--- title: "2. Choosing and diagnosing spatial folds" author: "Roozbeh Valavi" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: fig_caption: yes vignette: > %\VignetteIndexEntry{2. Choosing and diagnosing spatial folds} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## 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_summary` gives a one-call overview of fold sizes, class prevalence, environmental novelty, nearest-neighbour distances, and automatically detected warnings. - `cv_similarity` checks whether testing folds contain novel environmental conditions relative to their training data. - `cv_distance` checks whether the test-to-train distances in the folds resemble the distances the model will face when predicting over the target area. - `cv_spatial_autocor` and `cv_block_size` help 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](https://doi.org/10.1111/2041-210X.13107) ```{r echo=FALSE} options(scipen = 10) ``` ## Prepare data and folds Load the example presence-absence records and environmental covariates included with `blockCV`. ```{r message=FALSE, warning=FALSE} 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](tutorial_1.html) for a fuller explanation of the fold-creation functions. ```{r warning=FALSE, message=FALSE} 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. ```{r} summary_basic <- cv_summary(sb) summary_basic ``` When `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`. ```{r} summary_full <- cv_summary(cv = sb, x = pa_data, r = rasters, num_sample = 3000, seed = 6, progress = FALSE) summary_full$distances summary_full$novelty summary_full$warnings ``` 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. ```{r fig.height=4, fig.width=6} sim <- cv_similarity(cv = ecv, x = pa_data, r = rasters, method = "MESS", plot = FALSE, progress = FALSE) sim$plot sim$extrapolation sim$overall ``` Setting `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. ```{r fig.height=4.5, fig.width=6} 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. ```{r fig.height=4, fig.width=6} 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. ```{r warning=FALSE, message=FALSE, fig.height=7, fig.width=8} 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. ```{r} dist_sb$W dist_sb$distances ``` Distances are geographical by default. Set `space = "feature"` to compute distances in the covariate space of `r`, optionally scaled with `scale = TRUE`. ```{r fig.height=4, fig.width=6} 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. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7.2} sac_raster <- cv_spatial_autocor(r = rasters, num_sample = 3000, progress = FALSE, plot = FALSE) plot(sac_raster) ``` ```{r} sac_raster$range summary(sac_raster) ``` Alternatively, supply the sample points and a response, count, continuous, or residual column. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7.2} sac_points <- cv_spatial_autocor(x = pa_data, column = "occ", plot = FALSE) plot(sac_points) ``` `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. ```{r eval=FALSE} cv_block_size(x = pa_data, column = "occ", r = rasters, min_size = 2e5, max_size = 9e5) ``` ## 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](https://doi.org/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](https://doi.org/10.1111/2041-210X.13107)