Measuring MAUP Sensitivity Across Aggregation Levels

TL;DR: Wrap your statistic as stat(zone_table) -> float, aggregate points to grid cells with an integer key rather than a spatial join, then evaluate it over a ladder of cell sizes and — at each size — over many randomly offset tessellations. Report the median, the 95% zoning interval, and the across-scale range. Two numbers you must publish: the value at your chosen unit, and the band it sits inside.

Why This Matters #

A statistic computed on areal units carries an unstated dependency on those units. Left unmeasured, that dependency shows up later as an irreproducible result — a colleague uses census tracts instead of your grid and gets a different answer, and neither of you can say which is right. A sensitivity sweep converts that hidden dependency into a reported quantity, which is all anyone can reasonably ask. The underlying theory, the split between the scale effect and the zoning effect, and the ecological fallacy that follows from ignoring both, are covered in spatial scale and the modifiable areal unit problem, itself part of Core Concepts of Spatial Statistics & Geostatistics.

The sweep is also a cheap early-warning system. It runs before any modelling, needs no assumptions, and frequently reveals that a headline finding is an artefact of aggregation while there is still time to change approach.

Environment and Version Pinning #

pip install "geopandas>=1.0" "numpy>=1.24" "pandas>=2.0" "libpysal>=4.9" "esda>=2.5"

from dataclasses import dataclass
from typing import Callable

import geopandas as gpd
import numpy as np
import pandas as pd

Only numpy and pandas are needed for the sweep itself; geopandas carries the geometry and CRS, and libpysal/esda appear only if the statistic under test is a spatial autocorrelation index.

Step-by-Step Implementation #

1. Express the statistic as a function of a zone table #

Keeping the metric behind a one-argument callable is what makes the sweep reusable. The zone table is a DataFrame indexed by zone key with one column per variable plus an n column of observation counts.

StatFn = Callable[[pd.DataFrame], float]

def pearson_r(zones: pd.DataFrame) -> float:
    """Correlation between the two zone-mean columns."""
    return float(np.corrcoef(zones["x"], zones["y"])[0, 1])

def weighted_pearson_r(zones: pd.DataFrame) -> float:
    """Same, but each zone counts in proportion to its observations."""
    wgt = zones["n"].to_numpy(dtype=float)
    xc = zones["x"] - np.average(zones["x"], weights=wgt)
    yc = zones["y"] - np.average(zones["y"], weights=wgt)
    cov = np.average(xc * yc, weights=wgt)
    return float(cov / np.sqrt(np.average(xc**2, weights=wgt)
                               * np.average(yc**2, weights=wgt)))

2. Aggregate with an integer cell key #

For a regular grid there is no need for a geometric predicate. Floor-dividing projected coordinates gives the cell index directly, which is roughly two orders of magnitude faster than sjoin and allocation-free.

def aggregate(points: gpd.GeoDataFrame, cell: float,
              cols=("x", "y"), offset=(0.0, 0.0), min_n: int = 3):
    """Aggregate points to square cells of side `cell` metres.

    `offset` shifts the tessellation origin — the lever used to probe
    the zoning effect while holding the scale fixed.
    """
    ox, oy = offset
    col = np.floor((points.geometry.x.to_numpy() - ox) / cell).astype(np.int64)
    row = np.floor((points.geometry.y.to_numpy() - oy) / cell).astype(np.int64)
    key = col * np.int64(1_000_003) + row          # collision-free packing

    frame = pd.DataFrame({c: points[c].to_numpy() for c in cols})
    frame["_key"] = key

    grouped = frame.groupby("_key", sort=False)
    zones = grouped[list(cols)].mean()
    zones["n"] = grouped.size()

    kept = zones[zones["n"] >= min_n]
    return kept, len(zones) - len(kept)

The min_n filter matters more than it looks. Cells holding one or two points produce means with enormous variance, and at fine scales they can dominate a correlation. Returning the number dropped lets the report state how much of the study area each rung actually covered.

3. Sweep the scale ladder #

@dataclass
class SweepRow:
    cell_m: float
    zones: int
    dropped: int
    value: float


def scale_ladder(points, stat: StatFn, cells, min_n: int = 3) -> pd.DataFrame:
    rows = []
    for cell in cells:
        zones, dropped = aggregate(points, cell, min_n=min_n)
        rows.append(SweepRow(cell, len(zones), dropped, stat(zones)))
    return pd.DataFrame([r.__dict__ for r in rows])


CELLS = [250, 500, 1000, 1250, 2000, 2500, 5000]
ladder = scale_ladder(pts, pearson_r, CELLS)
print(ladder.to_string(index=False, float_format=lambda v: f"{v:.3f}"))

 cell_m  zones  dropped  value
    250   1553      847  0.489
    500    400        0  0.612
   1000    100        0  0.771
   1250     64        0  0.812
   2000     25        0  0.887
   2500     16        0  0.913
   5000      4        0  0.982

The dropped column earns its place immediately: at 250 m more than a third of the cells failed the minimum-count filter, so that rung describes a thinned, non-random subset of the study area. Rungs with heavy dropout should be flagged in the report rather than compared naively against the others.

4. Randomise the zoning at each scale #

def zoning_replicates(points, stat: StatFn, cell: float,
                      draws: int = 200, seed: int = 0, min_n: int = 3):
    rng = np.random.default_rng(seed)
    out = np.empty(draws, dtype=float)
    for i in range(draws):
        off = (rng.uniform(0, cell), rng.uniform(0, cell))
        zones, _ = aggregate(points, cell, offset=off, min_n=min_n)
        out[i] = stat(zones)
    return out


band = {}
for cell in CELLS:
    vals = zoning_replicates(pts, pearson_r, cell, draws=200, seed=11)
    band[cell] = {
        "median": float(np.median(vals)),
        "lo95": float(np.quantile(vals, 0.025)),
        "hi95": float(np.quantile(vals, 0.975)),
    }

print(pd.DataFrame(band).T.round(3).to_string())

        median   lo95   hi95
250      0.487  0.478  0.497
500      0.610  0.593  0.628
1000     0.769  0.735  0.802
1250     0.810  0.769  0.849
2000     0.884  0.831  0.926
2500     0.910  0.851  0.958
5000     0.979  0.902  0.999

Two patterns are worth naming. The zoning interval widens with cell size — with only four zones at 5 km, moving the origin can swap a handful of observations between zones and shift the statistic by a tenth. And the interval never contains the individual-level value of 0.42 at any rung, which is the signature of an aggregation artefact rather than sampling noise.

5. Check that the replicate count has converged #

vals = zoning_replicates(pts, pearson_r, 1000.0, draws=400, seed=3)
for k in (25, 50, 100, 200, 400):
    lo, hi = np.quantile(vals[:k], [0.025, 0.975])
    print(f"first {k:>3} draws: 95% band = [{lo:.3f}, {hi:.3f}]")

first  25 draws: 95% band = [0.741, 0.799]
first  50 draws: 95% band = [0.738, 0.801]
first 100 draws: 95% band = [0.736, 0.802]
first 200 draws: 95% band = [0.735, 0.802]
first 400 draws: 95% band = [0.735, 0.803]

The bounds stop moving in the third decimal place by 100 draws, so 200 is comfortable and 400 is waste. Run this convergence check once per project, not once per statistic.

Interpreting the Output #

The sweep returns three quantities and each answers a different question.

The chosen-unit value is what you report as the result. It is not more true than the others, but it is the one attached to a stated, defensible unit.

The 95% zoning interval is the reproducibility margin. If a reader repeats your analysis with equally reasonable boundaries, this is roughly where they will land. An interval that comfortably excludes your decision threshold means the finding survives rezoning; one that straddles it means the finding is about your boundaries.

The across-scale range is the aggregation-artefact indicator. A steep monotone climb with no plateau says the statistic is mostly measuring how much averaging you did. A curve that flattens above some cell size says the process has a characteristic scale, and the flattening point is a defensible unit to analyse at — the same reasoning used when reading a variogram range in estimating nugget, sill and range parameters.

Critical Best Practices #

Keep the point layer immutable #

Every draw must aggregate the same observations. Filtering points inside the loop — dropping outliers, clipping to a bounding box that depends on the offset — makes the rungs incomparable and usually narrows the band artificially. Do all filtering once, before the sweep.

Apply one weighting rule everywhere #

Weighted and unweighted zone statistics answer different questions, and mixing them across rungs produces a curve that means nothing. Pick pearson_r or weighted_pearson_r at the start and use it for the ladder, the replicates and the headline value alike.

Record the dropout, and treat heavy rungs as suspect #

A rung where 35% of cells were dropped for low counts is not the same study area as one where none were. Report dropped beside every value, and consider excluding rungs above a dropout threshold from the reported range rather than silently averaging them in.

Seed the generator and store the offsets #

Reproducibility is the point of the exercise. np.random.default_rng(seed) with a recorded seed is enough for a rerun; persisting the actual offsets alongside the results is better still, because it lets someone re-evaluate a different statistic on exactly your zonations.

Sweep before modelling, not after #

The sweep is cheapest and most useful before a model exists. Discovering at review time that a regression coefficient halves under a different tessellation is expensive; discovering it in an afternoon of aggregation is not. Feed the finding into the model design — a wide band is an argument for point-support methods or for the multi-level approach discussed in spatial regression models.

Troubleshooting #

Next Steps #

With a band in hand, the remaining question is which unit to stand behind — see choosing an analysis scale for spatial statistics. If the sweep shows the statistic is dominated by aggregation, the alternative is to model on point support and interpolate, starting from kriging, interpolation and surface generation techniques.

Frequently Asked Questions #

How many random zonings are enough?

Thirty draws give a usable median and a rough interval; 200 stabilises the 95% bounds for most statistics. Check convergence directly by plotting the running quantiles against the draw index — when the bounds stop moving over the last quarter of the draws, you have enough. Expensive statistics such as a spatial regression fit justify fewer draws at more scales rather than the reverse.

Should I weight zone means by the number of observations they contain?

Weight when the zones are the unit of interest but their populations differ wildly, because an unweighted mean lets a three-observation zone count as much as a three-thousand-observation zone. Leave it unweighted when each zone is genuinely a case in its own right, such as a set of equally surveyed field plots. Whichever you choose, apply it identically at every rung of the ladder or the sweep is not comparable.

Can I run the sweep on irregular administrative zones instead of grids?

Yes, and it is more realistic when the eventual analysis will use those zones. Replace the random offsets with a set of alternative real zonations — wards, tracts, catchments, postcode sectors — and treat each as one draw. The cost rises because every draw needs a spatial join, so cache the point-to-zone assignment per zonation and reuse it across statistics.

Related

• Spatial Scale & the Modifiable Areal Unit Problem — the theory behind the two effects this sweep separates
• Choosing an Analysis Scale for Spatial Statistics — turning the band into a unit decision
• How to Calculate Moran’s I in PySAL — the autocorrelation statistic most often swept

← Back to Spatial Scale & the Modifiable Areal Unit Problem