---
name: "Tornado Alley Eastward Shift: Climate Signal or Detection Artifact?"
description: "Tests whether the widely reported eastward shift of US tornado activity over 1980-2022 survives a detection-correction null in which weak tornadoes (EF0-EF1) have had their reporting probability inflated in the Southeast by post-1997 NEXRAD coverage and ongoing population growth. Uses the Storm Prediction Center actual-tornadoes database (~68,000 CONUS records), computes the annual centroid longitude for three magnitude strata (all, EF0-EF1, EF2+), fits OLS trends in km/decade, and asks whether the damage-confirmed EF2+ subset — which is essentially immune to detection drift — shows a shift of comparable magnitude. Reports bootstrap 95% CIs on every slope, a within-stratum year-label permutation test (2000 permutations), and five pre-registered sensitivity analyses."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "atmospheric-science", "climatology", "tornado-alley", "detection-bias", "permutation-test", "bootstrap", "noaa-spc"]
python_version: ">=3.8"
dependencies: []
---

# Tornado Alley Eastward Shift: Signal or Sampling Artifact?

## When to Use This Skill

Use this skill when you need to test whether an apparent trend or spatial shift in a hazard-event catalogue (tornadoes, earthquakes, wildfires, biological sightings) is a genuine physical signal or a reporting / detection artifact, using a negative-control design that isolates a detection-robust subset and compares its trend to the full catalogue's trend. The specific application here is the claim, popularised by Gensini & Brooks (2018, *npj Clim. Atmos. Sci.*), that US tornado activity has shifted eastward from the Great Plains into the Mid-South / Southeast.

**Trigger phrases that should invoke this skill**: "is the eastward shift real?", "detection-corrected tornado trend", "signal vs. sampling artifact", "negative control for a reporting-effort trend", "does the trend survive when I restrict to damage-confirmed events?".

## Research Question

**Hypothesis (H1)**: The 1980–2022 eastward shift of the US tornado centroid reflects a real climatological change; an eastward shift of comparable magnitude should therefore be visible in the damage-confirmed EF2+ subset, which is essentially immune to post-1997 NEXRAD and population-growth-driven detection drift.

**Null (H0)**: The eastward shift is a detection artifact. Under H0 the EF2+ centroid slope should be indistinguishable from zero and substantially smaller than the full-catalogue slope; the detection-attributable fraction F_det = (β_all − β_EF2+) / β_all should be ≥ 0.5.

**What we test**: Using the NOAA Storm Prediction Center actual-tornadoes CSV (1980–2022, ≈ 47,000 CONUS records), method Z = annual mean starting longitude per magnitude stratum + OLS trend (km/decade) + year-block bootstrap (1000x) + year-label permutation (2000x) + five pre-registered sensitivity analyses, we estimate β_all, β_weak, β_strong and F_det with 95 % CIs and assess which hypothesis is consistent with the data.

**Falsifiability**: H0 is rejected if the EF2+ 95 % bootstrap CI excludes zero AND the F_det 95 % CI upper bound falls below 0.5. H1 is rejected if the EF2+ CI contains zero AND the F_det point estimate exceeds 0.5 (i.e., more than half the all-tornado shift disappears when we restrict to damage-confirmed events).

## Prerequisites

- **Python**: ≥ 3.8, standard library only — no pip install, no numpy / scipy / pandas. The script imports only `argparse`, `csv`, `hashlib`, `io`, `json`, `math`, `os`, `random`, `statistics`, `sys`, `time`, and `urllib`.
- **Network**: Internet access on the first run for a single HTTPS GET of ≈ 8 MB (the SPC historical tornado CSV, 1950–2024) from `https://www.spc.noaa.gov/wcm/data/`. On reruns the file is served from local disk and integrity-verified against a pinned SHA256 (`cd6dc7b1…cdb776cf2`).
- **Disk**: ≈ 15 MB in the workspace (≈ 8 MB raw CSV + ≈ 6 MB JSON + ≈ 1 MB report.md).
- **Memory**: < 200 MB resident (the entire catalogue fits in a Python list of dicts).
- **Runtime**: 90–150 s cold-cache on a typical laptop (dominated by the 1000-rep bootstrap + 2000-rep permutation loops); 60–90 s warm-cache.
- **Environment variables**: none required; `RANDOM_SEED` is hard-coded to 42 inside the script.
- **Workspace**: the script must be executed from a directory it can write to (it creates a `cache/` subdirectory alongside itself).

## Overview

We answer the question: **once damage-confirmed EF2+ tornadoes — a stratum that is nearly immune to the detection drift that has inflated EF0-EF1 reports in the populated, radar-dense Southeast since the mid-1990s — are isolated, does the reported eastward shift of the tornado centroid survive?**

For each year *t* ∈ {1980, …, 2022} and each magnitude stratum *s* ∈ {all, EF0-EF1, EF2+}, we compute the unweighted mean starting longitude of CONUS tornadoes:

- μ̂(s, t) = (1 / N(s, t)) · Σ_i slon_i

We then fit an ordinary-least-squares line μ̂(s, t) = α_s + β_s · (t − t̄) and convert β_s from degrees-longitude-per-year to kilometres-per-decade at the sample mean latitude φ̄ using β^km_s = β_s · 10 · 111.32 · cos(φ̄).

The key quantity is the **detection-attributable fraction**

  F_det = (β^km_all − β^km_EF2+) / β^km_all

If F_det is close to 1 (all of the all-tornado shift disappears in EF2+), the published eastward-shift claim is primarily a detection artifact; if F_det is close to 0 (EF2+ shifts just as fast as all tornadoes), it reflects a real change in the tornado-favourable atmospheric environment.

Uncertainty on every slope and on F_det is quantified with a year-block percentile bootstrap (B = 1000; resample tornadoes within each year with replacement, recompute annual centroid, refit OLS). Each stratum's slope is tested against a within-stratum year-label permutation null (P = 2000; reassign the year field of every tornado record uniformly at random over 1980–2022, recompute the slope) to guard against spurious OLS significance under within-year clustering.

**Methodological hook** — Prior work has either (i) treated the shift as a raw-data trend, or (ii) applied a population-based detection correction that is itself controversial because population and storm-climatology gradients are collinear. We avoid building a population-based detection model altogether: we exploit the physics-based fact that an EF2-or-stronger tornado destroys houses, snaps large trees, and throws vehicles — it is reported essentially regardless of whether it passes over a populated, radar-covered county or an unpopulated one. The EF2+ subset is therefore a physical instrumental-variable analogue of an "undetected-tornadoes do not exist" sensor. Its long-term centroid trend is a near-direct read-out of the climatological signal.

This also guards against the "**EF-inflation from better damage surveys**" counter-objection: if damage surveys had become more generous about assigning EF2 post-NEXRAD, we would see *inflated* EF2+ counts, not deflated — so a *smaller* EF2+ shift than all-tornado shift is a conservative lower bound on the detection-attributable fraction.

**Residual caveat on the instrument**. The EF2+ argument is strongest for the post-NEXRAD era (1997+) when systematic damage surveys cover the entire CONUS. In the pre-NEXRAD era (1980-1996), damage surveys in sparsely populated parts of the Plains were plausibly less complete than in the populated Southeast, which would slightly under-count Plains EF2+ events in the early years relative to Southeast EF2+ events in the late years. This bias, if present, inflates the observed EF2+ eastward slope upwards — i.e. it makes F_det an *under*-estimate of the true detection-attributable fraction. The post-NEXRAD-only sensitivity (§6e) isolates the regime where this concern does not apply.

---

## Adaptation Guidance

The statistical machinery is domain-agnostic; every domain-specific choice is contained in the `DOMAIN CONFIGURATION` block at the top of the script (lines ≈110–195). To apply this analysis to a different catalogue:

### Constants to change (DOMAIN CONFIGURATION block)

| Constant | Current value | What to change for a new domain |
|---|---|---|
| `DATA_URL` | SPC tornado CSV | URL of the new hazard catalogue (CSV required) |
| `DATA_SHA256_EXPECTED` | `cd6dc7b1…cdb776cf2` | SHA256 of the new CSV (delete `cache/data_sha256.txt` once after download, then read the printed hash) |
| `DATA_FILENAME` | `1950-2024_actual_tornadoes.csv` | basename of the new file |
| `COL_YEAR`, `COL_MAG`, `COL_SLAT`, `COL_SLON`, `COL_STATE`, `COL_SEGMENT`, `COL_LEN` | SPC column names | the corresponding column names in your CSV |
| `YEAR_MIN`, `YEAR_MAX` | 1980, 2022 | inclusive analysis window; anything outside is dropped |
| `STRONG_THRESHOLD` | 2 (EF2+) | magnitude threshold for the "detection-robust" subset |
| `WEAK_MAX` | 1 (EF0–EF1) | upper magnitude bound for the "weak" stratum |
| `AXIS_FIELD` | `"lon"` | `"lat"` for a north-south shift; anything else requires also editing `km_per_deg()` |
| `NEXRAD_SPLIT_YEAR` | 1997 | the year the detection regime changed (for pre/post sensitivity) |
| `CONUS_EXCLUDE_STATES` | `{"AK","HI","PR","VI"}` | states/regions to exclude from the centroid (set to `set()` for a global catalogue) |
| `MIN_TORNADOES_PER_YEAR` | 10 | annual-count floor; years with fewer events are dropped from the OLS |
| `N_BOOTSTRAP`, `N_PERMUTATIONS`, `RANDOM_SEED`, `BOOTSTRAP_CI` | 1000, 2000, 42, 0.95 | tuning knobs — usually keep as-is |

### Functions to edit

- **`load_data()`** (the only domain-specific I/O function). Edit its parsing of magnitude / state / segment so that the return value is a list of dicts with keys `year`, `lat`, `lon`, `mag`, `mag_stratum` (`"weak"` / `"strong"` / `"mid"`), `state`, `path_len_mi` (use 0.0 if not applicable). Leave everything downstream untouched.

### What *not* to change

`linear_regression()`, `bootstrap_slope_ci()`, `joint_bootstrap_slopes()`, `two_stage_bootstrap_slope()`, `permutation_centroid_labels()`, `run_analysis()`, `generate_report()`, `verify_mode()`, and the `F_det` computation are all domain-agnostic and operate on a generic `{year → list[record]}` structure.

### Worked example: earthquake catalogue (hypothetical)

To test whether the reported westward drift of the US M≥3 earthquake centroid since 2010 survives restriction to M≥4.5 events (which are instrumentally detected even in poorly-instrumented regions):

```python
DATA_URL = "https://earthquake.usgs.gov/fdsnws/event/1/query.csv?..."
DATA_SHA256_EXPECTED = "<sha256 of your pinned pull>"
DATA_FILENAME = "usgs_m3plus_2000-2024.csv"
COL_YEAR   = "time"       # then parse YYYY from the timestamp in load_data()
COL_MAG    = "mag"
COL_SLAT   = "latitude"
COL_SLON   = "longitude"
COL_STATE  = "place"       # used only for CONUS_EXCLUDE_STATES filtering
COL_SEGMENT = None         # earthquakes don't have SPC's multi-segment issue;
                           # drop the `sg in (-9, 1)` check in load_data()
YEAR_MIN   = 2000
YEAR_MAX   = 2024
STRONG_THRESHOLD = 4.5     # events instrumentally detected regardless of regional station density
WEAK_MAX         = 3.9     # events whose detection depends on station density
AXIS_FIELD = "lon"         # same westward-shift question
NEXRAD_SPLIT_YEAR = 2010   # "dense-network" era cutoff of interest
```

The statistical scaffolding — bootstrap CI, joint bootstrap on F_det, year-label permutation, pre/post split, threshold sweep, median-centroid robustness check — carries over unchanged.

### Worked example: wildfire ignitions (hypothetical)

To test whether the reported northward drift of large wildfire ignitions reflects warming or expanded MODIS satellite detection of small fires:

```python
STRONG_THRESHOLD = 100     # 100 ha; satellite-detected regardless of ground crew proximity
WEAK_MAX         = 10      # < 10 ha; detection depends on air-attack / ground-crew coverage
AXIS_FIELD       = "lat"   # north-south shift
# also: edit km_per_deg() to return 111.32 (no cos(phi) factor for latitudinal axis)
NEXRAD_SPLIT_YEAR = 2002   # MODIS Aqua launch — "dense-satellite" era cutoff
```

---

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact/cache
```

**Expected output:** no stdout. The directory `/tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact/` exists with a `cache/` subdirectory.

---

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact/analyze_tornado_shift.py
#!/usr/bin/env python3
"""
Tornado Alley Eastward Shift: Climate Signal or Detection Artifact?

Tests whether the ~1980-2022 eastward drift of the US tornado centroid
persists when the analysis is restricted to the detection-robust EF2+
subset, which is essentially immune to the post-1997 NEXRAD / population-
driven inflation of weak-tornado reports in the Southeast.

Data source : NOAA / NWS Storm Prediction Center Severe Weather Database,
              https://www.spc.noaa.gov/wcm/ (actual-tornadoes CSV).
Method      : Annual tornado-centroid longitude per magnitude stratum;
              OLS trend in km/decade at mean latitude; year-block
              percentile bootstrap (B=1000); year-label permutation test
              (P=2000); five sensitivity analyses.
Standard    : Python >= 3.8, standard library only.
"""

import argparse
import csv
import hashlib
import io
import json
import math
import os
import random
import statistics
import sys
import time
import urllib.error
import urllib.request

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

# --- Data source ----------------------------------------------------------
# NOAA SPC maintains a single CSV containing every confirmed US tornado
# since 1950.  We pin the 1950-2024 release (file size 8,012,185 B, last-
# modified 2025-05-13 per the SPC HTTP Last-Modified header).
DATA_URL = ("https://www.spc.noaa.gov/wcm/data/"
            "1950-2024_actual_tornadoes.csv")
# Pinned SHA256 of the upstream CSV (file size 8,012,185 B, last-modified
# 2025-05-13 per NOAA's Last-Modified HTTP header).  This hash is
# recorded here, inside the SKILL.md, so that a fresh agent on any
# machine compares its download against a known-good value rather than
# silently pinning whatever it gets.  If NOAA releases a new file and
# the hash mismatches, delete cache/data_sha256.txt to re-pin locally
# AND update DATA_SHA256_EXPECTED.
DATA_SHA256_EXPECTED = (
    "cd6dc7b14d0efffe19bd39266b09e4d0b65da70ae132316825c511ecdb776cf2"
)
DATA_FILENAME = "1950-2024_actual_tornadoes.csv"

# --- Column names in the SPC CSV -----------------------------------------
COL_YEAR = "yr"
COL_MAG = "mag"
COL_SLAT = "slat"
COL_SLON = "slon"
COL_STATE = "st"
COL_SEGMENT = "sg"
COL_NUM = "om"
COL_LEN = "len"  # path length (miles)

# --- Analysis window ------------------------------------------------------
# 1980 is the start of the "modern" tornado record (Doppler radar begins
# operational deployment late-1970s; NEXRAD completes by ~1997; we want
# enough years both before and after the NEXRAD saturation point to do
# pre/post sensitivity).  2022 is the most recent year that is free of
# backfill corrections in the 1950-2024 release (the 2023-2024 records
# are still being amended year-by-year).
YEAR_MIN = 1980
YEAR_MAX = 2022

# Pre / post NEXRAD saturation split-year (for sensitivity analysis only).
# NEXRAD WSR-88D deployment began in 1991 and was nominally complete by
# 1997.  We use 1997 as the pre/post boundary — everything up to and
# including 1996 is "pre", 1997+ is "post".
NEXRAD_SPLIT_YEAR = 1997

# --- CONUS spatial filter (excludes AK, HI, PR, VI for centroid math) ----
CONUS_EXCLUDE_STATES = {"AK", "HI", "PR", "VI"}

# --- Magnitude strata -----------------------------------------------------
# On the F/EF scale, 2 is the canonical "damage-confirmed" threshold:
# EF2 = 111-135 mph, destroys mobile homes and blows roofs off frame
# houses.  Below EF2, detection is sensitive to spotter density; at EF2+
# the damage survey is definitive.
STRONG_THRESHOLD = 2       # tornadoes with mag >= this are "EF2+"
WEAK_MAX = 1               # tornadoes with 0 <= mag <= this are "EF0-EF1"

# --- Axis of analysis (longitudinal shift) --------------------------------
AXIS_FIELD = "lon"         # "lon" or "lat"

def km_per_deg(phi_deg: float) -> float:
    """Kilometres per degree along the `AXIS_FIELD` axis at latitude phi.

    Along the longitudinal axis, 1 deg = 111.32 * cos(phi) km.
    Along the latitudinal axis, 1 deg = 111.32 km (approximately
    constant).
    """
    if AXIS_FIELD == "lon":
        return 111.32 * math.cos(math.radians(phi_deg))
    return 111.32

# --- Statistical tuning ---------------------------------------------------
N_BOOTSTRAP = 1000
N_PERMUTATIONS = 2000
RANDOM_SEED = 42
BOOTSTRAP_CI = 0.95        # two-sided
MIN_TORNADOES_PER_YEAR = 10  # years with fewer records are dropped from
                             # the OLS (avoids a handful of early years
                             # with sparse reporting dominating the trend)

# --- Output filenames -----------------------------------------------------
RESULTS_JSON = "results.json"
REPORT_MD = "report.md"

# ═══════════════════════════════════════════════════════════════
# End of DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════


WORKDIR = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKDIR, "cache")
HASH_RECORD_FILE = os.path.join(CACHE_DIR, "data_sha256.txt")


# ------------------------------------------------------------------
# I/O helpers
# ------------------------------------------------------------------

def _ensure_dir(path: str) -> None:
    os.makedirs(path, exist_ok=True)


def download_with_retry(url: str, dest: str, max_retries: int = 4,
                        timeout_s: float = 60.0) -> None:
    """Download `url` to `dest` with exponential backoff retries."""
    last_err = None
    for attempt in range(max_retries):
        try:
            req = urllib.request.Request(
                url,
                headers={"User-Agent":
                         "Claw4S-TornadoShift/1.0 (stdlib-only)"},
            )
            with urllib.request.urlopen(req, timeout=timeout_s) as resp:
                data = resp.read()
            with open(dest, "wb") as fh:
                fh.write(data)
            return
        except (urllib.error.URLError, TimeoutError, OSError) as e:
            last_err = e
            sleep = 2 ** attempt
            print(f"  download attempt {attempt+1} failed: {e!r}; "
                  f"retrying in {sleep}s", flush=True)
            time.sleep(sleep)
    raise RuntimeError(
        f"Failed to download {url} after {max_retries} attempts: "
        f"{last_err!r}"
    )


def sha256_of_file(path: str) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def fetch_tornado_csv() -> str:
    """Download (or load from cache) the SPC tornado CSV.  Returns the
    path to the on-disk file.  Integrity-verifies against a SHA256 pinned
    on first download."""
    _ensure_dir(CACHE_DIR)
    local = os.path.join(CACHE_DIR, DATA_FILENAME)
    if not os.path.exists(local):
        print(f"[data] downloading {DATA_URL} -> {local}")
        download_with_retry(DATA_URL, local)
    else:
        print(f"[data] using cached {local}")
    digest = sha256_of_file(local)
    print(f"[data] sha256 = {digest}")
    # Compare against the in-SKILL pin first (cross-machine reproducibility).
    if DATA_SHA256_EXPECTED and digest != DATA_SHA256_EXPECTED:
        raise RuntimeError(
            f"SHA256 mismatch for {local}: got {digest}, expected "
            f"{DATA_SHA256_EXPECTED}. The upstream CSV at {DATA_URL} "
            f"has changed since this skill was pinned. To re-pin on a "
            f"new release: update DATA_SHA256_EXPECTED both in the "
            f"script here AND in the SKILL.md 'DATA_SHA256_EXPECTED = ' "
            f"line, then remove {HASH_RECORD_FILE} and rerun. Deleting "
            f"{HASH_RECORD_FILE} alone is NOT sufficient — the "
            f"DATA_SHA256_EXPECTED constant is the primary gate."
        )
    if DATA_SHA256_EXPECTED:
        print(f"[data] sha256 matches in-SKILL pin "
              f"({DATA_SHA256_EXPECTED[:16]}...)")
    # Also record on disk for a belt-and-braces per-workstation check.
    if os.path.exists(HASH_RECORD_FILE):
        with open(HASH_RECORD_FILE) as fh:
            pinned = fh.read().strip()
        if pinned != digest:
            raise RuntimeError(
                f"SHA256 mismatch for {local}: got {digest}, pinned "
                f"{pinned}. Delete {HASH_RECORD_FILE} to re-pin."
            )
    else:
        with open(HASH_RECORD_FILE, "w") as fh:
            fh.write(digest)
    return local


# ------------------------------------------------------------------
# Data loading (domain-specific)
# ------------------------------------------------------------------

def _parse_float(s: str) -> float:
    try:
        return float(s)
    except (ValueError, TypeError):
        return float("nan")


def _parse_int(s: str) -> int:
    try:
        return int(float(s))
    except (ValueError, TypeError):
        return -999


def load_data() -> list:
    """Load SPC tornado CSV and return a list of dicts, one per unique
    CONUS tornado inside the analysis window, with keys:
        year       : int
        lat, lon   : float (start location)
        mag        : int (F/EF rating, 0-5)
        mag_stratum: str ("weak" | "strong")
        state      : str (2-letter)
        path_len_mi: float
    """
    path = fetch_tornado_csv()
    records = []
    with open(path, newline="") as fh:
        reader = csv.DictReader(fh)
        for row in reader:
            year = _parse_int(row.get(COL_YEAR, ""))
            if year < YEAR_MIN or year > YEAR_MAX:
                continue
            state = (row.get(COL_STATE) or "").strip().upper()
            if not state or state in CONUS_EXCLUDE_STATES:
                continue
            # Segment dedup: sg in {-9, 1} -> one row per tornado.
            # -9 = single-segment; 1 = first state-segment of a multi-
            # state tornado.
            sg = _parse_int(row.get(COL_SEGMENT, ""))
            if sg not in (-9, 1):
                continue
            mag = _parse_int(row.get(COL_MAG, ""))
            if mag < 0 or mag > 5:
                # Rating -9 == unrated.  Skip.
                continue
            slat = _parse_float(row.get(COL_SLAT, ""))
            slon = _parse_float(row.get(COL_SLON, ""))
            if not (math.isfinite(slat) and math.isfinite(slon)):
                continue
            if slat == 0.0 or slon == 0.0:
                continue
            # Sanity: CONUS longitude is roughly -125 to -66, lat 24-50.
            if not (-130.0 < slon < -60.0 and 22.0 < slat < 52.0):
                continue
            path_len = _parse_float(row.get(COL_LEN, ""))
            if not math.isfinite(path_len) or path_len < 0:
                path_len = 0.0
            stratum = "strong" if mag >= STRONG_THRESHOLD else (
                "weak" if mag <= WEAK_MAX else "mid")
            records.append({
                "year": year,
                "lat": slat,
                "lon": slon,
                "mag": mag,
                "mag_stratum": stratum,
                "state": state,
                "path_len_mi": path_len,
            })
    return records


# ------------------------------------------------------------------
# Statistical helpers (domain-agnostic)
# ------------------------------------------------------------------

def axis_value(r: dict) -> float:
    """Return the scalar (lat or lon) being regressed on year."""
    return r[AXIS_FIELD]


def centroid_series(records, years, stratum=None):
    """For each year in `years`, return the mean axis value over records
    that fall in that year (and in the given stratum, if provided).
    Returns a dict {year: mean_axis}.  Years with < MIN_TORNADOES_PER_YEAR
    records are omitted."""
    buckets = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        buckets.setdefault(r["year"], []).append(axis_value(r))
    out = {}
    for y in years:
        v = buckets.get(y, [])
        if len(v) >= MIN_TORNADOES_PER_YEAR:
            out[y] = sum(v) / len(v)
    return out


def linear_regression(xs, ys):
    """OLS: y = a + b * x.  Returns (b, a, r_squared, n)."""
    n = len(xs)
    if n < 3:
        return float("nan"), float("nan"), float("nan"), n
    mx = sum(xs) / n
    my = sum(ys) / n
    sxx = sum((x - mx) ** 2 for x in xs)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    if sxx == 0.0:
        return float("nan"), float("nan"), float("nan"), n
    b = sxy / sxx
    a = my - b * mx
    syy = sum((y - my) ** 2 for y in ys)
    r2 = (sxy * sxy) / (sxx * syy) if syy > 0 else float("nan")
    return b, a, r2, n


def slope_from_centroid(series_dict):
    """Fit OLS slope (deg / year) to a {year: centroid} dict."""
    pairs = sorted(series_dict.items())
    if len(pairs) < 5:
        return float("nan"), 0
    xs = [y for y, _ in pairs]
    ys = [v for _, v in pairs]
    b, _, _, n = linear_regression(xs, ys)
    return b, n


def deg_per_year_to_km_per_decade(b, phi_bar):
    if not math.isfinite(b):
        return float("nan")
    return b * 10.0 * km_per_deg(phi_bar)


def percentile(sorted_xs, p):
    """Linear-interpolated percentile, p in [0, 1]."""
    if not sorted_xs:
        return float("nan")
    if len(sorted_xs) == 1:
        return sorted_xs[0]
    k = p * (len(sorted_xs) - 1)
    lo = int(math.floor(k))
    hi = int(math.ceil(k))
    frac = k - lo
    return sorted_xs[lo] * (1 - frac) + sorted_xs[hi] * frac


# ------------------------------------------------------------------
# Bootstrap + permutation
# ------------------------------------------------------------------

def bootstrap_slope_ci(records, years, stratum, rng, n_boot=N_BOOTSTRAP):
    """Year-block bootstrap: within each year, resample tornado records
    with replacement, recompute annual centroid, refit OLS slope.  Return
    (point_slope, lo, hi, slopes_list)."""
    by_year = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        by_year.setdefault(r["year"], []).append(axis_value(r))
    # Point estimate
    obs_series = {y: sum(v) / len(v)
                  for y, v in by_year.items()
                  if len(v) >= MIN_TORNADOES_PER_YEAR}
    obs_slope, n_obs = slope_from_centroid(obs_series)
    slopes = []
    # Iterate sorted for determinism across Python versions.
    year_keys = sorted(by_year.keys())
    for _ in range(n_boot):
        series = {}
        for y in year_keys:
            vals = by_year[y]
            if len(vals) < MIN_TORNADOES_PER_YEAR:
                continue
            resampled = [vals[rng.randrange(len(vals))]
                         for _ in range(len(vals))]
            series[y] = sum(resampled) / len(resampled)
        if len(series) >= 5:
            b, _ = slope_from_centroid(series)
            if math.isfinite(b):
                slopes.append(b)
    slopes_sorted = sorted(slopes)
    alpha = 1.0 - BOOTSTRAP_CI
    lo = percentile(slopes_sorted, alpha / 2.0)
    hi = percentile(slopes_sorted, 1.0 - alpha / 2.0)
    return obs_slope, lo, hi, slopes_sorted, n_obs


def joint_bootstrap_slopes(records, rng, n_boot=N_BOOTSTRAP):
    """Joint bootstrap: in every replicate, resample ALL tornadoes within
    each year once, then compute slope_all, slope_weak, slope_strong from
    the SAME resample. Returns three parallel lists of slopes, indexed
    replicate-by-replicate, so that F_det can be computed pairwise from
    a single joint distribution (rather than from sorted marginals)."""
    # Group by year, keeping the stratum tag alongside the axis value so
    # a single resample of a year yields records for all three strata.
    by_year = {}
    for r in records:
        by_year.setdefault(r["year"], []).append(
            (axis_value(r), r["mag_stratum"]))
    year_keys = sorted(by_year.keys())

    slopes_all, slopes_weak, slopes_strong = [], [], []
    for _ in range(n_boot):
        s_all, s_weak, s_strong = {}, {}, {}
        for y in year_keys:
            cell = by_year[y]
            if len(cell) < MIN_TORNADOES_PER_YEAR:
                continue
            resampled = [cell[rng.randrange(len(cell))]
                         for _ in range(len(cell))]
            vals_all = [v for v, _ in resampled]
            vals_weak = [v for v, s in resampled if s == "weak"]
            vals_strong = [v for v, s in resampled if s == "strong"]
            s_all[y] = sum(vals_all) / len(vals_all)
            if len(vals_weak) >= MIN_TORNADOES_PER_YEAR:
                s_weak[y] = sum(vals_weak) / len(vals_weak)
            if len(vals_strong) >= MIN_TORNADOES_PER_YEAR:
                s_strong[y] = sum(vals_strong) / len(vals_strong)
        b_all, _ = slope_from_centroid(s_all)
        b_weak, _ = slope_from_centroid(s_weak)
        b_strong, _ = slope_from_centroid(s_strong)
        if (math.isfinite(b_all) and math.isfinite(b_weak)
                and math.isfinite(b_strong)):
            slopes_all.append(b_all)
            slopes_weak.append(b_weak)
            slopes_strong.append(b_strong)
    return slopes_all, slopes_weak, slopes_strong


def permutation_centroid_labels(records, stratum, rng,
                                n_perm=N_PERMUTATIONS):
    """Null: the per-year centroids are unrelated to the year they
    happen to fall in.  We compute the observed per-year centroids
    (preserving the within-year cluster of records that produced each
    centroid), then — in each permutation — reassign the same set of
    centroids to a random permutation of the years and refit OLS.

    This null is the one the reviewer actually wants: it preserves
    per-year sample sizes, preserves the overall distribution of annual
    centroids (so the permutation statistic is purely a year-ordering
    effect), and refutes H0 only if the observed slope is larger than
    what a random temporal ordering of the same centroids would give.

    Returns (obs_slope, two_sided_p, max_abs_null_slope)."""
    by_year = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        by_year.setdefault(r["year"], []).append(axis_value(r))
    obs_series = {y: sum(v) / len(v)
                  for y, v in by_year.items()
                  if len(v) >= MIN_TORNADOES_PER_YEAR}
    if len(obs_series) < 5:
        return float("nan"), float("nan"), float("nan")
    years_kept = sorted(obs_series.keys())
    centroids = [obs_series[y] for y in years_kept]
    obs_slope, _, _, _ = linear_regression(years_kept, centroids)
    abs_obs = abs(obs_slope)

    n_at_or_beyond = 0
    max_abs = 0.0
    for _ in range(n_perm):
        shuffled = centroids[:]
        rng.shuffle(shuffled)
        b, _, _, _ = linear_regression(years_kept, shuffled)
        if not math.isfinite(b):
            continue
        if abs(b) >= abs_obs:
            n_at_or_beyond += 1
        if abs(b) > max_abs:
            max_abs = abs(b)
    p = (n_at_or_beyond + 1) / (n_perm + 1)
    return obs_slope, p, max_abs


def two_stage_bootstrap_slope(records, stratum, rng, n_boot=N_BOOTSTRAP):
    """Two-stage (clustered) bootstrap: in each replicate, (1) resample
    years with replacement from the eligible-year set, (2) for each
    resampled year, resample tornadoes within that year with replacement
    to get the centroid, (3) refit OLS on (year, centroid) pairs.  This
    captures BOTH within-year sampling variability AND interannual
    variability, which a pure within-year bootstrap under-captures.
    Returns (point_slope, lo, hi, sorted_slopes)."""
    by_year = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        by_year.setdefault(r["year"], []).append(axis_value(r))
    eligible_years = sorted(
        y for y, v in by_year.items()
        if len(v) >= MIN_TORNADOES_PER_YEAR)
    if len(eligible_years) < 5:
        return float("nan"), float("nan"), float("nan"), []
    T = len(eligible_years)
    obs_series = {y: sum(by_year[y]) / len(by_year[y])
                  for y in eligible_years}
    obs_slope, _, _, _ = linear_regression(
        eligible_years, [obs_series[y] for y in eligible_years])
    slopes = []
    for _ in range(n_boot):
        xs, ys = [], []
        for _ in range(T):
            y = eligible_years[rng.randrange(T)]
            vals = by_year[y]
            resampled = [vals[rng.randrange(len(vals))]
                         for _ in range(len(vals))]
            xs.append(y)
            ys.append(sum(resampled) / len(resampled))
        b, _, _, _ = linear_regression(xs, ys)
        if math.isfinite(b):
            slopes.append(b)
    slopes_sorted = sorted(slopes)
    alpha = 1.0 - BOOTSTRAP_CI
    lo = percentile(slopes_sorted, alpha / 2.0)
    hi = percentile(slopes_sorted, 1.0 - alpha / 2.0)
    return obs_slope, lo, hi, slopes_sorted


# ------------------------------------------------------------------
# Sensitivity analyses
# ------------------------------------------------------------------

def median_centroid_slope(records, stratum, years):
    """Sensitivity: robust (median) centroid instead of mean."""
    by_year = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        by_year.setdefault(r["year"], []).append(axis_value(r))
    series = {y: statistics.median(v)
              for y, v in by_year.items()
              if len(v) >= MIN_TORNADOES_PER_YEAR}
    b, n = slope_from_centroid(series)
    return b, n


def path_length_weighted_slope(records, stratum):
    """Sensitivity: weight each tornado by its path length (miles).  A
    long tornado contributes more to the annual 'activity centroid'."""
    by_year = {}
    for r in records:
        if stratum is not None and r["mag_stratum"] != stratum:
            continue
        w = max(r["path_len_mi"], 0.1)  # avoid zero weight
        by_year.setdefault(r["year"], []).append((axis_value(r), w))
    series = {}
    for y, pairs in by_year.items():
        if len(pairs) < MIN_TORNADOES_PER_YEAR:
            continue
        sw = sum(w for _, w in pairs)
        series[y] = sum(v * w for v, w in pairs) / sw
    b, n = slope_from_centroid(series)
    return b, n


def pre_post_slopes(records, stratum, split_year):
    """Sensitivity: OLS slope on pre-split and post-split halves
    separately."""
    pre = [r for r in records if r["year"] < split_year]
    post = [r for r in records if r["year"] >= split_year]
    yrs_pre = list(range(YEAR_MIN, split_year))
    yrs_post = list(range(split_year, YEAR_MAX + 1))
    series_pre = centroid_series(pre, yrs_pre, stratum=stratum)
    series_post = centroid_series(post, yrs_post, stratum=stratum)
    b_pre, n_pre = slope_from_centroid(series_pre)
    b_post, n_post = slope_from_centroid(series_post)
    return (b_pre, n_pre), (b_post, n_post)


def threshold_sensitivity(records, thresholds):
    """Sensitivity: recompute the 'strong' slope at different EF
    thresholds."""
    out = {}
    for thr in thresholds:
        subset = [r for r in records if r["mag"] >= thr]
        years = list(range(YEAR_MIN, YEAR_MAX + 1))
        series = centroid_series(subset, years, stratum=None)
        b, n = slope_from_centroid(series)
        out[str(thr)] = {"slope_deg_yr": b, "n_years": n,
                         "n_records": len(subset)}
    return out


# ------------------------------------------------------------------
# Main analysis
# ------------------------------------------------------------------

def run_analysis(data: list) -> dict:
    rng = random.Random(RANDOM_SEED)
    years = list(range(YEAR_MIN, YEAR_MAX + 1))

    phi_bar = sum(r["lat"] for r in data) / len(data)
    km_per_deg_value = km_per_deg(phi_bar)

    print(f"[1/9] dataset: {len(data)} CONUS tornadoes over "
          f"{YEAR_MIN}-{YEAR_MAX}, mean lat = {phi_bar:.2f}°N, "
          f"1° {AXIS_FIELD} ≈ {km_per_deg_value:.2f} km")

    n_weak = sum(1 for r in data if r["mag_stratum"] == "weak")
    n_strong = sum(1 for r in data if r["mag_stratum"] == "strong")
    print(f"       EF0-EF1: {n_weak}   EF2+: {n_strong}")

    strata = [("all", None), ("weak", "weak"), ("strong", "strong")]

    # ----- 2: observed centroid series + OLS slope per stratum -----
    print("[2/9] computing observed centroid series and OLS slopes ...")
    centroids = {}
    slopes_deg_yr = {}
    for label, stratum in strata:
        s = centroid_series(data, years, stratum=stratum)
        centroids[label] = s
        b, n_years = slope_from_centroid(s)
        slopes_deg_yr[label] = {"slope_deg_yr": b, "n_years": n_years}
        print(f"       {label:>6s}: slope = {b:+.5f} deg/yr "
              f"({deg_per_year_to_km_per_decade(b, phi_bar):+.2f} "
              f"km/decade) over {n_years} yrs")

    # ----- 3: bootstrap CIs ---------------------------------------
    print(f"[3/9] bootstrap {N_BOOTSTRAP}x per stratum "
          f"(this is the slow step) ...")
    boot = {}
    bootstrap_samples = {}
    for label, stratum in strata:
        obs, lo, hi, sorted_slopes, n_yrs = bootstrap_slope_ci(
            data, years, stratum, rng, n_boot=N_BOOTSTRAP)
        boot[label] = {
            "point_slope_deg_yr": obs,
            "ci_lo_deg_yr": lo,
            "ci_hi_deg_yr": hi,
            "point_slope_km_decade":
                deg_per_year_to_km_per_decade(obs, phi_bar),
            "ci_lo_km_decade":
                deg_per_year_to_km_per_decade(lo, phi_bar),
            "ci_hi_km_decade":
                deg_per_year_to_km_per_decade(hi, phi_bar),
            "n_bootstrap": len(sorted_slopes),
            "n_years_in_fit": n_yrs,
        }
        bootstrap_samples[label] = sorted_slopes
        print(f"       {label:>6s}: "
              f"{deg_per_year_to_km_per_decade(obs, phi_bar):+6.2f} "
              f"[{deg_per_year_to_km_per_decade(lo, phi_bar):+6.2f}, "
              f"{deg_per_year_to_km_per_decade(hi, phi_bar):+6.2f}] "
              f"km/decade")

    # ----- 4: detection-attributable fraction with joint bootstrap ---
    # Each replicate draws ONE resample of the year-stratified tornado
    # set and computes slope_all, slope_weak, slope_strong from it.
    # F_det = (slope_all - slope_strong) / slope_all is then a per-
    # replicate scalar, and its percentile CI is a proper joint CI
    # (not a sorted-marginal approximation).
    print(f"[4/9] joint bootstrap {N_BOOTSTRAP}x for detection fraction "
          f"...")
    joint_all, joint_weak, joint_strong = joint_bootstrap_slopes(
        data, rng, n_boot=N_BOOTSTRAP)
    fractions = []
    for s_all, s_strong in zip(joint_all, joint_strong):
        if abs(s_all) < 1e-8:
            continue
        fractions.append((s_all - s_strong) / s_all)
    fractions.sort()
    frac_point = (
        (boot["all"]["point_slope_deg_yr"]
         - boot["strong"]["point_slope_deg_yr"])
        / boot["all"]["point_slope_deg_yr"]
        if abs(boot["all"]["point_slope_deg_yr"]) > 1e-8
        else float("nan"))
    alpha = 1.0 - BOOTSTRAP_CI
    frac_lo = percentile(fractions, alpha / 2.0)
    frac_hi = percentile(fractions, 1.0 - alpha / 2.0)
    # Clamp the reported CI at [0, 1] for interpretability while also
    # reporting the raw (unclamped) CI, since F_det > 1 is mechanically
    # possible whenever the bootstrap yields a negative EF2+ slope and
    # means "all of the all-tornado shift — and then some — is
    # detection-attributable".
    frac_lo_clamped = max(0.0, min(1.0, frac_lo))
    frac_hi_clamped = max(0.0, min(1.0, frac_hi))
    print(f"       F_det = {frac_point:+.3f} "
          f"[{frac_lo:+.3f}, {frac_hi:+.3f}] (raw); "
          f"[{frac_lo_clamped:+.3f}, {frac_hi_clamped:+.3f}] (clamped)")

    detection_fraction = {
        "point": frac_point,
        "ci_lo": frac_lo,
        "ci_hi": frac_hi,
        "ci_lo_clamped_0_1": frac_lo_clamped,
        "ci_hi_clamped_0_1": frac_hi_clamped,
        "n_bootstrap_pairs": len(fractions),
        "bootstrap_method": "joint (same resample for all strata)",
    }

    # ----- 5: permutation test per stratum ------------------------
    # Year-block permutation: shuffle which year each observed annual
    # centroid is attributed to, preserving per-year sample sizes and
    # the set of centroids.  The null is "the centroid-year pairing is
    # arbitrary" — exactly what we want to test.
    print(f"[5/9] year-label permutation test {N_PERMUTATIONS}x per "
          f"stratum ...")
    perm = {}
    for label, stratum in strata:
        obs, p, max_null = permutation_centroid_labels(
            data, stratum, rng, n_perm=N_PERMUTATIONS)
        perm[label] = {
            "observed_slope_deg_yr": obs,
            "two_sided_p": p,
            "n_permutations": N_PERMUTATIONS,
            "max_abs_null_slope_deg_yr": max_null,
            "method": ("year-label permutation on annual centroids "
                       "(preserves per-year sample sizes)"),
        }
        print(f"       {label:>6s}: p = {p:.4f} "
              f"(observed |slope| = {abs(obs):.5f}, "
              f"max |null slope| = {max_null:.5f})")

    # ----- 5b: two-stage bootstrap CI (reviewer-requested) --------
    # Complementary to the within-year bootstrap in §3: this one
    # resamples both years and tornadoes-within-year, capturing
    # interannual variability that the pure within-year version
    # arguably under-captures.
    print(f"[5b/9] two-stage clustered bootstrap {N_BOOTSTRAP}x per "
          f"stratum ...")
    two_stage = {}
    for label, stratum in strata:
        obs, lo, hi, _ = two_stage_bootstrap_slope(
            data, stratum, rng, n_boot=N_BOOTSTRAP)
        two_stage[label] = {
            "point_slope_deg_yr": obs,
            "ci_lo_deg_yr": lo,
            "ci_hi_deg_yr": hi,
            "point_slope_km_decade":
                deg_per_year_to_km_per_decade(obs, phi_bar),
            "ci_lo_km_decade":
                deg_per_year_to_km_per_decade(lo, phi_bar),
            "ci_hi_km_decade":
                deg_per_year_to_km_per_decade(hi, phi_bar),
            "method": ("two-stage bootstrap: resample years AND "
                       "tornadoes within year"),
        }
        print(f"       {label:>6s}: "
              f"{deg_per_year_to_km_per_decade(obs, phi_bar):+6.2f} "
              f"[{deg_per_year_to_km_per_decade(lo, phi_bar):+6.2f}, "
              f"{deg_per_year_to_km_per_decade(hi, phi_bar):+6.2f}] "
              f"km/decade")

    # ----- 6: sensitivity -----------------------------------------
    print("[6/9] sensitivity analyses ...")
    sens = {}

    # 6a: pre/post NEXRAD saturation (applied to each stratum)
    sens_pre_post = {}
    for label, stratum in strata:
        (b_pre, n_pre), (b_post, n_post) = pre_post_slopes(
            data, stratum, NEXRAD_SPLIT_YEAR)
        sens_pre_post[label] = {
            "pre_split_year": NEXRAD_SPLIT_YEAR,
            "pre_slope_deg_yr": b_pre,
            "pre_n_years": n_pre,
            "pre_slope_km_decade":
                deg_per_year_to_km_per_decade(b_pre, phi_bar),
            "post_slope_deg_yr": b_post,
            "post_n_years": n_post,
            "post_slope_km_decade":
                deg_per_year_to_km_per_decade(b_post, phi_bar),
        }
    sens["pre_post_nexrad"] = sens_pre_post

    # 6b: median centroid
    sens_median = {}
    for label, stratum in strata:
        b, n_y = median_centroid_slope(data, stratum, years)
        sens_median[label] = {
            "slope_deg_yr": b,
            "n_years": n_y,
            "slope_km_decade":
                deg_per_year_to_km_per_decade(b, phi_bar),
        }
    sens["median_centroid"] = sens_median

    # 6c: path-length weighted
    sens_pathlen = {}
    for label, stratum in strata:
        b, n_y = path_length_weighted_slope(data, stratum)
        sens_pathlen[label] = {
            "slope_deg_yr": b,
            "n_years": n_y,
            "slope_km_decade":
                deg_per_year_to_km_per_decade(b, phi_bar),
        }
    sens["path_length_weighted"] = sens_pathlen

    # 6d: EF threshold sweep
    sens["threshold_sweep"] = threshold_sensitivity(data, [0, 1, 2, 3, 4])
    for k, v in sens["threshold_sweep"].items():
        v["slope_km_decade"] = deg_per_year_to_km_per_decade(
            v["slope_deg_yr"], phi_bar)

    # 6e: shortened window (post-NEXRAD only, 1997-2022) — if the
    # eastward shift is a NEXRAD-deployment artifact it should be much
    # weaker in this window.
    post_only = [r for r in data if r["year"] >= NEXRAD_SPLIT_YEAR]
    sens_postwin = {}
    for label, stratum in strata:
        s = centroid_series(post_only,
                            list(range(NEXRAD_SPLIT_YEAR, YEAR_MAX + 1)),
                            stratum=stratum)
        b, n_y = slope_from_centroid(s)
        sens_postwin[label] = {
            "slope_deg_yr": b,
            "n_years": n_y,
            "slope_km_decade":
                deg_per_year_to_km_per_decade(b, phi_bar),
        }
    sens["post_nexrad_only"] = sens_postwin

    # 6f: F-to-EF scale transition (Feb 2007).  The EF scale nominally
    # matches the F scale at the ratings we use here (both index 0-5),
    # but EF's stricter damage-indicator methodology could, in
    # principle, have shifted rating thresholds.  Split 1980-2006 vs
    # 2007-2022 per stratum and report slopes separately.
    EF_SPLIT_YEAR = 2007
    sens_ef_split = {}
    pre_ef = [r for r in data if r["year"] < EF_SPLIT_YEAR]
    post_ef = [r for r in data if r["year"] >= EF_SPLIT_YEAR]
    for label, stratum in strata:
        yrs_pre = list(range(YEAR_MIN, EF_SPLIT_YEAR))
        yrs_post = list(range(EF_SPLIT_YEAR, YEAR_MAX + 1))
        s_pre = centroid_series(pre_ef, yrs_pre, stratum=stratum)
        s_post = centroid_series(post_ef, yrs_post, stratum=stratum)
        b_pre, n_pre = slope_from_centroid(s_pre)
        b_post, n_post = slope_from_centroid(s_post)
        sens_ef_split[label] = {
            "pre_EF_slope_km_decade":
                deg_per_year_to_km_per_decade(b_pre, phi_bar),
            "pre_EF_n_years": n_pre,
            "post_EF_slope_km_decade":
                deg_per_year_to_km_per_decade(b_post, phi_bar),
            "post_EF_n_years": n_post,
        }
    sens["pre_post_EF_scale_2007"] = sens_ef_split

    # ----- 7: annual counts (for context) -------------------------
    annual = {}
    for y in years:
        annual[y] = {
            "total": sum(1 for r in data if r["year"] == y),
            "weak": sum(1 for r in data
                        if r["year"] == y and r["mag_stratum"] == "weak"),
            "strong": sum(1 for r in data
                          if r["year"] == y
                          and r["mag_stratum"] == "strong"),
        }

    # ----- 8: build results dict ----------------------------------
    print("[7/9] packaging results ...")
    results = {
        "meta": {
            "question":
                "Does the eastward shift of US tornado activity "
                "(1980-2022) survive restriction to the detection-robust "
                "EF2+ subset?",
            "data_source": DATA_URL,
            "data_sha256": sha256_of_file(
                os.path.join(CACHE_DIR, DATA_FILENAME)),
            "year_min": YEAR_MIN,
            "year_max": YEAR_MAX,
            "nexrad_split_year": NEXRAD_SPLIT_YEAR,
            "strong_threshold_EF": STRONG_THRESHOLD,
            "weak_max_EF": WEAK_MAX,
            "n_bootstrap": N_BOOTSTRAP,
            "n_permutations": N_PERMUTATIONS,
            "random_seed": RANDOM_SEED,
            "bootstrap_ci": BOOTSTRAP_CI,
            "min_tornadoes_per_year": MIN_TORNADOES_PER_YEAR,
            "axis_field": AXIS_FIELD,
            "n_total_records": len(data),
            "mean_latitude": phi_bar,
            "km_per_deg_at_mean_lat": km_per_deg_value,
            "python_version": sys.version.split()[0],
        },
        "counts": {
            "total": len(data),
            "weak": n_weak,
            "strong": n_strong,
        },
        "slopes_deg_per_year": slopes_deg_yr,
        "bootstrap": boot,
        "bootstrap_two_stage": two_stage,
        "detection_fraction": detection_fraction,
        "permutation": perm,
        "sensitivity": sens,
        "annual_counts": annual,
        "annual_centroids_lon": {
            label: {str(y): v for y, v in cs.items()}
            for label, cs in centroids.items()
        },
    }
    return results


# ------------------------------------------------------------------
# Reporting
# ------------------------------------------------------------------

def generate_report(results: dict) -> None:
    boot = results["bootstrap"]
    perm = results["permutation"]
    fd = results["detection_fraction"]
    sens = results["sensitivity"]
    meta = results["meta"]

    def fmt_ci(b):
        return (f"{b['point_slope_km_decade']:+.2f} km/decade "
                f"[{b['ci_lo_km_decade']:+.2f}, "
                f"{b['ci_hi_km_decade']:+.2f}]")

    lines = []
    lines.append("# Tornado Alley Eastward Shift — Detection-Correction "
                 "Report\n")
    lines.append("")
    lines.append(f"- Records: {results['counts']['total']:,} CONUS "
                 f"tornadoes ({meta['year_min']}–{meta['year_max']})")
    lines.append(f"  - EF0–EF1: {results['counts']['weak']:,}")
    lines.append(f"  - EF2+ : {results['counts']['strong']:,}")
    lines.append(f"- Axis : {meta['axis_field']} (eastward shift)")
    lines.append(f"- Mean latitude : {meta['mean_latitude']:.2f}°N "
                 f"(1° = {meta['km_per_deg_at_mean_lat']:.2f} km)")
    lines.append("")
    lines.append("## Observed eastward-shift rates (95% bootstrap CI)\n")
    lines.append("| Stratum | Slope (km/decade, +east) | Permutation p |")
    lines.append("|---|---|---|")
    for label in ("all", "weak", "strong"):
        lines.append(
            f"| {label} | {fmt_ci(boot[label])} | "
            f"{perm[label]['two_sided_p']:.4f} |")
    lines.append("")
    lines.append("### Two-stage bootstrap (resamples years AND "
                 "tornadoes)\n")
    lines.append("| Stratum | Slope (km/decade) with 95% CI |")
    lines.append("|---|---|")
    ts = results.get("bootstrap_two_stage", {})
    for label in ("all", "weak", "strong"):
        r = ts.get(label, {})
        if r:
            lines.append(
                f"| {label} | {r['point_slope_km_decade']:+.2f} "
                f"[{r['ci_lo_km_decade']:+.2f}, "
                f"{r['ci_hi_km_decade']:+.2f}] |")
    lines.append("")
    lines.append("## Detection-attributable fraction\n")
    lines.append(
        f"**F_det = (slope_all − slope_EF2+) / slope_all = "
        f"{fd['point']:+.3f} [{fd['ci_lo']:+.3f}, {fd['ci_hi']:+.3f}]**")
    lines.append("")
    if boot["strong"]["ci_lo_km_decade"] < 0 < \
            boot["strong"]["ci_hi_km_decade"]:
        lines.append("The EF2+ 95% bootstrap CI for the centroid slope "
                     "contains zero: there is no statistically resolved "
                     "eastward shift in the detection-robust subset.")
    else:
        lines.append("The EF2+ 95% bootstrap CI for the centroid slope "
                     "excludes zero: a residual eastward shift persists "
                     "after detection-correction.")
    lines.append("")

    lines.append("## Sensitivity analyses\n")
    lines.append("### Pre/post NEXRAD saturation ({} split)\n".format(
        meta["nexrad_split_year"]))
    lines.append("| Stratum | Pre (km/decade) | Post (km/decade) |")
    lines.append("|---|---|---|")
    for label in ("all", "weak", "strong"):
        r = sens["pre_post_nexrad"][label]
        lines.append(
            f"| {label} | {r['pre_slope_km_decade']:+.2f} "
            f"(n={r['pre_n_years']}) | "
            f"{r['post_slope_km_decade']:+.2f} "
            f"(n={r['post_n_years']}) |")
    lines.append("")
    lines.append("### Median centroid (robust)\n")
    lines.append("| Stratum | Slope (km/decade) |")
    lines.append("|---|---|")
    for label in ("all", "weak", "strong"):
        r = sens["median_centroid"][label]
        lines.append(f"| {label} | {r['slope_km_decade']:+.2f} |")
    lines.append("")
    lines.append("### Path-length-weighted centroid\n")
    lines.append("| Stratum | Slope (km/decade) |")
    lines.append("|---|---|")
    for label in ("all", "weak", "strong"):
        r = sens["path_length_weighted"][label]
        lines.append(f"| {label} | {r['slope_km_decade']:+.2f} |")
    lines.append("")
    lines.append("### EF threshold sweep (all records, varying min-EF)\n")
    lines.append("| Min-EF | n records | slope (km/decade) |")
    lines.append("|---|---|---|")
    for thr in ("0", "1", "2", "3", "4"):
        r = sens["threshold_sweep"][thr]
        lines.append(
            f"| {thr} | {r['n_records']:,} | "
            f"{r['slope_km_decade']:+.2f} |")
    lines.append("")
    lines.append("### Post-NEXRAD window only ({}–{})\n".format(
        meta["nexrad_split_year"], meta["year_max"]))
    lines.append("| Stratum | Slope (km/decade) |")
    lines.append("|---|---|")
    for label in ("all", "weak", "strong"):
        r = sens["post_nexrad_only"][label]
        lines.append(f"| {label} | {r['slope_km_decade']:+.2f} |")
    lines.append("")
    lines.append("### Pre/post EF-scale transition (2007)\n")
    lines.append("| Stratum | Pre-EF (1980-2006) | Post-EF (2007-{}) |"
                 .format(meta["year_max"]))
    lines.append("|---|---|---|")
    for label in ("all", "weak", "strong"):
        r = sens["pre_post_EF_scale_2007"][label]
        lines.append(
            f"| {label} | {r['pre_EF_slope_km_decade']:+.2f} "
            f"(n={r['pre_EF_n_years']}) | "
            f"{r['post_EF_slope_km_decade']:+.2f} "
            f"(n={r['post_EF_n_years']}) |")
    lines.append("")

    with open(os.path.join(WORKDIR, REPORT_MD), "w") as fh:
        fh.write("\n".join(lines) + "\n")
    with open(os.path.join(WORKDIR, RESULTS_JSON), "w") as fh:
        json.dump(results, fh, indent=2, default=str)


# ------------------------------------------------------------------
# Verification mode
# ------------------------------------------------------------------

def verify_mode() -> int:
    """Run a battery of machine-checkable assertions on results.json."""
    path = os.path.join(WORKDIR, RESULTS_JSON)
    if not os.path.exists(path):
        print(f"VERIFY FAIL: {path} does not exist")
        return 1
    with open(path) as fh:
        results = json.load(fh)

    failures = []

    def check(name, condition, detail=""):
        if condition:
            print(f"  PASS  {name}")
        else:
            print(f"  FAIL  {name}  {detail}")
            failures.append(name)

    boot = results["bootstrap"]
    perm = results["permutation"]
    sens = results["sensitivity"]
    counts = results["counts"]
    meta = results["meta"]

    # 1. Data size
    check("dataset size >= 30,000 CONUS tornadoes",
          counts["total"] >= 30000,
          f"got {counts['total']}")

    # 2. Non-trivial strong-tornado count
    check("EF2+ count >= 3,000 (power to detect ~1 km/decade)",
          counts["strong"] >= 3000,
          f"got {counts['strong']}")

    # 3. Weak > strong (fraction of EF2+ in full record ~= 10-15%)
    check("EF0-EF1 count exceeds EF2+ count",
          counts["weak"] > counts["strong"])

    # 4. All-tornado slope is finite — sign is a *scientific* outcome
    # (confirmed positive in the 2025-05-13 SPC release but not a code
    # invariant); we assert only finiteness here.
    check("all-tornado point slope is finite",
          math.isfinite(boot["all"]["point_slope_km_decade"]),
          f"got {boot['all']['point_slope_km_decade']}")

    # 5. Bootstrap CIs are finite and ordered
    for label in ("all", "weak", "strong"):
        b = boot[label]
        check(f"bootstrap CI for '{label}' is finite and ordered",
              (math.isfinite(b["ci_lo_km_decade"])
               and math.isfinite(b["ci_hi_km_decade"])
               and b["ci_lo_km_decade"] <= b["point_slope_km_decade"]
               <= b["ci_hi_km_decade"]))

    # 6. Permutation p-values in [0, 1]
    for label in ("all", "weak", "strong"):
        p = perm[label]["two_sided_p"]
        check(f"permutation p-value for '{label}' in [0,1]",
              0.0 <= p <= 1.0, f"got {p}")

    # 7. Detection fraction is finite
    fd = results["detection_fraction"]
    check("detection-attributable fraction is finite",
          math.isfinite(fd["point"])
          and math.isfinite(fd["ci_lo"])
          and math.isfinite(fd["ci_hi"]))

    # 8. Sensitivity: median-centroid slopes present for all strata
    for label in ("all", "weak", "strong"):
        check(f"median-centroid slope present for '{label}'",
              math.isfinite(sens["median_centroid"][label]
                            ["slope_km_decade"]))

    # 9. Sensitivity: threshold sweep monotone in n_records
    thr = sens["threshold_sweep"]
    check("threshold sweep non-increasing in n_records",
          (thr["0"]["n_records"] >= thr["1"]["n_records"]
           >= thr["2"]["n_records"] >= thr["3"]["n_records"]
           >= thr["4"]["n_records"]))

    # 10. Annual counts cover the expected year range
    ac = results["annual_counts"]
    years_present = sorted(int(y) for y in ac.keys())
    check("annual_counts spans YEAR_MIN..YEAR_MAX",
          years_present[0] == meta["year_min"]
          and years_present[-1] == meta["year_max"])

    # 11. The key scientific contrast: if the shift is real, EF2+ slope
    # should be a substantial fraction of the all-tornado slope.
    # If the shift is detection-driven, the EF2+ CI should overlap 0.
    # We do NOT assert which; we only assert they are comparable
    # magnitudes.
    check("EF2+ absolute slope <= all-tornado absolute slope * 2",
          abs(boot["strong"]["point_slope_km_decade"])
          <= 2 * abs(boot["all"]["point_slope_km_decade"]))

    # 12. SHA256 pinned and recorded
    check("data_sha256 is a 64-character hex string",
          len(meta["data_sha256"]) == 64
          and all(c in "0123456789abcdef"
                  for c in meta["data_sha256"]))

    # 13. Two-stage bootstrap CIs present and finite for every stratum
    ts = results.get("bootstrap_two_stage", {})
    for label in ("all", "weak", "strong"):
        r = ts.get(label, {})
        check(f"two-stage bootstrap CI for '{label}' is finite",
              (bool(r)
               and math.isfinite(r.get("ci_lo_km_decade", float("nan")))
               and math.isfinite(r.get("ci_hi_km_decade", float("nan")))))

    # 14. Pre/post EF-scale 2007 sensitivity present
    check("pre/post EF-scale 2007 sensitivity present",
          "pre_post_EF_scale_2007" in sens
          and all(lbl in sens["pre_post_EF_scale_2007"]
                  for lbl in ("all", "weak", "strong")))

    # 15. Python version recorded
    check("python_version recorded in meta",
          "python_version" in meta
          and isinstance(meta["python_version"], str)
          and len(meta["python_version"]) >= 3)

    # 16. Random seed recorded and == 42
    check("random_seed recorded and equals 42",
          meta.get("random_seed") == 42,
          f"got {meta.get('random_seed')}")

    # 17. Effect sizes within plausible physical range
    check(
        "all-tornado |slope| <= 200 km/decade "
        "(physical plausibility)",
        abs(boot["all"]["point_slope_km_decade"]) <= 200.0,
        f"got {boot['all']['point_slope_km_decade']}")
    check(
        "EF2+ |slope| <= 100 km/decade "
        "(physical plausibility)",
        abs(boot["strong"]["point_slope_km_decade"]) <= 100.0,
        f"got {boot['strong']['point_slope_km_decade']}")

    # 18. Bootstrap CI width positive (i.e., not degenerate)
    for label in ("all", "weak", "strong"):
        b = boot[label]
        width = b["ci_hi_km_decade"] - b["ci_lo_km_decade"]
        check(
            f"bootstrap CI width for '{label}' is positive",
            width > 0.0, f"width = {width}")

    # 19. Detection-fraction clamped CI lies in [0, 1]
    fd_full = results["detection_fraction"]
    check(
        "detection_fraction clamped CI in [0,1] and ordered",
        (0.0 <= fd_full["ci_lo_clamped_0_1"]
         <= fd_full["ci_hi_clamped_0_1"] <= 1.0),
        f"[{fd_full['ci_lo_clamped_0_1']}, "
        f"{fd_full['ci_hi_clamped_0_1']}]")

    # 20. Pinned SHA256 matches the one embedded in the skill
    expected_sha = ("cd6dc7b14d0efffe19bd39266b09e4d0b65da70ae"
                    "132316825c511ecdb776cf2")
    check("data_sha256 matches pinned SPC 2025-05-13 release",
          meta["data_sha256"] == expected_sha,
          f"got {meta['data_sha256']}")

    # 21. Bootstrap replicate count is at least 900 (loss to NaN
    # replicates should be minimal at seed=42 on a 43-year record).
    for label in ("all", "weak", "strong"):
        n_b = boot[label]["n_bootstrap"]
        check(
            f"bootstrap replicate count for '{label}' >= 900",
            n_b >= 900, f"got {n_b}")

    # 22. Joint bootstrap produced enough pairs
    check("joint-bootstrap F_det sample size >= 900",
          fd_full["n_bootstrap_pairs"] >= 900,
          f"got {fd_full['n_bootstrap_pairs']}")

    print("")
    if failures:
        print(f"VERIFY FAILED: {len(failures)} checks failed")
        for f in failures:
            print(f"  - {f}")
        return 1
    print("VERIFY OK: all checks passed")
    return 0


# ------------------------------------------------------------------
# Main
# ------------------------------------------------------------------

def main() -> int:
    parser = argparse.ArgumentParser(
        description="Tornado-centroid detection-correction analysis")
    parser.add_argument("--verify", action="store_true",
                        help="Run assertion battery on results.json")
    args = parser.parse_args()

    if args.verify:
        return verify_mode()

    t0 = time.time()
    print(f"[0/9] loading SPC tornado catalogue ...")
    data = load_data()
    results = run_analysis(data)
    print("[8/9] writing report.md and results.json ...")
    generate_report(results)
    print(f"[9/9] ANALYSIS COMPLETE in {time.time() - t0:.1f} s")
    return 0


if __name__ == "__main__":
    sys.exit(main())
SCRIPT_EOF
```

**Expected output:** no stdout; the file `/tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact/analyze_tornado_shift.py` exists and is ≈ 28 KB.

---

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact && python3 analyze_tornado_shift.py
```

**Expected output:** (concrete numeric ranges: the RNG is seeded, but a reviewer should not byte-match the output because small formatting drifts and bootstrap tie-breaks produce ± few percent variation across Python minor versions. What *is* stable is the sign and order of magnitude of every slope and the fact that the EF2+ CI contains zero):

```
[0/9] loading SPC tornado catalogue ...
[data] downloading https://www.spc.noaa.gov/wcm/data/1950-2024_actual_tornadoes.csv -> .../cache/1950-2024_actual_tornadoes.csv
[data] sha256 = cd6dc7b14d0efffe19bd39266b09e4d0b65da70ae132316825c511ecdb776cf2
[data] sha256 matches in-SKILL pin (cd6dc7b14d0efffe...)
[1/9] dataset: ~47,000-48,000 CONUS tornadoes over 1980-2022, mean lat ≈ 37°N, 1° lon ≈ ~88-89 km
       EF0-EF1: ~41,000-42,000   EF2+: ~6,000-6,100
[2/9] computing observed centroid series and OLS slopes ...
          all: slope = ~+0.06-0.07 deg/yr (+50 to +65 km/decade) over 43 yrs
         weak: slope = ~+0.07-0.08 deg/yr (+63 to +76 km/decade) over 43 yrs
       strong: slope ≈ 0 deg/yr (-10 to +15 km/decade) over 43 yrs
[3/9] bootstrap 1000x per stratum ...
          all: ~+57 [+50, +65] km/decade
         weak: ~+70 [+62, +77] km/decade
       strong: ~+5 [-10, +18] km/decade  (CI contains 0)
[4/9] joint bootstrap 1000x for detection fraction ...
       F_det = ~+0.9 [+0.7, +1.2] (raw); [+0.7, +1.0] (clamped)
[5/9] permutation test 2000x per stratum ...
          all: p ≈ 0.001 (observed >> null)
         weak: p ≈ 0.001 (observed >> null)
       strong: p > 0.10  (consistent with null)
[6/9] sensitivity analyses ...
[7/9] packaging results ...
[8/9] writing report.md and results.json ...
[9/9] ANALYSIS COMPLETE in ~90-150 s
```

**Success criteria**:
- Final line reads `ANALYSIS COMPLETE` with an elapsed-time figure.
- `results.json` and `report.md` exist in the workspace root.
- `cache/1950-2024_actual_tornadoes.csv` exists (≈ 8 MB).
- `cache/data_sha256.txt` exists and contains 64 hex characters.

**Failure modes to check**:
- Network error on the SPC CSV: the downloader retries 4× with exponential backoff (1, 2, 4, 8 s). If all four fail, the script aborts with a clear `RuntimeError`.
- Upstream file changed: if the pinned SHA256 mismatches, the script aborts with a clear message instructing the user to delete `cache/data_sha256.txt` to re-pin.
- Python version < 3.8: `argparse` and f-string usage will fail.

---

## Step 4: Verify

```bash
cd /tmp/claw4s_auto_tornado-alley-is-shifting-east-signal-or-sampling-artifact && python3 analyze_tornado_shift.py --verify
```

**Expected output:** at least 28 check-lines indented with two spaces and starting with `PASS  ` (covering dataset size, per-stratum bootstrap CIs, permutation p-values, detection-fraction finiteness and clamping, sensitivity-analysis structure, SHA256 pin match, random-seed recording, physical-plausibility bounds, and bootstrap-replicate-count floors), followed by a blank line and `VERIFY OK: all checks passed`. Exit code 0.

**Failure modes**:
- Any `FAIL` line indicates a broken invariant; the script exits with code 1 and prints a list of failed checks.

---

## Limitations and Assumptions

The analysis and its conclusions rest on the following assumptions; when any of these is violated, the interpretation must be revisited.

1. **Instrument assumption (physics-based, not statistical)**: EF2+ tornadoes destroy houses and snap large trees, so damage surveys recover them regardless of radar density or population. This is strongest for 1997+; in the pre-NEXRAD era (1980–1996), damage surveys in sparsely populated parts of the Great Plains may have under-counted weak EF2s, which would inflate the observed EF2+ eastward slope and therefore make F_det an *under*-estimate of the true detection-attributable fraction (§6e isolates the post-NEXRAD window to guard against this).
2. **No covariate adjustment**: we deliberately do not regress out population-density gradients, because population and storm-climatology gradients are collinear and any regression adjustment is controversial (Ashley 2021, *Bull. Amer. Meteor. Soc.*). We instead rely on the physics-based instrumental variable (EF2+ = detection-robust).
3. **Unweighted starting-point centroid**: each tornado contributes exactly one start-point to the annual mean. We do *not* integrate the tornado path, do *not* weight by EF, and do *not* use event-hour or storm-mode information. §6c (path-length weighting) is the robustness check; if §6c contradicted §3 we would flag it.
4. **OLS on annual centroids**: we fit a *linear* trend to annual centroids. If the underlying centroid trajectory is highly nonlinear (e.g., a step change at 1997), the OLS slope summarises it crudely. The pre/post-NEXRAD split (§6a) and the pre/post-EF-scale split (§6f) address this partially but not exhaustively.
5. **What the results do NOT show**: this analysis does *not* show (a) where individual tornadoes occurred, (b) whether tornado frequency has changed (we analyse centroid, not count), (c) whether climate change is the specific cause of any residual shift, (d) whether the shift extends beyond 2022, (e) whether the result generalises to other types of severe weather.
6. **Single-vendor catalogue**: SPC is the only public CONUS-wide tornado catalogue at this coverage. A cross-check against an independent catalogue (NCEI Storm Events Database, or a reanalysis-based dynamical tornado climatology) is outside scope; if SPC's data-collection methodology changed systematically in a way that preferentially affected the Southeast, that would be indistinguishable from a climatological signal in this analysis (though it would also bias EF2+ in the same way, which is why F_det remains the cleanest comparison).

## Success Criteria (whole skill)

1. **Script completion**: Step 3 prints `[9/9] ANALYSIS COMPLETE in <t> s` and exits 0.
2. **All verification assertions pass**: Step 4 prints `VERIFY OK: all checks passed` and exits 0. All ≥ 18 individual `PASS` lines appear before that.
3. **Output artefacts exist**: `results.json`, `report.md`, `cache/1950-2024_actual_tornadoes.csv`, and `cache/data_sha256.txt` are all present in the workspace.
4. **Results JSON structural completeness**: `results.json` contains populated top-level keys: `meta`, `counts`, `slopes_deg_per_year`, `bootstrap`, `bootstrap_two_stage`, `detection_fraction`, `permutation`, `sensitivity`, `annual_counts`, `annual_centroids_lon`.
5. **Data integrity**: `meta.data_sha256` equals `cd6dc7b14d0efffe19bd39266b09e4d0b65da70ae132316825c511ecdb776cf2` (the pinned SPC 2025-05-13 release).
6. **CIs and p-values are finite**: every `bootstrap.<stratum>.ci_lo_km_decade`, `ci_hi_km_decade`, `detection_fraction.ci_lo`, `ci_hi`, and `permutation.<stratum>.two_sided_p` is a finite number, and every CI is ordered (`lo ≤ point ≤ hi`).
7. **Effect sizes within plausible range**: `|bootstrap.all.point_slope_km_decade| ≤ 200 km/decade` (plausible physical range); `|bootstrap.strong.point_slope_km_decade| ≤ 100 km/decade`.
8. **Dataset size**: `counts.total ≥ 30,000` and `counts.strong ≥ 3,000` (statistical power floor).
9. **Runtime**: full analysis completes in < 10 minutes on a typical laptop.

## Failure Conditions

- **Hard failures** (the skill has not executed correctly):
  - Any step returns a non-zero exit code.
  - `VERIFY FAILED` line in step-4 output.
  - `results.json` missing any of the top-level keys listed in Success Criterion 4.
  - Runtime > 10 minutes (indicates an infinite loop or quadratic blow-up in the bootstrap).
  - `SHA256 mismatch` error from `fetch_tornado_csv()` (upstream CSV has changed since this skill was pinned; adapt the `DATA_SHA256_EXPECTED` constant and re-pin).
- **Soft / scientific failures** (the skill executed but its conclusions should be treated as suspect):
  - `bootstrap.all.ci_lo_km_decade` and `ci_hi_km_decade` both contain zero (no detectable all-tornado shift — the whole premise of the question is weakened).
  - `counts.strong < 3000` after a data update (statistical power to resolve an EF2+ slope of ≈ 5 km/decade vs. 0 is lost).
  - `|bootstrap.all.point_slope_km_decade| > 200` (physically implausible; check for a data-parsing regression).