Who’s “rich” in Rutherford County? With data from the U.S. Census Bureau’s latest American Community Survey and some code from the open-source R Project for Statistical Computing, you can compare median household income across Rutherford County, the greater Nashville area, and the whole state. MTSU’s School of Journalism and Strategic Media offers instruction in data journalism using R and other data tools.

Median household income estimates

Rutherford County, Tennessee (It’s a touch screen; try it out!)

Median household income estimates

Rutherford, Davidson and Williamson counties, Tennessee

Median household income estimates

Tennessee

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

# census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")

# Fetching ACS codebooks

DetailedTables <- load_variables(2022, "acs5", cache = TRUE)
SubjectTables <- load_variables(2022, "acs5/subject", cache = TRUE)
ProfileTables <- load_variables(2022, "acs5/profile", cache = TRUE)

# Specifying target variables

VariableList = 
  c(Estimate_ = "DP03_0062")

# Some potentially interesting ACS variables

# DP05_0037P - Pct. white alone
# DP03_0074P - Pct. with SNAP benefits in past 12 months
# 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
# DP02PR_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("Rutherford 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 = "darkblue") + 
  theme_minimal(base_size = 12.5) + 
  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)

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

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)

mapdata2 <- st_as_sf(mapdata2)

mapviewOptions(basemaps.color.shuffle = FALSE)

DivisionMap2 <- mapview(mapdata2,
                       zcol = "Estimate",
                       layer.name = "Estimate",
                       popup = popupTable(
                         mapdata2,
                         feature.id = FALSE,
                         row.numbers = FALSE,
                         zcol = c(
                           "State",
                           "County",
                           "Division",
                           "Estimate",
                           "Range")))

DivisionMap2

# Mapping the Whole state

filtereddata3 <- mydata

# Mapping multi-county data

mapdata3 <- filtereddata3 %>% 
  rename(Estimate = Estimate_E,
         Range = Estimate_M)

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

DivisionMap3