--- title: "1. blockCV introduction: how to create block cross-validation folds" author: "Roozbeh Valavi" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: fig_caption: yes vignette: > %\VignetteIndexEntry{1. blockCV introduction: how to create block cross-validation folds} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## Introduction The package `blockCV` offers a range of functions for generating train and test folds for **k-fold** and **leave-one-out (LOO)** cross-validation (CV). It allows for separation of data spatially, environmentally, by nearest-neighbour distance matching, or by user-defined groups. The `blockCV` package is suitable for the evaluation of a variety of spatial modelling applications, including classification of remote sensing imagery, soil mapping, and species distribution modelling (SDM). It also provides support for different SDM scenarios, including presence-absence and presence-background species data, rare and common species, and raster data for predictor variables. This tutorial demonstrates `cv_spatial`, `cv_cluster`, `cv_group`, `cv_buffer`, `cv_nndm`, and `cv_knndm`, plus tools for plotting folds. See Tutorial 2 (open with `vignette("tutorial_2", package = "blockCV")`) for fold-quality diagnostics and block-size selection. 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) ## Installation The `blockCV` is available in CRAN and the latest update can also be downloaded from GitHub. It is recommended to install the dependencies of the package. To install the package use: ```{r, eval=FALSE} # install stable version from CRAN install.packages("blockCV", dependencies = TRUE) # install latest update from GitHub remotes::install_github("rvalavi/blockCV", dependencies = TRUE) ``` ```{r message=TRUE, warning=TRUE} # loading the package library(blockCV) ``` ## Package data The package contains the raw format of the following data: - Raster covariates of Australia (`.tif`) - Simulated species data (`.csv`) These data are used to illustrate how the package is used. The raster data include several bioclimatic variables for Australia. The species data include presence-absence records (binary) of a simulated species. To load the package raster data use: ```{r, fig.height=5, fig.width=7.2, warning=FALSE, message=FALSE} library(sf) # working with spatial vector data library(terra) # working with spatial raster data library(tmap) # plotting maps # load raster data rasters <- terra::rast( list.files(system.file("extdata/au/", package = "blockCV"), full.names = TRUE) ) ``` The presence-absence species data include `243` presence points and `257` absence points. ```{r, fig.height=4.5, fig.width=7.1} # load species presence-absence data and convert to sf points <- read.csv(system.file("extdata/", "species.csv", package = "blockCV")) head(points) ``` The appropriate format of species data for the `blockCV` package is simple features (from the `sf` package). The data is provide in [GDA2020 / GA LCC](https://epsg.io/7845) coordinate reference system with `"EPSG:7845"` as defined by `crs = 7845`. We convert the `data.frame` to `sf` as follows: ```{r, fig.height=4.5, fig.width=7.1} pa_data <- sf::st_as_sf(points, coords = c("x", "y"), crs = 7845) ``` Let's plot the species data using [`tmap`](https://cran.r-project.org/package=tmap) package: ```{r, fig.height=4.5, fig.width=7.1} tm_shape(rasters[[1]]) + tm_raster( col.scale = tm_scale_continuous(values = gray.colors(10)), col.legend = tm_legend_hide() ) + tm_shape(pa_data) + tm_dots( fill = "occ", fill.scale = tm_scale_categorical(), size = 0.5, fill_alpha = 0.5 ) ``` ## Block cross-validation strategies The `blockCV` stores training and testing folds in a common list format: `folds_list`, a list of the indices of observations in each fold. For k-fold methods that assign each observation to one fold, the folds are also stored as `folds_ids`, a vector of fold numbers for each observation, and often as `biomod_table`, a matrix suitable for the `biomod2` package. Use the format that suits your modelling workflow. ### Spatial blocks The function `cv_spatial` creates spatial blocks/polygons then assigns blocks to the training and testing folds with *random*, *checkerboard pattern* or a *systematic* way (with the selection argument). When `selection = "random"`, the function tries to find evenly distributed records in training and testing folds. Spatial blocks can be defined either by `size` or number of rows and columns. Consistent with other functions, the distance (`size`) should be in **metres**, regardless of the unit of the reference system of the input data. When the input map has *geographic coordinate system* (i.e. decimal degrees), the block size is calculated based on dividing `size` by 111325 (the standard distance of a degree in metres, on the Equator) to change metre to degree. In reality, this value varies by a factor of the cosine of the latitude. So, an alternative sensible value could be `cos(mean(sf::st_bbox(x)[c(2,4)]) * pi/180) * 111325`. The `offset` argument can be used to shift the spatial position of the blocks in horizontal and vertical axes, respectively. This only works when the block have been built based on `size`, and the `extend` option allows user to enlarge the blocks ensuring all points fall inside the blocks (most effectve when `rows_cols` is used). The blocks argument allows users to define an external spatial polygon as blocking layer. Here are some spatial block settings: ```{r, results='hide', fig.keep='all', warning=FALSE, message=FALSE, fig.height=5, fig.width=7} sb1 <- cv_spatial(x = pa_data, column = "occ", # the response column (binary or multi-class) k = 5, # number of folds size = 350000, # size of the blocks in metres selection = "random", # random blocks-to-fold iteration = 50, # find evenly dispersed folds biomod2 = TRUE) # also create folds for biomod2 ``` The output object is an R S3 object and you can get its elements by a `$`. Explore `sb1$folds_ids`, `sb1$folds_list`, and `sb1$biomod_table` for the three types of generated folds from the `cv_spatial` object `sb1`. Use the one suitable for you modelling practice to evaluate your models. See the explanation of all other outputs/elements of the function in the help file of the function. The default spatial blocks are hexagons, and you can make this explicit with `hexagon = TRUE`. You can optionally add a raster layer (using `r` argument) to be used for creating blocks and in the background of the plot (a raster can also be added later only for visualising blocks using `cv_plot`). ```{r, warning=FALSE, message=FALSE, fig.height=5, fig.width=7} sb2 <- cv_spatial(x = pa_data, column = "occ", r = rasters, # optionally add a raster layer k = 5, size = 350000, hexagon = TRUE, # use hexagonal blocks selection = "random", progress = FALSE, # turn off progress bar for vignette iteration = 50, biomod2 = TRUE, report = TRUE, plot = TRUE) ``` The assignment of folds to each block can also be done in a systematic manner using `selection = "systematic"`, or a checkerboard pattern using `selection = "checkerboard"`. Square blocks are created with `hexagon = FALSE`, as shown below; checkerboard selection requires square blocks. The blocks can also be created by number of rows and columns when no `size` is supplied by e.g. `rows_cols = c(12, 10)`. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7} # systematic fold assignment # and also use row/column for creating blocks instead of size sb3 <- cv_spatial(x = pa_data, column = "occ", rows_cols = c(12, 10), hexagon = FALSE, selection = "systematic", report = TRUE, plot = TRUE) ``` The function's output report reveals that setting the selection to 'random' results in a more even distribution of presence/absence instances between the train and test folds compared to 'systematic'. This is because the random assignment process is repeated multiple times, controlled by the `iteration` parameter, to ensure that the folds are evenly distributed. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7} # checkerboard block to CV fold assignment sb4 <- cv_spatial(x = pa_data, column = "occ", size = 350000, hexagon = FALSE, selection = "checkerboard", report = TRUE, plot = TRUE) ``` Let's visualise the checkerboard blocks with `tmap`: ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7} tm_shape(sb4$blocks) + tm_fill( fill = "folds", fill.scale = tm_scale_categorical() ) ``` ### Spatial and environmental clustering The function `cv_cluster` uses *clustering* methods to specify sets of similar environmental conditions based on the input covariates. Species data corresponding to any of these groups or clusters are assigned to a fold. Alternatively, the clusters can be based on spatial coordinates of sample points (the `x` argument). Using spatial coordinate values for clustering: ```{r} # spatial clustering set.seed(6) scv <- cv_cluster(x = pa_data, column = "occ", # optional: counting number of train/test records k = 5, report = TRUE) ``` The default spatial clustering keeps each fold spatially compact, but the number of records or classes in each fold can be uneven. To prioritise more even folds, set `balance = TRUE`. In this case, `cv_cluster()` first creates `k * k_multiplier` smaller clusters and then assigns those clusters to `k` folds. ```{r} # balanced spatial clustering set.seed(6) bscv <- cv_cluster(x = pa_data, column = "occ", k = 5, balance = TRUE, k_multiplier = 3, report = TRUE) ``` ```{r warning=FALSE, message=FALSE, fig.height=4, fig.width=8} scv_plot <- cv_plot(scv, pa_data, combine_folds = TRUE) + ggplot2::labs(title = "Default clustering") bscv_plot <- cv_plot(bscv, pa_data, combine_folds = TRUE) + ggplot2::labs(title = "Balanced clustering") cowplot::plot_grid(scv_plot, bscv_plot, nrow = 1) ``` Increasing `k_multiplier` usually improves the balance of records or classes across folds because there are more small clusters to assign, but it reduces the spatial separation and compactness of the final folds. For species **presence-background** data (presences and background or pseudo-absence points), supply the response `column` and set `presence_bg = TRUE`. The fold balancing is then driven by the presence points (`1`s) only, so that the usually far more numerous background points do not dominate the split. This applies to the balancing in `cv_spatial`, `cv_cluster`, and `cv_knndm`. The example below uses the presence-background data supplied with the package in `species_pb.csv`. ```{r} # load species presence-background data and convert to sf points_pb <- read.csv(system.file("extdata/", "species_pb.csv", package = "blockCV")) pb_data <- sf::st_as_sf(points_pb, coords = c("x", "y"), crs = 7845) # balance folds on presences only (presence-background data) set.seed(6) pbcv <- cv_cluster(x = pb_data, column = "occ", k = 5, balance = TRUE, presence_bg = TRUE, # balance on presences (1s) only report = TRUE) ``` The clustering can be done in environmental space by supplying `r`. Notice, this could be an extreme case of cross-validation as the testing folds could possibly fall in novel environmental conditions compared with the training points. Tutorial 2 shows how to check this with `cv_similarity`. Note that the input raster layer should cover all the species points, otherwise an error will rise. The records with no raster value should be deleted prior to the analysis or a different raster be used. ```{r warning=FALSE, message=FALSE} # environmental clustering set.seed(6) ecv <- cv_cluster(x = pa_data, column = "occ", r = rasters, k = 5, scale = TRUE) ``` When `r` is supplied, all the input rasters are first centred and scaled to avoid one raster variable dominate the clusters using `scale = TRUE` option. By default, the clustering will be done based only on the values of the predictors at the sample points. In this case, and the number of the folds will be the same as `k`. If `raster_cluster = TRUE`, the clustering is done in the raster space. In this approach, the clusters will be consistent throughout the region and across species (in the same region). However, this may result in cluster(s) that cover none of the species records especially when species data is not dispersed throughout the region (or environmental ranges) or the number of clusters (`k` or folds) is high. Purely environmental clusters can be scattered across the study area, since points with similar covariate values may be far apart in space. The `spatial_weight` argument (a number in `[0, 1]`) adds a soft spatial-compactness pressure by blending the coordinates into the environmental clustering: `spatial_weight = 0` (the default) is pure environmental clustering, larger values pull the clusters to be more geographically compact, and `spatial_weight = 1` clusters on the coordinates alone. A moderate value such as `0.4` keeps the folds environmentally coherent while making them more spatially separated. ```{r warning=FALSE, message=FALSE} # spatially constrained environmental clustering set.seed(6) secv <- cv_cluster(x = pa_data, column = "occ", r = rasters, k = 5, scale = TRUE, spatial_weight = 0.4) # blend geography into the environmental clusters ``` ```{r warning=FALSE, message=FALSE, fig.height=4, fig.width=8} ecv_plot <- cv_plot(ecv, pa_data, combine_folds = TRUE) + ggplot2::labs(title = "Environmental clustering") secv_plot <- cv_plot(secv, pa_data, combine_folds = TRUE) + ggplot2::labs(title = "Spatially constrained (spatial_weight = 0.4)") cowplot::plot_grid(ecv_plot, secv_plot, nrow = 1) ``` ### Leave-group-out cross-validation The function `cv_group` creates folds from an existing grouping column, such as a site, plot, sampling campaign, or individual identifier. Records that share a group are always kept together, so the same group never appears in both the training and testing set of a fold. The package example data do not include a site identifier, so the next chunk creates an illustrative grouping column. In your own data, use the real grouping variable that defines the dependence structure you want to hold out. ```{r} set.seed(6) coords <- sf::st_coordinates(pa_data) pa_data$site <- paste0("site_", stats::kmeans(coords, centers = 10)$cluster) ``` With `k = NULL` (the default), `cv_group()` performs leave-group-out cross-validation: each group is left out once. ```{r} site_lgo <- cv_group(x = pa_data, group_col = "site", column = "occ", report = TRUE) ``` If there are many groups, you can merge them into fewer folds by setting `k`. With `balance = TRUE`, the function searches for a group-to-fold assignment with more even record or class counts, while still keeping whole groups together. ```{r} site_cv <- cv_group(x = pa_data, group_col = "site", column = "occ", k = 5, balance = TRUE, iteration = 50, seed = 6, report = TRUE) ``` ```{r warning=FALSE, message=FALSE, fig.height=4.5, fig.width=6} cv_plot(site_cv, pa_data, combine_folds = TRUE) + ggplot2::labs(title = "Grouped folds") ``` ### Buffering LOO (also known as Spatial LOO) The function `cv_buffer` generates spatially separated training and testing folds by considering buffers of specified distance around each observation point. This approach is a form of leave-one-out (LOO) cross-validation. Each fold is generated by excluding nearby observations around each testing point within the specified distance (ideally the range of spatial autocorrelation). In this method the test set never directly abuts a training set. Using buffering to create CV folds: ```{r results='hide', fig.keep='all'} bloo <- cv_buffer(x = pa_data, column = "occ", size = 350000) ``` When using species **presence-background** data (or presence and pseudo-absence), you need to supply the `column` and set `presence_bg = TRUE`. In this case, only presence points (1s) are considered as target points. For more information read the details section in the help of the function (i.e. `help(cv_buffer)`). For species **presence-absence** data and any other types of data (such as **continuous**, **counts**, and **multi-class** targets) keep `presence_bg = FALSE` (default). In this case, all sample points other than the target point within the buffer are excluded, and the training set comprises all points outside the buffer. ## Nearest Neighbour Distance Matching (NNDM) LOO The `cv_nndm` is a fast implementation of the Nearest Neighbour Distance Matching (NNDM) algorithm (Milà et al., 2022) in C++. Similar to `cv_buffer`, this is a variation of leave-one-out (LOO) cross-validation. It tries to match the nearest neighbour distance distribution function between the test and training data to the nearest neighbour distance distribution function between the target prediction and training points (Milà et al., 2022). ```{r fig.height=4, fig.width=7} nncv <- cv_nndm(x = pa_data, column = "occ", r = rasters, size = 350000, num_sample = 5000, sampling = "random", min_train = 0.1, plot = TRUE) ``` ## k-fold Nearest Neighbour Distance Matching (kNNDM) The `cv_knndm` function implements the k-fold version of NNDM (Linnenbrink et al., 2024). Unlike `cv_nndm`, it returns a smaller set of `k` folds, making it cheaper to fit models while still matching the nearest neighbour distance distribution between testing and training data to the prediction-to-training distribution. The original kNNDM grouping options are `clustering = "hierarchical"` and `clustering = "kmeans"`. The example below uses hierarchical clustering; use `clustering = "kmeans"` for the k-means alternative. For a quick vignette run, the examples use `num_sample = 5000` and `nk_len = 50`; for applied analyses, consider using the defaults or larger values, especially for large prediction areas or complex sampling patterns. ```{r warning=FALSE, message=FALSE, fig.height=4, fig.width=7} knn_hier <- cv_knndm(x = pa_data, column = "occ", r = rasters, k = 5, clustering = "hierarchical", num_sample = 5000, nk_len = 50, seed = 6, report = TRUE, plot = TRUE) ``` The default `clustering = "blocks"` option is specific to `blockCV`. It first groups sample points with spatial blocks (hexagonal by default, or square with `hexagon = FALSE`) and then allocates those groups to kNNDM folds. ```{r warning=FALSE, message=FALSE, fig.height=4, fig.width=7} knn_blocks <- cv_knndm(x = pa_data, column = "occ", r = rasters, k = 5, clustering = "blocks", keep_blocks = TRUE, num_sample = 5000, nk_len = 50, seed = 6, report = TRUE, plot = TRUE) ``` ```{r warning=FALSE, message=FALSE, fig.height=4.5, fig.width=7} cv_plot(knn_blocks, pa_data, combine_folds = TRUE) ``` ## Visualising the folds You can visualise the generate folds for all block cross-validation strategies. You can optionally add a raster layer as background map using `r` option. When `r` is supplied the plots might be slightly slower. The output of `cv_plot()` is a `ggplot` object, so it can be customised with standard `ggplot2` layers. For example, you can add a title, change the theme, or modify colours and scales. Let's plot spatial clustering folds created in previous section (using `cv_cluster`): ```{r warning=FALSE, message=FALSE, fig.height=6, fig.width=8} cv_plot(cv = scv, x = pa_data) ``` Optionally, all folds can be plotted in a single panel, with colours indicating different folds. This option is available for k-fold cross-validation methods: `cv_spatial`, `cv_cluster`, `cv_group`, and `cv_knndm`. ```{r warning=FALSE, message=FALSE, fig.height=4.5, fig.width=6} cv_plot(cv = scv, x = pa_data, combine_folds = TRUE) + ggplot2::labs( title = "Spatial clustering folds", x = "Longitude", y = "Latitude" ) ``` When the fold object was built with `presence_bg = TRUE`, `cv_plot()` draws the background points (`0`s) faded so the presences (`1`s) stand out. The background opacity is set with `bg_alpha`; it is capped at `points_alpha` so the background is never more prominent than the presences (set `bg_alpha = points_alpha` to disable the fading). ```{r warning=FALSE, message=FALSE, fig.height=4.5, fig.width=6} cv_plot(cv = pbcv, x = pb_data, combine_folds = TRUE, bg_alpha = 0.05) + ggplot2::labs(title = "Presence-background folds (background faded)") ``` When `cv_buffer` is used for plotting, only first 10 folds are shown. You can choose any set of CV folds for plotting. If `remove_na = FALSE` (default is `TRUE`), the `NA` in the legend shows the excluded points. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=8} cv_plot(cv = bloo, x = pa_data, num_plots = c(1, 50, 100)) # only show folds 1, 50 and 100 ``` If you do not supply `x` when plotting a `cv_spatial` object, only the spatial blocks are plotted. ```{r warning=FALSE, message=FALSE, fig.height=5, fig.width=7} cv_plot(cv = sb1, r = rasters, raster_colors = terrain.colors(10, alpha = 0.5), label_size = 4) ``` ## Next steps After creating folds, use Tutorial 2 (open with `vignette("tutorial_2", package = "blockCV")`) to check fold size balance, environmental novelty, nearest-neighbour distance matching, and candidate spatial block sizes. ## References: - 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. and Meyer H. (2022) Nearest Neighbour Distance Matching Leave-One-Out Cross-Validation for map validation. Methods in Ecology and Evolution. - 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)