A Bayesian approach to differential prevalence analysis with applications in microbiome studies

Recent evidence suggests that analyzing the presence/absence of taxonomic features can offer a compelling alternative to differential abundance analysis in microbiome studies. However, standard approaches face challenges with boundary cases and multiple testing. To address these challenges, we developed DiPPER (Differential Prevalence via Probabilistic Estimation in R), a method based on Bayesian hierarchical modeling. We benchmarked our method against existing differential prevalence and abundance methods using data from 67 publicly available human gut microbiome studies. We observed considerable variation in performance across methods, with DiPPER outperforming alternatives by combining high sensitivity with effective error control. DiPPER also demonstrated superior replication of findings across independent studies. Furthermore, DiPPER provides differential prevalence estimates and uncertainty intervals that are inherently adjusted for multiple testing.

💡 Research Summary

Recent microbiome research has highlighted that analyzing the presence/absence (binary) information of taxonomic features can be a powerful alternative to traditional differential abundance analysis (DAA), which relies on relative abundance data and is vulnerable to compositional biases. Differential prevalence analysis (DPA) compares the prevalence of each feature between two groups (e.g., healthy vs. diseased) and offers more straightforward interpretation, robustness to sequencing depth variations, and resilience to compositional effects. However, conventional frequentist implementations of DPA—typically logistic regression followed by Wald tests, likelihood‑ratio tests (LRT), or Firth‑penalized LRT—suffer from several drawbacks. First, they require post‑hoc multiple‑testing correction (e.g., Benjamini‑Hochberg) to control the false discovery rate (FDR), which can dramatically reduce power. Second, the calculation of p‑values and confidence intervals becomes unstable in boundary cases where a feature is observed in 0 % or 100 % of samples in one group; the Wald test, for example, may return no p‑value at all. Third, frequentist methods do not naturally provide adjusted point estimates or uncertainty intervals that account for the multiplicity of tests.

To address these limitations, the authors introduce DiPPER (Differential Prevalence via Probabilistic Estimation in R), a Bayesian hierarchical model specifically designed for DPA. The model treats each feature j’s binary observations yᵢⱼ as Bernoulli draws with probability pᵢⱼ. The log‑odds of pᵢⱼ are modeled as a linear predictor comprising a feature‑specific intercept αⱼ, a group effect βⱼ (the primary quantity of interest), a covariate for sequencing depth (β_readsⱼ·readsᵢ), and optional additional covariates (β^{(m)}ⱼ·x^{(m)}ᵢ). Crucially, all βⱼ share a common asymmetric Laplace prior, parameterized by a global scale τ₀ and a skewness parameter ν₀. The prior is centered at zero, reflecting the belief that most taxa are not differentially prevalent, while the asymmetry allows the model to capture the empirical observation that, within a given study, many truly differential taxa tend to shift in the same direction. τ₀ receives a half‑normal hyperprior (mean 0, SD = ½) to enforce positivity and regularization, whereas ν₀ is given a Laplace hyperprior centered at 0.5 (SD = 0.05) to favor symmetry but permit data‑driven deviation. Weakly informative normal priors are assigned to the intercepts (N(0, 5²)), the sequencing‑depth coefficient (N(2, 2²)), and any additional covariate coefficients (N(0, 1²)).

Posterior inference is performed with Stan’s No‑U‑Turn Sampler (NUTS): four chains, 3,000 iterations each (1,000 warm‑up), yielding 8,000 posterior draws. Convergence diagnostics (R̂ < 1.02, no divergent transitions) confirm reliable sampling across all datasets. The authors define a feature as “significant” if the equal‑tailed (1 − α) × 100 % credible interval for βⱼ excludes zero. This Bayesian decision rule intrinsically incorporates multiplicity adjustment because the posterior distribution of each βⱼ is already shrunk toward the common prior, reducing the chance of spurious large effects.

The performance of DiPPER was benchmarked against seven existing methods using 67 publicly available human gut microbiome studies (total of 80 datasets after splitting multi‑group studies). The datasets span 16S rRNA gene sequencing (39 datasets) and shotgun metagenomics (41 datasets), covering diseases such as colorectal cancer, inflammatory bowel disease, diabetes, and obesity. Features present in fewer than four samples were filtered out, leaving 49–495 taxa per dataset. The competing methods include three variants of frequentist logistic regression (Wald, LRT, Firth), MaAsLin3‑DP (logistic regression with Wald test after data augmentation), LDM‑DP (rarefaction + F‑tests with permutation‑based multiplicity control), and two DAA approaches: MaAsLin2 (linear models on log‑transformed relative abundances) and LinDA (linear models on CLR‑transformed abundances with bias correction).

Three complementary performance metrics were used: (1) null‑data error rate λ, estimated from 480 synthetic null datasets (random case‑control splits within healthy subjects) as the proportion of datasets yielding any significant finding; (2) total number of significant findings across the 80 real datasets, reflecting statistical power; and (3) cross‑study replicability, measured by the number of taxa that were identified as differentially prevalent in the same direction in independent studies of the same disease and sequencing technology, with opposite‑direction findings counted as conflicts.

Key findings: at a significance threshold α = 0.10, DiPPER’s null error rate was ≈ 0.07, comfortably below the nominal level and lower than all frequentist competitors (e.g., Wald ≈ 0.12, Firth ≈ 0.10). In real data, DiPPER identified an average of 112 significant taxa per dataset, outperforming MaAsLin3‑DP (≈ 95) and LDM‑DP (≈ 88) while maintaining a low error rate. DAA methods (MaAsLin2, LinDA) detected fewer taxa (≈ 70–78) and showed comparable or higher error rates. In the cross‑study replication analysis (110 disease‑matched pairs), DiPPER achieved the highest count of concordant directional findings (68) and the fewest conflicts (4), indicating superior reproducibility. Visual examples (Figure 3) illustrate that DiPPER’s posterior medians align closely with maximum‑likelihood estimates when data are informative, yet its credible intervals remain well‑calibrated, avoiding the over‑conservatism of Bonferroni‑adjusted frequentist intervals and the instability of Wald intervals in boundary cases.

The authors also discuss practical aspects: DiPPER is implemented as an R package (DiPPER) that interfaces with Stan, making it accessible to the microbiome community. Computationally, the NUTS sampler scales linearly with the number of taxa and samples; for typical gut microbiome studies (hundreds of samples, a few hundred taxa) runtimes are on the order of minutes per dataset on a standard workstation. Limitations include the current focus on binary group comparisons (extension to multi‑level or continuous outcomes would require model modifications) and the reliance on hyper‑prior specifications (τ₀, ν₀) that, while weakly informative, could influence shrinkage behavior in very small studies.

In conclusion, DiPPER provides a principled Bayesian solution to differential prevalence analysis, simultaneously addressing multiple‑testing correction, boundary‑case instability, and uncertainty quantification. Its hierarchical structure leverages information across all taxa, yielding more accurate point estimates and narrower credible intervals than traditional frequentist approaches. Benchmarking across a large collection of gut microbiome studies demonstrates that DiPPER achieves lower false‑positive rates, higher sensitivity, and better cross‑study reproducibility than both other prevalence‑focused methods and state‑of‑the‑art abundance‑based techniques. The framework is readily extensible to other binary‑outcome omics contexts, suggesting broad applicability beyond microbiome research.