# Use ACS Subject Tables to loand data into  R
# Subject Tables can be found in 
# https://www.census.gov/acs/www/data/data-tables-and-tools/subject-tables/
# Median Income in the Past 12 Months	median income, income	S1903
# Compare to race, poverty, income.
# Jin Man Lee, PhD, DePaul University
# 6/5/2020

# packages
# install.packages("tidyverse")
# install.packages("tidycensus")
# census_api_key("YOUR KEY GOES HERE")  # Census API key. Obtain one at http://api.census.gov/data/key_signup.html

library(tidyverse)
library(tidycensus)
library(ggplot2)

# variables
year <- 2018

# FIPS codes for 7 counties in CMAP region
# counties <- c("031", "043", "089", "093", "097", "111", "197")
counties <- c("031")


# load full ACS variables list, if needed
#load_variables(year, "acs5/subject", cache = TRUE) %>% 
#  View()


# get descriptive variable names for table S1903
varnames <- load_variables(year, "acs5/subject", cache = TRUE) %>% 
  # get only records from this table
  filter(str_sub(name, 1, 100) == "S1903") %>%
  # remove unnecessary label text
  mutate(label = sub('^Estimate!!', '', label)) %>% 
  # drop unnecessary columns and rename
  select(variable = name, label)


# get occupation data
S1903 <- get_acs(geography = "tract", table = "S1903", cache_table = TRUE, 
                 year = year, state = "17", county = counties, survey = "acs5")
# household numbers and median hhincome
hhnum     <- filter(S1903,variable=="S1903_C01_001") %>%
             rename(n_hh = estimate,n_hh_se=moe) %>%
             group_by(GEOID) %>%
             select(GEOID,NAME,n_hh,n_hh_se)
hhincome  <- filter(S1903,variable=="S1903_C03_001") %>%
             rename(mhhincome = estimate,hhincome_se=moe) %>%
             group_by(GEOID) %>%
             select(GEOID,mhhincome,hhincome_se) %>%
             left_join(hhnum,hhincome, by=c("GEOID"="GEOID")) %>%
             filter(mhhincome > 0 )
 
# 
complete.cases(hhincome)
# Graphical comparisons
ggplot(hhincome, aes(x=mhhincome,color="blue"))  + 
      geom_histogram(bins=25,color = "white",fill="darkblue")  

ggplot(hhincome, aes(x=mhhincome,color="blue"))  + 
  geom_histogram(bins=25,color = "white",fill="darkblue")  +
  geom_vline(aes(xintercept=mean(mhhincome)),
               color="green", linetype="dashed", size=1)

ggplot(hhincome, aes(x=mhhincome))  + 
  geom_histogram(aes(y=..density..), colour="blue", fill="white") +
  geom_density(alpha=.2, fill="#FF6666")