add: pandemic traffic fatalities

2 files changed, 264 insertions, 1 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"
+  )
+)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/LICENSE b/LICENSE
index a157816..d8f154b 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2024 hrbrmstr
+Copyright (c) 2024-2026 hrbrmstr
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: