### DEFINING NULL MODEL

mod.null <- '
x1 ~~ c(psi1, psi1)*x1
x2 ~~ c(psi2, psi2)*x2
x3 ~~ c(psi3, psi3)*x3
x4 ~~ c(psi4, psi4)*x4
x5 ~~ c(psi5, psi5)*x5
x6 ~~ c(psi6, psi6)*x6
x7 ~~ c(psi7, psi7)*x7
x8  ~~ c(psi8, psi8)*x8
x9  ~~ c(psi9, psi9)*x9
x10  ~~ c(psi10, psi10)*x10
x11 ~~ c(psi11, psi11)*x11
x12 ~~ c(psi12, psi12)*x12
x13 ~~ c(psi13, psi13)*x13
x14 ~~ c(psi14, psi14)*x14
x15 ~~ c(psi15, psi15)*x15
x16 ~~ c(psi16, psi16)*x16

x1 ~ c(tau1, tau1)*1
x2 ~ c(tau2, tau2)*1
x3 ~ c(tau3, tau3)*1
x4 ~ c(tau4, tau4)*1
x5 ~ c(tau5, tau5)*1
x6 ~ c(tau6, tau6)*1
x7 ~ c(tau7, tau7)*1
x8  ~ c(tau8, tau8)*1
x9  ~ c(tau9, tau9)*1
x10  ~ c(tau10, tau10)*1
x11 ~ c(tau11, tau11)*1
x12 ~ c(tau12, tau12)*1
x13 ~ c(tau13, tau13)*1
x14 ~ c(tau14, tau14)*1
x15 ~ c(tau15, tau15)*1
x16 ~ c(tau16, tau16)*1
'

### DEFINING BIFACTOR MODEL BASED ON PRIOR ANALYSES

bifactor.model<-'
g =~ x8 + x1 + x2 + x3 + x4 + x5 + x6 + x16 + x7 + x9 + x10 + x11 + x12 + x13 + x14 + x15
f1 =~ x3 + x5 + x4 + x1 + x2 + x6 + x16
f2 =~ x10 + x7 + x8 + x9
f3 =~ x12 + x14 + x15 + x11 + x13
'

### GENERATING AND EDITING SYNTAX FOR CONFIGURAL MODEL

bifactor.config.syntax <- measEq.syntax(configural.model = bifactor.model, data = CFA, ordered = TRUE, parameterization = "theta", estimator = "WLSMV", ID.fac = "UV", ID.cat = "Wu.Estabrook.2016", group = "Gender")
cat(as.character(bifactor.config.syntax))

mod.config <- '
## LOADINGS:

g =~ c(NA, NA)*x8 + c(lambda.1_1.g1, lambda.1_1.g2)*x8
g =~ c(NA, NA)*x1 + c(lambda.2_1.g1, lambda.2_1.g2)*x1
g =~ c(NA, NA)*x2 + c(lambda.3_1.g1, lambda.3_1.g2)*x2
g =~ c(NA, NA)*x3 + c(lambda.4_1.g1, lambda.4_1.g2)*x3
g =~ c(NA, NA)*x4 + c(lambda.5_1.g1, lambda.5_1.g2)*x4
g =~ c(NA, NA)*x5 + c(lambda.6_1.g1, lambda.6_1.g2)*x5
g =~ c(NA, NA)*x6 + c(lambda.7_1.g1, lambda.7_1.g2)*x6
g =~ c(NA, NA)*x16 + c(lambda.8_1.g1, lambda.8_1.g2)*x16
g =~ c(NA, NA)*x7 + c(lambda.9_1.g1, lambda.9_1.g2)*x7
g =~ c(NA, NA)*x9 + c(lambda.10_1.g1, lambda.10_1.g2)*x9
g =~ c(NA, NA)*x10 + c(lambda.11_1.g1, lambda.11_1.g2)*x10
g =~ c(NA, NA)*x11 + c(lambda.12_1.g1, lambda.12_1.g2)*x11
g =~ c(NA, NA)*x12 + c(lambda.13_1.g1, lambda.13_1.g2)*x12
g =~ c(NA, NA)*x13 + c(lambda.14_1.g1, lambda.14_1.g2)*x13
g =~ c(NA, NA)*x14 + c(lambda.15_1.g1, lambda.15_1.g2)*x14
g =~ c(NA, NA)*x15 + c(lambda.16_1.g1, lambda.16_1.g2)*x15
f1 =~ c(NA, NA)*x1 + c(lambda.2_2.g1, lambda.2_2.g2)*x1
f1 =~ c(NA, NA)*x2 + c(lambda.3_2.g1, lambda.3_2.g2)*x2
f1 =~ c(NA, NA)*x3 + c(lambda.4_2.g1, lambda.4_2.g2)*x3
f1 =~ c(NA, NA)*x4 + c(lambda.5_2.g1, lambda.5_2.g2)*x4
f1 =~ c(NA, NA)*x5 + c(lambda.6_2.g1, lambda.6_2.g2)*x5
f1 =~ c(NA, NA)*x6 + c(lambda.7_2.g1, lambda.7_2.g2)*x6
f1 =~ c(NA, NA)*x16 + c(lambda.8_2.g1, lambda.8_2.g2)*x16
f2 =~ c(NA, NA)*x8 + c(lambda.1_3.g1, lambda.1_3.g2)*x8
f2 =~ c(NA, NA)*x7 + c(lambda.9_3.g1, lambda.9_3.g2)*x7
f2 =~ c(NA, NA)*x9 + c(lambda.10_3.g1, lambda.10_3.g2)*x9
f2 =~ c(NA, NA)*x10 + c(lambda.11_3.g1, lambda.11_3.g2)*x10
f3 =~ c(NA, NA)*x11 + c(lambda.12_4.g1, lambda.12_4.g2)*x11
f3 =~ c(NA, NA)*x12 + c(lambda.13_4.g1, lambda.13_4.g2)*x12
f3 =~ c(NA, NA)*x13 + c(lambda.14_4.g1, lambda.14_4.g2)*x13
f3 =~ c(NA, NA)*x14 + c(lambda.15_4.g1, lambda.15_4.g2)*x14
f3 =~ c(NA, NA)*x15 + c(lambda.16_4.g1, lambda.16_4.g2)*x15

## THRESHOLDS:

x8 | c(NA, NA)*t1 + c(x8.thr1.g1, x8.thr1.g2)*t1
x8 | c(NA, NA)*t2 + c(x8.thr2.g1, x8.thr2.g2)*t2
x8 | c(NA, NA)*t3 + c(x8.thr3.g1, x8.thr3.g2)*t3
x8 | c(NA, NA)*t4 + c(x8.thr4.g1, x8.thr4.g2)*t4
x1 | c(NA, NA)*t1 + c(x1.thr1.g1, x1.thr1.g2)*t1
x1 | c(NA, NA)*t2 + c(x1.thr2.g1, x1.thr2.g2)*t2
x1 | c(NA, NA)*t3 + c(x1.thr3.g1, x1.thr3.g2)*t3
x1 | c(NA, NA)*t4 + c(x1.thr4.g1, x1.thr4.g2)*t4
x2 | c(NA, NA)*t1 + c(x2.thr1.g1, x2.thr1.g2)*t1
x2 | c(NA, NA)*t2 + c(x2.thr2.g1, x2.thr2.g2)*t2
x2 | c(NA, NA)*t3 + c(x2.thr3.g1, x2.thr3.g2)*t3
x2 | c(NA, NA)*t4 + c(x2.thr4.g1, x2.thr4.g2)*t4
x3 | c(NA, NA)*t1 + c(x3.thr1.g1, x3.thr1.g2)*t1
x3 | c(NA, NA)*t2 + c(x3.thr2.g1, x3.thr2.g2)*t2
x3 | c(NA, NA)*t3 + c(x3.thr3.g1, x3.thr3.g2)*t3
x3 | c(NA, NA)*t4 + c(x3.thr4.g1, x3.thr4.g2)*t4
x4 | c(NA, NA)*t1 + c(x4.thr1.g1, x4.thr1.g2)*t1
x4 | c(NA, NA)*t2 + c(x4.thr2.g1, x4.thr2.g2)*t2
x4 | c(NA, NA)*t3 + c(x4.thr3.g1, x4.thr3.g2)*t3
x4 | c(NA, NA)*t4 + c(x4.thr4.g1, x4.thr4.g2)*t4
x5 | c(NA, NA)*t1 + c(x5.thr1.g1, x5.thr1.g2)*t1
x5 | c(NA, NA)*t2 + c(x5.thr2.g1, x5.thr2.g2)*t2
x5 | c(NA, NA)*t3 + c(x5.thr3.g1, x5.thr3.g2)*t3
x5 | c(NA, NA)*t4 + c(x5.thr4.g1, x5.thr4.g2)*t4
x6 | c(NA, NA)*t1 + c(x6.thr1.g1, x6.thr1.g2)*t1
x6 | c(NA, NA)*t2 + c(x6.thr2.g1, x6.thr2.g2)*t2
x6 | c(NA, NA)*t3 + c(x6.thr3.g1, x6.thr3.g2)*t3
x6 | c(NA, NA)*t4 + c(x6.thr4.g1, x6.thr4.g2)*t4
x16 | c(NA, NA)*t1 + c(x16.thr1.g1, x16.thr1.g2)*t1
x16 | c(NA, NA)*t2 + c(x16.thr2.g1, x16.thr2.g2)*t2
x16 | c(NA, NA)*t3 + c(x16.thr3.g1, x16.thr3.g2)*t3
x16 | c(NA, NA)*t4 + c(x16.thr4.g1, x16.thr4.g2)*t4
x7 | c(NA, NA)*t1 + c(x7.thr1.g1, x7.thr1.g2)*t1
x7 | c(NA, NA)*t2 + c(x7.thr2.g1, x7.thr2.g2)*t2
x7 | c(NA, NA)*t3 + c(x7.thr3.g1, x7.thr3.g2)*t3
x7 | c(NA, NA)*t4 + c(x7.thr4.g1, x7.thr4.g2)*t4
x9 | c(NA, NA)*t1 + c(x9.thr1.g1, x9.thr1.g2)*t1
x9 | c(NA, NA)*t2 + c(x9.thr2.g1, x9.thr2.g2)*t2
x9 | c(NA, NA)*t3 + c(x9.thr3.g1, x9.thr3.g2)*t3
x9 | c(NA, NA)*t4 + c(x9.thr4.g1, x9.thr4.g2)*t4
x10 | c(NA, NA)*t1 + c(x10.thr1.g1, x10.thr1.g2)*t1
x10 | c(NA, NA)*t2 + c(x10.thr2.g1, x10.thr2.g2)*t2
x10 | c(NA, NA)*t3 + c(x10.thr3.g1, x10.thr3.g2)*t3
x10 | c(NA, NA)*t4 + c(x10.thr4.g1, x10.thr4.g2)*t4
x11 | c(NA, NA)*t1 + c(x11.thr1.g1, x11.thr1.g2)*t1
x11 | c(NA, NA)*t2 + c(x11.thr2.g1, x11.thr2.g2)*t2
x11 | c(NA, NA)*t3 + c(x11.thr3.g1, x11.thr3.g2)*t3
x11 | c(NA, NA)*t4 + c(x11.thr4.g1, x11.thr4.g2)*t4
x12 | c(NA, NA)*t1 + c(x12.thr1.g1, x12.thr1.g2)*t1
x12 | c(NA, NA)*t2 + c(x12.thr2.g1, x12.thr2.g2)*t2
x12 | c(NA, NA)*t3 + c(x12.thr3.g1, x12.thr3.g2)*t3
x12 | c(NA, NA)*t4 + c(x12.thr4.g1, x12.thr4.g2)*t4
x13 | c(NA, NA)*t1 + c(x13.thr1.g1, x13.thr1.g2)*t1
x13 | c(NA, NA)*t2 + c(x13.thr2.g1, x13.thr2.g2)*t2
x13 | c(NA, NA)*t3 + c(x13.thr3.g1, x13.thr3.g2)*t3
x13 | c(NA, NA)*t4 + c(x13.thr4.g1, x13.thr4.g2)*t4
x14 | c(NA, NA)*t1 + c(x14.thr1.g1, x14.thr1.g2)*t1
x14 | c(NA, NA)*t2 + c(x14.thr2.g1, x14.thr2.g2)*t2
x14 | c(NA, NA)*t3 + c(x14.thr3.g1, x14.thr3.g2)*t3
x14 | c(NA, NA)*t4 + c(x14.thr4.g1, x14.thr4.g2)*t4
x15 | c(NA, NA)*t1 + c(x15.thr1.g1, x15.thr1.g2)*t1
x15 | c(NA, NA)*t2 + c(x15.thr2.g1, x15.thr2.g2)*t2
x15 | c(NA, NA)*t3 + c(x15.thr3.g1, x15.thr3.g2)*t3
x15 | c(NA, NA)*t4 + c(x15.thr4.g1, x15.thr4.g2)*t4

## INTERCEPTS:

x8 ~ c(nu.1.g1, nu.1.g2)*1 + c(0, 0)*1
x1 ~ c(nu.2.g1, nu.2.g2)*1 + c(0, 0)*1
x2 ~ c(nu.3.g1, nu.3.g2)*1 + c(0, 0)*1
x3 ~ c(nu.4.g1, nu.4.g2)*1 + c(0, 0)*1
x4 ~ c(nu.5.g1, nu.5.g2)*1 + c(0, 0)*1
x5 ~ c(nu.6.g1, nu.6.g2)*1 + c(0, 0)*1
x6 ~ c(nu.7.g1, nu.7.g2)*1 + c(0, 0)*1
x16 ~ c(nu.8.g1, nu.8.g2)*1 + c(0, 0)*1
x7 ~ c(nu.9.g1, nu.9.g2)*1 + c(0, 0)*1
x9 ~ c(nu.10.g1, nu.10.g2)*1 + c(0, 0)*1
x10 ~ c(nu.11.g1, nu.11.g2)*1 + c(0, 0)*1
x11 ~ c(nu.12.g1, nu.12.g2)*1 + c(0, 0)*1
x12 ~ c(nu.13.g1, nu.13.g2)*1 + c(0, 0)*1
x13 ~ c(nu.14.g1, nu.14.g2)*1 + c(0, 0)*1
x14 ~ c(nu.15.g1, nu.15.g2)*1 + c(0, 0)*1
x15 ~ c(nu.16.g1, nu.16.g2)*1 + c(0, 0)*1

## UNIQUE-FACTOR VARIANCES:

x8 ~~ c(1, 1)*x8 + c(theta.1_1.g1, theta.1_1.g2)*x8
x1 ~~ c(1, 1)*x1 + c(theta.2_2.g1, theta.2_2.g2)*x1
x2 ~~ c(1, 1)*x2 + c(theta.3_3.g1, theta.3_3.g2)*x2
x3 ~~ c(1, 1)*x3 + c(theta.4_4.g1, theta.4_4.g2)*x3
x4 ~~ c(1, 1)*x4 + c(theta.5_5.g1, theta.5_5.g2)*x4
x5 ~~ c(1, 1)*x5 + c(theta.6_6.g1, theta.6_6.g2)*x5
x6 ~~ c(1, 1)*x6 + c(theta.7_7.g1, theta.7_7.g2)*x6
x16 ~~ c(1, 1)*x16 + c(theta.8_8.g1, theta.8_8.g2)*x16
x7 ~~ c(1, 1)*x7 + c(theta.9_9.g1, theta.9_9.g2)*x7
x9 ~~ c(1, 1)*x9 + c(theta.10_10.g1, theta.10_10.g2)*x9
x10 ~~ c(1, 1)*x10 + c(theta.11_11.g1, theta.11_11.g2)*x10
x11 ~~ c(1, 1)*x11 + c(theta.12_12.g1, theta.12_12.g2)*x11
x12 ~~ c(1, 1)*x12 + c(theta.13_13.g1, theta.13_13.g2)*x12
x13 ~~ c(1, 1)*x13 + c(theta.14_14.g1, theta.14_14.g2)*x13
x14 ~~ c(1, 1)*x14 + c(theta.15_15.g1, theta.15_15.g2)*x14
x15 ~~ c(1, 1)*x15 + c(theta.16_16.g1, theta.16_16.g2)*x15

## LATENT MEANS/INTERCEPTS:

g ~ c(alpha.1.g1, alpha.1.g2)*1 + c(NA, NA)*1
f1 ~ c(alpha.2.g1, alpha.2.g2)*1 + c(NA, NA)*1
f2 ~ c(alpha.3.g1, alpha.3.g2)*1 + c(NA, NA)*1
f3 ~ c(alpha.4.g1, alpha.4.g2)*1 + c(NA, NA)*1

## COMMON-FACTOR VARIANCES:

g ~~ c(1, 1)*g + c(psi.1_1.g1, psi.1_1.g2)*g
f1 ~~ c(1, 1)*f1 + c(psi.2_2.g1, psi.2_2.g2)*f1
f2 ~~ c(1, 1)*f2 + c(psi.3_3.g1, psi.3_3.g2)*f2
f3 ~~ c(1, 1)*f3 + c(psi.4_4.g1, psi.4_4.g2)*f3

## COMMON-FACTOR COVARIANCES:

g ~~ c(NA, NA)*f1 + c(psi.2_1.g1, psi.2_1.g2)*f1
g ~~ c(NA, NA)*f2 + c(psi.3_1.g1, psi.3_1.g2)*f2
g ~~ c(NA, NA)*f3 + c(psi.4_1.g1, psi.4_1.g2)*f3
f1 ~~ c(NA, NA)*f2 + c(psi.3_2.g1, psi.3_2.g2)*f2
f1 ~~ c(NA, NA)*f3 + c(psi.4_2.g1, psi.4_2.g2)*f3
f2 ~~ c(NA, NA)*f3 + c(psi.4_3.g1, psi.4_3.g2)*f3
'

### GENERATING AND EDITING SYNTAX FOR THRESHOLD-CONSTRAINED MODEL

bifactor.thresh.syntax <- measEq.syntax(configural.model = bifactor.model, data = CFA, ordered = TRUE, parameterization = "theta", estimator = "WLSMV", ID.fac = "UV", ID.cat = "Wu.Estabrook.2016", group = "Gender", group.equal = "thresholds")
cat(as.character(bifactor.thresh.syntax))

mod.thresh <- '
## LOADINGS:

g =~ c(NA, NA)*x8 + c(lambda.1_1.g1, lambda.1_1.g2)*x8
g =~ c(NA, NA)*x1 + c(lambda.2_1.g1, lambda.2_1.g2)*x1
g =~ c(NA, NA)*x2 + c(lambda.3_1.g1, lambda.3_1.g2)*x2
g =~ c(NA, NA)*x3 + c(lambda.4_1.g1, lambda.4_1.g2)*x3
g =~ c(NA, NA)*x4 + c(lambda.5_1.g1, lambda.5_1.g2)*x4
g =~ c(NA, NA)*x5 + c(lambda.6_1.g1, lambda.6_1.g2)*x5
g =~ c(NA, NA)*x6 + c(lambda.7_1.g1, lambda.7_1.g2)*x6
g =~ c(NA, NA)*x16 + c(lambda.8_1.g1, lambda.8_1.g2)*x16
g =~ c(NA, NA)*x7 + c(lambda.9_1.g1, lambda.9_1.g2)*x7
g =~ c(NA, NA)*x9 + c(lambda.10_1.g1, lambda.10_1.g2)*x9
g =~ c(NA, NA)*x10 + c(lambda.11_1.g1, lambda.11_1.g2)*x10
g =~ c(NA, NA)*x11 + c(lambda.12_1.g1, lambda.12_1.g2)*x11
g =~ c(NA, NA)*x12 + c(lambda.13_1.g1, lambda.13_1.g2)*x12
g =~ c(NA, NA)*x13 + c(lambda.14_1.g1, lambda.14_1.g2)*x13
g =~ c(NA, NA)*x14 + c(lambda.15_1.g1, lambda.15_1.g2)*x14
g =~ c(NA, NA)*x15 + c(lambda.16_1.g1, lambda.16_1.g2)*x15
f1 =~ c(NA, NA)*x1 + c(lambda.2_2.g1, lambda.2_2.g2)*x1
f1 =~ c(NA, NA)*x2 + c(lambda.3_2.g1, lambda.3_2.g2)*x2
f1 =~ c(NA, NA)*x3 + c(lambda.4_2.g1, lambda.4_2.g2)*x3
f1 =~ c(NA, NA)*x4 + c(lambda.5_2.g1, lambda.5_2.g2)*x4
f1 =~ c(NA, NA)*x5 + c(lambda.6_2.g1, lambda.6_2.g2)*x5
f1 =~ c(NA, NA)*x6 + c(lambda.7_2.g1, lambda.7_2.g2)*x6
f1 =~ c(NA, NA)*x16 + c(lambda.8_2.g1, lambda.8_2.g2)*x16
f2 =~ c(NA, NA)*x8 + c(lambda.1_3.g1, lambda.1_3.g2)*x8
f2 =~ c(NA, NA)*x7 + c(lambda.9_3.g1, lambda.9_3.g2)*x7
f2 =~ c(NA, NA)*x9 + c(lambda.10_3.g1, lambda.10_3.g2)*x9
f2 =~ c(NA, NA)*x10 + c(lambda.11_3.g1, lambda.11_3.g2)*x10
f3 =~ c(NA, NA)*x11 + c(lambda.12_4.g1, lambda.12_4.g2)*x11
f3 =~ c(NA, NA)*x12 + c(lambda.13_4.g1, lambda.13_4.g2)*x12
f3 =~ c(NA, NA)*x13 + c(lambda.14_4.g1, lambda.14_4.g2)*x13
f3 =~ c(NA, NA)*x14 + c(lambda.15_4.g1, lambda.15_4.g2)*x14
f3 =~ c(NA, NA)*x15 + c(lambda.16_4.g1, lambda.16_4.g2)*x15

## THRESHOLDS:

x8 | c(NA, NA)*t1 + c(x8.thr1, x8.thr1)*t1
x8 | c(NA, NA)*t2 + c(x8.thr2, x8.thr2)*t2
x8 | c(NA, NA)*t3 + c(x8.thr3, x8.thr3)*t3
x8 | c(NA, NA)*t4 + c(x8.thr4, x8.thr4)*t4
x1 | c(NA, NA)*t1 + c(x1.thr1, x1.thr1)*t1
x1 | c(NA, NA)*t2 + c(x1.thr2, x1.thr2)*t2
x1 | c(NA, NA)*t3 + c(x1.thr3, x1.thr3)*t3
x1 | c(NA, NA)*t4 + c(x1.thr4, x1.thr4)*t4
x2 | c(NA, NA)*t1 + c(x2.thr1, x2.thr1)*t1
x2 | c(NA, NA)*t2 + c(x2.thr2, x2.thr2)*t2
x2 | c(NA, NA)*t3 + c(x2.thr3, x2.thr3)*t3
x2 | c(NA, NA)*t4 + c(x2.thr4, x2.thr4)*t4
x3 | c(NA, NA)*t1 + c(x3.thr1, x3.thr1)*t1
x3 | c(NA, NA)*t2 + c(x3.thr2, x3.thr2)*t2
x3 | c(NA, NA)*t3 + c(x3.thr3, x3.thr3)*t3
x3 | c(NA, NA)*t4 + c(x3.thr4, x3.thr4)*t4
x4 | c(NA, NA)*t1 + c(x4.thr1, x4.thr1)*t1
x4 | c(NA, NA)*t2 + c(x4.thr2, x4.thr2)*t2
x4 | c(NA, NA)*t3 + c(x4.thr3, x4.thr3)*t3
x4 | c(NA, NA)*t4 + c(x4.thr4, x4.thr4)*t4
x5 | c(NA, NA)*t1 + c(x5.thr1, x5.thr1)*t1
x5 | c(NA, NA)*t2 + c(x5.thr2, x5.thr2)*t2
x5 | c(NA, NA)*t3 + c(x5.thr3, x5.thr3)*t3
x5 | c(NA, NA)*t4 + c(x5.thr4, x5.thr4)*t4
x6 | c(NA, NA)*t1 + c(x6.thr1, x6.thr1)*t1
x6 | c(NA, NA)*t2 + c(x6.thr2, x6.thr2)*t2
x6 | c(NA, NA)*t3 + c(x6.thr3, x6.thr3)*t3
x6 | c(NA, NA)*t4 + c(x6.thr4, x6.thr4)*t4
x16 | c(NA, NA)*t1 + c(x16.thr1, x16.thr1)*t1
x16 | c(NA, NA)*t2 + c(x16.thr2, x16.thr2)*t2
x16 | c(NA, NA)*t3 + c(x16.thr3, x16.thr3)*t3
x16 | c(NA, NA)*t4 + c(x16.thr4, x16.thr4)*t4
x7 | c(NA, NA)*t1 + c(x7.thr1, x7.thr1)*t1
x7 | c(NA, NA)*t2 + c(x7.thr2, x7.thr2)*t2
x7 | c(NA, NA)*t3 + c(x7.thr3, x7.thr3)*t3
x7 | c(NA, NA)*t4 + c(x7.thr4, x7.thr4)*t4
x9 | c(NA, NA)*t1 + c(x9.thr1, x9.thr1)*t1
x9 | c(NA, NA)*t2 + c(x9.thr2, x9.thr2)*t2
x9 | c(NA, NA)*t3 + c(x9.thr3, x9.thr3)*t3
x9 | c(NA, NA)*t4 + c(x9.thr4, x9.thr4)*t4
x10 | c(NA, NA)*t1 + c(x10.thr1, x10.thr1)*t1
x10 | c(NA, NA)*t2 + c(x10.thr2, x10.thr2)*t2
x10 | c(NA, NA)*t3 + c(x10.thr3, x10.thr3)*t3
x10 | c(NA, NA)*t4 + c(x10.thr4, x10.thr4)*t4
x11 | c(NA, NA)*t1 + c(x11.thr1, x11.thr1)*t1
x11 | c(NA, NA)*t2 + c(x11.thr2, x11.thr2)*t2
x11 | c(NA, NA)*t3 + c(x11.thr3, x11.thr3)*t3
x11 | c(NA, NA)*t4 + c(x11.thr4, x11.thr4)*t4
x12 | c(NA, NA)*t1 + c(x12.thr1, x12.thr1)*t1
x12 | c(NA, NA)*t2 + c(x12.thr2, x12.thr2)*t2
x12 | c(NA, NA)*t3 + c(x12.thr3, x12.thr3)*t3
x12 | c(NA, NA)*t4 + c(x12.thr4, x12.thr4)*t4
x13 | c(NA, NA)*t1 + c(x13.thr1, x13.thr1)*t1
x13 | c(NA, NA)*t2 + c(x13.thr2, x13.thr2)*t2
x13 | c(NA, NA)*t3 + c(x13.thr3, x13.thr3)*t3
x13 | c(NA, NA)*t4 + c(x13.thr4, x13.thr4)*t4
x14 | c(NA, NA)*t1 + c(x14.thr1, x14.thr1)*t1
x14 | c(NA, NA)*t2 + c(x14.thr2, x14.thr2)*t2
x14 | c(NA, NA)*t3 + c(x14.thr3, x14.thr3)*t3
x14 | c(NA, NA)*t4 + c(x14.thr4, x14.thr4)*t4
x15 | c(NA, NA)*t1 + c(x15.thr1, x15.thr1)*t1
x15 | c(NA, NA)*t2 + c(x15.thr2, x15.thr2)*t2
x15 | c(NA, NA)*t3 + c(x15.thr3, x15.thr3)*t3
x15 | c(NA, NA)*t4 + c(x15.thr4, x15.thr4)*t4

## INTERCEPTS:

x8 ~ c(nu.1.g1, nu.1.g2)*1 + c(0, NA)*1
x1 ~ c(nu.2.g1, nu.2.g2)*1 + c(0, NA)*1
x2 ~ c(nu.3.g1, nu.3.g2)*1 + c(0, NA)*1
x3 ~ c(nu.4.g1, nu.4.g2)*1 + c(0, NA)*1
x4 ~ c(nu.5.g1, nu.5.g2)*1 + c(0, NA)*1
x5 ~ c(nu.6.g1, nu.6.g2)*1 + c(0, NA)*1
x6 ~ c(nu.7.g1, nu.7.g2)*1 + c(0, NA)*1
x16 ~ c(nu.8.g1, nu.8.g2)*1 + c(0, NA)*1
x7 ~ c(nu.9.g1, nu.9.g2)*1 + c(0, NA)*1
x9 ~ c(nu.10.g1, nu.10.g2)*1 + c(0, NA)*1
x10 ~ c(nu.11.g1, nu.11.g2)*1 + c(0, NA)*1
x11 ~ c(nu.12.g1, nu.12.g2)*1 + c(0, NA)*1
x12 ~ c(nu.13.g1, nu.13.g2)*1 + c(0, NA)*1
x13 ~ c(nu.14.g1, nu.14.g2)*1 + c(0, NA)*1
x14 ~ c(nu.15.g1, nu.15.g2)*1 + c(0, NA)*1
x15 ~ c(nu.16.g1, nu.16.g2)*1 + c(0, NA)*1

## UNIQUE-FACTOR VARIANCES:

x8 ~~ c(1, NA)*x8 + c(theta.1_1.g1, theta.1_1.g2)*x8
x1 ~~ c(1, NA)*x1 + c(theta.2_2.g1, theta.2_2.g2)*x1
x2 ~~ c(1, NA)*x2 + c(theta.3_3.g1, theta.3_3.g2)*x2
x3 ~~ c(1, NA)*x3 + c(theta.4_4.g1, theta.4_4.g2)*x3
x4 ~~ c(1, NA)*x4 + c(theta.5_5.g1, theta.5_5.g2)*x4
x5 ~~ c(1, NA)*x5 + c(theta.6_6.g1, theta.6_6.g2)*x5
x6 ~~ c(1, NA)*x6 + c(theta.7_7.g1, theta.7_7.g2)*x6
x16 ~~ c(1, NA)*x16 + c(theta.8_8.g1, theta.8_8.g2)*x16
x7 ~~ c(1, NA)*x7 + c(theta.9_9.g1, theta.9_9.g2)*x7
x9 ~~ c(1, NA)*x9 + c(theta.10_10.g1, theta.10_10.g2)*x9
x10 ~~ c(1, NA)*x10 + c(theta.11_11.g1, theta.11_11.g2)*x10
x11 ~~ c(1, NA)*x11 + c(theta.12_12.g1, theta.12_12.g2)*x11
x12 ~~ c(1, NA)*x12 + c(theta.13_13.g1, theta.13_13.g2)*x12
x13 ~~ c(1, NA)*x13 + c(theta.14_14.g1, theta.14_14.g2)*x13
x14 ~~ c(1, NA)*x14 + c(theta.15_15.g1, theta.15_15.g2)*x14
x15 ~~ c(1, NA)*x15 + c(theta.16_16.g1, theta.16_16.g2)*x15

## LATENT MEANS/INTERCEPTS:

g ~ c(alpha.1.g1, alpha.1.g2)*1 + c(NA, NA)*1
f1 ~ c(alpha.2.g1, alpha.2.g2)*1 + c(NA, NA)*1
f2 ~ c(alpha.3.g1, alpha.3.g2)*1 + c(NA, NA)*1
f3 ~ c(alpha.4.g1, alpha.4.g2)*1 + c(NA, NA)*1

## COMMON-FACTOR VARIANCES:

g ~~ c(1, 1)*g + c(psi.1_1.g1, psi.1_1.g2)*g
f1 ~~ c(1, 1)*f1 + c(psi.2_2.g1, psi.2_2.g2)*f1
f2 ~~ c(1, 1)*f2 + c(psi.3_3.g1, psi.3_3.g2)*f2
f3 ~~ c(1, 1)*f3 + c(psi.4_4.g1, psi.4_4.g2)*f3

## COMMON-FACTOR COVARIANCES:

g ~~ c(NA, NA)*f1 + c(psi.2_1.g1, psi.2_1.g2)*f1
g ~~ c(NA, NA)*f2 + c(psi.3_1.g1, psi.3_1.g2)*f2
g ~~ c(NA, NA)*f3 + c(psi.4_1.g1, psi.4_1.g2)*f3
f1 ~~ c(NA, NA)*f2 + c(psi.3_2.g1, psi.3_2.g2)*f2
f1 ~~ c(NA, NA)*f3 + c(psi.4_2.g1, psi.4_2.g2)*f3
f2 ~~ c(NA, NA)*f3 + c(psi.4_3.g1, psi.4_3.g2)*f3
'

### FITTING CONFIGURAL AND THRESHOLD-CONSTRAINED MODELS

bifactor.config <- cfa(mod.config, data = CFA, group = "Gender",ordered = TRUE, estimator = "WLSMV", parameterization = "theta")
bifactor.thresh <- cfa(mod.thresh, data = CFA, group = "Gender",ordered = TRUE, estimator = "WLSMV", parameterization = "theta")

### RUNNING PERMUTATIONAL TEST FOR THRESHOLD INVARIANCE

myAFIs <- c("chisq","chisq.scaled","df","cfi.scaled","tli.scaled","rmsea.scaled","srmr")
moreAFIs <- NULL # c("gammaHat","gammaHat.scaled")

bifactor.thresh.test <- permuteMeasEq(nPermute = 1000, uncon = bifactor.config, con = bifactor.thresh, param = "thresholds", AFIs = myAFIs, moreAFIs = moreAFIs, null = fit.null, parallelType = "none", iseed = 3141593)