home / methods / hot spots

Getis-Ord Gi* & Moran's I

MESA detects where high or low values cluster using the Getis-Ord Gi* local statistic, reports global Moran's I as a diagnostic, and controls false discovery across the many local tests. Inference is by Monte-Carlo permutation and is bit-reproducible across platforms.

on this page

1Getis-Ord Gi*
2Spatial weights
3Inference
4Gi_Bin classification
5Global Moran's I
6Multiple testing
7Empirical-Bayes smoothing
8Numerics & edge cases
9Validation

01 the Gi* statistic

The Getis-Ord Gi* statistic¹² measures, for each location i, whether the values in its neighborhood sum to more (a hot spot) or less (a cold spot) than you would expect if the values were distributed at random. MESA computes the standardized form:

Getis-Ord Gi*\[ G_i^{*} = \frac{\sum_j w_{ij}\,x_j - \bar{X}\sum_j w_{ij}}{S\,\sqrt{\dfrac{n\sum_j w_{ij}^{2} - \left(\sum_j w_{ij}\right)^{2}}{n-1}}} \]

wᵢⱼspatial weight between i and j

xⱼattribute value at location j

X̄, Sglobal mean and (population-form) standard deviation of x

nnumber of features

The result is a z-score. Under the null of complete spatial randomness Gi* ~ N(0,1), and a two-tailed p-value is erfc(|z| / √2). Mean and variance use the population form (divided by n) because the (n−1) correction already appears in the denominator, and both are accumulated with Welford's numerically stable algorithm.

The choice between Gi* and Gi is whether the focal location is counted. By default the neighbor sum includes location i itself (its own weight wᵢᵢ), which is the standard Gi* used by ArcGIS and most applied work. Excluding the focal point gives the original Gi statistic of Getis & Ord (1992), computed leave-one-out with (n−1)/(n−2) scaling to match PySAL G_Local(star=False) and R spdep::localG.

02 spatial weights

Neighbors are defined by one of two models (contiguity is planned but not yet in the engine):

fixed distanceall features within a distance threshold are neighbors; symmetric, but can leave isolated features in sparse rural areas

k-nearest (k=8)the k closest features; the graph is symmetrized by union (if A→B or B→A, both are neighbors), following the PySAL convention

Weights are then either binary (1 for a neighbor) or, by default, row-standardized (wᵢⱼ = 1/nᵢ, so each feature's weights sum to 1 and Gi* becomes a weighted average of neighbor values, the ArcGIS default). Distance is Euclidean for projected coordinates or ellipsoidal geodesic on WGS84; MESA hard-rejects Euclidean distance on unprojected geographic coordinates, which would otherwise introduce silent, large distortion. Two documented presets bundle these choices: an ArcGIS-compatible preset (fixed distance, row-standardized, Gi*) and an epidemiology preset (k-nearest with k=8, binary weights, minimum neighbor and event counts).

03 inference

MESA offers an analytical z-score (using the randomization-form variance) and, by default, Monte-Carlo permutation inference. Permuting the observed values many times and recomputing Gi* builds an empirical null against which the observed statistic is judged, which avoids leaning on the normal approximation in the tails.

Two permutation schemes are available for the local statistic:

totalevery value is shuffled across all locations (MESA's default)

conditionalthe focal value xᵢ is held fixed and the other n−1 values are permuted (Anselin 1995; the LISA convention in PySAL)⁴

The permutation seed is fixed and the reductions are order-independent, so a given input and seed produce identical p-values on every platform and core count (see reproducibility). A permutation count of 999 is the documented recommendation; the analytical form is available as a fallback.

04 the Gi_Bin classification

MESA reduces each location's signed significance to a seven-class Gi_Bin, the same scheme ArcGIS uses. When a multiple-testing correction is applied, the adjusted p-value drives the classification.

bin	class	rule
+3 / −3	hot / cold, 99%	p ≤ 0.01
+2 / −2	hot / cold, 95%	p ≤ 0.05
+1 / −1	hot / cold, 90%	p ≤ 0.10
0	not significant	p > 0.10

The sign is the sign of the z-score: positive means local values are higher than expected (hot), negative means lower (cold). A location with undefined statistics (NaN z or p, from insufficient variance or no neighbors) is flagged as undefined rather than forced into a class.

05 global Moran's I

Global Moran's I³ summarizes autocorrelation across the whole study area in a single number:

Moran's I\[ I = \frac{n}{\sum_{ij} w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij}\,(x_i - \bar{X})(x_j - \bar{X})}{\sum_i (x_i - \bar{X})^{2}} \]

Its expectation under the null is E[I] = −1/(n−1), so values above that indicate positive autocorrelation (similar values cluster) and below it, negative. All sums use Kahan compensated summation, and inference is by total-randomization permutation.

Moran's I describes; it doesn't decide. MESA reports it as an informational diagnostic. Significant global autocorrelation does not say where clustering is, and its absence does not preclude meaningful local clusters. Use it to characterize the surface, then use Gi* or the scan to locate clusters.

06 multiple testing

Computing a local test at every location means thousands of hypotheses, so raw p-values overstate significance. By default MESA applies the Benjamini-Hochberg procedure⁵ to control the false discovery rate, following its established use for local statistics of spatial association⁶.

Benjamini-Hochberg adjusted p\[ \tilde{p}_{(i)} = \min\!\left( 1,\ \min_{j \ge i} \frac{n}{j}\, p_{(j)} \right), \qquad p_{(1)} \le \cdots \le p_{(n)} \]

A Bonferroni option (p̃ᵢ = min(1, n·pᵢ)) is available for family-wise error control when a conservative bound is wanted. The correction matches R stats::p.adjust to 10⁻⁶. Benjamini-Yekutieli is discussed in MESA's internal methodology notes but is not implemented in the shipped engine, so it is not offered here.

07 empirical-Bayes rate smoothing

Rates computed from small denominators are unstable: one case in a village of thirty reads as an enormous rate. For rate data MESA provides a global empirical-Bayes smoother⁷ (the estimator used by GeoDa and PySAL, assuming a Poisson-Gamma prior⁸), which shrinks each local rate toward the global mean in proportion to how unstable it is:

shrinkage\[ \begin{gathered} \text{smoothed}_i = w_i\,r_i + (1 - w_i)\,m, \qquad w_i = \frac{a}{a+b} \\ a = \text{prior variance of the rates}, \quad b = m / \text{population}_i, \quad m = \text{global rate} \end{gathered} \]

An area with a large population gets a weight near 1 (its own rate is trusted); a tiny population gets a weight near 0 (it is pulled toward the global rate). The prior variance is a population-weighted variance of the raw rates, matching R spdep::EBest. Smoothing is applied to rate inputs; it is a rate-preparation step, not a change to the Gi* or scan statistics themselves.

08 numerics & edge cases

Degenerate inputs are handled explicitly rather than left to crash or return a misleading result:

input	behavior
zero variance (all values identical)	returns z = 0, p = 1 with a warning, rather than dividing by zero
isolated feature (no neighbors)	returns z = 0, p = 1 with a warning
duplicate coordinates	distances of exactly 0 are handled without singular matrices
NaN or infinite values	propagate under IEEE-754 and are classified as undefined, not crashed on
empty dataset	returns gracefully with a numerical-instability warning

The analysis functions keep no global mutable state, so they are thread-safe and re-entrant, and results are bit-for-bit identical from 1 to 128 cores.

09 validation

The Gi* engine is cross-validated against PySAL esda.G_Local: z-scores agree within 10⁻⁴, FDR-adjusted p-values within 10⁻³, and the Gi_Bin classification matches at the 90/95/99% levels wherever the neighbor graphs agree. (Where they differ, it is a documented tie-breaking difference: PySAL resolves exact distance ties by KD-tree traversal order, which varies by architecture, while MESA breaks ties deterministically by feature id.)

The engine is also calibrated for statistical validity. Under simulated complete spatial randomness the empirical false-positive rate tracks the nominal level (≈ 0.05 at α = 0.05, and likewise at 0.01 and 0.001), asserted against 99.999% Wilson score intervals, which confirms the Type-I error is controlled. To check power, a 3×3 hot spot injected into a grid produces a true-positive rate that rises monotonically with effect size, exceeding 90% at large effect sizes while noise is correctly suppressed.

§ references

Getis, A., & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206.
Ord, J. K., & Getis, A. (1995). Local spatial autocorrelation statistics: distributional issues and an application. Geographical Analysis, 27(4), 286-306.
Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37(1/2), 17-23.
Anselin, L. (1995). Local indicators of spatial association (LISA). Geographical Analysis, 27(2), 93-115.
Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B, 57(1), 289-300.
Caldas de Castro, M., & Singer, B. H. (2006). Controlling the false discovery rate: a new application to account for multiple and dependent tests in local statistics of spatial association. Geographical Analysis, 38(2), 180-208.
Marshall, R. J. (1991). Mapping disease and mortality rates using empirical Bayes estimators. Journal of the Royal Statistical Society: Series C (Applied Statistics), 40(2), 283-294.
Clayton, D., & Kaldor, J. (1987). Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics, 43(3), 671-681.