feat(diagnostics): backfill power-scaling sensitivity across the stored fits (#381)

config/spellcheck/allow-en.txt

@@ -548,6 +548,7 @@ residualized
 resumable
 retau
 reto
+reweighting
 rglob
 rhat
 riclpm

notes/202607261700-psense-coverage-backfill.md

@@ -0,0 +1,62 @@
+<!-- SPDX-License-Identifier: CC-BY-4.0 -->
+
+> [!NOTE]
+> Drafted by an LLM-based AI tool (Claude Code/Opus 5).
+
+# Power-scaling sensitivity is now measured for 170 of 193 reporting fits (#381)
+
+**What changed (2026-07-26, #381):** every stored `--config reporting` fit whose trace carries the `log_prior` and `log_likelihood` groups now has a `psense_summary.csv` and `psense.png`, backfilled by `scripts/regenerate_psense.py` without refitting anything. Coverage goes from 86 to 170 of 193 fits. The remaining 23 need a refit and are listed below.
+
+## Why this was a gap rather than a clean result
+
+Power-scaling sensitivity was wired into the family pipelines in #408 and #416, so the emission code has been correct for some time — `run_psense` is called from 22 fit functions covering every family. What was missing was **artefacts for fits stored before those merged**. A report with no psense table showed no flags, and a reader had no way to tell "measured, no concern" from "never measured". #381 named that distinction as its central meta-finding, and it is the whole value of this pass: it converts a large block of unverified verdicts into verified ones, in both directions.
+
+## Why no refit was needed
+
+Power-scaling is importance reweighting over the draws already in hand, not resampling. Given a trace with `log_prior` and `log_likelihood`, `arviz_stats.psense_summary` reconstructs the sensitivity after the fact, and the numbers belong to exactly the posterior the published report was written from — which is a stronger guarantee than a fresh fit would give, since a refit would produce a different (if equivalent) chain.
+
+The reported-parameter set for each backfill comes from that fit's own `diagnostics.csv` rather than from re-deriving it off the spec. That file _is_ the fit's record of which parameters it reported, so the backfill covers exactly what the report shows, and it stays correct for a model whose spec has been edited since.
+
+## What the measurement found
+
+Across 1,647 parameters in 170 fits, 46 % carry some power-scaling flag. That headline number should **not** be read as "46 % of the suite is unsound". At n ≈ 54 a weakly-informative prior against a small sample will often show sensitivity, and most of the flagged parameters are nuisance or scale terms whose prior dependence changes no reported claim. The decision-relevant question is narrower: what happens to the terms the suite actually reports as estimands?
+
+Restricting to each family's headline causal term — 155 such terms across 77 models:
+
+| Family | term         |   ✓ | prior–data conflict | strong prior / weak likelihood |
+| ------ | ------------ | --: | ------------------: | -----------------------------: |
+| `itt`  | `tau`        |  28 |                  15 |                          **0** |
+| `lf`   | `b_grp_time` |  20 |                  11 |                             13 |
+| `hs`   | `beta`       |  19 |                  17 |                              2 |
+| `did`  | `tau_t2`     |   4 |                   4 |                              3 |
+| `mech` | `beta_mech`  |  10 |                   5 |                              0 |
+| `bx`   | `delta`      |   2 |                   2 |                              0 |
+
+**The ITT suite — the primary deliverable — has no headline `tau` in the serious class.** Fifteen show potential prior–data conflict, which at this sample size is unsurprising and informational; none is prior-driven. That is a genuinely reassuring result that was previously unmeasured for a third of the suite.
+
+**Five models have their one randomised causal quantity flagged `potential strong prior / weak likelihood`** — the class where the prior is doing more work than the data:
+
+- `lrp-rli-lf-005` and `lrp-rli-lf-006` at `b_grp_time[1]`, the clean randomised t2 contrast (phase coordinates are 0-based, so index 1 is t2). Both are floor-limited outcomes — phonetic spelling off-floor, and phoneme blending — so a thin likelihood is explicable rather than alarming, but it means the reported interval reflects the prior substantially and should be read that way.
+- `lrp-rli-did-001` (word reading), `lrp-rli-did-003` (phoneme blending) and `lrp-rli-did-013` at `tau_t2`, the clean randomised t2 contrast.
+
+These are precisely the refits #382 already proposes, and this pass supplies the measurement that motivates them rather than assuming it.
+
+Two other results worth recording. `lrp-rli-hs-001`'s slab scale `hs_c2` comes back at prior 1.12 against likelihood 0.028 — `potential strong prior / weak likelihood`, exactly the dependence #381 predicted for horseshoe rankings ("the reported P(|β| > 0.1) is a direct function of `tau0`/`slab_scale`"). And `rli-lcf` and `rli-surv` flag on 100 % of their measured parameters, on small parameter sets.
+
+## Still uncovered: 23 fits needing a refit
+
+Their traces predate the `log_prior` / `log_likelihood` wiring, so power-scaling cannot be reconstructed from them: the `med` family (18 of 19), all three `rli-mm`, `rlm-jc-001` and `rlm-mm-001`. `rlm-mm-001` is the one legitimate exemption #381 already records — its posterior has not converged, so a sensitivity diagnostic on it would not mean anything. The rest need a `--config reporting` refit, at which point the fit-time wiring emits psense on its own with no further change.
+
+`scripts/regenerate_psense.py` exits non-zero while any target is in this state, so a sweep cannot report success while leaving estimands unmeasured.
+
+## Not addressed here
+
+#381 also asks for `prior_pushforward.csv` across every estimand-reporting family, and for an indicator-scale prior-predictive check on the measurement/CFA families. That is **not** the same shape of job. `prior_pushforward` pushes a treatment term through the items-scale average marginal effect and is ITT-shaped by construction — it takes the arm indicator and a trial count. Families whose estimand is a curve (`mech`), a ranking (`hs`), a latent change (`lcsm`) or a correlation (`lcf`) have no single "estimand scale" to push a prior onto, so extending it needs a per-family decision about what the estimand-scale check even is, not a call-site added. 51 of 193 fits currently carry one. The outcome-scale prior-predictive check, by contrast, is near-universal already (169 of 193).
+
+## Reproducing
+
+```bash
+python scripts/regenerate_psense.py all --dry-run
+```
+
+Drop `--dry-run` to write. `--force` recomputes where a summary already exists.

scripts/regenerate_psense.py

@@ -0,0 +1,205 @@
+#!/usr/bin/env python
+# Copyright (c) 2026 Down Syndrome Education International and contributors
+# SPDX-License-Identifier: AGPL-3.0-or-later
+
+"""Backfill ``psense_summary.csv`` / ``psense.png`` over existing fit output dirs.
+
+Power-scaling prior sensitivity (#381) is emitted at fit time by every family
+pipeline, but that wiring landed in #408 / #416 — so every fit stored *before*
+those merged carries no psense artefacts, and its report shows no flags because
+none were **measured**, not because the estimand was measured clean. That
+distinction matters most for exactly the deliverables #381 names as the most
+prior-dependent: horseshoe rankings, HSGP mechanism curves, the EiV slope.
+
+Power-scaling is importance reweighting over the draws already in hand, **not** a
+refit, so the gap closes without resampling anything: a stored trace carrying the
+``log_prior`` and ``log_likelihood`` groups can be measured after the fact and the
+numbers belong to the same posterior the published report was written from.
+
+The reported-parameter set comes from each fit's own ``diagnostics.csv`` rather
+than from re-deriving it off the spec. That file *is* the fit's record of what it
+reported, so the backfill covers exactly the parameters the report shows, and it
+stays correct for a fit whose spec has since been edited.
+
+Targets mirror ``regenerate_key_findings.py``:
+
+    regenerate_psense.py all                        # every fit output dir
+    regenerate_psense.py lrp-rli-hs-001             # one model, all its -<config> dirs
+    regenerate_psense.py lrp-rli-hs-001-reporting   # one specific fit dir
+
+Fits whose trace predates the ``log_prior`` wiring cannot be repaired here and are
+reported as skips with that reason — they need a refit, which is a separate call.
+"""
+
+from __future__ import annotations
+
+import argparse
+import warnings
+from pathlib import Path
+
+import pandas as pd
+from rich.console import Console
+
+from language_reading_predictors import paths as _paths
+from language_reading_predictors.statistical_models.diagnostics import (
+    psense_artifacts,
+)
+
+_console = Console()
+
+
+def _subdirs(root: Path) -> list[Path]:
+    """Published fit directories, excluding in-flight output transactions.
+
+    ``StatisticalFitContext.reset_output_dir`` stages each run in a *hidden* sibling
+    (``.<id>-<config>.staging-XXXX``) and promotes it only on success, so a dotted
+    directory is either a run in progress or an abandoned one. Backfilling into it
+    writes artefacts that are about to be discarded — or, worse, races a live fit.
+    """
+    if not root.is_dir():
+        return []
+    return sorted(
+        d for d in root.iterdir() if d.is_dir() and not d.name.startswith(".")
+    )
+
+
+def resolve_targets(target: str) -> list[Path]:
+    """Fit output dirs for the requested target (statistical models only)."""
+    root = _paths.stat_models_dir()
+    if target == "all":
+        return _subdirs(root)
+    return [
+        d for d in _subdirs(root) if d.name == target or d.name.startswith(f"{target}-")
+    ]
+
+
+def reported_var_names(fit_dir: Path) -> list[str]:
+    """Base parameter names the fit reported, in the order it reported them.
+
+    ``diagnostics.csv`` carries one row per *element* (``beta[L]``, ``beta[R]``,
+    ...); ``psense_summary`` takes the base name and expands it itself. Order is
+    preserved because it drives the psense figure's row layout.
+    """
+    path = fit_dir / "diagnostics.csv"
+    if not path.exists():
+        return []
+    index = pd.read_csv(path, index_col=0).index
+    seen: dict[str, None] = {}
+    for label in index:
+        seen.setdefault(str(label).split("[")[0].strip(), None)
+    return list(seen)
+
+
+def _trace_groups(idata) -> set[str]:
+    return {g.rstrip("/").split("/")[-1] for g in idata.groups}
+
+
+def backfill(fit_dir: Path, *, force: bool, dry_run: bool) -> tuple[str, str]:
+    """Return ``(status, detail)`` for one fit directory."""
+    summary = fit_dir / "psense_summary.csv"
+    if summary.exists() and not force:
+        return "present", "psense_summary.csv already written"
+    if not (fit_dir / "trace.nc").exists():
+        return "skipped", "no trace.nc"
+
+    var_names = reported_var_names(fit_dir)
+    if not var_names:
+        return "skipped", "no diagnostics.csv, so no reported-parameter set to measure"
+
+    import arviz as az
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        idata = az.from_netcdf(fit_dir / "trace.nc")
+    groups = _trace_groups(idata)
+    missing = {"log_prior", "log_likelihood"} - groups
+    if missing:
+        return (
+            "needs refit",
+            f"trace predates the {'/'.join(sorted(missing))} wiring; "
+            "power-scaling cannot be reconstructed without resampling",
+        )
+
+    # Only measure what the posterior actually carries: a spec edit since the fit can
+    # leave a name in diagnostics.csv that this trace never sampled, and psense would
+    # fail the whole model over one absent term.
+    present = set(az.extract(idata, group="posterior").data_vars)
+    measurable = [n for n in var_names if n in present]
+    dropped = [n for n in var_names if n not in present]
+    if not measurable:
+        return "skipped", f"none of {var_names} are in the posterior group"
+
+    if dry_run:
+        detail = f"would measure {len(measurable)} parameters"
+        if dropped:
+            detail += f" (not in posterior: {', '.join(dropped)})"
+        return "would write", detail
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        frame = psense_artifacts(idata, str(fit_dir), measurable)
+    if frame is None:
+        return "failed", "psense_summary could not be computed (see warning above)"
+
+    detail = f"{len(frame)} rows"
+    if "diagnosis" in frame.columns:
+        flagged = int((frame["diagnosis"].astype(str).str.strip() != "✓").sum())
+        detail += f", {flagged} flagged"
+    if dropped:
+        detail += f" (not in posterior: {', '.join(dropped)})"
+    return "written", detail
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "target", help="'all', a model id, or a fit dir name (<id>-<config>)"
+    )
+    parser.add_argument(
+        "--force",
+        action="store_true",
+        help="Recompute even where psense_summary.csv already exists",
+    )
+    parser.add_argument(
+        "--dry-run",
+        action="store_true",
+        help="Report what would be written without touching any file",
+    )
+    parser.add_argument(
+        "--output-dir",
+        default=None,
+        help="Output root override (takes precedence over DSE_LRP_OUTPUT_DIR)",
+    )
+    args = parser.parse_args()
+    if args.output_dir:
+        _paths.set_output_root(args.output_dir)
+    _console.print(f"Output root: {_paths.describe_output_root()}")
+
+    targets = resolve_targets(args.target)
+    if not targets:
+        raise SystemExit(f"No fit output directories matched {args.target!r}.")
+
+    tally: dict[str, int] = {}
+    for fit_dir in targets:
+        status, detail = backfill(fit_dir, force=args.force, dry_run=args.dry_run)
+        tally[status] = tally.get(status, 0) + 1
+        colour = {
+            "written": "green",
+            "would write": "cyan",
+            "present": "dim",
+            "skipped": "yellow",
+            "needs refit": "yellow",
+            "failed": "red",
+        }.get(status, "white")
+        _console.print(f"[{colour}]{status:11}[/{colour}] {fit_dir.name}: {detail}")
+
+    _console.print()
+    _console.print(", ".join(f"{k}: {v}" for k, v in sorted(tally.items())))
+    # A fit that needs a refit is a real coverage gap, not a no-op: surface it in the
+    # exit status so a sweep cannot look clean while leaving estimands unmeasured.
+    if tally.get("failed") or tally.get("needs refit"):
+        raise SystemExit(1)
+
+
+if __name__ == "__main__":
+    main()

src/language_reading_predictors/statistical_models/diagnostics.py

@@ -856,22 +856,26 @@ def _psense_layout(trace, var_names: list[str]) -> tuple[dict, dict]:
     return {"figure_kwargs": {"figsize": (11.0, height)}}, rc
 
 
-def run_psense(
-    context: StatisticalFitContext,
-    *,
+def psense_artifacts(
+    trace,
+    out: str,
     var_names: list[str],
-) -> None:
-    """Power-scaling prior/likelihood sensitivity (issue #125 Area 1, secondary).
-
-    Writes ``psense_summary.csv`` and ``psense.png`` for the named parameters
-    (usually the causal term). Requires the ``log_prior`` and ``log_likelihood``
-    groups added by :func:`compute_log_likelihood_and_loo`. Guarded — a missing
-    group or an API mismatch degrades to a warning (psense is recommended-but-
-    secondary at this n). Kallioinen et al. 2024.
+):
+    """Write ``psense_summary.csv`` and ``psense.png`` for ``var_names``.
+
+    Split out of :func:`run_psense` (#381) so the regeneration script can produce
+    exactly the fit-time artefacts from a stored trace, with no second
+    implementation to drift. Power-scaling is importance reweighting over the draws
+    already in hand, **not** a refit, so any fit whose trace carries the
+    ``log_prior`` and ``log_likelihood`` groups can be measured after the fact.
+
+    Returns the summary frame, or ``None`` when psense could not be computed (a
+    missing group, an API mismatch). The caller decides whether that is fatal: at
+    fit time it is a warning, because psense is recommended-but-secondary at this
+    n; a regeneration run reports it as a skip with its reason.
     """
-    out = context.output_dir
     summary_path = os.path.join(out, "psense_summary.csv")
-    context.tables.pop("psense_summary", None)
+    df = None
     try:
         os.unlink(summary_path)
     except FileNotFoundError:
@@ -884,7 +888,7 @@ def run_psense(
         # the FULL trace (thin_for_plots' own contract: "numeric summaries always
         # use the full trace"); the thinned view is only for the figure below
         # (issue #270 item 2).
-        s = azs.psense_summary(context.trace, var_names=var_names)
+        s = azs.psense_summary(trace, var_names=var_names)
         if hasattr(s, "to_dataframe"):
             df = s.to_dataframe()
         else:
@@ -900,8 +904,8 @@ def run_psense(
         df.to_csv(temporary_path)
         os.replace(temporary_path, summary_path)
         temporary_path = None
-        context.tables["psense_summary"] = df
     except Exception as exc:  # pragma: no cover
+        df = None
         if temporary_path is not None:
             try:
                 os.unlink(temporary_path)
@@ -911,7 +915,7 @@ def run_psense(
 
     import arviz_plots as azp
 
-    tr, plot_var_names = _psense_plot_view(thin_for_plots(context.trace), var_names)
+    tr, plot_var_names = _psense_plot_view(thin_for_plots(trace), var_names)
     plot_kwargs, rc = _psense_layout(tr, plot_var_names)
     with az.rc_context(rc):
         _save_pc(
@@ -920,6 +924,26 @@ def run_psense(
             "psense.png",
             title="Prior/likelihood power-scaling sensitivity",
         )
+    return df
+
+
+def run_psense(
+    context: StatisticalFitContext,
+    *,
+    var_names: list[str],
+) -> None:
+    """Power-scaling prior/likelihood sensitivity (issue #125 Area 1, secondary).
+
+    Writes ``psense_summary.csv`` and ``psense.png`` for the named parameters
+    (usually the causal term). Requires the ``log_prior`` and ``log_likelihood``
+    groups added by :func:`compute_log_likelihood_and_loo`. Guarded — a missing
+    group or an API mismatch degrades to a warning (psense is recommended-but-
+    secondary at this n). Kallioinen et al. 2024.
+    """
+    context.tables.pop("psense_summary", None)
+    df = psense_artifacts(context.trace, context.output_dir, var_names)
+    if df is not None:
+        context.tables["psense_summary"] = df
 
 
 def sample_posterior_predictive(