Hi,

I am running a trivariate ACE Cholesky model where trait 1 and 2 are positively correlated, but trait 3 is negatively correlated with 1 and 2. The problem is that when I run the model, the upperbound CI's for the phenotypic correlation due to a, c and e is over 1 for some estimates (e.g. MZM.h2[1,3] ubound= 1.14034072)

I have 3 other similar analyses with different traits, with similar correlations and the same problem.

Any ideas how to fix this?

Thank you very much!

# -------------------------------------------------------------------------------------------------

# Specify Cholesky ACE Decomposition

# -------------------------------------------------------------------------------------------------

# To create Labels for Lower Triangular Matrices

aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

cLabs <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

eLabs <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

# Matrices to store a, c, and e Path Coefficients

pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=4, labels=aLabs,name="a" )

pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,F,F,T,F,T), values=c(1,0,0,1,0,1),labels=cLabs, lbound=0, ubound=NA, name="c" )

pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=4, labels=eLabs, name="e" )

# Algebra for expected Means

MeanM <-mxMatrix(type="Full", nrow=1, ncol=ntv, free=T, values=1 , labels=c("Mm11", "Mm21", "Mm31","Mm11", "Mm21", "Mm31"), name="ExpMeanM")

MeanF <-mxMatrix(type="Full", nrow=1, ncol=ntv, free=T, values=1, labels=c("Fm11", "Fm21","Fm31", "Fm11", "Fm21","Fm31"), name="ExpMeanF")

MeanOS <-mxMatrix(type="Full", nrow=1, ncol=ntv, free=T, values=1, labels=c("Mm11", "Mm21", "Mm31","Fm11", "Fm21","Fm31"), name="ExpMeanOS")

# Matrices generated to hold A, C, and E computed Variance Components

covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )

covE <- mxAlgebra( expression=e %*% t(e), name="E" )

# Algebra to compute standardized variance components

covP <- mxAlgebra( expression=A+C+E, name="V" )

StA <- mxAlgebra( expression=A/V, name="h2")

StC <- mxAlgebra( expression=C/V, name="c2")

StE <- mxAlgebra( expression=E/V, name="e2")

# Algebra to compute Phenotypic, A, C & E correlations

matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")

rph <- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")

rA <- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )

rC <- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )

rE <- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )

covMZM <- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),

cbind(A+C , A+C+E)), name="ExpCovMZM" )

covDZM <- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),

cbind(0.5%x%A+C , A+C+E)), name="ExpCovDZM" )

covMZF <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),

cbind(A+C , A+C+E)), name="ExpCovMZF" )

covDZF <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C),

cbind(0.5%x%A+C , A+C+E)), name="ExpCovDZF" )

covDZOS <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C),

cbind(0.5%x%A+C , A+C+E)), name="ExpCovDZOS" )

# Data objects for Multiple Groups

dataMZM <-mxData(observed=mzmData, type="raw")

dataDZM <-mxData(observed=dzmData, type="raw")

dataMZF <-mxData(observed=mzfData, type="raw")

dataDZF <-mxData(observed=dzfData, type="raw")

dataDZOS <-mxData(observed=dzosData, type="raw")

# Objective objects for Multiple Groups

objMZM <- mxExpectationNormal(covariance="ExpCovMZM", means="ExpMeanM", dimnames=selVars)

objDZM <- mxExpectationNormal(covariance="ExpCovDZM", means="ExpMeanM", dimnames=selVars)

objMZF <- mxExpectationNormal(covariance="ExpCovMZF", means="ExpMeanF", dimnames=selVars)

objDZF <- mxExpectationNormal(covariance="ExpCovDZF", means="ExpMeanF", dimnames=selVars)

objDZOS <- mxExpectationNormal(covariance="ExpCovDZOS", means="ExpMeanOS", dimnames=selVars)

fitFunction <- mxFitFunctionML()

# Combine Groups

pars <-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE,matI, rph, rA, rC, rE,fitFunction)

modelMZM <-mxModel(pars, MeanM, covMZM, dataMZM, objMZM, name="MZM")

modelDZM <-mxModel(pars, MeanM, covDZM, dataDZM, objDZM, name="DZM")

modelMZF <-mxModel(pars, MeanF, covMZF, dataMZF, objMZF, name="MZF")

modelDZF <-mxModel(pars, MeanF, covDZF, dataDZF, objDZF, name="DZF")

modelDZOS <-mxModel(pars, MeanM, MeanF, MeanOS, covDZOS, dataDZOS, objDZOS,name="DZOS")

minus2ll<- mxAlgebra( expression=MZM.objective + DZM.objective + MZF.objective + DZF.objective + DZOS.objective, name="m2LL" )

obj <- mxFitFunctionAlgebra( "m2LL" )

ci1 <- mxCI (c ('MZM.h2[1,1]','MZM.h2[2,2]','MZM.h2[3,3]','MZM.c2[1,1]','MZM.c2[2,2]','MZM.c2[3,3]',

'MZM.e2[1,1]','MZM.e2[2,2]','MZM.e2[3,3]')) #a2, c2, e2

ci2 <- mxCI (c ('MZM.Rph[1,2]', 'MZM.Rph[1,3]','MZM.Rph[2,3]')) # rPh

ci3 <- mxCI (c ('MZM.Ra[1,2]', 'MZM.Ra[1,3]','MZM.Ra[3,2]',

'MZM.Rc[1,2]','MZM.Rc[1,3]','MZM.Rc[3,2]',

'MZM.Re[1,2]','MZM.Re[1,3]','MZM.Re[3,2]')) #rA, rC, rE

ci5 <-mxCI (c('MZM.h2[1,2]','MZM.h2[1,3]','MZM.h2[2,3]','MZM.c2[1,2]','MZM.c2[1,3]','MZM.c2[2,3]',

'MZM.e2[1,2]','MZM.e2[1,3]','MZM.e2[2,3]')) # rph due to A,C,E (%)

CholAceModel<- mxModel( "CholACE", pars, modelMZM, modelDZM, modelMZF, modelDZF, modelDZOS, minus2ll, obj,ci1,ci2,ci3,ci5)

# ---------------------------------------------------------------------------------------------------------------------

# RUN Cholesky Decomposition ACE Model

CholAceFit <- mxRun(CholAceModel, intervals=T)

(CholAceSum <- summary(CholAceFit))

OUT PUT:

Summary of CholACE

free parameters:

name matrix row col Estimate Std.Error A lbound ubound

1 a11 a 1 1 4.205001e+00 0.4512848

2 a21 a 2 1 1.912439e+00 0.3111661

3 a31 a 3 1 -1.294262e+00 0.2762099

4 a22 a 2 2 1.370645e+00 0.4425790

5 a32 a 3 2 -9.997835e-01 0.4921637

6 a33 a 3 3 -1.409512e+00 0.8254633

7 c11 c 1 1 1.015753e+00 1.3711636 0

8 c22 c 2 2 -4.529173e-13 1.0553574 0!

9 c33 c 3 3 1.000354e+00 0.8782109 ! 0

10 e11 e 1 1 5.105396e+00 0.1506051

11 e21 e 2 1 5.557607e-01 0.1884338

12 e31 e 3 1 -1.812631e-01 0.1744390

13 e22 e 2 2 3.473323e+00 0.1324603

14 e32 e 3 2 -6.184976e-01 0.1714142

15 e33 e 3 3 3.132411e+00 0.1334209

16 Mm11 MZM.ExpMeanM 1 rEOE1 -1.706490e-02 0.2070563

17 Mm21 MZM.ExpMeanM 1 rNeuro1 -2.181607e-02 0.2019051

18 Mm31 MZM.ExpMeanM 1 rExtra1 -4.209456e-02 0.2002149

19 Fm11 MZF.ExpMeanF 1 rEOE1 3.887317e-02 0.1791824

20 Fm21 MZF.ExpMeanF 1 rNeuro1 -1.335770e-03 0.1586503

21 Fm31 MZF.ExpMeanF 1 rExtra1 2.600313e-02 0.1569875

confidence intervals:

lbound estimate ubound note

MZM.h2[1,1] 2.237774e-01 3.948746e-01 0.48066135

MZM.h2[2,2] 1.655729e-01 3.091246e-01 0.40500087

MZM.h2[3,3] 9.902937e-02 2.933637e-01 0.45343557

MZM.c2[1,1] 1.622543e-21 2.304111e-02 0.14867725

MZM.c2[2,2] 1.145429e-26 1.145429e-26 0.10657173 !!!

MZM.c2[3,3] 0.000000e+00 6.297918e-02 0.22569071

MZM.e2[1,1] 5.193389e-01 5.820843e-01 0.65408871

MZM.e2[2,2] 5.949989e-01 6.908754e-01 0.79387528

MZM.e2[3,3] 5.465646e-01 6.436571e-01 0.75100534

MZM.Rph[1,2] 3.335501e-01 3.841715e-01 0.43199250

MZM.Rph[1,3] -2.935708e-01 -2.387245e-01 -0.18176681

MZM.Rph[2,3] -4.114417e-01 -3.612849e-01 -0.30900866

MZM.Ra[1,2] 6.118084e-01 8.128044e-01 1.00000000

MZM.Ra[1,3] -1.000000e+00 -5.994647e-01 -0.34234354

MZM.Ra[3,2] -1.000000e+00 -7.570034e-01 -0.48455014

MZM.Rc[1,2] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.Rc[1,3] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.Rc[3,2] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.Re[1,2] 5.435395e-02 1.579986e-01 0.25785609

MZM.Re[1,3] -1.626099e-01 -5.667960e-02 0.05022950

MZM.Re[3,2] -3.019891e-01 -1.999256e-01 -0.09445959

MZM.h2[1,2] 5.677666e-01 7.391920e-01 0.90848317

MZM.h2[1,3] 5.834506e-01 8.546718e-01 1.14034072

MZM.h2[2,3] 4.318567e-01 6.309834e-01 0.82191707

MZM.c2[1,2] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.c2[1,3] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.c2[2,3] 0.000000e+00 0.000000e+00 0.00000000 !!!

MZM.e2[1,2] 9.151570e-02 2.608080e-01 0.43223321

MZM.e2[1,3] -1.403418e-01 1.453282e-01 0.41655169

MZM.e2[2,3] 1.780829e-01 3.690166e-01 0.56814336

observed statistics: 5178

estimated parameters: 21

degrees of freedom: 5157

fit value ( -2lnL units ): 31565.78

number of observations: 1446

Information Criteria:

| df Penalty | Parameters Penalty | Sample-Size Adjusted

AIC: 21251.785 31607.78 NA

BIC: -5959.417 31718.59 31651.88

Your MxAlgebra with Mxname "h2" is defined as

`A/V`

, which is the additive-genetic covariance matrix elementwise divided by the total phenotypic covariance matrix. The diagonals of`h2`

are the narrow-sense heritabilities of the three traits. The off-diagonals of`h2`

can be greater than 1, as they're neither correlations nor proportions. I don't know of any meaningful interpretation of those off-diagonals, other than the rather literal "ratio of two traits' genetic covariance to their total phenotypic covariance." Edit: Actually, they'd be interpretable as bivariate heritabilities under other circumstances; see Prof. Neale's post below.There seems to be some confusion in the OpenMx user base about this. See two recent threads, http://openmx.psyc.virginia.edu/thread/3969 and http://openmx.psyc.virginia.edu/thread/4134 .

Thank you Rob. Re " I don't know of any meaningful interpretation of those off-diagonals, other than the rather literal "ratio of two traits' genetic covariance to their total phenotypic covariance." My understanding is that those off-diagonals can be interpreted as the proportion of phenotypic correlations between two traits due to their correlating genetic/shared/nonshared environments. These add up to 1(which are usually presented as percentage), and I have never seen CIs for this statistics to be over 1 or under -1, in the published articles where they report this stat. So I don't understand why my analysis is showing a CI over 1. Am I doing something wrong? Any ideas how I can fix this so that the CI is not over 1 or under -1. Maybe this happens to other people, but no one talks about it?

OK, now that I think about it, that interpretation may be correct. Nonetheless, the values of the off-diagonals are not bounded to the interval [0,1]. Suppose we have two phenotypes with a genetic covariance of -5.0, a shared-environmental covariance of 1.5, and a non-shared environmental covariance of 1.0. The phenotypic covariance will then be -5 + 1.5 + 1 = -2.5. The off-diagonal for additive genetics would then be -5/-2.5 = 2.0; for the shared environment, 1.5/-2.5 = -0.6; for the non-shared environment, 1/-2.5 = -0.4. Nothing about this scenario is implausible, and yet not a single one of the off-diagonals falls within the interval [0,1]. In your case, two of the three phenotypic correlations are negative; consider the substantive interpretability of a proportion of a whole, when the "whole" can be negative.

If the phenotypic correlation were positive, things might "look saner," but there's still no guarantee that all those off-diagonals will fall in [0,1]. Suppose we have two phenotypes with a genetic covariance of 5.0, a shared-environmental covariance of 1.5, and a non-shared environmental covariance of -1.0. The phenotypic covariance will then be 5 + 1.5 - 1 = 5.5. The off-diagonal for additive genetics would then be 5/5.5 = approx 0.909; for the shared environment, 1.5/5.5 = approx 0.272; for the non-shared environment, -1/5.5 = approx -0.182.

I think it's the negative correlations that are really throwing you off. It might be possible to force all of your phenotypes to be positively correlated if you reflect ("reverse-score") trait 3. Then, in the resulting paper/talk, make sure your reader/audience understands that trait 3 is reflected, e.g. if trait 3 is extraversion, then when reflected, it would be introversion.

The background of proportion is really in a ratio, the proportion of two positive quantities. When the quantities can contain mixtures of positive and negative components, the proportion (sometimes known as "bivariate heritability" - another problematic quantity when signs are mixed) becomes meaningless. In these circumstances I prefer to describe the contribution of each component to the phenotypic correlation. So a phenotypic correlation might be made up of .2 due to A, .2 due to C and -.3 due to E, giving a phenotypic correlation of .1

Thank you both for the responses. It makes more sense now.