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

ggplot(headache, aes(x=pain_score,fill = risk)) + 
  geom_histogram(bins=10) + 
  facet_grid(gender ~ treatment) + 
  theme_classic()

bxp <- ggboxplot(
  headache, x = "treatment", y = "pain_score", 
  color = "risk", palette = "jco", facet.by = "gender"
  )+
  geom_jitter(aes(color = treatment))
bxp

index <- c(1:72)
ggplot(headache, aes(index,y=pain_score,color = treatment,size = risk, alpha = 0.2)) + 
  geom_point() + 
  scale_color_manual(values = c('#FF1493','#0000FF','#00ff00')) +
    theme(plot.margin=unit(c(1,2,0,0),"cm")) +
  theme(legend.position=c(1,1), legend.justification=c(0,1))+
    facet_wrap(~gender)

## Warning: Using size for a discrete variable is not advised.

on on this grpah we can see scores that will be shown in outlier

headache %>%
  group_by(gender, risk, treatment) %>%
  get_summary_stats(pain_score, type = "full",
  show = c("n","mean","sd","max","min","se","q1","q3","median"))

## # A tibble: 12 x 13
##    gender risk  treatment variable       n  mean    sd   max   min    se    q1
##    <fct>  <fct> <fct>     <chr>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 male   high  X         pain_score     6  92.7  5.12 100    86.3  2.09  88.9
##  2 male   high  Y         pain_score     6  82.3  5.00  91.2  77.5  2.04  79.0
##  3 male   high  Z         pain_score     6  79.7  4.05  85.1  74.4  1.65  76.6
##  4 male   low   X         pain_score     6  76.1  3.86  81.2  70.8  1.57  73.6
##  5 male   low   Y         pain_score     6  73.1  4.76  80.7  67.9  1.94  69.5
##  6 male   low   Z         pain_score     6  74.5  4.89  80.4  68.3  2.00  70.7
##  7 female high  X         pain_score     6  78.9  5.32  82.7  68.4  2.17  79.1
##  8 female high  Y         pain_score     6  81.2  4.62  86.6  73.1  1.89  80.0
##  9 female high  Z         pain_score     6  81.0  3.98  87.1  75.0  1.63  79.8
## 10 female low   X         pain_score     6  74.2  3.69  78.1  68.6  1.51  72.0
## 11 female low   Y         pain_score     6  68.4  4.08  74.7  63.7  1.67  65.2
## 12 female low   Z         pain_score     6  69.8  2.72  73.1  65.4  1.11  68.9
## # ... with 2 more variables: q3 <dbl>, median <dbl>

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       FALSE

Normality

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

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.801

there is p-value under 0.05 almost everywhere exept where we found outlier in group gender:female, risk:high, treatment:X where p-value is 0.00869, and on the qqplot we can see that not all points are falling in normality spectrum in this group

model  <- lm(pain_score ~ gender*risk*treatment, data = headache)
# Create a QQ plot of residuals
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.398

p-value is not significant we can assume normality in relationship between pain_score and interactions between gender,risk, treatment

Homogeneity of variance

bartlett.test(pain_score ~ interaction(gender,treatment,risk), data=headache)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pain_score by interaction(gender, treatment, risk)
## Bartlett's K-squared = 3, df = 11, p-value = 1

since samples are normal I’m using Bartlett test to test if variances are equal, p-value is not significant so they are equal

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.198

there is significant p-value in interaction gender:ris and risk:treatment, so we have to do post-hoc test

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

# Group the data by gender and 
# fit simple two-way interaction 
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.087

the only statsitcally signifact interaction is risk:treatment(p<0.05,p=0.0014) when it comes to males, now we have to look for which level of risk is the most sensitive for treatment, we can rejcet the rest of effects becasue we proved that there is no other significance explenation of pain_score

Main effects

treatment.effect <- headache %>%
  group_by(risk, gender) %>%
  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()

treatment.effect

## # A tibble: 4 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 female high  treatment     2    60  0.52 0.597     ""      0.017
## 3 male   low   treatment     2    60  0.66 0.521     ""      0.022
## 4 female low   treatment     2    60  2.83 0.067     ""      0.086

main effect that is significant(p-value<0.05, p = 0.000006) is effect of treatment for males at high risk of migraine, that suggest that type of treatment is signifficant only with males with high risk of migraine #### 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      ns

there is significance difference between treatment in pair X and Y(p-value<0.05) and in pair X and Z (p-value<0.05), there is no significance difference between Y an Z treatments so those can be used exchangable(only if you are male with high risk of migraine) we can confirm that running tukey honest significant difference test

result <- aov(formula = pain_score ~ gender*risk*treatment, data = headache)
TukeyHSD(result,"treatment")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pain_score ~ gender * risk * treatment, data = headache)
## 
## $treatment
##        diff   lwr   upr p adj
## Y-X -4.1986 -7.25 -1.15 0.004
## Z-X -4.2152 -7.27 -1.16 0.004
## Z-Y -0.0166 -3.07  3.04 1.000

###Now we can perform two-way interactions and pairwise comparisons for another significant interaction

Two-way interactions

# Group the data by gender and 
# fit simple two-way interaction 
headache %>%
  group_by(treatment) %>%
  anova_test(pain_score ~ risk*gender, error = model)

## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

## # A tibble: 9 x 8
##   treatment Effect        DFn   DFd      F            p `p<.05`   ges
## * <fct>     <chr>       <dbl> <dbl>  <dbl>        <dbl> <chr>   <dbl>
## 1 X         risk            1    60 35.5   0.000000145  "*"     0.372
## 2 X         gender          1    60 19.3   0.0000466    "*"     0.243
## 3 X         risk:gender     1    60 11.1   0.001        "*"     0.156
## 4 Y         risk            1    60 37.6   0.0000000744 "*"     0.385
## 5 Y         gender          1    60  2.74  0.103        ""      0.044
## 6 Y         risk:gender     1    60  1.01  0.319        ""      0.017
## 7 Z         risk            1    60 21.1   0.0000233    "*"     0.26 
## 8 Z         gender          1    60  0.856 0.359        ""      0.014
## 9 Z         risk:gender     1    60  2.82  0.098        ""      0.045

we can see there is only one siginificance level for this interaction(p-value =0.001 ) and it appers with ‘X’ treatment

Main effects

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

## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

gender.effect

## # A tibble: 6 x 9
##   risk  treatment Effect   DFn   DFd      F        p `p<.05`   ges
## * <fct> <fct>     <chr>  <dbl> <dbl>  <dbl>    <dbl> <chr>   <dbl>
## 1 high  X         gender     1    10 21.2   0.000971 "*"     0.68 
## 2 low   X         gender     1    10  0.757 0.405    ""      0.07 
## 3 high  Y         gender     1    10  0.176 0.684    ""      0.017
## 4 low   Y         gender     1    10  3.48  0.092    ""      0.258
## 5 high  Z         gender     1    10  0.341 0.572    ""      0.033
## 6 low   Z         gender     1    10  4.19  0.068    ""      0.295

this effect show us that we werent wrong in previous time because we can see that there is significant difference only when the risk is high pairwise comparison is unnessecary here because we have only two groups(male,female) so we know that the difference is significant in those groups

Summary

pwc <- pwc %>% add_xy_position(x="treatment")
pwc.filtered <-filter(pwc,gender == "male", risk == "high")
bxp + 
  stat_pvalue_manual(   pwc.filtered) +
  labs(
    subtitle=get_test_label(res.aov,detailed=TRUE),
    caption=get_pwc_label(pwc)
  )

A three-way ANOVA was conducted to examine the effects of gender, risk and treatment on migraine headache episode pain_score.

Residual analysis was performed to test for the assumptions of the two-way ANOVA. Outliers were assessed by box plot method, normality was assessed using Shapiro-Wilk’s normality test and homogeneity of variances was assessed by Bartlett test

There was one extreme outlier, residuals were normally distributed (p > 0.05) and there was homogeneity of variances (p > 0.05).

There was a statistically significant interaction between gender and risk of migraine on pain score, F(2, 60) = 0.141, p = 0.002 and between risk and treatment, F(2, 60) = 0.713, p = 0.023

Consequently, an analysis of simple main effects for treatment and gender was performed with statistical significance receiving a Bonferroni adjustment. There was a statistically significant difference in mean of pain score(when it comes to effect for treatment) for males with high risk of migraine (F = 14.77, p = 0.0000061) And with effect for gender there was significant difference in mean of pain socre for treatment X with high risk of migriane (F = 21.214, p = 0.0000061)

All pairwise comparisons were analyzed between the different treatment groups organized by risk and gender There was a significant difference of pain score between groups X Y and X Z (p < 0.05)

Taking into consideration whole raport and summing up conclusions for males in high risk of migraine the best treatment is X(looking on graph and taking into consideration second interaction)