---
title: "MaAsLin 3 Tutorial"
author:
- name: William Nickols
email: willnickols@g.harvard.edu
- name: Jacob Nearing
email: nearing@broadinstitute.org
- name: Sagun Maharjan
email: smaharjan@hsph.harvard.edu
output: html_document
vignette: >
%\VignetteIndexEntry{Tutorial}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## Introduction ##
[**MaAsLin 3**](http://huttenhower.sph.harvard.edu/MaAsLin3) is the next
generation of **M**a**A**s**L**in (**M**icrobiome **M**ultivariable
**A**ssociations with **L**inear Models). This comprehensive R package
efficiently determines multivariable associations between clinical
metadata and microbial meta-omics features. Relative to MaAsLin 2,
MaAsLin 3 introduces the ability to quantify and test for **both
abundance and prevalence associations** while **accounting for
compositionality**. By incorporating generalized linear models, MaAsLin
3 accommodates most modern epidemiological study designs including
cross-sectional and longitudinal studies.
If you use the MaAsLin 3 software, please cite our manuscript:
> William A. Nickols, Thomas Kuntz, Jiaxian Shen, Sagun Maharjan,
Himel Mallick, Eric A. Franzosa, Kelsey N. Thompson, Jacob T. Nearing,
Curtis Huttenhower. MaAsLin 3: Refining and extending generalized
multivariable linear models for meta-omic association discovery.
bioRxiv 2024.12.13.628459; doi: https://doi.org/10.1101/2024.12.13.628459
### Support ###
The user manual can be found
[here](https://github.com/biobakery/maaslin3). If using vignettes,
users should start with the `maaslin3_tutorial.Rmd` vignette and then
refer to the `maaslin3_manual.Rmd` vignette as necessary. If you have questions,
please direct them to the [MaAsLin
Forum](https://forum.biobakery.org/c/Downstream-analysis-and-statistics/maaslin)
## Contents
* [1. Installing R](#installing-r)
* [2. Installing MaAsLin 3](#installing-maaslin-3)
* [3. Microbiome association detection with MaAsLin
3](#microbiome-association-detection-with-maaslin-3)
* [3.1 MaAsLin 3 input](#maaslin-3-input)
* [3.2 Running MaAsLin 3](#running-maaslin-3)
* [3.3 MaAsLin 3 output](#maaslin-3-output)
* [4. Advanced Topics](#advanced-topics)
* [4.1 Absolute abundance](#absolute-abundance)
* [4.2 Random Effects (repeated sampling)](#random-effects)
* [4.3 Interactions (differences in differences)](#interactions)
* [4.4 Level contrasts](#level-contrasts)
* [4.5 Group-wise differences](#group-wise-differences)
* [4.6 Contrast tests](#contrast-tests)
* [5. Command line](#command-line)
## 1. Installing R {#installing-r}
[R](https://www.r-project.org/) is a programming language specializing
in statistical computing and graphics. You can use R just the same as
any other programming languages, but it is most useful for statistical
analyses, with well-established packages for common tasks such as
[linear
modeling](https://cran.r-project.org/web/packages/lme4/index.html),
['omic data analysis](https://bioconductor.org/), [machine
learning](http://topepo.github.io/caret/index.html), and
[visualization](https://ggplot2.tidyverse.org/).
#### Installing R for the first time
You can download and install the free R software environment
[here](https://cloud.r-project.org/). Note that you should download the
latest release - this will ensure the R version you have is compatible
with MaAsLin 3.
#### (Optional) the RStudio IDE
[RStudio](https://rstudio.com/products/rstudio/) is a freely available
IDE (integrated development environment) for R. It is a "wrapper" around
R with some additional functionalities that make programming in R a bit
easier. Feel free to download RStudio and explore its interface - but it
is not required for this tutorial.
#### Important: the correct R version
If you already have R installed, then it is possible that the R version
you have does not support MaAsLin 3. The easiest way to check this is to
launch R/RStudio, and in the console ("Console" pane in RStudio), type
in the following command (without the `>` symbol):
```{R, eval=FALSE, cache = FALSE}
sessionInfo()
```
The first line of output message should indicate your current R version.
For example:
```
> sessionInfo()
R version 4.4.1 (2024-06-14)
```
For MaAsLin 3 to install, you will need R >= 4.4. If your version is
older than that, please refer to section [Installing R for the first
time](#installing-r-for-the-first-time) to download the latest R. Note
that RStudio users will need to have RStudio "switch" to this later
version once it is installed. This should happen automatically for
Windows and Mac OS users when you relaunch RStudio. For Linux users you
might need to bind the correct R executable.
Either way, once you have the correct version installed, launch the
software and use `sessionInfo()` to make sure that you indeed have R >=
4.4.
## 2. Installing MaAsLin 3 {#installing-maaslin-3}
The latest version of MaAsLin 3 can be installed from BiocManager.
```{r, eval=FALSE, cache = FALSE}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("biobakery/maaslin3")
```
To compile vignettes yourself, specify `dependencies = TRUE`.
```{r, eval=TRUE, cache = FALSE, echo=FALSE}
for (lib in c('maaslin3', 'dplyr', 'ggplot2', 'knitr', 'kableExtra')) {
suppressPackageStartupMessages(require(lib, character.only = TRUE))
}
```
## 3. Microbiome association detection with MaAsLin 3
To run MaAsLin 3, the user must supply a table of per-sample feature
abundances (with zeros still included), a table of per-sample metadata,
and a model specifying how the metadata should relate to the feature
prevalence (how likely the feature is to be present or absent) and
abundance (how much of the feature is there if it's there). MaAsLin 3
will return a table of associations including an effect size and p-value
for each feature-metadatum association and a folder of visuals including
a summary plot and diagnostic plots for significant associations.
### 3.1 MaAsLin 3 input {#maaslin-3-input}
MaAsLin3 requires two input files.
1. Feature abundance data frame
* Formatted with features as columns and samples as rows.
* The transpose of this format is also okay.
* Possible features include taxonomy or genes. These can be relative
abundances or counts.
* This can be a filepath to a tab-delimited file.
2. Metadata data frame
* Formatted with variables as columns and samples as rows.
* The transpose of this format is also okay.
* Possible metadata include gender or age.
* This can be a filepath to a tab-delimited file.
The data file can contain samples not included in the metadata file
(along with the reverse case). For both cases, those samples not
included in both files will be removed from the analysis. Also, the
samples do not need to be in the same order in the two files.
Samples with NA values for the metadata will be excluded when fitting the
models. It is recommended to drop these samples in advance under a
missing-at-random assumption or use multiple imputation to run
the complete dataset.
Example input files can be found in the ``inst/extdata`` folder of the
MaAsLin 3 source. The files provided were generated from the Human
Microbiome Project 2 (HMP2) data which can be downloaded from
https://ibdmdb.org/.
* ``HMP2_taxonomy.tsv``: a tab-delimited file with samples as rows and
species as columns. It is a subset of the full HMP2 taxonomy that
includes just some of the the species abundances.
* ``HMP2_metadata.tsv``: a tab-delimited file with samples as rows and
metadata as columns. It is a subset of the full HMP2 metadata that
includes just some of the fields.
```{R, cache = FALSE}
# Read abundance table
taxa_table_name <- system.file("extdata", "HMP2_taxonomy.tsv",
package = "maaslin3")
taxa_table <- read.csv(taxa_table_name, sep = '\t', row.names = 1)
# Read metadata table
metadata_name <- system.file("extdata", "HMP2_metadata.tsv",
package = "maaslin3")
metadata <- read.csv(metadata_name, sep = '\t', row.names = 1)
# Factor the categorical variables to test IBD against healthy controls
metadata$diagnosis <-
factor(metadata$diagnosis, levels = c('nonIBD', 'UC', 'CD'))
metadata$dysbiosis_state <-
factor(metadata$dysbiosis_state, levels =
c('none', 'dysbiosis_UC', 'dysbiosis_CD'))
metadata$antibiotics <-
factor(metadata$antibiotics, levels = c('No', 'Yes'))
taxa_table[1:5, 1:5]
metadata[1:5, 1:5]
```
### 3.2 Running MaAsLin 3 {#running-maaslin-3}
MaAsLin 3 is run by specifying the
abundance table (`input_data`), the metadata table (`input_metadata`),
the output directory (`output`), and a model. The model can come from a
formula or vectors of terms. In either case, variable names should not
have spaces or unusual characters.
* *Formula*: The `formula` parameter should be set to any formula that
satisfies the `lme4` specifications: fixed effects, random effects
([section 4.2](#random-effects)), interaction terms ([section
4.3](#interactions)), polynomial terms, and more can all be included. If
categorical variables are included as fixed effects, each level will be
tested against the first factor level. In addition, ordered predictors
([section 4.4](#level-contrasts)) and group predictors ([section
4.5](#group-wise-differences)) can be included by including
`group(variable_name)` and `ordered(variable_name)` in the formula. When
the data involve paired or matched case-control samples, conditional
logistic regression can be run by specifying the grouping variable with
`strata(variable_name)`.
* *Vectors*: Alternatively, a vector of variable names can be supplied
to the parameters `fixed_effects`, `random_effects`, `group_effects`,
`ordered_effects`, and `strata_effects`. Categorical variables
should either be ordered as
factors beforehand, or `reference` should be provided as a string of
'variable,reference' semi-colon delimited for multiple variables (e.g.,
`variable_1,reference_1;variable_2,reference_2`). NOTE: adding a space
between the variable and level might result in the wrong reference level
being used.
**Because MaAsLin 3 identifies prevalence (presence/absence)
associations, sample read depth (number of reads) should be included as
a covariate if available. Deeper sequencing will likely increase feature
detection in a way that could spuriously correlate with metadata of
interest when read depth is not included in the model.**
#### Running MaAsLin 3 on HMP2 data
The following command runs MaAsLin 3 on the HMP2 data, running a
multivariable
regression model to test for the association between microbial species
abundance and prevalence
versus IBD diagnosis and [dysbiosis
scores](https://www.nature.com/articles/s41586-019-1237-9) while
controlling for antibiotic usage, age, and the sampling depth (reads).
The abundance associations come from fitting multivariate linear models
to the log base 2 relative abundances of features in samples in which
the feature is present. The prevalence associations come from fitting
multivariate logistic models to the presence/absence data for a feature
across all samples. Outputs are directed to a folder called
`hmp2_output` under the current working directory (`output =
"hmp2_output"`).
To show one of each type of plot (see below), the maximum number of
significant associations to plot has been increased with `max_pngs=100`
(default `30`). All other parameters beyond the first four are the
default parameters but are still included for clarity:
* Total sum scaling (`normalization = 'TSS'`) with a log transformation
(`'transform = 'LOG'`, base 2) afterwards will be used. These are almost
always the recommended choices (and all MaAsLin 3 evaluations were
performed with these options), but other normalizations and
transformations are allowed (see `normalization` and `transform` in
`?maaslin3`).
* A data augmentation procedure is used to avoid linear separability in
the logistic regressions (`augment = TRUE`). This is almost always
recommended, though it can be turned off (see `augment` in the manual).
* Continuous metadatum variables are z-score standardized (subtract the
mean, divide by the standard deviation) with `standardize = TRUE` so
that the coefficients will be on the same scale (improving
visualization).
* A nominal FDR level of 0.1 (`max_significance = 0.1`) determines which
associations are significant, but this can also be changed (see
`max_significance` in the manual).
* For the abundance coefficients, to account for compositionality in
relative abundance data, significance is determined by comparing against
the median coefficient for a metadatum across the features
(`median_comparison_abundance`). Since prevalence coefficients do not
have the same compositional properties, they are instead compared
against 0 (`median_comparison_prevalence`). See below for a discussion of
these parameters.
* One CPU is used for fitting (`cores = 1`).
The abundance, prevalence, and variance filtering parameters are not
included, and they are 0 by default. **Different from other differential
abundance tools, low prevalence features need not be filtered out**
since the prevalence modeling in MaAsLin 3 already accounts for high
proportions of zeros. However, filtering low prevalence features might improve
power.
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
set.seed(1)
fit_out <- maaslin3(input_data = taxa_table,
input_metadata = metadata,
output = 'hmp2_output',
formula = '~ diagnosis + dysbiosis_state +
antibiotics + age + reads',
normalization = 'TSS',
transform = 'LOG',
augment = TRUE,
standardize = TRUE,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 10,
cores = 1)
```
By default, many lines of information are printed while the models fit,
but the verbosity can be changed with the `verbosity` argument (e.g.,
`verbosity = 'WARN'`).
MaAsLin 3 can also take as input a `SummarizedExperiment` or
`TreeSummarizedExperiment` object from BioConductor:
```{R, echo = TRUE, results = 'hide', warning = FALSE, eval = FALSE}
se <- SummarizedExperiment(
assays = list(taxa_table = t(taxa_table)),
colData = metadata
)
fit_out <- maaslin3(input_data = se,
output = 'hmp2_output',
formula = '~ diagnosis + dysbiosis_state +
antibiotics + age + reads',
normalization = 'TSS',
transform = 'LOG',
augment = TRUE,
standardize = TRUE,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 10,
cores = 1)
```
#### Median comparisons
When `median_comparison_abundance` or `median_comparison_prevalence` are
on, the coefficients for a metadatum will be tested against the median
coefficient for that metadatum (median across the features). Otherwise,
the coefficients will be tested against 0. For abundance associations,
this is designed to account for compositionality, the idea that
**testing against zero on the relative scale is not equivalent to
testing against zero on the absolute scale**. The user manual contains a
more extensive discussion of when to use `median_comparison_abundance`
and `median_comparison_prevalence`, but in summary:
* `median_comparison_abundance` is `TRUE` by default and should be used
to make inference on the absolute scale when using relative abundance
data. When `median_comparison_abundance` is `TRUE`, only the p-values
and q-values change; the coefficients returned are still the relative
abundance coefficients unless `subtract_median` is `TRUE` in which case
the median will be subtracted from the coefficients.
* `median_comparison_abundance` should be `FALSE` when (1) testing for
significant relative associations, (2) testing for absolute abundance
associations under the assumption that the total absolute abundance is
not changing, or (3) testing for significant absolute associations when
supplying spike-in or total abundances with `unscaled_abundance`.
* `median_comparison_prevalence` is `FALSE` by default.
### 3.3 MaAsLin 3 output {#maaslin-3-output}
#### Significant associations
The main output from MaAsLin 3 is the list of significant associations
in `significant_results.tsv`. This lists all associations with joint
or individual q-values
that pass MaAsLin 3's significance
threshold, ordered by increasing individual q-value. The format is:
* `feature` and `metadata` are the feature and metadata names.
* `value` and `name` are the value of the metadata and variable name
from the model.
* `coef` and `stderr` are the fit coefficient and standard error from
the model. In abundance models, a one-unit change in the metadatum
variable corresponds to a $2^{\textrm{coef}}$ fold change in the
relative abundance of the feature. In prevalence models, a one-unit
change in the metadatum variable corresponds to a $\textrm{coef}$ change
in the log-odds of a feature being present.
* `null_hypothesis` is the value of the null hypothesis against which
the coefficient is tested (zero or the per-metadatum median).
* `pval_individual` and `qval_individual` are the p-value and false
discovery rate (FDR) corrected q-value of the individual association.
FDR correction is performed over all p-values without errors in the abundance
and prevalence modeling together.
* `pval_joint` and `qval_joint` are the p-value and q-value of the joint
prevalence and abundance association. The p-value comes from plugging in
the minimum of the association's abundance and prevalence p-values into
the Beta(1,2) CDF. It is interpreted as the probability that either the
abundance or prevalence association would be as extreme as observed if
there was neither an abundance nor prevalence association between the
feature and metadatum.
* `model` specifies whether the association is abundance or prevalence.
* `N` and `N_not_zero` are the total number of data points and the total
number of non-zero data points for the feature.
```{r, echo = FALSE, cache = FALSE}
signif_results <- read.csv('hmp2_output/significant_results.tsv',
sep='\t')
head(signif_results, 20) %>%
dplyr::mutate_if(is.numeric, .funs = function(x){(
as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
#### Full output file structure
MaAsLin 3 generates two types of output files explained below: data and
visualizations. In addition, the object returned from `maaslin3`
(`fit_out` above) contains all the data and results (see
`?maaslin_fit`).
1. Data output files
* ``significant_results.tsv``
* Described above.
* ``all_results.tsv``
* Same format as above but for all associations and with an extra column
listing any errors from the model fitting.
* ``features``
* This folder includes the filtered, normalized, and transformed
versions of the input feature table.
* These steps are performed sequentially in the above order.
* If an option is set such that a step does not change the data, the
resulting table will still be output.
* ``models_linear.rds`` and ``models_logistic.rds``
* These files contain a list with every model fit object (`linear` for
linear models, `logistic` for logistic models).
* It will only be generated if `save_models` is set to TRUE.
* ``residuals_linear.rds`` and ``residuals_logstic.rds``
* These files contain a data frame with residuals for each feature.
* ``fitted_linear.rds`` and ``fitted_logistic.rds``
* These files contain a data frame with fitted values for each feature.
* ``ranef_linear.rds`` and ``ranef_logistic.rds``
* These files contain a data frame with extracted random effects for
each feature (when random effects are specified).
* ``maaslin3.log``
* This file contains all log information for the run.
* It includes all settings, warnings, errors, and steps run.
2. Visualization output files
* ``summary_plot.pdf``
* This file contain a combined coefficient plot and heatmap of the most
significant associations. In the heatmap, one star indicates the
individual q-value is below the parameter `max_significance`, and two
stars indicate the individual q-value is below `max_significance`
divided by 10. (Make sure to set `include=T` to view the plots when
knitting the Rmd.)
```{r, out.width='100%', echo=FALSE, cache = FALSE, include=FALSE, eval=FALSE}
# Rename summary plot to avoid knitting issues later
quiet_out <- file.rename('hmp2_output/figures/summary_plot.png',
'hmp2_output/figures/summary_plot_first.png')
knitr::include_graphics("hmp2_output/figures/summary_plot_first.png")
```
* ``association_plots/[metadatum]/[association]/
[metadatum]_[feature]_[association].png``
* A plot is generated for each significant association up to `max_pngs`.
* Scatter plots are used for continuous metadata abundance associations.
* Box plots are used for categorical data abundance associations.
* Box plots are used for continuous data prevalence associations.
* Grids are used for categorical data prevalence associations.
* Data points plotted are after filtering, normalization, and
transformation so that the scale in the plot is the scale that was used
in fitting.
* **To see the association plots, set `max_pngs = 250` above and set
`eval=TRUE` and `include=TRUE` in this chunk**
```{r, echo=FALSE, fig.show='hold',include=FALSE, eval=FALSE}
prefix <- "hmp2_output/figures/association_plots"
plot_vec <-
c("/age/linear/age_Enterocloster_clostridioformis_linear.png",
"/dysbiosis_state/linear/dysbiosis_state_Escherichia_coli_linear.png",
"/age/logistic/age_Bifidobacterium_longum_logistic.png",
paste0("/dysbiosis_state/logistic/",
"dysbiosis_state_Faecalibacterium_prausnitzii_logistic.png"
))
knitr::include_graphics(paste0(prefix, plot_vec))
```
At the top right of each association plot is the name of the significant
association in the results file, the FDR corrected q-value for the
individual association, the number of samples in the dataset, and the
number of samples with non-zero abundances for the feature. In the plots
with categorical metadata variables, the reference category is on the
left, and the significant q-values and coefficients in the top right are
in the order of the values specified above. Because the displayed
coefficients correspond to the full fit model with (possibly) scaled
metadata variables, the marginal association plotted might not match the
coefficient displayed. However, the plots are intended to provide an
interpretable visual while usually agreeing with the full model.
#### Diagnostics
There are a few common issues to check for in the results:
1. When warnings or errors are thrown during the fitting process, they
are recorded in the `error` column of `all_results.tsv`. Often, these
indicate model fitting failures or poor fits that should not be trusted,
but sometimes the warnings can be benign, and the model fit might still
be reasonable. Users should check associations of interest if they
produce errors.
* With `warn_prevalence=TRUE`, a note will be produced in the error column
indicating when prevalence associations are likely induced by abundance
associations. These should be checked visually with the diagnosic visuals.
2. Despite the error checking, significant results could still result
from poor model fits. These can usually be diagnosed with the visuals in
the `association_plots` directory.
* Any significant abundance associations with a categorical variable
should usually have **at least 10 observations in each category**.
* Significant prevalence associations with categorical variables should
also have **at least 10 samples in which the feature was present and at
least 10 samples in which it was absent for each significant category**.
* Significant abundance associations with continuous metadata should be
checked visually for influential outliers.
3. The q-values are FDR corrected over all abundance and prevalence
relationships, so it may be preferable to FDR correct just
the p-values from the variables of interest. This can reduce false
positives when there are many significant but uninteresting associations
(e.g., many read depth associations).
4. There are also a few rules of thumb to keep in mind.
* Models should ideally have about 10 times as many samples (all samples
for logistic fits, non-zero samples for linear fits) as covariate terms
(all continuous variables plus all categorical variable levels).
* Coefficients (effect sizes) larger than about 15 in absolute value are
usually suspect unless very small unstandardized predictors are being
included. (A coefficient of 15 corresponds to a fold change >30000!). If
you encounter such coefficients, check that (1) no error was thrown, (2)
the diagnostics look reasonable, (3) a sufficient number of samples were
used in fitting, (4) the q-value is significant, (5) the metadata are
not highly collinear, and (6) the random effects are plausible.
#### Replotting
Once `maaslin3` has been run once, `maaslin_plot_results_from_output`
can be run to
(re-)create the plots. This allows the user to plot the associations
even without having the R object returned by `maaslin_fit` or `maaslin3`
(e.g., if fitting models through the command line). It is recommended to
fit models with simple variables names that are robust to formula
plotting and then convert these into proper names for plotting.
Likewise, `heatmap_vars` and `coef_plot_vars` can be specified at any point,
but it is usually easier to see how the names
come out first and then choose which metadata variables will go in the
coefficient plot and which will go in the heatmap afterwards with
`maaslin_plot_results_from_output`.
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
# This section is necessary for updating the
# summary plot and the association plots
# Rename results file with clean titles
all_results <- read.csv('hmp2_output/all_results.tsv', sep='\t')
all_results <- all_results %>%
mutate(metadata = case_when(metadata == 'age' ~ 'Age',
metadata == 'antibiotics' ~ 'Abx',
metadata == 'diagnosis' ~ 'Diagnosis',
metadata == 'dysbiosis_state' ~ 'Dysbiosis',
metadata == 'reads' ~ 'Read depth'),
value = case_when(value == 'dysbiosis_CD' ~ 'CD',
value == 'dysbiosis_UC' ~ 'UC',
value == 'Yes' ~ 'Used', # Antibiotics
value == 'age' ~ 'Age',
value == 'reads' ~ 'Read depth',
TRUE ~ value),
feature = gsub('_', ' ', feature) %>%
gsub(pattern = 'sp ', replacement = 'sp. '))
# Write results
write.table(all_results, 'hmp2_output/all_results.tsv', sep='\t')
# Set the new heatmap and coefficient plot variables and order them
heatmap_vars = c('Dysbiosis UC', 'Diagnosis UC',
'Abx Used', 'Age', 'Read depth')
coef_plot_vars = c('Dysbiosis CD', 'Diagnosis CD')
# This section is necessary for updating the association plots
taxa_table_copy <- taxa_table
colnames(taxa_table_copy) <- gsub('_', ' ', colnames(taxa_table_copy)) %>%
gsub(pattern = 'sp ', replacement = 'sp. ')
# Rename the features in the norm transformed data file
data_transformed <-
read.csv('hmp2_output/features/data_transformed.tsv', sep='\t')
colnames(data_transformed) <-
gsub('_', ' ', colnames(data_transformed)) %>%
gsub(pattern = 'sp ', replacement = 'sp. ')
write.table(data_transformed,
'hmp2_output/features/data_transformed.tsv',
sep='\t', row.names = FALSE)
# Rename the metadata like in the outputs table
metadata_copy <- metadata
colnames(metadata_copy) <-
case_when(colnames(metadata_copy) == 'age' ~ 'Age',
colnames(metadata_copy) == 'antibiotics' ~ 'Abx',
colnames(metadata_copy) == 'diagnosis' ~ 'Diagnosis',
colnames(metadata_copy) == 'dysbiosis_state' ~ 'Dysbiosis',
colnames(metadata_copy) == 'reads' ~ 'Read depth',
TRUE ~ colnames(metadata_copy))
metadata_copy <- metadata_copy %>%
mutate(Dysbiosis = case_when(Dysbiosis == 'dysbiosis_UC' ~ 'UC',
Dysbiosis == 'dysbiosis_CD' ~ 'CD',
Dysbiosis == 'none' ~ 'None') %>%
factor(levels = c('None', 'UC', 'CD')),
Abx = case_when(Abx == 'Yes' ~ 'Used',
Abx == 'No' ~ 'Not used') %>%
factor(levels = c('Not used', 'Used')),
Diagnosis = case_when(Diagnosis == 'nonIBD' ~ 'non-IBD',
TRUE ~ Diagnosis) %>%
factor(levels = c('non-IBD', 'UC', 'CD')))
# Recreate the plots
scatter_plots <- maaslin_plot_results_from_output(
output = 'hmp2_output',
metadata = metadata_copy,
normalization = "TSS",
transform = "LOG",
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_significance = 0.1,
max_pngs = 20)
```
In the new summary plot below, we can see that the feature names are
cleaned up, the metadata names are cleaned up, the set of metadata
variables used in the coefficient plot is different, and the metadata
used in the heatmap is reordered.
```{r, out.width='100%', echo=FALSE, cache = FALSE, include=FALSE}
knitr::include_graphics("hmp2_output/figures/summary_plot.png")
```
In the association plots, the taxa and metadata have been renamed to be
consistent with the results file from earlier.
```{r, echo=FALSE, fig.show='hold', cache = FALSE, include=FALSE, eval=FALSE}
prefix <- "hmp2_output/figures/association_plots"
plot_vec <- c("/Age/linear/Age_Enterocloster clostridioformis_linear.png",
"/Dysbiosis/linear/Dysbiosis_Escherichia coli_linear.png",
"/Age/logistic/Age_Bifidobacterium longum_logistic.png",
paste0("/Dysbiosis/logistic/",
"Dysbiosis_Faecalibacterium prausnitzii_logistic.png"
))
knitr::include_graphics(paste0(prefix, plot_vec))
```
## 4. Advanced Topics {#advanced-topics}
### 4.1 Absolute abundance {#absolute-abundance}
Most microbiome methodologies have historically focused on relative
abundances (proportions out of 1). However, some experimental protocols
can enable estimation of absolute abundances (cell count/concentration).
MaAsLin 3 can be used with two types of absolute abundance estimation:
spike-ins and total abundance scaling. In a spike-in procedure, a small,
known quantity of a microbe that otherwise would not be present in the
sample is added, and the sequencing procedure is carried out as usual.
Then, the absolute abundance of a microbe already in the community is
estimated as:
$$\textrm{Absolute abundance other microbe}=\frac{\textrm{Relative
abundance other microbe}}{\textrm{Relative abundance spike-in
microbe}}\cdot (\textrm{Absolute abundance spike-in microbe})$$
Alternatively, the total microbial abundance of a sample can be
determined (e.g., with qPCR of a marker gene or by cell counting). Then,
the absolute abundance of a microbe in the community is estimated as:
$$\textrm{Absolute abundance microbe}=(\textrm{Total absolute
abundance})\cdot(\textrm{Relative abundance microbe})$$
#### Spike-in
The following section shows a synthetic spike-in procedure with
simulated data from SparseDOSSA2. The abundance table is like any other
abundance table input to MaAsLin 3 except that 'Feature101' is the
spike-in (which must be present in every sample). The scaling factor
data frame has 'Feature101' as its only column name with the absolute
abundance of the spike-in for each sample (rownames). If the same
quantity of the spike-in was added to all samples, this column could be
entirely 1 (or any other arbitrary value). When using a spike-in
procedure, the scaling factor data frame should have a single column
named the same as one of the features in the abundance table, which will
be used as the spike-in.
```{R, cache = FALSE}
# Abundance table
taxa_table_name <- system.file("extdata", "abundance_spike_in_ex.tsv",
package = "maaslin3")
spike_in_taxa_table <- read.csv(taxa_table_name, sep = '\t', row.names = 1)
# Metadata table
metadata_name <- system.file("extdata", "metadata_spike_in_ex.tsv",
package = "maaslin3")
spike_in_metadata <- read.csv(metadata_name, sep = '\t', row.names = 1)
for (col in c('Metadata_1', 'Metadata_2', 'Metadata_5')) {
spike_in_metadata[,col] <- factor(spike_in_metadata[,col])
}
# Spike-in table
unscaled_name <- system.file("extdata",
"scaling_factors_spike_in_ex.tsv",
package = "maaslin3")
spike_in_unscaled <- read.csv(unscaled_name, sep = '\t', row.names = 1)
spike_in_taxa_table[c(1:5, 101),1:5]
spike_in_metadata[1:5,]
spike_in_unscaled[1:5, , drop=FALSE]
```
The following code fits the absolute abundance model. Note that the
normalization must be set to `TSS` since the procedure involves scaling
the relative abundance ratios. Also, `median_comparison_abundance` is
now set to `FALSE` since we want to test the absolute coefficients
against zero. The data frame of absolute abundances is included as the
`unscaled_abundance` parameter, and the spike-in strategy will
automatically be detected based on the column name.
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
fit_out <- maaslin3(
input_data = spike_in_taxa_table,
input_metadata = spike_in_metadata,
output = 'spike_in_demo',
formula = '~ Metadata_1 + Metadata_2 + Metadata_3 +
Metadata_4 + Metadata_5',
normalization = 'TSS',
transform = 'LOG',
median_comparison_abundance = FALSE,
unscaled_abundance = spike_in_unscaled)
```
The absolute abundance coefficients can be inspected:
```{R, cache = FALSE}
rownames(fit_out$fit_data_abundance$results) <- NULL
head(fit_out$fit_data_abundance$results, 20) %>%
dplyr::mutate_if(is.numeric, .funs =
function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
#### Total abundance scaling
The following section shows a synthetic total abundance scaling
procedure with simulated data from SparseDOSSA2. The abundance table is
like any other abundance table input to MaAsLin 3 without any extra
features. The scaling factor data frame has 'total' as its only column
name with the total absolute abundance for each sample (rownames).
```{r, cache = FALSE}
# Abundance table
taxa_table_name <- system.file("extdata", "abundance_total_ex.tsv",
package = "maaslin3")
total_scaling_taxa_table <- read.csv(taxa_table_name, sep = '\t', row.names = 1)
# Metadata table
metadata_name <- system.file("extdata", "metadata_total_ex.tsv",
package = "maaslin3")
total_scaling_metadata <- read.csv(metadata_name, sep = '\t', row.names = 1)
for (col in c('Metadata_1', 'Metadata_3', 'Metadata_5')) {
spike_in_metadata[,col] <- factor(spike_in_metadata[,col])
}
# Total abundance table
unscaled_name <- system.file("extdata", "scaling_factors_total_ex.tsv",
package = "maaslin3")
total_scaling_unscaled <- read.csv(unscaled_name, sep = '\t', row.names = 1)
total_scaling_taxa_table[1:5, 1:5]
total_scaling_metadata[1:5,]
total_scaling_unscaled[1:5, , drop=FALSE]
```
The following code fits the absolute abundance model. As before, `TSS`
must be set to `TRUE`, `median_comparison_abundance` is set to `FALSE`,
and the data frame of absolute abundances is included as the
`unscaled_abundance` parameter. The total abundance scaling procedure
will be detected based on the column name.
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
fit_out <- maaslin3(
input_data = total_scaling_taxa_table,
input_metadata = total_scaling_metadata,
output = 'total_scaling_demo',
formula = '~ Metadata_1 + Metadata_2 + Metadata_3 +
Metadata_4 + Metadata_5',
normalization = 'TSS',
transform = 'LOG',
median_comparison_abundance = FALSE,
unscaled_abundance = total_scaling_unscaled)
```
The absolute abundance coefficients can be inspected:
```{R, cache = FALSE}
rownames(fit_out$fit_data_abundance$results) <- NULL
head(fit_out$fit_data_abundance$results, n = 20) %>%
dplyr::mutate_if(is.numeric,
.funs = function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
#### Computationally estimated absolute abundance
When using a purely computational approach to estimate absolute abundance
(e.g., Nishijima et al. Cell 2024 with the function `MLP`), the abundance
table can be estimated
and then run with the `normalization = 'NONE'` to avoid scaling the results.
### 4.2 Random effects (clustered samples) {#random-effects}
Certain studies have a natural "grouping" of sample observations, such
as by subject in longitudinal designs or by family in family designs. It
is important for statistic analysis to address the non-independence
between samples belonging to the same group. MaAsLin 3 provides a simple
interface for this with random effects, where the user can specify the
grouping variable to run a [mixed effect
model](https://cran.r-project.org/web/packages/lme4/index.html). This
grouping variable can be provided with the `random_effects` parameter or
specified in the model with `(1|grouping_variable)`. This will add in a
"random intercept," a per-group adjustment to the model intercept. More
complicated random effects can be specified in the formula in accordance
with `lme4` formula parsing. For example, HMP2 has a longitudinal design
where the same subject (column `participant_id`) can have multiple
samples. We thus use `participant_id` as a random effect grouping
variable in the model.
```{R, echo = TRUE, results = 'hide', warning = FALSE, messages = FALSE}
# Subset to only CD cases for time; taxa are subset automatically
fit_out <- maaslin3(
input_data = taxa_table,
input_metadata = metadata[metadata$diagnosis == 'CD',],
output = 'random_effects_output',
formula = '~ dysbiosis_state + antibiotics +
age + reads + (1|participant_id)',
plot_summary_plot = FALSE,
plot_associations = FALSE)
```
```{r, echo = FALSE, cache = FALSE}
signif_results <-
read.csv('random_effects_output/significant_results.tsv',
sep='\t')
head(signif_results, n = 20) %>%
dplyr::mutate_if(is.numeric, .funs =
function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
Note that no `participant_id` terms are included in the outputs; the
random intercepts are just used to control for grouping. If you are
interested in testing the effect of time in a longitudinal study, the
time point variable should be included as a fixed effect in your MaAsLin
3 formula.
#### Alternatives to random effects
Random effects are not the solution to every correlated data problem. In
particular, **logistic random intercept models may produce poor model fits
when there are fewer than 4 average observations per random effect group.** The
fit coefficients can be large in magnitude and very significant while
modeling the data poorly. By default, a warning will be added to the `error`
column of the outputs in these scenarios. This warning can be disregarded
if the user is confident the model is fit correctly, but recommendations
are provided in the table below:
| Scenario |$\geq4$ average observations per group|$<4$ average observations\
per group, metadata not collinear with the groups |$<4$ average observations\
per group, metadata collinear with the groups|
|----|----|-----|------|
| Example | Subjects sampled $\geq4$ times | Subjects sampled once before \
and once after treatment | Subjects have fixed disease status \
and are sampled twice (disease \
status is constant per-subject, so collinear with the grouping) |
| Action | Use random effects with `small_random_effects=FALSE` \
(default) | Use random effects for abundance but fixed effects for \
prevalence by setting `small_random_effects=TRUE` | Subset or average data to \
one observation per group |
| Reason | Random effects are robust and \
increase power | Prevalence random effects are not robust but fixed effects \
are robust and control for grouping | Prevalence random effects are \
not robust and fixed effects cannot be used when groups are \
collinear with the metadata |
### 4.3 Interactions (differences in differences) {#interactions}
One of the benefits of MaAsLin 3's formula mode is the ability to
include interaction terms. Mathematically, an interaction term
corresponds to the product of two columns in the design matrix. When a
continuous variable is interacted with a categorical term, the
interaction term corresponds to the change in the continuous variable's
slope between the categories. For two categorical variables interacted,
see below; it is better explained through an example. In the formula,
interactions can be specified with the `:` symbol to include only the
interaction term(s) or the `*` symbol to include both the interaction
term and the non-interacted terms (i.e., `a*b` adds the terms `a`, `b`,
and `a:b`).
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
metadata <- metadata %>%
mutate(dysbiosis_general = ifelse(dysbiosis_state != 'none',
'dysbiosis', 'none')) %>%
mutate(dysbiosis_general = factor(dysbiosis_general, levels =
c('none', 'dysbiosis')))
fit_out <- maaslin3(
input_data = taxa_table,
input_metadata = metadata,
output = 'interaction_output',
formula = '~ diagnosis + diagnosis:dysbiosis_general +
antibiotics + age + reads')
```
```{R, cache = FALSE}
full_results <- rbind(fit_out$fit_data_abundance$results,
fit_out$fit_data_prevalence$results)
full_results <- full_results %>%
dplyr::arrange(qval_joint) %>%
dplyr::filter(metadata == "diagnosis")
rownames(full_results) <- NULL
head(full_results, n = 20) %>%
dplyr::mutate_if(is.numeric,
.funs = function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
In the model above, we have converted `dysbiosis_CD` and `dysbiosis_UC`
into `dysbiosis` in the variable `dysbiosis_general` and then specified
an interaction between `diagnosis` and `dysbiosis_general` with
`diagnosis:dysbiosis_general`. Since `diagnosis` has 3 levels itself
(`nonIBD`, `UC`, `CD`), this will produce five terms (`name` column) for
each feature:
* `diagnosisCD` is the difference between non-dysbiotic CD and non-IBD
* `diagnosisUC` is the difference between non-dysbiotic UC and non-IBD
* `diagnosisCD:dysbiosis_generaldysbiosis` (the interaction of
`diagnosisCD` and `dysbiosis_generaldysbiosis`) is the difference
between non-dysbiotic CD and dysbiotic CD
* `diagnosisnonIBD:dysbiosis_generaldysbiosis` (the interaction of
`diagnosisnonIBD` and `dysbiosis_generaldysbiosis`) is the difference
between dysbiotic and non-dysbiotic non-IBD. Note that these
coefficients are all NA since no subjects labeled non-IBD were marked as
being dysbiotic.
* `diagnosisUC:dysbiosis_generaldysbiosis` (the interaction of
`diagnosisUC` and `dysbiosis_generaldysbiosis`) is the difference
between non-dysbiotic UC and dysbiotic UC
Since we have effectively fit the original model in a different way,
this output table is equivalent to the original output table but with
new names in the columns `metadata`, `value`, and `name` based on the
interaction terms.
### 4.4 Level contrasts {#level-contrasts}
Another feature of MaAsLin 3 is the ability to test for
level-versus-level differences in ordered predictors with contrast
tests. Ordered predictors are categorical variables with a natural
ordering such as cancer stage, consumption frequency of a dietary
factor, or dosage group. Here, we assess how microbial abundances and
prevalences are associated with eating red meat by including
`ordered(red_meat)` in the formula. This will perform a contrast test
for whether there is a difference between each pair of subsequent levels
(e.g., "Yesterday, 3 or more times" versus "Yesterday, 1 to 2 times")
rather than whether there are differences between the levels and the
baseline (e.g., "Yesterday, 3 or more times" versus "Not in the last 7
days"). The coefficient, standard error, and p-value all correspond to
the difference between the level in `value` and the previous level.
Ordered predictors should only be included as fixed effects (i.e., no
ordered predictors as random effects etc.).
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
# Put the red meat consumption responses in order
metadata <- metadata %>%
mutate(red_meat = ifelse(
red_meat == 'No, I did not consume these products in the last 7 days',
'Not in the last 7 days',
red_meat) %>%
factor(levels = c('Not in the last 7 days',
'Within the past 4 to 7 days',
'Within the past 2 to 3 days',
'Yesterday, 1 to 2 times',
'Yesterday, 3 or more times'))
)
# Create the model with only non-IBD subjects
fit_out <- maaslin3(
input_data = taxa_table,
input_metadata = metadata[metadata$diagnosis == 'nonIBD',],
output = 'ordered_outputs',
formula = '~ ordered(red_meat) + antibiotics + age + reads',
plot_summary_plot = TRUE,
plot_associations = TRUE,
heatmap_vars = c('red_meat Within the past 4 to 7 days',
'red_meat Within the past 2 to 3 days',
'red_meat Yesterday, 1 to 2 times',
'red_meat Yesterday, 3 or more times'),
max_pngs = 30)
```
```{r, out.width='100%', echo=FALSE, cache = FALSE, include=FALSE}
knitr::include_graphics("ordered_outputs/figures/summary_plot.png")
```
If we look at the resulting heatmap, we can identify the *Alistipes
shahii* prevalence association as potentially interesting since it
increases in prevalence with more meat consumption in all but one
level-versus-level comparison.
```{r, echo = FALSE, cache = FALSE}
full_results <- rbind(fit_out$fit_data_abundance$results,
fit_out$fit_data_prevalence$results)
full_results <- full_results %>%
dplyr::filter(metadata %in% c("red_meat") &
feature == 'Alistipes_shahii' &
model == 'logistic') %>%
dplyr::mutate(
value =
factor(
value,
levels =
c('No, I did not consume these products in the last 7 days',
'Within the past 4 to 7 days',
'Within the past 2 to 3 days',
'Yesterday, 1 to 2 times',
'Yesterday, 3 or more times'))) %>%
dplyr::arrange(value)
rownames(full_results) <- NULL
head(full_results, n = 20) %>%
dplyr::mutate_if(is.numeric,
.funs = function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "200px")
```
### 4.5 Group-wise differences {#group-wise-differences}
MaAsLin 3 can also test
for group-wise differences in categorical predictors using an ANOVA or
ANOVA-like procedure. Group-wise predictors are categorical variables
with or without an ordering such as race, country, or consumption
frequency of a dietary factor. When performing a group-wise difference
test, we are not comparing any two particular levels; rather, we are
investigating whether abundances and prevalences are the same for all
levels. Therefore, no coefficient is returned, only a p-value. Here, we
test whether there are differences in microbial abundances and
prevalences across people with different levels of red meat consumption
by including `group(red_meat)`. Group-wise predictors should only be
included as fixed effects (i.e., no group-wise predictors as random
effects etc.).
```{R, echo = T, results = 'hide', warning = FALSE, cache = FALSE}
fit_out <- maaslin3(
input_data = taxa_table,
input_metadata = metadata[metadata$diagnosis == 'nonIBD',],
output = 'group_outputs',
formula = '~ group(red_meat) + antibiotics + age + reads',
plot_summary_plot = TRUE,
plot_associations = TRUE,
heatmap_vars = c('red_meat Within the past 4 to 7 days',
'red_meat Within the past 2 to 3 days',
'red_meat Yesterday, 1 to 2 times',
'red_meat Yesterday, 3 or more times'),
max_pngs = 200)
```
If we look at the same *Alistipes shahii* association from before, we
see that no coefficient or standard error is returned for a group
predictor, but the q-value is very low. This corroborates the
observation earlier that there are differences in *Alistipes shahii*
prevalence with red meat consumption.
```{r, echo = FALSE, cache = FALSE}
full_results <- rbind(fit_out$fit_data_abundance$results,
fit_out$fit_data_prevalence$results)
full_results <- full_results %>%
dplyr::filter(metadata %in% c("red_meat") &
feature == 'Alistipes_shahii' &
model == 'logistic') %>%
dplyr::mutate(
value =
factor(
value,
levels =
c('No, I did not consume these products in the last 7 days',
'Within the past 4 to 7 days',
'Within the past 2 to 3 days',
'Yesterday, 1 to 2 times',
'Yesterday, 3 or more times'))) %>%
dplyr::arrange(value)
rownames(full_results) <- NULL
head(full_results, n = 20) %>%
dplyr::mutate_if(is.numeric,
.funs = function(x){(as.character(signif(x, 3)))}) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "100px")
```
### 4.6 Contrast tests {#contrast-tests}
We can also perform contrast testing in MaAsLin 3 to check whether a linear
combination of the coefficients is significantly different from some constant
on the right hand size (RHS). Specify a contrast matrix
with column names matching the coefficient names and rows with the contrasts
of interest. For each row $C^T$, the hypothesis $H_0:C^T\vec{\beta}=RHS$
will be tested. If no RHS is specified, a median test will be performed,
and the median test RHS will be returned.
Below, we test whether the UC and CD diagnosis coefficients are the same
and whether the UC and CD dysbiosis coefficients are the same. The contrast
matrix can easily be extended to test any groups pairwise.
```{R, echo = TRUE, results = 'hide', warning = FALSE, cache = FALSE}
fit_out <- maaslin3(input_data = taxa_table,
input_metadata = metadata,
output = 'contrast_test_output',
formula = '~ diagnosis + dysbiosis_state +
antibiotics + age + reads',
plot_summary_plot = FALSE,
plot_associations = FALSE,
cores = 1)
```
```{R, echo = TRUE, warning = FALSE, cache = FALSE}
contrast_mat <- matrix(c(1, -1, 0, 0, 0, 0, 1, -1),
ncol = 4, nrow = 2, byrow = TRUE)
colnames(contrast_mat) <- c("diagnosisUC",
"diagnosisCD",
"dysbiosis_statedysbiosis_UC",
"dysbiosis_statedysbiosis_CD")
rownames(contrast_mat) <- c("diagnosis_test", "dysbiosis_test")
contrast_mat
contrast_out <- maaslin_contrast_test(maaslin3_fit = fit_out,
contrast_mat = contrast_mat)
head(contrast_out$fit_data_abundance$results, n = 20) %>%
knitr::kable() %>%
kableExtra::kable_styling("striped", position = 'center') %>%
kableExtra::scroll_box(width = "800px", height = "400px")
```
## 5. Command line #### {#command-line}
MaAsLin 3 can also be run with a command line interface. For example,
the first HMP2 analysis can be performed with the following command (the
slashes may need to be removed):
```{r, engine = 'bash', eval = FALSE, cache = FALSE}
./R/maaslin3.R \
inst/extdata/HMP2_taxonomy.tsv \
inst/extdata/HMP2_metadata.tsv \
command_line_output \
--formula='~ diagnosis + dysbiosis_state + antibiotics + age + reads' \
--reference='diagnosis,nonIBD;dysbiosis_state,none;antibiotics,No'
```
* Make sure to provide the full path to the MaAsLin3 executable (i.e.
`./R/maaslin3.R`).
* In the demo command:
* ``inst/extdata/HMP2_taxonomy.tsv`` is the path to your data (or
features) file
* ``inst/extdata/HMP2_metadata.tsv`` is the path to your metadata file
* ``command_line_output`` is the path to the folder to write the output
## Tool comparison ##
| Tool | Models differential abundance | \
Models differential prevalence | Models allowed | Accounts for\
compositionality | Can use experimental absolute abundance information |
|:---:|:---:|:---:|:---:|:---:|:---:|
| MaAsLin 3 | Yes | Yes, for any model type | Any lme4 model\
(general mixed models), level-vs-level differences, group-wise differences, \
feature-specific covariates, contrast tests | Yes | Yes |
| MaAsLin 2 | Yes | Only insofar as it affects abundance | Fixed and random \
effects (subset of general mixed models) | No | No |
| ANCOM-BC2 | Yes | Only when completely absent in a group | Any lme4 model \
(general mixed models), level-vs-level differences, group-wise differences, \
contrast tests | Yes | No |
| ALDEx2 | Yes | Only insofar as it affects abundance | Fixed \
effects | Yes | No |
| LinDA | Yes | Only insofar as it affects abundance | Any lme4 model \
(general mixed models), group-wise differences, contrast tests | Yes | No |
```{R}
sessionInfo()
```
```{R}
# Clean-up
unlink('hmp2_output', recursive = TRUE)
unlink('spike_in_demo', recursive = TRUE)
unlink('total_scaling_demo', recursive = TRUE)
unlink('random_effects_output', recursive = TRUE)
unlink('interaction_output', recursive = TRUE)
unlink('ordered_outputs', recursive = TRUE)
unlink('group_outputs', recursive = TRUE)
unlink('contrast_test_output', recursive = TRUE)
```