4.8 Exercises

Part: 4 — Spatial Cluster Analysis (Chapter 4.8 of 8)

Dataset: diabetes_northeast.shp (217 counties)

Focus: Five solved exercises spanning Chapters 4.2, 4.4, and 4.5

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

---

Table of Contents

1. Overview

2. Python Setup

   • 2.1 Gi* vs. Local Moran agreement

   • 2.2 Scan-statistic window size sensitivity

   • 2.3 The permutation count where Gi* FDR starts flagging counties

   • 2.4 Agglomerative vs. SKATER: which is contiguous?

   • 2.5 A minimum-population floor for SKATER, and the island

3. Summary and Conclusions

4. Resources

---

1. Overview

Five solved exercises, each revisiting a specific method from Chapters 4.2–4.5 with a fresh, self-contained computation.

#	Question	Chapter revisited
1	Do Gi* and Local Moran flag the same counties for Dia_pct?	4.2
2	What changes in the scan statistic between 10% and 50% maximum window sizes?	4.4
3	At what permutation count does Gi* + FDR start flagging anything?	4.2
4	Which regionalisation method is actually contiguous?	4.5
5	What happens to the island under a minimum-population floor?	4.5

---

2. Python Setup

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
from spatialstats.cluster import (
    getis_ord_gi, local_moran_clusters, spatial_scan, skater, agglomerative,
    agreement_matrix, partition_agreement,
)

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__}")

spatialstats : version 0.2.0

2.1 Gi* vs. Local Moran agreement

Do Gi and Local Moran flag the same counties for Dia_pct? Quantify with agreement_matrix, for both the hot/HH and cold/LL sides.*

gi_star = getis_ord_gi(counties["Dia_pct"], w, star=True, inference="analytic", correction="fdr")
lisa = local_moran_clusters(counties["Dia_pct"], w, permutations=999, seed=1, correction="none")

masks = {
    "Gi* hot": gi_star.hot, "Gi* cold": gi_star.cold,
    "LISA HH": lisa.labels == "HH", "LISA LL": lisa.labels == "LL",
}
print({k: int(v.sum()) for k, v in masks.items()})
agreement_matrix(masks).round(3)

{'Gi* hot': 7, 'Gi* cold': 3, 'LISA HH': 13, 'LISA LL': 11}

	Gi* hot	Gi* cold	LISA HH	LISA LL
Gi* hot	1.000	0.000	0.538	0.000
Gi* cold	0.000	1.000	0.000	0.273
LISA HH	0.538	0.000	1.000	0.000
LISA LL	0.000	0.273	0.000	1.000

Answer: Gi* hot vs. LISA HH: Jaccard 0.54; Gi* cold vs. LISA LL: Jaccard 0.27 (both reproduced from Chapter 4.2). Neither pair is identical, and the cold/LL pair agrees noticeably less than the hot/HH pair — LISA's uncorrected significance screen (30 counties total) is simply more permissive than Gi's FDR-corrected one (10 counties total), so LISA's LL set is the larger of the two and Gi's cold set is close to a subset of it rather than a near-identical match.

2.2 Scan-statistic window size sensitivity

Re-run the scan with maximum window sizes of 10% and 50% of the population; what changes?

coords = np.c_[counties.geometry.centroid.x, counties.geometry.centroid.y]
cases = counties["Dia_count"].to_numpy()
pop = counties["POP_Total"].to_numpy()

for frac in (0.10, 0.50):
    res = spatial_scan(coords, cases, population=pop, max_pop_fraction=frac, permutations=999, seed=1, n_clusters=3)
    print(f"\nmax_pop_fraction = {frac}:")
    print(res.clusters[["n_locations", "rr", "p_value", "significant"]].round(4).to_string(index=False))

max_pop_fraction = 0.1:
 n_locations     rr  p_value  significant
           3 1.1668    0.001         True
           9 1.1623    0.001         True
          47 1.1270    0.001         True

max_pop_fraction = 0.5:
 n_locations     rr  p_value  significant
           3 1.1668    0.001         True
          95 1.1120    0.001         True
           9 1.1623    0.001         True

Answer: the most-likely (3-county) and one secondary (9-county) cluster are identical under both settings — they are small enough to fit within a 10% population cap regardless. What changes is the second secondary cluster: capped at 10%, it tops out at 47 counties with a relative risk of 1.127; allowed up to 50%, it grows to 95 counties with a slightly lower relative risk of 1.112. A larger allowed window can absorb more territory, but the extra counties it picks up are, on average, less elevated than the core — diluting the relative risk even as the window's total case excess grows. Capping the window size trades a smaller, more/tightly elevated cluster for a larger, more diffuse one; which is more useful depends on whether the application wants "where is risk most concentrated" (small cap) or "how far does the elevated-risk region extend" (larger cap).

2.3 The permutation count where Gi* FDR starts flagging counties

Show why 999 permutations plus FDR flags nothing, and find the number of permutations at which it starts to.

for m in (999, 2999, 4999, 6999, 8999, 9999):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        res = getis_ord_gi(counties["Dia_pct"], w, star=True, inference="permutation",
                           permutations=m, correction="fdr", seed=1)
    n_flagged = int((res.hot | res.cold).sum())
    print(f"m={m:5d}  (p-floor={2/(m+1):.5f})  ->  {n_flagged} flagged")

m=  999  (p-floor=0.00200)  ->  0 flagged
m= 2999  (p-floor=0.00067)  ->  0 flagged

m= 4999  (p-floor=0.00040)  ->  0 flagged

m= 6999  (p-floor=0.00029)  ->  4 flagged

m= 8999  (p-floor=0.00022)  ->  3 flagged

m= 9999  (p-floor=0.00020)  ->  1 flagged

Answer: exactly Chapter 4.2's finding — at $m=999$ the smallest attainable p-value ($2/1000=0.002$) is not small enough for any county to survive FDR correction across 217 simultaneous tests. Zooming in around the transition:

for m in (4200, 4400, 4600, 4800, 4900, 4950, 4999, 5000, 5500, 6000, 6500):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        res = getis_ord_gi(counties["Dia_pct"], w, star=True, inference="permutation",
                           permutations=m, correction="fdr", seed=1)
    n_flagged = int((res.hot | res.cold).sum())
    print(f"m={m:5d}  ->  {n_flagged} flagged")

m= 4200  ->  0 flagged

m= 4400  ->  4 flagged

m= 4600  ->  3 flagged

m= 4800  ->  0 flagged

m= 4900  ->  3 flagged

m= 4950  ->  4 flagged

m= 4999  ->  0 flagged

m= 5000  ->  1 flagged

m= 5500  ->  1 flagged

m= 6000  ->  1 flagged

m= 6500  ->  4 flagged

The transition is not a single clean threshold — the flagged count fluctuates between 0 and 4 in the $m\approx4200$–$6500$ range before settling into consistently nonzero territory. This is genuine Monte Carlo noise: which specific counties' p-values happen to land just below or just above the FDR cutoff at a given permutation count depends on the exact random draws, not on a smooth deterministic function of $m$. The practical answer is therefore a range, roughly 4,000–7,000 permutations for this dataset and seed, rather than one exact number — consistent with (and a bit more precise than) the library's own built-in warning, which estimates a conservative minimum of about 8,680 permutations.

2.4 Agglomerative vs. SKATER: which is contiguous?

attributes = counties[["Dia_pct", "Obesity", "SVI"]]
regions_skater = skater(attributes, w, n_clusters=6)
regions_agg = agglomerative(attributes, w, n_clusters=6)

print(f"SKATER        : R²={regions_skater.r2:.4f}, contiguous={regions_skater.is_contiguous(w)}")
print(f"Agglomerative  : R²={regions_agg.r2:.4f}, contiguous={regions_agg.is_contiguous(w)}")

agreement = partition_agreement(regions_skater.labels, regions_agg.labels)
print(f"partition agreement: adjusted Rand={agreement['adjusted_rand']:.3f}, NMI={agreement['nmi']:.3f}")

SKATER        : R²=0.2723, contiguous=True
Agglomerative  : R²=0.3496, contiguous=False
partition agreement: adjusted Rand=0.672, NMI=0.670

Answer: SKATER is contiguous; agglomerative is not (Chapter 4.5, Section 4.1) — despite agglomerative's higher $R^2$ (it is not constrained by a tree structure, so it can find a strictly better attribute-space fit). SKATER's minimum-spanning-tree construction makes contiguity a structural guarantee; agglomerative's connectivity= constraint is a preference that scikit-learn's implementation can silently override when the neighbour graph has more than one connected component, as this data's isolated island ensures it does. When strict contiguity is a hard requirement (as it usually is for redistricting-style applications), SKATER is the safer choice; when the best possible attribute fit matters more than contiguity, agglomerative can outperform it — but its output must always be checked with is_contiguous, never assumed.

2.5 A minimum-population floor for SKATER, and the island

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    regions_floor = skater(attributes, w, n_clusters=6, min_size=5,
                           floor=counties["POP_Total"], floor_value=500_000)
    for wmsg in caught:
        print("WARNING:", str(wmsg.message))

print(f"\nregion sizes: {sorted(regions_floor.sizes, reverse=True)}")
region_pops = pd.Series(counties["POP_Total"].to_numpy()).groupby(regions_floor.labels).sum()
print(f"region populations: {sorted(region_pops.round(0).astype(int).tolist())}")

WARNING: 1 disconnected component(s) of the neighbour graph (e.g. islands; sizes [1]) cannot meet min_size/floor on their own. They cannot be joined to another region, so each forms its own region and violates the constraint. Link them to a neighbour in `w` if that is not acceptable.

region sizes: [np.int64(136), np.int64(46), np.int64(20), np.int64(8), np.int64(6), np.int64(1)]
region populations: [11268, 2660249, 3753068, 7367630, 18753114, 23579439]

Answer: the island (Nantucket) still ends up as its own one-county region, still below the population floor (its own population is far under 500,000), and the warning fires again, naming the same reason as before: a disconnected component of the neighbour graph can never be merged into another region by skater's tree-cutting procedure, no matter what min_size or floor_value is requested — the constraint can only be enforced on cuts within a connected component, and the island is not part of any tree that has other counties to cut into it. The floor and minimum size do take effect everywhere else (every other region comfortably exceeds both), which confirms the constraint mechanism itself works correctly — it is specifically inapplicable to a graph component with no other member to share a region with. The only way to bring the island into compliance is the same one noted in Chapter 4.5: link it to a neighbour in w explicitly before regionalising.

---

3. Summary and Conclusions

This chapter closed Part 4 with five self-contained, freshly computed exercises.

What we did

1. Quantified Gi*-vs-LISA agreement on both the hot/HH and cold/LL sides, finding the hot side agrees more (Jaccard 0.54) than the cold side (0.27), traced to LISA's more permissive uncorrected screening.
2. Compared scan-statistic results at 10% vs. 50% maximum window fractions: the tightest clusters are identical, but the largest allowed cluster grows from 47 to 95 counties with a correspondingly diluted relative risk.
3. Located the (noisy, not sharp) permutation-count range — roughly 4,000–7,000 — where Gi* + FDR first starts flagging counties for this dataset, and explained why the transition is not a single clean threshold.
4. Confirmed SKATER is contiguous by construction while constrained agglomerative clustering is not, despite the latter's higher attribute-space $R^2$.
5. Showed a minimum-population floor is correctly enforced everywhere except the one disconnected island component, which structurally cannot comply regardless of the constraint's specific value.

Part 4 in review

Across eight chapters, Part 4 built and validated the full spatial cluster-detection toolkit: statistical hot spots (Gi*/Gi, benchmarked against Part 2's LISA), a population-adjusted scan statistic (calibrated and power-checked with real Monte Carlo studies), density-based point clustering (with a concretely demonstrated duplicate-location hazard), and spatially constrained regionalisation (with a concretely demonstrated contiguity gap in a commonly used method). Chapter 4.6 showed these methods should be compared, not treated as interchangeable, and Chapter 4.7 previewed the two most direct next steps: Part 5 (population-adjusted smoothing for the small-area instability problem) and Part 9 (formal space-time extensions of this Part's purely spatial, single-snapshot methods).

---

4. Resources

Textbooks

Reference	Notes
Fischer, M. M. & Getis, A. (Eds.) (2010). Handbook of Applied Spatial Analysis. Springer.	Comprehensive reference for every method exercised in this chapter

Key papers

Paper	Contribution
Benjamini, Y. & Hochberg, Y. (1995). Controlling the false discovery rate. JRSS B, 57, 289–300.	The FDR correction whose resolution limit is explored in Exercise 2.3

In this library

Object	Role
spatialstats.cluster.*	Every function used across Part 4, exercised fresh in this closing chapter
Part 5	Population-adjusted smoothing — the direct continuation flagged in Chapter 4.7
Part 9	Space-time cluster detection — the other direct continuation flagged in Chapter 4.7