Structural Equation Modeling using lavaan in R: Introduction to Stacked Models and Factorial Measurement Invariance Testing

Table of Contents

  1. Read in datafile
  2. Stacked Models in lavaan
  3. Factorial Measurement Invariance Testing in lavaan

Made for Jonathan Butner’s Structural Equation Modeling Class, Fall 2017, University of Utah. This handout will focus on implementing stacked models in lavaan, which allow us to test a model for two different groups (for example, control vs. intervention). This syntax imports the X variable, 192 person dataset called HW9 2017.csv, which is an comma separated data file. The data are for 10 items on the SF36 scale collected at 2 timepoints (20 total observed items for the scale), as well as a condition variable where 0 is control and 1 is intervention. There is also a diagnosis variable with 3 diagnoses. There is missing data coded as 99 in the original datafile. See http://lavaan.ugent.be/tutorial/groups.html for more on stacked models in lavaan.

Data Input

First, we will read in the datafile and clean up the datafile a little bit, as well as load required packages.

require(lavaan)  #for doing the SEM
require(dplyr)  #for subsetting data quickly if needed

# make sure to set your working directory to where the file
# is located, e.g. setwd('D:/Downloads/Handout 1 -
# Regression') you can do this in R Studio using menus by
# going to the Session menu - 'Set Working Directory' - To
# Source File Location
semdat <- read.csv("HW9 2017.csv", na.strings = 99)
# note that this data is in free format so we use read.table

# fix first column name
names(semdat)[1] <- c("condition")

# create dummy variables for diagnosis (for demonstration - i
# won't show any models with diagnosis)
semdat$d1[semdat$diag == 2] <- 0
semdat$d1[semdat$diag == 1] <- 1
semdat$d2[semdat$diag == 2] <- 1
semdat$d2[semdat$diag == 1] <- 0

head(semdat)  #take a look at the beginning of our datafile
##   condition diag sf3 sf4 sf5 sf6 sf7 sf8 sf9 sf10 sf11 sf12 sf3_1 sf4_1
## 1         1    2   1   1   1   1   1   1   1    1    1    2     1     1
## 2         1    2   2   2   2   2   2   3   1    1    2    3     2     2
## 3         1    3   1   1   1   1   1   1   1    1    1    2     1     2
## 4         0    2   1   1   1   1   1   2   1    2    3    2     1     1
## 5         1    2   1   2   2   1   1   2   1    2    3    3     1     1
## 6         1    3   1   1   1   1   1   1   1    1    1    1    NA    NA
##   sf5_1 sf6_1 sf7_1 sf8_1 sf9_1 sf10_1 sf11_1 sf12_1 d1 d2
## 1     1     1     1     1     1      1      1      1  0  1
## 2     2     2     2     1     3      3      2      2  0  1
## 3     1     2     3     2     1      1      1      3 NA NA
## 4     4     4     4     3     4      3      3      3  0  1
## 5     2     1     2     3     1      2      3      3  0  1
## 6    NA    NA    NA    NA    NA     NA     NA     NA NA NA

Stacked Models in Lavaan

Your group variable should be specified using the group = 'variablename' argument. Notice that lavaan, unlike Mplus, does not automatically assume equality of lambdas (regression parameters) and taus (intercepts) across groups when generating stacked models. Note that here we add an additional argument called group.equal = c("loadings"), which equates lambdas across groups. Think of this as assuming that there is no interaction with group (effect is same across groups). You can also extend this to include group.equal = c('loadings', 'intercepts', 'residuals'), amongst other arguments depending on which equality constraints you want to include in your model. Note that we can compare models where we have equal parameters across groups to a model where the parameters are different across groups. Note you can do the same with intercepts. YOu would then do a chi-sq difference test where you compare a model with equality constrants to a model where. This concept is called factorial measurement invariance testing (see the next section for how to do this). The bit of code specifying missing = "fiml" implements full information maximum likelihood estimation to handle missing data, which assumes that your data is missing at random (MAR) or missing completely at random (MCAR).

sem.model.equality <- "physfun =~ sf3+sf4+sf5+sf6+sf7+sf8+sf9+sf10+sf11+sf12"
sem.fit.equality <- sem(sem.model.equality, data = semdat, group = "condition",
    group.equal = c("loadings"), missing = "fiml")  #missing = fiml implements full
# information ML for missingness
summary(sem.fit.equality, fit.measures = TRUE)
## lavaan (0.5-23.1097) converged normally after  73 iterations
##
##   Number of observations per group
##   1                                                101
##   0                                                 91
##
##   Number of missing patterns per group
##   1                                                  1
##   0                                                  1
##
##   Estimator                                         ML
##   Minimum Function Test Statistic              295.364
##   Degrees of freedom                                79
##   P-value (Chi-square)                           0.000
##
## Chi-square for each group:
##
##   1                                            186.914
##   0                                            108.450
##
## Model test baseline model:
##
##   Minimum Function Test Statistic              802.850
##   Degrees of freedom                                90
##   P-value                                        0.000
##
## User model versus baseline model:
##
##   Comparative Fit Index (CFI)                    0.696
##   Tucker-Lewis Index (TLI)                       0.654
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)              -1387.924
##   Loglikelihood unrestricted model (H1)      -1240.242
##
##   Number of free parameters                         51
##   Akaike (AIC)                                2877.848
##   Bayesian (BIC)                              3043.980
##   Sample-size adjusted Bayesian (BIC)         2882.428
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.169
##   90 Percent Confidence Interval          0.149  0.190
##   P-value RMSEA <= 0.05                          0.000
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.099
##
## Parameter Estimates:
##
##   Information                                 Observed
##   Standard Errors                             Standard
##
##
## Group 1 [1]:
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   physfun =~
##     sf3               1.000
##     sf4     (.p2.)    5.427    1.377    3.941    0.000
##     sf5     (.p3.)    5.955    1.530    3.893    0.000
##     sf6     (.p4.)    3.330    0.952    3.498    0.000
##     sf7     (.p5.)    6.114    1.612    3.793    0.000
##     sf8     (.p6.)    4.679    1.305    3.587    0.000
##     sf9     (.p7.)    3.383    0.916    3.695    0.000
##     sf10    (.p8.)    4.235    1.127    3.758    0.000
##     sf11    (.p9.)    6.192    1.626    3.808    0.000
##     sf12    (.10.)    4.373    1.241    3.523    0.000
##
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               1.050    0.032   32.428    0.000
##    .sf4               1.366    0.067   20.489    0.000
##    .sf5               1.515    0.072   21.010    0.000
##    .sf6               1.149    0.055   21.019    0.000
##    .sf7               1.416    0.076   18.544    0.000
##    .sf8               1.663    0.076   21.822    0.000
##    .sf9               1.109    0.048   22.897    0.000
##    .sf10              1.238    0.053   23.296    0.000
##    .sf11              1.604    0.072   22.383    0.000
##    .sf12              2.327    0.073   32.077    0.000
##     physfun           0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               0.098    0.014    7.015    0.000
##    .sf4               0.223    0.037    6.073    0.000
##    .sf5               0.252    0.042    6.054    0.000
##    .sf6               0.216    0.034    6.453    0.000
##    .sf7               0.301    0.050    6.034    0.000
##    .sf8               0.418    0.062    6.746    0.000
##    .sf9               0.149    0.023    6.380    0.000
##    .sf10              0.147    0.025    6.001    0.000
##    .sf11              0.224    0.040    5.589    0.000
##    .sf12              0.384    0.057    6.696    0.000
##     physfun           0.008    0.004    1.936    0.053
##
##
## Group 2 [0]:
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   physfun =~
##     sf3               1.000
##     sf4     (.p2.)    5.427    1.377    3.941    0.000
##     sf5     (.p3.)    5.955    1.530    3.893    0.000
##     sf6     (.p4.)    3.330    0.952    3.498    0.000
##     sf7     (.p5.)    6.114    1.612    3.793    0.000
##     sf8     (.p6.)    4.679    1.305    3.587    0.000
##     sf9     (.p7.)    3.383    0.916    3.695    0.000
##     sf10    (.p8.)    4.235    1.127    3.758    0.000
##     sf11    (.p9.)    6.192    1.626    3.808    0.000
##     sf12    (.10.)    4.373    1.241    3.523    0.000
##
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               1.022    0.022   46.600    0.000
##    .sf4               1.297    0.056   23.220    0.000
##    .sf5               1.396    0.064   21.672    0.000
##    .sf6               1.132    0.055   20.714    0.000
##    .sf7               1.440    0.079   18.298    0.000
##    .sf8               1.747    0.072   24.220    0.000
##    .sf9               1.088    0.043   25.389    0.000
##    .sf10              1.165    0.049   23.804    0.000
##    .sf11              1.527    0.074   20.704    0.000
##    .sf12              2.396    0.073   32.749    0.000
##     physfun           0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               0.039    0.006    6.591    0.000
##    .sf4               0.135    0.026    5.266    0.000
##    .sf5               0.198    0.036    5.530    0.000
##    .sf6               0.216    0.034    6.426    0.000
##    .sf7               0.374    0.062    6.079    0.000
##    .sf8               0.363    0.058    6.275    0.000
##    .sf9               0.109    0.018    6.081    0.000
##    .sf10              0.127    0.022    5.825    0.000
##    .sf11              0.302    0.051    5.944    0.000
##    .sf12              0.390    0.060    6.458    0.000
##     physfun           0.005    0.003    1.940    0.052

You don’t have to constrain everything at once, either, if you have a theoretical reason not to. You can constrain individaul parameters across groups if you don’t want to constrain all loadings, for example. I do this here, where I constrain just certain indicators to be the same across groups. Though I don’t show it in this example, you can also fix parameters in one group, but not the other group through syntax such as c(1,NA)*sf4. This will fix the loading of the sf4 item to 1 in the first group but not the second group, where NA tells lavaan to freely estimate the parameter for that group. Note that this syntax for fixing and equating variables can be extended to more than two groups simply by increasing the number of items in the list. For example, if you wanted to equate parameters across three groups in the example below you would use c(v6,v6,v6)*sf6. This syntax can be extended to any parameters in your model that you want to fix or equate, not just loadings of observed variables on latent factors.

sem.model.partequality <- "physfun =~ c(v3,v3)*sf3+sf4+sf5+c(v6,v6)*sf6+sf7+sf8+sf9+sf10+sf11+sf12"
sem.fit.partequality <- sem(sem.model.partequality, data = semdat,
    group = "condition", missing = "fiml")  #missing = fiml implements full
# information ML for missingness
summary(sem.fit.partequality, fit.measures = TRUE)
## lavaan (0.5-23.1097) converged normally after 109 iterations
##
##   Number of observations per group
##   1                                                101
##   0                                                 91
##
##   Number of missing patterns per group
##   1                                                  1
##   0                                                  1
##
##   Estimator                                         ML
##   Minimum Function Test Statistic              288.803
##   Degrees of freedom                                71
##   P-value (Chi-square)                           0.000
##
## Chi-square for each group:
##
##   1                                            184.417
##   0                                            104.386
##
## Model test baseline model:
##
##   Minimum Function Test Statistic              802.850
##   Degrees of freedom                                90
##   P-value                                        0.000
##
## User model versus baseline model:
##
##   Comparative Fit Index (CFI)                    0.694
##   Tucker-Lewis Index (TLI)                       0.613
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)              -1384.643
##   Loglikelihood unrestricted model (H1)      -1240.242
##
##   Number of free parameters                         59
##   Akaike (AIC)                                2887.286
##   Bayesian (BIC)                              3079.479
##   Sample-size adjusted Bayesian (BIC)         2892.585
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.179
##   90 Percent Confidence Interval          0.158  0.200
##   P-value RMSEA <= 0.05                          0.000
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.091
##
## Parameter Estimates:
##
##   Information                                 Observed
##   Standard Errors                             Standard
##
##
## Group 1 [1]:
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   physfun =~
##     sf3       (v3)    1.000
##     sf4               4.443    1.276    3.481    0.000
##     sf5               4.932    1.427    3.456    0.001
##     sf6       (v6)    3.250    0.936    3.472    0.001
##     sf7               5.238    1.518    3.450    0.001
##     sf8               3.396    1.130    3.005    0.003
##     sf9               2.955    0.909    3.250    0.001
##     sf10              3.905    1.158    3.372    0.001
##     sf11              5.534    1.600    3.458    0.001
##     sf12              3.840    1.232    3.116    0.002
##
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               1.050    0.033   31.969    0.000
##    .sf4               1.366    0.065   20.937    0.000
##    .sf5               1.515    0.071   21.377    0.000
##    .sf6               1.149    0.056   20.425    0.000
##    .sf7               1.416    0.076   18.691    0.000
##    .sf8               1.663    0.073   22.745    0.000
##    .sf9               1.109    0.048   22.957    0.000
##    .sf10              1.238    0.054   22.710    0.000
##    .sf11              1.604    0.073   22.013    0.000
##    .sf12              2.327    0.073   31.915    0.000
##     physfun           0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               0.099    0.014    6.994    0.000
##    .sf4               0.228    0.038    6.067    0.000
##    .sf5               0.258    0.042    6.070    0.000
##    .sf6               0.211    0.034    6.296    0.000
##    .sf7               0.298    0.050    5.914    0.000
##    .sf8               0.422    0.062    6.821    0.000
##    .sf9               0.146    0.023    6.242    0.000
##    .sf10              0.143    0.025    5.835    0.000
##    .sf11              0.222    0.041    5.392    0.000
##    .sf12              0.385    0.058    6.664    0.000
##     physfun           0.010    0.006    1.796    0.073
##
##
## Group 2 [0]:
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   physfun =~
##     sf3       (v3)    1.000
##     sf4               6.849    2.161    3.169    0.002
##     sf5               7.606    2.427    3.134    0.002
##     sf6       (v6)    3.250    0.936    3.472    0.001
##     sf7               7.416    2.473    2.998    0.003
##     sf8               6.880    2.353    2.924    0.003
##     sf9               3.915    1.340    2.922    0.003
##     sf10              4.409    1.524    2.893    0.004
##     sf11              6.604    2.317    2.850    0.004
##     sf12              5.034    1.922    2.619    0.009
##
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               1.022    0.022   47.330    0.000
##    .sf4               1.297    0.057   22.714    0.000
##    .sf5               1.396    0.066   21.238    0.000
##    .sf6               1.132    0.053   21.504    0.000
##    .sf7               1.440    0.080   18.094    0.000
##    .sf8               1.747    0.075   23.147    0.000
##    .sf9               1.088    0.043   25.313    0.000
##    .sf10              1.165    0.047   24.626    0.000
##    .sf11              1.527    0.072   21.276    0.000
##    .sf12              2.396    0.073   32.956    0.000
##     physfun           0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .sf3               0.039    0.006    6.606    0.000
##    .sf4               0.127    0.026    4.903    0.000
##    .sf5               0.184    0.036    5.165    0.000
##    .sf6               0.214    0.033    6.507    0.000
##    .sf7               0.377    0.063    5.974    0.000
##    .sf8               0.347    0.057    6.100    0.000
##    .sf9               0.113    0.019    6.056    0.000
##    .sf10              0.133    0.023    5.900    0.000
##    .sf11              0.311    0.052    5.978    0.000
##    .sf12              0.389    0.060    6.441    0.000
##     physfun           0.004    0.002    1.624    0.104

Model Comparison Using lavaan

Note that models that are compared using most fit statistics must be nested in order for the tests to be valid. Nested models are models that contain at least all of the same exact observed variables contained in the less complicated model. For nested models, you can only free or fix parameter(s), not do both. The underlying data matrix must be the exact same between models. The code below compares the reduced model with more df (paths between physical functioning and observed indicators equated to be the same) to the more saturated model with less df (paths between physical functioning and indicators allowed to differ, and thus estimated differently for both conditions). Note that you can compare models directly in order to have an explicit test of factorial measurement invariance. If we retain the null hypothesis, our less saturated model with more parameters fixed is a better fit ot the data, and we can claim that the loadings are equivalent across groups (we have achieved one level of factorial measurement invariance).

First, we’ll fit our model where we allow the lambdas to be different across conditions, which is the default in lavaan.

sem.model.diff <- "physfun =~ sf3+sf4+sf5+sf6+sf7+sf8+sf9+sf10+sf11+sf12"
sem.fit.diff <- sem(sem.model.diff, data = semdat, group = "condition",
    missing = "fiml")
# missing = fiml implements full information ML for
# missingness

Now, we’ll compare the models (the model we just made to the model earlier where we equated lambdas - the paths from our latent variables to our observed variables - to be the same across groups) using a chi-squared difference of fit test. This is considered a test of what’s typically referred to as weak invariance.

anova(sem.fit.equality, sem.fit.diff)  #model fit comparison
## Chi Square Difference Test
##
##                  Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## sem.fit.diff     70 2888.7 3084.2 288.26
## sem.fit.equality 79 2877.8 3044.0 295.36     7.1058       9     0.6261

Since we retain the null hypothesis using the chi-squared difference of fit test, this means that the less saturated model with less parameters estimated (in this case, when the groups have factor loadings equated) fits the data better than the more free model with those factor loadings allowed to differ across groups. This means that we have met weak invariance. You could now keep building on this model, next testing whether intercepts are the same across groups (all the way up to means and residuals constrained to the same across groups, etc.), which would be a test of strong invariance. Note in these models you would keep the factor loadings constrained to be the same since you already know this reflects your data more accurately.