Group-wise comparisons with negative binomial

Read more on negative binomials

Load example data:

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

# Probiotics intervention example data 
data(peerj32) # Source: https://peerj.com/articles/32/
pseq <- peerj32$phyloseq # Rename the example data

Visually compare Akkermansia abundance between time point 1 and 2

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

Test statistical significance with negative binomial:

library(MASS)
library(tidyr)

# Analyse specific taxa
tax <- "Akkermansia"

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

# Negative binomial test with group and gender included
res <- glm.nb(signal ~ group + sex, data = meta(pseq))

# Show the results
print(coef(summary(res)))
##                Estimate Std. Error   z value      Pr(>|z|)
## (Intercept)   6.2466556  0.2244993 27.824837 2.172082e-170
## groupPlacebo  0.3416809  0.2540709  1.344825  1.786817e-01
## sexmale      -0.7876635  0.2624965 -3.000663  2.693926e-03