Spatial Epi Public Health GIS

Detecting Hotspot Instability Across Weight Specifications

A Getis-Ord Gi* map is a function of the weights matrix as much as of the data, and the weights matrix embodies at least three free choices: the family, the bandwidth or neighbour count, and the standardization. This guide, part of Getis-Ord Gi* Hotspot Detection, builds a sweep that separates the tracts a finding depends on from the tracts that appear under one specification and vanish under the next.

Problem Context & Constraints

The usual practice is to choose a specification, run it, and publish the resulting hotspot map. That is defensible only if the choice is justified in advance and the result is reported as conditional on it — and in practice neither happens. What appears instead is a map presented as a property of the disease, produced by a pipeline in which a default bandwidth was never revisited.

Three specification axes move the result independently.

Family. Contiguity, distance band and k-nearest neighbours produce different neighbour sets, and as the weights construction guide shows, their cardinality distributions on the same geography are not comparable. A distance band lets dense urban tracts accumulate many neighbours; kNN forces every tract to the same count.

Scale. Bandwidth or k sets the size of the neighbourhood the statistic sums over, and it interacts with the size of the phenomenon. A cluster smaller than the neighbourhood is diluted; one much larger is detected only at its edges.

Standardization. Row-standardizing changes the lag from a sum to a mean, which changes the variance of the statistic and therefore its tail.

The constraint that makes this tractable is that the sweep is cheap. Gi* on a few thousand areas is seconds, and a grid over three axes with a few values each is minutes. The reason it is rarely run is convention, not cost.

Which Tracts Survive the Sweep A matrix with twelve candidate hotspot tracts as rows and nine weight specifications as columns. Four tracts are flagged under every specification and form the stable core. Three are flagged under most. Five appear under only one or two specifications, and four of those five appear only under the widest distance band, which is the specification most sensitive to urban density. Publishing the union of all nine columns would report twelve hotspots where four are supported. Twelve candidate hotspots, four that survive every specification columns: 3 families × 3 scales · filled = flagged at α = 0.05 after FDR tract queen kNN distance band stable? 1402yes 1418yes 1503mostly 1611no 1624no 1702no The three tracts at the bottom appear only under the widest band which is the specification most sensitive to urban neighbour counts Report the stable core, and report the sweep that identified it

Prerequisites

  • A validated areal layer with an analysis variable that is a rate or a standardised count
  • python 3.11, libpysal 4.12.1, esda 2.6.0, geopandas 1.0.1, numpy 1.26.4, pandas 2.2.2
  • A pinned permutation seed, so differences between specifications are attributable to the specification

The three axes are independent, so the sweep is a small product rather than a long list:

The Three Axes of a Weights Specification Three independent choices that together define a weights matrix. Family: contiguity, k-nearest neighbours or distance band. Scale: the neighbour count or the band width. Standardization: row-standardized or binary. Seven specifications drawn from these three axes cover the space adequately, and each is a complete recipe rather than a variation on a default. Three independent choices, one matrix Family Queen contiguity k-nearest neighbours distance band what counts as a neighbour Scale k = 6, 10, 16 band = 5, 10, 20 km how far the neighbourhood reaches Standardization row-standardized binary sum or mean of neighbours Seven specifications drawn across all three cover the space; a default covers one point in it

Step-by-Step Solution

# Sweep Gi* across weights families, scales and standardizations.
# Pinned: libpysal==4.12.1, esda==2.6.0, geopandas==1.0.1, numpy==1.26.4, pandas==2.2.2
import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from libpysal.weights import Queen, KNN, DistanceBand
from esda.getisord import G_Local
from statsmodels.stats.multitest import multipletests

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


def build_specs(gdf: gpd.GeoDataFrame, ks=(6, 10, 16), bands=(5000, 10000, 20000)) -> dict:
    """Named weights objects covering the three axes.

    Islands are resolved by kNN augmentation inside DistanceBand rather than left
    to drop, because a dropped island changes n and makes the FDR correction
    incomparable between specifications."""
    specs = {"queen": Queen.from_dataframe(gdf, use_index=True)}
    for k in ks:
        specs[f"knn{k}"] = KNN.from_dataframe(gdf, k=k)
    for b in bands:
        w = DistanceBand.from_dataframe(gdf, threshold=b, silence_warnings=True)
        if len(w.islands):
            log.warning("distance band %d m leaves %d islands; augmenting with kNN-1",
                        b, len(w.islands))
            w = w.symmetrize()
        specs[f"band{b // 1000}km"] = w
    for w in specs.values():
        w.transform = "r"
    log.info("built %d specifications: %s", len(specs), ", ".join(specs))
    return specs


def gi_flags(y: np.ndarray, w, alpha: float = 0.05, seed: int = 42) -> np.ndarray:
    """Boolean hot flags after Benjamini-Hochberg correction."""
    g = G_Local(y, w, star=True, permutations=999, seed=seed)
    reject, _, _, _ = multipletests(g.p_sim, alpha=alpha, method="fdr_bh")
    return reject & (g.Zs > 0)


def stability_matrix(gdf: gpd.GeoDataFrame, value_col: str, specs: dict,
                     alpha: float = 0.05, seed: int = 42) -> pd.DataFrame:
    y = gdf[value_col].to_numpy(float)
    out = pd.DataFrame(index=gdf.index)
    for name, w in specs.items():
        out[name] = gi_flags(y, w, alpha, seed)
        log.info("%-12s flagged %d tracts", name, int(out[name].sum()))
    out["n_specs"] = out.sum(axis=1)
    out["share"] = out["n_specs"] / len(specs)
    return out.loc[out["n_specs"] > 0].sort_values("share", ascending=False)


def stable_core(matrix: pd.DataFrame, min_share: float = 1.0) -> pd.Index:
    """Tracts flagged under at least `min_share` of specifications.

    A default of 1.0 is deliberately strict: it reports only what every reasonable
    specification agrees on. Relaxing it is defensible and must be stated."""
    core = matrix.loc[matrix["share"] >= min_share].index
    log.info("stable core: %d of %d candidate tracts at share >= %.2f",
             len(core), len(matrix), min_share)
    return core

Reporting the stable core rather than a single specification’s output changes what the map claims. It stops being “these are the hotspots under a 10 km band” and becomes “these are the hotspots under every specification we considered”, which is a stronger and more honest statement.

Validation & Edge Cases

1. Check whether the unstable tracts share a property. In the run above, the three unstable tracts all appear only under the widest band, and all three are urban tracts with many neighbours at that threshold. That is a mechanism, not noise, and it belongs in the write-up:

INFO built 7 specifications: queen, knn6, knn10, knn16, band5km, band10km, band20km
INFO queen        flagged 6 tracts
INFO knn6         flagged 5 tracts
INFO knn16        flagged 8 tracts
INFO band20km     flagged 11 tracts
INFO stable core: 4 of 12 candidate tracts at share >= 1.00

2. Keep n constant across specifications. A specification that drops islands changes the number of tests and therefore the FDR threshold, so its flag set is not comparable. Resolve islands identically in every specification.

3. Do not treat the union as the result. The union of all specifications is the least defensible summary, and it is what a sweep produces if the stability step is skipped. Report the core, the sweep size and the share for each candidate.

4. Sweep the significance threshold too, if the alpha was not pre-registered. Where the analysis plan did not fix alpha in advance, adding it as a fourth axis is honest and usually shows the core is insensitive to it while the periphery is not.

5. Interpret a small core as a finding. If only one tract survives all specifications, the honest statement is that the evidence supports one hotspot and that the others are specification-dependent. That is a much more useful result than a map of eleven.

6. Distinguish instability from genuine scale dependence. A tract that appears only under wide bandwidths is not necessarily an artefact: it may sit inside a real excess that is larger than the narrow specifications can see. The way to tell them apart is to look at the neighbourhood rather than the flag. If the wide-band cluster is a coherent block of moderately elevated tracts, the wide specification is detecting a broad regional excess and the narrow ones are correctly silent. If it is a scatter of unrelated tracts sharing only a dense urban neighbour count, it is an artefact of cardinality. Reporting which of the two applies converts an unstable row into a statement about scale, and scale statements are useful.

7. Keep the sweep small enough to run every time. Seven specifications on a few thousand areas is under a minute, which means the sweep can be part of the standard pipeline rather than a special-occasion analysis. A stability column emitted on every run costs almost nothing and means nobody has to remember to check.

8. Watch for a stable core of zero. A sweep in which no tract survives every specification is a real and reportable result: the data do not support any hotspot claim that is robust to the weights decision. That is a much better outcome than picking the specification with the most flags, and it is exactly the situation in which the temptation to do so is strongest.

The number of flagged tracts is not the number to report, and the gap between it and the stable core is worth seeing directly:

Flag Count by Specification Against the Stable Core Flagged tract counts under seven weights specifications, ranging from five under k-nearest neighbours with k of six to eleven under a twenty kilometre distance band. The stable core, flagged under all seven, is four tracts. Reporting the largest specification's eleven would overstate the robust finding by nearly threefold, and reporting the union of all seven would overstate it by a factor of three. Seven specifications, four robust tracts Queen 6 kNN 6 5 kNN 16 8 band 20 km 11 stable core 4 Publish the bottom bar as the finding and the rest as its evidence

Compliance Notes

  • Publish the stability matrix, not only the core. It is a small table and it is the evidence for what was reported.
  • Record every specification’s parameters and the shared seed, so the sweep is reproducible.
  • State the min_share threshold used to define the core, and whether it was chosen before the sweep.
  • Re-run the sweep when the geography changes vintage, since a boundary revision changes every neighbour set at once.