library(tidyverse)
library(lme4)
library(immr)Group-Level Predictors
Abstract
In the examples considered so far, all predictor variables are at the observation level: each observation has its own value. Some predictors, however, take the same value for every observation in a group: they are defined at the group level. School size, sector, or mean socioeconomic status are properties of schools, not of individual pupils. This guide develops the use of group-level predictors in linear mixed effects models using the MathAchieve dataset, explains why group-level predictors cannot be assigned random slopes, and shows how they interact with observation-level predictors to produce cross-level interactions.
The data
Two related datasets from a study of mathematical achievement in US high schools.
mathach_df <- mathachieve
mathachschool_df <- mathachieveschool
mathach_df# A tibble: 7,185 × 6
pupil school minority sex ses mathach
<fct> <fct> <fct> <fct> <dbl> <dbl>
1 p1 s1224 No Female -1.53 5.88
2 p2 s1224 No Female -0.588 19.7
3 p3 s1224 No Male -0.528 20.3
4 p4 s1224 No Male -0.668 8.78
5 p5 s1224 No Male -0.158 17.9
6 p6 s1224 No Male 0.022 4.58
7 p7 s1224 No Female -0.618 -2.83
8 p8 s1224 No Male -0.998 0.523
9 p9 s1224 No Female -0.888 1.53
10 p10 s1224 No Male -0.458 21.5
# ℹ 7,175 more rowsmathachschool_df# A tibble: 160 × 7
school size sector pracad disclim himinty meanses
<fct> <dbl> <fct> <dbl> <dbl> <fct> <dbl>
1 s1224 842 Public 0.35 1.60 0 -0.428
2 s1288 1855 Public 0.27 0.174 0 0.128
3 s1296 1719 Public 0.32 -0.137 1 -0.42
4 s1308 716 Catholic 0.96 -0.622 0 0.534
5 s1317 455 Catholic 0.95 -1.69 1 0.351
6 s1358 1430 Public 0.25 1.54 0 -0.014
7 s1374 2400 Public 0.5 2.02 0 -0.007
8 s1433 899 Catholic 0.96 -0.321 0 0.718
9 s1436 185 Catholic 1 -1.14 0 0.569
10 s1461 1672 Public 0.78 2.10 0 0.683
# ℹ 150 more rowsThe mathach_df data records individual pupils with their mathematical achievement score (mathach), home socioeconomic status (ses), sex, and minority status. The mathachschool_df data records school-level properties: size, sector (public/private), proportion on an academic track (pracad), discrimination climate (disclim), proportion of minority students (himinty), and mean SES of the school (meanses).
Join the two datasets by school.
mathach_df2 <- left_join(mathach_df, mathachschool_df, by = "school")
mathach_df2# A tibble: 7,185 × 12
pupil school minority sex ses mathach size sector pracad disclim
<fct> <fct> <fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <dbl>
1 p1 s1224 No Female -1.53 5.88 842 Public 0.35 1.60
2 p2 s1224 No Female -0.588 19.7 842 Public 0.35 1.60
3 p3 s1224 No Male -0.528 20.3 842 Public 0.35 1.60
4 p4 s1224 No Male -0.668 8.78 842 Public 0.35 1.60
5 p5 s1224 No Male -0.158 17.9 842 Public 0.35 1.60
6 p6 s1224 No Male 0.022 4.58 842 Public 0.35 1.60
7 p7 s1224 No Female -0.618 -2.83 842 Public 0.35 1.60
8 p8 s1224 No Male -0.998 0.523 842 Public 0.35 1.60
9 p9 s1224 No Female -0.888 1.53 842 Public 0.35 1.60
10 p10 s1224 No Male -0.458 21.5 842 Public 0.35 1.60
# ℹ 7,175 more rows
# ℹ 2 more variables: himinty <fct>, meanses <dbl>A standard two-level model
Begin with a model in which ses (a pupil-level predictor) has both a fixed effect and a random slope across schools.
M_base <- lmer(mathach ~ ses + (ses | school), data = mathach_df2)
fixef(M_base)(Intercept) ses
12.66502 2.39381 VarCorr(M_base) Groups Name Std.Dev. Corr
school (Intercept) 2.1974
ses 0.6426 -0.109
Residual 6.0688 The fixed effect of ses gives the population-average slope. The random effects allow this slope to differ across schools, and the standard deviation of the random slopes quantifies that variation.
Why group-level predictors cannot be random slopes
Consider adding pracad — the proportion of students on an academic track — as a predictor. pracad is a school-level variable: every pupil in the same school has the same value.
Attempting to give pracad a random slope across schools fails.
lmer(mathach ~ pracad + (pracad | school), data = mathach_df2)boundary (singular) fit: see help('isSingular')Linear mixed model fit by REML ['lmerMod']
Formula: mathach ~ pracad + (pracad | school)
Data: mathach_df2
REML criterion at convergence: 47016.62
Random effects:
Groups Name Std.Dev. Corr
school (Intercept) 2.1966
pracad 0.2594 -1.00
Residual 6.2564
Number of obs: 7185, groups: school, 160
Fixed Effects:
(Intercept) pracad
8.391 8.219
optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings The model is not identified because the variation in pracad is entirely between schools, not within them. A random slope for pracad would require pracad to vary within schools, but it does not: each school has a single value.
Group-level predictor as a fixed effect
The correct way to include pracad is as a fixed effect, alongside a random intercept for schools.
M_glp <- lmer(mathach ~ pracad + (1 | school), data = mathach_df2)
fixef(M_glp)(Intercept) pracad
8.384325 8.232278 VarCorr(M_glp) Groups Name Std.Dev.
school (Intercept) 2.0614
Residual 6.2565 The fixed effect of pracad tells us how the average mathach in a school changes as that school’s academic track proportion increases. The random intercept captures residual between-school variation that pracad does not explain.
Cross-level interactions
A more informative question is whether the effect of ses on mathach depends on school-level characteristics. For example, does the ses effect differ between schools with a high proportion of minority students (himinty = 1) and those without?
This is a cross-level interaction: a group-level predictor interacting with an observation-level predictor.
Formally, the model allows the school-specific slope and intercept for ses to depend on himinty:
\[ \begin{aligned} y_i &\sim N(\mu_i, \sigma^2), \\ \mu_i &= \beta_{[s_i]0} + \beta_{[s_i]1} x_i, \\ \boldsymbol{\beta}_j &\sim N(\boldsymbol{b}_0 + \boldsymbol{b}_1 z_j, \Sigma), \end{aligned} \]
where \(x_i\) is pupil \(i\)’s ses and \(z_j\) is school \(j\)’s himinty. Expanding this shows that the fixed effects include an intercept, a ses term, a himinty1 term, and a ses:himinty1 interaction term.
M_glp2 <- lmer(
mathach ~ ses + himinty + ses:himinty + (ses | school),
data = mathach_df2
)
fixef(M_glp2) (Intercept) ses himinty1 ses:himinty1
13.1542443 2.5436697 -1.8593332 -0.5760956 The interaction coefficient ses:himinty1 tells us how the slope of ses differs between high-minority schools (himinty = 1) and low-minority schools (himinty = 0). If it is negative, the effect of socioeconomic background on mathematics achievement is smaller in schools with a high proportion of minority students.
VarCorr(M_glp2) Groups Name Std.Dev. Corr
school (Intercept) 2.05931
ses 0.60327 -0.285
Residual 6.06828 The residual variance in random slopes after accounting for himinty can be compared to that from M_base to assess how much of the school-to-school variation in the ses slope is explained by himinty.
Combining pupil-level and school-level predictors
There is no constraint on combining multiple observation-level and group-level predictors in the same model.
M_combined <- lmer(
mathach ~ ses + pracad + himinty + ses:himinty + (ses | school),
data = mathach_df2
)
fixef(M_combined) (Intercept) ses pracad himinty1 ses:himinty1
10.1198458 2.4317382 5.7305137 -1.6886875 -0.5445438 Each school-level predictor enters as a fixed effect. The random intercept for school captures residual between-school variation not explained by the included school-level predictors. The random slope for ses captures residual variation in the pupil-level effect across schools, beyond what himinty and other predictors account for.
Confidence intervals
The default confint uses profile likelihood across all parameters, including the variance components. For this model the profile for the school-level random effects correlation is non-monotonic near the boundary of the parameter space, which causes the profiling algorithm to fail with a cascade of warnings. Since the discussion here concerns only the fixed effects, we restrict to those using parm = "beta_".
confint(M_glp2, parm = "beta_")Computing profile confidence intervals ... 2.5 % 97.5 %
(Intercept) 12.739047 13.56895697
ses 2.269340 2.82193860
himinty1 -2.653663 -1.06338996
ses:himinty1 -1.071993 -0.07801305The confidence intervals for the fixed effects include the cross-level interaction. Wide intervals on the interaction term relative to its estimate suggest limited evidence for a differential ses effect between school types.