Improving Bayesian Local Spatial Models in Large Data Sets

TL;DR

This paper introduces a Bayesian three-step local modeling approach that enhances inference and prediction in large, non-stationary spatial data by effectively managing the bias-variance trade-off through smaller, more manageable regions.

Contribution

A novel Bayesian three-step method enabling smaller regional models without increasing variance, improving inference and prediction in large spatial datasets.

Findings

Improved inference accuracy demonstrated in simulated wind speed data.

Effective propagation of uncertainty across modeling steps.

Enhanced prediction performance over traditional methods.

Abstract

Environmental processes resolved at a sufficiently small scale in space and time will inevitably display non-stationary behavior. Such processes are both challenging to model and computationally expensive when the data size is large. Instead of modeling the global non-stationarity explicitly, local models can be applied to disjoint regions of the domain. The choice of the size of these regions is dictated by a bias-variance trade-off; large regions will have smaller variance and larger bias, whereas small regions will have higher variance and smaller bias. From both the modeling and computational point of view, small regions are preferable to better accommodate the non-stationarity. However, in practice, large regions are necessary to control the variance. We propose a novel Bayesian three-step approach that allows for smaller regions without compromising the increase of the variance that would follow. We are able to propagate the uncertainty from one step to the next without issues caused by reusing the data. The improvement in inference also results in improved prediction, as our simulated example shows. We illustrate this new approach on a data set of simulated high-resolution wind speed data over Saudi Arabia.

Equations

$$y_r(t)$$

$$x_r(t)$$

$${\mathbf{y}}_r\mid {\mathbf{x}}_r,{\mathbf{\theta }}_r$$

$${\mathbf{x}}_r\mid {\mathbf{\theta }}_r$$

$${\mathbf{\theta }}_r$$

$$\pi ({\mathbf{x}}_r,{\mathbf{\theta }}_r\mid {\mathbf{y}}_r)\propto \pi ({\mathbf{\theta }}_r)\pi ({\mathbf{x}}_r\mid {\mathbf{\theta }}_r)\prod _{i=1}^{N_r}\pi \{y_r(i)\mid x_r(i),{\mathbf{\theta }}_r\}.$$

$$\pi \{x_r(i)\mid {\mathbf{y}}_r\}$$

$$\pi \{{\theta }_r(j)\mid {\mathbf{y}}_r\}$$

$$\pi ({\mathbf{\theta }}_r\mid {\mathbf{y}}_r)\approx {\frac{\pi (\mathbf{y}\mid \mathbf{x},\mathbf{\theta })\pi (\mathbf{x}\mid \mathbf{\theta })\pi (\mathbf{\theta })}{\tilde{\pi }(\mathbf{x}\mid \mathbf{y},\mathbf{\theta })}|}_{x=x^{\ast }(\mathbf{\theta })}=\tilde{\pi }({\mathbf{\theta }}_r\mid {\mathbf{y}}_r),$$

$$\tilde{\pi }\{x_r(i)\mid {\mathbf{y}}_r\}$$

$$\tilde{\pi }\{{\theta }_r(j)\mid {\mathbf{y}}_r\}$$

$${\tilde{\pi }}_{\text{mbox}smooth}\{x_r(i)\mid {\mathbf{y}}_r\}$$

$${\tilde{\pi }}_{\text{mbox}smooth}\{{\theta }_r(j)\mid {\mathbf{y}}_r\}$$

$${\mathbf{y}}_r\mid {\mathbf{x}}_r,{\varphi }_r\sim {\mathcal{N}}_T({\mathbf{x}}_r,{\tau }^{-1}{\text{I}}_T),$$

$$\pi ({\mathbf{x}}_r\mid {\varphi }_r)\sim {\mathcal{N}}_T(0,{\mathbf{Q}}_{x,r}^{-1}),$$

$${\Delta }^2{u}_r=u_r-2u_{r+1}+u_{r+2}\sim \mathcal{N}(0,{\tau }_u^{-1}),r=1,\dots ,R-2,$$

$$\text{mbox}{KL}_r=\int \tilde{\pi }({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\varphi }_r)\log \left\{\frac{\tilde{\pi }({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\varphi }_r)}{\pi ({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\varphi }_r)}\right\}\mathrm{d}{\mathbf{x}}_r,$$

$$\pi ({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\varphi }_r)$$

$$\displaystyle\propto\mbox{exp}\bigg{(}-\frac{1}{2}\mathbf{x}_{r}^{\top}\mathbf{Q}_{x,r}\mathbf{x}_{r}\bigg{)}\times\mbox{exp}\bigg{\{}-\frac{1}{2}\tau(\mathbf{x}_{r}^{\top}\mathbf{x}_{r}-2\mathbf{y}_{r}^{\top}\mathbf{x}_{r})\bigg{\}}$$

$$\displaystyle=\mbox{exp}\bigg{\{}-\frac{1}{2}\mathbf{x}_{r}^{\top}(\mathbf{Q}_{x,r}+\tau\mbox{I})\mathbf{x}_{r}+\tau\mathbf{y}_{r}^{\top}\mathbf{x}_{r}\bigg{\}}$$

$$\displaystyle=\mbox{exp}\bigg{\{}-\frac{1}{2}\mathbf{x}^{\top}\mathbf{P}_{r}\mathbf{x}_{r}+\mathbf{b}_{r}^{\top}\mathbf{x}_{r}\bigg{\}},$$

$$\mbox{KL}_{r}=\frac{1}{2}\bigg{\{}\mbox{log}\frac{|\bm{\Sigma}_{1,r}|}{|\bm{\Sigma}_{0,r}|}-T+\mbox{tr}(\bm{\Sigma}_{1,r}^{-1}\bm{\Sigma}_{0,r})+(\bm{\mu}_{1,r}-\bm{\mu}_{0,r})^{\top}\bm{\Sigma}_{1,r}^{-1}(\bm{\mu}_{1,r}-\bm{\mu}_{0,r})\bigg{\}},$$

$$\text{mbox}{CPO}_r(t)=\pi \{y_r(t)|{\mathbf{y}}_r(-t)\}=\int \pi \{y_r(t)|{\mathbf{y}}_r(-t),{\theta }_r\}\pi \{{\theta }_r|{\mathbf{y}}_r(-t)\}\mathrm{d}{\theta }_r,$$

$$y_r({\mathbf{s}}_i)={\mathbf{z}}_r({\mathbf{s}}_i{)}^{\top }{\mathbf{\beta }}_r+u_r({\mathbf{s}}_i)+{ϵ}_r({\mathbf{s}}_i),i=1,\dots ,N_r,$$

$${\mathbf{y}}_r|{\mathbf{\beta }}_r,{\mathbf{u}}_r,{\tau }_{ϵ,r}\sim {\mathcal{N}}_{N_r}({\mathbf{Z}}_r{\mathbf{\beta }}_r+{\mathbf{u}}_r,{\tau }_{ϵ,r}^{-1}{\mathbf{I}}_{N_r}),$$

$$\text{mbox}cov\{u_r({\mathbf{s}}_1),u_r({\mathbf{s}}_2)\}=C_r(h)={\sigma }_{u,r}^2\frac{1}{\Gamma ({\nu }_r)2^{{\nu }_r-1}}({\kappa }_{r}h{)}^{{\nu }_r}{\mathcal{K}}_{{\nu }_r}({\kappa }_{r}h),$$

$${\mathbf{\theta }}_r=({\theta }_{1,r},{\theta }_{2,r},{\theta }_{3,r}{)}^{\top }=\{\log ({\tau }_{ϵ,r}),\log ({\tau }_{u,r}),\log ({\rho }_r){\}}^{\top }.$$

$${\tilde{\pi }}_{\text{mbox}smooth}\{x_r(i)\mid {\mathbf{y}}_r\}$$

$$\text{mbox}log\{\tilde{\pi }({\varphi }_r)\}$$

$$\text{mbox}log\{{\tilde{\pi }}_{\text{mbox}smooth}({\varphi }_r)\}$$

$$\pi ({\mathbf{y}}_r\mid {\varphi }_r)=\frac{\pi ({\mathbf{y}}_r,{\mathbf{x}}_r\mid {\varphi }_r)}{\pi ({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\varphi }_r)},$$

$$\pi ({\mathbf{x}}_r\mid {\varphi }_r)\sim {\mathcal{N}}_T(0,{\mathbf{Q}}_{x,r}),$$

$${\mathbf{Q}}_{x,r}=\left(\begin{array}{cccccc}1& -{\varphi }_r& & & & \\ {\varphi }_r& 1+{\varphi }_r^2& -{\varphi }_r& & & \\ & \cdots & \cdots & \cdots & & \\ & & & -{\varphi }_r& 1+{\varphi }_r^2& -{\varphi }_r\\ & & & & -{\varphi }_r& 1\end{array}\right).$$

$$\pi ({\mathbf{x}}_r,{\varphi }_r\mid {\mathbf{y}}_r)$$

$$\displaystyle\propto\pi(\phi_{r})|\mathbf{Q}_{x,r}|^{1/2}\tau^{1/2}\mbox{exp}\bigg{[}-\frac{1}{2}\bigg{\{}\mathbf{x}_{r}^{\top}\mathbf{Q}_{x,r}\mathbf{x}_{r}+\tau(\mathbf{y}_{r}-\mathbf{x}_{r})^{\top}(\mathbf{y}_{r}-\mathbf{x}_{r})\bigg{\}}\bigg{]}.$$

$$\pi ({\mathbf{x}}_r\mid {\mathbf{y}}_r,{\theta }_r)$$

Full text

Improving Bayesian Local Spatial Models in Large Data Sets

Amanda Lenzi111 Statistics Program, King Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi Arabia., Stefano Castruccio222 Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA.

This publication is based on research supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No: OSR-2018-CRG7-3742., Håvard Rue1, and Marc G. Genton1

March 10, 2024

Abstract

Environmental processes resolved at a sufficiently small scale in space and time inevitably display non-stationary behavior. Such processes are both challenging to model and computationally expensive when the data size is large. Instead of modeling the global non-stationarity explicitly, local models can be applied to disjoint regions of the domain. The choice of the size of these regions is dictated by a bias-variance trade-off; large regions will have smaller variance and larger bias, whereas small regions will have higher variance and smaller bias. From both the modeling and computational point of view, small regions are preferable to better accommodate the non-stationarity. However, in practice, large regions are necessary to control the variance. We propose a novel Bayesian three-step approach that allows for smaller regions without compromising the increase of the variance that would follow. We are able to propagate the uncertainty from one step to the next without issues caused by reusing the data. The improvement in inference also results in improved prediction, as our simulated example shows. We illustrate this new approach on a data set of simulated high-resolution wind speed data over Saudi Arabia.

Keywords: Integrated nested Laplace approximation, latent processes, local models, spatial models, wind speed

Short title: Local Spatial Models

1 Introduction

The rising popularity of statistical methods for environmental data calls for the development of new methods that are able to capture the underlying varying dependencies and to provide computationally efficient inference for the ever increasing amount of data. Traditional geostatistical approaches are not only computationally intensive but are also based on stationarity assumptions, which is convenient but too restrictive and rarely realistic. For instance, wind at sufficiently small temporal resolution (e.g., hourly or sub-hourly) tends to be more variable over complex terrain than over flat surfaces due to geographical features creating eddies. Additionally, failing to account for how physical processes such as weather patterns vary over time or space can lead to an unrealistic assessment of the dependence, and hence suboptimal inference and prediction.

Traditionally, methods have focused on characterizing the spatial and spatio-temporal non-stationarity explicitly via the covariance function. The deformation method in Sampson and Guttorp, (1992) constructs a non-stationary covariance structure from a stationary structure by re-scaling the spatial distance, which was subsequently extended to the Bayesian context in Damian et al., (2001) and Schmidt and O’Hagan, (2003). Another class of non-stationary methods is built on the process convolution or kernel smoothing method, introduced by Higdon, (1998), which uses a spatially varying kernel and a white noise process to create the covariance structure. Other well-known approaches to model non-stationarity include representing the covariance function as a linear combination of basis functions and modelling the covariance matrix of the random coefficients (Nychka et al.,, 2002), and to account for the effect of covariate information directly in the covariance function (Schmidt et al.,, 2011; Neto et al.,, 2014). For a review on the existing literature on non-stationary methods, see Risser, (2016).

Although all of the above methods produce valid models, their computational burden for inference and prediction can be unfeasible for large data sets. Indeed, for evaluating a Gaussian likelihood in a data set of size $n$, $O(n^{2})$ entries need to be stored and $O(n^{3})$ flops need to be computed for the log-determinant and matrix factorization. This task is feasible in modern computers only when $n$ is at most a few tens of thousands of points. Additionally, evaluating a non-stationary model implies inference on a larger parameter space, which requires an exponentially increasing number of likelihood evaluations for frequentist inference or posterior sampling (Edwards et al.,, 2019). To address the difficulties in computation for large data sets, Nychka et al., (2018) used a multi-resolution representation of Gaussian processes to represent non-stationarity based on windowed estimates of the covariance function under the assumption of local stationarity, and successfully used this idea to emulate fields from climate models. Kuusela and Stein, (2018) proposed modelling Argo profiling float data using locally stationary Gaussian process regression, where parameter estimation and prediction were carried out in a moving window. Other works related to moving window methods have been developed and applied in Hammerling et al., (2012) and Tadić et al., (2015) to model remote sensing data.

The seminal work of Lindgren et al., (2011) predicated avoiding modeling the covariance function altogether and modeled the data via a Stochastic Partial Differential Equation (SPDE) instead. By considering a spatial field as a solution of an SPDE, and describing the covariance function only implicitly, inference is of the order $O(n^{3/2})$ (Rue et al.,, 2017), thus allowing inference on considerably larger data sets than covariance-based methods. The computational benefits arise from the precision matrix (inverse covariance matrix) resulting from the approximate stochastic weak solutions of the SPDE, which has a Markovian structure where only close neighbours are non-zero (Rue and Held,, 2005). By spatially varying the coefficients in the SPDEs, it is also possible to construct a variety of non-stationary models. Bolin and Lindgren, (2011) developed such a method for global ozone mapping, whereas Bakka et al., (2019) defined a continuous solution to an SPDE with spatially varying coefficients for solving problems that involve a physical barrier to spatial correlation. By combining the SPDE representation of a stationary Matérn field with the deformation method, Hildeman et al., (2019) modelled non-stationarity in significant wave heights. Locally non-stationary fields were considered in Fuglstad et al., 2015a by letting the coefficients in the SPDE vary with position, and further discussed and generalized for spatially varying marginal standard deviations and correlation structure in Fuglstad et al., 2015b . More recently, Fuglstad and Castruccio, (2020) formulated a global SPDE model with locally varying coefficients with a change of structure across land and ocean. Another application of the SPDE approach to model non-stationarity is to include covariates directly into the model parameters; see Ingebrigtsen et al., (2014) for an application to annual precipitation in Norway.

The aim of this paper is to develop a new method for modelling large data sets with spatial dependence that not only improves local models in terms of inference and prediction, but is also computationally affordable. As a motivating example, we use the high-resolution simulated wind data from a computer model displayed in Figure 1(a). We partition this data into several small disjoint subsets of the data, which we call ‘regions’, as shown in Figure 1(b). Modeling and predicting such variable over a large region present several challenges. First, the data structure at this high resolution is very complex, with details and features that are difficult to capture with a single model. As a consequence, the assumption of stationarity for the entire region is inappropriate. Second, because of the large number of locations, we need a method that is computationally efficient. We show that our method is able to address not only the modeling challenges arising from the inherent non-stationarity of hourly wind, but also the computational issues that are implied by the large data size.

When choosing the size of these regions, we face the conflicting issue of bias-variance trade-off in parameter estimation. Ideally, one wants to choose regions that accurately capture the features in the data (low variance), but also have high predictive out-of-sample skills (low bias). Indeed, small regions reduce the model bias and allow fast computations, at the expense of low accuracy (high variance) in the parameter estimation. Large regions instead allow a control of the variance but also imply a sub-optimal characterization of the dependence structure, hence a bias.

We propose a novel three-step approach, which simultaneously allows for small regions and low variance. The key is to allow small regions to model the local dependence, and correct the estimated parameter distribution with a smoothing step that borrows strength from neighboring regions. The smoothing step is performed so that it accounts for the uncertainty of the parameter estimates from the first step. The resulting smoothed distribution represents the adjusted uncertainty of the local parameters, which is then used for re-fitting the models. Allowing this adjusted uncertainty to be used as a new prior would imply the incorrect premise of the model being influenced by the data twice, hence our approach restricts the information propagation by including it as the new posterior estimates instead. Crucially, the approach we propose is computationally fast and scalable to massive spatial data sets, as it can be fully parallelized across regions. We start with a simple example where the new posterior is the mode of the distribution from the smoothing step. Then, using the wind data in Figure 1(a), we show that it is possible improve the predictive performances by also allowing the uncertainty to propagate from one step to the next.

Our three-step approach is best exemplified by considering a toy data set, where each region consists of an autoregressive process of order one, AR(1). We simulate $R$ time series from this model, where each time series contains $T$ observations, $\mathbf{y}_{r}=\{y_{r}(1),\ldots,y_{r}(T)\}^{\top}$. For each $r$, the observations $\mathbf{y}_{r}$ are assumed to be conditionally independent, given the latent Gaussian random field $\mathbf{x}_{r}=\{x_{r}(1),\ldots,x_{r}(T)\}^{\top}$ and the hyperparameter $\phi_{r}$:

[TABLE]

where $t=2,\ldots,T$ is an index for time, $|\phi_{r}|<1$ and $\tau$ is the fixed precision (known and the same for all time series). Figure 2 shows the different values of $\phi_{r}$ used to simulate $R=100$ time series from (1), where $\phi_{r}$ changes according to a series of sine squared (black squares in Figure 2). For each time series, we set $T=50$ and two different values for the precision: $\tau=2$ and $\tau=1$ in Figure 2(a) and 2(b), respectively. In the first step, we estimate local models for each time series (red circles in Figure 2). In the second step, we apply a correction on the parameters’ estimates from the first step, based on information from neighbouring regions (blue triangles in Figure 2). The third step consists of re-fitting the model in (1) to each time series, propagating the information from the adjusted posterior estimates from the second step back into the analysis. Figure 2 shows that our correction improves the parameter estimates substantially not only for the more extreme case where $\tau=1$ in panel (b), but also when $\tau=2$ in panel (a). More details on this example will be provided in Section 3.

The remainder of this paper is organized as follows. In Section 2 we provide an overview of the proposed methodology. Further details of our approach using the AR(1) example are given in Section 3. The application to the wind speed data in Figure 1 is presented in Section 4. A comprehensive discussion and conclusions are provided in Section 5.

2 Overview of the proposed methodology

2.1 Background

We consider a non-stationary and possibly very large data set, and a partition of the domain into regions where the assumption of stationarity is plausible, defined as $\Omega_{r},r=1,\ldots,R$, where each observation is associated with exactly one $\Omega_{r}$. Each region contains $N_{r}$ observations, $\mathbf{y}_{r}=\{y_{r}(1),\ldots,y_{r}(N_{r})\}^{\top}$. For each $\Omega_{r}$, consider the following hierarchical structure:

[TABLE]

where $\mathbf{x}_{r}=\{x_{r}(1),\ldots,x_{r}(N_{r})\}^{\top}$ is the vector of the latent field that describes the underlying spatial dependence structure, $\bm{\theta}_{r}$ is the $m$-dimensional vector of hyperparameters and $\pi$ is a generic distribution. The observations $\mathbf{y}_{r}$ are assumed to be conditionally independent, given $\mathbf{x}_{r}$ and $\bm{\theta}_{r}$. The resulting joint posterior distribution of $\mathbf{x}_{r}$ and $\bm{\theta}_{r}$ is given by

[TABLE]

Our main goal is to extract the posterior marginal distributions for the elements of the latent field, $\pi\{x_{r}(i)\mid\mathbf{y}_{r}\}$ and hyperparameters, $\pi\{\theta_{r}(j)\mid\mathbf{y}_{r}\}$, and use them to obtain predictive distributions at unsampled locations. Calculation of these univariate posterior distributions requires integrating with respect to $\mathbf{x}_{r}$ and $\bm{\theta}_{r}$:

[TABLE]

where $\bm{\theta}_{r}(-j)$ is the vector of all but the $j$-th hyperparameter component omitted. When the integrals in (3) cannot be found analytically, approximations are typically obtained via simulation-based methods such as MCMC. Alternatively, Rue et al., (2009) proposed an approximate Bayesian inference approach that has become increasingly popular in the last decade. Approximations for $\pi(x_{r}(i)\mid\mathbf{y}_{r},\bm{\theta}_{r})$ and $\pi(\bm{\theta}_{r}\mid\mathbf{y}_{r})$ are obtained via a Laplace approximation (see Rue et al., (2017) for a comprehensive review on this approximation). The posterior $\pi(\bm{\theta}_{r}\mid\mathbf{y}_{r})$ is computed as

[TABLE]

where $\tilde{\pi}(\mathbf{x}\mid\mathbf{y},\bm{\theta})$ is a Laplace approximation, and $x^{*}(\bm{\theta})$ is the mode of $\mathbf{x}$ for a specific value of $\bm{\theta}$. Similarly we obtain $\tilde{\pi}(x_{r}(i)\mid\mathbf{y}_{r},\bm{\theta}_{r})$, the approximation of $\pi(x_{r}(i)\mid\mathbf{y}_{r},\bm{\theta}_{r})$. These are then used to construct the following nested approximations

[TABLE]

2.2 Improving the local estimates

We propose a new method for improving the estimation of $\tilde{\pi}\{\bm{\theta}_{r}\mid\mathbf{y}_{r}\}$ in (4) and hence also improving the estimated $\tilde{\pi}\{x_{r}(i)\mid\mathbf{y}_{r}\}$, for $i=1,\ldots,N_{r}$. Since each region is selected to be small enough to approximate the local non-stationarity well, the resulting parameters’ estimates are likely to have a large variance, and smoothing across the regions is used to reduce it.

The method is based on two extra steps in the estimation procedure from the previous section. In Step 2 we apply a correction to the posteriors $\tilde{\pi}(\bm{\theta}_{r}\mid\mathbf{y}_{r})$ by smoothing the mode of this distribution across $r$. In Section 3, we show a one dimensional example with a smoothing spline, while in Section 4.3 we describe the two dimensional case with a spatial model. We denote by $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\bm{\theta}_{r}\mid\mathbf{y})$ the resulting smoothed distribution for region $r$ in Step 2 of our approach, where $\mathbf{y}$ is the combined data sets from all regions, i.e., $\mathbf{y}=(\mathbf{y}_{1}^{\top},\ldots,\mathbf{y}_{R}^{\top})^{\top}$. In Step 3, the correction from Step 2 is propagated back into the analysis as the posterior for each region:

[TABLE]

where $\tilde{\pi}({x}_{r}(i)\mid\mathbf{y}_{r},\bm{\theta}_{r})$ is obtained by plugging values of $\bm{\theta}_{r}$ from $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\bm{\theta}_{r}\mid\mathbf{y}_{r})$ obtained in Step 2. Step 3 is very computationally efficient, since the posteriors for the hyperparameters have already been estimated, and as in Step 1 the models for each region can be fully parallelized. Also, as the posterior marginals in (5) are the basis to derive the predictive distributions, the proposed correction will also have a direct impact in prediction performance.

Here, the vector $\bm{\theta}_{r}$ contains the hyperparameters that need to be smoothed, while the ones that do not require the smoothing are included in $\mathbf{x}_{r}$. In practice, it is more important to smooth hyperparameters that have a higher variability and are harder to estimate.

Our approach has a crucial difference compared to Empirical Bayes methods. The key is to account for the information from the smoothing in Step 2 directly into the posterior distribution in Step 3, as opposed to introducing it through priors as in Empirical Bayes methods. By doing so, we prevent the estimation in Step 3 to be influenced by the likelihood of the data that was already used in Step 1, and thus avoiding using the data twice. Moreover, our approach allows for uncertainty propagation from Step 2 to Step 3.

3 Simulation with spatially varying AR(1) process

3.1 Model description

In the Introduction we briefly introduced our method on a simulated example (see Figure 2) using the AR(1) model in (1). Here, we provide all the details about the methodology in light of the steps proposed in the previous section. For the ease of exposition, we fix the precision $\tau$ in (1), so that for each region $r$ only the hyperparameter $\phi_{r}$ needs to be estimated. No covariates or additional random effects have been included in (1), but the steps below can be easily adapted to account for them.

The model is a special case of the hierarchical framework proposed in (2). Indeed for the first equation of the hierarchy, the likelihood of the data $\mathbf{y}_{r}$ given the latent field $\mathbf{x}_{r}$ and the hyperparameter $\phi_{r}$ is given by

[TABLE]

where $\textbf{I}_{T}$ is the $T\times T$ identity matrix and $\tau$ is the fixed precision, while $\mathcal{N}_{T}$ is a $T$-dimensional normal distribution. For the latent process $\mathbf{x}_{r}$, we assume that the marginal distribution of $x_{r}(1)$ is Gaussian with mean zero and variance $1/(1-\phi_{r}^{2})$ to have a stationary process. The joint distribution can be written as

[TABLE]

where $\mathbf{Q}_{x,r}$ is the tridiagonal precision matrix of an AR(1) process.

The three steps of our approach can be summarized as follows:

Step 1: The model fitted to each region. Fit the AR(1) model in (1) with fixed known $\tau$ to each time series, $\mathbf{y}_{r}$, separately. Following the notation in Section 2, we define the variance-stabilizing transformation \bm{\theta}_{r}=\theta_{r}=\mbox{log}\Big{(}\frac{1+\phi_{r}}{1-\phi_{r}}\Big{)}, where $\theta_{r}$ has a normal prior with mean zero and precision 0.15, independent across $r$. We then obtain the posterior marginal distributions for the latent field and for the hyperparameter $\theta_{r}$, which we denote by $\tilde{\pi}\{x_{r}(t)\mid\mathbf{y}_{r}\}$ and $\tilde{\pi}(\theta_{r}\mid\mathbf{y}_{r})$, respectively, for $t=1,\ldots,T$ and $r=1,\ldots,R$. Inference is performed using the R-INLA package (Rue et al.,, 2009).

Step 2: Smoothing the hyper-parameter. As in Lindgren and Rue, (2008), we assume a continuous spline on a discrete set of knots with a second order random walk RW(2). We denote by ${\hat{\theta}_{r}}$ the mode for $\tilde{\pi}(\theta_{r}\mid\mathbf{y}_{r})$ from Step 1, and we assume a normal distribution: ${\hat{\theta}_{r}}\sim\mathcal{N}(u_{r},\tau_{\theta;r}^{-1})$, where $\tau_{\theta;r}$ is the precision and is such that $\mbox{log}(\tau_{\theta;r})=\mbox{log}(1/\widehat{\mbox{sd}}_{r}^{2})$, where $\widehat{\mbox{sd}}_{r}$ is the estimated standard deviation of the posterior distribution $\tilde{\pi}(\theta_{r}\mid\mathbf{y}_{r})$. The vector ${\bm{u}}=({u_{1}},\ldots,{u_{R}})^{\top}$ is assumed to have independent second-order increments:

[TABLE]

where $\tau_{u}$ is the precision parameter and can be used to control the degree of smoothing across regions. Section 3.2 discusses a method for choosing the optimal value of $\tau_{u}$.

Step 3: Re-fit the model to each region using the estimated mode. For each region $r$, we assume that the posterior distribution for the hyperparameters, namely $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\theta_{r}\mid\mathbf{y}_{r})$, is a point mass concentrated at the mode of $\hat{\theta}_{r}$ from Step 2. Our choice was dictated by ease of exposition, and in the wind data application in Section 4.3 we will show a more general approach with integration points and weights instead of just the mode. The marginal posterior for the latent process $\mathbf{x}_{r}$ is then obtained from the first equation in (5). Because here there are no hyperparameters that need to be re-estimated in this example, re-fitting the model is equivalent to updating the posterior for $\mathbf{x}_{r}$ given the data under the smoothing from Step 2.

Step 3 implies a change of the original posterior in Step 1, and hence a change in the prior of the model. While retrieving the appropriate prior is not relevant for our method, it is still however possible, and in the Appendix we show the steps to do so. Figure 3 shows (a) the log posterior distributions, (b) log likelihood function, and (c) log prior distributions from Step 1 (solid red) and Step 2 (dashed blue), for $\phi_{r}=0.88$. The log prior distributions were obtained simply by subtracting the log likelihood from the log posterior distributions, and the vertical line represents the true value. The proposed smoothing in Step 2 concentrates the posterior (and consequently the prior) considerably closer to the true value $\phi_{r}$ than a standard approach with no smoothing. Similar results can be observed for other choices of $\phi_{r}$; the mean absolute error across $r$ of the estimated mode posterior distributions from Steps 1 and 2 are $0.23$ and $0.08$, whereas the mean absolute error of the estimated mode priors for these steps are $0.61$ and $0.09$, respectively.

3.2 Sensitivity of prediction to smoothing

There are different approaches to control the degree of smoothness in Step 2. This can be, for instance, dictated by the case study and prior knowledge. Here, we present one possible method, which is based on two metrics: the first focuses on the departure of the estimated posterior against the exact simulated distribution, and the second is based on cross-validation.

To assess the improved accuracy in capturing the true distribution of the latent process $\mathbf{x}_{r}$, for each value of $\phi_{r}$, we calculate the Kullback-Leibler Divergence (KL), a widely used metric for comparing two probability distributions. The departure from the true posterior $\pi(\mathbf{x}_{r}\mid\mathbf{y}_{r},\phi_{r})$ is defined as

[TABLE]

where $\phi_{r}$ in $\tilde{\pi}$ is calculated using either Step 1 or Step 2 of our approach. A small KLr indicates a small departure from the target posterior, and a zero KLr indicates that the two distributions are the same.

The data are simulated from a known model, and the posterior distribution of the latent process $\pi(\mathbf{x}_{r}\mid\mathbf{y}_{r},\phi_{r})$ can be easily obtained from the joint distribution $\pi(\mathbf{x}_{r},\phi_{r}\mid\mathbf{y}_{r})$:

[TABLE]

where, $\mathbf{P}_{r}=\mathbf{Q}_{x,r}+\tau\mbox{I}$ and $\mathbf{b}=\mathbf{y}_{r}^{\top}\tau$. This implies that $\pi(\mathbf{x}_{r}\mid\mathbf{y}_{r},\phi_{r})\sim\mathcal{N}_{T}(\bm{\mu}_{0,r},\bm{\Sigma}_{0,r})$, with $\bm{\mu}_{0,r}=\mathbf{P}_{r}^{-1}\mathbf{b}_{r}$ and $\bm{\Sigma}_{0,r}=\mathbf{P}_{r}^{-1}$. We also assume that the approximated posterior in (7) is normal, i.e., $\pi_{\mbox{\scriptsize{appr}}}(\mathbf{x}_{r}\mid\mathbf{y}_{r},\phi_{r})\sim\mathcal{N}_{T}(\bm{\mu}_{1,r},\bm{\Sigma}_{1,r})$, and we obtain $\bm{\mu}_{1,r}$ and $\bm{\Sigma}_{1,r}$ based on the sample mean vector and covariance matrix from 10,000 posterior samples. The KL divergence expression in (7) can be simplified in the case of two multivariate Gaussian distributions. Indeed, if the target distribution is $\mathcal{N}_{T}(\bm{\mu}_{0,r},\bm{\Sigma}_{0,r})$ and the approximation is $\mathcal{N}_{T}(\bm{\mu}_{1,r},\bm{\Sigma}_{1,r})$, we have

[TABLE]

where $|\bm{\Sigma}|$ denotes the determinant of $\bm{\Sigma}$. Since the KL changes across different orders of magnitudes, we opted for a variance stabilizing estimator, the Expected Mean Log KL (EMLKL) divergence, defined as $\mbox{EMLKL}=\exp\left\{\frac{1}{R}\sum_{r=1}^{R}\log(\mbox{KL}_{r})\right\}$

Therefore, we assess the impact of smoothing on the prediction skills of the estimated process. We use the conditional predictive ordinate (CPO) for leave-one-out cross-validation, defined as

[TABLE]

where $\mathbf{y}_{r}(-t)$ represents the vector of observations $\mathbf{y}_{r}$ with the $t$-th component omitted. In other words, $\mbox{CPO}_{r}(t)$ is calculated by first obtaining the predictive distribution at $t$ given all but the $t$-th observation in the time series, and then evaluating it at the actual withheld value $y_{r}(t)$. The CPO can be interpreted as a continuous equivalent of the posterior probability that the observation is predicted from the model, so larger values are preferable. The CPO can be computed efficiently without re-running the model $R\times T$ times (Held et al.,, 2010). The CPOs are then aggregated in an overall score for comparing different models by averaging across time and regions. As with the KL, we propose the Expected Mean Log Conditional Predictive Ordinate (EMLCPO), defined as \mbox{EMLCPO}=\mbox{exp}\bigg{[}\frac{1}{RT}\sum_{r=1}^{R}\sum_{t=1}^{T}\mbox{log}\{\mbox{CPO}_{r}(t)\}\bigg{]}, with models having relatively higher values of EMLCPO, showing a better fit.

We compare the EMLKL and the EMLCPO based on six different degrees of smoothing by changing the values of $\tau_{u}$ in (6): $\mbox{log}(\tau_{u})=\{-5,-1,3,7,11,15\}$, with lower values of $\mbox{log}(\tau_{u})$ indicating less smoothing. Here, $\mbox{log}(\tau_{u})=15$ results in a constant value across the regions (complete smoothing), so no larger values are considered. Figure 4 shows the results based on (a) KL and on (b) CPO according to the various degrees of smoothing. The first value in the $x$-axis, ‘no smooth’, corresponds to the estimates directly from Step 1 of our approach. According to the EMLKL (panel (a), left $y$-axis) and the EMLCPO (panel (b)), the best fit occurs when $\mbox{log}(\tau_{u})=7$ and $\mbox{log}(\tau_{u})=-5$, respectively. For the EMLKL there is a minimal difference between the $\mbox{log}(\tau_{u})=3$ and $\mbox{log}(\tau_{u})=7$, and the right $y$-axis highlights how the first choice results in less variable KL divergences. Both scores show that there is a clear improvement against a model with no smoothing for $\mbox{log}(\tau_{u})=\{-5,-1,3,7\}$. After $\mbox{log}(\tau_{u})=7$ the posteriors are oversmoothed and this worsens the fit compared to no smoothing (high EMLKL and low EMLCPO values). Evidence from this numerical study suggests that smoothing almost always improves the estimation of the latent process and prediction. The overall agreement between EMLKL and EMLCPO is essential, as, in a real application, the actual underlying distribution is unknown, and a cross-validation metric, such as the EMLCPO, would be used for choosing the optimal degree of smoothness.

Smoothing does not just improve the prediction and decrease the bias, but also results in less variable estimates. Figure 4a (right $y$-axis) shows the spread of KLr for the different amounts of smoothing, displayed as a boxplot. It is readily apparent that optimal smoothing results in more stable estimates by decreasing the variance across regions.

4 Application to the WRF data set

In this section, we apply the approaches in Section 2 to model and predict a sizeable simulated wind speed data set in Saudi Arabia. The predictive ability at the sub-grid scale is of interest for statistical downscaling. Interpolated wind from the numerical simulation can be used as a baseline to build a mathematical relationship (e.g., pattern scaling) from in-situ ground wind data at the same location, hence allowing to generate more accurate, observation-driven wind maps

We apply our method to a spatial data set of simulated wind speed detailed in Section 4.1. In Section 4.2 we present the local model that is fitted to each region and in Section 4.3 we explain in details each step of our approach and present the results.

4.1 The WRF data set

We focus on a simulation generated by Yip, (2018) from the Weather Research and Forecasting (WRF) model, which is a state-of-the-art Numerical Weather Prediction model, developed at the National Center for Atmospheric Research and also recently used in wind energy assessment in Tagle et al., (2020). Mesoscale numerical models such as WRF rely on large-scale atmospheric phenomena or meteorological reanalysis to provide boundary conditions and solve physical equations driving the real processes on a fine scale. The boundary conditions used to simulate the WRF data are obtained from the Modern-Era Retrospective analysis for Research and Applications (MERRA, Rienecker et al., (2011)), a reanalysis product developed at NASA’s Global Modeling and Assimilation Office, using the Goddard Earth Observing System Version 5 general circulation model, together with satellite and surface observations through a data assimilation system.

Each simulation corresponds to hourly data of the zonal and meridional ($U$ and $V$) wind components on a regular grid of $769\times 659$ points in space (5-km resolution) bounded by 5-$35^{\circ}$N and 30-$65^{\circ}$E during the 2009-2014 period, at 2 meters above ground level. The full data set comprises of 506,771 spatial locations. We select data that fall inside Saudi Arabia from 03/06/2010 between 14:00 and 15:00 local time, when wind speeds tend to peak, resulting in 84,494 points in space. The $U$ and $V$ components are converted into wind speed: $\sqrt{U^{2}+V^{2}}$. Figure 1a shows the map of the wind field.

We first partition the domain into $R$ regions small enough so that the assumption of stationarity is plausible. The disjoint subsets are obtained using the $k$-means clustering method, which minimizes the sum of squares from points to the assigned region centers (Hartigan and Wong,, 1979). It is, in principle, possible to provide a more formal assessment of stationarity and use it as a metric for clustering, for example, by fitting directional variograms to each region. However, because these estimates are corrected with a smoothing step, the clustering method is less critical. Our partition results into $R=2000$ regions, (see Figure 1b), with the smallest region containing 26 locations ($\approx 28\times 28$ km) and the largest 62 ($\approx 48\times 48$ km).

4.2 The spatial model

The distribution of wind speed is bounded below by zero and is significantly right-skewed. Therefore, wind speed cannot be directly modeled with the Gaussian distribution. Common transformations for normalizing wind speed data include logarithmic transformation and square-root transformation (Taylor et al.,, 2009). Haslett and Raftery, (1989) showed that square-root transformation is well suited for wind data, as the resulting transformed wind speed resembles the Gaussian distribution. Hence, for each region $r$ we model the square-root transformed wind speed $y_{r}$ at sampling locations $\mathbf{s}=(\mathbf{s}_{1},\ldots,\mathbf{s}_{N_{r}})$ with a latent Gaussian model, a special case of the hierarchical framework proposed in (2). For each region $r$, we assume

[TABLE]

where $\mathbf{z}_{r}$ is a $p$-dimensional vector of covariates, and $\bm{\beta}_{r}$ is the vector of the linear coefficients. Here, $\{\epsilon_{r}(\mathbf{s}_{1}),\ldots,\epsilon_{r}(\mathbf{s}_{N_{r}})\}\sim\mathcal{N}_{N_{r}}(0,\tau_{\epsilon,r}^{-1}I_{N_{r}})$ is the iid random noise that accounts for the model uncertainty. The aforementioned model can be written in the vector form

[TABLE]

where $\mathbf{y}_{r}=\{y(\mathbf{s}_{1}),\ldots,y(\mathbf{s}_{N_{r}})\}^{\top}$ is the observation vector and the $N_{r}\times P$ design matrix is $\mathbf{Z}_{r}=\{\mathbf{z}_{r}(\mathbf{s}_{1}),\ldots,\mathbf{z}_{r}(\mathbf{s}_{N_{r}})\}^{\top}$. We consider $p=2$, thus two covariates: elevation and distance to the coast. In terms of the hierarchical framework in Section 2.1, (8) is the first equation, i.e., the data level, in (2).

The spatial field $u_{r}(\mathbf{s}_{i})$ is assumed to be Gaussian and isotropic, with a covariance described by the Matérn function, a widely popular choice in spatial statistics. For two locations $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ at distance $h=\|\mathbf{s}_{1}-\mathbf{s}_{2}\|$, the Matérn covariance is defined as (Stein,, 1999)

[TABLE]

where $\sigma^{2}_{u,r}=1/\tau_{u,r}$ is the marginal variance and ${\cal{K}}_{\nu_{r}}$ is the modified Bessel function of the second kind of order $\nu_{r}>0$. The popularity of the Matérn is mainly attributable to the control of the number of mean square derivatives of the underlying process through the parameter $\nu_{r}$. The range is controlled by $\kappa_{r}>0$ and $\rho_{r}=\sqrt{8\nu_{r}}/\kappa_{r}$ represents the distance at which the spatial correlation is approximately $0.13$, and we set $\nu_{r}=1$.

The vector of hyperparameters to be estimated is given by the precision of the data, the precision of the latent process, and its range, so that

[TABLE]

The linear coefficients $\bm{\beta}_{r}$ in (8) are less variable, so they are not included in the vector of hyperparameters to be smoothed.

We provide a joint distribution for the range $\rho_{r}$ and the variance $\sigma^{2}_{u,r}$ using the concept of the Penalized Complexity (PC) prior that was recently introduced by Simpson et al., (2017). PC develops priors that allow shrinkage towards a base model, which is assumed to be the reference. The prior is then built by allowing a control of the KL divergence from the base to the actual model. Following Fuglstad et al., (2019), we assume a base model with infinite range and precision, i.e., a constant, and we assign PC priors to $\rho_{r}$ and $\tau_{u,r}$ that are able to control the tail probabilities: $\mbox{P}(\sigma^{2}_{u,r}>\sigma^{2}_{0,r})=\alpha_{1}$ and $\mbox{P}(\rho_{r}<\rho_{0,r})=\alpha_{2}$. We choose $\alpha_{1}=\alpha_{2}=0.01$, $\rho_{0,r}$ to be the $20\%$ of the range of the observations and $\sigma^{2}_{0,r}$ the variance estimated from the data at region $r$. In other words, we assume a prior that bounds the variance to be larger than that estimated from the data with a $1\%$ chance, and the range to be below $20\%$ of the range of the observations with a $1\%$ chance. For $r=1,\ldots,R$, we assume a vague Gamma prior with parameters 1 and 0.00005 for $\tau_{\epsilon,r}$ and a vague Gaussian prior $\mathcal{N}(0,1000)$ for $\bm{\beta}_{r}$. The priors are also assumed to be independent across components. The R-INLA package is used for model fitting and predictions (Rue et al.,, 2009).

4.3 Results

We now detail our approach with the data and the model described in the previous sections. The three steps are described as follows:

Step 1: The model fitted to each region.

We fit the model outlined in Section 4.2 to each of the $R=2000$ regions in Figure 1b separately, and obtain estimates of the posterior distribution for the $k$-th element of $\bm{\theta}_{r}$ for $k=1,2,3$, which we denote by $\tilde{\pi}(\theta_{k,r}\mid\mathbf{y}_{r})$. We denote as $\hat{\theta}_{k,r}$ the mode of $\tilde{\pi}(\theta_{k,r}\mid\mathbf{y}_{r})$, while the posterior standard deviation is denoted as $\widehat{\mbox{sd}}_{k,r}$. We show the results for $\theta_{3,r}=\log(\rho_{r})$, since the range is the hardest parameter to identify, and hence the most variable across regions. Figure 5a shows the maps of $\hat{\theta}_{3,r}$. Many regions have a considerably higher estimated posterior mode than the neighboring regions, hence smoothing is necessary. Figure 5b shows the map of the posterior standard deviation $\widehat{\mbox{sd}}_{3,r}$, and it is apparent how the locations with large range values correspond to the ones with low posterior variance. The high variance in the estimates of Figure 5 a is a consequence of the small region size needed to accommodate the non-stationarity. The region size is another tuning parameter of the method, and cross-validation could have been used to choose the optimal region size.

Step 2: Smoothing the hyperparameters.

The modes $\hat{\theta}_{k,r}$ from Step 1 are smoothed here independently across $k$ for simplicity and are normalized by subtracting the mean and dividing by its standard deviation computed across $r=1,\ldots,R$, for each $k$. For a general Matérn in region $r$ with known smoothness $\nu_{r}$, under infill asymptotics, only $\sigma_{r}^{2}\kappa_{r}^{2\nu_{r}}$ can be consistently estimated (Zhang,, 2004). As the likelihood and the posteriors can be concentrated around non-linear manifolds of the parameter space, modelling the hyperparameters in the log-scale alleviates the problem of smoothing them separately. With an abuse of notation, we now refer to $\bm{\theta}_{r}$ and their components as their normalized version. We assume an additive model for smoothing: $\hat{\theta}_{k,r}(\mathbf{s}_{c})=u_{k}(\mathbf{s}_{c})+\varepsilon_{k,r}(\mathbf{s}_{c})$, where the locations $\mathbf{s}_{c}$ are the centroids of each region $r$. The process $u_{k}(\mathbf{s}_{c})$ is assumed to be Gaussian and modeled with the Matérn covariance in (9), with marginal variance $\tilde{\sigma}^{2}_{u,k}$, range $\tilde{\rho}_{k}$, and the iid noise is $\varepsilon_{k,r}\sim\mathcal{N}(0,\tilde{\tau}^{-1}_{\epsilon;k,r})$, for $k=1,2$ and 3.

We assume $\tilde{\tau}_{\epsilon;k,r}$ to be fixed at the value of $1/\widehat{\mbox{sd}}_{k,r}^{2}$, $r=1,\ldots,R$, from Step 1. This ensures that the same degree of smoothness is applied to all three additive models, i.e., the hyperparameters with a larger standard deviation will be smoothed more than the ones with a smaller standard deviation. Here, $\tilde{\rho}_{k}$ is fixed to half of the domain of the study region. A choice of considerably different values, such as the size of the domain, would result in oversmoothing. The choice of $\tilde{\tau}_{u}=1/\tilde{\sigma}^{2}_{u,k}$ is performed via cross-validation and will be discussed later. Because $\hat{\theta}_{k,r},k=1,2,3$, are at the same scale after normalization, we can use the same smoothness and therefore $\tilde{\tau}_{u}$ will not be strongly dependent on $k$. We use six equally spaced values for log($\tilde{\tau}_{u}$), varying from $-7.5$ to $5$. The fitted values from the smoothing are then transformed back from the normalized to the original scale. Figure 5c is an example of the estimated posterior mode of $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\theta_{3,r}\mid\mathbf{y}_{r})$ with $\mbox{log}(\tilde{\tau}_{u})=-5$.

Step 3: Re-fit the model to each region using integration points.

In the AR(1) simulated example in Section 3, the smoothed hyperparameter posterior was assumed to be a point mass concentrated at the smoothed posterior mode from Step 2, so that calculation of $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}\{x_{r}(i)\mid\mathbf{y}_{r}\}$ in (5) was trivial. In this application, we propose a more articulated method which numerically approximates the integral in the first equation of (5).

We use the Gauss–Hermite quadrature, a numerical scheme to approximate integrals of the form $\int e^{-\xi^{2}}f(\xi)\mathrm{d}\xi\approx\sum_{l=1}^{L}f\{\xi^{(l)}\}\Delta^{(l)}$ for a fixed $L$. The abscissas for the quadrature of order $L$, which are given by the roots of the Hermite polynomials $\xi^{(l)}$, and the weights $\Delta^{(l)}$, both have a closed form expression (Abramowitz and Stegun,, 1964).

We operate under the assumption that $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\bm{\theta}_{r}\mid\mathbf{y})$ can be well approximated by a product of marginal normal distributions $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\bm{\theta}_{r}\mid\mathbf{y})\approx\prod_{k=1}^{3}\mathcal{N}(\mu_{k,r},\sigma^{2}_{k,r})$, where $\mu_{k,r}$ and $\sigma^{2}_{k,r}$ are the posterior mean and variance of $\tilde{\pi}_{\text{smooth}}({\theta}_{k,r}\mid\mathbf{y})$, respectively. The independence implied by the product is made for convenience, although empirically, we found a relatively low correlation between the components of $\bm{\theta}_{r}$. Also, because of the log-scale, the posteriors can be well approximated by a Gaussian distribution, and the first expression in (5) becomes

[TABLE]

where the latent field $\mathbf{x}_{r}=(\mathbf{u}^{\top}_{r},\bm{\beta}^{\top}_{r})^{\top}$ contains the linear coefficients and the spatial process in (8). Using a change of variables, we obtain $\xi_{k,r}=\frac{\theta_{k,r}-\mu_{k,r}}{\sqrt{2}\sigma_{k,r}}\Leftrightarrow\theta_{k,r}=\mu_{k,r}+\sqrt{2}\xi_{k,r}\sigma_{k,r}$. For this case study, $L=5$ integration points in each of the three dimensions provide an approximation that is sufficiently accurate. Thus, the required number of configurations to evaluate the integral is $L^{3}=5^{3}=125$. Since each configuration can be evaluated independently, the computations can be easily parallelized.

4.4 Choice of the smoothing parameter

There is no true underlying model here, so the EMLKL in Section 3.2 is not applicable and we only focus on the cross-validation score EMLCPO. We compare the leave-one-out predictive performance using the different degrees of smoothing, as explained previously in Section 3. Figure 6 shows this comparison: lower values of $\mbox{log}(\tilde{\tau}_{u})$ indicate more smoothing than higher values. The highest value corresponds to the results obtained directly from Step 1. The EMLCPO value attains its maximum at $\mbox{log}(\tilde{\tau}_{u})=-5$, and any of the smoothing levels improves the original estimates from Step 1. Differently from the AR(1) case in Section 3, where at some point the smoothing becomes excessive and the scores progressively deteriorate, here the performance is significantly improved even for a large smoothing. We also compare the predictive performances of the integration method against the approach using only the mode as in the AR(1) case. The Gauss-Hermite integration shows marginal improvement, especially for low degrees of smoothing. For higher degrees of smoothing, the estimated posterior distribution is more narrow, and the effect of the integration is less apparent.

5 Discussion

In this work, we developed a new three-step approach for analyzing large data sets with spatial dependence that improves local models in terms of inference and prediction. The method is scalable to extremely large spatial data sets and can properly propagate the uncertainty across steps in a Bayesian framework. In Step 1, the domain is partitioned into regions, and local models are fit to each region. The size of these regions is a bias-variance trade-off; larger regions will have a smaller variance and more substantial bias, whereas smaller regions will have higher variance and lower bias. We choose to use smaller regions, thus allowing the capture of local non-stationarities, followed by a correction for the high variance, based on borrowing information from neighboring regions in Step 2, while accounting for the uncertainty of the parameter estimates from Step 1. Finally, in Step 3, the model is re-fitted to each region, propagating the uncertainty from the smoothing back into the analysis as the new posterior, thus avoiding problems of using the data twice. The approach allows flexible modeling of complex dependence structures, but is at the same time computationally affordable, as the proposed adjustment is amenable to full parallelization across regions.

In both the AR(1) simulated data and the application, the improvement from our method compared to fitting local models to each region is apparent. Indeed, the smoothing adjustment allows us to better recover the actual posterior distribution in the simulation study, and most importantly, it enables a superior predictive skill. The smoothing can be chosen to achieve the best possible advantage over the uncorrected model. Ad-hoc sensitivity analysis shows that our method is robust concerning the smoothing technique, with improved results for a wide range of smoothings. The existing methods for the nonstationary case involve highly complex model fitting strategies, and only a few software is available, see, e.g., Gramacy, (2007) and Risser and Calder, (2017). However, a detailed comparison with other methods was not feasible due to the large size of our simulated wind data set, which could not be presently handled by either of these packages.

Our method is general and can be applied to many settings: space, time, space/time, and different domains, as long as a partition is provided. It relies on local models defined through a hierarchical latent process framework, a class large enough to allow a wide range of applications. If better local models are provided, our method can still be used to correct the variance of the estimated parameters.

A limitation of this approach lies in the assumption of a domain partition. For some applications such as wind, the regions imply a discontinuity at the border, and hence prediction at unsampled locations at the border might be suboptimal. Partition-based approaches that do not imply independence across blocks of the partition with a globally valid model are available. Guinness and Stein, (2013) and Castruccio and Guinness, (2017) proposed an evolutionary spectrum model to capture high-frequency temperature across day and night regimes. If the nature of the problem suggests a change in spatial dependence dictated by some geographical features such as mountain range as proposed in Jeong et al., (2018), then this strategy could be naturally employed with appropriate likelihood approximation. For our domain and wind, however, the partition must be provided by a clustering scheme such as the $k$-means.

An application of our model to spatio-temporal data is possible. Still, it would likely require additional approximations and a careful choice of the regions as the data size and the hyperparameter space will be considerably larger.

Appendix: Retrieving the priors

The re-fitting procedure in Step 3 of our approach uses the information from Step 2 as the new posterior distribution. We show how to retrieve the prior distribution that corresponds to the posterior for the toy example in Section 3.

For each $\phi_{r}$ and corresponding data $\mathbf{y}_{r}$, with $r=1,\ldots,R$, let $\pi(\mathbf{y}_{r}\mid\phi_{r})$ be the likelihood of observing data $\mathbf{y}_{r}$ given the hyperparameter $\phi_{r}$. We denote by $\tilde{\pi}(\phi_{r}\mid\mathbf{y}_{r})$ and $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\phi_{r}\mid\mathbf{y}_{r})$ the posterior distributions from Steps 1 and 3, respectively, and $\tilde{\pi}(\phi_{r})$ and $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\phi_{r})$ are the corresponding priors.

From Bayes’ Theorem, the prior distributions are given by:

[TABLE]

where $A$ is the normalizing constant.

Recall that the posteriors $\tilde{\pi}(\phi_{r}\mid\mathbf{y}_{r})$ and $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\phi_{r}\mid\mathbf{y}_{r})$ in the right hand sides of (10), are readily available from Steps 1 and 2, respectively. Therefore, to evaluate $\tilde{\pi}(\phi_{r})$ and $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\phi_{r})$, what remains to be computed is the likelihood term $\pi(\mathbf{y}_{r}\mid\phi_{r})$, which is the same in both equations given in (10). To compute this term, we start by writing:

[TABLE]

and then compute (11) in two steps:

The joint distribution $\pi(\mathbf{y}_{r},\mathbf{x}_{r}\mid\phi_{r}):$

We assume that the marginal distribution of $x_{r}(1)$ is Gaussian with mean zero and variance $1/(1-\phi_{r}^{2})$. Then, we can express the joint distribution of $\mathbf{x}_{r}$, $\pi(\mathbf{x}_{r}\mid\phi_{r})=\pi\{x_{r}(1)\}\pi\{x_{r}(2)\mid x_{r}(1)\},\ldots,\pi\{x_{r}(T)\mid x_{r}(T-1)\}$, as

[TABLE]

where $\mathbf{Q}_{x,r}$ is the tridiagonal precision matrix of an AR(1) process

[TABLE]

It follows that the joint posterior distribution is

[TABLE] 2. 2.

The conditional distribution $\pi(\mathbf{x}_{r}\mid\mathbf{y}_{r},\phi_{r}):$

We use the fact that the conditional distribution of $\mathbf{x}_{r}$ is just the joint distribution between $\mathbf{x}_{r}$ and $\mathbf{y}_{r}$, without the terms that do not depend on $\mathbf{x}_{r}$ since $\mathbf{y}_{r}$ and $\phi_{r}$ are fixed:

[TABLE]

Using the canonical form of the multivariate Gaussian distribution, we can write:

[TABLE]

where, $\mathbf{P}_{r}=\mathbf{Q}_{x,r}+\tau\mbox{I}$ and $\mathbf{b}_{r}=\mathbf{y}_{r}^{\top}\tau$. It follows that:

[TABLE]

Finally, from (13) and (14), we can write $\pi(\mathbf{y}_{r}\mid\phi_{r})$ in (11) evaluated at $\mathbf{x}_{r}=0$ as

[TABLE]

Next, from the posteriors $\tilde{\pi}(\phi_{r}\mid\mathbf{y}_{r})$ and $\tilde{\pi}_{\mbox{\scriptsize{smooth}}}(\phi_{r}\mid\mathbf{y}_{r})$ on the right hand side of (10) that are computed in Steps 1 and 2, respectively, together with the likelihood term in (15), we can obtain the corresponding priors in (10). The right hand side plot of Figure 3 shows these exact scaled log prior distributions.