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..9e1b64a --- /dev/null +++ b/examples/improve_panel/skill_ladder.py @@ -0,0 +1,180 @@ +#!/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 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: + + 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)