Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions src/language_reading_predictors/statistical_models/pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand Down
81 changes: 80 additions & 1 deletion src/language_reading_predictors/statistical_models/reporting.py
Original file line number Diff line number Diff line change
Expand Up @@ -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).

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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).

Expand Down Expand Up @@ -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
Comment on lines +3348 to +3351


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)
Expand Down Expand Up @@ -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")
)
Expand Down
40 changes: 40 additions & 0 deletions tests/statistical_models/test_key_findings.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
48 changes: 48 additions & 0 deletions tests/statistical_models/test_reporting.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
evidence_label,
favoured_direction,
joint_treatment_marginals,
level_prior_pushforward,
level_t2_marginal_effect,
longitudinal_conditional_slopes,
longitudinal_factor_correlations,
Expand Down Expand Up @@ -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.
Expand Down