Model including groups with both binary and continuous/ordinal data

7 replies [Last post]
ebejer's picture
Offline
Joined: 03/18/2010

Hi,

I have heard that it is now possible to run a model estimating genetic and environmental effects with twin data that includes separate variables with binary and continuous (or ordinal?) coding.

I have been running a thresholds model with age and sex as covariates (which I will attach) and I would like to include the same behavioural measures from different data sets in the one model.

Two different scales have been used to collect this information (binary and continuousish), and I want to equate the A, C and E paths across datasets with the differing scales. I have been told this will necessitate the inclusion of a group for each dataset. But I am also unsure how to do this.

Is this possible and can you help? - and forgive me if I am missing some important technical detail.

Best wishes
Jane

AttachmentSize
twinModels.R9.21 KB
neale's picture
Offline
Joined: 07/31/2009
Aha

Set E to 1 for continuous variable also - and estimate an SD parameter for the continuous variable. So you pre- and post-multiply the continuous expCov by diag matrix with sd on diagonal. Then you have the option to equate A and C parameters without nonlinear constraint.

ebejer's picture
Offline
Joined: 03/18/2010
Thanks Mike, I'll start

Thanks Mike, I'll start testing models, you may here from me again soon.
j

neale's picture
Offline
Joined: 07/31/2009
Can be done but there's likely a problem.

Jane

Since the scales are different, there is likely to be differences in mean and variance between the two. Possibly, what you want to know is whether the proportions of A C and E are the same for the two scales?

The approach you take depends on whether the two scales were obtained from the same or different pairs of twins (i.e. the same sample for both scales or different samples). In the former case, you would not increase the number of groups, but make a multivariate model. In the latter - different samples - case you would add groups. Equating the proportions of variance could be done at least two ways. The simplest would be to use algebras to compute, e.g., A/(A+C+E) for the two measures, and to add a non-linear constraint that equates the two (all three could be cbind()'ed on each side of the ==, but since they are proportions, two would be sufficient). In the case of measures on the same sample, you would set up a multivariate analysis (say with a Cholesky model) and equate the two diagonal elements of the (now 2x2) A/(A+C+E) matrix and the C/(A+C+E) matrix. Extracting the elements from the diagonal could be done with A[1,1] type notation, so you might try A[1,1]/(A[1,1]+C[1,1}+E[1,1] == A[2,2]/(A[2,2]+C[2,2]+E[2,2]. I think. If it doesn't work then, well then we'll learn something making it do so :).

HTH
Mike

ebejer's picture
Offline
Joined: 03/18/2010
I am using data from two

I am using data from two separate samples so I need to include separate calculation groups. And I do want to compare the genetic and environmental effects across the approximately same variables in the two samples. One variable represents a normal distribution (continuous) and the other is represented with a binary code.

Unfortunately I haven't been able to find a model running both continuous and ordinal /binary data including separate groups. So what I have done is put together a mishmash of script that I though might work - but I don't really understand how the matrices will combine and the use of certain code. I will attach my script.

There is a matrix used to specify the values for each a, c and e path (cholSVnv) I'm not sure why - and I don't understand the abbreviation. I'm also unsure where and how to add additional groups. I thought perhaps as two sets of mz and dz models (one for each sample). And I'm not really up to interpreting the meaning of your comments yet Mike, sorry.

Can you please give me a tip on (1) how to add a calculation group, (2) whether I need a lower cholesky for 2 vars to represent the a, c and e paths (if that's what they are), and (3) what I need to specify in the objective function to include both variables in the model and (4) the order of how the matrices combine to provide the output.

perhaps I'm asking too much?

Jane

AttachmentSize
testJointModel2.R 4.62 KB
neale's picture
Offline
Joined: 07/31/2009
Not too much but I'm a bit pressed for time just now

Jane

If you have two univariate scripts, one for the ordinal, and one for the continuous data, I would be tempted to try combining them with something like the following

jointACE <- mxModel(aceOrd, aceCont, mxAlgebraObjective(aceOrd.Objective+aceCont.Objective))

at least this would get you started. In fact I think you are on the right track with your script as is, but you need an mxAlgebraObjective to combine the objective functions from all four groups. Something like

mxAlgebra( expression=MZ1.objective+DZ1.objective+MZ2.objective+DZ2.objective, name="minus2sumloglikelihood" ),
mxAlgebraObjective("minus2sumloglikelihood")

Adding the mxConstraint would be needed to equate the heritabilities, per previous post.

ebejer's picture
Offline
Joined: 03/18/2010
Hello again, I have found a

Hello again,

I have found a model that seems to combine the two separate datasets, however it returns a code red and there are two SD's that are quite large - I have played around with the starting values and the threshold and what I have in the script appears to work best. I am wanting to get some feedback on the 'reasonableness' of the model I have specified and what I might need to change to change the code red status and reduce the standard errors. Any suggestions will be greatly appreciated (I'll attach the output and script below).

Thanks,
Jane

AttachmentSize
output.R 4.49 KB
comparingData.R 6.44 KB
smedland's picture
Offline
Joined: 08/04/2009
We have a version of the

We have a version of the script that works now. One data set has continuous data and the other has a binary measure of the same trait. To identify the ordinal model we have fixed e at 1 and estimated a and c.
The only way I can think of to 'equate' the standardised A,C and E estimates is through the mxConstraint function.
This does work and the results look sensible. But because we've used mxConstraint we lose the SE.
We can still get CI but it would be nice to have the se as well.

Can anyone think of an alturnate apporach?

thanks
Sarah (on behalf of Jane)