Region table and statewide map

Estimates by Region
(2019–2023 ACS 5-Year ACS)
Region Percent (%) Margin of Error Count
Davidson 47.3 0.7 493,272
Doughnut 35.4 1.4 715,926
Non-doughnut 16.9 2.3 179,896

Percent with a bachelor’s degree or higher
Tennessee, 2019-2023. Source: American Community Survey.

Code

# ============================================================
# 0. INSTALL AND LOAD REQUIRED PACKAGES
# ============================================================

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("leaflet.extras2")) install.packages("leaflet.extras2")
if (!require("gt")) install.packages("gt")
if (!require("gtExtras")) install.packages("gtExtras")
if (!require("plotly")) install.packages("plotly")

library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(leaflet.extras2)
library(gt)
library(gtExtras)
library(plotly)

# ============================================================
# 1. CENSUS API KEY
# ============================================================

# census_api_key("_______")

# ============================================================
# 2. LOAD ACS CODEBOOKS
# ============================================================

DetailedTables <- load_variables(2023, "acs5", cache = TRUE)
SubjectTables  <- load_variables(2023, "acs5/subject", cache = TRUE)
ProfileTables  <- load_variables(2023, "acs5/profile", cache = TRUE)

# ============================================================
# 3. DEFINE VARIABLES OF INTEREST
# ============================================================

VariableList =
  c(
    Count_   = "DP02_0059",
    Percent_ = "DP02_0068P"
  )

# ============================================================
# 4. FETCH COUNTY SUBDIVISION DATA (TENNESSEE)
# ============================================================

mydata <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = VariableList,
  year = 2023,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
)

# ============================================================
# 5. CLEAN AND REFORMAT GEOGRAPHIC NAMES
# ============================================================

mydata <- mydata %>%
  separate_wider_delim(
    NAME,
    delim  = ", ",
    names  = c("Division", "County", "State")
  ) %>%
  mutate(County = str_remove(County, " County"))

# ============================================================
# NOTE: See Practical R parts 1 & 2 for sections 6 through 11
# of this script
# ============================================================

# ============================================================
# 12. REGION-LEVEL AGGREGATION AND TABLE
# ============================================================
CountyLevelData <- get_acs(
  geography = "county",
  state = "TN",
  variables = VariableList,
  year = 2023,
  survey = "acs5",
  output = "wide",
  geometry = FALSE
)

SelectedCounties <- CountyLevelData %>% 
  filter(NAME %in% c(
    "Cannon County, Tennessee","Cheatham County, Tennessee",
    "Davidson County, Tennessee","Dickson County, Tennessee",
    "Hickman County, Tennessee","Macon County, Tennessee",
    "Maury County, Tennessee","Robertson County, Tennessee",
    "Rutherford County, Tennessee","Smith County, Tennessee",
    "Sumner County, Tennessee","Trousdale County, Tennessee",
    "Williamson County, Tennessee","Wilson County, Tennessee")
  ) %>%
  select(NAME, Percent_E, Percent_M, Count_E) %>%
  rename(
    County  = NAME,
    Percent = Percent_E,
    Range   = Percent_M,
    Count   = Count_E
  ) %>%
  arrange(desc(Percent))

SelectedCounties <- SelectedCounties %>%
  mutate(
    CountyClean = str_remove(County, ", Tennessee"),
    Region = case_when(
      CountyClean == "Davidson County" ~ "Davidson",
      CountyClean %in% c(
        "Rutherford County","Williamson County","Cheatham County",
        "Robertson County","Sumner County","Wilson County"
      ) ~ "Doughnut",
      TRUE ~ "Non-doughnut"
    )
  )

RegionTableData <- SelectedCounties %>%
  group_by(Region) %>%
  summarise(
    Percent = mean(Percent, na.rm = TRUE),
    Range   = mean(Range, na.rm = TRUE),
    Count   = sum(Count, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(Percent))

RegionTable <- RegionTableData %>%
  gt() %>%
  gt_theme_espn() %>%
  cols_label(
    Percent = "Percent (%)",
    Range   = "Margin of Error",
    Count   = "Count"
  ) %>%
  fmt_number(columns = c(Percent, Range), decimals = 1) %>%
  fmt_number(columns = c(Count), decimals = 0, use_seps = TRUE) %>%
  tab_header(
    title = "Estimates by Region",
    subtitle = "(2019–2023 ACS 5-Year ACS)"
  )

RegionTable

# ============================================================
# 13. MAP ENTIRE STATE
# ============================================================

mapdata3 <- mydata %>%
  rename(
    Percent   = Percent_E,
    PercentEM = Percent_M,
    Count     = Count_E,
    CountEM   = Count_M
  ) %>%
  st_as_sf() %>%
  st_transform(4326)

# Get Tennessee county boundaries for overlay
county_borders <- get_acs(
  geography = "county",
  state = "TN",
  variables = "B01001_001",  # dummy population variable
  year = 2023,
  geometry = TRUE
) %>%
  st_as_sf() %>%
  st_transform(4326)

# Viridis color palette
pal <- colorNumeric(
  palette = "viridis",
  domain = mapdata3$Percent
)

# Build leaflet map with county borders
DivisionMap3 <- leaflet(mapdata3) %>%
  # Base maps
  addProviderTiles("CartoDB.Positron", group = "Positron (Light)") %>%
  addProviderTiles("Esri.WorldStreetMap", group = "ESRI Streets") %>%
  addProviderTiles("Esri.WorldImagery", group = "ESRI Satellite") %>%
  # Main polygons
  addPolygons(
    fillColor = ~pal(Percent),
    fillOpacity = 0.7,
    color = "#333333",
    weight = 1,
    group = "Data Layer",
    popup = ~paste0(
      "<b>", Division, "</b><br>",
      "County: ", County, "<br>",
      "Percent: ", Percent, "%<br>",
      "Percent MOE: ±", PercentEM, "<br>",
      "Count: ", Count, "<br>",
      "Count MOE: ±", CountEM
    )
  ) %>%
  # County borders overlay
  addPolylines(
    data = county_borders,
    color = "black",
    weight = 1,
    opacity = 0.8,
    group = "County Borders"
  ) %>%
  # Legend
  addLegend(
    pal = pal,
    values = mapdata3$Percent,
    position = "topright",
    title = "Percent"
  ) %>%
  addLayersControl(
    baseGroups = c(
      "Positron (Light)",
      "ESRI Streets",
      "ESRI Satellite"
    ),
    overlayGroups = c("Data Layer", "County Borders"),
    options = layersControlOptions(
      collapsed = TRUE,
      position = "bottomleft"   # ?? Move control here
    )
  )

DivisionMap3