Spatial Epi Public Health GIS

Kernel Density & Spatial Relative Risk Surfaces

This guide is part of the Disease Clustering & Spatial Statistical Modeling section. It covers continuous surfaces estimated from point data — the method to use when cases are geocoded to addresses and the question is where risk is elevated, without first forcing the points into administrative areas.

Concept & Epidemiological Alignment

A kernel density estimate places a small smooth bump on every case point and sums them, producing a continuous surface of case intensity. It answers where cases are concentrated, and on its own it is almost useless for epidemiology, because cases are concentrated where people are. A density map of any disease in a state looks like a population map of that state.

The epidemiologically meaningful quantity is the ratio of two densities: the density of cases divided by the density of the population at risk. Where the population at risk is represented by controls — a random sample of the population, or cases of an unrelated condition — the log of that ratio is the spatial log relative risk surface, ρ(s)=log⁡{f(s)/g(s)}\rho(s) = \log\{f(s)/g(s)\}, where ff is the case density and gg the control density. It is zero where cases are distributed like controls and positive where cases are over-represented.

This approach sits between the areal methods and the point-pattern methods elsewhere in the section. Unlike Global & Local Moran’s I Implementation, it needs no areal units and so avoids the modifiable areal unit problem. Unlike K-Function & Point Pattern Analysis, which asks whether clustering exists at some scale, it asks where — and produces a map a program can act on.

Its assumptions are specific. The controls must represent the population at risk, so their spatial distribution matches the at-risk population; a control series drawn from one hospital’s catchment does not. The bandwidth must suit the scale of the effect; the ratio is exquisitely sensitive to it. And the surface near the study-area boundary is biased unless edge-corrected, because the kernel’s mass partly falls outside the region.

Surveillance use differs from research use in one important respect. A research study chooses its controls once and analyses a fixed dataset; a surveillance program re-estimates the surface as cases accrue, often weekly. That makes the control series a standing asset rather than a one-off sample — it must be refreshed as the population changes, geocoded with the same pipeline as incoming cases, and versioned so that a surface from March can be reproduced in December. Programs that treat the control set as disposable end up comparing surfaces built against different denominators and reading the difference as a change in risk.

Method Choices

Choice Options Consequence
Kernel Gaussian, quartic Shape matters little; bandwidth matters a lot
Bandwidth Fixed, adaptive, jointly optimal Fixed over-smooths cities or under-smooths countryside
Denominator Controls, population raster Controls allow significance testing; rasters give a rate-like surface
Edge correction None, uniform, Jones–Diggle Uncorrected surfaces dip at the boundary
Significance Asymptotic, Monte Carlo tolerance contours Monte Carlo is more robust with few cases
Density Shows Where Cases Are; the Ratio Shows Where Risk Is Three panels along a transect. The case density peaks in the city centre and at a smaller town, closely following population. The control density follows the same population pattern. Their log ratio is flat near zero across most of the transect, including the city centre, and rises only at one location near the edge of the city where cases are over-represented relative to controls. A case density map alone would have highlighted the city centre, which has no elevated risk at all. The case peak is the city; the risk peak is somewhere else case density control density log relative risk elevated risk city centre: many cases, no excess

Spatial Data Prerequisites

  • Geometry: case and control points as geocoded coordinates with match quality; a study-region polygon for edge correction.
  • Coordinate reference system: a projected metric system, since bandwidths are distances.
  • Sample size: the ratio is unstable with fewer than about fifty cases; below that, prefer an areal method or a global test.
  • Control design: controls sampled from the same population and period as cases, with the same geocoding process, so geocoding error affects both alike.
  • Positional accuracy: a bandwidth smaller than the geocode positional error resolves noise, not risk.

Production Implementation

The routine below computes case and control kernel densities on a grid with a shared bandwidth, applies a uniform edge correction, forms the log relative risk, and builds Monte Carlo tolerance contours by random relabelling.

# Spatial log relative risk with edge correction and Monte Carlo tolerance.
# Pinned: numpy==1.26.4, scipy==1.13.1, geopandas==0.14.4, shapely==2.0.4
import logging

import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree
from shapely.geometry import Point

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

SEED = 20260918


def grid_over(region: gpd.GeoDataFrame, step: float) -> np.ndarray:
    minx, miny, maxx, maxy = region.total_bounds
    xs = np.arange(minx + step / 2, maxx, step)
    ys = np.arange(miny + step / 2, maxy, step)
    xx, yy = np.meshgrid(xs, ys)
    g = np.c_[xx.ravel(), yy.ravel()]
    poly = region.union_all()
    inside = np.array([poly.contains(Point(p)) for p in g])
    return g[inside]


def gaussian_kde(points: np.ndarray, grid: np.ndarray, h: float) -> np.ndarray:
    """Normalised Gaussian kernel density; neighbours beyond 4h ignored for speed."""
    tree = cKDTree(points)
    dens = np.zeros(len(grid))
    for i, idx in enumerate(tree.query_ball_point(grid, 4 * h)):
        if idx:
            d2 = ((points[idx] - grid[i]) ** 2).sum(1)
            dens[i] = np.exp(-d2 / (2 * h * h)).sum()
    return dens / (len(points) * 2 * np.pi * h * h)


def edge_weights(grid: np.ndarray, region, h: float, n_mc: int = 64, seed: int = SEED):
    """Share of each grid point's kernel mass inside the region (uniform correction)."""
    rng = np.random.default_rng(seed)
    offs = rng.normal(0, h, size=(n_mc, 2))
    poly = region.union_all()
    w = np.array([np.mean([poly.contains(Point(p + o)) for o in offs]) for p in grid])
    return np.clip(w, 0.05, 1.0)


def log_relative_risk(cases: np.ndarray, controls: np.ndarray, grid: np.ndarray,
                      region, h: float, n_sim: int = 199, seed: int = SEED) -> dict:
    ew = edge_weights(grid, region, h)
    f = gaussian_kde(cases, grid, h) / ew
    g = gaussian_kde(controls, grid, h) / ew
    floor = 1e-12
    rho = np.log(np.maximum(f, floor)) - np.log(np.maximum(g, floor))

    # Monte Carlo: relabel pooled points at random, recompute, rank the observed.
    rng = np.random.default_rng(seed)
    pooled = np.vstack([cases, controls])
    n1 = len(cases)
    exceed = np.zeros(len(grid))
    for _ in range(n_sim):
        perm = rng.permutation(len(pooled))
        fs = gaussian_kde(pooled[perm[:n1]], grid, h) / ew
        gs = gaussian_kde(pooled[perm[n1:]], grid, h) / ew
        rs = np.log(np.maximum(fs, floor)) - np.log(np.maximum(gs, floor))
        exceed += rs >= rho
    p_upper = (exceed + 1) / (n_sim + 1)
    log.info("relative risk: %d grid cells, bandwidth %.0f m, %d cells with p<0.05",
             len(grid), h, int((p_upper < 0.05).sum()))
    return {"rho": rho, "p_upper": p_upper, "edge_weight": ew}

Using one bandwidth for both densities is a deliberate choice: different bandwidths for cases and controls make the ratio reflect the difference in smoothing rather than the difference in distribution. The Monte Carlo relabelling tests the null hypothesis that cases and controls come from the same spatial distribution, and the resulting cell-wise p-values define tolerance contours — the regions where the surface is significantly elevated — which are what a map should outline.

The edge correction divides each density by the fraction of the kernel’s mass that falls inside the region. Without it both densities dip near the boundary; the dips partly cancel in the ratio, but not exactly, and the uncorrected surface shows spurious gradients along coastlines and state lines.

Bandwidth Decides What the Surface Can Show The same case and control data smoothed at three bandwidths. At a very small bandwidth the log relative risk surface is fragmented into many small spikes around individual cases, most of them noise. At an intermediate bandwidth a single coherent area of elevated risk appears with a clear tolerance contour. At a very large bandwidth the elevation is spread across half the study area and diluted below significance. The appropriate bandwidth depends on the spatial scale of the effect, which is why selection is a topic of its own. Too small is noise, too large is nothing h = 200 m spikes at single cases h = 1.5 km one area, significant h = 10 km diluted, not significant dashed outline: tolerance contour at p < 0.05

Parameter Selection & Tuning

Bandwidth is the dominant choice. Data-driven selectors that jointly minimise the error of the ratio exist and are the defensible default; the details and their failure modes are in Choosing a Kernel Bandwidth for Case Density Maps.

Fixed versus adaptive. A fixed bandwidth suits a study area of roughly uniform density; where density varies strongly, an adaptive bandwidth that shrinks in cities and widens in the countryside is better, as covered in Adaptive Bandwidth KDE for Uneven Population Density.

Grid resolution should be a fraction of the bandwidth — a quarter or less — so the surface is sampled finely enough to draw smooth contours. Finer adds computation without information.

Simulation count for tolerance contours: 199 is a practical minimum for a 0.05 threshold; 999 gives more stable contours and is worth it for a published map.

Control ratio. Two to four controls per case is usually adequate; beyond that the control density is already well estimated and more controls add little.

Kernel shape. The Gaussian and quartic kernels give nearly identical ratios at matched effective bandwidths; the practical difference is that the quartic kernel has finite support, so each grid cell depends only on nearby points, which makes computation faster and keeps a distant outlier from contributing. Choose on computational grounds and match bandwidths when comparing with published work that used the other shape — a quartic bandwidth corresponds to a Gaussian bandwidth roughly 2.5 times smaller.

Temporal window. A surface pooled over five years mixes periods in which risk may have moved. For surveillance, estimate on a rolling window long enough to accumulate the minimum case count and short enough to reflect current risk, and state the window alongside the bandwidth, since both control the effective resolution of the map.

Edge Cases & Failure Modes

Controls that do not represent the population. A control series from a single clinic or a convenience sample reproduces that source’s catchment, and the ratio then maps the difference between case and control catchments. Validate the control distribution against census population before use.

Differential geocoding. If cases geocode at rooftop level and controls at street level, the smaller-scale structure of the two densities differs for a reason unrelated to risk. Use the same geocoder and the same quality filter for both.

Too few cases. With twenty cases the relative risk surface is dominated by the kernel shape around each point. Report a global test instead of a map.

Multiple testing across the grid. Cell-wise p-values are not independent and a 0.05 threshold over thousands of cells will outline some noise. Tolerance contours are exploratory; the significance procedure is covered in Testing Significance of a Relative Risk Surface.

Clustered controls. Controls sampled by household or by clinic arrive in clumps, and the control density then has spikes that the case density lacks, producing apparent troughs of risk around each clump. Sample controls independently at the individual level, or thin clumped controls to one per location before estimation.

Changing boundaries of the study region. Edge correction depends on the region polygon. A surface recomputed after the region’s boundary is edited — a county added, a coastline refined — changes near the edges for that reason alone. Hold the region polygon fixed per vintage and record it with the surface.

Disclosure. A kernel surface at a small bandwidth reveals the locations of individual cases — each spike is a person. Publication bandwidths and output resolution must satisfy disclosure review regardless of what the analysis used internally.

A Narrow Kernel Is a Map of Addresses Two versions of a case density surface over the same small area. At an analytic bandwidth of one hundred metres each case appears as its own distinct peak, and the peak locations are the geocoded addresses to within the kernel's width, which is a disclosure. At a publication bandwidth chosen by disclosure review, the peaks merge into a smooth surface from which no individual location can be recovered. The analysis may use the narrow kernel internally; the published product may not. Analytic bandwidth and publication bandwidth are different decisions h = 100 m: one peak per case peaks mark addresses publication bandwidth no individual recoverable publication bandwidth is set by disclosure review, not by the analysis

Compliance & Audit Controls

Record the case and control definitions and sampling, the geocoder and quality filter applied to both, the kernel, the bandwidth and how it was selected, the edge-correction method, the grid resolution, the simulation count and seed, and the significance threshold for tolerance contours.

Kernel surfaces are among the most disclosive outputs in spatial epidemiology because they are derived directly from point locations. Treat the internal analytic surface as protected data. The published surface should use a bandwidth and grid resolution approved under Privacy-Preserving Spatial Analytics and should never be released alongside the point count per cell.

Production Implementation Checklist

Frequently Asked Questions

Why not just map case density? Because case density mostly reproduces population density. A cluster of cases in a city centre is expected; what matters is whether cases are over-represented relative to the population there, which only the ratio of case to control densities shows.

Can I use a population raster instead of controls? Yes, which gives a rate-like surface. You lose the relabelling test, since there are no control points to permute, and must use a different significance approach such as simulating cases from the population surface under the null.

What if cases and controls are geocoded differently? Fix it before analysis. Differential positional error changes the fine-scale structure of one density and not the other, producing spurious relative risk at scales below the error.

How do I report a relative risk surface? Map the log relative risk with tolerance contours outlining significant areas, state the bandwidth and edge correction, and describe the contours as exploratory unless a global test has confirmed a departure from the null.

Is a kernel surface safe to publish? Only at a bandwidth and resolution approved by disclosure review. At small bandwidths each case appears as a distinct peak at its address.