Modeling spatial data using local likelihood estimation and a Matérn to SAR translation

Let $f(\mbox{\boldmath$x$})$ be a mean zero Gaussian process on $\mbox{\boldmath$x$}\in{\mathbb{R}}^{2}$ with covariance function $k(\mbox{\boldmath$x$},\mbox{\boldmath$x$}^{\prime})\in{\mathbb{R}}$. The Matérn family of stationary covariance models is important because of its flexibility and the interpretability of its parameters. The Matérn covariance function with a unit range parameter is

$$C(\,d\,|\,\nu,\sigma^{2}\,)=\sigma^{2}{\frac{2^{1-\nu}}{\Gamma(\nu)}}{(d)}^{\nu}\mathcal{K}_{\nu}{(d)},$$

where $d$ is the Euclidean distance between $x$ and $\mbox{\boldmath$x$}^{\prime}$, $\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind of order $\nu$, and $\Gamma(\cdot)$ is the gamma function. $\sigma^{2}$ is the spatial process variance, and $\nu$ is the smoothness parameter which controls the mean square differentiability of the process. The isotropic covariance function with range parameter $\kappa$ is given by

$$k(\mbox{\boldmath$x$},\mbox{\boldmath$x$}^{\prime})=C(\kappa d\,|\,\nu,\sigma^{2}).$$

In contrast to modeling a covariance function for a process of continuous spatial variation, the SAR model parameterizes the precision matrix for the process on a discrete lattice. For the following development we denote by $\mathbf{y}$ a Gaussian process on an infinite regular lattice in ${\mathbb{R}}^{2}$.

Denote by $y_{i,j}$ the element of $y$ at lattice location $(i,j)$. A simple isotropic SAR model can be written using lattice notation as

$$\begin{array}[]{c|c|c}0&-1&0\\ \hline\cr-1&4+\kappa_{S}^{2}&-1\\ \hline\cr 0&-1&0\\ \end{array}$$

(1)

which visually illustrates a set of decorrelating weights on the random vector $y$. To be precise, we interpret (1) as implying that the following equation holds

$$(4+\kappa_{S}^{2})y_{ij}-(y_{i-1,j}+y_{i+1,j}+y_{i,j-1}+y_{i,j+1})=e_{ij},$$

(2)

for a mean zero unit variance normally distributed white noise vector $e$ with components registered to the spatial grid $\{e_{ij}\}_{i,j}$.

Here $\kappa_{S}>0$ is suggestive of a range parameter controlling the correlation length scale dependence of the field and is similar, but not identical, to $\kappa$ for the continuous Matérn case above. For the model in (1), one can populate a matrix $B$ using (2) such that $B\mbox{\boldmath$y$}=\mbox{\boldmath$e$}$. Furthermore, it is clear that if $\kappa_{S}>0$, the diagonal dominance of $B$ guarantees its invertibility. With $\mbox{\boldmath$y$}=B^{-1}\mbox{\boldmath$e$}$, the covariance matrix for $y$ is $B^{-1}B^{-T}$ and thus $y$ has precision matrix $Q=B^{T}B$. Moreover, it is straightforward to see that the precision matrix implied by the SAR model is sparse. The sparsity property makes the SAR model amenable to modeling large data sets because the sparse precision matrix can be used in place of a dense covariance matrix for likelihood estimation and simulation. Following the ideas from [20] one can iterate the spatial autoregressive weights to obtain higher order models that approximate smoother processes. For example, if $BB\mbox{\boldmath$y$}=\mbox{\boldmath$e$}$, this implies a SAR model extending to second-order nearest neighbors and gives the precision matrix: $Q_{2}=(BB)^{T}(BB)$. The SAR model detailed here is a special case of a GMRF. For a given row of the precision matrix, the nonzero, off-diagonal entries index the neighbors that determine the Markov property. Citing the order of neighborhoods can cause some confusion depending on whether one is referring to SAR weights , or precision, matrix. For the first-order SAR described above, the nonzero elements in $Q$ will include second-order neighbors. A first-order SAR will be a GMRF based on second-order neighbors and the weights will depend on $B$. Two additional points concerning the SAR model are important. First, for a given dataset $y$ is not defined on an infinite lattice, but rather a (typically) rectangular one; this stencil of (1) should be modified at the boundaries of the domain. The center value of the stencil should be the $\kappa_{S}^{2}$ minus the sum of the weights of its non-zero neighbors, although a simpler option is to artificially extend the lattice a few nodes outwards to reduce boundary effects. Second, the value of $\kappa_{S}$ is one-to-one functionm for the marginal variance of the process.

A version of the SAR model that exhibits approximate stationarity and also geometric anisotropy will be detailed in Section 3.

Lindgren and Rue [20] developed an approximation of Gaussian random fields with Matérn covariance functions using GMRFs with particular SAR structures. The connection is established through a stochastic partial differential equation (SPDE) formulation in that a Gaussian field $u(\mathbf{s})$ with stationary Matérn covariance is a solution to the SPDE

$$(\kappa^{2}-\Delta)^{\nicefrac{{\alpha}}{{2}}}u(\mathbf{s})=\mathcal{W}(\mathbf{s})$$

where $\alpha=\nu+\frac{d}{2}$, $\kappa>0,\nu>0$, $\mathbf{s}\in\mathbb{R}^{d}$, $d=1$ or $2$, and $\mathcal{W}(s)$ is a generalized white noise process with zero mean and variance $\sigma^{2}$. As in the Matérn model, $\nu$ controls the smoothness of realizations of the Gaussian field. Fixing $\nu=1$ and $d=2$, [20] showed that the SAR covariance structure obtained by discretizing the pseudodifferential operator $(\kappa^{2}-\Delta)$ approximates a Matérn covariance structure with range $\kappa\approx\kappa_{S}$; this relationship is approximate. Similar results can be obtained for different smoothness parameters $\nu$ by convolving the finite difference stencil in (1) with itself $\nu$ times, as detailed in the previous section for $\nu=1$ and $\nu=2$.

The connection between the isotropic Matern family and a SAR relies on the approximation of a discretized Laplacian operator with finite differences of the fields on a lattice. To develop an accurate statistical model, it is important to quantify the error in such an approximation and improve its calibration over the limiting expression suggested in [20]. In this section we provide numerical evidence to show that an accurate calibration is possible if restricted to specific ranges of the covariance parameters.

Our calibration setup is as follows: given a Matérn range parameter $\kappa$, we estimate the value of $\kappa_{S}$ in the SAR model which gives the best approximation to the Matérn correlation function. We consider the smoothness of the Matérn model fixed at $\nu=1$ and $\nu=2$, and with unit marginal variance for all models. The first step is to fix the Matérn range parameter and evaluate a Matérn correlation matrix for the process evaluated on the lattice grid. Then, we find the optimal $\kappa_{S}$ from the SAR model by computing its equivalent correlation matrix (the standardized inverse precision matrix) and minimizing distance between the Matérn and SAR correlation matrices.

It is known that the SAR covariance model suffers from edge effects. To avoid the interference of edge effects, we quantify the difference between the two correlation matrices by only comparing the correlations from the central lattice point under both models. For an $N\times N$ lattice of locations, with $N$ odd, let $\boldsymbol{\sigma}_{\kappa}$ denote the vector of correlations between the center point in this lattice and all other locations based on the Matérn covariance function with range $\kappa$. Let $\boldsymbol{\sigma}_{\kappa S}$ be the analogous correlation vector for the SAR model with range parameter $\kappa_{S}$. We then find

$$\hat{\kappa}_{S}(\kappa)=\underset{\kappa_{S}}{\operatorname{argmin}}\|\boldsymbol{\sigma}_{\kappa}-\boldsymbol{\sigma}_{\kappa S}\|_{2}.$$

If the relationship proposed by [20] holds we would expect $\hat{\kappa}_{S}(\kappa)\approx\kappa$.

$\kappa^{-1}$ is varied over the interval $[1,20]$ with $N=51$ and lattice points having unit spacing. The approximation results are summarized in Figure 1. In 1(a) $\hat{\kappa}_{S}(\kappa)^{-1}$ is plotted as a function of $\kappa^{-1}$. Orange corresponds to the $\nu=1$ case, cyan corresponds to $\nu=2$, and the solid black line shows the theoretical relationship, $\kappa_{S}^{-1}=\kappa^{-1}$ from [20]. From this experiment we conclude that, at least at this level of discretization, it is important not to rely on the analytic formula to translate between $\kappa$ and $\kappa_{S}$ parameters.

The relative error of using the SAR correlation with $\kappa_{S}$ value derived from the numerical experiment is shown in 1(b). The $\ell_{2}$ distance measure used in the optimization of the model correlation matrices is used to quantify the resulting model error, normalized by $\|\boldsymbol{\sigma}_{\kappa S}\|$. The relative error incurred when using $\frac{1}{\kappa_{S}}=\frac{1}{\kappa}$ as shown by the dashed lines can be on the order of $20\%$ or worse for higher correlation ranges. However, using the calibrated range (solid lines) gives relative errors on the order of $5\%$.

The Matérn family can be extended to include geometric anisotropy by defining a distance measure based on a linear scaling and rotation of the coordinates. Let $A=D^{-1}U^{T}$ be a $2\times 2$ matrix where $U$ is a rotation matrix parameterized by angle $\theta$

$$U=\begin{bmatrix}\cos(\theta)&-\sin(\theta)\\ \sin(\theta)&\cos(\theta)\end{bmatrix},$$

and

$$D=\begin{bmatrix}\xi_{1}&0\\ 0&\xi_{2}\end{bmatrix}$$

is a diagonal matrix scaling the $s_{1}$ and $s_{2}$ coordinate axes separately. Then the pairwise Mahalanobis distance between two locations $\mbox{\boldmath$s$},\mbox{\boldmath$s$}^{\prime}$ is defined as $d=\|A\mbox{\boldmath$s$}-A\mbox{\boldmath$s$}^{\prime}\|$ which is used as the argument to the isotropic Matérn covariance function. A useful interpretation of this form is that if one transforms the coordinates according the linear transformation $A$ then the resulting field will be isotropic.

The SAR model can also be extended to incorporate geometric anisotropy. Let $H$ denote a $2\times 2$ symmetric positive definite anisotropy matrix and modify the Laplacian in the pseudodifferential operator as follows

$$(\kappa^{2}-\nabla\cdot H\nabla)^{\nicefrac{{\alpha}}{{2}}}u(\mathbf{s})=\mathcal{W}(\mathbf{s}).$$

(3)

To avoid potential ambiguity it is helpful to identify the Laplacian operator above for two dimensions in its expanded form as

$$\nabla\cdot H\nabla\equiv H_{1,1}\frac{\partial^{2}}{\partial^{2}s_{1}}+2H_{2,1}\frac{\partial^{2}}{\partial s_{1}\partial s_{2}}+H_{2,2}\frac{\partial^{2}}{\partial^{2}s_{2}}.$$

From this expression, a first-order finite difference discretization of the anisotropic SPDE at (3) gives the following stencil for filling the rows of the $B$ matrix,

$$\begin{array}[]{c|c|c}\frac{2H_{12}}{h_{s_{1}}h_{s_{2}}}&-\frac{H_{22}}{h_{s_{2}}^{2}}&-\frac{2H_{12}}{h_{s_{1}}h_{s_{2}}}\\ \hline\cr-\frac{H_{11}}{h_{s_{1}}^{2}}&\;\;\kappa^{2}+\frac{2H_{11}}{h_{s_{1}}^{2}}+\frac{2H_{22}}{h_{s_{2}}^{2}}&-\frac{H_{11}}{h_{s_{1}}^{2}}\\ \hline\cr-\frac{2H_{12}}{h_{s_{1}}h_{s_{2}}}&-\frac{H_{22}}{h_{s_{2}}^{2}}&\frac{2H_{12}}{h_{s_{1}}h_{s_{2}}}\\ \end{array}$$

(4)

where $h_{s_{1}}$ and $h_{s_{2}}$ are the grid spacings along the x-axis and y-axis. This is just a reparameterization of the results in Appendix A of [20] but facilitates the practical translation between these models. Moreover, setting $h_{s_{1}}=h_{s_{2}}=1$, $H_{12}=H_{21}=0$, and $H_{11}=H_{22}=1$ yields the first-order isotropic model from (1).

Finally, we connect the role of $H$ in the SPDE formulation to the anisotropic model for the Matérn. Under the linear transformation $A=D^{-1}U^{T}$ from Section 2.1, let $\mbox{\boldmath$s$}^{*}=A^{-1}\mbox{\boldmath$s$}$, and let $u$ be an isotropic field solution to the SPDE with Laplacian $\nabla\cdot\nabla$. Furthermore, set $u^{*}(\mbox{\boldmath$s$})=u(A^{-1}\mbox{\boldmath$s$})=u(\mbox{\boldmath$s$}^{*})$. Then from elementary properties of the gradient

$$\nabla u^{*}(\mbox{\boldmath$s$})=\nabla\left(u(A^{-1}\mbox{\boldmath$s$})\right)=\left.A^{-1}\nabla u(A^{-1}\mbox{\boldmath$s$})\right|_{\mbox{\boldmath$s$}=A\mbox{\boldmath$s$}^{*}}=\left.A^{-1}\nabla u(\mbox{\boldmath$s$}^{*})\right|_{\mbox{\boldmath$s$}^{*}=A^{-1}\mbox{\boldmath$s$}}$$

and so we have

$$\nabla\cdot\nabla u^{*}(\mbox{\boldmath$s$})=(A^{-1}\nabla)\cdot A^{-1}\nabla u(\mbox{\boldmath$s$}^{*})=\left.\nabla\cdot A^{-T}A^{-1}\nabla u(\mbox{\boldmath$s$}^{*})\right|_{\mbox{\boldmath$s$}^{*}=A^{-1}\mbox{\boldmath$s$}}.$$

From this expression we identify $H=A^{-T}A^{-1}$. From Section 2.1, if $u$ is an isotropic field then $u^{*}$ will be anisotropic with coordinates transformed by $A^{-1}$. Moreover, $u^{*}$ will also be the solution to the SPDE with $H=A^{-T}A^{-1}$. This connection provides guidance how to interpret $H$. Finally, note that if $A$ is a pure rotation then $H=I$ and isotropy is preserved.

In the climate data analysis below, we find it necessary to include geometric anisotropy in the covariance model. For this reason, we also investigate how the presence of geometric anisotropy affects the numerical correspondence established in 2.4. In this experiment, we encode fixed values for $\xi_{s_{1}}$ and $\xi_{s_{2}}$ in a Matérn correlation matrix such that the length scale ratio $\xi_{1}\colon\xi_{2}=4\colon 1$, which is consistent with the anisotropic estimates in the data analysis. In particular, we let $\xi_{1}=1,\cdots,20$ and $\xi_{2}=4\xi_{s_{1}}$. Then, the optimal eigenvalues of the SAR anisotropy matrix $H$ are found. The experiment was repeated with rotation angles $\theta=0^{\circ},10^{\circ},\cdots,90^{\circ}$, where the rotation angle is assumed to be known and fixed in the eigendecomposition of the $H$ matrix. The anisotropic parameter translation results for $\xi_{s_{1}}$ and $\xi_{s_{2}}$ are shown in panels (a) and (b) of Fig 2, respectively, and the relative error of approximation is shown in panel (c). Overall, the behavior of the approximation is similar to the isotropic case: the estimated eigenvalues of the SAR anisotropy matrix $H$ are smaller than the eigenvalues fixed in the Matérn anisotropy $D^{2}$ matrix.

The approximation results may be slightly affected by the rotation angle and oblateness of the geometric anisotropy, but the effect is negligible in practice. From these results, we have ascertained a numerical translation among the anisotropy parameters. We can use these results to translate locally estimated anisotropic Matérn model parameters into SAR parameters with improved accuracy over the conjectured analytic relationship.

A non-stationary SAR model can be constructed by allowing the parameters $\kappa,H$, and $\sigma^{2}$ in the generating SPDE to vary over space. Let

$${\cal{L}}(\mbox{\boldmath$s$})=H_{1,1}(\mbox{\boldmath$s$})\frac{\partial^{2}}{\partial s_{1}^{2}}+2H_{2,1}(\mbox{\boldmath$s$})\frac{\partial^{2}}{\partial s_{1}\partial s_{2}}+H_{2,2}(\mbox{\boldmath$s$})\frac{\partial^{2}}{\partial s_{2}^{2}}.$$

The SPDE can then be written as

$$(\kappa^{2}(\mathbf{s})-{\cal{L}}(\mbox{\boldmath$s$}))^{\nicefrac{{\alpha}}{{2}}}u(\mathbf{s})=\mathcal{W}(\mathbf{s})$$

where $\kappa(\mathbf{s})>0$, $\mathcal{W}(s)\sim\operatorname{WN}(0,\sigma^{2}(\mathbf{s}))$, and $\sigma^{2}(\mathbf{s})>0$. Furthermore, we specialize to a spatially varying linear transformation of the coordinates, $A(\mbox{\boldmath$s$})$, and so $H(\mbox{\boldmath$s$})=A^{-T}(\mbox{\boldmath$s$})A^{-1}(\mbox{\boldmath$s$})$. Note that $A(\mbox{\boldmath$s$})$ varying in space is equivalent to specifying spatial fields for $\theta$, $\xi_{s_{1}}$ and $\xi_{s_{2}}$ within $U$ and $D$.

Discretizing this equation results in a valid GMRF that is non-stationary. In particular, the autoregressive $B$ matrix from Section 2.2 could have different elements in each row based on the variation in $H(\mbox{\boldmath$s$})$ or $\kappa(\mbox{\boldmath$s$})$. However, $B$ will still be a sparse matrix and $Q=B^{T}B$ will always be positive-definite. The process variance can also be allowed to vary in the same way as with the non-stationary Matérn model, but this must be done balancing the identifiability of $\kappa$ and $H$ and the fact that edge effects may introduce spurious variation in the GMRF variance. Our approach is to first construct the precision matrix and then use sparse matrix methods to solve for the diagonal elements of the covariance matrix. The rows of $B$ are then weighted so that this new version gives a GMRF with constant marginal variance. With this normalization of the SAR model $\sigma(\mbox{\boldmath$s$})^{2}$ can be introduced to capture explicit spatial variation in the process marginal variance.

Estimating a non-stationarity model can be challenging due to the increased number of covariance parameters. When enough data is available, however, local estimation can give insight into what type of non-stationarity is present. Moreover, we illustrate in this section that a modest number of replicated fields results in stable local covariance estimates.

Local estimation is usually accompanied by the assumption of approximate local stationarity. For this work, we define local stationarity and the local likelihood estimation technique for a Gaussian process with stationary Matérn covariance as follows. First, divide the region of interest $\mathcal{D}$ into $M$ possibly overlapping subregions, or windows, $\mathcal{D}_{1},\mathcal{D}_{2},\cdots,\mathcal{D}_{M}$. Then the assumption of approximate local stationarity is that we can model the data $\mbox{\boldmath$y$}_{i}$ within the subregion $\mathcal{D}_{i}$ using a stationary Gaussian process $Y_{i}$ defined using the following specification:

$$Y_{i}(\mathbf{s})=\mu_{i}(\mathbf{s})+Z_{i}(\mathbf{s})+\varepsilon_{i}(\mathbf{s})$$

(5)

where $\varepsilon_{i}$ is mean zero spatial white noise with variance $\tau_{i}^{2}$, and $Z_{i}$ is a mean zero Gaussian process with anisotropic but stationary Matérn covariance function, and $\mu_{i}(\mathbf{s})$ is a fixed mean function. Let $G_{i}=\text{Var}\,\mbox{\boldmath$y$}_{i}$ be the spatial covariance matrix of $\mbox{\boldmath$y$}_{i}$ for the $i$th region. Then the local log likelihood based on $p$ independent replicates $\{\mbox{\boldmath$y$}_{ij}\}_{j=1}^{p}$ for the $i$th region is, up to a constant,

$$\frac{p}{2}\log|G_{i}|^{-1}-\sum_{j=1}^{p}\frac{1}{2}(\mbox{\boldmath$y$}_{ij}-\boldsymbol{\mu}_{i})^{T}G_{i}^{-1}(\mbox{\boldmath$y$}_{ij}-\boldsymbol{\mu}_{i})$$

(6)

where $\boldsymbol{\mu}_{i}$ is the mean function $\mu_{i}(\mathbf{s})$ evaluated at the locations of $\mbox{\boldmath$y$}_{i}$. The likelihood is approximate because we are assuming stationarity within each data window.

After partitioning the data, finding each local likelihood estimate is an embarrassingly parallel task, which makes it a viable strategy for large data sets using many computational cores. In fact, in our application the parallelization is efficient to the point that we take the subregions to be an exhaustive set of moving windows centered at every grid point. We assign these estimated parameters to the location of the center of the subregion $\mathcal{D}_{i}$, and after translating into the SAR parameterization these parameters specify the row of $B$, the SAR matrix, at this location. This assignment is, of course, predicated on the assumption that over the region there is little variation in these parameters. This issue will be discussed in more detail in the last section.

Given that the SAR model also gives a specification of the covariance it may seem indirect that the local estimates focus on the Matérn model, and then the estimates are transformed into the SAR representation. An alternative would be to estimate the SAR version directly from local likelihood windowing. There are several reasons for the two-step procedure. Fitting the covariance model directly avoids any boundary effects that would come about by applying the SAR to a small window. Furthermore, the Matérn parameters are easier to interpret and will be more simple to model in a hierarchical statistical framework.

A practical issue for a local approach, especially in the context of determining covariance parameters, is whether the number of replicates and the size of the window are adequate for robust estimation of parameters. Although choosing a data adaptive window is beyond the scope of this work, it is important to identify the conditions under which parameter estimates will be accurate. Moreover, it is also useful to understand the benefits of replicate spatial fields in estimating a covariance model. In particular, the hope is that replication makes it possible to estimate correlation ranges that are much larger than the local window size.

We perform a Monte Carlo experiment with four factors: window size ranging between a $5\times 5$ grid and a $33\times 33$ grid, the Matérn range parameter being multiples of $1,2,3,$ and $4$ times the window size, the Matérn smoothness parameter taking on values $1$ and $2$, and the number of replicates ranging between 5 and 60. Thus the full factorial design is $11\times 4\times 2\times 9$ (window size $\times$ range parameter $\times$ smoothness parameter $\times$ replicates). For each combination, replicates with given range and smoothness were simulated with given window size, and the Matérn range was estimated using maximum likelihood. This was done 100 times for each combination, and statistics were assembled from the 100 independent maximum likelihood estimates (MLEs) for the range parameter. The main quantity of interest is the percent error of the estimate, and so percent error surfaces as a function of replicate number and window size are summarized in Figure 3. In this experiment, the range parameter was varied based on the window size which may seem unusual. However, the motivation was to address the computational requirements of the problem: given a computational budget to accommodate windows of a specific size, what size range parameter can be accurately estimated? Note that with constraints on the window size, accuracy can also be improved by increasing the number of replicates.

The surfaces in Figure 3 can be used as guidelines to decide how many replicates are necessary and what window size should be used to achieve a specific estimation error tolerance, given something is known about the size of the range to be estimated. The results are encouraging: only a small number of replicates ($>10$) are needed with a window size of $>10$ to estimate a range of 10. In the extreme case, a Matérn range four times the size of the window might be estimated to within $10\%$ error if 30 replicates are available and using a window size of 10 or greater. Using these guidelines, we can be more confident that local moving window likelihood estimation is a viable technique if enough data is used.

We experiment with moving window local MLEs using window sizes between $8\times 8$ and $15\times 15$. Among these choices there was little change in the estimates and subsequent analysis uses a $9\times 9$ window. The estimation was performed on the NCAR Cheyenne supercomputer [4] using the R programming language [26] with the Rmpi [33] and fields packages [6]. The details of the parallel implementation are the same as in [23]. Since the fields were standardized, the constraint $\sigma^{2}=1-\tau^{2}$ was included.

The estimates for the spatially varying parameters are shown in Figure 4. The variance and nugget are shown in (a) and (b). Panel (c) shows the geometric mean of $\xi_{s_{1}}$ and $\xi_{s_{2}}$ as a measure of the average correlation range, and this also agrees with the range in the isotropic case. Finally in panel (d), the estimated anisotropy matrix $A(\mathbf{s}_{i})$ is depicted by glyphs indicating the range and departure from isotropy. The large signal to noise ratio $\sigma^{2}/\tau^{2}$ (not shown) and the evident transition in the covariance structure between land and ocean suggests that the non-stationarity in the second-order structure of the data is being accurately estimated. Based on the coastlines in some regions, we hypothesize that the addition of a land/ocean covariate may also be useful.