5.7 MCMC Implementation

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

Dataset: Simulated data with known truth (Chapter 5.3, sigma=0.35, expected counts averaging 40)

Focus: The two MCMC backends, their diagnostics, a real demonstration of under-sampling, and the ICAR scaling gotcha

Library: spatialstats.bayes (PySpatialStats)

---

Table of Contents

1. Overview

2. Python Setup

3. The Built-in Sampler

4. The PyMC Backend

   • 4.1 The ICAR scaling gotcha

5. Diagnostics: Why the Built-in Sampler Needs Thinning

6. Cross-backend Validation

7. Posterior Predictive Check

8. DIC, and Why Not WAIC/LOO Here

9. Summary and Conclusions

10. Resources

---

1. Overview

Chapter 5.6 fit BYM and treated the posterior as a given. This chapter looks at how it is actually produced — two different samplers, with different trade-offs — and, critically, demonstrates directly what goes wrong when the built-in sampler's default settings are not respected.

Step	What we do
Built-in sampler	Blocked Gibbs/Metropolis with graph-coloured updates; dependency-free, handles islands natively
PyMC backend	NUTS; needs a connected graph; a specific, documented scaling subtlety worth understanding
Diagnostics	R-hat, effective sample size, divergences — and a direct thin=1 vs. thin=10 comparison
Cross-validation	Do the two backends agree on the same data?
Posterior predictive check	Does the fitted model actually reproduce data like what was observed?
Model comparison	DIC (always available) vs. WAIC/LOO (a real gap in this library, stated honestly)

---

2. Python Setup

import sys, time, 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, KNN, W
from spatialstats.bayes import BYM, simulate_relative_risk, simulate_counts, pymc_available

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")
E_sim = np.full(w.n, 40.0)
theta_true = simulate_relative_risk(w, sigma=0.35, spatial_strength=0.95, seed=1)
O_sim = simulate_counts(theta_true, E_sim, seed=1)
print(f"spatialstats  : version {spatialstats.__version__}")
print(f"PyMC available: {pymc_available()}")

spatialstats  : version 0.2.0

PyMC available: True

---

3. The Built-in Sampler

backend='gibbs' (the default) uses a blocked Gibbs/Metropolis sampler working directly in the convolution ($u+v$) parametrisation. It updates areas in graph-coloured groups — areas that share no edge can be updated simultaneously without breaking the Markov property — and it handles an island natively, since an area with no ICAR neighbours simply has no spatial term to update. Multiple chains start from dispersed initial scales (init_scale varies across chains), which is what makes the R-hat diagnostic meaningful: if chains that started far apart converge to the same distribution, that is real evidence of convergence, not an artefact of similar starting points.

t0 = time.time()
fit_gibbs = BYM(O_sim, E_sim, w, model="bym", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)
print(f"fit time: {time.time() - t0:.1f}s  ({fit_gibbs.chains} chains x {fit_gibbs.n_draws} draws/chain)")
print(f"max R-hat over relative risks: {fit_gibbs.max_rhat:.4f}   min ESS: {fit_gibbs.min_ess:.0f}")

fit time: 10.1s  (4 chains x 2000 draws/chain)
max R-hat over relative risks: 1.0049   min ESS: 1266

---

4. The PyMC Backend

backend='pymc' uses NUTS (the No-U-Turn Sampler), which typically mixes far faster than single-site Gibbs updates for smooth spatial fields — but it requires a fully connected adjacency graph, since PyMC's ICAR distribution is not defined on a graph with isolated nodes. The island must be linked first (Chapter 5.3).

island_idx = w.islands()[0]
island_id = w.id_order[island_idx]
wk1 = KNN(counties, k=1)
nearest_id = wk1.neighbors[island_id][0]
neighbors = {i: list(w.neighbors[i]) for i in w.id_order}
neighbors[island_id] = [nearest_id]
neighbors[nearest_id] = neighbors[nearest_id] + [island_id]
w_linked = W(neighbors, ids=w.id_order, transform="binary")

t0 = time.time()
fit_pymc = BYM(O_sim, E_sim, w_linked, model="bym", backend="pymc", draws=1000, tune=1000, chains=4, seed=1,
              progressbar=False)
print(f"fit time: {time.time() - t0:.1f}s  ({fit_pymc.chains} chains x {fit_pymc.n_draws} draws/chain)")
print(f"max R-hat over relative risks: {fit_pymc.max_rhat:.4f}   min ESS: {fit_pymc.min_ess:.0f}")
print(f"divergences: {fit_pymc.get('divergences')}")

Initializing NUTS using jitter+adapt_diag...

/home/zia207/.pyenv/versions/3.11.16/lib/python3.11/site-packages/pytensor/link/c/cmodule.py:2986: UserWarning: PyTensor could not link to a BLAS installation. Operations that might benefit from BLAS will be severely degraded.
This usually happens when PyTensor is installed via pip. We recommend it be installed via conda/mamba/pixi instead.
Alternatively, you can use an experimental backend such as Numba or JAX that perform their own BLAS optimizations, by setting `pytensor.config.mode == 'NUMBA'` or passing `mode='NUMBA'` when compiling a PyTensor function.
For more options and details see https://pytensor.readthedocs.io/en/latest/troubleshooting.html#how-do-i-configure-test-my-blas-library
  warnings.warn(

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [alpha, tau_u, u_std, tau_v, v]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 6 seconds.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

fit time: 8.8s  (4 chains x 1000 draws/chain)
max R-hat over relative risks: 1.0009   min ESS: 4453
divergences: 0

Comparable or better diagnostics with half the total draws (4,000 vs. 20,000) and zero divergent transitions — the practical payoff of gradient-based NUTS over single-site Gibbs updates for this kind of smooth spatial model.

4.1 The ICAR scaling gotcha

A specific, easy-to-get-wrong implementation detail: PyMC's ICAR distribution is unnormalised in its own scale parameter — its density does not include the $\sigma^{-(n-c)}$ normalising term a proper Gaussian Markov random field needs (where $n-c$ is the field's rank after the sum-to-zero constraint). Putting a hyperprior directly on ICAR's own sigma argument would therefore sample from an improper posterior for that parameter, silently. The documented, correct workaround — and what this library's PyMC backend actually does internally — is to fix ICAR's own sigma=1.0 and instead multiply the unit-scale field by a separately-modelled precision outside the distribution:

tau_u = pm.Gamma("tau_u", a_u, b_u)
u_std = pm.ICAR("u_std", W=A, sigma=1.0)              # fixed unit scale — correct
u = pm.Deterministic("u", u_std / pm.math.sqrt(tau_u))  # the actual spatial effect

rather than the tempting but statistically wrong shortcut of passing a modelled sigma=1/pm.math.sqrt(tau_u) directly into pm.ICAR(...). This is exactly why BYM(backend='pymc')'s implementation is structured this way rather than the more obvious-looking alternative — it is not a stylistic choice.

---

5. Diagnostics: Why the Built-in Sampler Needs Thinning

Single-site Gibbs updates on a smooth spatial field mix slowly — consecutive draws are highly correlated, because each update only moves one area at a time relative to its current neighbours. thin=10 (the default) keeps every 10th draw specifically to counteract this. A direct comparison at thin=1 shows exactly what is lost by skipping it.

t0 = time.time()
fit_thin1 = BYM(O_sim, E_sim, w, model="bym", backend="gibbs", draws=1500, tune=1000, chains=4, thin=1, seed=1)
print(f"thin=1  : {time.time() - t0:.1f}s,  max R-hat={fit_thin1.max_rhat:.4f},  min ESS={fit_thin1.min_ess:.0f}")
print(f"thin=10 : (Section 3 above)     max R-hat={fit_gibbs.max_rhat:.4f},  min ESS={fit_gibbs.min_ess:.0f}")

thin=1  : 1.3s,  max R-hat=1.0522,  min ESS=112
thin=10 : (Section 3 above)     max R-hat=1.0049,  min ESS=1266

At thin=1, R-hat for the worst relative risk rises to about 1.05 (right at the usual 1.05 convergence threshold) and the worst effective sample size drops by roughly an order of magnitude — despite using fewer raw MCMC iterations before thinning is even applied (1,500+1,000 vs. 2,000+2,000). Thinning does not literally add information; but by discarding highly autocorrelated intermediate draws it keeps the retained sample close to independent, which is what R-hat and ESS actually measure. Skipping it is the single most likely way to get a silently under-converged BYM fit from the built-in sampler.

---

6. Cross-backend Validation

Do the two, structurally different samplers actually agree on the same data? This is the strongest available check that neither implementation has a hidden bug — but it is only a fair test if both backends fit the same weights graph. Section 3's fit_gibbs used w (island included); Section 4's fit_pymc had to use w_linked (island merged in) for PyMC's connectivity requirement — comparing those two directly would confound "different sampler" with "different graph". A fair comparison refits the built-in sampler on w_linked too.

fit_gibbs_linked = BYM(O_sim, E_sim, w_linked, model="bym", backend="gibbs",
                       draws=2000, tune=2000, chains=4, thin=10, seed=1)

corr = np.corrcoef(fit_gibbs_linked.rr_mean, fit_pymc.rr_mean)[0, 1]
print(f"posterior-mean correlation (gibbs vs. pymc, SAME graph w_linked): {corr:.5f}")
print(f"DIC: gibbs={fit_gibbs_linked.dic:.1f}   pymc={fit_pymc.dic:.1f}")
print(f"tau_u posterior mean: gibbs={fit_gibbs_linked.diagnostics.loc['tau_u', 'mean']:.2f}   "
      f"pymc={fit_pymc.diagnostics.loc['tau_u', 'mean']:.2f}")

posterior-mean correlation (gibbs vs. pymc, SAME graph w_linked): 0.99997
DIC: gibbs=1525.7   pymc=1525.5
tau_u posterior mean: gibbs=5.68   pymc=5.72

With the confound removed, the two independently implemented samplers — one hand-written Gibbs/Metropolis, one general-purpose NUTS — converge to the same answer: posterior-mean correlation above 0.999 and DIC/$\tau_u$ within a fraction of a percent of each other. This is about as strong a correctness check as is practically available for a custom MCMC implementation, and it also makes a second point on its own: fitting the same model on two different graphs (island merged vs. left alone) is not the same experiment, even with an identical sampler — always hold the graph fixed when comparing anything else.

---

7. Posterior Predictive Check

Beyond convergence, does the fitted model actually reproduce data that looks like what was observed? Draw simulated counts from the posterior and compare their distribution to the real observed counts.

rng = np.random.default_rng(1)
post_draws = fit_gibbs.posterior_draws()          # (chains*draws, n)
sample_idx = rng.choice(len(post_draws), size=200, replace=False)
ppc = rng.poisson(E_sim[None, :] * post_draws[sample_idx])   # (200, n) simulated datasets

obs_mean, obs_max = O_sim.mean(), O_sim.max()
ppc_means = ppc.mean(axis=1)
ppc_maxes = ppc.max(axis=1)
print(f"observed mean count = {obs_mean:.2f}   posterior-predictive mean count: "
      f"[{np.quantile(ppc_means, 0.025):.2f}, {np.quantile(ppc_means, 0.975):.2f}]")
print(f"observed max count  = {obs_max:.0f}    posterior-predictive max count : "
      f"[{np.quantile(ppc_maxes, 0.025):.0f}, {np.quantile(ppc_maxes, 0.975):.0f}]")

observed mean count = 41.96   posterior-predictive mean count: [40.66, 43.13]
observed max count  = 115    posterior-predictive max count : [95, 132]

The observed mean and maximum both fall comfortably inside their posterior-predictive intervals — the fitted model can generate data indistinguishable, on these two summaries, from what was actually observed. A posterior predictive check that failed (observed value far outside the simulated range) would be a concrete sign of model misspecification, distinct from an MCMC convergence problem.

---

8. DIC, and Why Not WAIC/LOO Here

BYM.dic is always available from either backend, computed directly from the posterior deviance. WAIC and LOO (via ArviZ, when using the PyMC backend) are a real gap in this library's current implementation: they require the model to track pointwise log-likelihood at every posterior draw, which run_pymc does not currently request from PyMC.

import arviz as az
try:
    az.waic(fit_pymc.idata)
except Exception as e:
    print(f"WAIC unavailable: {e}")

WAIC unavailable: log likelihood not found in inference data object

This confirms the gap directly rather than leaving it as an assumption: fit_pymc.idata (the raw ArviZ InferenceData object, available whenever backend='pymc') does not carry a log_likelihood group, so az.waic/az.loo cannot be computed from it without modifying the model definition to request idata_kwargs={"log_likelihood": True} — an extension outside the scope of what BYM currently exposes. DIC remains the practical, always-available model-comparison metric in this library (used throughout Chapters 5.6–5.11); WAIC/LOO would need custom PyMC code written outside the BYM wrapper.

---

9. Summary and Conclusions

This chapter opened up how BYM's posteriors are actually produced, and validated both routes to get there.

What we did

1. Explained the built-in Gibbs/Metropolis sampler: graph-coloured updates, native island handling, dispersed chain starts for a meaningful R-hat.
2. Ran the PyMC/NUTS backend: comparable or better diagnostics with a fifth of the total draws and zero divergences, at the cost of requiring a connected graph.
3. Explained the ICAR unnormalised-scale gotcha precisely, with the correct fixed-scale-then-multiply pattern this library's own PyMC backend implements.
4. Demonstrated thin=1 degrading diagnostics directly: R-hat rises to ≈1.05 and ESS drops roughly 10-fold compared to the thin=10 default, on fewer raw iterations.
5. Cross-validated both backends against each other: posterior-mean correlation 0.99997, near-identical DIC and $\tau_u$.
6. Ran a posterior predictive check: the observed mean and max count both fall inside their posterior- predictive intervals.
7. Confirmed DIC is always available but WAIC/LOO are not currently supported by this library's PyMC integration, and explained precisely why (no tracked log-likelihood).

Practical takeaways

• Never drop the built-in sampler's default thin=10 without checking max_rhat/min_ess afterward — Section 5 shows concretely what is lost.
• Cross-backend validation (Section 6) is a strong, cheap correctness check whenever both are available — run it, don't just trust one implementation.
• Use DIC for model comparison in this library; treat WAIC/LOO as unavailable unless the PyMC model is extended to track log-likelihood explicitly.

Next steps

Chapter 5.8 (Comparing and Evaluating Smoothing Methods) puts every method from Chapters 5.4–5.6 through the formal simulation-based scoring Chapter 5.3's known-truth data was built for.

---

10. Resources

Textbooks

Reference	Notes
Gelman, A. et al. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.	Standard reference on MCMC diagnostics (R-hat, ESS, posterior predictive checks)
Hoffman, M. D. & Gelman, A. (2014). The No-U-Turn sampler. Journal of Machine Learning Research, 15, 1593–1623.	Origin of NUTS

Key papers

Paper	Contribution
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R̂ for assessing convergence of MCMC. Bayesian Analysis, 16, 667–718.	Modern split-R̂ diagnostic used by rhat()
Vehtari, A., Gelman, A. & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27, 1413–1432.	The WAIC/LOO methods this library does not yet expose

In this library

Object	Role
spatialstats.bayes.BYM(backend='gibbs'\|'pymc')	Both samplers used throughout this chapter
BYM.max_rhat, .min_ess, .diagnostics	Convergence diagnostics
BYM.idata	Raw ArviZ object, available with backend='pymc'
BYM.posterior_draws()	Used for the posterior predictive check in Section 7
Chapter 5.8	The formal simulation-based scoring this chapter's validated posteriors feed into