# Code Red Error message NPSOL

8 posts / 0 new
Offline
Joined: 06/02/2010 - 22:22
Code Red Error message NPSOL

I am running a series of univariate models before I combine my variables into a trivariate cholesky. I am running a full ACE model as well as nested AE, CE and E models and my AE is best fitting (not sure if this is important but I am only requesting CIs for my AE model since it is best-fitting). The following error messages consistently pop up under the output for AE. The output for the full ACE and other nested models does not show this error message.

1: In model 'univAE' while computing the lower bound for 'univAE.stPathEst_a[1,1]' NPSOL 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)
2: In model 'univAE' while computing the lower bound for 'univAE.stPathEst_e[1,1]' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).
3: In model 'univAE' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).

I saw a previous post for the same error message and Dr. Neale stated this may mean the model is underidentified (as evidenced by identical -2LL). So, I changed the start values for all the models and came up with more of the same error messages. The best option (I think) is when I keep the start value for the AE .1 and the rest at 0. This gives me different -2LL and only one error message (the same one):

1: In model 'univAE' while computing the lower bound for 'univAE.stCovComp_E[1,1]' NPSOL 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)
2: In model 'univAE' NPSOL returned a non-zero status code 1. The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).

My (very limited beginner) understanding is that the code greens are ok, and that some code reds are ok as well. But I do not know how to distinguish whether this code red is ok for me to proceed with the rest of my analysis? Any help would be appreciated.

Thank you!

Offline
Joined: 07/31/2009 - 15:12
While NPSOL 6 errors are

While NPSOL 6 errors are generally something to worry about, they are of little concern for confidence intervals. We expect poor convergence for the bounds of confidence intervals; they're not minima of the objective function, so they don't technically converge. However, the change in -2LL when you change starting values indicates that at least one set of starting values yield a local rather than global minimum. You need to try multiple (more than two) sets of starting values and find the sets that yield the lowest -2LL. If you post your code, we'd have a better chance of diagnosing the problem.

Offline
Joined: 06/02/2010 - 22:22

Thanks for your message, it is very helpful. My code is below, the starting values are all at 0 for this version of the code.

require(OpenMx)
source("GenEpiHelperFunctions.R")

# -----------------------------------------------------------------------

require(psych)
summary (FamData)
describe(FamData)
Famvars <- colnames (FamData)
Vars <- c ('CB')
nv <- 1
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('CB1', 'CB2')
ntv <- nv*2
mzData <- subset(FamData, zyg==1, selVars)
dzData <- subset(FamData, zyg==2, selVars)

# -----------------------------------------------------------------------

describe(mzData)
colMeans(mzData,na.rm=TRUE)
cov(mzData,use="complete")
describe(dzData)
colMeans(dzData,na.rm=TRUE)
cov(dzData,use="complete")

# -----------------------------------------------------------------------

univTwinSatModel <- mxModel("univTwinSat",
mxModel("MZ",
mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="CholMZ" ),
mxAlgebra( expression=CholMZ %% t(CholMZ), name="expCovMZ" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=20, name="expMeanMZ" ),
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars),
# Algebra's needed for equality constraints
mxAlgebra( expression=expMeanMZ[1,1:nv], name="expMeanMZt1"),
mxAlgebra( expression=expMeanMZ[1,(nv+1):ntv], name="expMeanMZt2"),
mxAlgebra( expression=t(diag2vec(expCovMZ)), name="expVarMZ"),
mxAlgebra( expression=expVarMZ[1,1:nv], name="expVarMZt1"),
mxAlgebra( expression=expVarMZ[1,(nv+1):ntv], name="expVarMZt2")
),
mxModel("DZ",
mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="CholDZ" ),
mxAlgebra( expression=CholDZ %
% t(CholDZ), name="expCovDZ" ),
mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=20, name="expMeanDZ" ),
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars),
# Algebra's needed for equality constraints
mxAlgebra( expression=expMeanDZ[1,1:nv], name="expMeanDZt1"),
mxAlgebra( expression=expMeanDZ[1,(nv+1):ntv], name="expMeanDZt2"),
mxAlgebra( expression=t(diag2vec(expCovDZ)), name="expVarDZ"),
mxAlgebra( expression=expVarDZ[1,1:nv], name="expVarDZt1"),
mxAlgebra( expression=expVarDZ[1,(nv+1):ntv], name="expVarDZt2")
),
mxAlgebra( MZ.objective + DZ.objective, name="minus2sumll" ),
mxAlgebraObjective("minus2sumll")
)

univTwinSatFit <- mxRun(univTwinSatModel)
univTwinSatSumm <- summary(univTwinSatFit)
univTwinSatSumm

# -----------------------------------------------------------------------

parameterSpecifications(univTwinSatFit)
expectedMeansCovariances(univTwinSatFit)
tableFitStatistics(univTwinSatFit)

# -----------------------------------------------------------------------

univACEModel <- mxModel("univACE",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ), #X
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ), #Y
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ), #Z
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %% t(a), name="A" ),
mxAlgebra( expression=c %
% t(c), name="C" ),
mxAlgebra( expression=e %% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I
V)), name="iSD"),
## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, label="mean", name="Mean" ),
mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
# Algebra for expected variance/covariance matrix in MZ
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, note use of 0.5, converted to 1*1 matrix
mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
),
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars )
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumll" ),
mxAlgebraObjective("minus2sumll")
)

univACEFit <- mxRun(univACEModel)
univACESumm <- summary(univACEFit)
univACESumm

# -----------------------------------------------------------------------

parameterSpecifications(univACEFit)
expectedMeansCovariances(univACEFit)
tableFitStatistics(univACEFit)

# Generate Table of Parameter Estimates using mxEval

pathEstimatesACE <- mxEval(cbind(ACE.a,ACE.c,ACE.e), univACEFit)
rownames(pathEstimatesACE) <- 'pathEstimates'
colnames(pathEstimatesACE) <- c('a','c','e')
pathEstimatesACE

varComponentsACE <- mxEval(cbind(ACE.A/ACE.V,ACE.C/ACE.V,ACE.E/ACE.V), univACEFit)
rownames(varComponentsACE) <- 'varComponents'
colnames(varComponentsACE) <- c('a^2','c^2','e^2')
varComponentsACE

# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices

ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %% ACE.a","ACE.iSD %% ACE.c","ACE.iSD %*% ACE.e")
ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","iSD","stPathEst_a","stPathEst_c","stPathEst_e")
formatOutputMatrices(univACEFit,ACEpathMatrices,ACEpathLabels,Vars,4)

ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V")
ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E")
formatOutputMatrices(univACEFit,ACEcovMatrices,ACEcovLabels,Vars,4)

# -----------------------------------------------------------------------

univAEModel <- mxModel(univACEFit, name="univAE",
mxModel(univACEFit$ACE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="c11", name="c" ) # drop c at 0 ) ) univAEFit <- mxRun(univAEModel) univAESumm <- summary(univAEFit) # Create mxAlgebras for standardized estimates (e.g., A11) stCovComp_A <- mxAlgebra((ACE.A/ACE.V), name="stCovComp_A") stCovComp_E <- mxAlgebra((ACE.E/ACE.V), name="stCovComp_E") # Add estimates to model and ask the model to report CIs when it runs CImodel <- mxModel(univAEFit, stCovComp_A, stCovComp_E, mxCI(c("stCovComp_A", "stCovComp_E"))) CImodelfit<-mxRun(CImodel,intervals=TRUE) CIAEsumm<- summary(CImodelfit) CIAEsumm # Creat mxAlgebras for individaul standardized path estimates (e.g., a11) stPathEst_a <- mxAlgebra((ACE.iSD %% ACE.a), name="stPathEst_a") stPathEst_e <- mxAlgebra((ACE.iSD % % ACE.e), name="stPathEst_e") # Add estimates to model and ask the model to report CIs when it runs CImodel <- mxModel(univAEFit, stPathEst_a, stPathEst_e, mxCI(c("stPathEst_a", "stPathEst_e"))) CImodelfit<-mxRun(CImodel,intervals=TRUE) CIAEsumm<- summary(CImodelfit) CIAEsumm # Fit CE model # ----------------------------------------------------------------------- univCEModel <- mxModel(univACEFit, name="univCE", mxModel(univACEFit$ACE,
mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) # drop a at 0
)
)
univCEFit <- mxRun(univCEModel)
univCESumm <- summary(univCEFit)

# -----------------------------------------------------------------------

univEModel <- mxModel(univAEFit, name="univE",
mxModel(univAEFit$ACE, mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) # drop a at 0 ) ) univEFit <- mxRun(univEModel) univESumm <- summary(univEFit) # Print Comparative Fit Statistics # ----------------------------------------------------------------------- univACENested <- list(univAEFit, univCEFit, univEFit) tableFitStatistics(univACEFit,univACENested) Thanks again! Offline Joined: 06/02/2010 - 22:22 As an addition to my previous As an addition to my previous reply. I managed to eliminate the code red by changing the value from 0 to .1 in this part of the script. This seemed to have solved my error problem, but because I am a beginner in Mx, can I trust my results? # Fit AE model # ----------------------------------------------------------------------- univAEModel <- mxModel(univACEFit, name="univAE", mxModel(univACEFit$ACE,
mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=.1, label="c11", name="c" ) # drop c at 0
)
)
univAEFit <- mxRun(univAEModel)
univAESumm <- summary(univAEFit)

# Create mxAlgebras for standardized estimates (e.g., A11)

stCovComp_A <- mxAlgebra((ACE.A/ACE.V), name="stCovComp_A")
stCovComp_E <- mxAlgebra((ACE.E/ACE.V), name="stCovComp_E")

# Add estimates to model and ask the model to report CIs when it runs

CImodel <- mxModel(univAEFit, stCovComp_A, stCovComp_E, mxCI(c("stCovComp_A", "stCovComp_E")))
CImodelfit<-mxRun(CImodel,intervals=TRUE)
CIAEsumm<- summary(CImodelfit)
CIAEsumm

# Creat mxAlgebras for individaul standardized path estimates (e.g., a11)

stPathEst_a <- mxAlgebra((ACE.iSD %% ACE.a), name="stPathEst_a")
stPathEst_e <- mxAlgebra((ACE.iSD %
% ACE.e), name="stPathEst_e")

# Add estimates to model and ask the model to report CIs when it runs

CImodel <- mxModel(univAEFit, stPathEst_a, stPathEst_e, mxCI(c("stPathEst_a", "stPathEst_e")))
CImodelfit<-mxRun(CImodel,intervals=TRUE)
CIAEsumm<- summary(CImodelfit)
CIAEsumm

Offline
Joined: 07/31/2009 - 14:25
oops, no: you are now not

oops, no: you are now not setting C to 0, but to .1

mxMatrix( type="Lower", nrow=1, ncol=1, free=F, values=.1, label="c11", name="c" ) # drop c at 0

That's not dropping C, but just fixing it.

Offline
Joined: 06/02/2010 - 22:22
Ok, thanks, that is very

Ok, thanks, that is very helpful. Any suggestions on the NPSOL error messages or this one:
1: In model 'univAE' while computing the upper bound for 'univAE.stCovComp_E[1,1]' NPSOL 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)

It sounds like it may not be anything to worry about, although it says code red?

Offline
Joined: 07/31/2009 - 15:10
A code 6 (status RED) is a

A code 6 (status RED) is a warning that the optimizer thinks it might have gotten stuck in a local minimum.

Oftentimes, the solution is correct, but there's a little bit of wiggle room at the bottom of the likelihood well, and the optimizer isn't confident enough to call it a global minimum.

I usually recommend that people change their starting values and run the model again. If you get the same result from several sets of starting values, you can be pretty confident that this is the global maximum.