Getting Started with mlmoderator

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Overview

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 report

The built-in dataset

library(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).

Step 1: Center variables

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.33844

Step 2: Fit the model

mod <- 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.005

Step 3: Probe simple slopes

mlm_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.191

By 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.952

Step 4: Johnson–Neyman interval

The 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.

Step 5: Plot the interaction

mlm_plot() returns a ggplot object, so it is fully customisable.

mlm_plot(mod, pred = "ses_c", modx = "climate_c")

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")

Step 6: Full summary report

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
#> 
#> ========================================

Tips

Always center your predictors before fitting the model. This makes the simple slopes interpretable and reduces multicollinearity.
For cross-level interactions, group-mean-center the level-1 predictor and grand-mean-center the level-2 moderator.
The mlm_plot() output is a ggplot — you can add theme_*(), labs(), or any other ggplot2 layer on top.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.