Geostatistical Modeling Spatial Statistics with Python

Getis-Ord Gi* vs Local Moran's I for Hot Spots

TL;DR: Use esda.getisord.G_Local(y, w, star=True) when the question is where is the intensity high or low — Gi* answers it directly with a signed z-score. Use esda.moran.Moran_Local(y, w) when spatial outliers matter, because Local Moran’s I splits significant units into high-high, low-low, high-low and low-high, and Gi* structurally cannot see the last two. Run both on one weights object; disagreement is information, not error.

Why This Matters

Both statistics get called “hot-spot analysis” and both produce a map of significant units, which makes them look interchangeable until a stakeholder asks why an obviously extreme location is not highlighted. They are testing different null hypotheses. Choosing the wrong one either hides the anomalies you were hired to find or clutters an intensity map with single-unit spikes. This choice sits inside hot spot analysis and rests on the neighbour definitions covered in spatial weight matrices.

The formulas make the difference plain. Gi* compares a neighbourhood sum against the global total:

Gi∗=∑jwijxj−xˉ∑jwijS[n∑jwij2−(∑jwij)2n−1]1/2,G_i^* = \frac{\sum_j w_{ij} x_j - \bar{x}\sum_j w_{ij}}{S \left[\dfrac{n\sum_j w_{ij}^2 - \left(\sum_j w_{ij}\right)^2}{n-1}\right]^{1/2}} ,

with jj including ii itself. Local Moran’s I instead multiplies the focal deviation by the neighbouring deviations:

Ii=xi−xˉm2∑j≠iwij(xj−xˉ),m2=∑i(xi−xˉ)2n.I_i = \frac{x_i - \bar{x}}{m_2} \sum_{j \ne i} w_{ij} (x_j - \bar{x}), \qquad m_2 = \frac{\sum_i (x_i - \bar{x})^2}{n} .

A sum cannot change sign when one term is high and the rest are low; a product can. That single algebraic fact is the whole of the difference.

How the two statistics classify the same four neighbourhoods Four three-by-three neighbourhoods. In the first, a high centre sits among high neighbours: both statistics agree it is a hot spot. In the second, a low centre sits among low neighbours: Gi star reports a cold spot and Local Moran's I reports low-low. In the third, a high centre sits among low neighbours: Gi star returns an unremarkable score while Local Moran's I flags a high-low outlier. In the fourth, a low centre among high neighbours: Gi star again returns nothing while Local Moran's I flags low-high. The same neighbourhood, two verdicts shaded cells are high values, open cells are low high among high Gi*: hot spot LISA: high-high they agree low among low Gi*: cold spot LISA: low-low they agree high among low Gi*: not significant LISA: high-low outlier only LISA sees it low among high Gi*: not significant LISA: low-high outlier only LISA sees it why Gi* sums the values LISA multiplies deviations only a product A sum cannot flip sign when one term opposes the rest — which is exactly why Gi* has no outlier category

Environment and Version Pinning

bash
pip install "geopandas>=1.0" "libpysal>=4.9" "esda>=2.5" "numpy>=1.24" "pandas>=2.0"
python
import geopandas as gpd
import numpy as np
import pandas as pd
import libpysal
from esda.getisord import G_Local
from esda.moran import Moran_Local
from shapely.geometry import box

Step-by-Step Implementation

1. Build one weights object for both statistics

Any divergence between the two maps must be attributable to the statistics, not the neighbourhoods. Build w once.

python
cells = [box(x, y, x + 1, y + 1) for y in range(20) for x in range(20)]
gdf = gpd.GeoDataFrame(geometry=cells, crs="EPSG:32618").reset_index(drop=True)

rng = np.random.default_rng(11)
cx = gdf.geometry.centroid.x.to_numpy()
cy = gdf.geometry.centroid.y.to_numpy()

# Two broad clusters plus deliberately planted single-cell spikes.
field = (np.exp(-((cx - 5) ** 2 + (cy - 5) ** 2) / 12)
         + np.exp(-((cx - 15) ** 2 + (cy - 15) ** 2) / 12))
gdf["rate"] = 20 + 40 * field + rng.normal(0, 2.0, len(gdf))
gdf.loc[[87, 152, 233], "rate"] += 45      # high values in low neighbourhoods
gdf.loc[[46, 118], "rate"] -= 28           # low values in high neighbourhoods

w = libpysal.weights.Queen.from_dataframe(gdf, use_index=False)
w.transform = "R"
y = gdf["rate"].to_numpy()

2. Run Gi* with the self-inclusive star

python
gi = G_Local(y, w, star=True, permutations=999, seed=7)
gdf["gi_z"] = gi.Zs
gdf["gi_p"] = gi.p_sim

ALPHA = 0.05
gdf["gi_class"] = "not significant"
gdf.loc[(gdf.gi_p < ALPHA) & (gdf.gi_z > 0), "gi_class"] = "hot spot"
gdf.loc[(gdf.gi_p < ALPHA) & (gdf.gi_z < 0), "gi_class"] = "cold spot"
print(gdf["gi_class"].value_counts())
text
not significant    312
hot spot            47
cold spot           41

3. Run Local Moran’s I on the identical inputs

python
lm = Moran_Local(y, w, permutations=999, seed=7)
LABELS = {0: "not significant", 1: "high-high", 2: "low-high",
          3: "low-low", 4: "high-low"}
q = np.where(lm.p_sim < ALPHA, lm.q, 0)
gdf["lisa_class"] = pd.Series(q).map(LABELS).to_numpy()
print(gdf["lisa_class"].value_counts())
text
not significant    306
high-high           44
low-low             39
high-low             7
low-high             4

4. Cross-tabulate to find the disagreements

python
xtab = pd.crosstab(gdf["gi_class"], gdf["lisa_class"])
print(xtab.to_string())
text
lisa_class       high-high  high-low  low-high  low-low  not significant
gi_class
cold spot                0         0         3       38                0
hot spot                43         1         0        0                3
not significant          1         6         1        1              304

The interesting cell is the bottom row: six units that Local Moran’s I flags as high-low outliers, Gi* leaves unclassified. Those are the deliberately planted spikes. Gi* is not wrong — a lone high value genuinely does not make its neighbourhood’s total unusual — but if the analysis is about detecting anomalies, its silence is a failure to answer the question asked.

Note also the three units where Gi* says hot and Local Moran’s I says nothing: neighbourhoods with consistently mild elevation that sum to a significant total without any single deviation being large enough to produce a significant product.

Where the two classifications agree and diverge A matrix with Getis-Ord Gi star classes as rows and Local Moran's I classes as columns. The diagonal cells for hot spot with high-high and cold spot with low-low are heavily populated and shaded as agreement. Two off-diagonal regions are highlighted: Local Moran's I outliers that Gi star calls not significant, and Gi star hot spots that Local Moran's I calls not significant. Agreement is the diagonal; the value is in the corners high-highlow-low high-lowlow-highnot sig. Local Moran’s I Getis-Ord Gi* hot spotcold spotnot significant 43 38 010 00 11 3 6 1 outliers invisible to Gi* 3 mild but consistent elevation Gi* sums up 304 Publishing the cross-tabulation alongside the map turns a methodological choice into a visible, checkable one

Interpreting the Output

gi.Zs is a signed z-score: positive means the neighbourhood sums higher than a random neighbourhood of the same size would, negative means lower. It has no outlier category by construction, so the map divides into hot, cold and unremarkable.

lm.q is a quadrant code, and it is only meaningful where lm.p_sim is significant — which is why the masking step above sets non-significant units to 0 before mapping labels. Quadrants 1 and 3 are clusters; 2 and 4 are spatial outliers, the units whose value opposes their surroundings.

Where each statistic draws its significance boundary A plane with each cell's own standardised value on the horizontal axis and its neighbours' mean on the vertical axis, divided into the four Local Moran quadrants. The Getis-Ord Gi star cut-offs appear as two dashed, almost horizontal lines, because Gi star is a sum dominated by the eight neighbours. The Local Moran cut-offs appear as hyperbolas, because the statistic is a product. In the high-low and low-high quadrants the hyperbola runs closer to the axis than the dashed line does, leaving two shaded lens-shaped regions that Local Moran's I calls significant while Gi star does not; six high-low units and one low-high unit from the cross-tabulation sit inside them. A sum draws straight cut-offs; a product draws hyperbolas the same 400 cells, plotted as own value against neighbours’ mean −2 +2 +1 −1 the cell’s own standardised value z spatial lag Wz: neighbours’ mean 6 units here – Gi* stays silent 1 unit, same story low-high high-high 43 units – both agree low-low 38 units – both agree high-low Reading the plane Gi* cut-off: a sum → near-flat line LISA cut-off: a product → hyperbola significant to LISA, invisible to Gi* With eight queen neighbours, Gi* is dominated by the neighbours’ mean, so its cut-off barely tilts with z. The product z × Wz, though, grows large even where Wz stays close to zero. Those lenses are what a sum can never reach: they hold 6 of the 7 high-low units and 1 of the 4 low-high units from the cross-tabulation.

Read the two together. Agreement on the clusters is expected and reassuring. Where Gi* is silent and Local Moran’s I flags an outlier, you have found a place that differs from its context — often a data error, a boundary artefact, or a genuinely singular site. Where Gi* flags a hot spot and Local Moran’s I does not, you have found a broad, gentle elevation rather than a sharp cluster. Neither pattern is noise; they are different phenomena.

Critical Best Practices

Use star=True for hot-spot maps

Without the star, the focal unit is excluded from its own neighbourhood sum, and the map answers “how unusual are my neighbours” rather than “how unusual is this place”. The non-star form has legitimate uses in spillover analysis, but a hot-spot map built on it will look subtly wrong to anyone who checks a specific location.

Correct for multiple comparisons, and say which correction

Testing 400 units at α=0.05\alpha = 0.05 yields about 20 false positives before any real signal. Apply an FDR correction, or tighten to 0.01. The extent of the coloured area on the map changes materially with this choice, so it belongs in the caption, not in a methods appendix.

Keep the weights identical, and justify them once

Both statistics inherit their notion of “local” entirely from w. A distance band of 2 km and a k-nearest of 8 will produce different maps from the same data with either statistic. Choose the neighbourhood from the process scale — see spatial weight matrices — and apply it to both.

Standardise rates before comparing regions of unequal population

Both statistics are sensitive to the variance instability of raw rates computed over small denominators. A cell with three cases out of five people will dominate any neighbourhood it sits in. Use an empirical Bayes smoothed rate, or restrict the analysis to units meeting a minimum denominator.

Report both when the audience is mixed

The cross-tabulation above takes one line of code and pre-empts the question “why isn’t this location highlighted”. Publishing both classifications, with the disagreements named, is more defensible than picking one and hoping nobody checks.

Troubleshooting

Symptom Likely cause Fix
Gi* highlights nothing at all Variable has little spatial structure, or weights are too broad Check global Moran’s I first; shrink the distance band or reduce k
Every unit is significant No multiple-comparison correction with many units Apply FDR, or use permutations=9999 with α=0.01\alpha = 0.01
lm.q values look wrong Quadrants read without masking on p_sim Set non-significant units to 0 before mapping labels
Hot spots trace population density, not the phenomenon Raw counts used instead of rates Convert to a rate, then smooth small-denominator units
G_Local raises on islands in w Units with no neighbours Rebuild with KNN, or drop islands and note the exclusion
Results shift between runs Permutation seed not fixed Pass seed= to both estimators for reproducible pseudo p-values

Next Steps

For the full Gi* workflow including mapping and rate smoothing, see Getis-Ord Gi* hot spot analysis in Python; for the cluster-and-outlier map in depth, see computing Local Moran’s I (LISA) in Python. When hot spots need to be tracked over time rather than mapped once, continue to emerging hot spot analysis over time in Python.

Frequently Asked Questions

Why does Getis-Ord Gi never identify spatial outliers?*

Gi* is a ratio of a neighbourhood sum to the global sum, so it only measures whether a neighbourhood totals more or less than expected. A single very high value surrounded by low values produces a middling sum, and Gi* returns an unremarkable z-score. Local Moran’s I multiplies the unit’s own deviation by its neighbours’ average deviation, so opposite signs give a negative product that is flagged as a high-low or low-high outlier. Detecting outliers requires that cross-product; a sum cannot do it.

Should star be True or False in G_Local?

Use star=True for hot-spot mapping. The star variant includes the unit’s own value in its neighbourhood sum, which is what makes Gi* a statement about a place including itself rather than about its surroundings only. The non-star Gi excludes the focal value and answers a different question — how unusual are my neighbours — which is rarely what a hot-spot map is meant to show.

Do I need to correct for multiple comparisons?

Yes, for any map with more than a handful of units. Both statistics test every unit, so at a 0.05 threshold roughly one in twenty units will look significant by chance alone. Apply a false discovery rate correction to the pseudo p-values, or use the conditional permutation p-values with a stricter threshold such as 0.01. Report which correction was used, because the visible extent of hot spots depends heavily on it.

Related

← Back to Hot Spot Analysis

Related

Emerging Hot Spot Analysis Over Time in Python Getis-Ord Gi* Hot Spot Analysis in Python