SEM Fit and Modification
PSY 597 Week 8
Michael Hallquist
12 Mar 2019
1 Global fit
Measures of global fit in SEM provide information about how well the model fits the data. Importantly, these statistics attempt to quantify the overall recovery of the observed data without typically considering specific components of fit or misfit in each element of the mean and covariance structure.
Global fit statistics can be divided into absolute and comparative fit indices. Absolute indices are often a function of the test statistic \(T\), which quantifies global fit to the population covariance structure (see model chi-square goodness of fit below). Absolute fit indices can also be a function of the model residuals.
Comparative fit indices compare a candidate model (specified by you) against a baseline model, which is a minimal model containing only variances for observed endogenous variables, but not covariances among them. Thus, the baseline model represents the view that there are no meaningful relationships among variables. Comparative fit indices describe how much better your model fits the data compared to this independence representation.
Global fit statistics can also be divided into tests of goodness versus badness of model fit, so it’s important to remember what are acceptable ranges for each.
1.1 Matrix expression of CFA
\[ \boldsymbol{\Sigma} = \boldsymbol{\Lambda}_y \boldsymbol{\Psi} \boldsymbol{\Lambda}_y' + \boldsymbol{\Theta}_\varepsilon \]
Reminders about SEM parameterization:
- The number of observed variables (sometimes called ‘observations’) is \(k\).
- The sample covariance matrix, estimated in a sample of \(N\) individuals, is \(\boldsymbol{S_{XX}}\).
- The population covariance matrix (unknown) is denoted \(\boldsymbol{\Sigma}\).
- The model-implied covariance matrix, given \(q\) free (estimated) parameters, is \(\boldsymbol{\Sigma}({\hat{\boldsymbol{\theta}}})\).
- The number of unique elements in a covariance matrix is: \(p = k (k+1) / 2\).
- SEM tries to minimize the discrepancy between \(\boldsymbol{S_{XX}}\) and \(\boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})\) according to an objective (fit) function: \(F = (\textbf{s} - \hat{\boldsymbol{\sigma}}(\theta))' \boldsymbol{W}^{-1} (\textbf{s} - \hat{\boldsymbol{\sigma}}(\theta))\).
- In standard ML estimation, the objective function is: \(\hat{F} = \textrm{log} | \boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}}) | + \textrm{tr} | \boldsymbol{S_{XX}} \boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})^{-1} | - \textrm{log} | \boldsymbol{S_{XX}} | - k\).
1.2 Model chi-square goodness of fit
The model \(\chi^2\) test is the most common global fit index in SEM and is a component of several other fit indices. It is a test of the null hypothesis that the model-implied covariance matrix \(\boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})\) equals the population covariance \(\boldsymbol{\Sigma}\). Therefore, we would prefer not to reject the null hypothesis since doing so would be indicative of misfit.
In standard ML, model \(\chi^2\) is computed as the product of \(N - 1\) with the minimum value of \(\hat{F}\) estimated in the optimization algorithm (often expectation-maximization [EM]). The minimum value, denoted \(f\), represents the fit of the model to the data given a vector of parameters with the greatest sample likelihood. The model \(\chi^2\) is then computed:
\[ T = (N-1) f \] which is (asymptotically) distributed as \(\chi^2\) with \(df = p - q\). Thus, we can test the null hypothesis of perfect model fit in the population given value \(T\) and \(df\). People often test this hypothesis at \(\alpha = .05\), but as we will discuss below, the \(\chi^2\) test can be ‘significant’ in cases where there is only trivial misfit, especially in large samples.
1.2.1 Assumptions of conventional model chi-square test
- Observed variables are distributed as multivariate normal.
- \(N\) is sufficiently large (approx. 150+ observations).
- None of the tested parameters is at a boundary or invalid (e.g., variance of zero).
1.2.2 Relationship to likelihood ratio test (LRT)
The model \(\chi^2\) is a specific case of a likelihood ratio test (LRT) in which two models are compared according to differences in their log likelihoods (LL) and degrees of freedom. In the case of the model \(\chi^2\) test, we are comparing the fit of our proposed SEM against the perfect fit of a saturated model, which estimates all variances and covariances individually (i.e., using \(p\) parameters).
More generally, LRTs can be used to test fit differences in nested models, where model A is considered nested in model B if the free parameters of A are a subset of the parameters in B.
In this situation, the difference in model fit can be computed as the difference in model \(\chi^2\) values, which is also distributed according to a \(\chi^2\) distribution. The degrees of freedom for a nested LRT is the difference in the individual degrees of freedom for each model:
\[ \begin{align*} \Delta \chi^2 &= \chi^2_\textrm{Fewer} - \chi^2_\textrm{More} \\ \Delta df &= df_\textrm{Fewer} - df_\textrm{More} \end{align*} \] where \(Fewer\) refers to the model with fewer parameters that is nested in the \(More\) model.
1.3 Comparing candidate models against a saturated model
Let’s walk through an example that highlights model nesting and chi-square model comparisons.
Say that we are interested in how well fear of global warming (FoGWar) and concern for animal welfare (CoAWe) predict donations to charities involved in environmental protection. We measure FoGWar and CoAWe using two self-report scales with four items each. We measure donations in terms of total dollars to five large charities most involved in environmental protection (i.e., observed variable). We predict that latent standing (factor scores) on both FoGWar and CoAWe will be associated with greater donations.
1.3.1 Observed correlations
#simulate data along the hypothesized lines above
demo.model <- '
l_FoGWar =~ .8*FoGWar1 + .8*FoGWar2 + .8*FoGWar3 + .8*FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ .8*CoAWe5 + .8*CoAWe6 + .8*CoAWe7 + .8*CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_FoGWar =~ .2*CoAWe5 #small cross-loading
donations ~ 0.3*l_FoGWar + 0.3*l_CoAWe
l_FoGWar ~~ 0.4*l_CoAWe
#equal (unit) residual variance
FoGWar1 ~~ 1*FoGWar1
FoGWar2 ~~ 1*FoGWar2
FoGWar3 ~~ 1*FoGWar3
FoGWar4 ~~ 1*FoGWar4
CoAWe5 ~~ 1*CoAWe5
CoAWe6 ~~ 1*CoAWe6
CoAWe7 ~~ 1*CoAWe7
CoAWe8 ~~ 1*CoAWe8
'
# generate data; note, standardized lv is default
simData <- lavaan::simulateData(demo.model, sample.nobs=500)
#round(cor(simData), 2)
ggcorrplot::ggcorrplot(cor(simData), type="lower")
1.3.2 Fit a saturated model
Here is the lavaan model syntax for the saturated model.
#pairwise combinations of observed variables
pairwise <- combn(x = names(simData), m = 2)
msat.syntax <- paste(apply(pairwise, 2, function(pair) { paste(pair, collapse="~~")}), collapse="\n")
msat <- sem(msat.syntax, simData)
cat(msat.syntax)## FoGWar1~~FoGWar2
## FoGWar1~~FoGWar3
## FoGWar1~~FoGWar4
## FoGWar1~~CoAWe5
## FoGWar1~~CoAWe6
## FoGWar1~~CoAWe7
## FoGWar1~~CoAWe8
## FoGWar1~~donations
## FoGWar2~~FoGWar3
## FoGWar2~~FoGWar4
## FoGWar2~~CoAWe5
## FoGWar2~~CoAWe6
## FoGWar2~~CoAWe7
## FoGWar2~~CoAWe8
## FoGWar2~~donations
## FoGWar3~~FoGWar4
## FoGWar3~~CoAWe5
## FoGWar3~~CoAWe6
## FoGWar3~~CoAWe7
## FoGWar3~~CoAWe8
## FoGWar3~~donations
## FoGWar4~~CoAWe5
## FoGWar4~~CoAWe6
## FoGWar4~~CoAWe7
## FoGWar4~~CoAWe8
## FoGWar4~~donations
## CoAWe5~~CoAWe6
## CoAWe5~~CoAWe7
## CoAWe5~~CoAWe8
## CoAWe5~~donations
## CoAWe6~~CoAWe7
## CoAWe6~~CoAWe8
## CoAWe6~~donations
## CoAWe7~~CoAWe8
## CoAWe7~~donations
## CoAWe8~~donationsFirst, let’s check the parameterization of the saturated model. Is it identified? How many degrees of freedom?
inspect(msat)## $lambda
##
## FoGWar1
## FoGWar2
## FoGWar3
## FoGWar4
## CoAWe5
## CoAWe6
## CoAWe7
## CoAWe8
## donations
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe5 CoAWe6 CoAWe7 CoAWe8 dontns
## FoGWar1 37
## FoGWar2 1 38
## FoGWar3 2 9 39
## FoGWar4 3 10 16 40
## CoAWe5 4 11 17 22 41
## CoAWe6 5 12 18 23 27 42
## CoAWe7 6 13 19 24 28 31 43
## CoAWe8 7 14 20 25 29 32 34 44
## donations 8 15 21 26 30 33 35 36 45
##
## $psi
## <0 x 0 matrix>Next, let’s examine fit and parameter values:
summary(msat, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 41 iterations
##
## Optimization method NLMINB
## Number of free parameters 45
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 0.000
## Degrees of freedom 0
## Minimum Function Value 0.0000000000000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 945.926
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7025.937
## Loglikelihood unrestricted model (H1) -7025.937
##
## Number of free parameters 45
## Akaike (AIC) 14141.875
## Bayesian (BIC) 14331.532
## Sample-size adjusted Bayesian (BIC) 14188.700
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.000
## P-value RMSEA <= 0.05 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## FoGWar1 ~~
## FoGWar2 0.581 0.080 7.291 0.000 0.581 0.345
## FoGWar3 0.559 0.079 7.079 0.000 0.559 0.334
## FoGWar4 0.588 0.079 7.419 0.000 0.588 0.352
## CoAWe5 0.299 0.079 3.765 0.000 0.299 0.171
## CoAWe6 0.263 0.073 3.580 0.000 0.263 0.162
## CoAWe7 0.212 0.080 2.646 0.008 0.212 0.119
## CoAWe8 0.240 0.075 3.195 0.001 0.240 0.144
## donations 0.231 0.064 3.613 0.000 0.231 0.164
## FoGWar2 ~~
## FoGWar3 0.727 0.082 8.835 0.000 0.727 0.430
## FoGWar4 0.668 0.081 8.233 0.000 0.668 0.396
## CoAWe5 0.468 0.082 5.732 0.000 0.468 0.265
## CoAWe6 0.380 0.075 5.065 0.000 0.380 0.233
## CoAWe7 0.319 0.082 3.905 0.000 0.319 0.177
## CoAWe8 0.416 0.077 5.393 0.000 0.416 0.249
## donations 0.368 0.066 5.596 0.000 0.368 0.258
## FoGWar3 ~~
## FoGWar4 0.669 0.081 8.283 0.000 0.669 0.399
## CoAWe5 0.386 0.080 4.816 0.000 0.386 0.221
## CoAWe6 0.118 0.073 1.613 0.107 0.118 0.072
## CoAWe7 0.269 0.081 3.323 0.001 0.269 0.150
## CoAWe8 0.289 0.075 3.825 0.000 0.289 0.174
## donations 0.224 0.064 3.499 0.000 0.224 0.158
## FoGWar4 ~~
## CoAWe5 0.308 0.079 3.871 0.000 0.308 0.176
## CoAWe6 0.236 0.073 3.212 0.001 0.236 0.145
## CoAWe7 0.129 0.080 1.614 0.106 0.129 0.072
## CoAWe8 0.274 0.075 3.634 0.000 0.274 0.165
## donations 0.312 0.065 4.822 0.000 0.312 0.221
## CoAWe5 ~~
## CoAWe6 0.748 0.083 9.025 0.000 0.748 0.441
## CoAWe7 0.873 0.092 9.475 0.000 0.873 0.468
## CoAWe8 0.795 0.085 9.308 0.000 0.795 0.458
## donations 0.349 0.068 5.144 0.000 0.349 0.236
## CoAWe6 ~~
## CoAWe7 0.749 0.084 8.885 0.000 0.749 0.433
## CoAWe8 0.616 0.077 7.987 0.000 0.616 0.382
## donations 0.405 0.064 6.339 0.000 0.405 0.296
## CoAWe7 ~~
## CoAWe8 0.760 0.086 8.826 0.000 0.760 0.430
## donations 0.389 0.070 5.583 0.000 0.389 0.258
## CoAWe8 ~~
## donations 0.365 0.065 5.632 0.000 0.365 0.260
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## FoGWar1 1.671 0.106 15.811 0.000 1.671 1.000
## FoGWar2 1.701 0.108 15.811 0.000 1.701 1.000
## FoGWar3 1.678 0.106 15.811 0.000 1.678 1.000
## FoGWar4 1.675 0.106 15.811 0.000 1.675 1.000
## CoAWe5 1.829 0.116 15.811 0.000 1.829 1.000
## CoAWe6 1.573 0.099 15.811 0.000 1.573 1.000
## CoAWe7 1.902 0.120 15.811 0.000 1.902 1.000
## CoAWe8 1.647 0.104 15.811 0.000 1.647 1.000
## donations 1.195 0.076 15.811 0.000 1.195 1.000semPaths_default(msat, layout="circle")
Note that the \(LL\) for our model (notated H0) is identical to the saturated model implicitly fit by lavaan (notated H1). This is because we explicitly fit a saturated model, though any just identified model would have equivalent fit. So, yes, the saturated model is just identified (\(df = 0\)), and it has 45 free parameters.
1.3.3 Fit the true SEM model that we used to generate the data
Note that in real life, we never know the ‘true’ model. To make this example more like real data, I’ve simulated the data with a small cross loading of CoAWe5 (the first item on the concern for animal welfare scale) onto the fear of global warming factor. Here, we include this in the SEM to verify that we have a good fit to the data when we use a good model (the groundtruth in this case).
msem.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_FoGWar =~ CoAWe5 #allow small cross-loading
donations ~ l_FoGWar + l_CoAWe #predict criterion based on SEM
l_FoGWar ~~ l_CoAWe
'
msem <- sem(msem.syntax, simData)
semPaths_default(msem)
What about identification and parameterization?
How many free parameters do we have?
inspect(msem)## $lambda
## l_FGWr l_CoAW dontns
## FoGWar1 0 0 0
## FoGWar2 1 0 0
## FoGWar3 2 0 0
## FoGWar4 3 0 0
## CoAWe5 7 0 0
## CoAWe6 0 4 0
## CoAWe7 0 5 0
## CoAWe8 0 6 0
## donations 0 0 0
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe5 CoAWe6 CoAWe7 CoAWe8 dontns
## FoGWar1 11
## FoGWar2 0 12
## FoGWar3 0 0 13
## FoGWar4 0 0 0 14
## CoAWe5 0 0 0 0 15
## CoAWe6 0 0 0 0 0 16
## CoAWe7 0 0 0 0 0 0 17
## CoAWe8 0 0 0 0 0 0 0 18
## donations 0 0 0 0 0 0 0 0 0
##
## $psi
## l_FGWr l_CoAW dontns
## l_FoGWar 20
## l_CoAWe 10 21
## donations 0 0 19
##
## $beta
## l_FGWr l_CoAW dontns
## l_FoGWar 0 0 0
## l_CoAWe 0 0 0
## donations 8 9 0summary(msem, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 35 iterations
##
## Optimization method NLMINB
## Number of free parameters 21
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 33.905
## Degrees of freedom 24
## P-value (Chi-square) 0.086
##
## Model test baseline model:
##
## Minimum Function Test Statistic 945.926
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.989
## Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7042.890
## Loglikelihood unrestricted model (H1) -7025.937
##
## Number of free parameters 21
## Akaike (AIC) 14127.780
## Bayesian (BIC) 14216.287
## Sample-size adjusted Bayesian (BIC) 14149.632
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.029
## 90 Percent Confidence Interval 0.000 0.050
## P-value RMSEA <= 0.05 0.954
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.030
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.688 0.532
## FoGWar2 1.288 0.142 9.069 0.000 0.886 0.679
## FoGWar3 1.183 0.134 8.814 0.000 0.814 0.628
## FoGWar4 1.152 0.132 8.711 0.000 0.792 0.612
## l_CoAWe =~
## CoAWe5 1.000 0.897 0.663
## CoAWe6 0.889 0.090 9.891 0.000 0.797 0.636
## CoAWe7 1.030 0.102 10.134 0.000 0.924 0.670
## CoAWe8 0.920 0.092 9.950 0.000 0.825 0.643
## l_FoGWar =~
## CoAWe5 0.151 0.103 1.463 0.144 0.104 0.077
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.328 0.094 3.491 0.000 0.226 0.206
## l_CoAWe 0.374 0.071 5.273 0.000 0.335 0.307
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.245 0.050 4.884 0.000 0.397 0.397
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.198 0.089 13.466 0.000 1.198 0.717
## .FoGWar2 0.916 0.086 10.641 0.000 0.916 0.539
## .FoGWar3 1.016 0.086 11.876 0.000 1.016 0.606
## .FoGWar4 1.047 0.086 12.203 0.000 1.047 0.625
## .CoAWe5 0.940 0.083 11.340 0.000 0.940 0.514
## .CoAWe6 0.938 0.075 12.466 0.000 0.938 0.596
## .CoAWe7 1.049 0.089 11.767 0.000 1.049 0.551
## .CoAWe8 0.966 0.078 12.327 0.000 0.966 0.587
## .donations 0.971 0.065 14.956 0.000 0.971 0.813
## l_FoGWar 0.473 0.087 5.423 0.000 1.000 1.000
## l_CoAWe 0.804 0.127 6.348 0.000 1.000 1.000We see strong evidence of global fit: we do not reject the model \(\chi^2\), CFI is > .95, RMSEA is < .06 (and its confidence interval does not include .06), and SRMR is < .08. Thus, at a global level, this model characterizes the data well.
1.3.4 Compare the saturated model with the candidate model
Let’s compare the saturated model against the specified two-factor SEM. Technically, all SEM models are nested in the saturated model.
anova(msat, msem)| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| msat | 0 | 14142 | 14332 | 0.0 | NA | NA | NA |
| msem | 24 | 14128 | 14216 | 33.9 | 33.9 | 24 | 0.086 |
Notice how the LRT exactly equals the model chi-square test of the two-factor SEM specified above. This is because the model chi-square test is a nested model comparison against the saturated model. An important point is that our candidate model has 24 fewer parameters than the saturated model but does not fit significantly worse. That’s evidence of a good fit:parsimony tradeoff.
1.3.5 Compare a slightly simpler nested model
What about a nested model in which we remove one of the measurement parameters? In particular, the cross-loading is an interpretive nuisance (not simple structure). Let’s remove the cross-loading. We are fixing this path to be exactly zero, which means that we are estimating one less parameter, but the rest of the model structure is the same. Thus, the simpler (reduced) model is nested in the more complex model.
mnested.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
donations ~ l_FoGWar + l_CoAWe #predict criterion based on SEM
l_FoGWar ~~ l_CoAWe
'
mnested <- sem(mnested.syntax, simData)
summary(mnested, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 34 iterations
##
## Optimization method NLMINB
## Number of free parameters 20
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 36.005
## Degrees of freedom 25
## P-value (Chi-square) 0.072
##
## Model test baseline model:
##
## Minimum Function Test Statistic 945.926
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.988
## Tucker-Lewis Index (TLI) 0.983
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7043.940
## Loglikelihood unrestricted model (H1) -7025.937
##
## Number of free parameters 20
## Akaike (AIC) 14127.880
## Bayesian (BIC) 14212.172
## Sample-size adjusted Bayesian (BIC) 14148.690
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.030
## 90 Percent Confidence Interval 0.000 0.050
## P-value RMSEA <= 0.05 0.951
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.031
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.688 0.533
## FoGWar2 1.289 0.142 9.073 0.000 0.887 0.680
## FoGWar3 1.177 0.134 8.798 0.000 0.810 0.625
## FoGWar4 1.150 0.132 8.710 0.000 0.792 0.612
## l_CoAWe =~
## CoAWe5 1.000 0.952 0.704
## CoAWe6 0.832 0.074 11.288 0.000 0.793 0.632
## CoAWe7 0.960 0.082 11.666 0.000 0.914 0.663
## CoAWe8 0.866 0.076 11.420 0.000 0.824 0.642
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.318 0.095 3.365 0.001 0.219 0.200
## l_CoAWe 0.353 0.067 5.294 0.000 0.336 0.307
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.279 0.049 5.660 0.000 0.426 0.426
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.197 0.089 13.452 0.000 1.197 0.716
## .FoGWar2 0.913 0.086 10.594 0.000 0.913 0.537
## .FoGWar3 1.022 0.086 11.921 0.000 1.022 0.609
## .FoGWar4 1.048 0.086 12.195 0.000 1.048 0.626
## .CoAWe5 0.922 0.083 11.084 0.000 0.922 0.504
## .CoAWe6 0.945 0.075 12.636 0.000 0.945 0.601
## .CoAWe7 1.067 0.088 12.054 0.000 1.067 0.561
## .CoAWe8 0.968 0.078 12.452 0.000 0.968 0.587
## .donations 0.971 0.065 14.957 0.000 0.971 0.813
## l_FoGWar 0.474 0.087 5.428 0.000 1.000 1.000
## l_CoAWe 0.907 0.116 7.807 0.000 1.000 1.000Note that the model chi-square is still non-significant, suggesting good fit relative to the population covariance matrix. But is the fit worse than the model that includes the cross-loading?
anova(msem, mnested)| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| msem | 24 | 14128 | 14216 | 33.9 | NA | NA | NA |
| mnested | 25 | 14128 | 14212 | 36.0 | 2.1 | 1 | 0.147 |
No, we do not see a significant degradation in fit by excluding the cross-loading parameter, suggesting we should retain the model with simple structure. That is, the LRT \(p\) value of .15 suggests that the model with a freed cross-loading parameter is not significantly better than the model that the model with simple structure.
1.4 Challenges of using the model chi-square test
Although the model \(\chi^2\) is a common and useful test of global fit, it is not without flaws. Foremost among these is that the power to detect trivial misfit increases toward 1.0 as the sample size \(N\) increases. That is, in large samples, we will reject the null hypothesis of good fit even if the model is excellent. This problem applies to comparison among nested SEMs using the LRT (difference in \(\chi^2\)) approach, too.
A related problem is that overemphasis on model \(\chi^2\) may lead some researchers to prefer small samples in which the null hypothesis (of good fit) is not rejected. However, in small samples, parameters will be poorly estimated, and it is the parameters, not fit tests, that are usually of substantive interest. Thus, there is a risk of encouraging bad scientific practices (collecting small samples to avoid rejection the model \(\chi^2\) test).
A second problem is that the model statistic \(T\) does not follow a \(\chi^2\) distribution if the data are not multivariate normal or when \(N\) is small (see assumptions above). There have been additional developments in the SEM literature for estimators that are robust to non-normality, especially the Satorra-Bentler correction to the \(\chi^2\)
A third problem is that in small samples, even when other assumptions are met, the model \(\chi^2\) test can reject the null hypothesis even when the model fits well… so yeah, you can’t win (reject with small or large \(N\))! On top of that, in small samples, the test can fail to detect serious misfit.
1.5 Practical fit indices
As Kline notes, the model \(\chi^2\) tests a hypothesis of exact fit and is typically framed in the binary terms of accept or reject. Some of the problems with this the model \(\chi^2\) test noted above, along with its formulation as a null hypothesis significance test, have promoted the development of alternative measures of fit. These measures, sometimes called practical or approximate fit indices, seek to provide guidance on the extent to which the model fits the data reasonably well.
1.5.1 \(\chi^2/df\) ratio
The ratio of misfit, quantified by \(\chi^2\), to the model degrees of freedom, \(df\), is an absolute fit index. The intuition is that for every trivial free parameter removed from the model, the \(\chi^2\) should increase by about 1.0. That is, the expected value of \(\chi^2\) for a correct model is \(df\). Thus, in larger, more complex models, this ratio may be more indicative of fit than the \(\chi^2\). The \(\chi^2\), which is proportionate to the sample likelihood, is sensitive the case where there are many trivial deviations between the sample covariance matrix and the model-implied covariance matrix. Like \(R^2\) in regression, penalties for these deviations can never be negative, so the model \(\chi^2\) can become quite large in models with many free parameters.
The \(\chi^2/df\) ratio is intended to help correct for such problems. Rules of thumb for this ratio in well-fitting models are < 5, < 4, or sometimes even < 2. As the ratio becomes higher (6, 10, 50, etc.), we increasingly believe the model is a poor description of the data. Importantly, \(\chi^2/df\) is a badness of fit index such that we become concerned at higher values, but there is no lower bound that is ideal. Also, because of the sample size dependence of \(\chi^2\), the \(\chi^2/df\) ratio is also dependent on \(N\), leading it to be more of an ad hoc heuristic than a go-to fit index.
1.5.2 GFI and AGFI
The Goodness-of-Fit Index (GFI) and Adjusted Goodness-of-Fit Index (AGFI) provide an absolute metric for how well the model describes the data. The GFI is conceptually related to the proportion of variance explained by the model, as in regression. Thus, it ranges between 0 (bad) and 1 (perfect).
The AGFI extends this metric by adjusting for the \(df\) of the model, downweighting goodness in proportion to the number of free parameters. This is analogous to adjusted \(R^2\) in regression, where a model with many predictors is penalized in its overall explanatory power because of the positive variance explained by each.
Although these are reasonable metrics, they are not often used in practice.
1.5.3 RMR and SRMR
The root mean residual (RMR) and standardized root mean residual (SRMR) provide an estimate of the average misfit for each estimated versus observed variance/covariance parameter. That is, the RMR quantifies the square root of the average of the squared residuals comparing the model-implied covariance matrix and observed covariance matrix:
\[\sqrt{\sum_p (\boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})_p - \boldsymbol{S}_p)^2 } \]
where \(p\) is my sloppy notation for the pth unique element of the covariance matrix. Note that this gives an intuitive estimate of the level of average misfit in the model-implied covariance matrix in the units of the original data. The RMR is most interpretable when all variables are on a similar scale. RMR is a badness of fit index in that large values are indicative of misfit.
One virtue of RMR is that it is not sensitive to variables with small unique (residual) variances. The weighting procedure of many other fit measures (e.g., CFI or GFI) can misstate the badness of fit in this case.
The SRMR solves the problem of scaling the metric of RMR. Because RMR is in the original units, it is not comparable across studies and is hard to interpret when the variances of variables are disparate. The SRMR standardizes the average residuals in terms of the magnitude of elements in \(\boldsymbol{S_{XX}}\). Thus, it ranges between 0 and 1, where 1 represents residuals that are as large, on average, as a given element of the sample covariance matrix.
The rule of thumb for SRMR is < .08.
1.5.4 RMSEA
The root mean squared error of approximation (RMSEA) is a badness of fit index that quantifies level of approximate fit (as opposed to perfect fit in model \(\chi^2\)). The RMSEA seeks to examine whether a model exceeds a reasonable level of close fit to the data, which would indicate misfit. It is based on the fact that the null hypothesis test (model \(T\)) follows a \(\chi^2\) distribution, whereas evidence for the alternative follows a non-central \(\chi^2\) distribution, whose non-centrality is governed by a parameter \(\lambda\).
Conceptually, RMSEA penalizes more complex models by dividing the non-centrality term by the degrees of freedom. This is a form of a parsimony-adjusted index such that simpler models are preferred. Furthermore, RMSEA is in the same units as SRMR, namely the size of the average element in \(\boldsymbol{S_{XX}}\). Thus, it ranges from 0-1. Values above .10 indicate serious misfit. Values below .06 are considered good.
Many people encourage the best practice of reporting a 90% CI on RMSEA to see how close the model is to the boundary of misfit. Some programs also report a test on the null hypothesis that RMSEA is less than or equal to .05. One hopes not to reject this hypothesis since low RMSEA is a good thing.
Mechanics of RMSEA
Close fit is defined as:
\[ \hat{\Delta}_M = \textrm{max}(0, \chi^2_M - df_M) \]
When the \(\chi^2_M\) is less than \(df_M\), then this metric is zero, indicating no departure from close fit. The term \(\hat{\Delta}_M\) is the noncentrality parameter in a \(\chi^2\) distribution that examines how well the model fits the population covariance matrix. When \(\hat{\Delta}_M = 0\), this reduces to a standard \(\chi^2\) testing close fit. But in general, the model will not perfectly fit the population, and its deviation according to RMSEA is:
\[ \hat{\varepsilon} = \sqrt{\frac{\hat{\Delta}_M}{df_M (N-1) }} \] Thus, the intuition is that the deviation becomes smaller when a) \(df_M\) is larger (i.e., more parsimonious model) and b) the sample size increases. In larger samples, the correction for \(df_M\) will be much smaller than the correction for \(N\), indicating a weaker penalty for complex models in large samples
1.5.5 TLI/NNFI
The Tucker-Lewis index (TLI), which is also called the non-normed fit index (NNFI), ranges between 0 and 1, with larger values indicating better relative fit. The typical rule of thumb is to judge a model as good if TLI > .95. It is related to the \(\chi^2/df\) ratio described above, but instead of relying on rules of thumb, it compares the ratio of your model relative to a ‘true model’ (where the expected \(\chi^2/df_M\) ratio of such a model is 1) and the baseline (independence) model.
That is, the TLI represents the distance (in terms of fit) between your model and the independence model as a proportion of the distance between the independence model and the saturated model. Because the independence model is lousy by definition and the saturated model is perfect (fit-wise) by definition, the distance between these models (e.g., in terms of an LRT) represents a benchmark for how good (or bad) fit can be. Hence, if we can achieve 90 or 95% of this distance using a candidate model, but with many fewer parameters, we would believe our model fits well.
More formally, the TLI is defined as:
\[ TLI = 1 - \frac{(\chi^2_k - df_k) / (\chi^2_0 - df_0)}{df_k/df_0} \] where \(k\) is the candidate model and \(0\) is the independence model. Note that the TLI can technically go below 0 or above 1, which is confusing and unlikely, but nevertheless, know that this metrics is not formally bounded. If you see a value below 0, your model fits very badly (worse than independence), or paradoxically, the independence model fits well. If the value is above 1, your model fits very well, or the independence model is especially bad.
1.5.6 NFI
The normed fit index (NFI) is related to TLI in that it quantifies the distance between a candidate model and the baseline. It does not use the \(\chi^2/df\) ratios, but instead is based on the log-likelihood of the models or model \(\chi^2\). Similar to TLI, it represents the improvement in fit of a candidate model against the baseline (independence) as a proportion of the baseline fit.
NFI ranges between 0 and 1, and values > .95 are considered ‘good.’
\[ \textrm{NFI} = \frac{\chi^2_0 - \chi^2_k}{\chi^2_0} \]
Note that NFI does not penalize for complexity, which may encourage a researcher to overfit the data. Likewise, it is affected by sample size, with small \(N\) underestimating NFI.
In an attempt to handle this underestimation problem, Bollen introduced the incremental fit index (IFI), but this metric introduced several new problems and is not recommended.
1.5.7 CFI
The comparative fit index (CFI) is goodness of fit index that ranges between 0 and 1, with values above .95 indicating ‘good’ fit. It is conceptually related to the TLI in that it expresses the fit of a candidate model against an independence model. It is related to RMSEA in that both rely on the idea that the evidence for a candidate model follows a noncentral \(\chi^2\) distribution. More specifically, CFI considers how much the noncentrality parameter is reduced as one moves from a baseline model to a candidate model (in terms of fit).
The CFI is quantified in terms of the relative difference between fit (model \(\chi^2\)) and degrees of freedom for both candidate and baseline models:
\[ \textrm{CFI} = \frac{\textrm{max}(\chi^2_0 - df_0, 0) - \textrm{max}(\chi^2_k - df_k, 0)}{\textrm{max}(\chi^2_0 - df_0, 0)} \] By taking the max of each operand, CFI is formally bounded between 0 and 1. Importantly, CFI is not affected by \(N\) and has performed well in simulation studies relative to similar indices such as TLI. Finally, note that there is a very similar measure called relative noncentrality index (RNI) that differs from CFI primarily in that it can go above 1.0 in the same way as TLI.
1.6 Challenges of practical fit indices
Many, if not most, of the practical indices described above cannot be interpreted blindly according to the rules of thumb above. Indeed, these rules of thumb were primarily developed for models containing only covariance structure, but not means. We’ll get to mean structure later in the class, but suffice it to say that most fit indices struggle to balance misfit in the mean structure and covariance structure. Is a model that fits one well, but not the other, considered ‘good?’ For example, SRMR typically ignores mean structure altogether (since it is based on comparisons of covariance elements), so a model with poor representation of means could have a low SRMR.
There is also a challenge of specifying an appropriate baseline model in more complex SEMs, especially ones that include mean structure. Likewise, some models do not have a proper saturated model, such as models with specific patterns of missing data. The short version is to interpret global fit indices with caution and always to examine the fit of each part of the model (residuals are useful).
One may be tempted to modify a model according to modification indices in order to improve global fit, but atheoretically (fit) motivated modifications are unlikely to replicate in an independent sample. As West notes, model modification can often shift the epistemological perspective from confirmatory (a priori model) to (partially) exploratory. Exercise caution!
1.7 tl;dr summary of global fit statistics
Kline suggests reporting the following minimal set of fit indices (and I agree):
- Model \(\chi^2\) along with its \(df\) and \(p\).
- standardized root mean residual (SRMR): .08 or less.
- root mean squared error of approximation (RMSEA), along with 90% CI: .06 or less
- comparative fit index (CFI): .95 or greater.
Altogether, many of the other fit indices that have been developed for SEM try to overcome dependence on \(N\). Here’s the tl;dr synopsis from the West et al. chapter.
2 Local fit
This brings us to the topic of ‘local fit’, by which I mean how well each part of the model captures the data. For example, in a two-factor model, does each factor fit the data well in isolation? Things typically only get worse when we consider a model built from poorly fitting elements. That is, global fit statistics can be:
- sensitive to misfit in a single bad component amidst many good components
- insensitive to a bad component when many are good, and/or
- sensitive to mild misfit in many reasonable components when these are considered together.
Tomarken and Waller (2003, J Abnormal Psych) argue for the importance of testing meaningful partitions of the model in order to verify good fit in each. Most commonly, one can separately test the measurement and structural models. But in larger models (esp. multi-factor or hierarchical factor), testing individual components of measurement is often important. Likewise, testing a sub-part of the SEM (e.g., both measurement and structural components for self-report data, but not the informant data) is useful.
They note the following problems with models that appear to fit well according to global statistics:
- There are many equivalent models that fit equally well
- There may be several alternative models that could fit better
- Important variables may be omitted (that would affect fit)
- Fit of lower-order components may be poor
- An ill-fitting partition (component) may be masked within a well-fitting overall model
- Global stats may be insensitive to particular kinds of misspecification (esp. in mean structure models)
- One may make post hoc modifications to the model that reduce the validity and reproducibility of the results.
Tomarken and Waller also offer the following advice for avoiding overreliance on global fit statistics:
- Acknowledge plausible equivalent models and design studies to rule out these alternatives
- Compare the fit of a candidate model against nonequivalent alternatives (e.g., one-factor versus three-factor)
- Acknowledge potentially important variables that were omitted from the model or not measured. How might these alter the conclusions?
- Report and evaluate lower-order model components
- Test components of the model to verify fit in each (this can proceed hierarchically)
- Design studies with sensitivity to detect non-trivial misspecifications and include measures to flag trivial misspecifications (I think this is quite hard in practice)
- Clearly distinguish between a priori models and post hoc modifications, or models generated from a search process.
2.1 Testing fit of components
Let’s return to the example above about fear of global warming (FoGWar) and concern for animal welfare (CoAWe). How might we test local fit? First, let me mess up the model a bit to make it complex so we can encounter fit challenges. :-)
In addition, let’s make the structural model a bit more complex by testing the hypothesis that the effects of FoGWar and CoAWe on worldwide governmental spending on the environment (observed, in USD) are mediated by charitable donations. That is, we believe the environmental NGOs (charities) influence governments through lobbying and the development of legislation. Thus, if these NGOs have a larger budget, they can effect greater change in governmental spending.
demo.model <- '
l_FoGWar =~ .8*FoGWar1 + .8*FoGWar2 + .8*FoGWar3 + .6*FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ .6*CoAWe5 + .8*CoAWe6 + .8*CoAWe7 + .8*CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_FoGWar =~ .6*CoAWe5 #cross-loading
l_CoAWe =~ .6*FoGWar4 #cross-loading
donations ~ 0.3*l_FoGWar + 0.3*l_CoAWe
gspending ~ 0.5*donations + 0.3*l_FoGWar + 0.3*l_CoAWe
l_FoGWar ~~ 0.4*l_CoAWe
#equal (unit) residual variance
FoGWar1 ~~ 1*FoGWar1
FoGWar2 ~~ 1*FoGWar2
FoGWar3 ~~ 1*FoGWar3
FoGWar4 ~~ 1*FoGWar4
CoAWe5 ~~ 1*CoAWe5
CoAWe6 ~~ 1*CoAWe6
CoAWe7 ~~ 1*CoAWe7
CoAWe8 ~~ 1*CoAWe8
'
# generate data; note, standardized lv is default
simData <- lavaan::simulateData(demo.model, sample.nobs=500)
ggcorrplot::ggcorrplot(cor(simData), type = "lower")
#round(cor(simData), 2)Let’s test the same model as above, but in the new dataset, and with the addition of gspending. We will start with the model having simple structure (no cross-loadings for FoGWar or CoAWe).
msem.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
donations ~ l_FoGWar + l_CoAWe #predict criterion based on SEM
gspending ~ donations
l_FoGWar ~~ l_CoAWe
'
msem <- sem(msem.syntax, simData)
semPaths_default(msem)
inspect(msem)## $lambda
## l_FGWr l_CoAW dontns gspndn
## FoGWar1 0 0 0 0
## FoGWar2 1 0 0 0
## FoGWar3 2 0 0 0
## FoGWar4 3 0 0 0
## CoAWe5 0 0 0 0
## CoAWe6 0 4 0 0
## CoAWe7 0 5 0 0
## CoAWe8 0 6 0 0
## donations 0 0 0 0
## gspending 0 0 0 0
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe5 CoAWe6 CoAWe7 CoAWe8 dontns
## FoGWar1 11
## FoGWar2 0 12
## FoGWar3 0 0 13
## FoGWar4 0 0 0 14
## CoAWe5 0 0 0 0 15
## CoAWe6 0 0 0 0 0 16
## CoAWe7 0 0 0 0 0 0 17
## CoAWe8 0 0 0 0 0 0 0 18
## donations 0 0 0 0 0 0 0 0 0
## gspending 0 0 0 0 0 0 0 0 0
## gspndn
## FoGWar1
## FoGWar2
## FoGWar3
## FoGWar4
## CoAWe5
## CoAWe6
## CoAWe7
## CoAWe8
## donations
## gspending 0
##
## $psi
## l_FGWr l_CoAW dontns gspndn
## l_FoGWar 21
## l_CoAWe 10 22
## donations 0 0 19
## gspending 0 0 0 20
##
## $beta
## l_FGWr l_CoAW dontns gspndn
## l_FoGWar 0 0 0 0
## l_CoAWe 0 0 0 0
## donations 7 8 0 0
## gspending 0 0 9 0Any comments or concerns so far? We see that the latent factors for fear of global warming and concern for animal welfare predict individual donations, which in turn predict governmental spending. This fit with our hypotheses.
summary(msem, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 37 iterations
##
## Optimization method NLMINB
## Number of free parameters 22
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 153.328
## Degrees of freedom 33
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 1405.036
## Degrees of freedom 45
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.912
## Tucker-Lewis Index (TLI) 0.879
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7899.614
## Loglikelihood unrestricted model (H1) -7822.950
##
## Number of free parameters 22
## Akaike (AIC) 15843.227
## Bayesian (BIC) 15935.949
## Sample-size adjusted Bayesian (BIC) 15866.119
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.085
## 90 Percent Confidence Interval 0.072 0.099
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.077
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.774 0.584
## FoGWar2 0.948 0.103 9.252 0.000 0.734 0.542
## FoGWar3 0.903 0.099 9.141 0.000 0.699 0.533
## FoGWar4 1.451 0.130 11.199 0.000 1.123 0.745
## l_CoAWe =~
## CoAWe5 1.000 1.099 0.761
## CoAWe6 0.682 0.060 11.394 0.000 0.750 0.578
## CoAWe7 0.715 0.058 12.267 0.000 0.786 0.626
## CoAWe8 0.623 0.059 10.539 0.000 0.685 0.532
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.326 0.194 1.683 0.092 0.253 0.222
## l_CoAWe 0.272 0.136 2.005 0.045 0.299 0.263
## gspending ~
## donations 0.766 0.045 17.159 0.000 0.766 0.609
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.712 0.081 8.837 0.000 0.836 0.836
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.156 0.085 13.594 0.000 1.156 0.659
## .FoGWar2 1.299 0.092 14.047 0.000 1.299 0.707
## .FoGWar3 1.233 0.087 14.130 0.000 1.233 0.716
## .FoGWar4 1.010 0.098 10.296 0.000 1.010 0.445
## .CoAWe5 0.878 0.086 10.221 0.000 0.878 0.421
## .CoAWe6 1.124 0.081 13.832 0.000 1.124 0.666
## .CoAWe7 0.960 0.073 13.245 0.000 0.960 0.608
## .CoAWe8 1.187 0.083 14.259 0.000 1.187 0.717
## .donations 1.010 0.067 14.984 0.000 1.010 0.783
## .gspending 1.285 0.081 15.811 0.000 1.285 0.629
## l_FoGWar 0.599 0.094 6.365 0.000 1.000 1.000
## l_CoAWe 1.208 0.136 8.856 0.000 1.000 1.000In terms of global fit, we see that only SRMR (barely) suggests good fit. RMSEA is in the meh range (~ .08), but not terrible (> .10). CFI is not acceptable at .91. Even if these statistics were acceptable, it is usually worthwhile to examine the fit of individual components, such as the measurement model or specific factors.
2.2 Test FoGWar alone
What could be wrong? Let’s examine the FoGWar factor in isolation.
m1.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
'
m1 <- sem(m1.syntax, simData)
inspect(m1)## $lambda
## l_FGWr
## FoGWar1 0
## FoGWar2 1
## FoGWar3 2
## FoGWar4 3
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4
## FoGWar1 4
## FoGWar2 0 5
## FoGWar3 0 0 6
## FoGWar4 0 0 0 7
##
## $psi
## l_FGWr
## l_FoGWar 8semPaths_default(m1)
summary(m1, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 27 iterations
##
## Optimization method NLMINB
## Number of free parameters 8
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 1.751
## Degrees of freedom 2
## P-value (Chi-square) 0.417
##
## Model test baseline model:
##
## Minimum Function Test Statistic 330.007
## Degrees of freedom 6
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.002
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3307.663
## Loglikelihood unrestricted model (H1) -3306.787
##
## Number of free parameters 8
## Akaike (AIC) 6631.325
## Bayesian (BIC) 6665.042
## Sample-size adjusted Bayesian (BIC) 6639.650
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.085
## P-value RMSEA <= 0.05 0.740
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.012
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.863 0.651
## FoGWar2 0.906 0.100 9.013 0.000 0.781 0.576
## FoGWar3 0.881 0.097 9.038 0.000 0.760 0.579
## FoGWar4 1.116 0.118 9.440 0.000 0.963 0.639
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.011 0.094 10.753 0.000 1.011 0.576
## .FoGWar2 1.228 0.099 12.434 0.000 1.228 0.668
## .FoGWar3 1.144 0.092 12.382 0.000 1.144 0.665
## .FoGWar4 1.345 0.121 11.073 0.000 1.345 0.592
## l_FoGWar 0.745 0.114 6.532 0.000 1.000 1.000No evidence of misfit. All looks good.
2.3 Test CoAWe alone
m2.syntax <- '
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
'
m2 <- sem(m2.syntax, simData)
inspect(m2)## $lambda
## l_CoAW
## CoAWe5 0
## CoAWe6 1
## CoAWe7 2
## CoAWe8 3
##
## $theta
## CoAWe5 CoAWe6 CoAWe7 CoAWe8
## CoAWe5 4
## CoAWe6 0 5
## CoAWe7 0 0 6
## CoAWe8 0 0 0 7
##
## $psi
## l_CoAW
## l_CoAWe 8semPaths_default(m2)
summary(m2, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 25 iterations
##
## Optimization method NLMINB
## Number of free parameters 8
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 0.987
## Degrees of freedom 2
## P-value (Chi-square) 0.611
##
## Model test baseline model:
##
## Minimum Function Test Statistic 380.404
## Degrees of freedom 6
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.008
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3202.985
## Loglikelihood unrestricted model (H1) -3202.491
##
## Number of free parameters 8
## Akaike (AIC) 6421.969
## Bayesian (BIC) 6455.686
## Sample-size adjusted Bayesian (BIC) 6430.294
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.072
## P-value RMSEA <= 0.05 0.854
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.008
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_CoAWe =~
## CoAWe5 1.000 0.957 0.663
## CoAWe6 0.833 0.084 9.906 0.000 0.797 0.614
## CoAWe7 0.874 0.085 10.286 0.000 0.836 0.665
## CoAWe8 0.804 0.082 9.752 0.000 0.769 0.598
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .CoAWe5 1.170 0.106 11.050 0.000 1.170 0.561
## .CoAWe6 1.052 0.086 12.163 0.000 1.052 0.623
## .CoAWe7 0.880 0.080 10.979 0.000 0.880 0.557
## .CoAWe8 1.065 0.085 12.467 0.000 1.065 0.643
## l_CoAWe 0.916 0.133 6.889 0.000 1.000 1.000Everything is excellent here. No cause for alarm.
2.4 Measurement model altogether
Let’s see if part of the misfit reflects the difficulty of modeling FoGWar and CoAWe together. Note that the fits of each factor alone don’t consider the possibility that there is important variation among indicators across factors, such as cross-loading or important residual covariance terms. By zooming in on a single factor, we can’t know anything about covariances across factors, either at the indicator level (residual covariance) or factor level (factor correlation).
m3.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_FoGWar ~~ l_CoAWe
'
m3 <- sem(m3.syntax, simData)
inspect(m3)## $lambda
## l_FGWr l_CoAW
## FoGWar1 0 0
## FoGWar2 1 0
## FoGWar3 2 0
## FoGWar4 3 0
## CoAWe5 0 0
## CoAWe6 0 4
## CoAWe7 0 5
## CoAWe8 0 6
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe5 CoAWe6 CoAWe7 CoAWe8
## FoGWar1 8
## FoGWar2 0 9
## FoGWar3 0 0 10
## FoGWar4 0 0 0 11
## CoAWe5 0 0 0 0 12
## CoAWe6 0 0 0 0 0 13
## CoAWe7 0 0 0 0 0 0 14
## CoAWe8 0 0 0 0 0 0 0 15
##
## $psi
## l_FGWr l_CoAW
## l_FoGWar 16
## l_CoAWe 7 17semPaths_default(m3)
summary(m3, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 32 iterations
##
## Optimization method NLMINB
## Number of free parameters 17
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 73.440
## Degrees of freedom 19
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 994.023
## Degrees of freedom 28
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.944
## Tucker-Lewis Index (TLI) 0.917
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -6404.193
## Loglikelihood unrestricted model (H1) -6367.472
##
## Number of free parameters 17
## Akaike (AIC) 12842.385
## Bayesian (BIC) 12914.033
## Sample-size adjusted Bayesian (BIC) 12860.074
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.076
## 90 Percent Confidence Interval 0.058 0.094
## P-value RMSEA <= 0.05 0.010
##
## 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|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.781 0.589
## FoGWar2 0.944 0.102 9.251 0.000 0.737 0.544
## FoGWar3 0.904 0.098 9.182 0.000 0.706 0.538
## FoGWar4 1.425 0.128 11.101 0.000 1.113 0.738
## l_CoAWe =~
## CoAWe5 1.000 1.105 0.765
## CoAWe6 0.677 0.060 11.252 0.000 0.748 0.576
## CoAWe7 0.710 0.059 12.109 0.000 0.785 0.625
## CoAWe8 0.615 0.059 10.373 0.000 0.679 0.528
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.721 0.081 8.860 0.000 0.836 0.836
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.146 0.085 13.432 0.000 1.146 0.653
## .FoGWar2 1.295 0.093 13.949 0.000 1.295 0.705
## .FoGWar3 1.223 0.087 14.005 0.000 1.223 0.711
## .FoGWar4 1.034 0.100 10.330 0.000 1.034 0.455
## .CoAWe5 0.865 0.088 9.877 0.000 0.865 0.415
## .CoAWe6 1.128 0.082 13.768 0.000 1.128 0.669
## .CoAWe7 0.963 0.073 13.151 0.000 0.963 0.610
## .CoAWe8 1.195 0.084 14.230 0.000 1.195 0.721
## l_FoGWar 0.610 0.095 6.391 0.000 1.000 1.000
## l_CoAWe 1.221 0.138 8.838 0.000 1.000 1.000Yep, something is not quite right in the measurement model. Note the CFI is slightly below .95 and RMSEA is greater than .06. The model \(\chi^2\) is also quite high relative to the \(df\).
Let’s examine the modification indices. Here, I will only consider \(\chi^2\) improvements of > 2 points, and I arrange the MIs in descending order (print first 8 only):
modificationindices(m3) %>% filter(mi > 2) %>% arrange(desc(mi)) %>% head(n=8)| lhs | op | rhs | mi | epc | sepc.lv | sepc.all | sepc.nox |
|---|---|---|---|---|---|---|---|
| l_FoGWar | =~ | CoAWe5 | 27.59 | 1.640 | 1.280 | 0.886 | 0.886 |
| l_CoAWe | =~ | FoGWar4 | 20.29 | 1.084 | 1.198 | 0.795 | 0.795 |
| FoGWar2 | ~~ | CoAWe5 | 16.26 | 0.251 | 0.251 | 0.237 | 0.237 |
| l_FoGWar | =~ | CoAWe8 | 11.38 | -0.778 | -0.608 | -0.472 | -0.472 |
| FoGWar2 | ~~ | FoGWar4 | 10.34 | -0.249 | -0.249 | -0.215 | -0.215 |
| FoGWar2 | ~~ | CoAWe8 | 10.24 | -0.198 | -0.198 | -0.159 | -0.159 |
| CoAWe7 | ~~ | CoAWe8 | 9.89 | 0.182 | 0.182 | 0.169 | 0.169 |
| FoGWar4 | ~~ | CoAWe7 | 7.16 | 0.160 | 0.160 | 0.160 | 0.160 |
The cross-loading of CoAWe5 onto l_FoGWar is clearly the most troubling problem in our model. If we include this cross-loading, the model \(\chi^2\) would improve by 27.59 points. We could also consider residual correlations among some of the FoGWar and CoAWe items if these were substantively motivated (e.g., specific similarities in content or wording). Whether we should dump CoAWe5 altogether (bad item) or allow it to cross-load is a substantive consideration that would depend on the goal of the analysis. Dropping the item would not cause problems with identification and would ease interpretability by maintaining simple structure. Here are the fits for the alternatives.
2.4.1 Drop CoAWe5
m4.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 3 items
l_FoGWar ~~ l_CoAWe
'
m4 <- sem(m4.syntax, simData)
inspect(m4)## $lambda
## l_FGWr l_CoAW
## FoGWar1 0 0
## FoGWar2 1 0
## FoGWar3 2 0
## FoGWar4 3 0
## CoAWe6 0 0
## CoAWe7 0 4
## CoAWe8 0 5
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8
## FoGWar1 7
## FoGWar2 0 8
## FoGWar3 0 0 9
## FoGWar4 0 0 0 10
## CoAWe6 0 0 0 0 11
## CoAWe7 0 0 0 0 0 12
## CoAWe8 0 0 0 0 0 0 13
##
## $psi
## l_FGWr l_CoAW
## l_FoGWar 14
## l_CoAWe 6 15semPaths_default(m4)
summary(m4, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 36 iterations
##
## Optimization method NLMINB
## Number of free parameters 15
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 39.972
## Degrees of freedom 13
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 701.420
## Degrees of freedom 21
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.960
## Tucker-Lewis Index (TLI) 0.936
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -5640.465
## Loglikelihood unrestricted model (H1) -5620.479
##
## Number of free parameters 15
## Akaike (AIC) 11310.930
## Bayesian (BIC) 11374.150
## Sample-size adjusted Bayesian (BIC) 11326.539
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.064
## 90 Percent Confidence Interval 0.042 0.088
## P-value RMSEA <= 0.05 0.133
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.042
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.784 0.592
## FoGWar2 0.905 0.104 8.721 0.000 0.709 0.523
## FoGWar3 0.906 0.101 8.934 0.000 0.710 0.541
## FoGWar4 1.436 0.137 10.449 0.000 1.126 0.747
## l_CoAWe =~
## CoAWe6 1.000 0.793 0.610
## CoAWe7 1.087 0.113 9.605 0.000 0.861 0.686
## CoAWe8 0.937 0.105 8.966 0.000 0.743 0.577
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.440 0.061 7.194 0.000 0.708 0.708
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.141 0.088 12.909 0.000 1.141 0.650
## .FoGWar2 1.336 0.097 13.817 0.000 1.336 0.727
## .FoGWar3 1.217 0.089 13.609 0.000 1.217 0.707
## .FoGWar4 1.005 0.110 9.096 0.000 1.005 0.442
## .CoAWe6 1.059 0.089 11.923 0.000 1.059 0.628
## .CoAWe7 0.837 0.084 9.998 0.000 0.837 0.530
## .CoAWe8 1.105 0.088 12.568 0.000 1.105 0.667
## l_FoGWar 0.614 0.099 6.230 0.000 1.000 1.000
## l_CoAWe 0.628 0.101 6.191 0.000 1.000 1.000We now see reasonable, though not perfect, fit on all practical measures. Thus, it appears that dropping this item largely resolves woes in the measurement model component.
2.4.2 Allow CoAWe5 cross-loading
m5.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe5 + CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 3 items
l_FoGWar =~ CoAWe5 #allow cross-loading
l_FoGWar ~~ l_CoAWe
'
m5 <- sem(m5.syntax, simData)
inspect(m5)## $lambda
## l_FGWr l_CoAW
## FoGWar1 0 0
## FoGWar2 1 0
## FoGWar3 2 0
## FoGWar4 3 0
## CoAWe5 7 0
## CoAWe6 0 4
## CoAWe7 0 5
## CoAWe8 0 6
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe5 CoAWe6 CoAWe7 CoAWe8
## FoGWar1 9
## FoGWar2 0 10
## FoGWar3 0 0 11
## FoGWar4 0 0 0 12
## CoAWe5 0 0 0 0 13
## CoAWe6 0 0 0 0 0 14
## CoAWe7 0 0 0 0 0 0 15
## CoAWe8 0 0 0 0 0 0 0 16
##
## $psi
## l_FGWr l_CoAW
## l_FoGWar 17
## l_CoAWe 8 18semPaths_default(m5)
summary(m5, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 44 iterations
##
## Optimization method NLMINB
## Number of free parameters 18
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 47.988
## Degrees of freedom 18
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 994.023
## Degrees of freedom 28
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.969
## Tucker-Lewis Index (TLI) 0.952
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -6391.466
## Loglikelihood unrestricted model (H1) -6367.472
##
## Number of free parameters 18
## Akaike (AIC) 12818.933
## Bayesian (BIC) 12894.796
## Sample-size adjusted Bayesian (BIC) 12837.663
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.058
## 90 Percent Confidence Interval 0.038 0.078
## P-value RMSEA <= 0.05 0.238
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.039
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.792 0.598
## FoGWar2 0.962 0.100 9.585 0.000 0.762 0.562
## FoGWar3 0.900 0.096 9.348 0.000 0.713 0.543
## FoGWar4 1.366 0.122 11.156 0.000 1.082 0.718
## l_CoAWe =~
## CoAWe5 1.000 0.487 0.337
## CoAWe6 1.642 0.384 4.271 0.000 0.800 0.616
## CoAWe7 1.742 0.407 4.284 0.000 0.849 0.675
## CoAWe8 1.539 0.362 4.249 0.000 0.750 0.583
## l_FoGWar =~
## CoAWe5 0.850 0.146 5.818 0.000 0.673 0.466
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.269 0.067 4.045 0.000 0.697 0.697
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.128 0.084 13.351 0.000 1.128 0.642
## .FoGWar2 1.257 0.091 13.779 0.000 1.257 0.684
## .FoGWar3 1.213 0.087 13.976 0.000 1.213 0.705
## .FoGWar4 1.101 0.100 11.002 0.000 1.101 0.485
## .CoAWe5 0.938 0.079 11.869 0.000 0.938 0.450
## .CoAWe6 1.047 0.084 12.387 0.000 1.047 0.620
## .CoAWe7 0.858 0.078 11.043 0.000 0.858 0.544
## .CoAWe8 1.093 0.084 12.952 0.000 1.093 0.660
## l_FoGWar 0.628 0.096 6.513 0.000 1.000 1.000
## l_CoAWe 0.237 0.105 2.264 0.024 1.000 1.000Allowing for this cross-loading also solves the fit problems with our measurement model. What to do?
For now, let’s enforce simple structure and drop CoAWe5 from further consideration.
2.5 Testing the structural model
Perhaps we also have problems in the structural model relating FoGWar and CoAWe to governmental spending and charitable donations. The informal way to handle this might be to retain factor scores from a well-fitting measurement model and conduct a simple path analysis to identify misfit. This somewhat flawed approach is still conceptually revealing since it reminds us that we are trying to assess structural relationships among variables, both latent and observed.
2.5.1 Using factor scores directly
Let’s extract the factor scores from the measurement model in which we dropped CoAWe5 (the problematic item).
m6.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 3 items
l_FoGWar ~~ l_CoAWe
'
m6 <- sem(m6.syntax, simData)
#summary(m6, fit.measures=TRUE)
fscores <- lavPredict(m6)
simData_augment <- bind_cols(simData, data.frame(fscores))
mstruct.syntax <- '
donations ~ l_FoGWar + l_CoAWe #predict criterion based on saved factors scores
gspending ~ donations #donations to NGOs predict governmental spending
'
mstruct <- sem(mstruct.syntax, simData_augment)
semPaths_default(mstruct)
summary(mstruct, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 16 iterations
##
## Optimization method NLMINB
## Number of free parameters 5
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 66.921
## Degrees of freedom 2
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 393.685
## Degrees of freedom 5
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.833
## Tucker-Lewis Index (TLI) 0.582
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1497.602
## Loglikelihood unrestricted model (H1) -1464.141
##
## Number of free parameters 5
## Akaike (AIC) 3005.204
## Bayesian (BIC) 3026.277
## Sample-size adjusted Bayesian (BIC) 3010.407
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.255
## 90 Percent Confidence Interval 0.205 0.309
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.109
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.421 0.126 3.349 0.001 0.421 0.253
## l_CoAWe 0.305 0.128 2.391 0.017 0.305 0.181
## gspending ~
## donations 0.766 0.045 17.159 0.000 0.766 0.609
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .donations 1.066 0.067 15.811 0.000 1.066 0.827
## .gspending 1.285 0.081 15.811 0.000 1.285 0.629We see from this analysis that the structural model (using saved factor scores) fits quite badly. These global fit statistics would never be considered acceptable and suggest that we cannot trust the structural parameter estimates. This is problematic since our hypotheses focused primarily on how fear of global warming and concern for animal welfare predict individual charitable contributions, which in turn influences governmental spending. We only have 2 \(df\) to play with, so what can we do?!
modificationIndices(mstruct)| lhs | op | rhs | mi | epc | sepc.lv | sepc.all | sepc.nox | |
|---|---|---|---|---|---|---|---|---|
| 9 | donations | ~~ | gspending | 62.46 | -0.993 | -0.993 | -0.849 | -0.849 |
| 10 | donations | ~ | gspending | 62.46 | -0.773 | -0.773 | -0.973 | -0.973 |
| 11 | gspending | ~ | l_FoGWar | 60.05 | 0.631 | 0.631 | 0.301 | 0.441 |
| 12 | gspending | ~ | l_CoAWe | 52.36 | 0.594 | 0.594 | 0.279 | 0.416 |
| 14 | l_FoGWar | ~ | gspending | 7.70 | 0.041 | 0.041 | 0.085 | 0.085 |
| 17 | l_CoAWe | ~ | gspending | 1.62 | 0.018 | 0.018 | 0.039 | 0.039 |
Watch out for those first two… That would yield a non-recursive model (bidirectional relationships between donations and gspending), which is rarely a good approach in SEM. But we do see a couple of more plausible candidates. As our path diagram reveals, we have assumed that the effects of FoGWar and CoAWe on governmental spending are entirely mediated by donations to charitable organizations. What if there are other reasons for the charity-government spending association? There could be unmeasured mediators, but our model enforces zero residual covariance of l_FoGWar and l_CoAWe with gspending after accounting for effects on charitable donations (donations).
Let’s consider the possibility that there is partial, but not full, mediation. To keep the structural model in overidentified territory, we can free just the largest MI and examine global fit.
mstruct2.syntax <- '
donations ~ l_FoGWar + l_CoAWe #predict criterion based on saved factors scores
gspending ~ donations #donations to NGOs predict governmental spending
gspending ~ l_CoAWe #allow direct effect of CoAWe on governmental spending
'
mstruct2 <- sem(mstruct2.syntax, simData_augment)
semPaths_default(mstruct2)
summary(mstruct2, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 19 iterations
##
## Optimization method NLMINB
## Number of free parameters 6
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 11.616
## Degrees of freedom 1
## P-value (Chi-square) 0.001
##
## Model test baseline model:
##
## Minimum Function Test Statistic 393.685
## Degrees of freedom 5
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.973
## Tucker-Lewis Index (TLI) 0.863
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1469.949
## Loglikelihood unrestricted model (H1) -1464.141
##
## Number of free parameters 6
## Akaike (AIC) 2951.899
## Bayesian (BIC) 2977.186
## Sample-size adjusted Bayesian (BIC) 2958.142
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.146
## 90 Percent Confidence Interval 0.079 0.226
## P-value RMSEA <= 0.05 0.011
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.019
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.421 0.126 3.349 0.001 0.421 0.253
## l_CoAWe 0.305 0.128 2.391 0.017 0.305 0.181
## gspending ~
## donations 0.628 0.046 13.661 0.000 0.628 0.499
## l_CoAWe 0.594 0.078 7.647 0.000 0.594 0.279
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .donations 1.066 0.067 15.811 0.000 1.066 0.827
## .gspending 1.151 0.073 15.811 0.000 1.151 0.563We have moved in the right direction, with acceptable values of SRMR and CFI, but not TLI or RMSEA. Let’s examine the just-identified model in which both mediators have a direct effect:
mstruct3.syntax <- '
donations ~ l_FoGWar + l_CoAWe #predict criterion based on saved factors scores
gspending ~ donations #donations to NGOs predict governmental spending
gspending ~ l_FoGWar + l_CoAWe #allow direct effects of both FoGWar and CoAWe on governmental spending
'
mstruct3 <- sem(mstruct3.syntax, simData_augment)
semPaths_default(mstruct3)
summary(mstruct3, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 22 iterations
##
## Optimization method NLMINB
## Number of free parameters 7
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 0.000
## Degrees of freedom 0
## Minimum Function Value 0.0000000000000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 393.685
## Degrees of freedom 5
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1464.141
## Loglikelihood unrestricted model (H1) -1464.141
##
## Number of free parameters 7
## Akaike (AIC) 2942.283
## Bayesian (BIC) 2971.785
## Sample-size adjusted Bayesian (BIC) 2949.567
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.000
## P-value RMSEA <= 0.05 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## donations ~
## l_FoGWar 0.421 0.126 3.349 0.001 0.421 0.253
## l_CoAWe 0.305 0.128 2.391 0.017 0.305 0.181
## gspending ~
## donations 0.605 0.046 13.161 0.000 0.605 0.480
## l_FoGWar 0.448 0.131 3.428 0.001 0.448 0.214
## l_CoAWe 0.227 0.132 1.719 0.086 0.227 0.107
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .donations 1.066 0.067 15.811 0.000 1.066 0.827
## .gspending 1.124 0.071 15.811 0.000 1.124 0.551By definition, we have achieved perfect fit, but notice that both FoGWar and CoAWe have strong relationships with governmental spending after accounting for charitable donations. We should not ignore these findings if we want the overall SEM to fit well and to trust the mediated effects of interest. Thus, we should strongly consider the saturated structural model.
2.5.2 Formal testing of the structural model within the full SEM
As I mentioned above, extracting the factors scores is intuitive, and in some ways sensible, but it does not follow the SEM principle of trying to estimate all important effects, both measurement and structural, in a single model. In particular, we lose any information about covariances among the indicators of CoAWe and FoGWar with the outcomes of charitable donations and governemental spending. If there are strong residual associations after accounting for the factor model, we would entirely miss these by extracting fator scores.
Thus, to formally test the local fit of our candidate structural model (full mediation), we can leave the measurement model alone (drop CoAWe5, as above), but test an alternative in which we saturate the structural model only. That is, we can consider how well a model fits when l_FoGWar, l_CoAWe, donations, and gspending are allowed to correlate freely with each other.
This leads to a tricky situation, however, because in SEM (at least conventional LISREL all-y notation), there is no matrix that represents the correlations among observed and latent variables. Rather, \(\boldsymbol{\Theta}_\varepsilon\) contains the (residual) covariance structure for the observed variables and \(\boldsymbol{\Psi}\) contains the covariance structure of the latent variables.
When we regress observed variables (e.g., gspending) on latent variables (e.g., l_FoGWar) above, lavaan converts gspending into a single-indicator latent variable without our awareness. This allows the residual variance of gspending to be incorporated into \(\boldsymbol{\Psi}\) while the regression coefficient is in \(\boldsymbol{B}\). Here, we need to create the single-indicator latent variable ourselves. This is done by defining a latent variable with a fixed 1.0 loading onto the observed variable. Note that this does not increase our number of free parameters because the loading is fixed. It is more of a trick to get the matrix math to work properly.
mstruct_saturated.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 3 items
l_gspending =~ 1*gspending #could omit 1* since that is the default; left here for clarity
l_donations =~ 1*donations
#saturate structural model with undirected relationships
#4 choose 2 = 6 parameters
l_donations ~~ l_FoGWar
l_donations ~~ l_CoAWe
l_donations ~~ l_gspending
l_gspending ~~ l_FoGWar
l_gspending ~~ l_CoAWe
l_FoGWar ~~ l_CoAWe
'
mstruct_saturated <- sem(mstruct_saturated.syntax, simData)
semPaths_default(mstruct_saturated, sizeMan=10)
inspect(mstruct_saturated)## $lambda
## l_FGWr l_CoAW l_gspn l_dntn
## FoGWar1 0 0 0 0
## FoGWar2 1 0 0 0
## FoGWar3 2 0 0 0
## FoGWar4 3 0 0 0
## CoAWe6 0 0 0 0
## CoAWe7 0 4 0 0
## CoAWe8 0 5 0 0
## gspending 0 0 0 0
## donations 0 0 0 0
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 12
## FoGWar2 0 13
## FoGWar3 0 0 14
## FoGWar4 0 0 0 15
## CoAWe6 0 0 0 0 16
## CoAWe7 0 0 0 0 0 17
## CoAWe8 0 0 0 0 0 0 18
## gspending 0 0 0 0 0 0 0 0
## donations 0 0 0 0 0 0 0 0 0
##
## $psi
## l_FGWr l_CoAW l_gspn l_dntn
## l_FoGWar 19
## l_CoAWe 11 20
## l_gspending 9 10 21
## l_donations 6 7 8 22Notice how the $psi matrix is constructed.
What about fit?
summary(mstruct_saturated, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 38 iterations
##
## Optimization method NLMINB
## Number of free parameters 22
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 48.342
## Degrees of freedom 23
## P-value (Chi-square) 0.002
##
## Model test baseline model:
##
## Minimum Function Test Statistic 1103.651
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.976
## Tucker-Lewis Index (TLI) 0.963
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7104.518
## Loglikelihood unrestricted model (H1) -7080.347
##
## Number of free parameters 22
## Akaike (AIC) 14253.036
## Bayesian (BIC) 14345.758
## Sample-size adjusted Bayesian (BIC) 14275.928
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.047
## 90 Percent Confidence Interval 0.028 0.065
## P-value RMSEA <= 0.05 0.580
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.035
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.774 0.584
## FoGWar2 0.932 0.104 8.952 0.000 0.722 0.532
## FoGWar3 0.930 0.102 9.148 0.000 0.720 0.549
## FoGWar4 1.444 0.134 10.787 0.000 1.118 0.742
## l_CoAWe =~
## CoAWe6 1.000 0.794 0.611
## CoAWe7 1.083 0.109 9.977 0.000 0.860 0.684
## CoAWe8 0.936 0.102 9.197 0.000 0.743 0.578
## l_gspending =~
## gspending 1.000 1.429 1.000
## l_donations =~
## donations 1.000 1.136 1.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_donations 0.388 0.056 6.935 0.000 0.442 0.442
## l_CoAWe ~~
## l_donations 0.378 0.057 6.573 0.000 0.419 0.419
## l_gspending ~~
## l_donations 0.988 0.085 11.628 0.000 0.609 0.609
## l_FoGWar ~~
## l_gspending 0.606 0.076 7.956 0.000 0.548 0.548
## l_CoAWe ~~
## l_gspending 0.568 0.076 7.454 0.000 0.501 0.501
## l_FoGWar ~~
## l_CoAWe 0.434 0.060 7.267 0.000 0.706 0.706
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.157 0.087 13.325 0.000 1.157 0.659
## .FoGWar2 1.317 0.095 13.930 0.000 1.317 0.717
## .FoGWar3 1.203 0.087 13.755 0.000 1.203 0.699
## .FoGWar4 1.023 0.103 9.882 0.000 1.023 0.450
## .CoAWe6 1.056 0.087 12.194 0.000 1.056 0.626
## .CoAWe7 0.840 0.081 10.432 0.000 0.840 0.532
## .CoAWe8 1.103 0.086 12.798 0.000 1.103 0.666
## .gspending 0.000 0.000 0.000
## .donations 0.000 0.000 0.000
## l_FoGWar 0.599 0.096 6.259 0.000 1.000 1.000
## l_CoAWe 0.630 0.100 6.317 0.000 1.000 1.000
## l_gspending 2.042 0.129 15.811 0.000 1.000 1.000
## l_donations 1.289 0.082 15.811 0.000 1.000 1.000It’s not perfect, but looks pretty good – all practical fit indices are in acceptable ranges. Thus, in the full SEM, saturating the structural model and fixing up the measurement model yields an acceptable solution.
The crucial test is how this model fits with our a priori hypothesis of mediation. We can use an LRT to test this because the restricted model is nested within the saturated (structural) model. To get the LRT right (i.e., exactly comparable models), we should use the single-indicator approach to the original model explicitly, rather than allowing lavaan to change the parameters behind the scenes.
msem.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_gspending =~ 1*gspending #could omit 1* since that is the default; left here for clarity
l_donations =~ 1*donations
l_donations ~ l_FoGWar + l_CoAWe #predict criterion based on SEM
l_gspending ~ l_donations #donations predict governmental spending
l_FoGWar ~~ l_CoAWe #correlated factors
'
msem <- sem(msem.syntax, simData)
semPaths_default(msem)
inspect(msem)## $lambda
## l_FGWr l_CoAW l_gspn l_dntn
## FoGWar1 0 0 0 0
## FoGWar2 1 0 0 0
## FoGWar3 2 0 0 0
## FoGWar4 3 0 0 0
## CoAWe6 0 0 0 0
## CoAWe7 0 4 0 0
## CoAWe8 0 5 0 0
## gspending 0 0 0 0
## donations 0 0 0 0
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 10
## FoGWar2 0 11
## FoGWar3 0 0 12
## FoGWar4 0 0 0 13
## CoAWe6 0 0 0 0 14
## CoAWe7 0 0 0 0 0 15
## CoAWe8 0 0 0 0 0 0 16
## gspending 0 0 0 0 0 0 0 0
## donations 0 0 0 0 0 0 0 0 0
##
## $psi
## l_FGWr l_CoAW l_gspn l_dntn
## l_FoGWar 17
## l_CoAWe 9 18
## l_gspending 0 0 19
## l_donations 0 0 0 20
##
## $beta
## l_FGWr l_CoAW l_gspn l_dntn
## l_FoGWar 0 0 0 0
## l_CoAWe 0 0 0 0
## l_gspending 0 0 0 8
## l_donations 6 7 0 0anova(mstruct_saturated, msem)| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| mstruct_saturated | 23 | 14253 | 14346 | 48.3 | NA | NA | NA |
| msem | 25 | 14316 | 14400 | 115.3 | 67 | 2 | 0 |
We see that the LRT is highly significant and the AIC and BIC for the saturated model is much better than the restricted mediation model. Thus, the hypothesis that charitable donations fully mediates the relationship of FoGWar and CoAWe with governmental spending is not supported.
What about the alternative that the mediation can be partial, and that FoGWar and CoAWe also have direct effects?
msem_partial.syntax <- '
l_FoGWar =~ FoGWar1 + FoGWar2 + FoGWar3 + FoGWar4 #definition of factor FoGWar with loadings on 4 items
l_CoAWe =~ CoAWe6 + CoAWe7 + CoAWe8 #definition of factor CoAWe with loadings on 4 items
l_gspending =~ 1*gspending #could omit 1* since that is the default; left here for clarity
l_donations =~ 1*donations
l_donations ~ l_FoGWar + l_CoAWe #predict criterion based on SEM
l_gspending ~ l_donations + l_FoGWar + l_CoAWe #all three predict spending
l_FoGWar ~~ l_CoAWe #correlated factors
'
msem_partial <- sem(msem_partial.syntax, simData)
semPaths_default(msem_partial)
inspect(msem_partial)## $lambda
## l_FGWr l_CoAW l_gspn l_dntn
## FoGWar1 0 0 0 0
## FoGWar2 1 0 0 0
## FoGWar3 2 0 0 0
## FoGWar4 3 0 0 0
## CoAWe6 0 0 0 0
## CoAWe7 0 4 0 0
## CoAWe8 0 5 0 0
## gspending 0 0 0 0
## donations 0 0 0 0
##
## $theta
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 12
## FoGWar2 0 13
## FoGWar3 0 0 14
## FoGWar4 0 0 0 15
## CoAWe6 0 0 0 0 16
## CoAWe7 0 0 0 0 0 17
## CoAWe8 0 0 0 0 0 0 18
## gspending 0 0 0 0 0 0 0 0
## donations 0 0 0 0 0 0 0 0 0
##
## $psi
## l_FGWr l_CoAW l_gspn l_dntn
## l_FoGWar 19
## l_CoAWe 11 20
## l_gspending 0 0 21
## l_donations 0 0 0 22
##
## $beta
## l_FGWr l_CoAW l_gspn l_dntn
## l_FoGWar 0 0 0 0
## l_CoAWe 0 0 0 0
## l_gspending 9 10 0 8
## l_donations 6 7 0 0summary(msem_partial, fit.measures=TRUE, standardized=TRUE)## lavaan 0.6-3 ended normally after 49 iterations
##
## Optimization method NLMINB
## Number of free parameters 22
##
## Number of observations 500
##
## Estimator ML
## Model Fit Test Statistic 48.342
## Degrees of freedom 23
## P-value (Chi-square) 0.002
##
## Model test baseline model:
##
## Minimum Function Test Statistic 1103.651
## Degrees of freedom 36
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.976
## Tucker-Lewis Index (TLI) 0.963
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -7104.518
## Loglikelihood unrestricted model (H1) -7080.347
##
## Number of free parameters 22
## Akaike (AIC) 14253.036
## Bayesian (BIC) 14345.758
## Sample-size adjusted Bayesian (BIC) 14275.928
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.047
## 90 Percent Confidence Interval 0.028 0.065
## P-value RMSEA <= 0.05 0.580
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.035
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar =~
## FoGWar1 1.000 0.774 0.584
## FoGWar2 0.932 0.104 8.952 0.000 0.722 0.532
## FoGWar3 0.930 0.102 9.148 0.000 0.720 0.549
## FoGWar4 1.444 0.134 10.787 0.000 1.118 0.742
## l_CoAWe =~
## CoAWe6 1.000 0.794 0.611
## CoAWe7 1.083 0.109 9.977 0.000 0.860 0.684
## CoAWe8 0.936 0.102 9.197 0.000 0.743 0.578
## l_gspending =~
## gspending 1.000 1.429 1.000
## l_donations =~
## donations 1.000 1.136 1.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_donations ~
## l_FoGWar 0.427 0.132 3.223 0.001 0.291 0.291
## l_CoAWe 0.305 0.131 2.336 0.020 0.214 0.214
## l_gspending ~
## l_donations 0.551 0.049 11.142 0.000 0.438 0.438
## l_FoGWar 0.481 0.139 3.452 0.001 0.261 0.261
## l_CoAWe 0.240 0.135 1.781 0.075 0.133 0.133
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## l_FoGWar ~~
## l_CoAWe 0.434 0.060 7.267 0.000 0.706 0.706
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .FoGWar1 1.157 0.087 13.325 0.000 1.157 0.659
## .FoGWar2 1.317 0.095 13.930 0.000 1.317 0.717
## .FoGWar3 1.203 0.087 13.755 0.000 1.203 0.699
## .FoGWar4 1.023 0.103 9.882 0.000 1.023 0.450
## .CoAWe6 1.056 0.087 12.194 0.000 1.056 0.626
## .CoAWe7 0.840 0.081 10.432 0.000 0.840 0.532
## .CoAWe8 1.103 0.086 12.798 0.000 1.103 0.666
## .gspending 0.000 0.000 0.000
## .donations 0.000 0.000 0.000
## l_FoGWar 0.599 0.096 6.259 0.000 1.000 1.000
## l_CoAWe 0.630 0.100 6.317 0.000 1.000 1.000
## .l_gspending 1.070 0.072 14.797 0.000 0.524 0.524
## .l_donations 1.008 0.068 14.744 0.000 0.782 0.782Hmm, looks pretty gosh-darned similar to the model with a saturated structural component.
anova(msem_partial, mstruct_saturated)## Warning in lavTestLRT(object = new("lavaan", version = "0.6.3", call =
## lavaan::lavaan(model = msem_partial.syntax, : lavaan WARNING: some models
## have the same degrees of freedom| Df | AIC | BIC | Chisq | Chisq diff | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|
| msem_partial | 23 | 14253 | 14346 | 48.3 | NA | NA | NA |
| mstruct_saturated | 23 | 14253 | 14346 | 48.3 | 0 | 0 | NA |
Ignore the p-value for a second and just notice the identical df and AIC. Likewise, here are the log-likelihood estimates for each:
Partial mediation model:
logLik(msem_partial)## 'log Lik.' -7105 (df=22)Saturated structural model:
logLik(mstruct_saturated)## 'log Lik.' -7105 (df=22)These are equivalent models, a topic we will spend more time on next week. For now, know that replacing a covariance term with a regression term yields identical fit. And recall that covariance is undirected. Thus, in our candidate model that allows for partial mediation, we believe that the effects are directed, but because there are strong associations among all 4 variables (FoGWar, CoAWe, gspending, and donations), we cannot constrain the model to be more specific. Therefore, there is no statistical evidence for our view of partial mediation relative to a model in which all structural relationships are undirected.
2.6 Enough of fit statistics already, how do I know where my model fails?
As noted above, modification indices provide direct feedback about potential omitted effects in the model. That is, they help to identify relationships (measurement or structural) that, if included, would improve model fit (in \(\chi^2\) terms).
However, examination of residuals and standardized residuals can provide additional, complementary, guidance about sources of misfit in the model. Recall that the goal of SEM is to provide a good representation of the covariance matrix with the fewest parameters possible. Thus, residuals are conceptualized in terms of observed vs. predicted values for specific covariance terms, as in RMR and SRMR.
2.6.1 Raw Residuals
We can obtain raw residuals using the resid() function. In our earlier SEM with a saturated structural model and no CoAWe5, we would get:
resid(mstruct_saturated)## $type
## [1] "raw"
##
## $cov
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 0.000
## FoGWar2 0.130 0.000
## FoGWar3 0.063 0.102 0.000
## FoGWar4 -0.015 -0.098 -0.057 0.000
## CoAWe6 -0.096 0.008 -0.037 0.124 0.000
## CoAWe7 -0.052 -0.083 -0.049 0.155 -0.040 0.000
## CoAWe8 -0.093 -0.158 -0.111 0.059 0.034 0.019 0.000
## gspending -0.068 0.035 0.052 -0.005 0.007 0.015 -0.031 0.000
## donations -0.041 0.013 -0.055 0.048 0.007 -0.015 0.015 0.000 0.000Note that we get a vector of means (all zero), even though we didn’t fit mean structure. The values of each of these reflect the difference between the model-implied covariance matrix, \(\boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})\), and the sample covariance matrix, \(\boldsymbol{S_{XX}}\).
The model-implied covariance is given by fitted():
sigma_hat <- fitted(mstruct_saturated)$cov
print(sigma_hat)## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 1.756
## FoGWar2 0.559 1.838
## FoGWar3 0.557 0.520 1.721
## FoGWar4 0.865 0.807 0.805 2.272
## CoAWe6 0.434 0.405 0.404 0.627 1.687
## CoAWe7 0.470 0.438 0.437 0.679 0.683 1.579
## CoAWe8 0.407 0.379 0.378 0.587 0.590 0.639 1.656
## gspending 0.606 0.565 0.564 0.876 0.568 0.615 0.532 2.042
## donations 0.388 0.362 0.361 0.561 0.378 0.409 0.354 0.988 1.289And the sample covariance is given by cov():
s_xx <- cov(dplyr::select(simData, -CoAWe5))
print(s_xx)## FoGWar1 FoGWar2 FoGWar3 FoGWar4 CoAWe6 CoAWe7 CoAWe8 donations
## FoGWar1 1.759 0.690 0.622 0.852 0.339 0.419 0.314 0.348
## FoGWar2 0.690 1.842 0.623 0.710 0.414 0.356 0.222 0.376
## FoGWar3 0.622 0.623 1.725 0.750 0.368 0.389 0.267 0.307
## FoGWar4 0.852 0.710 0.750 2.277 0.753 0.835 0.647 0.610
## CoAWe6 0.339 0.414 0.368 0.753 1.690 0.644 0.626 0.385
## CoAWe7 0.419 0.356 0.389 0.835 0.644 1.582 0.660 0.395
## CoAWe8 0.314 0.222 0.267 0.647 0.626 0.660 1.659 0.370
## donations 0.348 0.376 0.307 0.610 0.385 0.395 0.370 1.292
## gspending 0.540 0.602 0.618 0.873 0.577 0.632 0.502 0.990
## gspending
## FoGWar1 0.540
## FoGWar2 0.602
## FoGWar3 0.618
## FoGWar4 0.873
## CoAWe6 0.577
## CoAWe7 0.632
## CoAWe8 0.502
## donations 0.990
## gspending 2.046Notice how the column order is not necessarily the same between these? This can be an irritation. Here’s a bit of code to order the model-implied covariance matrix identically to the sample covariance matrix.
print(sigma_hat[rownames(s_xx), colnames(s_xx)])## FoGWar1 FoGWar2 FoGWar3 FoGWar4 CoAWe6 CoAWe7 CoAWe8 donations
## FoGWar1 1.756 0.559 0.557 0.865 0.434 0.470 0.407 0.388
## FoGWar2 0.559 1.838 0.520 0.807 0.405 0.438 0.379 0.362
## FoGWar3 0.557 0.520 1.721 0.805 0.404 0.437 0.378 0.361
## FoGWar4 0.865 0.807 0.805 2.272 0.627 0.679 0.587 0.561
## CoAWe6 0.434 0.405 0.404 0.627 1.687 0.683 0.590 0.378
## CoAWe7 0.470 0.438 0.437 0.679 0.683 1.579 0.639 0.409
## CoAWe8 0.407 0.379 0.378 0.587 0.590 0.639 1.656 0.354
## donations 0.388 0.362 0.361 0.561 0.378 0.409 0.354 1.289
## gspending 0.606 0.565 0.564 0.876 0.568 0.615 0.532 0.988
## gspending
## FoGWar1 0.606
## FoGWar2 0.565
## FoGWar3 0.564
## FoGWar4 0.876
## CoAWe6 0.568
## CoAWe7 0.615
## CoAWe8 0.532
## donations 0.988
## gspending 2.0422.6.2 Standardized residuals
The raw residuals are always worth examination, but sometimes if variables vary in scale, or the variances are far from 1, the covariances are hard to interpret. How much is too much in the residuals?
This is where standardized residuals play a useful role in diagnosing misfit. Standardization is performed by dividing the raw residuals the square root of the difference between the model-implied and observed statistics. Long story short, they are in z-score terms, but if you see NA values, you can safely ignore them. We simply modify our resid() call slightly to obtain these:
resid(mstruct_saturated, type="standardized")## $type
## [1] "standardized"
##
## $cov
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 0.000
## FoGWar2 2.723 0.000
## FoGWar3 1.443 2.091 0.000
## FoGWar4 -0.512 -3.157 -1.887 0.000
## CoAWe6 -1.839 0.148 -0.686 2.595 0.000
## CoAWe7 -1.120 -1.665 -1.031 3.675 -1.882 0.000
## CoAWe8 -1.716 -2.731 -2.013 1.195 1.014 0.757 0.000
## gspending -1.555 0.691 1.085 -0.136 0.178 0.458 -0.716 0.000
## donations -1.037 0.298 -1.295 1.455 0.180 -0.513 0.394 0.000 0.000Large positive values (e.g., > 2) indicate that the model is underestimating the covariance term, whereas large negative values (e.g., < -2) indicate overestimation. This follows from regression, where residuals are \(Y - \hat{Y}\). Thus, standardized residuals given us direct information about which covariance terms are most misestimated.
When examining standardized residuals, remember that, as in null hypothesis significance testing, there is a chance of a false positive. That is, if you have 20 standardized residuals, you expect at least one of them to be ‘significant’ (\(|z| > 1.96\)) by chance even if the model is correct.
2.6.3 Residuals in correlation units
We can also simply convert raw residuals into correlational units (standardize variances to 1) if that is more intuitive. These units can give us an sense of misfit on a common scale, where we might consider values \(|r| > .10\) as troubling a priori (just a rule of thumb).
resid(mstruct_saturated, type="cor")## $type
## [1] "cor.bollen"
##
## $cov
## FoGWr1 FoGWr2 FoGWr3 FoGWr4 CoAWe6 CoAWe7 CoAWe8 gspndn dontns
## FoGWar1 0.000
## FoGWar2 0.072 0.000
## FoGWar3 0.036 0.058 0.000
## FoGWar4 -0.007 -0.048 -0.029 0.000
## CoAWe6 -0.056 0.005 -0.022 0.064 0.000
## CoAWe7 -0.031 -0.049 -0.030 0.082 -0.024 0.000
## CoAWe8 -0.054 -0.091 -0.066 0.030 0.021 0.012 0.000
## gspending -0.036 0.018 0.028 -0.002 0.004 0.008 -0.017 0.000
## donations -0.027 0.009 -0.037 0.028 0.004 -0.010 0.010 0.000 0.000Here’s a handy trick for finding only those residual correlations above a given threshold, here |r| > .05:
rmat <- resid(mstruct_saturated, type="cor")$cov
rmat[upper.tri(rmat)] <- NA #make sure that we don't double count the upper and lower triangles of the symmetric matrix
rmat %>% reshape2::melt() %>% filter(abs(value) > .05)| Var1 | Var2 | value |
|---|---|---|
| FoGWar2 | FoGWar1 | 0.072 |
| CoAWe6 | FoGWar1 | -0.056 |
| CoAWe8 | FoGWar1 | -0.054 |
| FoGWar3 | FoGWar2 | 0.058 |
| CoAWe8 | FoGWar2 | -0.091 |
| CoAWe8 | FoGWar3 | -0.066 |
| CoAWe6 | FoGWar4 | 0.064 |
| CoAWe7 | FoGWar4 | 0.082 |
Finally, examination of residuals can help to guide your search through the modification indices to see which may be both justifiable and most related to covariance misestimation that is central to your substantive question.
For more details on residuals, see: https://www.statmodel.com/download/StandardizedResiduals.pdf.