Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes
Justin D. Silverman Kimberly Roche Zachary C. Holmes Lawrence A. David Sayan Mukherjee


1,2]Justin D. Silverman
1]Kimberly Roche
5]Zachary C. Holmes
1,4,5,6]Lawrence A. David
1,3,4,6]Sayan Mukherjee



[1]Program in Computational Biology and Bioinformatics, Duke University, Durham, NC 27708
[2]Medical Scientist Training Program, Duke University, Durham, NC 27708
[3]Departments of Statistical Science, Mathematics, Computer Science, Biostatistics & Bioinformatics, Duke University, Durham, NC 27708
[4]Center for Genomic and Computational Biology, Duke University, Durham, NC 27708
[5]Department of Molecular Genetics and Microbiology, Duke University, Durham, NC 27708
[6]Denotes Co-corresponding Authors


[
    [
    December 30, 2023
=====================








   Bayesian multinomial logistic-normal (MLN) models are popular for the analysis of sequence count data (e.g., microbiome or gene expression data) due to their ability to model multivariate count data with complex covariance structure. However, existing implementations of MLN models are limited to handling small data sets due to the non-conjugacy of the multinomial and logistic-normal distributions. We introduce MLN models which can be written as marginally latent matrix-t process (LTP) models. Marginally LTP models describe a flexible class of generalized linear regression, non-linear regression, and time series models. We develop inference schemes for Marginally LTP models and, through application to MLN models, demonstrate that our inference schemes are both highly accurate and often 4-5 orders of magnitude faster than MCMC.




  
Keywords: Bayesian Statistics; Multivariate Analysis; Count Data; Microbiome; Gene Expression


§ INTRODUCTION


In this article we develop methods for efficient inference of Bayesian multinomial-logistic normal (MLN) models. MLN models are used for the analysis of compositions measured through multivariate counting. In contrast to multinomial Dirichlet models, multinomial logistic-normal models permit both positive and negative covariation between multinomial categories <cit.>. While multinomial logistic-normal topic models have been used in natural language processing for some time <cit.>, more recently these models have been adopted for regression and time-series modeling of microbiome data <cit.>. 

Yet, inference in MLN models is difficult due to lack of conjugacy between the multinomial and the logistic normal. Early work with MLN models used Metropolis within Gibbs samplers <cit.> and could scale to just a small number multinomial categories (i.e., less than 5). Recently, Pólya–Gamma data augmentation was proposed as a means of inference in MLN regression by augmenting Pólya–Gamma random variables between the multinomial and logistic normal components of a model and sampling each variable as a separate Gibbs sampling step <cit.>. Numerous authors have found this approach too computationally intensive to scale to large multinomial models and have instead developed augmentation methods based on a stick-breaking representation of the multinomial  <cit.>. However, this stick breaking representation does not maintain the logistic-normal form of the model and is sensitive to the labeling of multinomial categories <cit.>. Most recently, 
several authors <cit.> have shown that Hamiltonian Monte Carlo (HMC) provides for a more efficient and scalable approach to inference in MLN models. In particular, <cit.> used a HMC within a Gibbs sampler whereas both <cit.>, and <cit.> found that the No-U-Turn-Sampling algorithm provided by the Stan Modeling language <cit.>, provided more scalable inference. However, both these approaches are still limited in the number of categories or samples that they can handle. <cit.> analyzed approximately 800 samples each with only 10 multinomial categories; <cit.> analyzed 36 multinomial categories but had to run their model over the dataset using a sliding window of 60 samples at a time; and <cit.> analyzed 166 samples and 2662 categories but had to impose low rank structure on the logistic normal model for computational tractability. In this work we show that our inference methods scale to hundreds to thousands of categories and samples and permit inference for a wide variety of models including non-linear regression models (as in <cit.>), dynamic linear models (as in <cit.>), and linear regression models (as in <cit.>).


The key ideas developed in this article are three-fold. First, we demonstrate that a number of popular multivariate Bayesian models are special cases of what we term marginally latent matrix-t process (Marginally LTP) models. As the name suggests, Marginally LTP models share a canonical marginal form which is a generalization of Student-t processes <cit.> to multivariate problems and alternative observation distributions. Second, through applications to MLN linear models, we demonstrate that the LTP representation can provide efficient inference schemes for Marginally LTP models. In particular our inference is based on collapsing a target model to its LTP form, obtaining samples from that LTP, and then uncollapsing the samples (i.e. sampling form a particular posterior conditional) so as to produce samples from the original target model. Finally, we show that for multinomial-logistic normal models, a Laplace approximation to the LTP form can provide efficient and accurate inference for the original target model.  


The layout of this article is as follows. In Section 2 we introduce the MLN linear models which motivate subsequent sections. In Section 3 we introduce the class of Marginally LTP models which encompasses MLN linear models. In Section 4 we develop inference methods for Marginally LTP Models. In Section 5 we demonstrate our approaches through an extensive simulation studies of MLN linear models. In Section 6 we demonstrate the utility of MLN linear Models through application to a microbiome sequence count study. Finally, we close with a discussion in Section 7. 


§ MULTINOMIAL LOGISTIC-NORMAL (MLN) LINEAR MODELS
 

Our primary motivation in this work was to develop efficient inference for a class of models we refer to as multinomial logistic-normal (MLN) linear models. Consider a dataset consisting of paired observations Y and X where Y_· j denotes the j-th D-variate count vector with a paired Q dimensional covariate vector X_· j. For example, in the analysis of microbiome data we may consider Y_· j to be the number of D different bacterial taxa sequenced in the j-th sample and X_· j to be Q different covariates describing the j-th sample. We denote the number of samples as N. In this setting we define the class of multinomial logistic normal linear models as 

    Y_· j   ∼Multinomial(π_· j)  
    π_· j   = ALR_D^-1(η_· j) 
    η_· j   ∼ N(Λ X_· j, Σ) 
    Λ   ∼ N(Θ, Σ, Γ)
    Σ   ∼ IW(Ξ, υ).

where ALR^-1_D is the inverse additive log-ratio transform  given by <cit.>

    ALR_D^-1(η_· j)=(e^η_1j/1+∑^D-1_i=1e^η_ij,⋯,e^η_(d-1)j/1+∑^D-1_i=1e^η_ij, 1/1+∑^D-1_i=1e^η_ij).

ALR^-1_D is also known as the multivariate-logit and is identical to the softmax transform where one component is set to zero for identifiability. In this way Equations (<ref>) and (<ref>) is equivalent to the multivariate-logit parameterized multinomial commonly used in multinomial regression models. Importantly, the pull-back of the multivariate normal measure through the inverse ALR transform is the logistic-normal measure (<cit.>) and we therefore refer to this class of models as multinomial logistic-normal linear models. The definition of ALR^-1 implies that η is a real valued matrix of dimension (D-1) × N, Σ is a covariance matrix of dimension (D-1) × (D-1) and Γ is a covariance matrix of dimension Q × Q. Additionally we note that the term Λ∼ N(Θ, Σ, Γ) denotes a matrix-normal distribution defined by (⃗Λ) ∼ N((⃗Θ), Γ⊗Σ) where ⊗ denotes the Kronecker product and $⃗ denotes the vectorization operation that converts a matrix into a column vector.

MLN linear models provide a means of inferring log-linear interactions between compositions measured through multivariate counting and covariates. Each elementΛ_ikdescribes the linear effect of thek-th covariate on thei-th log-ratio coordinate. In contrast to standard multinomial logistic regression, the above model assumes that each observation is subject to extra-multinomial variation[This is reflected in the fact that Equation (<ref>) reads as η_· j∼ N(Λ X_· j, Σ) rather than η_· j = Λ X_· j. ] and uncertainty regarding the covariation between the rows ofηandΛ(Equation (<ref>)).  Additinoally, the prior specified in Equations (<ref>) and (<ref>) produce shrinkage of estimates ofΛin settings whereD>N.  

We will show in the next section that MLN linear models are one example of a larger class of models which we term marginally latent matrix-t process (LTP) models. We will show that for MLN linear models this canonical form is characterized by marginalizing overΛandΣso as to replace Equations (<ref>)-(<ref>) with a stochastic process we term a matrix-t process. In Section <ref> we build efficient inference methods for MLN linear models by developing tools for inference of Marginally LTP models. 


§ MODELING OVERVIEW

In this section we will introduce marginally latent matrix-t process (Marginally LTP) models as a flexible class of models capable of describing a wide variety of linear regression, non-linear regression, and time-series models. To build Marginally LTP models we first describe matrix-t processes, then build latent matrix-t processes (LTP), and then generalize this class to Marginally LTP models. 



 §.§ Matrix-Normal and Matrix-T, Distributions and Processes

To build the class of Marginally LTP models we first review matrix-normal distributions and processes as well as matrix-t distributions and processes highlighting properties we will make use of in this article.



  
Matrix-Normal Distribution The matrix-normal distribution is a generalization of the multivariate normal distribution to random matrices. We describe a randomm×nmatrixXas being distributed matrix-normalY∼N(M, U, V)if(⃗Y) ∼N((⃗M), V ⊗U)whereUis am×mcovariance matrix andVis an ×ncovariance matrix. 



  
Matrix-Normal Process We define a stochastic process𝖸as a matrix-normal process on the set𝒳 = 𝒳^(1) ×𝒳^(2)and denoted𝖸 ∼𝖦𝖯(𝖬, 𝖪, 𝖠)if𝖸evaluated on any two finite subsets𝐱^(1)=(x_1^(1), …, x_P^(1)) ∈𝒳^(1)and𝐱^(2)=(x_1^(2), …, x_N^(2))∈𝒳^(2)is distributed asY∼N(M, K, A)whereM_ij = 𝖬(x^(1)_i, x^(2)_j),K_ij = 𝖪(x^(1)_i, x^(1)_j),A_ij = 𝖠(x^(2)_i, x^(2)_j)for matrix function𝖬and kernel functions𝖪and𝖠. The requirement that𝖪and𝖠be kernel functions implies that the matricesKandAare covariance matricies (i.e., they are symmetric positive definite). 




  
Matrix-t Distribution  The matrix-t distribution is a generalization of the multivariate-t distribution to random matrices. Like the multivariate-t, the matrix-t can be defined constructively through its relationship to the matrix-normal and inverse Wishart distributions. LetΣdenote a random covariance matrix such thatΣ∼IW(Ξ, υ)whereΞrepresents a positive semi-definite scale matrix andυ> 0. Also suppose thatX∼N(0, I, V). IfCC^T=Σthen the distribution ofY=CXis denoted as matrix-t such thatY∼T(υ, 0, Ξ, V). For a random matrixη∼T(υ, B, K, A)the log density ofηmay be written

    log T_P× N(η|υ, B, K, A)  = logΓ_P(υ+N+P-1/2) -  logΓ_P(υ+P-1/2) -NP/2logπ
    
        -N/2log|K| - p/2log|A| 
        -υ + N + P - 1/2log| I_p + K^-1[η-B]A^-1[η - B]^T |

whereΓ_a(b)refers to the multivariate gamma function. These results follows directly from <cit.>. 



  
Matrix-t Process Through analogy to our definition of matrix normal processes, we define a matrix-t process through its relationship to the matrix-t distribution. We define a stochastic process𝖸 ∼𝖳𝖯(υ, 𝖡, 𝖪, 𝖠)defined on the set𝒳 = 𝒳^(1) ×𝒳^(2)as a matrix-t process if𝖸evaluated on any two finite subsets𝐱^(1)=(x_1^(1), …, x_P(1)) ∈𝒳^(1)and𝐱^(2)=(x_1^(2), …, x_N^(2))∈𝒳^(2)is distributed asY∼T(υ, B, K, A)whereυis a scalar strictly greater than zero,B_ij=𝖡(x^(1)_i, x^(2)_j),K_ij = 𝖪(x^(1)_i, x^(1)_j), andA_ij=𝖠(x^(2)_i, x^(2)_j)for matrix function𝖡, and kernel functions𝖪and𝖠. Matrix-t processes can be alternatively seen as a multivariate generalization of Student-t processes which have found widespread use in statistical analysis <cit.>.



 §.§ Latent Matrix-t Processes (LTPs)

To allow us to generalize matrix-t processes to a more flexible set of data types, for example count data, we now define LTPs as a generalization of a matrix-t processes. In essence, the following definition generalizes matrix-t processes by defining a stochastic process𝖸as a hierarchical process formed by a process𝖥having parameters that, with appropriate transformationϕ, follow a matrix-t process. LTPs allows us to flexibly model complex patterns of covariation in a wide variety of data types. For example, with LTPs we can develop multinomial models with a more flexible covariance structure than the multinomial-Dirichlet, which has strong independence assumptions between multinomial categories <cit.>. The final change compared to our definition of matrix-t processes, is that we now explicitly denote dependence on model hyperparameters which we collectively refer to asδ.  

Latent Matrix-t Process We define a stochastic process 𝖸 as a latent matrix-t process 𝖸∼𝖫𝖳𝖯(𝖥, ϕ, υ,𝖡, 𝖪, 𝖠, δ) on the set 𝒳 = 𝒳^(1)×𝒳^(2) if 𝖸 evaluated on any P dimensional finite subset x^(1)∈𝒳^(1) and any N dimensional finite subset x^(2)∈𝒳^(2) is distributed 

    Y    ∼ f(π, δ) 
    π   = ϕ^-1(η) 
    η   ∼ T(υ, B(δ), K(δ), A(δ)).

where η denotes a P × N real valued matrix, B(δ) a P× N dimensional real valued matrix function of parameters δ defined by [B(δ)]_ij = 𝖡(x^(1)_i, x^(2)_j, δ), K(δ) is a P× P covariance matrix defined as [K(δ)]_ij = 𝖪(x^(1)_i, x^(1)_j, δ),  A(δ) is an N× N covariance matrix defined as [A(δ)]_ij = 𝖠(x^(2)_i, x^(2)_j, δ), υ is a scalar subject to υ>0, π is an element of a space Π defined via the one-to-one mapping ϕ^-1:ℛ^P× N→Π, f denotes a likelihood function with parameters π and δ which is itself an evaluation of the process 𝖥 evaluated on a finite subset of the set Π. 















Before proceeding further we motivate the generality of the above definition with two examples.



  §.§.§ LTPs for analysis of Generalized Linear Models
 
Our first example demonstrates how LTPs can be used to infer Bayesian generalized linear models. 
Consider the following univariate exponential generalized linear model (GLM) for observed dataYand covariatesXY_j    ∼Exponential(λ_j) 
    λ_j    = -1/η X_· j
    η   ∼ N(0, Σ) 
    Σ   ∼ IW(Ξ, u).

Noting thatη= N(0, Σ)can alternatively be written asη∼N(0, Σ, 1)and noting the definition of the matrix-t, we can stateηmarginalizing overΣasη∼T(u, 0, Ξ, 1). Thus, marginalizing overΣin the above model results in an LTP with parametersf = ∏_i=1^N Exponential(Y_i |λ_i), υ=u, B = 0, K = Ξ, A = 1, δ= ∅where∅denotes the empty set andϕ^-1is defined by the relationλ_j = [ϕ^-1(η)]_j = -1/ηX_·j. This result demonstrates that inference ofηin the above univariate exponential GLM can instead be replaced by inference of an LTP with the above parameters.   
 


  §.§.§ LTPs for Non-Linear Modeling of Microbiome Time-Series

Motivated by the non-linear time-series model in <cit.>, our second example illustrates how LTPs can be used to construct a non-linear multinomial time-series model for the analysis of microbiome data.
LetYdenote aD×Tmatrix of counts such thatY_·jrepresents the sequence count vector ofDtaxa measured in the sample from thet-th time-point. This data could be modeled as 

    Y_· t   ∼Multinomial(π_· t) 
    π_· t   =ALR^-1(η_· t) 
    η   ∼ T(υ, B, K, A) 
    
        B_· t   = 0_d-1
    
        K_i,j   = κ^2exp(-γ^-2[d_phylo(i,j)]^2) 
    
        A_t,s   = α^2 exp ( -ρ ^-2 (t-s)^2)

whereA_t,sdenotes a radial basis function kernel describing the similarity between time-pointstandsandd_phylo(i,j)represents the phylogenetic (evolutionary) distance between microbial log-ratiosiandjas introduced in <cit.>. Intuitively, this specification models a microbial time-series as a composition evolving through time subject to non-linear dynamics and with the prior assumption that evolutionary more similar taxa are likely to behave more similarly. Additionally, the termB_·t = 0_D-1can be interpreted as saying we have little prior knowledge regarding which taxa are the most abundant at any time-point[When all ALR coordinates are equal to zero this corresponds to π_i=1/D for all taxa i.]. Additionally the relationship between the matrix-t, the inverse Wishart and the matrix-normal implies that this model can also be seen as a multinomial, log-ratio transformed matrix-normal (i.e., matrix-variate version of a logistic-normal) with uncertainty in the covariance matricies of the matrix-normal described by an inverse Wishart distribution. In this way the specification ofKandAis likely less sensitive to miss-specification than if we hadη∼N(B, K, A). In particular, this implies that the scalarυrepresents our uncertainty in our specification ofKandA. 

Finally, this non-linear time-series model can be seen to represent an LTP with the following specifications:f = ∏_t=1^T Multinomial(Y_·t |π_·t), 
    ϕ^-1 = ALR^-1,
    υ= υ,
    B_·t = 0_D-1,
    K_i,j = κ^2exp(-γ^-2[d_phylo(i,j)]^2),A_t,s = α^2 exp( -ρ^-2 (t-s)^2),
    δ= { κ, γ, α, ρ} which together also imply that𝐱^(1) = {1, …, D-1} ∈𝒳^(1),𝐱^(2) = {1, …, T} ∈𝒳^(2), and correspondingly thatP = D-1andN=T. 
















 §.§ Marginally LTP Models

We next generalize our definition of LTPs to a larger class which we term Marginally LTP models. This definition is straight forward, we define Marginally LTP models as those models which have a marginal that is an LTP. 

Marginally LTP models If a model described by the joint distribution p(η, Ψ, Y) may be written as p(Ψ|η, Y)   p(η,Y)  where p(η,Y) is a latent matrix-t process, we refer to p(η, Ψ, Y) as a Marginally LTP model and p(η, Y) as the model's collapsed representation. 


In the next two subsections we demonstrate that Marginally LTP models provide a rich class of models. We give two examples of Marginally LTP models, the first represents a generalization of the multinomial logistic normal linear models introduced in Section <ref>, the second represents a flexible class of models for inference in multivariate non-Gaussian time-series. 



  §.§.§ Generalized Multivariate Conjugate Linear (GMCL) Models

As in Section <ref>, let us considerYto representNindependentD-variate measurements and considerXto representNsets ofQ-dimensional covariates. We define generalized multivariate conjugate linear (GMCL) models as 

    Y_· j   ∼ f(π_· j)  
    π_· j   = ϕ^-1(η_· j)
    η_· j   ∼ N(Λ X_· j, Σ) 
    Λ   ∼ N(Θ, Σ, Γ) 
    Σ   ∼ IW(Ξ, υ) .

We may describe the joint density of this model asp(Λ, Σ, η, Y)which can be factored asp(Λ, Σ|η, Y)   p(η, Y). Therefore, in parallel to the definition of Marginally LTP models we may equateΨ= {Λ, Σ}.  In Appendix <ref> we show thatp(η, Y)is LTP matrix-t with parameters

    B    = Θ X
    
        K    = Ξ
    
        A   =I_N+X^TΓ X

and where parametersf,ϕ, andυare defined flexibly as in Equations (<ref>) and (<ref>). Finally this model specification implies that{Θ, Γ, Ξ} ∈δ. This result demonstrates that all GMCL models are Marginally LTP models. Further, by lettingfdenote the multinomial distribution andϕ^-1denote the inverse ALR transform, we see that the multinomial logistic-normal linear models introduced in Section <ref> are a special case of GMCL models. 




  §.§.§ Generalized Multivariate Dynamic Linear Models (GMDLMs)

 
We develop a flexible class of multivariate time-series models for non-Gaussian observations. We term this class of models generalized multivariate dynamic linear models (GMDLMs). GMDLMs represent an extension of the multivariate dynamic linear models introduced in <cit.> and developed further in <cit.> to non-Gaussian observations. Using the notation from <cit.> we letη_t^Tdenote a row-vector (i.e., the transpose of thet-th column ofη). We define the GMDLM as

    Y_· j   ∼ f(π_· j)  
    π_· j   = ϕ^-1(η_· j) 
    η_t^T    = F_t^TΘ_t + ν_t^T,   ν_t ∼ N(0, γ_t Σ) 
    Θ_t    = G_t Θ_t-1 + Ω_t,   Ω_t ∼ N(0, W_t, Σ) 
    Θ_0    ∼ N(M_0, C_0, Σ)  
    Σ   ∼ IW(Ξ, υ)

whereΘ_trepresents aQ ×Pmatrix describing the state of the time-series at timet,G_tdenotes theQ×Qstate transition matrix at timet,F_tdenotes a1×Qvector describing a linear model relating the latent space to the parametersη_t,Σis aP ×Pcovariance matrix specifying the covariation between thePdimensions of the time-series,W_tis aQ×Qcovariance matrix describing the covariation of the perturbations affecting latent states, andγ_tis a  scalar allowing an analyst to weight the importance of select observations (γ_tis typically equal to1).  

The joint model for the GMDLM can be writtenp(Θ, Σ, η, Y)which can be factored asp(Θ, Σ|η, Y)   p(η, Y). To parallel the definition of Marginally LTP models, here we haveΨ= {Θ, Σ}. In Appendix <ref> we demonstrate thatp(η, Y)is an LTP with parameters 

    B    = [                ; α_1   ⋯ α_t   ⋯ α_T;                     ]
    α_t    = (F_t^T𝒢_t:1M_0)^T 
    
        K    = Ξ
    
        A_t,t-k   = 
          γ_t + F_t^T[
            W_t + ∑_ℓ=t^2𝒢_t:ℓW_ℓ-1𝒢^T_ℓ:t + 𝒢_t:1C_0𝒢_1:t^T
            ]F_t  if  k=0
    
            F_t^T[𝒢_t:t-k+1  W_t-k +
        ∑_ℓ=t-k^2 𝒢_t:ℓ W_ℓ-1G^T_ℓ:t-k
        + 𝒢_t:1 C_0G^T_1:t-k]F_t-k if  k>0

where we have introduced𝒢_t:ℓas a short hand notation for the productG_t⋯G_ℓ. Additionally here LTP parametersυ,f, andϕare defined flexibly as in Equations (<ref>) and (<ref>) and we have hyperparameters{ Ξ, M_0, C_0, W_1, …, W_T, γ_1, …, γ_T, G_1, …, G_T, F_1,…, F_T }∈δ. This result demonstrates that GMDLMs are a special case of Marginally LTP models. 



§ INFERENCE IN MARGINALLY LTP MODELS

Our overarching goal was to develop efficient and accurate posterior inference for multinomial logistic-normal linear models which are a special case of Marginally LTP models. In this section we demonstrate how the canonical LTP form of Marginally LTP Models can be exploited for efficient inference of this larger model class. We consider two types of parameters,ηwhich are distributed matrix-t andΨwhich are marginalized out of a model to produce a LTP form. The inference scheme we introduce will involve two steps: 1) samplingηconditioned on observed data (p(η|Y)), 2) samplingΨconditioned onηand observed data (p(Ψ|η,  Y)). In Section <ref> we introduce this sampling scheme which we refer to as the collapse-uncollapse (CU) sampler. In Section <ref> we further build on the CU sampler by introducing a Laplace approximation as a means of samplingp(η|Y). In Sections <ref> we discuss efficient means of samplingp(Ψ|η, Y)for the GMCL model and GMDLM models introduced in the previous section.  



 §.§ The Collapse-Uncollapse (CU) Sampler

Consider the task of sampling from the posterior distribution of a Marginally LTP model with joint densityp(Ψ, η, Y). The corresponding posterior density can be decomposed as

    p(η, Ψ| Y) = p(Ψ|η, Y)p(η, Y)/p(Y).

This decomposition implies the following: given a model described by the joint distributionp(η, Ψ, Y), we may sample from the posterior by first sampling from the posterior of the collapsed (LTP) modelp(η,Y)and then given that sample ofηand the observedYwe may then sampleΨfrom the conditionalp(Ψ|η, Y). Together the sample ofηandΨthen represents a single sample from the posterior of the Marginally LTP model,p(Ψ, η|Y)(Algorithm <ref>). 



The CU Sampler for Marginally LTP Models therefore requires two features for effcient inference. First, we require an efficient means of producing samples from the collapsed (LTP) formp(η|Y). As we will show in Section <ref>, samplingp(η|Y)can be more efficient than samplingp(Ψ, η|Y)just by virtue of the fact that the former has fewer dimensions. Therefore the CU sampler alone can be more efficient than sampling the full (uncollapsed) model. Still, in Section <ref> we develop a Laplace approximation forp(η|Y)which can further improve efficiency. Second, we require an efficient means of sampling from the posterior conditionalp(Ψ|η, Y). In Section <ref> we discuss efficient means of samplingp(Ψ|η, Y)for the GMCL and GMDLM models.  

Our rationale for focusing on the CU sampler for inference in Marginally LTP models is as follows. We expect that many Marginally LTP models (such as those introduced in Section <ref>) have partial conjugacy. Methods such as Metropolis-within-Gibbs have been popular methods of exploiting this conjugacy <cit.>. Yet, by embedding MCMC steps within a Gibbs sampler techniques such as adaptation <cit.> or approximate methods such as Laplace approximations may not make sense. In contrast, the CU sampler allows the non-conjugate sampling to occur up front (in the sampling ofp(η|Y)) so that such techniques can be used. Moreover, after multiple samples ofηhave been produced, uncollapsing the model can be done in parallel for each sample ofη. Therefore, the CU sampler may be advantageous as it permits the use of adaptive or approximate methods for sampling the non-conjugate model components and permits a degree of parallelism not allowed by Metropolis-within-Gibbs.



 §.§ Laplace Approximation for the Collapsed Form

Samplingp(η|Y)is often the major computational bottleneck when inferring Marginally LTP models via the CU sampler. To accelerate inference we developed a Laplace approximation for the densityp(η|Y). This approximation is defined asq(η|Y) = N( η̂, H^-1(η̂))whereη̂denotes the MAP estimate ofp(η|Y)andH^-1(η̂)denotes the inverse Hessian matrix oflogp(η|Y)evaluated at the pointη̂whereη̂represents the maximum a posteriori (MAP) estimate. That isη̂is defined as the solution to the following optimization problem

    η̂ = _η∈ℛ^P× N[-log p(η| Y)].


We hypothesized that such a Laplace approximation would provide an accurate approximation to an LTP posterior based on the following observations. First, the matrix-normal can provide a good approximation for the matrix-t for suitably largeυ<cit.> as it is both globally symmetric and log-convex about the MAP estimate. Second, all exponential family likelihoods are log-convex with respect to their natural parameters <cit.>. Notably, the ALR parameters represent the natural parameters of the multinomial distribution. Together these features suggested that the posterior of a LTP could be approximated as globally symmetric and log-convex like the multivariate normal. 

Developing an efficient Laplace approximation for LTP models required closed form solutions for the gradient and Hessian of LTPs. To develop these tools note that by Bayes rule we may write 

    -log p(η| Y)    ∝ - log f(Y |ϕ^-1(η))  + p(η) .

By linearity of the derivative operator we may write the gradient and Hessian of-logp(η|Y)as

    -dlog p(η| Y)/d(⃗η)   = -d log f(Y |ϕ^-1(η))/d(⃗η)- dlog p(η)/d(⃗η)
    
        -d^2log p(η|Y)/d(⃗η)d(⃗η)^T   = -d^2 log f(Y |ϕ^-1(η))/d(⃗η)d(⃗η)^T- d^2log p(η)/d(⃗η)d(⃗η)^T.

Thus we find that calculating the gradient and Hessian of LTPs reduces to calculating the gradient and hessian oflogf(Y |ϕ^-1(η))and the matrix-t densitylogp(η|X)separately. The additive structure of the gradient and Hessian are central to generalizing this approach to a variety of different observation distributionsfand transformationsϕ^-1. In Appendix <ref> we provide the gradient and Hessian for the matrix-t density. With these results, to derive a flexible class of multinomial logistic-normal models we only need to provide the gradient and Hessian for the logit-parameterized multinomial which we state in Appendix <ref>. 



 §.§ Efficient Sampling from Posterior Conditionals

The second step of the CU sampler involves sampling from the densityp(Ψ|η, Y). While the density ofp(Ψ|η, Y)is specific to the particular Marginally LTP model, we develop efficient means of sampling from this density for the GMCL and GMDLM models in Appendices <ref> and <ref> respectively. In particular, for both the GMCL and GMDLM models we make use of the fact thatΨis conditionally independent ofYgivenη, that isp(Ψ|η, Y) = p(Ψ|η). This conditional independence also reduces sampling from the conditionals to computing the posterior distribution of standard Bayesian multivariate linear regression for the GMCL model and conjugate multivariate dynamic linear models for the GMDLM model. That is, for both the GMCL and GMDLM models sampling the conditionals reduces to posterior inference for equivalent Bayesian Gaussian models that have been well described previously. 



 §.§ Software for Multinomial Logistic-Normal Linear Models

For inference of multinomial logistic normal linear models we developed the R package stray<cit.>. Stray implements the CU sampler with Laplace approximation described above using optimized C++ code. Estimation ofη̂is performed using the L-BFGS optimizer which we have found provides efficient and stable numerical results. 

Additionally all code required to reproduce the results of the next two sections, including the alternative implementations of multinomial logistic-normal linear models discussed <ref> is available as a GitHub repository at https://github.com/jsilve24/stray_paper_code. 


§ SIMULATIONS

We performed a series of simulation studies to evaluate both the CU sampler and the Laplace Approximation in terms of accuracy and efficiency of posterior inference in multinomial logistic-normal linear models. To evaluate the utility of the collapse-uncollapse sampler we compared Hamiltonian Monte Carlo (HMC) of the full model (HMC Uncollapsed) to the CU sampler where sampling of the collapsed (LTP) form was performed using HMC (HMC Collapsed). Both HMC implementations were inferred using the highly optimized  No-U-Turn-Sampler provided in the Stan modeling language <cit.> which has been frequently for the analysis of MLN models <cit.>. To further evaluate the utility of the Laplace approximation to the collapsed form in the CU sampler (LA Collapsed), we used the stray software package described in Section <ref>. Finally, to compare LA Collapsed to an alternative scheme for approximate inference, we included two mean-field automatic-differentiation Variational Bayes (VB) implementations <cit.>. The first was a VB approximation to the full form (VB Uncollapsed), the second was a VB approximation to just the collapsed form of the CU sampler (VB Collapsed). As VB Uncollapsed was unstable in practice, often resulting in error during optimization, only the results form VB Collapsed could be shown below. 

In order to compare these implementations, we created a series of simulations based on the likelihood model described in Equations (<ref>)-(<ref>). We identified three key parameters,  the sample size (N),  the observation dimension (D), and the number of model covariates (Q) which we varied in order to test each implementation over a wide variety of conditions. We choose the tuple{N=100, D=30, Q=5}as our base condition and independently varied each simulation parameter from that base condition (Nfrom 10 to 1000,Dfrom 3 to 500, andQfrom 2 to 500). Importantly, since we hypothesized that the Laplace and variational approximations would be challenged at high levels of zero values, we developed these simulations such that the proportion of zero counts also increased with increasingDandQ(Figure <ref>). Additionally, to account for the stochastic nature of the simulations, three simulations were performed for each tuple{N, D, Q}. For each simulation, each of the five implementations were fit and allowed a maximum of 48 hours to produce 2000 samples. Further details of the simulation and model fitting procedure can be found in Appendix <ref>. 

To quantify the accuracy and efficiency of each implementation we defined the following performance metrics. As a measure of efficiency, we calculated the average number of seconds needed for the implementation to produce one independent sample from the target posterior (i.e., Seconds per Effective Sample - SpES). Specifying independent samples is important as HMC samplers produce autocorrelated samples. To quantify the accuracy of point estimates from each implementation (i.e., either the posterior mean or MAP) we defined the root mean squared error of the point estimate forΛfrom its true simulated value.  Notably, given finiteNwe do not expect that any implementation will be able to perfectly reconstruct the true simulated value forΛ; rather, this metric provides a means of comparing the relative performance of each implementation. Finally, to quantify the accuracy of uncertainty quantification from each implementation we compared posterior intervals against those of the HMC Collapsed model which was taken as a gold standard. In particular, we define the root mean squared error of standard deviations as the average difference between the estimated posterior standard deviations,sd(Λ_ij), compared to the estimates produced by HMC Collapsed.

Overall, we hypothesized that LA Collapsed would be more efficient than the other implementations but would have low accuracy of uncertainty quantification when the percent of zeros in the observed data was higher than 30%. The later portion of this hypothesis is based on the observation that all exponential family distributions parameterized by their natural parameters (e.g., a multinomial parameterized by log-ratio coordinates) is globally log-convex <cit.> whereas the matrix-t distribution is only log-convex near the mean. We therefore hypothesized that a weaker likelihood would lead the Laplace approximation to have thinner tails than the true posterior. In practice LA Collapses provided nearly identical estimates of posterior uncertainty to both HMC implementations up to over 90% data sparsity. Additionally, LA Collapsed provided nearly identical point estimates to both HMC implementation over the full spectrum of simulations. Finally, LA Collapsed was often up to 5 orders of magnitude faster than HMC and often 1-2 orders of magnitude faster than VB. 



 §.§ Computational Efficiency


Overall, the CU sampler with a Laplace approximation (LA Collapsed) provides the most efficient inference across all tested conditions. More specifically LA Collapsed displays speed ups of between 1 to 5 orders of magnitude in comparison to HMC Collapsed and Uncollapsed and often between 1-2 order of magnitude compared to VB Collapsed. Notably, HMC Uncollapsed fails to complete sampling within 48 hours forD>100. 

Beyond the high efficiency of LA Collapsed, our results also demonstrate that the CU sampler can improve inference in HMC without the use of approximate inference methods. These results likely stem from the smaller number of dimensions in HMC for the collapsed versus uncollapsed implementations. Most noticeably, the collapsed representation completely removes dependency onQfrom HMC run-times asΛis marginalized out of the collapsed representation. However, for largeNthe HMC Uncollapsed is more efficient than HMC Collapsed. This result may reflect that the heavy tails of the matrix-t distribution produce a more challenging geometry for HMC than the expanded  matrix normal and inverse Wishart forms. Such a result has been well described previously for both univariate and multivariate-t distributions <cit.>. 




 §.§ Point Estimation


Overall point estimation using LA Collapsed (i.e., MAP estimates) is nearly identical to point estimation using either HMC Collapsed or HMC Uncollapsed (i.e., mean estimates). In contrast point estimation using VB Collapsed can produce substantially larger errors, especially for large values ofD. Overall these results demonstrate that the CU sampler maintains accuracy in point estimation and that MAP estimation provides an excellent approximation to the mean in multinomial logistic normal linear models. 




 §.§ Uncertainty Quantification

Beyond accuracy of point estimates, we also wanted to study the accuracy of estimates of uncertainty from each implementation. We consider the HMC Collapsed implementation to be the gold standard which we base our metric RMSE of standard deviations on. Except for values ofQgreater or equal to 250 (where the proportion of zero values is>90%), the uncertainty estimates of LA Collapsed are identical to those of both HMC implementations. Yet, at larger values ofQ, when sparsity is>90%, we see differences not only between LA Collapsed and HMC but between the two HMC implementations themselves. There are two possible explanations for this. First, that stray has a slightly better point estimation accuracy in these same largeQsimulations could point to the fact that stray is correct and instead HMC estimates of uncertainty are incorrect due to the often small effective sample size for largeQ. Alternatively, this could support our previous hypothesis that the Laplace approximation has higher error in uncertainty quantification with higher data sparsity. Given the ergodicity of HMC it seems more likely that the Laplace approximation is in error in these regions of high sparsity. Yet, that the approximation only begins to show error when sparsity is>90% is notable. Beyond LA Collapsed and the HMC implementations, VB Collapsed consistently demonstrates higher error in uncertainty quantification as compared to the other implementations. 

Finally, to provide context regarding the size of the differences in uncertainty quantification, we provide direct visualizations of posterior intervals for all four implementations in Figure <ref> and <ref>. These two simulations were chosen to highlight a case in which LA Collapsed is highly accurate (<ref>) in terms of uncertainty quantification and a case in which it differs from HMC estimates (<ref>).  Notably, visualization of posterior intervals consistently demonstrated that the posterior mean was centered symmetrically in the 95% credible regions. This symmetry suggested that our metric RMSE of standard deviations captures much of the discrepancies in uncertainty quantification without higher order moments. Additionally, for context we include a fifth implementation, PCLM (pseudo-count augmented linear model). The PCLM uses a pseudo-count based estimate ofηwhich ignores the multinomial count variation. Such approximations are common in the analysis of microbiome sequence count data <cit.>. Unsurprisingly, this PCLM implementation demonstrates substantially higher error rates than any of the other implementations (Figure <ref>). 










§ REAL DATA ANALYSIS

Crohn's Disease (CD) is a type of inflammatory bowel disease that has been linked to aberrant immune response to intestinal microbiota <cit.>. To demonstrate that LA Collapsed (from the R package stray) provides an accurate and efficient means of modeling real microbiome data, we reanalyzed a previously published study comparing microbial composition in the terminal ileum of subjects with CD to healthy controls. Only LA Collapsed could efficiently scale to this data size (49 taxa, 250 samples, 4 covariates). To allow us to compare to alternative implementations we randomly subset the data to contain 83 samples. On this subset HMC Uncollapsed and VB Collapsed repeatedly failed to run due to numerical instability. In addition, LA Collapsed produced posterior estimates nearly identical to HMC Collapsed but more than 1000 times faster. 



Using the four model implementations introduced in Section <ref>, a Bayesian multinomial logistic normal linear model was fit to investigate the relationship between bacterial composition and CD. For both the full data-set and the subset, our regression model was defined for thej-th sample by the covariate vector 
    x_j=[1,x_j(CD),x_j(Inflamed),x_j(Age)]^T
 wherex_j(CD)is a binary variable denoting whether thej-th sample was from a patient with CD or a healthy control,x_j(Inflamed)a binary variable denoting inflammation at time of sample collection,x_j(Age)denoting age of the subject, and the preceding 1 represents a constant intercept. A detailed description of our prior assumptions is given in Appendix <ref> and results of posterior predictive checks are shown in Figure <ref>.

Even though all four implementations were initialized identically, both the HMC Uncollapsed and VB Collapsed implementations repeatedly resulted in errors due to numerical instability. Thus only LA Collapsed and the HMC Collapsed implementations could fit this model for even the subset dataset. Whereas the HMC Collapsed model took approximately 30 minutes, LA Collapsed took only 3 seconds. Thus LA Collapsed is over 1000 times faster than HMC Collapsed on real data.  Additionally, posterior estimates ofΛproduced by both the HMC Collapsed and LA Collapsed implementations are nearly identical (Figure <ref>). These results demonstrate that in real data scenarios LA Collapsed can provide efficient and accurate posterior inference. 

By modeling the full dataset we found the centered log-ratio (CLR) coordinates corresponding to 12 genera to be associated with CD status (95% credible interval not covering 0; Figure <ref>). These results are in general agreement with prior analyses <cit.>. As in prior analyses, we find a substantial increase in the abundance of proteobacteria in CD versus healthy controls. Similarly, we find that the families Pasteurellaceae and Enterobacteriacaeae, Gemellaceaem, and Fusobacteriaceae are highly enriched and that the class Clostridia are depleted in CD. Notably, Fusobacteria has been independently suggested as a marker of IBD <cit.>. These findings serve to validate our results and build confidence in our methods. 

In contrast, our results differ from prior analyses of this data in certain respects. We find that the family Peptostreptococcaceae is likely decreased in CD versus healthy controls and we find no association for Veillonellaceae. Three factors support our results. First, our analysis accounts for  count variation and compositional constraints whereas prior analyses have not. Notably, the handing of count variation and compositional constraints can have substantial impact on conclusions in the analysis of sequence count data <cit.>. Second, Peptostreptococcaceae has been found to be decreased in CD based on the analysis of independent data <cit.>. Third, in visualizing the count data for Peptostreptococcaceae and Veillonellaceae (Figure <ref>) we find no visual difference in Veillonellaceae but a notable difference in Peptostreptococcaceae. Therefore, we conclude that our approach has revealed novel associations in this data and excluded potentially spurious conclusions.


§ CONCLUSION


In this work we have developed efficient inference methods for the analysis of multinomial logistic-normal linear models through the use of a marginally latent matrix-t Process (LTP) representation. Through the use of the CU sampler and Laplace approximation (referred to as LA Collapsed) we have demonstrated up to 5 orders of magnitude increases in computational efficiency compared to HMC while preserving accuracy of point estimation and uncertainty quantification.  Beyond the inference of multinomial logistic-normal linear models, we find that our approach generalizes to a larger class of models including non-linear regression (LTP) models and time-series models (GMDLMs). Importantly, the key computational advance we have introduced relates to the marginal representation of multinomial logistic-normal linear models which is identical to that of multinomial logistic-normal GMDLMs. For this reason we hypothesize that the LA Collapsed sampler should share similar gains in efficiency for these other multinomial logistic-normal models. Yet, the performance of our Laplace approximation under observation distributions beyond the log-ratio parameterized multinomial is more uncertain. Notably, we hypothesize that our results could generalize to other exponential family distributions parameterized by natural parameters due to the fact that such distributions are globally log-concave. Yet, we expect that there are other observation distributions where a Laplace approximation to the LTP form may be sub-optimal. Rather than resorting to MCMC, we suggest that methods of particle refinement of the initial Laplace approximation (e.g., parallel MCMC steps for each sample from the LTP form, or sequential importance resampling) may be more efficient. We believe such extensions are prime areas for future work. 

Until this point we have not considered the presence of unknown hyperparameters in the LTP form (i.e. we have consideredδas given). Yet for a number of Marginally LTP models we expect estimation of such hyperparameters will be of interest. For example, within the GMDLM model we anticipate researchers may want to allow the termsW_tto be subject to their own stochastic model. This would in turn require that some portion ofδis unknown. Unfortunately, the multitude of ways in which such unknown hyper-parameters may arise is beyond the scope of this work. However, we note that our early efforts at extending multinomial logistic-normal linear models to handle variance components with unknown scale suggests that such hyperparameters can be learned efficiently during MAP estimation ofη̂. Alternatively, when the number of unknown parameters is small, such parameters could be choosen by cross-validation as is often done for choosing unknown parameters in kernel functions when using Gaussian process models. Still, we believe that future work incorporating unknown hyper-parameters into the LTP form as required by specific Marginally LTP models would be impactful. 


§ ACKNOWLEDGEMENTS


 JDS was supported in part by the Duke University Medical Scientist Training Program (GM007171). JDS and LAD were supported in part by the Global Probiotics Council, a Searle Scholars Award, the Hartwell Foundation, an Alfred P. Sloan Research Fellowship, the Translational Research Institute through Cooperative Agreement NNX16AO69A, the Damon Runyon Cancer Research Foundation, the Hartwell Foundation, and NIH 1R01DK116187-01. SM and KR would like to acknowledge the support of grants NSF IIS-1546331, NSF DMS-1418261, NSF IIS-1320357, NSF DMS-1045153, NSF DMS1613261, and NSF DEB-1840223. The authors thank Rachel Silverman, Liat Shenhav, and Shaobo Han for their helpful comments. 


jasa

§ GENERALIZED CONJUGATE MULTIVARIATE LINEAR (GMCL) MODELS

Here we demonstrate that the GMCL models defined in Equations (<ref>)-(<ref>) are Marginally LTP models by deriving their collapsed (LTP) form. Additionally, we demonstrate that uncollapsing the LTP form can be done efficiently. 



 §.§ Derivation of Collapsed Form

Our goal in this section is to derive a collapsed (or marginalized) form of the GMCL model defined in Equations (<ref>)-(<ref>) that will have a posterior distribution of the formp(η|Y,X). To do this, note that Equations (<ref>)-(<ref>) can alternatively be written as 

    η   = Λ X+E^η   E^η∼ N(0, Σ, I_N) 
    Λ   = Θ + E^Λ   E^Λ∼  N(0, Σ, Γ) 
    Σ   ∼ IW(Ξ, υ).


We will first marginalizeΛin the above model by using the following affine transformation properties of the matrix normal distribution. Given anm ×mmatrixA, ap ×mmatrixB, andn ×rmatrixC, and a random matrixX∼N(M, U, V); then for a random matrixZ = A + BXCwe haveZ ∼N(A+BMC, BUB^T, C^TVC)<cit.>. We note that Equations (<ref>) and (<ref>) can be collapsed into the following form

    η   = Θ X + E^Λ X + E^η   E^η∼ N(0, Σ, I_N)    E^Λ∼  N(0, Σ, Γ) 
       = Θ X + E^*    E^⋆∼ N(0, Σ, I_N + X^T Γ X).

Thus we may rewrite Equations (<ref>)-(<ref>) as 

    η   = Θ X + E^*    E^⋆∼ N(0, Σ, I_N + X^T Γ X) 
    Σ   ∼ IW(Ξ, υ).

By using the definition of the matrix-t given in Section <ref> we can marginalize overΣin Equations (<ref>) and (<ref>) to get 

    η∼ T(υ, Θ X, Ξ, I_N+X^TΓ X).

Finally, incorporating equations (<ref>) and (<ref>) allow us to write the marginalized form of GMCL models,p(η|Y, X), as

    Y    ∼ f(π)  
    π   = ϕ^-1(η) 
    η   ∼ T(υ, B, K, A)

whereB=ΘX,K=Ξ, andA=I_N + X^TΓX. 



 §.§ Efficient Form for Uncollapsing

Here we demonstrate that for GMCL models, the conditional posteriorp(Λ, Σ|η, Y, X)can be computed and sampled efficiently. AsΛandΣare conditionally independent ofYgivenηin GMCL models, we may write 

    p(Λ, Σ|η, Y, X) = p(Λ, Σ|η, X) = p(Λ|Σ, η, X)  p(Σ|η, X).
 
The right hand side of the above equation represents the posterior of a multivariate conjugate linear model that can be sampled efficiently using the following relations <cit.>:

    υ_N    = υ+N 
    Γ_N    = (XX^T+Γ^-1)^-1
    Λ_N    = (η X^T+ΘΓ^-1)Γ_N 
    Ξ_N    = Ξ + (η - Λ_N X)(η - Λ_N X)^T + (Λ_N - Θ)Γ^-1(Λ_N - Θ)^T 
    
        p(Σ | η, X)    = IW(Ξ_N, υ_N)
    
        p(Λ | Σ, η, X)    = N(Λ_N, Σ, Γ_N).


§ GENERALIZED MULTIVARIATE DYNAMIC LINEAR MODEL (GMDLM)

Here we demonstrate that the GMDLM defined in Equations (<ref>)-(<ref>) is a marginally LTP model and show how the theory and tools we develop improves inference in this model class. Additionally we provide a recursive procedure for uncollapsing an LTP to a GMDLM (i.e., an efficient recursive approach to samplingp(Θ, Σ|η, Y). 



 §.§ Derivation of Collapsed Form

Here we derive the marginal distribution ofηin terms of quantitiesF_t,G_t,W_t,Σ,M_0andC_0. As all densities involved are multivariate or matrix-variate normal, the result must also be multivariate or matrix-variate normal and thus fully described by the mean and covariance ofη. To derive the mean and covariance we first derive a useful alternative representation ofη^T_twith respect toΘ_t-k-1for some positive integerk < t. 

Substituting Equation (<ref>) into Equation (<ref>) allowsη_t^Tbe expressed with respect toΘ_t-1as 

    η^T_t = F_t^T G_t Θ_t-1 + F_t^T Ω_t + ν_t^T.

Repeated substitution ofΘ_t-kleads to the following form forη_t^Tin terms ofΘ_t-k-1η_t^T = F_t^T𝒢_t:t-kΘ_t-k-1 + F_t^TΩ_t + ∑_ℓ=t^t-k-1 F_t^T𝒢_t:ℓΩ_ℓ-1 + ν_t^T

where𝒢_t:t-kis shorthand forG_tG_t-1⋯G_t-k.

Note that ifX∼N(0, I, I),AA^T = UandBB^T=Vthen forZ = M+ AXB^Twe have the marginal distribution ofZ(e.g. marginalizing overX) given byZ ∼N(M, U, V). This identity in combination with (<ref>) allows us to marginalize over the random variablesΩ_t, …, Ω_1, ν_tin Equation (<ref>) giving

    η^T_t ∼ N(F_t^T𝒢_t:1M_0, γ_t + F_t^T[
        W_t + ∑_ℓ=t^2𝒢_t:ℓW_ℓ-1𝒢^T_ℓ:t + 𝒢_t:1C_0𝒢_1:t^T
        ]F_t, Σ) .


Next we calculateCov(η_t^T, η_t-k^T). In parallel to Equation (<ref>) we may writeη_t-k^Tas 

    η_t-k^T = F_t-k^TG_t-kΘ_t-k-1 + F_t-k^TΩ_t-k + ν_t-k^T.

Using Equation (<ref>) and (<ref>) along with the factCov(AX_1 + BX_2, Y) = A   Cov(X_1,Y) + B   Cov(X_2, Y)and thatCov(Θ_s, ν_ℓ) = Cov(Θ_s, Ω_ℓ) = Cov(Ω_ℓ, ν_s)=0for allsandℓ, we can write

    Cov(η_t^T, η_t-k^T)    = F_t^T𝒢_t:t-kVar(θ_t-k-1)G_t-k^TF_t-k + F_t^T𝒢_t:t-k+1Var(Ω_t-k)F_t-k

whereVar(Θ_t-k-1)can be written recursively as 

    Var(Θ_t-k-1) = G_t-k-1Var(Θ_t-k-2)G_t-k-1^T + Var(Ω_t-k-1)

whereVar(Ω_t-k-1) = Σ⊗W_t-k-1.  Combining this recursive form with equation (<ref>) gives 

    Cov(η_t^T, η_t-k^T) = F_t^T𝒢_t:t-k+1 (Σ⊗ W_t-k)F_t-k +
        ∑_ℓ=t-k^2 F_t^T𝒢_t:ℓ(Σ⊗ W_ℓ-1)G^T_ℓ:t-kF_t-k
    
        + F_t^T𝒢_t:1(Σ⊗ C_0)G^T_1:t-kF_t-k.


Together Equations (<ref>) and (<ref>) characterize the marginal distribution ofη_t^Tin terms ofF_t,G_t,W_t,Σ,M_0andC_0. Noting that ifX∼N(M, U, V)thenX^T∼N(M^T, V, U),  if follows that

    η   ∼ N(B, Σ, A) 
    
        B    = [                ; α_1   ⋯ α_t   ⋯ α_T;                     ]
    α_t    = (F_t^T𝒢_t:1M_0)^T 
    
        A_t,t-k   = 
          γ_t + F_t^T[
            W_t + ∑_ℓ=t^2𝒢_t:ℓW_ℓ-1𝒢^T_ℓ:t + 𝒢_t:1C_0𝒢_1:t^T
            ]F_t  if  k=0
    
            F_t^T[𝒢_t:t-k+1  W_t-k +
        ∑_ℓ=t-k^2 𝒢_t:ℓ W_ℓ-1G^T_ℓ:t-k
        + 𝒢_t:1 C_0G^T_1:t-k]F_t-k if  k>0

Finally, using the  marginalization property of the matrix normal and the inverse Wishart used in our definition of the matrix-t distribution and incorporating Equations (<ref>), (<ref>) and (<ref>) it follows

    Y    ∼ f(π)
    π   = ϕ^-1(η) 
    η   ∼ T(υ, B, Ξ, A).

This result demonstrates that the GMDLM model is a marginally LTP model. 




 §.§ Efficient Form for Uncollapsing

Here we provide an efficient means of sampling from the conditional densityp(Θ, Σ|η, Y)for the GMDLM. First we recognize thatΘis conditionally independent ofYgivenη. Therefore, our task simplifies to sampling fromp(Θ, Σ|η). The problem is identical to the standard filtering and simulation smoothing problem solved by <cit.>. Again, the problem is defined by the following model (which we will refer to as the MDLM model) 

    η_t^T    = F_t^TΘ_t + ν_t^T,   ν_t ∼ N(0, γ_t Σ) 
    Θ_t    = G_t Θ_t-1 + Ω_t,   Ω_t ∼ N(0, W_t, Σ) 
    Θ_0    ∼ N(M_0, C_0, Σ)
    Σ   ∼ IW(Ξ, υ).


Following <cit.>, below we restate the filtering and retrospective recursions needed to sample fromp(Θ, Σ|η). Note that all densities in this subsection are conditional on the parametersF_t,G_t,W_t,Σ,M_0andC_0but that this dependence has been suppressed for notational simplicity. Let us introduceυ_tandΞ_tas filtering parameters at stept. Further, we defineυ_0=υandΞ_0=Ξ. As a final piece of notation we introduceH^T_tas a shorthand for the set{η^T_t, …, η^T_1}§.§.§ Filtering Recursions for MDLM Model

(1) Posterior att-1:

    p(Σ| H^T_t-1)    ∼ IW(Ξ_t-1, υ_t-1) 
    
        p(Θ_t-1|Σ, H^T_t-1)    ∼ N(M_t-1, C_t-1, Σ)

(2) Prior att:

    A_t    = G_tM_t-1
    
        R_t    = G_tC_t-1G_t^T+W_t 
    
        p(Σ| H^T_t-1)    ∼ IW(Ξ_t-1, υ_t-1) 
    
        p(Θ_t |Σ, H^T_t-1)    ∼ N(A_t, R_t, Σ)

(3) One-step ahead forecast att:

    f_t^T    = F_t^TA_t 
    
        q_t    = γ_t + F_t^TR_tF_t 
    
        p(Σ| H^T_t-1)    ∼ IW(Ξ_t-1, υ_t-1) 
    
        p(η_t |Σ, H^T_t-1)    ∼ N(f_t, q_tΣ)

(4) Posterior att:

    e_t^T    = η_t^T-f_t^T  
    
        S_t    = R_tF_t/q_t
    
        M_t    = A_t + S_te_t^T 
     
        C_t    = R_t - q_tS_tS_t^T 
    υ_t    = υ_t-1 + 1 
    Ξ_t    = Ξ_t-1 + e_t e_t^T/q_t
    
        p(Σ| H^T_t-1)    ∼ IW(Ξ_t, υ_t) 
    
        p(Θ_t |Σ, H^T_t)    ∼ N(m_t, C_t, Σ)

Equation (<ref>) differs slightly from the presentation in <cit.> as the parameterization of the inverse-Wishart we employ throughout this paper differs from that source. Throughout this work we use the following parmeterization for a random matrixΣ∼IW(Ξ, υ):





    p(Σ) ∝| Σ|^-(P + υ + 1)/2exp( -1/2tr( ΞΣ^-1) ).


  §.§.§ Simulation Smoothing Recursion

The recursions provided here follow directly from <cit.>

(1) SampleΣ∼IW(Ξ_T, υ_T)and thenΘ_T ∼N(M_t, C_t, Σ). 

(2) For each timetfromT-1to0, samplep(Θ_t | Θ_t+1, H_T^T) ∼N(M_t^*, C_t^*, Σ)where 

    Z_t   =C_tG_t+1^TR_t+1^-1
    
    M_t^*    = M_t + Z_t(Θ_t+1-a_t+1) 
    
    C_t^*    = C_t - Z_t R_t+1Z_t^T.


§ GRADIENT AND HESSIAN CALCULATIONS FOR THE MATRIX-T DISTRIBUTION

Here we are concerned with calculating the gradient and Hessian of  

    log p(η| Y) ∝ -υ + N + P-1/2log| I_P + K^-1(η - B)A^-1(η-B)^T|.
 
LettingS = I_P + K^-1(η- B)A^-1(η-B)^Twe concern ourselves with calculating the quantitiesd log|S|/d(⃗η^T)andd log|S|/(⃗dη)(⃗dη)^T. We will use the following identity from matrix calculus  <cit.>: 

    dlog |S|    = Tr(S^-1 dS)  
    
    dS    = d(I_P + K^-1(η - B)A^-1(η-B)^T) 
       = d(K^-1(η A^-1η^T  -η A^-1 B^T - B A^-1η^T)) 
        = K^-1(dη A^-1η^T + η A^-1 dη^T - dη A^-1B^T-B A^-1 dη^T) 
        = K^-1(dη(A^-1η^T-A^-1B^T) + (η A^-1 - B A^-1)dη^T) 
        = K^-1(dη C + C^T dη^T)

where in the last line we have defined theN ×PmatrixC = A^-1(η^T - B^T). Further simplifying and using the identitiesTr(A) = Tr(A^T)andTr(AB)= Tr(BA)for matricesAandBwe get 

    dlog|S|    = Tr(S^-1K^-1(dη C + C^T dη^T)) 
       = Tr(S^-1K^-1dη C) + Tr(S^-1K^-1C^T dη^T) 
       = Tr(CS^-1 K^-1 dη) + Tr(CK^-1S^-Tdη) 
       = Tr(C(S^-1K^-1 + K^-1S^-T)dη) 
       =(⃗[C(S^-1K^-1+K^-1S^-T)]^T)^T(⃗dη) 
    
    dlog|S|   =(⃗(S^-1K^-1+K^-1S^-T)C^T)^T(⃗dη) 
    dlog|S|/(⃗dη)   = (⃗(S^-1K^-1+K^-1S^-T)C^T)^T


The HessianH=d^2log|S|/(⃗dη)(⃗dη)^Tcan then be calculated from equation (<ref>) by taking the differential again and manipulating the result into the following canonical formd^2log|S| = (⃗dη)^TH(⃗dη). In particular we make use of the following identities(⃗ABC) = (C^T⊗A)(⃗B)andd(S^-1) = -S^-1dS S^-1. We also make use of the vec-transposition matrix defined byT_m,n(⃗A) = (⃗A^T)whereAis anm×nmatrix andT_m,nis anmn×mnpermutation matrix. The vec-transposition matrix also satisfies the following propertiesT_m,n = T^T_n,m = T^-1_n,m. 

    d^2log|S|    = (⃗(S^-1K^-1+K^-1S^-T)dC^T+d(S^-1)K^-1C^T + K^-1d(S^-T)C^T)^T(⃗dη) 
       = [(⃗(S^-1K^-1+K^-1S^-T)dC^T)^T+(⃗d(S^-1)K^-1C^T)^T + (⃗K^-1d(S^-T)C^T)^T](⃗dη) 
       = [#1 + #2 + #3 ](⃗dη) 
    #1     = (⃗(S^-1K^-1+K^-1S^-T)dη A^-1)^T
        = ((A^-1⊗(S^-1K^-1+K^-1S^-T)) (⃗dη))^T 
        = (⃗dη)^T(A^-1⊗ (S^-1K^-1+K^-1S^-T))^T 
    #2    = -(⃗S^-1dSS^-1K^-1C^T)^T 
       =-((CK^-1S^-T⊗ S^-1)(⃗dS))^T 
       =-(⃗dS)^T(S^-1K^-1C^T ⊗ S^-T) 
    (⃗dS)^T    = (⃗K^-1(dη C + C^T dη^T))^T 
       =(⃗K^-1dη C)^T + (⃗K^-1C^Tdη^T)^T 
       =((C^T⊗ K^-1)(⃗dη))^T + ((I_D-1⊗ K^-1C^T)(⃗dη^T))^T 
       =(⃗dη)^T(C⊗ K^-1) + (⃗dη^T)^T(I_P⊗ CK^-1) 
    #2   = [-(⃗dη)^T(C⊗ K^-1) - (⃗dη^T)^T(I_P⊗ CK^-1)](S^-1K^-1C^T⊗ S^-T) 
       = -(⃗dη)^T(CS^-1K^-1C^T⊗ K^-1 S^-T) - (⃗dη^T)^T(S^-1K^-1C^T⊗ CK^-1 S^-T) 
       = -(⃗dη)^T(CS^-1K^-1C^T⊗ K^-1 S^-T) - (⃗dη)^TT_N, P(S^-1K^-1C^T⊗ CK^-1 S^-T) 
    #3    = (⃗K^-1d(S^-T)C^T)^T
       =-(⃗K^-1S^-TdS^TS^-TC^T) 
       =-((CS^-1⊗ K^-1S^-T)(⃗dS^T))^T 
       =-(⃗dS^T)^T(S^-TC^T⊗ S^-1K^-1) 
    (⃗dS^T)^T   = (⃗(dη C + C^T dη^T)^TK^-1)^T 
       = ((K^-1C^T⊗ I_P)(⃗dη))^T + ((K^-1⊗ C^T)(⃗dη^T))^T 
    #3    = [-(⃗dη)^T(CK^-1⊗ I_P)-(⃗dη^T)^T(K^-1⊗ C)](S^-TC^T ⊗ S^-1K^-1) 
       = -(⃗dη)^T(CK^-1S^-TC^T⊗ S^-1K^-1)-(⃗dη)^TT_N,D-1(K^-1S^-TC^T⊗ CS^-1K^-1) 
    
    d^2log|S|    = (⃗dη)^T [(A^-1⊗ (S^-1K^-1 + K^-1S^-T))^T - (CS^-1K^-1C^T⊗ K^-1S^-T)  
            -(CK^-1S^-TC^T) ⊗ S^-1K^-1)  
            -T_N, D-1((S^-1K^-1C^T⊗ CK^-1S^-T)+(K^-1S^-TC^T⊗ CS^-1K^-1))] (⃗dη)

Summarizing the above results together we obtain

    S    = I_P + K^-1(η-B)A^-1(η-B)^T 
    
    C    = A^-1(η -B)^T
    
    R    = S^-1K^-1
    dlog |S|/(⃗dη)   = (⃗(R+R^T)C^T)^T 
    
    L    = (CRC^T⊗ R^T) 
    d^2log|S|/(⃗dη)(⃗dη)^T   = (A^-1⊗ (R+R^T)) - (L + L^T) - T_N, D-1[(RC^T⊗ CR^T)+(R^TC^T⊗ C R)].

The following computational trick makes evaluation of this Hessian far more computationally efficient. We may quickly calculateT_m,nX=X^*for anm×mmatrixXhaving already computedXby noting that fori∈1 …mandj∈1…nwe can writeX^*_(i-1)n+j, · = X_(j-1)m+i, ·whereX_l,·denotes thel-th row of the matrixX. 


§ GRADIENTS AND HESSIANS FOR  THE LOG-RATIO PARAMETERIZED MULTINOMIAL

Unfortunately we cannot provide a general form for the gradient and Hessian of all possible likelihoodsf(Y |ϕ^-1(η)). For the purposes of this article, here we derive the gradient and Hessian for the case wherefis multinomial andϕ^-1is the inverse ALR transform:

    ∑_j logMultinomial(Y_· j|ALR_D^-1(η_· j))

which for notational simplicity we will refer to asg. Thus our goal is to find efficient forms for calculatingg,dg/d(⃗η)andd^2g/d(⃗η)d(⃗η)^T. Using the fact thatlogMultinomial(Y_·j |π_·j) ∝Y_1jlogπ_1j+…+Y_Djlogπ_Djand Equation (<ref>) we can write 

    g = ∑_j=1^N(∑_i=1^D-1η_ijY_ij-n_jlog(1+∑_i=1^D-1e^η_ij) ).

Differentiating with respect toη_ijgives 

    dg/dη_ij =  Y_ij - n_j e^η_ij/1+∑_i e^η_ij.

Differentiating again with respect toη_kℓgives

    d^2g/dη_ijdη_kℓ = 
    
    -n_j( e^η_ij/1+∑_i e^η_ij - e^2η_ij/(1+∑_i e^η_ij)^2)    if ℓ=j, i=k
    
    n_j ( e^η_ije^η_kj/(1+∑_i e^η_ij)^2)    if ℓ=j, i≠ k
    
    0    if ℓ≠ j.

These results directly imply the following matrix forms. 

    O    = expη
    
        m    = 1_N + O^T 1_D-1
    ρ   = (⃗O)⊘(⃗1_D-1m^T) 
    
        n    = 1^T_DY 
    
        g    = - (⃗η)^T(⃗Y_/D·) - n⊙log(m) 
    dg/d(⃗η)   = ((⃗Y_/D·) - (⃗1_Dn)⊙ρ)^T 
    
        W^(j)   = n_j(ρ_(j)ρ_(j)^T-diag(ρ_(j))) 
    d^2g/d(⃗η)d(⃗η)^T   = diag(W^(1), …, W^(N))

whereexpXandlogXrefers to the element-wise exponentiation and logarithm of a matrixX,⊙and⊘refer to element-wise product and division respectively,Y_/D·refers to the firstD-1rows of the matrixY,ρ_(j)denotes elements(j-1)(D-1)+1toj(D-1)in the vectorρ, anddiag(X_1, …, X_D)refers to a block diagonal matrix where thei-th block isX_i. 



§ ACCELERATED MATRIX-T GRADIENTS VIA SYLVESTER'S DETERMINANT IDENTITY

Sylvester's determinant identity states that for matricesAandBof sizem×nandn×mrespectively,|I_m + AB| = |I_n+BA|. This relationship can be used to speed up calculation of the log-likelihood and gradient of the matrix-tdistribution whenN < Pas the the log determinant or inverse of the matrixScan dominate computational time. To take advantage of this speed up we note that we can replace the relations given in Appendix <ref> with 

    S    = I_N + A^-1(η - B)^TK^-1 (η - B) 
    
        C    = K^-1(η - B)
    
        R    = S^-1A^-1
    dlog|S|/d(η)   = (C(R+R^T))^T.

While this result can greatly accelerate inference for matrix-t gradients whenP ≫N, this result provides only minimal improvement for calculating the corresponding Hessian terms. Therefore we suggest that, for simplicity, the Hessian form provided in Appendix <ref> be used even ifP ≫N. 



§ SIMULATIONS AND MODEL FITTING

To compare the performance of the multiple multinomial logistic-normal linear model implementations described in Section <ref> over a range of sample sizes (N), observation dimensions (D), and covariate dimensions (Q), we created the following simulation scheme. For each evaluated triple (N,D,Q), three simulated data-sets were created based on the multinomial logistic-normal linear model
with the following specified likelihood:

    Y_· j   = ALR_D^-1(η_· j) 
    η_· j   = N(Λ X_· j, Σ) 
    Λ   ∼ N(Θ, Σ, Γ) 
    Σ   ∼ IW(Ξ, υ)

withυ=D+10,Ξ=I_D-1. AdditionallyX,Λ, andΣwere simulated as 

    Λ   ∼ N(0, I, I) 
    Σ   ∼ IW(Ξ, υ) 
    
        X    ∼ N(0, I, I).

The percent of zero counts naturally increased with largeDor largeQrelative to other parameters. We took advantage of this behavior to study the performance of all implementations in sparse data regimes. 


All implementations were compiled and run using gcc version 6.2.0, R version 3.4.2, and Intel(R) Math Kernel Library version 2019 where possible. All replicates of the simulated count data were supplied to the various implementations 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.


§ PRIORS FOR CROHN'S DISEASE DATA


 §.§ Data Preprocessing

Sequence count data was obtained from the R package MicrobeDS (github.com/twbattaglia/MicrobeDS). Only samples from the terminal illeum from healthy donors and patients with Crohn's Disease, who had no recent history of steroids, antibiotics, or biologics were included in analysis. Samples with a sequencing depth below 5000 counts were excluded from analyses. Only families seen with at least 3 counts in at least 10% of samples were retained for subsequent analyses. 



 §.§ Priors

The regression model required that 4 hyper-parametersΓ,Θ,Ξ, andυbe specified.  We setΘto aD×Qmatrix of zeros representing our prior assumption that, on average, there was no association between each covarariate and microbial composition. We specifiedΓ=I_Qto constrain associations between microbial composition and covariates to remain small. We specifiedυ=D+3andΞ_ij = (υ-D)ifi==jandΞ_ij = (υ-D)/2ifi ≠jto reflect our weak prior assumption that the log absolute abundance of each taxa is uncorrelated <cit.>. 



§ SUPPLEMENTARY FIGURES