knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

### -------------
### packages
###

library(gganimate) # Used to create an animated plot
# install.packages('gifski') necessary for animation to compile
library(ggrepel) # Used to create cleveland dot plot
library(here)
library(kableExtra) # Used to create tables
library(patchwork) # Used to arrange plots as one image
library(sf) # Used to manipulate the spatial data for maps
library(tidycensus) # Used to pull census data 
library(tidyverse)
library(tigris) # Used to pull spatial features
library(UpSetR) # For plotting missing data
library(visdat) # For plotting missing data

### ------------------------
### helper functions + theme
###

summary_table <- function(data,
                          orient = "cols") { # cols or rows
  
  # This logic handles the orientation of the resulting table
  # "cols" means variable names are in the header
  # "rows" means variable names are their own column
  x = "Variable"
  y = ".value"
  
  if(orient == "cols") {
    x = ".value"
    y = "Statistic"
  }
  
  data %>%
    summarise(# calculates all other summary statistics
      across(.cols = everything(),
             .fns = list(Minimum = min, Mean = mean, `Standard Deviation` = sd, # Names the functions in the output
                         Median = median, Max = max), # "Name" = function
             .names = "{.col}_{.fn}", 
             na.rm = T)) %>%  # Passes na.rm = T to all of the listed functions
    pivot_longer(cols = everything(),
                 names_to = c(x, y),
                 names_pattern = "(.*)_(.*)") %>%
    mutate(across(where(is.numeric), round, digits = 2))
  
}

# Creates table to adjust 2014 census dollar values to 2019 values
# old code
# inflation_table <- blscrapeR::inflation_adjust(base_year = 2019) %>%
#   mutate(year = as.numeric(year)) %>%
#   filter(year %in% c(2014:2020))
# for CUSRxxxxxxx and CUURxxxxxxx (pre-1947) BLS variables

base_year = 2019

inflation_temp <- blscrapeR::bls_api("CUSR0000SA0", startyear = 2014) %>%
  subset(period!="M13" & 
           period!="S01" & 
           period!="S02" & 
           period!="S03") %>%
  mutate(date=as.Date(paste(year, period,"01",sep="-"),"%Y-M%m-%d"), 
         year=format(date, '%Y')) %>% 
  select(date, period, year, value) 

inflation_table <- 
  inflation_temp %>%
  group_by(year) %>%
  summarize(avg_cpi = mean(value)) %>%
  mutate(
    adj_value = avg_cpi / as.numeric(which(year==as.character(base_year))),
    base_year = as.character(base_year),
    pct_increase = (1-adj_value) * -100,
    adj_value = round(adj_value, 2)) %>%
  filter(year %in% c(2014:2020))

# Defines acs5_pull function
# Creates list of name-code pairs 
# Defines final variable calculations
# NOTE: To change the variables you can find codes using tidycensus::load_variables(year, "acs5")
source(here("census_tools.R"))

# Custom theme
custom_theme <- theme_minimal() +
  theme(plot.title = element_text(size = 16),
        plot.subtitle = element_text(size = 12),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12))

### ----------------
### research data
###

# Ingests and cleans the HUD data
# Creates the object hud_data_raw in env
# Source: https://www.huduser.gov/portal/datasets/fmr.html
# Methodology: https://www.huduser.gov/portal/sites/default/files/pdf/fmr-overview.pdf
# NOTE: some years have more ZIP codes than others in this data
source(here("hud_data.R")) # creates hud_data_raw object in env

# Zillow Observed Rent Index (ZORI) data
# Source: https://www.zillow.com/research/data/
# NOTE: Not every ZIP code in Chicago is included in this data
zori_raw <- read_csv(here("Data", "ZORI_AllHomesPlusMultifamily_SSA_ZIP.csv")) %>%
  filter(RegionName %in% hud_data_raw$ZIP) # 92 ZIP Codes

Introduction

Housing is a social determinant of health. Affordable housing is related to better health outcomes for tenants, freeing money up for better nutrition, daycare and healthcare, and other necesities for individual wellbeing. Evictions are related to poor health outcomes. Housing instability has been found to lead to higher stress levels, increased rates of depression and anxiety, and increased likelihood of living in lower quality housing or homelessness. Because of this, it is important to understand current metrics for what constitutes affordable housing. This project explores housing affordability through federal government metrics used to calculate housing subsidies known as Fair Market Rent (FMR).

Fair Market Rent is determined by a combination of a census survey of rents paid by recent movers, rent inflation adjustment, and forecasted expected growth in the rental housing market. FMR includes contract rent and utilities, ensuring that a full cost of base living necessities are captured in the metric. FMR calculations are made using American Community Survey data from two years prior to the calculation and are primarily used to set public housing flat rates and Section 8 housing choice voucher payment standards. Housing choice voucher recipients and public housing tenants make up approximately 3 million, or 2.5%, of United States households. Because of the diversity of contexts in which it is utilized, HUD FMR calculations can serve as a lens into the effectiveness of housing markets across the country in addressing issues of housing provision. This project aims to describe public data sources that can be used in analyzing housing affordability and some simple indicators that can be applied at various scales using Chicago as a case study.

Research Questions

  • What public data exists that can be used to develop indicators for housing affordability?

  • Are the U.S. Department of Housing and Urban Development’s Fair Market Rent calculations representative of actual rents in Chicago?

  • How have Fair Market Rent calculations and Zillow Observed Rent Index values changed over the period of 2014 to 2019 in Chicago?

  • Do socio-economic characteristics such as rent burden, gini index values, and educational attainment have an association with housing affordability?

Data

To perform this analysis, we needed data from three sources:

  • Fair Market Rent (FMR) data from the United States Department of Housing and Urban Development (HUD)

  • Observed rent data provided by Zillow

  • Socio-economic data provided by the United States Census in their American Community Survey data product

  • Geographic data describing the boundaries of ZIP codes within our study area provided by the Census Bureau

All of these data needed to be available at the same geographic level, aggregated to the ZIP code level. For this project we focused our efforts on ZIP codes within the City of Chicago. HUD FMR data is available through the HUD Policy Development & Research Website, Zillow data is available through their website’s research data page, and socio-economic data was available through the Census data API which we accessed using Kyle Walker’s TidyCensus package.

Methods

The following section describes methods used to wrangle and analyze the public data described above. All of the following methods were performed using R and various open source packages developed by the R community.

Data Wrangling

Housing Market Data

Zillow’s ZORI data provides average observed rents without regard for the number of bedrooms within. HUD provides FMR calculations for each bedroom size. To compare these two datasets we averaged the FMR for each bedroom size (studio through 4 bedroom) together to create an average FMR for the ZIP code. This method was chosen for its relative simplicity, with proposed alternatives including a weighted average with weights tied to bedroom size or the selection of one bedroom size as the standard for comparison.

HUD FMR data was collected at the ZIP code level, with calculations being made for each year in the study period. Zillow data was provided for each month within the study period. To compare data across the same unit of time, we transformed the ZORI data and averaged the monthly values to approximate a year average for each ZIP code. All dollar values were adjusted to 2019 values after calculations were completed.

HUD data was clean and complete, meaning that no rows contained any missing values and represented 487 ZIP codes within the Chicago metropolitan region across the 7 year study period. Zillow data was relatively more limited, sharing only 92 ZIP codes with the HUD data and missing many values within those ZIP codes for certain time periods within our study period. Missingness is discussed further in the Conclusions section.

### --------------
### hud data
###

hud_data <- hud_data_raw %>%
  # Create a variable to compare with ZORI values
  mutate(br_avg = (br_0 + br_1 + br_2 + br_3 + br_4) / 5) %>% 
  select(!c(br_0:br_4)) %>%
  # saw some duplicates on the year-zip level
  distinct()

### ------------
### zori data
###

zori_staging <- zori_raw %>%
  filter(RegionName %in% hud_data_raw$ZIP) %>% # 92 ZIP Codes
  pivot_longer(cols = "2014-01":"2020-12",
               names_to = c("year","month"),
               names_sep = "-",
               values_to = "rent",
               names_transform = list(year = as.double)) %>%
  select("RegionName","year","month","rent") 

zori <- zori_staging %>%
  group_by(RegionName, year) %>%
  summarize(rent_avg = mean(rent, na.rm = TRUE)) %>%
  ungroup()

### -----------
### market data
###

data <- zori %>%
  left_join(hud_data, by = c("RegionName" = "ZIP", "year" = "year")) %>%
  rename(hud_avg = br_avg, zori_avg = rent_avg) %>%
  # Adjust dollar values to 2019 dollars
  mutate(across(c(hud_avg, zori_avg),
                ~ .x / (inflation_table %>% filter(year == year) %>% pull(adj_value))))
  
# Output: table
data %>% 
  select("Zillow Observed Rent Index Value" = zori_avg,
         "Fair Market Rent Value" = hud_avg) %>%
  drop_na() %>%
  summary_table() %>% 
  # Display as currency
  mutate(across(.cols = 2:3, 
                .fns = ~ scales::dollar(.x))) %>%
  # basic table arrangements, label aesthetics
  kable(escape = F, # Tells cell_spec HTML to process instead of print
        align = "lrr",
        caption = "Housing Market Data Summary",
        format.args = list(big.mark = ",")) %>%
  # full-table aesthetics
  kable_styling(font_size = 20,
                html_font = "Lucida",
                full_width = F) %>%
  # header
  row_spec(0, bold = T, background = "#e5e5e5")
Housing Market Data Summary
Statistic Zillow Observed Rent Index Value Fair Market Rent Value
Minimum $17.98 $20.84
Mean $37.82 $33.52
Standard Deviation $8.74 $6.44
Median $36.25 $32.64
Max $71.82 $46.92

Socio-Economic Data

To place the housing market data into better context we combined it with socio-economic data provided by the US Census Bureau’s American Community Survey. ZIP code-level data is provided only in their 5 year estimate data products, which limited our ability to thoroughly analyze changes over 1 year time periods. To accommodate this limitation we pulled the data for two non-overlapping 5 year estimates, 2014 and 2019, and used them to supplement our findings. This is the recommended method for comparing ACS estimates according to the Census Bureau. 2014 dollar values were adjusted to 2019 values to account for inflation.

### -----------------------
### Census data
###

zips <- data %>%
  pull(RegionName) %>%
  unique()

acs_2014 <- acs5_pull(2014, zips) %>%
  # Including zips in the function call is supposed to prevent us needing
  # to filter by zip, but it's not working.
  mutate(ZIP = str_extract(NAME, "[0-9]{5}?")) %>%
  filter(ZIP %in% zips)

acs_2019 <- acs5_pull(2019, zips) %>%
  mutate(ZIP = str_extract(NAME, "[0-9]{5}?")) %>%
  filter(ZIP %in% zips)

data_2014 <- data %>%
  filter(year == 2014) %>%
  left_join(acs_2014, by = c("RegionName" = "ZIP", "year")) %>%
  mutate(median_hh_income =  median_hh_income / (inflation_table %>% filter(year == year) %>% pull(adj_value)))

data_2019 <- data %>%
  filter(year == 2019) %>%
  left_join(acs_2019, by = c("RegionName" = "ZIP", "year"))

# Output: tables
data_2014 %>%
  select("Gini Index" = gini_index,
         "Rent Burdened Households" = percent_unaffordable,
         "Median Household Income (USD)" = median_hh_income,
         "Black Population" = percent_black_pop) %>%
  drop_na() %>%
  summary_table() %>%
    # Display as currency
  mutate(`Median Household Income (USD)` = scales::dollar(`Median Household Income (USD)`), 
         across(.cols = c(3, 5),
                .fns = scales::percent,
                accuracy = 1)) %>%
  # basic table arrangements, label aesthetics
  kable(escape = F, # Tells cell_spec HTML to process instead of print
        align = "lrrrr",
        caption = "2014 Census Data Summary[note]",
        format.args = list(big.mark = ",")) %>%
  # full-table aesthetics
  kable_styling(font_size = 20,
                html_font = "Lucida",
                full_width = F) %>%
  # header
  row_spec(0, bold = T, background = "#e5e5e5") %>%
  add_footnote("Dollars in 2019 values")
2014 Census Data Summarya
Statistic Gini Index Rent Burdened Households Median Household Income (USD) Black Population
Minimum 0.32 30% $529.32 1%
Mean 0.45 47% $1,659.76 18%
Standard Deviation 0.06 8% $625.05 27%
Median 0.45 46% $1,677.72 6%
Max 0.58 65% $3,864.76 97%
a Dollars in 2019 values 
data_2019 %>%
  select("Gini Index" = gini_index,
         "Rent Burdened Households" = percent_unaffordable,
         "Median Household Income (USD)" = median_hh_income,
         "Black Population" = percent_black_pop) %>%
  drop_na() %>%
  summary_table() %>%
    # Display as currency
  mutate(`Median Household Income (USD)` = scales::dollar(`Median Household Income (USD)`),
         across(.cols = c(3, 5),
                .fns = scales::percent,
                accuracy = 1)) %>%
  # basic table arrangements, label aesthetics
  kable(escape = F, # Tells cell_spec HTML to process instead of print
        align = "lrrrr",
        caption = "2019 Census Data Summary",
        format.args = list(big.mark = ",")) %>%
  # full-table aesthetics
  kable_styling(font_size = 20,
                html_font = "Lucida",
                full_width = F) %>%
  # header
  row_spec(0, bold = T, background = "#e5e5e5")
2019 Census Data Summary
Statistic Gini Index Rent Burdened Households Median Household Income (USD) Black Population
Minimum 0.31 22% $23,429 1%
Mean 0.45 44% $82,324 17%
Standard Deviation 0.06 8% $33,092 26%
Median 0.46 44% $83,028 6%
Max 0.64 61% $191,528 96%

Spatial Data

There may be a spatial component that would be necessary to understand the dynamics present within our data. To explore this we used the tigris package to acquire spatial data that would allow us to map the research data to real space. Spatial data was then joined to our research data for the purpose of visualizing the changes in the data over time in a way that observes geographic differences. The following map shows the location and area of ZIP codes included in the research data.

### -----------------------------
### Spatial Data
### 

# Transforming to this CRS makes the map look "normal"
crs <- 3435 # coordinate reference system: IL NAD83 East

# Get Cook County geometry to put zctas in context
cook_county <- counties(state = "IL", cb = T, progress_bar = F) %>%
  filter(NAME == "Cook") %>% # Chicago + suburbs are almost entirely within Cook County
  select(NAME, geometry) %>%
  st_transform(crs)

# Create filter for ZIP geom pull
starts_with <- data %>% 
  pull(RegionName) %>% # Vectorize the column
  str_extract(pattern = "[0-9]{2}") %>%
  unique() # c("53", "60") 

### Get ZCTA geometry
tigris_cache_dir(here())
zcta <- zctas(cb = F, starts_with = starts_with, 
              progress_bar = F, refresh = T,
              keep_zipped_shapefile = T) %>%
  # changed from GEOID10. Was the variable or API updated?
  select(GEOID20, geometry) %>%
  st_transform(crs) %>% # Without this transformation to ILNAD83 the map looks weird
  st_filter(cook_county) # There are ZIP codes that fall outside of Cook county
                         # and make the map harder to digest. st_filter removes them

# Output: Map
data %>%
  left_join(zcta, by = c("RegionName" = "GEOID20")) %>%
  ggplot(aes(geometry = geometry)) +
  geom_sf(data = cook_county) + # County boundary
  geom_sf(aes(fill = 'red')) + # ZCTA boundaries
  labs(title = "Map of ZIP Codes in Research Data",
       subtitle = "ZIP Codes in red, Cook County in grey") +
  custom_theme +
  guides(fill = FALSE) # remove legend

Analysis

Pulling Out the Magnifying Glass

The general trends described above provide insight into Chicago as a whole, opening up questions about whether these trends are experienced equally across Chicago ZIP codes.

An article written by Deborah Thrope in the Journal of Affordable Housing titled Achieving Housing Choice and Mobility in the Voucher Program outlined the inadequacies of the Housing Choice Voucher program in enabling tenant mobility and prescribed a collection of policy recommendations to improve outcomes. Among them was a recommendation that HUD use Small Area FMR (SAFMR) values to better account for local market conditions.

The following animation visualizes the difference between observed rent and FMR values over time for ZIP codes (Note: SAFMR and ZIP Code are used interchangeably in this report) where data was available. Please note that because of the way the difference was calculated it is not immediately apparent whether changes are due to a shift in HUD data, Zillow data, or both.

### -----------------------------
### Year-over-Year difference map
### Animated

map_anim <- data %>%
  mutate(diff = zori_avg - hud_avg) %>% 
  # changed from GEOID10 here as well
  left_join(zcta, by = c("RegionName" = "GEOID20")) %>%
  ggplot(aes(geometry = geometry)) +
  geom_sf(data = cook_county) + # County boundary
  geom_sf(aes(fill = diff)) + # ZCTA map
  # Adds the year value to the title through interpolation
  labs(title = "Difference in FMR Verus Zillow Rents for Year: {current_frame}", 
       subtitle = "in 2019 dollars", 
       fill = "ZORI - HUD Value ($$)") + 
  scale_fill_viridis_c() +
  custom_theme +
  transition_manual(year) # transition_time doesn't work well with sf for some reason

animate(map_anim, fps = 1, nframes = 6)

The following distribution plots help us understand that the change in the difference between ZORI values and HUD FMR is an effect of FMR increasing, not ZORI decreasing. What does this mean for housing affordability? This could indicate that HUD voucher recipients are being provided a wider range of places to move to as HUD FMR begins to match average rents.

### --------------------------
### ZIP/Rent distribution plot
### Animated

dist_anim <- data %>% 
  group_by(RegionName, year) %>%
  summarize(ZORI = mean(zori_avg, na.rm = T),
            HUD = mean(hud_avg, na.rm = T)) %>%
  ggplot() +
  geom_density(aes(x = ZORI), fill = "blue", alpha = .5) +
  geom_density(aes(x = HUD), fill = "red", alpha = .5) +
  labs(title = "Distribution of Rent Values by ZIP Code, {current_frame}",
       subtitle = "blue is ZORI, red is HUD",
       x = "Rent", 
       y = "Density") +
  custom_theme +
  transition_manual(year)

animate(dist_anim, fps = 1, nframes = 6)

The above animations show that not every SAFMR value is equally responsive to market conditions, even when larger trends indicate progress towards intended outcomes. ZIP Codes in some areas saw relatively little change in the difference between their observed rents and calculated FMR values, while others saw drastic changes over the same time period - shifting from a large positive value to a negative value. This may indicate inefficiencies in the data and methods used to calculate SAFMR since an ideal model of FMR calculation would result in a relatively homogeneous set of difference values across all areas. Further exploration of these differences could take the form of a case study of areas where the difference value was positive and saw little change across the study period compared with those where the difference value was close to negative and saw little change across the study period. The following section attempts to explore a methodological avenue for doing such a case study.

Examining the Edges: Largest Differences in FMR Calculations and Observed Rents

The following section performs a brief case study of the areas with the largest differences at the beginning of our study period, 2014. One notable theme in this set was gentrification. Five of the ten ZIP codes in this case study contained areas that were known to be gentrifying or at risk of gentrification in this study performed by the Nathalie P. Voorhees Center at UIC - Logan Square, Pilsen, Cabrini Green, Chinatown, and East Garfield Park. Expanding on this case study could provide more explanatory or descriptive factors in gentrification research.

To add some context to this issue we combined our market data with socio-economic data and examined the trend in income changes over time as compared to rent for the ZIP codes with the largest difference in ZORI and FMR values in 2014. Because the Census Bureau recommends not comparing overlapping 5-year survey estimates, and because ZIP code-level data is unavailable in 1-year survey data, this approach is unable to provide a granular time-series analysis of the relationships between socio-economic characteristics and the market data. The following table lists those ZIP codes and their associated town or neighborhood. Label values were defined using Google Maps and assigning names based on which names appeared within the ZIP code boundaries.

# Defines labels for ZIP codes
neighborhoods <-
  list(
    "60045" = "Lake Forest, IL",
    "60602" = "Downtown, Washington",
    "60610" = "Near Northside, Cabrini Green",
    "60126" = "Elmhurst, IL",
    "60642" = "River West",
    "60647" = "Logan Square",
    "60616" = "Chinatown, Armour Square, Bridgeport",
    "60525" = "La Grange, IL",
    "60612" = "East Garfield Park",
    "60608" = "Pilsen",
    "60538" = "Aurora, IL",
    "60601" = "Downtown, Lake",
    "60611" = "Streeterville",
    "60654" = "River North"
  )

 data %>%
  mutate(diff = (zori_avg - hud_avg) / hud_avg) %>%
  ungroup() %>%
  filter(year == 2014) %>%
  arrange(desc(diff)) %>%
   slice(head(row_number(), 10)) %>%
   select(RegionName) %>%
   mutate(Neighborhood = as.character(neighborhoods[RegionName])) %>%
   # basic table arrangements, label aesthetics
   kable(escape = F, # Tells cell_spec HTML to process instead of print
         align = "lr",
         caption = "ZIP Codes and Associated Neighborhoods") %>%
   # full-table aesthetics
   kable_styling(font_size = 20,
                 html_font = "Lucida",
                 full_width = F) %>%
   # header
   row_spec(0, bold = T, background = "#e5e5e5")
ZIP Codes and Associated Neighborhoods
RegionName Neighborhood
60045 Lake Forest, IL
60616 Chinatown, Armour Square, Bridgeport
60525 La Grange, IL
60642 River West
60602 Downtown, Washington
60126 Elmhurst, IL
60610 Near Northside, Cabrini Green
60538 Aurora, IL
60608 Pilsen
60647 Logan Square

The following plots visualize the change in the selected ZIP codes between 2014 and 2019 in terms of observed minus fair market rent values and rent burden. Rent burden is defined by HUD as a tenant’s gross rent being greater than 30% of their income. Any rent above 50% of tenant income is considered to be a severe rent burden. The lines between the plotted values are colored based on whether the proportion of households in the ZIP code area were majority or minority rent burdened.

### ------------------
### Cleveland Dot Plot 
### 

cleveland_data <-
  data %>%
  mutate(diff = (zori_avg - hud_avg) / hud_avg) %>%
  ungroup() %>%
  filter(RegionName %in% names(neighborhoods) & year %in% c(2014,2019)) %>%
  pivot_longer(cols = c(zori_avg:hud_avg),
               names_to = "source",
               values_to = "rent") %>%
  mutate(neighborhood = neighborhoods[RegionName])

acs_to_year_mapping <- list(
  "2014" = data_2014,
  "2019" = data_2019
)

cleveland_plot <- function(input_data, yr, acs) {
  plot_data <-
    input_data %>%
    left_join(acs[[as.character(yr)]], by = c("RegionName", "year")) %>%
    filter(year == yr) %>%
    mutate(burden = cut(percent_unaffordable,
                        breaks=c(0, 0.5, 1),
                        labels=c("Minority","Majority")))
  
  ggplot(data = plot_data,
         mapping = aes(x = RegionName, y = rent)) +
    # line showing the difference between the HUD and Zillow indices
    geom_line(aes(group = RegionName, color = burden), size = 5) +
    # shapes representing either HUD or Zillow index
    geom_point(aes(shape = source), size = 5) +
    # labeling the difference by percentage between indices
    geom_text(
      data = plot_data %>% filter(source == "zori_avg"),
      mapping = aes(label = paste0("", scales::percent(x = percent_unaffordable, 
                                                       accuracy = 0.01)),
                    color = burden),
      size = 3.5,
      hjust = -0.25
    ) +
    # labeling the lines
    geom_text(
      data = plot_data %>% filter(source == "hud_avg"),
      mapping = aes(label = neighborhood),
      size = 3.5,
      hjust = 0,
      nudge_x = 0.5
    ) +
    scale_y_continuous(labels = scales::dollar_format(),
                       limits = c(800, 3200)) +
    coord_flip() +
    labs(
      title = paste0("Fair and Market Rate Differences by ZIP, ", yr),
      subtitle = "HUD versus Zillow with Proportion of Population Rent Burdened",
      x = "ZIP",
      y = "Rent (in 2019 $$)",
      shape = "Rent Value",
      color = "Percent Rent Burdened"
    ) +
    # Proper labels with aesthetics
    scale_color_manual(values = c("Minority"="#00ba38","Majority"="#F8766D")) +
    scale_shape_manual(values = c("circle","square"), labels = c("HUD", "Zillow")) +
    custom_theme
}

cleveland_plot(cleveland_data, 2014, acs_to_year_mapping) / cleveland_plot(cleveland_data, 2019, acs_to_year_mapping)

As expected, all of the neighborhoods saw an increase in their observed rent and subsidy values. The two Chicago neighborhoods that were majority rent burdened saw rent burdened households become the minority, though the shift was only by a few points. Aurora stayed majority rent burdened between the two study periods, indicating that households in Aurora continue to face issues related to affordable housing. It’s important to note that the separation between majority and minority labeling of the population as burdened is at 50%. This means that even though the historically marginalized neighborhoods, like Pilsen and Garfield Park, shifted into being minority rent burdened in 2019, they are still over 40% rent burdened!

Socio-Economic Contexts

To further explore the data, we put together correlation plots describing the relationship between the difference values and three census variables: Median Household Income, Gini Index Coefficient, and Percent Black Population. Note that no correlation tests were performed in this section, this is a sequence of plots without underlying statistics.

### -------------------
### coor plots
### hh income

plot1 <- data_2014 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = median_hh_income, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Household Income, 2014",
       subtitle = "dollars adjusted to 2019 values",
       x = "Median Household Income",
       y = "ZORI - FMR Values (in USD)") +
  scale_x_continuous(limits = c(0, 190000)) + 
  custom_theme

plot2 <- data_2019 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = median_hh_income, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Household Income, 2019",
      x = "Median Household Income",
      y = "ZORI - FMR Values (in USD)") +
  scale_x_continuous(limits = c(0, 190000)) + 
  custom_theme

plot1 / plot2

Correlation between difference values and median household income were positive for both years. This makes sense since areas with higher household incomes are likely to have newer, more expensive housing stock.

### -------------------
### coor plots
### gini index

plot3 <- data_2014 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = gini_index, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Gini Index, 2014",
       subtitle = "dollars adjusted to 2019 values",
       x = "Gini index",
       y = "ZORI - FMR Values (in USD)") +
  scale_x_continuous(limits = c(0.3, .65)) + 
  custom_theme

plot4 <- data_2019 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = gini_index, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Gini Index, 2019",
      x = "Gini Index",
      y = "ZORI - FMR Values (in USD)") +
  scale_x_continuous(limits = c(0.3, .65)) + 
  custom_theme

plot3 / plot4

These plots demonstrate a positive correlation between Gini Index Coefficient values and difference values. Though the relationship was less strong in 2019 than 2014. The difference value indicates a lack of adequate housing subsidy for low income housholds. This correlational change makes sense given the decrease in difference values between the two time periods. There are multiple possible explanations for this change ranging from a decrease in income inequality, a shift in the makeup of subsidized households in the area, or a displacement of the lowest income households from the study areas. Further research on this phenomenon seems warranted in this circumstance.

### -------------------
### coor plots
### share black pop

plot5 <- data_2014 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = percent_black_pop, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Share Black Population, 2014",
       subtitle = "dollars adjusted to 2019 values",
       x = "Black Population",
       y = "") +
  scale_x_continuous(limits = c(0, 1),
                     labels = scales::percent) + 
  custom_theme

plot6 <- data_2019 %>%
  mutate(diff = zori_avg - hud_avg) %>%
  ggplot(aes(x = percent_black_pop, y = diff)) +
  geom_jitter() +
  geom_smooth(method = "lm") +
  labs(title = "Observed Minus Fair Market Rent Compared with Share Black Population, 2019",
      x = "Gini Index",
      y = "") +
  scale_x_continuous(limits = c(0, 1),
                     labels = scales::percent) + 
  custom_theme

plot5 / plot6

These plots show that the relation between difference values and the proportion of Black population in the area is negative. The majority of ZIP codes in the study areas have a sub-25% proportion of Black residents. Those that are majority Black are associated with lower or negative difference values.

Overall, these findings point to three possible independent variables that could be used in a regression analysis of this data. Further investigation of these relationships using statistical tests would be warranted given these exploratory findings.

Conclusions

The above analysis provides an analytic approach to exploring the deficiencies of the federal government in addressing issues of access to affordable housing. Our research found that household rent subsidies vary between ZIP Code areas, though generally Fair Market Rent calculations have been catching up with actual observed rents. Areas with the widest differences between observed and calculated fair market values were spaces that have been described as gentrifying by other researchers. Visual correlations were found between measures of segregation and inequality and the difference between FMR and observed rent values.

A strength of this analysis is that, because of the public-facing nature of the data, our methods are able to be applied in other municipalities or across municipalities. We hope that this project serves as a proof of concept for a broader analysis of housing subsidies in the United States.

Policy Implications

Further research and literature review is warranted given the findings of this project. Initial findings demonstrate a need for more efficient calculations of Fair Market Rent to expand the quantity and quality of rental units available to low-income tenants. This research also pointed to possible relationships between rental housing subsidy availability and gentrification, implying a need for more cooperation between federal and local housing authorities in addressing issues related to gentrification and segregation.

Further Questions

Given the results of our correlation plots, possible relationships between those values and housing subsidies could be tested using statistical tests and in-depth time-series regression analysis. Further, research in urban planning may benefit from comparing measures of gentrification with housing subsidy data to better understand dynamics of displacement, disinvestment, and dispossession in urban environments.

Limitations

Our data was limited to a subset of ZIP codes within Cook County, and because of this our findings may not be generalizable to the entire City of Chicago or Cook County as a whole. ZORI data had higher degrees of missingness in earlier time periods than later ones, indicating possible hurdles to overcome in a time-series analysis.

Missingness in data

zori_raw %>%
  select(c("2014-01":"2021-02")) %>%
  vis_miss() +
  geom_vline(xintercept = seq(from = 0, to = 95, by = 12), size = 1) +
  annotate(
    geom = "text",
    x = seq(from = 6, to = 88, by = 12),
    y = 95,
    label = c(2014:2020),
    size = 8,
    angle = 0
  ) +
  custom_theme +
  theme(axis.text.x = element_text(size = 10, angle = 90),
        axis.title = element_text(size = 16),
        title = element_text(size = 40)) +
  theme(axis.text.x = element_blank()) +
  annotate(
    geom = "text",
    x = seq(from = 1.5, to = 87, by = 3),
    y = -2.5,
    label = c(rep(x = c("Q1","Q2","Q3","Q4"), time = length(c(2014:2020)), each = 1),"Q1"),
    size = 4,
    angle = 90
  ) +
  labs(title = "Missing Monthly Data from Zillow",
       subtitle = "Split by Quarter and Year")

We see here that missing data is much more of a concern for 2014 and 2015, especially in the upper range of rents, but less of an issue for the subsequent years. This missing data is not at random, and may be due to a lack of sufficient quantity of listings in earlier years. 2014 had far more missing values than any subsequent year.

After the join between HUD and Zillow data, we find that Zillow has data on some years that HUD doesn’t, but it was not enough to cause concern with only 3 values missing in the three most recent years.

ggplot(data %>%
         group_by(year) %>%
         summarise(n_missing = sum(is.na(hud_avg)))) +
  geom_bar(aes(x = year, y = n_missing), stat = "identity") +
  scale_x_continuous(breaks = data$year) +
  scale_y_continuous(limits = c(0, 5)) +
  labs(x = "Year",
       y = "# of Missing Values",
       title = "Missing Data from HUD") +
  custom_theme