You are here

Simplex models with thresholds - CIs?

7 posts / 0 new
Last post
Charlotte's picture
Offline
Joined: 07/02/2012 - 11:13
Simplex models with thresholds - CIs?

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

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
CIs

In general, you'll be better off using SLSQP for confidence intervals.

Charlotte's picture
Offline
Joined: 07/02/2012 - 11:13
How to implement Rob's suggestion?

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
try this

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 changing plan to this:

plan <- mxComputeSequence(
    steps=list(
    mxComputeGradientDescent(engine="CSOLNP"),
    mxComputeConfidenceInterval(plan=mxComputeGradientDescent(engine="SLSQP"),constraintType="ineq"),
    mxComputeStandardError(),
    mxComputeHessianQuality(),
    mxComputeReportDeriv()#,
    #mxComputeReportExpectation()
))

Also, consider using mxTryHardOrdinal() instead of mxTryHard(). The former is just a wrapper for the latter, but the default values of its arguments are tailored specifically for ordinal analyses.

Charlotte's picture
Offline
Joined: 07/02/2012 - 11:13
Thank you!

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
you're welcome

Glad I could be of help.

I should mention that the biggest difference between the defaults for mxTryHard() and mxTryHardOrdinal() is that argument finetuneGradient defaults to TRUE for mxTryHard() but FALSE for mxTryHardOrdinal(). Like the help page says, "finetuneGradient=FALSE is 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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
I second Joshua's suggestion

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

plan <- mxComputeSequence(steps=list(
    mxComputeGradientDescent(engine="CSOLNP"),
    mxComputeConfidenceInterval(
        plan=mxComputeGradientDescent(engine="SLSQP")),
    mxComputeStandardError(),
        mxComputeHessianQuality()
    mxComputeReportDeriv(),
        mxComputeReportExpectation()
))

Then, you would just put plan into 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.