Variance-Component Tests For Distributed Epistasis

Why Test a Shared Interaction Variance?

A pairwise interaction analysis can include many genotype-product terms. If there are too many terms to estimate precisely, testing each coefficient separately may have little power: each interaction could be small, yet the collection could still contribute meaningful variation. A variance-component model instead asks whether the interaction effects, considered as a group, vary around zero enough to explain trait variation beyond the fixed covariates and main effects.

This is a model-based test of distributed epistasis. It does not identify which pair has an effect, and a positive result does not establish a biological mechanism. The central question is whether the data support a nonzero variance for a specified collection of interaction coefficients.

Random-Effect Model for Pairwise Products

Let $y$ be the $n$-vector of a quantitative trait, and let $X$ contain an intercept, covariates, and the main-effect genotype terms required by the model. Let $Z$ be an $n \times q$ matrix whose columns are the specified pairwise genotype products. For example, a column could contain the product of the two coded genotypes at a particular locus pair. The model treats the columns of $Z$ as specified predictors while treating their coefficients as random.

Definition: The Gaussian variance-component model for distributed interaction effects is $$ y = X\beta + Zb + \varepsilon,\qquad b \sim N(0,\tau I_q),\qquad \varepsilon \sim N(0,\sigma^2 I_n), $$ where $b$ and $\varepsilon$ are independent. The fixed-effect coefficients are $\beta$, $\tau \ge 0$ is the common variance of the interaction coefficients, and $\sigma^2 > 0$ is the residual variance. Conditional on the genotype design matrices, integrating over $b$ gives $$ y \sim N(X\beta,V_\tau),\qquad V_\tau=\sigma^2 I_n+\tau ZZ'. $$

The random-effects assumption says that, conditional on the included fixed effects and genotype products, the interaction coefficients have mean zero, share a variance, and are independent of one another. This is a working distributional model, not a claim that the particular interaction coefficients are literally random draws in a biological process. With dimensionless genotype coding, $b_j$ is measured in trait units and $\tau$ in squared trait units. Changing the scale of the genotype products changes the interpretation and numerical scale of $\tau$.

The null and alternative are

H_0:\tau=0 \qquad\text{versus}\qquad H_a:\tau>0.

At $H_0$, the interaction contribution $Zb$ has zero variance and the model reduces to the fixed-effects model. Because a variance cannot be negative, $H_0$ is on the boundary of the allowed parameter space. That boundary is why standard tests for an unconstrained coefficient cannot simply be applied.

Worked Example: Setting Up a Distributed Interaction Model

Suppose a hypothetical study has a quantitative trait measured in 500 participants and specifies 40 pairwise genotype products. The fixed-effect matrix $X$ includes an intercept, age, ancestry covariates, and the two locus main effects for the specified pairs. The matrix $Z$ contains the 40 product columns.

The model is $y=X\beta+Zb+\varepsilon$, with $b_j$ independently distributed as $N(0,\tau)$. If, for example, the trait is measured in milligrams per deciliter and the genotype products are dimensionless, each interaction coefficient is measured in milligrams per deciliter and $\tau$ in squared milligrams per deciliter.

The test is $H_0:\tau=0$ versus $H_a:\tau>0$. It asks whether the 40 effects collectively have variation beyond residual variation, conditional on the terms in $X$. It does not test 40 separate null hypotheses, nor does it say which of the 40 coefficients is responsible if the null is rejected.

Score Test at Zero Interaction Variance

A score test evaluates the slope of the likelihood as $\tau$ moves away from zero, rather than fitting a separate unrestricted coefficient for every interaction. For a transparent derivation, first suppose $\sigma^2$ is known. Define $R=ZZ'$, so $V_\tau=\sigma^2 I_n+\tau R$. Under the null, let

M=I_n-X(X'X)^{-1}X'.

Here $M$ projects away the fitted fixed effects; assume $X$ has full column rank. The residual vector at the null is $My$. The restricted log likelihood, up to constants that do not depend on $\tau$, can be written using $P=V_\tau^{-1}-V_\tau^{-1}X(X'V_\tau^{-1}X)^{-1}X'V_\tau^{-1}$. Its derivative at $\tau=0$ is the score

U_\tau(0)=\frac{1}{2}\left[y'P R P y-\operatorname{tr}(P R)\right], \qquad P=\frac{M}{\sigma^2}\text{ at }\tau=0.

The two terms have a direct interpretation: the first measures interaction-related residual variation, and the trace term is its expectation under the null. A large positive score indicates more residual variation in directions associated with $Z$ than the null model predicts. The score itself is centered, so it can be negative; a negative value is not evidence for a negative variance.

Justification: For a Gaussian model with covariance $V_\tau$, the derivative of the log likelihood with respect to $\tau$ is one half the quadratic form in $V_\tau^{-1}(\partial V_\tau/\partial\tau)V_\tau^{-1}$, minus one half the trace of $V_\tau^{-1}(\partial V_\tau/\partial\tau)$. Here $\partial V_\tau/\partial\tau=R$. When fixed effects are estimated, replacing the inverse covariance by the residualized matrix $P$ gives the restricted-likelihood score above. At $\tau=0$, $P=M/\sigma^2$, yielding the stated formula.

The score's null distribution is generally not an ordinary chi-square distribution. To see why, let $\lambda_1,\ldots,\lambda_{n-p}$ be the eigenvalues, including any zeros, of $MRM/\sigma^2$ on the residual subspace, where $p$ is the number of columns of $X$. Under $H_0$, the residuals are Gaussian with variance $\sigma^2$. Orthogonally rotating them into the eigenvector directions gives

Q=y'P R P y \ \overset{H_0}{\sim}\ \sum_{i=1}^{n-p}\lambda_i\chi^2_{1,i}, \qquad U_\tau(0)=\frac{1}{2}\left(Q-\sum_{i=1}^{n-p}\lambda_i\right),

where the $\chi^2_{1,i}$ variables are independent chi-square random variables with one degree of freedom. Unequal eigenvalues produce a weighted sum, so calibration must account for the design and the interaction columns. When the positive eigenvalues are all equal, the sum reduces to a scaled chi-square distribution. Otherwise, software can evaluate or approximate the weighted-sum tail probability. If nuisance parameters such as $\sigma^2$ are estimated, the null calibration must account for that estimation, for example through a restricted-likelihood procedure or a parametric bootstrap under the fitted null.

Worked Example: A Score Test With Two Interaction Directions

Consider hypothetical, residualized data with known $\sigma^2=1$. Suppose the two nonzero eigenvalues of $MRM/\sigma^2$ are both 2, and the squared residual coordinates in those two directions are 3 and 2. All other residual coordinates have zero interaction weight.

The observed quadratic form is $$ Q=2(3)+2(2)=10. $$ The trace term is $2+2=4$, so the score is $$ U_\tau(0)=\frac{1}{2}(10-4)=3. $$

Because the two weights are equal, the null distribution of $Q$ is $2\chi^2_2$. Thus the upper-tail p-value is $$ P(2\chi^2_2\ge 10)=P(\chi^2_2\ge 5)=e^{-5/2}=0.0821. $$ At a 0.05 significance level, this result does not reject $H_0$. It is not evidence that every interaction effect is zero; it means this score test does not provide sufficient evidence for a positive shared variance at that threshold.

Likelihood-Ratio Testing and the Boundary

A likelihood-ratio test compares the largest likelihood under $\tau\ge 0$ with the likelihood under $\tau=0$. In practice, the fixed effects and other nuisance parameters are estimated as part of the fit. The statistic is

\mathrm{LRT}=2\{\ell(\widehat{\tau},\widehat{\beta},\widehat{\sigma}^2) -\ell(0,\widehat{\beta}_0,\widehat{\sigma}_0^2)\}, \qquad \widehat{\tau}\ge 0.

For a single variance component under standard large-sample conditions, the null reference distribution is a 50:50 mixture: half its probability is at zero, and half follows a $\chi^2_1$ distribution. This differs from the ordinary $\chi^2_1$ reference distribution because estimates of $\tau$ cannot cross below zero. The mixture result is an asymptotic result, not a guarantee for a small or poorly identified design. If there are multiple variance components or additional irregularities, the reference distribution may differ; calibration should match the fitted model.

Worked Example: A One-Contrast Likelihood-Ratio Calculation

For a small hypothetical illustration of the likelihood calculation, suppose one residual contrast $Y^*$ has variance $1+\tau$, with known residual variance 1, and the observed contrast is $Y^*=2$. Under the alternative, the log likelihood, omitting constants, is $$ \ell(\tau)=-\frac{1}{2}\left\{\log(1+\tau)+\frac{4}{1+\tau}\right\}. $$

The likelihood is maximized at $\widehat{\tau}=3$, since the fitted variance $1+\widehat{\tau}$ equals the observed squared contrast, 4. The log-likelihood difference from $\tau=0$ is $$ \ell(3)-\ell(0) =-\frac{1}{2}(\log 4+1)-\left(-2\right) =0.8069. $$ Therefore, $\mathrm{LRT}=2(0.8069)=1.614$.

In this deliberately simple model, the null is $Y^*\sim N(0,1)$. The likelihood-ratio statistic is at least as large as the observed value exactly when $|Y^*|\ge 2$. Its exact null p-value is therefore $P(\chi^2_1\ge4)=0.0455$. This calculation uses the exact one-contrast null distribution; it should not be confused with the large-sample mixture approximation. The example illustrates likelihood mechanics, not the usual data-rich setting for a distributed interaction analysis.

Assumptions, Interpretation, and Common Pitfall

The model and its test are meaningful only relative to the specified analysis. The Gaussian derivations above assume independent residual errors with common variance, a correctly specified fixed-effect mean, and genotype-product columns treated as fixed conditional on the observed genotypes. The random interaction coefficients have mean zero and a shared variance; the model also assumes their covariance matrix is $\tau I_q$. If interaction effects are expected to be correlated or to have different variances, that is a different variance-component model. For binary outcomes, a corresponding generalized mixed model requires an appropriate likelihood and calibration rather than using the Gaussian formulas unchanged.

The fixed effects should include the intended covariates and the relevant genotype main effects so that the variance test addresses interaction variation conditional on those terms. The coding and scaling of genotype products matter: they determine what a unit of interaction coefficient means and can affect the fitted variance. A significant test supports variation attributable to the specified interaction-effect model; it does not show that each pair contributes, locate a causal pair, or imply a biological interaction.

Common mistake: “There are 40 interaction columns, so the likelihood-ratio statistic should be compared with $\chi^2_{40}$.” This treats the model as if it added 40 unconstrained fixed coefficients.
Correction: The random-effects model adds one variance parameter, $\tau$, constrained to be nonnegative. Under standard large-sample conditions its likelihood-ratio reference is the boundary mixture for one variance component, not $\chi^2_{40}$. The score test instead has a design-dependent weighted chi-square null distribution.

A variance-component test is useful when the scientific target is a collection of interactions rather than a single prespecified pair. It can combine weak signals without requiring each effect to pass an individual test, but the result is necessarily aggregate. Report the interaction set, genotype coding, fixed-effect adjustment, variance-component estimate, test statistic, and calibration method so that the conclusion remains tied to the model actually fitted.

Check Your Understanding

Use the model assumptions and test definitions to answer the following.

In the random-effect model, what does $\tau=0$ imply about the interaction contribution to the covariance of $y$?
Why is testing $H_0:\tau=0$ a boundary problem?
When the eigenvalues in the score-test null distribution are unequal, why is an ordinary chi-square reference generally inappropriate?
A significant variance-component test is obtained for a set of 25 genotype products. What can be concluded, and what cannot be concluded, about the individual pairs?
Why does treating 25 interaction coefficients as random effects not automatically make the likelihood-ratio reference distribution $\chi^2_{25}$?