Efficient Bayesian Additive Regression Models For Microbiome and Gene Expression Studies

Venue: NeurIPS_cc_2024_Workshop_BDU
Published: 2024-08-30 22:32:13

Abstract: Analyzing sequence count data, such as microbiome or gene expression data, poses a challenge. Researchers often want to estimate linear and non-linear effects of covariates on microbial or gene composition. Bayesian multinomial logistic-normal (MLN) models have gained popularity due to their ability to account for the count compositional nature of these data. However, inferring MLN models can be challenging. Recently, we developed a computationally efficient and accurate approach to inferring MLN models with a Marginally Latent Matrix-T Process (MLTP) form. Our approach is based on a particle filter with marginal Laplace approximation -- called the \textit{Collapse-Uncollapse} (CU) sampler. Here, we introduce a new class of MLN Additive Gaussian Process models (\textit{MultiAddGPs}) for additive deconvolution of overlapping linear and non-linear effects. MultiAddGPs can be efficiently inferred using the CU sampler. Furthermore, we develop an efficient approach to hyperparameter selection via Maximum Marginal Likelihood estimation. We validate our approach via both simulated and real data studies. The code is now available: https://cran.r-project.org/web/packages/fido/index.html

Main paper: 
 introduction: 
 High-throughput sequencing (HTS) technologies in microbial ecology generate multivariate count data to characterize and analyze microbial communities from a variety of habitats, such as human body sites, soil, and water. These studies often involve multiple individuals resulting in a count table where each column is a vector of counts repressing the number of observed microbes in each of D different classes (e.g., genera).Host-associated (e.g., human and animal) microbiota are dynamic: community composition can change on daily time scales in response to various factors including diet, disease, and drugs ( [3] , [11] ). Disentangling the impact of host factors on microbiota is therefore complicated by multiple overlapping factors that can lead to linear and non-linear changes in composition over time.Existing methods for analyzing this data are not suitable as they either ignore the dual count and compositional nature of these data by assuming the data is Gaussian distributed, or they struggle with scalability in larger sample sizes and diverse taxa. (e.g., [1] , [12] ). More recently, we showed that a wide variety of Bayesian MLN models including generalized linear models and generalized Gaussian process regression models share a common marginal form called a Latent Matrix-T Process (LTP) [13] . We proved that a Laplace approximation to this marginal form was extremely accurate leading to an efficient and accurate approximate inference procedure called the Collapse-Uncollapse (CU) sampler. We showed this approach is often 4-5 orders of magnitude faster than HMC-based methods with little to no error in posterior calculations [13] .Despite these advances, there is still a lack of tools to disentangle the effects of multiple factors on microbiota. Our previous work with Marginal Latent Matrix-T Process (MLTPs) was limited to factors with linear effects on microbial composition (generalized linear models) or a single factor with a non-linear effect (generalized Gaussian process regression models). Additionally, we did not address hyperparameter selection, which is crucial in both linear and non-linear modeling (e.g., kernel parameter selection). Although alternative frameworks maybe more flexible (e.g., 2), they assume the data is transformed Gaussian, ignoring some of the critical complexities of the measurement process. This article addresses the limitations of prior methods and develops a flexible, and computationally efficient approach to disentangling both linear and non-linear effects on microbiota. As in [2] , our approach is based on a class of additive Gaussian process regression models. Unlike [2] , we do not assume the data is transformed Gaussian and instead prove that Bayesian MLN Additive Gaussian Process Models are also MLTPs. Using those results, we extend the CU sampler to this class of models. Moreover, we extend the results of [13] by developing an efficient maximum marginal likelihood estimator for model hyperparameters. Finally, we validate our approach through both simulated and real data studies. Through those studies we find our methods improve inferential performance: posterior intervals from our models often cover the truth 10 to even 100 times more often than those of transformed Gaussian models. Motivating by problems in microbiota analysis, this work provides the first flexible and computationally efficient approach to disentangling linear and non-linear processes using Bayesian MLN models. 
methods: 
 To facilitate additive linear and nonlinear modeling within a Bayesian MLN framework, this article introduces Multinomial Logistic Normal Additive Gaussian Process Models (MultiAddGPs).In this section, we first present the model and demonstrate that MultiAddGPs are a specific type of MLTP model (a comprehensive review of MLTP models and Collapse-Uncollapsed (CU) sampler are provided in Appendix A). Following this, we outline the process for conducting posterior inference using an extended CU sampler in MultiAddGPs. Additionally, we propose a computationally efficient approach for estimating hyperparameters in MLTP models using maximum marginal likelihood.Let Y •n denote a D vector of observed data, X •n a Q 0 vector of covariates to model linearly, and each Z (k∈{1,...,K})•n a Q k vector of covariates to be modeled with distinct non-linear functions. The MultiAddGPs models have the following form:EQUATIONEQUATIONEQUATIONwith priorsB ∼ N (Θ (0) , Σ, Γ (0) ), f (k) ∼ GP(Θ (k) , Σ, Γ (k) ), and Σ ∼ InvWishart(Ξ, ζ). As in the Appendix A, ϕ denotes any log-ratio transform from S D to R D-1 . Σ is a D -1 × D -1covariance matrix. For the matrix-normal prior on the linear term, Θ (0) is the mean matrix andΓ (0) is a Q 0 ×Q 0 covariance matrix representing covariance in the parameters of the Q 0 covariates.The terms Θ (k) and Γ (k) in the K matrix-normal process priors echo their linear counterparts but are functions (e.g., mean and kernel functions) rather than fixed dimensional matrices.Additionally, we permit hyperparameters Ω in any part of the model including but not limited to hyperparameters in the kernel functions (e.g., bandwidth parameters) or hyperparameters in the mean functions. As we will show through simulated and real data analyses in Appendix B and Appendix C, this is a very flexible form of model that can be used in a wide range of additive linear and non-linear modeling tasks.Posterior Estimation in MultiAddGPs. We use MLTP theory to sample from the posterior of MultiAddGPs: p(H, B, f (1) , . . . , f (K) , Σ | Y). To connect to that theory we first prove that MultiAddGPs are Bayesian Multinomial Logistic Normal MLTPs (MLN-MLTPs) with Φ = {B, f (1) , . . . , given in [13] . We obtain posterior samples of the collapsed form using a Laplace approximation as described in Appendix A.f (K) , Σ}, M = Θ (0) X + K k=1 Θ (k) (Z (k) ); V = Ξ; A = X T Γ (0) X + K k=1 Γ (k) (Z (k) ) + I; ν = ζ (Up to this point we have not distinguished between the set of points n ∈ {1, . . . , N } at which we have observed data Y and the potentially different set n * ∈ {1, . . . , N * } at which we want to predict the functions F, f 1 , . . . , f K . In what follows we use the symbols F, f 1 , . . . , and f K to denote the evaluation of corresponding infinite dimensional functions at the set of evaluation points {1, . . . , N * }, i.e., they are each D × N * -dimensional random matrices. In contrast all other random matricies (e.g., H) represent their corresponding infinite dimensional analogues evaluated at the set of observed points ({1, . . . , N }).Sampling from the uncollapsed form starts by obtaining samples from p(F, Σ | H). To sample from this distribution note that, conditioning on H and marginalizing over f (1) , . . . , f (k) in the MultiAddGP model, we obtain matrix-normal regression model:H •n ∼ N (F •n , Σ)with prior: ζ ). This is the same matrix-normal process regression model discussed in [13] : samples from p(F, Σ | H) can be obtained via methods described in Appendix C of that article.F ∼ GP Θ (0) X + K k=1 Θ (k) (Z (k) ), Σ, X T Γ (0) X + K k=1 Γ (k) (Z (k) ) Σ ∼ InvWishart(Ξ,Finally, conditioned on samples of F and Σ, we can obtain samples of each f (k) . Inspired by the backfitting algorithm used for estimation in generalized additive models [5] , we developed a backsampling algorithm for this task. The intuition of backsampling algorithm is iteratively sam-pling from p(B | F), p(f 1 | F, B), p(f 2 | F, B, f 1 ), . . ., and then p(f (K) | F, B, f 1 , . . . , f (K-1) ).For example, to sample p(f (k) | F, B, f (1) , . . . , f (k-1) ) we can use the following result. Letf * = F -BX -f (1) -• • • -f (k-1)then we have the following model for f * : (k) where Γ * = K j=k+1 Γ (j) . This model has a closed-form conditional:f * ∼ N f (k) , Σ, Γ * f (k) ∼ N Θ (k) , Σ, Γf (k) | Σ, f * ∼ N f * Σ - * + Θ (k) Γ -(k) Γ - * + Γ -(k) -1 , Σ, Γ - * + Γ -(k) -1where Γ - * and Γ -(k) are short-hand for (Γ * ) -1 and (Γ (k) ) -1 respectively. We describe the full procedure of the algorithm and provide pseudo-code in the Appendix D and E.Model Identification. Identifiability is a common problem in function decomposition models such as generalized additive models [5] . This problem is often addressed by imposing sum-to-zero constraints on the functions (e.g., f (k) dZ (k) = 0) [5] , yet others have proposed modifying kernel functions to impose identification [7] . For simplicity, this article uses the former for simulated and real data analyses: posterior samples of each f (k) are centered as f (k) -1 • mean(f (k)A first-order Laplace approximation to this integral is given by: where ĤΩ denotes the MAP estimate for H H H evaluated at the observed set {1, . . . , N } which depends on the hyperparameters Ω. While the MML estimate already provides a degree of penalization [16] , when specifying weakly identified MultiAddGPs it can be helpful to add an additional penalty term (λ) to the marginal likelihood. As we demonstrate in Appendix C and F, we often choose this penalty term as the log density of a prior for the hyperparameters. Due to MML is often non-convex especially when estimating bandwidth parameters in kernel functions, we generally recommend Bayesian optimization procedures. These procedures are well-suited for optimizing non-convex functions and can reduce the number of function evaluations compared to more typical zero-order optimizers (e.g., the L-BFGS optimizer).log p(H H H, Y Y Y | Ω)dH H H ≈ (D -1)N 2 log(2π) + log p( ĤΩ , Y | Ω) - 1 2 log(-|∇ 2 [vec( ĤΩ )]|) (5)Here, we demonstrate how to efficiently compute and sample from the conditional posterior p(F, B, f (1) , . . . , f (K) , Σ|H, Y, X, Z). Since F and Σ are conditionally independent of Y given H, and B, f (1) , . . . , f (K) are conditionally independent of H given F, by applying the chain rule, we can rewrite the equation as:p(F, B, f (1) , . . . , f (K) , Σ|H, Y, X, Z) = p(B, f (1) , . . . , f (K) |F, Σ, X, Z)p(F|Σ, H, X, Z)p(Σ|H, X, Z)The second and third parts on the right-hand side of the equation represent the posterior of a multivariate conjugate linear model, which can be sampled efficiently from Appendix C of [13] .To sample from the first part of the equation, we developed backsampling algorithm. The idea is motivated by the back-fitting algorithm in the Generalized Additive Model. Specifically, given the samples from F and Σ, we draw sample iteratively from p(B|F), p(f (1) |F, B), . . . , p(f (K) |F, B, f (1) , . . . , f (K-1) ). Starting with B, defineB * = F - K j=1 Θ (j) (Z (j)), then we can write:B * ∼ M N (BX, Σ, Γ * ) B ∼ N (Θ (0) , Σ, Γ (0) )where Γ * = K j=1 Γ (j) , As the above model is a matrix conjugate linear model, we can sample from its closed form:B|B * , Σ ∼ M N ((B * Γ - * X T + Θ (0) Γ -(0) )(XΓ - * X T + Γ -(0) ) -1 , Σ, (XΓ - * X T + Γ -(0) ) -1 )where Γ - * and Γ -(0) are short-hand for (Γ * ) -1 and (Γ (0) ) -1 respectively.We then use a similar process for f (k) . Definef * = F -BX - k-1 i=1 f (i) - K j=k+1 Θ (j) (Z (j) ).Then we can use a similar process to sample for f (k) :f * ∼ M N (f (k) , Σ, Γ * ) f (k) ∼ N (Θ (k) , Σ, Γ (k) )where Γ * = K j=k+1 Γ (j) and we can sample from its closed-form conditional distribution:f (k) | Σ, f * ∼ N f * Σ - * + Θ (k) Γ -(k) Γ - * + Γ -(k) -1 , Σ, Γ - * + Γ -(k) -1where Γ - * and Γ -(k) are short-hand for (Γ * ) -1 and (Γ (k) ) -1 respectively. Finally, we setf (K) = F -BX - K-1 k=1 f (k) 
empirical result: 
 We conduct a series of simulations with varying numbers of taxa D and samples N to evaluate the performance of MultiAddGPs against the transformed Gaussian model, which ignores the counting uncertainty. We implemented this model by setting H •n = ϕ(Y •n + 0.5) 
conclusion & future work: 
 We have introduced MultiAddGPs, a Bayesian Multinomial Logistic-Normal additive linear and non-linear regression model aimed at addressing statistical challenges in analyzing microbiota data.By incorporating recent advancements in Marginally Latent Matrix-t Processes (MLTPs), we've developed efficient inference methods and a framework for hyperparameter selection through Maximum Marginal Likelihood (MML) and penalized MML, implemented in the fido R package.Our approach has proven effective in both simulated and real data analyses, which offer significant flexibility in modeling. However, challenges persist in hyperparameter selection, particularly for large numbers or complex landscapes. Further research and improvement are likely needed when selecting larger numbers of hyperparameters simultaneously or when the optimization landscape is particularly challenging. Furthermore, while the Laplace Approximation in the Collapse-Uncollapse (CU) sampler performs well with high sequencing depths, it may prove inadequate for low-count scenarios. In this case, other approaches are also needed. 
 
 Appendix: Notation. In this article we denote matrix and vector dimensions with unbolded uppercase letters (e.g., N ), matricies and matrix-valued functions using bold uppercase symbols (e.g., X), vectors and vector-valued functions with bold lowercase symbols (e.g., x x x), and scalars and scalar functions as unbolted lowercase symbols (e.g., x). For matrices, we index specific rows as X d• and columns as X •n .We describe the class of Marginally Latent Matrix-t Processes (MLTPs) by sequentially generalizing from Matrix-t Processes to Latent Matrix-t Processes (LTPs) and finally MLTPs. We then describe the subset of Bayesian Multinomial Logistic-Normal MLTPs (MLN-MLTPs) before reviewing the inference of this class of models.Just as Gaussian processes can be defined based on the marginal properties of the multivariate normal, Matrix Normal and Matrix-t processes can be defined by the marginal properties of the Matrix Normal and Matrix-t distributions [13] . Matrix-t processes generalize student-T processes and Gaussian processes [13] . 2) is a matrix-t process if Y Y Y evaluated on any two finite subsets X (1) ⊂ W (1) and X (2) ⊂ W (2) is a random matrix Y of dimension |X (1) | × |X (2) | that follows a matrix-t distribution: Y ∼ T (ν, M, V, A). Here ν is scalar value strictly greater than one,Definition A.1 (Matrix-t Process). A stochastic process Y Y Y ∼ T P T P T P (ν, M M M , V V V , A A A) defined on the set W = W (1) × W (EQUATIONj ∈ X (2) ) is the matrix function representing the mean, andEQUATIONi ) andEQUATIONj ) are kernel functions. Latent Matrix-t Processes (LTP) generalize Matrix-t processes. Y Y Y is said to be an LTP if Y Y Y ∼ f (Π Π Π, λ) Π Π Π = ϕ -1 (H H H) H H H ∼ T P (ν, M M M , V V V , A A A).where f is any distribution depending on parameters Π Π Π and hyperparameters λ and ϕ is a knownA stochastic process Y Y Y is a MLTP if it has a marginal that is an LTP. More formally, a modeldescribed by a joint distribution p(Y Y Y , H H H, Φ Φ Φ) is an MLTP if it has a marginal p(Y Y Y , H H H) that is anLTP. [13] showed that a wide variety of linear, dynamic linear, and non-linear regression models are MLTP. In Appendix D, we show that our proposed class of generalized additive Gaussian process regression models are MLTP as well.Bayesian Multinomial Logistic Normal MLTPs (MLN-MLTPs) are a subtype of MLTPs that are particularly useful for the analysis of microbiome data. In MLN-MLTPs, the distributionf is a product multinomial: p(Y •1 , . . . , Y •N ) ∼ N n=1 Multinomial(Π •n ) and the transform ϕ is an invertible log-ratio transform from the D-dimensional simplex to D -1 dimentional real-space: H •n ∈ D-1 = ϕ(Π •n ∈ § D ).Canonically we used the following Additive Log-Ratio (ALR) transform which takes the D-th taxa as a reference:EQUATION)We choose this transform for computational reasons discussed in [13] and note that the choice of transformation makes little difference at the inference stages as posterior samples taken with respect to one transform can easily be transformed to samples consistent with any other choice of transform [4] . For context, note that the ALR D transform given above is the inverse of the softmax transformation.The p(H, Φ | Y ) = p(Φ | H, Y )p(H | Y )with the uncollapsed form as the first term and the posterior of the collapsed form as the second.As the collapsed form is rarely conjugate, techniques such as MCMC can be used to obtain samples from it. Then, conditioned on those samples, the uncollapsed from can be sampled to obtain samples from the posterior. Especially when Φ is high-dimensional, this Collapse-Uncollapse sampler can be much more efficient than using MCMC or Gibbs-based samplers [13] . Still, the most substantial enhancements occur when approximations to the collapsed form are considered.We have developed a Laplace approximation for the collapsed forms of MLTPs which takes the form:p(vec(H) | Y ) ≈ N (vec( Ĥ), ∇ 2 [vec( Ĥ)])where Ĥ denotes the Maximum A Posteriori (MAP) estimate of the collapsed form and We simulated a suite of longitudinal studies of microbiota with varying numbers of taxa D ∈ {3, . . . , 100} and samples N ∈ {20, . . . , 1000}. In Figure 1 we show a small simulation D = 4 and N = 600 so that results can be easily visualized. We simulated microbial composition influenced by batch effects, daily periodicity (e.g., circadian rhythm [6] ), and longer-term non-linear trends.For inference, we specify a MultiAddGP modelF n = b 0 + b 1 x (batch) n + f (periodic) (t n ) + f (trend) (t n )as follows. For covariates, we set (periodic) and f (trend) were given matrix-normal process priors with mean functionEQUATIONZ (trend) •n = t n . For priors we set B = [b 0 = 2.7; b 1 = 1]. Both fΘ (k) = 0. A periodic kernel K per (t, t ′ ) = σ per exp - 2 sin 2 π|t-t ′ | p ρ 2 perwas used in the prior for f (periodic) and a squared exponentialK se (t, t ′ ) = σ 2 se exp(-(t-t ′ ) 2 2ρ 2 se )was used for f (trend) . Hyperparameters Ψ = {σ per , ρ per , p, σ se , ρ se } were selected using MML estimation.For comparison, we created a nearly identical model that ignored uncertainty due to counting and assumed the data was transformed Gaussian. We implemented this model by setting In Figure 2 we show that our findings generalize as N and D increase. As D and N increase, it became increasingly difficult to simulate data where there was a clear enough non-linear trend to warrant additive modeling. For this reason, we replace the non-linear trend f (trend) in these simulations with a linear temporal trend and a linear kernel in the corresponding matrix-normal The first row shows the percentage of zero in each of the simulated data. In the second row, the y-axis represents the coverage ratio, and the x-axis represents the number of taxa and observations, respectively. The color represents each decomposed component, with the yellow representing the cumulative latent function F . The coverage ratio is calculated as the difference between the percentage of each method's interval that covers the true function, A coverage ratio above 1 indicates that MultiAddGP has a higher coverage percentage than NAddGP, while a ratio below 1 suggests the opposite. A higher coverage ratio signifies better performance of MultiAddGP compared to NAddGP.H •n = ϕ(Y •n + 0.process prior. As the posterior of these high-dimensional models cannot be easily visualized, we instead quantify model performance based on the coverage of posterior 95% intervals with respect to the true function decomposition. Since both MultiAddGP and NAddGP are Bayesian models we do not expect that the posterior of either method covers the truth with 95% probability: as a result we focus on the ratio between the probability of coverage from the MultiAddGP model and that of the NAddGP model. Positive values of this coverage ratio indicate that the MultiAddGP model covers the truth more often than the NAddGP model. In all simulations, we find the coverage ratio is positive.To initiate the simulation, a random seed is set to ensure reproducibility, and key parameters are defined: the number of dimensions (D) and the number of samples (N ). A sequence of time points, ranging from 0 to N -1, is generated to represent the temporal aspect of the data. These time points are standardized by subtracting the mean and dividing by the standard deviation.Samples are randomly assigned to two batches, labeled Batch1 and Batch2, in equal proportions, introducing a batch effect into the simulation.A design matrix is constructed using the batch labels and standardized time points. Hyperparameters for the covariance structures are predefined, guiding the simulation of various Gaussian processes. Specifically, an inverse-Wishart distribution is used to simulate the covariance matrix (Σ), with Σ set as a matrix with 0.9 off-diagonal elements and 1.5 diagonal elements.The following steps involve modeling different effects using Gaussian processes with distinct kernels. The daily effect, f (daily) (t n ), is captured using a periodic kernel, given by:K(t, t ′ ) = σ exp   - 2 sin 2 π|t-t ′ | p ρ 2  with σ = 4, ρ = 5, and p = 1. The hourly effect, f (hourly) (t n ), is modeled using a linear kernel:K(t, t ′ ) = σ 2 b + σ 2 v (t -c)(t ′ -c)where σ b = 20 and σ v = 1. These effects, combined with a batch effect modeled linearly with coefficients β 0 = 2.7 and β 1 = 3, form the overall effect matrix F:F .n = b 0 + b 1 x (batch) n + f (hourly) (t n ) + f (daily) (t n )The combined effects and batch influences are incorporated into the model, with Gaussian noise added to simulate realistic variability. The resultant effect matrix H is transformed using the additive log-ratio (ALR) inverse function to produce the multinomial probabilities Π:Π .n = ϕ -1 (H)Finally, the response variable Y is generated from a multinomial distribution:Y .n ∼ Multinomial(Π .n )For each evaluated combination of (D, N ), three simulated data were created and the coverage ratio in Figure 2 was calculated by the average of these three simulated data. All implements were compiled and run using gcc version 9.1.0 and R version 4.3.2. All replicates of the simulated count data were supplied to the various implementions independently and the models were fit on identical hardware, allotted 64GB RAM, 4 cores, and restricted to a 48-hour upper limit on run-time.We demonstrate MultiAddGP on a previously published longitudinal study of four artificial gut models [12] . There are two key complexities of this data. First, in order to capture both daily and sub-daily dynamics, the data consists of both daily sampling (1 month of daily samples) with a short period of hourly sample in the middle of the longer monthly period (5 days of hourly samples). As a result, the time-series is irregularly spaced and therefore ideal for Gaussian Process methods as opposed to more standard time-series methods. Second, between experiment days 11 and 13 feed lines supplying fresh media to Vessels 1 and 2 clogged leading to partial starvation in those vessels but leaving Vessels 3 and 4 unaffected. Previously we analyzed these data using Bayesian MLN dynamic linear models however our approach was limited. Our prior approach was too computationally intensive and we were limited to only modeling the dynamics of the vessels using a random walk. In particular, we were not able to isolate the effects of the starvation event from the rest of the time-series. Still, that random walk model took over 4 hours to infer in the absence of hyperparameters (which had to be hand-tuned for computational tractability).We demonstrate that our MultiAddGP framework overcomes these limitations by estimating the model, and hyperparameters, and accurately partitioning the effects of the starvation event in under 20 minutes. Our results are consistent with our prior analysis yet also reveal novel biology.Data preprocessing following the analysis of [12] resulted in 10 bacterial families and 537 samples (97.7% of sequence counts were retained). We analyze these data with the following MultiAddGPY •n ∼ Multinomial(Π •n ) Π •n = ALR -1 D (H •n ) H •n ∼ N (F •n , Σ) F = 4 v=1 f (base,v) (t n ) + I[v ∈ {1, 2}]f (disrupt,v) (t) f (base,v) ∼ GP(Θ (base,i) , Σ, Γ (v) ⊙ Γ (base) ) f (disrupt,v) ∼ GP(Θ (disrupt,i) , Σ, Γ (v) ⊙ Γ (disrupt) ) Σ ∼ InvWishart(Ξ, ζ)where Ξ was chosen to reflect a prior belief that more evolutionary similar taxa would behave more similarly [13] , ζ = D + 5 reflecting weak prior knowledge in this assumption, and ⊙ represents elementwise multiplication of kernels. The kernel Γ (v) was a block diagonal kernel with all non-zero elements equal to 1: it was used to model conditional independence between the vessels.Γ (base) was modeled as a squared exponential kernel to capture long-term non-linear trends.In contrast, Γ (distrup) was modeled using a non-stationary rational quadratic kernel. This was achieved by warping t with the sigmoid function, such that the majority of variation occurs after days 10 to 13 to capture shorter-term effects of the starvation event in Vessels 1 and 2. However, due to the use of the warp function kernel, it became difficult to locate the identified region of hyperparameters. Therefore, hyperparameters Ω = {σ base , ρ base , σ disrupt , ρ disrupt } were selected using penalized MML. Full details on penalization, kernel, and mean functions are provided in Appendix F.Consistent with prior reports, we find that the Rikenellaceae family demonstrates a clear and substantial decrease during the initial starvation period followed by a slow recovery to baseline [12, 14] . This finding is particularly interesting considering that, within living hosts, decreases in Rikenellaceae are observed in obese individuals [8, 9, 15] . Our results could either highlight the limitations of artificial gut systems for obesity research or suggest that the relationship between Rikenellaceae and obesity is modulated by host factors (e.g., the host immune system).Still, our model reveals novel biological features of these systems that were not captured in our prior models. In particular, we find that there are a number of taxa including the Synergistaceae, Ruminococcaceae, Fusobacteriaceae, and most notably the Rikenellaceae who appear to demonstrate a distinct overcoorection phenomena. That is, for each of these taxa, our model finds strong evidence that their relative abundance in the community increases (or decreases) drastically during the initial starvation event followed by a later overcorrection where their abundance ultimately decreases (or increases) from pre-starvation levels. These types of over-corrections are predicted by ecological predator-prey models which can lead to damped oscillations in response to exogenous stimuli [10] . 1) -( 4), are a type of MLTPs models with parameters Φ = {B, f (1) , . . . , f (K) , Σ} and a collapsed form p(Y , H) defined by parametersM = Θ (0) X + K k=1 Θ (k) (Z (k) ) V = Ξ A = X T Γ (0) X + K k=1 Γ (k) (Z (k) ) + I. ν = ζProof. We prove this theorem by showing that the marginal of MultiAddGP models is an LTP.By definition 2 from [13] , if p(H, Y ) is an LTP, then p(H, Y, F, f (1) , . . . , f (K) , Σ) is a MLTP model. To begin with, we note that equation ( 4) with its priors can alternatively be written as:EQUATIONEQUATIONEQUATIONEQUATIONUsing this form in combination with the affine transformation property of the matrix normal distribution, it is straightforward to marginalize over B and f (k) producing the following form:F = Θ (0) X + K k=1 Θ (k) (Z (k) ) + E B + E f (k) E B ∼ N (0, Σ, Γ (0) ) E f (k) ∼ N (0, Σ, Γ (k) ) F = Θ (0) X + K k=1 Θ (k) (Z (k) ) + E F E F ∼ N (0, Σ, X T Γ (0) X + K k=1 Γ (k) (Z (k) )Thus we may rewrite Equations ( 7) - (11) asEQUATIONEQUATIONEQUATIONBy using the definition of the matrix-t process given in section A.1, we can marginalize over F and Σ in Equations ( 12) and ( 14) to getH ∼ T P (ν, Θ (0) X + K k=1 Θ (k) (Z (k) ), Ξ, X T Γ (0) X + K k=1 Γ (k) (Z (k) ) + I)Finally, incorporating Equations ( 1) and ( 2) allows us to write the marginalized form of theY Y Y ∼ f (Π Π Π, λ) Π Π Π = ϕ -1 (H H H) H H H ∼ T P (ν, M M M , V V V , A A A).where f is multinomial distribution, ϕ is an invertible log-ratio transformation, M = Θ(0) X + K k=1 Θ (k) (Z (k) ), V = Ξ, and A = X T Γ (0) X + K k=1 Γ (k) (Z (k) ) + I, ν = ζ D.In this section, we first present the pseudo-code for the Back Sampler (BS), which efficiently samples B and f k for k ∈ 1, . . . , K. Following this, we provide the full pseudo-code for the extended Collapse-Uncollapsed (CU) sampler designed for MultiAddGPs models.1: Input: {Y, X, Z} are data observation, {F, Σ} are samples from CU sampler , Λ = {Θ (0) , . . . , Θ (k) , Γ (0) , . . . , Γ (k) } is a set of prior input 2: Output: S samples of the form (B, f (k) , k ∈ {1, . . . , K}) 3: for s = 1 to S do 4:B * = F -K j=1 Θ (j) (Z (j) ) 5:Γ * = K j=1 Γ (j) 6:Sample B|B * , Σ ∼ M N ((B * Γ - * X T + Θ (0) Γ -( 0) )(XΓ - * X T + Γ -(0) ) -1 , Σ, (XΓ - * X T + Γ -(0) ) -1 ) where Γ - * and Γ -(0) are short-hand for (Γ * ) -1 and (Γ (0) ) -1 respectively. 7:for j = 1 to K do 8:if j = 1, ..., k -1 then 9:f * = F -BX - k-1 i=1 f (i) - K j=k+1 Θ (j) (Z (j) ) 10: Γ * = K j=k+1 Γ (j) 11: Sample f (k) | Σ, f * ∼ N f * Σ - * + Θ (k) Γ -(k) Γ - * + Γ -(k) -1 , Σ, Γ - * + Γ -(k) -1where Γ - * and Γ -(k) are short-hand for (Γ * ) -1 and (Γ (k) ) -1 respectively. 12: else 13:Samplef (K) = F -BX - K-1 k=1 f (k) 14:end if 15:end for 16: end for 17: return B, f (k) Algorithm 2 The Collapse-Uncollapse (CU) Sampler for MultiAddGPs Models Before fitting the model, we standardize all continuous covariates so that their means are zero and their standard deviations are one. This standardization helps the definition of generally applicable priors for the kernel parameters. The selection of kernels and the specification of prior distributions for the kernel parameters in the analysis of artificial gut data are detailed below.For kernel selection, we employ a squared exponential kernel to model Γ (base) , which is intended to capture long-term non-linear trends. To model the disruption effects, we utilize a non-stationary rational quadratic kernel. This choice is based on the realistic assumption that the target variable might exhibit rapid changes near or after the starvation period. To do that, we incorporate an input-warping approach from [2] , defining a bijective mapping w : (-∞, +∞) → (-c, c) for the continuous time t as: w(t) = 2c × (-0.5 + 1 1 + e -a×(t-b) )where a, b and c are predefined parameters : a controls the size of the effective time window, b controls its starvation location, and c controls the maximum range effect of target variable. In our study, we set our a = 2.5, b = -1.2, and c = 2.5 to reflect the starvation occurring between days 10 -13, with the effect window extending to the end of the study. The non-stationary rational quadratic kernel is then defined as:K wrq (t, t ′ ) = σ 2 wrq (1 + (w(t) -w(t) ′ ) 2 2aρ 2 wrq ) -aIn terms of the prior setting for the hyperparameters, we use the inverse gamma prior for the length scales (ρ se , ρ wrq ) and the magnitude parameters (σ se , σ wrq ). The parameters of the inverse 



Reviews: review: title: title: The paper titled "Efficient Bayesian Additive Regression Models For Microbiome and Gene Expression Studies" presents a framework for analyzing complex sequence count data, particularly within microbiome and gene expression contexts.
main_review: main_review: The authors address key challenges inherent in complex sequence count data, particularly its compositional and count-based nature, which traditional models often overlook by treating the data as Gaussian. By introducing the MultiAddGP models within a Bayesian Multinomial Logistic-Normal (MLN) framework, the paper extends upon prior methods like generalized linear models and Gaussian process models, allowing for both linear and non-linear effects to be accounted for simultaneously.

Strengths:
-The CU sampler speeds up inference, making it much faster than traditional methods like HMC, improving scalability for large datasets.
-The MultiAddGP model can capture both linear and non-linear effects, allowing for more accurate modeling of complex microbiome and gene expression data, which often have overlapping factors.

Limitations:
-The model's reliance on the Laplace approximation may not perform well in cases with low-sequencing depth, requiring alternative methods in such situations.
-While powerful for microbiome and gene expression data, the model’s focus on compositional count data may make it less applicable for simpler datasets.

NOTE: There is a link in the paper which is breaking the anonymity rule.
title: title: This is a nice paper that improves the MLTP approach, which implies a complex covariance structure, by adding the deconvoluted overlapping linear and non-linear effects. The earlier approach of MLTP either can have a non-linear effect or multiple disjoint linear effects. The motivation of the paper is to improve the estimation of the linear and non-linear effects of covariates on microbial or gene compositions.
main_review: main_review: 1. My biggest problem with the paper is how the authors clearly reveal their identity by mentioning “we” and referring to their previous work in the reference.
2. The paper could be of more value if it were compared to mixed models since those also include the additive concept. And it could be that Mixed Models are simpler and better suited.
3. I am a bit unsure if the simulated data is generated in a fair manner, but at least it shows that some logic is working in their method.
4. The paper sounds really powerful until I read the algorithm, seeing that we actually have observable variables signaling clearly what is linear and non-linear. This makes the whole paper less useful as, in many biological applications, we do not have access to this level of information, and we merely have one observation, which is the result of collided linear and nonlinear covariates. I may be misunderstanding this last point here but if it is sensible I am hoping to be useful for the reader of this comment.
5. Finally, I think this work is absolutely valuable and sound and it should be accepted.