Differential prevalence analysis with DiPPER
13 October 2026
Package
DiPPER 0.99.1
1 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 (Giliberti et al. 2022; Karwowska et al. 2025). 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 (Pelto et al. 2025; Nickols et al. 2026).
DiPPER (Differential Prevalence via Probabilistic Estimation 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.
1.1 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 (Firth 1993)), 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”.
1.2 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 (Nickols et al. 2026) and LDM (Hu, Lane, and Satten 2021) packages, and two differential abundance methods (MaAsLin2 (Mallick et al. 2021) and LinDA (Zhou et al. 2022)).
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.
1.3 More information
A detailed description of DiPPER, along with comprehensive benchmarking results, can be found in our pre-print (Pelto et al. 2026).
2 Installation
Installing DiPPER is generally straightforward, but it requires dependencies from both Bioconductor (for handling microbiome data structures) and Stan (for Bayesian inference).
2.1 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.
# 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()2.2 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.
# Install BiocManager if it is not already installed
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("DiPPER")3 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).
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 (hereFat) is the variable of interest (another term can be chosen withvar.of.interest). - As our data consist of sequencing counts, it is recommended to control for
sequencing depth by setting
read.depth = TRUE.
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
)We first print the details of the DiPPER run:
print(fit)
#> DiPPER Model Fit
#> ----------------
#> Model formula: ~Fat + XOS + log10_read_depth
#> Hierarchical prior: Asymmetric Laplace (tau 0.67 [0.54, 0.86], nu 0.50 [0.44, 0.55])
#> Features/taxa: 76 retained, 186 removed by filtering
#> Samples: 40
#> Posterior draws: 4000 (4 chains x 1000 iterations)
#> MCMC diagnostics: No divergent transitions or max treedepth hits
#> Convergence: max R-hat 1.006, min bulk ESS 2252, min tail ESS 2000
#> MCMC sampling time: 20.02 secondsNext, 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’).
# Extract the results
res <- summary(fit)
# Print the first seven rows
res[1:7, ]
#> taxon log_or lwr upr significant
#> 1 D_5__ASF356 -2.344015 -4.1605614 -0.8367085 TRUE
#> 2 D_5__Clostridium sensu stricto 1 -2.268650 -3.9415148 -0.8509062 TRUE
#> 3 D_5__Faecalibaculum -4.263514 -8.5262620 -1.9273439 TRUE
#> 4 D_5__Turicibacter -3.992722 -8.1047578 -1.7648139 TRUE
#> 5 D_5__Enterobacter 2.376496 0.9234157 4.0046487 TRUE
#> 6 D_5__Romboutsia -3.312785 -7.2800959 -1.1661719 TRUE
#> 7 D_5__Sellimonas 3.259861 0.9675345 7.2432577 TRUE
#> pseudo_q
#> 1 0.000125
#> 2 0.000125
#> 3 0.000125
#> 4 0.000125
#> 5 0.001000
#> 6 0.001000
#> 7 0.001000The 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.
plot(fit)
4 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) 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.
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.
VatanenT_2016_subset |>
SummarizedExperiment::colData() |>
head()
#> DataFrame with 6 rows and 4 columns
#> subject_id antibiotics age_point gender
#> <factor> <factor> <factor> <factor>
#> G80490 E002338 no 1.5 months female
#> G80455 E002473 no 7 months female
#> G80498 E002473 no 1.5 months female
#> G80621 E002681 yes 7 months female
#> G80522 E004781 no 1.5 months female
#> G80574 E004781 yes 7 months femaleNote that the abundance data are now proportions (relative abundances):
vatanen_assay <- SummarizedExperiment::assay(VatanenT_2016_subset,
"relative_abundance")
vatanen_assay[c("Bifidobacterium", "Escherichia", "Veillonella"), 1:5]
#> G80490 G80455 G80498 G80621 G80522
#> Bifidobacterium 0.0542750 0.1427629 0.0000000 0.0014688 0.0389856
#> Escherichia 0.0002381 0.5788488 0.7626766 0.6043643 0.0228112
#> Veillonella 0.6600033 0.1537560 0.0280686 0.1524510 0.0000000As 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.
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
)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.
# Print the details of the run
print(fit_long)
#> DiPPER Model Fit
#> ----------------
#> Model formula: ~age_point + antibiotics + gender + (1 | subject_id)
#> Hierarchical prior: Asymmetric Laplace (tau 0.46 [0.31, 0.66], nu 0.28 [0.19, 0.37])
#> Features/taxa: 58 retained, 151 removed by filtering
#> Samples: 90
#> Posterior draws: 4000 (4 chains x 1000 iterations)
#> MCMC diagnostics: No divergent transitions or max treedepth hits
#> Convergence: max R-hat 1.007, min bulk ESS 1810, min tail ESS 2234
#> MCMC sampling time: 153.82 seconds
# Extract the results for age_point
res_long <- summary(fit_long)
res_long[1:7, ]
#> taxon log_or lwr upr significant pseudo_q
#> 1 Anaerostipes 2.654603 1.2037014 4.742133 TRUE 0.000125
#> 2 Blautia 3.674634 1.9290329 6.417722 TRUE 0.000125
#> 3 Faecalibacterium 3.257079 1.5625498 5.839684 TRUE 0.000125
#> 4 Fusicatenibacter 3.625883 1.2731047 8.450943 TRUE 0.000125
#> 5 Intestinibacter 2.227003 0.8953533 4.070835 TRUE 0.000125
#> 6 Roseburia 2.908464 1.1798282 5.601842 TRUE 0.000500
#> 7 Actinomyces 2.583856 0.9678193 4.992881 TRUE 0.001000
# Create a forest plot for significant taxa
plot(fit_long)
5 Input data and other settings
5.1 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.
5.2 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 = TRUEcomputes read depth as the sample totals from the provided sequencing count data. Use this with non-rarefied sequencing counts. (Requiresdata.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 = FALSEleaves read depth out of the model. Use this with rarefied data, or when sequencing depth infromation is not available. (The default for everydata.typeother 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.
5.3 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.
6 Session Information
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] DiPPER_0.99.1 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] SummarizedExperiment_1.43.0 gtable_0.3.6
#> [3] xfun_0.61 bslib_0.12.0
#> [5] ggplot2_4.0.3 processx_3.9.0
#> [7] Biobase_2.73.2 lattice_0.23-1
#> [9] callr_3.8.0 yulab.utils_0.2.5
#> [11] vctrs_0.7.3 tools_4.6.1
#> [13] ps_1.9.3 generics_0.1.4
#> [15] stats4_4.6.1 parallel_4.6.1
#> [17] tibble_3.3.1 pkgconfig_2.0.3
#> [19] Matrix_1.7-6 RColorBrewer_1.1-3
#> [21] S7_0.2.2 S4Vectors_0.51.10
#> [23] lifecycle_1.0.5 compiler_4.6.1
#> [25] farver_2.1.2 treeio_1.37.1
#> [27] Biostrings_2.81.9 tinytex_0.61
#> [29] Seqinfo_1.3.2 codetools_0.2-20
#> [31] htmltools_0.5.9 sass_0.4.10
#> [33] lazyeval_0.2.3 yaml_2.3.12
#> [35] tidyr_1.3.2 pillar_1.11.1
#> [37] crayon_1.5.3 jquerylib_0.1.4
#> [39] BiocParallel_1.47.0 SingleCellExperiment_1.35.2
#> [41] cachem_1.1.0 DelayedArray_0.39.8
#> [43] magick_2.9.1 abind_1.4-8
#> [45] nlme_3.1-171 tidyselect_1.2.1
#> [47] digest_0.6.39 purrr_1.2.2
#> [49] dplyr_1.2.1 bookdown_0.48
#> [51] TreeSummarizedExperiment_2.21.0 fastmap_1.2.0
#> [53] grid_4.6.1 cli_3.6.6
#> [55] SparseArray_1.13.4 magrittr_2.0.5
#> [57] S4Arrays_1.13.2 dichromat_2.0-1
#> [59] ape_5.8-1 withr_3.0.3
#> [61] rappdirs_0.3.4 scales_1.4.0
#> [63] rmarkdown_2.32 XVector_0.53.0
#> [65] matrixStats_1.5.0 otel_0.2.0
#> [67] instantiate_0.2.4 evaluate_1.0.5
#> [69] knitr_1.52 GenomicRanges_1.65.4
#> [71] IRanges_2.47.5 rlang_1.3.0
#> [73] Rcpp_1.1.2 glue_1.8.1
#> [75] tidytree_0.4.8 BiocManager_1.30.27
#> [77] BiocGenerics_0.59.12 jsonlite_2.0.0
#> [79] R6_2.6.1 MatrixGenerics_1.25.0
#> [81] fs_2.1.0