Hierarchical Models and Partial Pooling

Eight high schools ran randomised trials of SAT coaching. School A estimated an effect of 28 points with a standard error of 15; school C estimated \(-3\) with a standard error of 16 (Rubin, 1981, Estimation in Parallel Randomized Experiments, Journal of Educational Statistics 6(4)). Taken separately, school A's coaching looks remarkable. Pooled into one precision-weighted average, every school gets \(7.7 \pm 4.1\). Neither answer is right. The first ignores that the eight programmes were similar; the second pretends they were identical. The spread of the eight estimates is in fact smaller than their standard errors alone would produce, and a hierarchical model is the machinery for turning that observation into an estimate.

The model and the shrinkage it implies

Let \(y_j\) be group \(j\)'s observed estimate with known standard error \(\sigma_j\), and \(\theta_j\) its true effect:

\(\mu\) is the population mean effect and \(\tau\) the between-group standard deviation. The groups are treated as exchangeable: before seeing data, nothing distinguishes school A from school C.

Conditional on \(\mu\) and \(\tau\), the posterior for each group is a precision-weighted compromise:

\(B_j\) is the shrinkage factor. Noisy groups (large \(\sigma_j\)) shrink hard toward the mean; precise groups keep their own estimate. \(\tau\) interpolates between the extremes: \(\tau \to 0\) gives complete pooling, \(\tau \to \infty\) gives no pooling.

Put numbers on school A. If \(\tau = 10\), then \(B_A = 225/325 = 0.69\) and the estimate moves from 28 to about 14. If \(\tau = 5\), \(B_A = 0.9\) and it moves to about 10. The whole question becomes how large \(\tau\) is, and a full Bayesian analysis does not pick one value: it averages over the posterior of \(\tau\), so the uncertainty about how much to pool is carried into every group's interval.

Why not estimate tau and plug it in

Empirical Bayes estimates \(\tau\) from the data and conditions on it. For eight schools the method-of-moments estimate is exactly zero, because the observed spread is below what sampling error predicts: the weighted sum of squared deviations from the pooled mean is 4.7, against an expected value of 7 under \(\tau = 0\). Plugging in \(\hat{\tau} = 0\) collapses every school to 7.7 with intervals that ignore any possibility of real variation. The posterior for \(\tau\) instead peaks near zero but keeps non-negligible mass at values of 10 or more, and the resulting intervals for school A stay wide. This is the standard argument for full Bayes over plug-in estimates when the number of groups is small, and it is also where the prior on \(\tau\) starts to matter.

That prior was itself contested. The once-common inverse-gamma\((\epsilon, \epsilon)\) prior on \(\tau^2\) is not uninformative: when the data are compatible with \(\tau\) near zero, inferences depend strongly on \(\epsilon\). Gelman recommended uniform or half-\(t\) (including half-Cauchy) priors on \(\tau\) itself instead (Gelman, 2006, Prior distributions for variance parameters in hierarchical models, Bayesian Analysis 1(3)), and weakly informative half-normal priors are now a common default. The priors-and-conjugacy concept on this track explains why "uninformative" is always relative to a scale.

Centred versus non-centred parameterisation

The model as written is centred: \(\theta_j\) is sampled directly from \(\mathcal{N}(\mu, \tau^2)\). When \(\tau\) is small, every \(\theta_j\) is forced into a narrow band around \(\mu\); when \(\tau\) is large they range freely. The joint posterior of \((\theta, \log\tau)\) is a funnel whose neck is orders of magnitude narrower than its mouth, and Hamiltonian Monte Carlo with one step size cannot explore both. The sampler reports divergent transitions and, worse, under-samples small \(\tau\), biasing the answer (Betancourt and Girolami, 2013, Hamiltonian Monte Carlo for Hierarchical Models, arXiv:1312.0906).

The non-centred version samples a standardised offset and rescales it:

Now \(\eta\) and \(\tau\) are a priori independent, and the funnel disappears from the space the sampler moves in. It is the same model with a different geometry.

Non-centring is not universally better, which is the part often missed. When the data strongly inform each group (small \(\sigma_j\) relative to \(\tau\), or many observations per group), the likelihood pins each \(\theta_j\) and the non-centred form creates strong posterior dependence between \(\eta_j\) and \(\tau\) instead. The two parameterisations are complementary: centred wins when data dominate, non-centred when the prior dominates (Papaspiliopoulos, Roberts and Sköld, 2007, A General Framework for the Parametrization of Hierarchical Models, Statistical Science 22(1)). Models with some data-rich and some data-poor groups can mix the two per group.

When it breaks

Few groups means the prior on tau is doing real work. With eight groups there are effectively eight data points about \(\tau\). Report sensitivity to the \(\tau\) prior rather than presenting one posterior as the answer.

Exchangeability is an assumption about the grouping. If schools differ systematically, say in baseline scores that predict coaching effect, shrinking all toward one mean biases the distinctive ones. The fix is to model \(\theta_j \sim \mathcal{N}(X_j\beta, \tau^2)\) with group-level predictors, so groups shrink toward their own conditional mean.

Shrinkage punishes genuine outliers. A group that truly is different, with little data, will be pulled toward the population. That is the right behaviour on average and the wrong one for that group, and a decision about that specific group deserves its own scrutiny.

Normal group effects hide heavy tails. If a few groups have very large effects, a normal population distribution over-shrinks them; a Student-\(t\) group distribution is the usual robustification, at the cost of harder sampling.