Spatial Epi Public Health GIS

CUSUM Aberration Detection by Small Area

A scan statistic asks whether any region has an excess right now. A cumulative-sum detector asks a different question, one area at a time: has this area’s level shifted, and if so, when did the shift begin? This guide, part of Space-Time Cluster Detection, covers running CUSUM per small area, choosing its two parameters, and handling the fact that a hundred detectors alarm a hundred times as often.

Problem Context & Constraints

CUSUM accumulates the signed distance between what was observed and what was expected, resetting whenever the accumulation goes negative. That accumulation is what gives it a property no single-period test has: it detects a sustained small shift that never produces a single striking day. A ten-percent rise that persists for six weeks is invisible to a threshold on daily counts and obvious to CUSUM by week three.

Its cost is that it is univariate and blind to geography. Run independently in each area, it has no notion that two adjacent areas alarming together is more meaningful than two distant ones, and it will never find a cluster that straddles a boundary without either area alarming alone. That is precisely complementary to a scan statistic, which sees the geography and struggles with slow shifts.

Two constraints govern its use in surveillance. It needs a baseline expectation per area, which small areas with a handful of cases per month cannot supply reliably. And it produces one detector per area, so a system covering 900 tracts at a nominal one-in-500-periods false-alarm rate expects roughly two false alarms per period across the map.

What CUSUM Sees That a Threshold Does Not Weekly counts in one area sit near their expected level, then rise by about ten percent from week nine onward. No single week exceeds a conventional threshold. The CUSUM statistic, plotted below, accumulates the small weekly excesses, crosses its decision interval at week twelve, and signals. A threshold on the weekly count alone would not have signalled at all. No single week is remarkable; the accumulation is weekly count against expectation expected CUSUM statistic decision interval h shift begins signal, 3 weeks later The three-week lag is the price of detecting a shift no single week reveals

Prerequisites

  • Counts by area and period with a stable baseline expectation per area — a seasonal model, a trailing mean, or an expected count from a population denominator
  • At least 24 baseline periods per area before the detector is trusted
  • python 3.11, numpy 1.26.4, pandas 2.2.2
  • A decision about which areas are too small to run at all

The decision interval has to be calibrated against each area’s own baseline, because the same value behaves very differently at different expected counts:

In-Control Run Length at h = 5 by Area Size Average run length before a false signal, at a fixed decision interval of five, for areas of different expected weekly counts. An area expecting two cases signals falsely every 118 periods on average. At five expected cases it is 291. At twenty it is 486. At eighty it is 812. A single decision interval therefore gives the sparse areas a false-alarm rate seven times that of the dense ones. One decision interval, seven-fold difference in false alarms 2 expected/week 118 periods 5 expected/week 291 20 expected/week 486 80 expected/week 812 Calibrate h per area, or accept that sparse areas dominate the alarm stream

Step-by-Step Solution

# Per-area CUSUM with a Poisson reference value and cross-area alarm control.
# 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("aberration.cusum")


def poisson_reference(mu0: float, shift: float = 1.5) -> float:
    """Reference value k for detecting a `shift`-fold increase in a Poisson mean.

    k is the midpoint on the log scale between the in-control mean and the
    out-of-control mean. Choosing `shift` is choosing what size of change the
    detector is tuned for; it is not a tuning knob to be swept for effect."""
    mu1 = shift * mu0
    return (mu1 - mu0) / np.log(mu1 / mu0)


def cusum_series(obs: np.ndarray, expected: np.ndarray, shift: float = 1.5,
                 h: float = 5.0) -> pd.DataFrame:
    """Upper one-sided CUSUM. Returns the statistic, the signal flag and the
    inferred start of the shift (the last reset before the signal)."""
    s = np.zeros(len(obs))
    start = np.zeros(len(obs), dtype=int)
    last_reset = 0
    for t in range(len(obs)):
        k = poisson_reference(max(expected[t], 1e-9), shift)
        prev = s[t - 1] if t else 0.0
        s[t] = max(0.0, prev + obs[t] - k)
        if s[t] == 0.0:
            last_reset = t
        start[t] = last_reset
    return pd.DataFrame({"cusum": s, "signal": s >= h, "shift_start": start})


def run_all_areas(counts: pd.DataFrame, min_baseline_mean: float = 2.0,
                  shift: float = 1.5, h: float = 5.0) -> pd.DataFrame:
    """CUSUM in every area, skipping areas too sparse to support a baseline."""
    out = []
    skipped = 0
    for area, g in counts.groupby("area_id"):
        g = g.sort_values("period")
        if g["expected"].mean() < min_baseline_mean:
            skipped += 1
            continue
        res = cusum_series(g["observed"].to_numpy(float), g["expected"].to_numpy(float),
                           shift, h)
        res["area_id"] = area
        res["period"] = g["period"].to_numpy()
        out.append(res)
    log.info("ran CUSUM in %d areas, skipped %d with mean expected < %.1f",
             len(out), skipped, min_baseline_mean)
    return pd.concat(out, ignore_index=True) if out else pd.DataFrame()


def contiguous_signal_groups(signals: pd.DataFrame, neighbours: dict) -> list[set]:
    """Group simultaneously signalling areas into contiguous sets.

    This is the cheap repair for CUSUM's geographic blindness: two adjacent areas
    signalling in the same period is far stronger evidence than two isolated ones,
    and grouping them turns a list of alarms into a candidate cluster."""
    active = set(signals.loc[signals["signal"], "area_id"])
    groups, seen = [], set()
    for a in active:
        if a in seen:
            continue
        stack, comp = [a], set()
        while stack:
            x = stack.pop()
            if x in comp:
                continue
            comp.add(x)
            seen.add(x)
            stack.extend(n for n in neighbours.get(x, []) if n in active and n not in comp)
        groups.append(comp)
    groups.sort(key=len, reverse=True)
    log.info("%d signalling areas formed %d contiguous groups (largest %d)",
             len(active), len(groups), len(groups[0]) if groups else 0)
    return groups

The last function is what makes per-area CUSUM usable as a surveillance system rather than as 900 unrelated detectors. A contiguous group of five signalling tracts is a candidate cluster worth investigating; five scattered singletons across a county are, at the usual thresholds, roughly what chance produces.

Validation & Edge Cases

1. Calibrate h on quiet data, per area size. The decision interval controls the in-control average run length, and the same h gives very different run lengths in an area expecting two cases a week and one expecting forty. Calibrate by simulation from each area’s own baseline:

INFO ran CUSUM in 704 areas, skipped 204 with mean expected < 2.0
INFO h=5.0 gives in-control ARL 486 periods (median across areas, IQR 291-812)
INFO 9 signalling areas formed 3 contiguous groups (largest 5)

2. Reset after an investigated signal. A CUSUM that signalled and was investigated must be reset, or it will signal again next period from the same accumulation. Failing to reset is the most common cause of a detector that alarms every week forever.

3. Watch the inferred shift start. The last reset before a signal estimates when the change began, and it is one of CUSUM’s most useful outputs for an investigation — often more useful than the signal itself. Report it.

4. Do not run CUSUM on areas that cannot support it. An area expecting 0.3 cases per week produces a CUSUM dominated by whether a single case arrived. The min_baseline_mean guard is not optional, and the skipped areas should be covered by the scan statistic instead.

5. Treat a step change in the baseline as a specification problem. If a laboratory comes online mid-series, every area it serves shifts, and CUSUM will faithfully signal for all of them. Annotate known changes and re-baseline rather than acting on the alarms.

6. Consider a two-sided detector where under-reporting matters. The upper CUSUM above detects increases. A lower CUSUM detects sustained decreases, which in surveillance usually means a reporting failure rather than falling incidence — a laboratory that stopped submitting, a feed that silently broke. Running both is nearly free and turns the detector into a data-quality monitor as well as an epidemiological one.

7. Match the period to the reporting cadence. Running a daily CUSUM on weekly-batched data produces six zero days and one spike per week, and the accumulator will signal on the batching pattern. Aggregate to the reporting period first, and where different streams report at different cadences, run separate detectors rather than a single mixed one.

Grouping simultaneous signals by contiguity is what converts a list of alarms into something an investigation can scope:

Nine Signalling Areas, Three Groups Nine areas signalling in the same period, arranged on a grid. Five of them form a contiguous block in the upper left, three form a smaller adjacent pair and single in the lower right, and one signals alone in the centre. The five-area group is a candidate cluster worth investigating; the isolated signal is what chance produces at this threshold across this many detectors. Contiguity turns nine alarms into one investigation group of 5 isolated group of 2 Investigate the group of five; the isolated signal is what 704 detectors produce by chance

Compliance Notes

  • Log shift, h, the baseline model and the reset history per area. A signal is only reproducible if the accumulator’s state history is.
  • Record skipped areas explicitly. A map showing signals from 704 areas, with 204 silently absent, understates coverage in exactly the sparse places most likely to need it.
  • Keep the alarm history, since the false-alarm rate is a property of the deployed system that can only be measured from the full series.
  • Apply disclosure controls to signalling-area lists before circulation, as with any small-area output.