Latent Network Estimation and Variable Selection for Compositional Data Via Variational EM

2.1. Prior on Covariate Effects B

Here, we describe the prior on the covariate effects, which enables selection of the important associations between X and other potentially related factors M. We consider the effects of the covariates M on each column of Z separately, which means that we will be able to update the columns of B independently of each other. Here B is a q × p matrix, where each column of B represents the vector of regression coefficients for the q covariates of M on the jth column of Z. We use a spike-and-slab prior on each element of the matrix B, which shrinks the effects of features that do not influence Z to zero. Remember that we are looking at the columns of Z one at a time, and can thus say that any entry from the jth column, Z_i,j, comes from a $\text{Normal}\left(B_{0j}+M_i{B}_j,{\sigma }_j^{*}\right)$, where ${\sigma }_j^{*}$ is the standard deviation of the jth column of Z, found by using the properties of the multivariate normal distribution shown in Equation (1). The prior on B is as follows:

for j = 1, …, p and k = 1, …, q, and with δ₀ a point mass at 0, indicating that when γ_k,j is 0, B_k,j is exactly 0. Here, θ_γj is the probability of a variable being relevant in B_j. Notice that the mixture prior (2) allows each variable to be relevant for individual responses (Richardson, Bottolo, and Rosenthal 2010; Stingo et al. 2010), as opposed to spike-and-slab constructions that select variables as relevant to either all or none of the responses (Brown, Vannucci, and Fearn 1998). We also put a non-informative prior on each element of B₀, that is, B_0j ∝ 1.

2.2. Prior on Precision Matrix Ω

Next we introduce the prior on the precision matrix Ω, which allows us to learn a sparse association network. We consider the prior of Wang (2015) in the formulation proposed by Li and McCormick (2019)

where ν₀ and ν₁ are fixed standard deviations, that assume small and large values, respectively, δ_i,j is a latent variable indicating whether or not an edge is present between nodes i and j, and τ is a scaling parameter, with a hyperprior Gamma(a_τ, b_τ) that allows us to adaptively learn the standard deviations. The original prior of Wang (2015) is obtained by setting τ = 1. Additional complexity can be added to the prior on τ to include existing knowledge about variable associations, as shown in Li and McCormick (2019). The mixture of normals on the off-diagonal precision matrix entries enables the selection of interactions, represented by edges in a network, since nonzero precision matrix entries reflect conditional dependence relationships (Dempster 1972). Here, entries reflecting conditional independence relations do not equal exactly zero, but get shrunk to close to zero. The diagonal entries are drawn from a common exponential prior. The final term in Equation (3) expresses a constraint to the space of positive definite matrices M⁺. This prior is particularly advantageous in our model, as it allows for efficient estimation via the EM algorithm and leads to less bias in estimation of the off-diagonal precision matrix elements than the graphical LASSO, as shown by Li and McCormick (2019).

We complete the modeling of the precision matrix Ω by setting the prior on the graph structure, assuming independent Bernoulli distributions on the inclusion of each edge as follows:

3.1. VI Step

Here, we use a VI step to estimate the regression coefficients B by updating the free parameters μ_k,j, σ_k,j, and ϕ_k,j, where μ_k,j and σ_k,j are the mean and variance, respectively, of B_k,j when γ_k,j 1, and ϕ = k,j is interpreted as the probability of γ_k,j = 1, resulting in B_k,j ≠ 0. Following the work of Titsias and Lázaro-Gredilla (2011) and Carbonetto and Stephens (2012) the free parameters can then be updated as

Updating each column of B can then be done independently and, when resources are available, these updates can be done in parallel. While updating a column of B_j, each component of μ_j, σ_j, ϕ_j is updated given all other components. This component-wise update of μ_j, σ_j, ϕ_j is repeated until ELBO(μ_j, σ_j, ϕ_j) has converged. Once all μ and ϕ have been updated, the individual elements of B are assigned $\mathbb{E}\left(B_{k,j}\right)={\mu }_{k,j}{\varphi }_{k,j}$. Once the SINC algorithm has converged, as common in variational spike-and-slab literature (Huang, Wang, and Liang 2016; Miao et al. 2020), we set B_k,j = μ_k,j if ϕ_k,j > 0.5 and B_k,j = 0 if ϕ_k,j ≤ 0.5. The threshold of 0.5 equivalent to selecting the median model of Barbieri and Berger (2004) and can be adjusted to include or exclude more covariates, but the threshold of 0.5 is the most commonly used.

3.2. E Step

In this step, we focus on updates to the edge inclusion parameter δ_i,j. For the first step we take the expectation of the posterior distribution, treating δ as the latent variable. We define Q(Θ| Θ^(t), ξ^(t)) as $E_{\delta \mid {\Omega }^{(t)},\pi ,X}\left(\text{log\hspace{0.17em}}p\left(Z,\Omega ,B,B_0,\pi \mid X,M\right)\mid {\Omega }^{(t)},\pi ,X\right)$. Following the results shown in Li and McCormick (2019), the E step can be broken into two steps:

where a_i,j = p(ω_i,j | δ_i,j = 1)π and b_i,j = p(ω_i,j | δ_i,j = 0)(1 − π), and (i, j) is the (i, j)th entry of the precision matrix, where i and j ∈ {1, …, P}.

3.3. M Step

The remainder of the unknown parameters can be found by maximizing the posterior distribution with regards to each of the parameters we are interested in. Here, we first update the column-wise centering parameters B_0j independently as

Next, we update the precision matrix, Ω. Following Li and McCormick (2019) and Wang (2015) the conditional distribution of each column of Ω can be found in closed form. For this, let

Then, the conditional distributions are

We can then do a column-by-column update as

The point estimates of θ_γ and π are also updated as follows:

If using the adaptive scale parameter, τ, then an additional update is done by setting

Finally, the matrix of latent variables can be estimated by finding a point estimate for each entry of the matrix. This is done by updating each row of the matrix independent of the others. As shown in Yang, Chen, and Chen (2017), the objective function to optimize with respect to Z is

where $\tilde{\Gamma }$ is the log-gamma function, and $s\left(x_i\right)={\sum }_{j=1}^p{x}_{i,j}$. To accomplish optimization of each Z_i we use the limited-memory quasi-Newton (L-BFGS) algorithm, which is a quasi-newton gradient descent method that makes use of the inverse gradient to direct where to search through the variable space.

For posterior inference, we iterate through the VI step, which iterates between updating ϕ_k,j, μ_k,j, and σ_k,j, and the E and M steps, which update Ω one column at a time, until the algorithm has converged. For both the VI and the M steps, we run each of those steps until the respective parameter estimates have converged. We determine the algorithm to have converged if the ELBO changes by less than a predefined tolerance from one iteration to the next. To obtain a selected network and set of covariates based on these posterior estimates, we select edges i, j with $p_{i,j}^{*}$ in Equation (9) ≥ 0.50, and covariate associations with ${\varphi }_{k,j}^{*}$ in Equation (8) ≥ 0.50. In practice, both the $p_{i,j}^{*}$ and ϕ_k,j values, which reflect the posterior probabilities for the selection of edges and covariates, tend to converge to values close to 0 or close to 1. A similar trend has been noted by Kook et al. (2021), who reported that the variational parameters for the marginal posterior probabilities of inclusion tended to become more widely separated as the algorithm converges.

4.1. Simulation Setup

Simulated data were generated with the following steps. First, the covariates, M, were generated from a normal distribution MVNorm(0, I_p), and subsequently scaled. The values of the regression coefficients B, related to M, were then sampled. Each element of B, B_k,j, was assigned either a random value between [−1, −0.5] with probability 0.1, a value in [0.5, 1] with probability 0.1, or else 0. Each B_0j was then sampled from the interval [6, 8] with probability 0.2, and from [2, 4] with probability 0.8. This allowed for some variables to have larger counts and others to be sparser, as common in microbiome data. Simulated counts were then sampled by first drawing Z from MVNorm(MB + B₀, Ω), with Ω as described below, and then assigning α as exp(Z). Finally, h_i was sampled as a random draw from Dirichlet(α_i), and X_i drawn from a Multinomial(h_i, nint(N_i)), with N_i generated from a Normal(3000, 250), allowing for samples to have different numbers of total counts, where nint() represents the nearest integer function. We set p = 100, q = 50 and n = 300.

To explore performance for a range of possible network structures, we simulated a variety of configurations. These networks were created using the R package huge (Jiang et al. 2019). For this simulation, we used a band, cluster, hub, and random graph structure. An example of what these networks look like can be seen in Figure 1. Band graphs and random graphs are common test cases for network learning, while the hub and cluster graphs capture some aspects of biological networks, such as highly connected nodes and community structure (Girvan and Newman 2002). The probability of an edge in the network was set to 0.025 for the random graph and 0.30 for each cluster in the cluster graph. The bandwidth in the graph was set to 3 and the number of hubs in the hub graph was set to 3. The precision matrix Ω used to generate the simulated data was also constructed using the function huge.generate of the R package huge Jiang et al. (2019). Parameters v and u of huge.generate, which control the off-diagonal elements of the precision matrix and magnitude of the partial correlations, were set to 1 and 0.0001, respectively.

The results we report below for our proposed model were obtained with the following hyperparameter settings. Fairly non-informative priors were set by choosing a_γ = b_γ = 2 in the prior probability of inclusion (2) of each covariate. The same setting was used for the hyperparameters a_π and b_π in Equation (4), which determines the prior probability of an edge being included. The standard deviation of the prior on the selected regression coefficients, ν_B in Equation (2), was set to 1. Following guidelines given by Wang (2015) and Li and McCormick (2019), we set ν₁ = 10, the standard deviation in prior (3) on the off-diagonal precision matrix entries corresponding to selected edges, and fit the model across a grid of values, ranging from 0.0001 to 0.1, of ν₀, the prior standard deviation of off-diagonal elements of the precision matrix corresponding to nonselected edges. We then chose the final model by using the ν₀ value that gave sparsity closest to 0.10. This sparsity level was selected arbitrarily, and did not result in any specific advantage for our method, as none of the simulation networks had sparsity equal to 0.10. Finally, the rate parameter λ of Equation (3), which appears in the prior on diagonal elements of the precision matrix, was set to 150. We also compare the model when the scaling parameter, τ, is learned, and use a Gamma(2,2) prior to do so. We comment on the sensitivity of the results to parameter choices in Section 4.3.

4.2. Simulation Results

We compare the performance of our method to several existing alternative approaches in terms of accuracy in network estimation and covariate selection. For comparison, we used mLDM (Yang, Chen, and Chen 2017), which is specifically designed for estimating networks of compositional data while controlling for covariates, and mSSL-DPE (Deshpande, Ročková, and George 2019), which we applied to the centered log ratio (CLR) transformed version of the simulated count data, a common method to account for the compositionality (Fang et al. 2017; Kurtz et al. 2015). For network estimation, we also considered SpiecEasi (Kurtz et al. 2015), which applies graphical LASSO to the CLR-transformed data. These methods were applied by using the default selection criteria in their respective R packages. To more precisely characterize factors contributing to the network estimation performance of the SINC method, we apply a version of SINC with the covariate effects B constrained to be 0. A comparison of the results from this constrained version of SINC to those of SpiecEasi reflects the performance advantages arising from differences in the network estimation procedure, while comparison to the full unconstrained SINC method provides quantitative insight on the benefit of simultaneous estimation of covariate effects on network recovery. Similarly, for the comparison of variable selection accuracy, we apply a constrained version of SINC with Ω fixed to the identity matrix. Comparison of the results from this approach to those of the full unconstrained SINC method illustrates the added value of accounting for the residual covariance in estimation of B.

We report results in terms of true positive rate (TPR), false positive rate (FPR), F1 score, and Matthew’s correlation coefficient (MCC). For edge selection, we also report the area under the curve (AUC). This was calculated, for the SINC method, over a grid of ν₀ values, and for the mLDM and mSSL-DPE methods by using the LASSO penalization parameter for the coefficients associated with the best selected graph, and then varying the graph penalization parameter over a grid of values. The AUC for SpiecEasi was calculated by varying the penalization parameter. Tables 1 and 2 show the results for network estimation and variable selection, respectively. From Table 1, we can see that mSSL-DPE and SpiecEasi are generally not competitive in terms of their performance, with low F1, MCC, and AUC values. This is likely because mSSL-DPE was not designed for compositional data, and SpiecEasi is not able to account for the effects of the covariates on the counts. We also see that mLDM performs better than the other two methods but is still outperformed by the proposed model, which does better in all F1, MCC, and AUC scores across all of the network structures except the hub structure. Finally, we see that when the proposed model does not control for additional covariates, the network estimation scores decrease and are comparable to the other methods. Across all methods, SINC while learning τ performed best in all network structures in terms of F1, MCC, and AUC. The performance metrics for SINC are pretty similar across all network types, though the best performance is achieved on the random graph, while the hub and cluster settings are more challenging. From Table 2, we can see that mSSL-DPE performs quite well in selecting the covariates of interest. In fact, its performance is very close to the proposed model, SINC, on all metrics. Similarly, the proposed model, while holding the estimated network and precision matrix fixed performs well for coefficient estimation, but does not do as well as mSSL-DPE and the full version of the proposed model. Additionally, we did not see any significant difference in performance in the full SINC models when τ is fixed or learned. SpiecEasi is not included in Table 2, as it is not able to select or adjust for relevant covariates, which is a limitation of the method.

We did not compare our model to MCMC approaches because of the computational complexity resulting from a lack of conjugacy. We did, however, experiment by using a Monte Carlo draw to update Ω at each iteration of SINC, and found that the point estimates of SINC without a Monte Carlo step were very close to the mean of the Monte Carlo draws.

4.3. Influence of Parameters

An advantage of using spike-and-slab priors for covariate and network edge selection, over penalized methods, is given by the flexible level of sparsity induced on the regression coefficients and the precision matrix entries. For example, Li and McCormick (2019) showed that $d_{i,j}^{*}$ from Equation (10) is comparable to the penalty parameter, λ, in the graphical LASSO (Dempster 1972). However, $d_{i,j}^{*}$ is unique to each edge and is adaptively learned from the data. In Figure 2, we show this advantage over penalized methods by plotting the estimated coefficients and precision matrix values, for a smaller simulation scenario with p = 10, q = 15, and n = 5000, using SINC (with τ fixed at 1) and the penalization based method mLDM, while varying the sparsity-inducing parameters. The top-left plot shows the estimated B coefficients by the proposed model when increasing the variance parameter ν_B of the spike-and-slab prior in Equation (2). The top-right plot shows the estimated B coefficients via mLDM when increasing the LASSO penalty parameter. The bottom-left plot shows the estimated off-diagonal values of the precision matrix Ω when increasing the variance parameter ν₀ in Equation (3) and the bottom-right plot shows the estimated off-diagonal estimates of the precision matrix when increasing the graphical LASSO penalty parameter. In all plots, red lines correspond to true associations in the simulated data, and black lines correspond to coefficients representing no underlying association. The flat trend in the red lines of the plots related to SINC shows that the estimated covariate effects and precision matrix entries corresponding to true associations are stable, while for mLDM, depicted at bottom, they get shrunk to zero as the penalty parameters increase.

We conclude this section by providing some comments on the sensitivity of the results to the choice of the hyperparameters. As shown in Figure 2, with sufficient data, the estimates of Ω and B are stable for increasing values of ν₀, in the prior of Equation (3), and ν_B, in the prior of Equation (2), respectively. These parameters, however, affect the sparsity of the selection. In particular, as ν₀ increases, holding all other parameters constant, the selected network becomes sparser. Similarly, as ν₁, which appears in the prior of Equation (3), increases, holding ν₀ constant, the network sparsity increases. In recent work using this type of mixture prior, Li and McCormick (2019) and Ročková and George (2014) have suggested holding ν₁ constant while varying ν₀. It should be apparent, then, that increasing ν_B while holding all other parameters constant decreases the number of selected coefficients, as increasing ν_B is analogous to increasing ν₁. The remainder of the hyperparameters influence sparsity as well, but to a lesser extent. For example, changing a_π and b_π in the prior given in Equation (4) to put more weight on larger values of π results in sparser networks. Similarly, selecting a_γ and b_γ in the prior of Equation (2) to reflect a stronger prior belief in larger θ_γ values results in an increase in the number of selected coefficients. Since π and θ_γ are both updated at each iteration of the SINC algorithm, selecting relatively non-informative priors, such as the ones used in the simulations of a_γ = b_γ = 2, allows the sparsity levels to be primarily controlled by ν₀ and ν_B. Alternative choices that would also be appropriate include a_γ = b_γ = 1, a more non-informative setting corresponding to a uniform prior on the unit interval, or a_γ = 1 and b_γ = p, which would more strongly favor sparsity, as discussed in Ročková and George (2014). We found that our variable selection results were not overly sensitive to the choice of these parameters. Increasing λ also increases the network sparsity because it changes the scale of the estimated Ω values by making them smaller. Appropriate λ values need to be selected based on the scale of the data that is being used.

5.1. MOMS-PI Results

Figure 3(a) shows the adjacency matrix of the microbial network inferred from the MOMS-PI data, with filled boxes representing selected edges, together with a plot of the number of edges for each OTU in 3(b). A network diagram of the inferred microbial network is shown in in Figure 3(c), with node sizes representing the degree of the nodes, so the larger a node, the more edges that node has with other nodes. In these plots, OTUs are grouped based on their phylum (Firmicutes, Actinobacteria, Bacteroidetes, Proteobacteria, Fusobacteria, and TM7). Looking at the adjacency matrix and network representation, we notice that Actinobacteria have few shared edges with Bacteroidetes, Proteobacteria and Fusobacteria, instead sharing the majority of their inter-phylum edges with Firmicutes, while the other phyla (Firmicutes, Bacteroidetes, Proteobacteria and Fusobacteria) show no trend in inter-phylum edges. We also notice that within the Firmicutes subnetwork, OTUs of the genus Lactobacillus (OTU 1 through 26 in the adjacency matrix) form their own subnetwork with very few inter-genus connections. Also, from the node degrees plot we can see that many Firmicutes have large numbers of edges. Indeed, when looking at the most connected nodes of the inferred microbial network, we found that 5 of the 6 most connected OTUs belong to the Firmicutes phylum.

Next, we show the inference on the microbe-cytokine association network. Figure 4 shows the adjacency matrix of the selected microbe-cytokine associations, with microbes colored based on phylum. We observe clear patterns of association, with both cytokines that show relationships to many OTUs, and OTUs that show relationships with several cytokines. In particular, we see that the cytokines IP-10, IL-1b, IL-17A, FGF basic, and IL-8 have the most associations with OTU abundances. When looking at which OTUs have the most associations with cytokines, we found that 6 of the 10 most connected are Lactobacillus. This is also seen in Figure 4, where many of the first 26 microbes (columns) have several cytokine associations. Lactobacillus has previously been shown to be largely influenced by cytokines (Valenti et al. 2018).

We also compare the results from SINC to those from SpiecEasi (Kurtz et al. 2015) and the B-constrained version of SINC, and found that the two methods that do not control for covariates shared 22 edges that were not selected by the proposed model. Two of these edges can be seen in Figure 5, which shows the network of a subset of three microbes and three cytokines estimated by SINC, as well as the network of the same subset of microbes estimated by SpiecEasi and a variant of SINC with the B coefficients not estimated. We hypothesize that SINC did not select the same edges as the other two methods, that is, the edge between OTUs 14 and 19 and the edge between OTUs 19 and 20, because edges selected by methods that do not control for cytokines may incorrectly determine an edge when microbe pairs have a mutual association with a cytokine. This can be seen in Figure 5, where OTUs 14,19, and 20 all have an association with IP-10. This illustrates the ability of our model to discover covariate effects and a sparse network accounting for these effects.