You are here

Twin Correlations for Scalar-Sex Limited Model?

4 posts / 0 new

Log in or register to post comments

Last post

Fri, 03/25/2016 - 17:41

#1

Jeremy's picture

Jeremy

Offline

Joined: 03/25/2016 - 13:07

Twin Correlations for Scalar-Sex Limited Model?

New poster here, and new-ish user to twin analyses. I am trying to run an ACE model that I believe is called "scalar sex limited" -- I want to allow the total variance to differ across sex, but I want to constrain the A, C, and E to be the same so that I am reporting just one heritability estimate. Here are the key parts of my syntax:

twinACEModel <- mxModel("twinACE",
mxModel("ACE",
# ace path coefficients for males and females
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="am11", name="am"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="cm11", name="cm"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.5, label="em11", name="em"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="am11", name="af"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.3, label="cm11", name="cf"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.5, label="em11", name="ef"),
mxMatrix(type="Diag", 2, 2, labels=rep("k",2), free=TRUE, values=1, name="scalar"),
 
# ACE variance components for males and females
mxAlgebra(am %*% t(am), name="Am"),
mxAlgebra(cm %*% t(cm), name="Cm"),
mxAlgebra(em %*% t(em), name="Em"),
mxAlgebra(af %*% t(af), name="Af"),
mxAlgebra(cf %*% t(cf), name="Cf"),
mxAlgebra(ef %*% t(ef), name="Ef"),

I'm correcting for age and sex.

mxModel("DZm",
mxData(observed=dzData , type="raw"),
 
 mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age"  ), 
 mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ),  
 
mxAlgebra( expression= ACE.expMean + ACE.betasex %x% obs_sex + ACE.betaage %x% obs_age, name="expMean"),
 
mxFIMLObjective(covariance="ACE.expCovDZm", means="ACE.expMean", dimnames=selVars)
),

I'm trying to figure out how to get the twin correlations. I've seen some other posts where people get the twin correlations through saturated models -- is there something wrong with trying to grab them from my ACE model? E.g., just asking for confidence intervals for expCovDZm, expCovDZf, expCovMZm, and expCovMZf? (If my model has worked right, I believe expCovDZm should = expCovDZf, and same for MZ).

Thanks.

Fri, 03/25/2016 - 20:48

#2

Jeremy's picture

Jeremy

Offline

Joined: 03/25/2016 - 13:07

Scratch that

I'm realizing that expCovDZm will not equal expCovDZf since this isn't standardized, right?
Do I need to add something to the syntax I've pasted at the bottom like this:

mxAlgebra((Am+Cm)/Vm, name="CorMZm"),
mxAlgebra(((.5*Am)+Cm)/Vm, name="CorDZm"),
mxAlgebra((Af+Cf)/Vf, name="CorMZf"),
mxAlgebra(((.5*Af)+Cf)/Vf, name="CorDZf"),

And then CorMZm will equal CorMZf, and CorDZm will equal CorDZf, and this is an acceptable way to calculate the twin correlations? Is there another way I should be doing it?

twinACEModel <- mxModel("twinACE",
mxModel("ACE",
# ace path coefficients for males and females
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="am11", name="am"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="cm11", name="cm"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="em11", name="em"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="am11", name="af"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="cm11", name="cf"),
mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="em11", name="ef"),
mxMatrix(type="Diag", 2, 2, labels=rep("k",2), free=TRUE, values=1, name="scalar"),
 
# ACE variance components for males and females
mxAlgebra(am %*% t(am), name="Am"),
mxAlgebra(cm %*% t(cm), name="Cm"),
mxAlgebra(em %*% t(em), name="Em"),
mxAlgebra(af %*% t(af), name="Af"),
mxAlgebra(cf %*% t(cf), name="Cf"),
mxAlgebra(ef %*% t(ef), name="Ef"),
 
# Declare a matrix for the definition variable regression parameters, called beta
    mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("bage_caster"), name="betaage"),
    mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("bsex_caster"), name="betasex"),
 
# Total variation
mxAlgebra(Am+Cm+Em, name="Vm"),
mxAlgebra(Af+Cf+Ef, name="Vf"),
 
# Standard deviations for males and females
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*Vm)), name="isdm"),
mxAlgebra( expression=solve(sqrt(I*Vf)), name="isdf"),
 
# Standardized path estimates
mxAlgebra(am %*% isdm, name="stam"),
mxAlgebra(cm %*% isdm, name="stcm"),
mxAlgebra(em %*% isdm, name="stem"),
mxAlgebra(af %*% isdf, name="staf"),
mxAlgebra(cf %*% isdf, name="stcf"),
mxAlgebra(ef %*% isdf, name="stef"),
 
# Standardized variance components
mxAlgebra(Am/Vm, name="h2m"),
mxAlgebra(Cm/Vm, name="c2m"),
mxAlgebra(Em/Vm, name="e2m"),
mxAlgebra(Af/Vf, name="h2f"),
mxAlgebra(Cf/Vf, name="c2f"),
mxAlgebra(Ef/Vf, name="e2f"),

Wed, 03/30/2016 - 17:48

(Reply to #2) #3

AdminRobK's picture

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

Only one phenotype?

I'm realizing that expCovDZm will not equal expCovDZf since this isn't standardized, right?
Do I need to add something to the syntax I've pasted at the bottom like this:
mxAlgebra((Am+Cm)/Vm, name="CorMZm"),
mxAlgebra(((.5*Am)+Cm)/Vm, name="CorDZm"),
mxAlgebra((Af+Cf)/Vf, name="CorMZf"),
mxAlgebra(((.5*Af)+Cf)/Vf, name="CorDZf"),

I'm assuming from your syntax that you are only working with one phenotype. If that's the case, then yes, those four algebras will give you the correlations you want. They would be the model-expected correlations from the ACE model, though. You may want/need to report correlations from a saturated model instead, or in addition. For the saturated model, cov2cor() can be useful for turning a covariance matrix into a correlation matrix.

Sat, 04/02/2016 - 13:41

(Reply to #3) #4

Jeremy's picture

Jeremy

Offline

Joined: 03/25/2016 - 13:07

Thank you. This was very

Thank you. This was very helpful.