---
name: "Disproportionality Metric Benchmark Against FDA Label Mentions"
description: "Compares four pharmacovigilance disproportionality metrics (PRR, ROR, IC, EBGM) on their ability to predict whether adverse-event terms appear in current FDA drug labels. Uses openFDA FAERS counts and SPL label text for ~50 high-volume drugs, reports AUC/PPV with bootstrap CIs, permutation-tested significance, and sensitivity to cell-size filters."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "pharmacovigilance", "faers", "disproportionality", "prr", "ror", "ebgm", "ic", "biomedical"]
python_version: ">=3.8"
dependencies: []
---

# Disproportionality Metric Benchmark Against FDA Label Mentions

## When to Use This Skill

Use this skill when you need to investigate how well pharmacovigilance disproportionality
metrics (PRR, ROR, IC, EBGM) computed from the FDA Adverse Event Reporting System (FAERS)
align with adverse events explicitly listed on current FDA drug labels — an operational
proxy for "known" adverse reactions. The analysis yields AUC, PPV at low-recall thresholds,
bootstrap confidence intervals, a permutation null, a negative-control / falsification
check on shuffled labels, and sensitivity analyses across cell-size filters, label section
sets, and indication overlap.

Agents should invoke this skill when a user asks any of:
- "Which disproportionality metric is best at recovering known adverse reactions?"
- "Do PRR, ROR, IC, and EBGM agree in practice on FAERS?"
- "How sensitive are pharmacovigilance signal-detection rankings to cell-size or label
  scope choices?"

## Research Question

We test whether disproportionality metrics computed from spontaneous-report data can
recover whether an adverse-event term is listed on the current FDA label. Specifically:

> **Hypothesis (H1):** For a drug-event pair (D, E), the disproportionality scores
> PRR(D, E), ROR(D, E), IC(D, E), EBGM(D, E) are independent of whether E appears in
> D's current FDA label (null), versus predictive (alternative, AUC > 0.5) or
> anti-predictive (alternative, AUC < 0.5).

This is falsifiable (permutation test on shuffled labels collapses AUC to 0.5) and
specific (uses openFDA FAERS counts + openFDA SPL labels for the top 50 drugs by report
volume). The analysis reports point estimates with bootstrap 95% CIs, a BH-FDR-corrected
permutation p-value per metric, pairwise paired-bootstrap AUC differences between
metrics, and a sensitivity sweep over the operational definition of "positive" (label
section set, indication overlap).

## Prerequisites

- **Python:** 3.8 or newer. Standard library only — no `pip install` required, no
  third-party packages imported.
- **Network:** HTTPS access to `https://api.fda.gov`. First run issues ~150 API calls
  downloading ~1–2 MB of JSON. Retries with exponential backoff on transient errors.
- **Disk:** ~5 MB for the workspace including the on-disk cache.
- **Runtime:** 5–10 minutes on first run (network-bound); <5 seconds from cache.
- **Environment variables:** None required. The script respects the `PYTHONHASHSEED`
  environment convention but does not rely on it for determinism.
- **Random seed:** all stochastic steps (bootstrap, permutation) use
  `random.Random(RANDOM_SEED=42)`. Same cache + same code ⇒ identical `results.json`.

## Adaptation Guidance

To adapt this skill to a new pharmacovigilance comparison or even a different tabular
"signal-versus-ground-truth" problem, modify only the `# DOMAIN CONFIGURATION` block of the
analysis script:

- **Change the data source**: update `API_EVENT_URL`, `API_LABEL_URL`, `DRUG_COUNT_FIELD`,
  `EVENT_COUNT_FIELD`, and `LABEL_DRUG_FIELD` to point to a different registry. The helper
  functions `_fetch_json()`, `load_data()`, and `_extract_label_text()` encapsulate all
  openFDA-specific parsing.
- **Change metrics**: the functions `compute_prr()`, `compute_ror()`, `compute_ic()`,
  `compute_ebgm()` each take the 2×2 cell counts `(a, b, c, d)` and are independent; add a
  new metric by writing a sibling function and registering it in the `METRICS` dict.
- **Change the ground truth**: `_is_event_in_label()` defines the positive-class criterion
  (substring match on MedDRA preferred term in one or more label sections, controlled by
  `LABEL_SECTIONS`). Replace with any function mapping a drug-event pair to a 0/1 label.
- **What stays the same**: `run_analysis()` is domain-agnostic — it computes AUC (via the
  Mann-Whitney U identity), PPV at a list of recall thresholds, bootstrap CIs, a permutation
  null, and a sensitivity sweep over cell-size thresholds. None of it references
  pharmacovigilance terminology.
- **Tuning knobs** in the CONFIGURATION block: `N_DRUGS`, `N_EVENTS_PER_DRUG`,
  `MIN_CELL_SIZE_MAIN`, `SENSITIVITY_MIN_CELL_SIZES`, `BOOTSTRAP_N`, `PERMUTATION_N`,
  `RECALL_THRESHOLDS`, `RANDOM_SEED`.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g/cache
```

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

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g/analyze.py
#!/usr/bin/env python3
"""
Disproportionality Metric Benchmark Against FDA Label Mentions.

For a set of high-volume drugs in openFDA FAERS, this script:
  1. Downloads per-drug report totals, per-drug top events with joint counts,
     global per-event totals, the overall report corpus size, and current-label
     adverse-reactions text from openFDA.
  2. Builds 2x2 contingency cells (a, b, c, d) for each drug-event pair.
  3. Computes four disproportionality metrics: PRR, ROR, IC, EBGM.
  4. Labels each pair positive iff the event's MedDRA preferred term appears in
     any configured label section.
  5. Benchmarks each metric's separability of positives from negatives using
     AUC (Mann-Whitney U) and PPV at operational recall thresholds.
  6. Estimates 95 percent bootstrap confidence intervals for each AUC.
  7. Runs a permutation null (shuffled labels) to test significance.
  8. Runs a sensitivity sweep over minimum cell-size filters.

Python 3.8+ standard library only. No pip install.

Data source: openFDA (FAERS + SPL). https://open.fda.gov/apis/
"""

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

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

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")

API_EVENT_URL = "https://api.fda.gov/drug/event.json"
API_LABEL_URL = "https://api.fda.gov/drug/label.json"

DRUG_COUNT_FIELD = "patient.drug.openfda.generic_name.exact"
EVENT_COUNT_FIELD = "patient.reaction.reactionmeddrapt.exact"
LABEL_DRUG_FIELD = "openfda.generic_name.exact"

# Label sections scanned for positive-class evidence.
LABEL_SECTIONS = ["adverse_reactions", "warnings_and_cautions", "boxed_warning"]
# Label sections used to build an "indication exclusion" filter: events whose
# text appears in these sections are considered likely indications (reasons
# for treatment) rather than adverse reactions, and excluded in a sensitivity
# analysis.
INDICATION_SECTIONS = ["indications_and_usage"]

# Number of top drugs (by total FAERS report count) to study.
N_DRUGS = 50
# Top events per drug (by joint report count with that drug). openFDA caps
# the `count` aggregation at limit < 1000; we use 150 to get a useful range
# of joint-count cell sizes from small (a~5) to large (a~50,000) while
# keeping the bootstrap permutation loops tractable in pure-Python stdlib.
N_EVENTS_PER_DRUG = 150
# Global top events used to build the marginal dictionary.
N_EVENTS_GLOBAL = 999
# Events with fewer than this many global reports are dropped from analysis
# (they produce unstable disproportionality ratios).
MIN_EVENT_GLOBAL_COUNT = 100
# A drug-event pair is included in the main analysis only if a >= this.
MIN_CELL_SIZE_MAIN = 3
# Sensitivity sweep over minimum cell sizes (main analysis uses MIN_CELL_SIZE_MAIN).
# Thresholds span the observed joint-count range (openFDA returns only a drug's top
# events, so raw a values are large: hundreds to tens of thousands). Values
# <= min(a) keep the full table; higher values filter progressively more pairs.
SENSITIVITY_MIN_CELL_SIZES = [1, 3, 1000, 5000, 10000, 50000]

# Bootstrap and permutation configuration.
BOOTSTRAP_N = 500
PERMUTATION_N = 1000
RECALL_THRESHOLDS = [0.01, 0.05, 0.10, 0.20]
RANDOM_SEED = 42

# Noise-robust heuristic tokens that often appear in labels but are rarely
# useful matches. When the MedDRA term equals one of these, we skip it.
GENERIC_TOKEN_SKIP = {
    "DEATH", "DRUG INEFFECTIVE", "OFF LABEL USE", "CONDITION AGGRAVATED",
    "NO ADVERSE EVENT", "PAIN", "UNEVALUABLE EVENT", "PRODUCT QUALITY ISSUE",
    "MALAISE"
}

# API backoff and pacing.
REQUEST_DELAY_SECONDS = 0.25
MAX_RETRIES = 4
INITIAL_BACKOFF = 1.5

# ═══════════════════════════════════════════════════════════════
# HELPERS - statistical utilities and I/O
# ═══════════════════════════════════════════════════════════════

def sha256_hex(data_bytes):
    return hashlib.sha256(data_bytes).hexdigest()


def _url_encode_term(term):
    return urllib.parse.quote('"' + term + '"', safe="")


def _fetch_raw(url):
    last_err = None
    for attempt in range(MAX_RETRIES):
        try:
            req = urllib.request.Request(url)
            req.add_header("User-Agent", "Claw4S-Research/1.0 (disproportionality benchmark)")
            with urllib.request.urlopen(req, timeout=60) as resp:
                return resp.read()
        except urllib.error.HTTPError as e:
            # openFDA returns 404 when a search matches zero documents.
            if e.code == 404:
                return b'{"results": [], "meta": {"results": {"total": 0}}}'
            last_err = e
        except (urllib.error.URLError, OSError) as e:
            last_err = e
        wait = INITIAL_BACKOFF * (2 ** attempt)
        time.sleep(wait)
    raise RuntimeError(f"Failed to fetch {url} after {MAX_RETRIES} attempts: {last_err}")


def _fetch_json(url, cache_name):
    """Fetch a URL with on-disk caching and SHA256 verification."""
    cache_file = os.path.join(CACHE_DIR, cache_name + ".json")
    hash_file = os.path.join(CACHE_DIR, cache_name + ".sha256")
    if os.path.exists(cache_file) and os.path.exists(hash_file):
        with open(cache_file, "rb") as fh:
            data = fh.read()
        with open(hash_file, "r") as fh:
            expected = fh.read().strip()
        if sha256_hex(data) == expected:
            return json.loads(data.decode("utf-8"))
    raw = _fetch_raw(url)
    with open(cache_file, "wb") as fh:
        fh.write(raw)
    with open(hash_file, "w") as fh:
        fh.write(sha256_hex(raw))
    time.sleep(REQUEST_DELAY_SECONDS)
    return json.loads(raw.decode("utf-8"))


def _safe_filename(name):
    return re.sub(r"[^A-Za-z0-9_-]+", "_", name)[:80]


# ------- Disproportionality metrics (2x2 table: a, b, c, d) -------
# a = reports with drug D and event E
# b = reports with drug D and NOT event E
# c = reports without drug D and event E
# d = reports without drug D and NOT event E

def compute_prr(a, b, c, d):
    """Proportional Reporting Ratio (Evans 2001) with Haldane-Anscombe +0.5."""
    a_, b_, c_, d_ = a + 0.5, b + 0.5, c + 0.5, d + 0.5
    return (a_ / (a_ + b_)) / (c_ / (c_ + d_))


def compute_ror(a, b, c, d):
    """Reporting Odds Ratio (van Puijenbroek 2002) with +0.5."""
    a_, b_, c_, d_ = a + 0.5, b + 0.5, c + 0.5, d + 0.5
    return (a_ * d_) / (b_ * c_)


def compute_ic(a, b, c, d):
    """BCPNN Information Component (Bate 1998), uniform prior shrinkage.

    IC = log2((a + 0.5) / (E + 0.5))
    where E = (a+b)(a+c) / N is the expected count under independence.
    """
    n = a + b + c + d
    if n == 0:
        return 0.0
    expected = (a + b) * (a + c) / n
    return math.log2((a + 0.5) / (expected + 0.5))


def compute_ebgm(a, b, c, d, prior):
    """Empirical Bayes Geometric Mean (DuMouchel 1999), method-of-moments prior.

    The prior is a single gamma(alpha, beta) estimated from the observed
    distribution of RR = a / E across all pairs. Posterior mean under a
    conjugate gamma-Poisson is (a + alpha) / (E + beta).
    """
    n = a + b + c + d
    if n == 0:
        return 0.0
    expected = (a + b) * (a + c) / n
    alpha, beta = prior
    # Guard against numerical pathologies.
    if expected < 1e-9:
        expected = 1e-9
    return (a + alpha) / (expected + beta)


def _fit_gamma_mom(values):
    """Fit shape=alpha, rate=beta of a gamma via method of moments.

    For a Gamma(alpha, beta) with E[X] = alpha/beta and Var[X] = alpha/beta^2:
      beta  = mean / var
      alpha = mean * beta
    """
    if not values:
        return (0.5, 0.5)
    mean = sum(values) / len(values)
    if len(values) < 2 or mean <= 0:
        return (0.5, 0.5)
    var = sum((v - mean) ** 2 for v in values) / (len(values) - 1)
    if var <= 0:
        return (0.5, 0.5)
    beta = mean / var
    alpha = mean * beta
    # Clamp to sensible range to avoid runaway priors.
    alpha = min(max(alpha, 0.01), 100.0)
    beta = min(max(beta, 0.01), 100.0)
    return (alpha, beta)


def _mann_whitney_auc(positive_scores, negative_scores):
    """AUC via Mann-Whitney U, counting ties as half. O(n log n)."""
    if not positive_scores or not negative_scores:
        return float("nan")
    n_pos = len(positive_scores)
    n_neg = len(negative_scores)
    combined = [(s, 1) for s in positive_scores] + [(s, 0) for s in negative_scores]
    combined.sort(key=lambda x: x[0])
    # Rank with tie-averaging.
    ranks = [0.0] * len(combined)
    i = 0
    while i < len(combined):
        j = i
        while j + 1 < len(combined) and combined[j + 1][0] == combined[i][0]:
            j += 1
        avg_rank = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[k] = avg_rank
        i = j + 1
    r_pos = sum(r for r, (_, lbl) in zip(ranks, combined) if lbl == 1)
    u = r_pos - n_pos * (n_pos + 1) / 2.0
    return u / (n_pos * n_neg)


def auc_roc(scores, labels):
    pos = [s for s, l in zip(scores, labels) if l == 1]
    neg = [s for s, l in zip(scores, labels) if l == 0]
    return _mann_whitney_auc(pos, neg)


def ppv_at_recall(scores, labels, target_recall):
    """At the smallest top-k that achieves recall >= target_recall, return PPV."""
    total_pos = sum(labels)
    if total_pos == 0:
        return float("nan")
    order = sorted(range(len(scores)), key=lambda i: -scores[i])
    cum_tp = 0
    for rank, idx in enumerate(order, start=1):
        if labels[idx] == 1:
            cum_tp += 1
        if cum_tp / total_pos >= target_recall:
            return cum_tp / rank
    return cum_tp / len(order) if order else float("nan")


def bootstrap_auc_ci(scores, labels, n_iter, rng, alpha=0.05):
    n = len(scores)
    if n < 10:
        point = auc_roc(scores, labels)
        return (point, point, point)
    aucs = []
    for _ in range(n_iter):
        idx = [rng.randrange(n) for _ in range(n)]
        boot_scores = [scores[i] for i in idx]
        boot_labels = [labels[i] for i in idx]
        if sum(boot_labels) == 0 or sum(boot_labels) == n:
            continue
        aucs.append(auc_roc(boot_scores, boot_labels))
    aucs.sort()
    if not aucs:
        point = auc_roc(scores, labels)
        return (point, point, point)
    lo = aucs[int(alpha / 2 * len(aucs))]
    hi = aucs[int((1 - alpha / 2) * len(aucs)) - 1]
    return (auc_roc(scores, labels), lo, hi)


def permutation_auc_pvalue(scores, labels, n_iter, rng, two_sided=True):
    """Permutation p-value for AUC.

    If two_sided (default) the statistic is |AUC - 0.5|, so the test rejects
    the null of "score uninformative about label" when the metric has
    discrimination in *either* direction (predictive or anti-predictive).
    If not two_sided the original one-sided AUC >= observed test is used.
    """
    observed = auc_roc(scores, labels)
    if not (0 <= observed <= 1):
        return float("nan")
    target = abs(observed - 0.5) if two_sided else observed
    n_ge = 0
    shuffled = list(labels)
    for _ in range(n_iter):
        rng.shuffle(shuffled)
        perm_auc = auc_roc(scores, shuffled)
        stat = abs(perm_auc - 0.5) if two_sided else perm_auc
        if stat >= target:
            n_ge += 1
    return (n_ge + 1) / (n_iter + 1)


def paired_bootstrap_auc_diff(scores_a, scores_b, labels, n_iter, rng, alpha=0.05):
    """Paired bootstrap for AUC_a - AUC_b on the same drug-event pairs."""
    n = len(scores_a)
    assert n == len(scores_b) == len(labels)
    diffs = []
    for _ in range(n_iter):
        idx = [rng.randrange(n) for _ in range(n)]
        la = [labels[i] for i in idx]
        if sum(la) == 0 or sum(la) == n:
            continue
        sa = [scores_a[i] for i in idx]
        sb = [scores_b[i] for i in idx]
        diffs.append(auc_roc(sa, la) - auc_roc(sb, la))
    diffs.sort()
    if not diffs:
        obs = auc_roc(scores_a, labels) - auc_roc(scores_b, labels)
        return (obs, obs, obs)
    lo = diffs[int(alpha / 2 * len(diffs))]
    hi = diffs[int((1 - alpha / 2) * len(diffs)) - 1]
    obs = auc_roc(scores_a, labels) - auc_roc(scores_b, labels)
    return (obs, lo, hi)


def benjamini_hochberg(pvals):
    """Return BH-FDR-adjusted q-values for a list of p-values (same order)."""
    n = len(pvals)
    if n == 0:
        return []
    order = sorted(range(n), key=lambda i: pvals[i])
    ranked = [pvals[i] for i in order]
    # Monotone from the end.
    adj = [0.0] * n
    running_min = 1.0
    for k in range(n - 1, -1, -1):
        val = ranked[k] * n / (k + 1)
        if val < running_min:
            running_min = val
        adj[k] = min(1.0, running_min)
    out = [0.0] * n
    for i, orig_idx in enumerate(order):
        out[orig_idx] = adj[i]
    return out


def spearman_rho(x, y):
    """Spearman rank correlation on two equal-length sequences."""
    n = len(x)
    if n < 3:
        return float("nan")
    def rank(v):
        order = sorted(range(n), key=lambda i: v[i])
        r = [0.0] * n
        i = 0
        while i < n:
            j = i
            while j + 1 < n and v[order[j + 1]] == v[order[i]]:
                j += 1
            avg = (i + j) / 2.0 + 1.0
            for k in range(i, j + 1):
                r[order[k]] = avg
            i = j + 1
        return r
    rx = rank(x)
    ry = rank(y)
    mx = sum(rx) / n
    my = sum(ry) / n
    num = sum((rx[i] - mx) * (ry[i] - my) for i in range(n))
    denom_x = math.sqrt(sum((rx[i] - mx) ** 2 for i in range(n)))
    denom_y = math.sqrt(sum((ry[i] - my) ** 2 for i in range(n)))
    if denom_x == 0 or denom_y == 0:
        return float("nan")
    return num / (denom_x * denom_y)


# ═══════════════════════════════════════════════════════════════
# DATA LOADING - domain-specific openFDA calls
# ═══════════════════════════════════════════════════════════════

def _normalize_for_match(text):
    text = text.lower()
    text = re.sub(r"[^a-z0-9 ]+", " ", text)
    text = re.sub(r"\s+", " ", text).strip()
    return text


def _extract_label_text(label_record, sections):
    parts = []
    for sec in sections:
        val = label_record.get(sec)
        if val:
            if isinstance(val, list):
                parts.extend(str(v) for v in val)
            else:
                parts.append(str(val))
    return _normalize_for_match(" ".join(parts))


def _is_event_in_label(event_term, label_text_normalized):
    norm = _normalize_for_match(event_term)
    if not norm:
        return False
    # Require a whole-token or whole-phrase match to avoid matching short
    # substrings inside unrelated words (e.g. "pain" inside "painting").
    pattern = r"\b" + re.escape(norm) + r"\b"
    return re.search(pattern, label_text_normalized) is not None


def load_data():
    print("  [load] fetching top drugs by report volume")
    url = f"{API_EVENT_URL}?count={DRUG_COUNT_FIELD}&limit={N_DRUGS}"
    resp = _fetch_json(url, "drugs_top")
    drugs = [r["term"] for r in resp.get("results", [])][:N_DRUGS]
    drug_totals = {r["term"]: int(r["count"]) for r in resp.get("results", [])}
    print(f"  [load] got {len(drugs)} drugs")

    print("  [load] fetching global event marginals")
    url = f"{API_EVENT_URL}?count={EVENT_COUNT_FIELD}&limit={N_EVENTS_GLOBAL}"
    resp = _fetch_json(url, "events_global")
    event_totals = {r["term"]: int(r["count"]) for r in resp.get("results", [])
                    if int(r["count"]) >= MIN_EVENT_GLOBAL_COUNT}
    print(f"  [load] kept {len(event_totals)} events with count >= {MIN_EVENT_GLOBAL_COUNT}")

    print("  [load] fetching total report count")
    resp = _fetch_json(f"{API_EVENT_URL}?limit=1", "total")
    total_reports = int(resp.get("meta", {}).get("results", {}).get("total", 0))
    event_api_last_updated = resp.get("meta", {}).get("last_updated", "")
    print(f"  [load] total FAERS reports: {total_reports:,}")
    print(f"  [load] event API last_updated: {event_api_last_updated}")

    label_api_last_updated = ""
    try:
        meta_resp = _fetch_json(f"{API_LABEL_URL}?limit=1", "label_meta")
        label_api_last_updated = meta_resp.get("meta", {}).get("last_updated", "")
        print(f"  [load] label API last_updated: {label_api_last_updated}")
    except Exception:
        pass

    print(f"  [load] fetching per-drug top events (N={len(drugs)} drugs)")
    pair_counts = {}  # (drug, event) -> joint count
    for i, drug in enumerate(drugs, start=1):
        search = f"search={DRUG_COUNT_FIELD}:{_url_encode_term(drug)}"
        url = f"{API_EVENT_URL}?{search}&count={EVENT_COUNT_FIELD}&limit={N_EVENTS_PER_DRUG}"
        resp = _fetch_json(url, f"drug_{_safe_filename(drug)}_events")
        for r in resp.get("results", []):
            event = r["term"]
            if event in event_totals:
                pair_counts[(drug, event)] = int(r["count"])
        if i % 10 == 0:
            print(f"    [load] {i}/{len(drugs)} drugs fetched")
    print(f"  [load] got {len(pair_counts)} drug-event pairs")

    print("  [load] fetching labels for each drug")
    labels = {}
    indications = {}
    for i, drug in enumerate(drugs, start=1):
        search = f"search={LABEL_DRUG_FIELD}:{_url_encode_term(drug)}"
        url = f"{API_LABEL_URL}?{search}&limit=1"
        resp = _fetch_json(url, f"label_{_safe_filename(drug)}")
        results = resp.get("results", [])
        if results:
            labels[drug] = _extract_label_text(results[0], LABEL_SECTIONS)
            indications[drug] = _extract_label_text(results[0], INDICATION_SECTIONS)
        else:
            labels[drug] = ""
            indications[drug] = ""
        if i % 10 == 0:
            print(f"    [load] {i}/{len(drugs)} labels fetched")
    coverage = sum(1 for d in drugs if labels.get(d))
    print(f"  [load] label coverage: {coverage}/{len(drugs)}")

    return {
        "drugs": drugs,
        "drug_totals": drug_totals,
        "event_totals": event_totals,
        "total_reports": total_reports,
        "pair_counts": pair_counts,
        "labels": labels,
        "indications": indications,
        "event_api_last_updated": event_api_last_updated,
        "label_api_last_updated": label_api_last_updated,
    }


# ═══════════════════════════════════════════════════════════════
# ANALYSIS - domain-agnostic statistical method
# ═══════════════════════════════════════════════════════════════

def _build_table(data, min_cell_size):
    """Build rows: {drug, event, a, b, c, d, label}"""
    rows = []
    skipped_generic = 0
    skipped_small = 0
    skipped_nolabel = 0
    a_gt_ndrug = 0
    a_gt_nevent = 0
    min_a = float("inf")
    max_a = 0
    n_total = data["total_reports"]
    for (drug, event), a in data["pair_counts"].items():
        if event in GENERIC_TOKEN_SKIP:
            skipped_generic += 1
            continue
        if a < min_cell_size:
            skipped_small += 1
            continue
        n_drug = data["drug_totals"].get(drug, 0)
        n_event = data["event_totals"].get(event, 0)
        label_text = data["labels"].get(drug, "")
        if not label_text:
            skipped_nolabel += 1
            continue
        # Integrity: marginals and joint come from separate openFDA aggregations;
        # log the frequency of cross-query inconsistency rather than silently clamp.
        if a > n_drug:
            a_gt_ndrug += 1
        if a > n_event:
            a_gt_nevent += 1
        b = max(0, n_drug - a)
        c = max(0, n_event - a)
        d = max(0, n_total - a - b - c)
        is_positive = 1 if _is_event_in_label(event, label_text) else 0
        if a < min_a:
            min_a = a
        if a > max_a:
            max_a = a
        rows.append({
            "drug": drug,
            "event": event,
            "a": a, "b": b, "c": c, "d": d,
            "label": is_positive,
        })
    return rows, {
        "skipped_generic": skipped_generic,
        "skipped_small": skipped_small,
        "skipped_nolabel": skipped_nolabel,
        "rows_kept": len(rows),
        "a_exceeds_n_drug": a_gt_ndrug,
        "a_exceeds_n_event": a_gt_nevent,
        "min_joint_count_a": (min_a if min_a != float("inf") else 0),
        "max_joint_count_a": max_a,
    }


def _score_table(rows, prior):
    metrics = {"PRR": [], "ROR": [], "IC": [], "EBGM": []}
    labels = []
    for r in rows:
        a, b, c, d = r["a"], r["b"], r["c"], r["d"]
        metrics["PRR"].append(compute_prr(a, b, c, d))
        metrics["ROR"].append(compute_ror(a, b, c, d))
        metrics["IC"].append(compute_ic(a, b, c, d))
        metrics["EBGM"].append(compute_ebgm(a, b, c, d, prior))
        labels.append(r["label"])
    return metrics, labels


def run_analysis(data):
    rng = random.Random(RANDOM_SEED)
    results = {
        "analysis": "disproportionality_vs_label_mention",
        "config": {
            "n_drugs": N_DRUGS,
            "n_events_per_drug": N_EVENTS_PER_DRUG,
            "min_event_global_count": MIN_EVENT_GLOBAL_COUNT,
            "min_cell_size_main": MIN_CELL_SIZE_MAIN,
            "sensitivity_min_cell_sizes": SENSITIVITY_MIN_CELL_SIZES,
            "label_sections": LABEL_SECTIONS,
            "bootstrap_n": BOOTSTRAP_N,
            "permutation_n": PERMUTATION_N,
            "recall_thresholds": RECALL_THRESHOLDS,
            "random_seed": RANDOM_SEED,
        },
        "data_summary": {
            "total_reports": data["total_reports"],
            "n_drugs": len(data["drugs"]),
            "drugs": data["drugs"],
            "n_pairs_raw": len(data["pair_counts"]),
            "label_coverage": sum(1 for d in data["drugs"] if data["labels"].get(d)),
            "event_api_last_updated": data.get("event_api_last_updated", ""),
            "label_api_last_updated": data.get("label_api_last_updated", ""),
        },
    }

    # Fit the EBGM prior from all (drug, event) pairs using RR = a / E.
    rr_values = []
    n_total = data["total_reports"]
    for (drug, event), a in data["pair_counts"].items():
        n_drug = data["drug_totals"].get(drug, 0)
        n_event = data["event_totals"].get(event, 0)
        if n_drug == 0 or n_event == 0 or n_total == 0:
            continue
        b = max(0, n_drug - a)
        c = max(0, n_event - a)
        expected = (a + b) * (a + c) / n_total
        if expected > 0:
            rr_values.append(a / expected)
    prior = _fit_gamma_mom(rr_values)
    results["ebgm_prior"] = {"alpha": prior[0], "beta": prior[1],
                              "rr_mean": sum(rr_values) / len(rr_values) if rr_values else 0,
                              "n_pairs_for_prior": len(rr_values)}

    # Main benchmark at MIN_CELL_SIZE_MAIN.
    print(f"  [analysis] building main table at min_cell_size={MIN_CELL_SIZE_MAIN}")
    rows, skip_stats = _build_table(data, MIN_CELL_SIZE_MAIN)
    results["main_table_counts"] = {
        "n_rows": len(rows),
        "n_positive": sum(r["label"] for r in rows),
        "n_negative": sum(1 for r in rows if r["label"] == 0),
        **skip_stats,
    }
    print(f"  [analysis] {len(rows)} pairs, positive={sum(r['label'] for r in rows)}")

    if len(rows) < 50 or sum(r["label"] for r in rows) < 5:
        raise RuntimeError("Insufficient data for benchmark — check API or filters")

    metric_scores, labels = _score_table(rows, prior)

    # AUC and PPV per metric (point estimate + bootstrap CI + permutation p).
    results["main_benchmark"] = {}
    main_pvals = []
    main_pval_names = []
    for name, scores in metric_scores.items():
        point, lo, hi = bootstrap_auc_ci(scores, labels, BOOTSTRAP_N, rng)
        pval = permutation_auc_pvalue(scores, labels, PERMUTATION_N, rng, two_sided=True)
        ppv = {}
        for r in RECALL_THRESHOLDS:
            ppv[f"ppv_at_recall_{int(r*100)}pct"] = ppv_at_recall(scores, labels, r)
        results["main_benchmark"][name] = {
            "auc": point,
            "auc_95ci": [lo, hi],
            "permutation_p_two_sided": pval,
            "ppv_at_recall": ppv,
            "n_rows": len(rows),
        }
        main_pvals.append(pval)
        main_pval_names.append(name)
        print(f"  [analysis] {name}: AUC={point:.4f} [{lo:.4f}, {hi:.4f}] p_perm_2s={pval:.4g}")

    # BH-FDR correction across the four main-benchmark permutation p-values.
    if main_pvals:
        qvals = benjamini_hochberg(main_pvals)
        for name, q in zip(main_pval_names, qvals):
            results["main_benchmark"][name]["permutation_q_bh"] = q

    # Paired AUC differences between metrics.
    print("  [analysis] pairwise paired-bootstrap AUC differences")
    diff_results = {}
    metric_names = list(metric_scores.keys())
    for i in range(len(metric_names)):
        for j in range(i + 1, len(metric_names)):
            a_name, b_name = metric_names[i], metric_names[j]
            obs, lo, hi = paired_bootstrap_auc_diff(
                metric_scores[a_name], metric_scores[b_name], labels, BOOTSTRAP_N, rng
            )
            diff_results[f"{a_name}_minus_{b_name}"] = {
                "diff": obs, "ci_95": [lo, hi],
            }
    results["pairwise_auc_differences"] = diff_results

    # Spearman correlation between metric score vectors.
    print("  [analysis] Spearman correlation between metric scores")
    corr_matrix = {}
    for i in range(len(metric_names)):
        for j in range(i + 1, len(metric_names)):
            a_name, b_name = metric_names[i], metric_names[j]
            rho = spearman_rho(metric_scores[a_name], metric_scores[b_name])
            corr_matrix[f"{a_name}_vs_{b_name}"] = rho
    results["metric_score_spearman"] = corr_matrix

    # Sensitivity sweep over min cell size. Bootstrap 95% CIs at each setting.
    print("  [analysis] sensitivity sweep over min_cell_size")
    sweep = {}
    for mcs in SENSITIVITY_MIN_CELL_SIZES:
        rows_mcs, _skip = _build_table(data, mcs)
        if len(rows_mcs) < 30 or sum(r["label"] for r in rows_mcs) < 3:
            sweep[str(mcs)] = {"n_rows": len(rows_mcs), "n_positive": sum(r["label"] for r in rows_mcs),
                               "note": "skipped — insufficient data"}
            continue
        scores_sw, labels_sw = _score_table(rows_mcs, prior)
        entry = {
            "n_rows": len(rows_mcs),
            "n_positive": sum(r["label"] for r in rows_mcs),
            "auc": {},
            "auc_95ci": {},
        }
        for name, scores in scores_sw.items():
            pt, lo, hi = bootstrap_auc_ci(scores, labels_sw, BOOTSTRAP_N, rng)
            entry["auc"][name] = pt
            entry["auc_95ci"][name] = [lo, hi]
        sweep[str(mcs)] = entry
    results["sensitivity_min_cell_size"] = sweep

    # Sensitivity sweep over label-section set.
    print("  [analysis] sensitivity sweep over label section set")
    section_sweep = {}
    section_configs = {
        "adverse_reactions_only": ["adverse_reactions"],
        "warnings_only": ["warnings_and_cautions"],
        "boxed_warning_only": ["boxed_warning"],
        "ar_plus_warnings": ["adverse_reactions", "warnings_and_cautions"],
        "all_three": LABEL_SECTIONS,
    }
    # Pre-extract label text under each config from cache already loaded.
    # Because load_data normalized with LABEL_SECTIONS only, we re-read raw labels here.
    raw_label_cache = {}
    for drug in data["drugs"]:
        fname = os.path.join(CACHE_DIR, f"label_{_safe_filename(drug)}.json")
        try:
            with open(fname, "r") as fh:
                payload = json.load(fh)
            recs = payload.get("results", [])
            raw_label_cache[drug] = recs[0] if recs else None
        except (FileNotFoundError, json.JSONDecodeError):
            raw_label_cache[drug] = None

    for cfg_name, sections in section_configs.items():
        data_cfg = {
            **data,
            "labels": {
                d: _extract_label_text(raw_label_cache[d], sections) if raw_label_cache.get(d) else ""
                for d in data["drugs"]
            },
        }
        rows_cfg, _ = _build_table(data_cfg, MIN_CELL_SIZE_MAIN)
        if len(rows_cfg) < 30 or sum(r["label"] for r in rows_cfg) < 3:
            section_sweep[cfg_name] = {"n_rows": len(rows_cfg),
                                         "n_positive": sum(r["label"] for r in rows_cfg),
                                         "note": "skipped — insufficient data"}
            continue
        scores_cfg, labels_cfg = _score_table(rows_cfg, prior)
        entry = {
            "n_rows": len(rows_cfg),
            "n_positive": sum(r["label"] for r in rows_cfg),
            "auc": {},
            "auc_95ci": {},
            "permutation_p_two_sided": {},
        }
        for n in metric_names:
            pt, lo, hi = bootstrap_auc_ci(scores_cfg[n], labels_cfg, BOOTSTRAP_N, rng)
            entry["auc"][n] = pt
            entry["auc_95ci"][n] = [lo, hi]
            entry["permutation_p_two_sided"][n] = permutation_auc_pvalue(
                scores_cfg[n], labels_cfg, PERMUTATION_N, rng, two_sided=True)
        section_sweep[cfg_name] = entry
    results["sensitivity_label_section"] = section_sweep

    # Sensitivity: exclude pairs where the event appears in indications_and_usage.
    # These are often the drug's indication rather than an adverse reaction, and
    # may dominate disproportionate-reporting signals, inverting AUCs.
    print("  [analysis] sensitivity: exclude pairs overlapping indications_and_usage")
    indication_overlap_sweep = {}
    indication_texts = data.get("indications", {})

    # Two configs: (a) drop overlapping pairs entirely; (b) drop them only from
    # the negative class (overlap pairs are likely mislabeled positives).
    def _build_filtered(data, min_cell_size, mode):
        rows_local, _ = _build_table(data, min_cell_size)
        filtered = []
        for r in rows_local:
            ind_text = indication_texts.get(r["drug"], "")
            is_ind = _is_event_in_label(r["event"], ind_text) if ind_text else False
            if mode == "drop_overlap" and is_ind:
                continue
            if mode == "drop_overlap_negatives_only" and is_ind and r["label"] == 0:
                continue
            filtered.append(r)
        return filtered

    for mode in ("drop_overlap", "drop_overlap_negatives_only"):
        rows_f = _build_filtered(data, MIN_CELL_SIZE_MAIN, mode)
        if len(rows_f) < 30 or sum(r["label"] for r in rows_f) < 3:
            indication_overlap_sweep[mode] = {"n_rows": len(rows_f),
                                                "n_positive": sum(r["label"] for r in rows_f),
                                                "note": "skipped — insufficient data"}
            continue
        scores_f, labels_f = _score_table(rows_f, prior)
        entry = {
            "n_rows": len(rows_f),
            "n_positive": sum(r["label"] for r in rows_f),
            "auc": {},
            "auc_95ci": {},
        }
        for name in metric_names:
            pt, lo, hi = bootstrap_auc_ci(scores_f[name], labels_f, BOOTSTRAP_N, rng)
            entry["auc"][name] = pt
            entry["auc_95ci"][name] = [lo, hi]
        indication_overlap_sweep[mode] = entry
    results["sensitivity_indication_exclusion"] = indication_overlap_sweep

    # Negative control / falsification: if metrics carry label information, a single
    # shuffle should destroy it (AUC -> ~0.5). Average over NEG_CTRL_ITERS shuffles
    # using the seeded RNG so we get a stable estimate. We assert that the mean
    # |AUC_shuffled - 0.5| < 0.05 in verification.
    print("  [analysis] negative-control / falsification check (shuffled labels)")
    neg_ctrl = {}
    NEG_CTRL_ITERS = 50
    shuffled_labels = list(labels)
    for name, scores in metric_scores.items():
        aucs = []
        for _ in range(NEG_CTRL_ITERS):
            rng.shuffle(shuffled_labels)
            aucs.append(auc_roc(scores, shuffled_labels))
        mean_auc = sum(aucs) / len(aucs)
        max_dev = max(abs(a - 0.5) for a in aucs)
        neg_ctrl[name] = {
            "mean_shuffled_auc": mean_auc,
            "mean_abs_deviation_from_0_5": abs(mean_auc - 0.5),
            "max_abs_deviation_from_0_5": max_dev,
            "n_shuffles": NEG_CTRL_ITERS,
        }
        print(f"  [analysis] {name} neg-ctrl: mean shuffled AUC={mean_auc:.4f} "
              f"(|AUC-0.5| mean={abs(mean_auc-0.5):.4f}, max={max_dev:.4f})")
    results["negative_control"] = neg_ctrl

    # Robustness: report whether the ranking of metrics by AUC is preserved across
    # the primary and each sensitivity configuration. We rank metrics by AUC in
    # each configuration and compute the Spearman rank correlation with the main
    # ranking; values near 1.0 indicate robustness.
    main_ranked = sorted(metric_names, key=lambda n: -results["main_benchmark"][n]["auc"])
    main_rank_map = {n: i for i, n in enumerate(main_ranked)}
    ranking_consistency = {}
    def _rank_map(aucs):
        order = sorted(metric_names, key=lambda n: -aucs[n])
        return {n: i for i, n in enumerate(order)}
    for mcs, entry in results["sensitivity_min_cell_size"].items():
        if "auc" in entry:
            rm = _rank_map(entry["auc"])
            # Kendall-ish agreement: count pairs with same relative order.
            agree = 0
            tot = 0
            for i in range(len(metric_names)):
                for j in range(i + 1, len(metric_names)):
                    tot += 1
                    a, b = metric_names[i], metric_names[j]
                    if (main_rank_map[a] < main_rank_map[b]) == (rm[a] < rm[b]):
                        agree += 1
            ranking_consistency[f"min_cell_size={mcs}"] = agree / tot if tot else 0.0
    results["ranking_consistency_vs_main"] = ranking_consistency

    # Explicit per-result limitations block, echoed in stdout and in report.md so
    # any downstream reader is confronted with the key caveats alongside numbers.
    results["limitations"] = [
        "Ground truth: label mention is an imperfect proxy for 'known' ADRs — labels "
        "omit many recognised reactions and include boilerplate, and FAERS is affected "
        "by confounding, notoriety bias, and indication bias (so AUC < 0.5 is plausible).",
        "FAERS bias: spontaneous-reporting data is subject to under-reporting, "
        "stimulated reporting after media attention, and selection bias (reports tend to "
        "come from patients who suffered events).",
        "Indication confounding: events that are indications for the drug (reasons for "
        "treatment) appear disproportionately in reports but are not labelled as adverse "
        "reactions, which inflates false positives; we quantify this in the indication-"
        "exclusion sensitivity analysis but cannot eliminate it without a curated set.",
        "openFDA pagination: the per-drug event query is truncated at "
        f"N_EVENTS_PER_DRUG={N_EVENTS_PER_DRUG}, so low-count events are missing and the "
        "min-cell-size sensitivity sweep cannot probe very small joint counts.",
        "Marginal consistency: drug and event marginals come from separate openFDA "
        "count aggregations; small cross-query inconsistencies (logged as a_exceeds_* "
        "fields) are possible when openFDA's index updates between queries.",
        "The analysis does not adjust for drug-event co-medication, temporal trends, "
        "or age/sex stratification; the 2x2 table collapses these into single counts.",
    ]

    return results


# ═══════════════════════════════════════════════════════════════
# REPORT - results.json + report.md
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open(RESULTS_FILE, "w") as fh:
        json.dump(results, fh, indent=2, default=str)

    lines = []
    lines.append("# Disproportionality Metric Benchmark Against FDA Label Mentions\n")
    ds = results["data_summary"]
    cfg = results["config"]
    lines.append(f"**Total FAERS reports:** {ds['total_reports']:,}\n")
    lines.append(f"**Drugs studied:** {ds['n_drugs']}, "
                 f"label coverage: {ds['label_coverage']}/{ds['n_drugs']}\n")
    lines.append(f"**Raw drug-event pairs:** {ds['n_pairs_raw']:,}\n")

    mt = results["main_table_counts"]
    lines.append(f"\n**Main table at min_cell_size={cfg['min_cell_size_main']}:** "
                 f"{mt['n_rows']:,} pairs — {mt['n_positive']} positive, "
                 f"{mt['n_negative']} negative "
                 f"(skipped: generic={mt['skipped_generic']}, "
                 f"small={mt['skipped_small']}, no-label-coverage={mt['skipped_nolabel']})\n")

    prior = results["ebgm_prior"]
    lines.append(f"\n**EBGM gamma prior (method of moments):** alpha={prior['alpha']:.4g}, "
                 f"beta={prior['beta']:.4g}, fit on {prior['n_pairs_for_prior']:,} pairs.\n")

    lines.append("\n## Main benchmark\n")
    lines.append("| Metric | AUC | 95% CI | Perm p (2-sided) | BH q | PPV@1% | PPV@5% | PPV@10% | PPV@20% |")
    lines.append("|--------|----:|:------:|:----------------:|:----:|------:|------:|-------:|-------:|")
    for name, entry in results["main_benchmark"].items():
        lo, hi = entry["auc_95ci"]
        ppv = entry["ppv_at_recall"]
        pval = entry["permutation_p_two_sided"]
        qval = entry.get("permutation_q_bh", float("nan"))
        lines.append(
            f"| {name} | {entry['auc']:.4f} | [{lo:.4f}, {hi:.4f}] | "
            f"{pval:.4g} | {qval:.4g} | "
            f"{ppv['ppv_at_recall_1pct']:.4f} | {ppv['ppv_at_recall_5pct']:.4f} | "
            f"{ppv['ppv_at_recall_10pct']:.4f} | {ppv['ppv_at_recall_20pct']:.4f} |"
        )

    lines.append("\n## Pairwise paired-bootstrap AUC differences\n")
    lines.append("| Contrast | Difference | 95% CI |")
    lines.append("|----------|-----------:|:------:|")
    for name, entry in results["pairwise_auc_differences"].items():
        lo, hi = entry["ci_95"]
        lines.append(f"| {name} | {entry['diff']:+.4f} | [{lo:+.4f}, {hi:+.4f}] |")

    lines.append("\n## Metric score Spearman correlations\n")
    lines.append("| Pair | Spearman rho |")
    lines.append("|------|-------------:|")
    for name, rho in results["metric_score_spearman"].items():
        lines.append(f"| {name} | {rho:.4f} |")

    lines.append("\n## Sensitivity: minimum cell size\n")
    lines.append("| min_cell_size | n_rows | n_positive | AUC PRR | AUC ROR | AUC IC | AUC EBGM |")
    lines.append("|--------------:|------:|-----------:|-------:|-------:|------:|--------:|")
    for mcs, entry in results["sensitivity_min_cell_size"].items():
        if "note" in entry:
            lines.append(f"| {mcs} | {entry['n_rows']} | {entry['n_positive']} | — | — | — | — |")
        else:
            a = entry["auc"]
            lines.append(
                f"| {mcs} | {entry['n_rows']} | {entry['n_positive']} | "
                f"{a['PRR']:.4f} | {a['ROR']:.4f} | {a['IC']:.4f} | {a['EBGM']:.4f} |"
            )

    lines.append("\n## Sensitivity: label section set\n")
    lines.append("| Sections | n_rows | n_positive | AUC PRR [95% CI] | AUC ROR [95% CI] | AUC IC [95% CI] | AUC EBGM [95% CI] |")
    lines.append("|----------|-------:|-----------:|:----------------:|:----------------:|:---------------:|:-----------------:|")
    for cfg_name, entry in results["sensitivity_label_section"].items():
        if "note" in entry:
            lines.append(f"| {cfg_name} | {entry['n_rows']} | {entry['n_positive']} | — | — | — | — |")
        else:
            a = entry["auc"]
            ci = entry.get("auc_95ci", {})
            def _fmt(name):
                c = ci.get(name, [None, None])
                return (f"{a[name]:.3f} [{c[0]:.3f},{c[1]:.3f}]"
                        if c[0] is not None and c[1] is not None else f"{a[name]:.3f}")
            lines.append(
                f"| {cfg_name} | {entry['n_rows']} | {entry['n_positive']} | "
                f"{_fmt('PRR')} | {_fmt('ROR')} | {_fmt('IC')} | {_fmt('EBGM')} |"
            )

    if "sensitivity_indication_exclusion" in results:
        lines.append("\n## Sensitivity: indication overlap exclusion\n")
        lines.append("| Mode | n_rows | n_positive | AUC PRR [95% CI] | AUC ROR [95% CI] | AUC IC [95% CI] | AUC EBGM [95% CI] |")
        lines.append("|------|-------:|-----------:|:----------------:|:----------------:|:---------------:|:-----------------:|")
        for mode, entry in results["sensitivity_indication_exclusion"].items():
            if "note" in entry:
                lines.append(f"| {mode} | {entry['n_rows']} | {entry['n_positive']} | — | — | — | — |")
            else:
                a = entry["auc"]; ci = entry.get("auc_95ci", {})
                def _fmt2(name):
                    c = ci.get(name, [None, None])
                    return (f"{a[name]:.3f} [{c[0]:.3f},{c[1]:.3f}]"
                            if c[0] is not None and c[1] is not None else f"{a[name]:.3f}")
                lines.append(
                    f"| {mode} | {entry['n_rows']} | {entry['n_positive']} | "
                    f"{_fmt2('PRR')} | {_fmt2('ROR')} | {_fmt2('IC')} | {_fmt2('EBGM')} |"
                )

    if "negative_control" in results:
        lines.append("\n## Negative control (shuffled labels)\n")
        lines.append("Sanity check — under repeatedly shuffled labels, AUC should collapse "
                     "to 0.5 for every metric.\n")
        lines.append("| Metric | mean shuffled AUC | mean |AUC-0.5| | max |AUC-0.5| |")
        lines.append("|--------|------------------:|---------------:|---------------:|")
        for name, e in results["negative_control"].items():
            lines.append(f"| {name} | {e['mean_shuffled_auc']:.4f} | "
                         f"{e['mean_abs_deviation_from_0_5']:.4f} | "
                         f"{e['max_abs_deviation_from_0_5']:.4f} |")

    if "ranking_consistency_vs_main" in results:
        lines.append("\n## Metric ranking consistency across min-cell-size sweep\n")
        lines.append("Pairwise agreement with the main-config ranking (1.0 = identical).\n")
        lines.append("| Config | Pairwise agreement |")
        lines.append("|--------|------------------:|")
        for cfg, v in results["ranking_consistency_vs_main"].items():
            lines.append(f"| {cfg} | {v:.3f} |")

    if "limitations" in results:
        lines.append("\n## Limitations and failure modes\n")
        for i, lim in enumerate(results["limitations"], start=1):
            lines.append(f"{i}. {lim}")

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


# ═══════════════════════════════════════════════════════════════
# VERIFICATION
# ═══════════════════════════════════════════════════════════════

def verify_results():
    print("\n=== VERIFICATION MODE ===\n")
    if not os.path.exists(RESULTS_FILE):
        print("FAIL: results.json not found")
        sys.exit(1)
    with open(RESULTS_FILE, "r") as fh:
        results = json.load(fh)

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

    ds = results.get("data_summary", {})
    check("total FAERS reports >= 1,000,000",
          ds.get("total_reports", 0) >= 1_000_000)
    check("studied at least 20 drugs",
          ds.get("n_drugs", 0) >= 20)
    check("label coverage >= 60% of drugs",
          ds.get("label_coverage", 0) >= 0.6 * ds.get("n_drugs", 1))

    mt = results.get("main_table_counts", {})
    check("main table has at least 200 drug-event pairs",
          mt.get("n_rows", 0) >= 200)
    check("main table has both positive and negative examples",
          mt.get("n_positive", 0) > 10 and mt.get("n_negative", 0) > 10)

    mb = results.get("main_benchmark", {})
    check("all four metrics benchmarked (PRR, ROR, IC, EBGM)",
          set(mb.keys()) == {"PRR", "ROR", "IC", "EBGM"})
    for metric in ("PRR", "ROR", "IC", "EBGM"):
        entry = mb.get(metric, {})
        check(f"{metric} AUC in (0, 1)",
              0 < entry.get("auc", -1) < 1)
        check(f"{metric} AUC 95% CI has finite bounds",
              (entry.get("auc_95ci", [None, None])[0] is not None
               and entry["auc_95ci"][1] is not None
               and entry["auc_95ci"][0] <= entry["auc"] <= entry["auc_95ci"][1] + 1e-9))
    check("all main permutation p-values computed (two-sided)",
          all("permutation_p_two_sided" in mb.get(m, {}) for m in ("PRR", "ROR", "IC", "EBGM")))

    # PRR and ROR are mathematically closely related — their rank correlation should be near 1.
    corr = results.get("metric_score_spearman", {}).get("PRR_vs_ROR", 0)
    check("PRR and ROR Spearman >= 0.95 (sanity — same family)",
          corr >= 0.95)

    sweep = results.get("sensitivity_min_cell_size", {})
    check("sensitivity sweep includes at least 4 thresholds",
          sum(1 for v in sweep.values() if "auc" in v) >= 4)

    sec_sweep = results.get("sensitivity_label_section", {})
    check("label-section sweep includes at least 3 configs with AUCs",
          sum(1 for v in sec_sweep.values() if "auc" in v) >= 3)

    diffs = results.get("pairwise_auc_differences", {})
    check("all 6 pairwise AUC differences present",
          len(diffs) == 6)

    ind_sweep = results.get("sensitivity_indication_exclusion", {})
    check("indication-exclusion sensitivity computed for both modes",
          sum(1 for v in ind_sweep.values() if "auc" in v) >= 2)

    # At least one of the label-section sensitivity configs should have an AUC
    # that substantively departs from 0.5 in either direction (evidence that
    # disproportionality carries SOME information about label content).
    sec_aucs = []
    for v in sec_sweep.values():
        if "auc" in v:
            sec_aucs.extend(v["auc"].values())
    if sec_aucs:
        check("at least one label-section AUC |0.5 - AUC| >= 0.05",
              any(abs(a - 0.5) >= 0.05 for a in sec_aucs))
    else:
        check("at least one label-section AUC |0.5 - AUC| >= 0.05", False)

    check("report.md exists and > 800 bytes",
          os.path.exists(REPORT_FILE) and os.path.getsize(REPORT_FILE) > 800)

    # Negative control / falsification: shuffled-label AUC should be near 0.5.
    # Effect-size plausibility: observed AUC deviations should be bounded (<= 0.2),
    # i.e. no unrealistic "near-perfect" AUCs given the noisy ground truth.
    neg = results.get("negative_control", {})
    check("negative control computed for all four metrics",
          set(neg.keys()) == {"PRR", "ROR", "IC", "EBGM"})
    if neg:
        max_dev = max(e.get("mean_abs_deviation_from_0_5", 1.0) for e in neg.values())
        check("negative control mean |AUC - 0.5| < 0.05 for all metrics",
              max_dev < 0.05)

    # Effect-size / CI plausibility sanity checks.
    mb = results.get("main_benchmark", {})
    for metric in ("PRR", "ROR", "IC", "EBGM"):
        entry = mb.get(metric, {})
        auc = entry.get("auc", 0.5)
        check(f"{metric} |AUC - 0.5| <= 0.5 (bounded effect size)",
              abs(auc - 0.5) <= 0.5)
        lo, hi = entry.get("auc_95ci", [None, None])
        if lo is not None and hi is not None:
            check(f"{metric} CI width > 0.001 (reality check, not zero-width)",
                  (hi - lo) > 0.001)
            check(f"{metric} CI width < 0.5 (not uselessly wide)",
                  (hi - lo) < 0.5)

    # Ranking consistency: the metric ordering should be reasonably stable across
    # min-cell-size configurations (pharmacovigilance metrics are mathematically
    # close cousins, so large reshuffles would indicate numerical instability).
    rc = results.get("ranking_consistency_vs_main", {})
    if rc:
        mean_rc = sum(rc.values()) / len(rc)
        check("mean pairwise ranking agreement across min-cell-size sweep >= 0.6",
              mean_rc >= 0.6)

    # Key scientific finding from main benchmark: at least one metric has a CI
    # that does NOT cover 0.5 (i.e. we can distinguish from chance on at least
    # one metric). This is a falsifiable claim tied to our Research Question.
    mb = results.get("main_benchmark", {})
    distinguishable = any(
        (m.get("auc_95ci", [0.5, 0.5])[1] < 0.5
         or m.get("auc_95ci", [0.5, 0.5])[0] > 0.5)
        for m in mb.values()
    )
    check("main benchmark: at least one metric's 95% CI excludes 0.5",
          distinguishable)

    # Limitations block present.
    check("results.json contains a limitations block with >= 4 items",
          isinstance(results.get("limitations", []), list)
          and len(results.get("limitations", [])) >= 4)

    # Determinism sanity: every main-benchmark AUC 95% CI contains the point.
    all_contained = all(
        e.get("auc_95ci", [None, None])[0] is not None
        and e["auc_95ci"][0] - 1e-6 <= e["auc"] <= e["auc_95ci"][1] + 1e-6
        for e in mb.values()
    )
    check("all main benchmark point AUCs lie within their bootstrap CIs",
          all_contained)

    print(f"\n  Results: {passed}/{total} checks passed")
    if passed != total:
        print("  VERIFICATION FAILED")
        sys.exit(1)
    print("  VERIFICATION PASSED")
    print("  ALL CHECKS PASSED")


# ═══════════════════════════════════════════════════════════════
# MAIN
# ═══════════════════════════════════════════════════════════════

def main():
    if "--verify" in sys.argv:
        verify_results()
        return
    os.makedirs(CACHE_DIR, exist_ok=True)

    print("[1/4] Loading openFDA data (drugs, events, labels)")
    try:
        data = load_data()
    except RuntimeError as e:
        print(f"ERROR: data load failed — openFDA unreachable after retries: {e}",
              file=sys.stderr)
        print("This analysis REQUIRES network access to https://api.fda.gov. "
              "Retry when connectivity is restored.", file=sys.stderr)
        sys.exit(2)
    except Exception as e:
        print(f"ERROR: unexpected failure during data load: {e.__class__.__name__}: {e}",
              file=sys.stderr)
        sys.exit(3)

    print("[2/4] Running analysis")
    try:
        results = run_analysis(data)
    except RuntimeError as e:
        print(f"ERROR: analysis aborted — {e}", file=sys.stderr)
        print("Most likely cause: insufficient data after filters. "
              "Consider reducing MIN_CELL_SIZE_MAIN or MIN_EVENT_GLOBAL_COUNT.",
              file=sys.stderr)
        sys.exit(4)

    print("[3/4] Writing outputs")
    results["generated_at_utc"] = datetime.datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ")
    generate_report(results)
    print(f"  wrote {RESULTS_FILE}")
    print(f"  wrote {REPORT_FILE}")

    print("[4/4] Done")
    print("\n--- Known limitations and failure modes ---")
    for i, lim in enumerate(results.get("limitations", []), start=1):
        print(f"  ({i}) {lim}")
    print("-------------------------------------------\n")
    print("ANALYSIS COMPLETE")


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

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

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g && python3 analyze.py
```

**Expected output (values are approximate; exact numbers depend on openFDA snapshot):**

```
[1/4] Loading openFDA data (drugs, events, labels)
  [load] fetching top drugs by report volume
  [load] got 50 drugs
  [load] fetching global event marginals
  [load] kept ~900-999 events with count >= 100
  [load] fetching total report count
  [load] total FAERS reports: ~20,000,000
  [load] fetching per-drug top events (N=50 drugs)
    [load] 10/50 drugs fetched
    ...
  [load] got ~4000-8000 drug-event pairs
  [load] fetching labels for each drug
    [load] 10/50 drugs fetched
    ...
  [load] label coverage: ~30-50/50
[2/4] Running analysis
  [analysis] building main table at min_cell_size=3
  [analysis] ~2000-5000 pairs, positive=~500-1500
  [analysis] PRR: AUC=0.35-0.50 [lo, hi] p_perm<0.01 (two-sided)
  [analysis] ROR: AUC=0.35-0.50 [lo, hi] p_perm<0.01 (two-sided)
  [analysis] IC:  AUC=0.35-0.50 [lo, hi] p_perm<0.01 (two-sided)
  [analysis] EBGM:AUC=0.35-0.50 [lo, hi] p_perm<0.01 (two-sided)
  [analysis] pairwise paired-bootstrap AUC differences
  [analysis] Spearman correlation between metric scores
  [analysis] sensitivity sweep over min_cell_size
  [analysis] sensitivity sweep over label section set
  [analysis] sensitivity: exclude pairs overlapping indications_and_usage
  [analysis] negative-control / falsification check (shuffled labels)
  [analysis] PRR neg-ctrl: mean shuffled AUC~0.500 (|AUC-0.5| mean<0.01, max<0.05)
  ...
[3/4] Writing outputs
  wrote /tmp/.../results.json
  wrote /tmp/.../report.md
[4/4] Done
--- Known limitations and failure modes ---
  (1) Ground truth: label mention is an imperfect proxy ...
  ...
-------------------------------------------

ANALYSIS COMPLETE
```

**Expected files created:**
- `/tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g/results.json`
- `/tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g/report.md`
- `/tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g/cache/*.json` + `*.sha256`

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_disproportionality-metric-comparison-with-label-changes-as-g && python3 analyze.py --verify
```

**Expected output:**

```
=== VERIFICATION MODE ===

  [1] PASS: total FAERS reports >= 1,000,000
  [2] PASS: studied at least 20 drugs
  ... (~36 checks total covering data volume, per-metric AUC/CI sanity,
       Spearman consistency, sensitivity coverage, negative control,
       effect-size plausibility, CI-width bounds, and ranking consistency)

  Results: 36/36 checks passed
  VERIFICATION PASSED
  ALL CHECKS PASSED
```

## Success Criteria

1. All verification checks pass (expect >= 30 checks after improvements).
2. `results.json` contains a benchmark for four metrics (PRR, ROR, IC, EBGM) with AUCs,
   bootstrap 95% CIs, and two-sided permutation p-values.
3. `report.md` reproduces the benchmark, sensitivity, negative-control, ranking-
   consistency and limitations tables/sections.
4. The analysis uses >= 20 drugs and >= 200 drug-event pairs.
5. Sensitivity analyses (cell size, label section, indication exclusion) are populated with
   at least two configurations each.
6. At least one label-section configuration shows |AUC - 0.5| >= 0.05, i.e. metrics carry
   information about label content in at least one operational definition of the gold
   standard.
7. Negative-control shuffled-label AUC is within 0.05 of 0.5 for every metric
   (falsification check passes: metric discrimination is driven by the labels, not by
   incidental ordering).
8. Main-benchmark effect-size plausibility: every |AUC - 0.5| <= 0.5, every
   bootstrap CI width is in (0.001, 0.5), and every point AUC lies within its own CI.
9. Script stdout terminates with `ANALYSIS COMPLETE`; `--verify` stdout terminates with
   `ALL CHECKS PASSED`.

## Failure Conditions

1. Any verification check fails.
2. openFDA is unreachable after retries (network error or rate limit): script exits with
   non-zero status and a clear stderr message, instead of producing corrupt output.
3. Fewer than 20 drugs return usable event counts.
4. Label coverage drops below 60% of studied drugs (unexpected openFDA data gap).
5. Script requires any non-stdlib dependency or interactive input.
6. Negative-control mean |AUC - 0.5| >= 0.05 for any metric (indicates a bug in the
   label-permutation loop or broken AUC computation).
7. A main-benchmark bootstrap CI has zero width or exceeds 0.5 width (indicates
   degenerate bootstrap or severely under-powered sample).