Spatial Epi Public Health GIS

Choosing an IDW Power for Sparse Monitor Networks

This guide, part of Interpolating Environmental Exposure Surfaces, addresses the situation every small health department eventually faces: too few monitors to fit a variogram, a deadline, and a default power of 2 sitting unexamined in the call.

Problem Context & Constraints

Inverse distance weighting predicts at a location as a weighted average of observations, with weights proportional to 1/dp1/d^{p}. The exponent pp is the entire model. It is not a tuning detail; it decides whether the surface is nearly the global mean or a mosaic of nearest-neighbour patches, and the two produce different health effect estimates from the same monitors.

The default of p=2p = 2 has no physical justification for environmental concentrations. It is inherited from an early implementation and survives because it is the default. For a regionally smooth pollutant, p=2p = 2 is usually too aggressive: it concentrates weight on the closest site, reproduces its value nearly exactly, and generates the characteristic bullseye pattern. For a pollutant with genuine short-range structure, p=2p = 2 may be too flat, smearing a real gradient across a county.

Two constraints shape the fix. The exponent must be selected out of sample, because in-sample error falls monotonically as pp rises — at the limit, the surface passes exactly through every observation and the training error is zero while the prediction error is at its worst. And the search must be bounded by behaviour, not taste: below roughly p=0.5p = 0.5 the surface is indistinguishable from the mean, above roughly p=6p = 6 it is indistinguishable from nearest-neighbour assignment, and anywhere in between the difference is measurable.

What the Exponent Actually Controls Three profiles across the same transect between two monitors with different values. At a low exponent the profile is nearly flat at the average of the two sites, so local contrast disappears. At a moderate exponent the profile is a smooth gradient between the two values, which is the behaviour most regional pollutants call for. At a high exponent the profile is almost a step, holding each monitor's own value until it flips abruptly at the midpoint, which is nearest-neighbour assignment in disguise. The three panels share the same axis so the transition is directly comparable. The same two monitors, three exponents, three different exposures p = 0.5 — near the mean contrast erased p = 2 — smooth gradient usually the right shape, rarely the right value p = 6 — patches nearest monitor with extra steps transect between two monitors · vertical axis is concentration In-sample error prefers the right-hand panel; out-of-sample error almost never does which is why the exponent must be chosen by holding sites out

Prerequisites

  • python 3.11
  • numpy 1.26.4, geopandas 0.14.4
  • Monitor sites as points with a stable site identifier, in a projected coordinate reference system in metres
  • Site-level values already reduced to one averaging window with a completeness rule applied
  • At least eight sites — below that, cross-validation over a power grid is noise and a nearest-monitor assignment with an explicit caveat is more honest

Step-by-Step Solution

Sweep the exponent over a bounded grid, score each value by leave-one-site-out error, and select the minimiser — with a tie-break toward the smaller exponent, because a flatter surface fails more gracefully than a patchy one when the network is sparse.

# Select an IDW exponent by leave-one-site-out cross-validation.
# 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.idw")

POWERS = np.round(np.arange(0.5, 6.01, 0.25), 2)


def idw_predict(train_xy, train_z, target_xy, power, eps=1.0):
    """IDW prediction. eps (metres) avoids a division by zero at a co-located site."""
    d = np.sqrt(((target_xy[:, None, :] - train_xy[None, :, :]) ** 2).sum(-1))
    d = np.maximum(d, eps)
    w = d ** (-power)
    return (w * train_z[None, :]).sum(1) / w.sum(1)


def select_power(sites: gpd.GeoDataFrame, value_col: str,
                 powers=POWERS) -> dict:
    """Leave-one-site-out sweep over the exponent grid."""
    if sites.crs is None or sites.crs.is_geographic:
        raise ValueError("project to metres before weighting by distance")
    sites = sites.sort_values("site_id").reset_index(drop=True)   # deterministic
    xy = np.c_[sites.geometry.x.values, sites.geometry.y.values]
    z = sites[value_col].to_numpy(float)
    n = len(sites)
    if n < 8:
        raise ValueError(f"{n} sites is too few to select a power by cross-validation")

    scores = {}
    for p in powers:
        resid = np.empty(n)
        for i in range(n):
            m = np.ones(n, dtype=bool)
            m[i] = False
            resid[i] = z[i] - idw_predict(xy[m], z[m], xy[i:i + 1], p)[0]
        scores[float(p)] = {
            "rmse": float(np.sqrt((resid ** 2).mean())),
            "bias": float(resid.mean()),
            "mae": float(np.abs(resid).mean()),
        }

    best_rmse = min(s["rmse"] for s in scores.values())
    # Tie-break toward the flatter surface: accept any power within 2% of the best.
    admissible = [p for p, s in scores.items() if s["rmse"] <= best_rmse * 1.02]
    chosen = min(admissible)
    flat_rmse = float(np.sqrt(((z - z.mean()) ** 2).mean()))
    log.info("idw power chosen=%.2f rmse=%.4f (global-mean rmse=%.4f) n_sites=%d",
             chosen, scores[chosen]["rmse"], flat_rmse, n)
    if scores[chosen]["rmse"] >= flat_rmse:
        log.warning("IDW does not beat the global mean — the network carries no "
                    "usable spatial information at this scale")
    return {"power": chosen, "scores": scores, "global_mean_rmse": flat_rmse,
            "n_sites": n}

The comparison against the global-mean error is the part most implementations omit and the part a reviewer will ask about first. If leave-one-site-out error for the best exponent is no better than simply assigning everyone the state average, the interpolation is decoration: it varies across the map, but none of that variation is supported by the data. Saying so in the log is what stops the surface from being used as though it meant something.

The tie-break rule matters for sparse networks specifically. The error curve over pp is usually shallow near its minimum, and picking the exact argmin can land on a high exponent whose advantage is within noise. Taking the smallest exponent within a small tolerance of the best error buys robustness at negligible cost.

Reading the Exponent Sweep A curve of leave-one-site-out root mean squared error against the inverse distance weighting exponent. The curve falls from a high value at a very low exponent, reaches a broad shallow minimum, then rises again as the exponent approaches nearest-neighbour behaviour. A horizontal dashed line marks the error obtained by assigning every site the global mean; the curve must fall clearly below it for the interpolation to be worth anything. A shaded tolerance band two per cent above the minimum shows the range of exponents that are statistically indistinguishable, and the chosen exponent is the smallest one inside that band. A shallow minimum means several exponents are equally defensible error of assigning the global mean within 2% of the minimum chosen: smallest in band 0.5 1.75 3.5 6 inverse distance weighting exponent If the curve never crosses the dashed line, do not ship the surface held-out error

Validation & Edge Cases

Co-located or near-co-located sites. Two monitors a few metres apart make the weight at one of them blow up, and the epsilon floor in the code is what keeps the prediction finite. Check the nearest-neighbour distance distribution first; a pair under about 100 m should usually be averaged into one pseudo-site instead.

A monitor with an extreme value. Inverse distance weighting has no mechanism to discount an outlier — a site reading three times its neighbours dominates a wide area at any exponent above about 1.5. Screen for outliers before the sweep, and if one is genuine rather than erroneous, that is itself a signal that the field is not smooth and the method is wrong.

A directional network. Sites strung along a highway or a valley floor produce a surface that is well informed along the string and pure extrapolation perpendicular to it, and inverse distance weighting cannot express that asymmetry. The predictions off-axis will look no less confident than the on-axis ones, because there is no variance to tell you otherwise.

The sweep picks the top of the grid. An exponent selected at the boundary of the search — 6 in this implementation — means the data prefer nearest-neighbour assignment. Do not extend the grid; use nearest monitor explicitly, which at least makes the assumption visible to a reader.

INFO exposure.idw idw power chosen=1.75 rmse=1.3820 (global-mean rmse=2.4471) n_sites=23

A sweep run separately per year. Refitting the exponent independently for each year of a multi-year cohort produces a surface whose smoothness changes between years, and a change in smoothness moves every exposure value even when the underlying field did not move. Select one exponent on pooled data across the study period, confirm it is inside the tolerance band in each individual year, and refit only when a year falls clearly outside it. Record which years were checked, because a reviewer comparing your 2019 and 2024 exposures needs to know that the smoothing was held constant rather than re-optimised.

A grid resolution finer than the network can support. Rendering a 1.75-exponent surface onto a 30 m grid does not add information; it adds pixels. The visual precision of the output invites readers to interpret block-level differences that come entirely from the interpolation geometry. Choose a cell size on the order of a fraction of the typical inter-site distance, and state that number next to the map.

No variance, ever. This is the defining limitation. Inverse distance weighting gives a number with no stated precision, so the honest downstream treatment is to use it for ranking and screening rather than for a per-unit dose-response estimate, and to say in the methods that no prediction uncertainty is available. Where the analysis needs that uncertainty, the network must be large enough for the geostatistical route described in the parent guide.

The one-exponent-per-cohort rule is easiest to see as a timeline of what refitting each year does to an address whose surroundings never changed:

Refitting Each Year Moves Exposures That Did Not Move Assigned exposure at a single residence over five years. With the exponent refitted independently each year, the value jumps up and down because the smoothness of the surface changed between fits, not because the field changed. With one exponent selected on pooled data and checked against each year, the series is smooth and its movement reflects changes in the monitor values only. Hold the exponent fixed across the cohort period exponent refitted yearly one pooled exponent 2020 2022 2024 assigned exposure at one address whose surroundings did not change exposure

Compliance Notes

Record the exponent, the search grid, the selection rule including the tolerance used for the tie-break, the leave-one-site-out error and bias, the global-mean reference error, the epsilon distance floor and the site count. The tolerance rule in particular must be written down: it is a discretionary choice that changes the shipped surface, and a reviewer who cannot see it cannot reproduce your result even with the same data.

Two statements belong in the methods section of anything built on this surface. First, that the exponent was selected by cross-validation rather than left at a default. Second, that the surface carries no prediction variance, with the consequence that measurement-error correction of the kind described in Exposure Measurement Error & Uncertainty is not available for it.