add: pandemic traffic fatalities

1 files changed, 263 insertions, 0 deletions

diff --git a/2026/2026-02-21-pandemic-paradox.R b/2026/2026-02-21-pandemic-paradox.R
new file mode 100644
index 0000000..b4a6e61
--- /dev/null
+++ b/2026/2026-02-21-pandemic-paradox.R
@@ -0,0 +1,263 @@
+library(tidyverse)
+library(stringi)
+library(hrbrthemes) # it's on sourcehut, codeberg, & (ugh) github
+library(tidycensus)
+library(lme4)
+library(lmerTest)
+library(ggtext)
+library(grid)
+library(gtable)
+
+# ingest and enrich! ------------------------------------------------------
+
+### NOTE ###
+# You will likely need to get your own curl since I'm not bothering
+# to figure out if the cookie "expires" or not.
+
+r"(
+curl \
+  --url 'https://cdan.dot.gov/SASJobExecution/' \
+  --cookie $'sas-ingress-nginx=31f69ac20404464b36c87be9bc93b2ee|4149f60c1217d8889a1d2f7717687849; eb4b13ae-059e-4c4e-86c0-ac97a3c6a64b=MTc3MTY2OTQ2NnxEWDhFQVFMX2dBQUJEQUVNQUFBcF80QUFBUU4xY21rZ0wxTkJVMHB2WWtWNFpXTjFkR2x2Ymk4X2MzTnZYMmQxWlhOMFBYUnlkV1U9fLPpx1XJW8Uux-80R0q7_osNqAg-OhaSmds7iXDeFwrG; JSESSIONID=2ed667f7-76a7-45b6-ae88-b700ad4c8aac; 2ed667f7-76a7-45b6-ae88-b700ad4c8aac=MTc3MTY2OTkyNXxOd3dBTkV0S1NrUkNNbFZRUkVKT01qSlBRVWxOU0ZwV1RFOVVRVmRCV0V0S1VscE9SRlJhVVVwR00xWXlNMDFHVEZGWk0wMUhTMEU9fI9lJDLP04nVQzC1ian8kzKDtUCFJVegsL3D1kbtIRUS; ak_bmsc=4C6092472E95F0A25B89C1880E772754~000000000000000000000000000000~YAAQnZQZuKYZ2GecAQAAsp+6fx6/3te2YSVq4dN/bx0rdqjfyFAPWqwGJsbiCXJ1GBWcajgQcRM5KYebtPC8uSpZ+erSJThjdbtmlAghJluQQEW/Mysp36uwBvwC4boJ5n9egsNgzn3wdqQG2MXe7WLftPcsa/P+W7D/730yLJ5D8amYBvpw78GXr3AG9MylU61rcByZTF2fSaTTjZhfOCuLoieK4M+EwWersCZQ++M2/Yvwx1HpuxZquZ0b/G4u2kQ4MSZowNkxO6wX9+PnqQD3g6v3OWXLbt0mLCtf193Pz8wXI4PfoJrLua5KDnf+JX7DDAzTA3hja1loyE87zAJm9wTue9ddq77vRYbF9eZ7Ic0sjqByK4HqUW1HFJ5AC2Ded/hkp6ISC2fHz7/eJ2L7OaMRBM+iJ89HQWqHwZRRUw==; NHTSA=\u0021eKlbc2T2GyI7o+bI3sRetaVg0a8n9EVtxMnMLRs7MLFoSGUIX8xSs9KT1kSwv4GrVDzvbvZ6JL72OrQ=; bm_sv=EE12702B2B24BD62FAC0CB636C242FB3~YAAQbabQF/H4bVucAQAA0Ku/fx7fwxilxzzPNR8vMgQ2nW+1lWbzn5H60WJMoPftu4j2g/hbg0jgTdz3B5na9+d3/40gkXdWTGQfCiRemWkB/gnsMRg0JvBd5H8QLiP52UA58iz5FOdfs/lVwg8zsytu04+3dbaYjiJv/EMbqyTzFDeN8Vk2rUJO6hCwgzYEMx581x4jdc1ZP7jx+lqru2Er7EtPYeT4aX1k2/ol2wpBgRVeiA7vGSJkWlPHZA==~1' \
+  --output '/Users/hrbrmstr/Data/stsi/{FIPS}.xlsx' \
+  --header 'accept: text/html,application/xhtml+xml,application/xml;q=0.9,image/avif,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3;q=0.7' \
+  --header 'accept-language: en-US,en;q=0.9' \
+  --header 'cache-control: max-age=0' \
+  --header 'content-type: application/x-www-form-urlencoded' \
+  --header 'origin: https://cdan.dot.gov' \
+  --header 'priority: u=0, i' \
+  --header 'referer: https://cdan.dot.gov/SASJobExecution/' \
+  --header 'sec-ch-ua: "Chromium";v="145", "Not:A-Brand";v="99"' \
+  --header 'sec-ch-ua-mobile: ?0' \
+  --header 'sec-ch-ua-platform: "macOS"' \
+  --header 'sec-fetch-dest: iframe' \
+  --header 'sec-fetch-mode: navigate' \
+  --header 'sec-fetch-site: same-origin' \
+  --header 'sec-fetch-user: ?1' \
+  --header 'upgrade-insecure-requests: 1' \
+  --header 'user-agent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/145.0.0.0 Safari/537.36' \
+  --data-raw '_PROGRAM=%2FPublic%2FOTRA%2FApps%2FSTSI%2FSTSIReport&STATE={FIPS}&COUNTY=NOT_SELECTED&OUTPUTDEST=EXCEL&ExcelPGM=1'
+)" -> template
+
+walk(1:56, \(FIPS) {
+  system(
+    glue::glue(template)
+  )
+})
+
+fips_codes |>
+  as_tibble() |>
+  distinct(state_name, state_code) |>
+  mutate(
+    numeric_code = as.numeric(state_code)
+  ) |>
+  filter(
+    numeric_code <= 56
+  ) -> fips_xlate
+
+# see if states are missing and ingest the ones that aren't
+
+system("file /Users/hrbrmstr/Data/stsi/*", intern=TRUE) |>
+  stri_split_fixed(": ", simplify = TRUE) |>
+  as.data.frame.matrix() |>
+  setNames(c("file", "type")) |>
+  filter(
+    stri_detect_fixed(type, "Excel")
+  ) |>
+  select(-type) |>
+  mutate(
+    numeric_code = basename(file) |>
+      tools::file_path_sans_ext() |>
+      as.numeric()
+  ) |>
+  left_join(
+    fips_xlate,
+    by = "numeric_code"
+  ) |>
+  mutate(
+    data = map(file, readxl::read_excel, range = "C4:L5") |> # we only need one line of counts
+      map(gather, year, fatalities)
+  ) |>
+  unnest(data) |>
+  select(
+    -file,
+    -numeric_code,
+    -state_code
+  ) |>
+  mutate(
+    state_name = ifelse(is.na(state_name), "Puerto Rico", state_name)
+  ) -> stsi_state_fatalities
+
+# stats! ------------------------------------------------------------------
+
+# Pre-pandemic baseline per state (2014-2019)
+stsi_state_fatalities |>
+  filter(
+    year %in% as.character(2014:2019)
+  ) |>
+  group_by(state_name) |>
+  summarise(
+    baseline_mean = mean(fatalities),
+    .groups = "drop"
+  ) -> baseline
+
+# Get 2020 and 2021 values
+stsi_state_fatalities |>
+  filter(
+    year == "2020"
+  ) |>
+  select(
+    state_name,
+    fat_2020 = fatalities
+  ) -> y2020
+
+stsi_state_fatalities |>
+  filter(
+    year == "2021"
+  ) |>
+  select(
+    state_name,
+    fat_2021 = fatalities
+  ) -> y2021
+
+baseline |>
+  left_join(y2020, by = "state_name") |>
+  left_join(y2021, by = "state_name") |>
+  mutate(
+    diff_2020 = fat_2020 - baseline_mean,
+    diff_2021 = fat_2021 - baseline_mean,
+    pct_change_2020 = diff_2020 / baseline_mean * 100,
+    pct_change_2021 = diff_2021 / baseline_mean * 100
+  ) -> comparison
+
+# 👀
+summary(comparison$pct_change_2020)
+summary(comparison$pct_change_2021)
+
+# Paired approach: one-sample t-test on the differences
+# H0: mean difference = 0 (no change from baseline)
+t_2020 <- t.test(comparison$diff_2020, mu = 0, alternative = "less")  # testing for decrease
+t_2021 <- t.test(comparison$diff_2021, mu = 0, alternative = "less")
+
+# Non-parametric alternative (doesn't assume normality of differences)
+w_2020 <- wilcox.test(comparison$diff_2020, mu = 0, alternative = "less")
+w_2021 <- wilcox.test(comparison$diff_2021, mu = 0, alternative = "less")
+
+# 👀
+t_2020
+t_2021
+w_2020
+w_2021
+
+stsi_state_fatalities |>
+  mutate(
+    year_num = as.numeric(year),
+    is_2020 = as.integer(year == "2020"),
+    is_2021 = as.integer(year == "2021"),
+    is_pandemic = as.integer(year %in% c("2020", "2021"))
+  ) -> model_df
+
+# Model with state random intercepts, year trend, and pandemic indicators
+fit <- lmer(fatalities ~ year_num + is_2020 + is_2021 + (1 | state_name), data = model_df)
+
+# 👀
+summary(fit)
+
+# Ttest pandemic as a single block
+fit2 <- lmer(fatalities ~ year_num + is_pandemic + (1 | state_name), data = model_df)
+
+# 👀
+summary(fit2)
+
+# plots! ------------------------------------------------------------------
+
+# we'll order states by most increasing to most decreasing
+
+comparison |>
+  arrange(desc(pct_change_2021)) |>
+  pull(state_name) -> state_order
+
+r"(States ordered by highest increasing fatality count to most decreasing (excluding slight increases in 2021 for 2 states).
+Records-based Studies use Police Accident Reports (PAR) as the basis for all qualifying cases.
+Cases are coded solely from information obtained via the crash documents.
+Data source: <https://cdan.dot.gov/STSI/stsi.htm>
+Not all states report data into STSI.
+Note free Y scales)" -> caption_text
+
+subtitle_text <- r"(In 2020, the median state saw a +6.9% increase over baseline. Mean difference was +54.5 fatalities per state (t-test p = 1.0, Wilcoxon p = 1.0). Only 12 of 52 states saw any decrease at all, and those were mostly small-population states (Hawaii, Alaska, North Dakota, etc.). 2021 was even worse. The median state saw +17.3% increase. Mean difference was +137.6 fatalities per state. p = 1.0 across the board. Only 5 states decreased. <span style='font-weight:bold; color:#BA8E23'>Highlighted region</span> is 2020-2021.)"
+
+stsi_state_fatalities |>
+  mutate(
+    state_name = factor(state_name, levels = state_order),
+    year_as_date = as.Date(sprintf("%s-01-01", year))
+  ) |>
+  ggplot() +
+  annotate(
+    geom = "rect",
+    xmin = as.Date("2020-01-01"),
+    xmax = as.Date("2021-01-01"),
+    ymin = -Inf, ymax = Inf,
+    fill = alpha("gold", 1/3)
+  ) +
+  geom_line(aes(year_as_date, fatalities, group = state_name)) +
+  scale_x_date(
+    breaks = as.Date(c("2014-01-01", "2023-01-01")),
+    date_labels = "%Y"
+  ) +
+  scale_y_comma() +
+  facet_wrap(~state_name, scales = "free_y") +
+  labs(
+    x = NULL, y = NULL,
+    title = "The Pandemic Paradoxically Made Driving Deadlier",
+    subtitle = subtitle_text,
+    caption = NULL  # <-- kill the caption
+  ) +
+  theme_ipsum_gs(grid = "", plot_title_size = 24, subtitle_size = 16, axis_text_size = 8, strip_text_size = 13) +
+  theme(
+    strip.text.x.top = element_textbox_simple(
+      height = unit(2.3, "lines"),
+      margin = margin(b = 2/3, unit="lines"),
+      lineheight = 0.9,
+      valign = 0,
+      vjust = 0
+    ),
+    plot.subtitle = element_textbox_simple(
+      margin = margin(t = 1/2, b = 1.5, unit = "lines")
+    ),
+    panel.spacing.y = unit(1, "lines")
+  )
+
+grid::grid.text(
+  caption_text,
+  x = 0.976, y = 0.08,
+  just = c("right", "center"),
+  gp = gpar(
+    fontsize = 9,
+    col = "grey40",
+    lineheight = 0.9,
+    fontfamily = "Goldman Sans"
+  )
+)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+