Spatial Epi Public Health GIS

Comparing Scan Statistics and Local Moran’s I

This guide, part of Cluster Detection Power & Benchmarking, puts the two most widely used local cluster methods against each other on identical synthetic outbreaks so the choice between them rests on evidence rather than habit.

Problem Context & Constraints

Health departments usually run one of two local methods. The spatial scan statistic searches over many candidate windows — typically circles of increasing radius around each centroid — and reports the window with the highest likelihood ratio, tested by Monte Carlo replication of the whole search. Local Moran’s I computes, for every area, how strongly its value co-varies with its neighbours’ values, tested by conditional permutation, and labels each area high-high, low-low, high-low or low-high.

They are often presented as interchangeable, and they are not. The scan statistic asks is there a region where risk is elevated, and where is the most likely one? It controls the family-wise error across all windows by construction, and it returns a small number of clusters with a single p-value each. Local Moran’s I asks which areas resemble their neighbours unusually strongly? It returns a label for every area, each with its own p-value, and without correction for many tests it flags areas by chance in proportion to the number of areas.

Comparing them fairly takes three commitments. Use the same outbreaks and the same null data for both, so differences reflect the methods rather than the inputs. Score on more than detection — location accuracy and false alarm rate matter as much. And feed each method the input it expects: the scan statistic works on counts and population, while Local Moran’s I works on a rate, and it should be a smoothed rate, since raw rates in small areas swamp the statistic with noise.

Two Methods, One Outbreak The same planted outbreak over a grid of areas, analysed two ways. The scan statistic returns a single circular window covering the outbreak plus some surrounding areas, with one p-value for the whole window. Local Moran's I returns a label for each area: most outbreak areas are labelled high-high, one edge area inside the outbreak is missed, and two areas far from the outbreak are flagged by chance because each area is tested separately. One window versus a label for every area Scan statistic one window p = 0.004 covers the outbreak and some neighbours Local Moran's I HH HH HH miss HH LH per-area labels, chance flags elsewhere neither output is wrong; they answer different questions

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4
  • libpysal 4.10, esda 2.5.1 for Local Moran’s I
  • A scan implementation: SaTScan 10.1 driven by parameter files, or a Python Poisson circular scan for fast benchmarking loops
  • A scenario library of planted outbreaks, as built in Simulating Outbreaks to Measure Detection Power

Step-by-Step Solution

Wrap each method behind the same interface — take counts and population, return the set of flagged areas and a signal decision — then score both against the truth for every replicate.

# Head-to-head benchmark: Poisson circular scan vs Local Moran's I.
# Pinned: numpy==1.26.4, pandas==2.2.2, geopandas==0.14.4, libpysal==4.10, esda==2.5.1
import logging

import numpy as np
import pandas as pd
from libpysal.weights import Queen
from esda.moran import Moran_Local

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

ALPHA = 0.05


def eb_rate(cases: np.ndarray, pop: np.ndarray) -> np.ndarray:
    """Global empirical Bayes smoothed rate, the input Local Moran's I needs."""
    r = cases / pop
    m = cases.sum() / pop.sum()
    v = max((pop * (r - m) ** 2).sum() / pop.sum() - m / pop.mean(), 0.0)
    shrink = v / (v + m / pop)
    return m + shrink * (r - m)


def lisa_flags(cases, pop, w, perms: int = 999, seed: int = 1) -> dict:
    """High-high areas from Local Moran's I, optionally FDR-corrected."""
    lm = Moran_Local(eb_rate(cases, pop), w, permutations=perms, seed=seed)
    p = lm.p_sim
    order = np.argsort(p)
    thresh = ALPHA * np.arange(1, len(p) + 1) / len(p)
    passed = np.zeros(len(p), dtype=bool)
    k = np.where(p[order] <= thresh)[0]
    if k.size:
        passed[order[: k.max() + 1]] = True
    hh = (lm.q == 1) & passed
    return {"flagged": set(np.flatnonzero(hh)), "signal": bool(hh.any())}


def scan_flags(cases, pop, centroids, max_frac: float = 0.25,
               reps: int = 999, seed: int = 1) -> dict:
    """Poisson circular scan: most likely cluster and its Monte Carlo p-value."""
    rng = np.random.default_rng(seed)
    d = np.linalg.norm(centroids[:, None, :] - centroids[None, :, :], axis=2)
    order = np.argsort(d, axis=1)
    total_c, total_p = cases.sum(), pop.sum()

    def best(c):
        top, window = 0.0, None
        for i in range(len(c)):
            cc = np.cumsum(c[order[i]])
            pp = np.cumsum(pop[order[i]])
            ok = pp <= max_frac * total_p
            e = total_c * pp / total_p
            with np.errstate(divide="ignore", invalid="ignore"):
                llr = np.where((cc > e) & ok,
                               cc * np.log(cc / e) + (total_c - cc) * np.log((total_c - cc) / (total_c - e)),
                               0.0)
            j = int(np.nanargmax(llr))
            if llr[j] > top:
                top, window = float(llr[j]), order[i][: j + 1]
        return top, window

    obs, window = best(cases)
    null = [best(rng.multinomial(total_c, pop / total_p))[0] for _ in range(reps)]
    p = (1 + sum(n >= obs for n in null)) / (reps + 1)
    flagged = set(window.tolist()) if (window is not None and p <= ALPHA) else set()
    return {"flagged": flagged, "signal": p <= ALPHA, "p": p}


def score(result: dict, truth: set) -> dict:
    """Detection, location precision and recall against the planted areas."""
    f = result["flagged"]
    tp = len(f & truth)
    return {
        "signal": result["signal"],
        "hit": tp > 0,
        "precision": tp / len(f) if f else np.nan,
        "recall": tp / len(truth) if truth else np.nan,
        "false_areas": len(f - truth),
    }


def compare(areas, scenarios, n_reps: int = 200) -> pd.DataFrame:
    """Run both methods on identical simulated data and collect scores."""
    w = Queen.from_dataframe(areas, use_index=False)
    w.transform = "r"
    xy = np.column_stack([areas.geometry.centroid.x, areas.geometry.centroid.y])
    pop = areas["population"].to_numpy(dtype=float)
    rows = []
    for sc in scenarios:
        rng = np.random.default_rng(sc["seed"])
        for r in range(n_reps):
            lam = pop * sc["rr"]
            cases = rng.multinomial(sc["total"], lam / lam.sum())
            for name, res in (("scan", scan_flags(cases, pop, xy, seed=r)),
                              ("lisa", lisa_flags(cases, pop, w, seed=r))):
                rows.append({"scenario": sc["name"], "method": name, "rep": r,
                             **score(res, sc["truth"])})
    out = pd.DataFrame(rows)
    log.info("scored %d method-replicates across %d scenarios", len(out), len(scenarios))
    return out

The null scenarios — where rr is all ones and truth is empty — are essential. They give each method’s false alarm rate, and a method that signals on 30 percent of null replicates has not earned its power on the outbreak scenarios.

The scan’s Monte Carlo loop dominates run time. For a benchmark of several hundred replicates across a dozen scenarios, run it with fewer replications (199) during development and the full 999 only for the reported results; the ranking of methods rarely changes, but the p-value resolution does.

A Typical Head-to-Head Scorecard A table-like scorecard with rows for compact outbreak power, corridor outbreak power, location precision, and false alarm rate on null data, and columns for the scan statistic and Local Moran's I with false discovery rate correction. The scan statistic is stronger on compact outbreaks and holds its false alarm rate at the nominal five percent, but its location precision is lower because circles include extra areas. Local Moran's I is stronger on corridor outbreaks and more precise in location, but without correction its false alarm rate would be far above nominal. Each method wins somewhere Scan Local Moran's I + FDR Power, compact outbreak 0.86 0.71 Power, corridor outbreak 0.48 0.69 Location precision 0.58 0.81 False alarms on null data 0.05 0.07 illustrative results; bold marks the stronger method on each row

Validation & Edge Cases

Unfair inputs. Running Local Moran’s I on raw rates while the scan uses counts and population hands the scan an advantage it did not earn: small-population areas with one or two cases dominate raw-rate LISA. Smooth the rate first, as the code does, and state it.

Different definitions of detection. The scan signals once per map; Local Moran’s I signals if any area is high-high and significant. Without FDR correction the second definition fires on almost every null map of a few hundred areas, which makes its power meaningless. Apply the correction — see Applying FDR Correction to Local Moran’s I Results — or compare at matched false alarm rates.

Weights choice. Local Moran’s I results depend on the neighbour definition; queen contiguity, k-nearest neighbours and distance bands give different power, especially on irregular geographies. Benchmark at least two weight specifications and report the spread.

Scan window limits. A maximum window of 50 percent of population, the SaTScan default, lets the scan return huge windows that swallow the outbreak and half the state. For surveillance, 10 to 25 percent is usually more useful, and the limit must be recorded alongside the result.

Location scoring for tiny outbreaks. When the planted outbreak is a single area, precision and recall are coarse — either 0 or 1 for recall. Score location by distance from the flagged centroid to the true centre as well, so near misses count differently from far misses.

Replicate agreement. The same seed feeding both methods is what makes the comparison paired. Analyse the results as paired: count replicates where one method detects and the other does not, and test that difference directly rather than comparing two independent power estimates with overlapping intervals.

INFO scan_vs_lisa scored 4800 method-replicates across 12 scenarios
INFO scan_vs_lisa paired discordance, corridor scenarios: lisa-only 211, scan-only 58, McNemar p < 0.001
Which Method for Which Question Three boxes connected to a decision. If the question is whether any cluster exists in the jurisdiction this period, with strict false alarm control, use the scan statistic. If the question is which specific areas behave like their neighbours, for mapping and targeting, use Local Moran's I with false discovery rate correction. If both matter, run both and treat agreement as the strongest evidence and disagreement as a prompt to look closer. Pick by the question, not by habit Is there a cluster? strict alarm control compact outbreaks scan statistic Which areas? mapping and targeting irregular shapes Local Moran's I + FDR Both matter run both methods agreement = strongest triangulate disagreement between methods is a prompt to look closer, not a tie to break

Compliance Notes

Record the version of each method’s implementation, the scan’s window shape and maximum size, the weights specification and rate smoothing used for Local Moran’s I, the multiple-testing correction, the replication counts, and the paired discordance results.

When a program adopts one method on the strength of this benchmark, keep the benchmark report with the surveillance protocol. The choice of method is a choice about which outbreaks the program is most likely to miss, and that is a decision reviewers and partner agencies are entitled to see.