Heritability from Mixed-Model Variance Components

What a Mixed-Model Heritability Estimate Measures

A mixed model separates variation associated with genomic resemblance among individuals from variation left in the residual component. Heritability summarizes the relative size of those variance components. The calculation is simple once the model and the scale of its relationship matrix are clear; the interpretation depends on both.

For a quantitative trait, consider the model

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

Here \(\mathbf{y}\) is the vector of phenotypes, \(\mathbf{X}\boldsymbol{\beta}\) is the fixed-effects mean, \(\mathbf{g}\) is a random genetic effect, and \(\mathbf{e}\) is residual variation. The matrix \(\mathbf{K}\) is a genomic relationship matrix. Its entries describe relative genetic resemblance, and its scaling determines how \(\sigma_g^2\) translates into genetic variance. The variance components are usually estimated by restricted maximum likelihood (REML), as in the earlier material on REML estimation of variance components.

Definition: When \(\mathbf{K}\) is scaled so that its average diagonal entry is one, the model-based heritability is \(h^2_{\mathrm{model}}=\sigma_g^2/(\sigma_g^2+\sigma_e^2)\). It is the fraction of the modeled random-component variance assigned to \(\mathbf{g}\), on the phenotype scale and for the specified model and relationship matrix.

The assumptions matter. The genetic effect has covariance \(\sigma_g^2\mathbf{K}\); the residuals have common variance \(\sigma_e^2\) and no residual covariance; and the genetic and residual effects are independent. The relationship matrix must be positive semidefinite and appropriately constructed for the genetic component being modeled. Fixed effects describe the conditional mean, not another random variance component in this ratio. If the phenotype’s variance is summarized across a population that varies in fixed covariates, that total variance can include additional variation not represented in the denominator above.

Deriving the Ratio and Accounting for Matrix Scaling

The model’s covariance is \(\operatorname{Var}(\mathbf{y}\mid\mathbf{X})=\sigma_g^2\mathbf{K}+\sigma_e^2\mathbf{I}\). For individual \(i\), the diagonal entry is \(\sigma_g^2K_{ii}+\sigma_e^2\). If the average diagonal of \(\mathbf{K}\) is one, the average marginal variance across individuals is therefore \(\sigma_g^2+\sigma_e^2\). The modeled genetic share is the genetic component divided by this total.

\frac{\text{average genetic variance}}{\text{average total random-component variance}} = \frac{\sigma_g^2\left(\frac{1}{n}\sum_{i=1}^n K_{ii}\right)} {\sigma_g^2\left(\frac{1}{n}\sum_{i=1}^n K_{ii}\right)+\sigma_e^2} \quad\longrightarrow\quad \frac{\sigma_g^2}{\sigma_g^2+\sigma_e^2} \quad\text{if}\quad \frac{1}{n}\sum_{i=1}^n K_{ii}=1.

This also shows why matrix scaling is not a cosmetic choice. If the average diagonal is \(d\ne1\), the average genetic variance is \(d\sigma_g^2\), not \(\sigma_g^2\). The corresponding average-variance ratio is \(d\sigma_g^2/(d\sigma_g^2+\sigma_e^2)\). Alternatively, rescaling \(\mathbf{K}\) changes the numerical estimate of \(\sigma_g^2\); the product that determines genetic variance must be interpreted on the chosen scale.

Justification: The ratio follows from the variance decomposition in the model: independent genetic and residual components have variances that add. It is not a proportion of the observed phenotype explained by every possible genetic mechanism; it is a proportion of the model’s specified random-component variance.

Worked Example: Heritability With a Standard-Scaled Relationship Matrix

Suppose a hypothetical height analysis uses a relationship matrix with average diagonal one. REML estimates the genetic variance as \(18\ \text{cm}^2\) and the residual variance as \(32\ \text{cm}^2\). Assume the fitted model’s covariance and residual assumptions are reasonable for this illustration.

The total random-component variance estimate is \(18+32=50\ \text{cm}^2\). Thus

\hat h^2_{\mathrm{model}} = \frac{\hat\sigma_g^2}{\hat\sigma_g^2+\hat\sigma_e^2} = \frac{18\ \text{cm}^2}{18\ \text{cm}^2+32\ \text{cm}^2} = \frac{18}{50} = 0.360.

The estimated modeled genetic component accounts for 36.0% of the estimated genetic-plus-residual variance in this model. The units cancel, so heritability is unitless. This estimate does not say that 36.0% of an individual’s height is genetic, nor does it imply that changing genetic makeup would cause a particular height change.

Worked Example: Correcting for Relationship-Matrix Scale

Consider a hypothetical trait measured in standardized units. Suppose \(\hat\sigma_g^2=0.40\) and \(\hat\sigma_e^2=0.50\), while the relationship matrix has average diagonal \(d=1.25\). The variance attributable to the genetic component averaged over individuals is \(1.25(0.40)=0.50\) standardized-units squared.

The scale-adjusted estimate is

\hat h^2_{\mathrm{model}} = \frac{d\hat\sigma_g^2}{d\hat\sigma_g^2+\hat\sigma_e^2} = \frac{1.25(0.40)}{1.25(0.40)+0.50} = \frac{0.50}{1.00} = 0.500.

Using \(0.40/(0.40+0.50)=0.444\) would incorrectly treat the genetic variance parameter as the average genetic variance despite the matrix’s diagonal scale. The example assumes the residual variance is on the same phenotype scale and the average diagonal is the relevant scaling summary. In practice, inspect how \(\mathbf{K}\) was constructed and scaled before calculating or reporting a ratio.

Relation to SNP Heritability

“SNP heritability” usually refers to the proportion of phenotypic variance attributable to additive effects of a defined set of measured or imputed SNPs, under a specified genetic model and phenotype scale. A genomic relationship matrix built from those SNPs can encode their aggregate similarity. In that setting, a variance-component estimate based on \(\mathbf{K}\) is often an estimate of the additive variance tagged by that marker set, and the resulting ratio can be called SNP heritability.

The two terms are not automatically interchangeable. A variance-component ratio describes the genetic component represented by the particular matrix in the fitted model. Its target depends on the variants used to construct \(\mathbf{K}\), their weighting and scaling, the sample, and the assumptions about effects and covariance. It does not automatically include additive effects poorly tagged by those markers, rare or unmeasured variants, or non-additive genetic effects. Conversely, if the relationship matrix is constructed from a particular SNP set and the assumptions align, the model-based ratio may be the method used to estimate that set’s SNP heritability.

Heritability is also population- and environment-dependent: the ratio can change when allele frequencies, linkage disequilibrium, environmental variation, or the analyzed population changes. It is not a fixed percentage of a trait that applies to every population. For binary traits, the scale requires special care: a variance ratio calculated on an observed zero-one scale is not the same as a ratio on an assumed latent liability scale.

Worked Example: Distinguishing Two Heritability Targets

Suppose a hypothetical analysis of a quantitative trait uses a matrix constructed from common genotyped SNPs and estimates \(\hat\sigma_g^2=12\) and \(\hat\sigma_e^2=28\), with average diagonal one. A separate hypothetical pedigree analysis reports a heritability estimate of \(0.42\). The estimates differ, but they need not be estimates of exactly the same quantity.

For the genomic relationship model,

\hat h^2_{\mathrm{model}} = \frac{12}{12+28} = \frac{12}{40} = 0.300.

This is the estimated share assigned to the common-SNP relationship component in this model. The pedigree estimate may reflect a broader genetic similarity structure and relies on different assumptions about how resemblance arises. A difference of \(0.120\) between the reported numbers does not, by itself, show that either method is wrong. Before comparing them, specify the population, phenotype scale, genetic effects represented, and variance denominator for each estimate. If both analyses were designed to target the same genetic component on the same scale, then their uncertainty and modeling assumptions would also need to be considered.

Diagnostics, Boundaries, and Limitations

A calculated ratio is only as informative as its variance estimates and model. Review the relationship matrix’s construction, scaling, and eigenstructure; verify that phenotype and genotype records align to the same individuals; and examine whether residual independence and common variance are plausible. Systematic residual patterns or unmodeled group structure can make the decomposition misleading. A ratio within the interval from zero to one is expected when both variance estimates are nonnegative and the denominator is positive. An estimate outside that interval usually signals a calculation, scaling, or model-parameterization problem.

Variance estimates can lie on the boundary. If REML gives \(\hat\sigma_g^2=0\) and \(\hat\sigma_e^2>0\), the plug-in ratio is \(0\). This is a boundary estimate, not proof that the population genetic variance is exactly zero. It can reflect limited information, an unsuitable relationship matrix, or genuine weak signal; uncertainty and model diagnostics are needed for interpretation. If both estimated components are zero, the ratio is undefined because its denominator is zero.

The ratio is an estimate, not a known constant. Uncertainty in both variance components carries into uncertainty in heritability, so a point estimate alone can conceal substantial imprecision. REML variance estimates, their covariance, and an appropriate interval or other uncertainty summary should be considered when inferential precision matters. Also, do not interpret the model-based ratio as a causal effect or as the percentage of a trait that could be changed by intervention.

Alternative approaches may target different quantities or use different information. Pedigree methods infer genetic variance from resemblance among relatives; SNP-based methods use marker-derived relatedness to estimate variance tagged by a specified marker set; and single-variant association results describe effects at individual loci rather than the aggregate variance ratio. A meaningful comparison requires matching genetic target, scale, population, and assumptions rather than comparing labels alone.

Common Mistakes

Ignoring matrix scaling. The incorrect reasoning is that \(\sigma_g^2\) always equals average genetic variance. It does so under the relevant normalization; otherwise use the matrix scale, such as \(d\sigma_g^2\) for average diagonal \(d\).
Calling every ratio SNP heritability. The correction is to name the genetic component represented by the fitted relationship matrix and the marker set used to build it.
Reading the estimate as an individual-level statement. A ratio summarizes variance in a specified population and model; it is not the fraction of one person’s phenotype caused by genes.
Treating a boundary estimate as proof of no genetic contribution. A zero estimate is a fitted boundary value; it does not establish that the underlying population parameter is exactly zero.

Key takeaway: Calculate model-based heritability as genetic variance divided by total modeled random-component variance, using the relationship matrix’s actual scale. Interpret the result as specific to the genetic component, population, phenotype scale, and assumptions of the fitted model.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A model uses a relationship matrix with average diagonal one and estimates genetic variance \(9\ \text{units}^2\) and residual variance \(21\ \text{units}^2\). What is the estimated model-based heritability?
A matrix has average diagonal \(0.80\), with \(\hat\sigma_g^2=0.50\) and \(\hat\sigma_e^2=0.40\). Calculate the average-variance ratio.
Why does a model-based estimate from a common-SNP relationship matrix not necessarily equal heritability from a pedigree study?
A fitted model returns \(\hat\sigma_g^2=0\) and \(\hat\sigma_e^2=10\). What is the plug-in ratio, and what does it not establish?

Solutions. (1) \(9/(9+21)=9/30=0.300\), or 30.0% of modeled random-component variance. (2) The average genetic variance is \(0.80(0.50)=0.40\), so the ratio is \(0.40/(0.40+0.40)=0.500\). (3) The methods can represent different genetic components and depend on different assumptions, samples, and scales; matching their labels does not make their targets identical. (4) The plug-in ratio is \(0/(0+10)=0\). This boundary estimate does not prove that the population genetic variance is exactly zero.