model <- '
level: 1
fw =~ y1 + y2 + y3
fw ~ x1 + x2 + x3
level: 2
fb =~ y1 + y2 + y3
fb ~ w1 + w2
'Multilevel SEM
If the data is clustered, one way to handle the clustering is to use a multilevel modeling approach. In the SEM framework, this leads to multilevel SEM. The multilevel capabilities of lavaan are still limited, but you can fit a two-level SEM with random intercepts (and, since version 0.7, also random slopes; see the last section below).
Multilevel SEM model syntax
To fit a two-level SEM, you must specify a model for both levels, as follows:
This model syntax contains two blocks, one for level 1, and one for level 2. Within each block, you can specify a model just like in the single-level case. To fit this model, using a toy dataset Demo.twolevel that is part of the lavaan package, you need to add the cluster= argument to the sem/lavaan function call:
fit <- sem(model = model, data = Demo.twolevel, cluster = "cluster")The output looks similar to a multigroup SEM output, but where the two groups are now the within and the between level respectively.
summary(fit)lavaan 0.7-2 ended normally after 36 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 20
Number of observations 2500
Number of clusters [cluster] 200
Model Test User Model:
Test statistic 8.092
Degrees of freedom 10
P-value (Chi-square) 0.620
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Level 1 [within]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fw =~
y1 1.000
y2 0.774 0.034 22.671 0.000
y3 0.734 0.033 22.355 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fw ~
x1 0.510 0.023 22.037 0.000
x2 0.407 0.022 18.273 0.000
x3 0.205 0.021 9.740 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.986 0.046 21.591 0.000
.y2 1.066 0.039 27.271 0.000
.y3 1.011 0.037 27.662 0.000
.fw 0.546 0.040 13.539 0.000
Level 2 [cluster]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fb =~
y1 1.000
y2 0.717 0.052 13.824 0.000
y3 0.587 0.048 12.329 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fb ~
w1 0.165 0.079 2.093 0.036
w2 0.131 0.076 1.715 0.086
Intercepts:
Estimate Std.Err z-value P(>|z|)
.y1 0.024 0.075 0.327 0.743
.y2 -0.016 0.060 -0.269 0.788
.y3 -0.042 0.054 -0.777 0.437
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.058 0.047 1.213 0.225
.y2 0.120 0.031 3.825 0.000
.y3 0.149 0.028 5.319 0.000
.fb 0.899 0.118 7.592 0.000After fitting the model, you can inspect the intra-class correlations:
lavInspect(fit, "icc") y1 y2 y3 x1 x2 x3
0.331 0.263 0.232 0.000 0.000 0.000 The see the unrestricted (h1) within and between means and covariances, you can use
lavInspect(fit, "h1")$within
$within$cov
y1 y2 y3 x1 x2 x3
y1 2.000
y2 0.789 1.674
y3 0.749 0.564 1.557
x1 0.489 0.393 0.376 0.982
x2 0.416 0.322 0.299 0.001 1.011
x3 0.221 0.160 0.155 -0.006 0.008 1.045
$within$mean
y1 y2 y3 x1 x2 x3
0.000 0.000 0.000 -0.007 -0.003 0.020
$cluster
$cluster$cov
y1 y2 y3 w1 w2
y1 0.992
y2 0.668 0.598
y3 0.548 0.391 0.469
w1 0.125 0.119 0.036 0.870
w2 0.086 0.057 0.130 -0.128 0.931
$cluster$mean
y1 y2 y3 w1 w2
0.020 -0.019 -0.045 0.052 -0.091 Important notes
note that in
level: 1the colon follows thelevelkeyword; if you typelevel 1:, you will get an erroryou must specify a model for each level; the following syntax is not allowed and will produce an error:
model <- ' level: 1 fw =~ y1 + y2 + y3 fw ~ x1 + x2 + x3 level: 2 'if you do not have a model in mind for level 2, you can specify a saturated level by adding all variances and covariances of the endogenous variables (here: y1, y2 and y3):
model <- ' level: 1 fw =~ y1 + y2 + y3 fw ~ x1 + x2 + x3 level: 2 y1 ~~ y1 + y2 + y3 y2 ~~ y2 + y3 y3 ~~ y3 '
Convergence issues and solutions
By default, the current version of lavaan (0.7) uses a quasi-Newton procedure to maximize the loglikelihood of the data given the model (just like in the single-level case). For most model and data combinations, this will work fine (and fast). However, every now and then, you may experience convergence issues.
Non-convergence is typically a sign that something is not quite right with either your model, or your data. Typical settings are: a small number of clusters, in combination with (almost) no variance of an endogenous variable at the between level.
However, if you believe nothing is wrong, you may want to try another optimization procedure. The current version of lavaan allows for using the Expectation Maximization (EM) algorithm as an alternative. To switch to the EM algorithm, you can use:
fit <- sem(model = model, data = Demo.twolevel, cluster = "cluster",
verbose = TRUE, optim.method = "em")Since version 0.7, the EM algorithm is accelerated (using SQUAREM by default), so it usually reaches a solution in far fewer iterations than before. Still, on hard model/data combinations, many iterations may be needed. The settings of the EM algorithm are controlled by the em.args argument, which takes a list of options. To change the stopping criteria, you can use:
fit <- sem(model = model, data = Demo.twolevel, cluster = "cluster",
verbose = TRUE, optim.method = "em",
em.args = list(max_iter = 20000,
fx_tol = 1e-08, dx_tol = 1e-04))The fx_tol element is used to monitor the change in loglikelihood between the current step and the previous step. If this change is smaller than fx_tol, the algorithm stops. The dx_tol element is used to monitor the (unscaled) gradient. When a solution is reached, all elements of the gradient should be near zero. When the largest gradient element is smaller than dx_tol, the algorithm stops. The max_iter element sets the maximum number of EM iterations.
(Note: in older versions of lavaan, these settings were provided as separate em.iter.max, em.fx.tol and em.dx.tol arguments; they are now collected in the em.args list. A related em.h1.args list controls the EM algorithm used for the unrestricted (h1) model.)
A word of caution: the EM algorithm can always be forced to ‘converge’ (perhaps after changing the stopping criteria), but that does not mean you have a model/dataset combination that deserves to converge.
Random slopes
Up to now, only the intercepts were allowed to vary across clusters (random intercepts). Since version 0.7, lavaan also supports random slopes: regression coefficients at the within level that are allowed to vary across clusters. To turn a within-level regression coefficient into a random slope, you wrap the corresponding predictor in the rv() modifier. The (quoted) argument of rv() is the name of a latent variable that lives at the between level, and that represents the cluster-specific slope. At level 2, this new latent variable can be given a variance (the slope variance), an intercept (the average or ‘fixed’ effect of the slope), covariances with other between-level variables, and it can be regressed on between-level covariates.
In the example below (using the Demo.twolevel dataset again), we let the within-level regressions of the latent variable fw on the covariates x1 and x2 vary across clusters. The random slopes are named s1 and s2, and are regressed on the between-level covariates w1 and w2:
model <- '
level: 1
fw =~ y1 + y2 + y3
fw ~ rv("s1")*x1 + rv("s2")*x2 + x3
level: 2
fb =~ y1 + y2 + y3
fb ~ w1 + w2
# the random slopes are latent variables at the between level
s1 + s2 ~ w1 + w2
'
fit <- sem(model = model, data = Demo.twolevel, cluster = "cluster")
summary(fit)lavaan 0.7-2 ended normally after 55 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 26
Number of observations 2500
Number of clusters [cluster] 200
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Level 1 [within]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fw =~
y1 1.000
y2 0.773 0.034 22.643 0.000
y3 0.733 0.033 22.353 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fw ~
x1 0.000
x2 0.000
x3 0.204 0.021 9.706 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.985 0.046 21.555 0.000
.y2 1.067 0.039 27.274 0.000
.y3 1.011 0.037 27.676 0.000
.fw 0.540 0.041 13.041 0.000
Level 2 [cluster]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fb =~
y1 1.000
y2 0.717 0.052 13.812 0.000
y3 0.587 0.048 12.318 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fb ~
w1 0.163 0.079 2.074 0.038
w2 0.130 0.076 1.700 0.089
s1 ~
w1 0.010 0.022 0.448 0.654
w2 0.009 0.024 0.376 0.707
s2 ~
w1 -0.041 0.022 -1.846 0.065
w2 -0.025 0.023 -1.126 0.260
Intercepts:
Estimate Std.Err z-value P(>|z|)
.y1 0.025 0.075 0.332 0.740
.y2 -0.016 0.060 -0.264 0.791
.y3 -0.042 0.054 -0.772 0.440
.s1 0.511 0.024 21.609 0.000
.s2 0.406 0.023 17.895 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.058 0.047 1.213 0.225
.y2 0.120 0.031 3.824 0.000
.y3 0.149 0.028 5.318 0.000
.fb 0.899 0.118 7.588 0.000
.s1 0.003 0.008 0.363 0.716
.s2 0.002 0.008 0.267 0.790In the output, the average (fixed) effect of each random slope shows up as the intercept of s1 and s2 at the between level, while their (residual) variances describe how strongly the slopes vary from cluster to cluster.
The predictor that carries a random slope may also be a latent within-level variable, for example:
model <- '
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
'In this case, the random slope multiplies two random quantities, and the marginal loglikelihood is no longer available in closed form; it is then computed by (Gauss-Hermite) quadrature over the latent-covariate slopes (the number of quadrature points per dimension can be set with the integration.ngh option).
The current implementation of random slopes is limited to two-level models, a single group, continuous data, and estimator = "ML"; observed covariates that carry a random slope must be within-only (level-1) variables. Missing data can be handled with missing = "ml". Note that no chi-square test statistic (or the fit indices derived from it) is available for models with random slopes; fitMeasures() will only report the loglikelihood-based measures (such as AIC and BIC).