---
title: "Differential prevalence analysis with DiPPER"
author:
  - name: Juho Pelto
    affiliation: University of Turku
    email: jepelt@utu.fi
  - name: Leo Lahti
    affiliation: University of Turku
    email: leo.lahti@iki.fi
package: DiPPER
output:
  BiocStyle::html_document:
    toc_float: true
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Differential Prevalence Analysis with DiPPER}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# Checks whether CmdStan is available and the model-fitting chunks can be
# evaluated only when CmdStan is present. Otherwise the pre-computed fits in
# inst/extdata are used.
cmdstan_ok <- !identical(Sys.getenv("DIPPER_NO_CMDSTAN"), "true") &&
    requireNamespace("instantiate", quietly = TRUE) &&
    instantiate::stan_cmdstan_exists()

cache_file <- function(x) {
    system.file("extdata", x, package = "DiPPER")
}

have_cache <- nzchar(cache_file("fit_hintikka.rds")) &&
    nzchar(cache_file("fit_vatanen.rds"))

# Run the DiPPER models?
run_fit <- cmdstan_ok

# Can results be shown?
show_res <- cmdstan_ok || have_cache
```

# Introduction

Recent research suggests that presence/absence information of taxonomic
features (e.g., species or genera) alone, rather than complete abundance
profiles, can be highly informative when analyzing associations between
microbiomes and external variables such as health outcomes
[@Giliberti2022HostTaxab; @Karwowska2025EffectsData]. In particular,
*differential prevalence analysis (DPA)*, which in its simplest form compares
feature prevalences between two groups (e.g., healthy vs. non-healthy
subjects), represents a compelling complementary or even alternative approach
to traditional differential abundance analysis
[@Pelto2025ElementaryAnalysis; @Nickols2026MaAsLinDiscovery].

DiPPER (**Di**fferential **P**revalence via **P**robabilistic **E**stimation
in **R**) is a Bayesian hierarchical model-based approach designed for
differential prevalence analysis. In addition to the main variable of interest
(e.g., host disease status), it can incorporate other *covariates* (e.g., age,
sex, BMI). It can also be used to analyze *longitudinal* (or repeated measures)
microbiome data.

## Advantages of DiPPER

While differential prevalence analysis can be performed using approaches based
on standard (frequentist) logistic regression, such approaches suffer certain
limitations that DiPPER aims to solve:

* **Robustness in boundary cases:** Boundary cases (or quasi-complete
  separation), where a taxon is completely absent or fully present in one
  study group, are common in microbiome data. Standard approaches may have
  difficulties in these cases. For instance, the default logistic regression
  implementation in R (based on the Wald test) yields non-sensible *p*-values
  and infinite parameter estimates and confidence intervals in boundary cases.
  DiPPER, however, remains robust in these cases, providing finite estimates
  and uncertainty intervals.

* **Applicability to longitudinal data:** The aforementioned boundary cases
  can be even more problematic in longitudinal data modeling. While there are
  robust versions of frequentist logistic regression (e.g., Firth's penalized
  regression [@Firth1993BiasEstimates]), they can be difficult to generalize
  to longitudinal settings. DiPPER, however, provides a single, unified
  robust framework for both cross-sectional and longitudinal analyses.

* **Inherent multiplicity adjustment:** Standard approaches typically calculate
  an unadjusted *p*-value for each feature and subsequently apply a post-hoc
  adjustment (e.g., the Benjamini-Hochberg procedure) to control for multiple
  testing (while leaving the corresponding regression coefficients and
  confidence intervals unadjusted). In contrast, multiplicity adjustment in
  DiPPER occurs naturally via hierarchical shrinkage. The differential
  prevalence parameters are assumed to originate from a common (Asymmetric)
  Laplace prior distribution. This shrinks the estimates toward zero (and toward
  each other), providing effectively multiplicity-adjusted estimates and
  intervals that can directly be used to assign "statistical significance".

## Benchmarking results against alternative methods

We have benchmarked DiPPER against alternative methods using real data from
57 human gut microbiome studies. The performance was compared against three
versions of frequentist logistic regression, DPA implementations in the
MaAsLin3 [@Nickols2026MaAsLinDiscovery] and LDM [@Hu2021AMicrobiome] packages,
and two differential abundance methods (MaAsLin2
[@Mallick2021MultivariableStudies] and LinDA [@Zhou2022LinDA:Data]).

DiPPER demonstrated high sensitivity in detecting potentially differentially
prevalent features while maintaining a well-calibrated family-wise error rate
under the global null hypothesis. Moreover, it outperformed all
alternatives in the replication of findings across independent studies.

## More information

A detailed description of DiPPER, along with comprehensive benchmarking
results, can be found in our pre-print [@Pelto2026DiPPER:Studies].

# Installation

Installing DiPPER is generally straightforward, but it requires dependencies
from both Bioconductor (for handling microbiome data structures) and Stan (for
Bayesian inference).

## Stan backend (cmdstanr)

Unless you already have `cmdstanr` installed, you will need to set it up first.
While configuring the Stan backend can occasionally require some
troubleshooting depending on your operating system, following these steps
should generally ensure a smooth setup.

First, install the `cmdstanr` R package from the Stan R-universe repository.

```{r install_stan, eval=FALSE}
# Install cmdstanr package
install.packages("cmdstanr",
                 repos = c("https://stan-dev.r-universe.dev",
                           getOption("repos")))

# Check C++ toolchain and install the backend
cmdstanr::check_cmdstan_toolchain(fix = TRUE)
cmdstanr::install_cmdstan()
```

## Bioconductor dependencies and DiPPER

Once CmdStan is ready, you can install DiPPER.

*Note: DiPPER compiles its Stan models during installation, which may take a
few minutes and produce verbose output. This is expected, and it is done only
once: afterwards the models are ready to use.*

```{r install_bioc, eval=FALSE}
# Install BiocManager if it is not already installed
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("DiPPER")
```


# Quick start

Below is a minimal workflow using the built-in `tse_hintikka` dataset, which
compares the gut microbiomes of rats (N = 20 + 20) on a high-fat versus low-fat
diet (`Fat`) while controlling for xylo-oligosaccharide (`XOS`)
supplementation. The dataset originates from
[Hintikka et al. (2021)](https://www.mdpi.com/1660-4601/18/8/4049).

To run DiPPER, we need to provide

- TSE object (`tse = tse_hintikka`)
- Name of the assay (`assay.type = "counts"`)
- Type of data in the assay. As our data consist of sequencing counts we set
  `data.type = "counts"`.
- Model formula: `formula = ~ Fat + XOS`. By default, the first term (here
  `Fat`) is the *variable of interest* (another term can be chosen with
  `var.of.interest`).
- As our data consist of sequencing counts, it is recommended to control for
  sequencing depth by setting `read.depth = TRUE`.

```{r quick_start, eval=run_fit}
library(DiPPER)

# Load dataset (TreeSummarizedExperiment object)
data("tse_hintikka")

# Run DiPPER
fit <- dipper(
    tse = tse_hintikka,
    assay.type = "counts",
    data.type = "counts",
    formula = ~ Fat + XOS,
    read.depth = TRUE
)
```

```{r quick_start_cached, eval=!run_fit && have_cache, echo=FALSE}
# CmdStan is not available on this machine: use the pre-computed fit instead
library(DiPPER)
fit <- readRDS(cache_file("fit_hintikka.rds"))
```

We first print the details of the DiPPER run:

```{r quick_print, eval=show_res}
print(fit)
```

Next, we extract the summarized results as a `data.frame` and print the first
seven rows (the most "significant" genera). Note that the data frame contains
the results only for the variable of interest, which in this case is the `Fat`
variable indicating the diet, low-fat ('Low', the reference level) or high-fat
('High').

```{r quick_results, eval=show_res}
# Extract the results
res <- summary(fit)

# Print the first seven rows
res[1:7, ]
```

The summarized results `data.frame` contains the following columns:

* **`taxon`**: Taxon name (here: genus).
* **`log_or`**: Differential prevalence estimate *for the variable of interest*
  (here: `Fat`), expressed as a log(odds ratio).
* **`lwr`**, **`upr`**: Lower and upper limits of the 95% credible interval for
  the differential prevalence.
* **`significant`**: A logical value indicating whether the 95% credible
  interval excludes zero. This roughly corresponds to traditional "statistical
  significance" based on multiplicity-adjusted *p*-values.
* **`pseudo_q`**: Indicates the highest "significance level" at which the
  result would be considered "significant" (i.e., the credible interval
  excluding zero). It thus bears some resemblance to a multiplicity-adjusted
  two-sided *p*-value (*q*-value), but it does not have any strict statistical
  interpretation.

Lastly, create a forest plot (with 95% credible intervals) showing the
significant genera (determined by the 95% credible intervals). Here, the odds
ratio (OR) values > 1 indicate higher prevalence in the high-fat diet group.

```{r quick_plot, eval=show_res, out.width="100%"}
plot(fit)
```

# More complex example

This example demonstrates a more complex case, where we have longitudinal data.

We use a built-in subset of
the [Vatanen et al. (2016)](https://pmc.ncbi.nlm.nih.gov/articles/PMC4950857/)
dataset. The dataset includes gut microbiome samples from 79 infants at one
or two time points (1.5 and/or 7 months of age). We want to study how
microbe prevalences depend on age point (1.5 months vs 7 months)
while controlling for gender and antibiotics use.

```{r load_vatanen}
library(DiPPER)

# Load dataset (TreeSummarizedExperiment object, agglomerated to genus level)
data("VatanenT_2016_subset")
```

The metadata includes the variables `subject_id` (indicating the infant),
`age_point` (1.5 months or 7 months), `antibiotics` (use of antibiotics:
yes/no), and `gender` (male/female). Note that there are 1-2 measurements per
infant in the data.

```{r vatanen_metadata}
VatanenT_2016_subset |>
    SummarizedExperiment::colData() |>
    head()
```

Note that the abundance data are now proportions (relative abundances):

```{r vatanen_abundance}
vatanen_assay <- SummarizedExperiment::assay(VatanenT_2016_subset,
                                             "relative_abundance")
vatanen_assay[c("Bifidobacterium", "Escherichia", "Veillonella"), 1:5]
```

As we want to study the effect of age point (1.5 months vs 7 months) on the
microbe prevalences, we run DiPPER by setting `age_point` as the first variable
in the formula, making it the variable of interest. As we control for
`antibiotics` use and `gender`, we put them as the second and third variables.
Lastly, as some infants have two measurements, we must
account for the longitudinality by specifying `subject_id` (indicating infant)
as a random intercept term: `(1 | subject_id)`. Our formula is thus
`~ age_point + antibiotics + gender + (1 | subject_id)`.

Furthermore, because the abundance data are proportions (and not counts), we
must set `data.type = "relabundance"`. Moreover, they carry no information about
sequencing depth, so we set `read.depth = FALSE`.

```{r assay_vatanen, eval=run_fit}
fit_long <- dipper(
    tse = VatanenT_2016_subset,
    assay.type = "relative_abundance", # Name of the assay
    data.type = "relabundance",
    formula = ~ age_point + antibiotics + gender + (1 | subject_id),
    read.depth = FALSE
)
```

```{r assay_vatanen_cached, eval=!run_fit && have_cache, echo=FALSE}
# CmdStan is not available on this machine: use the pre-computed fit instead
fit_long <- readRDS(cache_file("fit_vatanen.rds"))
```

Finally, we extract the results (for `age_point`) and create a forest plot
showing only the significant taxa. As expected, the prevalence of several taxa
increases with infant age, although the prevalence of *Staphylococcus*
decreases.

```{r plot_vatanen, eval=show_res, out.width="100%"}
# Print the details of the run
print(fit_long)

# Extract the results for age_point
res_long <- summary(fit_long)
res_long[1:7, ]

# Create a forest plot for significant taxa
plot(fit_long)
```

# Input data and other settings

## Type of abundance data

The type of input abundance data, i.e., `data.type`, must be provided. This is
to enable `dipper()` to derive presence/absence information, and optionally, to
calculate sequencing depth correctly. The following values are accepted:

- **`"counts"`**: sequencing read counts (non-negative integers). The most
  typical format for 16S and shotgun data. It is preferable to use this option
  if you have read count data available, as these data can be used to
  automatically compute and control for the sequencing depth (read depth) in
  the analysis.
- **`"relabundance"`**: relative abundances, that is, proportions between 0
  and 1 whose sample sums are at most 1. Tools such as MetaPhlAn report these
  directly.
- **`"abundance"`**: any other non-negative abundances, for example, CPM or TPM
  values.
- **`"pa"`**: presence/absence data expressed as 0s and 1s. Use this whenever
  you want to define presence and absence values yourself.

Transformed values, such as CLR- or log-transformed assays, are **not**
supported, as DiPPER cannot derive presence/absence information from them.


## Controlling for sequencing (read) depth

Whether a feature/taxon is observed in a sample depends partly on how deeply
the sample was sequenced. Therefore, it is advisable to control for the
sequencing/read depth in the analysis whenever possible. DiPPER can compute the
read depth itself when sequencing count data are provided. The user can also
provide the read depths as a metadata column.

Whether read depth is controlled for is determined by the `read.depth`
argument. For sequencing count data (`data.type = counts`) there is
no default, because DiPPER cannot detect how the data have been
preprocessed.

- **`read.depth = TRUE`** computes read depth as the sample totals from the
  provided sequencing count data. Use this with non-rarefied sequencing counts.
  (Requires `data.type = "counts"`.)
- **`read.depth = "column_name"`** takes the read depths from a metadata
  column. Use this whenever you have computed the read depths yourself.
- **`read.depth = FALSE`** leaves read depth out of the model. Use this
  with rarefied data, or when sequencing depth infromation is not available.
  (The default for every `data.type` other than `"counts"`.)

When `read.depth = TRUE` or a metadata column is provided, `dipper()` applies a
log10 transformation to the read depths and adds the variable to the model.
The variable can be seen in the `print()` output as the `log10_read_depth` term
in the model formula.


## Symmetric or asymmetric hierarchical prior

By default DiPPER uses a hierarchical Asymmetric Laplace prior for the
differential prevalence parameters of interest (`symmetric = FALSE`). This
means that the prior is allowed to be asymmetric and the degree of
asymmetry is estimated from the data. Optionally, the user can set
`symmetric = TRUE` when the prior is fixed to be symmetric around zero.

In most cases, the symmetric and asymmetric option yield rather similar results
but the symmetric option generally provides slightly more conservative results.

In the case of strong asymmetry the prior pulls the differential prevalence
estimates not only towards zero, but also towards each other. Consequently, in
such cases, DiPPER induces a small bias to the estimates (e.g. all estimates
may have the same sign). While this bias may seem suspicious, and it may yield
a few additional false positives, it should also lead to clearly higher power
to detect true effects.

It is generally recommended to first run DiPPER with the default setting
`symmetric = FALSE`. The estimated values of the hyperparameters of the prior
can then be read from the `Hierarchical prior` line of the `print()` output,
where `tau` is the prior scale and `nu` indicates the asymmetry (`nu` is
constrained between 0 and 1, with 0.5 meaning symmetry). Now, if `nu` is
clearly away from 0.5 (e.g. below 0.30 or above 0.70) **and** `tau` is at the
same time small (e.g. < 0.10), the asymmetry of the prior may have a somewhat
strong effect on the results. **If** in such a case the user suspects that the
results are driven by some technical bias in the data, they may consider using
the symmetric version.

# Session Information

```{r session_info}
sessionInfo()
```