Recently I was trying a trivariate Cholesky model with the latest version of OpenMx. However, when I estimated the standardized CI of parameters, it reported “Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings, : NLOPT fatal error -1”. Then I followed the advices of Rob replied to another topic that change the optimizer from “SLSQP” to “NPSOL”. Although it run smoothly, the results contained minus value and bigger than 1 that were very strange, such as:
CholACE.StandardizedA[1,1] 4.835575e-02 0.4772359455 8.043944e-01
CholACE.StandardizedA[2,1] -3.001322e+03 1.6422971828 1.889364e+03
CholACE.StandardizedA[3,1] -4.236242e+02 -0.2316164474 1.479296e+00
The genetic correlation seemed also unreasonable:
CholACE.CorA[1,1] 1.000000e+00 1.0000000000 1.000000e+00 !!!
CholACE.CorA[2,1] -9.999988e-01 -0.9999978673 9.999902e-01
CholACE.CorA[3,1] -9.999983e-01 0.9999635762 9.999960e-01
I also checked my dataset which contained no singular value. I totally have no idea about this situation. Could anyone help me find the problems of my codes (as follow)? Many thanks!
meVals <- c(0.5,0.5,0.5) # start value for means meanLabs <- paste(Vars,"mean",sep="") # create labels for means paVal <- .6 # start value for path coefficient paLBo <- .0001 # start value for lower bounds paValD <- vech(diag(paVal,nv,nv)) # start values for diagonal of covariance matrix paLBoD <- diag(paLBo,nv,nv) # lower bounds for diagonal of covariance matrix paLBoD[lower.tri(paLBoD)] <- -10 # lower bounds for below diagonal elements paLBoD[upper.tri(paLBoD)] <- NA # lower bounds for above diagonal elements # Create Labels for Lower Triangular Matrices aLabs <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") cLabs <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") eLabs <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paValD, labels=aLabs, lbound=paLBoD, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paValD, labels=cLabs, lbound=paLBoD, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paValD, labels=eLabs, lbound=paLBoD, name="e" ) # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=meVals, labels=meanLabs, name="Mean" ) meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ) covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) # Standeraized parameters and CI sta <- mxAlgebra(iSD %*% a, name="sta") stc <- mxAlgebra(iSD %*% c, name="stc") ste <- mxAlgebra(iSD %*% e, name="ste") StandardizedA <- mxAlgebra(A/V, name="StandardizedA") StandardizedC <- mxAlgebra(C/V, name="StandardizedC") StandardizedE <- mxAlgebra(E/V, name="StandardizedE") CorA <- mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA") CorC <- mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC") CorE <- mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE") CorP <- mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP") CI <- mxCI(c('sta', 'stc', 'ste','StandardizedA', 'StandardizedC', 'StandardizedE','CorA', 'CorC', 'CorE', 'CorP')) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG, meanT, sta, stc, ste, StandardizedA, StandardizedC, StandardizedE, CorA, CorC, CorE, CorP) modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, funML, name="DZ" ) fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") ) CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, fitML, CI)