3-way Anova

Outliers

tab <- headache %>%
  group_by(gender, risk, treatment) %>%
  identify_outliers(pain_score)

tab$is.extreme <- cell_spec(tab$is.extreme, color =  ifelse(tab$is.extreme == TRUE, "red", "white"))

tab %>%
  kbl(caption = "Outliers identification", escape = F) %>% 
  kable_material_dark("striped")

It can be seen that, the data contain one extreme outlier (id = 57, female at high risk of migraine taking drug X) In further analysis I have decided to keep the outlier and see the results of the analysis. Then, analysis could be repeated without the outlier.

Normality

Check normality assumption by analyzing the model residuals. QQ plot and Shapiro-Wilk test of normality are used.

model  <- lm(pain_score ~ gender*risk*treatment, data = headache)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

In the QQ plot, as all the points fall approximately along the reference line, we can assume normality.

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 x 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.982   0.398

This conclusion is supported by the Shapiro-Wilk test. The p-value is not significant (p = 0.4), so we can assume normality.

Check normality assumption by groups. Computing Shapiro-Wilk test for each combinations of factor levels.

sw <- headache %>%
  group_by(gender, risk, treatment) %>%
  shapiro_test(pain_score)

sw$p <- cell_spec(round(sw$p, 4), color =  ifelse(sw$p < 0.05, "red", "white"))

sw %>% kbl(caption = "Outliers identification", escape = F) %>% 
  kable_material_dark("striped")

The pain scores were normally distributed (p > 0.05) except for one group (female at high risk of migraine taking drug X, p = 0.0086), as assessed by Shapiro-Wilk’s test of normality.

ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
  facet_grid(gender + risk ~ treatment, labeller = "label_both")

All the points fall approximately along the reference line, except for one group (female at high risk of migraine taking drug X), where we already identified an extreme outlier. #### Homogeneity of variance This can be checked using the Levene’s test:

headache %>% levene_test(pain_score ~ gender*risk*treatment) %>% 
  kbl(caption = "Levene's test") %>% 
  kable_material_dark()

The Levene’s test is not significant (p > 0.05). Therefore, we can assume the homogeneity of variances in the different groups.

plot(lm(pain_score~gender*risk*treatment, data=headache), 1)

Points 29, 57 and 62 are leverage points which violate the homogenity of variances assumption. However the homogenity was proved in Levene’s test.

Two-way interactions

We want to evaluate the effect of risk*treatment interaction on pain_score at each level of gender.

model  <- lm(pain_score ~ gender*risk*treatment, data = headache)
headache %>%
  group_by(gender) %>%
  anova_test(pain_score ~ risk*treatment, error = model) %>% 
  kbl(caption = "Two-way interactions") %>% 
  kable_material_dark()

For males, this result suggests that the effect of treatment on “pain_score” depends on one’s “risk” of migraine. In other words, the risk moderates the effect of the type of treatment on pain_score.

Main effects

Group the data by gender and risk and analyze the simple simple main effects of treatment on pain_score:

treatment.effect <- headache %>%
  group_by(gender, risk) %>%
  anova_test(pain_score ~ treatment, error = model)

treatment.effect %>% kbl(caption = "Treatment effect") %>% 
  kable_material_dark() %>%
  row_spec(c(1,2), color ="red", bold = T)

This analysis indicates that, the type of treatment taken has a statistically significant effect on pain_score in males who are at high risk.

In other words, the mean pain_score in the treatment X, Y and Z groups was statistically significantly different for males who at high risk, but not for males at low risk.

Pairwise comparisons

Compare the different treatments by gender and risk variables:

pwc <- headache %>%
  group_by(gender, risk) %>%
  emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>%
  select(-df, -statistic, -p) # Remove details
# Show comparison results for male at high risk
pwc %>% filter(gender == "male", risk == "high") %>%
  kbl(caption = "Pairwise comparisons") %>% 
  kable_material_dark()

There is a significant difference between treatments X and Y as well as between X and Z, however the difference between treatments Z and Y is statistically insignificant.