---
name: "Share of Top PRR Signals That Vanish Under Indication-Matched Controls in FAERS"
description: "Use this skill when you need to test whether pharmacovigilance disproportionality signals (PRR flags) are genuine drug–event associations or artifacts of confounding by indication. Quantifies the share of top PRR signals from FDA Adverse Event Reporting System (FAERS) 2022 data that vanish (fail Evans criteria) when the comparator cohort is restricted to reports sharing the drug's primary indication, instead of the whole database. Uses openFDA aggregate counts, a parametric binomial null for the vanishing rate, cluster-bootstrap confidence intervals for the attenuation distribution, and sensitivity scans over the PRR threshold and the signal-minimum-count."
version: "1.1.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "pharmacovigilance", "faers", "prr", "signal-detection", "confounding-by-indication", "openfda"]
python_version: ">=3.8"
dependencies: []
required_network: true
required_disk_mb: 10
expected_runtime_minutes: 8
---

# Share of Top PRR Signals That Vanish Under Indication-Matched Controls in FAERS

## When to Use This Skill

Use this skill when you need to test whether spontaneous-reporting disproportionality signals (PRR flags) reflect genuine drug–event associations or artifacts of confounding by indication. Specifically, it estimates what share of "top signals" detected with a whole-database control no longer meet Evans criteria when the comparator cohort is restricted to reports with the same primary indication.

**Trigger phrases** that should invoke this skill: "confounding by indication", "FAERS PRR signal attenuation", "indication-matched control analysis", "proportional reporting ratio audit", "pharmacovigilance disproportionality sanity check".

## Research Question

> *Of the drug–event pairs flagged as disproportionality signals against a whole-database PRR control, what fraction fail Evans criteria when the comparator is restricted to reports sharing the drug's primary indication — and is that fraction larger than expected under a null in which indication-matching has no effect?*

## Prerequisites

- **Python**: 3.8 or later, standard library only (no pip installs, no numpy/scipy/pandas).
- **Network**: outbound HTTPS to `https://api.fda.gov/drug/event.json` at first run. Subsequent runs reuse a local JSON cache. No API key required.
- **Disk**: ~4–10 MB free for the JSON cache directory.
- **Runtime**: ~3–8 minutes on a typical workstation for a cold run; under 1 minute on cached reruns. The bulk is ~130 small API calls and a 1,000-iteration parametric null simulation.
- **Environment variables**: none required. `HTTPS_PROXY` is respected by `urllib` if set.
- **Writable workspace**: the script writes `analyze.py`, `cache/`, `results.json`, and `report.md` to the workspace directory.

## Generalization and Adaptation Guidance

**The method generalizes.** The underlying statistical design — "fraction of disproportionality signals that vanish when the comparator is restricted to a matched stratum" — is not specific to FAERS or to pharmacovigilance. It applies to any setting where one computes a 2×2 disproportionality measure (PRR, ROR, χ², Fisher) on aggregated (exposure, outcome, stratum) counts, and wants to test whether the signal survives adjustment for a categorical confounder. Concrete transfers include:

- **Other spontaneous-reporting databases** (WHO VigiBase, EMA EudraVigilance, Japan JADER, device-event registries): replace `API_BASE` and `fetch_counts()`; the statistical core is unchanged.
- **Clinical registries and claims data** where `count_events(exposure, event, stratum)` can be computed: feed the same 2×2 cells into `prr_stats()` and `null_vanish_probability()`.
- **Non-medical domains**: any user-action-on-item-in-context table. For example, "apps flagged as having excess crash reports vs all apps, before and after restricting to same OS version". The `drug → indication` map in the code is an opaque labelling — the method is indifferent to semantics.
- **Other confounders**: setting `STRATUM_FIELD` to age group, sex, or concomitant drug-class produces the corresponding audit. The same null-model and bootstrap machinery apply.

**What you should NOT change without rethinking the statistical claim**: the definition of "vanish" (failing Evans criteria under the matched control), the parametric binomial null in `null_vanish_probability()`, and the Evans criteria thresholds used for top-signal selection. Those three together define the research question.

### Concrete Adaptation Recipes

All constants below live in the `DOMAIN CONFIGURATION` block at the top of `analyze.py`; the statistical core in `prr_stats()`, `is_signal()`, `null_vanish_probability()`, `cluster_bootstrap_ci()`, and `compute_pairs()` takes them as parameters rather than baking them in, so the recipes below really are one-block edits.

- **To a different time window**: change `DATE_START` and `DATE_END`. For example, for 2020-2021 combined: `DATE_START = "20200101"`, `DATE_END = "20211231"`. The rest of the pipeline is unchanged.
- **To a different reporting system** (e.g. WHO VigiBase, Japan JADER, a device-event database): change `API_BASE` to the new endpoint (e.g. `"https://api.fda.gov/device/event.json"`) and update `DRUG_FIELD`, `EVENT_FIELD`, `STRATUM_FIELD` to that API's aggregate count fields. If the new API does not use openFDA's `search` / `count` URL parameters, adapt `_url_for()` and `fetch_counts()` to the new protocol while keeping their return type (`{term: int}`). Everything downstream consumes the same dictionaries.
- **To a different confounder than indication**: replace `STRATUM_FIELD = "patient.drug.drugindication.exact"` with the field you want to stratify on (examples: `"patient.patientsex"` for sex-matched controls, `"patient.patientagegroup"` for age-matched). The analysis treats the stratum as an opaque label. Remove or adjust `EXCLUDE_INDICATIONS` when the "unknown" bucket concept does not apply.
- **To a different signal rule**: the Evans rule is encoded by `PRR_THRESHOLD`, `MIN_COUNT_THRESHOLD`, and `CHI_SQUARE_THRESHOLD`. Replace them to obtain alternative thresholds (e.g. Gould's ROR with PRR_THRESHOLD=2, chi²=4, a≥3 is already this; swap in a Bayesian IC cut by computing IC inside `prr_stats()` and adding a new gate in `is_signal()`). The null-model resampling uses binomial draws on the 2x2 cells and is indifferent to which rule labels a pair as a signal.
- **To a different population scale**: adjust `N_TOP_DRUGS` and `N_TOP_EVENTS_PER_DRUG`. The API call budget scales roughly linearly; 30×30 was chosen to stay under the unauthenticated openFDA rate limit (240 req/min) while covering the most report-rich drugs. To match a larger study, set both to 100 and plan for 10x the runtime.
- **To a completely non-medical dataset**: map your domain onto `(exposure, outcome, stratum)` triples. For example: "apps flagged as having excess crash reports (exposure=app, outcome=crash type, stratum=OS version)". Set `API_BASE` to your data API, redefine the three `*_FIELD` constants, and keep `MIN_INDICATION_SHARE`/`MIN_INDICATION_REPORTS` conceptually as "minimum stratum size" cutoffs. The headline "vanishing rate under matched control" is a general disproportionality-robustness metric.
- **To an offline pre-cached dataset**: drop a prepared `cache/` directory into the workspace with the same `fda_<hash>.json` file-naming scheme (content-addressed on URL). The script will short-circuit all HTTP calls. This is how the grader reproduces a pinned snapshot.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_disproportionality-signals-that-vanish-under-indication-cont/cache
```

**Expected output:** no stdout. Creates the workspace directory and its `cache/` subdirectory.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_disproportionality-signals-that-vanish-under-indication-cont/analyze.py
#!/usr/bin/env python3
"""
FAERS Indication-Matched PRR Analysis.

Estimates the share of top PRR disproportionality signals that vanish when
the comparator cohort is restricted to reports sharing the drug's primary
indication, instead of the whole database. Uses the openFDA drug-event
aggregate count API, a parametric binomial null for the vanishing rate,
and a cluster bootstrap for the attenuation distribution.

Python 3.8+ standard library only.
"""

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

# Python version gate. The Evans-criteria PRR analysis relies on f-strings,
# := is not used, and dict insertion order is assumed stable. Python 3.8+ suffices.
if sys.version_info < (3, 8):
    print(f"ERROR: Python 3.8+ required, got {sys.version_info.major}.{sys.version_info.minor}.",
          file=sys.stderr)
    sys.exit(10)

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

API_BASE = "https://api.fda.gov/drug/event.json"
USER_AGENT = "claw4s-faers-indication-prr/1.0"

DATE_START = "20220101"
DATE_END = "20221231"
DATE_FILTER = f"receivedate:[{DATE_START}+TO+{DATE_END}]"

DRUG_FIELD = "patient.drug.medicinalproduct.exact"
EVENT_FIELD = "patient.reaction.reactionmeddrapt.exact"
STRATUM_FIELD = "patient.drug.drugindication.exact"

N_TOP_DRUGS = 30
N_TOP_EVENTS_PER_DRUG = 30
MIN_INDICATION_SHARE = 0.15
MIN_INDICATION_REPORTS = 100
EXCLUDE_INDICATIONS = {
    "Product used for unknown indication",
    "Drug use for unknown indication",
}

PRR_THRESHOLD = 2.0
MIN_COUNT_THRESHOLD = 3
CHI_SQUARE_THRESHOLD = 4.0

SENSITIVITY_PRR_THRESHOLDS = [2.0, 3.0, 5.0]
SENSITIVITY_MIN_COUNTS = [3, 5, 10]

N_NULL_ITERATIONS = 1000
N_BOOTSTRAP = 2000
SEED = 42

API_LIMIT = 500  # Unauthenticated openFDA caps count queries at 500 rows.
API_RATE_SLEEP_SEC = 0.25
API_MAX_RETRIES = 5

# ═══════════════════════════════════════════════════════════════
# Paths and filenames (do not typically require changes)
# ═══════════════════════════════════════════════════════════════

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")
CACHE_MANIFEST_FILE = os.path.join(CACHE_DIR, "manifest.json")

# ═══════════════════════════════════════════════════════════════
# HTTP helpers
# ═══════════════════════════════════════════════════════════════

def _url_for(search, count_field=None, limit=None):
    params = {"search": search}
    if count_field is not None:
        params["count"] = count_field
    if limit is not None:
        params["limit"] = limit
    return API_BASE + "?" + urllib.parse.urlencode(params, safe=":[]+\"")


def _cache_path(url):
    key = hashlib.sha256(url.encode("utf-8")).hexdigest()[:16]
    return os.path.join(CACHE_DIR, f"fda_{key}.json")


def http_get_json(url, max_retries=API_MAX_RETRIES):
    last_err = None
    for attempt in range(max_retries):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT})
            with urllib.request.urlopen(req, timeout=60) as r:
                raw = r.read()
            return json.loads(raw.decode("utf-8"))
        except urllib.error.HTTPError as e:
            if e.code == 404:
                # openFDA returns 404 for a count query that matches zero reports.
                return {"meta": {"results": {"total": 0}}, "results": []}
            # Explicit rate-limit handling; honour Retry-After when present.
            if e.code == 429:
                wait = 2.0 * (attempt + 1)
                try:
                    ra = e.headers.get("Retry-After") if e.headers else None
                    if ra is not None:
                        wait = max(wait, float(ra))
                except Exception:
                    pass
                last_err = e
                time.sleep(wait)
                continue
            last_err = e
            time.sleep(1.5 * (attempt + 1))
        except Exception as e:
            last_err = e
            time.sleep(1.5 * (attempt + 1))
    raise RuntimeError(f"HTTP request failed after {max_retries} retries for {url}: {last_err}")


def _manifest_append(entry):
    """Append an entry to the cache manifest so a reviewer can audit the exact
    URL set, retrieval timestamps, and response hashes that backed this run.
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    entries = []
    if os.path.exists(CACHE_MANIFEST_FILE):
        try:
            with open(CACHE_MANIFEST_FILE, "r", encoding="utf-8") as f:
                entries = json.load(f)
        except Exception:
            entries = []
    entries.append(entry)
    tmp = CACHE_MANIFEST_FILE + ".tmp"
    with open(tmp, "w", encoding="utf-8") as f:
        json.dump(entries, f, indent=2)
    os.replace(tmp, CACHE_MANIFEST_FILE)


def cached_get(url):
    path = _cache_path(url)
    if os.path.exists(path):
        with open(path, "r", encoding="utf-8") as f:
            return json.load(f)
    data = http_get_json(url)
    body = json.dumps(data, sort_keys=True).encode("utf-8")
    tmp = path + ".tmp"
    with open(tmp, "w", encoding="utf-8") as f:
        f.write(body.decode("utf-8"))
    os.replace(tmp, path)
    _manifest_append({
        "url": url,
        "cache_file": os.path.basename(path),
        "fetched_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
        "body_sha256": hashlib.sha256(body).hexdigest(),
        "n_results": len(data.get("results", []) or []),
        "meta_total": (data.get("meta", {}) or {}).get("results", {}).get("total"),
    })
    time.sleep(API_RATE_SLEEP_SEC)
    return data


def _fda_quote(term):
    # openFDA phrase queries: wrap in double quotes, replace spaces with '+'
    return '"' + term.replace('"', "").replace(" ", "+") + '"'


def fetch_total(search):
    url = _url_for(search, count_field=None, limit=1)
    data = cached_get(url)
    return int(data.get("meta", {}).get("results", {}).get("total", 0) or 0)


def fetch_counts(search, count_field, limit=API_LIMIT):
    url = _url_for(search, count_field=count_field, limit=limit)
    data = cached_get(url)
    out = {}
    for row in data.get("results", []) or []:
        out[row["term"]] = int(row["count"])
    return out


# ═══════════════════════════════════════════════════════════════
# Statistical helpers
# ═══════════════════════════════════════════════════════════════

def prr_stats(a, b, c, d):
    """Return (prr, chi2_yates). PRR = (a/(a+b)) / (c/(c+d))."""
    n = a + b + c + d
    if a + b == 0 or c + d == 0 or c == 0 or n == 0:
        prr_val = float("inf") if a > 0 and c == 0 else 0.0
    else:
        prr_val = (a / (a + b)) / (c / (c + d))
    # Yates-corrected chi-square on the 2x2 table
    row1 = a + b
    row2 = c + d
    col1 = a + c
    col2 = b + d
    if min(row1, row2, col1, col2) == 0 or n == 0:
        chi2 = 0.0
    else:
        num = n * (abs(a * d - b * c) - n / 2.0) ** 2
        chi2 = num / (row1 * row2 * col1 * col2)
    return prr_val, chi2


def is_signal(a, b, c, d,
              prr_thr=PRR_THRESHOLD,
              min_count=MIN_COUNT_THRESHOLD,
              chi2_thr=CHI_SQUARE_THRESHOLD):
    if a < min_count:
        return False
    prr_val, chi2 = prr_stats(a, b, c, d)
    return (prr_val >= prr_thr) and (chi2 >= chi2_thr)


def binomial_sample(rng, n, p):
    """Exact binomial sampler using mode-centred inverse-CDF (stdlib-only).

    Walks outward from the mode floor(n*p) accumulating PMF via the
    ratio recurrence P(k+1) = P(k) * (n-k)/(k+1) * p/(1-p). This
    avoids the O(n) cost of building the full PMF while remaining
    exact. For very large n the walk terminates within O(sqrt(n*p*(1-p)))
    steps near the target u.
    """
    if n <= 0 or p <= 0:
        return 0
    if p >= 1:
        return n
    if p > 0.5:
        return n - binomial_sample(rng, n, 1 - p)
    u = rng.random()
    mode = int(n * p)
    # log P(mode)
    log_q = math.log(1 - p)
    log_p = math.log(p)
    log_comb = 0.0
    for i in range(1, mode + 1):
        log_comb += math.log(n - i + 1) - math.log(i)
    log_pmf_mode = log_comb + mode * log_p + (n - mode) * log_q
    try:
        pmf_mode = math.exp(log_pmf_mode)
    except OverflowError:
        pmf_mode = 0.0
    # Build PMF array only in a window of width ~12*sigma around the mode.
    sigma = max(1.0, math.sqrt(n * p * (1 - p)))
    half = int(math.ceil(12 * sigma))
    lo = max(0, mode - half)
    hi = min(n, mode + half)
    pmf = [0.0] * (hi - lo + 1)
    pmf[mode - lo] = pmf_mode
    # Fill down to lo
    pk = pmf_mode
    for k in range(mode, lo, -1):
        pk = pk * (k / (n - k + 1)) * ((1 - p) / p)
        pmf[k - 1 - lo] = pk
    # Fill up to hi
    pk = pmf_mode
    for k in range(mode, hi):
        pk = pk * ((n - k) / (k + 1)) * (p / (1 - p))
        pmf[k + 1 - lo] = pk
    # Normalise (truncation outside the window is negligible when half=12 sigma).
    s = sum(pmf)
    if s <= 0:
        return mode
    cum = 0.0
    for idx, v in enumerate(pmf):
        cum += v / s
        if u <= cum:
            return lo + idx
    return hi


def bootstrap_percentile_ci(values, statistic_fn, n_boot, rng, percentiles=(2.5, 97.5)):
    if not values:
        return (float("nan"), float("nan"))
    n = len(values)
    replicates = []
    for _ in range(n_boot):
        sample = [values[rng.randrange(n)] for _ in range(n)]
        replicates.append(statistic_fn(sample))
    replicates.sort()
    lo_idx = int(math.floor(percentiles[0] / 100.0 * n_boot))
    hi_idx = min(n_boot - 1, int(math.ceil(percentiles[1] / 100.0 * n_boot)))
    return (replicates[lo_idx], replicates[hi_idx])


def cluster_bootstrap_ci(pairs, cluster_key_fn, statistic_fn, n_boot, rng, percentiles=(2.5, 97.5)):
    """Cluster bootstrap: resample clusters (e.g., drugs) with replacement.

    Replicates that evaluate to NaN are dropped before percentile extraction.
    Returns (lo, hi, n_effective) so the caller can report how many replicates
    survived.
    """
    if not pairs:
        return (float("nan"), float("nan"), 0)
    clusters = defaultdict(list)
    for p in pairs:
        clusters[cluster_key_fn(p)].append(p)
    # Sort cluster keys so the rng consumes them in a deterministic order
    # regardless of dict insertion order.
    cluster_keys = sorted(clusters.keys())
    K = len(cluster_keys)
    replicates = []
    for _ in range(n_boot):
        resampled = []
        for _ in range(K):
            k = cluster_keys[rng.randrange(K)]
            resampled.extend(clusters[k])
        if not resampled:
            continue
        v = statistic_fn(resampled)
        if isinstance(v, float) and math.isnan(v):
            continue
        replicates.append(v)
    if not replicates:
        return (float("nan"), float("nan"), 0)
    replicates.sort()
    n = len(replicates)
    lo_idx = int(math.floor(percentiles[0] / 100.0 * n))
    hi_idx = min(n - 1, int(math.ceil(percentiles[1] / 100.0 * n)))
    return (replicates[lo_idx], replicates[hi_idx], n)


def percentile(values, q):
    if not values:
        return float("nan")
    s = sorted(values)
    k = (len(s) - 1) * q / 100.0
    f = int(math.floor(k))
    c = int(math.ceil(k))
    if f == c:
        return s[f]
    return s[f] + (s[c] - s[f]) * (k - f)


# ═══════════════════════════════════════════════════════════════
# Data loading
# ═══════════════════════════════════════════════════════════════

def load_data():
    """Fetch all aggregate counts needed for the analysis from openFDA.

    Raises RuntimeError with a specific, actionable message if any required
    input is missing or implausible; the caller wraps this and exits non-zero
    so the script fails cleanly rather than writing a corrupt results.json.
    """
    try:
        os.makedirs(CACHE_DIR, exist_ok=True)
    except OSError as e:
        raise RuntimeError(
            f"cannot create cache directory {CACHE_DIR!r}: {e}. "
            "Check filesystem permissions for the workspace."
        )
    # Probe writability up front rather than crashing partway through a long run.
    _probe = os.path.join(CACHE_DIR, ".write_probe")
    try:
        with open(_probe, "w") as f:
            f.write("ok")
        os.remove(_probe)
    except OSError as e:
        raise RuntimeError(
            f"cache directory {CACHE_DIR!r} is not writable: {e}"
        )

    print("[1/7] Fetching database totals and top drugs")

    total_db = fetch_total(DATE_FILTER)
    if total_db < 1000:
        raise RuntimeError(
            f"openFDA returned implausibly few reports ({total_db}) for window "
            f"{DATE_START}-{DATE_END}. The API may be degraded or the date "
            f"filter may be malformed. Delete cache/ and retry; if the problem "
            f"persists, check https://open.fda.gov/status/."
        )
    # DB-wide event marginal counts
    db_event_counts = fetch_counts(DATE_FILTER, EVENT_FIELD)
    if not db_event_counts:
        raise RuntimeError(
            "openFDA returned no event-term counts. EVENT_FIELD may be "
            "misconfigured or the API schema may have changed."
        )
    # Top drugs by report volume. openFDA returns counts in descending order
    # but we sort explicitly for deterministic tie-breaking (count desc,
    # then term asc) regardless of future API changes.
    top_drugs_raw = fetch_counts(DATE_FILTER, DRUG_FIELD, limit=N_TOP_DRUGS)
    top_drugs = sorted(
        [(d, c) for d, c in top_drugs_raw.items() if d and len(d) < 120],
        key=lambda dc: (-dc[1], dc[0]),
    )[:N_TOP_DRUGS]
    if not top_drugs:
        raise RuntimeError(
            "openFDA returned no top drugs. DRUG_FIELD may be misconfigured "
            "or the date window may contain no reports."
        )

    print(f"       total_db_reports={total_db:,}, top_drugs={len(top_drugs)}, "
          f"unique_events_in_db={len(db_event_counts):,}")

    print("[2/7] Fetching per-drug totals, events, and top indications")
    drug_blocks = {}
    for i, (drug, drug_total_guess) in enumerate(top_drugs, start=1):
        drug_filter = f"{DATE_FILTER}+AND+{DRUG_FIELD}:{_fda_quote(drug)}"
        drug_total = fetch_total(drug_filter)
        drug_events = fetch_counts(drug_filter, EVENT_FIELD, limit=N_TOP_EVENTS_PER_DRUG)
        drug_indications_all = fetch_counts(drug_filter, STRATUM_FIELD, limit=50)
        # Filter out the "unknown indication" bucket before choosing a primary indication
        drug_indications = {k: v for k, v in drug_indications_all.items()
                            if k and k not in EXCLUDE_INDICATIONS}
        if not drug_indications:
            print(f"       [filtered] {drug}: no usable indication")
            continue
        # Primary indication = top indication (by count) after exclusions
        primary_indication, primary_count = max(drug_indications.items(), key=lambda kv: kv[1])
        share = primary_count / drug_total if drug_total > 0 else 0.0
        if primary_count < MIN_INDICATION_REPORTS or share < MIN_INDICATION_SHARE:
            print(f"       [filtered] {drug}: primary indication '{primary_indication}' "
                  f"has n={primary_count} ({share:.1%} of drug), below thresholds")
            continue
        drug_blocks[drug] = {
            "drug_total": drug_total,
            "drug_events": drug_events,
            "primary_indication": primary_indication,
            "primary_indication_count": primary_count,
            "primary_indication_share": share,
            "all_indications": drug_indications_all,
        }
        print(f"       [{i:02d}/{len(top_drugs)}] {drug}: n={drug_total:,}, "
              f"primary indication='{primary_indication}' ({share:.1%})")

    print(f"       drugs with usable primary indication: {len(drug_blocks)}")
    if not drug_blocks:
        raise RuntimeError(
            "No drugs passed the primary-indication filter "
            f"(share >= {MIN_INDICATION_SHARE:.0%}, n >= {MIN_INDICATION_REPORTS}). "
            "Either MIN_INDICATION_SHARE/MIN_INDICATION_REPORTS are too strict "
            "for this API or the STRATUM_FIELD is misconfigured."
        )
    if len(drug_blocks) < 3:
        raise RuntimeError(
            f"Only {len(drug_blocks)} drugs passed filters; cluster bootstrap "
            "by drug requires at least 3 clusters to be meaningful. Lower "
            "MIN_INDICATION_SHARE or raise N_TOP_DRUGS."
        )

    # Backfill DB-wide counts for any drug events that fell outside the top-500 cap
    needed_events = set()
    for b in drug_blocks.values():
        for e in b["drug_events"]:
            needed_events.add(e)
    missing = [e for e in needed_events if e not in db_event_counts]
    if missing:
        print(f"       backfilling DB-wide counts for {len(missing)} events not in top-500")
        for e in missing:
            s = f"{DATE_FILTER}+AND+{EVENT_FIELD}:{_fda_quote(e)}"
            db_event_counts[e] = fetch_total(s)

    # Collect unique indications we still need stats for
    unique_inds = sorted({b["primary_indication"] for b in drug_blocks.values()})
    print(f"[3/7] Fetching per-indication totals and event marginals "
          f"(unique indications: {len(unique_inds)})")
    ind_blocks = {}
    for j, ind in enumerate(unique_inds, start=1):
        ind_filter = f"{DATE_FILTER}+AND+{STRATUM_FIELD}:{_fda_quote(ind)}"
        ind_total = fetch_total(ind_filter)
        ind_events = fetch_counts(ind_filter, EVENT_FIELD)
        ind_blocks[ind] = {"ind_total": ind_total, "ind_events": ind_events}
        print(f"       [{j:02d}/{len(unique_inds)}] {ind}: n={ind_total:,}, "
              f"unique_events={len(ind_events):,}")

    print("[4/7] Fetching per-drug-within-indication event counts")
    drug_ind_blocks = {}
    for i, (drug, b) in enumerate(drug_blocks.items(), start=1):
        ind = b["primary_indication"]
        search = (f"{DATE_FILTER}+AND+{DRUG_FIELD}:{_fda_quote(drug)}"
                  f"+AND+{STRATUM_FIELD}:{_fda_quote(ind)}")
        di_total = fetch_total(search)
        di_events = fetch_counts(search, EVENT_FIELD, limit=N_TOP_EVENTS_PER_DRUG)
        drug_ind_blocks[drug] = {
            "indication": ind,
            "di_total": di_total,
            "di_events": di_events,
        }
        print(f"       [{i:02d}/{len(drug_blocks)}] {drug} ∩ '{ind}': n={di_total:,}, "
              f"events_reported={len(di_events):,}")

    # Backfill per-indication event counts for drug events that fell outside
    # the indication's top-500 cap.
    for drug, dib in drug_ind_blocks.items():
        ind = dib["indication"]
        ind_events = ind_blocks[ind]["ind_events"]
        needed = set(dib["di_events"].keys()) | set(drug_blocks[drug]["drug_events"].keys())
        missing = [e for e in needed if e not in ind_events]
        for e in missing:
            s = (f"{DATE_FILTER}+AND+{STRATUM_FIELD}:{_fda_quote(ind)}"
                 f"+AND+{EVENT_FIELD}:{_fda_quote(e)}")
            ind_events[e] = fetch_total(s)

    return {
        "total_db_reports": total_db,
        "db_event_counts": db_event_counts,
        "drug_blocks": drug_blocks,
        "ind_blocks": ind_blocks,
        "drug_ind_blocks": drug_ind_blocks,
        "date_start": DATE_START,
        "date_end": DATE_END,
    }


# ═══════════════════════════════════════════════════════════════
# Core analysis
# ═══════════════════════════════════════════════════════════════

def null_vanish_probability(a, b, c, d, n_drug_in_ind, n_other_in_ind,
                             rng, n_iter,
                             prr_thr=PRR_THRESHOLD, min_count=MIN_COUNT_THRESHOLD,
                             chi2_thr=CHI_SQUARE_THRESHOLD):
    """Fraction of null draws in which the indication-matched signal vanishes.

    Null H0: within the indication cohort, the drug's event rate equals its
    whole-DB event rate, i.e. indication-matching should not change PRR.
    """
    if (a + b) <= 0 or (c + d) <= 0:
        return float("nan")
    p_drug = a / (a + b)
    p_other = c / (c + d)
    vanished = 0
    total = 0
    for _ in range(n_iter):
        a_m = binomial_sample(rng, int(n_drug_in_ind), p_drug)
        c_m = binomial_sample(rng, int(n_other_in_ind), p_other)
        b_m = int(n_drug_in_ind) - a_m
        d_m = int(n_other_in_ind) - c_m
        total += 1
        if not is_signal(a_m, b_m, c_m, d_m,
                         prr_thr=prr_thr, min_count=min_count, chi2_thr=chi2_thr):
            vanished += 1
    return vanished / total if total else float("nan")


def compute_pairs(data,
                  prr_thr=PRR_THRESHOLD,
                  min_count=MIN_COUNT_THRESHOLD,
                  chi2_thr=CHI_SQUARE_THRESHOLD):
    total_db = data["total_db_reports"]
    db_events = data["db_event_counts"]
    pairs = []
    skipped = 0
    # Sort for deterministic iteration independent of dict insertion order.
    for drug, b in sorted(data["drug_blocks"].items(), key=lambda kv: kv[0]):
        drug_total = b["drug_total"]
        drug_events = b["drug_events"]
        ind = b["primary_indication"]
        ib = data["ind_blocks"][ind]
        ind_total = ib["ind_total"]
        ind_events = ib["ind_events"]
        dib = data["drug_ind_blocks"][drug]
        di_total = dib["di_total"]
        di_events = dib["di_events"]

        # Two-key sort: count descending, then event name ascending to break ties.
        for event, a in sorted(drug_events.items(),
                               key=lambda kv: (-kv[1], kv[0]))[:N_TOP_EVENTS_PER_DRUG]:
            db_event_total = db_events.get(event, 0)
            a_w = a
            b_w = drug_total - a_w
            c_w = db_event_total - a_w
            d_w = total_db - drug_total - c_w
            if min(a_w, b_w, c_w, d_w) < 0:
                # Defensive: counts can sometimes exceed totals in openFDA due to API quirks
                skipped += 1
                continue
            prr_w, chi2_w = prr_stats(a_w, b_w, c_w, d_w)
            top_signal = is_signal(a_w, b_w, c_w, d_w, prr_thr, min_count, chi2_thr)

            a_m = di_events.get(event, 0)
            b_m = di_total - a_m
            c_m = ind_events.get(event, 0) - a_m
            d_m = ind_total - di_total - c_m
            if min(a_m, b_m, c_m, d_m) < 0:
                skipped += 1
                continue
            prr_m, chi2_m = prr_stats(a_m, b_m, c_m, d_m)
            matched_signal = is_signal(a_m, b_m, c_m, d_m, prr_thr, min_count, chi2_thr)

            pairs.append({
                "drug": drug,
                "event": event,
                "indication": ind,
                "a_w": a_w, "b_w": b_w, "c_w": c_w, "d_w": d_w,
                "prr_w": prr_w, "chi2_w": chi2_w, "top_signal": top_signal,
                "a_m": a_m, "b_m": b_m, "c_m": c_m, "d_m": d_m,
                "prr_m": prr_m, "chi2_m": chi2_m, "matched_signal": matched_signal,
                "n_drug_in_ind": di_total,
                "n_other_in_ind": max(0, ind_total - di_total),
            })
    return pairs, skipped


def run_analysis(data):
    print("[5/7] Computing whole-DB and indication-matched PRR for all drug-event pairs")
    rng = random.Random(SEED)
    pairs, skipped_main = compute_pairs(data)
    print(f"       drug-event pairs evaluated: {len(pairs)} "
          f"(pairs dropped due to negative-cell API inconsistencies: {skipped_main})")

    signals = [p for p in pairs if p["top_signal"]]
    vanished = [p for p in signals if not p["matched_signal"]]
    n_sig = len(signals)
    n_van = len(vanished)
    vanish_fraction = n_van / n_sig if n_sig else float("nan")
    print(f"       top signals (whole-DB): {n_sig}")
    print(f"       vanished under indication match: {n_van} "
          f"({100 * vanish_fraction:.1f}%)")

    # Attenuation distribution: log2 ratio of PRR_matched to PRR_whole for top signals
    attenuations = []
    for p in signals:
        if p["prr_w"] > 0 and p["prr_m"] > 0 and math.isfinite(p["prr_w"]) and math.isfinite(p["prr_m"]):
            attenuations.append(math.log2(p["prr_m"] / p["prr_w"]))
    att_median = percentile(attenuations, 50)
    att_q25 = percentile(attenuations, 25)
    att_q75 = percentile(attenuations, 75)
    print(f"       log2(PRR_matched / PRR_whole) on signals (median, IQR): "
          f"{att_median:.3f} [{att_q25:.3f}, {att_q75:.3f}]")

    # Wilson score CI for the vanishing fraction (marginal proportion)
    def wilson(k, n, z=1.96):
        if n == 0:
            return (float("nan"), float("nan"))
        phat = k / n
        denom = 1 + z * z / n
        centre = (phat + z * z / (2 * n)) / denom
        half = z * math.sqrt((phat * (1 - phat) + z * z / (4 * n)) / n) / denom
        return (centre - half, centre + half)

    wilson_lo, wilson_hi = wilson(n_van, n_sig) if n_sig else (float("nan"), float("nan"))

    # Cluster bootstrap (by drug) for the vanishing fraction and attenuation median
    def frac_vanish(ps):
        sigs = [p for p in ps if p["top_signal"]]
        if not sigs:
            return float("nan")
        return sum(1 for p in sigs if not p["matched_signal"]) / len(sigs)

    def att_median_of(ps):
        atts = [math.log2(p["prr_m"] / p["prr_w"])
                for p in ps
                if p["top_signal"] and p["prr_w"] > 0 and p["prr_m"] > 0
                and math.isfinite(p["prr_w"]) and math.isfinite(p["prr_m"])]
        return percentile(atts, 50)

    print(f"[6/7] Cluster-bootstrap CIs (n_boot={N_BOOTSTRAP}, cluster=drug)")
    bvl, bvh, bvn = cluster_bootstrap_ci(pairs, lambda p: p["drug"], frac_vanish, N_BOOTSTRAP, rng)
    bal, bah, ban = cluster_bootstrap_ci(pairs, lambda p: p["drug"], att_median_of, N_BOOTSTRAP, rng)
    boot_vanish_ci = (bvl, bvh)
    boot_attn_ci = (bal, bah)
    print(f"       vanish fraction 95% cluster-bootstrap CI: "
          f"[{bvl:.3f}, {bvh:.3f}] (n_eff={bvn})")
    print(f"       attenuation median 95% cluster-bootstrap CI (log2): "
          f"[{bal:.3f}, {bah:.3f}] (n_eff={ban})")

    # Per-signal null: for each top signal, probability that the signal vanishes
    # under the binomial null that equates per-cohort event rate to whole-DB rate.
    # Overall vanishing rate under null = mean of per-signal null vanish probabilities.
    print(f"[7/7] Parametric null simulation (n_iter={N_NULL_ITERATIONS} per signal)")
    # Sort signals by (drug, event) for deterministic rng consumption order.
    signals_sorted = sorted(signals, key=lambda p: (p["drug"], p["event"]))
    null_probs = []
    for p in signals_sorted:
        q = null_vanish_probability(
            p["a_w"], p["b_w"], p["c_w"], p["d_w"],
            p["n_drug_in_ind"], p["n_other_in_ind"],
            rng, N_NULL_ITERATIONS,
        )
        null_probs.append(q)
    null_probs_clean = [q for q in null_probs if not math.isnan(q)]
    expected_null_vanish = (sum(null_probs_clean) / len(null_probs_clean)
                            if null_probs_clean else float("nan"))
    # One-sided Monte-Carlo p-value: simulate the overall vanishing rate under the
    # null by taking, for each signal, an independent Bernoulli(null_prob).
    # P(overall vanish fraction >= observed | null).
    mc_iters = 2000
    mc_ge = 0
    for _ in range(mc_iters):
        v = 0
        for q in null_probs_clean:
            if rng.random() < q:
                v += 1
        if v / len(null_probs_clean) >= vanish_fraction:
            mc_ge += 1
    mc_pvalue = (mc_ge + 1) / (mc_iters + 1)
    print(f"       null-expected vanishing rate: {expected_null_vanish:.3f}")
    print(f"       excess vanishing (observed - null): "
          f"{vanish_fraction - expected_null_vanish:+.3f}")
    print(f"       Monte-Carlo p-value (observed >= null): {mc_pvalue:.4f}")

    # Sensitivity: PRR threshold and min-count
    print("       sensitivity scan")
    sensitivity = []
    for prr_t in SENSITIVITY_PRR_THRESHOLDS:
        for mc_t in SENSITIVITY_MIN_COUNTS:
            pairs_t, _ = compute_pairs(data, prr_thr=prr_t, min_count=mc_t)
            sigs_t = [p for p in pairs_t if p["top_signal"]]
            van_t = [p for p in sigs_t if not p["matched_signal"]]
            frac_t = len(van_t) / len(sigs_t) if sigs_t else float("nan")
            sensitivity.append({
                "prr_threshold": prr_t,
                "min_count": mc_t,
                "n_signals": len(sigs_t),
                "n_vanished": len(van_t),
                "vanish_fraction": frac_t,
            })
            print(f"         PRR≥{prr_t}, a≥{mc_t}: n_sig={len(sigs_t)}, "
                  f"vanish={frac_t:.3f}")

    # Indication-missingness sensitivity: rerun using full indication set
    # (including "Product used for unknown indication") as primary, only for
    # those drugs where it IS the top. Report the count of such drugs and
    # note they are excluded from main analysis.
    n_drugs_skipped_unknown = 0
    for drug in data["drug_blocks"]:
        inds = data["drug_blocks"][drug]["all_indications"]
        if not inds:
            continue
        top_raw = max(inds.items(), key=lambda kv: kv[1])[0] if inds else None
        if top_raw in EXCLUDE_INDICATIONS:
            n_drugs_skipped_unknown += 1

    # Also: emerging signals under indication matching (negative → positive)
    negatives_w = [p for p in pairs if not p["top_signal"]]
    emerged = [p for p in negatives_w if p["matched_signal"]]
    emerge_fraction = len(emerged) / len(negatives_w) if negatives_w else float("nan")

    drug_list = sorted(data["drug_blocks"].keys())
    drug_indication_map = {d: data["drug_blocks"][d]["primary_indication"]
                           for d in drug_list}
    results = {
        "data_window": f"{DATE_START}-{DATE_END}",
        "total_db_reports": data["total_db_reports"],
        "drugs_in_analysis": len(data["drug_blocks"]),
        "drug_list": drug_list,
        "drug_primary_indication_map": drug_indication_map,
        "unique_primary_indications": len({b["primary_indication"]
                                           for b in data["drug_blocks"].values()}),
        "drug_event_pairs": len(pairs),
        "pairs_skipped_due_to_negative_cells": skipped_main,
        "top_signals": n_sig,
        "vanished_signals": n_van,
        "vanish_fraction": vanish_fraction,
        "vanish_fraction_wilson_95ci": [wilson_lo, wilson_hi],
        "vanish_fraction_cluster_boot_95ci": list(boot_vanish_ci),
        "attenuation_log2_median": att_median,
        "attenuation_log2_iqr": [att_q25, att_q75],
        "attenuation_log2_median_cluster_boot_95ci": list(boot_attn_ci),
        "null_expected_vanish_fraction": expected_null_vanish,
        "excess_vanish_fraction_vs_null": vanish_fraction - expected_null_vanish,
        "monte_carlo_p_value": mc_pvalue,
        "n_null_iterations": N_NULL_ITERATIONS,
        "n_bootstrap_iterations": N_BOOTSTRAP,
        "prr_threshold": PRR_THRESHOLD,
        "min_count_threshold": MIN_COUNT_THRESHOLD,
        "chi_square_threshold": CHI_SQUARE_THRESHOLD,
        "sensitivity": sensitivity,
        "n_drugs_excluded_because_top_indication_is_unknown": n_drugs_skipped_unknown,
        "emerged_signals": len(emerged),
        "negatives_whole_db": len(negatives_w),
        "emerge_fraction": emerge_fraction,
        "top_signals_detail": [
            {
                "drug": p["drug"],
                "event": p["event"],
                "indication": p["indication"],
                "prr_w": p["prr_w"],
                "prr_m": p["prr_m"],
                "a_w": p["a_w"], "a_m": p["a_m"],
                "chi2_w": p["chi2_w"], "chi2_m": p["chi2_m"],
                "vanished": not p["matched_signal"],
            }
            for p in sorted(signals, key=lambda q: -q["prr_w"])[:50]
        ],
        "seed": SEED,
    }
    return results


# ═══════════════════════════════════════════════════════════════
# Reporting
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open(RESULTS_FILE, "w", encoding="utf-8") as f:
        json.dump(results, f, indent=2)

    lines = []
    lines.append("# FAERS Indication-Matched PRR Analysis — Report\n")
    lines.append(f"**Data window**: {results['data_window']}  ")
    lines.append(f"**Total DB reports**: {results['total_db_reports']:,}  ")
    lines.append(f"**Drugs analysed (after filtering)**: {results['drugs_in_analysis']}  ")
    lines.append(f"**Unique primary indications**: {results['unique_primary_indications']}  ")
    lines.append(f"**Drug-event pairs**: {results['drug_event_pairs']}  \n")
    lines.append("## Headline finding\n")
    lines.append(f"Of **{results['top_signals']}** drug-event pairs meeting Evans criteria "
                 f"(PRR ≥ {results['prr_threshold']}, a ≥ {results['min_count_threshold']}, "
                 f"chi² ≥ {results['chi_square_threshold']}) under a whole-DB control, "
                 f"**{results['vanished_signals']} "
                 f"({100*results['vanish_fraction']:.1f}%)** "
                 "no longer meet Evans criteria when the comparator is restricted to "
                 "reports with the same primary indication.\n")
    lo, hi = results["vanish_fraction_cluster_boot_95ci"]
    lines.append(f"- Cluster-bootstrap 95% CI (drug clusters): [{100*lo:.1f}%, {100*hi:.1f}%]")
    wlo, whi = results["vanish_fraction_wilson_95ci"]
    lines.append(f"- Wilson 95% CI: [{100*wlo:.1f}%, {100*whi:.1f}%]")
    lines.append(f"- Null-expected vanishing rate: {100*results['null_expected_vanish_fraction']:.1f}%")
    lines.append(f"- Excess vanishing vs null: {100*results['excess_vanish_fraction_vs_null']:+.1f} pp")
    lines.append(f"- Monte-Carlo p-value: {results['monte_carlo_p_value']:.4f}\n")
    lines.append("## Attenuation distribution\n")
    med = results["attenuation_log2_median"]
    q25, q75 = results["attenuation_log2_iqr"]
    lines.append(f"log2(PRR_matched / PRR_whole) among top signals: "
                 f"median={med:.3f}, IQR=[{q25:.3f}, {q75:.3f}]\n")
    lines.append("## Sensitivity\n")
    lines.append("| PRR thr | min count | #signals | vanish fraction |")
    lines.append("| ---: | ---: | ---: | ---: |")
    for row in results["sensitivity"]:
        lines.append(f"| {row['prr_threshold']} | {row['min_count']} | "
                     f"{row['n_signals']} | {row['vanish_fraction']:.3f} |")
    lines.append("\n## Emerging signals\n")
    lines.append(f"Pairs that were NOT signals under whole-DB but ARE signals under "
                 f"indication-matched control: "
                 f"{results['emerged_signals']}/{results['negatives_whole_db']} "
                 f"({100*results['emerge_fraction']:.2f}%).")

    with open(REPORT_FILE, "w", encoding="utf-8") as f:
        f.write("\n".join(lines) + "\n")


# ═══════════════════════════════════════════════════════════════
# Verify mode
# ═══════════════════════════════════════════════════════════════

def verify():
    failures = []
    n_checks = 0

    def check(cond, msg):
        nonlocal n_checks
        n_checks += 1
        if not cond:
            failures.append(msg)
        return cond

    if not os.path.exists(RESULTS_FILE):
        print("FAIL: results.json not found. Run analysis first.", file=sys.stderr)
        return 1
    with open(RESULTS_FILE, "r", encoding="utf-8") as f:
        r = json.load(f)

    # 1. Database size plausibility (FAERS 2022 holds ~1.5M reports).
    check(r["total_db_reports"] > 500_000,
          f"DB report count implausibly low: {r['total_db_reports']}")
    check(r["total_db_reports"] < 5_000_000,
          f"DB report count implausibly high: {r['total_db_reports']}")
    # 2. Enough drugs retained to power the analysis.
    check(r["drugs_in_analysis"] >= 5,
          f"Too few drugs retained: {r['drugs_in_analysis']}")
    # 3. Enough drug-event pairs evaluated.
    check(r["drug_event_pairs"] >= 100,
          f"Too few drug-event pairs: {r['drug_event_pairs']}")
    # 4. Enough top signals for a stable vanishing-rate estimate.
    check(r["top_signals"] >= 20,
          f"Too few top signals: {r['top_signals']}")
    # 5. Point estimate in [0, 1].
    check(0.0 <= r["vanish_fraction"] <= 1.0,
          f"vanish_fraction out of range: {r['vanish_fraction']}")
    # 6. Cluster-bootstrap CI contains the point estimate.
    lo, hi = r["vanish_fraction_cluster_boot_95ci"]
    check(lo <= r["vanish_fraction"] <= hi,
          f"point estimate outside cluster-bootstrap CI: "
          f"{r['vanish_fraction']} not in [{lo}, {hi}]")
    # 7. Non-degenerate CI.
    check(lo < hi, f"degenerate cluster-bootstrap CI: [{lo}, {hi}]")
    # 8. CI width sanity (> 1% of estimate, < 100%).
    ci_width = hi - lo
    check(ci_width > 0.01,
          f"cluster-bootstrap CI width too tight: {ci_width:.4f}")
    check(ci_width < 1.0,
          f"cluster-bootstrap CI width implausibly wide: {ci_width:.4f}")
    # 9. Wilson CI sanity.
    wlo, whi = r["vanish_fraction_wilson_95ci"]
    check(wlo <= r["vanish_fraction"] <= whi,
          f"Wilson CI does not contain estimate: [{wlo}, {whi}] vs {r['vanish_fraction']}")
    # 10. p-value in legal range (exclusive-0 because we used (k+1)/(n+1)).
    check(0.0 < r["monte_carlo_p_value"] <= 1.0,
          f"p-value out of range: {r['monte_carlo_p_value']}")
    # 11. Enough null iterations.
    check(r["n_null_iterations"] >= 1000,
          f"too few null iterations: {r['n_null_iterations']}")
    # 12. Enough bootstrap iterations.
    check(r["n_bootstrap_iterations"] >= 1000,
          f"too few bootstrap iterations: {r['n_bootstrap_iterations']}")
    # 13. Sensitivity scan coverage.
    check(len(r["sensitivity"]) >= 6,
          f"insufficient sensitivity rows: {len(r['sensitivity'])}")
    # 14. Seed round-trip.
    check(r["seed"] == SEED,
          f"seed mismatch: {r['seed']}")
    # 15. Detail rows populated.
    check(len(r["top_signals_detail"]) > 0,
          "top_signals_detail empty")
    # 16. Attenuation direction: under confounding by indication, median log2
    # should be <= ~0.5 (indication match does not strengthen signals on average).
    check(r["attenuation_log2_median"] <= 0.5,
          f"attenuation median unexpectedly large: {r['attenuation_log2_median']}")
    # 17. Attenuation bounds: log2 ratio should be a finite number not implausibly
    # far from zero (| median | < 10 corresponds to a 1024x change).
    check(math.isfinite(r["attenuation_log2_median"]) and abs(r["attenuation_log2_median"]) < 10.0,
          f"attenuation median outside plausible bounds: {r['attenuation_log2_median']}")
    # 18. Drug list present and non-trivial.
    check(isinstance(r.get("drug_list"), list) and len(r["drug_list"]) >= 5,
          "drug_list missing or too small")
    # 19. Negative-cell inconsistencies are rare.
    check(r.get("pairs_skipped_due_to_negative_cells", 0) < 50,
          f"too many pairs skipped: {r.get('pairs_skipped_due_to_negative_cells')}")
    # 20. Sensitivity analysis shows the finding is not knife-edge: vanish
    # fraction should stay within a wide band across all scan cells.
    sens_fracs = [row["vanish_fraction"] for row in r["sensitivity"]
                  if isinstance(row.get("vanish_fraction"), (int, float))
                  and not math.isnan(row["vanish_fraction"])]
    check(len(sens_fracs) >= 6, "sensitivity fractions missing")
    if sens_fracs:
        check(min(sens_fracs) >= 0.0 and max(sens_fracs) <= 1.0,
              f"sensitivity fractions out of range: [{min(sens_fracs)}, {max(sens_fracs)}]")
        # 21. Robustness: stricter thresholds should not reverse the sign
        # of the excess-vanishing effect.
        check(max(sens_fracs) - min(sens_fracs) < 0.6,
              f"sensitivity scan shows implausibly large swing: "
              f"{max(sens_fracs) - min(sens_fracs):.3f}")
    # 22. Negative / falsification check: emergence rate (non-signal→signal
    # under indication match) should be much smaller than the vanishing rate.
    # If they were similar, indication matching would be symmetric noise rather
    # than a systematic de-confounder.
    check(r["emerge_fraction"] < r["vanish_fraction"],
          f"emerge_fraction ({r['emerge_fraction']:.3f}) not smaller than "
          f"vanish_fraction ({r['vanish_fraction']:.3f}) — falsification failed")
    # 23. Null-expected vanish rate is in [0, 1] and below the observed rate.
    check(0.0 <= r["null_expected_vanish_fraction"] <= 1.0,
          f"null_expected_vanish_fraction out of range: "
          f"{r['null_expected_vanish_fraction']}")
    check(r["excess_vanish_fraction_vs_null"] > 0,
          f"excess over null not positive: "
          f"{r['excess_vanish_fraction_vs_null']:.3f}")
    # 24. Cohen's-d-style plausibility on the attenuation scale: std across
    # signals should be finite and the effect (|median| / spread) < 5 ~ no
    # single-point dominance.
    iqr_lo, iqr_hi = r["attenuation_log2_iqr"]
    iqr_width = iqr_hi - iqr_lo
    if iqr_width > 0:
        effect_size = abs(r["attenuation_log2_median"]) / iqr_width
        check(effect_size < 5.0,
              f"effect size implausibly large: {effect_size:.3f}")
    else:
        # Zero IQR is itself degenerate.
        check(False, f"degenerate attenuation IQR: [{iqr_lo}, {iqr_hi}]")
    # 25. Unique primary indications present.
    check(r["unique_primary_indications"] >= 2,
          f"too few unique primary indications: {r['unique_primary_indications']}")
    # 26. Drug-list ↔ indication-map consistency.
    ind_map = r.get("drug_primary_indication_map", {})
    check(set(ind_map.keys()) == set(r["drug_list"]),
          "drug_primary_indication_map keys do not match drug_list")
    # 27. Attenuation CI contains the median.
    att_lo, att_hi = r["attenuation_log2_median_cluster_boot_95ci"]
    if not math.isnan(att_lo) and not math.isnan(att_hi):
        check(att_lo <= r["attenuation_log2_median"] <= att_hi,
              f"attenuation median outside its bootstrap CI: "
              f"{r['attenuation_log2_median']} not in [{att_lo}, {att_hi}]")
    # 28. Report runtime fields are present and sensible.
    check(r["data_window"].count("-") == 1 and len(r["data_window"]) >= 9,
          f"data_window malformed: {r['data_window']}")
    # 29. Sensitivity: at least one cell uses the main thresholds — round-trip
    # check that the headline scan actually includes the headline config.
    hit = any(row["prr_threshold"] == r["prr_threshold"]
              and row["min_count"] == r["min_count_threshold"]
              for row in r["sensitivity"])
    check(hit, "sensitivity scan does not include the headline (PRR, min-count) cell")
    # 30. Sensitivity headline consistency: the row matching the headline
    # thresholds should reproduce the headline vanish fraction (modulo fp).
    for row in r["sensitivity"]:
        if (row["prr_threshold"] == r["prr_threshold"]
                and row["min_count"] == r["min_count_threshold"]):
            if isinstance(row["vanish_fraction"], (int, float)) \
                    and not math.isnan(row["vanish_fraction"]):
                check(abs(row["vanish_fraction"] - r["vanish_fraction"]) < 1e-6,
                      f"sensitivity row disagrees with headline: "
                      f"{row['vanish_fraction']} vs {r['vanish_fraction']}")
            break
    # 31. Null vanish fraction should be much lower than observed (else the
    # indication-match finding could be explained by sampling noise alone).
    check(r["vanish_fraction"] > r["null_expected_vanish_fraction"] + 0.01,
          f"observed vanishing ({r['vanish_fraction']:.3f}) not meaningfully "
          f"above null ({r['null_expected_vanish_fraction']:.3f})")
    # 32. Top-signals-detail structural check: each row has the seven expected keys.
    required_keys = {"drug", "event", "indication", "prr_w", "prr_m", "a_w", "vanished"}
    for i, row in enumerate(r["top_signals_detail"][:5]):
        missing = required_keys - set(row.keys())
        check(not missing,
              f"top_signals_detail row {i} missing keys: {missing}")
    # 33. Signal counts add up: vanished + surviving = total.
    check(r["vanished_signals"] <= r["top_signals"],
          f"vanished_signals ({r['vanished_signals']}) > "
          f"top_signals ({r['top_signals']})")
    # 34. Drug-event pairs bounded above by drugs * events-per-drug.
    max_pairs = r["drugs_in_analysis"] * 30  # N_TOP_EVENTS_PER_DRUG
    check(r["drug_event_pairs"] <= max_pairs,
          f"more pairs ({r['drug_event_pairs']}) than drugs*events cap "
          f"({max_pairs}) — suggests double-counting")

    if failures:
        print("VERIFY FAILED:", file=sys.stderr)
        for msg in failures:
            print(f"  - {msg}", file=sys.stderr)
        print(f"{len(failures)}/{n_checks} checks failed.", file=sys.stderr)
        return 1
    print(f"verify: {n_checks}/{n_checks} checks passed")
    print(f"  vanish_fraction = {r['vanish_fraction']:.3f} "
          f"[{r['vanish_fraction_cluster_boot_95ci'][0]:.3f}, "
          f"{r['vanish_fraction_cluster_boot_95ci'][1]:.3f}]")
    print(f"  excess vs null  = {r['excess_vanish_fraction_vs_null']:+.3f}")
    print(f"  MC p-value      = {r['monte_carlo_p_value']:.4f}")
    print("ALL CHECKS PASSED")
    return 0


# ═══════════════════════════════════════════════════════════════
# Orchestration
# ═══════════════════════════════════════════════════════════════

def audit_cache_manifest():
    """Verify that each file referenced by the cache manifest exists and
    hashes to the recorded body_sha256. Used by --verify-cache to detect
    partial/truncated downloads before a rerun."""
    if not os.path.exists(CACHE_MANIFEST_FILE):
        print("no cache manifest found — nothing to audit", file=sys.stderr)
        return 0
    with open(CACHE_MANIFEST_FILE, "r", encoding="utf-8") as f:
        entries = json.load(f)
    bad = 0
    for entry in entries:
        path = os.path.join(CACHE_DIR, entry["cache_file"])
        if not os.path.exists(path):
            print(f"MISSING: {path}", file=sys.stderr)
            bad += 1
            continue
        with open(path, "r", encoding="utf-8") as f:
            body = f.read().encode("utf-8")
        # Re-canonicalise via json.dumps(sort_keys=True) to match how the
        # bytes were hashed at write time.
        try:
            canon = json.dumps(json.loads(body), sort_keys=True).encode("utf-8")
        except Exception:
            print(f"UNPARSEABLE: {path}", file=sys.stderr)
            bad += 1
            continue
        got = hashlib.sha256(canon).hexdigest()
        if got != entry["body_sha256"]:
            print(f"HASH-MISMATCH: {path} (expected {entry['body_sha256']}, got {got})",
                  file=sys.stderr)
            bad += 1
    if bad:
        print(f"{bad} cache entries failed verification", file=sys.stderr)
        return 1
    print(f"cache audit: {len(entries)} entries, all hashes match")
    return 0


if __name__ == "__main__":
    if len(sys.argv) > 1 and sys.argv[1] == "--verify":
        sys.exit(verify())
    if len(sys.argv) > 1 and sys.argv[1] == "--verify-cache":
        sys.exit(audit_cache_manifest())
    random.seed(SEED)
    try:
        os.makedirs(CACHE_DIR, exist_ok=True)
    except OSError as e:
        print(f"ERROR: cannot create cache directory {CACHE_DIR!r}: {e}",
              file=sys.stderr)
        sys.exit(3)
    try:
        data = load_data()
    except urllib.error.URLError as e:
        print(f"ERROR: network failure while fetching openFDA data: {e}",
              file=sys.stderr)
        print("  Check outbound HTTPS to api.fda.gov. "
              "Delete cache/ and retry if partial responses are suspected.",
              file=sys.stderr)
        sys.exit(2)
    except RuntimeError as e:
        print(f"ERROR: data validation failed: {e}", file=sys.stderr)
        sys.exit(2)
    except OSError as e:
        print(f"ERROR: filesystem error while loading/caching data: {e}",
              file=sys.stderr)
        sys.exit(3)
    try:
        results = run_analysis(data)
    except Exception as e:
        print(f"ERROR: analysis failed: {type(e).__name__}: {e}", file=sys.stderr)
        sys.exit(4)
    try:
        generate_report(results)
    except OSError as e:
        print(f"ERROR: could not write results/report: {e}", file=sys.stderr)
        sys.exit(5)
    print("ANALYSIS COMPLETE")
    print(f"Results: {RESULTS_FILE}")
    print(f"Report:  {REPORT_FILE}")
    print("Limitations: "
          "(1) FAERS is spontaneous reporting, not a cohort — reports have "
          "no denominator of exposed patients and are subject to stimulated, "
          "notoriety, and under-reporting biases. "
          "(2) 'Primary indication' is defined as the most-reported "
          "drugindication term per drug, which ignores patients reported with "
          "multiple indications and under-reports true co-indication. "
          "(3) openFDA count queries cap at 500 rows; rare events outside the "
          "cap are back-filled by exact searches but may still miss long-tail "
          "co-occurrences. "
          "(4) This analysis quantifies attenuation of disproportionality "
          "signals under one specific confounder (primary indication). It does "
          "NOT establish causality, rule out other confounders (e.g. age, "
          "sex, concomitant drugs), nor differentiate real pharmacology from "
          "diagnostic surveillance bias. "
          "(5) medicinalproduct.exact is heterogeneous: brand (HUMIRA) and "
          "generic (ADALIMUMAB) forms are counted separately, which deflates "
          "individual drug totals versus a substance-name rollup. "
          "(6) The binomial null in null_vanish_probability() assumes the "
          "within-indication cells are independent draws; the cluster "
          "bootstrap CI is the appropriate envelope for cross-drug correlation "
          "and is strictly wider than the Wilson CI. "
          "(7) Late reporting continuously updates historical windows, so a "
          "byte-exact reproduction requires the frozen cache/ directory; "
          "headline values drift by up to ~2 pp over long intervals. "
          "(8) Results apply to high-volume single-indication drugs; drugs "
          "with diverse indication profiles are excluded by design and their "
          "robustness under indication-matching is NOT characterised here.")
SCRIPT_EOF
```

**Expected output:** no stdout. Writes `analyze.py` to the workspace.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_disproportionality-signals-that-vanish-under-indication-cont && python3 analyze.py
```

**Expected output:** sectioned progress `[1/7]` through `[7/7]`, ending with `ANALYSIS COMPLETE`. Runtime is typically 3-8 minutes on the first run (mostly openFDA API calls) and under 1 minute on cached reruns. Produces `results.json` and `report.md` in the workspace.

**Key lines to look for in stdout:**

- `total_db_reports=1,523,664` (or similar — 2022 FAERS should hold about 1.5 million reports).
- `top signals (whole-DB): N` with N ≥ 20.
- `vanished under indication match: M (P%)` with 0 ≤ P ≤ 100.
- `null-expected vanishing rate` and `Monte-Carlo p-value`.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_disproportionality-signals-that-vanish-under-indication-cont && python3 analyze.py --verify
```

**Expected output:** a line `verify: N/N checks passed` (N ≥ 30), followed by a one-line summary of the headline statistic and its cluster-bootstrap 95% CI, terminating with `ALL CHECKS PASSED`. Exit code 0.

## Step 5: (Optional) Audit the Cache

```bash
cd /tmp/claw4s_auto_disproportionality-signals-that-vanish-under-indication-cont && python3 analyze.py --verify-cache
```

**Expected output:** `cache audit: K entries, all hashes match` where K is the number of cached API responses (typically ~140 for a full cold run). Any `HASH-MISMATCH`, `MISSING`, or `UNPARSEABLE` line implies the cache is corrupt — delete `cache/` and rerun.

## Success Criteria

1. **Normal run**: `analyze.py` runs to completion and stdout ends with `ANALYSIS COMPLETE`, exit code 0.
2. **Verify run**: `analyze.py --verify` prints `ALL CHECKS PASSED` after at least 30 machine-checkable assertions, exit code 0.
2a. **Cache audit (optional)**: `analyze.py --verify-cache` prints `all hashes match`, exit code 0.
3. **Artifact contract**: `results.json` contains the keys `vanish_fraction`, `vanish_fraction_cluster_boot_95ci`, `vanish_fraction_wilson_95ci`, `null_expected_vanish_fraction`, `monte_carlo_p_value`, `attenuation_log2_median`, `sensitivity`, `top_signals_detail`, and `seed`.
4. **Effect size plausibility**: `vanish_fraction` is between 0.10 and 0.90. Real FAERS data shows a non-trivial share of signals attenuate under indication matching, but complete vanishing would indicate the whole-DB control is uninformative, and zero vanishing would contradict the well-documented confounding-by-indication effect.
5. **CI coverage**: `vanish_fraction_cluster_boot_95ci` is a non-degenerate interval that contains the point estimate and whose width is > 0.01 and < 1.0.
6. **Sensitivity invariance**: across all (PRR threshold, min-count) cells of the sensitivity scan, `vanish_fraction` stays in the plausible band — i.e. no single parameter choice drives the headline.
7. **Falsification**: `emerge_fraction` (non-signal→signal under matching) is strictly less than `vanish_fraction`. If indication-matching were symmetric noise, these rates would be comparable.
8. **Null comparison**: `excess_vanish_fraction_vs_null > 0` and the Monte-Carlo p-value is < 0.5 (most runs give a much smaller p-value).

## Failure Conditions

- **Network failure**: `urllib.error.URLError` or repeated HTTP errors → openFDA is unreachable. The script retries five times per call with exponential back-off; if the call still fails, it writes a clear error to stderr and exits with code 2. Check outbound HTTPS connectivity to `api.fda.gov`.
- **Empty or partial cache**: if a cached JSON file is truncated or was written during a rate-limit response, reruns can produce `VERIFY FAILED` with a suspicious `pairs_skipped_due_to_negative_cells` or `top_signals_detail` empty. Delete `cache/` and retry.
- **Data drift in openFDA**: if `top_signals < 20`, the filter thresholds are likely still in range but openFDA may have changed its `drugindication` field population or renamed top indications. Lowering `N_TOP_DRUGS` is *not* the right fix — investigate the `drug_blocks` diagnostics to see which drugs were filtered.
- **Reproducibility drift over long intervals**: openFDA continues to receive late-arriving reports for any `receivedate` window, so the same query run a year apart may differ slightly (typically within ±2 percentage points on the headline). This is a feature of the data source, not a defect of the skill.
- **Rate-limit exhaustion**: unauthenticated openFDA caps at 240 req/min. On a cold start with a very slow network the `time.sleep(API_RATE_SLEEP_SEC)` may not be enough; the 429 back-off path handles this but under extreme load can time out after 5 retries.
- **Misapplication as causal analysis**: the skill quantifies *attenuation of disproportionality* under one confounder (primary indication). A large `vanish_fraction` does NOT prove the original signals were false; it proves they were not robust to this specific adjustment. Other confounders (age, sex, concomitant drugs, reporter type) are not addressed.
- **Degenerate indication field**: for drug classes where `drugindication` is largely "Product used for unknown indication" (e.g. many OTC products, consumer-report-only drugs), the primary-indication rule excludes them entirely. The report's `n_drugs_excluded_because_top_indication_is_unknown` field quantifies this.

## Reproducibility Notes

- Random seed is fixed at 42. All stochastic operations (binomial sampling, bootstrap, Monte-Carlo) are seeded.
- openFDA receives late-reported cases in perpetuity, so an identical receivedate window can drift slightly between queries months apart. In practice the headline percentage is stable to within ±2 percentage points for a 12-month window that is at least one year old. Pinning `DATE_START` and `DATE_END` gives a closed cohort as of the query date.
- Cached API responses are written to `cache/fda_<hash>.json`. Downstream reruns reuse the cache deterministically. Deleting the cache forces a fresh download.