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.

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”.

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 et al. 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.

More information

A detailed description of DiPPER, along with comprehensive benchmarking results, can be found in our pre-print (Pelto et al. 2026).

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.

# 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.

# 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).

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.

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 seconds

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’).

# 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.001000

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.

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) 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     female

Note 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.0000000

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.

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)

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

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> 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               Seqinfo_1.3.2                  
#> [29] codetools_0.2-20                htmltools_0.5.9                
#> [31] sys_3.4.3                       buildtools_1.0.0               
#> [33] sass_0.4.10                     lazyeval_0.2.3                 
#> [35] yaml_2.3.12                     tidyr_1.3.2                    
#> [37] pillar_1.11.1                   crayon_1.5.3                   
#> [39] jquerylib_0.1.4                 BiocParallel_1.47.0            
#> [41] SingleCellExperiment_1.35.2     cachem_1.1.0                   
#> [43] DelayedArray_0.39.8             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                     maketools_1.3.2                
#> [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                 ape_5.8-1                      
#> [59] withr_3.0.3                     rappdirs_0.3.4                 
#> [61] scales_1.4.0                    rmarkdown_2.32                 
#> [63] XVector_0.53.0                  matrixStats_1.5.0              
#> [65] otel_0.2.0                      instantiate_0.2.4              
#> [67] evaluate_1.0.5                  knitr_1.52                     
#> [69] GenomicRanges_1.65.4            IRanges_2.47.5                 
#> [71] rlang_1.3.0                     Rcpp_1.1.2                     
#> [73] glue_1.8.1                      tidytree_0.4.8                 
#> [75] BiocManager_1.30.27             BiocGenerics_0.59.12           
#> [77] jsonlite_2.0.0                  R6_2.6.1                       
#> [79] MatrixGenerics_1.25.0           fs_2.1.0

Firth, David. 1993. “Bias Reduction of Maximum Likelihood Estimates.” Biometrika 80 (1): 27. https://doi.org/10.2307/2336755.

Giliberti, Renato, Sara Cavaliere, Italia Elisa Mauriello, Danilo Ercolini, and Edoardo Pasolli. 2022. “Host phenotype classification from human microbiome data is mainly driven by the presence of microbial taxa.” PLoS Computational Biology 18 (4). https://doi.org/10.1371/JOURNAL.PCBI.1010066.

Hu, Yi Juan, Andrea Lane, and Glen A. Satten. 2021. “A rarefaction-based extension of the LDM for testing presence–absence associations in the microbiome.” Bioinformatics 37 (12): 1652. https://doi.org/10.1093/BIOINFORMATICS/BTAB012.

Karwowska, Zuzanna, Oliver Aasmets, Tõnu Esko, et al. 2025. “Effects of data transformation and model selection on feature importance in microbiome classification data.” Microbiome 13 (1): 1–14. https://doi.org/10.1186/S40168-024-01996-6/FIGURES/4.

Mallick, Himel, Ali Rahnavard, Lauren J. McIver, et al. 2021. “Multivariable association discovery in population-scale meta-omics studies.” PLOS Computational Biology 17 (11): e1009442. https://doi.org/10.1371/journal.pcbi.1009442.

Nickols, William A., Thomas Kuntz, Jiaxian Shen, et al. 2026. “MaAsLin 3: refining and extending generalized multivariable linear models for meta-omic association discovery.” Nature Methods 2026, January, 1–11. https://doi.org/10.1038/s41592-025-02923-9.

Pelto, Juho, Kari Auranen, Janne V. Kujala, and Leo Lahti. 2026. DiPPER: A Bayesian approach to differential prevalence analysis with applications in microbiome studies. February. https://arxiv.org/abs/2602.05938v2.

Pelto, Juho, Kari Auranen, Janne V Kujala, and Leo Lahti. 2025. “Elementary methods provide more replicable results in microbial differential abundance analysis.” Briefings in Bioinformatics 26 (2): 130. https://doi.org/10.1093/BIB/BBAF130.

Zhou, Huijuan, Kejun He, Jun Chen, and Xianyang Zhang. 2022. “LinDA: linear models for differential abundance analysis of microbiome compositional data.” Genome Biology 23 (1): 95. https://doi.org/10.1186/s13059-022-02655-5.

Differential prevalence analysis with DiPPER