Spatial Epi Public Health GIS

Estimating a Spatial Relative Risk Surface

This guide, part of Kernel Density & Spatial Relative Risk Surfaces, runs the whole method once, from a case list to a published map, and concentrates on the steps that decide whether the result means anything: the control series and the matching of the two point sets.

Problem Context & Constraints

The mathematics of a relative risk surface is two kernel densities and a logarithm. The epidemiology is almost entirely in the controls. The surface compares where cases are with where the population at risk is, and the control points are the only representation of that population the method sees. If the controls are drawn from the wrong population, or geocoded differently, the surface faithfully maps the difference — which is not risk.

Three control strategies are common in surveillance. Random residential addresses sampled in proportion to population from an address-point file are the cleanest representation of the at-risk population, as long as the address file covers the whole study area evenly. Cases of an unrelated condition from the same surveillance system share the cases’ reporting and geocoding pathway, which cancels those biases, but carry their own geography if the control condition has any spatial structure. Population-weighted random points within census blocks are easy to generate and ignore within-block distribution, which is acceptable at bandwidths well above block size.

Three constraints follow. Controls must match the at-risk population in time and, where the outcome is age-specific, in age. Both series must pass through the same geocoding and quality filter, so positional error affects them alike. And the case count must be adequate — below roughly fifty cases the surface is dominated by kernel shape and should not be mapped.

Where the Controls Come From Decides What the Map Shows Three control strategies. Random residential addresses drawn in proportion to population represent the at-risk population cleanly but depend on the address file's even coverage. Cases of an unrelated condition share the reporting and geocoding pathway of the cases, cancelling those biases, but import any geography the control condition has. Population-weighted random points in census blocks are simple and ignore within-block distribution, acceptable when the bandwidth is well above block size. Each carries a stated assumption. Every control strategy makes one assumption; state which Random addresses sampled by population assumes: address file covers evenly usually the best default Unrelated-condition cases same reporting pathway assumes: control condition has no geography cancels reporting bias Points in blocks weighted by block population assumes: bandwidth well above block size simple, coarse validate the chosen controls against census population before estimating anything

Prerequisites

  • python 3.11
  • numpy 1.26.4, scipy 1.13.1, geopandas 0.14.4, pandas 2.2.2
  • The relative risk routine from the parent guide
  • Geocoded cases with match type; an address-point file or block populations for controls
  • A study-region polygon in a projected metric coordinate reference system

Step-by-Step Solution

Draw controls, filter both series identically, check the controls against population, estimate the surface and tolerance contours, and export contours rather than raw peaks.

# End-to-end relative risk surface: controls, matching, estimation, export.
# Pinned: numpy==1.26.4, scipy==1.13.1, geopandas==0.14.4, pandas==2.2.2
import logging

import numpy as np
import pandas as pd
import geopandas as gpd

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

SEED = 20260918
GOOD_MATCH = {"rooftop", "parcel"}
MIN_CASES = 50


def draw_controls(addresses: gpd.GeoDataFrame, n: int, seed: int = SEED) -> gpd.GeoDataFrame:
    """Simple random sample of residential address points — each address equally likely,
    so the sample is proportional to where households are."""
    rng = np.random.default_rng(seed)
    idx = rng.choice(len(addresses), size=n, replace=False)
    ctrl = addresses.iloc[np.sort(idx)].copy()
    ctrl["geocode_match"] = "rooftop"      # address points ARE rooftop locations
    log.info("drew %d controls from %d residential addresses", n, len(addresses))
    return ctrl


def match_filters(cases: gpd.GeoDataFrame, controls: gpd.GeoDataFrame):
    """Apply the SAME quality filter to both series."""
    c = cases[cases["geocode_match"].isin(GOOD_MATCH)]
    k = controls[controls["geocode_match"].isin(GOOD_MATCH)]
    log.info("after identical filtering: cases %d of %d, controls %d of %d",
             len(c), len(cases), len(k), len(controls))
    if len(c) < MIN_CASES:
        raise ValueError(f"{len(c)} cases after filtering — too few for a surface")
    return c, k


def control_check(controls: gpd.GeoDataFrame, blocks: gpd.GeoDataFrame,
                  pop_col: str = "population") -> float:
    """Correlation between control counts and population across blocks."""
    j = gpd.sjoin(controls[["geometry"]], blocks[["geoid", "geometry"]], predicate="within")
    counts = j.groupby("geoid").size().reindex(blocks["geoid"]).fillna(0)
    r = float(np.corrcoef(counts.to_numpy(), blocks[pop_col].to_numpy())[0, 1])
    log.info("control count vs block population correlation r=%.3f", r)
    if r < 0.8:
        log.warning("controls track population poorly — the surface will map the control design")
    return r


def run(cases, controls, blocks, region, h, relative_risk_fn, grid_fn, step):
    c, k = match_filters(cases, controls)
    control_check(k, blocks)
    grid = grid_fn(region, step)
    xy = lambda g: np.c_[g.geometry.x.to_numpy(), g.geometry.y.to_numpy()]
    res = relative_risk_fn(xy(c), xy(k), grid, region, h)
    out = gpd.GeoDataFrame({"rho": res["rho"], "p_upper": res["p_upper"]},
                           geometry=gpd.points_from_xy(grid[:, 0], grid[:, 1]),
                           crs=region.crs)
    sig = out[out["p_upper"] < 0.05]
    log.info("surface: %d cells, %d inside tolerance contour, max rho=%.2f",
             len(out), len(sig), float(out["rho"].max()))
    return out

The control-population check is the step that protects the whole analysis. If controls drawn from the address file correlate weakly with block population — because the address file under-covers rural routes or new subdivisions — the surface will show elevated risk wherever controls are sparse, which is a map of the address file’s gaps. A correlation above about 0.8 at block level is a reasonable minimum; below it, investigate the control source before estimating anything.

Filtering both series with the same match-type rule is equally important and easy to forget, because cases usually go through a quality filter in the surveillance pipeline while controls, drawn from an address file, are assumed to be perfect. If cases are filtered to rooftop matches and controls include everything, rural cases are under-represented relative to rural controls and the surface shows spurious low risk in the countryside.

Filter Both Series the Same Way Two runs on the same data. In the first, cases are filtered to rooftop geocodes, which removes a larger share of rural cases because rural addresses geocode worse, while controls are kept unfiltered. The resulting surface shows reduced risk across rural areas, entirely because rural cases were dropped. In the second, the same filter is applied to both series, rural cases and controls are reduced in proportion, and the rural surface returns to no difference. An asymmetric filter paints the countryside as low risk Cases filtered, controls not rural cases lost: 22% rural controls lost: 0% spurious low rural risk Same filter on both rural cases lost: 22% rural controls lost: 22% no rural difference controls from an address file must be put through the case pipeline's geocoder and filter

Validation & Edge Cases

Controls sampled once and reused for years. Population moves; a control set from five years ago misrepresents today’s at-risk population in growing areas. Redraw controls when the analysis period changes materially.

Institutional addresses. Nursing homes and prisons appear once in an address file and may hold hundreds of residents. Weight controls by household size, or handle group quarters separately, so the control density reflects people, not addresses.

Age-specific outcomes. For a paediatric outcome, controls should be drawn in proportion to the child population, not all residents. An all-ages control series maps where children live relative to adults.

Temporal mismatch. Cases accumulated over three years and controls drawn from this year’s address file differ by any construction in between. Align the address vintage with the middle of the case period.

Control-to-case ratio drifting across runs. A weekly surface re-estimated as cases accrue while the control set stays fixed changes its ratio of controls to cases every week. The log relative risk is invariant to that ratio in expectation, but its variance is not, so contours tighten or loosen from week to week for a reason unrelated to risk. Hold the ratio roughly constant by redrawing controls in proportion to the case count, or report the ratio alongside every surface.

Cases missing from the address file. Some case addresses — new builds, informal housing — are absent from the address-point file used for controls. Their areas have cases but systematically too few controls, which reads as elevated risk. Compare the share of cases whose address appears in the file across areas; a pronounced deficit in one area is a warning sign about that area’s surface.

Too-narrow tolerance contours. Small contours around one or two cases are usually noise. Require a contour to contain a minimum number of cases before reporting it.

INFO rr.run after identical filtering: cases 412 of 468, controls 1528 of 1640
INFO rr.run control count vs block population correlation r=0.884
INFO rr.run surface: 21406 cells, 318 inside tolerance contour, max rho=0.94

Before the map leaves the analysis environment, convert it to the form a program uses: contours, not a raster of raw values. A contour of the significant area with the maximum relative risk inside it is actionable and much less disclosive than a fine raster whose peaks mark individual addresses.

From Analytic Raster to Actionable Contour The analytic output is a fine raster of log relative risk and cell-wise p-values, which stays inside the protected environment. The exported product is a set of tolerance-contour polygons, each carrying the number of cases inside it rounded to a band, the maximum relative risk, and the bandwidth used. Contours containing fewer than a minimum number of cases are dropped. The export is what a program acts on and what disclosure review evaluates. Export contours, keep the raster inside analytic raster rho and p per cell protected filter contours drop those with too few cases contour product banded case count max relative risk bandwidth the contour product is what disclosure review evaluates and what a program acts on

Compliance Notes

Record the control strategy and source, the sampling seed, the geocoder and filter applied to both series, the control-population correlation, the bandwidth and its selection method, the edge correction, the simulation count and seed, and the minimum-case rule for reporting a contour.

Controls drawn from a residential address file are themselves locations of real households, even though those households are not cases. The control sample should be held inside the protected environment, never published, and discarded or retained under the same rules as the case file, since a published control set alongside a published case surface can narrow individual locations.