Banner

5.11 Hand-offs and Exercises¶

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

Dataset: diabetes_northeast.shp (217 counties); simulated data with known truth

Focus: What carries forward to Parts 6 and 9, and five solved exercises spanning this Part

Library: spatialstats.bayes, spatialstats.gwmodel (PySpatialStats)

Table of Contents¶

  1. Overview

  2. Python Setup

  3. Hand-offs to Parts 6 and 9

  4. Exercises

    • 4.1 A genuinely rare outcome: EB vs. raw SMR

    • 4.2 The weight on the data as E → ∞

    • 4.3 Rerunning BYM with thin=1

    • 4.4 Adding Obesity as a covariate

    • 4.5 Comparing model='bym', 'icar', and 'iid' by DIC

  5. Summary and Conclusions

  6. Resources

1. Overview¶

Part 5 closes with two things: an honest account of what its methods cannot do (and where the library that can do it lives), and five solved exercises that revisit specific claims made across Chapters 5.2–5.9 with fresh, independent checks.

Step What we do
Hand-offs Why GWGLM (Part 6) cannot replace this Part's expected-count-adjusted models, and why space-time BYM (Part 9's natural next step) is not implemented
Exercises Five solved problems spanning Chapters 5.2 (EB vs. raw), 5.5 (the shrinkage weight), 5.7 (thin=1), 5.6 (covariates), and DIC-based model comparison

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, eb_gamma_poisson, BYM, simulate_relative_risk, simulate_counts, evaluate_smoothing,
)

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

spatialstats : version 0.2.0

3. Hand-offs to Parts 6 and 9¶

To Part 6 (geographically weighted models): GWGLM(family='poisson') fits a locally varying Poisson regression — a natural-looking way to let disease risk's relationship with covariates vary smoothly over space, extending Part 3's spatial regression the way Part 3, Chapter 3.6 previewed. But GWGLM accepts no offset argument (confirmed directly: its constructor takes family, bandwidth, and standard GWR arguments, with no place for $\log E_i$), which means it cannot reproduce this Part's $O_i \sim \text{Poisson}(E_i\theta_i)$ structure — it can only model raw counts as a function of covariates, not counts relative to an appropriately-varying population base. GWGLM is not a drop-in replacement for BYM; it answers a different question (how does a covariate relationship vary locally) and needs its own population-adjustment handled separately if disease-mapping-style analysis is the goal.

To Part 9 (space-time methods): every model in this Part is a single cross-section — one snapshot in time. Chapter 4.7 already previewed the natural next question (are the same areas hot spots year after year, or does the pattern move?) with a simple two-snapshot Gi* comparison. A proper space-time BYM extension — jointly modelling spatial and temporal random effects, letting strength be borrowed across both dimensions at once — is not implemented in this library; Part 9's own tools are the place to look for it.

4. Exercises¶

4.1 A genuinely rare outcome: EB vs. raw SMR¶

Simulate a rare outcome and show EB beats the SMR in RMSE.

In [2]:
theta_rare_true = simulate_relative_risk(w, sigma=0.35, spatial_strength=0.95, seed=5)
E_rare = np.full(w.n, 8.0)     # genuinely rare: expected count of 8, not this Part's usual 40+
O_rare = simulate_counts(theta_rare_true, E_rare, seed=5)
print(f"simulated rare counts: min={O_rare.min()}, max={O_rare.max()}, mean={O_rare.mean():.1f}")

raw_rare = smr(O_rare, E_rare)
eb_rare = eb_gamma_poisson(O_rare, E_rare, method="mle")
ev = evaluate_smoothing(theta_rare_true, {"Raw SMR": raw_rare, "EB gamma-Poisson": eb_rare.estimate})
ev.round(4)

simulated rare counts: min=1, max=22, mean=8.1
Out[2]:
rmse mae bias spearman
Raw SMR 0.3639 0.2885 -0.0467 0.6962
EB gamma-Poisson 0.2695 0.2096 -0.0467 0.6962

Answer: at this genuinely rare-outcome scale (mean count ≈8), EB's RMSE (0.270) is about 26% lower than the raw SMR's (0.364) — a far larger gap than Chapter 5.5 found on the real, non-sparse diabetes data (where EB's improvement over raw was modest). The smaller and rarer the counts, the more empirical Bayes has to offer.

4.2 The weight on the data as E → ∞¶

Show that the weight on the data tends to 1 as E → ∞.

In [3]:
a, b = eb_rare.get("shape"), eb_rare.get("rate")
print(f"fitted prior: a={a:.3f}, b={b:.3f}")
for E_test in (1, 10, 100, 1_000, 10_000, 100_000, 10_000_000):
    weight = E_test / (E_test + b)
    print(f"E={E_test:>12,.0f}: weight on the data = {weight:.6f}")

fitted prior: a=7.928, b=7.811
E=           1: weight on the data = 0.113490
E=          10: weight on the data = 0.561439
E=         100: weight on the data = 0.927546
E=       1,000: weight on the data = 0.992249
E=      10,000: weight on the data = 0.999219
E=     100,000: weight on the data = 0.999922
E=  10,000,000: weight on the data = 0.999999

Answer: $w = E/(E+b)$ rises monotonically from 0.11 at $E=1$ to essentially 1.000000 by $E=10^7$ — directly verifying the formula's asymptotic behaviour. Interpreted plainly: as a county's population (and hence its expected count) grows without bound, its raw SMR eventually becomes precise enough that the prior contributes nothing at all, and EB's estimate converges to the raw SMR exactly — exactly the "large counties barely move" behaviour observed empirically back in Chapter 5.1's very first look.

4.3 Rerunning BYM with thin=1¶

This was worked out directly in Chapter 5.7, Section 5, on the same simulated data used throughout Chapters 5.7–5.8: at thin=1 (1,500 draws, 1,000 tuning), the worst R-hat across all 217 relative risks rose to about 1.05 (right at the standard convergence threshold) and the worst effective sample size dropped to about 112 — roughly a 10-fold degradation from the thin=10 default's R-hat of 1.005 and ESS of 1,266, despite thin=1 using fewer total iterations before thinning. What goes wrong: the built-in sampler's single-site updates are highly autocorrelated from one raw iteration to the next for a smooth spatial field, so without thinning to counteract it, the retained draws are far from independent — R-hat and ESS (which assume approximately independent draws) both correctly flag this as unreliable convergence evidence, even though the sampler itself has not diverged or produced wrong values on average.

4.4 Adding Obesity as a covariate¶

Also already worked out directly, in Chapter 5.6, Section 5, on the real diabetes data: adding Obesity gives a precisely estimated positive coefficient (0.0247, 95% CI [0.021, 0.029] — comfortably excluding 0), and the spatial fraction drops from 0.47 to 0.37. The smoothed relative risks themselves shift modestly in the direction obesity would predict for each county, but the more diagnostic change is the spatial fraction: about a fifth of what looked like unexplained spatial structure in the covariate-free model is, once obesity enters, explained by obesity's own (substantially spatially patterned) distribution rather than left for the random effects to absorb.

4.5 Comparing model='bym', 'icar', and 'iid' by DIC¶

In [4]:
fit_bym = BYM(O, E, w, model="bym", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)
fit_icar = BYM(O, E, w, model="icar", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)
fit_iid = BYM(O, E, w, model="iid", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)

print(f"BYM  (spatial + unstructured): DIC = {fit_bym.dic:.2f}")
print(f"ICAR (spatial only)         : DIC = {fit_icar.dic:.2f}")
print(f"IID  (unstructured only)    : DIC = {fit_iid.dic:.2f}")

BYM  (spatial + unstructured): DIC = 2806.34
ICAR (spatial only)         : DIC = 2861.73
IID  (unstructured only)    : DIC = 2805.70

Answer: on the real (non-sparse) diabetes data, IID has the lowest DIC (2805.7), essentially tied with BYM (2806.3), while pure ICAR is clearly worst (2861.7). This is a genuinely informative, slightly surprising result: it does not mean diabetes prevalence has no spatial structure — Part 2 and this Part's own Moran's I checks (Chapter 5.4) confirm real, significant autocorrelation. It means that, given how well-measured these particular counts already are (minimum $E \approx 327$, Chapter 5.3), the extra flexibility of a combined spatial-plus-unstructured model does not pay for its added complexity in DIC's effective-parameter accounting — a simpler unstructured model fits nearly as well once the data are this precise. The scenario looks different for genuinely sparse counts: Chapter 5.8's simulated, $E=40$ data (Section 4.1's even sparser $E=8$ scenario more so still) is exactly where BYM's spatial borrowing is expected to earn a clearer DIC advantage over IID alone — this real-data result is a reminder that DIC-preferred model complexity depends on how much the data can already tell you on their own, not a fixed ranking of the three model families.

5. Summary and Conclusions¶

This chapter closed Part 5 with its scope boundaries and five independently verified exercises.

What we did¶

  1. Confirmed GWGLM has no offset support and so cannot replace this Part's expected-count-adjusted models; noted space-time BYM is not implemented and belongs to Part 9.
  2. Exercise 1: at a genuinely rare-outcome scale, EB's RMSE beat raw SMR's by about 26% — a far bigger gap than on the real, well-measured diabetes data.
  3. Exercise 2: verified the shrinkage weight $E/(E+b)$ rises from 0.11 to ≈1.000000 as $E$ grows from 1 to $10^7$.
  4. Exercise 3: restated Chapter 5.7's direct demonstration that thin=1 degrades R-hat and ESS by roughly an order of magnitude relative to the thin=10 default.
  5. Exercise 4: restated Chapter 5.6's finding that adding Obesity drops the spatial fraction from 0.47 to 0.37.
  6. Exercise 5: compared bym/icar/iid by DIC on real data, finding IID and BYM essentially tied and both ahead of pure ICAR — a genuinely informative result about how much spatial complexity well-measured data can support, not a fixed ranking.

Part 5 in review¶

Across eleven chapters, Part 5 took Part 2 and Part 4's implicit assumption — that every county's reported rate is equally trustworthy — and replaced it with a full, validated toolkit for when it is not: the mechanism (Poisson variance, Chapter 5.2), two families of remedies (empirical Bayes, Chapter 5.5; fully Bayesian BYM, Chapters 5.6–5.7), formal validation against a known truth (Chapter 5.8), per-area interpretation (Chapter 5.9), and the reporting responsibilities that come with publishing any of it (Chapter 5.10). Part 6 (geographically weighted models) picks up the "one global relationship" limitation Part 3 flagged, with the important caveat (Section 3) that it does not inherit this Part's population-adjustment machinery automatically; Part 9 (space-time methods) is where the single-snapshot limitation shared by every model in this Part gets lifted.

6. Resources¶

Textbooks¶

Reference Notes
Lawson, A. B. (2018). Bayesian Disease Mapping (3rd ed.). CRC Press. Chapter 10 covers space-time extensions, the Part 9 hand-off
Fotheringham, A. S., Brunsdon, C. & Charlton, M. (2002). Geographically Weighted Regression. Wiley. The Part 6 hand-off, and GWGLM's scope

Key papers¶

Paper Contribution
Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19, 2555–2567. Foundational space-time disease-mapping model, the Part 9 target

In this library¶

Object Role
spatialstats.bayes.* Every object exercised across Part 5, summarised one final time in this chapter
spatialstats.gwmodel.GWGLM The Part 6 model whose offset limitation this chapter states explicitly
Part 6 Geographically weighted models — the next Part
Part 9 Space-time methods — where this Part's single-snapshot limitation is lifted