County-Level Analysis
For my county level analysis, I chose Hughes County, South Dakota which is home to the State Capital Pierre. I actually used to live in Pierre, SD back when I worked at KCCR Radio in 2025 which is why I was curious to learn more information about it. I chose to look at the median age for the county as that was always something that fascinated me when I tried to make friends and find people my age. Compared to the neighboring counties, Hughes County at 39.5 was middle of the road with Buffalo County, Dewey County, and Lyman County all having younger ages. Meanwhile, Hyde County, Stanley County, and Sully County all had larger median ages. It did surprise me to see that Hughes County was not the youngest county in the area overall considering it is the biggest city in the area by a wide margin, but with the surrounding counties being ranch areas, perhaps that should have been something I expected going in. The map for this analysis is displayed below.
Within-County Analysis
For the within-county analysis, therer was significant variation within the county itself. The City of Pierre had the youngest age at 38.5 while North Hughes and Crow Creek each at ages at 60+ years old. It seemed as though the further east you go, the older people got in the county as a whole. The Pierre City and West Hughes areas each had the youngest populations at 39.4 and 38.5 for median age respectively. In addition to my previous points, I also noticed that the further people got away from the City of Pierre proper, the older their average age usually ended up being. This reflects the expectation that as you get further from the city, the age usually goes up for an average citizen. I think what largely explains the age discrepancies is that Pierre is home to young state government employees looking to get their first big job in life while going out further is where the rich ranchers live in solitude away from the hustle and bustle of the city. It’s also overall a better place and quieter place to raise a family and set down roots. This is demonstrated in the map below.
Reflection
Before starting my analysis, I expected to see mostly what was shown with the within-county analysis when it came to median age, since the bigger cities are where the majority of a younger population resides. I was somewhat shocked by the age discrepancy in the county analysis with the fact that some neighboring counties had younger populations than Hughes County, which I thought would for sure have the youngest. In addition, I thought that areas with more ranches and open terrain would equal older people, but that was not always the case as previously shown. One thing that was hidden in a big way with the subcounty map was that the City of Pierre does not take up that much of Hughes County, and it is obvious to see that it is a very densely populated area. There is a lot more open space than what might previously think in a county of that size. If I was a journalist, I would likely want to write a story about how older folks who live in Crow Creek or North Hughes deal with the commutes to Pierre to get groceries/supplies and why they chose to settle down in such an isolated area? I would also want to compare the average day in the life of someone who is in the city compared to someone who is in those other areas to demonstrate the wildly different lifestyles and people within the same county.
Code
Below is the first code for the County Map that gave me Hughes County and its neighbors. The second code after that is what allowed me to access the subcounty map for the within-county analysis.
############################################################
# RUN ONCE TO SAVE YOUR CENSUS API KEY
#
# Run this section once, then delete it, so that the
# script begins with STEP 1, below
############################################################
# if (!require("tidycensus")) install.packages("tidycensus")
# library(tidycensus)
#
# census_api_key(
# "PASTE_YOUR_CENSUS_API_KEY_BETWEEN_THESE_QUOTES",
# install = TRUE,
# overwrite = TRUE
# )
#
# readRenviron("~/.Renviron")
############################################################
# 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 <- "DP05_0018"
# 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 = "SD",
year = 2024,
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(
"Hughes County, South Dakota",
"Hyde County, South Dakota",
"Sully County, South Dakota",
"Buffalo County, South Dakota",
"Lyman County, South Dakota",
"Dewey County, South Dakota",
"Stanley County, South Dakota"
)
)
############################################################
# 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
############################################################
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 = ","
)
)
CountyMap <- AreaMap
CountyMap
############################################################
# 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
)############################################################
# RUN ONCE TO SAVE YOUR CENSUS API KEY
#
# Run this section once, then delete it, so that the
# script begins with STEP 1, below
############################################################
# if (!require("tidycensus")) install.packages("tidycensus")
# library(tidycensus)
#
# census_api_key(
# "PASTE_YOUR_CENSUS_API_KEY_BETWEEN_THESE_QUOTES",
# install = TRUE,
# overwrite = TRUE
# )
#
# readRenviron("~/.Renviron")
############################################################
# 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 <- "DP05_0018"
# 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 = "SD",
year = 2024,
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(
str_detect(
Area,
"Hughes"
)
)
############################################################
# 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
############################################################
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 = ","
)
)
SubCountyMap <- AreaMap
SubCountyMap
############################################################
# 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
)