SEM madness with the Wheaton data
Dr. Thomas Pollet, Northumbria University (thomas.pollet@northumbria.ac.uk)
2022-07-25 | disclaimer
Hypothesised model
Data
require(lavaan)## Loading required package: lavaan## This is lavaan 0.6-11
## lavaan is FREE software! Please report any bugs.lower <- '
11.834,
6.947, 9.364,
6.819, 5.091, 12.532,
4.783, 5.028, 7.495, 9.986,
-3.839, -3.889, -3.841, -3.625, 9.610,
-21.899, -18.831, -21.748, -18.775, 35.522, 450.288 '
# convert to a full symmetric covariance matrix with names
wheaton.cov <- getCov(lower, names=c("anomia67","powerless67", "anomia71",
"powerless71","education","sei"))Model
require(lavaan)
# the model
wheaton.model = '
# measurement model
ses =~ education + sei
alien67 =~ anomia67 + powerless67
alien71 =~ anomia71 + powerless71
# structural model
alien71 ~ aa*alien67 + ses
alien67 ~ sa*ses
# correlated residuals
anomia67 ~~ anomia71
powerless67 ~~ powerless71
# Indirect effect
IndirectEffect := sa*aa
'Fit the model.
require(lavaan)
alienation <- sem(wheaton.model, sample.cov=wheaton.cov, sample.nobs=932)
summary(alienation, fit.measures=T)## lavaan 0.6-11 ended normally after 84 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 17
##
## Number of observations 932
##
## Model Test User Model:
##
## Test statistic 4.735
## Degrees of freedom 4
## P-value (Chi-square) 0.316
##
## Model Test Baseline Model:
##
## Test statistic 2133.722
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 0.999
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15213.274
## Loglikelihood unrestricted model (H1) -15210.906
##
## Akaike (AIC) 30460.548
## Bayesian (BIC) 30542.783
## Sample-size adjusted Bayesian (BIC) 30488.792
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.014
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.053
## P-value RMSEA <= 0.05 0.930
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.007
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ses =~
## education 1.000
## sei 5.219 0.422 12.364 0.000
## alien67 =~
## anomia67 1.000
## powerless67 0.979 0.062 15.895 0.000
## alien71 =~
## anomia71 1.000
## powerless71 0.922 0.059 15.498 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## alien71 ~
## alien67 (aa) 0.607 0.051 11.898 0.000
## ses -0.227 0.052 -4.334 0.000
## alien67 ~
## ses (sa) -0.575 0.056 -10.195 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .anomia67 ~~
## .anomia71 1.623 0.314 5.176 0.000
## .powerless67 ~~
## .powerless71 0.339 0.261 1.298 0.194
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .education 2.801 0.507 5.525 0.000
## .sei 264.597 18.126 14.597 0.000
## .anomia67 4.731 0.453 10.441 0.000
## .powerless67 2.563 0.403 6.359 0.000
## .anomia71 4.399 0.515 8.542 0.000
## .powerless71 3.070 0.434 7.070 0.000
## ses 6.798 0.649 10.475 0.000
## .alien67 4.841 0.467 10.359 0.000
## .alien71 4.083 0.404 10.104 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## IndirectEffect -0.349 0.041 -8.538 0.000Plot
require(semPlot)
require(qgraph)
semPaths(alienation,style = "lisrel",what="std", curvePivot = TRUE)
Model without the arrow?
This model drops the arrow.
require(lavaan)
# the model
wheaton.model2 = '
# measurement model
ses =~ education + sei
alien67 =~ anomia67 + powerless67
alien71 =~ anomia71 + powerless71
# structural model
alien71 ~ ses
alien67 ~ sa*ses
# correlated residuals
anomia67 ~~ anomia71
powerless67 ~~ powerless71
'You’ll notice that the fit statistics are identical?!
Why is that? We have modelled the covariance now instead of the arrow! In short, we cannot differentiate between double-headed arrow and directed arrow in terms of fit! It is a decision by the researcher to have one or the other but they are equivalent models – Read more in the sources cited.
require(lavaan)
alienation_without <- sem(wheaton.model2, sample.cov=wheaton.cov, sample.nobs=932)
summary(alienation_without, fit.measures=T)## lavaan 0.6-11 ended normally after 82 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 17
##
## Number of observations 932
##
## Model Test User Model:
##
## Test statistic 4.735
## Degrees of freedom 4
## P-value (Chi-square) 0.316
##
## Model Test Baseline Model:
##
## Test statistic 2133.722
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 0.999
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15213.274
## Loglikelihood unrestricted model (H1) -15210.906
##
## Akaike (AIC) 30460.548
## Bayesian (BIC) 30542.783
## Sample-size adjusted Bayesian (BIC) 30488.792
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.014
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.053
## P-value RMSEA <= 0.05 0.930
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.007
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ses =~
## education 1.000
## sei 5.219 0.422 12.364 0.000
## alien67 =~
## anomia67 1.000
## powerless67 0.979 0.062 15.895 0.000
## alien71 =~
## anomia71 1.000
## powerless71 0.922 0.059 15.498 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## alien71 ~
## ses -0.576 0.058 -9.956 0.000
## alien67 ~
## ses (sa) -0.575 0.056 -10.195 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .anomia67 ~~
## .anomia71 1.623 0.314 5.176 0.000
## .powerless67 ~~
## .powerless71 0.339 0.261 1.298 0.194
## .alien67 ~~
## .alien71 2.939 0.351 8.380 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .education 2.801 0.507 5.525 0.000
## .sei 264.597 18.126 14.597 0.000
## .anomia67 4.731 0.453 10.441 0.000
## .powerless67 2.563 0.403 6.359 0.000
## .anomia71 4.399 0.515 8.542 0.000
## .powerless71 3.070 0.434 7.070 0.000
## ses 6.798 0.649 10.475 0.000
## .alien67 4.841 0.467 10.359 0.000
## .alien71 5.867 0.554 10.585 0.000require(semPlot)
require(qgraph)
semPaths(alienation_without,style = "lisrel",what="std", curvePivot = TRUE)
Model without covariance between latent constructs
require(lavaan)
# the model
wheaton.model3 = '
# measurement model
ses =~ education + sei
alien67 =~ anomia67 + powerless67
alien71 =~ anomia71 + powerless71
# structural model
alien71 ~ ses
alien67 ~ sa*ses
# Orthogonal (not correlated)
alien71 ~~ 0*alien67
# correlated residuals
anomia67 ~~ anomia71
powerless67 ~~ powerless71
'require(lavaan)
alienation_without2 <- sem(wheaton.model3, sample.cov=wheaton.cov, sample.nobs=932)
summary(alienation_without2, fit.measures=T)## lavaan 0.6-11 ended normally after 72 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 932
##
## Model Test User Model:
##
## Test statistic 113.883
## Degrees of freedom 5
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 2133.722
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.949
## Tucker-Lewis Index (TLI) 0.846
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15267.848
## Loglikelihood unrestricted model (H1) -15210.906
##
## Akaike (AIC) 30567.696
## Bayesian (BIC) 30645.093
## Sample-size adjusted Bayesian (BIC) 30594.278
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.153
## 90 Percent confidence interval - lower 0.129
## 90 Percent confidence interval - upper 0.178
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.061
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ses =~
## education 1.000
## sei 5.821 0.411 14.172 0.000
## alien67 =~
## anomia67 1.000
## powerless67 0.943 0.055 17.235 0.000
## alien71 =~
## anomia71 1.000
## powerless71 0.895 0.053 16.944 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## alien71 ~
## ses -0.926 0.068 -13.703 0.000
## alien67 ~
## ses (sa) -0.908 0.066 -13.676 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .alien67 ~~
## .alien71 0.000
## .anomia67 ~~
## .anomia71 1.755 0.254 6.905 0.000
## .powerless67 ~~
## .powerless71 0.811 0.199 4.083 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .education 4.731 0.346 13.689 0.000
## .sei 284.832 16.584 17.175 0.000
## .anomia67 4.595 0.412 11.153 0.000
## .powerless67 3.039 0.344 8.839 0.000
## .anomia71 4.336 0.460 9.435 0.000
## .powerless71 3.459 0.369 9.371 0.000
## ses 4.869 0.470 10.354 0.000
## .alien67 2.995 0.349 8.576 0.000
## .alien71 3.797 0.416 9.118 0.000require(semPlot)
require(qgraph)
semPaths(alienation_without2,style = "lisrel",what="std", curvePivot = TRUE)
This model is clearly inferior in terms of fit and also in absolute terms no longer a good fit in RMSEA.
Fit the model.
require(lavaan)
alienation_no_corr_resid <- sem(wheaton.model_no_corr_resid, sample.cov=wheaton.cov, sample.nobs=932)
summary(alienation_no_corr_resid, fit.measures=T)## lavaan 0.6-11 ended normally after 73 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
##
## Number of observations 932
##
## Model Test User Model:
##
## Test statistic 71.546
## Degrees of freedom 6
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 2133.722
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.969
## Tucker-Lewis Index (TLI) 0.923
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15246.680
## Loglikelihood unrestricted model (H1) -15210.906
##
## Akaike (AIC) 30523.359
## Bayesian (BIC) 30595.919
## Sample-size adjusted Bayesian (BIC) 30548.281
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.108
## 90 Percent confidence interval - lower 0.087
## 90 Percent confidence interval - upper 0.131
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.021
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ses =~
## education 1.000
## sei 5.329 0.430 12.406 0.000
## alien67 =~
## anomia67 1.000
## powerless67 0.889 0.041 21.422 0.000
## alien71 =~
## anomia71 1.000
## powerless71 0.849 0.040 21.252 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## alien71 ~
## alien67 (aa) 0.705 0.054 13.170 0.000
## ses -0.174 0.054 -3.232 0.001
## alien67 ~
## ses (sa) -0.614 0.056 -10.879 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .education 2.941 0.499 5.894 0.000
## .sei 260.713 18.212 14.315 0.000
## .anomia67 4.011 0.343 11.708 0.000
## .powerless67 3.188 0.271 11.763 0.000
## .anomia71 3.697 0.373 9.917 0.000
## .powerless71 3.621 0.292 12.417 0.000
## ses 6.659 0.640 10.404 0.000
## .alien67 5.301 0.472 11.235 0.000
## .alien71 3.737 0.387 9.658 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## IndirectEffect -0.433 0.048 -9.100 0.000You can tell that this model is substantially worse. But what if we just remove the powerlessness correlated residuals?
require(lavaan)
# the model
wheaton.model_only_corr_resid = '
# measurement model
ses =~ education + sei
alien67 =~ anomia67 + powerless67
alien71 =~ anomia71 + powerless71
# structural model
alien71 ~ aa*alien67 + ses
alien67 ~ sa*ses
# correlated residuals (anomia)
anomia67 ~~ anomia71
# Indirect effect
IndirectEffect := sa*aa
'Fit the model.
Nearly identical AIC but some improvement in BIC.
require(lavaan)
alienation_only_corr_resid <- sem(wheaton.model_only_corr_resid, sample.cov=wheaton.cov, sample.nobs=932)
summary(alienation_only_corr_resid, fit.measures=T)## lavaan 0.6-11 ended normally after 78 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 932
##
## Model Test User Model:
##
## Test statistic 6.338
## Degrees of freedom 5
## P-value (Chi-square) 0.275
##
## Model Test Baseline Model:
##
## Test statistic 2133.722
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.999
## Tucker-Lewis Index (TLI) 0.998
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -15214.075
## Loglikelihood unrestricted model (H1) -15210.906
##
## Akaike (AIC) 30460.150
## Bayesian (BIC) 30537.548
## Sample-size adjusted Bayesian (BIC) 30486.733
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.017
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.051
## P-value RMSEA <= 0.05 0.944
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.011
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ses =~
## education 1.000
## sei 5.163 0.421 12.259 0.000
## alien67 =~
## anomia67 1.000
## powerless67 1.027 0.053 19.330 0.000
## alien71 =~
## anomia71 1.000
## powerless71 0.971 0.049 19.657 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## alien71 ~
## alien67 (aa) 0.617 0.050 12.427 0.000
## ses -0.211 0.049 -4.294 0.000
## alien67 ~
## ses (sa) -0.550 0.053 -10.298 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .anomia67 ~~
## .anomia71 1.885 0.240 7.868 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .education 2.727 0.515 5.290 0.000
## .sei 266.610 18.164 14.678 0.000
## .anomia67 5.060 0.370 13.658 0.000
## .powerless67 2.212 0.317 6.976 0.000
## .anomia71 4.807 0.395 12.181 0.000
## .powerless71 2.680 0.329 8.138 0.000
## ses 6.873 0.657 10.463 0.000
## .alien67 4.700 0.432 10.869 0.000
## .alien71 3.862 0.343 11.262 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## IndirectEffect -0.339 0.040 -8.564 0.000Session info.
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] qgraph_1.9.2 semPlot_1.1.5 lavaan_0.6-11
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-157 RColorBrewer_1.1-3 mi_1.1
## [4] tools_4.2.1 backports_1.4.1 bslib_0.3.1
## [7] utf8_1.2.2 R6_2.5.1 rpart_4.1.16
## [10] Hmisc_4.7-0 DBI_1.1.3 colorspace_2.0-3
## [13] nnet_7.3-17 tidyselect_1.1.2 gridExtra_2.3
## [16] mnormt_2.1.0 compiler_4.2.1 fdrtool_1.2.17
## [19] cli_3.3.0 htmlTable_2.4.1 sass_0.4.1
## [22] scales_1.2.0 checkmate_2.1.0 psych_2.2.5
## [25] pbapply_1.5-0 sem_3.1-15 stringr_1.4.0
## [28] digest_0.6.29 pbivnorm_0.6.0 foreign_0.8-82
## [31] minqa_1.2.4 rmarkdown_2.14 base64enc_0.1-3
## [34] jpeg_0.1-9 pkgconfig_2.0.3 htmltools_0.5.2
## [37] lme4_1.1-29 lisrelToR_0.1.5 highr_0.9
## [40] fastmap_1.1.0 htmlwidgets_1.5.4 rlang_1.0.3
## [43] rstudioapi_0.13 jquerylib_0.1.4 generics_0.1.2
## [46] jsonlite_1.8.0 gtools_3.9.3 dplyr_1.0.9
## [49] zip_2.2.0 magrittr_2.0.3 OpenMx_2.20.6
## [52] Formula_1.2-4 interp_1.1-2 Matrix_1.4-1
## [55] Rcpp_1.0.8.3 munsell_0.5.0 fansi_1.0.3
## [58] abind_1.4-5 rockchalk_1.8.152 lifecycle_1.0.1
## [61] stringi_1.7.6 yaml_2.3.5 carData_3.0-5
## [64] MASS_7.3-57 plyr_1.8.7 grid_4.2.1
## [67] parallel_4.2.1 crayon_1.5.1 deldir_1.0-6
## [70] lattice_0.20-45 kutils_1.70 splines_4.2.1
## [73] knitr_1.39 pillar_1.7.0 igraph_1.3.4
## [76] boot_1.3-28 corpcor_1.6.10 reshape2_1.4.4
## [79] stats4_4.2.1 XML_3.99-0.10 glue_1.6.2
## [82] evaluate_0.15 latticeExtra_0.6-30 RcppParallel_5.1.5
## [85] data.table_1.14.2 png_0.1-7 vctrs_0.4.1
## [88] nloptr_2.0.3 gtable_0.3.0 purrr_0.3.4
## [91] assertthat_0.2.1 ggplot2_3.3.6 xfun_0.31
## [94] openxlsx_4.2.5 xtable_1.8-4 coda_0.19-4
## [97] glasso_1.11 survival_3.3-1 tibble_3.1.7
## [100] arm_1.12-2 cluster_2.1.3 ellipsis_0.3.2