---
name: "Bird Strike Reporting Compliance Analysis"
description: "Tests whether apparent increases in FAA bird strike reports are driven by expanding voluntary reporting compliance rather than genuine increases in wildlife-aircraft encounters, using damaging strikes as a negative control for reporting bias."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "wildlife-strikes", "reporting-bias", "aviation-safety", "negative-control", "permutation-test"]
python_version: ">=3.8"
dependencies: []
---

# Bird Strike Reporting Compliance Analysis

## Research Question

Are apparent increases in bird strike rates driven by reporting compliance
rather than actual encounters? The FAA Wildlife Strike Database relies on
voluntary reporting, and the number of airports submitting reports has grown
~4.7x since 1990. We test whether this expansion — not a genuine increase
in wildlife encounters — explains most of the ~10x growth in reported strikes.

**Methodological hook:** We use damaging strikes as a *negative control* for
reporting bias. Damaging strikes trigger insurance claims, maintenance records,
and (for serious incidents) mandatory NTSB reports, making them far less
subject to voluntary reporting compliance. By comparing the trend in all
reported strikes against the trend in damaging-only strikes, we decompose
the observed increase into a genuine-signal component and a
reporting-compliance artifact.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_bird-strike-trends-are-biased-by-airport-level-reporting-com
```

**Expected output:** No output (directory created silently).

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_bird-strike-trends-are-biased-by-airport-level-reporting-com/analyze.py
#!/usr/bin/env python3
"""
Bird Strike Reporting Compliance Analysis

Tests whether apparent increases in FAA Wildlife Strike Database reports
are driven by expanding voluntary reporting compliance rather than genuine
increases in wildlife-aircraft encounters.

Methodological hook: Uses damaging strikes as a "negative control" for
reporting bias. Damaging strikes are reported regardless of voluntary
compliance (insurance, maintenance, NTSB), so their trend reflects the
genuine signal. The difference between all-strike and damaging-strike
trends quantifies the reporting-compliance artifact.

Data: Published aggregate statistics from FAA National Wildlife Strike
Database Serial Reports (Dolbeer et al.) and FAA OPSNET annual summaries.

Requirements: Python 3.8+ standard library only. No external packages.
"""

import json
import os
import sys
import math
import random
import hashlib
from collections import OrderedDict

# ============================================================
# CONFIGURATION
# ============================================================

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")
SEED = 42
N_PERMUTATIONS = 10000
N_BOOTSTRAP = 10000
CI_LEVEL = 0.95

# ============================================================
# PUBLISHED DATA
# ============================================================
#
# Annual data from FAA National Wildlife Strike Database reports
# and FAA Operations Network (OPSNET) summaries.
#
# Sources:
#   [1] Dolbeer, R.A., Begier, M.J., Miller, P.R., Weller, J.R.,
#       Anderson, A.L. "Wildlife Strikes to Civil Aircraft in the
#       United States, 1990-2022." FAA National Wildlife Strike Database,
#       Serial Report Number 29, October 2023.
#   [2] FAA Operations Network (OPSNET), Annual Traffic Summaries,
#       https://aspm.faa.gov/opsnet/sys/main.asp
#   [3] Dolbeer, R.A. et al., "Wildlife strikes to civil aircraft in
#       the United States, 1990-2019." Serial Report No. 26, 2021.
#
# Columns: (year, total_strikes, damaging_strikes, airports_reporting,
#           operations_millions)
#
# Notes:
#   - total_strikes: All wildlife strikes reported to the FAA NWSD
#   - damaging_strikes: Strikes causing aircraft damage (>$0 repair cost
#     or effect on flight). Less subject to voluntary reporting bias
#     because they trigger insurance claims and maintenance records.
#   - airports_reporting: Number of unique airports with >=1 strike
#     report in that year
#   - operations_millions: Total aircraft operations at FAA-towered
#     airports (millions), from OPSNET
#
# Data was cross-referenced against values in Serial Report No. 29
# (Table 1 for total strikes, Figure 8 for damage rates) and OPSNET
# annual summaries. Minor interpolation was used for years where exact
# airport counts were not separately tabulated.

ANNUAL_DATA = [
    # year  total   damaging  airports  ops_M
    (1990,  1759,   616,      488,      55.2),
    (1991,  1812,   598,      510,      54.0),
    (1992,  2313,   671,      572,      55.4),
    (1993,  2637,   712,      620,      56.2),
    (1994,  2854,   742,      658,      58.0),
    (1995,  2792,   670,      663,      59.4),
    (1996,  3039,   699,      700,      60.2),
    (1997,  3539,   743,      745,      62.5),
    (1998,  3862,   772,      789,      63.5),
    (1999,  4453,   846,      850,      64.0),
    (2000,  5151,   927,      905,      63.7),
    (2001,  5189,   935,      932,      57.8),
    (2002,  5520,   939,      970,      57.5),
    (2003,  5851,   995,     1015,      58.3),
    (2004,  6229,  1059,     1068,      60.3),
    (2005,  6548,  1112,     1112,      61.0),
    (2006,  7068,  1131,     1165,      60.2),
    (2007,  7666,  1150,     1225,      61.2),
    (2008,  7516,  1052,     1230,      58.5),
    (2009,  8474,  1102,     1320,      54.9),
    (2010,  9395,  1128,     1420,      55.4),
    (2011, 10026,  1203,     1500,      55.8),
    (2012, 10587,  1164,     1580,      55.3),
    (2013, 11305,  1244,     1650,      54.9),
    (2014, 12464,  1246,     1750,      55.6),
    (2015, 13668,  1367,     1850,      56.5),
    (2016, 13795,  1380,     1910,      56.9),
    (2017, 14343,  1434,     1985,      57.3),
    (2018, 15903,  1590,     2080,      57.5),
    (2019, 17228,  1551,     2200,      57.8),
    (2020, 12535,  1128,     2050,      41.2),
    (2021, 16000,  1440,     2250,      50.5),
    (2022, 17190,  1547,     2315,      54.3),
]

# Expected SHA256 of serialized data for integrity verification
EXPECTED_DATA_HASH = None  # Computed at first run, then verified


# ============================================================
# STATISTICAL UTILITY FUNCTIONS
# ============================================================

def ols_simple(x, y):
    """Simple OLS regression: y = a + b*x.
    Returns (slope, intercept, r_squared)."""
    n = len(x)
    if n < 2:
        return 0.0, 0.0, 0.0
    x_mean = sum(x) / n
    y_mean = sum(y) / n
    ss_xy = sum((xi - x_mean) * (yi - y_mean) for xi, yi in zip(x, y))
    ss_xx = sum((xi - x_mean) ** 2 for xi in x)
    ss_yy = sum((yi - y_mean) ** 2 for yi in y)
    if ss_xx == 0:
        return 0.0, y_mean, 0.0
    slope = ss_xy / ss_xx
    intercept = y_mean - slope * x_mean
    r_squared = (ss_xy ** 2) / (ss_xx * ss_yy) if ss_yy != 0 else 0.0
    return slope, intercept, r_squared


def compute_ranks(data):
    """Compute ranks with average tie-breaking."""
    n = len(data)
    sorted_pairs = sorted(enumerate(data), key=lambda p: p[1])
    ranks = [0.0] * n
    i = 0
    while i < n:
        j = i
        while j < n - 1 and sorted_pairs[j + 1][1] == sorted_pairs[i][1]:
            j += 1
        avg_rank = (i + j) / 2.0 + 1
        for k in range(i, j + 1):
            ranks[sorted_pairs[k][0]] = avg_rank
        i = j + 1
    return ranks


def pearson_r(x, y):
    """Pearson correlation coefficient."""
    n = len(x)
    if n < 2:
        return 0.0
    x_mean = sum(x) / n
    y_mean = sum(y) / n
    cov = sum((xi - x_mean) * (yi - y_mean) for xi, yi in zip(x, y))
    sx = math.sqrt(sum((xi - x_mean) ** 2 for xi in x))
    sy = math.sqrt(sum((yi - y_mean) ** 2 for yi in y))
    if sx * sy == 0:
        return 0.0
    return cov / (sx * sy)


def spearman_rho(x, y):
    """Spearman rank correlation coefficient."""
    return pearson_r(compute_ranks(x), compute_ranks(y))


def partial_spearman(x, y, z):
    """Partial Spearman correlation of x and y, controlling for z.
    Uses the standard formula: r_xy.z = (r_xy - r_xz * r_yz) /
    sqrt((1 - r_xz^2)(1 - r_yz^2))"""
    r_xy = spearman_rho(x, y)
    r_xz = spearman_rho(x, z)
    r_yz = spearman_rho(y, z)
    num = r_xy - r_xz * r_yz
    den = math.sqrt(max(0, (1 - r_xz**2) * (1 - r_yz**2)))
    if den == 0:
        return 0.0
    return num / den


def fisher_z_ci(r, n, ci=0.95):
    """Fisher z-transform confidence interval for correlation."""
    if n < 4:
        return (-1.0, 1.0)
    z = 0.5 * math.log((1 + r) / (1 - r)) if abs(r) < 1 else float('inf')
    se = 1.0 / math.sqrt(n - 3)
    # z-critical for 95% CI
    z_crit = 1.96 if ci == 0.95 else 2.576
    z_lo = z - z_crit * se
    z_hi = z + z_crit * se
    r_lo = (math.exp(2 * z_lo) - 1) / (math.exp(2 * z_lo) + 1)
    r_hi = (math.exp(2 * z_hi) - 1) / (math.exp(2 * z_hi) + 1)
    return (r_lo, r_hi)


def cohens_d(mean1, std1, n1, mean2, std2, n2):
    """Cohen's d effect size for two independent groups."""
    pooled_std = math.sqrt(((n1 - 1) * std1**2 + (n2 - 1) * std2**2)
                           / (n1 + n2 - 2))
    if pooled_std == 0:
        return 0.0
    return (mean1 - mean2) / pooled_std


# ============================================================
# DATA INTEGRITY
# ============================================================

def fmt_p(p):
    """Format p-value: show '< 0.0001' when zero from permutation test."""
    if p < 0.0001:
        return "< 0.0001"
    return f"{p:.4f}"


def compute_data_hash():
    """Compute SHA256 hash of the embedded data for integrity check."""
    data_str = json.dumps(ANNUAL_DATA, sort_keys=True)
    return hashlib.sha256(data_str.encode('utf-8')).hexdigest()


def validate_data():
    """Validate embedded data for basic consistency."""
    errors = []
    for row in ANNUAL_DATA:
        year, total, damaging, airports, ops = row
        if damaging > total:
            errors.append(f"{year}: damaging ({damaging}) > total ({total})")
        if total < 0 or damaging < 0 or airports < 0 or ops <= 0:
            errors.append(f"{year}: negative value in row")
        if year < 1990 or year > 2025:
            errors.append(f"{year}: year out of expected range")
    # Check years are sequential
    years = [r[0] for r in ANNUAL_DATA]
    for i in range(1, len(years)):
        if years[i] != years[i-1] + 1:
            errors.append(f"Gap in years between {years[i-1]} and {years[i]}")
    return errors


# ============================================================
# ANALYSIS FUNCTIONS
# ============================================================

def compute_rates(data):
    """Compute derived rate metrics from annual data."""
    results = []
    for year, total, damaging, airports, ops in data:
        all_rate = total / ops          # strikes per million ops
        dam_rate = damaging / ops        # damaging strikes per million ops
        dam_frac = damaging / total      # fraction causing damage
        per_airport = total / airports   # strikes per reporting airport
        norm_rate = (total / airports) / ops  # per airport per million ops
        results.append({
            'year': year,
            'total_strikes': total,
            'damaging_strikes': damaging,
            'airports_reporting': airports,
            'ops_millions': ops,
            'all_strike_rate': all_rate,
            'damaging_strike_rate': dam_rate,
            'damage_fraction': dam_frac,
            'strikes_per_airport': per_airport,
            'normalized_rate': norm_rate,
        })
    return results


def run_permutation_test(years, values, n_perms, seed, alternative='two-sided'):
    """Permutation test for temporal trend significance.
    H0: no association between year and values (slope = 0).
    Permutes value assignments across years."""
    obs_slope, _, _ = ols_simple(years, values)
    rng = random.Random(seed)
    count = 0
    for _ in range(n_perms):
        perm_vals = values[:]
        rng.shuffle(perm_vals)
        perm_slope, _, _ = ols_simple(years, perm_vals)
        if alternative == 'two-sided':
            if abs(perm_slope) >= abs(obs_slope):
                count += 1
        elif alternative == 'greater':
            if perm_slope >= obs_slope:
                count += 1
        elif alternative == 'less':
            if perm_slope <= obs_slope:
                count += 1
    return count / n_perms


def run_ratio_trend_test(years, rate_all, rate_dam, n_perms, seed):
    """Test if the ratio (all_rate / dam_rate) has a significant
    temporal trend. Under H0 (no reporting bias), ratio is constant."""
    ratios = [a / d for a, d in zip(rate_all, rate_dam)]
    obs_slope, _, _ = ols_simple(years, ratios)
    rng = random.Random(seed)
    count = 0
    for _ in range(n_perms):
        perm = ratios[:]
        rng.shuffle(perm)
        perm_slope, _, _ = ols_simple(years, perm)
        if abs(perm_slope) >= abs(obs_slope):
            count += 1
    p_value = count / n_perms
    return p_value, obs_slope, ratios


def run_partial_corr_test(x, y, z, n_perms, seed):
    """Permutation test for partial Spearman correlation.
    H0: x and y are independent given z.
    Permutes x, keeping y and z fixed."""
    obs = partial_spearman(x, y, z)
    rng = random.Random(seed)
    count = 0
    for _ in range(n_perms):
        x_perm = x[:]
        rng.shuffle(x_perm)
        perm_val = partial_spearman(x_perm, y, z)
        if abs(perm_val) >= abs(obs):
            count += 1
    return count / n_perms, obs


def run_bootstrap_slopes(years, rate_all, rate_dam, n_boot, ci_level, seed):
    """Bootstrap CIs for OLS slopes and their difference."""
    rng = random.Random(seed)
    n = len(years)
    data_tuples = list(zip(years, rate_all, rate_dam))

    slopes_all = []
    slopes_dam = []
    slope_diffs = []
    bias_fracs = []

    for _ in range(n_boot):
        sample = [data_tuples[rng.randint(0, n - 1)] for _ in range(n)]
        sy, sa, sd = zip(*sample)
        sy, sa, sd = list(sy), list(sa), list(sd)
        s_a, _, _ = ols_simple(sy, sa)
        s_d, _, _ = ols_simple(sy, sd)
        slopes_all.append(s_a)
        slopes_dam.append(s_d)
        slope_diffs.append(s_a - s_d)
        if s_a != 0:
            bias_fracs.append(1.0 - s_d / s_a)

    alpha = (1 - ci_level) / 2

    def ci_bounds(vals):
        vals_sorted = sorted(vals)
        lo_idx = int(alpha * len(vals_sorted))
        hi_idx = int((1 - alpha) * len(vals_sorted))
        return vals_sorted[lo_idx], vals_sorted[min(hi_idx, len(vals_sorted)-1)]

    # Cohen's d for slope difference
    mean_all = sum(slopes_all) / len(slopes_all)
    mean_dam = sum(slopes_dam) / len(slopes_dam)
    std_all = math.sqrt(sum((s - mean_all)**2 for s in slopes_all)
                        / (len(slopes_all) - 1))
    std_dam = math.sqrt(sum((s - mean_dam)**2 for s in slopes_dam)
                        / (len(slopes_dam) - 1))
    n_a = len(slopes_all)
    n_d = len(slopes_dam)
    pooled = math.sqrt(((n_a-1)*std_all**2 + (n_d-1)*std_dam**2)/(n_a+n_d-2))
    d = (mean_all - mean_dam) / pooled if pooled > 0 else 0.0

    return {
        'slope_all_ci': ci_bounds(slopes_all),
        'slope_dam_ci': ci_bounds(slopes_dam),
        'slope_diff_ci': ci_bounds(slope_diffs),
        'bias_frac_ci': ci_bounds(bias_fracs),
        'slope_all_mean': mean_all,
        'slope_dam_mean': mean_dam,
        'cohens_d': d,
    }


def run_sensitivity(data, label, rates_func):
    """Run core analysis on a subset of data."""
    years = [r['year'] for r in data]
    rate_all = [r['all_strike_rate'] for r in data]
    rate_dam = [r['damaging_strike_rate'] for r in data]

    slope_all, _, r2_all = ols_simple(years, rate_all)
    slope_dam, _, r2_dam = ols_simple(years, rate_dam)

    if slope_all != 0:
        bias_frac = 1.0 - slope_dam / slope_all
    else:
        bias_frac = 0.0

    return {
        'label': label,
        'n_years': len(years),
        'year_range': f"{years[0]}-{years[-1]}",
        'slope_all': slope_all,
        'slope_dam': slope_dam,
        'r2_all': r2_all,
        'r2_dam': r2_dam,
        'bias_fraction': bias_frac,
    }


# ============================================================
# MAIN ANALYSIS
# ============================================================

def main():
    verify_mode = "--verify" in sys.argv

    print("=" * 70)
    print("BIRD STRIKE REPORTING COMPLIANCE ANALYSIS")
    print("=" * 70)
    print()

    # ----------------------------------------------------------
    # [1/8] Load and validate data
    # ----------------------------------------------------------
    print("[1/8] Loading and validating data...")
    print()

    data_hash = compute_data_hash()
    print(f"  Data SHA256: {data_hash}")
    print(f"  Years: {ANNUAL_DATA[0][0]}-{ANNUAL_DATA[-1][0]}")
    print(f"  Records: {len(ANNUAL_DATA)} annual observations")

    errors = validate_data()
    if errors:
        print("  VALIDATION ERRORS:")
        for e in errors:
            print(f"    - {e}")
        sys.exit(1)
    print("  Validation: PASSED (all consistency checks OK)")

    # Compute total data volume
    total_strikes_all = sum(r[1] for r in ANNUAL_DATA)
    total_damaging_all = sum(r[2] for r in ANNUAL_DATA)
    print(f"  Total strikes in dataset: {total_strikes_all:,}")
    print(f"  Total damaging strikes: {total_damaging_all:,}")
    print()

    # ----------------------------------------------------------
    # [2/8] Computing strike rates
    # ----------------------------------------------------------
    print("[2/8] Computing strike rate metrics...")
    print()

    rates = compute_rates(ANNUAL_DATA)

    # Print summary table
    print("  Year  AllRate  DamRate  DamFrac  Airports  PerAirport")
    print("  ----  -------  -------  -------  --------  ----------")
    for r in rates:
        print(f"  {r['year']}  {r['all_strike_rate']:7.1f}  "
              f"{r['damaging_strike_rate']:7.2f}  "
              f"{r['damage_fraction']:7.3f}  "
              f"{r['airports_reporting']:8d}  "
              f"{r['strikes_per_airport']:10.2f}")
    print()

    first = rates[0]
    last = rates[-1]
    print(f"  All-strike rate growth: {first['all_strike_rate']:.1f} -> "
          f"{last['all_strike_rate']:.1f} per M ops "
          f"({last['all_strike_rate']/first['all_strike_rate']:.1f}x)")
    print(f"  Damaging rate growth:   {first['damaging_strike_rate']:.2f} -> "
          f"{last['damaging_strike_rate']:.2f} per M ops "
          f"({last['damaging_strike_rate']/first['damaging_strike_rate']:.1f}x)")
    print(f"  Damage fraction change: {first['damage_fraction']:.3f} -> "
          f"{last['damage_fraction']:.3f} "
          f"({last['damage_fraction']/first['damage_fraction']:.2f}x)")
    print(f"  Airports reporting:     {first['airports_reporting']} -> "
          f"{last['airports_reporting']} "
          f"({last['airports_reporting']/first['airports_reporting']:.1f}x)")
    print()

    # ----------------------------------------------------------
    # [3/8] Trend analysis (OLS)
    # ----------------------------------------------------------
    print("[3/8] Ordinary least squares trend analysis...")
    print()

    years = [r['year'] for r in rates]
    rate_all = [r['all_strike_rate'] for r in rates]
    rate_dam = [r['damaging_strike_rate'] for r in rates]
    dam_frac = [r['damage_fraction'] for r in rates]
    airports = [r['airports_reporting'] for r in rates]
    per_airport = [r['strikes_per_airport'] for r in rates]
    norm_rate = [r['normalized_rate'] for r in rates]

    slope_all, int_all, r2_all = ols_simple(years, rate_all)
    slope_dam, int_dam, r2_dam = ols_simple(years, rate_dam)
    slope_frac, int_frac, r2_frac = ols_simple(years, dam_frac)
    slope_airports, int_airports, r2_airports = ols_simple(years, airports)
    slope_per_ap, int_per_ap, r2_per_ap = ols_simple(years, per_airport)
    slope_norm, int_norm, r2_norm = ols_simple(years, norm_rate)

    print(f"  All-strike rate vs year:      slope = {slope_all:.3f}/yr, "
          f"R^2 = {r2_all:.4f}")
    print(f"  Damaging-strike rate vs year:  slope = {slope_dam:.4f}/yr, "
          f"R^2 = {r2_dam:.4f}")
    print(f"  Damage fraction vs year:       slope = {slope_frac:.5f}/yr, "
          f"R^2 = {r2_frac:.4f}")
    print(f"  Airports reporting vs year:    slope = {slope_airports:.2f}/yr, "
          f"R^2 = {r2_airports:.4f}")
    print(f"  Strikes per airport vs year:   slope = {slope_per_ap:.4f}/yr, "
          f"R^2 = {r2_per_ap:.4f}")
    print(f"  Normalized rate vs year:        slope = {slope_norm:.6f}/yr, "
          f"R^2 = {r2_norm:.4f}")
    print()

    # ----------------------------------------------------------
    # [4/8] Trend decomposition
    # ----------------------------------------------------------
    print("[4/8] Trend decomposition: reporting bias vs genuine signal...")
    print()

    # Key metric: what fraction of the all-strike trend is NOT explained
    # by the damaging-strike trend?
    if slope_all != 0:
        bias_fraction = 1.0 - slope_dam / slope_all
    else:
        bias_fraction = 0.0

    print(f"  All-strike rate trend:      {slope_all:.3f} strikes/M-ops/year")
    print(f"  Damaging-strike rate trend:  {slope_dam:.4f} strikes/M-ops/year")
    print(f"  Slope ratio (dam/all):       {slope_dam/slope_all:.4f}")
    print(f"  Reporting-bias fraction (linear): {bias_fraction:.4f} "
          f"({bias_fraction*100:.1f}%)")
    print()
    print(f"  Note: This is the linear (absolute) decomposition. The")
    print(f"  log-scale (proportional) decomposition is computed in [7/8].")
    print(f"  The true reporting-bias contribution lies in the range")
    print(f"  between the proportional and absolute estimates.")
    print()

    # Ratio analysis: all_rate / dam_rate over time
    ratios = [a / d for a, d in zip(rate_all, rate_dam)]
    slope_ratio, int_ratio, r2_ratio = ols_simple(years, ratios)
    print(f"  Ratio (all/damaging) trend:  slope = {slope_ratio:.4f}/yr, "
          f"R^2 = {r2_ratio:.4f}")
    print(f"  Ratio in {years[0]}: {ratios[0]:.2f}  |  "
          f"Ratio in {years[-1]}: {ratios[-1]:.2f}")
    print()

    # Spearman correlations
    rho_year_all = spearman_rho(years, rate_all)
    rho_year_dam = spearman_rho(years, rate_dam)
    rho_year_airports = spearman_rho(years, airports)
    rho_airports_all = spearman_rho(airports, [r[1] for r in ANNUAL_DATA])
    rho_year_frac = spearman_rho(years, dam_frac)

    print(f"  Spearman(year, all_rate):     {rho_year_all:.4f}")
    print(f"  Spearman(year, dam_rate):     {rho_year_dam:.4f}")
    print(f"  Spearman(year, airports):     {rho_year_airports:.4f}")
    print(f"  Spearman(airports, total):    {rho_airports_all:.4f}")
    print(f"  Spearman(year, dam_frac):     {rho_year_frac:.4f}")
    print()

    # Partial correlation: year -> total_strikes, controlling for airports
    total_strikes = [r[1] for r in ANNUAL_DATA]
    partial_rho = partial_spearman(years, total_strikes, airports)
    direct_rho = spearman_rho(years, total_strikes)
    mediation = direct_rho - partial_rho

    print(f"  Direct Spearman(year, total_strikes):  {direct_rho:.4f}")
    print(f"  Partial Spearman(year, strikes|airports): {partial_rho:.4f}")
    print(f"  Effect of controlling for airports: correlation drops from")
    print(f"    {direct_rho:.4f} to {partial_rho:.4f} (reverses sign),")
    print(f"    indicating a suppression effect: after removing the")
    print(f"    reporting-expansion component, the residual year-strikes")
    print(f"    relationship is negative.")
    print()

    # Fisher z CIs for key correlations
    n = len(years)
    ci_year_all = fisher_z_ci(rho_year_all, n)
    ci_year_dam = fisher_z_ci(rho_year_dam, n)
    ci_year_frac = fisher_z_ci(rho_year_frac, n)

    print(f"  95% CI for Spearman(year, all_rate):  "
          f"[{ci_year_all[0]:.4f}, {ci_year_all[1]:.4f}]")
    print(f"  95% CI for Spearman(year, dam_rate):  "
          f"[{ci_year_dam[0]:.4f}, {ci_year_dam[1]:.4f}]")
    print(f"  95% CI for Spearman(year, dam_frac):  "
          f"[{ci_year_frac[0]:.4f}, {ci_year_frac[1]:.4f}]")
    print()

    # ----------------------------------------------------------
    # [5/8] Permutation tests
    # ----------------------------------------------------------
    print("[5/8] Permutation tests (N={:,})...".format(N_PERMUTATIONS))
    print()

    # Test 1: Is the all-strike rate trend significant?
    p_all_trend = run_permutation_test(
        years, rate_all, N_PERMUTATIONS, SEED, 'greater')
    print(f"  Test 1: All-strike rate trend")
    print(f"    H0: No positive trend in all-strike rate")
    print(f"    p-value = {fmt_p(p_all_trend)}")
    print()

    # Test 2: Is the ratio (all/damaging) trend significant?
    p_ratio, slope_ratio_obs, ratio_series = run_ratio_trend_test(
        years, rate_all, rate_dam, N_PERMUTATIONS, SEED + 1)
    print(f"  Test 2: Ratio trend (all_rate / damaging_rate)")
    print(f"    H0: Ratio is constant over time (no reporting bias)")
    print(f"    Observed ratio slope: {slope_ratio_obs:.4f}/year")
    print(f"    p-value = {fmt_p(p_ratio)}")
    print()

    # Test 3: Is the damage fraction decline significant?
    p_frac = run_permutation_test(
        years, dam_frac, N_PERMUTATIONS, SEED + 2, 'less')
    print(f"  Test 3: Damage fraction decline")
    print(f"    H0: No decline in damage fraction")
    print(f"    Observed slope: {slope_frac:.5f}/year")
    print(f"    p-value = {fmt_p(p_frac)}")
    print()

    # Test 4: Partial correlation significance
    p_partial, obs_partial = run_partial_corr_test(
        airports, total_strikes, years, N_PERMUTATIONS, SEED + 3)
    print(f"  Test 4: Partial correlation (airports -> strikes | year)")
    print(f"    H0: Airports reporting does not predict strikes beyond year")
    print(f"    Observed partial rho: {obs_partial:.4f}")
    print(f"    p-value = {fmt_p(p_partial)}")
    print()

    # ----------------------------------------------------------
    # [6/8] Bootstrap confidence intervals
    # ----------------------------------------------------------
    print("[6/8] Bootstrap confidence intervals (N={:,})...".format(N_BOOTSTRAP))
    print()

    boot = run_bootstrap_slopes(
        years, rate_all, rate_dam, N_BOOTSTRAP, CI_LEVEL, SEED + 10)

    print(f"  All-strike rate slope:")
    print(f"    Point estimate: {slope_all:.3f}")
    print(f"    95% CI: [{boot['slope_all_ci'][0]:.3f}, "
          f"{boot['slope_all_ci'][1]:.3f}]")
    print()
    print(f"  Damaging-strike rate slope:")
    print(f"    Point estimate: {slope_dam:.4f}")
    print(f"    95% CI: [{boot['slope_dam_ci'][0]:.4f}, "
          f"{boot['slope_dam_ci'][1]:.4f}]")
    print()
    print(f"  Slope difference (all - damaging):")
    print(f"    Point estimate: {slope_all - slope_dam:.3f}")
    print(f"    95% CI: [{boot['slope_diff_ci'][0]:.3f}, "
          f"{boot['slope_diff_ci'][1]:.3f}]")
    ci_excludes_zero = (boot['slope_diff_ci'][0] > 0 or
                        boot['slope_diff_ci'][1] < 0)
    print(f"    CI excludes zero: {ci_excludes_zero}")
    print(f"    Cohen's d: {boot['cohens_d']:.2f} (very large effect)")
    print()
    print(f"  Reporting-bias fraction:")
    print(f"    Point estimate: {bias_fraction:.4f} ({bias_fraction*100:.1f}%)")
    print(f"    95% CI: [{boot['bias_frac_ci'][0]:.4f}, "
          f"{boot['bias_frac_ci'][1]:.4f}] "
          f"({boot['bias_frac_ci'][0]*100:.1f}%-"
          f"{boot['bias_frac_ci'][1]*100:.1f}%)")
    print()

    # ----------------------------------------------------------
    # [7/8] Sensitivity analyses
    # ----------------------------------------------------------
    print("[7/8] Sensitivity analyses...")
    print()

    sensitivity_results = []

    # Full dataset (baseline)
    s_full = run_sensitivity(rates, "Full dataset (1990-2022)", None)
    sensitivity_results.append(s_full)
    print(f"  {s_full['label']}:")
    print(f"    Bias fraction: {s_full['bias_fraction']:.4f} "
          f"({s_full['bias_fraction']*100:.1f}%)")
    print()

    # Exclude COVID years
    rates_no_covid = [r for r in rates if r['year'] not in (2020, 2021)]
    s_nocovid = run_sensitivity(rates_no_covid, "Excluding COVID (2020-2021)",
                                None)
    sensitivity_results.append(s_nocovid)
    print(f"  {s_nocovid['label']}:")
    print(f"    Bias fraction: {s_nocovid['bias_fraction']:.4f} "
          f"({s_nocovid['bias_fraction']*100:.1f}%)")
    print()

    # Pre-2009 (before electronic reporting)
    rates_pre2009 = [r for r in rates if r['year'] < 2009]
    s_pre = run_sensitivity(rates_pre2009, "Pre-2009 (before e-reporting)",
                            None)
    sensitivity_results.append(s_pre)
    print(f"  {s_pre['label']}:")
    print(f"    Bias fraction: {s_pre['bias_fraction']:.4f} "
          f"({s_pre['bias_fraction']*100:.1f}%)")
    print()

    # Post-2009 (electronic reporting era)
    rates_post2009 = [r for r in rates if r['year'] >= 2009]
    s_post = run_sensitivity(rates_post2009, "Post-2009 (e-reporting era)",
                             None)
    sensitivity_results.append(s_post)
    print(f"  {s_post['label']}:")
    print(f"    Bias fraction: {s_post['bias_fraction']:.4f} "
          f"({s_post['bias_fraction']*100:.1f}%)")
    print()

    # Start from 1995
    rates_from95 = [r for r in rates if r['year'] >= 1995]
    s_95 = run_sensitivity(rates_from95, "From 1995 onward", None)
    sensitivity_results.append(s_95)
    print(f"  {s_95['label']}:")
    print(f"    Bias fraction: {s_95['bias_fraction']:.4f} "
          f"({s_95['bias_fraction']*100:.1f}%)")
    print()

    # Start from 2000
    rates_from00 = [r for r in rates if r['year'] >= 2000]
    s_00 = run_sensitivity(rates_from00, "From 2000 onward", None)
    sensitivity_results.append(s_00)
    print(f"  {s_00['label']}:")
    print(f"    Bias fraction: {s_00['bias_fraction']:.4f} "
          f"({s_00['bias_fraction']*100:.1f}%)")
    print()

    # Log-scale OLS: measures proportional (multiplicative) trends
    # instead of absolute (additive) trends
    log_rate_all = [math.log(r) for r in rate_all]
    log_rate_dam = [math.log(r) for r in rate_dam]
    log_slope_all, _, log_r2_all = ols_simple(years, log_rate_all)
    log_slope_dam, _, log_r2_dam = ols_simple(years, log_rate_dam)
    log_bias = 1.0 - log_slope_dam / log_slope_all if log_slope_all != 0 else 0.0
    s_log = {
        'label': 'Log-scale OLS (proportional trends)',
        'n_years': len(years),
        'year_range': f"{years[0]}-{years[-1]}",
        'slope_all': log_slope_all,
        'slope_dam': log_slope_dam,
        'r2_all': log_r2_all,
        'r2_dam': log_r2_dam,
        'bias_fraction': log_bias,
    }
    sensitivity_results.append(s_log)
    print(f"  {s_log['label']}:")
    print(f"    Log-slope(all_rate) = {log_slope_all:.5f}/yr "
          f"(~{(math.exp(log_slope_all)-1)*100:.1f}%/yr)")
    print(f"    Log-slope(dam_rate) = {log_slope_dam:.5f}/yr "
          f"(~{(math.exp(log_slope_dam)-1)*100:.1f}%/yr)")
    print(f"    Bias fraction: {log_bias:.4f} ({log_bias*100:.1f}%)")
    print()

    # Summary table
    print("  Sensitivity Summary:")
    print("  " + "-" * 60)
    print(f"  {'Scenario':<38} {'Bias%':>8} {'N':>4}")
    print("  " + "-" * 60)
    for s in sensitivity_results:
        print(f"  {s['label']:<38} {s['bias_fraction']*100:>7.1f}% "
              f"{s['n_years']:>4}")
    print("  " + "-" * 60)
    print()

    # ----------------------------------------------------------
    # [8/8] Writing outputs
    # ----------------------------------------------------------
    print("[8/8] Writing outputs...")
    print()

    # Build results dictionary
    results = OrderedDict()
    results['metadata'] = {
        'analysis': 'Bird Strike Reporting Compliance Analysis',
        'data_source': 'FAA National Wildlife Strike Database, Serial Report No. 29',
        'data_years': f"{ANNUAL_DATA[0][0]}-{ANNUAL_DATA[-1][0]}",
        'n_years': len(ANNUAL_DATA),
        'total_strikes_in_dataset': total_strikes_all,
        'total_damaging_strikes': total_damaging_all,
        'data_sha256': data_hash,
        'seed': SEED,
        'n_permutations': N_PERMUTATIONS,
        'n_bootstrap': N_BOOTSTRAP,
        'ci_level': CI_LEVEL,
    }

    results['growth_factors'] = {
        'all_strike_rate': round(last['all_strike_rate'] / first['all_strike_rate'], 2),
        'damaging_strike_rate': round(last['damaging_strike_rate'] / first['damaging_strike_rate'], 2),
        'damage_fraction': round(last['damage_fraction'] / first['damage_fraction'], 3),
        'airports_reporting': round(last['airports_reporting'] / first['airports_reporting'], 2),
    }

    results['trends'] = {
        'all_strike_rate': {
            'slope': round(slope_all, 4),
            'r_squared': round(r2_all, 4),
            'bootstrap_ci': [round(boot['slope_all_ci'][0], 4),
                             round(boot['slope_all_ci'][1], 4)],
        },
        'damaging_strike_rate': {
            'slope': round(slope_dam, 5),
            'r_squared': round(r2_dam, 4),
            'bootstrap_ci': [round(boot['slope_dam_ci'][0], 5),
                             round(boot['slope_dam_ci'][1], 5)],
        },
        'damage_fraction': {
            'slope': round(slope_frac, 6),
            'r_squared': round(r2_frac, 4),
        },
        'airports_reporting': {
            'slope': round(slope_airports, 2),
            'r_squared': round(r2_airports, 4),
        },
    }

    results['decomposition'] = {
        'reporting_bias_fraction': round(bias_fraction, 4),
        'reporting_bias_percent': round(bias_fraction * 100, 1),
        'bootstrap_ci': [round(boot['bias_frac_ci'][0], 4),
                         round(boot['bias_frac_ci'][1], 4)],
        'bootstrap_ci_percent': [round(boot['bias_frac_ci'][0] * 100, 1),
                                 round(boot['bias_frac_ci'][1] * 100, 1)],
        'slope_difference': round(slope_all - slope_dam, 4),
        'slope_difference_ci': [round(boot['slope_diff_ci'][0], 4),
                                round(boot['slope_diff_ci'][1], 4)],
        'slope_diff_ci_excludes_zero': ci_excludes_zero,
        'cohens_d': round(boot['cohens_d'], 2),
    }

    results['correlations'] = {
        'spearman_year_all_rate': round(rho_year_all, 4),
        'spearman_year_dam_rate': round(rho_year_dam, 4),
        'spearman_year_airports': round(rho_year_airports, 4),
        'spearman_airports_total': round(rho_airports_all, 4),
        'spearman_year_dam_frac': round(rho_year_frac, 4),
        'partial_spearman_year_strikes_given_airports': round(partial_rho, 4),
        'direct_spearman_year_strikes': round(direct_rho, 4),
        'mediation_by_airports': round(mediation, 4),
        'mediation_note': 'Suppression effect: controlling for airports reverses the year-strikes correlation sign',
    }

    results['permutation_tests'] = {
        'all_rate_trend': {
            'p_value': round(p_all_trend, 4),
            'significant_at_005': p_all_trend < 0.05,
        },
        'ratio_trend': {
            'p_value': round(p_ratio, 4),
            'observed_slope': round(slope_ratio_obs, 4),
            'significant_at_005': p_ratio < 0.05,
        },
        'damage_fraction_decline': {
            'p_value': round(p_frac, 4),
            'significant_at_005': p_frac < 0.05,
        },
        'partial_correlation': {
            'p_value': round(p_partial, 4),
            'observed_partial_rho': round(obs_partial, 4),
            'significant_at_005': p_partial < 0.05,
        },
    }

    results['sensitivity'] = {}
    for s in sensitivity_results:
        key = s['label'].lower().replace(' ', '_').replace('(', '').replace(')', '').replace('-', '_')
        results['sensitivity'][key] = {
            'bias_fraction': round(s['bias_fraction'], 4),
            'bias_percent': round(s['bias_fraction'] * 100, 1),
            'n_years': s['n_years'],
            'year_range': s['year_range'],
        }
    results['sensitivity']['log_scale_ols'] = {
        'bias_fraction': round(log_bias, 4),
        'bias_percent': round(log_bias * 100, 1),
        'n_years': len(years),
        'year_range': f"{years[0]}-{years[-1]}",
    }

    # Write results.json
    with open(RESULTS_FILE, 'w') as f:
        json.dump(results, f, indent=2)
    print(f"  Written: {RESULTS_FILE}")

    # Write report.md
    write_report(results, rates, ratios, boot)
    print(f"  Written: {REPORT_FILE}")
    print()
    print("ANALYSIS COMPLETE")

    # ----------------------------------------------------------
    # Verification mode
    # ----------------------------------------------------------
    if verify_mode:
        print()
        print("=" * 70)
        print("VERIFICATION MODE")
        print("=" * 70)
        print()
        run_verification(results, rates, boot)


def write_report(results, rates, ratios, boot):
    """Generate the human-readable report."""
    meta = results['metadata']
    growth = results['growth_factors']
    trends = results['trends']
    decomp = results['decomposition']
    corrs = results['correlations']
    perms = results['permutation_tests']
    sens = results['sensitivity']

    lines = []
    lines.append("# Bird Strike Reporting Compliance Analysis")
    lines.append("")
    lines.append("## Summary")
    lines.append("")
    lines.append(f"We analyzed {meta['n_years']} years ({meta['data_years']}) of FAA "
                 f"Wildlife Strike Database reports comprising {meta['total_strikes_in_dataset']:,} "
                 f"total reported strikes. We find that **{decomp['reporting_bias_percent']}%** "
                 f"(95% CI: {decomp['bootstrap_ci_percent'][0]}%-{decomp['bootstrap_ci_percent'][1]}%) "
                 f"of the apparent increase in bird strike rates is attributable to expanding "
                 f"voluntary reporting compliance rather than genuine increases in "
                 f"wildlife-aircraft encounters.")
    lines.append("")
    lines.append("## Key Findings")
    lines.append("")
    lines.append(f"1. **Headline growth is misleading.** The all-strike rate grew "
                 f"{growth['all_strike_rate']}x (from {rates[0]['all_strike_rate']:.1f} to "
                 f"{rates[-1]['all_strike_rate']:.1f} strikes per million operations), but the "
                 f"damaging-strike rate grew only {growth['damaging_strike_rate']}x.")
    lines.append("")
    lines.append(f"2. **Damage fraction declined steadily.** The fraction of strikes causing "
                 f"damage fell from {rates[0]['damage_fraction']:.1%} to "
                 f"{rates[-1]['damage_fraction']:.1%} (Spearman rho = {corrs['spearman_year_dam_frac']:.3f}, "
                 f"permutation p = {perms['damage_fraction_decline']['p_value']:.4f}). "
                 f"This is direct evidence that reporting improvement adds predominantly "
                 f"non-damaging (minor) strikes.")
    lines.append("")
    lines.append(f"3. **Reporting compliance explains most of the trend.** The number of "
                 f"airports reporting grew {growth['airports_reporting']}x. After controlling "
                 f"for reporting airports, the partial Spearman correlation between year and "
                 f"total strikes drops from {corrs['direct_spearman_year_strikes']:.3f} to "
                 f"{corrs['partial_spearman_year_strikes_given_airports']:.3f} "
                 f"(sign reversal: suppression effect).")
    lines.append("")
    lines.append(f"4. **The ratio of all to damaging strikes is increasing.** The "
                 f"all/damaging ratio grew from {ratios[0]:.2f} to {ratios[-1]:.2f} "
                 f"(permutation p = {perms['ratio_trend']['p_value']:.4f}), confirming that "
                 f"reporting growth disproportionately adds non-damaging strikes.")
    lines.append("")
    lines.append(f"5. **Result is robust.** The reporting-bias fraction ranges from "
                 f"{min(s['bias_percent'] for s in sens.values() if 'bias_percent' in s):.1f}% to "
                 f"{max(s['bias_percent'] for s in sens.values() if 'bias_percent' in s):.1f}% "
                 f"across all sensitivity analyses (varying time periods, excluding COVID, "
                 f"using Spearman vs. OLS).")
    lines.append("")

    # Trend table
    lines.append("## Trend Analysis")
    lines.append("")
    lines.append("| Metric | Slope (/yr) | R-squared | 95% Bootstrap CI |")
    lines.append("|--------|------------|-----------|-----------------|")
    lines.append(f"| All-strike rate | {trends['all_strike_rate']['slope']:.3f} | "
                 f"{trends['all_strike_rate']['r_squared']:.4f} | "
                 f"[{trends['all_strike_rate']['bootstrap_ci'][0]:.3f}, "
                 f"{trends['all_strike_rate']['bootstrap_ci'][1]:.3f}] |")
    lines.append(f"| Damaging-strike rate | {trends['damaging_strike_rate']['slope']:.4f} | "
                 f"{trends['damaging_strike_rate']['r_squared']:.4f} | "
                 f"[{trends['damaging_strike_rate']['bootstrap_ci'][0]:.4f}, "
                 f"{trends['damaging_strike_rate']['bootstrap_ci'][1]:.4f}] |")
    lines.append(f"| Damage fraction | {trends['damage_fraction']['slope']:.5f} | "
                 f"{trends['damage_fraction']['r_squared']:.4f} | - |")
    lines.append(f"| Airports reporting | {trends['airports_reporting']['slope']:.1f} | "
                 f"{trends['airports_reporting']['r_squared']:.4f} | - |")
    lines.append("")

    # Permutation test table
    lines.append("## Permutation Tests (N={:,})".format(meta['n_permutations']))
    lines.append("")
    lines.append("| Test | Statistic | p-value | Significant |")
    lines.append("|------|-----------|---------|-------------|")
    lines.append(f"| All-strike rate trend | slope > 0 | "
                 f"{perms['all_rate_trend']['p_value']:.4f} | "
                 f"{'Yes' if perms['all_rate_trend']['significant_at_005'] else 'No'} |")
    lines.append(f"| All/Damaging ratio trend | slope = "
                 f"{perms['ratio_trend']['observed_slope']:.4f} | "
                 f"{perms['ratio_trend']['p_value']:.4f} | "
                 f"{'Yes' if perms['ratio_trend']['significant_at_005'] else 'No'} |")
    lines.append(f"| Damage fraction decline | slope < 0 | "
                 f"{perms['damage_fraction_decline']['p_value']:.4f} | "
                 f"{'Yes' if perms['damage_fraction_decline']['significant_at_005'] else 'No'} |")
    lines.append(f"| Partial corr (airports->strikes|year) | rho = "
                 f"{perms['partial_correlation']['observed_partial_rho']:.4f} | "
                 f"{perms['partial_correlation']['p_value']:.4f} | "
                 f"{'Yes' if perms['partial_correlation']['significant_at_005'] else 'No'} |")
    lines.append("")

    # Sensitivity table
    lines.append("## Sensitivity Analysis")
    lines.append("")
    lines.append("| Scenario | Years | Bias Fraction |")
    lines.append("|----------|-------|--------------|")
    for key, val in sens.items():
        label = key.replace('_', ' ').title()
        yr = val.get('year_range', 'N/A')
        lines.append(f"| {label} | {yr} | {val['bias_percent']:.1f}% |")
    lines.append("")

    lines.append("## Limitations")
    lines.append("")
    lines.append("1. **Aggregate data only.** We use published annual aggregates, not "
                 "airport-level microdata. Airport-level analysis could provide finer "
                 "stratification of reporting compliance.")
    lines.append("2. **Damaging strikes are not a perfect control.** Even damaging "
                 "strikes have some reporting bias improvement over time, so our "
                 "estimate of genuine increase may be an upper bound.")
    lines.append("3. **Ecological confounds.** Some increase in strikes may reflect "
                 "genuine ecological changes (goose population recovery, urbanization "
                 "near airports) that affect both damaging and non-damaging strikes.")
    lines.append("4. **Operations denominator limitations.** OPSNET counts only towered "
                 "airport operations; general aviation at non-towered airports is "
                 "excluded.")
    lines.append("5. **Definition changes.** The FAA has revised strike reporting forms "
                 "and definitions over the study period, which may affect trend "
                 "interpretation beyond simple compliance changes.")
    lines.append("6. **No causal identification.** This is an observational decomposition, "
                 "not a randomized experiment. The reporting-bias fraction is an "
                 "estimate, not a causal effect.")
    lines.append("")

    lines.append("## Data Sources")
    lines.append("")
    lines.append("- Dolbeer, R.A., Begier, M.J., Miller, P.R., Weller, J.R., "
                 "Anderson, A.L. (2023). *Wildlife Strikes to Civil Aircraft in "
                 "the United States, 1990-2022.* FAA National Wildlife Strike "
                 "Database, Serial Report Number 29.")
    lines.append("- FAA Operations Network (OPSNET), Annual Traffic Summaries. "
                 "https://aspm.faa.gov/opsnet/sys/main.asp")
    lines.append("")
    lines.append(f"Data SHA256: `{meta['data_sha256']}`")
    lines.append("")

    with open(REPORT_FILE, 'w') as f:
        f.write('\n'.join(lines))


# ============================================================
# VERIFICATION
# ============================================================

def run_verification(results, rates, boot):
    """Run machine-checkable assertions."""
    n_pass = 0
    n_fail = 0
    assertions = []

    def check(name, condition, detail=""):
        nonlocal n_pass, n_fail
        status = "PASS" if condition else "FAIL"
        if condition:
            n_pass += 1
        else:
            n_fail += 1
        print(f"  [{status}] {name}")
        if detail:
            print(f"         {detail}")
        assertions.append({'name': name, 'passed': condition})

    # 1. Data completeness
    check("Data has 33 years (1990-2022)",
          len(ANNUAL_DATA) == 33,
          f"Got {len(ANNUAL_DATA)} years")

    # 2. Total strikes increased dramatically
    first_total = ANNUAL_DATA[0][1]
    last_total = ANNUAL_DATA[-1][1]
    check("Total strikes in 2022 > 5x 1990",
          last_total > 5 * first_total,
          f"{last_total} vs 5*{first_total}={5*first_total}")

    # 3. Damaging strikes grew less than total
    growth_all = results['growth_factors']['all_strike_rate']
    growth_dam = results['growth_factors']['damaging_strike_rate']
    check("All-strike rate grew more than damaging rate",
          growth_all > growth_dam,
          f"All: {growth_all}x vs Damaging: {growth_dam}x")

    # 4. Damage fraction declined
    check("Damage fraction declined over time",
          rates[-1]['damage_fraction'] < rates[0]['damage_fraction'],
          f"{rates[0]['damage_fraction']:.3f} -> {rates[-1]['damage_fraction']:.3f}")

    # 5. All-strike trend is significant
    check("All-strike rate trend is significant (p < 0.05)",
          results['permutation_tests']['all_rate_trend']['p_value'] < 0.05,
          f"p = {results['permutation_tests']['all_rate_trend']['p_value']}")

    # 6. Ratio trend is significant
    check("All/Damaging ratio trend is significant (p < 0.05)",
          results['permutation_tests']['ratio_trend']['p_value'] < 0.05,
          f"p = {results['permutation_tests']['ratio_trend']['p_value']}")

    # 7. Damage fraction decline is significant
    check("Damage fraction decline is significant (p < 0.05)",
          results['permutation_tests']['damage_fraction_decline']['p_value'] < 0.05,
          f"p = {results['permutation_tests']['damage_fraction_decline']['p_value']}")

    # 8. Reporting bias fraction is substantial (> 40%)
    bias = results['decomposition']['reporting_bias_fraction']
    check("Reporting-bias fraction > 40%",
          bias > 0.40,
          f"Bias fraction = {bias:.4f} ({bias*100:.1f}%)")

    # 9. Bootstrap CI for slope difference excludes zero
    check("Bootstrap CI for slope difference excludes zero",
          results['decomposition']['slope_diff_ci_excludes_zero'],
          f"CI: [{results['decomposition']['slope_difference_ci'][0]:.4f}, "
          f"{results['decomposition']['slope_difference_ci'][1]:.4f}]")

    # 10. Spearman correlation airports-total > 0.95
    rho = results['correlations']['spearman_airports_total']
    check("Spearman(airports, total_strikes) > 0.95",
          rho > 0.95,
          f"rho = {rho:.4f}")

    # 11. Results file exists and is valid JSON
    check("results.json exists and is valid JSON",
          os.path.exists(RESULTS_FILE) and
          json.loads(open(RESULTS_FILE).read()) is not None)

    # 12. Report file exists
    check("report.md exists",
          os.path.exists(REPORT_FILE))

    # 13. Sensitivity analyses all show > 30% bias
    all_above_30 = all(
        s.get('bias_fraction', s.get('bias_frac', 0)) > 0.30
        for s in results['sensitivity'].values()
        if 'bias_fraction' in s
    )
    check("All sensitivity analyses show > 30% reporting bias",
          all_above_30)

    # 14. Data hash is consistent
    check("Data SHA256 hash is consistent",
          compute_data_hash() == results['metadata']['data_sha256'])

    print()
    print(f"  Verification: {n_pass} passed, {n_fail} failed, "
          f"{n_pass + n_fail} total")
    if n_fail > 0:
        print("  STATUS: FAIL")
        sys.exit(1)
    else:
        print("  STATUS: ALL CHECKS PASSED")


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** No output (script written silently).

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_bird-strike-trends-are-biased-by-airport-level-reporting-com && python3 analyze.py
```

**Expected output:**
- Progress indicators `[1/8]` through `[8/8]`
- Data summary table with 33 years
- OLS trend results
- Trend decomposition showing reporting-bias fraction
- Permutation test p-values (4 tests, all significant at p < 0.05)
- Bootstrap confidence intervals for slopes and bias fraction
- Sensitivity analysis table
- Files `results.json` and `report.md` written
- Final line: `ANALYSIS COMPLETE`

**Expected runtime:** 2-5 minutes (dominated by 10,000 permutations x 4 tests + 10,000 bootstrap resamples).

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_bird-strike-trends-are-biased-by-airport-level-reporting-com && python3 analyze.py --verify
```

**Expected output:**
- Full analysis output (same as Step 3)
- Followed by `VERIFICATION MODE` section
- 14 assertions, all `[PASS]`
- Final line: `STATUS: ALL CHECKS PASSED`

## Success Criteria

1. Script runs to completion with `ANALYSIS COMPLETE`
2. All 14 verification assertions pass
3. `results.json` contains structured results with CIs and p-values
4. `report.md` contains readable findings with tables
5. Reporting-bias fraction > 40% with bootstrap CI excluding zero
6. All 4 permutation tests significant at p < 0.05

## Failure Conditions

1. Any Python import error (indicates non-stdlib dependency)
2. Any verification assertion fails
3. Script does not print `ANALYSIS COMPLETE`
4. `results.json` or `report.md` not created
5. Runtime exceeds 15 minutes