diff --git a/docs/spec.md b/docs/spec.md index 2299263..75fcd1d 100644 --- a/docs/spec.md +++ b/docs/spec.md @@ -320,11 +320,28 @@ makes the interval tight, because a score and its recalibration move together wi is the interval the pitch shows in section 3; it replaces the earlier `(interval pending §4.5)` marker. The gain is emitted wherever the ladder runs (overall, and gated-pass subgroups). -**What is still deferred** is the interval on a recalibrated metric *level* that refits the -correction inside every bootstrap replicate — and, for the cross-fitted rungs, inside every -replicate *and* every fold — so that it carries the variance of the correction itself. The paired -gain holds the cross-fit mapping fixed; the full refit-in-replicate interval is a design surface of -its own and remains backlog rather than a silent omission. +**The fixed-mapping interval under-covers, and a refit-in-replicate interval is available.** The +paired interval above holds the cross-fit mapping fixed, so it carries none of the variance of +fitting the correction. In simulation (clustered logit random effect, true slope 0.7 or 0.5, +G = 10–50 clusters of 5, `forecastbench/sim_refit_origin.py` in the project record) its empirical +coverage at a nominal 0.95 was 0.74–0.93, with the misses almost all on the high side: the +cross-fitted gain is a little pessimistic relative to the gain a full-sample fit delivers, and the +percentile bootstrap does not widen for that. + +`uncertainty.gain_interval: refit` selects the alternative: every cluster-bootstrap replicate +re-runs the grouped K-fold cross-fit on the resampled rows and takes *that* replicate's paired +gain (`ci_method = paired_cluster_bootstrap_refit`; the point estimate is unchanged). Fold +membership inside a replicate is keyed on the **original** cluster id, so a cluster drawn twice has +both copies held out together. That detail is load-bearing: relabelling copies as distinct clusters +lets a copy in a train fold score its twin and biases every replicate gain upward (coverage then +looks fine for the wrong reason). With original-id folds the same simulation gave coverage +0.90–0.97 in every cell, at the cost of lower power (fewer intervals exclude zero) and B × K extra +rung fits per model and stratum. + +The default stays `fixed` in this version so that frozen golden outputs remain reproducible; the +report and manifest name which interval was used. New cohorts should declare `refit`. Still +deferred: an interval on a recalibrated metric *level* (not the gain) that refits inside every +replicate and every fold. ### 4.6 Every metric at every rung, with two integrity checks (S3.1) diff --git a/src/predval/evaluate.py b/src/predval/evaluate.py index 2268207..d9a91a4 100644 --- a/src/predval/evaluate.py +++ b/src/predval/evaluate.py @@ -517,6 +517,7 @@ def evaluate(cohort: Cohort, predictions: Predictions) -> Evaluation: "seed": unc.seed, "ci_level": ci, "show_naive_ci": unc.show_naive_ci, + "gain_interval": unc.gain_interval, "clustered": clustered, "clustering_field": spec.clustering.field if spec.clustering else None, # Human-readable noun for the clustering unit ("slide", "region", ...), threaded into @@ -712,10 +713,27 @@ def row( # out-of-fold Brier, with a paired slide bootstrap (item 1). This is the interval the pitch # shows in section 3, replacing the deferred "(interval pending §4.5)" marker. Emitted for # every rung whose cross-fit mapping is available; ci_method is paired_cluster_bootstrap. + # With `uncertainty.gain_interval: refit` the correction is re-cross-fitted inside every + # replicate instead (§4.5); same point estimate, ci_method paired_cluster_bootstrap_refit. + n_folds = min(5, n_clusters) for lr_rung, cf_pred in lrep.crossfit_predictions.items(): - gain, gain_ci = U.paired_brier_gain_interval( - y, p, cf_pred, groups, ci, n_boot=unc.n_boot, seed=unc.seed - ) + if unc.gain_interval == "refit": + gain, gain_ci = U.paired_brier_gain_refit_interval( + y, + p, + cf_pred, + groups, + ci, + n_boot=unc.n_boot, + seed=unc.seed, + crossfit=lambda yy, pp, gg, _r=lr_rung: R._crossfit_predictions( + _r, yy, pp, gg, n_folds + )[0], + ) + else: + gain, gain_ci = U.paired_brier_gain_interval( + y, p, cf_pred, groups, ci, n_boot=unc.n_boot, seed=unc.seed + ) rows.append(row("paired_gain_brier", gain, gain_ci, rung=lr_rung, fit_mode=R.CROSSFIT)) # rung3 non-monotonicity diagnostics: ALWAYS recorded in the artefact (S4.1 item 2), so a diff --git a/src/predval/schema.py b/src/predval/schema.py index ccb24ad..d83041e 100644 --- a/src/predval/schema.py +++ b/src/predval/schema.py @@ -89,6 +89,13 @@ class UncertaintySpec(StrictModel): #: unit-of-analysis exhibit. It is labelled incorrect; the point is to make the cost of #: ignoring clustering concrete. show_naive_ci: bool = True + #: How the paired cross-fit Brier gain interval is built (docs/spec.md section 4.5). + #: ``fixed`` holds the out-of-fold correction fixed and resamples clusters of its per-row loss + #: difference; it is fast but under-covers (about 0.75-0.90 at a nominal 0.95 in simulation) + #: because it carries no variance from the correction itself. ``refit`` re-runs the grouped + #: cross-fit inside every cluster-bootstrap replicate, with every copy of a resampled cluster + #: held out together, and recovers nominal coverage at the price of B x K extra fits. + gain_interval: Literal["fixed", "refit"] = "fixed" class CompletenessSpec(StrictModel): diff --git a/src/predval/uncertainty.py b/src/predval/uncertainty.py index 32e389a..7630e45 100644 --- a/src/predval/uncertainty.py +++ b/src/predval/uncertainty.py @@ -12,6 +12,7 @@ from __future__ import annotations +from collections.abc import Callable from dataclasses import dataclass import numpy as np @@ -127,6 +128,51 @@ def paired_brier_gain_interval( return point, percentile_interval(samples, ci_level, "paired_cluster_bootstrap") +def paired_brier_gain_refit_interval( + y: np.ndarray, + p0: np.ndarray, + p_r: np.ndarray, + groups: np.ndarray, + ci_level: float, + *, + n_boot: int, + seed: int, + crossfit: Callable[[np.ndarray, np.ndarray, np.ndarray], np.ndarray | None], +) -> tuple[float, Interval]: + """Refit-in-replicate interval for the paired Brier gain (docs/spec.md section 4.5). + + The point estimate is the same as :func:`paired_brier_gain_interval`: the gain of the given + out-of-fold correction ``p_r`` over ``p0`` on the rows as they are. The interval differs: each + cluster-bootstrap replicate re-runs ``crossfit`` on the resampled rows and takes *that* + replicate's paired gain, so the variance of fitting the correction is carried, not conditioned + away. Fold membership inside a replicate is keyed on the ORIGINAL cluster id: a cluster drawn + twice has both copies held out together, which is what stops a copy in a train fold from scoring + its twin in the eval fold -- that leak would push the replicate gains upward and the interval + would look well-covered for the wrong reason. + + ``crossfit(y, p, groups)`` returns out-of-fold corrected probabilities for the rows it is given, + or None when the fit degenerates; such replicates are dropped (the percentile interval ignores + undefined draws). Costs B grouped cross-fits per rung; in simulation at G = 10-50 clusters this + interval covers at 0.90-0.97 where the fixed-mapping one covers at 0.74-0.93. + """ + d = (p0 - y) ** 2 - (p_r - y) ** 2 + point = float(np.mean(d)) if d.size else float("nan") + clusters = cluster_indices(groups) + if len(clusters) < 2 or not np.isfinite(point): + return point, Interval(float("nan"), float("nan"), "paired_cluster_bootstrap_refit") + rng = np.random.default_rng(seed) + samples = np.full(n_boot, np.nan) + for b in range(n_boot): + pick = rng.integers(0, len(clusters), len(clusters)) + idx = np.concatenate([clusters[c] for c in pick]) + rep_groups = np.concatenate([np.full(clusters[c].size, c) for c in pick]) + rep_pred = crossfit(y[idx], p0[idx], rep_groups) + if rep_pred is None: + continue + samples[b] = np.mean((p0[idx] - y[idx]) ** 2 - (rep_pred - y[idx]) ** 2) + return point, percentile_interval(samples, ci_level, "paired_cluster_bootstrap_refit") + + def brier_skill_interval( y: np.ndarray, p: np.ndarray, diff --git a/tests/golden/gusto/manifest.json b/tests/golden/gusto/manifest.json index 5292fea..796f71b 100644 --- a/tests/golden/gusto/manifest.json +++ b/tests/golden/gusto/manifest.json @@ -38,6 +38,7 @@ "clustered": true, "clustering_field": "regl", "clustering_name": "region", + "gain_interval": "fixed", "n_boot": 2000, "seed": 1337, "show_naive_ci": true, diff --git a/tests/golden/improve/manifest.json b/tests/golden/improve/manifest.json index 7c752cd..6683104 100644 --- a/tests/golden/improve/manifest.json +++ b/tests/golden/improve/manifest.json @@ -145,6 +145,7 @@ "clustered": true, "clustering_field": "patient", "clustering_name": "patient", + "gain_interval": "fixed", "n_boot": 2000, "seed": 1337, "show_naive_ci": true, diff --git a/tests/golden/pcam/manifest.json b/tests/golden/pcam/manifest.json index 8b7921c..39fabbd 100644 --- a/tests/golden/pcam/manifest.json +++ b/tests/golden/pcam/manifest.json @@ -119,6 +119,7 @@ "clustered": true, "clustering_field": "wsi", "clustering_name": "slide", + "gain_interval": "fixed", "n_boot": 2000, "seed": 1337, "show_naive_ci": true, diff --git a/tests/test_evaluate.py b/tests/test_evaluate.py index 1a82155..bd79909 100644 --- a/tests/test_evaluate.py +++ b/tests/test_evaluate.py @@ -523,3 +523,35 @@ def test_no_clustering_declared_is_flagged(tmp_path: Path) -> None: res = evaluate(load_cohort(tmp_path / "cohort.yaml"), preds) assert any(f.code == "no_clustering_declared" for f in res.flags) assert res.manifest["uncertainty"]["clustered"] is False + + +def test_refit_gain_interval_is_selectable_and_keeps_the_point(tmp_path) -> None: + """`uncertainty.gain_interval: refit` swaps the paired-gain ci_method, records the choice in the + manifest, and leaves every point estimate identical to the fixed-mapping run (§4.5).""" + (tmp_path / "f").mkdir() + (tmp_path / "r").mkdir() + fixed = evaluate(*build(tmp_path / "f")) + refit = evaluate( + *build( + tmp_path / "r", + spec_extra={"uncertainty": {"n_boot": 60, "seed": 1337, "gain_interval": "refit"}}, + ) + ) + assert fixed.manifest["uncertainty"]["gain_interval"] == "fixed" + assert refit.manifest["uncertainty"]["gain_interval"] == "refit" + + def gains(res): + g = res.metrics[res.metrics["metric"] == "paired_gain_brier"] + return g.set_index( + ["model_id", "subset", "stratum_kind", "subgroup_level", "rung"] + ).sort_index() + + gf, gr = gains(fixed), gains(refit) + assert (gf["ci_method"] == "paired_cluster_bootstrap").all() + assert (gr["ci_method"] == "paired_cluster_bootstrap_refit").all() + assert gf.index.equals(gr.index) + np.testing.assert_allclose(gf["value"].to_numpy(), gr["value"].to_numpy()) + # the rest of the table is untouched by the choice + rest_f = fixed.metrics[fixed.metrics["metric"] != "paired_gain_brier"].reset_index(drop=True) + rest_r = refit.metrics[refit.metrics["metric"] != "paired_gain_brier"].reset_index(drop=True) + pd.testing.assert_frame_equal(rest_f, rest_r) diff --git a/tests/test_uncertainty.py b/tests/test_uncertainty.py index f14818e..9d3a09f 100644 --- a/tests/test_uncertainty.py +++ b/tests/test_uncertainty.py @@ -202,6 +202,60 @@ def test_paired_gain_is_tighter_than_an_unpaired_difference() -> None: assert (paired.high - paired.low) < 0.5 * (unpaired.high - unpaired.low) +def _rung2_crossfit(y: np.ndarray, p: np.ndarray, groups: np.ndarray): + from predval.recalibrate import _crossfit_predictions + + return _crossfit_predictions("rung2", y, p, groups, min(5, np.unique(groups).size))[0] + + +def test_refit_gain_interval_shares_the_point_and_is_deterministic() -> None: + y, p0, _, g = _recalibrated_pair() + p_r = _rung2_crossfit(y, p0, g) + point, ci = U.paired_brier_gain_refit_interval( + y, p0, p_r, g, 0.95, n_boot=100, seed=1337, crossfit=_rung2_crossfit + ) + fixed_point, _ = U.paired_brier_gain_interval(y, p0, p_r, g, 0.95, n_boot=100, seed=1337) + assert ci.method == "paired_cluster_bootstrap_refit" + assert point == pytest.approx(fixed_point), "refit changes the interval, never the point" + assert ci.low <= ci.high and np.isfinite(ci.low) + p2, ci2 = U.paired_brier_gain_refit_interval( + y, p0, p_r, g, 0.95, n_boot=100, seed=1337, crossfit=_rung2_crossfit + ) + assert (point, ci.low, ci.high) == (p2, ci2.low, ci2.high) + + +def test_refit_gain_interval_holds_every_copy_of_a_cluster_out_together() -> None: + """A cluster drawn twice must land in one fold with both copies: the crossfit callable sees + group ids that are the ORIGINAL cluster labels, so duplicates never straddle train and eval.""" + y, p0, _, g = _recalibrated_pair(n_clusters=6, per=40) + p_r = _rung2_crossfit(y, p0, g) + seen: list[int] = [] + + def spy(yy, pp, gg): + seen.append(np.unique(gg).size) + return _rung2_crossfit(yy, pp, gg) + + U.paired_brier_gain_refit_interval(y, p0, p_r, g, 0.95, n_boot=30, seed=7, crossfit=spy) + # With 6 clusters drawn with replacement, some replicates must contain a duplicate -> fewer + # distinct ids than clusters. If copies were relabelled fresh this would always equal 6. + assert min(seen) < 6 and max(seen) <= 6 + + +def test_refit_gain_interval_drops_degenerate_replicates() -> None: + y, p0, _, g = _recalibrated_pair() + p_r = _rung2_crossfit(y, p0, g) + calls = {"n": 0} + + def flaky(yy, pp, gg): + calls["n"] += 1 + return None if calls["n"] % 2 else _rung2_crossfit(yy, pp, gg) + + point, ci = U.paired_brier_gain_refit_interval( + y, p0, p_r, g, 0.95, n_boot=40, seed=1, crossfit=flaky + ) + assert np.isfinite(point) and np.isfinite(ci.low) and np.isfinite(ci.high) + + def test_paired_gain_degrades_to_nan_with_one_cluster() -> None: y, p0, p_r, _ = _recalibrated_pair(n_clusters=1, per=200) g = np.full(y.size, "only")