Multiple pairwise-comparison between groups

Load example data:

# Load libraries
library(microbiome)
library(ggplot2)
library(dplyr)

data("dietswap")
pseq <- dietswap # Rename the example data

Visually compare Akkermansia abundance between “DI,” “ED” and “HE”

p <- boxplot_abundance(pseq, x = "group", y = "Akkermansia") + scale_y_log10()
print(p)

Test statistical significance by Kruskal-Wallis

library(tidyr)

# Analyse specific taxa
tax <- "Akkermansia"

# Pick the signal (abundance) for this tax
sample_data(pseq)$signal <- get_sample(pseq, tax)


attach(meta(pseq))
res<- pairwise.wilcox.test(signal, group, p.adjust.method = 'BH')
resK<- kruskal.test(signal, group)
print(resK)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  signal and group
## Kruskal-Wallis chi-squared = 2.3967, df = 2, p-value = 0.3017

From the output of the Kruskal-Wallis test, if there is a significant difference between groups, but we don’t know which pairs of groups are different, then we shall move to:

Test statistical significance with common post hoc as follows:

1. Pairwise Wilcoxon tests with multiple testing correction, p.adjust.method is Benjamini-Hochberg
library(tidyr)

# Analyse specific taxa
tax <- "Akkermansia"

# Pick the signal (abundance) for this tax
sample_data(pseq)$signal <- get_sample(pseq, tax)


attach(meta(pseq))
res<- pairwise.wilcox.test(signal, group, p.adjust.method = 'BH')

# Show the results
print(res)
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  signal and group 
## 
##    DI   ED  
## ED 0.41 -   
## HE 0.45 0.58
## 
## P value adjustment method: BH

2. Dunn’s test, p.adjust.method is Benjamini-Hochberg

library(dunn.test)

# Analyse specific taxa
tax <- "Akkermansia"

# Pick the signal (abundance) for this tax
sample_data(pseq)$signal <- get_sample(pseq, tax)


attach(meta(pseq))
res<- dunn.test::dunn.test(signal, group, method = 'bh')  
##   Kruskal-Wallis rank sum test
## 
## data: signal and group
## Kruskal-Wallis chi-squared = 2.3967, df = 2, p-value = 0.3
## 
## 
##                          Comparison of signal by group                         
##                              (Benjamini-Hochberg)                              
## Col Mean-|
## Row Mean |         DI         ED
## ---------+----------------------
##       ED |   1.524102
##          |     0.1912
##          |
##       HE |   1.011341  -0.518075
##          |     0.2339     0.3022
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

3. Pairwise Multiple Comparison of Mean Ranks (PMCMR)

library(PMCMRplus)

# Analyse specific taxa
tax <- "Akkermansia"

# Pick the signal (abundance) for this tax
sample_data(pseq)$signal <- get_sample(pseq, tax)


attach(meta(pseq))
res<- PMCMRplus::kwAllPairsNemenyiTest(signal, group,dist='Chisq',p.adjust.methods=  'BH')  

# Show the results
print(res)
##    DI   ED  
## ED 0.31 -   
## HE 0.60 0.87

4. Conover-Iman Test of Multiple Comparisons Using Rank Sums, p.adjust.method is Benjamini-Hochberg

library(conover.test)

# Analyse specific taxa
tax <- "Akkermansia"

# Pick the signal (abundance) for this tax
sample_data(pseq)$signal <- get_sample(pseq, tax)


attach(meta(pseq))
res<- conover.test::conover.test(signal, group, method = 'bh')
##   Kruskal-Wallis rank sum test
## 
## data: signal and group
## Kruskal-Wallis chi-squared = 2.3967, df = 2, p-value = 0.3
## 
## 
##                          Comparison of signal by group                         
##                              (Benjamini-Hochberg)                              
## Col Mean-|
## Row Mean |         DI         ED
## ---------+----------------------
##       ED |   1.525485
##          |     0.1929
##          |
##       HE |   1.012258  -0.518545
##          |     0.2344     0.3023
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2