This map and chart show how the percentage of people aged 5 and older who speak a language other than English at home has changed between 2019 and 2024 across county subdivisions in and around Davidson County, Tennessee. The map highlights both increases and decreases, while the chart focuses on increases and indicates whether those changes are statistically significant at the 90% confidence level.

Map

Graphic

Code

# Load packages
library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

# Census API key
census_api_key("04d318cbd6d7417e56fec7f3873d8302c28fe21d", install = FALSE)

# Years + variables
year_old <- 2019
year_new <- 2024

Variable_old <- c(Estimate = "DP02_0113P")
Variable_new <- c(Estimate = "DP02_0114P")

# Fetch data
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)

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)

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

# Filter counties
search_terms <- c("Davidson", "Rutherford", "Williamson", 
                  "Cheatham", "Robertson", "Sumner", "Wilson")

or_pattern <- paste(search_terms, collapse = "|")

filtereddata <- mydata_joined %>%
  filter(str_detect(Area, regex(or_pattern, ignore_case = TRUE)))

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

# Map data
mapdata <- st_as_sf(filtereddata) %>%
  st_transform(4326)

# Popup
mapdata$popup <- paste0(
  "<strong>", mapdata$Area, "</strong><br/>",
  "2019: ", round(mapdata$Estimate_old, 1), "%<br/>",
  "2024: ", round(mapdata$Estimate_new, 1), "%<br/>",
  "<b>Change:</b> ", round(mapdata$Diff, 1), "%<br/>",
  "<em>", mapdata$SigDirection, "</em>"
)

# Map
bins <- c(-20, -15, -10, -5, 0, 5, 10, 15, 20)

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

DivisionMap <- leaflet(mapdata) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor = ~pal(Diff),
    fillOpacity = 0.7,
    color = "black",
    weight = 1,
    popup = ~popup
  ) %>%
  addLegend(
    pal = pal,
    values = ~Diff,
    title = "Change (2019–2024)"
  )

# Plot
plotdf <- filtereddata %>%
  st_drop_geometry() %>%
  filter(Diff >= 0) %>%
  mutate(
    y_ordered = reorder(Area, Diff),
    sig_class = ifelse(SigDirection == "Significant increase",
                       "Significant", "Not significant")
  )

mygraph <- plot_ly(
  data = plotdf,
  x = ~Diff,
  y = ~as.character(y_ordered),
  type = "scatter",
  mode = "markers",
  color = ~sig_class
)