Data journalism with R
 |
 |
Who goes hungry if SNAP gets cut?
With R code, journalists can get neighborhood-level American Community Survey data for Davidson County, the Nashville area, and the whole state. For example, these maps show percentages of households receiving Supplemental Nutrition Assistance Program benefits.
MTSU’s School of Journalism and Strategic Media and its new, fully online, one-year Digital Media Master’s degree both offer instruction in R-based data journalism and data analysis.
Percent of homes with SNAP assistance
Davidson County, 2019-2023
Percent of homes with SNAP assistance
Davidson, Rutherford and Williamson counties, 2019-2023
Percent of homes with SNAP assistance
Tennessee, 2019-2023
Code
The R code shown below will generate each of the maps shown above. For a tutorial on how to analyze American Community Survey data using R, See: Data Wizardry with R.
# Installing and loading 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")
library(tidyverse)
library(tidycensus)
library(sf)
library(mapview)
library(leaflet)
library(leafpop)
# Transmitting API key obtained from:
# https://api.census.gov/data/key_signup.html
census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")
# Fetching 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)
# Specifying target variables
VariableList =
c(Estimate_ = "DP03_0074P",
Households_ = "DP02_0001",
Population_ = "DP02_0124")
# Some potentially interesting ACS variables
# DP02_0124 - Total population
# DP05_0037P - Pct. white alone
# DP03_0074P - Pct. HH with SNAP benefits in past 12 months
# DP02_0001 - Total households
# DP03_0072P - Pct. with cash public assistance
# DP03_0128P - Pct. all people below poverty level
# DP03_0099P - Pct. no health insurance coverage
# DP04_0142P - GRAPI 35% or more
# DP04_0115P - Sel. Mo. Owner Costs 35% or more
# DP03_0092 - Median earnings for workers
# DP03_0093 - Median earnings, full-time year-round male workers
# DP03_0094 - Median earnings, full-time year-round female workers
# DP03_0062 - Median household income
# DP02_0068P - Pct. 25 or older with a bachelor's degree or higher
# ---------- County subdivision analysis ----------
# Fetching data
mydata <- get_acs(
geography = "county subdivision",
state = "TN",
variables = VariableList,
year = 2023,
survey = "acs5",
output = "wide",
geometry = TRUE)
# Reformatting data
mydata <-
separate_wider_delim(mydata,
NAME,
delim = ", ",
names = c("Division", "County", "State"))
# Filtering data
filtereddata <- mydata %>%
filter(County %in% c("Davidson County"))
# Plotting data
mygraph <- ggplot(filtereddata, aes(x = Estimate_E, y = reorder(Division, Estimate_E))) +
geom_errorbarh(aes(xmin = Estimate_E - Estimate_M, xmax = Estimate_E + Estimate_M)) +
geom_point(size = 3, color = "#099d91") +
theme_minimal(base_size = 12.5) +
theme(axis.text.y = element_text(size = 8)) +
labs(title = "Estimates by area",
subtitle = "County subdivisions. Brackets show error margins.",
x = "2019-2023 ACS estimate",
y = "")
mygraph
# Mapping data
mapdata <- filtereddata %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata <- st_as_sf(mapdata)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap <- mapview(mapdata,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households")))
DivisionMap
# Filtering for Nashville-area counties
filtereddata2 <- mydata %>%
filter(County %in% c("Rutherford County",
"Davidson County",
"Williamson County"))
# Mapping multi-county data
mapdata2 <- filtereddata2 %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata2 <- st_as_sf(mapdata2)
# Fetching county borders
Counties <- get_acs(
geography = "county",
state = c("TN"),
county = c("Rutherford", "Davidson", "Williamson"),
variables = VariableList,
year = 2023,
survey = "acs5",
output = "wide",
geometry = TRUE)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap2 <- mapview(
Counties,
zcol = NULL,
legend = FALSE,
col.regions = "transparent",
alpha.regions = 0,
lwd = 1.5,
color = "black"
) + mapview(
mapdata2,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata2,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households"
)
)
)
DivisionMap2
# Mapping the Whole state
filtereddata3 <- mydata
# Mapping multi-county data
mapdata3 <- filtereddata3 %>%
rename(Estimate = Estimate_E,
Range = Estimate_M,
Population = Population_E,
Households = Households_E)
mapdata3 <- st_as_sf(mapdata3)
mapviewOptions(basemaps.color.shuffle = FALSE)
DivisionMap3 <- mapview(mapdata3,
zcol = "Estimate",
layer.name = "Estimate",
popup = popupTable(
mapdata3,
feature.id = FALSE,
row.numbers = FALSE,
zcol = c(
"State",
"County",
"Division",
"Estimate",
"Range",
"Population",
"Households")))
DivisionMap3