Accurate estimation of heritability in genome wide studies using random effects models

1 INTRODUCTION

Many complex traits display a remarkable gap between the overall genetic variance estimated by population studies (such as twin studies) and the variance explained by specific genetic variants identified by genome wide association studies (GWASs). This gap has been coined the ‘dark matter’ of heritability (Manolio et al., 2008). Suggested explanations for this gap include gene–gene and gene-environment interactions (Manolio et al., 2008; Visscher et al., 2008), unidentified by the univariate GWAS scheme, though some of these explanations have been refuted by studies (Hill et al., 2008).

Height, being a trait that is easy to measure but genetically complex, has attracted special attention. Heritability of height, defined as the portion of height variability due to genetic factors, was estimated at around 80% (Macgregor et al., 2006; Visscher, 2008; Visscher et al., 2008) while the tens of specific single nucleotide polymorphisms (SNPs) identified by large-scale GWASs account for only 10% of the heritability (Gudbjartsson et al., 2008; Lettre et al., 2008; Weedon et al., 2008). Extrapolating from these data, under specific statistical assumptions, Goldstein (2009) estimated the number of causal SNPs affecting height at 93 000, which is at the high end of reasonable values.

(Yang et al., 2010) recently suggested the use of random effects models as a multivariate approach for estimating the heritability of height directly from GWAS data. In their model, each individual's height is affected by a genetic random effect, which is correlated across individuals by virtue of sharing some of the genetic variants affecting height, and an environmental random effect, which is uncorrelated across individuals. Consistent with the common definition, they define heritability as the ratio of the variance of the genetic random effect to the total variance.

Using a straightforward maximum-likelihood approach for estimating heritability in this model requires knowledge of the identity of the causal SNPs, and hence the covariance matrix. This is of course not available, so Yang et al. (2010) approximate the genetic correlation between every pair of individuals across the set of causal SNPs by the genetic correlation across the set of all genotyped SNPs. Using this approximation they estimate the heritability of height at 45%. By adding some further assumptions about the nature and distribution of the causative SNPs they obtain heritability estimates of 56% and 80%.

While the random effects approach proposed by Yang et al. (2010) is innovative and appropriate, the concept of using an estimated correlation matrix instead of the actual correlation matrix is a questionable heuristic. The very large number of SNPs used for estimating the genetic correlations—most of them likely not causative—might mask out the correlations on the set of causal SNPs. As we show below, this indeed leads to inaccurate and suboptimal estimation of heritability.

Instead, we propose here to treat the identity of the causal SNPs as missing data, and find true maximum-likelihood estimates of our parameters of interest in this setting. Towards this end, we develop an expectation-maximization (EM) algorithm for this problem. Because of the exponentially large data space, we employ Markov chain Monte Carlo (MCMC) methods for the estimation of the E-step. This approach is known as Markov-chain expectation-maximization (MCEM; Wei et al., 1990). We also develop an efficient optimization approach for the M-step, based on spectral decomposition of the estimated covariance matrix.

The resulting algorithm is computationally intensive, yet tractable for reasonably sized problems, as we demonstrate. As a by-product of our approach, we also get maximum-likelihood estimates of the portion of genotyped SNPs that are causative in the random effect model. This represents independently useful information, as it indicates whether the estimated heritability is likely due to a small number of large genetic effects, or a large number of smaller effects. For example, this knowledge may be useful in planning further studies to identify these specific factors.

We empirically demonstrate using simulations that our method, in addition to being theoretically sound, does indeed lead to significantly more accurate estimation of heritability than the heuristic approach of Yang et al. (2010).

2 METHODS AND SIMULATIONS

2.1 Random effect models of genetic causality

We first describe the random effects model proposed by Yang et al. (2010). Assume we have n individuals with their quantitative phenotype (say, height) denoted by y_j , j=1…n, and we know there are m genetic factors affecting the phenotype, as described by the following model:

where z_ij is the value of the i-th causative SNP of the j-th individual after standardization, i.e. where x_ij is the actual genotype and p_i is the minor allele frequency (MAF), so E(z_ij)=0 and V(z_ij)=1. u_i~N(0,σ_u²) is the effect of the i-th SNP on height. e_j~N(0,σ_e²) is the environmental effect which is assumed to be i.i.d. across individuals. μ is the mean term, that is typically assumed to be known and equal to 0 (i.e. it has been subtracted), and m is the number of causal SNPs.

Defining g_j=∑_i=1^mz_iju_i as the genetic random effect of individual j, the genetic variance is given by var(g_j)=σ_g²=mσ_u², and

From here it follows that y~N(0,Gσ_g²+Iσ_e²) where and Z is the genotypes matrix. G can be interpreted as the genetic correlation matrix. Heritability, or the share of the genetic variance in the overall variance, is given by . All that is needed for maximum-likelihood estimation of h² using standard random effects regression approaches (Searle et al., 1992) is the genetic correlation matrix G and the phenotypes vector y.

2.2 Estimation approach of Yang et al. (2010)

Yang et al. (2010) acknowledge the fact that the identity of the causal SNPs is unknown, and therefore G cannot be directly obtained. Assume that we have performed genome-wide genotyping on the n individuals, with a total of N genotyped SNPs on each. Yang et al. (2010)) propose to use the overall genetic correlations matrix computed on the entire genotyped data W (instead of just the ‘causative’ matrix Z) as an estimate of G. Yang et al. (2010) further improve this estimate in two ways.

First, Yang et al. (2010) correct the sampling bias along the diagonal of the matrix by defining A not as the cross-product of the genotype matrix, but by defining instead:

Second, Yang et al. (2010) account for linkage disequilibrium (LD) between genotyped SNPs and causal SNPs by replacing A with A*=β(A−I)+I where β is a parameter which is estimated from the genotyped data, using simulations which are based on assumptions as to the number of causal SNPs and their MAF. This can be thought of as regressing G on A and replacing A with the conditional expectation of G given A.

After obtaining the A* matrix, Yang et al. (2010) use restricted maximum likelihood (REML) for obtaining estimates of the variance components and hence heritability in the model y~N(0,A*σ_g²+Iσ_e²) instead of the real model y~N(0,Gσ_g²+Iσ_e²) for which the variance–covariance matrix is unknown. As we discuss later, this heuristic is statistically unbased and introduces substantial error to the estimation. We thus seek a more rigorous estimation approach.

2.3 Missing data representation of the problem

We wish to treat the identity of the causal SNPs as missing data, and therefore assume each SNP is causal with probability p. It is convenient to think about the set of causal SNPs as a bit mask represented by the number c with that bit mask as a binary representation. So 0 represents no causal SNPs while 2^N−1 means all SNPs are causal (where N is the number of genotyped SNPs), and the probability of observing a causal SNPs set c is p(c)=p^#c(1−p)^N−#c where #c is the number of causal SNPs in c. We denote by G(c) the genetic correlation matrix based on the causal SNPs in c. Given c, the log of the full likelihood can now be written as:

Since we do not know c we use EM (Dempster et al., 1977) to find maximum-likelihood estimates of our parameters of interest.

2.4 E-Step

The E-step takes the form:

From Bayes' theorem it follows that for all j:

where ϕ(y;Σ) is the multivariate normal density of the observed vector y given the variance-covariance matrix Σ.

2.5 MCMC estimation of the E-step

Clearly, it is not feasible to sum over 2^N summands for even modest values of N. We therefore use a MCMC scheme to estimate the conditional expectation of the log-likelihood. Given M samples from the conditional distribution of c: P(c|y;σ_g^2(t),σ_e^2(t),p^(t)) labeled c₁,…,c_M we estimate the E-step with the empirical mean :

To sample the conditional distribution of c, we specify the Metropolis algorithm in a manner similar to McCulloch, (1990). We specify the candidate distribution h_c(c) from which potential values are drawn and the acceptance function that gives the probability of accepting the new value.

Let c be the last draw from the conditional distribution of c with (b₁,b₂,…,b_N) being its binary representation. We generate a new value b_k^* for the k-th component of c using the candidate distribution h_c. We accept c*=(b₁,…,b_k^*,…,b_N) as the new value with the probability A_k(c,c*) given by:

Taking h_c(c)=p^#c(1−p)^1−#c the second term takes a particularly convenient form:

which is simply the ratio of the likelihoods.

2.6 M-step

Rewriting the empirical mean resulting from the MCMC:

indicates that the optimization in the M-step is done separately for p and for the variance components. For p we wish to solve:

Differentiating w.r.t p yields:

which is simply the average of causal SNPs proportions in the samples.

For the variance components, it is beneficial to rewrite the equation as a function of (r,σ_e²) with instead of (σ_g²,σ_e²). We then get:

Differentiating w.r.t. σ_e² yields:

Notice that plugging the last result into the E-step, the last term of the E-step expression becomes constant:

so the optimization problem for r becomes:

which has no closed form solution. However, given the spectral decomposition of each G(j) it is easy to solve numerically. G(j) is a correlation matrix and therefore has n real and non-negative eigenvalues. Let {λ_i(j)}_i=1ⁿ denote the eigenvalues of a matrix G(j) associated with eigenvectors {v_i(j)}_i=1ⁿ, then for all a,b:

and:

It is useful to denote w_i(j)=[y^Tv_i(j)]².

Hence, the maximized function is easy to compute when the spectral decomposition of each G(j) is known by setting a=r,b=1 and solving the following maximization problem:

Recall that the densities ϕ(y;G(j),σ_g²,σ_e²) are computed during the MCMC process for each G(j) matrix used. Computing the normal likelihood involves inverting the variance–covariance matrix G(j)σ_g²+Iσ_e² which takes O(n³) operations. Using the previous results, we can easily compute the log likelihood given the spectral decomposition by setting a=σ_g² and b=σ_e²:

Since obtaining the spectral decomposition is also done in O(n³), we can compute the spectral decomposition and use it for the likelihood calculations as well as for constructing the optimization problem for r with hardly any additional complexity cost. Doing so, the complexity of the optimization of r is negligible compared to the time complexity of the spectral decomposition and the MCMC burn-in.

3 RESULTS

We first study the convergence properties of the MCEM algorithm. As can be seen in Figure 1, convergence is achieved after a relatively small number of iterations, after which the estimators are almost identical, independent of their starting points. Clearly, some level of variance is expected due to the stochastic nature of the MCMC approximation of the E-step, and could therefore be reduced by either taking more samples or averaging over several consecutive MCEM iterations. The number of EM iterations required for convergence did not vary greatly even when the proportion of the causal SNPs was low (data not shown). We have also found that the lengths of burn-in periods, as determined by the Geweke (1992) algorithm varied greatly between simulations and between consecutive iterations of the MCEM algorithm. In general, we have found that burn-in periods were longer when the proportion of causal SNPs was high, and shorter when the estimators were closer to their true value.

The estimators of the proportion of causal SNPs converged much slower, but did display convergence, as can be seen in Figure 2. Starting from the 10 starting points that spread uniformly in the range (0.1%, 1%), all 10 estimators were in the range (0.4%,0.7%) within 7 iterations.

We noticed that the in-sample variances, and hence the in-sample heritability, varied greatly between simulations, especially when the number of causal SNPs was small. Thus, to compare the methods we first calculated the difference between the in-sample heritability (as estimated using the actual genetic correlation matrix G generated in our simulation) and the heritability estimates obtained by the Yang et al. (2010) method and our method. An advantage of this variability of in-sample heritability is that we, in fact, compared the methods over wide range of heritability values ranging from ~0.2 to ~0.8, and not just for 0.5.

When comparing the heritability estimators obtained after 5 MCEM iterations to the estimators obtained by the Yang et al. (2010) method (Fig. 3), we found that in sets A and B, where the initial estimators of the proportion of causal SNPs were set close to the in-sample proportion, our method significantly reduced the variance of the heritability estimation errors (p<0.03 in both cases using Bartlett's test (Snedecor and Cochran, 1989)).

It is important to note that while the performance of the method of Yang et al. (2010) deteriorates as the proportion of causal SNPs decreases, our method is more robust to such differences. The heritability estimation errors of Yang et al. (2010) in scenario A, where the proportion of causal SNPs is high, display a significantly lower variance than the heritability estimation errors of the same method in scenario B, when the proportion of causal SNPs is low (p<10⁻³). This is reasonable, as a lower proportion of causal SNPs means that the approximation of G used by the method of Yang et al. (2010) is less appropriate. Our method explicitly models the causal SNPs, and therefore we expected it to be less affected by this problem. Indeed, while our method does show higher estimation error variance in scenario B, the difference between the variances in scenarios A and B is not significant.

When initializing the estimate of the proportion of causal SNPs at a very different value, as in scenarios C and D, we saw a significant reduction in heritability estimation error from using our method only for scenario D(p<0.03). We suspect this is due to slow convergence of the estimators of p in scenario C, as discussed later.

It is also interesting to note that the estimation errors of the environmental random effects variance were smaller than the estimation errors of the genetic random effects variance. This was true for both methods and for all simulation scenarios, as well as for the convergence studies we conducted earlier. We suggest that this is because the environmental random effects are uncorrelated while the genetic random effects are correlated (with G as the correlation matrix) and thus provide a more difficult estimation problem.

4 DISCUSSION

The idea of using random effects models for analyzing genetic data is appealing and powerful. Even when only tag SNPs are genotyped, the high LD between tag SNPs and ungenotyped SNPs allows us to expect it to yield accurate results. However, the method suggested by Yang et al. (2010) might carry some hidden pitfalls. Replacing the correct genetic correlation matrix by a different matrix estimated from the data as if the latter matrix were the correct matrix is unfounded statistically. Indeed, using the A* matrix as the variance–covariance matrix introduces two additional sources of variance to the estimation. One is the error originating from the element-wise difference between A* and G (the noise from the non-causal variants) and one stems from the built-in variance of the estimated β. The former in particular is highly sensitive to the proportion of causal SNPs, as our simulations have demonstrated: with few causal SNPs, the noise in A* overwhelms the signal and estimation accuracy in the approach of Yang et al. (2010) deteriorates significantly. Ignoring both sources of variances would lead to inaccurate estimation and incorrect inference (overly optimistic confidence intervals).

Our approach avoids these issues by performing maximum likelihood-estimation of the relevant parameters: the genetic and environmental effects (variances) and the proportion of causal SNPs. To intuitively understand how this is accomplished, it is beneficial to think of the Metropolis algorithm used in the MCMC evaluation of the E-step as a ‘soft’ simulated annealing optimization scheme where we always switch to a SNP combination with higher likelihood, but switch to a lower likelihood combination with a probability that depends on the likelihood ratio. We thus expect the MCMC to visit often high-likelihood areas of the search space, corresponding to combinations of candidate causal SNPs that lead to good estimation. The M-step combines all these good combinations to estimate the parameters.

Convergence of EM in general, and MCEM in particular, to the true maximum-likelihood estimates is hard to guarantee and test (Boyles, 1983; Wu, 1983). Our simulations present evidence that the estimators of the variance components and heritability do, in fact, converge quickly to the neighborhood of the true value of the parameters, and provide some evidence that the estimator of the proportion of causal SNPs also converges. Some variance of the estimators is of course due to the Monte-Carlo nature of the MCMC approximation of the E-step and can be reduced either by taking more samples from the MCMC process, or by averaging consecutive EM iterations. More extensive research and simulations would be required to establish clear criteria for how long to run the algorithm in order to obtain optimal tradeoff between accuracy and computation time.

Putting aside the theoretical difficulties with convergence of EM, maximum-likelihood estimates possess many desirable qualities, which facilitate accurate estimation and statistical inference (Lehmann and Casella, 1998). In particular they are consistent and asymptotically normal. Moreover, once the process converges we can easily extract the information matrix (Meilijson, 1989), with which we can construct asymptotically accurate confidence intervals and carry out proper statistical inference.

Running times of the MCEM estimation are governed by two main factors: the number of individuals and the number of genotyped SNPs. Our current simulations have dealt with a toy example of 700 individuals and 10 000 SNPs. Current running times of the MCEM algorithm would render it difficult, but not impossible, to estimate the heritability when using a realistic number of genotyped SNPs. Since the lion's share of running time is ‘wasted’ on burn-in iterations (over 90% in most runs) before the MCMC process converges (as indicated by the convergence diagnostics of Geweke, 1992), we believe that smarter initialization of the MCMC process could reduce running times by an order of magnitude or more. Increasing the number of individuals is more problematic, as the time-complexity is cubic in the number of individuals, since every MCMC step involves computing the likelihood, and hence the inversion of a n×n matrix. We have shown that with even a small number of individuals our method produces very accurate estimates. Moreover, as noted earlier, the main bottleneck of the estimation is the eigenvector–eigenvalue decomposition, which is done in each MCMC iteration. The time complexity of this decomposition can be reduced from O(n³) to O(n²) using parallel algorithms (see, for example, Bientinesi et al., 2006). We therefore believe that using the MCEM method for heritability estimation in realistic databases is definitely feasible.

It is also appealing to think of replacing the ‘soft’ EM scheme we propose by a ‘hard’ EM, where in each M-step the parameters are estimated using the combination of SNPs with the highest likelihood encountered. In this case the Metropolis algorithm represents a simulated annealing optimization scheme. This setup has the advantage of also offering a possible identification of the set of causal SNPs, in addition to parameter estimation. We also expect such an approach to dramatically reduce running times. However, the hard approach is not a proper EM algorithm and is not as theoretically based as the standard EM approach. Therefore, it does not possess the favorable properties of EM estimators. We aim to study this alternative as means of reducing running times. Another alternative is to adopt a Bayesian approach altogether, as suggested by Guan,Y. and Stephens,M. (submitted for publication). The two approaches are not necessarily exclusive, as one might use the Bayesian estimate as a good starting point for the MCEM algorithm to reduce running times.

In this article, we have concentrated on the application of random effects linear regression models for estimating heritability of quantitative traits like height. While this application area is clearly important, many of the phenotypes of interest in current research are dichotomous, like disease status. The extension to the linear regression approach that would apply in such cases would be to apply random-effects generalized linear models (Breslow and Clayton, 1993), such as logistic regression. This direction presents significant complications because of fundamental differences in the statistical setup between linear random effects models and their generalized counterparts and the complexities involved in estimation methods of the latter (for example a MCEM algorithm by McCulloch, 1990). We intend to pursue this direction further in future work.