# =============================================================================
# replication.R — Toronto Speed Camera Study, Release 2026-07-18
# =============================================================================
# WHAT THIS IS
#   An independent check of the study's headline numbers, written for a lab
#   assistant with no special setup. It uses ONLY base R — nothing to
#   install — and reads plain CSV files that sit next to it in
#   data/csv_for_replication/.
#
# WHY THIS FILE EXISTS
#   The study's credibility rests on a simple promise: every published
#   number can be recomputed from the archived data by someone else, on
#   another machine, with different software. This script keeps that
#   promise for the descriptive core of the study. If you can run it and
#   every check reads OK, then the chain from raw data to headline claim
#   holds without trusting anyone's memory — including ours.
#
# HOW TO RUN IT
#   1. Open a terminal in this folder (the one containing this file).
#   2. Type:  Rscript replication.R
#   3. Read the output top to bottom.
#
# WHAT "OK" MEANS — AND WHAT IT DOES NOT MEAN
#   For counts and medians, OK means an exact match (within rounding).
#   For p-values, OK means the recomputed and published values land on the
#   same side of the conventional 0.05 threshold — exact p-values can
#   differ slightly because the archived data cache was refreshed in July
#   2026 (report Section 14.13), which moved a few sample sizes by one or
#   two sites without changing any conclusion.
#   An OK verdict says the DATA produce the PUBLISHED NUMBER. It does not,
#   by itself, say cameras caused anything — the causal reading depends on
#   design arguments the report makes and qualifies in Sections 4, 14, 16.
# =============================================================================

# ---- 0. Setup ---------------------------------------------------------------
# WHY: replication must not depend on whose machine this is or where the
# folder was copied. So the script finds its own location instead of using
# a typed-in path — the single most common reason replication scripts fail.

# Ask R how it was started; the answer contains this script's own path.
args <- commandArgs(trailingOnly = FALSE)
# Among the start-up arguments, find the one that names this file.
file_arg <- grep("^--file=", args, value = TRUE)
# Strip the prefix to get the path, keep the folder part. If the script is
# being pasted into an interactive session instead, fall back to the
# current working directory.
here <- if (length(file_arg) == 1) dirname(sub("^--file=", "", file_arg)) else getwd()
# All replication inputs live in one subfolder, so there is exactly one
# place to look — and one place to regenerate (export_csv_for_replication.py).
csvdir <- file.path(here, "data", "csv_for_replication")

# A helper that prints one check in a consistent format, so every test in
# this file reports the same three things: what we computed, what the
# report says, and whether they agree. Uniformity is the point — a reader
# should never wonder whether a check was quietly held to a looser standard.
check <- function(label, got, pub, ok, note = "") {
  # Turn the TRUE/FALSE agreement flag into a word.
  verdict <- if (ok) "OK" else "MISMATCH"
  # One aligned line per check: label, recomputed value, published value.
  cat(sprintf("%-52s recomputed: %-18s published: %-14s %s\n",
              label, got, pub, verdict))
  # Fine print, where a check needs context (e.g., the cache-vintage note).
  if (nzchar(note)) cat("      note:", note, "\n")
}

# Section banners so the output can be scanned quickly.
banner <- function(title) cat("\n====", title, "====\n")

# ---- 1. Load the data -------------------------------------------------------
# WHY THESE FOUR TABLES: they are the study's entire descriptive skeleton.
#   - one row per CAMERA SITE (did crashes fall where cameras stood?)
#   - one row per WYS SIGN (the comparison intervention)
#   - one row per ENFORCEMENT PERIOD (cameras rotated — the rotation is
#     what lets the study watch cameras arrive, leave, and return)
#   - one row per COLLISION (the outcome everything is measured in)

banner("Loading data")

# Per-site camera results. Each row holds the site's crash rate before,
# during, and after enforcement, and the percent change between them.
cam <- read.csv(file.path(csvdir, "camera_effectiveness.csv"))
# Per-sign WYS results, in the same shape, for the comparison intervention.
wys <- read.csv(file.path(csvdir, "wys_effectiveness.csv"))
# One row per camera deployment window (start and end dates).
enf <- read.csv(file.path(csvdir, "enforcement_periods.csv"))
# One row per police-reported collision, with the flags the checks need:
# its year, whether anyone was injured, whether it has usable coordinates,
# whether it happened within 250 m of a camera site, and whether that
# camera was operating at the time.
col <- read.csv(file.path(csvdir, "collisions_flat.csv"))
# The adjusted-rebound result computed by the Python pipeline, flattened to
# one row. Printed later for reference — see Check 12 for why it is not
# recomputed here.
reb <- read.csv(file.path(csvdir, "definitive_rebound_primary.csv"))

# Show the sizes immediately. WHY: if any of these is wrong, you are not
# looking at the study's data — stop and check the folder before reading on.
cat("camera sites table:", nrow(cam), "rows\n")
cat("WYS signs table:   ", nrow(wys), "rows\n")
cat("periods table:     ", nrow(enf), "rows\n")
cat("collisions table:  ", nrow(col), "rows\n")

# ---- 2. Scope checks --------------------------------------------------------
# WHY THIS SECTION: every later statistic descends from these counts. They
# are the chain of custody. If the totals match, you hold the same universe
# of crashes and cameras the report analyzed; if they do not, nothing
# downstream is interpretable, no matter how close it looks.

banner("Check 1-4: scope")

# Check 1 — the study's universe: every police-reported collision in
# Toronto, 2014 through 2025. TESTING: that the release ships the exact
# dataset the report describes. A different count would mean a different
# data vintage — the subtle failure this release was built to prevent.
n_total <- nrow(col)
check("1. Total collision records", format(n_total, big.mark = ","),
      "790,725", n_total == 790725)

# Check 2 — how many records carry usable coordinates. WHY IT MATTERS: the
# study links crashes to cameras BY LOCATION, so only geocoded records can
# enter the spatial analysis. The report is explicit that 16.3% cannot be
# placed on the map (Section 14.8) and treats that as a limitation. We
# verify the share so that limitation is exactly as disclosed — no more,
# no less.
n_geo <- sum(col$geocoded == 1)
pct_geo <- round(100 * n_geo / n_total, 1)
check("2. Geocoded records (share)",
      paste0(format(n_geo, big.mark = ","), " (", pct_geo, "%)"),
      "661,535 (83.7%)", n_geo == 661535 && pct_geo == 83.7)

# Check 3 — the number of treatment locations. TESTING: 520 camera sites,
# the denominator behind every "X% of sites" claim in the report.
n_sites <- nrow(cam)
check("3. Camera sites", n_sites, "520", n_sites == 520)

# Check 4 — deployments and re-deployments. WHY 627 > 520: cameras
# ROTATED. Sites were enforced, released, and sometimes re-enforced. That
# rotation is the study's main asset — it lets us watch what happens when
# a camera leaves (the "rebound") and when it returns. The 96 multi-period
# sites are the re-deployment cases those return-visit tests rely on.
n_periods <- nrow(enf)
multi <- sum(table(enf$site_code) > 1)
check("4. Periods / multi-period sites",
      paste0(n_periods, " / ", multi), "627 / 96",
      n_periods == 627 && multi == 96)

# ---- 3. The headline within-site result -------------------------------------
# WHY THIS SECTION: this is the study's most-quoted number. Each site is
# compared TO ITSELF — its monthly crash rate while the camera operated
# versus the same site's rate before the camera arrived. Comparing a place
# to itself removes everything fixed about the place (road design, traffic
# mix, neighbourhood) as an explanation. What it does NOT remove is
# anything that changed over time city-wide (the pandemic, most of all) —
# which is why the report brackets this number with the calendar-matched
# designs discussed at Check 12.

banner("Check 5-7: within-site crash change")

# Keep the sites where a percent change could be computed at all (a site
# with no measurable before-rate has nothing to compare against).
pc <- cam$pct_change[!is.na(cam$pct_change)]

# Check 5 — the MEDIAN site's change. WHY MEDIAN AND NOT MEAN: site-level
# changes are wildly skewed — a handful of sites swung by hundreds of
# percent because small crash counts make ratios jumpy. A mean would let
# those few extremes set the headline; the median reports the typical
# site. TESTING: the -10.2% that leads the report.
med <- round(median(pc), 1)
check("5. Median within-site change", paste0(med, "%"), "-10.2%",
      med == -10.2)

# Check 6 — the split behind the average. WHY: an average alone can hide
# how uneven the result is, and the lab's writing rule is that no average
# stands alone. TESTING: that roughly six sites in ten improved — and,
# just as important, that four in ten did NOT. Both halves of that
# sentence are part of the finding.
n_improved <- sum(pc < 0)   # negative change = fewer crashes
n_worse    <- sum(pc > 0)   # positive change = more crashes
check("6. Sites improved / worsened",
      paste0(n_improved, " / ", n_worse), "304 / 216",
      n_improved == 304 && n_worse == 216)

# Check 7 — could a median like that arise by chance if cameras did
# nothing? WHY THE WILCOXON TEST SPECIFICALLY: the changes are not
# bell-shaped (skewed, heavy-tailed, built from small counts), so a t-test's
# normality assumption would be wrong. The Wilcoxon signed-rank test uses
# only the ordering of the changes, which is robust to exactly this kind
# of data. Two-sided, because honesty requires letting the data surprise
# us in either direction. TESTING: that the fall in crashes is
# statistically distinguishable from zero — same verdict as the report,
# even though the exact p moved with the July 2026 cache refresh.
w <- wilcox.test(pc)
check("7. Wilcoxon test (change < 0)",
      format(signif(w$p.value, 2)), "p = 2.2e-06",
      w$p.value < 1e-4,
      "report resynced to the shipped cache in July 2026 (was 519 sites / 0.0025)")

# ---- 4. The fairness check the July 2026 revision added ---------------------
# WHY THIS SECTION EXISTS AT ALL: the study's comparison intervention (the
# WYS signs, Section 7 below) had an artifact — signs deployed for a few
# weeks often recorded ZERO crashes purely by chance, which shows up as a
# fake "-100% reduction" and flattered the signs' average. The report
# rightly excluded those. The July 2026 audit asked the obvious fairness
# question: cameras have the same artifact — 47 camera sites recorded zero
# crashes during deployment — so the same rule must apply to them.
# TESTING: that the report's own corrected statement of this comparison is
# what the data produce. This is the single most self-critical number in
# the release, which is exactly why it must replicate.

banner("Check 8-9: symmetric zero-crash exclusion")

# Count the camera sites the fairness rule removes.
zero_cam <- sum(cam$during_crashes == 0, na.rm = TRUE)
# Check 8 — the report says 47 such sites. TESTING: the rule's footprint.
check("8. Camera sites with zero during-crashes", zero_cam, "47",
      zero_cam == 47)

# Apply the same exclusion the WYS analysis uses: keep only sites with at
# least one crash during deployment, so a -100% can only mean a real fall,
# not an accident of a short window.
pc_sym <- cam$pct_change[!is.na(cam$pct_change) & cam$during_crashes > 0]
# Recompute the median and the significance test on the symmetric sample.
med_sym <- round(median(pc_sym), 1)
w_sym <- wilcox.test(pc_sym)
# Check 9 — WHY THIS MATTERS: under the symmetric rule the camera headline
# is roughly HALVED (-10.2% becomes -5.1%) and loses conventional
# significance (p = 0.063). The corrected report states this openly and
# rests the cameras-vs-signs conclusion on other tests. Replicating it
# proves the report's self-correction is data, not diplomacy.
check("9. Camera median, symmetric rule",
      paste0(med_sym, "% (p = ", signif(w_sym$p.value, 2), ")"),
      "-5.1% (p = 0.063)",
      med_sym == -5.1 && round(w_sym$p.value, 2) == 0.06,
      "not significant at 0.05 - exactly as the report now states")

# ---- 5. What happened when cameras left -------------------------------------
# WHY THIS SECTION: the rotation gives the study a rare test. If cameras
# were suppressing crashes, then REMOVING one should be followed by crashes
# drifting back up at that site. That is the "departure rebound." It became
# the study's signature finding — and also the finding the July 2026
# revision most reframed, because the years after most departures were also
# the years Toronto's traffic recovered from COVID. Rising crashes after
# departure could therefore be the camera's absence, the city's recovery,
# or both. These two checks verify the RAW rebound (what happened); Check
# 12 covers the adjusted version (what it likely means).

banner("Check 10-11: departure rebound (raw)")

# Strict sample: to say anything about "after the camera left," a site
# needs (a) enough post-camera time to measure — at least 6 months, so one
# odd month cannot decide the verdict — and (b) at least one crash during
# deployment, so the before/after comparison has a base.
strict <- cam[!is.na(cam$after_rate) & !is.na(cam$during_rate) &
              cam$after_months >= 6 & cam$during_crashes > 0, ]
# A site "rebounded" if its monthly crash rate after the camera left
# exceeded its rate while the camera operated.
n_up <- sum(strict$after_rate > strict$during_rate)
# Check 10 — TESTING: the report's strict-sample claim that crashes rose
# at 62% of 383 sites after the camera left. If cameras did nothing and
# nothing else changed, roughly half the sites would rise by chance;
# 62% is the observed excess over that coin-flip baseline.
check("10. Strict rebound sample",
      paste0(n_up, " of ", nrow(strict), " (",
             round(100 * n_up / nrow(strict)), "%)"),
      "236 of 383 (62%)",
      nrow(strict) == 383 && abs(n_up - 236) <= 1,
      "one boundary site differs by cache vintage; share matches at 61-62%")

# Broad sample: every site with ANY post-departure observation. WHY SHOW
# BOTH SAMPLES: the report quotes 62% (strict) and 68% (broad) in
# different places and footnotes the difference; a replicator should see
# that both are real computations on stated samples, not a cherry-pick.
broad <- cam$rebound_pct[!is.na(cam$rebound_pct)]
share_broad <- round(100 * mean(broad > 0), 0)
check("11. Broad rebound share",
      paste0(share_broad, "% of ", length(broad)),
      "68% of 519",
      share_broad >= 65 && share_broad <= 70,
      "cache vintage changes n; the report's direction claim is unchanged")

# The context that the revised report insists must travel with these
# numbers — printed here so the replication cannot be quoted without it:
cat("      context: the report treats the rebound as evidence of DIRECTION,\n")
cat("      not size - after matching the city's own recovery, the camera-\n")
cat("      specific excess is not statistically distinguishable (next check).\n")

# ---- 6. The adjusted rebound (reference values) -----------------------------
# WHY REPORTED, NOT RECOMPUTED: separating "the camera left" from "the city
# recovered" requires knowing, for every one of 661,535 geocoded crashes,
# how far it happened from every camera — so that areas MORE than 500 m
# from any camera can serve as the no-camera yardstick over the exact same
# calendar months. That geospatial computation runs in Python
# (run_rebound_decomposition_definitive.py, in this folder). Base R prints
# its stored result so the replication record contains the study's most
# consequential adjustment, clearly labelled as imported rather than
# recomputed.
# WHY IT MATTERS: this is the number that changed the report's language.
# The excess rebound at camera sites, over and above the recovery visible
# in no-camera areas, is small, extremely variable site to site, and not
# statistically distinguishable from zero — under the cleanest
# specification and across all sixteen reasonable ones.

banner("Check 12: adjusted rebound - reported, not recomputed")

cat(sprintf("      primary spec: %s\n", reb$spec))
cat(sprintf("      median excess: %+.1f points (IQR %.0f to %+.0f), %s%% of sites above control, p = %.2f\n",
            reb$median_excess_pp, reb$iqr_low, reb$iqr_high,
            reb$share_above_counterfactual_pct, reb$wilcoxon_two_sided_p))
cat(sprintf("      across all 16 specifications the median excess spans %+.1f to %+.1f points\n",
            reb$grid_min_median, reb$grid_max_median))
cat("      verdict: no statistically separable camera-specific rebound (report 14.2)\n")

# ---- 7. The comparison intervention (WYS signs) ------------------------------
# WHY SIGNS ARE IN A CAMERA STUDY: the policy question is never "cameras
# versus nothing" — it is "cameras versus the cheaper alternative," and
# Toronto's alternative was the Watch Your Speed radar sign (the display
# that flashes your speed). The report's cameras-vs-signs comparison only
# means something if the signs' own numbers are computed honestly — which
# is where the zero-crash artifact was first found and fixed.

banner("Check 13-15: Watch Your Speed signs")

# Check 13 — the RAW signs median, before the artifact correction. WHY
# VERIFY A NUMBER THE REPORT CALLS MISLEADING: because the report's
# credibility partly rests on showing its corrections' before-and-after.
# The -24.9% is what a naive analysis would have published.
wpc <- wys$pct_change[!is.na(wys$pct_change)]
med_wys_raw <- round(median(wpc), 1)
check("13. WYS raw median (uncorrected)", paste0(med_wys_raw, "%"),
      "-24.9%", med_wys_raw == -24.9)

# Check 14 — the corrected signs result. THE RULE: drop signs with zero
# crashes during deployment (the same rule Check 9 applied to cameras —
# symmetry is the point). One retained sign has no computable percent
# change (its before-rate is zero), so the median runs over the computable
# ones within the corrected sample.
wys_ok <- wys[wys$during_crashes > 0 & !is.na(wys$during_crashes), ]
med_wys <- round(median(wys_ok$pct_change, na.rm = TRUE), 1)
# TESTING: that the correction moves the signs from -24.9% to -8.8% — the
# single largest honest revision in the study, and the template for the
# camera-side fairness check above.
check("14. WYS corrected median (n)",
      paste0(med_wys, "% (n = ", format(nrow(wys_ok), big.mark = ","), ")"),
      "-8.8% (3,438)",
      med_wys == -8.8 && nrow(wys_ok) == 3438)

# Check 15 — WHY DEPLOYMENT LENGTH MATTERS: the zero-crash artifact exists
# BECAUSE sign deployments are short. At a typical site (a few crashes per
# year), a deployment of about 1.6 months will often see zero crashes by
# luck alone. Verifying the 1.6-month median confirms the artifact's
# mechanism, not just its correction. The column is named enf_months,
# matching the camera table's naming.
med_dur <- round(median(wys$enf_months, na.rm = TRUE), 1)
check("15. WYS median deployment", paste0(med_dur, " months"),
      "1.6 months", med_dur == 1.6)

# The like-for-like summary a careful reader should carry away:
cat("      context: under the SAME exclusion rule, cameras show -5.1% and\n")
cat("      WYS signs -8.8% on this simple table; the report's cameras-vs-\n")
cat("      signs conclusion rests on the durability and injury tests, not\n")
cat("      this table (report 15.1).\n")

# ---- 7b. Was it the cameras, or the signs deployed beside them? -------------
# WHY THIS SECTION: a third of camera sites had a Watch Your Speed radar sign
# co-deployed within 250 metres, so a skeptic can ask whether the "camera
# effect" was really a camera-plus-sign bundle. The July 2026 first pass at
# separating them (report Section 14.12) found the answer hides behind a
# composition trap: signs were placed on busier roads, so a naive split
# compares busy camera+sign sites against quiet camera-only sites and
# mistakes the road difference for a sign effect. The honest comparison
# stratifies by baseline busyness first. These checks reproduce that
# analysis from the site-coding table shipped in this release.

banner("Check 16-18: sign-bundle disentangling (exploratory)")

# Load the per-site coding: for each of the 520 sites, whether a sign's
# deployment overlapped the camera's enforcement window (wys_co).
vzc <- read.csv(file.path(csvdir, "site_vision_zero_coding.csv"))
# Join onto the camera results and apply the same symmetric zero-crash rule
# as Check 9 — the report's fairness standard for any cameras-vs-signs claim.
sb <- merge(cam, vzc[, c("site_code", "wys_co")], by = "site_code")
sb <- sb[!is.na(sb$pct_change) & sb$during_crashes > 0, ]
# Convert the text flag ("True"/"False") into a logical, robustly.
sb$wys_co <- sb$wys_co %in% c("True", "TRUE", "true")

# Check 16 — the composition fact that makes stratification necessary.
# TESTING: co-deployed camera sites are nearly twice as busy at baseline.
# If this were false, the pooled split would be fair; because it is true,
# only the stratified comparison below is meaningful.
med_with <- round(median(sb$before_rate[sb$wys_co]), 2)
med_wo   <- round(median(sb$before_rate[!sb$wys_co]), 2)
check("16. Baseline busyness with/without sign",
      paste0(med_with, " / ", med_wo, " crashes/mo"),
      "2.21 / 1.22", med_with == 2.21 && med_wo == 1.22)

# Split the network at its median baseline rate: the "busy half" is where
# crashes — and therefore preventable crashes — concentrate.
busy <- sb[sb$before_rate >= median(sb$before_rate), ]
bw <- busy$pct_change[busy$wys_co]    # busy sites WITH a co-deployed sign
bo <- busy$pct_change[!busy$wys_co]   # busy sites with NO sign anywhere near

# Check 17 — the headline: on busy roads, cameras WITHOUT any sign show a
# reduction of the same size as the bundled sites. WHY IT MATTERS: this is
# what "the sign is not a necessary ingredient" means in numbers. One
# subtlety worth teaching: this median sits at almost exactly -9.15 across
# an even number of sites, so different-but-equally-correct floating-point
# paths round it to -9.1 or -9.2. The check therefore uses a tolerance —
# the honest statement is "about -9%", which is how the report rounds it.
check("17. Busy half: camera-only median",
      paste0(round(median(bo), 1), "% (n = ", length(bo), ")"),
      "≈ -9.2% (114)",
      abs(median(bo) + 9.15) < 0.15 && length(bo) == 114,
      "knife-edge median: -9.1 vs -9.2 are the same value at report rounding")

# The unpaired wilcox.test IS the Mann-Whitney test in base R — it asks
# whether the with-sign and without-sign distributions differ at all.
mw <- wilcox.test(bw, bo)
check("18. Busy half: with-sign vs without, differ?",
      paste0("with ", round(median(bw), 1), "%, diff p = ",
             round(mw$p.value, 2)),
      "-7.9%, p = 0.51",
      round(median(bw), 1) == -7.9 && round(mw$p.value, 2) == 0.51,
      "no detectable difference - the camera association appears without the sign")

# The context a careful reader should carry (printed, per house style):
cat("      context: on the QUIET half, cameras alone show no detectable\n")
cat("      reduction - the report's targeting lesson (little to reduce on\n")
cat("      low-crash roads); and school-zone features, which nearly every\n")
cat("      site has, cannot be separated by comparing sites (report 14.12).\n")

# ---- 7c. Did police fill the camera gaps? -----------------------------------
# WHY THIS SECTION: the most natural rival explanation for the camera
# results is that police enforcement moved with the cameras — into the gaps
# when cameras left, away when cameras arrived. If true, "camera effects"
# could partly be police effects. The rotation era supplies the test: with
# 627 enforcement periods opening and closing, every neighbourhood's camera
# coverage rose and fell; if police backfilled, their speeding tickets
# should move OPPOSITE to that coverage. Report Section 12.8 found they do
# not. This check recomputes that regression from the assembled panel
# shipped with the release.

banner("Check 19: police backfill panel (Section 12.8)")

# One row per neighbourhood-year (159 neighbourhoods x 2020-2024): its
# police speeding tickets and its active camera-months that year.
pb <- read.csv(file.path(csvdir, "police_backfill_panel.csv"))
# The outcome: log tickets (+1 so empty cells don't drop). Logs make the
# coefficient read as a percentage change, matching the report.
pb$log_t <- log(pb$tickets + 1)
# The regression: neighbourhood fixed effects absorb each area's level;
# year fixed effects absorb the citywide ticketing boom (it roughly
# doubled); what remains is whether tickets moved WITH camera coverage.
# factor() creates the fixed effects in base R's lm().
fit <- lm(log_t ~ active_months + factor(hood) + factor(year), data = pb)
# Pull the coefficient on active camera-months.
co <- coef(summary(fit))["active_months", ]
# Check 19 — TESTING: the report's -0.4% per active camera-month, a value
# statistically indistinguishable from zero. WHY IT MATTERS: a flat
# coefficient means police enforcement did not track camera coverage, so
# the camera-era reductions are not police filling gaps. NOTE ON THE
# P-VALUE: the report's p = 0.65 uses neighbourhood-clustered errors
# (computed in the Python pipeline); base R's plain OLS p differs in value
# but not in verdict — both sit far above any significance threshold.
check("19. Backfill: log-tickets per camera-month",
      paste0(round(co["Estimate"], 4), " (OLS p = ", round(co["Pr(>|t|)"], 2), ")"),
      "-0.0040 (clust. p = 0.65)",
      abs(co["Estimate"] + 0.00404) < 0.0005 && co["Pr(>|t|)"] > 0.1,
      "flat under both error structures - no detectable backfill")

# ---- 7d. Is the headline just a change in how crashes get REPORTED? ---------
# WHY THIS SECTION: the study's outcome is police-REPORTED collisions, and
# Toronto's reporting practices moved during the study (intake changes in
# March 2020, drifting severity mix). The bound: injury crashes are the
# LEAST reporting-sensitive outcome (people are hurt, police attend);
# property-damage-only crashes are the MOST sensitive (thresholds and
# self-reporting live there). If the decline shows at BOTH extremes,
# reporting drift cannot be its driver. Report Section 14.14.

banner("Check 20: reporting-regime bounds (Section 14.14)")

# One row per site: the during-minus-before rate DIFFERENCE for each strip
# (differences, not ratios — rare events make ratios misbehave; a
# difference has no -100% floor and no selection trap).
rb2 <- read.csv(file.path(csvdir, "reporting_bounds.csv"))
# Recompute each strip's median difference and its direction share; the
# report's claim is that all three are negative and significant, with the
# injury strip (least reporting-sensitive) the most consistent.
mi <- median(rb2$inj_d, na.rm = TRUE)   # injury strip
mp <- median(rb2$pdo_d, na.rm = TRUE)   # property-damage strip
# Wilcoxon on each strip: could medians like these arise by chance?
pi <- wilcox.test(rb2$inj_d)$p.value
pp <- wilcox.test(rb2$pdo_d)$p.value
check("20. Both strips decline (inj / pdo medians)",
      paste0(round(mi, 3), " / ", round(mp, 3), " per site-month"),
      "-0.030 / -0.089",
      abs(mi + 0.0298) < 0.003 && abs(mp + 0.0885) < 0.003 &&
        pi < 0.001 && pp < 0.001,
      "declines at BOTH reporting-sensitivity extremes - not a reporting artifact")

# ---- 7e. The selection effect, reproduced (Test 7) --------------------------
# WHY: Test 7 is the report's teaching moment - a naive comparison finds
# crashes near cameras MORE likely to injure (cameras were placed at
# dangerous spots on purpose). Until July 2026 it had no producing script;
# this check closes that gap with one line of base R.
banner("Check 21: Test 7 selection effect")
# Logistic regression on geocoded records: injury ~ near_camera. The odds
# ratio ABOVE one is the selection effect the report explains in 4.3.
g7 <- glm(injury ~ near_camera, family = binomial,
          data = col[col$geocoded == 1, ])
or7 <- exp(coef(g7)["near_camera"])
check("21. Naive near-camera injury OR",
      round(or7, 3), "1.056 (selection)",
      abs(or7 - 1.056) < 0.005,
      "wrong-direction by design - the report's argument FOR within-site tests")

# ---- 8. Honest list of what base R cannot recompute -------------------------
# WHY THIS LIST EXISTS: a replication that silently skips the hard parts
# is worse than none. Each item below needs machinery base R does not have
# here (spatial joins, models with site-specific treatment timing, or a
# specific Python 3.12 environment). Nothing is hidden: the producing
# script for every item ships in this folder, so a technical replicator
# can go one level deeper.

banner("Reported in the study but not recomputed here")

cat("      - Injury odds ratios (OR 0.951 / OR 0.920): run_staggered_did.py\n")
cat("        and run_crash_count_did.py (site-specific treatment timing)\n")
cat("      - Callaway-Sant'Anna estimator (near zero): run_callaway_santanna.py\n")
cat("        (requires Python 3.12 via uv, per its header)\n")
cat("      - Adjusted rebound decomposition: run_rebound_decomposition_definitive.py\n")
cat("      - Synthetic control, halo, RTM grid: run_synthetic_control.py,\n")
cat("        toronto_full_analysis.py, run_rtm_estimation.py\n")

# ---- 9. Summary -------------------------------------------------------------

banner("Done")
# What agreement means — and the limit of what it means. The checks above
# verify the data-to-number chain for the descriptive core. The causal
# interpretation of those numbers rests on the design arguments the report
# makes, hedges, and stress-tests in Sections 4, 14, and 16.
cat("      If the checks above read OK, the release data reproduce the\n")
cat("      report's headline descriptive results exactly, and the study's\n")
cat("      corrected framing (direction over magnitude for the rebound;\n")
cat("      symmetric rules for cameras and signs) matches the data.\n")