Testing Epistasis Under Linkage Disequilibrium

Why Test Epistasis Under Linkage Disequilibrium?

A two-locus interaction asks whether the association between one genetic marker and a phenotype varies with the genotype at another marker. The challenge is that markers are often correlated because of linkage disequilibrium (LD): some genotype combinations are much more common than others. This dependence affects how precisely an interaction can be estimated and which comparisons the data support.

A test can establish evidence for an interaction in a specified statistical model. It does not, by itself, establish that the two measured markers are causal variants or that a particular molecular mechanism connects them. The conclusion is conditional on how genotypes and the phenotype were coded, what other variables were included, and the assumptions used to calculate uncertainty.

Definition: In a model with two genetic predictors, a statistical interaction is a departure from the model’s specified no-interaction form, represented by a term involving both predictors. Linkage disequilibrium is non-independence of alleles or genotypes at different loci in a population. LD is dependence between predictors; it is not itself evidence that their effects interact.

A Two-Locus Regression Test

Let \(G_1\) and \(G_2\) be the observed genotype dosages at two biallelic markers, coded as 0, 1, or 2 copies of a designated allele. For a quantitative phenotype \(Y\), one common model uses centered dosages \(X_1=G_1-\bar{G}_1\) and \(X_2=G_2-\bar{G}_2\):

Y_i=\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\beta_{12}X_{1i}X_{2i}+\varepsilon_i.

The main effects \(X_1\) and \(X_2\) belong in the model along with their product. The interaction coefficient \(\beta_{12}\) describes how the fitted slope for one marker changes with the other marker. Specifically, for fixed \(X_2\), the fitted slope with respect to \(X_1\) is \(\beta_1+\beta_{12}X_2\). If the dosages have been standardized instead of merely centered, the coefficient’s units and numerical value change accordingly.

For the usual exact small-sample \(t\) test, assume independent individuals, a correctly specified conditional mean, full rank of the model matrix, and errors that are independent and normally distributed with mean zero and common variance given the modeled predictors. The full-rank condition means the product column cannot be written exactly as a combination of the intercept and the two main-effect columns. Under these assumptions, test \(H_0:\beta_{12}=0\) against an alternative such as \(H_a:\beta_{12}\ne 0\). If the errors are not normal, large-sample tests may still be used with suitable standard errors, but their justification is different.

LD does not automatically bias the interaction estimate. If the model is correctly specified, the predictors have full rank, and the conditional mean assumption holds, correlated predictors can still be included together. However, their dependence can make some coefficient estimates less precise or unstable. With strong LD, some genotype combinations may be rare or absent, which can leave little information about the interaction or make it impossible to estimate.

Worked Example: Testing an Interaction Coefficient

Suppose a hypothetical study of 240 unrelated individuals measures a quantitative phenotype in score units and two centered marker dosages standardized to have unit variance. The markers have genotype correlation \(r=0.65\). A regression including the intercept, both main effects, and their product estimates \(\hat{\beta}_{12}=0.80\) score units per product unit, with \(SE=0.30\). Assume the stated linear model and exact \(t\)-test conditions are reasonable, and that the four regression columns have full rank.

There are four fitted parameters, so the residual degrees of freedom are \(df=240-4=236\). For \(H_0:\beta_{12}=0\), the test statistic is:

t=\frac{\hat{\beta}_{12}-0}{SE} =\frac{0.80}{0.30}=2.6667.

The two-sided p-value is \(2P(T_{236}\ge 2.6667)\), approximately \(0.008191\). At significance level \(\alpha=0.05\), this is evidence against a zero interaction coefficient in this model. The p-value is calculated assuming the null hypothesis and the model assumptions; it is not the probability that the null hypothesis is true.

Using the approximate 95% critical value \(t_{0.975,236}=1.970\), a confidence interval is:

0.80\pm 1.970(0.30) =0.80\pm0.591 \approx(0.209,\ 1.391).

Under repeated sampling and the model assumptions, this procedure has 95% coverage for the interaction coefficient. In this example, the interval is in score units per product unit and excludes zero. The estimate describes a change in the fitted marker slope across values of the other centered, standardized marker. It does not show that either marker causes the phenotype.

What LD Changes: Information and Identification

The interaction term is estimated from how the outcome varies across joint marker values, beyond the patterns accounted for by the intercept and main effects. If LD concentrates observations in only a few combinations, comparisons involving rare combinations carry high uncertainty. If a combination is absent, the model may not have enough independent information to identify all the intended effects.

A useful diagnostic is to inspect genotype-combination counts, not only the correlation between marker dosages. In a regression analysis, also check that the design matrix has full rank and that the interaction estimate is not driven by a handful of observations. A large standard error, sensitivity to reasonable coding choices, or a near-singular model matrix signals weak information. The sign or magnitude of LD alone does not determine whether an interaction coefficient is estimable or precise.

Worked Example: Sparse Genotype Combinations Under Strong LD

Consider a hypothetical study of 100 individuals with two binary genotype indicators \(X_1,X_2\), each coded 0 or 1. The joint counts are 45 for \((0,0)\), 5 for \((1,0)\), 5 for \((0,1)\), and 45 for \((1,1)\). Each indicator has mean 0.5 and variance 0.25. Their covariance is \(45/100-(0.5)(0.5)=0.20\), so their correlation is \(0.20/0.25=0.80\). The genotype indicators are strongly correlated, but all four combinations are present.

With a binary-by-binary model containing both main effects and the product, the product coefficient is the difference-in-differences of the four cell means:

\hat{\beta}_{12} =\bar{Y}_{11}-\bar{Y}_{10}-\bar{Y}_{01}+\bar{Y}_{00}.

Suppose the hypothetical cell means are \(\bar{Y}_{00}=0\), \(\bar{Y}_{10}=0.2\), \(\bar{Y}_{01}=-0.1\), and \(\bar{Y}_{11}=0.9\) score units. Then:

\hat{\beta}_{12}=0.9-0.2-(-0.1)+0=0.80\text{ score units}.

Assume independent normal residuals with known standard deviation 1 score unit in every cell. Since the cell means are based on independent individuals, the variance of this contrast is the sum of the four cell-mean variances:

\operatorname{Var}(\hat{\beta}_{12}) =1\left(\frac{1}{45}+\frac{1}{5}+\frac{1}{5}+\frac{1}{45}\right) =0.4444, \qquad SE=\sqrt{0.4444}=0.6667.

The known-variance normal test statistic is \(0.80/0.6667=1.20\), giving a two-sided p-value of approximately \(0.230\). Thus, this hypothetical estimate is not statistically significant at \(\alpha=0.05\), despite being 0.80 score units. The two discordant genotype combinations each contain only five individuals, so the contrast is imprecise.

For comparison, if the four combinations each had 25 individuals, with the same assumed residual standard deviation, the standard error would be \(\sqrt{4/25}=0.40\) score units. The estimate’s precision depends on the joint genotype distribution, not simply on the total sample size. If any cell were empty, this saturated binary-marker model would not identify all four cell means from the data.

Confounding, Interpretation, and Model Scale

LD is not the same as confounding. LD describes dependence between genetic predictors. Confounding arises when a variable related to both the predictors and the phenotype is not adequately accounted for, leading the fitted relationship to reflect more than the association of interest. For example, ancestry differences can be associated with allele frequencies and with a phenotype for environmental or other genetic reasons. A marker interaction can then partly reflect population structure rather than a within-population interaction.

The interaction coefficient is conditional on the model’s included terms and covariates. A diagnostic review should therefore ask whether the model includes the relevant main effects, whether important design or ancestry variables have been handled appropriately, and whether the result is robust to defensible specifications. LD may also make it difficult to distinguish which of two correlated markers carries a signal. A significant interaction between measured markers does not prove that those markers themselves interact biologically; they may tag other variants.

Interaction also depends on the scale of the outcome model. For a quantitative outcome, the product coefficient in linear regression measures departure from additivity of the fitted mean on the outcome scale. For a binary outcome, a logistic model describes log odds. Its interaction coefficient measures departure from additivity on the log-odds scale, not necessarily on the risk scale. A zero interaction on one scale need not be zero on another.

Worked Example: Risk-Scale and Log-Odds-Scale Interaction

Suppose hypothetical outcome risks for two binary genotype indicators are 0.10 for \((0,0)\), 0.20 for \((1,0)\), 0.20 for \((0,1)\), and 0.30 for \((1,1)\). On the risk scale, the interaction contrast is:

$$ 0.30-0.20-0.20+0.10=0. $$

Thus, these risks show no interaction as a difference-in-differences on the additive risk scale. In a logistic model, the corresponding interaction contrast is the same difference-in-differences applied to the cell logits:

\operatorname{logit}(0.30)-\operatorname{logit}(0.20) -\operatorname{logit}(0.20)+\operatorname{logit}(0.10).

Using \(\operatorname{logit}(p)=\log\{p/(1-p)\}\), the four logits are approximately \(-0.8473\), \(-1.3863\), \(-1.3863\), and \(-2.1972\). Therefore the log-odds interaction is:

-0.8473-(-1.3863)-(-1.3863)+(-2.1972) \approx -0.272.

The interaction is zero on the additive risk scale but not on the log-odds scale. These are not contradictory results: they answer different questions. A report of an interaction test should name the model and scale rather than presenting “interaction” as an outcome-independent property.

A Short Justification of the Test

Under the linear model assumptions, the estimated interaction coefficient is centered at its population model coefficient. Under \(H_0:\beta_{12}=0\), its estimation error divided by its standard error follows a \(t\) distribution with the residual degrees of freedom. The test compares the observed standardized estimate with that null distribution; a two-sided p-value counts outcomes at least as far from zero in either direction.

The product term must contribute information not already contained in the intercept and main effects. If it does not, the design matrix is rank-deficient and the coefficient cannot be uniquely estimated. If it contributes only a small amount of independent information, its standard error can be large. In this way, the usual test calculation connects model identification and precision directly to the observed genotype combinations.

Common Mistakes and Practical Checks

Leaving out a main effect. Testing a product term without both marker main effects changes the model and can make the product pick up ordinary main-effect patterns. Include both main effects when testing their interaction.
Calling LD an interaction. Correlated genotypes do not imply that one marker’s effect changes with the other. The interaction is a coefficient or contrast in a specified outcome model.
Treating a non-significant result as proof of no interaction. Sparse genotype combinations or a small sample can yield low power. Report the estimate and uncertainty, and inspect the joint genotype counts.
Calling the marker interaction causal epistasis. The fitted coefficient is a statistical association conditional on the model. LD with unmeasured variants, population structure, or other model limitations may affect its interpretation.
Ignoring the scale. An interaction on the log-odds scale does not have the same meaning as an interaction on the risk scale. State the phenotype model and the scale of the coefficient.

A careful analysis therefore reports the genotype coding, model terms, estimate and standard error, test and degrees of freedom (when applicable), p-value, and a contextual interpretation. It also describes the degree of LD and checks whether the observed joint genotypes provide enough support for the comparison the interaction test requires.

Check Your Understanding

Use the model assumptions and examples above to answer these questions.

Why should a two-locus product term ordinarily be fitted alongside both marker main effects?
In the regression example, what does the estimate \(\hat{\beta}_{12}=0.80\) describe, and why does it not establish a causal biological interaction?
How can strong LD reduce precision even when the total sample size is large?
What is the distinction between LD and confounding? Give one possible source of confounding in a genetic study.
How can the same binary-outcome data show no interaction on the risk scale but a nonzero interaction in a logistic model?