Freely estimate genetic correlation DOS twins

9 replies [Last post]
Charlotte's picture
Offline
Joined: 07/02/2012

Hi,

I'm adapting the BIVARIATE Boulder-skript atm and I'm trying to freely estimate the paths connecting the As of the two phenotypes for the DOS-twins. So it should be 0.5 for the DZ twins, but freely estimated for the DOS twins.

With this skrip, I get ONE estimate for both paths:

pathRados <- mxMatrix (type="Full", nrow=1, ncol=1, free=TRUE, values=0.5, label="rados", lbound=-1, ubound=1, name="ra")

expCovDOS <- mxAlgebra(name = "expCovDOS",
expression = rbind (cbind(Am+Cm+Em, (ra%x%(am%*%t(af)))+cm%*%t(cf)),
cbind((ra%x%(af%*%t(am)))+cf%*%t(cm), Af+Cf+Ef))

But I'd like to have two estimates:

Estimate 1: The path connecting the A of phenotype 1, first twin, with the A of phenotype 1, second twin.
Estimate 2: The path connecting the A of phenotype 2, first twin, with the A of phenotype 2, second twin.

I couldn't find an elegant solution yet - do you have any ideas or perhaps even a working skript?

Any help would be very much appreciated!

Regards
Charlotte

neale's picture
Offline
Joined: 07/31/2009
Diagonal matrix?

Charlotte

I think you want a diagonal rados matrix with the ra's on the diagonal. I would recommend:

am%*% rados %*% af

Beware of a couple of gotchas. First, probably not sensible to have the elements of rados outside the range of 0 to .5 (.5 is what it would be for DZ twins of same sex, i.e. same genetic factors in both sexes). Second, make sure am & af parameters are bounded to be positive. Otherwise, the algorithm may 'cheat' by making say am positive and af negative, thereby estimating a negative rados correlation. And third, don't use the Cholesky decomposition, it does weird things like change the fit if you change the order of the variables. See http://www.vipbg.vcu.edu/vipbg/Articles/twinres-2006-sexlimgxe.pdf for why.

Charlotte's picture
Offline
Joined: 07/02/2012
Thank you so much for the

Thank you so much for the fast response, this is very helpful!

In your article, you suggest a way to change the Cholesky model in order to make it applicable to the "nonscalar sex limitation"-case (p.486). Do you have - by any chance - a working openMx script on this?

Regards
Charlotte

neale's picture
Offline
Joined: 07/31/2009
Sorry

I don't think so but will scout around.

Charlotte's picture
Offline
Joined: 07/02/2012
Hi Michael,I've changed the

Hi Michael,

I've changed the skript according to the correlational approach (scalar & non-scalar sex limitations) that you've suggested. I'm not sure whether it is correct now. My ACE-estimates are quite different from the ones I got with the Cholesky approach (e.g. A going from 9% to 22%) - is that possible?

##### A (here males & mf, in my skript I also have females & fm, not included here)
mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.5,0.5), labels = c("am1","am2"), name = "am")

mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,0.5,0.5,1), labels=c("fixed","ram","ram","fixed"), lbound=-1, ubound=1, name="rAM")

mxAlgebra(expression = am %*% rAM %*% am, name = "Am")

and for DOSmf:
mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=TRUE, values=c(0.5,0.5,0.5,0.5), labels=c("rAmAf","rSmAf","rAmSf","rSmSf"),lbound=-1, ubound=1, name="rAmf")

##### C (same for males and females)
mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.5,0.5), labels = c("cm1","cm2"), name = "cm")

mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,0.5,0.5,1), labels=c("fixed","rc","rc","fixed"), lbound=-1, ubound=1, name="rC")

mxAlgebra(expression = cm %*% rC %*% cm, name = "Cm")

#### E
mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.5,0.5), labels = c("em1","em2"), name = "em")

mxAlgebra(expression = em %*% em, name = "Em")

I'm not sure about how to estimate the covariance of A. I don't know how to take into account that it is different for DOS twins?
I guess that this is wrong:
mxAlgebra(expression = Am+Cm+Em, name = "Vm" [=total variance])? (given that mxAlgebra(expression = am %*% rAM %*% am, name = "Am"), thus not including the DOS-matrices)
And the same problem hold for the variance as it has additional parameters (rSmAf and rAmSf) for the DOS twins?

Also, I am a bit confused as I thought that "everyone" uses Cholesky for bivariate models even if there are quantitative sex differences? (although I'm not very familiar with bivariate models yet) But you are right, changing the order of the variables does change ACE estimates. I would like to know "how wrong" it is to use Cholesky in this case - how different are the estimates from the actual estimates?

Regards
Charlotte

Charlotte's picture
Offline
Joined: 07/02/2012
NANs

Hi Mike,

The model is running now and I get ACE-estimates that are very close to what I would expect based on the correlations.
I ran this bivariate model on a couple of phenotypes. I am a bit confused as for some phenotypes (not for all), I get NaNs for the standard errors. For those estimates that are related to the C-component I could imaging that it’s because C had almost no impact on the phenotype. However, for the others (I think it’s always the genetic correlations), I don’t understand why the program has trouble calculating the std. errors. Is something going wrong here? You told me that NaN could be a sign of underidentification? A part of the script and the estimates for one phenotype are attached.

Also, I was trying to test for the significance of quantitative and qualitative sex differences. I wanted to test for QUANTITATIVE differences as follows (labels taken from figure 6 of your paper):
a1m = a1f (& the same for C and E)
a2m = a2f (& the same for C and E)
rgm=rgf (& the same for C and E)
And I wanted to keep the four DOS correlations free (rgm1f1, rgm1f2, rgm2f1, rgm2f2) – is that right?

For the qualitative sex differences, I wanted to equate this:
rgm1f1 = 0.5
rgm2f2 = 0.5
but then I don’t know what to do with the other two (rgm1f2, rgm2f1). I guess that these can only be equated to 0.5*rg in the absence of quantitative sex differences?

I would be very grateful for any advice!

Regards
Charlotte

# Matrices a, c, and e to store a, c, and e Path Coefficients
pathAm <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.4,1.2), labels = c("am1","am2"), name = "am", lbound=0)
pathCm <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.2,0.3), labels = c("cm1","cm2"), name = "cm", lbound=0)
pathEm <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.8,1.1), labels = c("em1","em2"), name = "em", lbound=0)

pathAf <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.4,0.9), labels = c("af1","af2"), name = "af", lbound=0)
pathCf <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.3,0.3), labels = c("cf1","cf2"), name = "cf", lbound=0)
pathEf <- mxMatrix(type = "Diag", nrow = nphen, ncol = nphen, free=TRUE, values=c(0.7,1.1), labels = c("ef1","ef2"), name = "ef", lbound=0)

pathRAM <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,-0.6,-0.6,1), labels=c("fix","ram","ram","fix"), lbound=-1, ubound=1, name="rAM")
pathRAF <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,-0.4,-0.4,1), labels=c("fix","raf","raf","fix"), lbound=-1, ubound=1, name="rAF")
pathRAmf <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=TRUE, values=c(0.2,-0.2,0.01,0.2), labels=c("rAmAf","rSmAf","rAmSf","rSmSf"),lbound=c(0,-1,-1,0), ubound=1, name="rAmf")
pathRAfm <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=TRUE, values=c(0.2,0.01,-0.2,0.2), labels=c("rAmAf","rAmSf","rSmAf","rSmSf"),lbound=c(0,-1,-1,0), ubound=1, name="rAfm")

pathRC <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,-1,-1,1), labels=c("fix","rc","rc","fix"), lbound=-1, ubound=1, name="rC")
pathRE <- mxMatrix (type="Full", nrow= nphen, ncol= nphen, free=c(FALSE,TRUE,TRUE,FALSE), values=c(1,-0.2,-0.2,1), labels=c("fix","re","re","fix"), lbound=-1, ubound=1, name="rE")

free parameters:
name matrix row col Estimate Std.Error
1 am1 MZMmodel.am 1 1 1.719317e-01 0.131521455
2 am2 MZMmodel.am 2 2 1.062506e+00 0.067926413
3 cm1 MZMmodel.cm 1 1 1.321520e-01 0.126093474
4 cm2 MZMmodel.cm 2 2 1.161182e-01 0.347575873
5 em1 MZMmodel.em 1 1 4.568589e-01 0.020430192
6 em2 MZMmodel.em 2 2 1.084936e+00 0.046684280
7 af1 MZMmodel.af 1 1 2.413808e-01 0.015959201
8 af2 MZMmodel.af 2 2 8.656119e-01 0.094059102
9 cf1 MZMmodel.cf 1 1 -5.793705e-19 NaN
10 cf2 MZMmodel.cf 2 2 2.592262e-01 0.245674655
11 ef1 MZMmodel.ef 1 1 3.920726e-01 0.009685523
12 ef2 MZMmodel.ef 2 2 1.044822e+00 0.029209402
13 ram MZMmodel.rAM 2 1 7.965951e-01 0.525799393
14 raf MZMmodel.rAF 2 1 6.214928e-01 NaN
15 rAmAf MZMmodel.rAmf 1 1 5.795599e-01 0.205931625
16 rSmAf MZMmodel.rAmf 2 1 2.418333e-01 0.092369742
17 rAmSf MZMmodel.rAmf 1 2 -1.372497e-01 NaN
18 rSmSf MZMmodel.rAmf 2 2 2.937841e-01 0.116614314
19 rc MZMmodel.rC 2 1 1.000000e+00 NaN
20 re MZMmodel.rE 2 1 1.798579e-01 0.030778450
21 mean1m MZMmodel.Mm 1 1 2.185487e+00 0.033377761
22 mean2m MZMmodel.Mm 1 2 1.911841e+00 0.104139629
23 mean1f MZMmodel.Mf 1 1 2.218440e+00 0.032525170
24 mean2f MZMmodel.Mf 1 2 1.720874e+00 0.101020062
25 batt MZMmodel.b 1 1 -1.942780e-03 0.001023805
26 bsp MZMmodel.b 1 2 -1.306742e-02 0.003183240

Charlotte's picture
Offline
Joined: 07/02/2012
... I'm getting slightly

... I'm getting slightly different CIs for the [1,2] versus [2,1], altough this matrix should be symmetrical (and it is for the point estimates).
Is this related to the problem above?

MZMmodel.emCI[1,1] 6.907786e-01 0.816132558 0.9171619
MZMmodel.emCI[1,2] 2.158782e-01 0.356576871 0.5053817
MZMmodel.emCI[2,1] 2.297038e-01 0.356576871 0.5053817
MZMmodel.emCI[2,2] 4.252797e-01 0.507476963 0.6050545

Cheers
Charlotte

neale's picture
Offline
Joined: 07/31/2009
Order

That's a bit weird. I can only think it is a numerical precision thing related to the previous state of the system when the CI is estimated. If you ask for them in reverse order, are the two lower CI's reversed?

In the end, CI's are only approximate so perhaps you are bumping up against the limited accuracy.

Charlotte's picture
Offline
Joined: 07/02/2012
Hi Mike, Thank you! I'll

Hi Mike,

Thank you! I'll check the reversed order!
So you don't see a problem with the NANs (see above)?

Cheers
Charlotte

Charlotte's picture
Offline
Joined: 07/02/2012
That would be very helpful,

That would be very helpful, thanks! =)