Package {mlmoderator}

Type:	Package
Title:	Probing, Plotting, and Interpreting Multilevel Interaction Effects
Version:	0.3.0
Description:	Provides a workflow for probing, plotting, and checking cross-level interaction effects in two-level mixed-effects models fitted with 'lme4' (Bates et al., 2015) <doi:10.18637/jss.v067.i01>. Implements simple slopes analysis following Aiken and West (1991, ISBN:9780761907121), Johnson-Neyman intervals following Johnson and Fay (1950) <doi:10.1007/BF02288864> and Bauer and Curran (2005) <doi:10.1207/s15327906mbr4003_5>, and grand- or group-mean centering as described in Enders and Tofighi (2007) <doi:10.1037/1082-989X.12.2.121>. Tests and intervals use Satterthwaite degrees of freedom via 'lmerTest' (Kuznetsova et al., 2017) <doi:10.18637/jss.v082.i13> by default, with Kenward-Roger and between-cluster alternatives. Also provides confidence and new-cluster prediction intervals for simple slopes in random-slope models, contour plots of predicted outcomes over the predictor-by-moderator space, and leave-one-cluster-out influence diagnostics for the interaction.
License:	MIT + file LICENSE
Encoding:	UTF-8
Language:	en-US
LazyData:	true
Depends:	R (≥ 4.1.0)
Imports:	lme4, lmerTest (≥ 3.1-0), ggplot2, stats, utils, grDevices, rlang
Suggests:	testthat (≥ 3.0.0), dplyr, tibble, knitr, rmarkdown, patchwork, pbkrtest, emmeans, mlmRev
Config/testthat/edition:	3
VignetteBuilder:	knitr
URL:	https://github.com/subirhait/mlmoderator
BugReports:	https://github.com/subirhait/mlmoderator/issues
NeedsCompilation:	no
Config/roxygen2/version:	8.0.0
Packaged:	2026-09-29 23:11:00 UTC; subir
Author:	Subir Hait [aut, cre]
Maintainer:	Subir Hait <haitsubi@msu.edu>
Repository:	CRAN
Date/Publication:	2026-09-30 15:00:39 UTC

mlmoderator: Probing and Visualizing Multilevel Interaction Effects

Description

'mlmoderator' provides a unified workflow for estimating, probing, and visualizing multilevel moderation effects from mixed-effects models fitted with 'lme4::lmer()'. It supports simple slopes analysis, Johnson—Neyman intervals, publication-ready interaction plots, and grand- or group-mean centering. Tests and intervals use Satterthwaite degrees of freedom by default (see the 'df_method' argument of [mlm_probe()]).

## Core functions

| Function | Description | |—|—| | 'mlm_center()' | Grand- or group-mean center variables | | 'mlm_probe()' | Compute simple slopes at selected moderator values | | 'mlm_jn()' | Johnson-Neyman significance regions | | 'mlm_plot()' | Publication-ready interaction plot | | 'mlm_summary()' | Consolidated moderation report | | 'mlm_variance_decomp()' | Confidence vs. new-cluster prediction intervals for slopes | | 'mlm_surface()' | Slope surface heatmap over the full predictor x moderator space | | 'mlm_sensitivity()' | Leave-one-cluster-out influence on the interaction |

## Typical workflow

“'r library(mlmoderator) library(lme4)

data(school_data)

# 1. Center variables dat <- mlm_center(school_data, vars = "ses", cluster = "school", type = "group")

# 2. Fit model mod <- lmer(math ~ ses * climate + gender + (1 + ses | school), data = dat)

# 3. Probe interaction mlm_probe(mod, pred = "ses", modx = "climate")

# 4. Johnson-Neyman interval mlm_jn(mod, pred = "ses", modx = "climate")

# 5. Plot mlm_plot(mod, pred = "ses", modx = "climate")

# 6. Summary report mlm_summary(mod, pred = "ses", modx = "climate") “'

Author(s)

Maintainer: Subir Hait haitsubi@msu.edu (ORCID)

Authors:

• Subir Hait haitsubi@msu.edu (ORCID)

See Also

Useful links:

• https://github.com/subirhait/mlmoderator

• Report bugs at https://github.com/subirhait/mlmoderator/issues

Center variables for multilevel modeling

Description

Performs grand-mean centering, group-mean centering, or both (within-between decomposition) on one or more variables in a data frame. Group-mean centering is the standard preparation for cross-level interaction models.

Usage

mlm_center(
  data,
  vars,
  cluster = NULL,
  type = c("grand", "group", "both"),
  suffix_within = "_within",
  suffix_between = "_between"
)

Arguments

Value

The input data frame with new centered columns appended. Original columns are not modified.

Examples

data(school_data)

# Grand-mean center SES
d1 <- mlm_center(school_data, vars = "ses", type = "grand")
head(d1[, c("ses", "ses_c")])

# Group-mean center SES within schools
d2 <- mlm_center(school_data, vars = "ses", cluster = "school", type = "group")
head(d2[, c("ses", "ses_c")])

# Within-between decomposition
d3 <- mlm_center(school_data, vars = "ses", cluster = "school", type = "both")
head(d3[, c("ses", "ses_within", "ses_between")])

Johnson—Neyman interval for multilevel two-way interactions

Description

Computes the Johnson—Neyman (JN) interval: the region(s) of the moderator ('modx') where the simple slope of 'pred' transitions between statistical significance and non-significance. Useful for identifying *exactly* at which moderator values an effect becomes significant.

Usage

mlm_jn(
  model,
  pred,
  modx,
  alpha = 0.05,
  modx.range = NULL,
  grid = 200L,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual")
)

Arguments

Details

A Johnson-Neyman boundary is a moderator value w at which |t(w)| = t_{1-\alpha/2,\,df(w)}, where t(w) = (\hat\beta_1 + \hat\beta_3 w) / SE(w). With a constant number of degrees of freedom ('df_method = "between"' or '"residual"') the boundaries are the real roots of a quadratic in w (Bauer and Curran, 2005), and these closed-form roots are returned. With Satterthwaite or Kenward-Roger degrees of freedom, which change with w, the boundaries are found by root-finding ('stats::uniroot()') on the exact criterion, starting from brackets located on the grid.

Only boundaries inside 'modx.range' are returned in 'jn_bounds'. The closed-form roots outside that range are reported in 'jn_bounds_all' when available, because a boundary outside the observed data is an extrapolation.

Value

An object of class 'mlm_jn' with components: * 'jn_bounds': numeric vector of moderator values inside 'modx.range' where the slope crosses the significance threshold. 'NA' if none. * 'jn_bounds_all': for constant-df methods, all real roots of the Johnson-Neyman quadratic (possibly outside the data); 'NULL' otherwise. * 'slopes_df': data frame of slope estimates and significance across the grid. * 'pred', 'modx', 'alpha', 'modx.range', 'df_method'.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
jn <- mlm_jn(mod, pred = "x", modx = "m")
print(jn)

Publication-ready interaction plot for multilevel models

Description

Creates a 'ggplot2'-based interaction plot showing predicted values of the outcome across levels of 'pred', with separate lines for each selected value of 'modx'. Confidence bands and raw data overlay are optional.

Usage

mlm_plot(
  model,
  pred,
  modx,
  modx.values = c("mean-sd", "quartiles", "tertiles", "custom"),
  at = NULL,
  interval = TRUE,
  conf.level = 0.95,
  points = FALSE,
  point_alpha = 0.3,
  colors = NULL,
  line_size = 1,
  x_label = NULL,
  y_label = NULL,
  legend_title = NULL,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual"),
  modx.level = c("auto", "cluster", "observation")
)

Arguments

Value

A 'ggplot' object.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200),
  x   = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
mlm_plot(mod, pred = "x", modx = "m")
mlm_plot(mod, pred = "x", modx = "m", modx.values = "quartiles")

Probe simple slopes from a multilevel interaction

Description

Computes simple slopes of a focal predictor ('pred') at selected values of a moderator ('modx') from a two-level mixed-effects model fitted with [lme4::lmer()]. Returns estimates, standard errors, *t*-values, *p*-values, and confidence intervals in a tidy data frame.

Usage

mlm_probe(
  model,
  pred,
  modx,
  modx.values = c("mean-sd", "quartiles", "tertiles", "custom"),
  at = NULL,
  conf.level = 0.95,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual"),
  modx.level = c("auto", "cluster", "observation")
)

Arguments

Value

An object of class 'mlm_probe' (a list) with components: * 'slopes': a data frame with columns 'modx_value', 'slope', 'se', 't', 'df', 'p', 'ci_lower', 'ci_upper'. * 'df_method': the degrees-of-freedom method used. * 'pred', 'modx': names of the predictor and moderator. * 'modx.values': the strategy used. * 'conf.level': the confidence level. * 'model': the original model (stored for downstream use).

References

Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. *Biometrics Bulletin*, 2(6), 110–114.

Snijders, T. A. B., & Bosker, R. J. (2012). *Multilevel analysis* (2nd ed.). Sage.

Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. *Journal of Statistical Software*, 82(13), 1–26. doi:10.18637/jss.v082.i13

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
mlm_probe(mod, pred = "x", modx = "m")
mlm_probe(mod, pred = "x", modx = "m", modx.values = "quartiles")
mlm_probe(mod, pred = "x", modx = "m", modx.values = "custom", at = c(-1, 0, 1))

Leave-one-cluster-out influence diagnostics for a cross-level interaction

Description

Refits the model once for each cluster, leaving that cluster out, and records how the interaction coefficient changes. This answers a specific question: *is the interaction driven by a small number of clusters?*

Usage

mlm_sensitivity(
  model,
  pred,
  modx,
  alpha = 0.05,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual"),
  verbose = FALSE,
  icc_range,
  icc_grid,
  loco = TRUE
)

Arguments

Details

Each refit uses [stats::update()] on the data the model was fitted to (via [lme4::getData()]), so it keeps the original formula, estimation method (REML or ML), weights, and control settings. Influence is measured by a DFBETA-type statistic,

\mathrm{DFBETA}_j = (\hat\beta_{3(-j)} - \hat\beta_3) / SE(\hat\beta_3),

the change in the interaction when cluster j is omitted, in units of the full-data standard error. Clusters with |\mathrm{DFBETA}_j| > 2/\sqrt{J} are flagged, following the usual size-adjusted cutoff for DFBETAS (Belsley, Kuh, and Welsch, 1980). The cutoff is a screening heuristic, not a test.

**Changes in 0.3.0.** Earlier versions also reported an "ICC-shift" analysis that rescaled the interaction's standard error by a design-effect ratio for hypothetical intraclass correlations. That formula applies to means under a random-intercept model, not to a cross-level interaction in a random-slope model, and the accompanying adjusted Johnson-Neyman boundary was not correct, so the analysis has been removed. The earlier leave-one-cluster-out refits used maximum likelihood even when the original model used REML, and refitted from the model frame, which failed silently for formulas with transformed variables; both are fixed.

**Scope.** These are influence diagnostics. They do not address unmeasured confounding of the interaction.

Value

An object of class 'mlm_sensitivity' with components: * 'observed': list with the full-data interaction estimate, SE, df, p-value, and number of clusters. * 'loco': data frame with one row per cluster: 'cluster', 'n_obs', 'b_int', 'se_int', 'df', 'p_int', 'b_change', 'dfbeta', 'influential', 'same_sign', 'sig', and 'error' (the error message when a refit failed, otherwise 'NA'). * 'cutoff': the DFBETA flagging cutoff 2/\sqrt{J}. * 'summary': list with the number of failed refits, the proportion of successful refits keeping the sign and significance of the interaction, and the range of refitted estimates. * Metadata: 'pred', 'modx', 'alpha', 'df_method', 'grp_name', 'int_term'.

References

Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). *Regression diagnostics*. Wiley.

Examples

set.seed(42)
n_j <- 20; n_i <- 10
dat_small <- data.frame(
  y   = rnorm(n_j * n_i),
  x   = rnorm(n_j * n_i),
  m   = rep(rnorm(n_j), each = n_i),
  grp = factor(rep(seq_len(n_j), each = n_i))
)
dat_small$y <- dat_small$y + 0.5 * dat_small$x * dat_small$m
mod_small <- lme4::lmer(y ~ x * m + (1 | grp), data = dat_small)

sens <- mlm_sensitivity(mod_small, pred = "x", modx = "m",
                        df_method = "between")
sens
plot(sens)

Summary table for a multilevel moderation effect

Description

Returns a consolidated summary of the moderation effect: the focal interaction coefficient, simple slopes at selected moderator values, and (optionally) the Johnson—Neyman interval. Designed for quick reporting and results sections.

Usage

mlm_summary(
  model,
  pred,
  modx,
  modx.values = c("mean-sd", "quartiles", "tertiles", "custom"),
  at = NULL,
  conf.level = 0.95,
  jn = TRUE,
  alpha = 0.05,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual"),
  modx.level = c("auto", "cluster", "observation")
)

Arguments

Value

An object of class 'mlm_summary' (a list) with components: * 'interaction': one-row data frame for the interaction term. * 'simple_slopes': data frame from 'mlm_probe()'. * 'jn': output of 'mlm_jn()' (or 'NULL' if 'jn = FALSE'). * Other metadata.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
mlm_summary(mod, pred = "x", modx = "m")

Contour plot of predicted outcomes over the predictor x moderator space

Description

Plots iso-outcome contour lines of \hat{Y}(x, w) over the full joint space of 'pred' (x-axis) and 'modx' (y-axis). This is the most direct geometric representation of a two-way interaction:

Usage

mlm_surface(
  model,
  pred,
  modx,
  grid = 80L,
  n_contours = 10L,
  fill = TRUE,
  probe_lines = TRUE,
  x_label = NULL,
  y_label = NULL,
  legend_title = NULL
)

Arguments

Details

* **No interaction**: contour lines are perfectly straight and parallel — the effect of 'pred' does not depend on 'modx'. * **Positive interaction**: contour lines fan outward (rotate clockwise) — higher 'modx' steepens the 'pred' slope. * **Negative interaction**: contour lines fan inward (rotate counter-clockwise).

The degree of non-parallelism among contours is a direct visual index of interaction strength: the larger \beta_3, the more the lines rotate.

An optional overlay draws the three standard simple-slope evaluation lines (mean – 1 SD of 'modx') as horizontal reference lines, connecting the plot to 'mlm_probe()' output.

Predicted values are computed from fixed effects only ('re.form = NA'), with all covariates held at their means or reference levels. The surface therefore represents the population-average predicted outcome, not any specific cluster.

**Reading the plot:** Pick any contour line. Its slope in the 'pred' direction tells you how fast the outcome changes with 'pred' at that 'modx' value. If the contour slopes steeply up-right, 'pred' has a strong positive effect there. If contours become more horizontal as 'modx' increases, the 'pred' effect is weakening. If they rotate from positive to flat to negative, you have a sign-changing interaction — and the 'modx' value where they are perfectly horizontal is the Johnson-Neyman boundary.

Value

A 'ggplot' object.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
mlm_surface(mod, pred = "x", modx = "m")
mlm_surface(mod, pred = "x", modx = "m", fill = FALSE, n_contours = 15)

Confidence and prediction intervals for simple slopes in a random-slope model

Description

In a model with a random slope for 'pred', two different questions can be asked about the slope at a given moderator value w:

Usage

mlm_variance_decomp(
  model,
  pred,
  modx,
  modx.values = c("mean-sd", "quartiles", "tertiles", "custom"),
  at = NULL,
  conf.level = 0.95,
  df_method = c("satterthwaite", "kenward-roger", "between", "residual"),
  modx.level = c("auto", "cluster", "observation")
)

Arguments

Details

1. **How precisely is the average slope known?** The confidence interval for \beta_1 + \beta_3 w reflects only estimation uncertainty in the fixed effects (as in [mlm_probe()]). 2. **What slope should be expected in a new cluster?** The prediction interval adds the residual random-slope variance \tau_{11}:

\hat\beta_1 + \hat\beta_3 w \pm t_{df}\sqrt{SE^2(w) + \hat\tau_{11}}.

Because 'modx' is in the model, \tau_{11} is the slope variance that remains *after* the moderator: a significant interaction does not mean the moderator explains the between-cluster slope differences, and a large \tau_{11} relative to the average slope means individual clusters can have slopes far from it.

The prediction interval treats \hat\tau_{11} as known and assumes normally distributed random slopes, so it is approximate and too narrow when the number of clusters is small. The random slope is taken from the grouping factor whose random-effects terms include 'pred'.

Versions before 0.3.0 reported 'pct_random', the ratio \tau_{11} / (\tau_{11} + SE^2). It mixed a population variance with a sampling variance, so it grew with sample size without any change in the data-generating process; it has been removed.

Value

An object of class 'mlm_variance_decomp' (a list) with: * 'decomp': data frame with columns 'modx_value', 'slope', 'se_fixed' (SE of the average slope), 'tau11' (random-slope variance), 'tau11_sd' (its square root), 'se_total' (\sqrt{SE^2 + \tau_{11}}), 'df', 'ci_lower', 'ci_upper' (confidence interval for the average slope), and 'pi_lower', 'pi_upper' (prediction interval for a new cluster). * 'tau11', 'tau11_sd': the random-slope variance and SD. * 'has_random_slope': logical; does the model include a random slope for 'pred'? * Metadata: 'pred', 'modx', 'conf.level', 'df_method', 'grp_name'.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 + x | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
vd <- mlm_variance_decomp(mod, pred = "x", modx = "m")
print(vd)

Plot a Johnson-Neyman interval from a multilevel model

Description

Creates a ggplot2 figure showing the simple slope of 'pred' across the full range of 'modx', with shading indicating regions of significance. A vertical dashed line marks each Johnson-Neyman boundary.

Usage

## S3 method for class 'mlm_jn'
plot(
  x,
  x_label = NULL,
  y_label = NULL,
  sig_color = "#2166AC",
  nonsig_color = "#D6604D",
  ...
)

Arguments

Value

A 'ggplot' object.

Examples

set.seed(1)
dat <- data.frame(
  y   = rnorm(200), x = rnorm(200),
  m   = rep(rnorm(20), each = 10),
  grp = factor(rep(1:20, each = 10))
)
dat$y <- dat$y + dat$x * dat$m
mod <- lme4::lmer(y ~ x * m + (1 | grp), data = dat,
                  control = lme4::lmerControl(optimizer = "bobyqa"))
jn <- mlm_jn(mod, pred = "x", modx = "m")
plot(jn)

Plot leave-one-cluster-out influence for a cross-level interaction

Description

Plots the DFBETA of each cluster, ordered by size, with the flagging cutoff \pm 2/\sqrt{J} as dashed lines.

Usage

## S3 method for class 'mlm_sensitivity'
plot(x, ...)

Arguments

Value

A ggplot object.

Plot the variance decomposition of simple slopes

Description

Shows the simple slope at each moderator value as a point with two interval layers: an inner confidence interval (fixed-effect uncertainty only) and an outer prediction interval (fixed + random-slope variance). The gap between the two intervals is the contribution of random-slope heterogeneity.

Usage

## S3 method for class 'mlm_variance_decomp'
plot(x, x_label = NULL, y_label = NULL, ...)

Arguments

Value

A 'ggplot' object.

Simulated school achievement dataset

Description

A simulated two-level dataset with students nested within schools, designed to illustrate multilevel moderation analysis. The true data-generating model includes a cross-level interaction between student socioeconomic status and school climate.

Usage

school_data

Format

A data frame with 3,000 rows and 6 variables:

school

A factor indicating the school identifier (1–100).

student

An integer indicating the student identifier (1–3000).

math

A numeric mathematics achievement score.

ses

A numeric student socioeconomic status variable (standardized).

climate

A numeric school climate rating (standardized level-2 variable).

gender

A factor indicating student gender with levels '"female"' and '"male"'.

Details

The data were generated from the model

math_{ij} = 50 + 1.5\,ses_{ij} + 0.8\,climate_j + 0.5\,ses_{ij} climate_j + u_{0j} + u_{1j} ses_{ij} + e_{ij}.

The level-2 random effects were generated as u_{0j} \sim N(0, 9) and u_{1j} \sim N(0, 0.25), and the level-1 residuals were generated as e_{ij} \sim N(0, 25).

Examples

data(school_data)
head(school_data)
str(school_data)


library(lme4)
mod <- lmer(math ~ ses * climate + gender + (1 + ses | school),
            data = school_data)
mlm_probe(mod, pred = "ses", modx = "climate")