if (FALSE) {
data(HolzingerSwineford1939, package = "lavaan")
HS.model <- ' visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9 '
## fit single-group model
fit <- bcfa(HS.model, data = HolzingerSwineford1939,
n.chains = 2, burnin = 1000, sample = 500)
## fit multigroup model
fitg <- bcfa(HS.model, data = HolzingerSwineford1939,
n.chains = 2, burnin = 1000, sample = 500, group = "school")
## Use fit.measures as a shortcut for global fitMeasures only
## - Note that indices calculated from the "df" are only appropriate under
## noninformative priors, such that pD approximates the number of estimated
## parameters counted under ML estimation; incremental fit indices
## introduce further complications)
AFIs <- ppmc(fit, thin = 10, fit.measures = c("srmr","chisq","rmsea","cfi"))
summary(AFIs) # summarize the whole vector in a data.frame
hist(AFIs, element = "rmsea") # only plot one discrepancy function at a time
plot(AFIs, element = "srmr")
## define a list of custom discrepancy functions
## - (global) fit measures
## - (local) standardized residuals
discFUN <- list(global = function(fit) {
fitMeasures(fit, fit.measures = c("cfi","rmsea","srmr","chisq"))
},
std.cov.resid = function(fit) lavResiduals(fit, zstat = FALSE,
summary = FALSE)$cov,
std.mean.resid = function(fit) lavResiduals(fit, zstat = FALSE,
summary = FALSE)$mean)
out1g <- ppmc(fit, discFUN = discFUN)
## summarize first discrepancy by default (fit indices)
summary(out1g)
## some model-implied correlations look systematically over/underestimated
summary(out1g, discFUN = "std.cov.resid", central.tendency = "EAP")
hist(out1g, discFUN = "std.cov.resid", element = c(1, 7))
plot(out1g, discFUN = "std.cov.resid", element = c("x1","x7"))
## For ease of investigation, optionally export summary as a data.frame,
## sorted by size of average residual
summary(out1g, discFUN = "std.cov.resid", central.tendency = "EAP",
to.data.frame = TRUE, sort.by = "EAP")
## or sorted by size of PPP
summary(out1g, discFUN = "std.cov.resid", central.tendency = "EAP",
to.data.frame = TRUE, sort.by = "PPP_sim_LessThan_obs")
## define a list of custom discrepancy functions for multiple groups
## (return each group's numeric output using a different function)
disc2g <- list(global = function(fit) {
fitMeasures(fit, fit.measures = c("cfi","rmsea","mfi","srmr","chisq"))
},
cor.resid1 = function(fit) lavResiduals(fit, zstat = FALSE,
type = "cor.bollen",
summary = FALSE)[[1]]$cov,
cor.resid2 = function(fit) lavResiduals(fit, zstat = FALSE,
type = "cor.bollen",
summary = FALSE)[[2]]$cov)
out2g <- ppmc(fitg, discFUN = disc2g, thin = 2)
## some residuals look like a bigger problem in one group than another
pairs(out2g, discFUN = "cor.resid1", horInd = 1:3, verInd = 7:9) # group 1
pairs(out2g, discFUN = "cor.resid2", horInd = 1:3, verInd = 7:9) # group 2
## print all to file: must be a LARGE picture. First group 1 ...
png("cor.resid1.png", width = 1600, height = 1200)
pairs(out2g, discFUN = "cor.resid1")
dev.off()
## ... then group 2
png("cor.resid2.png", width = 1600, height = 1200)
pairs(out2g, discFUN = "cor.resid2")
dev.off()
}
Run the code above in your browser using DataLab