From 6c26b9a8123f8d59500759f4cfb2010f0ae9b39e Mon Sep 17 00:00:00 2001 From: Ethan Buckley Date: Mon, 27 Jul 2026 09:59:21 +0100 Subject: [PATCH] feat(level-factors): estimand-scale prior pushforward + surfaced t2 psense warning (#389 finding 3) The level family lacked the estimand-scale prior_pushforward.csv the ITT/gain families emit, and its key-findings headline stayed unqualified even when power-scaling flagged the t2 term for prior-data conflict (the warning was hidden in a collapsed prior section). - reporting.level_t2_marginal_effect gains group="posterior"|"prior" so the same t2 net-out transform runs on either group (byte-identical default); - reporting.level_prior_pushforward pushes the prior on the t2 group contrast through that transform, returning the same items-scale schema as prior_pushforward; fit_level_factors writes prior_pushforward.csv beside the posterior ROPE card (the level model has a single t2 term, so the pushforward is well-defined here, unlike the curve / matrix families #381 could not cover); - _kf_build_level_factors surfaces a caution bullet right after the headline when psense_summary.csv flags b_grp_time[1], via a new _kf_psense_diagnosis reader (graceful no-op when psense is absent or clear). Tests: level_prior_pushforward reads the prior group and returns the schema; the KF warning appears when the t2 term is flagged, is absent when clear, and the golden level case is unchanged. ruff / format:check / spellcheck clean; lf-001 dev fit emits prior_pushforward.csv. Addresses #389 finding 3 (bar the flagged-outcome prior-grid sensitivity, which is refit-gated). Co-Authored-By: Claude Opus 4.8 --- .../statistical_models/pipeline.py | 17 ++++ .../statistical_models/reporting.py | 81 ++++++++++++++++++- tests/statistical_models/test_key_findings.py | 40 +++++++++ tests/statistical_models/test_reporting.py | 48 +++++++++++ 4 files changed, 185 insertions(+), 1 deletion(-) diff --git a/src/language_reading_predictors/statistical_models/pipeline.py b/src/language_reading_predictors/statistical_models/pipeline.py index 854cdbd8..7eea9454 100644 --- a/src/language_reading_predictors/statistical_models/pipeline.py +++ b/src/language_reading_predictors/statistical_models/pipeline.py @@ -5877,6 +5877,23 @@ def fit_level_factors(spec: ModelSpec, config: str = "dev") -> StatisticalFitCon f"{spec.outcome_symbol}, delta={delta_prob:g})" ) items = ame_prob * n_marg + # Estimand-scale prior pushforward for the t2 term (#389 finding 3): the + # prior-predictive counterpart of this card, pushed through the same t2 + # net-out transform, so the level family emits the check the ITT / gain + # families do (it previously lacked one). Same n_trials scaling as the AME. + pf = _report.level_prior_pushforward( + ctx.trace, + phase=built.prepared.phase, + G=built.prepared.G, + n_trials=n_marg, + ability=ability, + ci_prob=ctx.reporting.ci_prob, + ) + pd.DataFrame([pf]).to_csv( + os.path.join(ctx.output_dir, "prior_pushforward.csv"), index=False + ) + ctx.tables["prior_pushforward"] = pd.DataFrame([pf]) + meta_extra["prior_pushforward"] = pf rope_s = _report.rope_card( contrast_draws, items, delta=delta, ci_prob=ctx.reporting.ci_prob ) diff --git a/src/language_reading_predictors/statistical_models/reporting.py b/src/language_reading_predictors/statistical_models/reporting.py index f87ed48a..e783c516 100644 --- a/src/language_reading_predictors/statistical_models/reporting.py +++ b/src/language_reading_predictors/statistical_models/reporting.py @@ -2961,6 +2961,7 @@ def level_t2_marginal_effect( interaction_term: str = "gamma_grp_ability", ability: np.ndarray | None = None, eta_name: str = "eta", + group: str = "posterior", ) -> tuple[np.ndarray, np.ndarray]: """Per-draw t2 randomised contrast and its items-scale AME (LRPLF, #127). @@ -2988,7 +2989,10 @@ def level_t2_marginal_effect( standardised ability covariate aligned with ``eta``'s ``obs_id`` axis (pass ``None`` when the model has no group×ability term). """ - posterior = trace.posterior + # ``group`` selects the posterior (the estimate) or the prior (the estimand-scale + # prior-predictive pushforward, #389 finding 3); both carry eta / contrast / + # interaction, so the same net-out transform applies to either. + posterior = getattr(trace, group) phase = np.asarray(phase) G = np.asarray(G, dtype=float) mask = phase == t2_phase @@ -3036,6 +3040,45 @@ def level_t2_marginal_effect( return contrast_draws, ame_prob +def level_prior_pushforward( + trace: xr.DataTree, + *, + phase: np.ndarray, + G: np.ndarray, + n_trials: int, + ability: np.ndarray | None = None, + ci_prob: float = 0.95, +) -> dict[str, float]: + """Push the **prior** on the t2 group contrast through the items-scale AME (#389 finding 3). + + The estimand-scale prior-predictive check for the level family, the counterpart + of :func:`prior_pushforward` for the ITT/gain families: before seeing data, what + does the prior on the t2 group term imply for the items-scale average marginal + effect? It runs :func:`level_t2_marginal_effect` on the persisted ``prior`` group, + so the prior is pushed through the *same* t2 net-out transform as the posterior + estimate (the level model's per-timepoint group vector + group x ability term + means the plain ``eta - term*G`` ITT pushforward does not apply). Returns the same + schema as :func:`prior_pushforward` so ``config.json`` and any consumer treat the + two families' prior checks uniformly. + """ + contrast_draws, ame_prob = level_t2_marginal_effect( + trace, phase=phase, G=G, ability=ability, group="prior" + ) + items = ame_prob * float(n_trials) + lo_q, hi_q = (1 - ci_prob) / 2, 1 - (1 - ci_prob) / 2 + return { + "prior_logit_median": float(np.median(contrast_draws)), + "prior_logit_lo": float(np.quantile(contrast_draws, lo_q)), + "prior_logit_hi": float(np.quantile(contrast_draws, hi_q)), + "prior_items_median": float(np.median(items)), + "prior_items_lo50": float(np.quantile(items, 0.25)), + "prior_items_hi50": float(np.quantile(items, 0.75)), + "prior_items_lo": float(np.quantile(items, lo_q)), + "prior_items_hi": float(np.quantile(items, hi_q)), + "n_trials": int(n_trials), + } + + def horseshoe_ranking(trace: xr.DataTree, *, delta: float = 0.1) -> pd.DataFrame: """Per-predictor ranking from a horseshoe fit (LRPHS, #116 Phase E). @@ -3286,6 +3329,28 @@ def _kf_csv_row(output_dir, name: str) -> dict | None: return df.iloc[0].to_dict() +def _kf_psense_diagnosis(output_dir, term: str) -> str | None: + """Power-scaling diagnosis for ``term`` from ``psense_summary.csv`` (#389 finding 3). + + Returns the ``diagnosis`` string (e.g. "potential prior-data conflict") when the + parameter is flagged, or ``None`` when the file or row is absent or the parameter + is clear — so a caller can surface a warning beside the headline without breaking + fits that never ran power-scaling.""" + path = os.path.join(str(output_dir), "psense_summary.csv") + if not os.path.exists(path): + return None + try: + df = pd.read_csv(path, index_col=0) + except Exception: + return None + if "diagnosis" not in df.columns or term not in df.index: + return None + diag = str(df.loc[term, "diagnosis"]).strip() + if not diag or diag.lower() in {"-", "nan", "none"}: + return None + return diag + + def _kf_csv(output_dir, name: str) -> pd.DataFrame | None: """Read one fit CSV, returning ``None`` when it is absent or empty.""" path = os.path.join(str(output_dir), name) @@ -3620,6 +3685,20 @@ def _kf_build_level_factors(output_dir, config: Mapping) -> list[dict[str, str]] rope, outcome_label, "at the end of the randomised period (t2)" ) sentences.append(_kf_sentence(headline, "headline")) + # Surface the t2 power-scaling verdict beside the headline (#389 finding 3): the + # level headline was previously unqualified even when psense flagged the t2 term + # for prior-data conflict, with the warning hidden in a collapsed prior section. + _t2_psense = _kf_psense_diagnosis(output_dir, "b_grp_time[1]") + if _t2_psense: + sentences.append( + _kf_sentence( + "Caution: a power-scaling check flags this t2 estimate as " + f"prior-sensitive ({_t2_psense}) — the prior meaningfully shapes it, " + "so read the headline alongside the prior-sensitivity section rather " + "than as a purely data-driven result.", + "warning", + ) + ) sentences.append( _kf_sentence(_kf_direction_words(rope["pd"], is_rd=is_rd), "confidence") ) diff --git a/tests/statistical_models/test_key_findings.py b/tests/statistical_models/test_key_findings.py index 8377691c..d3bf1859 100644 --- a/tests/statistical_models/test_key_findings.py +++ b/tests/statistical_models/test_key_findings.py @@ -489,6 +489,46 @@ def test_level_factors_golden_sentences(tmp_path): assert "at the end of the randomised period (t2)" in texts[0] assert "Only this t2 comparison is randomised" in texts[3] assert "crossed over" in texts[3] + # No psense_summary.csv → no caution bullet (base case is four sentences). + assert [s["kind"] for s in payload["sentences"]] == [ + "headline", "confidence", "rope", "causal" + ] + + +def test_level_factors_surfaces_t2_psense_warning(tmp_path): + # #389 finding 3: when power-scaling flags the t2 term (b_grp_time[1]), a caution + # bullet is surfaced right after the headline instead of being hidden in a + # collapsed prior section. psense_summary.csv is index-keyed by parameter. + d = _setup_dir(tmp_path, "level_factors") + _write_csv(d, "rope_summary.csv", _rope_row()) + pd.DataFrame( + { + "prior": [0.083], + "likelihood": [0.053], + "diagnosis": ["potential prior-data conflict"], + }, + index=["b_grp_time[1]"], + ).to_csv(d / "psense_summary.csv") + payload = generate_key_findings(d) + assert payload["status"] == "ok" + assert [s["kind"] for s in payload["sentences"]] == [ + "headline", "warning", "confidence", "rope", "causal" + ] + warn = payload["sentences"][1]["text"] + assert "prior-sensitive" in warn + assert "potential prior-data conflict" in warn + + +def test_level_factors_no_warning_when_t2_psense_clear(tmp_path): + # A clear ("-") diagnosis for the t2 term produces no caution bullet. + d = _setup_dir(tmp_path, "level_factors") + _write_csv(d, "rope_summary.csv", _rope_row()) + pd.DataFrame( + {"prior": [0.01], "likelihood": [0.02], "diagnosis": ["-"]}, + index=["b_grp_time[1]"], + ).to_csv(d / "psense_summary.csv") + payload = generate_key_findings(d) + assert "warning" not in [s["kind"] for s in payload["sentences"]] def test_did_golden_sentences(tmp_path): diff --git a/tests/statistical_models/test_reporting.py b/tests/statistical_models/test_reporting.py index ad620024..fd5136c9 100644 --- a/tests/statistical_models/test_reporting.py +++ b/tests/statistical_models/test_reporting.py @@ -25,6 +25,7 @@ evidence_label, favoured_direction, joint_treatment_marginals, + level_prior_pushforward, level_t2_marginal_effect, longitudinal_conditional_slopes, longitudinal_factor_correlations, @@ -1411,6 +1412,53 @@ def test_level_t2_marginal_effect_requires_t2_rows(): ) +def test_level_prior_pushforward_uses_prior_group_and_reports_items_schema(): + # #389 finding 3: the level-family estimand-scale prior pushforward runs the t2 + # net-out transform on the *prior* group (group="prior") and returns the same + # items-scale schema as the ITT/gain prior_pushforward. The logit summary must + # match b_grp_time[t2] read from the prior group (proving it reads the prior, not + # the posterior), and the items summary is that AME scaled by n_trials. + n_chain, n_draw, n_obs, n_phase = 1, 8, 6, 4 + rng = np.random.default_rng(7) + prior = xr.Dataset( + { + "eta": (("chain", "draw", "obs_id"), rng.normal(0.0, 1.0, (n_chain, n_draw, n_obs))), + "b_grp_time": (("chain", "draw", "phase"), rng.normal(0.2, 0.5, (n_chain, n_draw, n_phase))), + "gamma_grp_ability": (("chain", "draw"), rng.normal(0.0, 0.1, (n_chain, n_draw))), + }, + coords={ + "chain": np.arange(n_chain), "draw": np.arange(n_draw), + "obs_id": np.arange(n_obs), "phase": np.arange(n_phase), + }, + ) + # A distinct posterior group: if the pushforward wrongly read it, the logit + # summary below would not match the prior's b_grp_time[t2]. + posterior = prior.copy(deep=True) + posterior["b_grp_time"] = posterior["b_grp_time"] + 100.0 + trace = SimpleNamespace(prior=prior, posterior=posterior) + + phase = np.array([0, 1, 2, 3, 1, 0]) + G = np.array([1.0, 1.0, 0.0, 1.0, 0.0, 1.0]) + ability = np.array([0.5, -1.0, 0.2, 0.3, 0.8, -0.4]) + + pf = level_prior_pushforward( + trace, phase=phase, G=G, n_trials=79, ability=ability, ci_prob=0.95 + ) + assert set(pf) == { + "prior_logit_median", "prior_logit_lo", "prior_logit_hi", + "prior_items_median", "prior_items_lo50", "prior_items_hi50", + "prior_items_lo", "prior_items_hi", "n_trials", + } + assert pf["n_trials"] == 79 + # Read the prior group's b_grp_time[t2] directly and confirm the logit summary + # matches it (so the pushforward used the prior, not the +100 posterior). + contrast, _ = level_t2_marginal_effect( + trace, phase=phase, G=G, ability=ability, group="prior" + ) + assert pf["prior_logit_median"] == pytest.approx(float(np.median(contrast))) + assert pf["prior_logit_lo"] < pf["prior_logit_hi"] + + def test_proportion_at_zero_ppc(): prepared = SimpleNamespace(post_counts={"N": np.array([0.0, 0.0, 1.0, 3.0])}) # Replicated y_post (chain, draw, obs): zero-fractions 3/4 and 1/4.