Hi,
I am trying to calculate the heritability for some fat-related traits. I have fitted the ACE model and some submodels (CE and AE) and then learned that it was recommended to use mxModelAverage()
to gain the model-averaged point estimates.
Here is part of the script I am using to fit ACE model:
modelbase <- mxModel("base", mxMatrix(type="Full",nrow=1,ncol=ncv,free=TRUE, values=0,labels=c('Bage','BSUBTOT_MASS'), name="B" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, labels="mean", name="meanG" ), # Matrices declared to store a, c, and e Path Coefficients mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=sqrt(svace), label="a11",lbound=0.0001, name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=sqrt(svace), label="c11",lbound=0.0001, name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=sqrt(svace), label="e11",lbound=0.0001, name="e" ), # Matrices generated to hold A, C, and E components and total Variance mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), mxAlgebra( expression=A+C+E, name="V" ), mxAlgebra( expression=A/V, name="h2" ), mxAlgebra( expression=C/V, name="c2" ), mxAlgebra( expression=E/V, name="e2" ))
And the code I'm using to calculate the model-averaged point estimates is as follows:
mxModelAverage(reference=c("base.h2[1,1]","base.c2[1,1]","base.e2[1,1]"), include="onlyFree", type="AICc", SE=T, models=list(AceFit,AEfit, CEfit))
I learned that, for the purpose of statistical inferences, it's recommended to set include="onlyFree"
. However, in this way, the sum of model-averaged point estimates of h2, c2 and e2 is not equal to 1. The result is showed as follows:
$`Model-Average Estimates` Estimate SE base.h2[1,1] 0.6400199 0.14768808 base.c2[1,1] 0.2225559 0.11887596 base.e2[1,1] 0.2259699 0.02830278 $`Model-wise Estimates` ACE AE CE base.h2[1,1] 0.5502440 0.7759045 NA base.c2[1,1] 0.2225497 NA 0.6412221 base.e2[1,1] 0.2272063 0.2240955 0.3587779 $`Model-wise Sampling Variances` ACE AE CE base.h2[1,1] 0.0154727014 0.000742769 NA base.c2[1,1] 0.0141290903 NA 0.001150016 base.e2[1,1] 0.0008354342 0.000742769 0.001150016 $`Akaike-Weights Table` model AICc delta AkaikeWeight inConfidenceSet 1 ACE 995.5385 0.0000000 6.021586e-01 * 2 AE 996.3675 0.8289793 3.978325e-01 * 3 CE 1017.7799 22.2413237 8.913921e-06
Thus, what do you think can be done to fix it? Using include="all"
at the cost of violating normal sampling distribution?
Additionally, according to the code of the ACE model, it seems also feasible to calculate model-averaged point estimates of path coefficients a, c, e and variance of A, C, E. However, this will also result in different estimates of heritability. In this case, what parameters do you recommend to be calculated model-averaged point estimates?
Many thanks in advance!
I think it makes the most sense to report model-averaged estimates for the raw variance components--that is, the unstandardized variances of A, C, and E.
Many thanks for your suggestion!
Just to check, if I calculate the model-averaged point estimates of A, C, and E, then I can only get the point estimate of heritability, rather than the CI of it. Is it right? Thank you.
Yes, that's true, if you use
include="onlyFree"
andrefAsBlock=FALSE
. If you instead useinclude="all"
andrefAsBlock=TRUE
, you'll get the joint sampling covariance matrix of the model-averaged A, C, and E variance components. You can then use that joint covariance matrix to calculate the standard error of the heritability coefficient, using the delta method.Per Chebyshev's Inequality, ±5SEs will give you a confidence interval with coverage probability of at least 96%.