---
title: "Model Comparison with Blocked Cross-Validation"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Model Comparison with Blocked Cross-Validation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  dev = "svglite",
  fig.ext = "svg"
)
library(BORG)
library(ggplot2)

theme_borg <- theme_minimal(base_size = 12) +
  theme(
    panel.background = element_rect(fill = "transparent", colour = NA),
    plot.background = element_rect(fill = "transparent", colour = NA),
    legend.background = element_rect(fill = "transparent", colour = NA)
  )
theme_set(theme_borg)
```

# Introduction

Model comparison in ecology and environmental science usually means fitting
several candidate models to a dataset, evaluating each with cross-validation,
and picking the one with the best metric. The logic is sound, but the
implementation depends on one assumption that spatial data violates:
observations are independent.

When observations are spatially autocorrelated, random k-fold CV scatters
nearby points across training and test sets. The training set already
"knows" much of what the test set contains because neighbouring observations
leak information across the fold boundary. Every model benefits from this
leak, but not equally. A complex model with many parameters can exploit
spatial autocorrelation in the predictors more effectively than a simple one.
The result is that metric inflation varies by model complexity, which means
random CV can reverse the ranking of candidate models compared to honest,
spatially blocked CV.

This is not a theoretical concern. Roberts et al. (2017) demonstrated
inflation exceeding 50% in realistic species distribution modelling
scenarios. Wadoux et al. (2021) showed that model selection based on random
CV favours overly complex models that perform poorly at genuinely new
locations. The mechanism is straightforward: a model that memorises local
spatial patterns scores well when tested on nearby points but fails when
extrapolating to unseen regions.

The practical consequence for model selection is subtle but important. If
random CV inflates all models by the same amount, the ranking stays the same
and the modeller simply over-reports absolute performance. In practice,
however, the inflation is differential. Models with more parameters or more
flexible functional forms tend to benefit more from spatial leakage because
they have the capacity to learn short-range spatial patterns that blocked CV
would exclude from the training signal. A linear model with two predictors
might see 10% inflation. A random forest with 50 predictors might see 40%.
The modeller picks the random forest, reports excellent cross-validated
performance, and publishes. The random forest then fails at new locations
because the signal it learned was partly spatial artefact.

BORG provides a complete toolkit for fair model comparison under spatial
dependence. This vignette walks through each component: quantifying inflation,
ranking models on blocked folds, attaching confidence intervals, testing
against a null model, and diagnosing per-fold variation. We use simulated
data throughout so the ground truth is known and every result can be verified.

# Simulating spatial data

We start with `borg_simulate()`, which generates a dataset with controlled
spatial autocorrelation, a known number of informative predictors, and a
known true R-squared. Having the ground truth lets us check whether our
comparison pipeline recovers the correct signal.

```{r simulate}
sim <- borg_simulate(
  n = 500,
  type = "spatial",
  n_predictors = 5,
  signal_strength = 0.3,
  autocorrelation = 0.6,
  seed = 42
)
```

The simulated dataset contains coordinate columns `x` and `y`, five
predictors (`x1` through `x5`), and a response column. The `signal_strength`
of 0.3 means the predictors explain 30% of the variance in the response under
ideal conditions. The `autocorrelation` of 0.6 produces moderate spatial
structure in both predictors and residuals.

```{r inspect-sim}
d <- sim$data
cat("Dimensions:", nrow(d), "x", ncol(d), "\n")
cat("True R-squared:", round(sim$true_r2, 3), "\n")
cat("Coordinate columns:", paste(sim$coords, collapse = ", "), "\n")
```

The true R-squared is the upper bound for any honest evaluation. Any
cross-validated R-squared that substantially exceeds this value is a sign
of information leakage.

Before comparing models, we diagnose the spatial structure. `borg_diagnose()`
computes Moran's I, estimates the autocorrelation range, and recommends a
CV strategy.

```{r diagnose}
diag <- borg_diagnose(
  d,
  coords = sim$coords,
  target = "y",
  verbose = TRUE
)
```

With autocorrelation at 0.6, the diagnosis should flag moderate to strong
spatial dependence and recommend blocked CV. The diagnosis object contains
the estimated autocorrelation range, which determines the minimum distance
between training and test observations in the blocked CV scheme. Folds are
constructed so that no test observation is closer to any training observation
than this range, eliminating the information leak that plagues random CV.

We use this diagnosis throughout the rest of the vignette to ensure
consistent fold generation. In a real analysis, you would compute the
diagnosis once and pass it to all downstream functions to guarantee that
every comparison uses the same blocking structure.

# Quantifying inflation

The first step in any model comparison is to measure how much random CV
inflates performance estimates on your data. `borg_compare_cv()` runs the
same model under both random and blocked CV, repeats each evaluation
multiple times for stability, and reports the difference.

```{r compare-cv-rmse}
comp_rmse <- borg_compare_cv(
  data = d,
  formula = y ~ x1 + x2 + x3 + x4 + x5,
  coords = sim$coords,
  metric = "rmse",
  v = 5,
  repeats = 10,
  seed = 42,
  verbose = FALSE
)
comp_rmse
```

The output shows the mean RMSE under random and blocked CV, the inflation
percentage, and a paired t-test p-value. For error metrics like RMSE, the
inflation manifests as deflation: random CV reports a lower (apparently
better) error than the honest blocked estimate.

We can repeat the comparison with R-squared to see inflation from the other
direction.

```{r compare-cv-rsq}
comp_rsq <- borg_compare_cv(
  data = d,
  formula = y ~ x1 + x2 + x3 + x4 + x5,
  coords = sim$coords,
  metric = "rsq",
  v = 5,
  repeats = 10,
  seed = 42,
  verbose = FALSE
)
comp_rsq
```

For R-squared, random CV inflates the estimate directly. The blocked CV
R-squared should sit closer to (and often below) the true R-squared of
`r round(sim$true_r2, 3)`. The random CV R-squared will typically exceed
the true value, sometimes substantially.

The `plot()` method for `borg_comparison` objects supports three
visualisation types. The boxplot is the clearest for showing separation
between strategies.

```{r plot-boxplot, fig.cap = "Distribution of RMSE across 10 repeats of random vs blocked CV."}
plot(comp_rmse, type = "boxplot")
```

The density view shows the full distribution of estimates across repeats.
Well-separated densities indicate that the choice of CV strategy materially
affects the reported metric.

```{r plot-density, fig.cap = "Density of R-squared estimates under random and blocked CV."}
plot(comp_rsq, type = "density")
```

The paired plot links each repeat, making it easy to see whether the
inflation is consistent or driven by a few extreme repeats.

```{r plot-paired, fig.cap = "Paired comparison: each point is one repeat of the full CV."}
plot(comp_rsq, type = "paired")
```

The inflation percentage is the key quantity for reporting. If blocked CV
gives RMSE = 1.2 and random CV gives RMSE = 1.0, the inflation (deflation
for error metrics) is (1.2 - 1.0) / 1.2 = 16.7%. Values below 5% suggest
that spatial structure is weak enough to tolerate random CV. Values above
15% indicate that random CV produces misleading estimates.

The p-value from the paired t-test tells us whether the difference between
strategies is statistically significant. Because each repeat generates a
pair of estimates (one random, one blocked) on the same data, the pairing
controls for repeat-to-repeat variation in fold assignment. A significant
result (p < 0.05) is standard evidence that random CV is not appropriate
for this dataset. An insignificant result does not prove that random CV is
safe, only that the test lacked power to detect the difference at the
current number of repeats.

The magnitude of inflation depends on three factors: the strength of
autocorrelation, the density of observations relative to the autocorrelation
range, and the complexity of the model. Strongly autocorrelated data with
dense sampling and complex models produce the largest inflation. Weakly
autocorrelated data with sparse sampling and simple models produce the
smallest. Our simulated data with autocorrelation 0.6 and 500 observations
is a moderate scenario. Real-world ecological datasets often have stronger
autocorrelation and denser sampling, so the inflation we see here is
conservative.

# Fair model ranking

Knowing that inflation exists, the next step is to compare candidate models
on a level playing field. `borg_compare_models()` evaluates multiple
formulas on the same set of blocked folds, ensuring that differences in
metric reflect genuine predictive ability rather than differential leakage.

We define three candidate models of increasing complexity: a simple model
with two predictors, a medium model with three, and the full model with
all five.

```{r define-models}
models <- list(
  simple = y ~ x1 + x2,
  medium = y ~ x1 + x2 + x3,
  full   = y ~ x1 + x2 + x3 + x4 + x5
)
```

We generate blocked folds once and pass them to `borg_compare_models()`.
Using the same folds for all models eliminates fold assignment as a source
of variation.

```{r blocked-folds}
cv_blocked <- borg_cv(
  d,
  coords = sim$coords,
  target = "y",
  v = 5,
  verbose = FALSE
)
```

```{r compare-models-blocked}
rank_blocked <- borg_compare_models(
  data = d,
  folds = cv_blocked,
  models = models,
  metric = "rmse"
)
rank_blocked
```

Now we repeat the comparison with random folds to show how rankings can
shift.

```{r random-folds}
n <- nrow(d)
random_fold_ids <- sample(rep(1:5, length.out = n))
cv_random <- lapply(1:5, function(k) {
  list(
    train = which(random_fold_ids != k),
    test = which(random_fold_ids == k)
  )
})
```

```{r compare-models-random}
rank_random <- borg_compare_models(
  data = d,
  folds = cv_random,
  models = models,
  metric = "rmse"
)
rank_random
```

When the inflation is differential, rankings can reverse. A complex model
that overfits spatial structure may rank first under random CV but last
under blocked CV. The blocked ranking is the one that predicts how models
will perform at new locations outside the study area.

The `autoplot()` method produces a bar chart with error bars for quick
visual comparison.

```{r autoplot-models, fig.cap = "Model ranking under blocked spatial CV. Error bars show +/- 1 SD across folds."}
autoplot(rank_blocked)
```

When the best model's error bar overlaps with the second-best model's error
bar, the two are statistically indistinguishable. In that case, prefer the
simpler model. More parameters mean more opportunities for overfitting to
spatial structure that will not be present at prediction locations.

The comparison between random and blocked rankings is itself informative.
If both strategies produce the same ranking, the relative ordering of
models is robust to the choice of CV strategy even if the absolute metrics
differ. If the rankings diverge, the models that benefit most from random CV
are the ones most likely to disappoint at new locations. We can formalise
this by computing the rank correlation (Kendall's tau or Spearman's rho)
between the two rankings. A rank correlation below 0.7 across more than
three candidate models is a strong signal that random CV is producing
misleading model selection.

One subtlety: `borg_compare_models()` evaluates each model on the same
folds, but it refits the model on each fold's training set. This means
the comparison is between the models' generalisable predictive performance,
not between their fits to the full dataset. A model that fits the full
dataset well but generalises poorly will rank low. This is the correct
behaviour for prediction-oriented model comparison, but it differs from
AIC or BIC comparisons, which evaluate the fit to the full dataset without
a train-test split.

# Confidence intervals on metrics

A point estimate of RMSE or R-squared is incomplete without a measure of
uncertainty. Standard bootstrap assumes independent observations, which
underestimates variance when data are spatially correlated. `borg_bootstrap()`
uses spatial block bootstrap: it partitions the study area into spatial
blocks, resamples blocks with replacement, and refits the model on each
bootstrap sample.

We compute block bootstrap confidence intervals for the full model.

```{r bootstrap-full}
fit_full <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = d)

boot_full <- borg_bootstrap(
  model = fit_full,
  data = d,
  target = "y",
  coords = sim$coords,
  metric = "rmse",
  n_boot = 200,
  n_blocks = 10,
  conf_level = 0.95,
  seed = 42
)
boot_full
```

The output reports the point estimate, 95% confidence interval, bootstrap
standard error, and bias. The CI width reflects the true uncertainty in the
metric after accounting for spatial dependence. A standard (non-block)
bootstrap on the same data would produce narrower intervals that understate
the uncertainty.

We repeat for the simple model.

```{r bootstrap-simple}
fit_simple <- lm(y ~ x1 + x2, data = d)

boot_simple <- borg_bootstrap(
  model = fit_simple,
  data = d,
  target = "y",
  coords = sim$coords,
  metric = "rmse",
  n_boot = 200,
  n_blocks = 10,
  conf_level = 0.95,
  seed = 42
)
boot_simple
```

When comparing two models, the question is whether their confidence
intervals overlap. Non-overlapping intervals provide strong evidence that
one model genuinely outperforms the other. Overlapping intervals mean the
difference is within the range of sampling variability, and neither model
can be declared superior with confidence.

```{r autoplot-bootstrap, fig.cap = "Block bootstrap distribution of RMSE for the full model. Dashed lines mark the 95% CI."}
autoplot(boot_full)
```

The `autoplot()` method shows the bootstrap distribution with the point
estimate and CI bounds. A symmetric, roughly normal distribution indicates
that the bootstrap is well-behaved. Heavy skewness or multimodality suggests
that the metric is unstable across spatial blocks, which is itself a
diagnostic signal.

A useful heuristic: if the bootstrap SE exceeds 20% of the point estimate,
the metric is too noisy to support fine-grained model comparisons. In that
case, focus on whether models are clearly separated rather than on small
differences in rank.

The number of spatial blocks (`n_blocks`) controls the trade-off between
preserving spatial structure and having enough resampling units for stable
estimates. Too few blocks (e.g., 3) means each block is large and spatially
heterogeneous, which weakens the block bootstrap's ability to capture
within-block correlation. Too many blocks (e.g., 50) means each block
contains few observations, making the resampled datasets unrepresentative.
Ten blocks is a reasonable default for datasets with 200-1000 observations.
For larger datasets, 15-20 blocks work well.

The bias estimate deserves attention. Block bootstrap bias measures the
systematic difference between the bootstrap estimate and the point estimate.
A large bias (more than 10% of the point estimate) suggests that the
spatial structure interacts with the metric in a way that the point estimate
does not capture. This happens most often with R-squared on small datasets,
where the point estimate can be optimistically biased. In that case, report
the bias-corrected estimate alongside the CI.

# Is the signal real?

Before interpreting model comparisons, we should verify that the best model
captures a genuine signal rather than fitting noise. `borg_null_test()` runs
a permutation test: it shuffles the response variable, re-evaluates the
model under the same blocked CV scheme, and builds a null distribution of
metric values. If the observed metric falls far enough from the null
distribution, the model has detected real structure.

```{r null-test}
null_full <- borg_null_test(
  data = d,
  folds = cv_blocked,
  formula = y ~ x1 + x2 + x3 + x4 + x5,
  metric = "rmse",
  n_null = 99,
  seed = 42
)
null_full
```

The output reports the empirical metric, the null mean and standard
deviation, a z-score, and a p-value. For RMSE (lower is better), a
significant result means the observed RMSE is substantially below the null
distribution. With a signal strength of 0.3 and 500 observations, we
expect a clear rejection.

```{r autoplot-null, fig.cap = "Null distribution of RMSE from 99 permutations. The red line marks the observed value."}
autoplot(null_full)
```

The `autoplot()` method overlays the empirical value on the null
distribution. When the observed value sits in the tail, the model is
capturing real structure. When it sits within the bulk of the null
distribution, the model has no predictive skill beyond chance.

We can also test the simple model to see whether two predictors are
sufficient to detect the signal.

```{r null-test-simple}
null_simple <- borg_null_test(
  data = d,
  folds = cv_blocked,
  formula = y ~ x1 + x2,
  metric = "rmse",
  n_null = 99,
  seed = 42
)
null_simple
```

If the simple model also rejects the null, both models detect real signal.
The model comparison then becomes about which captures more of it. If the
simple model fails to reject, the additional predictors in the full model
are necessary for detecting the signal at all.

A note on the number of permutations: 99 gives a minimum possible p-value
of 0.01, which suffices for exploratory analysis. For publication, use
999 permutations (minimum p = 0.001). Going beyond 999 rarely changes
conclusions but increases computation time linearly.

The permutation test is conservative in the presence of spatial
autocorrelation. Shuffling the response variable destroys the spatial
structure of the response but preserves the spatial structure of the
predictors. If the predictors themselves are spatially autocorrelated (which
they usually are), the null model still captures some spatial pattern through
the predictors alone. This means the null distribution is shifted towards
better performance compared to a truly skill-less model, making it harder
to reject the null. A significant result under these conditions is strong
evidence of genuine predictive skill.

The z-score provides a continuous measure of how far the observed metric
deviates from the null. A z-score of -3 (for RMSE, where lower is better)
means the observed error is 3 standard deviations below the null mean. In
ecological modelling, z-scores beyond +/- 2 are typical for models with
real signal, and z-scores beyond +/- 4 indicate strong effects. A z-score
near zero means the model is no better than random assignment of the
response to locations.

# Per-fold performance

Aggregate metrics hide spatial variation in model quality. A model might
perform well on average but fail badly in a particular region of the study
area. `borg_fold_performance()` returns the metric for each fold separately,
along with the spatial centroid of each test set.

```{r fold-perf}
fold_perf <- borg_fold_performance(
  data = d,
  folds = cv_blocked,
  formula = y ~ x1 + x2 + x3 + x4 + x5,
  coords = sim$coords,
  metric = "rmse"
)
fold_perf
```

The output is a data frame with one row per fold, showing the metric value,
training and test set sizes, and the spatial centroid of the test set. Large
variation across folds indicates that the model performs unevenly across the
study area. This can happen when environmental gradients are steeper in some
regions, when sampling density is uneven, or when the model's functional
form is a poor match for local patterns.

```{r autoplot-folds, fig.cap = "Per-fold RMSE. Variation across folds reflects spatial heterogeneity in model performance."}
autoplot(fold_perf)
```

The `autoplot()` method highlights folds with unusually high or low metric
values. Folds with high RMSE correspond to regions where the model
struggles. If coordinate centroids are available, these can be mapped to
identify the problematic areas.

```{r fold-variation}
fold_range <- range(fold_perf$value)
fold_cv <- sd(fold_perf$value) / mean(fold_perf$value)
cat("RMSE range across folds:", round(fold_range[1], 3),
    "to", round(fold_range[2], 3), "\n")
cat("Coefficient of variation:", round(fold_cv * 100, 1), "%\n")
```

A coefficient of variation below 15% suggests fairly uniform performance.
Above 30%, the aggregate metric is a poor summary and the modeller should
investigate whether a different model specification or additional predictors
could reduce the spatial heterogeneity.

When one fold is substantially worse than the others, check whether its
test region is at the edge of the study area (extrapolation), has lower
sampling density, or lacks coverage of a key environmental gradient. These
diagnostics often reveal more about the modelling problem than the aggregate
metric alone.

Per-fold analysis also helps with model selection. Compare the fold
performance profiles of two candidate models. If model A beats model B on
four out of five folds but loses badly on the fifth, model A has a
vulnerability in a specific spatial region. Model B might be the safer choice
if that region matters for the application. The aggregate metric would
favour model A, but the per-fold breakdown reveals the risk.

Spatial centroids in the output allow geographic mapping of performance.
Plotting metric values at centroid locations produces a performance map
that shows where predictions are reliable and where they are not. This is
especially valuable for species distribution models, where conservation
decisions depend on the spatial pattern of prediction quality, not just
the average.

# Full comparison pipeline

The individual tools combine into a complete model comparison workflow. Here
we run the full pipeline from simulation to final table, using all five
components.

**Step 1: Generate data and diagnose structure.**

```{r pipeline-sim}
sim2 <- borg_simulate(
  n = 400,
  type = "spatial",
  n_predictors = 5,
  signal_strength = 0.25,
  autocorrelation = 0.7,
  seed = 99
)
d2 <- sim2$data
```

```{r pipeline-diagnose}
diag2 <- borg_diagnose(
  d2,
  coords = sim2$coords,
  target = "y",
  verbose = FALSE
)
```

**Step 2: Quantify inflation to justify blocked CV.**

```{r pipeline-inflation}
comp2 <- borg_compare_cv(
  data = d2,
  formula = y ~ x1 + x2 + x3 + x4 + x5,
  coords = sim2$coords,
  metric = "rsq",
  v = 5,
  repeats = 10,
  seed = 99,
  verbose = FALSE
)
cat("Inflation:", round(abs(comp2$inflation$estimate) * 100, 1), "%\n")
cat("P-value:", format(comp2$p_value, digits = 3), "\n")
```

With autocorrelation at 0.7, the inflation should be substantial, providing
clear justification for using blocked CV in all subsequent comparisons.

**Step 3: Generate blocked folds and rank candidate models.**

```{r pipeline-folds}
cv2 <- borg_cv(
  d2,
  coords = sim2$coords,
  target = "y",
  v = 5,
  verbose = FALSE
)
```

```{r pipeline-rank}
models2 <- list(
  intercept = y ~ 1,
  simple    = y ~ x1 + x2,
  medium    = y ~ x1 + x2 + x3,
  full      = y ~ x1 + x2 + x3 + x4 + x5
)

rank2 <- borg_compare_models(
  data = d2,
  folds = cv2,
  models = models2,
  metric = "rmse"
)
rank2
```

The intercept-only model serves as a baseline. If the best candidate model
barely beats the intercept, the signal is too weak to support reliable
predictions at new locations.

**Step 4: Test the top model against the null.**

```{r pipeline-null}
top_formula <- models2[[rank2$model[1]]]

null2 <- borg_null_test(
  data = d2,
  folds = cv2,
  formula = top_formula,
  metric = "rmse",
  n_null = 99,
  seed = 99
)
null2
```

**Step 5: Attach confidence intervals to the top model.**

```{r pipeline-bootstrap}
fit2 <- lm(top_formula, data = d2)

boot2 <- borg_bootstrap(
  model = fit2,
  data = d2,
  target = "y",
  coords = sim2$coords,
  metric = "rmse",
  n_boot = 200,
  n_blocks = 10,
  seed = 99
)
boot2
```

**Step 6: Check per-fold variation.**

```{r pipeline-folds-perf}
fold2 <- borg_fold_performance(
  data = d2,
  folds = cv2,
  formula = top_formula,
  coords = sim2$coords,
  metric = "rmse"
)
fold2
```

The full pipeline produces five pieces of evidence:

1. The inflation test justifies the use of blocked CV.
2. The model ranking identifies the best candidate on fair folds.
3. The null test confirms that the best model captures real signal.
4. The bootstrap CI quantifies uncertainty in the reported metric.
5. The per-fold breakdown reveals spatial heterogeneity in performance.

Together, these form a complete reporting package. A reviewer can verify
that the comparison accounted for spatial dependence, that the best model
genuinely outperforms alternatives, that the signal is real, that the
reported uncertainty is honest, and that performance is reasonably uniform
across the study area.

The order matters. Running the inflation test first establishes whether
blocked CV is necessary. If it is not (inflation < 5%), the remaining steps
can use random CV without loss of validity. Running the model ranking before
the null test ensures that we test the best candidate, not an arbitrary
one. Running the bootstrap after the ranking avoids computing expensive CIs
for models that will be discarded. And running the per-fold analysis last
provides a final check on the selected model's spatial reliability.

In practice, most ecological datasets with coordinates will show significant
inflation. The pipeline above takes a few minutes on a dataset with 500
observations and five predictors. For larger datasets or more complex models
(random forests, gradient boosting), computation time scales with the number
of repeats and permutations, so start with the defaults (10 repeats, 99
permutations, 200 bootstrap replicates) and increase for publication.

# Practical guidance

## How many repeats for `borg_compare_cv()`

Ten repeats (the default) are sufficient for screening. For publication,
use 20-30 repeats. Beyond 30, the paired t-test p-value stabilises and
additional repeats add computation without changing conclusions. Each
repeat generates new fold assignments, so the variation across repeats
reflects sensitivity to the specific partitioning of space.

## How many permutations for `borg_null_test()`

Use 99 for exploratory work, 999 for publication. The minimum achievable
p-value is 1/(n_null + 1), so 99 permutations can detect p < 0.01 and
999 permutations can detect p < 0.001. If computation is cheap (linear
models, small datasets), default to 999.

## How many bootstrap replicates

Two hundred replicates (the default) give stable 95% CIs for most metrics.
If the CI bounds change by more than 5% when you increase to 500, the
metric is inherently unstable and you should interpret it cautiously. For
99% CIs, use at least 500 replicates.

## When to report which metric

RMSE is the default choice for regression: it is in the units of the
response and penalises large errors. MAE is more robust to outliers and
appropriate when large errors are no worse than moderate ones. R-squared
is useful for comparing across datasets with different response scales
but is sensitive to the variance of the response and can be misleading
when the signal is weak.

For model comparison, use the same metric throughout. Switching metrics
between the inflation test, the ranking, and the bootstrap CI makes the
results incoherent.

## When model comparison is misleading

Model comparison assumes that all candidates are trained and evaluated on
the same data under the same conditions. This breaks down in several
situations.

First, when models require different preprocessing. A model that needs
standardised predictors must have the standardisation computed inside each
fold, using only training data. If standardisation is done before splitting,
the test set statistics leak into training, and the comparison is unfair.
BORG handles fold generation but not preprocessing, so the user must ensure
that any feature transformations respect the fold boundaries.

Second, when models have fundamentally different response types. Comparing
a regression model to a classification model requires choosing a metric
that makes sense for both, which usually means the comparison is not
meaningful. Stick to models that solve the same task.

Third, when the dataset is too small for the number of folds. With 5-fold
CV and 100 observations, each test fold has 20 points. Metrics computed on
20 points have high variance, making fine-grained comparisons unreliable.
Reducing to 3-fold CV helps but also reduces the training set size per fold,
which can hurt complex models disproportionately.

Fourth, when the candidate models are not nested. Comparing `y ~ x1 + x2`
to `y ~ x3 + x4` is valid in principle, but the interpretation changes. With
nested models, the question is whether additional predictors improve
prediction. With non-nested models, the question is which set of predictors
is more useful, and the answer may depend on which spatial region we care
about.

Finally, when the study area will change at prediction time. Blocked CV
estimates performance at new locations within the study area. If the model
will be applied to a different region with different environmental
conditions, even blocked CV is optimistic. In that case, consider leave-
region-out CV or spatial transferability analysis (see the `borg_transferability()`
function).

## Reporting for publications

A methods section should include:

1. The CV strategy used (e.g., "5-fold spatial block CV") and the
   justification (e.g., "Moran's I = 0.45, p < 0.001; inflation test
   showed 23% R-squared inflation under random CV").
2. The number of repeats and the metric.
3. The model ranking table with mean and SD across folds.
4. A null test p-value for the selected model.
5. Block bootstrap 95% CI for the reported metric.

This gives reviewers everything they need to assess whether the model
comparison was fair and the reported performance is honest.

A results section should present the model ranking table, the best model's
RMSE or R-squared with its 95% block bootstrap CI, the null test p-value,
and a statement about per-fold variation (e.g., "Per-fold RMSE ranged from
0.85 to 1.14, with a coefficient of variation of 12%, indicating
approximately uniform performance across the study area"). Include the
inflation estimate in the methods section as justification for the CV
strategy, not in the results.

## Number of folds

Five folds is the standard choice and works well for datasets with 200+
observations. Three folds are appropriate for small datasets (50-200
observations) where each fold needs a reasonable number of test points. Ten
folds increase variance across folds because each test set is smaller, but
they reduce bias because the training set is larger (90% vs 80% of the
data). For spatial data, the number of folds interacts with the block size:
more folds require more spatial blocks, which may not be achievable if the
study area is small relative to the autocorrelation range.

## Minimum sample size

Blocked CV reduces effective sample size because spatially correlated
observations within a block are not independent. Use `borg_power()` to
estimate the effective sample size and check whether the dataset supports
the comparison.

```{r power}
pwr <- borg_power(
  data = d,
  coords = sim$coords,
  target = "y",
  alpha = 0.05,
  power = 0.80
)
pwr
```

If the power analysis indicates that the effective sample size is too small,
consider reducing the number of folds (from 5 to 3), using a simpler model
with fewer parameters, or collecting additional data in undersampled
regions of the study area.

# References

Lahiri, S. N. (2003). *Resampling Methods for Dependent Data*. Springer.

Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G.,
Hauenstein, S., Lahoz-Monfort, J. J., Schroder, B., Thuiller, W., Warton, D. I.,
Wintle, B. A., Hartig, F., & Dormann, C. F. (2017). Cross-validation strategies
for data with temporal, spatial, hierarchical, or phylogenetic structure.
*Ecography*, 40(8), 913-929. doi:10.1111/ecog.02881

Wadoux, A. M. J.-C., Heuvelink, G. B. M., de Bruin, S., & Brus, D. J. (2021).
Spatial cross-validation is not the right way to evaluate map accuracy.
*Ecological Modelling*, 457, 109692. doi:10.1016/j.ecolmodel.2021.109692
