Spatial Autocorrelation Explained

Problem statement

Tobler's first law — "near things are more related than distant things" — is usually quoted as an observation. It is better understood as a warning, because almost every statistical method you were taught assumes the opposite.

Ordinary regression assumes independent residuals. Confidence intervals assume independent observations. A t-test on 400 spatial units with strongly correlated values is not really working with 400 independent pieces of information — it might be closer to 40 — and its p-values are correspondingly too small.

Spatial autocorrelation is how you measure that dependence. Moran's I is the standard statistic, and it has a property that surprises people:

A perfect checkerboard — maximally dispersed values
  Rook weights (4 edge neighbours)    I = -0.9980   p = 0.0010
  Queen weights (8 neighbours)        I = -0.0386   p = 0.0700

Identical data. One weights definition says "the most dispersed pattern possible". The other says "indistinguishable from random". Neither is a bug.

Quick answer

Moran's I compares each value to the average of its neighbours:

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

w = Queen.from_dataframe(zones, use_index=True)
w.transform = "r"                                  # row-standardise

moran = Moran(zones["value"].to_numpy(), w, permutations=999)
print(f"I = {moran.I:+.4f}  E[I] = {moran.EI:+.4f}  "
      f"z = {moran.z_sim:+.2f}  p = {moran.p_sim:.4f}")

I = +0.9046  E[I] = -0.0025  z = +34.52  p = 0.0010

Value of I	Meaning
near +1	strong clustering — similar values adjacent
near E[I] = −1/(n−1)	no spatial pattern
near −1	dispersion — dissimilar values adjacent (a checkerboard)

E[I] is not zero. Under the null hypothesis it is −1/(n−1), which for 400 units is −0.0025. Comparing I to zero rather than to E[I] is a small error at n=400 and a large one at n=20.

Moran's I from minus one through the expected value near zero to plus one, with clustered, random and checkerboard patterns illustrated. — The midpoint is E[I] = −1/(n−1), not zero. With small n that distinction matters.

Step-by-step solution

1. Understand what the statistic actually compares

Moran's I relates each unit's value to the weighted mean of its neighbours' values. Precisely, with row-standardised weights it is the regression slope of the spatial lag on the standardised value:

import numpy as np

values = zones["value"].to_numpy()
z = (values - values.mean()) / values.std()
lag = w.sparse @ z                       # weighted mean of neighbours' z-scores

print(f"Moran's I             {moran.I:+.6f}")
print(f"OLS slope of lag on z {np.polyfit(z, lag, 1)[0]:+.6f}")
print(f"z'Wz / z'z            {(z @ lag) / (z @ z):+.6f}")
print(f"correlation(z, lag)   {np.corrcoef(z, lag)[0, 1]:+.6f}")

Moran's I             +0.904632
OLS slope of lag on z +0.904632
z'Wz / z'z            +0.904632
correlation(z, lag)   +0.958867

The first three agree exactly. The correlation does not — it is 0.959 against a slope of 0.905, because the spatial lag is an average and therefore less variable than the values it averages (sd(lag)/sd(z) = 0.943 here, and 0.958867 × 0.943 = 0.904632).

So I is the slope of the Moran scatterplot, not its correlation coefficient. The distinction is worth keeping straight: averaging always shrinks variance, so I is systematically smaller in magnitude than the correlation, and quoting one for the other overstates the pattern.

2. Recognise that the weights are half the analysis

The w object defines what "neighbour" means, and different definitions answer different questions. On the same clustered surface:

  Queen     mean nb  7.4  I = +0.9046  p = 0.0010
  Rook      mean nb  3.8  I = +0.9152  p = 0.0010
  KNN k=4   mean nb  4.0  I = +0.9167  p = 0.0010
  KNN k=8   mean nb  8.0  I = +0.9005  p = 0.0010

For a smooth field the choice hardly matters — every definition sees the same gradient. But for a pattern whose structure operates at the scale of the neighbourhood definition, it decides the answer entirely:

A perfect checkerboard
  Rook  (4 edge neighbours)  I = -0.9980  z = -28.28  p = 0.0010
  Queen (8 neighbours)       I = -0.0386  z =  -1.39  p = 0.0700

Rook weights see only the four edge-sharing neighbours, which on a checkerboard are all the opposite colour — perfect dispersion. Queen weights add the four diagonals, which are all the same colour, and the two effects cancel almost exactly.

Choose the weights from the process, and say which you used. A Moran's I quoted without its weights specification is not reproducible.

3. Row-standardise, and know what it does

w.transform = "r"

Row standardisation makes each unit's weights sum to 1, so the spatial lag is a weighted mean of neighbours rather than a sum. Without it, a unit with twelve neighbours contributes three times as much as one with four, purely because of how the boundaries fall.

This matters most for irregular administrative units, where neighbour counts vary enormously. For a regular grid it changes little.

4. Use the permutation p-value

moran = Moran(values, w, permutations=999)
print(f"analytical p = {moran.p_norm:.4f}   permutation p = {moran.p_sim:.4f}")

analytical p = 0.0000   permutation p = 0.0010

The analytical p-value assumes normality. The permutation ("conditional randomisation") p-value shuffles the values across the units many times and asks how often a random arrangement produces an I as extreme as the observed one. It makes no distributional assumption and is the right default.

With 999 permutations the smallest achievable p is 0.001 — so p = 0.0010 means "never matched in 999 tries", not "p equals exactly one in a thousand".

5. Interpret a non-significant result correctly

  random     I = +0.0260  E[I] = -0.0025  z = +1.08  p = 0.1450

That is the correct answer for data with no spatial structure. Note that I is not zero — it is +0.026, because any finite random arrangement has some incidental pattern. The z-score of 1.08 is what tells you it is unremarkable.

Reporting I alone, without the z-score or p-value, invites the reader to see structure in the second decimal place.

A checkerboard grid with rook neighbours all the opposite colour and queen neighbours mixed, producing Moran's I of minus 0.998 and minus 0.039 respectively. — Rook sees four opposite-coloured neighbours. Queen adds four of the same colour, and they cancel.

Code examples

Example 1 — global Moran's I with the choices made explicit

import geopandas as gpd
import numpy as np
import pandas as pd
from esda.moran import Moran
from libpysal.weights import KNN, Queen, Rook


WEIGHT_BUILDERS = {
    "queen": lambda gdf, **kw: Queen.from_dataframe(gdf, use_index=True),
    "rook": lambda gdf, **kw: Rook.from_dataframe(gdf, use_index=True),
    "knn4": lambda gdf, **kw: KNN.from_dataframe(gdf, k=4),
    "knn8": lambda gdf, **kw: KNN.from_dataframe(gdf, k=8),
}


def global_moran(zones, column, scheme="queen", *, permutations=999, transform="r"):
    w = WEIGHT_BUILDERS[scheme](zones)
    w.transform = transform
    if w.islands:
        print(f"  warning: {len(w.islands)} unit(s) with no neighbours were excluded")

    values = zones[column].to_numpy()
    moran = Moran(values, w, permutations=permutations)

    print(f"{scheme:6} n={w.n} mean_nb={w.mean_neighbors:4.1f}  "
          f"I={moran.I:+.4f} E[I]={moran.EI:+.4f} z={moran.z_sim:+7.2f} p={moran.p_sim:.4f}")
    return moran, w


def weights_sensitivity(zones, column):
    rows = []
    for scheme in WEIGHT_BUILDERS:
        moran, w = global_moran(zones, column, scheme)
        rows.append({"scheme": scheme, "mean_nb": round(w.mean_neighbors, 1),
                     "I": round(moran.I, 4), "p": moran.p_sim})
    frame = pd.DataFrame(rows)
    print(f"\nI ranges {frame['I'].min():+.4f} to {frame['I'].max():+.4f} across weights")
    return frame


weights_sensitivity(zones, "clustered")

queen  n=400 mean_nb= 7.4  I=+0.9046 E[I]=-0.0025 z= +34.52 p=0.0010
rook   n=400 mean_nb= 3.8  I=+0.9152 E[I]=-0.0025 z= +34.19 p=0.0010
knn4   n=400 mean_nb= 4.0  I=+0.9167 E[I]=-0.0025 z= +36.02 p=0.0010
knn8   n=400 mean_nb= 8.0  I=+0.9005 E[I]=-0.0025 z= +40.13 p=0.0010

I ranges +0.9005 to +0.9167 across weights

A spread of 0.016 across four definitions is a robust result. Run the same function on the checkerboard and the spread is 0.96 — which is the signal that the weights, not the data, are driving the answer.

Example 2 — the Moran scatterplot, which is where the interpretation lives

import matplotlib.pyplot as plt


def moran_scatter(zones, column, w, moran, ax=None):
    values = zones[column].to_numpy()
    z = (values - values.mean()) / values.std()
    lag = w.sparse @ z

    ax = ax or plt.subplots(figsize=(6, 6))[1]
    ax.axhline(0, color="#94a3b8", lw=1)
    ax.axvline(0, color="#94a3b8", lw=1)
    ax.scatter(z, lag, s=14, alpha=0.6, color="#0ea5e9")

    line = np.linspace(z.min(), z.max(), 10)
    ax.plot(line, moran.I * line, color="#ef4444", lw=2,
            label=f"slope = I = {moran.I:.3f}")

    for x, y, label in [(0.6, 0.9, "HH"), (-0.9, 0.9, "LH"),
                        (-0.9, -0.9, "LL"), (0.6, -0.9, "HL")]:
        ax.text(x * z.max(), y * lag.max(), label, fontsize=13,
                weight="bold", color="#64748b")

    ax.set_xlabel(f"{column} (standardised)")
    ax.set_ylabel("spatial lag of neighbours")
    ax.legend()
    return ax


moran, w = global_moran(zones, "clustered")
moran_scatter(zones, "clustered", w, moran)

The four quadrants are the whole vocabulary of local spatial statistics:

HH — high value, high neighbours: a hotspot
LL — low value, low neighbours: a coldspot
HL and LH — spatial outliers: a unit unlike its surroundings

The global I is the slope of the fitted line — one number covering four distinct kinds of behaviour, which is exactly why local statistics exist.

Example 3 — from global to local

from esda.moran import Moran_Local

def local_moran(zones, column, w, *, permutations=999, alpha=0.05):
    lisa = Moran_Local(zones[column].to_numpy(), w, permutations=permutations)

    labels = np.array(["not significant", "HH", "LH", "LL", "HL"])
    quadrant = np.where(lisa.p_sim <= alpha, lisa.q, 0)

    out = zones.copy()
    out["lisa_I"] = lisa.Is
    out["lisa_p"] = lisa.p_sim
    out["cluster"] = labels[quadrant]

    print(out["cluster"].value_counts().to_string())
    return out


result = local_moran(zones, "clustered", w)

cluster
not significant    173
LL                 116
HH                 111

One hundred and eleven high-high units and one hundred and sixteen low-low — the two halves of a smooth gradient. There are no significant LH or HL units at all, which is itself informative: this surface is smooth, so no unit is unlike its surroundings.

On real data the spatial outliers (LH and HL) are often the most interesting rows in the table. A low-value unit surrounded by high-value ones is precisely the anomaly a global statistic averages away.

One caution: this runs a hypothesis test for every unit, so with 400 units and α = 0.05 you would expect about 20 false positives by chance. The multiple-comparisons correction matters, and Moran_Local exposes p_z_sim and FDR-adjusted options for it.

Explanation

Why E[I] is negative

Under the null hypothesis of no spatial pattern, the expected value of Moran's I is:

E[I] = −1 / (n − 1)

Not zero. The reason is that I is computed from deviations about the sample mean, and each unit's own value contributes to that mean — so the deviations are very slightly negatively correlated with each other by construction.

For 400 units, E[I] = −0.0025 and the distinction is cosmetic. For 20 units it is −0.053, and an observed I of −0.04 is positive autocorrelation relative to the null, despite being a negative number. Always compare against moran.EI, never against zero.

Why the weights are not a technicality

The definition of "neighbour" is the model of the spatial process. Contiguity says influence stops at a shared boundary; k-nearest-neighbours says every unit has the same number of influences regardless of geography; a distance band says influence has a range in metres.

The checkerboard result shows what happens when the definition and the process operate at different scales. And there is a second, quieter problem: contiguity weights on administrative units inherit all the arbitrariness of those units. Redraw the boundaries and the adjacency graph changes, and so does I.

The practical discipline: choose the definition from the process, test two or three alternatives, and report both the choice and the spread.

Why autocorrelation invalidates ordinary statistics

Independent observations each contribute one unit of information. Correlated observations do not — if you know a unit's neighbours, you already know most of what that unit will say.

The result is an effective sample size smaller than n. Standard errors computed as if all n observations were independent are too small, confidence intervals are too narrow, and p-values are too optimistic. It is not a small effect: with I near 0.9, most of the apparent information in 400 units is redundant.

This is why spatial regression exists. If a regression's residuals show significant Moran's I, the model is missing spatial structure, and the fix is a spatial lag or spatial error specification rather than tighter confidence intervals.

A Moran scatterplot with four quadrants labelled HH, LH, LL and HL, and the fitted slope equal to Moran's I. — The global I is one slope through four different behaviours. Local statistics split them apart again.

Why a global statistic can hide everything interesting

Moran's I is one number for the whole map. A study area with a strong hotspot in one corner and strong dispersion elsewhere can average to an I near zero — "no spatial pattern" — while containing two very strong ones.

That is the entire motivation for local indicators (LISA) and Getis-Ord Gi*: compute the statistic per unit, and map where the pattern actually is. Use the global statistic to answer "is there structure here at all", and local statistics to answer "where".

Edge cases or notes

Islands — units with no neighbours — get zero weight and are effectively dropped, with a warning. See spatial weights warn about islands.
use_index=True is required by current libpysal versions when building weights from a GeoDataFrame; without it you get a deprecation warning and index-alignment risk.
Row-standardisation changes I slightly and changes the interpretation from "sum of neighbours" to "mean of neighbours". Row-standardise for irregular units.
999 permutations gives a minimum p of 0.001. For stricter thresholds use more permutations.
Rates need care. Moran's I on a rate with wildly varying denominators picks up small-population noise. Use esda.moran.Moran_Rate, which applies an Empirical Bayes correction.
Moran's I on regression residuals is the standard diagnostic for whether a model needs a spatial specification.
Local Moran runs n tests. Correct for multiple comparisons, or expect roughly 5% false positives.
Geary's C is an alternative that is more sensitive to local differences; the two often disagree at the margins, which is informative.

Internal links

How to calculate Moran's I in Python — the practical implementation
Spatial weights warn about islands or Moran's I returns nan — the commonest failure
How to find hotspots with Getis-Ord Gi* in Python — locating what the global statistic averages away
The modifiable areal unit problem explained — why the zones change the answer too
How to bin points into hexagons in Python — a neighbour structure with no diagonal ambiguity
Spatial join cardinality explained — the adjacency relationships underneath
Kernel density explained — the continuous alternative to zone-based statistics
How to aggregate spatial data by region in GeoPandas — producing the zone values in the first place

FAQ

What does Moran's I actually measure?

The correlation between each unit's value and the weighted average of its neighbours' values. Near +1 means similar values cluster; near −1 means they alternate.

Why is the expected value not zero?

Because I is computed from deviations about the sample mean. E[I] = −1/(n−1), which is negligible at n=400 and substantial at n=20. Always compare against moran.EI.

Which spatial weights should I use?

Choose from the process — contiguity for administrative units where boundaries matter, k-nearest for uneven unit sizes, a distance band when influence has a physical range. Test two or three and report the spread.

Why did my Moran's I change completely with different weights?

Because the pattern operates at the scale of the neighbourhood definition. A checkerboard gives I = −0.998 with rook weights and −0.039 with queen. If your result is that sensitive, the weights are the finding.

Should I use the analytical or permutation p-value?

The permutation one (p_sim). It makes no normality assumption. With 999 permutations the minimum reportable p is 0.001.

My Moran's I is significant. Now what?

If it was on raw values, map the local statistics to find where. If it was on regression residuals, your model is missing spatial structure and needs a spatial lag or error specification.

Can Moran's I be used on point data?

Not directly — it needs values attached to units with a neighbour structure. Aggregate the points to zones or a hex grid first, accepting the aggregation caveats.

Spatial Autocorrelation Explained #

Problem statement #

Quick answer #

Step-by-step solution #

1. Understand what the statistic actually compares #

2. Recognise that the weights are half the analysis #

3. Row-standardise, and know what it does #

4. Use the permutation p-value #

5. Interpret a non-significant result correctly #

Code examples #

Example 1 — global Moran's I with the choices made explicit #

Example 2 — the Moran scatterplot, which is where the interpretation lives #

Example 3 — from global to local #

Explanation #

Why E[I] is negative #

Why the weights are not a technicality #

Why autocorrelation invalidates ordinary statistics #

Why a global statistic can hide everything interesting #

Edge cases or notes #

Internal links #

FAQ #

What does Moran's I actually measure? #

Why is the expected value not zero? #

Which spatial weights should I use? #

Why did my Moran's I change completely with different weights? #

Should I use the analytical or permutation p-value? #

My Moran's I is significant. Now what? #

Can Moran's I be used on point data? #