How to run spatial block cross-validation¶

You need an honest out-of-sample comparison between spatial models, and az.loo will not give you one — leave-one-out assumes observations are exchangeable given the parameters, which spatial dependence breaks. spatial_kfold holds out spatially contiguous blocks instead, refits on the remainder, and scores the held-out fold under the full-data joint implied by the model [Roberts et al., 2017].

The scoring rule and its derivation are in Architecture.

!pip install git+https://github.com/ljwolf/geovalidate/ --no-deps

Collecting git+https://github.com/ljwolf/geovalidate/
  Cloning https://github.com/ljwolf/geovalidate/ to /tmp/pip-req-build-amecmuag
  Running command git clone --filter=blob:none --quiet https://github.com/ljwolf/geovalidate/ /tmp/pip-req-build-amecmuag
  Resolved https://github.com/ljwolf/geovalidate/ to commit ab72fa3336a6bcecefd4b69184a06307df1fe24e
  Installing build dependencies ... ?25l-
 \
 |
 done
?25h  Getting requirements to build wheel ... ?25ldone
?25h  Preparing metadata (pyproject.toml) ... ?25ldone
?25hBuilding wheels for collected packages: geovalidate
  Building wheel for geovalidate (pyproject.toml) ... ?25ldone
?25h  Created wheel for geovalidate: filename=geovalidate-0.1.0-py3-none-any.whl size=68551 sha256=bb3a33ec9bf9b56be56045d60a3297154d9059b5bda77c8c57689639fa672d91
  Stored in directory: /tmp/pip-ephem-wheel-cache-krn42y_f/wheels/07/e1/1c/505deea87953881f24f62e5ffe6dc847b905492d550844705a
Successfully built geovalidate
Installing collected packages: geovalidate
Successfully installed geovalidate-0.1.0

Set up data with a known answer¶

We simulate from a SAR data-generating process on a regular grid so the “correct” model is known in advance. A well-behaved CV criterion should prefer the SAR over the misspecified OLS, SLX, and SEM alternatives.

import numpy as np
import pandas as pd
from libpysal.graph import Graph

from neighbayes.dgp import simulate_sar
from neighbayes.diagnostics import spatial_kfold
from neighbayes.models import OLS, SAR, SEM, SLX

SEED = 0
GRID = 12  # 12 x 12 = 144 cells
RHO_TRUE = 0.6
BETA_TRUE = np.array([1.0, 2.0, -1.5])

gdf = simulate_sar(
    n=GRID,
    rho=RHO_TRUE,
    beta=BETA_TRUE,
    sigma=1.0,
    seed=SEED,
    create_gdf=True,
    geometry_type="polygon",
)

W = Graph.build_contiguity(gdf, rook=True).transform("r")
print(f"n = {len(gdf)}, true rho = {RHO_TRUE}, columns = {list(gdf.columns)}")

n = 144, true rho = 0.6, columns = ['y', 'X_0', 'X_1', 'X_2', 'geometry']

Fit the models you want to compare¶

We fit four candidates — only the SAR matches the DGP — using short MCMC runs sufficient for a pedagogical example.

FIT_KW = dict(draws=400, tune=400, chains=2, random_seed=SEED, progressbar=False)
FORMULA = "y ~ X_1 + X_2"

models = {
    "OLS": OLS(formula=FORMULA, data=gdf, W=W),
    "SLX": SLX(formula=FORMULA, data=gdf, W=W),
    "SAR": SAR(formula=FORMULA, data=gdf, W=W, logdet_method="eigenvalue"),
    "SEM": SEM(formula=FORMULA, data=gdf, W=W, logdet_method="eigenvalue"),
}
for name, m in models.items():
    m.fit(**FIT_KW)
    print(f"  fit {name}")

  fit OLS
  fit SLX
  fit SAR
  fit SEM

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 2 chains for 400 tune and 400 draw iterations (800 + 800 draws total) took 5 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 2 chains for 400 tune and 400 draw iterations (800 + 800 draws total) took 5 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

Score them out of sample¶

spatial_kfold derives 5 spatial blocks from the cell centroids via KMeans, refits each model on each train fold, and accumulates the held-out elpd. Higher elpd is better.

CV_KW = dict(draws=200, tune=200, chains=1, random_seed=SEED, progressbar=False)

results = {}
for name, m in models.items():
    results[name] = spatial_kfold(m, geometry=gdf.geometry, n_blocks=5, **CV_KW)

summary = pd.DataFrame(
    {
        "elpd": [r.elpd for r in results.values()],
        "se": [r.se for r in results.values()],
        "n_folds": [r.n_folds for r in results.values()],
    },
    index=list(results.keys()),
).sort_values("elpd", ascending=False)
summary["delta_elpd"] = summary["elpd"] - summary["elpd"].max()
summary

With data simulated from a SAR DGP, spatial block CV should rank SAR on top, with the misspecified alternatives trailing behind. The delta_elpd column shows the elpd difference relative to the best model — values close to zero indicate near-equivalent predictive performance, while large negative values indicate worse out-of-block prediction.

Swap in a different fold design¶

The built-in KMeans blocking is convenient but coarse. For a richer menu of spatial cross-validation designs — Hilbert space-filling curves, DGGS-cell stratification, clustered hold-outs, and exclusion-buffered leave-one-out — spatial_kfold accepts any sklearn-style splitter via the splitter= keyword. The geovalidate package provides a comprehensive set of these splitters.

Install with pip install geovalidate.

Below we re-run CV using HilbertKFold, which interleaves observations along a Hilbert curve so every fold covers the entire study area — typically a more balanced design than KMeans blocking.

from geovalidate import HilbertKFold

splitter = HilbertKFold(n_splits=5, random_state=SEED)

results_hilbert = {
    name: spatial_kfold(m, splitter=splitter, geometry=gdf.geometry, **CV_KW)
    for name, m in models.items()
}

summary_hilbert = pd.DataFrame(
    {
        "elpd": [r.elpd for r in results_hilbert.values()],
        "se": [r.se for r in results_hilbert.values()],
        "method": [r.method for r in results_hilbert.values()],
    },
    index=list(results_hilbert.keys()),
).sort_values("elpd", ascending=False)
summary_hilbert["delta_elpd"] = summary_hilbert["elpd"] - summary_hilbert["elpd"].max()
summary_hilbert

The method column now records the splitter class name (HilbertKFold) instead of "kmeans". Any other geovalidate splitter slots in the same way — for example LeaveClusterOut for leave-one-region-out CV, or BallKFold for spatially-exclusive folds with an explicit exclusion radius.

See also¶

• How to run Bayesian LM specification tests — the in-sample counterpart, cheaper and often enough

• Architecture — the held-out scoring rule

• Supported Models — the candidate models compared here