How to do a two-way ANOVA test with interaction (in Python, using Statsmodels)
Task
When we analyze the impact that two factors have on a response variable, we often consider the possible relationship between the two factors. That is, does their interaction term affect the response variable? A two-way ANOVA test with interaction can answer that question.
Related tasks:
- How to do a one-way analysis of variance (ANOVA)
- How to do a two-way ANOVA test without interaction
- How to compare two nested linear models using ANOVA
- How to conduct a mixed designs ANOVA
- How to conduct a repeated measures ANOVA
- How to perform an analysis of covariance (ANCOVA)
Solution
We’re going to use R’s esoph
dataset, about esophageal cancer cases.
We will focus on the impact of age group (agegp
) and alcohol consumption (alcgp
)
on the number of cases of the cancer (ncases
). We ask, does the interaction
between these two factors affect the number of cases?
First, we load in the dataset. (See how to quickly load some sample data.)
1
2
3
from rdatasets import data
data = data('esoph')
data.head()
agegp | alcgp | tobgp | ncases | ncontrols | |
---|---|---|---|---|---|
0 | 25-34 | 0-39g/day | 0-9g/day | 0 | 40 |
1 | 25-34 | 0-39g/day | 10-19 | 0 | 10 |
2 | 25-34 | 0-39g/day | 20-29 | 0 | 6 |
3 | 25-34 | 0-39g/day | 30+ | 0 | 5 |
4 | 25-34 | 40-79 | 0-9g/day | 0 | 27 |
Next, we create a model that includes the response variable we care about, plus the two categorical variables we will be testing, as well as their interaction term.
1
2
3
4
import statsmodels.api as sm
from statsmodels.formula.api import ols
# C(...) means the variable is categorical, and : means multiplication.
model = ols('ncases ~ C(alcgp) + C(agegp) + C(alcgp):C(agegp)', data = data).fit()
A two-way ANOVA with interaction tests the following three null hypotheses.
- There is no interaction between the two categorical variables. (If we reject this we do not test the other two hypotheses.)
- The mean response is the same across all groups of the first factor.
(In our example, that says the mean
ncases
is the same for all age groups.) - The mean response is the same across all groups of the second factor.
(In our example, that says the mean
ncases
is the same for all alcohol consumption groups.)
We choose a value, $0 \le \alpha \le 1$, as the Type I Error Rate. Let’s let $\alpha=0.05$ here.
1
sm.stats.anova_lm(model, typ=2)
sum_sq | df | F | PR(>F) | |
---|---|---|---|---|
C(alcgp) | 52.695287 | 3.0 | 4.723387 | 4.862447e-03 |
C(agegp) | 267.026108 | 5.0 | 14.361068 | 2.021935e-09 |
C(alcgp):C(agegp) | 107.557743 | 15.0 | 1.928206 | 3.632710e-02 |
Residual | 238.000000 | 64.0 | NaN | NaN |
The $p$-value for the interaction of age group and alcohol consumption is in the third row, final column, $3.63271\times10^{-2}$. It is less than $\alpha$, so we can reject the null hypothesis that age group and alcohol consumption do not interact with regards to the number of esophageal cancer cases. That is, we have reason to believe that their interaction does effect the outcome.
As we stated when we listed the hypotheses to test, since we rejected the first null hypothesis, we will not test the other two. In the case where you failed to reject the first null hypothesis, you could consider each $p$-value in the first two rows of the above table, one for each of the two factors.
Content last modified on 24 July 2023.
See a problem? Tell us or edit the source.
Contributed by Elizabeth Czarniak (CZARNIA_ELIZ@bentley.edu)