Hello,
I have completed a standard ACE model of a phenotype with no errors. I now want to permute the zygosity labels of the twins and perform 10,000 permutations of the ACE model. I will then use this distribution of permitted A estimates to estimate the significance of my observed A.
When I run the analysis with the observed data I receive no warnings. However, when I run the permutations I usually receive about 20 of the warnings shown at the bottom of this message. Based upon the info in this forum https://openmx.ssri.psu.edu/thread/754 it seems like this error is most likely from when the permutations are operating on the extremes of the model and are safe to ignore.
However, I am not sure how to confirm this is the case with so many models to work with.
Thanks!
Brittany
In model 'ACE' 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)
Brittany,
When I've done simulation analyses with large numbers of models, I've kept track of the optimization status. If you have a model called
yourRunModel
that has been run, thenyourRunModel$output$status$code
gives the status code as a number. Zero means it converged fine. Six is the code you mentioned. Given that you're doing 10,000 permutations and only about 20 are giving status 6, I'd just compute the confidence interval after dropping the models that have status code 6. A more strict criteria would be dropping all models that have a non-zero status codes.Cheers,
Mike
Thanks for the great suggestions!
might be worth adding some code to mxRun() them again: often the code 6 will resolve on a second or third try.
Or better yet use
mxTryHard
Overall, I agree with Bates and Hunter. As I see it, you basically have four options:
1. During the permutation procedure, simply don't store the result when the optimizer status code is 6 (Red); or, just throw out code-6 results when the procedure is over. It probably won't matter much if it's only 20 out of 10,000.
2. By similar reasoning, just ignore the code 6's and use those results when calculating your permutation p-value. Again, it might not make much difference if it's just 20 out of 10,000.
3. You could re-write your permutation procedure's code so that, when the optimizer status code is 6, the result is not stored and the loop index is decremented by 1. For instance, if iteration number 4756 results in code 6, the procedure would just re-shuffle the zygosity labels and "try again" for iteration number 4756.
4. You could use
mxTryHard()
(with the appropriate arguments) instead ofmxRun()
, which should greatly cut down on the number of code 6's that occur.With regard to option 3 above, the loop index in R can be modified within a loop but any changes to it are thrown out when it changes to the next iteration.
Yeah. I should have specified to use a
while(i<=10000)
loop for this purpose.Thanks for laying those options out! I am interested in trying mxTryHard. If I am just running the ACE matrix tutorial script (http://openmx.psyc.virginia.edu/docs/openmx/latest/_static/demo/UnivariateTwinAnalysis_MatrixRaw.R) my understanding is I would substitute mxRun() to mxTryHard() in the two locations below(towards the end of the script) correct? Thanks!
Run Model
twinACEFit <- mxRun(twinACEModel, intervals=T)
twinACESum <- summary(twinACEFit)
twinACESum
Fit ACE Model with RawData and Matrices Input
-----------------------------------------------------------------------------
twinACEFit <- mxRun(twinACEModel)
Run ACE Model
-----------------------------------------------------------------------------
Basically, yes. But look at the documentation for
mxTryHard()
. Its default arguments might not be what you want. For instance,greenOK=TRUE
is probably fine for your purpose. Also, you might want to use fewerextraTries
than 10.Great, thanks!