Dear all,
I have fitted a simplex model with six time points. It's a threshold model with 2 thresholds. I fixed the thresholds and freely estimated the means. I used the CSOLNP optimizer.
I get stable estimates, they are in accordance with what I would expect based on the univariate results, and the gradients look ok.
As a next step, I wanted to calculate confidence intervals, but I only get this:
confidence intervals:
lbound estimate ubound note
atm21 NA 0.4871751 NA !!!
I had many problems with fitting threshold models in openMx before, so I was wondering whether calculating those CIs is a realistic plan anyway, and whether there are any tricks to get the model with CIs run properly?
Many thanks in advance!
Kind regards
Charlotte
Printer-friendly version
In general, you'll be better off using SLSQP for confidence intervals.
Dear Joshua and Rob,
Thanks a lot for your fast replies!
The script is now running with SLSQP, let's see what happens (I don't know whether this optimizer can deal with the thresholds?).
I would also like to try Rob's suggestion, but I'm not sure how to implement it into my script (see below)?
Best regards,
Charlotte
Combine groups
pars <- list(THm, add, pathAtm, pathAim, pathCtm, pathCim, pathEm, covAm, covCm, covEm, varM, matI, corAm, corCm, corEm, corVm, fitFunction )
modelMZM <- mxModel( pars, meanMM, CovMZM, ThresMM, dataMZM, objMZM, name = "MZM" )
modelDZM <- mxModel( pars, meanMM, CovDZM, ThresMM, dataDZM, objDZM, name = "DZM" )
Objective
minus2ll <- mxAlgebra( expression=MZM.objective + DZM.objective, name = "m2LL" )
obj <- mxFitFunctionAlgebra( "m2LL" )
Conf <- mxCI( reference="atm21", interval=0.99 )
plan <- mxComputeSequence(steps=list(
mxComputeGradientDescent(engine="CSOLNP"),
mxComputeConfidenceInterval(plan=mxComputeGradientDescent(engine="SLSQP")),
mxComputeStandardError(),
mxComputeHessianQuality(),
mxComputeReportDeriv(),
mxComputeReportExpectation()
))
SimplexACEModel <- mxModel( "SimplexACE", pars, modelMZM, modelDZM, minus2ll, obj, relativeAm, relativeCm, relativeEm, Conf, plan )
------------------------------------------------------------------------------
RUN SIMPLEX ACE MODEL
mxOption( NULL,"Default optimizer","CSOLNP")
mxOption( NULL,"mvnAbsEps",1.e-6 )
mxOption( NULL,"mvnMaxPointsC",5000 )
mxOption( NULL,"Major iterations",2000 )
SimplexACEFit <- mxTryHard(SimplexACEModel, extraTries=5, intervals=T)
Error in .local(.Object, ...) :
could not find function "mxComputeReportExpectation"
Calls: mxComputeSequence -> new -> initialize -> initialize -> .local
Execution halted
Hmm, I forgot that
mxComputeReportExpectation()isn't in the most recent version of OpenMx. It WILL be in the next version released, though, so try changingplanto this:Also, consider using
mxTryHardOrdinal()instead ofmxTryHard(). The former is just a wrapper for the latter, but the default values of its arguments are tailored specifically for ordinal analyses.Dear Rob,
Thanks again for your fast reply! I have adapted the script according to your suggestions, I will keep you updated on whether it works!
I will also use mxTryHardOrdinal for threshold models from now on, thank you for this advice!
Best regards
Charlotte
Glad I could be of help.
I should mention that the biggest difference between the defaults for
mxTryHard()andmxTryHardOrdinal()is that argumentfinetuneGradientdefaults toTRUEformxTryHard()butFALSEformxTryHardOrdinal(). Like the help page says, "finetuneGradient=FALSEis recommended for analyses involving thresholds."Also, for the benefit of any lurkers reading this thread: it's quite possible to just use SLSQP as the default optimizer and not worry about a custom compute plan, in which case OpenMx will automatically use SLSQP for finding point estimates and confidence limits.
I second Joshua's suggestion of not using CSOLNP for confidence intervals, which are still one of that optimizer's weak points. If you still want to use CSOLNP to find point estimates, you could try a custom compute plan, which would look something like
Then, you would just put
planinto your MxModel, just like you would an MxMatrix or something like that.If you post the script you're using, I could tailor the compute-plan syntax to it.