---
title: "Post-Fit Diagnostics"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Post-Fit Diagnostics}
  %\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)
```

# Setup: simulating spatial data and fitting a model

Before running any diagnostics, we need a dataset with known spatial structure and
a fitted model we can interrogate. The `borg_simulate()` function generates
synthetic datasets with controlled dependencies, which makes it possible to know
in advance how much metric inflation to expect. Throughout this vignette we use a
spatial simulation with moderate autocorrelation and moderate signal strength. The
simulation places 300 observations on a spatial grid, generates five autocorrelated
predictors, and constructs a response that depends linearly on those predictors
plus spatially correlated noise.

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

d <- sim$data
str(d)
```

Because the coordinate columns are named `x` and `y` and the response variable is
also named `y`, the resulting data frame renames the response to `y.1`. This is a
deliberate design choice in `borg_simulate()` to avoid column name collisions. We
fit a linear model using only the five predictors, then generate blocked
cross-validation folds that respect the spatial dependency structure. These folds
will be reused across most of the diagnostics that follow.

```{r fit-model}
model <- lm(y.1 ~ x1 + x2 + x3 + x4 + x5, data = d)

# Blocked spatial CV folds
cv <- borg_cv(d, coords = c("x", "y"), target = "y.1", v = 5)
```

The true R-squared for this simulation is `r round(sim$true_r2, 3)`. Any metric
estimate that substantially exceeds this value is evidence of inflation caused by
leakage through the cross-validation scheme.

# Quick diagnostic check

The fastest way to audit a dataset for leakage risks is `borg_check()`. This
single function runs the full BORG pipeline internally: it diagnoses dependency
structure, generates blocked CV folds, validates the split, and optionally
inspects a fitted model. The result is a tidy data frame where each row represents
one detected risk, with columns for the risk type, severity level, a plain-language
description, the number of affected observations, and the source object that
triggered the detection.

```{r quick-check}
chk <- borg_check(d, target = "y.1", coords = c("x", "y"))
chk
```

The severity column classifies each risk into one of three categories:
`"hard_inflation"` indicates a problem that will inflate metrics substantially,
such as test observations appearing in the training set or a feature that directly
encodes the response. `"soft_inflation"` flags a subtler issue like spatial
autocorrelation that leaks through random CV folds but may not cause dramatic
inflation in every case. `"info"` marks observations that the system flagged for
awareness but that do not necessarily compromise evaluation.

Zero rows means all tests passed.
No guarantee of zero leakage, but the check covers everything BORG can detect:
spatial and temporal autocorrelation, clustered observations, index overlap between
train and test sets, target leakage via proxy features, and preprocessing leakage
from shared normalisation statistics.

You can also pass a fitted model to `borg_check()` to run model-level inspections
in the same call. This checks whether the model object carries preprocessing
artefacts or formula terms that reference the target through a proxy variable.

```{r check-with-model}
chk_model <- borg_check(d, model = model, target = "y.1", coords = c("x", "y"))
chk_model
```

The attributes of the returned object contain the full diagnosis, the generated CV
folds, and the raw risk list. You can extract these for more detailed inspection
if the summary table raises questions.

```{r check-attributes}
# Access the underlying BorgDiagnosis
diag <- attr(chk, "diagnosis")
diag@dependency_type
diag@severity
diag@recommended_cv
```

For non-stopping diagnostics on a model object alone (without running the full
pipeline), use `borg_inspect()`. Unlike `borg()`, which throws errors on hard
violations, `borg_inspect()` returns a `BorgRisk` object that you can examine at
your own pace without interrupting a script.

# Residual diagnostics

A model might pass all pre-fit leakage checks and still produce biased predictions
if it fails to capture the spatial process in the data. `borg_check_residuals()`
tests whether the residuals from a fitted model retain spatial autocorrelation.
If they do, the model has left spatial signal on the table, and predictions at
new locations will be systematically biased toward the values of their neighbours
in the training set.

The function accepts either a fitted model object (anything with a `residuals()`
method) or a raw numeric vector of residuals. When given a model, it extracts
residuals automatically. It then computes Moran's I on those residuals using an
inverse-distance weight matrix derived from the spatial coordinates. The p-value
comes from a randomization test that permutes the residuals across locations.

```{r residual-check}
res_check <- borg_check_residuals(model, d, coords = c("x", "y"))
res_check
```

Three quantities are reported. First, Moran's I itself: a value near zero
indicates no spatial pattern in the residuals, positive values indicate clustering
of similar residuals, and negative values (rare in practice) indicate spatial
repulsion. Second, the p-value from the randomization test. Third, an overall
assessment that translates the test result into a plain-language verdict: `"clean"`
means the residuals show no statistically significant spatial structure,
`"mild"` means Moran's I is significant but below 0.3, and `"strong"` means
Moran's I exceeds 0.3, indicating that the model is missing a spatial effect of
practical importance.

The 0.3 threshold is conservative. Even a Moran's I of 0.15 can bias predictions
in interpolation-heavy workflows, so the threshold exists to avoid alarm fatigue
when the residual signal is too weak to act on, not to certify that the spatial
structure is absent or harmless.

The `autoplot()` method for the residual check object produces a variogram plot. A
variogram shows semivariance (half the mean squared difference between residual
values) as a function of the distance between observation pairs. If the model has
fully captured the spatial process, the variogram should be flat: the semivariance
should not increase with lag distance. A variogram that rises steeply from the
origin and then levels off at a sill indicates residual spatial structure with a
characteristic range.

```{r residual-variogram, fig.height = 4}
autoplot(res_check)
```

The dashed horizontal line marks the sill (maximum semivariance). Point sizes
encode the number of pairs at each distance bin, so you can judge how much data
supports each point on the curve. Thin bins at long lags are noisy and should not
drive interpretation.

When the assessment is `"mild"` or `"strong"`, adding spatial covariates
(coordinates, distance to features, a spatial smooth) to the model formula is
usually the first corrective step. After refitting, run `borg_check_residuals()`
again to confirm the spatial signal has been absorbed.

# Calibration

Performance metrics like RMSE and R-squared tell you how accurate a model is on
average, but they say nothing about whether the model's uncertainty estimates are
trustworthy. A model might report narrow prediction intervals that are right on
average but wrong for specific subsets of the input space. Calibration diagnostics
answer the question: when the model says it is 80% confident, is it actually
correct 80% of the time?

`borg_calibration()` supports both classification and regression problems. For
classification, it bins predicted probabilities and compares each bin's mean
predicted probability against the observed event frequency. For regression, it
checks whether the residual distribution matches the assumed noise model by
computing observed coverage at multiple nominal levels.

```{r calibration-regression}
preds <- predict(model, newdata = d)
cal <- borg_calibration(preds, d$y.1, type = "regression", n_bins = 10)
cal
```

Several metrics appear in the result. The Expected Calibration Error (ECE) is the
weighted average of the absolute difference between predicted and observed
frequencies across bins. The Maximum Calibration Error (MCE) captures the
worst-case bin. The reliability slope and intercept come from a weighted regression
of observed on predicted bin values: a perfectly calibrated model has slope 1 and
intercept 0. A slope below 1 means the model is overconfident (its predictions are
too spread out), while a slope above 1 means the model is underconfident.

Visualize these metrics with `autoplot()`.

```{r reliability-diagram, fig.height = 4.5}
autoplot(cal, type = "reliability")
```

The diagonal dashed line represents perfect calibration. Points above the line
indicate the model underestimates the true values in that bin; points below
indicate overestimation. Point sizes reflect the number of observations per bin,
and the shaded ribbon shows approximate confidence intervals.

For regression models, `borg_calibration()` also computes a coverage curve. This
shows, for each nominal coverage level (10%, 20%, ..., 90%), how many observations
actually fall within the corresponding prediction interval.

```{r coverage-plot, fig.height = 4.5}
autoplot(cal, type = "coverage")
```

Coverage curves that track the diagonal closely indicate well-calibrated
estimates. Curves consistently below the diagonal indicate
overconfidence: the model's intervals are too narrow and miss more observations
than they should.

Binning matters. `borg_calibration()` offers three options via the
`strategy` argument. `"uniform"` creates equal-width bins across the prediction
range, which works well when predictions are roughly uniformly distributed.
`"quantile"` creates equal-count bins, ensuring each bin has enough observations
for a reliable estimate. `"spatial"` uses k-means on the coordinates to define
bins, preventing calibration curves from being inflated by spatial autocorrelation
among predictions within the same spatial neighbourhood.

```{r calibration-spatial}
# Spatial-aware binning
coord_mat <- d[, c("x", "y")]
cal_spatial <- borg_calibration(
  preds, d$y.1,
  type = "regression",
  n_bins = 8,
  coords = coord_mat,
  strategy = "spatial"
)
cat("ECE (uniform):", round(cal$ece, 4), "\n")
cat("ECE (spatial):", round(cal_spatial$ece, 4), "\n")
```

When the spatial ECE substantially exceeds the uniform ECE, the calibration
curve was being flattered by autocorrelation: nearby predictions are similar,
so same-bin groupings artificially reduce the apparent calibration error.

# Null hypothesis testing

A low RMSE or high R-squared does not, by itself, prove that the model has learned
a real signal. If the predictors happen to be spatially structured and the response
is also spatially structured, blocked CV can still produce decent-looking metrics
even when the predictors carry no information about the response. This is
especially common in ecology and environmental science, where "everything is
correlated with everything" through shared spatial gradients.

`borg_null_test()` constructs a null distribution by permuting the response
variable and re-evaluating the model under the same blocked CV scheme. Each
permutation breaks the link between predictors and response while preserving the
spatial structure of the predictors. Should the empirical metric fall far into the
tail of the null distribution, the signal is real; if it sits near the centre, the
model's apparent performance could be explained by spatial confounding alone.

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

Reported values include the empirical metric, the null distribution's mean and
standard deviation, a z-score, and a one-sided p-value. The p-value is computed as
the proportion of null permutations that produce a metric as extreme as (or more
extreme than) the observed metric, with a finite-sample correction that adds one
to both numerator and denominator.

A significant result (p < 0.05) means the model performs better than random
expectation given the CV scheme. A non-significant result does not necessarily mean
the predictors are useless; it may mean the sample size is too small or the blocked
CV is too aggressive for the effect size. In that case, the power analysis section
below will help you determine whether the experiment was adequately powered.

`autoplot()` produces a histogram of the null distribution with the empirical
value marked as a vertical red line.

```{r null-plot, fig.height = 4}
autoplot(nt)
```

A useful convention is to report the null test alongside the primary metrics in any
manuscript or report. Reviewers who see "R-squared = 0.35, permutation p = 0.01
(99 permutations)" immediately know the signal survived a spatial confounding
check. The number of permutations controls the resolution of the p-value: 99
permutations give a minimum p of 0.01, 999 give 0.001. For publication, 999
permutations is standard.

Be aware that the null test assumes the model-fitting function is deterministic
(or at least produces stable results across permutations). For stochastic algorithms like random forests, set the seed inside `fit_fun` or use a
wrapper that fixes internal randomness. Otherwise the null distribution will mix
model variance with permutation variance, inflating the apparent significance.

# Model comparison

Once you know the signal is real, the next question is which model captures it
best. BORG provides two functions for model comparison: `borg_compare_cv()`, which
compares random versus blocked cross-validation to quantify metric inflation, and
`borg_compare_models()`, which evaluates multiple models under the same blocked CV
folds for fair side-by-side ranking.

## Quantifying inflation with `borg_compare_cv()`

`borg_compare_cv()` runs both random and blocked cross-validation on the same data
and model, repeating each scheme multiple times for stable estimates. It then
computes the inflation: the percentage by which random CV underestimates the error
(for error metrics like RMSE) or overestimates the performance (for metrics like
R-squared). A paired t-test assesses whether the difference is statistically
significant.

```{r compare-cv}
comp <- borg_compare_cv(
  data    = d,
  formula = y.1 ~ x1 + x2 + x3 + x4 + x5,
  coords  = c("x", "y"),
  metric  = "rmse",
  repeats = 10,
  seed    = 123,
  verbose = FALSE
)
comp
```

The inflation estimate tells you how much your reported metrics would change if
you switched from random to blocked CV. Put concretely, an inflation of 15% on
RMSE means random CV underestimates the true prediction error by 15%, which is
more than enough to mislead model selection decisions, overstate the model's value
in a manuscript, and give reviewers a false sense of how well the predictions will
generalise to new spatial locations outside the training domain.

The base R `plot()` method offers three views of the comparison.

```{r compare-cv-plot, fig.height = 4.5}
plot(comp, type = "boxplot")
```

The boxplot shows the distribution of metric values across repeats for each CV
scheme. Non-overlapping boxes indicate inflation that is both statistically
significant and practically meaningful. The density view (`type = "density"`)
overlays the two distributions for a more direct visual comparison, and the paired
view (`type = "paired"`) shows each repeat as a point in random-vs-blocked space,
with perfect agreement on the diagonal.

## Ranking models with `borg_compare_models()`

When you have several candidate models, `borg_compare_models()` evaluates all of
them using the same blocked CV folds. This ensures differences in metrics reflect
genuine predictive differences rather than variation in fold assignments.

```{r compare-models}
model_list <- list(
  simple   = y.1 ~ x1,
  additive = y.1 ~ x1 + x2 + x3,
  full     = y.1 ~ x1 + x2 + x3 + x4 + x5
)

ranking <- borg_compare_models(d, cv, models = model_list, metric = "rmse")
ranking
```

The result is a data frame sorted by mean metric value (ascending for error
metrics, descending for performance metrics), with a rank column. The standard
deviation across folds gives a rough sense of how stable each model's performance
is. Stability counts. A model with slightly better mean RMSE but much higher
variability across folds may not be the safest choice.

```{r compare-models-plot, fig.height = 4}
autoplot(ranking)
```

The bar chart shows mean metric values with error bars spanning one standard
deviation. Labels on the bars give exact values. The models are ordered by rank,
so the best performer appears at the top.

You can pass pre-fitted model objects instead of formulas if your models require
non-standard fitting procedures. In that case, omit `fit_fun` and let
`borg_compare_models()` call each model's own `predict()` method on the test
folds.

# Conformal prediction

Classical prediction intervals assume a parametric noise distribution (usually
Gaussian) and rely on asymptotic theory. Conformal prediction produces
distribution-free intervals with finite-sample coverage guarantees: if you
request 90% intervals, at least 90% of future observations will fall inside them,
with no distributional assumptions required. The only assumption is
exchangeability.

Spatial data violates exchangeability. Nearby observations are more similar than
distant ones, so residuals computed on a random calibration set underestimate the
true prediction error at new locations. `borg_conformal()` fixes this by computing
nonconformity scores from spatially blocked calibration residuals. Because the
calibration residuals come from predictions on held-out spatial blocks, they
reflect the model's actual error on genuinely separated data.

```{r conformal-split}
conf <- borg_conformal(
  model  = model,
  data   = d,
  target = "y.1",
  coords = c("x", "y"),
  alpha  = 0.10,       # 90% intervals
  method = "split"
)
conf
```

The print method reports the nominal and empirical coverage, the number of
calibration scores, the interval half-width (the conformal quantile q-hat), and
summary statistics of the interval widths. Should the empirical coverage fall
substantially below the nominal level, the spatial blocking was too coarse or the
model's errors are spatially heterogeneous.

Three conformal methods are available. `"split"` divides the data into training and
calibration blocks, computes scores on the calibration blocks, and generates
intervals. It is fast but uses only part of the data. `"block_jackknife"` refits
the model leaving out each block in turn, producing more calibration scores at the
cost of repeated model fitting. `"block_cv"` uses residuals from pre-computed
blocked CV folds, which is the most natural integration with the rest of the BORG
pipeline.

```{r conformal-blockcv}
conf_cv <- borg_conformal(
  model  = model,
  data   = d,
  target = "y.1",
  coords = c("x", "y"),
  alpha  = 0.10,
  method = "block_cv",
  folds  = cv
)
conf_cv
```

The `autoplot()` method offers two views. The default shows prediction intervals
sorted by predicted value, which makes it easy to spot whether interval width
varies systematically with the prediction magnitude.

```{r conformal-intervals, fig.height = 4.5}
autoplot(conf_cv, type = "intervals")
```

When coordinates are available, the map view shows interval width at each location,
revealing whether uncertainty is spatially clustered. Wider intervals in certain
regions suggest the model has less information there, perhaps due to sparse
training data or stronger local noise.

```{r conformal-map, fig.height = 5}
autoplot(conf_cv, type = "map")
```

To generate intervals for new locations, pass a `new` data frame. The conformal
quantile is calibrated on the existing data and applied to predictions at the new
points.

```{r conformal-new}
new_locs <- data.frame(
  x  = runif(50, 0, 1),
  y  = runif(50, 0, 1),
  x1 = rnorm(50), x2 = rnorm(50), x3 = rnorm(50),
  x4 = rnorm(50), x5 = rnorm(50)
)
pred_new <- borg_conformal(
  model  = model,
  data   = d,
  new    = new_locs,
  target = "y.1",
  coords = c("x", "y"),
  alpha  = 0.10
)
head(pred_new)
```

The intervals have a constant width (equal to 2 * q-hat) because split conformal
prediction uses the same quantile for all predictions. Adaptive intervals that widen or narrow based on local difficulty require normalised conformal
prediction or locally-weighted variants, though BORG does not
currently implement those.

# Bootstrap uncertainty

A single cross-validation run produces a point estimate of model performance:
one RMSE, one R-squared. But how precise is that estimate? Repeating the
entire experiment with new data would yield a different number.
`borg_bootstrap()` answers this by constructing confidence intervals through
spatial block bootstrap resampling.

Standard bootstrap assumes independent observations. For spatial data, resampling
individual observations produces bootstrap samples where nearby points are
duplicated together, understating the true sampling variability. Block bootstrap
instead resamples entire spatial blocks (defined by k-means clustering on
coordinates), preserving within-block correlation while generating valid
resamples of the data.

```{r bootstrap}
boot <- borg_bootstrap(
  model   = model,
  data    = d,
  target  = "y.1",
  coords  = c("x", "y"),
  formula = y.1 ~ x1 + x2 + x3 + x4 + x5,
  metric  = "rmse",
  n_boot  = 200,
  n_blocks = 10,
  conf_level = 0.95,
  seed    = 42
)
boot
```

The output reports the point estimate (from leave-one-block-out CV on the original
data), the 95% confidence interval, bootstrap standard error, and bootstrap bias.
Bias quantifies the systematic difference between the bootstrap mean and the point
estimate; large bias suggests the point estimate is optimistic or pessimistic in a
consistent direction.

```{r bootstrap-plot, fig.height = 4}
autoplot(boot)
```

The histogram shows the distribution of bootstrap metric values, with the solid
red line marking the point estimate and dashed blue lines marking the confidence
interval bounds. Ideally the distribution is symmetric and centred on the point
estimate. Heavy skew or multimodality suggests instability, often because a small
number of blocks dominate the resampling.

Block count matters. The `n_blocks` parameter controls resampling granularity:
too few blocks (< 5) means each resample is nearly identical to the original data,
underestimating uncertainty, while too many blocks risks splitting spatially
coherent groups into fragments that no longer preserve within-block correlation.
Start with 10 blocks for datasets of 200-500 observations and scale up for
larger ones.

When you already have blocked CV folds from `borg_cv()`, pass them via the `folds`
argument so the bootstrap uses the same block structure for both the point estimate
and the resampling.

# Power analysis

Switching from random to blocked cross-validation reduces the effective sample
size. Observations within the same block are not independent, so the amount of
information in 300 spatially clustered observations is less than in 300
independent ones. This reduction is captured by the design effect (DEFF): the
ratio of the actual sample size to the effective sample size. A DEFF of 3 means
300 observations carry as much information as 100 independent ones.

`borg_power()` computes the design effect from the dependency structure, estimates
the effective sample size, and reports how much statistical power you lose by
switching to blocked CV. It also computes the minimum detectable effect size at
your target power level (default 80%), so you know the smallest signal your
experiment can reliably detect.

```{r power}
pw <- borg_power(
  data   = d,
  coords = c("x", "y"),
  target = "y.1",
  alpha  = 0.05,
  power  = 0.80
)
pw
```

Printed results include the actual sample size, the effective sample size after blocking,
the design effect, and the minimum detectable effect size under both random and
blocked CV. The effect size ratio tells you how much larger a signal needs to be
for blocked CV to detect it at the same power. An effect size ratio of 1.5 means
you need a 50% stronger signal with blocked CV compared to random CV.

When you supply a specific effect size via the `effect_size` argument, the function
computes the actual power under both CV schemes for that effect size, plus the
absolute power loss from switching to blocked CV.

```{r power-effect}
pw2 <- borg_power(
  data        = d,
  coords      = c("x", "y"),
  target      = "y.1",
  effect_size = 0.5
)
cat("Power (random CV):", round(pw2$power_random * 100, 1), "%\n")
cat("Power (blocked CV):", round(pw2$power_blocked * 100, 1), "%\n")
cat("Power loss:", round(pw2$power_loss * 100, 1), "percentage points\n")
```

The design effect is computed differently depending on the dependency type. For
spatial data, BORG uses Moran's I and the effective sample size from the spatial
diagnosis. For temporal data, it uses the lag-1 autocorrelation with the formula
DEFF = (1 + rho) / (1 - rho). For clustered data, it uses the intraclass
correlation coefficient (ICC) and the mean cluster size with the formula
DEFF = 1 + (m - 1) * ICC. For mixed dependencies, the components are combined
multiplicatively, which is conservative but safe.

A common situation in ecological studies is that the dataset technically passes all
leakage checks but has a DEFF above 5, leaving insufficient power under blocked CV
to detect the expected effect. Falling back on random CV is not the answer, because
random CV gives invalid inference. Better options: collect more data, reduce the
spatial grain by aggregating observations, or simplify the model. The
`recommendation` field in the output provides guidance tailored to the specific
diagnosis.

Plan before you collect. Power analysis is most useful at the design stage, before
data collection.
Simulate plausible data with `borg_simulate()`, run `borg_power()`, and check
whether the planned experiment will have sufficient power under blocked CV. This
prevents a frustrating sequence: collecting expensive field data, discovering that
blocked CV is required, then finding out that the dataset is too small.

# Per-fold performance

Aggregated metrics like mean RMSE hide spatial variation in model performance.
A model might predict well in data-dense regions and poorly in sparsely sampled
areas, but the mean metric masks this pattern entirely, which is why per-fold
inspection matters. `borg_fold_performance()` evaluates the model on each fold
separately and records the metric value along with the spatial centroid of each
fold's test set.

```{r fold-performance}
fold_perf <- borg_fold_performance(
  data    = d,
  folds   = cv,
  formula = y.1 ~ x1 + x2 + x3 + x4 + x5,
  coords  = c("x", "y"),
  metric  = "rmse"
)
fold_perf
```

The result is a data frame with one row per fold, showing the metric value, the
training and test set sizes, and the spatial centroid of the test set. If one fold
has RMSE twice as high as the rest, that fold's spatial region deserves closer
inspection.

```{r fold-map, fig.height = 5}
autoplot(fold_perf)
```

When coordinates are available, `autoplot()` produces a spatial performance map
where each point represents a fold's test set centroid, coloured by metric value
and sized by the number of test observations. Fold labels help identify which
region is problematic, and hovering (in interactive contexts) shows the metric
value and fold index. Without coordinates, the fallback is a bar chart of
per-fold metrics with a dashed line at the mean.

Spatial variation in performance often points to nonstationarity: the relationship
between predictors and response changes across space, and the model cannot capture
this shift with a single set of coefficients. If one region consistently performs
poorly, consider fitting region-specific models, adding spatial interaction terms,
or including covariates that capture the local process.

# Diagnostic workflow

The tools described above are most effective when applied in a fixed order. This
checklist summarises the recommended diagnostic workflow for any spatial modelling
project.

**Step 1: Detect and block.** Run `borg_check()` on the raw data to identify
dependency structure and generate blocked folds. If risks are detected, address
them before proceeding. This takes one line of code and should be the first thing
you do after loading data.

**Step 2: Quantify inflation.** Run `borg_compare_cv()` to measure how much random
CV would inflate your metrics. If the inflation is negligible (< 5%), spatial
blocking is still correct but the practical consequences are small. If inflation
exceeds 10%, random CV results should not be reported.

**Step 3: Test significance.** Run `borg_null_test()` with the blocked CV folds to
confirm the model's performance exceeds random expectation. A non-significant
result stops the workflow: there is no point optimizing a model that has not
demonstrated a signal.

**Step 4: Check residuals.** Run `borg_check_residuals()` to verify the model has
absorbed the spatial structure. If residuals show significant autocorrelation,
revise the model (add spatial terms) and return to Step 3.

**Step 5: Assess calibration.** Run `borg_calibration()` with spatial-aware binning
to check whether the model's uncertainty estimates are trustworthy. Poor
calibration does not affect point prediction accuracy but invalidates any
confidence or prediction intervals.

**Step 6: Quantify uncertainty.** Run `borg_bootstrap()` for confidence intervals
on your primary metric, and `borg_conformal()` for prediction intervals at new
locations. Report these intervals alongside point estimates.

**Step 7: Check power.** Run `borg_power()` to assess whether the dataset has
sufficient effective sample size for the detected effect. Report the design effect
and effective sample size alongside p-values.

**Step 8: Map performance.** Run `borg_fold_performance()` to identify spatial
regions where the model struggles. This step is optional but highly informative for
applied work where predictions will be used across the full study area.

# Practical guidance

## When to run which diagnostic

Not every project needs every diagnostic. Here is a decision guide.

For manuscripts reporting spatial model results, Steps 1-3 are needed
at minimum (detect, quantify inflation, null test) plus Step 6 for uncertainty
estimates. Reviewers in ecology and remote sensing increasingly expect spatial CV
comparisons and permutation p-values as standard evidence of valid evaluation.

Predictive models headed for deployment demand Steps 4-5 (residuals and
calibration). A deployed model that produces uncalibrated uncertainty
estimates will mislead downstream decisions, even if its point predictions are
accurate.

At the design stage, before data collection, Step 7 (power analysis)
alone may be sufficient. Simulate plausible data, check that blocked CV would leave
enough effective sample size, and adjust the sampling design before going to the
field.

## Common pitfalls

**Using too few null permutations.** With `n_null = 19`, the minimum achievable
p-value is 0.05, which means you cannot distinguish a marginally significant result
from a highly significant one or report the fine-grained p-values that reviewers
expect. Resolution matters. Use `n_null = 99` for exploratory work and
`n_null = 999` for publication, since for linear models each permutation is cheap
and the cost scales linearly.

**Ignoring the design effect.** A dataset of 500 observations with a design effect
of 10 behaves like 50 independent observations, and reporting standard errors
computed from n = 500 overstates precision by a factor of sqrt(10). Always report
the effective sample size when using blocked CV.

**Treating calibration as optional.** Calibration errors compound through
downstream decisions. If a species distribution model reports 90% habitat
suitability at a location and the true value is 60%, conservation resources get
allocated to the wrong areas. The question is not whether the model predicts well
on average, but whether you can trust what the model says about its own confidence
at any given point.

**Over-blocking.** Aggressive blocking (many large blocks) reduces the effective
sample size and can make the model appear worse than it truly is. Block at the
scale of the dependency, not larger. BORG's `borg_diagnose()` estimates the
autocorrelation range and recommends a block size, but verify that recommendation
against domain knowledge about the spatial process driving the autocorrelation.

**Comparing models on different folds.** Always use `borg_compare_models()` with a
shared fold set rather than generating new folds for each model. Different fold
assignments introduce fold-level variance that confounds the comparison, making it
impossible to tell whether performance differences come from the models or from the
fold assignments.

**Reporting random CV alongside blocked CV.** This sounds transparent but is
misleading. Reporting both "R-squared = 0.65 (random CV)" and
"R-squared = 0.45 (blocked CV)", readers fixate on the higher number.
`borg_compare_cv()` is for your private understanding of the inflation magnitude;
the only metric you should report publicly is the one from the valid (blocked)
evaluation scheme.

## Integration with external frameworks

BORG's blocked CV folds can be exported to `rsample`, `caret`, or `mlr3` format
via the `output` argument of `borg_cv()`. This means you can use BORG for
diagnosis and fold generation while running the actual model fitting in your
preferred framework. The diagnostics in this vignette work with any model that has
`predict()` and `residuals()` methods, which covers `lm`, `glm`, `ranger`,
`xgboost`, and most packages in the R modelling ecosystem.
