library(lavaan)

dat5 <- read.csv("example5c.csv", header = FALSE)
colnames(dat5) <- c("u1", "u2", "u3", "u4", "u5", "u6", "u7", "u8", "g")
for(i in 1:8) dat5[,i] <- ordered(dat5[,i])
dat5[,9] <- factor(dat5[,9], labels = c("male", "female"))

configural5 <- "
f1 =~ c(1, 1)*u1 + u2 + u3 + u4
f2 =~ c(1, 1)*u5 + u6 + u7 + u8
u1 | c(t11, t11)*t1 + c(t12, t12)*t2 + t3 + t4
u2 | c(t21, t21)*t1 + t2 + t3 + t4
u3 | c(t31, t31)*t1 + t2 + t3 + t4
u4 | c(t41, t41)*t1 + t2 + t3 + t4
u5 | c(t51, t51)*t1 + c(t52, t52)*t2 + t3 + t4
u6 | c(t61, t61)*t1 + t2 + t3 + t4
u7 | c(t71, t71)*t1 + t2 + t3 + t4
u8 | c(t81, t81)*t1 + t2 + t3 + t4
f1 ~~ NA*f1
f2 ~~ NA*f2
f1 ~~ NA*f2
f1 ~ c(0, NA)*1
f2 ~ c(0, NA)*1
u1 ~~ c(1, NA)*u1
u2 ~~ c(1, NA)*u2
u3 ~~ c(1, NA)*u3
u4 ~~ c(1, NA)*u4
u5 ~~ c(1, NA)*u5
u6 ~~ c(1, NA)*u6
u7 ~~ c(1, NA)*u7
u8 ~~ c(1, NA)*u8
"

outConfigural5 <- cfa(configural5, data = dat5, group = "g", parameterization="theta", estimator="wlsmv")

weak5 <- "
f1 =~ c(1, 1)*u1 + c(f21, f21)*u2 + c(f31, f31)*u3 + c(f41, f41)*u4
f2 =~ c(1, 1)*u5 + c(f62, f62)*u6 + c(f72, f72)*u7 + c(f82, f82)*u8
u1 | c(t11, t11)*t1 + c(t12, t12)*t2 + t3 + t4
u2 | c(t21, t21)*t1 + t2 + t3 + t4
u3 | c(t31, t31)*t1 + t2 + t3 + t4
u4 | c(t41, t41)*t1 + t2 + t3 + t4
u5 | c(t51, t51)*t1 + c(t52, t52)*t2 + t3 + t4
u6 | c(t61, t61)*t1 + t2 + t3 + t4
u7 | c(t71, t71)*t1 + t2 + t3 + t4
u8 | c(t81, t81)*t1 + t2 + t3 + t4
f1 ~~ NA*f1
f2 ~~ NA*f2
f1 ~~ NA*f2
f1 ~ c(0, NA)*1
f2 ~ c(0, NA)*1
u1 ~~ c(1, NA)*u1
u2 ~~ c(1, NA)*u2
u3 ~~ c(1, NA)*u3
u4 ~~ c(1, NA)*u4
u5 ~~ c(1, NA)*u5
u6 ~~ c(1, NA)*u6
u7 ~~ c(1, NA)*u7
u8 ~~ c(1, NA)*u8
"

outWeak5 <- cfa(weak5, data = dat5, group = "g", parameterization="theta", estimator="wlsmv")

strong5 <- "
f1 =~ c(1, 1)*u1 + c(f21, f21)*u2 + c(f31, f31)*u3 + c(f41, f41)*u4
f2 =~ c(1, 1)*u5 + c(f62, f62)*u6 + c(f72, f72)*u7 + c(f82, f82)*u8
u1 | c(t11, t11)*t1 + c(t12, t12)*t2 + c(t13, t13)*t3 + c(t14, t14)*t4
u2 | c(t21, t21)*t1 + c(t22, t22)*t2 + c(t23, t23)*t3 + c(t24, t24)*t4
u3 | c(t31, t31)*t1 + c(t32, t32)*t2 + c(t33, t33)*t3 + c(t34, t34)*t4
u4 | c(t41, t41)*t1 + c(t42, t42)*t2 + c(t43, t43)*t3 + c(t44, t44)*t4
u5 | c(t51, t51)*t1 + c(t52, t52)*t2 + c(t53, t53)*t3 + c(t54, t54)*t4
u6 | c(t61, t61)*t1 + c(t62, t62)*t2 + c(t63, t63)*t3 + c(t64, t64)*t4
u7 | c(t71, t71)*t1 + c(t72, t72)*t2 + c(t73, t73)*t3 + c(t74, t74)*t4
u8 | c(t81, t81)*t1 + c(t82, t82)*t2 + c(t83, t83)*t3 + c(t84, t84)*t4
f1 ~~ NA*f1
f2 ~~ NA*f2
f1 ~~ NA*f2
f1 ~ c(0, NA)*1
f2 ~ c(0, NA)*1
u1 ~~ c(1, NA)*u1
u2 ~~ c(1, NA)*u2
u3 ~~ c(1, NA)*u3
u4 ~~ c(1, NA)*u4
u5 ~~ c(1, NA)*u5
u6 ~~ c(1, NA)*u6
u7 ~~ c(1, NA)*u7
u8 ~~ c(1, NA)*u8
"

outStrong5 <- cfa(strong5, data = dat5, group = "g", parameterization="theta", estimator="wlsmv")

strict5 <- "
f1 =~ c(1, 1)*u1 + c(f21, f21)*u2 + c(f31, f31)*u3 + c(f41, f41)*u4
f2 =~ c(1, 1)*u5 + c(f62, f62)*u6 + c(f72, f72)*u7 + c(f82, f82)*u8
u1 | c(t11, t11)*t1 + c(t12, t12)*t2 + c(t13, t13)*t3 + c(t14, t14)*t4
u2 | c(t21, t21)*t1 + c(t22, t22)*t2 + c(t23, t23)*t3 + c(t24, t24)*t4
u3 | c(t31, t31)*t1 + c(t32, t32)*t2 + c(t33, t33)*t3 + c(t34, t34)*t4
u4 | c(t41, t41)*t1 + c(t42, t42)*t2 + c(t43, t43)*t3 + c(t44, t44)*t4
u5 | c(t51, t51)*t1 + c(t52, t52)*t2 + c(t53, t53)*t3 + c(t54, t54)*t4
u6 | c(t61, t61)*t1 + c(t62, t62)*t2 + c(t63, t63)*t3 + c(t64, t64)*t4
u7 | c(t71, t71)*t1 + c(t72, t72)*t2 + c(t73, t73)*t3 + c(t74, t74)*t4
u8 | c(t81, t81)*t1 + c(t82, t82)*t2 + c(t83, t83)*t3 + c(t84, t84)*t4
f1 ~~ NA*f1
f2 ~~ NA*f2
f1 ~~ NA*f2
f1 ~ c(0, NA)*1
f2 ~ c(0, NA)*1
u1 ~~ c(1, 1)*u1
u2 ~~ c(1, 1)*u2
u3 ~~ c(1, 1)*u3
u4 ~~ c(1, 1)*u4
u5 ~~ c(1, 1)*u5
u6 ~~ c(1, 1)*u6
u7 ~~ c(1, 1)*u7
u8 ~~ c(1, 1)*u8
"

outStrict5 <- cfa(strict5, data = dat5, group = "g", parameterization="theta", estimator="wlsmv")

summary(outConfigural5, fit = TRUE)
summary(outWeak5, fit = TRUE)
summary(outStrong5, fit = TRUE)
summary(outStrict5, fit = TRUE)

anova(outConfigural5, outWeak5)
anova(outWeak5, outStrong5)
anova(outStrong5, outStrict5)