Multivariate binary-continuous twin model

7 posts / 0 new

Fri, 03/17/2023 - 14:06

Zeynep_N

Offline

Joined: 10/03/2019 - 12:58

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

Fri, 03/17/2023 - 17:51

AdminHunter

Offline

Joined: 03/01/2013 - 11:03

?omxRunCI

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 from confint(). I think the latter will only give CIs for the parameters, not algebras.

Fri, 03/17/2023 - 20:53

(Reply to #2) #3

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

To briefly elaborate

To briefly elaborate on Hunter's post:

You can get "better" standard errors, and therefore, "better" standard-error-based Wald-type confidence intervals, if you pass your MxModel object through imxRobustSE().
Look into mxBootstrap() and related functions for yet another way to get confidence intervals.

Sun, 03/19/2023 - 14:50

(Reply to #3) #4

Zeynep_N

Offline

Joined: 10/03/2019 - 12:58

Thanks both for your speedy

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

Mon, 04/03/2023 - 12:51

AdminHunter

Offline

Joined: 03/01/2013 - 11:03

Odd

>>> As far as I understand, the function omxRunCI() only works when a model has been run through mxRun().

The mxTryHard() function is just a wrapper around mxRun() that wiggles start values around a few times. So, you can use omxRunCI() after mxTryHard() the same as you would after mxRun().

>>> 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 from mxTryHard(). 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.

Tue, 04/04/2023 - 09:13

(Reply to #5) #6

Zeynep_N

Offline

Joined: 10/03/2019 - 12:58

MxRun still gives me errors

Thanks so much for your helpful reply, I really appreciate your time!

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 request omxRunCI() 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

Fri, 04/14/2023 - 11:01

(Reply to #6) #7

AdminNeale

Offline

Joined: 03/01/2013 - 14:09

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 bAgeage and bSexsex 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

Main menu

Navigation

You are here

Multivariate binary-continuous twin model

Main menu

User login

Navigation

You are here

Multivariate binary-continuous twin model

Search form