Spatial Epi Public Health GIS

Choosing a Kernel Bandwidth for Case Density Maps

This guide, part of Kernel Density & Spatial Relative Risk Surfaces, deals with the single parameter that dominates every kernel-based map: the bandwidth.

Problem Context & Constraints

The kernel shape barely matters; the bandwidth decides almost everything. Too small and every case becomes its own peak, the map is noise, and the individual addresses are visible. Too large and a real concentration is spread across a county and diluted below detection. Between those extremes the map changes character continuously, and two analysts choosing bandwidths by eye from the same data will publish different maps and draw different conclusions.

Several selectors exist and they disagree, predictably. Rule-of-thumb selectors derive a bandwidth from the spread of the points assuming a roughly Gaussian distribution; case data are never Gaussian, so these over-smooth anything with local structure. Least-squares cross-validation chooses the bandwidth that best predicts left-out points; it adapts to local structure but is notoriously unstable, often selecting a very small value when there are tied or near-duplicate points. Oversmoothing selectors give an upper bound — the largest bandwidth consistent with the data’s spread — useful as a sanity ceiling rather than a choice.

For relative risk surfaces there is an extra subtlety: the case and control densities share one bandwidth, and the bandwidth that is optimal for each density separately is not optimal for their ratio. A jointly optimal bandwidth, chosen to minimise the error of the log ratio, is the right target.

Two floors bound any choice from below. The geocode positional error: a bandwidth smaller than the uncertainty in the point locations resolves noise. And disclosure: at small bandwidths each case is a visible peak, so the published bandwidth is constrained by review regardless of what the analysis prefers.

Four Selectors, Four Answers, Two Floors A horizontal bandwidth axis with the choices of four selectors on the same case data marked. Least-squares cross-validation picks a very small bandwidth driven by near-duplicate points. The jointly optimal relative-risk selector picks an intermediate value. A rule-of-thumb Gaussian selector picks a larger value that over-smooths local structure. The oversmoothing bound sits at the far right as a ceiling. Two shaded floors on the left mark the geocode positional error and the disclosure minimum; the cross-validation choice falls below both and is inadmissible. Selectors disagree; the floors rule some of them out positional error disclosure floor least-squares CV joint relative-risk Gaussian rule of thumb oversmoothing ceiling 200 m 12 km bandwidth on a log scale for one case-control dataset report the selector, the value, and which floor was binding

Prerequisites

  • python 3.11
  • numpy 1.26.4, scipy 1.13.1
  • Case and control coordinates in a projected metric coordinate reference system
  • The geocode positional error distribution from the geocoding stage
  • The disclosure review’s minimum publication bandwidth

Step-by-Step Solution

Compute a candidate from each selector, choose the jointly optimal value for relative risk, apply the floors, and check the result’s stability by resampling.

# Bandwidth selection for kernel density and relative risk, with floors.
# Pinned: numpy==1.26.4, scipy==1.13.1
import logging

import numpy as np
from scipy.optimize import minimize_scalar
from scipy.spatial.distance import pdist

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

SEED = 20260918


def scott(points: np.ndarray) -> float:
    """Gaussian rule of thumb for 2-D data."""
    n = len(points)
    sigma = np.sqrt(np.mean(points.var(axis=0, ddof=1)))
    return float(sigma * n ** (-1.0 / 6.0))


def oversmooth(points: np.ndarray) -> float:
    """Upper bound on a sensible fixed bandwidth (Terrell's oversmoothing rule)."""
    n = len(points)
    sigma = np.sqrt(np.mean(points.var(axis=0, ddof=1)))
    return float(1.08 * sigma * n ** (-1.0 / 6.0) * 1.7)


def lscv(points: np.ndarray, lo: float, hi: float) -> float:
    """Least-squares cross-validation for a 2-D Gaussian kernel."""
    d2 = pdist(points, "sqeuclidean")
    n = len(points)

    def score(h):
        a = np.exp(-d2 / (4 * h * h)).sum() * 2 / (4 * np.pi * h * h * n * n)
        b = np.exp(-d2 / (2 * h * h)).sum() * 2 / (2 * np.pi * h * h * n * (n - 1))
        return 1 / (4 * np.pi * h * h * n) + a - 2 * b

    r = minimize_scalar(score, bounds=(lo, hi), method="bounded")
    return float(r.x)


def joint_rr(cases: np.ndarray, controls: np.ndarray, lo: float, hi: float) -> float:
    """Shared bandwidth minimising a cross-validated error of the log ratio.

    Leave-one-out log-density of each point under its own group versus the other,
    averaged — a practical proxy for the integrated error of the ratio.
    """
    def loo_log_density(pts, h, other=None):
        tgt = pts if other is None else other
        d2 = ((tgt[:, None, :] - pts[None, :, :]) ** 2).sum(-1)
        k = np.exp(-d2 / (2 * h * h))
        if other is None:
            np.fill_diagonal(k, 0.0)
            denom = len(pts) - 1
        else:
            denom = len(pts)
        return np.log(np.maximum(k.sum(1) / (denom * 2 * np.pi * h * h), 1e-300))

    def score(h):
        lc = loo_log_density(cases, h).mean()
        lk = loo_log_density(controls, h).mean()
        return -(lc + lk)

    r = minimize_scalar(score, bounds=(lo, hi), method="bounded")
    return float(r.x)


def choose(cases, controls, pos_error_p90: float, disclosure_min: float) -> dict:
    pooled = np.vstack([cases, controls])
    cand = {"scott": scott(pooled), "oversmooth": oversmooth(pooled)}
    lo, hi = max(10.0, pos_error_p90 / 4), cand["oversmooth"]
    cand["lscv"] = lscv(cases, lo, hi)
    cand["joint_rr"] = joint_rr(cases, controls, lo, hi)
    floor = max(pos_error_p90, disclosure_min)
    chosen = max(cand["joint_rr"], floor)
    binding = ("selector" if chosen == cand["joint_rr"] else
               "positional error" if floor == pos_error_p90 else "disclosure")
    log.info("candidates (m): %s; floor %.0f; chosen %.0f (binding: %s)",
             {k: round(v) for k, v in cand.items()}, floor, chosen, binding)
    return {"candidates": cand, "chosen_m": chosen, "binding": binding}


def stability(cases, controls, selector, n_boot: int = 50, seed: int = SEED) -> dict:
    """Does the selected bandwidth survive resampling of the data?"""
    rng = np.random.default_rng(seed)
    hs = []
    for _ in range(n_boot):
        c = cases[rng.integers(0, len(cases), len(cases))]
        k = controls[rng.integers(0, len(controls), len(controls))]
        hs.append(selector(c, k))
    hs = np.asarray(hs)
    out = {"median": float(np.median(hs)), "cv": float(hs.std() / hs.mean())}
    log.info("bandwidth stability: median %.0f m, CV %.2f over %d resamples",
             out["median"], out["cv"], n_boot)
    return out

The binding field in the result is what should appear in the methods text. It says whether the bandwidth reflects the data (the selector chose it), the data quality (positional error forced it up), or the release constraint (disclosure forced it up). Those are three different stories about the resulting map, and a reader cannot tell them apart from the number alone.

The stability check matters because data-driven selectors can be fragile. A bandwidth whose value swings by half across bootstrap resamples is not well determined by the data, and the map should be presented as one of several plausible smoothings rather than as the answer.

Is the Bandwidth Determined by the Data? Bootstrap distributions of the selected bandwidth for two datasets. For a dataset with four hundred cases the distribution is narrow, with a coefficient of variation near a tenth, so the selected value is well determined. For a dataset with sixty cases the distribution is wide and bimodal, with some resamples choosing very small bandwidths driven by near-duplicate points, so no single value is supported and the map should be shown at more than one bandwidth. A bandwidth that jumps between resamples is not an answer 400 cases: CV 0.1 well determined 60 cases: bimodal, CV 0.6 show several smoothings selected bandwidth across bootstrap resamples

Validation & Edge Cases

Duplicate coordinates. Multiple cases at one address — a household, an apartment building — create zero distances that pull cross-validation toward a tiny bandwidth. Jitter duplicates by a fraction of the positional error for selection only, or use the joint selector, which is less sensitive.

Very uneven density. No fixed bandwidth suits both a dense city and the surrounding countryside. If the selected value is plainly too large for the city and too small for the countryside, move to an adaptive bandwidth.

Different case and control counts. A selector applied to controls alone, when controls outnumber cases fourfold, picks a smaller bandwidth than the cases can support. Always select on the ratio or on the pooled data, never on controls alone.

Selecting after looking at maps. Trying several bandwidths and keeping the one whose map “shows the cluster” is the kernel equivalent of choosing a buffer radius by its result. Select by the declared method first; show alternatives as sensitivity.

Anisotropic spread. Cases strung along a river valley or a transit corridor spread much more in one direction than the other, and a circular kernel with one bandwidth over-smooths across the corridor while under-smoothing along it. Where the pattern is clearly directional, consider a bandwidth matrix rather than a scalar, or rotate coordinates to the corridor axis before selection, and say so in the methods.

Bandwidth across time. A weekly surface whose bandwidth is re-selected every week changes smoothness with the case count. Fix the bandwidth for a surveillance series, or select it on a pooled baseline period.

INFO kde.bandwidth candidates (m): {'scott': 4210, 'oversmooth': 7730, 'lscv': 180, 'joint_rr': 1460}; floor 350; chosen 1460 (binding: selector)
INFO kde.bandwidth bandwidth stability: median 1420 m, CV 0.14 over 50 resamples
Which Constraint Chose the Bandwidth Three cases for the final bandwidth. When the selector's value is above both floors, the map reflects the data's own structure. When the positional error floor binds, the map is limited by geocode quality, and finer structure may exist that the data cannot show. When the disclosure floor binds, the map is limited by the release constraint, and the analytic surface used internally may show more than the published one. The methods text should say which case applies. Say which constraint set the number selector binds map reflects the data's own structure positional error binds limited by geocode quality; finer structure unseen disclosure binds limited by release rules; internal surface is finer the same published number means three different things

Compliance Notes

Record every candidate bandwidth with its selector, the positional error quantile used as a floor, the disclosure minimum, the chosen value and which constraint was binding, and the bootstrap stability result. For a surveillance series, record the fixed bandwidth and the baseline period it was selected on.

The disclosure floor deserves its own entry in the release record: it is a decision by the disclosure reviewer, not the analyst, and changing it changes what the published map can reveal. Where the analytic bandwidth is smaller than the published one, the internal surface is a protected artifact and must not leak into presentations or appendices.