From 5d26be95dcbbf75a3502ed1ec9a39f10ea2e7dee Mon Sep 17 00:00:00 2001 From: aRealGem <203440826+aRealGem@users.noreply.github.com> Date: Sun, 4 Oct 2026 23:35:00 +0000 Subject: [PATCH 1/2] examples: add the IMPROVE third-party predictor panel, and fix the coverage-delta column it broke Before: predval's neoepitope evidence was one model. examples/improve/ showed that IMPROVE's own random forest ships probabilities anchored ~15x above the 2.67% base rate, with the whole fault in the intercept. Read strictly, that is a statement about one training pipeline. The report's per-model coverage-delta column, meanwhile, had never executed on any fixture: GUSTO, PCam and examples/improve all had uniform coverage across their rosters, and the column renders only when coverage differs. After: examples/improve_panel/ holds the patients, peptides, labels, clustering and both prespecified subgroups fixed and swaps in four immunogenicity predictors published by four other groups (DeepImmuno, DeepHLApan, DeepNetBim, iTTCA-RF), with IMPROVE's best variant kept as the in-house reference. One of them scores only 9-mers, so this is also the first fixture in the repo with non-uniform coverage -- which crashed the report renderer, and that is fixed here with regression tests covering the branch. The fault was in coverage_delta's pivot: dropna=False fills the full cross product of its index levels, inventing cells no stratum ever produced (stratum_kind="overall" crossed with a subgroup level), so a per-model lookup returned a Series rather than a scalar and np.isfinite raised. Those invented cells are now dropped at the source, since they carry neither subset, and the report pins all five keys when it reads the primary table's own cell. A missing delta now renders "n/a" rather than "0.0000", because the two are different statements. How: the fixture joins the release's benchmark_comparison table to its own patient-partitioned cross-validation output. Mut_peptide + HLA_allele is not a unique key -- 167 peptide-HLA pairs recur across patients and joining on them silently expands 17,520 rows to 17,884 with spurious label disagreements. Adding Expression, which is measured per patient, makes it exact, and build_fixture.py asserts uniqueness, the row count, label agreement and the [0, 1] bound rather than trusting any of them. Predictors whose published output is a binding or presentation score rather than an immunogenicity probability are excluded by name in cohort.yaml, as part of the prespecification; holding a presentation probability to the immunogenicity base rate would manufacture a miscalibration the model never asserted. skill_ladder.py is a sidecar, outside predval's contract like budget_yield.py: the report quotes Brier skill at the admissible rung, which is parsimony-first by design, so this scores every rung the ladder can fit and lets a reader check whether a skill result is a property of the scores or of the selection rule. It reports three outcomes rather than two, because "no rung reaches reliably positive skill" is a weaker claim than "every rung is reliably negative". No library behaviour changes for a uniform-coverage roster: the frozen GUSTO and PCam golden files still match byte for byte. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0138yQfFhvnVRHd1XhTEz54w --- .gitignore | 5 + examples/improve_panel/README.md | 165 +++++++++++++++++++ examples/improve_panel/build_fixture.py | 174 ++++++++++++++++++++ examples/improve_panel/cohort.yaml | 117 ++++++++++++++ examples/improve_panel/run_evaluation.py | 198 +++++++++++++++++++++++ examples/improve_panel/skill_ladder.py | 179 ++++++++++++++++++++ src/predval/evaluate.py | 6 + src/predval/report.py | 19 ++- tests/test_report.py | 43 ++++- 9 files changed, 900 insertions(+), 6 deletions(-) create mode 100644 examples/improve_panel/README.md create mode 100644 examples/improve_panel/build_fixture.py create mode 100644 examples/improve_panel/cohort.yaml create mode 100644 examples/improve_panel/run_evaluation.py create mode 100644 examples/improve_panel/skill_ladder.py diff --git a/.gitignore b/.gitignore index fb80953..389b642 100644 --- a/.gitignore +++ b/.gitignore @@ -34,3 +34,8 @@ validation/predrupdate/gusto_estimates.csv examples/improve/cohort.parquet examples/improve/predictions.parquet examples/improve/out/ + +# And for the third-party predictor panel on the same cohort. +examples/improve_panel/cohort.parquet +examples/improve_panel/predictions.parquet +examples/improve_panel/out/ diff --git a/examples/improve_panel/README.md b/examples/improve_panel/README.md new file mode 100644 index 0000000..be5b936 --- /dev/null +++ b/examples/improve_panel/README.md @@ -0,0 +1,165 @@ +# IMPROVE neoepitope immunogenicity, third-party predictor panel — predval example cohort + +The companion to [`examples/improve/`](../improve/). Same 17,520 candidate neoepitopes, same 70 +patients, same experimentally determined T-cell recognition labels — scored by immunogenicity +predictors published by **other groups** instead of by IMPROVE's own random forest. + +## The question this fixture asks + +`examples/improve/` found that IMPROVE's published probabilities sit roughly 15x above the base +rate, with the whole fault in the intercept, repairable by a single global shift. On its own that +is a statement about one training pipeline: train on a class-balanced set, deploy at natural +prevalence, and this is what you get. + +The obvious objection is that it is one model. So this fixture holds everything else fixed — the +patients, the peptides, the assay, the clustering, the two prespecified subgroups — and swaps in +four independently developed and independently published predictors from four other groups. If +the fault is a pipeline artefact, it should look different in each. If it is what the field does, +it should look the same in all of them. + +**This is explicitly not independent replication.** The patients, labels and assay are shared +with `examples/improve/`, so nothing here tests whether the finding survives a new patient +sample. It tests breadth across *predictors*, not breadth across *cohorts*. Replication needs a +second cohort that publishes per-peptide scores with a patient-level clustering unit at natural +prevalence, which is an acquisition problem, not a code problem. + +## Provenance + +Source: the predictor comparison table released with Borch et al. 2024, *IMPROVE: a feature model +to predict neoepitope immunogenicity through broad-scale validation of T cell recognition*, +Frontiers in Immunology — [github.com/SRHgroup/IMPROVE_paper](https://github.com/SRHgroup/IMPROVE_paper), +`data.zip` → `data/benchmark_comparison/neoepitopes_predictions.tsv`, joined to that same +release's patient-partitioned cross-validation output. + +The join is the part worth checking, and `build_fixture.py` documents it in full. In short: +`Mut_peptide` + `HLA_allele` is **not** a unique key — 167 peptide-HLA pairs recur across +patients, and joining on them produces a silent many-to-many expansion (17,884 rows from 17,520) +with apparent label disagreements that are really the cross-product of those duplicates. Adding +`Expression`, which is measured per patient, makes it exact: 17,520 of 17,520 rows, both sides +unique, `validate="1:1"`, zero label disagreements. The builder asserts all of that rather than +trusting it, and asserts the `[0, 1]` bound after rescaling. + +## Roster + +Admitted, because every member's published output is nominally an immunogenicity **probability**: + +- `deepimmuno` — DeepImmuno. +- `deephlapan_immuno` — DeepHLApan's *immunogenic* score, not its binding score. +- `deepnetbim_immuno` — DeepNetBim's immunogenicity probability. 9-mers only, so it covers 10,561 + of 17,520 peptides; see the coverage note below. +- `ittca_rf` — iTTCA-RF, published on a 0–100 scale and divided by 100 in the builder. +- `improve_tme_incl` — IMPROVE's best-discriminating variant, carried over from + `examples/improve/` as the in-house reference point. It is **not** a competitor on this table: + it was cross-validated on this very cohort, so it reads with a home advantage the other four + do not have. + +Excluded, and the exclusion list is part of the prespecification because scoring these as +immunogenicity probabilities would be a strawman rather than a finding: + +- **MHCflurry** (affinity, processing, presentation) and **HLAthena** (MSiC / MSiE / MSiCE) + predict antigen *presentation*. A presentation probability of 0.5 is not a claim that the + peptide is immunogenic, so holding it to the immunogenicity base rate would manufacture a + miscalibration the model never asserted. +- **MixMHCpred**, **RankEL**, **RankBA** — binding predictors, and the ranks are percentiles. +- **PRIME** — binding-weighted, not published on a probability scale. +- **IEDB immunogenicity** — runs from −0.914 to 0.639 on this cohort. It is a score, not a + probability, and predval would reject it at the `[0, 1]` bound anyway. + +The distinction being enforced: a calibration critique applies only to a model that claims to +emit a probability of the event the labels record. Anything else is a ranking tool, and +`examples/improve/`'s own result is that the ranking is the part that works. + +## Coverage, and why the floor was left where it is + +`coverage.min_fraction` stays at 0.95 knowing that DeepNetBim breaches it. Scoring only 9-mers is +a real property of that model, and the right response is a coverage flag plus an honest +common-subset comparison — not a relaxed floor, and not imputation. A model that did not score a +peptide gets no row for it. + +The excluded 40% is far above the informative-selection threshold, so the report also cautions +that selection into the common subset matters here: the intersection is 9-mers only, which is the +easiest stratum for every binding-derived feature. **Read each model's full-subset column for its +own behaviour, and the common-subset column only for model-to-model comparison.** + +Building this fixture is also what exposed a real defect in predval itself: the per-model +coverage-delta column in the report renders only when coverage differs across the roster, and +every prior fixture (GUSTO, PCam, `examples/improve/`) had uniform coverage, so that branch had +never executed. DeepNetBim is the first partial-coverage model in the repo's history and it +crashed the renderer. Fixed, with regression tests covering non-uniform coverage, in the same +change as this directory. + +## Unit of analysis and prespecified subgroups + +Identical to `examples/improve/`, carried over unchanged so the two fixtures are readable side by +side, and neither subgroup was chosen after seeing a result on this roster. + +The cluster is the **patient**. Peptides from one patient share that patient's HLA genotype, +mutational process, tumour transcriptome and microenvironment — true of these predictors' inputs +exactly as it was of IMPROVE's. Subgroups are `tumour_cohort` (the release's three accrual +groups) and `peptide_length` (9-mer versus everything else). + +## The skill-ladder sidecar + +`skill_ladder.py` exists to answer the first objection a reviewer should raise. + +predval's report quotes one Brier skill score per model, at the **admissible** rung — the lowest +rung whose cross-fitted paired gain interval excludes zero. That rule is parsimony-first by +design: do not buy a more complex correction than you can prove you need. But it means the +headline number depends on a selection rule, and on this roster most models come out at or below +zero skill. Maybe a higher rung would have rescued them and the rule hid it. + +So the sidecar drops the rule and scores **every** rung the ladder can fit — rung0 as published, +then cross-fitted rung1, rung2 and rung3 — each with predval's own `brier_skill_interval` and the +same patient-clustered bootstrap, B and seed. A rung the ladder cannot fit for a model is +reported unavailable with its reason, never imputed. If no rung reaches reliably positive skill, +that is a property of the scores rather than of the selection rule, and that is the claim this +script lets a reader check. It reports three outcomes, not two: "reliably positive at some rung", +"reliably negative at every fitted rung", and "no rung reaches reliably positive skill" — the +last is not the same claim as the middle one, and collapsing them would overstate the result in +exactly the direction this fixture argues. + +Like `budget_yield.py` in `examples/improve/`, it is a sidecar and deliberately not part of +predval's contract. It reuses predval's private `_crossfit_predictions` rather than +reimplementing it, so the corrections scored here are byte-for-byte the ones the report's own +rung columns are scored on. + +## Reproducing + +```sh +# 1. fetch the release and extract both archives +git clone https://github.com/SRHgroup/IMPROVE_paper.git +unzip -q IMPROVE_paper/data.zip -x '__MACOSX/*' -d improve_data +unzip -q IMPROVE_paper/results.zip -x '__MACOSX/*' -d improve_results + +# 2. build the fixture (writes cohort.parquet and predictions.parquet here). +# Two explicit paths, because the two tables come from two different archives. +uv run python examples/improve_panel/build_fixture.py \ + --panel-tsv improve_data/data/benchmark_comparison/neoepitopes_predictions.tsv \ + --cv-txt 'improve_results/results copy/5_fold_CV/TME_included/pred_df_TME_included.txt' + +# 3. the main evaluation -- B=2000 patient-clustered bootstrap, several minutes +uv run python examples/improve_panel/run_evaluation.py + +# 4. the skill-at-every-rung sidecar +uv run python examples/improve_panel/skill_ladder.py +``` + +`run_evaluation.py` writes `out/report.html`, `out/findings.json` and the parquet artefacts; +`skill_ladder.py` writes `out/skill_ladder.parquet`. + +## Limits worth stating plainly + +- **Not a replication.** Shared patients, labels and assay with `examples/improve/`. Breadth + across predictors only. +- **`improve_tme_incl` is not a competitor.** It was cross-validated on this cohort; the other + four were trained elsewhere and are being applied to it cold. That asymmetry is the normal + condition for anyone selecting antigens with an off-the-shelf tool, which is why it is on the + table at all — but it is not a like-for-like ranking of model quality, and reading it as one + would be wrong. +- **Discrimination near chance for the third-party models is itself a finding about this cohort, + not only about the tools.** A cohort on which nothing discriminates cannot establish that a + tool is uninformative in general. +- Three tumour types, none colorectal. Nothing here licenses a claim about MSS CRC. +- The labels are a screen: "not immunogenic" means "not detected under the assays applied to that + patient". A calibration result is a statement about agreement with these labels. +- A calibration test can fail to reject "calibrated". It can never establish it. diff --git a/examples/improve_panel/build_fixture.py b/examples/improve_panel/build_fixture.py new file mode 100644 index 0000000..2cd8317 --- /dev/null +++ b/examples/improve_panel/build_fixture.py @@ -0,0 +1,174 @@ +#!/usr/bin/env python +"""Build the IMPROVE third-party predictor panel fixture. + +Same 70-patient cohort as examples/improve/, scored by independently published predictors from +other groups instead of by IMPROVE's own random forest. The question this fixture asks is +whether the miscalibration found in examples/improve/ is a property of that one model's +training pipeline or a property of how the field ships probabilities. + + uv run python examples/improve_panel/build_fixture.py --src + +Writes cohort.parquet and predictions.parquet next to cohort.yaml. No network calls. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import pandas as pd + +HERE = Path(__file__).resolve().parent + +#: model_id -> (column in the released panel table, divisor to put it on the [0, 1] probability +#: scale). Only predictors whose published output is nominally an immunogenicity PROBABILITY are +#: admitted; see cohort.yaml for the exclusions and the reasoning behind them. +PANEL = { + "deepimmuno": ("DeepImmuno", 1.0), + "deephlapan_immuno": ("DeepHLApan immunogenic score", 1.0), + "deepnetbim_immuno": ("DeepNetBim immunogenicity probability", 1.0), + "ittca_rf": ("iTTCA-RF probability", 100.0), +} + +#: The in-house reference point, carried over from examples/improve/ so the panel is read +#: against a model whose calibration this repo has already characterised. +IMPROVE_REFERENCE = "improve_tme_incl" + + +def _join_key(frame: pd.DataFrame) -> pd.Series: + """The only key that joins the panel table to the cross-validation tables 1:1. + + `Mut_peptide` + `HLA_allele` is NOT unique: 167 peptide-HLA pairs recur across patients, and + joining on them silently produces a many-to-many expansion (17,884 rows from 17,520) with + apparent label disagreements that are really the cross-product of those duplicates. Adding + `Expression` disambiguates, because expression is measured per patient. The panel table + rounds it to three decimals while the CV tables carry full precision, so both sides are + formatted to three decimals before comparison; at that precision the join is exactly 1:1 + over all 17,520 rows with zero label disagreements, which `main` asserts rather than trusts. + """ + return ( + frame["Mut_peptide"].astype(str) + + "|" + + frame["HLA_allele"].astype(str) + + "|" + + frame["Expression"].map(lambda v: f"{float(v):.3f}") + ) + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + # Two explicit file paths rather than one archive root: the two tables come from two + # different zips in the release, and where each lands depends on how they were extracted. + ap.add_argument( + "--panel-tsv", + type=Path, + required=True, + help="data.zip -> data/benchmark_comparison/neoepitopes_predictions.tsv", + ) + ap.add_argument( + "--cv-txt", + type=Path, + required=True, + help="results.zip -> 5_fold_CV/TME_included/pred_df_TME_included.txt", + ) + ap.add_argument("--out", type=Path, default=HERE) + args = ap.parse_args(argv) + + panel = pd.read_csv(args.panel_tsv, sep="\t") + cv = pd.read_csv(args.cv_txt, sep=r"\s+") + + panel["join_key"] = _join_key(panel) + cv["join_key"] = _join_key(cv) + if not panel["join_key"].is_unique or not cv["join_key"].is_unique: + raise SystemExit("join key is not unique on both sides; refusing to guess the pairing") + + merged = panel.merge( + cv[["join_key", "Patient", "cohort", "Partition", "response", "prediction_rf"]], + on="join_key", + how="inner", + validate="1:1", + suffixes=("_panel", ""), + ) + if len(merged) != len(cv): + raise SystemExit(f"join dropped rows: {len(merged)} of {len(cv)}") + if (merged["response_panel"] != merged["response"]).any(): + raise SystemExit("labels disagree between the panel table and the CV table") + + merged["subject_id"] = ( + merged["Patient"].astype(str) + + "|" + + merged["Mut_peptide"].astype(str) + + "|" + + merged["HLA_allele"].astype(str) + + "|" + + merged["Expression"].map(lambda v: f"{float(v):.3f}") + ) + if not merged["subject_id"].is_unique: + raise SystemExit("subject_id is not unique") + + pep_len = merged["Mut_peptide"].str.len() + cohort = pd.DataFrame( + { + "subject_id": merged["subject_id"], + "patient": merged["Patient"].astype(str), + "tumour_cohort": merged["cohort"].astype(str), + "peptide_length": pep_len.map(lambda n: "9mer" if n == 9 else "non_9mer"), + "cv_partition": merged["Partition"].astype(int), + "mut_peptide": merged["Mut_peptide"].astype(str), + "hla_allele": merged["HLA_allele"].astype(str), + "response": merged["response"].astype(int), + } + ).reset_index(drop=True) + + rows = [ + pd.DataFrame( + { + "subject_id": merged["subject_id"], + "model_id": IMPROVE_REFERENCE, + "predicted": merged["prediction_rf"].astype(float), + } + ) + ] + for model_id, (column, divisor) in PANEL.items(): + score = pd.to_numeric(merged[column], errors="coerce") / divisor + keep = score.notna() + out_of_range = keep & ((score < 0.0) | (score > 1.0)) + if out_of_range.any(): + raise SystemExit( + f"{model_id}: {int(out_of_range.sum())} scores outside [0, 1] after scaling " + f"by {divisor}; the published scale is not what this mapping assumes" + ) + # A model that did not score a peptide gets NO ROW. predval measures coverage against + # the declared floor and refuses to read "no opinion" as "confidently negative". + rows.append( + pd.DataFrame( + { + "subject_id": merged.loc[keep, "subject_id"], + "model_id": model_id, + "predicted": score[keep].astype(float), + } + ) + ) + + predictions = pd.concat(rows, ignore_index=True) + + args.out.mkdir(parents=True, exist_ok=True) + cohort.to_parquet(args.out / "cohort.parquet", index=False) + predictions.to_parquet(args.out / "predictions.parquet", index=False) + + prevalence = float(cohort["response"].mean()) + print(f"cohort: {len(cohort):,} peptides, {cohort['patient'].nunique()} patients") + print(f"outcome: {int(cohort['response'].sum())} immunogenic ({prevalence:.3%} prevalence)") + print(f"predictions: {len(predictions):,} rows, {predictions['model_id'].nunique()} models") + for model_id, grp in predictions.groupby("model_id"): + mean_p = float(grp["predicted"].mean()) + print( + f" {model_id:22} n={len(grp):6,} coverage={len(grp) / len(cohort):.4f} " + f"mean={mean_p:.4f} ratio={mean_p / prevalence:5.1f}x " + f"range=[{grp['predicted'].min():.4f}, {grp['predicted'].max():.4f}]" + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/improve_panel/cohort.yaml b/examples/improve_panel/cohort.yaml new file mode 100644 index 0000000..b700a5c --- /dev/null +++ b/examples/improve_panel/cohort.yaml @@ -0,0 +1,117 @@ +# IMPROVE cohort, third-party predictor panel. +# +# The same 17,520 candidate neoepitopes from 70 patients as examples/improve/, with the same +# experimentally determined T-cell recognition outcomes (467 immunogenic, 2.666% prevalence), +# scored by independently published immunogenicity predictors from other groups. +# +# examples/improve/ established that IMPROVE's own random forest ships probabilities anchored +# roughly 15x above the base rate, with the whole fault in the intercept. That result on its own +# is a statement about one training pipeline. This fixture asks the question it provokes: is +# that one model's artefact, or is it what the field does? Same cohort, same labels, same +# clustering, different model families from different groups. +# +# This is explicitly NOT independent replication. The patients, the labels and the assay are +# shared with examples/improve/, so nothing here tests whether the finding survives a new +# patient sample. It tests breadth across predictors, not breadth across cohorts. +# +# Source: the predictor table released with Borch et al. 2024 (Front Immunol) at +# github.com/SRHgroup/IMPROVE_paper, data.zip, +# data/benchmark_comparison/neoepitopes_predictions.tsv, joined to that release's own +# patient-partitioned cross-validation table. The join key and why it is the only correct one +# are documented in build_fixture.py. +# +# Regenerate the parquet files with: +# uv run python examples/improve_panel/build_fixture.py --src + +cohort_id: improve-neoepitope-thirdparty-panel +version: 0 +subject_key: subject_id +data: cohort.parquet + +outcome: + type: binary + field: response + positive_label: 1 + +# Identical to examples/improve/: the cluster is the patient, because peptides from one patient +# share that patient's HLA genotype, mutational process, tumour transcriptome and +# microenvironment. That is true of these predictors' inputs exactly as it was of IMPROVE's. +unit_noun: peptide + +clustering: + field: patient + name: patient + +# 0.95 is deliberately left where examples/improve/ set it, knowing one roster member will +# breach it. DeepNetBim scores only 9-mers and therefore covers 10,561 of 17,520 peptides, and +# the right response to that is a coverage flag and an honest common-subset comparison, not a +# relaxed floor and not imputation. The excluded 40% is far above common_warn_frac, so the +# report will also caution that selection into the common subset is informative here: the +# intersection is 9-mers only, which is the easiest stratum for every binding-derived feature. +# Read the full-subset column for each model's own behaviour and the common-subset column only +# for model-to-model comparison. +coverage: + min_fraction: 0.95 + compare_on: both + +completeness: + require_outcome: true + on_violation: drop_and_report + +# The same two prespecified subgroups as examples/improve/, carried over unchanged so the two +# fixtures are readable side by side. Neither was chosen after seeing a result on this roster. +subgroups: + - name: tumour_cohort + field: tumour_cohort + - name: peptide_length + field: peptide_length + +thresholds: [0.5] + +uncertainty: + n_boot: 2000 + seed: 1337 + ci_level: 0.95 + show_naive_ci: true + +# Roster: 5 declared, 5 expected present. +# +# ADMITTED -- every member's published output is nominally an immunogenicity probability: +# deepimmuno DeepImmuno, a deep model trained on immunogenicity assay data. +# deephlapan_immuno DeepHLApan's immunogenic score (its binding score is a separate output +# and is NOT in this roster). +# deepnetbim_immuno DeepNetBim's immunogenicity probability (likewise distinct from its +# binding output). 9-mers only; see the coverage note above. +# ittca_rf iTTCA-RF's probability, published on a 0-100 scale and divided by 100 in +# build_fixture.py. Its observed floor on this cohort is exactly 0.500 -- +# it never expresses a probability below a coin flip. +# improve_tme_incl IMPROVE's best-discriminating variant, carried from examples/improve/ as +# the in-house reference point this repo has already characterised. +# +# EXCLUDED, and this list is part of the prespecification because scoring these as +# immunogenicity probabilities would be a strawman rather than a finding: +# MHCflurry (MHC affinity, processing, presentation) -- predicts antigen PRESENTATION. A +# presentation probability of 0.5 is not a claim that the peptide is immunogenic, so holding +# it to the immunogenicity base rate would manufacture a miscalibration that the model never +# asserted. The same reasoning excludes HLAthena's MSiC / MSiE / MSiCE. +# MixMHCpred, RankEL, RankBA -- binding predictors, and the ranks are percentiles, not +# probabilities. +# PRIME (score and rank) -- binding-weighted and not published on a probability scale. +# IEDB immunogenicity -- ranges from -0.914 to 0.639 on this cohort; it is a score, not a +# probability, and predval would reject it at the [0, 1] bound anyway. +# +# The distinction being enforced: a calibration critique only applies to a model that claims to +# emit a probability of the event the labels record. Anything else is a ranking tool, and +# examples/improve/'s own result is that the ranking is the part that works. +expected_models: + - improve_tme_incl + - deepimmuno + - deephlapan_immuno + - deepnetbim_immuno + - ittca_rf +on_missing_model: fail + +# No member is a blend of several separately trained models. Same reasoning as +# examples/improve/: a single trained model's internal bagging is part of its fitting procedure, +# not a post-hoc blend over a roster. +ensemble_members: [] diff --git a/examples/improve_panel/run_evaluation.py b/examples/improve_panel/run_evaluation.py new file mode 100644 index 0000000..73d2084 --- /dev/null +++ b/examples/improve_panel/run_evaluation.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python +"""Evaluate the IMPROVE third-party predictor panel end to end and write the artefacts. + + uv run python examples/improve_panel/run_evaluation.py [--out DIR] [--n-boot N] + +Writes metrics.parquet, fragility.parquet, coverage.parquet, manifest.json, report.html and +findings.json, then prints the summary: prevalence anchoring, discrimination, the +unit-of-analysis exhibit over patients, calibration slope and intercept, the recalibration +ladder (as-published versus cross-fitted) and fragility to dropping any single patient. +At the declared B=2000 over 7 models this takes a few minutes; pass --n-boot for a smoke run. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import pandas as pd + +from predval import evaluate, load_cohort, load_predictions, write_findings, write_report +from predval.metrics import optimism + +HERE = Path(__file__).resolve().parent +W = 22 # model ids here are shorter than examples/improve/'s + + +def _overall(metrics: pd.DataFrame, metric: str, ci_method: str | None = None) -> pd.DataFrame: + sel = ( + (metrics["stratum_kind"] == "overall") + & (metrics["subset"] == "common") + & (metrics["metric"] == metric) + ) + if ci_method is not None: + sel &= metrics["ci_method"] == ci_method + return metrics[sel] + + +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, help="override uncertainty.n_boot") + args = ap.parse_args(argv) + + cohort = load_cohort(HERE / "cohort.yaml") + predictions = load_predictions(HERE / "predictions.parquet") + + if args.n_boot is not None: + unc = cohort.spec.uncertainty.model_copy(update={"n_boot": args.n_boot}) + spec = cohort.spec.model_copy(update={"uncertainty": unc}) + cohort = type(cohort)(**{**cohort.__dict__, "spec": spec}) + + result = evaluate(cohort, predictions) + paths = result.write(args.out) + report_path = write_report(result, args.out / "report.html") + findings_path = write_findings(result, args.out / "findings.json") + + print(f"cohort: {cohort.n_subjects:,} peptides") + print(f"models: {len(predictions.model_ids)}") + print(f"metrics: {len(result.metrics):,} rows -> {paths['metrics']}") + print(f"report: {report_path}") + print(f"findings: {findings_path}") + print() + + print("FLAGS") + for flag in result.flags: + print(f" [{flag.severity:7}] {flag.code}: {flag.message}") + print() + + # ---- prevalence anchoring: the first thing to look at under 2.7% positives ------------- + frame = predictions.frame + prev = float(cohort.table[cohort.spec.outcome.field].mean()) + print(f"Prevalence anchoring -- observed immunogenic rate {prev:.4%}") + print(f" {'model':{W}} {'mean pred':>10} {'ratio':>7} {'min':>8} {'max':>8}") + for model, grp in frame.groupby("model_id"): + mean_p = float(grp["predicted"].mean()) + print( + f" {model:{W}} {mean_p:10.4f} {mean_p / prev:6.1f}x " + f"{grp['predicted'].min():8.4f} {grp['predicted'].max():8.4f}" + ) + print() + + auroc = _overall(result.metrics, "auroc", "cluster_bootstrap").sort_values( + "value", ascending=False + ) + print("AUROC by model (common subset, patient-clustered bootstrap 95% CI)") + for _, r in auroc.iterrows(): + print(f" {r['model_id']:{W}} {r['value']:.4f} [{r['ci_low']:.4f}, {r['ci_high']:.4f}]") + print() + + ap_rows = _overall(result.metrics, "average_precision", "cluster_bootstrap").sort_values( + "value", ascending=False + ) + if len(ap_rows): + print(f"Average precision (no-skill baseline = prevalence = {prev:.4f})") + for _, r in ap_rows.iterrows(): + print( + f" {r['model_id']:{W}} {r['value']:.4f} [{r['ci_low']:.4f}, {r['ci_high']:.4f}]" + ) + print() + + naive = _overall(result.metrics, "auroc") + print("Unit-of-analysis exhibit -- AUROC interval width, patient vs naive per-peptide") + for model, grp in naive.groupby("model_id"): + widths = {r["ci_method"]: r["ci_high"] - r["ci_low"] for _, r in grp.iterrows()} + if "naive_row_bootstrap" in widths and widths["naive_row_bootstrap"] > 0: + ratio = widths["cluster_bootstrap"] / widths["naive_row_bootstrap"] + print( + f" {model:{W}} patient {widths['cluster_bootstrap']:.4f} " + f"naive {widths['naive_row_bootstrap']:.4f} ratio {ratio:5.1f}x" + ) + print() + + for metric, blurb in ( + ("calibration_slope", "1.0 = correct spread; <1 = over-spread/over-confident"), + ("calibration_intercept", "0.0 = correctly anchored; <0 = systematically too high"), + ): + rows = _overall(result.metrics, metric, "cluster_robust_t").sort_values("value") + if not len(rows): + continue + print(f"{metric.replace('_', ' ').capitalize()} ({blurb})") + for _, r in rows.iterrows(): + interval = f"[{r['ci_low']:8.3f}, {r['ci_high']:8.3f}]" + print(f" {r['model_id']:{W}} {r['value']:8.3f} {interval}") + print() + + lad = _overall(result.metrics, "brier").pivot_table( + index="model_id", columns=["rung", "fit_mode"], values="value" + ) + + def cell(model: str, rung: str, mode: str) -> float: + col = (rung, mode) + return float(lad.loc[model, col]) if col in lad.columns else float("nan") + + def cf_gain(model: str) -> float: + r0 = cell(model, "rung0", "apparent") + cf = min(cell(model, r, "crossfit") for r in ("rung1", "rung2", "rung3")) + return r0 - cf + + order = sorted(lad.index, key=cf_gain, reverse=True) + print("Recalibration ladder -- Brier, as-published vs cross-fitted by held-out patients") + print(f" reference: a constant predictor at the prevalence scores {prev * (1 - prev):.4f}") + print(" optim_r2 > 0 == the apparent rung2 fit flattered itself versus held-out patients") + print( + f" {'model':{W}} {'rung0':>8} {'r1_cf':>8} {'r2_cf':>8} {'r3_cf':>8} " + f"{'gain_cf':>8} {'optim_r2':>9}" + ) + for model in order: + r0 = cell(model, "rung0", "apparent") + r1_cf = cell(model, "rung1", "crossfit") + r2_cf = cell(model, "rung2", "crossfit") + r3_cf = cell(model, "rung3", "crossfit") + gain = r0 - min(r1_cf, r2_cf, r3_cf) + optim = optimism("brier", cell(model, "rung2", "apparent"), r2_cf) + print( + f" {model:{W}} {r0:8.4f} {r1_cf:8.4f} {r2_cf:8.4f} {r3_cf:8.4f} " + f"{gain:8.4f} {optim:+9.4f}" + ) + print() + + gains = _overall(result.metrics, "paired_gain_brier") + if len(gains): + print("Paired Brier gain, cross-fitted (CI excluding 0 == a demonstrated repair)") + for model, grp in gains.groupby("model_id"): + parts = [] + for _, r in grp.sort_values("rung").iterrows(): + mark = " " + if pd.notna(r["ci_low"]) and pd.notna(r["ci_high"]): + if r["ci_low"] > 0: + mark = "+" + elif r["ci_high"] < 0: + mark = "-" + parts.append(f"{r['rung']}{mark}{r['value']:+.5f}") + print(f" {model:{W}} " + " ".join(parts)) + print(" + = repair demonstrated - = recalibration reliably HARMS blank = unproven") + print() + + frag = result.fragility[ + (result.fragility["stratum_kind"] == "overall") + & (result.fragility["subset"] == "common") + & (result.fragility["metric"] == "auroc") + ].sort_values("max_abs_delta", ascending=False) + print("Fragility -- largest AUROC change from dropping any single patient") + for _, r in frag.head(7).iterrows(): + print(f" {r['model_id']:{W}} {r['max_abs_delta']:.4f} (patient {r['culprit_cluster']})") + + boundary = _overall(result.metrics, "boundary_count") + boundary = boundary[boundary["value"] > 0] + if len(boundary): + print() + print("Data quality -- predictions exactly at 0 or 1 (eps-clip load)") + for _, r in boundary.iterrows(): + print(f" {r['model_id']:{W}} {int(r['value']):5d}") + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/improve_panel/skill_ladder.py b/examples/improve_panel/skill_ladder.py new file mode 100644 index 0000000..dbcb9c5 --- /dev/null +++ b/examples/improve_panel/skill_ladder.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python +"""Brier skill at EVERY rung, for the panel: is negative skill a model fact or a rung choice? + + uv run python examples/improve_panel/skill_ladder.py [--out DIR] [--n-boot N] + +predval's report quotes one Brier skill score per model, at the *admissible* rung -- the lowest +rung whose cross-fitted paired gain interval excludes zero. That rule is parsimony-first by +design: do not buy a more complex correction than you can prove you need. It also means the +headline number depends on a selection rule, and on this roster four of five models come out with +NEGATIVE skill, i.e. worse than quoting the cohort base rate to every patient. A reviewer's first +objection is the right one: maybe a higher rung would have rescued them and the rule hid it. + +So this sidecar drops the selection rule and scores every rung the ladder could fit: + + rung0 the published probabilities, uncorrected + rung1 cross-fitted intercept shift (level) + rung2 cross-fitted intercept + slope (level and spread) + rung3 cross-fitted restricted cubic spline + +Each is scored with predval's own `brier_skill_interval`, anchored at 0 = always predict the +observed prevalence and 1 = perfect, with the same patient-clustered bootstrap, B and seed the +cohort spec declares. A rung the ladder cannot fit for a model is reported unavailable with its +reason, never imputed. + +If a model's interval stays below zero at every rung, "negative skill" is a property of the +scores, not of the admissible-rung rule. That is the claim this script exists to let a reader +check, and it is the claim predval's own report is not structured to make. + +Sidecar, not part of predval's contract: the library reports the admissible rung on purpose. +Writes skill_ladder.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 corrections scored here are +# byte-for-byte those the main report's rung1-3 columns are 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 +from predval.uncertainty import brier_skill_interval + +HERE = Path(__file__).resolve().parent +RUNGS = ("rung1", "rung2", "rung3") +MAX_FOLDS = 5 +W = 22 + + +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, help="override uncertainty.n_boot") + 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"} + ) + + rows = [] + for model_id, grp in predictions.frame.groupby("model_id"): + # inner join: a peptide this model did not score contributes no row, and the model's own + # prevalence reference is recomputed on the rows it does cover (see the coverage note in + # cohort.yaml -- the full and common subsets are different questions). + 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(MAX_FOLDS, np.unique(groups).size) + + corrected: dict[str, np.ndarray | None] = {"rung0": p} + reasons: dict[str, str | None] = {"rung0": None} + for rung in RUNGS: + cf, reason = _crossfit_predictions(rung, y, p, groups, k_folds) + corrected[rung] = cf + reasons[rung] = reason + + for rung, q in corrected.items(): + if q is None: + rows.append( + { + "model_id": model_id, + "rung": rung, + "n": int(len(merged)), + "available": False, + "unavailable_reason": reasons[rung], + "bss": float("nan"), + "bss_ci_low": float("nan"), + "bss_ci_high": float("nan"), + "brier": float("nan"), + "prevalence": float(np.mean(y)), + "n_boot": n_boot, + "seed": seed, + "ci_level": ci, + } + ) + continue + bss, interval = brier_skill_interval(y, q, groups, ci, n_boot=n_boot, seed=seed) + rows.append( + { + "model_id": model_id, + "rung": rung, + "n": int(len(merged)), + "available": True, + "unavailable_reason": None, + "bss": float(bss), + "bss_ci_low": float(interval.low), + "bss_ci_high": float(interval.high), + "brier": float(np.mean((q - y) ** 2)), + "prevalence": float(np.mean(y)), + "n_boot": n_boot, + "seed": seed, + "ci_level": ci, + } + ) + + out = pd.DataFrame(rows) + args.out.mkdir(parents=True, exist_ok=True) + path = args.out / "skill_ladder.parquet" + out.to_parquet(path, index=False) + + print(f"skill ladder -> {path}") + print(f"B={n_boot} patient-clustered resamples, seed {seed}, {ci:.0%} intervals") + print("BSS: 0 = always predict this model's own observed prevalence, 1 = perfect.") + print("A rung whose whole interval is below 0 is reliably WORSE than quoting the base rate.\n") + + for model_id, grp in out.groupby("model_id"): + prev = float(grp["prevalence"].iloc[0]) + n = int(grp["n"].iloc[0]) + print( + f"{model_id} (n={n:,}, prevalence {prev:.4%}, reference Brier {prev * (1 - prev):.4f})" + ) + for rung in ("rung0", *RUNGS): + r = grp[grp["rung"] == rung].iloc[0] + if not bool(r["available"]): + print(f" {rung:{8}} unavailable: {r['unavailable_reason']}") + continue + mark = " " + if np.isfinite(r["bss_ci_high"]) and r["bss_ci_high"] < 0: + mark = "-" # reliably worse than the base rate + elif np.isfinite(r["bss_ci_low"]) and r["bss_ci_low"] > 0: + mark = "+" # reliably better + interval = f"[{100 * r['bss_ci_low']:+8.2f}%, {100 * r['bss_ci_high']:+8.2f}%]" + print( + f" {rung:{8}} brier {r['brier']:.4f} BSS {100 * r['bss']:+8.2f}% " + f"{interval} {mark}" + ) + ok = grp[grp["available"]] + if len(ok): + # Three outcomes, not two. "No rung is reliably positive" is the finding; it is NOT + # the same claim as "every rung is reliably negative", and collapsing the two would + # overstate the result in exactly the direction this fixture is arguing. + if (ok["bss_ci_low"] > 0).any(): + verdict = "reliably POSITIVE skill at some rung" + elif (ok["bss_ci_high"] < 0).all(): + verdict = "reliably NEGATIVE skill at every fitted rung" + else: + verdict = "no rung reaches reliably positive skill (best rung straddles zero)" + print(f" -> {verdict}") + print() + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/predval/evaluate.py b/src/predval/evaluate.py index 2268207..42bc9bf 100644 --- a/src/predval/evaluate.py +++ b/src/predval/evaluate.py @@ -100,6 +100,12 @@ def coverage_delta(metrics: pd.DataFrame) -> pd.DataFrame: ).reset_index() if "full" not in wide or "common" not in wide: return pd.DataFrame(columns=[*keys, "cov_delta"]) + # `dropna=False` keeps a cell whose `common` value is missing -- that absence is the thing + # this column exists to show. It also fills in the full cross product of the index levels, + # which invents cells that never existed (stratum_kind="overall" crossed with a + # subgroup_level, say). Those have neither subset present, so drop exactly those: the result + # then carries one row per cell the evaluation actually computed. + wide = wide[wide["full"].notna() | wide["common"].notna()].reset_index(drop=True) wide["cov_delta"] = wide["full"] - wide["common"] return wide[[*keys, "cov_delta"]] diff --git a/src/predval/report.py b/src/predval/report.py index b51d605..6ceac8c 100644 --- a/src/predval/report.py +++ b/src/predval/report.py @@ -230,9 +230,18 @@ def _primary_rows(metrics: pd.DataFrame, models: list[str], show_cov_delta: bool & (metrics["stratum_kind"] == "overall") ] cov = coverage_delta(metrics) - cov_auroc = cov[(cov["stratum_kind"] == "overall") & (cov["metric"] == "auroc")].set_index( - "model_id" - )["cov_delta"] + # The cell this column annotates is the primary table's own: overall stratum, no subgroup, + # threshold-free AUROC. Pinning all five keys is what makes the lookup one row per model; + # selecting on stratum_kind and metric alone also catches the per-subgroup AUROC cells and + # yields a Series per model. + cov_auroc = cov[ + (cov["stratum_kind"] == "overall") + & (cov["metric"] == "auroc") + & cov["subgroup_name"].isna() + & cov["subgroup_level"].isna() + & cov["threshold"].isna() + ] + cov_auroc = dict(zip(cov_auroc["model_id"], cov_auroc["cov_delta"], strict=True)) out = [] for model in models: @@ -254,7 +263,9 @@ def cell(metric: str, method: str, sub: pd.DataFrame = sub) -> str: } if show_cov_delta: d = cov_auroc.get(model, float("nan")) - row["cov_delta"] = _f(d, 4) if np.isfinite(d) else "0.0000" + # "n/a", not "0.0000": a missing delta means one of the two subsets had no AUROC to + # compare, which is not the same statement as "restriction changed nothing". + row["cov_delta"] = _f(d, 4) if np.isfinite(d) else "n/a" out.append(row) return out diff --git a/tests/test_report.py b/tests/test_report.py index 5e9ff43..5d0b71a 100644 --- a/tests/test_report.py +++ b/tests/test_report.py @@ -40,6 +40,7 @@ def _build( clustering_name=_UNSET, single_model=False, ensemble_members=None, + drop_for_weak=0, ): """A small clustered cohort: 2 models (or 1, if single_model), a 2-level subgroup (4 clusters each -> gated).""" @@ -70,8 +71,13 @@ def _build( } ) - pd.DataFrame(rows).to_parquet(tmp_path / "cohort.parquet", index=False) - pd.DataFrame(preds).to_parquet(tmp_path / "predictions.parquet", index=False) + cohort_df = pd.DataFrame(rows) + pred_df = pd.DataFrame(preds) + if drop_for_weak: + victims = set(cohort_df["subject_id"].head(drop_for_weak)) + pred_df = pred_df[~((pred_df["model_id"] == "weak") & pred_df["subject_id"].isin(victims))] + cohort_df.to_parquet(tmp_path / "cohort.parquet", index=False) + pred_df.to_parquet(tmp_path / "predictions.parquet", index=False) clustering = {"field": "wsi"} if clustering_name is not _UNSET and clustering_name is not None: clustering["name"] = clustering_name @@ -233,6 +239,39 @@ def test_provenance_carries_hashes_seed_and_versions(html, evaluation) -> None: # ----------------------------------------------------------------------------- determinism +def test_partial_coverage_renders_a_per_model_cov_delta_column(tmp_path: Path) -> None: + """The cov_delta column only renders when coverage differs across the roster, so every + uniform-coverage fixture leaves it untested. One model scoring a subset is what exercises it: + the column must appear once per model, as a single number, not a Series repr.""" + ev = _build(tmp_path, drop_for_weak=int(N_CLUSTERS * PER_CLUSTER * 0.2)) + ctx = build_context(ev.metrics, ev.fragility, ev.coverage, ev.manifest) + assert not ctx["identical_coverage"], "the fixture must have non-uniform coverage" + rows = {r["model"]: r["cov_delta"] for r in ctx["primary_rows"]} + assert set(rows) == {"good", "weak"} + for model, cell in rows.items(): + assert isinstance(cell, str), model + assert "\n" not in cell and "dtype" not in cell, f"{model}: {cell!r} is not a scalar cell" + assert cell == "n/a" or float(cell) == pytest.approx(float(cell)) + html = render_evaluation(ev) + assert "AUROC cov-delta" in html + + +def test_coverage_delta_has_one_row_per_computed_cell(tmp_path: Path) -> None: + """The pivot that builds cov_delta fills the full cross product of its index levels, which + invents cells no stratum ever produced (overall crossed with a subgroup level). Those are + dropped, so the keys stay unique per computed cell.""" + from predval import coverage_delta + + ev = _build(tmp_path, drop_for_weak=int(N_CLUSTERS * PER_CLUSTER * 0.2)) + cov = coverage_delta(ev.metrics) + keys = ["model_id", "stratum_kind", "subgroup_name", "subgroup_level", "metric", "threshold"] + assert not cov.duplicated(keys, keep=False).any() + overall = cov[(cov["stratum_kind"] == "overall") & cov["subgroup_name"].isna()] + assert len(overall), "the overall stratum must survive the drop" + invented = cov[(cov["stratum_kind"] == "overall") & cov["subgroup_level"].notna()] + assert not len(invented), "overall must never carry a subgroup level" + + def test_render_is_deterministic(evaluation) -> None: a = render_evaluation(evaluation) b = render_evaluation(evaluation) From e2e23a871d207d284ba6f06b1c17a9d952008f54 Mon Sep 17 00:00:00 2001 From: aRealGem <203440826+aRealGem@users.noreply.github.com> Date: Sun, 4 Oct 2026 23:41:56 +0000 Subject: [PATCH 2/2] examples/improve_panel: state the admissible-rung skill result as the B=2000 run found it Three of five models are reliably negative at the admissible rung and a fourth straddles zero; the docstring had said four of five were negative. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0138yQfFhvnVRHd1XhTEz54w --- examples/improve_panel/skill_ladder.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/improve_panel/skill_ladder.py b/examples/improve_panel/skill_ladder.py index dbcb9c5..9e1b64a 100644 --- a/examples/improve_panel/skill_ladder.py +++ b/examples/improve_panel/skill_ladder.py @@ -6,9 +6,10 @@ predval's report quotes one Brier skill score per model, at the *admissible* rung -- the lowest rung whose cross-fitted paired gain interval excludes zero. That rule is parsimony-first by design: do not buy a more complex correction than you can prove you need. It also means the -headline number depends on a selection rule, and on this roster four of five models come out with -NEGATIVE skill, i.e. worse than quoting the cohort base rate to every patient. A reviewer's first -objection is the right one: maybe a higher rung would have rescued them and the rule hid it. +headline number depends on a selection rule, and on this roster three of five models come out +there with reliably NEGATIVE skill, i.e. worse than quoting the cohort base rate to every +patient, and a fourth is indistinguishable from zero. A reviewer's first objection is the right +one: maybe a higher rung would have rescued them and the rule hid it. So this sidecar drops the selection rule and scores every rung the ladder could fit: