# ---------------------------------------------------------------------------
# double_jeopardy_map.R (tidycensus version)
#
# Recreates "Double jeopardy: shrinking and aging" -- a 2x2 bivariate county
# map of the 50 states + DC:
#
# Population Growing / Declining population 2025 vs. 2000
# Age Older / Younger 2025 residents 65+ vs. under 15
#
# "Double jeopardy" = Declining AND Older (dark red).
#
# Data, all through tidycensus:
# 2000 get_decennial(year = 2000, sumfile = "sf1") -- Census 2000 counts
# (counties, plus CT towns and CO places for the two fixes below)
# 2025 get_estimates(vintage = 2025, product = "characteristics",
# breakdown = "AGEGROUP") -- July 1, 2025 PEP by age
# Boundaries: tigris cartographic counties (20m), Alaska/Hawaii shifted.
# Metro status: OMB July 2023 delineation file (not in tidycensus).
#
# Geography harmonization (so 2000 and 2025 compare like-for-like):
# * Virginia independent cities merged into their surrounding county
# (BEA-style combinations); Bedford city and Clifton Forge folded in too.
# * Alaska boroughs split/redrawn since 2000 recombined to stable units.
# * Shannon SD (46113) -> Oglala Lakota (46102).
# * Connecticut: 2025 data exist only by planning region, so 2000 population
# is rebuilt for planning regions by summing 2000 towns.
# * Broomfield CO didn't exist in 2000; its 2000 population is taken from
# Broomfield city (2000 census place). Parent counties are left as-is;
# all four grew strongly, so their classification doesn't change.
#
# Requires: tidycensus, tigris, sf, dplyr, tidyr, ggplot2, patchwork, readxl
# (optional: ragg, systemfonts for the Aptos font)
# Census API key (free: https://api.census.gov/data/key_signup.html):
# tidycensus::census_api_key("YOUR_KEY", install = TRUE); readRenviron("~/.Renviron")
# ---------------------------------------------------------------------------
suppressPackageStartupMessages({
library(tidycensus); library(tigris); library(sf)
library(dplyr); library(tidyr); library(ggplot2); library(patchwork)
})
options(tigris_use_cache = TRUE)
# --- settings --------------------------------------------------------------
out_png <- "double_jeopardy_map.png"
out_csv <- "double_jeopardy_counties.csv"
vintage <- 2025 # PEP vintage; set 2024 if 2025 isn't in your tidycensus
cache_dir <- "data_census" # every download is saved here as .rds
refresh <- FALSE # TRUE to ignore the cache and pull everything fresh
url_metro <- "https://www2.census.gov/programs-surveys/metro-micro/geographies/reference-files/2023/delineation-files/list1_2023.xlsx"
pal <- c("Younger Growing" = "#ECE6DA",
"Younger Declining" = "#9FB3C4",
"Older Growing" = "#E1A563",
"Older Declining" = "#6B1E2C")
font <- if (requireNamespace("systemfonts", quietly = TRUE) &&
"Aptos" %in% systemfonts::system_fonts()$family) "Aptos" else "sans"
# --- cache -------------------------------------------------------------------
# cached("name", expr): return data_census/name.rds if it exists, otherwise
# evaluate expr, save the result, and return it. Failed pulls aren't saved,
# so a vintage tidycensus doesn't have yet is retried on the next run.
dir.create(cache_dir, showWarnings = FALSE, recursive = TRUE)
cached <- function(name, expr) {
f <- file.path(cache_dir, paste0(name, ".rds"))
if (!refresh && file.exists(f)) {
message(" cache: ", basename(f))
return(readRDS(f))
}
x <- expr # evaluated only on a cache miss
saveRDS(x, f)
x
}
if (!nzchar(Sys.getenv("CENSUS_API_KEY")))
stop("No Census API key found. Run once:\n",
' tidycensus::census_api_key("YOUR_KEY", install = TRUE)\n',
' readRenviron("~/.Renviron")', call. = FALSE)
# 50 states + DC (drops PR and territories)
keep_states <- unique(fips_codes$state_code[fips_codes$state_code <= "56"])
# --- harmonization ---------------------------------------------------------
# County FIPS -> harmonized unit id. Anything not listed maps to itself.
unit_map <- c(
# Virginia independent cities -> surrounding county
"51540" = "51003", # Charlottesville -> Albemarle
"51580" = "51005", "51560" = "51005", # Covington, Clifton Forge -> Alleghany
"51790" = "51015", "51820" = "51015", # Staunton, Waynesboro -> Augusta
"51515" = "51019", # Bedford city -> Bedford
"51680" = "51031", # Lynchburg -> Campbell
"51640" = "51035", # Galax -> Carroll
"51570" = "51053", "51730" = "51053", # Colonial Hts, Petersburg -> Dinwiddie
"51600" = "51059", "51610" = "51059", # Fairfax city, Falls Church -> Fairfax
"51840" = "51069", # Winchester -> Frederick
"51595" = "51081", # Emporia -> Greensville
"51690" = "51089", # Martinsville -> Henry
"51830" = "51095", # Williamsburg -> James City
"51750" = "51121", # Radford -> Montgomery
"51590" = "51143", # Danville -> Pittsylvania
"51670" = "51149", # Hopewell -> Prince George
"51683" = "51153", "51685" = "51153", # Manassas, Manassas Park -> Prince William
"51775" = "51161", # Salem -> Roanoke County
"51530" = "51163", "51678" = "51163", # Buena Vista, Lexington -> Rockbridge
"51660" = "51165", # Harrisonburg -> Rockingham
"51620" = "51175", # Franklin city -> Southampton
"51630" = "51177", # Fredericksburg -> Spotsylvania
"51520" = "51191", # Bristol -> Washington
"51720" = "51195", # Norton -> Wise
"51735" = "51199", # Poquoson -> York
# Alaska: recombine areas redrawn since 2000
"02230" = "02232", "02105" = "02232", # Skagway + Hoonah-Angoon
"02063" = "02261", "02066" = "02261", # Chugach + Copper River (Valdez-Cordova)
"02130" = "02SE", "02198" = "02SE", "02195" = "02SE", # Ketchikan / PoW /
"02275" = "02SE", "02201" = "02SE", "02280" = "02SE", # Petersburg / Wrangell
"02270" = "02158", # Wade Hampton -> Kusilvak (rename)
"46113" = "46102" # Shannon -> Oglala Lakota (rename)
)
to_unit <- function(geoid) ifelse(geoid %in% names(unit_map), unit_map[geoid], geoid)
# --- 2000 population (Census 2000, SF1) ------------------------------------
p00 <- cached("dec2000_county",
get_decennial(geography = "county", variables = "P001001",
year = 2000, sumfile = "sf1")) |>
filter(substr(GEOID, 1, 2) %in% keep_states, substr(GEOID, 1, 2) != "09") |>
select(GEOID, pop2000 = value)
# Connecticut: 2000 towns -> planning regions (towns nest in both systems)
ct_towns <- cached("dec2000_ct_towns",
get_decennial(geography = "county subdivision", state = "CT",
variables = "P001001", year = 2000, sumfile = "sf1")) |>
mutate(COUSUBFP = substr(GEOID, 6, 10))
ct_xw <- cached("ct_town_region_xwalk",
county_subdivisions("CT", year = 2024, cb = TRUE) |>
st_drop_geometry() |> distinct(COUSUBFP, COUNTYFP))
ct00 <- inner_join(ct_towns, ct_xw, by = "COUSUBFP") |>
group_by(GEOID = paste0("09", COUNTYFP)) |>
summarise(pop2000 = sum(value), .groups = "drop")
message(sprintf("CT: %d of %d towns matched to planning regions",
sum(ct_towns$COUSUBFP %in% ct_xw$COUSUBFP), nrow(ct_towns)))
# Broomfield CO: county created 2001 from Broomfield city
broomfield <- cached("dec2000_co_places",
get_decennial(geography = "place", state = "CO",
variables = "P001001", year = 2000, sumfile = "sf1")) |>
filter(grepl("^Broomfield", NAME)) |>
transmute(GEOID = "08014", pop2000 = value)
p00 <- bind_rows(p00, ct00, broomfield) |>
group_by(unit = to_unit(GEOID)) |>
summarise(pop2000 = sum(pop2000), .groups = "drop")
# --- 2025 population and age structure (PEP) ------------------------------
# tidycensus adds each PEP product as the Census Bureau posts it, so the age
# file can lag the totals by a release. Try the requested vintage first and
# step back one year if tidycensus doesn't have it yet.
try_vintages <- function(f, vintages) {
for (v in vintages) {
res <- tryCatch(f(v), error = function(e) {
message(sprintf(" vintage %d not available: %s", v, conditionMessage(e))); NULL })
if (!is.null(res) && nrow(res) > 0) return(list(data = res, vintage = v))
}
stop("No vintage available among: ", paste(vintages, collapse = ", "), call. = FALSE)
}
# Total population (drives Growing / Declining)
tot <- try_vintages(function(v)
cached(sprintf("pep_v%d_county_pop", v),
get_estimates(geography = "county", variables = "POPESTIMATE",
vintage = v, year = v)),
c(vintage, vintage - 1))
pop_vintage <- tot$vintage
# Age groups (drives Older / Younger)
age <- try_vintages(function(v)
cached(sprintf("pep_v%d_county_age", v),
get_estimates(geography = "county", product = "characteristics",
breakdown = "AGEGROUP", breakdown_labels = FALSE,
vintage = v, year = v)),
c(vintage, vintage - 1))
age_vintage <- age$vintage
message(sprintf("Using population from Vintage %d, age structure from Vintage %d",
pop_vintage, age_vintage))
e25 <- age$data |>
filter(substr(GEOID, 1, 2) %in% keep_states) |>
mutate(AGEGROUP = as.integer(as.character(AGEGROUP)))
# If tidycensus also returns a SEX column, keep only the both-sexes total
if ("SEX" %in% names(e25))
e25 <- filter(e25, as.character(SEX) %in% c("0", "Both sexes"))
# AGEGROUP codes: 0 = total, 1-3 = 0-14, 14-18 = 65+
p_age <- e25 |>
group_by(unit = to_unit(GEOID)) |>
summarise(u15 = sum(value[AGEGROUP %in% 1:3]),
o65 = sum(value[AGEGROUP %in% 14:18]),
.groups = "drop")
p_tot <- tot$data |>
filter(substr(GEOID, 1, 2) %in% keep_states) |>
group_by(unit = to_unit(GEOID)) |>
summarise(pop2025 = sum(value), .groups = "drop")
p25 <- full_join(p_tot, p_age, by = "unit")
# --- classify --------------------------------------------------------------
dat <- full_join(p25, p00, by = "unit") |>
mutate(growth = if_else(pop2025 < pop2000, "Declining", "Growing"),
age = if_else(o65 > u15, "Older", "Younger"),
cls = if_else(is.na(growth) | is.na(age), NA_character_,
paste(age, growth)),
pct_change = round(100 * (pop2025 / pop2000 - 1), 1))
unmatched <- filter(dat, is.na(cls))
if (nrow(unmatched) > 0)
warning(nrow(unmatched), " unit(s) missing a year: ",
paste(unmatched$unit, collapse = ", "))
# Metro status (a merged unit is metro if any piece is)
metro <- tryCatch({
f <- file.path(cache_dir, basename(url_metro))
if (refresh || !file.exists(f)) download.file(url_metro, f, mode = "wb", quiet = TRUE)
d <- readxl::read_excel(f, skip = 2)
d <- d[!is.na(d$`FIPS County Code`), ]
d |>
filter(`Metropolitan/Micropolitan Statistical Area` == "Metropolitan Statistical Area") |>
transmute(unit = to_unit(paste0(`FIPS State Code`, `FIPS County Code`))) |>
distinct() |> mutate(metro = TRUE)
}, error = function(e) { warning("metro delineation failed: ", conditionMessage(e)); NULL })
dat <- if (is.null(metro)) mutate(dat, metro = NA) else
left_join(dat, metro, by = "unit") |> mutate(metro = coalesce(metro, FALSE))
n_units <- sum(!is.na(dat$cls))
n_dj <- sum(dat$cls == "Older Declining", na.rm = TRUE)
n_dj_nm <- sum(dat$cls == "Older Declining" & !dat$metro, na.rm = TRUE)
message(sprintf("%d of %d units in double jeopardy (%.0f%%), %d nonmetro",
n_dj, n_units, 100 * n_dj / n_units, n_dj_nm))
print(count(dat, cls))
write.csv(arrange(dat, unit), out_csv, row.names = FALSE)
# --- geometry ---------------------------------------------------------------
# get_estimates() returns a plain tibble for the 2020s characteristics files
# even with geometry = TRUE, so boundaries come from tigris (same source
# tidycensus uses). 2024 cartographic counties include CT planning regions.
geo_res <- "5m" # "500k" = most detail, "20m" = coarsest
# The dissolve + shift is the slow step, so the finished shapes are cached too
# (one file per resolution; changing unit_map needs refresh = TRUE).
geo <- cached(paste0("geo_units_", geo_res),
counties(cb = TRUE, resolution = geo_res, year = 2024) |>
filter(STATEFP %in% keep_states) |>
mutate(unit = to_unit(GEOID)) |>
group_by(unit) |>
summarise(.groups = "drop") |> # dissolves merged VA/AK pieces
shift_geometry()) |> # Alaska & Hawaii inset, Albers
left_join(dat, by = "unit")
# State outlines, same resolution and projection as the counties
state_lines <- cached(paste0("geo_states_", geo_res),
states(cb = TRUE, resolution = geo_res, year = 2024) |>
filter(STATEFP %in% keep_states) |>
shift_geometry())
# Move Hawaii to sit below Louisiana. shift_geometry() only offers "below"
# (Alaska & Hawaii in the lower left) or "outside", so Hawaii is nudged
# afterward: centered on Louisiana's x, its top hi_gap below LA's coast.
hi_gap <- 75000 # meters (map units are EPSG:5070-style Albers)
la_bb <- st_bbox(filter(state_lines, STATEFP == "22"))
hi_bb <- st_bbox(filter(state_lines, STATEFP == "15"))
hi_shift <- unname(c((la_bb["xmin"] + la_bb["xmax"]) / 2 - (hi_bb["xmin"] + hi_bb["xmax"]) / 2,
(la_bb["ymin"] - hi_gap) - hi_bb["ymax"]))
move_hawaii <- function(x, is_hi) {
g <- st_geometry(x)
crs <- st_crs(g)
g[is_hi] <- st_set_crs(g[is_hi] + hi_shift, crs)
st_geometry(x) <- g
x
}
geo <- move_hawaii(geo, substr(geo$unit, 1, 2) == "15")
state_lines <- move_hawaii(state_lines, state_lines$STATEFP == "15")
# --- map ---------------------------------------------------------------------
map <- ggplot(geo) +
geom_sf(aes(fill = cls), colour = NA) +
geom_sf(data = state_lines, fill = NA, colour = "grey55", linewidth = 0.15) +
scale_fill_manual(values = pal, na.value = "grey85", guide = "none") +
theme_void(base_family = font)
leg_df <- expand.grid(growth = factor(c("Growing", "Declining"), c("Growing", "Declining")),
age = factor(c("Younger", "Older"), c("Younger", "Older"))) |>
mutate(cls = paste(age, growth))
legend <- ggplot(leg_df, aes(growth, age, fill = cls)) +
geom_tile(colour = "white", linewidth = 1.2) +
annotate("text", x = 2, y = 2, label = "Double\njeopardy",
colour = "white", size = 2.6, family = font, lineheight = 0.9) +
scale_fill_manual(values = pal, guide = "none") +
coord_equal(expand = FALSE) +
labs(x = "Population", y = "Age") +
theme_minimal(base_size = 7.5, base_family = font) +
theme(panel.grid = element_blank(),
axis.title = element_text(face = "bold", size = 8),
axis.text.y = element_text(angle = 90, hjust = 0.5),
plot.background = element_rect(fill = "white", colour = NA))
subtitle <- sprintf(paste0(
"Counties that lost population since 2000 AND now have more residents 65 and older than under 15. ",
"%s counties (%.0f%% of the country, dark red)\nare in double jeopardy on a 21st-century basis",
"%s"),
format(n_dj, big.mark = ","), 100 * n_dj / n_units,
if (is.null(metro)) "." else sprintf(" \u2014 %s of them nonmetro.", format(n_dj_nm, big.mark = ",")))
# Legend sits in its own column right of the map, so it can't cover Florida.
# The spacers push it toward the lower right; change heights to move it up.
legend_col <- plot_spacer() / legend / plot_spacer() +
plot_layout(heights = c(1.6, 1, 0.3))
p <- (map | legend_col) +
plot_layout(widths = c(1, 0.14)) +
plot_annotation(
title = "Double jeopardy: shrinking and aging",
subtitle = subtitle,
caption = sprintf("Source: U.S. Census 2000 & Vintage %d estimates; %d age structure. Broomfield CO = 2000 city; VA cities merged.",
pop_vintage, age_vintage),
theme = theme(
text = element_text(family = font),
plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
plot.subtitle = element_text(size = 9.5, colour = "grey30", hjust = 0.5,
lineheight = 1.2, margin = margin(4, 0, 6, 0)),
plot.caption = element_text(size = 7.5, colour = "grey40", hjust = 1),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(12, 12, 8, 12)))
dev <- if (requireNamespace("ragg", quietly = TRUE)) ragg::agg_png else "png"
ggsave(out_png, p, width = 13, height = 8, dpi = 300, bg = "white", device = dev)
message("Saved ", out_png, " and ", out_csv)
print(p)
# ---------------------------------------------------------------------------
# Notes
# * 2000 is the April 1 census count, not a July 1 estimate: tidycensus
# doesn't reach the 2000-2010 intercensal estimates. The difference is
# ~1% of growth and flips only counties within a hair of zero change.
# * Ties (o65 == u15) count as Younger; pop2025 == pop2000 counts as Growing.
# * Vintage 2025 PEP support needs a recent tidycensus
# (install.packages("tidycensus") or remotes::install_github("walkerke/tidycensus")).
# * Expect counts near, not exactly on, the original's 1,233 / 1,035.
# ---------------------------------------------------------------------------