Nested data needs a multilevel model

When your data come in groups (students inside classrooms, patients inside hospitals, repeated measurements inside the same person), ordinary regression quietly assumes something false: that every row is independent. This post walks through why that matters, how a multilevel (mixed-effects) model fixes it, and how to fit one in R with lme4.

The problem: rows that know each other

Imagine 3,000 students spread across 100 classrooms. Students in the same classroom share a teacher, a timetable, and often a neighbourhood. Their test scores will be more alike than scores drawn at random from the whole population.

Ordinary least squares (OLS) does not know that. It treats 3,000 rows as 3,000 independent pieces of information. In reality, part of the information is duplicated within each classroom, so OLS standard errors come out too small, confidence intervals too narrow, and p-values too optimistic. Classroom-level predictors such as "teacher has a master's degree" are hit hardest, because there are really only 100 independent units for that variable, not 3,000.

Measuring the similarity: the intraclass correlation

The standard way to quantify "how alike are members of the same group" is the intraclass correlation (ICC). Start with an empty (intercept-only) multilevel model:

y_ij = γ00 + u_j + e_ij

• u_j is classroom j's deviation from the grand mean, with variance τ²
• e_ij is student i's deviation from their classroom mean, with variance σ²

Then

ICC = τ² / (τ² + σ²)

It is the share of total variance that lies between groups, and it also equals the expected correlation between two observations from the same group.

From ICC to the design effect

Survey statisticians have long used the design effect to translate the ICC into lost precision. Following Kish (1965), for groups of average size m:

deff ≈ 1 + (m − 1) × ICC

A purely hypothetical illustration: 100 classrooms of 30 students and an ICC of 0.2 give deff = 1 + 29 × 0.2 = 6.8. The effective sample size is about 3000 / 6.8 ≈ 441, and naive standard errors would be roughly 1/√6.8 ≈ 0.38 of what they should be. Even a modest ICC of 0.05 gives a design effect of 2.45 with classes of 30. Small ICCs matter when groups are large.

Three ways to handle groups

Gelman and Hill (2007) frame the choice as three options:

1. Complete pooling. Ignore groups and fit one regression. Simple, but it erases between-group variation (which is often the very thing you care about) and gives overconfident standard errors.
2. No pooling. Fit each group separately, or add a dummy variable per group. Small groups get noisy estimates, and group dummies are perfectly collinear with any group-level predictor, so you cannot estimate the effect of, say, school type at the same time.
3. Partial pooling (multilevel model). Each group gets its own parameters, but those parameters are drawn from a common distribution. Groups with little data are "shrunk" towards the overall mean; groups with lots of data mostly keep their own estimates.

Partial pooling is the core idea. It is a data-driven compromise between "all groups are the same" and "all groups are unrelated".

Random intercepts and random slopes

A random-intercept model lets each group have its own baseline:

y_ij = γ00 + γ10 x_ij + u0j + e_ij

A random-slope model also lets the effect of a predictor vary by group:

y_ij = γ00 + (γ10 + u1j) x_ij + u0j + e_ij

For example, the link between hours of homework and test scores might be steeper in some classrooms than in others. The variance of u1j tells you how much that slope varies, and the covariance between u0j and u1j tells you whether high-baseline groups also have steeper slopes.

Fitting it in R with lme4

The lme4 package (Bates, Mächler, Bolker & Walker, 2015) is the workhorse for this in R. It ships with the sleepstudy dataset: reaction times for 18 participants measured over 10 days of sleep restriction, so measurements are nested in people.

library(lme4)

# Complete pooling: one line for everyone
m_ols <- lm(Reaction ~ Days, data = sleepstudy)

# Random intercept per participant
m_ri <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

# Random intercept and random slope for Days
m_rs <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

summary(m_rs)

Things worth comparing yourself:

• The standard error of Days in m_ols versus m_rs. The OLS version treats 180 rows as independent.
• The Subject variance components in summary(m_ri). Plug them into the ICC formula above.
• coef(m_rs)$Subject against per-subject lm() fits. You will see shrinkage: participants with extreme individual slopes are pulled towards the average.

lmer() uses REML (restricted maximum likelihood) by default, which gives less biased variance estimates than plain ML. If you compare models that differ in their fixed effects with a likelihood-ratio test, refit with REML = FALSE first. anova() in lme4 does this automatically for merMod objects.

Common mistakes

Too few groups. Variance components are estimated from the number of groups, not the number of rows. Maas and Hox (2005) found that with very few groups the standard errors of the group-level variances were noticeably biased. Five schools with 1,000 students each is a weak basis for a random effect.

Ignoring the design and just clustering standard errors. Cluster-robust standard errors can fix inference for fixed effects, which is a reasonable choice in some settings. But they do not give you estimates of between-group variation, group-specific predictions, or random slopes. Pick the tool that answers your question.

Adding random slopes for everything. Models with many random slopes and correlations often fail to converge or produce singular fits (isSingular() returns TRUE). Start from a theoretically motivated structure and simplify when the data cannot support it.

Forgetting centring. With a random slope on x, the intercept variance depends on where x = 0. Centring predictors, either at the grand mean or within each group, changes what the intercept means. Within-group centring also separates within-group and between-group effects, which are often different things entirely.

Treating the ICC as a test of whether you "need" a multilevel model. A small ICC can still produce a large design effect when groups are big, and the grouping is part of the design whatever the ICC turns out to be.

When to reach for a multilevel model

Use one when:

• observations are grouped, and the grouping comes from how the data were collected (sampling, randomisation, or repeated measurement);
• you care about variation between groups, or about predictors measured at the group level;
• some groups are small and you want stable estimates for them.

If you want a longer, education-focused walk-through of the same ideas, including a worked example with classroom data, the original article on nested data and multilevel models goes into more detail.

References

• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48.
• Gelman, A., & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
• Kish, L. (1965). Survey Sampling. Wiley.
• Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 86–92.
• Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical Linear Models: Applications and Data Analysis Methods (2nd ed.). Sage.