Intro

For this R script, I chose Knox County, which is in the area where I grew up. For this analysis, I decided to look into the Poverty Rate in 2021, specifically around the time COVID-19 had hit all areas, which caused alot of poverty to surrounding counties.

County-Level Analysis

Aiming near the bottom, Knox County’s estimate rate was at a 12.7 Poverty rate percentage and is at the same poverty rate percentage as Hamilton County in Chattanooga. Surrounding areas such as Anderson County (15.4%) and Roane County (13.7%) hit higher in the poverty rate as well.

Within-County Analysis

The county-level map revealed noticeable variation in poverty rates across Tennessee. When I examined Knox County, one result that stood out was that the county’s poverty rate was lower than many rural counties across the state. Given that Knox County contains the city of Knoxville and several suburban communities, this result was not surprising. However, it highlighted how poverty can vary significantly even within a county that is generally considered economically stable compared to many other parts of Tennessee.

Many of the highest poverty rates were concentrated in rural counties across West Tennessee and parts of East Tennessee, while more urban and suburban counties tended to have lower poverty rates. Knox County fell somewhere in the middle-to-lower range compared to the counties with the highest levels of poverty. This suggests that larger metropolitan areas may provide greater access to employment, education, healthcare, and community resources that can help reduce poverty rates.

To better understand Knox County’s results, I compared the poverty rate to other socioeconomic indicators such as median household income and educational attainment. Counties with higher incomes and larger percentages of residents holding bachelor’s degrees generally tended to have lower poverty rates. Knox County follows this pattern because its relatively strong economy, major employers, and educational institutions, including the University of Tennessee, help support higher incomes and create employment opportunities for residents.

Code

############################################################
# STEP 1: Load Packages
############################################################

if (!require("tidycensus")) install.packages("tidycensus")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("sf")) install.packages("sf")
if (!require("leaflet")) install.packages("leaflet")
if (!require("htmlwidgets")) install.packages("htmlwidgets")
if (!require("plotly")) install.packages("plotly")

library(tidycensus)
library(tidyverse)
library(kableExtra)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

############################################################
# STEP 2: Verify Census API Key
############################################################

if (!nzchar(Sys.getenv("CENSUS_API_KEY"))) {
  
  stop(
    paste(
      "\nNo Census API key was found.",
      "\n\nTo obtain a free Census API key, visit:",
      "\nhttps://api.census.gov/data/key_signup.html",
      "\n\nThen run the one-time setup code at the",
      "\ntop of this script to store your key.",
      "\n"
    )
  )
  
}

############################################################
# STEP 3: Explore 2024 Codebooks and Choose an ACS Variable
############################################################

ProfileTables <- load_variables(2024, "acs5/profile")

# DetailedTables <- load_variables(2024, "acs5")
# SubjectTables  <- load_variables(2024, "acs5/subject")

# Replace the variable code below with the ACS variable
# you want to download.

my_variable <- "DP03_0128P"

# Example variables:
#
# DP05_0001  = Total population
# DP03_0062  = Median Household Income
# DP03_0128P = Poverty Rate
# DP02_0068P = Bachelor's Degree or Higher Rate
# DP05_0018  = Median Age
# DP04_0046P = Homeownership Rate
# DP04_0134  = Median Gross Rent
# DP04_0142P = Pct. renters paying 35%+ HH income in rent (GRAPI)
# DP05_0090P = Hispanic or Latino Population Percentage
# DP02_0094P = Foreign-Born Percentage
# DP03_0009P = Unemployment Rate
# DP03_0025  = Average Commute Time

############################################################
# STEP 4: Download ACS Data
############################################################

# Common ACS geography options:
#
# "state"
# "county"
# "county subdivision"
# "tract"
# "block group"
# "zcta"
#
# Most of this script works for any ACS geography.
# Some geographies may require additional filtering
# or display adjustments.

mydata <- get_acs(
  geography = "county",
  state = "TN",
  year = 2021,
  survey = "acs5",
  variables = my_variable,
  geometry = TRUE
)

############################################################
# STEP 5: Keep Needed Variables and Geometry
############################################################

mydata <- mydata |>
  select(
    Area = NAME,
    Variable = variable,
    Estimate = estimate,
    Margin_of_Error = moe,
    geometry
  )

############################################################
# STEP 5B: Transform Coordinates for Leaflet
############################################################

mydata <- st_transform(
  mydata,
  crs = 4326
)

############################################################
# STEP 6: Optional Filtering
############################################################

# Uncomment and edit the code below if you want to keep only
# selected areas.
#
# To uncomment:
# 1. Select the lines with your mouse.
# 2. Press Ctrl+Shift+C (Windows) or Cmd+Shift+C (Mac).  

# mydata <- mydata |>
#   filter(
#     Area %in% c(
#       "Knox County, Tennessee"
#       "Anderson County, Tennessee"
#       
#     )
#   )

############################################################
# STEP 7: Sort Results from Highest to Lowest
############################################################

mydata <- mydata |>
  arrange(desc(Estimate))

############################################################
# STEP 8: Display Results as a Table
############################################################

ResultsTable <- mydata |>
  st_drop_geometry() |>
  kbl(
    col.names = c(
      "Area",
      "Variable",
      "Estimate",
      "Margin of Error"
    ),
    format.args = list(big.mark = ",")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

ResultsTable

############################################################
# STEP 9: Create Pop-Up Content
############################################################

mydata$popup <- paste0(
  "<strong>", mydata$Area, "</strong><br/>",
  "Estimate: ",
  format(mydata$Estimate, big.mark = ","),
  "<br/>",
  "Plus/Minus: ",
  format(mydata$Margin_of_Error, big.mark = ",")
)

############################################################
# STEP 10: Build Color Palette
############################################################

qs <- unique(
  quantile(
    mydata$Estimate,
    probs = seq(0, 1, length.out = 6),
    na.rm = TRUE
  )
)

pal <- colorBin(
  palette = "Blues",
  domain = mydata$Estimate,
  bins = qs,
  pretty = FALSE
)

############################################################
# STEP 11: Create Choropleth Map
############################################################

AreaMap1 <- leaflet(mydata) |>
  addProviderTiles(
    providers$Esri.WorldTopoMap
  ) |>
  addPolygons(
    fillColor = ~pal(Estimate),
    fillOpacity = 0.5,
    color = "black",
    weight = 1,
    popup = ~popup
  ) |>
  addLegend(
    pal = pal,
    values = ~Estimate,
    title = "Estimate",
    labFormat = labelFormat(
      big.mark = ","
    )
  )

AreaMap1

############################################################
# STEP 12: Create Plotly Dot Plot
############################################################

plotdf <- mydata |>
  st_drop_geometry() |>
  mutate(
    point_color = pal(Estimate),
    y_ordered = reorder(
      Area,
      Estimate
    ),
    hover_text = paste0(
      Area,
      "
Estimate: ",
      format(Estimate, big.mark = ","),
      "
Plus/Minus: ",
      format(Margin_of_Error, big.mark = ",")
    )
  )

DotPlot <- plot_ly(
  data = plotdf,
  x = ~Estimate,
  y = ~as.character(y_ordered),
  type = "scatter",
  mode = "markers",
  showlegend = FALSE,
  marker = list(
    color = ~point_color,
    size = 8,
    line = list(
      color = "rgba(120,120,120,0.9)",
      width = 0.5
    )
  ),
  error_x = list(
    type = "data",
    array = ~Margin_of_Error,
    arrayminus = ~Margin_of_Error,
    color = "rgba(0,0,0,0.65)",
    thickness = 1
  ),
  text = ~hover_text,
  hovertemplate = "%{text}"
) |>
  layout(
    title = list(
      text = paste0(
        "ACS Estimates by Area",
        "
Error bars show ACS margins of error."
      )
    ),
    xaxis = list(
      title = "Estimate",
      tickformat = ",.0f",
      automargin = TRUE
    ),
    yaxis = list(
      title = "",
      automargin = TRUE,
      categoryorder = "array",
      categoryarray = levels(plotdf$y_ordered)
    ),
    margin = list(
      l = 200,
      r = 20,
      b = 60,
      t = 60,
      pad = 2
    )
  )

DotPlot

############################################################
# STEP 13: Export Dot Plot (Optional)
############################################################

saveWidget(
  widget = as_widget(DotPlot),
  file = "ACSGraph.html",
  selfcontained = TRUE
)

############################################################
# STEP 14: Export Map (Optional)
############################################################

saveWidget(
  widget = AreaMap,
  file = "ACSMap.html",
  selfcontained = TRUE
)
############################################################
# STEP 1: Load Packages
############################################################

if (!require("tidycensus")) install.packages("tidycensus")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("kableExtra")) install.packages("kableExtra")
if (!require("sf")) install.packages("sf")
if (!require("leaflet")) install.packages("leaflet")
if (!require("htmlwidgets")) install.packages("htmlwidgets")
if (!require("plotly")) install.packages("plotly")

library(tidycensus)
library(tidyverse)
library(kableExtra)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

############################################################
# STEP 2: Verify Census API Key
############################################################

if (!nzchar(Sys.getenv("CENSUS_API_KEY"))) {
  
  stop(
    paste(
      "\nNo Census API key was found.",
      "\n\nTo obtain a free Census API key, visit:",
      "\nhttps://api.census.gov/data/key_signup.html",
      "\n\nThen run the one-time setup code at the",
      "\ntop of this script to store your key.",
      "\n"
    )
  )
  
}

############################################################
# STEP 3: Explore 2021 Codebooks and Choose an ACS Variable
############################################################

ProfileTables <- load_variables(2021, "acs5/profile")

# DetailedTables <- load_variables(2024, "acs5")
# SubjectTables  <- load_variables(2024, "acs5/subject")

# Replace the variable code below with the ACS variable
# you want to download.

my_variable <- "DP02_0128P"

# Example variables:
#
# DP05_0001  = Total population
# DP03_0062  = Median Household Income
# DP03_0128P = Poverty Rate
# DP02_0068P = Bachelor's Degree or Higher Rate
# DP05_0018  = Median Age
# DP04_0046P = Homeownership Rate
# DP04_0134  = Median Gross Rent
# DP04_0142P = Pct. renters paying 35%+ HH income in rent (GRAPI)
# DP05_0090P = Hispanic or Latino Population Percentage
# DP02_0094P = Foreign-Born Percentage
# DP03_0009P = Unemployment Rate
# DP03_0025  = Average Commute Time

############################################################
# STEP 4: Download ACS Data
############################################################

# Common ACS geography options:
#
# "state"
# "county"
# "county subdivision"
# "tract"
# "block group"
# "zcta"
#
# Most of this script works for any ACS geography.
# Some geographies may require additional filtering
# or display adjustments.

mydata <- get_acs(
  geography = "county subdivision",
  state = "TN",
  year = 2021,
  survey = "acs5",
  variables = my_variable,
  geometry = TRUE
)

############################################################
# STEP 5: Keep Needed Variables and Geometry
############################################################

mydata <- mydata |>
  select(
    Area = NAME,
    Variable = variable,
    Estimate = estimate,
    Margin_of_Error = moe,
    geometry
  )

############################################################
# STEP 5B: Transform Coordinates for Leaflet
############################################################

mydata <- st_transform(
  mydata,
  crs = 4326
)

############################################################
# STEP 6: Optional Filtering
############################################################

# OPTION A: Filter for specific counties.
#
# Uncomment to keep only county subdivisions whose
# names contain the specified county names. To search
# for just a single county, use just the county name
# without the | character.

mydata <- mydata |>
  filter(
    str_detect(
      Area,
      "Knox"
    )
  )

############################################################
# STEP 7: Sort Results from Highest to Lowest
############################################################

mydata <- mydata |>
  arrange(desc(Estimate))

############################################################
# STEP 8: Display Results as a Table
############################################################

ResultsTable <- mydata |>
  st_drop_geometry() |>
  kbl(
    col.names = c(
      "Area",
      "Variable",
      "Estimate",
      "Margin of Error"
    ),
    format.args = list(big.mark = ",")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

ResultsTable

############################################################
# STEP 9: Create Pop-Up Content
############################################################

mydata$popup <- paste0(
  "", mydata$Area, "
",
"Estimate: ",
format(mydata$Estimate, big.mark = ","),
"
",
"Plus/Minus: ",
format(mydata$Margin_of_Error, big.mark = ",")
)

############################################################
# STEP 10: Build Color Palette
############################################################

qs <- unique(
  quantile(
    mydata$Estimate,
    probs = seq(0, 1, length.out = 6),
    na.rm = TRUE
  )
)

pal <- colorBin(
  palette = "Blues",
  domain = mydata$Estimate,
  bins = qs,
  pretty = FALSE
)

############################################################
# STEP 11: Create Choropleth Map
############################################################

AreaMap <- leaflet(mydata) |>
  addProviderTiles(
    providers$Esri.WorldTopoMap
  ) |>
  addPolygons(
    fillColor = ~pal(Estimate),
    fillOpacity = 0.5,
    color = "black",
    weight = 1,
    popup = ~popup
  ) |>
  addLegend(
    pal = pal,
    values = ~Estimate,
    title = "Estimate",
    labFormat = labelFormat(
      big.mark = ","
    )
  )

AreaMap

############################################################
# STEP 12: Create Plotly Dot Plot
############################################################

plotdf <- mydata |>
  st_drop_geometry() |>
  mutate(
    point_color = pal(Estimate),
    
    y_ordered = reorder(
      Area,
      Estimate
    ),
    
    hover_text = paste0(
      Area,
      "
Estimate: ",
      format(Estimate, big.mark = ","),
      "
Plus/Minus: ",
      format(Margin_of_Error, big.mark = ",")
    )
  )

DotPlot <- plot_ly(
  data = plotdf,
  x = ~Estimate,
  y = ~as.character(y_ordered),
  type = "scatter",
  mode = "markers",
  showlegend = FALSE,
  marker = list(
    color = ~point_color,
    size = 8,
    line = list(
      color = "rgba(120,120,120,0.9)",
      width = 0.5
    )
  ),
  error_x = list(
    type = "data",
    array = ~Margin_of_Error,
    arrayminus = ~Margin_of_Error,
    color = "rgba(0,0,0,0.65)",
    thickness = 1
  ),
  text = ~hover_text,
  hovertemplate = "%{text}"
) |>
  layout(
    title = list(
      text = paste0(
        "ACS Estimates by State",
        "
Error bars show ACS margins of error."
      )
    ),
    xaxis = list(
      title = "Estimate",
      tickformat = ",.0f",
      automargin = TRUE
    ),
    yaxis = list(
      title = "",
      automargin = TRUE,
      categoryorder = "array",
      categoryarray = levels(plotdf$y_ordered)
    ),
    margin = list(
      l = 200,
      r = 20,
      b = 60,
      t = 60,
      pad = 2
    )
  )

DotPlot

############################################################
# STEP 13: Export Dot Plot (Optional)
############################################################

saveWidget(
  widget = as_widget(DotPlot),
  file = "ACSGraph.html",
  selfcontained = TRUE
)

############################################################
# STEP 14: Export Map (Optional)
############################################################

saveWidget(
  widget = AreaMap,
  file = "ACSMap.html",
  selfcontained = TRUE
)