Twin Correlations for Scalar-Sex Limited Model?
Jeremy
Joined: 03/25/2016
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.
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"),
Log in or register to post comments
In reply to Scratch that by Jeremy
Only one phenotype?
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.Log in or register to post comments
In reply to Only one phenotype? by AdminRobK
Thank you. This was very
Thank you. This was very helpful.
Log in or register to post comments