How to Evaluate a Spatial Model Honestly
Problem statement
A single score is not an evaluation of a map. It is one number summarising a model's performance on the data it was validated against, and every map has parts where it is much better and much worse than that number.
Measured on a real problem, the same model reported three ways:
random 5-fold CV, inside the sampled area +0.555 sd 0.019
spatial block CV +0.213 sd 0.106
true error over the whole unsampled map +0.557 sd 0.014
the model applied to the adjacent half -0.771
All four are correct answers to different questions. Reporting one without saying which question it answers is the most common failure in applied spatial modelling.
Quick answer
An honest evaluation has five parts:
report = {
"scheme": "spatial block CV, 10 blocks, 5 folds, 5 layouts",
"metric": {"r2": 0.213, "r2_sd": 0.106, "mae": 0.161},
"deployment": "prediction within the sampled extent",
"applicability": "78% of the map is inside the area of applicability",
"residuals": {"morans_i": 0.08, "bias": -0.004},
}The scheme and the deployment are what make the metric interpretable. Without them, R² = 0.213 is not information.
Step-by-step solution
1. State the validation scheme and the deployment task
These are the two facts that let a reader interpret everything else. R² = 0.555, random 5-fold, for gap-filling within the sampled extent is a result; R² = 0.555 is not.
2. Report the spread, not a point estimate
Block CV over ten different block layouts gave a standard deviation of 0.106, ranging from +0.061 to +0.429. Quoting one of those to three decimal places is false precision by a wide margin.
Repeat over layouts, over seeds, over train/test splits, and report the mean and spread.
3. Use metrics in the variable's units
R² compares against predicting the mean, which for a spatially structured variable is a weak baseline. It comes out high for mediocre models and is hard to act on.
MAE and RMSE in the variable's own units — "±0.16 NDVI", "±54 m of elevation" — are directly interpretable and comparable across studies.
Report bias separately: a model can have a good RMSE and a systematic offset.
4. Map the residuals
A single error number hides where the errors are. Two checks:
- Spatial autocorrelation of residuals. Structured residuals mean an unmodelled spatial process. Moran's I on the residuals is the standard test.
- Error against covariate values. Systematically larger errors at high elevation or steep slope means the model fits one part of the range and not another.
5. Report the area of applicability
The share of the map where the covariates resemble the training data. Outside it, the model is extrapolating and the metric does not apply at all.
A map that is 60% inside its area of applicability with an RMSE of 0.1 is a very different product from one that is 98% inside with the same RMSE.
Code examples
Example 1 — a full evaluation report
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.model_selection import GroupKFold, KFold
def evaluate(model_factory, X, y, coords, feature_names,
deployment="prediction within the sampled extent",
n_blocks=10, repeats=5):
"""Scheme, metrics with spread, residual diagnostics, applicability."""
def cv(splitter):
predicted = np.empty(len(y))
for train, test in splitter:
predicted[test] = model_factory().fit(X[train], y[train]) \
.predict(X[test])
return predicted
random_pred = cv(KFold(5, shuffle=True, random_state=0).split(X))
block_scores, block_pred = [], None
for seed in range(repeats):
blocks = KMeans(n_blocks, random_state=seed,
n_init=10).fit_predict(coords)
predicted = cv(GroupKFold(5).split(X, y, blocks))
block_scores.append(r2_score(y, predicted))
if block_pred is None:
block_pred = predicted
block_scores = np.array(block_scores)
residual = y - block_pred
print(f" deployment task : {deployment}")
print(f" random 5-fold : R2 {r2_score(y, random_pred):+.3f}, "
f"MAE {mean_absolute_error(y, random_pred):.4f}")
print(f" block CV : R2 {block_scores.mean():+.3f} "
f"± {block_scores.std():.3f} over {repeats} layouts "
f"(range {block_scores.min():+.3f} to {block_scores.max():+.3f})")
print(f" MAE (block) : {mean_absolute_error(y, block_pred):.4f} "
f"in units of the target")
print(f" bias : {residual.mean():+.4f}")
print(f" target sd : {y.std():.4f} (the baseline RMSE)")
if block_scores.std() > 0.05:
print(" ! block CV is unstable — report the range, not a point value")
return {"random_r2": float(r2_score(y, random_pred)),
"block_r2_mean": float(block_scores.mean()),
"block_r2_sd": float(block_scores.std()),
"mae": float(mean_absolute_error(y, block_pred)),
"bias": float(residual.mean()), "deployment": deployment}Printing the target's standard deviation next to the MAE gives the reader the baseline. An MAE of 0.16 against a target sd of 0.19 is a weak model; against a sd of 1.9 it is a strong one.
Example 2 — residual diagnostics
import numpy as np
from scipy.spatial import cKDTree
def residual_diagnostics(residual, coords, X=None, feature_names=None,
k=8):
"""Is the error spatially structured, or related to a covariate?"""
tree = cKDTree(coords)
_, neighbours = tree.query(coords, k=k + 1)
lag = residual[neighbours[:, 1:]].mean(axis=1)
centred = residual - residual.mean()
lag_centred = lag - lag.mean()
morans = float((centred * lag_centred).sum() /
np.sqrt((centred ** 2).sum() * (lag_centred ** 2).sum()))
print(f" residual mean {residual.mean():+.4f}, sd {residual.std():.4f}")
print(f" spatial autocorrelation of residuals (k={k}): {morans:+.3f}")
if abs(morans) > 0.2:
print(" ! residuals are spatially structured — an unmodelled "
"spatial process remains")
if X is not None and feature_names is not None:
print(" error against each covariate (correlation with |residual|):")
for i, name in enumerate(feature_names):
correlation = float(np.corrcoef(X[:, i], np.abs(residual))[0, 1])
flag = " <- heteroscedastic" if abs(correlation) > 0.3 else ""
print(f" {name:16} {correlation:+.3f}{flag}")
return {"morans_i": morans, "bias": float(residual.mean())}The correlation between a covariate and the absolute residual is the heteroscedasticity check. A model that is accurate on gentle slopes and poor on steep ones has a real limitation that RMSE averages away.
Example 3 — the report that ships with the map
import json
def evaluation_report(metrics, diagnostics, applicability_share,
scheme, deployment, features, n_train, path=None):
"""Everything a reader needs, in one serialisable object."""
report = {
"deployment_task": deployment,
"validation_scheme": scheme,
"n_training_samples": int(n_train),
"features": list(features),
"metrics": {
"r2_mean": round(metrics["block_r2_mean"], 3),
"r2_sd": round(metrics["block_r2_sd"], 3),
"mae": round(metrics["mae"], 4),
"bias": round(metrics["bias"], 4),
},
"residuals": {
"spatial_autocorrelation": round(diagnostics["morans_i"], 3),
},
"area_of_applicability": {
"share_inside": round(applicability_share, 3),
"note": "outside this share the model is extrapolating and the "
"metrics above do not apply",
},
"caveats": [
"metrics are estimated at the training locations",
"error grows with distance from the training data",
f"the model cannot predict outside the training target range",
],
}
text = json.dumps(report, indent=2)
print(text)
if path:
with open(path, "w") as handle:
handle.write(text)
return reportThe caveats list is the part that survives summarisation. When someone reduces the model to a single number in a slide, the caveats are what they should have carried with it — and writing them down at least makes their omission a choice.
Explanation
Why R² is a weak metric here
R² is one minus the ratio of the model's squared error to the variance of the target. It measures improvement over predicting the mean.
For a spatially structured variable, predicting the mean is a very weak baseline, so R² comes out high for models that are not very useful. It also has no units, so it cannot be checked against a requirement — "is ±0.16 NDVI good enough?" is answerable; "is 0.55 good enough?" is not.
Report MAE or RMSE in the target's units, with the target's standard deviation for context.
Why the spread matters as much as the mean
Block CV's standard deviation was 0.106 across ten block layouts on identical data and an identical model — a range of 0.37 R².
Two papers reporting 0.21 and 0.43 for the same method on the same data are not disagreeing about anything except which blocks the clustering produced.
Repeating and reporting the spread is what makes the number reproducible. It also, usefully, tells the reader how much to trust it.
Why residual autocorrelation matters
If the residuals are spatially structured — nearby errors similar — the model has missed a spatial process. Something varies smoothly across the area that the covariates do not capture.
That is actionable in a way an error number is not: it points at a missing covariate, or at a case for adding a spatial term such as kriging the residuals.
Residuals that look like noise mean the covariates have captured the structure, which is the goal.
Why applicability belongs in the metric report
The metric was estimated at the training locations. Wherever the map's covariates fall outside the training range, the metric does not describe the error — it describes the error somewhere else.
A model with an RMSE of 0.1 covering 98% of its map inside the area of applicability is a usable product. The same RMSE covering 55% is a research result with a large asterisk.
Reporting both turns a number into a statement about a map.
Edge cases or notes
- State the scheme and the deployment task with every metric.
- Report the spread over repeats, not a point estimate.
- Use the target's units; give its standard deviation as the baseline.
- Report bias separately from RMSE.
- Test residual autocorrelation — structure means a missing covariate.
- Check error against covariate values for heteroscedasticity.
- Report the applicability share.
- Write the caveats down, because they are the first thing dropped in a summary.
Internal links
- Spatial leakage explained: why random cross-validation lies — why the scheme must be stated
- How to do spatial cross-validation in Python — producing the numbers
- Extrapolation in space explained: the area of applicability — the applicability share
- My model scores 0.95 in testing and fails in the field — what a bare metric hides
- Spatial autocorrelation explained — the residual test
- How to calculate Moran's I in Python — computing it properly
- My prediction raster is striped, blocky or full of NoData — what the map shows that metrics do not
- How to cross-validate an interpolated surface — the same argument for interpolation
FAQ
What should I report for a spatial model?
The validation scheme, the deployment task, a metric in the target's units with its spread, residual diagnostics, and the share of the map inside the area of applicability.
Should I use R² or RMSE?
RMSE or MAE, in the target's units, with the target's standard deviation as context. R² compares against predicting the mean, which is a weak baseline for spatial data.
Why report the spread of the score?
Because block CV over ten block layouts gave a standard deviation of 0.106 on identical data. A single number to three decimals is false precision.
What does spatially autocorrelated residual mean?
That the model has missed a spatial process. Something varies smoothly across the area that the covariates do not capture — which points at a missing covariate.
How do I check for heteroscedasticity?
Correlate each covariate with the absolute residual. A strong correlation means the model is accurate over part of the range and not the rest.
Why report the area of applicability with the metric?
Because the metric was estimated at the training locations. Outside the applicability area it describes the error somewhere else entirely.
What if my model has a good RMSE and a bias?
Report both. A systematic offset is a different problem from scatter, and RMSE conflates them.