Hierarchical Models for Clustered Assignment

Python · PyMC · ArviZ  ·  Project STAR, 3,743 students in 79 schools

What a Cluster-Robust Standard Error Cannot Tell You

The clustered-designs notebook established the problem on Project STAR: 3,743 students across 79 schools, where children in the same school share a teacher, a building and a peer group. Ignoring that correlation gives a standard error of 1.042; clustering by school gives 1.850, and a known-truth simulation showed that naive intervals cover only 54% of the time. The cluster-robust standard error is the right fix, and it is a patch: it corrects one number and tells you nothing else. It cannot say how much schools differ, cannot produce an estimate for any individual school, and cannot address whether small classes help equally everywhere — which for a policy rolled out school by school is the question that decides whether it travels.

readi=αs[i]+τsmalli+εi,αsN(μα, σα2)\text{read}_i = \alpha_{s[i]} + \tau\,\text{small}_i + \varepsilon_i, \qquad \alpha_s \sim \mathcal{N}(\mu_\alpha,\ \sigma_\alpha^2)

A hierarchical model puts the school in the model rather than in a variance correction. School intercepts are drawn from a common distribution whose spread σα\sigma_\alpha is estimated, and the intraclass correlation follows directly as σα2/(σα2+σy2)\sigma_\alpha^2/(\sigma_\alpha^2+\sigma_y^2) — with a full posterior rather than a point proxy. The result is 0.214 with a 95% interval of [0.160, 0.279], so "strong clustering" finally acquires error bars.

Does the Effect Itself Vary by School?

The version above still forces the small-class benefit to be identical in every school. Adding a varying slope tests that, and it is where the notebook stops confirming and starts contradicting.

readi=αs[i]+τs[i]smalli+εi,τsN(μτ, στ2)\text{read}_i = \alpha_{s[i]} + \tau_{s[i]}\,\text{small}_i + \varepsilon_i, \qquad \tau_s \sim \mathcal{N}(\mu_\tau,\ \sigma_\tau^2)

The between-school spread of the effect comes out at στ=11.9\sigma_\tau = 11.9 points, with 95% interval [9.24, 14.85] and Pr(στ<5)=0.000\Pr(\sigma_\tau < 5) = 0.000. That is larger than the average effect itself (6.57). School-level effects run from −16.0 to +37.5, and 18 of 79 schools have intervals excluding zero where chance alone would give about four. Small classes do not help equally everywhere — and nothing in the published analysis could have revealed that, because its entire treatment of clustering was a correction to one standard error.

One caveat belongs with that finding. The median school contributes only 43 students, so any individual school's effect is poorly determined on its own — which is exactly why the model partially pools them. στ\sigma_\tau is identified from the spread across schools rather than from any single one, and a spread this large is hard to manufacture from noise at 79 schools. The defensible reading is that the effect genuinely varies; the indefensible one is to trust any particular school's number.

quantityposteriorreading
ICC0.214, 95% [0.160, 0.279]the point proxy, now with an interval
μτ average effect6.57, 95% [3.33, 9.92]
στ spread across schools11.89, 95% [9.24, 14.85]larger than the average effect itself
school-level effects−16.0 to +37.518 of 79 exclude zero; chance would give ~4

The Width Comparison, Which Reverses the Moral

Now the result that reverses the expected moral. Line the four analyses up by interval width rather than by sophistication:

analysiseffect95% intervalwidth
OLS, naive SE5.82[3.78, 7.86]4.08 known to under-cover
OLS, cluster-robust SE5.82[2.19, 9.44]7.25
hierarchical, varying intercept6.54[4.63, 8.38]3.74 tightest — and misspecified
hierarchical, + varying slope6.57[3.33, 9.92]6.59

The varying-intercept model is the tightest of the four — width 3.74 against cluster-robust's 7.25 — and that tightness is not an improvement. It is bought by assuming one common treatment effect, an assumption the varying-slope model then demolishes. Allow the effect to vary and the interval widens back to 6.59, close to cluster-robust's.

So the cluster-robust standard error was right to be wide. Making no assumption about how schools differ, it quietly absorbed heterogeneity that the simpler hierarchical model wished away. The usual framing — robust standard errors as a crude patch, hierarchical models as the principled alternative — gets this backwards. A hierarchical model beats a robust standard error only when its structure is correct; here the first version was tighter and wrong, and it took a second, richer model to notice. A robust standard error is the more conservative choice precisely because it declines to model what it does not know.

A smaller point worth noting: the hierarchical estimates (6.54, 6.57) sit above the OLS estimate of 5.82. The hierarchical models estimate the effect within schools, holding each school's level fixed, while OLS pools across them. The gap says small classes were not distributed evenly across schools of differing baseline achievement — a between-school component that the within-school comparison removes.

Where this sits

The honest summary is that these tools answer different numbers of questions rather than competing on one. Cluster-robust inference is the right choice when the effect is all you want and you distrust your model, since it assumes nothing about how schools vary. The hierarchical model asks for a distributional assumption and returns the ICC with uncertainty, seventy-nine school-level effects, and a testable constant-effect assumption — but it is only as good as that assumption, which is exactly what the width comparison demonstrates. This is the same partial-pooling machinery as the Electric Company example, scaled from four grades to seventy-nine schools and applied to grouping that was a nuisance rather than a design feature.

Notebook

Downloads

References