Why Linkage Disequilibrium Matters for Interaction Tests
A pairwise interaction test asks whether the joint contribution of two marker predictors improves a specified model beyond their main effects. Its answer depends not only on the phenotype and the model, but also on the observed distribution of genotypes. When markers are in linkage disequilibrium (LD), their predictors are statistically dependent. Their product can then resemble one or more main-effect predictors, or another interaction term. The resulting design matrix may contain little information for distinguishing the corresponding coefficients.
This is a problem of estimability and precision, not by itself evidence for or against biological interaction. A fitted model may have a numerically unstable interaction estimate even when it is correctly specified, because the observed marker combinations do not provide independent variation in the relevant columns. Conversely, LD does not automatically make every interaction test unstable: the relevant question is how much independent variation remains in the particular interaction column after accounting for the other predictors.
The haplotype definition describes dependence between alleles; a regression design matrix contains individual-level genotype columns, often coded as dosages 0, 1, or 2. These are related but not interchangeable descriptions. The examples below state their coding explicitly. In a hierarchical regression model, the interaction term is included along with both corresponding main effects. For a quantitative outcome, a simple model is \(Y_i=\beta_0+\beta_Ax_{Ai}+\beta_Bx_{Bi}+\beta_{AB}x_{Ai}x_{Bi}+\epsilon_i\). The interaction coefficient describes the extra product-term contribution on the model's outcome scale, conditional on the included terms; LD affects how well it can be estimated, not the definition of that coefficient.
Exact and Near Aliasing
Two design-matrix columns are exactly aliased if one is a linear combination of the other columns. Then the model cannot uniquely estimate all the associated coefficients: different coefficient values produce the same fitted values. Near aliasing means the columns are not exactly dependent, but one column is close to a linear combination of the others. The coefficients are then identifiable in principle, but small changes in the data can produce large changes in their estimates.
LD can create either situation. If two marker predictors are strongly correlated, a product involving one marker can track a main effect or a different product. The pattern depends on allele frequencies, genotype coding, and the joint distribution of the markers. For example, multiplying two similar nonnegative dosage columns can yield a product that varies much like either dosage. With several markers, interaction terms can also be correlated with each other because they share predictors and because the markers themselves are dependent.
A marginal correlation is a useful warning, but it is not a complete diagnosis. The interaction coefficient is estimated jointly with all other included predictors. The more direct question is whether the interaction column still varies after its linear association with the rest of the design matrix has been removed.
Worked Example: LD and an Interaction Column That Resembles Main Effects
Consider a hypothetical sample of 100 individuals with two binary marker indicators \(A\) and \(B\), coded 0 or 1. The joint counts are 49 with \((A,B)=(0,0)\), 1 with \((0,1)\), 1 with \((1,0)\), and 49 with \((1,1)\). Let the product predictor be \(W=AB\).
Both marker indicators have sample proportions \(0.50\), and the proportion with \(A=B=1\) is \(0.49\). Thus the empirical LD measure is \(D=0.49-(0.50)(0.50)=0.24\), and \(r^2=0.24^2/[(0.50)(0.50)(0.50)(0.50)]=0.9216\). The marker indicators are highly correlated: their covariance is \(0.49-(0.50)(0.50)=0.24\), so their correlation is \(0.24/\sqrt{(0.25)(0.25)}=0.96\).
The interaction column is also highly correlated with each main effect. For example, \(E(AW)=0.49\), \(E(A)=0.50\), and \(E(W)=0.49\), so \(\operatorname{Cov}(A,W)=0.49-(0.50)(0.49)=0.245\). Since \(\operatorname{Var}(A)=0.25\) and \(\operatorname{Var}(W)=0.49(0.51)=0.2499\), the correlation is \(0.245/\sqrt{(0.25)(0.2499)}\approx0.9802\).
However, the joint model includes both \(A\) and \(B\), so residualizing \(W\) against just one main effect is not enough. Regress \(W\) on an intercept, \(A\), and \(B\). By symmetry, the fitted coefficients are intercept \(-0.01\) and slopes \(0.50\) for both \(A\) and \(B\). The residuals for the four patterns \(00,01,10,11\) are \(0.01,-0.49,-0.49,0.01\), respectively. Their mean squared value is \((98(0.01^2)+2(0.49^2))/100=0.0049\), while the variance of \(W\) is \(0.2499\). Only \(0.0049/0.2499\approx0.01961\), or about 1.96%, of the interaction column's variation remains after adjustment for the intercept and main effects. This is near aliasing: the product is not exactly redundant here, but its independent information is limited.
Residualizing the Interaction: A Design-Matrix Justification
The role of the remaining interaction variation follows directly from least squares. Let \(y\) be the outcome vector, \(w\) the interaction column being tested, and \(Z\) the matrix containing an intercept and all other model columns, including the relevant main effects and any other included interactions. Assume \(Z\) has full column rank. Define \(P_Z=Z(Z^\mathsf{T}Z)^{-1}Z^\mathsf{T}\) and \(M_Z=I-P_Z\). The matrix \(M_Z\) removes the part of a vector explained by \(Z\). Consider the model \(y=Z\gamma+w\beta+\epsilon\).
For any fixed \(\beta\), the least-squares residual after fitting \(\gamma\) is \(M_Z(y-w\beta)\), because \(M_Z\) projects onto the space orthogonal to the columns of \(Z\).
Minimize \([M_Z(y-w\beta)]^\mathsf{T}[M_Z(y-w\beta)]\). Differentiating with respect to \(\beta\) gives \(-2w^\mathsf{T}M_Z(y-w\beta)=0\), since \(M_Z\) is symmetric and idempotent.
When \(w^\mathsf{T}M_Zw>0\), solving the equation gives \(\hat{\beta}=w^\mathsf{T}M_Zy/(w^\mathsf{T}M_Zw)\). If \(w^\mathsf{T}M_Zw=0\), then \(w\) lies in the span of \(Z\) and the separate interaction coefficient is not identifiable.
For the linear model with independent errors of common variance \(\sigma^2\), the variance of this coefficient estimate, conditional on the design, is \(\operatorname{Var}(\hat{\beta})=\sigma^2/(w^\mathsf{T}M_Zw)\). Indeed, \(\hat{\beta}-\beta=w^\mathsf{T}M_Z\epsilon/(w^\mathsf{T}M_Zw)\). The variance of its numerator is \(\sigma^2 w^\mathsf{T}M_ZM_Zw=\sigma^2 w^\mathsf{T}M_Zw\), because \(M_Z^2=M_Z\). Dividing by the squared denominator gives the stated result. This formula assumes the stated linear model and error variance conditions; normal errors are additionally needed for exact normal-theory inference in a finite sample.
Thus the residualized sum of squares \(w^\mathsf{T}M_Zw\) is a direct measure of information for the target coefficient in this model. A small positive value means a large standard error relative to the outcome noise. A value of zero means exact aliasing. For generalized linear models, the corresponding local information is weighted by the model's information matrix, but the same geometric concern remains: the target column must contribute information not already represented by the other columns.
Worked Example: Translating Residual Variation Into Precision
Use the hypothetical binary-marker table above in a quantitative-trait model with columns \(1,A,B,W\). Suppose the linear-model errors are independent with common variance \(\sigma^2=1\) trait-unit squared. From the residualization calculation, \(w^\mathsf{T}M_Zw=100(0.0049)=0.49\).
The conditional standard error of the interaction estimate is therefore \(SE(\hat{\beta}_{AB})=\sqrt{1/0.49}=1.4286\) trait units per product-predictor unit. If the interaction column had retained all its original centered sum of squares, that sum would have been \(100(0.2499)=24.99\), corresponding to a standard error of \(\sqrt{1/24.99}\approx0.2001\) under the same error-variance assumption. The ratio is \(1.4286/0.2001\approx7.14=\sqrt{51}\). This comparison isolates the loss of design information; in a real fitted analysis, the residual variance is estimated and the reported standard error also reflects that estimate.
This calculation does not imply that the interaction coefficient is large, small, or statistically significant. It shows that the design provides relatively little precision for separating the product term from the main effects. A confidence interval or test must still use the fitted coefficient and its standard error, with the model assumptions and outcome scale stated.
Dosage Coding and Centering
With additive genotype coding, the marker predictors are often dosages \(G_A,G_B\in\{0,1,2\}\), and the product is \(G_AG_B\). LD can make the dosage columns strongly correlated, and the product can correlate with them as well. The strength and direction of these relationships depend on the dosage distribution. Centering changes the numerical form of the product: \((G_A-\bar{G}_A)(G_B-\bar{G}_B)\) is not the same column as \(G_AG_B\). When an intercept and both main effects are included, these two product columns differ only by a linear combination of the intercept and main effects, so they span the same model space. Centering can improve interpretation or numerical behavior, but does not create independent genotype combinations that were absent from the sample.
Worked Example: A Dosage Product Correlated With a Main Effect
Suppose a hypothetical sample has 10 individuals, with marker dosages equal at the two loci for every individual: five have \((G_A,G_B)=(0,0)\), four have \((1,1)\), and one has \((2,2)\). Write the common dosage as \(G\), so the product column is \(G_AG_B=G^2\).
The sample moments are \(E(G)=0.6\), \(E(G^2)=0.8\), \(E(G^3)=1.2\), and \(E(G^4)=2.0\). Therefore \(\operatorname{Var}(G)=0.8-0.6^2=0.44\), \(\operatorname{Var}(G^2)=2.0-0.8^2=1.36\), and \(\operatorname{Cov}(G,G^2)=1.2-(0.6)(0.8)=0.72\). The correlation is \(0.72/\sqrt{(0.44)(1.36)}\approx0.9305\).
The interaction product is strongly correlated with each dosage main effect in this constructed sample. Yet it is not exactly a linear combination of an intercept and \(G\), because the three observed dosage levels give \(G^2\) values \(0,1,4\), which do not lie on a single straight line as a function of \(G\). The interaction is therefore distinguishable in principle, although its coefficient may be imprecise. This example also shows why the coding and observed genotype distribution must be examined rather than assuming that LD implies exact aliasing.
When Separate Coefficients Cannot Be Estimated
Exact aliasing is easiest to see algebraically. If a marker coding makes two columns identical, or makes an interaction column identical to a main-effect column, the model can estimate only a combination of their coefficients. No choice of test statistic can recover separate coefficients from a design matrix that does not distinguish them. This is different from a large standard error: exact aliasing is non-identifiability, whereas near aliasing retains some information.
Worked Example: Perfectly Dependent Binary Marker Predictors
In a hypothetical sample, suppose two binary marker indicators always agree: every individual has either \((A,B)=(0,0)\) or \((1,1)\). In the hierarchical model, the interaction column is \(W=AB\). For every observed individual, \(B=A\) and \(AB=A\). The design-matrix columns for \(A\), \(B\), and \(AB\) are therefore identical.
The linear predictor becomes \(\beta_0+\beta_AA+\beta_BB+\beta_{AB}AB =\beta_0+(\beta_A+\beta_B+\beta_{AB})A\). Only the sum \(\beta_A+\beta_B+\beta_{AB}\) can be informed by these observations. For example, adding 1 to \(\beta_{AB}\) and subtracting 1 from \(\beta_A\) leaves every fitted value unchanged. The individual coefficients are not separately estimable. This is exact aliasing, not merely a large standard error or a nonsignificant interaction test.
Practical Diagnosis and Interpretation
A useful diagnostic sequence starts with the genotype predictors and then examines the actual model matrix. Review marker-genotype correlations and the distribution of joint genotypes; calculate correlations among product terms and between products and main effects; then assess each target interaction after residualizing it against the full set of other model columns. Matrix rank, small eigenvalues of the information matrix, condition numbers, or a large variance inflation factor can flag dependence. A high condition number is scale-dependent, so predictors should be placed on sensible, comparable scales before using it as a diagnostic.
The important object is the model as fitted. If an analysis includes covariates or multiple interaction terms, a target interaction must be assessed against those columns too. Pairwise correlations alone can miss a dependence involving several columns. Conversely, a high pairwise correlation is a warning rather than a proof that the target coefficient is not estimable; the residualized column determines whether independent information remains.
When estimates are unstable, report the limitation rather than interpreting a large coefficient as strong evidence. Consider whether the scientific question requires separate effects for highly dependent markers, and describe the coding and relevant genotype dependence. Reparameterizing or centering can make a model easier to compute or interpret, but it cannot resolve a true rank deficiency. LD-induced aliasing is therefore a reason to qualify what the data can distinguish, not a reason to infer causation or biological interaction from the LD pattern itself.
Common Mistake
- Wrong reasoning: “The two markers are in high LD, so their interaction must be unestimable.” Correction: High LD can produce near aliasing, but whether a particular interaction is estimable depends on its independent variation after adjustment for the full model matrix. Check rank and residualized information.
- Wrong reasoning: “The interaction column has a high correlation with a main effect, so the interaction coefficient equals that main effect.” Correction: Correlation does not make coefficients equivalent. Exact aliasing requires a linear dependence among columns; near aliasing instead inflates uncertainty.
- Wrong reasoning: “Centering the dosages removes LD and fixes the interaction test.” Correction: Centering can alter the product column by combinations of main-effect columns, but it does not change the underlying joint genotype information or cure exact dependence.
- Wrong reasoning: “An unstable interaction estimate shows that there is no biological interaction.” Correction: Instability describes limited statistical separation of model terms in the observed design. It is not evidence that an interaction is absent.
Check Your Understanding
Use the definitions and calculations above to answer these questions.
- What distinguishes exact aliasing from near aliasing in a pairwise interaction model?
- In the 100-individual binary-marker example, why is the residual variation of \(W=AB\) after adjustment for both main effects more informative than its marginal correlation with \(A\)?
- Under the stated linear-model assumptions, what happens to the interaction coefficient's variance as \(w^\mathsf{T}M_Zw\) becomes small?
- Why can centering both dosages change the product column without creating new independent genotype information?
- If two binary marker predictors always agree and their product equals the same indicator, which combination of coefficients is estimable from the observed design?