Practical R 3

Percent with a Bachelors degree or higher
Rutherford County, 2019-2023. Source: American Community Survey
Percent with a Bachelor’s degree or higher
Nashville MSA, 2019-2023. Source: American Community Survey
County	Percent (%)	Percent MOE	Count	Count MOE
Estimates by County
(2019–2023 ACS 5-Year ACS)
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
Region Table and Statewide Map
Region	Percent (%)	Margin of Error	Count
Estimates by Region
(2019–2023 ACS 5-Year ACS)
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
# ============================================================
# 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("mapview"))
  install.packages("mapview")
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(mapview)
library(leaflet)
library(leafpop)
library(gt)
library(gtExtras)
library(plotly)

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

# NOTE: Replace PasteYourAPIKeyBetweenTheseQuoteMarks with your
# actual API key. If you don't the script won't work.

# ============================================================

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"))

# ============================================================
# 6. FILTER FOR A SINGLE COUNTY
# ============================================================

filtereddata <- mydata %>%
  filter(County %in% c("Rutherford"))

# ============================================================
# 7. MAP DATA FOR SINGLE COUNTY SUBDIVISIONS
# ============================================================

mapdata <- filtereddata %>%
  rename(
    Percent = Percent_E,
    PercentEM = Percent_M,
    Count = Count_E,
    CountEM = Count_M
  ) %>%
  st_as_sf()

mapviewOptions(basemaps.color.shuffle = FALSE)

DivisionMap <- mapview(
  mapdata,
  zcol = "Percent",
  layer.name = "Percent",
  popup = popupTable(
    mapdata,
    feature.id  = FALSE,
    row.numbers = FALSE,
    zcol = c("State", "County", "Division", "Percent",
             "PercentEM", "Count", "CountEM")
  )
)

DivisionMap

# ============================================================
# 8. INTERACTIVE PLOTLY GRAPH OF ESTIMATES WITH ERROR BARS
# ============================================================

mygraph <- plot_ly(
  data = filtereddata,
  x = ~Percent_E,
  y = ~reorder(Division, Percent_E),
  type = 'scatter',
  mode = 'markers',
  error_x = list(
    type = "data",
    array = ~Percent_M,
    visible = TRUE
  ),
  marker = list(color = "#099d91", size = 10)
) %>%
  layout(
    title = "Estimates by Area (Interactive)",
    xaxis = list(title = "2019–2023 ACS Estimate"),
    yaxis = list(title = "", automargin = TRUE)
  )

mygraph

# ============================================================
# 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

# ============================================================
# 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: 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
Practical R 3

Gracie Stanley

2025-12-10
