Reproducing Spatial Information Explained

1 Method Overview & Reproduction Scope

1.1 Core Idea

A spatial prediction model is usually judged by RMSE or R²: how far its predictions are from the observations, averaged over the study area. Those numbers do not say whether the model has captured the geography of the response. Two models with the same RMSE can leave very different residuals — one spatially random, the other still organised into regional clusters, gradients or zones of unequal spread. Spatial information explained (SIE) measures that difference directly.

SIE treats spatial information as a measurable quantity: the mutual information I(Y; S) between a variable and its location, i.e. how much knowing where an observation is reduces the uncertainty about its value. A model that learns the spatial structure of the response leaves residuals that carry little of it. SIE is the relative reduction from the response to the residuals, E = 1 − I(ε; S) / I(Y; S), computed from out-of-fold residuals of spatial block cross-validation and aggregated over several grid resolutions.

• Research problem: validation metrics quantify error, not explained spatial structure; residual diagnostics such as Moran's I detect remaining dependence but not its share of the total; and existing spatial validation metrics are computed under random cross-validation, where nearby observations leak between training and test sets.
• One-line contribution: a normalised metric of explained spatial information — block-CV residuals + mutual information of the response and the residuals with location + a softmax-weighted multi-scale aggregate that needs no scale to be chosen.
• Why it matters here: across ten models, SIE rankings differ from accuracy rankings (Spearman ρ = −0.38 with RMSE, 0.53 with R²): random forest has the second-lowest RMSE but ranks seventh in E. Random cross-validation overstates E for all ten models, by 2.1–17.0 % (mean 9.9 %). This tutorial regenerates all of it from the published code.

1.2 Method Logic

SIE is computed in five steps. The notation follows the paper's Section 2 and the equation numbers are the paper's; every step is exercised by the pipeline in §3.

Step 1 — residuals from spatial block cross-validation

The study area is divided into contiguous blocks that are randomly assigned to k folds. Each fold is predicted by a model trained on the other folds, and the held-out predictions are concatenated into Ŷ_OOF, aligned with Y. Residuals from a model fitted on all the data would overstate SIE, because spatial autocorrelation lets training points inform predictions at nearby locations.

Step 2 — discretization

Y and ε are each cut into n_b equal-frequency bins at their empirical quantiles, which makes the entropy robust to skewed distributions. Location is discretized jointly: the bounding box of the data is split into a regular g × g grid and every observation takes the label of its cell. Empty cells carry no probability and are dropped; n_s is the number of occupied cells.

Step 3 — spatial information as mutual information

Spatial information is the reduction in uncertainty about Y that comes from knowing its spatial unit, estimated from the discrete joint distribution in nats:

The same estimator applied to the residuals gives I(ε; S), the spatial information the model leaves behind. Mutual information reacts to any difference between the distribution inside a cell and the overall distribution — in mean, spread, skewness or tails — so it detects, for example, residuals whose variance changes across space even when their mean and autocorrelation do not.

Step 4 — SIE at one grid resolution

E = 1: the residuals hold no spatial information left; E = 0: the model has not reduced it (for a no-model baseline ε = Y, so E = 0 exactly). Values below zero, when residuals carry more spatial information than the response, are set to zero and signal misspecification. Mutual information is not additive, so E is a relative reduction, not a share of a fixed quantity split between model and residuals.

Step 5 — multi-scale SIE with information weights

Spatial dependence acts at several scales, so E_r is evaluated on a set of grids {r₁, …, r_K} and combined with weights from the softmax of the response's spatial information at each scale (τ = 1):

The weights depend on Y only, so every model is scored with the same weights and no scale is picked by looking at model performance — picking the resolution that maximises E_r would bias SIE upwards.

1.3 The Two Experiments — Purpose, Results, Interpretation

The paper validates SIE on a simulation with known amounts of spatial learning, then applies it to ten machine learning models on real data. This tutorial re-runs the case study; the simulation is shown from the published article for reference.

Experiment 1 — simulation with a controlled learning ladder

In the paper, E is zero both without spatial structure and when structure is present but not learned, and rises monotonically with the learned proportion of the signal in every configuration. Full learning stops at about 0.80 rather than 1: with residuals that are pure noise, the plug-in estimate of I(ε; S) is still positive — the finite-sample noise floor of the estimator, close to the leading bias term — so the ceiling of E depends on the bins, the grid and the sample size. The simulation is not part of this tutorial; its design and results are summarised in the SIE article.

Experiment 2 — C4 grass area across Australia, ten models

Demo run · further-analyses.R The data of the case study, drawn from c4-data.csv: the response Y and the three predictors at 0.5°, AC4/AC3 ratio (X₁), mean annual temperature (X₂) and precipitation (X₃). Data for 2019 from Luo et al. (2024) and the CRU-JRA climate data.

The reading in one sentence: the three most spatially informative models — KNN, the neural network and XGBoost, all at E ≈ 0.50 — are not the three most accurate ones. SVR and random forest have the lowest RMSE and the highest R², yet random forest explains less spatial information than six other models — its averaged trees smooth the response surface — while KNN reaches the top E with a mid-range RMSE. The regularised linear models are last on both. Under random CV every model gains E, most of all the flexible learners that can exploit spatial proximity (GBM +17.0 %), least the linear ones (+2.1 %).

1.4 Reproduction Scope

Everything tagged Demo run regenerates from this folder with Rscript run-all.R on R 4.6.0. The ten models' SIE and R² match reference-results.csv of the SIE repository to machine precision (largest difference 4.4 × 10⁻¹⁶).

How to use this page

• To reproduce the demo: run Rscript run-all.R. Every figure below regenerates in under a minute.
• To use SIE on your own model: §2.3 gives the data contract — a response, out-of-fold predictions from spatial block CV, and coordinates — and §4.1 the porting steps. sie() itself needs nothing but base R.

2 Setup, Structure & Data Contract

2.1 Environment & Dependencies

The metric is one base-R file; the packages below are only needed to fit the paper's ten models and draw the figures.

# The metric: nothing to install, sie.R uses base R only
# The case study: modelling engines and plotting helpers (from CRAN)
install.packages(c("tidymodels", "bonsai", "ranger", "xgboost", "lightgbm",
                   "kernlab", "glmnet", "kknn", "earth", "nnet",
                   "patchwork", "ggrepel"))

2.2 Project Structure & Configuration

One entry script runs three scripts in order. The four files from the SIE repository are used exactly as published; the tutorial adds the worked example, the further analyses and the runner.

sie/
├── run-all.R                  # entry point: steps 1-5, ~45 s
├── sie.R                      # sie(): the metric, base R      (SIE repository)
├── example.R                  # ten models, SIE, figure, check (SIE repository)
├── c4-data.csv                # 2,480 cells: lon, lat, Y, 3 predictors (SIE repository)
├── reference-results.csv      # E and R2 of the reference run  (SIE repository)
├── worked-example.R           # SIE by hand on a 4 x 4 grid
├── further-analyses.R         # further analyses: 5 figures + results/*.csv
├── figures/                   # sie-vs-r2.jpg + 5 fig-*.png    (generated)
├── results/                   # CSV tables, runtimes, session  (generated)
├── README.md
└── LICENSE                    # GPL-3.0

Settings

The settings sit at the top of example.R; every value is the paper's. The block grid of the folds is set inside its fold code.

N_FOLDS <- 5                  # cross-validation folds
N_BINS  <- 20                 # equal-frequency bins for Y and residuals
SCALES  <- c(4, 6, 8, 10)     # grids of 900, 600, 450 and 360 km cells
SEED    <- 42

# Spatial block CV: a 10 x 10 block grid over the bounding box of the data
xb <- cut(S[, 1], breaks = seq(min(S[, 1]), max(S[, 1]), length.out = 11),
          include.lowest = TRUE, labels = FALSE)

2.3 Data Contract

SIE itself needs three aligned vectors: the observed response, out-of-fold predictions from a spatially separated validation, and coordinates. The case study starts one step earlier, from a table with the predictors.

3 Reproduction Pipeline, Results & Validation

3.1 Pipeline Execution

cd sie                       # this folder
Rscript run-all.R            # full pipeline, ~45 s

Rscript worked-example.R     # step 1 alone, base R, instant
Rscript example.R            # steps 2-4 alone: the script of the SIE repository

== Step 1    worked-example.R: SIE by hand on a 4 x 4 grid ===========
Response Y (rows i = 1..4, columns j = 1..4):
     [,1] [,2] [,3] [,4]
[1,]  1.0  2.6  3.2  4.4
[2,]  1.4  2.2  3.6  4.0
[3,]  1.2  2.4  3.0  4.6
[4,]  1.6  2.0  3.4  4.2

One grid resolution (2 x 2 cells), n_bins = 2:
  model A: I(Y;S) = 0.69315  I(e;S) = 0.00000  E = 1.00000  (by hand 1.00000)
  model C: I(Y;S) = 0.69315  I(e;S) = 0.13081  E = 0.81128  (by hand 0.81128)
  model B: I(Y;S) = 0.69315  I(e;S) = 0.69315  E = 0.00000  (by hand 0.00000)
  ln 2 = 0.69315;  (3/4) ln 3 - ln 2 = 0.13081

Too fine a grid (4 x 4, one observation per cell), model A:
  I(Y;S) = 0.69315  I(e;S) = 0.69315  E = 0.00000   bias term 0.469

Worked example: all three E values equal the hand calculation.

== Step 2-4  example.R: ten models, multi-scale SIE, figure, check ===

Block CV, KNN:
Spatial information explained (SIE)
  observations: 2480 | bins: 20 | scales: 4x4, 6x6, 8x8, 10x10
  multi-scale E = 0.501
 n_grid cells I_Y_S I_e_S   E_r  bias weight
      4    15 0.670 0.258 0.615 0.054  0.199
      6    29 0.829 0.358 0.568 0.107  0.234
      8    49 0.971 0.491 0.495 0.184  0.269
     10    72 1.071 0.664 0.379 0.272  0.298

Multi-scale SIE (E) and R2 of the ten models:
  [... 20 rows: E, E_se, R2, R2_se per model and design -- drawn in figures/sie-vs-r2.jpg ...]

Figure written: figures/sie-vs-r2.jpg
Largest difference from the reference: E 4.44e-16, R2 5.55e-16 -- identical to the reference results

== Step 5    further-analyses.R: further analyses, redrawn ===========
Block CV accuracy: RMSE 8.26 (SVR) to 9.75 (Lasso); R2 0.696 (Lasso) to 0.789 (SVR)
Spearman rho of E with RMSE -0.38, with R2 0.53
Weights 900-360 km: 0.199 0.234 0.269 0.298; linear or uniform weights change E by at most 0.013
Random CV overstates E for 10 of 10 models: 2.1% to 17.0% (mean 9.9%)
Sensitivity, 42 settings: Spearman rho 0.879 to 0.988 (median 0.927)
Figures written: figures/fig-data.png, fig-sie-overview.png, fig-accuracy.png, fig-block-vs-random.png, fig-sensitivity.png

Finished in 45.6 s

3.2 Core Analytical Steps

Five steps: purpose → code → output → result → how to read it. Step 1 is the worked example; steps 2–4 are the four parts of example.R; step 5 redraws the paper's further analyses.

source("sie.R")
i <- rep(1:4, times = 4); j <- rep(1:4, each = 4)   # row, column
S <- cbind(x = j, y = i)
Y <- j + o[cbind(i, j)]                              # gradient + perturbation

Y_hat <- list(A = j,                                 # the gradient, learned
              C = Y - e_C[cbind(i, j)],              # part of it
              B = rep(0, 16))                        # nothing
for (m in names(Y_hat))
  print(sie_scale(Y, Y_hat[[m]], S, n_bins = 2, n_grid = 2))

Response Y (rows i = 1..4, columns j = 1..4):
     [,1] [,2] [,3] [,4]
[1,]  1.0  2.6  3.2  4.4
[2,]  1.4  2.2  3.6  4.0
[3,]  1.2  2.4  3.0  4.6
[4,]  1.6  2.0  3.4  4.2

One grid resolution (2 x 2 cells), n_bins = 2:
  model A: I(Y;S) = 0.69315  I(e;S) = 0.00000  E = 1.00000  (by hand 1.00000)
  model C: I(Y;S) = 0.69315  I(e;S) = 0.13081  E = 0.81128  (by hand 0.81128)
  model B: I(Y;S) = 0.69315  I(e;S) = 0.69315  E = 0.00000  (by hand 0.00000)
  ln 2 = 0.69315;  (3/4) ln 3 - ln 2 = 0.13081

Too fine a grid (4 x 4, one observation per cell), model A:
  I(Y;S) = 0.69315  I(e;S) = 0.69315  E = 0.00000   bias term 0.469

Worked example: all three E values equal the hand calculation.

block_id <- as.integer(interaction(xb, yb, drop = TRUE))   # 72 occupied blocks
fold_of_block <- sample(rep(seq_len(N_FOLDS),
                            length.out = length(unique(block_id))))
fold_block <- fold_of_block[block_id]                      # fold of every cell

rec <- recipe(C4_grass_area_pct ~ ., data = df_model) |>
  step_normalize(all_predictors())
models <- list(
  random_forest = rand_forest(trees = 500) |> set_engine("ranger") |> set_mode("regression"),
  xgboost = boost_tree(trees = 500, learn_rate = 0.1, tree_depth = 3) |>
    set_engine("xgboost") |> set_mode("regression"),
  knn = nearest_neighbor(neighbors = 10) |> set_engine("kknn") |> set_mode("regression")
  # ... lightgbm, svr, elastic_net, lasso, gbm, mars: see example.R
)
oof_block  <- oof_cv(folds_block,  "block")    # fit_resamples(), predictions by row
oof_random <- oof_cv(folds_random, "random")
oof_block$neural_net <- oof_nn(fold_block, "block", SEED + 1)   # 64 units, 5 starts

Step 3 · example.R, step 2

Multi-scale SIE of every model

What this step does

Applies sie() to the out-of-fold predictions of each model under each design, with 20 bins and the four grids. E is computed from all 2,480 residuals together; its standard error comes from the E values computed within each of the five folds. R² is computed within each held-out fold and averaged. The full output of sie() is printed for KNN, the top model.

Code

R — condensed from example.R

E <- sie(Y, oof_block$knn, S, n_bins = N_BINS, scales = SCALES)
E$E          # multi-scale SIE
E$profile    # one row per grid: cells, I_Y_S, I_e_S, E_r, bias, weight

# all models, both designs: E, its fold SE, R2 (fold mean) and its SE
res <- rbind(evaluate(oof_block,  fold_block,  "block"),
             evaluate(oof_random, fold_random, "random"))

Example output

Console output

Block CV, KNN:
Spatial information explained (SIE)
  observations: 2480 | bins: 20 | scales: 4x4, 6x6, 8x8, 10x10
  multi-scale E = 0.501
 n_grid cells I_Y_S I_e_S   E_r  bias weight
      4    15 0.670 0.258 0.615 0.054  0.199
      6    29 0.829 0.358 0.568 0.107  0.234
      8    49 0.971 0.491 0.495 0.184  0.269
     10    72 1.071 0.664 0.379 0.272  0.298

Demo run · further-analyses.R SIE under block CV, as in the paper's Fig. 6: (a) total and residual spatial information at 360 km; (b) scale-specific E_r across the four grids; (c) multi-scale E with one standard error across folds — from 0.501 (KNN) to 0.333 (lasso), while I(ε; S) ranges from 0.646 nats (XGBoost) to 0.840 (lasso).

How to read this result

• Read the profile row by row. At 900 km the response has 0.670 nats of spatial information and KNN leaves 0.258 of it, E_r = 0.615. At 360 km there is more to explain (1.071 nats) and more is left (0.664), E_r = 0.379: fine-grained dependence is harder to absorb, and E_r falls with resolution for all ten models.
• The weights are nearly flat (0.199–0.298) and identical for every model, because they come from Y alone. Spatial information is spread across the four scales rather than concentrated at one, so the aggregate is insensitive to the weighting rule: linear or uniform weights change E by less than 0.013 and leave the ranking unchanged (printed by Step 5).
• The bias column is the guard rail. At 360 km it is 0.272 nats, about a quarter of I(Y; S) — the paper's limit for the finest grid. At 1.071 / log 20, location accounts for about 36 % of the entropy of C4 grass cover at this grain.

Step 5 · further-analyses.R

The paper's further analyses, redrawn

What this step does

Takes the out-of-fold predictions that example.R left in the session and redraws the paper's further analyses without refitting anything: the data maps, the mutual information components and scale profiles of Step 3, RMSE and R² under both designs, the block-vs-random comparison of E, and the sensitivity of E to the discretization over 6 bin counts × 7 grids. Tables of every quantity go to results/; the console prints a short summary.

Code

R — condensed from further-analyses.R

# RMSE and R2 = 1 - SSE/SST within each held-out fold, both designs
acc <- rbind(fold_stat(oof_block,  fold_block,  rmse, "Block CV"),
             fold_stat(oof_block,  fold_block,  r2,   "Block CV"),
             fold_stat(oof_random, fold_random, r2,   "Random CV"))

# 42 settings from the same block CV residuals: no refitting
NB <- c(5, 10, 15, 20, 25, 30)
NG <- c(4, 6, 8, 10, 12, 16, 20)                  # 900 ... 180 km
single <- sie_scale(Y, oof_block[[m]], S, n_bins = nb, n_grid = ng)  # per model, nb, ng
rho <- cor(e, ref_rank, method = "spearman")      # vs multi-scale E at n_b = 20

Example output

Demo run · further-analyses.R Accuracy of the ten models: (a) RMSE and (b) R² = 1 − SSE/SST, means over the five held-out folds with one standard error, under block CV and random CV (as in the paper's Figs. 8 and 11).

Demo run · further-analyses.R Block CV E against random CV E, with one standard error; all ten points lie above the diagonal (paper Fig. 14b).

Demo run · further-analyses.R Sensitivity to the discretization, from the fixed block CV residuals (paper Fig. 10): (a) multi-scale E across bin counts; (b) single-scale E_r across grids at n_b = 20, the four aggregated resolutions shaded; (c) Spearman ρ between the model ranking at each of the 42 settings and the reference ranking (multi-scale E, n_b = 20).

How to read this result

• Magnitude moves, order holds. Refining either the bins or the grid lowers E for every model — mean E_r falls from 0.552 at 900 km to 0.181 at 180 km — because the bias term grows (0.05 to 0.95 nats) faster than I(Y; S) (0.67 to 1.58). The ranking barely changes: the leading group stays on top in 40 of 42 settings, the linear models stay last in all of them. Compare E across models only at identical settings and sample sizes.
• Random CV is optimistic in both metrics, for the same reason. Mean E rises from 0.443 to 0.489, mean RMSE falls from 9.12 to 8.15 and mean R² rises from 0.734 to 0.815. Randomly held-out cells have training neighbours that reproduce local structure the predictors cannot, so residuals look less spatial. The flexible learners profit most (GBM +17.0 %, LightGBM +14.9 %), the linear models least (+2.1 %).
• Accuracy ranks the models differently. SVR has the lowest RMSE (8.26) and the highest R² (0.789), random forest the second-lowest RMSE, while the leading E group — KNN, the neural network and XGBoost — sits mid-range in accuracy. Across the ten models Spearman ρ is −0.38 between E and RMSE and 0.53 between E and R².
• Multi-scale is not single-scale. At 360 km alone the mean E_r is 0.333, at the coarsest grid well above E; the aggregate integrates the grains instead of depending on one.

3.3 Reading the Results Together

• SIE separates models. E runs from 0.333 (lasso) to 0.501 (KNN): the instance-based, smooth nonlinear and tree-ensemble models between 0.459 and 0.501, the linear and spline models between 0.333 and 0.368.
• It is not accuracy in disguise. Spearman ρ = −0.38 with RMSE and 0.53 with R²; random forest is second in RMSE and seventh in E.
• Block CV keeps it honest. Random CV overstates E for all ten models, by 2.1–17.0 % (mean 9.9 %), most for the models that exploit spatial proximity.
• The ranking is robust to the weighting rule (linear or uniform weights change E by less than 0.013) and to 42 discretization settings (ρ ≥ 0.879).

4 Adaptation, Writing & Reproducibility

4.1 Port to Your Domain

SIE applies to any spatial prediction with a continuous response — soil properties, air quality, vegetation cover, land surface temperature, urban indicators — and to any model, from kriging to a deep network: it only reads predictions. The port is one function call after your own block cross-validation:

source("sie.R")                       # base R, no packages

# y_hat: out-of-fold predictions of YOUR model from spatial block CV,
# in the same row order as y and coords
E <- sie(y = mydata$y, y_hat = y_hat,
         coords = mydata[, c("x", "y")],
         n_bins = 20, scales = c(4, 6, 8, 10))
E$E                                    # multi-scale SIE
E$profile                              # check: bias small relative to I_Y_S
max(E$profile$bias / E$profile$I_Y_S)  # the case study stays at about 0.25

What to edit

1. The folds. Build spatial blocks for your study area and keep every model on the same folds. The case study's 10 × 10 block grid over the bounding box is a simple default; blocks larger than the range of spatial dependence of the residuals give a stricter test.
2. The models. Collect the out-of-fold predictions of each model you compare, aligned with the data rows. Fix or nest-tune hyperparameters inside the training folds; avoid settings that make the fits unstable, since unstructured noise also changes E.
3. n_bins. 20 suits a few thousand observations. Fewer observations or many tied values call for fewer bins.
4. scales. From a coarse grid with a handful of occupied cells to the finest grid at which the bias term stays small relative to I(Y; S). Fix them before looking at the models, and do not choose the scale that makes a favoured model look best.
5. Nothing else: tau stays at 1.

Choosing grids and bins with the bias rule

• Compute the profile of the response first. sie(y, y_hat = rep(0, length(y)), coords) returns I(Y; S), occupied cells and the bias term for every candidate grid without any model; drop grids whose bias exceeds about a quarter of I(Y; S).
• Mind the occupied cells. Irregular study areas leave many grid cells empty — 72 of 100 in the Australian case — and only occupied cells count in the bias term.
• Same settings for every model. E values are comparable only under identical bins, grids and sample size; report them with those settings.

4.2 Write the Paper

The published study separates the metric from its evidence: define spatial information and SIE, check the metric where the answer is known, then use it to compare real models along a dimension accuracy does not cover.

Manuscript skeleton

1  Introduction        - spatial validation relies on error metrics; no measure
                          of how much spatial structure a model explains, and none
                          designed for block cross-validation (aims)
2  SIE metric          - 2.1 spatial information I(Y;S)
                          2.2 E = 1 - I(e;S)/I(Y;S); discretization and bias
                          2.3 residuals from spatial block CV
                          2.4 multi-scale SIE with softmax weights
3  Simulation          - controlled learning ladder; boundary cases; robustness
4  Case study          - data; ten models; block and random CV; SIE settings
5  Results             - SIE per model; SIE vs accuracy; single vs multi-scale;
                          sensitivity; block vs random CV
6  Discussion          - what mutual information captures; joint reading with
                          accuracy; position among validation metrics; limits
7  Conclusions         - 6-8 sentences answering the aims

4.3 Final Reproducibility Package

What is in the download, and what is not

References & credits

Every figure on this page was generated by the scripts in this folder from the code and data of the SIE repository.

1. Song Y (2026). Spatial information explained by prediction models under block cross-validation. ISPRS Journal of Photogrammetry and Remote Sensing 242:984–1000. doi:10.1016/j.isprsjprs.2026.09.034 (source of the method)
2. SIE code, data and example — github.com/yongzesong/SIE (sie.R, example.R, c4-data.csv and the reference results, as published with the paper)
3. Luo X, Zhou H, Satriawan TW, Tian J, Zhao R, Keenan TF, Griffith DM, Sitch S, Smith NG, Still CJ (2024). Mapping the global distribution of C4 vegetation using observations and optimality theory. Nature Communications 15:1219. doi:10.1038/s41467-024-45606-3 (C4 grass area and the AC4/AC3 ratio)
4. Harris I, Osborn TJ, Jones P, Lister D (2020). Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Scientific Data 7:109. doi:10.1038/s41597-020-0453-3; Kobayashi S, et al. (2015). The JRA-55 reanalysis. Journal of the Meteorological Society of Japan 93:5–48. doi:10.2151/jmsj.2015-001 (temperature and precipitation)
5. The SIE method article — yongzesong.com/sie/ (definition, properties, evaluation and related methods)