Spatial Epi Public Health GIS

Adaptive Bandwidth KDE for Uneven Population Density

This guide, part of Kernel Density & Spatial Relative Risk Surfaces, handles the study area that no single bandwidth fits: a state with a dense metropolitan core and a sparse rural hinterland.

Problem Context & Constraints

A fixed bandwidth imposes one scale of smoothing everywhere. In the city, where points are hundreds of metres apart, a bandwidth chosen to stabilise the rural estimate blurs neighbourhood-scale variation into a district average. In the countryside, where points are kilometres apart, a bandwidth chosen to resolve the city leaves isolated spikes around individual cases with nothing between them. Any fixed value is wrong for one of the two, and the usual compromise is wrong for both.

Adaptive kernels give each point its own bandwidth, small where points are dense and large where they are sparse. The classic scheme sets each point’s bandwidth inversely proportional to the square root of a pilot density estimate at that point, scaled by a global bandwidth. In the city, where the pilot density is high, bandwidths shrink; in the countryside they grow. The resulting surface has roughly constant precision across the study area rather than constant smoothing.

Three constraints shape a sound implementation. The pilot density must itself be stable, so it is estimated with a fixed bandwidth that is deliberately generous. Bandwidths must be bounded, because the inverse-square-root rule sends bandwidths to extreme values at the edges of the point cloud. And for a relative risk surface, cases and controls must use one adaptive scheme derived from the pooled data, or the ratio reflects the difference in adaptation rather than risk.

Fixed Footprints Versus Adaptive Footprints Two strips crossing from a dense city on the left to sparse countryside on the right, with points and their kernel footprints drawn as circles. Under a fixed bandwidth every circle is the same size: in the city they overlap heavily and blur local detail, and in the countryside they barely touch, leaving isolated spikes. Under an adaptive bandwidth the circles are small in the city and large in the countryside, overlapping to a similar degree everywhere, which gives the surface roughly uniform precision. Adapt the kernel so every point sees a similar number of neighbours fixed city: blurred countryside: isolated spikes adaptive city: detail kept countryside: stable

Prerequisites

  • python 3.11
  • numpy 1.26.4, scipy 1.13.1
  • Case and control coordinates in a projected metric coordinate reference system
  • A global bandwidth from Choosing a Kernel Bandwidth for Case Density Maps
  • Lower and upper bounds for per-point bandwidths — the lower at least the geocode positional error

Step-by-Step Solution

Estimate a pilot density on the pooled points, derive per-point bandwidths by the inverse-square-root rule, clip them to the bounds, and evaluate both densities with those same bandwidths.

# Abramson adaptive kernel density shared by cases and controls.
# Pinned: numpy==1.26.4, scipy==1.13.1
import logging

import numpy as np
from scipy.spatial import cKDTree

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


def fixed_kde_at(points: np.ndarray, at: np.ndarray, h: float) -> np.ndarray:
    tree = cKDTree(points)
    out = np.zeros(len(at))
    for i, idx in enumerate(tree.query_ball_point(at, 4 * h)):
        if idx:
            d2 = ((points[idx] - at[i]) ** 2).sum(1)
            out[i] = np.exp(-d2 / (2 * h * h)).sum()
    return out / (len(points) * 2 * np.pi * h * h)


def abramson_bandwidths(pooled: np.ndarray, h0: float, h_min: float, h_max: float,
                        pilot_factor: float = 1.5) -> np.ndarray:
    """h_i = h0 * (pilot(x_i) / geometric_mean(pilot)) ** -0.5, clipped."""
    pilot = fixed_kde_at(pooled, pooled, h0 * pilot_factor)   # generous pilot
    pilot = np.maximum(pilot, pilot[pilot > 0].min())
    gmean = np.exp(np.mean(np.log(pilot)))
    h = h0 * (pilot / gmean) ** -0.5
    clipped_lo = int((h < h_min).sum())
    clipped_hi = int((h > h_max).sum())
    h = np.clip(h, h_min, h_max)
    log.info("adaptive bandwidths: median %.0f m, range %.0f-%.0f m, clipped low=%d high=%d",
             np.median(h), h.min(), h.max(), clipped_lo, clipped_hi)
    return h


def adaptive_kde(points: np.ndarray, h_points: np.ndarray, grid: np.ndarray) -> np.ndarray:
    """Sample-point adaptive estimator: each point carries its own bandwidth."""
    dens = np.zeros(len(grid))
    tree = cKDTree(grid)
    for p, h in zip(points, h_points):
        idx = tree.query_ball_point(p, 4 * h)
        if idx:
            d2 = ((grid[idx] - p) ** 2).sum(1)
            dens[idx] += np.exp(-d2 / (2 * h * h)) / (2 * np.pi * h * h)
    return dens / len(points)


def adaptive_log_rr(cases, controls, grid, h0, h_min, h_max):
    pooled = np.vstack([cases, controls])
    h_all = abramson_bandwidths(pooled, h0, h_min, h_max)
    h_cases, h_ctrl = h_all[:len(cases)], h_all[len(cases):]
    f = adaptive_kde(cases, h_cases, grid)
    g = adaptive_kde(controls, h_ctrl, grid)
    rho = np.log(np.maximum(f, 1e-300)) - np.log(np.maximum(g, 1e-300))
    return rho, h_all

Deriving the bandwidths from the pooled points is what keeps the ratio honest. If cases and controls each derived their own adaptive bandwidths, a region where cases happen to be denser than controls would give cases smaller bandwidths there, sharpening their density relative to the controls’ and producing an apparent excess from the adaptation alone. One pooled pilot gives both densities identical smoothing at every location.

The clipping bounds are not optional refinements. Without a lower bound, points in the densest cluster receive bandwidths below the geocode positional error and the surface resolves noise; without an upper bound, the most isolated rural points receive bandwidths of tens of kilometres and smear a single case across a county.

The Adaptation Rule and Its Bounds Per-point bandwidth plotted against pilot density on logarithmic axes. The inverse square root rule is a straight descending line. A horizontal cap at the upper bound flattens it at low pilot densities, preventing isolated rural points from receiving enormous bandwidths. A horizontal floor at the lower bound flattens it at high densities, preventing points in the densest cluster from receiving bandwidths below the geocode positional error. The counts of points clipped at each bound are reported. The rule is a line; the bounds keep its ends sensible upper bound lower bound h ∝ pilot^-1/2 rural urban core pilot density at the point bandwidth

Validation & Edge Cases

Pilot too narrow. A pilot estimated with a small bandwidth is noisy, and its noise propagates into the adaptive bandwidths, producing a surface with more texture than the data support. Use a pilot bandwidth somewhat larger than the global one.

Many points clipped. If a large share of points hits either bound, the bounds rather than the data are shaping the surface. Revisit the global bandwidth or the bounds, and report the clipped counts.

Edge correction with variable bandwidths. The uniform edge correction from the parent guide assumes one bandwidth. With adaptive bandwidths, each point’s kernel extends a different distance outside the region; compute the correction per point rather than per grid cell.

Testing with relabelling. Monte Carlo relabelling must keep the pooled adaptive bandwidths attached to locations, not to labels. Relabel which points are cases while each point keeps its bandwidth, so the null distribution reflects only the labelling.

Adaptive versus fixed as a sensitivity check. The adaptive surface should not be the only one examined. Estimate the fixed-bandwidth surface at the global bandwidth as well, and compare their tolerance contours. A contour that appears under both is robust to the smoothing scheme; a contour that appears only under adaptive smoothing is usually a small dense cluster the fixed kernel blurred, which may be real, and should be examined with the case count inside it before it is reported.

Computational cost. The sample-point estimator evaluates each point’s kernel over the grid within four of its own bandwidths, so large rural bandwidths touch many grid cells. For state-wide grids, evaluate on a coarser grid in low-density regions, or tile the computation, rather than letting a few isolated points dominate runtime.

Controls much denser than cases. With four controls per case the pooled pilot is driven mostly by controls, which is appropriate — it represents the population — but means case clusters that are dense only in cases do not shrink their own bandwidths. That is the intended behaviour: adaptation should follow the population at risk, not the outcome.

Communicating variable smoothing. A reader assumes a map is equally smooth everywhere. State that smoothing varies with population density, and consider showing the bandwidth surface beside the risk surface.

INFO kde.adaptive adaptive bandwidths: median 1180 m, range 350-6000 m, clipped low=214 high=37

The clipped counts in that line are worth reading as a pair. Many points at the lower bound is expected in a dense city and usually fine. Many points at the upper bound means the rural cases are so sparse that even a six-kilometre kernel isolates them, and the rural part of the surface should be presented with that caveat or not at all.

Show the Bandwidth Beside the Risk Two maps of the same study area side by side. The left shows the adaptive log relative risk with a tolerance contour in the city. The right shows the bandwidth used at each location, small in the city and large in the countryside. Read together, they tell a viewer that the city contour is resolved at neighbourhood scale while any rural pattern is resolved only at several kilometres, which prevents over-reading rural detail that the smoothing cannot support. Variable smoothing should be visible, not implied log relative risk contour in the city bandwidth used ~400 m ~5 km fine in city, coarse outside publish the pair so rural detail is not over-read

Compliance Notes

Record the global bandwidth and how it was selected, the pilot bandwidth factor, the lower and upper bounds with their justification, the counts of points clipped at each bound, the edge-correction method, and confirmation that cases and controls shared one pooled adaptive scheme.

Adaptive kernels raise a specific disclosure point: bandwidths shrink where points are dense, and in a small cluster of cases the lower bound — not the global bandwidth — sets how precisely those cases can be located from the surface. The disclosure review should evaluate the published surface at the lower bound, since that is where it is most revealing.