A complete workflow with real data: SES and school context in High School and Beyond

This vignette analyses a classic cross-level interaction with public data: does the within-school relationship between students’ socioeconomic status (SES) and mathematics achievement depend on the average SES of the school? The data are the 1982 High School and Beyond subsample analysed by Raudenbush and Bryk (2002): 7,185 students in 160 schools.

Model

cses is student SES centred at the school mean, so its coefficient is a purely within-school slope. meanses is the school mean SES, a cluster-level moderator. Following Raudenbush and Bryk, the model also lets the SES slope differ by school sector.

library(mlmoderator)
library(lme4)
data("Hsb82", package = "mlmRev")

fit <- lmer(mAch ~ cses * meanses + cses * sector + (1 + cses | school),
            data = Hsb82)

Probing the interaction

mlm_summary(fit, pred = "cses", modx = "meanses")
#> 
#> ========================================
#>   Multilevel Moderation Summary
#> ========================================
#> Focal predictor : cses 
#> Moderator       : meanses 
#> Confidence level: 0.95 
#> df method       : satterthwaite 
#> 
#> --- Interaction Term ---
#>   cses:meanses                  b =   1.039  SE =  0.299  t(160.4) =  3.477  p = < .001  [0.449, 1.629]
#> 
#> --- Simple Slopes ---
#>   meanses        slope        se         t        df         p  95% CI lo  95% CI hi 
#> ------------------------------------------------------------------------------------ 
#>   -0.420         2.508     0.169    14.847     131.5    < .001     2.174     2.843
#>   -0.006         2.939     0.155    18.951     139.2    < .001     2.632     3.245
#>   0.408          3.369     0.224    15.042     156.9    < .001     2.926     3.811
#> 
#> --- Johnson-Neyman Region ---
#>   Slope of 'cses' is significant across the entire range of 'meanses'.
#> 
#> ========================================

Because meanses is constant within schools, the default probe points (mean and ±1 SD) are computed across the 160 schools rather than across students, so large schools do not dominate them.

Tests use Satterthwaite degrees of freedom by default. Inference for a cross-level interaction draws its information from the clusters, so the relevant degrees of freedom are on the order of the number of schools, not the number of students:

sapply(c("satterthwaite", "kenward-roger", "between", "residual"),
       function(m) {
         r <- mlm_summary(fit, "cses", "meanses", jn = FALSE,
                          df_method = m)$interaction
         round(c(SE = r$se, df = r$df, p = r$p), 5)
       })
#>    satterthwaite kenward-roger   between   residual
#> SE       0.29888       0.29945   0.29888    0.29888
#> df     160.43742     159.06791 157.00000 7179.00000
#> p        0.00065       0.00067   0.00066    0.00051

With 160 schools the choice hardly matters here. It matters a great deal with few clusters: in a simulation with 10 clusters and no true interaction, the "residual" rule (students minus fixed effects, the default before mlmoderator 0.3.0) rejected in 10.0% of 400 replications at the nominal 5% level, while Satterthwaite rejected in 5.5% and the between-cluster rule in 5.75%.

Plots

mlm_plot(fit, pred = "cses", modx = "meanses",
         x_label = "Student SES (school-centred)",
         y_label = "Mathematics achievement",
         legend_title = "School mean SES")

plot(mlm_jn(fit, pred = "cses", modx = "meanses"))

How much do school slopes vary beyond the moderator?

The interaction describes how the average SES slope changes with school SES. Individual schools still differ around that average:

vd <- mlm_variance_decomp(fit, pred = "cses", modx = "meanses")
vd
#> 
#> ========================================
#>   Slope Variance Decomposition
#> ========================================
#> Focal predictor : cses 
#> Moderator       : meanses 
#> Cluster grouping: school 
#> Confidence level: 0.95 
#> 
#> Random-slope SD (tau11):  0.318
#> Random-slope Var (tau11): 0.101
#> 
#> --- Per-moderator-value intervals ---
#> 
#>  meanses slope    SE    df 95% CI (average slope) 95% PI (new cluster)
#>   -0.420 2.508 0.169 131.5         [2.174, 2.843]       [1.796, 3.221]
#>   -0.006 2.939 0.155 139.2         [2.632, 3.245]       [2.239, 3.638]
#>    0.408 3.369 0.224 156.9         [2.926, 3.811]       [2.601, 4.137]
#> 
#> The prediction interval treats tau11 as known; it is approximate
#> with few clusters.
plot(vd)

The confidence intervals describe the average slope at each value of meanses; the prediction intervals describe the slope to expect in a new school with that mean SES. They are wider because school mean SES explains only part of the between-school variation in SES slopes.

Is the interaction driven by a few schools?

sens <- mlm_sensitivity(fit, pred = "cses", modx = "meanses",
                        df_method = "between")
sens
#> 
#> --- Leave-one-cluster-out influence: mlm_sensitivity ---
#> Interaction     : cses:meanses 
#> Clusters        : 160 ( school )
#> df method       : between 
#> 
#> Full data       : b = 1.039, SE = 0.299, df = 157.0, p = < .001
#> Refitted range  : 0.898 to 1.120
#> Same sign and significance decision in 100.0% of refits
#> Flagged (|DFBETA| > 0.158): 18 cluster(s)
#>   2277         DFBETA = -0.472  b without it = 0.898
#>   2655         DFBETA =  0.270  b without it = 1.120
#>   1461         DFBETA = -0.267  b without it = 0.960
#>   2639         DFBETA = -0.261  b without it = 0.961
#>   8800         DFBETA =  0.251  b without it = 1.114
plot(sens)

Dropping any single school leaves the interaction positive and significant. A few schools move the estimate by more than the screening cutoff and are worth inspecting, but none changes the conclusion.