PUMS Housing Affordability Analysis

TidyCensus reference:

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

See, particularly, chapters 9 and 10.

Rent burden in the Nashville area

About 40 percent of Nashville-area renter households struggle to pay their leases, the latest Census Bureau data show.

In some neighborhoods, like North Nashville and Antioch, the figure ranges closer to half.

---

---

The map shows the proportion of renter households spending 35 percent or more of their monthly income on a lease in each of 11 Census-defined areas across Davidson, Rutherford and Williamson counties. Financial experts generally say rent should consume no more than 30 percent of a household budget. The 35-percent or more category is the highest rent-burden category available in published summaries of the Census Bureau’s annual American Community Survey results. The figures on the map come from the single-year 2024 ACS, the most recent ACS that the bureau has released.

This figure shows each area’s estimate along with the estimate’s error margin. The error margins for all 11 areas overlap, meaning that random sampling error could account for the difference between any pair of areas. But the error margins for Antioch and North Nashville are the only two that fall wholly above the 40 percent mark.

---

---

Some types of renter households might be especially likely to struggle financially. For example, couples who rent - whether married or cohabitating - can sometimes afford their leases more often than single renters who live alone, especially among renter households headed by someone who is white.

Once again, though, random sampling error could explain some of the differences. In the second of the two figures below, groups of renters with non-overlapping error margins have statistically significant differences. For example, renter households headed by someone who is white and who is living with at least one other person, whether a spouse, a partner or someone else, are significantly less likely to be rent burdened that renter households headed by someone who is Black and living alone or without a spouse or partner.

---

---

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:

---

---

Combining ZIP codes and PUMAs

It might be helpful to overlay the ZORI map with the PUMA boundaries:

---

---

ACS / 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")
if (!require("survey")) install.packages("survey")
if (!require("srvyr")) install.packages("srvyr")
library(tidycensus)
library(tidyverse)
library(tigris)
library(sf)
library(leaflet)
library(kableExtra)
library(plotly)
library(survey)
library(srvyr)


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

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

# ============================================================
# 3. LOAD ACS CODEBOOKS
# ============================================================

# DetailedTables <- load_variables(2024, "acs1")
# SubjectTables  <- load_variables(2024, "acs1/subject")
ProfileTables  <- load_variables(2024, "acs1/profile")

# ============================================================
# 4. DEFINE VARIABLES OF INTEREST
# ============================================================

VariableList =
  c(
    UnitCount = "DP04_0002",
    RentCount = "DP04_0136",
    Rent35Count = "DP04_0142",
    Rent35Pct = "DP04_0142P"
  )

# =========================================================
# 5. Pull renter data
# =========================================================

CountyLevelData <- get_acs(
  geography = "county",
  state = "TN",
  variables = VariableList,
  year = 2024,
  survey = "acs1",
  output = "wide",
  geometry = TRUE
)

mydata <- get_acs(
  geography = "public use microdata area",
  state = "TN",
  variables = VariableList,
  year = 2024,
  survey = "acs1",
  output = "wide",
  geometry = TRUE
)

# ============================================================
# 6. FILTER FOR SELECTED COUNTIES & PROJECT TO WGS84
# ============================================================
counties <- c("Davidson County",
              "Rutherford County",
              "Williamson County",
              "Murfreesboro city")

CountyData <- mydata %>%
  filter(str_detect(NAME, paste(counties, collapse = "|"))) %>%
  st_transform(4326)

# ============================================================
# 7. CREATE POPUP TEXT
# ============================================================

CountyData <- CountyData %>%
  mutate(
    popup_text = paste0(
      "<strong>", NAME, "</strong><br><br>",
      
      "Total occupied housing units: ",
      format(UnitCountE, big.mark = ","), "<br>",
      
      "Renter-occupied units: ",
      format(RentCountE, big.mark = ","), "<br>",
      
      "Renters paying 35%+ of income for rent: ",
      format(Rent35CountE, big.mark = ","), "<br>",
      
      "Percent paying 35%+ of income for rent: ",
      round(Rent35CountPE, 1), "%<br>",
      
      "Margin of error: ±",
      round(Rent35CountPM, 1),
      " percentage points"
    )
  )

# ============================================================
# 8. CREATE COLOR PALETTE
# ============================================================

pal <- colorNumeric(
  palette = "Blues",
  domain = CountyData$Rent35CountPE,
  na.color = "#D3D3D3"
)

# ============================================================
# 9. CREATE INTERACTIVE MAP
# ============================================================

Map <- leaflet(CountyData) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  addPolygons(
    fillColor = ~pal(Rent35CountPE),
    fillOpacity = 0.75,
    color = "black",
    weight = 1,
    opacity = 1,
    popup = ~popup_text,
    highlightOptions = highlightOptions(
      weight = 3,
      color = "yellow",
      bringToFront = TRUE
    )
  ) %>%
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~Rent35CountPE,
    title = "% of Renters Paying<br>35%+ of Income for Rent",
    opacity = 0.8
  )

Map

# ============================================================
# 10. INTERACTIVE PLOTLY GRAPH OF ESTIMATES WITH ERROR BARS
# ============================================================

PlotData <- CountyData %>%
  st_drop_geometry() %>%
  mutate(ShortName = str_remove(NAME, " PUMA; Tennessee"))

mygraph <- plot_ly(
  data = PlotData,
  x = ~ Rent35CountPE,
  y = ~ reorder(ShortName, Rent35CountPE),
  type = "scatter",
  mode = "markers",
  
  error_x = list(
    type = "data",
    array = ~ Rent35CountPM,
    visible = TRUE,
    color = "black",
    thickness = 1.5,
    width = 4
  ),
  
  marker = list(
    color = "#2171B5",
    # medium blue from the Blues palette
    size = 10
  ),
  
  hovertemplate = paste(
    "<b>%{y}</b><br>",
    "Estimate: %{x:.1f}%<br>",
    "MOE: ±%{customdata:.1f}<br>",
    "<extra></extra>"
  ),
  
  customdata = ~ Rent35CountPM,
  
  showlegend = FALSE
) %>%
  
  layout(
    title = "Renters Paying 35%+ of Income for Rent",
    xaxis = list(title = "Percent of Renters", ticksuffix = "%"),
    yaxis = list(title = "", automargin = TRUE)
  )

mygraph

# =========================================================
# 11. Get PUMS Codebooks
# =========================================================

data(pums_variables)

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


# =========================================================
# 12. Pull renter data
# =========================================================

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

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

# =========================================================
# 13. Filter Nashville renters
# =========================================================

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

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

# =========================================================
# 14. Rent burden by race and household type
# =========================================================

NASH_hh_recoded <- NASH_Renters %>%
  mutate(
    race_ethnicity = case_when(
      HISP != "01" ~ "Hispanic",
      HISP == "01" & RAC1P == "1" ~ "White",
      HISP == "01" & RAC1P == "2" ~ "Black",
      TRUE ~ "Other"),
    
    HHT2Type = case_when(
      HHT2 %in% c("01", "02", "03", "04") ~ "Couple renting",
      HHT2 %in% c("05", "09") ~ "Lone renter",
      TRUE ~ "Other renter"
    )
  )

NASH_hh_summary <- NASH_hh_recoded %>%
  filter(race_ethnicity != "Other") %>%
  group_by(race_ethnicity, HHT2Type) %>%
  summarize(
    prop_burdened = sum(WGTP[GRPIP >= 35]) / sum(WGTP)
  )

# =========================================================
# 15. DISPLAY SUMMARY TABLE
# =========================================================

NASH_hh_table <- NASH_hh_summary %>%
  mutate(
    PercentBurdened = round(100 * prop_burdened, 1)
  ) %>%
  select(
    race_ethnicity,
    HHT2Type,
    PercentBurdened
  )

Table1 <- NASH_hh_table %>%
  kbl(
    col.names = c(
      "Race / Ethnicity",
      "Household Type",
      "Percent Burdened"
    ),
    caption = "Percent of Renters Paying 35%+ of Income for Housing"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  )

Table1

# =========================================================
# 16. DISPLAY GRAPHIC
# =========================================================

PlotData <- NASH_hh_summary %>%
  mutate(
    PercentBurdened = 100 * prop_burdened
  )

# Blues palette matching the map

BluePalette <- c(
  "#08306B",  # dark blue
  "#4292C6",  # medium blue
  "#C6DBEF"   # light blue
)

hh_graph <- plot_ly(
  data = PlotData,
  x = ~race_ethnicity,
  y = ~PercentBurdened,
  color = ~HHT2Type,
  colors = BluePalette,
  type = "bar",
  
  hovertemplate = paste(
    "<b>%{x}</b><br>",
    "Household Type: %{fullData.name}<br>",
    "Percent Burdened: %{y:.1f}%<br>",
    "<extra></extra>"
  )
) %>%
  layout(
    barmode = "group",
    title = "Rent Burden by Race/Ethnicity and Household Type",
    yaxis = list(
      title = "Percent Paying 35%+ of Income for Housing",
      ticksuffix = "%"
    ),
    xaxis = list(
      title = ""
    )
  )

hh_graph

# =========================================================
# 17. Pull renter data again, with replicate weights
# =========================================================

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

tn_hh_replicates <- get_pums(
  variables = vars,
  state = "TN",
  survey = "acs1",
  variables_filter = list(SPORDER = 1),
  year = 2024,
  recode = TRUE,
  rep_weights = "housing"
)

# =========================================================
# 18. Apply replicate weights
# =========================================================

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

NASH_hh_svy <- tn_hh_replicates %>%
  to_survey(type = "housing", 
            design = "rep_weights") %>%
  filter(TEN == "3") %>% 
  filter(PUMA %in% PUMA_list)

class(NASH_hh_svy)

# =========================================================
# 19. RECODE VARIABLES IN SURVEY OBJECT
# =========================================================

NASH_hh_svy_recoded <- NASH_hh_svy %>%
  
  mutate(
    race_ethnicity = case_when(
      HISP != "01" ~ "Hispanic",
      HISP == "01" & RAC1P == "1" ~ "White",
      HISP == "01" & RAC1P == "2" ~ "Black",
      TRUE ~ "Other"
    ),
    
    HHT2Type = case_when(
      HHT2 %in% c("01", "02", "03", "04") ~ "Couple renting",
      HHT2 %in% c("05", "09") ~ "Lone renter",
      TRUE ~ "Other renter"
    ),
    
    Burdened = GRPIP >= 35
  ) %>%
  
  filter(race_ethnicity != "Other")

# =========================================================
# 20. ESTIMATE RENT BURDEN WITH MOEs
# =========================================================

NASH_hh_MOE <- NASH_hh_svy_recoded %>%
  
  group_by(race_ethnicity, HHT2Type) %>%
  
  summarize(
    prop_burdened = survey_mean(
      Burdened,
      vartype = "se"
    )
  ) %>%
  
  mutate(
    PercentBurdened = 100 * prop_burdened,
    
    MOE90 = 1.645 * prop_burdened_se,
    
    MOE90Percent = 100 * MOE90,
    
    LowerCI = PercentBurdened - MOE90Percent,
    
    UpperCI = PercentBurdened + MOE90Percent
  )

# =========================================================
# 21. DISPLAY TABLE WITH MOEs
# =========================================================

Table2 <- NASH_hh_MOE %>%
  
  transmute(
    race_ethnicity,
    HHT2Type,
    Estimate = round(PercentBurdened, 1),
    MOE = round(MOE90Percent, 1)
  ) %>%
  
  kbl(
    col.names = c(
      "Race / Ethnicity",
      "Household Type",
      "Estimate (%)",
      "MOE (±)"
    ),
    caption = "Rent Burden Estimates with ACS 90% MOEs"
  ) %>%
  
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  )

Table2

# =========================================================
# 22. PLOTLY CHART WITH ERROR BARS
# =========================================================

BluePalette <- c(
  "#08306B",
  "#4292C6",
  "#C6DBEF"
)

hh_graph_moe <- plot_ly(
  data = NASH_hh_MOE,
  
  x = ~PercentBurdened,
  
  y = ~paste(race_ethnicity, HHT2Type, sep = " | "),
  
  color = ~HHT2Type,
  
  colors = BluePalette,
  
  type = "scatter",
  
  mode = "markers",
  
  error_x = list(
    type = "data",
    array = ~MOE90Percent,
    visible = TRUE,
    color = "black",
    thickness = 1.5,
    width = 4
  ),
  
  marker = list(size = 10),
  
  hovertemplate = paste(
    "<b>%{y}</b><br>",
    "Estimate: %{x:.1f}%<br>",
    "MOE: ±%{customdata:.1f}%<br>",
    "<extra></extra>"
  ),
  
  customdata = ~MOE90Percent
) %>%
  
  layout(
    title =
      "Rent Burden by Race/Ethnicity and Household Type",
    
    xaxis = list(
      title = "Percent Paying 35%+ of Income for Housing",
      ticksuffix = "%"
    ),
    
    yaxis = list(
      title = "",
      automargin = TRUE
    )
  )

hh_graph_moe

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 18A. Download Nashville-Area PUMA Boundaries
# =============================================================================

PUMA_shapes <- get_acs(
  geography = "public use microdata area",
  state = "TN",
  variables = "DP04_0002",
  year = 2024,
  geometry = TRUE
) |>
  st_transform(4326)

PUMA_shapes <- PUMA_shapes |>
  filter(
    GEOID %in% c(
      "4702401",
      "4702402",
      "4702403",
      "4702404",
      "4702405",
      "4702406"
    )
  )

# =============================================================================
# STEP 18B. Create PUMA Popups
# =============================================================================

PUMA_shapes <- PUMA_shapes |>
  mutate(
    popup_text =
      paste0(
        "<strong>",
        NAME,
        "</strong><br>",
        "PUMA Code: ",
        substr(GEOID, 3, 7)
      )
  )

# =============================================================================
# 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 19A. Build ZIP Map with PUMA Overlay
# =============================================================================

ZORImap_PUMA <- leaflet() |>
  
  addProviderTiles(
    providers$CartoDB.Positron
  ) |>
  
  addPolygons(
    data = zip_map,
    fillColor = ~pal(latest_zori),
    fillOpacity = 0.75,
    color = "#777777",
    weight = 1,
    smoothFactor = 0.3,
    popup = ~popup,
    
    highlightOptions = highlightOptions(
      weight = 3,
      color = "black",
      bringToFront = TRUE
    )
  ) |>
  
  addPolygons(
    data = PUMA_shapes,
    fill = FALSE,
    color = "#B22222",
    weight = 5,
    opacity = 1,
    dashArray = "8,6",
    popup = ~popup_text
  ) |>
  
  addLegend(
    position = "bottomright",
    pal = pal,
    values = zip_map$latest_zori,
    title = "Latest ZORI ($)",
    opacity = 1
  )

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

ZORImap

# =============================================================================
# STEP 20A. Display ZIP + PUMA Overlay Map
# =============================================================================

ZORImap_PUMA