Change in Non-English Speaking Households around Nashville, TN

The percentage change of people, 5 years or older, around the Nashville area hasn’t seen that drastic of an increase, contrary to popular belief. The most significant increase has been 15% in Robertson County, while the least significant change has been around 6% in Rutherford County.

Change in pct. who speak non-English at home, 2019-2024

Code:

# ----------------------------------------------------------
# Step 1: Install required packages (if missing)
# ----------------------------------------------------------

if (!require("tidyverse"))
  install.packages("tidyverse")
if (!require("tidycensus"))
  install.packages("tidycensus")
if (!require("sf"))
  install.packages("sf")
if (!require("leaflet"))
  install.packages("leaflet")
if (!require("htmlwidgets"))
  install.packages("htmlwidgets")
if (!require("plotly"))
  install.packages("plotly")

# ----------------------------------------------------------
# Step 2: Load libraries
# ----------------------------------------------------------

library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

# ----------------------------------------------------------
# Step 3: Transmit Census API key
# ----------------------------------------------------------

census_api_key

# ----------------------------------------------------------
# Step 4: Define analysis years, variables, labels, and codebooks
# ----------------------------------------------------------

year_old <- 2019  # Avoid overlapping ACS datasets
year_new <- 2024  # Avoid overlapping ACS datasets

# ---- Load codebooks for both years (for inspection / validation) ----

# Detailed tables
Old_DetailedTables <- load_variables(year_old, "acs5", cache = TRUE)
New_DetailedTables <- load_variables(year_new, "acs5", cache = TRUE)

# Subject tables
Old_SubjectTables  <- load_variables(year_old, "acs5/subject", cache = TRUE)
New_SubjectTables  <- load_variables(year_new, "acs5/subject", cache = TRUE)

# Profile tables
Old_ProfileTables  <- load_variables(year_old, "acs5/profile", cache = TRUE)
New_ProfileTables  <- load_variables(year_new, "acs5/profile", cache = TRUE)

# ---- Enter variable name(s); names often change year to year ----
Variable_old <- c(Estimate = "DP02_0113P")
Variable_new <- c(Estimate = "DP02_0114P")

# Human-readable label for variable
target_label <- "Pct. speaking non-English at home" # <<< Edit as needed

# ----------------------------------------------------------
# Step 5: Fetch ACS data – newer year
# ----------------------------------------------------------

mydata_new <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = Variable_new,
  year = year_new,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
) %>%
  rename(Area = NAME) %>%
  transmute(
    GEOID,
    Area,
    geometry,
    Estimate_new = EstimateE,
    MoE_new      = EstimateM
  )

# ----------------------------------------------------------
# Step 6: Fetch ACS data – older year
# ----------------------------------------------------------

mydata_old <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = Variable_old,
  year = year_old,
  survey = "acs5",
  output = "wide",
  geometry = FALSE
) %>%
  rename(Area = NAME) %>%
  transmute(
    GEOID,
    Area,
    Estimate_old = EstimateE,
    MoE_old      = EstimateM
  )

# ----------------------------------------------------------
# Step 7: Join datasets
# ----------------------------------------------------------

mydata_joined <- mydata_new %>%
  left_join(mydata_old, by = c("GEOID", "Area"))

# ----------------------------------------------------------
# Step 8: Flexible filtering by Area
# ----------------------------------------------------------

search_terms <- c("Rutherford", "Davidson", "Williamson", "Robertson", "Sumner", "Wilson", "Cheatham")
# Use c("ALL") to include all areas

if (length(search_terms) == 1 && toupper(search_terms) == "ALL") {
  filtereddata <- mydata_joined
} else {
  or_pattern <- paste(
    stringr::str_replace_all(search_terms, "([\\W])", "\\\\\\1"),
    collapse = "|"
  )
  
  filtereddata <- mydata_joined %>%
    filter(
      stringr::str_detect(
        Area %>% replace_na(""),
        regex(or_pattern, ignore_case = TRUE)
      )
    )
}

# ----------------------------------------------------------
# Step 9: Compute change and statistical significance
# ----------------------------------------------------------

filtereddata <- filtereddata %>%
  mutate(
    Diff     = Estimate_new - Estimate_old,
    Diff_MOE = tidycensus::moe_sum(MoE_new, MoE_old),
    Sig_90   = tidycensus::significance(
      Estimate_new, MoE_new,
      Estimate_old, MoE_old
    ),
    SigDirection = case_when(
      Sig_90 & Diff > 0 ~ "Significant increase",
      Sig_90 & Diff < 0 ~ "Significant decrease",
      TRUE             ~ "No statistically significant change"
    )
  )

# ----------------------------------------------------------
# Step 10: Prepare spatial data
# ----------------------------------------------------------

mapdata <- filtereddata %>%
  st_as_sf() %>%
  st_transform(4326)

# ----------------------------------------------------------
# Step 11: Ranked change dot plot (significance encoded by color)
# ----------------------------------------------------------

plotdf <- filtereddata %>%
  st_drop_geometry() %>%
  mutate(
    y_ordered = reorder(Area, Diff),
    sig_class = case_when(
      SigDirection == "Significant increase" ~ "Increase (significant)",
      SigDirection == "Significant decrease" ~ "Decrease (significant)",
      TRUE ~ "Not statistically significant"
    ),
    hover_text = paste0(
      "Area: ", Area,
      "<br>Change (", year_old, " to ", year_new, "): ", scales::comma(Diff),
      "<br>Significance: ", sig_class,
      "<br>", year_old, ": ", scales::comma(Estimate_old),
      "<br>", year_new, ": ", scales::comma(Estimate_new)
    )
  )

mygraph <- plotly::plot_ly(
  data = plotdf,
  x = ~Diff,
  y = ~as.character(y_ordered),
  type = "scatter",
  mode = "markers",
  color = ~sig_class,
  colors = c(
    "Increase (significant)" = "#b2182b",
    "Decrease (significant)" = "#2166ac",
    "Not statistically significant" = "#bdbdbd"
  ),
  marker = list(size = 9),
  text = ~hover_text,
  hovertemplate = "%{text}<extra></extra>"
) %>%
  layout(
    title = list(
      text = paste0(
        "Estimated Change in ", target_label,
        " (", year_old, " to ", year_new, ")",
        "<br><sup>County subdivisions. Color indicates statistical significance (90%).</sup>"
      )
    ),
    xaxis = list(title = "Estimated change", tickformat = ",.0f"),
    yaxis = list(title = "", automargin = TRUE)
  )

mygraph

# ----------------------------------------------------------
# Step 12: Leaflet popup content (no MOE for change)
# ----------------------------------------------------------

mapdata$popup <- paste0(
  "<strong>", mapdata$Area, "</strong><br/>",
  "<hr>",
  year_old, ": ", format(mapdata$Estimate_old, big.mark = ","),
  " (plus/minus ", format(mapdata$MoE_old, big.mark = ","), ")<br/>",
  year_new, ": ", format(mapdata$Estimate_new, big.mark = ","),
  " (plus/minus ", format(mapdata$MoE_new, big.mark = ","), ")<br/>",
  "<b>Change:</b> ", format(mapdata$Diff, big.mark = ","), "<br/>",
  "<em>", mapdata$SigDirection, " (90%)</em>"
)

# ----------------------------------------------------------
# Step 13: Leaflet map of change
# ----------------------------------------------------------

max_abs <- max(abs(mapdata$Diff), na.rm = TRUE)
bins <- pretty(c(-max_abs, max_abs), n = 6)

pal <- colorBin(
  palette = "RdBu",
  domain  = mapdata$Diff,
  bins    = bins,
  reverse = TRUE
)

DivisionMap <- leaflet(mapdata) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor   = ~pal(Diff),
    fillOpacity = 0.6,
    color       = "black",
    weight      = 1,
    popup       = ~popup
  ) %>%
  addLegend(
    pal    = pal,
    values = ~Diff,
    title  = paste0(
      "Change (", year_old, " to ", year_new, ")",
      "<br/><small>", target_label, "</small>"
    ),
    labFormat = labelFormat(big.mark = ",")
  )

DivisionMap

# ----------------------------------------------------------
# Step 14: Export outputs
# ----------------------------------------------------------

saveWidget(as_widget(mygraph), "ACSGraph_Diff.html", selfcontained = TRUE)
saveWidget(DivisionMap, "ACSMap_Diff.html", selfcontained = TRUE)