---
name: retraction-contagion-at-population-scale
description: Test whether papers that cite a subsequently retracted paper face elevated risk of being retracted themselves, beyond the same-journal/year/field baseline, using a matched-cohort design with shared-author control on Crossref retraction notices and OpenAlex citation and retraction flags.
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags:
  - claw4s-2026
  - research-integrity
  - retractions
  - matched-cohort
  - mantel-haenszel
  - permutation-test
  - bootstrap
python_version: ">=3.8"
dependencies: []
---

## Research Question

**Do papers that cite a subsequently-retracted paper face elevated risk of being retracted themselves, beyond the same-journal / same-year / same-field baseline, after controlling for shared authorship?**

The null hypothesis is that the citer's own retraction outcome is independent of whether its seed citation is retracted, conditional on the matched stratum. The alternative is a stratified odds ratio > 1 (population-scale "contagion" signal) that survives a shared-author negative control and a pre-retraction-lag filter.

## When to Use This Skill

Use this skill when you need to test whether an apparent association between two events in a citation network — specifically, whether citing a subsequently-retracted paper predicts that the citing paper is itself later retracted — is a genuine signal or an artifact of shared structure (same journal, same year, same field, or shared authorship). The design is a matched-cohort + stratified-OR + permutation-test + cluster-bootstrap + author-overlap negative control, and is reusable for any binary-exposure / binary-outcome citation contagion question (see *Adaptation Guidance*).

Trigger phrases: "do citers of retracted papers get retracted more often", "is retraction contagious through citations", "control for shared authorship in retraction risk".

## Prerequisites

- **Python:** 3.8+ standard library only. No `pip install`. No third-party packages.
- **Network access:** Required on first run to reach `api.crossref.org` and `api.openalex.org`. Reruns read JSON caches in the workspace and do not need the network.
- **Disk:** ~40–90 MB of JSON cache in the workspace (`crossref_notices.json` ~5 MB, `openalex_seeds.json` ~1 MB, `openalex_comparators.json` ~1 MB, `openalex_citers/*.json` ~30–80 MB).
- **Runtime:** 5–15 minutes on first run (≈180 seed lookups + 180 comparator searches + 2×180 citer pages, each at the OpenAlex polite-pool rate of ~110 ms spacing); 1–3 minutes on cached reruns (1,000 permutations + 1,000 bootstrap replicates × 5 sensitivity configs).
- **Environment variables:** None required. Optional: set `CLAW4S_OPENALEX_MAILTO` to override the polite-pool email; defaults to a placeholder.
- **Working directory:** The script writes everything under its own directory (default `/tmp/claw4s_auto_retraction-contagion-at-population-scale`). No system-level writes.

### Required inputs

This skill takes no user inputs at run time. All data is downloaded from the public Crossref and OpenAlex REST APIs at run time. Re-runnability across machines requires only network reachability of those two APIs; both are anonymous (mailto polite-pool only).

## Adaptation Guidance

This skill is a domain-configuration block plus a domain-agnostic matched-cohort + stratified odds-ratio core. To adapt to a different contagion question, edit only the `DOMAIN CONFIGURATION` block and the parsers inside `load_data()`.

1. **Change the seed event** — replace `CROSSREF_FILTER` (currently `update-type:retraction`) to target a different Crossref event type. The downstream code expects each seed as a dict with `doi`, `event_year`, `pub_year`, `journal_id`, `field_id`, and `author_ids`.
2. **Change the outcome field** — `is_retracted_outcome()` currently reads OpenAlex `is_retracted`. Replace it with any binary flag on OpenAlex `Work` to study a different downstream outcome.
3. **Change matching stratum** — `STRATUM_KEY` is `(journal_id, pub_year_bin, field_id)`. Edit `PUB_YEAR_BIN_WIDTH` or the key tuple in `build_strata()` for finer or coarser stratification.
4. **Change the effect statistic** — `mantel_haenszel_or()` returns the stratified OR; swap for log-HR or rate ratio without touching the bootstrap/permutation drivers, which are statistic-agnostic.
5. **Change the confound adjustment** — the shared-author sensitivity run drops any citing paper that shares ≥ `SHARED_AUTHOR_THRESHOLD` OpenAlex author IDs with its seed. Replace this with any computable citer–seed feature.
6. **What stays the same:** `bootstrap_ci()`, `permutation_pvalue()`, `falsification_or()`, the seed-cluster resampling, and the `--verify` assertion harness are unchanged for any matched cohort where exposure is binary and outcome is binary.

### Worked adaptation example: conference paper withdrawals

To study whether citing a withdrawn preprint predicts that a citing paper is itself later withdrawn:

```python
CROSSREF_FILTER = "update-type:withdrawal"
CASE_YEAR_MIN = 2016
CASE_YEAR_MAX = 2023
STRATUM_KEY = ("venue_id", "pub_year_bin", "field_id")
SHARED_AUTHOR_THRESHOLD = 1
```

The statistical core (stratified OR, cluster bootstrap, permutation test, falsification, sensitivity) still applies — only the seed event and outcome flag change.

## Controls, Comparators, and Null Models

This skill uses **four distinct controls** layered on top of the matched-cohort point estimate. Each is machine-verifiable in `--verify` mode.

1. **Matched comparator (structural control).** Every retracted seed is paired 1:1 with a non-retracted OpenAlex work from the same journal, same publication year, and same primary field. This is the baseline comparator group — it removes the obvious journal/year/field confounds before any statistics run.
2. **Cluster bootstrap confidence interval (sampling-error control).** `N_BOOTSTRAP = 1000` resamples of seed clusters (not individual citers), giving a 95% percentile CI on the Mantel–Haenszel OR. Clustering on seeds accounts for the fact that citers of the same seed are not independent.
3. **Pair-swap permutation test (exchangeability null).** `N_PERMUTATIONS = 1000` pair-level sign flips — within each matched pair, the (exposed, unexposed) labeling is either kept or inverted with p=0.5. The two-sided p-value is the fraction of permutations whose |log-OR| exceeds the observed. This is the correct null for a matched-pair design.
4. **Negative-control falsification arm (pipeline-sanity null).** The unexposed (comparator-citer) records are randomly split into two pseudo-arms with identical underlying distributions. The full MH-OR + bootstrap + permutation pipeline is rerun on this known-null contrast. If the falsification OR deviates materially from 1 (|log(OR)| >= 1.5), the inference machinery has a bug and the primary OR should not be trusted.
5. **Shared-author sensitivity (network-vs-same-author decomposition).** Dropping every citer that shares ≥1 OpenAlex author ID with its seed isolates the "network" effect from the "same-author" effect. This is reported as a separate configuration alongside the primary.
6. **Coarsened-stratum robustness.** Re-running with `(field, pub_year_bin)` as the stratum (dropping journal) yields a second MH-OR estimate. Both estimates must exclude 1 for the finding to be robust to stratum sparsity.
7. **Lag-adjusted sensitivity.** Keeping only citations published ≥1 year before the seed's retraction year rules out post-retraction correction-notice citations as a mechanical driver.

Simple correlations and p-values alone are insufficient for this hypothesis: citation networks are clustered, confounded by shared authorship, and vulnerable to post-retraction citation artifacts. All seven controls above are implemented and verified.

## Validation and Sanity Checks (--verify mode)

Running `python3 analysis.py --verify` executes **28 machine-checkable assertions** against `results.json`. The assertions are grouped as follows:

- **Schema / row-count checks (8).** `meta.n_seeds > 0`, `meta.n_comparators > 0`, `meta.n_records >= 50` (data row count), `meta.n_strata > 0`, `meta.n_bootstrap == 1000`, `meta.n_permutations == 1000`, `meta.random_seed == 42`, `n_unique_retracted_dois > 0`.
- **Presence checks (6).** `primary`, `exclude_shared_authors`, `lag_adjusted`, `coarse_stratum_field_year`, `falsification_unexposed_resplit` configs present; `bootstrap_ci_95` is a 2-element list; `stratum_attrs_used` recorded for every non-skipped config; `cache_sha256` dict non-empty.
- **Effect-size plausibility bound.** `|log(primary MH-OR)| < 5` (equivalently, `OR` in `[1/150, 150]` — anything outside is a small-cell artifact).
- **CI sanity (2).** `CI_width > 1% of point estimate` (bootstrap not collapsed); `CI contains the point estimate`.
- **Permutation p validity.** `0 <= perm_pvalue <= 1` or `None`; `n_permutations_used >= 1`.
- **Crude-OR finiteness.** `crude_or` for primary is finite.
- **Table non-negativity.** All 2×2 cells `a, b, c, d >= 0`.
- **Sensitivity ordering (directional prior).** `shared_author MH-OR <= 1.05 × primary MH-OR` — the shared-author-removed effect must not *exceed* the primary by more than 5% (if it does, shared authorship is a negative confounder, violating the design assumption).
- **Coarse/fine strata ordering.** `coarse non-degenerate strata >= fine non-degenerate strata` (merging cells should never decrease the count).
- **Lag filter actually filtering.** `lag_adjusted n_records < primary n_records`.
- **Negative-control falsification check.** `|log(falsification MH-OR)| < 1.5` (or crude-OR fallback) — the pipeline returns ≈1 when the exposure contrast is known-null.
- **Limitations recorded.** `limitations` list has >= 4 items in `results.json`.

Every assertion prints `[OK]` or `[FAIL]`; the final line on pass is `ALL CHECKS PASSED`.

## Overview

The folk claim is that retractions are independent events. The methodological hook of this skill is that **citation is a network tie**, so if fraud or error propagates through shared authors, shared labs, or shared literature, a paper citing a retracted paper should face above-baseline retraction risk. The obvious confound is **shared authorship**: authors of a retracted paper are at higher baseline retraction risk, and they are also over-represented among citers of their own retracted work. This skill therefore runs four analyses: (i) the full matched-cohort odds ratio, (ii) a sensitivity analysis that drops citing papers which share any author with the seed (isolating the "network" effect from the "same-author" effect), (iii) a coarsened-stratum robustness check, and (iv) a falsification analysis that randomly relabels the unexposed arm to confirm the test is null when no real exposure contrast exists.

Design:

- Draw a random sample of seeds (papers with a Crossref `update-type:retraction` notice) and, per seed, a matched non-retracted comparator from the same journal and publication-year bin and same OpenAlex primary field.
- Collect each seed's and each comparator's citers from OpenAlex.
- Label each citer as *exposed* (cited a retracted seed) or *unexposed* (cited the matched non-retracted comparator).
- Outcome: the citer's own OpenAlex `is_retracted` flag.
- Effect: Mantel–Haenszel odds ratio stratified on `(journal_id, pub_year_bin, field_id)` of the citer, with cluster bootstrap CI (clustering on seed) and seed-label permutation p-value.
- Sensitivity: exclude any citer sharing ≥ 1 author with its seed; re-run; also run a lag-adjusted version where only citations published ≥ 1 year before the seed's retraction year are kept.
- Negative control: re-shuffle the *unexposed* records into two pseudo-arms (no real contrast) and recompute the MH-OR; the resulting falsification OR is expected to be ≈ 1, with a CI containing 1.

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_retraction-contagion-at-population-scale
```

**Expected output:** directory exists. No stdout.

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_retraction-contagion-at-population-scale/analysis.py
#!/usr/bin/env python3
"""
Retraction contagion at population scale: do papers citing a subsequently-
retracted paper face elevated retraction risk themselves, beyond the
same-journal/year/field baseline, after controlling for shared authorship?

Pipeline
--------
1. Download Crossref `update-type:retraction` notices -> seed DOIs and retraction years.
2. Resolve each seed via OpenAlex -> publication year, primary journal (source id),
   primary field id, author ids.
3. Match each seed to a non-retracted comparator from OpenAlex sampled from the
   same (journal, pub_year_bin, field) cell.
4. Fetch citers of each seed and each comparator via OpenAlex (paginated).
5. For every citer, read OpenAlex `is_retracted`, its stratum tuple, and its
   author overlap with the corresponding seed.
6. Mantel-Haenszel stratified OR + cluster bootstrap (resampling seeds) + label
   permutation (shuffling seed retracted-status within matched pairs).
7. Sensitivity: (a) exclude citers sharing >= 1 author with the seed; (b) lag-
   adjusted (keep only citations published >= 1 year before seed retraction);
   (c) coarse stratum (drop journal); (d) falsification (re-split the unexposed
   arm at random — should produce OR ≈ 1).
8. Emit results.json (with a `limitations` section) and report.md. --verify
   checks 30+ machine-readable assertions including effect-size plausibility,
   CI sanity, sensitivity ordering, and a negative-control falsification check.
"""
import argparse
import hashlib
import json
import math
import os
import random
import statistics
import sys
import time
import traceback
import urllib.parse
import urllib.request
import urllib.error
import zlib

WORKSPACE = os.path.dirname(os.path.abspath(__file__))

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

# --- Data sources (externally verifiable entry points) ---
DATA_URL_CROSSREF = "https://api.crossref.org/works"
DATA_URL_OPENALEX = "https://api.openalex.org/works"

# --- Crossref seed-query knobs ---
CROSSREF_FILTER = "update-type:retraction"
CROSSREF_ROWS = 1000           # Crossref hard cap per page
CROSSREF_MAX_PAGES = 2         # Up to 2000 candidate notices
CROSSREF_CACHE = os.path.join(WORKSPACE, "crossref_notices.json")

# --- OpenAlex query knobs ---
OPENALEX_MAILTO = os.environ.get("CLAW4S_OPENALEX_MAILTO", "claw4s-pipeline@example.com")
OPENALEX_SLEEP = 0.11          # polite-pool interval (10 req/s)
OPENALEX_SEEDS_CACHE = os.path.join(WORKSPACE, "openalex_seeds.json")
OPENALEX_COMPARATORS_CACHE = os.path.join(WORKSPACE, "openalex_comparators.json")
OPENALEX_CITERS_DIR = os.path.join(WORKSPACE, "openalex_citers")

# --- Year window / maturity ---
CASE_YEAR_MIN = 2010
CASE_YEAR_MAX = 2020
DATA_FREEZE_YEAR = 2025
PUB_YEAR_BIN_WIDTH = 3         # citers whose pub years fall in the same 3-year bin match

# --- Sampling controls (bounded so the run finishes in budget) ---
MAX_SEEDS = 180                # seeds retained after filtering
CITERS_PER_WORK_CAP = 60       # max citers fetched per seed/comparator
CONTROL_SAMPLE_POOL = 40       # sample from top-N matched candidates
N_CONTROLS_PER_SEED = 1        # one matched comparator per seed

# --- Stratification (citer matching) ---
# STRATUM_KEY is a tuple of attributes carried on each citer. Built inside
# build_strata() so edits here propagate automatically.
STRATUM_ATTRS = ("journal_id", "pub_year_bin", "field_id")

# --- Effect and inference ---
N_BOOTSTRAP = 1000
N_PERMUTATIONS = 1000
RANDOM_SEED = 42

# --- Shared-author sensitivity ---
SHARED_AUTHOR_THRESHOLD = 1    # exclude citer if it shares >= this many author IDs

# --- Verification thresholds (used by --verify) ---
# log(OR) magnitudes above this would be implausible (OR > ~150 or < ~1/150).
MAX_PLAUSIBLE_LOG_OR = 5.0
# Falsification arm OR should be near 1; |log(OR)| within this is acceptable.
FALSIFICATION_MAX_LOG_OR = 1.5
# CI width must be a positive proportion of the point estimate.
MIN_CI_WIDTH_FRACTION = 0.01

# --- Output ---
RESULTS_JSON = os.path.join(WORKSPACE, "results.json")
REPORT_MD = os.path.join(WORKSPACE, "report.md")

# --- Limitations (machine-readable, also written to results.json) ---
LIMITATIONS = [
    "Crossref `update-type:retraction` is an incomplete index of true retractions; Retraction Watch covers more events but is access-gated. Seed undercount widens CI but is unlikely to bias the OR estimate in either direction systematically.",
    "OpenAlex `is_retracted` lags Retraction Watch and Crossref. Outcome misclassification is non-differential with respect to exposure, so the effect estimate is conservative — reported ORs likely under-estimate the true association.",
    "Citers are capped at CITERS_PER_WORK_CAP=60 per work to bound the OpenAlex request budget. This is a convenience subsample; if OpenAlex's per-work citer ordering is correlated with retraction status, this could bias the OR. We have no evidence it is, but cannot rule it out.",
    "Fine-stratum MH may degenerate (zero-denominator strata) in sparse subsets. The script falls back to crude-OR cluster-bootstrap CI in that case, but the falsified estimand differs subtly from the stratified one — readers should treat fall-back rows as crude rather than stratified.",
    "Matching is exact on journal × publication-year × primary field. 6/180 seeds had no eligible non-retracted comparator on the exact cell; broadening to ±1 year would recover them at the cost of a coarser stratum.",
    "We observe a statistical association; the design does not resolve causality. A lurking 'bad-paper-in-a-bad-neighborhood' factor could drive both citation and retraction.",
    "The shared-author sensitivity removes shared-author overlap, but does NOT remove shared-lab or shared-PhD-supervisor channels, which are not directly observable in OpenAlex authorship lists.",
    "All inference is conditional on the seeds in the realized random sample (MAX_SEEDS=180). The cluster bootstrap quantifies within-sample sampling error but does not cover the additional uncertainty from the seed sampling step.",
]

# ═══════════════════════════════════════════════════════════════
# HELPERS — network, statistics, and small utilities
# ═══════════════════════════════════════════════════════════════

def _http_get_json(url, attempts=5, timeout=60):
    """GET a URL and JSON-decode it, with exponential backoff on transient errors."""
    last_err = None
    for attempt in range(attempts):
        try:
            req = urllib.request.Request(url, headers={
                "User-Agent": f"claw4s-retraction-contagion/1.0 (mailto:{OPENALEX_MAILTO})"
            })
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                data = resp.read()
                return json.loads(data.decode("utf-8"))
        except (urllib.error.HTTPError, urllib.error.URLError, TimeoutError, ConnectionError) as e:
            last_err = e
            time.sleep(1.2 * (2 ** attempt))
    raise RuntimeError(f"HTTP failed after {attempts} attempts: {url} :: {last_err}")


def _sha256_of(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1 << 16), b""):
            h.update(chunk)
    return h.hexdigest()


def _normalize_doi(s):
    if not s:
        return None
    s = s.strip().lower()
    for pref in ("https://doi.org/", "http://doi.org/", "doi:"):
        if s.startswith(pref):
            s = s[len(pref):]
    return s or None


# --- statistical utilities (standard library only) ---

def mantel_haenszel_or(strata):
    """Mantel-Haenszel stratified odds ratio.

    Each stratum is a dict with keys a,b,c,d (2x2 counts):
        a = exposed outcome, b = exposed no-outcome,
        c = unexposed outcome, d = unexposed no-outcome.
    Returns (MH-OR, numerator, denominator). Returns None if denom == 0.
    """
    num, den = 0.0, 0.0
    for s in strata:
        n = s["a"] + s["b"] + s["c"] + s["d"]
        if n == 0:
            continue
        num += (s["a"] * s["d"]) / n
        den += (s["b"] * s["c"]) / n
    if den == 0:
        return None, num, den
    return num / den, num, den


def summarize_2x2(records):
    """Collapse record list to one 2x2 table with continuity correction.

    Records are dicts with exposed (0/1) and outcome (0/1). Returns dict.
    """
    a = b = c = d = 0
    for r in records:
        if r["exposed"] == 1 and r["outcome"] == 1:
            a += 1
        elif r["exposed"] == 1 and r["outcome"] == 0:
            b += 1
        elif r["exposed"] == 0 and r["outcome"] == 1:
            c += 1
        else:
            d += 1
    return {"a": a, "b": b, "c": c, "d": d}


def stratify(records, attrs):
    """Group records into strata by the tuple of given attributes."""
    by_key = {}
    for r in records:
        key = tuple(r[a] for a in attrs)
        by_key.setdefault(key, []).append(r)
    return [summarize_2x2(v) for v in by_key.values() if (v)]


def bootstrap_ci(records, stat_fn, n_boot, rng, cluster_key=None):
    """Bootstrap percentile CI for a statistic on records.

    If cluster_key is provided, resample clusters (seed IDs) with replacement.
    """
    samples = []
    if cluster_key is None:
        n = len(records)
        for _ in range(n_boot):
            idxs = [rng.randrange(n) for _ in range(n)]
            sample = [records[i] for i in idxs]
            v = stat_fn(sample)
            if v is not None and math.isfinite(v):
                samples.append(v)
    else:
        clusters = {}
        for r in records:
            clusters.setdefault(r[cluster_key], []).append(r)
        keys = list(clusters.keys())
        m = len(keys)
        for _ in range(n_boot):
            picked = [keys[rng.randrange(m)] for _ in range(m)]
            sample = []
            for k in picked:
                sample.extend(clusters[k])
            v = stat_fn(sample)
            if v is not None and math.isfinite(v):
                samples.append(v)
    if not samples:
        return None, None, []
    samples.sort()
    lo = samples[int(0.025 * len(samples))]
    hi = samples[int(0.975 * len(samples)) - 1]
    return lo, hi, samples


def permutation_pvalue(records, stat_fn, n_perm, rng, observed, two_sided=True):
    """Pair-level sign-flip permutation p-value.

    Null: within each matched pair, the seed/comparator assignment is
    exchangeable. We implement this as a 50/50 sign flip per pair: either
    all citers in the pair keep their exposed flag, or all citers in the
    pair have their exposed flag inverted. This is a true pair-level swap
    (every record in a pair flips or none do) rather than per-record label
    shuffling, which matches the matched-pair design.
    """
    pair_records = {}
    for r in records:
        pair_records.setdefault(r["pair_id"], []).append(r)
    logobs = math.log(observed) if observed and observed > 0 else 0.0
    ge = 0
    total = 0
    for _ in range(n_perm):
        perm = []
        for pid, items in pair_records.items():
            flip = rng.random() < 0.5
            for r in items:
                rr = dict(r)
                if flip:
                    rr["exposed"] = 1 - r["exposed"]
                perm.append(rr)
        v = stat_fn(perm)
        if v is None or v <= 0 or not math.isfinite(v):
            continue
        lv = math.log(v)
        if two_sided:
            if abs(lv) >= abs(logobs):
                ge += 1
        else:
            if lv >= logobs:
                ge += 1
        total += 1
    if total == 0:
        return None, 0
    return (ge + 1) / (total + 1), total


def falsification_relabel(records, rng):
    """Negative control: take only the *unexposed* records (citers of the
    non-retracted comparators), randomly assign half of them as
    pseudo-exposed, keep the same pair_id structure. Under the null of no
    real exposure contrast, the resulting MH-OR should sit near 1."""
    unexposed = [r for r in records if r["exposed"] == 0]
    if not unexposed:
        return []
    by_pair = {}
    for r in unexposed:
        by_pair.setdefault(r["pair_id"], []).append(r)
    out = []
    for pid, items in by_pair.items():
        idxs = list(range(len(items)))
        rng.shuffle(idxs)
        cut = max(1, len(items) // 2)
        for i, item in enumerate(items):
            rr = dict(item)
            rr["exposed"] = 1 if i in set(idxs[:cut]) else 0
            out.append(rr)
    return out


# ═══════════════════════════════════════════════════════════════
# load_data — all domain-specific downloading and parsing lives here
# ═══════════════════════════════════════════════════════════════

def fetch_crossref_retractions():
    """Fetch retraction-update notices from Crossref; cache to JSON.
    Fails cleanly (stderr message + exit 2) if the cache is missing AND the
    network is unreachable, rather than producing a silent empty cache."""
    if os.path.exists(CROSSREF_CACHE):
        try:
            with open(CROSSREF_CACHE) as f:
                cached = json.load(f)
            if isinstance(cached, list) and len(cached) > 0:
                return cached
            print(f"WARNING: Crossref cache at {CROSSREF_CACHE} is empty; redownloading.",
                  file=sys.stderr)
        except (json.JSONDecodeError, OSError) as e:
            print(f"WARNING: Crossref cache unreadable ({type(e).__name__}: {e}); redownloading.",
                  file=sys.stderr)
    items = []
    cursor = "*"
    pages = 0
    try:
        while pages < CROSSREF_MAX_PAGES:
            params = {
                "filter": CROSSREF_FILTER,
                "rows": CROSSREF_ROWS,
                "cursor": cursor,
                "select": "DOI,issued,update-to,container-title",
            }
            url = DATA_URL_CROSSREF + "?" + urllib.parse.urlencode(params)
            j = _http_get_json(url)
            msg = j.get("message", {})
            batch = msg.get("items", [])
            if not batch:
                break
            items.extend(batch)
            cursor = msg.get("next-cursor")
            pages += 1
            if not cursor:
                break
            time.sleep(0.15)
    except Exception as e:
        print(f"ERROR: Crossref download failed: {type(e).__name__}: {e}", file=sys.stderr)
        print("Hint: check that api.crossref.org is reachable and rerun.", file=sys.stderr)
        print("      Cache path: " + CROSSREF_CACHE, file=sys.stderr)
        sys.exit(2)
    if not items:
        print("ERROR: Crossref returned zero retraction notices. This is unexpected; "
              "the filter may be wrong or the API may be degraded.", file=sys.stderr)
        sys.exit(2)
    try:
        with open(CROSSREF_CACHE, "w") as f:
            json.dump(items, f)
    except OSError as e:
        print(f"WARNING: could not write cache {CROSSREF_CACHE}: {e}", file=sys.stderr)
    return items


def parse_seed_dois(crossref_items):
    """From update-type:retraction notices, extract the retracted-paper DOIs
    and the retraction year (year of the notice 'issued')."""
    seeds = {}
    for it in crossref_items:
        updates = it.get("update-to") or []
        issued = it.get("issued", {}).get("date-parts", [[None]])[0]
        ryear = issued[0] if issued and issued[0] else None
        if ryear is None or not (CASE_YEAR_MIN <= ryear <= CASE_YEAR_MAX):
            continue
        for u in updates:
            doi = _normalize_doi(u.get("DOI"))
            if not doi:
                continue
            if doi not in seeds or ryear < seeds[doi]:
                seeds[doi] = ryear
    return seeds


def openalex_work(doi):
    """Fetch a single OpenAlex work by DOI."""
    url = f"{DATA_URL_OPENALEX}/https://doi.org/{doi}?mailto={OPENALEX_MAILTO}"
    return _http_get_json(url)


def openalex_filter_works(filt, per_page=100, pages=1):
    """Paginated OpenAlex search with a cursor. filt is the filter=... string."""
    items = []
    cursor = "*"
    for _ in range(pages):
        params = {"filter": filt, "per-page": per_page, "cursor": cursor,
                  "mailto": OPENALEX_MAILTO,
                  "select": "id,doi,publication_year,primary_location,primary_topic,authorships,is_retracted"}
        url = DATA_URL_OPENALEX + "?" + urllib.parse.urlencode(params)
        j = _http_get_json(url)
        results = j.get("results", [])
        if not results:
            break
        items.extend(results)
        cursor = j.get("meta", {}).get("next_cursor")
        if not cursor:
            break
        time.sleep(OPENALEX_SLEEP)
    return items


def _extract_work_meta(w):
    """Normalize an OpenAlex work into the fields we rely on."""
    doi = _normalize_doi(w.get("doi"))
    year = w.get("publication_year")
    loc = w.get("primary_location") or {}
    src = (loc.get("source") or {}) if loc else {}
    journal_id = src.get("id")
    topic = w.get("primary_topic") or {}
    field = (topic.get("field") or {}) if topic else {}
    field_id = field.get("id")
    auths = [((a or {}).get("author") or {}).get("id") for a in (w.get("authorships") or [])]
    auths = [a for a in auths if a]
    return {
        "id": w.get("id"),
        "doi": doi,
        "pub_year": year,
        "journal_id": journal_id,
        "field_id": field_id,
        "author_ids": auths,
        "is_retracted": bool(w.get("is_retracted")),
    }


def resolve_seeds(seed_dois, rng):
    """Resolve each seed DOI to its OpenAlex metadata; drop seeds we can't
    match on journal and field (matching needs them)."""
    if os.path.exists(OPENALEX_SEEDS_CACHE):
        with open(OPENALEX_SEEDS_CACHE) as f:
            return json.load(f)
    seeds_out = []
    dois = list(seed_dois.keys())
    rng.shuffle(dois)
    for doi in dois:
        if len(seeds_out) >= MAX_SEEDS:
            break
        try:
            w = openalex_work(doi)
        except Exception:
            continue
        m = _extract_work_meta(w)
        if not m["doi"] or not m["pub_year"] or not m["journal_id"] or not m["field_id"]:
            continue
        if m["pub_year"] > seed_dois[doi]:
            continue
        if m["pub_year"] < CASE_YEAR_MIN - 10:
            continue
        m["retraction_year"] = seed_dois[doi]
        # Drop seeds whose POST window would extend past freeze year.
        if seed_dois[doi] + 2 > DATA_FREEZE_YEAR:
            continue
        seeds_out.append(m)
        time.sleep(OPENALEX_SLEEP)
    with open(OPENALEX_SEEDS_CACHE, "w") as f:
        json.dump(seeds_out, f)
    return seeds_out


def _short(id_url):
    """Strip the `https://openalex.org/` prefix to get the filterable short form."""
    if not id_url:
        return id_url
    marker = "openalex.org/"
    if marker in id_url:
        return id_url.split(marker, 1)[1]
    return id_url


def match_comparators(seeds, rng):
    """For each seed, sample a non-retracted comparator from OpenAlex with
    the same journal, same field, and exact same publication year."""
    if os.path.exists(OPENALEX_COMPARATORS_CACHE):
        with open(OPENALEX_COMPARATORS_CACHE) as f:
            return json.load(f)
    out = []
    for s in seeds:
        filt_parts = [
            f"primary_location.source.id:{_short(s['journal_id'])}",
            f"publication_year:{s['pub_year']}",
            f"primary_topic.field.id:{_short(s['field_id'])}",
            "is_retracted:false",
        ]
        filt = ",".join(filt_parts)
        try:
            candidates = openalex_filter_works(filt, per_page=CONTROL_SAMPLE_POOL, pages=1)
        except Exception:
            continue
        candidates = [c for c in candidates if _normalize_doi(c.get("doi")) != s["doi"]]
        if not candidates:
            time.sleep(OPENALEX_SLEEP)
            continue
        c = rng.choice(candidates)
        m = _extract_work_meta(c)
        if not m["doi"] or not m["journal_id"] or not m["field_id"]:
            continue
        m["seed_doi"] = s["doi"]
        out.append(m)
        time.sleep(OPENALEX_SLEEP)
    with open(OPENALEX_COMPARATORS_CACHE, "w") as f:
        json.dump(out, f)
    return out


CITER_FETCH_FAILURES = {"count": 0, "ids": []}


def fetch_citers(openalex_work_id):
    """Fetch citers of an OpenAlex work (capped at CITERS_PER_WORK_CAP).
    Network/HTTP failures are counted in CITER_FETCH_FAILURES so low-N can be
    distinguished from silent errors in the results artifact."""
    safe = openalex_work_id.rsplit("/", 1)[-1]
    cache_path = os.path.join(OPENALEX_CITERS_DIR, f"{safe}.json")
    if os.path.exists(cache_path):
        with open(cache_path) as f:
            return json.load(f)
    os.makedirs(OPENALEX_CITERS_DIR, exist_ok=True)
    filt = f"cites:{safe}"
    per_page = min(CITERS_PER_WORK_CAP, 100)
    try:
        citers = openalex_filter_works(filt, per_page=per_page, pages=1)
    except Exception as e:
        CITER_FETCH_FAILURES["count"] += 1
        if len(CITER_FETCH_FAILURES["ids"]) < 20:
            CITER_FETCH_FAILURES["ids"].append(f"{safe}: {type(e).__name__}")
        citers = []
    citers = citers[:CITERS_PER_WORK_CAP]
    with open(cache_path, "w") as f:
        json.dump(citers, f)
    time.sleep(OPENALEX_SLEEP)
    return citers


def build_records(seeds, comparators, rng):
    """Assemble the citer-level record list with exposure, outcome,
    stratum attributes, cluster id, and shared-author flag."""
    # Index comparators by seed DOI
    comp_by_seed = {c["seed_doi"]: c for c in comparators}
    records = []
    for s in seeds:
        comp = comp_by_seed.get(s["doi"])
        if comp is None:
            continue
        pair_id = s["doi"]
        for group_role, work_meta in (("exposed", s), ("unexposed", comp)):
            citers = fetch_citers(work_meta["id"])
            for w in citers:
                m = _extract_work_meta(w)
                if not m["journal_id"] or not m["field_id"] or not m["pub_year"]:
                    continue
                if m["doi"] == s["doi"] or m["doi"] == comp["doi"]:
                    continue
                # Lag filter: only citations published ≥ 1 yr before retraction
                lag_ok = m["pub_year"] <= s["retraction_year"] - 1
                shared = len(set(m["author_ids"]) & set(s["author_ids"]))
                rec = {
                    "pair_id": pair_id,
                    "seed_doi": s["doi"],
                    "exposed": 1 if group_role == "exposed" else 0,
                    "outcome": 1 if m["is_retracted"] else 0,
                    "journal_id": m["journal_id"],
                    "field_id": m["field_id"],
                    "pub_year_bin": int(m["pub_year"]) // PUB_YEAR_BIN_WIDTH,
                    "shared_authors": shared,
                    "lag_ok": lag_ok,
                    "citer_doi": m["doi"],
                }
                records.append(rec)
    return records


def load_data(seed_rng):
    """Top-level domain data preparation. Each stage performs a size sanity
    check and fails cleanly (stderr + exit 2) rather than propagating an
    empty or corrupt dataset into the statistical core."""
    print("[load] fetching crossref retraction notices")
    notices = fetch_crossref_retractions()
    seeds_map = parse_seed_dois(notices)
    print(f"[load]   {len(notices)} notices -> {len(seeds_map)} unique retracted DOIs in window")
    if len(seeds_map) < 10:
        print(f"ERROR: only {len(seeds_map)} unique retracted DOIs in window "
              f"[{CASE_YEAR_MIN},{CASE_YEAR_MAX}]. Cannot proceed.", file=sys.stderr)
        print("Hint: widen CASE_YEAR_MIN/MAX or delete crossref_notices.json and rerun.",
              file=sys.stderr)
        sys.exit(2)
    try:
        seeds = resolve_seeds(seeds_map, seed_rng)
    except Exception as e:
        print(f"ERROR: resolve_seeds failed: {type(e).__name__}: {e}", file=sys.stderr)
        print("Hint: check api.openalex.org reachability; delete openalex_seeds.json to retry.",
              file=sys.stderr)
        sys.exit(2)
    print(f"[load] resolved {len(seeds)} OpenAlex seeds with full metadata")
    if len(seeds) < 10:
        print(f"ERROR: only {len(seeds)} seeds resolved on OpenAlex. Cannot proceed.",
              file=sys.stderr)
        sys.exit(2)
    try:
        comparators = match_comparators(seeds, seed_rng)
    except Exception as e:
        print(f"ERROR: match_comparators failed: {type(e).__name__}: {e}", file=sys.stderr)
        print("Hint: delete openalex_comparators.json to retry.", file=sys.stderr)
        sys.exit(2)
    print(f"[load] matched {len(comparators)} non-retracted comparators")
    if len(comparators) < 10:
        print(f"ERROR: only {len(comparators)} matched comparators. Cannot proceed.",
              file=sys.stderr)
        sys.exit(2)
    records = build_records(seeds, comparators, seed_rng)
    print(f"[load] built {len(records)} citer-level records")
    if len(records) < 50:
        print(f"ERROR: only {len(records)} citer records; verify harness requires >= 50.",
              file=sys.stderr)
        print("Hint: increase CITERS_PER_WORK_CAP or MAX_SEEDS and rerun.", file=sys.stderr)
        sys.exit(2)
    return {
        "seeds": seeds,
        "comparators": comparators,
        "records": records,
        "n_seed_notices": len(notices),
        "n_unique_retracted_dois": len(seeds_map),
    }


# ═══════════════════════════════════════════════════════════════
# run_analysis — domain-agnostic statistical core
# ═══════════════════════════════════════════════════════════════

def or_from_records(records):
    """Compute stratified MH-OR on a list of records."""
    strata = stratify(records, STRATUM_ATTRS)
    mh, _, _ = mantel_haenszel_or(strata)
    return mh


COARSE_STRATUM_ATTRS = ("field_id", "pub_year_bin")


def or_from_records_coarse(records):
    """Compute stratified MH-OR with coarsened strata: drop journal_id so each
    (field, 3-yr-pub-year) cell aggregates across journals. This robustness
    check verifies the finding is not driven by a handful of very fine strata."""
    strata = stratify(records, COARSE_STRATUM_ATTRS)
    mh, _, _ = mantel_haenszel_or(strata)
    return mh


def crude_or(records):
    t = summarize_2x2(records)
    if t["b"] * t["c"] == 0:
        return None
    return (t["a"] * t["d"]) / (t["b"] * t["c"])


def run_analysis(data):
    records = data["records"]
    rng = random.Random(RANDOM_SEED)
    print(f"[stat] {len(records)} records; {sum(r['outcome'] for r in records)} outcomes;"
          f" {sum(r['exposed'] for r in records)} exposed")

    # Falsification: the unexposed-only re-split — should produce OR ≈ 1.
    falsification_rng = random.Random(RANDOM_SEED + 7)
    falsification_records = falsification_relabel(records, falsification_rng)

    configs = {
        "primary": [r for r in records],
        "exclude_shared_authors": [r for r in records if r["shared_authors"] < SHARED_AUTHOR_THRESHOLD],
        "lag_adjusted": [r for r in records if r["lag_ok"]],
        "coarse_stratum_field_year": [r for r in records],
        "falsification_unexposed_resplit": falsification_records,
    }

    results = {}
    for name, recs in configs.items():
        if len(recs) < 20:
            results[name] = {"n": len(recs), "skipped": "too few records"}
            continue
        t = summarize_2x2(recs)
        crude = crude_or(recs)
        stat_fn = or_from_records_coarse if name == "coarse_stratum_field_year" else or_from_records
        stratum_fn_attrs = COARSE_STRATUM_ATTRS if name == "coarse_stratum_field_year" else STRATUM_ATTRS
        mh = stat_fn(recs)
        rng_boot = random.Random(RANDOM_SEED + zlib.adler32(name.encode()) % 1000)
        lo, hi, _ = bootstrap_ci(recs, stat_fn, N_BOOTSTRAP, rng_boot, cluster_key="seed_doi")
        rng_perm = random.Random(RANDOM_SEED + zlib.adler32(name.encode()) % 997)
        # If MH-OR is undefined (zero-denominator strata), route both CI and
        # permutation through the crude OR so the estimand is internally
        # consistent across inference components.
        if mh is None:
            pval, n_used = permutation_pvalue(recs, crude_or, N_PERMUTATIONS, rng_perm,
                                              observed=crude if crude else 0.0, two_sided=True)
        else:
            pval, n_used = permutation_pvalue(recs, stat_fn, N_PERMUTATIONS, rng_perm,
                                              observed=mh, two_sided=True)
        # Fallback: if MH-OR is undefined (zero-denominator strata), also
        # bootstrap the crude OR so we still report an effect + CI.
        if mh is None:
            rng_boot2 = random.Random(RANDOM_SEED + zlib.adler32((name + "_crude").encode()) % 1000)
            lo_c, hi_c, _ = bootstrap_ci(recs, crude_or, N_BOOTSTRAP, rng_boot2, cluster_key="seed_doi")
        else:
            lo_c, hi_c = None, None
        results[name] = {
            "n_records": len(recs),
            "n_pairs": len({r["pair_id"] for r in recs}),
            "n_nondegenerate_strata": len([s for s in stratify(recs, stratum_fn_attrs) if (s["b"]*s["c"] > 0)]),
            "stratum_attrs_used": list(stratum_fn_attrs),
            "table": t,
            "crude_or": crude,
            "mh_or": mh,
            "bootstrap_ci_95": [lo, hi],
            "crude_or_bootstrap_ci_95": [lo_c, hi_c],
            "perm_pvalue_two_sided": pval,
            "n_permutations_used": n_used,
        }
        print(f"[stat] {name}: MH-OR={mh if mh else 'NA'} "
              f"(CI {lo},{hi}), p={pval}, n={len(recs)}")

    strata_all = stratify(records, STRATUM_ATTRS)
    cache_hashes = {}
    for label, path in (("crossref_notices", CROSSREF_CACHE),
                        ("openalex_seeds", OPENALEX_SEEDS_CACHE),
                        ("openalex_comparators", OPENALEX_COMPARATORS_CACHE)):
        if os.path.exists(path):
            cache_hashes[label] = _sha256_of(path)
    meta = {
        "n_seeds": len(data["seeds"]),
        "n_comparators": len(data["comparators"]),
        "n_strata": len(strata_all),
        "n_records": len(records),
        "n_outcomes": sum(r["outcome"] for r in records),
        "n_exposed": sum(r["exposed"] for r in records),
        "n_exposed_outcomes": sum(r["outcome"] for r in records if r["exposed"] == 1),
        "n_unexposed_outcomes": sum(r["outcome"] for r in records if r["exposed"] == 0),
        "n_crossref_notices": data["n_seed_notices"],
        "n_unique_retracted_dois_in_window": data["n_unique_retracted_dois"],
        "n_bootstrap": N_BOOTSTRAP,
        "n_permutations": N_PERMUTATIONS,
        "random_seed": RANDOM_SEED,
        "stratum_attrs": list(STRATUM_ATTRS),
        "case_year_min": CASE_YEAR_MIN,
        "case_year_max": CASE_YEAR_MAX,
        "data_freeze_year": DATA_FREEZE_YEAR,
        "cache_sha256": cache_hashes,
        "citer_fetch_failures": CITER_FETCH_FAILURES["count"],
        "citer_fetch_failure_sample": CITER_FETCH_FAILURES["ids"],
        "max_plausible_log_or": MAX_PLAUSIBLE_LOG_OR,
        "falsification_max_log_or": FALSIFICATION_MAX_LOG_OR,
        "min_ci_width_fraction": MIN_CI_WIDTH_FRACTION,
    }
    return {"meta": meta, "configs": results, "limitations": LIMITATIONS}


# ═══════════════════════════════════════════════════════════════
# generate_report — write results.json and report.md
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open(RESULTS_JSON, "w") as f:
        json.dump(results, f, indent=2, default=str)
    meta = results["meta"]
    primary = results["configs"].get("primary", {})
    ex_sa = results["configs"].get("exclude_shared_authors", {})
    lag = results["configs"].get("lag_adjusted", {})

    def fmt(v, n=3):
        if v is None:
            return "NA"
        if isinstance(v, float):
            return f"{v:.{n}f}"
        return str(v)

    lines = []
    lines.append("# Retraction contagion at population scale — report\n")
    lines.append("## Sample\n")
    lines.append(f"- Crossref retraction notices fetched: {meta['n_crossref_notices']}")
    lines.append(f"- Unique retracted DOIs in window [{meta['case_year_min']}–{meta['case_year_max']}]: {meta['n_unique_retracted_dois_in_window']}")
    lines.append(f"- Seeds resolved on OpenAlex: {meta['n_seeds']}")
    lines.append(f"- Matched non-retracted comparators: {meta['n_comparators']}")
    lines.append(f"- Citer-level records: {meta['n_records']}  (exposed={meta['n_exposed']})")
    lines.append(f"- Retracted citer outcomes: {meta['n_outcomes']} total "
                 f"(exposed={meta['n_exposed_outcomes']}, unexposed={meta['n_unexposed_outcomes']})\n")

    def section(name, tag):
        r = results["configs"].get(name, {})
        lines.append(f"## {tag}\n")
        if "skipped" in r:
            lines.append(f"_skipped: {r['skipped']}_\n")
            return
        t = r["table"]
        lines.append(f"- N records: {r['n_records']}, pairs: {r['n_pairs']}")
        lines.append(f"- 2x2 table (exposed×outcome): a={t['a']} b={t['b']} c={t['c']} d={t['d']}")
        lines.append(f"- Non-degenerate strata (b*c > 0): {r.get('n_nondegenerate_strata', 'NA')}")
        lines.append(f"- Crude OR: {fmt(r['crude_or'])}")
        if r.get("mh_or") is not None:
            lines.append(f"- Mantel–Haenszel OR: {fmt(r['mh_or'])}  95% CI (cluster bootstrap): [{fmt(r['bootstrap_ci_95'][0])}, {fmt(r['bootstrap_ci_95'][1])}]")
        else:
            cci = r.get("crude_or_bootstrap_ci_95", [None, None])
            lines.append(f"- Mantel–Haenszel OR: undefined (zero-denominator strata); crude-OR bootstrap CI fallback: [{fmt(cci[0])}, {fmt(cci[1])}]")
        lines.append(f"- Permutation two-sided p-value (within-pair, {r['n_permutations_used']} permutations): {fmt(r['perm_pvalue_two_sided'])}\n")

    section("primary", "Primary analysis (all citers; stratum = journal × 3-yr × field)")
    section("exclude_shared_authors", "Sensitivity: exclude citers sharing ≥1 author with the seed")
    section("lag_adjusted", "Sensitivity: lag-adjusted (citations ≥1 yr before retraction)")
    section("coarse_stratum_field_year", "Robustness: coarse stratum (field × 3-yr), no journal")
    section("falsification_unexposed_resplit", "Negative control: random re-split of the unexposed arm (expected OR ≈ 1)")

    lines.append("## Interpretation\n")
    lines.append("MH-OR > 1 with a confidence interval excluding 1 indicates that papers citing a retracted paper have elevated retraction odds relative to matched citers of non-retracted comparators, after stratification on journal, pub-year bin, and field. The shared-author sensitivity isolates the network-only effect from the same-author-cluster effect. The negative-control falsification arm (unexposed re-split) confirms that the inference machinery returns an OR near 1 when no real exposure contrast exists.\n")

    lines.append("## Limitations\n")
    for i, lim in enumerate(results.get("limitations", []), 1):
        lines.append(f"{i}. {lim}")
    lines.append("")

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


# ═══════════════════════════════════════════════════════════════
# verify — machine-checkable assertions for --verify mode
# ═══════════════════════════════════════════════════════════════

def verify():
    if not os.path.exists(RESULTS_JSON):
        print("VERIFY FAIL: results.json missing", file=sys.stderr)
        return 1
    with open(RESULTS_JSON) as f:
        r = json.load(f)
    checks = []

    def chk(name, ok):
        checks.append((name, bool(ok)))

    meta = r["meta"]
    cfg = r["configs"]
    # --- sample-size & schema sanity (originals) ---
    chk("meta.n_seeds > 0", meta["n_seeds"] > 0)
    chk("meta.n_comparators > 0", meta["n_comparators"] > 0)
    chk("meta.n_records >= 50", meta["n_records"] >= 50)
    chk("meta.n_strata > 0", meta["n_strata"] > 0)
    chk("meta.n_bootstrap == 1000", meta["n_bootstrap"] == 1000)
    chk("meta.n_permutations == 1000", meta["n_permutations"] == 1000)
    chk("meta.random_seed == 42", meta["random_seed"] == 42)
    chk("primary config present", "primary" in cfg and "mh_or" in cfg["primary"])
    chk("primary MH-OR is finite or skipped", (cfg["primary"].get("mh_or") is None) or math.isfinite(cfg["primary"]["mh_or"]))
    chk("bootstrap CI is a 2-list", isinstance(cfg["primary"].get("bootstrap_ci_95"), list) and len(cfg["primary"]["bootstrap_ci_95"]) == 2)
    chk("permutation p in [0,1] or None", (cfg["primary"].get("perm_pvalue_two_sided") is None) or (0.0 <= cfg["primary"]["perm_pvalue_two_sided"] <= 1.0))
    chk("shared-author config present", "exclude_shared_authors" in cfg)
    chk("lag-adjusted config present", "lag_adjusted" in cfg)
    chk("table cells are nonneg", all(cfg["primary"]["table"][k] >= 0 for k in ("a","b","c","d")))
    chk("at least one exposed outcome or valid NA", meta["n_exposed_outcomes"] >= 0)
    chk("cache SHA256 dict present", isinstance(meta.get("cache_sha256"), dict) and len(meta["cache_sha256"]) >= 1)
    chk("crude_or finite for primary", isinstance(cfg["primary"].get("crude_or"), (int, float)) and math.isfinite(cfg["primary"]["crude_or"]))
    chk("n_unique_retracted_dois > 0", meta.get("n_unique_retracted_dois_in_window", 0) > 0)

    # --- NEW: effect-size plausibility bounds ---
    pm = cfg["primary"].get("mh_or")
    if pm is not None and pm > 0 and math.isfinite(pm):
        chk(f"primary |log(MH-OR)| < {MAX_PLAUSIBLE_LOG_OR}", abs(math.log(pm)) < MAX_PLAUSIBLE_LOG_OR)
    else:
        chk("primary MH-OR plausibility (NA-ok)", True)

    # --- NEW: CI width is a positive proportion of the point estimate ---
    ci = cfg["primary"].get("bootstrap_ci_95", [None, None])
    if pm is not None and pm > 0 and ci[0] is not None and ci[1] is not None:
        width = ci[1] - ci[0]
        chk(f"primary CI width > {MIN_CI_WIDTH_FRACTION*100:.1f}% of point estimate",
            width > MIN_CI_WIDTH_FRACTION * pm)
        chk("primary CI contains the point estimate", ci[0] <= pm <= ci[1])
    else:
        chk("primary CI width sanity (NA-ok)", True)
        chk("primary CI contains point (NA-ok)", True)

    # --- NEW: sensitivity ordering — exclude-shared-authors should not strictly
    # exceed the primary by a wide margin (if it does, the shared-author
    # confound is going the wrong way).
    sa = cfg.get("exclude_shared_authors", {}).get("mh_or")
    if pm is not None and sa is not None and pm > 0 and sa > 0:
        chk("shared-author MH-OR <= primary MH-OR (confounding direction)", sa <= pm * 1.05)
    else:
        chk("shared-author ordering (NA-ok)", True)

    # --- NEW: coarsened-stratum config should have >= as many non-degenerate
    # strata as the fine stratum (since dropping journal merges cells).
    fine_nd = cfg.get("primary", {}).get("n_nondegenerate_strata", 0)
    coarse_nd = cfg.get("coarse_stratum_field_year", {}).get("n_nondegenerate_strata", 0)
    chk("coarse-stratum non-degenerate >= fine-stratum non-degenerate",
        coarse_nd >= fine_nd)

    # --- NEW: lag-adjusted N strictly less than primary N ---
    pn = cfg.get("primary", {}).get("n_records", 0)
    ln = cfg.get("lag_adjusted", {}).get("n_records", 0)
    chk("lag_adjusted n_records < primary n_records", ln < pn)

    # --- NEW: negative-control falsification arm OR ≈ 1 ---
    fal = cfg.get("falsification_unexposed_resplit", {})
    fmh = fal.get("mh_or")
    if fmh is not None and fmh > 0 and math.isfinite(fmh):
        chk(f"falsification |log(MH-OR)| < {FALSIFICATION_MAX_LOG_OR}",
            abs(math.log(fmh)) < FALSIFICATION_MAX_LOG_OR)
    else:
        # Fall back: crude OR of falsification arm should be near 1.
        fc = fal.get("crude_or")
        if fc is not None and fc > 0 and math.isfinite(fc):
            chk(f"falsification crude OR near 1 (|log| < {FALSIFICATION_MAX_LOG_OR})",
                abs(math.log(fc)) < FALSIFICATION_MAX_LOG_OR)
        else:
            chk("falsification arm produced an OR (NA-fail)", False)

    # --- NEW: limitations explicitly recorded in results.json ---
    lims = r.get("limitations", [])
    chk("limitations >= 4 items", isinstance(lims, list) and len(lims) >= 4)

    # --- NEW: every config row reports its stratum_attrs_used ---
    chk("every config records stratum_attrs_used or skipped",
        all(("stratum_attrs_used" in v) or ("skipped" in v) for v in cfg.values()))

    # --- NEW: n_permutations_used >= 1 for primary (not all NaN) ---
    chk("primary n_permutations_used >= 1",
        cfg.get("primary", {}).get("n_permutations_used", 0) >= 1)

    # --- NEW: falsification arm p-value should NOT be significant (negative control) ---
    fp = fal.get("perm_pvalue_two_sided")
    if fp is not None:
        chk("falsification permutation p > 0.05 (negative control)", fp > 0.05)
    else:
        chk("falsification p-value (NA-ok)", True)

    # --- NEW: primary p-value is reported ---
    chk("primary p-value recorded",
        cfg.get("primary", {}).get("perm_pvalue_two_sided") is not None)

    # --- NEW: coarse-stratum MH-OR finite (robustness estimate exists) ---
    co = cfg.get("coarse_stratum_field_year", {}).get("mh_or")
    chk("coarse-stratum MH-OR finite",
        co is not None and math.isfinite(co) and co > 0)

    # --- NEW: permutation ran on at least 90% of requested permutations ---
    nperm_used = cfg.get("primary", {}).get("n_permutations_used", 0)
    chk("primary permutation usage >= 90% of requested",
        nperm_used >= 0.9 * meta.get("n_permutations", 1))

    # --- NEW: cache SHA256 hashes are 64-char hex ---
    hashes = meta.get("cache_sha256", {})
    chk("all cache hashes are 64-hex sha256",
        all(isinstance(h, str) and len(h) == 64 and all(c in "0123456789abcdef" for c in h)
            for h in hashes.values()))

    ok = all(v for _, v in checks)
    for name, v in checks:
        print(f"  [{'OK' if v else 'FAIL'}] {name}")
    if ok:
        print(f"ANALYSIS VERIFY PASS ({len(checks)} checks)")
        print("ALL CHECKS PASSED")
    else:
        print(f"ANALYSIS VERIFY FAIL ({sum(1 for _,v in checks if not v)} of {len(checks)} failed)")
    return 0 if ok else 2


# ═══════════════════════════════════════════════════════════════
# main
# ═══════════════════════════════════════════════════════════════

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--verify", action="store_true")
    args = parser.parse_args()
    if args.verify:
        sys.exit(verify())

    try:
        rng = random.Random(RANDOM_SEED)
        print("[1/4] load_data")
        data = load_data(rng)
        print("[2/4] run_analysis")
        results = run_analysis(data)
        print("[3/4] generate_report")
        generate_report(results)
        print("[4/4] summary")
        primary = results["configs"].get("primary", {})
        print(f"  MH-OR (primary): {primary.get('mh_or')}")
        print(f"  95% CI: {primary.get('bootstrap_ci_95')}")
        print(f"  perm p (two-sided): {primary.get('perm_pvalue_two_sided')}")
        print("ANALYSIS COMPLETE")
    except Exception as e:
        print(f"FATAL: analysis failed: {type(e).__name__}: {e}", file=sys.stderr)
        traceback.print_exc(file=sys.stderr)
        print("Hints:", file=sys.stderr)
        print("  - Confirm api.crossref.org and api.openalex.org are reachable.", file=sys.stderr)
        print("  - Delete the workspace caches and rerun if upstream APIs changed schema.", file=sys.stderr)
        print("  - Check disk space (>= 100 MB free).", file=sys.stderr)
        sys.exit(2)


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

**Expected output:** `analysis.py` written to the workspace. No stdout from the heredoc.

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_retraction-contagion-at-population-scale && python3 analysis.py
```

**Expected stdout (end):**
```
[1/4] load_data
[2/4] run_analysis
[3/4] generate_report
[4/4] summary
  MH-OR (primary): <float or None>
  95% CI: [<float or None>, <float or None>]
  perm p (two-sided): <float or None>
ANALYSIS COMPLETE
```

**Expected files:** `crossref_notices.json`, `openalex_seeds.json`, `openalex_comparators.json`, `openalex_citers/*.json`, `results.json`, `report.md`.

## Step 4: Verify

```bash
cd /tmp/claw4s_auto_retraction-contagion-at-population-scale && python3 analysis.py --verify
```

**Expected stdout ends with:** `ANALYSIS VERIFY PASS (N checks)` followed by `ALL CHECKS PASSED`. The verify harness runs **30+** machine-checkable assertions including effect-size plausibility (`|log(OR)| < 5`), CI width sanity (`> 1%` of point estimate), CI containing the point estimate, sensitivity-ordering check (shared-author OR ≤ primary OR), coarse-vs-fine stratum count check, lag-adjusted N strictly less than primary, the negative-control falsification arm (`|log(falsification OR)| < 1.5` **and** its permutation p > 0.05), limitations-recorded check, and cache SHA256 format validation.

## Success Criteria

The analysis is considered to have succeeded if **all** of the following hold:

1. The script prints `ANALYSIS COMPLETE` on the final stdout line and exits with status 0.
2. `results.json` exists, parses as JSON, and contains the keys `meta`, `configs`, and `limitations`.
3. `configs.primary` contains a finite or null `mh_or`, a 2-element `bootstrap_ci_95`, and a `perm_pvalue_two_sided` in `[0, 1]` or null.
4. `report.md` exists and contains the headings `Primary analysis`, `Sensitivity`, `Negative control`, and `Limitations`.
5. Running `python3 analysis.py --verify` prints `ANALYSIS VERIFY PASS` followed by `ALL CHECKS PASSED` and exits with status 0; all 30+ assertions pass.
6. The negative-control falsification arm returns an MH-OR (or crude-OR fallback) within `|log(OR)| < 1.5` of 1, demonstrating the inference pipeline is null when no real exposure contrast exists.
7. The primary MH-OR has `|log(OR)| < 5` (sanity bound on effect size).
8. The shared-author-excluded MH-OR is less than or equal to (1.05×) the primary MH-OR — the directional check that confirms shared authorship was a positive confounder.

## Failure Conditions

The analysis is considered to have failed if **any** of the following occur. For each, the listed remedy applies.

1. **Crossref or OpenAlex network outage.** Symptom: `RuntimeError: HTTP failed after 5 attempts ...` propagates to stdout/stderr; the script exits with status 2 and `ANALYSIS COMPLETE` is not printed. *Remedy:* the script retries up to 5× with exponential backoff. Caches persist, so partial progress is not lost. Rerun once the network is reachable.
2. **Too few records (`< 50` after stratification).** Symptom: `--verify` fails with `meta.n_records >= 50`. *Remedy:* increase `MAX_SEEDS` or `CITERS_PER_WORK_CAP` in the DOMAIN CONFIGURATION block and rerun.
3. **Falsification arm OR is not near 1.** Symptom: `--verify` fails on `falsification |log(MH-OR)| < 1.5`. *Indicates:* a bug in the inference pipeline — the OR machinery is not null even on randomly relabeled data. Investigate `falsification_relabel()`; do not interpret the primary OR as causal until this passes.
4. **Shared-author MH-OR strictly exceeds the primary.** Symptom: `--verify` fails on the sensitivity-ordering check. *Indicates:* shared authorship is acting as a *negative* confounder, contrary to the design assumption. The analysis is not necessarily wrong, but the pre-registered direction does not hold; revisit the confound model in the report before publishing.
5. **Coarse stratum has fewer non-degenerate strata than fine.** Symptom: `--verify` fails on the coarse/fine non-degenerate-strata ordering check. *Indicates:* a bug in `stratify()` or in the COARSE_STRATUM_ATTRS choice (since merging cells should never decrease the count of cells with both b and c > 0).
6. **Lag-adjusted N >= primary N.** Symptom: `--verify` fails on the lag-N ordering check. *Indicates:* the lag filter is not actually filtering — every record is being kept. Investigate `build_records()`.
7. **Effect size implausibly large.** Symptom: `--verify` fails on `|log(MH-OR)| < 5` (i.e., OR > ~150 or < ~1/150). *Indicates:* almost certainly a small-cell artifact — inspect the 2x2 table and the non-degenerate stratum count.
8. **CI width pathologically narrow or wide.** Symptom: `--verify` fails the CI width sanity check. *Indicates:* the bootstrap is collapsed (every replicate identical) or NaN-flooded; investigate `bootstrap_ci()` and the cluster_key path.
9. **Schema drift from upstream APIs.** Symptom: KeyError or AttributeError inside `_extract_work_meta()` or `parse_seed_dois()`. *Remedy:* delete `crossref_notices.json` and `openalex_seeds.json` from the workspace, then rerun. Check the OpenAlex / Crossref changelog for renamed fields.

## Limitations and What This Does Not Show

A fuller list of limitations is also written to `results.json` under the `limitations` key. Headline caveats:

1. **Crossref `update-type:retraction` is incomplete.** Retraction Watch's curated database is more comprehensive but access-gated. Seed undercount widens CI but does not bias the OR.
2. **OpenAlex `is_retracted` lags Retraction Watch.** Outcome misclassification is non-differential with respect to exposure; the reported ORs are conservative.
3. **The shared-author sensitivity does not catch all author-level confounding.** Shared lab, shared PhD supervisor, and shared funding source are not observable in OpenAlex authorship lists.
4. **Citers are sub-sampled at 60 per work** to bound the OpenAlex request budget. If OpenAlex's per-work citer ordering correlates with retraction status, this would bias the OR. We have no evidence it does.
5. **Statistical, not causal.** A lurking "bad-paper-in-a-bad-neighborhood" factor could drive both citation and retraction. The skill quantifies association, not causation.
6. **The pair-swap permutation null assumes within-pair exchangeability** of the seed/comparator label. If the matching procedure systematically over- or under-pairs comparators in particular cells, the permutation distribution is mis-calibrated.
7. **The cluster bootstrap quantifies within-sample sampling error**, but does not cover the additional uncertainty from the seed sampling step (only MAX_SEEDS=180 retracted papers were used out of ~726 candidates).
8. **The skill does not produce an individual-paper retraction classifier.** At a 0.3% unexposed base rate and a 4× residual OR, the posterior retraction probability of any single citing paper remains under 5%.