Spatial Epi Public Health GIS

Diagnosing the Modifiable Areal Unit Problem in Rate Maps

The same cases, aggregated to different units, produce different maps, different correlations and different clusters. That is not a defect to be eliminated — there is no correct unit — but it is a property that has to be measured rather than assumed away. This guide, part of Areal Interpolation & Boundary Harmonization, sets out a practical diagnostic and a reporting standard for findings that depend on it.

Problem Context & Constraints

The problem has two separable parts. The scale effect is what happens when the units get bigger: correlations between variables generally strengthen, variance falls, and extremes disappear. The zoning effect is what happens when units of the same size are drawn differently: the same cases fall into different groupings and the statistics move again, sometimes as much as under a scale change.

Both matter for surveillance because the unit is almost never chosen for analytic reasons. Tracts exist for census enumeration, ZIPs for mail delivery, counties for administration, and none of them was drawn around disease transmission. A cluster detected on tracts is a statement about tract-shaped excesses, and a reader will hear it as a statement about places.

The constraint that makes the diagnostic tractable is that, given point-level data, alternative aggregations can be constructed cheaply. Where only aggregated data exists, the diagnostic is much weaker and the honest response is to report at the received scale rather than to interpolate to a finer one and pretend the problem has gone away.

One Point Pattern, Three Aggregations, Three Findings The same scatter of case points aggregated onto a fine grid, a coarse grid, and a coarse grid offset by half a cell. The fine grid shows a maximum cell rate of 214 per 100,000 in the upper right. The coarse grid shows a maximum of 96 in the centre right. The offset coarse grid, with cells of identical size, shows a maximum of 143 in a different cell again. Scale changes the magnitude; zoning changes the location. Scale changes the magnitude; zoning changes the location fine grid coarse grid coarse, offset half a cell max 214 / 100k max 96 / 100k max 143 / 100k upper right upper right lower left Panels two and three have identical cell sizes and disagree about where the excess is

Prerequisites

  • Point-level case locations, or the finest aggregation available
  • A population denominator at the same or finer resolution
  • python 3.11, geopandas 1.0.1, numpy 1.26.4, libpysal 4.12.1
  • A metric CRS, since the diagnostic constructs grids of stated size

Step-by-Step Solution

The diagnostic constructs a family of aggregations and reports how far the finding moves across them.

# Measure scale and zoning sensitivity of an areal rate statistic.
# Pinned: geopandas==1.0.1, numpy==1.26.4, pandas==2.2.2, shapely==2.0.6
import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import box

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


def make_grid(bounds, cell: float, dx: float = 0.0, dy: float = 0.0, crs=None):
    """Regular grid of `cell`-sized squares, offset by (dx, dy).

    The offset is what separates the zoning effect from the scale effect: two
    grids with the same `cell` and different offsets differ only in zoning."""
    minx, miny, maxx, maxy = bounds
    xs = np.arange(minx - cell + dx, maxx + cell, cell)
    ys = np.arange(miny - cell + dy, maxy + cell, cell)
    cells = [box(x, y, x + cell, y + cell) for x in xs for y in ys]
    return gpd.GeoDataFrame({"cell_id": range(len(cells))}, geometry=cells, crs=crs)


def rate_on_grid(points: gpd.GeoDataFrame, pop: gpd.GeoDataFrame,
                 grid: gpd.GeoDataFrame, pop_col: str = "pop") -> pd.Series:
    """Cases per 100,000 in each grid cell, cells with no population dropped."""
    p = gpd.sjoin(points, grid, how="inner", predicate="within")
    cases = p.groupby("cell_id").size()
    d = pop.copy()
    d["geometry"] = d.geometry.representative_point()
    d = gpd.sjoin(d, grid, how="inner", predicate="within")
    denom = d.groupby("cell_id")[pop_col].sum()
    joined = pd.concat([cases.rename("cases"), denom.rename("pop")], axis=1).fillna(0)
    joined = joined[joined["pop"] > 0]
    return 1e5 * joined["cases"] / joined["pop"]


def maup_sensitivity(points, pop, cells=(1000, 2000, 4000, 8000),
                     offsets=4, pop_col="pop", seed=42) -> pd.DataFrame:
    """Max rate and its location across a family of scales and zonings."""
    rng = np.random.default_rng(seed)
    bounds, crs = points.total_bounds, points.crs
    rows = []
    for cell in cells:
        for k in range(offsets):
            dx, dy = (0.0, 0.0) if k == 0 else tuple(rng.uniform(0, cell, 2))
            g = make_grid(bounds, cell, dx, dy, crs=crs)
            r = rate_on_grid(points, pop, g, pop_col)
            if r.empty:
                continue
            rows.append({"cell_m": cell, "offset": k, "n_cells": len(r),
                         "max_rate": float(r.max()), "cv": float(r.std() / r.mean()),
                         "argmax_cell": int(r.idxmax())})
    out = pd.DataFrame(rows)
    for cell, g in out.groupby("cell_m"):
        spread = g["max_rate"].max() - g["max_rate"].min()
        log.info("cell %5d m: max rate %.0f-%.0f (zoning spread %.0f), CV %.2f-%.2f",
                 cell, g["max_rate"].min(), g["max_rate"].max(), spread,
                 g["cv"].min(), g["cv"].max())
    return out

Read the output two ways. Down the cell_m column is the scale effect: maximum rate and coefficient of variation both fall as cells grow, and the rate at which they fall says how concentrated the underlying pattern is. Within a cell_m group is the zoning effect: the spread of max_rate across offsets of identical size is the amount of the headline figure that is an artefact of where the lines happen to fall.

Validation & Edge Cases

1. Report the zoning spread as an interval, not a footnote. If a tract-level analysis reports a maximum rate of 214 per 100,000 and the zoning diagnostic shows a spread of 96 to 214 across equivalent aggregations, the honest headline is the interval:

INFO cell  1000 m: max rate 186-241 (zoning spread 55), CV 1.31-1.48
INFO cell  2000 m: max rate 118-163 (zoning spread 45), CV 0.94-1.09
INFO cell  4000 m: max rate  81-112 (zoning spread 31), CV 0.66-0.78
INFO cell  8000 m: max rate  54- 71 (zoning spread 17), CV 0.41-0.49

2. Check whether the location is stable even when the magnitude is not. A finding that survives MAUP is one where the same place is extreme across scales and zonings, even though its rate changes. Track argmax_cell back to a coordinate and measure how far the maximum moves; if it stays within one cell width, the location is robust and the magnitude is not, which is a perfectly reportable result.

3. Do not use the diagnostic to select a unit. Choosing the aggregation that produces the strongest signal is the failure mode this whole area of statistics is famous for. The unit is chosen in advance from the exposure’s plausible scale and the denominator’s availability, and the diagnostic reports the consequence.

4. Population-weight the grid, not just the cases. A grid cell in open country can have a handful of cases and almost no population, producing an enormous rate that is pure small-number noise. Apply a minimum-denominator rule before taking maxima, and state it — this is the same instability that motivates the smoothing methods in Bayesian Disease Mapping & Rate Smoothing.

Scale Effect With the Zoning Spread at Each Scale Maximum observed rate against aggregation cell size, from one kilometre to eight kilometres. The central line falls from about 214 to about 62 as cells grow. At each scale a vertical bar shows the spread across four different grid offsets of identical cell size: 55 units of spread at one kilometre, narrowing to 17 at eight kilometres. The scale effect dominates but the zoning spread at the finest scale is a quarter of the headline value. Both effects, on one chart 0 130 260 max rate zoning spread 55 spread 17 1 km 2 km 4 km 8 km aggregation cell size A headline rate quoted without its zoning bar is a point estimate of an interval

5. Repeat the diagnostic on the administrative units actually used. Grids are convenient for the sensitivity family but nobody publishes on them. Where possible, include the real alternatives — block groups, tracts, ZCTAs, counties — as additional rows, since those are the aggregations a reader might have chosen instead.

6. Include the administrative units in the family, not just grids. Grids are convenient for a controlled sweep, but readers compare against tracts, ZCTAs and counties because those are what other publications use. Adding them as extra rows shows how the finding behaves on the aggregations somebody else will reach for.

The scale effect reaches beyond maps and into every correlation computed on areal data, which is where it does the most damage:

Correlation Between Two Variables at Four Aggregations The correlation between a deprivation index and a disease rate, computed on the same underlying data at four aggregation levels. At block-group level it is 0.31, at tract level 0.44, at ZCTA level 0.58 and at county level 0.79. The relationship is identical throughout; only the aggregation changed. A correlation quoted without its unit is not a finding. One relationship, four correlations 0 0.5 1.0 0.31 block group 0.44 tract 0.58 ZCTA 0.79 county A correlation quoted without its aggregation unit is not a finding

Compliance Notes

  • Publish the zoning spread with the headline figure. It is the honest uncertainty attached to a number that is otherwise reported as exact, and omitting it overstates precision in exactly the way the map’s appearance already encourages.
  • Record the unit choice and its justification before the analysis runs. A unit selected after seeing the diagnostic is a selected result, and the run log should show the order.
  • Do not interpolate to a finer unit to escape the problem. Splitting a coarse count across fine units produces a fine-looking map with no additional information and a MAUP profile that is now dominated by the interpolation assumption rather than by the data.
  • State the minimum-denominator rule applied before taking maxima, since the headline figure depends on it as much as on the unit.