Spatial Epi Public Health GIS

Diagnosing Over-Smoothing in a Disease Map

Smoothing is meant to remove small-number noise, and a smoothed map always looks better than a raw one. That is exactly the problem: a map that has smoothed away a real localised excess looks just as good as one that has not. This guide, part of Bayesian Disease Mapping & Rate Smoothing, gives three diagnostics that distinguish the two.

Problem Context & Constraints

Over-smoothing has no visual signature. Both an appropriately smoothed map and an over-smoothed one show fewer extremes than the raw map, smoother gradients, and no isolated spikes. The difference is only visible against something external: the residuals, a held-out prediction, or a signal you planted yourself.

Three checks cover the ground, and they fail in different ways.

Residual structure. If the smoother is right, the standardised residuals should be spatially unstructured. Leftover spatial autocorrelation in the residuals means the model has not captured the geography; leftover anti-correlation — adjacent residuals of opposite sign — is the specific signature of over-smoothing, because a suppressed peak leaves a positive residual surrounded by negatives.

Predictive score. Leave-one-out predictive density scores the model on data it did not see. An over-smoothed model predicts the held-out area’s count too close to its neighbours’ level, and the score falls.

Planted-signal recovery. Inject a known excess into a known area, refit, and measure how much of it survives. This is the most direct diagnostic and the one most rarely run, and it is what turns “the model may suppress outbreaks” into a number.

How Much of a Planted Excess Survives A twofold relative-risk excess is injected into a single tract and the model refitted. Under light smoothing 84 percent of the excess is recovered in the posterior mean. Under the fitted BYM2 configuration 61 percent is recovered. Under the legacy gamma prior only 19 percent survives, so an outbreak of that size in that tract would be invisible on the published map. The same test run in a densely populated tract recovers 93 percent under all three. A 2× excess planted in a sparse tract, and what came back Light smoothing 84% Fitted BYM2 61% Legacy gamma prior 19% Same test, dense tract 93% 0% 50% 100% share of the planted excess recovered in the posterior mean The last row is why the test must be run in a sparse tract, where the risk is

Prerequisites

  • A fitted smoothing model with posterior draws or, for empirical Bayes, the smoothed estimates and shrinkage weights
  • The raw observed and expected counts
  • A weights matrix for the residual autocorrelation test
  • python 3.11, numpy 1.26.4, pandas 2.2.2, libpysal 4.12.1, esda 2.6.0

The residual signature is directional, and knowing which direction means what is the whole of the first diagnostic:

Reading the Sign of Residual Autocorrelation Three residual patterns and their interpretations. Positive residual autocorrelation means adjacent residuals share a sign, indicating spatial structure the model has not captured and therefore under-smoothing or a missing covariate. Near-zero residual autocorrelation is the target. Negative residual autocorrelation means adjacent residuals alternate in sign, the signature of a suppressed peak surrounded by compensating deficits, which is over-smoothing. The sign says which way the model is wrong I > 0 structure not captured I ≈ 0 the target I < 0 a peak was suppressed Most implementations test only the left case, which is the one that is not over-smoothing

Step-by-Step Solution

# Three over-smoothing diagnostics: residual structure, LOO score, planted signal.
# Pinned: numpy==1.26.4, pandas==2.2.2, libpysal==4.12.1, esda==2.6.0, scipy==1.13.1
import logging
import numpy as np
import pandas as pd
from esda.moran import Moran
from scipy.stats import poisson

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


def residual_structure(observed, expected, fitted_rr, w, seed: int = 42) -> dict:
    """Moran's I on standardised residuals.

    Positive I means unmodelled spatial structure remains. Negative I is the
    over-smoothing signature: a suppressed peak leaves a positive residual ringed
    by negatives, which reads as spatial anti-correlation."""
    mu = expected * fitted_rr
    resid = (observed - mu) / np.sqrt(np.maximum(mu, 1e-9))
    mi = Moran(resid, w, permutations=999, seed=seed)
    verdict = ("under-smoothed / unmodelled structure" if mi.I > 0 and mi.p_sim < 0.05
               else "over-smoothed" if mi.I < 0 and mi.p_sim < 0.05
               else "no residual structure detected")
    log.info("residual Moran's I = %.4f (p=%.4f) -> %s", mi.I, mi.p_sim, verdict)
    return {"I": float(mi.I), "p": float(mi.p_sim), "verdict": verdict}


def loo_predictive_score(observed, expected, draws) -> float:
    """Mean log pointwise predictive density, computed from posterior draws.

    Higher is better. Comparing this across smoothing configurations is the
    cheapest defensible way to say one is over-smoothed relative to another."""
    lppd = []
    for i in range(len(observed)):
        mu = expected[i] * draws[:, i]
        lppd.append(np.log(np.mean(poisson.pmf(observed[i], mu))))
    score = float(np.mean(lppd))
    log.info("mean log pointwise predictive density: %.4f", score)
    return score


def planted_signal_recovery(observed, expected, target_idx: int, factor: float,
                            refit, seed: int = 42) -> dict:
    """Inject a `factor`-fold excess into one area and measure what survives.

    `refit` takes (observed, expected) and returns fitted relative risks. Run this
    in a SPARSE area: recovery in a dense area is always high and tells you nothing
    about the case you are worried about."""
    y = observed.copy().astype(float)
    added = expected[target_idx] * (factor - 1.0)
    y[target_idx] += added
    rr = refit(y, expected)
    baseline = refit(observed, expected)[target_idx]
    recovered = (rr[target_idx] - baseline) / (factor - 1.0)
    log.info("planted %.1fx in area %d (expected %.1f): recovered %.0f%% of the excess",
             factor, target_idx, expected[target_idx], 100 * recovered)
    return {"area": int(target_idx), "factor": factor,
            "expected": float(expected[target_idx]), "recovered_fraction": float(recovered)}


def recovery_profile(observed, expected, refit, factor=2.0, n=12, seed=42) -> pd.DataFrame:
    """Recovery across the expected-count range, which is the publishable summary."""
    rng = np.random.default_rng(seed)
    order = np.argsort(expected)
    picks = order[np.linspace(0, len(order) - 1, n).astype(int)]
    rows = [planted_signal_recovery(observed, expected, int(i), factor, refit) for i in picks]
    out = pd.DataFrame(rows).sort_values("expected")
    log.info("recovery ranges from %.0f%% (expected %.1f) to %.0f%% (expected %.1f)",
             100 * out["recovered_fraction"].iloc[0], out["expected"].iloc[0],
             100 * out["recovered_fraction"].iloc[-1], out["expected"].iloc[-1])
    return out

The recovery profile is the artefact worth publishing. It converts a vague caveat — “smoothing may suppress local excesses” — into a curve showing exactly how large an excess has to be, in an area of a given size, before the published map would show it.

Validation & Edge Cases

1. Run the residual test both ways. Most implementations check only for positive residual autocorrelation. The negative direction is the one that indicates over-smoothing and it is routinely ignored:

INFO residual Moran's I = -0.0871 (p=0.0140) -> over-smoothed
INFO mean log pointwise predictive density: -2.4118   [fitted BYM2]
INFO mean log pointwise predictive density: -2.6902   [legacy gamma prior]
INFO planted 2.0x in area 41 (expected 3.2): recovered 19% of the excess
INFO recovery ranges from 14% (expected 1.8) to 96% (expected 412.0)

2. Plant the signal in a sparse area, and in several. Recovery in a well-populated area is near-perfect under any configuration. A single test in a randomly chosen area will usually be reassuring and uninformative.

3. Compare predictive scores only between models on identical data. The score is not interpretable in absolute terms; a difference of 0.28 in mean log predictive density between two configurations on the same counts is meaningful, the value itself is not.

4. Distinguish over-smoothing from a missing covariate. Positive residual autocorrelation can mean either. Adding a plausible spatially varying covariate and re-testing separates them: if the residual structure disappears, the model was under-specified rather than under-smoothed.

5. Report the recovery profile with the map, not only in the methods. A reader deciding whether to act on the absence of a hotspot needs to know that an excess below a certain size in an area below a certain population would not have appeared.

6. Run the diagnostics on the published configuration, not a convenient one. It is easy to validate a fast approximation and then publish a fuller model, or to validate on a subset and publish on the whole state. Both leave the actual published map untested, and the difference between configurations is exactly what the diagnostics exist to measure.

7. Repeat the recovery test annually. The population denominators shift, the case counts change scale, and a configuration whose recovery profile was acceptable three years ago may no longer be. Because the test is a simulation on top of an existing fit, re-running it costs a fraction of the original analysis, and the resulting profile is the natural place to notice that a surveillance programme’s effective sensitivity has drifted downward as its areas have grown or its counts have fallen. Where the profile has moved materially, that is worth reporting in its own right rather than silently re-baselining.

The recovery profile is more useful as a curve than as a single number, because it tells a reader exactly which excesses the map could not have shown:

Detection Floor by Area Size The fold increase required for an excess to be visible on the published smoothed map, plotted against the area's expected count. An area expecting two cases needs a fivefold excess before the map would show it. At ten expected cases the threshold falls to about 2.2, at fifty to 1.4 and at two hundred to 1.15. The curve is the map's sensitivity, and it is steep exactly where surveillance most wants sensitivity. How large an excess has to be before the map shows it 1× 3× 5× 2 expected — needs 5× 10 expected — needs 2.2× 200 expected — needs 1.15× 2 30 400 expected cases in the area (log scale)

Compliance Notes

  • Publish the recovery profile alongside any smoothed surveillance map used for resource allocation. It is the map’s sensitivity, and a map without a stated sensitivity cannot support a decision not to act.
  • Record the diagnostic results in the run metadata, including the residual Moran’s I and its sign, so a later reviewer can see they were run.
  • Do not tune the smoothing parameter to maximise recovery. Recovery and noise suppression trade off directly, and optimising one produces a map that is unstable in the other direction.
  • Treat a negative residual Moran’s I as a blocking finding, not a footnote: it means the published map is suppressing structure that the data contains.