A Flexible Zero-Inflated Poisson-Gamma Model with application to microbiome read counts

For taxon count $W_{ijk}$, we decompose the mean of Poisson distribution in eq (1) into a true abundance level $\lambda_{ijk}$ and a multiplicative factor $U_{ijk}$ and consider the following hierarchical Zero-Inflated Poisson-Gamma (ZIPG) model:

where $\lambda_{ijk}$ represents the true abundance level for taxon $k$ on the $j$th observation from subject $i$, and $p_{k}$ denotes the zero-inflation parameter describing the probability of true zero occurrence of taxon $k$. $U_{ijk}$ follows Gamma distribution with the same rate parameter and shape parameter $\theta_{ik}^{-1}$. This factor does not change the average abundance level given the fact that $E(U_{ijk})=1$. On the other hand, $Var(U_{ijk})=\theta_{ik}$ allows extra variation of the observed abundance $W_{ijk}$ around the average level $\lambda_{ijk}$, a phenomenon described as the deviation of observed abundance to unobserved true abundance. It is of note that the variability term, $\theta_{ik}$, remains the same for all measurements (i.e., $j=1,...,n_{i}$) across individual $i$. Thus it reflects the stability of taxon $k$ in the microbial system of individual $i$. This is motivated by the assumption about metagenomic sequencing bias being taxon-specific but not sample-specific made in literature (McLaren et al., 2019). Finally, we assume the zero-inflation parameter $p_{k}$ is only taxon-specific and is common across samples ($j$) and individuals ($i$). This is because many experimental factors in sequencing can introduce the measurement of zeros (Silverman et al., 2020) and hence it is less appealing to link zero inflation parameter $p_{k}$ to other covariates (e.g., biological or environmental conditions) possessed by individuals or samples. We have checked the sensitivity of ZIPG under model misspecification when this assumption is violated by comparing the performance of the current ZIPG model (3) to a full ZIPG model (denoted as ZIPG-full), which replaces $p_{k}$ by $p_{ijk}$ and links $p_{ijk}$ to covariates. Results in Section 4.4 indicates that the current ZIPG model tends to have better model-fitting performance than ZIPG-full, which is consistent with previous empirical conclusions that modeling zero inflation in sequence count data should be careful (with respect to underlying zero generating process) and numerical evidence tends to favor simpler models (Silverman et al., 2020).

To further explore the new ZIPG framework, let $W_{ijk}^{\rm Pois}\sim{\rm Poisson}(\lambda_{ijk})$ denote as the random variable generated from the Poisson distribution and $W_{ijk}^{\rm PG}$ as the random variable generated from the Poisson-Gamma part in eq (3). We have

Thus, the mean of the Poisson-Gamma distribution is the same as the regular Poisson distribution, but its variance is multiplied by $(1+\lambda_{ijk}\theta_{ik})$ to account for the over-dispersion caused by the multiplicative measurement error factor $U_{ijk}$. We also observe the similar phenomenon of using a more sophisticated hierarchical model to account for over-dispersion in microbiome data analysis, such as the Beta-Binomial distribution (Martin et al., 2020) and Dirichlet-multinomial distribution (La Rosa et al., 2012). In this paper, we refer to $\lambda_{ijk}$ as the abundance mean parameter and $\theta_{ik}$ as the abundance dispersion parameter.

A critical task in microbiome research is to explore the relationship between taxa abundances and covariates of interest. Compared to noisy observed counts $W_{ijk}$, it is more interesting to investigate the association between key parameters (i.e., $\lambda_{ijk}$ and $\theta_{ik}$) of the underlying taxa abundance distribution and covariates of interest. To this end, we connect the mean parameter and dispersion parameter with different sets of covariates, respectively. In the repeated measures or longitudinal study design considered in the current paper, some covariates vary across different samples within the same subject, such as dietary intake, and we refer to them as “time-dependent” covariates. Other covariates, which do not change during the study, such as the treatment group assigned at the beginning of the study, are referred to as “time-independent” covariates. Since the dispersion parameter $\theta_{ik}$ describes the deviation of short-term abundance from long-term mean abundance $\lambda_{ijk}$, we propose to link it to time-independent covariates, supported by the evidence of microbiome stability perturbation in Morgan et al. (2012) and Couch et al. (2021). For the mean abundance $\lambda_{ijk}$, it can be linked to either time-dependent covariates or time-independent covariates, or both. Therefore, we define the following link functions:

where $\boldsymbol{X_{ij}}\in\mathbb{R}^{d_{1}}$ is a vector of covariates associated with $\lambda_{ijk}$, $\boldsymbol{X^{*}_{i}}\in\mathbb{R}^{d_{2}}$ include covariates associated with $\theta_{ik}$, $\boldsymbol{\beta_{k}}=(\beta_{1},...,\beta_{d_{1}})^{T}$ and $\boldsymbol{\beta_{k}^{*}}=(\beta_{1}^{*},...,\beta_{d_{2}}^{*})^{T}$ are regression coefficients of interest. We allow overlapped covariates in two models. For ease of presentation, we term the two models in (4) as ZIPG mean model and ZIPG dispersion model, respectively. The $\log(M_{ij})$ term in the mean model accounts for the effect of sequencing depth variation on mean abundances. The same offset is used in previous ZINB models (Zhang and Yi, 2020; Robinson et al., 2010), and the log of median-of-ratios is another possible candidate for offset used in literature (Xu et al., 2021; Love et al., 2014). While most existing ZINB-based methods (2) models $\mu_{ijk}$ but treat $\alpha_{k}$ as a nuisance parameter, our approach allows additional time-independent covariates linked to the dispersion of abundance, which would lead to better model fit and more powerful association analysis as will be shown later in this paper. According to previous discussions, we suggest including time-independent covariates such as demographic and lifestyle-related variables in $\boldsymbol{X^{*}_{i}}$. All covariates of interest, regardless of time-dependent or not, shall be included in $\boldsymbol{X_{ij}}$. By testing the coefficients $\boldsymbol{\beta_{k}}$ and $\boldsymbol{\beta^{*}_{k}}$, we can detect differential abundance and differential variability impacted by physiological status or host environment, respectively. Finally, if we do not include any covariates in $\boldsymbol{X^{*}_{i}}$, our ZIPG model will degenerate to ZINB (i.e., model (2)) with $\theta_{ik}=\alpha_{k}$ for any subject. Throughout this paper we choose a logarithmic link function $g(x)=g^{*}(x):=\log(x)$ to ensure $\lambda_{ijk}>0$ and $\theta_{ik}>0$.

The proposed hierarchical ZIPG model has several key advantages. First, to handle over-dispersion in microbiome count data, we decompose abundances into the long-term true abundance and its individual-level variation through a multiplicative factor. Second, we allow different sets of variables associated with the mean and the variation of abundance and provide explanations for variations in the individual-level microbial system. Thus, the proposed method is not only able to test the change of the mean abundance but also the microbiome stability affected by physiological status or host environment, which cannot be achieved by existing ZIP models. In addition, parameter $\theta_{ik}$ also controls the skewness and kurtosis of the Gamma distribution. That is, we allow the higher-order moments (or the shape of the distribution) to be linked to covariates, which is another feature that is typically missed in existing models.

Given covariates $\boldsymbol{X}$ and $\boldsymbol{X^{*}}$, observed count data $\boldsymbol{W}$ and sequencing depth $\boldsymbol{M}$, we write the log-likelihood of $\boldsymbol{\Omega}$ as follows:

where $\lambda_{ij}$ and $\theta_{i}$ are functions of $\boldsymbol{\Omega}$ defined in eq (4) and $\Gamma(x)=\int_{0}^{+\infty}t^{x-1}\exp(-t)dt$ is the Gamma function. The log-likelihood eq (5) is non-concave in $\boldsymbol{\Omega}$ (see Section 1.1 of Supplements). In practice, we found that directly maximizing eq (5) can cause trouble in distinguishing zeros from the Poisson-Gamma part and the other zero-inflation part of the ZIPG model, leading to an unreasonably low estimator of $\gamma$. A similar phenomenon has been observed for pscl, with more discussions provided in Section 4. Therefore, we use the EM algorithm for a reliable estimator of $\boldsymbol{\Omega}$.

Let $z_{ij}$ be the latent variable, where $z_{ij}=1$ indicates $W_{ij}$ is generated from zero-inflated part with probability $p=\exp(\gamma)/(\exp(\gamma)+1)$, and $z_{ij}=0$ indicates $W_{ij}$ is generated from the Poisson-Gamma distribution with probability $1-p$. The log-likelihood with complete data $\{W_{ij},M_{ij},X_{ij},X_{ij}^{*},z_{ij}\}$ for $i=1,...,n$ and $j=1,...,n_{i}$ is written as:

The detailed procedure of the EM algorithm is provided in Algorithm 1. We first initialize $\boldsymbol{\Omega}^{(0)}$ by the results of zeroinfl in pscl with $\boldsymbol{\beta^{*}}=\boldsymbol{0}$. $p_{ij}^{(0)}$ is adjusted to the proportion of observed zeros of $\boldsymbol{W}$ to avoid the local maximum at the start point. Then we can repeat the E-step and M-step until convergence or the maximum number of iterations $t_{\rm max}$ is reached. For the $t$-th iteration, in M-step, we update $\boldsymbol{\Omega}$ by maximizing eq (6) using BFGS in R function optim (Broyden (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970)) given latent variable $\boldsymbol{z}^{(t-1)}$. The gradient of eq (6) is applied to improve the computational efficiency (see Section 1.2 of Supplements).

In E-step, we update latent variables $z_{ij}$ by their conditional expectations given $\boldsymbol{\Omega}^{(t)}$ estimated from M-step:

Substituting latent variables $z_{ij}^{(t)}=z(W_{ij},\boldsymbol{\Omega}^{(t)})$ into eq (6), we can write the expectation of log-likelihood on a single observation as a new function

then $Q(\boldsymbol{\Omega}\mid\boldsymbol{\Omega}^{(t)})=\sum_{ij}Q_{ij}(\boldsymbol{\Omega}\mid\boldsymbol{\Omega}^{(t)})$ is the quantity maximized in M-step equivalently.

According to Theorem 6 in Wu (1983), since $\partial Q(\boldsymbol{\Omega}\mid\boldsymbol{\Omega}^{*})/\partial\boldsymbol{\Omega}$ is continuous in $\boldsymbol{\Omega}$ and $\boldsymbol{\Omega}^{*}$, and our algorithm ensures $\partial Q(\boldsymbol{\Omega}\mid\boldsymbol{\Omega}^{(t)})/\partial\boldsymbol{\Omega}\mid_{\boldsymbol{\Omega}=\boldsymbol{\Omega}^{(t+1)}}=\textbf{0}$ in each iteration, then EM estimator will converge to a stationary point of $L(\boldsymbol{\Omega}\mid\boldsymbol{W})$. Under mild conditions, EM Algorithm converges to the local maximum (Wu, 1983, Theorem 3). Based on our suggested initialization, namely the MLE assuming $\boldsymbol{\beta^{*}}=0$ with adjusted $p^{(0)}$, numerical studies suggest that our estimators are nearly unbiased.

In this paper, proving the asymptotic normality of the EM estimator is less of our interest, and despite its lack of rigorous theoretical justification, the EM estimator has been essentially treated as MLE in literature. Thus, we treat EM estimator $\boldsymbol{\hat{\Omega}}$ as MLE and construct the Wald test statistics and confidence interval based on the asymptotic normality of MLE. Further, we carefully consider the potential practical issues and evaluate different bootstrap methods and interval construction strategies with extensive numerical studies.

Consider the null hypothesis $H_{0}:\boldsymbol{A\boldsymbol{\Omega}=b}$ as the general form to test arbitrary subsets and linear combinations of parameters within $\boldsymbol{\Omega}$, where $\boldsymbol{A}\in\mathbb{R}^{r\times(d_{1}+d_{2}+3)}$ has full rank, $r<d_{1}+d_{2}+3$, and $\boldsymbol{b}\in\mathbb{R}^{r}$. The classical Wald test statistic can be constructed as

where $\boldsymbol{V}$ is the asymptotic covariance matrix of $\boldsymbol{\Omega}$, $N$ is the sample size, $T_{\rm Wald}$ is asymptotic $\chi^{2}_{r}$ under the null hypothesis. In practice, we found that directly using the inverse of observed information derived from the last M-step might underestimate the variance matrix (see Section 3.1 of the Supplements). Thus, we propose to use nonparametric bootstrap to estimate the covariance matrix and construct the test statistic as in Algorithm 2.

The 100$(1-\alpha)\%$ confidence interval for any single parameter $\beta\in\boldsymbol{\Omega}$ can also be obtained through nonparametric bootstrap Wald test as $\left(\hat{\beta}-z_{\alpha/2}{\rm SD}(\beta^{b}),\hat{\beta}+z_{\alpha/2}{\rm SD}(\beta^{b})\right)$ based on the standard error (SD) of $\beta^{b}\in\boldsymbol{\Omega^{b}}$ from $B$ bootstrap samples. When the sample size is small or when the collected data is unbalanced (e.g., Romero data in Section 5), parametric bootstrap could be applied for more robust results as in Martin et al. (2020). Thus, we also developed the Wald test based on parametric bootstrap (i.e., ZIPG-pbWald). That is, we simulate bootstrap samples from the ZIPG model with its parameters estimated under $H_{0}$. A detailed algorithm is provided in Section 2 of the Supplements.

There are multiple ways to construct test statistics and confidence intervals. We conducted extensive simulations studies evaluating the performance of (1) the Wald test without bootstrap (ZIPG-Wald) and the likelihood ratio test (ZIPG-LRT); (2) nonparametric and parametric bootstrap through different construction strategies, such as normality-based/quantile-based/$\rm BC_{a}$(Efron and Tibshirani, 1994) intervals; (3) resampling schemes for ZIPG-bWald, such as resampling based on measurements or subjects. Simulation results suggest that bootstrap-based Wald tests (i.e., ZIPG-bWald and ZIPG-pbWald) with normality-based confidence intervals are desired, and resampling based on measurements yields satisfactory results. More discussions are provided in Section 4.5.

We simulated $n=20$ subjects with $m=\{5,25\}$ measurements for each subject, with a total sample size of $N=\{100,500\}$. For $i=1,...,n$ and $j=1,...,m$, we generated covariates $\boldsymbol{X_{ij}}=(X_{i,1},X_{ij,2})$ associated with the mean parameter $\lambda_{ij}$ and $\boldsymbol{X^{*}_{i}}=X_{i,1}$ associated with the dispersion parameter $\theta_{i}$, where $X_{i,1}$ was a time-independent indicator sampled from a Bernoulli distribution with equal probabilities and served as the group indicator for each subject, such as pregnant or non-pregnant, alcohol drinker or non-alcohol drinker. $X_{ij,2}$ was a time-dependent longitudinal measurement generated from $X_{ij,2}=X_{i,2}+\epsilon_{ij}$, where $X_{i,2}\sim\mathcal{N}(0,1)$ and $\epsilon_{ij}\sim\mathcal{N}(0,0.1)$, representing covariates with variation in different measurements, such as calorie intake. Sequencing depths $M_{ij}$’s were generated based on the empirical distribution of observed sequencing depths in the Romero data set (Romero et al., 2014) analyzed later in Section 5. Finally, the observed count data $\boldsymbol{W}$ was generated based on models (3)-(4) with $\boldsymbol{\beta}$ and $\boldsymbol{\beta^{*}}$ specified as below.

To imitate real-world microbial data, we set $(\beta_{0},\beta^{*}_{0})=(-4.23,0.6)$, guided by the ZIPG estimates in Romero data. We investigate performance of ZIPG under different model parameter configurations (Table 2). In each simulation setting, we explored three different values of the zero-inflation proportion $p=\{0.3,0.5,0.7\}$, which is equivalent to $\gamma=\{-0.847,0,0.847\}$. For power analysis, we set $p=0.5$ and vary the parameter of interest under each null hypothesis from 0 to 1.8. We compare ZIPG-bWald, ZIPG-pbWald, pscl, NBZIMM, and PG for the inference of $\beta_{1}$ and only ZIPG-bWald, ZIPG-pbWald, and PG for the inference on $\beta^{*}_{1}$, because inference of the dispersion parameter is not applicable in pscl and NBZIMM. Results are presented based on $L=1,000$ Monte Carlo replicates for each scenario. The performances of ZIPG-bWald and ZIPG-pbWald are evaluated using $B=200$ bootstrap samples, which are numerically sufficient and stable based on our experience.

To demonstrate the advantage of ZIPG regarding point estimators and confidence intervals, we report the average bias (e.g., $\{\sum_{l}(\hat{\beta}_{l}-\beta)\}/L$) with its standard error over $l=1,\ldots,L=1000$ Monte Carlo replicates, average bootstrap standard error (avg-SE, e.g., $\{\sum_{l}\widehat{SE}(\hat{\beta}_{l})\}/L$), root mean squared error (RMSE, e.g., $\sqrt{\{\sum_{l}(\hat{\beta}_{l}-\beta)^{2})\}/L}$), and converge rate of confidence intervals (CR) for the settings with $\beta_{1}=0,\beta^{*}_{1}=1$ (Figure 2(a)) and $\beta_{1}=1,\beta^{*}_{1}=0$ (Figure 2(b)) with $p=\{0.5,0.7\}$ and $N=500$. Results under other settings are similar and hence not reported.

In Table 3, we observe that ZIPG-bWald has the smallest bias of $\beta_{1}$, $\beta_{1}^{*}$ and $\gamma$ among all methods when $\beta_{1}^{*}=1$. In addition, ZIPG is often more efficient than the other two methods, providing a smaller RMSE. For $\beta_{1}$ and $\beta_{1}^{*}$, ZIPG-bWald always maintains a valid confidence interval with its coverage rate close to nominal level 0.95, whereas pscl and NBZIMM provide underestimated confidence intervals in most cases, and PG estimate $\beta_{1}^{*}$ with strong bias and provide a more conservative CI for $\beta_{1}$. Additional results with setting $\beta_{1}=1,\beta^{*}_{1}=0$ corresponds to Figure 2(b) are presented in Section 3.4 of the Supplements.

While ZIPG assumes that all differential variability comes from the Poisson-Gamma part (i.e., $\theta$) and the true zero proportion $p$ is only taxon-specific, one may wonder how the model fits when these assumptions are violated. Here, we evaluate ZIPG under the misspecified model, in which the true zero proportion $p$ is also associated with covariates (denoted as “ZIPG-full”). Of note, ZIPG-full is equivalent to Omnibus (Chen et al., 2018) from a hypothesis testing perspective, whereas the Omnibus test does not provide point/interval estimation and hypothesis test for each parameter separately. We consider the group covariates $X_{1}$ only and sample size $N=500$. We use a logistic link ${\rm logit}(p)=\gamma_{0}+\gamma_{1}X_{1}$ with $(\beta_{0},\beta_{0}^{*},\gamma_{0})=(-4.23,0.6,-0.847)$, and then we report the performance of ZIPG with $\gamma_{1}$ increasing from 0 to 2.5, which is equivalent to increasing $p$ from $0.30$ to $0.84$ in the group with $X_{1}=1$.

For two model fittings, i.e., ZIPG and ZIPG-full, we report the proportion of simulations suggesting better BIC from ZIPG (Figure 4). Over 1000 replicated simulations in each setting, ZIPG has a smaller BIC than ZIPG-full in most cases. Moreover, we also use Kolmogorov-Smirnov test (Kolmogorov, 1933; Smirnoff, 1939) to compare the ZIPG predicted distribution with the simulated observations: 100% replicates report insignificant differences between the two distributions ($p>0.05$) , suggesting that ZIPG-predicted distribution has no difference to the observed samples.

We also evaluate ZIPG’s performance under other misspecified models: (1) Poisson-Gamma without zero inflation and (2) zero-inflated Beta-Binomial model (see Section 3.5 of the Supplements). In both scenarios, ZIPG preserves nominal type I error and retains its superior power in detecting differential abundance/variability, especially for large-sample cases.

We conduct additional simulations to compare multiple ways of constructing test statistics and confidence intervals. For hypothesis testing, we investigate different test statistics for the proposed ZIPG (Section 3.1 of the Supplements), including the Wald test without bootstrap (ZIPG-Wald) and the likelihood ratio test (ZIPG-LRT). Simulation results suggest that ZIPG with bootstrap-based Wald tests (i.e., ZIPG-bWald and ZIPG-pbWald) are desired. For confidence intervals, we compare the coverage rate for both nonparametric and parametric bootstrap through different construction strategies, such as normality-based/quantile-based/$\rm BC_{a}$(Efron and Tibshirani, 1994) intervals. Results show that the normality-based confidence interval often has close-to-nominal coverage with a low computational cost (Section 3.2 of the Supplements).

We further evaluate different resampling schemes for ZIPG-bWald, such as resampling based on measurements or subjects. Given the total sample size $N=200$, results suggest no significant difference between the two strategies (Section 3.6 of the Supplements).

Additional simulations about how measurement times affect ZIPG performance are provided in Section 3.6 of the Supplements. We consider the following two scenarios: (1) when the numbers of measurements per subject are very small (i.e., $m=2$) and (2) when the numbers of measurements per subject are unequal. Results show that the proposed ZIPG method is valid for both cases above.

We first present numbers of identified taxa regarding the covariates of interest in the two studies (Figure 5). For ZIPG, we show the set of taxa associated with $\boldsymbol{X}$ in the mean model (denoted as ZIPG $\beta$) and $\boldsymbol{X^{*}}$ in the dispersion model (denoted as ZIPG $\beta^{*}$), separately. In Romero, we observe that pregnant subjects are clustered under age 35, while non-pregnant subjects are collected in a much wider range of ages. Owing to the unbalanced sample that pregnant women have fewer measurements and are often younger than non-pregnant women, we use parametric bootstrap (i.e., ZIPG-pbWald) for more stable results. ZIPG identified 18 taxa with differential abundance ($H_{0}:\beta_{1}=0$) and 17 taxa with differential variability ($H_{0}:\beta_{1}^{*}=0$) associated with pregnancy, while 13 taxa are associated with pregnant status with both abundance and variability (Figure 5(a)). Most of the taxa found by pscl and DESeq2 are also identified by ZIPG, while ZIPG also identified 3 additional taxa with significant differential variability, which are not detected by other methods. Further, ZIPG identified totally additionally 6 taxa associated to age with differential abundance ($H_{0}:\beta_{2}=0$) or differential variability ($H_{0}:\beta_{2}^{*}=0$) , compared to other methods (Figure 5(b)).

In Dietary, ZIPG is performed using nonparametric bootstrap Wald test (i.e., ZIPG-bWald), because the data is balanced and the sample size is sufficient. ZIPG identified 33 taxa with differential abundance ($H_{0}:\beta_{1}=0$) and 24 taxa with differential variability ($H_{0}:\beta_{1}^{*}=0$) associated to the alcohol intake (Figure 5(c)). Compared to other methods, ZIPG discovered 5 extra taxa only associated with differential variability and 1 extra taxon associated with both differential abundance and differential variability.

We further use the Omnibus test (Chen et al., 2018) as verification for ZIPG detected taxa. The Omnibus test links covariates of interest to all three parameters in the Negative-Binomial distribution and rejects the null hypothesis if any covariate is associated with any of the parameters. Thus, taxa identified by both ZIPG and Omnibus are less likely to be false positives. As expected, in Romero, all taxa identified by ZIPG are also detected by Omnibus. In Dietary, 35 out of 41 taxa identified by ZIPG are also detected by Omnibus. Though the Omnibus test detected more taxa than ZIPG, it is worth pointing out that the Omnibus test can not distinguish differential abundance and differential variability and does not provide point/interval estimation as ZIPG does. Details of taxa detected by each method and estimation results for those taxa regarding parameters of interest are shown in Section 4.3 of the Supplements.

To visualize the differential variability tested by ZIPG (i.e., $H_{0}:\beta^{*}=0$), we further analyze the result of Bifidobacteriaceae and Lactobacillus.vaginalis from Romero, and Burkholderiales bacterium and Alistipes indistinctus from Dietary as examples. Other taxa with differential variability identified by ZIPG have similar conclusions.

First, we compare the goodness of fit of results from fitted models to the empirical distribution of the relative abundance. The log of relative abundance observed in real data (i.e., Bifidobacteriaceae from Romero) is compared to the predicted distribution according to the estimated model by ZIPG, pscl, and DESeq2 (Figure 6). For boxplots, we generated samples from each predicted distribution with a quintuple of the observed sample sizes for a better visualization at the tail (e.g., Figure 6(6(a))). All three methods can estimate the median of two groups (pregnant and non-pregnant) accurately, but pscl and DESeq2 cannot distinguish the overdispersion between the two groups, as the box length for the pregnant and non-pregnant groups are similar. On the contrary, ZIPG identified the differential variability of this taxon associated with the factor pregnant with $p=0.00256$ regarding the null hypothesis $H_{0}:\beta_{1}^{*}=0$, and thus its simulated data matches the real data better, providing a shorter interquartile range with a long tail in the pregnant group. We also present the empirical cumulative distribution functions (ECDF) of the log of the relative abundance for the pregnant group, using sampled data from fitted ZIPG, pscl, and DESeq2 models (Figure 6(6(b))). It has been shown that ZIPG also fits the real data better than other methods. Quantile-quantile plots for real data versus predicted distribution also show that ZIPG models the entire distribution better (Section 4.2 of the Supplements).

For Lactobacillus.vaginalis in Romero, we present the relative abundance based on the simulated distributions of the fluctuation factor $U\sim{\rm Gamma}(\theta_{i}^{-1},\theta_{i})$ from ZIPG and pscl at four representative ages (Figure 7). ZIPG identified that the differential variability of this taxon is associated to age with $p$ ¡ 0.001 regarding $H_{0}:\beta_{2}^{*}=0$. Accordingly, we observe that the shape of the fitted distribution changes with the increase of age. However, pscl is not able to model the change in the entire shape of distribution as the average relative abundance remained the same in this group. The ECDF plot of the pregnant group again shows that ZIPG fitted model is closer to the empirical distribution.

In Dietary, we present the box plots for Burkholderiales bacterium and Alistipes indistinctus to show the differential variability between the alcohol and non-alcohol drinkers (Figure 8). For Burkholderiales bacterium, the differential abundance in the two groups (i.e., alcohol and non-alcohol drinkers) are similar, but the overdispersion in alcohol drinkers is obviously larger than that in teetotal subjects. ZIPG can distinguish the differential variability in two groups with $p=0.00293$ ($H_{0}:\beta_{1}^{*}=0$). ZIPG also provides a shorter interquartile in the non-alcohol drinker group, which is consistent with the raw data. On the contrary, pscl failed to detect any differential abundance ($p=0.636$), while DESeq2 provided $p=0.049$ which is significant but much larger than the p-value of ZIPG. For Alistipes indistinctus, though pscl and DESeq2 identified the differential abundance between two groups, both of them did not approximate the data from a distributional perspective because of the ignorance of differential variability. In contrast, ZIPG can distinguish it with $p$ = 8.54e-6 ($H_{0}:\beta_{1}^{*}=0$), and provide a long box for the non-alcohol group showing their small overdispersion and a median more closer to real data.