Spatial Epi Public Health GIS

Testing Significance of a Relative Risk Surface

This guide, part of Kernel Density & Spatial Relative Risk Surfaces, deals with the inferential step that turns a smooth map into a claim: deciding which areas of elevated relative risk are real enough to report.

Problem Context & Constraints

The standard output of a relative risk analysis is a surface with tolerance contours: regions where the observed log relative risk exceeds, say, 95% of its values under random relabelling of cases and controls. Each grid cell gets its own Monte Carlo p-value, and the contour outlines cells below 0.05.

The difficulty is that there are thousands of cells. Even when cases and controls come from exactly the same spatial distribution, around 5% of cells will fall below 0.05 by chance, and because neighbouring cells share most of their kernel mass, those chance exceedances cluster into plausible-looking contours. A surface built from pure noise will, most of the time, show one or two outlined regions.

Three constraints shape a defensible procedure. Test globally first: a single statistic asking whether there is any spatial variation in risk at all, with a proper null distribution. If it is not significant, the local contours are exploratory at best. Control multiplicity locally, either by a false discovery rate adapted to spatial dependence or by requiring contours to exceed a threshold derived from the maximum statistic under the null. And require substance: a contour enclosing two cases is not a finding whatever its p-value.

Contours From Pure Noise A relative risk surface estimated from simulated data in which cases and controls were drawn from the same distribution, so there is no true variation in risk. Cell-wise tolerance contours at the 0.05 level nevertheless outline two small regions, because roughly five per cent of cells cross the threshold by chance and neighbouring cells are correlated. The same surface tested globally returns a non-significant result, and a maximum-statistic threshold outlines nothing, which is the correct answer. No true variation in risk — and still two contours cell-wise p < 0.05 two regions outlined by chance global test, then max-statistic global p = 0.61 nothing outlined the correct answer for noise simulated data with identical case and control distributions

Prerequisites

  • python 3.11
  • numpy 1.26.4, scipy 1.13.1
  • The relative risk estimation routine and a grid, from the parent guide
  • Case and control coordinates, with the chosen bandwidth
  • A simulation budget: at least 999 relabellings for a published map

Step-by-Step Solution

Compute a global statistic — the integrated squared log relative risk — and its null distribution from relabellings; then use the same relabellings to derive a maximum-statistic threshold for local contours; finally filter contours by case count.

# Global and multiplicity-aware local tests for a log relative risk surface.
# Pinned: numpy==1.26.4, scipy==1.13.1
import logging

import numpy as np
from scipy import ndimage

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

SEED = 20260918


def global_statistic(rho: np.ndarray, cell_area: float) -> float:
    """Integrated squared log relative risk: zero under a constant ratio."""
    return float(np.sum(rho ** 2) * cell_area)


def relabel_surfaces(cases, controls, grid, estimate_rho, n_sim: int, seed: int = SEED):
    """Yield null surfaces from random relabelling of the pooled points."""
    rng = np.random.default_rng(seed)
    pooled = np.vstack([cases, controls])
    n1 = len(cases)
    for _ in range(n_sim):
        perm = rng.permutation(len(pooled))
        yield estimate_rho(pooled[perm[:n1]], pooled[perm[n1:]], grid)


def test_surface(cases, controls, grid, cell_area, estimate_rho,
                 n_sim: int = 999, alpha: float = 0.05) -> dict:
    rho = estimate_rho(cases, controls, grid)
    t_obs = global_statistic(rho, cell_area)
    t_null, max_null = [], []
    for rs in relabel_surfaces(cases, controls, grid, estimate_rho, n_sim):
        t_null.append(global_statistic(rs, cell_area))
        max_null.append(rs.max())
    t_null = np.asarray(t_null)
    p_global = (1 + np.sum(t_null >= t_obs)) / (n_sim + 1)
    # Max-statistic threshold: controls the chance of ANY false cell under the null.
    threshold = float(np.quantile(max_null, 1 - alpha))
    local = rho >= threshold
    log.info("global test: T=%.3f p=%.4f; max-stat threshold rho>=%.3f outlines %d cells",
             t_obs, p_global, threshold, int(local.sum()))
    return {"rho": rho, "p_global": float(p_global), "threshold": threshold,
            "significant_cells": local}


def filter_regions(grid_shape, significant_cells, case_counts, min_cases: int = 5):
    """Keep only contiguous significant regions containing enough cases."""
    mask = significant_cells.reshape(grid_shape)
    lab, n = ndimage.label(mask)
    keep = np.zeros_like(mask)
    kept = 0
    for r in range(1, n + 1):
        cells = lab == r
        if case_counts.reshape(grid_shape)[cells].sum() >= min_cases:
            keep |= cells
            kept += 1
    log.info("regions: %d outlined, %d kept with >= %d cases", n, kept, min_cases)
    return keep.ravel()

Reusing the same relabellings for the global test and the local threshold keeps the computation affordable: the expensive part is re-estimating the surface, and each relabelled surface feeds both statistics. The maximum-statistic threshold is conservative — it controls the probability of outlining any cell at all under the null — which is appropriate for a map whose contours will be read as findings. A false discovery rate approach is less conservative and suits exploratory screening where some false contours are acceptable.

The region filter applies the substance requirement. A contour that encloses two cases may be statistically extreme because the local control density happens to be low, but it is not something a program can investigate, and publishing it carries disclosure risk. A minimum case count per region removes those.

Where the Two Thresholds Sit Two null distributions from the same relabellings. The cell-wise distribution of log relative risk at a single grid cell is narrow, and its ninety-fifth percentile is low. The distribution of the maximum log relative risk anywhere on the grid is shifted well to the right, because the maximum over thousands of correlated cells is almost always larger than any single cell's value. The ninety-fifth percentile of the maximum is the threshold that controls false contours across the whole map. The map-wide threshold is set by the maximum, not by one cell one cell under the null cell-wise 95% maximum over the grid map-wide 95% log relative risk under random relabelling cells between the two thresholds are the ones cell-wise contours wrongly outline

Validation & Edge Cases

Global test not significant but local contours present. Report the global result as the headline. Local contours may be shown as exploratory, clearly labelled, but should not be described as clusters.

Too few simulations. The maximum-statistic threshold is an upper quantile of a maximum, and with 99 relabellings it is estimated from five or so values. Use at least 999 for a published threshold.

Adaptive bandwidths in relabelling. When the surface uses adaptive bandwidths, relabelling must keep each point’s bandwidth attached to its location, not its label, so that the null reflects only case status. Re-deriving bandwidths per relabelling is also valid but costlier; mixing the two approaches is not.

Edge regions. Near the boundary the edge-corrected density has higher variance, so the maximum often occurs at the edge under the null. That inflates the map-wide threshold for the interior. Where boundary cells dominate the null maximum, compute the threshold over interior cells only and report boundary contours separately.

Repeated testing in surveillance. A surface re-estimated and re-tested every week accumulates chances for a false contour just as a daily scan statistic does. The map-wide threshold controls false contours within one map, not across a year of maps. For prospective use, calibrate the threshold against the history of weekly surfaces, following the approach in Calibrating Alert Thresholds Against Historical Signals.

Low-risk regions. The procedure above tests for elevated risk. Regions of unusually low risk are tested with the lower tail of the null and the minimum rather than the maximum statistic. Report them separately, because they answer a different question and are more often driven by under-ascertainment than by genuinely lower risk.

Different bandwidths across relabellings. Re-selecting the bandwidth inside each relabelling adds variability that the observed surface did not have. Fix the bandwidth from the observed data and hold it constant across simulations.

INFO rr.test global test: T=18.412 p=0.0030; max-stat threshold rho>=0.871 outlines 142 cells
INFO rr.test regions: 3 outlined, 2 kept with >= 5 cases

Reporting what the contours mean. A significant contour says that, within it, cases are more concentrated relative to controls than random relabelling produces. It does not locate a source, estimate the magnitude of risk precisely, or establish causation; the log relative risk inside the contour is itself a smoothed estimate with its own uncertainty. The report should state the global p-value, the threshold method, the minimum-case rule and the bandwidth alongside the map.

From Surface to Reportable Regions A three-stage sequence. The global test asks whether any spatial variation in risk exists; if not, the map is labelled exploratory and no clusters are claimed. If so, the map-wide maximum-statistic threshold outlines cells that exceed chance across the whole grid. Contiguous outlined regions are then kept only if they contain a minimum number of cases. The regions that survive all three stages are the ones reported, with the global p-value, threshold method and case-count rule stated. Three gates between a surface and a reported cluster 1 · global test any variation at all? 2 · map-wide threshold maximum-statistic 3 · case-count filter enough cases to act on not significant: label map exploratory report surviving regions with the global p-value, threshold method and case rule

Compliance Notes

Record the global statistic and its p-value, the number of relabellings and the seed, the threshold method and value, whether edge cells were treated separately, the minimum-case rule and the number of regions outlined and kept.

A significant contour is often the trigger for public health action and sometimes for public communication, so the inferential record behind it should be as complete as for any published statistic. The regions reported publicly should also pass disclosure review: a small contour with a minimum of five cases is a statement that at least five people with the condition live in a defined small area, and whether that can be published is a separate decision from whether it is significant.