## You are here

The OpenMx website will be down for maintenance from 9 AM EDT on Tuesday, September 17th, and is expected to return by the end of the day on Wednesday, September 18th. During this period, the backend will be updated and the website will get a refreshed look.

# Multivariate model with positive and negative correlations

6 posts / 0 new
Offline
Joined: 11/26/2015 - 11:14
Multivariate model with positive and negative correlations

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")
dataMZF <-mxData(observed=mzfData, 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

Offline
Joined: 01/24/2014 - 12:15
MxAlgebra "h2"

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.

Offline
Joined: 11/26/2015 - 11:14
Thank you Rob. Re " I don't

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?

Offline
Joined: 01/24/2014 - 12:15
My understanding is that
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.

Offline
Joined: 03/01/2013 - 14:09
Proportion is not the right term

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

Offline
Joined: 11/26/2015 - 11:14
Thank you both for the

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