examples: add the IMPROVE neoepitope immunogenicity cohort

.gitignore

@@ -28,3 +28,9 @@ examples/gusto/out/
 validation/predrupdate/gusto_y_p.csv
 validation/predrupdate/synpm_estimates.csv
 validation/predrupdate/gusto_estimates.csv
+
+# Same reasoning for the IMPROVE neoepitope fixture: build_fixture.py regenerates both parquets
+# deterministically from the released IMPROVE_paper results archive.
+examples/improve/cohort.parquet
+examples/improve/predictions.parquet
+examples/improve/out/

examples/improve/README.md

@@ -0,0 +1,120 @@
+# IMPROVE neoepitope immunogenicity — predval example cohort
+
+This directory applies predval to a problem outside the clinical-risk-score setting it was
+built in: deciding which candidate neoepitopes from a patient's tumour are worth putting into a
+personalised vaccine or a TCR-T product.
+
+## Why this cohort
+
+Every published neoepitope immunogenicity predictor reports **ranking** metrics — AUROC, AUPRC,
+top-N recall, fraction-ranked. A literature check (2026-10-04) found no published work reporting
+**calibration** for peptide-level immunogenicity predictors: whether a peptide scored 0.8
+validates 80% of the time. The nearest work is adjacent rather than overlapping and is cited
+below.
+
+Calibration is not a cosmetic concern here, because selection is budget-constrained. A
+personalised vaccine encodes on the order of 20 neoepitopes; a TCR-T product carries 3-5. Those
+slots are the scarce resource, and a probability that over-promises by a large factor misprices
+every one of them.
+
+## Provenance
+
+Source: Borch et al. 2024, *IMPROVE: a feature model to predict neoepitope immunogenicity
+through broad-scale validation of T cell recognition*, Frontiers in Immunology. Data release:
+[github.com/SRHgroup/IMPROVE_paper](https://github.com/SRHgroup/IMPROVE_paper), `results.zip`.
+
+- 17,520 candidate neoepitopes, 70 patients, three accrual cohorts (bladder 6,237 peptides / 24
+  patients; melanoma 5,921 / 26; a mixed "Basket" group 5,362 / 20).
+- 467 peptides are experimentally immunogenic — **2.666% prevalence**.
+- Predictions are the paper's own **out-of-fold** scores from its 5-fold cross-validation. That
+  CV was partitioned **by patient**: every patient sits in exactly one partition, verified in
+  `build_fixture.py`. No score here comes from a model that had seen another peptide from the
+  same patient, which is what makes the file admissible as a predictions-only validation input
+  rather than an apparent-fit artefact.
+- Seven models: the IMPROVE random forest in its three released feature configurations
+  (`Simple`; `TME_excluded`, adding cellular prevalence and a prioritisation score;
+  `TME_included`, further adding cytolytic activity, HLA expression and MCP-counter estimates),
+  each in full and in the `wo_prime` ablation that drops the PRIME feature, plus the paper's
+  `nnalign` comparator.
+
+Nothing in this directory re-trains, re-tunes or re-ranks anything. predval reads probabilities
+and labels and nothing else.
+
+## Unit of analysis
+
+The cluster is the **patient**, not the peptide. Peptides from one patient share that patient's
+HLA genotype, mutational process, tumour transcriptome and microenvironment — the dominant
+sources of variance in every feature these models consume. Mean cluster size is 250 peptides,
+an order of magnitude larger than the PCam fixture's slides, so the cost of getting this wrong
+is correspondingly larger. Every interval in the report resamples whole patients.
+
+## Prespecified subgroups
+
+Both were declared before any metric was computed on this cohort.
+
+- `tumour_cohort` — the release's own three accrual groups. A model that holds in melanoma, the
+  indication these predictors were overwhelmingly trained and tuned on, and fails elsewhere is
+  precisely the failure mode that matters for carrying a verdict to another disease. This is the
+  closest available proxy for that transfer question, and its result bounds how far anything
+  here may be carried. MSS colorectal cancer is not in this cohort at all.
+- `peptide_length` — 9-mers versus everything else. The binding predictors that supply the
+  dominant features are trained overwhelmingly on 9-mers and are known to degrade off that
+  length, so the off-length stratum is where a calibration fault would be expected *a priori*.
+
+## Reproducing
+
+```sh
+# 1. fetch the release and extract results.zip
+git clone https://github.com/SRHgroup/IMPROVE_paper.git
+unzip -q IMPROVE_paper/results.zip -x '__MACOSX/*' -d improve_results
+
+# 2. build the fixture (writes cohort.parquet and predictions.parquet here)
+uv run python examples/improve/build_fixture.py --src 'improve_results/results copy'
+
+# 3. the main evaluation -- B=2000 patient-clustered bootstrap, several minutes
+uv run python examples/improve/run_evaluation.py
+
+# 4. the budgeted-yield sidecar
+uv run python examples/improve/budget_yield.py
+```
+
+`run_evaluation.py` writes `out/report.html`, `out/findings.json` and the parquet artefacts.
+`budget_yield.py` writes `out/budget_yield.parquet`.
+
+## The budgeted-yield sidecar
+
+`budget_yield.py` is deliberately *not* part of predval's contract. predval reports calibration;
+the mapping from calibration to a decision is cohort-specific, and this is that mapping for this
+cohort. For each model and each budget k ∈ {3, 5, 10, 20} it takes every patient's own top-k
+candidates by score and reports, per patient and averaged with a patient-level bootstrap:
+
+- `observed` — immunogenic peptides actually in that top-k;
+- `promised` — the sum of the model's published probabilities over the *same* top-k set, which is
+  what those probabilities assert the yield will be;
+- `recalibrated` — the same sum after predval's cross-fitted rung1 (intercept-only) correction,
+  fitted on other patients.
+
+Because all three columns score the identical top-k set, `promised − observed` isolates the
+miscalibration cost in units of antigens. It is not a ranking error.
+
+## Prior art this sits next to
+
+- **ImmUQBench** (Qayyum et al. 2026, *Oxford Open Immunology* iqag003, PMID 41859687) reports
+  ECE, Brier and NLL for uncertainty-quantification methods on **whole-protein** antigen
+  immunogenicity. Complementary scope: no peptide-HLA level, no patient clustering, no budget.
+- **arXiv 2604.13254** — calibrated conformal abstention for TCR-pMHC **binding**.
+- **CaliPPer** (arXiv 2606.07258) — distance-aware Platt recalibration, largely for binding and
+  presentation rather than immunogenicity.
+- **TESLA** (Wells et al. 2020, *Cell*) and the shortcut-bias analysis in Zhang et al. 2026
+  (*Cell Genomics*, PMC13347938) report ranking metrics only.
+
+## Limits worth stating plainly
+
+- Three tumour types, none of them colorectal. Nothing here licenses a claim about MSS CRC.
+- 70 patients is a small number of independent units for interval estimation, which is exactly
+  why the patient-clustered intervals are reported and the naive per-peptide ones are shown
+  beside them labelled incorrect.
+- The experimental labels are themselves a screen: "not immunogenic" means "not detected under
+  the assays applied to that patient", which is not the same as "not immunogenic". A calibration
+  result is a statement about agreement with these labels.
+- A calibration test can fail to reject "calibrated". It can never establish it.

examples/improve/budget_yield.py

@@ -0,0 +1,167 @@
+#!/usr/bin/env python
+"""Budgeted-yield sidecar for the IMPROVE fixture: what a slot budget actually buys.
+
+    uv run python examples/improve/budget_yield.py [--out DIR] [--n-boot N]
+
+Calibration is a property of probabilities; antigen selection is a property of a *budget*. A
+personalised vaccine encodes on the order of 20 neoepitopes and a TCR-T product carries 3-5, so
+the operative question is not "is peptide X immunogenic" but "of the k slots I am going to fill
+for this patient, how many carry a real target, and how many did the score promise?"
+
+For each model and each budget k this script takes every patient's own top-k candidates by
+score and reports three numbers per patient, averaged over patients with a patient-level
+bootstrap interval:
+
+  observed   -- immunogenic peptides actually in the top-k (the truth)
+  promised   -- the sum of the model's own published probabilities over that same top-k, which
+                is what those probabilities assert the yield will be
+  recalibrated -- the same sum after a cross-fitted intercept correction (predval's rung1),
+                fitted on other patients and applied here
+
+`promised` minus `observed` is the slot budget's miscalibration cost, in units of antigens. It
+is not a ranking error: the same top-k set is used for all three columns.
+
+This is a sidecar, not part of predval's contract: predval reports calibration, and the mapping
+from calibration to a decision is cohort-specific. It writes budget_yield.parquet and prints a
+summary. No network calls.
+"""
+
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+from predval import load_cohort, load_predictions
+
+# predval's own cross-fitting, reused rather than reimplemented so the correction applied here is
+# byte-for-byte the one the main report's rung1 column is scored on. Private because it is not a
+# contract surface; this example is allowed to reach for it, library code is not.
+from predval.recalibrate import _crossfit_predictions
+
+HERE = Path(__file__).resolve().parent
+BUDGETS = (3, 5, 10, 20)
+W = 26
+
+
+def _per_patient_topk(
+    frame: pd.DataFrame, k: int, score: str, extra: tuple[str, ...]
+) -> pd.DataFrame:
+    """Each patient's own top-k rows by `score`, ties broken deterministically by subject_id."""
+    ordered = frame.sort_values([score, "subject_id"], ascending=[False, True])
+    top = ordered.groupby("patient", sort=True).head(k)
+    cols = ["patient", "response", score, *extra]
+    return top[cols]
+
+
+def _cluster_bootstrap(
+    per_patient: pd.DataFrame, columns: list[str], n_boot: int, seed: int, ci: float
+) -> dict[str, tuple[float, float]]:
+    """Resample whole patients with replacement; return a percentile interval per column."""
+    rng = np.random.default_rng(seed)
+    values = per_patient[columns].to_numpy(dtype=float)
+    n = len(per_patient)
+    draws = np.empty((n_boot, len(columns)), dtype=float)
+    for b in range(n_boot):
+        idx = rng.integers(0, n, size=n)
+        draws[b] = values[idx].mean(axis=0)
+    lo = (1.0 - ci) / 2.0 * 100.0
+    return {
+        col: (float(np.percentile(draws[:, j], lo)), float(np.percentile(draws[:, j], 100.0 - lo)))
+        for j, col in enumerate(columns)
+    }
+
+
+def main(argv: list[str] | None = None) -> int:
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--out", type=Path, default=HERE / "out")
+    ap.add_argument("--n-boot", type=int, default=None)
+    args = ap.parse_args(argv)
+
+    cohort = load_cohort(HERE / "cohort.yaml")
+    predictions = load_predictions(HERE / "predictions.parquet")
+    spec = cohort.spec
+    n_boot = args.n_boot if args.n_boot is not None else spec.uncertainty.n_boot
+    seed = spec.uncertainty.seed
+    ci = spec.uncertainty.ci_level
+
+    table = cohort.table[[spec.subject_key, "patient", spec.outcome.field]].rename(
+        columns={spec.subject_key: "subject_id", spec.outcome.field: "response"}
+    )
+    prevalence = float(table["response"].mean())
+
+    rows = []
+    for model_id, grp in predictions.frame.groupby("model_id"):
+        merged = table.merge(
+            grp[["subject_id", "predicted"]], on="subject_id", how="inner", validate="1:1"
+        )
+        y = merged["response"].to_numpy(dtype=float)
+        p = merged["predicted"].to_numpy(dtype=float)
+        groups = merged["patient"].to_numpy()
+        k_folds = min(5, np.unique(groups).size)
+        cf, reason = _crossfit_predictions("rung1", y, p, groups, k_folds)
+        if cf is None:
+            raise SystemExit(f"{model_id}: rung1 cross-fit unavailable ({reason})")
+        merged["recalibrated"] = cf
+
+        for k in BUDGETS:
+            top = _per_patient_topk(merged, k, "predicted", ("recalibrated",))
+            per_patient = top.groupby("patient", sort=True).agg(
+                observed=("response", "sum"),
+                promised=("predicted", "sum"),
+                recalibrated=("recalibrated", "sum"),
+                slots=("response", "size"),
+            )
+            cols = ["observed", "promised", "recalibrated"]
+            intervals = _cluster_bootstrap(per_patient, cols, n_boot, seed, ci)
+            means = {c: float(per_patient[c].mean()) for c in cols}
+            rows.append(
+                {
+                    "model_id": model_id,
+                    "budget_k": k,
+                    "n_patients": int(len(per_patient)),
+                    "slots_median": float(per_patient["slots"].median()),
+                    **means,
+                    **{f"{c}_ci_low": intervals[c][0] for c in cols},
+                    **{f"{c}_ci_high": intervals[c][1] for c in cols},
+                    "patients_with_zero_hits": int((per_patient["observed"] == 0).sum()),
+                    "random_baseline": k * prevalence,
+                    "n_boot": n_boot,
+                    "seed": seed,
+                    "ci_level": ci,
+                }
+            )
+
+    out = pd.DataFrame(rows)
+    args.out.mkdir(parents=True, exist_ok=True)
+    path = args.out / "budget_yield.parquet"
+    out.to_parquet(path, index=False)
+
+    print(f"budget yield -> {path}")
+    print(
+        f"cohort prevalence {prevalence:.4%}; {out['n_patients'].iloc[0]} patients; "
+        f"B={n_boot} patient-level resamples, seed {seed}\n"
+    )
+    for k in BUDGETS:
+        sub = out[out["budget_k"] == k].sort_values("observed", ascending=False)
+        print(f"budget k={k} slots per patient  (random pick would yield {k * prevalence:.3f})")
+        print(
+            f"  {'model':{W}} {'observed':>24} {'promised':>10} {'recalib':>9} "
+            f"{'over by':>8} {'0-hit pts':>10}"
+        )
+        for _, r in sub.iterrows():
+            print(
+                f"  {r['model_id']:{W}} {r['observed']:8.3f} "
+                f"[{r['observed_ci_low']:.3f}, {r['observed_ci_high']:.3f}]"
+                f" {r['promised']:10.3f} {r['recalibrated']:9.3f}"
+                f" {r['promised'] / max(r['observed'], 1e-9):7.1f}x"
+                f" {r['patients_with_zero_hits']:6d}/{r['n_patients']}"
+            )
+        print()
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())