In this section, we formally define a time-lagged effect and its corresponding time lag. For each of the n subjects, we collect repeated measures on p covariates, X 1 , … , X p , each of which is measured T times, denoted by X j = ( X j , 1 , … , X j , T ) ′ , j = 1 , … , p . In a high-dimensional setting, we have p ≫ n . In addition, we measure a response Y , recorded at the final time point T for each subject, thus denoted by Y T . The relationship between the repeated measures of the covariates X 1 , … , X p and the response Y T is modeled through a linear model or a generalized linear model in which g E Y T | X 1 , … , X p = β 0 + X 1 ′ β 1 + ⋯ + X p ′ β p , (1) where β j = ( β j , 1 , … , β j , T ) ′ are the coefficients for X j = ( X j , 1 , … , X j , T ) ′ and g ( ⋅ ) is a link function depending on the distribution of the response. For example, g ( ⋅ ) can be the identity function when Y T follows a normal distribution and the logit function when Y T follows a Bernoulli distribution.

In this work, we focus on identifying and estimating time-lagged associations between the covariates and the response. We start with the definition of the time-lagged effect of a covariate on the response. To facilitate the definition, we illustrate the underlying relationship between two covariates and the response over time in Figure 1. In Figure 1, we lay out the repeated measures of two covariates X 1 and X 2 over time as well as the measurement of the response Y at the final time point T . We also include the hypothetical measurements of the response Y at previous time points t = 1 , … , T − 1 to facilitate the interpretation of the time-lagged effect, although these measurements were not collected in practice (indicated by the dotted circles instead of solid circles in the diagram). The longitudinal dependence between the repeated measures of a covariate or between the repeated measures of the response is shown as dashed arrows, indicating the causal effect from a previous time point to the next. For simplicity of illustration, we did not include confounders for the repeated measures of covariates or the response in Figure 1, although adding them does not really change the interpretation of time-lagged effects.

In Figure 1, on the one hand, we say that X 1 has an instantaneous effect (with no time lag or lag 0) on Y , as shown by the solid directional arrows from X 1 , t to Y t at every time point t = 1 , … , T , as if Y 1 , … , Y T − 1 were measured. In this case, the instantaneous effect of X 1 , T on Y T , represented by β 1 , T in model 1, is nonzero. However, given that we do not observe the responses at the previous time points, Y 1 , … , Y T − 1 , we need to assess β 1,1 , … , β 1 , T − 1 in model 1, the effects of X 1,1 , … , X 1 , T − 1 on Y T . From the figure, it is seen that there are indirect associations between X 1,1 , … , X 1 , T − 1 and Y T through the longitudinal dependence between the repeated measures of either X 1 or Y . Even if we include all repeated measures of X 1 in model 1, the effects of X 1,1 , … , X 1 , T − 1 on Y T , represented by β 1,1 , … , β 1 , T − 1 , remain nonzero in model 1 due to the longitudinal dependence between the hypothetical repeated measures of the response. In summary, if X 1 has an instantaneous effect on Y , the coefficients in model 1) possess the following pattern: β 1,1 ≠ 0 , … , β 1 , T − 1 ≠ 0 , β 1 , T ≠ 0 .

On the other hand, X 2 has a 1-lagged effect on Y (time-lagged effect with lag 1), i.e., X 2 , t has a direct effect on Y t + 1 for any given time t = 1 , … , T − 1 , indicated by the solid directional arrows from X 2 , t to Y t + 1 for t = 1 , … , T − 1 . In this case, the 1-lagged effect of X 2 , T − 1 on Y T , represented by β 2 , T − 1 , is obviously nonzero. We still need to assess β 2,1 , … , β 2 , T − 2 and β 2 , T in model 1. Similar to X 1 , β 2,1 , … , β 2 , T − 2 are also nonzero in the model due to the indirect associations between X 2,1 , … , X 2 , T − 2 and Y T caused by the longitudinal dependence between the hypothetical repeated measures of the response. However, β 2 , T in 1 is zero, because the indirect association between X 2 , T and Y T is through X 2 , T − 1 , which has been included in the model. In summary, if X 2 has a 1-lagged effect on Y , the coefficients in model 1 possess the following pattern: β 2,1 ≠ 0 , … , β 2 , T − 1 ≠ 0 , β 2 , T = 0 .

Based on the above discussion, we formally define the time-lagged effect of a covariate X j on Y T with lag d through the following sparsity pattern of the corresponding coefficients β j = ( β j , 1 , … , β j , T ) ′ : β j , 1 ≠ 0 , … , β j , T − d ≠ 0 , β j , T − d + 1 = 0 , … , β j , T = 0 . (2)

In other words, if X j has a d -lagged effect on Y T , the first T − d coefficients are nonzero and the last d coefficients are zero, as it takes d time intervals for X j to impact Y T .

The time-lagged effect of microbiome on host status is not uncommon in practice. For example, our real data analysis identifies several microbial taxa in the gut that have a variety of lagged associations with zebrafish parasite worm burden. In particular, abundances of genus Chitinibacter were found to have a 29-day-lagged association with parasite worm burden, whereas abundances of genus Mycobacterium were found to have an instantaneous association with parasite worm burden. More details of such findings can be found in the real data section.

Due to the correspondence between the time-lagged effect and the sparsity pattern of the grouped coefficients in (2), identifying and estimating time-lagged effects can be regarded as a grouped variable selection problem. To identify these effects, we aim to pinpoint which time points have non-zero coefficient values, indicating a true association with the response. This process parallels traditional variable selection, where the objective is to determine which covariates contribute to the model. By considering the measurements of a covariate at different time points as a group, we treat the identification of time-lagged effects as a grouped variable selection problem, often solved using group penalization approaches, ensuring that only the relevant time points are retained in the model.

In the group penalization framework, we estimate the parameters β = ( β 1 ′ , … , β p ′ ) ′ in model 1 by minimizing a objective function of the following form: Q β | X , Y = L β | X , Y + P β , where L ( β | X , Y ) is the negative log-likelihood function depending on the underlying model 1 and P ( β ) is a penalty function, with observations X = ( ( X 1,1 ′ , … , X 1 , p ′ ) ′ , … , ( X n , 1 ′ , … , X n , p ′ ) ′ ) ′ and Y = ( Y 1 , T , … , Y n , T ) ′ of sample size n .

In this subsection, we briefly review existing grouped variable selection methods as well as their corresponding penalty functions P ( β ) by classifying them into two categories—group-level/bi-level selection methods and overlapping group selection methods.

To mimic the real data in our simulation, we make use of the microbiome data from a longitudinal study sampling the zebrafish microbiome (Hammer et al., 2024). The fecal microbiome of zebrafish was analyzed across three separate days ( T = 3 ) for 38 taxa ( p = 38 ) . After filtering for samples present in each of the 3 days, we are left with an initial sample size of 21 ( n = 21 ) . Measurements are then transformed using the centered log-ratio transformation (Aitchison, 1982) and serve as the covariates X 1 , … , X p . Details of this study and the data can be found in the real data section.

We simulate the response with two distributional settings, the normal distribution and the Poisson distribution. With the normal distribution, we simulate the response using a linear model: Y = X 1 ′ β 1 + ⋯ + X p ′ β p + ϵ where ϵ ∼ N ( 0,1 ) . With the Poisson distribution, we set μ = exp ( X 1 ′ β 1 + ⋯ + X p ′ β p ) , and simulate the response using a Poisson model: Y ∼ Poisson ( μ ) .

For longitudinal data sampled at three time points, there are four possible lags, lag 0, 1, 2, and 3, where lag 3 corresponds to no association. Of the 38 taxa, we set a sparse signal of six true associations, two instances of each of lag 0, 1, and 2. We randomly assign the six true associations among the 38 taxa. Additionally, we test three magnitudes of the signal of the true associations, small, medium, and large. For the linear case these are β j , t ≡ c when β j , t ≠ 0 , where c ∈ { 0.1 , 0.5 , 2 } . For the Poisson case these are β j , t ≡ c when β j , t ≠ 0 , where c ∈ { 0.01 , 0.05 , 0.1 } . These values were chosen to simulate counts in a similar range to those present in the real data.

The original sample size of the zebrafish study ( n = 21 ) represents a small, but not unusual, sample size. We also want to test other sample sizes. To do this, we resample from the original 21 samples, and add noise from N ( 0,1 ) , resulting in a medium ( n = 50 ) and a large ( n = 100 ) sample size.

We compare seven grouped variable selection methods in two categories: group-level/bi-level selection methods and overlapping group selection methods. For the former, we applied group lasso (GrLasso), group exponential lasso (GEL), composite MCP (cMCP), and group bridge (GrBridge) to the simulated data. For the latter, we applied overlapping group lasso (O-GrLasso), overlapping group MCP (O-GrMCP), and overlapping group SCAD (O-GrSCAD).

We compare the performance of the seven grouped variable selection methods with various simulation settings: sample size ( 21,50,100 ) and signal magnitude (small, medium, large), for both linear model and Poisson model. The simulation results are summarized from 100 replicates. In the Poisson model, the GrBridge method did not converge in half of the simulation replicates in the medium-signal case when n = 100 , and never converged in the large-signal case. Results are shown only for the cases in which the algorithm did converge.

We apply our group penalization framework to the real dataset from Hammer et al. (2024), which originally studied the role of the gut microbiome in mediating parasitic infection of Pseudocapillaria tomentosa in zebrafish. In our application, we make use of the data collected from longitudinally sampled zebrafish fecal samples (sampled on days 0, 3, and 32) to find time-lagged associations between the abundances of microbial taxa in the zebrafish gut and the parasite burden on zebrafish.

After day 0 and before day 3, half of the tanks of the zebrafish in the study were given an antibiotic and the other half were not. After day 3, in each group of zebrafish (antibiotic and control), roughly half of them were exposed to the parasite P. tomentosa. All zebrafish were sacrificed to assess intestinal histopathology on day 32, and the parasite burden on zebrafish was measured. In other words, microbiome data from the zebrafish gut were collected (a) prior to antibiotic exposure (day 0), (b) just prior to parasite exposure but after antibiotic exposure (day 3), and (c) 29 days post-parasite exposure (day 32). Figure 6 further explains the experimental design.

Since we use the final parasite burden as our response Y (day 32), our analysis only focuses on the half of fish that were exposed to the parasite. After filtering for die-off, 21 zebrafish remain in the parasite exposed group, 9 of which were given an antibiotic and 12 were not. In addition, we examine taxa at the genus level and apply an inclusion threshold to analyze genera present in at least 30% of samples, resulting in p = 38 genera being included in our analysis.

In our analysis, the response Y is the final parasite burden on zebrafish, measured by the counts of the parasite P. tomentosa via intestinal histopathology. Therefore, we use a Poisson regression model for Y , with the link function in (1) chosen as the logarithm function. In addition, we apply a centered log-ratio transformation on the genus abundances collected on each day, and treat them as the repeated measures of the genera X 1 , … , X p where X j = ( X j , 0 , X j , 3 , X j , 32 ) ′ for j = 1 , … , p .

We use these data to find time-lagged associations between final parasite burden of the zebrafish and genus abundances measured on day 0, 3, and 32. Since our samples are split into two groups, those exposed to antibiotics and those not exposed, we also take potential effects of antibiotic exposure into account in our modeling. To this end, we include antibiotic exposure as a main effect ( A , an indicator function for antibiotic exposure) as well as its interaction effects with the longitudinal genus abundances ( A X 1 , … , A X p ) , where A X j = ( A X j , 3 , A X j , 32 ) ′ for j = 1 , … , p . It is noteworthy that there is no interaction between A and X j , 0 as the zebrafish are exposed to the antibiotic after day 0. In summary, the model in this data analysis can be written as: Y | A , X 1 , … , X p ∼ Poisson μ , where  log μ = α 0 + A α 1 + X 1 ′ β 1 + ⋯ + X p ′ β p + A X 1 ′ γ 1 + ⋯ + A X p ′ γ p , (3) where β j = ( β j , 0 , β j , 3 , β j , 32 ) ′ and γ j = ( γ j , 3 , γ j , 32 ) ′ for j = 1 , … , p .

We do not penalize the main effects α 0 and α 1 , and penalize the main effects β 1 , … , β p and the interaction effects γ 1 , … , γ p using either group-level/bi-level selection methods or overlapping group selection methods. For group-level/bi-level selection methods, a group of coefficients is formed for each taxon, ( β j , 0 , β j , 3 , β j , 32 , γ j , 3 , γ j , 32 ) ′ for j = 1 , … , p . Group-level/bi-level selection methods then select from these groups and bi-level selection methods further select individual parameters in each group. For overlapping group selection methods, we construct 6 groups for each taxon j : ∅ , β j , 0 , β j , 0 , β j , 3 ′ , β j , 0 , β j , 3 , β j , 32 ′ , β j , 0 , β j , 3 , γ j , 3 ′ , β j , 0 , β j , 3 , β j , 32 , γ j , 3 , γ j , 32 ′ , Corresponding to different patterns of association with different lags. The first empty set corresponds to no association, the next three correspond to main effects with lags 32, 29, and 0 days, the last two correspond to main and interaction effects with lags 29, and 0 days. Note that the above construction of groups enforces the main effect to be included in the model before the corresponding interaction effects, whereas the bi-level selection methods have no such restriction.

We tabulate in Table 2 the list of genera that are identified to be associated with parasite burden, together with their lags of such associations. For main effects, there are three possible lags—0, 29, and 32 days; for interaction effects, there are two possible lags—0 and 29 days, as there is no interaction between the antibiotic treatment and the microbial abundance on day 0. A lag of 0 days indicates an instantaneous effect, a lag of 29 or 32 days implies there is a lag for the microbial abundance to affect the parasite burden, either after antibiotic exposure or before antibiotic exposure. From the results in Table 2, we can draw the following conclusions.

First, the results highlight the importance of identifying time-lagged associations. Five of the eleven identified genera have a time-lagged effect, which would not have been found if we were not using the proposed approaches.

Second, the overlapping group selection methods perform more similarly to each other, with O-GrLasso and O-GrSCAD producing the same results, while the bi-level selection methods do not have many similarities amongst themselves. GrBridge did not even identify any associated taxa. The most commonly identified genus was Candidatus Accumulibacter which was present in every method except O-GrMCP and GrBridge.

Third, we see a variety of identified lags. Candidatus Accumulibacter was identified with an instantaneous effect (lag of 0 days) by O-GrLasso/O-GrSCAD and GEL, Chitinibacter was found to have a lag of 29 days, including a 29-day lagged antibiotic interaction effect by all overlapping methods. Hyphomicrobium had a lag of 29 days by O-GrLasso and O-GrSCAD, and a lag of 32 days by O-GrMCP and cMCP.

Table 3 further presents the estimated coefficients for the identified genera in each model. Table 3 demonstrates a key difference between the overlapping group selection methods and the bi-level selection methods, namely, whether or not the sparsity pattern in (2) is met. For example, Candidatus Accumulibacter has a lag of 0 days for O-GrLasso and O-GrSCAD, as well as GEL, but GEL does not find an association on day 3, violating the sparsity pattern in (2). Additionally, Candidatus Accumulibacter is identified by cMCP to have an association with worm burden, but only the interaction effect on day 32, with no effects on earlier days and no main effects.

The original paper (Hammer et al., 2024) including these zebrafish data conducted an analysis on the mediating role played by the gut microbiome on the relationship between gastrointestinal metabolites and parasitic infection outcome. This section highlights how our findings can further inform the longitudinal aspect of the relationship between the gut microbiome and parasite burden. Many of the genera the mediation study comments on as interesting are also found by our approach, and we expand upon these genera below.

Hammer et al. (2024) identified taxa in the Pseudomonas and Mycobacterium genera to be mediators in the relationship of the important Vitamin E metabolite γ -tocopherol and final worm burden. These genera are also present in our results. We found that abundance of Pseudomonas has a time-lagged and negative association with worm burden, where notably, the association only exists at the first time point (pre-antibiotic and parasite exposure). We also find that Mycobacterium has an instantaneous and negative association with worm burden and see that it gets stronger as the time is closer to the final time point.

Our findings also unite previous research which used these data by uncovering temporal changes in microbial association that help explain microbiota connections to parasite infection burden. For instance, results from Hammer et al. (2024) show that salicylaldehyde, which may be partially controlled by Pelomonas, has particularly strong effects on egg larvation and development. Our finding here that there exists an association between Pelomonas and infection burden on day 3 is both novel and important because this is the time fish hosts were exposed to parasite eggs and based on experimental evidence from stemming from this earlier work it is expected that salicylaldehyde-related inhibition of helminth maturation would be most pronounced at this time point in the study. Thus, this association from O-GrMCP points to time-dependent activity of Pelomonas that could help to explain and unite a connection between Pelomonas and salicylaldehyde to in vivo and in vitro results that implicate salicylaldehyde as an anthelmintic agent.

Additionally, results of applying these methods highlight a possible route by which gut microbes might regulate host intestinal structure to limit helminth infection, by regulating tight junction integrity. Results from both O-GrMCP and GEL point to day-3 associations in Cetobacterium are inversely related to helminth parasite infection. Prior work using the zebrafish model has shown that taxa within Cetobacterium synthesize vitamin B12 which mechanistically enhances gut barrier tight junction integrity to prevent microbial pathogen infection and improve gut microbiome structure stability (Qi et al., 2023). These findings are also relevant to nematode infection, where it has been shown using other infection models that early parasite exposure results in loss of epithelial barrier integrity as a result of changes in tight-junction related protein expression (Fernández-Blanco et al., 2015). Similar tight-junction regulating activity could underscore our results which point to early time points of Cetobacterium relative abundance negatively associating with parasite infection burden, potentially as a result of vitamin B12 biosynthesis.

Overall, we identify similar taxa associating with parasite worm burden as in Hammer et al. (2024), but our results provide important nuance to these findings by revealing time-dependent microbial associations with infection burden. For instance, our finding that Pelomonas is inversely associated with infection parasite burden on day 3 could help explain distinct results from Hammer et al. (2024) that revealed a connection between Pelomonas with salicylaldehyde, and salicylaldehyde to egg larvation. Furthermore, the connection between Cetobacterium and infection burden also uncovers a testable hypothesis regarding the relationship between microbes to metabolites and parasite infection that was not previously identified in the Hammer et al. (2024) analysis, and points to an additional route by which microbes could regulate helminth parasite infection. Together, clarifying the time during which a microbe might produce potent anthelmintic products or influence host response to infection can elucidate insights into the activity of microbes across the time range of parasite infection, possibly imparting new ways the gut microbiome may be harnessed to combat helminth parasite infection.

To validate the microbial taxa that were identified to be associated with parasitic infection, we conduct a validation analysis using the additional measurements of metabolites in the zebrafish study (see Figure 6). We focus our analysis on the two metabolites that were found to be linked to parasite infection and whose effect on infection burden is mediated by members of the gut microbiome (Hammer et al., 2024), namely, salicylaldehyde and γ -tocopherol. We refer to Section 4.4 for detailed discussion on the role of these two metabolites in the microbiome-parasite ecosystem. Instead of worm burden used as a response, we now use each of the two metabolite measurements as the response and investigate the (potentially time-lagged) effects of microbial taxa on these metabolites. This additional analysis serves as the validation of the identified microbiome-parasite associations.

In the validation analysis, the response Y is the final metabolite measurement on day 32. Therefore, we use a Gaussian regression model for Y , with the link function in (1) chosen as the identity function. Metabolite measurements were log transformed. Since parasite burden is not included in this analysis, we include all groups of fish used in the initial study, i.e., using both the parasite and control groups. After filtering for die-off, we have 59 overall fish.

We use these data to find time-lagged associations between genus abundances measured on day 0, 3, 32 and final metabolite levels of the zebrafish measured on day 32. Since our samples are split into four groups, those exposed and not exposed to antibiotics, as well as those exposed and not exposed to parasites, we take potential effects of antibiotic and parasite exposure into account in our modeling. Compared to the model in (3), we include parasite exposure as an additional main effects ( P as an indicator function for parasite exposure). Thus, the model for the validation analysis can be written as: Y | A , P , X 1 , … , X p ∼ N μ , σ 2 , where  μ = α 0 + A α 1 + P α 2 + X 1 ′ β 1 + ⋯ + X p ′ β p + A X 1 ′ γ 1 + ⋯ + A X p ′ γ p , (4) where β j = ( β j , 0 , β j , 3 , β j , 32 ) ′ and γ j = ( γ j , 3 , γ j , 32 ) ′ for j = 1 , … , p .

Similar to Section 4.2, we use bi-level/group-level selection methods and overlapping group selection methods to estimate the coefficients in (4). We do not penalize the main effects α 0 , α 1 , and α 2 , and penalize the main effects β 1 , … , β p and the interaction effects γ 1 , … , γ p . The same groups of coefficients are formed as in Section 4.2 for group-level/bi-level selection methods and overlapping group selecion methods. The results of this validation analysis are presented in Tables 4–7. Tables 4, 5 show the time lags of microbial effects on salicylaldehyde and γ -tocopherol, respectively; Tables 6, 7 report the coefficient estimates of the above associations.