From 507e8df55634a4dbe76bca6e7ce4949bbb1a80e9 Mon Sep 17 00:00:00 2001 From: Frank Buckley Date: Tue, 28 Jul 2026 23:21:33 +0100 Subject: [PATCH 1/2] feat(diagnostics): close the psense refit gap and diagnose the last two holdouts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Refits the 21 fits whose traces predated the log_prior/log_likelihood wiring, so power-scaling could not be reconstructed from their draws. Coverage goes from 170 to 192 of 194 reporting fits. The refits overwrote 19 published mediation and measurement fits, so the first check was whether anything moved. Nothing did: the mediation headlines reproduce the values already quoted in #404 and the findings notes (med-066 NIE_B -0.030 words, med-074 proportion mediated 0.7%), and the three rli-mm models reproduce their documented gate failures to the digit against the 2026-07-21 findings note. That exact reproduction across a two-week gap and the dse-research-utils v0.8.0 bump is a reproducibility check the suite had not previously run. Also wires psense into fit_rlm_joint_growth, which #416 missed. That family is compute_loo=False, so the sampling stage does not attach the groups psense needs and they have to be requested explicitly. The refit shows the real obstacle is not an intractable likelihood: the model draws pm.LKJCorr("measure_corr_chol"), deliberately chosen over LKJCholeskyCov, and PyMC stores the transformed value variable as "measure_corr_chol_cholesky" while the free RV keeps the base name. compute_log_prior and compute_log_likelihood both require an exact name match, so both refuse. That is a naming seam, plausibly fixable, and reconciling it reaches beyond psense coverage — left as a follow-up. The value of the wiring is that the absence is now diagnosed with a specific error rather than silent. rlm-mm-001 remains the one true exemption (non-converged posterior). Note corrected: it had listed rlm-jc-001 as a refit gap when it was a code gap. Refs #381, #383 Co-Authored-By: Claude Opus 5 --- .../202607261700-psense-coverage-backfill.md | 32 ++++++++++++++++--- .../statistical_models/pipeline.py | 9 ++++++ 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/notes/202607261700-psense-coverage-backfill.md b/notes/202607261700-psense-coverage-backfill.md index a099edb7..61e9f616 100644 --- a/notes/202607261700-psense-coverage-backfill.md +++ b/notes/202607261700-psense-coverage-backfill.md @@ -3,9 +3,9 @@ > [!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) +# Power-scaling sensitivity is now measured for 192 of 194 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. +**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 86 → 170. A second pass the same day refitted the 21 fits whose traces predated that wiring, taking coverage to **192 of 194**; see the closing section for what that changed (nothing) and for the two that remain. ## Why this was a gap rather than a clean result @@ -43,12 +43,36 @@ These are precisely the refits #382 already proposes, and this pass supplies the 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 +## The remaining 23, resolved (2026-07-26, later the same day) -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. +Their traces predate the `log_prior` / `log_likelihood` wiring, so power-scaling could not be reconstructed from them. **Twenty-one were refitted at `--config reporting` and now carry psense; coverage is 192 of 194.** The two that remain are not refit gaps, and this section's original claim that they were is corrected below. `scripts/regenerate_psense.py` exits non-zero while any target is in this state, so a sweep cannot report success while leaving estimands unmeasured. +### The refits changed no published result + +Nineteen mediation and measurement fits were overwritten, so the first thing checked was whether anything moved. Nothing did. The mediation headlines reproduce the values already quoted in #404 and the findings notes — `med-066` NIE_B −0.030 words, `med-075` −0.033, `med-074` proportion mediated 0.7 % — and the three `rli-mm` models reproduce their **gate failures to the digit** against `notes/202607210914-findings-measurement.md`: `mm-001` R-hat 1.0195 / ESS 354 / 1 divergence, `mm-002` 1.0478 / 64 / 0, `mm-101` 1.0213 / 260 / 57. Those failures are pre-existing and documented, their structural legs already on HOLD; the exact reproduction across a two-week gap and a `dse-research-utils` v0.8.0 bump is itself a reproducibility check the suite had not previously run. + +One consequence for reading the new numbers: the three `rli-mm` fits flag ~96 % of their parameters (168/176, 170/179, 169/176), which corroborates #383's diagnosis that `HalfNormal(1)` on loadings-and-residuals puts ~32 % of prior mass on loadings > 1. But those posteriors **fail the gate**, so the same reasoning that exempts `rlm-mm-001` applies: treat this as corroboration of #383, not as an independent measurement. + +### `rlm-jc-001` was a code gap, not a refit gap — and is now a diagnosed one + +This note originally listed `rlm-jc-001` among the fits needing a refit. That was wrong. `fit_rlm_joint_growth` calls `_run_sampling_and_loo(ctx, compute_loo=False)` and never called `compute_log_likelihood_and_prior` **or** `run_psense` at all — #416 added that wiring to `fit_mediation` and `fit_correlated_factor` but not to this family — so a refit alone would have burned sampling time and produced nothing. + +The wiring is now added, and the refit shows the real obstacle, which is neither what this note assumed nor what was predicted before running it. Both groups are skipped with: + +``` +exact match required for all data variable names, but +['eta_cell', 'sigma_subject', 'measure_corr_chol_cholesky', 'z_subject', 'kappa'] != +['eta_cell', 'sigma_subject', 'kappa', 'measure_corr_chol', 'z_subject'] +``` + +The model draws `pm.LKJCorr("measure_corr_chol", ...)` — chosen deliberately over `LKJCholeskyCov` to avoid its unidentified nuisance sd scales — and PyMC stores the **transformed value variable** in the posterior as `measure_corr_chol_cholesky` while the model's free RV keeps the base name. `compute_log_prior` and `compute_log_likelihood` both require an exact name match, so both refuse. It is a naming seam, not an intractable likelihood, and it blocks the prior group as well as the likelihood group. + +That makes it plausibly fixable rather than an intrinsic exemption, and the value of adding the wiring is precisely that the absence is now **diagnosed with a specific error** instead of silent. Fixing it means reconciling the LKJCorr value-variable name across the `compute_log_*` seam, which reaches beyond psense coverage and is left as a follow-up. + +`rlm-mm-001` remains the one true exemption #381 already records: its posterior has not converged, so a sensitivity diagnostic on it would not mean anything. + ## 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). diff --git a/src/language_reading_predictors/statistical_models/pipeline.py b/src/language_reading_predictors/statistical_models/pipeline.py index aab60fe5..25444fea 100644 --- a/src/language_reading_predictors/statistical_models/pipeline.py +++ b/src/language_reading_predictors/statistical_models/pipeline.py @@ -8835,6 +8835,15 @@ def fit_rlm_joint_growth(spec: ModelSpec, config: str = "dev") -> StatisticalFit diag_vars = ["eta_cell", "sigma_subject", "kappa", "measure_corr_pairs"] section_header("Summary diagnostics") _diag.summary_diagnostics(ctx, var_names=diag_vars) + # Power-scaling prior sensitivity on the reported parameters (#381). This family + # is ``compute_loo=False``, so the groups psense needs are not attached by the + # sampling stage and have to be requested here. ``strict=False`` because this is + # the model whose LKJ-Cholesky likelihood ``pm.compute_log_likelihood`` cannot + # evaluate — the very case that argument exists for. When it cannot, the fit + # degrades to a warning and simply gets no psense, which is a *measured and + # declined* exemption rather than the silent absence it was before. + _diag.compute_log_likelihood_and_prior(ctx, strict=False) + _diag.run_psense(ctx, var_names=diag_vars) _run_ppc(ctx, var_names=[f"score_{m}" for m in measure_syms]) From 8e05647aa31e3d0a246974f8fdfec268af9a18a0 Mon Sep 17 00:00:00 2001 From: Frank Buckley Date: Fri, 31 Jul 2026 10:18:38 +0100 Subject: [PATCH 2/2] docs(notes): draft the upstream PyMC report for the LKJCorr value-variable naming bug MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds notes/assets/draft-pymc-issue-get-untransformed-name.md, the reviewed-before- posting draft of a pymc-devs/pymc bug report for the seam diagnosed in #453: get_untransformed_name strips a fixed three underscore-separated components, which mangles any transform whose name itself contains an underscore. LKJCorr's default CholeskyCorrTransform ("cholesky_corr") and LogExpM1 ("log_exp_m1") both break, making compute_log_prior/compute_log_likelihood unusable on models containing them — the reason lrp-rlm-jc-001 cannot emit psense (#381, #452). Every code block in the draft was executed against the installed PyMC 6.1.0, the error messages are verbatim output, and the v6.2.0 tag was checked byte-identical. The suggested fix direction (derive the untransformed name from node.op.name, which the same loop already passes to model_free_rv, instead of parsing value.name) was verified empirically on 6.1.0. Two allow-list additions (myvar, xarray) so the draft passes spellcheck. Refs #453 Co-Authored-By: Claude Fable 5 --- config/spellcheck/allow-en.txt | 2 + ...draft-pymc-issue-get-untransformed-name.md | 123 ++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 notes/assets/draft-pymc-issue-get-untransformed-name.md diff --git a/config/spellcheck/allow-en.txt b/config/spellcheck/allow-en.txt index 7e17c642..a56902b8 100644 --- a/config/spellcheck/allow-en.txt +++ b/config/spellcheck/allow-en.txt @@ -474,6 +474,7 @@ mumedu mumedupost mumocc mypy +myvar ncol ndarray networkx @@ -620,6 +621,7 @@ woco workstream worktree worktrees +xarray xcheck yarcewr yarclet diff --git a/notes/assets/draft-pymc-issue-get-untransformed-name.md b/notes/assets/draft-pymc-issue-get-untransformed-name.md new file mode 100644 index 00000000..ade814bc --- /dev/null +++ b/notes/assets/draft-pymc-issue-get-untransformed-name.md @@ -0,0 +1,123 @@ + + + +**Title:** `compute_log_prior` / `compute_log_likelihood` fail for any transform whose name contains an underscore (`LKJCorr`, `LogExpM1`) + +--- + +> [!NOTE] +> Drafted by an LLM-based AI tool (Claude Code/Opus 5, revised by Claude Code/Fable 5). + +### Describe the issue + +`get_untransformed_name` recovers a variable's untransformed name by dropping a **fixed three** trailing underscore-separated components — one for the transform name and two for the `__` marker: + +```python +# pymc/util.py +def get_transformed_name(name, transform): + return f"{name}_{transform.name}__" + +def get_untransformed_name(name): + if not is_transformed_name(name): + raise ValueError(f"{name} does not appear to be a transformed name") + return "_".join(name.split("_")[:-3]) +``` + +That is only correct when `transform.name` itself contains no underscore. Two shipped transforms violate this: + +| Transform class | `name` | value var | `get_untransformed_name` | | +| --------------------------------------------------- | --------------- | ---------------------- | ------------------------ | --- | +| `CholeskyCorrTransform` (the default for `LKJCorr`) | `cholesky_corr` | `corr_cholesky_corr__` | `corr_cholesky` | ❌ | +| `LogExpM1` | `log_exp_m1` | `x_log_exp_m1__` | `x_log_exp` | ❌ | + +Every other shipped transform round-trips correctly (`log`, `logodds`, `ordered`, `sumto1`, `zerosum`, `cholesky-cov`, `cholesky-cov-packed`). + +The user-visible consequence is that **`compute_log_prior` and `compute_log_likelihood` cannot be run at all** on any model containing an `LKJCorr` (or `LogExpM1`-transformed) variable. `compute_log_density` builds its elemwise function over `remove_value_transforms(model).value_vars` but subsets the posterior by `[rv.name for rv in model.free_RVs]`; `remove_value_transforms` → `change_value_transforms` renames each value variable through `get_untransformed_name(value.name)`, so the two name sets no longer agree and `xarray` raises. + +This blocks any downstream consumer that needs those groups. In our case it is power-scaling prior sensitivity (`arviz_stats.psense`), which requires both `log_prior` and `log_likelihood`, so a model using `LKJCorr` for a correlation matrix cannot be checked for prior sensitivity at all. + +### Reproducible code example + +A single-variable model with no data is enough: + +```python +import pymc as pm + +with pm.Model() as m: + pm.LKJCorr("corr", n=3, eta=2.0) + idata = pm.sample(draws=10, tune=10, chains=1, progressbar=False) + +pm.compute_log_prior(idata, model=m) +``` + +The same failure with `LogExpM1`: + +```python +import pymc as pm +from pymc.distributions import transforms as tr + +with pm.Model() as m: + pm.HalfNormal("x", 1.0, default_transform=tr.log_exp_m1) + idata = pm.sample(draws=10, tune=10, chains=1, progressbar=False) + +pm.compute_log_prior(idata, model=m) +``` + +And the underlying round-trip, without sampling: + +```python +from pymc.util import get_transformed_name, get_untransformed_name + +class T: + def __init__(self, name): + self.name = name + +get_untransformed_name(get_transformed_name("myvar", T("log"))) # 'myvar' OK +get_untransformed_name(get_transformed_name("myvar", T("cholesky-cov"))) # 'myvar' OK +get_untransformed_name(get_transformed_name("myvar", T("cholesky_corr"))) # 'myvar_cholesky' BUG +get_untransformed_name(get_transformed_name("myvar", T("log_exp_m1"))) # 'myvar_log_exp' BUG +``` + +### Error message + +``` +ValueError: exact match required for all data variable names, but ['corr_cholesky'] != ['corr']: {'corr', 'corr_cholesky'} are not in both. +``` + +and for the `LogExpM1` case: + +``` +ValueError: exact match required for all data variable names, but ['x_log_exp'] != ['x']: {'x_log_exp', 'x'} are not in both. +``` + +### PyMC version information + +``` +PyMC: 6.1.0 +PyTensor: 3.1.2 +ArviZ: 1.2.0 +xarray: 2026.7.0 +Python: 3.14.6 +OS: Darwin 25.5.0 (arm64) +Installed via conda-forge +``` + +`get_untransformed_name` and `CholeskyCorrTransform.name` are unchanged at the `v6.2.0` tag, so this is present on the current release as well as 6.1.0. + +### Context for the issue + +A possible fix, offered tentatively since I do not know the constraints on `get_untransformed_name` as public API: the one call site that causes this failure does not need to parse the name at all. In `change_value_transforms` (`pymc/model/transform/conditioning.py`) the loop recovers the untransformed name by string surgery on the value variable: + +```python +try: + untransformed_name = get_untransformed_name(value.name) +except ValueError: + untransformed_name = value.name +``` + +but the same iteration already has that name directly — a few lines later it passes `node.op.name` (the model variable's untransformed name; `rv.name` from `rv, value = node.inputs` carries it too, verified on 6.1.0) into `model_free_rv` when reconstructing the node. Deriving `untransformed_name` from `node.op.name` instead of parsing `value.name` would fix every underscore-named transform at once and remove the implicit "transform names must not contain underscores" constraint, without touching `get_untransformed_name` itself. + +If `get_untransformed_name` should also be fixed for its other callers, note that it cannot be made correct from the string alone in general — a variable genuinely named `myvar_cholesky` transformed by a hypothetical `corr` transform is indistinguishable from `myvar` transformed by `cholesky_corr` — so it would need the transform passed in, or the naming constraint documented and enforced (e.g. a check that `transform.name` contains no underscore, which would at least turn a silent mismatch into a clear error at definition time).