Should the Davidson County Doughnut Expand its Resources for Non-English Speakers?

Data from the American Community Survey has shown how over the past 5 years, more counties than not have shown an increase in non-English speakers. These individuals do also speak English, but it may not be the most accurate, this data set looks for those who mainly speak languages other than English at home. While less counties have a significant increase than those that do, the increase can still be seen.

The question of if more resources should be allocated to help these bilingual individuals has been long debated between activists and lawmakers. If these numbers continue to follow the trend that has been seen over the past 5 years, these resource allocations will become necessary.

Change in Percentage of People 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("b2b546278d54eb1ced62b21ec507abf728af01c2", install = FALSE)

# ----------------------------------------------------------
# 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", "Cheatham", "Robertson", "Sumner", "Wilson")
# 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)