Login

News

Multivariate model with positive and negative correlations

Posted on
No user picture. Ell Joined: 11/26/2015
Forums
Behavioral Genetics Models
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

Replied on Fri, 05/13/2016 - 15:48
Picture of user. AdminRobK Joined: 01/24/2014

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 .

Log in or register to post comments
Replied on Wed, 05/18/2016 - 10:23
No user picture. Ell Joined: 11/26/2015

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?
Log in or register to post comments
Replied on Wed, 05/18/2016 - 11:24
Picture of user. AdminRobK Joined: 01/24/2014

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.

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.

Log in or register to post comments
Replied on Mon, 05/23/2016 - 15:23
Picture of user. AdminNeale Joined: 03/01/2013

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

Log in or register to post comments
  Login

News