AdvancedIndependently reviewed, not yet spot-checked by a human

Interaction terms and subgroup analysis

What each of the four coefficients in an interaction model means, why exp(interaction) is a ratio of ratios, how to recover the stratum estimates from one model instead of splitting the data, and the mistake a subgroup forest plot invites: reading a difference in significance as a significant difference.

What an interaction term is for

The primary result of a trial is an average effect. Having read that number, the next clinical question is almost always: who does this drug work best for? What about older patients? What about the ones who were high risk to begin with?

That question has a formal name — interaction, or effect modification in the epidemiological literature. It asks whether the treatment effect changes with some characteristic of the patient.

The way to answer it is not to cut the data into pieces, fit each piece separately, and see whose p-value is smaller. That is the error this page exists to take apart, and it is one of the most common statistical misreadings in the clinical literature:

“significant in women, not significant in men” = “the treatment effect differs between women and men”

That equals sign does not hold. The left-hand side is two tests, each against a null of its own. The right-hand side is one test of the two against each other. The strata differ in sample size and in event count, and that alone is enough to push one of them across 0.05 while the other stays above it — even when nothing underneath varies by stratum at all. This page puts both readings of the same data side by side so the gap is visible at its actual size.

The example on this page

This page uses medicaldata::indo_rct, the trial behind the randomised controlled trial chapter: rectal indomethacin to prevent post-ERCP pancreatitis (PEP). There are 602 participants and 79 PEP events.

The overall effect, from an unadjusted logistic regression: 27/295 (9.2%) in the indomethacin arm against 52/307 (16.9%) on placebo, OR 0.49 (0.30–0.81), p 0.005.

This dataset was chosen for this page because it ships with a whole set of pre-specified subgroup variables, and on several of them the stratum estimates look far apart while the interaction test gives no support at all for a difference. That combination is the hardest thing on this topic to teach and the most important thing to see.

library(medicaldata)
data(indo_rct, package = "medicaldata")

d <- indo_rct
d$y    <- as.integer(d$outcome == "1_yes")          # 1 = PEP occurred
d$tx   <- as.integer(d$rx == "1_indomethacin")      # 1 = indomethacin
d$male <- as.integer(d$gender == "2_male")          # 0 = female (reference)

# One model with an interaction, four coefficients
fit <- glm(y ~ tx * male, data = d, family = binomial)
summary(fit)

# Stratum ORs: recover them from the one model, do not split and refit
exp(coef(fit)["tx"])                                 # female stratum
exp(coef(fit)["tx"] + coef(fit)["tx:male"])          # male stratum

# The male stratum's CI needs the COVARIANCE of the two coefficients,
# not just the two standard errors added up
V  <- vcov(fit)
b  <- unname(coef(fit)["tx"] + coef(fit)["tx:male"])
se <- sqrt(V["tx", "tx"] + V["tx:male", "tx:male"] + 2 * V["tx", "tx:male"])
exp(c(OR = b, lcl = b - 1.96 * se, ucl = b + 1.96 * se))

# P for interaction: likelihood ratio test, not Wald
fit0 <- glm(y ~ tx + male, data = d, family = binomial)
anova(fit0, fit, test = "LRT")

# The same thing for every subgroup variable
for (v in c("gender", "sod", "pep", "psphinc", "precut", "train", "bsphinc")) {
  f  <- droplevels(d[[v]])
  m0 <- glm(y ~ tx + f, data = d, family = binomial)
  m1 <- glm(y ~ tx * f, data = d, family = binomial)
  cat(v, anova(m0, m1, test = "LRT")[["Pr(>Chi)"]][2], "\n")
}

Verified against R 4.6.0 with medicaldata 0.2.0. Every interaction test uses anova(..., test = "LRT"); the last section says why.

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

d = sm.datasets.get_rdataset("indo_rct", "medicaldata").data
d["y"]    = (d["outcome"] == "1_yes").astype(int)
d["tx"]   = (d["rx"] == "1_indomethacin").astype(int)
d["male"] = (d["gender"] == "2_male").astype(int)

fit = smf.logit("y ~ tx * male", data=d).fit()
print(fit.summary())

# Stratum ORs
b_tx, b_int = fit.params["tx"], fit.params["tx:male"]
print(np.exp(b_tx), np.exp(b_tx + b_int))

# The male stratum's CI: a contrast vector pulls in the covariance for you
c = np.zeros(len(fit.params))
c[list(fit.params.index).index("tx")] = 1
c[list(fit.params.index).index("tx:male")] = 1
se = np.sqrt(c @ fit.cov_params().values @ c)
print(np.exp([b_tx + b_int - 1.96 * se, b_tx + b_int + 1.96 * se]))

# P for interaction: build the LRT by hand
fit0 = smf.logit("y ~ tx + male", data=d).fit()
lr = 2 * (fit.llf - fit0.llf)
print(lr, stats.chi2.sf(lr, df=1))

statsmodels has no ready-made LRT interface, so the code below subtracts the two log-likelihoods and looks the statistic up in a chi-square distribution.

What the four coefficients mean

Written out, glm(y ~ tx * male) is this model:

\log \frac{p}{1-p} = \beta_0 + \beta_1 \cdot \text{tx} + \beta_2 \cdot \text{male} + \beta_3 \cdot (\text{tx} \times \text{male})

and these are the four coefficients it produces:

Term	What it is	Estimate (log scale)	SE	Exponentiated	p
`(Intercept)`	Log odds of PEP in the reference cell: placebo, female	-1.5569	0.1678	0.2108	< 0.001
`tx`	Log odds ratio for indomethacin vs placebo AMONG FEMALES (the reference stratum)	-0.7897	0.2880	0.4540	0.006
`male`	Log odds ratio for male vs female AMONG PLACEBO (the reference treatment arm)	-0.1777	0.3986	0.8372	0.656
`tx:male`	Difference between the male and the female log odds ratio for treatment	0.3927	0.6111	1.4809	0.521

One row at a time:

The intercept. Set both tx and male to 0 and the intercept is all that is left. Exponentiated, 0.2108 is the odds of PEP among women on placebo — and you can check it by hand, because that cell holds 43 events and 204 non-events, and one divided by the other is exactly this number.
The tx coefficient is the treatment effect in women, not the treatment effect overall. Exponentiated it is 0.4540, the OR in the female stratum. This is the single most misread cell on the page.
The male coefficient is the male-versus-female contrast among patients on placebo, not the overall risk in men. Exponentiated it is 0.8372.
Only tx:male is the interaction: the treatment log OR in men minus the treatment log OR in women. It is the difference between two effects, not either effect itself.

Adding an interaction term changes what the main effects mean.

Without the interaction term, the tx coefficient is “the average treatment effect adjusted for male”. Add tx:male and the very same tx coefficient becomes “the treatment effect when male = 0” — the effect in the reference stratum.

Two practical consequences:

Change the reference level and the main effect changes, while the model fit, the fitted values and the p-value for the interaction all stay exactly where they were. Someone reporting that “the main effect shrank once we added the interaction” has usually not discovered anything; they have changed what the coefficient is conditional on.
A continuous modifier makes this worse. In glm(y ~ tx * age) the tx coefficient is the treatment effect at age zero, and there is nobody aged zero in the data. Centre the modifier — subtract its mean, or a clinically meaningful reference value — before fitting, and the main effect becomes readable again.

exp(interaction) is a ratio of ratios

Expand the definition and the treatment OR in each stratum is:

\mathrm{OR}_{\text{female}} = \exp(\beta_1), \qquad \mathrm{OR}_{\text{male}} = \exp(\beta_1 + \beta_3)

Divide one by the other and $\beta_1$ cancels:

\frac{\mathrm{OR}_{\text{male}}}{\mathrm{OR}_{\text{female}}} = \frac{\exp(\beta_1 + \beta_3)}{\exp(\beta_1)} = \exp(\beta_3)

So exponentiating the interaction term gives a ratio of odds ratios. The script that makes the figures computes both sides and compares them:

Route	Value
Male stratum OR ÷ female stratum OR = 0.6723 ÷ 0.4540	1.480909
exp(interaction coefficient) = exp(0.3927)	1.480909
Absolute difference between the two	0

The identity is algebra, not a coincidence in this dataset — the script holds it with stopifnot(), so nothing is produced if it ever fails.

Knowing it is a ratio of ratios settles how to read it. Its null value is 1, not 0; it is tested and given a confidence interval on the log scale; and like any other ratio it has to be read with that interval attached. Here it is 1.48 (0.45–4.91).

That interval has to be converted back into ORs before its width means anything, because it is a ratio and not an effect. Multiply each end by the female-stratum OR of 0.45 and you get the range of male-stratum ORs compatible with these data.

The lower bound, 0.447, corresponds to a male-stratum OR of 0.20 — lower than the 0.45 seen in women, meaning men would be getting more protection, not less. The upper bound, 4.906, corresponds to a male-stratum OR of 2.23 — past 1, meaning treatment could be tending towards harm in men.

One interval holding both “men benefit more” and “men are harmed” is an interval that rules out nothing. No evidence was detected that the treatment effect differs by sex.

A ratio of ratios is interaction on the multiplicative scale, and only on that scale.

The same set of effects can answer differently on the multiplicative scale (OR, RR) and the additive scale (risk difference), and that is not a contradiction. As long as the two strata have different baseline risks, a relative effect that is constant across strata forces the absolute benefit to differ across them, and the reverse holds too. The single exception is a treatment effect that is exactly null.

The two scales are asking different questions: the multiplicative scale asks whether the relative effect is constant, the additive scale asks whether the absolute benefit is.

Clinical decisions — whether to give this drug to this patient — turn on absolute benefit. A P for interaction on the OR scale is therefore not enough on its own to support a claim that some group is more worth treating. Answering that needs the absolute risks and risk differences in each stratum; see Marginal estimates and G-computation.

Recovering the stratum estimates from the model

With the four coefficients in hand, the stratum ORs follow straight from the definition: the reference stratum is exp(second coefficient), the other is exp(second coefficient + interaction).

The confidence intervals are less direct. The male stratum’s log OR is a sum of two coefficients, so its variance is

\mathrm{Var}(\beta_1 + \beta_3) = \mathrm{Var}(\beta_1) + \mathrm{Var}(\beta_3) + 2\,\mathrm{Cov}(\beta_1, \beta_3)

Drop the covariance term and the interval comes out wrong, in a direction set by its sign. In this model it is -0.083, so dropping it would make the interval too wide. vcov(fit) in the R code above is what fetches it.

So why not just split the data and fit each stratum separately?

In this particular example you would get the same numbers. The script computes both routes and compares them with stopifnot() — on the ORs and on the confidence intervals: the ORs agree to eight decimal places, and the largest gap between the two sets of interval bounds is 4.8e-5.

The reason is that this model is saturated: two binary variables give four cells, the model has four parameters, and it is therefore a rewriting of the four cell counts with nothing constrained.

That stops being true the moment any covariate enters. glm(y ~ tx * male + age + risk) assumes the effects of age and risk are shared between men and women; splitting the data lets each stratum estimate its own. The first uses the whole dataset to estimate the covariate effects and is more efficient; the second assumes less and pays for that in sample size.

Three practical reasons make “recover it from the model” the default:

P for interaction can only come from a model that contains the interaction. Splitting the data cannot produce the test at all.
Covariate effects are estimated once, from everybody, so the stratum estimates are more precise.
Splitting shrinks the sample twice over, and subgroup analysis is already operating at the edge of its power.

How to read a subgroup forest plot

This trial has 7 pre-specified binary subgroup variables. Plotting the treatment OR in every stratum alongside each variable’s P for interaction gives the figure that appears in the paper:

Subgroup forest plot. The top row is the overall effect, OR 0.49 (0.30–0.81), drawn as a red diamond sitting to the left of OR = 1. Below it are 7 subgroup variables, each with two stratum rows, 14 rows in all. Each row carries the events over patients for the indomethacin arm and the placebo arm on the left, and that stratum's OR with its 95% confidence interval on the right; the P for interaction appears only on the variable heading row, in the far right column. Two vertical lines cross the plot: a solid one at OR = 1 and a dashed red one at the overall effect. All 14 stratum point estimates fall to the left of OR = 1, and 6 of their confidence intervals cross 1. The interval on the "Yes" row for precut sphincterotomy ends in an arrowhead on the left, meaning it runs off the axis. — The axis is clipped between OR 0.05 and 4; an interval that runs past either end is drawn with an arrowhead, so the true interval is longer than what you see. The right-hand column holds each variable's P for interaction, not the p-values of the individual strata — what goes in that column decides what conclusion a reader walks away with.Plotting script `figures/scripts/B2-06-interaction.R`

The same numbers as a table:

Subgroup	Stratum	n	Events	OR	95% CI	Stratum p	P for interaction
Sex	Female	476	63	0.45	0.26–0.80	0.006	0.52
	Male	126	16	0.67	0.23–1.93	0.461
Sphincter of Oddi dysfunction	No	107	16	0.37	0.11–1.24	0.108	0.60
	Yes	495	63	0.53	0.31–0.91	0.022
Previous post-ERCP pancreatitis	No	506	56	0.54	0.30–0.96	0.037	0.49
	Yes	96	23	0.36	0.13–0.98	0.046
Pancreatic sphincterotomy	No	259	32	0.33	0.14–0.77	0.010	0.23
	Yes	343	47	0.63	0.33–1.17	0.142
Precut sphincterotomy	No	570	74	0.52	0.31–0.86	0.011	0.49
	Yes	32	5	0.23	0.02–2.36	0.217
Trainee involvement	No	319	31	0.49	0.22–1.08	0.076	0.95
	Yes	283	48	0.47	0.25–0.90	0.023
Biliary sphincterotomy	No	258	33	0.31	0.13–0.72	0.006	0.16
	Yes	344	46	0.66	0.35–1.23	0.192

The order to read the figure in:

Start with the overall row. It is the reference point for every stratum. Here it is OR 0.49.
Then each variable’s P for interaction, which is what says whether there is evidence that the effect varies with that characteristic. Across the 8 subgroup variables on this page, P for interaction runs from 0.16 to 0.95; the number below 0.05 is 0.
Only then the stratum estimates and intervals — and what you are looking for is whether they sit consistently near the overall effect, not which of them happens to clear 1.
While you are there, read the n and the event count in every row. The “Yes” stratum for precut sphincterotomy holds 32 patients and 5 events, and its interval of 0.02–2.36 spans everything from near-complete prevention to more than doubled risk. A stratum estimate like that carries almost no information, and on the plot it looks exactly as prominent as every other row.

A difference in significance is not a significant difference

Now take that second-to-last column seriously for a moment. Of the 7 binary subgroup variables, 6 produce the pattern “p below 0.05 in one stratum, above it in the other”. Judged on stratum p-values alone, this one trial would support 6 different stories of the form “the drug only works in such-and-such patients”.

Not one of the interaction tests on the same data reaches significance. This figure puts the two readings next to each other:

Two panels side by side. Panel A on the left shows three subgroup variables (sex, pancreatic sphincterotomy, biliary sphincterotomy), two rows each, on a logarithmic axis of the stratum odds ratio with a vertical line at OR = 1. Behind each variable is a shaded band marking the range of odds ratios that both stratum confidence intervals contain; all three bands are wide. Strata with p below 0.05 are drawn as filled blue squares and the rest as open brown squares, with the stratum p-value printed at the right of each row. All three variables show the same shape: the upper row, the filled one, has an interval lying entirely left of 1, while the lower row, the open one, has an interval whose right end passes 1 — and the two intervals overlap heavily. Panel B on the right shows the same three variables, one row each, as the ratio of the two stratum odds ratios with its 95% confidence interval, drawn as a red diamond. All three intervals contain 1. Each row is labelled at the right with the ratio, its confidence interval and the P for interaction. — Two readings of one set of numbers. On the left it looks like two different effects; on the right the difference itself has been estimated, and every interval contains 1 and is too wide to rule out either direction.Plotting script `figures/scripts/B2-06-interaction.R`

Take sex, the two rows in the left panel, on their own:

Stratum	n	OR	95% CI	Stratum p
Female	476	0.45	0.26–0.80	0.006
Male	126	0.67	0.23–1.93	0.461

The point estimates, 0.45 and 0.67, do look some distance apart. But:

The female stratum’s confidence interval sits entirely inside the male one. The overlap runs 0.26–0.80, which is 100% of the narrower of the two intervals. Every OR compatible with the data in women is also compatible with the data in men.
The strata differ in size by a factor of 3.78 (476 against 126), and in events by a factor of 3.94 (63 against 16). The male interval is far wider than the female one, and that on its own is enough to carry it across 1.
Estimate the difference itself and the ratio of the two ORs is 1.48 (0.45–4.91), with a P for interaction of 0.52.

The one sentence this page exists for: not significant does not mean not there.

A P for interaction of 0.52 is written up as “no evidence was detected that the treatment effect differs by sex”.

It must not be written as “there was no interaction”, “the effect was the same in men and women”, or “sex does not modify the effect of treatment”. Each of those asserts an equality, and subgroup analyses almost never have the power to support an equality. A widely quoted rule of thumb: detecting an interaction as large as the main effect itself needs roughly four times the sample size that detecting the main effect needs. A trial powered for its primary outcome is therefore, almost by construction, underpowered for its subgroups.

By the same token, the fact that 6 of the variables came out significant in one stratum and not the other must not become “the drug only works in some patients”. The right-hand panel of the figure above is the evidence against that sentence: every interval on the difference contains 1.

Claiming that two strata really do respond alike requires a pre-specified equivalence margin, not a p-value larger than 0.05.

How many tests did you actually run

There is one more reason the stratum p-values on such a plot cannot be taken at face value: there is never just one of them.

This page tested 8 subgroup variables. Suppose the treatment effect varies with none of those characteristics, and each test carries its own 0.05 chance of a false alarm. The probability that at least one p-value comes in under 0.05 is then

1 - (1 - \alpha)^{k}

which works out to 33.7% here — more than one chance in three. And that counts only the interaction tests; go by stratum p-values instead and the number of tests is multiplied by the number of strata, pushing the probability higher still.

Real papers usually report more subgroups than this, so the price is correspondingly larger. For the mechanism, the correction methods, and why pre-specification is the only defence that actually works, see Multiple comparisons and subgroup analyses.

When there are more strata, the stratum estimates break first

This trial carries one more subgroup variable, with 4 levels: the enrolling site. It was deliberately kept out of the forest plot above, and the reason is worth a section:

Site	n	Events	OR	95% CI
1_UM	164	36	0.41	0.19–0.91
2_IU	413	41	0.55	0.28–1.07
3_UK	22	2	1.22	0.07–22.40
4_Case	3	0	not estimable	—

The last site enrolled 3 patients with 0 events, so an entire row of its 2×2 table is zero and the odds ratio simply does not exist. The site above it has 22 patients and 2 events; an OR can be computed, but its interval of 0.07–22.40 has an upper bound 336 times its lower bound, which carries no information either.

The interaction test, meanwhile, still runs, and it is not dragged down by that stratum: it is an overall test on 3 degrees of freedom, with a P for interaction of 0.89.

But one of those degrees of freedom is empty. With 3 patients and 0 events, the interaction parameter belonging to 4_Case is not identifiable in the likelihood — nothing in the data can pin down its value. R reports 3 degrees of freedom because that is the nominal parameter count, not because the data support that many parameters. The effective degrees of freedom are fewer than 3.

It changes nothing here: a P for interaction of 0.89 is far from significant on 3 degrees of freedom and would still be far from it on one fewer. But for a test that landed near the threshold this would decide the conclusion. When you meet an overall test on a many-level subgroup variable, count the events in each level before you decide whether to trust the df.

These two things fail separately, and always in the same order: the stratum estimates go first, the test keeps running.

What that produces is a table that looks complete — every row has numbers in it, the last column has a respectable-looking p-value — while several of those rows are estimating nothing. A reader who does not check the n and the event count cannot tell which.

Three practical rules:

Report n and the event count for every stratum. It is the only thing that lets a reader judge which rows are worth reading.
Do not plot a point estimate for a stratum with too few events, or at least mark it as not estimable, rather than letting it occupy the same visual weight as the rest.
Test many-level variables as a whole rather than comparing levels pairwise — pairwise comparisons multiply the multiplicity problem again. Enrolling site in particular is a poor post-hoc subgroup: differences between sites mix case-mix, operator experience and local process, and are not one interpretable modifier.

LRT or Wald

There are two routes to an interaction test:

The Wald test: read the p-value off the interaction row of summary().
The likelihood ratio test (LRT): compare the likelihoods of the model with the interaction and the model without it, which is anova(fit0, fit, test = "LRT").

For sex the two barely differ: Wald 0.5205 against LRT 0.5222. This page uses the LRT throughout, for three reasons:

A multi-level variable leaves no choice. Enrolling site has 4 levels and its interaction is 3 coefficients; the question is whether those 3 coefficients are simultaneously zero. Wald gives you one p-value per coefficient, and reading them one at a time is another round of multiple testing.
Wald depends on the parameterisation; the LRT does not. Change the reference level or recode the variable and the Wald p-value moves. The LRT does not.
Wald misbehaves in small samples and sparse tables, and it misbehaves by understating significance — the extreme version of this is the separation case on the logistic regression page, where the standard error explodes and the Wald p-value approaches 1. Subgroup analysis is exactly where sparse cells show up.

Common misuses

Misuse	Why it is wrong
Comparing the two strata’s p-values and concluding the effects differ	Those are two tests against a null each; “do they differ” needs the interaction test
Reading P for interaction above 0.05 as “the effect is the same in every group”	Underpowered tests cannot support an equality; write “not detected”
Reading the main effect of an interaction model as the overall effect	It is the effect in the reference stratum, and it moves when the reference moves
Reading the main effect of an uncentred continuous modifier	It is the effect at zero, which may lie outside the data
Splitting the data in two and calling that an interaction analysis	It cannot produce a P for interaction, and covariate effects get estimated twice
Forgetting the covariance when building the other stratum’s CI	The variance of a sum contains 2×Cov; dropping it gets the interval wrong, in a direction set by the sign of Cov
Reporting a subgroup that was chosen after seeing the data	Post-selection p-values are invalid; it is hypothesis-generating at best
Not disclosing how many subgroups were examined	Readers cannot judge the scale of the multiplicity problem
Discussing clinical decisions from a multiplicative interaction alone	Decisions turn on absolute benefit, which needs the risk-difference scale
Plotting a point estimate for a stratum with almost no events	The estimate may not exist, yet it looks as prominent as any other row
Comparing levels of a many-level variable pairwise	Use one overall test with more than one degree of freedom
Slicing by enrolling site after the fact	Site mixes case-mix, operator and process, and is not one interpretable modifier
Using a Wald test for a many-level interaction	Wald tests one coefficient at a time and depends on the parameterisation

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-06-interaction.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

On the subgroup forest plot the odds ratio is significant among women and not among men. Which p-value decides whether the treatment effect differs by sex?

Show the answer and why

Correct answer: 0.522 - the interaction test, which asks directly whether the two strata differ enough

The interaction p-value is 0.522, nowhere near significance: these data did not detect a difference in treatment effect between the sexes. 0.006 and 0.461 are the within-stratum p-values, and "a star on one side and not the other" is not evidence of interaction - it usually reflects unequal stratum sizes, since the larger stratum has more power. Answering a between-stratum question with within-stratum p-values assembles a test that was never performed out of two that address something else.

This trial tested eight subgroup variables, and six of them show one significant stratum and one not. What does that tell you?

Show the answer and why

Correct answer: That the discordant pattern appears for 6 variables, and such a pattern is close to inevitable

Six of eight variables show the pattern, while the number with a significant interaction test is 0 - and the contrast between six and zero is the point: whenever strata differ in size, one star and one blank is close to inevitable, yet it reads as though something has been found. 8 is simply how many variables were tested. What makes subgroup forest plots dangerous is that they render an inevitable visual pattern as though it were a result.

Which number estimates how much the treatment effect actually differs between the sexes?

Show the answer and why

Correct answer: 1.48 - the ratio of the two strata's odds ratios, which estimates the difference itself

The estimate of "how much they differ" is the ratio of the two odds ratios, 1.48, which asks the same question as the interaction test. 0.45 and 0.67 are the strata's own odds ratios: either one alone is that stratum's effect and not a difference between them - and on a ratio scale differences are taken by dividing, not subtracting. The interval around 1.48 is wide, which is the real trouble with subgroup analysis: not whether there is a difference, but that the study was never able to answer the question.

Chapters that use this method

Randomised controlled trial

Sources and licences

This page is original writing