Cholesky with latent traits

Hello, I want to estimate a modified CP model, where I have a Cholesky model with three latent traits (see picture attached).
The d- and e-variables are ordinal (d = 2 categories; e = 3 categories). The s-variable is continuous and in order to integrate it into the LISREL-matrices I define a single factor measurement model with the residual variance set to 0 and the factor loading to 1.
For the ordinal variables I fix the means to zero and the residual variances to 1. I followed the approach proposed [here](https://openmx.ssri.psu.edu/thread/4127) to fix the residual variances of the categorical variables to one instead of the variance of the underlying continuous variable. So I can set the thresholds free.
Do I need an additional constraint for the second order factors?
I get the following error:
<%
All fit attempts resulted in errors - check starting values or model specification
%>
since I have tried to estimate the model with different starting values and I am using mxTryHardOrdinal() I suppose that the specification and not the starting values is the cause of the error.
I used WLS because of the high number of categorical variables: I tried the estimation with FIML but after waiting a long time without results I stopped the process. Since the s-variable has a relative high number of NAs, this may be an issue?

Here is my code:
<%
rm(list=ls())
library(OpenMx)
library(dplyr)

# Import Data
data <- read.csv(file = "data_wide.csv",
header = TRUE)

# Number of variables
m <- 6 # lat. end. V.
n <- 18 # lat. ex. V.
p <- 32 # Ind. of lat. end. V.
q <- 0 # Ind. of lat. ex. V.

# Vectors of variables
enLV <- c("E","D","S")
mlab <- c(paste0(enLV,1),paste0(enLV,2))
exLV <- c("A","C","E")
nlab <- c(paste0(exLV[1],enLV,1),paste0(exLV[2],enLV,1),paste0(exLV[3],enLV,1),
paste0(exLV[1],enLV,2),paste0(exLV[2],enLV,2),paste0(exLV[3],enLV,2))
plab <- c(paste0("e",1:5,"t",1),paste0("d",1:10,"t",1),"st1",paste0("e",1:5,"t",2),paste0("d",1:10,"t",2),"st2")

del <- c(plab[c(6:15,22:31)])
ext <- c(plab[c(1:5,17:21)])

mzData <- subset(data, zyg==1, plab)
dzData <- subset(data, zyg==2, plab)
# define categorical items as mxfactors
# 1. externalizing items (3 levels)
mzData[ext] <- mxFactor(mzData[ext], levels = 1:3)
dzData[ext] <- mxFactor(dzData[ext], levels = 1:3)

# 2. delinquency items (2 levels)
mzData[del] <- mxFactor(mzData[del], levels = 0:1)
dzData[del] <- mxFactor(dzData[del], levels = 0:1)

# Define the expected covariance matrix using hte LISREL approach
pathA <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE,
free = c(TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,
TRUE,TRUE,TRUE),
values = c(.5,0,0,
0,.5,0,
0,0,.5),
labels = c("a11",NA,NA,
"a21","a22",NA,
"a31","a32","a33"),
lbound = c(.0001,NA,NA,
-10,.0001,NA,
-10,-10,.0001),
name = "A")
pathC <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE,
free = c(TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,
TRUE,TRUE,TRUE),
values = c(.5,0,0,
0,.5,0,
0,0,.5),
labels = c("c11",NA,NA,
"c21","a22",NA,
"c31","c32","c33"),
lbound = c(.0001,NA,NA,
-10,.0001,NA,
-10,-10,.0001),
name = "C")
pathE <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE,
free = c(TRUE,FALSE,FALSE,
TRUE,TRUE,FALSE,
TRUE,TRUE,TRUE),
values = c(1,0,0,
0,1,0,
0,0,1),
labels = c("e11",NA,NA,
"e21","e22",NA,
"e31","e32","e33"),
lbound = c(.0001,NA,NA,
-10,.0001,NA,
-10,-10,.0001),
name = "E")
Null <- mxMatrix(type = "Zero", nrow = 3, ncol = 9, name = "Zeros")
Gamma <- mxAlgebra(rbind(cbind(A,C,E,Zeros),
cbind(Zeros,A,C,E)), name = "GA")

Beta <- mxMatrix(type = "Zero", nrow = m, ncol = m, name = "B")
Id <- mxMatrix(type = "Iden", nrow = m, ncol = m, name = "I")

LaExt <- mxMatrix(type = "Full", nrow = 5, ncol = 1, byrow = TRUE,
free = c(FALSE,TRUE,TRUE,TRUE,TRUE),
values = c(1,.8,.8,.8,.8),
labels = c("la_ex1","la_ex2","la_ex3","la_ex4","la_ex5"),
name = "LExt")
NullExt <- mxMatrix(type = "Zero", nrow = 16-5, ncol = 1, name = "NExt")

LaDel <- mxMatrix(type = "Full", nrow = 11, ncol = 1, byrow = TRUE,
free = c(FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),
values = c(1,.8,.8,.8,.8,.8,.8,.8,.8,.8,0),
labels = c("la_de1","la_de2","la_de3","la_de4","la_de5","la_de6","la_de7","la_de8","la_de9","la_de10", NA),
name = "LDel")
NullDel <- mxMatrix(type = "Zero", nrow = 5, ncol = 1, name = "NDel")

LaSt <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE,
values = 1, labels = "la_st", name = "LSt")
NullSt <- mxMatrix(type = "Zero", nrow = 16-1, ncol = 1, name = "NSt")

NullFull <- mxMatrix(type = "Zero", nrow = 16, ncol = 1, name = "NFull")

LaY <- mxAlgebra(expression = cbind(rbind(LExt, NExt, NFull),
rbind(NDel, LDel, NFull),
rbind(NSt, LSt, NFull),
rbind(NFull, LExt, NExt),
rbind(NFull, NDel, LDel),
rbind(NFull, NSt, LSt)),
name = "LY")
errorval <- rep(c(rep(0,15),0),2)
errorlab <- rep(c("eex1","eex2","eex3","eex4","eex5",
"ed1","ed2","ed3","ed4","ed5","ed6","ed7","ed8","ed9","ed10",
"es"),2)
Thee <- mxMatrix(type = "Diag", nrow = 32, ncol = 32, byrow = FALSE,
free = FALSE, values = errorval, lbound = 0,
labels = errorlab, name = "THE")

Psi <- mxMatrix(type = "Diag", nrow = 6, ncol = 6, free = FALSE, values = 0, lbound = 0,
labels = rep(c("eEx","eD","eS"),2), name = "PS")

PhiVar <- mxMatrix(type = "Iden", ncol = 9, nrow = 9, name = "PHV")
PhiCMZ <- mxMatrix(type = "Diag", ncol = 9, nrow = 9,
values = c(1,1,1,1,1,1,0,0,0), name = "PHCMZ")
PhiCDZ <- mxMatrix(type = "Diag", ncol = 9, nrow = 9,
values = c(.5,.5,.5,1,1,1,0,0,0), name = "PHCDZ")

PhiMZ <- mxAlgebra(rbind(cbind(PHV,PHCMZ),
cbind(PHCMZ,PHV)), name = "PHMZ")
PhiDZ <- mxAlgebra(rbind(cbind(PHV,PHCDZ),
cbind(PHCDZ,PHV)), name = "PHDZ")

eCovMZ <- mxAlgebra(expression= LY%*%solve(I-B)%*%(GA%*%PHMZ%*%t(GA)+PS)%*%t(solve(I-B))%*%t(LY)+THE, name="expCovMZ") # Formula: Bollen (1989)
eCovDZ <- mxAlgebra(expression= LY%*%solve(I-B)%*%(GA%*%PHDZ%*%t(GA)+PS)%*%t(solve(I-B))%*%t(LY)+THE, name="expCovDZ")
#eCovMZ <- mxAlgebra(expression= LY%*%I%*%(GA%*%PHMZ%*%t(GA)+PS)%*%I%*%t(LY)+THE, name="expCovMZ")
#eCovDZ <- mxAlgebra(expression= LY%*%I%*%(GA%*%PHDZ%*%t(GA)+PS)%*%I%*%t(LY)+THE, name="expCovDZ")

# Defining the mean matrix
# mean of isei
mean <- mxMatrix(type = "Full", nrow = 1, ncol = 32,
free = rep(c(rep(FALSE,15),TRUE),2), values = rep(c(rep(0,15),40),2),
labels = rep(c(rep(NA,15),"mst"),2),
name = "expMean")

# Defining the threshold matrix
# Example of logic of premultiplication of thinG with inc ! -> Ensures ordering of thresholds for items with >1 thresholds
#start <- matrix(c(rep(-1.5,30),rep(1,30)),nrow = 2, ncol = 30, byrow = TRUE)
#one <- matrix(c(1,0,1,1), nrow = 2, ncol = 2, byrow = TRUE)
#one%*%start

thDim <- plab[c(1:5, 17:21, 6:15, 22:31)]
thLab <- c(rep(paste0("th1","e",1:5),2),rep(paste0("th1","d",1:10),2),rep(paste0("th2","e",1:5),2),rep(NA,20))
thFree <- c(rep(TRUE,40),rep(FALSE,20))
thVal <- c(rep(-1.5,30),rep(1,10),rep(NA,20))
thBound <- c(rep(-3,30),rep(0.001,10),rep(NA,20))
thinG <- mxMatrix( type="Full", nrow=2, ncol=30, byrow = TRUE, free=thFree,
values=thVal, lbound=thBound,
labels=thLab, name="thinG" ) # Incremental matrix: 1 row = start values of 1. threshold; next rows = increments departing from upper row (lbound necessary so that there is at least a 0.0001 increase! nad no decrease to insure the order!)
inc <- mxMatrix( type="Lower", nrow=2, ncol=2, free=FALSE, values=1, name="inc" )
threG <- mxAlgebra( expression= inc %*% thinG, name="expTh" )

# Define data object
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Define expectation objects
expMZ <- mxExpectationNormal(covariance="expCovMZ", means="expMean",
thresholds="expTh",
dimnames=plab,
threshnames = thDim)
expDZ <- mxExpectationNormal(covariance="expCovDZ", means="expMean",
thresholds="expTh",
dimnames=plab,
threshnames = thDim)
# Fit function (WLS)
fitfun <- mxFitFunctionWLS()

# parameters
pars <- list(pathA, pathC, pathE, Null, Gamma, Beta, Id,
LaExt, NullExt, LaDel, NullDel, LaSt, NullSt, NullFull, LaY,
Thee, Psi, mean, thinG, inc, threG, PhiVar)

# group specific model objects
modelMZ <- mxModel(pars, PhiMZ, PhiCMZ, eCovMZ, dataMZ, expMZ, fitfun, name="MZ")
modelDZ <- mxModel(pars, PhiDZ, PhiCDZ, eCovDZ, dataDZ, expDZ, fitfun, name="DZ")
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )

# overall model object
modelACE <- mxModel( "ACE", pars, modelMZ, modelDZ, multi)

# run model
mxOption( NULL , 'Default optimizer' , 'SLSQP' )
set.seed(1)
fitACE <- mxTryHardOrdinal(modelACE, extraTries = 10, exhaustive = TRUE)
sumACE <- summary( fitACE )
sumACE

%>

Here is the output:
<%
Summary of ACE

Starting values are not feasible. Consider mxTryHard()

free parameters:
name matrix row col Estimate lbound ubound
1 a11 A 1 1 0.66719606 1e-04
2 a21 A 2 1 0.24639404 -10
3 a31 A 3 1 0.01989468 -10
4 a22 A 2 2 0.58462127 1e-04
5 a32 A 3 2 0.23233684 -10
6 a33 A 3 3 0.20895600 1e-04
7 c11 C 1 1 0.38163935 1e-04
8 c21 C 2 1 -0.23787426 -10
9 c31 C 3 1 0.05621162 -10
10 c32 C 3 2 0.09830437 -10
11 c33 C 3 3 0.65105809 1e-04
12 e11 E 1 1 0.78303864 1e-04
13 e21 E 2 1 -0.17010759 -10
14 e31 E 3 1 0.24441680 -10
15 e22 E 2 2 0.81938448 1e-04
16 e32 E 3 2 -0.15604409 -10
17 e33 E 3 3 1.10636594 1e-04
18 la_ex2 LExt 2 1 0.92273536
19 la_ex3 LExt 3 1 0.88631576
20 la_ex4 LExt 4 1 1.03505741
21 la_ex5 LExt 5 1 0.95816739
22 la_de2 LDel 2 1 0.54168182
23 la_de3 LDel 3 1 0.65248564
24 la_de4 LDel 4 1 1.02872796
25 la_de5 LDel 5 1 0.47255096
26 la_de6 LDel 6 1 0.58154706
27 la_de7 LDel 7 1 0.92095349
28 la_de8 LDel 8 1 0.83431902
29 la_de9 LDel 9 1 0.75920953
30 la_de10 LDel 10 1 0.74839136
31 mst expMean 1 16 48.11579465
32 th1e1 thinG 1 1 -1.72455967 -3
33 th2e1 thinG 2 1 0.88925538 0.001
34 th1e2 thinG 1 2 -1.31678786 -3
35 th2e2 thinG 2 2 1.29860129 0.001
36 th1e3 thinG 1 3 -1.97153396 -3
37 th2e3 thinG 2 3 0.81570670 0.001
38 th1e4 thinG 1 4 -1.31957163 -3
39 th2e4 thinG 2 4 1.23111741 0.001
40 th1e5 thinG 1 5 -1.09611506 -3
41 th2e5 thinG 2 5 1.07165598 0.001
42 th1d1 thinG 1 11 -0.91010120 -3
43 th1d2 thinG 1 12 -1.38949672 -3
44 th1d3 thinG 1 13 -1.66416074 -3
45 th1d4 thinG 1 14 -0.97191166 -3
46 th1d5 thinG 1 15 -1.64255135 -3
47 th1d6 thinG 1 16 -1.15400035 -3
48 th1d7 thinG 1 17 -1.11629872 -3
49 th1d8 thinG 1 18 -1.49572992 -3
50 th1d9 thinG 1 19 -1.97838685 -3
51 th1d10 thinG 1 20 -1.98501499 -3

Model Statistics:
| Parameters | Degrees of Freedom | Fit ( units)
Model: 51 62399 NA
Saturated: NA NA NA
Independence: NA NA NA
Number of observations/statistics: 2042/62450

timestamp: 2020-11-13 08:57:16
Wall clock time: 2.551282 secs
OpenMx version number: 2.18.1
Need help? See help(mxSummary)
%>