Critical Window Variable Selection for Mixtures: Estimating the Impact of Multiple Air Pollutants on Stillbirth

Throughout their lives, humans are simultaneously exposed to multiple contaminants that may adversely impact their health. Ambient air pollutants represent a major source of potentially hazardous exposures. Numerous past studies have estimated associations between ambient exposures and multiple adverse health outcomes, with some findings consistent enough to suggest a causal relationship (Brunekreef and Holgate, 2002; Kampa and Castanas, 2008; Stieb et al., 2012). However, the majority of these studies relied on statistical models that examined the role of a single pollutant on a health outcome. Therefore, the findings may represent an incomplete understanding of the impact that pollution mixtures have on health. New statistical methods are needed to address this issue, though a number of challenges in their development have been previously identified, particularly high correlation between exposures due to shared emission sources and meteorological drivers (Dominici et al., 2010).

Recently, there has been a shift away from single pollutant approaches towards the development of multipollutant methods (e.g., Davalos et al. (2017); Ferrari and Dunson (2019, 2020); Antonelli et al. (2020); Reich et al. (2020)). However, the majority of these methods were not designed for the analysis of time-varying mixtures, a feature that introduces additional modeling complications such as strong temporal correlation between exposures. Methods that examine the role of time-varying pollutants on health are particularly important for identifying critical periods of exposure as is often of interest when analyzing pregnancy outcomes, though other health outcomes have also been explored in this context (Warren et al., 2012a, b, 2016, 2020a, 2020b; Chang et al., 2015; Wilson et al., 2017).

Methods that seek to identify critical windows of exposure most often rely on distributed lag models (DLMs) in some form to quantify the relationship between exposure over time and health, with the majority of methods developed for a single or small number of pollutants. Using a Gaussian process, Warren et al. (2012b) introduced a model which allowed for two pollutants to be analyzed jointly while accounting for their cross-covariance. Warren et al. (2012a) extended this model using a semiparametric Bayesian approach which more generally considered four pollutants. However, both frameworks became computationally demanding as the number of pollutants increased. Chen, Mukherjee and Berrocal (2019) proposed a method to explore the impact of two time-varying exposures, including interactions, within a DLM framework. Liu et al. (2018a) extended Bayesian kernel machine regression (BKMR) to accommodate time-varying mixtures and further introduced a computationally feasible framework using mean field variational Bayesian inference (Liu et al., 2018b). However, the authors noted that this method is better suited for a smaller number of shared exposure time periods. More recently, Wilson et al. (2019) introduced BKMR-DLM, a method designed to more flexibly and generally handle multiple exposure time periods while exploring nonlinear and interaction associations between pollutants and the outcome.

Weighted quantile sum (WQS) regression (Carrico et al., 2015) was developed to investigate the impact of multiple non-time-varying exposures on an outcome. It introduces regularization via a set of pollutant-specific weights that are forced to sum to one to stabilize model fitting when considering a large number of contaminants. Extending WQS to the setting of multiple exposure periods, Bello et al. (2017) proposed an algorithm that relied on a DLM, though Gennings et al. (2020) later suggested that it may be more appropriate for a smaller number of pollutants. As a result, Gennings et al. (2020) introduced lagged WQS (LWQS) regression which proposed to fit WQS regression separately at each exposure period to estimate the time-varying weight parameters and obtain the weighted average exposures across time. It then used a DLM to estimate the association between the estimated weighted average exposures and the outcome. However, this two-stage algorithm may be ignoring sources of uncertainty during model fitting (e.g., weight parameters are assumed known in the second stage DLM), may be inefficient since weights are estimated separately at each exposure period (i.e., ignoring potential correlation between weights in nearby exposure periods), and ignores the impact of interactions between exposures.

In this work, we extend the recently developed method for more accurately identifying critical exposure windows in a single pollutant setting (Critical Window Variable Selection (CWVS)) (Warren et al., 2020a) to the multipollutant setting by introducing CWVS for Mixtures (CWVSmix). Similar to LWQS, CWVSmix introduces sum-to-one constraints on time-varying mixture weights in order to introduce regularization and accommodate a large number of pollutants. However, unlike LWQS, CWVSmix is designed to jointly estimate these weights and the estimated mixtures’ impacts on health while employing hierarchical Bayesian variable selection techniques to more formally identify critical windows of exposure and the important pollutants involved in the time-varying mixtures. Additionally, first order interactions between pollutants on a given exposure period are also considered within the framework.

We design a simulation study to compare the performances of CWVSmix and other competing options (including LWQS) with respect to (i) accuracy of the identified critical exposure periods, (ii) estimation of the time-varying pollutant-specific weights, and (iii) risk parameter estimation. We then apply the methodology to investigate the impact that multiple gas and particulate matter ambient air pollutants, estimated using a novel data fusion model (Senthilkumar et al., 2019), have on stillbirth risk in New Jersey (NJ), 2005-2014, with a focus on investigating differences between various race/ethnicity groups.

We analyze live birth and fetal death records from the Division of Family Health Services in the NJ Department of Health in 2005-2014. All births/deaths included in our study are singletons with a clinically estimated gestational age of at least 20 weeks, no birth defects, and a conception date in 2005-2013 (based on gestational age and date of birth/death) to avoid the fixed-cohort bias (Strand, Barnett and Tong, 2011). A stillbirth was defined as a fetal death occurring on or after 20 weeks of gestation (e.g., Faiz et al. (2012); Green et al. (2015)). From these records, we extracted information on maternal educational attainment ($<$ high school, high school, $>$ high school), maternal age category ($<25$, $[25,30)$, $[30,35)$, $>35$), maternal race/ethnicity (non-Hispanic Black, Hispanic, non-Hispanic White), maternal tobacco use during pregnancy (yes, no), maternal residence at delivery (latitude, longitude), sex of the fetus (male, female), and the season and year of conception.

Due to the large number of live births in each year of analysis and the low prevalence of stillbirth, we opt for a case-control study design similar to Warren et al. (2020a) who investigated the impact of air pollution on very preterm birth. Additionally, we analyze each maternal race/ethnicity group separately to better understand the unique associations corresponding to each group instead of assuming common patterns of risk. Therefore, for each eligible stillbirth in the dataset, we select five controls (i.e., live births) while matching only on race/ethnicity. In Tables S1-S3 of the Supplement, we display summary information for each of the three race/ethnicity groups by stillbirth status.

We obtained daily estimates of ambient air pollution concentrations from 12 pollutants that covered a 12 by 12 kilometer grid across NJ from 2005-2014 with no spatial or temporal missingness. The pollutants include 24-hour average particulate matter with an aerodynamic diameter less than or equal to 2.5 and 10 microns (PM_2.5 and PM₁₀, respectively) and PM_2.5 constituents nitrate (NO${}_{3}^{-}$), ammonium (NH${}_{4}^{+}$), sulfate (SO${}_{4}^{2-}$), elemental carbon (EC), and organic carbon (OC); 1-hour maximum carbon monoxide (CO), nitrogen dioxide (NO₂), nitrogen oxides (NO${}_{\text{x}}$), and sulfur dioxide (SO₂); and 8-hour maximum ozone (O₃). These estimates come from a recently introduced spatiotemporal data fusion model that combines measured ambient concentrations with deterministic output from the Community Multiscale Air Quality model (Senthilkumar et al., 2019). The latitude/longitude of the residence at delivery for each birth in the study was linked with the closest 12 by 12 kilometer grid cell centroid and weekly exposures for each pollutant (through gestational week 20) were calculated based on the date of birth and gestational age. This study was approved by the institutional review boards at Rowan and Yale Universities.

We introduce CWVSmix, a model for (i) identifying critical windows of exposure to mixtures of time-varying pollutants, (ii) estimating the time-varying relative importance of each individual pollutant and their first order interactions within the mixture, and (iii) quantifying the impact that the identified mixtures have on a health outcome. CWVSmix incorporates hierarchical Bayesian variable selection techniques for defining critical windows and identifying important pollutants within a time-varying mixture, and uses Gaussian processes to allow for smoothness/regularization when estimating a potentially high-dimensional set of parameters that vary across time.

Specifically, we extend CWVS (Warren et al., 2020a) from the single pollutant setting to accommodate multiple exposures and their interactions such that

where $Y_{i}$ is the binary adverse health outcome of interest observed for study participant $i$; $n$ is the total number of people in the study; $p_{i}$ is the probability of experiencing the adverse health outcome; $\textbf{x}_{i}$ is a vector of person-specific covariates/confounders with accompanying regression parameters $\boldsymbol{\beta}$; $\text{z}_{ij}\left(t\right)$ represents exposure to the $j^{\text{th}}$ pollutant ($j=1,\ldots,q$), corresponding to person $i$’s spatial location, averaged over exposure period $t$ ($t=1,\ldots,m$) corresponding to the relevant calendar date range specific to person $i$’s timing of exposure (e.g., calendar dates that align with the selected gestational week); $\lambda_{j}\left(t\right)\in\left[0,1\right]$ are time-varying weight parameters that describe the main effect of pollutant $j$ in the mixture corresponding to exposure period $t$; $\widetilde{\lambda}_{jk}\left(t\right)\in\left[0,1\right)$ describes the interaction effect between pollutants $j$ and $k$ on exposure period $t$; and $\alpha\left(t\right)$ describes the impact that the mixture of pollutants defined by

has on the probability of developing the adverse health outcome.

We design a simulation study to determine how CWVSmix performs in comparison to competing approaches with respect to (i) identification of the true critical window set, (ii) estimation of the important pollutants and interactions that comprise the time-varying mixtures, and (iii) estimation of time-varying mixture risk parameters.

For the study, we simulate data from the model in (3.1) where $\textbf{x}_{i}^{\text{T}}\boldsymbol{\beta}\equiv 0$ for all $i$ since estimation of the covariates/confounders is not the primary interest. The intercept is also set to zero so that approximately half of the simulated outcomes result in the adverse event, similar to the case-control design of our stillbirth application in NJ (Section 5). To investigate realistic settings, we simulate data that closely resemble data from the stillbirth analysis. Specifically, from (3.1) we set the number of exposure periods to 20 ($m=20$) and the number of pollutants to 5 ($q=5$). The vector of exposures specific to person $i$ across all exposure periods and pollutants,

is randomly selected from the complete NJ dataset (i.e., including all race/ ethnicities). An entire exposure profile is selected from a real individual and assigned to a simulated individual in order to maintain the temporal correlation in exposures across time and the cross-covariance among pollutants; both of which would be difficult to mimic otherwise. For each simulated dataset, we use this process to create full exposure profiles for $n=2,534$ individuals, where around half of these people ($1,267$) develop the adverse outcome in each dataset. This matches the number of stillbirths seen in the non-Hispanic Black stillbirth analysis from the real data application.

From (3.1), we modify $\lambda_{j}\left(t\right)$, $\widetilde{\lambda}_{jk}\left(t\right)$, and $\alpha\left(t\right)$ to simulate data from different scenarios, with an interest in determining how the competing methods perform as the complexity of the time-varying mixtures increases. In Settings 1, 2, 3, 4, and 5, we assume that 1, 2, 3, 4, and 5 pollutants, respectively, are important drivers of the risk due to mixtures (i.e., those pollutants where $\lambda_{j}\left(t\right)>0$). We include Setting 1 (a single important pollutant) to determine how the methods perform if mixtures are not important features of the analysis and a single pollutant approach is optimal. Within each of these settings, we include sub-settings that determine the exact composition of the time-varying mixture (i.e., choice of $\lambda_{j}\left(t\right)$ and $\widetilde{\lambda}_{jk}\left(t\right)$) and control the level of smoothness of these parameters across exposure period.

In Sub-Setting A, we assume that the same pollutants and interactions are impacting the outcome at each exposure period and at the same magnitude during each period (i.e., $\boldsymbol{\lambda}\left(t\right)$ is the same for each $t$ included in the critical window set). This represents complete smoothness of the weights across time. In Sub-Setting B, we allow for an intermediate level of smoothness by once again assuming that the same set of pollutants and interactions are important at each exposure period but allowing their relative contributions to the mixtures to change across time (i.e., only the non-zero entries of $\boldsymbol{\lambda}\left(t\right)$ change for each $t$ in the critical window set). In Sub-Setting C, we assume that the set of important pollutants/interactions and their relative contributions to the mixtures, can entirely vary across exposure periods, with only the total number of important main effects consistent across time (i.e., the zero and non-zero entries of $\boldsymbol{\lambda}\left(t\right)$ can change for each $t$ in the critical set, with only the total number of non-zero $\lambda_{j}\left(t\right)$ parameters consistent across time). This results in a complete lack of smoothness in these parameters across time.

The interaction pattern changes with each simulated dataset and is determined by the specific simulation setting. For example, in Setting 2 where only two main effect pollutants are non-zero, only the interaction between those two pollutants is eligible to be non-zero (i.e., strict hierarchy for interactions). Additionally, each eligible interaction term has a 50% chance of being active in a given simulated dataset. This allows there to be variability in the interaction pattern across the simulated datasets, resulting in a more realistic setup where only a subset of potential interactions are active. In Sub-Setting C, this selection process happens separately at each critical week so that the important interactions are varying across time. Across all settings, we simulate the non-zero weight parameters using a $\text{Dirichlet}\left(1,\ldots,1\right)$ distribution to ensure that they sum to one.

One simulated realization of the weight parameters across all simulation study settings is shown in Figure 1 where the true critical window set includes five exposure periods that we arbitrarily define as the first five weeks of exposure for plotting purposes. We note that simulating data from Setting 1B is not possible with only one important pollutant at each critical week since we can not vary the weights among this single pollutant (i.e., the non-zero weight must always be equal to one). From Figure 1, it is clear how the complexity of the mixtures increase as we move from Settings 1 to 5 and how the weights become less smooth as we move from Sub-Settings A to C.

For each simulated dataset, we randomly select how many exposure periods will be part of the critical window set (i.e., between one and seven; e.g., $m_{0}$), the first exposure period where the risk becomes elevated (i.e., between one and $m$ when possible; e.g., $t_{0}$), and define the critical window set as the consecutive $m_{0}$ exposure periods beginning on $t_{0}$. For the exposure periods selected in the critical window set, the corresponding $\alpha\left(t\right)$ parameters are fixed at $0.23$, which is the sum of the posterior means of $\alpha\left(t\right)|\gamma\left(t\right)=1$ for all exposure periods identified in the critical window set from the non-Hispanic Black stillbirth NJ CWVSmix data analysis; providing realistic effect sizes for simulation. The remaining $\alpha\left(t\right)$ parameters are set to zero. The result is simulated datasets where the true critical window sets differ by length and location during the exposure period, with consecutively defined non-zero effect sizes during this period. A similar simulation framework was used in Warren et al. (2020a).

In addition to fitting CWVSmix to each simulated dataset, we also explore two comparable alternative approaches for identifying critical exposure windows when multiple time-varying exposures are available. Both methods represent modifications of (3.1) and are given as:

• •

  Equal Weight CWVS (EW): All pollutants and interactions are given the same weight at each exposure period (i.e., $\lambda_{j}\left(t\right)=\widetilde{\lambda}_{jk}\left(t\right)=2/q(q+1)$ for all $t$, $j$, and $k$) and the original CWVS method is used to analyze the averaged time-varying exposures. Therefore, (3.1) becomes

  $$\sum_{t=1}^{m}\left[\frac{2}{q(q+1)}\left\{\sum_{j=1}^{q}\text{z}_{ij}\left(t\right)+\sum_{j=1}^{q-1}\sum_{k=j+1}^{q}\text{z}_{ij}\left(t\right)\text{z}_{ik}\left(t\right)\right\}\right]\alpha\left(t\right).$$

• •

  LWQS CWVS: At each individual exposure period, WQS regression is used to estimate $\boldsymbol{\lambda}\left(t\right)$ as $\widecheck{\boldsymbol{\lambda}}\left(t\right)$ and the original CWVS method is used to analyze the resulting weighted time-varying exposures. Therefore, (3.1) becomes

  $$\sum_{t=1}^{m}\left\{\sum_{j=1}^{q}\widecheck{\lambda}_{j}\left(t\right)\text{z}_{ij}\left(t\right)+\sum_{j=1}^{q-1}\sum_{k=j+1}^{q}\widecheck{\widetilde{\lambda}}_{jk}\left(t\right)\text{z}_{ij}\left(t\right)\text{z}_{ik}\left(t\right)\right\}\alpha\left(t\right).$$

EW serves as a baseline method to determine the impact that ignoring the time-varying weighting has on various aspects of analysis. In our version of LWQS, we opt for CWVS instead of a traditional DLM that does not include a variable selection component in order to more fairly compare findings across methods. Comparing CWVSmix with LWQS in this study allows us to determine the benefits of jointly estimating the weights and risk parameters within a single framework that accounts for temporal correlation between parameters and variable selection of the weights (CWVSmix) as opposed to a more ad hoc two-stage model fitting approach (LWQS).

We analyze the impact that exposure to multiple ambient air pollutants has on the risk of stillbirth using the data from NJ, 2005-2014, previously described in Section 2. Stillbirth is defined as the death of a baby before or during a delivery occurring after 20 completed weeks of gestation and impacts around 1 in 160 births in the United States (Hoyert and Gregory, 2016). Some known maternal risk factors include Black race, 35 years of age and older, low socioeconomic status, and smoking. Past studies have investigated the link between stillbirth and exposure to various indoor and ambient air pollution with somewhat mixed findings, though the evidence suggests further studies are warranted; see Pope et al. (2010), Siddika et al. (2016), and Bekkar et al. (2020) for recent reviews.

First, we conduct single pollutant analyses using CWVS, where each of the 12 pollutants in the study are analyzed individually and critical windows of exposure are identified. These analyses are conducted for each pollutant and race/ethnicity (i.e., non-Hispanic Black, Hispanic, and non-Hispanic White) separately. Important individual pollutants identified by the single pollutant analyses are then used in the multipollutant/mixture approaches for critical window identification by applying CWVSmix and the competing methods from Section 4 to the data.

We also include an additional competing method which uses principal component analyses (PCAs) to reduce the dimensionality of the exposures prior to analysis. Specifically, at a selected gestational week, we apply PCA to the standardized exposures from all pollutants and interactions (across all individuals in the study) and extract factor loadings from the first principal component (i.e., pollutant/interaction- and gestational week-specific weights). Using these factor loadings, we construct the factor score (i.e., weighted exposure) for each individual. This process is repeated separately at each gestational week, resulting in 20 weighted exposures for each individual. We then use CWVS to analyze the impact that these standardized PCA weighted exposures have on stillbirth risk. This method differs from CWVSmix and LWQS because the weights are not forced to sum to one, and the weight values are determined by correlation between pollutants and do not consider associations with the stillbirth outcome.

For all models, the outcome $Y_{i}$ is an indicator of a stillbirth for pregnancy $i$ and we select the first $m=20$ exposure weeks since each included pregnancy had a gestational age of at least 20 weeks. In total, we consider average weekly exposures to the $q=12$ pollutants described in Section 2 (i.e., CO, EC, NH${}_{4}^{-}$, NO₂, NO${}_{3}^{-}$, NO${}_{\text{x}}$, O₃, OC, PM₁₀, PM_2.5, SO₂, and SO${}_{4}^{2-}$). The covariates included in $\textbf{x}_{i}$ include categorical variables for the year and season of conception, an indicator of maternal tobacco use during pregnancy (yes vs. no), maternal age category ([25, 30) vs. $<$25, [30, 35) vs. $<$ 25, $\geq$ 35 vs. $<$25), maternal educational attainment (high school vs. $<$ high school, $>$ high school vs. $<$ high school), sex of the fetus (male vs. female), and the latitude/longitude of the maternal residence at delivery to account for unexplained spatial variation in stillbirth risk.

Average weekly air pollution exposures were scaled by subtracting off the median value and dividing by the interquartile range (IQR) on a specific gestational week of pregnancy, resulting in an IQR interpretation for the risk parameters. Specifically, an IQR increase from the median exposure level of each pollutant during gestational week $t$ corresponds to an odds ratio of $\exp\left\{\alpha\left(t\right)\right\}$ with respect to stillbirth risk for EW, LWQS, and CWVSmix (see Section 3.3 for more information). For PCA, because the weight parameters do not sum to one, the corresponding effect is given as

where $\lambda_{j}^{(\text{pca})}$ are the PCA main effect weights (similar definition holds for the interaction parameters). In Figures S1-S3 of the Supplement, we display the IQRs for each pollutant and race/ethnicity across every gestational week for the purposes of interpretation.

For CWVSmix and each application of CWVS, we collect 10,000 posterior samples from the joint posterior distribution after discarding the first 10,000 during a burn-in period and thinning the next 100,000 by a factor of 10 to reduce posterior autocorrelation. Trace plots of individual parameters and the Geweke diagnostic (Geweke et al., 1991) were used to assess convergence of each model, with no obvious signs of non-convergence observed. Prior distributions and model fitting details match those from Section 4.1.

In this work, we introduced CWVSmix, a method to simultaneously identify critical windows of exposure to multiple time-varying exposures, determine the time-varying relative importance of individual pollutants/interactions within the mixture, and estimate risk associations for exposure mixtures. Through a simulation study, we showed that CWVSmix offers the best balance of performance across each of these areas with respect to the competing methods. In particular, LWQS, which attempts to estimate these mixture weights and associations at each exposure time period separately, is shown to produce critical window identification accuracy, estimates of mixture weights, and estimates of risk that are always less optimal than CWVSmix. This suggests that jointly estimating the weights and risk parameters and accounting for temporal correlation in exposures while identifying critical windows may be important features of an analysis. Additionally, the variable selection procedure introduced for the weight parameters allows for a more intuitive explanation of the drivers of risk during a specific exposure period.

By introducing weights that sum to one and are between zero and one, we assume that (i) pollutants impact the risk in the same direction on a given exposure period and (ii) main and interaction effects are in the same direction. For (i), we do not enforce that the global risk parameter must be positive or negative, which avoids bias in these estimates when that assumption is violated. For (ii), we anticipate that if two main effects have adverse (or protective) impacts on risk, then their interaction will either intensify this relationship or be negligible. The benefit of this formulation is the interpretability of the findings when describing the composition of a harmful mixture, similar to source apportionment methods, and that it can effectively collapse to a single pollutant model given that a main effect weight can be exactly equal to one.

In our application to stillbirth risk due to multiple ambient air pollution exposures, we found elevated risk of stillbirth for non-Hispanic Black mothers who experience higher concentrations of NH${}_{4}^{+}$, NO${}_{3}^{-}$, NO${}_{\text{x}}$, PM_2.5, and SO${}_{4}^{2-}$. Previous work in this area has led to mixed findings regarding the harmful pollutants and exposure periods, with several works relying on single pollutant approaches fit across multiple exposure periods. Using CWVSmix, we can now more accurately identify the critical exposure periods while accounting for exposure mixtures. Additionally, use of the fused air pollution estimates gives us the opportunity to investigate less commonly studied pollutants.

The differences observed between the competing methods in the data application must be interpreted in the context of the simulation study conclusions when determining which results are likely most reliable. Findings in Table 1 from the simulation study, which was designed to resemble the NJ stillbirth data application, indicate that CWVSmix is more likely to accurately identify the true critical window set, yield improved estimates of the weight parameters, and yield improved estimates of the risk parameters when the proportion of important pollutants/interactions is relatively small. The variable selection findings from the NJ stillbirth analysis in Table 3 suggest that the data application most closely resembles Settings 1 and 2 from the simulation study, where CWVSmix outperformed the competing methods across all metrics (Table 1). Additionally, the protective association estimated using LWQS in the Hispanic data application (but not CWVSmix) is likely due to the lack of smoothness in the estimated weights, and provides some additional support to suggest that CWVSmix may yield more reliable results in this case. Because WQS is fit separately at each exposure period in the LWQS algorithm, the model remains unadjusted for the other time-varying exposures when the weight parameters are estimated. This may also contribute to the findings observed in the simulation study and NJ data applications.

The ambient air pollution estimates from the data fusion model represent a major strength of this work. Typically, only a handful of regularly monitored pollutants are able to be included in a particular analysis, depending on the timeframe and spatial location of the study. In our case, we had excellent space-time coverage of our estimates, allowing us to investigate associations that are less common, with increased sample sizes. However, for some of the pollutants (e.g., PM_2.5 constituents) the measurements in space-time are limited and the model-derived estimates further away from the measurements may be more biased than those near the measurements. In addition, the SO₂ estimates have substantial uncertainty due to the plume nature of these pollutant fields that is difficult to capture with either monitoring or modeling.

As is typically the case when linking environmental exposures to birth records data, we used residence at delivery to define the spatial location of data linkage. This may introduce measurement error for some participants in the study, though Warren et al. (2017) showed that critical window identification is relatively robust to this type of error. In (3.6), we assume independence a priori between the weight parameters within an exposure time period due to the potentially large dimension of this matrix (i.e., $q(q+1)/2$ by $q(q+1)/2$) as the number of pollutants increases. Future methods work that could reduce this dimensionality, possibly through incorporation of auxiliary information about the chemicals in a mixture, may be able to accommodate this correlation (Reich et al., 2020). Future work may also want to consider interactions between pollutants on different exposure periods, which will require further methodological extensions.

When determining which method to use in future multipollutant applications, we recommend designing a simulation study based on important features of the analysis dataset and comparing the competing approaches. Additionally, computing times may be important to consider based on the size of the dataset. In Table S7 of the Supporting Information, we show the average computing times in minutes for each of the methods in the simulation study. EW, which relies on CWVS, is clearly the fastest algorithm, followed by CWVSmix and LWQS. The LWQS run times are impacted by the number of bootstrap samples requested and the number of exposure periods where the WQS algorithm must be applied. Overall, CWVSmix offers an intuitive and interpretable method for analyzing the impact of multiple time-varying exposures on an adverse health outcome.