Model checking for the educational testing example#

Part 2 · Stage 6 · 🔍 Model Checking & Comparison · Lesson 044 of 144 · intermediate

◀ Previous · Graphical posterior predictive checks · Next · Measures of predictive accuracy ▶ · ↑ Section

Important

✨ AI-generated content. This page was written with the assistance of an AI language model and is provided as a learning aid. Despite careful review, it may still contain mistakes, omissions, or out-of-date information. Whether you are new to the topic, a team lead, or a senior practitioner, treat it as a starting point rather than an authoritative reference: read it critically and independently verify anything you act on (code, commands, figures, and factual claims) against official documentation and primary sources before relying on it.

Checking eight schools#

The hierarchical model of Stage 5 gave a satisfying answer: school A’s dramatic 28-point effect shrinks toward the pooled mean, and \(\tau\) sits near zero. Satisfying is not the same as correct. Step three demands that we ask whether data generated by this fitted model would look like the eight estimates actually reported.

Generating replicated schools#

For each posterior draw \((\mu^{(s)}, \tau^{(s)}, \theta^{(s)})\), simulate a complete replicated experiment — eight new estimates, each with the same known standard errors:

\[y^{\text{rep}(s)}_j \sim \mathrm{N}\bigl(\theta_j^{(s)},\; \sigma_j^2\bigr), \qquad j = 1, \dots, 8 .\]

Now choose test quantities the model does not automatically match. The largest of the eight effects is the natural one, since the largest observed value (28, from school A) is what motivated the analysis in the first place. The smallest and the range (largest minus smallest) probe the spread.

import numpy as np
yrep = ...                       # shape (S, 8) replicated datasets
for name, T in [("max", np.max), ("min", np.min),
                ("range", lambda d: d.max() - d.min())]:
    p_B = (np.array([T(r) for r in yrep]) >= T(y)).mean()
    print(f"{name:6s}  observed {T(y):6.1f}   p_B = {p_B:.2f}")

What the checks say#

The observed maximum, minimum and range all fall comfortably inside their replicated distributions: these \(p\)-values are unremarkable, and there is no evidence that the hierarchical model misfits the data in these respects. Specifically, a largest-of-eight value of 28 is entirely ordinary for eight noisy estimates drawn around a common mean with these standard errors. That is the whole point: the “striking” result that provoked the study is exactly what the model predicts should happen by chance.

What a check cannot do#

Passing does not make the model true. With \(J = 8\) groups the checks have little power: many models would survive them, including ones with quite different \(\tau\) priors. And the check tests the likelihood and hierarchy, not the assumption that the eight schools are exchangeable or that \(\sigma_j\) are truly known. Model checking bounds what the data can refute; it does not certify the assumptions the data cannot see. The honest next step is sensitivity analysis — refit under a \(t\) population distribution, or a different prior on \(\tau\), and see whether any conclusion you would act on moves.

Hint