# Bayesian Hierarchical Compositional Models for Analysing Longitudinal Abundance Data from Microbiome Studies

Authors: I. Creus Martí; A. Moya; F. J. Santonja
Journal: Complexity
 (2022)

---

## Abstract

Gut microbiome plays a significant role in defining the health status of subjects, and recent studies highlight the importance of using time series strategies to analyse microbiome dynamics. In this paper, we develop a Bayesian model for microbiota longitudinal data, based on Dirichlet distribution with time-varying parameters, that take into account the compositional paradigm and consider principal balances. The proposed model can be effective for predicting the future dynamics of a microbial community in the short term and for analysing the microbial interactions using the value of the estimated parameters. The usefulness of the proposed model is illustrated with six different datasets, and a comparison study with four alternative models is provided.


---

## Body

## 1. Introduction

The microbiome is defined in [1] as “the ecological community of commensal, symbiotic, and pathogenic microorganisms that literally share our body space.” Recent studies suggest that gut microbiome plays a significant role in defining the health status of subjects. In fact, human gut microbiome influences brain function and behavior [2] and is related to several diseases such as cancer [3], HIV [4], Alzheimer [5], type 2 diabetes [6], cardiovascular disease [7], irritable bowel syndrome (IBS) [8], Crohn’s disease [9], chronic fatigue syndrome (CFS) [10], obesity [11], and malnutrition [12], among others. What is more, in [13], the paramount interest of gut microbiome, pointing out the importance of microbes in regulating most if not all human bodily functions, has been underscored. As a result, monitoring human microbiome will probably lead us to improve human health, as has been mentioned in [14].Given the importance of the microbiome, a concerted effort is being made to design models that enable us to extract information on bacterial dynamics. In [15], Martí et al. explain that temporal microbiota stability is linked to health status. In [16], Dinleyici et al. observed that a disruption in the stable state of the microbiome is related to a disease, but the stable state is reached again when the subject is cured. These facts highlight the importance of using time series strategies to study the microbiome. In that respect, longitudinal studies are an appropriate approach to tackle microbiome modeling, as they can help to gain a better understanding of microbiota regularities, bacterial interactions, and microbiome response to perturbations. Furthermore, this longitudinal approach allows us to predict future microbiome dynamics and analyse not only microbiome stability but also the time it needs to recover and reach a new stable state, if possible.At present, there are studies where the microbiome is manipulated in humans in order to treat certain illnesses [17]. This manipulation is carried out by fecal microbiota transplantation, probiotic supplementation, or a combination of both strategies [17]. Nowadays, there are also nutrition-based approaches to treatment [18]. For instance, the work detailed in [19] points out that the administration of Lactobacillus delbrueckii and antibiotics to treat bacterial vaginosis can provide a cure, [20, 21] indicate that fecal microbiota transplantation is a promising therapy for C. difficile infection, and [22] reports that fecal microbiota transplantation in HIV patients attenuates HIV-associated dysbiosis. Taking this research approach into account, modeling bacterial dynamics should be a research priority because these models will help us to design better interventions, as has been mentioned in [17].One line of action in longitudinal analysis to model the abundance of each microbial entity at each time point is to use generalized Lotka–Volterra equations [23–26]. We must mention that these models have usually been used in a deterministic way; however, [27] takes stochastic processes into account, using multiplicative noise with normal distribution. These approaches consider pair-wise interactions, and consequently, they do not take into account the effect that a third microbe may have on an interacting pair of microbes. In [28], more limitations of using the Lotka–Volterra structure have been pointed out. In [29, 30], a nonparametric approach with an additive structure has been proposed. The advantage of these models is the absence of microbial relationship specifications, while the drawback is the fact that the additivity of relationships, allowed in this approach, may not be realistic for complex microbial communities. Other proposals consist of using state space models [31] or linear mixed models [32, 33]. In [34], and references therein, more proposal to analyse microbiota longitudinal data is detailed.Recent investigations consider compositional data analysis [35] is also a valid approach to analyse microbiome high-throughput sequencing data, and therefore, microbiome datasets can be converted to relative abundance values, or normalized counts, prior to analysis (see, for instance, [36] or Chapter 3 in [34]). Fundamentally, this approach uses the log-ratio transformation techniques to convert microbiome data, thus removing the compositional constraints to make the standard techniques suitable for analysis. This transformation focuses on removing the compositional constraint: all microbial relative abundances sum to one. Microbiome data can be considered as a compositional dataset due to limitations of the techniques used to collect the data [36]. These limitations are based on the fact that datasets collected by 16S rRNA gene amplimers have an arbitrary total imposed by the instrument. This total is known as sequencing depth or library size. This implies that an increase in the abundance of one taxon results in a decrease of some of the other taxa because the sum of all of the raw abundances cannot exceed the sequencing depth. As the increase of a taxon produces a decrease in some of the other taxa, the correlation between the taxa is bogus. The work presented in [37] specified that the correlations between compositional components are spurious. This fact was described under various headings: the constant-sum problem, the closure problem, the negative bias problem, or the null correlation difficulty [38]. These characteristics of compositional data require that the analysis of this type of data must be conducted using appropriate methodology.There are approaches that analyse microbiome taxa taking its compositional nature into account. In [39], Kynčlová et al. model the isometric log-ratio transformation of the data using Gaussian distribution and a VAR structure, while in [40] phylogenetic isometric log-ratio transformation, to allow off-the-shelf statistical tools to be optimally applied to microbiota datasets, is introduced. This approach takes into account phylogenies to construct the ILR transform and redefine data considering coordinates termed balances. Rivera-Pinto et al. in [41], take also into account balances. In this case, to identify microbial signatures and consider them as predictive covariates of a phenotype of interest, Joseph et al. describe microbial dynamics using generalized Lotka–Volterra (gLV) equations and compositional data analysis (CODA) [42]. In [43, 44], Bayesian CODA methods for microbiome time series which model the counting process and perform inference and apply log-ratio transformations to the model parameters are detailed.Our proposal is based on modeling a normalized transformation of the observed counts. We assume that normalized transformation of the observed counts can be explained by an autoregressive structure considering Dirichlet distribution with time-varying parameters that take into account principal balances. Principal balances are a compositional tool that enables us to reduce the number of parameters of the model by selecting the balances that maximize the variance [45]. As balances with more variability will better explain the variability of the data, this formulation allows us to encapsulate information and decrease the number of parameters to be estimated. Our proposal is in line with [46], though this one is not used with microbiome datasets. We must mention that the Dirichlet distribution has been considered inappropriate for modeling compositional data due to its strong implied independence structure, which can be deduced from its definition by a set of independent, gamma-distributed, random variables, with equal scale parameter. This fact makes it inappropriate to consider this probability distribution for modeling compositional data [35]. However, it has also been found useful when used as a conditional distribution (see, for instance, [47, 48]). Our approach is presented as an alternative for predicting the future dynamics of a microbial community in the short term and for analysing the microbial interactions using the value of the estimated parameters.


## 2. Some Basic Principles of CODa

The first steps towards developing a suitable methodology to analyse compositional data focused on projecting the sample space of compositional data (the simplex) to real space using log-ratio transformations [35]. The additive log-ratio transformation (alr) and centered log-ratio transformation (clr) were defined in [35] and are defined by the following expressions:(1)alrx=logx1xD,logx2xD,…,logxD−1xD,clrx=logx1gx,logx2gx,…,logxDgx,where x=x1,x2,…,yD denotes a D-composition and g⋅ denotes the geometric mean. Note that when projecting the simplex to the real space, the clr transformation preserves the distances [49]; see equation (2), where dax,y denotes the Aitchison distance between two vectors x and y and de⋅,⋅ denotes the Euclidean distance.(2)da2x,y=12D∑i=1D∑j=1Dlogxixj−logyiyj2=de2clrx,clry.


### 2.1. Balances

Balances are defined in [50]. Using the notation applied in [45], they are expressed as follows:(3)b=r⋅sr+s⋅lngx+gx−,where x+ and x− are two nonoverlapping groups of parts of a D-composition, r and s are the number of parts in x+ and x−, respectively, and g⋅ denotes the geometric mean.Balances are a compositional tool that makes it possible to study the relationship between the components of the numerator and the components of the denominator of the balance. Balances that are near zero indicate similarity between the two groups; however, the greater the absolute value of the balance, the lager is the dissimilarity between them.Depending on which components we put inx+ and which we put in x− we will have different balances. In this study, we want to choose the components of x+ and the components of x− that maximize the variance of the balance because balances with more variability will better explain the variability of the microbiome data. Principal balances are the strategy selected to tackle this issue.


### 2.2. Principal Balances

Principal balances act as the Principal Component Analysis (PCA) for compositional data. In other words, they allow us to reduce the dimension with a minimum loss of information. In addition, the numerator and the denominator of the principal balances are selected so as to maximize the variance of the balance [45].In order to obtain the principal balances, Ward’s clustering method [51] is used (see details in [52]). We apply this method because it enhances the proportionality between groups of parts and reduces the computation time [52].


## 2.1. Balances

Balances are defined in [50]. Using the notation applied in [45], they are expressed as follows:(3)b=r⋅sr+s⋅lngx+gx−,where x+ and x− are two nonoverlapping groups of parts of a D-composition, r and s are the number of parts in x+ and x−, respectively, and g⋅ denotes the geometric mean.Balances are a compositional tool that makes it possible to study the relationship between the components of the numerator and the components of the denominator of the balance. Balances that are near zero indicate similarity between the two groups; however, the greater the absolute value of the balance, the lager is the dissimilarity between them.Depending on which components we put inx+ and which we put in x− we will have different balances. In this study, we want to choose the components of x+ and the components of x− that maximize the variance of the balance because balances with more variability will better explain the variability of the microbiome data. Principal balances are the strategy selected to tackle this issue.


## 2.2. Principal Balances

Principal balances act as the Principal Component Analysis (PCA) for compositional data. In other words, they allow us to reduce the dimension with a minimum loss of information. In addition, the numerator and the denominator of the principal balances are selected so as to maximize the variance of the balance [45].In order to obtain the principal balances, Ward’s clustering method [51] is used (see details in [52]). We apply this method because it enhances the proportionality between groups of parts and reduces the computation time [52].


## 3. Model

### 3.1. Selecting Principal Balances

Letyt=y1t,y2t,…,yDt be the relative abundances at time t, so that yit∈0,1 and ∑i=1Dyit=1. Note that yit denotes the relative abundance of taxon i at time point t. Let T be the number of time points available in the data. In order to implement Ward’s clustering method, we denote as Si the vector that contains the clr transformation of the relative abundance of species i at all time points available in the data. As a result,(4)Si=logyi1gy1,logyi2gy2,…,logyiTgyTwithi=1,2,…,D.First of all, we calculate the matrixM=mijD−1xD−1wheremij=deSi,Sj. This matrix is proportional to the variation matrix due to the following equation defined by [52].(5)Varlogzizj=1Tda2zi,zj=1Tde2Si,Sj.In this equation,zi=yi1,yi2,…,yiT denotes the vector that contains the relative abundance of species i at all time points available in the data. Before going into further detail, it is important to point out that Ward’s method distance between two clusters B and C is defined by equation (6), where C¯ and B¯ denote the centroids (i.e., barycenters) of C and B, respectively.(6)dC,B=2C‖BC+BdeC¯,B¯.Second, we apply Ward’s method, because according to [52], the variance of a balance is proportional to Ward’s method distance between the numerator and the denominator of the balance. Ward’s method is hierarchical, and it detects the smallest entry in the matrix M (which functions as the variation matrix), and the corresponding parts are merged to form a group. The matrix M is recalculated at each iteration using equation (6) and, the final stage of the process, the last two remaining groups are fused into one, which gives the balance with the largest variance. We have implemented this process using the hclust function included in the fastcluster package in R [53].During the hierarchical process, we generate a dendrogram that plots the groups created in each iteration of the method. As a result, the groups merged in the dendrogram give us the numerator and the denominator of the principal balances. In conclusion, through this process, we have applied the new approach to Ward’s clustering method in order to obtain the principal balances, which are denoted asPBal1,PBal2,…,PBalD−1.Finally, we select the principal balances for which the sum of the percentage of variance is higher than80%. We choose this percentage because we consider that taking account of the balances that describe at least 80% of the variability is enough for our model to capture the variability of the data [54]. We denote the selected principal balances as SPBal1,SPBal2…,SPBalM, where M denotes the total amount of SPBal. We must mention that a detailed example of the process described in this section to obtain the selected principal balances is provided in Section 5.


### 3.2. Proposed Model (PM)

We assume thatyt follows a Dirichlet distribution with nonnegative parameters αt=α1t,α2t…,αDt. Reference [46] presents the ability of the Dirichlet distribution with an autoregressive structure that uses the alr or the clr transformation as a link function and imposes a parameter that is equal to the sum of the Dirichlet parameters. Our proposal is in line with this because it uses the Dirichlet distribution and an autoregressive structure; however, we use a Bayesian framework with the logarithm as a link function, we introduce the principal balances in the autoregressive structure, and we do not need to estimate an extra parameter that is equal to the sum of the Dirichlet parameters.We model each transformed parameterlogαit as a linear function of the selected principal balances:(7)logαit=ai0+∑j=0Maij⋅SPBaljt−1.Note thatSPBalj,t−1 are calculated using the values of y^i,t−1, and the parameters of the model are aij. The Gaussian distribution with zero mean and standard deviation σij is considered as a prior distribution of aij. As [55] suggests, σij is assigned the uniform distribution in the interval 0,5. Note that αit=expai0+∑j=0Maij⋅SPBaljt−1 and then Eyit can be estimated by Eyit=αit/∑i=1DαitAs we have mentioned before, the numerator and the denominator of the principal balances are selected so as to maximize the variance of the balance. As balances with more variability will better explain the variability of the data, principal balances enable us to encapsulate information and decrease the number of parameters to be estimated. Balances provide information about the similarity between the components of the numerator and the elements of the denominator: balances that are approximately zero indicate that the two groups are similar; however, the larger the absolute value of the balances, the greater the dissimilarity between them.


## 3.1. Selecting Principal Balances

Letyt=y1t,y2t,…,yDt be the relative abundances at time t, so that yit∈0,1 and ∑i=1Dyit=1. Note that yit denotes the relative abundance of taxon i at time point t. Let T be the number of time points available in the data. In order to implement Ward’s clustering method, we denote as Si the vector that contains the clr transformation of the relative abundance of species i at all time points available in the data. As a result,(4)Si=logyi1gy1,logyi2gy2,…,logyiTgyTwithi=1,2,…,D.First of all, we calculate the matrixM=mijD−1xD−1wheremij=deSi,Sj. This matrix is proportional to the variation matrix due to the following equation defined by [52].(5)Varlogzizj=1Tda2zi,zj=1Tde2Si,Sj.In this equation,zi=yi1,yi2,…,yiT denotes the vector that contains the relative abundance of species i at all time points available in the data. Before going into further detail, it is important to point out that Ward’s method distance between two clusters B and C is defined by equation (6), where C¯ and B¯ denote the centroids (i.e., barycenters) of C and B, respectively.(6)dC,B=2C‖BC+BdeC¯,B¯.Second, we apply Ward’s method, because according to [52], the variance of a balance is proportional to Ward’s method distance between the numerator and the denominator of the balance. Ward’s method is hierarchical, and it detects the smallest entry in the matrix M (which functions as the variation matrix), and the corresponding parts are merged to form a group. The matrix M is recalculated at each iteration using equation (6) and, the final stage of the process, the last two remaining groups are fused into one, which gives the balance with the largest variance. We have implemented this process using the hclust function included in the fastcluster package in R [53].During the hierarchical process, we generate a dendrogram that plots the groups created in each iteration of the method. As a result, the groups merged in the dendrogram give us the numerator and the denominator of the principal balances. In conclusion, through this process, we have applied the new approach to Ward’s clustering method in order to obtain the principal balances, which are denoted asPBal1,PBal2,…,PBalD−1.Finally, we select the principal balances for which the sum of the percentage of variance is higher than80%. We choose this percentage because we consider that taking account of the balances that describe at least 80% of the variability is enough for our model to capture the variability of the data [54]. We denote the selected principal balances as SPBal1,SPBal2…,SPBalM, where M denotes the total amount of SPBal. We must mention that a detailed example of the process described in this section to obtain the selected principal balances is provided in Section 5.


## 3.2. Proposed Model (PM)

We assume thatyt follows a Dirichlet distribution with nonnegative parameters αt=α1t,α2t…,αDt. Reference [46] presents the ability of the Dirichlet distribution with an autoregressive structure that uses the alr or the clr transformation as a link function and imposes a parameter that is equal to the sum of the Dirichlet parameters. Our proposal is in line with this because it uses the Dirichlet distribution and an autoregressive structure; however, we use a Bayesian framework with the logarithm as a link function, we introduce the principal balances in the autoregressive structure, and we do not need to estimate an extra parameter that is equal to the sum of the Dirichlet parameters.We model each transformed parameterlogαit as a linear function of the selected principal balances:(7)logαit=ai0+∑j=0Maij⋅SPBaljt−1.Note thatSPBalj,t−1 are calculated using the values of y^i,t−1, and the parameters of the model are aij. The Gaussian distribution with zero mean and standard deviation σij is considered as a prior distribution of aij. As [55] suggests, σij is assigned the uniform distribution in the interval 0,5. Note that αit=expai0+∑j=0Maij⋅SPBaljt−1 and then Eyit can be estimated by Eyit=αit/∑i=1DαitAs we have mentioned before, the numerator and the denominator of the principal balances are selected so as to maximize the variance of the balance. As balances with more variability will better explain the variability of the data, principal balances enable us to encapsulate information and decrease the number of parameters to be estimated. Balances provide information about the similarity between the components of the numerator and the elements of the denominator: balances that are approximately zero indicate that the two groups are similar; however, the larger the absolute value of the balances, the greater the dissimilarity between them.


## 4. Datasets

We demonstrated the ability of our model by using both publicly available and simulated datasets. We studied the gut microbiome time series data of two healthy subjects, one female and one male. In both subjects, the samples were collected daily. The data of the female subject were extracted from [56] and were analysed on different taxonomic levels. We used Subject C in [57] as the male subject, taking into account that this subject is obese class I according to his body mass index.We simulated three databases with fifteen, twenty, and sixty microbial taxa, respectively. Following the scheme given by [58], we generated the interaction matrix using the algorithm proposed by [59], and we generated the initial abundances using the Poisson distribution. With these two tools, we simulated the data using the generalized Lotka–Volterra structure. We carried out the simulation using the R package seqtime [58]. Focusing on technical details, to generate the interaction matrix we set the clique size at 10, the diagonal values at −1, the interaction connectance at 0.01, the positive edge percentage at 64%, and the maximal absolute off-diagonal interaction strength at 1.Microbial datasets are zero-inflated [60]; nevertheless, zeros in a compositional data set are below detection limit values [61]. For this reason, in order to analyse the publicly available datasets, we selected the bacterial taxa that did not have zeros at any time point, and we gathered the information of the nonselected taxa in the Other collection.


## 5. Results

First of all, we will explain in detail how we applied Ward’s clustering method to obtain the principal balances and present a comparative study on selecting the principal balances. Second, we analyse the information that estimated parameters provide on bacterial dynamics. Third, we offer a study on the presence of zeros at credible intervals. After this, we examine the graphic results and test the results when modeling datasets with a large number of microbial taxa. Finally, we compare the proposed model with other approaches.


### 5.1. Applying Ward’s Clustering Method and Selecting the Principal Balances

In order to introduce the results in an accessible way, we first present the outcomes obtained with the publicly available datasets. Moreover, we start by analysing the female dataset at family level, because it has aspects in common with the other two datasets and, consequently, it is easily comparable.By applying Ward’s clustering method as it is described in Section 3.1, we obtain the dendrogram shown in Figure1. As it is presented by [62], the horizontal bars of the dendrogram describe the groups of parts defined at each iteration of the process; for this reason, the groups merged in the dendrogram give us the numerator and the denominator of the principal balances. Moreover, the sum of all vertical bars represents the total variance of the sample; therefore, a long vertical bar implies a balance that explains a good deal of the cc variance and a short vertical bar means that balance has a little variability in the sample. In Figure 1, we can also see the principal balances described by the dendogram and their percentage of variance (Table1).Figure 1
Female dataset at family taxonomic level. Dendogram obtained using Ward’s method for obtaining the principal balances. The table shows the expression of the principal balances and their percentage of variance.Table 1
The expression of the principal balances and their percentage of variance.
PBalExpression% varPBal11⋅1/1+1⋅lnx3/x69.00PBal22⋅1/2+1⋅lnx3⋅x61/2/x715.38PBal31⋅1/1+1⋅lnx5/x92.15PBal42⋅1/2+1⋅lnx5⋅x91/2/x47.62PBal53⋅1/3+1⋅lnx5⋅x9⋅x41/3/x110.62PBal61⋅1/1+1⋅lnx2/x85.31PBal74⋅2/4+2⋅lnx5⋅x9⋅x4⋅x11/4/x2⋅x81/214.08PBal86⋅3/6+3⋅lnx5⋅x9⋅x4⋅x1⋅x2⋅x81/6/x7⋅x3⋅x61/335.81As explained in Section3.1, we select the principal balances that accounted for at least 80% of the variability. Consequently, we pick PBal8, PBal2, PBal7, PBal5, and PBal1 because they explain a cumulative variance of 84.90% (see Table 2). Note that we have written the principal balances in descending order by percentage of variance and we denote them as SPBal1,SPBal2…SPBal5, respectively. For comparative purposes, we have estimated the proposed model using different selected principal balances to compare the results when extracting different percentages of variance of the data. In Table 2, we can observe that the formulation with SPBal1,SPBal2…SPBal5 is the one with smallest value of 2002 DIC, 2004 DIC, and RMSD when fitting; consequently, this formulation is the most suitable option. In fact, it is underscored in Section 5.2 that SPBal5 is particularly interesting for explaining the dynamics of the microbial community.Table 2
Female dataset at family taxonomic level. Results obtained when estimating the proposed model using different selected principal balances. Additionally, the table shows the deviance information criterion (DIC) defined in [63], the DIC defined in [64], the root mean square deviation (RMSD) of the results obtained fitting the model, and the RMSD of the results obtained using the model for prediction.
Proposed model different SPBal% varD¯pD (2002)DIC (2002)pD (2004)DIC (2004)RMSD (fitting)RMSD (predicting)SPBal135.81−1820.7616.93−1803.8317.60−1803.150.02990.0408SPBal1, SPBal251.19−1827.1124.16−1802.9431.55−1795.560.02920.0409SPBal1, SPBal2, SPBal365.27−1826.6131.08−1795.5240.64−1785.960.02880.0442SPBal1,SPBal2,…,SPBal475.89−1832.6537.68−1794.9645.90−1786.740.02820.0452SPBal1,SPBal2,…,SPBal584.90−1872.2345.48−1826.7564.49−1807.740.02640.0498Table 3
Estimated parametersaij and credible interval obtained by applying the proposed model to the female dataset at family taxonomic level. The table also contains the percentage of influence that each SPBalj has on each family and the percentage of sway that each family has on the bacterial community is represented in bold.
Family taxonomic levelaijInterceptSPBal1SPBal2SPBal3SPBal4SPBal5o_bacteroidales; f_−0.03660.42550.16660.0145−0.3427−0.3303(−1.06, 0.96)(0.13, 0.69)(−0.11, 0.46)(−0.33,0.37)(−0.76, 0.04)(−0.69, 0.009)0.58%33.25%13.02%1.13%26.78%25.81%Other1.10760.7105−0.07530.09130.2463−0.1231(0.25, 1.95)(0.44, 0.96)(−0.31, 0.13)(−0.20, 0.39)(−0.04, 0.53)(−0.47, 0.18)17.50%57.00%6.04%7.32%19.76%9.88%f_erysipelotrichaceae−0.5757−0.04970.1025−0.01910.01410.0055(−2.05, 0.61)(−0.40, 0.25)(−0.32, 0.54)(−0.50, 0.47)(−0.53, 0.56)(−0.50, 0.46)9.10%26.03%53.69%10.01%7.39%2.88%o_clostridiales; f_−0.87050.1802−0.03460.0038−0.02810.0198(−2.44, 0.38)(−0.15, 0.54)(−0.46, 0.35)(−0.45, 0.48)(−0.56, 0.50)(−0.50, 0.50)13.75%67.62%12.98%1.43%10.54%7.43%f_lachnospiraceae0.48470.3356−0.01990.0123−0.0225−0.5353(−0.33, 1.47)(0.04, 0.61)(−0.27, 0.21)(−0.29, 0.31)(−0.36, 0.30)(−0.88, −0.19)7,66%36,26%2,15%1,33%2,43%57,83%f_ruminococcaceae0.2661−0.20230.1230−0.21090.0594−0.7166(−0.71, 1.44)(−0.52, 0.06)(−0.19, 0.44)(−0.56, 0.10)(−0.30, 0.43)(−1.11, −0.30)4.20%15.42%9.37%16.07%4.53%54.61%f_veillonellaceae−0.20570.2451−0.1049−0.0496−0.0238−0.3804(−1.26, 0.75)(−0.05, 0.54)(−0.45, 0.16)(−0.42, 0.30)(−0.41, 0.34)(−0.76, −0.02)3,25%30,49%13,05%6,17%2,96%47,33%f_desulfovibrionaceae−0.63600.29930.03630.0866−0.1100−0.0535(−2.03, 0.45)(−0.03, 0.63)(−0.34, 0.40)(−0.32, 0.55)(−0.62, 0.36)(−0.52, 0.35)10,05%51,10%6,20%14,79%18,78%9,13%f_bacteroidaceae2.14600.2476−0.0018−0.1890−0.0426−0.5470(1.30, 3.04)(−0.008, 0.49)(−0.25, 0.22)(−0.47, 0.08)(−0.33, 0.25)(−0.87, −0.24)33,91%24,090,1818,394,1453,21


### 5.2. Estimated Parameters

Once we have selected the principal balances, we can estimate the proposed model. We are going to do so using Bayesian analysis because it allows us to simulate the posterior distribution of the model parameters using Markov Chain Monte Carlo (MCMC), as implemented in the free statistical software JAGS [65]. We fixed a burn-in period of 10000 iterations and an adaptation period of 1000 iterations. We kept one posterior sample in 100 iterations until a set of 20000 iterations was obtained, except for the first and second simulation datasets where just 2000 iterations were needed. We implemented this process using the runjags package in R [66]. Note that in Table S1 (see supplementary material), we can find the potential scale reduction factor defined by [67] and obtained using the gelman.diag function of the Coda package [68]. Table S1 shows that the potential scale reduction factor belongs to the interval 1,1.1 and is smaller than its upper 95% confidence limit.In Table3 we can find the estimated parameters, denoted in Section 3.2 as aij. Note that in each equation in the model (see equation (7)), we need to estimate six parameters, the intercept, and the five parameters associated with the selected principal balances. The estimated value of ai0 represents the abundance of the family i in the community; in fact, tiny intercepts are related to bacterial families of low abundance, whereas large intercepts indicate taxa of high abundance. In addition, the absolute value of the estimated aijwithj≠0 represents the weight that the balance SPBalj has in describing the relative abundance of the family i. As a result, we can say that the absolute value of aijwithj≠0 represents the influence of the balance SPBalj on the relative abundance of the family i.Analysing the estimated parameters allows us to extract some information about the bacterial dynamics. Focusing on intercepts described in Table3, we observe that the families o_Clostridiales; f_, f_Desulfovibrionaceae, and f_Erysipelotrichaceae have a low intercept, while the intercepts of f_Bacteroidaceae and Other are meaningful. A negative intercept indicates families of low abundance, whereas large positive intercepts indicate abundant ones. We define the percentage of sway that each family has on the bacterial community as ai0/∑i=1Dai0. This percentage of sway provides information about the abundance of a family taking into account the other families; consequently, it enables us to compare the abundance of the different families. When the intercept is positive, a high percentage indicates high abundance and vice versa. However, when the intercept is negative, high percentage is related to low abundance and vice versa. In Table 3, the percentage of sway is written in bold, and we can observe that f_Bacteroidaceae and Other, which have positive intercepts and a percentage of sway of 33.91% and 17.50%, respectively, are the two most abundant families. In addition, o_Clostridiales; f_ has a negative intercept and a percentage of 13.75%; it is therefore the least abundant family. The percentage of sway can help to gain a better understanding of the different dynamics of microbes because low percentage indicates small families, and high percentage is related to the extreme families: it indicates the most abundant if the intercept is positive or the least abundant if the percentage is negative.As we have mentioned before, the absolute value of the estimated parametersaijwithj≠0 represents the weight that the balance has in describing the relative abundance of the family, and they therefore provide information about the influence that each balance has on the family. Centering our attention on the absolute value of the estimated aij, we can observe that SPBal5 is the one that most influences f_Ruminococcaceae, f_Lachnospiraceae, f_Veillonellaceae, and f_Bacteroidaceae. Consequently, the families that belong to the SPBal5 are particularly interesting to explain the dynamics of these families. Note that SPBal5,t−1=1⋅1/1+1⋅lny3,t−1/y6,t−1; therefore, this balance explains the relationship between the components y3 and y6 (which are f_Erysipelotrichaceae and f_Ruminococcaceae respectively) at time point t−1. A balance value close to zero can be interpreted as a high similarity between the numerator elements and the denominator components; however, the larger the absolute value of the balance, the greater the dissimilarity between them. Furthermore, in Table 3, we can also observe that the coefficient of SPBal1 has a significant absolute value in all families. As a result, SPBal1, which is the selected principal balance with the highest percentage of variance, is also the one that has a significant influence on most of the bacterial families.Families with low abundances will have low parameter values; as a result, in order to compare the influence of oneSPBalj0 in different families, it is interesting to study the percentage of influence aij/∑j=1Maij. Observing Table 3, it seems that SPBal1 does not have an special influence on f_Erysipelotrichaceae because aij<0.1, while in the rest of the families aij>0.1, however, if we observe the value of all aij at f_Erysipelotrichaceae, we can observe that all the values are small and, in comparison, SPBal1 has a significant influence. This is the reason why it is interesting to study the percentage of influence present in Table 3. In f_Erysipelotrichaceae, the SPBal1 has a 26.03% percentage of influence although the selected principal balance that most affects this family is SPBal2 with an influence of 53.69%. It would be a mistake to think that SPBal1 has more influence on f_Ruminococcaceae than on f_Erysipelotrichaceae because the absolute value of the estimated parameter in f_Ruminococcaceae is greater than the value in f_Erysipelotrichaceae; in fact, the percentage of influence of SPBal1 on f_Ruminococcaceae is 15.42%, which is smaller than the percentage of influence of the selected principal balance on f_Erysipelotrichaceae.In TableS2 (see supplementary material), we present the estimated parameters obtained with the female dataset at genus level as it is the publicly available dataset with most parameters. Focusing on intercepts, we observe that g_Bacteroides and Other are the two genera that most influence the bacterial community, and g_Dorea is the genus with the smallest intercept. Studying the percentage of influence that each selected principal balance has on each genus, we can observe in Table S2 that SPBal5 is the one that most influences eight genera, and SPBal4 is the one that most influences six genera. This reflects the importance of SPBal5 and SPBal4, which are the balances with the least percentage of variance. Consequently, taking into account the principal balances that accumulate at least 80% of the variability is paramount in order to collect enough information. In addition, in Table S3, we can observe the estimated parameters obtained by analysing the male dataset, where both the selected principal balances with the lowest percentage of variance are also essential for describing the microbiome dynamics of the male subject.


### 5.3. Credible Intervals

The credible intervals present in Tables3, S2, and S3 (see supplementary material) are noteworthy because most of them contain zero. In order to study this fact, we analyse the percentage of zeros included at the credible intervals and the number of parameters present in each equation of the model for all databases. We scrutinise these two aspects with the proposed model and the four approaches used for comparison; see details in Section 5.5.Table4 collects this analysis and points out two overall trends: (a) as the number of microbial entities present in the data grows, the number of parameters that contain zero in his credible interval rises; (b) the model with the fewest parameters in each equation of the model is the one with the fewest parameters containing the zero in its credible interval.Table 4
Percentage of credible intervals that do not contain zero and the number of parameters in each equation of the model for all databases.S1, S2, and S3 denote the first, second, and third simulations. Note that the AR(1)-clr model does not achieve convergence in the F. genus dataset.
ModelInformationF. familiesD=9F. genusD=18MaleD=9S1 D=15S2 D=20S3 D=60AR (1)-alr% zeros12.5%6.86%0%3.80%1.57%0%Param/equation9189152060alr-alr-Bal% zeros41.66%39.21%4.16%11.90%12.28%0%Param/equation333333PM% zeros9.54%8.33%4.44%7.69%0.1%0%Param/equation665131315AR (1)-clr% zeros2.22%0%4.16%4.76%0%Param/equation101910162161alr-Bal% zeros43.75%35.29%31.25%14.28%15.78%0%Param/equation222222In all models, the expected value of each microbial entity depends on its corresponding autoregressive expression. In Table4, we have denoted the number of parameters in each autoregressive expression as Param/equation. As the autoregressive expressions have the same number of parameters for all bacteria, the estimate parameters of the less abundant microbial entities must necessarily be small. In conclusion, small families will inevitably have small parameters, and as a result of working with small parameters, the credible intervals contain zero. As we are working with compositional data, the more microbial entities we have, the smaller these microbial entities are. As a result, the number of small parameters increases in datasets with a lot of microbial entities. To sum up, a zero present in a credible interval does not mean that the parameter is irrelevant; the zeros appear due to working with compositional data.


### 5.4. Graphic Results and Modeling Datasets with a Large Number of Microbial Entities

Graphic results show the ability of our model to both explain and predict the dynamics of microbiome communities. In order to predict, we estimate the model, and we denote asaijk the value of the parameter aij at the iteration k of the Markov Chain, with k=1,…,J. In order to obtain predictions using the model, we apply the following reasoning:(8)y^itk=αit^k∑i=1Dαit^kwhereαit^k=expai0k+∑j=0Maijk⋅SPBaljt−1.Finally,y^it=meany^it1,…,y^itK denotes the predictions obtained with the model. In Figures 2 and 3 we can observe how many time points we have used to estimate the model and to make predictions for each dataset.Figure 2
Model estimation and prediction at family taxonomic level. Black points represent the data, the black line represents the expected value obtained with the proposed model, and dashed lines represent the value ofEyit±2Varyit. The vertical black line represents when the estimation of the model finishes and the prediction starts. (a) Female: the families of bacteria are, from left to right, o_bacteroidales;f_, other, f_erysipelotrichaceae, o_clostridiales;f_, f_lachnospiraceae, f_ruminococcaceae, f_veillonellaceae, f_desulfovibrionaceae, and f_bacteroidaceae. (b) Male: the families of bacteria are, from left to right, f_alcaligenaceae, f_bacteroidaceae, f_lachnospiraceae, f_porphyromonadaceae, other, f_rikenellaceae, f_ruminococcaceae, f_veillonellaceae, and f_prevotellaceae.
(a)(b)Figure 3
Model estimation and prediction of the female dataset at genus taxonomic level. Black points represent the data, the black line represents the expected value obtained with the proposed model, and dashed lines represent the value ofEyit±2Varyit. The vertical black line represents when the estimation of the model finishes and the prediction starts. The genera are, from left to right, o_bacteroidales;f_;g_, other, o_clostridiales;f_;g_, f_lachnospiraceae;g_, g_blautia, f_lachnospiraceae;g_clostridium, g_coprococcus, g_dorea, g_eubacterium, g_roseburia, g_ruminococcus, f_ruminococcaceae;g_, f_ruminococcaceae;g_clostridium, g_phascolarctobacterium, g_desulfovibrio, f_erysipelotrichaceae;g_clostridium, g_akkermansia, and g_bacteroides.First, we analyse the graphic results obtained with the female dataset at family taxonomic level because it is straightforwardly comparable with the male dataset and the genus study, as mentioned in Section5.1. Figure 2(a) indicates that the two families with the highest intercept values (f_Bacteroidaceae and Other, mentioned in Section 5.2) are also the families with the most relative abundance and variability. What is more, o_Clostridiales, f_Desulfovibrionaceae, and f_Erysipelotrichaceae, which are the families with the lowest intercepts, are the families with the lowest abundance. In order to measure the goodness of fit of the estimating and predicting part, we can observe the RMSD presented in Table 5.Table 5
Deviance information criterion (DIC) defined in [63] and the DIC defined in [64]. The root mean square deviation (RMSD) of the results obtained fitting the model and the RMSD of the results obtained using the model for prediction.
DataModelD¯pD (2002)DIC (2002)pD (2004)DIC (2004)RMSD (fitting)RMSD (predicting)Female familiesD=9PM−1.87⋅1034.54⋅101−1.82⋅1036.44⋅101−1.80⋅1030.02640.0498AR(1)-alr2.98⋅1042.09⋅1045.07⋅1049.36⋅1089.37⋅1080.02450.0398AR(1)-clr3.63⋅1042.34⋅1045.98⋅1041.03⋅1091.03⋅1090.02450.0395alr-Bal5.81⋅1044.48⋅1041.02⋅1053.55⋅1093.55⋅1090.03460.0414alr-alr-Bal4.95⋅1043.25⋅1048.20⋅1042.15⋅1092.15⋅1090.03150.0399Female generaD=18PM−5.25⋅1038.96⋅101−5.16⋅1039.57⋅101−5.15⋅1030.01740.0280AR(1)-alr4.39⋅1052.80⋅1057.19⋅1053.37⋅10103.37⋅10100.01720.0368alr-Bal8.31⋅1052.20⋅1051.05⋅1067.49⋅10107.49⋅10100.02710.0216alr-alr-Bal4.38⋅1051.41⋅1055.79⋅1052.49⋅10102.49⋅10100.02080.0207MaleD=9PM−5.19⋅1023.34⋅101−4.86⋅1025.77⋅101−4.61⋅1020.01110.0458AR(1)-alr2.96⋅1033.42⋅1036.38⋅1032.72⋅1072.72⋅1070.00980.0278AR(1)-clr6.09⋅1036.49⋅1031.25⋅1046.67⋅1076.67⋅1070.00890.0289alr-Bal1.22⋅1031.37⋅1032.60⋅1036.60⋅1066.60⋅1060.01440.0309alr-alr-Bal1.56⋅1031.70⋅1033.26⋅1038.43⋅1068.43⋅1060.01420.0306Simulated dataD=15PM−1.40⋅1031.23⋅102−1.28⋅1031.93⋅102−1.21⋅1030.02830.0401AR(1)-alr−8.59⋅1027.81⋅102−7.74⋅1015.40⋅1045.32⋅1040.02360.0398AR(1)-clr−3.67⋅1001.24⋅1038.77⋅1001.02⋅1061.02⋅1060.02390.0413alr-Bal−1.22⋅1036.94⋅101−1.15⋅1032.19⋅1039.72⋅1020.03390.0382alr-alr-Bal−1.20⋅1031.07⋅102−1.10⋅1033.13⋅1031.92⋅1030.03290.0385Simulated dataD=20PM−2.20⋅1032.16⋅102−1.98⋅1033.74⋅102−1.83⋅1030.01900.0310AR(1)-alr−1.11⋅1031.63⋅1035.13⋅1022.02⋅1052.01⋅1050.01440.0327AR(1)-clr6.63⋅1003.38⋅1034.05⋅1031.55⋅1071.55⋅1070.01470.0331alr-Bal−2.02⋅1037.67⋅101−1.94⋅1031.22⋅103−7.96⋅1020.02570.0308alr-alr-Bal−2.00⋅1031.10⋅102−1.89⋅1031.56⋅103−4.43⋅1020.02520.0309Simulated dataD=60PM−3.50⋅1033.51⋅102−3.15⋅1031.60⋅1041.25⋅1040.00010.0143AR(1)-alr5.12⋅1055.16⋅1051.02⋅1061.92⋅10111.92⋅10110.00090.0142AR(1)-clr4.91⋅1054.95⋅1059.86⋅1051.89⋅10111.89⋅10110.00090.0124alr-Bal−2.24⋅1034.86⋅102−1.75⋅1031.49⋅1051.47⋅1050.00520.0060alr-alr-Bal−1.89⋅1038.64⋅102−1.03⋅1037.79⋅1057.77⋅1050.00500.0061In Figure2 and Table 5, we can compare the results of the female and male datasets. Analysing Figure 2, we conclude that the families of the male have less variability than the families of the female. This leads to a decrease in the RMSD value obtained with the male dataset not only in the fitting part but also in the predicting part. Note that this drop is illustrated at Table 5.Hitherto, we have analysed the datasets at family taxonomic level. In Figure2, we can see that both female and male datasets have D=9 families. When we increase the number D, we increase the number of equation (7) that our model has (see Section 3.2). In order to test our model with a higher value of D, we analyse the female dataset at genus taxonomic level. In this case, D=18 and the dimension of the parameter estimation matrix aij is 18x6, as we can observe at Table S2. We must mention that the dimensions of aij with the female and male dataset at family taxonomic level are 9x6 and 9x5, respectively (see Tables 5 and S3).In Figure3 and Table 5, we can observe the results of the estimation and prediction when analysing the data at the genus taxonomic level. Figure 3 shows that g_Bacteroides, g_Akkermansia, and Other are the genera with the most variability. Indeed, we can observe in Table S2 that g_Bacteroides and Other are the two genera with the highest intercept values, and g_Akkermansia is especially influenced by SPBal1 and SPBal4.Using simulated data, we have also managed to test our model whenD=20 and D=60. In Table 4, we can see the results obtained, and we can conclude that the model achieves modeling and predicting in short-term microbiome dynamics with higher values of D. In fact, Table 5 shows that the RMSD of the fitting part of the data from the third simulation (with D=60) is smaller than the RMSD of the data from the second simulation (with D=20).


### 5.5. Comparison

For comparative purposes, we also show the results obtained with alternative models that differ from the proposed model due to the autoregressive component and the reparameterization of the Dirichlet parameters.Reference [46] highlights the ability of the Dirichlet distribution with an autoregressive structure that uses the alr or the clr transformation as a link function and imposes a parameter that is equal to the sum of the Dirichlet parameters. On the one hand, we proposed the AR(1)-alr and AR(1)-clr models that take this construction and use it in a Bayesian framework while writing in the autoregressive structured an AR expression with the alr or clr transformation, respectively. On the other hand, [41] emphasizes the importance of balances to identify microbial signatures; for this reason, we present alr-Bal and alr-alr-Bal models, which are two alternative that maintain the general construction of [46] but introduce balances in the autoregressive structure and use it in a Bayesian framework. The following alternatives have been implemented: a modification of the alr-alr-Bal model using a balance whose numerator does not include the microbial entity i, the alr-alr-Bal model but using the clr transformation to accompany the ai1 parameter and a balance that compares all microbial entities except i with all the microbial entities, and an alteration of this last model applying the clr to transform the Dirichlet parameters. However, these alternatives did not lead us to obtain relevant fitting results, and consequently, they are not included in this article.


#### 5.5.1. Comparison with AR(1)-alr Model

In this case, in order to linkαt=α1t,α2t…,αDt with yt−1,yt−2,yt−3,…,y1, we propose:(9)logαitαDt≔μit,with(10)μit=ai0+∑j=1D−1aij⋅logyjt−1yDt−1wherei=1,…,D−1.Note thatτt=α1t+α2t+…+αDt is the concentration parameter of the Dirichlet distribution and must also be estimated.


#### 5.5.2. Comparison with AR(1)-clr Model

This consideration is similar to the AR(1)-alr model. In this case, the centered log-ratioclr transformation is used. Note that this transformation scales each component vector by the geometric mean, g⋅.(11)logαitgαt:=μit,μit=ai0+∑j=1D−1aij⋅logyjt−1gyt−1wherei=1,…,D.In this alternative formulation,τt must also be estimated.


#### 5.5.3. Comparison with alr-Bal Model

This model has the same structure as the AR(1)-alr model, but in this case:(12)μit=ai0+ai1⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.


#### 5.5.4. Comparison with alr-alr-Bal Model

Following the structure of AR(1)-alr, this model presents an alternative formulation of the autoregressive component, combining the alr transformation and the balance used in the alr-Bal model.(13)μit=ai0+ai1⋅logyi,t−1yD,t−1+ai2⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.We compare the proposed model with these four alternatives to assess the ability of principal balances compared to AR models or other balance options. In fact, analysing Table5, we can thoroughly compare the five methods. The proposed model is the one that achieves the lowest DIC values not only at the 2004 DIC but also at the 2002 DIC. If we focus on the 2004 DIC, we observe that the proposed model has by far the smallest DIC value; in addition, alr-Bal and alr-alr-Bal models obtain similar results, but they enhance the values obtained with the AR methods (except in female datasets). In spite of the fact that AR methods describe similar behaviours, in general, we can observe an improvement in the AR(1)-alr model. This dynamic is also matched at the 2002 DIC values; however, the values of the 2002 DIC are smaller than the values of the 2004 DIC in all models.We can find the smallest 2002 DIC and 2004 DIC in the female dataset analysed at genus taxonomic level. The proposed model minimizes the 2002 and 2004 DICs with the values of−5.16⋅103 and −5.15⋅103, respectively. On the other hand, we observe the highest values of both DICs in the simulation with D=60. The AR(1)-clr model maximizes the 2002 DIC with a value of 9.86⋅105, and AR(1)-alr maximizes the 2004 DIC with the value of 1.92⋅1011. Note that in the female dataset analysed at genus taxonomic level, we can find the biggest difference between the DIC values obtained with the proposed model and the DICs achieved with other datasets. We must mention that the AR(1)-alr model did not achieve convergence with the dataset at genus taxonomic level. In conclusion, the difference between the proposed model and the other models is reflected far more clearly in this dataset.In Table5, we also can study the RMSD values of the fitting and the predicting parts. Notice that the values of the fitting part are smaller than the values of the predicting part in all models. Focusing on the fitting part, we can observe that alr-Bal and alr-alr-Bal have larger RMSE values than the proposed model in all datasets; however, AR(1)-alr and AR(1)-clr have in general smaller values than the proposed model. Despite this, the proposed model has advantages over these models: it reduces the number of parameters and this enhances the estimation of large datasets. This is reflected in the simulated dataset with D=60 where the proposed model has the smallest RMSE value, improving on the AR(1)-alr and AR(1)-clr RMSE values. Shifting the focus from fitting to predicting, we can observe that the AR(1)-alr or the AR(1)-clr has a larger RMSE value in all the datasets with D>9 except for the simulated dataset with D=60, where the difference between the RMSE values of the proposed model and the AR(1)-alr or AR(1)-clr is small. As a result, in this dataset, the models with AR structure and the proposed model have a similar prediction capacity; however, the proposed model describes the bacterial dynamics better, as we have seen before with the RMSE values of the fitting part. The same is true of the alr-Bal and alr-alr-Bal models; they have a smaller RMSE value in the predicting part than the proposed model, but they have the worst RSME value in the fitting part, so these models are the ones that worst describe the bacterial dynamics.In conclusion, the proposed model is the one with the smallest DIC values, and this leads to us conclude that the proposed model (PM) is the most suitable approach to model microbiome dynamics. We must take into account that the proposed model does not have the smallest RMSE value in general; in these cases, however, the proposed model still has an advantage over the other models: it better describes the bacterial dynamics of large datasets.


## 5.1. Applying Ward’s Clustering Method and Selecting the Principal Balances

In order to introduce the results in an accessible way, we first present the outcomes obtained with the publicly available datasets. Moreover, we start by analysing the female dataset at family level, because it has aspects in common with the other two datasets and, consequently, it is easily comparable.By applying Ward’s clustering method as it is described in Section 3.1, we obtain the dendrogram shown in Figure1. As it is presented by [62], the horizontal bars of the dendrogram describe the groups of parts defined at each iteration of the process; for this reason, the groups merged in the dendrogram give us the numerator and the denominator of the principal balances. Moreover, the sum of all vertical bars represents the total variance of the sample; therefore, a long vertical bar implies a balance that explains a good deal of the cc variance and a short vertical bar means that balance has a little variability in the sample. In Figure 1, we can also see the principal balances described by the dendogram and their percentage of variance (Table1).Figure 1
Female dataset at family taxonomic level. Dendogram obtained using Ward’s method for obtaining the principal balances. The table shows the expression of the principal balances and their percentage of variance.Table 1
The expression of the principal balances and their percentage of variance.
PBalExpression% varPBal11⋅1/1+1⋅lnx3/x69.00PBal22⋅1/2+1⋅lnx3⋅x61/2/x715.38PBal31⋅1/1+1⋅lnx5/x92.15PBal42⋅1/2+1⋅lnx5⋅x91/2/x47.62PBal53⋅1/3+1⋅lnx5⋅x9⋅x41/3/x110.62PBal61⋅1/1+1⋅lnx2/x85.31PBal74⋅2/4+2⋅lnx5⋅x9⋅x4⋅x11/4/x2⋅x81/214.08PBal86⋅3/6+3⋅lnx5⋅x9⋅x4⋅x1⋅x2⋅x81/6/x7⋅x3⋅x61/335.81As explained in Section3.1, we select the principal balances that accounted for at least 80% of the variability. Consequently, we pick PBal8, PBal2, PBal7, PBal5, and PBal1 because they explain a cumulative variance of 84.90% (see Table 2). Note that we have written the principal balances in descending order by percentage of variance and we denote them as SPBal1,SPBal2…SPBal5, respectively. For comparative purposes, we have estimated the proposed model using different selected principal balances to compare the results when extracting different percentages of variance of the data. In Table 2, we can observe that the formulation with SPBal1,SPBal2…SPBal5 is the one with smallest value of 2002 DIC, 2004 DIC, and RMSD when fitting; consequently, this formulation is the most suitable option. In fact, it is underscored in Section 5.2 that SPBal5 is particularly interesting for explaining the dynamics of the microbial community.Table 2
Female dataset at family taxonomic level. Results obtained when estimating the proposed model using different selected principal balances. Additionally, the table shows the deviance information criterion (DIC) defined in [63], the DIC defined in [64], the root mean square deviation (RMSD) of the results obtained fitting the model, and the RMSD of the results obtained using the model for prediction.
Proposed model different SPBal% varD¯pD (2002)DIC (2002)pD (2004)DIC (2004)RMSD (fitting)RMSD (predicting)SPBal135.81−1820.7616.93−1803.8317.60−1803.150.02990.0408SPBal1, SPBal251.19−1827.1124.16−1802.9431.55−1795.560.02920.0409SPBal1, SPBal2, SPBal365.27−1826.6131.08−1795.5240.64−1785.960.02880.0442SPBal1,SPBal2,…,SPBal475.89−1832.6537.68−1794.9645.90−1786.740.02820.0452SPBal1,SPBal2,…,SPBal584.90−1872.2345.48−1826.7564.49−1807.740.02640.0498Table 3
Estimated parametersaij and credible interval obtained by applying the proposed model to the female dataset at family taxonomic level. The table also contains the percentage of influence that each SPBalj has on each family and the percentage of sway that each family has on the bacterial community is represented in bold.
Family taxonomic levelaijInterceptSPBal1SPBal2SPBal3SPBal4SPBal5o_bacteroidales; f_−0.03660.42550.16660.0145−0.3427−0.3303(−1.06, 0.96)(0.13, 0.69)(−0.11, 0.46)(−0.33,0.37)(−0.76, 0.04)(−0.69, 0.009)0.58%33.25%13.02%1.13%26.78%25.81%Other1.10760.7105−0.07530.09130.2463−0.1231(0.25, 1.95)(0.44, 0.96)(−0.31, 0.13)(−0.20, 0.39)(−0.04, 0.53)(−0.47, 0.18)17.50%57.00%6.04%7.32%19.76%9.88%f_erysipelotrichaceae−0.5757−0.04970.1025−0.01910.01410.0055(−2.05, 0.61)(−0.40, 0.25)(−0.32, 0.54)(−0.50, 0.47)(−0.53, 0.56)(−0.50, 0.46)9.10%26.03%53.69%10.01%7.39%2.88%o_clostridiales; f_−0.87050.1802−0.03460.0038−0.02810.0198(−2.44, 0.38)(−0.15, 0.54)(−0.46, 0.35)(−0.45, 0.48)(−0.56, 0.50)(−0.50, 0.50)13.75%67.62%12.98%1.43%10.54%7.43%f_lachnospiraceae0.48470.3356−0.01990.0123−0.0225−0.5353(−0.33, 1.47)(0.04, 0.61)(−0.27, 0.21)(−0.29, 0.31)(−0.36, 0.30)(−0.88, −0.19)7,66%36,26%2,15%1,33%2,43%57,83%f_ruminococcaceae0.2661−0.20230.1230−0.21090.0594−0.7166(−0.71, 1.44)(−0.52, 0.06)(−0.19, 0.44)(−0.56, 0.10)(−0.30, 0.43)(−1.11, −0.30)4.20%15.42%9.37%16.07%4.53%54.61%f_veillonellaceae−0.20570.2451−0.1049−0.0496−0.0238−0.3804(−1.26, 0.75)(−0.05, 0.54)(−0.45, 0.16)(−0.42, 0.30)(−0.41, 0.34)(−0.76, −0.02)3,25%30,49%13,05%6,17%2,96%47,33%f_desulfovibrionaceae−0.63600.29930.03630.0866−0.1100−0.0535(−2.03, 0.45)(−0.03, 0.63)(−0.34, 0.40)(−0.32, 0.55)(−0.62, 0.36)(−0.52, 0.35)10,05%51,10%6,20%14,79%18,78%9,13%f_bacteroidaceae2.14600.2476−0.0018−0.1890−0.0426−0.5470(1.30, 3.04)(−0.008, 0.49)(−0.25, 0.22)(−0.47, 0.08)(−0.33, 0.25)(−0.87, −0.24)33,91%24,090,1818,394,1453,21


## 5.2. Estimated Parameters

Once we have selected the principal balances, we can estimate the proposed model. We are going to do so using Bayesian analysis because it allows us to simulate the posterior distribution of the model parameters using Markov Chain Monte Carlo (MCMC), as implemented in the free statistical software JAGS [65]. We fixed a burn-in period of 10000 iterations and an adaptation period of 1000 iterations. We kept one posterior sample in 100 iterations until a set of 20000 iterations was obtained, except for the first and second simulation datasets where just 2000 iterations were needed. We implemented this process using the runjags package in R [66]. Note that in Table S1 (see supplementary material), we can find the potential scale reduction factor defined by [67] and obtained using the gelman.diag function of the Coda package [68]. Table S1 shows that the potential scale reduction factor belongs to the interval 1,1.1 and is smaller than its upper 95% confidence limit.In Table3 we can find the estimated parameters, denoted in Section 3.2 as aij. Note that in each equation in the model (see equation (7)), we need to estimate six parameters, the intercept, and the five parameters associated with the selected principal balances. The estimated value of ai0 represents the abundance of the family i in the community; in fact, tiny intercepts are related to bacterial families of low abundance, whereas large intercepts indicate taxa of high abundance. In addition, the absolute value of the estimated aijwithj≠0 represents the weight that the balance SPBalj has in describing the relative abundance of the family i. As a result, we can say that the absolute value of aijwithj≠0 represents the influence of the balance SPBalj on the relative abundance of the family i.Analysing the estimated parameters allows us to extract some information about the bacterial dynamics. Focusing on intercepts described in Table3, we observe that the families o_Clostridiales; f_, f_Desulfovibrionaceae, and f_Erysipelotrichaceae have a low intercept, while the intercepts of f_Bacteroidaceae and Other are meaningful. A negative intercept indicates families of low abundance, whereas large positive intercepts indicate abundant ones. We define the percentage of sway that each family has on the bacterial community as ai0/∑i=1Dai0. This percentage of sway provides information about the abundance of a family taking into account the other families; consequently, it enables us to compare the abundance of the different families. When the intercept is positive, a high percentage indicates high abundance and vice versa. However, when the intercept is negative, high percentage is related to low abundance and vice versa. In Table 3, the percentage of sway is written in bold, and we can observe that f_Bacteroidaceae and Other, which have positive intercepts and a percentage of sway of 33.91% and 17.50%, respectively, are the two most abundant families. In addition, o_Clostridiales; f_ has a negative intercept and a percentage of 13.75%; it is therefore the least abundant family. The percentage of sway can help to gain a better understanding of the different dynamics of microbes because low percentage indicates small families, and high percentage is related to the extreme families: it indicates the most abundant if the intercept is positive or the least abundant if the percentage is negative.As we have mentioned before, the absolute value of the estimated parametersaijwithj≠0 represents the weight that the balance has in describing the relative abundance of the family, and they therefore provide information about the influence that each balance has on the family. Centering our attention on the absolute value of the estimated aij, we can observe that SPBal5 is the one that most influences f_Ruminococcaceae, f_Lachnospiraceae, f_Veillonellaceae, and f_Bacteroidaceae. Consequently, the families that belong to the SPBal5 are particularly interesting to explain the dynamics of these families. Note that SPBal5,t−1=1⋅1/1+1⋅lny3,t−1/y6,t−1; therefore, this balance explains the relationship between the components y3 and y6 (which are f_Erysipelotrichaceae and f_Ruminococcaceae respectively) at time point t−1. A balance value close to zero can be interpreted as a high similarity between the numerator elements and the denominator components; however, the larger the absolute value of the balance, the greater the dissimilarity between them. Furthermore, in Table 3, we can also observe that the coefficient of SPBal1 has a significant absolute value in all families. As a result, SPBal1, which is the selected principal balance with the highest percentage of variance, is also the one that has a significant influence on most of the bacterial families.Families with low abundances will have low parameter values; as a result, in order to compare the influence of oneSPBalj0 in different families, it is interesting to study the percentage of influence aij/∑j=1Maij. Observing Table 3, it seems that SPBal1 does not have an special influence on f_Erysipelotrichaceae because aij<0.1, while in the rest of the families aij>0.1, however, if we observe the value of all aij at f_Erysipelotrichaceae, we can observe that all the values are small and, in comparison, SPBal1 has a significant influence. This is the reason why it is interesting to study the percentage of influence present in Table 3. In f_Erysipelotrichaceae, the SPBal1 has a 26.03% percentage of influence although the selected principal balance that most affects this family is SPBal2 with an influence of 53.69%. It would be a mistake to think that SPBal1 has more influence on f_Ruminococcaceae than on f_Erysipelotrichaceae because the absolute value of the estimated parameter in f_Ruminococcaceae is greater than the value in f_Erysipelotrichaceae; in fact, the percentage of influence of SPBal1 on f_Ruminococcaceae is 15.42%, which is smaller than the percentage of influence of the selected principal balance on f_Erysipelotrichaceae.In TableS2 (see supplementary material), we present the estimated parameters obtained with the female dataset at genus level as it is the publicly available dataset with most parameters. Focusing on intercepts, we observe that g_Bacteroides and Other are the two genera that most influence the bacterial community, and g_Dorea is the genus with the smallest intercept. Studying the percentage of influence that each selected principal balance has on each genus, we can observe in Table S2 that SPBal5 is the one that most influences eight genera, and SPBal4 is the one that most influences six genera. This reflects the importance of SPBal5 and SPBal4, which are the balances with the least percentage of variance. Consequently, taking into account the principal balances that accumulate at least 80% of the variability is paramount in order to collect enough information. In addition, in Table S3, we can observe the estimated parameters obtained by analysing the male dataset, where both the selected principal balances with the lowest percentage of variance are also essential for describing the microbiome dynamics of the male subject.


## 5.3. Credible Intervals

The credible intervals present in Tables3, S2, and S3 (see supplementary material) are noteworthy because most of them contain zero. In order to study this fact, we analyse the percentage of zeros included at the credible intervals and the number of parameters present in each equation of the model for all databases. We scrutinise these two aspects with the proposed model and the four approaches used for comparison; see details in Section 5.5.Table4 collects this analysis and points out two overall trends: (a) as the number of microbial entities present in the data grows, the number of parameters that contain zero in his credible interval rises; (b) the model with the fewest parameters in each equation of the model is the one with the fewest parameters containing the zero in its credible interval.Table 4
Percentage of credible intervals that do not contain zero and the number of parameters in each equation of the model for all databases.S1, S2, and S3 denote the first, second, and third simulations. Note that the AR(1)-clr model does not achieve convergence in the F. genus dataset.
ModelInformationF. familiesD=9F. genusD=18MaleD=9S1 D=15S2 D=20S3 D=60AR (1)-alr% zeros12.5%6.86%0%3.80%1.57%0%Param/equation9189152060alr-alr-Bal% zeros41.66%39.21%4.16%11.90%12.28%0%Param/equation333333PM% zeros9.54%8.33%4.44%7.69%0.1%0%Param/equation665131315AR (1)-clr% zeros2.22%0%4.16%4.76%0%Param/equation101910162161alr-Bal% zeros43.75%35.29%31.25%14.28%15.78%0%Param/equation222222In all models, the expected value of each microbial entity depends on its corresponding autoregressive expression. In Table4, we have denoted the number of parameters in each autoregressive expression as Param/equation. As the autoregressive expressions have the same number of parameters for all bacteria, the estimate parameters of the less abundant microbial entities must necessarily be small. In conclusion, small families will inevitably have small parameters, and as a result of working with small parameters, the credible intervals contain zero. As we are working with compositional data, the more microbial entities we have, the smaller these microbial entities are. As a result, the number of small parameters increases in datasets with a lot of microbial entities. To sum up, a zero present in a credible interval does not mean that the parameter is irrelevant; the zeros appear due to working with compositional data.


## 5.4. Graphic Results and Modeling Datasets with a Large Number of Microbial Entities

Graphic results show the ability of our model to both explain and predict the dynamics of microbiome communities. In order to predict, we estimate the model, and we denote asaijk the value of the parameter aij at the iteration k of the Markov Chain, with k=1,…,J. In order to obtain predictions using the model, we apply the following reasoning:(8)y^itk=αit^k∑i=1Dαit^kwhereαit^k=expai0k+∑j=0Maijk⋅SPBaljt−1.Finally,y^it=meany^it1,…,y^itK denotes the predictions obtained with the model. In Figures 2 and 3 we can observe how many time points we have used to estimate the model and to make predictions for each dataset.Figure 2
Model estimation and prediction at family taxonomic level. Black points represent the data, the black line represents the expected value obtained with the proposed model, and dashed lines represent the value ofEyit±2Varyit. The vertical black line represents when the estimation of the model finishes and the prediction starts. (a) Female: the families of bacteria are, from left to right, o_bacteroidales;f_, other, f_erysipelotrichaceae, o_clostridiales;f_, f_lachnospiraceae, f_ruminococcaceae, f_veillonellaceae, f_desulfovibrionaceae, and f_bacteroidaceae. (b) Male: the families of bacteria are, from left to right, f_alcaligenaceae, f_bacteroidaceae, f_lachnospiraceae, f_porphyromonadaceae, other, f_rikenellaceae, f_ruminococcaceae, f_veillonellaceae, and f_prevotellaceae.
(a)(b)Figure 3
Model estimation and prediction of the female dataset at genus taxonomic level. Black points represent the data, the black line represents the expected value obtained with the proposed model, and dashed lines represent the value ofEyit±2Varyit. The vertical black line represents when the estimation of the model finishes and the prediction starts. The genera are, from left to right, o_bacteroidales;f_;g_, other, o_clostridiales;f_;g_, f_lachnospiraceae;g_, g_blautia, f_lachnospiraceae;g_clostridium, g_coprococcus, g_dorea, g_eubacterium, g_roseburia, g_ruminococcus, f_ruminococcaceae;g_, f_ruminococcaceae;g_clostridium, g_phascolarctobacterium, g_desulfovibrio, f_erysipelotrichaceae;g_clostridium, g_akkermansia, and g_bacteroides.First, we analyse the graphic results obtained with the female dataset at family taxonomic level because it is straightforwardly comparable with the male dataset and the genus study, as mentioned in Section5.1. Figure 2(a) indicates that the two families with the highest intercept values (f_Bacteroidaceae and Other, mentioned in Section 5.2) are also the families with the most relative abundance and variability. What is more, o_Clostridiales, f_Desulfovibrionaceae, and f_Erysipelotrichaceae, which are the families with the lowest intercepts, are the families with the lowest abundance. In order to measure the goodness of fit of the estimating and predicting part, we can observe the RMSD presented in Table 5.Table 5
Deviance information criterion (DIC) defined in [63] and the DIC defined in [64]. The root mean square deviation (RMSD) of the results obtained fitting the model and the RMSD of the results obtained using the model for prediction.
DataModelD¯pD (2002)DIC (2002)pD (2004)DIC (2004)RMSD (fitting)RMSD (predicting)Female familiesD=9PM−1.87⋅1034.54⋅101−1.82⋅1036.44⋅101−1.80⋅1030.02640.0498AR(1)-alr2.98⋅1042.09⋅1045.07⋅1049.36⋅1089.37⋅1080.02450.0398AR(1)-clr3.63⋅1042.34⋅1045.98⋅1041.03⋅1091.03⋅1090.02450.0395alr-Bal5.81⋅1044.48⋅1041.02⋅1053.55⋅1093.55⋅1090.03460.0414alr-alr-Bal4.95⋅1043.25⋅1048.20⋅1042.15⋅1092.15⋅1090.03150.0399Female generaD=18PM−5.25⋅1038.96⋅101−5.16⋅1039.57⋅101−5.15⋅1030.01740.0280AR(1)-alr4.39⋅1052.80⋅1057.19⋅1053.37⋅10103.37⋅10100.01720.0368alr-Bal8.31⋅1052.20⋅1051.05⋅1067.49⋅10107.49⋅10100.02710.0216alr-alr-Bal4.38⋅1051.41⋅1055.79⋅1052.49⋅10102.49⋅10100.02080.0207MaleD=9PM−5.19⋅1023.34⋅101−4.86⋅1025.77⋅101−4.61⋅1020.01110.0458AR(1)-alr2.96⋅1033.42⋅1036.38⋅1032.72⋅1072.72⋅1070.00980.0278AR(1)-clr6.09⋅1036.49⋅1031.25⋅1046.67⋅1076.67⋅1070.00890.0289alr-Bal1.22⋅1031.37⋅1032.60⋅1036.60⋅1066.60⋅1060.01440.0309alr-alr-Bal1.56⋅1031.70⋅1033.26⋅1038.43⋅1068.43⋅1060.01420.0306Simulated dataD=15PM−1.40⋅1031.23⋅102−1.28⋅1031.93⋅102−1.21⋅1030.02830.0401AR(1)-alr−8.59⋅1027.81⋅102−7.74⋅1015.40⋅1045.32⋅1040.02360.0398AR(1)-clr−3.67⋅1001.24⋅1038.77⋅1001.02⋅1061.02⋅1060.02390.0413alr-Bal−1.22⋅1036.94⋅101−1.15⋅1032.19⋅1039.72⋅1020.03390.0382alr-alr-Bal−1.20⋅1031.07⋅102−1.10⋅1033.13⋅1031.92⋅1030.03290.0385Simulated dataD=20PM−2.20⋅1032.16⋅102−1.98⋅1033.74⋅102−1.83⋅1030.01900.0310AR(1)-alr−1.11⋅1031.63⋅1035.13⋅1022.02⋅1052.01⋅1050.01440.0327AR(1)-clr6.63⋅1003.38⋅1034.05⋅1031.55⋅1071.55⋅1070.01470.0331alr-Bal−2.02⋅1037.67⋅101−1.94⋅1031.22⋅103−7.96⋅1020.02570.0308alr-alr-Bal−2.00⋅1031.10⋅102−1.89⋅1031.56⋅103−4.43⋅1020.02520.0309Simulated dataD=60PM−3.50⋅1033.51⋅102−3.15⋅1031.60⋅1041.25⋅1040.00010.0143AR(1)-alr5.12⋅1055.16⋅1051.02⋅1061.92⋅10111.92⋅10110.00090.0142AR(1)-clr4.91⋅1054.95⋅1059.86⋅1051.89⋅10111.89⋅10110.00090.0124alr-Bal−2.24⋅1034.86⋅102−1.75⋅1031.49⋅1051.47⋅1050.00520.0060alr-alr-Bal−1.89⋅1038.64⋅102−1.03⋅1037.79⋅1057.77⋅1050.00500.0061In Figure2 and Table 5, we can compare the results of the female and male datasets. Analysing Figure 2, we conclude that the families of the male have less variability than the families of the female. This leads to a decrease in the RMSD value obtained with the male dataset not only in the fitting part but also in the predicting part. Note that this drop is illustrated at Table 5.Hitherto, we have analysed the datasets at family taxonomic level. In Figure2, we can see that both female and male datasets have D=9 families. When we increase the number D, we increase the number of equation (7) that our model has (see Section 3.2). In order to test our model with a higher value of D, we analyse the female dataset at genus taxonomic level. In this case, D=18 and the dimension of the parameter estimation matrix aij is 18x6, as we can observe at Table S2. We must mention that the dimensions of aij with the female and male dataset at family taxonomic level are 9x6 and 9x5, respectively (see Tables 5 and S3).In Figure3 and Table 5, we can observe the results of the estimation and prediction when analysing the data at the genus taxonomic level. Figure 3 shows that g_Bacteroides, g_Akkermansia, and Other are the genera with the most variability. Indeed, we can observe in Table S2 that g_Bacteroides and Other are the two genera with the highest intercept values, and g_Akkermansia is especially influenced by SPBal1 and SPBal4.Using simulated data, we have also managed to test our model whenD=20 and D=60. In Table 4, we can see the results obtained, and we can conclude that the model achieves modeling and predicting in short-term microbiome dynamics with higher values of D. In fact, Table 5 shows that the RMSD of the fitting part of the data from the third simulation (with D=60) is smaller than the RMSD of the data from the second simulation (with D=20).


## 5.5. Comparison

For comparative purposes, we also show the results obtained with alternative models that differ from the proposed model due to the autoregressive component and the reparameterization of the Dirichlet parameters.Reference [46] highlights the ability of the Dirichlet distribution with an autoregressive structure that uses the alr or the clr transformation as a link function and imposes a parameter that is equal to the sum of the Dirichlet parameters. On the one hand, we proposed the AR(1)-alr and AR(1)-clr models that take this construction and use it in a Bayesian framework while writing in the autoregressive structured an AR expression with the alr or clr transformation, respectively. On the other hand, [41] emphasizes the importance of balances to identify microbial signatures; for this reason, we present alr-Bal and alr-alr-Bal models, which are two alternative that maintain the general construction of [46] but introduce balances in the autoregressive structure and use it in a Bayesian framework. The following alternatives have been implemented: a modification of the alr-alr-Bal model using a balance whose numerator does not include the microbial entity i, the alr-alr-Bal model but using the clr transformation to accompany the ai1 parameter and a balance that compares all microbial entities except i with all the microbial entities, and an alteration of this last model applying the clr to transform the Dirichlet parameters. However, these alternatives did not lead us to obtain relevant fitting results, and consequently, they are not included in this article.


### 5.5.1. Comparison with AR(1)-alr Model

In this case, in order to linkαt=α1t,α2t…,αDt with yt−1,yt−2,yt−3,…,y1, we propose:(9)logαitαDt≔μit,with(10)μit=ai0+∑j=1D−1aij⋅logyjt−1yDt−1wherei=1,…,D−1.Note thatτt=α1t+α2t+…+αDt is the concentration parameter of the Dirichlet distribution and must also be estimated.


### 5.5.2. Comparison with AR(1)-clr Model

This consideration is similar to the AR(1)-alr model. In this case, the centered log-ratioclr transformation is used. Note that this transformation scales each component vector by the geometric mean, g⋅.(11)logαitgαt:=μit,μit=ai0+∑j=1D−1aij⋅logyjt−1gyt−1wherei=1,…,D.In this alternative formulation,τt must also be estimated.


### 5.5.3. Comparison with alr-Bal Model

This model has the same structure as the AR(1)-alr model, but in this case:(12)μit=ai0+ai1⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.


### 5.5.4. Comparison with alr-alr-Bal Model

Following the structure of AR(1)-alr, this model presents an alternative formulation of the autoregressive component, combining the alr transformation and the balance used in the alr-Bal model.(13)μit=ai0+ai1⋅logyi,t−1yD,t−1+ai2⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.We compare the proposed model with these four alternatives to assess the ability of principal balances compared to AR models or other balance options. In fact, analysing Table5, we can thoroughly compare the five methods. The proposed model is the one that achieves the lowest DIC values not only at the 2004 DIC but also at the 2002 DIC. If we focus on the 2004 DIC, we observe that the proposed model has by far the smallest DIC value; in addition, alr-Bal and alr-alr-Bal models obtain similar results, but they enhance the values obtained with the AR methods (except in female datasets). In spite of the fact that AR methods describe similar behaviours, in general, we can observe an improvement in the AR(1)-alr model. This dynamic is also matched at the 2002 DIC values; however, the values of the 2002 DIC are smaller than the values of the 2004 DIC in all models.We can find the smallest 2002 DIC and 2004 DIC in the female dataset analysed at genus taxonomic level. The proposed model minimizes the 2002 and 2004 DICs with the values of−5.16⋅103 and −5.15⋅103, respectively. On the other hand, we observe the highest values of both DICs in the simulation with D=60. The AR(1)-clr model maximizes the 2002 DIC with a value of 9.86⋅105, and AR(1)-alr maximizes the 2004 DIC with the value of 1.92⋅1011. Note that in the female dataset analysed at genus taxonomic level, we can find the biggest difference between the DIC values obtained with the proposed model and the DICs achieved with other datasets. We must mention that the AR(1)-alr model did not achieve convergence with the dataset at genus taxonomic level. In conclusion, the difference between the proposed model and the other models is reflected far more clearly in this dataset.In Table5, we also can study the RMSD values of the fitting and the predicting parts. Notice that the values of the fitting part are smaller than the values of the predicting part in all models. Focusing on the fitting part, we can observe that alr-Bal and alr-alr-Bal have larger RMSE values than the proposed model in all datasets; however, AR(1)-alr and AR(1)-clr have in general smaller values than the proposed model. Despite this, the proposed model has advantages over these models: it reduces the number of parameters and this enhances the estimation of large datasets. This is reflected in the simulated dataset with D=60 where the proposed model has the smallest RMSE value, improving on the AR(1)-alr and AR(1)-clr RMSE values. Shifting the focus from fitting to predicting, we can observe that the AR(1)-alr or the AR(1)-clr has a larger RMSE value in all the datasets with D>9 except for the simulated dataset with D=60, where the difference between the RMSE values of the proposed model and the AR(1)-alr or AR(1)-clr is small. As a result, in this dataset, the models with AR structure and the proposed model have a similar prediction capacity; however, the proposed model describes the bacterial dynamics better, as we have seen before with the RMSE values of the fitting part. The same is true of the alr-Bal and alr-alr-Bal models; they have a smaller RMSE value in the predicting part than the proposed model, but they have the worst RSME value in the fitting part, so these models are the ones that worst describe the bacterial dynamics.In conclusion, the proposed model is the one with the smallest DIC values, and this leads to us conclude that the proposed model (PM) is the most suitable approach to model microbiome dynamics. We must take into account that the proposed model does not have the smallest RMSE value in general; in these cases, however, the proposed model still has an advantage over the other models: it better describes the bacterial dynamics of large datasets.


## 5.5.1. Comparison with AR(1)-alr Model

In this case, in order to linkαt=α1t,α2t…,αDt with yt−1,yt−2,yt−3,…,y1, we propose:(9)logαitαDt≔μit,with(10)μit=ai0+∑j=1D−1aij⋅logyjt−1yDt−1wherei=1,…,D−1.Note thatτt=α1t+α2t+…+αDt is the concentration parameter of the Dirichlet distribution and must also be estimated.


## 5.5.2. Comparison with AR(1)-clr Model

This consideration is similar to the AR(1)-alr model. In this case, the centered log-ratioclr transformation is used. Note that this transformation scales each component vector by the geometric mean, g⋅.(11)logαitgαt:=μit,μit=ai0+∑j=1D−1aij⋅logyjt−1gyt−1wherei=1,…,D.In this alternative formulation,τt must also be estimated.


## 5.5.3. Comparison with alr-Bal Model

This model has the same structure as the AR(1)-alr model, but in this case:(12)μit=ai0+ai1⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.


## 5.5.4. Comparison with alr-alr-Bal Model

Following the structure of AR(1)-alr, this model presents an alternative formulation of the autoregressive component, combining the alr transformation and the balance used in the alr-Bal model.(13)μit=ai0+ai1⋅logyi,t−1yD,t−1+ai2⋅loggy1,t−1,y2,t−1,…,yD−1,t−1gyD,t−1.Note thatlogg⋅/g⋅ is a balance expression and τt must also be estimated.We compare the proposed model with these four alternatives to assess the ability of principal balances compared to AR models or other balance options. In fact, analysing Table5, we can thoroughly compare the five methods. The proposed model is the one that achieves the lowest DIC values not only at the 2004 DIC but also at the 2002 DIC. If we focus on the 2004 DIC, we observe that the proposed model has by far the smallest DIC value; in addition, alr-Bal and alr-alr-Bal models obtain similar results, but they enhance the values obtained with the AR methods (except in female datasets). In spite of the fact that AR methods describe similar behaviours, in general, we can observe an improvement in the AR(1)-alr model. This dynamic is also matched at the 2002 DIC values; however, the values of the 2002 DIC are smaller than the values of the 2004 DIC in all models.We can find the smallest 2002 DIC and 2004 DIC in the female dataset analysed at genus taxonomic level. The proposed model minimizes the 2002 and 2004 DICs with the values of−5.16⋅103 and −5.15⋅103, respectively. On the other hand, we observe the highest values of both DICs in the simulation with D=60. The AR(1)-clr model maximizes the 2002 DIC with a value of 9.86⋅105, and AR(1)-alr maximizes the 2004 DIC with the value of 1.92⋅1011. Note that in the female dataset analysed at genus taxonomic level, we can find the biggest difference between the DIC values obtained with the proposed model and the DICs achieved with other datasets. We must mention that the AR(1)-alr model did not achieve convergence with the dataset at genus taxonomic level. In conclusion, the difference between the proposed model and the other models is reflected far more clearly in this dataset.In Table5, we also can study the RMSD values of the fitting and the predicting parts. Notice that the values of the fitting part are smaller than the values of the predicting part in all models. Focusing on the fitting part, we can observe that alr-Bal and alr-alr-Bal have larger RMSE values than the proposed model in all datasets; however, AR(1)-alr and AR(1)-clr have in general smaller values than the proposed model. Despite this, the proposed model has advantages over these models: it reduces the number of parameters and this enhances the estimation of large datasets. This is reflected in the simulated dataset with D=60 where the proposed model has the smallest RMSE value, improving on the AR(1)-alr and AR(1)-clr RMSE values. Shifting the focus from fitting to predicting, we can observe that the AR(1)-alr or the AR(1)-clr has a larger RMSE value in all the datasets with D>9 except for the simulated dataset with D=60, where the difference between the RMSE values of the proposed model and the AR(1)-alr or AR(1)-clr is small. As a result, in this dataset, the models with AR structure and the proposed model have a similar prediction capacity; however, the proposed model describes the bacterial dynamics better, as we have seen before with the RMSE values of the fitting part. The same is true of the alr-Bal and alr-alr-Bal models; they have a smaller RMSE value in the predicting part than the proposed model, but they have the worst RSME value in the fitting part, so these models are the ones that worst describe the bacterial dynamics.In conclusion, the proposed model is the one with the smallest DIC values, and this leads to us conclude that the proposed model (PM) is the most suitable approach to model microbiome dynamics. We must take into account that the proposed model does not have the smallest RMSE value in general; in these cases, however, the proposed model still has an advantage over the other models: it better describes the bacterial dynamics of large datasets.


## 6. Conclusion

In this paper, we develop a compositional autoregressive model in Bayesian framework that is able to describe microbiome time series dynamics. In order to model the relative abundance of microbiome longitudinal data, we use the Dirichlet distribution with time-varying parameters. We propose that the logarithm of the Dirichlet parameters can be explained by an autoregressive structure, which contains the principal balances that account for at least80% of the variability. A comparison study selecting principal balances that account for a lower percentage of variance is provided in this article. Even though we have developed the model for the analysis of microbiome longitudinal data, it can be adapted for other scenarios defined by longitudinal data with compositional nature.Principal balances are a compositional tool that enables us to extract the balances that maximize the variance. As balances with more variability will better explain the variability of the data, this formulation allows us to reduce the number of parameters to be estimated without removing or grouping microbial taxa. A direct consequence of reducing the number of parameters is being able to model datasets with a higher number of microbial taxa. In fact, we demonstrated the ability of our model with six datasets, three of which are publicly available and present the data of a healthy female and a healthy male, while the rest are simulated. Two of the publicly available datasets have 9 families to model and the other has 18 genera. The simulated datasets have 15, 20, and 60 microbial taxa to model, respectively.In conclusion, we have tested the ability of our model with datasets that have different numbers of microbial taxa, different taxonomic levels, healthy people of different sexes, and datasets that are publicly available and simulated. This paper presents the results obtained by analysing these datasets using the proposed model and four alternative models with which we compare the proposed model. These alternative models differ from the proposed model because they change the reparameterization of the Dirichlet parameters using the alr and the clr transformations, and they modify the autoregressive component using an AR model structure or other balance options. After comparing them, we conclude that models with balances are the ones that worst describe the bacterial dynamics but obtain the best prediction results. Models with AR structure describe the data better than the proposed model, except for big datasets, where they also present the worst prediction results. Finally, the proposed model better describes the bacterial dynamics of big datasets due to reducing the number of parameters and obtains the lowest DIC values. In addition, it can be effective for short-term prediction of the future dynamics of a microbial community. This approach can be an interesting alternative to predict the future dynamics of a microbial community.Additionally, the proposed model allows us to analyse the microbial interactions using the value of the estimated parameters. Indeed, the estimated value of the intercept represents the abundance of the microbial taxa in the community; as a matter of fact, small intercepts are related to bacterial families of low abundance, whereas large intercepts indicate taxa of high abundance. In addition, the absolute value of the rest of the estimated parameters represents the weight that the selected principal balance has in describing the relative abundance of the family; as a result, we can say that the absolute value ofaijwithj≠0 represents the influence of the balance SPBalj at the relative abundance of the family i. We must mention that balances provide information about the similarity between their numerator and denominator components; in fact, balances that are near zero indicate that the two groups are similar, and the higher the absolute value of the balance, the larger is the dissimilarity between them.Nevertheless, our proposal also has some limitations. Even if we choose the method for obtaining the principal balances that requires least computation time, estimating the models in big datasets still requires a significant amount of time. As future development, in order to improve the usefulness of the model, we propose to extend it to consider external covariates. The incorporation of covariates affecting microbiota dynamics, such as antibiotics or diet variations, may enable us to study the effectiveness of some medical treatments. Moreover, including random effects will allow us to study different subjects at the same time, such as healthy and unhealthy individuals, for example.



---

*Source: 4907527-2022-08-30.xml*