Adding thresholds to model using umxThresholdMatrix

Error: The thresholds matrix associated with the expectation function in m
odel 'MZ' contains column names and the expectation function has specified
 non-identical threshnames.

library(umx) 
library(dplyr)
library(ggplot2) 
 
data(docData)
 
mzData  = subset(docData, zygosity %in% c("MZFF", "MZMM")) %>%
  # remove all with B in the name
  select(-contains("B"), -zygosity) %>%
   mutate(varA1_T1 =  cut_number(varA1_T1, 4, ordered_result=T),
  varA1_T2 = cut_number(varA1_T2, 4, ordered_result = T)
  )
 
# same for dz
 
dzData  = subset(docData, zygosity %in% c("DZFF", "DZMM")) %>%
  # remove all with B in the name
  select(-contains("B"), -zygosity) %>%
   mutate(varA1_T1 =  cut_number(varA1_T1, 4, ordered_result=T),
  varA1_T2 = cut_number(varA1_T2, 4, ordered_result = T)
  )
 
 
top = mxModel("top",
  umxMatrix("BE", type = "Full",nrow=3,  ncol = 3,
    labels = c(NA, "g2", "b1", 
      "g1", NA, "b2", 
      NA, NA, NA),
    free = c(FALSE, FALSE, TRUE,
      TRUE, FALSE, TRUE,
      FALSE, FALSE, FALSE)),
  umxMatrix('I', type='Iden', nrow= 3,ncol= 3 ),
  umxMatrix('F', type='Full', nrow=5, ncol=6, free=FALSE,
    values=c(1,0,0,0,0,0,
      0,1,0,0,0,0,
      0,0,1,0,0,0,
      0,0,0,1,0,0,
      0,0,0,0,1,0)),
  umxMatrix('LY', type='Full',nrow=3, ncol = 3, free = FALSE, values = diag(3),
    labels = NA),
  umxMatrix('A', type='Symm', nrow=3,ncol = 3,
    labels=c("ab2","abraas",NA,
      "abraas","as2",NA,
      NA,NA,"x2"),
    free=c(TRUE,TRUE,FALSE,
      TRUE,TRUE,FALSE,
      FALSE,FALSE,TRUE)),
  umxMatrix('C', type='Symm',nrow=3, ncol = 3,
    labels =c("cb2","cbrccs",NA,
      "cbrccs","cs2",NA,
      NA,NA,NA),
    free=c(TRUE,TRUE,FALSE,
      TRUE,TRUE,FALSE,
      FALSE,FALSE,FALSE)),
  umxMatrix('E', type='Symm', nrow=3, ncol = 3,
    labels =c("eb2",NA,NA,
      NA,"es2",NA,
      NA,NA,NA),
    free= c(TRUE,FALSE,FALSE,
      FALSE,TRUE,FALSE,
      FALSE,FALSE,FALSE)),
  umxMatrix('dzmu', type='Full', nrow=1, ncol=6, free=TRUE, values= 0, labels=c('muPh1','muPh2','muPS1','muPh1','muPh2','muPS1')),
  mxAlgebra('mzmu', expression = dzmu%*%t(F)),
  mxAlgebra('A_'  , expression = solve(I - BE)%&%A),
  mxAlgebra('C_'  , expression = solve(I - BE)%&%C),
  mxAlgebra('E_'  , expression = solve(I - BE)%&%E),
  mxAlgebra('SPh' , expression = A_ + C_ + E_),
  mxAlgebra('Smz_', expression = rbind(
    cbind(SPh, A_+C_),
    cbind(A_+C_, SPh))),
  mxAlgebra('Sdz', expression=rbind(
    cbind(SPh, .5%x%A_+C_),
    cbind(.5%x%A_+C_, SPh))),
  mxAlgebra('Smz', expression= F%&%Smz_)
)
 
tymz <-   umxThresholdMatrix(mzData, fullVarNames = c("varA1_T1", "varA1_T2"),
                             sep = "_T", method ="allFree")
tydz <-   umxThresholdMatrix(dzData, fullVarNames = c("varA1_T1", "varA1_T2"), 
                             sep = "_T", method ="allFree")
 
MZ = mxModel("MZ",
  mxData(mzData[1:6], type="raw"),
  mxExpectationNormal(covariance = "top.Smz",means = "top.mzmu", dimnames =colnames(mzData[1:5]), thresholds="threshMat"),
  tymz,
  mxFitFunctionML()
)
 
DZ = mxModel("DZ",
  mxData(dzData[1:6], type="raw"),
  mxExpectationNormal(covariance = "top.Sdz",means = "top.dzmu", dimnames =colnames(mzData[1:6]),thresholds="threshMat"),
  tydz,
  mxFitFunctionML()
)
 
 
model = mxModel("mrdoc", top, MZ, DZ, mxFitFunctionMultigroup(c("MZ","DZ") ))
model <- omxAssignFirstParameters(model)
 
out <- mxRun(model)
summary(out)