Cross-validation for Geospatial Data: Estimating Generalization Performance in Geostatistical Problems
Authors: Jing Wang
Laurel Hopkins
Tyler Hallman
W. Douglas Robinson
Rebecca Hutchinson
Venue: Accepted by TMLR ()

Abstract: Geostatistical learning problems are frequently characterized by spatial autocorrelation in the input features and/or the potential for covariate shift at test time. These realities violate the classical assumption of independent, identically distributed data, upon which most cross-validation algorithms rely in order to estimate the generalization performance of a model. In this paper, we present a theoretical criterion for unbiased cross-validation estimators in the geospatial setting. We also introduce a new cross-validation algorithm to evaluate models, inspired by the challenges of geospatial problems. We apply a framework for categorizing problems into different types of geospatial scenarios to help practitioners select an appropriate cross-validation strategy. Our empirical analyses compare cross-validation algorithms on both simulated and several real datasets to develop recommendations for a variety of geospatial settings. This paper aims to draw attention to some challenges that arise in model evaluation for geospatial problems and to provide guidance for users.

Published in Transactions on Machine Learning Research (09/2023)
Cross-validation for Geospatial Data: Estimating
Generalization Performance in Geostatistical Problems
Jing Wang
wangii@@onggonsttteeedt
School of Electrical Engineering and Computer Science
Oregon State University
Corvellis, OR 97331-5501, USA
Laurel M. Hopkins
hopkilaax@oregonsaaaae edu
School of Electrical Engineering and Computer Science
Oregon State University
Corvallis, OR 97331-5501, USA
Tyler A. Hallman
f. t.halman hattimaanraaaaggr.ccco
School of Natural Sciences Bongor University Bongor LL57 2DG, UK
W. Douglas Robinson Department of Pisheries, Wildlife, and Conservation Sciences Oregon State University Corvallis, OR 97331-5501, USA
douglas, www.rrooinnnooooooossstieoe edu
Rebecca A. Hutchinson
ral mah@oregonssaae.com edu
School of Electricul Engineering and Computer Science Department Pisheries, Wildlife, and Conservation Sciences Oregon State University Corvallis, OR 97331-5501, USA
Reviewed on OpenReview: https: o/pp::////eeeeveee// net/foruse YyJ
Abstract
Geostatistical learning problems are frequently characterized by spatial autocorrelation in the input features and/or the potential for covariate shift at test time. These realities violate the classical assumption of independent, identically distributed data, upon which most cross-validation algorithuns rely order estimate the generalization performance of a model. In this paper, we present a theoretical criterion for unbiased cross-validation estimators in the geospatial setting. We also introduce new cross-validation algorithm to evaluate models, inspired by the challenges of geospatial problems. We apply framework for categorizing problems into different types of geospatial scenarios to help practitioners select an appropriate cross-validation strategy. Our empirical analyses compare cross-validation algorithis on both simulated and several real datasets to develop recommendations for a variety of geospatial settings. This paper aims to draw attention to some challenges that arisc in model evaluation for geospatial problems and to provide guidance for users.
1 Introduction
A variety of goospatial problems benefit from machine learning (ML) analysos. Of note, many problems with socictal and environmental importance have geospatial properties. For example, machine learning models support crop production in agricultural applications (Dadsetan al. 2021 predictions of poverty

---

Published in Transactions on Machine Learning Research (09/2023)
(Jean et al. 2016), analysis of forest fire risk (Yang et al. 2021), and creation of land cover maps (Bjorck
et al. 2021 across space. As specific motivating example, consider the problem of species distribution modeling (SDM). Species distribution models capture relationships between species and their habitats based on observations at a set of locations and the associated environmental variables at those locations (Elith & Leathwick 2009]. SDMs are critical tools for ecologists and natural resource managers seeking to model the current, potential, and future distributions of threatened species (e.g.. |Wilson et al., 2021 Hallman et all 2021). Machine learning approaches like random forests, boosted regression trees, and maximum entropy models have been shown to perform well on this task (Valavi et al. 2021); these methods are flexible and powerful but not explicitly spatial models. They are typically fit with only environmental features and not spatial coordinates as imputs. Some key common characteristics of the geospatial problems we focus on in this paper are: 1) they fit models to and predict the response variable al geolocated data points; 2) despite the spatial nature of the data, practitioners may choose to apply models that do not explicitly account for space, instcad favoring casy-to-use, flexible, non-spatial techniques; 3) the intended test region and its environmental features may be known or estimable.
Just in non-spatial domains, geospatial applications rely on estimates of generalization performance for model evaluation and selection. Classical vross-validation (CV) techniques estimate generalization perfor- mance by dividing the training set into folds; by default, they assign data points to folds uniformly at random. Each fold takes a turn acting as a validation fold on which performance is measured after training a model on the other folds. These measurements estimate how well a model will perform when applied to a new, unseen dataset (i.e., the test dataset). These performance estimates may have stand-alone value; e.g., natural resource manager may only wish to proceed with some conservation action based a sufficiently high-quality SDM. Performance estimates may also serve model selection goals; c.g., ccologist may want to select the best of set of modeling approaches for further study.
Geospatial problems pose unique challenges to standard CV approaches. In particular, theoretical results showing the unbiasedness of the standard CV estimators assume independent, identically distributed (iid) data (Arlot & Celisse 2010 Hollimann et al. 2021). However, spalial aulocorrelation among features (i.e., feature values at geographically proximal points are more similar than feature values at geographically distant points) in geospatial contexts induces a correlation structure among CV folds, data points in the training and validation folds may violate the assumption of independence. Furthermore, covariate shift is common in geospatial problems (i.e,, the distribution of features may change across space), so data points in the training and validation folds may violate the assumption of being identically distributed. In a given problem, either or both of these challenges may be present to varying degrecs. For example, in one SDM application in Sec. 7.3.1 below, we applied CV to models built with a training set comprised of bird observations at two national parks to estimate performance on a third, geographically distinct, national park. Prediction points in the test park are not spatially autocorrelated with those in the training parks, though the points within the training parks are spatially autocorrelated with each other. In addition, the feature distributions at the test park differ from that of the training parks. Therefore, this problem does not meet the iid assumption required to guarantee unbiased performance estimates from standard CV techniques. Some extensions to these techniques have been developed lo address spatial autocorrelation and covariale shift (reviewed below), but an approach to address both simultaneously is still lacking (Hoffimann et al. 2021).
In this paper, we report on an investigation into how spatial autocorrelation and covariate shift combine in different geospatial learning settings to affect generalization performance estimates from a variely of CV strategies. Our specific contributions are:
We prove theoretical criterion for unbiasedness of CV estimators in set with spatial autocor- relation and/or covariate shift.
We explore criteria for sorting problems into geospatial scenarios to aid in selection of an appropriate CV approach for a given analysis.
We propose a new CV algorithm to address spatial autocorrelation and covariate shift simultaneously.
We provide simulated data experiments measuring the bias of several CV estimators, SDM examples demonstrating four geospatial scenarios, and further empirical analyses of our proposed algorithm.
2

---

Published in Transactions on Machine Learning Research (09/2023)
The rest of the paper is organized follows. Section 2 gives background on cross-validation and reviews related work. Section B defines the geospatial problem setting formally. Section proves criterion for unbiased CV estimators in geospatial settings. Section 5 explores a framework for categorizing geospatial problems into four discrete scenarios, based on whether they display spatial autocorrclation and/or covariate shift. Section introduces a new cross- vooss-vaiiittion procedure aimed at geospatial settings in which both spatial autocorrelation and covariate shift are at play. Section reports on three sets of experiments evaluate set of CV procedures. Section g
G
discusses the insights and recommendations we can glean from this work far. Section conchudes with summary and some ideas for future research.
2 Background & Related Work
First, let be clear about how we quantify generalization performance. Test error and risk are two common statistics in model evaluation (Hastic et al. 2001). Test error (Errr) is the expected loss over test samples, given a training set T of n examples. Risk (Rt) is the expected test error over training sets from the same population. The critical difference between these is whether an expectation is taken over training sets (risk), as opposed to conditioning on a single training set (test error). Hereafter, tr and te subscripts denote training and test, respectively, and i and subscripts index training and test samples, respectively. We use X and x for features and y and y for response variables. Then for general loss function C, lest error and risk are defined as
Err.") (u) = Er, [Cys ((x;; T)] R( ET,T, [C( [xx;; T))]
(1)
where 7; = (x), 353 is test sample randomly drawn from the joint distribution Px T {xi, yi Hi is training set randomly drawn from the joint distribution Px...err and (xx;: is the prediction of yj given Xj and T.
Cross-validation assesses a model's predictive capabilities. Its target is sometimes misunderstood be an estimate of test error, but under the iid assumption, actually estimates the risk (Bates et al., 2023). Hastie et al. 2001 conclude that the estimation of test error for particular training set is not easy in general, but that validatai may provide rcasonable estimates of risk. Both risk and test error may be of interest to practitioners, in lerms of expected performance on testing sets drawn from some population. Risk estimates speak to a model's expected performance when fit to training sets drawn from some population, whereas test error estimates speak to a model's expected performance when fit to a particular training set. Risk is amenable to theoretical analysis and sinulation experiments, and this the quantity of focus for those sections below. However, when only a single training set available, as in our empirical analyses, only test error is estimable. The discrepancies between the risk and test error perspectives are challenging aspect of research cross-validation.
Standard techniques for cross-validation put data points into folds uniformly at random. K-fold cross- validation (KFCV) randomly divides a training set T into K non-overlapping folds and iteratively holds out one fold at a time, training a model on the remainder (training folds T-x) and measuring error on the held-out fold (validation fold Tk). The average of these model errors across folds is the estimate of risk. Leave-onc-out vross-valiittion (LOOCV) is n-fold CV, in which the validation fold is single data point, Ti, and the training fold contains the remaining data, T-.. The corresponding risk estimates for these approaches are
A(n) 1
RKFC - -EK Eien, Cly, (x): 71
(()) Roocv = - En L(y,, (x:
The dias-variance trade-off in KFCV controlled by the user-defined parameter K (Alot & Celisse 2010 In general, bias decrcases while variance increases as K increascs. To see this, note that LOOCV builds n models, each of which is trained on n 1 samples, while KFCV builds K models, each of which is trained (K - 1)n/ samples. Since the size of the training folds under LOOCV is closer to the full dataset size
3

---

Published in Transactions on Machine Learning Research (09/2023)
than in KFCV, LOOCV incurs less bias. However, LOOCV is more computationally expensive than KFCV when n > K. For iid data, these standard CV estimators are unbiased; their random resampling mechanism mimics the way that a new sample would be collected from the population.
In addition to these standard approaches, a few spatial cross-validation procedures have been developed to deal with non-independent data. The issue motivating these techniques that when nearby points are placed in different CV folds, spatial autocorrelation among the features at those points makes those folds dependent. This violates the assumptions of the methods, and it could optimistically bias generalization performance estimates if information crosses fold boundaries. Inspired by h-block CV (Burman et al., 1994), there are two general families of spatial CV approachos. Block cross-validation (BICV) groups geographically proximal points into the same block (RuB & Brenning, 2010 Roberts et al. 2017). Buffered cross-validation (BFCV) places a buffer between training and validation folds and removes points within the buffer so that they are used in neither set (Le Rest et al. 2014; PPohjankukka et al.] 2017/. Block size and buffer size are user- defined hyperparameters. One suggested practice is to set them to the modian of the spatial autocorrclation ranges of the features. Spatial CV is implemented in some software packages (Brenning) 2012 Muscarella et al. 2014 Valavi et al. 2018) and has been used in previous studies (Zurell et al. 2020 Seo et al. 2021
Valavi et al. 2023) to increase the average distance between training and validation points. One potential issue that may arise with spatial CV techniques is the introduction of covariate shift between folds. By Tobler's first law of geography (Tobler) 1970). things that are close together are more similar than things that are far apart. When creating spatially separated cross-validation folds, it is possible to construct folds with differing feature spaces, and in some cases, this could result in pessimistically biased error estimates. Outside of spatial dependence structure, prior work has also considered the impacts of correlation structure cross-validation through the lens of latent random effects rather than autocorrelation (Rabinowiez & Rosset 2020).
There is also prior work on CV procedures that do not assume that the training and test sets are identically distributed. Covariate shift ariscs frequently in goospatial learning problems when the training feature distribution is different from the feature distribution under which predictions are desired. For example, when considering whether to translocate a threatened species to a currently unoccupied region, natural rosource managers may need to construct habitat suitability in the new area. Even when predicting at points that interpolate spatially between the training locations, recall that spatial CV approaches for breaking up spatial autocorrelation may actually induce covariate shift. CV algorithms to account for covariate shift have been developed in non-spatial settings, most notably including ampprrancc-weightaee IWCV assumes that the training and test. data are independent but from different distributions. It develops an asymptotically unbiased estimator by rectifying the loss function with the ratio between the test and training probability densities. IWCV works well when the support of the test distribution is within that of the the training distribution but poorly with minimal overlap between the two, and it does not consider dependent samples. Meyer & Pebesma
SDM with data from a currently occupied region and predict
vross-validation (IWCV, Sugiyama et al. (2007)).
(2021) suggest an alternative for dealing with covariate shift, by using a dissimilarity index to define an appropriate area of applicability, outside of which cross-validation estimates of model performance do not apply.
anportance-weigheed CV procedures (both in prior work and our proposed approach) require density ratio estimates. Density ratios can be estimated efficiently by variety of algorithms, among which Relative racoedtrained Least-Squares Importance Fitting (RuLSIF, Yamada et al. 22011) outstanding in theory and practice. Its predecessor, unconstrained Least-Squares Importance Fitting (uLSIF, Kanamori et al. (2009) provides a closed-form solution by solving a system of lincar equations and thus solution is always numerical stable. In importance estimation and covariate shift adaptation tasks, uLSIF is empirically more accurate and faster than other competing importance estimation algorithms (Kanamori et al. (2009)). RuL.IF shares several advantages with uLSIF, like an analytical solution, numerical stability, and robustness, but RuLSIF achieves even better non- por-maraett convergence. For these reasons, we estimated density ratios with RuLSIF in all experiments below.
When users are unaware of the nuances and applicability of these methods, both standard and spatial CV algorithms may be misused. A few recent studies have compared the CV algorithms above. Roberts et al. 22117) Meyer et al. (2019). and Ploton et al. (2020) conclude that spatial CV is less biased than
4

---

Published in Transactions on Machine Learning Research (09/2023)
non-spatial CV. However, Hoffimann et al. (2021) show that blocking cross-validation overestimate error. Wadoux et al. (2021) consider sampling probability and conclude that for clustered data KFCV underestimates error but spatial CV algorithms severely overestimate it. Such contradictory conclusions may stem from a lack of clarity on the applicability of a variety CV strategics. This paper seeks to further the discussion by examining their utility across different geospatial scenarios. Mill et al. (2022) specifically considered geographical interpolation and extrapolation and developed the Nearest Neighbour Distance Matching (NNDM) LOOCV, and compared it with no-spatial LOOCV and buffered LOOCV. While we investigate various CV methods under the K-fold setting which is more common in practice.
Finally, by way of background, we remind readers of the concept of scmivariogram range (or variogram range) to decide spatial dependence and some vocabulary that will be used below. Semivariogram range is a widely used concept to determine spatial dependence in geostatistics. It is defined as the lag distance where semivariance levels off at its maximum, indicating that pairs of observations at least as far apart as the range are no longer spatially dependent (i.c., they become spatially independent) (Cressic 2015).
3 The Geospatial Problem Setting
Goospatial problems pertain to data with both spatial coordinates (locations) and variables measured those locations (features); here, we formalize this notion. Our notation draws on that used by Hoffimann ct al. 2021]. A spatial random variable is a stochastic process Z : D x Ω - R, where D C Rd is a spatial domain, and Ω is a sample space. For example, an environmental feature vector x; = (2r
(1) .......""! (m)) is
collection of processes (e-g., temperature, precipitation) a specifie location E D. Once we obtain feature values of X1 from an observation, a realization of the processes is produced by fixing w E Ω. In gcostatistical learning, we seek to predict the response variable in the test set Yte = {uy 1== given training data {Xer: Yurr = {x1, Mi collected from locations in spatial domain Drr and test. set input features Nte = [ collected from locations in spatial domain Dec- Assuming the underlying spatial processes are stationary, {xr are from the distribution Px., and are spatially autocorrelated. We use P and p for distributions and probability density functions, respectively. Similarly, } are from Px.. and are spatially autocorrelated. If training and test data have the same spatial-temporal extent, then Dtr and Dte are considered the same; otherwise, they are different.
Geospatial problems may exhibit different kinds of data shift, which happens when the joint distribution of features and responses differs between the training and test sets (Quinonero-Canded Candela et al. 2009). Covariate shift, one type of data shift, occurs when Px.. * Px.. while Py..ix.r = Py..IX..- Concept shift, the most intractable form of data shift, happens when the relationship of features to response differs across training and test phases: Px # l.,:
; we do not consider it in this paper. Another type of data shift is possible but has not been formally defined before: we refer to the case where the training and test sets have different variogram ranges range shift. When Dir and Dte are different, the problem is likely to involve covariate shift and/or range shift.
4 A Criterion for Unbiased CV
Rabinowicz & Rosset (2020) prove that when the features of the training and test sets are iid and a latent variable induces correlation structure in the response variable, the CV risk estimate (under squared error loss) is unbiased when the joint distribution of test and training sets is the same as the joint distribution of validation and training folds, for all folds. To extend this result, we consider autocorrelation as the mechanism creating dependence, with potential covariate shift, and general loss function. With the geospatial settings defined in Sec. 3 we further assume the training and test features have the same domain: X1, E 2' C R" and the function f f + is unchanged between training and test sets (i.c., no concept shift). Therefore, the domains of the response variables are also the same: Di, E V C RR. Notations with ir and le subscripts (e.g., Tor and Tie) demonstrate different prediction goals, rather than actual datasets. Specifically, Ter and Th denote training data and test data, respoctively; but does not imply that T is the whole datasct. As defined in See. 2, T is the training set which is further divided into validation fold T3 and training fold T With these assumptions in place, we can state the main theorem.
5

---

Published in Transactions on Machine Learning Research (09/2023)
Theorem 1. If Xe Px,Ix at Ak E I,...,K. K, then cross-volidation is asymptotically unbiased estimator of the risk R(a).
Proof. Since p X(rt
we have
p(Xte =
p(X,, Xx))/((X))))
(3)
This restatement of the antecedent of the theorem applics the features X, but to evaluate the bias of the risk estimator, we need lo work with the full datasets (Tar, Tae), which include the response variables y- To do this, multiply the LHS by p(yte|XXte and the RHS by p(yx|Xa)) Since is assumed constant, these quantities are equal (i.c., p(Yte (Xe) = p(ytr/X,r)) P(yx/Xx)).
We focus first on the LHS. Multiplying by p(ytal (10) gives
p(Xee X p(yce Xte)
p(Xor)
Since y conditionally independent of all other variables given its corresponding X, we can condition on additional variables (i.e., p(ycelXte) p(ycelXie Xtr)). Replacing p(ycelXte) with p(ycelXte: Xtr) gives
p(Xes Xar) pler Xar)
p(X+r)
Next, multiply the numerator and denominator by p(yer Xtr), except that by the same reasoning above, p(yte(Xirr) = p(yo) Xtr, Xte: Bhe), so we write the multiplier as plyteltee P(y
This gives
p(Xtes Xer) p(yrol Xte: Xec) -p(yer (ytor Xte,Xer)
(((rr) p(y)r(X(rr
plyte, Xte, Xtx) ppytrl/tao Xeer Xtr)
plytr: X(r)
plyter Xtes Ste.
p(ytr; Xer)
p([Tee Ter)
p(Tir)
This quantity, which now includes the response variables as well as features, is p(Ti
The same reasoning applics as we multiply the RHS of Eqn. by plyx Xx).. The condensed process is:
p(xx, X X- p(yx|(a)
p(X-+)
p(X,, x-+) x) p(yl) XR. X b) ply - k/yk, X-k)
p(X-p 2)
plyk. X X p(y -k Yk. X X_k)
ply k. X-a)
plyR, X k. x-x)
ply x-))
((T),T_K)
=P(TN/T_A).

---

Published in Transactions on Machine Learning Research (09/2023)
Combining the manipulated LHS and RHS, we conclude that
(4)
Now we are ready to show unbiasedness for the LOO setting, For this setting, Eqn. 4 is written as p[T) = p(T;|T-.), recalling that 7j is a single intended test instance outside of the full training set T, T; is a single-instance fold, and T_ is the data from the remaining folds. As is typical, we assume that T-, is distributed as and of size n, ignoring the bias from the different sizes of T-; and Т; this gives p(T,) (T) I(T) (TT[[T_),, which is needed for step (1) below. We use shorthand L for Clys, (x;;;:T))) and Lj for Cyy, U3 (x;:
By [R
Loocvl
Σ x==
[Ca]
[ / p(x), y vi(T C.
Ei ET
n1
(1) 1 Σ1 ET [1 J p(xs, Ms (T= C
n
Ei= ET_,2; [c,]
El R(a-1) "33 R((n)
n
All of these claims also hold for KFCV, with more bookkeeping required to account for varying fold sizes.
Thm. D states that to achieve unbiasedness, the distributional relationship between training and lest features must be preserved by the split. into cross-validation folds. To see the intuition, consider the example of a multivariate Gaussian distribution over Xte and Xtr, and one over X₃ and X-R- To meet the condition of Thm. [ the means, variances, and covariances of the distributions should match. In particular, the off- diagonal elements of the covariance matrices should be the same (i.e., cov(Xte, Xxr) = cou(X X-x)) This means that if there is spatial (in)dependence between the test and training sets' features, then there should be spatial (in)dependence between the validation and training folds' features. In addition, note that the mean and variances of X,e and X2 should match for all folds.
5 Framework for Selecting CV Method
In this section, we present framework for assessing the relevant characteristics of geospatial problems and linking them to the choice of CV procedure. The framework uses two criteria to assess whether spatial dependence and/or covariate shift are at play, which sorts the problem into one of four scenarios (Tab. 1. This discrete categorization facilitates our exploration in this paper of which CV strategy is appropriate for each scenario. However, both properties truly exist on spectra rather than being simply present versus absent, and it may be useful in future work to explore other characterizations.
5.1 Characteristics of Spatial Problems
First, we characterize the spatial antocorrelation between training and testing sets. Using idcas from spatial statistics, we compare the nearest. distance (d) between training and test samples, and the semivariogram range (r) of the training features. If d 2 r, treat the training and test sets displaying spatially independence (right column of Tab. 1 otherwise, they are spatially dependent.
Then, we apply the Cramér-von Mises two-sample test (abbr. Cramér test; Anderson (1662) to assess covariate shift between training and test features. The Cramér test is a multivariate, distribution free test with a null hypothesis of no covariate shift; i.c., Ho: two samples (c.g., X,, and XGa) are distributod identically. When the p-value (p) is smaller than a predefined significance level a, we reject Hy and treat covariate shift as present (bottom row of Tab. [] herwiee we treat it as absent.
7

---

Published in Transactions on Machine Learning Research (09/202))
Table 1: Goostatistical scenarios determined by semivariogram range and the Cramér test. We consider
spatial (in)dependence between training and testing sets comparing the nearest distance (d) bet
training and test samples with the semivariogram range (r) of the training features; and covariate shift by comparing the p-value (p) of Cramér test and the user defined significance level
Spatial Dependence
Spatial Independence
Scenario SD:
Scenario SI: d2 r and p 2
No Covariate Shift
d <<r and p a
Scenario SD CS:
Scenario SI + CS:
Covariate Shift
d r and < a
d 2 r and p < a
Table 2: Summary of how cross-validattin methods approach spatial aulocorrelation and covariate shift. The first four methods are discussed above; the last method (IBCV) is proposed in the next section. Note that for BLCV and BFCV, the extent to which spatial dependence is ameliorated and covariate shift introduced is moderated by tuning parameters (block or buffer size) and the pattern of spatial autocorrelation.
Method
Approach to Spatial Auto- correlation
Approach to Covariate Shift (CS)
Intended
Scenario
KFCV
Preserves spatial relationships of training set when creating folds
None
SD
(Standard K-Fold)
IWCV
Preserves spatial relationships of training set when creating folds
Density ratio estimates
SD + CS
(Importance- Weighted)
BLCV
Reduces spatial dependence
May introduce CS among folds
(BLocking)
across folds
BFCV
Removes spatial dependence
May introduce CS among SI folds Density ratio estimates
(BuFfered)
across folds
IBCV
Removes spatial dependence
ST CS
(Importance- Weighted
across folds
Bullered)

---

Published in Transactions on Machine Learning Research (09/2023)
5.2 Characteristics of CV Algorithms
The data splitting mechanism of CV procedure aflects the spatial dependence between training and vali- dation folds. Random partitioning methods, like standard k-fold cross-validate valsdavliittion (KFCV) and importance- weighted cross-validation (IWCV), preserve the spatial relationships in the training set, so if the training samples are spatially autocorrelated, then the training and validation folds will be spatially dependent. Blocking cross-validation (BI.CV) reduces spatial dependence across folds, though some dependency may remain among points around the splitting boundaries. Intuitively, with more blocks, more spatial depen- dence remains. Buffered cross-validation (BFCV) makes training and validation folds spatially independent by choosing a buffer least large the spatial autocorrelation, at the cost of losing some data points to the buffer.
CV algorithms also differ in their ability to deal with covariate shift. IWCV is the only one among them that explicitly addresses covariate shift. The spatial partitioning methods, like BLCV and BFCV, have no mochanism to rectify covariate shift but may actually introduce covariate shift when splitting the training set into folds. When spalially disparate folds imposed over spatially autocorrelated features, the folds may be sampling different parts of feature space. The extent to which this issue arises is function of the degree and pattern of spatial autocorrelation as well as tuning parameters of the approaches like the buffer or block size. In contrast, random partitioning like KFCV is not likely to introduce covariate shift across the folds because it samples uniformly at random across geographic space, which translates to sampling uniformly at random across feature space, even for autocorrelated features.
The approaches that the CV algorithms take to spatial autocorrelation and covariate shift are summarized in Tab. [ The last line of the table introduces our inspiration for new method, which we propose in the next section.
6 Proposed Method
Algorithm LOOIBCV
Input: training set {Ti}% = {xi Mi and density ratios {wi)
Parameters: buffer size r
Output: estimated error Err
for = 1 to do
2: Remove data points within the buffer area based on the longitude (long) and the latitude (lat) of T,
[Long long + r, lat - lat + r].
Fit model f on the remaining data T-i-r.
4: Calculate density ratio weighted loss on the validation fold T3 Err: wi L(y,, (xix):--T))).
5: end for
6: Return the estimated error: Err = 1 E. Erry
Inspired by scenarios in which training and test features are spatially independent and covariate shift is present (SI+ + CS), propose Importance Weighted Buffered Cross- Validation (IBCV; Alg Alg.[), which inherits the advantages of both BFCV and IWCV It isolates training and validation folds by adding buffer area between them and applies density ratio weighting to correct the bias from covariate shift. Modified from the IWCV estimator Sugiyama et al. 2007), the K-fold variant and the leave-one-out variant of IBCV are
((a) 1 E me
K
Eieteh Pio(Xi) C(yss (x);;TT Pee(Ni)
KKBBO
K
(m) FLOOBBV
1 Σn (x)( L(vi, Dr (x;;T
psr
where T-k-p and T-,-r are the training fold, n& is the size of the validation fold, i.e., the 4thh fold, and (x) is the density ratio weight. The buffer size r is decided by the range of spatial autocorrelation, e.g., the largest semivariogram range of all features. In geostatistics, the variogram/semivrriggaan range indicates
px-(xi)

---

Published in Transactions on Machine Learning Research (09/2023)
the threshold distance at which two observations are longer spatially dependent, i.c., they are spatially independent. With the same settings and assumptions of Thm. 1 we show that IBCV is asymptotically unbiased.
Proposition 2. IBCV is asymptotically unbiased: Ey [((m) [B = R(() when n - oc.
Proof. We demonstrate the claim for LOOIBCV; also valid for KIBCV with more bookkeeping for the folds. Step (1) below holds because T-,-, and Т, are independent. Step (2) holds because T-i-r and T3 are independent. We use shorthand L, for L(yi, De (x:;TT-->)) and Lj for (115)1)(55) (y::T---)).
Er LooIBovl
p(n)
E Er
(ee Ci] Per (25)
vi
[//
/
Pte
E Er
( Per (x1, yi) Ledxidyi
Per (xi)
1
E ET 1 //
Pio L
(2)
Σ ii ET_.... [C,]
Σ p(n-11nnr) ~3 R(())
71
Potential users IBCV should weigh a few caveats. First, IBCV may not be suitable for small datasets. Removing buffer points will have stronger effect for smaller datasets, potentially leading a pessimistic bias. Second, IBCV should be expected struggle with severe covariate shift just other importance- weighed methods do. Such methods perform poorly when the support of Px., has little overlap with the support of Px..
7 Experiments
In this section, we report on empirical evaluations of our proposed framework and on the performance of several CV estimators on simulated and real datasets. First, we used simulated datasets to estimate the bias of each CV algorithm in terms risk, following up on Thm. II This was possible because we could create multiple simulated landscapes under identical conditions. Second, we created example applications by subsetting from dataset of bird observations in Oregon, USA to create a diverse of geospatial scenarios. Finally, to further explore performance of the proposed IBCV algorithm in Scenario SI CS, we also applied the algorithms to datasets comprising bird observations in Alaska and California. housing pricos. On these real datasots, it was not possible to estimate risk directly, since we had only one landscape per analysis, so we measured test error as a proxy. To this end, the real datasets were each divided into a training set and a test set, each with known values for the features and response variables. The ground truth test error was the classification error rate or RMSE computed by fitting the model to the training set and applying
the test set. The estimated test errors were computed via each of the vross-validation procedures on the training set only and compared the ground truth test error. The subsections below describe each analysis, with more experimental details in the appendix. Code and data are available https://githnn http://cchh..aaoaaaaooaaaaaaaaaa
7.1
Simulation Experiments
7.1.1
Simulated Data
The simulation settings were inspired by the four main dataset scenarios discussod above (Tab. 1 with subscenarios for some cases and four autocorrelation ranges per subscenario. To construct synthetic land- scapes, simulated two independent, identical, spatial processes Zn and Z2 with underlying features (1)
10

---

Published in Transactions on Machine Learning Research (09/2023)
(a) r
(b)
(c) r 8
(d)
r
12
Figure Simulated spatial processes (such as those used for Z1 and Z2) with mcan=0, variance=1, and four variogram ranges v: (a) no, (b) mild, (c) moderate, and (d) strong spatial autocorrelation.
10
Feature 2
Feature 2
Feature
-3-
2.5
Feature
Feature 1
Feature
(a) Inside
(b) Partial
(c) Outside
Figure Joint probability density functions of the features in Scenarios SI + CSi (a) and SI + CSp (b). Cases like (c) are beyond the scope of this paper. These examples Blue and orange points denote training and test set feature values, respectively. The inner, middle, and outer contours illustrate the 68%, 95%, and 99.7% quantiles.
(i.e., no spatial autocorrelation)
and (2) 60 x 60 grid (Fig. D. We defined spatial random process by mean function m(.) and variogram function 7(h). Three parameters are used to describe variograms: the sill(s) is the total vari- ance; the nugget(n) is the small-scale variability of the data; the range(r) is the lag or distance where the variogram levels off to the sill. Pairs of points at least as far apart as the range are not considered spatially autocorrelated. In the simulation, we adopted the Ganssian variogram below:
7(h) (s (s- expl (cpp--13/72)) +n,
(5)
where h = |4 - 221 is the distance between any pair of two points 2,, 2,, 2y sampled from a contimuous two- dimensional spatial domain. Tab. specifics the parameters of the spatial processes for cach scenario. In Scenario SI CS, these variants differed based on whether the support of Px,, was contained within (Scenario SI + CSI, Fig, (h) partially overlapped with (Scenario SI + CSp, Fig. 2b) the support of Px... We consider extreme covariate shift, such that the test distribution is entirely outside the training distribution (Fig. 22), beyond the scope of this paper. Below, we treat SI + CSi as the main representative of the SI + CS scenario; we treat. ST + CSp and SI + RS as more seannniinn/ppecialzeeeee cases. The response variable y was generated as:
y xx) +x +1(1) 2)
11

---

Published in Transactions on Machine Learning Research (09/2023)
Table 3: Parameters of the spatial processes E1, Zg in each scenario: spatially dependent (SD), spatially independent (SI), spatially dependent with covariate shift (SD + CS), spatially independent with covariate shift (SI + CS: CSi inside, CSp- partially overlapped configurations), and spatially independent with range We set n = 0 for all scenarios. Note that indicates spatial autocorrelation the resolution of the 60 x 60 landscape is 1 grid cell.
shift (SI RS). The hyperparameters 8, r and n refer to Eqn. [
Intended Scenario
Training Set
Test Set
SD, SI, SD + CS
m(.) = 0, s = 1, r = 1, 4, 8, 12
m(.) = 0, 8 = 1, 1, 4, 8, 12
SI + CSi
m(-) = 0, s= 1, r = 1, 4,
8, 12
m(.) = 0,, 8 = 0.5, r 1, 4, 8, 12
SI + CSp
m(-) = 0, 8,8== 1, r= = 1, 4, 8, 12
m(.) =
1, 4, 8, 12
SI + RS
m(. = 0, s- 1, ---==0 12
m(.)) = =0,8= 10=11.-- 1- 1,
4,
60
60
60
Latiude
40
Latitude
40
Lastude
Laciture
40
20
20
20
40
60
20
40
60
20
40
60
100
120
140
160
Longitude
Longitude
Longitude
Longitude
(a) Scenario SD
(b) Scenario SD + CS
(c) Scenarios SI, SI CS and SI RS
Figure 3: Examples of spatial sampling in simulations from each scenario. Training points are blue; testing points are orange. Note that (c) shows points far apart in geographic space.
where (1) and 22) are the (potentially spatially autocorrelated) features, and c is an i.i.d. Gaussian error term. For every subscenario and range value combination, we repeated this procedure 100 times to allow us to estimate the risk.
Every dataset contained 1800 training and 500 testing samples, but the scenarios varied in the geographic sampling of these grid locations. In Scenario SD, training and test locations were uniformly sampled from the same spatial domain such that test points were interspersed among training points (Fig, Ba). In Scenario SD + CS, they were sampled non-uniformly (e.g., with spatial pattern) from the same spatial domain, so the average distance between training and test points was larger (Fig. 3bp. In Scenarios SI, SI + RS and SI + CS, training and test points were uniformly sampled from different spatial domains (Fig. 3c For the case of S1 + RS, the training landscapes had autocorrelation range 12, while the testing landscapes had different ranges. For SI + CS, the generating distributions changed from training to test landscapes, either concentrating the test features within smaller region of the support of the training landscapes (CSi) creating partial overlap between the densities (CSp). Consistent with our proposed framework, in Scenarios SD and SD+ + CS, the minimum distance between training points and test points is smaller than the variogram range (i.e., d < r in Fig. 3a and Fig. 3b), and in Scenarios SI and SI + CS, the opposite is true (i.c., 2 r in Fig. 3c The sampling methods and the underlying Gaussian processes are summarized Fig. and Tab.
Despite designing these simulations around each scenario in Tab. I some sampled datasets did not produce the expected result when testing for covariate shift (Tab. 4p. To understand why, the key thing to consider is the role of spatial autocorrelation and its potential to induce covariate shift in samples from these simulation processes. Even when Z, and Z2 are independent spatial processes, their underlying distributtioss present stronger spatial pattern as the range increases. With stronger spatial pattern, even a geographically uniform sample from the processes presents more skewed representation of the fcature space. In Scenario SI, where there are two separate realizations of these processes, the likelihood of observing covariate shift between them increases with the range of the underlying spatial autocorrelation. In Scenario SI + RS, range shift
12

---

Published in Transactions on Machine Learning Research (09/2023)
Table 4: Proportion of simulations that reject the Cramér test (a = 0.01). We summarize the major characteristics of each scenario based on the above statistical test results from two perspectives: spatially dependent (SD) or spatially independent (SI), and with or without covariate shift (CS). The omitted cases are those that do not fit the scenario definitions (i.c., SD does not occur with = 1 and RS does not occur
for r 12 since the training landscape also has r = 12.
Intended Scenario
SD
SI
SD
+ CS
SI + CSi
SI + CSp
SI + RS
F' = 1
0.00
0.00
1.00
1.00
0.17
1
0.01
0.00
0.86
1.00
1.00
0.03
r = 8
0.01
0.10
0.98
1.00
1.00
0.27
r = 12
0
0.33
1.00
1.00
1.00
exacerbates these differences. In Scemario SD + CS, when spatial sampling is conducted on a single landscape with no spatial autocorrelation (r = 1), almost all samples may be identically distributed. However, when spatial sampling is conducted on the distribution with an increasing spatial autocorrelation (r = 4, 8, 12), the likelihood of covariate shift increases, since similar feature values are more clustered and more often sampled together (i.e.; sampling Fig. 3b on a landscape like Fig. [d)) Minor discrepancies in Tab. 4 (i.e., rejecting 1% of simulations in some SD scenarios) are simply result of the random sampling and the statistical testing framework. In the results below, we excluded all simulations that were not sorted as intended by the Cramér test, to allow appropriate calculations of risk.
7.1.2 Simulation Experiment Design
Our primary goal was to measure the bias of each CV method in each scenario, to inform recommendations about which approach to recommend for different types of datasets. Toward this end, we defined the following model class for all simulation experiments:
g -= Bo + (1x(1) + 82x
We chose this linear model, which is missing the interaction term of the data-generating model, to mimic the usual case of model misspecification in real applications. We measured performance with root-mean-squaet squared error (RMSE). To evaluate the framework and CV methods, we estimated the risk (Eqn. 2 by averaging the test crrors of these models across the simulations in each subscenario. The bias of cach CV algorithm is the difference wetween the average CV estimate (over the simulations) and the est risk. The exception to this is Scenario SD+ + CS. Since the distributions of the training and test features change from one simulation to the next, we are not able to compute the expectations necessary to the target risk by averaging across simulations. In this case, we used the absolute bias as a proxy, defined as the average absolute difference between the CV estimates and the test error on a per-simulation basis: n E ICVest. - testerror,), where n is the number of simulations sorted into the scenario.
We also compared methods focusing on test error instead of risk, for two reasons. First, even though cross- validation actually estimates risk, practitioners are often interested in test error, given that it conditions on the training set at hand. Second, in the real data analyses below, we can calculate test error but not the true empirical risk, since again we have only a single dataset rather than the simulated replicates. Looking at both in simulated studies may provide insight into the differences between these viewpoints.
Some of the CV methods had hyperparameters to set. We set k = 9 folds for all cross-validation methods. The standard KFCV and importance- weightted IWCV split the training set randomly (Fig. [p). Blocking BLCV partitioned the training samples into blocks and then put these blocks into 9 folds (Fig. Ap). We set the block size of BLCV = 2, 8, 12 grids respectively when r 1,4, 8, 12, such that the block size would mimic the spatial aulocorrelation range. This is consistent with the motivation of BLCV, which is to break up spatial autocorrelation. Note that BLCV is sensitive to the block size, and different settings of this hyperparamcter produce different results (Fig. in Appendix). In particular, spatial autocorrelation, the average BLCV estimate of test error tends to increase with increasing block size, inter-block covariate shift increases (Fig. in Appendix). It is not trivial in general to set block
the presence of
13

---

Published in Transactions on Machine Learning Research (09/2023)
60
60
60
48
40
40
Latitude
Latitude
36
Latitude
24
20
20
12
0-
20 Longitude
0
12
36 Longitude
24
48
60
40
Longitude
(a) KFCV and IWCV
(b) BLCV (bs=12)
(c) BFCV and IBCV
Figure 4: Visualization of how different CV methods split 1800 training data points a 60 x 60 landscape into 9 folds.
(a) Each color represents fold. (b) An example of BLCV with block size of 12 grid cells, each assigned to one of the 9 folds. We used block sizes of 2, 4, 8, and 12. (c) BFCV and IBCV were based on a grid with 20x20 cells. This example shows fold 5 as the validation fold, where training samples in its surrounding buffer region (buffer size 12 here) have been removed.
size optimally (Valavi et al. 2018). Buffered BFCV and importance weighted buffered IBCV use a similar grid to BLCV, but with block size of 20, which produces 9 blocks on the 60 60 landscape, each of which becomes its own fold (Fig. 4 We removed all training points within buffer of each validation fold in that the
turn. The size of the buffer was equal to the spatial autocorrelation range of the simulation minimum distance between the training and validation folds was equal or greater than the range, and thus the validation fold could be considered spatially independent. We chose a block size of 20 instead of varying it in order to reduce the amount of data falling into the buffer; with multiple validation folds, the overall bufler area would increase, resulting in fewer training data points. Finally, IWCV and IBCV required density ratio estimates. We applied the Relative unconstrained Least-Squares Importance Fitting (RuLSTFF) method to estimate density ratios (Yaaada et al., 2011), with = 0 and 50 kernels.
7.1.3 Simulation Results
Here, we present results for the four main scenarios: SD, SI, SD+CS, and S1+CSi (denoted SI+CS below). We consider SS+RS and SI+CSp to be auxiliary cases of CS with SI+RS being somewhat uncommon and CI-CSp being arbitrarily difficult (as we move from SI+CSi through SI+CSp to SI+CSo, the covariate shift becomes increasingly unmanageable). Measurements of bias in the risk estimates of each CV method across the simulated datasets for each scenario are presented graphically in Fig. and numerically in Tab. We examine the results for each scenario below.
In Scenario SD, where Xtr and Xtc are spatially dependent and distributed identically, KFCV is least biased (Tab. 5) Its internal random partitioning mechanism maintains the same spatial relationship between training and validation folds as between training and test sets. The estimates for all methods show increased variation as the range of spatial autocorrelation increases (Fig. 5p). The bias incurred by IWCV reflects estimation error for the density ratios; if the density ratios are set (as correct for the SD populations), KFCV and IWCV are equivalent. Similarly, IBCV would be equivalent to BFCV with density ratios of 1, so the discrepancies reflect the density ratio estimation process. The other three spatial CV variants alleviate the spatial dependence between points across fold boundaries; in this case, this is unhelpful because while the folds are decorrelated with each other, the training and test sets are correlated. These methods may also introduce some covariate shift among the folds. In both ways, they overestimate errors.
In Scenario SI, where X,r and X,e are spatially independent and distributed identically, BFCV shows promise (Tab. Fig. 5p). Some spatial dependence remains along the boundaries of folds, which may cause under- estimation of error, but spatial separation of the folds can also introduce covariate shift, which may cause
14

---

Published in Transactions on Machine Learning Research (09/2023)
1.5
1.5
1.0
1.0
0.5
0.5
Bias
Bias
0.0
0.0
-0.5
-0.5
-1.0
-1.0
-1.5
-1.5
8
12
4
12
Variogram range
Variogram range
(a) Scenario SD
(b) Scenario SI
1.5
1.5
1.0
1.0
Absolute Bias
0.5
0.5
Bias
0.0
0.0
-0.5
-0.5
-1.0
-1.0
-1.5
-1.5
4
8
12
4
8
12
Variogram range
Variogram range
(c) Scenario SD CS
(d) Scenario SI + CS
CV algorithm:
KFCV
BLCV
BFCV
IWCV
IBCV
Figure Biases of CV estimates scenarios with various characteristies: spatially dependent. (SD), spa- tially independent (SI), spatially dependent with covariate shift (SD + CS), and spatially independent with covariate shift (SI CS). Circles inside the boxes display the mean values of biases. The black dash lines illustrate bias. (a), (b) and (d) show bias of the average CV estimate to the risk. Since the feature distributions change across simulations in Scenario SD + CS, (c) plots the absolute bias, i.e., the absolute value of CV estimate minus test error in each simulation.
15

---

Published in Transactions on Machine Learning Research (09/2023)
Table 5: Risk perspective: Bias of the average CV estimate (9 folds) with respect to risk, across varying degrees of spatial autocorrelation (i.e,, training set range 1, 4, 8, 12). Bias is calculated the mean CV estimate minus the risk except for Scenario SD + CS, where bias is calculated as the average of absolute differences between CV estimates and test errors. The least biases in each row of the same scenario block are bold.
KFCV BLCV
BFCV
IWCV
IBCV
KFCV
BLCV
BFCV
IWCV
IBCV
Scenario SD
Scenario SI
0.0022
0.0023
0.0024
-0.0445
-0.0443
0.0016
0.0017
0.0018
-0.0435
-0.0433
0.0103
0.0174
0.0280
-0.0517
-0.0354
0.0023
0.0094
0.0200
-0.0611
-0.0449
-0.0024
0.0030
0.0691
-0.0605
0.0062
-0.0299
0.0046
0.0388
-0.1106
-0.0482
-0.0010
0.0778
0.1602
-0.0550
0.0956
-0.0438
0.0347
0.1123
-0.1778
-0.0477
Scenario SD +
Scenario S1 + CS
0.4977
0.4978
0.4979
-0.0052
-0.0052
0.1043
0.1068
0.0973
0.1651
0.1453
0.4883
0.4954
0.5060
-0.0018
0.0059
0.1938
0.1861
0.1853
0.2524
0.1872
0.4603
0.4955
0.5320
-0.0652
-0.0347
12
0.2472
0.2492
0,2583
0.2472
0.1870
0.4429
0.5217
0.6041
-0.0646
0.0070
Table 6: Test error perspective: Proportion of sorted simulations in which each CV algorithu is the closest estimate to the test error, for varying degrees of spatial autocorrelation (i.e., when the range of training set 4, 8, 12). The highest numbers in each row of the same scenario block are in bold.
KFCV BLCV
BFCV
IWCV
IBCV
KFCV
BLCV
BFCV
IWCV
IBCV
Scenario SD
Scenario SI
0.23
0.20
0.25
0.11
0.21
0.25
0.24
0.20
0.16
0.15
0.20
0.07
0.35
0.22
0.15
0.17
0.11
0.36
0.25
0.11
8
0.23
0.21
0.15
0.26
0.14
0.12
0.10
0.42
0.23
0.12
12
0.37
0.17
0.03
0.30
0.13
0.19
0.10
0.39
0.28
0.03
Scenario SD
CS
Scenario SI + CS
0.00
0.00
0.00
0.57
0.43
0.16
0.09
0.45
0.08
0.21
0.00
0.00
0.00
0.55
0.45
0.17
0.09
0.31
0.07
0.36
0.00
0.00
0.00
0.19
0.81
0.10
0.15
0.25
0.15
0.35
0.01
0.00
0.00
0.32
0.67
16

---

Published in Transactions on Machine Learning Research (09/2023)
overestimation of error, especially when the range is large. Setting the block size based the spatial autocorrelation range appears to be reasonable heuristic in these results, balancing these two potential influences. The upward bias of BFCV increases with range, more substantially than BLCV; this may reflect more covariate shift incurred by BFCV than BLCV and/or the data loss incurred by buffering. Again, the differences belween KFCV and IWCV and between BFCV and IBCV are attributable to the density ratio estimates.
In Scenario SD + CS, where Xxr and Xte are spatially dependent and distributed differently, we used the average of absolute bias to the test error (|Rev Errrl) across simulations compare the CV algorithms (shown by the Y-axis of Fig. E()) The results omit the case where r 1 (no covariate shift). BFCV performed best with mild and moderate spatial autocorrelation, and IBCV performs best with strong spatial autocorrelation, but the results are similar across all methods (Tab. Fig. 5P)
For Scenario SI + CS, where Xrr and Xte are spatially independent and distributed differently, the CV algorithms with density ratio weightings have clear advantage (Tab. 5 Fig. F(). As may be expected, the improvements of IBCV over IWCV increase with variogram range; ie., spatial autocorrelation increases, the importance of the buffering component of the CV method increases best on the auxiliary SI+CS cases in the appendix as well, with the caveat that the overall quality of the estimates from all methods decreases with increasing amounts of covariate shift..
well. IBCV usually performs
Comparing CV methods in terms of test error instead of risk yields slightly different results. Tab. G lists the proportion of simulations for which each CV method's estimate was closest to the test error for that simulation (i.e., how often each CV method was best KFCV was least biased regardless of range, the best method from test error perspective varied, with BFCV and IBCV sometimes performing better on a larger proportion of simulations. Results from Scenario SI show similar trend in that KFCV is best for low r while spatial CV approaches perform better as range increases, but the test error analysis favors BFCV while the risk analysis favors BLCV. Results for the scenarios with covariate shift showed more similarity between the two perspectives, in that importance weighting and buflering methods provided best performance.
per-simulation basis)- For Scenario SD, while
7.2 Example Application
In this section, we provide a real-world example of how different CV methods may affect geospatial analyses. We divided a dataset on bird species distributions in Oregon for training and testing in four distinct ways to set up the different scenarios described above. We were interested to see whether the CV methods that best estimated test performance matched the results from the simulation study above. However, note that in real applications like this, we have only one dataset (in contrast to 100 simulated datasets), we only compare estimated and true test errors (Eqn. 1 rather than risks (Eqn. 2 This is more similar to the analysis presented in Tab. than Tab. 5
7.2.1 Oregon 2020 Data
The Oregon 2020 project collected bird species observations by both citizen scientists (a.k.a. community scientists) and professional ornithologists across the state of Oregon by 202 (Robinson et al.] 2020b. The project was designed to rigorously sample both the geographic space and habitat variety of Oregon. We focused on two common species in the Oregon 2020 dataset: the Hermit Warbler (HEWA) and the Western Tanager (WETA). Our analysis task was to build species distribution models (SDMs), which learn and prodict whether bird species was obscrved or not across the region of analysis, basod on environmental features (Elith & Leathwick 2009). SDMs are valuable tools for scientists and natural resource managers (Phillips et al.) 2004). At each bird observation location, assembled environmental features to characterize the surrounding habitat. The four features were vegetation indexes computed from remote sensing data using the Tasscled Cap transformation of LandSat imagery (Crist & Ciconel 1984): the mean values of the Tasscled Cap angle (TCA), brightness (TCB), greenness (TCG), and wetness (TCW) in a 600-m radius area around
Project website: https://oregon2222..omm
17

---

Published in Transactions on Machine Learning Research (09/2023)
46
46
Lattude 88
45
Latitude
45
Latttude
45
Lattude
45
44
44
44
44
43
43
42
42
42
42
-125.0
-122.5 -12.00 -17.55
-125.0
-12..5 1122.0-12.0-----10 -117.5
-125.0
-122.5
-120.0
117.55
-125.0 -122.5
120.00 -117.5
Longitude
Longitude
Longitude
Longitude
(a) Scenario SD
(b) Scenario Si
(c) Scenario SD CS
(d) Scenario SI+
Figure 6: Training and test sets for the four scenarios in the HEWA1800 dataset. Training points are blue and testing points are orange. Subplots (a), (c) and (d) share the same training set.
Table 7: Feature statistics for the Oregon 2020 datasets: Range lists the maximal semivariogram range of training features in degrees. Cramér lists p-value of the Cramér two-sample test whose mull hypothesis is that training and test samples are from the same distribution.
Dataset
Range
Cramér
Dataset
Range
Cramér
HEWA 1000 SD
0.28
0.7123
HEWA 1000 SI
0.33
0.0430
HEWA1 180 SD
0.30
0.6184
HEWA 1800 ST
0.32
0.1389
WETA1800 SD
0.27
0.0999
WETA 1800 SI
0.28
0.0649
HEWA 1000 SD+CCS
0.28
0.0000
HEWATONM ST+CS
0.28
0.0000
HEWA 1800 SD+CS
0.30
0.0000
HEWA1800 SI+CS
0.30
0.0000
WETA1800 SD+CS
0.27
0.0000
WETA1100 SI+CS
0.27
0.0000
each survey site. Previous work found these four variables to be predictive of bird species in Oregon (Hopkins lt al. 2022).
We created training and testing datasets for three combinations of species and dataset size, each with four geographic layouts to set up the four scenarios of Tab. [ We assembled datasets with either 1000 or 1800 training observations and in each case tested on 500 held-out observations; datasets are named by their species abbreviation and training sample size (i.e., HEWA1000, HEWA1100 and WETA1100). In all cases, we randomly removed samples from the majority class to produce balanced binary classification tasks, in order to control for any effects of imbalanced data on the classifiers building the SDMs. In Scenarios SD and SD CS, the training and test samples were spatially proximal (Fig. a and c). In Scenarios S1 and SI + CS, the minimum distance between training and test points is larger than the maximum semivariogram range among the four features (Fig, 6 and d). The p-values of the Cramér test in Scenarios SD and SI for all three datasets are greater than = 0.01, so we did not reject the mull hypothesis; the training features were distributed the same as the test features. The p-values are all zeros in Scenarios SD + CS and SI + CS, so we rejected the null hypothesis; covariate shift was present. In summary, each dataset meets the criteria for the intended evaluation scenario.
7.2.2 Oregon 2020 Experimental Design
As in the simulation experiments, we sought to determine the best CV approach for each SDM. Whereas in simulation we fit a single regression, we explored five classification models for cach SDM: Ridge classifier (Ridge), Linear SVM (LSVM), K-Nearest Neighbors (KNN), Random Forest (RF), and Naive Bayes (NB), and we compared their test errors with the CV error estimates. We used the default hyperparameters from the scikit-learn Python package (Pedregosa et al., 2011) for all models. While peak-to-peak comparison with tuned hyperparameters is critical for comparing state-of-the-art models, an unbiased CV procedure should be able to estimate model performance correctly whether it is high- or low-performing. We evaluated models with classification error rate: the proportion of samples misclassified.
An example showing the splits produced by the CV methods for the HEWA1800 dataset given in Fig. 7 The CV hyperparameters for HFWA1800 were set follows; the other analyses followed an analogous
18

---

Published in Transactions on Machine Learning Research (09/2023)
46.0
46.0
46.0
45.5
45.5
45.5
Latiude
45.0
Lattude
45.0
Latiude
45.0
44.5
44.5
44.5
44.0
44.0
44.0
43.5
43.5
43.5
-125
-124
-123
-122
-125
-124
-12
-122
-125
1124
-123
-122
Longitude
Longitude
Longitude
(a) KFCV, IWCV
(b) BLCV
(c) BFCV, IBCV
Figure Splitting strategies of five CV algorithms HEWA1800, Scenario SD dataset when block size 0.30 degrees: eight training folds (green points) and validation fold (red points) in each subplot. Compared with (b), training samples in buffer regions have been removed in (c).
procedure. To calculate ranges of the features of training sets, we fitted Matérn variogram functions with the lag class estimated by Scott's rule (Mälicke 2021). The maximal ranges among the four features are 0.30, 0.32, 0.30 and 0.30 degrees in Scenario SD, ST, SD+CS, ST+CS, respectively (Tab. 7 We set buffer size equal to the maximum range over all features to make training and validation folds spatially independent. We fine-Luned the block size among values equal to the range, double the range, and triple the range, specifically [0.30, 0.60, 0.90] for SD, SD+CS, and SI-CS. For SI, tripling the range produces fewer than 9 folds, so we only selected among [0.32, 0.64]. The Relative unconstrained Least-Squarcs Importance Fitting (RuLSIF) method was applied to estimate density ratios Yamada et al. 2011 again with 0 and the number of kernels = 50.
7.2.3 Oregon 2020 Analysis Results
This analysis highlights the importance of the hyperparameters of the CV approaches. In Scenario SD, the test errors generated by the CV methods usually overestimate the true test error (21/25 cases in Tab, g and the best method varies across classifiers. In this case, setting the hyperparameters based on the range of the features works well; selecting the best hyperparameters post hoc matches this heuristic. In Sconario SI, there is variation the best. CV approach again, and CV approaches underestimate the true error for all classifiers. Here, setting the hyperparameters for the spatially explicit methods optimally would change the results for KNN and RF, such that BLCV would be the best choice for both. In Scenario SD CS, all methods underestimate test error for all classifiers when setting the hyperparameters based on feature ranges, and KFCV performs best for four out of five classifiers. With optimal hyperparameter selection, spatially explicit approaches would be better in three of these cases. In Scenario SI + CS, IWCV is the best CV approach for four out five classifiers, and IBCV for the fifth. With optimal tuning, three of these would shift from IWCV to IBCV. Without hyperparameter tuning, most test errors are underestimated (24/25 cases); with tuning, this is still often the case (11/15 cases).
The results for the other Oregon 2020 datasets (HEWA1000, WETA1800) had similar variability to HEWA1800 (Appx. Tab. 13 Appx. Tab. 14). In all scenarios, test errors for WETA1800 were roughly double those of the HEWA analysis; in these regions, HEWA is more of a habitat specialist (clcarly prefer- ring some habitats over others) than WETA (generally widespread), so the environmental features may be less informative overall for the WETA distribution. In Scenario SD, analysis of the other datasets provided less support for spatially explicit CV approaches, and for WETA1800, setting hyperparameters a priori versus post hoc showed some differences. For Scenario SI, the analyses were more consistent in supporting spatially explicit CV approaches; BLCCV and BFCV were usually provided the strongest. estimates. Again
19

---

Published in Transactions on Machine Learning Research (09/202))
Table 8: HEWA1800: model classification error rates (targets) and 9-fold CV estimates thereof (best estimates in each column in bold). BLCV-best, BFCV-best and IBCV-best estimates are selected from the best ones in the Appx. Tab. 15 for peak-to-peak comparison. BLCV-range, BFCV-range, and IBCV- range set the tuning parameters priori based on the range of spatial autocorrelation in the features. A dash means that setting the hyperparameters based on the range gives the best value (i.e., the methods are equivalent).
Model
Test error (target)
KFCV
IWCV
BLCV
BFCV -range
IBCV
BLCV -best
BFCV
IBCV
-range
-range
-best
-best
SD
Ridge 0.1640
0.1844
0.1849
0.1786
0.1780
0.1781
LSVM
0.1640
0.1850
0.1855
0.1752
0.1842
0.1846
KNN
0.2020
0.2067
0.2070
0.2084
0.2035
0.2031
RF
0.2040
0.1972
0.1971
0.2057
0.1945
0.1938
NB
0.1720
0.1878
0.1884
0.1810
0.1757
0.1757
SI
Ridge
0.2000
0.1800
0.1797
0.1711
0.1842
0.1839
LSVM
0.2040
0.1828
0.1825
0.1805
0.1989
0.1986
KNN
0,2260
0.1861
0.1858
0.1753
0.1697
0.1695
0.2117
0.2028
0.2025
RF 0.2580
0.1894
0.1892
0.1700
0.1950
0.1947
0.2215
0.2114
0.2111
NB
0.2140
0.1883
0.1881
0.1892
0.1853
0.1850
SD
CS
Ridge
0.2140
0.1844
0.1842
0.1786
0.1780
0.1778
0.2083
0.2120
0.2117
LSVM
0.2060
0.1850
0.1848
0.1752
0.1842
0.1840
0.2092
0.2081
0.2078
KNN
0.2120
0.2067
0.2064
0.2084
0.2035
0.2032
RF
0.2000
0.1972
0.1970
0.2057
0.1945
0.1942
NB
0.2160
0.1878
0.1875
0.1810
0.1757
0.1755
0.2131
0.2114
0.2111
SI + CS
Ridge
0.2400
0.1844
0.2196
0.1786
0.1780
0.2085
0.2083
0.2120
0.2397
LSVM
0.2460
0.1850
0.2207
0.1752
0.1842
0.2167
0.2111
0.2148
0.2421
KNN
0.2420
0.2067
0.2477
0.2084
0.2035
0.2370
0.2439
0.2458
RF
0.2380
0.1972
0.2341
0.2057
0.1945
0.2239
0.2438
0.2446
NB
0.2660
0.1878
0.2248
0.1810
0.1757
0.2059
0.2131
0.2114
0.2403
20

---

Published in Transactions on Machine Learning Research (09/2023)
62
42
61
60
Latitude
Latitude
Latitude
38
59
36
58
57
34
56
158
156
154
152
124
122
120
118
116
114
124
122
120
118
116
114
Longitude
Longitude
Longitude
(a) Alaska's Birds
(b) Housing by Latitude
(c) Housing by Bay
Figure 8: Alaskan National Parks and California Housing dataset spatial layouts. Training points are blue and testing points are orange.
though, the default hyperparameter settings were not the best. The Scenario SD + CS analyses were variable again, with each CV method being the best estimator in at least one case and hyperparameter settings being influential. Finally, in Sccnario SI + CS, there was strong support for importance weighting: 10/10 cases were best estimated IWCV or IBCV. With optimal hyperparameters, IBCV would be best in all cases.
7.3 Further Analysis of Scenario SI + CS
We analyzed two more datasets to further explore the performance of IBCV in particular. We performed classification for another species distribution modeling problem using bird data from Alaska, and we explored a regression problem on California housing prices, with two different spatial divisions of training and test points. Both datasets, described below, fell into the SI + CS scenario.
7.3.1 Alaskan National Park Data
This dataset records observations of 15 bird species across three national parks in Alaska, USA Amundson et al.] 2018). These surveys were conducted during May-June, 2004-2006 in Katmai National Park and Preserve (KATM) and Lake Clark National Park and Preserve (LACL), and during May-June, 2008 in Aniakchak National Monument and Preserve (ANIA). Wo set up balanced binary classification task for the most prevalent bird species, Golden-crowned Sparrow (GCSP). The features were the proportions of eight land cover classes: water, wetland (wetland vegetation), shrub (shrub vegetation). dshrub (dwarf shrub or herbaceous vegetation), dec (deciduous forest), mixed (mixed deciduous and evergreen forest), spruce (evergreen forest), baresnow (bareground or perennial ice and snow) within 50-m 150-m radius circular area around each survey site. The training set consisted of 710 data points (355 non-detections and 355 detections) collected from KATM and LACL in 2004-2006. The test set contained 134 data points (67 non- detections and 67 detections) collected from ANIA in 2008. As shown in Fig (c), the nearest distance of training and test samples is 2.01 degrees, greater than the maximal semivariogram range 0.82 degree. The p-valuc of the Cramér test
7.3.2 California Housing Data
This dataser Paace & Barry 1997 is derived from the 1990 census. The most common task with this dataset is to predict the median house value. We included six numeric features in our models: housing
2Lin the dataset:
Sink the dataset: http:///www.uaagg...ooooaaaaaaaa
21

---

Published in Transactions on Machine Learning Research (09/2023)
Table 9: S1 + CS real datasets: test error (targets) and 10-fold CV estimates thereof (best estimates each column bold). BLCV-best, BFCV-best and IBCV-best estimates are selected from the best ones in the Appx. Tab. 15 for a peak-to-peak comparison. BLCV-range, BFCV-range, and IBCV-range set the tuning parameters priori based on the range of spatial autocorrelation in the features. A dash line means setting the tuning parameters based on the range gives the best value (i.e., the methods are equivalent).
Model
Test error (target)
KFCV
IWCV
BLCV
BFCV
IBCV
BLCV -best
BFCV
IBCV
-range
-range
-range
-best
-best
Alaskan
birds
Ridge
0. 1.1866
0.3225
0.2526
0.2777
0.2881
0.2374
0.2696
0.2555
0.2345
LSVM
0.1866
0.3211
0.2488
0.2723
0.2791
0.2244
0.2714
KNN
0.3284
0.3113
0.2780
0.3030
0.3251
0.2542
0.3279
0.3254
RF
0.3657
0.3169
0.2911
0.3176
0.3343
0.2998
0.3460
0.3554
0.3483
NB
0.4030
0.3521
0.2938
0.3230
0.3597
0.2748
0.3522
0.3663
0.3667
Housing:
by latitude
Linear
0.5134
0.7102
0.4298
0.7024
0.7062
0.4286
0.5269
KRR
0.5931
0.7226
0.4350
0.7135
0.7184
0.4332
0.5594
SVR
0.3470
0.6496
0.4162
0.6430
0.6564
0.4178
0.4650
KNN
0.5246
0.7248
0.4483
0.7132
0.7389
0.4558
0.7370
0.5358
RF
0.4767
0.6825
0.4274
0.6728
0.7008
0.4318
0.4877
Housing:
by bay
Lincar
0. 1
0.6949
0.8121
0.6900
0.7199
0.8518
0.7170
0.7330
0.7793
KRR
7813
0.7063
0.8349
0.7006
0.7309
0.8722
0.7316
0.7394
0.7988
SVR
0. 0.6994
0.6295
0.7310
0.6262
0.6592
0.7744
0.6407
0.6687
0.6838
KNN
7830
0.7037
0.8232
0.6960
0.7430
0.8821
0.7143
0.7530
0.7859
RF
0.7176
0.6609
0.7679
0.6504
0.6966
0.8191
0.6702
0.7050
0.7336
median age, total rooms, total bedrooms, population, households, median income. We split the data into training and test sets in two ways: (a) by Latitude, where the test set includes 551 points above 39.5 N and the training sot includes 10968 points below 39 N; and (b) by Bay, whore the test set consists of 981 points within San Francisco Bay Area and the training set consists of 9775 points whose distance at least 0.3 degree from any test sample (Fig, S For cases (a) and (b), the nearest distances between training and test samples are 0.5 and 0.3 degrees, respectively, which are both greater than their maximal semivariogram
ranges of 0.23 and 0.25 degrees, respectively. The p-values of the Cramné tests are zeros for both cases.
7.3.3 SI + CS Experiments
For each analysis, we fit five models: Ridge, LSVM, KNN, RF and NB for the birds, and Linear Regression (Linear), Kernel Ridge Regression (KRR), Support Vector Regression (SVR), K-nearest Neighbors (KNN), and Random Forest (RF) for housing prices. We measured classification error rate and RMSE for the birds and prices, respectively. We split the training data with each of the CV methods for each case and measured performance on the held out test sets to serve as the target error for the CV estimators.
The results provide some support for the morits IBCV, but the importance of hyperparameter selection is evident. Overall, IBCV with the rasge-based hyperparameter heuristic was best. in 4/15 cases, but if optimal hyperparameters could be selected, it would be best in 10/15 cases (Tab. gh. When IBCV is not closest the target test error, there is still support for both the importance weighting and spatial buffering components of the approach; the other competitive methods are IWCV and BFCV.
8 Discussion
Overall, our investigation into cross-validation strategies for geospatial problems suggests that care war- ranted when designing evaluation schemes. Across the simulated and empirical data analyses, there is
22

---

Published in Transactions on Machine Learning Research (09/2023)
substantial variation in the results from different CV approaches. For example, the random forest trained on housing data in southern California and tested in northern California was evaluated as 43% above the true test error if using standard KFCV and 9% below the true test error with the proposed IBCV. Inaccurate assessments of model quality may cause problems in a variety of applications, especially when policies are based upon model predictions. These issues propagate into model selection as well; in some cases, the five models being compared had completely different rankings under different CV approaches (e.g., the Alaskan birds dataset under KFCV vs. IBCV).
While important, the task of comparing cross-validation strategies is challenging, particularly due to the distinction between risk and test error. Theoretical analyses of CV algorithms prove that they can produce asymplotically unbiased estimates of risk. Simulation experiments can bear this out, since we can sample landscapes from identical distributions and compute true and estimated risk from these replicated analyses. However, in real data experiments, we can only compute one true test error and one estimate from each CV algorithm. Our simulation results show that the CV estimator that is least biased overall not always a clear winner in terms of test error (Tab. 5 vs. Tab. G). This is a fundamental challenge, especially when a modeler has single training set to analyze and seeks performance estimates that condition on that particular dataset (i.c., test error).
In this work, we explored a discrete framework for sorting analyses into four scenarios. We characterized spatial autocorrelation and covariate shift as present versus absent, by thresholding continuous measures of these attributes. Naturally, this approach has some limitations, and it may be the case that some of the variation in our results comes from varying degrees of spatial autocorrelation and covariate shift among cases within each of the four discrete scenarios. For spatial autocorrelation, we applied strict distance check to the closest pair of points spanning the training and test sets, but in future work, we will explore more holistic characterization. For covariate shift, we thresholded the p-value from the Cramér test. In our experience, the Cramér test was very sensitive; in the simulation experiments, many p-values were zero, a few greater than 0.1, and only a small portion in between, so we expect that our results would not change substantially across values of a. In practice, for intermediate p-values (e.g., between 0.01 and 0.1), a practitioner might consider validation methods that do and do not address covariate shift. In future work, we will consider more flexible uses of the Cramér test as well as alternative metrics for characterizing covariate shift.
Our proposed algorithm combining the spatial buffering and importance weighting strategies, IBCV, shows some promise for challenging geospatial scenarios. In the simulation experiments where covariate shift was present and the test set was spatially distinct from the training set (SI + CS), IBCV was the best choice of CV method most of the time, both in terms of the bias of the risk estimates and the proportion of cases in which the test, error was closest to the truth. The results were more variable on the empirical datasets. IBCV was sometimes the best method with the default hyperparameter settings, but it would have been a more clear favorite if the hyperparameters were set optimally. In 17/20 cases of empirical SI + CS datascts presented above (and 6/10 in the appendix), IBCV yielded better estimates when using hyperparameters that differed from the default settings based on the range of the spatial autocorrelation.
The challenge of hyperparameter tuning was not unique to IBCV. The existing spatial CV methods (BLCV and BFCV) were also often improved by setting hyperparameters to values other than the defaults. Since the optimal hyperparameter results in the experiments above were basod on post hoc inspection, it is most appropriate to compare methods based on the default hyperparameler settings. The improvements made under the best hyperparameter settings indicate how much potential method may have, given a robust method for hyperparameter tuning. Currently though, setting these parameters optimally is non-trivial. Nested cross-validation procedures may be an option, but it is undesirable to introduce a hyperparameter search for cross-validation methods that are typically wrapped around other modeling algorithms, which may themselves require hyperparameter tuning. To realize the potential of spatial cross- cross-validation methods, continued research is needed on methods to setting block and buffer sizes appropriately.
Overall, we can make a few recommendations to practitioners based this study.
Scenario SD. When the training and intended testing datasets are interspersed spatially, random partitioning of the training set (KFCV) is appropriate; there is no covariate shift if both sets are evenly distributed across the landscape, and whatever spatial dependence structure exists between
23

---

Published in Transactions on Machine Learning Research (09/2023)
training and testing data is replicated between CV folds by random partitioning. If modeler chooses spatial partitioning instead, it could actually introduce pessimistic bias via covariate shift. especially for long-range autocorrelation.
Scenario SI. When training and test data are geographically separated, spatial partitioning methods like blocking (BLCV) and buffering (BFCV) are good choices, though as discussed above, perfor- mance sensitive to the size of the blocks or buffer.
Scenario SD + CS. The best CV approach was perhaps most variable for this scenario. Importance weighting (IWCV) can be valuable when covariate shift is at play, but spatial partitioning methods sometimes performed better than IWCV. This variation may be driven in part by the relative strengths of the two aspects of the problem (i.e., autocorrelation and covariate shift).
Scenario SI + CS. Our proposed approach of importance-weightted buffered valss-validation (IBCV) can address both of autocorrelation and covariate shift simultaneously. Again, performance was sensitive to the buffer size, and further work on setting this hyperparameter appropriately will help realize the potential of this approach.
Recent papers have come to mixed conclusions and recommendations about the value of spatial CV strategies, and our analysis yields points of both agreement and disagreement with the ongoing discussion. For example, Ploton et al. (2020) argue, with forest mapping example, that spatial CV should be the norm for ccological and biological studies of spatially autocorrelated processes. Our results do provide evidence for spatial CV strategies in some contexts, though we find that there are some spatial settings (e.g., Scenario SD) where standard k-fold cross-validation is appropriate (#1 above). Wadoux et al. (2021) disagree with Ploton et al. and find instead that spatial CV can be substantially pessimistic in its estimates. Our study finds evidence to support this perspective as well, particular for Scenario SD. Our agreements with these two disagreeing studies simply highlights one of our main messages: the appropriate evaluation strategy for given analysis is context dependent, particularly in relation to the predictive goals of the study, Another recent paper by Hoffimann et al. (2021) was more similar to ours; those authors also found that spatial blocking CV can produce biased estimates but that standard k-fold CV was also not adequate in the presence of covariate shift. A key difference ween their work and ours is that their analysis was limited to analysis of test error, whereas our simulation study looked at bias in both risk estimate and test error estimates. Our multi-faceted analysis provides additional perspective, since test error is more commonly assessed in practice, but theoretical guarantees about CV pertain risk estimates.
We expect this area of research to be of ongoing importance for applications like species distribution modeling. The ultimate goal of many SDMs is to make useful predictions about species adaptations changing climatic conditions, including novel 'no-analog" climate scenarios, which clearly present covariate shift. Such models sometimes predict that species will move into new spatial areas that they do not currently occupy but that are accessilee/ddjacent to their current distributions (e.g., shifts pole-ward and/or up in elevation). Alternately, natural resource managers may consider translocation of species to new spatial areas that are inaccessible to the species by its own dispersal means, as a climate mitigation strategy for threatened taxa. Beyond SDMs, remotely sensed data, are enabling the application of machine learning tools to variety of global monitoring applications (Lacoste et al. 2021), and as these methodologies grow, evaluation methods that are appropriate for spatial prediction tasks are ever more critical. Ongoing work in this field may serve biodiversity conservation and sustainability goals by providing more accurate model assessments for these challenging scenarios.
9 Conclusions
This paper has presenled an extensive investigation, on both simulated and empirical datasets, into cross- validation algorithms for geospatial applications. In these applications, the presence of spatial autocorre- lation and/or covariate shift. make it challenging to obtain unbiased estimates of model's generalization performance. Our theoretical results confirm that in the presence of spatial autocorrelation, the dependence structure between the training and intended test data should match the dependence structure of the CV
24

---

Published in Transactions on Machine Learning Research (09/2023)
folds into which the training data are divided. That is, rather intuitively, the CV folds should be designed to replicate the spatial structure of the intended application. For problems that include both spatial au- tocorrelation and covariate shift, we developed importance-weighted buffered cross validation which both removes data points within a buffer around test folds to break up the spatial dependence and uses density ratios to address covariale shift. IBCV shows promise for challenging spatial problems, though like other spatial CV approaches (based on blocking and buffering), the resulting performance estimates are sensitive to hyperparameters that are not always easy to set properly.
We see several directions for future work. While IBCV shows promise, further developments to aid hyperpa- rameter selection, and/or improvements to the method itself, may be fruitful. The framework developed in Tab. | served this investigation adequately, but to make more specific recommendations to practitioners, a more nuanced tool for characterizing spatial autocorrelation and covariate shift (beyond four discrete quad- rants) is needed. Finally, we look forward to deploying these methods more broadly, ecological, economic applications and beyond.
10 Acknowledgements
We thank John Kilbride for his work in preprocessing and collecting the remote sensing data for the Oregon bird species distribution experiments, Tom Dietterich for comments on an early version the manuscript, and three anonymous reviewers for comments that improved the paper. This research was supported in part by the National Science Foundation (NSF) under Grant No. III-2046678 (JW, RAH), the United States Department of Agriculture National Institute of Food and Agriculture (USDA-NIFA) award No. 2021-67021-35344 (AgAID AI Institute; JW, RAH), the National Acronautics and Space Administration (NASA) under Future Investigators in NASA Earth and Space Science and Technology (FINESST) Grant No. SONSSC2KK664 (LMH), and the Bob and Phyllis Mace professorship (WDR).
References
Courtney L. Amundson, Colleen M. Handel, Daniel R. Ruthrauff, T. Lee Tibbitts, and Robert E. Gill Jr. Data for montane-breeding bird distribution and abundance across national parks of southwestern alaska, 2004-2008. Tn U.S. Geological Survey data release, 2018. doi:
Theodore W Anderson. On the distribution of the two-sample cramer-von mises criterion. The Annals of Mathematical Statistics, 33(3) 1962. doi: 10.1214/a0ms//117700l URL https:///tt...... org/
10.
Sylvain Arlot and Alain Celisse. A survey of cross-validation procodures for model selection. surveys, 4:40 79, 2010. doi: 10.1214/09-SS05)
Statistics
Stephen Bates, Trevor Hastie, and Robert Tibshirani. Cross-validation: what does it estimate and how well does do it? Journal of the American Statistical Association, pp. 1 1-22, 2023. doi: 00.1000/01221 621499.
2023.2197686.
Johan Bjorck, Brendan H. Rappazzo, Qinru Shi, Carric Brown-Lima, Jennifer Dean, Angela Fuller, and Carla Gomes. Accelerating ecological sciences from above: Spatial contrastive learning for remote sens- ing. Proceedings of the AAAT Conference on Arlificial Intelligence, 35(17):14711-443) May 2021. URL https://ojs.. aaai. org/index.aor http:/iee..pp/////aaaaiii...o phtp://At//atticl///eee/////
Alexander Brenning. Spatial cros--validation and bootstrap for the assessment. of prediction rules re- mote sensing: The r package sperrorest. In 2012 IEEE International Geoscience and Remote Sensing Symposium, pp. 5372-5375, 2012. doi: 10.1109 10.1109/IGARSS. /IA///GARSS..0.. 211....21126332332
Prabir Burman, Edmond Chow, and Deborah Nolan. A cross-validator, method for dependent data. Biometrika, 81(2):351-335 358, 1994. doi: 10. httg://iotttttaeeeee (ttommt//11..35.// 18.25.51.
Noel Cressie. Statisties for spatial data. John Wiley & Sons, 2015.
25

---

Published in Transactions on Machine Learning Research (09/2023)
Eric P. Crist and Richard C. Cicone. A physically-based transformation of thematic mapper data the tm tasseled cap. IEEE Transactions on Geoscience and Remole Sensing, GF-22:256-2D 26222 1984.
Saba Dadsetan, Gisele Rose, Naira Hovakimyan, and Jennifer Hobbs. Detection and prediction of nutrient deficiency stress using longitudinal aerial imagery. Proceedings of the AAAT Conference Artificial Intelligence, 35(17):14729-441)
May 2021. URL https:///os aaai.
ptp/AAATTarticlll.
view/17730/17537
Jane Elith and John R. Leathwick. Species distribution models: Ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systemalics, 40(1):677 697, 2009. doi: 10.1146) anmurev.ccolsys. (eeyy..10308.110 110308.120159. URL attps://doi.orrs 11t.::1.1114//mmueeeee/e www.sse.....ol....o..oo. 110308. 1
Allan Gutjahr, Bryan Bullard, and Sean Hatch. General joint conditional simulations using a fast fourier transform method. Mathematical Geology, 2333)336,,,,,, 1997. doi: 10.1007 A007/BF00766641.
Tyler A. Hallman, W. Douglas Robinson, Jenna R. Curtis, and Edward R. Alverson. Building better baseline to estimate 160 years of avian population change and create historically informed conservation targets. Conservation Biology, 35(4):1256-13 2021. ISSN 15231739. doi: 10.1111/cobi 13676
Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning: Data mining, inference, and prediction. Springer, 2001.
Júlio Hoffimann. Geostats.jl high-performanee geostatistics in julia. Journal of Open Source Software, 3
(24):692, Apr. 2018. ISSN 2475-9066. doi: 10.21105/joss.00692 0069.. URL http://dx. http://dx.do..oorggg. wwtddxxdoo...o www. ww..oo..///00/.. http://t......../11......... 21105/
joss. 00692
Júlio Hoffimann, Maciel Zortea, Breno De Carvalho, and Bianca Zadrozny. Geostatistical learning: Chal-
lenges and opportunities. Frontiers in Applied Mathematics and Statistics, 7, 2021. ISSN 2297-4687,
doi: 10.3389/famm /fams. http://aa..222.8888/8 2021.689393. URL https://www.tt http://wwww.rrottirriinnn owg//rtiiee///........ 3389///0............mm.. .000:/ 2021.
689393
Laurel Hopkins, Tyler Hallman, John Kilbride, W Douglas Robinson, and Rebecca A Hutchinson. A
comparison of remotely sensod environmental predictors for avian distributions. Landscape Ecology, 37()): 997-1016, 2022. ISSN 1572-9761. doi: 10.1007//110900-022-01444
URL https:://oo: http::///ooi.rr/1..1110 ort./o..1000000000101...........
610980-002-014066-
Neal Jean, Marshall Burke, Michael Xie, W. Matthew Davis, David B Lobell, and Stefano Ermon. Combining satellite imagery and machine learning to predict poverty. Science, 353((331) 2016. doi: 10.1126/
science.aaf7894.
Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance cstimation. J. Mach. Learn. Res., 1445, 2009.
Alexandre Lacoste, Evan D Sherwin, Hannah R Kerner, Hamed Alemohammad, Björn Lütjens, Jeremy A Trvin, David Dao, Alex Chang, Mehmet Gunturkun, Alexandre Drouin, Pau Rodriguez, and David Vazquez. Toward foundation models for earth monitoring: Proposal for climate change bench- mark. In NewrIPS 2021 Workshop on Tackling Climate Change with Machine Learning, 2021. URL attps://wwe
Kévin Le Rest, David Pinaud, Pascal Monestiez, Joël Chadoeuf, and Vincent Bretagnolle. Spatial leave-one- out crossvaalddation for variable selection in the presence of spatial autocorrelation. Global ecology and biogeography, 23(7):811-82) 2014. doi: 10.1111/geb.121
Mirko Mälicke. Scikit-gstat 1.0: A scipy flavoured geostatistical variogram estimation toolbox written in python. Geoscientific Model Development Discussions, pp. 43, 2021. doi: 10.5184/gmm-22221-174,
Hanna Meyer and Edzer Pebesma. Predicting into unknown space? Estimating the area of applicability of spatial prediction models. Methods in Ecology and Evolution, 12(9):1620-163) 1633, 2021. ISSN 2041210X. doi: 10.1111/2041-200x 13650..
26

---

Published in Transactions on Machine Learning Research (09/2023)
Hanna Meyer, Christoph Reudenbach, Stephan Wöllauer, and Thomas Nauss. Importance of spatial pre- dictor variable selection in machine learning applications - Moving from data reproduction to spatial prediction. Ecological Modelling, (SSSepttmberr:008B1i). 2019. ISSN 03043800. doi: 10.1016/j.ecommodel.
2019.108815. URL https: //tp:://doi.. org/1/......... 1016/ eco.model......ee.....ooo..o.... 2019. 108815
Carles Milå, Jorge Mateu, Edzer Pebesma, and Hanna Meyer. Nearest neighbour distance matching leave-one-out cross-validation for map validation. Methods in Ecology and Evolution, 13(():1304-133) 1316, 2022. doi:
13851. URL https://besjoouraals.co www.orrnnaaas...iileiiirrrrrrrra
wiley. http:.commooaaaa/ aas///aab/00.... 11111111111111111110000 13851
Robert Muscarella, Peter J. Galante, Mariano Soley-Guardia, Robert A. Boria, Jamie M. Kass, María Uriarte, and Robert P. Anderson. ENMeval: An R. package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models. Meth- ods in ecology and evolution, 6(1):119-120, 2014. doi: https://bbl.org//00111111
12261. URL
attps://bessouurnals onlinelibrary wiley wwtp:/wwwwww..l................e wttp:///o//b////..111/////////// 10tp://000.1111111-1111010000000 12221
R. Kelley Pace and Ronald P. Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33: 291 291-297, 1997.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
Steven J Phillips, Miroslav Dudík, and Robert Schapire. A maximum entropy approach to species distri- bution modeling. In Proceedings of the twenty-first international conference on Machine learning, pp. 83, 2004.
Pierre Ploton, Frédéric Mortier, Maxime Réjon-Méchain, Nicolas Barbier, Nicolas Picard, Vivien Rossi, Carsten Dormann, Guillaune Cornu, Gaëlle Viennois, Nicolas Bayol, Alexei Lyapustin, Sylvie Gourlet- Fleury, and Pélissier Raphaël. Spatial validation reveals poor predictive performance of large-scale eco- logical mapping models. Nature communications, 11, 2020. doi: 101108//41447--00--111
Jonne Pohjankukka, Tapio Pahikkala, Paavo Nevalainen, and Jukka Heikkonen. Estimating the prediction performance of spatial models via spatial k-fold cross validation. International Journal of Geographical Information Science, 31110):0001-2019, 2017. doi:
Joaquin Quinonero-Canddela, Masashi Sugiyama, Anton Schwaighofer, and Neil D. Lawrence. Dalasel Shifi in Machine Learning. MIT Press, 2009.
Assaf Rabinowicz and Saharon Rosset. Cross-validation for correlated data. Journal of the American Statistical Association, pp. 14, 2020. doi: 200.//0001114992200.10...
David R Roberts, Volker Bahn, Simone Ciuti, Mark S. Boyce, Jane Elith, Gurutzeta Guillera- Severin Hauenstein, José J. Lahoz-Monfort, Boris Schröder, Wilfried Thuiller, David I. Warton, Brendan A. Wintle, Florian Hartig, and Carsten F. Dormann. Cross-validatiin strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography, 40(8):913 929, 2017. doi: (0.111//ccg.00288/
William D. Robinson, Tyler A. Hallman, and Jenna R. Curtis. Benchmarking the avian diversity of oregon in an era of rapid change. Northwestern Naturalist, 101:180 - 193, 2020.
Georg Ruß and Alexander Brenning. Data mining in precision agriculture: management of spatial in- formation. In International Conference on Information Processing and Management of Uncertainty in Khowledge-Based Systems, pp. 350-359. Springer, 2010. doi:
36.
Eugene Seo, Rebecca A Hutchinson, Xiao Fu, Chelsea Li, Tyler A Hallman, John Kilbride, and W Douglas Robinson. Stateconet: Statistical ecology neural networks for species distribution modeling. Proceedings of the AAAI Conference on Artificial Intelligence, 35(17):513-57) May 2021. URL htttss://js. htops://ojs.aa... aaai. brgginndex.php/AAATT/eae hptp/A/s//arrrraaaaaaa///2222553
27

---

Published in Transactions on Machine Learning Research (09/2023)
Masashi Sugiyama, Matthias Krauledat, and Klaus-Robert Müller. Covariate shift adaptation by importance weighted cross validation. Journal of Machine Learning Research, 8:985-1005, May 2007. URL http: http://////pper//////////aaaaaaa htmml
Waldo Tobler. A computer movie simulating urban growth in the detroit region. Economic Geography, 46: 234-240, 1970.
Roozbeh Valavi, Jane Elith, José J. package for gencrating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods in Ecology and Evolulion, 10(2):225 232, 2018. doi: 10.1111200412100. 13107.
and Gurutzeta Guillera-Arritaa BLOCKCV: An R
Roozbeh Valavi, Gurutzeta Guillera-rrroita, Jose J. Lahoz-Monfort and Jane Elith. Predictive performance of presenco-only species distribution models: a benchmark study with reproducible code. Ecological Mono- gruphs, 2021. doi: 10. 10.1002/ecm) 1486. 1486.
Roozbeh Valavi, Jane Elith, José J Lahoz-Monfort and Gurutzeta Guillera-Arroita Flexible species distribu- tion modelling methods perform well on spatially separated testing data. Global Ecology and Biogeography, 2023.
Alexandre M.J.-C. Wadoux, Gerard B. M. Heuvelink, Sytze de Bruin, and Dick J. Brus. Spatial cross- validation is not the right way to evaluate map accuracy. Ecological Modelling, 457:109692, 2021.
J. Keaton Wilson, Nicolas Casajus, Rebecca A. Hutchinson, Kent P. McFarland, Jeremy T. Kerr, Dominique Berteaux, Maxim Larrivée, and Kathleen L. Prudic. Climate Change and Local Host Availability Drive the Northern Range Boundary in the Rapid Expansion of a Specialist Insect Herbivore, Papilio cresphontes. Frontiers in Ecology and Evolution, (MMrchh) 2021. ISSN 2296701X. doi: 10.3389 2008899fevo.222..55
Makoto Yamada, Taiji Suzuki, Takafumi Kanamori, Hirotaka Hachiya, and Masashi Sugiyama. Relative ransity-ratio estimation for robust distribution comparison. Newral Computation, 25:1324 1370, 2011.
Suwei Yang, Massimo Lupascu, and Kuldeep S. Meel. Predicting forest. fire using remote sensing data and machine learning. Proceedings of the AAAI Conference on Artificial Intelligence, 35(17): May 2021. URL https://oja www.aaar..aaaaaaaaaaaaaaaa htsp/indexx.coe http://wwwwppp..a/tttttaaaa/////
Damaris Zurell, Niklaus E Zimmermann, Helge Gross, Andri Baltensweiler, Thomas Sattler, and Rafacl O Wüest.. Testing species assemblage predictions from stacked and joint species distribution models. Journal of Biogeography, 47(1):101-113, 2020. doi: 10. (0.11111jbbi.331
28

---

Published in Transactions on Machine Learning Research (09/2023)
Appendix
A Implementation Details
A.1 Additional Simulation Details
We simulated Gaussian random fields via the fast Fourier transform method (Gutjahr et al. 1997 Hoffimann 2018). Tab. summarizes the joint distributions of the two spatial features for training and test sets, and the sampling strategies we use to extract training and test data points in each scenario.
A.2 Density Ratio Estimation Details
We applied the Relative unconstrained Least-Squares Importance Fitting (RuLSIF) method is applied to es- timate density ratios (amaada et al. 2011]. Given the definition, set alpha 0, and the mumber of kernels
= 50. The algorithm automatically selects an optimal value pair in the range of sigma = [0.1, 0.5, 1, 2, 3, 4) and lambda = [0.1, 0.3, 0.5, 0.8, 1]. We fitted Matérn variogram functions with the lag class estimated by Scott's 's rule Mälicke 2021 to calculate ranges of the features of training sets.
29

---

Published in Transactions on Machine Learning Research (09/2023)
1.08
Error estimate
1.06
1.04
1.02
1.00
0.98
2
4
8
12
20
Block size
SAC Range:
12
Figure 9: The average BLCV estimates across 100 simulations in Scenario SD, ST, ST CSi, ST + CSp.
Table 10: Simulation: BLCV estimates with block sizes 2, 4, 8, 12, 20 grids. Best estimales in each column in bold.
Scenario SD
Scenario SI
4
8
12
4
8
12
Risk
0.9912
0.9888
0.9823
0.9994
1.0170
1.0425
BLCV2
1.0053
0.9918
0.9851
1.0053
0.9918
0.9851
BLCV4
1.0088
0.9987
0.9899
1.0088
0.9987
0.9899
BLCV8
1.0137
1.0222
1.0331
1.0137
1.0222
1.0331
BLCV12
1.0164
1.0313
1.0601
1.0164
1.0313
1.0601
BLCV20
1.0167
1.0415
1.0900
1.0167
1.0415
1.0900
Scenario SI + CSi
Scenario SI + CSp
r
4
8
12
8
12
Risk
0.5134
0.5267
0.5384
2.0107
2.0046
1.9702
BLCV2
1.0053
0.9918
0.9851
1.0053
0.9918
0.9851
BLCV4
1.0088
0.9987
0.9899
1.0088
0.9987
0.9899
BLCV8
1.0137
1.0222
1.0331
1.0137
1.0222
1.0331
BLCV12
1.0164
1.0313
1.0601
1.0164
1.0313
1.0601
BLCV20
1.0167
1.0415
1.0900
1.0167
1.0415
1.0900
B Additional Results
B.1 Effects of Block Size in Simulation
When there was spatial autocorrelation (i.e., r 1 ) in the simulated datascts, we fine-tuned the buffer size
= 2, 4, 8, 12, 20 grids. In general, the error estimate (the magnitude, not the bias) increased with increasing block size (Fig. 9 and Tab. 10p. This finding is consistent with the proposition that covariate shift is increased with larger blocks.
30

---

Published in Transactions on Machine Learning Research (09/2023)
1.5
1.5
1.0
1.0
0.5
0.5
Bias
Bias
0.0
0.0
-0.5
-0.5
-1.0
-1.0
-1.5
-1.5
4
8
12
4
8
12
Variogram range (a) Scenario SI + RS
Variogram range (b) Scenario SI + CSp
CV algorithm:
KFCV
BLCV
BFCV
IWCV
IBCV
Figure 10: Biases of CV estimates in scenarios with various characteristics: spatially independent with range shift (ST RS), spatially independent with covariate shift (SI CSp: partially overlapped). Circles inside the boxes display the mean values of biases. The black dash lines illustrate no bias. Both show bias of the average CV estimate to the risk.
Table 11: Risk perspective: Bias of the average CV estimate (over 9 folds) of risk. The least biased estimates in each row are bold.
CV Algo.
KFCV BLCV
BFCV
IWCV
IBCV
KFCV
BLCV
BFCV
IWCV
IBCV
r
Scenario +
Scenario + CSp
1
-0.0523
0.0250
0.1117
-0.1365
0.0086
-0.9935
-0.9934
-0.9933
-0.7157
-0.7157
4
-0.0463
0.0328
0.1143
-0.1256
0.0197
-1.0090
-1.0019
-0.9914
-0.7714
-0.7497
8
-0.0633
0.0151
0.0994
-0.1390
0.0093
-1.0177
-0.9824
-0.9460
-0.0024
-0.7121
12
-0.9889
-0.9101
-0.8277
-0.7805
-0.5922
B.2 Special Cases of SI + CS Simulation Experiments
We treated Scenario ST + RS, where training and test sets are spatially independent with different ranges, and Scenario SI + CSp, where training and test sets are spatially independent and their feature distribution are partially overlapped, special cases of Scenario SI + CS in our simulation experiments. Results for these scenarios are shown in Fig. 10 and listed in Tabs. and 12
Table 12: Test error perspective: Proportion of simulations in which a certain CV algorithm is the closest estimate to the test error. Highest numbers the average row bold.
KFCV
BLOV
BFCV
IWCV
IBCV
KFCV
BLCV
BFCV
IWCV
IBCV
r
Scenario SI + RS
Scenario SI CSp:
1
0.29
0.23
0.25
0.16
0.07
0.00
0.00
0.01
0.56
0.43
4
0.22
0.20
0.27
0.20
0.12
0.00
0.00
0.01
0.13
0.86
8
0.08
0.29
0.32
0.25
0.07
0.00
0.02
0.05
0.04
0.89
12
0.00
0.00
0.16
0.05
0.79
31

---

Published in Transactions on Machine Learning Research (09/2023)
B.3
Feature Densities for the Empirical Datasets
Fig.
1
16
illustrate the distributions of training and test features.
4
10.0
3
Densiy
Density
Density
7.5
Density
2
5.0
1
2.5
2.5
0.0
0.00
0.25
0.50
0.75
0.2
0.4
0.6
0.0
0.1
0.3
-0.3
-0.2
-0.1
TOW
training
test
Figure 11: HEWA1000 SI + CS: probability density estimation by kernel smoothing of four features.
4
10.0
10.0
Density
Density
6
Density
7.5
Density
7.5
5.0
2.5
0
0.0
0.0
0.00
0.25
0.50
0.75
0.2
0.4
0.6
0.0
0.1
0.2
0.3
-0.3
-0,2
D..1
D.0
TCA
TCB
TCG
TOW
training
test
Figure 12: HEWA1800 SI + CS: probability density estimation by kernel smoothing of four features.
3
Density
Density
Density
7.5
Density
6
5.0
2
0.25
ICA
TCB
TCW
training
test
Figure 13: WETA1800 SI + CS: probability density estimation by kernel smoothing of four features.
32

---

Published in Transactions on Machine Learning Research (09/2023)
2.0
Density
Density
Density P
2-
1.0
0.5
0.0
0.00
0.25
0,50
0.75
1.00
0.00
0.25
0,50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
Lousing_mddia__Pgg
total_rooms
totro_bedrooms
Density
Density
Density
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
population
households
median income
_ecdian_income
training
test
Figure 15: Housing (a) by Latitude: probability density estimation by kernel smoothing of features.
25
2.0
1.0
20
Donaty
Deneity
Donaty
0.8
Denaly
0.6
1.0
10
0.4
10
02
0.0
0.0
0.2
0.4
0.6
0.6
0.8
0..
0.0
1.0
1.5
0.5 dshrub
1.0
water
wetland
shrub
10
15.0
12.5
Donalty
10.0
Dereity
Dersity
7.5
5.0
2.5
0.0
spruce
baresnow
training
test
Figure 14: Alaskan birds: probability density estimation by kernel smoothing of eight leatures. Features dec, mixed and spruce of test samples are all zeros so their density curves are skipped.
33

---

Published in Transactions on Machine Learning Research (09/2023)
1.5
Density
Density
Density
1.0
2
0.0
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
housing median age
rooms
bedrooms
2.5
2.5
Density
Density
2.0
1.5
1.0
0.5
0.0
0.0
0.25
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.75 median_ _icilaemiccoeeeee
0.50
population
households
training
test
Figure 16: Housing (b) by Bay: probability density estimation by kernel smoothing of six features.
34

---

Published in Transactions on Machine Learning Research (09/2023)
Table 13: HEWA1000: model classification test error rates (targets) and 9-fold CV estimates thereof (best estimates in each column bold). BLCV-best, BFCV-best and IBCV-best estimates are selected from the best ones in the Appx. Tab. 15 for peak-to-peak comparison. BLCV-range, BFCV-range, and IBCV- range set the tuning paramcters priori based the range of spatial autocorrelation in the features. A dash means that setting the hyperparameters based on the range gives the best value (i.e., the methods are equivalent).
Model
Test error (target)
KFCV
IWCV
BLCV -range
BFCV -range
IBCV -range
BLCV -best
BFCV
IBCV -best
-best
SD
Ridge
0.1700
0.1709
0.1706
0.1664
0.1783
0.1780
LSVM
1720
0.1709
0.1706
0.1678
0.1775
0.1773
KNN
0.1740
0.1779
0.1777
0.1905
0.1909
0.1906
RF
1710
0.1910
0.1907
0.1964
0.1989
0.1986
NB
0.1700
0.1729
0.1727
0.1663
0.1678
0.1676
ST
Ridge
0.2320
0.1890
0.1888
0.2007
0.2173
0.2170
0.2357
0.2411
0.2408
LSVM
0.2280
0.1900
0.1898
0.2047
0.2876
0.2872
0.2406
0.2396
0.2394
KNN
0.2440
0.2180
0.2178
0.2292
0.2628
0,2625
0.2402
0.2383
0.2381
RF
0.2520
0.2120
0.2118
0.2607
0.2692
0.2689
0.2365
0.2363
NB
0.2440
0.1950
0.1948
0.2043
0.2110
0.2107
0.2424
0.2441
0.2438
SD + CS
Ridge
0.2140
0.1709
0.2400
0.1644
0.1783
0.2526
0.1976
0.1997
LSVM
0.2040
0.1709
0.2430
0.1678
0.1775
0.2489
0.2055
0.2035
KNN
0.2080
0.1779
0.2533
0.1905
0.1909
0.2649
0.2186
RF
0. 0.1840
0.1910
0.2706
0.1964
0.1989
0.2789
NB
0.2160
0.1729
0.2469
0.1663
0.1678
0.2393
0.1984
0.2037
SI
CS
Ridge
0.2420
0.1709
0.2239
0.1644
0.1783
0.2370
0.1976
0.1997
0.2384
LSVM
0.2540
0.1709
0.2245
0.1678
0.1775
0.2334
0.2055
0.2035
0.2544
KNN
0.2440
0.1779
0.2370
0.1905
0.1909
0.2484
0.2425
0.2406
RF 0.2540
0.1910
0.2507
0.1964
0.1989
0.2572
0.2322
0.2453
NB
0.2640
0.1729
0.2289
0.1663
0.1678
0,2222
0.1984
0.2037
0.2546
B.4 Additional Results from the Oregon 2020 Study
Tabs. and 14 show the analagous Oregon 2020 results from the main text, but for the HEWA1000 and
WETA1800 datasets, respectively.
35

---

Published in Transactions on Machine Learning Research (09/202))
Table 14: WETAI100: model classification test error rates (targets) and 9-fold CV estimates thereof (best estimates in each column in bold). BLCV-best, BFCV-best and IBCV-best estimates are selected from the best ones in the Appx. Tab. 15 for a peak-to-peak comparison. BLCVVrange, BFCV-range, and TBCV- range set the tuning paramcters a priori based on the range of spatial autocorrelation in the features. A dash means that setting the hyperparameters based on the range gives the best value (i.e., the methods are equivalent).
Model
Test error (target)
KFCV
IWCV
BLCV -range
BFCV -range
IBCV
BLCV -best
BFCV -best
IBCV -best
-range
SD
Ridge
0.3560
0.3850
0.3828
0.3915
0.3986
0.3963
0.3845
0.3885
0.3863
LSVM
0.3600
0.3822
0.3800
0.3862
0.3896
0.3873
0.3819
0.3799
0.3777
KNN
0.4040
0.4028
0.4004
0.4057
0.4094
0.4070
RF
0.4020
0.4022
0.3999
0.4020
0.4133
0.4109
0.4105
0.4082
NB
0.3620
0.3783
0.3762
0.3815
0.3788
0.3766
SI
Ridge
0.4320
0.4050
0.4043
0.4503
0.4810
0,4802
0.4151
0.4382
0.4375
LSVM
0.4120
0.3894
0.3888
0.4416
0.5029
0.5020
0.4050
0.4226
0.4219
KNN
0,4380
0.4011
0.4004
0.4220
0.4279
0.4271
0.4245
RF
0.4660
0.3983
0.3976
0.4205
0.4312
0.4304
NB
0.4240
0.3872
0.3866
0.4275
0.4268
0.4260
SD + CS
Ridge
0.3980
0.3850
0.3828
0.3915
0.3986
0.3963
LSVM
0.3960
0.3822
0.3801
0.3862
0.3896
0.3873
KNN
0.3680
0.4028
0.4005
0.4057
0.4094
0.4070
0.3846
RF
0.3320
0.4022
0.3999
0.4020
0.4133
0.4109
0.3790
0.4105
0.4082
NB
0.4020
0.3783
0.3762
0.3815
0.3788
0.3767
0.3923
0.3941
0.3920
SI + CS
Ridge
0.5100
0.3850
0.4299
0.3915
0.3986
0.4458
0.4253
0.4435
0.4872
LSVM
0.5100
0.3822
0.4255
0.3862
0.3896
0.4278
0.4272
0.4408
0.4803
KNN
0.4400
0.4028
0.4410
0.4057
0.4094
0.4398
0.4199
RF
0.4280
0.4022
0.4389
0.4020
0.4133
0.4451
0.4120
0.4212
NB
0.4800
0.3783
0.4175
0.3815
0.3788
0.4104
0.4330
0.4318
0.4868
36

---

Published in Transactions on Machine Learning Research (09/2023)
B.5
Empirical Results by Block Sizes
Tab.
17 show the CV estimates for each hyperparameter block size setting. We pick the best block size given the target risks or test errors although it is impossible to get these target model errors in real applications. Here we are more interested in the best estimates of each CV algorithm and their peak-to-peak comparisons.
10
5
16
Table 15: Oregon's birds: BLCV, BFCV and IBCV estimates with different block sizes. Best estimates of each CV algorithm in cach column bold.
Classifier
Ridge
LSVM
KNN
RF
NB
Scenario
SD, HEWA1000
Test error
0.1700
0.1720
0.1740
0.1740
0.1700
BLCV0.28
0.1644
0.1678
0.1905
0.1964
0.1663
BLCV0.56
0.1976
0.2055
0.2425
0.2272
0.1984
BLCV0.84
0.1868
0.1867
0.2186
0.2322
0.1808
BFCV0.28
0.1783
0.1775
0.1909
0.1989
0.1678
BFCVV0.66
0.1997
0.2035
0.2406
0,2453
0.2037
BFCV0.84
0.1889
0.1902
0.2235
0.2159
0.1848
IBCV0.28
0.1780
0.1773
0.1906
0.1986
0.1676
IBCV0.56
0.1994
0.2032
0.2403
0.2449
0.2034
IBCV0.84
0.1886
0.1899
0.2231
0,2156
0.1845
Scenario
SD, KIEWA1800
Test error
0.1640
0.1640
0,2020
0,2040
0.1720
BLCV0.30
0.1786
0.1752
0.2084
0.2057
0.1810
BLCV0.60
0.2081
0.2111
0.2439
0.2162
0.2131
BLCV0.90
0.2083
0.2092
0,2512
0,2438
0.2098
BFCV0.30
0.1780
0.1842
0.2035
0.1945
0.1757
BFCV0.60
0.2120
0.2148
0.2494
0.2446
0.2114
BFCV0.90
0.2072
0.2081
0.2458
0.2472
0.2084
IBCV0.30
0.1781
0.1846
0.2031
0.1938
0.1757
IBCV0.60
0.2107
0.2135
0.2478
0.2431
0.2099
IBCV0.90
0.2057
0.2066
0.2439
0.2440
0.2070
Scenario
SD, WETA1800
Test error
0.3560
0.3600
0.4040
0.4020
0.3620
BLCV0.27
0.3915
0.3862
0.4057
0.4020
0.3815
BLCVO.54
0.3845
0.3819
0.3846
0.3790
0.3923
BLCV0.81
0.4253
0.4272
0.3915
0.4120
0.4330
BFCV0.27
0.3986
0.3896
0.4094
0.4133
0.3788
BECV0.54
0.3885
0.3799
0.4199
0.4105
0.3941
BFCV0.81
0.4435
0.4408
0.4160
0.4212
0.4318
IBCV0.27
0.3963
0.3873
0.4070
0.4109
0.3766
IBCV0.54
0.3863
0.3777
0.4175
0.4082
0.3919
IBCV0.81
0.4409
0.4383
0.4136
0.4189
0.4294
Scenario
SI, HEWA
Test error
0.2320
0.2280
0.2440
0.2520
0.2440
BLCV0.33
0.2007
0.2047
0.2292
0.2607
0.2043
BLCV0.66
0.2357
0.2406
0.2402
0.2343
0.2424
BFCV0.33
0.2173
0.2876
0.2628
0,2692
0.2110
BFCV0.66
0.2411
0.2396
0.2383
0.2365
0.2441
IBCV0.33
0.2170
0.2872
0.2625
0.2689
0.2107
IBCV0.66
0.2408
0.2394
0.2381
0.2363
0.2438
Scenario SI, HEWA1800
Test error
0.2000
0.2040
0.2260
0.2580
0.2140
37

---

Published in Transactions on Machine Learning Research (09/2023)
BLCVO.32
0.1711
0.1805
0.1753
0.1700
0.1892
BLCV0064
0.2711
0.2628
0.2117
0.2215
0.2701
BFCV0.32
0.1842
0.1989
0.1697
0.1950
0.1853
BFCV0.64
0.2733
0.2694
0.2028
0.2114
0.2735
IBCV0.32
0.1839
0.1986
0.1695
0.1947
0.1850
IBCV0.64
0.2730
0.2691
0.2025
0.2111
0.2732
Scenario
SI, WETA1800
Test error
0.4320
0.4120
0.4380
0.4660
0.4240
BLCV0.28
0.4503
0.4416
0.4220
0.4205
0.4275
BLCV0.56
0.4851
0.4648
0.4245
0.4135
0.4611
BLCV0.84
0.4151
0.4050
0.3691
0.3773
0.3819
BFCV0.28
0.4810
0.5029
0.4279
0.4312
0.4268
BFCV0.56
0.5221
0.4917
0.4167
0.3953
0.4653
BFCV0.84
0.4382
0.4226
0.3836
0.3598
0.3895
IBCV0.25
0.4802
0.5020
0.4271
0.4304
0.4260
IBCV0.56
0.5212
0.4909
0.4161
0.3946
0.4646
IBCV0.84
0.4375
0.4219
0.3829
0.3592
0.3889
Scenario SD + CS,
HEWA
Test error
0.2140
0.2040
0.2080
0.1840
0.2160
BLCV0.28
0.1644
0.1678
0.1905
0.1964
0.1663
BLCV0.56
0.1976
0.2055
0.2425
0.2272
0.1984
BLCV0.84
0.1868
0.1867
0.2186
0.2322
0.1808
BFCV0.28
0.1783
0.1775
0.1909
0.1989
0.1678
BFCV0.56
0.1997
0.2035
0.2406
0.2453
0.2037
BFCV0.84
0.1889
0.1902
0.2235
0.2159
0.1848
IBCV0.28
0.2526
0.2489
0.2649
0.2789
0.2393
IBCV0.56
0.2724
0.2767
0.3246
0.3327
0.2791
IBCV0.84
0.2649
0.2669
0.3085
0.2991
0.2614
Scenario SD + CS,
HEWA1800
Test error
0.2140
0.2040
0.2080
0.1840
0.2160
BLCV0.30
0.1786
0.1752
0.2084
0.2057
0.1810
BLCV0.60
0.2081
0,2111
0.2439
0.2162
0.2131
BLCV0.90
0.2083
0.2092
0.2512
0.2438
0.2098
BFCV0.30
0.1780
0.1842
0.2035
0.1945
0.1757
BFCV0.60
0.2120
0.2148
0.2494
0.2446
0.2114
BFCVO.90
0.2072
0.2081
0.2458
0.2472
0.2084
IBCV0.30
0.1778
0.1840
0.2032
0.1942
0.1755
IBCV0.60
0.2117
0.2145
0.2491
0.2443
0.2111
IBCV0.90
0.2069
0.2078
0.2454
0.2468
0.2081
Scenario
+ CS,
WETA1800
Test error
0.3980
0.3960
0.3680
0.3320
0.4020
BLCV0.27
0.3915
0.3862
0.4057
0.4020
0.3815
BLCV0.54
0.3845
0.3819
0.3846
0.3790
0.3923
BLCV0.81
0.4253
0.4272
0.3915
0.4120
0.4330
BFCV0.27
0.3986
0.3896
0.4094
0.4133
0.3788
BFCV0.54
0.3885
0.3799
0.4199
0.4105
0.3941
BFCV0.81
0.4435
0.4408
0.4160
0.4212
0.4318
IBCV0.27
0.3963
0.3873
0.4070
0.4109
0.3767
IBCV0.54
0.3863
0.3778
0.4175
0.4082
0.3920
IBCV0.81
0.4410
0.4384
0.4137
0.4189
0.4295
Scenario SI
CS,
H0WA1000
Test error
0.2420
0.2540
0.2440
0.2540
0.2640
BLCV0.28
0.1644
0.1678
0.1905
0.1964
0.1663
BLCV0.56
0.1976
0.2055
0.2425
0.2272
0.1984
38

---

Published in Transactions on Machine Learning Research (09/2023)
BLCV0.84
0.1868
0.1867
0.2186
0.2322
0.1808
BFCV0.28
0.1783
0.1775
0.1909
0.1989
0.1678
BFCV0.56
0.1997
0.2035
0.2406
0.2453
0.2037
BFCV0.84
0.1889
0.1902
0.2235
0.2159
0.1848
IBCV0.28
0.2370
0.2334
0.2484
0.2572
0.2222
IBCV0.56
0.2494
0.2544
0.2994
0.3093
0.2546
IBCV0.84
0.2384
0.2405
0.2818
0.2746
0.2343
Scenario ST
CS,
HEWA 1800
Test error
0.2400
0.2460
0.2420
0.2380
0.2660
BLCV0.30
0.1786
0.1752
0.2084
0.2057
0.1810
BLCVV.60
0.2081
0.2111
0.2439
0.2162
0.2131
BLCV0.90
0.2083
0.2092
0.2512
0.2438
0.2098
BFCV0.30
0.1780
0.1842
0.2035
0.1945
0.1757
BFCV0.60
0.2120
0.2148
0.2494
0.2446
0.2114
BFCV0.90
0.2072
0.2081
0.2458
0.2472
0.2084
IBCV0.30
0.2085
0.2167
0.2370
0.2239
0.2059
IBCV0.60
0.2397
0.2421
0,2833
0.2734
0.2385
IBCV0.90
0.2394
0.2385
0.2799
0.2744
0.2403
Scenario SI + CS,
WETA1800
Test error
0.5100
0.5100
0.4400
0.4280
0.4800
BLCV0.27
0.3915
0.3862
0.4057
0.4020
0.3815
BLCV0.54
0.3845
0.3819
0.3846
0.3790
0.3923
BLCV0.81
0.4253
0.4272
0.3915
0.4120
0.4330
BFCV0.27
0.3986
0.3896
0.4094
0.4133
0.3788
BFCV0.54
0.3885
0.3799
0.4199
0.4105
0.3941
BFCV0.81
0.4435
0.4408
0.4160
0.4212
0.4318
IBCV0.27
0.4458
0.4278
0.4398
0.4451
0.4104
IBCV0.54
0.4469
0.4304
0.4750
0.4680
0.4510
IBCVO.81
0.4872
0.4803
0.4690
0.4756
0.4868
Table 16: Alaskan birds: BLCV, BFCV and IBCV estimates with block sizes = 0.82, 0.9, 1.0, 1.1 degrees. Best estimates of each CV algorithm in each column in bold.
Classifier
Ridge
LSVM
KNN
RF
NB
Test error
0.1866
0.1866
0.3284
0.3657
0.4030
BLCV0.82
0.2777
0.2723
0.3030
0.3176
0.3230
BLCVO.9
0.2696
0.2714
0.2896
0.3460
0.3193
BLCV1.0
0.3119
0.3081
0.3279
0.3252
0.3522
BLCV1.1
0.2945
0.3117
0.3216
0.3268
0.3351
BFCV0.82
0.2881
0.2791
0.3251
0.3343
0.3597
BFCV0.9
0.2555
0.3252
0.2764
0.3424
0.3332
BFCVL.0
0.3122
0.3088
0.3582
0.3554
0.3279
BFCV1.1
0.3142
0.3464
0.3713
0.3203
0.3663
IBCV0.82
0.2374
0.2244
0.2542
0.2998
0.2748
IBCV0.9
0.2444
0.2391
0.3254
0.3483
0.3667
IBCV1.
0.2485
0.2594
0.3118
0.3021
0.3031
IBCV1.1
0.2345
0.2456
0.2882
0.2800
0.2782
39

---

Published in Transactions on Machine Learning Research (09/2023)
Table 17: Housing: BLCV, BFCV and IBCV estimates with block sizes 0.23, 0.25, 0.5, 1.0, 1.5, 2.0, degrees. Best estimates of each CV algorithm in each column in bold.
Regressor
Lincar
KRR
SVR
KNN
RF
(a)
By Latitude
Test error
0.5134
0.5931
0.3469
0.5246
0.4771
BLCV0.23
0.7024
0.7135
0.6430
0.7132
0.6728
BLCV0.5
0.7561
0.7589
0.7049
0.7731
0.7269
BLCV1.0
0.7461
0.7643
0.6784
0.7586
0.7185
BLCVL5
0.7938
0.8018
0.7296
0.8151
0.8018
BLCV2.0
0.7153
0.7270
0.6671
0.7278
0.7129
BLCV2.5
0.7419
0.7656
0.6847
0.7765
0.7557
BFCV0.23
0.7062
0.7184
0.6564
0.7389
0.7010
BFCV0.5
0.7793
0.7854
0.7285
0.8042
0.7614
BFCV10
0.7537
0.7741
0.6885
0.7757
0.7346
BFCVL.5
0.8062
0.8125
0.7420
0.8344
0.8125
BFCV2.0
0.7201
0.7323
0.6714
0.7370
0.7189
BFCV2.5
0.7435
0.7671
0.6862
0.7798
0.7557
ICV00.23
0.4286
0.4332
0.4178
0.4558
0.4321
IBCV0.5
0.5532
0.5479
0.5504
0.5955
0.5590
IBCV1.0
0.4896
0.5094
0.4650
0.5113
0.4877
IBCV1.5
0.6938
0.6999
0.6651
0.6777
0.6824
IBCV2.0
0.5269
0.5594
0.5284
0.5358
0.5596
IBCV2.5
0.5317
0.5584
0.5094
0.5510
0.5471
(b) By Bay
Test error
0.7441
0.7813
0.6933
0.7830
0.7159
BLCVO.25
0.6900
0.7006
0.6262
0.6960
0.6510
BLCV0.5
0.7053
0.7123
0.6407
0.7143
0.6702
BLCV1.0
0.6906
0.6960
0.6139
0.6946
0.6563
BLCV155
0.6884
0.6980
0.6258
0.7045
0.6640
BLCV2.0
0.6618
0.6758
0.5861
0.6670
0.6258
BLCV2.5
0.7170
0.7316
0.6216
0.6998
0.6687
BFCV0.25
0.7199
0.7309
0.6592
0.7430
0.6970
BFCV0.5
0.7330
0.7394
0.6687
0.7530
0.7050
BFCV1.0
0.6988
0.7037
0.6239
0.7095
0.6627
BFCV1.5
0.6980
0.7095
0.6349
0.7177
0.6749
BFCV2.0
0.6676
0.6816
0.5920
0.6754
0.6339
BFCV2.5
0.7284
0.7391
0.6350
0.7146
0.6798
IBCV0.25
0.8518
0.8722
0.7744
0.8821
0.8196
IBCV0.5
0.8529
0.8671
0.7785
0.8790
0.8199
IBCVL.0
0.7969
0.8070
0.7076
0.8081
0.7502
IBCV1.5
0.7959
0.8161
0.7216
0.8212
0.7678
IBCV2.0
0.7793
0.7988
0.6838
0.7859
0.7336
IBCV2.5
0.8359
0.8507
0.7191
0.8088
0.7663
40


Reviews:
[{"review_id": "KDCsTu02ks", "reviewer": "Reviewer_nnud", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n*The IBCV procedure is interesting and provides a useful tool for geospatial cross-validation in particular for the SI+CS condition.\n\n*The experimentation is extensive providing thorough synthetic and real data experimentation.\n\nWeaknesses \n\n* The proposed IBCV is rather straightforward combining IWCV and BFCV to do both. As such, the technical innovation of the paper is rather limited.\n\n* The procedure relies on the tuning of hyperparameters and their careful tuning appear to be a challenge.\n\n* The overall recommendations and results are somewhat non-surprising making the contribution limited although well executed and extensive.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper considers cross-validation in the context of geospatial data in which cross-validation procedures are challenged by spatial dependence and covariate shifts between training and test splits violating the common assumption of i.i.d. (independence and identically distributed) observations. The authors propose the IBCV and importance weighted buffered approach to CV and make extensive evaluations of various geospatial CV procedures on synthetically generated and real datasets. Apart from introducing the IBCV procedure the key contribution is extensive comparison of CV procedures systematically evaluating spatial dependency as opposed to spatial independence (SD vs SI) and covariate shift as opposed to no covariate shift (no CS vs. CS) in which the IBCV is favorable in the SI+CS scenario. Spatial dependence is assessed using spatial autocorrelation in training and test data whereas covariate shift is assessed using a non-parametric Cramér test. Importantly, the extensive analyses provided are used to make general recommendations in regards to best practices for geospatial cross-validation.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "It would be good to provide a clear recommendation how parameters should be tuned and perhaps also provide the reader with a recipe for recommendations according to the four scenarios SD/SI, no CS/CS in a box summarizing the key take-aways and guidelines of the paper.\n\nThe proof of theorem 1 is very condensed and difficult to follow. I suggest the steps in the derivation be elaborated and terms combined explicitly to derive P(T_te|T_Tr) and p(T_k|T_-k).\n\nThe importance weighting relies on some density estimation procedure for p_tr(x) and p_te(x). Presently the authors consider Relative unconstrained Least-Squares Importance Fitting (RuLSIF) to quantify these density ratios. It would be good to further discuss this choice as opposed to many alternatives providing density estimates as well as the impact of choice of procedure here chosen which I would expect can substantially impact results.\n\nIt is very unclear to me how the target test error is known for the real datasets considered as this has to be quantified by some cross-validation means as the generative process, noise structure etc. is unknown. Please clarify how the test error target value is quantified and known.\n\n\nMinor comment: The paper in general reads well but please correct the following sentence:\nresponse variables are the also same-> response variables are also the same\n\n", "broader_impact_concerns": "There is in general no major ethical concerns of the executed research. The authors could emphasize on the importance of correct performance assessment not leading to erroneous conclusions, however, they do a good job in the introduction discussing and establishing this, while in the discussion discussing the broader importance of geo-spatial cross-validation.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Hftpb3ZFyD", "reviewer": "Reviewer_T2Qn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- Addresses an important and non-trivial problem and delivers concrete algorithms for practical use\n- Very comprehensive experimentation\n- Gives concrete recommendations for practitioners\n\nWeaknesses\n - Overly long\n - The new algorithm for the SI+CS scenario is fairly straightforward generalisation of previous methods\n - Paper mentions code is provided in the Supplement, but I could not find it. The value of the paper somewhat depends on the quality of the code, since robust and reliable tools are needed for practitioners to actually use these techniques", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The manuscript provides an overview of challenges in using cross-validation for geospatial problems, characterising different conditions systematically and explaining what kind of algorithms can be used in condition to learn unbiased risk estimates. In addition, the authors provide a new cross-validation variant applicable for one of the cases for which no previous unbiased methods were available. The authors provide comprehensive experimental evaluation of the methods on both synthetic and real data. It is likely that the paper will be interesting for practitioners needing CV for geospatial applications, and further theoretical development building on this is also possible.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The paper is overall well written and easy to read, but I do not think the length (21.5 pages before references) is properly justified by the content. It would probably be difficult to compress this all the way down to the recommended 12 pages, but something like 16 pages would make the paper more readable and could be achieved by moving some of the experimental material to a supplement and by compressing the verbose presentation style. I admit this is to an extent a personal preference, so this is not to interpreted as direct request for changing the paper.\n\nI have some trouble following Theorem 1 and would like to see it clarified. The main claim sounds reasonable, but it is difficult to follow what exactly happens in the proof since you are multiplying both sides with different constants (that I presume are actually the same under your assumptions) and are quite liberally introducing new terms. You should at least justify what happened on the second and third lines of the three-line equation right below Eq. (3) since especially how the term $p(y_{tr}|X_{te},X_{tr},y_{te})$ appears is not really obvious. Overall, I am having a bit of trouble understanding why the first part of the proof is even needed -- I would have been perfectly happy to start from Eq. (4) as the underlying assumption, and it feels like a more straightforward proof could also be done. If $p(X|X)$ and $p(y|X)$ are both the same in the two cases then quite obviously the same holds for the data sets.\n\nThe characterisation of the four different cases in Table 1 is nice and likely to be interesting as a general guide for practitioners, but I find hard thresholding based on a p-value quite dangerous since p-values are not really that reliable indicators. Wouldn't it make more sense for a practitioner to try out both methods when the p-value is close to the threshold? Some discussion on the possible problems of p-value thresholding should be added.", "broader_impact_concerns": "No concerns. The paper makes contributions that reduce the risk of incorrect modelling choices and has only positive impacts.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bi93lIx1Ht", "reviewer": "Reviewer_SgZ4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I think the topic is very important, and I am excited about the goals outlined in the paper.  \nFrom my view the literature review is well done, and helps provide useful perspective for the paper.  \nThe attempts at modeling this through simulations and on real data are fairly well done, and support the conclusions and help justify the motivation.  \n\nAs such I am supportive of accepting this paper in TMLR.  \n\nSome weakness are described via the suggestions below.  ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper provides a framework for thinking about how to perform cross validation on spatial regression problems.  Such problems aim to predict a quantity based on a combination of spatial location and feature values observed at that location.  An example application considered is species distribution modeling where one predicts the spatial density of a species (e.g., birds) from point-wise observations.  In these scenarios, observations are spatially correlated (two nearby observations are likely to record similar values) so if observations are dense enough a purely spatial model (e.g., based on Nadaraya-Watson Kernel Regression or K-Nearest Neighbor) could be used; but this is not the goal, the goal is often to understand the relationship between the features and the prediction.  This means if one is to just use standard cross-validation, on NWKR model may do very well without even observing the features.  \nTo address this, the community has proposed a number of methods to split data for cross validation via spatial blocks or with spatial buffers.  This paper attempts to explain and unify these approaches, provide guidance on when to use each, and introduces a new combination of ideas based on both importance weighting and buffering, that works best in some scenarios.  ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1.  It is written for a more statistical audience, less for a ML one.  However, similar problems and challenges occur in many subareas of ML.  A few that come to mind (in additional to spatial problems) are learning on temporal data and scientific ML where the parameter space acts like spatial parameters.  \n\n2.  The main theorem and derivation is hard to follow, and I think not given in enough detail.  I am particularly thinking about the derivation right after equation (3).  I could not follow all steps.  More steps could be added, and they should each be explained.  If this is the central technical result, it should be explicitly clear.  \n\n3.  The discrete decomposition into 4 cases (Table 1) is useful for presentation, but as discussed many other places, these 4 cases are not so clear.  The use of a hypothesis test with user-required parameters seems problematic for making this rigorous.  It is suggested these sorts of tests could be used to algorithmically decide what type of cross-validation is warranted.  I would either (a) note that these tests are only for exposition as they are introduced, or (b) lean into the usefulness of the tests, and argue that they *should* be used for choosing how to perform cross-validation.  \n\n4.  The very heavy use of acronyms made the paper hard to read, especially for a reader who had not used similar ones frequently in the past.  I suggest to more times in the paper remind the reader what they stand for.  \n\n5.  I took me several reads of the paper to really understand the issue at hand, and how the feature values and spatial values were meant to interact in the prediction.  Not until I looked closely at the bird example did I start to understand this (I think I described it correctly above).  I suggest the authors re-work the introductory settings to describe these issues, and better highlight them in the synthetic experiment.  If the prediction should be purely spatial, then it seems there should be no need to spatial splitting or buffering in cross-validation.  If the prediction should be purely based on features, again it should not matter as long as the spatial coordinates are not used.  It’s only when these should and do interact that this is an interesting challenge.  \n (Or if I am totally mis-understanding the issue .. again more clarification, I think, would be helpful).  \n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers support acceptance without much controversy. Having looked at the reviews and the updated version of the manuscript, I'm satisfied with the reply by the authors and support acceptance.
Decision: Accept (poster)