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

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.

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

I don't think so but will scout around.

That would be very helpful, thanks! =)

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

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

... 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

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.

Hi Mike,

Thank you! I'll check the reversed order!

So you don't see a problem with the NANs (see above)?

Cheers

Charlotte