Practical R Part 2

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

---

Estimates by County
(2019–2023 ACS 5-Year ACS)
County	Percent (%)	Percent MOE	Count	Count MOE
Williamson	61.8	1.2	167,620	266
Davidson	47.3	0.7	493,272	113
Wilson	37.2	1.2	105,888	312
Rutherford	34.4	1.1	223,304	85
Sumner	32.6	1.2	138,609	121
Maury	26.5	1.5	73,091	272
Cheatham	25.1	2.3	29,614	147
Dickson	21.4	1.9	38,559	203
Robertson	21.1	1.6	50,891	165
Smith	17.9	2.3	14,180	77
Cannon	16.9	3.2	10,397	97
Hickman	13.6	2.2	18,044	107
Macon	11.8	2.3	17,135	251
Trousdale	10.1	2.7	8,490	285

# ============================================================
# 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("e8b8b7b0e43cfa561bd5bcc3c19f634a7406836f")

# ============================================================
# 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: Sections 6 through 8 of this script have been omitted.
# See "Practical R (Part 1)" for those sections.
# ============================================================

# ============================================================
# 9. FILTER DATA FOR MULTI-COUNTY AREA
# ============================================================

CountyNameList <- c(
  "Cannon","Cheatham","Davidson","Dickson","Hickman",
  "Macon","Maury","Robertson","Rutherford","Smith",
  "Sumner","Trousdale","Williamson","Wilson"
)

filtereddata2 <- mydata %>%
  filter(County %in% CountyNameList)

# ============================================================
# 10. MAP MULTI-COUNTY AREA (WITH COUNTY BORDERS)
# ============================================================

mapdata2 <- filtereddata2 %>%
  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) %>%
  mutate(County = str_remove(NAME, " County, Tennessee")) %>%
  filter(County %in% CountyNameList)

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

# Build leaflet map with county borders
DivisionMap2 <- leaflet(mapdata2) %>%
  # 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 = mapdata2$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"   # 👈 Control position
    )
  )

DivisionMap2

# ============================================================
# 11. COUNTY-LEVEL TABLE FOR SELECTED COUNTIES
# ============================================================

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

SelectedCounties <- CountyLevelData %>% 
  mutate(NAME = str_remove(NAME, " County, Tennessee")) %>%
  filter(NAME %in% CountyNameList) %>%
  select(NAME, Percent_E, Percent_M, Count_E, Count_M) %>%
  rename(
    County  = NAME,
    Percent = Percent_E,
    PctMOE  = Percent_M,
    Count   = Count_E,
    CountMOE = Count_M
  ) %>%
  arrange(desc(Percent))

CountyTable <- SelectedCounties %>%
  gt() %>%
  gt_theme_espn() %>%
  cols_label(
    Percent = "Percent (%)",
    PctMOE = "Percent MOE",
    Count = "Count",
    CountMOE = "Count MOE"
  ) %>%
  fmt_number(columns = c(Percent, PctMOE), decimals = 1) %>%
  fmt_number(columns = c(Count), decimals = 0, use_seps = TRUE) %>%
  tab_header(
    title = "Estimates by County",
    subtitle = "(2019–2023 ACS 5-Year ACS)"
  )

CountyTable