Using random slopes with lavaan

Beginning with v0.7, lavaan can now fit multilevel models with random slopes (aside from that, lavaan can now estimate multilevel categorical CFA models too!).

1. Make a toy dataset

Here’s an example dataset with a random slope:

library(MASS)
set.seed(456)
J <- 100
nJ <- 30
tot <- J * nJ
gl <- rep(1:J, each = nJ)
x1 <- rnorm(tot)
x2 <- rnorm(tot)
xx <- mvrnorm(J, 
        Sigma = matrix(c(.5, .05, .05, .10), 2),
        mu = c(0, 0))
u2 <- rep(xx[, 1], each = nJ)
u2s <- rep(xx[, 2], each = nJ)
u1 <- rnorm(tot)
y <- .75 * x1 + .5 * x2 + u2 + u1 + x1 * u2s
dd = data.frame(y, gl, x1, x2)
dd$one <- 1

2. Fit a model

Test first if the random slope is warranted (just using a likelihood ratio test):

library(lme4)
t0 <- lmer(y ~ x1 + x2 + (1|gl), REML = F)
t1 <- lmer(y ~ x1 + x2 + (x1|gl), REML = F) 
anova(t1, t0) #random slope is warranted

Data: NULL
Models:
t0: y ~ x1 + x2 + (1 | gl)
t1: y ~ x1 + x2 + (x1 | gl)
   npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
t0    5 9118.1 9148.2 -4554.1    9108.1                         
t1    7 8980.2 9022.2 -4483.1    8966.2 141.93  2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(t1) #result will differ since MBSE ne robust SE

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: y ~ x1 + x2 + (x1 | gl)

      AIC       BIC    logLik -2*log(L)  df.resid 
   8980.2    9022.2   -4483.1    8966.2      2993 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0774 -0.6407 -0.0063  0.6476  3.1642 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 gl       (Intercept) 0.51117  0.7150        
          x1          0.09659  0.3108   0.28 
 Residual             1.01730  1.0086        
Number of obs: 3000, groups:  gl, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.03411    0.07389  -0.462
x1           0.79348    0.03665  21.648
x2           0.50553    0.01888  26.774

Correlation of Fixed Effects:
   (Intr) x1    
x1  0.223       
x2 -0.009 -0.005

Fit a model using maximum likelihood with robust standard errors:

library(WeMix)
d1 <- mix(y ~ x1 + x2 + (x1|gl),
          data = dd,
          weights = c('one', 'one'))
summary(d1) #this will be the same

Call:
mix(formula = y ~ x1 + x2 + (x1 | gl), data = dd, weights = c("one", 
    "one"))

Variance terms:
 Level    Group        Name Variance Std. Error Std.Dev. Corr1
     2       gl (Intercept)  0.51117    0.07274   0.7150      
     2       gl          x1  0.09659    0.02219   0.3108  0.28
     1 Residual              1.01730    0.02582   1.0086      
Groups:
 Level Group n size mean wgt sum wgt
     2    gl    100        1     100
     1   Obs   3000        1    3000

Fixed Effects:
            Estimate Std. Error t value
(Intercept) -0.03411    0.07427  -0.459
x1           0.79348    0.03683  21.542
x2           0.50553    0.01615  31.293

lnl= -4483.10 
Intraclass Correlation= 0.374 

We can compare this to the following model fit using lavaan.

3. Fit using lavaan

NOTE how the syntax is specified. Looks a little different.

library(lavaan)
model <- '
level: 1
  y ~ rv("sx1") * x1 + x2

level:2
  sx1 ~~ sx1 #variance for T11
  y ~~ sx1 #covariance T11 and T00
  y ~~ y #T00

  y ~ 1
  sx1 ~ 1
'

t3 <- sem(model, data = dd, cluster = 'gl')

To get the output:

summary(t3, standardized = TRUE)

lavaan 0.7-2 ended normally after 34 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         7

  Number of observations                          3000
  Number of clusters [gl]                          100


Parameter Estimates:

  Standard errors                             Standard
  Information                                 Observed
  Observed information based on                Hessian


Level 1 [within]:

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  y ~                                                                   
    x1                0.000                               0.000    0.000
    x2                0.506    0.019   26.752    0.000    0.506    0.449

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .y                 1.017    0.027   37.428    0.000    1.017    0.799


Level 2 [gl]:

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  sx1 ~~                                                                
   .y                 0.061    0.028    2.210    0.027    0.198    0.277

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .y                -0.034    0.074   -0.462    0.644   -0.034   -0.048
    sx1               0.793    0.037   21.646    0.000    2.553    2.553

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    sx1               0.097    0.019    5.154    0.000    1.000    1.000
   .y                 0.511    0.077    6.620    0.000    0.511    1.000

NOTE where to find the slope effects and the random variances. Consult the previous output to compare. Might not be in the place you expect. The fixed effect coefficient for x1 is at level 2 (B = 0.793, SE = 0.037).

What about fitting a simpler model?

I show in another manuscript that fitting a random intercept model with robust standard errors can be an acceptable alternative to fitting a more complicated random slope model (only if you are not interested in the random effects).

The standard errors looks similar to the RS model when using the RI model with robust standard errors.

mod_simp <- '
level: 1
  y ~ x1 + x2

level: 2
  y ~ 1'

t4a <- sem(mod_simp, data = dd, cluster = 'gl',
           estimator = "ML") #not robust
t4b <- sem(mod_simp, data = dd, cluster = 'gl',
          estimator = "MLR") #robust

library(modelsummary)
modelsummary(list("RI/not rob" = t4a,
                  "RI /rob" = t4b),
             coef_omit = '~~', 
             stars = TRUE)

— END