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