GREML Estimation of SNP Heritability from Individual Genotypes

Why GREML Uses Genomic Similarity

A SNP set can be informative about a trait even when no single variant explains much of its variation. GREML (genomic-relatedness-based restricted maximum likelihood) uses individual-level genotypes to estimate the additive variance tagged by a specified set of SNPs. Its key idea is that, under an additive genetic model, individuals who are more genetically similar at those SNPs tend to have more similar genetic values. The method estimates how much phenotypic covariance follows that pattern.

The target is therefore tied to the markers used to construct the genomic relationship matrix (GRM). It is not automatically the trait’s total genetic variance: genetic effects not tagged by the chosen SNP set need not contribute to the estimate.

The GRM and the GREML Model

Suppose \(n\) individuals have genotypes at \(m\) SNPs in the specified set. Let \(x_{ij}\) be the allele dosage for individual \(i\) at SNP \(j\), and let \(p_j\) be the allele frequency used for centering and scaling. A common construction standardizes each dosage and averages the resulting cross-products:

z_{ij}=\frac{x_{ij}-2p_j}{\sqrt{2p_j(1-p_j)}}, \qquad K=\frac{1}{m}ZZ^{\mathsf T}.

Here \(Z\) is the \(n\times m\) matrix of standardized dosages, and \(K\) is an \(n\times n\) matrix. The entry \(K_{ik}\) summarizes genetic similarity between individuals \(i\) and \(k\) across the SNP set. The diagonal entries describe each individual’s genetic similarity to themselves; they are often near, but not necessarily exactly, one. Centering, scaling, and any genotype handling should be applied consistently with the chosen GRM convention.

GREML relates this matrix to the covariance of a quantitative phenotype. Let \(y\) be the \(n\)-vector of observed phenotypes and \(X\) the design matrix for fixed effects, such as an intercept and prespecified covariates. The standard single-component model is:

y=X\beta+g+e, \qquad g\sim N(0,\sigma_g^2K), \qquad e\sim N(0,\sigma_e^2I).

The genetic random effect \(g\) represents additive effects tagged by the SNP set. The residual \(e\) includes remaining variation not represented by that component, as well as individual-level noise. In this model, \(g\) and \(e\) are independent, and the residuals have equal variance and no residual covariance between different individuals. Conditional on the fixed effects, the phenotype covariance matrix is therefore:

V=\operatorname{Var}(y\mid X)=\sigma_g^2K+\sigma_e^2I.

Definition: In the single-component GREML model, \(\sigma_g^2\) is the additive genetic variance attributed to the specified SNP set under the model, and \(\sigma_e^2\) is the residual variance. The model-based SNP heritability is \(h^2_{\mathrm{SNP}}=\sigma_g^2/(\sigma_g^2+\sigma_e^2)\), when the GRM scaling and variance components are defined consistently.

Worked Example: Computing a Small GRM

Consider a hypothetical set of three individuals and two SNPs. Suppose the standardized genotype matrix is:

Z= \begin{pmatrix} 1&0\\ 0&1\\ -1&-1 \end{pmatrix}.

Using \(K=ZZ^{\mathsf T}/m\) with \(m=2\), first calculate the cross-product matrix:

ZZ^{\mathsf T}= \begin{pmatrix} 1&0&-1\\ 0&1&-1\\ -1&-1&2 \end{pmatrix}, \qquad K= \begin{pmatrix} 0.5&0&-0.5\\ 0&0.5&-0.5\\ -0.5&-0.5&1 \end{pmatrix}.

For example, \(K_{13}=(1\cdot -1+0\cdot -1)/2=-0.5\), while \(K_{12}=(1\cdot0+0\cdot1)/2=0\). These entries are similarity summaries under this standardization; they are not probabilities or direct measures of a causal effect. The matrix would enter \(V\) in the GREML model, so its off-diagonal pattern determines the model’s genetic covariance between distinct individuals.

How REML Estimates the Variance Components

The unknown quantities of interest are \(\sigma_g^2\) and \(\sigma_e^2\), constrained to be nonnegative. Restricted maximum likelihood (REML) estimates them while accounting for the fixed effects in \(X\). Conceptually, REML evaluates the covariance of phenotype contrasts that remove those fixed effects, rather than treating the estimated fixed effects as if they were known.

For a precise expression, let \(A\) be an \(n\times(n-q)\) matrix whose columns are orthonormal and satisfy \(A^{\mathsf T}X=0\), where \(q\) is the rank of \(X\). The transformed data \(A^{\mathsf T}y\) no longer contain the fixed-effect mean. Its covariance is \(A^{\mathsf T}VA\). Up to the displayed normalizing constant, the restricted log-likelihood is:

\ell_R(\sigma_g^2,\sigma_e^2) =-\frac{1}{2}\left[ (n-q)\log(2\pi) +\log|A^{\mathsf T}VA| +y^{\mathsf T}A(A^{\mathsf T}VA)^{-1}A^{\mathsf T}y \right].

REML chooses nonnegative variance components that maximize this expression. In practice, software uses numerical optimization and may use equivalent formulations for computation. The result depends on the phenotype, fixed-effect specification, GRM, and model assumptions; it is not obtained by simply averaging pairwise phenotypic similarities.

Why the covariance pattern identifies the component: For distinct individuals \(i\ne k\), the model gives \(\operatorname{Cov}(y_i,y_k\mid X)=\sigma_g^2K_{ik}\), because the independent residual term contributes no off-diagonal covariance. On the diagonal, \(\operatorname{Var}(y_i\mid X)=\sigma_g^2K_{ii}+\sigma_e^2\). Thus the model attributes covariance that tracks \(K\) to the genetic component, while the residual component contributes to individual variance.

Worked Example: Turning REML Estimates Into SNP Heritability

Suppose a hypothetical GREML analysis of a quantitative trait uses a specified SNP set, an intercept and age as fixed effects, and returns \(\hat{\sigma}_g^2=0.36\) trait-units squared and \(\hat{\sigma}_e^2=0.64\) trait-units squared. The fitted total variance is \(0.36+0.64=1.00\) trait-units squared. The plug-in heritability estimate is:

\hat{h}^2_{\mathrm{SNP}} =\frac{\hat{\sigma}_g^2}{\hat{\sigma}_g^2+\hat{\sigma}_e^2} =\frac{0.36}{0.36+0.64} =\frac{0.36}{1.00} =0.36.

Under this model, the specified SNP set accounts for an estimated 36% of the modeled phenotypic variance in the analyzed population and trait scale. This is a variance proportion, so it has no physical units. It is not a statement that 36% of any individual’s trait value is genetic, nor does it establish that the tagged variants cause the trait.

Why the Ratio Is the Relevant Summary

The genetic variance component is measured in squared trait units, just like the residual variance. Dividing by their sum expresses the genetic component relative to the total variance represented by the model. The ratio is dimensionless and, when the components are nonnegative and the denominator is positive, lies between zero and one.

This ratio relies on the model’s variance decomposition. It does not imply that the phenotype is literally divided into a genetic part and a residual part for each person. It summarizes population variation under the fitted model. Changing the population, phenotype scale, fixed effects, or SNP set can change the estimate.

Worked Example: Comparing Estimates From Two SNP Sets

Consider two hypothetical GREML fits to the same trait in the same cohort, with the same fixed effects and compatible phenotype scaling, but with different SNP sets and corresponding GRMs. For Set A, suppose the fitted components are \(\hat{\sigma}_{g,A}^2=0.24\) and \(\hat{\sigma}_{e,A}^2=0.76\). For Set B, suppose they are \(\hat{\sigma}_{g,B}^2=0.31\) and \(\hat{\sigma}_{e,B}^2=0.69\). The calculations are:

\hat{h}^2_{\mathrm{SNP},A} =\frac{0.24}{0.24+0.76} =0.24, \qquad \hat{h}^2_{\mathrm{SNP},B} =\frac{0.31}{0.31+0.69} =0.31.

In these hypothetical fits, Set B has a larger estimated tagged variance proportion: 31% rather than 24%. This comparison does not show that each additional SNP in Set B has an independent effect or that the difference is statistically reliable. The sets define different estimands and can differ in how well they capture additive genetic similarity. A comparison should therefore report the SNP sets and modeling choices, not only the two ratios.

Assumptions and Interpretation

A GREML estimate is meaningful to the extent that its covariance model is appropriate. Important assumptions and checks include:

Appropriate GRM: The genotype data, SNP selection, allele-frequency estimates, and scaling must produce a GRM that represents the intended SNP set. The variance-component scale depends on the GRM convention.
Fixed effects adequately address mean structure: The columns of \(X\) should include relevant design variables and prespecified covariates. Unmodeled differences in means can distort covariance-based inference.
Residual covariance is adequately modeled: The simple model assumes residual variance \(\sigma_e^2\) for every individual and no residual covariance across individuals. Violations can be mistaken for, or obscure, genetic covariance.
Genetic and residual components are independent: The model assigns distinct covariance structures to \(g\) and \(e\). If the assumed decomposition is poor, the fitted components may not have their intended interpretation.

These assumptions do not make the result a universal property of a trait. GREML estimates variance tagged by the chosen markers in the analyzed population, under the specified model. It does not, by itself, separate direct genetic effects from every source of genetic similarity or establish causality.

Common Mistakes

Wrong: “The estimate is the fraction of an individual’s trait caused by these SNPs.” Correction: It is a modeled proportion of phenotypic variance across individuals, not an individual-level causal fraction.
Wrong: “A GREML estimate is total heritability.” Correction: It concerns additive variance tagged by the specified SNP set, under the chosen GRM and model.
Wrong: “A larger estimate proves that one SNP set has more causal variants.” Correction: Set differences can reflect marker coverage and the covariance structure they capture; the estimate alone does not identify causal variants.
Wrong: “The GRM diagonal is always exactly one.” Correction: Its diagonal depends on the data and scaling convention and is often, but not necessarily, near one.

Key takeaway: GREML estimates SNP-set additive variance by fitting a covariance pattern proportional to a genotype-based GRM. Its heritability estimate is the fitted genetic variance divided by fitted total variance, interpreted for the specified SNP set, trait scale, population, and model.

Check Your Understanding

Answer each question using the model and interpretation developed above.

In the single-component GREML model, what role does the GRM play in the phenotype covariance matrix?
If the fitted genetic and residual variances are 0.18 and 0.42 trait-units squared, calculate the SNP heritability estimate.
Why can the diagonal entries of a GRM differ from one?
What does an estimated SNP heritability of 0.30 describe, and what does it not establish?
Why might two different SNP sets yield different GREML estimates for the same trait and cohort?