Tutorials › Statistical Genetics › Hierarchical Regression Models For Genotype Interactions

Home  /  Tutorials  /  Gene–Gene Interaction

Gene–Gene Interaction · Tutorial 820 of 989

Hierarchical Regression Models For Genotype Interactions

Build hierarchical regression models for two loci and interpret nested likelihood comparisons without confusing an interaction effect with omitted main effects.

Advanced 12 min read

What You'll Learn

  • Specify additive two-locus regression models with both locus main effects and their product term
  • Distinguish a main-effects model, a hierarchical interaction model, and an interaction-only model
  • Interpret interaction coefficients conditional on the model’s coding and scale
  • Identify when two fitted models are nested and comparable by likelihood
  • Calculate a nested-model likelihood difference and explain what it does and does not establish

Why Model Hierarchy Matters

A two-locus interaction model is not just a regression with a product term added wherever convenient. Its specification determines which genetic effects are being held fixed, which are allowed to vary, and what a likelihood comparison can tell us. A central modeling convention is hierarchy: when an interaction between two predictors is included, retain the corresponding lower-order main effects. This makes the interaction a departure from a model that already allows each locus to have its own effect.

Let \(G_1\) and \(G_2\) denote genotype dosages at two loci, coded as 0, 1, or 2 copies of a specified allele. The product \(G_1G_2\) is one possible interaction predictor. Its coefficient describes how the association with one dosage changes across values of the other dosage, on the model’s chosen scale. The interaction is statistical and coding-dependent; it does not by itself establish a biological mechanism.

Definition: A two-locus regression model is hierarchical with respect to a specified interaction if it includes the lower-order terms for both loci whenever it includes their interaction. With dosage coding and one product term, the model contains \(G_1\), \(G_2\), and \(G_1G_2\), together with an intercept and any prespecified covariates.

Writing the Nested Models

For a quantitative trait \(Y_i\), one common specification is a linear regression with normally distributed errors. Let \(C_i\) be a vector of covariates, such as age or ancestry variables. The main-effects model and the hierarchical interaction model are:

$$ Y_i = \beta_0 + C_i^\mathsf{T}\eta + \beta_1G_{1i} + \beta_2G_{2i} + \varepsilon_i, \qquad \varepsilon_i \sim N(0,\sigma^2), $$
$$ Y_i = \beta_0 + C_i^\mathsf{T}\eta + \beta_1G_{1i} + \beta_2G_{2i} + \gamma G_{1i}G_{2i} + \varepsilon_i, \qquad \varepsilon_i \sim N(0,\sigma^2). $$

Here, \(\gamma\) is the change in the \(G_1\) slope for a one-unit increase in \(G_2\), and equivalently the change in the \(G_2\) slope for a one-unit increase in \(G_1\). This interpretation follows from the conditional mean: the \(G_1\) slope at \(G_2=g_2\) is \(\beta_1+\gamma g_2\). The coefficient is measured in units of \(Y\) per dosage unit squared for the product term.

For a binary trait, the same predictor structure can be used in logistic regression. The left side then models the log odds rather than the mean on the original outcome scale:

$$ \log\left(\frac{P(Y_i=1\mid G_{1i},G_{2i},C_i)} {1-P(Y_i=1\mid G_{1i},G_{2i},C_i)}\right) = \beta_0 + C_i^\mathsf{T}\eta + \beta_1G_{1i} + \beta_2G_{2i} + \gamma G_{1i}G_{2i}. $$

The product coefficient is now an interaction on the log-odds scale. It should not be described as a constant change in risk. Predicted probabilities can be calculated from the fitted linear predictor, but the interaction’s interpretation depends on the link and coding used.

Model assumptions: Likelihood comparisons require that the candidate models use the same outcome observations and the same likelihood framework. The regression likelihood must also be appropriate for the outcome and sampling design. For the Gaussian model above, this includes independent errors with the stated normal distribution and common variance; for logistic regression, the binary outcomes are modeled as independent Bernoulli observations conditional on predictors. Model terms must be estimable: exact collinearity, empty genotype combinations, or other rank deficiencies can prevent unique estimation.

Worked Example: Specifying a Hierarchical Binary-Trait Model

Worked Example: Specifying a Hierarchical Binary-Trait Model

Suppose a hypothetical analysis uses allele dosages \(G_1\) and \(G_2\) and a binary outcome \(Y\). A fitted logistic model has intercept \(-2.00\), locus coefficients \(0.30\) and \(0.20\), and product coefficient \(0.25\). For clarity, assume there are no additional covariates. What does the model predict for a person with \(G_1=1\) and \(G_2=2\)?

The hierarchical model is \(\operatorname{logit}\{P(Y=1\mid G_1,G_2)\}=-2.00+0.30G_1+0.20G_2+0.25G_1G_2\). Both main effects remain in the equation alongside the product term.

Substituting the dosages gives a fitted log odds of \(-2.00+(0.30)(1)+(0.20)(2)+(0.25)(1)(2)=-0.80\). The corresponding odds are \(\exp(-0.80)=0.4493\), and the probability is \(\exp(-0.80)/\{1+\exp(-0.80)\}=0.3100\).

Thus, under this hypothetical fitted model, the predicted probability is about 0.310 for the specified dosage combination. It is not a causal effect estimate. The product coefficient is a log-odds-scale departure from additivity in the two dosage terms, conditional on the model.

Nested Models and Likelihood Comparison

A model is nested within another when the smaller model can be obtained from the larger one by fixing one or more of its parameters to specified values, usually zero. For example, the main-effects model above is nested in the hierarchical interaction model by setting \(\gamma=0\). This makes their comparison directly relevant to the additional product term, provided both models use the same outcome, observations, covariates, and likelihood.

Let \(\ell_{\text{small}}\) and \(\ell_{\text{large}}\) be the maximized log likelihoods for nested models. The larger model contains all parameter settings available to the smaller one, so its maximum likelihood cannot be lower. This gives a simple justification for the likelihood comparison:

$$ \{\text{parameter values in the small model}\} \subseteq \{\text{parameter values in the large model}\} \quad\Longrightarrow\quad \ell_{\text{large}} \geq \ell_{\text{small}}. $$

The likelihood-ratio statistic is the change in deviance, \(2(\ell_{\text{large}}-\ell_{\text{small}})\). Its size describes how much the maximized likelihood improves when the additional parameters are allowed. Under regularity conditions, its null distribution is approximately chi-squared with degrees of freedom equal to the difference in the number of independently estimated parameters. These conditions include a correctly specified likelihood, identifiable parameters, and a sufficiently large sample; boundary parameters or sparse cells can invalidate the usual approximation.

This comparison is only meaningful when the models are genuinely nested and fitted on the same data. A larger log likelihood alone is not enough to justify a model with more terms: adding parameters cannot reduce the maximum likelihood, even when the apparent gain is due to sampling noise. Likelihood comparisons also do not decide whether a genetic interaction is biologically meaningful.

Worked Example: Comparing Main Effects With a Hierarchical Interaction

In a hypothetical Gaussian-trait analysis of the same 500 individuals, the main-effects model has maximized log likelihood \(-202.0\), and the model adding \(G_1G_2\) has maximized log likelihood \(-198.5\). All covariates and other modeling choices are identical. Calculate and interpret the likelihood improvement.

The main-effects model is nested within the hierarchical model because it is recovered by setting \(\gamma=0\). The larger model adds one parameter. The change in deviance is \(2\{(-198.5)-(-202.0)\}=2(3.5)=7.00\).

The fitted interaction model has the higher maximized log likelihood, as it must for nested models, and the improvement corresponds to a 7.00-unit decrease in deviance for one additional parameter. This is a likelihood comparison of the hierarchical models. It does not compare an interaction against a model that omits the locus main effects, nor does this arithmetic alone establish biological interaction.

What Changes Under Alternative Hierarchies?

An interaction-only model includes \(G_1G_2\) but excludes \(G_1\) and \(G_2\). That is a different restriction: it assumes that both locus main effects are zero on the model’s linear-predictor scale while allowing a product effect. It is not the usual hierarchical interaction model. A better fit for this model than for an intercept-only model would not show that the product term captures interaction independently of the missing main effects.

The hierarchy principle applies to the terms actually used to represent the predictors. If a locus is represented by multiple genotype indicators or includes a dominance term, an interaction involving that representation should retain its corresponding lower-order terms. In a categorical two-locus model, the full interaction allows the effect of one locus to vary across categories of the other; its number of interaction parameters depends on the number of levels and the reference coding. It is therefore important to state the coding and count the added parameters rather than assume every two-locus interaction has one degree of freedom.

Worked Example: Why an Interaction-Only Comparison Is Different

Consider a hypothetical nested comparison using the same binary-trait observations and logistic likelihood. The intercept-only model has log likelihood \(-120.0\), an interaction-only model containing \(G_1G_2\) has log likelihood \(-115.0\), a main-effects model containing \(G_1\) and \(G_2\) has log likelihood \(-110.0\), and the hierarchical model containing both main effects and the product has log likelihood \(-108.5\).

Comparing the hierarchical model with the main-effects model isolates the addition of the product term: the deviance change is \(2\{(-108.5)-(-110.0)\}=3.00\), for one added parameter. Comparing the hierarchical model with the interaction-only model instead adds two main-effect parameters, giving \(2\{(-108.5)-(-115.0)\}=13.0\). That second comparison asks about adding both omitted main effects, not about adding an interaction to a main-effects model.

The interaction-only model has a higher log likelihood than the intercept-only model, but this fact does not separate interaction from main effects. The main-effects model has a still higher likelihood, and the hierarchical model has the highest likelihood among these candidates. The models have different restrictions, so their comparison must be interpreted in terms of the parameters released, not by labeling any likelihood gain “epistasis.”

A Practical Comparison Workflow

Before reporting a model comparison, make the hierarchy and the target contrast explicit. The following sequence prevents a common mismatch between the scientific question and the models being compared.

1
Define the outcome model and coding.
State whether the outcome is quantitative or binary, specify the link or error distribution, and define how each locus and interaction is represented.
2
Fit the main-effects model.
Include both locus terms and the same covariates intended for the interaction model.
3
Add the interaction hierarchically.
Retain both locus main effects while adding the prespecified interaction terms.
4
Check comparability and estimability.
Confirm the models use the same observations and likelihood, that the smaller model is nested, and that the added terms can be estimated from the observed genotype combinations.
5
Describe the comparison precisely.
Report the log likelihoods, added parameter count, and likelihood improvement, then state which restrictions are being compared. Do not interpret the improvement as proof of a biological mechanism.

Common Mistakes and Interpretation

A frequent error is to compare the full hierarchical model with an interaction-only model and call the difference an interaction test. Those models differ in their main effects as well as their overall structure. To assess what the product terms add beyond the locus effects, compare the hierarchical model with the otherwise identical main-effects model.

Another error is to remove a main effect because its individual coefficient is not statistically significant, while retaining the interaction containing that locus. This changes the model’s assumptions and can make the product coefficient difficult to interpret. Hierarchy is a modeling principle, not a claim that every main-effect coefficient must be nonzero or individually significant.

Finally, a product term is tied to its coding and scale. A dosage product in a linear model describes departure from additivity in the conditional mean under that coding; in logistic regression it describes departure from additivity in log odds. Neither interpretation automatically transfers to the risk scale or identifies a biological interaction.

Key takeaway: To assess a two-locus interaction in a hierarchical regression, compare a model with both locus main effects to the same model with the prespecified interaction terms added. The likelihood comparison concerns those added terms only when the models are nested, use the same data and likelihood, and retain the lower-order effects.

Check Your Understanding

Answer each question by identifying the model restrictions involved.

  1. Write the dosage-coded mean model for a quantitative trait that includes both locus main effects and their product interaction.
  2. In a logistic regression, what scale does the coefficient of \(G_1G_2\) describe?
  3. Why must the main-effects and interaction models use the same observations for a direct likelihood comparison?
  4. An interaction-only model fits better than an intercept-only model. Why does that not isolate the interaction beyond the locus main effects?
  5. If a full categorical two-locus model adds four interaction parameters to a main-effects model, what parameter-count difference should be used when describing their nested comparison?

Previous

Next

Back to Tutorials
Clinical Biostats