Estimating Dominance Variance From Genomic Data

Why Estimate Dominance Variance?

A genetic model that includes only additive effects can miss systematic differences between heterozygotes and the midpoint of the two homozygote means. A dominance variance component quantifies the variation associated with those nonadditive, within-locus contrasts, as represented by the SNPs and model being analyzed. Estimating it is not simply a matter of counting heterozygotes: it requires a covariance model for individuals, a clearly stated target, and an assessment of uncertainty.

The target is also narrower than “all biological dominance.” A SNP-based estimate describes the dominance variation captured by the chosen marker set and its covariance model in the population and phenotype scale under study. It need not include untyped causal variants, and it can be difficult to distinguish from other sources of variation when their covariance patterns are similar.

From a SNP Contrast to a Variance Component

Let \(g\) denote the number of copies of a specified allele at a biallelic SNP, so \(g\in\{0,1,2\}\). Let \(p\) be that allele's frequency and \(q=1-p\). Under Hardy–Weinberg equilibrium, the genotype probabilities are \(q^2\), \(2pq\), and \(p^2\). An orthogonal dominance contrast under this model assigns values \(-2p^2\), \(2pq\), and \(-2q^2\) to genotypes \(0\), \(1\), and \(2\), respectively. Its mean is zero, its variance is \(4p^2q^2\), and it is uncorrelated with the centered additive genotype under the stated equilibrium assumptions. Dividing by \(2pq\) gives a contrast with variance one when \(0<p<1\).

These facts help define the genomic signal, but an estimated dominance variance is usually a component in a model of the joint outcomes of many individuals, rather than a separate estimate at each SNP added together without qualification. The model uses a dominance covariance structure derived from SNP genotypes. For the purposes of estimation, denote that structure by \(K_D\); it is treated as a known, positive semidefinite matrix whose average diagonal is normalized to one. This normalization puts the dominance variance component on the phenotype-variance scale.

The value of \(\sigma_D^2\) depends on the marker set, coding and normalization, population, phenotype scale, and other components included in the model. It is not automatically the variance due to every dominance effect in the genome. In particular, a component attributed to measured SNPs may capture tagged effects as well as effects directly associated with those SNPs.

A Mixed Model for Estimation

For individual-level phenotypes, a common model separates fixed covariates from random genetic and residual effects:

Here \(y\) is the vector of observed phenotypes, \(Xb\) represents fixed effects such as an intercept and specified covariates, \(u_A\) is the additive genetic effect, and \(u_D\) is the modeled dominance effect. The matrices \(K_A\) and \(K_D\) are the corresponding covariance structures; \(I\) is the identity matrix, and \(\sigma_e^2\) is residual variance. The target of this tutorial is \(\sigma_D^2\), not the construction of the matrices.

This model is used with explicit assumptions. The fixed-effects design is correctly specified for the covariates included; the covariance matrices are treated as known for the analysis; and the random effects are mean-zero multivariate normal with the stated covariances. The usual component model also assumes that \(u_A\), \(u_D\), and \(e\) are mutually independent, conditional on the model, so that their covariance contributions add. If there are other relevant sources of covariance, omitting them can affect the fitted components. Normality is a working likelihood assumption for standard likelihood-based inference; it should not be mistaken for proof that the biological effects themselves are normally distributed.

The phenotype covariance implied by the model is

Restricted maximum likelihood (REML) estimates variance components while accounting for the fixed effects. Under the model, it finds component values that make the observed phenotype covariance plausible after accounting for the fitted fixed-effect space. Estimation is constrained so that variance components cannot be negative. Consequently, a fitted dominance component may lie at the boundary, \(\hat{\sigma}_D^2=0\).

Standard Errors and Uncertainty

A variance-component estimate is incomplete without an uncertainty assessment. For an interior REML estimate, a model-based standard error can be obtained from the estimated information matrix: in a regular setting, the approximate covariance matrix of the component estimates is the inverse of the information matrix, and the standard error for a component is the square root of its corresponding diagonal entry. Equivalently, local likelihood curvature indicates how sharply the data distinguish nearby component values. A flatter likelihood implies less precise estimation.

This approximation needs care. Variance components are constrained to be nonnegative, and the null value is on the boundary. Near zero, the sampling distribution can be skewed and need not be well approximated by a symmetric normal distribution. A Wald interval of the form estimate plus or minus a critical value times its standard error can then extend below zero and misrepresent uncertainty. Profile-likelihood intervals or a parametric bootstrap can better reflect the constraint, although their validity still depends on the model and fitting procedure.

A proportion such as \(\sigma_D^2/(\sigma_A^2+\sigma_D^2+\sigma_e^2)\) is also estimated, not known. Its uncertainty depends on the uncertainty in every component and on their covariances. For an interior estimate, the delta method approximates that uncertainty by applying the gradient of the proportion to the estimated covariance matrix of the component estimates. Treating the component estimates as independent without justification can produce an incorrect standard error.

Why Dominance Estimates Can Be Imprecise

The precision of \(\hat{\sigma}_D^2\) depends on whether the data contain covariance patterns that distinguish \(K_D\) from the other modeled components. If \(K_D\) is similar to \(K_A\), for example, changes in their component values can produce similar fitted covariance matrices. The likelihood then has little curvature in some combinations of component values. Increasing the sample size can help, but it does not automatically solve weak separation between covariance structures.

Allele frequency also affects the information in dominance contrasts. Under Hardy–Weinberg equilibrium, the unstandardized contrast has variance \(4p^2q^2\), which is small when one allele is rare. Standardizing the contrast to unit variance rescales its values, but does not create observations of rare genotype classes that are absent or scarce. Estimates can consequently be sensitive to allele-frequency estimation, genotype counts, marker selection, and the assumed population structure.

Finally, measurement noise, phenotype scale, and misspecified covariance can all affect a fitted component. A dominance estimate should therefore be reported with its scale, normalization, uncertainty method, and the modeled covariates and variance components. A small estimate with a large standard error is not strong evidence that the true component is exactly zero; it may indicate that the available data do not distinguish it precisely from other components.

A Short Justification: Information Depends on Covariance Patterns

The role of the covariance structure can be made explicit. Let \(X\) be the fixed-effects design and \(V\) the phenotype covariance matrix. Define

For a Gaussian REML model with known covariance matrices and an interior parameter value, the expected information for variance components \(\theta_i\) and \(\theta_j\), whose covariance matrices are \(K_i\) and \(K_j\), is

To see why this expression matters, consider the REML log likelihood, up to constants,

Differentiating with respect to \(\theta_i\), for which \(\partial V/\partial\theta_i=K_i\), gives the score

Taking the expected negative derivative of this score with respect to \(\theta_j\), under the fitted Gaussian covariance model, yields \(\mathcal{I}_{ij}=\frac{1}{2}\operatorname{tr}(PK_iPK_j)\). Thus information about the dominance component depends on how its covariance pattern remains distinguishable after accounting for fixed effects and the other covariance components. Similar patterns make the information matrix difficult to invert precisely; consequently, the component standard error can be large.

Common Mistakes and Practical Interpretation

• Calling the estimate total dominance variance. The component pertains to the modeled SNPs, coding, normalization, population, and phenotype scale. It need not equal variance from every causal dominance effect.
• Interpreting a variance proportion as an individual-level causal fraction. A fitted share describes a model-based population variance decomposition; it does not describe the fraction of one person's phenotype caused by dominance.
• Reporting the estimate without its uncertainty. A point estimate alone does not show whether the data distinguish dominance variance from other components. Report an appropriate standard error or interval and explain the method.
• Using a symmetric interval mechanically at zero. A negative lower limit for a nonnegative variance component is a warning that the normal approximation may be unsuitable. Consider a boundary-aware method.
• Assuming standardization guarantees precision. Unit-variance coding sets a scale; it cannot supply rare genotype classes or separate highly similar covariance patterns.