The Nashville area and surrounding counties to Davidson County have seen some statistically significant increases since 2019 in people who speak a language other than English at home. A district in Smyrna had the largest increase, with a 16.3 percentage point difference from 2019 to 2024.
The map below shows the breakdown of percentage change in various districts in the mid-state, with blue showing fewer people who speak a language other than English at home and red showing the opposite. A graph shows any statistically significant changes in districts.
# ----------------------------------------------------------
# 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("e96c46601fca77f6f3fcb0f72b673a75aed0ff2a", 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. with broadband" # <<< 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)