diff options
Diffstat (limited to '2026')
| -rw-r--r-- | 2026/2026-02-21-pandemic-paradox.R | 263 |
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" + ) +) + + + + + + + + + + + + + + + + + + + + + + + + |