Whole-genome regression methods are being increasingly used for the analysis and prediction of complex traits and diseases. In human genetics, these methods are commonly used for inferences about genetic parameters, such as the amount of genetic variance among individuals or the proportion of phenotypic variance that can be explained by regression on molecular markers. This is so even though some of the assumptions commonly adopted for data analysis are at odds with important quantitative genetic concepts. In this article we develop theory that leads to a precise definition of parameters arising in high dimensional genomic regressions; we focus on the so-called genomic heritability: the proportion of variance of a trait that can be explained (in the population) by a linear regression on a set of markers. We propose a definition of this parameter that is framed within the classical quantitative genetics theory and show that the genomic heritability and the trait heritability parameters are equal only when all causal variants are typed. Further, we discuss how the genomic variance and genomic heritability, defined as quantitative genetic parameters, relate to parameters of statistical models commonly used for inferences, and indicate potential inferential problems that are assessed further using simulations. When a large proportion of the markers used in the analysis are in LE with QTL the likelihood function can be misspecified. This can induce a sizable finite-sample bias and, possibly, lack of consistency of likelihood (or Bayesian) estimates. This situation can be encountered if the individuals in the sample are distantly related and linkage disequilibrium spans over short regions. This bias does not negate the use of whole-genome regression models as predictive machines; however, our results indicate that caution is needed when using marker-based regressions for inferences about population parameters such as the genomic heritability.
Whole-genome regression (WGR) methods [1] are becoming increasingly used for analysis and prediction of complex traits, quantitative or categorical. These methods were first developed for prediction in plant and animal breeding (e. g. , [2,3]). More recently, there has been an increased interest in using WGR methods for inferring the proportion of variance that can be explained by a linear regression on a marker panel, or ‘genomic heritability’ [4–6]. Prediction and inference are two different problems, and a model that may yield good (e. g. , unbiased and precise) estimates of parameters of interest may have a relatively poor prediction performance, and vice versa. Most of the methodological research in WGR methods was developed in animal breeding with a focus on prediction. Unfortunately, little is known about the inferential properties of estimates derived from WGRs models. For example, it is unclear whether the commonly used likelihood-based (or Bayesian) estimators of variance components or of genomic heritability estimate population parameters consistently [7]. Before copious marker information became available, genetic analysis (e. g. , estimation of heritability) was mainly based on mixed effects linear models applied to family data [8]. In the so-called infinitesimal model, relatedness due to kinship is assessed using pedigrees, and a central element of model specification is the assumption that genotypic values result from the small and additive effects of alleles at a large number of loci. A number of studies have investigated the quality of fit of the infinitesimal model to experimental [9,10] and to simulated family data [11]. Most of these studies have concluded that the additive infinitesimal model is a useful abstraction, except in situations involving a few large-effect non-additive loci. Therefore, at least at some operational level, with family information the distinction between the model that generated the data and the one used for analysis has not seemed critical. The availability of genotype information on a large number of loci has made possible to assess kinship among nominally unrelated individuals [9–13]. In this setting, and due to imperfect linkage disequilibrium (LD) between markers and quantitative trait loci (QTL), the patterns of allele sharing at markers and at causal loci may be very different [6]. Hence, the distinction between the data generating process and the model used for data analysis, or instrumental model, must be made clearly: in the instrumental model marker genotype information is used in lieu of the causal genotypes that are at the basis of the classical model of quantitative genetic theory. Thus, clarifying the link between the parameters of the instrumental model (e. g. , the genomic or SNP variance) and classical quantitative genetic parameters (e. g. , the genetic variance) is essential. Yang et al. (2010) [4] proposed using the G-BLUP method [2], a particular class of WGR, applied to data involving distantly related individuals, for estimation of the proportion of variance accounted for by a multiple-linear regression on common SNP. The proportion of unexplained genetic variance can be interpreted as ‘missing heritability’, which conceptually can be attributed to imperfect LD between markers and QTL. Using a WGR approach Yang et al. (2010) found that approximately half of the heritability of human height was captured by common SNP. Other studies, e. g. , [6], have corroborated Yang’s results using both simulated and real data. More recently, WGRs have been used for estimation in scenarios where all causal variants are assumed to be included in the marker panel, and various suggestions have been made with the purpose of obtaining inferences of genomic heritability that resemble more closely those based on pedigrees [5]. In the literature on genomic analysis of complex traits published so far [4–6,14], genetic parameters have been defined based on the statistical assumptions of the instrumental model used for data analysis. This is so despite the fact that there is a key difference in the way genotypes and effects are treated in statistical models and in quantitative genetics theory. In the latter, inter-individual differences in genetic values are attributed to subject-to-subject differences on allele content at QTL [15–17]; therefore genetic variance stems from variation at QTL genotypes. In this framework, at any given time in a population, the effects of alleles on a trait (e. g. , the average effects of allele substitution) are fixed quantities, e. g. , [16] pp. 112–113. On the other hand, in the instrumental regression models, genotypes are treated as fixed and variation stems from uncertainty about marker effects (the so called ‘variance of marker effects’). This key difference in the treatment of genotypes and of their effects has important consequences that we further explore in this article. An important contribution of this paper is to establish theory aiming at a precise definition of parameters arising in regression models using genomic data (markers, sequence) as explanatory variables. Our approach is framed within the classical quantitative genetics paradigm. We discuss how these “instrumental model parameters” relate to “structural parameters” of an underlying conceptual QTL model. We also present stylized cases that shed light into the interpretation of the parameters of the instrumental model. Finally, we discuss potential estimation problems and provide a limited set of simulations where some statistical properties of likelihood-based estimates are assessed.
In standard quantitative genetic theory [15–17] additive genetic values are linear functions of allele content at QTL. Concepts such as the additive effect of an allele in a population or narrow sense (trait) heritability are defined with reference to this framework. The standard quantitative genetic model (hereinafter referred to as QTL-model) assumes that a trait of interest measured on individual i (yi; i = 1, …, n) is affected by alleles at q QTL. Hereinafter, for ease of presentation, we assume that all loci are bi-allelic and that genotype codes (zij; j = 1,2, …, q) and phenotypes have been centered, so that E (yi) = 0 and E (zij) = 0 for all individuals and loci. The genetic value of an individual is defined as the expected phenotypic value given QTL genotypes, gi = E (yi|zi), where zi = {zi1, …, ziq} is a vector of genotype codes observed at the ith individual at each of the q QTL. The conditional expectation function maps from genotypes to expected phenotypic value. The genetic variance of a trait in the conceptual population is simply the average (over individuals) squared deviation of the genetic values from the population mean. In our setting, because E (yi) = 0, the genetic variance becomes σg2 = EzEy|z yi|zi2 = Ezgi2. Clearly, genetic variation stems from inter-individual differences in allele content at QTL (this is what confers individuals different genetic values). In the previous section we argued that the statistical parameter σu2 agrees with the population parameter σg2 only under highly simplified scenarios. Another important question is that of whether estimators derived from potentially misspecified likelihoods (either in likelihood-based or Bayesian estimates) can estimate parameters consistently, meaning that that the estimator converges to the true value of the parameter asymptotically. This will happen, for instance, if the bias and variance of the estimator go to zero as sample size goes to infinity, e. g. , [24]). We focus our study on the G-BLUP procedure due to its widespread use and relative simplicity. In this method, phenotypes are regressed on markers using a linear model of the form yi = ∑jxijbj+εi where bj~iidN0, σb2 and εj~iidN0, σε2, with marker effects and residuals mutually independent. The model implies the following marginal distribution of phenotypes y~N0, Gσu2+Iσε2 where N (. , .) stands for a multivariate normal distribution, G is a genomic-relationship matrix and σu2 is a variance parameter. Maximizing the likelihood (σu2 and σε2 are the unknown parameters) function associated with above-expression yields maximum likelihood estimates of variance components and of the proportion of variance explained by the model: hu2 = σu2σu2+σε2. Maximization of the likelihood function is facilitated using the eigenvalue decomposition of the G matrix; further details about this are given in S3 Text. Since G is computed using markers that are not necessarily QTL or that are in imperfect LD with QTL, the (co) variance patterns of additive effects and, consequently, the likelihood function, can be misspecified [6]. This potentially leads to large finite-sample bias and to inconsistency of estimates, meaning that the genetic parameters may not be well estimated even in large samples. To assess some statistical properties of maximum likelihood estimators (MLE) we conducted two simulations. In both, phenotypes were generated according to the additive QTL model of (Eq 1), yi = α′zi + δi, with QTL effects sampled from a zero-mean normal distribution and from an independent normal distribution with IID residuals, that is δi~N (0, σδ 2). Variance parameters were chosen to generate a trait heritability of 0. 5. In our first simulation markers and QTL were generated according to a stylized LD pattern. Our second simulation uses real human genotypes.
Above, we discussed conceptual problems emerging when genetic parameters are defined based on assumptions of statistical models that are not in line with fundamental quantitative genetic concepts. A second problem arises when estimation of these parameters is based on markers and not on QTL-genotypes. Under regularity conditions, maximum likelihood estimates are asymptotically unbiased ([24]). However, consistency cannot be guaranteed unless the likelihood is correctly specified and, even if the estimator happens to be consistent, misspecified likelihoods can induce sizable finite-sample bias. The proportion of allele sharing at any given set of loci can be viewed as a random variable with expected value given by twice the kinship coefficient between the individuals (derived from the full pedigree and with additional assumptions such as absence of mutation and no selection) and random variation given by the effects of Mendelian sampling [32]. If a large number of markers segregate independently from the QTL the proportion of allele sharing at markers and at QTL can be very different (e. g. , [6]). This is particularly important for pairs of individuals sharing short chromosome segments (e. g. , distantly related individuals from populations with large effective population size). Under these conditions a likelihood, constructed based on proportions of allele sharing at markers, can be largely misspecified and consistency may not hold. Of course, this does not necessarily imply that estimates would be inconsistent or have noticeable finite-sample biases; it just poses a caveat: one needs to be careful with use of estimates derived from models based on misspecified likelihoods, as it is most probably the case for marker-based model for whole-genome inference. There are at least two cases where a likelihood function based on markers will not be largely misspecified. The first one is when patterns of allele sharing at markers and at QTL are very similar; this can occur if markers are in tight LD with QTL, and also with data from close relatives. A second case is when the components of genetic values that cannot be explained by markers (ξi in expression (4) ) are IID, in which case the likelihood will correctly represent the (co) variance structure of the data. In these two cases the parameters of the instrumental model (e. g. , hu2) could be inferred consistently. Unfortunately, in general there is no good reason to believe that the ξi′s are IID and therefore the likelihood function will typically be misspecified. Goddard et al. [33] suggested a way of computing the marker-based genomic relationship matrix which, the authors argue, would have the property that the expected value of the proportion of allele sharing at QTL (GQTL) given the realized proportions of allele sharing at markers (GMRK) is GMRK, that is E[GQTL|GMRK] = GMRK. However, even if this property were to hold, this does not imply that estimates of variance components derived using GMRK would be unbiased. Indeed, in each sample GQTL can differ from GMRK; if the differences in the proportion of allele sharing at markers and at QTL are large enough, misspecification of the likelihood function due to the fact that we use GMRK instead of GQTL, can induce a systematic bias. Recently Jiang et al. [34] examined large sample properties of REML estimators of variance parameters of a marker-based regression (σε2, σu2) and concluded that, under certain conditions, the REML estimators can converge in probability to the true value of the statistical parameters, this being the case even if the likelihood is misspecified. The proof of this result is based on three key assumptions: (a) the model used to estimate parameters, what we here call the instrumental model, includes all QTL genotypes plus a number of markers with no effect, (b) all the covariates used in the model (i. e. , both QTLs and markers) are mutually independent, and (c) the number of QLT is not too small relative to the number of markers. In this specific setting the authors prove consistency of the REML estimator of the error variance and the convergence in probability of the REML estimator of the statistical parameter σu2. This result, which the authors verified in simulations, indicates some robustness of the REML estimator. On the other hand, in our study we detected several cases where misspecified likelihoods produced sizable finite-sample bias, this occurring in settings where ML estimates derived from correctly specified likelihoods were seemingly unbiased. However, there are important differences between the study of Jiang et al. [34] and ours. Firstly, Jiang et al. [34] focused on estimation of the statistical parameters of the instrumental model (σε2, σu2) while here we addressed estimation of parameters defined from a quantitative genetics perspective (σε2, σg2) and, as we discussed here, σu2 and σg2 coincide only under very stylized conditions (e. g. , complete LE between loci, a condition assumed by Jiang et al. [34] which we did not adopt here). Secondly, in an attempt to resemble what one encounters with modern genomic data, our simulations used a ratio between the number of QTL relative to the total number of loci that is relatively small (200/50K in simulation 1, and 5/305 in simulation 2). Importantly, Jiang et al. [34] indicate that the asymptotic results presented in their study involve approximations that may not hold when the ratio of number of QTL relative to the total number of loci is close to zero. Thirdly, our simulations incorporated LD (either LD blocks in simulation 1 or LD patterns realized in real human genotypes in simulation 2); therefore, in our case loci were not in mutual LE as assumed in Jiang et al. [34]. Fourthly, in the simulation used in [34] the ratio between the total number of loci (both marker and QTL) and sample size was 10; in our simulations this ratio was 50 in simulation 1 and approximately 60 in simulation 2. We adopted these settings because with modern genomic data the ratio between the number of markers and sample size is expected to be large. In the light of the evidence presented here it seems clear that further studies are needed to characterize the finite sample properties of REML estimators derived using misspecified likelihoods in settings that resemble the conditions encountered in the analysis of real genomes. Simulation studies using real human genotypes, including the one presented here, indicate that estimates of hu2 based on the G-BLUP model incorporating both markers and all QTL genotypes are lower than the trait heritability of the trait (h2). This reinforces the idea that GBLUP may not lead to unbiased estimates of the genomic heritability (hg2) because when all QTL are included in the marker panel hg2 = h2. Speed et al. [5] argued that a main reason for incorrect estimation is that LD is ignored in the computation of the G matrix. These authors suggested an alternative method for computing G that resulted in estimates of genomic heritability closer to the simulated trait heritability. However, in a follow up discussion [35,36] presented alternative simulations scenarios where the method of Speed et al. [5] yielded biased estimates. All in all, this suggests that the appropriate choice of method used for computing G may depend on the genetic architecture of the trait, a feature that is typically unknown, even if attention is restricted to additive gene action only. From our perspective, the main problem does not reside in the manner the G matrix is computed but, rather, in the use of massive numbers of markers that are in LE with QTL. The use of such data increases the sampling variance of estimates and may cause large finite-sample bias and lack of consistency. We have argued that the instrumental parameter hu2 may not provide a good representation of the genomic heritability (hg2); an important reason for the difference between these two parameters is that while LD plays a role in the determination of σg2 associations among marker genotypes are ignored in the definition of σu2. An alternative approach that accounts for multi-locus LD, could be based on the sample-variance of the true genomic values [37] that is σ~g2 = n-1∑i = 1n xi' β-x¯i' β2 where x¯i represents the average genotype. Of course, marker effects are unknown; however as suggested by Sorensen et al. [37] one could infer σ~g2 in a Bayesian setting by evaluating σ~g2 using samples from the posterior distribution of marker effects. The parameter σ~g2 accounts for LD between markers in a very specific manner: σ~g2 becomes equal to the genomic variance parameter, σg2 = β' Σxβ (see expression 7), with Σx replaced by its method-of-moments estimator X′ Xn-1. Indeed, when genotypes are centered, σ~g2 = n-1∑i = 1n xi' β2 = n-1∑i = 1nβxi xi' β = n-1β' ∑i = 1nxi xi' β = β' n-1X' Xβ. Complex traits are possibly affected by large numbers of small-effect QTL and the analysis of such traits requires fitting a large number of variants concurrently using a WGR approach such as the one proposed by [1]. Close relatives share long chromosome segments and, under these circumstances, the patterns of allele sharing at markers and at QTL are very similar. This leads to high prediction accuracy and very small bias in genomic heritability estimates. When markers and QTL co-segregate, variable selection does not seem to be needed [6]. On the other hand with distantly related individuals, the addition of large numbers of markers that are in LE with QTL can lead to incorrect specification of genomic relationships and this can result in potential inconsistencies of estimates of genomic heritability. We do not question the use of WGR for analysis of complex traits and as a prediction machine. Rather, we warn about problems arising when these methods are used for inferences. In our opinion this problem has been overlooked and oversimplified, and further research is needed to understand if and under what circumstances WGRs such as the G-BLUP can be used to correctly assess the true proportion of variance that can be explained by a regression on markers in the population.
Whole-genome regression (WGR) methods are being increasingly used for inferring the proportion of variance that can be explained by a linear regression on a massive number of markers, called ‘genomic heritability. ’ However, the statistical assumptions involved in WGRs are somewhat at odds with important quantitative genetics concepts. We argue and show that the parameters of the statistical model used for data analysis typically bear a tenuous relationship with the quantitative genetic parameters of interest. We also study, using simulations, the extent of bias of likelihood-based estimates. We conclude that under certain circumstances estimates can have a sizable finite-sample bias; therefore, caution needs to be exercised when interpreting parameter estimates derived from WGR models.
Genomic Heritability: What Is It?