2. Choosing and diagnosing spatial folds
Introduction
This tutorial shows how to check whether a cross-validation design is
appropriate after the folds have been created. It covers
cv_summary, cv_similarity,
cv_distance, cv_spatial_autocor, and
cv_block_size.
The diagnostics answer different questions:
cv_summarygives a one-call overview of fold sizes, class prevalence, environmental novelty, nearest-neighbour distances, and automatically detected warnings.cv_similaritychecks whether testing folds contain novel environmental conditions relative to their training data.cv_distancechecks whether the test-to-train distances in the folds resemble the distances the model will face when predicting over the target area.cv_spatial_autocorandcv_block_sizehelp choose candidate spatial block sizes.
Please cite blockCV by: Valavi R, Elith J,
Lahoz-Monfort JJ, Guillera-Arroita G. blockCV: An R package for
generating spatially or environmentally separated folds for k-fold
cross-validation of species distribution models. Methods Ecol Evol.
2019; 10:225–232. doi:
10.1111/2041-210X.13107
Prepare data and folds
Load the example presence-absence records and environmental
covariates included with blockCV.
library(blockCV)
library(sf)
library(terra)
points <- read.csv(system.file("extdata/", "species.csv", package = "blockCV"))
pa_data <- sf::st_as_sf(points, coords = c("x", "y"), crs = 7845)
rasters <- terra::rast(
list.files(system.file("extdata/au/", package = "blockCV"), full.names = TRUE)
)Create several fold designs to compare. See Tutorial 1 for a fuller explanation of the fold-creation functions.
sb <- cv_spatial(x = pa_data,
column = "occ",
r = rasters,
k = 5,
size = 350000,
selection = "random",
iteration = 50,
progress = FALSE,
report = FALSE,
plot = FALSE)
set.seed(6)
bscv <- cv_cluster(x = pa_data,
column = "occ",
k = 5,
balance = TRUE,
k_multiplier = 3,
report = FALSE)
set.seed(6)
ecv <- cv_cluster(x = pa_data,
column = "occ",
r = rasters,
k = 5,
scale = TRUE)
bloo <- cv_buffer(x = pa_data,
column = "occ",
size = 350000,
progress = FALSE,
report = FALSE)
knn_blocks <- cv_knndm(x = pa_data,
column = "occ",
r = rasters,
k = 5,
clustering = "blocks",
keep_blocks = TRUE,
num_sample = 3000,
nk_len = 40,
seed = 6,
report = FALSE,
plot = FALSE)Fold-quality summary
The cv_summary function is the quickest way to inspect a
fold design. Without a prediction domain it reports the structural
summary from the cv object: number of folds, leave-one-out
status, train/test counts, and fold warnings.
## blockCV fold-quality summary
## Folds: 5
##
## Fold sizes and class prevalence (train/test counts):
## train_0 train_1 test_0 test_1
## 1 211 210 46 33
## 2 186 202 71 41
## 3 203 182 54 61
## 4 220 205 37 38
## 5 208 173 49 70
##
## Warnings:
## noneWhen r, pred_points, or
model_domain is supplied, cv_summary also adds
the distance diagnostic. When r is supplied, it also adds
the environmental novelty diagnostic. These are the same data-oriented
computations used by cv_distance and
cv_similarity.
summary_full <- cv_summary(cv = sb,
x = pa_data,
r = rasters,
num_sample = 3000,
seed = 6,
progress = FALSE)
summary_full$distances## fold n_test min q1 median q3 max pct_below_pred
## 1 1 79 8558 90160 134000 176000 258700 98.7
## 2 2 112 8558 76910 112900 159200 285600 98.2
## 3 3 115 12100 66010 104100 154200 311500 95.7
## 4 4 75 8558 56290 107200 175400 537200 82.7
## 5 5 119 8558 77020 108300 145700 367500 86.6## fold n_test pct_novel min median limiting_var
## 1 1 79 1.3 -11.33 18.53 bio_12
## 2 2 112 11.6 -7.02 4.64 bio_4
## 3 3 115 3.5 -5.15 12.99 bio_15
## 4 4 75 0.0 0.94 27.29 <NA>
## 5 5 119 2.5 -0.84 30.45 bio_12## fold type
## 1 1 high_leakage
## 2 2 high_leakage
## 3 3 high_leakage
## message
## 1 99% of test points nearer to train than the median prediction distance (optimistic/leaky fold)
## 2 98% of test points nearer to train than the median prediction distance (optimistic/leaky fold)
## 3 96% of test points nearer to train than the median prediction distance (optimistic/leaky fold)The warnings are returned as data, not as R warning conditions, so
they can be filtered or saved in an analysis pipeline. For k-fold
designs they can flag empty or tiny test folds, single-class test folds,
classes missing from training, severe imbalance, and high leakage. For
leave-one-out designs such as cv_buffer and
cv_nndm, structural fold-size warnings are skipped because
each fold holds a single test point by design.
Environmental similarity diagnostics
The cv_similarity function evaluates environmental
similarity between training and testing folds. It supports
method = "MESS", method = "L1", and
method = "L2". Negative values indicate environmental
novelty: at least one testing point is less similar to its fold’s
training data than the selected method’s reference.
The function returns a cv_similarity object. The plot is
stored in $plot, the per-fold table in
$extrapolation, and the overall novelty rate in
$overall. Set plot = FALSE when you want to
save or customise the plot object instead of drawing it immediately.
sim <- cv_similarity(cv = ecv,
x = pa_data,
r = rasters,
method = "MESS",
plot = FALSE,
progress = FALSE)
sim$plot
## fold n_test pct_novel min median limiting_var
## 1 1 79 17.7 -68.43 9.98 bio_12
## 2 2 51 82.4 -31.91 -12.88 bio_4
## 3 3 78 1.3 -4.64 25.12 bio_5
## 4 4 139 0.0 0.55 9.42 <NA>
## 5 5 153 60.1 -22.82 -1.83 bio_4## [1] 29.8Setting type = "map" shows where extrapolation occurs in
geographical space. Testing points are coloured by their similarity
value (here MESS), so novel points are easy to locate.
cv_similarity(cv = ecv,
x = pa_data,
r = rasters,
method = "MESS",
type = "map",
plot = TRUE,
progress = FALSE)
The distance-based L1 and L2 scores use a
random raster sample as their baseline, so set seed when
reproducibility matters. Use num_plots to show a subset of
folds.
cv_similarity(cv = ecv,
x = pa_data,
r = rasters,
method = "L2",
num_plots = 1:3,
num_sample = 3000,
seed = 6,
plot = TRUE,
progress = FALSE)
Nearest-neighbour distance diagnostics
The cv_distance function compares the nearest-neighbour
distance distribution of a fold design with the distance distribution
expected across the prediction domain. It returns a
cv_distance object with the plot in $plot,
per-fold distance summaries in $distances, and
Wasserstein-1 distances in $W. Lower Wasserstein-1 values
indicate a closer match to the prediction domain.
distance_theme <- ggplot2::theme(
legend.position = "bottom",
plot.subtitle = ggplot2::element_text(size = 7)
)
dist_sb <- cv_distance(cv = sb,
x = pa_data,
r = rasters,
num_sample = 3000,
num_random = 5,
seed = 6,
plot = FALSE)
dist_bscv <- cv_distance(cv = bscv,
x = pa_data,
r = rasters,
num_sample = 3000,
num_random = 5,
seed = 6,
plot = FALSE)
dist_bloo <- cv_distance(cv = bloo,
x = pa_data,
r = rasters,
add_random = FALSE,
num_sample = 3000,
seed = 6,
plot = FALSE)
dist_knn <- cv_distance(cv = knn_blocks,
x = pa_data,
r = rasters,
num_sample = 3000,
num_random = 5,
seed = 6,
plot = FALSE)
cowplot::plot_grid(
dist_sb$plot +
ggplot2::labs(title = "Spatial blocks") +
distance_theme,
dist_bscv$plot +
ggplot2::labs(title = "Spatial clustering (balanced)") +
distance_theme,
dist_bloo$plot +
ggplot2::labs(title = "Buffering LOO") +
distance_theme,
dist_knn$plot +
ggplot2::labs(title = "kNNDM blocks") +
distance_theme,
ncol = 2
)
The W vector summarises the area between each curve and
the prediction curve. The distances table shows the spread
of test-to-nearest-train distances within each fold. The
pct_below_pred column reports the percentage of test points
that are closer to training data than the median prediction distance;
high values indicate folds that may be easier than the prediction
task.
## CV LOO Random
## 179895.2 263821.6 259662.7## fold n_test min q1 median q3 max pct_below_pred
## 1 1 79 8558 90160 134000 176000 258700 98.7
## 2 2 112 8558 76910 112900 159200 285600 98.2
## 3 3 115 12100 66010 104100 154200 311500 95.7
## 4 4 75 8558 56290 107200 175400 537200 82.7
## 5 5 119 8558 77020 108300 145700 367500 86.6Distances are geographical by default. Set
space = "feature" to compute distances in the covariate
space of r, optionally scaled with
scale = TRUE.
dist_feature <- cv_distance(cv = sb,
x = pa_data,
r = rasters,
space = "feature",
num_sample = 3000,
seed = 6,
plot = FALSE)
dist_feature$plot
Choosing a candidate block size
To support a first choice of block size before model fitting,
cv_spatial_autocor estimates the effective range of spatial
autocorrelation in the response or continuous predictor rasters. This
tool is a guide, not an absolute rule.
Roberts et al. (2017) recommend choosing spatial blocks with
reference to the autocorrelation range of model residuals. The response
and predictor variograms shown here are useful exploratory guides before
a model has been fitted, but they are not residual autocorrelation
ranges and can under- or over-estimate the block size needed for
residual-based guidance. To estimate residual autocorrelation with
cv_spatial_autocor, fit the model first, add the residuals
as a column in the sample points object, and pass that column name to
column.
When only r is supplied, the function fits variograms to
the continuous raster covariates and reports the median range.
sac_raster <- cv_spatial_autocor(r = rasters,
num_sample = 3000,
progress = FALSE,
plot = FALSE)
plot(sac_raster)
## [1] 1651696## Summary statistics of spatial autocorrelation ranges of all input layers:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 972128 1300689 1651696 6631027 6982034 22248586
## layers range
## 3 bio_4 972128.3
## 2 bio_15 1410209.4
## 4 bio_5 1893182.7
## 1 bio_12 22248586.3Alternatively, supply the sample points and a response, count, continuous, or residual column.
cv_block_size provides an interactive Shiny app for
manually exploring spatial block sizes. The app cannot be shown inside
the vignette because it requires an interactive R session.
References:
Elith, J., Kearney, M., & Phillips, S. (2010) The art of modelling range-shifting species. Methods in Ecology and Evolution, 1(4), 330-342.
Hiemstra, P. H., Pebesma, E. J., Twenhöfel, C. J., & Heuvelink, G. B. (2009) Real-time automatic interpolation of ambient gamma dose rates from the Dutch radioactivity monitoring network. Computers & Geosciences, 35(8), 1711-1721.
Linnenbrink, J., Milà, C., Ludwig, M., & Meyer, H. (2024) kNNDM CV: k-fold nearest neighbour distance matching cross-validation for map accuracy estimation. Geoscientific Model Development, 17(15), 5897-5912. doi: 10.5194/gmd-17-5897-2024
Milà, C., Mateu, J., Pebesma, E., & Meyer, H. (2022) Nearest neighbour distance matching leave-one-out cross-validation for map validation. Methods in Ecology and Evolution, 13(6), 1304-1316.
O’Sullivan, D., & Unwin, D. J. (2010) Geographic Information Analysis (2nd ed.). John Wiley & Sons.
Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J. J., Schroeder, B., Thuiller, W., Warton, D. I., Wintle, B. A., Hartig, F., & Dormann, C. F. (2017) Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography, 40, 913-929.
Valavi R, Elith J, Lahoz-Monfort JJ, Guillera-Arroita G. (2019) blockCV: An R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods Ecol Evol. 10:225-232. doi: 10.1111/2041-210X.13107