SEM madness with the Wheaton data
Dr. Thomas Pollet, Northumbria University (thomas.pollet@northumbria.ac.uk)
2026-03-24 | disclaimer
Hypothesised model
Data
require(lavaan)## Loading required package: lavaan## This is lavaan 0.6-21
## 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"))## Warning: 'getCov' is deprecated.
## Use 'lav_getcov' instead.
## See help("Deprecated")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-21 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 (SABIC) 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 H_0: RMSEA <= 0.050 0.930
## P-value H_0: RMSEA >= 0.080 0.001
##
## 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-21 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 (SABIC) 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 H_0: RMSEA <= 0.050 0.930
## P-value H_0: RMSEA >= 0.080 0.001
##
## 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-21 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 (SABIC) 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 H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.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-21 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 (SABIC) 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 H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 0.983
##
## 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-21 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 (SABIC) 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 H_0: RMSEA <= 0.050 0.944
## P-value H_0: RMSEA >= 0.080 0.000
##
## 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.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.3
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/London
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] qgraph_1.9.8 semPlot_1.1.6 lavaan_0.6-21
##
## loaded via a namespace (and not attached):
## [1] psych_2.6.1 tidyselect_1.2.1 dplyr_1.2.0
## [4] farver_2.1.2 S7_0.2.1 fastmap_1.2.0
## [7] OpenMx_2.22.7 XML_3.99-0.18 digest_0.6.39
## [10] rpart_4.1.24 mi_1.1 lifecycle_1.0.5
## [13] cluster_2.1.8.1 magrittr_2.0.4 compiler_4.5.2
## [16] rlang_1.1.7 Hmisc_5.2-4 sass_0.4.10
## [19] tools_4.5.2 igraph_2.1.4 yaml_2.3.12
## [22] data.table_1.18.2.1 knitr_1.51 htmlwidgets_1.6.4
## [25] mnormt_2.1.1 plyr_1.8.9 RColorBrewer_1.1-3
## [28] abind_1.4-8 foreign_0.8-90 nnet_7.3-20
## [31] grid_4.5.2 stats4_4.5.2 colorspace_2.1-2
## [34] xtable_1.8-4 ggplot2_4.0.2 gtools_3.9.5
## [37] scales_1.4.0 MASS_7.3-65 cli_3.6.5
## [40] rmarkdown_2.30 reformulas_0.4.4 generics_0.1.4
## [43] RcppParallel_5.1.10 otel_0.2.0 rstudioapi_0.17.1
## [46] reshape2_1.4.5 pbapply_1.7-4 minqa_1.2.8
## [49] cachem_1.1.0 stringr_1.6.0 splines_4.5.2
## [52] parallel_4.5.2 base64enc_0.1-6 vctrs_0.7.1
## [55] boot_1.3-32 Matrix_1.7-4 jsonlite_2.0.0
## [58] carData_3.0-6 glasso_1.11 Formula_1.2-5
## [61] htmlTable_2.4.3 jpeg_0.1-11 jquerylib_0.1.4
## [64] glue_1.8.0 nloptr_2.2.1 stringi_1.8.7
## [67] sem_3.1-16 gtable_0.3.6 quadprog_1.5-8
## [70] lme4_1.1-38 tibble_3.3.1 lisrelToR_0.3
## [73] pillar_1.11.1 htmltools_0.5.9 R6_2.6.1
## [76] Rdpack_2.6.6 evaluate_1.0.5 pbivnorm_0.6.0
## [79] lattice_0.22-7 rbibutils_2.4.1 png_0.1-8
## [82] backports_1.5.0 rockchalk_1.8.157 kutils_1.73
## [85] openxlsx_4.2.8 arm_1.14-4 corpcor_1.6.10
## [88] bslib_0.10.0 fdrtool_1.2.18 Rcpp_1.1.1
## [91] zip_2.3.3 coda_0.19-4.1 gridExtra_2.3
## [94] nlme_3.1-168 checkmate_2.3.3 xfun_0.56
## [97] pkgconfig_2.0.3