Random Slopes in lavaan

Multilevel structural equation modeling

lavaan

Estimating two-level random slope models with lavaan, compared against lme4 and Mplus under three centering schemes. Includes a benchmark of the Gauss–Hermite quadrature node count.

Author

Mark Lai

Published

September 18, 2026

This post shows how to fit two-level random slope models in lavaan using the rv() syntax (Rosseel 2012), and compares results against lme4 and Mplus under three centering schemes: uncentered, observed mean centered, and latent mean centered. The last section benchmarks the Gauss–Hermite (GH) quadrature that lavaan uses when the random-slope predictor is a latent variable.

Data

Code

library(mlmRev)
library(lme4)
library(lavaan)
library(MplusAutomation)
data(Hsb82, package = "mlmRev")
d <- na.omit(Hsb82[, c("mAch", "ses", "cses", "meanses", "school")])

We use the Hsb82 data from the mlmRev package: 7185 students nested in 160 schools. The outcome is math achievement (mAch), the predictor is socioeconomic status (ses), and the data already contain meanses (cluster mean) and cses (= ses − meanses, the cluster-mean-centered variable).

Code

# Shared helpers for extracting estimates from lavaan / Mplus and building
# the comparison tables.
get_lav_se <- function(pe, lhs, op, rhs, level = NULL) {
  i <- which(pe$lhs == lhs & pe$op == op & pe$rhs == rhs)
  if (!is.null(level)) i <- i[pe$level[i] == level]
  pe$se[i]
}
mp_load <- function(inp, out) {
  if (!file.exists(out)) return(NULL)
  suppressMessages(mplusModel(inp_file = inp, out_file = out, read = TRUE))
}
mp_par <- function(m) m$parameters$unstandardized
mp_get <- function(p, hdr, prm, lvl, se_col = "se") {
  i <- which(p$paramHeader == hdr & p$param == prm & p$BetweenWithin == lvl)
  list(est = p$est[i], se = p[[se_col]][i])
}
# covariance between two params, regardless of which is the Mplus header
mp_cov <- function(p, a, b, lvl, se_col = "se") {
  i <- which(p$BetweenWithin == lvl &
    ((p$paramHeader == paste0(a, ".WITH") & p$param == b) |
     (p$paramHeader == paste0(b, ".WITH") & p$param == a)))
  list(est = p$est[i], se = p[[se_col]][i])
}
# pull a 6-row (est, se) vector from a Mplus results table via spec vectors
mplus_row <- function(p, specs, se_col = "se") {
  f <- function(s) if (s[1] == "cov") mp_cov(p, s[2], s[3], s[4], se_col)
                   else               mp_get(p, s[2], s[3], s[4], se_col)
  r <- lapply(specs, f)
  list(est = sapply(r, `[[`, "est"), se = sapply(r, `[[`, "se"))
}
build_cmp <- function(names, lavaan, lavaan_se, lme4 = NULL, lme4_se = NULL,
                      mplus = NULL, mplus_se = NULL) {
  # strip names so cbind does not turn them into row names
  args <- list(names = unname(names), lavaan = unname(lavaan),
               lavaan_se = unname(lavaan_se), lme4 = unname(lme4),
               lme4_se = unname(lme4_se), mplus = unname(mplus),
               mplus_se = unname(mplus_se))
  out <- data.frame(Parameter = args$names)
  if (!is.null(args$lme4))  out <- cbind(out, lme4 = args$lme4, se_lme4 = args$lme4_se)
  out <- cbind(out, lavaan = args$lavaan, se_lavaan = args$lavaan_se)
  if (!is.null(args$mplus)) out <- cbind(out, mplus = args$mplus, se_mplus = args$mplus_se)
  rownames(out) <- NULL
  out
}

Model and centering options

The two-level random slope model is (Raudenbush and Bryk 2002)

\[ \begin{aligned} \text{Level 1:} \quad y_{ij} &= \beta_{0j} + \beta_{1j}\,x_{ij} + r_{ij}, \quad r_{ij} \sim N(0,\,\sigma^2) \\[4pt] \text{Level 2:} \quad \boldsymbol{\beta}_j &= \boldsymbol{\beta} + \mathbf{u}_j, \quad \mathbf{u}_j \sim N(\mathbf{0},\, \mathbf{T}) \end{aligned} \]

The level-1 predictor decomposes as $x_{ij} = \bar{x}_j + w_{ij}$ (school mean + within-school deviation). The centering scheme determines what the random slope $\beta_{1j}$ multiplies and how the between-level coefficient is interpreted (Lüdtke et al. 2008).

Uncentered

The random slope applies to the raw $x_{ij}$. The school mean $\bar{x}_j$ enters at level 2; its coefficient is the contextual contrast—the difference between the between-level and within-level effects:

\[y_{ij} = \beta_{00} + \beta_{01}\,\bar{x}_j + \beta_{10}\,x_{ij} + u_{0j} + u_{1j}\,x_{ij} + r_{ij}\]

Within effect: $\beta_{10} + u_{1j}$
Between effect: $\beta_{10} + \beta_{01} + u_{1j}$
Contextual contrast: $\beta_{01}$ = between − within

Observed mean centered

The random slope applies to the within deviation $w_{ij} = x_{ij} - \bar{x}_j$ (centered by the observed school mean). The school mean $\bar{x}_j$ enters at level 2 with its own (fixed) coefficient—the pure between effect:

\[y_{ij} = \beta_{00} + \beta_{1b}\,\bar{x}_j + \beta_{1w}\,(x_{ij} - \bar{x}_j) + u_{0j} + u_{1j}\,(x_{ij} - \bar{x}_j) + r_{ij}\]

Within effect: $\beta_{1w} + u_{1j}$
Between effect: $\beta_{1b}$ (fixed)

Latent mean centered

As in observed mean centering, but the school mean is a latent variable $B_j$ (correcting for unreliability of the observed school mean due to finite cluster size). The random slope applies to $x_{ij} - B_j$, and $B_j$ enters as a between-level predictor with its own (fixed) slope:

\[y_{ij} = \beta_{00} + \beta_{1b}\,B_j + \beta_{1w}\,(x_{ij} - B_j) + u_{0j} + u_{1j}\,(x_{ij} - B_j) + r_{ij}\]

Within effect: $\beta_{1w} + u_{1j}$
Between effect: $\beta_{1b}$ (fixed, on the latent mean)

In all three cases the predictor is observed, so the integral over the random effects has a closed-form solution—no numerical quadrature is needed.

Uncentered

The random slope is on the raw ses, and meanses enters at level 2 as the contextual predictor.

Code

fit_u_lme4 <- lmer(mAch ~ ses + meanses + (ses | school), data = d, REML = FALSE)

# lavaan: meanses is between-only, so it must appear at level 1
model_u <- '
  level: 1
    mAch ~ rv("b1") * ses + meanses
  level: 2
    mAch ~~ b1
'
fit_u_lav <- suppressWarnings(
  sem(model_u, data = d, cluster = "school", estimator = "ML")
)

cat("Log-likelihood:  lme4 =", round(as.numeric(logLik(fit_u_lme4)), 3),
    "  lavaan =", round(as.numeric(logLik(fit_u_lav)), 3), "\n")

Log-likelihood:  lme4 = -23278.46   lavaan = -23278.46

Code

pe <- parameterEstimates(fit_u_lav, se = TRUE)
cf <- coef(fit_u_lav)
fe <- fixef(fit_u_lme4)
vc <- as.matrix(VarCorr(fit_u_lme4)$school)
se_lme4 <- summary(fit_u_lme4)$coefficients[, "Std. Error"]
cmp_u <- build_cmp(
  names = c("Slope (ses)", "Contextual (meanses)", "Rand. intercept var",
            "Rand. slope var", "Int–slope cov", "Residual var"),
  lavaan = c(cf["b1~1.l2"], cf["mAch~meanses"], cf["mAch~~mAch.l2"],
             cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]),
  lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","meanses"),
                get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2),
                get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1)),
  lme4 = c(fe["ses"], fe["meanses"], vc[1,1], vc[2,2], vc[1,2], sigma(fit_u_lme4)^2),
  lme4_se = c(se_lme4["ses"], se_lme4["meanses"], NA, NA, NA, NA)
)
m <- mp_load("mplus/model_rs_uncentered.inp", "mplus/model_rs_uncentered.out")
if (!is.null(m)) {
  r <- mplus_row(mp_par(m), list(
    c("get", "Means", "B1", "Between"),
    c("get", "MACH.ON", "MEANSES", "Between"),
    c("get", "Residual.Variances", "MACH", "Between"),
    c("get", "Variances", "B1", "Between"),
    c("cov", "MACH", "B1", "Between"),
    c("get", "Residual.Variances", "MACH", "Within")))
  cmp_u <- cbind(cmp_u, mplus = r$est, se_mplus = r$se)
}
knitr::kable(cmp_u, digits = 4)

Parameter	lme4	se_lme4	lavaan	se_lavaan	mplus	se_mplus
Slope (ses)	2.1902	0.1214	2.1902	0.1214	2.191	0.121
Contextual (meanses)	3.7779	0.3802	3.7779	0.3820	3.782	0.382
Rand. intercept var	2.6488	NA	2.6488	0.4029	2.651	0.403
Rand. slope var	0.4367	NA	0.4367	0.2313	0.438	0.231
Int–slope cov	-0.2347	NA	-0.2347	0.2293	-0.234	0.230
Residual var	36.7970	NA	36.7970	0.6280	36.796	0.628

Note

The meanses coefficient (3.78) is the contextual contrast: the difference between the between-level effect (5.97) and the within-level effect (2.19).

Equivalent Mplus input (ML estimator):

VARIABLE:
  NAMES ARE mAch ses meanses school;
  CLUSTER IS school;
  BETWEEN IS meanses;
  WITHIN IS ses;
ANALYSIS:
  TYPE = TWOLEVEL RANDOM;
  ESTIMATOR = ML;
MODEL:
  %WITHIN%
    b1 | mAch ON ses;
  %BETWEEN%
    mAch ON meanses;
    mAch WITH b1;

Observed mean centered

The random slope is on cses (centered by the observed school mean), and meanses enters at level 2 as the between-level predictor.

Code

fit_o_lme4 <- lmer(mAch ~ cses + meanses + (cses | school), data = d, REML = FALSE)

# lavaan
model_o <- '
  level: 1
    mAch ~ rv("b1") * cses + meanses
  level: 2
    mAch ~~ b1
'
fit_o_lav <- suppressWarnings(
  sem(model_o, data = d, cluster = "school", estimator = "ML")
)

cat("Log-likelihood:  lme4 =", round(as.numeric(logLik(fit_o_lme4)), 3),
    "  lavaan =", round(as.numeric(logLik(fit_o_lav)), 3), "\n")

Log-likelihood:  lme4 = -23276.59   lavaan = -23276.59

Code

pe <- parameterEstimates(fit_o_lav, se = TRUE)
cf <- coef(fit_o_lav)
fe <- fixef(fit_o_lme4)
vc <- as.matrix(VarCorr(fit_o_lme4)$school)
se_lme4 <- summary(fit_o_lme4)$coefficients[, "Std. Error"]
cmp_o <- build_cmp(
  names = c("Within slope (cses)", "Between slope (meanses)", "Rand. intercept var",
            "Rand. slope var", "Int–slope cov", "Residual var"),
  lavaan = c(cf["b1~1.l2"], cf["mAch~meanses"], cf["mAch~~mAch.l2"],
             cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]),
  lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","meanses"),
                get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2),
                get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1)),
  lme4 = c(fe["cses"], fe["meanses"], vc[1,1], vc[2,2], vc[1,2], sigma(fit_o_lme4)^2),
  lme4_se = c(se_lme4["cses"], se_lme4["meanses"], NA, NA, NA, NA)
)
m <- mp_load("mplus/model_rs_cses.inp", "mplus/model_rs_cses.out")
if (!is.null(m)) {
  r <- mplus_row(mp_par(m), list(
    c("get", "Means", "B1", "Between"),
    c("get", "MACH.ON", "MEANSES", "Between"),
    c("get", "Residual.Variances", "MACH", "Between"),
    c("get", "Variances", "B1", "Between"),
    c("cov", "MACH", "B1", "Between"),
    c("get", "Residual.Variances", "MACH", "Within")))
  cmp_o <- cbind(cmp_o, mplus = r$est, se_mplus = r$se)
}
knitr::kable(cmp_o, digits = 4)

Parameter	lme4	se_lme4	lavaan	se_lavaan	mplus	se_mplus
Within slope (cses)	2.1912	0.1276	2.1912	0.1276	2.191	0.128
Between slope (meanses)	5.8959	0.3577	5.8959	0.3588	5.895	0.359
Rand. intercept var	2.6478	NA	2.6478	0.3966	2.648	0.397
Rand. slope var	0.6701	NA	0.6701	0.2767	0.699	0.284
Int–slope cov	-0.2638	NA	-0.2638	0.2398	-0.261	0.243
Residual var	36.7133	NA	36.7133	0.6260	36.706	0.626

Here the meanses coefficient (5.9) is the pure between effect—the association between the school-mean SES and the school-mean outcome, net of the within-school association.

Equivalent Mplus input (using the pre-centered cses variable):

VARIABLE:
  NAMES ARE mAch cses meanses school;
  CLUSTER IS school;
  BETWEEN IS meanses;
  WITHIN IS cses;
ANALYSIS:
  TYPE = TWOLEVEL RANDOM;
  ESTIMATOR = ML;
MODEL:
  %WITHIN%
    b1 | mAch ON cses;
  %BETWEEN%
    mAch ON meanses;
    mAch WITH b1;

Latent mean centered

In lavaan, placing ses at both levels implements latent mean centering: the level-2 ses is a latent variable (with its own variance), the random slope applies to the within deviation from this latent mean, and the between-level slope is a fixed effect on the latent mean.

Code

model_l <- '
  level: 1
    mAch ~ rv("b1") * ses
  level: 2
    mAch ~ ses
    b1 ~~ ses
    mAch ~~ b1
'
# 7 GH nodes give the same LL as the default of 21, at 3x speed
fit_l_lav <- sem(model_l, data = d, cluster = "school", estimator = "ML",
                 integration.ngh = 7)
summary(fit_l_lav, fit.measures = FALSE)

lavaan 0.7-2 ended normally after 77 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        10

  Number of observations                          7185
  Number of clusters [school]                      160


Parameter Estimates:

  Standard errors                             Standard
  Information                                 Observed
  Observed information based on                Hessian


Level 1 [within]:

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  mAch ~                                              
    ses               0.000                           

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .mAch             36.711    0.626   58.644    0.000


Level 2 [school]:

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  mAch ~                                              
    ses               6.111    0.381   16.023    0.000

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  b1 ~~                                               
    ses               0.050    0.054    0.927    0.354
   .mAch             -0.268    0.237   -1.133    0.257

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .mAch             12.684    0.148   85.602    0.000
    ses              -0.007    0.033   -0.203    0.839
    b1                2.195    0.128   17.177    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .mAch              2.479    0.396    6.254    0.000
    ses               0.161    0.019    8.392    0.000
    b1                0.672    0.277    2.426    0.015

The Mplus equivalent requires ESTIMATOR = BAYES, which implements latent mean centering for random predictors. The actual Mplus input is

VARIABLE:
  NAMES ARE mAch ses school;
  CLUSTER IS school;
  MISSING ARE ALL (99);
ANALYSIS:
  TYPE = TWOLEVEL RANDOM;
  ESTIMATOR = BAYES;
  BCONVERGENCE = .005;
  BITER = 100000(50000);
  ALGORITHM = GIBBS(RW);
  CHAIN = 4;
  PROCESS = 4;
MODEL:
  %WITHIN%
    b1 | mAch ON ses;
  %BETWEEN%
    mAch ON ses;
    b1 WITH ses;
    mAch WITH b1;
    ses;

Code

pe <- parameterEstimates(fit_l_lav, se = TRUE)
cf <- coef(fit_l_lav)
cmp_l <- build_cmp(
  names = c("Within slope", "Between slope", "Rand. intercept var",
            "Rand. slope var", "Int–slope cov", "Residual var"),
  lavaan = c(cf["b1~1.l2"], cf["mAch~ses.l2"], cf["mAch~~mAch.l2"],
             cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]),
  lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","ses",2),
                get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2),
                get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1))
)
m <- mp_load("mplus/model_rs_bayes.inp", "mplus/model_rs_bayes.out")
if (!is.null(m)) {
  r <- mplus_row(mp_par(m), list(
    c("get", "Means", "B1", "Between"),
    c("get", "MACH.ON", "SES", "Between"),
    c("get", "Residual.Variances", "MACH", "Between"),
    c("get", "Variances", "B1", "Between"),
    c("cov", "MACH", "B1", "Between"),
    c("get", "Residual.Variances", "MACH", "Within")), se_col = "posterior_sd")
  cmp_l <- cbind(cmp_l, mplus = r$est, se_mplus = r$se)
}
knitr::kable(cmp_l, digits = 4)

Parameter	lavaan	se_lavaan	mplus	se_mplus
Within slope	2.1948	0.1278	2.197	0.132
Between slope	6.1106	0.3814	6.097	0.394
Rand. intercept var	2.4789	0.3964	2.736	0.423
Rand. slope var	0.6722	0.2771	0.805	0.289
Int–slope cov	-0.2682	0.2368	-0.287	0.253
Residual var	36.7107	0.6260	36.718	0.626

Code

cat("\nlavaan log-likelihood:", round(as.numeric(logLik(fit_l_lav)), 3), "\n")

lavaan log-likelihood: -30796.65

Note

Mplus BAYES reports posterior means and standard deviations (not ML estimates and SEs). The fixed effects (slopes) agree well; the random-effects variance components can differ due to the tendency of the posterior median/mean in the Bayesian posterior to be larger than the mode (which is similar to the MLE).

Code

# The Mplus ML "hybrid": random slope on the whole observed ses, but a latent
# between-level ses. Loaded only to quote its estimates below.
m <- mp_load("mplus/model_rs_hybrid.inp", "mplus/model_rs_hybrid.out")
hyb_between <- if (!is.null(m)) mp_get(mp_par(m), "MACH.ON", "SES", "Between")$est else NA
hyb_t11     <- if (!is.null(m)) mp_get(mp_par(m), "Variances", "B1", "Between")$est else NA

Warning

ESTIMATOR = ML in Mplus is not true latent mean centering. With the ML estimator, Mplus applies the random slope to the whole observed predictor (uncentered) while treating the between-level predictor as a latent variable, and prints the warning “the random regression predictor variable on the WITHIN level refers to the whole observed variable.” This hybrid—a combination of latent mean centering and the uncentered method—was used in multilevel mediation models with random slopes (Preacher et al. 2010). True latent mean centering (random slope on $x_{ij} - B_j$) requires ESTIMATOR = BAYES (Asparouhov and Muthén 2018). The two are not equivalent: on Hsb82 the hybrid estimates a between slope of 4.00 and a slope variance $T_{11}$ of 0.41, versus 6.11 and 0.67 for true latent mean centering.

Comparing the three schemes

Code

cf_u <- coef(fit_u_lav); cf_o <- coef(fit_o_lav); cf_l <- coef(fit_l_lav)
tbl <- data.frame(
  Scheme  = c("Uncentered", "Observed mean", "Latent mean"),
  within  = c(cf_u["b1~1.l2"], cf_o["b1~1.l2"], cf_l["b1~1.l2"]),
  between = c(cf_u["mAch~meanses"], cf_o["mAch~meanses"], cf_l["mAch~ses.l2"]),
  T00 = c(cf_u["mAch~~mAch.l2"], cf_o["mAch~~mAch.l2"], cf_l["mAch~~mAch.l2"]),
  T11 = c(cf_u["b1~~b1.l2"], cf_o["b1~~b1.l2"], cf_l["b1~~b1.l2"]),
  ll  = c(as.numeric(logLik(fit_u_lav)), as.numeric(logLik(fit_o_lav)),
          as.numeric(logLik(fit_l_lav)))
)
tbl$within  <- round(tbl$within, 3)
tbl$between <- round(tbl$between, 3)
tbl$T00     <- round(tbl$T00, 3)
tbl$T11     <- round(tbl$T11, 3)
tbl$ll      <- round(tbl$ll, 1)
colnames(tbl) <- c("Scheme", "Within slope", "Between slope", "T00", "T11", "logLik")
knitr::kable(tbl)

Scheme	Within slope	Between slope	T00	T11	logLik
Uncentered	2.190	3.778	2.649	0.437	-23278.5
Observed mean	2.191	5.896	2.648	0.670	-23276.6
Latent mean	2.195	6.111	2.479	0.672	-30796.6

The within slope is stable across schemes (~2.19). The between-level coefficient has a different meaning in each scheme: in the uncentered model it is the contextual contrast (between − within = 3.78); in the observed mean centered model it is the pure between effect (5.9); in the latent mean centered model it is the pure between effect on the latent mean (6.11). The random-effects covariance $\mathbf{T}$ changes substantially because the random intercept absorbs different amounts of between-level variation.

Gauss–Hermite quadrature: when it is needed and how many nodes

When the random-slope predictor is a latent variable (e.g., a factor measured by multiple indicators), the closed-form integration is no longer available (Rockwood 2020). lavaan then uses Gauss–Hermite quadrature over the nonlinear random effects. The integration.ngh option controls the number of 1-D nodes; for $q$ latent random slopes the total number of evaluation points is $n_{\text{nodes}}^{\,q}$.

Each evaluation point is a closed-form conditional log-likelihood, so the total cost scales linearly in the number of nodes.

Benchmark

I benchmarked the node count on three two-level random slope models, each fit over a range of integration.ngh values and compared against a high-node reference. The tables are read from a saved results file (gh_benchmark_results.rds); the simulation that produced it—including the slow 2-D Model B—is in the appendix. Model A, a single random slope whose predictor is split into a between and a within part so that the level-2 predictor is latent ($G = 100$ clusters of $n_j = 20$), is the reference case:

ngh	nodes	Δ LL	time (s)	speedup vs 21
1	1	43.4	0.6	21.1
2	2	1.87	1.4	9.6
3	3	1.53	2.4	5.4
5	5	0.559	3.1	4.1
7	7	0.149	3.6	3.6
12	12	0.00668	6.7	1.9
21	21	0.000373	13.0	1.0

Key observations (Model A):

1 node (evaluation at the prior mean) is catastrophically wrong.
7 nodes give |ΔLL| ≈ 0.15—typically below the threshold that affects parameter estimates at the second decimal.
12 nodes give |ΔLL| < 0.01, which is negligible in practice.
21 nodes (the current default) is already overkill for most applications.

The pattern holds across all three models. The table below gives the total number of quadrature evaluation points needed to reach a given $|\Delta LL|$ (for the 2-D Model B this is an $n \times n$ grid, i.e. $n^2$ points):

model	$q_{nl}$	points for Δ LL < 0.1	points for Δ LL < 0.01
A	1	12	12
C	1	7	12
B	2	16	36

Only Model B behaves differently, and only because its integral is 2-D: an $n$-node setting means $n^2$ evaluation points. It is still accurate to $|\Delta LL| < 0.1$ with a small grid, but it is the one case where the node count has a real runtime effect—which is exactly where the trade-off above matters.

Recommendation

Setting integration.ngh = 12 in lavaan options provides an excellent accuracy/speed trade-off for most applications. Combined with vectorisation improvements (which reduce the per-node constant by ~1.5× for the objective), the net speedup for a typical latent-covariate random slope model is ~3×.

Appendix: Gauss–Hermite quadrature simulation

The benchmark tables above based on the simulation below: three two-level random slope models, each fit over a range of integration.ngh values and compared to a high-node reference. The simulation code is shown below

Code

# gh_benchmark.R---regenerate gh_benchmark_results.rds
suppressPackageStartupMessages(library(lavaan))
mvrnorm1 <- function(n, mu, Sigma) {
  L <- chol(t(Sigma)); Z <- matrix(rnorm(length(mu) * n), ncol = length(mu))
  t(t(Z %*% L) + mu)
}
# fit at a node count; mute only the benign transient negative-variance warning
fit_gh <- function(model, data, cluster, ngh) {
  t0 <- Sys.time()
  f <- withCallingHandlers(
    sem(model, data, cluster = cluster, estimator = "ML",
        integration.ngh = ngh, se = "none", verbose = FALSE),
    warning = function(w) if (grepl("ov variances are negative", conditionMessage(w), fixed = TRUE))
      invokeRestart("muffleWarning"))
  c(ll = as.numeric(logLik(f)), time = as.numeric(difftime(Sys.time(), t0, units = "secs")))
}
gh_sweep <- function(model, data, cluster, nghs, ref_ngh, qnl) {
  ref <- fit_gh(model, data, cluster, ref_ngh)
  rows <- lapply(nghs, function(g) { r <- fit_gh(model, data, cluster, g)
    c(ngh = g, nodes = if (qnl > 1) g^2 else g,
      ll = r[["ll"]], abs_dll = abs(r[["ll"]] - ref[["ll"]]), time = r[["time"]]) })
  res <- as.data.frame(do.call(rbind, rows))
  res$speedup <- res$time[res$ngh == max(nghs)] / res$time
  res
}
# Model A: simple 1 random slope, split observed predictor (q.nl = 1)
gen_a <- function(G, n) {
  N <- G*n; x<-y<-id<-numeric(N); idx<-1
  for (j in 1:G) { xj<-rnorm(1,0,1); bj<-rnorm(1,0.5,0.3); uj<-rnorm(1,0,0.5)
    xw<-rnorm(n,0,1); ey<-rnorm(n,0,1)
    x[idx:(idx+n-1)] <- xj + xw; y[idx:(idx+n-1)] <- 0.3*xj + bj*xw + uj + ey
    id[idx:(idx+n-1)] <- j; idx<-idx+n }
  data.frame(x=x, y=y, id=id)
}
model_A <- '
  level: 1
    y ~ rv("s1") * x
  level: 2
    y ~ x
    s1 ~~ x
    y ~~ s1
'
# Model C: latent regression, lavaan tutorial model (q.nl = 1)
gen_c <- function(G, n) {
  N <- G*n; x1<-x2<-x3<-y1<-y2<-y3<-id<-numeric(N); idx<-1
  lamx<-c(1.0,0.8,0.7); lamy<-c(1.0,0.9,0.8); th<-0.5
  gx<-c(1.0,0.8,0.7); gy<-c(1.0,0.9,0.8)
  for (j in 1:G) { s1<-rnorm(1,0.5,0.3); fxb<-rnorm(1,0,1)
    fyb<-0.5*fxb+0.3*s1+rnorm(1,0,sqrt(0.3))
    for (i in 1:n) { fxw<-rnorm(1,0,1); fyw<-s1*fxw+rnorm(1,0,sqrt(0.5))
      x1[idx]<-gx[1]*fxb+lamx[1]*fxw+rnorm(1,0,sqrt(th))
      x2[idx]<-gx[2]*fxb+lamx[2]*fxw+rnorm(1,0,sqrt(th))
      x3[idx]<-gx[3]*fxb+lamx[3]*fxw+rnorm(1,0,sqrt(th))
      y1[idx]<-gy[1]*fyb+lamy[1]*fyw+rnorm(1,0,sqrt(th))
      y2[idx]<-gy[2]*fyb+lamy[2]*fyw+rnorm(1,0,sqrt(th))
      y3[idx]<-gy[3]*fyb+lamy[3]*fyw+rnorm(1,0,sqrt(th))
      id[idx]<-j; idx<-idx+1 } }
  data.frame(x1=x1,x2=x2,x3=x3,y1=y1,y2=y2,y3=y3,id=id)
}
model_C <- '
  level: 1
    fxw =~ x1 + x2 + x3
    fyw =~ y1 + y2 + y3
    fyw ~ rv("s1") * fxw
  level: 2
    fxb =~ x1 + x2 + x3
    fyb =~ y1 + y2 + y3
    fyb ~ fxb
    fyb ~~ s1
'
# Model B: two correlated latent random slopes (q.nl = 2)
gen_b <- function(G, n) {
  N <- G*n
  mus<-c(0.5,0.3); Sigmas<-matrix(c(0.09,0.02,0.02,0.04),2)
  Sigb<-matrix(c(1,0.3,0.3,1),2); Sigw<-matrix(c(1,0.3,0.3,1),2)
  lamx<-c(1.0,0.8,0.7); lamy<-c(1.0,0.9,0.8); th<-0.5
  gx<-c(1.0,0.8,0.7); gy<-c(1.0,0.9,0.8)
  x1a<-x1b<-x1c<-x2a<-x2b<-x2c<-y1<-y2<-y3<-id<-numeric(N); idx<-1
  for (j in 1:G) { s<-mvrnorm1(1,mus,Sigmas); s1j<-s[1]; s2j<-s[2]
    fb<-mvrnorm1(1,c(0,0),Sigb); b1j<-fb[1]; b2j<-fb[2]
    fyb<-0.4*b1j+0.3*b2j+0.2*s1j+0.1*s2j+rnorm(1,0,sqrt(0.3))
    for (i in 1:n) { fw<-mvrnorm1(1,c(0,0),Sigw); w1<-fw[1]; w2<-fw[2]
      fyw<-s1j*w1+s2j*w2+rnorm(1,0,sqrt(0.5))
      x1a[idx]<-gx[1]*b1j+lamx[1]*w1+rnorm(1,0,sqrt(th))
      x1b[idx]<-gx[2]*b1j+lamx[2]*w1+rnorm(1,0,sqrt(th))
      x1c[idx]<-gx[3]*b1j+lamx[3]*w1+rnorm(1,0,sqrt(th))
      x2a[idx]<-gx[1]*b2j+lamx[1]*w2+rnorm(1,0,sqrt(th))
      x2b[idx]<-gx[2]*b2j+lamx[2]*w2+rnorm(1,0,sqrt(th))
      x2c[idx]<-gx[3]*b2j+lamx[3]*w2+rnorm(1,0,sqrt(th))
      y1[idx]<-gy[1]*fyb+lamy[1]*fyw+rnorm(1,0,sqrt(th))
      y2[idx]<-gy[2]*fyb+lamy[2]*fyw+rnorm(1,0,sqrt(th))
      y3[idx]<-gy[3]*fyb+lamy[3]*fyw+rnorm(1,0,sqrt(th))
      id[idx]<-j; idx<-idx+1 } }
  data.frame(x1a=x1a,x1b=x1b,x1c=x1c,x2a=x2a,x2b=x2b,x2c=x2c,y1=y1,y2=y2,y3=y3,id=id)
}
model_B <- '
  level: 1
    fxw1 =~ x1a + x1b + x1c
    fxw2 =~ x2a + x2b + x2c
    fyw  =~ y1 + y2 + y3
    fyw ~ rv("s1") * fxw1 + rv("s2") * fxw2
  level: 2
    fxb1 =~ x1a + x1b + x1c
    fxb2 =~ x2a + x2b + x2c
    fyb  =~ y1 + y2 + y3
    fyb ~ fxb1 + fxb2
    fyb ~~ s1 + s2
    s1 ~~ s2
'
set.seed(42); dA <- gen_a(100, 20); ra <- gh_sweep(model_A, dA, "id", c(1,2,3,5,7,12,21), 50, 1)
set.seed(44); dC <- gen_c(100, 10); rc <- gh_sweep(model_C, dC, "id", c(1,3,7,12,21), 30, 1)
set.seed(43); dB <- gen_b(30, 8);  rb <- gh_sweep(model_B, dB, "id", c(2,3,4,5,6), 7, 2)
saveRDS(list(A = ra, C = rc, B = rb, meta = list(
  A = list(label = "Simple 1 random slope (split observed predictor)", qnl = 1, N = 2000, ref = 50),
  C = list(label = "Latent regression (lavaan tutorial model)",         qnl = 1, N = 1000, ref = 30),
  B = list(label = "Two correlated latent random slopes",               qnl = 2, N = 240,  ref = 7)
)), "gh_benchmark_results.rds")

Code

ghres <- readRDS("gh_benchmark_results.rds")
tab_gh <- function(r) data.frame(
  ngh = r$ngh, nodes = r$nodes,
  dll = formatC(r$abs_dll, format = "g", digits = 3),
  time = round(r$time, 1), speedup = round(r$speedup, 1))
lab <- function(k) { m <- ghres$meta[[k]]
  sprintf("Model %s -- %s (q.nl = %d, N = %d, reference %d nodes)",
          k, m$label, m$qnl, m$N, if (m$qnl > 1) m$ref^2 else m$ref) }
cat(lab("A"), "\n"); knitr::kable(tab_gh(ghres$A),
  col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 21"))

Model A -- Simple 1 random slope (split observed predictor) (q.nl = 1, N = 2000, reference 50 nodes)

ngh	nodes	Δ LL	time (s)	speedup vs 21
1	1	43.4	0.6	21.1
2	2	1.87	1.4	9.6
3	3	1.53	2.4	5.4
5	5	0.559	3.1	4.1
7	7	0.149	3.6	3.6
12	12	0.00668	6.7	1.9
21	21	0.000373	13.0	1.0

Code

cat(lab("C"), "\n"); knitr::kable(tab_gh(ghres$C),
  col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 21"))

Model C -- Latent regression (lavaan tutorial model) (q.nl = 1, N = 1000, reference 30 nodes)

ngh	nodes	Δ LL	time (s)	speedup vs 21
1	1	10.4	2.7	16.1
3	3	0.43	8.6	5.0
7	7	0.016	15.0	2.9
12	12	0.00121	28.1	1.5
21	21	7.39e-08	43.3	1.0

Code

cat(lab("B"), "\n"); knitr::kable(tab_gh(ghres$B),
  col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 6"))

Model B -- Two correlated latent random slopes (q.nl = 2, N = 240, reference 49 nodes)

ngh	nodes	Δ LL	time (s)	speedup vs 6
2	4	0.187	12.0	5.7
3	9	0.152	21.6	3.2
4	16	0.0455	29.9	2.3
5	25	0.0162	59.8	1.1
6	36	0.00144	68.6	1.0

Appendix: Reproduction in OpenMx

OpenMx fits the same three schemes with a hand-written RAM parameterisation: a byClus submodel holds the cluster-level random effects, joined to the student-level data by cluster id. All three are reproduced to numerical precision against lme4 (schemes 1–2) and Mplus (the hybrid); the hybrid’s non-significant intercept–slope covariance is the only parameter that differs noticeably from Mplus.

Code

library(OpenMx)

To take full advantage of multiple cores, use:
  mxOption(key='Number of Threads', value=parallel::detectCores()) #now
  Sys.setenv(OMP_NUM_THREADS=parallel::detectCores()) #before library(OpenMx)

Attaching package: 'OpenMx'

The following objects are masked from 'package:Matrix':

    %&%, expm

Code

od <- d
od$cid <- as.integer(factor(d$school))
od_b <- data.frame(cid = sort(unique(od$cid)))

# free-parameter estimates from a fitted OpenMx model + its submodels
omx_est <- function(model) {
  out <- NULL
  for (mm in model@matrices) {
    vals <- slot(mm, "values"); free <- slot(mm, "free"); lab <- slot(mm, "labels")
    nm <- slot(mm, "name"); if (is.null(vals) || is.null(free)) next
    idx <- which(free, arr.ind = TRUE)
    for (k in seq_len(nrow(idx))) {
      i <- idx[k, 1]; j <- idx[k, 2]
      out[if (!is.null(lab) && !is.na(lab[i, j])) lab[i, j] else paste0(nm, "_", i, "_", j)] <- vals[i, j]
    }
  }
  for (s in model@submodels) out <- c(out, omx_est(s))
  out
}
fit3 <- function(m) { f <- mxRun(m, silent = TRUE); stopifnot(f@output$status$code == 0)
  est <- omx_est(f); attr(est, "logLik") <- as.numeric(logLik(f)); est }

# (1) uncentered: random slope on raw ses, observed meanses at L2
e_unc <- fit3(mxModel("unc",
  mxModel("byClus", type = "RAM", latentVars = c("u0", "v"),
    data = mxData(od_b, "raw", primaryKey = "cid"),
    mxPath(from = c("u0", "v"), arrows = 2, free = TRUE, values = c(2.6, 0.44),
           labels = c("varInt", "varSlope")),
    mxPath(from = "u0", to = "v", arrows = 2, free = TRUE, values = -0.25, labels = "covIntSlope")),
  type = "RAM", manifestVars = c("y", "ses", "meanses"),
  data = mxData(data.frame(cid = od$cid, y = od$mAch, ses = od$ses, meanses = od$meanses), "raw", sort = FALSE),
  mxPath(from = "one", to = c("y", "ses", "meanses"), free = TRUE,
         values = c(mean(od$mAch), mean(od$ses), mean(od$meanses)), labels = c("meanY", "meanX", "meanB")),
  mxPath(from = c("y", "ses", "meanses"), arrows = 2, free = TRUE,
         values = c(36, var(od$ses), var(od$meanses)), labels = c("residY", "residX", "residB")),
  mxPath(from = "ses", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"),
  mxPath(from = "meanses", to = "y", free = TRUE, values = 3.8, labels = "contextual"),
  mxPath(from = c("byClus.u0", "byClus.v"), to = c("y", "y"), free = FALSE,
         values = c(1, NA), labels = c(NA, "data.ses"), joinKey = "cid")))

# (2) observed mean centered: random slope on cses, meanses at L2
e_obs <- fit3(mxModel("obs",
  mxModel("byClus", type = "RAM", latentVars = c("u0", "v"),
    data = mxData(od_b, "raw", primaryKey = "cid"),
    mxPath(from = c("u0", "v"), arrows = 2, free = TRUE, values = c(2.65, 0.67),
           labels = c("varInt", "varSlope")),
    mxPath(from = "u0", to = "v", arrows = 2, free = TRUE, values = -0.26, labels = "covIntSlope")),
  type = "RAM", manifestVars = c("y", "cses", "meanses"),
  data = mxData(data.frame(cid = od$cid, y = od$mAch, cses = od$cses, meanses = od$meanses), "raw", sort = FALSE),
  mxPath(from = "one", to = c("y", "cses", "meanses"), free = TRUE,
         values = c(mean(od$mAch), mean(od$cses), mean(od$meanses)), labels = c("meanY", "meanXw", "meanB")),
  mxPath(from = c("y", "cses", "meanses"), arrows = 2, free = TRUE,
         values = c(36, var(od$cses), var(od$meanses)), labels = c("residY", "residXw", "residB")),
  mxPath(from = "cses", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"),
  mxPath(from = "meanses", to = "y", free = TRUE, values = 5.9, labels = "betweenSlope"),
  mxPath(from = c("byClus.u0", "byClus.v"), to = c("y", "y"), free = FALSE,
         values = c(1, NA), labels = c(NA, "data.cses"), joinKey = "cid")))

# (3) Mplus hybrid: random slope on raw ses, between-level ses is latent B
e_hyb <- fit3(mxModel("hyb",
  mxModel("byClus", type = "RAM", latentVars = c("B", "intc", "v"),
    data = mxData(od_b, "raw", primaryKey = "cid"),
    mxPath(from = c("B", "intc", "v"), arrows = 2, free = TRUE,
           values = c(0.16, 2.5, 0.44), labels = c("varB", "varIntR", "varSlope")),
    mxPath(from = "B", to = "intc", free = TRUE, values = 4.0, labels = "contextual"),
    mxPath(from = "B", to = "v", arrows = 2, free = TRUE, values = 0.03, labels = "covSlopeB"),
    mxPath(from = "intc", to = "v", arrows = 2, free = TRUE, values = -0.36, labels = "covIntcSlope")),
  type = "RAM", manifestVars = c("y", "x"),
  data = mxData(data.frame(cid = od$cid, y = od$mAch, x = od$ses), "raw", sort = FALSE),
  mxPath(from = "one", to = c("y", "x"), free = TRUE, values = c(mean(od$mAch), mean(od$ses)),
         labels = c("meanY", "meanX")),
  mxPath(from = c("y", "x"), arrows = 2, free = TRUE, values = c(36, var(od$ses)),
         labels = c("residY", "residX")),
  mxPath(from = "x", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"),
  mxPath(from = c("byClus.intc", "byClus.v", "byClus.B"), to = c("y", "y", "x"),
         free = FALSE, values = c(1, NA, 1), labels = c(NA, "data.x", NA), joinKey = "cid")))

# ---- references -----------------------------------------------------------
fl_u <- lmer(mAch ~ ses + meanses + (ses | school), data = od, REML = FALSE)
fl_o <- lmer(mAch ~ cses + meanses + (cses | school), data = od, REML = FALSE)
vc_u <- as.matrix(VarCorr(fl_u)$school); vc_o <- as.matrix(VarCorr(fl_o)$school)
mph <- mp_load("mplus/model_rs_hybrid.inp", "mplus/model_rs_hybrid.out")
mg <- function(hdr, prm, lvl) {
  pp <- mp_par(mph)
  pp$est[pp$paramHeader == hdr & pp$param == prm & pp$BetweenWithin == lvl]
}

cat("### (1) Uncentered\n")

### (1) Uncentered

Code

t1 <- round(data.frame(
  Within = c(e_unc[["withinSlope"]], fixef(fl_u)["ses"]),
  Contextual = c(e_unc[["contextual"]], fixef(fl_u)["meanses"]),
  T00 = c(e_unc[["varInt"]], vc_u[1, 1]), T11 = c(e_unc[["varSlope"]], vc_u[2, 2]),
  T01 = c(e_unc[["covIntSlope"]], vc_u[1, 2]), Resid = c(e_unc[["residY"]], sigma(fl_u)^2)),
  3); rownames(t1) <- c("OpenMx", "lme4")
knitr::kable(t1, row.names = TRUE)

	Within	Contextual	T00	T11	T01	Resid
OpenMx	2.19	3.778	2.649	0.437	-0.235	36.797
lme4	2.19	3.778	2.649	0.437	-0.235	36.797

Code

cat("### (2) Observed mean centered\n")

### (2) Observed mean centered

Code

t2 <- round(data.frame(
  Within = c(e_obs[["withinSlope"]], fixef(fl_o)["cses"]),
  Between = c(e_obs[["betweenSlope"]], fixef(fl_o)["meanses"]),
  T00 = c(e_obs[["varInt"]], vc_o[1, 1]), T11 = c(e_obs[["varSlope"]], vc_o[2, 2]),
  T01 = c(e_obs[["covIntSlope"]], vc_o[1, 2]), Resid = c(e_obs[["residY"]], sigma(fl_o)^2)),
  3); rownames(t2) <- c("OpenMx", "lme4")
knitr::kable(t2, row.names = TRUE)

	Within	Between	T00	T11	T01	Resid
OpenMx	2.191	5.896	2.648	0.67	-0.264	36.713
lme4	2.191	5.896	2.648	0.67	-0.264	36.713

Code

cat("### (3) Mplus hybrid\n")

### (3) Mplus hybrid

Code

t3 <- round(data.frame(
  Within = c(e_hyb[["withinSlope"]], mg("Means", "B1", "Between")),
  Between = c(e_hyb[["contextual"]], mg("MACH.ON", "SES", "Between")),
  varB = c(e_hyb[["varB"]], mg("Variances", "SES", "Between")),
  T00 = c(e_hyb[["varIntR"]], mg("Residual.Variances", "MACH", "Between")),
  T11 = c(e_hyb[["varSlope"]], mg("Variances", "B1", "Between")),
  T01 = c(e_hyb[["covIntcSlope"]] - e_hyb[["contextual"]] * e_hyb[["covSlopeB"]],
          mg("B1.WITH", "MACH", "Between")),
  Resid = c(e_hyb[["residY"]], mg("Residual.Variances", "MACH", "Within")),
  LogLik = c(round(as.numeric(attr(e_hyb, "logLik")), 1), round(mph$summaries$LL, 1))),
  3); rownames(t3) <- c("OpenMx", "Mplus")
knitr::kable(t3, row.names = TRUE)

	Within	Between	varB	T00	T11	T01	Resid	LogLik
OpenMx	2.193	4.013	0.16	2.503	0.438	-0.360	36.792	-30798.8
Mplus	2.197	4.004	0.16	2.500	0.414	-0.256	36.801	-30798.8

The two log-likelihoods agree to within optimizer precision (-30798.8): the OpenMx and Mplus hybrid models are the same model fit to the same data by maximum likelihood.

Asparouhov, Thomas, and Bengt Muthén. 2018. “Latent Variable Centering of Predictors and Mediators in Multilevel and Time-Series Models.” Structural Equation Modeling: A Multidisciplinary Journal 26 (1): 119–42. https://doi.org/10.1080/10705511.2018.1511375.

Lüdtke, Oliver, Herbert W. Marsh, Alexander Robitzsch, Ulrich Trautwein, Tihomir Asparouhov, and Bengt Muthén. 2008. “The Multilevel Latent Covariate Model: A New, More Reliable Approach to Group-Level Effects in Contextual Studies.” Psychological Methods 13 (3): 203–29. https://doi.org/10.1037/a0012869.

Preacher, Kristopher J., Michael J. Zyphur, and Zhen Zhang. 2010. “A General Multilevel SEM Framework for Assessing Multilevel Mediation.” Psychological Methods 15 (3): 209–33. https://doi.org/10.1037/a0020141.

Raudenbush, Stephen W., and Anthony S. Bryk. 2002. Hierarchical Linear Models: Applications and Design Strategies. 2nd ed. Sage.

Rockwood, Nicholas J. 2020. “Maximum Likelihood Estimation of Multilevel Structural Equation Models with Random Slopes for Latent Covariates.” Psychometrika 85 (2): 275–300. https://doi.org/10.1007/s11336-020-09702-9.

Rosseel, Yves. 2012. “lavaan: An R Package for Structural Equation Modeling.” Journal of Statistical Software 48 (2): 1–36. https://doi.org/10.18637/jss.v048.i02.

--- title: "Random Slopes in lavaan" author: - Mark Lai bibliography: references.bib categories: - Multilevel structural equation modeling - lavaan date: 2026-09-18 description: | Estimating two-level random slope models with lavaan, compared against lme4 and Mplus under three centering schemes. Includes a benchmark of the Gauss--Hermite quadrature node count. --- This post shows how to fit two-level random slope models in `lavaan` using the `rv()` syntax (Rosseel 2012), and compares results against `lme4` and Mplus under three centering schemes: uncentered, observed mean centered, and latent mean centered. The last section benchmarks the Gauss--Hermite (GH) quadrature that `lavaan` uses when the random-slope predictor is a latent variable. ## Data ::: cell ``` {.r .cell-code} library(mlmRev) library(lme4) library(lavaan) library(MplusAutomation) data(Hsb82, package = "mlmRev") d <- na.omit(Hsb82[, c("mAch", "ses", "cses", "meanses", "school")]) ``` ::: We use the `Hsb82` data from the `mlmRev` package: 7185 students nested in 160 schools. The outcome is math achievement (`mAch`), the predictor is socioeconomic status (`ses`), and the data already contain `meanses` (cluster mean) and `cses` (= `ses` − `meanses`, the cluster-mean-centered variable). ::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} # Shared helpers for extracting estimates from lavaan / Mplus and building # the comparison tables. get_lav_se <- function(pe, lhs, op, rhs, level = NULL) { i <- which(pe$lhs == lhs & pe$op == op & pe$rhs == rhs) if (!is.null(level)) i <- i[pe$level[i] == level] pe$se[i] } mp_load <- function(inp, out) { if (!file.exists(out)) return(NULL) suppressMessages(mplusModel(inp_file = inp, out_file = out, read = TRUE)) } mp_par <- function(m) m$parameters$unstandardized mp_get <- function(p, hdr, prm, lvl, se_col = "se") { i <- which(p$paramHeader == hdr & p$param == prm & p$BetweenWithin == lvl) list(est = p$est[i], se = p[[se_col]][i]) } # covariance between two params, regardless of which is the Mplus header mp_cov <- function(p, a, b, lvl, se_col = "se") { i <- which(p$BetweenWithin == lvl & ((p$paramHeader == paste0(a, ".WITH") & p$param == b) | (p$paramHeader == paste0(b, ".WITH") & p$param == a))) list(est = p$est[i], se = p[[se_col]][i]) } # pull a 6-row (est, se) vector from a Mplus results table via spec vectors mplus_row <- function(p, specs, se_col = "se") { f <- function(s) if (s[1] == "cov") mp_cov(p, s[2], s[3], s[4], se_col) else mp_get(p, s[2], s[3], s[4], se_col) r <- lapply(specs, f) list(est = sapply(r, `[[`, "est"), se = sapply(r, `[[`, "se")) } build_cmp <- function(names, lavaan, lavaan_se, lme4 = NULL, lme4_se = NULL, mplus = NULL, mplus_se = NULL) { # strip names so cbind does not turn them into row names args <- list(names = unname(names), lavaan = unname(lavaan), lavaan_se = unname(lavaan_se), lme4 = unname(lme4), lme4_se = unname(lme4_se), mplus = unname(mplus), mplus_se = unname(mplus_se)) out <- data.frame(Parameter = args$names) if (!is.null(args$lme4)) out <- cbind(out, lme4 = args$lme4, se_lme4 = args$lme4_se) out <- cbind(out, lavaan = args$lavaan, se_lavaan = args$lavaan_se) if (!is.null(args$mplus)) out <- cbind(out, mplus = args$mplus, se_mplus = args$mplus_se) rownames(out) <- NULL out } ``` </details> ::: ## Model and centering options The two-level random slope model is (Raudenbush and Bryk 2002) $$ \begin{aligned} \text{Level 1:} \quad y_{ij} &= \beta_{0j} + \beta_{1j}\,x_{ij} + r_{ij}, \quad r_{ij} \sim N(0,\,\sigma^2) \\[4pt] \text{Level 2:} \quad \boldsymbol{\beta}_j &= \boldsymbol{\beta} + \mathbf{u}_j, \quad \mathbf{u}_j \sim N(\mathbf{0},\, \mathbf{T}) \end{aligned} $$ The level-1 predictor decomposes as $x_{ij} = \bar{x}_j + w_{ij}$ (school mean + within-school deviation). The centering scheme determines what the random slope $\beta_{1j}$ multiplies and how the between-level coefficient is interpreted (Lüdtke et al. 2008). ### Uncentered The random slope applies to the raw $x_{ij}$. The school mean $\bar{x}_j$ enters at level 2; its coefficient is the **contextual contrast**---the difference between the between-level and within-level effects: $$y_{ij} = \beta_{00} + \beta_{01}\,\bar{x}_j + \beta_{10}\,x_{ij} + u_{0j} + u_{1j}\,x_{ij} + r_{ij}$$ - Within effect: $\beta_{10} + u_{1j}$ - Between effect: $\beta_{10} + \beta_{01} + u_{1j}$ - Contextual contrast: $\beta_{01}$ = between − within ### Observed mean centered The random slope applies to the within deviation $w_{ij} = x_{ij} - \bar{x}_j$ (centered by the **observed** school mean). The school mean $\bar{x}_j$ enters at level 2 with its own (fixed) coefficient---the **pure between effect**: $$y_{ij} = \beta_{00} + \beta_{1b}\,\bar{x}_j + \beta_{1w}\,(x_{ij} - \bar{x}_j) + u_{0j} + u_{1j}\,(x_{ij} - \bar{x}_j) + r_{ij}$$ - Within effect: $\beta_{1w} + u_{1j}$ - Between effect: $\beta_{1b}$ (fixed) ### Latent mean centered As in observed mean centering, but the school mean is a **latent** variable $B_j$ (correcting for unreliability of the observed school mean due to finite cluster size). The random slope applies to $x_{ij} - B_j$, and $B_j$ enters as a between-level predictor with its own (fixed) slope: $$y_{ij} = \beta_{00} + \beta_{1b}\,B_j + \beta_{1w}\,(x_{ij} - B_j) + u_{0j} + u_{1j}\,(x_{ij} - B_j) + r_{ij}$$ - Within effect: $\beta_{1w} + u_{1j}$ - Between effect: $\beta_{1b}$ (fixed, on the latent mean) In all three cases the predictor is observed, so the integral over the random effects has a closed-form solution---no numerical quadrature is needed. ## Uncentered The random slope is on the raw `ses`, and `meanses` enters at level 2 as the contextual predictor. :::: cell ``` {.r .cell-code} fit_u_lme4 <- lmer(mAch ~ ses + meanses + (ses | school), data = d, REML = FALSE) # lavaan: meanses is between-only, so it must appear at level 1 model_u <- ' level: 1 mAch ~ rv("b1") * ses + meanses level: 2 mAch ~~ b1 ' fit_u_lav <- suppressWarnings( sem(model_u, data = d, cluster = "school", estimator = "ML") ) cat("Log-likelihood: lme4 =", round(as.numeric(logLik(fit_u_lme4)), 3), " lavaan =", round(as.numeric(logLik(fit_u_lav)), 3), "\n") ``` ::: {.cell-output .cell-output-stdout} Log-likelihood: lme4 = -23278.46 lavaan = -23278.46 ::: :::: :::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} pe <- parameterEstimates(fit_u_lav, se = TRUE) cf <- coef(fit_u_lav) fe <- fixef(fit_u_lme4) vc <- as.matrix(VarCorr(fit_u_lme4)$school) se_lme4 <- summary(fit_u_lme4)$coefficients[, "Std. Error"] cmp_u <- build_cmp( names = c("Slope (ses)", "Contextual (meanses)", "Rand. intercept var", "Rand. slope var", "Int–slope cov", "Residual var"), lavaan = c(cf["b1~1.l2"], cf["mAch~meanses"], cf["mAch~~mAch.l2"], cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]), lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","meanses"), get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2), get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1)), lme4 = c(fe["ses"], fe["meanses"], vc[1,1], vc[2,2], vc[1,2], sigma(fit_u_lme4)^2), lme4_se = c(se_lme4["ses"], se_lme4["meanses"], NA, NA, NA, NA) ) m <- mp_load("mplus/model_rs_uncentered.inp", "mplus/model_rs_uncentered.out") if (!is.null(m)) { r <- mplus_row(mp_par(m), list( c("get", "Means", "B1", "Between"), c("get", "MACH.ON", "MEANSES", "Between"), c("get", "Residual.Variances", "MACH", "Between"), c("get", "Variances", "B1", "Between"), c("cov", "MACH", "B1", "Between"), c("get", "Residual.Variances", "MACH", "Within"))) cmp_u <- cbind(cmp_u, mplus = r$est, se_mplus = r$se) } knitr::kable(cmp_u, digits = 4) ``` </details> ::: cell-output-display ----------------------------------------------------------------------------------- Parameter lme4 se_lme4 lavaan se_lavaan mplus se_mplus --------------------- --------- --------- --------- ----------- -------- ---------- Slope (ses) 2.1902 0.1214 2.1902 0.1214 2.191 0.121 Contextual (meanses) 3.7779 0.3802 3.7779 0.3820 3.782 0.382 Rand. intercept var 2.6488 NA 2.6488 0.4029 2.651 0.403 Rand. slope var 0.4367 NA 0.4367 0.2313 0.438 0.231 Int--slope cov -0.2347 NA -0.2347 0.2293 -0.234 0.230 Residual var 36.7970 NA 36.7970 0.6280 36.796 0.628 ----------------------------------------------------------------------------------- ::: :::: ::: {} > **Note** > > The `meanses` coefficient (3.78) is the **contextual contrast**: the > difference between the between-level effect (5.97) and the > within-level effect (2.19). ::: Equivalent Mplus input (ML estimator): VARIABLE: NAMES ARE mAch ses meanses school; CLUSTER IS school; BETWEEN IS meanses; WITHIN IS ses; ANALYSIS: TYPE = TWOLEVEL RANDOM; ESTIMATOR = ML; MODEL: %WITHIN% b1 | mAch ON ses; %BETWEEN% mAch ON meanses; mAch WITH b1; ## Observed mean centered The random slope is on `cses` (centered by the observed school mean), and `meanses` enters at level 2 as the between-level predictor. :::: cell ``` {.r .cell-code} fit_o_lme4 <- lmer(mAch ~ cses + meanses + (cses | school), data = d, REML = FALSE) # lavaan model_o <- ' level: 1 mAch ~ rv("b1") * cses + meanses level: 2 mAch ~~ b1 ' fit_o_lav <- suppressWarnings( sem(model_o, data = d, cluster = "school", estimator = "ML") ) cat("Log-likelihood: lme4 =", round(as.numeric(logLik(fit_o_lme4)), 3), " lavaan =", round(as.numeric(logLik(fit_o_lav)), 3), "\n") ``` ::: {.cell-output .cell-output-stdout} Log-likelihood: lme4 = -23276.59 lavaan = -23276.59 ::: :::: :::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} pe <- parameterEstimates(fit_o_lav, se = TRUE) cf <- coef(fit_o_lav) fe <- fixef(fit_o_lme4) vc <- as.matrix(VarCorr(fit_o_lme4)$school) se_lme4 <- summary(fit_o_lme4)$coefficients[, "Std. Error"] cmp_o <- build_cmp( names = c("Within slope (cses)", "Between slope (meanses)", "Rand. intercept var", "Rand. slope var", "Int–slope cov", "Residual var"), lavaan = c(cf["b1~1.l2"], cf["mAch~meanses"], cf["mAch~~mAch.l2"], cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]), lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","meanses"), get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2), get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1)), lme4 = c(fe["cses"], fe["meanses"], vc[1,1], vc[2,2], vc[1,2], sigma(fit_o_lme4)^2), lme4_se = c(se_lme4["cses"], se_lme4["meanses"], NA, NA, NA, NA) ) m <- mp_load("mplus/model_rs_cses.inp", "mplus/model_rs_cses.out") if (!is.null(m)) { r <- mplus_row(mp_par(m), list( c("get", "Means", "B1", "Between"), c("get", "MACH.ON", "MEANSES", "Between"), c("get", "Residual.Variances", "MACH", "Between"), c("get", "Variances", "B1", "Between"), c("cov", "MACH", "B1", "Between"), c("get", "Residual.Variances", "MACH", "Within"))) cmp_o <- cbind(cmp_o, mplus = r$est, se_mplus = r$se) } knitr::kable(cmp_o, digits = 4) ``` </details> ::: cell-output-display ------------------------------------------------------------------------------------- Parameter lme4 se_lme4 lavaan se_lavaan mplus se_mplus ----------------------- --------- --------- --------- ----------- -------- ---------- Within slope (cses) 2.1912 0.1276 2.1912 0.1276 2.191 0.128 Between slope (meanses) 5.8959 0.3577 5.8959 0.3588 5.895 0.359 Rand. intercept var 2.6478 NA 2.6478 0.3966 2.648 0.397 Rand. slope var 0.6701 NA 0.6701 0.2767 0.699 0.284 Int--slope cov -0.2638 NA -0.2638 0.2398 -0.261 0.243 Residual var 36.7133 NA 36.7133 0.6260 36.706 0.626 ------------------------------------------------------------------------------------- ::: :::: Here the `meanses` coefficient (5.9) is the **pure between effect**---the association between the school-mean SES and the school-mean outcome, net of the within-school association. Equivalent Mplus input (using the pre-centered `cses` variable): VARIABLE: NAMES ARE mAch cses meanses school; CLUSTER IS school; BETWEEN IS meanses; WITHIN IS cses; ANALYSIS: TYPE = TWOLEVEL RANDOM; ESTIMATOR = ML; MODEL: %WITHIN% b1 | mAch ON cses; %BETWEEN% mAch ON meanses; mAch WITH b1; ## Latent mean centered In `lavaan`, placing `ses` at both levels implements latent mean centering: the level-2 `ses` is a latent variable (with its own variance), the random slope applies to the within deviation from this latent mean, and the between-level slope is a fixed effect on the latent mean. :::: cell ``` {.r .cell-code} model_l <- ' level: 1 mAch ~ rv("b1") * ses level: 2 mAch ~ ses b1 ~~ ses mAch ~~ b1 ' # 7 GH nodes give the same LL as the default of 21, at 3x speed fit_l_lav <- sem(model_l, data = d, cluster = "school", estimator = "ML", integration.ngh = 7) summary(fit_l_lav, fit.measures = FALSE) ``` ::: {.cell-output .cell-output-stdout} lavaan 0.7-2 ended normally after 77 iterations Estimator ML Optimization method NLMINB Number of model parameters 10 Number of observations 7185 Number of clusters [school] 160 Parameter Estimates: Standard errors Standard Information Observed Observed information based on Hessian Level 1 [within]: Regressions: Estimate Std.Err z-value P(>|z|) mAch ~ ses 0.000 Variances: Estimate Std.Err z-value P(>|z|) .mAch 36.711 0.626 58.644 0.000 Level 2 [school]: Regressions: Estimate Std.Err z-value P(>|z|) mAch ~ ses 6.111 0.381 16.023 0.000 Covariances: Estimate Std.Err z-value P(>|z|) b1 ~~ ses 0.050 0.054 0.927 0.354 .mAch -0.268 0.237 -1.133 0.257 Intercepts: Estimate Std.Err z-value P(>|z|) .mAch 12.684 0.148 85.602 0.000 ses -0.007 0.033 -0.203 0.839 b1 2.195 0.128 17.177 0.000 Variances: Estimate Std.Err z-value P(>|z|) .mAch 2.479 0.396 6.254 0.000 ses 0.161 0.019 8.392 0.000 b1 0.672 0.277 2.426 0.015 ::: :::: The Mplus equivalent requires `ESTIMATOR = BAYES`, which implements latent mean centering for random predictors. The actual Mplus input is VARIABLE: NAMES ARE mAch ses school; CLUSTER IS school; MISSING ARE ALL (99); ANALYSIS: TYPE = TWOLEVEL RANDOM; ESTIMATOR = BAYES; BCONVERGENCE = .005; BITER = 100000(50000); ALGORITHM = GIBBS(RW); CHAIN = 4; PROCESS = 4; MODEL: %WITHIN% b1 | mAch ON ses; %BETWEEN% mAch ON ses; b1 WITH ses; mAch WITH b1; ses; ::::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} pe <- parameterEstimates(fit_l_lav, se = TRUE) cf <- coef(fit_l_lav) cmp_l <- build_cmp( names = c("Within slope", "Between slope", "Rand. intercept var", "Rand. slope var", "Int–slope cov", "Residual var"), lavaan = c(cf["b1~1.l2"], cf["mAch~ses.l2"], cf["mAch~~mAch.l2"], cf["b1~~b1.l2"], cf["b1~~mAch.l2"], cf["mAch~~mAch"]), lavaan_se = c(get_lav_se(pe, "b1","~1","",2), get_lav_se(pe, "mAch","~","ses",2), get_lav_se(pe, "mAch","~~","mAch",2), get_lav_se(pe, "b1","~~","b1",2), get_lav_se(pe, "b1","~~","mAch",2), get_lav_se(pe, "mAch","~~","mAch",1)) ) m <- mp_load("mplus/model_rs_bayes.inp", "mplus/model_rs_bayes.out") if (!is.null(m)) { r <- mplus_row(mp_par(m), list( c("get", "Means", "B1", "Between"), c("get", "MACH.ON", "SES", "Between"), c("get", "Residual.Variances", "MACH", "Between"), c("get", "Variances", "B1", "Between"), c("cov", "MACH", "B1", "Between"), c("get", "Residual.Variances", "MACH", "Within")), se_col = "posterior_sd") cmp_l <- cbind(cmp_l, mplus = r$est, se_mplus = r$se) } knitr::kable(cmp_l, digits = 4) ``` </details> ::: cell-output-display Parameter lavaan se_lavaan mplus se_mplus --------------------- --------- ----------- -------- ---------- Within slope 2.1948 0.1278 2.197 0.132 Between slope 6.1106 0.3814 6.097 0.394 Rand. intercept var 2.4789 0.3964 2.736 0.423 Rand. slope var 0.6722 0.2771 0.805 0.289 Int--slope cov -0.2682 0.2368 -0.287 0.253 Residual var 36.7107 0.6260 36.718 0.626 ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cat("\nlavaan log-likelihood:", round(as.numeric(logLik(fit_l_lav)), 3), "\n") ``` </details> ::: {.cell-output .cell-output-stdout} lavaan log-likelihood: -30796.65 ::: ::::: ::: {} > **Note** > > Mplus BAYES reports posterior means and standard deviations (not ML > estimates and SEs). The fixed effects (slopes) agree well; the > random-effects variance components can differ due to the tendency of > the posterior median/mean in the Bayesian posterior to be larger than > the mode (which is similar to the MLE). ::: ::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} # The Mplus ML "hybrid": random slope on the whole observed ses, but a latent # between-level ses. Loaded only to quote its estimates below. m <- mp_load("mplus/model_rs_hybrid.inp", "mplus/model_rs_hybrid.out") hyb_between <- if (!is.null(m)) mp_get(mp_par(m), "MACH.ON", "SES", "Between")$est else NA hyb_t11 <- if (!is.null(m)) mp_get(mp_par(m), "Variances", "B1", "Between")$est else NA ``` </details> ::: ::: {} > **Warning** > > **`ESTIMATOR = ML` in Mplus is not true latent mean centering.** With > the ML estimator, Mplus applies the random slope to the *whole > observed* predictor (uncentered) while treating the between-level > predictor as a latent variable, and prints the warning *"the random > regression predictor variable on the WITHIN level refers to the whole > observed variable."* This *hybrid*---a combination of latent mean > centering and the uncentered method---was used in multilevel mediation > models with random slopes (Preacher et al. 2010). True latent mean > centering (random slope on $x_{ij} - B_j$) requires > `ESTIMATOR = BAYES` (Asparouhov and Muthén 2018). The two are not > equivalent: on `Hsb82` the hybrid estimates a between slope of 4.00 > and a slope variance $T_{11}$ of 0.41, versus 6.11 and 0.67 for true > latent mean centering. ::: ## Comparing the three schemes :::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cf_u <- coef(fit_u_lav); cf_o <- coef(fit_o_lav); cf_l <- coef(fit_l_lav) tbl <- data.frame( Scheme = c("Uncentered", "Observed mean", "Latent mean"), within = c(cf_u["b1~1.l2"], cf_o["b1~1.l2"], cf_l["b1~1.l2"]), between = c(cf_u["mAch~meanses"], cf_o["mAch~meanses"], cf_l["mAch~ses.l2"]), T00 = c(cf_u["mAch~~mAch.l2"], cf_o["mAch~~mAch.l2"], cf_l["mAch~~mAch.l2"]), T11 = c(cf_u["b1~~b1.l2"], cf_o["b1~~b1.l2"], cf_l["b1~~b1.l2"]), ll = c(as.numeric(logLik(fit_u_lav)), as.numeric(logLik(fit_o_lav)), as.numeric(logLik(fit_l_lav))) ) tbl$within <- round(tbl$within, 3) tbl$between <- round(tbl$between, 3) tbl$T00 <- round(tbl$T00, 3) tbl$T11 <- round(tbl$T11, 3) tbl$ll <- round(tbl$ll, 1) colnames(tbl) <- c("Scheme", "Within slope", "Between slope", "T00", "T11", "logLik") knitr::kable(tbl) ``` </details> ::: cell-output-display Scheme Within slope Between slope T00 T11 logLik --------------- -------------- --------------- ------- ------- ---------- Uncentered 2.190 3.778 2.649 0.437 -23278.5 Observed mean 2.191 5.896 2.648 0.670 -23276.6 Latent mean 2.195 6.111 2.479 0.672 -30796.6 ::: :::: The within slope is stable across schemes (\~2.19). The between-level coefficient has a different meaning in each scheme: in the uncentered model it is the contextual contrast (between − within = 3.78); in the observed mean centered model it is the pure between effect (5.9); in the latent mean centered model it is the pure between effect on the latent mean (6.11). The random-effects covariance $\mathbf{T}$ changes substantially because the random intercept absorbs different amounts of between-level variation. ## Gauss--Hermite quadrature: when it is needed and how many nodes When the random-slope predictor is a **latent variable** (e.g., a factor measured by multiple indicators), the closed-form integration is no longer available (Rockwood 2020). `lavaan` then uses Gauss--Hermite quadrature over the nonlinear random effects. The `integration.ngh` option controls the number of 1-D nodes; for $q$ latent random slopes the total number of evaluation points is $n_{\text{nodes}}^{\,q}$. Each evaluation point is a closed-form conditional log-likelihood, so the total cost scales **linearly** in the number of nodes. ### Benchmark I benchmarked the node count on three two-level random slope models, each fit over a range of `integration.ngh` values and compared against a high-node reference. The tables are read from a saved results file (`gh_benchmark_results.rds`); the simulation that produced it---including the slow 2-D Model B---is in the appendix. **Model A**, a single random slope whose predictor is split into a between and a within part so that the level-2 predictor is latent ($G = 100$ clusters of $n_j = 20$), is the reference case: :::: cell ::: cell-output-display ngh nodes Δ LL time (s) speedup vs 21 ----- ------- ---------- ---------- --------------- 1 1 43.4 0.6 21.1 2 2 1.87 1.4 9.6 3 3 1.53 2.4 5.4 5 5 0.559 3.1 4.1 7 7 0.149 3.6 3.6 12 12 0.00668 6.7 1.9 21 21 0.000373 13.0 1.0 ::: :::: Key observations (Model A): - **1 node** (evaluation at the prior mean) is catastrophically wrong. - **7 nodes** give \|ΔLL\| ≈ 0.15---typically below the threshold that affects parameter estimates at the second decimal. - **12 nodes** give \|ΔLL\| \< 0.01, which is negligible in practice. - **21 nodes** (the current default) is already overkill for most applications. The pattern holds across all three models. The table below gives the total number of quadrature evaluation points needed to reach a given $|\Delta LL|$ (for the 2-D Model B this is an $n \times n$ grid, i.e. $n^2$ points): :::: cell ::: cell-output-display model $q_{nl}$ points for Δ LL \< 0.1 points for Δ LL \< 0.01 ------- ---------- ------------------------ ------------------------- A 1 12 12 C 1 7 12 B 2 16 36 ::: :::: Only Model B behaves differently, and only because its integral is 2-D: an $n$-node setting means $n^2$ evaluation points. It is still accurate to $|\Delta LL| < 0.1$ with a small grid, but it is the one case where the node count has a real runtime effect---which is exactly where the trade-off above matters. ### Recommendation Setting `integration.ngh = 12` in `lavaan` options provides an excellent accuracy/speed trade-off for most applications. Combined with vectorisation improvements (which reduce the per-node constant by \~1.5× for the objective), the net speedup for a typical latent-covariate random slope model is \~3×. # Appendix: Gauss--Hermite quadrature simulation The benchmark tables above based on the simulation below: three two-level random slope models, each fit over a range of `integration.ngh` values and compared to a high-node reference. The simulation code is shown below ::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} # gh_benchmark.R---regenerate gh_benchmark_results.rds suppressPackageStartupMessages(library(lavaan)) mvrnorm1 <- function(n, mu, Sigma) { L <- chol(t(Sigma)); Z <- matrix(rnorm(length(mu) * n), ncol = length(mu)) t(t(Z %*% L) + mu) } # fit at a node count; mute only the benign transient negative-variance warning fit_gh <- function(model, data, cluster, ngh) { t0 <- Sys.time() f <- withCallingHandlers( sem(model, data, cluster = cluster, estimator = "ML", integration.ngh = ngh, se = "none", verbose = FALSE), warning = function(w) if (grepl("ov variances are negative", conditionMessage(w), fixed = TRUE)) invokeRestart("muffleWarning")) c(ll = as.numeric(logLik(f)), time = as.numeric(difftime(Sys.time(), t0, units = "secs"))) } gh_sweep <- function(model, data, cluster, nghs, ref_ngh, qnl) { ref <- fit_gh(model, data, cluster, ref_ngh) rows <- lapply(nghs, function(g) { r <- fit_gh(model, data, cluster, g) c(ngh = g, nodes = if (qnl > 1) g^2 else g, ll = r[["ll"]], abs_dll = abs(r[["ll"]] - ref[["ll"]]), time = r[["time"]]) }) res <- as.data.frame(do.call(rbind, rows)) res$speedup <- res$time[res$ngh == max(nghs)] / res$time res } # Model A: simple 1 random slope, split observed predictor (q.nl = 1) gen_a <- function(G, n) { N <- G*n; x<-y<-id<-numeric(N); idx<-1 for (j in 1:G) { xj<-rnorm(1,0,1); bj<-rnorm(1,0.5,0.3); uj<-rnorm(1,0,0.5) xw<-rnorm(n,0,1); ey<-rnorm(n,0,1) x[idx:(idx+n-1)] <- xj + xw; y[idx:(idx+n-1)] <- 0.3*xj + bj*xw + uj + ey id[idx:(idx+n-1)] <- j; idx<-idx+n } data.frame(x=x, y=y, id=id) } model_A <- ' level: 1 y ~ rv("s1") * x level: 2 y ~ x s1 ~~ x y ~~ s1 ' # Model C: latent regression, lavaan tutorial model (q.nl = 1) gen_c <- function(G, n) { N <- G*n; x1<-x2<-x3<-y1<-y2<-y3<-id<-numeric(N); idx<-1 lamx<-c(1.0,0.8,0.7); lamy<-c(1.0,0.9,0.8); th<-0.5 gx<-c(1.0,0.8,0.7); gy<-c(1.0,0.9,0.8) for (j in 1:G) { s1<-rnorm(1,0.5,0.3); fxb<-rnorm(1,0,1) fyb<-0.5*fxb+0.3*s1+rnorm(1,0,sqrt(0.3)) for (i in 1:n) { fxw<-rnorm(1,0,1); fyw<-s1*fxw+rnorm(1,0,sqrt(0.5)) x1[idx]<-gx[1]*fxb+lamx[1]*fxw+rnorm(1,0,sqrt(th)) x2[idx]<-gx[2]*fxb+lamx[2]*fxw+rnorm(1,0,sqrt(th)) x3[idx]<-gx[3]*fxb+lamx[3]*fxw+rnorm(1,0,sqrt(th)) y1[idx]<-gy[1]*fyb+lamy[1]*fyw+rnorm(1,0,sqrt(th)) y2[idx]<-gy[2]*fyb+lamy[2]*fyw+rnorm(1,0,sqrt(th)) y3[idx]<-gy[3]*fyb+lamy[3]*fyw+rnorm(1,0,sqrt(th)) id[idx]<-j; idx<-idx+1 } } data.frame(x1=x1,x2=x2,x3=x3,y1=y1,y2=y2,y3=y3,id=id) } model_C <- ' level: 1 fxw =~ x1 + x2 + x3 fyw =~ y1 + y2 + y3 fyw ~ rv("s1") * fxw level: 2 fxb =~ x1 + x2 + x3 fyb =~ y1 + y2 + y3 fyb ~ fxb fyb ~~ s1 ' # Model B: two correlated latent random slopes (q.nl = 2) gen_b <- function(G, n) { N <- G*n mus<-c(0.5,0.3); Sigmas<-matrix(c(0.09,0.02,0.02,0.04),2) Sigb<-matrix(c(1,0.3,0.3,1),2); Sigw<-matrix(c(1,0.3,0.3,1),2) lamx<-c(1.0,0.8,0.7); lamy<-c(1.0,0.9,0.8); th<-0.5 gx<-c(1.0,0.8,0.7); gy<-c(1.0,0.9,0.8) x1a<-x1b<-x1c<-x2a<-x2b<-x2c<-y1<-y2<-y3<-id<-numeric(N); idx<-1 for (j in 1:G) { s<-mvrnorm1(1,mus,Sigmas); s1j<-s[1]; s2j<-s[2] fb<-mvrnorm1(1,c(0,0),Sigb); b1j<-fb[1]; b2j<-fb[2] fyb<-0.4*b1j+0.3*b2j+0.2*s1j+0.1*s2j+rnorm(1,0,sqrt(0.3)) for (i in 1:n) { fw<-mvrnorm1(1,c(0,0),Sigw); w1<-fw[1]; w2<-fw[2] fyw<-s1j*w1+s2j*w2+rnorm(1,0,sqrt(0.5)) x1a[idx]<-gx[1]*b1j+lamx[1]*w1+rnorm(1,0,sqrt(th)) x1b[idx]<-gx[2]*b1j+lamx[2]*w1+rnorm(1,0,sqrt(th)) x1c[idx]<-gx[3]*b1j+lamx[3]*w1+rnorm(1,0,sqrt(th)) x2a[idx]<-gx[1]*b2j+lamx[1]*w2+rnorm(1,0,sqrt(th)) x2b[idx]<-gx[2]*b2j+lamx[2]*w2+rnorm(1,0,sqrt(th)) x2c[idx]<-gx[3]*b2j+lamx[3]*w2+rnorm(1,0,sqrt(th)) y1[idx]<-gy[1]*fyb+lamy[1]*fyw+rnorm(1,0,sqrt(th)) y2[idx]<-gy[2]*fyb+lamy[2]*fyw+rnorm(1,0,sqrt(th)) y3[idx]<-gy[3]*fyb+lamy[3]*fyw+rnorm(1,0,sqrt(th)) id[idx]<-j; idx<-idx+1 } } data.frame(x1a=x1a,x1b=x1b,x1c=x1c,x2a=x2a,x2b=x2b,x2c=x2c,y1=y1,y2=y2,y3=y3,id=id) } model_B <- ' level: 1 fxw1 =~ x1a + x1b + x1c fxw2 =~ x2a + x2b + x2c fyw =~ y1 + y2 + y3 fyw ~ rv("s1") * fxw1 + rv("s2") * fxw2 level: 2 fxb1 =~ x1a + x1b + x1c fxb2 =~ x2a + x2b + x2c fyb =~ y1 + y2 + y3 fyb ~ fxb1 + fxb2 fyb ~~ s1 + s2 s1 ~~ s2 ' set.seed(42); dA <- gen_a(100, 20); ra <- gh_sweep(model_A, dA, "id", c(1,2,3,5,7,12,21), 50, 1) set.seed(44); dC <- gen_c(100, 10); rc <- gh_sweep(model_C, dC, "id", c(1,3,7,12,21), 30, 1) set.seed(43); dB <- gen_b(30, 8); rb <- gh_sweep(model_B, dB, "id", c(2,3,4,5,6), 7, 2) saveRDS(list(A = ra, C = rc, B = rb, meta = list( A = list(label = "Simple 1 random slope (split observed predictor)", qnl = 1, N = 2000, ref = 50), C = list(label = "Latent regression (lavaan tutorial model)", qnl = 1, N = 1000, ref = 30), B = list(label = "Two correlated latent random slopes", qnl = 2, N = 240, ref = 7) )), "gh_benchmark_results.rds") ``` </details> ::: ::::::::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} ghres <- readRDS("gh_benchmark_results.rds") tab_gh <- function(r) data.frame( ngh = r$ngh, nodes = r$nodes, dll = formatC(r$abs_dll, format = "g", digits = 3), time = round(r$time, 1), speedup = round(r$speedup, 1)) lab <- function(k) { m <- ghres$meta[[k]] sprintf("Model %s -- %s (q.nl = %d, N = %d, reference %d nodes)", k, m$label, m$qnl, m$N, if (m$qnl > 1) m$ref^2 else m$ref) } cat(lab("A"), "\n"); knitr::kable(tab_gh(ghres$A), col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 21")) ``` </details> ::: {.cell-output .cell-output-stdout} Model A -- Simple 1 random slope (split observed predictor) (q.nl = 1, N = 2000, reference 50 nodes) ::: ::: cell-output-display ngh nodes Δ LL time (s) speedup vs 21 ----- ------- ---------- ---------- --------------- 1 1 43.4 0.6 21.1 2 2 1.87 1.4 9.6 3 3 1.53 2.4 5.4 5 5 0.559 3.1 4.1 7 7 0.149 3.6 3.6 12 12 0.00668 6.7 1.9 21 21 0.000373 13.0 1.0 ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cat(lab("C"), "\n"); knitr::kable(tab_gh(ghres$C), col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 21")) ``` </details> ::: {.cell-output .cell-output-stdout} Model C -- Latent regression (lavaan tutorial model) (q.nl = 1, N = 1000, reference 30 nodes) ::: ::: cell-output-display ngh nodes Δ LL time (s) speedup vs 21 ----- ------- ---------- ---------- --------------- 1 1 10.4 2.7 16.1 3 3 0.43 8.6 5.0 7 7 0.016 15.0 2.9 12 12 0.00121 28.1 1.5 21 21 7.39e-08 43.3 1.0 ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cat(lab("B"), "\n"); knitr::kable(tab_gh(ghres$B), col.names = c("ngh", "nodes", "Δ LL", "time (s)", "speedup vs 6")) ``` </details> ::: {.cell-output .cell-output-stdout} Model B -- Two correlated latent random slopes (q.nl = 2, N = 240, reference 49 nodes) ::: ::: cell-output-display ngh nodes Δ LL time (s) speedup vs 6 ----- ------- --------- ---------- -------------- 2 4 0.187 12.0 5.7 3 9 0.152 21.6 3.2 4 16 0.0455 29.9 2.3 5 25 0.0162 59.8 1.1 6 36 0.00144 68.6 1.0 ::: ::::::::: # Appendix: Reproduction in OpenMx `OpenMx` fits the same three schemes with a hand-written `RAM` parameterisation: a `byClus` submodel holds the cluster-level random effects, joined to the student-level data by cluster id. All three are reproduced to numerical precision against `lme4` (schemes 1--2) and Mplus (the hybrid); the hybrid's non-significant intercept--slope covariance is the only parameter that differs noticeably from Mplus. :::::::::::: cell `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} library(OpenMx) ``` </details> ::: {.cell-output .cell-output-stderr} To take full advantage of multiple cores, use: mxOption(key='Number of Threads', value=parallel::detectCores()) #now Sys.setenv(OMP_NUM_THREADS=parallel::detectCores()) #before library(OpenMx) ::: ::: {.cell-output .cell-output-stderr} Attaching package: 'OpenMx' ::: ::: {.cell-output .cell-output-stderr} The following objects are masked from 'package:Matrix': %&%, expm ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} od <- d od$cid <- as.integer(factor(d$school)) od_b <- data.frame(cid = sort(unique(od$cid))) # free-parameter estimates from a fitted OpenMx model + its submodels omx_est <- function(model) { out <- NULL for (mm in model@matrices) { vals <- slot(mm, "values"); free <- slot(mm, "free"); lab <- slot(mm, "labels") nm <- slot(mm, "name"); if (is.null(vals) || is.null(free)) next idx <- which(free, arr.ind = TRUE) for (k in seq_len(nrow(idx))) { i <- idx[k, 1]; j <- idx[k, 2] out[if (!is.null(lab) && !is.na(lab[i, j])) lab[i, j] else paste0(nm, "_", i, "_", j)] <- vals[i, j] } } for (s in model@submodels) out <- c(out, omx_est(s)) out } fit3 <- function(m) { f <- mxRun(m, silent = TRUE); stopifnot(f@output$status$code == 0) est <- omx_est(f); attr(est, "logLik") <- as.numeric(logLik(f)); est } # (1) uncentered: random slope on raw ses, observed meanses at L2 e_unc <- fit3(mxModel("unc", mxModel("byClus", type = "RAM", latentVars = c("u0", "v"), data = mxData(od_b, "raw", primaryKey = "cid"), mxPath(from = c("u0", "v"), arrows = 2, free = TRUE, values = c(2.6, 0.44), labels = c("varInt", "varSlope")), mxPath(from = "u0", to = "v", arrows = 2, free = TRUE, values = -0.25, labels = "covIntSlope")), type = "RAM", manifestVars = c("y", "ses", "meanses"), data = mxData(data.frame(cid = od$cid, y = od$mAch, ses = od$ses, meanses = od$meanses), "raw", sort = FALSE), mxPath(from = "one", to = c("y", "ses", "meanses"), free = TRUE, values = c(mean(od$mAch), mean(od$ses), mean(od$meanses)), labels = c("meanY", "meanX", "meanB")), mxPath(from = c("y", "ses", "meanses"), arrows = 2, free = TRUE, values = c(36, var(od$ses), var(od$meanses)), labels = c("residY", "residX", "residB")), mxPath(from = "ses", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"), mxPath(from = "meanses", to = "y", free = TRUE, values = 3.8, labels = "contextual"), mxPath(from = c("byClus.u0", "byClus.v"), to = c("y", "y"), free = FALSE, values = c(1, NA), labels = c(NA, "data.ses"), joinKey = "cid"))) # (2) observed mean centered: random slope on cses, meanses at L2 e_obs <- fit3(mxModel("obs", mxModel("byClus", type = "RAM", latentVars = c("u0", "v"), data = mxData(od_b, "raw", primaryKey = "cid"), mxPath(from = c("u0", "v"), arrows = 2, free = TRUE, values = c(2.65, 0.67), labels = c("varInt", "varSlope")), mxPath(from = "u0", to = "v", arrows = 2, free = TRUE, values = -0.26, labels = "covIntSlope")), type = "RAM", manifestVars = c("y", "cses", "meanses"), data = mxData(data.frame(cid = od$cid, y = od$mAch, cses = od$cses, meanses = od$meanses), "raw", sort = FALSE), mxPath(from = "one", to = c("y", "cses", "meanses"), free = TRUE, values = c(mean(od$mAch), mean(od$cses), mean(od$meanses)), labels = c("meanY", "meanXw", "meanB")), mxPath(from = c("y", "cses", "meanses"), arrows = 2, free = TRUE, values = c(36, var(od$cses), var(od$meanses)), labels = c("residY", "residXw", "residB")), mxPath(from = "cses", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"), mxPath(from = "meanses", to = "y", free = TRUE, values = 5.9, labels = "betweenSlope"), mxPath(from = c("byClus.u0", "byClus.v"), to = c("y", "y"), free = FALSE, values = c(1, NA), labels = c(NA, "data.cses"), joinKey = "cid"))) # (3) Mplus hybrid: random slope on raw ses, between-level ses is latent B e_hyb <- fit3(mxModel("hyb", mxModel("byClus", type = "RAM", latentVars = c("B", "intc", "v"), data = mxData(od_b, "raw", primaryKey = "cid"), mxPath(from = c("B", "intc", "v"), arrows = 2, free = TRUE, values = c(0.16, 2.5, 0.44), labels = c("varB", "varIntR", "varSlope")), mxPath(from = "B", to = "intc", free = TRUE, values = 4.0, labels = "contextual"), mxPath(from = "B", to = "v", arrows = 2, free = TRUE, values = 0.03, labels = "covSlopeB"), mxPath(from = "intc", to = "v", arrows = 2, free = TRUE, values = -0.36, labels = "covIntcSlope")), type = "RAM", manifestVars = c("y", "x"), data = mxData(data.frame(cid = od$cid, y = od$mAch, x = od$ses), "raw", sort = FALSE), mxPath(from = "one", to = c("y", "x"), free = TRUE, values = c(mean(od$mAch), mean(od$ses)), labels = c("meanY", "meanX")), mxPath(from = c("y", "x"), arrows = 2, free = TRUE, values = c(36, var(od$ses)), labels = c("residY", "residX")), mxPath(from = "x", to = "y", free = TRUE, values = 2.2, labels = "withinSlope"), mxPath(from = c("byClus.intc", "byClus.v", "byClus.B"), to = c("y", "y", "x"), free = FALSE, values = c(1, NA, 1), labels = c(NA, "data.x", NA), joinKey = "cid"))) # ---- references ----------------------------------------------------------- fl_u <- lmer(mAch ~ ses + meanses + (ses | school), data = od, REML = FALSE) fl_o <- lmer(mAch ~ cses + meanses + (cses | school), data = od, REML = FALSE) vc_u <- as.matrix(VarCorr(fl_u)$school); vc_o <- as.matrix(VarCorr(fl_o)$school) mph <- mp_load("mplus/model_rs_hybrid.inp", "mplus/model_rs_hybrid.out") mg <- function(hdr, prm, lvl) { pp <- mp_par(mph) pp$est[pp$paramHeader == hdr & pp$param == prm & pp$BetweenWithin == lvl] } cat("### (1) Uncentered\n") ``` </details> ::: {.cell-output .cell-output-stdout} ### (1) Uncentered ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} t1 <- round(data.frame( Within = c(e_unc[["withinSlope"]], fixef(fl_u)["ses"]), Contextual = c(e_unc[["contextual"]], fixef(fl_u)["meanses"]), T00 = c(e_unc[["varInt"]], vc_u[1, 1]), T11 = c(e_unc[["varSlope"]], vc_u[2, 2]), T01 = c(e_unc[["covIntSlope"]], vc_u[1, 2]), Resid = c(e_unc[["residY"]], sigma(fl_u)^2)), 3); rownames(t1) <- c("OpenMx", "lme4") knitr::kable(t1, row.names = TRUE) ``` </details> ::: cell-output-display Within Contextual T00 T11 T01 Resid -------- -------- ------------ ------- ------- -------- -------- OpenMx 2.19 3.778 2.649 0.437 -0.235 36.797 lme4 2.19 3.778 2.649 0.437 -0.235 36.797 ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cat("### (2) Observed mean centered\n") ``` </details> ::: {.cell-output .cell-output-stdout} ### (2) Observed mean centered ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} t2 <- round(data.frame( Within = c(e_obs[["withinSlope"]], fixef(fl_o)["cses"]), Between = c(e_obs[["betweenSlope"]], fixef(fl_o)["meanses"]), T00 = c(e_obs[["varInt"]], vc_o[1, 1]), T11 = c(e_obs[["varSlope"]], vc_o[2, 2]), T01 = c(e_obs[["covIntSlope"]], vc_o[1, 2]), Resid = c(e_obs[["residY"]], sigma(fl_o)^2)), 3); rownames(t2) <- c("OpenMx", "lme4") knitr::kable(t2, row.names = TRUE) ``` </details> ::: cell-output-display Within Between T00 T11 T01 Resid -------- -------- --------- ------- ------ -------- -------- OpenMx 2.191 5.896 2.648 0.67 -0.264 36.713 lme4 2.191 5.896 2.648 0.67 -0.264 36.713 ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} cat("### (3) Mplus hybrid\n") ``` </details> ::: {.cell-output .cell-output-stdout} ### (3) Mplus hybrid ::: `<details class="code-fold"> <summary>`{=html}Code`</summary>`{=html} ``` {.r .cell-code} t3 <- round(data.frame( Within = c(e_hyb[["withinSlope"]], mg("Means", "B1", "Between")), Between = c(e_hyb[["contextual"]], mg("MACH.ON", "SES", "Between")), varB = c(e_hyb[["varB"]], mg("Variances", "SES", "Between")), T00 = c(e_hyb[["varIntR"]], mg("Residual.Variances", "MACH", "Between")), T11 = c(e_hyb[["varSlope"]], mg("Variances", "B1", "Between")), T01 = c(e_hyb[["covIntcSlope"]] - e_hyb[["contextual"]] * e_hyb[["covSlopeB"]], mg("B1.WITH", "MACH", "Between")), Resid = c(e_hyb[["residY"]], mg("Residual.Variances", "MACH", "Within")), LogLik = c(round(as.numeric(attr(e_hyb, "logLik")), 1), round(mph$summaries$LL, 1))), 3); rownames(t3) <- c("OpenMx", "Mplus") knitr::kable(t3, row.names = TRUE) ``` </details> ::: cell-output-display Within Between varB T00 T11 T01 Resid LogLik -------- -------- --------- ------ ------- ------- -------- -------- ---------- OpenMx 2.193 4.013 0.16 2.503 0.438 -0.360 36.792 -30798.8 Mplus 2.197 4.004 0.16 2.500 0.414 -0.256 36.801 -30798.8 ::: :::::::::::: The two log-likelihoods agree to within optimizer precision (-30798.8): the OpenMx and Mplus hybrid models are the same model fit to the same data by maximum likelihood. ::::::::: {#refs .references .csl-bib-body .hanging-indent} ::: {#ref-asparouhov2018 .csl-entry} Asparouhov, Thomas, and Bengt Muthén. 2018. "Latent Variable Centering of Predictors and Mediators in Multilevel and Time-Series Models." *Structural Equation Modeling: A Multidisciplinary Journal* 26 (1): 119--42. <https://doi.org/10.1080/10705511.2018.1511375>. ::: ::: {#ref-ludtke2008 .csl-entry} Lüdtke, Oliver, Herbert W. Marsh, Alexander Robitzsch, Ulrich Trautwein, Tihomir Asparouhov, and Bengt Muthén. 2008. "The Multilevel Latent Covariate Model: A New, More Reliable Approach to Group-Level Effects in Contextual Studies." *Psychological Methods* 13 (3): 203--29. <https://doi.org/10.1037/a0012869>. ::: ::: {#ref-preacher2010 .csl-entry} Preacher, Kristopher J., Michael J. Zyphur, and Zhen Zhang. 2010. "A General Multilevel SEM Framework for Assessing Multilevel Mediation." *Psychological Methods* 15 (3): 209--33. <https://doi.org/10.1037/a0020141>. ::: ::: {#ref-raudenbush2002 .csl-entry} Raudenbush, Stephen W., and Anthony S. Bryk. 2002. *Hierarchical Linear Models: Applications and Design Strategies*. 2nd ed. Sage. ::: ::: {#ref-rockwood2020 .csl-entry} Rockwood, Nicholas J. 2020. "Maximum Likelihood Estimation of Multilevel Structural Equation Models with Random Slopes for Latent Covariates." *Psychometrika* 85 (2): 275--300. <https://doi.org/10.1007/s11336-020-09702-9>. ::: ::: {#ref-rosseel2012 .csl-entry} Rosseel, Yves. 2012. "[lavaan]{.nocase}: An R Package for Structural Equation Modeling." *Journal of Statistical Software* 48 (2): 1--36. <https://doi.org/10.18637/jss.v048.i02>. ::: :::::::::