Mixedmodels.Rmd

---
title: "Comparisons of microbiome community composition"
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}  
-->


## Mixed models for univariate comparisons


Load example data:

```{r boxplot-example_w, warning=FALSE, message=FALSE}
# Load libraries
library(microbiome)
library(ggplot2)
library(dplyr)
library(IRanges)

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


Abundance boxplot

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


## Linear model comparison with random effect subject term

Test individual taxonomic group

```{r comparisons-lmer_w, message=FALSE, warning=FALSE}
# Get sample metadata
dfs <- meta(pseq)

# Add abundance as the signal to model
dfs$signal <- abundances(pseq)["Akkermansia", rownames(dfs)]

# Paired comparison
# with fixed group effect and random subject effect
library(lme4)
out <- lmer(signal ~ group + (1|subject), data = dfs)
out0 <- lmer(signal ~ (1|subject), data = dfs)
comp <- anova(out0, out)
pv <- comp[["Pr(>Chisq)"]][[2]]
print(pv)
```