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
headache %>% sample_n_by(gender, risk, treatment, size = 1)## # A tibble: 12 x 5
## id gender risk treatment pain_score
## <int> <fct> <fct> <fct> <dbl>
## 1 21 male high X 92.7
## 2 30 male high Y 80.1
## 3 33 male high Z 81.3
## 4 2 male low X 76.8
## 5 8 male low Y 80.7
## 6 18 male low Z 74.7
## 7 57 female high X 68.4
## 8 65 female high Y 82.1
## 9 70 female high Z 80.4
## 10 42 female low X 78.1
## 11 48 female low Y 64.1
## 12 49 female low Z 69.0headache %>%
group_by(gender, risk, treatment) %>%
get_summary_stats(pain_score, type = "mean_sd")## # A tibble: 12 x 7
## gender risk treatment variable n mean sd
## <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 male high X pain_score 6 92.7 5.12
## 2 male high Y pain_score 6 82.3 5.00
## 3 male high Z pain_score 6 79.7 4.05
## 4 male low X pain_score 6 76.1 3.86
## 5 male low Y pain_score 6 73.1 4.76
## 6 male low Z pain_score 6 74.5 4.89
## 7 female high X pain_score 6 78.9 5.32
## 8 female high Y pain_score 6 81.2 4.62
## 9 female high Z pain_score 6 81.0 3.98
## 10 female low X pain_score 6 74.2 3.69
## 11 female low Y pain_score 6 68.4 4.08
## 12 female low Z pain_score 6 69.8 2.72effect of the treatment type - focal variable
gender, risk - moderator variables
Effect of the treatment will depend on moderator variables.
bplot <- ggplot(headache, aes(x = treatment, y = pain_score, color = risk)) +
geom_boxplot() +
facet_wrap(~gender)
bplot
X treatment among men has significantly bigger pain score, with the high risk, than the others.
pain score is generally higher when the risk is higher.
Female’s pain score have some outliers, when the risk is high.
Assumptions
Outliers
headache %>% group_by(gender, risk, treatment) %>% identify_outliers(pain_score)## # A tibble: 4 x 7
## gender risk treatment id pain_score is.outlier is.extreme
## <fct> <fct> <fct> <int> <dbl> <lgl> <lgl>
## 1 female high X 57 68.4 TRUE TRUE
## 2 female high Y 62 73.1 TRUE FALSE
## 3 female high Z 67 75.0 TRUE FALSE
## 4 female high Z 71 87.1 TRUE FALSE1 Extreme outlier - id = 57
Normality
Model Residuals
model <- lm(pain_score ~ gender*risk*treatment, data = headache)
ggqqplot(residuals(model))
shapiro_test(residuals(model))## # A tibble: 1 x 3
## variable statistic p.value
## <chr> <dbl> <dbl>
## 1 residuals(model) 0.982 0.398Normality can be assumed since in the QQ plot, all the points are approximately compatible with the reference line.
By groups
headache %>% group_by(gender, risk, treatment) %>% shapiro_test(pain_score)## # A tibble: 12 x 6
## gender risk treatment variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 male high X pain_score 0.958 0.808
## 2 male high Y pain_score 0.902 0.384
## 3 male high Z pain_score 0.955 0.784
## 4 male low X pain_score 0.982 0.962
## 5 male low Y pain_score 0.920 0.507
## 6 male low Z pain_score 0.924 0.535
## 7 female high X pain_score 0.714 0.00869
## 8 female high Y pain_score 0.939 0.654
## 9 female high Z pain_score 0.971 0.901
## 10 female low X pain_score 0.933 0.600
## 11 female low Y pain_score 0.927 0.555
## 12 female low Z pain_score 0.958 0.801ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
facet_grid(gender + risk ~ treatment, labeller = "label_both")
The pain scores are normally distributed - the p-value is bigger than 0.05, for every group except from one - female, drug X, high risk.
In the QQ plot, all the points corresponds to the reference lines except from one that has been defined as an extreme outlier.
Homogeneity of variance
headache %>% levene_test(pain_score ~ gender*risk*treatment)## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 11 60 0.179 0.998The Levene’s test is not significant (p > 0.05). Therefore, we can assume the homogeneity of variances in the different groups.
Anova
res_aov <- headache %>% anova_test(pain_score ~ gender*risk*treatment)## Coefficient covariances computed by hccm()res_aov## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gender 1 60 16.196 0.000163000000000 * 0.213
## 2 risk 1 60 92.699 0.000000000000088 * 0.607
## 3 treatment 2 60 7.318 0.001000000000000 * 0.196
## 4 gender:risk 1 60 0.141 0.708000000000000 0.002
## 5 gender:treatment 2 60 3.338 0.042000000000000 * 0.100
## 6 risk:treatment 2 60 0.713 0.494000000000000 0.023
## 7 gender:risk:treatment 2 60 7.406 0.001000000000000 * 0.198There is a statistically significant three-way interaction between gender, risk and treatment.
p = 0.001
F(2, 60) = 7.41
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
model <- lm(pain_score ~ gender*risk*treatment, data = headache)
headache %>% group_by(gender) %>% anova_test(pain_score ~ risk*treatment, error = model)## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()## # A tibble: 6 x 8
## gender Effect DFn DFd F p `p<.05` ges
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male risk 1 60 50.0 0.00000000187 "*" 0.455
## 2 male treatment 2 60 10.2 0.000157 "*" 0.253
## 3 male risk:treatment 2 60 5.25 0.008 "*" 0.149
## 4 female risk 1 60 42.8 0.000000015 "*" 0.416
## 5 female treatment 2 60 0.482 0.62 "" 0.016
## 6 female risk:treatment 2 60 2.87 0.065 "" 0.087For males - there is a statistically significant simple two-way interaction between risk and treatment
p = 0.008
F(2, 60) = 5.25
For Females - there is NOT a statistically significant simple two-way interaction between risk and treatment p = 0.065
F(2, 60) = 2.87
Main effects
effects <- headache %>% group_by(gender, risk) %>% anova_test(pain_score ~ treatment, error = model)## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()effects[1:2,]## # A tibble: 2 x 9
## gender risk Effect DFn DFd F p `p<.05` ges
## <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male high treatment 2 60 14.8 0.0000061 "*" 0.33
## 2 male low treatment 2 60 0.66 0.521 "" 0.022For males with high risk - There is a statistically significant simple simple main effect of treatment
p < 0.0001
F(2, 60) = 14.8
For males with low risk - There is NOT a statistically significant simple simple main effect of treatment
p = 0.521
F(2, 60) = 0.66
Pairwise comparisons
pwc <- headache %>% group_by(gender, risk) %>% emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>% select(-df, -statistic, -p)
pwc %>% filter(gender == "male", risk == "high")## # A tibble: 3 x 8
## gender risk term .y. group1 group2 p.adj p.adj.signif
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 male high treatment pain_score X Y 0.000386 ***
## 2 male high treatment pain_score X Z 0.00000942 ****
## 3 male high treatment pain_score Y Z 0.897 nsget_emmeans(pwc) %>% filter(gender == "male", risk == "high")## # A tibble: 3 x 9
## gender risk treatment emmean se df conf.low conf.high method
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 male high X 92.7 1.80 60 89.1 96.3 Emmeans test
## 2 male high Y 82.3 1.80 60 78.7 85.9 Emmeans test
## 3 male high Z 79.7 1.80 60 76.1 83.3 Emmeans testFor males with high risk:
Treatment X and Y - there is a statistically significant mean difference of 10.4
Treatment X and Z - there is a statistically significant mean difference of 13.1
Treatment Y and Z - there is NOT a statistically significant mean difference (p.adj = 0.897)
Final Boxplots & p-values
pwc <- pwc %>% add_xy_position(x = "treatment")
pwc.filtered <- pwc %>% filter(gender == "male", risk == "high")
bplot +
stat_pvalue_manual(
pwc.filtered, color = "risk", linetype = "risk", hide.ns = TRUE,
tip.length = 0, step.increase = 0.1, step.group.by = "gender"
) +
labs(
subtitle = get_test_label(res_aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)