SEM madness with the Wheaton data
Dr. Thomas Pollet, Northumbria University (thomas.pollet@northumbria.ac.uk)
2020-10-08 | disclaimer
Hypothesised model
Data
require(lavaan)## Loading required package: lavaan## This is lavaan 0.6-7## lavaan is BETA 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-7 ended normally after 84 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free 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-7 ended normally after 82 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free 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-7 ended normally after 72 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free 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-7 ended normally after 73 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free 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-7 ended normally after 78 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free 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.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] qgraph_1.6.5 semPlot_1.1.2 lavaan_0.6-7
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-149 RColorBrewer_1.1-2 mi_1.0
## [4] tools_4.0.2 backports_1.1.10 R6_2.4.1
## [7] d3Network_0.5.2.1 rpart_4.1-15 Hmisc_4.4-1
## [10] colorspace_1.4-1 nnet_7.3-14 tidyselect_1.1.0
## [13] gridExtra_2.3 mnormt_2.0.2 compiler_4.0.2
## [16] fdrtool_1.2.15 htmlTable_2.1.0 regsem_1.6.2
## [19] scales_1.1.1 checkmate_2.0.0 psych_2.0.9
## [22] pbapply_1.4-3 sem_3.1-11 stringr_1.4.0
## [25] digest_0.6.25 pbivnorm_0.6.0 foreign_0.8-80
## [28] minqa_1.2.4 rmarkdown_2.4 base64enc_0.1-3
## [31] jpeg_0.1-8.1 pkgconfig_2.0.3 htmltools_0.5.0
## [34] lme4_1.1-23 lisrelToR_0.1.4 htmlwidgets_1.5.2
## [37] rlang_0.4.7 huge_1.3.4.1 rstudioapi_0.11
## [40] generics_0.0.2 gtools_3.8.2 dplyr_1.0.2
## [43] zip_2.1.1 magrittr_1.5 OpenMx_2.18.1
## [46] Formula_1.2-3 Matrix_1.2-18 Rcpp_1.0.5
## [49] munsell_0.5.0 abind_1.4-5 rockchalk_1.8.144
## [52] lifecycle_0.2.0 whisker_0.4 stringi_1.5.3
## [55] yaml_2.2.1 carData_3.0-4 MASS_7.3-53
## [58] plyr_1.8.6 matrixcalc_1.0-3 grid_4.0.2
## [61] parallel_4.0.2 crayon_1.3.4 lattice_0.20-41
## [64] kutils_1.70 splines_4.0.2 tmvnsim_1.0-2
## [67] knitr_1.30 pillar_1.4.6 igraph_1.2.5
## [70] rjson_0.2.20 boot_1.3-25 corpcor_1.6.9
## [73] BDgraph_2.63 reshape2_1.4.4 stats4_4.0.2
## [76] XML_3.99-0.5 glue_1.4.2 evaluate_0.14
## [79] latticeExtra_0.6-29 data.table_1.13.0 png_0.1-7
## [82] vctrs_0.3.4 nloptr_1.2.2.2 gtable_0.3.0
## [85] purrr_0.3.4 ggplot2_3.3.2 xfun_0.18
## [88] openxlsx_4.2.2 xtable_1.8-4 coda_0.19-4
## [91] Rsolnp_1.16 glasso_1.11 survival_3.2-7
## [94] truncnorm_1.0-8 tibble_3.0.3 arm_1.11-2
## [97] cluster_2.1.0 statmod_1.4.34 ellipsis_0.3.1