# principles and practice of sem, 5th ed.
# rex b. kline, guilford press, 2023
# chapter 22, table 22.1, analysis 2
# two-factor cfa model of divergent thinking
# analyzed simultaneously over chinese and american samples
# to avoid the problem that some R packages share the same name
# for different functions, all functions are specified next as
# package::function
# which prevents masking, or the default hiding of
# a function with a redundant name from a package used next
date()
v <- R.Version()
print(paste0(v$language, " version ", v$major, ".",
v$minor, " (", v$year, "-", v$month, "-", v$day, ")"))
library(lavaan)
# get citation information
citation("lavaan", auto = TRUE)
# name subtests
subtest.name <- c("LM1", "LM2", "LM3", "RW1", "RW2")
# chinese, n = 316, originality
chineseLower.cor <- '
1.000
.715 1.000
.631 .738 1.000
.377 .408 .319 1.000
.452 .495 .391 .445 1.000 '
# name the variables and convert to full correlation matrix
chinese.cor <- lavaan::getCov(chineseLower.cor, names = subtest.name)
# add the standard deviations and convert to covariances
chinese.sd <- c(1.532, 1.878, 2.012, 1.002, 1.055)
chinese.cov <- lavaan::cor2cov(chinese.cor, sds = chinese.sd)
# create mean vector
chinese.mean = c(1.24, 1.47, 1.86, .77, .89)
# american, n = 302, originality
americanLower.cor <- '
1.000
.519 1.000
.543 .639 1.000
.280 .321 .317 1.000
.242 .380 .300 .665 1.000 '
# name the variables and convert to full correlation matrix
american.cor <- lavaan::getCov(americanLower.cor, names = subtest.name)
# add the standard deviations and convert to covariances
american.sd <- c(1.094, 1.350, 1.494, 2.257, 1.541)
american.cov <- lavaan::cor2cov(american.cor, sds = american.sd)
# create mean vector
american.mean = c(.88, 1.35, 1.27, 1.71, 1.37)
# combine covariances, group sizes, and means
# into single list objects
combined.cov <- list(chinese = chinese.cov, american = american.cov)
combined.n <- list(chinese = 316, american = 302)
combined.mean <- list(chinese = chinese.mean, american = american.mean)
# specify configural invariance model
configural.model <- '
# turn off (NA) default reference variable scaling for
# first indicator of each factor
# label loadings and intercepts for both factors
# over both groups
LineMeaning =~ NA*LM1 + c(a1,a2)*LM1 + c(b1,b2)*LM2 + c(c1,c2)*LM3
LM1 ~ c(d1,d2)*1
LM2 ~ c(e1,e2)*1
LM3 ~ c(f1,f2)*1
RealWorld =~ NA*RW1 + c(g1,g2)*RW1 + c(h1,h2)*RW2
RW1 ~ c(i1,i2)*1
RW2 ~ c(j1,j2)*1
# free factor means in both groups
LineMeaning ~ c(m1,m2)*1
RealWorld ~ c(q1,q2)*1 '
# specify configural invariance model constraints
# all constraints are for effects coding only
# average loading is 1, average intercept is zero
# there are no cross-group equality constraints
configural.constraints <- '
# LineMeaning, chinese, loadings and intercepts
3 - a1 - b1 - c1 == 0
0 - d1 - e1 - f1 == 0
# LineMeaning, american, loadings and intercepts
3 - a2 - b2 - c2 == 0
0 - d2 - e2 - f2 == 0
# RealWorld, chinese, loadings and intercepts
2 - g1 - h1 == 0
0 - i1 - j1 == 0
# Real World, american, loadings and intercepts
2 - g2 - h2 == 0
0 - i2 - j2 == 0 '
# specify weak invariance model constraints
# when loadings are constrained to equality over groups,
# release the effects coding scaling constraint that the
# average loading is 1.0 in one group (here, american)
weak.constraints <- '
# LineMeaning, chinese, loadings and intercepts
3 - a1 - b1 - c1 == 0
0 - d1 - e1 - f1 == 0
# LineMeaning, american, loadings and intercepts
# effects coding scaling constraint for
# loadings is released in this group by "commenting
# out" the original specification
# 3 - a2 - b2 - c2 == 0
0 - d2 - e2 - f2 == 0
# RealWorld, chinese, loadings and intercepts
2 - g1 - h1 == 0
0 - i1 - j1 == 0
# RealWorld, american, loadings and intercepts
# release constraint for 2nd factor in this group
# 2 - g2 - h2 == 0
0 - i2 - j2 == 0
# add equality constraints for loadings
a1 == a2
b1 == b2
c1 == c2
g1 == g2
h1 == h2 '
# specify constraints for strong invariance model
# when intercepts are also constrained to equality over groups,
# release the effects coding scaling constraint that the
# average loading is 0 in one group (here, american)
strong.constraints <- '
# LineMeaning, chinese, loadings and intercepts
3 - a1 - b1 - c1 == 0
0 - d1 - e1 - f1 == 0
# LineMeaning, american, loadings and intercepts
# both effects coding constraints are released in this group
# 3 - a2 - b2 - c2 == 0
# 0 - d2 - e2 - f2 == 0
# RealWorld, chinese, loadings and intercepts
2 - g1 - h1 == 0
0 - i1 - j1 == 0
# RealWorld, american, loadings and intercepts
# both effects coding constraints are released in this group
# 2 - g2 - h2 == 0
# 0 - i2 - j2 == 0
# equality constraints for loadings
a1 == a2
b1 == b2
c1 == c2
g1 == g2
h1 == h2
# add equality constraints for intecepts
d1 == d2
e1 == e2
f1 == f2
i1 == i2
j1 == j2 '
# specify constraints for partial strong invariance model
# release equality constraints for two indicators
# with the highest modification indices
# LM2 and RW1
partialStrong.constraints <- '
# LineMeaning, chinese, loadings and intercepts
3 - a1 - b1 - c1 == 0
0 - d1 - e1 - f1 == 0
# LineMeaning, american, loadings and intercepts
# both effects coding constraints are released in this group
# 3 - a2 - b2 - c2 == 0
# 0 - d2 - e2 - f2 == 0
# RealWorld, chinese, loadings and intercepts
2 - g1 - h1 == 0
0 - i1 - j1 == 0
# # RealWorld, american, loadings and intercepts
# both effects coding constraints are released in this group
# 2 - g2 - h2 == 0
# 0 - i2 - j2 == 0
# equality constraints for loadings
a1 == a2
b1 == b2
c1 == c2
g1 == g2
h1 == h2
# add equality constraints for intecepts
d1 == d2
# release equality constraint on LM2 intercept
# e1 == e2
f1 == f2
# release equality constraint on RW1 intercept
# i1 == i2
j1 == j2 '
# fit models to data
configural <- lavaan::cfa(configural.model, sample.cov = combined.cov,
sample.mean = combined.mean, sample.nobs = combined.n,
meanstructure = TRUE, constraints = configural.constraints)
weak <- lavaan::cfa(configural.model, sample.cov = combined.cov,
sample.mean = combined.mean, sample.nobs = combined.n,
meanstructure = TRUE, constraints = weak.constraints)
strong <- lavaan::cfa(configural.model, sample.cov = combined.cov,
sample.mean = combined.mean, sample.nobs = combined.n,
meanstructure = TRUE, constraints = strong.constraints)
partialStrong <- lavaan::cfa(configural.model, sample.cov = combined.cov,
sample.mean = combined.mean, sample.nobs = combined.n,
meanstructure = TRUE, constraints = partialStrong.constraints)
# generate chi-square difference tests
anova(configural, weak)
anova(weak, strong)
anova(weak, partialStrong)
# configural model estimates, predicted covariances
# and means, and residuals
lavaan::summary(configural, standardized = TRUE, fit.measures = TRUE)
lavaan::fitted(configural)
lavaan::residuals(configural, type = "raw")
lavaan::residuals(configural, type = "standardized")
lavaan::residuals(configural, type = "cor.bollen")
# weak model estimates, predicted covariances
# and means, and residuals
lavaan::summary(weak, standardized = TRUE, fit.measures = TRUE)
lavaan::fitted(weak)
lavaan::residuals(weak, type = "raw")
lavaan::residuals(weak, type = "standardized")
lavaan::residuals(weak, type = "cor.bollen")
# strong model estimates, predicted covariances
# and means, residuals, and modification indices
# for equality-constrained parameters
lavaan::summary(strong, standardized = TRUE, fit.measures = TRUE)
lavaan::fitted (strong)
# observed means and standard deviations
chinese.mean
chinese.sd
american.mean
american.sd
lavaan::residuals(strong, type = "raw")
lavaan::residuals(strong, type = "standardized")
lavaan::residuals(strong, type = "cor.bollen")
# equality constraints for intercepts
# are labeled d, e, f, i, j
lavaan::lavTestScore(weak, univariate = TRUE)
# partial strong model estimates, predicted covariances
# and means, residuals, and modification indices
# for equality-constrained parameters
lavaan::summary(partialStrong, standardized = TRUE, fit.measures = TRUE)
# observed means, chinese, american
chinese.mean
american.mean
# predicted covariances and means
lavaan::fitted (partialStrong)
lavaan::residuals(partialStrong, type = "raw")
lavaan::residuals(partialStrong, type = "standardized")
lavaan::residuals(partialStrong, type = "cor.bollen")
lavaan::lavTestScore(partialStrong, univariate = TRUE)