mlmoderator provides a unified workflow for probing,
plotting, and interpreting multilevel interaction effects from
mixed-effects models fitted with lme4::lmer().
The core workflow is:
mod <- lmer(math ~ ses * climate + (1 + ses | school), data = dat)
mlm_probe(mod, pred = "ses", modx = "climate") # simple slopes
mlm_jn(mod, pred = "ses", modx = "climate") # Johnson-Neyman region
mlm_plot(mod, pred = "ses", modx = "climate") # interaction plot
mlm_summary(mod, pred = "ses", modx = "climate") # full reportlibrary(mlmoderator)
library(lme4)
data(school_data)
head(school_data)
#> school student math ses climate gender
#> 1 1 1 61.06 -0.005 1.371 female
#> 2 1 2 60.54 0.760 1.371 male
#> 3 1 3 60.47 0.039 1.371 male
#> 4 1 4 59.82 0.735 1.371 female
#> 5 1 5 59.16 -0.146 1.371 male
#> 6 1 6 49.71 -0.058 1.371 female
str(school_data)
#> 'data.frame': 3000 obs. of 6 variables:
#> $ school : Factor w/ 100 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
#> $ student: int 1 2 3 4 5 6 7 8 9 10 ...
#> $ math : num 61.1 60.5 60.5 59.8 59.2 ...
#> $ ses : num -0.005 0.76 0.039 0.735 -0.146 ...
#> $ climate: num 1.37 1.37 1.37 1.37 1.37 ...
#> $ gender : Factor w/ 2 levels "female","male": 1 2 2 1 2 1 2 2 1 1 ...school_data contains 3,000 students nested in 100
schools. The outcome (math) was generated with a true
cross-level interaction between student SES (ses) and
school climate (climate).
Before fitting, it is good practice to center predictors.
mlm_center() supports grand-mean centering, group-mean
(within-cluster) centering, and a full within–between decomposition.
# Group-mean center SES within school
dat <- mlm_center(school_data,
vars = "ses",
cluster = "school",
type = "group")
# Grand-mean center school climate (level-2 variable)
dat <- mlm_center(dat,
vars = "climate",
type = "grand")
head(dat[, c("school", "ses", "ses_c", "climate", "climate_c")])
#> school ses ses_c climate climate_c
#> 1 1 -0.005 0.01513333 1.371 1.33844
#> 2 1 0.760 0.78013333 1.371 1.33844
#> 3 1 0.039 0.05913333 1.371 1.33844
#> 4 1 0.735 0.75513333 1.371 1.33844
#> 5 1 -0.146 -0.12586667 1.371 1.33844
#> 6 1 -0.058 -0.03786667 1.371 1.33844mod <- lmer(math ~ ses_c * climate_c + gender + (1 + ses_c | school),
data = dat, REML = TRUE)
summary(mod)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: math ~ ses_c * climate_c + gender + (1 + ses_c | school)
#> Data: dat
#>
#> REML criterion at convergence: 18414.9
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.7607 -0.6591 -0.0070 0.6580 3.1325
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> school (Intercept) 7.841 2.8002
#> ses_c 0.168 0.4099 -0.28
#> Residual 24.916 4.9916
#> Number of obs: 3000, groups: school, 100
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 49.84009 0.30871 161.446
#> ses_c 1.42623 0.10133 14.075
#> climate_c 1.01234 0.28421 3.562
#> gendermale 0.07486 0.18608 0.402
#> ses_c:climate_c 0.46554 0.09932 4.687
#>
#> Correlation of Fixed Effects:
#> (Intr) ses_c clmt_c gndrml
#> ses_c -0.105
#> climate_c -0.002 0.000
#> gendermale -0.300 0.012 0.007
#> ss_c:clmt_c 0.001 -0.009 -0.105 -0.005mlm_probe() computes the simple slope of the focal
predictor (ses_c) at selected values of the moderator
(climate_c).
probe_out <- mlm_probe(mod, pred = "ses_c", modx = "climate_c")
probe_out
#>
#> --- Simple Slopes: mlm_probe ---
#> Focal predictor : ses_c
#> Moderator : climate_c
#> Confidence level: 0.95
#>
#> climate_c slope se t df p ci_lower ci_upper
#> -1.036 0.944 0.145 6.504 2995 < .001 0.659 1.228
#> 0.000 1.426 0.101 14.075 2995 < .001 1.228 1.625
#> 1.036 1.909 0.144 13.277 2995 < .001 1.627 2.191By default, moderator values are set to −1 SD, Mean, and +1 SD. Other
options include "quartiles", "tertiles", or
"custom" values via at.
mlm_probe(mod, pred = "ses_c", modx = "climate_c",
modx.values = "quartiles")
#>
#> --- Simple Slopes: mlm_probe ---
#> Focal predictor : ses_c
#> Moderator : climate_c
#> Confidence level: 0.95
#>
#> climate_c slope se t df p ci_lower ci_upper
#> -0.649 1.124 0.121 9.318 2995 < .001 0.887 1.360
#> 0.057 1.453 0.101 14.324 2995 < .001 1.254 1.652
#> 0.629 1.719 0.119 14.502 2995 < .001 1.487 1.952The Johnson–Neyman (JN) interval identifies the exact moderator
value(s) at which the simple slope of ses_c crosses the
significance threshold.
jn_out <- mlm_jn(mod, pred = "ses_c", modx = "climate_c")
jn_out
#>
#> --- Johnson-Neyman Interval: mlm_jn ---
#> Focal predictor : ses_c
#> Moderator : climate_c
#> Alpha : 0.05
#> Moderator range : -3.026 to 2.254
#>
#> Johnson-Neyman boundary/boundaries ( climate_c ):
#> -2.088
#>
#> Interpretation: The simple slope of 'ses_c' is statistically
#> significant (p < 0.05) for values of 'climate_c'
#> between -2.07 and 2.254.mlm_plot() returns a ggplot object, so it
is fully customisable.
Add confidence bands, raw data points, or custom labels:
mlm_plot(mod,
pred = "ses_c",
modx = "climate_c",
points = TRUE,
point_alpha = 0.2,
x_label = "Student SES (group-mean centred)",
y_label = "Mathematics Achievement",
legend_title = "School Climate")
mlm_summary() combines the interaction coefficient,
simple slopes, and JN region into a single printable object.
mlm_summary(mod, pred = "ses_c", modx = "climate_c")
#>
#> ========================================
#> Multilevel Moderation Summary
#> ========================================
#> Focal predictor : ses_c
#> Moderator : climate_c
#> Confidence level: 0.95
#>
#> --- Interaction Term ---
#> ses_c:climate_c b = 0.466 SE = 0.099 t(2995) = 4.687 p = < .001 [0.271, 0.66]
#>
#> --- Simple Slopes ---
#> climate_c slope se t df p 95% CI lo 95% CI hi
#> ------------------------------------------------------------------------------------
#> -1.036 0.944 0.145 6.504 2995.0 < .001 0.659 1.228
#> 0.000 1.426 0.101 14.075 2995.0 < .001 1.228 1.625
#> 1.036 1.909 0.144 13.277 2995.0 < .001 1.627 2.191
#>
#> --- Johnson-Neyman Region ---
#> JN boundary/boundaries for 'climate_c': -2.088
#>
#> ========================================mlm_plot() output is a ggplot — you
can add theme_*(), labs(), or any other
ggplot2 layer on top.