Spatial Epi Public Health GIS

Bayesian Disease Mapping & Rate Smoothing

A map of raw small-area rates is dominated by the places with the fewest people. One extra case in a tract of 400 residents produces a rate that towers over everything around it, and the resulting map shows population sparsity rather than disease. This guide is part of Disease Clustering & Spatial Statistical Modeling, and it covers the smoothing methods that fix this, what they cost, and how to report the result.

Concept & Epidemiological Alignment

The problem is variance, not bias. A rate estimated from a small denominator is unbiased and enormously variable, so the extremes of a raw rate map are almost entirely the small-denominator units, in both directions. Ranking areas by raw rate ranks them approximately by how few people live in them.

Smoothing methods address this by borrowing strength: each area’s estimate is pulled toward a mean, and the amount of pull depends on how little information the area itself carries. An area with 200 expected cases barely moves; an area with 0.4 expected cases moves almost all the way. The methods differ in what mean they pull toward.

Global empirical Bayes pulls toward the study-wide rate. It is the simplest, it needs no neighbour graph, and it removes small-number noise effectively, but it discards geography entirely — a sparse tract in a genuinely high-rate region is pulled toward the global mean, away from its neighbours.

Local empirical Bayes pulls toward the mean of each area’s neighbours, which respects geography at the cost of needing a weights matrix, built as in Spatial Weights Matrix Construction.

Fully Bayesian hierarchical models, of which BYM2 is the current standard, estimate the amount of smoothing from the data rather than assuming it, and separate the variation into a spatially structured component and an unstructured one. That separation is the reason to prefer them: it answers whether the residual variation is geographic at all.

Where the Extremes on a Raw Rate Map Come From Estimated rate plotted against the expected case count for each area. Raw rates form a funnel: areas with fewer than about five expected cases scatter from near zero to more than four times the mean, while areas with fifty or more expected cases sit tightly around it. After empirical Bayes smoothing the funnel closes, with the sparse areas pulled strongly toward the mean and the well-populated ones essentially unchanged. The extremes of the raw map were the small denominators, not the high risks. The funnel is the whole problem 0 2× 4× rate / mean raw rates smoothed 1 10 50 200 expected cases in the area (log scale) Areas above 50 expected cases barely move — smoothing only touches what was unreliable

Method Selection

Method Pulls toward Needs Best when
Global empirical Bayes study-wide rate nothing but counts quick stabilisation, no spatial structure expected
Local empirical Bayes neighbours’ mean a weights matrix spatial structure expected, speed matters
BYM2 a mixture of structured and unstructured weights matrix, MCMC or INLA the default for publication-quality maps
Poisson-gamma hierarchical study-wide rate, with uncertainty MCMC or conjugate update when full posteriors are needed without spatial terms

BYM2 is the recommended default because of its parameterisation. It expresses the total random-effect variance as one parameter and the proportion of it that is spatially structured as another, so the model reports directly how much of the residual variation is geographic. Earlier BYM formulations confounded the two, and their smoothing parameters were not interpretable across studies.

Spatial Data Prerequisites

  • Expected counts, not just population. Age-standardise before smoothing, or the model attributes an age-structure difference to spatial risk.
  • A connected weights matrix. An island produces an improper conditional distribution and either fails or silently degenerates; handle it as in Handling Island Polygons in Spatial Weights.
  • A stable geography for the whole period, harmonized if a boundary vintage changed.
  • Counts, not rates, as the response. The Poisson likelihood is what makes the shrinkage depend on the denominator; feeding it rates discards that.

Production Implementation

# Global and local empirical Bayes smoothing of small-area rates.
# Pinned: python 3.11, numpy==1.26.4, pandas==2.2.2, libpysal==4.12.1, esda==2.6.0
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.smooth")


def global_empirical_bayes(observed: np.ndarray, expected: np.ndarray) -> dict:
    """Poisson-gamma empirical Bayes shrinkage toward the study-wide relative risk.

    Returns the smoothed relative risk and the shrinkage weight per area, because
    the weight is what tells a reader how much of each estimate came from the
    area's own data rather than from its neighbours."""
    if np.any(expected <= 0):
        raise ValueError("expected counts must be positive; an area with zero "
                         "expected cases cannot be smoothed, only excluded")
    theta = observed.sum() / expected.sum()                 # study-wide RR
    # Method-of-moments estimate of the gamma prior's shape.
    s2 = np.average((observed / expected - theta) ** 2, weights=expected)
    v = theta / expected.mean()
    alpha = theta ** 2 / max(s2 - v, 1e-9)
    beta = alpha / theta

    shrunk = (observed + alpha) / (expected + beta)
    weight = expected / (expected + beta)                   # 1 = all own data, 0 = all prior
    log.info("global EB: theta=%.3f, prior alpha=%.2f beta=%.2f, "
             "median shrinkage weight %.2f", theta, alpha, beta, float(np.median(weight)))
    return {"rr": shrunk, "own_data_weight": weight, "theta": float(theta)}


def local_empirical_bayes(observed: np.ndarray, expected: np.ndarray, w) -> dict:
    """Shrink each area toward the mean of its neighbours rather than the global mean.

    `w` is a libpysal W. Islands must be resolved beforehand: an area with no
    neighbours has no local prior and would silently fall back to its raw rate."""
    if len(w.islands):
        raise ValueError(f"{len(w.islands)} island(s) present; resolve before smoothing")
    ids = list(w.id_order)
    idx = {k: i for i, k in enumerate(ids)}
    rr = np.empty(len(ids))
    wt = np.empty(len(ids))
    for k in ids:
        i = idx[k]
        nb = [idx[j] for j in w.neighbors[k]] + [i]         # include self
        o_nb, e_nb = observed[nb], expected[nb]
        theta_i = o_nb.sum() / e_nb.sum()
        s2 = np.average((o_nb / e_nb - theta_i) ** 2, weights=e_nb)
        v = theta_i / e_nb.mean()
        beta = theta_i / max(s2 - v, 1e-9)
        rr[i] = (observed[i] + theta_i * beta) / (expected[i] + beta)
        wt[i] = expected[i] / (expected[i] + beta)
    log.info("local EB: median shrinkage weight %.2f, min %.2f", float(np.median(wt)), float(wt.min()))
    return {"rr": rr, "own_data_weight": wt}


def shrinkage_report(expected: np.ndarray, weight: np.ndarray) -> pd.DataFrame:
    """How much of the map is essentially prior. This belongs in the publication."""
    band = pd.cut(expected, [0, 1, 5, 20, 100, np.inf],
                  labels=["<1", "1-5", "5-20", "20-100", "100+"])
    out = (pd.DataFrame({"band": band, "w": weight})
           .groupby("band", observed=True)["w"].agg(["size", "median"]))
    log.info("shrinkage by expected-count band:\n%s", out.to_string())
    return out

own_data_weight is the field most implementations discard and the one a reader most needs. An area whose smoothed rate is 90% prior is not an estimate of that area’s risk; it is an estimate of its neighbourhood’s risk displayed at that area’s location, and a map that does not distinguish the two invites exactly the over-reading it was meant to prevent.

Parameter Selection & Tuning

  • Choose the prior’s target deliberately. Global versus local is not a technical detail: global smoothing removes spatial structure, and using it before a Moran’s I test will suppress the very signal the test is looking for.
  • For BYM2, put a penalised-complexity prior on the total variance and on the spatial proportion, with the prior mass concentrated on “no effect”. The default is to shrink toward simplicity unless the data insists otherwise.
  • Do not tune the smoothing to make the map look right. The amount of shrinkage is determined by the data under the model; adjusting it to sharpen a feature is fitting the picture.
  • Set an exclusion floor. Areas with an expected count below roughly 0.5 contribute almost nothing and receive almost pure prior; consider aggregating them rather than mapping them.

Edge Cases & Failure Modes

Smoothing then testing. Running a cluster test on smoothed values is circular: the smoothing has already imposed spatial correlation, and the test will find it. Test on raw counts with an appropriate model, and smooth only for display.

Zero expected counts. An area with no population at risk cannot have a rate. Exclude it explicitly rather than letting a division produce infinity or a silent NaN.

Over-smoothing a true isolated excess. A genuine single-tract outbreak in a sparse area is exactly what shrinkage suppresses. This is the honest cost of the method, and it is why smoothed maps are a poor primary detector — pair them with the scan statistics in Space-Time Cluster Detection.

Mis-specified expected counts. If the expectation omits a covariate that varies spatially, the model attributes that variation to spatial risk and produces a confident, smooth, wrong map.

What Smoothing Removes and What It Keeps Two maps of the same county. The raw-rate map shows scattered extreme high and low tracts with no coherent pattern, all of them sparsely populated. The smoothed map shows a coherent elevated band across the centre-east and no isolated extremes. The scattered highs in the first map were single extra cases in tracts of a few hundred people; the band in the second was present in both and is what survives. Scattered extremes vanish; a coherent band survives raw rates smoothed relative risk 4 extremes, all in tracts under 600 people one coherent band, present in both Publish the shrinkage weight beside the map, or the second panel reads as more certain than it is

Compliance & Audit Controls

  • Publish the shrinkage weight, or an equivalent uncertainty measure, per area. A smoothed map without it presents modelled values with the visual authority of measured ones.
  • Record the expected-count model — the standard population, the age bands, the covariates — since the smoothed result depends on it at least as much as on the smoothing method.
  • State that the map is modelled in the legend, not the footnote.
  • Do not use smoothed values as a suppression input. Disclosure decisions are about the underlying counts, and a smoothed value can be safely above a threshold while the count behind it is one.
  • Pin the seed and sampler settings for fully Bayesian fits, and record convergence diagnostics with the output rather than only inspecting them.

Reporting a Smoothed Map Without Overclaiming

A smoothed map is a model output that looks like a measurement, and most of the trouble it causes downstream comes from that resemblance rather than from anything wrong in the fit. Three reporting habits close most of the gap.

Show the shrinkage. Every area’s estimate is a weighted blend of its own data and its neighbourhood’s, and the weight is already computed. Publishing it — as a second panel, as hatching on the main map, or simply as a column in the accompanying table — lets a reader see which parts of the map are supported by observation and which are largely interpolation. Without it, a tract whose value is nine parts prior is indistinguishable from one whose value is nine parts data, and readers will treat both as findings.

Show the interval, not only the point. A posterior mean is a summary of a distribution, and for sparse areas that distribution is wide. Where a single map is required, prefer the exceedance probability treatment described in the companion guide, because it folds the uncertainty into the quantity being coloured rather than discarding it. Where a table accompanies the map, give the credible interval for every area rather than only for the extremes.

Say what the map is for. A smoothed map answers “what is the underlying level of risk across this region”, and it answers it well. It does not answer “where should we investigate this week”, because the shrinkage that stabilises the level also suppresses the isolated excess an investigation would target. Stating the intended use in the caption prevents the map being repurposed as a detector by someone who was not in the room when it was made.

Be explicit about what changed since the last release. A smoothed series can move because incidence moved, because the population estimates were revised, because the expected-count model was refitted, or because the smoothing parameters changed. Only the first is epidemiological. Where a map is published on a schedule, carry a short change log alongside it recording which of the four applies, since a reader comparing this quarter to last has no other way to tell.

One further habit is worth adopting even though it costs a little effort: publish the raw counts and expected counts alongside the smoothed values. It seems redundant, and it is the single most effective defence against a smoothed figure being quoted years later as an observation. Anyone who needs to check the map can, and anyone who wants to fit a different model has what they need.

The size of the effect is worth quantifying before choosing a method, because it decides how much of the map smoothing will actually touch:

Share of Areas Dominated by Their Prior Share of areas whose smoothed estimate draws more than half its weight from the prior rather than from their own data, at four aggregation levels. At block-group level 71 percent of areas are prior-dominated, at tract level 44 percent, at ZCTA level 19 percent and at county level 4 percent. The coarser the unit, the less smoothing does, which is the same statement as saying smoothing exists to repair small denominators. How much of the map is mostly prior Block group 71% Census tract 44% ZCTA 19% County 4% At county level smoothing changes almost nothing — which is a reason to use it there and to be careful below

Implementation Checklist

FAQ

Does smoothing hide real outbreaks? It can, and that is its known cost. Shrinkage suppresses isolated excesses in sparse areas, which is exactly what a single-tract outbreak looks like. Use smoothed maps for describing level and scan statistics for detecting change.

Global or local empirical Bayes? Local, unless there is a specific reason to remove spatial structure. Global smoothing is appropriate when the goal is a stable ranking without geographic assumptions, and inappropriate before any spatial test.

Why BYM2 rather than the original BYM? Because its two parameters are separately interpretable: total random-effect standard deviation and the proportion of that variance which is spatially structured. The original parameterisation confounded them, so its estimates were not comparable across studies.

Can I smooth and then compute Moran’s I? No. Smoothing imposes spatial correlation by construction, so the test will confirm what the smoother assumed. Test on the raw counts using a model that accounts for the varying denominators.

Comparing Smoothed Maps Across Places and Times

Two smoothed maps are only comparable when the quantities behind them are. Three conditions have to hold, and each fails routinely in practice.

The expected-count model must be the same. A map standardised to the 2000 US standard population and one standardised to a local population are on different scales, and their relative risks cannot be placed on a shared legend. Where a comparison across jurisdictions is required, agree the standard first; where a comparison across years within one programme is required, hold it fixed even when a newer standard becomes available, and note the version.

The smoothing must be the same. Two counties fitted separately, each with its own variance estimated from its own data, will have shrunk by different amounts, so a tract at the same relative risk in each has a different amount of evidence behind it. Fitting jointly, with a shared hierarchy, solves this at the cost of assuming the two counties belong to one population — which is a substantive assumption and should be stated rather than assumed by convenience.

The geography must be the same vintage. This is the condition most often violated silently in a time series, and the remedy is harmonization before fitting rather than after.

Where those conditions cannot all be met, the honest presentation is separate maps with separate legends and an explicit note that cross-panel comparison is not supported. That is less satisfying than a single continuous surface across a region, and it is a great deal better than a surface whose apparent gradients are artefacts of three different standardisations meeting at a state line.