Spatial Epi Public Health GIS

Cross-Validating an Exposure Surface With Held-Out Monitors

This guide, part of Interpolating Environmental Exposure Surfaces, fixes a specific and very common overstatement: an exposure surface reported with a cross-validated fit of 0.85 that in practice predicts poorly everywhere a monitor is not.

Problem Context & Constraints

Cross-validation assumes the held-out observations are independent of the training set. Spatial data violate that assumption by construction — that is the whole premise of interpolating them. Shuffle a monitor network into five random folds and most held-out sites still have a training-set neighbour a few kilometres away, so the model is asked to interpolate a short distance, which it does well. The reported error is the error of a task the surface is never actually asked to perform.

The task it is asked to perform is prediction at residences, and residences are not distributed like monitors. Regulatory networks are concentrated in population centres and near known sources; large rural areas have no site within the correlation range. The relevant question is therefore how well the surface predicts at a location with no nearby monitor, and only a validation design that removes the neighbours can answer it.

Three constraints follow. The holdout unit must be the site, not the observation — holding out a single day from a site that contributes 365 rows measures nothing. Any tuning that looks at the data must happen inside the fold, including variable selection, variogram fitting and exponent selection; a model whose variogram was fitted on all sites has already seen the held-out one. And the error must be reported in the units of the exposure, not only as a correlation, because a correlation of 0.8 with a two-microgram bias is a different product from a correlation of 0.8 without one.

Random Folds Versus Spatially Blocked Folds Two maps of the same monitor network. On the left, random assignment to folds: the held-out sites are scattered among training sites, so each one has a training neighbour within a few kilometres and the model only has to interpolate a short distance. On the right, spatially blocked assignment: an entire block of the study area is held out at once, so the nearest training site to a held-out site is far away. Below each map, the error the design reports is described — the left reports the error of short interpolation, the right reports the error that applies in a region with no monitor. The fold design decides which question the number answers Random folds every held-out site keeps a close neighbour reports: error of interpolating a few kilometres Spatially blocked folds held-out block nearest training site is far from every holdout reports: error where there is no monitor training site held-out site

Prerequisites

  • python 3.11
  • numpy 1.26.4, geopandas 0.14.4, scipy 1.13.1
  • A surface-fitting function that takes a training frame and target geometries and returns predictions — the kriging routine from the parent guide satisfies this
  • Monitor sites in a projected coordinate reference system in metres, with a stable site identifier
  • An estimate of the correlation range, used to size the blocks

Step-by-Step Solution

Implement both designs and report both. Leave-one-site-out is the right default for small networks; spatial blocking is what you report when the study area contains large unmonitored regions.

# Leave-one-site-out and spatially blocked cross-validation for an exposure surface.
# Pinned: numpy==1.26.4, geopandas==0.14.4
import logging

import numpy as np
import geopandas as gpd

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("exposure.cv")


def _scores(obs, pred):
    resid = obs - pred
    ss_res = float((resid ** 2).sum())
    ss_tot = float(((obs - obs.mean()) ** 2).sum())
    return {"n": int(len(obs)),
            "rmse": float(np.sqrt((resid ** 2).mean())),
            "mae": float(np.abs(resid).mean()),
            "bias": float(resid.mean()),
            "r2": float(1.0 - ss_res / ss_tot) if ss_tot > 0 else float("nan")}


def loo_cv(sites: gpd.GeoDataFrame, value_col: str, fit_predict) -> dict:
    """Leave-one-site-out. fit_predict(train_gdf, value_col, target_gdf) -> array."""
    sites = sites.sort_values("site_id").reset_index(drop=True)
    obs = sites[value_col].to_numpy(float)
    pred = np.empty(len(sites))
    for i in range(len(sites)):
        train = sites.drop(index=i)
        # The whole model is refitted inside the fold: variogram, trend, selection.
        pred[i] = fit_predict(train, value_col, sites.iloc[[i]])[0]
    s = _scores(obs, pred)
    log.info("loo n=%d rmse=%.4f bias=%+.4f r2=%.3f", s["n"], s["rmse"], s["bias"], s["r2"])
    return s


def blocked_cv(sites: gpd.GeoDataFrame, value_col: str, fit_predict,
               block_m: float, seed: int = 20260918) -> dict:
    """Spatially blocked folds: whole square blocks of the extent are held out."""
    sites = sites.sort_values("site_id").reset_index(drop=True)
    x = sites.geometry.x.to_numpy()
    y = sites.geometry.y.to_numpy()
    bx = np.floor((x - x.min()) / block_m).astype(int)
    by = np.floor((y - y.min()) / block_m).astype(int)
    block_id = bx * 10_000 + by
    blocks = np.unique(block_id)
    rng = np.random.default_rng(seed)
    rng.shuffle(blocks)

    obs, pred, used = [], [], 0
    for b in blocks:
        hold = block_id == b
        if hold.sum() == 0 or (~hold).sum() < 10:
            continue                      # too little left to fit anything
        p = fit_predict(sites[~hold], value_col, sites[hold])
        obs.extend(sites.loc[hold, value_col].to_numpy(float))
        pred.extend(np.asarray(p, dtype=float))
        used += 1
    s = _scores(np.asarray(obs), np.asarray(pred))
    s.update({"blocks_used": used, "block_m": block_m, "seed": seed})
    log.info("blocked block=%.0fm blocks=%d n=%d rmse=%.4f bias=%+.4f r2=%.3f",
             block_m, used, s["n"], s["rmse"], s["bias"], s["r2"])
    return s

The block size is the parameter that decides what the blocked design measures. Blocks smaller than the correlation range leave neighbours across the block edge and the result drifts back toward the random-fold number. Blocks around one to two times the fitted range remove the local information and report the error that applies inside a network gap of that size. Choose the block size to match the size of the actual unmonitored regions in your study area, and state it — a blocked error without its block size is uninterpretable.

The seeded shuffle exists so the fold assignment is reproducible. Block cross-validation is otherwise deterministic given the block size and origin, but the shuffle order affects nothing except the log, which is exactly why it should be fixed rather than left to chance.

Reported Error as a Function of Block Size A curve showing held-out error against the size of the spatial block removed. At a block size of zero the design is equivalent to random folds and reports the lowest, most optimistic error. As the block grows past the fitted correlation range the error climbs steeply, because held-out sites lose every informative neighbour. Beyond roughly twice the range the curve flattens, because the prediction has already fallen back to the regional mean and removing more sites cannot make it worse. A marker shows where the typical unmonitored gap in the study area falls, which is the block size that should be reported. Pick the block size from your real gaps, then report that number fitted range typical gap in this study area report the error here random-fold equivalent error of the regional mean 0 size of the block held out held-out error

Validation & Edge Cases

Selection leakage. The most frequent defect in published exposure models is choosing predictors on all the data and then cross-validating the chosen model. The fix is mechanical: the fold loop must call a function that performs selection, not a pre-selected formula. If refitting inside every fold is too slow to run, that is a performance problem, not a licence to leak.

A fold that cannot be fitted. Blocking sometimes leaves too few training sites for a variogram. Skip the block, count it, and report how many were skipped — a blocked result computed over a third of the blocks is not the same statistic as one computed over all of them.

A network that changes between years. Sites open and close, so the fold structure is not stable across a multi-year surface. Validating only the years with the fullest network reports the best case and hides the worst. Run the design separately for every year the surface is applied to and keep the per-year table; a year in which four sites closed will show it, and that year’s exposure values deserve wider uncertainty in the health model.

Too few sites for blocking at all. With fewer than roughly twenty-five sites, most blocks leave a training set too small to refit, and the blocked statistic is computed from a handful of folds with wide sampling error. In that situation report leave-one-site-out only, say explicitly that blocked validation was not feasible, and treat any use of the surface in unmonitored regions as unvalidated rather than as validated by the optimistic number.

Bias hiding inside a good correlation. Report bias separately and always. A surface that tracks the spatial pattern well but sits systematically low will produce a correlation near 0.9 and an exposure-response slope that is wrong by the ratio of the scales.

Error that varies across the region. A single number hides the thing you most need to know. Map the residuals; if they are large in one part of the state, the assumption of one covariance model for the whole area has failed, and the remedy is in Handling Non-Stationary Variograms Across a State.

INFO exposure.cv loo n=48 rmse=1.1802 bias=-0.0431 r2=0.742
INFO exposure.cv blocked block=150000m blocks=9 n=48 rmse=2.0617 bias=-0.2085 r2=0.213

Those two lines describe the same surface. Reporting only the first would claim the model explains three quarters of the variation; the second says that in a region the size of a real network gap it explains about a fifth. Both are true, and a study that uses the surface in rural counties is governed by the second.

Where the residual map belongs in the workflow is easiest to see laid against the validation outputs themselves:

A Pooled Statistic Can Hide a Regional Failure Held-out residuals at each monitoring site drawn as upward or downward marks across a state. Across the western two thirds residuals are small and mixed in sign. In the eastern third they are consistently positive, meaning the surface under-predicts there. The pooled root mean squared error looks acceptable because the eastern sites are a minority, but the map shows a structured regional bias that a single number conceals. Map the residuals; the pooled number averages the problem away east: all under-predicted small, mixed residuals pooled error looks acceptable; the eastern third needs its own model mark length is the size of the held-out residual; upward means under-predicted

Compliance Notes

The validation record is part of the surface’s provenance, not a separate analysis artifact. Store, alongside the surface: the design used, the block size and seed where applicable, the number of folds attempted and completed, the site count, and root mean squared error, mean absolute error, bias and fit for each design. Where the study applies the surface across several years, store the statistics per year rather than pooled, because a surface that degrades in later years is a finding that pooling hides.

State the validation design in the methods text of any product built on the surface. “Cross-validated R2R^2 = 0.74” is not a sufficient statement — the same surface honestly yields 0.21 under a design that matches how it is used. The distinction is the difference between a defensible claim and one that collapses under review.