Hot Spot Analysis with Getis-Ord Statistics

Hot spot analysis answers a question that global spatial statistics cannot: where on the map do high or low values concentrate strongly enough to be unlikely under spatial randomness? Epidemiologists use it to locate disease clusters, crime analysts to target patrol resources, and environmental scientists to flag pollution intensification zones. The Getis-Ord family of local statistics, $G_i$ and $G_i^*$, provides the rigorous answer by attaching a z-score and a p-value to every observation. This page develops the estimator, its dependence on the spatial weight matrix, permutation inference, and the multiple-comparison correction that separates real hot spots from statistical noise, all within the broader context of Core Concepts of Spatial Statistics & Geostatistics.

Prerequisites #

• Python 3.9+
• numpy>=1.22, geopandas>=0.13, libpysal>=4.8, esda>=2.5
• A GeoDataFrame in a projected (metric) CRS — distance-band neighbourhoods require metres, not degrees
• A numeric attribute column with no NaN values and non-zero variance
• A validated spatial weight matrix with a single connected component and no unintended islands

Mathematical Core #

The Getis-Ord statistic quantifies whether the sum of attribute values in a location’s neighbourhood departs from what the global distribution would predict. For the starred variant $G_i^*$, which includes the focal observation $i$ in its own neighbourhood, the standardised (z-score) form is:

$$G_i^* = \frac{\sum_{j=1}^{n} w_{ij}\, x_j - \bar{x} \sum_{j=1}^{n} w_{ij}}{S \sqrt{\dfrac{n \sum_{j=1}^{n} w_{ij}^2 - \left(\sum_{j=1}^{n} w_{ij}\right)^2}{n-1}}}$$

where $x_j$ is the attribute value at location $j$, $w_{ij}$ are elements of a binary weight matrix (with $w_{ii} = 1$ for $G_i^*$), $n$ is the number of observations, $\bar{x}$ is the global mean, and $S$ is the global standard deviation:

$$\bar{x} = \frac{1}{n}\sum_{j=1}^{n} x_j, \qquad S = \sqrt{\frac{1}{n}\sum_{j=1}^{n} x_j^2 - \bar{x}^2}$$

The numerator is the difference between the observed local sum and its expected value under a random spatial arrangement of the same data. The denominator is the standard error of that sum. Because the result is already standardised, $G_i^*$ is read directly as a z-score: large positive values mark hot spots (clusters of high values), large negative values mark cold spots (clusters of low values), and values near zero indicate no local concentration.

The non-starred $G_i$ excludes the focal unit ($w_{ii} = 0$ and the summations run over $j \neq i$), so it measures the neighbourhood context around $i$ without $i$’s own value. In practice $G_i^*$ is the default for hot spot mapping: a census tract that is itself a pollution maximum should count toward whether it lies inside a hot spot.

The Role of the Weight Matrix #

Getis-Ord statistics are defined on binary weights, unlike Moran’s I which is typically computed on row-standardised weights. The reason is structural: $G_i^*$ is a sum of neighbour values compared to a global expectation, not a weighted average. Row-standardising the matrix would divide every neighbourhood sum by its neighbour count, collapsing the very magnitude information the statistic is built to detect. Keep the transform binary and, for $G_i^*$, set the diagonal $w_{ii} = 1$ so each unit includes itself.

The neighbourhood definition remains a modelling decision. A fixed distance band expresses a hypothesis that influence operates within a known radius; k-nearest-neighbour weights guarantee a constant neighbourhood size and avoid isolated units in sparse point data. Whichever topology you choose, the same attribute data will produce a different hot spot map under a 2 km band than under an 8 km band. Selecting and validating that structure is the subject of the spatial weight matrices guide.

Permutation Inference #

The analytical z-score above assumes the local sum is approximately normal, which holds poorly for skewed data or small neighbourhoods. Conditional permutation inference sidesteps that assumption. For each location $i$, the values at all other locations are randomly reshuffled many times while $i$’s neighbourhood cardinality is held fixed, producing an empirical reference distribution of $G_i^*$ under the null of no local association. The pseudo p-value is:

$$p_{sim} = \frac{m + 1}{M + 1}$$

where $M$ is the number of permutations and $m$ is the count of simulated statistics at least as extreme as the observed one. With $M = 999$ the smallest resolvable p-value is $0.001$. Permutation inference is distribution-free and is the recommended default for applied hot spot mapping.

Multiple-Comparison Correction #

A hot spot map runs one hypothesis test per unit. Across $n = 3{,}000$ tracts, an uncorrected threshold of $\alpha = 0.05$ would flag roughly 150 units as significant purely by chance. Two corrections are common. The conservative approach divides $\alpha$ by the number of tests (a Bonferroni-style bound), controlling the family-wise error rate but discarding many true signals. The false discovery rate (FDR) approach of Benjamini and Hochberg is the practical middle ground: sort the $p$ p-values ascending, find the largest rank $k$ for which $p_{(k)} \le \frac{k}{n}\alpha$, and treat all units up to that rank as significant. FDR controls the expected proportion of false discoveries and typically retains far more real hot spots than a family-wise bound.

Contrast with Spatial Autocorrelation Metrics #

Getis-Ord and Local Moran’s I are complementary local statistics, and it is worth being precise about what each detects. Local Moran’s I, one of the spatial autocorrelation metrics, is a local correlation and returns four labels: high-high and low-low clusters, plus high-low and low-high spatial outliers. Getis-Ord $G_i^*$ is a local sum relative to the global mean and returns only two kinds of significant unit — hot spots and cold spots. It cannot flag a spatial outlier because a low value surrounded by highs still contributes a large neighbourhood sum and reads as part of a hot spot. Use $G_i^*$ when the question is “where are the intense high or low clusters”; use LISA when you also need to isolate observations that contradict their neighbourhood.

Annotated Implementation #

The esda library implements both variants through the G_Local class. The star=True argument switches from $G_i$ to $G_i^*$ by including the self-weight.

1. Data Preparation and Binary Weights #

import geopandas as gpd
import numpy as np
import libpysal
from esda.getisord import G_Local

# Load and reproject to a metric CRS (UTM Zone 18N)
gdf = gpd.read_file("incidents.gpkg").to_crs("EPSG:26918")
gdf = gdf.reset_index(drop=True)          # align rows with weight indices

# Build a fixed distance-band weight matrix with BINARY weights.
# binary=True is essential — Getis-Ord uses sums, not weighted averages.
w = libpysal.weights.DistanceBand.from_dataframe(
    gdf, threshold=3000.0, binary=True, silence_warnings=True
)

DistanceBand yields binary weights by default; contiguity weights (Queen/Rook) are also valid for areal data. Do not call w.transform = "r" — row-standardisation is wrong for this statistic.

2. Computing Gi* with Permutation Inference #

y = gdf["incident_rate"].values.astype(float)

# star=True includes the focal unit (Gi*); permutations drive inference.
gi_star = G_Local(y, w, star=True, permutations=999, seed=42)

# Key attributes on the fitted object:
#   gi_star.Zs      -> analytical z-scores per unit
#   gi_star.p_sim   -> permutation pseudo p-values
#   gi_star.Gs      -> raw Gi* values (before standardisation)
print("z-score range:", gi_star.Zs.min().round(2), "to", gi_star.Zs.max().round(2))

The G_Local signature is G_Local(y, w, transform='R', permutations=999, star=False, ...). Passing star=True sets the self-weight internally so you do not modify w yourself. The Zs array holds the standardised statistic; p_sim holds the permutation p-values used for classification.

3. False Discovery Rate Correction #

def fdr_threshold(p_values, alpha=0.05):
    """Benjamini-Hochberg critical p-value for a family of tests."""
    p_sorted = np.sort(p_values)
    n = len(p_values)
    ranks = np.arange(1, n + 1)
    below = p_sorted <= (ranks / n) * alpha
    if not below.any():
        return 0.0                       # nothing survives correction
    return p_sorted[below].max()

p_crit = fdr_threshold(gi_star.p_sim, alpha=0.05)
print(f"FDR-corrected p threshold: {p_crit:.4f}")

Every unit with p_sim <= p_crit is treated as significant. When no unit clears the bar, p_crit is 0.0 and the map correctly shows no hot spots.

4. Classification into Hot and Cold Spots #

def classify_getis_ord(z_scores, p_values, p_crit):
    """Return labels: 'Hot Spot', 'Cold Spot', or 'Not Significant'."""
    labels = np.full(len(z_scores), "Not Significant", dtype=object)
    significant = p_values <= p_crit
    labels[significant & (z_scores > 0)] = "Hot Spot"
    labels[significant & (z_scores < 0)] = "Cold Spot"
    return labels

gdf["gi_class"] = classify_getis_ord(gi_star.Zs, gi_star.p_sim, p_crit)
print(gdf["gi_class"].value_counts())

The sign of the z-score sets the direction (hot vs cold); the corrected p-value sets significance. This two-part rule is the whole of hot spot classification.

Diagnostic and Confidence Configuration #

Analysts commonly report hot spots at graded confidence levels rather than a single threshold, mirroring the three-band convention (90%, 95%, 99%) used in mapping software.

def confidence_band(z, p):
    """Signed confidence band from z-sign and permutation p-value."""
    az = abs(z)
    if p <= 0.01 and az >= 2.58:
        band = 99
    elif p <= 0.05 and az >= 1.96:
        band = 95
    elif p <= 0.10 and az >= 1.65:
        band = 90
    else:
        return 0
    return band if z > 0 else -band       # positive = hot, negative = cold

gdf["gi_band"] = [confidence_band(z, p)
                  for z, p in zip(gi_star.Zs, gi_star.p_sim)]

A value of +99 means a 99% confidence hot spot; -90 means a 90% confidence cold spot; 0 means not significant. Reserve the 99% band for decisions with real consequences (intervention targeting) and treat the 90% band as exploratory.

Output Interpretation #

Read the fitted result through three lenses:

• Sign of Zs: positive marks a hot spot (a neighbourhood summing above the global expectation), negative marks a cold spot. The sign is the direction of the anomaly.
• Magnitude of Zs: a z of 3.2 is a stronger, more spatially coherent cluster than a z of 2.0. Magnitude conveys intensity, not just presence.
• Corrected p_sim: only units clearing the FDR threshold should be coloured on the map. Everything else is “not significant” regardless of z-sign.

Good signs include spatially contiguous blocks of hot or cold units — real processes produce clustered hot spots, not scattered singletons. Warning signs include isolated significant units surrounded by non-significant neighbours (often an edge effect or an artefact of an over-tight distance band), and a suspiciously high fraction of significant units, which usually means multiple-comparison correction was skipped. Because hot spot geometry is sensitive to non-stationary trend in the mean, cross-check against the assumptions covered in stationarity and trend analysis: a global gradient can manufacture a hot spot at one edge and a cold spot at the other with no local clustering at all.

Production Considerations #

Complexity. G_Local with permutations is roughly $O(M \cdot n \cdot \bar{k})$, where $M$ is the permutation count and $\bar{k}$ the mean neighbour count. Runtime scales linearly in permutations, so exploratory runs can use permutations=99 and publication runs permutations=9999.

Sparsity. Always keep the weight matrix sparse (w.sparse, CSR). Never materialise the dense $n \times n$ form for large $n$; at $n = 50{,}000$ the dense float64 array alone would need roughly 20 GB.

Reproducibility. Pass an explicit seed to G_Local so the permutation null — and therefore every reported p-value and map — is deterministic across runs. Pin esda and libpysal versions, since permutation internals and default transforms have shifted between releases.

Edge effects. Units near the study-area boundary have truncated neighbourhoods and unstable statistics. Where possible, buffer the analysis with data beyond the reporting extent, or flag boundary units in the output so downstream consumers treat them cautiously.

Troubleshooting #

Next Steps #

For a copy-pasteable end-to-end workflow — building binary weights, running G_Local(y, w, star=True), extracting Zs and p_sim, and mapping graded confidence bands — follow Getis-Ord Gi* Hot Spot Analysis in Python. To contrast intensity mapping with local correlation and outlier detection, work through the spatial autocorrelation metrics guide, and when your goal shifts from finding hot spots to grouping units into homogeneous regions, see spatial clustering and regionalization.

Frequently Asked Questions #

What is the difference between Getis-Ord Gi and Gi?*

The $G_i$ statistic sums the values of an observation’s neighbours only, excluding the observation itself. $G_i^*$ additionally includes the focal observation by setting the self-weight to one. $G_i^*$ is the standard choice for hot spot mapping because a location that is itself an extreme value should contribute to whether it sits in a hot spot; $G_i$ is used when you specifically want neighbourhood context without the focal value.

How does hot spot analysis differ from Local Moran’s I (LISA)?

Local Moran’s I measures local correlation and distinguishes four categories including spatial outliers (high-low and low-high). Getis-Ord $G_i^*$ measures the local sum of values relative to the global mean and identifies only hot spots and cold spots, not outliers. Use $G_i^*$ to answer where the intense clusters of high or low values are; use LISA when you also need to flag observations that differ sharply from their neighbours.

Why do I need multiple-comparison correction for Gi?*

$G_i^*$ runs one hypothesis test per spatial unit, so a map of 3,000 units performs 3,000 tests. At an uncorrected 0.05 threshold roughly 150 units would appear significant by chance alone. A false discovery rate correction such as Benjamini-Hochberg, or a simpler Bonferroni-style bound, rescales the threshold so the map reflects genuine spatial structure rather than the expected volume of false positives.

Should Gi use binary or row-standardised weights?*

Getis-Ord statistics are defined on binary weights because $G_i^*$ is a sum of neighbour values, not a weighted average. Row-standardisation would rescale that sum and change its interpretation. Build the weight matrix with a binary transform, include the self-weight for $G_i^*$, and reserve row-standardisation for Moran’s I and spatial regression.

Related

• Getis-Ord Gi* Hot Spot Analysis in Python — the full step-by-step esda implementation with mapping
• Spatial Autocorrelation Metrics — Moran’s I, Geary’s C, and LISA, the local-correlation counterparts to Gi*
• Spatial Weight Matrices — building and validating the binary weights that Getis-Ord requires
• Spatial Clustering & Regionalization — grouping units into homogeneous regions once hot spots are identified

← Back to Core Concepts of Spatial Statistics & Geostatistics