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
Printer-friendly version
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 ofh2are the narrow-sense heritabilities of the three traits. The off-diagonals ofh2can 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.