Interpreting bivariate ACE output

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! 

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!

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?

2. 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!