Chapter 11 Week8_2: Lavaan Lab 8 Estimation Methods
In this lab, we will learn how to estimate parameters in CFA/SR models.
Load up the lavaan library:
library(lavaan)11.1 PART I: Hypothetical Example
11.1.1 One-factor CFA model
A made-up sample covariance matrix with n = 200:
n = 200
S_3fac = matrix(c(5, 2, 3.5, 2, 3, 2, 3.5, 2, 6), 3, 3,
dimnames = list(c('Y1', 'Y2', 'Y3'), c('Y1', 'Y2', 'Y3')))
S_3fac## Y1 Y2 Y3
## Y1 5.0 2 3.5
## Y2 2.0 3 2.0
## Y3 3.5 2 6.0Fit a one-factor CFA to the sample covariance matrix:
one_fac_syntax <- "
eta =~ Y1 + Y2 + Y3
"Request Unweighted Least Squares (ULS):
one_fac_fit2 <- lavaan::sem(one_fac_syntax,
sample.cov = S_3fac,
sample.nobs = n,
estimator = "ULS",
fixed.x = FALSE)
summary(one_fac_fit2, standardized = T)## lavaan 0.6-12 ended normally after 23 iterations
##
## Estimator ULS
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta =~
## Y1 1.000 1.871 0.837
## Y2 0.571 0.023 24.496 0.000 1.069 0.617
## Y3 1.000 0.050 19.950 0.000 1.871 0.764
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Y1 1.500 0.202 7.423 0.000 1.500 0.300
## .Y2 1.857 0.094 19.749 0.000 1.857 0.619
## .Y3 2.500 0.202 12.372 0.000 2.500 0.417
## eta 3.500 0.189 18.497 0.000 1.000 1.000Sigma:
fitted(one_fac_fit2)$cov## Y1 Y2 Y3
## Y1 5.0
## Y2 2.0 3.0
## Y3 3.5 2.0 6.0Sigma = fitted(one_fac_fit2)$cov
diff = Sigma[colnames(S_3fac), colnames(S_3fac)] - S_3fac
round(diff,3)## Y1 Y2 Y3
## Y1 0 0 0
## Y2 0 0 0
## Y3 0 0 0all zeros. Meaning that Sigma = S.
11.2 PART II: ULS on the Positive Affect Example
Let’s read this dataset in:
cfaData<- read.csv("cfaInclassData.csv", header = T)Fit a two-factor CFA model:
fixedIndTwoFacSyntax <- "
#Factor Specification
posAffect =~ glad + happy + cheerful
satisfaction =~ satisfied + content + comfortable
"two_fac_fit_uls <- lavaan::sem(fixedIndTwoFacSyntax,
data = cfaData,
fixed.x = FALSE,
estimator = "ULS")
summary(two_fac_fit_uls, standardized = T)## lavaan 0.6-12 ended normally after 33 iterations
##
## Estimator ULS
## Optimization method NLMINB
## Number of model parameters 13
##
## Number of observations 1000
##
## Model Test User Model:
##
## Test statistic 0.837
## Degrees of freedom 8
## P-value (Unknown) NA
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## posAffect =~
## glad 1.000 0.696 0.708
## happy 1.067 0.066 16.205 0.000 0.743 0.760
## cheerful 1.105 0.068 16.157 0.000 0.769 0.781
## satisfaction =~
## satisfied 1.000 0.775 0.768
## content 1.073 0.064 16.845 0.000 0.832 0.768
## comfortable 0.906 0.053 17.124 0.000 0.702 0.738
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## posAffect ~~
## satisfaction 0.264 0.017 15.795 0.000 0.489 0.489
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .glad 0.482 0.050 9.597 0.000 0.482 0.499
## .happy 0.402 0.054 7.488 0.000 0.402 0.422
## .cheerful 0.377 0.056 6.762 0.000 0.377 0.390
## .satisfied 0.417 0.055 7.615 0.000 0.417 0.410
## .content 0.482 0.060 8.106 0.000 0.482 0.411
## .comfortable 0.411 0.049 8.322 0.000 0.411 0.455
## posAffect 0.484 0.039 12.395 0.000 1.000 1.000
## satisfaction 0.600 0.045 13.455 0.000 1.000 1.000Sigma:
S = cov(cfaData[,-1])
Sigma = fitted(two_fac_fit_uls)$cov[colnames(S), colnames(S)]
diff = Sigma - S
round(diff,3)## glad cheerful happy satisfied content comfortable
## glad 0.000 -0.002 0.004 -0.014 0.002 0.009
## cheerful -0.002 0.000 -0.002 0.009 -0.004 0.003
## happy 0.004 -0.002 0.000 -0.003 -0.011 0.011
## satisfied -0.014 0.009 -0.003 0.000 0.009 -0.007
## content 0.002 -0.004 -0.011 0.009 0.000 -0.003
## comfortable 0.009 0.003 0.011 -0.007 -0.003 0.000print(paste0("The difference between S and Sigma ranged between ", round(min(diff),4), " and ", round(max(diff),4), "."))## [1] "The difference between S and Sigma ranged between -0.0143 and 0.0113."Sigma is not the same as S, but close.
The sum of squared differences is:
sum(diff[lower.tri(diff,diag = T)]^2)## [1] 0.000837587411.3 PART III: Calculate ULS test statistic manually
ULS test statistic is calculated as:
T_uls = (1000-1)*sum(diff[lower.tri(diff,diag = T)]^2)
T_uls## [1] 0.8367498One can also obtain vectors of S and Sigma first:
s = lav_matrix_vech(S)
sigma = lav_matrix_vech(Sigma)and calculate the ULS test statistic:
T_uls = (1000-1)*sum((s - sigma)^2)
T_uls## [1] 0.8367498- No p-value for ULS test statistic
- as there is no known distribution for this test statistic
- i.e., no suitable method for model fit evaluation
11.4 PART IV: ML vs ULS vs WLS
11.4.1 ML Estimation
two_fac_fit_ml <- lavaan::sem(fixedIndTwoFacSyntax,
data = cfaData,
fixed.x = FALSE,
estimator = "ML")11.4.2 WLS Estimation
two_fac_fit_wls <- lavaan::sem(fixedIndTwoFacSyntax,
data = cfaData,
fixed.x = FALSE,
estimator = "WLS")11.4.3 Compare the parameter estimates
coefTable = parameterEstimates(two_fac_fit_ml)[,1:3]
coefTable = cbind(coefTable,
ML = parameterEstimates(two_fac_fit_ml)$est,
ULS = parameterEstimates(two_fac_fit_uls)$est,
WLS = parameterEstimates(two_fac_fit_wls)$est)
coefTable## lhs op rhs ML ULS WLS
## 1 posAffect =~ glad 1.0000000 1.0000000 1.0000000
## 2 posAffect =~ happy 1.0665969 1.0674339 1.0645600
## 3 posAffect =~ cheerful 1.1121981 1.1045607 1.1150756
## 4 satisfaction =~ satisfied 1.0000000 1.0000000 1.0000000
## 5 satisfaction =~ content 1.0683189 1.0732316 1.0719435
## 6 satisfaction =~ comfortable 0.9175967 0.9059151 0.9151164
## 7 glad ~~ glad 0.4838829 0.4823845 0.4835721
## 8 happy ~~ happy 0.4048385 0.4021235 0.4052179
## 9 cheerful ~~ cheerful 0.3711009 0.3772155 0.3643077
## 10 satisfied ~~ satisfied 0.4185958 0.4165368 0.4182162
## 11 content ~~ content 0.4909589 0.4823048 0.4902739
## 12 comfortable ~~ comfortable 0.4003673 0.4114721 0.4034071
## 13 posAffect ~~ posAffect 0.4816239 0.4840887 0.4812181
## 14 satisfaction ~~ satisfaction 0.5973737 0.6004497 0.5964489
## 15 posAffect ~~ satisfaction 0.2617756 0.2638580 0.261915911.4.4 Compare the standard errors:
seTable = parameterEstimates(two_fac_fit_ml)[,1:3]
seTable = cbind(seTable,
ML = parameterEstimates(two_fac_fit_ml)$se,
ULS = parameterEstimates(two_fac_fit_uls)$se,
WLS = parameterEstimates(two_fac_fit_wls)$se)
seTable## lhs op rhs ML ULS WLS
## 1 posAffect =~ glad 0.00000000 0.00000000 0.00000000
## 2 posAffect =~ happy 0.05528238 0.06587199 0.05328971
## 3 posAffect =~ cheerful 0.05715937 0.06836322 0.05484498
## 4 satisfaction =~ satisfied 0.00000000 0.00000000 0.00000000
## 5 satisfaction =~ content 0.05204954 0.06371370 0.05010033
## 6 satisfaction =~ comfortable 0.04512163 0.05290389 0.04329804
## 7 glad ~~ glad 0.02906697 0.05026166 0.02803266
## 8 happy ~~ happy 0.02813623 0.05369964 0.02823333
## 9 cheerful ~~ cheerful 0.02853712 0.05578216 0.02731647
## 10 satisfied ~~ satisfied 0.02921942 0.05470277 0.02802346
## 11 content ~~ content 0.03376166 0.05950183 0.03304768
## 12 comfortable ~~ comfortable 0.02614297 0.04944382 0.02655468
## 13 posAffect ~~ posAffect 0.04210374 0.03905424 0.04141331
## 14 satisfaction ~~ satisfaction 0.04708792 0.04462502 0.04883263
## 15 posAffect ~~ satisfaction 0.02545399 0.01670551 0.0271884811.5 PART V: Improper Solutions
Going back to the 1-factor toy example…
Suppose we have a new covariance matrix now and a super large sample size (n = 200000000):
n = 200000000
S_3fac_new = matrix(c(5, 1, 3.5, 1, 3, 2, 3.5, 2, 6), 3, 3,
dimnames = list(c('Y1', 'Y2', 'Y3'), c('Y1', 'Y2', 'Y3')))
S_3fac_new## Y1 Y2 Y3
## Y1 5.0 1 3.5
## Y2 1.0 3 2.0
## Y3 3.5 2 6.0The one-factor syntax is the same:
one_fac_syntax <- "
eta =~ Y1 + Y2 + Y3
"ML Estimation:
one_fac_fit_new <- lavaan::sem(one_fac_syntax,
sample.cov = S_3fac_new,
sample.nobs = n,
estimator = "ML",
fixed.x = FALSE)## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negative- lavaan WARNING: some estimated ov variances are negative
- negative residual variances
- changing the estimation from ML to ULS doesn’t help
summary(one_fac_fit_new, standardized = T)## lavaan 0.6-12 ended normally after 28 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 2e+08
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta =~
## Y1 1.000 1.323 0.592
## Y2 0.571 0.000 6181.151 0.000 0.756 0.436
## Y3 2.000 0.000 4714.045 0.000 2.646 1.080
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Y1 3.250 0.000 7039.895 0.000 3.250 0.650
## .Y2 2.429 0.000 9150.327 0.000 2.429 0.810
## .Y3 -1.000 0.001 -760.286 0.000 -1.000 -0.167
## eta 1.750 0.001 3486.948 0.000 1.000 1.000Label and constraint sig3 to be larger than 0:
one_fac_syntax_const <- "
eta =~ Y1 + Y2 + Y3
Y3~~sig3*Y3
# constraints
sig3 > 0
"ML Estimation:
one_fac_fit_new2 <- lavaan::sem(one_fac_syntax_const,
sample.cov = S_3fac_new,
sample.nobs = n,
fixed.x = FALSE)## Warning in lav_object_post_check(object): lavaan WARNING: some estimated ov variances are negativesummary(one_fac_fit_new2, standardized = T)## lavaan 0.6-12 ended normally after 63 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
## Number of inequality constraints 1
##
## Number of observations 2e+08
##
## Model Test User Model:
##
## Test statistic 806452.704
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## eta =~
## Y1 1.000 1.429 0.639
## Y2 0.571 0.000 6357.073 0.000 0.816 0.471
## Y3 1.714 0.000 11748.540 0.000 2.449 1.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Y3 (sig3) -0.000 -0.000 -0.000
## .Y1 2.958 0.000 10000.000 0.000 2.958 0.592
## .Y2 2.333 0.000 10000.000 0.000 2.333 0.778
## eta 2.042 0.000 5065.023 0.000 1.000 1.000
##
## Constraints:
## |Slack|
## sig3 - 0 0.000fitted(one_fac_fit_new)$cov # perfect## Y1 Y2 Y3
## Y1 5.0
## Y2 1.0 3.0
## Y3 3.5 2.0 6.0fitted(one_fac_fit_new2)$cov # compromised, at the cost of model fit ## Y1 Y2 Y3
## Y1 5.000
## Y2 1.167 3.000
## Y3 3.500 2.000 6.000S_3fac_new## Y1 Y2 Y3
## Y1 5.0 1 3.5
## Y2 1.0 3 2.0
## Y3 3.5 2 6.0