A general unified framework to assess the sampling variance of heritability estimates using pedigree or marker-based relationships.

Heritability is a population parameter of importance in evolution, plant and animal breeding, and human medical genetics. It can be estimated using pedigree designs and, more recently, using relationships estimated from markers. We derive the sampling variance of the estimate of heritability for a wide range of experimental designs, assuming that estimation is by maximum likelihood and that the resemblance between relatives is solely due to additive genetic variation. We show that well-known results for balanced designs are special cases of a more general unified framework. For pedigree designs, the sampling variance is inversely proportional to the variance of relationship in the pedigree and it is proportional to 1/N, whereas for population samples it is approximately proportional to 1/N(2), where N is the sample size. Variation in relatedness is a key parameter in the quantification of the sampling variance of heritability. Consequently, the sampling variance is high for populations with large recent effective population size (e.g., humans) because this causes low variation in relationship. However, even using human population samples, low sampling variance is possible with high N.

HERITABILITY (h2), the proportion of phenotypic variation that is explained by additive genetic variation, is an important parameter in plant and animal breeding, evolutionary genetics, and human and medical genetics. It is central in quantifying the role of genetics in complex traits, predicting response to selection in natural and artificial breeding programs, and determining the limits of trait or disease prediction using information from relatives or DNA markers. Traditionally, the estimation of heritability is from pedigree data, by modeling the observed resemblance between relatives (Falconer and Mackay 1996; Lynch and Walsh 1998). More recently, genetic variation has been estimated using genetic marker information (Ritland 2000; Thomas 2005; Visscher ; Yang ; Robinson ; Berenos ). These designs estimate the genetic variance explained by the markers, which may be less than the additive genetic variance (Yang ), but in this article we refer to the parameter estimated as the heritability regardless of whether it is estimated from relationships defined by pedigree or by markers. In general, designs to estimate heritability can be grouped by their use of (i) the expected identity-by-descent (IBD) sharing between relatives, i.e., using pedigree relationships, (ii) marker-based estimated IBD relationships between relatives for known pedigree relationships, and (iii) marker-based estimated genomic relationship matrices for unknown pedigree relationships. For a review of these designs with a particular focus on human populations, see Vinkhuyzen .

Even with large sample sizes, the standard error of heritability estimates is often disappointingly large and it varies greatly between experimental designs. Therefore it is important to calculate the expected standard error before committing resources to collecting the data. Given a particular experimental design and the population value of h2, its sampling variance can be determined using a number of methods. After the data have been collected, the (asymptotic) sampling variance of the estimate can be derived from the analysis, for example, from mean squares in balanced designs, from the information matrix when using maximum likelihood or from the posterior density in Bayesian analysis. Prior to collecting data on phenotypes, the sampling variance can be predicted using statistical theory, typically for balanced designs, or obtained from computer simulation for more complex pedigree structures. In this study, we provide a single framework for calculating the asymptotic sampling variance of the heritability across a wide range of designs, for a class of models with two random variables and when analysis is by maximum likelihood (ML). We derive the sampling variance using the expected value of the information matrix. We show that previous results are special cases of the general framework and that the variance in relationships in the sample is a key parameter in all experimental designs.

Model and Assumptions

We assume a linear model with no fixed effects (or fixed effects that have been adjusted for without error) and two random components, a genetic effect (g), and a residual effect (e). There are N individuals, each with a single observation, y,where y, g, and e are vectors of length N of the phenotypic observations, genetic value, and residuals, respectively. G is the genetic relationship matrix (GRM), either from pedigree relationships, in which case it is the usual numerator relationship matrix (twice the kinship matrix), or derived from SNP similarity (Vanraden 2008; Stranden and Garrick 2009; Yang ). The genetic, residual, and total variances are σg2, σe2, and σ2, respectively. The N × N covariance matrix of all observations (V) iswhere h2 = σg2/(σg2 + σe2) = σg2/σ2, the heritability.

General Formula for Sampling Variance

We can decompose the symmetric GRM aswith TT′ = T′T = I and T−1 = T′ because T is orthogonal and D a diagonal matrix containing eigenvalues (λ) of G. Inference on h2 from data y does not change upon a linear transformation of y. We can therefore transform y by using the eigenvectors of G, which for the simple model used here are also eigenvectors of V (Thompson and Shaw 1990, 1992; Lippert ; Blangero ; Raffa and Thompson 2014).

withThe log likelihood with respect to h2 and σ2 isas shown previously (Thompson and Shaw 1990; Raffa and Thompson 2014). Equation 1 is very similar to that in Blangero , but with added parameter σ2. Elements of the (Fisher) information matrix (F) are obtained by taking the second derivative of (1) taken at the maximum with respect to h2 and σ2, and then the negative value of its expectation over y*, usingThe derivation of the first element of F (F11) is given here. The other two elements are derived analogously,and soThe resulting elements of the 2 × 2 matrix F arewith constants a and b,andThese elements are similar to those presented in Thompson and Atkins (1994), who parameterized the likelihood in a genetic and residual variance component, whereas we have parameterized in heritability and phenotypic variance. Thompson and Atkins do not have the factor and have λi2 and λi in the equations above where we have (λ – 1)2 and (λ – 1), respectively, the difference due to the choice of parameters in the model. In the article that developed the method of estimation of variance component in linear mixed models using restricted maximum likelihood (Patterson and Thompson 1971), the authors presented both the log likelihood and the information matrix in terms of eigenvalues of the covariance matrix.

The asymptotic sampling (co)variance for the estimates of heritability and phenotypic variance are from F−1. Therefore, the asymptotic sampling variance of the estimate of the heritability isHence, under the assumptions given, this is a completely general expression for the asymptotic sampling variance of an estimate of heritability and depends only on the eigenvalues of the GRM, the population value of heritability, and the experimental sample size.

Special Cases

With additional assumptions or for balanced designs, terms for a and b simplify and simple solutions for the sampling variance of can be derived. We go through a number of these special cases in this section that encompass pedigree and marker-based GRM.

Phenotypic variance (σ2) known

In many applications, the sampling variance of the total phenotypic variance is small or known before the experiment is conducted, and therefore it is useful to consider the sampling variance of heritability under the assumption that the phenotypic variance is known without error. For example, Blangero assume that σ2 is known in their derivations of the expected likelihood-ratio-test statistic (ELRT). If we assume here that the phenotypic variance is known without error then the resulting sampling variance of the estimate of heritability isThis expression is smaller than that in (2); hence assuming that phenotypic variance is known when it is not will lead to an underestimate of the sampling variance of heritability. This underestimate will be small when b2/N is small relative to the term a.

For a small heritability, a → b → andAssuming that the phenotypic variance is known and h2 is small giveswhich is close to (4) because the mean eigenvalue will be 1 in the absence of inbreeding when the GRM is from pedigree identity-by-descent and very close to 1 when the GRM is estimated from SNP data (Janss ). Hence, when the population value of heritability is small, its sampling variance is only a function of the variation in relatedness and sample size.

All

Equation 4 is also the result for when all λ are close to 1, such that their variance approaches zero. This situation can occur when the GRM is created from population SNP data on unrelated individuals in a population with a large effective population size. However, as we derive below, the variance of eigenvalues depends both on experimental sample size and effective population size, and so these parameters affect the sampling variance of heritability. In particular, the variance in eigenvalues is proportional to experimental sample size, so the larger the sample size the wider the spread around a mean value of 1.

Pairs of relatives with relationship r

If there are m pairs of relatives of the same degree r, then 2m = N and there are m eigenvalues λ1 with value 1 + r and m eigenvalue λ2 with value 1 − r (Searle 1982; Blangero ). Let ρ = rh2. ThenandFor pairs of monozygotic (MZ) twins (r = 1), Equation 5 becomes var() = (1 – ρ2)2/m. For pairs of full-sibs (r = ), the sampling variance is 4(1 – ρ2)2/m. For bivariate normality, the sampling variance of a correlation coefficient between two variates with population value ρ is ∼(1 – ρ2)2/N (e.g., Lynch and Walsh 1998, p. 819), so consistent with Equation 5.

Balanced design of multiple families

For m families with n individuals of relationship r, there are (n − 1) eigenvalues of (1 – r) and 1 eigenvalue of (1 + r(n – 1)) per family. This follows from known results on eigenvalues for symmetrical matrices that can be written as cI + dJ, with c and d constants (Searle 1982). Substituting these eigenvalues into the equation for parameters a and b givesandThis is consistent with the intraclass correlation sampling variance (e.g., Falconer and Mackay 1996, p. 180), apart from having m in the denominator [the least-squares derivation has (m − 1) instead]. Although we have assumed no fixed effects, in practice at least a mean would be included in the model and this absorbs one degree of freedom from the comparison of families. The least-squares formula takes account of this but ML estimation ignores it. Assuming that the phenotypic variance is known givessmaller than (5) by a factor of 1/(1 + ρ2(n − 1)). For large half-sib families, this term can be substantial.

Twin design

In human populations, the classical twin design is common for estimating genetic and nongenetic variance components. Let N = 2mM + 2mD, with mM and mD the number of MZ and dizygotic (DZ) pairs, respectively. In total, there are four different eigenvalues: 2, 0, 3/2, and 1/2 (Blangero ), with multiplicity mM, mM, mD, and mD. Let c = mM/(mM + mD), the proportion of all twin pairs that are MZ pairs. Using Equation 5, a − b2/N = NT, withandThis analysis assumes that there are no common environmental effects so the sampling variance is not appropriate for the usual practice of estimation of heritability using maximum likelihood fitting both an additive genetic and common environmental component (Neale and Cardon 1992).

Within-family estimation using realized relationships estimates from markers

Full-sibs have an expected pedigree relationship of 0.5 but the actual amount of the genome shared varies around 0.5 and this realized relationship can be estimated using genetic markers and used to estimate heritability (Visscher , 2007; Hemani ). These relationships can be estimated using identity-by-descent calculations conditional on observed marker genotypes. For full-sibs and half-sibs in human populations, the standard deviation of realized relationships is ∼0.04 and 0.03, around the expected value of and , respectively. For a comprehensive theory on the variance of realized relationships, see Hill and Weir (2011). A feature of this design is that common environmental factors that vary between families do not bias the heritability estimate. Visscher derived an approximate sampling variance of the estimate of heritability from multiple families with two full-sibs each. Hill (2013) derived the sampling variance of the estimate of genetic variation using REML for the general case of f families each of size n and expected relationship θ (twice the kinship coefficient). We can use the same general framework as developed here to approximate the sampling variance from within-family estimation. The difference between this design and those previously discussed is that the GRM is not fixed. That is, the eigenvalues of the GRM are themselves random variables and to derive the sampling variance of the estimate of heritability we need to first derive the expected value of the elements of the Information matrix over repeated samples. We provide details of an approximation in Appendix A. It results inThis equation shows that the sampling variance reduces by the square of the sample size per family (n), essentially because every individual adds a contrast with all other family members in the sample. As detailed in Appendix A, this approximation breaks down when h2 and n are large.

Random sampling from the population

One design to estimate the amount of additive genetic variation captured by SNPs is to take a random sample of individuals from the population, derive a GRM from SNP similarity, and estimate variance components from (residual) maximum likelihood (Yang ). In this sampling scheme, individuals are not sampled or ascertained based upon particular pedigree relationships, and any pedigree relationship, if known, is not taken into account in the analysis. The sampled individuals are related to some extent, even if very distantly, because the population size is finite. In human populations, this sampling scheme corresponds to sampling individuals who are conventionally unrelated. As for the case of realized relationships within families, the GRM is not fixed. We approximate E(a) and E(b2) in Appendix B. The resulting sampling variance of the estimate of heritability iswhere v(θ) is the variance of relatedness in the population, which is a function of effective population size (Goddard 2009; Goddard ). Analogous to the within-family design, the sampling variance is inversely proportional to the square of the sample size, rather than by 1/N in pedigree designs. Rijsdijk and Sham (2002) derived the same result (parameterized as the noncentrality-parameter, NCP, of the test statistic for heritability) for QTL linkage mapping in pedigrees, assuming that the variance in relatedness is small. Equation 8 was previously derived for SNP-based estimation of variance components from linear regression theory, assuming that the phenotypic variance is known without error (Vinkhuyzen ; Visscher ).

Statistical Power

The interest in this study is not about hypothesis testing but about quantifying the sampling variance of the estimate of heritability. For a detailed treatment on statistical power in variance component estimation using (restricted) maximum likelihood we refer to previous publications (Self and Liang 1987; Shaw 1987; Thompson and Shaw 1990; Almasy and Blangero 1998; Williams and Blangero 1999; Rijsdijk ; Purcell ; Raffa and Thompson 2014). Here we briefly consider the expected value of two test statistics that have been used for hypothesis testing in variance component estimation, the Wald test, and the likelihood-ratio-test statistic.

The Wald test is based on /var(), which under the null hypothesis that h2 = 0, follows a χ2 distribution. However, if h2 > 0, the Wald test statistic follows approximately a noncentral χ2 with noncentrality parameter (NCPW)If the estimation of phenotypic variance is ignored, thenAlternatively, the null hypothesis that h2 = 0 can be tested with a likelihood-ratio test. Blangero and colleagues (Blangero ) presented a very simple equation for the ELRT statistic to test the null hypothesis of h2 = 0,Equation 11 converges to Equation 10 when h2(λi − 1) → 0. For pairs of relatives with relationship r, NCPW = Nr2h4 / (1 – r2h4) and NCPLRT = −N ln(1 − r2h4). These expressions are equivalent when r2h4 → 0. When the true parameter is far from the one being tested under the null, these expressions can give quite different values. Raffa and Thompson (2014) give an analysis based on asymmetrical confidence intervals for the heritability.

Numerical Examples

Figure 1 shows the approximation to the standard error of an estimate of heritability as a function of the population value, experimental sample size, and design. Four different designs were used: a pedigree design of unrelated full-sib pairs, a pedigree design with MZ and DZ twins pairs with a ratio of 1:2 MZ and DZ pairs, a within-family design using full-sib pairs, and a population design using nominally unrelated individuals. In the last two designs, GRM are estimated with SNP data. These designs are less powerful than the pedigree-based experimental designs, but make fewer assumptions. At N = 10,000 the sampling variance of the population design approaches that of the pedigree designs, and at N = 100,000 it becomes the most powerful design. Sample sizes of 100,000 are realistic in human population and even larger samples sizes are expected in the next few years. Therefore, strong inference on heritability can be drawn using random samples from the population, while not having to make assumptions about the resemblance between relatives due to common environmental factors. The within-family design, which is the most robust with respect to assumptions of the model, remains inaccurate even when the analysis is on 50,000 full-sib pairs. However, in species such as fish with huge full-sib family sizes, accurate estimation could be achieved (Odegard and Meuwissen 2012; Hill 2013).

Figure 1

Standard error of estimates of heritability from different experimental designs in human populations, as a function of the population value of the heritability (x-axis), experimental sample size, and experimental design. For the within-family design (Within-family estimation using realized relationships estimates from markers), the variance in realized relationships was assumed to be 0.0392. For the population design (Random sampling from the population), the variance is relatedness was approximated assuming Ne = 10,000, a genome length of 35 M, and an average chromosome length of 1 M (Goddard 2009).

Figure 2 shows results for the population design for species with different Ne values of 1000, 10,000, and 100,000. It shows the increase in sampling variation with increasing effective population size, which is due to the decrease in the variation in relatedness. For the within-family design the sampling variance of heritability does not depend on the effective population size.

Figure 2

Standard error of the estimate of heritability from random samples of individuals from populations with different effective size and SNP-derived relationship matrices. For each population, a genome length of 35 M and an average chromosome length of 1 M was assumed (Goddard 2009).

Discussion

We have presented a general framework to quantify the sampling variance of heritability as a function of its population value, the sample size (N), and experimental design. Figure 1 shows that the sampling variance is relatively insensitive to the true value of h2 except when h2 → 1. The results recapitulate results from balanced designs and show that for pedigree designs, the sampling variance tends to be proportional to 1/N. In contrast, for designs that use genetic markers to estimate relatedness within families or estimate relatedness among randomly sampled individuals, the sampling variance is proportional to 1/N2. Consequently, very large samples of “unrelated” individuals are powerful for estimating h2. The key feature of the experimental design is the variation in relatedness. This is small within families of full-sibs and consequently the sampling variance of h2 is large.

There are a number of limitations to our study. First, we have assumed that the parameter whose sampling variance we derive is the same in different experimental designs. Even in the absence of confounding factors such as common environmental effect or nonadditive genetic factors, this is not necessarily the case. For the pedigree and within-family design, the parameter given our model assumptions is the narrow-sense heritability. But for the population design it is the proportion of phenotypic variance captured by genetic markers. If these markers are not sufficiently correlated with the genetic variants that cumulatively contribute to the total narrow sense heritability, then the use of a marker-based GRM will estimate additive genetic variation that is less than the total additive genetic variance. This can occur if the properties of the markers used to create the GRM are different from the segregating causal variants, for example, if the GRM is based upon common SNPs and the causal variants have lower heterozygosity, leading to loss of information due to imperfect linkage disequilibrium (Yang ). Although a “marker heritability” is conditional on the markers used to estimate relatedness, it is a valid population parameter with predictable sampling properties (as shown in this study). In human populations, it has been used to address the question of “missing heritability” from genome-wide association studies (Yang ).

Second, we assume that all resemblance between relatives is due to additive genetic covariance, so that there are only two random effects in the model. Additional random effects, for example, common environmental effects, make the covariance matrix V more complicated and generally not diagonalizable. When there are additional variance components, the residual variance as used in this study is partitioned in two or more components. These additional components are also estimated with error and will have a sampling covariance with the estimate of heritability. We suspect that having additional variance components in the model will tend to increase the sampling variance of the heritability, except for some balanced designs. However, we have not investigated general properties for designs with multiple random effects. With more than two variance components, computer simulation might be an efficient way to quantify the sampling variance of heritability and the proportion of variance due to additional random effects.

A third assumption is that estimation is by maximum likelihood or, alternatively, that fixed effects and covariates have been adjusted for without error. In practice, researchers tend to use least squares for balanced designs and restricted maximum likelihood (REML) or Bayesian methods for unbalanced designs. The difference in sampling variance between ML and REML is small when there are few fixed effects relative to the sample size, as, for example, in human genetic applications, but larger in situations where there are many fixed effects (e.g., in livestock applications).

Recently, Raffa and Thompson (2014) extended the work of Blangero by deriving approximations to the ELRT and confidence intervals of the heritability estimate using Taylor series expansions of the expected likelihood-ratio test with respect to the distribution of the eigenvalues of a given pedigree. Their simplest approximation can be expressed as an approximate sampling variance of the estimate of heritability as 2/[(N − 1)var(λ)] ≈ 2/(N var(λ)). This expression is the same as our special cases and All . The authors show that this approximation is not accurate when the assumptions break down, in particular when eigenvalues are not closely distributed around the mean of 1, and provide a better approximation using the logarithm of the eigenvalues (Raffa and Thompson 2014). They also show that confidence intervals of the estimates of heritability are not symmetrical when the variance in eigenvalues is large and that Wald statistic-based confidence intervals can be too narrow, implying that the use of the derived standard errors in our study to construct a confidence interval can be anticonservative. Although the derivations from Raffa and Thompson were for a pedigree design, they should also apply to other experimental designs, such as those where GRMs are estimated from marker data.

In conclusion, we have proposed a general unified framework to assess the sampling variance of the estimate of heritability using pedigree or marker-based relationships and have quantified how the sampling variance depends on sample size and the variation in relatedness.