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 asA/V
, which is the A matrix elementwise-divided by the phenotypic covariance matrix. Only the diagonal elements ofStandardizedA
--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] ofStandardizedA
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
andCorA
(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
, andste
. It sounds like you're specifically interested in the [2,1] element ofsta
andste
.I'm not sure what you mean here. If you're summing standardized variance components (from the diagonals of
StandardizedA
,StandardizedC
, andStandardizedE
), 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