Banner

5.9 Shrinkage Assessment and Interpretation¶

Part: 5 — Disease Mapping and Bayesian Smoothing (Chapter 5.9 of 11)

Dataset: diabetes_northeast.shp (217 counties)

Focus: Which areas move the most under smoothing, why, and a numerical caveat in the library's own shrinkage ratio

Library: spatialstats.bayes (PySpatialStats)

Table of Contents¶

  1. Overview

  2. Python Setup

  3. Two Kinds of Shrinkage

  4. The Shrinkage Ratio, and a Numerical Caveat

  5. A More Reliable Diagnostic: Absolute Change

  6. Which Areas Move Most?

  7. Summary and Conclusions

  8. Resources

1. Overview¶

Chapter 5.8 validated that smoothing improves estimates on average. This chapter asks the complementary, area-by-area question: for any specific county, how much did smoothing change its estimate, and can that be explained? Along the way, a real numerical caveat in this library's own shrinkage() ratio surfaces — worth knowing before trusting it blindly.

Step What we do
Two kinds of shrinkage Toward the global mean (EB, Chapter 5.5) vs. toward the neighbourhood (BYM, Chapters 5.6–5.7)
The shrinkage ratio BYM.shrinkage() / shrinkage_summary — and where its ratio-based definition becomes numerically unstable
A more reliable diagnostic Plain absolute change, and why it tells the more trustworthy story here
Which areas move most Small population, extreme raw values, and isolation — checked directly, not assumed

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
from spatialstats.bayes import expected_counts, smr, BYM, shrinkage_summary

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, transform="binary")
O = counties["Dia_count"].to_numpy()
POP = counties["POP_Total"].to_numpy()
E = expected_counts(POP, observed=O)
raw = smr(O, E)
fit = BYM(O, E, w, model="bym", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)
print(f"spatialstats : version {spatialstats.__version__}")

spatialstats : version 0.2.0

3. Two Kinds of Shrinkage¶

Chapter 5.5's empirical Bayes methods pull every county toward a single reference point — the global mean (or, for local EB, the neighbourhood mean, computed once and fixed). BYM's random effects instead let every county's implicit reference point emerge jointly from the whole model: the spatial effect $u_i$ pulls toward whatever the current posterior estimate of the local neighbourhood level is, which itself depends on every other county's data simultaneously. The practical difference: EB's shrinkage target is a simple, interpretable summary statistic; BYM's is a moving target implicit in the fitted model, more flexible but less directly interpretable in isolation.

4. The Shrinkage Ratio, and a Numerical Caveat¶

BYM.shrinkage() reports, per area, the fraction of the distance between the raw SMR and the overall mean risk that smoothing removed: $s_i = (\text{raw}_i - \hat\theta_i) / (\text{raw}_i - \bar g)$, where $\bar g = \sum O / \sum E$. This is a natural definition — but look at what it actually returns on real data before trusting it as a ranking.

In [2]:
shrink_ratio = fit.shrinkage()
ranked = pd.DataFrame({
    "COUNTY": counties["COUNTY"], "E": E, "raw_SMR": raw, "shrinkage_ratio": shrink_ratio,
}).sort_values("shrinkage_ratio", ascending=False)
ranked.head(5)
Out[2]:
COUNTY E raw_SMR shrinkage_ratio
15 Lincoln County 2536.484255 0.998232 0.634582
209 Grand Isle County 521.687415 1.002516 0.430475
177 Montour County 1342.515294 0.990678 0.298089
107 Otsego County 4392.215860 0.992438 0.233117
33 Nantucket County 830.006136 0.807223 0.122870

The counties with the largest shrinkage ratio are not the smallest-$E$ or most-extreme-raw counties at all — they are counties whose raw SMR happens to sit very close to 1.0 (the overall mean). The reason is purely arithmetic: $s_i$'s denominator is $\text{raw}_i - \bar g$, which is tiny whenever a county's raw rate is already near average — so even a small absolute change in that county's estimate produces a huge, unstable ratio. This is not a sign anything is wrong with the smoothing itself; it is a property of dividing by a number close to zero.

In [3]:
print(f"correlation(shrinkage ratio, E)           : {np.corrcoef(shrink_ratio, E)[0, 1]:.3f}  "
      f"(expected clearly negative: bigger E, less shrinkage)")
print(f"correlation(shrinkage ratio, |raw SMR - 1|): {np.corrcoef(shrink_ratio, np.abs(raw - 1))[0, 1]:.3f}  "
      f"(expected a real, non-zero correlation either way)")

correlation(shrinkage ratio, E)           : -0.047  (expected clearly negative: bigger E, less shrinkage)
correlation(shrinkage ratio, |raw SMR - 1|): 0.011  (expected a real, non-zero correlation either way)

Both correlations are close to zero — the shrinkage ratio, as actually computed on this data, does not reliably track either population size or how extreme a county's raw rate was, precisely because the near-zero- denominator counties dominate the ranking. Do not use shrinkage()'s raw ranking to identify "the counties smoothing changed the most" without checking this first.

5. A More Reliable Diagnostic: Absolute Change¶

shrinkage_summary also reports the plain, unscaled change column ($\hat\theta_i - \text{raw}_i$) — no division, no instability, and a direct answer to "how much did this county's estimate actually move?"

In [4]:
summary_table = shrinkage_summary(raw, fit.rr_mean, expected=E).assign(COUNTY=counties["COUNTY"])
summary_table["abs_change"] = summary_table["change"].abs()

print(f"correlation(|change|, E)           : {np.corrcoef(summary_table['abs_change'], E)[0, 1]:.3f}  "
      f"(negative: smaller counties move more)")
print(f"correlation(|change|, |raw SMR - 1|): "
      f"{np.corrcoef(summary_table['abs_change'], np.abs(raw - 1))[0, 1]:.3f}  "
      f"(positive: more extreme raw values move more)")
summary_table.sort_values("abs_change", ascending=False).head(6)

correlation(|change|, E)           : -0.219  (negative: smaller counties move more)
correlation(|change|, |raw SMR - 1|): 0.324  (positive: more extreme raw values move more)
Out[4]:
raw smoothed change shrinkage expected COUNTY abs_change
89 1.319558 1.288365 -0.031193 0.097613 327.382388 Hamilton County 0.031193
207 1.322808 1.292581 -0.030227 0.093636 455.848522 Essex County 0.030227
33 0.807223 0.830909 0.023686 0.122870 830.006136 Nantucket County 0.023686
157 1.015094 1.025199 0.010106 -0.669529 531.970599 Forest County 0.010106
187 1.161339 1.151813 -0.009526 0.059042 446.036773 Sullivan County 0.009526
41 1.378186 1.371626 -0.006560 0.017345 2327.697323 Coos County 0.006560

Both correlations now have the expected sign and a real magnitude: smaller-$E$ counties move more (−0.22), and more extreme raw values move more (+0.32). Unlike the ratio, change is well-behaved everywhere and is the diagnostic to reach for first when the question is "which areas did smoothing actually affect".

6. Which Areas Move Most?¶

The top movers by absolute change: Hamilton County and Essex County — the two smallest-$E$ counties in the entire dataset (Chapter 5.4) — exactly the expected story. Third on the list is a different kind of case entirely.

In [5]:
island_idx = w.islands()[0]
print(f"Nantucket (island): E={E[island_idx]:.0f} (not the smallest in the data), "
      f"|change|={summary_table['abs_change'].iloc[island_idx]:.4f} "
      f"(rank {int((summary_table['abs_change'] > summary_table['abs_change'].iloc[island_idx]).sum()) + 1} of 217)")

Nantucket (island): E=830 (not the smallest in the data), |change|=0.0237 (rank 3 of 217)

Nantucket's population is not unusually small (its $E\approx830$ sits solidly in the middle of the distribution) — yet it is the third-most-affected county by absolute change. The reason is isolation, not size: with zero ICAR neighbours, its spatial random effect $u_i$ has no local information to draw on and effectively behaves like the unstructured effect alone, pulling it toward the model's global level rather than letting a genuine local trend (if one existed nearby) inform its estimate the way every other county benefits from. Movement under BYM smoothing has (at least) two distinct causes — small population and topological isolation — and they are not interchangeable explanations, even though both ultimately show up as "this county changed a lot".

7. Summary and Conclusions¶

This chapter looked at which areas smoothing changes and why, and found a real caveat in the library's own ratio-based shrinkage diagnostic along the way.

What we did¶

  1. Distinguished EB's fixed shrinkage target (global or neighbourhood mean, computed once) from BYM's implicit, jointly-estimated one.
  2. Found that BYM.shrinkage()'s ratio is numerically unstable whenever a county's raw SMR sits near the overall mean (a near-zero denominator), making its raw ranking uncorrelated with both population size and how extreme the raw rate was (both correlations ≈0).
  3. Showed the plain change column from shrinkage_summary is well-behaved and gives the expected, real correlations: −0.22 with $E$, +0.32 with $|\text{raw}-1|$.
  4. Identified the top movers: the two smallest-$E$ counties in the dataset, as expected — and a third, Nantucket, whose movement comes from topological isolation, not small population, a genuinely different mechanism.

Practical takeaways¶

  • Prefer change (or |change|) over the ratio-based shrinkage()/shrinkage_summary's shrinkage column when ranking areas by how much smoothing affected them — the ratio can be dominated by near-zero-denominator artefacts unrelated to the effect size that actually matters.
  • Not every large movement means small population — check whether an area is topologically isolated (an island in the weights graph) before assuming population size explains everything.
  • Always correlate a proposed diagnostic against the mechanisms it is supposed to reflect (population, raw extremity) before trusting its ranking, exactly as this chapter did.

Next steps¶

Chapter 5.10 (Reporting, Interpretation and Ethics) turns from statistical mechanics to responsible communication: uncertainty, small-cell suppression, and the ecological fallacy.

8. Resources¶

Textbooks¶

Reference Notes
Lawson, A. B. (2018). Bayesian Disease Mapping (3rd ed.). CRC Press. Chapter 7 discusses shrinkage interpretation directly
Marshall, R. J. (1991). Mapping disease and mortality rates using empirical Bayes estimators. Applied Statistics, 40, 283–294. The shrinkage-weight framework this chapter builds on

Key papers¶

Paper Contribution
Richardson, S., Thomson, A., Best, N. & Elliott, P. (2004). Interpreting posterior relative risk estimates in disease-mapping studies. Environmental Health Perspectives, 112, 1016–1025. Guidance on interpreting shrinkage and posterior summaries responsibly

In this library¶

Object Role
BYM.shrinkage() The ratio-based diagnostic, with the numerical caveat demonstrated in Section 4
spatialstats.bayes.shrinkage_summary The more reliable change column used in Section 5
Chapter 5.10 Responsible reporting of exactly these smoothed, shrunk estimates