| Attachment | Size |
|---|---|
 MVModel_ACE_OpenMx_ZN.R | 44.23 KB |
Hi,
I am currently working on a multivariate correlated factors twin model with combined binary and continuous traits. I have six variables in total, the first two being binary data and the remaining four being continuous. I have been struggling for a while now in extracting meaningful output, specifically for the aetiological correlations. I receive genetic correlations with corresponding confidence interval lower bounds that are higher than the estimate itself (even though status is ok and model has converged).
For example:
confidence intervals:
lbound estimate ubound note
MZ.Acor[3,2] 0.4411063 0.4411063 0.999999 !!!
CI details:
parameter side value fit diagnostic statusCode method a11
1 MZ.Acor[3,2] lower 0.9285487 11548.08 success OK neale-miller-1997 0.6466693
2 MZ.Acor[3,2] upper 0.9999990 11548.08 success OK neale-miller-1997 0.6227424
This is the output to suggest the model has converged:
> AceFit$output$status$code
[1] 0
In case helpful, I am using the SLSQP optimiser and mxversion is 2.20.6
I have searched the openmx forums for different ways to adjust my script and troubleshoot (including adjusting 'function precision', mxTryHard/MxTryHardOrdinal with extratries) but have had no luck.
I am attaching my script here, though I cannot share the data. I appreciate this is quite a big model, so any tips and keen eyes on my script will be hugely appreciated!
Thanks in advance,
Zeynep
Printer-friendly version
I'd recommend separating your thinking about the model convergence from the confidence interval convergence. The
mxTryHard()/mxTryHardOrdinal()functions are best suited to getting your model to converge on the best parameter estimates. This is literally a separate optimization step from getting the best confidence intervals.The function
omxRunCI()is designed just for running the confidence interval optimization. There's more information on it here at the help page. The big thing to try is switching optimizers for the CIs. The best optimizer for your parameter estimates might not be the best for your CIs. The different optimizers have rather different strategies for getting CIs.If that all fails and you are willing to dance with the devil, then you could try standard-error-based Wald-type confidence intervals (gasp!). Standard errors for arbitrary algebras are available from
mxSE()and Wald-type confidence intervals fromconfint(). I think the latter will only give CIs for the parameters, not algebras.To briefly elaborate on Hunter's post:
imxRobustSE().mxBootstrap()and related functions for yet another way to get confidence intervals.Thanks both for your speedy response, this is extremely helpful.
I have re-run my model with the 'CSOLNP' optimiser (without any other adjustments, still using mxTryHard) and am now receiving more meaningful output:
confidence intervals:
lbound estimate ubound note
MZ.Acor[3,2] 0.04128288 0.4411438 1
CI details:
parameter side value fit diagnostic statusCode method
1 MZ.Acor[3,2] lower 0.04128288 11548.07 success OK neale-miller-1997
2 MZ.Acor[3,2] upper 1.00000000 11548.04 success OK neale-miller-1997
I am, however, curious to use omxRunCI. As far as I understand, the function omxRunCI() only works when a model has been run through mxRun(). I have tried to run my ACE model using mxRun() but have come across the following error:
Error: The thresholds matrix associated with the expectation function in model 'MZ' is not of the same length as the 'threshnames' argument provided by the expectation function. The 'threshnames' argument is of length 4 and the expected covariance matrix has 2 columns.
I think it comes down to the way my thresholds have been aligned, but I am unsure of how to fix this. Any suggestions?
Once I get past this, I can run the model using mxRun and request CIs seperately via omxRunCI(). I am thinking something along the lines of:
ACEModel <-mxModel( ACEModel, mxCI (c ('Acor[3,2]') ))
ACEfit <- mxRun(ACEModel, intervals=F)
omxRunCI(AceFit, verbose = T, optimizer = "CSOLNP") #optimiser TBC
On a related note, is there a particular optimiser that would work best in this type of situation based on your experience (e.g. CSOLNP/ NPSOL)?
Thank you again for your time!
Zeynep