Interaction effects in regression: how to read them

An interaction effect in regression means that the slope of one predictor is not the same at every value of another predictor. You fit it by adding the product of the two predictors to the model, and you read it by asking what the first predictor’s slope is at each level of the second. The main-effect coefficients alone will mislead you.

This post reads one interaction model line by line. The data are simulated to look like a dose-by-supplement tooth-length study: a numeric dose, a supplement coded 0 or 1, and a built-in interaction. They are not a real study and not a published dataset. Every number is a deterministic engine output or comes from a base R snippet shown below, and the ones that do not come from the engine are marked as recomputed.

What an interaction is

A model without an interaction says: the outcome rises by the same amount per unit of dose whichever supplement an animal gets, and the supplement shifts every prediction by a fixed amount. The two predictors are additive. The prediction lines for the two supplements are parallel.

A model with an interaction drops that constraint. Write the outcome as y, the dose as d and the supplement code as s:

y = b0 + b1 d + b2 s + b3 (d x s)

The last term is the product of the two predictors, and b3 is the interaction coefficient. When s is 0 the product vanishes and the dose slope is b1. When s is 1 the dose slope is b1 + b3. So b3 is the difference between the two dose slopes, and “interaction” is another word for “the slopes differ”. Brambor, Clark and Golder (2006) make the same point for applied work: the effect of a variable in an interaction model is a function of the other variable, and a single coefficient cannot summarise it.

The product is an ordinary column. The engine’s linear regression takes a list of predictors, so the example data carry a column dose_x_supp that is the dose multiplied by the supplement code, and the model is fitted with dose_mg, supp_vc and dose_x_supp as three plain predictors. In R, the formula y ~ d * s builds the same column for you.

The data and the model

The simulated study has 60 animals: 10 in each of six cells made by crossing two supplements with three daily doses of 0.5, 1 and 2 mg. Each animal is measured once and is its own independent unit, and the outcome, tooth_length, is continuous, so a linear regression fits. The post on multiple linear regression and when to switch explains how that choice is made, and the linear regression method reference lists what the engine computes.

The data were generated with a real interaction on purpose, so the example has something to find. The engine fitted all 60 rows with 56 residual degrees of freedom and got R-squared = 0.7855.

The table holds four coefficient tests with no multiplicity correction. The one the question is about is the last row, so read the others for what they mean, not as four separate findings.

Two engine checks deserve a look on this fit, and the post on linear regression assumptions walks through them. The Shapiro-Wilk check on the residuals found no evidence against normality (p = 0.184). The engine also noted five observations above the 4/n Cook’s distance screening threshold (largest 0.146), which is a prompt to look at those rows, not a reason to drop them.

Reading the coefficient table

Go from the bottom row up.

The product row, dose_x_supp = 3.920. This is the difference between the two dose slopes: the dose slope is 3.920 units per mg steeper in the group coded 1 than in the group coded 0. Its 95% interval, 1.507 to 6.333, excludes zero, and the t statistic is 3.254 on 56 degrees of freedom with p = 0.0019. The interval is the result to report. It says the data are compatible with a difference in slopes anywhere from modest to large, and the p-value only says that a difference this far from zero would be unusual if the slopes were equal and the model held (Greenland and colleagues 2016 list the usual misreadings).

The dose_mg row, 6.314. This is not the effect of dose. It is the effect of dose when supp_vc is 0, so it is the dose slope in the group coded 0 only. Read it as an average over both groups and you get the wrong answer. This is the most common error with interaction models, and Westreich and Greenland (2013) describe the same trap for covariate coefficients in general.

The supp_vc row, -2.761. This is the difference between the groups when the dose is 0. A dose of 0 is outside the study: the lowest dose is 0.5 mg. The coefficient is a difference at a point where there are no animals, and its interval, -5.953 to 0.4315, runs from a clear deficit to a small gain. It is a construction that anchors the lines, not a finding about the supplement.

The intercept, 11.64. The predicted outcome at dose 0 in the group coded 0. Again the dose is 0, so it is a short extrapolation, not an observed case.

Here d is the example’s data, read with d <- read.csv("data/interaction-effects-regression.csv"). The same fit in plain R is one call, and dose_mg * supp_vc builds the product for you:

fit <- lm(tooth_length ~ dose_mg * supp_vc, data = d)
summary(fit)$coefficients

                  Estimate Std. Error   t value     Pr(>|t|)
(Intercept)      11.637814  1.1267685 10.328487 1.398849e-14
dose_mg           6.314045  0.8517569  7.412966 7.188272e-10
supp_vc          -2.760630  1.5934913 -1.732441 8.869989e-02
dose_mg:supp_vc   3.919785  1.2045662  3.254105 1.931289e-03

The numbers match the engine’s coefficient table. The product row is labelled dose_mg:supp_vc there, and it is the same column as dose_x_supp.

Simple slopes: the numbers that answer the question

The useful quantity is the dose slope at each supplement, called a simple slope (Preacher, Curran and Bauer 2006 describe computational tools for probing interactions this way).

The slope for the group coded 1 is the sum of two coefficients, 6.314 + 3.920 = 10.23. The engine does not print that sum or its interval, so the interval was recomputed in base R from the fitted model’s covariance matrix, and you can reproduce it:

b <- coef(fit)
V <- vcov(fit)
slope1 <- b["dose_mg"] + b["dose_mg:supp_vc"]
se1 <- sqrt(V["dose_mg", "dose_mg"] + V["dose_mg:supp_vc", "dose_mg:supp_vc"] +
            2 * V["dose_mg", "dose_mg:supp_vc"])
slope1 + c(-1, 1) * qt(0.975, df = 56) * se1

# supplement difference at a chosen dose: a contrast of coefficients, with its interval
for (x in c(0.5, 1, 2)) {
  L <- c(0, 0, 1, x)  # supp_vc + x * dose_mg:supp_vc
  est <- sum(L * b); se <- sqrt(drop(t(L) %*% V %*% L))
  print(c(dose = x, est = est, lo = est - qt(0.975, 56) * se, hi = est + qt(0.975, 56) * se))
}

The interval cannot be built by adding the two coefficients’ own intervals. The variance of a sum is both variances plus twice their covariance, so standard errors do not add, and here the covariance is negative (correlation -0.71, recomputed), which narrows the interval further. Adding the intervals would give 6.11 to 14.35 (recomputed), far wider than the correct 8.53 to 11.94.

The same simple-slope idea works the other way round. The supplement difference depends on the dose: it is -2.761 + 3.920 x dose, which is -0.80 at 0.5 mg, 1.16 at 1 mg and 5.08 at 2 mg. These three values are recomputed from the coefficients, and the fitted lines cross near 0.7 mg (2.761 / 3.920 = 0.704, recomputed). The difference at 0.5 mg has a 95% interval of -3.00 to 1.40 (recomputed), so the data do not settle which group is ahead at the lowest dose; at 2 mg the group coded 1 is clearly ahead (5.08, 95% CI 2.57 to 7.59, recomputed). The loop in the code above gives the interval at any dose.

Reading the interaction plot

An interaction plot draws one line per group, with the dose on the horizontal axis and the predicted outcome on the vertical axis. Each line’s slope is a simple slope, so the plot and the table above must agree.

nd <- expand.grid(dose_mg = c(0.5, 1, 2), supp_vc = c(0, 1))
nd$pred <- predict(fit, nd)
plot(range(nd$dose_mg), range(nd$pred), type = "n",
     xlab = "Dose (mg)", ylab = "Predicted tooth length")
for (g in 0:1) lines(nd$dose_mg[nd$supp_vc == g], nd$pred[nd$supp_vc == g], lty = g + 1)

Here the line for the group coded 0 climbs about 6.3 units per mg and the line for the group coded 1 climbs about 10.2, so the lines start close together at 0.5 mg and fan apart by 2 mg. The observed cell means, recomputed in base R, tell the same story:

Three cautions apply to any such plot. Parallel lines are consistent with no interaction, but with few animals per cell they do not prove it, so look at the interval for the difference in slopes. Lines that cross mean the sign of the group difference may depend on the dose, so a single “supplement effect” would mislead; here the fitted lines cross near 0.7 mg, but the difference at the low end is too uncertain to say its sign has flipped. And a line is only as good as its straight-line assumption: with three dose levels the data can reveal a bend (the middle dose sitting off the line through the other two), but they say nothing about the shape between the studied doses or beyond them.

Why the main effects looked different before

Compare the model without the product. Fitting tooth_length ~ dose_mg + supp_vc gives a dose slope of 8.27 and a supplement coefficient of +1.81. Both are averages over the interaction. In this balanced design, 8.27 is exactly the mean of the two simple slopes, 6.31 and 10.23, and +1.81 is the interaction model’s group difference at the average dose of 1.17 mg. It differs in sign from the -2.76 above because that coefficient is the difference at dose 0. Both are recomputed in base R.

Neither fit is wrong as arithmetic. They answer different questions. The additive model reports one blended dose slope and one blended supplement shift, and the comparison of the two nested models shows the product term earns its place: the F test for adding it gives F = 10.59 with p = 0.0019 (recomputed in base R), the same p-value as the t test on its row, because one added term makes the two tests equivalent. If you present only the additive coefficients when the slopes really differ, you hide the structure in the data.

with(d, tapply(tooth_length, list(dose_mg, supp_vc), mean))
add <- lm(tooth_length ~ dose_mg + supp_vc, data = d)
coef(add)
anova(add, fit)

Centring and the overlap between columns

Because the main effects are read at zero, centring the predictors changes what they mean. If the dose is centred at 1 mg, the supplement coefficient becomes the group difference at 1 mg, which is the recomputed 1.16 above, and the dose coefficient stays the slope in the group coded 0. If you build the product column by hand, rebuild it from the centred dose; centring the dose alone leaves the supplement coefficient at dose 0. The interaction coefficient does not change. Schielzeth (2010) discusses centring as a way to make coefficients interpretable.

Centring also affects overlap. The engine flagged a variance inflation factor above 5 for the product column, which is expected: a product shares information with its components. A high VIF on an interaction term is not by itself a reason to distrust the fit, and the post on multicollinearity and VIF explains how to judge it. Centring lowers the overlap, but it does not change the fitted lines or the predictions.

Common mistakes

• Reading a main effect as an average effect. In a model with a product term it is the effect where the other predictor is zero.
• Reporting only the interaction p-value. Give the difference in slopes with its interval, then the simple slopes.
• Dropping a non-significant interaction without looking at the interval. A wide interval means the data cannot say whether the slopes differ.
• Interpreting a zero that is not in the data. The supplement coefficient here sits at a dose nobody received.
• Testing many interactions. Each extra interaction is another test. Decide which ones are of interest before looking at the data, or say that the search was exploratory.
• Calling an interaction a cause. It describes how slopes differ in these data. Whether changing the supplement would change the dose response is a causal question that the model alone cannot answer.

A design with two categorical factors, and no numeric predictor, is usually analysed with a two-way analysis of variance instead; the guide to two-way ANOVA covers that case, and choosing a statistical test helps you decide between them. For how to read a single coefficient before any interaction is added, see the post on interpreting regression coefficients.

How to report an interaction

State the model, the difference in slopes with its interval, and the simple slopes. For this example:

In a simulated dose-by-supplement study of 60 animals (10 per cell), the dose slope of tooth length was 6.31 units per mg (95% CI 4.61 to 8.02) in the group coded 0 and 10.23 units per mg (95% CI 8.53 to 11.94, computed from the model’s covariance matrix) in the group coded 1. The difference in slopes was 3.92 units per mg (95% CI 1.51 to 6.33), so the slopes differ in these data.

State whether the model was planned, give the sample size per cell, and say which numbers were computed outside the main software output. Do not write that the interaction “proves” a mechanism, and do not describe a non-significant interaction as showing the slopes are equal.

What the example run recomputes

The reproduce.R file in the example folder reads the bundled data and refits the three-predictor model with lm(), so it recomputes the estimates, standard errors, t statistics and p-values in the coefficient table above and prints the sensitivity power calculation; confint() on the same model gives the intervals, and it does not rerun the engine’s normality, influence or VIF checks. The simple slope for the group coded 1, the supplement differences by dose, the cell means and the additive-model comparison are not part of that file. They were recomputed separately in base R from the same data, and the code that produces them is shown in this post.

# Results

- Test: Linear regression
- n: 60
- statistic: 7.413
- p-value: < 0.001
- effect size: R-squared = 0.7855
- 95% CI (dose_mg): [4.608, 8.02]
- seed: 20261025

## Assumption checks
- normality (Shapiro-Wilk): passed=TRUE (no evidence against normality)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=60: Cohen's f2 = 0.1949
- to detect a conventional medium effect: 77 total to detect f2=0.15
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Linear regression  |  engine 0.31.8  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261025)
.data_path <- "data/interaction-effects-regression.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("tooth_length", "dose_mg", "supp_vc", "dose_x_supp")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(tooth_length ~ dose_mg + supp_vc + dose_x_supp, data = data))
print(result)

# Sensitivity: minimum detectable effect at n=60 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.f2.test(u = 3, v = 56, sig.level = 0.05, power = 0.8)
print(power)