School effects

This example fits a linear mixed model to the 250 students from 12 schools in the school scores data. Test score depends on study hours, and each school has its own random intercept. The model estimates how much schools differ and how much of each school's apparent effect is noise.

The page compares the mixed model with a no-pooling model that gives every school a separate fixed intercept. The mixed model pulls each school's intercept toward the overall intercept, and it pulls schools with fewer students further.

Code

use function Phit\fit_lmer;

$scores = DataFrame::fromQuery(StudentScore::query());

$mixed = fit_lmer('score ~ study_hours + (1 | school)', $scores, reml: true);
$noPooling = $scores |> lm('score ~ 0 + school + study_hours');

$mixed['random'][0]['covariance'];  // school variance
$mixed['random'][0]['effects'];     // one random intercept per school

use function Phit\fit_lmer;

$scores = DataFrame::fromQuery(StudentScore::query());

$mixed = fit_lmer('score ~ study_hours + (1 | school)', $scores, reml: true);
$noPooling = phplyr($scores)->lm('score ~ 0 + school + study_hours');

$mixed['random'][0]['covariance'];  // school variance
$mixed['random'][0]['effects'];     // one random intercept per school

Phit\fit_lmer() takes the phibble DataFrame directly. Its results are verified against R's lme4.

Estimation method REML Maximum likelihood

Fixed effects

Term	Estimate	Std. error	z	95% CI
(Intercept)	56.493	1.541	36.66	[53.473, 59.513]
study_hours	2.343	0.109	21.52	[2.129, 2.556]

Variance components

Group	Term	Variance	Std. dev.
school	(Intercept)	20.867	4.568
Residual		29.094	5.394

REML fit on 250 students in 12 schools. Log-likelihood = -791.66, AIC = 1,591.3. Schools account for 41.8% of the variance in scores after adjusting for study hours (the intraclass correlation).

Shrinkage of school intercepts

Red points are the no-pooling intercepts. Blue points are the mixed-model intercepts, which move toward the overall intercept (dashed line).

School	Students	No pooling	Partial pooling	Shrinkage
school-01	18	57.48	57.39	7.2%
school-02	24	58.41	58.28	5.5%
school-03	15	53.74	53.95	8.5%
school-04	22	50.43	50.76	6.0%
school-05	20	57.14	57.07	6.5%
school-06	27	59.51	59.34	4.9%
school-07	19	67.25	66.49	6.8%
school-08	21	50.65	50.99	6.2%
school-09	16	59.17	58.93	8.0%
school-10	25	56.76	56.73	5.3%
school-11	23	56.70	56.67	5.7%
school-12	20	50.98	51.32	6.5%

Shrinkage is σ²/n ÷ (τ² + σ²/n), where τ² is the school variance, σ² is the residual variance, and n is the number of students in the school. Each school's random intercept is its mean residual from the fixed part of the model, reduced by this share. Schools with fewer students are shrunk more. Shrinkage is small in these data because the school variance is large compared with σ²/n.

The no-pooling model fits its own study-hours slope, so the distance between its intercepts and the mixed-model intercepts is not exactly this share.

What if

School Study hours

Predicted score at a new school: 75.2
Population prediction: 75.2

Phit\predict() gives both predictions for a student who studies 8 hours. The first adds the school's random intercept, as lme4's predict() does by default. The second uses the fixed effects only, as predict() does with re.form = NA. A school outside the data has no estimated random intercept, so with allowNewLevels it gets the population prediction. lme4 does not compute standard errors for these predictions, and neither does Phit\predict().