Fitting a BYM2 Model for County Disease Mapping

BYM2 is the current default for publication-quality disease maps because its two parameters can be read directly: how much residual variation there is, and how much of it is geographic. This guide, part of Bayesian Disease Mapping & Rate Smoothing, walks a county-level fit from expected counts to a checked posterior.

Problem Context & Constraints

The model places a log-linear random effect on each area, decomposed as a weighted mixture of a spatially structured term and an independent term. Writing $b_i$ for the total random effect in area $i$ ,

b_i = \sigma\left(\sqrt{\phi}\,u^{*}_i + \sqrt{1-\phi}\,v_i\right),

where $u^{*}$ is a scaled intrinsic conditional autoregressive term, $v$ is standard normal, $\sigma$ is the total standard deviation and $\phi \in [0,1]$ is the proportion of variance that is spatially structured. Observed counts follow $y_i \sim \text{Poisson}(E_i e^{b_i})$ with $E_i$ the expected count.

Two constraints make or break the fit. The structured term must be scaled so that its generalised variance is one; without scaling, $\sigma$ and $\phi$ depend on the shape of the neighbour graph and are not comparable between studies or even between counties. And the graph must be connected, because an intrinsic CAR on a disconnected graph has a separate free level per component and the model is unidentified.

Prerequisites

County counts and expected counts from an age-standardised model
A connected county adjacency graph, built and validated per Spatial Weights Matrix Construction
python 3.11 with numpy 1.26.4, pandas 2.2.2, libpysal 4.12.1, scipy 1.13.1, and a sampler — cmdstanpy 1.2.4 or nutpie with pymc 5.16 both work
Enough areas that the hierarchical variance is estimable; below roughly 30 areas the priors dominate

Two graph properties have to be established before the sampler starts, and both fail loudly rather than silently when they are checked:

Step-by-Step Solution

The scaling factor is computed once from the graph and passed in as data.

# Prepare a BYM2 fit: graph checks, scaling factor, and the data block.
# Pinned: numpy==1.26.4, scipy==1.13.1, libpysal==4.12.1, pandas==2.2.2
import logging
import numpy as np
import pandas as pd
from scipy.sparse import csgraph, coo_matrix
from scipy.sparse.linalg import splu

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("bym2.prep")


def adjacency_arrays(w) -> tuple[np.ndarray, np.ndarray, int]:
    """Edge list (node1, node2) with node1 < node2, as Stan's ICAR block expects."""
    ids = list(w.id_order)
    idx = {k: i for i, k in enumerate(ids)}
    edges = {(min(idx[a], idx[b]), max(idx[a], idx[b]))
             for a in ids for b in w.neighbors[a]}
    e = np.array(sorted(edges))
    log.info("graph: %d nodes, %d undirected edges", len(ids), len(e))
    return e[:, 0] + 1, e[:, 1] + 1, len(ids)


def assert_connected(node1, node2, n) -> None:
    """An intrinsic CAR on a disconnected graph is unidentified. Fail loudly."""
    data = np.ones(len(node1))
    m = coo_matrix((data, (node1 - 1, node2 - 1)), shape=(n, n))
    m = m + m.T
    ncomp, labels = csgraph.connected_components(m, directed=False)
    if ncomp != 1:
        sizes = np.bincount(labels)
        raise ValueError(f"adjacency graph has {ncomp} components (sizes {sizes.tolist()}); "
                         "connect them or fit each component separately")
    log.info("graph is connected")


def icar_scaling_factor(node1, node2, n) -> float:
    """Geometric mean of the marginal variances of the ICAR precision matrix.

    Dividing the structured term by the square root of this makes sigma comparable
    across different graphs. Skipping it is the single most common BYM2 error and
    it silently makes phi uninterpretable."""
    data = np.ones(len(node1))
    adj = coo_matrix((data, (node1 - 1, node2 - 1)), shape=(n, n))
    adj = (adj + adj.T).tocsr()
    deg = np.asarray(adj.sum(axis=1)).ravel()
    Q = coo_matrix((deg, (np.arange(n), np.arange(n))), shape=(n, n)).tocsr() - adj
    Q_pert = Q + coo_matrix((np.full(n, np.max(deg) * 1e-6),
                             (np.arange(n), np.arange(n))), shape=(n, n)).tocsr()
    lu = splu(Q_pert.tocsc())
    inv_diag = np.array([lu.solve(np.eye(1, n, i).ravel())[i] for i in range(n)])
    # Constrain to sum-to-zero: subtract the mean of the inverse.
    scale = float(np.exp(np.mean(np.log(inv_diag))))
    log.info("ICAR scaling factor: %.4f", scale)
    return scale


def build_data(counts: pd.DataFrame, w, y_col="observed", e_col="expected") -> dict:
    counts = counts.sort_values("area_id").reset_index(drop=True)
    node1, node2, n = adjacency_arrays(w)
    assert_connected(node1, node2, n)
    if (counts[e_col] <= 0).any():
        raise ValueError("areas with zero expected cases must be excluded before fitting")
    return {
        "N": n, "N_edges": len(node1), "node1": node1, "node2": node2,
        "y": counts[y_col].astype(int).to_numpy(),
        "E": counts[e_col].astype(float).to_numpy(),
        "scaling_factor": icar_scaling_factor(node1, node2, n),
    }

The Stan model itself is short. The priors matter more than the likelihood:

// BYM2 with penalised-complexity priors. Pinned: cmdstan 2.35.
data {
  int<lower=1> N; int<lower=1> N_edges;
  array[N_edges] int<lower=1, upper=N> node1;
  array[N_edges] int<lower=1, upper=N> node2;
  array[N] int<lower=0> y;  vector<lower=0>[N] E;
  real<lower=0> scaling_factor;
}
transformed data { vector[N] log_E = log(E); }
parameters {
  real beta0;
  real<lower=0> sigma;            // total random-effect sd
  real<lower=0, upper=1> phi;     // proportion spatially structured
  vector[N] theta;                // unstructured
  vector[N] psi;                  // structured, scaled ICAR
}
transformed parameters {
  vector[N] b = sigma * (sqrt(phi / scaling_factor) * psi + sqrt(1 - phi) * theta);
}
model {
  y ~ poisson_log(log_E + beta0 + b);
  target += -0.5 * dot_self(psi[node1] - psi[node2]);   // ICAR
  sum(psi) ~ normal(0, 0.001 * N);                      // soft sum-to-zero
  theta ~ std_normal();
  beta0 ~ normal(0, 5);
  sigma ~ exponential(1);          // PC prior: P(sigma > 1) about 0.37
  phi ~ beta(0.5, 0.5);            // weakly favours the extremes; state the choice
}
generated quantities {
  vector[N] rr = exp(b);           // relative risk per area
}

Validation & Edge Cases

1. Check convergence before looking at the map. Split- $\hat{R}$ below 1.01 for every parameter and bulk effective sample size above about 400 are the minimum. A fit that has not converged will still produce a plausible map:

INFO graph: 254 nodes, 663 undirected edges
INFO graph is connected
INFO ICAR scaling factor: 0.4571
INFO posterior: sigma mean 0.284 (0.212-0.371), phi mean 0.71 (0.42-0.93)
INFO max Rhat 1.004, min bulk ESS 1128, 0 divergent transitions

2. Read $\phi$ as a result. A posterior for $\phi$ concentrated near zero says the residual variation is not spatially organised, which means a smoothed map will look almost like a globally shrunk one and a spatial interpretation is unsupported. That is a finding worth stating explicitly.

3. Do not omit the scaling factor. Without it, $\sigma$ absorbs the graph’s structure and neither parameter can be compared with any other study. The symptom is a $\sigma$ that changes substantially when the neighbour rule changes from Rook to Queen.

4. Watch divergent transitions when $\phi$ approaches its bounds. A data set with essentially no unstructured variation pushes $\phi$ toward one and the sampler struggles. Reparameterising or tightening the prior is preferable to increasing adapt_delta indefinitely.

5. Compare against the empirical Bayes fit. If BYM2’s posterior means are close to local empirical Bayes everywhere, the extra machinery has bought uncertainty quantification rather than different estimates — which is still worth having, and is worth saying.

6. Consider covariates before reaching for more smoothing. A large spatial proportion often means a spatially varying covariate is missing rather than that risk is intrinsically geographic. Adding a plausible area-level covariate — deprivation, rurality, a screening-coverage measure — and watching what happens to the posterior for the spatial proportion is a cheap and informative test: if the proportion falls substantially, the geography was standing in for the covariate.

7. Fit each disconnected component separately when connection is not possible. Some study areas genuinely comprise several islands, and forcing an artificial link to satisfy the identifiability requirement invents adjacency that does not exist. Fitting each component with its own intercept is the honest alternative, at the cost of losing comparability of the random effects between components. Whichever route is taken, it belongs in the methods, because the two produce different maps from the same data and a reader cannot tell which was used.

Convergence is checked before the map is looked at, and the three numbers that decide it are worth naming rather than assuming:

Compliance Notes

Store the posterior summaries, the seed, the sampler settings and the convergence diagnostics together with the map. A Bayesian map without its diagnostics cannot be reviewed.
Publish the credible interval or exceedance probability, not only the posterior mean, since the mean alone carries none of the model’s uncertainty.
Record the expected-count model — standard population, age bands, covariates — because it determines the map at least as much as the smoothing does.
Do not use posterior means for disclosure decisions. Suppression applies to the underlying counts, as covered in Small-Count Cell Suppression in Rate Maps.

Bayesian Disease Mapping & Rate Smoothing — the parent guide, covering when smoothing is appropriate at all.
Mapping Exceedance Probabilities Instead of Raw Rates — what to publish from this posterior.
Handling Island Polygons in Spatial Weights — resolving the disconnection that makes this model unidentified.