Creating a housing vulnerability index using tidycensus

housingindex

library(tidycensus)
library(tidyverse)
library(sf)

#Select variables
v18<-load_variables(2018,"acs5",cache=TRUE)
vars<-v18 %>%
  filter(substr(name,1,6) %in% c("B25064","B03002","B25106","B19013"))
#write_csv(vars,"data/census_vars.csv") #Select just those variables I want
vars_sel<-read_csv("data/census_vars.csv")

acsdata<-get_acs(geography="tract",state="GA",variables=vars_sel$name) %>%
  left_join(vars_sel %>% rename(variable=name))

#Select group total variables
acsdata_totals<-acsdata %>%
  filter(group_var=="-99999") %>%
  select(-group_var) %>%
  rename(group_var=variable,
         group_est=estimate) %>%
  select(GEOID,group_var,group_est)

#Select non-rate variables (income/rent)
acsdata_other<-acsdata %>%
  filter(group_var=="-88888") %>%
  select(GEOID,var_short,estimate) %>%
  spread(var_short,estimate)

#Calculate percentage
acsdata_pct<-acsdata %>%
  filter(substr(group_var,1,1)=="B") %>%
  left_join(acsdata_totals) %>%
  mutate(est_pct=round(estimate/group_est*100,1)) %>%
  select(GEOID,var_short,est_pct) %>%
  spread(var_short,est_pct) %>%
  mutate(housestress=ownhouse1+ownhouse2+ownhouse3+ownhouse4+ownhouse5+renthouse1+renthouse2+renthouse3+renthouse4,
         more2_pop=more2_pop1+more2_pop2+more2_pop3) %>%
  select(GEOID,afam_pop,amind_pop,asian_pop,hisp_pop,more2_pop,nhpi_pop,wht_pop,other_pop,
         housestress,renthouse_all)

acsdata_all<-acsdata_pct %>%
  left_join(acsdata_other)

#Read in job loss data
jobs<-st_read("data/jobloss.gpkg") %>%
  select(GEOID,X000) %>%
  rename(jobloss=X000) %>%
  inner_join(acsdata_all)

#Create index based on norms 
jobs_index<-jobs %>%
  st_set_geometry(NULL) %>%
  mutate(nonwht=100-wht_pop,
         medinc=max(medinc,na.rm=TRUE)-medinc) %>% 
  select(-afam_pop:-other_pop,-medrent) %>%
  gather(jobloss:nonwht,key=var,value=value) %>%
  filter(is.na(value)==FALSE) %>%
  group_by(var) %>%
  mutate(est_z=(value-mean(value))/sd(value)) %>%
  ungroup() %>%
  group_by(GEOID) %>% 
  summarise(index=sum(est_z))

jobs_final<-jobs %>% left_join(jobs_index)

st_write(jobs_final,"data/vulnindex.geojson")