Banner

2.6 Sensitivity, Pitfalls and Exercises¶

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

Dataset: diabetes_northeast.shp (217 counties)

Focus: How much of Chapters 2.1–2.5's conclusions survive changes in weights, permutations, and correction choices — plus five worked exercises

Library: spatialstats.weights, spatialstats.esda, spatialstats.cluster (PySpatialStats)

Table of Contents¶

  1. Overview

  2. Python Setup

  3. Sensitivity: Weights Definition

  4. Sensitivity: Permutations, Seed, and Transform

    • 4.1 Random seed

    • 4.2 Row-standardised vs. binary weights

  5. Pitfalls

  6. Exercises

    • 6.1 Spatial lag by hand and with W.lag

    • 6.2 Rank six weight definitions by Moran's I

    • 6.3 Compare LISA clusters for Obesity and Dia_pct

    • 6.4 Permutations needed before FDR flags anything

    • 6.5 Reproduce global Moran's I with libpysal/esda

  7. Summary and Conclusions

  8. Resources

1. Overview¶

Every chapter so far reported a specific number — I = 0.240, 29 significant LISA counties, 156 BB joins — under one particular set of choices (Queen weights, 999 permutations, a particular seed, fdr=False). This closing chapter asks the question a careful reader should already be asking: how much of that is really about the data, and how much is an artefact of the choices?

Step What we do
Weights sensitivity Repeat LISA under three weight definitions and measure agreement with flag_stability and agreement_matrix
Permutation/seed sensitivity Quantify how much significance counts move with the RNG seed and permutation budget
Transform sensitivity Row-standardised vs. binary weights: does Moran's I change?
Pitfalls A consolidated list of ways Part 2's numbers can mislead, with pointers to where each is fixed
Exercises Five solved exercises spanning Chapters 2.2–2.5

This is also where Part 2 formally hands off: Part 3 (spatial regression) asks what to do once dependence is confirmed; Part 4 (cluster analysis) gives the local statistics of Chapter 2.5 their full, calibrated treatment; Part 5 (disease mapping) addresses the small-population instability flagged in Section 5 below.

2. Python Setup¶

In [1]:
import sys, 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, Rook, KNN, DistanceBand, Delaunay, Gabriel, W
from spatialstats.esda import Moran, Moran_Local
from spatialstats.cluster import local_moran_clusters, flag_stability, agreement_matrix

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")
print(f"spatialstats : version {spatialstats.__version__}")
print(f"counties     : {counties.shape[0]}")

spatialstats : version 0.2.0
counties     : 217

3. Sensitivity: Weights Definition¶

Chapter 2.2, Section 12 already showed Moran's I moving from 0.116 to 0.273 across six weight definitions. Does the same hold for local significance — do the same counties get flagged regardless of which W is used?

In [2]:
weight_defs = {
    "Queen": Queen(counties),
    "KNN(6)": KNN(counties, k=6),
    "DistanceBand(100km)": DistanceBand(counties, threshold=100_000, binary=True),
}
masks = {}
for name, w in weight_defs.items():
    lmc = local_moran_clusters(counties["Dia_pct"], w, permutations=999, seed=1, correction="none")
    masks[name] = lmc.labels != "Not significant"
    print(f"{name:22s}: {masks[name].sum()} significant counties")

Queen                 : 30 significant counties

KNN(6)                : 39 significant counties

DistanceBand(100km)   : 59 significant counties
In [3]:
agreement_matrix(masks).round(3)
Out[3]:
Queen KNN(6) DistanceBand(100km)
Queen 1.000 0.50 0.328
KNN(6) 0.500 1.00 0.380
DistanceBand(100km) 0.328 0.38 1.000
In [4]:
stability = flag_stability(list(masks.values()))
counties_stab = counties.assign(stability=stability)
summary = pd.Series(stability).value_counts().sort_index(ascending=False)
print("share of the 3 weight definitions each county is flagged under:")
print(summary)
print(f"\nrobust (flagged under all 3): {int((stability == 1.0).sum())} counties")
print(f"flagged under exactly 1 of 3 (i.e. a single-definition artefact): {int(np.sum((stability > 0) & (stability < 1.0)))} counties")

share of the 3 weight definitions each county is flagged under:
1.000000     17
0.666667     21
0.333333     35
0.000000    144
Name: count, dtype: int64

robust (flagged under all 3): 17 counties
flagged under exactly 1 of 3 (i.e. a single-definition artefact): 56 counties

Agreement (Jaccard) between weight definitions ranges from 0.33 to 0.50 — moderate at best. Only 17 of the up to 59 counties ever flagged under any definition are flagged under all three; the majority (56 counties) are flagged under exactly one. The practical conclusion: a LISA cluster map answers "where does this particular neighbourhood definition find local structure", not an absolute, definition-free fact about the county. Report the weights definition alongside every cluster map, and treat single-definition-only flags as hypotheses, not findings.

4. Sensitivity: Permutations, Seed, and Transform¶

4.1 Random seed¶

Permutation p-values are themselves random variables — a different seed gives a slightly different answer at any fixed permutation count.

In [5]:
w_queen = Queen(counties)
seed_counts = []
for seed in range(1, 11):
    lmc = local_moran_clusters(counties["Dia_pct"], w_queen, permutations=999, seed=seed, correction="none")
    seed_counts.append(int((lmc.labels != "Not significant").sum()))
print(f"significant-county counts across 10 seeds (999 perms, uncorrected): {seed_counts}")
print(f"range: {min(seed_counts)}-{max(seed_counts)}  (mean {np.mean(seed_counts):.1f}, std {np.std(seed_counts):.2f})")

significant-county counts across 10 seeds (999 perms, uncorrected): [30, 32, 30, 29, 30, 31, 30, 32, 30, 31]
range: 29-32  (mean 30.5, std 0.92)

4.2 Row-standardised vs. binary weights¶

Section 8 of Chapter 2.2 explained what the transforms mean; here is whether they change the headline number.

In [6]:
w_row = Queen(counties, transform="row")
w_bin = Queen(counties, transform="binary")
i_row = Moran(counties["Dia_pct"], w_row, permutations=999, seed=1)
i_bin = Moran(counties["Dia_pct"], w_bin, permutations=999, seed=1)
print(f"Moran's I, row-standardised : {i_row.I:.4f}  (p_sim={i_row['p_sim']:.3f})")
print(f"Moran's I, binary           : {i_bin.I:.4f}  (p_sim={i_bin['p_sim']:.3f})")

Moran's I, row-standardised : 0.2398  (p_sim=0.001)
Moran's I, binary           : 0.2119  (p_sim=0.001)

The two do not agree (0.240 vs. 0.212) — row-standardisation and binary weighting are not just cosmetic rescalings of the same underlying quantity, because row-standardisation implicitly down-weights counties with many neighbours relative to counties with few. Both are legitimate choices; neither is "the" Moran's I for this dataset, which is why Chapter 2.2 flagged the transform as a modelling decision, not a technical default to ignore.

5. Pitfalls¶

A consolidated list, drawing on every chapter in this Part, of ways a Part-2 result can mislead if read without its context:

Pitfall Where it shows up Where it is addressed
Small-population instability. Rates from sparsely populated counties are noisier, and that noise itself can look like spatial "clustering" of extremes purely from sampling variance, not a real regional effect. Not directly visible in Dia_pct (a well-measured CDC rate), but present in any small-area count-based rate Part 5 (empirical Bayes and BYM smoothing)
Assumed global stationarity. Moran's I, Geary's C, and LISA all implicitly assume the strength of spatial dependence is constant across the study region — a single number (or a single per-county number computed the same way everywhere) cannot capture a relationship that itself varies geographically. Every statistic in this Part Part 6 (geographically weighted models)
Edge effects, MAUP, and projection. Counties on the study region's boundary have artificially truncated neighbourhoods (fewer real-world neighbours than an interior county); the areal units themselves are administrative, not natural (Chapter 2.1); and every distance-based weight in Chapter 2.2 depends on using a genuinely equal-area/equal-distance CRS. Chapters 2.1, 2.2 Part 0 §0.4–0.5 (CRS), Chapter 2.1 §4 (MAUP)
Clusters are not causes. A significant HH cluster says where similar high values sit together; it says nothing about why — shared drivers (Section 3 of Chapter 2.1), coincidental co-location, or genuine spatial contagion are all consistent with the same LISA map. Chapter 2.5 Part 3 (spatial regression is the tool for explaining clustering, not just detecting it)
FDR under-powered at moderate n. Chapter 2.5 showed FDR correction can report zero significant local clusters even with an unambiguous global signal, unless permutations are increased by an order of magnitude. Chapter 2.5 Chapter 2.5 §4.3; Part 4's inference='analytic' alternative
Weights-definition dependence. Every number in this Part changes, sometimes substantially, with the weights definition (Section 3 above; Chapter 2.2 §12). All chapters This chapter, Section 3

6. Exercises¶

6.1 Spatial lag by hand and with W.lag¶

Compute the spatial lag on Dia_pct under Queen weights by hand (a direct loop) and with W.lag, and confirm they agree. (Chapter 2.2 did this on the toy lattice; here it is redone on the real dataset used throughout this Part, to confirm the identity is not an artefact of the toy example.)

In [7]:
w = Queen(counties)
z = counties["Dia_pct"].to_numpy() - counties["Dia_pct"].mean()

lag_hand = np.zeros(len(z))
for idx in range(w.n):
    id_i = w.id_order[idx]
    nbrs = w.neighbors[id_i]
    if nbrs:
        lag_hand[idx] = np.mean([z[w.id_to_index[j]] for j in nbrs])

lag_builtin = w.lag(z)
print(f"max |hand - W.lag()| over 217 counties: {np.max(np.abs(lag_hand - lag_builtin)):.2e}")

max |hand - W.lag()| over 217 counties: 2.22e-16

Answer: agreement to machine precision (as in Chapter 2.2), confirming W.lag(z) is exactly the row-standardised neighbourhood mean, on real as well as toy data.

6.2 Rank six weight definitions by Moran's I¶

Using the six weight definitions from Chapter 2.2, Section 12, rank them by Moran's I on Dia_pct and explain the pattern.

In [8]:
six_defs = {
    "Queen": Queen(counties), "Rook": Rook(counties), "KNN(6)": KNN(counties, k=6),
    "DistanceBand(100km)": DistanceBand(counties, threshold=100_000, binary=True),
    "Delaunay": Delaunay(counties), "Gabriel": Gabriel(counties),
}
rows = [{"weights": name, "mean_neighbours": w.cardinalities().mean(),
         "Moran_I": Moran(counties["Dia_pct"], w, permutations=999, seed=1).I}
        for name, w in six_defs.items()]
ranked = pd.DataFrame(rows).sort_values("Moran_I", ascending=False).reset_index(drop=True)
ranked
Out[8]:
weights mean_neighbours Moran_I
0 Gabriel 4.497696 0.272861
1 Delaunay 5.870968 0.248073
2 Queen 5.188940 0.239826
3 Rook 5.152074 0.239754
4 KNN(6) 6.000000 0.231631
5 DistanceBand(100km) 14.608295 0.116121
In [9]:
corr = np.corrcoef(ranked["mean_neighbours"], ranked["Moran_I"])[0, 1]
print(f"correlation between mean neighbour count and Moran's I across the six definitions: {corr:.3f}")

correlation between mean neighbour count and Moran's I across the six definitions: -0.984

Answer: the ranking is (highest to lowest I) Gabriel > Delaunay > Rook ≈ Queen > KNN(6) > DistanceBand. The near-perfect negative correlation (−0.98) with mean neighbour count is the explanation: sparser, more-local weight definitions (Gabriel, with only 4.5 neighbours on average) give Moran's I less averaging-out room, so genuinely local similarity shows through more strongly; denser definitions (the 100 km distance band, with 14.6 neighbours) average each county against a wider, less locally-representative neighbourhood, diluting the same signal toward zero.

6.3 Compare LISA clusters for Obesity and Dia_pct¶

In [10]:
lm_dia = Moran_Local(counties["Dia_pct"], w, permutations=999, seed=1, fdr=False)
lm_obesity = Moran_Local(counties["Obesity"], w, permutations=999, seed=1, fdr=False)

print(f"Dia_pct  significant: {(lm_dia.labels != 'NS').sum()}")
print(f"Obesity  significant: {(lm_obesity.labels != 'NS').sum()}")

am = agreement_matrix({"Dia_pct": lm_dia.labels != "NS", "Obesity": lm_obesity.labels != "NS"})
print(f"\nJaccard agreement (significant/not, any label): {am.loc['Dia_pct', 'Obesity']:.3f}")

both_hh = int(((lm_dia.labels == "HH") & (lm_obesity.labels == "HH")).sum())
both_ll = int(((lm_dia.labels == "LL") & (lm_obesity.labels == "LL")).sum())
print(f"counties in the SAME cluster type (HH) for both variables: {both_hh}")
print(f"counties in the SAME cluster type (LL) for both variables: {both_ll}")

Dia_pct  significant: 29
Obesity  significant: 32

Jaccard agreement (significant/not, any label): 0.271
counties in the SAME cluster type (HH) for both variables: 8
counties in the SAME cluster type (LL) for both variables: 3

Answer: the two variables' LISA maps agree on which counties are flagged at all only moderately (Jaccard 0.271) — obesity and diabetes prevalence, while strongly correlated within the same county (Chapter 2.4: same-region $r=0.68$), have distinct local geographies. Where they do agree on cluster type, it is concentrated in the HH cluster (8 counties jointly high-high) more than LL (3 jointly low-low) — consistent with a shared driver in the highest-prevalence region but more divergent local structure elsewhere.

6.4 Permutations needed before FDR flags anything¶

This was worked out directly in Chapter 2.5, Section 4.3: with local_moran_clusters at correction='fdr', zero counties survive at 999, 2,999, 4,999, 6,999, or 8,999 permutations; the first county to survive appears at 9,999 permutations (one HH cluster). The library's own internal warning (triggered whenever permutations is too small for n at the requested alpha) estimates a slightly lower threshold of ≈8,680 permutations — consistent with 9,999 being the first round number past that estimated minimum. The reason, restated: FDR at $n=217$ needs the smallest attainable p-value, $2/(m{+}1)$, to be small enough that at least a handful of counties' true p-values can rank below the Benjamini–Hochberg threshold; increasing $m$ shrinks that floor directly, while the counties' actual p-values do not change with $m$ once $m$ is already reasonably large — so the fix is squarely about giving the correction enough resolution, not about the underlying statistic being weak.

6.5 Reproduce global Moran's I with libpysal/esda¶

In [11]:
try:
    import esda as pysal_esda
    lw = w.to_libpysal()
    mine = Moran(counties["Dia_pct"], w, permutations=999, seed=1)
    np.random.seed(1)
    theirs = pysal_esda.Moran(counties["Dia_pct"].to_numpy(), lw, permutations=999)
    print(f"{'':16s}{'this library':>16s}{'libpysal/esda':>16s}")
    print(f"{'I':16s}{mine.I:16.6f}{theirs.I:16.6f}")
    print(f"{'z (norm)':16s}{mine.z_score:16.4f}{theirs.z_norm:16.4f}")
    print(f"{'p_norm':16s}{mine['p_norm']:16.3e}{theirs.p_norm:16.3e}")
    print(f"{'p_sim':16s}{mine['p_sim']:16.4f}{theirs.p_sim:16.4f}")
except ImportError:
    print("libpysal/esda not installed; skipping cross-check (pip install libpysal esda)")

('WARNING: ', 33, ' is an island (no neighbors)')
                    this library   libpysal/esda
I                       0.239826        0.239826
z (norm)                  5.6107          5.6092
p_norm                 2.015e-08       2.032e-08
p_sim                     0.0010          0.0010

/home/zia207/Github/PySpatialStats/spatialstats/weights/w.py:225: UserWarning: The weights matrix is not fully connected: 
 There are 2 disconnected components.
 There is 1 island with id: 33.
  return LibW(self.neighbors, self.weights, id_order=self.id_order)

Answer: $I$ agrees to five decimal places (0.239826 vs. 0.239826), and the z-score and p-values agree closely (small residual differences come from slightly different finite-population variance corrections, not a disagreement about the substantive conclusion). This is the strongest possible validation available for Moran: an independent, widely-used reference implementation reproduces it on real data.

7. Summary and Conclusions¶

Part 2 built the full pipeline from "does this areal map look clustered?" to a validated, sensitivity-checked answer, and this closing chapter quantified how much of that answer is robust.

What we did¶

  1. Measured agreement between LISA cluster flags under three weight definitions: only 17 of up to 59 ever-flagged counties are robust to all three; most are single-definition artefacts.
  2. Quantified seed-to-seed variability in significant-county counts (range across 10 seeds, roughly ±3 around a mean of 30) and confirmed row-standardised and binary weights give materially different Moran's I (0.240 vs. 0.212).
  3. Consolidated six pitfalls spanning this Part, each with a pointer to the Part that resolves it.
  4. Solved five exercises: hand-computed spatial lag on real data, ranked six weight definitions and explained the pattern via mean neighbour count, compared Obesity and Dia_pct LISA maps, pinned down the exact permutation count needed for FDR to flag anything, and cross-validated Moran's I against libpysal/esda.

Practical takeaways¶

  • No single number from this Part — Moran's I, a LISA label, a join count — should be reported without its weights definition, permutation count, and correction choice; all three move the answer.
  • Robustness across multiple reasonable weight definitions is stronger evidence of real spatial structure than significance under any one definition alone.
  • The apparatus built in this Part (weights, global and local autocorrelation, visualisation) is now validated against an external reference implementation and ready to be the diagnostic layer for Part 3's spatial regression models.

Next steps¶

Part 3 (Spatial Regression) picks up exactly where Part 2 leaves off: given confirmed spatial dependence in Dia_pct, it asks whether that dependence reflects a substantive spillover (spatial lag model) or unmodelled spatial structure in the residuals (spatial error model), and how to choose between them formally.

8. Resources¶

Textbooks¶

Reference Notes
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic. Comprehensive treatment tying together weights, autocorrelation, and the regression models of Part 3
Haining, R. (2003). Spatial Data Analysis: Theory and Practice. Cambridge University Press. Strong treatment of sensitivity and robustness in exploratory spatial analysis
Fotheringham, A. S., Brunsdon, C. & Charlton, M. (2000). Quantitative Geography: Perspectives on Spatial Data Analysis. Sage. Accessible discussion of MAUP, stationarity, and their practical consequences

Key papers¶

Paper Contribution
Openshaw, S. & Taylor, P. J. (1979). A million or so correlation coefficients: three experiments on the modifiable areal unit problem. In Statistical Applications in the Spatial Sciences. Foundational empirical demonstration of weights/zoning sensitivity
Griffith, D. A. (1996). Some guidelines for specifying the geographic weights matrix contained in spatial statistical models. Practical guidance directly relevant to Section 3's sensitivity results

In this library¶

Object Role
spatialstats.cluster.flag_stability, agreement_matrix Cross-setting robustness checks, used throughout Section 3
spatialstats.weights.W.to_libpysal Cross-validation against the reference PySAL ecosystem (Exercise 6.5)
Part 3 spatialstats.regression — the next step once dependence is confirmed
Part 4 Full hot-spot, scan-statistic, and regionalisation toolkit, extending Chapter 2.5's LISA
Part 5 Empirical Bayes and BYM smoothing, addressing the small-population pitfall of Section 5
Part 6 Geographically weighted models, addressing the global-stationarity pitfall of Section 5