PUMS Housing Affordability Analysis

TidyCensus reference:

https://walker-data.com/census-r/index.html

See, particularly, chapters 9 and 10.

Sampling at the “household” level

TidyCensus extracts only individual-level PUMS data, meaning that each row represents an individual person. That’s great if you are investigating individual-level traits but problematic if (like us) you are investigating household-level traits. If a houshold includes only one individual, then the household occupies only one row in the dataset, and the household’s traits are the same as the household’s traits. But if a household includes two or more individuals, then the household occupies multiple rows in the dataset (one for each household occupant). This possibility poses at least two problems for a household-level analysis. First, the number of rows equals the number of people, not the number of households. So if you treat each row as a household - and especially if you involve weights - you will drastically overestimate the number of households. Furthemore, for multiple-occupant households, no single set of individual traits (like age, race, gender, etc.) can correctly describe the whole houshold.

One way to solve this problem is to filter the dataset for individuals designated as “1” on the SPORDER variable. SPORDER contains the “person number” for each household, and the Census Bureau’s questionnaire instructs the person completing the survey to designate “Person 1” as “the person living or staying here in whose name this house or apartment is owned, being bought, or rented,” and adds that, “If there is no such person, start with the name of any adult living or staying here.” (See: The 2024 ACS questionnaire.) So, filtering for rows where SPORDER equals 1 leaves you with one row per household, with individual-level traits that represent the “head” of the household.

In the script, the filtering is accomplished but this code, which fetches and then filters the data:

# =========================================================
# 9. Pull renter data
# =========================================================

TN_Renters <- get_pums(
  variables = vars,
  state = "TN",
  survey = "acs1",
  variables_filter = list(TEN = 3),
  year = 2024,
  recode = TRUE
)

TN_Renters <- TN_Renters |> 
  filter(SPORDER == 1)

Measuring “rent burden”

The GRPIP variable must be treated carefully. In the raw data, it is a numeric variable with values from 0 to 100 representing the percentage of household income spent on rent, plus a 101 category representing “topcoded” households for which a meaningful GRPIP value cannot be computed. Treating it numerically - as if “101” means a household is spending 101 percent of its income on rent - wouldn’t make sense. One also could question whether “100” households are truly spending 100 percent of their household income on rent, or how anyone could be getting by while spending north of 70, 80 or 90 percent of their income on rent. About the only workable scenario I can think of involves a renter with an income so high that even when 70-plus percent of it goes toward rent, the remainder is enough to cover other household needs, like food, transportation, clothing, medical care, etc.

Here is what the distribution of raw GRPIP values looks like. Note that the values are sorted alphabetically, not numerically, which is why the “101” category shows up between the 100 and 11 categories.

---

---

Given this ambiguity about what high GRPIP values actually mean, it’s probably wise to use GRPIP to group renters into four broad categories:

• The “Not burdened,” who spend less than 30 percent of their household income on rent.

• The “Cost burdened,” who spend 30-49% of their household income on rent.

• The “Severely burdened,” who spend 50% or more of their household income on rent.

• The “Extreme / top-coded,” who for a range of reasons that may not be similar or valid, show up as spending 100% or more on rent.

Below is the distribution of these four categories, with weights applied so that the category totals estimated totals for all renter households in Davidson County. Step 14 in the code shows how the categories were created.

---

---

Exploring rent burden across demographics

Continuous and categorical variables. The PUMS data contain a range of demograpic variables. Some, like “Age,” are numeric in nature, or “continuous.” They measure countable quantities, like the number of years that have passed since a person was born. Others, like “Race,” are “categorical.” They sensibly place people into categories based on a some trait that is either fully present or fully absent, not present in amounts ranging from none to however much is possible. In an extremely important sense, nearly nobody is wholly “white” or wholly “black” or wholly “Asian.” Nearly everyone has ancestors from two or more racial groups, and “race” itself is is largely a socially constructed concept. But people do tend to describe themselves as belonging to distinct racial categories. For that reason, “Race” is typically treated as a categorical variable, often with at least one category for people who consider themselves to be a blend of two or more races. It is also important to note that “Race” and “ethnicity” are not the same thing. Ethnicity has more to do with cultural traits, like one’s language, preferred food, musical tastes, and (sometimes) religion. Often, people of different races can belong to the same ethnic group, and vice versa.

Typically, continuous variables are summarized with average or a median, usually with some indication of the variable’s variability or distribution shape. Categorical variables, meanwhile, are usually summarized using group percentages. However, it might be wise for us to convert all continuous variables into categorical ones, then use percentages as summaries. The reason: Averages - and, to a lesser degree, medians - imply that midrange values generally reflect all of the varible’s values. With demographics like age and income, this isn’t always the case.

Households and individuals. Another thing to consider: “Rent burden” is a household-level variable, in that all members of a household are affected by rent burden in more or less the same way. But as discussed above, demographics like age and race are individual-level variables, in that members of the same household can be different ages or belong to different racial groups. One way to address the mismatch is to focus on the demographic characteristics of a each household’s “householder.” The Census defines the “householder” as the person whose name is on the household’s deed or lease. If two or more people fit that criteria, the Census goes with whoever the members of the household designate as the “householder” or, if no such designation was made, infers the identity of the householder from answers to other questions in the ACS survey. The approach works fine, as long as we don’t try to generalize from the householder to all members of the household.

Here are graphics of rent burden broken down by the householder’s age and race categories:

---

---

---

Mapping rent burden

Rent burden can vary by geography as well as by householder demographics. Here’s a look at how rent burden compares across Nashville’s six Census-defined Public Use Microdata Areas:

---

---

Mapping current rent

The Zillow Observed Rent Index, or ZORI, from Zillow.com offers a more up-to-date look at current rents. It also organizes data by ZIP code. The index tracks asking rents while controlling for changes in the quality of available rental stock. Data are available for most, but not all, ZIP codes that fall at least partly within Metro Nashville:

---

---

Next steps

• Explore regional and/or national comparisons.

• Determine the best way to handle error margins.

PUMS Analysis Script:

# =========================================================
# 1. Load required libraries
# =========================================================

if (!require("tidycensus")) install.packages("tidycensus")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("tigris")) install.packages("tigris")
if (!require("sf")) install.packages("sf")
if (!require("leaflet")) install.packages("leaflet")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("plotly")) install.packages("plotly")
library(tidycensus)
library(tidyverse)
library(tigris)
library(sf)
library(leaflet)
library(kableExtra)
library(plotly)


# =========================================================
# 2. (Optional) Set Census API key
# =========================================================

# census_api_key("YOUR_API_KEY_HERE", install = TRUE)
# readRenviron("~/.Renviron")


# =========================================================
# 3. Load Tennessee PUMAs (2020)
# =========================================================

options(tigris_use_cache = TRUE)

TN_PUMAs <- pumas(state = "TN", cb = TRUE, year = 2020)


# =========================================================
# 4. Define Nashville PUMAs
# =========================================================

PUMA_list <- c("02401","02402","02403","02404","02405","02406")

NASH_PUMAs <- TN_PUMAs %>%
  filter(PUMACE20 %in% PUMA_list)


# =========================================================
# 5. Prepare geometry for Leaflet
# =========================================================

NASH_PUMAs <- st_transform(NASH_PUMAs, 4326)


# =========================================================
# 6. Create Leaflet map
# =========================================================

PUMA_map <- leaflet(data = NASH_PUMAs) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng = -86.78, lat = 36.17, zoom = 10) %>%
  addPolygons(
    fillColor = "steelblue",
    fillOpacity = 0.5,
    color = "black",
    weight = 1,
    highlightOptions = highlightOptions(
      weight = 2,
      color = "white",
      bringToFront = TRUE
    ),
    popup = ~paste0(
      "<strong>PUMA GEOID:</strong> ", GEOID20, "<br>",
      "<strong>PUMA code:</strong> ", PUMACE20, "<br>",
      "<strong>Name:</strong> ", NAMELSAD20
    )
  )

PUMA_map


# =========================================================
# 7. Define affordability variables (UPDATED)
# =========================================================

vars <- c(
  "TEN","GRNTP","RNTP","GRPIP","HINCP",
  "NP","NOC","AGEP","HHLDRAGEP",   # ✅ ADDED
  "RAC1P","HISP","HHLDRRAC1P",
  "SCHL","ESR",
  "ELEP","GASP","WATP","FULP",
  "RMSP","BDSP",
  "PUMA","WGTP"
)


# =========================================================
# Burden color palette
# =========================================================

burden_colors <- c(
  "1 Not burdened" = "#B3D7FF",
  "2 Cost burdened (30–49%)" = "#66A3FF",
  "3 Severely burdened (50–99%)" = "#0066CC",
  "4 Extreme / top-coded (100%+)" = "#003366"
)


# =========================================================
# 8. Codebooks
# =========================================================

data(pums_variables)

pums_codebook_full <- pums_variables %>%
  filter(survey == "acs1", year == 2024)

pums_codebook_affordability <- pums_codebook_full %>%
  filter(var_code %in% vars)


# =========================================================
# 9. Pull renter data
# =========================================================

TN_Renters <- get_pums(
  variables = vars,
  state = "TN",
  survey = "acs1",
  variables_filter = list(TEN = 3),
  year = 2024,
  recode = TRUE
)

TN_Renters <- TN_Renters |> 
  filter(SPORDER == 1)


# =========================================================
# 10. Filter Nashville renters
# =========================================================

NASH_Renters <- TN_Renters %>%
  filter(PUMA %in% PUMA_list)


# =========================================================
# 11. Raw GRPIP distribution (Plotly)
# =========================================================

grpip_counts <- NASH_Renters %>%
  count(GRPIP)

GRPIP_raw_plot <- plot_ly(
  grpip_counts,
  x = ~GRPIP,
  y = ~n,
  type = "bar",
  marker = list(color = "steelblue"),
  text = ~n,
  hovertemplate = "GRPIP: %{x}<br>Count: %{y}<extra></extra>"
) %>%
  layout(
    title = "Raw Distribution of GRPIP (Sample Counts)",
    xaxis = list(title = "GRPIP (Raw Values)"),
    yaxis = list(title = "Count")
  )

GRPIP_raw_plot


# =========================================================
# 12. Convert to numeric
# =========================================================

NASH_Renters <- NASH_Renters %>%
  mutate(across(
    c(GRNTP, RNTP, GRPIP, HINCP,
      NP, NOC, AGEP, SCHL,
      ELEP, GASP, WATP, FULP,
      RMSP, BDSP),
    as.numeric
  ))


# =========================================================
# 13. NA Summary Table
# =========================================================

na_summary <- NASH_Renters %>%
  summarize(across(
    c(GRNTP, GRPIP, HINCP, SCHL),
    ~ mean(is.na(.))
  )) %>%
  pivot_longer(
    everything(),
    names_to = "Variable",
    values_to = "Pct_Missing"
  )

na_summary %>%
  kable(digits = 3, caption = "NA Proportion by Variable") %>%
  kable_styling(full_width = FALSE)


# =========================================================
# 14. Affordability measures + burden groups
# =========================================================

NASH_Renters <- NASH_Renters %>%
  mutate(
    annual_rent = GRNTP * 12,
    residual_income = HINCP - annual_rent,
    college = SCHL >= 21,
    burden_group = case_when(
      GRPIP < 30 ~ "1 Not burdened",
      GRPIP < 50 ~ "2 Cost burdened (30–49%)",
      GRPIP < 100 ~ "3 Severely burdened (50–99%)",
      GRPIP >= 100 ~ "4 Extreme / top-coded (100%+)",
      TRUE ~ NA_character_
    )
  )


# =========================================================
# 15. Housing stress profile table
# =========================================================

profile_table <- NASH_Renters %>%
  group_by(burden_group) %>%
  summarize(
    households = sum(WGTP, na.rm = TRUE),
    avg_income = weighted.mean(HINCP, WGTP, na.rm = TRUE),
    avg_rent = weighted.mean(GRNTP, WGTP, na.rm = TRUE),
    pct_college = weighted.mean(college, WGTP, na.rm = TRUE)
  )

profile_table %>%
  kable(digits = 1, caption = "Housing Stress Profile by Burden Group") %>%
  kable_styling(full_width = FALSE)

profile_table


# =========================================================
# 16. Burden distribution (Plotly, weighted + labeled)
# =========================================================

burden_counts <- NASH_Renters %>%
  group_by(burden_group) %>%
  summarize(weighted_n = sum(WGTP, na.rm = TRUE))

burden_group_graph <- plot_ly(
  burden_counts,
  x = ~burden_group,
  y = ~weighted_n,
  color = ~burden_group,
  colors = burden_colors,
  type = "bar",
  text = ~round(weighted_n),
  textposition = "outside",
  hovertemplate = "%{x}<br>Households: %{y}<extra></extra>"
) %>%
  layout(
    title = "Distribution of Housing Burden (Weighted)",
    xaxis = list(title = "Burden Group"),
    yaxis = list(title = "Estimated Households"),
    legend = list(traceorder = "normal")
  )

burden_group_graph

# =========================================================
# 17. Summary statistics table
# =========================================================

summary_stats <- NASH_Renters %>%
  mutate(
    cost_burdened = !str_detect(burden_group, "Not burdened"),
    severe_burden = str_detect(burden_group, "Severely")
  ) %>%
  summarize(
    pct_cost_burdened =
      weighted.mean(cost_burdened, WGTP, na.rm = TRUE),
    pct_severe_burden =
      weighted.mean(severe_burden, WGTP, na.rm = TRUE),
    avg_rent =
      weighted.mean(GRNTP, WGTP, na.rm = TRUE),
    avg_income =
      weighted.mean(HINCP, WGTP, na.rm = TRUE)
  )

summary_stats_table <- summary_stats %>%
  kable(digits = 3, caption = "Summary Statistics (Weighted)") %>%
  kable_styling(full_width = FALSE)

summary_stats_table


# =========================================================
# 18. Save dataset
# =========================================================

saveRDS(NASH_Renters, "nashville_renters_acs1_2024.rds")

# =========================================================
# 19. Age-by-burden table and interactive visualization
# =========================================================

# ---- Create age groups ----
NASH_Renters <- NASH_Renters %>%
  mutate(
    age_group = case_when(
      HHLDRAGEP < 30 ~ "29 and under",
      HHLDRAGEP < 45 ~ "30–44",
      HHLDRAGEP < 65 ~ "45–64",
      TRUE ~ "65+"
    )
  )


# ---- Create weighted age-by-burden table ----
age_burden_table <- NASH_Renters %>%
  filter(!is.na(age_group), !is.na(burden_group)) %>%
  group_by(age_group, burden_group) %>%
  summarize(
    households = sum(WGTP, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(age_group) %>%
  mutate(
    pct = households / sum(households)
  ) %>%
  ungroup()


# ---- Nicely formatted table ----
age_burden_table %>%
  mutate(
    households = round(households),
    pct = scales::percent(pct, accuracy = 0.1)
  ) %>%
  arrange(age_group, burden_group) %>%
  kable(
    caption = "Housing Burden by Age of Householder",
    col.names = c("Age Group", "Burden Group", "Households", "Percent")
  ) %>%
  kable_styling(full_width = FALSE)


# ---- Create label text (SUPPRESS Extreme group) ----
age_burden_table <- age_burden_table %>%
  mutate(
    label_text = if_else(
      str_detect(burden_group, "^4"),
      "",  # ✅ suppress label for Extreme group
      paste0(
        burden_group,
        "<br>",
        scales::comma(round(households)), " households",
        "<br>",
        scales::percent(pct, accuracy = 0.1)
      )
    )
  )


# ---- Plotly stacked bar chart ----

age_burden_plot <- plot_ly(
  age_burden_table,
  x = ~age_group,
  y = ~households,
  color = ~burden_group,
  colors = burden_colors,
  type = "bar",
  text = ~label_text,
  hovertemplate = ~paste0(
    burden_group,
    "<br>Households: ", scales::comma(round(households)),
    "<br>Percent: ", scales::percent(pct, accuracy = 0.1),
    "<extra></extra>"
  )
) %>%
  layout(
    barmode = "stack",
    title = "Housing Burden by Age of Householder",
    xaxis = list(title = "Age Group"),
    yaxis = list(title = "Estimated Number of Households"),
    legend = list(traceorder = "normal")
  )

age_burden_plot

# =========================================================
# 20. Race-by-burden table and interactive visualization
# =========================================================

# ---- Create race categories ----
NASH_Renters <- NASH_Renters %>%
  mutate(
    race_group = case_when(
      HHLDRRAC1P == 1 ~ "White alone",
      HHLDRRAC1P == 2 ~ "Black or African American alone",
      is.na(HHLDRRAC1P) ~ "N/A",
      TRUE ~ "Other race"
    )
  )


# ---- Create weighted race-by-burden table ----
race_burden_table <- NASH_Renters %>%
  filter(!is.na(race_group), !is.na(burden_group)) %>%
  group_by(race_group, burden_group) %>%
  summarize(
    households = sum(WGTP, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(race_group) %>%
  mutate(
    pct = households / sum(households)
  ) %>%
  ungroup()


# ---- Nicely formatted table ----
race_burden_table %>%
  mutate(
    households = round(households),
    pct = scales::percent(pct, accuracy = 0.1)
  ) %>%
  arrange(race_group, burden_group) %>%
  kable(
    caption = "Housing Burden by Race of Householder",
    col.names = c("Race Group", "Burden Group", "Households", "Percent")
  ) %>%
  kable_styling(full_width = FALSE)


# ---- Create label text (suppress Extreme group labels) ----
race_burden_table <- race_burden_table %>%
  mutate(
    label_text = if_else(
      str_detect(burden_group, "^4"),
      "",
      paste0(
        burden_group,
        "<br>",
        scales::comma(round(households)), " households",
        "<br>",
        scales::percent(pct, accuracy = 0.1)
      )
    )
  )


# ---- Plotly stacked bar chart ----

race_burden_plot <- plot_ly(
  race_burden_table,
  x = ~race_group,
  y = ~households,
  color = ~burden_group,
  colors = burden_colors,
  type = "bar",
  text = ~label_text,
  hovertemplate = ~paste0(
    burden_group,
    "<br>Households: ", scales::comma(round(households)),
    "<br>Percent: ", scales::percent(pct, accuracy = 0.1),
    "<extra></extra>"
  )
) %>%
  layout(
    barmode = "stack",
    title = "Housing Burden by Race of Householder",
    xaxis = list(title = "Race of Householder"),
    yaxis = list(title = "Estimated Number of Households"),
    legend = list(traceorder = "normal")
  )

race_burden_plot

# =========================================================
# 21. PUMA-by-burden table and interactive visualization
# =========================================================

# ---- Create weighted PUMA-by-burden table ----
puma_burden_table <- NASH_Renters %>%
  filter(!is.na(PUMA), !is.na(burden_group)) %>%
  group_by(PUMA, burden_group) %>%
  summarize(
    households = sum(WGTP, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(PUMA) %>%
  mutate(
    pct = households / sum(households)
  ) %>%
  ungroup()


# ---- Nicely formatted table ----
puma_burden_table %>%
  mutate(
    households = round(households),
    pct = scales::percent(pct, accuracy = 0.1)
  ) %>%
  arrange(PUMA, burden_group) %>%
  kable(
    caption = "Housing Burden by Nashville PUMA",
    col.names = c("PUMA", "Burden Group", "Households", "Percent")
  ) %>%
  kable_styling(full_width = FALSE)


# ---- Create label text (suppress Extreme group labels) ----
puma_burden_table <- puma_burden_table %>%
  mutate(
    label_text = if_else(
      str_detect(burden_group, "^4"),
      "",
      paste0(
        burden_group,
        "<br>",
        scales::comma(round(households)), " households",
        "<br>",
        scales::percent(pct, accuracy = 0.1)
      )
    )
  )


# ---- Plotly stacked bar chart ----

puma_burden_plot <- plot_ly(
  puma_burden_table,
  x = ~factor(PUMA, levels = PUMA_list),
  y = ~households,
  color = ~burden_group,
  colors = burden_colors,
  type = "bar",
  text = ~label_text,
  hovertemplate = ~paste0(
    "PUMA: ", PUMA,
    "<br>Burden Group: ", burden_group,
    "<br>Households: ", scales::comma(round(households)),
    "<br>Percent: ", scales::percent(pct, accuracy = 0.1),
    "<extra></extra>"
  )
) %>%
  layout(
    barmode = "stack",
    title = "Housing Burden by Nashville PUMA",
    xaxis = list(
      title = "PUMA",
      categoryorder = "array",
      categoryarray = PUMA_list
    ),
    yaxis = list(
      title = "Estimated Number of Households"
    ),
    legend = list(traceorder = "normal")
  )

puma_burden_plot


# =========================================================
# 22. Map housing burden by PUMA
# =========================================================

# ---- Calculate burden statistics by PUMA ----

puma_map_data <- NASH_Renters %>%
  filter(!is.na(burden_group)) %>%
  group_by(PUMA) %>%
  summarize(
    
    total_households =
      sum(WGTP, na.rm = TRUE),
    
    cost_households =
      sum(
        WGTP[burden_group == "2 Cost burdened (30–49%)"],
        na.rm = TRUE
      ),
    
    severe_households =
      sum(
        WGTP[burden_group == "3 Severely burdened (50–99%)"],
        na.rm = TRUE
      ),
    
    pct_cost =
      cost_households / total_households,
    
    pct_severe =
      severe_households / total_households,
    
    pct_burdened =
      (cost_households + severe_households) /
      total_households,
    
    .groups = "drop"
  ) %>%
  mutate(
    PUMA = as.character(PUMA)
  )


# ---- Join burden statistics to PUMA geometry ----

PUMA_map_sf <- NASH_PUMAs %>%
  mutate(
    PUMACE20 = as.character(PUMACE20)
  ) %>%
  left_join(
    puma_map_data,
    by = c("PUMACE20" = "PUMA")
  )


# ---- Create color palette ----

burden_pal <- colorNumeric(
  palette = c(
    "#B3D7FF",
    "#66A3FF",
    "#0066CC",
    "#003366"
  ),
  domain = PUMA_map_sf$pct_burdened,
  na.color = "lightgray"
)


# ---- Create popup text ----

PUMA_map_sf <- PUMA_map_sf %>%
  mutate(
    popup_text = paste0(
      
      "<strong>PUMA:</strong> ",
      PUMACE20,
      
      "<br><br>",
      
      "<strong>Cost burdened households:</strong> ",
      scales::comma(round(cost_households)),
      
      "<br><strong>Percent cost burdened:</strong> ",
      scales::percent(pct_cost, accuracy = 0.1),
      
      "<br><br>",
      
      "<strong>Severely burdened households:</strong> ",
      scales::comma(round(severe_households)),
      
      "<br><strong>Percent severely burdened:</strong> ",
      scales::percent(pct_severe, accuracy = 0.1),
      
      "<br><br>",
      
      "<strong>Combined burden rate:</strong> ",
      scales::percent(pct_burdened, accuracy = 0.1)
    )
  )

# ---- Create interactive Leaflet map ----

PUMA_burden_map <- leaflet(PUMA_map_sf) %>%
  
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  addPolygons(
    fillColor = ~burden_pal(pct_burdened),
    fillOpacity = 0.6,
    color = "black",
    weight = 1,
    
    highlightOptions = highlightOptions(
      weight = 2,
      color = "white",
      bringToFront = TRUE
    ),
    
    popup = ~popup_text
  ) %>%
  
  addLegend(
    position = "bottomright",
    pal = burden_pal,
    values = ~pct_burdened,
    title = "Cost + Severe Burden (%)",
    labFormat = labelFormat(
      suffix = "%",
      transform = function(x) x * 100
    )
  )

PUMA_burden_map

Zillow Data Analysis Script:

# =============================================================================
# ZORI by ZIP for Nashville, TN
# Download -> Clean -> Reshape -> Summarize -> Visualize -> Map
# =============================================================================

# =============================================================================
# STEP 1. Set Parameters
# =============================================================================

CITY <- "Nashville"
STATE <- "TN"

MONTH_WINDOW <- 25

ZILLOW_URL <- paste0(
  "https://files.zillowstatic.com/research/public_csvs/",
  "zori/Zip_zori_uc_sfrcondomfr_sm_month.csv?t=1773853995"
)

# =============================================================================
# STEP 2. Load Required Libraries
# =============================================================================

suppressPackageStartupMessages({
  library(readr)
  library(dplyr)
  library(janitor)
  library(lubridate)
  library(tidyr)
  library(plotly)
  library(scales)
  library(gt)
  library(sf)
  library(leaflet)
  library(RColorBrewer)
  library(htmltools)
  library(tigris)
  library(stringr)
})

# =============================================================================
# STEP 3. Download the Zillow ZORI Dataset
# =============================================================================

local_file <- tempfile(fileext = ".csv")

download.file(
  url = ZILLOW_URL,
  destfile = local_file,
  mode = "wb"
)

# =============================================================================
# STEP 4. Import and Clean Variable Names
# =============================================================================

zori <- read_csv(
  local_file,
  show_col_types = FALSE
) |>
  clean_names()

# =============================================================================
# STEP 5. Keep Nashville, Tennessee ZIP Codes Only
# =============================================================================

zori_nashville <- zori |>
  filter(
    state == STATE,
    city == CITY
  )

# =============================================================================
# STEP 6. Convert Monthly Columns from Wide to Long Format
# =============================================================================

zori_nashville_long <- zori_nashville |>
  pivot_longer(
    cols = matches(
      "^\\d{1,2}/\\d{1,2}/\\d{4}$|^x?\\d{4}_\\d{2}_\\d{2}$"
    ),
    names_to = "date",
    values_to = "zori"
  ) |>
  mutate(
    date = sub("^x", "", date),
    
    date = ifelse(
      grepl("/", date),
      as.character(mdy(date)),
      as.character(
        ymd(gsub("_", "-", date))
      )
    ),
    
    date = as.Date(date),
    
    region_name = str_pad(
      as.character(region_name),
      width = 5,
      side = "left",
      pad = "0"
    )
  ) |>
  arrange(region_name, date)

# =============================================================================
# STEP 7. Determine Analysis Window
# =============================================================================

most_recent_date <- max(
  zori_nashville_long$date,
  na.rm = TRUE
)

cutoff_date <- most_recent_date %m-% months(MONTH_WINDOW)

# =============================================================================
# STEP 8. Keep Only Recent Observations
# =============================================================================

zori_nashville_window <- zori_nashville_long |>
  filter(
    date >= cutoff_date,
    !is.na(zori)
  )

# =============================================================================
# STEP 9. Create ZIP-Level Summary Statistics
# =============================================================================

zip_summary <- zori_nashville_window |>
  group_by(region_name) |>
  arrange(date) |>
  summarize(
    latest_date = max(date),
    latest_zori = dplyr::last(zori),
    first_zori = dplyr::first(zori),
    dollar_change = latest_zori - first_zori,
    pct_change = 100 * (latest_zori - first_zori) / first_zori,
    .groups = "drop"
  ) |>
  rename(zip = region_name)

# =============================================================================
# STEP 10. Create Endpoint Labels for ZIP Codes
# =============================================================================

zip_labels <- zori_nashville_window |>
  group_by(region_name) |>
  filter(date == max(date, na.rm = TRUE)) |>
  ungroup()

# =============================================================================
# STEP 11. Create Interactive Plotly Line Chart
# =============================================================================

ZORIplot <- plot_ly()

ZORIplot <- ZORIplot |>
  add_trace(
    data = zori_nashville_window,
    x = ~date,
    y = ~zori,
    split = ~region_name,
    type = "scatter",
    mode = "lines",
    hoverinfo = "text",
    showlegend = FALSE,
    text = ~paste(
      "ZIP:", region_name,
      "<br>Date:", format(date, "%Y-%m-%d"),
      "<br>ZORI:", dollar(zori)
    )
  )

ZORIplot <- ZORIplot |>
  add_trace(
    data = zip_labels,
    x = ~date,
    y = ~zori,
    type = "scatter",
    mode = "text",
    text = ~region_name,
    textposition = "middle right",
    hoverinfo = "none",
    showlegend = FALSE
  )

ZORIplot <- ZORIplot |>
  layout(
    title = paste0(
      "ZORI Rent Trends by ZIP — ",
      CITY,
      ", ",
      STATE
    ),
    xaxis = list(
      title = "Date"
    ),
    yaxis = list(
      title = "Typical Rent (ZORI Estimate, $)"
    ),
    showlegend = FALSE,
    margin = list(
      r = 120
    )
  )

# =============================================================================
# STEP 12. Display Interactive Plotly Chart
# =============================================================================

ZORIplot

# =============================================================================
# STEP 13. Download ZIP Code Boundaries
# =============================================================================

options(tigris_use_cache = TRUE)

zip_shapes <- tigris::zctas(
  year = 2020,
  cb = TRUE
) |>
  st_transform(4326)

# =============================================================================
# STEP 14. Standardize ZIP Codes
# =============================================================================

zip_shapes <- zip_shapes |>
  mutate(
    zip = str_pad(
      as.character(ZCTA5CE20),
      width = 5,
      side = "left",
      pad = "0"
    )
  )

# =============================================================================
# STEP 15. Verify ZIP Matches
# =============================================================================

unmatched_zips <- anti_join(
  zip_summary,
  st_drop_geometry(zip_shapes),
  by = "zip"
)

print(unmatched_zips)

# =============================================================================
# STEP 16. Join ZORI Data to ZIP Boundaries
# =============================================================================

zip_map <- zip_shapes |>
  filter(zip %in% zip_summary$zip) |>
  left_join(
    zip_summary,
    by = "zip"
  )

cat(
  "Mapped",
  nrow(zip_map),
  "ZIP codes out of",
  nrow(zip_summary),
  "ZIP codes in Zillow data.\n"
)

# =============================================================================
# STEP 17. Build Popup Content
# =============================================================================

zip_map <- zip_map |>
  mutate(
    popup = paste0(
      "<b>ZIP Code:</b> ", zip,
      "<br><b>Latest ZORI:</b> ",
      dollar(latest_zori),
      "<br><b>Change (",
      MONTH_WINDOW,
      " months):</b> ",
      dollar(dollar_change),
      "<br><b>Percent Change:</b> ",
      round(pct_change, 1),
      "%",
      "<br><b>Latest Month:</b> ",
      format(latest_date, "%B %Y")
    )
  )

# =============================================================================
# STEP 18. Create Color Scale
# =============================================================================

pal <- colorNumeric(
  palette = colorRampPalette(
    c(
      "#B3D7FF",  # light blue
      "#80BFFF",
      "#4D94DB",
      "#1F5FAF",
      "#003366"   # dark navy blue
    )
  )(100),
  domain = zip_map$latest_zori,
  na.color = "#DDDDDD"
)

# =============================================================================
# STEP 19. Build Interactive Leaflet Map
# =============================================================================

ZORImap <- leaflet(zip_map) |>
  
  addProviderTiles(
    providers$CartoDB.Positron
  ) |>
  
  addPolygons(
    fillColor = ~pal(latest_zori),
    fillOpacity = 0.75,
    color = "#777777",
    weight = 1,
    smoothFactor = 0.3,
    popup = ~popup,
    
    highlightOptions = highlightOptions(
      weight = 3,
      color = "black",
      bringToFront = TRUE
    )
  ) |>
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~latest_zori,
    title = "Latest ZORI ($)",
    opacity = 1
  )

# =============================================================================
# STEP 20. Display Interactive Leaflet Map
# =============================================================================

ZORImap