# Example of Multiple Group Codes in lavaan ##################### Hard way ########################## library(lavaan) # Configural Invariance model.configural <- ' # No constraint across groups # Factor loadings visual =~ c(load11, load12)*x1 + c(load21, load22)*x2 + c(load31, load32)*x3 textual =~ c(load41, load42)*x4 + c(load51, load52)*x5 + c(load61, load62)*x6 speed =~ c(load71, load72)*x7 + c(load81, load82)*x8 + c(load91, load92)*x9 # Intercepts x1 ~ c(int11, int12)*1 x2 ~ c(int21, int22)*1 x3 ~ c(int31, int32)*1 x4 ~ c(int41, int42)*1 x5 ~ c(int51, int52)*1 x6 ~ c(int61, int62)*1 x7 ~ c(int71, int72)*1 x8 ~ c(int81, int82)*1 x9 ~ c(int91, int92)*1 # Factor variances visual ~~ c(1, 1)*visual textual ~~ c(1, 1)*textual speed ~~ c(1, 1)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, 0)*1 textual ~ c(0, 0)*1 speed ~ c(0, 0)*1 ' configural <- cfa(model.configural, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Weak Invariance model.weak <- ' # Constrain Factor Loading # Free the second group factor variances # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int12)*1 x2 ~ c(int21, int22)*1 x3 ~ c(int31, int32)*1 x4 ~ c(int41, int42)*1 x5 ~ c(int51, int52)*1 x6 ~ c(int61, int62)*1 x7 ~ c(int71, int72)*1 x8 ~ c(int81, int82)*1 x9 ~ c(int91, int92)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, 0)*1 textual ~ c(0, 0)*1 speed ~ c(0, 0)*1 ' weak <- cfa(model.weak, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Weak invariance testing anova(configural, weak) # Strong Invariance model.strong <- ' # Constrain Measurement Intercepts # Free the second group factor means # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' strong <- cfa(model.strong, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Strong invariance testing anova(weak, strong) # Factor Variance Invariance model.varin <- ' # Constrain All Factor Variance to be Equal (fix all as 1) # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, 1)*visual textual ~~ c(1, 1)*textual speed ~~ c(1, 1)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' varin <- cfa(model.varin, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Omnibus test of invariance anova(strong, varin) # Test factor 1 variance model.varin1 <- ' # Constrain Only first factor variance # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, 1)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' varin1 <- cfa(model.varin1, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Test of variance 1. Do similary for Factors 2 and 3 anova(strong, varin1) # Covariance Invariance model.covin <- ' # Constrain All Covariances to be equal # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov11)*textual visual ~~ c(cov21, cov21)*speed textual ~~ c(cov31, cov31)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' covin <- cfa(model.covin, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Covariance invariance testing anova(strong, covin) # Test covariance between visual and textual model.covin1 <- ' # Constrain Covariance 1 to be equal # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov11)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' covin1 <- cfa(model.covin1, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Covariance invariance testing anova(strong, covin1) # Test factor mean invariance model.meanin <- ' # Constrain All Factor Means to be Equal (Fix All Factor Means as 0) # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, 0)*1 textual ~ c(0, 0)*1 speed ~ c(0, 0)*1 ' meanin <- cfa(model.meanin, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Factor means invariance testing anova(strong, meanin) # Test Factor Mean 1 model.meanin1 <- ' # Constrain The First Factor Mean to Be Equal (as 0) # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual textual ~~ c(1, NA)*textual speed ~~ c(1, NA)*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, 0)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 ' meanin1 <- cfa(model.meanin1, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Test for Factor Mean 1 difference. Do similarly for the other means anova(strong, meanin1) # Test for correlation difference # Lavaan currently (0.5-*) does not good for phantom variable approach for comparing correlations. # Use nonlinear constraints for comparing correlations instead model.corin <- ' # Label all factor variances # The factor covaraince of the first group is correlation already because the factor variances are 1 # Change factor covariances of the second group to correlation # Make a difference in correlation # Factor loadings visual =~ c(load11, load11)*x1 + c(load21, load21)*x2 + c(load31, load31)*x3 textual =~ c(load41, load41)*x4 + c(load51, load51)*x5 + c(load61, load61)*x6 speed =~ c(load71, load71)*x7 + c(load81, load81)*x8 + c(load91, load91)*x9 # Intercepts x1 ~ c(int11, int11)*1 x2 ~ c(int21, int21)*1 x3 ~ c(int31, int31)*1 x4 ~ c(int41, int41)*1 x5 ~ c(int51, int51)*1 x6 ~ c(int61, int61)*1 x7 ~ c(int71, int71)*1 x8 ~ c(int81, int81)*1 x9 ~ c(int91, int91)*1 # Factor variances visual ~~ c(1, NA)*visual + label("","var12")*visual textual ~~ c(1, NA)*textual + label("","var22")*textual speed ~~ c(1, NA)*speed + label("","var32")*speed # Factor covariances visual ~~ c(cov11, cov12)*textual visual ~~ c(cov21, cov22)*speed textual ~~ c(cov31, cov32)*speed # Factor means visual ~ c(0, NA)*1 textual ~ c(0, NA)*1 speed ~ c(0, NA)*1 # Factor correlations cor11 := cov11 cor21 := cov21 cor31 := cov31 cor12 := cov12 / sqrt(var12 * var22) cor22 := cov22 / sqrt(var12 * var32) cor32 := cov32 / sqrt(var22 * var32) # Factor correlation difference diff1 := cor11 - cor12 diff2 := cor21 - cor22 diff3 := cor31 - cor32 ' corin <- cfa(model.corin, data=HolzingerSwineford1939, group="school", std.lv=TRUE) # Use bootstrap for testing factor correlation difference corinBoot <- cfa(model.corin, data=HolzingerSwineford1939, group="school", se="boot", std.lv=TRUE) parameterEstimates(corinBoot, boot.ci.type="bca.simple")