Worked Haseman–Elston Estimates from Genotype Similarity

From Pairwise Data to a Heritability Estimate

Haseman–Elston regression turns information about pairs of individuals into an estimate of genetic variance. The intuition is that, under an additive genetic model, individuals with greater genomic similarity tend to have more similar phenotypes. The regression quantifies that relationship: its slope estimates the genetic covariance associated with a one-unit increase in genomic relationship.

The calculation is easy to misread. A fitted slope is not itself heritability, and an intercept is not automatically the residual variance. To interpret either coefficient, first specify the pairwise model and the scale of the phenotype and genomic relationship.

Definition: For each pair of distinct individuals $i$ and $j$, let $r_{ij}$ be the product of their adjusted phenotypes and let $K_{ij}$ be their genomic relationship. The pairwise Haseman–Elston regression is the ordinary least-squares regression of $r_{ij}$ on $K_{ij}$, with an intercept: $$ r_{ij}=a+bK_{ij}+\text{regression error}. $$ Here $a$ is the fitted intercept and $b$ is the fitted slope. Under the standard additive model and a genomic relationship matrix scaled to have diagonal entries near one, $b$ estimates additive genetic variance, $\sigma_g^2$. If the model's phenotypic variance is $\sigma_p^2=\sigma_g^2+\sigma_e^2$, SNP heritability is $h^2_{\mathrm{SNP}}=\sigma_g^2/\sigma_p^2$.

Model and Calculation

The model underlying this interpretation is $y=X\beta+g+e$, where $y$ is the phenotype vector, $X\beta$ represents fixed effects such as an intercept and specified covariates, $g$ is the additive genetic effect captured by the SNP set, and $e$ is residual variation. Assume $g\sim N(0,\sigma_g^2K)$, $e\sim N(0,\sigma_e^2I)$, and that $g$ and $e$ are independent. The pairwise products should be formed from phenotypes adjusted for the relevant fixed effects, rather than from unadjusted phenotypes when those effects could create covariance.

For distinct individuals, residual errors are uncorrelated in this model. Consequently, after fixed effects are handled, the model gives $E(r_{ij}\mid K)$ approximately equal to $\sigma_g^2K_{ij}$. The word “approximately” matters in a finite sample: the regression is a moment estimator, and observed pairwise products will not fall exactly on their model expectation. An intercept is included to allow for a baseline pairwise product in the fitted data; under the idealized model that baseline is zero. It is not the diagonal residual variance.

To calculate the coefficients, treat each selected pair as one regression row. For $m$ pairs, write $x_\ell=K_{ij}$ and $z_\ell=r_{ij}$ for pair $\ell$. Ordinary least squares gives

\hat b=\frac{\sum_{\ell=1}^{m}(x_\ell-\bar{x})(z_\ell-\bar{z})}{\sum_{\ell=1}^{m}(x_\ell-\bar{x})^2}, \qquad \hat a=\bar{z}-\hat b\bar{x}.

This formula uses the pair-level means $\bar{x}$ and $\bar{z}$, not the mean phenotype or mean genotype. If the phenotype variance is known or estimated on the same scale, the heritability estimate is $\hat h^2_{\mathrm{SNP}}=\hat b/\hat\sigma_p^2$. In particular, when phenotypes are standardized to variance one, the estimated slope and estimated SNP heritability have the same numerical value. That equivalence depends on the scale; it is not a general identity.

Worked Example: A Positive Slope on a Standardized Scale

In a hypothetical larger cohort, suppose the phenotype has been standardized to phenotypic variance one. For transparent arithmetic, consider four selected, distinct pairs. The table gives each pair's genomic relationship and adjusted-phenotype product. The pair products are hypothetical observations, not population expectations.

Pair	$K_{ij}$	$r_{ij}$
1	0.00	0.001
2	0.10	0.005
3	0.20	0.009
4	0.30	0.013

The pair-level means are $\bar{x}=0.15$ and $\bar{z}=0.007$. The centered relationship values are $-0.15,-0.05,0.05,0.15$, and the centered products are $-0.006,-0.002,0.002,0.006$. Thus

\hat b=\frac{(-0.15)(-0.006)+(-0.05)(-0.002)+(0.05)(0.002)+(0.15)(0.006)} {(-0.15)^2+(-0.05)^2+(0.05)^2+(0.15)^2} =\frac{0.002}{0.05}=0.04.

The intercept is $\hat a=0.007-(0.04)(0.15)=0.001$. Since phenotypic variance is one, $\hat h^2_{\mathrm{SNP}}=0.04/1=0.04$. Under the model and scaling assumptions, the fitted slope estimates genetic variance of 0.04 phenotype-squared units per genomic-relationship unit, corresponding here to an SNP heritability estimate of 0.04. The small positive intercept is the fitted pair-product value at $K_{ij}=0$; it is not an estimate of $\sigma_e^2$.

Why the Slope Estimates Genetic Variance

The connection follows directly from the covariance implied by the model. After fixed effects have been accounted for, write the adjusted phenotype as $y_i=g_i+e_i$. For a distinct pair, expand its product and take its expectation conditional on the relationship:

E(y_i y_j\mid K_{ij}) =E\{(g_i+e_i)(g_j+e_j)\mid K_{ij}\} =\operatorname{Cov}(g_i,g_j\mid K_{ij}) +\operatorname{Cov}(e_i,e_j) =\sigma_g^2K_{ij}.

The cross-covariances between genetic and residual effects are zero by independence. The residual covariance is zero for $i\ne j$ under $e\sim N(0,\sigma_e^2I)$. The genetic covariance is $\sigma_g^2K_{ij}$ by $g\sim N(0,\sigma_g^2K)$. Therefore the model's expected pair product changes by $\sigma_g^2$ for a one-unit increase in $K_{ij}$, which is the slope interpretation.

The diagonal variance has a different role. If $K_{ii}=1$, then the model gives $\operatorname{Var}(y_i)=\sigma_g^2+\sigma_e^2$. Dividing the genetic variance by this total gives SNP heritability. If the relationship matrix uses a different diagonal scale, the genetic variance and heritability interpretation must be adjusted to that scale; one cannot silently apply the unit-diagonal formula.

Worked Example: A Negative Moment Estimate

Consider a second hypothetical calculation with three selected pairs. Suppose the adjusted phenotype has variance one in the source cohort, and the observed pair products and relationships are as follows:

Pair	$K_{ij}$	$r_{ij}$
1	0.10	0.20
2	0.20	0.10
3	0.30	0.05

Here $\bar{x}=0.20$ and $\bar{z}=0.1167$. The centered relationship values are $-0.10,0,0.10$; the centered products are $0.0833,-0.0167,-0.0667$. Substitution gives

\hat b=\frac{(-0.10)(0.0833)+(0)(-0.0167)+(0.10)(-0.0667)} {(-0.10)^2+0^2+(0.10)^2} =\frac{-0.0150}{0.0200}=-0.75.

Then $\hat a=0.1167-(-0.75)(0.20)=0.2667$, and the moment estimate of heritability is $\hat h^2_{\mathrm{SNP}}=-0.75/1=-0.75$. A negative estimate is not a claim that biological genetic variance is negative. It indicates that, in these observed pairs, the fitted product–relationship trend is negative. Unconstrained moment estimates can fall outside the parameter space; this small calculation should not be reported as evidence for a negative biological heritability.

Reading the Intercept and the Scale

The intercept is the fitted pairwise phenotype product when the genomic relationship is zero. Under the idealized covariance model for distinct individuals, its population target is zero. A nonzero fitted value can arise from sampling variation or from a baseline pairwise covariance not represented by the simple model. For example, unaccounted shared environment, residual dependence, or inadequate adjustment for fixed effects can contribute to such a pattern. The intercept alone does not identify which explanation applies.

The intercept also depends on how the predictor is coded. If $K_{ij}$ is replaced by $K_{ij}+c$, where $c$ is a constant, the fitted slope is unchanged and the intercept changes by $-c\hat b$. The fitted values remain the same. This is another reason not to assign the intercept the meaning of residual variance. Under the specified model, residual variance is a diagonal component of the phenotype covariance, whereas this regression uses products for distinct pairs.

For a nonstandardized phenotype, the slope is in phenotype-squared units because $K_{ij}$ is dimensionless. If the phenotype is measured in centimeters, for example, the slope has units of square centimeters per relationship unit. Heritability is dimensionless: divide the genetic-variance estimate by phenotypic variance expressed in the same squared units. When the variance is estimated from data, the numerator and denominator should refer to a compatible phenotype definition and scale.

Worked Example: A Nonzero Intercept Is Not Residual Variance

A third hypothetical dataset has phenotype variance one and four selected pair-level observations:

Pair	$K_{ij}$	$r_{ij}$
1	-0.10	0.02
2	0.00	0.03
3	0.20	0.05
4	0.30	0.06

The means are $\bar{x}=0.10$ and $\bar{z}=0.04$. The centered relationship values are $-0.20,-0.10,0.10,0.20$, and the centered products are $-0.02,-0.01,0.01,0.02$. Therefore

\hat b=\frac{0.004+0.001+0.001+0.004}{0.04+0.01+0.01+0.04} =\frac{0.010}{0.10}=0.10, \qquad \hat a=0.04-(0.10)(0.10)=0.03.

With phenotypic variance one, the estimated SNP heritability is $0.10/1=0.10$. The intercept of 0.03 is the fitted product at zero genomic relationship. It is not the remaining variance $1-0.10=0.90$, nor is it the residual variance. In fact, the model's residual variance estimate would require a variance-component calculation; it cannot be read off as this pair-regression intercept.

Practical Checks and Common Mistakes

A small hand calculation is useful for understanding the estimator, but actual Haseman–Elston analyses generally use many pairs. Pairwise observations are not independent when pairs share individuals, so treating every pair as an independent observation in an ordinary regression standard-error calculation can give misleading uncertainty. The coefficient arithmetic shown here explains the point estimate; it does not provide a valid uncertainty assessment.

Mistake: calling the slope heritability in every setting. The slope estimates genetic variance on the model's phenotype and GRM scales. Divide by compatible phenotypic variance to obtain heritability.
Mistake: calling the intercept residual variance. The intercept describes the fitted product at zero relationship. Residual variance belongs to the diagonal covariance model.
Mistake: treating a negative or greater-than-one estimate as a valid parameter value. An unconstrained moment estimate can leave the allowable range because of finite-sample variation or model mismatch. Report and investigate it; do not interpret it as a biologically possible negative variance.
Mistake: forming products from unadjusted phenotypes without considering fixed effects. Shared mean structure or covariates can create pairwise product patterns that are not the genetic covariance of interest.

The central interpretation is conditional on the additive covariance model, suitable phenotype adjustment, and a clearly stated GRM scale. A positive slope indicates that greater genomic similarity is associated with larger adjusted phenotype products in the fitted pairs; it does not establish that genetic similarity causes the phenotype pattern.

Key takeaway: In pairwise Haseman–Elston regression, the slope estimates additive genetic variance on the chosen scales. SNP heritability is that variance divided by phenotypic variance. The intercept is a baseline pair-product term, not residual variance, and small-sample moment estimates need not lie between zero and one.

Check Your Understanding

Use the model and coefficient definitions above to answer these questions.

If a fitted slope is 0.18 and the phenotype variance is 0.60 on the same scale, what is the corresponding SNP heritability estimate?
Under the stated model, why does the residual variance not contribute to the expected product for a pair of distinct individuals?
A regression has intercept 0.04. What does that value describe, and what does it not estimate?
What happens to the fitted slope and intercept if the same constant is added to every genomic relationship value?
Why can an unconstrained Haseman–Elston moment estimate be negative even though a variance component cannot be negative?