Optimal Spatial Prediction Using Ensemble Machine Learning

1 Introduction

Optimal prediction of a spatially indexed variable is a crucial task in many scientific disciplines. For example, environmental health applications concerning air pollution often involve predicting the spatial distribution of pollutants of interest, and many agricultural studies rely heavily on interpolated maps of various soil properties.

Numerous algorithmic approaches to spatial prediction have been proposed (see Cressie [1] and Schabenberger and Gotway [2] for reviews), but selecting the best approach for a given data set remains a difficult statistical problem. One particularly challenging aspect of spatial prediction is that location is often used as a surrogate for large sets of unmeasured spatially indexed covariates. In such instances, effective prediction algorithms capable of capturing local variation must make strong, mostly untestable assumptions about the underlying spatial structure of the sampled surface and can be prone to overfitting. Ensemble predictors that combine the output of multiple predictors can be a useful approach in these contexts, allowing one to consider multiple aggressive predictors. There have been some recent examples of the use of ensemble approaches in the spatial and spatiotemporal literature. For example, Zaier et al. [3] used ensembles of artificial neural networks to estimate the ice thickness of lakes and Chen and Wang [4] used stacked generalization to combine support vector machines classifying land-cover types in hyperspectral imagery. Ensembling techniques have also been used to make spatially indexed risk maps. For example, Rossi et al. [5] used logistic regression to combine a library of four base learners trained on a subset of the observed data to obtain landslide susceptibility forecasts for the central Umbrian region of Italy. Kleiber et al. [6] have developed a Bayesian model averaging technique for obtaining locally calibrated probabilistic precipitation forecasts by combining output from multiple deterministic models.

The Super Learner prediction algorithm is an ensemble approach that combines a user-supplied library of heterogeneous candidate learners in such a way as to minimize ν-fold cross-validated risk [7]. It is a generalization of the stacking algorithm first introduced by Wolpert [8] within the context of neural networks and later adapted by Breiman [9] to the context of variable subset regression. LeBlanc and Tibshirani [10] discuss stacking and its relationship to the model-mix algorithm of Stone [11] and the predictive sample-reuse method of Geisser [12]. The library on which Super Learner trains can include parametric and nonparametric models as well as mathematical models and other ensemble learners. These learners are then combined in an optimal way in the sense that the Super Learner predictor will perform asymptotically as well as or better than any single prediction algorithm in the library under consideration. Super Learner has been used successfully in nonspatial prediction (see for example Polley et al. [13]). In this paper, we review its optimality properties and discuss the assumptions necessary for these optimality properties to hold within the context of spatial prediction. We also present the results of a simulation study, demonstrating that Super Learner works well in practice under a variety of spatial sampling schemes and data-generating distributions. In addition, we apply Super Learner to a real world dataset, predicting water acidity for a set of 112 lakes in the Southeastern United States. We show Super Learner is a practical, data-driven, theoretically supported way to build an optimal spatial prediction algorithm from a large, heterogeneous set of predictors, protecting against both model misspecification and over-fitting.

2 Problem formulation

Consider a random spatial process indexed by location over a fixed, continuous, d-dimensional domain, {Y(s): s∈D⊂Rd}. For a particular set of distinct sampling points {S1,…,Sn}⊂D, We observe {(Si,Yi∗): i=1,…,n}, where Y∗=Y(Si)+ϵi and ϵi represents measurement error for the i^th observation. For all i, we assume E[Yi∗|Si=s]=Y(s). Our objective is to predict Y(s′) for unobserved locations s′⊂D. Thus, our parameter of interest is the spatial process itself. We do not make any assumptions about the functional form of the spatial process. We do, however, assume that one of the following is true: for all i, either

1. (Si,Yi∗) are independently and identically distributed (i.i.d.), or

2. (Si,Yi∗) are independent but not identically distributed, or

3. Yi∗ are independent given S₁, …, S_n; and E[Yi∗|S1,…,Sn]=E[Yi∗|Si]=Y(Si). This corresponds to a fixed design.

While these are reasonable assumptions for many spatial prediction problems, they are nontrivial and may not always be appropriate. For instance, instrumentation and calibration error within sensor networks can result in spatially structured measurement error that is not mean zero given S1,…,Sn. There has been an effort on the part of researchers to develop ways to adapt the cross-validation procedure so as to minimize the effects of this kind of measurement error when choosing parameters such as bandwidth in local linear regression or smoothing parameters for splines. Interested readers should consult Opsomer et al. [14] and Francisco-Fernandez and Opsomer [15] for overviews.

3 The Super Learner algorithm

Suppose we have observed n copies of the random variable O with true data-generating distribution P0∈M, where the statistical model M contains all possible data generating distributions for O. The empirical distribution for our sample is denoted P_n. Define a parameter Ψ: M→Ψ≡{Ψ[P]: P∈M} in terms of a risk function R as follows: ΨP=argminψ∈ΨRψ,P In this paper, we will limit our discussion to so-called linear risk functions, where R(ψ,P)=PL(ψ)=∫L(ψ)(o)dP(o) for some loss function L. For a discussion of nonlinear risk functions, see van der Laan and Dudoit [16]. We write our parameter of interest as ψ0=Ψ[P0]=argminψR(ψ,P0), a function of the true data generating distribution P₀. For many spatial prediction applications, the Mean-Squared Error (MSE) is an appropriate choice for the risk function R, but this needn’t necessarily be the case.

Define a library of J base learners of the parameter of interest ψ0, denoted {Ψˆj: Pn→Ψˆj[Pn]}j=1J. We make no restrictions on the functional form of the base learners. For example, within the context of spatial prediction, a library could consist of various Kriging and smoothing splines algorithms, Bayesian hierarchical models, mathematical models, machine learning algorithms, and other ensemble algorithms. We make a minimal assumption about the size of the library: it must be at most polynomial in sample size. Given this library of base learners, we consider a family of combining algorithms {Ψˆα=f({Ψˆj: j},α):α} indexed by a Euclidean vector α for some function f. One possible choice of combining family is the family of linear combinations, Ψˆα=∑j=1Jα(j)Ψˆj. If it is known that ψ0∈[0,1], one might instead consider the logistic family, log[Ψˆα/Ψˆα(1−Ψˆα)(1−Ψˆα)]=∑j=1Jα(j)log[Ψˆα/Ψˆα(1−Ψˆα)(1−Ψˆα)]. In either of these families, one can also constrain the values α can take. In this paper, we constrain ourselves to convex combinations, i. e. for all j,α(j)≥0 and ∑jα(j)=1.

Let {B_n}be a collection of length n binary vectors that define a random partition of the observed data into a training set {Oi: Bn(i)=0} and a validation set {Oi: Bn(i)=1}. The empirical probability distributions for the training and validation sets are denoted Pn,Bn0 and Pn,Bn1, respectively. The estimated risk of a particular estimator Ψˆ: Pn→Ψˆ[Pn] obtained via cross-validation is defined as

Given a particular class of candidate estimators indexed by α the cross-validation selector selects the candidate which minimizes the cross-validated risk under the empirical distribution P_n,

The Super Learner estimate of ψ0 is denoted Ψˆαn[Pn].

4 Cross-validation and spatial data

The theoretical results outlined above depend on the training and validation sets being independent. When this is not the case, there are generally no developed theoretical guarantees of the asymptotic performance of any cross-validation procedure [18]. Bernstein’s inequality, which van der Laan and Dudoit [16] use in developing their proof of the oracle inequality, has been extended to accommodate certain weak dependency structures, so it may be that there are ways to justify certain optimality properties of ν-fold cross-validation in these cases. There have also been some extensions to potentially useful fundamental theorems that accommodate other specific dependency structures. Lumley [19] proved an empirical process limit theorem for sparsely correlated data which can be extended to the multidimensional case. Jiang [20] provided probability bounds for uniform deviations in data with certain kinds of exponentially decaying one-dimensional dependence, although it is unclear how to extend these results to multidimensional dependency structures where sampling may be irregular. Neither of these extensions is immediately applicable to the general spatial case, where sampling may or may not be regular and the extent of spatial correlation cannot necessarily be assumed to be sparse. There has been some attention in the spatial literature to the use of cross-validation within the context of Kriging and selecting the best estimates for the parameters in a covariance function, most of it urging cautious and exploratory use [1, 21]. Todini [22] has investigated methods to provide accurate estimates of model-based Kriging error when the covariance structure has been selected via leave-one-out cross-validation, although this remains an open problem.

Recall from Section 2 above that our parameter of interest is the spatial process Y (s) and we have assumed E[Y∗|S=s]=Y(s). Even if Y (s) is a spatially dependent stochastic process such as a Gaussian random field, the true parameter of interest in most cases is not the full stochastic process, but rather the particular realization from which we have sampled. Conditioning on this realization removes all randomness associated with the stochastic process, and any remaining randomness comes from the sampling design and measurement error. So long as the data conform to one of the statistical models outlined above in Section 2, the optimality properties outlined above will apply.

5 Simulation study

We applied the Super Learner prediction algorithm to six data sets with known data generating distributions simulated on a grid of 128×128=16,384 points in [0,1]2⊂R2. Each spatial process was simulated once, hence samples of stochastic processes were taken from a common realization. All simulated processes were scaled to [−4,4] before sampling. See Figure 1 for an illustration.

The function f1(⋅) is a mean zero stationary Gaussian random field (GRF) with Matérn covariance function [23]

where h is a distance magnitude between two spatial locations, σ2 is a scaling parameter, ϕ>0 is a range parameter influencing the spatial extent of the covariance function and τ2 is a parameter capturing micro-scale variation and/or measurement error. Kν(⋅) is a modified Bessel function of the third order and ν>0 parametrizes the smoothness of the spatial covariation. Learners were given spatial location as covariates.

The function f2(⋅) is a smooth sinusoidal surface used as a test function in both Huang and Chen [24] and Gu [25], f2(s)=1+3sin(2π[s1−s2]−π). Learners were given spatial location as covariates.

The function f3(⋅) is a weighted nonlinear function of a spatiotemporal cyclone GRF and an exponential decay function of distances to a set of randomly chosen points in [−0.5,1.5]2⊂R2. In addition to spatial location, learners were given the distance to the nearest point as a covariate.

Figure 1:

The six spatial processes used in the simulation study. All surfaces were simulated once on the domain [0, 1]². Process values for all surfaces were scaled to [−4,4]⊂R.

The function f4(⋅) is defined by the piecewise function

where w is Beta distributed with non-centrality parameter 3 and shape parameters 4 and 1.5. Learners were given spatial location and w as covariates.

The function f5(⋅) is a sum of several surfaces on [0,1]⊂R2; a nonlinear function of a random partition of [0, 1]²; a piecewise smooth function; and w2∼uniform(−1,1). Learners were given spatial location, partition membership (w₁) and w₂ as covariates.

The function f6(⋅) is a weighted sum of a spatiotemporal GRF with five time-points, a distance decay function of a random set of points in [0, 1]², and a beta-distributed random variable with non-centrality parameter 0 and shape parameters both equal to 0.5. Learners were given spatial location, the five GRFs and the beta-distributed random variable as covariates.

6 Practical example: predicting lake acidity

We applied Super Learner to a lake acidity data set previously analyzed by Gu [25] and Huang and Chen [24]. Increases in water acidity are known to have a deleterious effect on lake ecology. Having an accurate estimate of the spatial distribution of lake acidity is an essential first step toward crafting effective regulatory interventions to control it. The data were sampled by the U.S. Environmental Protection Agency during the Fall of 1984 in the Blue Ridge region of the Southeastern United States [47], and consist of longitudes and latitudes (in degrees), calcium ion concentrations (in milligrams per liter), and pH values. The EPA used a systematic stratified sampling design which we treated as fixed here. Because only one sample per lake was collected, we assume some measurement error that is independent of lake pH, calcium ion concentration, and spatial location. The data are freely available in the R package gss [48]. We used the same nearly equal area projection as Gu [25] and Huang and Chen [24],

where xˉlat and xˉlon are the midpoints of the latitude and longitude ranges, respectively.

Let xi=(xi,1,xi,2) denote the i^th sampling location; w_i denote the calcium ion concentration observed at the i^th sampling location; and Yi∗ be the pH value observed at the i^th sampling location. We assume that E[Yi∗|Si=s]=Y(s), where Si=(xi,wi). Our objective is to learn the lake pH spatial process from the data.

The library used to predict lake acidity was similar in composition to the simulation library described in Subsection 5.1, with some important differences. We reduced the number of parameterizations for some of the algorithm classes in the library. We used one DSA learner, which used 10-fold cross-validation and considered polynomials of up to five terms (m=5), each term being at most a two-way interaction (k=2) with a maximum sum of powers p=3. We used a reduced number of parameterizations of GAM, GBM, TPS, GP, and SVM learners, as well. We also included screening algorithms that allowed us to train learners on specific subsets of covariates: x, w, log w, (x, w), and (x, log w). We considered the L₂ loss function, and the predictions from all base learners were truncated to the observed pH range in order to ensure a uniformly bounded loss.

Table 6 in Appendix B provides a detailed list of the library and shows performance results for each base learner as well as Super Learner itself. Figure 2 provides graphical representations of Super Learner’s pH predictions. The Super Learner predictions are a slightly smoothed version of the observed data, with attenuated predictions for the highest and lowest observations. We also cross-validated Super Learner (V=5) to get an unbiased estimate of its risk. Super Learner had nearly identical performance in terms of cross-validated risk to GAM (4 degrees, using × and log w) and both versions of Random Forest. The GBM predictors and OK (using log w) had slightly smaller cross-validated risks.

Many of the algorithms in the library performed better when given log w as opposed to w, but for those algorithms like GBM and Random Forest that were not attempting to fit some kind of polynomial trend, logging the calcium ion concentration made little difference in performance. This is not surprising, given the relationship between the activities of Ca2+ and H+ As expected, most algorithms had cross-validated risk estimates that were worse than their empirical risk estimates calculated from predictions made after training on the full data set. The Kriging algorithms, for instance, were all exact interpolators when trained on the full data, and thus had estimated empirical MSEs of 0, whereas their MSEs estimated via cross-validation ranged from 0.07 (FVU=0.46) to 0.11 (FVU=0.72). The Gaussian processes with RBF kernel had the most pronounced differences between the two risk estimates. For example, GP (RBF) trained on the covariates (x, w) had an empirical MSE of 0.01 (FVU=0.08) and a cross-validated MSE of 0.22 (FVU=1.46).

One should use caution in general when interpreting the weights returned by the Super Learner algorithm. It can be tempting to use them to try to draw conclusions about the underlying data generating process. However, these weights are not necessarily stable, and we would expect them to differ given a different sample, or a different partitioning during cross-validation. Performance guarantees are with respect to the ensemble, and not the ranking of individual algorithms whose predictions are being averaged.

Figure 2:

(a) A map of Super Learner’s pH predictions, and (b) a plot of Super Learner’s predictions as a function of the observed data. Super Learner mildly attenuated the pH values at either end of the range, but otherwise provided a fairly close fit to the data.

7 Conclusion/discussion

In this article, we have demonstrated the use of an ensemble learner for spatial prediction that uses cross-validation to optimally combine the predictions from multiple, heterogeneous base learners. We have reviewed important theoretical results giving performance bounds that imply Super Learner will perform asymptotically at least as well as the best candidate in the library. We discussed the assumptions required for these optimality properties hold. These assumptions are reasonable for many measurement error scenarios and commonly implemented spatial sampling designs, including various forms of stratified and random regular sampling.

Our simulations and practical data analysis used comparatively modest sample sizes. With the increased availability of massive environmental sensor networks and extremely large remotely sensed imagery, scalability of a spatial prediction method is often an important practical concern. The cross-validation step of Super Learner is a so-called ’embarrassingly parallel’ problem, and there are implementations of Super Learner that parallelize this step (see for instance the function CV.SuperLearner in the SuperLearner R package). There are also a number of versions of the Super Learner algorithm under active development that are intended specifically to cope with very large data sets. One is called h2oEnsemble [49], and is a stand-alone R package that makes use of a Java-based machine learning platform that can be used for massive cluster computing. Another version of Super Learner has been modified to work with streaming data [50].

In this paper, we have not addressed dependent sampling designs, where sampling at one point changes the probability of sampling at another point. This is an important area for future research. We also limited our scope to the case where measurement error is at least conditionally mean-zero. Spatially structured measurement error that is not conditionally mean zero is a common problem in many spatial prediction applications, and there have been a number of attempts to alter the cross-validation procedure to accommodate it [51, 15]. These proposed techniques generally require one to estimate the error correlation structure from the data or to know it a priori. How well these algorithms perform if the correlation extent is substantially underestimated is unknown. Ideally, it would be best to have a stronger theoretical understanding of how the degree of dependence between training and validation sets affects cross-validated risk estimates both asymptotically and in finite samples. This is an important future area for research.