Likelihood-Ratio And Score Tests For Epistasis

Why Test the Interaction Jointly?

A two-locus model typically includes effects for each locus as well as terms representing how their joint genotype pattern departs from the model with locus effects alone. The scientific question addressed here is whether those interaction terms, considered as a group, improve the model after the main effects and other prespecified covariates have been retained. This is a joint test: it is not the same as testing one selected interaction coefficient or comparing a full model with an intercept-only model.

The likelihood-ratio and score tests answer this question in different ways. The likelihood-ratio test compares the maximized likelihoods under the full and reduced models. The score test measures how strongly the likelihood would initially increase in an interaction direction, using estimates obtained under the reduced model. Both have an asymptotic chi-square reference distribution under regularity conditions.

Definition: Write the model parameters as $(\eta,\gamma)$, where $\eta$ collects nuisance parameters (including locus main effects and any adjustment covariates) and $\gamma$ is the vector of interaction coefficients. The joint null hypothesis of no interaction is $H_0:\gamma=\mathbf{0}$; the alternative is $H_a:\gamma\ne\mathbf{0}$. The reduced model imposes $H_0$, while the full model estimates $\eta$ and $\gamma$.

The interaction degrees of freedom are the number $q$ of linearly independent interaction parameters being tested, provided the design is identifiable. With a categorical genotype factor having three levels at each locus and a full set of locus main effects, the saturated interaction has $(3-1)(3-1)=4$ degrees of freedom. A model that tests only the product of two dosage variables instead has one interaction coefficient and therefore one degree of freedom. The coding and scientific hypothesis determine the test; the degrees of freedom must match the actual joint restriction.

Likelihood-Ratio Test

Let $\ell_1$ be the maximized log-likelihood in the full model and $\ell_0$ the maximized log-likelihood in the reduced model. The models must use the same observations, outcome distribution, main effects, covariates, and other assumptions; the full model adds the $q$ interaction parameters. Since the full model includes the reduced model as a special case, $\ell_1\geq\ell_0$.

Theorem: Under $H_0$, if the models are correctly specified and regular, the null parameter is identifiable and lies in the interior of the parameter space, and the sample is sufficiently large, the likelihood-ratio statistic $$ \Lambda=2(\ell_1-\ell_0) $$ has an approximate $\chi^2_q$ distribution. Its p-value is $P(\chi^2_q\geq\Lambda_{\mathrm{obs}})$.

The factor of two follows from the standard likelihood-ratio scaling: for a regular nested comparison, twice the log-likelihood improvement converges in distribution to a chi-square random variable with degrees of freedom equal to the number of independent restrictions. Here those restrictions are the $q$ interaction coefficients set to zero. The p-value measures how often a statistic at least this large would arise under the no-interaction model; it is not the probability that $H_0$ is true.

Worked Example: Quantitative Trait Likelihood-Ratio Test

Suppose a hypothetical study has $n=180$ unrelated participants with a quantitative trait measured in trait units. A Gaussian linear model includes the same intercept, locus main effects, and adjustment covariates in both fits. The full model adds a four-degree-of-freedom interaction for two three-level genotype factors. The fitted residual sums of squares are $RSS_0=360$ for the reduced model and $RSS_1=324$ for the full model. Assume independent normal errors with a common variance, and that both models are identifiable and full rank.

For a Gaussian model with an unknown variance estimated by maximum likelihood, the profiled log-likelihood difference gives $\Lambda=n\log(RSS_0/RSS_1)$. Therefore,

\Lambda=180\log\left(\frac{360}{324}\right) =180\log(1.11111) =18.9649.

The null interaction has $q=4$ degrees of freedom, so the asymptotic p-value is

P(\chi^2_4\geq18.9649) =e^{-18.9649/2}\left(1+\frac{18.9649}{2}\right) \approx 0.00080.

For four degrees of freedom, the displayed survival-function expression is the chi-square tail probability. The result is strong evidence against the joint no-interaction hypothesis in this hypothetical analysis. It does not identify which genotype contrast is responsible, and it does not establish a causal biological interaction. In a normal linear model, a partial $F$ test can also provide an exact finite-sample comparison under the stated normal, common-variance assumptions; the chi-square likelihood-ratio calibration is asymptotic.

Worked Example: Binary Trait Likelihood-Ratio Test

Consider a hypothetical case-control sample of $600$ unrelated participants. The outcome is binary, and a logistic regression models the probability of the outcome using an intercept, both loci's main effects, and prespecified covariates. The full model adds all four interaction terms for the two three-level genotype factors. Assume independent Bernoulli outcomes conditional on the modeled covariates, no separation, identifiable parameters, and adequate sample size for the asymptotic approximation. Suppose the reduced and full model deviances are $512.4$ and $503.8$, respectively. For the same data and likelihood, the deviance is $-2$ times the log-likelihood up to a common saturated-model constant.

Thus the difference in deviance is the likelihood-ratio statistic:

\Lambda=D_0-D_1=512.4-503.8=8.600.

There are four interaction restrictions, so the reference distribution is $\chi^2_4$. Its tail probability is

P(\chi^2_4\geq8.600) =e^{-8.600/2}\left(1+\frac{8.600}{2}\right) =e^{-4.300}(5.300) \approx 0.07191.

At a prespecified significance level of $0.05$, this result does not provide sufficient evidence to reject the joint no-interaction hypothesis. It is not proof that interaction is absent: the data and model yield a p-value of about $0.0719$ under the specified test. The comparison concerns interaction on the logistic model's log-odds scale; it does not, by itself, test interaction on every possible outcome scale.

Score Test and Its Adjustment for Nuisance Parameters

The score test evaluates the interaction direction at the reduced-model estimate, without fitting the full model. Let $U_\gamma$ be the score vector for $\gamma$, the derivative of the log-likelihood with respect to the interaction parameters. Let the Fisher information matrix, evaluated under the null, be partitioned as

I= \begin{pmatrix} I_{\eta\eta} & I_{\eta\gamma}\\ I_{\gamma\eta} & I_{\gamma\gamma} \end{pmatrix}.

Because nuisance parameters are estimated rather than known, the interaction score must be adjusted for its information shared with $\eta$. The efficient information for $\gamma$, conditional on estimating $\eta$, is the Schur complement

I_{\gamma\cdot\eta} =I_{\gamma\gamma} -I_{\gamma\eta}I_{\eta\eta}^{-1}I_{\eta\gamma}.

The inverse $I_{\eta\eta}^{-1}$ exists under the relevant full-rank regularity condition. The second term removes the part of the interaction information explained by nuisance-parameter directions. This adjustment is essential when the interaction score and main-effect or covariate scores are correlated.

Theorem: Under the same regularity conditions used for the likelihood-ratio test, the score statistic for $H_0:\gamma=\mathbf{0}$ is $$ S=U_\gamma^{T}I_{\gamma\cdot\eta}^{-1}U_\gamma. $$ Under $H_0$, $S$ is asymptotically distributed as $\chi^2_q$. Both $U_\gamma$ and $I_{\gamma\cdot\eta}$ are evaluated at the reduced-model estimate of $\eta$ and at $\gamma=\mathbf{0}$.

Why the Score Statistic Has This Form

A short local argument explains the adjustment and the quadratic form. Near the null estimate, a small change $(\delta_\eta,\delta_\gamma)$ changes the log-likelihood approximately by the score contribution minus half the information-weighted squared change. Under the null fit, the nuisance score is zero. For a fixed small interaction change $\delta_\gamma$, choose the nuisance change that maximizes this local approximation. Differentiating with respect to $\delta_\eta$ gives

0=I_{\eta\eta}\delta_\eta+I_{\eta\gamma}\delta_\gamma, \qquad \delta_\eta=-I_{\eta\eta}^{-1}I_{\eta\gamma}\delta_\gamma.

Substituting this nuisance adjustment into the quadratic information term leaves the interaction information $I_{\gamma\cdot\eta}$. The resulting local log-likelihood gain, maximized over the interaction change, is one half of $U_\gamma^T I_{\gamma\cdot\eta}^{-1}U_\gamma$. Multiplying that gain by two yields $S$. Under the null, the appropriately standardized efficient score is asymptotically a $q$-dimensional standard normal vector; its squared length therefore has a $\chi^2_q$ distribution. This argument also shows why using $I_{\gamma\gamma}$ alone is generally wrong when nuisance parameters are estimated.

Worked Example: Score Test for Four Binary-Trait Interaction Terms

Suppose a hypothetical logistic regression is fitted under the null with both locus main effects and all adjustment covariates retained. For four interaction parameters, assume the null-fit efficient score vector and efficient information matrix are

U_\gamma= \begin{pmatrix}3\\-2\\4\\1\end{pmatrix}, \qquad I_{\gamma\cdot\eta}= \begin{pmatrix} 9&0&0&0\\ 0&8&0&0\\ 0&0&16&0\\ 0&0&0&5 \end{pmatrix}.

These are stipulated hypothetical outputs from the null fit; the diagonal information matrix is positive definite, as required for its inverse. Its inverse has diagonal entries $1/9,1/8,1/16,1/5$. The statistic is

S=U_\gamma^T I_{\gamma\cdot\eta}^{-1}U_\gamma =\frac{3^2}{9}+\frac{(-2)^2}{8}+\frac{4^2}{16}+\frac{1^2}{5} =1+0.500+1+0.200 =2.700.

Using the four-degree-of-freedom reference distribution,

P(\chi^2_4\geq2.700) =e^{-2.700/2}\left(1+\frac{2.700}{2}\right) =e^{-1.350}(2.350) \approx0.6092.

This score test gives little evidence against the joint null in the hypothetical example. Its interpretation is conditional on the stated logistic model and its assumptions. Unlike the likelihood-ratio test, it requires only the reduced-model fit, although the likelihood and information calculations for that fit must still be valid.

Conditions, Interpretation, and Common Pitfalls

Both tests rely on a regular, identifiable model and a suitable large-sample approximation. In particular, the interaction parameters must be estimable rather than aliased with other model terms; the likelihood must be sufficiently smooth near the null; and the null parameter must not lie on a problematic boundary. In quantitative-trait Gaussian regression, the usual likelihood model assumes independent normally distributed errors with a common variance. In logistic regression, the usual model assumes independent Bernoulli outcomes conditional on the covariates and genotype terms. Dependence, severe sparsity, or separation can make ordinary asymptotic calibrations unreliable; a nominal chi-square p-value should not be treated as automatically trustworthy in those settings.

The reduced and full models must differ only by the interaction terms being tested. Retaining both locus main effects respects the hierarchical model comparison used to ask whether interaction adds information beyond those effects. Any adjustment covariates included in one fit should also be included in the other. The degrees of freedom equal the rank of the tested restrictions, not necessarily the number of columns written down if some columns are redundant or unidentifiable.

Testing one term instead of the joint null. A p-value for one interaction coefficient does not test whether all interaction coefficients are zero. Use the joint statistic with the corresponding $q$ degrees of freedom.
Removing main effects from the reduced model. Comparing a model with interaction and main effects against one with neither asks a different question. For a no-interaction test, keep the main effects in both models.
Using the wrong reference distribution. The degrees of freedom come from the number of independent interaction restrictions. For example, a four-term interaction block is tested against $\chi^2_4$, not $\chi^2_1$.
Ignoring nuisance-parameter adjustment in the score test. The correct information is the efficient information $I_{\gamma\cdot\eta}$, not generally the unadjusted block $I_{\gamma\gamma}$.
Overinterpreting a p-value. A small p-value is evidence against the specified joint null under the model; it is not the probability that the null is true, proof of causation, or identification of a biological mechanism.

In practice, the likelihood-ratio test is convenient when both models can be fitted reliably, while the score test is useful when fitting the full model is costly or difficult. In regular large samples they test the same local null and often give similar results, but their numerical values need not match exactly in a finite sample. Report the model scale, the joint null, the statistic, its degrees of freedom, the p-value, and important assumptions so that the tested claim is clear.

Check Your Understanding

Answer each question using the joint-test framework in this tutorial.

In a model with two three-level genotype factors and a full interaction block, how many interaction degrees of freedom are tested if the design is full rank? Explain the count.
Write the likelihood-ratio statistic for reduced and full log-likelihoods, and state its asymptotic null distribution.
Why does the score test use $I_{\gamma\cdot\eta}$ rather than simply $I_{\gamma\gamma}$ when nuisance parameters are estimated?
A logistic model has reduced deviance $420.0$ and full deviance $413.0$, with four interaction restrictions. What is the likelihood-ratio statistic, and which chi-square distribution supplies its p-value?
Why is a test of one interaction coefficient not generally a test of the joint hypothesis of no two-locus interaction?