Statistical Inference

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 %>%
  group_by(gender, risk, treatment) %>%
  get_summary_stats(pain_score, type = "common")

## # A tibble: 12 × 13
##    gender risk  treat…¹ varia…²     n   min   max median   iqr  mean    sd    se
##    <fct>  <fct> <fct>   <fct>   <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 male   high  X       pain_s…     6  86.3 100     93.4  6.40  92.7  5.12  2.09
##  2 male   high  Y       pain_s…     6  77.5  91.2   81.2  4.94  82.3  5.00  2.04
##  3 male   high  Z       pain_s…     6  74.4  85.1   80.4  5.38  79.7  4.05  1.65
##  4 male   low   X       pain_s…     6  70.8  81.2   75.9  5.10  76.1  3.86  1.57
##  5 male   low   Y       pain_s…     6  67.9  80.7   73.4  5.40  73.1  4.76  1.94
##  6 male   low   Z       pain_s…     6  68.3  80.4   74.8  7.24  74.5  4.89  2.00
##  7 female high  X       pain_s…     6  68.4  82.7   81.1  2.29  78.9  5.32  2.17
##  8 female high  Y       pain_s…     6  73.1  86.6   81.8  3.66  81.2  4.62  1.89
##  9 female high  Z       pain_s…     6  75.0  87.1   80.8  2.59  81.0  3.98  1.63
## 10 female low   X       pain_s…     6  68.6  78.1   74.6  5.05  74.2  3.69  1.51
## 11 female low   Y       pain_s…     6  63.7  74.7   68.7  4.56  68.4  4.08  1.67
## 12 female low   Z       pain_s…     6  65.4  73.1   69.6  2.78  69.8  2.72  1.11
## # … with 1 more variable: ci <dbl>, and abbreviated variable names ¹treatment,
## #   ²variable

ggplot(headache, aes(x=pain_score,fill = risk)) + 
  geom_histogram(bins=15) + 
  facet_grid(gender ~ treatment) + 
  labs(title = "Dependency of pain score in relation to other factors",
         x = "Pain score",
         y = "Number of observations") +
  theme_light()

ggboxplot(headache, x = "treatment", y = "pain_score", 
  color = "risk", palette = "aas", facet.by = "gender")+
  geom_jitter(aes(color = treatment)) +
  labs(title = "Dependency of pain score in relation to other factors",
         x = "Treatment",
         y = "Pain score") +
  theme_light()

Assumptions

Outliers

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

## # A tibble: 4 × 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

There are 4 outliers in the dataset, moreover one of them is extreme and is equal 68.4.

Normality

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

## # A tibble: 12 × 6
##    gender risk  treatment variable   statistic       p
##    <fct>  <fct> <fct>     <chr>          <dbl>   <dbl>
##  1 male   high  X         pain_score     0.958 0.808  
##  2 female high  X         pain_score     0.714 0.00869
##  3 male   low   X         pain_score     0.982 0.962  
##  4 female low   X         pain_score     0.933 0.600  
##  5 male   high  Y         pain_score     0.902 0.384  
##  6 female high  Y         pain_score     0.939 0.654  
##  7 male   low   Y         pain_score     0.920 0.507  
##  8 female low   Y         pain_score     0.927 0.555  
##  9 male   high  Z         pain_score     0.955 0.784  
## 10 female high  Z         pain_score     0.971 0.901  
## 11 male   low   Z         pain_score     0.924 0.535  
## 12 female low   Z         pain_score     0.958 0.801

ggqqplot(headache, "pain_score", ggtheme = theme_light()) +
  facet_grid(gender + risk ~ treatment)+ 
  labs(title = "Normality of pain score in relation to other factors")

Each of the group passed the normality test, except one - female with high risk treated with X, as their p-value is lower than 0.05 (p=0.00869).On the graph we can clearly see the outlier which greatly influences the normality of that specific group.

Homogeneity of variance

headache %>%
  levene_test(pain_score~gender*treatment*risk)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1    11    60     0.179 0.998

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

To test homogeneity of variance we used Levene Test. The p-value is > 0.05 (0.998), which is not significant. This means that, there is not significant difference between variances across groups. Therefore, we can assume the homogeneity of variances in the different treatment groups. On the plot, there is no significant relationships between residuals and fitted values. That means we can assume the homogeneity of variances.

Anova

results <- aov(data=headache,pain_score ~ gender*risk*treatment)
anova(results)

## Analysis of Variance Table
## 
## Response: pain_score
##                       Df Sum Sq Mean Sq F value            Pr(>F)    
## gender                 1    313     313   16.20           0.00016 ***
## risk                   1   1794    1794   92.70 0.000000000000088 ***
## treatment              2    283     142    7.32           0.00143 ** 
## gender:risk            1      3       3    0.14           0.70849    
## gender:treatment       2    129      65    3.34           0.04220 *  
## risk:treatment         2     28      14    0.71           0.49422    
## gender:risk:treatment  2    287     143    7.41           0.00133 ** 
## Residuals             60   1161      19                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rep<-report(results)

We can see that there is a 3 way correlation because p-value is less than 0.05. Furthermore, there is a significant relationship between pain_score and risk also pain_score and gender.

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

Interaction between pain_score and risk based on gender

model <- lm(pain_score ~ gender * risk, data = headache)
headache %>%
  group_by(gender) %>%
  anova_test(pain_score ~ risk, error = model)

## # A tibble: 2 × 8
##   gender Effect   DFn   DFd     F           p `p<.05`   ges
## * <fct>  <chr>  <dbl> <dbl> <dbl>       <dbl> <chr>   <dbl>
## 1 male   risk       1    68  34.9 0.000000124 *       0.339
## 2 female risk       1    68  29.8 0.000000719 *       0.305

Interaction between pain_score and treatment based on gender

model2 <- lm(pain_score ~ gender * treatment, data = headache)
headache %>%
  group_by(gender) %>%
  anova_test(pain_score ~ treatment, error = model2)

## # A tibble: 2 × 8
##   gender Effect      DFn   DFd     F     p `p<.05`   ges
## * <fct>  <chr>     <dbl> <dbl> <dbl> <dbl> <chr>   <dbl>
## 1 male   treatment     2    66 3.97  0.024 "*"     0.107
## 2 female treatment     2    66 0.188 0.829 ""      0.006

Main effects

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

## # A tibble: 4 × 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    68 10.3   0.000124 "*"     0.232
## 2 female high  treatment     2    68  0.363 0.697    ""      0.011
## 3 male   low   treatment     2    68  0.46  0.633    ""      0.013
## 4 female low   treatment     2    68  1.97  0.147    ""      0.055

It seems that for men in high risk group chosen type of treatment is very influential when it comes to describing pain score.

Pairwise comparisons

headache %>%
  pairwise_t_test(
    pain_score ~ treatment, 
    p.adjust.method = "bonferroni")

## # A tibble: 3 × 9
##   .y.        group1 group2    n1    n2      p p.signif p.adj p.adj.signif
## * <chr>      <chr>  <chr>  <int> <int>  <dbl> <chr>    <dbl> <chr>       
## 1 pain_score X      Y         24    24 0.0514 ns       0.154 ns          
## 2 pain_score X      Z         24    24 0.0505 ns       0.152 ns          
## 3 pain_score Y      Z         24    24 0.994  ns       1     ns

The difference between X-Y and X-Z is statistically significant (p value is close to 0.05, but compared to p value for Y-Z which is close to 1 we should treat them as significant). Difference between Y and Z is not significant.

Summary

3way ANOVA was performed to find the correlation of gender, risk and treatment on pain_score. There was one extreme outlier and there was homogeneity of variances (p > 0.05). There was a statistically significant interaction between gender and risk. So, to sum up, there exist a significant 3 way interaction between gender, risk and treatment.