A Case Study Competition Among Methods for Analyzing Large Spatial Data

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the “big data” era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low-rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online. Electronic Supplementary Material Supplementary materials for this article are available at 10.1007/s13253-018-00348-w.


Introduction
For decades, the Gaussian process (GP) has been the primary tool used for the analysis of geostatistical (point-referenced) spatial data (Schabenberger and Gotway, 2004;Cressie, 1993;Cressie and Wikle, 2015;Banerjee, Carlin and Gelfand, 2014).A spatial process Y (s) for s ∈ D ⊂ R 2 is said to follow a GP if any realization Y = (Y (s 1 ), . . ., Y (s N )) at the finite number of locations s 1 , . . ., s N follows an N -variate Gaussian distribution.More specifically, let µ(s) : D → R denote a mean function returning the mean at location s (typically assumed to be linear in covariates X(s) = (1, X 1 (s), . . ., X P (s)) ) and C(s 1 , s 2 ) : D 2 → R + denote a positive definite covariance function.Then, if Y (s) follows a spatial Gaussian process, Y has the density function, where µ = (µ(s 1 ), . . ., µ(s N )) is the mean vector and Σ = {C(s i , s j )} ij is the N × N covariance matrix governed by C(s i , s j ) (e.g. the Matérn covariance function).
From this definition, the appealing properties of the Gaussian distribution (e.g.Gaussian marginal and conditional distributions) have rendered the GP an indispensable tool for any spatial data analyst to perform such tasks as kriging (spatial prediction) and proper uncertainty quantification.
With the modern onset of larger and larger spatial datasets, however, the use of Gaussian processes for scientific discovery has been hindered by computational intractability.Specifically, evaluating the density in (1.1) requires O(N 3 ) operations and O(N 2 ) memory which can quickly overwhelm computing systems when N is only moderately large.Early solutions to this problem included factoring (1.1) into a series of conditional distributions (Vecchia, 1988;Stein, Chi and Welty, 2004), the use of pseudo-likelihoods (Varin, Reid and Firth, 2011;Eidsvik et al., 2014), modeling in the spectral domain (Fuentes, 2007) or using tapered covariance functions (Furrer, Genton and Nychka, 2006;Kaufman, Schervish and Nychka, 2008;Stein, 2013).Beginning in the late 2000's, several approaches based on low rank approximations to Gaussian processes were developed (or became popular) including discrete process convolutions (Higdon, 2002;Lemos and Sansó, 2009), fixed rank kriging (Cressie and Johannesson, 2008;Kang and Cressie, 2011;Katzfuss and Cressie, 2011), predictive processes (Banerjee et al., 2008;Finley et al., 2009), lattice kriging (Nychka et al., 2015) and stochastic partial differential equations (Lindgren, Rue and Lindström, 2011).Sun, Li and Genton (2012) and Bradley et al. (2016) provide exceptional reviews of these methods and demonstrate their effectiveness for modeling spatial data.
After several years of their use, however, scientists have started to observe shortcomings in many of the above methods for approximating GPs such as the propensity to oversmooth the data (Simpson, Lindgren and Rue, 2012;Stein, 2014) and even, for some of these methods, an upper limit on the size of the dataset that can be modeled.Hence, recent scientific research in this area has focused on the efficient use of modern computing platforms and the development of methods that are parallelizable.For example, Paciorek et al. (2015) show how (1.1) can be calculated using parallel computing while Katzfuss and Hammerling (2017) and Katzfuss (2017) develop a basisfunction approach that lends itself to distributed computing.Alternatively, Barbian and Assunc ¸ão (2017) and Guhaniyogi and Banerjee (2018) propose dividing the data into a large number of subsets, draw inference on the subsets in parallel and then combining the inferences.Datta et al. (2016a,b) build upon Vecchia (1988) by developing novel approaches to factoring (1.1) as a series of conditional distributions based only on nearest neighbors.
Given the plethora of choices to analyze large spatially correlated data, for this paper, we seek to not only provide an overview of modern methods to analyze massive spatial datasets, but also lightly compare the methods in a unique way.Specifically, this research implements the common task framework of Wikle et al. (2017) by describing the outcome of a friendly case study competition between various research groups across the globe who each implemented their own method to analyze the same spatial datasets (see the list of participating groups in Table 1).That is, several research groups were provided with two spatial datasets (one simulated and one real) with a portion of each dataset removed to validate predictions (research groups were not provided with the removed portion so that this study is "blinded").The simulated data represents a scenario where the Gaussian process assumption is valid (i.e., a correctly specified model), whereas the real dataset is a scenario when the model is potentially misspecified due to inherent non-stationarity or non-Gaussian errors.Each group then implemented their unique method and provided a prediction (and prediction interval or standard error) of the spatial process at the held out locations.The predictions were compared by a third party and are summarized herein.
The case study competition described herein is unique and novel in that, typically, comparisons/reviews of various methods is done by a single research group implementing each method.However, single research groups may be more or less acquainted with some methods leading to a possibly unfair comparison with those methods they are less familiar with.In contrast, for the comparison/competition here, each method was implemented by a research group with strong expertise in the method and who is wellversed in any possible intricacies associated with its use.Hence, in terms of scientific contributions, this paper (i) serves as a valuable review, (ii) discusses a unique case study comparison of spatial methods for large datasets, (iii) provides code to implement each method to practitioners (see supplementary materials) and (iv) establishes a framework for future studies to follow when comparing various analytical methods.
The remainder of this paper is organized as follows.Section 2 gives a brief back- ground on each method.Section 3 provides the setting for the comparison along with background on the datasets.Section 4 then summarizes the results of the comparison in terms of predictive accuracy, uncertainty quantification and computation time.Section 5 draws conclusions from this study and highlights future research areas for the analysis of massive spatial data.
2. Overview of Methods for Analyzing Large Spatial Data

Fixed Rank Kriging
Fixed Rank Kriging (FRK, Cressie andJohannesson, 2006, 2008) is built around the concept of a spatial random effects (SRE) model.In FRK, one models the process Y (s), s ∈ D, as where µ(s) is the mean function that is itself modeled as a linear combination of known covariates, w(s) is the SRE model, and ξ(s) is a fine-scale process, modeled to be spatially uncorrelated with variance σ 2 ξ .The process ξ(s) in (2.4) is designed to soak up variability in Y (s) not accounted for by w(s).
The primary assumption of FRK is that the spatial process w(•) can be decomposed into a linear combination of K basis functions h(s) = (h 1 (s), . . ., h K (s)) , s ∈ D, and K basis function coefficients θ = (θ 1 , . . ., θ K ) such that, (2.2) The use of K basis functions ensures that all estimation and prediction equations only contain inverses of matrices of size K × K, where K N .In practice, the set {h k (•)} in (2.2) is comprised of functions at R different resolutions such that (2.2) can also be written as where h rk (s) is the k th spatial basis function at the r th resolution with associated coefficient θ rk , and K r is the number of basis functions at the r th resolution, such that K = R r=1 K r is the total number of basis functions used.For this research, we used R = 3 resolutions of bisquare basis functions following Cressie and Johannesson (2008).
The coefficients θ = {θ rk } have Var(θ) = S(φ) with covariance parameters φ that need to be estimated.In this work, S(φ) is a block-diagonal matrix composed from R dense matrices, where the r th block has i, jth element exp(−d r (i, j)/φ r ) and d r (i, j) is the distance between the centroids of the i th and j th basis function at the r th resolution, and φ = (φ 1 , . . ., φ R ) are the spatial correlation parameters of the imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 exponential correlation function.Note that S(φ) can also be unstructured in which case K(K + 1)/2 parameters need to be estimated, however this case is not considered here.
There are several variants of FRK.In this work, we use the implementation by Zammit-Mangion and Cressie (2017) which comes in the form of the R package FRK, available from the Comprehensive R Archive Network (CRAN).In this paper we utlize v0.1.6 of that package.In FRK the model for Y (s), s ∈ D, is composed as in (2.1).FRK further assumes that recorded observations Y (s i ) are noisy readings of where for i = 1, . . ., N , ε(s i ) denotes independent and identically normally distributed measurement error with mean 0 and known measurement error variance σ 2 ε .More details on the implementation of FRK for this study are included in the supplementary materials.

LatticeKrig
LatticeKrig (LK, Nychka et al., 2015) uses nearly the same setup as is employed by FRK.Specifically, LK assumes the model (2.1) and (2.4) but omits the fine-scale process ξ(•).Further, for w(s), LK follows the multiresolution approach in (2.3), but LK uses a different structure and constraints than FRK.First, the marginal variance of each resolution h r (s)θ r where h r (s) = (h r1 (s), . . ., h rKr (s)) are the basis functions of the r th resolution with coefficients θ r = (θ r1 , . . ., θ rKr ) is constrained to be σ 2 w α r where σ 2 w , α r > 0 and R r=1 α r = 1.To further reduce the number of parameters, LK sets α r ∼ r −ν where ν is a single free parameter.
LatticeKrig obtains multiresolution radial basis functions by translating and scaling a radial function in the following manner.Let u rk for r = 1, . . ., R and k = 1, . . ., K r denote a regular grid of K r points on D corresponding to resolution r.For this article, LK defines where the distance is taken to be Euclidean because the spatial region in this case is of small geographic extent and θ r = 2 −r .Further, LK defines which are Wendland polynomials and are positive definite (an attractive property when the basis is used for interpolation).Finally, the basis functions in (2.6) are normalized at each resolution so that the process marginal variance at all s is σ 2 w α r .This reduces edge effects and makes for a better approximation to a stationary covariance function.
LatticeKrig assumes the coefficients at each resolution θ r = (θ r1 , . . ., θ rKr ) are independent (similar to the block diagonal structure used in FRK) and follow a multivariate normal distribution with covariance Q −1 r parameterized by a single parameter φ r .
Because the locations {u rk } Kr k=1 are prescribed to be a regular grid, LK uses a spatial autoregression/Markov random field (see Banerjee, Carlin and Gelfand, 2014, Section 4.4) structure for Q −1 r leading to sparsity and computational tractability.Furthermore, because Q r is sparse, LK can set K to be very large (as in this competition greater than N ) without much additional computational cost.The supplementary material to this article contains additional information about the implementation of LatticeKrig used in this case study.

Predictive Processes
For the predictive process (PP) approach, let s 1 , . . ., s K denote a set of "knot" locations well dispersed over the spatial domain D. Assume that the SREs (w(s)) in (2.1) follow a mean zero Gaussian process with covariance function C(s, s ) = σ 2 w ρ(s, s ) where ρ(•, •) is a postive definite correlation function.Under this Gaussian process assumption, the SREs w = (w(s 1 ), . . ., w(s K )) ∼ N (0, Σ w ) where Σ w is a K ×K covariance matrix with ij th element C(s i , s j ).The PP approach exploits the Gaussian process assumption for the SREs and replaces w(s) in (2.1) with where C(s, s ) = (C(s, s 1 ), . . ., C(s, s K )) .Note that (2.7) can be equivalently written as the basis function expression given above in (2.2) where the basis functions are C(s, s )Σ −1 w and w effectively plays the role of the basis coefficients.Finley et al. (2009) noted that the basis function expansion in (2.7) systematically underestimates the marginal variance σ 2 w from the original process.That is, To counterbalance this underestimation of the variance, Finley et al. (2009) use the structure in (2.4), where ξ(s) are spatially independent with distribution N (0, σ 2 w −C (s, s )Σ −1 w C (s, s )) such that Var( w(s) + ξ(s)) = σ 2 w as in the original parent process.As with FRK and LatticeKrig, the associated likelihood under (2.8)only requires calculating the inverse and determinant of a dense K × K matrix and diagonal N × N matrices which results in massive computational savings when K N and K is small.However, one advertised advantage of using the PP approach as opposed to FRK or LatticeKrig is that the PP basis functions are completely determined by the choice of covariance function C(•, •).Hence, the PP approach is unaltered even when considering modeling complexities such as anisotropy, non-stationarity or even multivariate processes.At the same time, however, when C(•, •) is governed by unknown parameters (which is nearly always the case) the PP basis functions need to be calculated iteratively rather than once as in FRK or LatticeKrig which will subsequently increase computation time.

Spatial Partitioning
(2.9) where By way of distinction, this approach is inherently different from the "divide and conquer" approach (Liang et al., 2013;Barbian and Assunc ¸ão, 2017).In the divide and conquer approach, the full dataset is subsampled, the model is fit to each subset and the results across subsamples are pooled.In contrast, the spatial partition approach uses all the data simultaneously in obtaining estimates, but the independence across regions facilitates computation.
The key to implementing the spatial partitioning approach is the choice of partition and the literature is replete with various options.A priori methods to define the spatial partitioning include partitioning the region into equal areas (Sang, Jun and Huang, 2011), partitioning based on centroid clustering (Knorr-Held and Raßer, 2000;Kim, Mallick and Holmes, 2005) and hierarchical clustering based on spatial gradients (Anderson, Lee and Dean, 2014;Heaton, Christensen and Terres, 2017).Alternatively, model-based approaches to spatial partitioning include treed regression (Konomi, Sang and Mallick, 2014) and mixture modeling (Neelon, Gelfand and Miranda, 2014) but these approaches typically require more computation.For this analysis, a couple of different partitioning schemes were considered, but each scheme resulted in approximately equivalent model fit to the training data.Hence, based on the results from the training data, for the competition below we used an equal area partition of approximately 6000 observations per subregion.

Covariance Tapering
The idea of covariance tapering is based on the fact that many entries in the covariance matrix Σ in (1.1) are close to zero and associated location pairs could be considered as essentially independent.Covariance tapering multiplies the covariance function C(s i , s j ) with a compactly supported covariance function, resulting in another positive definite covariance function but with compact support.From a theoretical perspective, covariance tapering (in the framework of infill-asymptotics) is using the concept of Gaussian equivalent measures and mis-specified covariance functions (see, e.g., Stein, 1999 and references therein).Subsequently, Furrer, Genton and Nychka (2006) have imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 assumed a second-order stationary and isotropic Matérn covariance to show asymptotic optimality for prediction under tapering.This idea has been extended to different covariance structures (Stein, 2013), non-Gaussian response (Hirano and Yajima, 2013) and multivariate and/or spatio-temporal setting (Furrer, Bachoc and Du, 2016).
From a computational aspect, the compact support of the resulting covariance function provides the computational savings needed by employing sparse matrix algorithms to efficiently solve systems of linear equations.More precisely, to evaluate density (1.1), a Cholesky factorization for Σ is performed followed by two solves of triangular systems.For typical spatial data settings, the solve algorithm is effectively linear in the number of observations.
For parameter estimation in the likelihood framework, one-and two-taper approaches exist (see Kaufman, Schervish and Nychka, 2008;Du, Zhang and Mandrekar, 2009;Wang and Loh, 2011;Bevilacqua et al., 2016, for relevant literature).To distinguish the two approaches, notice that the likelihood in (1.1) can be rewritten as where etr(A) = exp(trace(A)).In the one-taper setting, only the covariance is tapered such that Σ in (2.10) is replaced by Σ T where " " denotes the Hadamard product and T is the N × N tapering matrix.In the two-tapered approach both the covariance and empirical covariance are affected such that not only is Σ replaced by Σ T but The one-taper equation results in biased estimates of model parameters while the two-taper approach is based on estimating equations (and is, therefore, unbiased) but comes at the price of a severe loss of computational efficiency.If the one-taper biased estimates of model parameters are used for prediction, the biases may result in some loss of predictive accuracy (Furrer, Bachoc and Du, 2016).
Although tapering can be adapted to better take into account uneven densities of locations and complex anisotropies, we use a simple straight-forward approach for this competition.The implementation here relies almost exclusively on the R package spam (Furrer and Sain, 2010;Furrer, 2016).Alternatively to likelihood approaches and in view of computational costs, we have minimized the squared difference between an empirical covariance and parameterized covariance function.The gridded structure of the data is exploited and the empirical covariance is estimated for a specific set of locations only; and thus is close to classical variogram estimation and fitting (Cressie, 1993).

Multiresolution Approximations
The multi-resolution approximation (MRA) can be viewed as a combination of several previously described approaches.Similar to FRK or LatticeKrig, the MRA expresses the spatial process of interest w(s) in (2.1) as a weighted sum of compactly supported basis functions at different resolutions as in (2.3).In contrast to FRK or LatticeKrig, the MRA basis functions and the prior distribution of the corresponding weights are chosen using the predictive-process approach to automatically adapt to any given covariance imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 A toy example of simulated observations (black dots) with a covariance function C with increasing smoothness on a one-dimensional spatial domain D = [0, 1], together with a multi-resolution approximation (MRA) with M = 4 resolutions with 3 subregions per region (vertical lines) and r 0 = 2 basis functions per region.The basis functions and their weights (symbolized by the height of the functions) adjust to the changing smoothness, here increasing from left to right.
function C(•), and so the MRA can adjust flexibly to a desired spatial smoothness and dependence structure.Scalability of the MRA is ensured in that for increasing resolution, the number of basis functions increases while the support of each function (i.e., the part of the spatial domain in which it is nonzero) decreases.Decreasing support (and increasing sparsity of the covariance matrices of the corresponding weights) is achieved either by increasingly severe tapering of the covariance function (MRA-taper; Katzfuss and Gong 2017) or by recursively partitioning the spatial domain (MRAblock; Katzfuss, 2017).This can lead to (nearly) exact approximations with quasilinear computational complexity.
While the MRA-taper has some attractive smoothness properties, we focus here on the MRA-block which is based on a recursive partitioning of the domain D into smaller and smaller subregions up to some level M .Within each (sub-)region at each resolution, there is a small number, say r 0 , of basis functions.The resulting approximation of the process (including its variance and smoothness) in each region at resolution M is exact.In addition, it is feasible to compute and store the joint posterior covariance matrix (i.e., not just its inverse as with related approaches) for a large number of prediction locations as a product of two sparse matrices.Figure 1 illustrates the MRA basis functions in a toy example.
The MRA-block is designed to take full advantage of high-performance computing systems, in that inference is well suited for massively distributed computing, with limited communication overhead.The computational task is split into small parts by assigning a computational node to each region of the recursive partitioning.The nodes then deal in parallel with the basis functions corresponding to their assigned regions leading to a polylogarithmic computational complexity.For this project, we use M = 9 levels, partition each domain in 2 parts and set the number of basis function in each partition to r 0 = 64.

Nearest Neighbor Processes
The nearest neighbor Gaussian process (NNGP) developed in Datta et al. (2016a) and Datta et al. (2016c) is defined from the conditional specification of the joint distribution of the SREs in (2.1).Let w(s) in (2.1) follow a mean zero Gaussian process with C(s, s ) = σ 2 w ρ(s, s ) where ρ(•) is a positive definite correlation function.Factoring the joint distribution of w(s 1 ), . . ., w(s N ) into a series of conditional distributions yields that w(s 1 ) = 0 + η(s 1 ) and where w 1: ) and η's are independent, mean zero, normally distributed random variables.More compactly, (2.11) is equivalent to w = Aw + η where A = (a ij ) is a lower triangular matrix with zeroes along the diagonal and η = (η(s 1 ), . . ., η( ).This effectuates a joint distribution w ∼ N (0, Σ) where Furthermore, when predicting for any s / ∈ {s 1 , . . ., s N }, one can define similar to (2.11).
A sparse formulation of A ensures that evaluating the likelihood of w (and, hence, of Y ) will be computationally scalable.Because spatial covariances decrease with increasing distance, Vecchia (1988) demonstrated that replacing the conditional set w 1:(i−1) by the smaller set of m nearest neighbors (in terms of Euclidean distance) of s i provides an excellent approximation to the conditional density in (2.11).Datta et al. (2016a) demonstrated that this is equivalent to A having at-most m non-zero entries in each row and thereby corresponds to a proper probability distribution.Similarly, for prediction at a new location s, a sparse a(s) in (2.12) is constructed based on mnearest neighbors of s among s 1 , . . ., s N .The resulting Gaussian Process is referred to as the Nearest Neighbor Gaussian Process (NNGP) and computation primarily involves small m × m matrix operations.Generalizing the use of nearest neighbors from expedient likelihood evaluations as in Vecchia (1988) and Stein, Chi and Welty (2004) to the well defined NNGP on the entire domain enables fully Bayesian inference and coherent recovery of the latent SREs.
Using an NNGP prior for Y (s) − X (s)β, the model can be written as Y ∼ N (Xβ, Σ(φ)) where Σ is the NNGP covariance matrix derived from the full GP.A Bayesian specification is completed by specifying priors for the parameters β and φ.For this application, the covariance function C consists of an stationary exponential GP with variance σ 2 and range φ and a nugget process with variance σ 2 ε (see (2.4)).We assign a normal prior for β, inverse gamma priors for σ 2 w and σ 2 ε and a uniform prior for φ.A Gibbs sampler for the model involves conjugate updates for β and metropolis random walk updates for φ = (σ 2 w , σ 2 ε , φ) .Letting α = σ 2 ε /σ 2 w , the model can also be expressed as Y ∼ N (Xβ, σ 2 w R(φ, α)) where R is the NNGP matrix derived from C(φ) + αI, C(φ) being the correlation imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 matrix of the exponential GP.Fixing α and φ gives a conjugate Normal-Inverse Gamma posterior distribution for β and σ 2 w .Predictive distributions for y(s) at new locations can also be obtained as t-distributions.The fixed values of α and φ can be chosen from a grid-search by minimizing root mean square predictive error score based on K-fold cross validation.This hybrid approach departs from fully Bayesian philosophy by using hyper-parameter tuning.However, it offers a pragmatic solution for massive spatial datasets.We refer to this model as the conjugate NNGP model and the fully Bayesian approach described above as the response NNGP model.Detailed algorithms for both the models are provided in Finley et al. (2017).NNGP models for analyzing massive spatial data are available on CRAN as the R-package spNNGP (Finley, Datta and Banerjee, 2017).

Stochastic PDEs
The stochastic partial differential equation approach (SPDE) is based on the equivalence between Matérn covariance fields and stochastic PDEs, in combination with the Markov property that on 2-dimensional domains holds for integer valued smoothness parameters in the Matérn family.The starting point is a basis expansion for w(s) of the form (2.2), where the basis functions h k (s) are chosen to be piecewise linear on a triangulation of the domain (Lindgren, Rue and Lindström, 2011).The optimal joint distribution for the θ k coefficients is obtained through a finite element construction, which leads to a sparse inverse covariance matrix (precision) Q θ (φ).The precision matrix elements are polynomials in the precision and inverse range parameters (1/φ 2 σ and 1/φ r ), with sparse matrix coefficients that are determined solely by the choice of triangulation.This differs from the sequential Markov construction of the NNGP method which instead constructs a square-root free LDL Cholesky decomposition of its resulting precision matrix (in a reverse order permutation of the elements).
The spatial process is specified through a joint Gaussian model for z = (θ, β) with prior mean 0 and block-diagonal precision , where Q β = I • 10 −8 gives a vague prior for β.Introducing the sparse basis evaluation matrix H with elements H ij = h j (s i ) and covariate matrix X = X j (s i ), the observation model is then Y = Xβ + Hθ + ε.The design matrix for the joint vector z is A = (H, X), and ε is a zero mean observation noise vector with diagonal precision Q ε = I/σ 2 ε .Using the precision based equations for multivariate Normal distributions, the conditional precision and expectation for z are given by , where sparse Cholesky factorisation of Q z|y is used for the linear solve.The elements of z are automatically reordered to keep the Cholesky factors as sparse as possible.The resulting computational and storage cost for the posterior predictions and multivariate Gaussian likelihood of a spatial Gaussian Markov random field of this type with K basis functions is O(K 3/2 ).Since the direct solver does not take advantage of the stationarity of the model, the same prediction cost would apply to non-stationary models.For larger problems, more easily parallelizeable iterative sparse solvers (e.g.multigrid) can be applied, but for the relatively small size of the problem here, the straightforward implementation of a direct solver is likely preferable.The posterior covariance elements of Q −1 z|y corresponding to the non-zero structure of Q z|y imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 are obtained through Takahashi recursions as a post-processing step on the Cholesky factor of Q z|y (see Rue, Martino and Chopin, 2009).These elements are precisely the ones needed to compute the final predictive variances Var[µ(s 0 ) + w(s 0 ) + ε 0 | Y ] for each prediction location s 0 .The triangulation nodes were here chosen to coincide with the observation lattice, and in order to avoid unwanted boundary effects, the triangulation extends a short distance outside the domain.This extension has only a small effect on the computational cost, since the triangles are allowed to be larger than inside the domain of interest, and therefore the extension doesn't need as many nodes as in a regular lattice extension.In addition, the exponential covariance is a Matérn covariance with smoothness 0.5, and hence is not Markovian on R 2 .Where the LK method approaches this by using a sum of several Markovian components, the SPDE implementation in INLA (Rue et al., 2017) instead uses a parsimonious Markovian spectral approximation for a single field.The resulting model is a second order Markov random field on the coefficients {θ k }.For details of the approximation see the authors' response to the discussion of Lindgren, Rue and Lindström (2011).
The implementation of the SPDE method used here is based on the R package INLA (Rue et al., 2017), which is aimed at Bayesian inference for latent Gaussian models (in particular Bayesian generalised linear, additive, and mixed models) using integrated nested Laplace approximations (Rue, Martino and Chopin, 2009).The parameter optimization for φ = (φ r , φ σ , σ 2 ε ) uses general numerical log-likelihood derivatives, thus the full Bayesian inference was therefore turned off, leading to an empirical Bayes estimate of the covariance parameters.Most of the running time is still spent on parameter optimization, but using the same parameter estimation technique as for LK, in combination with a purely Gaussian implementation, substantively reduces the total running time even without specialized code for the derivatives.

Periodic Embedding
When the observation locations form a regular grid, and the model is stationary, methods that make use of the discrete Fourier transform (DFT), also known as spectral methods, can be statistically and computationally beneficial, since the DFT is an approximately decorrelating transform, and it can be computed quickly and with low memory burden using fast Fourier transform (FFT) algorithms.For spatially gridded data in two or higher dimensions-as opposed to time series data in one dimensionthere are two prominent issues to be addressed.The first is edge effects, and the second is missing values.By projecting onto trigonometric bases, spectral methods essentially assume that the process is periodic on the observation domain, which leads to bias in the estimates of the spectrum (Guyon, 1982;Dahlhaus and Künsch, 1987).Guinness and Fuentes (2017) and Guinness (2017) propose the use of small domain expansions and imputing data in a periodic fashion on the expanded lattice.Imputation-based methods also solve the second issue of missing values, since the missing observations can be imputed as well.
The methods presented here follow the iterative semiparametric approach in Guinness (2017).Guinness and Fuentes (2017) provides an alternative parametric approach.
imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 For this section, let N = (N 1 , N 2 ) give the dimensions of the observation grid (in the surface temperature dataset N = (300, 500)).Let τ denote an expansion factor, and let m = τ N denote the size of the expanded lattice.We use τ = 1.2 in all examples, so that m = (360, 600) in the surface temperature dataset.Let U be the vector of observations, and V be the vector of missing values on the grid of size m, making the full vector Y = (U , V ) .The discrete Fourier transform of the entire vector is The procedure is iterative.At iteration k, the spectrum f k is updated with where α is a smoothing kernel, and E k is expected value under the multivariate normal distribution with stationary covariance function where F m is the set of Fourier frequencies on a grid of size m.This is critical since it ensures that R k is periodic on the expanded grid.In practice, the expected value in (2.13) is replaced with |J(ν)| 2 computed using an imputed vector V , a conditional simulation of missing values given U under covariance function R k .This ensures that the imputed vector V is periodic on the expanded lattice and reduces edge effects.The iterative procedure can also be run with an intermediate parametric step in which the Whittle likelihood (Whittle, 1954) is used to estimate a parametric spectral density, which is used to filter the imputed data prior to smoothing the spectrum.See Guinness ( 2017) for details about more elaborate averaging schemes and monitoring for convergence of the iterative method.

Metakriging
Spatial metakriging is an approximate Bayesian method that is not tied to any specific model and is partly algorithmic in nature.In particular, any spatial model described above can be used to draw inference from subsets (as described below).From (1.1), let the N × N covariance matrix be determined by a set of covariance parameters φ such that Σ = Σ(φ) (e.g.φ could represent decay parameters from the Matérn covariance function) and µ(s) = X (s)β where X(s) is a set of known covariates with unknown coefficients β.Further, let the sampled locations S = {s 1 , ..., s N } be partitioned into sets {S 1 , ..., S K } such that S i ∩S j = ∅ for i = j and the corresponding partition of the data be given by {y k , X k }, for k = 1, 2, . . ., K, where each y k is n k × 1 and X k is n k × p. Assume that we are able to obtain posterior samples for Ω = {β, φ} from (1.1) applied independently to each of K subsets of the data in parallel on different cores.To be specific, assume that } is a collection of M posterior samples from p(Ω | y k ).We refer to each p(Ω | y k ) as a "subset posterior."The metakriging approach we outline below attempts to combine, optimally and meaningfully, these subset posteriors to arrive at a legitimate probability density.We refer to this as the "metaposterior".
Metakriging relies upon the unique geometric median (GM) of the subset posteriors (Minsker et al., 2014;Minsker, 2015).We envision the individual posterior densities p k ≡ p(Ω | y k ) to be residing on a Banach space H equipped with norm • ρ .The GM is defined as where y = (y 1 , y 2 , . . ., y K ) .The norm quantifies the distance between any two posterior densities π 1 (•) and π 2 (•) as The GM is unique.Further, the geometric median lies in the convex hull of the individual posteriors, so π * (Ω | y) is a legitimate probability density.Specifically, Computation of the geometric median π * ≡ π * (Ω | y) proceeds by employing the popular Weiszfeld's iterative algorithm that estimates α ρ,k (y) for every k from the subset posteriors p k .To further elucidate, we use a well known result that the geometric median π * satisfies, Since there is no apparent closed form solution for α ρ,k (y) that satisfies this equation, one needs to resort to the Weiszfeld iterative algorithm outlined in Minsker et al. (2014) to produce an empirical estimate of α ρ,k (y) for all k = 1, .., K. Guhaniyogi and Banerjee (2018) show that, for a large sample, π * (• | y) provides desirable approximation of the full posterior distribution in certain restrictive settings.It is, therefore, natural to approximate the posterior predictive distribution p(y(s 0 ) | y) by the subset posterior predictive distributions p(y(s 0 ) | y k ).Let {y(s 0 ) (j,k) } M j=1 , k = 1, . . ., K, be samples obtained from the posterior predictive distribution p(y(s 0 )|y k ) from the k-th subset posterior.Then, Therefore, the empirical posterior predictive distribution of the metaposterior is given by imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 1 y(s0) (j,k) , from which the posterior predictive median and the 95% posterior predictive interval for the unobserved y(s 0 ) are readily available.
One important ingredient of spatial metakriging (SMK) is partitioning the dataset into subsets.For this article, we adopt a random partitioning scheme that randomly divides data into K = 30 exhaustive and mutually exclusive subsets.The random partitioning scheme facilitates each subset to be a reasonable representative of the entire domain, so that each subset posterior acts as a "weak learner" of the full posterior.We have explored more sophisticated partitioning schemes and found similar predictive inference.
For the sake of definiteness, this article uses the Gaussian process model for each subset inference which may lead to higher run time.However, the metakriging approach lends much more scalability when any of the above models is employed in each subset.In fact, ongoing research in spatial metakriging includes distributed spatial kriging (DISK) (Guhaniyogi et al., 2017) which scales the modified predictive process to millions of observations.

Gapfill
The gapfill method (Gerber et al., 2018) differs from the other herein presented methods in that it is purely algorithmic, distribution-free, and, in particular, not based on Gaussian processes.Like other prediction methods popular within the satellite imaging community (see Gerber et al. 2018 andWeiss et al. 2014 for reviews), the gapfill method is attractive because of its low computational workload.A key aspect of gapfill is that it is designed for parallel processing, which allows the user to exploit computing resources at different scales including large servers.Parallelization is enabled by predicting each missing value separately based on only a subset of the data.
To predict the value Y (s 0 ) at location s 0 gapfill first selects a suitable subset A = {Y (s i ) : s i ∈ N (s 0 )}, where N (s 0 ) defines a spatial neighborhood around s 0 .Finding A is formalized with rules, which reassure that A is small but contains enough observed values to inform the prediction.In this study, we require A to have an extent of at least 5 × 5 pixels and to contain at least 25 non-missing values.Subsequently, the prediction of Y (s 0 ) is based on A and relies on sorting algorithms and quantile regression.Moreover, prediction intervals are constructed using permutation arguments (see Gerber et al. 2018 for more details on the prediction and uncertainty intervals).
The gapfill method was originally designed for spatio-temporal data, in which case the neighborhood N (s 0 ) is defined in terms of the spatial and temporal dimensions of the data.As a consequence, the implementation of gapfill in the R package gapfill (Gerber, 2017) requires multiple images to work properly.To mimic this situation, we shift the given images by one, two, and three pixels in both directions along the x and y-axes.Then the algorithm is applied to those 13 images in total (one original image and 12 images obtained through shifts of the original image).

Local Approximate Gaussian Processes
The local approximate Gaussian process (laGP, Gramacy and Apley, 2015) addresses the big-N problem in GP regression by taking a so-called transductive approach to learning, where the fitting scheme is tailored to the prediction problem (Vapnik, 1995) as opposed to the usual inductive approach of fitting first and predicting later conditional on the fit.A special case of laGP, based on nearest neighbors, is simple to describe.In order to predict at s, simply train a Gaussian process predictor on the nearest m neighbors to s; i.e., use the data subset Y m = {Y (s i ) : s i ∈ N m (s)}, where N m (s) are the m closest observed locations to s in terms of Euclidean distance.If the data-generating mechanism is not at odds with modeling assumptions (e.g., having a well-specified covariance structure), then one can choose m to be as large as possible, up to computational limitations, in order to obtain an accurate approximation.Observe that this use of nearest neighbors (NNs) for prediction is more akin to the classical statistical/machine learning variety, in contrast to their use in determining the global (inverse) covariance structure as described in Section 2.7.
Interestingly, NNs do not comprise an optimal data subset for prediction under the usual criteria such as mean-squared error.However, finding the best m of N !/(m!(N − m)!) possible choices represents a combinatorially huge search.The laGP method generalizes this so-called nearest neighbor prediction algorithm (whose modern form in spatial statistical literature is described by Emery 2009) by approximating that search with a greedy heuristic.First, start with a NN set Y m0 (s) = {Y (s i ) : where m 0 < m, and then for j = m 0 + 1, . . ., m successively choose s j to augment Y m0 building up a local design data set one point at a time according to one of several simple objective criteria related to mean-square prediction error.The idea is to repeat in this way until there are m observations in Y m (s).Gramacy and Apley's preferred variation targets s j which maximizes the reduction in predictive variance at s.To recognize a similar global design criterion called active learning Cohn (Cohn, 1996), they dubbed this criterion ALC.Qualitatively, these local ALC designs tend to have a cluster of neighbors and "satellite" points and have been shown to offer demonstrably better predictive properties than NN and even full-data alternatives especially when the data generating mechanism is at odds with the modeling assumptions.The reason is that local fitting offers a way to cope with a certain degree of non-stationarity which is common in many real data settings.
ALC search iterations and GP updating considerations as designs are built up, are carefully engineered to lead to a method whose computations are of O(N 3 ) complexity (i.e., the same as the simpler NN alternative).A relatively modest local design size of m = 50 typically works well.Moreover, calculations for each s are statistically independent of the next, which means that they can be trivially parallelized.Through a cascade of multi-core, multi-node and GPU parallelization, Gramacy, Niemi and Weiss (2014) and Gramacy and Haaland (2016) illustrated how N in the millions, in terms of both training and testing data sizes could be handled (and yield accurate predictors) with less than an hour of computing time.The laGP method has been packaged for R and is available on CRAN (Gramacy, 2016).Symmetric multi-core parallelization (via OpenMP) and multi-node automations (via the built-in parallel package) work out-of-the box.GPU extensions are provided in the source code but require custom Research groups participating in the competition along with their selected method (competitor).A disadvantage to local modeling in this fashion is that a global predictive covariance is unavailable.Indeed, the statistically independent nature of calculation is what makes the procedure computationally efficient and parallelizable.In fact, the resulting global predictive surface, over a continuum of predictive s-locations, need not even be smooth.However in most visual representations of predictive surfaces it can be difficult to distinguish between a genuinely smooth surface and what is plotted via the laGP predictive equations (see Figures 3 and 4 below).Finally, it is worth noting that although laGP is applied here in a spatial modeling setting (i.e., with two input variables), it was designed for computer simulation modeling and has been shown to work well in input dimension as high as ten.

The Competition
At the initial planning phase of this competition, we desired to compare a broad variety of approaches: from frequentist to Bayesian and from well-established to modern developments.In accordance with this plan, efforts were made to contact a variety of research groups with strong expertise in a method to analyze the datasets.After this outreach period, the research teams listed in Table 1 agreed to participate and implement their associated method.
Each group listed in Table 1 were provided with two training datasets: one real and one simulated.The simulated dataset then represented a case where the covariance function was specified correctly while the real dataset represented a scenario where the covariance function was misspecified.Both datasets consisted of observations on the same 500×300 grid ranging longitude values of −95.91153 to −91.28381 and latitude values of 34.29519 to 37.06811.The real dataset consisted of daytime land surface temperatures as measured by the Terra instrument onboard the MODIS satellite on August 4, 2016 (Level-3 data).The data was downloaded from the MODIS reprojection tool web interface (MRTweb) located at https://mrtweb.cr.usgs.gov/and is provided as supplementary material to this article.The latitude and longitude range, as well as the date, were chosen because of the sparse cloud cover over the region on this date (rather than by scientific interest in the date itself).Namely, only 1.1% of the Level-3 MODIS data were corrupted by cloud cover leaving 148,309/150,000 observed values to use for our purposes.
The simulated dataset was created by, first, fitting a Gaussian process model with constant mean, exponential covariance function and a nugget effect to a random sample of 2500 observations from the above MODIS data.The resulting parameter estimates were then used to simulate 150,000 observations on the same grid as the MODIS data.
In order to ensure a realistic analysis scenario, the missing data pattern on August 6, 2016 from the same MODIS satellite data product was used to separate each dataset into training and test sets.Each group independently wrote code (all of which is included as supplementary material to this article) that provided (i) a point prediction for each location in the test set, (ii) a 95% prediction interval for location in the test set or a corresponding standard error for the prediction, (iii) the average time required to implement the method per iteration and (iv) the total clock time needed to implement the method.In order to minimize the number of confounding factors in this competition, each group was instructed to use an exponential correlation function (if applicable to their chosen method) and a nugget variance.For the simulated data the groups were instructed to only use a constant mean (because this was how the data was originally simulated).However, for the satellite data, the groups used a linear effect for latitude and longitude so that the residual process more closely resembled the exponential correlation.The code from each team was then run on the Becker computing environment (256 GB of RAM and 2 Intel Xeon E5-2680 v4 2.40GHz CPUs with 14 cores each and 2 threads per coretotaling 56 possible threads for use in parallel computing) located at Brigham Young University (BYU).Each team's code was run individually and no other processes were simultaneously run so as to provide an accurate measure of computing time.
Each method was compared in terms of mean absolute error ), continuous rank probability score (CRPS; see Gneiting and Raftery 2007;Gneiting and Katzfuss 2014), interval score (INT; see Gneiting and Raftery 2007) and prediction interval coverage (CVG; the percent of intervals containing the true value).To calculate the CRPS, we assumed the associated predictive distribution was well approximated by a Gaussian distribution with mean centered at the predicted value and standard deviation equal to the predictive standard error.In cases where only a prediction interval was provided, the predictive standard error was taken as (U − L)/(2 × Φ −1 (0.975)) where U and L are the upper and lower ends of the interval, respectively.

Results for Simulated Data
The numerical results for the simulated data competition are displayed in Table 2 and the associated predicted surfaces for each method are shown in Figure 3. First, consider the predictive accuracy as measured by the MAE and RMSE in Table 2.In terms of predictive accuracy, each method performed extremely well with the best MAE being 0.61 while the worst was only 1.03.Similarly, the best RMSE was 0.83 compared to a worst RMSE of only 1.31.Considering the range of the simulated data was 53.80 − 33.91 = 19.89, a RMSE of 1.31 is quite accurate.
While all the methods performed well in terms of predictive accuracy, when considering uncertainty quantification (UQ) some of the methods fared better than others.For example, LatticeKrig, LAGP, metakriging, MRA, periodic embedding and NNGP all achieved near the nominal 95% coverage rate.In contrast, FRK, Gapfill, partitioning and PP achieved lower than nominal coverage while SPDE and tapering have higher than nominal coverage.Considering UQ further, Gapfill and PP have large interval  scores suggesting possible wide predictive intervals in addition to the penalty incurred from missing the true value.In this regard, it is important to keep in mind that LAGP, metakriging, MRA, NNGP and PP all can specify the "correct" exponential correlation function.Additionally, LK and SPDE have settings that can approximate the exponential correlation function well.In contrast, some methods such as FRK and Gapfill are less suited to model fields with exponential correlation functions, which may partially explain their relatively poor prediction or coverage performance in this instance.Finally, Figure 3 displays the predictive surfaces for each method on the simulated data.The visual inspection of the predictive surfaces provides interesting insights into the various features of each method.For example, because the Gapfill method was primarily designed for spatio-temporal data, we shifted the images to create "pseudo" datasets for the Gapfill algorithm.However, this shifting resulted in a "smeared" pattern in the predictive surface which we hypothesize would not occur in the space-time setting.Likewise, arguments by Simpson, Lindgren and Rue (2012) and Stein (2014) suggest that low rank methods oversmooth the data and such possible oversmoothing is seen in the predictive surfaces for FRK and PP.

Results for Real Data
The results for the real MODIS data are displayed in Table 3 and largely reiterate the results from the simulated data.Namely, each method performed very well in terms of predictive accuracy.The largest RMSE was only 2.52 which, when considered on the data range of 55.41 − 24.37 = 31.04, is very small.We note that, under the setup of the competition, some of the methods were forced to approximate a GP with isotropic exponential covariance function, which is the true covariance function of the simulated data, but most certainly not for the real data.Thus, the scores are lowest for those approximations that happened to result in a good fit to the data and not necessarily lowest for those methods that best approximated the exponential covariance.
The largest discrepancies among the competing methods is again in terms of uncertainty quantification.Lattice kriging, metakriging, MRA, NNGP and periodic embed-   ding again achieved near nominal coverage rates with small interval scores and CRPS.
The SPDE and tapering approaches did better in terms of coverage in that the empirical rates were near nominal (recall that the corresponding coverage rates were too high for the simulated data for these methods).In contrast, the coverage rates on the MODIS data for FRK, Gapfill, LAGP, partitioning and predictive processes were too small resulting in larger interval scores.Finally, visual inspections of the predictive surfaces for the MODIS data are shown in Figure 4. Notably the majority of the methods smooth out the predictions in the north-central region.This is to be expected because such predictions are considered "long-range" with very little (or no) observed data in this region (see Figure 2).Hence, predictions for this region rely more heavily on the overall mean surface rather than borrowing information from neighboring observations (of which there is none).Again, the "shifting" used for the Gapfill algorithm is again apparent in the predictive surface.As with the simulated data, we hypothesize that such "smeared" predictive surfaces for Gapfill would not occur under the spatio-temporal setting.

Conclusions
The contribution of this article was four-fold: (i) provide an overview of the plethora of methods available for analyzing large spatial datasets, (ii) provide a brief comparison of the methods by implementing a case study competition among research groups, (iii) make available the code to analyze the data to the broader scientific community and (iv) provide an example of the common task framework for future studies to follow when comparing various analytical methods.In terms of comparison, each of the methods performed very well in terms in predictive accuracy suggesting that any of the above methods are well suited to the task of prediction.However, the methods differed in terms of their ability to accurately quantify the uncertainty associated with the predictions.While we saw that some methods did consistently well in both predictive performance and nominal coverage on the simulated and real data, in general we can expect performance of any method to change with size of the dataset, measurement  error variance, and the nature of missingness.However, the data scenario's considered here are relatively representative of a typical spatial analysis such that our results can be used as a guide for practitioners.
At the outset of this study, run time and computation time for each method was of interest.However, because many of these methods are very young in their use and implementation, the variability across run time was too great to be used as a measure to compare the methods.For example, some methods are implemented in R while others are implemented in MATLAB.Still, others use R as a front end to call C-optimized functions.Hence, while we reported the run times in the results section, we provide these as more of an "off the shelf" run time estimate rather than an optimized run time.Until time allows for each method to be further developed and software becomes available comparing run times can be misleading.
Importantly, no effort was made to standardize the time spent on this project by each group.Some groups were able to quickly code up their analysis from existing R or MATLAB libraries.Others, however, had to spend more time writing code specific to this analysis.Undoubtedly, some groups likely spent more time running "in house" cross-validation studies to validate their model predictions prior to the final run on the BYU servers while others did not.Because of this difference, we note that some of the discrepancies in results seen here may be attributable to the amount of effort expended by each group.However, we still feel that the results displayed herein give valuable insight into the strengths and weaknesses of each method.
This study, while thorough, is non-comprehensive in that other methods for large spatial data (e.g.Sang and Huang, 2012;Stein et al., 2013;Kleiber and Nychka, 2015;Castrillon-Candás, Genton and Yokota, 2016;Sun and Stein, 2016;Litvinenko et al., 2017) were not included.Additionally, methods are sure to be developed in the future which are also viable for modeling large spatial data (see Ton et al., 2017;Taylor-Rodriguez et al., 2018).We made attempts to invite as many groups as possible to participate in this case study but, due to time and other constraining factors, not all groups were able to participate.However, in our opinion, the methods compared herein are representative of the most common methods for large spatial data at the time of writing.
We note that the data scenarios considered in this case study do not cover the spectrum of issues related to spatial data.That is, spatial data may exhibit anisotropy, nonstationarity, large and small range spatial dependence as well as various signal-to-noise ratios.Hence, we note that further practical distinctions between these various methods could be made depending on their applicability to these various spatial data scenarios.However, the comparison included here serves as a nice baseline case for method performance.Further research can develop case study competitions for these more complicated scenarios.
Notably, each method was compared only in terms of predictive accuracy.Further comparisons could include estimation of underlying model parameters.The difficulty in comparing estimation, however, is that not all the methods use the same model structure.For example, NNGP uses an exponential covariance while Gapfill does not require a specified covariance structure.Hence, we leave the comparison of the parameter estimates to a future study.
This comparison focused solely on spatial data.Hence, we stress that the results imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 found here are applicable only to the spatial setting.However, spatio-temporal data are often considerably larger and more complex than spatial data.Many of the above methods have extensions to the space time setting (e.g., Gapfill is built directly for spatio-temporal settings).Further research is needed to compare these methods in the spatio-temporal setting.
Let the spatial domain D = D d=1 D d where D 1 , . . ., D D are subregions that form a partition (i.e.D d1 D d2 = ∅ for all d 1 = d 2 ).The modeling approach based on imsart-generic ver.2014/10/16 file: LargeNCaseStudy.texdate: April 26, 2018 spatial partitioning is to assume conditional dependence between observations within a subregion and conditional independence between observations across subregions.More specifically, if Y d = {Y (s i ) : s i ∈ D d } where d = 1, . . ., D, then

FIG 2 .
FIG 2. The top row displays the (a) full and (b) training satellite datasets.The bottom row displays the (c) full and (d) training simulated data.
After the split, the training set for the MODIS data consisted of 105,569 observations leaving 42,740 observations in the test set.The training set for the simulated data also consisted of 105,569 observations but a test set size of 44,431 (the difference in test set size is contributed to missing data due to cloud cover in the original MODIS data).Research teams were provided with the training set and the locations of the test set (but not the actual observation in the test set).
Figure 2 displays the full datasets along with the corresponding training set provided to each research group.All datasets used in this article are provided as supplementary material to this article.

FIG 3 .
FIG 3. Predictions for the simulated data using each of the competing methods.

FIG 4 .
FIG 4. Predictions for the satellite data using each of the competing methods.
X d is a design matrix containing covariates associated with Y d , H d is a matrix of spatial basis functions (such as those used in predictive processes, fixed rank kriging or lattice kriging mentioned above) and Σ(φ d ) is the covariance matrix for subregion d governed by covariance parameters φ d (e.g.decay and smoothness parameters).Notice that, in (2.9) each subregion shares common β and θ parameters which allows smoothing across subregions (hence, Y d1 ⊥ ⊥ Y d2 for d 1 = d 2 conditional on the parameters β and θ).Further, the assumption of independence across subregions allows the likelihood for β, θ and φ d is to be computed in parallel thereby facilitating computation.

TABLE 2
Numerical scoring for each competing method on the simulated data.The best result of each score is bolded.

TABLE 3
Numerical scoring for each competing method on the satellite data.The best result of each score is bolded.