The Modifiable Areal Unit Problem Explained
Problem statement
You have 4,000 households with an income and a yes/no outcome. At the household level the relationship is modest:
print(np.corrcoef(points["income"], points["incident"])[0, 1].round(3))-0.251
Aggregate those same households into 20 km grid cells and report the mean income and the incident rate per cell:
print(np.corrcoef(cells["income"], cells["rate"])[0, 1].round(3))-0.921
Nothing was added. No new data, no new households, no modelling choice. The correlation went from −0.25 to −0.92 because of the size of the boxes the data was put in.
This is the modifiable areal unit problem (MAUP), and it is not a subtle statistical caveat. It is the reason two competent analysts can produce contradictory results from the same dataset and both be right.
Quick answer
MAUP has two independent halves, and you have to defend against both:
| Effect | What changes | Example |
|---|---|---|
| Scale effect | the size of the zones | 5 km cells vs 20 km cells |
| Zoning effect | the boundaries, at a fixed size | the same 20 km grid, shifted 5 km |
Measured on the same 4,000 points:
cell size zones corr
individual 4000 -0.251
5 km 393 -0.538
10 km 100 -0.763
20 km 25 -0.921
25 km 16 -0.934
20 km grid, shifted zones corr
offset (0.0, 0.0) 25 -0.921
offset (5.0, 5.0) 36 -0.845
offset (2.5, 7.5) 36 -0.748
The scale effect moved the correlation by 0.68. The zoning effect — same cell size, grid shifted a few kilometres — moved it by 0.17.
Step-by-step solution
1. Recognise when you are exposed
MAUP applies whenever you aggregate point or individual data into areas and then analyse the areas. That covers most of applied spatial analysis:
- rates per administrative unit (crime, disease, unemployment)
- choropleth maps of anything derived from individual records
- correlations, regressions or clustering computed on zone-level values
- "hotspots" defined by counting points per polygon
It does not apply when the areal unit is the thing you are studying. If your question is genuinely about municipalities as decision-making entities, the municipality is not an arbitrary container — it is the unit of analysis, and MAUP is not the issue. (Comparing across countries where "municipality" means different things is a different problem, and a real one.)
2. Understand why aggregation strengthens relationships
Averaging removes individual variation. Two households with the same income can have opposite outcomes; two zones with the same mean income have much more similar rates, because each zone's rate is an average over hundreds of households and the noise has cancelled.
What survives aggregation is the systematic part of the relationship. What disappears is the scatter. Correlation measures the ratio of the two, so it rises — often dramatically — as zones get larger:
5 km 393 zones -0.538
20 km 25 zones -0.921
The 20 km correlation is not "clearer". It is a correlation between different quantities — zone means, not households — and it should not be reported as if it described people.
3. Notice where the trend reverses
25 km 16 zones -0.934
50 km 4 zones -0.447
At 50 km there are four zones. With four observations, a correlation is close to meaningless, and it collapses. This is the other end of the trade-off: aggregating strengthens the signal until there are too few units for the statistic to mean anything.
Between those extremes there is no "correct" answer, only a curve. Which is the point.
4. Test the zoning effect explicitly
The scale effect is well known. The zoning effect is the one that catches people, because it needs no change of scale at all:
for offset in [(0, 0), (5, 0), (5, 5), (2.5, 7.5)]:
zones, corr = correlate_on_grid(points, cell=20, offset=offset)
print(f"offset {offset} zones {zones} corr {corr:.3f}")offset (0, 0) zones 25 corr -0.921
offset (5, 0) zones 30 corr -0.876
offset (5, 5) zones 36 corr -0.845
offset (2.5, 7.5) zones 36 corr -0.748
The same cell size, the same points, the grid nudged a few kilometres — and the correlation moves by 0.17. Nobody chose those boundaries to produce that result; they are arbitrary, and arbitrary choices are exactly what MAUP exploits.
5. Report the sensitivity, not a single number
The defensible response is not to find the "right" zoning. It is to show that your conclusion survives reasonable alternatives:
results = [correlate_on_grid(points, cell=c)[1] for c in (5, 10, 15, 20)]
print(f"corr across cell sizes 5–20 km: {min(results):.2f} to {max(results):.2f}")corr across cell sizes 5–20 km: -0.92 to -0.54
"Negative at every scale we tested, between −0.54 and −0.92" is an honest, useful finding. "−0.92" alone is a number chosen by the grid.
Code examples
Example 1 — measuring your own exposure to the scale effect
import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import box
def aggregate_to_grid(points, cell, *, offset=(0.0, 0.0), min_n=5,
value="income", outcome="incident"):
"""Aggregate points into a square grid and return one row per populated cell."""
minx, miny, maxx, maxy = points.total_bounds
xs = np.arange(minx - cell, maxx + cell, cell) + offset[0]
ys = np.arange(miny - cell, maxy + cell, cell) + offset[1]
cells = [box(x, y, x + cell, y + cell) for x in xs for y in ys]
grid = gpd.GeoDataFrame(geometry=cells, crs=points.crs).reset_index(names="zid")
joined = gpd.sjoin(points, grid, predicate="within")
zones = joined.groupby("zid").agg(
value=(value, "mean"),
rate=(outcome, "mean"),
n=(outcome, "size"),
)
return zones[zones["n"] >= min_n] # tiny zones give unstable rates
def scale_sweep(points, sizes, **kwargs):
rows = []
for cell in sizes:
zones = aggregate_to_grid(points, cell, **kwargs)
rows.append({
"cell": cell,
"zones": len(zones),
"median_n": int(zones["n"].median()),
"corr": round(np.corrcoef(zones["value"], zones["rate"])[0, 1], 3),
})
individual = np.corrcoef(points["income"], points["incident"])[0, 1]
print(f"individual level: {individual:.3f}\n")
return pd.DataFrame(rows)
print(scale_sweep(points, [5, 10, 12.5, 20, 25, 50]).to_string(index=False))individual level: -0.251
cell zones median_n corr
5.0 393 8 -0.538
10.0 100 38 -0.763
12.5 64 61 -0.829
20.0 25 159 -0.921
25.0 16 249 -0.934
50.0 4 999 -0.447
The median_n column is the honest one. At 25 km each "observation" is 249 households averaged together, which is why the correlation is so high and why it says nothing about any individual.
Example 2 — measuring the zoning effect
def zoning_sweep(points, cell, offsets, **kwargs):
rows = []
for offset in offsets:
zones = aggregate_to_grid(points, cell, offset=offset, **kwargs)
rows.append({
"offset": f"{offset[0]}, {offset[1]}",
"zones": len(zones),
"corr": round(np.corrcoef(zones["value"], zones["rate"])[0, 1], 3),
})
frame = pd.DataFrame(rows)
spread = frame["corr"].max() - frame["corr"].min()
print(frame.to_string(index=False))
print(f"\nspread from zoning alone: {spread:.3f}")
return frame
zoning_sweep(points, cell=20,
offsets=[(0, 0), (5, 0), (10, 0), (0, 5), (5, 5), (10, 10), (2.5, 7.5)]) offset zones corr
0, 0 25 -0.921
5, 0 30 -0.876
10, 0 30 -0.857
0, 5 30 -0.905
5, 5 36 -0.845
10, 10 36 -0.854
2.5, 7.5 36 -0.748
spread from zoning alone: 0.173
If your headline finding moves by less than this spread when you shift the grid, the finding is a property of the grid.
Example 3 — reporting a range instead of a point estimate
def maup_report(points, sizes=(5, 10, 15, 20), offsets=((0, 0), (5, 5), (2.5, 7.5))):
values = []
for cell in sizes:
for offset in offsets:
zones = aggregate_to_grid(points, cell, offset=offset)
if len(zones) < 10:
continue # too few units to correlate
values.append(np.corrcoef(zones["value"], zones["rate"])[0, 1])
values = np.array(values)
print(f"{len(values)} aggregations tested")
print(f" correlation range {values.min():.3f} to {values.max():.3f}")
print(f" median {np.median(values):.3f}")
print(f" sign consistent {bool(np.all(np.sign(values) == np.sign(values[0])))}")
return values
maup_report(points)12 aggregations tested
correlation range -0.921 to -0.522
median -0.807
sign consistent True
sign consistent: True is the claim that survives. The direction of the relationship is robust to every aggregation tested; the magnitude is not, and quoting one magnitude would be quoting the grid.
Explanation
Why this is not the same as the ecological fallacy
They are related and they are not identical.
MAUP is about the arbitrariness of the zones: different zone definitions produce different statistics.
The ecological fallacy is about inference: concluding something about individuals from a relationship measured between areas.
The example above contains both. The scale sweep is MAUP. Reading "−0.92" as "richer households have far fewer incidents" is the ecological fallacy, and the household-level figure of −0.25 is the correction.
You can have MAUP without the fallacy — comparing zones to zones, honestly, while still being sensitive to zone definition. Avoiding the fallacy does not protect you from MAUP.
Why "just use the smallest zones" does not work
It is the obvious response and it fails for three reasons:
- Small zones have unstable rates. A zone with five households has an incident rate that can only be 0, 0.2, 0.4… Its variance is dominated by sample size, not by anything real.
- Small zones are suppressed. Statistical agencies withhold values for small units to protect privacy, so the data is missing exactly where you wanted it.
- Small zones are still arbitrary. A 1 km grid is as modifiable as a 20 km one. The zoning effect does not disappear with scale.
The right size is set by the process you are studying: a walking-distance effect wants zones of a few hundred metres, a labour-market effect wants travel-to-work areas. That is a substantive argument, not a statistical one.
Why aggregation inflates correlation
Correlation is systematic variation divided by total variation. Averaging leaves the systematic part roughly intact while shrinking the random part by roughly the square root of the group size.
With a median of 159 households per zone at 20 km, the noise in each zone's rate is roughly one-twelfth of an individual's. The signal is unchanged. So the ratio — the correlation — rises, and it will keep rising until the number of zones becomes too small for the statistic to be estimated at all. That is the −0.447 at 50 km with four zones.
Why this matters for hotspots and clustering too
MAUP is usually discussed with correlation, but every zone-based statistic inherits it:
- Moran's I is computed on zones and depends on both their size and their adjacency structure.
- Getis-Ord hotspots shift with the zone definition — a hotspot can appear or vanish under a different grid.
- Choropleth classification compounds it: arbitrary zones plus arbitrary class breaks.
- Hexagonal binning does not solve it. Hexagons remove the directional bias of a square grid, not its modifiability.
The one genuine escape is not to aggregate: analyse the points themselves, with kernel density or a point-pattern method, and let the surface be continuous rather than boxed.
Edge cases or notes
- Aggregating to natural units does not fix it. Watersheds and travel-to-work areas are less arbitrary than a grid, but they are still one of many possible partitions.
- The direction of a relationship is usually more robust than its magnitude. Report the sign confidently and the size with a range.
- Rates need denominators from the same zones. A count per zone divided by a population from differently-shaped zones adds a second, worse error on top of MAUP.
- Weighted correlation is not a fix. Weighting zones by population makes results a little more stable, and they still move with the boundaries.
- Very small zones invite the small-numbers problem, where the highest and lowest rates are always the least populated units. Empirical Bayes smoothing addresses that, not MAUP.
- MAUP has no p-value. It is not something you test for and rule out; it is a property of the data structure that you either bound or ignore.
- Say which zones you used, always. "Correlation of −0.8 at the 10 km grid scale" is reproducible. "Correlation of −0.8" is not.
Internal links
- Spatial autocorrelation explained — another statistic that depends on the zone definition
- Kernel density explained — a continuous alternative to aggregating into zones
- How to bin points into hexagons in Python — a better grid, not an escape
- How to find hotspots with Getis-Ord Gi* in Python — hotspots move when the zones do
- Choropleth classification explained — arbitrary breaks on top of arbitrary zones
- How to aggregate spatial data by region in GeoPandas — the mechanics of the aggregation itself
- How to download administrative boundaries in Python — where the arbitrary zones usually come from
- How to count points in polygons with GeoPandas — the operation that starts the problem
FAQ
What is the modifiable areal unit problem in one sentence?
Statistics computed on aggregated areas depend on how the areas were drawn, so the same data produces different answers under different zone definitions.
What is the difference between the scale effect and the zoning effect?
Scale is the size of the zones; zoning is where their boundaries fall at a fixed size. Both change the result, and the zoning effect is the one people forget to test.
Does using smaller zones solve it?
No. Small zones have unstable rates, are often suppressed for privacy, and are just as arbitrary. Zone size should be chosen from the process you are studying.
Is MAUP the same as the ecological fallacy?
No. MAUP is about the arbitrariness of zones; the ecological fallacy is about inferring individual behaviour from zone-level relationships. They frequently occur together.
Do hexagons fix it?
They remove the directional bias of a square grid, which is worth having. They do not remove modifiability — a hexagonal grid can be shifted and resized just as freely.
How do I report results honestly?
Sweep several scales and several offsets, report the range and whether the sign is consistent, and always state which zones produced the headline number.
Can I avoid it entirely?
Only by not aggregating — analysing individual records, or using a continuous surface such as kernel density instead of counts per polygon.