Spatial Epi Public Health GIS

Handling Non-Stationary Variograms Across a State

This guide, part of Interpolating Environmental Exposure Surfaces, deals with the assumption that breaks first when an exposure surface is scaled from a metropolitan area to a whole state: that one covariance structure describes the entire region.

Problem Context & Constraints

Ordinary kriging assumes second-order stationarity — the same mean and the same distance-decay of covariance everywhere in the study area. A metropolitan network mostly satisfies it. A state network usually does not, for reasons that are physical rather than statistical: a coastal strip with marine air and a continental interior with wildfire influence genuinely have different variability; an urban core where concentrations change over a kilometre and a rural plain where they change over a hundred are genuinely different fields.

Fitting one variogram across that heterogeneity produces a single compromise model. The consequences are asymmetric and predictable. In the high-variability region the fitted range is too long, so the surface over-smooths and erases real contrast; in the low-variability region the fitted nugget is too high, so the surface pulls toward the mean and the prediction variance is overstated. Neither failure announces itself — both produce a plausible map.

The constraint that makes this hard is sample size. The obvious fix, fitting a separate variogram per region, divides an already-small network. A state with 50 sites split into three regions leaves fitting problems in each. Any remedy therefore has to buy locality without spending all the data, and the three that do are trend removal, regional sub-models with pooled structure, and a moving-window variogram.

One Compromise Model Versus Two Honest Ones A semivariogram plot holding three curves. The coastal region curve rises quickly to a low sill, indicating small total variance that decorrelates over a short distance. The interior region curve rises slowly to a much higher sill, indicating large variance decorrelating over a long distance. The pooled curve fitted to all sites together lies between them and matches neither. Annotations state the consequence in each region: over-smoothing where the pooled range is too long, and an inflated nugget where the pooled model borrows variance from the other region. The pooled model is wrong in both regions, in opposite directions coastal: short range, low sill interior: long range, high sill pooled fit on the coast the pooled range is too long so real contrast is smoothed away inland the pooled sill is too low so the variance band understates the gaps separation distance semivariance

Prerequisites

  • python 3.11
  • numpy 1.26.4, geopandas 0.14.4, scipy 1.13.1
  • A variogram fitting routine — the weighted least-squares fit from the parent guide
  • Monitor sites in a projected metric coordinate reference system with a stable identifier
  • A candidate regionalisation: an ecoregion layer, an airshed boundary, or a simple urban/rural split. It must be defined from physical reasoning, never from the residuals you are about to examine

Step-by-Step Solution

Detect first, then choose the least intrusive remedy that fixes what you detected. The detector below fits the variogram separately in each candidate region and compares the parameters against the pooled fit with a simple ratio test, and then fits a moving-window range to show how the structure varies continuously.

# Detect and describe non-stationarity in an exposure monitor network.
# Pinned: numpy==1.26.4, geopandas==0.14.4, scipy==1.13.1
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.nonstat")

MIN_SITES_PER_REGION = 20
RANGE_RATIO_FLAG = 2.0
SILL_RATIO_FLAG = 2.0


def regional_variograms(sites: gpd.GeoDataFrame, value_col: str,
                        region_col: str, fit_variogram) -> dict:
    """Fit the variogram pooled and per region; flag material disagreement."""
    sites = sites.sort_values("site_id").reset_index(drop=True)
    xy = np.c_[sites.geometry.x.values, sites.geometry.y.values]
    z = sites[value_col].to_numpy(float)
    pooled = fit_variogram(xy, z)
    out = {"pooled": {"nugget": pooled[0], "sill": pooled[1], "range_m": pooled[2] * 3.0},
           "regions": {}}

    for region, idx in sites.groupby(region_col).groups.items():
        i = np.asarray(sorted(idx))
        if len(i) < MIN_SITES_PER_REGION:
            log.warning("region %s has %d sites — too few to fit separately", region, len(i))
            continue
        try:
            n, s, r = fit_variogram(xy[i], z[i])
        except RuntimeError:
            log.warning("region %s variogram did not converge", region)
            continue
        out["regions"][str(region)] = {"n_sites": int(len(i)), "nugget": n,
                                       "sill": s, "range_m": r * 3.0}

    rs = [v["range_m"] for v in out["regions"].values()]
    ss = [v["sill"] for v in out["regions"].values()]
    if len(rs) >= 2:
        out["range_ratio"] = max(rs) / max(min(rs), 1e-9)
        out["sill_ratio"] = max(ss) / max(min(ss), 1e-9)
        out["non_stationary"] = (out["range_ratio"] > RANGE_RATIO_FLAG
                                 or out["sill_ratio"] > SILL_RATIO_FLAG)
        log.info("range_ratio=%.2f sill_ratio=%.2f non_stationary=%s",
                 out["range_ratio"], out["sill_ratio"], out["non_stationary"])
    return out


def detrend_by_position(sites: gpd.GeoDataFrame, value_col: str) -> gpd.GeoDataFrame:
    """Remove a first-order spatial trend; often this is all the fix that is needed."""
    x = sites.geometry.x.to_numpy()
    y = sites.geometry.y.to_numpy()
    A = np.c_[np.ones(len(x)), x - x.mean(), y - y.mean()]
    z = sites[value_col].to_numpy(float)
    coef, *_ = np.linalg.lstsq(A, z, rcond=None)
    resid = z - A @ coef
    out = sites.copy()
    out[f"{value_col}_resid"] = resid
    var_explained = 1.0 - resid.var() / z.var()
    log.info("trend removed: intercept=%.4f dx=%.3e dy=%.3e var_explained=%.3f",
             coef[0], coef[1], coef[2], var_explained)
    if var_explained < 0.05:
        log.info("trend is negligible — do not carry it into the model")
    return out

Run the detrend step before concluding the field is non-stationary. A large-scale gradient across a state — concentrations falling from a populated coast toward an empty interior — makes an empirical variogram climb without ever reaching a sill, which looks exactly like a very long range. Removing a first-order trend and refitting on the residuals frequently collapses that apparent range to something sensible, and universal kriging then handles the rest with the full network intact. This is the cheapest of the three remedies and should always be tried first.

If the ratios still flag after detrending, the structure genuinely differs by region. Fit per-region models where each region clears the site-count minimum, and for regions that do not, borrow the pooled range while fitting only the local sill — a partial pooling that spends far fewer degrees of freedom than a full local fit. Predict each location with the model of its own region, but always from the full set of sites within the neighbourhood, including those across the regional boundary. Restricting the neighbourhood to the region as well produces a visible seam along the boundary, and that seam will end up in a published map.

Three Remedies, Cheapest First Three options laid out in increasing order of cost in degrees of freedom. Trend removal keeps the entire network and models the large-scale gradient as a function of position, costing three parameters. Regional sub-models with a pooled range fit a separate sill per region while sharing the distance-decay, costing one parameter per region. A moving-window variogram fits the structure locally everywhere, costing the most and requiring the densest network. Under each option, the situation it is the right answer to is named, and a note warns that predicting with a region-restricted neighbourhood creates a visible seam on the boundary. Spend the least network you can to fix what you actually detected 1 · Remove a trend whole network retained three extra parameters right when: a smooth gradient crosses the whole state 2 · Regional sills range shared, sill per region one parameter per region right when: variance differs but decay is similar 3 · Moving window structure fitted locally needs a dense network right when: both range and sill vary continuously In all three: predict from every nearby site, including across the regional boundary a neighbourhood restricted to the region draws a seam that will reach a published map Validate the remedy by spatially blocked folds in each region separately, not pooled

Validation & Edge Cases

The regionalisation was chosen from the residuals. Defining regions by looking at where the model fits badly and then fitting those regions separately guarantees an improvement in fit and guarantees nothing about prediction. The boundary must come from physical reasoning — airshed, ecoregion, elevation band, urban extent — and be fixed before any comparison.

A region with too few sites. Below about twenty sites the local variogram is noise. Borrow the pooled range and fit only the sill, or fold the region into its neighbour and note that the surface there inherits a neighbouring structure.

Apparent non-stationarity caused by one site. A single anomalous monitor inside a small region moves its sill enough to trip the ratio test. Refit each region with each site left out in turn and check the ratios are stable; if one site drives the conclusion, the conclusion is about that site.

Seams at the boundary. After fitting regional models, plot a transect across the boundary. A step in the predicted surface means the neighbourhood was restricted along with the model. Predictions must draw on all nearby sites regardless of region.

Validation must also be regional. A pooled blocked-validation number can improve while one region gets worse. Report the held-out error per region, using the design in Cross-Validating an Exposure Surface With Held-Out Monitors.

INFO exposure.nonstat trend removed: intercept=8.4126 dx=-2.104e-06 dy=-7.712e-07 var_explained=0.284
INFO exposure.nonstat range_ratio=2.71 sill_ratio=3.36 non_stationary=True
WARNING exposure.nonstat region mountain has 11 sites — too few to fit separately

The seam test is quick to run and its failure is unmistakable once plotted as a transect crossing the regional boundary:

The Seam Test Across a Regional Boundary Predicted concentration along a line crossing from the coastal region into the interior region. When the prediction neighbourhood is restricted to each region's own sites, the profile jumps abruptly at the boundary, a seam that would appear as a straight edge on a published map. When each location uses its own region's variogram but draws neighbours from all nearby sites regardless of region, the profile changes smoothly across the boundary. A step at the boundary means the neighbourhood was restricted too regional boundary restricted neighbourhood: seam all nearby sites: smooth coastal region interior region predicted concentration along a transect crossing the boundary

Compliance Notes

Non-stationarity handling is a modelling decision with a visible effect on published exposures, so it must be recorded in full: the regionalisation layer and its vintage, the site count per region, the pooled and per-region variogram parameters, the ratio thresholds used to declare non-stationarity, whether a trend was removed and how much variance it explained, and the per-region held-out error before and after the remedy.

Record also what you chose not to do. An analyst who finds a range ratio of 2.7 and decides the pooled model is adequate anyway has made a defensible call, but only if the ratio and the reasoning are in the record. A surface whose documentation is silent on stationarity cannot be distinguished from one where the question was never asked.