{'fraction_non_alphanumeric': 0.026982593515364146, 'fraction_numerical': 0.009610006610057457, 'mean_word_length': 4.645584154626351, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A nonstationary spatial covariance model for data on graphs

Spatial data can exhibit dependence structures more complicated than can be represented using models that rely on the traditional assumptions of stationarity and isotropy. Several statistical methods have been developed to relax these assumptions. One in particular, the"spatial deformation approach"defines a transformation from the geographic space in which data are observed, to a latent space in which stationarity and isotropy are assumed to hold. Taking inspiration from this class of models, we develop a new model for spatially dependent data observed on graphs. Our method implies an embedding of the graph into Euclidean space wherein the covariance can be modeled using traditional covariance functions such as those from the Mat\'{e}rn family. This is done via a class of graph metrics compatible with such covariance functions. By estimating the edge weights which underlie these metrics, we can recover the"intrinsic distance"between nodes of a graph. We compare our model to existing methods for spatially dependent graph data, primarily conditional autoregressive (CAR) models and their variants and illustrate the advantages our approach has over traditional methods. We fit our model and competitors to bird abundance data for several species in North Carolina. We find that our model fits the data best, and provides insight into the interaction between species-specific spatial distributions and geography.


Introduction
Models for spatially indexed data traditionally assume that observations close to one another in space are more highly correlated than observations which are far apart. Spatial models that assume the correlation between any two observations can be represented by a function of only the distance in space between those observations are said to be stationary, in which case the spatial autocorrelation is the same in all parts of the spatial domain, and isotropic, spatial autocorrelation is the same in all directions [Cressie, 1993, Banerjee et al., 2003. While these assumptions allow for a simple and interpretable mathematical representation of spatial processes, such assumptions are often not appropriate for real world data [Guttorp et al., 1994].
As an example of this, refer to Figure 1. This plot was taken from the website of the eBird database, a project of the Cornell Lab of Ornithology that crowd-sources bird population data by allowing anyone to record and submit their birdwatching observations to the database. This is done in the form of user-submitted, location-indexed checklists containing the counts of each species observed by the bird watcher [Sullivan et al., 2009]. Figure 1 depicts the percentage of submitted checklists for each 20 by 20 kilometer gridded region within eastern North Carolina that contain an observation of at least one brown pelican, a seabird common along the coast of North Carolina. Dark purple rectangles correspond to regions where the bird is most commonly observed. It is readily apparent that the spatial distribution of the brown pelican throughout this region is highly dependent on the proximity to the coast, and as such, a model for brown pelican abundance which assumed that the spatial dependence between observations is the same for all locations in North Carolina and in all directions, would be woefully inadequate at representing the data. Consider for example how the frequency of observation along the northeast-southwest axis of the domain is relatively constant when compared to frequency along the northwest-southeast axis. Additionally, when similar data for other bird species within the same spatial domain are considered, one may find nonstationary dependence exhibiting a completely different spatial pattern. However, it is reasonable to imagine that for any given species (or more generally, for any nonstationary and anisotropic spatial process) there is a unique distance metric underlying the space in which data are observed, in which the covariance patterns are isotropic and stationary with respect to this metric rather than to the geographic distances between observed location. Using Figure 1 as an example, one could say that for a brown pelican, two locations 50 kilometers apart along the coast might be viewed as closer together than two locations with one on the coast and the other 50 kilometers apart inland, perpendicular to the coastline.
Many classes of spatial models have been developed for use with nonstationary and anisotropic data [Schmidt and Guttorp, 2020]. Among these methods is the "spatial defor-mation approach," which parameterizes spatial covariance in terms of an estimated metric space corresponding to a deformation of the Euclidean space in which data are observed [Sampson and Guttorp, 1992] and which leads to interpretation similar to the "intrinsic distances" described in the brown pelican example from above. The vast majority of these methods have been developed specifically for point referenced data. However, areal data is also common within the world of spatial statistics, either because some random variables are observed only at a regional level, or because raw point-referenced data might be aggregated and processed into areal or gridded regions out of convenience or due to scientific or policy interest in questions at the county, state, or national level. Such data are representable using graphs, with a node representing each region and edges connecting each pair of adjacent regions.
Conditional autoregressive (CAR) models, first introduced by Besag [1974], are commonly used for modeling spatially dependent data observed on a graph, but are potentially limited when faced with more complex spatial dependence patterns. Given a random vector y with y 1 , ..., y p being observations at each of the p nodes (or vertices) of some graph G = (V, E), where V contains the set of nodes, and E is set of edges between those nodes, a CAR model specifies a Markov random field in which each element of y has a conditional distribution determined by the values at neighboring nodes, for example where j ∼ i indicates that node j is adjacent to node i, W is a p × p adjacency or weights matrix (where w ij > 0 if j ∼ i else w ij = 0), κ controls the degree of spatial correlation and σ 2 is a variance parameter. This results in a marginal distribution for y of y ∼ N (0, σ 2 (I n − κW) −1 ).
The CAR model may be thought of as the spatial generalization of the first order autoregressive (AR1) model often used in time series modeling. There are many varieties of CAR model that have been developed for use in generalized mixed models [Besag et al., 1991], with multivariate data [Gelfand and Vounatsou, 2003], and which can utilize weighting schemes incorporating geographic and covariate information Hooten, 2013, Ejigu andWencheko, 2020]. Most of the literature regarding spatially dependent graph data is limited to various forms of autoregressive models [Ver Hoef et al., 2018a]. The majority of CAR models and their variants tend to to treat all first-order neighbors within a graph as equally proximate due to the use of graph adjacency structure alone their definition of spatial dependence. This will rarely result in effective representation of nonstationary spatial processes. In instances where more complex weighting schemes have been utilized, such as in Hanks and Hooten [2013], and Ejigu and Wencheko [2020], edge weights are defined as linear functions of a small set of environmental covariates. While such approaches allow for more flexibility than those using adjacency alone, the space of possible weighting matrices under those models is restricted by the column space defined by the collection of environmental covariates being used. There are also instances in which covariates associated with graph edges may not be available to researchers, in which case a model must be defined using only the structure of the graph itself. One computational advantage of the CAR modeling framework is that the induced precision matrix will generally be sparse due to the assumption of independence for each observation conditional on the values of its neighbors. This assumption can come with trade-offs in terms of predictive efficiency; Gramacy and Apley [2015], Katzfuss and Guinness [2021], Guinness [2018] all demonstrate the shortcomings of models which condition only on local or nearest neighbor structure (as is the case in CAR models) in comparison to those which also condition on more distant observations. For more discussion on CAR models see Ver Hoef et al. [2018b] and Lee [2011] for review papers which provide detail on the use, development, and interpretation of these models and their variants. To address the concerns regarding limitations of CAR models for use with nonstationary spatial processes, we propose a model for the spatial covariance of graph data using a different set of assumptions, which we will describe in the following sections. We will also illustrate that our proposed method has greater capacity to represent certain correlation structures that may be of interest to statisticians and ecologists, especially in contrast to the most commonly used CAR variants.
One approach to modeling spatial nonstationarity and anisotropy for point referenced data is the "spatial deformation" approach introduced by Sampson and Guttorp [1992]. Unlike CAR models and most other spatial models which characterize covariance as a function of a known and fixed orientation of the data within some space (whether discrete as with areal data or continuous with point referenced data), the spatial deformation approach assumes instead that covariance is determined by the distances between the data within some latent space which is estimated. Given some spatial process {Y (s) : s ∈ G} over a space G (in most applications G is R 2 and we refer to it as the "geographic-" or "G"-space), which is repeatedly observed at some finite set of locations s 1 , s 2 , ..., s n ∈ G at times t ∈ 1, ..., T , we can model the covariance between any two locations in G as where σ 2 is a scale parameter, ρ(·) is a valid correlation function, and δ(·) is a non-linear transformation mapping G to a latent space (referred to as the "deformation-" or "D"space) in which the stationary and isotropy assumptions hold. Covariance is thus a direct function of the distances between points in the D-space, which allows for nonstationarity in the covariance when viewed from the G-space. The challenge with the implementation of this approach lies in determining an appropriate framework to estimate the function δ(·).
One should define δ(·) such that the orientation of the points in the G-space is considered within the function, (proximity in the G-space should generally imply relative proximity in the D-space), and in a manner that avoids producing a D-space on which it is impossible to obtain a positive definite covariance matrix. In their 1992 article, Sampson and Guttorp estimate δ(·) using a two-step process in which they first estimate the latent space locations δ(s 1 ), ..., δ(s n ) via a nonmetric multi-dimensional scaling algorithm that is initialized using the geographic locations such that {δ(s i )} (0) = {s i } before converging to a set of locations in the D-space such that a stress criterion based on the variogram of the data within the D-space is minimized; they then use thin plate splines to obtain the full function δ(·) which maps from the geographic locations to the latent space locations over the entirety of the Gspace. Later articles further developing this method propose using a Gaussian-process prior for δ : R 2 → R 2 centered at the observed locations {s i } [Schmidt and O'Hagan, 2003]. The spatial deformation approach allows for considerable flexibility in the covariance structure of the data, and provides an intuitive interpretation of model results if one views the distances between observations in the D-space as the "intrinsic distances" between locations in the G-space.
In this article we present a covariance model for graph data that has similar interpretation to the deformation approach of Sampson and Guttorp. Fitting the model involves the estimation of a distance metric over our graph which in turn corresponds to an embedding of our graph into a high-dimensional Euclidean space. Covariance is then assumed to be a function of the distances between nodes in the embedding, similarly to how the spatial deformation approach assumes stationarity in the D-space, resulting in potentially nonstationary, anisotropic covariance structure in the G-space.
In the next section of this article we establish some of the necessary background information regarding covariance functions and graph metrics. In particular, we highlight a class of graph metrics which are compatible with most traditional covariance functions such as those from the Matérn family [Matérn, 1960]. In Section 3 we formally introduce our model, and provide some detail on its interpretation as well as its implementation. In Section 4 we proceed to demonstrate the model's performance in comparison to traditional autoregressive models through a case study using bird abundance data for several species in the state of North Carolina, and show that our model performs best in terms of WAIC among the models considered for the majority of species evaluated. As part of this case study, we also demonstrate the manner in which our model can provide valuable insight into the interaction between a spatial process and the underlying geography via the posterior distance matrix of the graph on which a process is observed. We conclude this article with a discussion of our model's strengths and limitations in comparison to existing methods, as well as potential extensions for future work.

Background
Suppose we are given a graph G = (V, E) and random variables {Y v : v ∈ V } observed at the nodes of G. In this article we propose a model for Σ vv = Cov(Y v , Y v ) which is described visually in Figure 2. The covariance matrix Σ is parameterized by W, a matrix of unknown positive edge weights, and σ 2 > 0, a scale parameter. The parameter W is an element of W G , the space of all possible edge weight matrices given G, and W is mapped to a distance matrix D via a graph metric ∆(·); D is in turn mapped to a symmetric positive definite matrix Σ using a covariance function C(·) which also takes the parameter σ 2 > 0 as input. In the following two subsections, we establish the background information necessary to justify our model and present it formally in Section 3. Figure 2: Visual representation of the model presented in this article. The edge weight parameter W is mapped to a distance matrix, which in conjunction with the scale parameter σ 2 is mapped to the covariance matrix Σ.

Covariance Functions and Distance
Let D p be the space of metric distance matrices between p points. That is, if D ∈ D p then D is a p × p matrix where all diagonal elements of D are zero, the off-diagonal entries are positive, D is symmetric, and all elements of D satisfy the triangle inequality (i.e. for any . In our model we define the p × p covariance matrix Σ as is a covariance function which maps distance matrices to S + p , the space of p × p symmetric positive definite matrices. To allow for flexibility in the covariance, we will treat the distance matrix of a graph as an unknown parameter to be estimated from the data. This in contrast to most spatial models, which take the distances between observations as known and fixed. Instead we will parameterize distance as a function of edge weights on the graph, which are estimated from the data. This approach allows us define covariance structures for graph data that are considerably more complex than ones obtained using the neighborhood structure of the graph alone.
Covariance functions defined only by the distances between locations, and not the locations themselves, are said to be stationary and isotropic. However, as shown by the deformation approach of Sampson and Guttorp [1992], such functions can be used to represent nonstationary and anisotropic spatial dependencies if the distances used as inputs are allowed to deviate from the fixed distances defined by the geographic space in which data are observed. Stationary and isotropic covariance functions are often defined such that they can be applied element-wise to a distance matrix. For D ∈ D p we can also define Σ by where c(·) is some function applied to a distance between two points. For Σ to be a valid covariance matrix it must be symmetric and positive definite, and thus c(·) must be chosen carefully with respect to D. There are many well known covariance functions that can be used to produce positive definite matrices using distances as inputs. Functions from the Matérn covariance family, originally developed by Matérn [1960] and given by are among the most commonly utilized within spatial statistics [Cressie, 1993, Banerjee et al., 2003]. Here the input d is a distance between two points, σ 2 is a scale parameter, τ is a range parameter, ν is a smoothness parameter, and K ν is a modified Bessel function of the second kind. The Matérn family includes the commonly used exponential and Gaussian (also known as the radial basis kernel or squared exponential) covariance functions [Banerjee et al., 2003].
If c(·) is a member of the Matérn family of covariance functions it is necessary for the the distances contained in D to be Euclidean in order to ensure that Σ is a positive definite matrix [Ver Hoef, 2018]. Here we define the space of Euclidean distance matrices D E p ⊂ D p such that if D ∈ D E p , there exists {x 1 , ..., x p } ⊂ R k for some k such that d ij = x i − x j . In our model, for a given graph G = (V, E) with p nodes we will require that D, the matrix containing pairwise distances between the nodes of G, satisfies D ∈ D E p . This also implies that there exists an embedding of G into a Euclidean space such that the graph distances are equal to the distances between nodes in the embedding.
The fact that many of the most common covariance functions (including those from the Matérn family) require Euclidean distances as inputs to guarantee positive definiteness is reasonably well known, and can be confirmed when one considers Bochner's theorem, which provides a necessary and sufficient condition for functions on Euclidean vectors to be positive definite (see Gihman and Skorohod [1974]). However, this requirement has at times been disregarded by researchers who were interested in utilizing common spatial modeling tools

Graph Metrics
To parameterize the covariance of random variables observed on a graph by applying the Matérn covariance function to a matrix of estimated distances, we must be careful that D is parameterized in a manner that ensures that it is Euclidean over the entirety of the parameter space. To do this, we must consider what graph metrics produce Euclidean distances and how their properties will influence model interpretation.
Given an undirected, connected, simple graph G = (V, E) with p nodes, we define a class of distance matrices parameterized by the edge weights matrix W of G. Many graph metrics are functions of W where w ij > 0 if node i is adjacent to node j and w ij = 0 otherwise [Jungnickel, 2012]. We define W G to be the space of all edge weights matrices possible given G, and note that all information regarding the edges of G is contained in each W ∈ W G based on whether or not an element of W is equal to zero. The distance matrix D for a graph G is given by D = ∆(W) where ∆(·) is the function for a particular metric which maps an edge weight matrix to a distance matrix.
For a given graph G and metric ∆(·) we define D ∆ G as the space of distance matrices which are possible to obtain under different edge weight configurations for G, that is We have established that for Σ = C(D) to be a valid covariance matrix with C(·) a Matérn family covariance function, D must be a Euclidean distance matrix. Thus, in choosing a graph metric ∆(·), we must ensure that D ∆ G ⊂ D E p . In doing so, we can specify a parameterization of Σ in terms of an unobserved edge weights matrix W ∈ W G with the property that Σ = C(∆(W)) ∈ S + p for all W ∈ W G . This also implies that because D = ∆(W) ⊂ D E p , W characterizes an embedding of G into Euclidean space. This embedding is conceptually similar to the mapping from the G-space to the D-space in the deformation approach. In Section 3 we establish that for appropriate choice of metric, our covariance model is identifiable, which implies that the embedding of G implied by W is unique up to isometry.
The most commonly used graph metric is the shortest path distance [Chebotarev, 2011] (denoted here as ∆ SP (·)) which defines the distance between nodes i and j as where P ij is the set of all possible paths between nodes i and j, but can intuitively be described as the sum of the weights along the "shortest" path connecting those nodes. However, this metric generally produces non-Euclidean distance matrices [Klein and Zhu, 1998] and as such Σ = C(D) is not guaranteed to be positive definite if C(·) is Matérn.
A commonly used alternative graph metric is the resistance distance [Chebotarev, 2011], denoted ∆ Ω (·), where the distance between nodes i and j is defined as where L + is the Moore-Penrose generalized inverse of the graph Laplacian, defined as L = diag(W1 p ) − W where diag(·) is a function taking a p-length vector and returning a p × p diagonal matrix, 1 p is a vector of ones, and {e i } 1:p are the standard basis vectors. This metric was introduced to the mathematical literature by Klein and Randić [1993] with roots in electrical physics and is based on the formula for calculating the effective resistance between nodes in a network of resistors.
Attractive properties of resistance distance include the fact that the distance between two nodes is reduced by the number of paths connecting said nodes, a property that shortest path distance does not possess [Klein and Randić, 1993], and that it is proportional to commute time distance, i.e. the expected hitting time of a random walk between two nodes [Chandra et al., 1996]. Use of resistance distance is widespread within ecological settings as both of the above properties characterize a type of flow-based graph connectivity which corresponds nicely with existing models for the movement of individual organisms as well as animal populations [Thiele et al., 2018, Peterson et al., 2019.
Because L + is always a positive semi-definite matrix, resistance distance may be viewed as a squared Euclidean (or Mahalanobis) distance. Klein and Zhu [1998] highlight that a related class of metrics of the form where m > 0 is an exponent (i.e {L + } 2 = L + L + ), are guaranteed to produce Euclidean distance matrices for any m > 0 and all W ∈ W G . Thus, all metrics of this form can be used in conjunction with any covariance function which requires Euclidean distances as inputs, such as those from the Matérn class. The square root of resistance distance is one of these metrics (when m = 1), but we are particularly interested in using the second metric of this form (when m = 2)-often referred to as the quasi-Euclidean metric-because it shares many of the same properties regarding graph connectivity with resistance distance and has been demonstrated to be well behaved in comparison to other metrics, especially when used with planar graphs which are typically seen in spatial and ecological settings Klein, 1996, Ivanciuc et al., 2001]. The quasi-Euclidean metric also scales linearly with the inverse of the edge weights, that is, cD = ∆ (2) (W/c) for c > 0. This property is convenient for computation and enables greater parameter interpretability. Subsection 3.2 contains a numeric illustration of how the quasi-Euclidean metric operates within our model.

Covariance Model
Given a simple, undirected graph G = (V, E) with p nodes, we make repeated observations These observations may be stored in the n×p data matrix Y, where y ij is the ith observation at the jth node of G. We define Σ, the between-node covariance of each y i , in terms of the intrinsic distances between the nodes of G which are in turn a function of the known structure of the graph (V, E) and unobserved edge weight parameter matrix W, which is estimated from the data. Our parameterization of Σ could be used in conjunction with any number of distributional assumptions for Y, including temporal dependence or non-normality, but to make presentation concrete we initially define our model using the assumption that {y i } 1:n are independent draws from a N (0 p , Σ) distribution. We will discuss possible extensions later in the article. We parameterize Σ as follows: Under this model, the data Y comes from a mean-zero matrix-normal distribution with column covariance Σ and row covariance I n . The parameter W is a p × p edge weight matrix, σ 2 is a scale parameter, and ρ ν (·) is the Matérn correlation function given below: Common practice is to take the smoothness parameter ν as fixed rather than as a parameter to be estimated. Within this article we will set ν = 5/2, which implies a spatial process that is smoother than those given by an exponential covariance function, but less smooth than those defined with Gaussian covariance; however, any value of ν > 0 would be valid in this model [Stein, 1999].
Note that the correlation function defined above does not include an explicit range or spatial-decay parameter. Equation 4, the standard form of the Matérn covariance function, contains the parameter τ whereas the correlation function we are using, given in Equation 9, does not. Despite the lack of a parameter τ , which directly scales the distances used as inputs in the correlation function, our parameterization of Σ implicitly controls for the rate of spatial decay thru the edge weight matrix W: Let D be the p × p distance matrix produced by applying the quasi-Euclidean metric to the edge weight matrix W. The third and fourth lines of Equation 8 show how D is obtained. Let ∆ (2) (·) be the function such that D = ∆ (2) (W). Recall that cD = ∆ (2) (W/c) for any positive constant c. The distances in our model are directly scaled by the magnitude of the edge weights, which are themselves unrestricted above zero rendering an additional parameter τ where correlation is a function of d/τ unnecessary. To interpret our model's estimate of the spatial decay in the covariance of Y, we recommend re-scaling the distance matrix D after fitting. If we define D s = D/max(D) and τ s = 1/max(D), then the covariance defined by Σ jk = ρ ν (d s jk /τ s ) is equal to the one defined Σ jk = ρ ν (d jk ), and τ s has equivalent interpretation to a range parameter in a covariance model for point-referenced data in which the geographic distances between observations were scaled to have a maximum distance of one.
As parameterized, our model is identifiable: Proposition 1 (Identifiability) Using the parameterization for Σ given in Equation 8, A proof of the model's identifiability is provided in the appendix of this article.

Interpretation
Rather than characterizing nonstationary and anisotropic spatial covariance via the use of a complex covariance function or kernel, our model uses a relatively simple covariance function, and pushes most of the modeling complexity into the estimation of latent distances between graph nodes, which are themselves functions of a set of edge weight parameters; nodes that are close to one another are more correlated, while nodes that are far apart are less correlated. The metric we utilize to define distances over the graph was chosen to ensure that the resulting covariance matrix will be positive definite for any combination of edge weight parameters, but what of the edge weight parameters themselves? Generally speaking, larger edge weight parameters indicate greater correlation between nodes in the regions of the graph where that edge is found, while smaller edge weight parameters indicate lower correlation in those regions. When edge weights approach zero, the corresponding edges are effectively omitted from the graph, which may be thought of as a type of automatic model selection for the design of the graph itself. Figure 3 contains an illustration of how edge weights impact covariance in our model.
It depicts four sets of edge weights for the same five-node graph and shows the resulting correlation matrix when a Matérn correlation function with ν = 5/2 is applied to the distance matrix produced by the edge weights. Subfigures 3a and 3b highlight how the scale of edge weights matter, with larger edge weights producing greater between-node correlation.
Subfigure 3c demonstrates how it is possible to obtain nonstationarity when certain edge weights are larger than others, with nodes along the path from 4 to 1 to 2 exhibiting greater intercorrelation than the rest of the graph. Subfigure 3d shows that it is possible under our model to obtain covariance structures where non-adjacent nodes (1 and 5) are more correlated (in effect closer together) than any other pair of first-order neighbors, which is possible due to the properties of the quasi-Euclidean metric, which characterizes the distances between nodes as shorter the greater the number of connecting paths between them.

Comparison to CAR models
The covariance model specified here is quite flexible, especially in contrast to the commonly used graph-based methods which treat all first-order neighbors as equally proximate, such as most CAR models. We illustrate this with a simple example. The correlation component of our model for Σ is determined fully by W ∈ W G , while σ 2 simply acts to scale the correlation matrix defined by W. Applying the quasi-Euclidean metric to W in turn characterizes D ∆ G , the space of distance matrices possible under the model, which in turn can be used to define S (C•∆) G ⊂ S + p , the subset of covariance matrices that can be obtained under our model for a given graph, metric, and covariance function. Fitting our model to data results in us finding the "best"Σ in S (C•∆) to approximate Σ as defined by maximizing likelihood, minimizing Bayesian loss, or some other approach. Just as with our model, other specifications of Σ using fewer than p(p − 1)/2 parameters (the number of free elements in a p × p covariance matrix) characterize a proper subset of S + p containing all covariance matrices that are possible to obtain under that model. Figure 4 depicts a graph G with five nodes, and two correlation matrices, Σ 1 and Σ 2 , which represent possible dependence structures for data observed on G; we note that neither . The matrix Σ 1 depicts a correlation structure in which one side of G (the path connecting nodes 4-1-2) exhibits greater correlation than the rest of the graph, while Σ 2 depicts a correlation structure in which non-adjacent nodes (1 and 5) are more correlated with one another than any pair of adjacent nodes in G. These two configurations of Σ were chosen to show types of graph nonstationarity which could be present in real world data, but are not meant to represent the full range of possible cases. In this example, we assume a multivariate normal, mean-zero distribution and find that best approximations of Σ 1 and Σ 2 from S (C•∆) G minimize Kullback-Leibler (KL) divergence to a greater degree than the equivalent best approximations from the CAR models we considered.
Specifically,Σ 1 andΣ 2 are obtained by identifying the optimal set of edge weights such that The most widely used CAR variant is the first-order CAR model, which we will refer to as the CAR1 model; it incorporates only the basic adjacency structure of G into its design and has a marginal distribution which defines where A is the fixed adjacency matrix of G. Let S CAR1 G be the space of covariance matrices possible to obtain under the CAR1 model. Given that the CAR1 model has only two parameters, it may be unfair to compare it to a more complex model such as the one proposed in this article, which has number of parameters equal to the number of edges in G (8 in the case of G from Figure 4) plus one. We thus propose a second more complex CAR model, which we refer to as the weighted CAR or CARw model, which represents the maximally flexible correlation structure possible under the CAR framework, (as characterized by this article,) in which we treat all non-zero elements of the weights matrix as parameters, defining Let S CARw G be the space of covariance matrices possible to obtain under this CAR model, which is specified using a weights matrix W for which edge weights are estimated individually, just as with the model described in Section 3. This type of autoregressive model is uncommon within the literature, though a similar specification appears in Smith et al. [2015]. In general, the limited number of existing CAR approaches utilizing a flexible weights matrix tend to define the elements of W as some function of environmental covariates, geographic distances, or the sizes of areal units rather than estimating individual edge weights directly [Hanks and Hooten, 2013, Wang et al., 2016, Ver Hoef et al., 2018b.
For each of the three models discussed above and the two correlation matrices from This was done by generating 100,000 samples from N 5 (0 5 , Σ) and obtaining a Bayes estimator for Σ using the simulated data (the details for obtaining this estimator are described in the following subsection). The Bayes estimator has been shown to concentrate about the parameter value which minimizes KL-divergence between a model and the true distribution as sample size increases [Zellner, 1988, Gelman et al., 2014. For this exercise we scaled each draw of Σ within our sampler to ensure that it was a correlation matrix; this was done to improve the fairness of the comparison between our model and the CAR models which will generally have heterogeneous diagonal variance if specified as in equations 11 and 12. KL-divergence indicates that the subset of covariance matrices defined by a particular model gets closer to the true data-generating distribution. As can be seen, our model is able to more closely approximate both Σ 1 and Σ 2 than the CAR variants we considered. The CAR1 model performs notably worse than our method, which is unsurprising given that it uses far fewer parameters and the KL-divergence does not account for model complexity. The CARw model has one more parameter than our model, yet still performs worse for both Σ 1 and Σ 2 with the gap in performance being more apparent for Σ 2 . This is due in part to the fact that our model has greater flexibility than most autoregressive models to produce covariance matrices for which non-adjacent nodes are more correlated with one another than to any of the intermediate nodes due how the quasi-Euclidean metric defines distance in the presence of multiple viable paths between locations.

Priors and Model Fitting
While one could obtain parameter estimates using maximum likelihood estimatation, the complexity of the relationship between edge weights and the covariance matrix makes max-  Table 1: KL(N (0,Σ i ), N (0, Σ i )) for each model and correlation matrix. The KL-divergence of the best performing model for each "true" correlation matrix is in bold.
imizing the likelihood challenging in comparison to Bayes estimation using a Markov chain Monte Carlo (MCMC) approach. We suggest priors in accordance with the parameter space of our model which requires that σ 2 and all edge weight parameters {w jk : j ∼ k} must be positive real numbers. If the nodes of the graph correspond to a known orientation in physical space (as is frequently the case when dealing with areal data) one can obtain the geographic distance matrix D Geo by calculating the physical distances between centroids of areal units. It is then possible to encode this geographic information into the prior for the edge weights, which may be helpful for model fitting, especially when dealing with a limited number of observations of the process over the graph. In such cases, we propose using the following priors: Here d Geo jk is the geographic distance between nodes j and k, and the prior expectation of the associated edge weight is a w /(d Geo jk b w ) which is smaller for adjacent nodes which are far apart than for adjacent nodes which are closer in physical space. In instances where no geographic information is available about the graph on which data is observed, it is appropriate to simply use i.i.d. gamma priors for the edge weights. The choice of inverse-gamma prior for σ 2 is natural when Y has a normal likelihood due to built-in conjugacy. The choice

Algorithm 1 MCMC Sampler for Covariance Model
Input: Y an n × p data matrix, and G a graph with p nodes Output: T posterior samples for model parameters σ 2 and W Initialize σ 2(0) and W (0) and set Σ (0) = c ν (σ 2(0) , W (0) ) for t = 1 to T do Update σ 2(t) by taking draw from full conditional π(σ 2(t) |Y, Generate u from U (0, 1) if u < min(1, α) then w cur jk ← w * jk and set w cur jk = w cur kj for all k ∈ K Σ cur ← Σ * W (t) ← W cur Σ (t) ← Σ cur to use independent gamma priors for all edge weights is somewhat simplistic, as it may be reasonable to assume some level of dependence between the weights of incident edges. The MCMC sampler we use to fit this model is provided in Algorithm 1, which describes how to obtain posterior samples for σ 2 and edge weight matrix W. In essence, we use the full conditional inverse-gamma distribution of σ 2 to perform a Gibbs update at each iteration, and perform a block Metropolis-Hastings update for the edge weights, in which new values for groups of edges are proposed and either accepted or rejected.
The model is complex and the sampler used is somewhat computationally expensive, but we have found that it mixes well and within reasonable time periods for most small and medium sized graphs (p < 500 nodes). The data for the application described in the following section had n = 36 observations and p = 100 nodes. It required approximately three hours to obtain 2500 iterations of the sampler using code written in R and run single threaded on the lead author's personal laptop equipped with an Intel Core i7-9750H processor. The code used for our analysis is included with the supplemental materials.

An Application Using NC eBird Data
Returning to the data discussed in the introduction of this article, and which motivated the development of this method, we now fit our model to bird abundance data within the state of North Carolina. As alluded to in our introduction, there are varying patterns of between-site dependence among species within this spatial domain, and we are interested in obtaining the species-specific intrinsic distances which characterize that spatial dependence. The data were retrieved from the eBird database, a project of the Cornell Lab of Ornithology that allows bird watchers across the globe to record and submit their bird watching observations in the form of checklists that contain the specimen count for each observed species as well as the observer's location and "effort," or time spent bird watching [Sullivan et al., 2009].
The model was fit to the data for several common species within the state of North Carolina during the time period from January 2018 to December 2020. The species considered were the northern cardinal, the Carolina wren, the mourning dove, the turkey vulture, the Canada goose, the laughing gull, and the mallard. In addition to being abundant throughout the state, these species were chosen as they represent a considerable range of genetic and habitat diversity.
For each species in the data set, observed counts were aggregated temporally by month and spatially by county. The data was stored in a 36 × 100 matrix X, in which x ij represents the combined number of individual birds from that species observed during month i = 1, ..., 36 within the borders of county j = 1, ..., 100 by all birdwatchers contributing to the database. In addition to the species-specific counts matrices, there is a 36 × 100 effort matrix T common to all species, where t ij contains the total time spent bird watching during month i within county j, as well as a geographic distance matrix D Geo , obtained by calculating the physical distances between the geographical centroid of all bird watching locations within each county. To account for non-uniformity in T across counties, we transform the observed counts to obtain Y, the 36×100 matrix of log-rates at which species were observed. We define y ij = log([x ij + .5]/t ij ). After subtracting the row means from Y to account for temporal variation in bird abundance, we find that the residuals given by Y − {Y1 100 1 100 /100} have histograms that are approximately bell-shaped for all seven species, leading us to the choice to model Y as coming from a normal distribution.
By treating the 100 counties of North Carolina as nodes, with 256 edges connecting each pair of adjacent counties, we obtain the graph representation of the state depicted in Figure   5.
For each species, we model Y using the following spatial random effects model with temporally-varying mean to account for seasonal changes in bird abundance: Here µ is an intercept term and α i is the temporally-varying mean effect for month i. The spatial random effects matrix U has independent rows and between-county covariance given us is Σ and the underlying distance matrix D, which can be viewed as containing the speciesspecific intrinsic distances between counties in North Carolina, with other model parameters primarily serving to enable an appropriate fit for Y. Note that the model as written does not include any environmental covariates, although it would be straightforward to include them in the mean of the linear model. It could also be of interest to include environmental covariates in the estimation of Σ, a point which we will discuss in Section 5.
We use the priors given in equation 13 for σ 2 and W and use the following priors for the remaining model parameters: To ensure identifiability of the mean effects we fix α 1 = 0. In our implementation of this model we set prior parameters a σ 2 = a w = a ψ = a α = a β = 2 and b σ 2 = b w = b ψ = b α = b β = 1. Because these priors result in fully conjugate conditional distributions, this model can be fit using the sampler provided in Algorithm 1 with the addition of Gibbs steps for updating the parameters in equation 15.

Interpretation of Results
To interpret the output of our model (which we subsequently refer to as the nonstationary graph covariance or NSGC model), it is useful to have a cursory understanding of the geography of North Carolina. The majority of physical geographic variability within the state of North Carolina tends to run on an east to west axis. Figure 6 depicts the four main geographic regions of the state, with the easternmost tidewater region, followed by the inner coastal plain, the Piedmont plateau in the center of the state, with the Blue Ridge mountains running through the westernmost part of North Carolina [NCpedia, 2015, Swick, 2016.
We are interested in how these geographic regions may manifest themselves within the intrinsic distances across the state for each species. Are there species for which distances along the eastern coast of the state are greater or smaller relative to the intrinsic distances from the coast moving inland? Are there species for which the mountain ranges in the west of the state significantly impact their spatial distribution while other species are unaffected?  and four line segments labeled "A", "B", "C", and "D", with these segments chosen to reflect various physical regions within the state, or the borders between them. For each of the seven bird species, there is a different distance associated with each of the four line segments; these distances are shorter, and the counties at their endpoints are "closer" to one another if there is high correlation in the bird counts between locations. Subfigures 7b and 7c show the posterior distributions of the intrinsic distances associated with two pairings of those segments. Both of these pairings corresponds to a geographic contrast that may be useful for understanding how these species spatial distributions interact with the physical properties of the environments they inhabit. We see that the intrinsic distances associated with each of the segments differ from species to species, with the posterior for the laughing gull being more divergent from the other six species in both subfigures. This divergence makes sense when one considers that the laughing gull it is the only seabird of the seven species, having a natural range of habitats with the least overlap among the species considered [Swick, 2016]. We see this in subfigure 7c, which shows that for a laughing gull, the distance associated with line segment "D," which runs along the coast, is relatively shorter when compared to the other species, while the distance associated with "C," which crosses the border from the tidewater region inland to the coastal plains, is relatively longer.
This can be interpreted to mean that locations along the coast are "closer" together for this species than the other six, and indicates that proximity to the ocean influences the betweencounty covariance for laughing gulls to a greater degree than for the other species. Subfigure 7b highlights the contrast between the distances associated with segment "A," which spans across a mountain range, and "B" which runs along the relatively flat Piedmont. Birds with posteriors concentrated at higher values along the x-axis of this plot, such as the Carolina wren, may be more inhibited in movement by the changes in elevation and habitat along segment "A" than the others.
We would also like to visualize the entire mean posterior distance matrix for each species.
Because this matrix corresponds to a high dimensional Euclidean embedding of the graph, it is challenging to represent with a single plot and without significant distortion. We have found that plots of the posterior mean edge weights are not straightforward to interpret due to issues of scale and dependence on the number of nodes within a given neighborhood.
Instead we instead present a visualization which highlights the contrast between the posterior mean intrinsic distance and the physical distance between each pair of adjacent counties; this provides us not only a sense of the relative distances between counties, but more significantly of the deformation of the graph's geographic orientation necessary to produce the intrinsic distances. We have done this by calculating z d jk , a score designed to represent whether the intrinsic distance between adjacent nodes j and k is smaller or greater than suggested by the physical distance between those nodes.
Here D(E) indicates the set of distances between nodes connected by edges of our graph, scales a set to have mean zero and standard deviation one. High values of z d jk indicate that nodes j and k are intrinsically "closer" together relative to their physical separation, while low values indicate that they are more distant. Figure 8 contains two plots depicting z d jk for all adjacent counties. Subfigure 8a contains a map of North Carolina for the Carolina wren, a species present throughout the state, but is less common in colder regions and at higher altitudes [Swick, 2016]. We note the band of dark blue line segments running approximately along longitude −82 • E, which indicates that observations of Carolina wren on either side of this border are less correlated than one would expect if they believed the data generating process to be stationary. Significantly, that border tracks along a similar line between the mountainous and plateau regions of the state, revealing a physical explanation for what we have observed from our model's output.
Subfigure 8b depicts a similar map for the laughing gull, which is found most abundantly along the coast of the state, and becomes progressively less common as one moves inland [Swick, 2016]. The edges along the coast are red and orange, indicating greater connectivity, with a blue band running parallel to the coast at approximately −77 • E, along the border between the coast and coastal plains regions of the state. This may be viewed as a sort of buffer-zone which inhibits movement or between site dependence from the coast inland.
There is another teal band along −79 • E, where the coastal plain meets the Piedmont plateau, which may correspond to a second buffer zone which inhibits further westward movement. It is worth noting that laughing gulls west of that longitude are extremely rare within the data.
For both species, these plots reveal significant interactions between their spatial distributions and the geography of where they were observed, despite the fact that no environmental covariates were included in this version the model, a point on which we will elaborate in the following section. Taken together, the plots (along with the plots not included for the other five species analyzed) reveal that there is a unique covariance structure and underlying distance between locations associated with each species.
An additional quantity of interest in many spatial data applications is the ratio σ 2 /(σ 2 +ψ) which may be thought of as the percentage of total variability attributable to the spatial process (in comparison to the nugget effect represented by ψ).

Fitted Model Comparison
We are interested in evaluating the performance of our model on this dataset in comparison to existing methods for similar datasets, many of which have already been discussed within this article. In addition to the NSGC model, we fit the following four models for Σ to the data for each of the seven bird species. Other than the specification of Σ, each model was fit in the same manner according to equations 14 and 15.
Naive: Σ = σ 2 ρ 5/2 (D Geo /τ 2 ) The notation here is as before, with A 4 denoting the adjacency matrix for G containing first thru fourth-order neighbors-that is to say a 4jk = 1 if there is a path with four or fewer edges connecting nodes j and k and a 4jk = 0 otherwise-and indicating the Hadamard product. The notation exp(−D G /τ 2 ) indicates that the exponential function is applied element-wise to the matrix −D G /τ 2 . The naive model applies a Matérn covariance to the geographic distances between county centroids, treating the data as though it were pointreferenced rather than areal. The CAR1 model is the basic CAR model described in section 3.3, while the SAR1 model is another common autoregressive model that tends to produce smoother covariances than the CAR1 [Ver Hoef et al., 2018a]. The CAR4Dr model is the best performing model from Ver Hoef et al.
[2018b], a review paper which discusses and compares sixteen different autoregressive models for spatial data on graphs. It utilizes a higher order neighborhood structure, which is in turn weighted by the geographic distances between areal regions. The CAR4Dr model also uses row-standardization for its weight matrix, a common practice in the autoregressive model literature [Ver Hoef et al., 2018a].
In fitting these four models the same priors given in equation 15 were used for all parameters not involved in the construction of Σ. Of the remaining parameters, and using the notation from equation 17 each instance of σ 2 was given a diffuse inverse-gamma prior, τ 2 was given a diffuse gamma prior, and κ a uniform prior on (-1,1). As before, the full conditional of σ 2 results in a Gibbs update, while the other two parameters can be sampled using a Metropolis-Hastings step.
Upon obtaining posterior samples from each model for each species, we can analyze them to assess model performance. The Watanabe-Akaike information criteria (WAIC, also referred to at times as the "widely applicable information criteria") was introduced by Watanabe and Opper [2010], and is designed to approximate the results of leave-one-out cross validation. It also tends to penalize model complexity to a greater degree than other common information criteria such as AIC and DIC. Gelman et al. [2014] discusses WAIC and other information criteria within a Bayesian context, indicating a preference for WAIC, finding that it is reasonably consistent with the results given by cross-validation while being considerably cheaper to compute. Table 3 contains the WAICs for each of the five models when fit to each of the seven bird species we analyzed. As can be seen, our method (NSGC) performs best by this criteria for five of the seven species considered despite considerably greater model complexity, which suggests that the spatial dependence present in the data for several species is more complex than can be represented using simpler methods. The two birds for which our model did not perform best according to WAIC were the mourning dove and the mallard. In the case of the Mourning Dove, the NSGC model had the highest log likelihood, but the model complexity penalty of WAIC lowered it into second place, indicating a relatively simple dependence structure for that species. Our model performed the worst of the five for the mallard data, and even had the lowest log-likelihood among the five models. This is likely due to the relative lack of spatial structure in the Mallard data as discussed in Section 4.2. We note as well that the range of WAIC values across all five models for the Mallard data is extremely narrow in comparison to the ranges present for other species; this makes sense considering that the choice of how to model spatial random effects is irrelevant in cases where no meaningful spatial effect is present.

Discussion
We now briefly discuss a few potential applications, extensions, and questions related to the ideas presented within this article. One of the most significant pieces of output from our model is the posterior distance matrix. By itself, it can be analyzed to better understand the relationships between nodes of a graph in terms of their latent proximity to one another; but as was the case in our application with the eBird data, there are instances in which one may fit this model repeatedly for each of several variables (unique bird species in the case of the eBird example) which are observed on the same graph. In such cases, an analysis of how the posterior distance matrices for each individual variable compare to one another  circumstances under which one may have a data set and a spatial graph structure on which it was observed, but no covariate information or intuition regarding an appropriate set of covariates to include in a model. Our method is able to account for the spatial dependence that may be attributable to those missing covariates, and a posteriori analysis of the distance matrix may enable one to learn about the relationship between covariates and covariance.
That said, in many instances we will be interested in the inclusion of environmental covariates in our model. We may wish to incorporate them into the mean via a linear model, and allow Σ to represent the covariance between observations after accounting for environmental covariates. Alternatively, we may be interested in incorporating environmental covariates directly into the specification of Σ itself. There are many instances of researchers including environmental covariates in their spatial covariance model as a method of representing covariate-based nonstationarity as discussed in Schmidt and Guttorp [2020]. Such approaches have also been used in CAR models as in Hanks and Hooten [2013] and Wang et al. [2016]. A reanalysis of the eBird data using an adaptation of our method which accounts for the relationship between edge weight parameters and environmental covariates could be fascinating as a tool understand what factors most influence spatial connectivity for each species.
Within this article we limited ourselves to cases with normally distributed data. There is however a wealth of literature for spatially dependant non-Gaussian models. The generalized CAR framework of Besag et al. [1991] could be straightforwardly adapted for use with our covariance model. Species observation data is often zero inflated; accounting for this as in Agarwal et al. [2002] could prove beneficial to model quality in applied settings.
Other potential extensions include joint species modeling using our methodology. One approach could be to perform factor analysis [Lawley andMaxwell, 1962, Harman, 1976] on the model parameters characterizing species-specific covariance enabling one to better understand which species interact with the spatial domain similarly in terms of covariance and intrinsic distances, while providing a natural approach to modelling the intrinsic distances for each species hierarchically.
The method described in this article provides a flexible and novel framework for modeling the covariance of data observed on a graph structure. Its development was inspired by the spatial deformation approach of Sampson and Guttorp [1992] and the posterior distance matrix produced by fitting our model may be interpreted as containing the intrinsic distances between the nodes of a graph. We have demonstrated that our model is equipped to represent nonstationary and anisotropic spatial dependence, and that there are contexts in which it performs better than advanced versions of autoregressive models which are the standard tools for modeling spatially dependent graph data. We also find the specification of between-node covariance in terms of intrinsic distances a useful framing for understanding and interpreting model results, as exemplified by our discussion of model output in section 4.2 in which we found that an analysis of the posterior distance matrix revealed connections between the covariance in the eBird dataset and the physical geography of the spatial domain on which it was observed.
We now need only to prove that W 1 = W 2 ⇒ D 1 = D 2 . The transformation from edge weights matrix to to distance matrix as defined by the quasi-Euclidean metric can be rewritten as follows: Here d {L + } 2 is defined to be the vector with elements equal to the diagonal entries of {L + } 2 and (·) • 1 2 denotes the Hadamard (element-wise) square root of a matrix. The Hadamard root is an invertible transformation (B 1 = B 2 ⇔ B Let A p be the space of p × p symmetric matrices such that ∀A ∈ A p , A1 p = 0 p . We note that for any graph with p nodes, both its Laplacian matrix L and {L + } 2 are elements of A p . It can be seen that N (∆) ∩ A p = {0 p×p }, which implies that ∀A ∈ A p , ∆(A) = 0 p×p ⇒ A = 0 p×p . We now show that this condition implies the injectivity of ∆ over A p .