4 posts / 0 new

Mon, 03/29/2010 - 22:02

pgseye

Offline

Joined: 10/13/2009 - 23:50

Have I done this right?

Attachment	Size
MultivariateTwinAnalysis_AL.R	10.83 KB
AxialLength.txt	12.24 KB

Hi Everyone,

I'm still rather green with all of this, but am learning a lot. I've tried to adapt a script to work with right and left eye data for the simple case of a single variable (axial length of the eye) and was wondering if someone would be kind enough to run it and comment on what I've done.

I've attached the script and the data (change .txt to .RData - it'd be nice if we could upload RData files directly).

There are 4 main queries which I've marked in the text with a ? in large lettering, so hopefully it'll be easy to find.

At the end of the twin order and zygosity constraint statements, I've added two for equating R and L data. I get a p value of 0 - so I can't assume equal means across R and L, right? I find this surprising, because if I run it as it originally was (across twin order and zygosity), I get a p value of 0.28. Maybe I've done something wrong here?
The main thrust of this is to use data from both eyes in the heritability estimation (cholesky decomposition), so here I've written algebras so that I can later set equality constraints on elements of the A and E matrices. Is this sound?

Running the models gives me a satisfactory fit with an AE model (p 0.68).

I've written a model to equate R and L additive and environmental elements and compared this to the saturated model.I get equal path coefficients, but a p value of 0.01, so again seems a poor fit (if I've done this right).

name matrix row col Estimate Std.Error
1 a11 ACE.a 1 1 0.79730983 2.121996e-314
2 a21 ACE.a 2 1 0.79730983 6.365987e-314
3 a22 ACE.a 2 2 -0.09772423 NaN
4 e11 ACE.e 1 1 0.31841590 NaN
5 e21 ACE.e 2 1 0.31841590 NaN
6 e22 ACE.e 2 2 -0.27562709 1.060998e-313
7 ACE.Mean 1 1 23.18632658 NaN
8 ACE.Mean 1 2 23.14464328 NaN
Trying to run the format output function on this, gives me errors, so I'm unsure how to convert the path coefficient to a variance component.

I'd be extremely grateful for any help. Support in my department is rather thin and so these forums are a main source of assistance for me.

Thanks in advance.

Paul

Tue, 03/30/2010 - 10:29

neale

Offline

Joined: 07/31/2009 - 15:14

Hi Paul I think your attempt

Hi Paul

I think your attempt to equate the genetic and environmental variance components between left and right is misguided. In this stretch of code:

?2

# Constrain elements of A and E matrices across R and L eyes to be the same    
    mxAlgebra( expression=a[1,1], name="AR"),
    mxAlgebra( expression=a[2,1], name="AL"),
    mxAlgebra( expression=e[1,1], name="ER"),
    mxAlgebra( expression=e[2,1], name="EL"),

you are only equating the elements in the first column of the lower triangular matrices a[] and e[]. However, element a[2,2] also contributes to the genetic variance of the second trait (on paper set up a lower triangular matrix [[a,0],[b,c]] and multiply it by its transpose. You should end up with b^2+c^2 for the genetic variance of the second trait). If you want to equate those genetic variances, it is really the results of the Cholesky multiplication that you want to equate. In the present case, since A=a %*% t(a), you want to constrain A[1,1] to equal A[2,2]. That would set the genetic variances equal; do the same for the environmental component(s) E[1,1] and E[2,2] and you will be testing whether heritability is the same for R and L.

This in hand, I think things will start to make more sense!

Thu, 04/01/2010 - 05:25

pgseye

Offline

Joined: 10/13/2009 - 23:50

Thanks for your help Mike -

Thanks for your help Mike - always appreciated.

When I try what you've suggested ie,

    mxAlgebra( expression=A[1,1], name="AR"),
    mxAlgebra( expression=A[2,2], name="AL"),
    mxAlgebra( expression=E[1,1], name="ER"),
    mxAlgebra( expression=E[2,2], name="EL"),

the model fails to converge:

> equateAEModelFit <- mxRun(equateAEModel)
Running equateAE
Warning message:
In model 'equateAE' NPSOL returned a non-zero status code 6. model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

name matrix row col Estimate Std.Error
1 a11 ACE.a 1 1 6.589639e+04 NaN
2 a21 ACE.a 2 1 6.532897e+04 NaN
3 a22 ACE.a 2 2 8.628986e+03 NaN
4 e11 ACE.e 1 1 1.465043e+02 NaN
5 e21 ACE.e 2 1 1.465037e+02 NaN
6 e22 ACE.e 2 2 -3.856362e-01 NaN
7 ACE.Mean 1 1 -7.049629e+04 NaN
8 ACE.Mean 1 2 -7.389087e+04 NaN

Did I need to tweak anything else?

Thanks

Edit: I think I figured out why I got the error. Thanks again

Thu, 04/01/2010 - 11:03

(Reply to #3) #4

neale

Offline

Joined: 07/31/2009 - 15:14

Although the optimizer is

Although the optimizer is indicating that it is not happy, it is possible that it has found the right answer anyway. Standard response to status code 6 should be to refit the model to see if the fit improves:

equateAEModelFit2 <- mxRun(equateAEModelFit)
and possibly if it improves to redo a third time etc.

Also useful in this context can be to change the starting values, although since you may have already run the un-equated AE model they may be pretty decent (and the chance is increased that the non-convergence message was a type I error saying there was a problem when there really want).

I note that your parameter estimates are extremely large. Sometimes overflow conditions can occur when variances etc. are in the thousands, and big differences in the scale of parameters can make optimization difficult because the Hessian matrix becomes ill-conditioned and difficult to invert accurately enough for the search to work properly.

Glad to see your edit, but thought these comments might help you and others with similar problems anyway...

Main menu

Navigation

You are here

Have I done this right?

?2

Main menu

User login

Navigation

You are here

Have I done this right?

?2

Search form