Answer

Which interaction terms belong in a regression or mixed model?

Inspired by a question on Cross Validated · November 3, 2026

regressionmultilevel modelsmultiple comparisonsmodel selection

The short answer

Say a model has four predictors, A, B, C and D, and you think A, B and C interact. Two separate rules are involved, and they are often mixed up.

• Keep everything underneath an interaction you include. A three-way term A × B × C needs the three main effects and the three two-way interactions among A, B and C. This is the principle of marginality (also called the hierarchy principle). In R, writing A * B * C adds all seven terms for you.
• You do not need every other interaction of the same order. Nothing in the mathematics requires the three other three-way terms, each of which involves D, just because A × B × C is present. D can enter as a main effect only. A model is not "incomplete" because it leaves out interactions you had no reason to expect.
• What decides whether the result can be trusted is how the term was chosen. An interaction specified before you saw the outcome can be tested at face value. One chosen because it looked striking in exploratory plots of the same data cannot, because its p-value no longer means what it claims.

None of this changes in a mixed-effects model. Random effects describe how observations cluster (by subject, school or site); the rules about which fixed-effect terms to include are the same as in ordinary regression. In lme4 the model would be written y ~ A * B * C + D + (1 | subject). See our page on fixed vs random effects for that side of the model.

Why the lower-order terms must stay

An interaction coefficient measures how one predictor's slope changes with another, relative to the main effects and lower interactions in the model. Remove those lower terms and the interaction is forced to soak up whatever they would have explained, so its estimate is biased and hard to interpret.

There is also a practical test of whether a model makes sense: its conclusions should not depend on arbitrary choices such as where zero sits on a scale. With the lower-order terms in place, shifting a continuous predictor (for example, centering it) changes the main-effect coefficients but leaves the highest-order interaction and the model's predictions unchanged. Drop them and the interaction estimate, its test and the fitted values all change when you merely move the origin. The same issue arises with categorical predictors, where the answer would depend on how the factor is coded.

Keeping the lower-order terms does not mean they must be significant, and it does not mean you should interpret them on their own: once a predictor is part of an interaction, its main effect is only the effect at zero (or at the reference level) of the other variables.

See it in R and Python

The real risk is selection, so the code simulates it. Each of 2,000 datasets has 100 observations and four standard-normal predictors, and the outcome depends only on their main effects: there are no true interactions of any kind. For each of the four possible three-way interactions, the code fits a regression containing that term, its lower-order terms and the remaining main effect, and records the interaction's p-value. It then compares three habits: testing A × B × C because it was planned, testing whichever of the four looked best, and testing the best one with a Bonferroni correction for the four candidates.

# No true interactions: y depends only on the four main effects
set.seed(2026)
n_sims <- 2000
n <- 100
triples <- list(c("A", "B", "C"), c("A", "B", "D"),
                c("A", "C", "D"), c("B", "C", "D"))

# p-value of the three-way term, with all its lower-order terms kept
p_three_way <- function(dat, v) {
  others <- setdiff(c("A", "B", "C", "D"), v)
  f <- as.formula(paste("y ~", paste(v, collapse = " * "), "+",
                        paste(others, collapse = " + ")))
  fit <- lm(f, data = dat)
  coef(summary(fit))[paste(v, collapse = ":"), "Pr(>|t|)"]
}

p_mat <- matrix(NA, n_sims, 4)
for (s in 1:n_sims) {
  dat <- data.frame(A = rnorm(n), B = rnorm(n), C = rnorm(n), D = rnorm(n))
  dat$y <- with(dat, 0.5 * A + 0.5 * B + 0.5 * C + 0.5 * D + rnorm(n))
  p_mat[s, ] <- sapply(triples, function(v) p_three_way(dat, v))
}

# Pre-specified: test A:B:C only
cat("False-positive rate, A:B:C chosen in advance:",
    round(mean(p_mat[, 1] < 0.05), 3), "\n")
# Data-driven: look at all four, keep the most 'interesting' one
cat("False-positive rate, best of four three-way terms:",
    round(mean(apply(p_mat, 1, min) < 0.05), 3), "\n")
# Bonferroni over the four candidate terms
cat("Best of four with Bonferroni (p < .0125):",
    round(mean(apply(p_mat, 1, min) < 0.05 / 4), 3), "\n")

# No true interactions: y depends only on the four main effects
import itertools
import numpy as np
from scipy import stats

rng = np.random.default_rng(2026)
n_sims, n = 2000, 100
triples = list(itertools.combinations(range(4), 3))  # ABC, ABD, ACD, BCD

def p_three_way(X, y, v):
    # Intercept, all four main effects, the three two-way terms, the three-way term
    i, j, k = v
    cols = [np.ones(n)] + [X[:, m] for m in range(4)]
    cols += [X[:, i] * X[:, j], X[:, i] * X[:, k], X[:, j] * X[:, k]]
    cols.append(X[:, i] * X[:, j] * X[:, k])
    M = np.column_stack(cols)
    beta, *_ = np.linalg.lstsq(M, y, rcond=None)
    resid = y - M @ beta
    df = n - M.shape[1]
    sigma2 = resid @ resid / df
    se = np.sqrt(sigma2 * np.linalg.inv(M.T @ M)[-1, -1])
    return 2 * stats.t.sf(abs(beta[-1] / se), df)

p_mat = np.empty((n_sims, 4))
for s in range(n_sims):
    X = rng.normal(size=(n, 4))
    y = 0.5 * X.sum(axis=1) + rng.normal(size=n)
    p_mat[s] = [p_three_way(X, y, v) for v in triples]

print("False-positive rate, A:B:C chosen in advance:",
      round(np.mean(p_mat[:, 0] < 0.05), 3))
print("False-positive rate, best of four three-way terms:",
      round(np.mean(p_mat.min(axis=1) < 0.05), 3))
print("Best of four with Bonferroni (p < .0125):",
      round(np.mean(p_mat.min(axis=1) < 0.05 / 4), 3))

False-positive rate, A:B:C chosen in advance: 0.056 
False-positive rate, best of four three-way terms: 0.186 
Best of four with Bonferroni (p < .0125): 0.053

These figures come from one seeded run of the R code. Python's random numbers differ from R's, so its rates differ slightly (0.056, 0.176 and 0.048) but show the same pattern.

What the simulation shows

• A planned interaction behaves as advertised. Testing A × B × C because you meant to gave a false-positive rate of 5.6%, within simulation noise of the nominal 5%.
• Choosing the most interesting of four more than triples the error. Keeping whichever three-way term looked best produced a significant result in 18.6% of datasets that contained no interaction at all. That is close to the 1 − 0.95⁴ ≈ 18.5% you would expect from four independent tests.
• Choosing by eye is the same act. Picking terms from exploratory plots is selection too, only less visible, and the number of interactions you implicitly scanned is usually larger than four.
• Accounting for the search restores the error rate. Applying a Bonferroni correction across the four candidates brought the rate back to 5.3%. The catch is that you can only correct for a search you can count.

What to do in practice

• If the interaction is a hypothesis, state it in advance (ideally in a preregistration or analysis plan), include it with its lower-order terms, and test it at face value. Other interactions do not need to be added.
• If the study is exploratory, say so and treat the interaction as a lead, not a finding. Report how many interactions you looked at, adjust for that number or report all of them, and present the result as something to confirm in new data.
• Split or replicate when you can. Choosing terms on one part of the data and testing them on a held-out part gives an honest p-value for the chosen model, at the price of a smaller sample for each step.
• Do not let significance tests pick terms one at a time. Adding or dropping interactions according to their p-values is a form of stepwise selection, with the same optimistic standard errors and unstable choices.
• Mind the sample size. Three-way interactions are estimated far less precisely than main effects. A study powered for main effects usually cannot detect a modest three-way interaction, so a nonsignificant one is weak evidence that it is absent.

For help choosing an analysis in the first place, try our test chooser, and the DASS blog has more on regression and model building.

How to report a three-way interaction in APA style (7th edition)

In the method section, say which interactions were planned and how any others were chosen, so readers can judge the p-values:

"The model included the A × B × C interaction, specified in advance, together with all lower-order terms and a main effect of D. Three further three-way interactions were examined in exploratory analyses, and their p-values were Bonferroni-adjusted for four tests."

In the results, give the coefficient with its uncertainty, replacing the placeholders with your own figures: "The A × B × C interaction was significant, b = x.xx, SE = x.xx, 95% CI [LL, UL], t(df) = x.xx, p = .xxx." For a mixed model, name the software and the degrees-of-freedom method (for example Satterthwaite) used for the t test. Put the full set of fixed effects in a table, and follow a significant interaction with simple slopes or a plot rather than interpreting the coefficient alone.

Related tools and guides

• APA 7 formatter for regression
• Multilevel power and MDES calculator
• Fixed vs random effects in mixed models
• Centering vs standardizing predictors
• Why stepwise regression is a bad idea
• Bonferroni correction: when should you use it?
• Principle of marginality (Wikipedia)

More answered questions

• Bagging, boosting and stacking: how do these ensemble methods differ?
• Training, validation and test sets: what is each one for?
• Fixed vs random effects in mixed models: what is the difference?

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.