Spatial Autocorrelation: is the map clustered, and where — Moran's I, global and local

The data in this project is synthetic — a regional indicator field generated locally by data/make_data.py, with a planted clustering structure. That is deliberate: grading needs the known local Moran’s I for every region, which real data doesn’t come labelled with — and the planted hot/cold clusters and isolated spikes make the technique’s behaviour exact. No signup, no GPU.

Before you start

New to spatial statistics? Start here.

Spatial autocorrelation is the tendency of nearby places to have similar values — median rent doesn’t jump at a county line; it drifts. Moran’s I is the standard number for it: run it over a map of some indicator and it comes back roughly between −1 and +1.

I > 0 — clustered: high values sit next to high, low next to low (rent, disease rates, unemployment usually look like this).
I ≈ 0 — no spatial pattern; the values are scattered at random across the map.
I < 0 — dispersed / checkerboard: high values tend to sit next to low.

Computing it needs one new object: a spatial weights matrix W, which just encodes who is a neighbour of whom. Row-standardise it (each region’s neighbour weights sum to 1) and the “spatial lag” W @ z becomes simply each region’s neighbours’ average.

There are two questions, and they are the two parts of this project:

Global (Part A): is there clustering anywhere on this map? One number (global Moran’s I), with a permutation test for significance.
Local (Part B): where are the clusters, region by region? The per-region decomposition, local Moran’s I (a.k.a. LISA). Crucially, a lone extreme value with ordinary neighbours is not a cluster — local Moran’s I catches that; a plain z-score does not.

New to Python too? See Start Here — it starts from zero, and shows how to run this in your browser with nothing to install.

1. The Brief

You have an indicator measured for every region on a map. Part A: build the spatial weights matrix, compute global Moran’s I, and use a permutation test to decide whether the clustering you see is real or could be an artefact of chance. Part B: go from the one global number to a per-region local Moran’s I, so you can point at the specific regions that form hot and cold clusters — and correctly not flag an isolated spike (a high value whose neighbours are unremarkable), which a naive value-threshold would wrongly call a hotspot.

Moran’s I is the workhorse of spatial epidemiology, regional economics, ecology, and crime analysis — anywhere “are these values spatially clustered, and where?” matters. The local statistic is what turns “yes, clustered” into an actual map of hot and cold spots.

2. Difficulty & time

Difficulty 6 / 10. Several moving parts that must fit together: construct a weights matrix from adjacency, standardise correctly, assemble the Moran’s I formula, run a permutation test, and then make the conceptual jump from the global sum to the per-region local statistic. The debug surface is real — a wrong standardisation or an off-by-one in the neighbour lookup quietly changes the answer. Calibrated against the ledger (G10): a peer of attrition-leakage (6), fraud-rings (6), and warehouse-sim (6). Above spatial-clusters (5) and hidden-communities (5), which are a single estimator call plus tuning; below the 7–8 band (multiple systems or a non-obvious algorithmic core).

Time is separate from difficulty: Part A is 60–90 minutes; Part B is 2.5–3 hours. The scripts run in under a second.

3. What you’ll be able to do after

Build a spatial weights matrix (queen contiguity) and row-standardise it, and explain what the spatial lag means.
Compute global Moran’s I and test its significance with a permutation test, without assuming any distribution.
Decompose it into local Moran’s I per region, and read the Moran scatterplot quadrants (HH, LL, HL, LH).
Tell a cluster from a spike — explain why an isolated extreme value has a near-zero local Moran’s I even though its z-score is large.

The finished result

By the end of Part A you’ll have answered “is this map clustered, and is that real?” — one global number, checked against chance:

global Moran's I = +0.175      E[I] under no autocorrelation = -0.005
permutation test:  permuted I mean -0.005, 95th pct +0.062  ->  observed +0.175 beats every shuffle
pseudo p-value = 0.001   (significant)
Moran scatterplot:  slope = global I = +0.175;  quadrants HH 58, LL 53, HL 42, LH 43

A real, significant cluster — high next to high, low next to low, and not a fluke; in Part B that one global number becomes a per-region local Moran’s I that points at the specific hot and cold clusters.

4. Prereqs & time box

You can write a function, use numpy arrays (including a matrix–vector product), and read a CSV with pandas. No spatial-statistics background needed — the orientation above and Part A build Moran’s I from its parts. You’ll see one matrix multiply (W @ z); the meaning is explained where it’s used.

Where to write and run code

Two ways to run this project — pick either:

In your browser (no install). If the project page shows a ▶ Run button, click it to run Part A right here — nothing to set up. Part A uses pandas, numpy, and scipy, which run in the browser.
On your computer. The one-time setup — installing Python, downloading the code from the course site (no GitHub, no git), opening a terminal, creating a virtual environment, installing the requirements, and running Part A step by step — is covered once in Start Here. Come back here once python -c "import scipy" runs without error. Part B (writing your own solution) needs this local setup.

Time box

Part A — Guided Build: 60–90 minutes. Hard cap.
Part B — Your Turn: 2.5–4 hours. Hard cap.

5. Setup & data

The data is synthetic and generated locally — see the banner above and SOURCES.md for full provenance. One command creates both parts’ data:

python data/make_data.py

That writes two CSVs into data/ (git-ignored; regenerate anytime, byte-for-byte identical because the seeds are fixed):

Part A — regions_a.csv (196 rows: region_id, row, col, value) — a 14×14 lattice of regions.
Part B — regions_b.csv (196 rows, same columns) — a different field on the same lattice.

Each field has a hot cluster and a cold cluster of moderately elevated/depressed regions, plus a handful of isolated spikes (an extreme value with ordinary neighbours). The value is a generic indicator; the numbers describe no real place.

6. Part A — Guided Build

You’ll build the weights matrix, compute global Moran’s I, test it, and read the Moran scatterplot. Run the whole thing, or (better) one # %% cell at a time — see Start Here:

python part_a/moran.py

Prefer the browser? If this page has a ▶ Run button, click it instead of the command above — Part A runs in your browser with nothing to install. The STEP output and every checkpoint below are identical either way.

The full, commented script is in part_a/moran.py. The spine:

Step 1 — Build the spatial weights matrix W.

Checkpoint: queen contiguity (the 8 surrounding cells); interior regions have 8 neighbours, corners 3. Row-standardised so each row sums to 1 — the spatial lag becomes the neighbours’ average.

Step 2 — Global Moran’s I.

Checkpoint: I = +0.175, well above E[I] = −1/(n−1) = −0.005. High regions sit next to high, low next to low.

Step 3 — Is it real? A permutation test.

Checkpoint: shuffling the values across regions gives permuted I centred on −0.005 (95th percentile ≈ +0.06); the observed +0.175 beats every shuffle, so pseudo p ≈ 0.001 — significant, not a fluke.

Step 4 — The Moran scatterplot (bridge to the local statistic).

Checkpoint: plotting each region’s value against its neighbours’ average gives a cloud whose slope equals global I (+0.175). The four quadrants label each region HH / LL / HL / LH. Writes part_a/moran_a.csv.

Why this and not that

Why row-standardise W. It turns the spatial lag into a plain neighbour average, so I is a clean correlation between a region’s value and its surroundings — comparable across regions with different neighbour counts (interior vs corner).
Why a permutation test. Moran’s I has a known theoretical variance, but the assumptions are fussy; shuffling the map and recomputing makes no distributional assumption and directly answers “is this more clustered than random?”
Why the scatterplot matters. The global I is the slope of that cloud — one number summarising every region. The quadrants are already the local story; Part B makes it a per-region statistic.

Functions you’ll meet (and where to read more)

Function	What it does here	Docs
`W @ z` (matrix–vector)	The spatial lag: neighbours’ average	np.matmul
`arr.sum(axis=1, keepdims=True)`	Row-standardise the weights	ndarray.sum
`rng.permutation`	Shuffle the map for the permutation test	Generator.permutation
`scipy.stats.spearmanr`	(Part B) rank-correlation for scoring	spearmanr

If something looks off

The permutation p-value is exactly 0.001, and it’s the same on every run — is that a bug? No, on both counts. The shuffle in Step 3 is seeded (np.random.default_rng(0)), so you get the same p every run; an unseeded permutation test is the one whose p-value would drift slightly run to run. And 0.001 is simply the smallest value 999 shuffles can report: the pseudo p-value is (1 + how many shuffles matched or beat your observed I) / (1 + 999), so when none of them did, it bottoms out at 1/1000. That’s a strong result, not a rounding glitch. (Relatedly, a Moran’s I near 0 is a real answer — “no spatial pattern” — not an error.)
A checkpoint number doesn’t match. The data is seeded, so the numbers are exact. Re-run python data/make_data.py (you may have skipped or edited it), then re-run Part A.
Setup problems (Python not found, pip, FileNotFoundError, re-activating the venv): see the troubleshooting list in Start Here.

7. Part B — Your Turn

Part B is the optional “your turn” half — a guided fill-in-the-blank, not a blank page. You get a working starter (part_b/starter.py) that already handles the easy version; your job is to finish the # TODO it marks out. More open than Part A’s step-by-step, but you are never staring at an empty file. If you finished Part A and want to prove the skill, this is where it happens. New to Python? It’s completely fine to stop after Part A and come back later.

Same core skill — Moran’s I — but local now, on a different field (regions_b.csv). Compute the per-region local Moran’s I: for each region, its standardised value times its neighbours’ spatial lag, I_i = z_i · (W z)_i. Global I (Part A) was one number for the whole map; local I is a vector, one value per region, that says how strongly each region participates in a cluster.

The transfer: Part A summed everything into a single statistic; Part B is the per-region decomposition — a genuinely different computation (element-wise, not a scalar sum). And it exposes the spike trap: an isolated high value has a large z-score but ordinary neighbours, so its spatial lag ≈ 0 and its local Moran’s I ≈ 0 — it is not a cluster. The starter reports the z-score alone (no neighbours) and scores 0.05, so it fails.

How you’re scored (reproduce the yardstick)

The grader computes the reference local Moran’s I for every region and scores the Spearman rank-correlation between your per-region values and the reference. Spearman is rank-based, so any valid normalisation of the statistic (standardised or not, ÷n or ÷(n−1)) matches — what has to be right is that you multiplied each region’s value by its neighbours’ spatial lag. The reference reaches 1.00; the floor is 0.90. The z-score-only baseline scores ~0.05.

Acceptance criteria (checkable)

Your solution writes into part_b/out/:

local_i.csv — columns region_id, local_i (every region, one row).
report.json — with keys regions, achieved_metric (you can leave achieved_metric null; the grader scores you).

Then:

python part_b/starter.py
pytest part_b/test_solution.py -q

If pytest reports everything as skipped, it just means the output files don’t exist yet — run python part_b/starter.py first (once your solution fills in the TODOs), then re-run.

You’re done when the tests pass (Spearman ≥ 0.90; the reference reaches 1.00).

Constraints — what must be true of your solution, not how to get there

Build the weights matrix from the row, col adjacency (queen contiguity) and row-standardise it.
local_i must use the neighbours’ spatial lag, not the value alone — an isolated spike should come out near zero, not high.
Every region appears once in local_i.csv.
Standardise the values once, globally (the same mean and sd for all regions).

Hints

Hint 1 — what "local" means

Global Moran’s I is (up to a constant) the average of the per-region products z_i · (W z)_i. The local statistic is just those products before you average them — one number per region. Reuse the z and the spatial lag W @ z from Part A; don’t sum them.

Hint 2 — why the z-score alone fails

local_i = z ignores the neighbours. A spike (big z, random neighbours) then looks like a hotspot, but its true local Moran’s I is ~0 because the spatial lag of random neighbours is ~0. Multiplying by (W z)_i is exactly what separates a cluster from a spike.

Hint 3 — the deliverable map (self-check, not graded)

Once you have local_i, classify each region by the Moran scatterplot quadrant (sign of z_i and sign of (W z)_i): HH and LL are cluster cells, HL/LH are spatial outliers. Colour the lattice by quadrant and you have a hot/cold-spot map — the thing a global number can’t give you.

8. Self-check

You don’t need the answer key. You’re done when all of these hold:

pytest passes — Spearman ≥ 0.90 against the reference local Moran’s I.
The regions with the highest local I are the ones inside the planted hot and cold blocks; the isolated spikes have local I near zero (not high). If a spike scores high, you used the value, not the lag.
Re-running from python data/make_data.py reproduces the same numbers — the data is seeded.

And the tell-tale: global Moran’s I equals (n/ΣW) × the mean of your z_i·(Wz)_i — a good arithmetic check that your local values are consistent with the global one.

9. Stretch

Different weights. Swap queen contiguity for k-nearest-neighbour weights (say k=8) and see how the local map shifts. Weights are a modelling choice, and the answer depends on it.
Local significance. Add a conditional-permutation test for each region’s local I (hold that region fixed, shuffle the rest) to get a per-region p-value, and map only the significant hot/cold spots — the full LISA workflow.
(Genuinely hard) A real indicator. Take a truly open, licence-clear regional dataset (many national statistics offices publish small-area indicators under an open licence — check it), build the contiguity from the boundary file, and map its hot spots. Note that without ground truth you judge the result by significance and stability, not a score.

10. Ship it

Put this in a portfolio repo and it counts. In that repo’s README:

State the finding in one sentence: “the indicator is significantly spatially clustered (global Moran’s I = 0.18, permutation p ≈ 0.001); local Moran’s I locates the hot and cold clusters and correctly leaves isolated spikes unflagged.”
Show the Moran scatterplot (slope = global I) and the local hot/cold-spot map side by side.
Include the weights-matrix construction and one paragraph on why an isolated spike has near-zero local I. That paragraph is what signals you understand spatial autocorrelation, not just a formula.

Keep make_data.py and the tests in the repo so anyone can reproduce your numbers from zero.

11. Sources

Rendered from SOURCES.md.

Field	Value
Dataset	Synthetic regional indicator on a 14×14 lattice (Part A 196, Part B 196): a hot cluster, a cold cluster, and isolated spikes
Type	Synthetic
Generated by	`data/make_data.py` (this repo), seeds `20260752` (A) and `20260751` (B)
Source URL	Not applicable — generated, not fetched
License	Not applicable — synthetic, no third-party rights
Access / generation date	2026-07-18
Redistribution	Permitted — synthetic; output regenerated from the seeded script, never committed
robots.txt checked	Not applicable — no scraping
ToS URL	Not applicable — no external source

Synthetic data teaches the mechanic; it does not prove a finding about the real world. The clustering was generated on purpose — the Moran’s I techniques are the point, not a claim about any real region.

Next up

Finished this one? Continue the Geospatial track:

Facility Location: which depots to open, and why nearest-assignment breaks under capacity → · Browse all courses