Nothing
R/knndm.R
In CAST: 'caret' Applications for Spatial-Temporal Models
Defines functions
distclust_distmat
knndm_feature
knndm_geo
check_knndm_feature
check_knndm_geo
knndm
Documented in
knndm
#' K-fold Nearest Neighbour Distance Matching
#' @description
#' This function implements the kNNDM algorithm for prediction-domain adaptive resampling
#' and returns the necessary indices to perform train-test splits or k-fold NNDM CV.
#'
#' @author Carles Milà and Jan Linnenbrink
#' @param tpoints sf or sfc point object, or data.frame if dist_space = "feature". Contains the training points samples.
#' @param modeldomain sf polygon object or SpatRaster defining the prediction area. Optional; alternative to predpoints (see Details).
#' @param predpoints sf or sfc point object, or data.frame if dist_space = "feature". Contains the target prediction points. Optional; alternative to modeldomain (see Details).
#' @param dist_space character. Either "geographical" or "feature".
#' @param k integer. Number of folds desired for CV. Defaults to 10.
#' @param maxp numeric. Maximum fold size allowed, defaults to 0.5, i.e. a single fold can hold a maximum of half of the training points.
#' @param clustering character. Possible values include "hierarchical" and "kmeans". See details.
#' @param linkf character. Only relevant if clustering = "hierarchical". Link function for agglomerative hierarchical clustering.
#' Defaults to "ward.D2". Check `stats::hclust` for other options.
#' @param samplesize numeric. How many points in the modeldomain should be sampled as prediction points?
#' Only required if modeldomain is used instead of predpoints.
#' @param sampling character. How to draw prediction points from the modeldomain? See `sf::st_sample`.
#' Only required if modeldomain is used instead of predpoints.
#' @param dist_fun character. Currently covers `euclidean` (default), `gower`, `mahalanobis` and `great_circle`.
#' `gower` and `mahalanobis` only work with `dist_space`="feature", while `great_circle` only works with `dist_space`="geographical".
#' `mahalanobis` takes into account correlation between predictor values. While `euclidean` and `mahalanobis` only work with numerical variables,
#' `gower` also works with mixed data including numerical and categorical variables.
#' For the geographical space, `great_circle` covers lon/lat coordinates, whereas `euclidean` only works with projected coordinates.
#' @param algorithm see \code{\link[FNN]{knnx.dist}} and \code{\link[FNN]{knnx.index}}
#' @param scale_vars boolean. Should variables be scaled? Only for `dist_space`="feature".
#' Calculating Gower distances already includes scaling, and manually rescale the data is redundant.
#' For other distances (Mahalanobis, Euclidean), scaling the data is important. Thus, TRUE by default.
#' @param test_prop numeric. The proportion of test data. NULL by default (i.e., no train/test split).
#' @param test_tolerance numeric. The allowed deviance from `test_prop`. The higher the tolerance,
#' the larger the possibility to obtain train/test splits that yield good approximations of the prediction situation.
#' @param nk_len integer. The number of fold configurations to test. By default 100.
#' Larger numbers increase computational times, but also might lead to better W statistics.
#' Useful for train/test splits, where a large number of configurations is discarded.
#' @param space deprecated. Use `dist_space` instead.
#' @param useMD deprecated. Use `dist_fun` instead.
#' @return An object of class \emph{knndm} consisting of a list of eight elements:
#' indx_train, indx_test (indices of the observations to use as
#' training/test data in each kNNDM CV iteration), Gij (distances for
#' G function construction between prediction and target points), Gj
#' (distances for G function construction during LOO CV), Gjstar (distances
#' for modified G function during kNNDM CV), clusters (list of cluster IDs),
#' W (Wasserstein statistic), and dist_space (stated by the user in the function call).
#'
#' @details
#' knndm is an implementation of prediction-domain adaptive validation.
#' It is a k-fold version of NNDM LOO CV which makes it more suitable for medium and large datasets.
#' It can be used for cross-validation and train / test splits (the latter is experimental).
#' Briefly, the algorithm tries to find a configuration such that the integral of the absolute differences (Wasserstein W statistic)
#' between the empirical nearest neighbour distance distribution function between the test and training data (Gj*),
#' and the empirical nearest neighbour distance distribution function between the prediction and training points (Gij),
#' is minimised. It does so by performing clustering of the training points' coordinates for different numbers of
#' clusters that range from k to N (number of observations), merging them into k final folds,
#' and selecting the configuration with the lowest W.
#'
#' When using `knndm` to split the data into training and test sets (experimental), the proportion of points belonging to the test set (`test_prop`) replaces the number of folds `k`.
#' Based on the `test_prop` , `minp` and `maxp` are calculated as the `test_prop` +/- `test_tolerance`.
#' Compared to k-fold CV, using knndm for train/test splits is less flexible and often results in larger NNDs between test and train locations
#' than between prediction and train locations. Hence, it is essential to plot the results of `knndm` and check how well the split can resemble the prediction situation.
#' Modifying the `test_prop` parameter, as well as increasing `test_prop` allow more flexible matching and can potentially improve the match.
#'
#' Using a projected CRS in `knndm` has large computational advantages since fast nearest neighbour search can be
#' done via the `FNN` package, while working with geographic coordinates requires computing the full
#' spherical distance matrices. As a clustering algorithm, `kmeans` can only be used for
#' projected CRS while `hierarchical` can work with both projected and geographical coordinates, though it requires
#' calculating the full distance matrix of the training points even for a projected CRS.
#'
#' In order to select between clustering algorithms and number of folds `k`, different `knndm` configurations can be run
#' and compared, being the one with a lower W statistic the one that offers a better match. W statistics between `knndm`
#' runs are comparable as long as `tpoints` and `predpoints` or `modeldomain` stay the same.
#'
#' Map validation using `knndm` should be used using `CAST::global_validation`, i.e. by stacking all out-of-sample
#' predictions and evaluating them all at once. The reasons behind this are 1) The resulting folds can be
#' unbalanced and 2) nearest neighbour functions are constructed and matched using all CV folds simultaneously.
#'
#' If training data points are very clustered with respect to the prediction area and the presented `knndm`
#' configuration still show signs of Gj* > Gij, there are several things that can be tried. First, increase
#' the `maxp` parameter; this may help to control for strong clustering (at the cost of having unbalanced folds).
#' Secondly, decrease the number of final folds `k`, which may help to have larger clusters.
#'
#' The `modeldomain` is either a sf polygon that defines the prediction area, or alternatively a SpatRaster out of which a polygon,
#' transformed into the CRS of the training points, is defined as the outline of all non-NA cells.
#' Then, the function takes a regular point sample (amount defined by `samplesize`) from the spatial extent.
#' As an alternative use `predpoints` instead of `modeldomain`, if you have already defined the prediction locations (e.g. raster pixel centroids).
#' When using either `modeldomain` or `predpoints`, we advise to plot the study area polygon and the training/prediction points as a previous step to ensure they are aligned.
#'
#' `knndm` can also be performed in the feature space by setting `dist_space` to "feature".
#' Euclidean distances, Gower distance or Mahalanobis distances can be used for distance calculation, but only Euclidean are tested.
#' In this case, nearest neighbour distances are calculated in n-dimensional feature space rather than in geographical space.
#' `tpoints` and `predpoints` can be data frames or sf objects containing the values of the features. Note that the names of `tpoints` and `predpoints` must be the same.
#' `predpoints` can also be missing, if `modeldomain` is of class SpatRaster. In this case, the values of of the SpatRaster will be extracted to the `predpoints`.
#' In the case of any categorical features, Gower distances will be used to calculate the Nearest Neighbour distances [Experimental]. If categorical
#' features are present, and `clustering` = "kmeans", K-Prototype clustering will be performed instead.
#'
#' @note
#' For spatial visualization of fold affiliation see examples.
#' @references
#' \itemize{
#' \item Linnenbrink, J., Milà, C., Ludwig, M., and Meyer, H. (2024): kNNDM: k-fold Nearest Neighbour Distance Matching Cross-Validation for map accuracy estimation. Geosci. Model Dev., 17, 5897–5912. https://doi.org/10.5194/gmd-17-5897-2024.
#' \item 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, 1304– 1316. https://doi.org/10.1111/2041-210X.13851.
#' }
#' @seealso \code{\link{geodist}}, \code{\link{nndm}}
#'
#' @export
#' @examples
#' ########################################################################
#' # Example 1: Simulated data - Randomly-distributed training points
#' ########################################################################
#'
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 random training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- sf::st_sample(simarea, 1000, type = "random")
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Run kNNDM for the whole domain, here the prediction points are known.
#' knndm_folds <- knndm(train_points, predpoints = pred_points, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' plot(knndm_folds, type = "simple") # For more accessible legend labels
#' plot(knndm_folds, type = "simple", stat = "density") # To visualize densities rather than ECDFs
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'
#' ########################################################################
#' # Example 2: Simulated data - Clustered training points
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 clustered training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- clustered_sample(simarea, 1000, 50, 5)
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Run kNNDM for the whole domain, here the prediction points are known.
#' knndm_folds <- knndm(train_points, predpoints = pred_points, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' plot(knndm_folds, type = "simple") # For more accessible legend labels
#' plot(knndm_folds, type = "simple", stat = "density") # To visualize densities rather than ECDFs
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'}
#' ########################################################################
#' # Example 3: Real- world example; using a modeldomain instead of previously
#' # sampled prediction locations
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(terra)
#' library(ggplot2)
#'
#' ### prepare sample data:
#' data(cookfarm)
#' dat <- aggregate(cookfarm[,c("DEM","TWI", "NDRE.M", "Easting", "Northing","VW")],
#' by=list(as.character(cookfarm$SOURCEID)),mean)
#' pts <- dat[,-1]
#' pts <- st_as_sf(pts,coords=c("Easting","Northing"))
#' st_crs(pts) <- 26911
#' studyArea <- rast(system.file("extdata","predictors_2012-03-25.tif",package="CAST"))
#' pts <- st_transform(pts, crs = st_crs(studyArea))
#' terra::plot(studyArea[["DEM"]])
#' terra::plot(vect(pts), add = T)
#'
#' knndm_folds <- knndm(pts, modeldomain=studyArea, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = pts, aes(col = folds))
#'
#' #use for cross-validation:
#' library(caret)
#' ctrl <- trainControl(method="cv",
#' index=knndm_folds$indx_train,
#' savePredictions='final')
#' model_knndm <- train(dat[,c("DEM","TWI", "NDRE.M")],
#' dat$VW,
#' method="rf",
#' trControl = ctrl)
#' global_validation(model_knndm)
#'}
#' ########################################################################
#' # Example 4: Simulated data - Train/test split with clustered training points
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 clustered training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- clustered_sample(simarea, 1000, 50, 5)
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Use kNNDM to split the data into 30% +- 10% test and 70% train
#' knndm_folds <- knndm(train_points, predpoints = pred_points, test_prop = 0.3, test_tolerance = 0.1)
#' # How many samples have been used for testing:
#' table(knndm_folds$clusters)
#' plot(knndm_folds)
#' # The train/test split could not represent the prediction situation well
#' # Increase tolerance to increase number of configurations tried, and thus to find a suitable split
#' knndm_folds <- knndm(train_points, predpoints = pred_points, test_prop = 0.3, test_tolerance = 0.2)
#' plot(knndm_folds)
#' table(knndm_folds$clusters)
#' # This resulted in better match of the prediction situation, but a 50/50 split
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'}
#'
#' ########################################################################
#' # Example 5: Real- world example; kNNDM in feature space
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(terra)
#' library(ggplot2)
#'
#'data(splotdata)
#'splotdata <- splotdata[splotdata$Country == "Chile",]
#'
#'predictors <- c("bio_1", "bio_4", "bio_5", "bio_6",
#' "bio_8", "bio_9", "bio_12", "bio_13",
#' "bio_14", "bio_15", "elev")
#'
#'trainDat <- sf::st_drop_geometry(splotdata)
#'predictors_sp <- terra::rast(system.file("extdata", "predictors_chile.tif",package="CAST"))
#'
#'
#' terra::plot(predictors_sp[["bio_1"]])
#' terra::plot(vect(splotdata), add = T)
#'
#'knndm_folds <- knndm(trainDat[,predictors], modeldomain = predictors_sp, dist_space = "feature",
#' clustering="kmeans", k=4, maxp=0.8)
#'plot(knndm_folds)
#'
#'}
knndm <- function(tpoints, modeldomain = NULL, predpoints = NULL,
dist_space = "geographical",
k = 10, maxp = 0.5,
clustering = "hierarchical", linkf = "ward.D2",
samplesize = 1000, sampling = "regular", dist_fun="euclidean",
algorithm="brute", scale_vars = TRUE,
space = NULL, useMD = NULL,
test_prop = NULL, test_tolerance = NULL,
nk_len = 100){
# Check for deprecated arguments
if (!is.null(space)) {
warning("Argument 'space' is deprecated. Please use 'dist_space' instead.",
call. = FALSE)
dist_space <- space
}
if (!is.null(useMD)) {
warning("Argument 'useMD' is deprecated. Please use 'dist_fun' instead.",
call. = FALSE)
}
if (dist_space == "geo") dist_space <- "geographical"
## Check that dist_space was correctly defined
if (!dist_space %in% c("geographical", "feature")) {
stop("dist_space must be one of 'geographical' or 'feature'")
}
if (!(dist_fun %in% c("euclidean", "mahalanobis", "gower", "great_circle"))) {
stop("dist_fun must be one of 'euclidean', 'mahalanobis', 'gower' or 'great_circle'")
}
if(dist_space == "time" && dist_fun != "euclidean") stop("Temporal space only supports euclidean distances.")
if(dist_space == "feature" && dist_fun == "great_circle") stop("Great-circle distances only work with in geographical space.")
if(dist_space == "geographical" && dist_fun %in% c("mahalanobis", "gower")) stop("Mahalanobis and Gower distances only work in feature space.")
# Issue a warning if train/test split is used
if(!is.null(test_prop)) {
warning("A train/test split will be returned, which is currently experimental.")
}
# Check that test_prop and test_tolerance are correctly specified and align parameters
minp <- NULL
if(!is.null(test_prop)) {
if(test_prop >= 1 | test_prop <= 0) {
stop("test_prop must be greater than 0 and smaller than 1")
}
if(is.null(test_tolerance)) test_tolerance <- 0.1
# adjust parameters
k <- 2
maxp <- test_prop + test_tolerance
minp <- test_prop - test_tolerance
if(maxp <= 0 || minp >= 1 || minp > maxp) {
stop("Misspecified tolerance. Resulted in infeasible minp/maxp values")
}
if(minp <= 0) {
warning("minp was set to 0.1")
minp <- 0.1
} else if(maxp >= 1) {
warning("maxp was set to 0.9")
maxp <- 0.9
}
if(maxp == 1/k) {
# adds some numerical tolerance to avoid maxp = 1/k
eps <- .Machine$double.eps^0.5
maxp <- maxp + eps
}
} else {
minp <- NULL
}
# Check nk_len
if(!nk_len%%1==0) {
stop("nk_len must be an integer")
}
# create sample points from modeldomain
if(is.null(predpoints)&!is.null(modeldomain)){
# Check modeldomain is indeed a sf/SpatRaster
if(!any(c("sfc", "sf", "SpatRaster") %in% class(modeldomain))){
stop("modeldomain must be a sf/sfc object or a 'SpatRaster' object.")
}
# If modeldomain is a SpatRaster, transform into polygon
if(any(class(modeldomain) == "SpatRaster")){
# save predictor stack for extraction if dist_space = "feature"
if(dist_space == "feature") {
predictor_stack <- modeldomain
}
modeldomain[!is.na(modeldomain)] <- 1
modeldomain <- terra::as.polygons(modeldomain, values = FALSE, na.all = TRUE) |>
sf::st_as_sf() |>
sf::st_union()
if(any(c("sfc", "sf") %in% class(tpoints))) {
modeldomain <- sf::st_transform(modeldomain, crs = sf::st_crs(tpoints))
}
}
# Check modeldomain is indeed a polygon sf
if(!any(class(sf::st_geometry(modeldomain)) %in% c("sfc_POLYGON", "sfc_MULTIPOLYGON"))){
stop("modeldomain must be a sf/sfc polygon object.")
}
# Check whether modeldomain has the same crs as tpoints
if(!identical(sf::st_crs(tpoints), sf::st_crs(modeldomain)) & dist_space == "geographical"){
stop("tpoints and modeldomain must have the same CRS")
}
# We sample
message(paste0(samplesize, " prediction points are sampled from the modeldomain"))
predpoints <- suppressMessages(sf::st_sample(x = modeldomain, size = samplesize, type = sampling))
sf::st_crs(predpoints) <- sf::st_crs(modeldomain)
if(dist_space == "feature") {
message("predictor values are extracted for prediction points")
predpoints <- terra::extract(predictor_stack, terra::vect(predpoints), ID=FALSE)
}
}else if(!is.null(predpoints) & dist_space == "geographical"){
if(!identical(sf::st_crs(tpoints), sf::st_crs(predpoints))){
stop("tpoints and predpoints must have the same CRS")
}
}
# Conditional preprocessing actions
if(dist_space == "geographical") {
if (any(class(tpoints) %in% "sfc")) {
tpoints <- sf::st_sf(geom = tpoints)
}
if (any(class(predpoints) %in% "sfc")) {
predpoints <- sf::st_sf(geom = predpoints)
}
if(is.na(sf::st_crs(tpoints))){
warning("Missing CRS in training or prediction points. Assuming projected CRS.")
islonglat <- FALSE
}else{
islonglat <- sf::st_is_longlat(tpoints)
}
} else if (dist_space == "feature") {
# drop geometry if tpoints / predpoints are of class sf
if(any(class(tpoints) %in% c("sf","sfc"))) {
tpoints <- sf::st_set_geometry(tpoints, NULL)
}
if(any(class(predpoints) %in% c("sf","sfc"))) {
predpoints <- sf::st_set_geometry(predpoints, NULL)
}
# get names of categorical variables
catVars <- names(tpoints)[vapply(tpoints, function(z) inherits(z, c("factor", "character")), logical(1))]
if (length(catVars) == 0) catVars <- NULL
if(!is.null(catVars)) {
message(paste0("variable(s) '", catVars, "' is (are) treated as categorical variables"))
}
# omit NAs
if(any(is.na(predpoints))) {
message("some prediction points contain NAs, which will be removed")
predpoints <- stats::na.omit(predpoints)
}
if(any(is.na(tpoints))) {
message("some training points contain NAs, which will be removed")
tpoints <- stats::na.omit(tpoints)
}
}
# kNNDM in the geographical / feature space
if(isTRUE(dist_space == "geographical")){
# prior checks
check_knndm_geo(tpoints = tpoints, predpoints = predpoints, dist_space = dist_space,
k = k, maxp = maxp, clustering = clustering, dist_fun = dist_fun, test_prop = test_prop, islonglat = islonglat)
# kNNDM in geographical space
knndm_res <- knndm_geo(tpoints = tpoints, predpoints = predpoints, k = k, maxp = maxp, minp = minp,
test_prop = test_prop, clustering = clustering, linkf = linkf, nk_len = nk_len,
dist_fun = dist_fun, dist_space = dist_space, algorithm = algorithm)
} else if (isTRUE(dist_space == "feature")) {
# prior checks
check_knndm_feature(tpoints = tpoints, predpoints = predpoints, dist_space = dist_space,
k = k, maxp = maxp, clustering = clustering, dist_fun = dist_fun, test_prop = test_prop, catVars = catVars)
# kNNDM in feature space
knndm_res <- knndm_feature(tpoints = tpoints, predpoints = predpoints, k = k, maxp = maxp, minp = minp,
test_prop = test_prop, clustering = clustering, linkf = linkf, nk_len = nk_len, dist_fun = dist_fun,
dist_space = dist_space, algorithm = algorithm, catVars = catVars, scale_vars = scale_vars)
}
# Output
knndm_res
}
# kNNDM checks
check_knndm_geo <- function(tpoints, predpoints, dist_space, k, maxp, clustering, islonglat, dist_fun, test_prop){
if(!identical(sf::st_crs(tpoints), sf::st_crs(predpoints))){
stop("tpoints and predpoints must have the same CRS")
}
if (!(clustering %in% c("kmeans", "hierarchical"))) {
stop("clustering must be one of `kmeans` or `hierarchical`")
}
if(is.null(test_prop)) {
if (!(maxp < 1 & maxp > 1/k)) {
stop("maxp must be strictly between 1/k and 1")
}
}
if(isTRUE(islonglat) & clustering == "kmeans"){
stop("kmeans works in the Euclidean space and therefore can only handle
projected coordinates. Please use hierarchical clustering or project your data.")
}
if(isTRUE(islonglat) && dist_fun != "great_circle") {
stop("Only great-circle distances are allowed for lon/lat coordinates. Please use 'great_circle' as 'dist_fun'.")
}
}
check_knndm_feature <- function(tpoints, predpoints, dist_space, k, maxp, clustering, catVars, dist_fun, test_prop){
if (!is.null(catVars) && dist_fun != "gower") {
stop("Only gower distances work with categorical features. Please use dist_fun = 'gower'")
}
if(is.null(test_prop)) {
if (!(maxp < 1 & maxp > 1/k)) {
stop("maxp must be strictly between 1/k and 1")
}
}
if(is.null(predpoints)) {
stop("predpoints with predictor data missing")
}
if(length(setdiff(names(tpoints), names(predpoints)))>0) {
stop("tpoints and predpoints need to contain the predictor data and have the same colnames.")
}
for (catvar in catVars) {
if (any(!unique(tpoints[,catvar]) %in% unique(predpoints[,catvar]))) {
stop(paste0("Some values of factor", catvar, "are only present in training / prediction points.
All factor values in the prediction points must be present in the training points."))
}
}
}
# kNNDM in the geographical space
knndm_geo <- function(tpoints, predpoints, k, maxp, minp, test_prop,
clustering, linkf, nk_len, dist_fun, dist_space, algorithm){
# Gj and Gij calculation
tcoords <- sf::st_coordinates(tpoints)[,1:2]
if(isTRUE(dist_fun == "great_circle")){
# For great-circle distance, we calculate the distance matrix here once and then use
# distclust_distmat later to avoid re-calculating the dist_mat when using compute_NND
distmat <- sf::st_distance(tpoints)
units(distmat) <- NULL
diag(distmat) <- NA
Gj <- apply(distmat, 1, function(x) min(x, na.rm=TRUE))
Gij <- sf::st_distance(predpoints, tpoints)
units(Gij) <- NULL
Gij <- apply(Gij, 1, min)
}else{
Gj <- compute_NND(tpoints, dist_space = dist_space, dist_fun = dist_fun, algorithm = algorithm)$dist
Gij <- compute_NND(tpoints, y = predpoints, dist_space = dist_space, dist_fun = dist_fun, algorithm = algorithm)$dist
}
# Check if Gj > Gij (warning suppressed regarding ties)
testks <- suppressWarnings(stats::ks.test(Gj, Gij, alternative = "great"))
if(testks$p.value >= 0.05){
if(!is.null(test_prop)) {
ntest <- floor(test_prop * nrow(tpoints))
ntrain <- nrow(tpoints) - ntest
# Create a vector: 1 = train, 2 = test
clusters <- c(rep(1, ntrain), rep(2, ntest))
# Shuffle randomly
clust <- sample(clusters, nrow(tpoints))
} else {
clust <- sample(rep(1:k, ceiling(nrow(tpoints)/k)), size = nrow(tpoints), replace=F)
}
if(isTRUE(dist_fun == "great_circle")){
Gjstar <- distclust_distmat(distmat, clust)
}else{
Gjstar <- cv_distances(tcoords, CVtest = clust, algorithm=algorithm, dist_fun = dist_fun)
}
k_final <- "random CV"
W_final <- twosamples::wass_stat(Gjstar, Gij)
message("Gij <= Gj; a random CV assignment is returned")
}else{
if(clustering == "hierarchical"){
# For hierarchical clustering we need to compute the full distance matrix,
# but we can integrate geographical distances
if(isTRUE(dist_fun == "euclidean")){
distmat <- sf::st_distance(tpoints)
}
hc <- stats::hclust(d = stats::as.dist(distmat), method = linkf)
}
# Build grid of number of clusters to try - we sample low numbers more intensively
clustgrid <- data.frame(nk = as.integer(round(exp(seq(log(k), log(nrow(tpoints)-2),
length.out = nk_len)))))
clustgrid$W <- NA
clustgrid <- clustgrid[!duplicated(clustgrid$nk),]
clustgroups <- list()
# Compute 1st PC for ordering clusters
pcacoords <- stats::prcomp(tcoords, center = TRUE, scale. = FALSE, rank = 1)
# We test each number of clusters
for(nk in clustgrid$nk){
# Create nk clusters
if(clustering == "hierarchical"){
clust_nk <- stats::cutree(hc, k=nk)
}else if(clustering == "kmeans"){
clust_nk <- stats::kmeans(tcoords, nk)$cluster
}
tabclust <- as.data.frame(table(clust_nk))
tabclust$clust_k <- NA
# compute cluster centroids and apply PC loadings to shuffle along the 1st dimension
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca <- matrix(colMeans(tcoords[clust_nk %in% x, , drop = FALSE]), nrow = 1)
colnames(centrpca) <- colnames(tcoords)
return(predict(pcacoords, centrpca))
},numeric(1))
tabclust$centrpca <- centr_tpoints
tabclust <- tabclust[order(tabclust$centrpca),]
# We don't merge big clusters
if(is.null(test_prop)) {
clust_i <- 1
for(i in 1:nrow(tabclust)){
if(tabclust$Freq[i] >= nrow(tpoints)/k){
tabclust$clust_k[i] <- clust_i
clust_i <- clust_i + 1
}
}
rm("clust_i")
}
# And we merge the remaining into k groups
clust_i <- setdiff(1:k, unique(tabclust$clust_k))
tabclust$clust_k[is.na(tabclust$clust_k)] <- rep(clust_i, ceiling(nk/length(clust_i)))[1:sum(is.na(tabclust$clust_k))]
tabclust2 <- data.frame(ID = 1:length(clust_nk), clust_nk = clust_nk)
tabclust2 <- merge(tabclust2, tabclust, by = "clust_nk")
tabclust2 <- tabclust2[order(tabclust2$ID),]
clust_k <- tabclust2$clust_k
# Check size of clust_k
if(is.null(test_prop)) {
prop_valid <- !(any(table(clust_k)/length(clust_k)>maxp))
} else {
# For train/test splits, only compute W if < maxp and > minp
# Calculate the proportion by group (train/test)
prop_1 <- mean(clust_k == 1)
prop_2 <- mean(clust_k == 2)
props <- c(prop_1, prop_2)
# Keep only groups within range minp–maxp
prop_valid <- props >= minp & props <= maxp
}
if(any(prop_valid)){
if(isTRUE(dist_fun == "great_circle")){
Gjstar_i <- distclust_distmat(distmat, clust_k)
}else{
Gjstar_i <- cv_distances(tcoords, CVtest = clust_k,algorithm=algorithm, dist_fun = dist_fun)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
# Compute W statistic if not exceeding maxp
if(!any(table(clust_k)/length(clust_k)>maxp)){
if(isTRUE(dist_fun == "great_circle")){
Gjstar_i <- distclust_distmat(distmat, clust_k)
}else{
Gjstar_i <- cv_distances(tcoords, CVtest = clust_k,algorithm=algorithm, dist_fun = dist_fun)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
}
# Final configuration
k_final <- clustgrid$nk[which.min(clustgrid$W)]
W_final <- min(clustgrid$W, na.rm=T)
clust <- clustgroups[[paste0("nk", k_final)]]
if(!is.null(test_prop) && is.null(clust)) {
stop("No valid train/test configurations found in the range test_prop +/- tolerance. Increase tolerance.")
}
if(isTRUE(dist_fun == "great_circle")){
Gjstar <- distclust_distmat(distmat, clust)
}else{
Gjstar <- cv_distances(tcoords, CVtest = clust,algorithm=algorithm, dist_fun = dist_fun)
}
}
# Output
if(is.null(test_prop)) {
cfolds <- CAST::CreateSpacetimeFolds(data.frame(clust=clust), spacevar = "clust", k = k)
} else {
# Assign train/test classes
deviation_1 <- abs((table(clust)[[1]] / length(clust)) - test_prop)
deviation_2 <- abs((table(clust)[[2]] / length(clust)) - test_prop)
if(deviation_1 > deviation_2) {
test_class <- 2
} else {
test_class <- 1
}
clust[clust == test_class] <- "test"
clust[clust != "test"] <- "train"
cfolds <- list("indexOut" = which(clust == "test"), "index" = which(clust == "train"))
}
res <- list(clusters = clust,
indx_train = cfolds$index, indx_test = cfolds$indexOut,
Gij = Gij, Gj = Gj, Gjstar = Gjstar,
W = W_final, method = clustering, q = k_final, dist_space = "geographical")
class(res) <- c("knndm", "list")
res
}
# kNNDM in the feature space
knndm_feature <- function(tpoints, predpoints, k, maxp, minp, test_prop,
clustering, linkf, nk_len, dist_fun, dist_space, algorithm, catVars, scale_vars) {
# rescale data (optional)
if(isTRUE(scale_vars)) {
if(is.null(catVars)) {
scale_attr <- attributes(scale(tpoints))
tpoints <- scale(tpoints) |> as.data.frame()
predpoints <- scale(predpoints,center=scale_attr$`scaled:center`,
scale=scale_attr$`scaled:scale`) |>
as.data.frame()
} else {
tpoints_cat <- tpoints[,catVars,drop=FALSE]
predpoints_cat <- predpoints[,catVars,drop=FALSE]
tpoints_num <- tpoints[,-which(names(tpoints)%in%catVars),drop=FALSE]
predpoints_num <- predpoints[,-which(names(predpoints)%in%catVars),drop=FALSE]
scale_attr <- attributes(scale(tpoints_num))
tpoints <- scale(tpoints_num) |> as.data.frame()
predpoints <- scale(predpoints_num,center=scale_attr$`scaled:center`,
scale=scale_attr$`scaled:scale`) |>
as.data.frame()
tpoints <- as.data.frame(cbind(tpoints, lapply(tpoints_cat, as.factor)))
predpoints <- as.data.frame(cbind(predpoints, lapply(predpoints_cat, as.factor)))
}
}
# Gj and Gij calculation
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
tpoints_mat <- as.matrix(tpoints)
predpoints_mat <- as.matrix(predpoints)
# use Mahalanobis distances
if (dim(tpoints_mat)[2] == 1) {
S <- matrix(stats::var(tpoints_mat), 1, 1)
tpoints_mat <- as.matrix(tpoints_mat, ncol = 1)
} else {
S <- stats::cov(tpoints_mat)
}
S_inv <- MASS::ginv(S)
# calculate distance matrix
n_rows <- nrow(tpoints_mat)
distmat <- vapply(seq_len(n_rows), function(i) {
vapply(seq_len(n_rows), function(j) {
diff <- tpoints_mat[i, ] - tpoints_mat[j, ]
sqrt(t(diff) %*% S_inv %*% diff)
}, numeric(1))
}, numeric(n_rows))
diag(distmat) <- NA
Gj <- apply(distmat, 1, min, na.rm=TRUE)
n_rows_p <- nrow(predpoints_mat)
n_rows_t <- nrow(tpoints_mat)
Gij <- vapply(seq_len(n_rows_p), function(i) {
min(vapply(seq_len(n_rows_t), function(j) {
diff <- predpoints_mat[i, ] - tpoints_mat[j, ]
sqrt(t(diff) %*% S_inv %*% diff)
}, numeric(1)))
}, numeric(1))
} else {
# use FNN with Euclidean distances if no categorical variables are present
Gj <- c(FNN::knn.dist(tpoints, k = 1, algorithm=algorithm))
Gij <- c(FNN::knnx.dist(query = predpoints, data = tpoints, k = 1, algorithm=algorithm))
}
} else {
# use Gower distances if categorical variables are present
Gj <- vapply(1:nrow(tpoints), function(i) gower::gower_topn(tpoints[i,], tpoints[-i,], n=1)$distance[[1]], numeric(1))
Gij <- c(gower::gower_topn(predpoints, tpoints, n = 1)$distance)
}
# Check if Gj > Gij (warning suppressed regarding ties)
testks <- suppressWarnings(stats::ks.test(Gj, Gij, alternative = "great"))
if(testks$p.value >= 0.05){
if(!is.null(test_prop)) {
ntest <- floor(test_prop * nrow(tpoints))
ntrain <- nrow(tpoints) - ntest
# Create a vector: 1 = train, 2 = test
clusters <- c(rep(1, ntrain), rep(2, ntest))
# Shuffle randomly
clust <- sample(clusters, nrow(tpoints))
} else {
clust <- sample(rep(1:k, ceiling(nrow(tpoints)/k)), size = nrow(tpoints), replace=F)
}
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
} else {
Gjstar <- cv_distances(tpoints, clust,algorithm=algorithm, dist_fun = dist_fun)
}
} else {
Gjstar <- cv_distances(tpoints, clust, dist_fun = dist_fun)
}
k_final <- "random CV"
W_final <- twosamples::wass_stat(Gjstar, Gij)
message("Gij <= Gj; a random CV assignment is returned")
}else{
if(clustering == "hierarchical"){
# calculate distance matrix which is needed for hierarchical clustering
if(is.null(catVars)) {
if(isTRUE(dist_fun == "euclidean")) {
# calculate distance matrix with Euclidean distances if no categorical variables are present
# for MD: distance matrix was already calculated
distmat <- stats::dist(tpoints, upper=TRUE, diag=TRUE) |> as.matrix()
diag(distmat) <- NA
}
} else {
# calculate distance matrix with Gower distances if categorical variables are present
distmat <- matrix(nrow=nrow(tpoints), ncol=nrow(tpoints))
for (i in 1:nrow(tpoints)){
trainDist <- gower::gower_dist(tpoints[i,], tpoints)
trainDist[i] <- NA
distmat[i,] <- trainDist
}
}
hc <- stats::hclust(d = stats::as.dist(distmat), method = linkf)
}
# Build grid of number of clusters to try - we sample low numbers more intensively
clustgrid <- data.frame(nk = as.integer(round(exp(seq(log(k), log(nrow(tpoints)-2),
length.out = nk_len)))))
clustgrid$W <- NA
clustgrid <- clustgrid[!duplicated(clustgrid$nk),]
clustgroups <- list()
# Compute 1st PC for ordering clusters
if(is.null(catVars)) {
pcacoords <- stats::prcomp(tpoints, center = TRUE, scale. = FALSE, rank = 1)
} else {
pcacoords <- PCAmixdata::PCAmix(X.quanti = tpoints[,!(names(tpoints) %in% catVars), drop=FALSE],
X.quali = tpoints[,names(tpoints) %in% catVars, drop=FALSE],
graph = FALSE)
}
# We test each number of clusters
for(nk in clustgrid$nk) {
# Create nk clusters
if(clustering == "hierarchical"){
clust_nk <- stats::cutree(hc, k=nk)
} else if(clustering == "kmeans"){
if(is.null(catVars)) {
clust_nk <- tryCatch(stats::kmeans(tpoints, nk)$cluster,
error=function(e) e)
} else {
# prototype clustering for mixed data sets
clust_nk <- tryCatch(clustMixType::kproto(tpoints, nk,verbose=FALSE)$cluster,
error=function(e) e)
}
}
if (!inherits(clust_nk,"error")){
tabclust <- as.data.frame(table(clust_nk))
tabclust$clust_k <- NA
# compute cluster centroids and apply PC loadings to shuffle along the 1st dimension
if(is.null(catVars)) {
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca <- matrix(colMeans(tpoints[clust_nk %in% x, , drop = FALSE]), nrow = 1)
colnames(centrpca) <- colnames(tpoints)
return(predict(pcacoords, centrpca))
},numeric(1))
} else {
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca_num <- matrix(apply(tpoints[clust_nk %in% x, !(names(tpoints) %in% catVars), drop=FALSE], 2, mean), nrow = 1)
centrpca_cat <- matrix(apply(tpoints[clust_nk %in% x, names(tpoints) %in% catVars, drop=FALSE], 2,
function(y) names(which.max(table(y)))), nrow = 1)
colnames(centrpca_num) <- colnames(tpoints[,!(names(tpoints) %in% catVars), drop=FALSE])
colnames(centrpca_cat) <- colnames(tpoints[,names(tpoints) %in% catVars, drop=FALSE])
return(predict(pcacoords, centrpca_num, centrpca_cat)[,1])
}, numeric(1))
}
tabclust$centrpca <- centr_tpoints
tabclust <- tabclust[order(tabclust$centrpca),]
# We don't merge big clusters
if(is.null(test_prop)) {
clust_i <- 1
for(i in 1:nrow(tabclust)){
if(tabclust$Freq[i] >= nrow(tpoints)/k){
tabclust$clust_k[i] <- clust_i
clust_i <- clust_i + 1
}
}
rm("clust_i")
}
# And we merge the remaining into k groups
clust_i <- setdiff(1:k, unique(tabclust$clust_k))
tabclust$clust_k[is.na(tabclust$clust_k)] <- rep(clust_i, ceiling(nk/length(clust_i)))[1:sum(is.na(tabclust$clust_k))]
tabclust2 <- data.frame(ID = 1:length(clust_nk), clust_nk = clust_nk)
tabclust2 <- merge(tabclust2, tabclust, by = "clust_nk")
tabclust2 <- tabclust2[order(tabclust2$ID),]
clust_k <- tabclust2$clust_k
# Check size of clust_k
if(is.null(test_prop)) {
prop_valid <- !(any(table(clust_k)/length(clust_k)>maxp))
} else {
# For train/test splits, only compute W if < maxp and > minp
# Calculate the proportion by group (train/test)
prop_1 <- mean(clust_k == 1)
prop_2 <- mean(clust_k == 2)
props <- c(prop_1, prop_2)
# Keep only groups within range minp–maxp
prop_valid <- props >= minp & props <= maxp
}
# Compute W statistic if size of clust_k is valid
if(any(prop_valid)){
if(clustering == "kmeans") {
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")){
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k, dist_fun = dist_fun)
} else {
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k,algorithm = algorithm, dist_fun = dist_fun)
}
} else {
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k, dist_fun = dist_fun)
}
} else {
Gjstar_i <- distclust_distmat(distmat, clust_k)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
}
}
# Final configuration
k_final <- clustgrid$nk[which.min(clustgrid$W)]
W_final <- min(clustgrid$W, na.rm=T)
clust <- clustgroups[[paste0("nk", k_final)]]
if(!is.null(test_prop) && is.null(clust)) {
stop("No valid train/test configurations found in the range test_prop +/- tolerance. Increase tolerance.")
}
if(clustering == "kmeans") {
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
} else {
Gjstar <- cv_distances(tpoints, CVtest = clust,algorithm=algorithm, dist_fun = dist_fun)
}
} else {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
}
} else {
Gjstar <- distclust_distmat(distmat, clust)
}
}
# Output
if(is.null(test_prop)) {
cfolds <- CAST::CreateSpacetimeFolds(data.frame(clust=clust), spacevar = "clust", k = k)
} else {
# Assign train/test classes
deviation_1 <- abs((table(clust)[[1]] / length(clust)) - test_prop)
deviation_2 <- abs((table(clust)[[2]] / length(clust)) - test_prop)
if(deviation_1 > deviation_2) {
test_class <- 2
} else {
test_class <- 1
}
clust[clust == test_class] <- "test"
clust[clust != "test"] <- "train"
cfolds <- list("indexOut" = which(clust == "test"), "index" = which(clust == "train"))
}
res <- list(clusters = clust,
indx_train = cfolds$index, indx_test = cfolds$indexOut,
Gij = Gij, Gj = Gj, Gjstar = Gjstar,
W = W_final, method = clustering, q = k_final, dist_space = "feature")
class(res) <- c("knndm", "list")
res
}
# Helper function: Compute out-of-fold NN distance based on a distance matrix (geographical coordinates / numerical variables)
distclust_distmat <- function(distm, folds){
alldist <- rep(NA, length(folds))
for(f in unique(folds)){
alldist[f == folds] <- apply(distm[f == folds, f != folds, drop=FALSE], 1, min)
}
alldist
}
Try the CAST package in your browser
Any scripts or data that you put into this service are public.
CAST documentation built on Sept. 9, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions distclust_distmat knndm_feature knndm_geo check_knndm_feature check_knndm_geo knndm
Documented in knndm
#' K-fold Nearest Neighbour Distance Matching
#' @description
#' This function implements the kNNDM algorithm for prediction-domain adaptive resampling
#' and returns the necessary indices to perform train-test splits or k-fold NNDM CV.
#'
#' @author Carles Milà and Jan Linnenbrink
#' @param tpoints sf or sfc point object, or data.frame if dist_space = "feature". Contains the training points samples.
#' @param modeldomain sf polygon object or SpatRaster defining the prediction area. Optional; alternative to predpoints (see Details).
#' @param predpoints sf or sfc point object, or data.frame if dist_space = "feature". Contains the target prediction points. Optional; alternative to modeldomain (see Details).
#' @param dist_space character. Either "geographical" or "feature".
#' @param k integer. Number of folds desired for CV. Defaults to 10.
#' @param maxp numeric. Maximum fold size allowed, defaults to 0.5, i.e. a single fold can hold a maximum of half of the training points.
#' @param clustering character. Possible values include "hierarchical" and "kmeans". See details.
#' @param linkf character. Only relevant if clustering = "hierarchical". Link function for agglomerative hierarchical clustering.
#' Defaults to "ward.D2". Check `stats::hclust` for other options.
#' @param samplesize numeric. How many points in the modeldomain should be sampled as prediction points?
#' Only required if modeldomain is used instead of predpoints.
#' @param sampling character. How to draw prediction points from the modeldomain? See `sf::st_sample`.
#' Only required if modeldomain is used instead of predpoints.
#' @param dist_fun character. Currently covers `euclidean` (default), `gower`, `mahalanobis` and `great_circle`.
#' `gower` and `mahalanobis` only work with `dist_space`="feature", while `great_circle` only works with `dist_space`="geographical".
#' `mahalanobis` takes into account correlation between predictor values. While `euclidean` and `mahalanobis` only work with numerical variables,
#' `gower` also works with mixed data including numerical and categorical variables.
#' For the geographical space, `great_circle` covers lon/lat coordinates, whereas `euclidean` only works with projected coordinates.
#' @param algorithm see \code{\link[FNN]{knnx.dist}} and \code{\link[FNN]{knnx.index}}
#' @param scale_vars boolean. Should variables be scaled? Only for `dist_space`="feature".
#' Calculating Gower distances already includes scaling, and manually rescale the data is redundant.
#' For other distances (Mahalanobis, Euclidean), scaling the data is important. Thus, TRUE by default.
#' @param test_prop numeric. The proportion of test data. NULL by default (i.e., no train/test split).
#' @param test_tolerance numeric. The allowed deviance from `test_prop`. The higher the tolerance,
#' the larger the possibility to obtain train/test splits that yield good approximations of the prediction situation.
#' @param nk_len integer. The number of fold configurations to test. By default 100.
#' Larger numbers increase computational times, but also might lead to better W statistics.
#' Useful for train/test splits, where a large number of configurations is discarded.
#' @param space deprecated. Use `dist_space` instead.
#' @param useMD deprecated. Use `dist_fun` instead.
#' @return An object of class \emph{knndm} consisting of a list of eight elements:
#' indx_train, indx_test (indices of the observations to use as
#' training/test data in each kNNDM CV iteration), Gij (distances for
#' G function construction between prediction and target points), Gj
#' (distances for G function construction during LOO CV), Gjstar (distances
#' for modified G function during kNNDM CV), clusters (list of cluster IDs),
#' W (Wasserstein statistic), and dist_space (stated by the user in the function call).
#'
#' @details
#' knndm is an implementation of prediction-domain adaptive validation.
#' It is a k-fold version of NNDM LOO CV which makes it more suitable for medium and large datasets.
#' It can be used for cross-validation and train / test splits (the latter is experimental).
#' Briefly, the algorithm tries to find a configuration such that the integral of the absolute differences (Wasserstein W statistic)
#' between the empirical nearest neighbour distance distribution function between the test and training data (Gj*),
#' and the empirical nearest neighbour distance distribution function between the prediction and training points (Gij),
#' is minimised. It does so by performing clustering of the training points' coordinates for different numbers of
#' clusters that range from k to N (number of observations), merging them into k final folds,
#' and selecting the configuration with the lowest W.
#'
#' When using `knndm` to split the data into training and test sets (experimental), the proportion of points belonging to the test set (`test_prop`) replaces the number of folds `k`.
#' Based on the `test_prop` , `minp` and `maxp` are calculated as the `test_prop` +/- `test_tolerance`.
#' Compared to k-fold CV, using knndm for train/test splits is less flexible and often results in larger NNDs between test and train locations
#' than between prediction and train locations. Hence, it is essential to plot the results of `knndm` and check how well the split can resemble the prediction situation.
#' Modifying the `test_prop` parameter, as well as increasing `test_prop` allow more flexible matching and can potentially improve the match.
#'
#' Using a projected CRS in `knndm` has large computational advantages since fast nearest neighbour search can be
#' done via the `FNN` package, while working with geographic coordinates requires computing the full
#' spherical distance matrices. As a clustering algorithm, `kmeans` can only be used for
#' projected CRS while `hierarchical` can work with both projected and geographical coordinates, though it requires
#' calculating the full distance matrix of the training points even for a projected CRS.
#'
#' In order to select between clustering algorithms and number of folds `k`, different `knndm` configurations can be run
#' and compared, being the one with a lower W statistic the one that offers a better match. W statistics between `knndm`
#' runs are comparable as long as `tpoints` and `predpoints` or `modeldomain` stay the same.
#'
#' Map validation using `knndm` should be used using `CAST::global_validation`, i.e. by stacking all out-of-sample
#' predictions and evaluating them all at once. The reasons behind this are 1) The resulting folds can be
#' unbalanced and 2) nearest neighbour functions are constructed and matched using all CV folds simultaneously.
#'
#' If training data points are very clustered with respect to the prediction area and the presented `knndm`
#' configuration still show signs of Gj* > Gij, there are several things that can be tried. First, increase
#' the `maxp` parameter; this may help to control for strong clustering (at the cost of having unbalanced folds).
#' Secondly, decrease the number of final folds `k`, which may help to have larger clusters.
#'
#' The `modeldomain` is either a sf polygon that defines the prediction area, or alternatively a SpatRaster out of which a polygon,
#' transformed into the CRS of the training points, is defined as the outline of all non-NA cells.
#' Then, the function takes a regular point sample (amount defined by `samplesize`) from the spatial extent.
#' As an alternative use `predpoints` instead of `modeldomain`, if you have already defined the prediction locations (e.g. raster pixel centroids).
#' When using either `modeldomain` or `predpoints`, we advise to plot the study area polygon and the training/prediction points as a previous step to ensure they are aligned.
#'
#' `knndm` can also be performed in the feature space by setting `dist_space` to "feature".
#' Euclidean distances, Gower distance or Mahalanobis distances can be used for distance calculation, but only Euclidean are tested.
#' In this case, nearest neighbour distances are calculated in n-dimensional feature space rather than in geographical space.
#' `tpoints` and `predpoints` can be data frames or sf objects containing the values of the features. Note that the names of `tpoints` and `predpoints` must be the same.
#' `predpoints` can also be missing, if `modeldomain` is of class SpatRaster. In this case, the values of of the SpatRaster will be extracted to the `predpoints`.
#' In the case of any categorical features, Gower distances will be used to calculate the Nearest Neighbour distances [Experimental]. If categorical
#' features are present, and `clustering` = "kmeans", K-Prototype clustering will be performed instead.
#'
#' @note
#' For spatial visualization of fold affiliation see examples.
#' @references
#' \itemize{
#' \item Linnenbrink, J., Milà, C., Ludwig, M., and Meyer, H. (2024): kNNDM: k-fold Nearest Neighbour Distance Matching Cross-Validation for map accuracy estimation. Geosci. Model Dev., 17, 5897–5912. https://doi.org/10.5194/gmd-17-5897-2024.
#' \item 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, 1304– 1316. https://doi.org/10.1111/2041-210X.13851.
#' }
#' @seealso \code{\link{geodist}}, \code{\link{nndm}}
#'
#' @export
#' @examples
#' ########################################################################
#' # Example 1: Simulated data - Randomly-distributed training points
#' ########################################################################
#'
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 random training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- sf::st_sample(simarea, 1000, type = "random")
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Run kNNDM for the whole domain, here the prediction points are known.
#' knndm_folds <- knndm(train_points, predpoints = pred_points, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' plot(knndm_folds, type = "simple") # For more accessible legend labels
#' plot(knndm_folds, type = "simple", stat = "density") # To visualize densities rather than ECDFs
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'
#' ########################################################################
#' # Example 2: Simulated data - Clustered training points
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 clustered training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- clustered_sample(simarea, 1000, 50, 5)
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Run kNNDM for the whole domain, here the prediction points are known.
#' knndm_folds <- knndm(train_points, predpoints = pred_points, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' plot(knndm_folds, type = "simple") # For more accessible legend labels
#' plot(knndm_folds, type = "simple", stat = "density") # To visualize densities rather than ECDFs
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'}
#' ########################################################################
#' # Example 3: Real- world example; using a modeldomain instead of previously
#' # sampled prediction locations
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(terra)
#' library(ggplot2)
#'
#' ### prepare sample data:
#' data(cookfarm)
#' dat <- aggregate(cookfarm[,c("DEM","TWI", "NDRE.M", "Easting", "Northing","VW")],
#' by=list(as.character(cookfarm$SOURCEID)),mean)
#' pts <- dat[,-1]
#' pts <- st_as_sf(pts,coords=c("Easting","Northing"))
#' st_crs(pts) <- 26911
#' studyArea <- rast(system.file("extdata","predictors_2012-03-25.tif",package="CAST"))
#' pts <- st_transform(pts, crs = st_crs(studyArea))
#' terra::plot(studyArea[["DEM"]])
#' terra::plot(vect(pts), add = T)
#'
#' knndm_folds <- knndm(pts, modeldomain=studyArea, k = 5)
#' knndm_folds
#' plot(knndm_folds)
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = pts, aes(col = folds))
#'
#' #use for cross-validation:
#' library(caret)
#' ctrl <- trainControl(method="cv",
#' index=knndm_folds$indx_train,
#' savePredictions='final')
#' model_knndm <- train(dat[,c("DEM","TWI", "NDRE.M")],
#' dat$VW,
#' method="rf",
#' trControl = ctrl)
#' global_validation(model_knndm)
#'}
#' ########################################################################
#' # Example 4: Simulated data - Train/test split with clustered training points
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(ggplot2)
#'
#' # Simulate 1000 clustered training points in a 100x100 square
#' set.seed(1234)
#' simarea <- list(matrix(c(0,0,0,100,100,100,100,0,0,0), ncol=2, byrow=TRUE))
#' simarea <- sf::st_polygon(simarea)
#' train_points <- clustered_sample(simarea, 1000, 50, 5)
#' pred_points <- sf::st_sample(simarea, 1000, type = "regular")
#' plot(simarea)
#' plot(pred_points, add = TRUE, col = "blue")
#' plot(train_points, add = TRUE, col = "red")
#'
#' # Use kNNDM to split the data into 30% +- 10% test and 70% train
#' knndm_folds <- knndm(train_points, predpoints = pred_points, test_prop = 0.3, test_tolerance = 0.1)
#' # How many samples have been used for testing:
#' table(knndm_folds$clusters)
#' plot(knndm_folds)
#' # The train/test split could not represent the prediction situation well
#' # Increase tolerance to increase number of configurations tried, and thus to find a suitable split
#' knndm_folds <- knndm(train_points, predpoints = pred_points, test_prop = 0.3, test_tolerance = 0.2)
#' plot(knndm_folds)
#' table(knndm_folds$clusters)
#' # This resulted in better match of the prediction situation, but a 50/50 split
#' folds <- as.character(knndm_folds$clusters)
#' ggplot() +
#' geom_sf(data = simarea, alpha = 0) +
#' geom_sf(data = train_points, aes(col = folds))
#'}
#'
#' ########################################################################
#' # Example 5: Real- world example; kNNDM in feature space
#' ########################################################################
#' \dontrun{
#' library(sf)
#' library(terra)
#' library(ggplot2)
#'
#'data(splotdata)
#'splotdata <- splotdata[splotdata$Country == "Chile",]
#'
#'predictors <- c("bio_1", "bio_4", "bio_5", "bio_6",
#' "bio_8", "bio_9", "bio_12", "bio_13",
#' "bio_14", "bio_15", "elev")
#'
#'trainDat <- sf::st_drop_geometry(splotdata)
#'predictors_sp <- terra::rast(system.file("extdata", "predictors_chile.tif",package="CAST"))
#'
#'
#' terra::plot(predictors_sp[["bio_1"]])
#' terra::plot(vect(splotdata), add = T)
#'
#'knndm_folds <- knndm(trainDat[,predictors], modeldomain = predictors_sp, dist_space = "feature",
#' clustering="kmeans", k=4, maxp=0.8)
#'plot(knndm_folds)
#'
#'}
knndm <- function(tpoints, modeldomain = NULL, predpoints = NULL,
dist_space = "geographical",
k = 10, maxp = 0.5,
clustering = "hierarchical", linkf = "ward.D2",
samplesize = 1000, sampling = "regular", dist_fun="euclidean",
algorithm="brute", scale_vars = TRUE,
space = NULL, useMD = NULL,
test_prop = NULL, test_tolerance = NULL,
nk_len = 100){
# Check for deprecated arguments
if (!is.null(space)) {
warning("Argument 'space' is deprecated. Please use 'dist_space' instead.",
call. = FALSE)
dist_space <- space
}
if (!is.null(useMD)) {
warning("Argument 'useMD' is deprecated. Please use 'dist_fun' instead.",
call. = FALSE)
}
if (dist_space == "geo") dist_space <- "geographical"
## Check that dist_space was correctly defined
if (!dist_space %in% c("geographical", "feature")) {
stop("dist_space must be one of 'geographical' or 'feature'")
}
if (!(dist_fun %in% c("euclidean", "mahalanobis", "gower", "great_circle"))) {
stop("dist_fun must be one of 'euclidean', 'mahalanobis', 'gower' or 'great_circle'")
}
if(dist_space == "time" && dist_fun != "euclidean") stop("Temporal space only supports euclidean distances.")
if(dist_space == "feature" && dist_fun == "great_circle") stop("Great-circle distances only work with in geographical space.")
if(dist_space == "geographical" && dist_fun %in% c("mahalanobis", "gower")) stop("Mahalanobis and Gower distances only work in feature space.")
# Issue a warning if train/test split is used
if(!is.null(test_prop)) {
warning("A train/test split will be returned, which is currently experimental.")
}
# Check that test_prop and test_tolerance are correctly specified and align parameters
minp <- NULL
if(!is.null(test_prop)) {
if(test_prop >= 1 | test_prop <= 0) {
stop("test_prop must be greater than 0 and smaller than 1")
}
if(is.null(test_tolerance)) test_tolerance <- 0.1
# adjust parameters
k <- 2
maxp <- test_prop + test_tolerance
minp <- test_prop - test_tolerance
if(maxp <= 0 || minp >= 1 || minp > maxp) {
stop("Misspecified tolerance. Resulted in infeasible minp/maxp values")
}
if(minp <= 0) {
warning("minp was set to 0.1")
minp <- 0.1
} else if(maxp >= 1) {
warning("maxp was set to 0.9")
maxp <- 0.9
}
if(maxp == 1/k) {
# adds some numerical tolerance to avoid maxp = 1/k
eps <- .Machine$double.eps^0.5
maxp <- maxp + eps
}
} else {
minp <- NULL
}
# Check nk_len
if(!nk_len%%1==0) {
stop("nk_len must be an integer")
}
# create sample points from modeldomain
if(is.null(predpoints)&!is.null(modeldomain)){
# Check modeldomain is indeed a sf/SpatRaster
if(!any(c("sfc", "sf", "SpatRaster") %in% class(modeldomain))){
stop("modeldomain must be a sf/sfc object or a 'SpatRaster' object.")
}
# If modeldomain is a SpatRaster, transform into polygon
if(any(class(modeldomain) == "SpatRaster")){
# save predictor stack for extraction if dist_space = "feature"
if(dist_space == "feature") {
predictor_stack <- modeldomain
}
modeldomain[!is.na(modeldomain)] <- 1
modeldomain <- terra::as.polygons(modeldomain, values = FALSE, na.all = TRUE) |>
sf::st_as_sf() |>
sf::st_union()
if(any(c("sfc", "sf") %in% class(tpoints))) {
modeldomain <- sf::st_transform(modeldomain, crs = sf::st_crs(tpoints))
}
}
# Check modeldomain is indeed a polygon sf
if(!any(class(sf::st_geometry(modeldomain)) %in% c("sfc_POLYGON", "sfc_MULTIPOLYGON"))){
stop("modeldomain must be a sf/sfc polygon object.")
}
# Check whether modeldomain has the same crs as tpoints
if(!identical(sf::st_crs(tpoints), sf::st_crs(modeldomain)) & dist_space == "geographical"){
stop("tpoints and modeldomain must have the same CRS")
}
# We sample
message(paste0(samplesize, " prediction points are sampled from the modeldomain"))
predpoints <- suppressMessages(sf::st_sample(x = modeldomain, size = samplesize, type = sampling))
sf::st_crs(predpoints) <- sf::st_crs(modeldomain)
if(dist_space == "feature") {
message("predictor values are extracted for prediction points")
predpoints <- terra::extract(predictor_stack, terra::vect(predpoints), ID=FALSE)
}
}else if(!is.null(predpoints) & dist_space == "geographical"){
if(!identical(sf::st_crs(tpoints), sf::st_crs(predpoints))){
stop("tpoints and predpoints must have the same CRS")
}
}
# Conditional preprocessing actions
if(dist_space == "geographical") {
if (any(class(tpoints) %in% "sfc")) {
tpoints <- sf::st_sf(geom = tpoints)
}
if (any(class(predpoints) %in% "sfc")) {
predpoints <- sf::st_sf(geom = predpoints)
}
if(is.na(sf::st_crs(tpoints))){
warning("Missing CRS in training or prediction points. Assuming projected CRS.")
islonglat <- FALSE
}else{
islonglat <- sf::st_is_longlat(tpoints)
}
} else if (dist_space == "feature") {
# drop geometry if tpoints / predpoints are of class sf
if(any(class(tpoints) %in% c("sf","sfc"))) {
tpoints <- sf::st_set_geometry(tpoints, NULL)
}
if(any(class(predpoints) %in% c("sf","sfc"))) {
predpoints <- sf::st_set_geometry(predpoints, NULL)
}
# get names of categorical variables
catVars <- names(tpoints)[vapply(tpoints, function(z) inherits(z, c("factor", "character")), logical(1))]
if (length(catVars) == 0) catVars <- NULL
if(!is.null(catVars)) {
message(paste0("variable(s) '", catVars, "' is (are) treated as categorical variables"))
}
# omit NAs
if(any(is.na(predpoints))) {
message("some prediction points contain NAs, which will be removed")
predpoints <- stats::na.omit(predpoints)
}
if(any(is.na(tpoints))) {
message("some training points contain NAs, which will be removed")
tpoints <- stats::na.omit(tpoints)
}
}
# kNNDM in the geographical / feature space
if(isTRUE(dist_space == "geographical")){
# prior checks
check_knndm_geo(tpoints = tpoints, predpoints = predpoints, dist_space = dist_space,
k = k, maxp = maxp, clustering = clustering, dist_fun = dist_fun, test_prop = test_prop, islonglat = islonglat)
# kNNDM in geographical space
knndm_res <- knndm_geo(tpoints = tpoints, predpoints = predpoints, k = k, maxp = maxp, minp = minp,
test_prop = test_prop, clustering = clustering, linkf = linkf, nk_len = nk_len,
dist_fun = dist_fun, dist_space = dist_space, algorithm = algorithm)
} else if (isTRUE(dist_space == "feature")) {
# prior checks
check_knndm_feature(tpoints = tpoints, predpoints = predpoints, dist_space = dist_space,
k = k, maxp = maxp, clustering = clustering, dist_fun = dist_fun, test_prop = test_prop, catVars = catVars)
# kNNDM in feature space
knndm_res <- knndm_feature(tpoints = tpoints, predpoints = predpoints, k = k, maxp = maxp, minp = minp,
test_prop = test_prop, clustering = clustering, linkf = linkf, nk_len = nk_len, dist_fun = dist_fun,
dist_space = dist_space, algorithm = algorithm, catVars = catVars, scale_vars = scale_vars)
}
# Output
knndm_res
}
# kNNDM checks
check_knndm_geo <- function(tpoints, predpoints, dist_space, k, maxp, clustering, islonglat, dist_fun, test_prop){
if(!identical(sf::st_crs(tpoints), sf::st_crs(predpoints))){
stop("tpoints and predpoints must have the same CRS")
}
if (!(clustering %in% c("kmeans", "hierarchical"))) {
stop("clustering must be one of `kmeans` or `hierarchical`")
}
if(is.null(test_prop)) {
if (!(maxp < 1 & maxp > 1/k)) {
stop("maxp must be strictly between 1/k and 1")
}
}
if(isTRUE(islonglat) & clustering == "kmeans"){
stop("kmeans works in the Euclidean space and therefore can only handle
projected coordinates. Please use hierarchical clustering or project your data.")
}
if(isTRUE(islonglat) && dist_fun != "great_circle") {
stop("Only great-circle distances are allowed for lon/lat coordinates. Please use 'great_circle' as 'dist_fun'.")
}
}
check_knndm_feature <- function(tpoints, predpoints, dist_space, k, maxp, clustering, catVars, dist_fun, test_prop){
if (!is.null(catVars) && dist_fun != "gower") {
stop("Only gower distances work with categorical features. Please use dist_fun = 'gower'")
}
if(is.null(test_prop)) {
if (!(maxp < 1 & maxp > 1/k)) {
stop("maxp must be strictly between 1/k and 1")
}
}
if(is.null(predpoints)) {
stop("predpoints with predictor data missing")
}
if(length(setdiff(names(tpoints), names(predpoints)))>0) {
stop("tpoints and predpoints need to contain the predictor data and have the same colnames.")
}
for (catvar in catVars) {
if (any(!unique(tpoints[,catvar]) %in% unique(predpoints[,catvar]))) {
stop(paste0("Some values of factor", catvar, "are only present in training / prediction points.
All factor values in the prediction points must be present in the training points."))
}
}
}
# kNNDM in the geographical space
knndm_geo <- function(tpoints, predpoints, k, maxp, minp, test_prop,
clustering, linkf, nk_len, dist_fun, dist_space, algorithm){
# Gj and Gij calculation
tcoords <- sf::st_coordinates(tpoints)[,1:2]
if(isTRUE(dist_fun == "great_circle")){
# For great-circle distance, we calculate the distance matrix here once and then use
# distclust_distmat later to avoid re-calculating the dist_mat when using compute_NND
distmat <- sf::st_distance(tpoints)
units(distmat) <- NULL
diag(distmat) <- NA
Gj <- apply(distmat, 1, function(x) min(x, na.rm=TRUE))
Gij <- sf::st_distance(predpoints, tpoints)
units(Gij) <- NULL
Gij <- apply(Gij, 1, min)
}else{
Gj <- compute_NND(tpoints, dist_space = dist_space, dist_fun = dist_fun, algorithm = algorithm)$dist
Gij <- compute_NND(tpoints, y = predpoints, dist_space = dist_space, dist_fun = dist_fun, algorithm = algorithm)$dist
}
# Check if Gj > Gij (warning suppressed regarding ties)
testks <- suppressWarnings(stats::ks.test(Gj, Gij, alternative = "great"))
if(testks$p.value >= 0.05){
if(!is.null(test_prop)) {
ntest <- floor(test_prop * nrow(tpoints))
ntrain <- nrow(tpoints) - ntest
# Create a vector: 1 = train, 2 = test
clusters <- c(rep(1, ntrain), rep(2, ntest))
# Shuffle randomly
clust <- sample(clusters, nrow(tpoints))
} else {
clust <- sample(rep(1:k, ceiling(nrow(tpoints)/k)), size = nrow(tpoints), replace=F)
}
if(isTRUE(dist_fun == "great_circle")){
Gjstar <- distclust_distmat(distmat, clust)
}else{
Gjstar <- cv_distances(tcoords, CVtest = clust, algorithm=algorithm, dist_fun = dist_fun)
}
k_final <- "random CV"
W_final <- twosamples::wass_stat(Gjstar, Gij)
message("Gij <= Gj; a random CV assignment is returned")
}else{
if(clustering == "hierarchical"){
# For hierarchical clustering we need to compute the full distance matrix,
# but we can integrate geographical distances
if(isTRUE(dist_fun == "euclidean")){
distmat <- sf::st_distance(tpoints)
}
hc <- stats::hclust(d = stats::as.dist(distmat), method = linkf)
}
# Build grid of number of clusters to try - we sample low numbers more intensively
clustgrid <- data.frame(nk = as.integer(round(exp(seq(log(k), log(nrow(tpoints)-2),
length.out = nk_len)))))
clustgrid$W <- NA
clustgrid <- clustgrid[!duplicated(clustgrid$nk),]
clustgroups <- list()
# Compute 1st PC for ordering clusters
pcacoords <- stats::prcomp(tcoords, center = TRUE, scale. = FALSE, rank = 1)
# We test each number of clusters
for(nk in clustgrid$nk){
# Create nk clusters
if(clustering == "hierarchical"){
clust_nk <- stats::cutree(hc, k=nk)
}else if(clustering == "kmeans"){
clust_nk <- stats::kmeans(tcoords, nk)$cluster
}
tabclust <- as.data.frame(table(clust_nk))
tabclust$clust_k <- NA
# compute cluster centroids and apply PC loadings to shuffle along the 1st dimension
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca <- matrix(colMeans(tcoords[clust_nk %in% x, , drop = FALSE]), nrow = 1)
colnames(centrpca) <- colnames(tcoords)
return(predict(pcacoords, centrpca))
},numeric(1))
tabclust$centrpca <- centr_tpoints
tabclust <- tabclust[order(tabclust$centrpca),]
# We don't merge big clusters
if(is.null(test_prop)) {
clust_i <- 1
for(i in 1:nrow(tabclust)){
if(tabclust$Freq[i] >= nrow(tpoints)/k){
tabclust$clust_k[i] <- clust_i
clust_i <- clust_i + 1
}
}
rm("clust_i")
}
# And we merge the remaining into k groups
clust_i <- setdiff(1:k, unique(tabclust$clust_k))
tabclust$clust_k[is.na(tabclust$clust_k)] <- rep(clust_i, ceiling(nk/length(clust_i)))[1:sum(is.na(tabclust$clust_k))]
tabclust2 <- data.frame(ID = 1:length(clust_nk), clust_nk = clust_nk)
tabclust2 <- merge(tabclust2, tabclust, by = "clust_nk")
tabclust2 <- tabclust2[order(tabclust2$ID),]
clust_k <- tabclust2$clust_k
# Check size of clust_k
if(is.null(test_prop)) {
prop_valid <- !(any(table(clust_k)/length(clust_k)>maxp))
} else {
# For train/test splits, only compute W if < maxp and > minp
# Calculate the proportion by group (train/test)
prop_1 <- mean(clust_k == 1)
prop_2 <- mean(clust_k == 2)
props <- c(prop_1, prop_2)
# Keep only groups within range minp–maxp
prop_valid <- props >= minp & props <= maxp
}
if(any(prop_valid)){
if(isTRUE(dist_fun == "great_circle")){
Gjstar_i <- distclust_distmat(distmat, clust_k)
}else{
Gjstar_i <- cv_distances(tcoords, CVtest = clust_k,algorithm=algorithm, dist_fun = dist_fun)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
# Compute W statistic if not exceeding maxp
if(!any(table(clust_k)/length(clust_k)>maxp)){
if(isTRUE(dist_fun == "great_circle")){
Gjstar_i <- distclust_distmat(distmat, clust_k)
}else{
Gjstar_i <- cv_distances(tcoords, CVtest = clust_k,algorithm=algorithm, dist_fun = dist_fun)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
}
# Final configuration
k_final <- clustgrid$nk[which.min(clustgrid$W)]
W_final <- min(clustgrid$W, na.rm=T)
clust <- clustgroups[[paste0("nk", k_final)]]
if(!is.null(test_prop) && is.null(clust)) {
stop("No valid train/test configurations found in the range test_prop +/- tolerance. Increase tolerance.")
}
if(isTRUE(dist_fun == "great_circle")){
Gjstar <- distclust_distmat(distmat, clust)
}else{
Gjstar <- cv_distances(tcoords, CVtest = clust,algorithm=algorithm, dist_fun = dist_fun)
}
}
# Output
if(is.null(test_prop)) {
cfolds <- CAST::CreateSpacetimeFolds(data.frame(clust=clust), spacevar = "clust", k = k)
} else {
# Assign train/test classes
deviation_1 <- abs((table(clust)[[1]] / length(clust)) - test_prop)
deviation_2 <- abs((table(clust)[[2]] / length(clust)) - test_prop)
if(deviation_1 > deviation_2) {
test_class <- 2
} else {
test_class <- 1
}
clust[clust == test_class] <- "test"
clust[clust != "test"] <- "train"
cfolds <- list("indexOut" = which(clust == "test"), "index" = which(clust == "train"))
}
res <- list(clusters = clust,
indx_train = cfolds$index, indx_test = cfolds$indexOut,
Gij = Gij, Gj = Gj, Gjstar = Gjstar,
W = W_final, method = clustering, q = k_final, dist_space = "geographical")
class(res) <- c("knndm", "list")
res
}
# kNNDM in the feature space
knndm_feature <- function(tpoints, predpoints, k, maxp, minp, test_prop,
clustering, linkf, nk_len, dist_fun, dist_space, algorithm, catVars, scale_vars) {
# rescale data (optional)
if(isTRUE(scale_vars)) {
if(is.null(catVars)) {
scale_attr <- attributes(scale(tpoints))
tpoints <- scale(tpoints) |> as.data.frame()
predpoints <- scale(predpoints,center=scale_attr$`scaled:center`,
scale=scale_attr$`scaled:scale`) |>
as.data.frame()
} else {
tpoints_cat <- tpoints[,catVars,drop=FALSE]
predpoints_cat <- predpoints[,catVars,drop=FALSE]
tpoints_num <- tpoints[,-which(names(tpoints)%in%catVars),drop=FALSE]
predpoints_num <- predpoints[,-which(names(predpoints)%in%catVars),drop=FALSE]
scale_attr <- attributes(scale(tpoints_num))
tpoints <- scale(tpoints_num) |> as.data.frame()
predpoints <- scale(predpoints_num,center=scale_attr$`scaled:center`,
scale=scale_attr$`scaled:scale`) |>
as.data.frame()
tpoints <- as.data.frame(cbind(tpoints, lapply(tpoints_cat, as.factor)))
predpoints <- as.data.frame(cbind(predpoints, lapply(predpoints_cat, as.factor)))
}
}
# Gj and Gij calculation
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
tpoints_mat <- as.matrix(tpoints)
predpoints_mat <- as.matrix(predpoints)
# use Mahalanobis distances
if (dim(tpoints_mat)[2] == 1) {
S <- matrix(stats::var(tpoints_mat), 1, 1)
tpoints_mat <- as.matrix(tpoints_mat, ncol = 1)
} else {
S <- stats::cov(tpoints_mat)
}
S_inv <- MASS::ginv(S)
# calculate distance matrix
n_rows <- nrow(tpoints_mat)
distmat <- vapply(seq_len(n_rows), function(i) {
vapply(seq_len(n_rows), function(j) {
diff <- tpoints_mat[i, ] - tpoints_mat[j, ]
sqrt(t(diff) %*% S_inv %*% diff)
}, numeric(1))
}, numeric(n_rows))
diag(distmat) <- NA
Gj <- apply(distmat, 1, min, na.rm=TRUE)
n_rows_p <- nrow(predpoints_mat)
n_rows_t <- nrow(tpoints_mat)
Gij <- vapply(seq_len(n_rows_p), function(i) {
min(vapply(seq_len(n_rows_t), function(j) {
diff <- predpoints_mat[i, ] - tpoints_mat[j, ]
sqrt(t(diff) %*% S_inv %*% diff)
}, numeric(1)))
}, numeric(1))
} else {
# use FNN with Euclidean distances if no categorical variables are present
Gj <- c(FNN::knn.dist(tpoints, k = 1, algorithm=algorithm))
Gij <- c(FNN::knnx.dist(query = predpoints, data = tpoints, k = 1, algorithm=algorithm))
}
} else {
# use Gower distances if categorical variables are present
Gj <- vapply(1:nrow(tpoints), function(i) gower::gower_topn(tpoints[i,], tpoints[-i,], n=1)$distance[[1]], numeric(1))
Gij <- c(gower::gower_topn(predpoints, tpoints, n = 1)$distance)
}
# Check if Gj > Gij (warning suppressed regarding ties)
testks <- suppressWarnings(stats::ks.test(Gj, Gij, alternative = "great"))
if(testks$p.value >= 0.05){
if(!is.null(test_prop)) {
ntest <- floor(test_prop * nrow(tpoints))
ntrain <- nrow(tpoints) - ntest
# Create a vector: 1 = train, 2 = test
clusters <- c(rep(1, ntrain), rep(2, ntest))
# Shuffle randomly
clust <- sample(clusters, nrow(tpoints))
} else {
clust <- sample(rep(1:k, ceiling(nrow(tpoints)/k)), size = nrow(tpoints), replace=F)
}
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
} else {
Gjstar <- cv_distances(tpoints, clust,algorithm=algorithm, dist_fun = dist_fun)
}
} else {
Gjstar <- cv_distances(tpoints, clust, dist_fun = dist_fun)
}
k_final <- "random CV"
W_final <- twosamples::wass_stat(Gjstar, Gij)
message("Gij <= Gj; a random CV assignment is returned")
}else{
if(clustering == "hierarchical"){
# calculate distance matrix which is needed for hierarchical clustering
if(is.null(catVars)) {
if(isTRUE(dist_fun == "euclidean")) {
# calculate distance matrix with Euclidean distances if no categorical variables are present
# for MD: distance matrix was already calculated
distmat <- stats::dist(tpoints, upper=TRUE, diag=TRUE) |> as.matrix()
diag(distmat) <- NA
}
} else {
# calculate distance matrix with Gower distances if categorical variables are present
distmat <- matrix(nrow=nrow(tpoints), ncol=nrow(tpoints))
for (i in 1:nrow(tpoints)){
trainDist <- gower::gower_dist(tpoints[i,], tpoints)
trainDist[i] <- NA
distmat[i,] <- trainDist
}
}
hc <- stats::hclust(d = stats::as.dist(distmat), method = linkf)
}
# Build grid of number of clusters to try - we sample low numbers more intensively
clustgrid <- data.frame(nk = as.integer(round(exp(seq(log(k), log(nrow(tpoints)-2),
length.out = nk_len)))))
clustgrid$W <- NA
clustgrid <- clustgrid[!duplicated(clustgrid$nk),]
clustgroups <- list()
# Compute 1st PC for ordering clusters
if(is.null(catVars)) {
pcacoords <- stats::prcomp(tpoints, center = TRUE, scale. = FALSE, rank = 1)
} else {
pcacoords <- PCAmixdata::PCAmix(X.quanti = tpoints[,!(names(tpoints) %in% catVars), drop=FALSE],
X.quali = tpoints[,names(tpoints) %in% catVars, drop=FALSE],
graph = FALSE)
}
# We test each number of clusters
for(nk in clustgrid$nk) {
# Create nk clusters
if(clustering == "hierarchical"){
clust_nk <- stats::cutree(hc, k=nk)
} else if(clustering == "kmeans"){
if(is.null(catVars)) {
clust_nk <- tryCatch(stats::kmeans(tpoints, nk)$cluster,
error=function(e) e)
} else {
# prototype clustering for mixed data sets
clust_nk <- tryCatch(clustMixType::kproto(tpoints, nk,verbose=FALSE)$cluster,
error=function(e) e)
}
}
if (!inherits(clust_nk,"error")){
tabclust <- as.data.frame(table(clust_nk))
tabclust$clust_k <- NA
# compute cluster centroids and apply PC loadings to shuffle along the 1st dimension
if(is.null(catVars)) {
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca <- matrix(colMeans(tpoints[clust_nk %in% x, , drop = FALSE]), nrow = 1)
colnames(centrpca) <- colnames(tpoints)
return(predict(pcacoords, centrpca))
},numeric(1))
} else {
centr_tpoints <- vapply(tabclust$clust_nk, function(x){
centrpca_num <- matrix(apply(tpoints[clust_nk %in% x, !(names(tpoints) %in% catVars), drop=FALSE], 2, mean), nrow = 1)
centrpca_cat <- matrix(apply(tpoints[clust_nk %in% x, names(tpoints) %in% catVars, drop=FALSE], 2,
function(y) names(which.max(table(y)))), nrow = 1)
colnames(centrpca_num) <- colnames(tpoints[,!(names(tpoints) %in% catVars), drop=FALSE])
colnames(centrpca_cat) <- colnames(tpoints[,names(tpoints) %in% catVars, drop=FALSE])
return(predict(pcacoords, centrpca_num, centrpca_cat)[,1])
}, numeric(1))
}
tabclust$centrpca <- centr_tpoints
tabclust <- tabclust[order(tabclust$centrpca),]
# We don't merge big clusters
if(is.null(test_prop)) {
clust_i <- 1
for(i in 1:nrow(tabclust)){
if(tabclust$Freq[i] >= nrow(tpoints)/k){
tabclust$clust_k[i] <- clust_i
clust_i <- clust_i + 1
}
}
rm("clust_i")
}
# And we merge the remaining into k groups
clust_i <- setdiff(1:k, unique(tabclust$clust_k))
tabclust$clust_k[is.na(tabclust$clust_k)] <- rep(clust_i, ceiling(nk/length(clust_i)))[1:sum(is.na(tabclust$clust_k))]
tabclust2 <- data.frame(ID = 1:length(clust_nk), clust_nk = clust_nk)
tabclust2 <- merge(tabclust2, tabclust, by = "clust_nk")
tabclust2 <- tabclust2[order(tabclust2$ID),]
clust_k <- tabclust2$clust_k
# Check size of clust_k
if(is.null(test_prop)) {
prop_valid <- !(any(table(clust_k)/length(clust_k)>maxp))
} else {
# For train/test splits, only compute W if < maxp and > minp
# Calculate the proportion by group (train/test)
prop_1 <- mean(clust_k == 1)
prop_2 <- mean(clust_k == 2)
props <- c(prop_1, prop_2)
# Keep only groups within range minp–maxp
prop_valid <- props >= minp & props <= maxp
}
# Compute W statistic if size of clust_k is valid
if(any(prop_valid)){
if(clustering == "kmeans") {
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")){
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k, dist_fun = dist_fun)
} else {
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k,algorithm = algorithm, dist_fun = dist_fun)
}
} else {
Gjstar_i <- cv_distances(tpoints, CVtest = clust_k, dist_fun = dist_fun)
}
} else {
Gjstar_i <- distclust_distmat(distmat, clust_k)
}
clustgrid$W[clustgrid$nk==nk] <- twosamples::wass_stat(Gjstar_i, Gij)
clustgroups[[paste0("nk", nk)]] <- clust_k
}
}
}
# Final configuration
k_final <- clustgrid$nk[which.min(clustgrid$W)]
W_final <- min(clustgrid$W, na.rm=T)
clust <- clustgroups[[paste0("nk", k_final)]]
if(!is.null(test_prop) && is.null(clust)) {
stop("No valid train/test configurations found in the range test_prop +/- tolerance. Increase tolerance.")
}
if(clustering == "kmeans") {
if(is.null(catVars)) {
if(isTRUE(dist_fun == "mahalanobis")) {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
} else {
Gjstar <- cv_distances(tpoints, CVtest = clust,algorithm=algorithm, dist_fun = dist_fun)
}
} else {
Gjstar <- cv_distances(tpoints, CVtest = clust, dist_fun = dist_fun)
}
} else {
Gjstar <- distclust_distmat(distmat, clust)
}
}
# Output
if(is.null(test_prop)) {
cfolds <- CAST::CreateSpacetimeFolds(data.frame(clust=clust), spacevar = "clust", k = k)
} else {
# Assign train/test classes
deviation_1 <- abs((table(clust)[[1]] / length(clust)) - test_prop)
deviation_2 <- abs((table(clust)[[2]] / length(clust)) - test_prop)
if(deviation_1 > deviation_2) {
test_class <- 2
} else {
test_class <- 1
}
clust[clust == test_class] <- "test"
clust[clust != "test"] <- "train"
cfolds <- list("indexOut" = which(clust == "test"), "index" = which(clust == "train"))
}
res <- list(clusters = clust,
indx_train = cfolds$index, indx_test = cfolds$indexOut,
Gij = Gij, Gj = Gj, Gjstar = Gjstar,
W = W_final, method = clustering, q = k_final, dist_space = "feature")
class(res) <- c("knndm", "list")
res
}
# Helper function: Compute out-of-fold NN distance based on a distance matrix (geographical coordinates / numerical variables)
distclust_distmat <- function(distm, folds){
alldist <- rep(NA, length(folds))
for(f in unique(folds)){
alldist[f == folds] <- apply(distm[f == folds, f != folds, drop=FALSE], 1, min)
}
alldist
}
Try the CAST package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.