Three honest answers and one that is out on its own

The same coefficient, four standard errors

The statistics course established that the school feeding programme was assigned to 24 schools and measured on 1,200 children. Here is what each way of handling that produces.

import pandas as pd
import statsmodels.formula.api as smf

naive = smf.ols("rate ~ feeding_programme", data=per_student).fit()
robust = naive.get_robustcov_results(cov_type="cluster",
                                     groups=per_student["school_id"])

mixed = smf.mixedlm("rate ~ feeding_programme", data=per_student,
                    groups=per_student["school_id"]).fit()

by_school = (per_student.groupby(["school_id", "feeding_programme"])["rate"]
             .mean().reset_index())
aggregated = smf.ols("rate ~ feeding_programme", data=by_school).fit()

library(lme4); library(sandwich); library(lmtest)

naive      <- lm(rate ~ feeding_programme, data = per_student)
robust     <- coeftest(naive, vcov = vcovCL, cluster = ~school_id)
mixed      <- lmer(rate ~ feeding_programme + (1 | school_id), data = per_student)
aggregated <- lm(rate ~ feeding_programme, data = by_school)

Approach	Coefficient	SE	t	n
Naive OLS	+4.93 pts	0.88	5.63	1,200
Cluster-robust SE	+4.93 pts	1.46	3.38	1,200
Random intercept	+5.10 pts	1.48	3.44	1,200
School-level OLS	+5.21 pts	1.53	3.41	24

Three of the four agree and one does not. The coefficients span 0.3 points; the three honest standard errors span 0.07; and the naive one is 40% smaller than any of them.

That is the shape to expect. Clustering rarely changes an estimate much and routinely changes its precision a lot, and the three corrections are three routes to the same destination rather than three competing answers.

What a random intercept adds

Cluster-robust standard errors fix the inference and say nothing about the structure. A random intercept estimates it.

print(mixed.summary())
print(f"between-school variance: {float(mixed.cov_re.iloc[0, 0]):.5f}")
print(f"residual variance:       {mixed.scale:.5f}")

VarCorr(mixed)

Model	Between-school variance	Residual variance	ICC
Intercept only	0.00142	0.02042	0.065
+ feeding programme	0.00081	0.02042	0.038

Feeding explains 42.8% of the between-school variance and none of the within-school variance, which is exactly what a school-level programme should do and is worth checking because it is a way of catching a mis-specified model.

The intercept-only ICC of 0.065 is the same quantity the statistics course computed by hand from mean squares and got 0.060. Two estimators, two slightly different answers, one conclusion — about six per cent of the variation in attendance is between schools, which with fifty children per school is enough to quadruple the variance of a naive comparison.

Choosing among the three

The three corrections are not interchangeable, and the choice is usually decided by the number of clusters and by what varies at which level.

Use	When	Cost
Aggregate to the cluster	Few clusters; the exposure is cluster-level	A big school counts the same as a small one
Cluster-robust SE	Many clusters (30+); you want the individual-level model	Unreliable below about 30 clusters
Random intercept	You want the variance components, or predictors at both levels	Assumes the random effect is uncorrelated with the predictors

With 24 clusters, aggregate or fit the random intercept and say which. The cluster-robust sandwich is the standard tool and it is the one that behaves worst here: its asymptotics need more groups than this design has, and it will run without warning you.

When the exposure varies only between clusters, all three converge, which the table above shows. They come apart when a predictor varies within clusters — and then the random intercept is doing work the other two cannot.

Predictors at two levels

This is where a multilevel model stops being a correction and starts being a model.

two_level = smf.mixedlm(
    "rate ~ feeding_programme + age_years + disability_reported",
    data=d, groups=d["school_id"]).fit()
print(two_level.summary().tables[1])

lmer(rate ~ feeding_programme + age_years + disability_reported +
       (1 | school_id), data = d)

feeding_programme varies only between schools. age_years and disability_reported vary within them. A single model can carry both, and the random intercept is what lets it: the within-school predictors are estimated from within-school variation and the between-school predictor from between-school variation, without either contaminating the other’s standard error.

The failure mode is putting a cluster-level exposure in a model with no cluster term. That is the naive row of the first table, and it is the default that every software package gives you.

Three levels, and when to stop

Children sit in households, households in villages, villages in districts.

# Two grouping levels: households nested within enumeration areas.
smf.mixedlm("outcome ~ x", data=survey,
            groups=survey["ea_id"], re_formula="1").fit()

lmer(outcome ~ x + (1 | ea_id / household_id), data = survey)

Add a level when something is assigned or measured at that level and you have enough units of it. Three levels with six districts at the top will not estimate a district variance worth reporting — the top level needs enough groups for the same reason the multiple-comparisons lesson needed enough tests.

Do not add a level for tidiness. A level with four groups adds a parameter that cannot be estimated and a false sense of having handled something.

Report it whole

Attendance and school feeding, multilevel analysis

  Linear mixed model, 1,200 students in 24 schools.
  rate ~ feeding_programme + (1 | school_id)

  Feeding programme    +5.10 points   95% CI +2.20 to +8.00

  Between-school SD    0.028      Residual SD  0.143
  ICC 0.038 with the programme term; 0.065 without it, so the programme
  accounts for 43% of the between-school variance.

  The naive single-level model gives +4.93 points with a standard error of
  0.88 rather than 1.48. It is not reported: it treats 50 children in one
  school as 50 independent observations.

  Cluster-robust standard errors on the single-level model (+4.93, SE 1.46)
  and school-level aggregation (+5.21, SE 1.53) give the same answer. With
  24 clusters the robust sandwich is at the edge of its assumptions and is
  reported as a check rather than as the headline.

  Observational. Schools were not randomised into the programme.

Reporting all three is four lines and it forecloses the obvious challenge — that the result depends on the method. It does not, and showing that is cheaper than arguing it.

What comes next

This lesson matched the model to how the programme was assigned. The next one matches it to how the sample was drawn, on a survey where the strata were deliberately unequal — and finds coefficients that cross the significance line in both directions when the weights go in.