The real spread of COVID-19, from the Wuhan cluster on 31 December 2019 to 31 March 2020. Black dots show where people live; each red dot is about 100 confirmed cases. Planes fly the first documented routes, and on 11 March the WHO declares a pandemic. At the end of March the scan settles on a window around the Alps. Data: Johns Hopkins CSSE; Italy, Protezione Civile. Earth: NASA Blue Marble.
Say we want to know where COVID-19 cases were unusually dense, and by how much. The black dots show where people live: that’s the population, b. The red dots are confirmed cases: we call them m. The scan slides a window over the map, at every size and position, and gives every window a score, Kulldorff’s log-likelihood ratio:
Here m(Y) and b(Y) are the shares of all cases and of all people inside the window Y, and the expected share is what the same rate everywhere would give.
In plain English, for the window it found:
just over 1% of all people live here
so we expect just over 1% of the cases
17% of the cases are here
13× more than expected → high score
The scan keeps the window with the highest score.
Scan statistics are used to find cancer clusters, disease outbreaks such as COVID-19, emerging public-health threats and street crime: hotspots of cases, where they are and how unusual. The common tools work on points. SaTScan is the most widely used; for regions, it uses one point per region, typically the centroid. pyscan is a fast, open-source scanner, Python on a C++ backend, by my co-author Michael Matheny, and it is what we use.
The gap
Real data come as counts per region
The problem. Public-health and crime maps almost never come with addresses — you get one total for a whole county or zip code, and some counties are the size of a small state. The fastest algorithms for spotting a hotspot need actual points on the map, so the standard shortcut squashes each county down to a single dot at its center. That dot erases the county’s real shape and size — and a disease outbreak that a health department finds late, or not at all, is measured in lives, not in decimals of statistical power.
Published · one rate per countyCentroid · one point per county
Valley Fever is a fungal infection. Its fungus lives in the soil of California's San Joaquin Valley (dashed), so we know where the hotspot should be. The state publishes one rate per county, with no addresses, for privacy. With one point per county, at its centroid, the scan's window misses most of the valley. Data: CDPH, 2014–2018.Why it misses: one point stands for the whole region. Kern County, the big Valley Fever county: a window over 46% of it sees zero cases with the centroid, and 48% with 50 points. A ZIP code in Henderson, Nevada: the centroid is inside, but near a thin neck; a window over 68% of the ZIP code sees zero cases, and 68% with points. Henrico County, Virginia, wraps around the City of Richmond, so its centroid is not even in Henrico; a window over 58% of the county sees zero cases, and 58% with points. A region is not a point. The grey box is one scan window.
Existing methods
Existing ways to scan regions are lossy, slow, or capped
The centroid is fast, but each region becomes one dot, and a window sees all of it or nothing. Buchin and colleagues compute exactly how much each window overlaps each region: accurate, but slow, and the window size must be given in advance. FlexScan searches connected groups of whole regions, but it is slow, and its code caps a cluster at 15 regions. What is missing is a fast and accurate way to run point-based scanners on region data.
The method
Replace each region with 20–50 random points
If your data are counts per region, don't shrink each region to a single point. Scatter 20–50 random points inside it, and share the counts among them. Same data, same scan, far more statistical power — and with pyscan, almost no extra time.
For each region, draw k points uniformly at random inside it, with k between 20 and 50. Each point carries 1/k of the region’s cases and people (Geom k): with k = 20, a county with 120 cases and 10,000 people becomes 20 points with 6 cases and 500 people each. The region keeps its spatial extent, and the scan algorithm is untouched: it is a preprocessing step.
120 cases 10,000 people6 cases each 500 people each
Before · 75 centroidsAfter · 20 points per county
How I solved it. Our fix is almost embarrassingly simple: instead of one dot at the center, scatter 20–50 points across the county’s actual shape and share its case count among them. The map gets its geography back, and the fast algorithms never notice the difference. The proof it works: given nothing but county totals of Valley Fever cases, the method redraws California’s San Joaquin Valley endemic region — a region fixed by soil ecology — far more faithfully than the centroid shortcut ever could, while scanning the entire continental-U.S. map — 3,711 region polygons spanning 3,108 counties — in 0.33 seconds, roughly 3,000× faster than FlexScan’s 1,109 s (scan step only; FlexScan additionally runs a significance test the other methods do not).
Speed
At national scale, the whole US scans in 0.33 s
The whole pipeline on the US map: the county borders, 50 points per county, a planted anomaly in red, and the scan trying windows of every size and position until it locks onto the anomaly. The moving windows are an illustration; the time is real. Satellite image: NASA Blue Marble.
Fifty times more points does not mean a much slower scan. pyscan’s time is set by a data-adaptive grid, not by the number of points:
Scan time, one discovery
1 point / region
50 points / region (ours)
Buchin et al.
FlexScan
USA (3,711 regions)
0.29 s
0.33 s
114 s
1,109 s
Arkansas (75 counties)
0.21 s
0.28 s
25.3 s
8.00 s
Scan step only, one discovery, same machine. The USA map is the county × congressional-district partition (3,711 regions). FlexScan's time includes a Monte Carlo significance test the others do not run; Buchin et al. sweeps nine window sizes.
Results
A synthetic test: plant a hotspot, see which input finds it
Centroid · found window PJD 0.5450 points · found window PJD 0.08
Arkansas, one typical trial. Inside the planted region (green), each point is a case with probability p = 0.5; outside, with q = 0.2. The difference p − q is the signal strength. On centroids the scan finds a window too big and off-centre; with 50 points it lands on the planted region. The Point Jaccard distance (PJD) measures how much the found and planted windows disagree: 0 is perfect, and below 0.2 is a good match.More points per region, better detection. Twenty trials at each signal strength on Arkansas; lower is better, and below the grey line is a good match. The centroid never gets there, and neither does one random point per region. With 5 points per region a good match comes at a signal of 0.45, with 10 at 0.30, and with 50 at just 0.15.Sampling wins on all six maps, from 29 to 3,108 regions: on each, 50 points per region crosses the good-match line first.
In Georgia, a planted rectangle grows from 1.7% of the state to 60.8%. At 4.1%, 50 points per region find it and the centroid does not; the centroid succeeds only at 60.8%, the largest size tested. Sampling finds planted regions 15× smaller than centroids can.Against the region methods on the same Arkansas test. FlexScan's own code caps a cluster at 15 regions; the planted region spans more than 15 counties, so it plateaus near 0.35. Buchin and colleagues' code takes one window size per run and cannot find the size; run with the authors' protocol, nine sizes around the true size, which is its best case, it gets just below the line, but 90× slower than ours. Ours has no cap and no size input, and goes well below the line.The finest maps need the fewest points. At a weak signal (0.15), the first k that reaches a good match is 3 for the US, 7 for New York City, 20 for California, 50 for Arkansas and Georgia, and 75 for Utah. A convergence analysis in the paper explains why: as a map is divided into more regions, fewer sample points per region are needed, and the error roughly follows 1/√k.
Real data
Valley Fever: within 0.036 of the best rectangle
50 points per county
The ground truth is the eight San Joaquin Valley counties, fixed by the soil, not by the counts. With 50 points per county, the scan's window sits on the valley: its Point Jaccard distance to the truth drops from 0.717 with the centroid to 0.448, within 0.036 of 0.412, the best any axis-aligned rectangle can do, because the valley runs diagonally.
Paper
Sampling for Region-Aggregated Spatial Scan Statistics
Accepted as a full paper at ACM SIGSPATIAL 2026 — 58 of 248 submissions. Joint work with Drew McClelland (Taptap Send), Michael Matheny (Meta), and Jeff M. Phillips.
We recommend this sampling-based conversion as the default way to apply point-based spatial scan statistics to region-aggregated data. The companion repository contains the Python experiments and figure-rendering scripts that reproduce every figure and runtime table in the paper, built on pyscan’s C++ backend.
