# CI estimation of Trivariate Cholesky model

CholACE.StandardizedA[1,1] 4.835575e-02 0.4772359455 8.043944e-01

CholACE.StandardizedA[2,1] -3.001322e+03 1.6422971828 1.889364e+03

CholACE.StandardizedA[3,1] -4.236242e+02 -0.2316164474 1.479296e+00

The genetic correlation seemed also unreasonable:

CholACE.CorA[1,1] 1.000000e+00 1.0000000000 1.000000e+00 !!!

CholACE.CorA[2,1] -9.999988e-01 -0.9999978673 9.999902e-01

CholACE.CorA[3,1] -9.999983e-01 0.9999635762 9.999960e-01

I also checked my dataset which contained no singular value. I totally have no idea about this situation. Could anyone help me find the problems of my codes (as follow)? Many thanks!

meVals <- c(0.5,0.5,0.5) # start value for means

meanLabs <- paste(Vars,"mean",sep="") # create labels for means

paVal <- .6 # start value for path coefficient

paLBo <- .0001 # start value for lower bounds

paValD <- vech(diag(paVal,nv,nv)) # start values for diagonal of covariance matrix

paLBoD <- diag(paLBo,nv,nv) # lower bounds for diagonal of covariance matrix

paLBoD[lower.tri(paLBoD)] <- -10 # lower bounds for below diagonal elements

paLBoD[upper.tri(paLBoD)] <- NA # lower bounds for above diagonal elements

```
```# Create Labels for Lower Triangular Matrices

aLabs <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")

cLabs <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")

eLabs <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")

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

pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

values=paValD, labels=aLabs, lbound=paLBoD, name="a" )

pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

values=paValD, labels=cLabs, lbound=paLBoD, name="c" )

pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,

values=paValD, labels=eLabs, lbound=paLBoD, name="e" )

# 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 total variances and standard deviations (diagonal only)

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

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

invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,

values=meVals, labels=meanLabs, name="Mean" )

meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" )

covMZ <- mxAlgebra( expression=

rbind( cbind(A+C+E , A+C),

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

covDZ <- mxAlgebra( expression=

rbind( cbind(A+C+E , 0.5%x%A+C),

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

# Standeraized parameters and CI

sta <- mxAlgebra(iSD %*% a, name="sta")

stc <- mxAlgebra(iSD %*% c, name="stc")

ste <- mxAlgebra(iSD %*% e, name="ste")

StandardizedA <- mxAlgebra(A/V, name="StandardizedA")

StandardizedC <- mxAlgebra(C/V, name="StandardizedC")

StandardizedE <- mxAlgebra(E/V, name="StandardizedE")

CorA <- mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA")

CorC <- mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC")

CorE <- mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE")

CorP <- mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP")

CI <- mxCI(c('sta', 'stc', 'ste','StandardizedA', 'StandardizedC', 'StandardizedE','CorA', 'CorC', 'CorE', 'CorP'))

# Data objects for Multiple Groups

dataMZ <- mxData( observed=mzData, type="raw" )

dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups

objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean",

dimnames=selVars )

objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean",

dimnames=selVars )

funML <- mxFitFunctionML()

`# Combine Groups`

pars <- list( pathA, pathC, pathE, covA, covC, covE, covP,

matI, invSD, meanG, meanT, sta, stc, ste,

StandardizedA, StandardizedC, StandardizedE,

CorA, CorC, CorE, CorP)

modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, funML, name="MZ" )

modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, funML, name="DZ" )

fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") )

CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, fitML, CI)

## Although it run smoothly, the

Your MxAlgebra

`StandardizedA`

is defined as`A/V`

, which is theAmatrix elementwise-divided by the phenotypic covariance matrix. Only the diagonal elements of`StandardizedA`

--elements [1,1], [2,2], and [3,3]--are readily interpretable; they are, of course, the standardized variance components (i.e., variance proportions). Elements [2,1] and [3,1] of`StandardizedA`

are not bounded between 0 and 1.If you're looking for the genetic correlations among the three phenotypes, you need to be looking at elements [2,1], [3,1], and [3,2] (or equivalently, [1,2], [1,3], and [2,3]) of

`CorA`

. The CI you report for the [1,1] element is correct because that's a diagonal element of a correlation matrix, and is therefore fixed at 1.0. The other two appear OK, in that they are not flagged with "!!!". How large is your sample? If it is relatively small, it's quite possible to get an extreme estimate of genetic correlations, and for a valid CI for the genetic correlations to span nearly the whole parameter space--see this post for a similar situation.Log in or register to post comments

In reply to Although it run smoothly, the by AdminRobK

## How to estimate the CIs of a31, a21 and a32?

And yes, my sample size is very small that contains only 20 pairs of MZ and 15 pairs of DZ. How does the sample size affect the genetic correlations? Is it nessesary to calculate the statistical power? Thanks~ :D

Log in or register to post comments

In reply to How to estimate the CIs of a31, a21 and a32? by Liz

## how could I estimate the CIs

What do you mean by "a31," "a21," and "a32?" It looks like you're requesting CIs for all the elements of

`StandardizedA`

and`CorA`

(and for that matter,`sta`

), so you should be seeing entries for the elements you want in the confidence-interval output. Right?That's tiny! With a sample that small, you're going to get very statistically imprecise estimates of the parameters of interest, so it's not that surprising that you're seeing extremely wide CIs for the genetic correlations.

Log in or register to post comments

In reply to how could I estimate the CIs by AdminRobK

## The genetic contribution (A) from phenotype A to phenotype B

Log in or register to post comments

In reply to The genetic contribution (A) from phenotype A to phenotype B by Liz

## Refers to the "a21" an "e21",

The standardized path coefficients of the Cholesky parameterization are in your MxAlgebras

`sta`

,`stc`

, and`ste`

. It sounds like you're specifically interested in the [2,1] element of`sta`

and`ste`

.I'm not sure what you mean here. If you're summing standardized variance components (from the diagonals of

`StandardizedA`

,`StandardizedC`

, and`StandardizedE`

), you're going to get a sum of 1.0 for sure.Log in or register to post comments

In reply to Refers to the "a21" an "e21", by AdminRobK

## Thanks

Log in or register to post comments

In reply to Refers to the "a21" an "e21", by AdminRobK

## The contribution of genes to phenotypic correlation

Log in or register to post comments

In reply to The contribution of genes to phenotypic correlation by Liz

## Formula?

Note that none of the path coefficients from the Cholesky factorization appear in the above formula; only the genetic correlation and the (raw) variance components do.

Source: Posthuma, D. (2009). Multivariate genetic analysis. In Kim, Y-K. (Ed.),

Handbook of Behavior Genetics(pp. 47-59). New York: Springer Science+Business Media, LLC.Log in or register to post comments

In reply to Formula? by AdminRobK

## Thanks

Log in or register to post comments

## Many thanks, Rob. I think I

Log in or register to post comments