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.
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 is estimated, and the intraclass correlation follows directly as — 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.
The between-school spread of the effect comes out at points, with 95% interval [9.24, 14.85] and . 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. 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.
| quantity | posterior | reading |
|---|---|---|
| ICC | 0.214, 95% [0.160, 0.279] | the point proxy, now with an interval |
| μτ average effect | 6.57, 95% [3.33, 9.92] | — |
| στ spread across schools | 11.89, 95% [9.24, 14.85] | larger than the average effect itself |
| school-level effects | −16.0 to +37.5 | 18 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:
| analysis | effect | 95% interval | width |
|---|---|---|---|
| OLS, naive SE | 5.82 | [3.78, 7.86] | 4.08 known to under-cover |
| OLS, cluster-robust SE | 5.82 | [2.19, 9.44] | 7.25 |
| hierarchical, varying intercept | 6.54 | [4.63, 8.38] | 3.74 tightest — and misspecified |
| hierarchical, + varying slope | 6.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
- Efron, B. & Morris, C. (1975). Data Analysis Using Stein's Estimator and its Generalizations. Journal of the American Statistical Association 70(350), 311–319. — shrinkage, and why it beats the separate estimates
- Rubin, D. B. (1981). Estimation in Parallel Randomized Experiments. Journal of Educational Statistics 6(4), 377. — parallel randomized experiments
- Gelman, A. (2006). Conservative prior distributions for variance parameters in hierarchical models. Canadian Journal of Statistics 34(3), 377–390. — priors for the variance parameters