Bag of DAGs: Flexible Nonstationary Modeling of Spatiotemporal Dependence

We propose a class of nonstationary processes that characterize varying directional associations in space and time for point-referenced data. Our construction is based on local mixtures of sparse directed acyclic graphs (DAGs). In stochastically choosing DAG edges from a"bag,"we account for uncertainty in directional correlation patterns across a domain. The resulting"bag of DAGs"processes (BAGs) lead to interpretable nonstationarity and scalability for large data due to sparsity of all DAGs in the bag. We are motivated by spatiotemporal modeling of air pollutants in which a directed edge in a DAG represents a prevailing wind direction causing some associated covariance in the pollutants. We outline Bayesian hierarchical models embedding the resulting nonstationary BAGs and illustrate inferential and performance gains of our methods compared to existing alternatives. We analyze fine particulate matter (PM2.5) in California with high-resolution data from low-cost air quality sensors. The code for all analyses is publicly available at https://github.com/jinbora0720/BAG.

Introduction
In the spread of aqueous or air pollutants, local dynamics of currents or winds influence the strength of spatial and temporal correlations. Figure 1 (left) illustrates an incident in California (CA) in which smoke from wildfires exhibits local patterns of spread in different directions. Realistic models of pollutant spread should allow correlations between spatial locations s and s to peak if one is downwind of the other after a certain temporal lag. In principle, one can construct rich models for any point-referenced data including pollutants via Gaussian processes (GPs). Let w(s) ∼ GP (0, C(·, ·|θ)) denote a zero-mean univariate  GP where s ∈ D ⊆ R d is a location in a region of interest D, and d is the dimension of coordinates (e.g., d = 2 if s = (longitude, latitude)). A covariance function C(·, ·|θ) : D × D → R specifies associations across locations, indexed by a parameters θ. Then, finite realizations of w(·) at an arbitrary set of locations S = {s 1 , . . . , s k } follow a multivariate normal distribution with mean zero and a covariance matrix C S whose (i, j) element is C(s i , s j | θ) = Cov(w(s i ), w(s j )).
Customarily GPs assume a stationary covariance function such that C(s, s + h) = C(s , s + h) = C(h) for any s, s ∈ D, omitting θ for brevity. The stationarity assumption, however, is overly restrictive in our applications of interest because it assigns the same correlation to any pairs of locations with the same lag h without regard to their respective and relative positions. To overcome the limitation, GPs can rely on nonstationary covariance functions (see Sampson 2010). Moreover, when spatially-and/or temporally-varying covariates are available, one can explicitly model covariance parameters using the covariates (Calder, 2008;Schmidt et al., 2011;Neto et al., 2014;Risser and Calder, 2015). Calder (2008) and Neto et al. (2014) use wind direction data to characterize nonstationarity in air pollutants. However, reliable covariate data (e.g., wind directions) that are recorded at a sufficiently high spatial and/or temporal resolution are seldom available. In this paper, we aim to construct nonstationary processes that do not require nonstationarity-informing covariates but take into account effects of winds or currents.
Another hurdle in nonstationary GPs is poor scalability to large data. Some computationally efficient nonstationary GP models use the basis graphical lasso (Krock et al., 2021a,b) or fit locally stationary processes for different regions of a domain, allowing for different covariance parameters (Fuentes, 2001;Kim et al., 2005;Lindgren et al., 2011;Kleiber and Nychka, 2012). Several extensions (Datta et al., 2016a;Katzfuss and Guinness, 2021) of Vecchia (1988) approximation of stationary GPs also induce nonstationary GPs. Vecchia-style approximations achieve scalability by making spatial conditional independence assumptions across locations or knots. Each knot is associated to a node in a sparse directed acyclic graphs (DAG), and the edges pointing to it populate its parent or conditioning set. One can construct a Nearest-Neighbor GP (NNGP; Datta et al. 2016a) by restricting the conditioning set to include a few closest neighbors by Euclidean distance or a Meshed GP (MGP; ) by fixing a patterned DAG over domain partitions. Further generalizations are available in Katzfuss and Guinness (2021). DAG-based methods are appealing because they can be implemented with any covariance functions, can lead to standalone processes that facilitate predictions at arbitrary locations, and can be embedded seamlessly within Bayesian hierarchies. Multiple extensions have been proposed for data that are spatiotemporal (Datta et al., 2016b), multivariate (Peruzzi and Dunson, 2020;Zhang and Banerjee, 2021), or non-Gaussian (Peruzzi and Dunson, 2022).
Conditional independence encoded via predetermined DAGs, however, may inhibit inferences on certain types of nonstationarity. In the middle and the right plots of Figure 1, fixed DAGs respectively from nearest neighbors and adjacent partitions disregard the local dynamics in the spread of wildfire smoke and air pollutants. The middle plot illustrates a simple scenario of an NNGP in which four nearest neighbors are connected by pink arrows to a location approximately at (120 • W, 37 • N). Three of these neighbors, however, are outside the smoke above the location of interest and may be undesirably non-informative despite proximity. Similarly, on the right plot, a Cubic MGP (Q-MGP) model fails to capture locally varying associations due to the fixed and repeated pattern in its DAG. These exemplify potentially incorrect assumptions about spatial and/or temporal dependence previous DAGbased models rest on, which can also result in unreliable estimates of covariance parameters and predictions.
In this paper, we propose a novel class of nonstationary processes characterizing varying directional dependence structures in space and time. Conditional independence between partitions of a domain is encoded in a sparse DAG whose directed edges are allowed to vary locally to capture different directional dependencies across the domain. We call this a Bag of directed Acyclic Graphs process (BAG). One may choose observed locations as nodes of DAGs instead of partitions; however, it is unwieldy to infer directions at the level of observed locations unless they lie on a grid. Therefore, we simplify the problem by partitioning the domain into disjoint regions, each of which becomes a node. Importantly, BAGs are valid standalone processes, leading to inferential advantages by incorporating parameter estimation and prediction at arbitrary locations into a coherent framework. BAGs can be employed as a prior in Bayesian latent spatial/spatiotemporal process models and can scale to large data (the number of observations n > 10 5 ). Unlike a recently proposed DAG-based nonstationary model of Gaussian data with replicates (Kidd and Katzfuss, 2021), BAGs do not rely on approximations in posterior sampling and can be easily extended to multivariate non-Gaussian settings following Peruzzi and Dunson (2022). Furthermore, BAGs induce interpretable nonstationarity.
Our construction of BAGs is partly motivated by causal inference in which DAGs have been used as systematic representations of causal relationships (Wright, 1934;Greenland et al., 1999;Shrier and Platt, 2008;Textor et al., 2016). We propose a related idea using a local mixture of directed edges. Directed edges in a DAG are inferred probabilistically and interpreted as directional associations across a domain. Mixture-related methods have been proposed for added flexibility (Griffin and Steel, 2006;Duan et al., 2007) or to generalize underlying stationary GPs to nonstationarity and beyond Gaussian assumptions (Gelfand et al., 2005;Rodríguez et al., 2010). Distinctively, our mixture is local and involves smaller regions of the domain one-at-a-time, instead of the whole domain at once. In the context of air pollution, local mixture components for each region represent different prevailing wind directions which may explain associated local covariance in the pollutants. Through stochastic selection of arrows, BAGs account for uncertainties of directional dependence patterns.
Since different prevailing wind directions could be selected for different regions of the domain, BAGs are flexible enough to mimic the potential volatility of wind directions.
The rest of the paper is structured as follows: Section 2 formulates BAGs and Gaussian BAGs (G-BAGs) as a special case, followed by an empirical investigation on nonstationarity of G-BAG induced covariance functions. We demonstrate the performance and inferential benefits of our proposed approach via simulation results (Section 3) and application results to air quality data in CA (Section 4). Section 5 concludes the paper with useful discussions and possible extensions. Supporting Material is separately available including detailed derivations of the properties of BAGs, estimation and prediction steps, computational costs, and more results from simulation studies and data applications.

Spatial Process Modeling Using BAGs
Consider a univariate spatial process, {w(s) ∈ R | s ∈ D ⊆ R d }. The following discussions are identically applied to spatiotemporal processes defined in D ⊆ R d+1 , so we may proceed with "spatial" processes for brevity. We first divide the domain D into M disjoint partitions, . . , s k } denote a fixed "reference set" set of locations. Here, S need not coincide or intersect with observations T = {t 1 , . . . , t n }, but T is a practical choice for S nonetheless.

Fixed DAG
Consider finite realizations of the latent process w(·) over S, w S = (w T 1 , . . . , w T M ) T where w i = (w(s i 1 ), . . . , w(s i k i )) T are realizations over partition S i for i = 1, . . . , M . After choosing an arbitrary order of the M partitions, the joint density p(w S ) = p(w 1 , . . . , w M ) can be rewritten as a product of conditional densities: (1) We can represent (1) as a DAG G S = (A, E) with nodes A = {a 1 , . . . , a M } and directed The k i random variables w i are collectively mapped to a single node a i . As an example, realizations of w(·) at all of the blue dots in Figure 2 are assigned to the node a 1 . A set of nodes [a i ] is a subset of A whose edges are directed to a i and is referred to as a parent set.
Motivated by Vecchia approximation, we can drop some edges in E and build a new DAG G * S that leads to a distinct joint density p * as a product of conditional densities with reduced conditioning sets Here, w [i] is a subset of {w 1 , . . . , w i−1 } collecting realizations at parent nodes of a i , w (Lauritzen, 1996). Approximating p with p * coupled with a DAG renders advantages: (a) the approximation is highly accurate in Kullback-Leibler divergence; (b) accuracy can be improved with a careful order of nodes (Guinness, 2018); (c) the DAG can be extended to arbitrary locations to create standalone stochastic processes (Datta et al., 2016a); and (d) DAGs can be designed for faster computations . In these cases, a pre-specified DAG is introduced only for computational reasons.

Unknown DAG
Taking a different perspective, we allow the DAG to be inferred from data. We account for uncertainty about dependency patterns across the spatial domain by considering finite mixtures of possible local arrows within each partition. A simplifying assumption is made that each node can have at most one parent from neighboring partitions in the spatial domain. This assumption may be relaxed by letting the number of parents be > 1, but at the cost of reduced clarity in exposition, computational tractability, and interpretability.
Taking toy examples in Figure 2, directed edges 1, 2, 3, 4 (left) or 1, 2, 3 (right) correspond to competing assumptions on the dependence structure relevant for the shaded partitions given the tessellations using rectangles or hexagons, respectively. Such tessellations of the domain are advantageous because directed edges lead to intuitive interpretations, e.g. the directed edge 1 on the left of Figure 2 is interpreted as the direction from west to east.
Moreover, stochastic search of directed edges will proceed analogously for all partitions due to the same neighbor structure and thus the same "bag." For node a i (equivalently, reference subset S i ), we introduce a latent membership variable z i that determines the direction along which dependence is allowed to flow and thus the parent node. After enumerating possible directions to each node as h = 1, . . . , K, we let p(z i = h) = π ih be the probability that a i receives the hth directed edge. We denote the resulting parent node as [a i |z i = h]. We include only half of the possible directions in Figure 2 to ensure acyclicity of inferred directed graphs.
Prior knowledge about the domain may be utilized to choose one direction over the other on each axis. For a spatiotemporal domain, we assume that each node receives an additional directed edge from the partition that covers the same spatial region at the previous time point, and a spatial and the temporal parent are collectively denoted as [a i |z i = h].
Using these assumptions, we define a joint densityp conditional on a specific configuration of directed edges z = (z 1 , . . . , z M ), modifying equation (2) as follows: . . , a i−1 }}. Since z 1 , . . . , z M are independent a priori, integrating them out gives For notational brevity, z is hereafter z∈bag of DAGs unless otherwise specified. Then proposition 1 holds whose proof is provided in Supporting Material S1.
Proposition 1.p(w S ) is a proper joint density.

BAGs
The above discussion focuses on finite dimensional w S . We extend our approach to a valid process over all locations in the domain. We label non-reference locations as U = D \ S, and We extend the DAG over S to a larger DAGG, with each of the U i 's mapped to a node b i ∈ B. The construction ofG is completed by assigning directed edges from nodes in A to those in B, ensuring acyclicity of G. There are several possible ways to place these directed edges. Assuming that S i = ∅ for all i, one can fix the edges as a i → b i ([b i ] = a i ), implying that local reference locations become the parent set for the non-reference locations in the same partition. Kolmogorov consistency conditions are then easily verified (see Appendix A in , proving the validity of the resulting process. However, a better alternative lets [b i |z i ] = a i ∪ [a i |z i ], which allows modeling at any non-reference locations to depend not only on the local reference locations but also on their parent locations learned from data in our approach. In this way, predictions are made based on inferred wind directions. We prove the Kolmogorov consistency conditions for our proposed stochastic process when there is randomness in choosing the DAG for both A and B.
First, we assume conditional independence of non-reference locations given the local reference locations and their parents so that The equations (3) and (4) suffice to describe the joint density of any finite subset of spatial We confirm that the collection of finite dimensional densities in equation (5) satisfy the Kolmogorov conditions (see Supporting Material S2). This implies that there exists a stochastic process associated with them, which we call a BAG. Therefore, our approach generates a valid spatial process via domain partitioning and local mixtures of DAGs.
When modeling spatial variation, we are interested in understanding the covariance between process realizations at different spatial locations. With BAGs, we find that if Ep(w(l) | z) = 0 for any l ∈ D, where Ep and E z indicate expectation with respect tõ p and z, respectively. This suggests that pairwise covariances between process realizations can be studied starting from the conditional covariances given a DAG configuration z. We further study BAGs in the Gaussian case and the induced nonstationary covariance.

Gaussian BAGs
Suppose the spatial process {w(s)} is a zero-centered GP, w(s) ∼ GP (0, C(·, ·|θ)) with a valid covariance function C parametrized by θ. Then equation (3) becomes where The latent membership variable z determines the spatial DAG which in turn defines H i|z i and R i|z i . We also findp(w S | = a j and zero otherwise, and a block-diagonal matrix R = diag(R 1|z 1 , . . . , R M |z M ). Given a bag of directed edges, it is always possible to in a lower triangular matrix H z with zero diagonals. This renders |I k − H z | = 1, and thus Furthermore, in the spatial only case, the ith block-row H z [i, :] of H z has at most one non-zero block for any i due to the assumption that all nodes inG S (z) have at most one parent. As a result, the undirected moral graph found by "marrying the parents" of nodes inG S (z) generates no additional edges. This means we can immediately identify the (i, j)th block ofC −1 z as non-zero if and only if there is a directed edge between a i and a j inG S (z) (either a i = [a j |z j ] or a j = [a i |z i ]). Hence, the sparsity structure of lower-triangular ) T by symmetry. In spatiotemporal cases, the undirected moral graph will include additional edges connecting a spatial parent and a temporal parent for each node inG S (z). As a result,C −1 z will have a non-zero block at (i, j) not only if there is a directed edge between a i and a j , but also if a i and a j have a common child node (one as the spatial directional parent and the other as the temporal parent).
Similarly, for non-reference locations, equation (4) becomes In vector form, the joint density of w U is multivariate Gaussian with mean H U |z w S and covariance matrix R U |z in which H U |z is a |U| × k matrix with |U| being the number of locations in U, and R U |z is a diagonal has non-zero blocks in columns corresponding to a i and a j , and the ith block matrix of R U |z is also a diagonal matrix whose elements are R u|z i ∀u ∈ U i . Following equation (5) and the subsequent discussion, we obtain a Gaussian BAG (G-BAG) such that w(s) ∼ G-BAG(0,C(·, ·|θ)). The new covariance functionC for G-BAG is given asC(l 1 , l 2 | θ) = zC (l 1 , l 2 | z, θ)p(z) for any two locations l 1 , l 2 ∈ D by equation (6).C(l 1 , l 2 | z, θ) is the induced covariance function of a GP one obtains from fixing a single DAG as z: since equations (7) and (8) are all Gaussian, the densityp(w L | z) is also Gaussian for any finite subset L ⊂ D, leading to a GP for any given z. The implied covariance function conditional on z is theñ whereC P,Q is a submatrix ofC z corresponding to locations in the sets P and Q, and 1(·) is the indicator function. HereC(·, ·|z, θ) is nonstationary, as is the marginal covariance functionC(·, ·|θ). Derivation of equation (9) is described in Supporting Material S3.

Bayesian Hierarchical Model with G-BAGs
Consider a general regression model where y(t) is a response variable, x(t) is a p-dimensional vector of point-referenced predictors, w(t) is a spatial process over domain D, and (t) ∼ N (0, τ 2 ) is a measurement error. Over the observation set T = {t 1 , . . . , t n }, equation (10) is expressed in vector form as y = Xβ + w + with y = (y(t 1 ), . . . , y(t n )) T , w and similarly defined, a n × p matrix X having x(t i ) T as its ith row, and ∼ N (0, τ 2 I n ). We assume w(·) ∼ G-BAG(0,C(·, ·|θ)) a priori as specified in Section 2.4. We complete our Bayesian hierarchical model by specifying priors for all unknowns. The joint posterior distribution for θ is proportional to the likelihood times the prior distributions: where U = T \S, and Cat, Dir, and IG denote Categorical, Dirich-let, and inverse Gamma distribution, respectively. The vector π i = (π i1 , . . . , π iK ) stores the probability for each possible value of z i , and α = (α, . . . , α) T is the K × 1 Dirichlet hyperparameter vector. The prior p(θ) is left unspecified as it depends on the choice of a base covariance function C.
Spatiotemporal G-BAGs can use any spatiotemporal covariance function as a base. As a flexible default, we focus on the nonseparable (stationary) spatiotemporal covariance function of Gneiting (2002). For space-time lag (h, u) ∈ R d+1 , the covariance function is where κ ∈ [0, 1] is a space-time interaction parameter, a > 0 and c > 0 are temporal and spatial decays, respectively, and | · |, · denote 1, 2 dimensional Euclidean distance. (11) reduces to a separable function when κ = 0. We discuss nonstationarity induced by our approach from this stationary base covariance in Section 2.6. A straightforward Markov chain Monte Carlo (MCMC) sampler to obtain posterior samples is provided in Supporting Material S4. We show that G-BAGs have computational complexity of order n at each iteration of the MCMC sampler in Supporting Material S5.

Nonstationarity of G-BAGs
Even with a stationary base covariance, G-BAGs induce directional nonstationarity. Here, we consider a simplified setting to elucidate the nonstationarity features ofC(·, ·|θ) among   (a) Covariance heat maps between a reference point (red dot around (0.5, 0.5, 0.333)) and the others. G-BAG nonstationary covariances are depicted with two partitioning schemes (first two rows), using the base stationary covariance on the bottom row.

Applications on Simulated Data
We conduct simulation studies to evaluate performance of G-BAGs in prediction and learning DAGs. We consider a scenario in which a G-BAG model is correctly specified, and one in which a fitted G-BAG model is misspecified relative to the data generating model. In both scenarios, our spatiotemporal domain is D = [0, 1] 3 with locations t = (t 1 , t 2 , t 3 ) ∈ D. For G-BAG models, we assume that a latent process w(t) follows G-BAG(0,C(·, ·|θ)) a priori with the base covariance function in equation (11) with θ = (a, c, κ, σ 2 ).
We compare our G-BAG approach to Q-MGP whose cubic meshes are considered as fixed  (Rue et al., 2009) is used to implement SPDE models, whose basic interface only supports separable spatiotemporal models. Hence, an autoregressive model of order 1 (AR (1)) is assumed for temporal correlations, while nonstationarity is imposed in spatially varying spatial range such that ρ(t) ∝ exp(ψ 1 +ψ 2 t 1 +ψ 3 t 2 +ψ 4 (t 1 −1)(t 2 −1)). The reference sets for G-BAG and MGP models are the observed locations which are 80% of the simulated locations. Twenty percent of data are held-out for validation. The number of threads is set at 10; different methods benefit from multiple threads by different degrees.
Parameters β and τ 2 are set at 2 and 0.01, respectively. We partition each of the spatial axes   This suggests that when directional associations are believed to shift rapidly over space and time, it is recommended to use G-BAGs that explicitly model directional dependence structures. Its computational time (in minutes), however, is longer than the competitors largely due to difference in coding languages and mixtures that bring heavier computational burden.
Although we believe the current running time is acceptable, we expect some improvement in computational time by fully optimizing the code as a future work. The coefficient β is estimated correctly by all models. The SPDE models considerably overestimate τ 2 , while G-BAG and Q-MGP produce τ 2 estimates close to the truth. The decay parameters a and c are identified by G-BAG using weakly informative priors. The nostationary SPDE barely improves the stationary SPDE in terms of prediction. This may imply that when directional associations are prevalent, the nonstationarity specified via a location-specific spatial range does not suffice.

Fitted G-BAG is misspecified
In the misspecified case, we generate data on a regular lattice of size 193 × 193 × 59 assuming a univariate model y(t) = w(t) + (t) with (t) ∼ N (0, τ 2 ) and τ 2 = 0.1. We create w(t) with north arrows fixed over 1 × 193 × 59 partitions to mimic a surface under constant wind directions. The following covariance function introduced in Gneiting (2002) is used as a base: ν is the smoothness parameter of space, and K ν is the modified Bessel function of the second kind of order ν. Fixing ν = 1.5, we generate 25 synthetic data sets with θ 3 = (5, 20, 1, 150) and another 25 with θ 4 = (10, 20, 1, 150). The latter mimics faster wind speed than the former. Figure  Prediction results for the 25 data sets with θ 3 are illustrated in Figure S1 in Supporting better predictive performance of Q-MGP than the G-BAG model can be attributable to (a) the data generating DAG being fixed and (b) two spatial parents in Q-MGP versus one spatial parent in G-BAG. It is encouraging that G-BAG has comparable performance to Q-MGP using less information from parent nodes, by choosing one "best" spatial parent based on inferred directions. We observe similar results with θ 4 which are omitted to save space.
Moreover, G-BAG enables inferences on directions along which correlations move across the domain. Figure S2 in Supporting Material S6.2 demonstrates that the true direction (N) is properly found among four directions in the bag by the highest posterior probability in over 80% of the space-time partitions (left). Across time, each spatial partition assigns N an average posterior probability of around 0.7 (right).

Air Quality Data Analysis
Among many ambient pollutants, fine particulate matter that is 2.5 microns or less in diameter (PM2.5) is the primary concern due to its abundance and association to adverse health effects. In particular, acute exposures to PM2.5 have been associated with detrimental effects on human health including heart/lung diseases and associated premature deaths to various health effects and its influence on daily activities, accurate prediction of PM2.5 in regions without monitoring sites is an imperative task. We focus on CA (Section 4.1) and South Korea (Supporting Material S7) whose residents suffer from acute exposures to high PM2.5 seasonally due to wildfires and long-range movement of air pollutants, respectively.
The spatiotemporal spread of PM2.5 is heavily affected by winds. A large body of literature has included wind-relevant variables in an attempt to predict PM2.5 (Wu et al., 2006;Calder, 2008;Wang and Ogawa, 2015;Preisler et al., 2015;Kleine Deters et al., 2017;Aguilera et al., 2020). However, there is a fundamental issue in summarizing wind direction at discrete times. The ambiguity of representative values for wind direction is evident in a tendency for different data sources to provide different summaries. For instance, EPA provides average wind direction, while the National Oceanic and Atmospheric Administration (NOAA) provides direction of the fastest wind. Moreover, due to the volatile nature of winds, naive averages of such wind directions over a given time period may rarely be meaningful.
Therefore, we use the G-BAG approach on the latent process of log transformed PM2.5 to overcome these difficulties and implicitly incorporate wind effects.

California, the United States
The year 2020 was the largest wildfire season in CA's modern history. During fire events, dramatically poor air quality is witnessed due to wildfire emissions in which PM2.5 is the primary pollutant (Liu et al., 2017). These environmental risks have made CA outstanding in its widespread use of low-cost sensors such as PurpleAir. In particular, more than half of the PurpleAir sensors in the United States are concentrated in CA as of February, 2020 (deSouza and Kinney, 2021). Therefore, we analyze daily PM2.5 levels in CA during fire seasons (August 1 to October 22, 2020) using EPA and PurpleAir measurements. Low-cost sensors require appropriate correction for improved comparability to regulatory monitors; e.g., using the Federal Reference Method (FRM) or Federal Equivalent Method (FEM) (Datta et al., 2020). Using recommended practice by EPA Evans et al., 2021), PurpleAir data are calibrated relative to FRM/FEM via: P M 2.5 = 5.75 + 0.52P A CF 1 − 0.09RH if P A CF 1 ≤ 343µg/m 3 and P M 2.5 = 2.97+0.46P A CF 1 +3.93×10 −4 P A 2 CF 1 otherwise. P M 2.5 is the calibrated value, P A CF 1 is the PurpleAir measurements from their correction factor labeled as CF = 1, and RH is relative humidity in percentage.
The model is y(t) = x(t)β + w(t) + (t) where x(t) is the Euclidean distance to the nearest fire at spatiotemporal location t, and y(t) is the mean centered log(PM2.5). We expect x(t) to capture elevated PM2.5 level due to proximity to wildfires. The covariate x(t) is standardized to have mean 0 and standard deviation 1. A total of 110,473 spatiotemporal locations are observed of which 20% are held out for validation. We additionally predict at 274,564 new locations. The (W, NW, N, NE) edges are chosen in the bag because CA lies within the area of prevailing westerlies. CA is partitioned by M = 320 rectangles, and each covers approximately 55 2 km 2 (0.5 • resolution). With 8,000 burn-in and 2 thinning, 1,000 posterior samples are analyzed. We compare G-BAG to Q-MGP and SPDE-nonstationary.
Detailed prior specifications are available in Supporting Material S8. The total run time for G-BAG, Q-MGP, and SPDE-nonstationary with 10 threads is roughly 22 hours, 32 minutes, and 4 hours, respectively. Table 2 summarizes the results in comparison to Q-MGP and SPDE-nonstationary. All parameters are similarly estimated by G-BAG and Q-MGP, and the effect of the distance to the nearest fire appears insignificant in that the 95% CIs of β from both models contain zero. Although the fitted β is significantly negative from SPDE-nonstationary, it is likely misled due to unexplained space-time variations after fitting spatiotemporal random effects (see Figure S4 in Supporting Material S8). Much larger τ 2 fitted by the SPDE model than that by G-BAG also corroborates this remaining variability. These suggest that the nonstationarity via location-specific spatial range is insufficient to fully characterize the cause and removal of air pollution affected by wind directions in CA over this time period.
From Table 2, G-BAG yields smaller errors than the other two models and better accuracy of uncertainty quantification in prediction across the whole period compared to Q-MGP. A detailed comparison is illustrated in Figure S5 in Supporting Material S8 by which we confirm that the G-BAG approach improves prediction by a large margin at each time point especially over the fixed-DAG approach. We achieved down to 60% reduction in prediction errors.
These results signify that the unknown DAG approach and construction of nonstationarity via the inferred directions help enhance performance in this application.
The predicted surfaces of log(PM2.5) and observations are presented in Figure 6. We confirm that G-BAG can produce plausible predicted values even in regions lacking data with properly increased prediction uncertainty ( Figure S6 in Supporting Material S8). Moreover, we find that W and N are the most prevailing directions causing associations in PM2.5 in CA from August to October, 2020. They are selected in 70% of the space-time partitions.
In particular, the bottom row of Figure 6 shows that the growth and reduction of areas with poor PM2.5 levels (> 35µg/m 3 ) progress from W and N directions, and the posterior mode of arrows captures the tendency well. Between latitude 36 • N and 40 • N and longitude 121 • W

Discussion
We have developed a class of nonstationary processes based on flexible selection of DAGs for large point-referenced data. The interpretability of our model lies in its construction of nonstationarity and recovery of true directions among directed edges in the bag. Our model can manage high-resolution data at individual household level (e.g. PurpleAir) which may also be at road (e.g. Google Street View vehicles; Guan et al. 2020)   made. Although we have focused on univariate applications motivated by PM2.5 analysis, multivariate extensions are immediate with a proper cross-covariance function. In addition, as GPs were used to derive G-BAGs from the more general BAG approach, one may use other processes. For instance, Student-t processes can derive T-BAGs, and such heavier-tail processes increase robustness to outliers.
Our model has broad applicability in several research areas. First, it can be used for inferring environmental risks at individual level which is then linked to various health outcomes. Such imputation is highly valuable because measurements of environmental risks and health outcomes are often misaligned. Second, it can simulate realistic processes of air pollutants under various wind scenarios and provide mechanistic insights of how air quality would behave with abnormal events such as fires. This would contribute to a body of literature assessing impacts of wildfires or smoke plumes on air quality, which currently relies on separate modeling approaches with or without fires (Wu et al., 2006), regressions without spatial effects (Preisler et al., 2015;Larsen et al., 2018) or inverse distance weighting interpolations (Aguilera et al., 2020). Finally, in a multivariate setting, the stochastic selection of DAGs may be applied not only for spatiotemporal dependence of multiple outcomes, but also for the dependence among the outcomes. The latter can be useful in community ecology applications, for example, where one could use BAGs to scalably infer the relationships among species, choosing from DAGs based on prior information about species' phylogeny.
Our current implementation of BAGs ensures acyclicity by choosing directed edges from a subset of all possible directions (e.g., either W or E, either N or S, etc.). In addition, we restrict each node in a DAG to have at most one spatial parent. These assumptions help reduce the computational burden caused by the stochastic search of the DAG space.
However, they may be too restrictive in some contexts. Therefore, work needs to be done to relax these assumptions in a coherent framework, keeping computational burden low.
Another interesting future research is to develop a more realistic model of the time dynamics, rather than fixing a directed edge from a previous time point to the current, for better predictions in longer time lags.
Supporting Materials for "Bag of DAGs: Flexible Nonstationary Modeling of Spatiotemporal Dependence" S1 Proof of Proposition 1 Proof. The joint density in equation (3) is proper because there exists an associated DAG G S (z) whose nodes are A and edges are identified by z. Then . . , M . Hence,p(w S ) is a proper joint density.

S2 Kolmogorov Consistency Conditions for BAG
Kolmogorov consistency conditions indicate consistency with (a) permutation of indices and (b) marginalization. First, we show thatp(w L ) =p(w Lπ ) for any permutation L π of L. Let U Lπ = L π \ S with the fixed reference set S. Given fixed partitions, because U L and U Lπ , S \ L and S \ L π share the same collection of elements, respectively.
Second, we show thatp(w L ) = p(w L∪{l 0 } )dw(l 0 ) for a new location l 0 ∈ D \ L. Take L 1 = L ∪ {l 0 }. We consider two cases: when l 0 ∈ S and when l 0 / ∈ S. If l 0 ∈ S, then w S , z) due to the conditional independence assumption at non-reference locations given the reference set. Note that a finite sum is interchangeable with an integral.

S3 Conditional Nonstationary Covariance Function
The conditional nonstationary covariance function in equation (9) is computed as follows.
Covp and Ep implicitly depend on θ. When l 1 = s i ∈ S and l 2 = s j ∈ S, the result is trivial.
When l 1 is not in S but l 2 is, i.e. l 1 ∈ U i and l 2 = s j ∈ S, When neither l 1 nor l 2 is in S, i.e. l 1 ∈ U i and l 2 ∈ U j ,

S4 Estimation and Prediction of G-BAGs
The straightforward Markov chain Monte Carlo (MCMC) sampler to obtain posterior samples with a general C is provided below: (a) Update β from for h = 1, . . . , K and i = 1, . . . , M .
The response y * i is a k i × 1 vector whose jth element is y(s i j ) if s i j ∈ S * i or 0 otherwise, leaving n i non-zero elements. Similarly, X i is a k i × p matrix with zeros at rows corresponding to locations outside and H u,i is a submatrix of H u|z j by choosing columns that correspond to w i given that (g) Finally, θ is updated from a Metropolis-Hastings step with target density (i) In the covariance function in equation (11), θ is (a, c, κ, σ 2 ).
(ii) We use the robust adaptive Metropolis algorithm proposed by Vihola (2012) for g(a, c, κ) with target acceptance rate of 0.234. A link function g maps (a, c, κ) to (−∞, ∞) range. Uniform priors are used.
(iii) A Gibbs step to update σ 2 can be easily derived with a prior σ 2 ∼ IG(a σ , b σ ).
Posterior predictive samples at an arbitrary location l ∈ D are obtained as follows: at each MCMC iteration, sample y(l)|y ∼ N (x(l) T β + w(l), τ 2 ) if l ∈ S ∪ T . If l / ∈ S ∪ T , then first sample w(l) from N (H l w [l] , R l ) conditioned on w S and z if desired, and then sample y(l)|y ∼ N (x(l) T β + w(l), τ 2 ).

S5 Computational Cost of G-BAGs
In this section, we show that G-BAGs have computational complexity of order n at each iteration of the MCMC sampler in Section S4. For explanatory purposes, we assume |S i | = for H i|h and R i|h , and R −1 i|h for density. Since C [i|h] is of size 2m × 2m and R i|h is of size m × m, their inversion leads to complexity O(9m 3 ).
Repeating for each i and h, the overall complexity regarding (a) becomes O(9M Km 3 ). Similarly, (b) requires C −1 [u|h] for H u|h and R u|h , and R −1 u|h for density. The Jm × Jm matrix C [u|h] is common for all u ∈ U i , and R u|h is a 1 × 1 matrix. Thus, the total floating point

S6.1 Fitted G-BAG is correctly specified
To fit a G-BAG model, a vague normal prior N (0, 100) is specified for β. IG priors are specified for nugget (τ 2 ) and spatial variability (σ 2 ) such that τ 2 ∼ IG(2, 0.1) and σ 2 ∼ IG(2, 1). For instance, a ∼ U nif (7.330, 14.667) corresponds to a temporal range between roughly 0.5 and 1 at which the correlation drops to 0.1. The prior distribution for probabilities to choose one of the four directions is given as π i ∼ Dir(α) with α = 1/4 for i = 1, . . . , 288. In Q-MGP, its default priors (σ 2 , τ 2 ∼ IG(2, 1), β ∼ N (0, 100)) are used except for the covariance decay parameters which have the common uniform prior whose lower bound is the minimum of lower bounds for a and c in the G-BAG model, and an upper bound is the maximum of upper bounds for a and c. Penalized Complexity (PC) priors  are used for SPDE model parameters. For AR(1), the autocorrelation parameter is assumed to be larger than 0.05 with prior probability 0.99, and the precision of white noise is larger than 1 with prior probability 0.01. Under stationarity, the spatial standard deviation is larger than 1.5 with prior probability 0.01, and the range parameter is smaller than 1.8 or 9 when the truth is θ 1 or θ 2 , respectively. In a nonstationary case, ψ's are modeled with i.i.d . N (0, 3) priors.

S7 Air Quality Analysis in South Korea
South Korea offers locational benefits to study flow of air quality because it is embedded in a region where air quality exhibits rapid changes in time and space and strongly impacted by transboundary pollution, dust, or smoke from China (Crawford et al., 2021). We analyze hourly PM2.5 levels in Korea between 8pm on February 4, 2020 and 7am on February 5, 2020.
This time of a day is chosen to minimize traffic effects that appear to significantly influence local air quality. AirKorea data provided by the Korean Ministry Of Environment ( N (0, τ 2 ) is assumed with the mean centered log(PM2.5) as y(t). We assume τ 2 ∼ IG(2, 0.1), σ 2 ∼ IG(2, 1) a ∼ U nif (0.052, 38), c ∼ U nif (0.001, 0.081), and κ ∼ U nif (0, 1) as prior distributions. The following four directions (SW, W, NW, N) are placed in a bag because these have been identified as the prevailing winds over Korea (Heo et al., 2010), also aligning with the directions of long-range movement from China. Analysis results are reported based on a posterior sample of size 1,000 obtained from 10,000 iterations after 8,000 burn-ins and saving every second sample in the subsequent 2,000 samples. The total run time with 10 threads is reported in minutes in Table S2.
In predicted surfaces of Figure S3, it is visible that the PM2.5 level is gradually cleared over night from west and northwest directions, which is well captured by the posterior mean direction of winds. This exemplifies that the cause and removal of air pollution is strongly affected by wind directions. Among four directions in the bag, W and N tend to be selected with higher posterior probabilities (less uncertainty) over Korea. Table S2 shows that G-BAG outperforms Q-MGP in terms of lowering prediction errors and quantifying prediction uncertainty. This suggests that allowing for flexible selection of DAGs can improve prediction over a fixed DAG when directional association is evident or believed to exist in space and time.  Figure S3: Observed (left column) and predicted (middle) log(PM2.5) by G-BAG in South Korea. On the right, the posterior mean direction of winds is overlaid on discretized predicted results for PM2.5 based on its standards by Korean MOE.

S8 Air Quality Analysis in CA
Same priors are used as in Section S7 except a, c ∼ U nif (0, 1000) and β ∼ N (0, 100). The G-BAG and Q-MGP models have the same priors for all parameters. The location-specific spatial range is assumed in the SPDE model as ρ(t) ∝ exp(ψ 1 + ψ 2 t 1 + ψ 3 t 2 ) where t 1 and t 2 are scaled eastings and northings in [0, 1], respectively. PC priors are given identically as in Section 3.1 except independent N (0, 1) priors for ψ 1 , ψ 2 , and ψ 3 . Comparing residuals of SPDE-nonstationary and those of G-BAG, subtracting fitted w from y, the former exhibit apparent spatial patterns unlike the latter ( Figure S4). These patterns are rarely explained by x(t)β either because the fitted β by SPDE-nonstationary (-0.038) is small compared to the residuals ranging from -6 to 6. We observe that the spatiotemporal random effects from SPDE-nonstationary tend to underestimate PM2.5 when the true PM2.5 is high and overestimate when it is low.   Figure S6: The 95% posterior predictive CI widths on November 13, 2020. The surfaces bear much resemblance at any time point.