2.4 Measuring Global Spatial Autocorrelation

Part: 2 — Spatial Autocorrelation (Chapter 2.4 of 6)

Dataset: diabetes_northeast.shp (217 counties); data_1998_2010_long_lbc.csv + US_COUNTY.shp (3,107 counties, 1998–2012)

Focus: Moran's I, Geary's C, General G, join counts, bivariate and differential Moran's I

Library: spatialstats.esda (PySpatialStats)

---

Table of Contents

1. Overview

2. Python Setup

3. Moran's I: Definition and Permutation Mechanics

   • 3.1 The permutation p-value and its floor

4. Geary's C

   • 4.1 A verified bug: the permutation p-value

5. Getis-Ord General G

6. Join Counts for Binary and Categorical Data

7. Bivariate Moran's I

8. Differential Moran's I: Change Over Time

9. Output Handling

10. Performance: Vectorised vs. Python-loop Statistics

11. Summary and Conclusions

12. Resources

---

1. Overview

Chapter 2.3's scatterplot slope becomes, here, a family of formal test statistics with null distributions, z-scores, and permutation p-values. Moran's I is the central one; the rest answer related but distinct questions.

Statistic	Question	Function
Moran's I	Overall linear association between a value and its neighbours' values	Moran
Geary's C	Overall dissimilarity between neighbouring pairs (locally sensitive alternative to Moran's I)	Geary
General G	Do high values cluster near high values specifically (not just "similar near similar")?	GeneralG
Join counts	The categorical analogue, for binary/two-class data	Join_Counts
Bivariate Moran's I	Does variable $x$ at $i$ predict variable $y$ at $i$'s neighbours?	MoranBV
Differential Moran's I	Is the change in a variable over time itself spatially autocorrelated?	Moran_Diff

Along the way this chapter surfaces a real, verified bug in this library's Geary permutation p-value, and a real ~90–100× performance gap between the vectorised Moran and the pure-Python-loop Geary/GeneralG at national scale — both practical hazards, not textbook footnotes.

---

2. Python Setup

import sys, time, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)


def find_repo_root(start: Path = Path.cwd()) -> Path:
    for p in [start, *start.parents]:
        if (p / "pyproject.toml").exists() and (p / "data").is_dir():
            return p
    raise FileNotFoundError("Run this notebook from inside the PySpatialStats repository.")


ROOT = find_repo_root()
DATA = ROOT / "data"
try:
    import spatialstats
except ImportError:
    sys.path.insert(0, str(ROOT))
    import spatialstats

from spatialstats.weights import Queen, higher_order
from spatialstats.esda import Moran, Geary, GeneralG, Join_Counts, MoranBV, Moran_Diff

plt.rcParams.update({
    "figure.dpi": 100, "savefig.dpi": 100, "font.size": 10, "axes.titlesize": 11,
    "axes.titleweight": "bold", "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "legend.frameon": False,
})
C = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
     "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

counties = gpd.read_file(DATA / "diabetes_northeast.shp")
w = Queen(counties)
print(f"spatialstats : version {spatialstats.__version__}")
print(f"counties     : {counties.shape[0]}")

spatialstats : version 0.2.0
counties     : 217

---

3. Moran's I: Definition and Permutation Mechanics

For attribute $y$, centred as $z_i = y_i - \bar y$, and row-standardised weights $W$:

$$ I \;=\; \frac{n}{S_0} \cdot \frac{\mathbf{z}^\top W \mathbf{z}}{\mathbf{z}^\top \mathbf{z}}, \qquad S_0 = \sum_{i,j} w_{ij}. $$

For row-standardised $W$, $S_0 = n$, so $I$ reduces to the correlation-like ratio $\mathbf{z}^\top W\mathbf{z} / \mathbf{z}^\top \mathbf{z}$ seen as a slope in Chapter 2.3. Under the null of no spatial association, $\mathbb{E}[I] = -1/(n-1)$ (slightly negative, not zero — a finite-sample artefact, negligible for large $n$), with a variance from the normality or randomisation assumption that gives an analytic z-score and p-value (p_norm).

res = Moran(counties["Dia_pct"], w, permutations=999, seed=1)
print(f"I           = {res.I:.4f}")
print(f"E[I]        = {res['expected_I']:.4f}   (= -1/(n-1) = {-1/(len(counties)-1):.4f})")
print(f"Var[I]      = {res['var_I']:.6f}")
print(f"z-score     = {res.z_score:.3f}")
print(f"p_norm      = {res['p_norm']:.2e}   (analytic, normal approximation)")
print(f"p_sim       = {res['p_sim']:.4f}   (999-permutation conditional randomisation)")

I           = 0.2398
E[I]        = -0.0046   (= -1/(n-1) = -0.0046)
Var[I]      = 0.001898
z-score     = 5.611
p_norm      = 2.02e-08   (analytic, normal approximation)
p_sim       = 0.0010   (999-permutation conditional randomisation)

3.1 The permutation p-value and its floor

p_sim does not come from a formula; it comes from conditional permutation: hold the observed values fixed, shuffle which county gets which value $m$ times, recompute $I$ each time, and count how many shuffled datasets produced an $I$ at least as extreme as the real one:

$$ p_\text{sim} = \frac{1 + \#\{|I_\text{sim}| \ge |I_\text{obs}|\}}{1 + m}. $$

The "+1" in numerator and denominator (the observed dataset counts as one of the $m+1$ possible arrangements) means the smallest attainable p-value is $1/(m+1)$, no matter how extreme $I$ is — with 999 permutations, that floor is $1/1000 = 0.001$.

for m in (99, 199, 999, 4999):
    res_m = Moran(counties["Dia_pct"], w, permutations=m, seed=1)
    print(f"m={m:5d} permutations  ->  p_sim = {res_m['p_sim']:.4f}   (floor = {1/(m+1):.4f})")

m=   99 permutations  ->  p_sim = 0.0100   (floor = 0.0100)
m=  199 permutations  ->  p_sim = 0.0050   (floor = 0.0050)
m=  999 permutations  ->  p_sim = 0.0010   (floor = 0.0010)

m= 4999 permutations  ->  p_sim = 0.0002   (floor = 0.0002)

Every run above hits its own floor: Dia_pct's autocorrelation is strong enough that not one of up to 4,999 random relabellings produced as extreme a statistic as the real data. Reporting "$p < 0.001$" is honest here; reporting the literal floor value as if it were a precisely estimated small probability is not — with 999 permutations, $p_\text{sim}=0.001$ and $p_\text{sim}=0.0001$ are indistinguishable, and only p_norm (which does not have this floor) can express how extreme the result is beyond that point.

---

4. Geary's C

Geary's C uses squared differences between neighbouring pairs rather than products of centred values:

$$ C = \frac{(n-1) \sum_{i,j} w_{ij} (y_i - y_j)^2}{2 \, S_0 \sum_i (y_i - \bar y)^2}. $$

$C < 1$ indicates positive autocorrelation (neighbours tend to be similar, so squared differences are smaller than under randomness); $C > 1$ indicates negative autocorrelation; $C \approx 1$ indicates none. Because it uses local pairwise differences rather than a single global product, Geary's C is more sensitive to local structure than Moran's I, at the cost of being less directly comparable to the classical correlation-coefficient interpretation.

gc = Geary(counties["Dia_pct"], w, permutations=999, seed=1)
print(f"C           = {gc.C:.4f}   (< 1 -> positive autocorrelation, consistent with Moran's I = {res.I:.3f})")
print(f"E[C]        = {gc['expected_C']:.1f}   (always 1, regardless of n or W)")
print(f"p_value     = {gc.p_value}")

C           = 0.7337   (< 1 -> positive autocorrelation, consistent with Moran's I = 0.240)
E[C]        = 1.0   (always 1, regardless of n or W)
p_value     = 1.0

4.1 A verified bug: the permutation p-value

C = 0.734 says the same thing as Moran's I: clear positive autocorrelation. But the reported permutation p-value is 1.0 — the opposite of what the effect size implies. This is not a data problem; it is a bug in how the permutation test is implemented for Geary, and it reproduces on every run.

rng = np.random.default_rng(1)
y = counties["Dia_pct"].to_numpy()

# Reproduce the actual permutation loop used inside Geary(), by hand, to see why p ends up at 1.0
from spatialstats.esda.global_ import _geary_C
obs_C = _geary_C(y, w)
sim_Cs = np.array([_geary_C(rng.permutation(y), w) for _ in range(999)])
print(f"observed C = {obs_C:.4f}")
print(f"permuted C: min={sim_Cs.min():.3f}, mean={sim_Cs.mean():.3f}, max={sim_Cs.max():.3f}")
print(f"share of permuted |C| >= |observed C|: {np.mean(np.abs(sim_Cs) >= np.abs(obs_C)):.4f}")

observed C = 0.7337
permuted C: min=0.870, mean=1.000, max=1.147
share of permuted |C| >= |observed C|: 1.0000

The permutation test applies a two-tailed, absolute-value extremity rule — |C_sim| >= |C_obs| — copied from Moran's I, where $z$ is centred at 0 and that rule is the natural one. Geary's C is centred at 1, not 0, and under positive autocorrelation it moves toward 0, i.e. its absolute value increases in the same direction as an even larger, purely random C would. Every one of the 999 permuted values above is farther from zero (in absolute terms) than the true, clustered value of 0.734 — an artefact of comparing absolute values around the wrong centre, not evidence against clustering. Do not report Geary().p_value as implemented here without correcting for this; report the value of $C$ itself (which is correctly computed) and rely on Moran's I or a corrected permutation scheme (measuring extremity as distance from 1, i.e. |C_sim - 1| >= |C_obs - 1|) for a trustworthy significance statement.

# A corrected permutation test, measuring extremity as distance from the true null centre (C = 1)
p_corrected = (1 + np.sum(np.abs(sim_Cs - 1) >= np.abs(obs_C - 1))) / (1 + len(sim_Cs))
print(f"corrected p-value (distance from 1): {p_corrected:.4f}   vs. library's p_value: {gc.p_value}")

corrected p-value (distance from 1): 0.0010   vs. library's p_value: 1.0

---

5. Getis-Ord General G

General G asks a sharper question than Moran's I: not just "are neighbours similar?" but specifically "do high values sit near high values?" (it is not symmetric to low-low clustering the way Moran's I is).

$$ G = \frac{\sum_{i \ne j} w_{ij}\, y_i y_j}{\sum_{i \ne j} y_i y_j}, \qquad y_i \ge 0 \text{ required.} $$

gg = GeneralG(counties["Dia_pct"], w, permutations=999, seed=1)
print(f"G           = {gg.G:.6f}")
print(f"p_value     = {gg.p_value}")

G           = 0.004647
p_value     = 0.007

Dia_pct is significant under General G ($p = 0.007$) but visibly less extreme than under Moran's I ($p_\text{sim}$ at its 0.001 floor) — a reminder that General G specifically tests high-value clustering, while the data also contains a substantial low-value (LL) cluster (Chapter 2.3, Section 5) that General G does not credit toward the same statistic.

---

6. Join Counts for Binary and Categorical Data

For a two-category variable, join counts classify every pair of neighbours as Black–Black (BB), White–White (WW), or Black–White (BW) — using "Black"/"White" only as historical PySAL terminology for the two categories. Fewer BW (mixed) joins than expected under random labelling is the categorical signature of positive autocorrelation. We build a binary indicator from Dia_pct (above vs. at-or-below the regional median).

above_median = (counties["Dia_pct"] > counties["Dia_pct"].median()).astype(int)
jc = Join_Counts(above_median, w, permutations=999, seed=1)
print(f"BB (high-high joins) = {jc['BB']:.0f}   p = {jc['p_BB']:.3f}")
print(f"WW (low-low joins)   = {jc['WW']:.0f}   p = {jc['p_WW']:.3f}")
print(f"BW (mixed joins)     = {jc['BW']:.0f}   p = {jc['p_BW']:.3f}")
print(f"total joins (should equal Queen's undirected edge count, {w.sparse.nnz // 2}): "
      f"{jc['BB'] + jc['WW'] + jc['BW']:.0f}")

BB (high-high joins) = 156   p = 0.022
WW (low-low joins)   = 148   p = 0.226
BW (mixed joins)     = 259   p = 0.974
total joins (should equal Queen's undirected edge count, 563): 563

BB joins (156) are significantly more common than a random above/below-median labelling would produce ($p=0.022$) — high-prevalence counties border each other more than chance. WW and BW do not reach significance individually here; with only two categories and a single realised map, join counts are typically a less powerful test than Moran's I on the original continuous variable, and are most useful when the variable is genuinely categorical (like UrbRural) rather than an artificially dichotomised continuous one.

---

7. Bivariate Moran's I

MoranBV(x, y, w) asks whether variable $x$ at region $i$ is associated with variable $y$ at $i$'s neighbours — a spatial cross-correlation, not a same-region correlation. Does obesity at county $i$ predict diabetes prevalence in the surrounding counties?

mbv = MoranBV(counties["Obesity"], counties["Dia_pct"], w, permutations=999, seed=1)
same_region = np.corrcoef(counties["Obesity"], counties["Dia_pct"])[0, 1]
print(f"MoranBV(Obesity_i, lag(Dia_pct)) = {mbv.I:.4f}   p = {mbv['p_value']:.4f}")
print(f"for comparison, ordinary (same-region) correlation(Obesity, Dia_pct) = {same_region:.4f}")

MoranBV(Obesity_i, lag(Dia_pct)) = 0.0542   p = 0.0010
for comparison, ordinary (same-region) correlation(Obesity, Dia_pct) = 0.6834

The cross-neighbour association (0.054) is far weaker than the same-region correlation (0.68) — obesity is a strong predictor of diabetes within the same county, as expected, but a county's obesity rate says comparatively little, on its own, about its neighbours' diabetes rate once the shared regional pattern in each variable separately is set aside. This distinguishes "these two variables share a spatial trend" (true — Chapter 2.1 showed both cluster regionally) from "one variable spatially spills over into the other" (weak evidence here).

---

8. Differential Moran's I: Change Over Time

Moran_Diff(y_t, y_t0, w) tests whether the change in an outcome between two periods is itself spatially autocorrelated — i.e., whether places that are near each other changed together, not just whether they started or ended similar. We use the age-adjusted mortality RATE in the long-format county panel, comparing 1998 against 2012 across all 3,107 US counties.

panel = pd.read_csv(DATA / "data_1998_2010_long_lbc.csv")
us_county = gpd.read_file(DATA / "US_COUNTY.shp")

rate_1998 = panel.loc[panel["Year"] == 1998, ["FIPS", "RATE"]].rename(columns={"RATE": "RATE_1998"})
rate_2012 = panel.loc[panel["Year"] == 2012, ["FIPS", "RATE"]].rename(columns={"RATE": "RATE_2012"})
national = us_county.merge(rate_1998, on="FIPS", how="inner").merge(rate_2012, on="FIPS", how="inner")
print(f"matched counties: {len(national)} of {len(us_county)}")

w_national = Queen(national)
print(w_national.summary())

matched counties: 3107 of 3108

Spatial Weights (W)
===================
  Observations        3107
  Transform           row
  Min neighbors       0
  Mean neighbors      5.837
  Max neighbors       14
  Islands             2
  Components          3
  Nonzero weights     18136

moran_1998 = Moran(national["RATE_1998"], w_national, permutations=99, seed=1)
moran_2012 = Moran(national["RATE_2012"], w_national, permutations=99, seed=1)
moran_diff = Moran_Diff(national["RATE_2012"], national["RATE_1998"], w_national, permutations=99, seed=1)

print(f"Moran's I, RATE 1998        : {moran_1998.I:.4f}  (p_sim={moran_1998['p_sim']:.3f})")
print(f"Moran's I, RATE 2012        : {moran_2012.I:.4f}  (p_sim={moran_2012['p_sim']:.3f})")
print(f"Moran's I, change 2012-1998 : {moran_diff.I:.4f}  (p_sim={moran_diff['p_sim']:.3f})")
print(f"mean change: {(national['RATE_2012'] - national['RATE_1998']).mean():.2f}  "
      f"(rates fell nationally over the period)")

Moran's I, RATE 1998        : 0.7310  (p_sim=0.010)
Moran's I, RATE 2012        : 0.7122  (p_sim=0.010)
Moran's I, change 2012-1998 : 0.6836  (p_sim=0.010)
mean change: -7.34  (rates fell nationally over the period)

The rate level is strongly autocorrelated in both years (I ≈ 0.71–0.73) — no surprise, geography changes slowly. What is more informative is that the change itself is also strongly autocorrelated (I = 0.684): counties did not fall (or rise) independently of their neighbours; the decline in mortality over 1998–2012 was itself a regional phenomenon, consistent with regionally-shared drivers (policy, healthcare access, economic conditions) rather than 3,107 independent local trends.

---

9. Output Handling

Every result above is a SpatialResult subclass with a uniform interface: .summary() for a formatted report, .stats for the full dictionary, .get(key) / result[key] for a single value, and .to_frame() for a one-row DataFrame (Part 0 §0.2 introduced this pattern; it holds throughout esda).

print(res.summary())
print()
res.to_frame()

Moran's I
=========
  I                    0.239826
  expected_I           -0.00462963
  var_I                0.0018983
  z_score              5.61069
  p_value              0.001
  p_norm               2.01516e-08
  p_sim                0.001
  permutations         999
  n                    217

	I	expected_I	var_I	z_score	p_value	p_norm	p_sim	permutations	n
0	0.239826	-0.00463	0.001898	5.610695	0.001	2.015161e-08	0.001	999	217

---

10. Performance: Vectorised vs. Python-loop Statistics

Moran computes its statistic with sparse matrix algebra (z @ (W @ z)); Geary and GeneralG iterate over every nonzero weight in a pure Python loop. At 217 counties the difference is invisible; at national scale (3,107 counties, mean ~5.8 Queen neighbours) it is not.

rates = national["RATE_2012"].to_numpy()
timings = {}
for name, fn in [("Moran",     lambda: Moran(rates, w_national, permutations=999, seed=1)),
                  ("Geary",     lambda: Geary(rates, w_national, permutations=999, seed=1)),
                  ("GeneralG",  lambda: GeneralG(rates, w_national, permutations=999, seed=1))]:
    t0 = time.time()
    fn()
    timings[name] = time.time() - t0

for name, secs in timings.items():
    print(f"{name:10s} (999 permutations, n=3107): {secs:6.2f} s")
print(f"\nGeary is {timings['Geary'] / timings['Moran']:.0f}x slower than Moran for the same permutation count.")

Moran      (999 permutations, n=3107):   0.12 s
Geary      (999 permutations, n=3107):  11.36 s
GeneralG   (999 permutations, n=3107):   9.96 s

Geary is 93x slower than Moran for the same permutation count.

The gap (roughly two orders of magnitude) comes entirely from the inner loop: Moran's vectorised update costs one sparse matrix–vector product per permutation, while Geary's and GeneralG's cost is proportional to the number of nonzero weights, evaluated in interpreted Python, per permutation. On a 217-county dataset this never matters; on a national or finer-grained (e.g. tract-level) dataset with several thousand permutations requested, budget for it, or reduce permutations for Geary/GeneralG specifically while keeping Moran at its default.

---

11. Summary and Conclusions

This chapter built the full global-autocorrelation toolkit and surfaced two concrete, verified implementation hazards along the way.

What we did

1. Defined Moran's I formally, connected it to Chapter 2.3's scatterplot slope, and explained the permutation p-value's floor of $1/(m+1)$.
2. Computed Geary's C (0.734, correctly indicating positive autocorrelation) and demonstrated a real bug: its reported permutation p-value (1.0) is backwards, because the extremity rule compares |C| to 0 instead of distance from Geary's true null centre of 1 — and showed a corrected version giving a sensible p-value.
3. Computed General G, finding significant high-value clustering, but a weaker signal than Moran's I — General G does not credit the data's substantial low-value cluster.
4. Ran join counts on a median-dichotomised Dia_pct, finding significant high-high clustering (BB) but not low-low or mixed joins individually.
5. Computed bivariate Moran's I between obesity and neighbouring diabetes rates, and contrasted the weak cross-neighbour association (0.054) with the strong same-region correlation (0.79).
6. Computed differential Moran's I on national mortality data (1998 vs. 2012), showing that the change in rates over time, not just their level, is itself spatially autocorrelated.
7. Reviewed the uniform .summary() / .stats / .to_frame() output interface.
8. Measured a ~90–100× runtime gap between vectorised (Moran) and Python-loop (Geary, GeneralG) statistics at national scale, and explained its source.

Practical takeaways

• Report Moran's I permutation p-values as "< floor" once they hit $1/(m+1)$, not as a precise probability; use p_norm if a finer-grained significance statement is needed.
• Do not trust Geary().p_value from this library as-is — the C statistic is correct, but its permutation p-value is not; use Moran's I, or the corrected permutation shown in Section 4.1, for inference.
• Match the number of permutations to the statistic: 999 is essentially free for Moran, but can cost seconds to tens of seconds for Geary/GeneralG at national scale.

Next steps

Chapter 2.5 (Local Spatial Autocorrelation) takes these global summaries apart into one statistic per county — Local Moran's I, with a careful look at two more permutation-inference subtleties that echo this chapter's Geary bug in spirit, if not in mechanism.

---

12. Resources

Textbooks

Reference	Notes
Cliff, A. D. & Ord, J. K. (1981). Spatial Processes: Models and Applications. Pion.	The definitive early treatment of Moran's I and Geary's C, including their moments
Getis, A. & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24, 189–206.	Origin of the General G statistic
Cliff, A. D. & Ord, J. K. (1973). Spatial Autocorrelation. Pion.	Join count statistics for categorical spatial data

Key papers

Paper	Contribution
Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37, 17–23.	Original Moran's I
Geary, R. C. (1954). The contiguity ratio and statistical mapping. The Incorporated Statistician, 5, 115–146.	Original Geary's C
Anselin, L., Syabri, I. & Smirnov, O. (2002). Visualizing multivariate spatial correlation with dynamically linked windows. In Spatial Analysis: Modelling in a GIS Environment.	Bivariate Moran's I and linked visualisation

In this library

Object	Role
spatialstats.esda.Moran	Global Moran's I: .I, .z_score, .p_norm, .p_sim, .summary(), .to_frame()
spatialstats.esda.Geary	Global Geary's C — trust .C, treat .p_value with the caveat in Section 4.1
spatialstats.esda.GeneralG	High-value-specific clustering
spatialstats.esda.Join_Counts	Categorical/binary analogue
spatialstats.esda.MoranBV, Moran_Diff	Bivariate and differential (change-over-time) Moran's I
Chapter 2.5	esda.Moran_Local, esda.Geary_Local, esda.G_Local — the pointwise versions of this chapter's statistics