uncertainty: refit-in-replicate interval for the paired cross-fit gain (spec §4.5)

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)
 

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

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

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,

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,

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,

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,

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)

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