diff --git a/src/language_reading_predictors/statistical_models/pipeline.py b/src/language_reading_predictors/statistical_models/pipeline.py index 854cdbd8..c646e2ec 100644 --- a/src/language_reading_predictors/statistical_models/pipeline.py +++ b/src/language_reading_predictors/statistical_models/pipeline.py @@ -5877,6 +5877,32 @@ 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. + # Guarded like every other secondary artefact here: the ``prior`` group is + # grafted on by the (itself guarded) ``_attach_prior_groups``, and reading it + # off a trace that has none raises AttributeError. This runs after sampling, + # so an unguarded failure would lose a completed reporting-tier fit over a + # prior check. Degrade to a warning and no CSV instead. + try: + 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, + ) + except Exception as exc: + rprint(f"[yellow]prior_pushforward skipped: {exc}[/yellow]") + else: + 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..965fdd72 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,52 @@ def _kf_csv_row(output_dir, name: str) -> dict | None: return df.iloc[0].to_dict() +# Values in psense_summary.csv's ``diagnosis`` column that mean "not flagged". +# "✓" is what arviz_stats actually writes for a clear parameter; the wording +# variants match sensitivity.tau_psense_status, which is the established +# convention for reading this column. The rest are defensive placeholders. +_PSENSE_CLEAR_MARKERS = frozenset( + { + "✓", + "-", + "nan", + "none", + "ok", + "no concern", + "no conflict", + "no prior-data conflict", + } +) + + +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. + + ``arviz_stats`` writes a **tick** for an unflagged parameter, not a blank, and that + is the single most common value in the stored suite (1117 of 2648 rows). Treating + it as a diagnosis would publish a "prior-sensitive" caution on a *clean* estimate — + so the clear markers are matched explicitly. Anything unrecognised is deliberately + treated as a flag: an unknown marker should over-warn, not go silent.""" + 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 _PSENSE_CLEAR_MARKERS: + 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 +3709,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..4f57d4b7 100644 --- a/tests/statistical_models/test_key_findings.py +++ b/tests/statistical_models/test_key_findings.py @@ -489,6 +489,66 @@ 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 + + +@pytest.mark.parametrize("clear_marker", ["✓", "-", ""]) +def test_level_factors_no_warning_when_t2_psense_clear(tmp_path, clear_marker): + # A clear diagnosis for the t2 term produces no caution bullet. + # + # "✓" is the case that matters: it is what arviz_stats actually writes for an + # unflagged parameter and the most common value in the stored suite, while "-" + # never appears in it. Treating the tick as a diagnosis published a + # "prior-sensitive" caution on six of the eleven level-factor reporting fits + # whose t2 term is in fact clear. + d = _setup_dir(tmp_path, "level_factors") + _write_csv(d, "rope_summary.csv", _rope_row()) + pd.DataFrame( + {"prior": [0.01], "likelihood": [0.02], "diagnosis": [clear_marker]}, + 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_level_factors_warns_on_an_unrecognised_psense_marker(tmp_path): + # An unknown marker is treated as a flag, not silently dropped: under-warning on + # a real prior-data conflict is the worse failure of the two. + d = _setup_dir(tmp_path, "level_factors") + _write_csv(d, "rope_summary.csv", _rope_row()) + pd.DataFrame( + {"prior": [0.9], "likelihood": [0.01], "diagnosis": ["some future verdict"]}, + index=["b_grp_time[1]"], + ).to_csv(d / "psense_summary.csv") + payload = generate_key_findings(d) + assert "warning" 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..8c0a546a 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,72 @@ 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_level_prior_pushforward_raises_without_a_prior_group(): + # ``_attach_prior_groups`` is itself guarded, so a trace can reach this point with + # no ``prior`` group. Reading one off then raises rather than returning empty — + # which is why fit_level_factors wraps the call. Pin the contract so the call + # site's guard is not later removed as redundant. + posterior = xr.Dataset( + {"eta": (("chain", "draw", "obs_id"), np.zeros((1, 4, 3)))}, + coords={"chain": [0], "draw": np.arange(4), "obs_id": np.arange(3)}, + ) + trace = SimpleNamespace(posterior=posterior) + with pytest.raises(AttributeError): + level_prior_pushforward( + trace, + phase=np.array([0, 1, 2]), + G=np.array([1.0, 0.0, 1.0]), + n_trials=79, + ) + + 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.