This Gist creates a plot of smoothed zip-code-tabulation-area-level median household income against the percent of non-White residents, percent of Black residents and percent of Hispanic residents. Expect to be depressed. But first, set up your Census API key following the instructions in the link that follows. See section 3.3. http://eglenn.scr…

depressing_plots

library(acs)
library(dplyr)
library(ggplot2)

zips <- acs::geo.make(zip = "*")

# Get median household income.
income <-
  acs::acs.fetch(
    endyear = 2014, span = 5,
    geography = zips, variable = "B19013_001"
  )@estimate
income_dat <-
  data.frame(
    zip = substr(rownames(income), 7, 11), income = unname(income),
    stringsAsFactors = FALSE, row.names = NULL
  ) %>%
  dplyr::filter(!is.na(income))

# Get ethnic profiles.
ethnicity <-
  acs::acs.fetch(
    endyear = 2014, span = 5,
    geography = zips, table.number = "B03002"
  )@estimate
ethnicity_dat <-
  data.frame(ethnicity) %>%
  dplyr::mutate(
    zip = substr(rownames(ethnicity), 7, 11),
    p_nonwhite = 1 - (B03002_003 + B03002_013) / B03002_001,
    p_black = (B03002_004 + B03002_014) / B03002_001,
    p_hispanic = B03002_012 / B03002_001
  )

demog_dat <- dplyr::left_join(income_dat, ethnicity_dat)

ggplot2::ggplot(demog_dat, aes(p_nonwhite, income)) + geom_smooth()
ggplot2::ggplot(demog_dat, aes(p_black, income)) + geom_smooth()
ggplot2::ggplot(demog_dat, aes(p_hispanic, income)) + geom_smooth()