Choosing Priors for Spatial Random Effects in Disease Maps

In a hierarchical disease map the prior on the random-effect variance is not a technicality — it is the parameter that decides how much smoothing happens, and with small-area counts it can move the map more than the data does. This guide, part of Bayesian Disease Mapping & Rate Smoothing, covers choosing that prior defensibly and demonstrating that the conclusion does not depend on it.

Problem Context & Constraints

The difficulty is that the variance parameter is weakly identified. With a few hundred areas and small counts, the likelihood is nearly flat over a wide range of plausible variances, so the posterior is substantially the prior. A gamma prior on the precision with shape and rate both 0.001 — for decades the default in published disease-mapping code — is not “uninformative”: it places most of its mass at very small variances, which shrinks aggressively, and its behaviour changes with the number of areas.

Penalised-complexity priors were designed for exactly this situation. The idea is to define a base model — here, no random effect at all — and penalise distance from it at a constant rate. For a standard deviation the resulting prior is exponential, and it is set by one interpretable statement: the probability that the standard deviation exceeds some value $U$ is $\alpha$ . Writing “I believe there is a 1% chance the residual relative risk varies by more than a factor of about three” is a defensible sentence; writing “gamma(0.001, 0.001)” is not.

Two constraints govern any choice. The prior must be stated on an interpretable scale — a standard deviation on the log-relative-risk scale, not a precision — and its consequence must be demonstrated, by refitting under alternatives and showing what changes.

Prerequisites

A working hierarchical fit, such as the BYM2 model in Fitting a BYM2 Model for County Disease Mapping
A statement of what magnitude of residual variation is epidemiologically plausible for the condition
python 3.11, numpy 1.26.4, pandas 2.2.2, and a sampler
Compute budget for at least three refits, since the sensitivity analysis is the deliverable

The translation from a sentence to a distribution is the step worth making explicit, because it is where the prior stops being arbitrary:

Step-by-Step Solution

Translate a substantive statement into a prior, then refit under alternatives.

# Turn an interpretable statement into a PC prior, and score prior sensitivity.
# Pinned: numpy==1.26.4, pandas==2.2.2, scipy==1.13.1
import logging
import numpy as np
import pandas as pd

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


def pc_rate_for_sd(u: float, alpha: float = 0.01) -> float:
    """Exponential rate lambda such that P(sigma > u) = alpha.

    `u` is stated on the log-relative-risk scale: u = log(3) says the residual
    relative risk is unlikely to vary by more than threefold. That sentence is
    what goes in the paper, and this function is the only translation step."""
    if not (0 < alpha < 1) or u <= 0:
        raise ValueError("need u > 0 and 0 < alpha < 1")
    lam = -np.log(alpha) / u
    log.info("PC prior: sigma ~ Exponential(%.3f)  [P(sigma > %.3f) = %.2f, "
             "i.e. residual RR spread beyond %.1fx is %.0f%% likely]",
             lam, u, alpha, np.exp(u), 100 * alpha)
    return lam


def prior_predictive_rr(lam: float, n_draw: int = 20000, seed: int = 42) -> dict:
    """What relative risks does this prior imply before seeing data?

    Checking this is what catches a prior that permits absurd maps. A prior whose
    97.5th percentile relative risk is 40 is not weakly informative; it is wrong."""
    rng = np.random.default_rng(seed)
    sigma = rng.exponential(1 / lam, n_draw)
    rr = np.exp(rng.normal(0, sigma))
    q = np.percentile(rr, [2.5, 50, 97.5])
    log.info("prior-predictive RR: 2.5%% %.2f, median %.2f, 97.5%% %.2f", *q)
    return {"q025": q[0], "median": q[1], "q975": q[2]}


def prior_sensitivity(fits: dict, threshold: float = 1.5, p_cut: float = 0.95) -> pd.DataFrame:
    """Compare fits under different priors on what they actually change.

    `fits` maps a prior label to an (n_draws, n_areas) array of relative risks.
    The comparison is on the DECISION, not on the parameter, because a posterior
    for sigma that shifts while the flagged set does not is not a problem."""
    rows, flagged = [], {}
    for name, draws in fits.items():
        p = (draws > threshold).mean(axis=0)
        f = set(np.flatnonzero(p >= p_cut))
        flagged[name] = f
        rows.append({"prior": name, "median_sigma_proxy": float(np.std(np.log(draws), axis=0).mean()),
                     "n_flagged": len(f)})
    base = next(iter(flagged))
    for r in rows:
        f, b = flagged[r["prior"]], flagged[base]
        r["jaccard_vs_base"] = len(f & b) / len(f | b) if (f | b) else 1.0
    out = pd.DataFrame(rows)
    log.info("prior sensitivity:\n%s", out.to_string(index=False))
    return out

The prior_sensitivity function scores agreement on the flagged set, not on the variance parameter, because that is what a reader acts on. A prior that changes $\sigma$ by 20% while leaving the same eleven counties flagged has not changed the finding; one that leaves $\sigma$ alone while flagging a different eleven has.

Validation & Edge Cases

1. Always run the prior predictive check. It takes seconds and it catches the most embarrassing errors — priors that permit relative risks of forty, or that force every area to be identical:

INFO PC prior: sigma ~ Exponential(4.192)  [P(sigma > 1.099) = 0.01, i.e. residual RR spread beyond 3.0x is 1% likely]
INFO prior-predictive RR: 2.5% 0.62, median 1.00, 97.5% 1.63
INFO prior sensitivity:
      prior  median_sigma_proxy  n_flagged  jaccard_vs_base
   pc_u_log3               0.281         11            1.000
   pc_u_log10              0.334         14            0.786
 gamma_0.001               0.108          3            0.273

2. Report the Jaccard agreement, not just the counts. Two priors that both flag eleven areas may flag different elevens. The overlap is the number that says whether the conclusion is prior-driven.

3. Do not use a uniform prior on the variance to appear neutral. A uniform prior on a variance is strongly informative on the standard deviation and vice versa, and neither is neutral. Neutrality on a scale-parameter is not achievable; interpretability is.

4. Beware very small numbers of areas. Below roughly thirty areas the prior dominates whatever it is, and the honest report says so rather than presenting a posterior as if the data had spoken.

5. Keep the prior fixed across a series. A prior re-chosen each year makes the time series incomparable, and any change should be treated like a change of case definition — announced, dated and applied retrospectively where possible.

6. Distinguish prior sensitivity from model sensitivity. A conclusion that survives three priors but collapses when the neighbour rule changes from Rook to Queen is not prior-robust in any useful sense — it is fragile, and the prior sweep gave false comfort. Run at least one alternative weights specification alongside the prior alternatives, and report the two sensitivities together. The combined table is barely longer than either alone and it is much harder to misread.

7. Keep the sweep cheap enough to actually run. Three full MCMC fits on a county-level model is minutes; on a tract-level model with several thousand areas it can be hours, and a sensitivity analysis that nobody runs because it is expensive is worse than one designed to be affordable. An integrated-nested-Laplace approximation fit is fast enough to sweep freely and close enough to the MCMC posterior for this purpose, so use it for the sweep and reserve full sampling for the final reported fit. State which was used for which.

The sensitivity result is best read as an overlap rather than as a count, because two priors can flag the same number of different areas:

Compliance Notes

State the prior as a sentence, with its numeric form beside it. “We assumed a 1% prior probability that residual relative risk varies by more than threefold ( $\sigma \sim \text{Exp}(4.19)$ )” is reviewable; a bare distribution is not.
Publish the sensitivity table. It is short, it is cheap to produce, and it is the only evidence that the map is not an artefact of a default.
Record the prior in the run signature alongside the seed and the expected-count model, since a re-run under a different prior is a different analysis.
Do not report the prior-driven areas as findings when the sensitivity analysis shows they appear only under one prior.

Bayesian Disease Mapping & Rate Smoothing — the parent guide covering when to smooth at all.
Fitting a BYM2 Model for County Disease Mapping — the model these priors go into.
Mapping Exceedance Probabilities Instead of Raw Rates — the decision quantity the sensitivity analysis should be run on.