---
name: "pm25-negative-control-confounding"
description: "Uses negative control outcomes (accidents, suicide) to quantify spatial confounding in cross-sectional PM2.5-mortality associations across US states"
version: "2.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "causal-inference", "negative-control", "air-pollution", "epidemiology", "spatial-confounding"]
python_version: ">=3.8"
dependencies: []
---

# Can Negative Control Outcomes Quantify Spatial Confounding in PM2.5-Mortality Studies?

## Research Question

Cross-sectional studies consistently find that higher PM2.5 air pollution is associated
with higher cardiovascular and respiratory mortality. But PM2.5 is spatially correlated
with poverty, smoking, diet, and healthcare access — factors that independently predict
mortality. How much of the PM2.5-mortality association is causal vs. confounded?

**Methodological hook:** We apply a *negative control outcome* design. If PM2.5 truly
causes mortality, it should predict cardiopulmonary deaths but NOT deaths from causes
with no biological pathway from air pollution (accidents, suicide). If PM2.5 predicts
negative controls too, this indicates spatial confounding. The ratio of the negative
control association to the target association quantifies the confounding fraction.

**Enhancements (v2.0):** We add pooled multi-year analysis (2015-2017) for more stable
estimates, partial Spearman correlations controlling for accident rate as an SES proxy,
one-sided differential association tests for greater statistical power, and quintile
dose-response analysis to demonstrate monotonic exposure-response gradients.

## Data Sources

1. **EPA Air Quality System** — Annual PM2.5 concentrations by monitor, aggregated to
   state level. Source: https://aqs.epa.gov/aqsweb/airdata/download_files.html
2. **CDC NCHS Leading Causes of Death** — State-level age-adjusted death rates by cause.
   Source: https://data.cdc.gov (dataset bi63-dtpu)

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_spatial-confounding-in-pm25-studies-negative-control-approac/data
```

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

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_spatial-confounding-in-pm25-studies-negative-control-approac/analyze.py
#!/usr/bin/env python3
"""
Negative Control Outcomes for Quantifying Spatial Confounding
in PM2.5-Mortality Studies — Version 2.0

Improvements over v1.0:
- Pooled multi-year analysis (2015-2017) for more stable estimates
- Partial Spearman correlations controlling for accident rate (SES proxy)
- One-sided differential association tests for greater statistical power
- Quintile dose-response analysis

Data sources:
  - EPA Air Quality System: annual PM2.5 concentrations by monitor
  - CDC NCHS: Leading causes of death by state (age-adjusted rates)

Python 3.8+ standard library only. No third-party packages.
"""

import sys
import os
import json
import csv
import io
import math
import random
import hashlib
import zipfile
import time
import urllib.request
import urllib.error
from collections import defaultdict

# ============================================================
# CONFIGURATION
# ============================================================
SEED = 42
random.seed(SEED)

N_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
ALPHA = 0.05

PRIMARY_YEAR = 2017
SENSITIVITY_YEARS = [2015, 2016]
ALL_YEARS = [2015, 2016, 2017]

DATA_DIR = "data"
RESULTS_FILE = "results.json"
REPORT_FILE = "report.md"

# Target outcomes: biologically plausible causal pathway from PM2.5
TARGET_CAUSES = ["Heart disease", "CLRD"]
# Negative controls: no plausible PM2.5 causal pathway
NEGATIVE_CONTROL_CAUSES = ["Unintentional injuries", "Suicide"]
ALL_CAUSES = TARGET_CAUSES + NEGATIVE_CONTROL_CAUSES

# Data source URLs
EPA_URL_TEMPLATE = (
    "https://aqs.epa.gov/aqsweb/airdata/annual_conc_by_monitor_{year}.zip"
)
CDC_URL = (
    "https://data.cdc.gov/api/views/bi63-dtpu/rows.csv?accessType=DOWNLOAD"
)

# Minimum PM2.5 monitors per state for inclusion
MIN_MONITORS_PER_STATE = 3

# Known SHA256 hashes for data integrity verification
EXPECTED_SHA256 = {
    "annual_conc_by_monitor_2017.zip": "7be35af3811004c0f87adb3917f507b3c94cb5ed68e0253672a4164a6133ea66",
    "annual_conc_by_monitor_2015.zip": "508f00e1679dfc8651e744a36a66fc7c7f5eac0d9ff330023d4b30437491a87f",
    "annual_conc_by_monitor_2016.zip": "8edb9665d0e8da818983f2f1417536e1465e7fb144055811767beeba7178c4d9",
    "cdc_leading_causes_death.csv": "1971e6664c76457a5c32e4c86aa103ba404b7be5d17ddb6e6dfc0075479c85d9",
}


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

def download_with_retry(url, filepath, max_retries=3, timeout=180):
    """Download a file from URL with retry logic."""
    os.makedirs(os.path.dirname(filepath) if os.path.dirname(filepath) else ".",
                exist_ok=True)
    for attempt in range(max_retries):
        try:
            req = urllib.request.Request(url, headers={
                "User-Agent": "Mozilla/5.0 (Python/research-script)"
            })
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                data = resp.read()
            with open(filepath, "wb") as f:
                f.write(data)
            size_mb = len(data) / (1024 * 1024)
            print(f"    Downloaded {size_mb:.1f} MB -> {filepath}")
            return True
        except Exception as e:
            print(f"    Attempt {attempt + 1}/{max_retries} failed: {e}")
            if attempt < max_retries - 1:
                wait = 2 ** (attempt + 1)
                print(f"    Retrying in {wait}s...")
                time.sleep(wait)
    return False


def sha256_file(filepath):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(filepath, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def ensure_data(url, filepath):
    """Download data if not cached; verify SHA256 against known hashes."""
    basename = os.path.basename(filepath)
    expected = EXPECTED_SHA256.get(basename)

    if os.path.exists(filepath):
        h = sha256_file(filepath)
        if expected and h != expected:
            print(f"    WARNING: SHA256 mismatch for {basename}")
            print(f"      Expected: {expected}")
            print(f"      Got:      {h}")
            print(f"    Data source may have been updated. Continuing.")
        else:
            print(f"    Cached: {filepath} (SHA256 verified: {h[:16]}...)")
        return h

    print(f"    Downloading: {url}")
    if not download_with_retry(url, filepath):
        raise RuntimeError(f"FATAL: Failed to download {url} after retries")
    h = sha256_file(filepath)
    if expected and h != expected:
        print(f"    WARNING: SHA256 mismatch for {basename}")
        print(f"      Expected: {expected}")
        print(f"      Got:      {h}")
        print(f"    Data source may have been updated. Continuing.")
    else:
        print(f"    SHA256 verified: {h}")
    return h


# ============================================================
# STATISTICAL FUNCTIONS (pure standard library)
# ============================================================

def mean_val(values):
    """Arithmetic mean."""
    return sum(values) / len(values)


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


def pearson_r(x, y):
    """Pearson product-moment correlation coefficient."""
    n = len(x)
    if n < 3:
        return 0.0
    mx, my = mean_val(x), mean_val(y)
    var_x = sum((xi - mx) ** 2 for xi in x) / (n - 1)
    var_y = sum((yi - my) ** 2 for yi in y) / (n - 1)
    sx = math.sqrt(var_x) if var_x > 0 else 0.0
    sy = math.sqrt(var_y) if var_y > 0 else 0.0
    if sx < 1e-15 or sy < 1e-15:
        return 0.0
    cov = sum((x[i] - mx) * (y[i] - my) for i in range(n)) / (n - 1)
    return max(-1.0, min(1.0, cov / (sx * sy)))


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


def fisher_z_ci(rho, n, alpha=0.05):
    """Fisher z-transform confidence interval for a correlation."""
    if n <= 3:
        return (-1.0, 1.0)
    rho_c = max(-0.9999, min(0.9999, rho))
    z = 0.5 * math.log((1.0 + rho_c) / (1.0 - rho_c))
    se = 1.0 / math.sqrt(n - 3)
    z_crit = 1.96 if alpha == 0.05 else 2.576
    return (math.tanh(z - z_crit * se), math.tanh(z + z_crit * se))


def permutation_test(x, y, n_perms=N_PERMUTATIONS, seed=SEED):
    """Two-sided permutation test for Spearman correlation.

    Returns (observed_rho, p_value).
    """
    rng = random.Random(seed)
    observed = spearman_rho(x, y)
    y_perm = list(y)
    n_extreme = 0
    for _ in range(n_perms):
        rng.shuffle(y_perm)
        if abs(spearman_rho(x, y_perm)) >= abs(observed):
            n_extreme += 1
    p_value = (n_extreme + 1) / (n_perms + 1)
    return observed, p_value


def bootstrap_ci(x, y, n_boot=N_BOOTSTRAP, alpha=0.05, seed=SEED):
    """Bootstrap percentile confidence interval for Spearman correlation."""
    rng = random.Random(seed)
    n = len(x)
    stats = []
    for _ in range(n_boot):
        idx = [rng.randint(0, n - 1) for _ in range(n)]
        stats.append(spearman_rho([x[i] for i in idx], [y[i] for i in idx]))
    stats.sort()
    lo = int(alpha / 2 * n_boot)
    hi = int((1 - alpha / 2) * n_boot) - 1
    return (stats[lo], stats[hi])


def bootstrap_ratio_ci(x, y_target, y_control,
                        n_boot=N_BOOTSTRAP, alpha=0.05, seed=SEED):
    """Bootstrap CI for confounding ratio: rho(x,control) / rho(x,target).

    Skips bootstrap samples where |rho_target| < 0.02 to avoid
    division-by-near-zero instability.
    """
    rng = random.Random(seed)
    n = len(x)
    ratios = []
    for _ in range(n_boot):
        idx = [rng.randint(0, n - 1) for _ in range(n)]
        xb = [x[i] for i in idx]
        ytb = [y_target[i] for i in idx]
        ycb = [y_control[i] for i in idx]
        rt = spearman_rho(xb, ytb)
        rc = spearman_rho(xb, ycb)
        if abs(rt) > 0.02:
            ratios.append(rc / rt)
    if len(ratios) < 100:
        return (float("nan"), float("nan"))
    ratios.sort()
    lo = int(alpha / 2 * len(ratios))
    hi = int((1 - alpha / 2) * len(ratios)) - 1
    return (ratios[lo], ratios[hi])


def bootstrap_difference_ci(x, y_target, y_control,
                             n_boot=N_BOOTSTRAP, alpha=0.05, seed=SEED):
    """Bootstrap CI for rho(x,target) - rho(x,control).

    Returns dict with two-sided CI, one-sided p-value, and one-sided CI lower bound.
    """
    rng = random.Random(seed + 7)  # offset seed for independence
    n = len(x)
    diffs = []
    for _ in range(n_boot):
        idx = [rng.randint(0, n - 1) for _ in range(n)]
        xb = [x[i] for i in idx]
        ytb = [y_target[i] for i in idx]
        ycb = [y_control[i] for i in idx]
        diffs.append(spearman_rho(xb, ytb) - spearman_rho(xb, ycb))
    diffs.sort()
    # Two-sided CI
    lo_2s = int(alpha / 2 * n_boot)
    hi_2s = int((1 - alpha / 2) * n_boot) - 1
    ci_two_sided = (diffs[lo_2s], diffs[hi_2s])
    # One-sided p-value: proportion of bootstrap diffs <= 0
    p_one_sided = sum(1 for d in diffs if d <= 0) / n_boot
    # One-sided 95% CI lower bound
    lo_1s = int(alpha * n_boot)
    ci_one_sided_lower = diffs[lo_1s]
    return {
        "ci_two_sided": ci_two_sided,
        "p_one_sided": p_one_sided,
        "ci_one_sided_lower": ci_one_sided_lower,
    }


# --- NEW v2.0: Partial Spearman correlation ---

def partial_spearman(x, y, z):
    """Partial Spearman correlation rho(x, y | z).

    Uses the standard formula:
    rho_xy.z = (rho_xy - rho_xz * rho_yz) / sqrt((1 - rho_xz^2)(1 - rho_yz^2))
    """
    rho_xy = spearman_rho(x, y)
    rho_xz = spearman_rho(x, z)
    rho_yz = spearman_rho(y, z)
    denom_sq = (1 - rho_xz ** 2) * (1 - rho_yz ** 2)
    if denom_sq <= 1e-10:
        return float("nan")
    return (rho_xy - rho_xz * rho_yz) / math.sqrt(denom_sq)


def bootstrap_partial_ci(x, y, z, n_boot=N_BOOTSTRAP, alpha=0.05, seed=SEED):
    """Bootstrap percentile CI for partial Spearman correlation."""
    rng = random.Random(seed + 11)
    n = len(x)
    stats = []
    for _ in range(n_boot):
        idx = [rng.randint(0, n - 1) for _ in range(n)]
        xb = [x[i] for i in idx]
        yb = [y[i] for i in idx]
        zb = [z[i] for i in idx]
        val = partial_spearman(xb, yb, zb)
        if not math.isnan(val):
            stats.append(val)
    if len(stats) < 100:
        return (float("nan"), float("nan"))
    stats.sort()
    lo = int(alpha / 2 * len(stats))
    hi = int((1 - alpha / 2) * len(stats)) - 1
    return (stats[lo], stats[hi])


def partial_permutation_test(x, y, z, n_perms=N_PERMUTATIONS, seed=SEED):
    """Permutation test for partial Spearman correlation.

    Shuffles y while keeping x and z fixed, testing H0: rho(x,y|z) = 0.
    """
    rng = random.Random(seed + 17)
    observed = partial_spearman(x, y, z)
    if math.isnan(observed):
        return observed, 1.0
    y_perm = list(y)
    n_extreme = 0
    for _ in range(n_perms):
        rng.shuffle(y_perm)
        perm_val = partial_spearman(x, y_perm, z)
        if not math.isnan(perm_val) and abs(perm_val) >= abs(observed):
            n_extreme += 1
    return observed, (n_extreme + 1) / (n_perms + 1)


# --- NEW v2.0: Quintile dose-response ---

def quintile_analysis(x, y):
    """Compute mean y for each quintile of x.

    Returns list of 5 dicts with quintile number, PM2.5 mean, mortality mean, n.
    """
    n = len(x)
    paired = sorted(zip(x, y), key=lambda p: p[0])
    q_size = n // 5
    remainder = n % 5
    quintiles = []
    idx = 0
    for q in range(5):
        size = q_size + (1 if q < remainder else 0)
        q_x = [paired[idx + i][0] for i in range(size)]
        q_y = [paired[idx + i][1] for i in range(size)]
        quintiles.append({
            "quintile": q + 1,
            "pm25_mean": round(mean_val(q_x), 3),
            "mort_mean": round(mean_val(q_y), 2),
            "n": size,
        })
        idx += size
    return quintiles


# ============================================================
# DATA LOADING
# ============================================================

def match_cause(cause_name, target):
    """Flexible cause-name matching to handle minor label variations."""
    if cause_name.strip() == target:
        return True
    cn = cause_name.lower().strip()
    tl = target.lower().strip()
    if tl == "heart disease":
        return "heart disease" in cn
    if tl == "clrd":
        return "clrd" in cn or "chronic lower respiratory" in cn
    if tl == "unintentional injuries":
        return "unintentional" in cn
    if tl == "suicide":
        return "suicide" in cn
    return False


def load_epa_pm25(year):
    """Download EPA PM2.5 monitoring data, aggregate to state-level means.

    Uses Parameter Code 88101 (PM2.5 FRM/FEM, 24-hour).
    Returns dict {state_name: mean_pm25_ug_m3}.
    """
    url = EPA_URL_TEMPLATE.format(year=year)
    zip_path = os.path.join(DATA_DIR, f"annual_conc_by_monitor_{year}.zip")
    ensure_data(url, zip_path)

    # Extract CSV from ZIP
    with zipfile.ZipFile(zip_path, "r") as zf:
        csv_files = [f for f in zf.namelist() if f.endswith(".csv")]
        if not csv_files:
            raise RuntimeError(f"No CSV files found in {zip_path}")
        csv_name = csv_files[0]
        extract_path = os.path.join(DATA_DIR, csv_name)
        if not os.path.exists(extract_path):
            zf.extract(csv_name, DATA_DIR)
            print(f"    Extracted: {csv_name}")

    # Parse CSV: filter for PM2.5 (88101), collect monitor means by state
    state_monitors = defaultdict(list)
    rows_total = 0
    rows_pm25 = 0

    with open(extract_path, "r", encoding="utf-8-sig") as f:
        reader = csv.DictReader(f)
        for row in reader:
            rows_total += 1
            if row.get("Parameter Code", "").strip() != "88101":
                continue
            state = row.get("State Name", "").strip()
            try:
                conc = float(row.get("Arithmetic Mean", "0"))
            except (ValueError, TypeError):
                continue
            if state and conc > 0:
                state_monitors[state].append(conc)
                rows_pm25 += 1

    print(f"    Parsed {rows_total:,} rows total, {rows_pm25:,} PM2.5 monitor records")

    # Aggregate: state mean, require minimum monitor count
    result = {}
    for state, values in state_monitors.items():
        if len(values) >= MIN_MONITORS_PER_STATE:
            result[state] = mean_val(values)

    print(f"    {len(result)} states with >= {MIN_MONITORS_PER_STATE} monitors")
    return result


def load_cdc_mortality():
    """Download CDC Leading Causes of Death data.

    Returns nested dict: {year: {cause: {state: age_adjusted_rate}}}.
    """
    csv_path = os.path.join(DATA_DIR, "cdc_leading_causes_death.csv")
    ensure_data(CDC_URL, csv_path)

    data = defaultdict(lambda: defaultdict(dict))
    rows_total = 0
    rows_used = 0

    with open(csv_path, "r", encoding="utf-8-sig") as f:
        reader = csv.DictReader(f)
        for row in reader:
            rows_total += 1
            state = row.get("State", "").strip()
            if not state or state == "United States":
                continue

            try:
                year = int(row.get("Year", "0"))
            except (ValueError, TypeError):
                continue

            cause_raw = row.get("Cause Name", "").strip()
            rate_str = row.get("Age-adjusted Death Rate", "").strip()

            try:
                rate = float(rate_str)
            except (ValueError, TypeError):
                continue

            for target_cause in ALL_CAUSES:
                if match_cause(cause_raw, target_cause):
                    data[year][target_cause][state] = rate
                    rows_used += 1
                    break

    print(f"    Parsed {rows_total:,} rows total, {rows_used:,} matched records")
    years = sorted(data.keys())
    if years:
        print(f"    Years available: {years[0]}-{years[-1]}")
    return data


# ============================================================
# ANALYSIS
# ============================================================

def analyze_association(pm25_data, mort_data, cause, year):
    """Compute PM2.5-mortality association for one cause-year pair.

    Returns dict with correlation, CI, p-value, and matched data,
    or None if fewer than 15 states match.
    """
    cause_data = mort_data.get(year, {}).get(cause, {})
    common = sorted(set(pm25_data.keys()) & set(cause_data.keys()))
    if len(common) < 15:
        return None

    x = [pm25_data[s] for s in common]
    y = [cause_data[s] for s in common]

    rho = spearman_rho(x, y)
    ci_fisher = fisher_z_ci(rho, len(common))
    _, p_perm = permutation_test(x, y)
    ci_boot = bootstrap_ci(x, y)

    return {
        "cause": cause,
        "year": year,
        "n_states": len(common),
        "states": common,
        "rho": rho,
        "fisher_ci": ci_fisher,
        "bootstrap_ci": ci_boot,
        "perm_p": p_perm,
        "pm25": x,
        "mort": y,
    }


def analyze_pooled(pm25_by_year, mort_data, cause, years):
    """Pooled multi-year analysis: average PM2.5 and mortality across years.

    Returns same structure as analyze_association or None.
    """
    # Find states present in all years
    state_sets = []
    for year in years:
        pm_states = set(pm25_by_year[year].keys())
        cause_states = set(mort_data.get(year, {}).get(cause, {}).keys())
        state_sets.append(pm_states & cause_states)
    if not state_sets:
        return None
    common = sorted(set.intersection(*state_sets))
    if len(common) < 15:
        return None

    # Average across years
    x = [mean_val([pm25_by_year[y][s] for y in years]) for s in common]
    y = [mean_val([mort_data[y][cause][s] for y in years]) for s in common]

    rho = spearman_rho(x, y)
    ci_fisher = fisher_z_ci(rho, len(common))
    _, p_perm = permutation_test(x, y, seed=SEED + 100)
    ci_boot = bootstrap_ci(x, y, seed=SEED + 100)

    return {
        "cause": cause,
        "years": years,
        "n_states": len(common),
        "states": common,
        "rho": rho,
        "fisher_ci": ci_fisher,
        "bootstrap_ci": ci_boot,
        "perm_p": p_perm,
        "pm25": x,
        "mort": y,
    }


# ============================================================
# REPORT
# ============================================================

def write_report(results):
    """Write human-readable Markdown report."""
    meta = results["metadata"]
    lines = [
        "# Negative Control Outcomes for Spatial Confounding "
        "in PM2.5-Mortality Studies (v2.0)\n",
        "## Summary\n",
        f"- **Primary year:** {meta['primary_year']}",
        f"- **Pooled years:** {meta['pooled_years']}",
        f"- **Sensitivity years:** {meta['sensitivity_years']}",
        f"- **Permutations:** {meta['n_permutations']}",
        f"- **Bootstrap samples:** {meta['n_bootstrap']}",
        f"- **Random seed:** {meta['random_seed']}",
        f"- **Significance level:** {meta['alpha']}\n",
        "## Primary Results\n",
        "| Cause | Role | Spearman rho | 95% Boot. CI | Perm. p | N |",
        "|---|---|---|---|---|---|",
    ]
    for cause in ALL_CAUSES:
        r = results["primary_results"].get(cause)
        if not r:
            continue
        role = "Target" if r["role"] == "target" else "Neg. Control"
        lines.append(
            f"| {cause} | {role} | {r['rho']:.4f} | "
            f"[{r['bootstrap_ci'][0]:.4f}, {r['bootstrap_ci'][1]:.4f}] | "
            f"{r['perm_p']:.4f} | {r['n_states']} |"
        )

    lines += [
        "\n## Confounding Ratios\n",
        "Confounding ratio = rho(PM2.5, negative control) / "
        "rho(PM2.5, target outcome).\n",
        "| Comparison | Ratio | 95% Boot. CI | rho_target | rho_control | N |",
        "|---|---|---|---|---|---|",
    ]
    for key, cr in results["confounding_ratios"].items():
        rv = cr["ratio"]
        r_str = f"{rv:.4f}" if rv is not None else "N/A"
        ci = cr["ci"]
        ci_str = (f"[{ci[0]:.4f}, {ci[1]:.4f}]"
                  if ci[0] is not None else "N/A")
        lines.append(
            f"| {key} | {r_str} | {ci_str} | "
            f"{cr['rho_target']:.4f} | {cr['rho_control']:.4f} | "
            f"{cr['n_states']} |"
        )

    lines += [
        "\n## Differential Association Test\n",
        "Tests whether rho(PM2.5, target) > rho(PM2.5, negative control).\n",
        "| Comparison | Diff | 95% 2-sided CI | Excl. 0? | 1-sided p | 1-sided lower |",
        "|---|---|---|---|---|---|",
    ]
    for key, dt in results.get("differential_tests", {}).items():
        ci = dt["ci_two_sided"]
        lines.append(
            f"| {key} | {dt['diff']:.4f} | "
            f"[{ci[0]:.4f}, {ci[1]:.4f}] | "
            f"{'YES' if dt['ci_excludes_zero'] else 'NO'} | "
            f"{dt['p_one_sided']:.4f} | {dt['ci_one_sided_lower']:.4f} |"
        )

    # Pooled analysis
    lines += [
        "\n## Pooled Multi-Year Analysis (2015-2017 averaged)\n",
        "| Cause | Role | Spearman rho | 95% Boot. CI | Perm. p | N |",
        "|---|---|---|---|---|---|",
    ]
    for cause in ALL_CAUSES:
        r = results.get("pooled_analysis", {}).get(cause)
        if not r:
            continue
        role = "Target" if r["role"] == "target" else "Neg. Control"
        lines.append(
            f"| {cause} | {role} | {r['rho']:.4f} | "
            f"[{r['bootstrap_ci'][0]:.4f}, {r['bootstrap_ci'][1]:.4f}] | "
            f"{r['perm_p']:.4f} | {r['n_states']} |"
        )

    # Partial correlations
    lines += [
        "\n## Partial Spearman Correlations (controlling for accident rate)\n",
        "| Target | Partial rho(PM2.5, target | accidents) | 95% Boot. CI | Perm. p | Raw rho |",
        "|---|---|---|---|---|",
    ]
    for key, pc in results.get("partial_correlations", {}).items():
        lines.append(
            f"| {key} | {pc['partial_rho']:.4f} | "
            f"[{pc['bootstrap_ci'][0]:.4f}, {pc['bootstrap_ci'][1]:.4f}] | "
            f"{pc['perm_p']:.4f} | {pc['raw_rho']:.4f} |"
        )

    # Quintile analysis
    lines += [
        "\n## Quintile Dose-Response (PM2.5 quintiles, 2017)\n",
    ]
    for cause in ALL_CAUSES:
        qa = results.get("quintile_analysis", {}).get(cause)
        if not qa:
            continue
        lines.append(f"\n**{cause}:**\n")
        lines.append("| Quintile | PM2.5 (ug/m3) | Death rate | N |")
        lines.append("|---|---|---|---|")
        for q in qa:
            lines.append(
                f"| Q{q['quintile']} | {q['pm25_mean']:.1f} | "
                f"{q['mort_mean']:.1f} | {q['n']} |"
            )

    lines += [
        "\n## Temporal Sensitivity\n",
        "| Year | Cause | rho | p | N |",
        "|---|---|---|---|---|",
    ]
    for yr_s, yr_d in sorted(results["sensitivity_temporal"].items()):
        for cause in ALL_CAUSES:
            if cause in yr_d:
                r = yr_d[cause]
                lines.append(
                    f"| {yr_s} | {cause} | {r['rho']:.4f} | "
                    f"{r['perm_p']:.4f} | {r['n_states']} |"
                )

    lines += [
        "\n## Geographic Sensitivity (excl. Alaska, Hawaii)\n",
        "| Cause | rho | p | N |",
        "|---|---|---|---|",
    ]
    for cause in ALL_CAUSES:
        if cause in results["sensitivity_geographic"]:
            r = results["sensitivity_geographic"][cause]
            lines.append(
                f"| {cause} | {r['rho']:.4f} | "
                f"{r['perm_p']:.4f} | {r['n_states']} |"
            )

    lines += [
        "\n## Interpretation\n",
        "If PM2.5 significantly predicts negative control outcomes (accidents, "
        "suicide), the cross-sectional PM2.5-mortality association is at least "
        "partly attributable to spatial confounding by socioeconomic and "
        "behavioral factors (poverty, smoking, diet, healthcare access) that "
        "correlate with both air pollution and all-cause mortality.\n",
        "The confounding ratio quantifies what fraction of the target "
        "association could be explained by such confounding. The partial "
        "correlation analysis provides a complementary check by showing whether "
        "the PM2.5-target association persists after adjusting for the negative "
        "control outcome as an SES proxy.\n",
        "## Limitations\n",
        "1. **Ecological fallacy:** State-level associations may not reflect "
        "individual-level effects.",
        "2. **Exposure misclassification:** State-mean PM2.5 from fixed "
        "monitors may not represent true population exposure.",
        "3. **Imperfect negative controls:** Accidents or suicide may have "
        "weak indirect links to air quality or shared upstream causes.",
        "4. **Limited spatial resolution:** State-level analysis masks "
        "within-state urban/rural heterogeneity.",
        "5. **Cross-sectional design:** Cannot establish temporal ordering "
        "or rule out reverse causation.",
        "6. **Unmeasured effect modification:** The confounding ratio assumes "
        "comparable confounding structure for target and control outcomes.",
        "7. **No direct covariate measurement:** Partial correlations use "
        "negative control outcomes as SES proxies rather than direct measures.",
    ]

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


# ============================================================
# MAIN
# ============================================================

def main():
    N = 11
    print("=" * 70)
    print("NEGATIVE CONTROL OUTCOMES FOR SPATIAL CONFOUNDING")
    print("IN PM2.5-MORTALITY STUDIES (v2.0)")
    print("=" * 70)

    os.makedirs(DATA_DIR, exist_ok=True)

    # ---- [1/11] Download EPA PM2.5 ----
    print(f"\n[1/{N}] Downloading EPA PM2.5 monitoring data...")
    pm25_by_year = {}
    for year in ALL_YEARS:
        print(f"  Year {year}:")
        pm25_by_year[year] = load_epa_pm25(year)

    # ---- [2/11] Download CDC mortality ----
    print(f"\n[2/{N}] Downloading CDC cause-specific mortality data...")
    mort_data = load_cdc_mortality()

    # ---- [3/11] Primary analysis ----
    print(f"\n[3/{N}] Primary analysis: PM2.5 vs cause-specific mortality "
          f"({PRIMARY_YEAR})...")
    primary = {}
    for cause in ALL_CAUSES:
        result = analyze_association(
            pm25_by_year[PRIMARY_YEAR], mort_data, cause, PRIMARY_YEAR
        )
        if result:
            primary[cause] = result
            role = "TARGET" if cause in TARGET_CAUSES else "NEG CTRL"
            print(f"  [{role}] {cause}:")
            print(f"    Spearman rho  = {result['rho']:.4f}")
            print(f"    Bootstrap CI  = [{result['bootstrap_ci'][0]:.4f}, "
                  f"{result['bootstrap_ci'][1]:.4f}]")
            print(f"    Fisher CI     = [{result['fisher_ci'][0]:.4f}, "
                  f"{result['fisher_ci'][1]:.4f}]")
            print(f"    Perm. p-value = {result['perm_p']:.4f} "
                  f"({N_PERMUTATIONS} permutations)")
            print(f"    N states      = {result['n_states']}")
        else:
            print(f"  {cause}: INSUFFICIENT DATA (< 15 matched states)")

    # ---- [4/11] Confounding ratios ----
    print(f"\n[4/{N}] Computing confounding ratios with bootstrap CIs...")
    conf_ratios = {}
    for target in TARGET_CAUSES:
        for control in NEGATIVE_CONTROL_CAUSES:
            if target not in primary or control not in primary:
                continue
            tr = primary[target]
            cr = primary[control]
            common = sorted(set(tr["states"]) & set(cr["states"]))
            if len(common) < 15:
                continue

            pm_m = [pm25_by_year[PRIMARY_YEAR][s] for s in common]
            tgt_m = [mort_data[PRIMARY_YEAR][target][s] for s in common]
            ctl_m = [mort_data[PRIMARY_YEAR][control][s] for s in common]

            rho_t = spearman_rho(pm_m, tgt_m)
            rho_c = spearman_rho(pm_m, ctl_m)

            if abs(rho_t) > 0.02:
                ratio = rho_c / rho_t
            else:
                ratio = float("nan")

            ratio_ci = bootstrap_ratio_ci(pm_m, tgt_m, ctl_m)

            key = f"{control} / {target}"
            conf_ratios[key] = {
                "ratio": ratio,
                "ci": ratio_ci,
                "rho_target": rho_t,
                "rho_control": rho_c,
                "n_states": len(common),
            }
            print(f"  {key}:")
            print(f"    Ratio         = {ratio:.4f}")
            ci_lo, ci_hi = ratio_ci
            if not math.isnan(ci_lo):
                print(f"    95% Boot. CI  = [{ci_lo:.4f}, {ci_hi:.4f}]")
            print(f"    rho(target)   = {rho_t:.4f}")
            print(f"    rho(control)  = {rho_c:.4f}")
            print(f"    N states      = {len(common)}")

    # ---- [5/11] Differential association test (two-sided + one-sided) ----
    print(f"\n[5/{N}] Differential association test: target vs negative control...")
    diff_tests = {}
    for target in TARGET_CAUSES:
        for control in NEGATIVE_CONTROL_CAUSES:
            if target not in primary or control not in primary:
                continue
            tr = primary[target]
            cr = primary[control]
            common = sorted(set(tr["states"]) & set(cr["states"]))
            if len(common) < 15:
                continue
            pm_m = [pm25_by_year[PRIMARY_YEAR][s] for s in common]
            tgt_m = [mort_data[PRIMARY_YEAR][target][s] for s in common]
            ctl_m = [mort_data[PRIMARY_YEAR][control][s] for s in common]
            rho_t = spearman_rho(pm_m, tgt_m)
            rho_c = spearman_rho(pm_m, ctl_m)
            obs_diff = rho_t - rho_c
            diff_result = bootstrap_difference_ci(pm_m, tgt_m, ctl_m)
            ci_2s = diff_result["ci_two_sided"]
            p_1s = diff_result["p_one_sided"]
            ci_1s_lo = diff_result["ci_one_sided_lower"]
            key = f"{target} - {control}"
            excludes_zero = (ci_2s[0] > 0) or (ci_2s[1] < 0)
            diff_tests[key] = {
                "diff": obs_diff,
                "ci_two_sided": ci_2s,
                "ci_excludes_zero": excludes_zero,
                "p_one_sided": p_1s,
                "ci_one_sided_lower": ci_1s_lo,
                "rho_target": rho_t,
                "rho_control": rho_c,
                "n_states": len(common),
            }
            print(f"  {key}:")
            print(f"    rho_target - rho_control = {obs_diff:.4f}")
            print(f"    95% 2-sided CI = [{ci_2s[0]:.4f}, {ci_2s[1]:.4f}]")
            print(f"    CI excludes 0: {'YES' if excludes_zero else 'NO'}")
            print(f"    1-sided p = {p_1s:.4f}")
            print(f"    1-sided 95% CI lower = {ci_1s_lo:.4f}")

    # ---- [6/11] Pooled multi-year analysis ----
    print(f"\n[6/{N}] Pooled multi-year analysis ({ALL_YEARS})...")
    pooled = {}
    for cause in ALL_CAUSES:
        result = analyze_pooled(pm25_by_year, mort_data, cause, ALL_YEARS)
        if result:
            pooled[cause] = result
            role = "TARGET" if cause in TARGET_CAUSES else "NEG CTRL"
            print(f"  [{role}] {cause}:")
            print(f"    Spearman rho  = {result['rho']:.4f}")
            print(f"    Bootstrap CI  = [{result['bootstrap_ci'][0]:.4f}, "
                  f"{result['bootstrap_ci'][1]:.4f}]")
            print(f"    Perm. p-value = {result['perm_p']:.4f}")
            print(f"    N states      = {result['n_states']}")

    # ---- [7/11] Partial Spearman correlations ----
    print(f"\n[7/{N}] Partial Spearman correlations "
          f"(controlling for accident rate as SES proxy)...")
    partial_results = {}
    # Use accident rate as SES proxy covariate
    control_cause = "Unintentional injuries"
    if control_cause in primary:
        for target in TARGET_CAUSES:
            if target not in primary:
                continue
            tr = primary[target]
            cr = primary[control_cause]
            common = sorted(set(tr["states"]) & set(cr["states"]))
            if len(common) < 15:
                continue
            pm_m = [pm25_by_year[PRIMARY_YEAR][s] for s in common]
            tgt_m = [mort_data[PRIMARY_YEAR][target][s] for s in common]
            ctl_m = [mort_data[PRIMARY_YEAR][control_cause][s] for s in common]
            p_rho = partial_spearman(pm_m, tgt_m, ctl_m)
            p_ci = bootstrap_partial_ci(pm_m, tgt_m, ctl_m)
            _, p_perm = partial_permutation_test(pm_m, tgt_m, ctl_m)
            raw_rho = spearman_rho(pm_m, tgt_m)
            key = target
            partial_results[key] = {
                "partial_rho": p_rho,
                "bootstrap_ci": p_ci,
                "perm_p": p_perm,
                "raw_rho": raw_rho,
                "controlling_for": control_cause,
                "n_states": len(common),
            }
            print(f"  {target} | {control_cause}:")
            print(f"    Partial rho   = {p_rho:.4f}")
            if not math.isnan(p_ci[0]):
                print(f"    Bootstrap CI  = [{p_ci[0]:.4f}, {p_ci[1]:.4f}]")
            print(f"    Perm. p-value = {p_perm:.4f}")
            print(f"    Raw rho       = {raw_rho:.4f}")
            print(f"    N states      = {len(common)}")

    # Also partial correlations controlling for suicide
    control_cause_2 = "Suicide"
    if control_cause_2 in primary:
        for target in TARGET_CAUSES:
            if target not in primary:
                continue
            tr = primary[target]
            cr = primary[control_cause_2]
            common = sorted(set(tr["states"]) & set(cr["states"]))
            if len(common) < 15:
                continue
            pm_m = [pm25_by_year[PRIMARY_YEAR][s] for s in common]
            tgt_m = [mort_data[PRIMARY_YEAR][target][s] for s in common]
            ctl_m = [mort_data[PRIMARY_YEAR][control_cause_2][s] for s in common]
            p_rho = partial_spearman(pm_m, tgt_m, ctl_m)
            p_ci = bootstrap_partial_ci(pm_m, tgt_m, ctl_m, seed=SEED + 23)
            _, p_perm = partial_permutation_test(pm_m, tgt_m, ctl_m, seed=SEED + 23)
            raw_rho = spearman_rho(pm_m, tgt_m)
            key = f"{target} | {control_cause_2}"
            partial_results[key] = {
                "partial_rho": p_rho,
                "bootstrap_ci": p_ci,
                "perm_p": p_perm,
                "raw_rho": raw_rho,
                "controlling_for": control_cause_2,
                "n_states": len(common),
            }
            print(f"  {target} | {control_cause_2}:")
            print(f"    Partial rho   = {p_rho:.4f}")
            if not math.isnan(p_ci[0]):
                print(f"    Bootstrap CI  = [{p_ci[0]:.4f}, {p_ci[1]:.4f}]")
            print(f"    Perm. p-value = {p_perm:.4f}")
            print(f"    Raw rho       = {raw_rho:.4f}")
            print(f"    N states      = {len(common)}")

    # ---- [8/11] Quintile dose-response ----
    print(f"\n[8/{N}] Quintile dose-response analysis ({PRIMARY_YEAR})...")
    quintile_results = {}
    for cause in ALL_CAUSES:
        if cause not in primary:
            continue
        r = primary[cause]
        quints = quintile_analysis(r["pm25"], r["mort"])
        quintile_results[cause] = quints
        role = "TARGET" if cause in TARGET_CAUSES else "NEG CTRL"
        print(f"  [{role}] {cause}:")
        for q in quints:
            print(f"    Q{q['quintile']}: PM2.5={q['pm25_mean']:.1f}, "
                  f"DeathRate={q['mort_mean']:.1f}, n={q['n']}")

    # ---- [9/11] Temporal sensitivity ----
    print(f"\n[9/{N}] Sensitivity analysis: temporal stability...")
    sensitivity_temporal = {}
    for year in SENSITIVITY_YEARS:
        print(f"  Year {year}:")
        yearly = {}
        for cause in ALL_CAUSES:
            result = analyze_association(
                pm25_by_year[year], mort_data, cause, year
            )
            if result:
                yearly[cause] = result
                role = "TARGET" if cause in TARGET_CAUSES else "NEG CTRL"
                print(f"    [{role}] {cause}: rho={result['rho']:.4f}, "
                      f"p={result['perm_p']:.4f}, n={result['n_states']}")
        sensitivity_temporal[year] = yearly

    # ---- [10/11] Geographic sensitivity ----
    print(f"\n[10/{N}] Sensitivity: excluding Alaska and Hawaii...")
    exclude = {"Alaska", "Hawaii"}
    pm25_excl = {k: v for k, v in pm25_by_year[PRIMARY_YEAR].items()
                 if k not in exclude}
    geo_results = {}
    for cause in ALL_CAUSES:
        result = analyze_association(pm25_excl, mort_data, cause, PRIMARY_YEAR)
        if result:
            geo_results[cause] = result
            orig = primary.get(cause, {})
            orig_rho = orig.get("rho", float("nan")) if isinstance(orig, dict) else float("nan")
            role = "TARGET" if cause in TARGET_CAUSES else "NEG CTRL"
            print(f"  [{role}] {cause}: rho={result['rho']:.4f} "
                  f"(was {orig_rho:.4f} with all states), n={result['n_states']}")

    # ---- [11/11] Write outputs ----
    print(f"\n[11/{N}] Writing results...")

    results = {
        "metadata": {
            "analysis": "PM2.5 Negative Control Confounding Analysis v2.0",
            "primary_year": PRIMARY_YEAR,
            "pooled_years": ALL_YEARS,
            "sensitivity_years": SENSITIVITY_YEARS,
            "n_permutations": N_PERMUTATIONS,
            "n_bootstrap": N_BOOTSTRAP,
            "random_seed": SEED,
            "alpha": ALPHA,
            "min_monitors_per_state": MIN_MONITORS_PER_STATE,
        },
        "primary_results": {},
        "confounding_ratios": {},
        "differential_tests": {},
        "pooled_analysis": {},
        "partial_correlations": {},
        "quintile_analysis": {},
        "sensitivity_temporal": {},
        "sensitivity_geographic": {},
    }

    for cause, r in primary.items():
        results["primary_results"][cause] = {
            "rho": round(r["rho"], 4),
            "fisher_ci": [round(v, 4) for v in r["fisher_ci"]],
            "bootstrap_ci": [round(v, 4) for v in r["bootstrap_ci"]],
            "perm_p": round(r["perm_p"], 4),
            "n_states": r["n_states"],
            "role": "target" if cause in TARGET_CAUSES else "negative_control",
        }

    for key, cr in conf_ratios.items():
        r_val = cr["ratio"]
        ci = cr["ci"]
        results["confounding_ratios"][key] = {
            "ratio": round(r_val, 4) if not math.isnan(r_val) else None,
            "ci": [
                round(v, 4) if not math.isnan(v) else None
                for v in ci
            ],
            "rho_target": round(cr["rho_target"], 4),
            "rho_control": round(cr["rho_control"], 4),
            "n_states": cr["n_states"],
        }

    for key, dt in diff_tests.items():
        ci_2s = dt["ci_two_sided"]
        results["differential_tests"][key] = {
            "diff": round(dt["diff"], 4),
            "ci_two_sided": [round(v, 4) for v in ci_2s],
            "ci_excludes_zero": dt["ci_excludes_zero"],
            "p_one_sided": round(dt["p_one_sided"], 4),
            "ci_one_sided_lower": round(dt["ci_one_sided_lower"], 4),
            "rho_target": round(dt["rho_target"], 4),
            "rho_control": round(dt["rho_control"], 4),
            "n_states": dt["n_states"],
        }

    for cause, r in pooled.items():
        results["pooled_analysis"][cause] = {
            "rho": round(r["rho"], 4),
            "fisher_ci": [round(v, 4) for v in r["fisher_ci"]],
            "bootstrap_ci": [round(v, 4) for v in r["bootstrap_ci"]],
            "perm_p": round(r["perm_p"], 4),
            "n_states": r["n_states"],
            "role": "target" if cause in TARGET_CAUSES else "negative_control",
        }

    for key, pc in partial_results.items():
        ci = pc["bootstrap_ci"]
        results["partial_correlations"][key] = {
            "partial_rho": round(pc["partial_rho"], 4),
            "bootstrap_ci": [
                round(v, 4) if not math.isnan(v) else None for v in ci
            ],
            "perm_p": round(pc["perm_p"], 4),
            "raw_rho": round(pc["raw_rho"], 4),
            "controlling_for": pc["controlling_for"],
            "n_states": pc["n_states"],
        }

    for cause, quints in quintile_results.items():
        results["quintile_analysis"][cause] = quints

    for year, yearly in sensitivity_temporal.items():
        results["sensitivity_temporal"][str(year)] = {}
        for cause, r in yearly.items():
            results["sensitivity_temporal"][str(year)][cause] = {
                "rho": round(r["rho"], 4),
                "perm_p": round(r["perm_p"], 4),
                "n_states": r["n_states"],
            }

    for cause, r in geo_results.items():
        results["sensitivity_geographic"][cause] = {
            "rho": round(r["rho"], 4),
            "perm_p": round(r["perm_p"], 4),
            "n_states": r["n_states"],
        }

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

    write_report(results)
    print(f"  Written: {REPORT_FILE}")

    print("\n" + "=" * 70)
    print("ANALYSIS COMPLETE")
    print("=" * 70)


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

def verify():
    """Run machine-checkable verification assertions on results."""
    print("=" * 70)
    print("VERIFICATION MODE")
    print("=" * 70)

    passed = 0
    total = 0

    def check(cond, desc):
        nonlocal passed, total
        total += 1
        status = "PASS" if cond else "FAIL"
        if cond:
            passed += 1
        print(f"  {status} [{total:02d}]: {desc}")

    # 1. results.json exists
    check(os.path.exists(RESULTS_FILE), "results.json exists")

    with open(RESULTS_FILE) as f:
        res = json.load(f)

    # 2-4. Required keys
    check("metadata" in res, "results.json has 'metadata'")
    check("primary_results" in res, "results.json has 'primary_results'")
    check("confounding_ratios" in res, "results.json has 'confounding_ratios'")

    # 5-6. Target and negative control presence
    targets = [k for k, v in res["primary_results"].items()
               if v.get("role") == "target"]
    controls = [k for k, v in res["primary_results"].items()
                if v.get("role") == "negative_control"]
    check(len(targets) >= 1,
          f">=1 target cause ({len(targets)} found: {targets})")
    check(len(controls) >= 1,
          f">=1 negative control ({len(controls)} found: {controls})")

    # 7. All correlations in [-1, 1]
    rhos = [v["rho"] for v in res["primary_results"].values()]
    check(all(-1.0 <= r <= 1.0 for r in rhos),
          f"All rho in [-1,1] (values: {[round(r,3) for r in rhos]})")

    # 8. All p-values in [0, 1]
    ps = [v["perm_p"] for v in res["primary_results"].values()]
    check(all(0 <= p <= 1 for p in ps),
          f"All p-values in [0,1] (values: {[round(p,4) for p in ps]})")

    # 9. Sufficient sample size
    ns = [v["n_states"] for v in res["primary_results"].values()]
    check(all(n >= 30 for n in ns),
          f">=30 states per analysis (min={min(ns)}, max={max(ns)})")

    # 10. Confounding ratios computed
    check(len(res["confounding_ratios"]) >= 1,
          f">=1 confounding ratio ({len(res['confounding_ratios'])} found)")

    # 11-14. Bootstrap CIs properly ordered for each primary result
    for cause, r in res["primary_results"].items():
        ci = r.get("bootstrap_ci", [0, 0])
        check(ci[0] <= ci[1],
              f"Bootstrap CI ordered for {cause}: [{ci[0]:.4f}, {ci[1]:.4f}]")

    # 15. Differential tests present
    check(len(res.get("differential_tests", {})) >= 1,
          f"Differential association tests present "
          f"({len(res.get('differential_tests', {}))} found)")

    # 16. One-sided p-values in differential tests
    one_sided_ps = [v.get("p_one_sided", -1)
                    for v in res.get("differential_tests", {}).values()]
    check(all(0 <= p <= 1 for p in one_sided_ps),
          f"One-sided p-values in [0,1] (values: {[round(p,4) for p in one_sided_ps]})")

    # 17. Temporal sensitivity present
    check(len(res.get("sensitivity_temporal", {})) >= 1,
          "Temporal sensitivity analysis present")

    # 18. Geographic sensitivity present
    check(len(res.get("sensitivity_geographic", {})) >= 1,
          "Geographic sensitivity analysis present")

    # 19. report.md exists
    check(os.path.exists(REPORT_FILE), "report.md exists")

    # 20-21. Metadata correctness
    meta = res["metadata"]
    check(meta.get("random_seed") == SEED,
          f"Random seed = {SEED}")
    check(meta.get("n_permutations") == N_PERMUTATIONS,
          f"N permutations = {N_PERMUTATIONS}")

    # 22. Pooled analysis present
    check(len(res.get("pooled_analysis", {})) >= 1,
          f"Pooled analysis present "
          f"({len(res.get('pooled_analysis', {}))} causes)")

    # 23. Pooled analysis has target results
    pooled_targets = [k for k, v in res.get("pooled_analysis", {}).items()
                      if v.get("role") == "target"]
    check(len(pooled_targets) >= 1,
          f"Pooled analysis has >=1 target ({pooled_targets})")

    # 24. Pooled correlations in [-1,1]
    pooled_rhos = [v["rho"] for v in res.get("pooled_analysis", {}).values()]
    check(all(-1.0 <= r <= 1.0 for r in pooled_rhos),
          f"Pooled rho in [-1,1] (values: {[round(r,3) for r in pooled_rhos]})")

    # 25. Partial correlations present
    check(len(res.get("partial_correlations", {})) >= 1,
          f"Partial correlations present "
          f"({len(res.get('partial_correlations', {}))} found)")

    # 26. Partial correlation values in [-1,1]
    partial_rhos = [v["partial_rho"]
                    for v in res.get("partial_correlations", {}).values()]
    check(all(-1.0 <= r <= 1.0 for r in partial_rhos),
          f"Partial rho in [-1,1] (values: {[round(r,3) for r in partial_rhos]})")

    # 27. Quintile analysis present
    check(len(res.get("quintile_analysis", {})) >= 1,
          f"Quintile analysis present "
          f"({len(res.get('quintile_analysis', {}))} causes)")

    # 28. Quintile analysis has 5 quintiles per cause
    for cause, quints in res.get("quintile_analysis", {}).items():
        check(len(quints) == 5,
              f"Quintile analysis: {cause} has 5 quintiles")
        break  # just check one

    print(f"\n{'=' * 70}")
    print(f"VERIFICATION: {passed}/{total} assertions passed")
    if passed == total:
        print("ALL CHECKS PASSED")
    else:
        print(f"FAILED: {total - passed} assertion(s) did not pass")
        sys.exit(1)


# ============================================================
# ENTRY POINT
# ============================================================

if __name__ == "__main__":
    if "--verify" in sys.argv:
        verify()
    else:
        main()
SCRIPT_EOF
```

**Expected output:** No output (script file created).

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_spatial-confounding-in-pm25-studies-negative-control-approac && python3 analyze.py
```

**Expected output:**
- Sectioned progress `[1/11]` through `[11/11]`
- For each cause: Spearman rho, 95% CIs, permutation p-value, N states
- Confounding ratios with bootstrap CIs
- Differential association tests with two-sided CIs and one-sided p-values
- Pooled multi-year analysis results
- Partial Spearman correlations controlling for accident rate
- Quintile dose-response tables
- Sensitivity results for alternate years and geographic subsets
- Final line: `ANALYSIS COMPLETE`
- Creates `results.json` and `report.md` in the workspace

**Expected runtime:** 5-20 minutes (data download + permutation/bootstrap tests).

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_spatial-confounding-in-pm25-studies-negative-control-approac && python3 analyze.py --verify
```

**Expected output:**
- 28+ PASS assertions checking result structure, value ranges, and completeness
- Final line: `ALL CHECKS PASSED`
- Exit code 0

## Expected Outputs

| File | Description |
|---|---|
| `data/annual_conc_by_monitor_*.zip` | Cached EPA PM2.5 monitoring data (ZIP) |
| `data/cdc_leading_causes_death.csv` | Cached CDC mortality data |
| `results.json` | Structured results: correlations, p-values, CIs, confounding ratios, partial correlations, quintile analysis |
| `report.md` | Human-readable Markdown report with tables and interpretation |

## Success Criteria

1. Script downloads EPA and CDC data without errors
2. Prints `ANALYSIS COMPLETE` with all 11 sections completing
3. `results.json` contains primary associations, confounding ratios, differential tests, pooled analysis, partial correlations, quintile analysis, and sensitivity analyses
4. `report.md` contains readable tables and findings
5. Verification passes all 28+ assertions
6. All data cached with SHA256 integrity verification
7. Correlations in [-1, 1], p-values in [0, 1], CIs properly ordered
8. At least 30 matched states per analysis

## Failure Conditions

- Python < 3.8
- Network failure preventing download from EPA AQS or CDC data.cdc.gov
- CDC dataset bi63-dtpu removed or restructured (fatal: analysis requires real mortality data)
- EPA annual_conc_by_monitor file format changed (fatal: need PM2.5 concentrations)
- Fewer than 15 matched states for any cause (analysis skipped for that cause)

## Generalization Guide

This negative control framework generalizes to any exposure-outcome association:

1. **Different exposures:** Replace PM2.5 with lead, noise, temperature, ozone, etc. Change the EPA data source accordingly.
2. **Different target outcomes:** Replace heart disease/CLRD with asthma hospitalizations, stroke, lung cancer, etc.
3. **Different negative controls:** Any outcome with no plausible causal pathway from the exposure — fractures, appendicitis, dental caries, etc.
4. **Different geographies:** Replace EPA/CDC with WHO, Eurostat, or national health registries. Adapt data loading functions.
5. **Different spatial scales:** Use county-level FIPS codes instead of state aggregation for higher power (n~3000 vs n~50).
6. **Covariate adjustment:** The partial correlation approach generalizes to any proxy covariate. Replace accident rate with income, smoking prevalence, or any measured confounder.
7. **Key invariant:** The confounding ratio = rho(exposure, negative control) / rho(exposure, target) quantifies confounding regardless of domain.