Diagnosing Spatial Autocorrelation with Global Moran's I in Python

TL;DR Answer

Build a row-standardised spatial weights matrix with libpysal (Queen contiguity for polygons, KNN for points), then pass your values and the weights to esda.Moran. Read three numbers: I (the statistic, near +1 for clustering, near 0 for randomness), p_sim (the permutation p-value), and z_sim (the standardised score). Run it on your target variable first to confirm autocorrelation exists, then run it again on model residuals — a significant result there signals leftover spatial structure and model misspecification. This is the core diagnostic behind Handling Spatial Autocorrelation.

from libpysal.weights import Queen
from esda.moran import Moran

w = Queen.from_dataframe(gdf)
w.transform = "r"                      # row-standardise

mi = Moran(gdf["residual"].to_numpy(), w, permutations=999)
print(f"I={mi.I:.3f}  p_sim={mi.p_sim:.4f}  z={mi.z_sim:.2f}")

Part of: Handling Spatial Autocorrelation

What Global Moran’s I Actually Measures

Global Moran’s I is a single scalar that summarises the degree of spatial autocorrelation across an entire dataset. Conceptually it is a correlation coefficient between each observation’s value and the average value of its neighbours — the spatial lag. When high values sit next to high values and low next to low, the statistic tends toward its maximum of +1. When the map is a checkerboard of opposites, it tends toward −1. When values are scattered with no geographic pattern, it hovers near its expected value under randomness, E[I] = −1 / (n − 1), which is close to zero for any realistic sample size.

The statistic is defined over a spatial weights matrix W whose entries w_ij encode which observations are neighbours and how strongly. Everything downstream depends on that matrix: change the neighbour definition and you change the answer. That is why building W correctly is the real work, and computing I is almost an afterthought.

Two questions motivate this diagnostic in a machine-learning workflow. First, does the target variable carry spatial structure worth exploiting? A strongly autocorrelated target means nearby observations are not independent, which is exactly the condition that breaks random cross-validation and inflates your scores. Second, and more importantly, do the residuals of a fitted model still carry spatial structure? If they do, the model has failed to absorb a signal that geography could explain, and your error estimates are optimistic.

Global Versus Local: Don’t Confuse the Two

This page covers the global statistic — one number for the whole map, used as a pass/fail diagnostic. It is deliberately distinct from local Moran’s I, which returns one value per observation and pinpoints individual hot spots, cold spots, and spatial outliers. The local decomposition is a feature-engineering tool, and it has its own walkthrough in Computing Local Moran’s I for Feature Engineering. Reach for the global statistic when you want a yes/no answer about whether autocorrelation exists; reach for the local statistic when you want to know where it lives and turn that into a model input.

Building the Spatial Weights Matrix

The weights matrix is the single most consequential choice in the whole procedure, and it depends on your geometry type.

• Polygons with shared borders — census tracts, counties, catchments — call for contiguity weights. libpysal.weights.Queen.from_dataframe(gdf) treats any two polygons that touch at an edge or a corner as neighbours; Rook requires a shared edge. Queen is the common default.
• Points — soil cores, weather stations, property sales — have no shared borders, so contiguity is undefined. Use KNN.from_dataframe(gdf, k=8) to fix a neighbour count, or DistanceBand to include every point within a radius. K-nearest-neighbours guarantees every observation has neighbours even in sparse regions, which avoids the disconnected-island problem a fixed radius can create.

Whatever the construction, row-standardise by setting w.transform = "r". This rescales each row so its weights sum to one, turning the spatial lag into a simple neighbour average. Without it, observations with many neighbours dominate the statistic purely because of their connectivity, not their values. Row standardisation is what makes I comparable across a study area with uneven neighbour counts.

Coordinates must be in a projected metric CRS before you compute KNN or distance-band weights, for the same reason K-means blocking needs it — Euclidean distance on raw latitude and longitude is distorted. If your data is still in degrees, reproject it first; the mechanics are covered in CRS Alignment and Projection Handling.

A Production-Ready morans_i_report()

The function below wraps weights construction, row standardisation, the permutation test, and interpretation into one typed call. It validates inputs, aligns the values array with the weight matrix, warns about observations with no neighbours (which silently bias the statistic), and returns a structured result you can log or assert against in a pipeline.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np
from esda.moran import Moran
from libpysal.weights import KNN, Queen, W

logger = logging.getLogger(__name__)


@dataclass
class MoranResult:
    I: float
    expected_I: float
    p_sim: float
    z_sim: float
    n: int
    significant: bool
    interpretation: str


def morans_i_report(
    values: np.ndarray,
    weights: W,
    *,
    alpha: float = 0.05,
    permutations: int = 999,
    seed: int | None = 42,
) -> MoranResult:
    """Compute global Moran's I with a permutation test and interpret it.

    Parameters
    ----------
    values : np.ndarray of shape (n_samples,)
        The variable to test — a target column or model residuals.
        Order must match the observation order used to build ``weights``.
    weights : libpysal.weights.W
        A spatial weights object. Row-standardise it (``weights.transform = "r"``)
        before calling for interpretable results.
    alpha : float, default=0.05
        Significance threshold for the permutation p-value.
    permutations : int, default=999
        Conditional randomisations used to build the reference distribution.
    seed : int or None, default=42
        Seed for reproducible permutations.

    Returns
    -------
    MoranResult
    """
    values = np.asarray(values, dtype=float).ravel()
    if values.ndim != 1:
        raise ValueError("values must be one-dimensional.")
    if values.shape[0] != weights.n:
        raise ValueError(
            f"values has {values.shape[0]} rows but weights describes "
            f"{weights.n} observations — order/length must match."
        )
    if np.isnan(values).any():
        raise ValueError("values contains NaNs; drop or impute before testing.")
    if weights.islands:
        logger.warning(
            "%d observation(s) have no neighbours (islands); they bias I. "
            "Consider KNN weights or dropping them.", len(weights.islands)
        )

    if seed is not None:
        np.random.seed(seed)

    mi = Moran(values, weights, permutations=permutations)
    significant = bool(mi.p_sim < alpha)

    if not significant:
        interp = "No significant spatial autocorrelation — values look spatially random."
    elif mi.I > mi.EI:
        interp = "Positive spatial autocorrelation — similar values cluster together."
    else:
        interp = "Negative spatial autocorrelation — neighbouring values are dissimilar."

    result = MoranResult(
        I=float(mi.I),
        expected_I=float(mi.EI),
        p_sim=float(mi.p_sim),
        z_sim=float(mi.z_sim),
        n=int(weights.n),
        significant=significant,
        interpretation=interp,
    )
    logger.info(
        "Moran's I=%.4f (E=%.4f) p_sim=%.4f z=%.2f — %s",
        result.I, result.expected_I, result.p_sim, result.z_sim, result.interpretation,
    )
    return result

Step-by-Step Walkthrough

1. Test the target variable first

Before modelling, confirm the target actually has spatial structure. If it does not, spatial cross-validation and spatial features buy you little.

import geopandas as gpd
from libpysal.weights import Queen

tracts = gpd.read_file("housing_by_tract.geojson").to_crs("EPSG:32617")

w = Queen.from_dataframe(tracts)
w.transform = "r"

target = morans_i_report(tracts["median_price"].to_numpy(), w)
print(target.I, target.p_sim, target.interpretation)

A typical output for house prices might be I=0.61 p_sim=0.001 — strong, significant positive autocorrelation. Nearby tracts have similar prices, exactly as expected, and this confirms that random data splitting would leak information.

2. Fit a model and extract residuals

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

feature_cols = ["income", "rooms", "age", "distance_cbd", "school_score"]
X = tracts[feature_cols].to_numpy()
y = tracts["median_price"].to_numpy()

model = GradientBoostingRegressor(random_state=42).fit(X, y)
tracts["residual"] = y - model.predict(X)

3. Test the residuals — the diagnostic that matters

res = morans_i_report(tracts["residual"].to_numpy(), w)
print(f"Residual Moran's I = {res.I:.3f}, p_sim = {res.p_sim:.4f}")
print(res.interpretation)

If the residuals return I=0.28 p_sim=0.001, the model has left behind significant positive autocorrelation. Its errors cluster geographically: whole neighbourhoods are systematically over- or under-predicted. The features you supplied did not capture a spatially structured driver — perhaps a regional trend or an omitted neighbourhood effect. If instead you see I=0.02 p_sim=0.41, the residuals are effectively random and the model has absorbed the spatial signal.

4. Act on a significant residual result

Leftover autocorrelation points to two common remedies. Removing a large-scale gradient is handled in Removing Spatial Trend with Geographic Detrending, and turning local structure into explicit inputs is covered by Computing Local Moran’s I for Feature Engineering. Either way, a positive residual I also warns that your reported error is optimistic, which is precisely why honest evaluation needs Spatial Cross-Validation Strategies rather than random folds.

Verification: Reading the Permutation p-value Correctly

The esda.Moran p-value is not analytic — it comes from a conditional permutation test. The observed values are shuffled across locations permutations times, I is recomputed for each shuffle, and p_sim is the fraction of permuted statistics as extreme as the observed one. With 999 permutations the smallest reportable p-value is 1/1000 = 0.001, so p_sim = 0.001 means “more extreme than every shuffle,” not literally zero.

mi = Moran(tracts["residual"].to_numpy(), w, permutations=999)

print(f"Observed I     : {mi.I:.4f}")
print(f"Expected I (EI): {mi.EI:.4f}")     # ≈ -1/(n-1), near zero
print(f"p_sim          : {mi.p_sim:.4f}")
print(f"z_sim          : {mi.z_sim:.2f}")  # std. deviations above/below EI

assert mi.p_sim >= 0.05, (
    "Residuals are spatially autocorrelated — model is misspecified. "
    "Add spatial features or detrend before trusting error estimates."
)

Interpret the trio together. A small p_sim tells you the pattern is unlikely under randomness; the sign of I − EI tells you which pattern — positive for clustering, negative for a dispersed checkerboard; and z_sim reports how many standard deviations the observed I sits from its randomised mean, giving you effect magnitude. A z-score of +8 with p_sim=0.001 is a far stronger clustering signal than a borderline z of +2. Wiring the assert into a training pipeline turns the diagnostic into an automatic gate: the run fails loudly whenever residuals retain spatial structure, the same defensive posture that Reducing Spatial Leakage in Model Training applies to the split itself.

FAQ

What is the difference between global and local Moran’s I?

Global Moran’s I returns a single number that summarises whether the entire dataset is spatially clustered, dispersed, or random — it answers “is there autocorrelation anywhere in this variable?” Local Moran’s I decomposes that global statistic into one value per observation, identifying which specific locations sit inside hot spots, cold spots, or spatial outliers. Use the global statistic as a diagnostic gate; use the local statistic to engineer per-location features, as shown in Computing Local Moran’s I for Feature Engineering.

How do I interpret a significant p_sim on model residuals?

A small p_sim (for example below 0.05) with a positive Moran’s I means the residuals are spatially clustered: locations where the model over-predicts sit next to other over-predictions, and likewise for under-predictions. That is a symptom of model misspecification — the model has missed a spatially structured signal such as a trend or an omitted regional covariate. Honest cross-validation error is likely higher than your current estimate, and the fix is usually a spatial feature or detrending, not more trees.

Which spatial weights should I use, Queen contiguity or KNN?

Use Queen or Rook contiguity when your observations are polygons that share borders, such as census tracts or administrative zones. Use K-nearest-neighbours or a distance band when observations are points, such as soil samples or sensor stations, because points have no shared edges. Row-standardise the weights so each observation’s neighbours sum to one, which makes Moran’s I comparable across units with different neighbour counts.

Related

• Removing Spatial Trend with Geographic Detrending
• Reducing Spatial Leakage in Model Training
• Computing Local Moran’s I for Feature Engineering
• Spatial Cross-Validation Strategies

Part of: Handling Spatial Autocorrelation — Training Geospatial Predictive Models in Python