diff --git a/METHODS.md b/METHODS.md index f441f0c7..86ef85d4 100644 --- a/METHODS.md +++ b/METHODS.md @@ -58,6 +58,8 @@ For the post-hoc P/N floor-transition estimands, a power-scaling conflict for `t **Evaluate.** Check convergence _before_ reading any estimate: R-hat ≤ 1.01, adequate effective sample size (ESS ≥ 400), and zero divergent transitions. A non-converged or divergent fit's posterior is not interpretable — fix the model, do not report it. The pipeline writes a `diagnostics_summary.json` whose pass/fail gate requires **0 divergences, BFMI ≥ 0.3, R-hat ≤ 1.01, ESS ≥ 400** (evaluated unrounded, over the model's free RVs plus the headline scalars) and the report renders it as a **pass/fail convergence banner first**, before any τ. This is only the sampling gate: `pareto_k.csv` separately flags importance-sampled LOO points above ArviZ's reliability threshold and makes that LOO score unreliable. In single-period ITT/joint models, where each point is one child, `scripts/influence_sensitivity.py` performs the required direct treatment-effect refit without all flagged children before a robustness claim. It standardises the full and leave-out posteriors over the same retained children, so the headline refit shift is not contaminated by changing the covariate-averaging population; the smaller composition shift and the total full-sample-to-leave-out shift are reported separately. This checks whether the treatment estimate is sensitive to those children; it does not repair or replace the unreliable LOO score. In repeated-measures random-intercept families, each point is a child-by-period row conditional on that child's fitted intercept; use an observation-level sensitivity or exact/moment-matched LOO for that same predictive target, not the whole-child ITT/joint runner. Prior- and posterior-predictive checks assess whether the working likelihood reproduces the outcome by arm, baseline band and distribution shape; a material shape flag remains a qualification even when convergence passes. The 1,000 prior draws are persisted onto `trace.nc` (the `prior` / `prior_predictive` / `log_prior` groups) so the report can show the prior-predictive check, the prior-vs-posterior overlay, the estimand-scale prior pushforward, and power-scaling sensitivity without refitting. +**Reading diagnostics — always via `sampling_quality`.** Never re-derive R-hat, ESS, BFMI or divergences from a trace by hand, in pipeline code or in a one-off script: call `statistical_models.sampling_quality.sampling_quality(trace, var_names=…)`, which returns the four signals unrounded and correctly coerced. Hand-rolled extraction has gone wrong twice for the same reason. `az.summary()` rounds to `rcParams["stats.round_to"]` (`"2g"`, two significant figures) unless passed `round_to="none"` — the **string**; `round_to=None` and `"auto"` both fall through to the rounded default, and omitting the argument returns a _string_-dtype frame that raises on float formatting. Rounding erases exactly the digits the gate turns on: across the whole gate-relevant range **every R-hat from 1.011 to 1.049 rounds to `1.0`**, so an `R-hat ≤ 1.01` test silently becomes `R-hat < 1.05` (dseinternational/research#65; found again in #440 in the exact-LOO-refit gate, which had been gating refits on rounded values since #438, and in a prototype script that reported four fits as `1.0000` when they were 1.0011–1.0022). ESS must be the minimum of `ess_bulk` and `ess_tail`, as the gate is defined, not `ess_bulk` alone. The helper deliberately does not decide pass/fail — call sites differ in which variables they gate over and in how they treat a missing BFMI — so those policies stay at the call site while the extraction stays in one place. + **Compare.** PSIS-LOO via ArviZ: prefer the higher-`elpd` model only when the difference clears its standard error (`elpd_diff` against `dse`). The interaction models are tested against their own no-interaction baselines as clean nested comparisons; `scripts/compare_statistical_models.py` collects the cross-model views. For the families with a per-child random intercept (mechanism, gain-/level-factors, DiD, dose-response, LCSM, growth) the pointwise unit is a child × phase/period row and the LOO is **conditional** (within-child, one-row-out) — elpd scores predicting a held-out row given that child's own fitted intercept, not new-child prediction — so it is comparable only across models sharing the same grouping on the same rows (which the nested comparators do). Two separate questions govern whether a nested contrast may be read, and the comparison CSVs record them in separate columns (#438). **Reliability** — `comparison_valid` — asks whether the number can be trusted at all: a comparison is refused whenever the models do not share ordered rows, or when importance sampling is unreliable (`pareto_k` above ArviZ's `good_k`). The HSGP-curve mechanism models each carry one or two such influential rows out of ~150, because the basis coefficients plus a child random intercept at this sample size make a single child-phase row pivotal for the curve near its own exposure value; the linear-mechanism pairs carry none. Those points are repaired by **exact refit** (`reloo`, at most `_RELOO_MAX_REFITS` per model), and `loo_method` / `n_exact_refits` record that the repair happened rather than leaving a reader to infer it from a silently-valid row. **Discrimination** — `elpd_verdict` — is the different question of whether a trustworthy number separates the models, and follows the standing `|elpd_diff| < 4` rule. A nested pair differing by one regularised coefficient sits near the floor of what predictive comparison can resolve at n ≈ 54, so **an inconclusive contrast is the expected outcome and is a pass, not a failed comparison**: it is recorded and reported as "LOO does not distinguish these models at this sample size", and the scientific claim rests on the coefficient's own posterior (median, 50 %/89 % intervals, direction probability, ROPE) in the usual way. Do not read the `|elpd_diff| / dse` ratio in that regime — when two models are near-identical pointwise the standard error is itself unreliable, and a ratio above 5 on a difference of ~1 elpd says the models agree, not that one wins. diff --git a/config/spellcheck/allow-en.txt b/config/spellcheck/allow-en.txt index fb3a3061..f9353518 100644 --- a/config/spellcheck/allow-en.txt +++ b/config/spellcheck/allow-en.txt @@ -84,6 +84,7 @@ Gottman Graphviz Grasman Gunn +Guttman HPDI HSGP Haldane @@ -106,6 +107,7 @@ Isager JCPP Jakulin Jarrold +Kadane Kapalková Kass Kostewicz @@ -131,7 +133,6 @@ LRPLF LRPMM LRPSURV LSAM -MAXJOBS Laan Lakens Lawson @@ -146,6 +147,7 @@ Loveall Lucas Lyster Lüdecke +MAXJOBS MCMC MCSE MIMIC @@ -219,6 +221,7 @@ RLOG RMSE RMSEA Raftery +Rasch Ratz Reichow Rhemtulla @@ -284,6 +287,7 @@ agespeak anchoring aptgram aptinfo +argmax asis assoc autoregressive @@ -308,6 +312,7 @@ brms bytree caveated celf +chokepoint cloglog coef coeffs @@ -410,8 +415,8 @@ importorskip incl invlogit isfinite -isna isinstance +isna issubset iterdir iterrows @@ -461,7 +466,6 @@ miscalibration monofont monofontoptions morphosyntax -mypy multicollinear multicollinearity multimethod @@ -469,22 +473,23 @@ multitrait mumedu mumedupost mumocc +mypy ncol ndarray networkx noconstruct nogp nohup -noqa noncentred nonlinear nonparametric nonword nonwords nonzero -nunique +noqa numchil numpy +nunique nutpie offfloor operationalisation @@ -494,6 +499,8 @@ operatorname optuna orcid othertimeread +overcover +overcoverage overdispersion overfit overfits @@ -524,11 +531,11 @@ recomputation redescribing reframing reloo +removesuffix reparameterisation reparameterise reparameterised reparameterising -removesuffix replite repr resid @@ -560,6 +567,7 @@ scikit scipy scrartcl scratchpad +secondaries sens shutil sigmas @@ -588,10 +596,13 @@ trivariate trog unadjustable unblockable -chokepoint +uncentred unclustered unconfounded unconverged +undercoverage +underdispersed +underdispersion unfittable unidimensionality univariates @@ -615,10 +626,3 @@ yarcsi ypre zproc ΔELPD -secondaries -uncentred -argmax -cloglog -disattenuated -disattenuation -surv diff --git a/notes/202607261405-binomial-exchangeability-item-difficulty-review.md b/notes/202607261405-binomial-exchangeability-item-difficulty-review.md new file mode 100644 index 00000000..cd67b1bd --- /dev/null +++ b/notes/202607261405-binomial-exchangeability-item-difficulty-review.md @@ -0,0 +1,107 @@ + + +# Review: item-difficulty non-exchangeability, the Beta-Binomial working likelihood, and the negative age association + +> [!NOTE] +> Drafted by an LLM-based AI tool (Claude Code/Fable 5). Records the 2026-07-26 review of whether ordered item difficulty (later test items are harder) undermines the Binomial/Beta-Binomial working likelihoods; the decision to keep the Beta-Binomial with targeted diagnostics and local sensitivities; and the analysis showing the negative word-reading age association is **not** a difficulty-ladder artefact. + +## The question + +Frank asked (2026-07-26): we use Binomial/Beta-Binomial likelihoods for many outcomes, but the items on these measures are likely not truly exchangeable — later words are more difficult than earlier words on word-reading and vocabulary tests. (1) Is this a concern? (2) Should we adjust the likelihood? (3) Might it explain the negative association with age for word reading — older children are more likely to be tackling more difficult words? + +## 1. What ordered item difficulty does and does not threaten + +**The sum score is the right data object regardless.** Under a Rasch-type item model — items of arbitrary, heterogeneous difficulty but equal discrimination — the total score is a sufficient statistic for the child's ability (Andersen 1977). Modelling totals therefore discards nothing; the only live question is whether the Beta-Binomial is an adequate distribution _for the total_. `METHODS.md` (Likelihood and priors) already frames it exactly this way — "a working model … not a literal claim that every test item is an exchangeable Bernoulli trial" — and mandates a likelihood sensitivity when posterior-predictive checks show material mismatch. This note is the first systematic execution of that commitment. + +**Variance: two opposing violations, and the Beta-Binomial can express only one.** Given a child's ability, a sum of unequal-difficulty items is Poisson-binomial, whose variance is _below_ the binomial with the same mean — a steep, near-Guttman ladder makes the score nearly deterministic given ability. Between-child heterogeneity beyond the linear predictor pushes the marginal variance _above_ binomial. The Beta-Binomial (and the logistic-normal-Binomial used by the within-wave levels designs) can only represent at-or-above-binomial variance: it absorbs the second violation and cannot express the first. The observable symptom of unexpressed conditional underdispersion is posterior-predictive **overcoverage** — 50 % / 90 % prediction bands covering more than 50 % / 90 % of observations. The stored `lrp-rli-itt-010` (word reading) reporting fit already shows the signature: its 50 % band covers **62.3 %** of observations (33/53). Follow-up A below pools this statistic across the whole fitted suite from the `ppc_summary.csv` every fit writes. + +**Mean shape: the true score-versus-ability curve is not one logistic.** The test characteristic curve — the sum of the item curves — is flatter than a single logistic and kinked where an item bank changes (the `EWRSWR` composite's early-word-reading → single-word-reading boundary; the documented distributional kink above ~25 words, `notes/202607241600-findings-word-reading-bands.md`). For the 170-item vocabulary scales the observed maxima are 82 (`R`) and 77 (`E`), so the logit link's ceiling curvature sits at a ceiling no child approaches: the effective ceiling is ladder steepness and discontinue rules, not `n_trials`. Stopping rules (documented for `P` in `measures.py`; the ROWPVT/EOWPVT/TROG manuals' basal/discontinue rules) additionally make the count a _curtailed_ sum — items beyond the stopping point are scored as failures without being attempted — not `n` attempted trials. + +**Interpretation: items are not equal-interval ability units.** Three items at the hard end of a ladder represent more latent progress than three at the easy end. Items-scale effects _within_ a measure remain well defined (both arms take the same instrument, and the reported marginal effect averages over the analysis rows), but items-scale magnitudes are not comparable _across_ measures, and this should remain a stated reporting caveat. + +**Why τ is largely protected.** Randomisation, the same instrument in both arms, and the row-averaged marginal effect protect the sign and approximate size of the treatment contrast; likelihood misspecification mainly distorts interval calibration and the fine detail of the items-scale conversion. The floor rule for `P`/`N` (`bernoulli_offfloor`) is already a local likelihood repair at exactly the point where the graded working model is worst. + +## 2. Decision: keep the Beta-Binomial; diagnose, then patch locally + +A wholesale move to ordinal or item-response likelihoods across the suite would cost a great deal and mostly relabel conclusions. The adopted policy is a ladder: + +1. **Diagnose from artefacts we already produce.** Pool the `ppc_summary.csv` 50 % / 90 % coverage across reporting fits by outcome (follow-up A; `scripts/ppc_coverage_sweep.py`). Systematic overcoverage on the laddered tests is the conditional-underdispersion signature; undercoverage is ordinary shape misfit. Small-denominator caveat: central intervals of a discrete count distribution overcover mechanically (a "50 %" interval on a 10-item score covers at least 50 %), so measures of similar length are compared with each other, not against an absolute 0.50. +2. **One principled link fix regardless of diagnostics: blending.** `B` is ten three-alternative forced-choice items, so the expected score cannot fall below chance (≈ 3.3 of 10), yet the standard logit link happily models sub-chance means. A guessing-floor link, mean = 1/3 + (2/3) · logit⁻¹(η), is mechanically justified (follow-up B prototypes it). +3. **If conditional underdispersion is material: the Conway–Maxwell-binomial.** The CMB (Kadane 2016) is the natural drop-in that can express _under_- as well as over-dispersion relative to the Binomial (its extra parameter ν: ν = 1 Binomial, ν > 1 underdispersed, ν < 1 overdispersed), and its normalising constant is a finite sum over 0…n, so a custom PyMC log-probability is cheap even at n = 170 (follow-up C prototypes it for `W`). +4. **If mean-shape misfit is material: a flexible-monotone own-baseline term.** The HSGP machinery and the mechanism-family knee tests are the in-house precedent. Section 3 below suggests low urgency for the age question specifically. +5. **Ordinal (continuation-ratio — the natural generative model for a ladder with discontinue rules) or full item-response modelling only as a last resort.** Item-level responses are absent from `rli_data_long.csv` (block totals only). Settled same day: the original scoring sheets are not readily available, closing this route — see §7. + +## 3. The negative word-reading age association is not a ladder artefact + +The suite-wide finding under test: conditional on own baseline, linear age is credibly negative for word reading — **−0.12** on the logit scale in the stacked gain-factors models (P(< 0) = 0.99; `notes/202607082140-statistical-models-full-reporting-fit.md`, `notes/202607111100-replite-full-statistical-fit.md`), pointing the same way in the single-window ITT (`lrp-rli-itt-010` stored reporting trace: `gamma_A` median **−0.08**, 89 % [−0.21, +0.05], P(< 0) = 0.83 — one transition, so wider), **−0.26** in the adjusted-association model `adj-065`, and **−0.34** in the errors-in-variables mechanism fit (`notes/202607141700-lrp228-eiv-mechanism.md`). The hypothesis: older children sit higher on the item ladder, where words are harder, so the association is a measurement artefact. + +**Stating the mechanism precisely already narrows it.** Two children with the same baseline score face the _same next words_ regardless of age, so item difficulty cannot reach the age coefficient directly. It can only leak in through **unmodelled baseline-curvature absorbed by age** via the age–baseline correlation (observed: corr(age, W at t1) = **+0.37** in the randomised-window sample, n = 53). + +**Real-data check (t1 → t2, n = 53).** Regressing W post-score on baseline + age + arm, then adding hinge terms at 25 and 30 items (the kink region), the age slope barely moves: **−1.06 → −1.02 items per SD of age** (SE ≈ 0.55) on the raw scale; **−0.066 → −0.028** (SE ≈ 0.14) on the empirical-logit scale; the hinge coefficients themselves are null (±0.4, SE ≈ 1.0). Only three children start above 25 items, so power above the kink is minimal — but the gradient-boosting step, which fits the baseline fully flexibly with trees, reaches the same negative age ranking for word-reading gain (`notes/202607021100-word-reading-reporting-tier-ranking.md`), which is the same test by another route. + +**Simulation under the null.** `scripts/age_artefact_check.py` generates Rasch-type data from three 79-item ladders — kinked (30 easy early-word-reading items, then a single-word-reading difficulty jump), smoothly graded, and homogeneous — calibrated to the observed baseline mean/SD, corr(age, ability) = 0.37, mean gain ≈ 3.4 items, and **zero true age effect on growth**, then fits the ITT-style regression (1,000 replicates, n = 53, seed 20260726): + +| ladder | age coefficient (empirical-logit scale), mean [2.5, 97.5 %] | P(coef ≤ −0.12) | +| ----------------- | ----------------------------------------------------------- | --------------- | +| kinked (EWR/SWR) | **+0.175** [−0.090, +0.459] | 0.02 | +| smoothly graded | **+0.184** [−0.073, +0.452] | 0.01 | +| homogeneous items | **+0.218** [−0.074, +0.524] | 0.01 | + +Every ladder produces a **positive** mean age coefficient under the null (≈ +0.4 to +0.5 items per SD of age on the raw scale). The dominant channel runs the other way from the hypothesis: with a noisy baseline, age is a proxy for true ability, and higher true ability means a higher post-score — an errors-in-baseline channel that biases the age coefficient _upwards_, which the difficulty ladder does nothing to reverse. The chance that this data-generating process yields the observed −0.12 is 1–2 %. + +**The suite corroborates the direction.** If the errors-in-baseline channel biases the naive age slope upwards, then models that correct or bypass it should report _more_ negative slopes — and they do, in exactly that order: −0.08 (single-window ITT) / −0.12 (stacked gain-factors) → −0.26 (`adj-065`) → −0.34 (the errors-in-variables fit). The negative age association also appears on the taught-word set (`TE`, a bespoke item set with no difficulty ladder: −0.10 in the quick check) and on grammar, not only the laddered reading tests. + +**Conclusion.** The negative age association for word reading is not a likelihood or item-difficulty artefact; if anything the ladder-plus-measurement-error structure _masks_ part of a more negative true association. What this analysis cannot do is separate developmental timing from trajectory selection — being older at the same score _means_ having grown more slowly historically, and no cross-sectional adjustment can distinguish "older children gain less" from "slower-trajectory children are older when they arrive at any given score". The existing reading in the notes ("older children in this cohort having already banked their fastest-moving years", `notes/202607082140-statistical-models-full-reporting-fit.md`) stands, now with the artefact explanation ruled out; it is also consistent with the RLI trial literature (younger children made more progress; `notes/202606281600-literature-review.md`). + +**Caveats.** The quick regressions use empirical-logit ordinary least squares with naive standard errors on a single randomised window; the simulation's ability distribution is moment-matched on a coarse grid. These are direction-and-magnitude checks, not replacement estimates; the registered models remain the authoritative quantities. + +## 4. Follow-up A — posterior-predictive coverage sweep (results) + +`scripts/ppc_coverage_sweep.py` pooled 380 `ppc_summary.csv` coverage rows from **190 stored reporting fits** (2026-07-26; three output directories lacked the files). Headline: 50 % prediction bands **overcover everywhere** — pooled coverage by outcome runs **0.65 (`T`) to 0.92** (the historical-cohort measures), with `W` at **0.76** across 66 models, `L` 0.71, `R` 0.68, `E` 0.75, `B` 0.74, `TE` 0.80, `N` 0.82. By family the cross-sectional ITT suite is mildest (**0.66**; `aligned` 0.63, `corr_factor` 0.58 closest to nominal) and the longitudinal random-intercept families largest (`did` 0.82, `lcsm` 0.85, `historical_growth` 0.85). The 90 % bands sit much closer to nominal (0.95–1.00). The off-floor mode has 2–31 group cells per level and is uninformative about dispersion. + +Interpretation requires the two mechanical inflators named in §2: equal-tailed intervals on a _discrete_ count overcover by construction (worst for short scales and for the floor-concentrated score distributions, where the modal low counts sit inside any central band), and the longitudinal families' in-sample checks condition on fitted child effects. The pointed question — is the remaining overcoverage evidence that the Beta-Binomial's variance floor binds (conditional underdispersion)? — is answered by follow-up C: **no**. A dispersion family free to go below binomial variance chooses overdispersion for `W`, and its coverage barely improves (0.62 → 0.60). The 50 %-band overcoverage is therefore substantially mechanical, not a likelihood defect; the sweep stays useful as a _relative_ monitor (a measure drifting away from its peers) after future full runs. + +## 5. Follow-up B — blending guessing-floor link prototype (results) + +`scripts/likelihood_sensitivity_prototypes.py blending-chance-floor` refits the `lrp-rli-itt-008` structure (n = 54, 10 items) twice with shared priors and sampler (4 chains × 2,000 draws, `nutpie`, `target_accept` 0.92, seed 20260726; both fits: R-hat ≤ 1.002, min ESS ≥ 4,137, 0 divergences — comfortably inside the suite gate), differing only in the link. Validation: the standard arm reproduces the stored reporting fit (τ **+0.44** here vs **+0.45** stored; items-scale effect **+0.96** vs **+0.99**). + +| model | τ (logit), median [89 %] | treatment effect (items), median [89 %] | P(> 0) | +| ------------------------------------------- | ------------------------ | --------------------------------------- | ------ | +| standard Beta-Binomial | +0.44 [+0.09, +0.81] | **+0.96** [+0.20, +1.78] | 0.98 | +| chance-floor (mean = 1/3 + 2/3 · logit⁻¹ η) | +0.40 [−0.14, +0.92] | **+0.49** [−0.16, +1.10] | 0.90 | + +Respecting the guessing floor roughly **halves the items-scale treatment effect and weakens the evidence from strong to moderate** — under the chance-floor link, movement in the sub-chance region is attributed to guessing noise rather than ability, and what remains is the contrast on the above-chance share. Against the adopted ROPE δ(B) = 1 item, the standard headline sits at δ while the chance-floor estimate sits at about δ/2. PSIS-LOO cannot separate the two links (elpd difference 1.0 with difference-SE 2.7 — a tie), so the data do not adjudicate; the argument for the chance-floor link is mechanistic (the instrument is three-alternative forced choice), not predictive. Side-observations: `gamma_own` rises to ≈ 1.0 under the chance-floor link (closer to the autoregressive prior's centre), and predictive coverage is identical between links (0.63 / 0.93), so this is purely an estimand-scale sensitivity, not a fit improvement. + +## 6. Follow-up C — word-reading Conway–Maxwell-binomial prototype (results) + +`scripts/likelihood_sensitivity_prototypes.py word-reading-cmb` refits the `lrp-rli-itt-010` structure (n = 53, 79 items; same shared priors and sampler; both fits: R-hat ≤ 1.002, min ESS ≥ 1,494, 0 divergences), swapping the Beta-Binomial for the Conway–Maxwell-binomial with ν ~ LogNormal(0, 0.5). Validation: the standard arm reproduces the stored reporting fit to three decimals (τ +0.354 both). + +| model | τ (logit), median [89 %] | treatment effect (items), median [89 %] | P(> 0) | dispersion | 50 % / 90 % coverage | +| ----------------------- | ------------------------ | --------------------------------------- | ------ | ------------------------- | -------------------- | +| standard Beta-Binomial | +0.35 [+0.12, +0.60] | **+2.38** [+0.77, +4.01] | 0.99 | κ = 48 [30, 76] | 0.62 / 0.94 | +| Conway–Maxwell-binomial | +0.20 [+0.07, +0.34] | **+2.49** [+0.96, +4.02] | 0.99 | **ν = 0.50 [0.37, 0.67]** | 0.60 / 0.89 | + +Two findings. First, ν is decisively **below 1**: given the freedom to express sub-binomial variance, the word-reading residuals demand *over*dispersion — the Beta-Binomial's direction — so there is no evidence its variance floor binds in the cross-sectional models. (Conditional underdispersion, if it exists, lives _within_ child; the stacked gain-factors family with its child intercepts would be the place to probe it, not these one-row-per-child fits.) Second, the **items-scale treatment effect is invariant** (+2.38 vs +2.49, P(> 0) = 0.99 both) even though the logit-scale τ halves — latent-scale coefficients are not comparable across dispersion families, which is exactly why the house reporting convention leads with the items-scale marginal effect. The headline `W` result is likelihood-robust. + +## 7. Would a Rasch model help? (assessed and closed) + +Asked directly by Frank after the follow-ups. Assessment: a Rasch (one-parameter item-response) layer is the only route to the things sum scores cannot identify — an interval-scale ability metric (resolving the items-are-not-equal-interval caveat), an empirically pinned test characteristic curve, per-wave measurement standard errors to systematise the errors-in-variables corrections, and differential-item-functioning checks (notably whether any standardised-vocabulary items overlap taught content, an arm-DIF validity question for the transfer claims). It would **not** change the headline results: the sum score is sufficient for ability under the Rasch model itself (Andersen 1977), the items-scale treatment effect proved invariant to a materially different count family (§6), and a Rasch treatment of the three-alternative forced-choice blending test reduces to the same guessing-floor idea as §5. Even with item data, calibration would be thin exactly where the ladder matters — 54 children × 4 waves, with the discontinue rules leaving the hard end of each test rarely or never attempted — and publisher difficulties from typically-developing norming samples are a poor substitute (population DIF is likely). + +**Closed 2026-07-26: Frank confirmed the original scoring sheets are not readily available**, so no item-level response matrix can be assembled without a substantial archival effort that the expected gains do not justify. The end state for this reanalysis is therefore the Beta-Binomial working likelihood, the items-scale reporting caveat, and the blending chance-floor sensitivity. If item-level records ever do surface, the value-ordered agenda recorded here is: (1) a one-off Rasch calibration of `W` to pin the test characteristic curve and re-express the ≈ 2.4-item effect in ability logits; (2) Rasch-based measurement error feeding the mechanism/mediation couplings; (3) arm-DIF screening on `R`/`E` — and nothing suite-wide. + +## Decisions + +Recorded 2026-07-26; owners in brackets. + +1. **Keep the Beta-Binomial as the default working likelihood** for bounded-count outcomes. The CMB probe found overdispersion (ν = 0.50), vindicating the family's direction; the items-scale treatment estimand is invariant to the swap. No suite-wide likelihood change. +2. **Adopt `scripts/ppc_coverage_sweep.py` as a standing diagnostic** after full reporting sweeps, read _relatively_ (a measure drifting from its peers), with the discreteness and in-sample-random-effect caveats stated in §4. [rerun with future sweeps] +3. **Blending (`B`) headline is link-sensitive and needs a stated sensitivity.** The ≈ 1-item, strong-evidence ITT effect halves to ≈ 0.5 items, moderate evidence, under the mechanically-motivated guessing-floor link, and LOO cannot separate the links. Whether to promote the chance-floor variant to a registered model (a `link=` option in `factories.py` plus guard tests) — or at minimum to caveat `B`'s items-scale effect and its ROPE reading in reports — is a substantive call. [Frank / education lead] +4. **No Conway–Maxwell-binomial adoption; the ordinal / item-response route is closed.** Settled 2026-07-26: Frank confirmed the original scoring sheets are not readily available, so no item-level data can be assembled for a Rasch or ordinal measurement model (§7). Revisit only if item-level records surface, per the value-ordered agenda in §7. [closed] +5. **Propose one reporting-guardrail sentence for `METHODS.md`**: items-scale effects are averages over the realised score distribution; items are not equal-interval ability units and are not comparable across measures. Not made in this change. [future docs change] +6. **The negative word-reading age association stands as substantive** (§3): not a difficulty-ladder or likelihood artefact; developmental-timing versus trajectory-selection remains unresolvable cross-sectionally. Attribution corrected: −0.12 (P = 0.99) is the stacked gain-factors figure; the single-window ITT posterior is −0.08 (P(< 0) = 0.83). + +## References + +- Andersen, E. B. (1977). Sufficient statistics and latent trait models. _Psychometrika_, 42(1), 69–81. DOI [10.1007/BF02293746](https://doi.org/10.1007/BF02293746) +- Kadane, J. B. (2016). Sums of possibly associated Bernoulli variables: the Conway–Maxwell-binomial distribution. _Bayesian Analysis_, 11(2), 403–420. DOI [10.1214/15-BA955](https://doi.org/10.1214/15-BA955) +- Burgoyne, K., Duff, F. J., Clarke, P. J., Buckley, S., Snowling, M. J., & Hulme, C. (2012). Efficacy of a reading and language intervention for children with Down syndrome: a randomized controlled trial. _Journal of Child Psychology and Psychiatry_, 53(10), 1044–1053. DOI [10.1111/j.1469-7610.2012.02557.x](https://doi.org/10.1111/j.1469-7610.2012.02557.x) diff --git a/scripts/age_artefact_check.py b/scripts/age_artefact_check.py new file mode 100644 index 00000000..942a6355 --- /dev/null +++ b/scripts/age_artefact_check.py @@ -0,0 +1,193 @@ +# Copyright (c) 2026 Down Syndrome Education International and contributors +# SPDX-License-Identifier: AGPL-3.0-or-later + +"""Can item-difficulty non-exchangeability explain the negative word-reading age slope? + +Companion analysis to +notes/202607261405-binomial-exchangeability-item-difficulty-review.md. The suite +finds a credibly negative linear-age precision term for word reading conditional on +own baseline (``gamma_A`` ~ -0.12 logit; adj-065 -0.26; the errors-in-variables fit +-0.34). The hypothesis under test: older children sit higher on the item ladder where +words are harder, so the slope is a measurement artefact. Two children with the same +baseline score face the same next words regardless of age, so the artefact can only +enter through unmodelled baseline curvature absorbed by age via the age-baseline +correlation. This script tests that channel twice: + +Part 1 (real data, randomised t1->t2 window): regress W post on baseline + age + arm +with the baseline entered (a) linearly and (b) with hinge terms at 25/30 items (the +documented kink region). If the negative age slope is curvature leak, it should +attenuate under (b). + +Part 2 (simulation under the null): Rasch-type 79-item ladders (kinked EWR/SWR-style, +smoothly graded, homogeneous), age correlated with baseline ability at the observed +level, and growth INDEPENDENT of age (true age effect = 0). The fitted age +coefficient measures how much spurious slope each ladder can generate. + +Result recorded in the note: every ladder yields a POSITIVE mean age coefficient +under the null (errors-in-baseline makes age a proxy for true ability), so the +observed negative slope cannot be a difficulty-ladder artefact. + +Usage: + python scripts/age_artefact_check.py [--reps 1000] [--seed 20260726] [--data data/rli_data_long.csv] +""" + +from __future__ import annotations + +import argparse + +import numpy as np +import pandas as pd + +N_ITEMS = 79 + + +def z(x): + x = np.asarray(x, float) + return (x - x.mean()) / x.std() + + +def elogit(k, n): + k = np.asarray(k, float) + return np.log((k + 0.5) / (n - k + 0.5)) + + +def ols(y, cols): + names = list(cols) + X = np.column_stack([np.ones(len(y))] + [np.asarray(cols[c], float) for c in names]) + b, *_ = np.linalg.lstsq(X, y, rcond=None) + resid = y - X @ b + dof = len(y) - X.shape[1] + s2 = resid @ resid / dof + cov = s2 * np.linalg.inv(X.T @ X) + se = np.sqrt(np.diag(cov)) + return {n_: (bi, si) for n_, bi, si in zip(["const"] + names, b, se, strict=True)} + + +def show(res, keys): + for k in keys: + b, s = res[k] + print(f" {k:>16}: {b:+7.3f} (se {s:.3f})") + + +def real_data_checks(data_path: str) -> tuple[int, float, float, float, float]: + df = pd.read_csv(data_path) + w = df.pivot_table(index="subject_id", columns="time", values="ewrswr") + a = df.pivot_table(index="subject_id", columns="time", values="age") + g = df.groupby("subject_id")["group"].first() + + d = pd.DataFrame({"W1": w[1], "W2": w[2], "age1": a[1], "G": 2 - g}).dropna() + n = len(d) + print(f"== Real data: W (ewrswr), t1->t2, n = {n}") + print(f" W1 mean {d.W1.mean():.2f} sd {d.W1.std():.2f}; W2 mean {d.W2.mean():.2f}; " + f"mean gain {(d.W2 - d.W1).mean():.2f}") + r_age_w1 = float(np.corrcoef(d.age1, d.W1)[0, 1]) + print(f" corr(age1, W1) = {r_age_w1:+.3f} corr(age1, W2-W1) = " + f"{np.corrcoef(d.age1, d.W2 - d.W1)[0, 1]:+.3f}") + + bands = pd.cut(d.W1, [-0.5, 5, 15, 25, 80], labels=["0-5", "6-15", "16-25", ">25"]) + tab = d.assign(gain=d.W2 - d.W1, band=bands).groupby("band", observed=True).agg( + n=("gain", "size"), mean_gain=("gain", "mean"), mean_age=("age1", "mean")) + print(" gain and age by baseline band:") + print(" " + tab.to_string().replace("\n", "\n ")) + + y_raw = np.asarray(d.W2, float) + y_el = elogit(d.W2, N_ITEMS) + h25 = np.clip(d.W1 - 25, 0, None) + h30 = np.clip(d.W1 - 30, 0, None) + + print("\n [raw scale] W2 ~ z(W1) + z(age) + G") + show(ols(y_raw, {"zW1": z(d.W1), "zage": z(d.age1), "G": d.G}), ["zW1", "zage", "G"]) + print(" [raw scale] + hinge(W1-25) + hinge(W1-30)") + show(ols(y_raw, {"zW1": z(d.W1), "h25": h25, "h30": h30, "zage": z(d.age1), "G": d.G}), + ["zW1", "h25", "h30", "zage", "G"]) + print(" [empirical logit] elogit(W2) ~ z(W1) + z(age) + G (~ITT structure)") + show(ols(y_el, {"zW1": z(d.W1), "zage": z(d.age1), "G": d.G}), ["zW1", "zage", "G"]) + print(" [empirical logit] + hinge terms") + show(ols(y_el, {"zW1": z(d.W1), "h25": h25, "h30": h30, "zage": z(d.age1), "G": d.G}), + ["zW1", "h25", "h30", "zage", "G"]) + + for col, nmax, label in [("rowpvt", 170, "R (ROWPVT, graded ladder)"), + ("b1extau", 24, "TE (taught words, no ladder)")]: + p = df.pivot_table(index="subject_id", columns="time", values=col) + dd = pd.DataFrame({"y1": p[1], "y2": p[2], "age1": a[1], "G": 2 - g}).dropna() + r = ols(elogit(dd.y2, nmax), {"zpre": z(dd.y1), "zage": z(dd.age1), "G": dd.G}) + b_, s_ = r["zage"] + print(f" {label}: n={len(dd)} empirical-logit age coef {b_:+.3f} (se {s_:.3f})") + + return n, float(d.W1.mean()), float(d.W1.std()), r_age_w1, float((d.W2 - d.W1).mean()) + + +def simulate(n, w1_mean, w1_sd, r_age_w1, gain_obs, reps, rng): + print("\n== Simulation: ladders with ZERO true age effect on growth " + f"(n = {n}, reps = {reps})") + + ladders = { + "kinked (EWR easy, SWR jump)": np.r_[np.linspace(-2.5, 0.5, 30), np.linspace(1.5, 5.0, 49)], + "smooth graded": np.linspace(-2.5, 5.0, 79), + "homogeneous items": np.full(N_ITEMS, 2.0), + } + + def score(theta, b): + p = 1 / (1 + np.exp(-(theta[:, None] - b[None, :]))) + return (rng.random(p.shape) < p).sum(axis=1) + + def calibrate(b): + best = None + for mu in np.linspace(-4, 2, 25): + for sd in np.linspace(0.5, 2.5, 9): + th = rng.normal(mu, sd, 4000) + k = score(th, b) + loss = (k.mean() - w1_mean) ** 2 + (k.std() - w1_sd) ** 2 + if best is None or loss < best[0]: + best = (loss, mu, sd) + return best[1], best[2] + + for name, b_items in ladders.items(): + mu_t, sd_t = calibrate(b_items) + d0 = None + for cand in np.linspace(0.05, 1.5, 30): + th = rng.normal(mu_t, sd_t, 4000) + gg = score(th + cand, b_items).mean() - score(th, b_items).mean() + if d0 is None or abs(gg - gain_obs) < d0[0]: + d0 = (abs(gg - gain_obs), cand) + delta0 = d0[1] + + coefs_lin, coefs_hinge, coefs_raw = [], [], [] + for _ in range(reps): + th1 = rng.normal(mu_t, sd_t, n) + th1z = (th1 - mu_t) / sd_t + age = r_age_w1 * th1z + np.sqrt(1 - r_age_w1**2) * rng.normal(size=n) + delta = rng.normal(delta0, 0.3, n) # independent of age: true effect 0 + k1 = score(th1, b_items) + k2 = score(th1 + delta, b_items) + yl = elogit(k2, N_ITEMS) + try: + coefs_lin.append(ols(yl, {"zW1": z(k1), "zage": z(age)})["zage"][0]) + hh = np.clip(k1 - 25, 0, None) + coefs_hinge.append(ols(yl, {"zW1": z(k1), "h25": hh, "zage": z(age)})["zage"][0]) + coefs_raw.append(ols(k2.astype(float), {"zW1": z(k1), "zage": z(age)})["zage"][0]) + except np.linalg.LinAlgError: + continue + cl, ch, cr = map(np.asarray, (coefs_lin, coefs_hinge, coefs_raw)) + print(f" {name}: theta~N({mu_t:.2f},{sd_t:.2f}), delta0={delta0:.2f}") + print(f" empirical-logit age coef, linear baseline : mean {cl.mean():+.4f} " + f"[2.5,97.5%: {np.percentile(cl, 2.5):+.3f}, {np.percentile(cl, 97.5):+.3f}] " + f"P(<=-0.12) = {(cl <= -0.12).mean():.2f}") + print(f" empirical-logit age coef, hinge baseline : mean {ch.mean():+.4f}") + print(f" raw-scale age coef : mean {cr.mean():+.4f} items") + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("--reps", type=int, default=1000) + ap.add_argument("--seed", type=int, default=20260726) + ap.add_argument("--data", default="data/rli_data_long.csv") + args = ap.parse_args() + + rng = np.random.default_rng(args.seed) + n, w1_mean, w1_sd, r_age_w1, gain_obs = real_data_checks(args.data) + simulate(n, w1_mean, w1_sd, r_age_w1, gain_obs, args.reps, rng) + + +if __name__ == "__main__": + main() diff --git a/scripts/design_analysis.py b/scripts/design_analysis.py index 64c48ce8..604238d2 100644 --- a/scripts/design_analysis.py +++ b/scripts/design_analysis.py @@ -34,7 +34,6 @@ import shutil import warnings -import arviz as az import matplotlib import matplotlib.pyplot as plt import numpy as np @@ -47,6 +46,9 @@ from language_reading_predictors.statistical_models.itt import resolve_itt_run_plan from language_reading_predictors.statistical_models.measures import MEASURES, ROPE_DELTA from language_reading_predictors.statistical_models.preprocessing import load_and_prepare +from language_reading_predictors.statistical_models.sampling_quality import ( + sampling_quality, +) matplotlib.use("Agg") warnings.filterwarnings("ignore") @@ -95,19 +97,13 @@ def _refit_convergence(model, idata) -> dict: would fall through to ``rcParams["stats.round_to"]`` at 2 sig figs and let a borderline R-hat pass); free RVs match the headline gate's coverage. """ + n_div = None try: free = [rv.name for rv in model.free_RVs] - summ = az.summary(idata, var_names=free, round_to="none", kind="diagnostics") - max_rhat = float(np.nanmax(summ["r_hat"].values)) - min_ess = float(np.nanmin(summ[["ess_bulk", "ess_tail"]].min(axis=1).values)) + signals = sampling_quality(idata, var_names=free) + max_rhat, min_ess, n_div = signals.max_rhat, signals.min_ess, signals.n_divergences except Exception: # pragma: no cover - defensive max_rhat, min_ess = float("nan"), float("nan") - n_div = None - try: - if "diverging" in idata.sample_stats: - n_div = int(np.asarray(idata.sample_stats["diverging"].values).sum()) - except Exception: # pragma: no cover - defensive - pass converged = bool( np.isfinite(max_rhat) and max_rhat <= 1.01 diff --git a/scripts/likelihood_sensitivity_prototypes.py b/scripts/likelihood_sensitivity_prototypes.py new file mode 100644 index 00000000..a543cd8c --- /dev/null +++ b/scripts/likelihood_sensitivity_prototypes.py @@ -0,0 +1,245 @@ +# Copyright (c) 2026 Down Syndrome Education International and contributors +# SPDX-License-Identifier: AGPL-3.0-or-later + +"""Prototype likelihood sensitivities for the bounded-count working model. + +Two one-off sensitivity fits accompanying +notes/202607261405-binomial-exchangeability-item-difficulty-review.md. Both mirror +the uniform DAG-faithful single-outcome ITT structure (``factories.build_itt_model``: +``eta = alpha + gamma_own * pre_logit + gamma_A * A_std + tau * G``), use the shared +priors from ``priors.py``, load rows with +``preprocessing.load_and_prepare(phase_mode="itt", outcomes=(symbol,))``, and differ +ONLY in the observation model — so any difference in the treatment summaries is +attributable to the likelihood alone. + +``blending-chance-floor`` + ITT B (``lrp-rli-itt-008`` structure). Standard Beta-Binomial versus a + guessing-floor link ``mu = 1/3 + (2/3) * sigmoid(eta)``: blending is + three-alternative forced choice, so the expected score cannot fall below chance + (about 3.3 of 10), yet the standard logit link admits sub-chance means. The + shared priors are kept identical in both arms so the comparison isolates the + link; note ``gamma_own``'s autoregressive centring at 1 was calibrated for the + identity relation between pre- and post-logit and is retained unchanged. + +``word-reading-cmb`` + ITT W (``lrp-rli-itt-010`` structure). Beta-Binomial versus Conway-Maxwell- + binomial (Kadane 2016, DOI 10.1214/15-BA955), which can express UNDER- as well + as over-dispersion relative to the Binomial (``nu`` > 1 underdispersed, = 1 + Binomial, < 1 overdispersed). Heterogeneous item difficulty makes the score + conditionally underdispersed (Poisson-binomial), which the Beta-Binomial cannot + represent — its variance floor is the Binomial. + +These are NOT registered models: no diagnostics gate, no report artefacts, and +lighter sampling than the reporting preset. Promotion path if adopted: a +``likelihood=``/link option in ``factories.py`` with guard tests, a registered +variant spec, and the usual pipeline artefacts. + +Run from the repo root (worktree checkouts need ``PYTHONPATH=src``): + python scripts/likelihood_sensitivity_prototypes.py blending-chance-floor + python scripts/likelihood_sensitivity_prototypes.py word-reading-cmb +""" + +from __future__ import annotations + +import argparse + +import arviz as az +import numpy as np +import pymc as pm +import pytensor.tensor as pt +from scipy.special import gammaln +from scipy.stats import betabinom as sp_betabinom + +from language_reading_predictors.statistical_models import priors as _priors +from language_reading_predictors.statistical_models.preprocessing import load_and_prepare +from language_reading_predictors.statistical_models.sampling_quality import ( + sampling_quality, +) + +EPSILON = 1e-6 +SEED = 20260726 +Q50 = (0.25, 0.75) +Q89 = (0.055, 0.945) + + +def _summary(x: np.ndarray, label: str) -> str: + med = np.median(x) + lo50, hi50 = np.quantile(x, Q50) + lo89, hi89 = np.quantile(x, Q89) + return (f"{label}: median {med:+.3f} 50% [{lo50:+.3f}, {hi50:+.3f}] " + f"89% [{lo89:+.3f}, {hi89:+.3f}] P(>0) = {(x > 0).mean():.3f}") + + +def _diagnostics(idata) -> str: + return sampling_quality(idata).summary_line() + + +def _fit(symbol: str, observation: str, draws: int, tune: int, chains: int, + target_accept: float, sampler: str): + """Fit one ITT-structure variant; ``observation`` in {bb, bb_chance_floor, cmb}.""" + prepared = load_and_prepare(phase_mode="itt", outcomes=(symbol,)) + post = prepared.post_counts[symbol] + keep = ~np.isnan(post) + pre = prepared.pre_logit[symbol][keep] + a_std = prepared.A_std[keep] + g = prepared.G.astype(float)[keep] + y = post[keep].astype(np.int64) + n_trials = prepared.n_trials[symbol] + n_obs = int(keep.sum()) + + with pm.Model(coords={"obs_id": np.arange(n_obs)}) as model: + alpha = _priors.alpha_prior().to_pymc("alpha") + tau = _priors.tau_prior().to_pymc("tau") + gamma_own = _priors.gamma_own_prior().to_pymc("gamma_own") + gamma_A = _priors.gamma_age_prior().to_pymc("gamma_A") + eta = alpha + gamma_own * pre + gamma_A * a_std + tau * g + + if observation in ("bb", "bb_chance_floor"): + kappa = _priors.kappa_prior().to_pymc("kappa") + sig = pm.math.sigmoid(eta) + mu_raw = 1.0 / 3.0 + (2.0 / 3.0) * sig if observation == "bb_chance_floor" else sig + mu = pm.math.clip(mu_raw, EPSILON, 1 - EPSILON) + pm.BetaBinomial("y_post", n=n_trials, alpha=mu * kappa, beta=(1 - mu) * kappa, + observed=y, dims="obs_id") + else: # cmb + nu = pm.LogNormal("nu", mu=0.0, sigma=0.5) + j = np.arange(n_trials + 1) + log_c = gammaln(n_trials + 1) - gammaln(j + 1) - gammaln(n_trials + 1 - j) + terms = nu * pt.constant(log_c)[None, :] + eta[:, None] * pt.constant(j.astype(float))[None, :] + m = pt.max(terms, axis=1, keepdims=True) + log_z = (m + pt.log(pt.sum(pt.exp(terms - m), axis=1, keepdims=True)))[:, 0] + loglik = nu * pt.constant(log_c[y]) + y * eta - log_z + pm.Potential("y_cmb", loglik.sum()) + + kwargs = dict(draws=draws, tune=tune, chains=chains, target_accept=target_accept, + random_seed=SEED, progressbar=False) + if sampler == "nutpie": + try: + idata = pm.sample(nuts_sampler="nutpie", **kwargs) + except Exception as exc: # pragma: no cover - environment-dependent + print(f" [nutpie unavailable for {observation} ({exc!r}); falling back to pymc]") + idata = pm.sample(**kwargs) + else: + idata = pm.sample(**kwargs) + if observation != "cmb": + pm.compute_log_likelihood(idata) + + data = {"pre": pre, "a_std": a_std, "g": g, "y": y, "n_trials": n_trials, "model": model} + return data, idata + + +def _draws(idata, names, thin_to=1500): + post = idata.posterior + flat = {k: post[k].values.reshape(-1) for k in names} + n = len(next(iter(flat.values()))) + idx = np.linspace(0, n - 1, min(thin_to, n)).astype(int) + return {k: v[idx] for k, v in flat.items()} + + +def _eta_by_arm(dr, data): + base = (dr["alpha"][:, None] + dr["gamma_own"][:, None] * data["pre"][None, :] + + dr["gamma_A"][:, None] * data["a_std"][None, :]) + return base + dr["tau"][:, None], base # (eta with G=1, eta with G=0) + + +def _cmb_pmf(eta, nu, n_trials): + """Predictive pmf per (draw, obs) over 0..n; eta (D,N), nu (D,).""" + j = np.arange(n_trials + 1) + log_c = gammaln(n_trials + 1) - gammaln(j + 1) - gammaln(n_trials + 1 - j) + logw = nu[:, None, None] * log_c[None, None, :] + eta[:, :, None] * j[None, None, :] + logw -= logw.max(axis=2, keepdims=True) + w = np.exp(logw) + return w / w.sum(axis=2, keepdims=True) + + +def _coverage_from_pmf(pmf_rows, y, level): + """Equal-tailed predictive-interval coverage from a per-row pmf over 0..n.""" + cdf = np.cumsum(pmf_rows, axis=1) + lo_q, hi_q = (1 - level) / 2, 1 - (1 - level) / 2 + lo = (cdf < lo_q).sum(axis=1) + hi = (cdf < hi_q).sum(axis=1) + return float(((y >= lo) & (y <= hi)).mean()) + + +def _ppc_coverage(observation, dr, data, thin_to=400): + idx = np.linspace(0, len(dr["tau"]) - 1, min(thin_to, len(dr["tau"]))).astype(int) + sub = {k: v[idx] for k, v in dr.items()} + eta1, eta0 = _eta_by_arm(sub, data) + eta = np.where(data["g"][None, :] == 1.0, eta1, eta0) + n_trials = data["n_trials"] + if observation == "cmb": + pmf = _cmb_pmf(eta, sub["nu"], n_trials).mean(axis=0) + else: + sig = 1 / (1 + np.exp(-eta)) + mu = 1 / 3 + (2 / 3) * sig if observation == "bb_chance_floor" else sig + kk = np.arange(n_trials + 1)[None, None, :] + a = (mu * sub["kappa"][:, None])[:, :, None] + b = ((1 - mu) * sub["kappa"][:, None])[:, :, None] + pmf = sp_betabinom.pmf(kk, n_trials, a, b).mean(axis=0) + return {lvl: _coverage_from_pmf(pmf, data["y"], lvl) for lvl in (0.5, 0.9)} + + +def _ame_items(observation, dr, data): + """Average marginal effect on the items scale, toggling G over all analysis rows.""" + eta1, eta0 = _eta_by_arm(dr, data) + n_trials = data["n_trials"] + if observation == "cmb": + m1 = (_cmb_pmf(eta1, dr["nu"], n_trials) * np.arange(n_trials + 1)).sum(axis=2) + m0 = (_cmb_pmf(eta0, dr["nu"], n_trials) * np.arange(n_trials + 1)).sum(axis=2) + return (m1 - m0).mean(axis=1) + sig1, sig0 = 1 / (1 + np.exp(-eta1)), 1 / (1 + np.exp(-eta0)) + if observation == "bb_chance_floor": + sig1, sig0 = 1 / 3 + (2 / 3) * sig1, 1 / 3 + (2 / 3) * sig0 + return n_trials * (sig1 - sig0).mean(axis=1) + + +def run_pair(symbol, variants, labels, args): + results = {} + for obs, label in zip(variants, labels, strict=True): + print(f"\n--- fitting {label} ({obs}) for outcome {symbol} ---") + data, idata = _fit(symbol, obs, args.draws, args.tune, args.chains, + args.target_accept, args.sampler) + print(f" n_obs = {len(data['y'])}, n_trials = {data['n_trials']}; {_diagnostics(idata)}") + names = ["alpha", "tau", "gamma_own", "gamma_A"] + (["nu"] if obs == "cmb" else ["kappa"]) + dr = _draws(idata, names) + print(" " + _summary(dr["tau"], "tau (logit)")) + print(" " + _summary(_ame_items(obs, dr, data), "AME (items)")) + print(" " + _summary(dr["gamma_A"], "gamma_A")) + print(" " + _summary(dr["gamma_own"], "gamma_own")) + disp = "nu" if obs == "cmb" else "kappa" + print(" " + _summary(dr[disp], disp)) + cov = _ppc_coverage(obs, dr, data) + print(f" predictive coverage: 50% band {cov[0.5]:.3f}, 90% band {cov[0.9]:.3f}") + results[label] = (data, idata) + both_ll = [id_ for (_, id_) in results.values() if hasattr(id_, "log_likelihood")] + if len(both_ll) == len(results) == 2: + try: + comp = az.compare({lbl: id_ for lbl, (_, id_) in results.items()}) + cols = [c for c in ("rank", "elpd_loo", "elpd_diff", "dse", "p_loo") if c in comp.columns] + print("\nPSIS-LOO comparison:") + print(comp[cols].to_string()) + except Exception as exc: # pragma: no cover - ArviZ version dependent + print(f"\n[PSIS-LOO comparison unavailable: {exc!r}]") + return results + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("task", choices=["blending-chance-floor", "word-reading-cmb"]) + ap.add_argument("--draws", type=int, default=2000) + ap.add_argument("--tune", type=int, default=1500) + ap.add_argument("--chains", type=int, default=4) + ap.add_argument("--target-accept", type=float, default=0.92) + ap.add_argument("--sampler", choices=["nutpie", "pymc"], default="nutpie") + args = ap.parse_args() + + if args.task == "blending-chance-floor": + run_pair("B", ["bb", "bb_chance_floor"], + ["standard Beta-Binomial", "chance-floor Beta-Binomial"], args) + else: + run_pair("W", ["bb", "cmb"], + ["standard Beta-Binomial", "Conway-Maxwell-binomial"], args) + + +if __name__ == "__main__": + main() diff --git a/scripts/ppc_coverage_sweep.py b/scripts/ppc_coverage_sweep.py new file mode 100644 index 00000000..5a77e6ac --- /dev/null +++ b/scripts/ppc_coverage_sweep.py @@ -0,0 +1,99 @@ +# Copyright (c) 2026 Down Syndrome Education International and contributors +# SPDX-License-Identifier: AGPL-3.0-or-later + +"""Pool posterior-predictive coverage (``ppc_summary.csv``) across fitted models. + +Motivation (notes/202607261405-binomial-exchangeability-item-difficulty-review.md): +heterogeneous item difficulty makes a bounded count *conditionally* underdispersed +relative to the Binomial (a Poisson-binomial has at-or-below-binomial variance given +ability), which the Beta-Binomial cannot express — its variance floor is the Binomial. +The observable symptom is predictive OVERcoverage: 50 % / 90 % prediction bands +covering more than 50 % / 90 % of observations. This script pools the per-fit +``ppc_summary.csv`` coverage rows by outcome symbol and family so the suite-level +pattern is visible at a glance, without any refitting. + +Small-denominator caveat: central intervals of a discrete count distribution +overcover mechanically (a "50 %" interval on a 10-item score covers at least 50 %), +so compare measures of similar length with each other rather than reading any single +row against an absolute nominal level. + +Usage: + python scripts/ppc_coverage_sweep.py + python scripts/ppc_coverage_sweep.py --models-dir /path/to/output/statistical_models/models + python scripts/ppc_coverage_sweep.py --config reporting --out output/statistical_models/comparison/ppc_coverage.csv +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +import pandas as pd + + +def collect(models_dir: Path, config: str) -> pd.DataFrame: + rows = [] + skipped = [] + for model_dir in sorted(models_dir.glob(f"*-{config}")): + ppc_path = model_dir / "ppc_summary.csv" + cfg_path = model_dir / "config.json" + if not ppc_path.exists() or not cfg_path.exists(): + skipped.append(model_dir.name) + continue + cfg = json.loads(cfg_path.read_text()) + ppc = pd.read_csv(ppc_path) + ppc["model_id"] = cfg.get("model_id", model_dir.name) + ppc["kind"] = cfg.get("kind", "?") + ppc["outcome_symbol"] = cfg.get("outcome_symbol") or "-" + rows.append(ppc) + if skipped: + print(f"[skipped {len(skipped)} dirs without ppc_summary.csv/config.json]") + if not rows: + raise SystemExit(f"No ppc_summary.csv found under {models_dir} for config {config!r}") + return pd.concat(rows, ignore_index=True) + + +def pooled(df: pd.DataFrame, keys: list[str]) -> pd.DataFrame: + grp = df.groupby(keys + ["level_pct"], dropna=False) + out = grp.agg( + n_models=("model_id", "nunique"), + n_total=("n_total", "sum"), + n_inside=("n_inside", "sum"), + mean_coverage=("coverage", "mean"), + ).reset_index() + out["pooled_coverage"] = out["n_inside"] / out["n_total"] + out["excess"] = out["pooled_coverage"] - out["level_pct"] / 100.0 + return out.sort_values(["level_pct", "excess"], ascending=[True, False]) + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument( + "--models-dir", + type=Path, + default=Path("output/statistical_models/models"), + help="Directory holding {model_id}-{config} output folders", + ) + ap.add_argument("--config", default="reporting", help="Config suffix to sweep (default: reporting)") + ap.add_argument("--out", type=Path, default=None, help="Optional CSV path for the long per-model table") + args = ap.parse_args() + + df = collect(args.models_dir, args.config) + print(f"Collected {len(df)} coverage rows from {df['model_id'].nunique()} models ({args.config}).") + + if args.out is not None: + args.out.parent.mkdir(parents=True, exist_ok=True) + df.to_csv(args.out, index=False) + print(f"Wrote long table to {args.out}") + + for mode, sub in df.groupby("mode"): + print(f"\n=== mode: {mode} — pooled coverage by outcome symbol ===") + tab = pooled(sub, ["outcome_symbol"]) + print(tab.to_string(index=False, float_format=lambda v: f"{v:0.3f}")) + print(f"\n=== mode: {mode} — pooled coverage by family (kind) ===") + print(pooled(sub, ["kind"]).to_string(index=False, float_format=lambda v: f"{v:0.3f}")) + + +if __name__ == "__main__": + main() diff --git a/src/language_reading_predictors/statistical_models/diagnostics.py b/src/language_reading_predictors/statistical_models/diagnostics.py index 34b8f37a..f2dc4145 100644 --- a/src/language_reading_predictors/statistical_models/diagnostics.py +++ b/src/language_reading_predictors/statistical_models/diagnostics.py @@ -60,6 +60,9 @@ from language_reading_predictors.statistical_models.context import ( StatisticalFitContext, ) +from language_reading_predictors.statistical_models.sampling_quality import ( + sampling_quality as _sampling_quality, +) from language_reading_predictors.statistical_models.plotting import ( save_plotcollection, save_styled_figure, @@ -571,22 +574,18 @@ def subfit_convergence(trace, *, label: str, var_names: list[str] | None = None) "n_divergences": None, } try: - # ``round_to="none"`` (the string) genuinely disables rounding; ``round_to=None`` - # falls through to ``rcParams["stats.round_to"]`` (2 sig figs) and would silently - # gate on rounded R-hat/ESS — the dseinternational/research#65 bug this check must - # not reproduce (it advertises "unrounded" signals above). - summ = az.summary( - trace, var_names=var_names, round_to="none", kind="diagnostics" - ) - max_rhat = float(summ["r_hat"].max()) - min_ess = float(min(summ["ess_bulk"].min(), summ["ess_tail"].min())) - n_div = int(np.asarray(trace.sample_stats["diverging"].values).sum()) - bfmi = _bfmi_per_chain(trace) - min_bfmi = ( - float(np.min(bfmi)) - if bfmi is not None and np.all(np.isfinite(bfmi)) - else None - ) + # Unrounded extraction lives in ``sampling_quality`` — see that module for the + # ``round_to="none"`` and coercion traps it exists to stop recurring. + signals = _sampling_quality(trace, var_names=var_names) + max_rhat = signals.max_rhat + min_ess = signals.min_ess + min_bfmi = signals.min_bfmi + if signals.n_divergences is None: + # No ``diverging`` in sample_stats: the gate cannot be evaluated, which is + # the "uncheckable" case (``converged=None``), not a failure. Previously the + # missing key raised and landed in the except branch below. + raise KeyError("sample_stats has no 'diverging' variable") + n_div = signals.n_divergences result.update( max_rhat=max_rhat, min_ess=min_ess, diff --git a/src/language_reading_predictors/statistical_models/loo_refit.py b/src/language_reading_predictors/statistical_models/loo_refit.py index e06bca6f..0ea8f8bd 100644 --- a/src/language_reading_predictors/statistical_models/loo_refit.py +++ b/src/language_reading_predictors/statistical_models/loo_refit.py @@ -35,7 +35,6 @@ from dataclasses import dataclass from typing import Any -import arviz as az import numpy as np import pymc as pm from arviz_stats.loo.wrapper import SamplingWrapper @@ -43,12 +42,14 @@ BFMI_THRESHOLD, ESS_THRESHOLD, RHAT_MAX, - _bfmi_per_chain, ) from language_reading_predictors.statistical_models import mechanism as _mechanism from language_reading_predictors.statistical_models.factories import _subset from language_reading_predictors.statistical_models.preprocessing import PreparedData +from language_reading_predictors.statistical_models.sampling_quality import ( + sampling_quality, +) __all__ = ["MechanismSamplingWrapper", "RefitPlan", "build_mechanism_wrapper"] @@ -195,20 +196,25 @@ def sample(self, modified_observed_data): return idata def _assert_refit_converged(self, idata) -> None: - """Fail the refit unless it clears the suite's sampling-quality thresholds.""" - summary = az.summary(idata, round_to=None) - max_rhat = float(np.nanmax(summary["r_hat"].to_numpy())) - min_ess = float(np.nanmin(summary["ess_bulk"].to_numpy())) - divergences = int(np.asarray(idata.sample_stats["diverging"].values).sum()) - # Use the shared per-chain helper: ``az.bfmi`` returns a DataTree in ArviZ 1.x, - # which cannot be coerced to an array. - bfmi = _bfmi_per_chain(idata) - min_bfmi = ( - float(np.min(bfmi)) if bfmi is not None and np.all(np.isfinite(bfmi)) else np.inf - ) + """Fail the refit unless it clears the suite's sampling-quality thresholds. + + Signals come from :func:`sampling_quality` so the refit gate reads them exactly + as the fit-time gate does. This previously called ``az.summary(round_to=None)``, + which rounds to two significant figures: every R-hat from 1.011 to 1.049 became + ``1.0`` and cleared the ``<= 1.01`` threshold, so the R-hat arm of this gate was + effectively ``< 1.05``. It also took ESS from ``ess_bulk`` alone, where the gate + takes the bulk/tail minimum. + """ + signals = sampling_quality(idata) + max_rhat = signals.max_rhat + min_ess = signals.min_ess + divergences = signals.n_divergences + # A missing BFMI is not treated as a failure here (unlike the fit-time sub-fit + # check), preserving this gate's original policy. + min_bfmi = np.inf if signals.min_bfmi is None else signals.min_bfmi failures = [] - if divergences > GATE_MAX_DIVERGENCES: + if divergences is None or divergences > GATE_MAX_DIVERGENCES: failures.append(f"{divergences} divergences") if not np.isfinite(max_rhat) or max_rhat > RHAT_MAX: failures.append(f"max R-hat {max_rhat:.4f}") diff --git a/src/language_reading_predictors/statistical_models/sampling_quality.py b/src/language_reading_predictors/statistical_models/sampling_quality.py new file mode 100644 index 00000000..77587cb4 --- /dev/null +++ b/src/language_reading_predictors/statistical_models/sampling_quality.py @@ -0,0 +1,106 @@ +# Copyright (c) 2026 Down Syndrome Education International and contributors +# SPDX-License-Identifier: AGPL-3.0-or-later + +"""One correct way to read sampling-quality signals off a trace. + +Every gate and every ad-hoc script needs the same four numbers — max R-hat, min ESS, +min per-chain BFMI, total divergences — and each one that re-derives them can get them +subtly wrong. Two traps, both observed in this repository: + +**Rounding.** ``az.summary()`` rounds to ``rcParams["stats.round_to"]`` (``"2g"`` — two +significant figures) unless passed ``round_to="none"``, *the string*. ``round_to=None`` +and ``"auto"`` both fall through to the rounded default, and omitting the argument +entirely returns a **string**-dtype frame that raises on ``float`` formatting. Rounding +erases exactly the digits the gate turns on: across the whole gate-relevant range every +R-hat from 1.011 to 1.049 rounds to ``1.0`` and clears an ``R-hat <= 1.01`` test it +should fail (dseinternational/research#65; recurred twice — in ``loo_refit`` and in a +one-off prototype script, issue #440). + +**Coercion.** ``trace.sample_stats["diverging"]`` is an xarray ``DataArray``; reduce it +via ``np.asarray(....values).sum()`` rather than relying on ``DataArray.__int__``, and +take BFMI from :func:`dse_research_utils.statistics.diagnostics._bfmi_per_chain`, since +``az.bfmi`` returns a ``DataTree`` in ArviZ 1.x that cannot be coerced to an array. + +This module extracts the numbers and nothing else. It deliberately does **not** decide +whether a fit passed: the call sites differ in which variables they gate over and in how +they treat a missing BFMI, and those policies stay where they are rather than being +silently homogenised here. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import arviz as az +import numpy as np +from dse_research_utils.statistics.diagnostics import _bfmi_per_chain + +__all__ = ["SamplingQuality", "sampling_quality"] + + +@dataclass(frozen=True) +class SamplingQuality: + """Unrounded sampling-quality signals for one trace.""" + + max_rhat: float + """Largest R-hat over the summarised variables (NaNs skipped).""" + min_ess: float + """Smallest of bulk-ESS and tail-ESS over the summarised variables.""" + min_bfmi: float | None + """Smallest per-chain BFMI, or ``None`` when it cannot be computed.""" + n_divergences: int | None + """Total divergent transitions, or ``None`` when ``sample_stats`` lacks them.""" + + def summary_line(self) -> str: + """One-line human-readable rendering for logs and prototype scripts.""" + bfmi = "n/a" if self.min_bfmi is None else f"{self.min_bfmi:.2f}" + div = "n/a" if self.n_divergences is None else str(self.n_divergences) + return ( + f"max R-hat {self.max_rhat:.4f}, min ESS {self.min_ess:.0f}, " + f"min BFMI {bfmi}, divergences {div}" + ) + + +def sampling_quality(trace, *, var_names: list[str] | None = None) -> SamplingQuality: + """Read the four sampling-quality signals off ``trace``, unrounded. + + Parameters + ---------- + trace + An ArviZ ``InferenceData`` (or DataTree-backed equivalent) with a ``posterior`` + group; ``sample_stats`` is used for divergences and BFMI when present. + var_names + Restrict the R-hat / ESS summary to these variables. ``None`` summarises + everything ArviZ reports for the trace, which includes deterministics — pass the + caller's curated gate variables when that matters. + + Returns + ------- + SamplingQuality + The extracted signals. Exceptions from ArviZ propagate; callers that must + tolerate a failed diagnostic calculation should catch them and decide what an + uncheckable fit means for them. + """ + # ``round_to="none"`` must be the string — see the module docstring. + summ = az.summary(trace, var_names=var_names, round_to="none", kind="diagnostics") + # pandas ``.max()`` / ``.min()`` skip NaN by default, so a constant or unsampled + # variable does not poison the reduction. + max_rhat = float(summ["r_hat"].max()) + min_ess = float(min(summ["ess_bulk"].min(), summ["ess_tail"].min())) + + n_div: int | None = None + sample_stats = getattr(trace, "sample_stats", None) + if sample_stats is not None and "diverging" in sample_stats: + n_div = int(np.asarray(sample_stats["diverging"].values).sum()) + + bfmi = _bfmi_per_chain(trace) + min_bfmi = ( + float(np.min(bfmi)) if bfmi is not None and np.all(np.isfinite(bfmi)) else None + ) + + return SamplingQuality( + max_rhat=max_rhat, + min_ess=min_ess, + min_bfmi=min_bfmi, + n_divergences=n_div, + ) diff --git a/tests/statistical_models/test_diagnostics.py b/tests/statistical_models/test_diagnostics.py index 9dda974c..cf9c6a6a 100644 --- a/tests/statistical_models/test_diagnostics.py +++ b/tests/statistical_models/test_diagnostics.py @@ -15,6 +15,9 @@ import xarray as xr from language_reading_predictors.statistical_models import diagnostics as diag +from language_reading_predictors.statistical_models import ( + sampling_quality as sampling_quality_mod, +) def test_run_psense_removes_stale_summary_when_recomputation_fails( @@ -610,7 +613,11 @@ def test_subfit_convergence_catches_bad_nuisance_parameter(): def test_subfit_convergence_flags_low_bfmi(monkeypatch): - monkeypatch.setattr(diag, "_bfmi_per_chain", lambda _trace: np.asarray([0.2, 0.8])) + # BFMI is now read by the shared ``sampling_quality`` extractor, so that is the + # seam to patch; ``diag._bfmi_per_chain`` is no longer on this call path. + monkeypatch.setattr( + sampling_quality_mod, "_bfmi_per_chain", lambda _trace: np.asarray([0.2, 0.8]) + ) result = diag.subfit_convergence( _synthetic_trace(0.0), label="low-bfmi", var_names=["tau"] ) diff --git a/tests/statistical_models/test_sampling_quality.py b/tests/statistical_models/test_sampling_quality.py new file mode 100644 index 00000000..55b705f8 --- /dev/null +++ b/tests/statistical_models/test_sampling_quality.py @@ -0,0 +1,162 @@ +# Copyright (c) 2026 Down Syndrome Education International and contributors +# SPDX-License-Identifier: AGPL-3.0-or-later + +"""Guards for the shared unrounded sampling-quality extraction (#440). + +The point of ``sampling_quality`` is that there is exactly one place where R-hat, ESS, +BFMI and divergences are read off a trace, and that place does not round. The critical +regression these tests pin is the ``round_to`` trap: ``az.summary`` rounds to two +significant figures unless passed the *string* ``"none"``, which turns every R-hat from +1.011 to 1.049 into ``1.0`` and silently clears an ``R-hat <= 1.01`` gate +(dseinternational/research#65; recurred in ``loo_refit`` and in a prototype script). +""" + +from __future__ import annotations + +import arviz as az +import numpy as np +import pytest +import xarray as xr + +from language_reading_predictors.statistical_models.sampling_quality import ( + SamplingQuality, + sampling_quality, +) + + +def _trace( + *, + n_chains: int = 4, + n_draws: int = 400, + divergences: int = 0, + seed: int = 0, + offset: float = 0.0, + with_energy: bool = True, + with_diverging: bool = True, + extra_vars: dict | None = None, +): + """A small DataTree trace with well-mixed chains and a known divergence count. + + ``offset`` shifts chain 0's mean to manufacture a deliberately poor R-hat. + """ + rng = np.random.default_rng(seed) + draws = rng.normal(size=(n_chains, n_draws)) + draws[0] += offset + + data = {"theta": (("chain", "draw"), draws)} + for name, values in (extra_vars or {}).items(): + data[name] = (("chain", "draw"), values) + posterior = xr.Dataset( + data, coords={"chain": np.arange(n_chains), "draw": np.arange(n_draws)} + ) + + stats: dict = {} + if with_diverging: + div = np.zeros((n_chains, n_draws), dtype=bool) + if divergences: + div.reshape(-1)[:divergences] = True + stats["diverging"] = (("chain", "draw"), div) + if with_energy: + stats["energy"] = (("chain", "draw"), rng.normal(size=(n_chains, n_draws))) + + groups = {"posterior": posterior} + if stats: + groups["sample_stats"] = xr.Dataset( + stats, coords={"chain": np.arange(n_chains), "draw": np.arange(n_draws)} + ) + return xr.DataTree.from_dict(groups) + + +# --- the regression this module exists for ---------------------------------------- + + +def test_reports_unrounded_rhat_where_rounding_would_hide_the_failure(): + # Tuned to land in the (1.01, 1.05) band — the range where two-significant-figure + # rounding collapses a gate failure onto a clean-looking 1.0. + trace = _trace(offset=0.4, seed=2) + signals = sampling_quality(trace) + rounded = float(az.summary(trace, kind="diagnostics", round_to=None)["r_hat"].max()) + + assert signals.max_rhat > 1.01, "fixture should fail an unrounded gate" + # The bug: the rounded value clears the very gate the true value fails. + assert rounded <= 1.01 + assert signals.max_rhat != rounded + + +@pytest.mark.parametrize("true_rhat", [1.011, 1.02, 1.049]) +def test_two_significant_figure_rounding_collapses_the_gate_band(true_rhat): + """Every R-hat in (1.01, 1.05) rounds to 1.0 — the gate becomes ``< 1.05``.""" + assert true_rhat > 1.01, "must fail an unrounded R-hat <= 1.01 gate" + assert float(f"{true_rhat:.2g}") == 1.0, "but passes once rounded to 2 sig figs" + + +def test_default_summary_is_string_typed(): + """Omitting ``round_to`` returns strings, which raise on float formatting.""" + col = az.summary(_trace(), kind="diagnostics")["r_hat"] + assert col.dtype == object or col.dtype.kind in {"U", "O"} + with pytest.raises((TypeError, ValueError)): + f"{col.max():.4f}" + # The helper is immune. + assert isinstance(sampling_quality(_trace()).max_rhat, float) + + +# --- extraction behaviour ---------------------------------------------------------- + + +def test_min_ess_takes_the_bulk_tail_minimum(): + trace = _trace() + summ = az.summary(trace, round_to="none", kind="diagnostics") + expected = min(float(summ["ess_bulk"].min()), float(summ["ess_tail"].min())) + assert sampling_quality(trace).min_ess == pytest.approx(expected) + + +def test_counts_divergences_via_array_coercion(): + assert sampling_quality(_trace(divergences=7)).n_divergences == 7 + assert sampling_quality(_trace(divergences=0)).n_divergences == 0 + + +def test_divergences_none_when_sample_stats_lacks_them(): + trace = _trace(with_diverging=False, with_energy=False) + assert sampling_quality(trace).n_divergences is None + + +def test_var_names_restricts_the_summary(): + rng = np.random.default_rng(1) + nuisance = rng.normal(size=(4, 400)) + nuisance[0] += 5.0 # badly mixed + trace = _trace(extra_vars={"nuisance": nuisance}) + + assert sampling_quality(trace, var_names=["theta"]).max_rhat < 1.01 + assert sampling_quality(trace).max_rhat > 1.01 + + +def test_nan_diagnostics_do_not_poison_the_reduction(): + """A constant (unsampled) variable yields NaN diagnostics; those are skipped.""" + trace = _trace(extra_vars={"constant": np.ones((4, 400))}) + assert np.isfinite(sampling_quality(trace).max_rhat) + + +def test_bfmi_absent_without_energy(): + assert sampling_quality(_trace(with_energy=False)).min_bfmi is None + + +def test_bfmi_finite_when_energy_present(): + min_bfmi = sampling_quality(_trace()).min_bfmi + assert min_bfmi is None or np.isfinite(min_bfmi) + + +# --- rendering --------------------------------------------------------------------- + + +def test_summary_line_renders_values_and_missing_alike(): + line = SamplingQuality( + max_rhat=1.0018, min_ess=4137.4, min_bfmi=None, n_divergences=None + ).summary_line() + assert "1.0018" in line, "R-hat must keep four decimals, not round to 1.00" + assert "4137" in line + assert line.count("n/a") == 2 + + full = SamplingQuality( + max_rhat=1.0, min_ess=8000.0, min_bfmi=0.91, n_divergences=3 + ).summary_line() + assert "0.91" in full and "divergences 3" in full