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

Basic statistics for groups by three factors

## # A tibble: 12 × 7
##    gender risk  treatment variable       n  mean    sd
##    <fct>  <fct> <fct>     <fct>      <dbl> <dbl> <dbl>
##  1 male   high  X         pain_score     6  92.7  5.12
##  2 male   low   X         pain_score     6  76.1  3.86
##  3 male   high  Y         pain_score     6  82.3  5.00
##  4 male   low   Y         pain_score     6  73.1  4.76
##  5 male   high  Z         pain_score     6  79.7  4.05
##  6 male   low   Z         pain_score     6  74.5  4.89
##  7 female high  X         pain_score     6  78.9  5.32
##  8 female low   X         pain_score     6  74.2  3.69
##  9 female high  Y         pain_score     6  81.2  4.62
## 10 female low   Y         pain_score     6  68.4  4.08
## 11 female high  Z         pain_score     6  81.0  3.98
## 12 female low   Z         pain_score     6  69.8  2.72

Assumptions

Outliers

Various outliers have been identified.

## # 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

The extreme one is for female no. 57 with high risk, treated with X.

Normality

Shapiro-Wilk test and a Quantile-Quantile plot.

## # A tibble: 1 × 3
##   variable   statistic      p
##   <chr>          <dbl>  <dbl>
## 1 pain_score     0.970 0.0847

p-val = 0.085 > 0.05, generally normal

Testing normality by all 3 factors with Shapiro-Wilk test:

## # 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 male   low   X         pain_score     0.982 0.962  
##  3 male   high  Y         pain_score     0.902 0.384  
##  4 male   low   Y         pain_score     0.920 0.507  
##  5 male   high  Z         pain_score     0.955 0.784  
##  6 male   low   Z         pain_score     0.924 0.535  
##  7 female high  X         pain_score     0.714 0.00869
##  8 female low   X         pain_score     0.933 0.600  
##  9 female high  Y         pain_score     0.939 0.654  
## 10 female low   Y         pain_score     0.927 0.555  
## 11 female high  Z         pain_score     0.971 0.901  
## 12 female low   Z         pain_score     0.958 0.801

Visualising normality violations by gender, risk and treatment: Above QQ plot and a grouped test confirm normality for all groups apart from females with high risk treated with X (p = 0.008) where the extreme outlier has been found.

Homogeneity of variance

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

P-value, according to the Levene’s test above, is almost one and thus not significant. The homogeneity of variance in three groups is preserved.

Also the plot of residuals versus fitted values proves the homogeneity of variance.

Anova

## ANOVA Table (type II tests)
## 
##                  Effect DFn DFd      F                 p p<.05   ges
## 1                  risk   1  60 92.699 0.000000000000088     * 0.607
## 2                gender   1  60 16.196 0.000163000000000     * 0.213
## 3             treatment   2  60  7.318 0.001000000000000     * 0.196
## 4           risk:gender   1  60  0.141 0.708000000000000       0.002
## 5        risk:treatment   2  60  0.713 0.494000000000000       0.023
## 6      gender:treatment   2  60  3.338 0.042000000000000     * 0.100
## 7 risk:gender:treatment   2  60  7.406 0.001000000000000     * 0.198

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

Anova performed with anova_test (between) and aov indicate significant (p = 0.1%, almost zero at F = 7.41) interaction between gender, treatment and risk.

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

The analysis for each case of gender.

## # A tibble: 6 × 8
##   gender Effect           DFn   DFd      F             p `p<.05`   ges
## * <fct>  <chr>          <dbl> <dbl>  <dbl>         <dbl> <chr>   <dbl>
## 1 male   treatment          2    60 10.2   0.000157      "*"     0.253
## 2 male   risk               1    60 50.0   0.00000000187 "*"     0.455
## 3 male   treatment:risk     2    60  5.25  0.008         "*"     0.149
## 4 female treatment          2    60  0.482 0.62          ""      0.016
## 5 female risk               1    60 42.8   0.000000015   "*"     0.416
## 6 female treatment:risk     2    60  2.87  0.065         ""      0.087

A significant (p = 0.8% at F(2, 60) = 5.25) interaction has been observed between treatment and risk for males.

Main effects

The analysis for each case of gender and risk.

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

A significant (p = 0 at F(2, 60) = 14.77) main effect of treatment for males with high risk has been unveiled. Hence their level of pain is prone to changes by a treatment.

Pairwise comparisons

Estimated Marginal Means with Bonferroni correction:

## # A tibble: 12 × 11
##    gender risk  term   .y.   group1 group2    df stati…¹       p   p.adj p.adj…²
##  * <chr>  <chr> <chr>  <chr> <chr>  <chr>  <dbl>   <dbl>   <dbl>   <dbl> <chr>  
##  1 female high  treat… pain… X      Y         60 -0.910  3.67e-1 1   e+0 ns     
##  2 female high  treat… pain… X      Z         60 -0.855  3.96e-1 1   e+0 ns     
##  3 female high  treat… pain… Y      Z         60  0.0552 9.56e-1 1   e+0 ns     
##  4 female low   treat… pain… X      Y         60  2.28   2.61e-2 7.82e-2 ns     
##  5 female low   treat… pain… X      Z         60  1.72   9.00e-2 2.70e-1 ns     
##  6 female low   treat… pain… Y      Z         60 -0.558  5.79e-1 1   e+0 ns     
##  7 male   high  treat… pain… X      Y         60  4.09   1.29e-4 3.86e-4 ***    
##  8 male   high  treat… pain… X      Z         60  5.14   3.14e-6 9.42e-6 ****   
##  9 male   high  treat… pain… Y      Z         60  1.05   2.99e-1 8.97e-1 ns     
## 10 male   low   treat… pain… X      Y         60  1.15   2.56e-1 7.68e-1 ns     
## 11 male   low   treat… pain… X      Z         60  0.628  5.32e-1 1   e+0 ns     
## 12 male   low   treat… pain… Y      Z         60 -0.519  6.06e-1 1   e+0 ns     
## # … with abbreviated variable names ¹statistic, ²p.adj.signif

Pairwise by treatment:

## # 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

Pairwise by risk:

## # A tibble: 1 × 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 high   low       36    36 1.24e-10 ****     1.24e-10 ****

Means between risk groups are significantly different (at p = 0).

Pairwise by gender:

## # A tibble: 1 × 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 male   female    36    36 0.0172 *        0.0172 *

Results

The normality and homoscedasticity of data has been verified positively (except for one outlier of a female no. 57). This analysis has shown that all three examined factors interact at a significance level of p = 0.001 (F(2, 60) = 7.41). It has also unveiled an interaction between treatment and risk for males (p = 0.8% at F(2, 60) = 5.25). Moreover, a significant main effect has been identified: treatment for males with high risk (p = 0 at F(2, 60) = 14.77).