3-way Anova
3-Way ANOVA
The three-way ANOVA is an extension of the two-way ANOVA for assessing whether there is an interaction effect between three independent categorical variables on a continuous outcome variable.
We’ll use the headache dataset [datarium package], which contains the measures of migraine headache episode pain score in 72 participants treated with three different treatments. The participants include 36 males and 36 females. Males and females were further subdivided into whether they were at low or high risk of migraine.
We want to understand how each independent variable (type of treatments, risk of migraine and gender) interact to predict the pain score.
Descriptive statistics
| gender | risk | treatment | variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| male | high | X | pain_score | 6 | 86.3 | 100.0 | 93.4 | 88.9 | 95.3 | 6.40 | 5.98 | 92.7 | 5.12 | 2.09 | 5.37 |
| male | high | Y | pain_score | 6 | 77.5 | 91.2 | 81.2 | 79.0 | 83.9 | 4.95 | 4.36 | 82.3 | 5.00 | 2.04 | 5.25 |
| male | high | Z | pain_score | 6 | 74.4 | 85.1 | 80.4 | 76.6 | 82.0 | 5.38 | 4.83 | 79.7 | 4.05 | 1.65 | 4.25 |
| male | low | X | pain_score | 6 | 70.8 | 81.2 | 75.9 | 73.6 | 78.7 | 5.10 | 4.60 | 76.1 | 3.85 | 1.57 | 4.04 |
| male | low | Y | pain_score | 6 | 67.9 | 80.7 | 73.4 | 69.5 | 74.9 | 5.40 | 5.30 | 73.1 | 4.76 | 1.95 | 5.00 |
| male | low | Z | pain_score | 6 | 68.3 | 80.4 | 74.8 | 70.7 | 77.9 | 7.24 | 7.06 | 74.5 | 4.89 | 2.00 | 5.13 |
| female | high | X | pain_score | 6 | 68.4 | 82.7 | 81.1 | 79.1 | 81.4 | 2.29 | 1.49 | 78.9 | 5.32 | 2.17 | 5.58 |
| female | high | Y | pain_score | 6 | 73.1 | 86.6 | 81.8 | 80.0 | 83.7 | 3.65 | 3.48 | 81.2 | 4.62 | 1.89 | 4.85 |
| female | high | Z | pain_score | 6 | 75.0 | 87.1 | 80.8 | 79.8 | 82.4 | 2.59 | 2.35 | 81.0 | 3.98 | 1.63 | 4.18 |
| female | low | X | pain_score | 6 | 68.6 | 78.1 | 74.6 | 72.0 | 77.1 | 5.05 | 4.82 | 74.2 | 3.69 | 1.51 | 3.87 |
| female | low | Y | pain_score | 6 | 63.7 | 74.7 | 68.7 | 65.2 | 69.8 | 4.56 | 4.40 | 68.4 | 4.08 | 1.67 | 4.28 |
| female | low | Z | pain_score | 6 | 65.4 | 73.1 | 69.6 | 68.9 | 71.6 | 2.78 | 2.46 | 69.8 | 2.72 | 1.11 | 2.85 |
Assumptions
On a graph we may notice that not every mean is the same, so we can expect to obtain ANOVA result that will indicate rejection of H0- “all means are the same”.
Outliers
Also on a graph we may notice some outliers in FHX, FHY and FHZ, which means we should take care of Max and/or Min values in Females with High risk. I will replace outliers with mean values, then present graph once again.
As seen above, we got rid of outliers in FHX and FHY, however there still seem to be some extreme values in FHZ. I am gonna leave it be, because our dataset is small (only 6 observations per group).
Normality
##
## Shapiro-Wilk normality test
##
## data: headache$pain_score
## W = 1, p-value = 0.08As we can see from both test and QQ plot, our data has normal distribution.
Homogeneity of variance
##
## Bartlett test of homogeneity of variances
##
## data: pain_score by interaction(gender, risk, treatment)
## Bartlett's K-squared = 3, df = 11, p-value = 1Since we have 12 samples and our data is normal we should use Bartlett test for k-samples homogenity of variances. Obtained P value is equal or close to 1, so we do not reject H0-“variances are equal”. I would also like to point out that p-value being equal 1 is very weird, even unnatural.
Anova
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 313.36 | 313.36 | 16.196 | 0.000 |
| risk | 1 | 1793.56 | 1793.56 | 92.699 | 0.000 |
| treatment | 2 | 283.17 | 141.58 | 7.318 | 0.001 |
| gender:risk | 1 | 2.73 | 2.73 | 0.141 | 0.708 |
| gender:treatment | 2 | 129.18 | 64.59 | 3.338 | 0.042 |
| risk:treatment | 2 | 27.59 | 13.80 | 0.713 | 0.494 |
| gender:risk:treatment | 2 | 286.60 | 143.30 | 7.406 | 0.001 |
| Residuals | 60 | 1160.89 | 19.35 | NA | NA |
As we can see, p-value for gender, risk, treatment and interaction between all 3 of them is lower than 0.05 At this point we reject H0-“means are equal”. Just as assumed at the beginning, by looking at a graph.
Post-hoc tests
If there is a significant 3-way interaction effect, you can decompose it into:
- Simple two-way interaction: run two-way interaction at each level of third variable,
- Simple simple main effect: run one-way model at each level of second variable,
- Simple simple pairwise comparisons: run pairwise or other post-hoc comparisons if necessary.
If you do not have a statistically significant three-way interaction, you need to determine whether you have any statistically significant two-way interaction from the ANOVA output. You can follow up a significant two-way interaction by simple main effects analyses and pairwise comparisons between groups if necessary.
Two-way interactions
Step 1.1: third variable - gender - male:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 968 | 968.1 | 45.02 | 0.000 |
| treatment | 2 | 394 | 196.9 | 9.15 | 0.001 |
| risk:treatment | 2 | 203 | 101.6 | 4.72 | 0.016 |
| Residuals | 30 | 645 | 21.5 | NA | NA |
Significant level of p-value for: risk, treatment, interaction risk-treatment.
Step 1.2: third variable - gender - female:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 828.2 | 828.16 | 48.168 | 0.000 |
| treatment | 2 | 18.6 | 9.32 | 0.542 | 0.587 |
| risk:treatment | 2 | 111.0 | 55.48 | 3.227 | 0.054 |
| Residuals | 30 | 515.8 | 17.19 | NA | NA |
Significant level of p-value for: risk.
Step 2.1: third variable - risk - low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 129 | 128.79 | 7.794 | 0.009 |
| treatment | 2 | 119 | 59.48 | 3.600 | 0.040 |
| gender:treatment | 2 | 16 | 8.02 | 0.485 | 0.620 |
| Residuals | 30 | 496 | 16.52 | NA | NA |
Significant level of p-value for: gender, treatment.
Step 2.2: third variable - risk - high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 187 | 187.3 | 8.45 | 0.007 |
| treatment | 2 | 192 | 95.9 | 4.33 | 0.022 |
| gender:treatment | 2 | 400 | 199.9 | 9.01 | 0.001 |
| Residuals | 30 | 665 | 22.2 | NA | NA |
Significant level of p-value for: gender, treatment, interaction gender-treatment.
Step 3.1: third variable - treatment - X:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 373 | 373.0 | 18.0 | 0.000 |
| risk | 1 | 687 | 686.7 | 33.1 | 0.000 |
| gender:risk | 1 | 215 | 215.2 | 10.4 | 0.004 |
| Residuals | 20 | 415 | 20.7 | NA | NA |
Significant level of p-value for: gender, risk and interaction gender-risk.
Step 3.2: third variable - treatment - Y:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 53.0 | 53.0 | 2.472 | 0.132 |
| risk | 1 | 727.1 | 727.1 | 33.929 | 0.000 |
| gender:risk | 1 | 19.6 | 19.6 | 0.913 | 0.351 |
| Residuals | 20 | 428.6 | 21.4 | NA | NA |
Significant level of p-value for: risk.
Step 3.3: third variable - treatment - Z:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 16.6 | 16.6 | 1.04 | 0.320 |
| risk | 1 | 407.4 | 407.4 | 25.64 | 0.000 |
| gender:risk | 1 | 54.6 | 54.6 | 3.43 | 0.079 |
| Residuals | 20 | 317.7 | 15.9 | NA | NA |
Significant level of p-value for: risk.
Main effects
Step 1.1: treatment changes, gender = male + risk = low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| treatment | 2 | 25.5 | 12.8 | 0.623 | 0.55 |
| Residuals | 15 | 307.3 | 20.5 | NA | NA |
No significant value.
Step 1.2: treatment changes, gender = female + risk = low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| treatment | 2 | 109 | 54.7 | 4.36 | 0.032 |
| Residuals | 15 | 188 | 12.6 | NA | NA |
Significant level of p-value.
Step 1.3: treatment changes, gender = male + risk = high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| treatment | 2 | 571 | 285.7 | 12.7 | 0.001 |
| Residuals | 15 | 338 | 22.5 | NA | NA |
Significant level of p-value.
Step 1.4: treatment changes, gender = female + risk = high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| treatment | 2 | 20.1 | 10.1 | 0.461 | 0.639 |
| Residuals | 15 | 327.4 | 21.8 | NA | NA |
No significant value.
Step 2.1: gender changes, treatment = X + risk = low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 10.8 | 10.8 | 0.757 | 0.405 |
| Residuals | 10 | 142.4 | 14.2 | NA | NA |
No significant value.
Step 2.2: gender changes, treatment = Y + risk = low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 68.5 | 68.5 | 3.48 | 0.092 |
| Residuals | 10 | 196.8 | 19.7 | NA | NA |
No significant value.
Step 2.3: gender changes, treatment = Z + risk = low:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 65.6 | 65.6 | 4.19 | 0.068 |
| Residuals | 10 | 156.5 | 15.6 | NA | NA |
No significant value.
Step 2.4: gender changes, treatment = X + risk = high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 577 | 577.4 | 21.2 | 0.001 |
| Residuals | 10 | 272 | 27.2 | NA | NA |
Significant level of p-value.
Step 2.5: gender changes, treatment = Y + risk = high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 4.08 | 4.08 | 0.176 | 0.684 |
| Residuals | 10 | 231.74 | 23.17 | NA | NA |
No significant value.
Step 2.6: gender changes, treatment = Z + risk = high:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| gender | 1 | 5.5 | 5.5 | 0.341 | 0.572 |
| Residuals | 10 | 161.2 | 16.1 | NA | NA |
No significant value.
Step 3.1: risk changes, treatment = X + gender = male:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 835 | 835.4 | 40.7 | 0 |
| Residuals | 10 | 205 | 20.5 | NA | NA |
Significant level of p-value.
Step 3.2: risk changes, treatment = Y + gender = male:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 254 | 254.1 | 10.6 | 0.009 |
| Residuals | 10 | 239 | 23.9 | NA | NA |
Significant level of p-value.
Step 3.3: risk changes, treatment = Z + gender = male:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 81.9 | 81.9 | 4.07 | 0.071 |
| Residuals | 10 | 201.4 | 20.1 | NA | NA |
No significant value.
Step 3.4: risk changes, treatment = X + gender = female:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 66.5 | 66.5 | 3.18 | 0.105 |
| Residuals | 10 | 209.4 | 20.9 | NA | NA |
No significant value.
Step 3.5: risk changes, treatment = Y + gender = female:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 493 | 493 | 25.9 | 0 |
| Residuals | 10 | 190 | 19 | NA | NA |
Significant level of p-value.
Step 3.6: risk changes, treatment = Z + gender = female:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| risk | 1 | 380 | 380.1 | 32.7 | 0 |
| Residuals | 10 | 116 | 11.6 | NA | NA |
Significant level of p-value.
Pairwise comparisons
| gender | risk | .y. | group1 | group2 | statistic | p | p.adj |
|---|---|---|---|---|---|---|---|
| male | high | pain_score | X | Y | 2.894 | 0.034 | 0.102 |
| male | high | pain_score | X | Z | 4.268 | 0.008 | 0.024 |
| male | high | pain_score | Y | Z | 2.861 | 0.035 | 0.106 |
| male | low | pain_score | X | Y | 1.135 | 0.308 | 0.924 |
| male | low | pain_score | X | Z | 0.817 | 0.451 | 1.000 |
| male | low | pain_score | Y | Z | -0.385 | 0.716 | 1.000 |
| female | high | pain_score | X | Y | -0.829 | 0.445 | 1.000 |
| female | high | pain_score | X | Z | -0.803 | 0.458 | 1.000 |
| female | high | pain_score | Y | Z | 0.048 | 0.964 | 1.000 |
| female | low | pain_score | X | Y | 1.921 | 0.113 | 0.339 |
| female | low | pain_score | X | Z | 2.939 | 0.032 | 0.097 |
| female | low | pain_score | Y | Z | -0.618 | 0.563 | 1.000 |