A Bayesian Functional Concurrent Zero-Inflated Dirichlet-Multinomial Regression Model with Application to Infant Microbiome

1. Problem Statement

The paper addresses significant analytical challenges in modeling longitudinal microbiome data, specifically within the context of infant development. Key challenges include:

• Compositionality: Microbiome sequencing yields relative abundances (counts summing to a fixed total), requiring models that respect the simplex constraint. Ignoring this leads to spurious associations.
• Zero-Inflation: Data contains excess zeros arising from two sources: structural zeros (the microbe is truly absent) and at-risk zeros (the microbe is present but not detected).
• Overdispersion and High Dimensionality: Data exhibits high variance and often involves hundreds of taxa.
• Longitudinal Dynamics: Existing models often assume constant effects over time or fail to handle repeated measures. They cannot capture time-varying (functional) effects of covariates on microbial composition.
• Limitations of Current Methods: Standard Generalized Linear Models (GLMs) often ignore compositionality. Dirichlet-Multinomial (DM) models handle compositionality but lack zero-inflation handling or repeated measures capabilities. Existing zero-inflated DM models do not account for longitudinal repeated measures or functional covariate effects.

2. Methodology: The FunC-ZIDM Model

The authors propose a Functional Concurrent Zero-Inflated Dirichlet-Multinomial (FunC-ZIDM) regression model. This is a Bayesian hierarchical model designed to handle all the aforementioned challenges simultaneously.

A. Statistical Framework

• Likelihood: The observed counts $z_i(t)$ for individual $i$ at time $t$ follow a Multinomial distribution with total counts $\dot{z}_i(t)$ and relative abundance vector $\psi_i(t)$.
• Zero-Inflation Mechanism: The relative abundances $\psi_i(t)$ are derived from a reparameterized Dirichlet distribution using zero-inflated Gamma variables.
  • A latent binary indicator $\eta_{ij}$ determines if taxon $j$ is "at-risk" of being zero.
  • If $\eta_{ij}=0$, the count is a structural zero. If $\eta_{ij}=1$, the count follows a Gamma distribution (potentially resulting in an at-risk zero if the draw is zero).
• Functional Regression: The concentration parameters $\gamma_{ij}(t)$ of the Gamma distribution are modeled via a log-linear regression with functional coefficients:
  $$\log[\gamma_{ij}(t)] = \beta_{j0}(t) + \sum_{p=1}^P \beta_{jp}(t)x_{ip}(t) + r_{ij}$$
  • $\beta_{jp}(t)$: Smooth functions of time representing the effect of covariate $p$ on taxon $j$.
  • $r_{ij}$: Individual- and taxon-specific random intercepts to account for repeated measures and between-subject correlation.
• Basis Expansion: The functional coefficients $\beta_{jp}(t)$ are modeled using B-spline basis expansions.
• Regularization: To handle high dimensionality and induce sparsity (shrinking non-active coefficients to zero), the authors employ a Regularized Horseshoe prior on the spline coefficients. This prevents overfitting and handles weak identifiability in high-dimensional settings.

B. Inference and Computation

• Posterior Sampling: Inference is performed using a Metropolis-Hastings within Gibbs algorithm.
• Efficiency: Auxiliary parameterizations (Makalic and Schmidt, 2015) are used for the horseshoe prior, and $\theta_j$ (at-risk probability) is integrated out to sample $\eta_{ij}$ efficiently.
• Outputs: The model provides posterior distributions for:
  • Time-varying relative abundances ($RA_j$).
  • Multiplicative differences in relative abundance given covariate changes.
  • $\alpha$-diversity: Calculated using Hill's diversity metric (specifically $l \in [0, 1)$ to avoid triviality), allowing for holistic assessment of diversity changes over time.

3. Key Contributions

1. Novel Model Architecture: First model to simultaneously handle compositionality, zero-inflation, repeated measures, and functional (time-varying) covariate effects in a unified Bayesian framework.
2. Scalability: The model is designed to scale to large compositional spaces (hundreds of taxa) via regularization.
3. Software Implementation: The authors provide the FunCZIDM R package and a Shiny app to facilitate implementation, simulation, and visualization of results.
4. Theoretical Rigor: The use of Regularized Horseshoe priors ensures robust shrinkage of non-significant functional effects, addressing the "curse of dimensionality" common in microbiome studies.

4. Results

A. Simulation Study

• Performance: The FunC-ZIDM model significantly outperformed alternatives (standard DM, Zero-Inflated DM without random effects, and functional DM without zero-inflation) in terms of:
  • Interval Coverage: Achieved nominal ~95% coverage for credible intervals, whereas competitors dropped to 60–87% due to bias from unmodeled zero-inflation.
  • Accuracy: Low Average Root Mean Squared Error (ARMSE) for active coefficients.
  • Robustness: Performance remained stable as the number of taxa increased (tested up to $J=1000$) and under varying levels of zero-inflation.
• Comparison: Ignoring zero-inflation led to biased estimates toward zero and constant functions. Ignoring compositionality led to spurious associations.

B. Application: Preterm Infant Microbiome

The model was applied to data from 58 preterm infants (922 samples, 16 taxa) to investigate associations with gestational age and breast milk percentage.

• Baseline Dynamics: Identified distinct temporal trajectories for dominant taxa (Bacilli, Clostridia, Gammaproteobacteria). For example, Clostridia abundance increased after day 30, while Gammaproteobacteria peaked early and declined.
• Gestational Age:
  • Higher gestational age was associated with higher Clostridia abundance and lower Gammaproteobacteria abundance.
  • $\alpha$-diversity was positively associated with gestational age (more mature infants had higher diversity).
• Breast Milk:
  • Higher percentages of breast milk (>50%) were positively associated with $\alpha$-diversity.
  • The model revealed non-constant effects: The impact of breast milk on Clostridia and Bacilli varied significantly over the first 50 days of life, a nuance missed by constant-effect models.
• Comparison to Previous Studies: Unlike the original analysis (La Rosa et al., 2014) which used mixed models, FunC-ZIDM provided time-varying insights and corrected for compositional constraints, revealing that the association between breast milk and Gammaproteobacteria was not significant when properly modeled, while Clostridia showed a positive association.

5. Significance

• Clinical Relevance: The model allows clinicians to identify specific time windows where dietary interventions (e.g., increasing breast milk) are most effective for specific microbial taxa or diversity metrics.
• Reproducibility: By capturing functional effects rather than averaging them over time, the model addresses a major source of irreproducibility in microbiome research.
• Methodological Advancement: It sets a new standard for longitudinal microbiome analysis by rigorously addressing the "four horsemen" of microbiome data challenges: compositionality, zero-inflation, overdispersion, and temporal dynamics.
• Flexibility: The framework can be extended to other continuous covariates (e.g., modeling how the effect of an exposure changes as a function of gestational age rather than time).

In conclusion, the FunC-ZIDM model provides a robust, scalable, and interpretable tool for understanding the dynamic development of the infant microbiome, offering superior statistical performance over existing methods and actionable insights for pediatric health.