#sem, openmx, lava, lavaan install.packages("lavaan") library(lavaan) ## browse help files for all functions in lavaan package help(package = lavaan) ex1 <- read.table("http://www.sunthud.com/lavaan1.csv", sep = ",", header = TRUE, na.strings = "-999999") ex2 <- read.table("http://www.sunthud.com/lavaan2.csv", sep = ",", header = TRUE, na.strings = "-999999") ex3 <- read.table("http://www.sunthud.com/lavaan3.csv", sep = ",", header = TRUE, na.strings = "-999999") ex4 <- read.table("http://www.sunthud.com/lavaan4.csv", sep = ",", header = TRUE, na.strings = "-999999") ex5 <- read.table("http://www.sunthud.com/lavaan5.csv", sep = ",", header = TRUE, na.strings = "-999999") ex6 <- read.table("http://www.sunthud.com/lavaan6.csv", sep = ",", header = TRUE, na.strings = "-999999") ################################################################################ ## example 1 -- Multiple Regression as Path Analysis ################################################################################ mod1 <- 'y ~ x1 + x2' # sem, cfa, growth fit1 <- sem(mod1, data = ex1, meanstructure = TRUE) summary(fit1, standardized=TRUE, fit = TRUE, rsquare = TRUE) ## compare to regression as a linear model m1 <- lm(y ~ x1 + x2, data = ex1) summary(m1) # free x fit1a <- sem(mod1, data = ex1, meanstructure = TRUE, fixed.x = FALSE) summary(fit1a, standardized=TRUE) ################################################################################ ## example 2 -- Path Analysis with Missing Data ################################################################################ ## take a look at how much missing data there is library(Amelia) missmap(ex2) Amelia::missmap(ex2) mod2 <- ' ## regressions y1 + y2 ~ x1 + x2 + x3 y3 ~ y1 + y2 + x2 ## correlated residual variances y1 ~~ y2 ' fit2 <- sem(mod2, data = ex2, meanstructure = TRUE, missing = "fiml") summary(fit2, fit.measures = TRUE, standardized = TRUE) ## Making an extra parameter mod2a <- ' ## regressions y1 ~ a1*x1 + a2*x2 + a3*x3 y2 ~ a4*x1 + a5*x2 + a6*x3 y3 ~ b1*y1 + b2*y2 + cp2*x2 ## correlated residual variances y1 ~~ y2 ## Define indirect effect ab11 := a1*b1 ab12 := a4*b2 ab21 := a2*b1 ab22 := a5*b2 ab31 := a3*b1 ab32 := a6*b2 ab1t := ab11 + ab12 ab2t := ab21 + ab22 ab3t := ab31 + ab32 ' fit2a <- sem(mod2a, data = ex2, meanstructure = TRUE, missing = "fiml") summary(fit2a, fit.measures = TRUE, standardized = TRUE) parameterEstimates(fit2a) fit2b <- sem(mod2a, data = ex2, meanstructure = TRUE, missing = "fiml", se = "boot", bootstrap = 200) parameterEstimates(fit2b, boot.ci.type = "bca.simple") library(semTools) fit2c <- sem.auxiliary(mod2a, data = ex2, meanstructure = TRUE, missing = "fiml", aux = c("z1", "z2")) summary(fit2c) ################################################################################ ## example 3 -- CFA using sufficient statistics ################################################################################ mymeans <- colMeans(ex3) mycov <- cov(ex3) nObs <- nrow(ex3) mod3 <- ' ## factor loadings Positive =~ great + cheerful + happy Negative =~ sad + down + unhappy ' # Marker variable method fit3 <- cfa(mod3, sample.cov = mycov, sample.mean = mymeans, sample.nobs = nObs, meanstructure = TRUE) summary(fit3, standardized = TRUE, fit.measures = TRUE) fitMeasures(fit3) coef(fit3) fitted.values(fit3) resid(fit3) vcov(fit3) standardizedSolution(fit3) modindices(fit3) # Fixed factor method fit3a <- cfa(mod3, sample.cov = mycov, sample.mean = mymeans, sample.nobs = nObs,meanstructure = TRUE, std.lv = TRUE) ################################################################################ ## example 4 -- Longitudinal CFA with measurement invariance constraints ################################################################################ mod4configural <- ' ## factor loadings Pos1 =~ great1 + cheerful1 + happy1 Pos2 =~ great2 + cheerful2 + happy2 ## correlated residual variances great1 ~~ great2 cheerful1 ~~ cheerful2 happy1 ~~ happy2 ' mod4weak <- ' ## factor loadings Pos1 =~ L1*great1 + L2*cheerful1 + L3*happy1 Pos2 =~ L1*great2 + L2*cheerful2 + L3*happy2 ## free factor variance at second time Pos2 ~~ NA*Pos2 ## correlated residual variances great1 ~~ great2 cheerful1 ~~ cheerful2 happy1 ~~ happy2 ' mod4strong <- ' ## factor loadings Pos1 =~ L1*great1 + L2*cheerful1 + L3*happy1 Pos2 =~ L1*great2 + L2*cheerful2 + L3*happy2 ## free factor variance at second time Pos2 ~~ NA*Pos2 ## correlated residual variances great1 ~~ great2 cheerful1 ~~ cheerful2 happy1 ~~ happy2 ## constrain intercept great1 ~ I1*1 cheerful1 ~ I2*1 happy1 ~ I3*1 great2 ~ I1*1 cheerful2 ~ I2*1 happy2 ~ I3*1 # free factor mean at second time Pos2 ~ NA*1 ' fit4c <- cfa(mod4configural, data = ex4, std.lv = TRUE, meanstructure = TRUE) fit4w <- cfa(mod4weak, data = ex4, std.lv = TRUE, meanstructure = TRUE) fit4s <- cfa(mod4strong, data = ex4, std.lv = TRUE, meanstructure = TRUE) summary(fit4c, standardized = TRUE, fit.measures = TRUE) summary(fit4w, standardized = TRUE, fit.measures = TRUE) summary(fit4s, standardized = TRUE, fit.measures = TRUE) anova(fit4c, fit4w) anova(fit4w, fit4s) predict(fit4c) mod4weaka <- ' ## factor loadings Pos1 =~ L1*great1 + L2*cheerful1 + L3*happy1 Pos2 =~ L1*great2 + L2*cheerful2 + L3*happy2 ## free factor variance at second time Pos2 ~~ NA*Pos2 ## correlated residual variances great1 ~~ great2 cheerful1 ~~ cheerful2 happy1 ~~ happy2 ## regression Pos2 ~ Pos1 ' fit4wa <- cfa(mod4weaka, data = ex4, std.lv = TRUE, meanstructure = TRUE) summary(fit4wa, standardized = TRUE, fit.measures = TRUE) # longInvariance mod4 <- ' Pos1 =~ great1 + cheerful1 + happy1 Pos2 =~ great2 + cheerful2 + happy2 ' # Create list of variables var1 <- c("great1", "cheerful1", "happy1") var2 <- c("great2", "cheerful2", "happy2") constrainedVar <- list(var1, var2) fit4 <- longInvariance(mod4, auto=1, constrainAuto=TRUE, varList=constrainedVar, data=ex4) fit4[[1]] fit4[[2]] fit4[[3]] fit4[[4]] ################################################################################ ## example 5 -- Multiple Group CFA ################################################################################ mod5configural <- ' Agency =~ Agency1 + Agency2 + Agency3 Intrinsic =~ Intrin1 + Intrin2 + Intrin3 Extrinsic =~ Extrin1 + Extrin2 + Extrin3 Positive =~ PosAFF1 + PosAFF2 + PosAFF3 ' fit5configural <- cfa(mod5configural, data = ex5, meanstructure = TRUE, group = "Sex") summary(fit5configural, standardized = TRUE, fit.measures = TRUE) ## instead of labeling parameters to constrain: c(L1, L1)*Agency1 ## lavaan can do this automatically with the argument "group.equal" ## weak invariance fit5weak <- sem(mod5configural, data = ex5, meanstructure = TRUE, group = "Sex", group.equal = "loadings") summary(fit5weak, standardized = TRUE, fit.measures = TRUE) ## strong invariance fit5strong <- sem(mod5configural, data = ex5, meanstructure = TRUE, group = "Sex", group.equal = c("loadings", "intercepts")) summary(fit5strong, standardized = TRUE, fit.measures = TRUE) # measurementInvariance function fit5 <- measurementInvariance(mod5configural, data = ex5, group = "Sex", std.lv = TRUE) summary(fit5[[2]]) ################################################################################ ## example 6 -- Latent Growth Curve model ################################################################################ mod6 <- ' ## initial status and shape of change factors intercept =~ 1*NegT1 + 1*NegT2 + 1*NegT3 + 1*NegT4 slope =~ 0*NegT1 + 1*NegT2 + 2*NegT3 + 3*NegT4 ## constraint residual variance to equality over time NegT1 ~~ th1*NegT1 NegT2 ~~ th1*NegT2 NegT3 ~~ th1*NegT3 NegT4 ~~ th1*NegT4 ' fit6 <- growth(mod6, data = ex6) summary(fit6, standardized = TRUE, fit.measures = TRUE) ## compare to multilevel regression model (HLM) ex6$id <- 1:nrow(ex6) ex6long <- rbind(data.frame(time = 1, id = ex6$id, neg = ex6$NegT1), data.frame(time = 2, id = ex6$id, neg = ex6$NegT2), data.frame(time = 3, id = ex6$id, neg = ex6$NegT3), data.frame(time = 4, id = ex6$id, neg = ex6$NegT4)) m6 <- lmer(neg ~ time + (time | id), data = ex6long, REML = FALSE) summary(m6)