REML Uncertainty and Boundary Estimates for SNP Heritability

Why Uncertainty and Boundaries Matter

A GREML estimate of SNP heritability is a ratio built from estimated genetic and residual variance components. Its uncertainty therefore depends on both component estimates and on their covariance. It also depends on the shape of the likelihood: when a variance component is estimated as zero, the estimate lies at the edge of the allowed parameter space, where ordinary symmetric approximations can be misleading.

The ratio has a specific meaning only after the genetic relationship matrix has been scaled. In this tutorial, \(K\) is scaled so that its mean diagonal is one. This convention makes the average marginal genetic variance equal to the genetic variance component and gives the usual variance-ratio interpretation of SNP heritability.

GREML Model and Heritability Scale

The marginal covariance matrix of \(y\) is \(V=\sigma_g^2K+\sigma_e^2I\). In general, individual \(i\) has marginal variance \(\sigma_g^2K_{ii}+\sigma_e^2\), not necessarily \(\sigma_g^2+\sigma_e^2\). The mean of those marginal variances is \(\sigma_g^2\operatorname{tr}(K)/n+\sigma_e^2\), which equals \(\sigma_g^2+\sigma_e^2\) under the stated scaling. If a different scale for \(K\) is used, the variance ratio must be interpreted on that scale or the matrix must first be normalized.

REML estimates variance components using information in the data after accounting for the fixed effects. In standard notation, let \(p\) be the number of columns of \(X\), and define

Apart from constants, the restricted log likelihood is

REML chooses variance components to maximize this likelihood subject to nonnegativity. This constrained maximization is important: a variance estimate cannot be negative, so the maximizing value can be exactly zero. The corresponding heritability estimate can be zero or, if the residual variance is zero, one.

From REML Curvature to Uncertainty in the Ratio

Write \(\theta=(\sigma_g^2,\sigma_e^2)'\), and let \(V_a\) denote the derivative of \(V\) with respect to component \(\theta_a\). For the model above, these derivatives are \(K\) and \(I\). The expected REML information for the variance components has entries

This information describes the local curvature of the restricted log likelihood. At an interior maximum, the inverse of the information, or the inverse observed information evaluated at the estimates, is commonly used as an approximate covariance matrix for \(\hat{\theta}\). Thus the diagonal entries give approximate variances for the component estimates, and the off-diagonal entry records their dependence. This is an asymptotic, model-based approximation, not a guarantee of accurate uncertainty in every finite sample.

To transfer that uncertainty to heritability, use the delta method: for a smooth function of estimates, approximate its variance by the gradient of the function multiplied on both sides of the covariance matrix. For \(h^2=a/(a+b)\), where \(a=\sigma_g^2\) and \(b=\sigma_e^2\), the gradient is

The first derivative is positive because increasing genetic variance while holding residual variance fixed increases the ratio. The second is negative because increasing residual variance while holding genetic variance fixed decreases it. If \(\widehat{\operatorname{Cov}}(\hat{\theta})\) is the estimated covariance matrix, then

The off-diagonal covariance matters: ignoring it can either overstate or understate the uncertainty in the ratio. The delta method also relies on the estimate being sufficiently far from a boundary for a local, approximately normal description to be useful.

Wald and Profile-Likelihood Intervals

A Wald interval takes the estimate plus or minus a normal critical value times its standard error. It is convenient and can work well when the likelihood is approximately quadratic and the estimate is well inside the permitted range. But \(h^2\) is constrained to \([0,1]\). A Wald interval can extend below zero or above one, and simply truncating its endpoints does not generally restore its stated coverage.

A profile-likelihood interval uses the likelihood itself rather than only its local curvature. One way to form it is to fix a candidate \(h\) and maximize the restricted likelihood over the total variance \(\tau=\sigma_g^2+\sigma_e^2\), while setting \(\sigma_g^2=h\tau\) and \(\sigma_e^2=(1-h)\tau\). The resulting profile log likelihood is denoted \(\ell_R^p(h)\). If \(\ell_{\max}\) is the unrestricted constrained maximum, the likelihood-ratio statistic for a candidate value is

For an interior parameter value, the usual approximate 95% likelihood-ratio confidence set contains values with \(\Lambda(h)\leq 3.84\), the 95th percentile of a chi-square distribution with one degree of freedom. Equivalently, the profile log likelihood is at least \(\ell_{\max}-1.92\) for values in that set. A profile interval naturally respects the allowed range \(0\leq h\leq1\), although the cutoff’s usual calibration requires care at a boundary.

Why Estimates Can Be Zero or One

The variance components are restricted to be nonnegative. If the likelihood would continue to improve as \(\sigma_g^2\) decreases toward zero, the constrained maximum is at \(\hat{\sigma}_g^2=0\), giving \(\hat{h}^2=0\) when the residual estimate is positive. This does not prove that the true genetic variance is zero. It means that, in this sample and model, the likelihood is maximized at the zero boundary; the data may still be compatible with positive values, as a one-sided profile interval can show.

At the other endpoint, if \(\hat{\sigma}_e^2=0\) while \(\hat{\sigma}_g^2>0\), then \(\hat{h}^2=1\). This is also a boundary fit, and it should prompt scrutiny of the model, data, and optimization. A zero residual component can arise from limited information or model misspecification as well as from an exceptionally strong fitted genetic covariance structure. It is not, by itself, evidence that all phenotypic variation is genuinely attributable to the measured SNP set.

Practical Interpretation and Common Pitfalls

Uncertainty in GREML SNP heritability is uncertainty about a model-based variance ratio for a specified phenotype, sample, and variant set. It does not directly quantify uncertainty about an individual’s genetic contribution, establish causation, or show that the estimated value will transfer unchanged to another population. The interval also does not account automatically for every possible source of model misspecification.

A useful reporting practice is to give the point estimate, the uncertainty method, the interval, and whether the fit was at a boundary. If Wald and profile intervals differ substantially, that difference is informative: it suggests that local curvature alone may not adequately describe the likelihood over the relevant parameter range. Numerical convergence and the fitted variance components should also be checked, especially when a component is estimated as zero or nearly zero.

• Mistake: Treating a zero genetic variance estimate as proof that genetic variance is absent. Correction: It is a constrained maximum; uncertainty and power determine which positive values remain plausible.
• Mistake: Reporting a Wald interval that crosses zero or one and assuming truncating it fixes the problem. Correction: Boundary-aware profile methods or calibrated procedures may be more suitable.
• Mistake: Applying a chi-square cutoff directly to a profile log likelihood. Correction: Compare the likelihood-ratio statistic \(\Lambda(h)\) with its cutoff, or compare the profile log likelihood with \(\ell_{\max}\) minus half that cutoff.
• Mistake: Interpreting \(\sigma_g^2+\sigma_e^2\) as the average phenotype variance for any \(K\). Correction: Check the scale of \(K\); the mean marginal variance is \(\sigma_g^2\operatorname{tr}(K)/n+\sigma_e^2\).