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!

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.

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")

# Prepare Data

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

FamData require(psych)

summary (FamData)

describe(FamData)

Famvars Vars nv selVars ntv mzData dzData

# Print Descriptive Statistics

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

describe(mzData)

colMeans(mzData,na.rm=TRUE)

cov(mzData,use="complete")

describe(dzData)

colMeans(dzData,na.rm=TRUE)

cov(dzData,use="complete")

# Fit Univariate Saturated Model

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

univTwinSatModel 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 univTwinSatSumm univTwinSatSumm

# Generate Saturated Output

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

parameterSpecifications(univTwinSatFit)

expectedMeansCovariances(univTwinSatFit)

tableFitStatistics(univTwinSatFit)

# Fit ACE Model with RawData and Matrices Input

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

univACEModel 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 univACESumm univACESumm

# Generate ACE Output

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

parameterSpecifications(univACEFit)

expectedMeansCovariances(univACEFit)

tableFitStatistics(univACEFit)

# Generate Table of Parameter Estimates using mxEval

pathEstimatesACE rownames(pathEstimatesACE) colnames(pathEstimatesACE) pathEstimatesACE

varComponentsACE rownames(varComponentsACE) colnames(varComponentsACE) varComponentsACE

# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices

ACEpathMatrices ACEpathLabels formatOutputMatrices(univACEFit,ACEpathMatrices,ACEpathLabels,Vars,4)

ACEcovMatrices ACEcovLabels formatOutputMatrices(univACEFit,ACEcovMatrices,ACEcovLabels,Vars,4)

# Fit AE model

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

univAEModel mxModel(univACEFit$ACE,

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

)

)

univAEFit univAESumm

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

stCovComp_A stCovComp_E

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

CImodel CImodelfit CIAEsumm CIAEsumm

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

stPathEst_a stPathEst_e

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

CImodel CImodelfit CIAEsumm CIAEsumm

# Fit CE model

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

univCEModel mxModel(univACEFit$ACE,

mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) # drop a at 0

)

)

univCEFit univCESumm

# Fit E model

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

univEModel mxModel(univAEFit$ACE,

mxMatrix( type="Lower", nrow=1, ncol=1, free=FALSE, values=0, label="a11", name="a" ) # drop a at 0

)

)

univEFit univESumm

# Print Comparative Fit Statistics

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

univACENested tableFitStatistics(univACEFit,univACENested)

Thanks again!

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$ACE,

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

)

)

univAEFit univAESumm

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

stCovComp_A stCovComp_E

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

CImodel CImodelfit CIAEsumm CIAEsumm

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

stPathEst_a stPathEst_e

CImodel CImodelfit CIAEsumm CIAEsumm

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.

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?

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.

There's more info on interpreting NPSOL status codes in this thread:

http://openmx.psyc.virginia.edu/thread/448

And in this one:

http://openmx.psyc.virginia.edu/thread/287

I think there's an overhaul of these error messages on the list of things to do, but it may be a little while before we get there.

How about trying to increase default the optimal accuracy for the case inform 6?