Multivariate binary-continuous twin model
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
?omxRunCI
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](https://search.r-project.org/CRAN/refmans/OpenMx/html/omxParallelCI.html) 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.Log in or register to post comments
In reply to ?omxRunCI by AdminHunter
To briefly elaborate
imxRobustSE()
.mxBootstrap()
and related functions for yet another way to get confidence intervals.Log in or register to post comments
In reply to To briefly elaborate by AdminRobK
Thanks both for your speedy
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
Log in or register to post comments
Odd
The
mxTryHard()
function is just a wrapper aroundmxRun()
that wiggles start values around a few times. So, you can useomxRunCI()
aftermxTryHard()
the same as you would aftermxRun()
.>>> I have tried to run my ACE model using mxRun() but have come across the following error:
I'm not sure why you'd be getting an error like that from
mxRun()
but nor frommxTryHard()
. My first guess is you aren't actually running the same model in both cases, and one of them has this error.>>> 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)?
I'll let other folks chime in with their experiences for this.
Log in or register to post comments
In reply to Odd by AdminHunter
MxRun still gives me errors
I have re-tried running this same model with mxRun but still getting the same error. I have no errors when I run the model with
mxTryHard()
. When I requestomxRunCI()
on the MxTryHard model I get the same error message:> omxRunCI(AceFit, verbose = T, optimizer = "CSOLNP") # doesnt work with mxTryhard
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 suspect that this may be due to the way my thresholds have been set up and potentially the way I arrange them in the threshnames= argument.
e.g. I set the thresholds up as a 1x2 matrix:
# Setting up thresholds - one threshold regardless of twin order and zygosity
Thresh <-mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values=c(3), lbound=c(-6), ubound=c(6),
labels=c("Tmz11", "Tmz22"), name="Th" )
Thre <-mxAlgebra( expression= Th+ (bAge%x%age)+ (bSex%x%sex) , name="expThre" )
Then when I use the threshnames argument I set these up as follows (which I assume gives me a 1x4 matrix):
# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars, thresholds="expThre", threshnames=c("BED1", "BED2", "TFEQ_Emo1", "TFEQ_Emo2" ))
objDZ <- mxExpectationNormal( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="expThre", threshnames=c("BED1", "BED2", "TFEQ_Emo1", "TFEQ_Emo2" ))
I have tried multiple ways to try and fix this e.g. trying a 1x4 matrix to set the thresholds initially but this then seems to mess up the mxAlgebra for object 'Thre' where I add covariates. Any useful pointers on how I can get round this? For context, I am now using the CSOLNP optimiser for consistency and given the now-reasonable looking CIs.
Many thanks,
Zeynep
Log in or register to post comments
In reply to MxRun still gives me errors by Zeynep_N
Threshnames order
It's difficult to debug your script without a sample dataset. I think that the threshnames should be the column names of the thresholds matrix. If you change the threshnames to only have the one ordered threshold variable, TFEQ_Emo1 and 2 I think, then maybe it will be ok. Assuming that BED1 and BED2 are not threshold variables.
threshnames=c("TFEQ_Emo1", "TFEQ_Emo2" )
However, you state you want two variables per person to be ordinal, which implies that there should be four threshnames. One thing I would want to check is that the names are in the right order in both the threshold matrix and the SelVars list. It is typical in twin studies to order variables within twin, so it goes T1var1 T1var2 T1var3 then T2var1 T2var2 T2Var3 - pretty much all ACE models in OpenMx follow this structure. So what I am thinking is
threshnames=c("BED1", "TFEQ_Emo1", "BED2", "TFEQ_Emo2" )
Also make sure that the covariates age and sex are doing the right thing. Sometimes it is useful to write the matrix algebra out by hand to ensure that everything is going in the right place. I'm presently slightly concerned that the covariates are only used on the ordinal variables (but maybe they were regressed out of the continuous variables). What is slightly concerning is that only one beta coefficient for age is being used for both variables BED and TFEQ_Emo1. It's unusual to want to constrain such regressions to be equal. If it was not intended, then note that what you really want is to end up with additions to the threshold matrices that have a different beta on age for BED and TFEQ_Emo. One thing that is easy is that these variables only have one threshold (we might be in the business of a Kronecker product of a column vector of ones times the bAge*age and bSex*sex bits, or equating all the regressions on all the thresholds within a variable). So I think some changes are needed to achieve the separate regressions on age and sex for BED and TFEQ_Emo.
betaA <-mxMatrix( type="Full", nrow=nth, ncol=2, free=T, values=-.1,
labels=c("BaTH"), name="bAge" )
betaS <-mxMatrix( type="Full", nrow=nth, ncol=2, free=T, values=.2,
labels=c("BsTH"), name="bSex" )
I think this makes bSex and bAge the same dimensions as ThMZ and ThDZ so the algebra should work. Sorry that this got a bit involved.
HTH
Log in or register to post comments