DiD review: resolve estimand, diagnostic, and maintainability gaps

[frankbuckley]

Note

Drafted by a LLM-based AI tool (Codex/GPT-5).

Summary

A critical review of the waitlist-crossover/DiD statistical-model family found that the core arm-by-wave model is implemented coherently and is statistically defensible for its principal randomised estimand. However, the family should not receive an unconditional methodological sign-off yet: two estimand/reporting defects, an empirically anchored prior that needs an explicit methodological decision, an invalid current model comparison, and several diagnostic/audit gaps require follow-up.

The model is more accurately described as a saturated longitudinal arm-by-wave mixed model than as a classical difference-in-differences model. The clean causal quantity is tau_t2, the absolute immediate-versus-waitlist arm contrast at t2. arm_gap_t1 is a baseline-balance quantity. arm_gap_t3 and delta_crossover = tau_t2 - arm_gap_t3 are post-crossover associations, not randomised causal effects. Dose coefficients are adjusted observational associations. These distinctions are generally present in the implementation and reports and must be preserved.

Priority findings

P1 — Make “new-child” predicted probabilities and scores target the same population

In src/language_reading_predictors/statistical_models/predicted_scores.py, prob_control, prob_intervention, and ame_prob are calculated from the linear predictor before the child random effect is replaced or integrated for a new child. The Beta-Binomial score simulation subsequently applies the new-child adjustment, whereas the Bernoulli path returns before that adjustment. The resulting predicted_scores.csv therefore mixes population-integrated score simulations with fitted-child conditional probability effects; off-floor Bernoulli output is entirely fitted-child conditional despite being described as a “new typical child” result.

Recommended action: define the target population explicitly and apply it consistently. For new-child output, integrate over the child-effect distribution before calculating marginal probabilities and scores, using quadrature or Monte Carlo. If both estimands are scientifically useful, emit separately named finite-sample conditional and new-child marginal results. Add regression tests for both Beta-Binomial and Bernoulli paths.

P1 — Decide how to handle the outcome-informed intercept prior

In src/language_reading_predictors/statistical_models/factories.py, the pooled t1 outcome is used to construct a logit-scale intercept anchor, and the same t1 observations then enter the fitted likelihood. This is an empirical-Bayes prior, not a fully independent prior. It also means the reported prior-predictive distribution is partly informed by the observed outcome data.

Recommended action: preferably replace the anchor with a genuinely external or weakly informative scale-aware prior. If the empirical anchor is retained, label it explicitly as empirical Bayes in the model recipe and report, explain its inferential consequences, and add a sensitivity fit using an independent prior. This may have limited practical influence on tau_t2, but that must be demonstrated rather than assumed.

P1 — Do not interpret the current DID-007 versus DID-107 PSIS-LOO comparison

The reviewed reporting fit for DID-007 failed the convergence gate with 14 divergences in 36,000 draws. The reviewed traces also had maximum Pareto-k values of approximately 1.135 for DID-007 and 1.152 for DID-107. Ordinary PSIS-LOO is unreliable at these values. The current comparison code in scripts/compare_statistical_models.py checks convergence and row identity before calling az.compare, but does not independently reject or prominently flag a comparison with unreliable Pareto-k diagnostics.

Recommended action: first resolve the DID-007 divergences. Then make the comparison gate inspect Pareto-k reliability and refuse to mark the comparison valid when importance sampling is unreliable. Use exact refits/reloo or grouped K-fold validation if the comparison remains scientifically important. Until then, the available evidence cannot establish whether the period-varying dose model improves predictive performance. The existing conditional-within-child interpretation of this LOO comparison should remain explicit.

P2 — Correct boundary behaviour in cell zero-rate posterior-predictive checks

In src/language_reading_predictors/statistical_models/reporting.py, the zero-rate check calculates the share of replicated zero rates greater than or equal to the observed value and flags either extreme tail. When the observed zero rate is exactly zero, that probability is necessarily one, so a well-fitting cell can be falsely flagged. The reviewed artefacts contain several such flags where the replicated median zero rate is also zero.

Recommended action: use a discrete two-sided tail probability, a mid-p-style diagnostic, or direct inclusion of the observation in a posterior-predictive central interval. Add boundary tests for observed zero rates of zero and one.

P2 — Extend prior-sensitivity checks to variant-defining parameters

The prior-sensitivity routing in src/language_reading_predictors/statistical_models/pipeline.py selects only tau_t2, beta_dose, or mu_dose. It therefore omits parameters that define important variants, including DID-007’s sigma_dose/beta_dose_phase and DID-013’s sigma_delta. The existing critical-review note also records tau_t2 prior sensitivity in seven DiD fits, with especially weak information for off-floor outcomes.

Recommended action: define a per-model list of scientifically central sensitivity parameters and run alternative-prior refits for each. Include the heterogeneous-effect scale and period-varying dose parameters where applicable. Treat DID-013 as a sensitivity analysis only: each waitlist deviation is informed by one t3 observation and must be separated from child heterogeneity, Beta-Binomial dispersion, and measurement noise; it cannot identify individual responders.

P3 — Make row exclusions auditable by reason

The generic subsetting logic in src/language_reading_predictors/statistical_models/factories.py increments dropped_rows for both missing-data exclusions and deliberate design exclusions. For example, the reviewed DID-001 configuration reports dropped_rows = 57, while its detailed contract distinguishes one loader exclusion, 53 design-excluded rows, and three factory missing-data exclusions. The detailed contract is useful, but the top-level field is misleading.

Recommended action: replace the overloaded count with typed, reason-specific counters and calculate any total from those components. Test that the components are mutually exclusive and reconcile to the source-row count.

Refactoring recommendations

• Introduce immutable, validated DidRunSettings and DidRunPlan types analogous to the ITT family. Resolve contradictory or unknown options before data loading and output reset, and make that plan the single input to model construction, reporting, influence refits, and metadata.
• Split build_did_model into arm-by-wave and dose-model factories. The two designs have different rows, estimands, and assumptions and should not be controlled by a growing collection of spec.extra flags.
• Move family-owned construction, estimand derivation, diagnostics, and audit logic from the very large factories.py and pipeline.py compatibility modules into a focused statistical_models/did.py module.
• Replace the untyped extras mapping with a DiD result/audit dataclass so required fields and meanings are visible to unfamiliar readers and checked statically.
• Remove the legacy delta reporting branch once stale artefact compatibility is no longer required.
• Consider using “waitlist-crossover arm-by-wave” as the primary reader-facing name, retaining “DiD” as a historical/search alias. This will reduce the risk that readers interpret tau_t2 as a classical change-in-changes coefficient.
• Add one small real-sampler parameter-recovery test for the production Beta-Binomial random-intercept model. Current statistical validation covers closed-form saturated Binomial behaviour and prior prediction, but not recovery in the actual hierarchical likelihood.

Verification completed during review

Focused DiD validation, factory, reporting, predicted-score, comparison, and key-findings tests passed, as did targeted Ruff checks. Those tests confirm the currently encoded behaviour but do not cover the defects above. Thirteen of the fourteen reviewed production fits passed the convergence gate; DID-007 did not. No source files were changed during the review.

Definition of done

• New-child probability and score estimands use a single explicitly documented target population, with Beta-Binomial and Bernoulli regression tests.
• The intercept prior is made independent of the fitted outcomes or is explicitly documented and sensitivity-tested as empirical Bayes.
• DID-007 passes the convergence gate before comparison, and the comparison layer rejects unreliable Pareto-k diagnostics or uses a reliable alternative.
• Cell zero-rate PPC checks behave correctly at zero and one.
• Prior sensitivity covers all variant-defining estimands and reports the existing tau_t2 sensitivity clearly.
• Row exclusions are recorded once and partitioned by reason.
• The DiD family has a validated run-plan boundary and clearer separation between arm-by-wave and dose models.
• A hierarchical parameter-recovery test exercises the production likelihood.
• Reader-facing documentation preserves the causal/associational distinctions above and avoids describing DID-013 as individual-response identification.

[ethanbuckley]

Note

Drafted by an LLM-based AI tool (Claude Code/Opus 4.8).

#390 progress + scope for the remaining findings

Shipped: the self-contained diagnostic fix — P2 cell zero-rate PPC boundary — in #403 (mid-p tail so a cell with observed zero-rate exactly 0 or 1 is no longer falsely flagged; boundary regression test added). The rest of #390 splits by what it needs, so I've kept it out of that PR:

Needs a methodological decision (Frank)

• P1 — outcome-informed intercept prior (empirical Bayes). The pooled t1 outcome anchors the logit intercept and the same t1 observations enter the likelihood, so the prior is empirical Bayes and the reported prior-predictive is partly data-informed. Decision: (a) replace with a genuinely external / weakly-informative scale-aware prior, or (b) retain the empirical anchor but label it explicitly as EB in the recipe/report and add an independent-prior sensitivity fit (a refit). Both options need a reporting-tier refit to land; (a) also changes the fitted intercept. My lean: (b) if the empirical anchor is deliberate and its influence on tau_t2 is demonstrably small, else (a).

Needs a refit / compute cycle

• P1 — DID-007 divergences. The reviewed DID-007 reporting fit failed the convergence gate (14 divergences) with max Pareto-k ≈ 1.14 (DID-107 ≈ 1.15). PSIS-LOO — and therefore the DID-007-vs-DID-107 elpd comparison — is unreliable until this is resolved (reparameterise / retune, then refit).
• P2 — variant-parameter sensitivity coverage. The psense routing selects only tau_t2/beta_dose/mu_dose, omitting DID-007's sigma_dose/beta_dose_phase and DID-013's sigma_delta. Wiring a per-model sensitivity list is easy engineering, but the alternative-prior refits it enables need a compute cycle — this overlaps Emit power-scaling (psense) and prior-pushforward artifacts for all estimand-reporting families #381 (emit psense/pushforward) / Prior-sensitivity refits for the prior-dependent estimands #382 (sensitivity refits), so best folded into that refit pass. (Keep DID-013 explicitly framed as a sensitivity analysis, not individual-response identification.)

Actionable engineering (follow-up PRs, not blocked)

• P1 — comparison Pareto-k gate. scripts/compare_statistical_models.py::_write_comparison already gates on the convergence gate + shared row identity before az.compare, but marks comparison_valid=True without checking Pareto-k reliability. Add a khat gate (refuse / flag comparison_valid=False when any model's max Pareto-k exceeds ~0.7). Self-contained; I left it out of fix(did): mid-p zero-rate PPC + Pareto-k comparison gate (#390 P2 + P1) #403 only to keep the arviz-loo khat handling on its own testable PR. This can land independently of resolving the DID-007 divergences (it would correctly mark the current comparison invalid).
• P1 — predicted_scores target population. predicted_scores.py mixes population-integrated score simulations (Beta-Binomial path) with fitted-child conditional probabilities (Bernoulli path returns before the new-child adjustment). Define one documented target population and apply it consistently (integrate over the child-effect distribution for new-child output), or emit separately named conditional vs new-child results; add Beta-Binomial + Bernoulli regression tests. Medium-sized, has a modeling-choice element worth confirming.

Overlaps #394 (keep separate — mechanical relocation, not behaviour)

• P3 — row-exclusion counters. dropped_rows is a shared PreparedData field consumed across preprocessing / pipeline metadata / reporting, so replacing it with typed reason-specific counters is a cross-family contract change, best done with the DiD family module extraction.
• The DidRunSettings/DidRunPlan typing, splitting build_did_model into arm-by-wave + dose factories, and moving DiD logic into statistical_models/did.py are the Complete the statistical pipeline family split and artefact lifecycle refactor #394 family-migration work.

Suggested next steps I can take without a refit or your ruling: the comparison Pareto-k gate and the predicted_scores population fix. The EB-prior decision and the DID-007 refit are yours/Frank's calls; the sensitivity-coverage and row-exclusion items ride the refit cycle / #394 respectively.

[ethanbuckley]

Note

Drafted by an LLM-based AI tool (Claude Code/Opus 4.8).

Update: the comparison Pareto-k gate (the actionable P1 engineering item above) is now included in #403 alongside the P2 zero-rate PPC fix — _loo_compare reads each fit's persisted pareto_k.csv and marks the comparison invalid (per-model elpd_loo, not deltas) when any max khat exceeds its good-k threshold, so the current DID-007/DID-107 comparison would be correctly flagged rather than presented as a valid ranking. Resolving DID-007's divergences remains a separate refit.

[ethanbuckley]

Note

Drafted by an LLM-based AI tool (Claude Code/Opus 4.8).

P1 predicted_scores target population — this is an estimand/reporting decision, not a clean fix

Investigating the P1 predicted_scores finding, the "mixed populations" turn out to be partly deliberate and guard-enforced, so I'm flagging it for a decision rather than changing it:

• PredictiveContrast.ame_prob / prob_control / prob_intervention are the sample-conditional headline AME. The docstring states "ame_prob … reproduce reporting._itt_ame_draws exactly — the guard test relies on this", and test_predicted_scores.py::test_ame_matches_itt_ame_draws asserts contrast.ame_prob == _itt_ame_draws(...) to atol=1e-12. So integrating the child effect into these (making them new-child marginal) would break that guard and make the AME diverge from the reported headline — not a bug fix.
• The Beta-Binomial score_* are new-child marginal (fresh u ~ Normal(0, sigma_child) per simulated child).
• So the artefact deliberately carries a conditional effect (matching the headline) and marginal score simulations. The genuine defects are (i) the module framing calls the whole thing a "new typical child" prediction, which is only true of the scores, and (ii) the off-floor Bernoulli path has no scores, so its only output is the conditional prob — yet it is presented as the "new typical child" result.

The decision (yours/Frank's)

• Option A — document the intentional split (no value change). State per-column that ame_prob/prob_* are the sample-conditional AME (reproducing the headline, guard-tested) while score_* are new-child marginal; fix the blanket "new typical child" framing and the off-floor caption to be column-precise. Minimal, preserves the guard test and all reported values.
• Option B — emit separately named new-child-marginal outputs. Keep the guard-tested conditional ame_prob, and add prob_*_new_child / ame_prob_new_child computed by integrating the child effect (quadrature/MC), so both estimands are explicit. Additive schema; the off-floor binary display would then show a genuine new-child-marginal rate. More work; needs your view on whether the marginal prob is scientifically wanted.

My recommendation: A if the conditional AME (reproducing the headline) is the intended primary quantity and only the framing needs fixing; B if a new-child-marginal event probability is genuinely wanted for the off-floor binary display. I can implement either on your ruling — I've deliberately not changed ame_prob, since it is a guard-tested contract and the choice is a scientific one.

With this, #390's self-contained non-refit engineering is complete (P2 zero-rate PPC + P1 comparison Pareto-k gate, both in #403); the remainder (this, the empirical-Bayes intercept prior, the DID-007 divergences, and the variant-parameter sensitivity coverage) all need either a decision or a refit.

[ethanbuckley]

Note

Drafted by an LLM-based AI tool (Claude Code/Opus 4.8).

Cross-link: the first P1 here (make the "new-child" predicted probabilities/scores target the same population) is resolved by #417. predicted_scores.py is shared across the randomised-contrast families, and #417 implements the recommended action directly: counterfactual_predictive_contrast now integrates the child random effect out (Gauss–Hermite quadrature) to produce a new-child, population-average effect for both the Beta-Binomial and the Bernoulli paths (the floor-rule path no longer returns before the new-child step), while the existing fitted-child conditional AME is kept and matches treatment_marginal.csv/rope_summary.csv. Both estimands are emitted under distinct quantity names with an explicit intercept_basis column, and there are regression tests for both likelihoods. So the DiD/waitlist-crossover reports get the fix automatically.

For the record, the two other landed items — P1 (Pareto-k LOO-reliability gating) and P2 (mid-p zero-rate PPC) — merged via #403. Remaining here: the empirical-Bayes intercept-prior decision (P1), extending prior-sensitivity to variant-defining parameters (P2), typed reason-specific dropped_rows counters (P3), and the DidRunSettings/did.py refactors.

[frankbuckley]

@ethanbuckley -- follow up -- option B with three conditions:

1. The sensitivity fit must use a genuinely independent prior, not a wider σ around the same anchor — widening σ leaves the mean data-dependent, which is the actual defect.
2. The report must state the prior-predictive limitation explicitly, not just the EB label. That's the cost you're accepting.
3. Apply the same label to growth and LCSM, or record why not — growth's anchor is the weaker case and will be found.