Negativebinomial.Rmd

---
title: "Negative binomial"
author: "Leo Lahti, Sudarshan Shetty et al."
bibliography: 
- bibliography.bib
output:
  BiocStyle::html_document:
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float: true
    self_contained: true
    thumbnails: true
    lightbox: true
    gallery: true
    use_bookdown: false
    highlight: haddock
---
<!--
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{microbiome tutorial - comparisons}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->


## Group-wise comparisons with negative binomial

[Read more on negative binomials](http://www.ats.ucla.edu/stat/r/dae/nbreg.htm)

Load example data:

```{r boxplot-example, warning=FALSE, message=FALSE}
# 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

```{r boxplot2, warning=FALSE, message=FALSE}
p <- boxplot_abundance(pseq, x = "time", y = "Akkermansia", line = "subject") + scale_y_log10()
print(p)
```

Test statistical significance with negative binomial:

```{r comparisons2b, message=FALSE, error=FALSE, warning=FALSE}
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)))
```