Multivariate spatial data consist of multiple variables that exhibit inherent associations among themselves as well as underlying spatial associations unique to each of them. While modeling each variable separately captures its spatial distribution independent of other variables, the resulting data analysis fails to account for associations among the variables, which can substantially impair prediction or interpolation [see, e.g., 9, 11, 21, 50]. Joint modeling of multivariate spatial data proceeds from vector-valued latent spatial stochastic processes or random fields, such as a multivariate Gaussian process. These are specified with matrix-valued cross-covariance functions [see, e.g., 24, 43, 33, and references therein] that model pairwise associations at distinct locations. While theoretical characterizations of cross-covariance functions are well established, modeling implications and practicability depend upon the specific application [see, e.g. 5, 34, 46, 32, 20, 44].

Multivariate spatial processes are customarily specified using a vector of mean functions and a cross-covariance matrix function. Let z ( s ) = ( z 1 ( s ) , … , z q ( s ) ) T be a q × 1 stochastic process, where each z i ( s ) is a real-valued random variable at location s ∈ D ⊆ ℜ d . The process is specified by its means E [ z i ( s ) ] = μ i ( s ) and the second-order covariance functions C i j ( s , s ′ ) = Cov { z i ( s ) , z j ( s ′ ) } for s , s ′ ∈ ℜ d and i , j = 1 , 2 , … , q. These covariances define the matrix-valued q × q cross-covariance function C ( s , s ′ ) = { C i j ( s , s ′ ) } with ( i , j )-th entry C i j ( s , s ′ ). Since the mean structure can be estimated by absorbing it into a regression component, i.e., z i ( s ) = μ i ( s ) + w i ( s ), where w i ( s ) has zero mean, we focus on constructing a valid cross-covariance function for a zero-mean latent process. From its definition, C ( s , s ′ ) need not be symmetric, but must satisfy C ( s , s ′ ) T = C ( s ′ , s ). Also, since var { ∑ i n a i T z ( s i ) } > 0 for any finite set of distinct locations s 1 , s 2 , … , s n ∈ D and any set of nonzero constant vectors a 1 , a 2 , … , a n ∈ ℜ q ∖ { 0 }, we have ∑ i = 1 n ∑ j = 1 n a i T C ( s i , s j ) a i > 0. See [24] for a review.

Spatial factor models [35, 40, 47, 53] are among the most widely used approaches to build multivariate spatial models when the number of spatially dependent variables is large. These models build upon the popular linear model of coregionalization (LMC) [26, 44], which specifies the multivariate random field as a linear combination of r univariate random fields. Choosing r ≪ q produces low-rank or spatial factor models. Computational benefits ensue; one needs specify only a small number r of univariate processes to ensure non-negative definiteness and, hence, yields large q. While computationally convenient, multivariate spatial analysis using LMC prohibits interpretation of spatial structures of each variable. For example, LMC endows each z j ( s ) with the same smoothness (the smoothness of the roughest latent process). This is implausible in most applications because different spatial variables typically exhibit very different degrees of smoothness. An alternate approach constructs cross-covariances by convolving univariate processes with kernel functions [48, 49, 36]. However, barring certain very special cases, the resulting cross-covariance functions are analytically intractable, hence less interpretable, and may require cumbersome numerical integration for estimating process parameters. Some of the aforementioned difficulties are obviated by a conditional approach developed in [12], where the univariate GPs are specified sequentially, conditional on the previous GPs assuming some ordering of the q variables. Other notable approaches include [1] who considered using latent dimensions to embed all the variables in a larger dimensional space and use standard covariance functions on this augmented space. However, this embedding restricts all pairwise-correlations among the different variables to be positive. [25] and [2] directly formulated multivariate Matérn cross-covariance functions where both the univariate covariance functions for each variable and the cross-covariance functions between each pair of variables are members of the Matérn family. The multivariate Matérn GP is appealing in terms of interpretability, allowing direct spatial inference for each variable via estimates of the parameters of the corresponding component Matérn GP.

While such methods have been applied to diverse data sets, they have been restricted to a moderate number of outcomes. This article focuses on the high-dimensional multivariate setting where a significantly larger number of outcomes (tens to hundreds of variables) can be measured at each spatial location. Such settings are becoming increasingly commonplace in the environmental and physical sciences where inference is sought on large numbers of dependent outcomes. [17] developed a novel class of “Graphical Gaussian Process” (GGP) models for scalable and interpretable analysis of high-dimensional multivariate spatial data. The underlying idea is to exploit conditional independence among the variables using an undirected graph with the components of the multivariate GP as the nodes. Absence of edges between nodes represent conditional independence among the corresponding pair of component processes given all the other nodes.

Graphical models are extensively used to represent the joint distribution of several variables for many types of non-spatial data. Applications in spatial settings have been concerned with reducing the computational burden for large n by replacing the complete graph between locations with nearest-neighbor graphs [14]. Current multivariate GP covariance functions do not lend themselves naturally to incorporating an inter-variable graphical model. [30] proposed a method for parsimonious joint analysis of such high-dimensional multivariate spatial data via a common basis function expansion and imposing sparse graphical models for estimating the covariances of the coefficient vectors. This approach, akin to similar approaches common in multivariate functional modeling [52, 54], rely on multiple replicate measurements of each variable at each location which may not be available in many multivariate spatial applications.

The approach of Dey et al. [17] constructs a high-dimensional multivariate model by adapting graphical Gaussian models to process-based settings and does not require replicate data. This process-level graphical model approach addresses some key properties of multivariate GPs that are deemed critical for handling highly multivariate data, including the retention of the flexibility and interpretation of spatial and non-spatial dependencies. The balance of the manuscript proceeds as follows. We offer a brief overview of inference from graphical models in spatial statistics followed by an elucidation of the GGP and its implementation in a fully Bayesian modeling framework. We illustrate with some simulation examples, offer an analysis for a highly multivariate environmental data, and conclude with some discussion and pointers for future research.

Undirected graphical models have an established history in spatial statistics in the context of specifying spatial dependencies via Markov random fields (MRF) for regionally aggregated or areal data [see, e.g., 6, 7, 10, 41, 22, 5, and references therein]. In areal modeling [37, 23, 28, 27] the nodes of the graph represent regions and the presence of an edge between two nodes indicates that the two regions are neighbors (or adjacent). Spatial models are constructed from this graph by assuming dependence between nodes with edges between them. For example, the popular conditional auto-regression (CAR) models assume conditional independence between two nonadjacent nodes given all other nodes.

Our current development departs from the multivariate areal setting in two aspects. First, we consider point-referenced settings where inference is desired at the process-level. Unlike areal settings where inference is limited to the fixed set of areas, point-referenced datasets allow inference about the variables over every conceivable location in an entire region. Second, unlike MRFs where the graph is on the set of areas specifying spatial dependencies, here the undirected graph posits inter-variable conditional independence relationships. Thus each node corresponds to the stochastic process describing the distribution of the respective variable over space.

In point-referenced settings, graphical models among locations based on spatial proximity have been used to specify nearest neighbor Gaussian Processes (NNGP) [14, 18] that offer a scalable solution to use of GP priors in Bayesian hierarchical models. However, akin to the MRF models for areal data, NNGP uses a graph among the locations to parsimoniously specify spatial dependence. Neither approaches model inter-variable relationships using graphs.

Recent developments for graphical modeling of multivariate functional or spatial data represents the multiple variable processes using a common univariate basis function expansion with multivariate (vector-valued) coefficients [52, 30]. Graphical modeling among the variables is induced by a series of graphical models, one for each of the vector-valued coefficients. As practically, basis functions are truncated to a finite number of coefficients, the representation reduces to the standard graphical model estimation for vector-valued data using graphical Lasso type techniques [19]. The advantage of the approach is that the graph is allowed to vary for each coefficient and thereby allowing conditional independence structures to be resolution-specific. However, these methods rely on replicate data on each variable at each location which is atypical in many spatial settings. Also, the assumption of a common univariate basis function expansion for the multivariate process may be inadequate.

We first present a simple alternative approach to build a multivariate spatial model that explicitly respects a given inter-variable graphical model. Let G = { V , E } be an undirected graph with V being the set of vertices and E the set of edges. A customary approach for developing probabilistic models on such undirected graphs is to specify full conditional distributions for each node given others and then deriving a joint density from the set of full conditional distributions. This is achieved using Brook’s Lemma [6, 37], which can be adapted to our setting as follows.

Let us consider a multivariate spatial process z ( s ) with q univariate spatial processes and let z i = ( z i ( s i 1 ) , z i ( s i 2 ) , … , z i n i ( s i n i ) ) T be an n i × 1 random vector corresponding to the realizations of the i-th variable, z i ( s ), for each i = 1 , 2 , … , q. We specify the following sequence of full conditional distributions for modeling z i ’s given all other z j ’s: (2.1) z i | z ( − i ) ∼ N ∑ j = 1 q A i j z j , Γ i , i = 1 , 2 , … , q , where z − ( i ) = { z 1 , z 2 , … , z q } ∖ { z i }, i.e., the collection of z j ’s for j = 1 , 2 , … , q but excluding z i , A i j ’s are fixed n i × n j matrices, A i i = O (the matrix of zeroes), and Γ i ’s are fixed positive definite matrices. Since each variable can exhibit its own spatial dependence structure, the Γ i varies by variable. Brook’s Lemma provides a straightforward method for deriving the joint density from (2.1) using the identity: (2.2) π ( z 1 , z 2 , … , z q ) = ∏ i = 1 q π ( z i | z ˜ 1 , … , z ˜ i − 1 , z i + 1 , … , z q ) π ( z ˜ i | z ˜ 1 , … , z ˜ i − 1 , z i + 1 , … , z q ) ) × π ( z ˜ 1 , z ˜ 2 , … , z ˜ q ) , where z ˜ = ( z ˜ 1 T , z ˜ 2 T , … , z ˜ q T ) T is any fixed point in the support of π ( z ) and we assume that the joint density π ( · ) > 0 over its entire support. The proof is a straightforward verification proceeding from the last element of the right hand side (i.e., the joint density on the right hand side). Note that

Rather than building a Gaussian graphical model whose nodes are finite-dimensional random vectors, GGP builds a graphical model whose nodes are the component Gaussian processes of a multivariate spatial GP and edges represent process-level conditional dependence between the incident nodes given all other nodes. To do so, we first clarify what it means for two spatial process to be conditionally independent. Let the q nodes of G represent q different spatial processes, { z i ( s ) : s ∈ D } for i = 1 , 2 , … , q. Two distinct processes z i ( · ) and z j ( · ) are conditionally independent given the remaining q − 2 processes, which we denote by z i ( · ) ⊥ z j ( · ) | z − ( i j ) ( · ), if the covariance between z i ( s ) and z j ( s ′ ) is zero for any pair of locations s , s ′ ∈ D conditional on full realizations of all of the remaining processes. A q-variate process z ( s ) is a GGP with respect to an inter-variable graph G = { V , E } if z i ( · ) and z j ( · ) are conditionally independent given the remaining processes whenever ( i , j ) ∉ E.

Any collection of q independent spatial processes is a GGP with respect to any graph G but is of limited use as it does not leverage inter-variable dependencies. Sparse graphical models are often used as parsimonious working models serving as a scalable approximation to a richer (or fuller) data generating model. Hence, a more relevant question is given a reasonable inter-variable graph G, possibly inferred from scientific knowledge, how to construct a GGP that respects G and provides best approximation to a multivariate GP with a given cross-covariance function. Thus, we wish to construct a multivariate spatial process such that it conforms to the posited graphical relationships among the q dependent variables and best approximates the spatial structures of a given q × q cross-covariance matrix C ( s , s ′ ) = ( C i j ( s , s ′ ) ). [17] proved that such an optimal GGP exists and is unique, and is one that preserves the marginal spatial covariance structure for each variable i as specified by C i i ( s , s ′ ) as well as cross-covariances C i j ( s , s ′ ) between any pairs of variables ( i , j ) ∈ E.

In this regard, we recall a seminal paper by Dempster [16] on covariance selection. The key result can be described as follows. Let G = { V , E } be any undirected graph with q nodes and let F = ( f i j ) be any positive definite covariance matrix whose elements f i j give the covariance between random variables at nodes i and j. Then, the best approximation of F (in terms of the Kullback-Liebler distance) among the class of covariance matrices that satisfy the conditional independence relationships specified by the Gaussian graphical model G is a unique positive definite matrix F ˜ = ( f ˜ i j ) such that: 1.

f ˜ i i = f i i ;

Covariance selection, as described above, is applicable to finite-dimensional inference, where, for example, we restrict attention to parameter estimation and random effects at a specified set of n spatial locations. Let C = ( C ( s i , s j ) ) be an n q × n q spatial covariance matrix, where each C ( s i , s j ) is a q × q submatrix evaluated using a valid cross-covariance function. Given an undirected graph that posits conditional independence relations among the q variables, we can compute the unique C ˜ preserving marginal variances and covariances and conforming to the graph using an iterative proportional scaling (IPS) algorithm [45, 51].

Consider the graph in Figure 1 specifying the conditional independence structure among 4 variables. We consider each variable being observed at 10 locations uniformly sampled from a ( 0 , 1 ) × ( 0 , 1 ) grid. We assume a multivariate Maten covariance structure [2] between the processes and compute the 40 × 40 cross-covariance matrix C. Now, Algorithm 1 lays out the steps to obtain the unique cross-covariance matrix C ˜ corresponding to the GGP that conforms to the graph in Figure 1. The resulting precision matrix is plotted in Figure 2 which has zero entries only in the ( 1 , 3 ) and ( 2 , 4 )-th block, i.e., when there are no edges present between the variable nodes.

However, this finite-dimensional application of covariance selection or IPS does not immediately extend to the process-level formulation of conditional dependence in a multivariate GP. We now illustrate how GGP of [17] extends covariance selection to infinite-dimensional framework that will allow predictive inference at arbitrary spatial locations over the entire domain of interest D.

We begin the construction with a given inter-variable graph G, a multivariate cross-covariance function C ( s , s ′ ) and with a finite set of locations L. Using covariance selection, we obtain a covariance matrix C ˜ : = C ˜ ( L ) such that we can define a multivariate spatial model on L as (3.1) w ˜ ( L ) ∼ N ( 0 , C ˜ ) The properties of covariance selection ensure that w ˜ i ( L ) ∼ N ( 0 , C i i ( L ) ), C o v ( w ˜ i ( L ) , w ˜ j ( L ) ) = C i j ( L ) for all ( i . j ) ∈ E, and C o v ( w ˜ i ( L ) , w ˜ j ( L ) | w − i j ( L ) = ( C ˜ − 1 ) i j = O for all ( i , j ) ∉ E. Thus on L, each component of w ˜ retains marginal covariances as specified by C, marginal cross-covariances for variable pairs included in the graph are also preserved, and conditional dependencies align with the graph.

To extend this finite-dimensional framework on L into a process-level graphical model for a multivariate GP on D, we leverage the property that a univariate GP w i ( s ) with covariance function C i i ( s ) can be constructed as the sum of two orthogonal GPs w i ∗ ( s ) and r i ( s ). Here w i ∗ ( s ) is a finite rank predictive process that can be represented as w i ∗ ( s ) = b i ( s ) T w ˜ i ( L ) where w ˜ i ( L ) is some random vector such that w ˜ i ( L ) = d w i ( L ) and b i ( s ) are such that the predictive process w i ∗ ( · ) agrees with w i ( · ) on L [3]. The second GP r i ( · ) is often referred to as the residual GP as it explains the residual covariance between realizations of the GP w i ( · ) beyond what can be explained by the predictive process. Choosing w ˜ i ( L ) to be the components of w ˜ ( L ) from (3.1) and r i ( · )’s to be independent across the variables i, one creates a multivariate GP w ( s ) = B ( s ) w ˜ ( L ) + r ( s ) where B ( s ) = ⊕ b i ( s ) and r ( s ) = ( r 1 ( s ) , … , r q ( s ) ) T .

Covariance selection and predictive process harmonize to ensure that the resulting process w ( · ) exactly or approximately satisfies all the three criteria for being an optimal multivariate process given a covariance function C and a graphical model G. The component processes satisfy the law w i ( s ) ∼ G P ( 0 , C i i ) thereby preserving the marginals. Cross-covariances for ( i , j ) ∈ E are also preserved exactly on L and approximately elsewhere when choosing L to be sufficiently representative of the domain. Finally, the conditional independence of w ˜ ( L ) with respect to G induced by covariance selection together with the component-wise independent residual processes r i ( · ) ensure that the resulting multivariate process w ( · ) conforms to process-level conditional independences as specified by G and is thus a GGP. This construction is referred to as stitching – the multiple processes can be envisioned as multiple layers of fabric which are connected to each other only at the locations L via the edges of G, serving as threads.

We demonstrate the practical use of GGP for modeling non-stationary spatial time-series data using daily PM2.5 measurements from monitoring stations in 11 northeastern states of the contiguous US, including Washington DC, for February 2020. The data are publicly available from the United States Environmental Protection Agency (EPA) website. Our current analysis focuses on n = 99 monitoring stations that recorded at least 20 days of data in both 2019 and 2020. From the National Center for Environmental Prediction’s (NCEP) North American Regional Reanalysis (NARR) database, we obtained daily values of meteorological factors (temperature, barometric pressure, wind-speed, and relative humidity) that are posited to influence PM2.5 levels. We used multilevel B-spline smoothing using the MBA package in R to impute daily values of these variables and merge them with the EPA data. To account for baseline levels and weekly periodicity, we included a 7-day moving average of the PM2.5 data for each station and day in 2020, centered around the same day of 2019, and subtracted the day-of-the-week specific means from the raw PM2.5 values. We then analyzed the resulting highly multivariate (28-dimensional) spatial data set over n = 99 locations.

After conducting exploratory analysis that revealed autocorrelations among pollutant processes on consecutive days (for both lag 1 and 2), even after adjusting for covariates, we employed a graphical Matérn GP with an A R ( 2 ) graph (Figure 12). The marginal parameters for day t are σ t t , ϕ t t , and τ t 2 , the autoregressive cross-covariances are denoted by b t − 1 , t and b t − 2 , t between day t and days t − 1 and t − 2, respectively. GGP enables us to model non-separability in auto-covariances across both space and time, as well as time-varying marginal spatial parameters and autoregressive coefficients.

The GGP requires only 53 cross-covariance parameters (27 parameters for lag 1 and 26 parameters for lag 2). As the largest clique size in an AR(2) graph is 3, the largest matrix one must deal with for the data at 99 stations is only 297 × 297. We present the estimates and credible intervals for lag 1 and lag 2 auto-correlation parameters r t − 1 , t and r t − 2 , t (normalized b t − 1 , t and b t − 2 , t ) obtained from GGP in Figures 13 and 14, respectively. These estimates display substantial variation across time, with the many spikes indicating high positive autocorrelation. Specifically, 95 % Bayesian credible intervals for 6 out of the 27 ( 22 %) of the lag 1 estimates and 4 out of the 26 ( 15 %) of the lag 2 estimates exclude 0, providing strong evidence in support of non-stationary autocorrelation across time. The presence of significant lag 2 autocorrelations justifies our choice of A R ( 2 ) graph. The Graphical Matérn GP also yields impressive predictive performance ( 93 % coverage) on hold-out data (Figure 15).

This article has developed and elucidated, with examples, the construction of multivariate Gaussian processes with associations among variables modeled by a valid spatial cross-covariance function while conditional independence between variables conforming to a posited undirected graph. The sparsity of the posited graph accrues substantial computational gains and makes the proposed approach especially attractive for datasets with very large numbers of spatially dependent outcomes. The algorithms implemented here have been especially designed to exploit the structure of decomposable graphs. Our examples demonstrate the algorithmic efficiency of chromatic Gibbs samplers used to update the latent process and the cross-covariance parameters and also the inferential efficiency, in terms of estimation and prediction, of the graphical Gaussian process model. Future avenues for research will include incorporating scalable spatial processes, such as the Nearest-Neighbor Gaussian Process [15, 13] for spatial and spatial-temporal processes and extensive comparisons with competing methods for multivariate outcomes, such as factor models, from computational and inferential standpoints.