Let y i t i k = ( y i t i k 1 , · · · , y i t i k p ) T denote observations of some transformed abundance data of a microbiome with p taxa from subject i at time t_k (1 ≤ i ≤ m, 1 ≤ k ≤ n_i ) so thatkit is appropriate to assume that y i t i k ∼ N p ( µ , Σ ) , where µ = E ( y i t i k ) and Σ = V a r ( y i t i k ) . The precision matrix is defined as Ω = Σ⁻¹. Then, the n_ip vector y i = ( y i 1 T , ⋯ , y i n i T ) T represents all the observations on subject i, and vector y = ( y 1 T , ⋯ , y m T ) T represents the observations on all the m subjects with n = ∑ i = 1 m n i . For the correlations between the observations, we assume that the observations from different subjects are independent, i.e., c o v ( y i 1 t i k 1 , y i 2 t i k 2 ) = 0 p × p for i 1 ≠ i 2 , k 1 ≥ 1 , k 2 ≥ 1 . For observations from the 1 same subject, we assume c o v ( y i t i k 1 , y i t i k 2 ) = D H i k 1 k 2 D where D = diag(σ ₁,···,σ_p ) with σ 1 2 , ⋯ , σ p 2 diagonal elements of Σ, while H i k 1   k 2 is the correlation matrix between Y i t i k 1 and Y i t i k 2 for which the following form is assumed:

The symbol in (Equation 1) stands for the Hadamard product of matrices Φ i k 1 k 2 and R. Here, R is the correlation matrix with respect to covariance matrix Σ, while matrix Φ i k 1 k 2 = ( ϕ i k 1 k 2 ) p × p defines the dampening rates at which the components of H i k 1 k 2 decrease as time goes from t i k 1 to t i k 2 . For example, ( Φ i k 1 k 2 ) 12 is the dampening rate of correlation c o r ( Y i t i k 1 1 , Y i t i k 1 2 ) to correlation c o r ( Y i t i k 1 1 , Y i t i k 2 2 ) . Theoretically, dampening rates can vary from taxon to taxon and depend on the time points as long as the resulting matrix H_{ik 1} k ₂ is positive definite. However, in this paper, for subject i, we assume that the components of H i k 1 k 2 have the same dampening rate. Furthermore, they depend on time points ( t i k 1 , t i k 2 ) only through the distance between t i k 1 to t i k 2 , i.e., ϕ i k 1 k 2 = g i ( | t i k 1 − t i k 2 | ) for some decreasing function 0 ≤ g_i (·) ≤ 1. Motivated by studies on the longitudinal regression model (Diggle et al., 2002), we assume function g_i (·) has the form of exp ( − τ i | t i k 1 − t i k 2 | p ) . For p = 0, we have ϕ i k 1 k 2 = exp ( − τ i ) , which is referred to as the uniform correlation and can be used to model the spatial correlation. For example, specimens may be collected at different body sites from the same subjects, for which the uniform correlation seems to be a reasonable assumption. On the other hand, the cases of p > 0 can be used to model the irregularly spaced temporal correlation, which typically decreases as the time span | t i k 1 − t i k 2 | increases. In particular, functions exp ( τ i | t i k 1 − t k 2 | ) and exp ( τ i | t i k 1 − t i k 2 | 2 ) have been used in the marginal regression model for low-dimensional longitudinal data. Here, the parameters τ_i s, which are referred to as autocorrelation parameters, measure the dampening rates that are shared by all the components of y i t , ( i = 1 , ⋯ , m ) .

Without loss of generality, we always employ the correlation function exp ( τ i | t i k 1 − t i k 2 | ) in the following and assume that the observations have been centered so that μ = 0. Let Σ _i denote the covariance matrix of the observation vector y _i . The density function of y is then given by

Since the number of unknown parameters in Ω is much larger than the sample size in the context of the gut microbiome, the maximum likelihood estimate of Ω is unidentifiable, and sparsity is typically assumed in the literature. To this end, penalized maximum likelihood estimation (MLE) is usually adopted, e.g., the SPIEC-EASI model in (Kurtz et al., 2015). The SPIEC-EASI pipeline employes the L ₁-penalty to achieve the sparsity of Ω for cross-sectional observations. Here, we adopt the same strategy for longitudinal data and use the minimizer of the following L ₁-penalized negative log-likelihood function as the estimate of network Ω

We refer to model (Equations 2, 3) the stationary Gaussian graphical model (SGGM). Here, stationarity stems from the fact that the same network Ω is shared by all the subjects and at all time points. If the data are independent observations, then (3) can be solved by the graphical LASSO algorithm (Friedman et al., 2008) or the neighborhood method (Meinshansen and Bühlmann, 2006), and SGGM is just reduced to the SPIEC-EASI model. However, since the data are longitudinal and can correlate to each other, the performance of the SPIEC-EASI pipeline is not guaranteed when solving (3).

Notably, in model (2, 3), we assume the same correlation dampening rate τ_i for all the taxa in the microbiome of subject i. This assumption is motivated by the characteristics of the gut microbiome, where the taxa are typically influenced by the same perturbation sources, such as diet change and disease development. However, even if this assumption is violated and different taxa have different dampening rates, model (2, 3) can still be used as a working model for identifying the network structure, in which τ_i can be regarded as the mean dampening rate of the whole microbiome. In such cases, models (2, 3) still outperform the SPIEC-EASI pipeline, and the latter ignores the correlation structure of longitudinal data. We demonstrate this point through simulation studies in Section 3.1.

In the following sections, we propose three algorithms to identify the network Ω in (3) based on different dampening rate τ_i models. An algorithm for the homogeneous SGGM is first considered, and two extensions are proposed that allow the algorithms to deal with cohorts of heterogeneous subjects. These algorithms integrate the graphical LASSO algorithm with other algorithms, e.g., the EM algorithm, to find the penalized maximum likelihood estimator of Ω.

In this section, we consider identifying Ω under the assumption τ 1 = ⋯ = τ m = τ .

Thus, it is assumed that correlations between observations at different time points dampen at the same rate for each subject in the cohort. From the density function (2), the log-likelihood function for y = ( y 1 T , ⋯ y m T ) T is given by

up to a constant. Here, ⊗ stands for the Kronecker product. We use the formulas Σ i =   Φ i ⊗ Σ and | Σ i | = | Φ i | p | Σ | n i . Note that with formula ( Φ i ⊗ Σ ) − 1 = Φ i − 1 ⊗ Σ − 1 , the last term in (Equation 4) can be rewritten as

where S i ( τ ) = ∑ j = 1 n i ∑ k = 1 n i ϕ i j k − y i t k y i t i j T and Φ i − 1 = ( ϕ i j k − ) n i × n i . By substituting (Equation 5) into (4), we have

where n = ∑ i = 1 m n i , S ¯ ( τ ) = 1 n ∑ i = 1 m S i ( τ ) . Here, we use S ¯ ( τ ) to emphasize that matrix S ¯ is a function of unknown parameter τ. With (Equation 6) in hand, the sparse network can be achieved by minimization the L ₁-penalized negative log-likelihood function (3), i.e.,

for given tuning parameter λ > 0. The minimization problem (Equation 7) can be solved through a block coordinate descent procedure. First note that for a given τ, the solution of Ω can be obtained through the following minimization:

which has the same form as the GGM for independent data when the empirical covariance matrix is given by S ¯ ( τ ) . Consequently, the graphical LASSO algorithm can be used to compute the sparse estimate of Ω in (Equation 8). On the other hand, given Ω, the minimization of (7) with respect to τ does not involve any L ₁ penalty term and consequently can be carried out through the maximization of the likelihood function (6) with respect to τ. The conventional Newton algorithm can be used in this step. This process continues until convergence is achieved. This algorithm will be referred to as homogeneous longitudinal graphical LASSO (LGLASSO), for which the details are summarized in the following table.

In the homogeneous SGGM, we assume that a single correlation parameter τ applies to all the subjects. In real data analysis, this parameter may vary across subjects. In this section, we consider network identification without assuming τ 1 = ⋯ = τ m . Instead, we assume that the parameters τ_i ’s are independent random variables from the exponential distribution τ i ∼ exp ( α ) . Consequently, the joint density function for { y i , τ i } i = 1 m  is

from which the likelihood function for Σ and α is given by

With (Equations 9, 10) in hand, the sparse estimate for the network can then be obtained by minimizing the following L ₁-penalized negative log-likelihood function:

where l n ( Ω , α ) = log ( L n ( Ω , α | y ) ) . Since no explicit form for l n ( Ω , α ) is available, the expectation-maximization (EM) algorithm is proposed here to find the solution to (11) (Dempster et al., 1977). Since we are considering the negative log-likelihood function in (11), the maximization in the EM algorithm will be replaced by the minimization. The correlation parameters τ = ( τ 1 , ⋯ , τ m ) will be taken as the so-called missing data. Recall in the first step of the EM algorithm that the conditional distribution of missing data τ given y, Σ = Σ₀ , α = α ₀ has to be derived from (2) and (9) as follows:

For (12), the expectation of the complete log-likelihood function for (y, τ) to (12) has to be computed. Given the joint density function (9) of (y, τ), the expectation of its logarithmic transformation can be shown to be

where S ¯ ( τ ) = 1 n ∑ i = 1 m S i ( τ i ) , S ( 0 ) = E g S ¯ ( τ ) = 1 n ∑ i = 1 m E g S i ( τ i ) in which S i ( · ) is defined in the previous section. In the second step of the EM algorithm, the minimum point of the Q function in (13) has to be computed. This, again, is implemented through a block coordinate descent algorithm. First, for fixed Ω, it is straightforward to show that the minimum of the Q function with respect to α is attained at

i.e., the reciprocal of the sample mean of the conditional expectation of τ_i with respect to density (12). Then, for a given α in (Equation 14), the minimization of (13) with respect to Ω is equivalent to

which can be solved through the graphical LASSO algorithm. The difficult part of this algorithm is to find the expectation E g S ¯ ( τ ) , which may not have an explicit form given that S ¯ ( τ ) is a nonlinear function of τ and the complex form of density function g(τ); therefore, this expectation will be computed through a Monte Carlo method. This algorithm will be referred to as the heterogeneous longitudinal graphical LASSO, for which the details are summarized in the following table. Given initial values α ₀ and Ω₀, tuning parameter λ and error tolerance e > 0, the following algorithm is applied.

In the previous section, we assumed that E τ I = 1 / α is constant across the subjects. In this section, we further relax this constraint, and α can be a function of the covariates. Specifically, we assume τ i ∼ e x p ( α i ) , where α_i has the following form:

Here, x i =   ( 1 , x i 1 , ⋯ , x i q ) T represents covariates such as sex and race, and α = ( α 0 , ⋯ , α q ) T represents unknown parameters. The model (Equations 9–15) and Algorithm 2 in the previous section can then be revised straightforwardly to accommodate the current regression model (Equation 16). Specifically, first replace the parameter α in (Equations 9, 10) by e x p ( α T x i ) . Then, the conditional distribution of missing data is given by

Based on (Equation 17), the L ₁-penalized likelihood estimation of the MIN (11) is then given by

For the same reason, we need to use the EM algorithm to solve the optimization problem (Equation 18). The conditional distribution of missing value { τ i , i = 1 , ⋯ , m } is given by (12), with α ₀ replaced by exp ( α 0 T x i ) , from which the Q function is given by

Given the penalized Q function (Equation 19), the current estimate of the network Ω and parameter α, defined as the minimum point of Q, can be computed by the block coordinate descent algorithm. However, unlike (14), in the current case, the estimate of α does not admit an explicit form. We have to use numerical methods such as Newton algorithms to find the minimizer of the Q function (19). Specifically, we solve the following problem by using the BFGS algorithms:

Using α ^ in (Equation 20), the estimate of the MIN Ω is given by the solution of Equation (15) through the graphical LASSO algorithm. The algorithm is summarized in the following.

Remark: (1) All three algorithms described above leverage the graphical LASSO algorithm to achieve efficiency even though graphical LASSO itself is devised for independent data. (2) Note that correlation τ_i for subject i is a random variable. The forecast of τ_i , τ ^ í is given by the expectation of distribution (12), which is one of the outputs in Algorithms 2 and 3 . (3) The proposed algorithms can generate a solution path for a given sequence of tuning parameters λ. To select the optimal network Ω from the candidate networks, model selection criteria can be used, e.g., Akaike information criterion (AIC), Bayesian information criterion (BIC), or cross-validation (CV). In the numerical studies in the next section, we use the extended BIC (EBIC) that is dedicated to the graphical model to select λ (Foygel and Drton, 2010).

In this section, we compare the proposed algorithms, which are referred to as longitudinal graphical LASSO (LGLASSO) algorithms, with other existing network selection methods. These methods include graphical LASSO (Friedman et al., 2008; Friedman et al., 2019), neighborhood (Meinshansen and Bühlmann, 2006), GGMselectC01 and GGMselect-LA algorithms (Giraud et al., 2012; Bouvier et al., 2022). The graphical LASSO and neighborhood algorithms have been used in the SPIEC-EASI pipeline in (Kurtz et al., 2015) to select MINs and both are based on the L ₁-penalty. On the other hand, the GGMselect algorithm family provides different ways to construct and select the candidate models, e.g., GGMselectC01 employs the estimation procedure in (Wille and Bühlmann, 2006) to construct the candidate models, while GGMselect-LA uses the Fisher random variable to define the criterion for network selection. We demonstrate that the proposed longitudinal graphical LASSO algorithms can outperform these existing algorithms for simulated high-dimensional longitudinal microbiome data.

Let TP,P,FP,N,TN be the numbers of true positive edges, real positive edges, false-positive edges, null edges, and true null edges, respectively. In the following, we use the true/false positive rate (TPR/FPR) to measure the performance of each algorithm. They are defined as TPR = TP/P, FPR = FP/N. Also, for the conventional indices sensitivity/specificity, we have sensitivity = TPR, specificity = TN/N

The simulations consist of three scenarios. In the first scenario, we consider cases where the data follow the SGGM in Sections 2.2 and 2.3. Recall that the SGGM assumes that all the taxa in the microbiome share a common dampening rate τ_i for subject i. We refer to such microbiomes as having a homogeneous community. In the second scenario, we consider the microbiomes that violate such homogeneity, i.e., having heterogeneous communities. In scenario three, we consider the left-censored microbiome data, which aims to test the robustness of the SGGM with respect to zero inflation. The zero-inflation phenomena are widely observed in 16S rRNA gene sequencing experiments and violate the assumptions of the SGGM. Both the homogeneous and heterogeneous versions of LGLASSO in Sections 2.2 and 2.3 are investigated in each scenario. We use the receiver operating characteristic (ROC) curve to show the superiority of our algorithms over the graphical LASSO and the neighborhood algorithms in all scenarios. In the case of the heterogeneous LGLASSO, we also use (TPR,FPR) to compare the performances of the algorithms in which the networks are selected through the extended BIC (EBIC). The methods of GGMselect-C01 and GGNselect-LA have their own model selection approaches for network selection. The EBIC for the Gaussian graphic model is given by

where l_n is the log-likelihood function, n is the sample size, G is the network of interest, | G | is the number of edges in G, p is the number of nodes, and T is the tuning parameter. In (Equation 21), we choose a typical value T = 2 for model selection.

In this section, a longitudinal data set from a cohort of children with cystic fibrosis was investigated using the homogeneous version of the proposed algorithm in Section 2.2. Specifically, stool samples from thirty-eight children were collected from children aged 6 months to 51 months old (Madan et al., 2012). The number of observations from each child ranged from 2 to 17. Each observation consisted of the abundance of 16,383 amplicon sequence variants (ASVs) of the 16S rRNA gene. These sequences were then collapsed to the genus level using the R package DADA2 (Callahan et al., 2016). The sequences that had no genus-level information were dropped. Then, all the taxa with a proportion of nonzero observations less than 10% were combined, which was referred to as the composite taxon. There were 83 total remaining taxa. The observations of zeros for each of these 83 taxa were replaced by the minimum abundance of that taxon divided by 10. The log-ratio transformation was then carried out to obtain the relative abundance, in which the composite taxon was used as the reference. Similar log-ratio transformations have been used and justified in empirical studies (Kurtz et al., 2015; Greenacre et al., 2021).

The application of the homogeneous model in Section 2.2 to the transformed data yielded the estimated network, which is displayed in Figure 7 . Based on the modularity maximization algorithm (Newman, 2006; Blondel et al., 2008), five communities were identified in the estimated network, which is listed in Table 2 .

To show that the estimated network can reveal the true structure of the underlying network, we investigated the relationship between the estimated network and the phylogenetic tree of the 82 taxa (the composite network was excluded here). The phylogenetic tree constructed from the same data set is presented in Figure 8 and demonstrates the evolutionary relationship among these taxa. Our underlying hypothesis is that microbial taxa that are closer in terms of evolutionary history also have more/stronger interactions in the human body. To validate this hypothesis, the null hypothesis is set as follows: the estimated network in Figure 7 is independent of the phylogenetic tree in Figure 8 . Correlation between the estimated network and phylogenetic tree is employed to test this hypothesis. In particular, we computed the distances between two taxa in the estimated network and the phylogenetic tree. Here, the distance between taxa A and B is defined as the length of the shortest path from A to B in the estimated network (phylogenetic tree). If no paths exist between two taxa in the network (phylogenetic tree), that pair will be excluded from the following computations. Let d ₁ and d ₂ be the distance vectors for all possible pairs of taxa from the estimated network and phylogenetic tree, respectively. The correlation between d ₁ and d ₂ is used to measure the relatedness between the network and phylogenetic tree. This correlation was determined to be r _{0 =} 0.333 for the tuning parameter selected by the EBIC.

To understand the significance of r0 against the null hypothesis, we use the permutation method to estimate the null distribution. Specifically, we keep the structure of the estimated network unchanged and permute the order of the 82 taxa m = 5000 times on the estimated network. L e t   d 1 ( i ) ( i = 1 , ⋯ , 5000 ) be the distance vectors of the network for the ith permutation. Then, the correlations r ( i ) = c o r ( d 1 ( i ) , d 2 ) , i = 1 , ⋯ , 5000 , which collectively depict the null distribution, can be derived. Given r ^{( i )}, the p value of r ₀ is smaller than 1/5000, which means that the correlation between the estimated network and phylogenetic tree is statistically significant, i.e., the data support the hypothesis that microbial taxa with closer evolution histories tend to have more/stronger interactions.

It should be noted that some related findings have been reported in the literature. In (Chaffron et al., 2010), the authors performed a global meta-analysis of previously sampled microbial lineages in the environment. They found that genomes from coexisting taxa tended to be more similar than expected by chance, both with respect to pathway content and genome size. The studies in (Eiler et al., 2012) also revealed that ecological coherence is often dependent on taxonomic relatedness. These studies employed coefficient-based methods such as Fisher’s exact test to infer the interaction of taxa. This can lead to a misleading conclusion. It is known that the correlation between two taxa, A and B, may be induced by their correlation with a third taxon C, even though A and B are independent if C is fixed. In the current study, the interaction between the taxa is defined based on the conditional correlation coefficient, which by its definition eliminates the possible spurious correlation between taxa A and B induced by taxon C. Therefore, by using abundance instead of cooccurrence information and boosted by the proposed methods, we reach a more convincing and robust conclusion than existing ones in the literature.

The results above are derived with the fixed tuning parameter selected by the EBIC. The correlations between the estimated MIN and phylogenetic tree are actually robust with respect to the tuning parameter. To demonstrate this point, let us consider the networks generated from other tuning parameter choices. Specifically, for each of the 40 tuning parameters ranging from λ = 7 to 13, the same homogeneous model and the permutation test are carried out. The estimated and permuted correlations are displayed in Figure 9 in the form of a boxplot. From left to right, the boxplots in Figure 9 correspond to the tuning parameters increasing from 7 to 13. The dots linked by the line represent the correlations between the estimated MINs and the phylogenetic tree, while others correspond to the correlations computed from the permuted MINs. The most prominent feature of Figure 9 is that for all 40 cases, the observed correlations between the MINs and the phylogenetic tree are positive; in the first half of the boxplots, the estimated correlations are also significant. It should be noted that even though the estimated correlations in the second half of the boxplots do not appear significant, it does not mean that the edges in these estimated networks do not reflect the true structure. Instead, the insignificance may stem from the fact that the second half of the networks are sparser. A sparse network will generate shorter vectors d ₁ , d ₂, which in turn increase the variability of the correlation estimates.

Finally, let us compare the results of the proposed algorithm with those of the SPIECEASI algorithm in (Kurtz et al., 2015). Though in the original form of the SPIEC-EASI algorithm, the centered log-ratio transformation was employed for the relative abundances of the taxa, we use the additive log-ratio transformation here. Note that the SPIEC-EASI algorithm includes both graphical LASSO and neighborhood algorithms. We restrict ourselves to networks with edges between 20 and 1300, which we believe cover all biologically meaningful cases. First, we compute the solution path for each algorithm. For the networks on each solution path, their correlations with the phylogenetic tree are calculated. Figure 10 displays the correlations for the three solution paths corresponding to the three algorithms. It is evident that for most of the solution paths, the networks generated through the proposed algorithm have higher correlations with the phylogenetic tree than those generated through the SPIEC-EASI algorithms. The correlations at the beginning parts of the three paths appear to be comparable. We attribute this to the fact that the networks at the beginning parts are much sparser, leading to a smaller sample size when computing the correlations. A smaller sample size can blur the comparison of different algorithms, as shown in Figure 9 . In other words, if we assume that the phylogenetic tree represents the true structure of the MIN, then the proposed algorithms have greater power in the identification of the MIN than that of the SPIEC-EASI algorithm.