Interpreting bivariate ACE output

Posted on
No user picture. Ilaria Joined: 07/22/2024

Hi all, 

I am very new to genetic modelling and using ACE models. I have used Hermina's R script, slightly adapted to my question. I have having issues with the saturated models with constrained means (see other question I posted) but also in interpreting the ACE findings. First of all, in the ACE output I get negative A estimates which I would have thought would have been unexpected given the MZ/DZ correlations. 

Secondly, I am not sure I get the output from the ACE model, is it the first column first row (of for example SA) indicating the proportion of variance explained by additive genetic factor in trait 1 (i.e.,  0.17%) and the second row for trait 2 (1.9%)? However, I am not sure what the second row indicates. For example what is the SA=0.0283 indicating? The proportion of variance explained by additive genetic effect for the association between trait 1 and trait 2? Please find attached my script and my Ro file. 

Any help in understanding why do I get these negative numbers and how to interpret these output would be appreciated. Thanks !

         A       A      C      C      E      E      SA      SA     SC     SC     SE     SE
US  0.0038 -0.0158 0.5092 0.2834 1.7550 0.5635  0.0017 -0.0190 0.2245 0.3410 0.7738 0.6780
US -0.0158  0.0655 0.2834 0.5856 0.5635 1.6622 -0.0190  0.0283 0.3410 0.2532 0.6780 0.7185
Replied on Fri, 10/04/2024 - 09:36
Picture of user. AdminNeale Joined: 03/01/2013

Hi

Note that in this script the A C and E matrices are constructed from a Cholesky decomposition, by multiplying lower triangular matrix a by its transpose.  That forces the A C and E matrices to be positive definite, which isn't ideal from a statistical point of view.  See https://vipbg.vcu.edu/vipbg/Articles/PM30569348.pdf for why symmetric A C and E matrices that aren't constrained to be positive definite has statistical advantages.

In the SA columns the standardized estimates are arranged as a matrix, so .2245 and .2532 are the proportions of variation of A for the two traits, and -.0190 is the genetic correlation.

In the summary of the results, the diagonal of the matrix with columns A and A contains the estimated variance components for the two traits, and the -.0158 is the genetic covariance between the traits.  The other contributions to the covariance, .5092 from C and .5635 from E are positive.  Thus this is an example where the "bivariate heritability" idea fails.  In general, I don't think bivariate heritability is a useful statistic because it isn't bounded between 0 and 1 if the contributions to covariation are of mixed signs.  In this case it turns out negative.  All that is quite OK, but the bivariate heritability is clearly a nonsense statistic.  I generally prefer to consider the contributions to covariance without tacking on the term "heritability" since this is a proportion statistic that only works when all components have the same sign.

HTH!

Replied on Fri, 10/04/2024 - 12:48
No user picture. Ilaria Joined: 07/22/2024

Thanks for the paper. I was aware of this issue - from this blog and that exact paper. However, I became aware of this after pre-registering the analyses plan. We discussed to use the variance based models instead but the idea for now I think would be to fit this only as sensitivity analysis and check that we are indeed getting the same results to the Cholesky. 

Regarding your response: 

Thus this is an example where the "bivariate heritability" idea fails.  In general, I don't think bivariate heritability is a useful statistic because it isn't bounded between 0 and 1 if the contributions to covariation are of mixed signs.  In this case it turns out negative.  All that is quite OK, but the bivariate heritability is clearly a nonsense statistic.  I generally prefer to consider the contributions to covariance without tacking on the term "heritability" since this is a proportion statistic that only works when all components have the same sign.

Sorry if I am being slow to understand this. I am trying to use this paper to better interpret this output. 

Do you then mean that in my example -0.0158 is the genetic covariance (bivariate heritability) and -0.0190 is the genetic correlation (which are then distinct concepts?) ? 

Thanks again for all your help with this! 

Replied on Fri, 10/18/2024 - 10:38
Picture of user. AdminNeale Joined: 03/01/2013

In reply to by Ilaria

Yes, -0.0158 is the  genetic covariance estimate (but no it is not the bivariate heritability, that statistic is nonexistent/nonsensical here due to different signs among the three contributions to covariance) and -0.0190 is the genetic correlation estimate.  

Replied on Thu, 10/17/2024 - 11:06
No user picture. Ilaria Joined: 07/22/2024

Hi, 

Just following up on this and whether someone knows whether I am getting this distinction right from my data. Also how likely is that the values are indeed negative? 

Thanks a lot!

Replied on Thu, 11/14/2024 - 08:19
Picture of user. AdminNeale Joined: 03/01/2013

In reply to by Ilaria

In general, model fitting results should jive with the way the data look.  In your case, the small negative A variance estimates would seem likely to be due to slightly higher DZ than MZ correlations.  Does that check out with your data?

I'm not sure which distinction you mean (answering the thread has made prior posts invisible for now) but perhaps I addressed it in the first paragraph?

Replied on Mon, 06/23/2025 - 10:00
No user picture. Ilaria Joined: 07/22/2024

Hi, 

Apologies for the follow-up message and double posting (I have added this output elsewhere too). I have reran the models with the corrected ACE script adapted from Hermine's code, using the bivariate and multivariate (when I have two repeated outcomes) using matrix-based models. 

This is the script for my bivariate ACE model: 

# Matrices
meanG <- mxMatrix("Full", 1, ntv, T, svMe, labels = labVars("mean", vars), name = "meanG")
covA <- mxMatrix("Symm", nv, nv, T, valDiag(svPa, nv), label = labLower("VA", nv), name = "VA")
covC <- mxMatrix("Symm", nv, nv, T, valDiag(svPa, nv), label = labLower("VC", nv), name = "VC")
covE <- mxMatrix("Symm", nv, nv, T, valDiag(svPa, nv), label = labLower("VE", nv), name = "VE")

# Total variance and expected covariances
covP <- mxAlgebra(VA + VC + VE, name = "V")
covMZ <- mxAlgebra(VA + VC, name = "cMZ")
covDZ <- mxAlgebra(0.5 %x% VA + VC, name = "cDZ")
expCovMZ <- mxAlgebra(rbind(cbind(V, cMZ), cbind(t(cMZ), V)), name = "expCovMZ")
expCovDZ <- mxAlgebra(rbind(cbind(V, cDZ), cbind(t(cDZ), V)), name = "expCovDZ")

# Data
dataMZ <- mxData(observed = mzData, type = "raw")
dataDZ <- mxData(observed = dzData, type = "raw")

# Expectation
expMZ <- mxExpectationNormal("expCovMZ", means = "meanG", dimnames = selVars)
expDZ <- mxExpectationNormal("expCovDZ", means = "meanG", dimnames = selVars)
funML <- mxFitFunctionML()

# Models for MZ and DZ
dzModel <- mxModel(list(meanG, covA, covC, covE, covP), covDZ, expCovDZ, dataDZ, expDZ, funML, name = "DZ")
mzModel <- mxModel(list(meanG, covA, covC, covE, covP), covMZ, expCovMZ, dataMZ, expMZ, funML, name = "MZ")

# Multigroup fit
multi <- mxFitFunctionMultigroup(c("MZ", "DZ"))

# Standardisation and Correlations
matI <- mxMatrix("Iden", nv, nv, name = "I")
invSD <- mxAlgebra(solve(sqrt(I * V)), name = "iSD")
corA <- mxAlgebra(solve(sqrt(I * VA)) %&% VA, name = "rA")
corC <- mxAlgebra(solve(sqrt(I * VC)) %&% VC, name = "rC")
corE <- mxAlgebra(solve(sqrt(I * VE)) %&% VE, name = "rE")
corP <- mxAlgebra(solve(sqrt(I * V)) %&% V, name = "rP")

# Extract variance components
estUV <- mxAlgebra(cbind(VA, VC, VE), name = "UV", dimnames = list(rep("UV", nv), rep(c("VA", "VC", "VE"), each = nv)))
estSV <- mxAlgebra(cbind(VA/V, VC/V, VE/V), name = "SV", dimnames = list(rep("SV", nv), rep(c("SA", "SC", "SE"), each = nv)))
estCOR <- mxAlgebra(cbind(rA, rC, rE, rP), name = "COR", dimnames = list(rep("COR", nv), rep(c("rA", "rC", "rE", "rP"), each = nv)))

# Confidence Intervals
ciACE <- mxCI(c(
 paste0("SV[1,", 1:(3 * nv), "]"),
 paste0("SV[2,", 1:(3 * nv), "]"),
 paste0("COR[1,", 1:(4 * nv), "]"),
 paste0("COR[2,", 1:(4 * nv), "]"),
 "rA[2,1]", "rC[2,1]", "rE[2,1]", "rP[2,1]"
))

# Assemble the full ACE model
modelACE <- mxModel("ACE_bivariate_model",
                   meanG, covA, covC, covE, covP,
                   mzModel, dzModel,
                   matI, invSD,
                   corA, corC, corE, corP,
                   estUV, estSV, estCOR,
                   multi, ciACE)

# Run
fitACE <- mxTryHard(modelACE, intervals = TRUE)

# Output
fitGofs(fitACE)
fitEstCis(fitACE)

I have a few questions regarding the model outputs and would appreciate your recommendations on the next steps:

  1. Constraining C to Zero
    I have received conflicting advice on whether to constrain the shared environment component (C) to zero. On one hand, I was advised to constrain C to improve interpretability and reduce the risk of the covariance estimate exceeding 1 — although I am not convinced this would have a meaningful impact. On the other hand, I have been cautioned that omitting C could inflate the estimates for A (additive genetics) and E (unique environment).

Model fit comparisons between the ACE and AE models do not show major differences. Do you have a preference between including or excluding C in such cases? Most importantly, could you recommend a paper or reference that could help guide or justify this decision?

  1. Standardising Covariances
    Should the covariance be standardised so that it is bounded between 0 and 1? My initial thought was that doing so would result in values equivalent to correlations. However, this paper suggests that covariances and correlations are distinct concepts and that both should be reported https://link.springer.com/article/10.1007/s10519-021-10046-y#Sec20 

I would really appreciate your thoughts on whether both measures should be included and how best to present them.

If there is no standard way to scale or standardise the A, C, and E covariance estimates, I am uncertain about the most appropriate way to report these results clearly and meaningfully.

Thank you all for your help!  

 

SA (95%CI) 

SC (95%CI) 

SE (95%CI) 

COV A 

COV C 

COV E 

Corr A 

Corr C 

Corr E 

Main analysis: Bivariate model (using a variance-based model) exploring the association between exposure and outcome age 21 

Exposure 

0·66 (0·48 to 0·85) 

 

-0·11 (-0·27 to 0·04) 

 

0·45 (0·39 to 0·51) 

 

·· 

 

·· 

 

·· 

 

·· 

 

·· 

 

·· 

 

Outcome age 21  

 

0·73 (0·48 to 0·98) 

 

-0·22 (-0·43 to -0·02) 

 

0·49 (0·42 to 0·58) 

 

·· 

 

·· 

 

·· 

 

·· 

 

·· 

 

·· 

 

Exposure to Outcome

·· 

·· 

 

·· 

 

1·06 (0·73 to 1·41) 

 

-0·26 (-0·56 to 0·02) 

 

0·20 (0·10 to 0·31) 

 

0·77 (0·57 to 0·97) 

 

NA (NA to NA) 

0·21 (0·11 to 0·32) 

 

Main analysis: Multivariate model (using a variance-based model) exploring the association between exposure and outcomes ages 21 and 26 

Exposure 

0·67 (0·48 to 0·85) 

-0·11 (-0·28 to 0·04) 

0·45 (0·39 to 0·51) 

·· 

 

·· 

·· 

·· 

·· 

·· 

Outcome age 21  

 

0·33 (0·05 to 0·62) 

-0·02 (-0·25 to 0·21) 

0·68 (0·59 to 0·78) 

·· 

 

·· 

·· 

·· 

·· 

·· 

Outcome age 26

0·42 (0·20 to 0·64) 

0·05 (-0·19 to 0·27) 

0·53 (0·45 to 0·63) 

·· 

 

·· 

·· 

·· 

·· 

·· 

Exposure to Outcome age 21

·· 

 

·· 

 

·· 

 

0·64 (-0·07 to 1·34) 

0·01 (-0·58 to 0·59) 

0·35 (0·13 to 0·58) 

0·42 (0·23 to 0 ·62) 

NA (NA to NA) 

 

0·08 (-0·04 to 0·20) 

Exposure to Outcome age 26

·· 

 

·· 

 

·· 

 

1·05 (0·35 to 1·77) 

-0·23 (-0·84 to 0·35) 

0·18 (-0·04 to 0·41) 

0·40 (0·17 to 0·63) 

NA (NA to NA) 

 

0·11 (-0·01 to 0·22) 

Outcome age 21 to Outcome age 26 

·· 

 

·· 

 

·· 

 

0·79 (0·37 to 1·20) 

-0·13 (-0·48 to 0·20) 

0·35 (0·22 to 0·49) 

0·94 (0·82 to 1·07) 

NA (NA to NA) 

 

0·45 (0·35 to 0·55)