Hi,
I am trying to run a Cholesky Decomposition with 3 binary variables and I keep getting the error message
“In model 'multiCholACEModel' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)”
I have tried changing the threshold starting values (the ones currently in there is from the saturated model) and reducing mvnRelEps, but the model still will not converge. I am not sure if this is due to errors in the script or other reasons.
Due to results from the univariate models, I went straight for the AE model.
Any guidance or tips on how to proceed would be greatly appreciated, I have included the script, but due to confidentiality agreements I cannot include any data.
Kind regards/ Lisa
#setting up for the analysis nv <- 3 #number of variables ntv <- nv*2 #number of variables for a pair nth <- 1 #number of thresholds thRows <- paste("th",1:nth,sep="") thFree <- matrix(c(T,T,T),nrow=nth,ncol=nv) thValues <- matrix(c(1.56,0.75,1.74),nrow=nth,ncol=nv) thLBound <- matrix(rep(c(0.001,rep(0.001,nth-1)),nv),nrow=nth,ncol=nv) #constricting the values to be positive thLabels=rep(c("pain","cmd","disability"),2) #Organising the data as ordinal Vars <- c('PAIN','CMD', 'DISABILITY') selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") mzDataOrd <- subset(selData, bestzyg==1, selVars) dzDataOrd <- subset(selData, bestzyg==2, selVars) mzDataOrd [,1]<- mxFactor(mzDataOrd[,1], levels = c(0:1)) mzDataOrd [,2]<- mxFactor(mzDataOrd[,2], levels = c(0:1)) dzDataOrd [,1]<- mxFactor(dzDataOrd[,1], levels = c(0:1)) dzDataOrd [,2]<- mxFactor(dzDataOrd[,2], levels = c(0:1)) mzDataOrd [,3]<- mxFactor(mzDataOrd[,3], levels = c(0:1)) mzDataOrd [,4]<- mxFactor(mzDataOrd[,4], levels = c(0:1)) dzDataOrd [,3]<- mxFactor(dzDataOrd[,3], levels = c(0:1)) dzDataOrd [,4]<- mxFactor(dzDataOrd[,4], levels = c(0:1)) mzDataOrd [,5]<- mxFactor(mzDataOrd[,5], levels = c(0:1)) mzDataOrd [,6]<- mxFactor(mzDataOrd[,6], levels = c(0:1)) dzDataOrd [,5]<- mxFactor(dzDataOrd[,5], levels = c(0:1)) dzDataOrd [,6]<- mxFactor(dzDataOrd[,6], levels = c(0:1)) summary(mzDataOrd) summary(dzDataOrd) describe(mzDataOrd) describe(dzDataOrd) # Fit Multivariate ACE Model with RawData and Matrices Input # ----------------------------------------------------------------------- cholSVnv <- c(.5,0,0,.5,0,.5) AFac <- c( "a11","a21","a31", "a22","a32", "a33") CFac <- c( "c11","c21","c31", "c22","c32", "c33") EFac <- c( "e11","e21","e31", "e22","e32", "e33") # Matrices a, c, and e to store a, c, and e path coefficients a <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholSVnv, labels=AFac, name="a" ) c <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=FALSE, values=0, labels=CFac, name="c" ) e <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholSVnv, labels=EFac, name="e" ) # Matrices A, C, and E compute variance components A <- mxAlgebra( expression=a %*% t(a), name="A" ) C <- mxAlgebra( expression=c %*% t(c), name="C" ) E <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) V <- mxAlgebra( expression=A+C+E, name="V" ) I <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") iSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Constraint on variance of Ordinal variables Var1 <- mxConstraint( expression=I*V==I, name="Var1") #if ordinal ## Note that the rest of the mxModel statements do not change for bi/multivariate case # Matrix & Algebra for expected means vector expMean<- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" ) Thre<- mxMatrix( type="Full", nrow=nth, ncol=nv, free=TRUE, values=thValues, lbound=thLBound, name="Thre" ) Inc<- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" ) ThreInc<- mxAlgebra( expression= Inc %*% Thre, name="ThreInc") expThre<- mxAlgebra( expression= cbind(ThreInc,ThreInc), dimnames=list(thRows,selVars), name="expThre" ) # Algebra for expected variance/covariance matrix in MZ expCovMZ<- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) # Algebra for expected variance/covariance matrix in DZ expCovDZ<- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) funML <- mxFitFunctionML() dataMZ <- mxData( observed=mzDataOrd, type="raw" ) expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre") dataDZ<- mxData( observed=dzDataOrd, type="raw" ) expDZ<- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre") pars <- list(a, c, e, A, C, E, V, I, iSD) modelMZ <- mxModel( expCovMZ, expMean, Thre, Inc, ThreInc, expThre, dataMZ, expMZ, funML, pars, name="MZ" ) modelDZ <- mxModel( expCovDZ, expMean, Thre, Inc, ThreInc, expThre, dataDZ, expDZ, funML, pars, name="DZ") multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) multiCholACEModel <- mxModel( "multiCholACEModel", modelMZ, modelDZ, multi) multiCholACEModel <- mxTryHard(multiCholACEModel) multiCholACEModelfit <- mxRun( multiCholACEModel, intervals=F ) summultiCholACEModel <- summary( multiCholACEModelfit ) summultiCholACEModel
Printer-friendly version
I corrected a non-ASCII quote in line 12 of your script, from
Vars <- c('PAIN','CMD’, 'DISABILITY')toVars <- c('PAIN','CMD', 'DISABILITY'). I could not run that line under Windows or Linux/GNU. Are you on a Mac?The algorithm that calculates the multivariate-normal probability integral only has so much numerical accuracy, and the inaccuracy is worse out into the tails of the distribution. Further, the numerical error is compounded by the finite-differences subroutine that numerically approximates the fitfunction derivatives. As a result, sometimes models of threshold traits cannot satisfy the first or second-order conditions for a local minimum, even when they've actually reached one. In such situations, the best you can do is try to find the smallest possible -2logL. I notice you're using
mxTryHard(). Some of the default argument values formxTryHard()are downright counterproductive for optimization with threshold traits. UsemxTryHardOrdinal()instead (or equivalently use different argument values to vanillamxTryHard()). Also, if you're not using the on-load default optimizer, CSOLNP, you should be.The issue of "inescapable status red with threshold variables" is quite common. Two relevant forums threads that readily come to my mind are this one and this one. You could find others if you search our forums.
Thank you so much, now that I changed over to mxTryHardOrdinal() it is finding a solution. I am on a PC; I must have accidentally added a space when copying over the script.
However, now that I got some results, I noticed that my model is estimating separate thresholds for MZ and DZ twins, even though they seem to have the same label. Am I doing something wrong in the mxFitFunctionMultigroup?
It's nothing related to
mxFitFunctionMultigroup(). Your syntax defines a character vectorthLabelsbut doesn't use it when you createThre. Consequently, each threshold exists in your MxModel's namespace as 'containingModelName.matrixName[row,col]', e.g. 'MZ.Thre[1,1]'.The issue of separate thresholds for MZs and DZs is that you don't have any labels in your
Thre:You could add labels.
Also consider that
umxACEwill do all of this automatically for you in one maintained and flexible function:Bates, T. C., Maes, H., & Neale, M. C. (2019). umx: Twin and Path-Based Structural Equation Modeling in R. Twin Res Hum Genet, 22(1), 27-41. doi:10.1017/thg.2019.2
Thank you so much for the help, I have added the labels in the OpenMx model and it now works, however since I am still getting code red error message, I thought I would give umx a go. When I fit a cholesky model it estimates the thresholds separately for twin 1 and twin 2 and I am getting the following parameters:
1 PAIN1_dev1
2 CMD1_dev1
3 DISABILITY1_dev1
4 PAIN2_dev1
5 CMD2_dev1
6 DISABILITY2_dev1
7 a_r1c1
8 a_r2c1
9 a_r3c1
10 a_r2c2
11 a_r3c2
12 a_r3c3
13 c_r1c1
14 c_r2c1
15 c_r3c1
16 c_r2c2
17 c_r3c2
18 c_r3c3
19 e_r1c1
20 e_r2c1
21 e_r3c1
22 e_r2c2
23 e_r3c2
24 e_r3c3
The extra thresholds are giving me problems when I try to fit an IP model. I tried using umxThresholdMatrix and it is telling me that the thresholds all have the same labels and same values.
> umxThresholdMatrix(dzDataOrd, selDVs = selVars, sep = "", method = c("auto"), threshMatName = "threshMat", droplevels = FALSE, verbose = FALSE)
[[1]]
LowerMatrix 'lowerOnes_for_thresh'
$labels: No labels assigned.
$values
[,1]
[1,] 1
$free: No free parameters.
$lbound: No lower bounds assigned.
$ubound: No upper bounds assigned.
[[2]]
FullMatrix 'deviations_for_thresh'
$labels
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 "PAIN_dev1" "CMD_dev1" "DISABILITY_dev1" "PAIN_dev1" "CMD_dev1" "DISABILITY_dev1"
$values
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 1.754467 1.367181 1.957871 1.754467 1.367181 1.957871
$free
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 TRUE TRUE TRUE TRUE TRUE TRUE
$lbound: No lower bounds assigned.
$ubound: No upper bounds assigned.
[[3]]
mxAlgebra 'threshMat'
$formula: lowerOnes_for_thresh %*% deviations_for_thresh
$result: (not yet computed) <0 x 0 matrix>
dimnames:
[[1]]
[1] "th_1"
[[2]]
[1] "PAIN1" "CMD1" "DISABILITY1" "PAIN2" "CMD2" "DISABILITY2"
>
> umxThresholdMatrix(mzDataOrd, selDVs = selVars, sep = "", method = c("auto"), threshMatName = "threshMat", droplevels = FALSE, verbose = FALSE)
[[1]]
LowerMatrix 'lowerOnes_for_thresh'
$labels: No labels assigned.
$values
[,1]
[1,] 1
$free: No free parameters.
$lbound: No lower bounds assigned.
$ubound: No upper bounds assigned.
[[2]]
FullMatrix 'deviations_for_thresh'
$labels
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 "PAIN_dev1" "CMD_dev1" "DISABILITY_dev1" "PAIN_dev1" "CMD_dev1" "DISABILITY_dev1"
$values
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 1.802673 1.333306 2.006859 1.802673 1.333306 2.006859
$free
PAIN1 CMD1 DISABILITY1 PAIN2 CMD2 DISABILITY2
dev_1 TRUE TRUE TRUE TRUE TRUE TRUE
$lbound: No lower bounds assigned.
$ubound: No upper bounds assigned.
[[3]]
mxAlgebra 'threshMat'
$formula: lowerOnes_for_thresh %*% deviations_for_thresh
$result: (not yet computed) <0 x 0 matrix>
dimnames:
[[1]]
[1] "th_1"
[[2]]
[1] "PAIN1" "CMD1" "DISABILITY1" "PAIN2" "CMD2" "DISABILITY2"
Am I missing something in my model? Below is the syntax I have been using.
nv <- 3 #number of variables
Vars <- c('PAIN','CMD','DISABILITY')
Vars
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")
selVars
mzDataOrd <- subset(selData, bestzyg==1, selVars)
dzDataOrd <- subset(selData, bestzyg==2, selVars)
mzDataOrd [,1]<- mxFactor(mzDataOrd[,1], levels = c(0:1))
mzDataOrd [,2]<- mxFactor(mzDataOrd[,2], levels = c(0:1))
dzDataOrd [,1]<- mxFactor(dzDataOrd[,1], levels = c(0:1))
dzDataOrd [,2]<- mxFactor(dzDataOrd[,2], levels = c(0:1))
mzDataOrd [,3]<- mxFactor(mzDataOrd[,3], levels = c(0:1))
mzDataOrd [,4]<- mxFactor(mzDataOrd[,4], levels = c(0:1))
dzDataOrd [,3]<- mxFactor(dzDataOrd[,3], levels = c(0:1))
dzDataOrd [,4]<- mxFactor(dzDataOrd[,4], levels = c(0:1))
mzDataOrd [,5]<- mxFactor(mzDataOrd[,5], levels = c(0:1))
mzDataOrd [,6]<- mxFactor(mzDataOrd[,6], levels = c(0:1))
dzDataOrd [,5]<- mxFactor(dzDataOrd[,5], levels = c(0:1))
dzDataOrd [,6]<- mxFactor(dzDataOrd[,6], levels = c(0:1))
summary(mzDataOrd)
summary(dzDataOrd)
umxThresholdMatrix(dzDataOrd, selDVs = selVars, sep = "", method = c("auto"), threshMatName = "threshMat", droplevels = FALSE, verbose = FALSE)
umxThresholdMatrix(mzDataOrd, selDVs = selVars, sep = "", method = c("auto"), threshMatName = "threshMat", droplevels = FALSE, verbose = FALSE)
Chol1 = umxACE("CholACEModel", selDVs = Vars, sep = "", mzData = mzDataOrd, dzData = dzDataOrd)
umxSummary(Chol1)
plot(Chol1)
parameters(Chol1)
Could you be a bit more explicit when you say "The extra thresholds are giving me problems when I try to fit an IP model. I tried using umxThresholdMatrix and it is telling me that the thresholds all have the same labels and same values"?
It is pretty much impossible to help when one doesn't know what exactly the problems are. However, good that you post the code which can sometimes give the intrepid support people clues:
Note how in the example part of ?umxThresholdMatrix its result is put into an object called tmp:
That object should be included in the model. Running it without a target object to put its output into simply dumps the results to the screen and does nothing else.
should, I think, include a tmp <- at the beginning, and then be used in the mxModel. But (and it's a big one) umxACE does all this for you by detecting which variables are ordinal and which are not.
Incidentally, it is possible to generate data that have similar properties to the real data, but which are definitely not the observed data and therefore can be shared to help others to debug. See ?mxGenerateData
Finally, if code 6 is the only concern, note that the message is often a false alarm caused by the precision of numerical integration making it difficult to find the exact solution - but the solution found may be extremely close to the correct one (no differences in 2nd decimal place of estimates or -2lnL). Rerunning from different starting values can increase confidence that the approximate solution is very close indeed to the actual one.
umxThresholdMatrix: you don't need to call it unless you're building a model from scratch.FYI:
umxThresholdMatrixbuilds an algebra and 2 matrices: one consisting of offsets (deviations_for_thresh) and the other a lower-1s matrixlowerOnes_for_threshwhich the algebra (threshMat) uses to add up the deviations to form thresholds). It thus handles the complex-for-humans job of building correct thresholds which are used within a model.Sorry if my post was unclear, I have not been able to get a model to run with generated data, but I have fitted a umx model with some public data from a VCU workshop and the same thing happened, the model seem to generate seperate thresholds for twin one and twin 2 and I am getting warning messages about these paths in the IP model, see output below: