library(psych)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%() masks psych::%+%()
## ✖ ggplot2::alpha() masks psych::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggplot2)
library(lavaan)## This is lavaan 0.6-20
## lavaan is FREE software! Please report any bugs.
##
## Attaching package: 'lavaan'
##
## The following object is masked from 'package:psych':
##
## cor2covlibrary(stats)
library(lavaanPlot)## Warning: package 'lavaanPlot' was built under R version 4.5.2region_revd_educ <- read.csv("251120_region_educ.csv", header = TRUE)table(region_revd_educ$region)##
## 1 2 3 4 5 6
## 4260 803 1503 788 1789 1166Monrovia <- region_revd_educ %>%
filter(region == "1")`
mean(Monrovia$wi_continuous)## [1] 159328.6Northwest<- region_revd_educ %>%
filter(region == "2")mean(Northwest$wi_continuous)## [1] -16697.58So_Central <- region_revd_educ %>%
filter(region == "3")mean(So_Central$wi_continuous)## [1] 90770.46So_Easta <- region_revd_educ %>%
filter(region == "4")mean(So_Easta$wi_continuous)## [1] 1185.193So_Eastb <- region_revd_educ %>%
filter(region == "5")mean(So_Eastb$wi_continuous)## [1] 17560.73No_Central <- region_revd_educ %>%
filter(region == "6")mean(No_Central$wi_continuous)## [1] 50960.07### Run the full modelpov.model <- '
health =~ h2o_source + share_toilet
educ =~ p1_educ + p2_educ
wealth =~ radio + watch + computer + car_truck + sewing_machine + cell_phone
wealth ~~ health + educ
wealth ~ health
wealth ~ educ'fit <- cfa(pov.model, std.lv = TRUE, ordered = c("share_toilet", "radio", "car_truck", "cell_phone", "watch", "computer", "sewing_machine"), data = region_revd_educ)## Warning: lavaan->lav_model_vcov():
## Could not compute standard errors! The information matrix could not be
## inverted. This may be a symptom that the model is not identified.## Warning: lavaan->lav_test_satorra_bentler():
## could not invert information matrix needed for robust test statistic## Warning: lavaan->lav_object_post_check():
## some estimated ov variances are negativesummary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-20 ended normally after 46 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 28
##
## Number of observations 10309
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 179.919 NA
## Degrees of freedom 30 30
## P-value (Chi-square) 0.000 NA
## Scaling correction factor NA
## Shift parameter NA
##
##
## Model Test Baseline Model:
##
## Test statistic 18471.120 15649.402
## Degrees of freedom 45 45
## P-value 0.000 0.000
## Scaling correction factor 1.181
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.992 NA
## Tucker-Lewis Index (TLI) 0.988 NA
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.022 NA
## 90 Percent confidence interval - lower 0.019 NA
## 90 Percent confidence interval - upper 0.025 NA
## P-value H_0: RMSEA <= 0.050 1.000 NA
## P-value H_0: RMSEA >= 0.080 0.000 NA
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.043 0.043
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health =~
## h2o_source 1.307 NA 1.307 0.106
## share_toilet -1.072 NA -1.072 -1.072
## educ =~
## p1_educ 1.665 NA 1.665 1.078
## p2_educ -0.846 NA -0.846 -0.417
## wealth =~
## radio 0.639 NA 0.726 0.726
## watch 0.650 NA 0.739 0.739
## computer 0.751 NA 0.853 0.853
## car_truck 0.729 NA 0.828 0.828
## sewing_machine 0.447 NA 0.508 0.508
## cell_phone 0.745 NA 0.847 0.847
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wealth ~
## health 0.271 NA 0.238 0.238
## educ 0.129 NA 0.114 0.114
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health ~~
## .wealth 0.290 NA 0.290 0.290
## educ ~~
## .wealth 0.130 NA 0.130 0.130
## health ~~
## educ 0.154 NA 0.154 0.154
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source 75.135 NA 75.135 6.095
## .p1_educ 1.708 NA 1.708 1.106
## .p2_educ 2.731 NA 2.731 1.346
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## share_toilt|t1 -0.430 NA -0.430 -0.430
## radio|t1 -0.640 NA -0.640 -0.640
## watch|t1 -0.324 NA -0.324 -0.324
## computer|t1 1.953 NA 1.953 1.953
## car_truck|t1 1.640 NA 1.640 1.640
## sewing_mchn|t1 1.643 NA 1.643 1.643
## cell_phone|t1 -0.135 NA -0.135 -0.135
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source 150.247 NA 150.247 0.989
## .share_toilet -0.150 -0.150 -0.150
## .p1_educ -0.387 NA -0.387 -0.162
## .p2_educ 3.399 NA 3.399 0.826
## .radio 0.473 0.473 0.473
## .watch 0.454 0.454 0.454
## .computer 0.272 0.272 0.272
## .car_truck 0.314 0.314 0.314
## .sewing_machine 0.742 0.742 0.742
## .cell_phone 0.283 0.283 0.283
## health 1.000 1.000 1.000
## educ 1.000 1.000 1.000
## .wealth 1.000 0.774 0.774### This code did not run. Uncomment it to see the errors.
#res <- resid(fit, type = "standardized")
#reslavaanPlot(model = fit, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"), coefs = TRUE)### Run the model just on So_Easta. Computer removed b/c it had only 1 level.pov.model.sea <- '
health =~ h2o_source + share_toilet
educ =~ p1_educ + p2_educ
wealth =~ radio + watch + car_truck + sewing_machine + cell_phone
wealth ~~ health + educ
wealth ~ health
wealth ~ educ'fit <- cfa(pov.model.sea, std.lv = TRUE, ordered = c("share_toilet", "radio", "car_truck", "cell_phone", "watch", "sewing_machine"), data = So_Easta)## Warning: lavaan->lav_model_vcov():
## Could not compute standard errors! The information matrix could not be
## inverted. This may be a symptom that the model is not identified.## Warning: lavaan->lav_test_satorra_bentler():
## could not invert information matrix needed for robust test statistic## Warning: lavaan->lav_object_post_check():
## some estimated ov variances are negativesummary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-20 ended normally after 2559 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 26
##
## Number of observations 788
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 153.592 NA
## Degrees of freedom 22 22
## P-value (Chi-square) 0.000 NA
## Scaling correction factor NA
## Shift parameter NA
##
##
## Model Test Baseline Model:
##
## Test statistic 4065.457 2480.989
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.648
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.967 NA
## Tucker-Lewis Index (TLI) 0.947 NA
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.087 NA
## 90 Percent confidence interval - lower 0.074 NA
## 90 Percent confidence interval - upper 0.100 NA
## P-value H_0: RMSEA <= 0.050 0.000 NA
## P-value H_0: RMSEA >= 0.080 0.828 NA
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.143 0.143
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health =~
## h2o_source 129.594 NA 129.594 17.088
## share_toilet 0.025 NA 0.025 0.025
## educ =~
## p1_educ 1.330 NA 1.330 0.988
## p2_educ -1.185 NA -1.185 -0.549
## wealth =~
## radio 10.030 NA 0.820 0.820
## watch 7.688 NA 0.628 0.628
## car_truck 6.260 NA 0.511 0.511
## sewing_machine 5.212 NA 0.426 0.426
## cell_phone 8.971 NA 0.733 0.733
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wealth ~
## health -0.918 NA -11.231 -11.231
## educ -0.371 NA -4.540 -4.540
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health ~~
## .wealth 0.919 NA 0.919 0.919
## educ ~~
## .wealth 0.385 NA 0.385 0.385
## health ~~
## educ -0.001 NA -0.001 -0.001
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source 75.100 NA 75.100 9.902
## .p1_educ 1.410 NA 1.410 1.047
## .p2_educ 2.841 NA 2.841 1.317
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## share_toilt|t1 -0.900 NA -0.900 -0.900
## radio|t1 -0.511 NA -0.511 -0.511
## watch|t1 0.070 NA 0.070 0.070
## car_truck|t1 2.276 NA 2.276 2.276
## sewing_mchn|t1 1.998 NA 1.998 1.998
## cell_phone|t1 0.479 NA 0.479 0.479
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source -16737.113 NA -16737.113 -290.983
## .share_toilet 0.999 0.999 0.999
## .p1_educ 0.044 NA 0.044 0.024
## .p2_educ 3.248 NA 3.248 0.698
## .radio 0.328 0.328 0.328
## .watch 0.605 0.605 0.605
## .car_truck 0.738 0.738 0.738
## .sewing_machine 0.819 0.819 0.819
## .cell_phone 0.463 0.463 0.463
## health 1.000 1.000 1.000
## educ 1.000 1.000 1.000
## .wealth 1.000 149.787 149.787### Run the model for So_Eastb. Computer removed b/c it had only 1 level.pov.model.seb <- '
health =~ h2o_source + share_toilet
educ =~ p1_educ + p2_educ
wealth =~ radio + watch + car_truck + sewing_machine + cell_phone
wealth ~~ health + educ
wealth ~ health
wealth ~ educ'fit <- cfa(pov.model.seb, std.lv = TRUE, ordered = c("share_toilet", "radio", "car_truck", "cell_phone", "watch", "sewing_machine"), data = So_Eastb)## Warning: lavaan->lav_samplestats_step2():
## correlation between variables car_truck and share_toilet is (nearly) 1.0## Warning: lavaan->lav_model_vcov():
## Could not compute standard errors! The information matrix could not be
## inverted. This may be a symptom that the model is not identified.## Warning: lavaan->lav_test_satorra_bentler():
## could not invert information matrix needed for robust test statistic## Warning: lavaan->lav_object_post_check():
## some estimated ov variances are negative## Warning: lavaan->lav_object_post_check():
## covariance matrix of latent variables is not positive definite ; use
## lavInspect(fit, "cov.lv") to investigate.summary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-20 ended normally after 170 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 26
##
## Number of observations 1789
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 140.424 NA
## Degrees of freedom 22 22
## P-value (Chi-square) 0.000 NA
## Scaling correction factor NA
## Shift parameter NA
##
##
## Model Test Baseline Model:
##
## Test statistic 2962.216 2195.213
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.355
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.960 NA
## Tucker-Lewis Index (TLI) 0.934 NA
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.055 NA
## 90 Percent confidence interval - lower 0.046 NA
## 90 Percent confidence interval - upper 0.064 NA
## P-value H_0: RMSEA <= 0.050 0.167 NA
## P-value H_0: RMSEA >= 0.080 0.000 NA
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.140 0.140
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health =~
## h2o_source 0.232 NA 0.232 0.049
## share_toilet -0.130 NA -0.130 -0.130
## educ =~
## p1_educ 1.848 NA 1.848 1.352
## p2_educ -0.846 NA -0.846 -0.389
## wealth =~
## radio 0.608 NA 0.809 0.809
## watch 0.563 NA 0.749 0.749
## car_truck 0.558 NA 0.743 0.743
## sewing_machine 0.525 NA 0.699 0.699
## cell_phone 0.639 NA 0.850 0.850
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## wealth ~
## health 0.084 NA 0.063 0.063
## educ 0.147 NA 0.110 0.110
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## health ~~
## .wealth 4.296 NA 4.296 4.296
## educ ~~
## .wealth 0.017 NA 0.017 0.017
## health ~~
## educ 0.729 NA 0.729 0.729
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source 74.670 NA 74.670 15.707
## .p1_educ 1.380 NA 1.380 1.010
## .p2_educ 2.874 NA 2.874 1.321
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## share_toilt|t1 -0.849 NA -0.849 -0.849
## radio|t1 -0.370 NA -0.370 -0.370
## watch|t1 -0.098 NA -0.098 -0.098
## car_truck|t1 2.029 NA 2.029 2.029
## sewing_mchn|t1 1.653 NA 1.653 1.653
## cell_phone|t1 0.638 NA 0.638 0.638
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .h2o_source 22.546 NA 22.546 0.998
## .share_toilet 0.983 0.983 0.983
## .p1_educ -1.547 NA -1.547 -0.829
## .p2_educ 4.015 NA 4.015 0.849
## .radio 0.345 0.345 0.345
## .watch 0.438 0.438 0.438
## .car_truck 0.448 0.448 0.448
## .sewing_machine 0.512 0.512 0.512
## .cell_phone 0.277 0.277 0.277
## health 1.000 1.000 1.000
## educ 1.000 1.000 1.000
## .wealth 1.000 0.564 0.564