Getis-Ord Gi* Hotspot Detection

This guide is part of Disease Clustering & Spatial Statistical Modeling, and covers the production deployment of the Getis-Ord Gi* local statistic for feature-level hotspot and coldspot detection in public health surveillance. Gi* computes a localized z-score and p-value for every spatial unit, letting an agency rank census tracts or block groups by the statistical intensity of high-incidence clustering rather than collapsing a region into a single global coefficient. Use it during active surveillance of vector-borne outbreaks, opioid overdose surges, or environmental exposure events, where intervention budgets must be steered to specific, defensible geographies.

Concept & Epidemiological Alignment

Gi* measures the degree to which a feature and its defined neighborhood deviate from a spatially random distribution of the analysis variable. The statistic sums the analysis values inside the neighborhood of feature i (including i itself when computed in star form) and compares that local sum to what would be expected if all values were randomly redistributed across the study area. A large positive z-score marks a feature embedded in a high-value neighborhood — a hotspot; a large negative z-score marks a low-value neighborhood — a coldspot.

Three assumptions must hold before Gi* is epidemiologically valid:

• The analysis variable is a rate or count comparable across units. Raw case counts conflate incidence with population. Aggregate to an age-standardized or population-normalized rate, or model counts against an expected baseline, so that a hotspot reflects elevated risk rather than elevated denominator.
• Spatial support is consistent. Mixing tracts, ZCTAs, and point geocodes in one run produces a neighborhood structure that is not interpretable. Resolve to a single areal support before building weights.
• Underlying population is roughly homogeneous within neighborhoods, or its heterogeneity is captured by the weights. Where population density varies sharply, fixed-distance neighborhoods over- or under-bound the at-risk population; an adaptive neighborhood is required (covered under Optimizing Bandwidth for Getis-Ord Gi* Heatmaps).

Gi* answers where clustering concentrates. It does not, on its own, distinguish a focal point source from a diffuse gradient, and it does not test for clusters of arbitrary geometry across a time window. Choose the method that matches the surveillance question.

Method Selection

A common production pattern runs a global Moran’s I test first to confirm that any clustering exists, then drills into Gi* only when global autocorrelation is significant — this avoids interpreting local noise as signal.

The Gi* z-score for each feature places it on a signed significance scale: high values surrounded by high neighbors push the statistic above the upper critical threshold (a hotspot), low values surrounded by low neighbors push it below the lower threshold (a coldspot), and features near zero remain statistically indistinguishable from spatial randomness.

Spatial Data Prerequisites

Gi* operates on a single layer of areal polygons (or polygon centroids) carrying one numeric analysis column. Before computation, satisfy every prerequisite below — each maps to a validation gate enforced in the pipeline that follows.

• Geometry type: clean, valid polygons (or their centroids). Run gdf.is_valid.all() and repair with make_valid(); drop or flag null geometries explicitly rather than silently.
• Projection: a distance-preserving projected CRS, never a geographic one. Project to UTM (e.g., EPSG:32618 for Zone 18N) or the appropriate State Plane zone so that contiguity and distance neighborhoods are metrically accurate. Enforce canonical projections per Coordinate Reference Systems for Public Health before any neighborhood is built.
• Privacy de-identification: raw patient coordinates must never enter the analytical environment. Aggregate to stable areal units — census tracts, block groups, or ZCTAs — to satisfy HIPAA Safe Harbor and the disclosure thresholds documented in Compliance Mapping Frameworks.
• Analysis variable: a population-normalized rate or a count paired with an expected baseline, with no nulls in the analysis column. Impute or exclude missing values with a logged, domain-justified rule.
• Minimum sample size: permutation inference is unstable below roughly 30 features; below that, neighborhoods overlap so heavily that z-scores lose discriminating power. Report the feature count in execution metadata.
• Stable identifier: a unique, sortable feature ID (FIPS, GEOID) used to order the frame deterministically before computation, so repeated runs produce byte-identical output.

The implementation that follows realizes this five-stage pipeline:

The Gi* Statistic

The Gi* formula evaluates the ratio of the weighted sum of values in a local neighborhood to the global sum of values, standardized by the expected mean and variance under spatial randomness:

$$G_i^* = \frac{\sum_{j=1}^{n} w_{ij} x_j - \bar{X} \sum_{j=1}^{n} w_{ij}}{S \sqrt{\frac{n \sum_{j=1}^{n} w_{ij}^2 - (\sum_{j=1}^{n} w_{ij})^2}{n-1}}}$$

Where $w_{ij}$ is the spatial weight between features i and j, $x_j$ is the attribute value at j, $\bar{X}$ is the sample mean, and $S$ is the sample standard deviation. The resulting z-score indicates cluster directionality (positive = hot, negative = cold). Because local statistics evaluate hundreds or thousands of hypotheses simultaneously, uncorrected p-values yield unacceptable false-positive rates. The Benjamini-Hochberg False Discovery Rate (FDR) procedure is the production standard for spatial epidemiology, balancing statistical power against type I error control.

Production Implementation

The following implementation enforces row-standardization, handles isolated geometries, applies FDR correction, and outputs a validated GeoDataFrame ready for GIS rendering or API consumption. Pinned versions: geopandas==0.14.*, libpysal==4.11.*, esda==2.5.*, statsmodels==0.14.*, numpy==1.26.*.

import geopandas as gpd
import numpy as np
import libpysal
from libpysal.weights import KNN, Queen
from esda.getisord import G_Local
from statsmodels.stats.multitest import fdrcorrection

def compute_getis_ord_gi_star(
    gdf: gpd.GeoDataFrame,
    incidence_col: str,
    id_col: str,
    k_neighbors: int = 4,
    fdr_alpha: float = 0.05,
    permutations: int = 999,
    weight_type: str = "knn",
    seed: int = 42
) -> gpd.GeoDataFrame:
    """
    Production-ready Getis-Ord Gi* hotspot detection pipeline.
    Assumes gdf is projected to a metric CRS and incidence_col contains
    aggregated case counts or standardized rates.
    """
    # 0. Deterministic ordering: sort by a stable ID so reruns are byte-identical
    gdf = gdf.sort_values(id_col).reset_index(drop=True)
    np.random.seed(seed)

    # 1. Spatial validation: drop null geometries and verify CRS
    gdf = gdf.dropna(subset=["geometry"]).copy()
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("Geographic/undefined CRS. Project to a metric CRS (UTM/State Plane) before computation.")
    if gdf[incidence_col].isna().any():
        raise ValueError(f"Null values in '{incidence_col}'. Impute or exclude with a logged rule first.")

    # 2. Construct spatial weights matrix
    if weight_type.lower() == "knn":
        w = KNN.from_dataframe(gdf, k=k_neighbors)
    elif weight_type.lower() == "queen":
        w = Queen.from_dataframe(gdf)
        if w.islands:
            # Fall back to KNN when Queen weights produce isolated polygons
            w = KNN.from_dataframe(gdf, k=k_neighbors)
    else:
        raise ValueError("Unsupported weight_type. Use 'knn' or 'queen'.")

    w.transform = "r"  # Row-standardize

    # 3. Compute Gi* statistic (star=True includes the focal feature in its own neighborhood)
    y = gdf[incidence_col].values.astype(float)
    gi_star = G_Local(y, w, star=True, permutations=permutations, seed=seed)

    # 4. Extract results and apply FDR correction
    z_scores = gi_star.Zs
    raw_p_values = gi_star.p_sim

    # Handle potential NaNs from permutations or isolated features
    valid_mask = ~np.isnan(z_scores)
    corrected_p = np.full_like(raw_p_values, np.nan)
    if valid_mask.any():
        _, corrected_p[valid_mask] = fdrcorrection(raw_p_values[valid_mask], alpha=fdr_alpha)

    # 5. Assemble output GeoDataFrame
    out_gdf = gdf.copy()
    out_gdf["gi_z_score"] = z_scores
    out_gdf["gi_p_raw"] = raw_p_values
    out_gdf["gi_p_fdr"] = corrected_p
    out_gdf["gi_significant"] = (corrected_p < fdr_alpha) & valid_mask
    out_gdf["gi_cluster_type"] = np.select(
        [
            (out_gdf["gi_significant"]) & (out_gdf["gi_z_score"] > 0),
            (out_gdf["gi_significant"]) & (out_gdf["gi_z_score"] < 0)
        ],
        ["hotspot", "coldspot"],
        default="not_significant"
    )

    return out_gdf

Parameter Selection & Tuning

Four parameters drive Gi* output; each must be chosen against the surveillance context and recorded for reproducibility.

• Weights type. Contiguity-based weights (Queen or Rook) are standard for complete areal coverage, where adjacency encodes a meaningful “neighbor” relationship. K-nearest-neighbor or distance-band weights are preferred when administrative boundaries are fragmented, when polygons vary wildly in size, or when modeling mobile populations whose exposure is not bounded by tract lines.
• Row-standardization. Setting w.transform = "r" so each feature’s weights sum to 1 is mandatory; it stabilizes variance across irregularly shaped polygons and prevents large-perimeter units from dominating the statistic.
• Neighborhood size (k or bandwidth). Overly dense neighborhoods dilute localized intensity into the regional mean; sparse neighborhoods amplify noise and edge effects. For continuous surveillance, use adaptive bandwidths that scale with population density rather than a fixed Euclidean radius — radius calibration and kernel decay functions are detailed in Optimizing Bandwidth for Getis-Ord Gi* Heatmaps.
• Significance threshold and FDR strategy. Apply Benjamini-Hochberg FDR at a documented alpha (commonly 0.05). Run a sensitivity analysis across several k values or contiguity thresholds and confirm that the set of significant hotspots is stable; clusters that appear only at a single parameter setting are candidates for review, not action.

Gi* and Local Moran’s I Are Not Interchangeable

The two statistics are computed from the same weights matrix on the same variable, they are reported on the same map, and they are routinely treated as alternative renderings of one result. They are not. They answer different questions, and the difference shows up on exactly the units a surveillance programme cares most about.

Gi* asks whether the sum of values in a neighbourhood, including the unit itself, is higher than you would expect if the values were shuffled at random. It is a statement about a neighbourhood’s total. Local Moran’s I asks whether a unit’s value covaries with its neighbours’ values — whether it is like them, or strikingly unlike them. It is a statement about similarity. A high value surrounded by high values scores strongly on both. Everything else diverges.

Three consequences follow, and each has bitten a real surveillance programme.

Gi* cannot report a spatial outlier. A tract with a sharp excess whose neighbours are ordinary produces a neighbourhood sum barely above expectation, so its z-score is unremarkable. Local Moran’s I classifies the same tract as high-low — a genuine anomaly worth a phone call. If the programme’s brief includes “find single-facility or single-institution outbreaks”, Gi* alone is the wrong instrument.

Gi* does not distinguish a cold cluster from an unremarkable one in the same way. It reports low-sum neighbourhoods as significantly cold, which is useful, but it has no notion of a low value sitting among high ones. In screening-coverage work that low-high case is often the operational finding: the one tract in a well-covered region that is being missed.

Their significance masks differ even where they agree. Because Gi* includes the focal unit in its own sum while local Moran’s I does not, the two statistics have different variances and therefore different tails on identical data. Reporting “significant at 0.05” from each without saying which is which invites a reviewer to compare counts that are not comparable.

The practical recommendation is to compute both, always, and to publish them as two layers rather than one. They cost the same weights matrix and one extra pass over the data. When they agree, the finding is robust to the choice of statistic — which is itself worth stating. When they disagree, the disagreement localises the type of anomaly, and that is more information than either layer carries alone.

Edge Cases & Failure Modes

• Island polygons. Contiguity weights produce features with zero neighbors (w.islands), which yield NaN z-scores that propagate through downstream joins. Detect islands via w.islands and fall back to KNN, or attach a domain-justified nearest-neighbor injection — never silently drop them.
• Zero-inflation. Sparse case data with many zero-count units flattens the variance term and can manufacture spurious coldspots. Use a population-offset rate or smooth small-area estimates before running Gi*, and flag units whose denominator falls below a reliability threshold.
• Transboundary CRS drift. Multi-jurisdiction studies that span UTM zones or mix State Plane zones introduce distance error at the seams. Reproject every input to one documented CRS and assert a single gdf.crs before building weights; see Coordinate Reference Systems for Public Health.
• Memory constraints for N > 50k. Dense permutation inference scales poorly past tens of thousands of features. Use sparse weight representations, reduce permutations only after confirming p-value stability, and consider tiling the study area into overlapping panes with a buffer to preserve edge neighborhoods.
• Edge effects. Features on the study boundary have truncated neighborhoods and biased statistics. Buffer the analysis extent or flag boundary units in execution metadata rather than reporting them as confirmed hotspots.

Compliance & Audit Controls

Federal audit and inter-agency data-sharing protocols require that any published hotspot be reproducible from logged inputs.

• Deterministic execution. Sort by a stable feature ID before computation and pin the permutation seed, so every rerun reproduces the same z-scores, p-values, and classifications.
• Configuration logging. Record weights type, k/bandwidth, fdr_alpha, permutation count, seed, CRS EPSG code, temporal window, and aggregation method in a metadata registry. Hash the input geometry file and the configuration block with SHA-256 and attach both hashes to every output.
• Output schema. Emit explicit, documented columns — gi_z_score, gi_p_raw, gi_p_fdr, gi_significant, gi_cluster_type — alongside the stable ID, and write ISO 19115 lineage metadata (projection, source, processing steps) with the file. Prefer serialization formats covered in Spatial Data Types & Formats so provenance survives interagency handoff.
• Validation cross-check. When the surveillance program also holds point-level data, cross-validate Gi* areal hotspots against K-Function & Point Pattern Analysis to confirm that detected clustering is not an artifact of population heterogeneity.

Reporting a hotspot so it can be acted on

A flagged tract is the beginning of an investigation, not its conclusion, and the fields that make it actionable are cheap to emit at the moment the statistic is computed. Alongside the z-score and the corrected p-value, record the neighbourhood that produced the score — the list of unit identifiers inside the band — and the count and expected count that went into it. An investigator who receives only a z-score has to reconstruct the neighbourhood by hand before they can decide whether the cluster spans a jurisdictional boundary or sits inside a single service area.

It is equally worth recording what the tract’s own value contributed relative to its neighbours. Because Gi* includes the focal unit in its sum, a high score can come from a genuinely elevated neighbourhood or from one very high unit dragging an otherwise ordinary neighbourhood upward. Those two situations call for different responses, and the ratio of the focal value to the neighbourhood mean distinguishes them in a single number that costs nothing to compute and is impossible to recover later.

Production Readiness Checklist

• All inputs reprojected to one documented metric CRS; crs.is_geographic asserted false
• Patient-level coordinates aggregated to areal units that clear disclosure thresholds
• Analysis column is a population-normalized rate or count with a baseline; no nulls
• Frame sorted by a stable ID and permutation seed pinned for byte-identical reruns
• Weights row-standardized (w.transform = "r"); islands detected and handled, not dropped
• FDR (Benjamini-Hochberg) applied at a documented alpha; sensitivity sweep over k logged
• Weights type, k/bandwidth, alpha, seed, CRS, and input hash recorded in the metadata registry
• Output written with explicit z/p/classification columns and ISO 19115 lineage metadata

Related Topics

• Disease Clustering & Spatial Statistical Modeling — parent overview of the end-to-end surveillance pipeline.
• Global & Local Moran’s I Implementation — autocorrelation testing and spatial-outlier detection that complements Gi*.
• K-Function & Point Pattern Analysis — second-order clustering for raw point geocodes.
• Spatial Scan Statistics Configuration — variable-radius space-time cluster scanning.
• Optimizing Bandwidth for Getis-Ord Gi* Heatmaps — adaptive neighborhood calibration for this method.