library("lavaan.survey") ###################### ## Application 1 ## ###################### data("pisa.be.2003") model.pisa <- " math =~ PV1MATH1 + PV1MATH2 + PV1MATH3 + PV1MATH4 neg.efficacy =~ ST31Q01 + ST31Q02 + ST31Q03 + ST31Q04 + ST31Q05 + ST31Q06 + ST31Q07 + ST31Q08 neg.selfconcept =~ ST32Q02 + ST32Q04 + ST32Q06 + ST32Q07 + ST32Q09 neg.selfconcept ~ neg.efficacy + ESCS + male neg.efficacy ~ neg.selfconcept + school.type + ESCS + male math ~ neg.selfconcept + neg.efficacy + school.type + ESCS + male " fit <- lavaan(model.pisa, data = pisa.be.2003, auto.var = TRUE, std.lv = TRUE, meanstructure = TRUE, int.ov.free = TRUE, estimator = "MLM") des.rep <- svrepdesign(ids = ~1, weights = ~W_FSTUWT, data = pisa.be.2003, repweights = "W_FSTR[0-9]+", type = "Fay", rho = 0.5) fit.surv <- lavaan.survey(lavaan.fit = fit, survey.design = des.rep) fit # Show fitmeasures results fit.surv # Show fitmeasures results library("plyr") # For arrange() head(arrange(modificationIndices(fit.surv)[, -c(7, 9)], mi.scaled, decreasing = TRUE)) ####### Table 2 ####### which.coefs <- grep("[a-zA-Z]~[a-zA-Z]", names(coef(fit))) which.coefs <- which.coefs[-grep("(ESCS|male)", names(coef(fit))[which.coefs])] est.fit <- coef(fit)[which.coefs] se.fit <- sqrt(diag(vcov(fit)))[which.coefs] est.fit.surv <- coef(fit.surv)[which.coefs] se.fit.surv <- sqrt(diag(vcov(fit.surv)))[which.coefs] creff <- se.fit.surv/se.fit pars <- strsplit(names(est.fit), "~") restab.2 <- cbind('Naive'=est.fit, 'PML'=est.fit.surv, 'Naive'=se.fit,'PML'= se.fit.surv, creff) ## Best viewed with # options(digits=4) restab.2 ###################### ## Application 2 ## ###################### data("ess4.gb") model.cfa <- "range =~ gvjbevn + gvhlthc + gvslvol + gvslvue + gvcldcr + gvpdlwk goals =~ sbprvpv + sbeqsoc + sbcwkfm" fit.cfa.ml <- lavaan(model.cfa, data = ess4.gb, estimator = "MLM", meanstructure = TRUE, int.ov.free = TRUE, auto.var = TRUE, auto.fix.first = TRUE, auto.cov.lv.x = TRUE) fit.cfa.ml des.gb <- svydesign(ids = ~psu, strata = ~stratval, weights = ~dweight, data = ess4.gb) fit.cfa.surv <- lavaan.survey(fit.cfa.ml, survey.design = des.gb) fit.cfa.surv fit.cfa.surv.wls <- lavaan.survey(fit.cfa.ml, survey.design = des.gb, estimator = "WLS") fit.cfa.surv.wls.yb <- lavaan.survey(fit.cfa.ml, survey.design = des.gb, estimator = "WLS", estimator.gamma = "Yuan-Bentler") des.gb.rep <- as.svrepdesign(des.gb, type = "JKn") fit.cfa.surv.rep <- lavaan.survey(fit.cfa.ml, survey.design = des.gb.rep) fit.cfa.surv.rep ###### Table 3 ###### get.restab <- function(fit, other.fit.list, interest="~~") { which.coefs <- grep(interest, names(coef(fit))) est.fit <- coef(fit)[which.coefs] se.fit <- sqrt(diag(vcov(fit)))[which.coefs] other.fit.res <- lapply(other.fit.list, function(fit.surv) { list(est = coef(fit.surv)[which.coefs], se = sqrt(diag(vcov(fit.surv)))[which.coefs], creff = sqrt(diag(vcov(fit.surv)))[which.coefs]/se.fit ) }) pars <- strsplit(names(est.fit), "~~") restab <- rbind('Naive'=c(est.fit, se.fit, rep(NA, 3))) for(ifit in seq(along=other.fit.res)) { restab <- rbind(restab, with(other.fit.res[[ifit]], c(est, se, creff)) ) } rownames(restab) <- c("ML, Naive", names(other.fit.list)) restab <- restab[, c(1,4,7,2,5,8,3,6,9)] colnames(restab) <- rep(c("Est.","s.e.","creff"), 3) restab } restab.3 <- get.restab(fit.cfa.ml, list("ML, Taylor"=fit.cfa.surv, "ML, jackknife"=fit.cfa.surv.rep, "WLS"=fit.cfa.surv.wls, "WLS, Y-B"=fit.cfa.surv.wls.yb), interest="(range|goals)~~(range|goals)") ## Best viewed with # options(digits=4) restab.3 ###################### ## Application 3 ## ###################### data("liss") model.liss <- " cs08 =~ 1 * cs08a247 cs09 =~ 1 * cs09b247 cs10 =~ 1 * cs10c247 cs11 =~ 1 * cs11d247 cs09 ~ cs08 cs10 ~ cs09 cs11 ~ cs10 cs08a247 ~~ vare * cs08a247 cs09b247 ~~ vare * cs09b247 cs10c247 ~~ vare * cs10c247 cs11d247 ~~ vare * cs11d247 cs08 ~~ vart08 * cs08 reliab.ratio := vart08 / (vart08 + vare) " fit.liss <- lavaan(model.liss, auto.var = TRUE, meanstructure = TRUE, int.ov.free = TRUE, data = liss) des.liss <- svydesign(ids = ~nohouse_encr, prob = ~1, data = liss) fit.liss.surv <- lavaan.survey(fit.liss, des.liss) fit.liss.surv parameterEstimates(fit.liss.surv)[24, ] # Running the imputations takes a long time, load the imputed datasets instead. load("liss.imp.rda") ## The loaded mice object was generated using # set.seed(20140221) # liss.imp <- mice(liss, m = 100, method = "norm", maxit = 100) library("mitools") library("mice") # For complete() liss.implist <- lapply(seq(liss.imp$m), function(im) complete(liss.imp, im)) liss.implist <- imputationList(liss.implist) des.liss.imp <- svydesign(ids = ~nohouse_encr, prob = ~1, data = liss.implist) # Estimation can take a while (9.4 secs elapsed on my machine) fit.liss.surv.mi <- lavaan.survey(fit.liss, des.liss.imp) fit.liss.surv.mi parameterEstimates(fit.liss.surv.mi)[24, ] ###### Figure 3 ###### par(mar=c(2,4,2,0)) percent.missing <- sapply(liss[,-(1:2)], function(x) 100*mean(is.na((x)))) year <- 2008:2011 plot(percent.missing ~ year, type="n", axes=FALSE, ylim=c(0,100), xlab="", ylab="% item nonrespondents", main="Attrition in the LISS panel 2008-2011") pts <- c(25, 50, 75) abline(h=pts, col=c("#44bb44", "orange", "#bb4444"), lwd=0.5) lines(year, percent.missing, type="b", pch=20) axis(2, at=c(0,pts,100), labels=paste(c(0,pts,100), "%", sep=""), las=2) axis(1, at = year, labels=as.character(year)) ###################### ## Application 4 ## ###################### # Cardinale BJ, Bennett DM, Nelson CE, Gross K (2009). "Does Productivity Drive # Diversity or Vice Versa? A Test of the Multivariate Productivity-Diversity # Hypothesis in Streams." Ecology, 90(5), 1227–1241. # Model and data from this example were obtained from Jarrett Byrnes' GitHub: # https://github.com/jebyrnes/Ecological-SEMs-in-lavaan # Note that I changed the variable names to be easier to understand, e.g. # logNutrient was logN, PatchDiversity was SA, etc. data("cardinale") model.card <- ' PatchDiversity ~ logNutrient + logNutrient2 + StreamDiversity Biomass ~ PatchDiversity + logNutrient O2Production ~ logNutrient + Biomass logNutrient ~~ logNutrient2' fit.card <- sem(model.card, data = cardinale, fixed.x = FALSE, estimator = "MLM") des.card <- svydesign(ids = ~Stream, probs = ~1, data = cardinale) fit.card.survey <- lavaan.survey(fit.card, des.card, estimator = "MLM") pval.pFsum(fit.card.survey, survey.design = des.card)