When One Model Is Not Enough: Hierarchical Models for Grouped Data

A dataset with 600 observations and 30 groups. The question is simple: does a student's score depend on a predictor, and does that relationship vary by school? Answering it with a single regression line throws away information. Answering it with 30 separate regressions throws away statistical power. Hierarchical models do neither.

Three models were fit to the same synthetic dataset using PyMC5 with NUTS. Here is what changed.

The Models

Three approaches, same data, same NUTS sampler, same 4 chains with 2,000 draws each. What varied was how group membership entered the model.

WAIC is the Widely Applicable Information Criterion. Lower values indicate better out-of-sample predictive fit. The full hierarchical model scores 396 WAIC points lower than the pooled model - a 19% improvement. That gap quantifies how much information group structure carries.

Why Not 30 Separate Regressions?

The obvious alternative: fit an independent regression for each of the 30 schools. This avoids the pooled model's naive assumption that all groups are identical, but it loses the ability to share information between groups.

The problem is sample size. Each school has only 20 observations. A separate regression for each school estimates its own slope and intercept from those 20 points. The hierarchical model does something different: it estimates school-level parameters while pooling information across schools. Schools with noisy data get pulled toward the group average. Schools with strong signals keep their individual estimates. This is partial pooling - it sits between the extremes of "one line for everyone" and "everyone for themselves."

The Intraclass Correlation Coefficient was 0.73, meaning 73% of the variance in the outcome comes from school-level differences. Grouping is not optional here - it is the structure of the data.

How The Model Learned

NUTS converged without issues. All R-hat values were at or below 1.01, and effective sample sizes ranged from 405 to 7,357. Zero divergent transitions across 8,000 posterior draws.

Recovery checks - comparing posterior estimates to the known true values used to generate the synthetic data - showed what the hierarchy accomplished:

• Intercept recovery: r = 0.99. The model recovered school-level baselines nearly perfectly. The hierarchical prior borrowed strength from the group mean without washing out individual differences.

• Slope recovery: r = 0.86. Slopes are harder to estimate from 20 observations each, and the model showed appropriate caution, shrinking uncertain school slopes toward the group mean. This is not a bug - it is the model correctly recognizing that some schools do not have enough data to support a confident individual slope estimate.

Posterior predictive checks confirmed the fit. The posterior predictive mean was 68.89, compared to an observed mean of 68.90 - a difference of 0.01.

Where Hierarchical Models Apply

Any dataset with repeated observations inside groups. The groups do not need to be schools.

Detection models measured across browser varieties. Treatment effects tracked across 20 clinics with 50 patients each. A/B tests where users are nested inside geographic regions. Server performance metrics grouped by datacenter. Survey responses clustered by interviewer. Any situation where observations share membership in a higher-level unit and that membership likely affects the outcome.

The test for whether a hierarchical structure is needed: run the pooled model, extract group-level residuals, and compute the Intraclass Correlation Coefficient. If ICC is above 0.1, the groups carry information. If ICC is above 0.4, flat regression is actively misleading.

The implementation uses PyMC5, a probabilistic programming library that builds on PyMC and PyTensor. Model specification follows the standard multilevel formula notation: a population-level slope and intercept, plus group-level offsets drawn from a shared prior distribution. NUTS handles the rest.

No special hardware required. The three models trained in sequence on a commodity CPU. The varying-intercept model with 34 parameters converged in under a minute. The full hierarchical model with 64 parameters converged in under two.