knndm is a k-fold version of NNDM LOO CV for medium and large datasets. Brielfy, the algorithm tries to
find a k-fold 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 during CV (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.
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 `ppoints` 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 a sf polygon that defines the prediction area. The function takes a regular point sample
(amount defined by `samplesize`) from the spatial extent. As an alternative use `ppoints` instead of
`modeldomain`, if you have already defined the prediction locations (e.g. raster pixel centroids).
When using either `modeldomain` or `ppoints`, we advise to plot the study area polygon and the
training/prediction points as a previous step to ensure they are aligned.