# Structural Equation Modeling in R Tutorial 7: Full Structural Equation Modeling

Made for Jonathan Butner’s Structural Equation Modeling Class, Fall 2017, University of Utah. This handout begins by showing how to import a matrix into R. Then, we will overview how to establish a measurement model in R using the lavaan package. As we go, I’ll demonstrate how to quickly and easily plot the results of your confirmatory factor analysis using semPlot. ast, I’ll demonstrate how to do basic model comparisons using lavaan objects, which will help to inform decisions related to which model fits your data better. This syntax imports the 14 variable, 307 person dataset called EXMP8MATRIX.txt, which is a free format text file. The data are for 13 observed variables comprising three theorized latent factors (Socioeconomic Status, Communication, and Conduct), as well as biological sex, dummy coded as 0 or 1, where 0 is female and 1 is male. There is no missing data.

## Data Input

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

``````library(lavaan)  #for doing the CFA
library(dplyr)  #for subsetting data quickly if needed

# make sure to set your working directory to where the file
# Regression') you can do this in R Studio using menus by
# going to the Session menu - 'Set Working Directory' - To
# Source File Location
# note that this data is in free format so we use read.table

names(semdat) <- c("Sex", "Pa_Ed", "Ma_Ed", "Pa_Occ", "Ma_Occ",
"Income", "Accept", "Listen", "Commun", "Openness", "Patience",
"Act_Out", "Agress", "Hostile")
head(semdat)  #take a look at the beginning of our datafile
``````
``````##   Sex Pa_Ed Ma_Ed Pa_Occ Ma_Occ Income Accept Listen Commun Openness
## 1   0 -0.28  0.54   1.27   0.06   1.58   1.46   1.88   2.86     0.48
## 2   0 -0.26  0.02  -0.67   0.20  -0.15   0.84   0.62  -1.06     0.22
## 3   0 -0.61 -0.01  -0.05  -0.85  -0.04   0.85   0.11  -0.07     0.36
## 4   0 -0.67 -2.26   0.45   0.16   0.43   0.12  -0.37  -0.57    -1.93
## 5   0 -0.35 -1.05   0.22   0.57  -0.63   0.14  -0.36   0.36     1.14
## 6   0  1.85 -0.90   0.97  -0.19   0.86  -0.47   1.23   0.99     0.46
##   Patience Act_Out Agress Hostile
## 1     2.00   -2.78  -2.71   -1.85
## 2     1.21   -3.07  -2.09   -1.40
## 3     0.07   -1.19  -0.52   -1.48
## 4    -0.43    0.46  -1.22   -0.26
## 5    -0.06   -2.39  -1.63    2.26
## 6    -0.11    0.37  -0.57   -1.47
``````

## Structural Equation Modeling Using lavaan: Measurement Model

Typically the first step in structural equation modeing is to establish what’s called a “measurement model”, a model which includes all of your observed variables that are going to be represented with latent variables. Constructing a measurement model allows to determine model fit related to the latent portion of your model. This way you can more precisely know where model misfit is most prevalent in your model.

Remember that `*` fixes variables to a particular value. Here we adopt the marker variable identification approach, where we fix the loading of one indicator in each latent to 1 in order to identify the model and give each latent factor a metric. This will not change model fit, just some of the loadings in the model. This makes it so that a one unit change in the latent factor is interpreted as a unit unit change in the scale of the observed variable selected as the marker variable (the variable loading we fix to 1). Note also that variances and disturbances (errors) are estimated by default.

``````sem.model.measurement <- "SES =~ 1*Pa_Ed + Ma_Ed + Pa_Occ + Ma_Occ + Income
#make 5 indicator latent socioeconomic status (SES) factor with parents education variable as the marker
COM =~ 1*Accept + Listen + Commun + Openness + Patience
#make 5 indicator latent communication factor with accept variable as the marker
Conduct =~ 1*Act_Out + Agress + Hostile
#make 3 indicator latent conduct factor with acting out variable as the marker"

sem.fit.measurement <- sem(sem.model.measurement, data = semdat)
summary(sem.fit.measurement, fit.measures = TRUE)
``````
``````## lavaan 0.6-5 ended normally after 34 iterations
##
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         29
##
##   Number of observations                           307
##
## Model Test User Model:
##
##   Test statistic                               120.948
##   Degrees of freedom                                62
##   P-value (Chi-square)                           0.000
##
## Model Test Baseline Model:
##
##   Test statistic                              1403.563
##   Degrees of freedom                                78
##   P-value                                        0.000
##
## User Model versus Baseline Model:
##
##   Comparative Fit Index (CFI)                    0.956
##   Tucker-Lewis Index (TLI)                       0.944
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)              -5835.118
##   Loglikelihood unrestricted model (H1)      -5774.644
##
##   Akaike (AIC)                               11728.236
##   Bayesian (BIC)                             11836.314
##   Sample-size adjusted Bayesian (BIC)        11744.339
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.056
##   90 Percent confidence interval - lower         0.041
##   90 Percent confidence interval - upper         0.070
##   P-value RMSEA <= 0.05                          0.251
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.049
##
## Parameter Estimates:
##
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard errors                             Standard
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   SES =~
##     Pa_Ed             1.000
##     Ma_Ed             1.457    0.144   10.144    0.000
##     Pa_Occ            1.245    0.130    9.559    0.000
##     Ma_Occ            0.896    0.101    8.893    0.000
##     Income            0.876    0.113    7.776    0.000
##   COM =~
##     Accept            1.000
##     Listen            0.811    0.067   12.035    0.000
##     Commun            1.132    0.095   11.937    0.000
##     Openness          0.602    0.070    8.609    0.000
##     Patience          0.737    0.079    9.298    0.000
##   Conduct =~
##     Act_Out           1.000
##     Agress            0.752    0.073   10.305    0.000
##     Hostile           0.893    0.082   10.842    0.000
##
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   SES ~~
##     COM               0.467    0.070    6.645    0.000
##     Conduct          -0.223    0.061   -3.659    0.000
##   COM ~~
##     Conduct          -0.487    0.088   -5.528    0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Pa_Ed             0.837    0.077   10.856    0.000
##    .Ma_Ed             0.642    0.079    8.108    0.000
##    .Pa_Occ            0.762    0.078    9.737    0.000
##    .Ma_Occ            0.580    0.055   10.609    0.000
##    .Income            0.950    0.084   11.361    0.000
##    .Accept            0.786    0.082    9.541    0.000
##    .Listen            0.480    0.052    9.317    0.000
##    .Commun            0.976    0.103    9.449    0.000
##    .Openness          0.876    0.076   11.480    0.000
##    .Patience          1.050    0.093   11.249    0.000
##    .Act_Out           0.723    0.108    6.697    0.000
##    .Agress            0.817    0.085    9.601    0.000
##    .Hostile           0.740    0.094    7.833    0.000
##     SES               0.510    0.093    5.488    0.000
##     COM               0.960    0.136    7.037    0.000
##     Conduct           1.231    0.172    7.139    0.000
``````
``````semPaths(sem.fit.measurement, "par", edge.label.cex = 1.2, fade = FALSE)  #plot our CFA
``````

In the figure, you can see the marker variables that we fixed to one as well as how each of the other observed variables loads onto each latent factor.From this model output, we can see that our fit is pretty good (CFI/TLI around .95, RMSEA approaching .05, SRMR < .05).

## Structural Equation Modeling Using lavaan: Full Model

Now, we will estimate a full model that includes predictive components in addition to measurement components. We will still use the marker variable approach, making the measurement model portion of our model unchanged.

``````sem.model.full <- "SES =~ 1*Pa_Ed + Ma_Ed + Pa_Occ + Ma_Occ + Income #make 5 indicator latent socioeconomic status (SES) factor with parents education variable as the marker
COM =~ 1*Accept + Listen + Commun + Openness + Patience; #make 5 indicator latent communication factor with accept variable as the marker
Conduct =~ 1*Act_Out + Agress + Hostile #make 3 indicator latent conduct factor with acting out variable as the marker
COM ~ Sex + SES #Regress COM on Sex and SES
Conduct ~ Sex + COM #Regress Conduct on Sex and COM
"

sem.fit.full <- sem(sem.model.full, data = semdat)
summary(sem.fit.full, fit.measures = TRUE)  #ask for model results
``````
``````## lavaan 0.6-5 ended normally after 32 iterations
##
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         30
##
##   Number of observations                           307
##
## Model Test User Model:
##
##   Test statistic                               143.535
##   Degrees of freedom                                74
##   P-value (Chi-square)                           0.000
##
## Model Test Baseline Model:
##
##   Test statistic                              1486.134
##   Degrees of freedom                                91
##   P-value                                        0.000
##
## User Model versus Baseline Model:
##
##   Comparative Fit Index (CFI)                    0.950
##   Tucker-Lewis Index (TLI)                       0.939
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)              -5805.126
##   Loglikelihood unrestricted model (H1)      -5733.359
##
##   Akaike (AIC)                               11670.252
##   Bayesian (BIC)                             11782.058
##   Sample-size adjusted Bayesian (BIC)        11686.911
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.055
##   90 Percent confidence interval - lower         0.042
##   90 Percent confidence interval - upper         0.069
##   P-value RMSEA <= 0.05                          0.247
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.053
##
## Parameter Estimates:
##
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard errors                             Standard
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   SES =~
##     Pa_Ed             1.000
##     Ma_Ed             1.479    0.146   10.154    0.000
##     Pa_Occ            1.246    0.131    9.492    0.000
##     Ma_Occ            0.907    0.102    8.910    0.000
##     Income            0.875    0.113    7.717    0.000
##   COM =~
##     Accept            1.000
##     Listen            0.807    0.066   12.160    0.000
##     Commun            1.131    0.093   12.104    0.000
##     Openness          0.604    0.069    8.756    0.000
##     Patience          0.740    0.078    9.462    0.000
##   Conduct =~
##     Act_Out           1.000
##     Agress            0.739    0.068   10.900    0.000
##     Hostile           0.877    0.074   11.782    0.000
##
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   COM ~
##     Sex              -0.311    0.105   -2.971    0.003
##     SES               0.955    0.119    7.994    0.000
##   Conduct ~
##     Sex               0.904    0.133    6.789    0.000
##     COM              -0.485    0.076   -6.377    0.000
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Pa_Ed             0.844    0.077   10.924    0.000
##    .Ma_Ed             0.626    0.078    8.026    0.000
##    .Pa_Occ            0.772    0.078    9.865    0.000
##    .Ma_Occ            0.575    0.054   10.611    0.000
##    .Income            0.956    0.084   11.405    0.000
##    .Accept            0.788    0.082    9.587    0.000
##    .Listen            0.487    0.052    9.437    0.000
##    .Commun            0.980    0.103    9.507    0.000
##    .Openness          0.874    0.076   11.482    0.000
##    .Patience          1.046    0.093   11.250    0.000
##    .Act_Out           0.695    0.100    6.982    0.000
##    .Agress            0.826    0.083    9.936    0.000
##    .Hostile           0.753    0.089    8.450    0.000
##     SES               0.503    0.092    5.453    0.000
##    .COM               0.498    0.083    6.026    0.000
##    .Conduct           0.803    0.119    6.741    0.000
``````
``````semPaths(sem.fit.full, "par", edge.label.cex = 1.2, fade = FALSE)  #plot our CFA. you can change layout with layout = argument. see ?semPaths() for more.
``````

Next, we’ll estimated an alternative, more saturated model that has one additional parameter estimated: the path between SES and conduct, with SES predicting Conduct. Essentially, we add an additional predictor in our regression, whereas that parameter was fixed to 0 in our previous model.

``````sem.model.full.1free <- "SES =~ 1*Pa_Ed + Ma_Ed + Pa_Occ + Ma_Occ + Income #make 5 indicator latent socioeconomic status (SES) factor with parents education variable as the marker
COM =~ 1*Accept + Listen + Commun + Openness + Patience; #make 5 indicator latent communication factor with accept variable as the marker
Conduct =~ 1*Act_Out + Agress + Hostile #make 3 indicator latent conduct factor with acting out variable as the marker
COM ~ Sex + SES #Regress COM on Sex and SES
Conduct ~ Sex + COM + SES #Regress Conduct on Sex, COM, and SES
"
sem.fit.full.1free <- sem(sem.model.full.1free, data = semdat)
summary(sem.fit.full.1free, fit.measures = TRUE)
``````
``````## lavaan 0.6-5 ended normally after 34 iterations
##
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         31
##
##   Number of observations                           307
##
## Model Test User Model:
##
##   Test statistic                               142.417
##   Degrees of freedom                                73
##   P-value (Chi-square)                           0.000
##
## Model Test Baseline Model:
##
##   Test statistic                              1486.134
##   Degrees of freedom                                91
##   P-value                                        0.000
##
## User Model versus Baseline Model:
##
##   Comparative Fit Index (CFI)                    0.950
##   Tucker-Lewis Index (TLI)                       0.938
##
## Loglikelihood and Information Criteria:
##
##   Loglikelihood user model (H0)              -5804.567
##   Loglikelihood unrestricted model (H1)      -5733.359
##
##   Akaike (AIC)                               11671.134
##   Bayesian (BIC)                             11786.666
##   Sample-size adjusted Bayesian (BIC)        11688.348
##
## Root Mean Square Error of Approximation:
##
##   RMSEA                                          0.056
##   90 Percent confidence interval - lower         0.042
##   90 Percent confidence interval - upper         0.069
##   P-value RMSEA <= 0.05                          0.236
##
## Standardized Root Mean Square Residual:
##
##   SRMR                                           0.054
##
## Parameter Estimates:
##
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard errors                             Standard
##
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   SES =~
##     Pa_Ed             1.000
##     Ma_Ed             1.480    0.146   10.162    0.000
##     Pa_Occ            1.242    0.131    9.475    0.000
##     Ma_Occ            0.909    0.102    8.924    0.000
##     Income            0.876    0.113    7.720    0.000
##   COM =~
##     Accept            1.000
##     Listen            0.810    0.066   12.221    0.000
##     Commun            1.130    0.093   12.109    0.000
##     Openness          0.600    0.069    8.709    0.000
##     Patience          0.738    0.078    9.448    0.000
##   Conduct =~
##     Act_Out           1.000
##     Agress            0.741    0.068   10.970    0.000
##     Hostile           0.875    0.074   11.822    0.000
##
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   COM ~
##     Sex              -0.309    0.105   -2.945    0.003
##     SES               0.948    0.119    7.944    0.000
##   Conduct ~
##     Sex               0.945    0.136    6.969    0.000
##     COM              -0.390    0.113   -3.459    0.001
##     SES              -0.165    0.153   -1.079    0.281
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Pa_Ed             0.844    0.077   10.928    0.000
##    .Ma_Ed             0.623    0.078    8.009    0.000
##    .Pa_Occ            0.777    0.078    9.902    0.000
##    .Ma_Occ            0.573    0.054   10.602    0.000
##    .Income            0.956    0.084   11.407    0.000
##    .Accept            0.785    0.082    9.544    0.000
##    .Listen            0.480    0.051    9.330    0.000
##    .Commun            0.979    0.103    9.474    0.000
##    .Openness          0.877    0.076   11.489    0.000
##    .Patience          1.048    0.093   11.247    0.000
##    .Act_Out           0.694    0.099    6.986    0.000
##    .Agress            0.822    0.083    9.915    0.000
##    .Hostile           0.757    0.089    8.508    0.000
##     SES               0.503    0.092    5.453    0.000
##    .COM               0.507    0.084    6.041    0.000
##    .Conduct           0.808    0.119    6.793    0.000
``````
``````semPaths(sem.fit.full.1free, "par", edge.label.cex = 1.2, fade = FALSE)  #plot our CFA
``````

Note how the effect of SES on Conduct is fairly weak (-.16 unstandardized, not significantly different from 0 ,p = 0.28). In addition, fit appears to have worsened (CFI/TLI are smaller, RMSEA/SRMR are larger), but we can explicitly quantitatively test this using the Chi-squared difference of fit test since the models are nested (more on that later), which we will do in the next section.

## 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 (regression between Conduct and SES fixed to 0) to the more saturated model with one less df (regression between Conduct and SES estimated).

``````anova(sem.fit.full, sem.fit.full.1free)  #model fit comparison
``````
``````## Chi-Squared Difference Test
##
##                    Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## sem.fit.full.1free 73 11671 11787 142.42
## sem.fit.full       74 11670 11782 143.54     1.1182       1     0.2903
``````

Since we retain the null hypothesis using the chi-squared difference of fit test, this means that the less saturated model with one less parameter estimated (in this case, when the regression between Conduct and SES is fixed to 0) fits the data better than the more free model with that regression coefficient estimated. Specifically, the previous model chi-square was 139.486, DF=73, the change between the two is 139.486-138.345=1.141, DF=73-72=1, cutoff value in a chi-squared distribution for 1 df is 3.84 at p=.05. we retain the null thus the new model wth an additional parameter estimated is not a significant improvement over the first more restricted model where that parameter was fixed to 0.