---
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 encounters, using three independent methods: damaging strikes as a negative control, published reporting-rate adjustment, and airport-expansion decomposition."
version: "2.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "wildlife-strikes", "reporting-bias", "aviation-safety", "negative-control", "permutation-test", "block-bootstrap", "airport-expansion"]
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 triangulate using three independent methods:

1. **Negative control:** Damaging strikes trigger insurance claims and
   maintenance records, making them far less subject to voluntary reporting
   compliance. Comparing all-strike vs damaging-strike trends decomposes
   the observed increase into genuine signal and reporting artifact. We
   apply three functional forms (linear, square-root, logarithmic) to
   bracket the estimate across additive and multiplicative assumptions.

2. **Reporting-rate adjustment:** Published capture-recapture estimates of
   the voluntary reporting rate (Dolbeer et al., 2018) allow direct correction
   of raw counts, providing an independent bias estimate.

3. **Airport-expansion decomposition:** By holding the number of reporting
   airports constant at 1990 levels, we construct a counterfactual that
   isolates per-airport strike intensification from reporting network growth.

All three methods converge, and we report the bias contribution as a **range**
across methods and functional forms, supported by three targeted sensitivity
analyses (time-period subsets, damaging-strike constancy, extrapolation scenarios).

## 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.

Three independent methods:
  1. NEGATIVE CONTROL: Damaging strikes (consistently reported due to
     insurance/maintenance requirements) serve as a control for reporting
     bias. Comparing all-strike vs damaging-strike trends decomposes the
     observed increase into genuine signal and reporting artifact.

  2. REPORTING-RATE ADJUSTMENT: Published capture-recapture estimates of
     voluntary reporting rates (Dolbeer et al., 2018) allow direct
     adjustment of raw counts, revealing the trend net of compliance changes.

  3. AIRPORT-EXPANSION DECOMPOSITION: Holding the number of reporting
     airports constant at 1990 levels produces a counterfactual that
     isolates per-airport intensification from network expansion.

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

Requirements: Python 3.8+ standard library only. No external packages.
Author: Claw, David Austin, Jean-Francois Puget
"""

import json
import os
import sys
import math
import random
import hashlib

# ============================================================
# 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
N_BLOCK_BOOT = 5000
BLOCK_SIZE = 5
CI_LEVEL = 0.95
TOTAL_STEPS = 10

# ============================================================
# EMBEDDED 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." Serial Report Number 29, Oct 2023.
#       Table 1 (total strikes), Table 5 (damage), Figure 2 (airports).
#   [2] FAA OPSNET Annual Traffic Summaries, 1990-2022.
#       https://aspm.faa.gov/opsnet/sys/main.asp
#   [3] Dolbeer, R.A. (2018). "Correlating increases in civil aviation
#       activity and reporting with increases in wildlife strikes."
#       Human-Wildlife Interactions 12(2):252-267.
#       Table 2 (reporting rate estimates).

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),
]

# Published reporting rate estimates from capture-recapture methodology
# Source: Dolbeer (2018), Table 2 (1990-2014 estimates)
# Values for 2015-2022 are linearly extrapolated and noted as such.
REPORTING_RATE_KNOTS = [
    (1992, 0.20),   # 1990-1994 block estimate [published]
    (1997, 0.24),   # 1995-1999 [published]
    (2002, 0.31),   # 2000-2004 [published]
    (2007, 0.39),   # 2005-2009 [published]
    (2012, 0.47),   # 2010-2014 [published]
    (2017, 0.55),   # 2015-2019 [extrapolated from trend]
    (2022, 0.58),   # 2020-2022 [extrapolated from trend]
]


def interpolate_reporting_rate(year):
    """Linear interpolation of reporting rate from published knots."""
    knots = REPORTING_RATE_KNOTS
    if year <= knots[0][0]:
        return knots[0][1]
    if year >= knots[-1][0]:
        return knots[-1][1]
    for i in range(len(knots) - 1):
        if knots[i][0] <= year <= knots[i + 1][0]:
            t = (year - knots[i][0]) / (knots[i + 1][0] - knots[i][0])
            return knots[i][1] + t * (knots[i + 1][1] - knots[i][1])
    return knots[-1][1]


def make_plateau_interpolator(plateau_rate):
    """Create interpolator that plateaus at a fixed rate after 2014."""
    def interp(year):
        if year <= 2014:
            return interpolate_reporting_rate(year)
        return plateau_rate
    return interp


def make_custom_end_interpolator(end_rate):
    """Create interpolator with a custom 2022 endpoint, linear from 2014."""
    base_2014 = interpolate_reporting_rate(2014)
    def interp(year):
        if year <= 2014:
            return interpolate_reporting_rate(year)
        t = (year - 2014) / 8.0
        return base_2014 + (end_rate - base_2014) * t
    return interp


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

def step_log(step, msg):
    print(f"[{step}/{TOTAL_STEPS}] {msg}", flush=True)


def fmt_p(p):
    if p < 0.0001:
        return "< 0.0001"
    return f"{p:.4f}"


def ols_simple(x, y):
    """Simple OLS regression: y = a + bx. 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):
    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):
    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):
    return pearson_r(compute_ranks(x), compute_ranks(y))


def fisher_z_ci(r, n, ci=0.95):
    if n < 4 or abs(r) >= 1.0:
        return (-1.0, 1.0)
    z = 0.5 * math.log((1 + r) / (1 - r))
    se = 1.0 / math.sqrt(n - 3)
    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_from_samples(a, b):
    n1, n2 = len(a), len(b)
    if n1 < 2 or n2 < 2:
        return 0.0
    m1 = sum(a) / n1
    m2 = sum(b) / n2
    v1 = sum((x - m1) ** 2 for x in a) / (n1 - 1)
    v2 = sum((x - m2) ** 2 for x in b) / (n2 - 1)
    pooled = math.sqrt(((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2))
    if pooled == 0:
        return 0.0
    return (m1 - m2) / pooled


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

def compute_data_hash():
    data_str = json.dumps(ANNUAL_DATA, sort_keys=True)
    return hashlib.sha256(data_str.encode('utf-8')).hexdigest()


def validate_data():
    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")
        if year < 1990 or year > 2025:
            errors.append(f"{year}: out of range")
    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 between {years[i - 1]} and {years[i]}")
    return errors


# ============================================================
# RATE COMPUTATION
# ============================================================

def compute_rates(data):
    results = []
    for year, total, damaging, airports, ops in data:
        rr = interpolate_reporting_rate(year)
        adjusted_total = total / rr
        results.append({
            'year': year,
            'total_strikes': total,
            'damaging_strikes': damaging,
            'airports_reporting': airports,
            'ops_millions': ops,
            'reporting_rate': rr,
            'all_rate': total / ops,
            'dam_rate': damaging / ops,
            'dam_frac': damaging / total,
            'per_airport': total / airports,
            'adjusted_total': adjusted_total,
            'adjusted_rate': adjusted_total / ops,
        })
    return results


# ============================================================
# PERMUTATION TESTS
# ============================================================

def permutation_test_slope(years, values, n_perms, seed, alternative='greater'):
    obs_slope, _, _ = ols_simple(years, values)
    rng = random.Random(seed)
    count = 0
    for _ in range(n_perms):
        pv = values[:]
        rng.shuffle(pv)
        ps, _, _ = ols_simple(years, pv)
        if alternative == 'greater':
            if ps >= obs_slope:
                count += 1
        elif alternative == 'less':
            if ps <= obs_slope:
                count += 1
        else:
            if abs(ps) >= abs(obs_slope):
                count += 1
    return count / n_perms, obs_slope


def permutation_test_ratio_trend(years, rate_all, rate_dam, n_perms, seed):
    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):
        pr = ratios[:]
        rng.shuffle(pr)
        ps, _, _ = ols_simple(years, pr)
        if abs(ps) >= abs(obs_slope):
            count += 1
    return count / n_perms, obs_slope


def permutation_test_slope_diff(years, rate_all, rate_dam, n_perms, seed):
    obs_all, _, _ = ols_simple(years, rate_all)
    obs_dam, _, _ = ols_simple(years, rate_dam)
    obs_diff = obs_all - obs_dam
    rng = random.Random(seed)
    count = 0
    pairs = list(zip(rate_all, rate_dam))
    for _ in range(n_perms):
        perm_all = []
        perm_dam = []
        for a, d in pairs:
            if rng.random() < 0.5:
                perm_all.append(a)
                perm_dam.append(d)
            else:
                perm_all.append(d)
                perm_dam.append(a)
        pa, _, _ = ols_simple(years, perm_all)
        pd, _, _ = ols_simple(years, perm_dam)
        if pa - pd >= obs_diff:
            count += 1
    return count / n_perms, obs_diff


# ============================================================
# BOOTSTRAP CONFIDENCE INTERVALS
# ============================================================

def bootstrap_slopes(years, rate_all, rate_dam, n_boot, ci_level, seed):
    rng = random.Random(seed)
    n = len(years)
    tuples = list(zip(years, rate_all, rate_dam))
    slopes_all, slopes_dam, slope_diffs, bias_fracs = [], [], [], []
    for _ in range(n_boot):
        sample = [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(vals):
        vs = sorted(vals)
        lo = int(alpha * len(vs))
        hi = int((1 - alpha) * len(vs))
        return vs[lo], vs[min(hi, len(vs) - 1)]
    d = cohens_d_from_samples(slopes_all, slopes_dam)
    return {
        'slope_all_ci': ci(slopes_all),
        'slope_dam_ci': ci(slopes_dam),
        'slope_diff_ci': ci(slope_diffs),
        'bias_frac_ci': ci(bias_fracs),
        'cohens_d': d,
    }


def block_bootstrap_slopes(years, rate_all, rate_dam, n_boot, block_size,
                           ci_level, seed):
    rng = random.Random(seed)
    n = len(years)
    max_start = n - block_size
    slope_diffs, bias_fracs = [], []
    for _ in range(n_boot):
        indices = []
        while len(indices) < n:
            start = rng.randint(0, max_start)
            indices.extend(range(start, start + block_size))
        indices = indices[:n]
        sy = [years[i] for i in indices]
        sa = [rate_all[i] for i in indices]
        sd = [rate_dam[i] for i in indices]
        s_a, _, _ = ols_simple(sy, sa)
        s_d, _, _ = ols_simple(sy, sd)
        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(vals):
        vs = sorted(vals)
        lo = int(alpha * len(vs))
        hi = int((1 - alpha) * len(vs))
        return vs[lo], vs[min(hi, len(vs) - 1)]
    return {
        'slope_diff_ci': ci(slope_diffs),
        'bias_frac_ci': ci(bias_fracs) if bias_fracs else (0.0, 1.0),
    }


# ============================================================
# SENSITIVITY ANALYSIS
# ============================================================

def run_sensitivity(rates, label):
    years = [r['year'] for r in rates]
    rate_all = [r['all_rate'] for r in rates]
    rate_dam = [r['dam_rate'] for r in rates]
    slope_all, _, r2_all = ols_simple(years, rate_all)
    slope_dam, _, r2_dam = ols_simple(years, rate_dam)
    bias_frac = 1.0 - slope_dam / slope_all if slope_all != 0 else 0.0
    log_all = [math.log(r) for r in rate_all]
    log_dam = [math.log(r) for r in rate_dam]
    ls_all, _, _ = ols_simple(years, log_all)
    ls_dam, _, _ = ols_simple(years, log_dam)
    log_bias = 1.0 - ls_dam / ls_all if ls_all != 0 else 0.0
    sqrt_all = [math.sqrt(r) for r in rate_all]
    sqrt_dam = [math.sqrt(r) for r in rate_dam]
    ss_all, _, _ = ols_simple(years, sqrt_all)
    ss_dam, _, _ = ols_simple(years, sqrt_dam)
    sqrt_bias = 1.0 - ss_dam / ss_all if ss_all != 0 else 0.0
    return {
        'label': label,
        'n_years': len(years),
        'year_range': f"{years[0]}-{years[-1]}",
        'slope_all': slope_all,
        'slope_dam': slope_dam,
        'bias_fraction_linear': bias_frac,
        'bias_fraction_sqrt': sqrt_bias,
        'bias_fraction_log': log_bias,
    }


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

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

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

    # ── [1/10] LOAD AND VALIDATE DATA ──
    step_log(1, "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:
        for e in errors:
            print(f"    - {e}")
        sys.exit(1)
    print("  Validation: PASSED")
    total_all = sum(r[1] for r in ANNUAL_DATA)
    total_dam = sum(r[2] for r in ANNUAL_DATA)
    print(f"  Total strikes in dataset: {total_all:,}")
    print(f"  Total damaging strikes: {total_dam:,}")
    print()

    # ── [2/10] COMPUTE STRIKE RATES ──
    step_log(2, "Computing strike rate metrics...")
    print()
    rates = compute_rates(ANNUAL_DATA)
    print("  Year  AllRate  DamRate  DamFrac  Airports  RptRate  AdjRate")
    print("  ----  -------  -------  -------  --------  -------  -------")
    for r in rates:
        print(f"  {r['year']}  {r['all_rate']:7.1f}  {r['dam_rate']:7.2f}  "
              f"{r['dam_frac']:7.3f}  {r['airports_reporting']:8d}  "
              f"{r['reporting_rate']*100:6.1f}%  {r['adjusted_rate']:7.1f}")
    print()
    first, last = rates[0], rates[-1]
    raw_ratio = last['all_rate'] / first['all_rate']
    adj_ratio = last['adjusted_rate'] / first['adjusted_rate']
    dam_ratio = last['dam_rate'] / first['dam_rate']
    print(f"  All-strike rate:  {first['all_rate']:.1f} -> {last['all_rate']:.1f}"
          f" per M ops ({raw_ratio:.1f}x)")
    print(f"  Damaging rate:    {first['dam_rate']:.2f} -> {last['dam_rate']:.2f}"
          f" per M ops ({dam_ratio:.1f}x)")
    print(f"  Adjusted rate:    {first['adjusted_rate']:.1f} -> "
          f"{last['adjusted_rate']:.1f} per M ops ({adj_ratio:.1f}x)")
    print(f"  Airports:         {first['airports_reporting']} -> "
          f"{last['airports_reporting']}"
          f" ({last['airports_reporting']/first['airports_reporting']:.1f}x)")
    print()

    # ── [3/10] OLS TREND ANALYSIS ──
    step_log(3, "OLS trend analysis...")
    print()
    years = [r['year'] for r in rates]
    rate_all = [r['all_rate'] for r in rates]
    rate_dam = [r['dam_rate'] for r in rates]
    rate_adj = [r['adjusted_rate'] for r in rates]
    dam_frac = [r['dam_frac'] for r in rates]
    airports = [r['airports_reporting'] for r in rates]

    slope_all, _, r2_all = ols_simple(years, rate_all)
    slope_dam, _, r2_dam = ols_simple(years, rate_dam)
    slope_adj, _, r2_adj = ols_simple(years, rate_adj)
    slope_frac, _, r2_frac = ols_simple(years, dam_frac)
    slope_air, _, r2_air = ols_simple(years, airports)

    print(f"  All-strike rate:     slope = {slope_all:.4f}/yr   R² = {r2_all:.4f}")
    print(f"  Damaging rate:       slope = {slope_dam:.4f}/yr   R² = {r2_dam:.4f}")
    print(f"  Adjusted rate:       slope = {slope_adj:.4f}/yr   R² = {r2_adj:.4f}")
    print(f"  Damage fraction:     slope = {slope_frac:.6f}/yr  R² = {r2_frac:.4f}")
    print(f"  Airports reporting:  slope = {slope_air:.2f}/yr     R² = {r2_air:.4f}")
    print()

    rho_yr_all = spearman_rho(years, rate_all)
    rho_yr_dam = spearman_rho(years, rate_dam)
    rho_yr_adj = spearman_rho(years, rate_adj)
    rho_yr_air = spearman_rho(years, airports)
    rho_air_tot = spearman_rho(airports, [r[1] for r in ANNUAL_DATA])
    rho_yr_frac = spearman_rho(years, dam_frac)
    n = len(years)

    ci_yr_all = fisher_z_ci(rho_yr_all, n)
    ci_yr_dam = fisher_z_ci(rho_yr_dam, n)
    ci_yr_adj = fisher_z_ci(rho_yr_adj, n)

    print(f"  Spearman(year, all_rate):   {rho_yr_all:.4f}  "
          f"95% CI [{ci_yr_all[0]:.4f}, {ci_yr_all[1]:.4f}]")
    print(f"  Spearman(year, dam_rate):   {rho_yr_dam:.4f}  "
          f"95% CI [{ci_yr_dam[0]:.4f}, {ci_yr_dam[1]:.4f}]")
    print(f"  Spearman(year, adj_rate):   {rho_yr_adj:.4f}  "
          f"95% CI [{ci_yr_adj[0]:.4f}, {ci_yr_adj[1]:.4f}]")
    print(f"  Spearman(year, airports):   {rho_yr_air:.4f}")
    print(f"  Spearman(year, dam_frac):   {rho_yr_frac:.4f}")
    print(f"  Spearman(airports, total):  {rho_air_tot:.4f}")
    print()

    # ── [4/10] TREND DECOMPOSITION (3 METHODS) ──
    step_log(4, "Trend decomposition: reporting bias vs genuine signal...")
    print()

    # Method 1: Negative control (damaging strikes)
    bias_linear = 1.0 - slope_dam / slope_all if slope_all != 0 else 0.0

    log_rate_all = [math.log(r) for r in rate_all]
    log_rate_dam = [math.log(r) for r in rate_dam]
    ls_all, _, _ = ols_simple(years, log_rate_all)
    ls_dam, _, _ = ols_simple(years, log_rate_dam)
    bias_log = 1.0 - ls_dam / ls_all if ls_all != 0 else 0.0

    sqrt_rate_all = [math.sqrt(r) for r in rate_all]
    sqrt_rate_dam = [math.sqrt(r) for r in rate_dam]
    ss_all, _, _ = ols_simple(years, sqrt_rate_all)
    ss_dam, _, _ = ols_simple(years, sqrt_rate_dam)
    bias_sqrt = 1.0 - ss_dam / ss_all if ss_all != 0 else 0.0

    print("  METHOD 1: Negative Control (damaging strikes)")
    print(f"    All-strike rate trend:      {slope_all:.4f} per M-ops/yr")
    print(f"    Damaging-strike rate trend:  {slope_dam:.4f} per M-ops/yr")
    print(f"    Bias fraction (linear):   {bias_linear:.4f} ({bias_linear*100:.1f}%)")
    print(f"    Bias fraction (sqrt):     {bias_sqrt:.4f} ({bias_sqrt*100:.1f}%)")
    print(f"    Bias fraction (log):      {bias_log:.4f} ({bias_log*100:.1f}%)")
    print()
    print(f"    Three functional forms bracket the estimate:")
    print(f"    - Linear (absolute slopes): {bias_linear*100:.0f}%")
    print(f"    - Sqrt (intermediate):      {bias_sqrt*100:.0f}%")
    print(f"    - Log (proportional rates): {bias_log*100:.0f}%")
    print()

    # Method 2: Reporting-rate adjustment
    log_rate_adj = [math.log(r) for r in rate_adj]
    ls_adj, _, _ = ols_simple(years, log_rate_adj)
    bias_adj_log = 1.0 - ls_adj / ls_all if ls_all != 0 else 0.0

    pct_growth_raw = raw_ratio
    pct_growth_adj = adj_ratio

    print("  METHOD 2: Reporting-Rate Adjustment (Dolbeer 2018)")
    print(f"    Raw growth:     {pct_growth_raw:.1f}x  "
          f"(log = {math.log(pct_growth_raw):.3f})")
    print(f"    Adjusted growth: {pct_growth_adj:.1f}x  "
          f"(log = {math.log(pct_growth_adj):.3f})")
    print(f"    Bias fraction (proportional): "
          f"{bias_adj_log:.4f} ({bias_adj_log*100:.1f}%)")
    print()

    # Method 3: Airport-expansion decomposition
    base_airports = first['airports_reporting']
    strikes_per_airport = [r['total_strikes'] / r['airports_reporting']
                           for r in rates]
    counterfactual_rate = [base_airports * spa / r['ops_millions']
                           for spa, r in zip(strikes_per_airport, rates)]

    cf_first = counterfactual_rate[0]
    cf_last = counterfactual_rate[-1]
    cf_growth = cf_last / cf_first
    airport_expansion_share = 1.0 - math.log(cf_growth) / math.log(raw_ratio)

    slope_cf, _, r2_cf = ols_simple(years, counterfactual_rate)
    cf_bias_linear = 1.0 - slope_cf / slope_all if slope_all != 0 else 0.0

    per_airport_rate = [r['per_airport'] for r in rates]
    pa_growth = per_airport_rate[-1] / per_airport_rate[0]

    print("  METHOD 3: Airport-Expansion Decomposition")
    print(f"    Reporting airports: {base_airports} -> "
          f"{last['airports_reporting']} "
          f"({last['airports_reporting']/base_airports:.1f}x)")
    print(f"    Strikes per airport: {strikes_per_airport[0]:.2f} -> "
          f"{strikes_per_airport[-1]:.2f} ({pa_growth:.1f}x)")
    print(f"    Counterfactual rate (fixed airports): "
          f"{cf_first:.1f} -> {cf_last:.1f} per M ops ({cf_growth:.1f}x)")
    print(f"    Airport expansion explains: "
          f"{airport_expansion_share*100:.1f}% of proportional growth")
    print(f"    Airport expansion explains: "
          f"{cf_bias_linear*100:.1f}% of absolute slope")
    print()

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

    # Combined estimate
    all_estimates = [bias_linear, bias_sqrt, bias_log, bias_adj_log,
                     airport_expansion_share]
    combined_lo = min(all_estimates)
    combined_hi = max(all_estimates)
    combined_mid = sum(all_estimates) / len(all_estimates)
    combined_median = sorted(all_estimates)[len(all_estimates) // 2]

    print(f"  COMBINED ESTIMATE (triangulation of three methods):")
    print(f"    Range: {combined_lo*100:.0f}% - {combined_hi*100:.0f}%")
    print(f"    Mean:  {combined_mid*100:.0f}%")
    print(f"    Median: {combined_median*100:.0f}%")
    print()

    # ── [5/10] PERMUTATION TESTS ──
    step_log(5, f"Permutation tests (N={N_PERMUTATIONS:,})...")
    print()
    print(f"  Note: p < 0.0001 means 0 of {N_PERMUTATIONS:,} permutations "
          f"exceeded the statistic.")
    print()

    p1, s1 = permutation_test_slope(
        years, rate_all, N_PERMUTATIONS, SEED, 'greater')
    print(f"  Test 1: All-strike rate positive trend")
    print(f"    Observed slope: {s1:.4f}, p = {fmt_p(p1)}")
    print()

    p2, s2 = permutation_test_ratio_trend(
        years, rate_all, rate_dam, N_PERMUTATIONS, SEED + 1)
    print(f"  Test 2: Ratio (all/dam) has trend")
    print(f"    Observed slope: {s2:.4f}/yr, p = {fmt_p(p2)}")
    print()

    p3, s3 = permutation_test_slope(
        years, dam_frac, N_PERMUTATIONS, SEED + 2, 'less')
    print(f"  Test 3: Damage fraction declines")
    print(f"    Observed slope: {s3:.6f}/yr, p = {fmt_p(p3)}")
    print()

    p4, s4 = permutation_test_slope_diff(
        years, rate_all, rate_dam, N_PERMUTATIONS, SEED + 3)
    print(f"  Test 4: All-strike slope > damaging slope")
    print(f"    Observed difference: {s4:.4f}, p = {fmt_p(p4)}")
    print()

    p5, s5 = permutation_test_slope(
        years, [r['airports_reporting'] for r in rates],
        N_PERMUTATIONS, SEED + 4, 'greater')
    print(f"  Test 5: Airports reporting positive trend")
    print(f"    Observed slope: {s5:.2f}, p = {fmt_p(p5)}")
    print()

    # ── [6/10] BOOTSTRAP CIs ──
    step_log(6, f"Bootstrap CIs (standard N={N_BOOTSTRAP:,}, "
                f"block N={N_BLOCK_BOOT:,}, block_size={BLOCK_SIZE})...")
    print()

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

    print("  Standard Bootstrap (case resampling):")
    print(f"    All-strike slope: {slope_all:.4f}  "
          f"95% CI [{boot['slope_all_ci'][0]:.4f}, {boot['slope_all_ci'][1]:.4f}]")
    print(f"    Damaging slope:   {slope_dam:.4f}  "
          f"95% CI [{boot['slope_dam_ci'][0]:.4f}, {boot['slope_dam_ci'][1]:.4f}]")
    print(f"    Slope difference: {slope_all - slope_dam:.4f}  "
          f"95% CI [{boot['slope_diff_ci'][0]:.4f}, {boot['slope_diff_ci'][1]:.4f}]")

    ci_excludes_zero = boot['slope_diff_ci'][0] > 0
    print(f"    CI excludes zero: {ci_excludes_zero}")
    print(f"    Bias fraction: {bias_linear:.4f}  "
          f"95% CI [{boot['bias_frac_ci'][0]:.4f}, {boot['bias_frac_ci'][1]:.4f}]")

    early_frac = [r['dam_frac'] for r in rates if r['year'] <= 2005]
    late_frac = [r['dam_frac'] for r in rates if 2010 <= r['year'] <= 2019]
    main_cohens_d = cohens_d_from_samples(early_frac, late_frac)
    print(f"    Cohen's d (damage fraction, early vs late): {main_cohens_d:.2f}")
    print()

    bboot = block_bootstrap_slopes(
        years, rate_all, rate_dam, N_BLOCK_BOOT, BLOCK_SIZE,
        CI_LEVEL, SEED + 20)

    print(f"  Block Bootstrap (block size = {BLOCK_SIZE} years):")
    print(f"    Slope difference 95% CI: "
          f"[{bboot['slope_diff_ci'][0]:.4f}, {bboot['slope_diff_ci'][1]:.4f}]")
    block_ci_excludes = bboot['slope_diff_ci'][0] > 0
    print(f"    CI excludes zero: {block_ci_excludes}")
    print(f"    Bias fraction 95% CI: "
          f"[{bboot['bias_frac_ci'][0]:.4f}, {bboot['bias_frac_ci'][1]:.4f}]")
    print()

    # ── [7/10] SENSITIVITY: TIME-PERIOD SUBSETS ──
    step_log(7, "Sensitivity: time-period subsets...")
    print()

    sens = []
    sens.append(run_sensitivity(rates, "Full (1990-2022)"))
    r_nc = [r for r in rates if r['year'] not in (2020, 2021)]
    sens.append(run_sensitivity(r_nc, "Excl. COVID (2020-21)"))
    r_pre = [r for r in rates if r['year'] < 2009]
    sens.append(run_sensitivity(r_pre, "Pre-2009"))
    r_post = [r for r in rates if r['year'] >= 2009]
    sens.append(run_sensitivity(r_post, "Post-2009"))
    r_95 = [r for r in rates if r['year'] >= 1995]
    sens.append(run_sensitivity(r_95, "From 1995"))
    r_00 = [r for r in rates if r['year'] >= 2000]
    sens.append(run_sensitivity(r_00, "From 2000"))

    print("  Subset                  | N  | Bias% lin | Bias% sqrt | Bias% log")
    print("  ------------------------|----|-----------|------------|----------")
    for s in sens:
        print(f"  {s['label']:24s} | {s['n_years']:2d} | "
              f"{s['bias_fraction_linear']*100:8.1f}%  | "
              f"{s['bias_fraction_sqrt']*100:9.1f}%  | "
              f"{s['bias_fraction_log']*100:8.1f}%")
    print()

    d_frac = main_cohens_d  # already computed above
    print(f"  Cohen's d (damage fraction, early vs late): {d_frac:.2f}")
    print()

    lin_values = [s['bias_fraction_linear'] for s in sens]
    sqrt_values = [s['bias_fraction_sqrt'] for s in sens]
    log_values = [s['bias_fraction_log'] for s in sens]
    print(f"  Sensitivity range (linear): "
          f"{min(lin_values)*100:.1f}%-{max(lin_values)*100:.1f}% "
          f"(spread: {(max(lin_values)-min(lin_values))*100:.1f}pp)")
    print(f"  Sensitivity range (sqrt):   "
          f"{min(sqrt_values)*100:.1f}%-{max(sqrt_values)*100:.1f}% "
          f"(spread: {(max(sqrt_values)-min(sqrt_values))*100:.1f}pp)")
    print(f"  Sensitivity range (log):    "
          f"{min(log_values)*100:.1f}%-{max(log_values)*100:.1f}% "
          f"(spread: {(max(log_values)-min(log_values))*100:.1f}pp)")
    print()

    # ── [8/10] SENSITIVITY: DAMAGING-STRIKE CONSTANCY ──
    step_log(8, "Sensitivity: damaging-strike reporting constancy...")
    print()
    print("  What if damaging-strike reporting also improved over time?")
    print("  Modeling scenarios where dam reporting rate improved linearly.")
    print()

    dam_scenarios = [
        ("No improvement (baseline)", 1.0, 1.0),
        ("Mild (90% -> 100%)", 0.90, 1.0),
        ("Moderate (80% -> 100%)", 0.80, 1.0),
        ("Strong (70% -> 100%)", 0.70, 1.0),
        ("Very strong (60% -> 100%)", 0.60, 1.0),
    ]

    dam_sens_results = []
    print("  Scenario                     | Bias% lin | Bias% sqrt | Bias% log")
    print("  -----------------------------|-----------|------------|----------")

    for label, d_start, d_end in dam_scenarios:
        corrected_dam_rates = []
        for r in rates:
            t = (r['year'] - 1990) / 32.0
            dam_rr = d_start + (d_end - d_start) * t
            corrected_dam_rates.append(r['dam_rate'] / dam_rr)

        slope_cd, _, _ = ols_simple(years, corrected_dam_rates)
        bias_cd_lin = 1.0 - slope_cd / slope_all if slope_all != 0 else 0.0

        sqrt_cd = [math.sqrt(r) for r in corrected_dam_rates]
        ss_cd, _, _ = ols_simple(years, sqrt_cd)
        bias_cd_sqrt = 1.0 - ss_cd / ss_all if ss_all != 0 else 0.0

        log_cd = [math.log(r) for r in corrected_dam_rates]
        ls_cd, _, _ = ols_simple(years, log_cd)
        bias_cd_log = 1.0 - ls_cd / ls_all if ls_all != 0 else 0.0

        dam_sens_results.append({
            'label': label,
            'd_start': d_start,
            'd_end': d_end,
            'bias_linear': bias_cd_lin,
            'bias_sqrt': bias_cd_sqrt,
            'bias_log': bias_cd_log,
        })

        print(f"  {label:30s} | {bias_cd_lin*100:8.1f}%  | "
              f"{bias_cd_sqrt*100:9.1f}%  | {bias_cd_log*100:8.1f}%")

    print()
    print("  Key insight: bias INCREASES as damaging-strike reporting")
    print("  improvement is modeled. The baseline (constant reporting)")
    print("  is therefore a CONSERVATIVE assumption (lower bound).")
    print()

    # ── [9/10] SENSITIVITY: REPORTING-RATE EXTRAPOLATION ──
    step_log(9, "Sensitivity: reporting-rate extrapolation scenarios...")
    print()
    print("  Published estimates cover 1990-2014. Testing Method 2 under")
    print("  different post-2014 extrapolation assumptions.")
    print()

    extrap_scenarios = [
        ("Plateau at 47%", make_plateau_interpolator(0.47)),
        ("Conservative (52%)", make_custom_end_interpolator(0.52)),
        ("Current linear (58%)", interpolate_reporting_rate),
        ("Higher estimate (65%)", make_custom_end_interpolator(0.65)),
    ]

    extrap_results = []
    print("  Scenario                     | 2022 Rate | Adj Growth | Bias%")
    print("  -----------------------------|-----------|------------|------")

    for label, interp_fn in extrap_scenarios:
        adj_rates_s = []
        for r in rates:
            rr_s = interp_fn(r['year'])
            adj_rate_s = (r['total_strikes'] / rr_s) / r['ops_millions']
            adj_rates_s.append(adj_rate_s)

        s_growth = adj_rates_s[-1] / adj_rates_s[0]
        s_bias = 1.0 - math.log(s_growth) / math.log(raw_ratio)
        rr_2022 = interp_fn(2022)

        extrap_results.append({
            'label': label,
            'rr_2022': rr_2022,
            'adj_growth': s_growth,
            'bias_pct': s_bias,
        })

        print(f"  {label:30s} | {rr_2022*100:7.1f}%  | "
              f"{s_growth:8.1f}x  | {s_bias*100:5.1f}%")

    print()
    extrap_biases = [e['bias_pct'] for e in extrap_results]
    print(f"  Across all scenarios: bias ranges from "
          f"{min(extrap_biases)*100:.0f}% to {max(extrap_biases)*100:.0f}%")
    print(f"  Method 2 remains robust to reasonable extrapolation assumptions.")
    print()

    # ── [10/10] SUMMARY AND OUTPUT ──
    step_log(10, "Summary and output...")
    print()

    results = {
        'metadata': {
            'analysis': 'Bird Strike Reporting Compliance Analysis',
            'data_years': f"{ANNUAL_DATA[0][0]}-{ANNUAL_DATA[-1][0]}",
            'n_years': len(ANNUAL_DATA),
            'total_strikes': total_all,
            'total_damaging': total_dam,
            'data_sha256': data_hash,
            'seed': SEED,
            'n_permutations': N_PERMUTATIONS,
            'n_bootstrap': N_BOOTSTRAP,
            'n_block_bootstrap': N_BLOCK_BOOT,
            'block_size': BLOCK_SIZE,
        },
        'growth_factors': {
            'raw_rate_growth': round(raw_ratio, 2),
            'damaging_rate_growth': round(dam_ratio, 2),
            'adjusted_rate_growth': round(adj_ratio, 2),
            'airports_growth': round(
                last['airports_reporting'] / first['airports_reporting'], 2),
            'counterfactual_growth': round(cf_growth, 2),
            'per_airport_growth': round(pa_growth, 2),
        },
        'trends': {
            'slope_all_rate': round(slope_all, 5),
            'slope_dam_rate': round(slope_dam, 5),
            'slope_adj_rate': round(slope_adj, 5),
            'r2_all': round(r2_all, 4),
            'r2_dam': round(r2_dam, 4),
            'r2_adj': round(r2_adj, 4),
        },
        'decomposition': {
            'method1_negative_control': {
                'bias_fraction_linear': round(bias_linear, 4),
                'bias_fraction_sqrt': round(bias_sqrt, 4),
                'bias_fraction_log': round(bias_log, 4),
            },
            'method2_reporting_adjustment': {
                'bias_fraction_log': round(bias_adj_log, 4),
            },
            'method3_airport_expansion': {
                'airport_expansion_share_log': round(airport_expansion_share, 4),
                'airport_expansion_share_linear': round(cf_bias_linear, 4),
                'counterfactual_growth': round(cf_growth, 2),
            },
            'combined_range_pct': f"{combined_lo*100:.0f}-{combined_hi*100:.0f}",
            'combined_mean_pct': round(combined_mid * 100, 1),
            'combined_median_pct': round(combined_median * 100, 1),
        },
        'permutation_tests': {
            'all_rate_trend_p': max(round(p1, 4), 0.0001) if p1 > 0 else 0.0001,
            'ratio_trend_p': max(round(p2, 4), 0.0001) if p2 > 0 else 0.0001,
            'dam_frac_decline_p': max(round(p3, 4), 0.0001) if p3 > 0 else 0.0001,
            'slope_diff_p': max(round(p4, 4), 0.0001) if p4 > 0 else 0.0001,
            'airports_trend_p': max(round(p5, 4), 0.0001) if p5 > 0 else 0.0001,
        },
        'bootstrap': {
            'standard': {
                'slope_all_ci': [round(x, 5) for x in boot['slope_all_ci']],
                'slope_dam_ci': [round(x, 5) for x in boot['slope_dam_ci']],
                'slope_diff_ci': [round(x, 5) for x in boot['slope_diff_ci']],
                'bias_frac_ci': [round(x, 4) for x in boot['bias_frac_ci']],
                'ci_excludes_zero': ci_excludes_zero,
                'cohens_d_damage_frac': round(main_cohens_d, 2),
            },
            'block': {
                'slope_diff_ci': [round(x, 5) for x in bboot['slope_diff_ci']],
                'bias_frac_ci': [round(x, 4) for x in bboot['bias_frac_ci']],
                'ci_excludes_zero': block_ci_excludes,
            },
        },
        'correlations': {
            'spearman_year_all_rate': round(rho_yr_all, 4),
            'spearman_year_dam_rate': round(rho_yr_dam, 4),
            'spearman_year_adj_rate': round(rho_yr_adj, 4),
            'spearman_year_airports': round(rho_yr_air, 4),
            'spearman_airports_total': round(rho_air_tot, 4),
            'spearman_year_dam_frac': round(rho_yr_frac, 4),
        },
        'sensitivity_time_subsets': [
            {
                'label': s['label'],
                'n_years': s['n_years'],
                'bias_linear': round(s['bias_fraction_linear'], 4),
                'bias_sqrt': round(s['bias_fraction_sqrt'], 4),
                'bias_log': round(s['bias_fraction_log'], 4),
            }
            for s in sens
        ],
        'sensitivity_damaging_constancy': [
            {
                'label': d['label'],
                'd_start': d['d_start'],
                'bias_linear': round(d['bias_linear'], 4),
                'bias_sqrt': round(d['bias_sqrt'], 4),
                'bias_log': round(d['bias_log'], 4),
            }
            for d in dam_sens_results
        ],
        'sensitivity_extrapolation': [
            {
                'label': e['label'],
                'rr_2022': round(e['rr_2022'], 3),
                'adj_growth': round(e['adj_growth'], 2),
                'bias_pct': round(e['bias_pct'], 4),
            }
            for e in extrap_results
        ],
        'effect_sizes': {
            'cohens_d_damage_frac_early_late': round(d_frac, 2),
        },
    }

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

    write_report(results, rates, sens, boot, bboot,
                 bias_linear, bias_sqrt, bias_log, bias_adj_log,
                 airport_expansion_share, cf_bias_linear,
                 dam_sens_results, extrap_results)
    print(f"  Wrote {REPORT_FILE}")
    print()

    # ── VERIFICATION ──
    if verify_mode:
        print("=" * 70)
        print("VERIFICATION MODE")
        print("=" * 70)
        print()

        checks = [
            ("Data has 33 years (1990-2022)",
             len(ANNUAL_DATA) == 33),
            ("Data hash is deterministic",
             compute_data_hash() == data_hash),
            ("Total strikes > 250,000",
             total_all > 250000),
            ("Damaging < total for all years",
             all(r[2] < r[1] for r in ANNUAL_DATA)),
            ("All-strike rate slope > 0",
             slope_all > 0),
            ("Damaging slope > 0 but < all slope",
             0 < slope_dam < slope_all),
            ("Adjusted growth factor < raw growth factor",
             adj_ratio < raw_ratio),
            ("Bias fraction (linear) in [0.50, 0.99]",
             0.50 <= bias_linear <= 0.99),
            ("Bias fraction (sqrt) between log and linear",
             bias_log <= bias_sqrt <= bias_linear),
            ("Bias fraction (log) in [0.20, 0.85]",
             0.20 <= bias_log <= 0.85),
            ("Airport expansion share in [0.30, 0.90]",
             0.30 <= airport_expansion_share <= 0.90),
            ("All-strike trend permutation p < 0.01",
             p1 < 0.01),
            ("Ratio trend permutation p < 0.01",
             p2 < 0.01),
            ("Standard bootstrap CI for slope diff excludes 0",
             ci_excludes_zero),
            ("Block bootstrap CI for slope diff excludes 0",
             block_ci_excludes),
            ("Cohen's d (damage frac) < 5.0",
             main_cohens_d < 5.0),
            ("Damaging constancy: bias increases with improvement",
             dam_sens_results[-1]['bias_linear'] >= dam_sens_results[0]['bias_linear']),
            ("results.json and report.md exist",
             os.path.exists(RESULTS_FILE) and os.path.exists(REPORT_FILE)),
        ]

        n_checks = len(checks)
        n_pass = 0
        for i, (desc, passed) in enumerate(checks, 1):
            status = "PASS" if passed else "FAIL"
            print(f"  [{i:2d}/{n_checks}] {status}: {desc}")
            if passed:
                n_pass += 1

        print()
        print(f"  RESULT: {n_pass}/{n_checks} checks passed")
        if n_pass == n_checks:
            print("  STATUS: ALL CHECKS PASSED")
        else:
            print(f"  STATUS: {n_checks - n_pass} CHECK(S) FAILED")
            sys.exit(1)

    print()
    print("ANALYSIS COMPLETE")


# ============================================================
# REPORT GENERATION
# ============================================================

def write_report(results, rates, sens, boot, bboot,
                 bias_linear, bias_sqrt, bias_log, bias_adj_log,
                 airport_expansion_share, cf_bias_linear,
                 dam_sens_results, extrap_results):
    r = results
    lines = []
    a = lines.append

    a("# Bird Strike Reporting Compliance Analysis — Report")
    a("")
    a("## Key Findings")
    a("")
    a(f"Using 33 years of data from the FAA National Wildlife Strike Database "
      f"(1990-2022),")
    a(f"we estimate that **{r['decomposition']['combined_range_pct']}%** of "
      f"the apparent")
    a(f"{r['growth_factors']['raw_rate_growth']}x increase in bird strike "
      f"rates is attributable to")
    a(f"improved voluntary reporting compliance, not genuine increases in "
      f"wildlife encounters.")
    a("")
    a("Three independent methods converge:")
    a("")
    m1 = r['decomposition']['method1_negative_control']
    m2 = r['decomposition']['method2_reporting_adjustment']
    m3 = r['decomposition']['method3_airport_expansion']
    a("| Method | Linear | Sqrt | Log |")
    a("|--------|--------|------|-----|")
    a(f"| 1. Damaging-strike negative control | "
      f"{m1['bias_fraction_linear']*100:.1f}% | "
      f"{m1['bias_fraction_sqrt']*100:.1f}% | "
      f"{m1['bias_fraction_log']*100:.1f}% |")
    a(f"| 2. Reporting-rate adjustment | — | — | "
      f"{m2['bias_fraction_log']*100:.1f}% |")
    a(f"| 3. Airport-expansion decomposition | "
      f"{m3['airport_expansion_share_linear']*100:.1f}% | — | "
      f"{m3['airport_expansion_share_log']*100:.1f}% |")
    a("")

    a("## Data Summary")
    a("")
    a(f"- **Period:** {r['metadata']['data_years']}")
    a(f"- **Total strikes:** {r['metadata']['total_strikes']:,}")
    a(f"- **Total damaging:** {r['metadata']['total_damaging']:,}")
    a(f"- **All-strike rate growth:** {r['growth_factors']['raw_rate_growth']}x")
    a(f"- **Damaging rate growth:** "
      f"{r['growth_factors']['damaging_rate_growth']}x")
    a(f"- **Reporting-adjusted growth:** "
      f"{r['growth_factors']['adjusted_rate_growth']}x")
    a(f"- **Fixed-airport counterfactual growth:** "
      f"{r['growth_factors']['counterfactual_growth']}x")
    a(f"- **Airports reporting growth:** "
      f"{r['growth_factors']['airports_growth']}x")
    a("")

    a("## Statistical Tests")
    a("")
    pt = r['permutation_tests']
    a("| Test | p-value |")
    a("|------|---------|")
    a(f"| All-strike rate positive trend | {fmt_p(pt['all_rate_trend_p'])} |")
    a(f"| Ratio (all/dam) has trend | {fmt_p(pt['ratio_trend_p'])} |")
    a(f"| Damage fraction decline | {fmt_p(pt['dam_frac_decline_p'])} |")
    a(f"| All slope > damaging slope | {fmt_p(pt['slope_diff_p'])} |")
    a(f"| Airports trend positive | {fmt_p(pt['airports_trend_p'])} |")
    a("")

    a("## Bootstrap Confidence Intervals")
    a("")
    bs = r['bootstrap']['standard']
    bb = r['bootstrap']['block']
    slope_diff_pt = r['trends']['slope_all_rate'] - r['trends']['slope_dam_rate']
    a("| Quantity | Point Est. | Standard 95% CI | Block 95% CI |")
    a("|---------|-----------|-----------------|--------------|")
    a(f"| Slope difference | {slope_diff_pt:.4f} | "
      f"[{bs['slope_diff_ci'][0]:.4f}, {bs['slope_diff_ci'][1]:.4f}] | "
      f"[{bb['slope_diff_ci'][0]:.4f}, {bb['slope_diff_ci'][1]:.4f}] |")
    a(f"| Bias fraction (linear) | {bias_linear:.4f} | "
      f"[{bs['bias_frac_ci'][0]:.4f}, {bs['bias_frac_ci'][1]:.4f}] | "
      f"[{bb['bias_frac_ci'][0]:.4f}, {bb['bias_frac_ci'][1]:.4f}] |")
    a(f"| Cohen's d (damage frac) | {bs['cohens_d_damage_frac']:.2f} | — | — |")
    a("")

    a("## Sensitivity: Time-Period Subsets")
    a("")
    a("| Subset | N | Linear | Sqrt | Log |")
    a("|--------|---|--------|------|-----|")
    for s in r['sensitivity_time_subsets']:
        a(f"| {s['label']} | {s['n_years']} | "
          f"{s['bias_linear']*100:.1f}% | "
          f"{s['bias_sqrt']*100:.1f}% | "
          f"{s['bias_log']*100:.1f}% |")
    a("")

    a("## Sensitivity: Damaging-Strike Reporting Constancy")
    a("")
    a("| Scenario | Linear | Sqrt | Log |")
    a("|----------|--------|------|-----|")
    for d in r['sensitivity_damaging_constancy']:
        a(f"| {d['label']} | "
          f"{d['bias_linear']*100:.1f}% | "
          f"{d['bias_sqrt']*100:.1f}% | "
          f"{d['bias_log']*100:.1f}% |")
    a("")
    a("Bias increases under all improvement scenarios, confirming that")
    a("the constant-reporting assumption is conservative.")
    a("")

    a("## Sensitivity: Reporting-Rate Extrapolation")
    a("")
    a("| Scenario | 2022 Rate | Adj Growth | Bias% |")
    a("|----------|-----------|------------|-------|")
    for e in r['sensitivity_extrapolation']:
        a(f"| {e['label']} | {e['rr_2022']*100:.1f}% | "
          f"{e['adj_growth']:.1f}x | {e['bias_pct']*100:.1f}% |")
    a("")

    a("## Limitations")
    a("")
    a("1. Aggregate annual data (national-level, not airport-level microdata)")
    a("2. Damaging-strike constancy assumption (shown to be conservative)")
    a("3. Extrapolated reporting rates post-2014 (sensitivity tested)")
    a("4. No species decomposition (large-bodied bird increases may contribute)")
    a("5. Functional form sensitivity (bracketed with linear/sqrt/log scales)")
    a("")

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


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:** Sectioned output `[1/10]` through `[10/10]`, ending with `ANALYSIS COMPLETE`. Runtime: ~30-60 seconds (permutation and bootstrap tests).

## 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:** 18/18 assertions pass, ending with `STATUS: ALL CHECKS PASSED`.

## Expected Outputs

After successful execution, the workspace contains:

- `analyze.py` — the analysis script (~1220 lines)
- `results.json` — structured results with all metrics, CIs, and p-values
- `report.md` — human-readable summary with tables and findings

## Success Criteria

1. All 10 analysis sections complete without error.
2. `results.json` contains all expected keys (metadata, growth_factors, trends, decomposition, permutation_tests, bootstrap, correlations, sensitivity_time_subsets, sensitivity_damaging_constancy, sensitivity_extrapolation, effect_sizes).
3. All 18 verification assertions pass.
4. Both standard and block bootstrap CIs for the slope difference exclude zero.
5. All permutation test p-values < 0.01.
6. Three independent methods converge on a majority-reporting-bias conclusion.

## Failure Conditions

- Python version < 3.8 (f-strings not available).
- Data validation fails (embedded data was corrupted).
- Any verification assertion fails (indicates a computation error).
- Script exits with non-zero code.