How meaningful are heritability estimates of liability?

Abstract

It is commonly acknowledged that estimates of heritability from classical twin studies have many potential shortcomings. Despite this, in the post-GWAS era, these heritability estimates have come to be a continual source of interest and controversy. While the heritability estimates of a quantitative trait are subject to a number of biases, in this article we will argue that the standard statistical approach to estimating the heritability of a binary trait relies on some additional untestable assumptions which, if violated, can lead to badly biased estimates. The ACE liability threshold model assumes at its heart that each individual has an underlying liability or propensity to acquire the binary trait (e.g., disease), and that this unobservable liability is multivariate normally distributed. We investigated a number of different scenarios violating this assumption such as the existence of a single causal diallelic gene and the existence of a dichotomous exposure. For each scenario, we found that substantial asymptotic biases can occur, which no increase in sample size can remove. Asymptotic biases as much as four times larger than the true value were observed, and numerous cases also showed large negative biases. Additionally, regions of low bias occurred for specific parameter combinations. Using simulations, we also investigated the situation where all of the assumptions of the ACE liability model are met. We found that commonly used sample sizes can lead to biased heritability estimates. Thus, even if we are willing to accept the meaningfulness of the liability construct, heritability estimates under the ACE liability threshold model may not accurately reflect the heritability of this construct. The points made in this paper should be kept in mind when considering the meaningfulness of a reported heritability estimate for any specific disease.

Appendix

Consider the function $f({\hat{\mathbf{p}}}) = \arg \max_{\varvec{\uptheta} } \left\{ {\log L(\varvec{\uptheta} |{\hat{\mathbf{p}}},{\mathbf{n}})} \right\},$ which relates the estimated proportions of twin pairs, with (0,1, or 2) affected in each pair, to the maximum likelihood parameter estimates (i.e., $\hat {\varvec{\uptheta }}_{\text{ML}} = f({\hat{\mathbf{p}}})$). As indicated by the notation, $f( \cdot )$ is not a function of ${\mathbf{n}}$ because the model is saturated. By the law of large numbers, the estimated proportions converge in probability to the true proportions: ${\hat{\mathbf{p}}} \xrightarrow{P} {\mathbf{p}}^{\text{true}}$. Then, assuming f is continuous at ${\mathbf{p}}_{\text{o}}$, by Slutsky’s Continuity Theorem (Ferguson 1996) we have, $\hat {\varvec{\uptheta }}_{\text{ML}} \xrightarrow{P} f({\mathbf{p}}^{\text{true}} )$. Therefore, even in the case where the model is falsely specified, the maximum likelihood estimate will become very close to a function of the true proportions: $f({\mathbf{p}}^{\text{true}} )$. If the model is correctly specified, then the function of the true proportions will equal the true parameter values: $f({\mathbf{p}}^{\text{true}} ) = \varvec{\uptheta}_{\text{o}}$, otherwise it may be a biased estimate of the parameter values: $f({\mathbf{p}}^{\text{true}} ) \ne \varvec{\uptheta}_{\text{o}}$. In all of the scenarios discussed, it is easy to calculate the “true” proportions (${\mathbf{p}}^{\text{true}}$) using numerical integration. It is also possible to calculate the function $f({\mathbf{p}}^{\text{true}} )$ using numerical optimization techniques combined with numerical integration. Thus, the large sample properties of the estimator under the ACE liability model may be investigated by statistical theory without the use of simulations. We calculated the bivariate normal integral using the method of Genz and Bretz (2002) as implemented in the R package “mvtnorm,” and we performed numerical integration over the t distribution using the adaptive quadrature (Piessens et al. 1983) “integrate” function in R (R Team 2010). We performed numerical optimization using the “optim” function in R (R Team 2010). Also, we calculated the thresholds using the R function “uniroot” (R Team 2010).

As an example of calculating ${\mathbf{p}}^{\text{true}}$, consider the second scenario with a dichotomous common exposure where the exposure has frequency $\gamma$. Let $E \in \{ 0,1\}$ represent the exposure status of a pair. Consider the probability that neither of the twins is affected:

$$\begin{aligned} {\mathbf{p}}_{0j}^{\text{true}} &= \gamma P\left( {Z_{1} < \tau_{j} ,Z_{2} < \tau_{j} |E = 1} \right) \\&\quad + (1 - \gamma )P\left( {Z_{1} < \tau_{j} ,Z_{2} < \tau_{j} |E = 0} \right).\end{aligned} $$

(3)

Similar formulas to (3) can be written for ${\mathbf{p}}_{1j}^{\text{true}}$ and ${\mathbf{p}}_{2j}^{\text{true}}$, so the probability of 0, 1 and 2 affected in a twin pair can be calculated. If we allow $\beta = \sqrt {\sigma_{\text{C}}^{2} /\gamma \left( {1 - \gamma } \right)}$ to represent the effect size for the exposure, then variance of the common environment component is ${\text{Var(}}\beta E )= \sigma_{C}^{2}$. Furthermore, conditional on E we have

$$ \left[ \begin{gathered} Z_{1} \hfill \\ Z_{2} \hfill \\ \end{gathered} \right] = {\text{MVN}}\left( {\left[ {\begin{array}{*{20}c} {\beta E} \\ {\beta E} \\ \end{array} } \right],\left[ {\begin{array}{*{20}c} 1 & {\phi \sigma_{\text{A}}^{2} } \\ {\phi \sigma_{\text{A}}^{2} } & 1 \\ \end{array} } \right]} \right). $$

(4)

Thus, Eq. (4) may be used in conjunction with Eq. (3) to calculate ${\mathbf{p}}^{\text{true}}$. Very similar equations may be created for the other scenarios discussed in the paper.