The Modifiable Areal Unit Problem Explained

Problem statement

Take Leeds' 2,607 Census 2021 output areas — the smallest census geography, about 120 households each — and ask how closely two things go together: the share of households with no car, and the share of residents in bad or very bad health. At output-area level the relationship is moderate:

print(np.corrcoef(oas["pct_no_car"], oas["pct_bad_health"])[0, 1].round(3))

0.43

Add the same output areas up into 5 km grid cells, recompute both percentages from the summed counts, and ask again:

cells = aggregate_to_grid(oas, 5000)
print(np.corrcoef(cells["pct_no_car"], cells["pct_bad_health"])[0, 1].round(3))

0.671

Nothing was added. No new data, no new households, no modelling choice. The correlation changed because of the size of the boxes the data was put in — and shift that grid, or change its size again, and it changes again, in either direction.

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:

Measured on the same output areas (see The data used on this page):

for cell in (1000, 2000, 3000, 5000, 8000):
    cells = aggregate_to_grid(oas, cell)
    print(f"{cell / 1000:>4.0f} km  {len(cells):4} zones  corr {np.corrcoef(cells['pct_no_car'], cells['pct_bad_health'])[0, 1]:.3f}")
for offset in ((0, 0), (4000, 0), (4000, 4000)):
    cells = aggregate_to_grid(oas, 8000, offset=offset)
    print(f"8 km grid shifted {offset}  {len(cells):3} zones  corr "
          f"{np.corrcoef(cells['pct_no_car'], cells['pct_bad_health'])[0, 1]:.3f}")

   1 km   396 zones  corr 0.503
   2 km   136 zones  corr 0.454
   3 km    72 zones  corr 0.429
   5 km    32 zones  corr 0.671
   8 km    16 zones  corr 0.360
8 km grid shifted (0, 0)   16 zones  corr 0.360
8 km grid shifted (4000, 0)   16 zones  corr 0.180
8 km grid shifted (4000, 4000)   17 zones  corr 0.139

The scale effect and the zoning effect each moved the correlation by more than the output-area value itself is worth quoting to. The sign never changed; the size did, and not in a straight line.

When you are exposed, and testing scale and zoning

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.

Note that the "fine" level here is already an aggregate. Household-level census records are not published; an output area of about 120 households is the finest public unit, so every census analysis starts one step up the scale ladder.

2. Understand what aggregation does to a relationship

Averaging removes local variation. Two neighbouring output areas can differ sharply; the 5 km cells that contain them are each an average over dozens of output areas, so the idiosyncratic part cancels and what is shared survives. That often strengthens a correlation — here from 0.43 to 0.67 at 5 km.

But it is not a law. Aggregation also mixes different places: a 5 km cell can hold a car-free inner-city estate and a commuter village, and whether their combination strengthens or weakens the relationship depends on where the cell boundaries fall. That is why the scale sweep above is not monotonic — 0.50 at 1 km, 0.45 at 2 km, 0.43 at 3 km, 0.67 at 5 km.

A 5 km correlation is not "clearer". It is a correlation between different quantities — cell totals, not neighbourhoods — and it should not be reported as if it described people.

3. Notice where it breaks down

At 8 km there are 16 or 17 cells. With that few observations a correlation is barely estimable, and it moved from 0.36 to 0.14 when the grid was shifted by half a cell. At 12 km, with 11 cells, it is −0.05: the sign of the relationship has flipped. This is the other end of the trade-off: aggregating strengthens the shared signal until there are too few units for the statistic to mean anything — and then it can say the opposite.

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 (Example 2): the same cell size, the same output areas, the grid nudged by half a cell — and the correlation moves. 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 (Example 3). "Positive at every scale and offset from 1 to 5 km, between 0.41 and 0.67, and unreliable beyond that" is an honest, useful finding. "0.67" alone is a number chosen by the grid.

Car access and health across Leeds output areas, regrouped

The data used on this page

Every example runs on real data: the 2,607 Census 2021 output areas of Leeds, located by the population-weighted centroids ONS publishes with the boundaries, and two census tables from Nomis — car or van availability (TS045) and general health (TS037). aggregate_to_grid adds output areas up into square cells and recomputes both percentages from the summed counts, which is how zone-level census figures are made.

import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query")
rows, offset = [], 0
while True:                                   # the service pages at 2,000 records
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "OA21CD,BNG_E,BNG_N",
                                     "returnGeometry": "false", "f": "json", "resultOffset": offset}, timeout=120).json()
    rows += [f["attributes"] for f in page["features"]]
    if not page.get("exceededTransferLimit"):
        break
    offset += len(page["features"])
oas = pd.DataFrame(rows).set_index("OA21CD").sort_index()


def census(dataset, dimension, categories):
    """One Nomis table for the output areas, 400 codes per request."""
    parts = [pd.read_csv(f"https://www.nomisweb.co.uk/api/v01/dataset/{dataset}.data.csv?date=latest"
                         f"&geography={','.join(oas.index[i:i + 400])}&{dimension}={categories}"
                         f"&measures=20100&select=geography_code,{dimension},obs_value")
             for i in range(0, len(oas), 400)]
    return pd.concat(parts).pivot_table(index="GEOGRAPHY_CODE", columns=dimension.upper(), values="OBS_VALUE")


cars = census("NM_2063_1", "c2021_cars_5", "0,1")           # 0 all households, 1 no car or van
health = census("NM_2055_1", "c2021_health_6", "0,4,5")     # 0 all residents, 4 bad, 5 very bad
oas["households"], oas["no_car"] = cars[0], cars[1]
oas["residents"], oas["bad_health"] = health[0], health[4] + health[5]
oas["pct_no_car"] = 100 * oas["no_car"] / oas["households"]
oas["pct_bad_health"] = 100 * oas["bad_health"] / oas["residents"]


def aggregate_to_grid(units, cell, *, offset=(0, 0), min_households=50):
    """Add units up into square cells of `cell` metres and recompute the percentages."""
    gx = (units["BNG_E"] - offset[0]) // cell
    gy = (units["BNG_N"] - offset[1]) // cell
    grid = units.groupby([gx, gy])[["households", "no_car", "residents", "bad_health"]].sum()
    grid = grid[grid["households"] >= min_households]           # a few houses give unstable rates
    grid["n_oas"] = units.groupby([gx, gy]).size()
    grid["pct_no_car"] = 100 * grid["no_car"] / grid["households"]
    grid["pct_bad_health"] = 100 * grid["bad_health"] / grid["residents"]
    return grid


print(f"{len(oas):,} output areas, {int(oas['households'].sum()):,} households, "
      f"{int(oas['residents'].sum()):,} residents; median {oas['households'].median():.0f} households each")

2,607 output areas, 341,381 households, 811,904 residents; median 127 households each

Example 1 — measuring your own exposure to the scale effect

def scale_sweep(units, sizes):
    rows = [{"cell": "output area", "zones": len(units), "median_oas": 1,
             "corr": round(np.corrcoef(units["pct_no_car"], units["pct_bad_health"])[0, 1], 3)}]
    for cell in sizes:
        zones = aggregate_to_grid(units, cell)
        rows.append({"cell": f"{cell / 1000:g} km", "zones": len(zones), "median_oas": int(zones["n_oas"].median()),
                     "corr": round(np.corrcoef(zones["pct_no_car"], zones["pct_bad_health"])[0, 1], 3)})
    return pd.DataFrame(rows)


print(scale_sweep(oas, [1000, 2000, 3000, 5000, 8000, 12000]).to_string(index=False))

       cell  zones  median_oas   corr
output area   2607           1  0.430
       1 km    396           4  0.503
       2 km    136          11  0.454
       3 km     72          19  0.429
       5 km     32          40  0.671
       8 km     16          61  0.360
      12 km     11          69 -0.053

The median_oas column is the honest one. At 8 km each "observation" is dozens of output areas added together, which is why the correlation can say nothing about any single neighbourhood, and at 12 km there are too few cells left to estimate one at all.

Example 2 — measuring the zoning effect

def zoning_sweep(units, cell, offsets):
    rows = []
    for offset in offsets:
        zones = aggregate_to_grid(units, cell, offset=offset)
        rows.append({"offset": f"{offset[0]:g}, {offset[1]:g}", "zones": len(zones),
                     "corr": round(np.corrcoef(zones["pct_no_car"], zones["pct_bad_health"])[0, 1], 3)})
    frame = pd.DataFrame(rows)
    print(frame.to_string(index=False))
    print(f"\nspread from zoning alone: {frame['corr'].max() - frame['corr'].min():.3f}")
    return frame


zoning_sweep(oas, cell=3000,
             offsets=[(0, 0), (1000, 0), (1500, 0), (0, 1500), (1500, 1500), (750, 2250)])

    offset  zones  corr
      0, 0     72 0.429
   1000, 0     72 0.480
   1500, 0     71 0.518
   0, 1500     74 0.413
1500, 1500     70 0.540
 750, 2250     70 0.560

spread from zoning alone: 0.147

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(units, sizes=(1000, 2000, 3000, 5000), offsets=((0, 0), (0.5, 0.5), (0.25, 0.75))):
    values = []
    for cell in sizes:
        for fx, fy in offsets:                   # offsets as fractions of the cell
            zones = aggregate_to_grid(units, cell, offset=(fx * cell, fy * cell))
            if len(zones) < 10:
                continue                         # too few units to correlate
            values.append(np.corrcoef(zones["pct_no_car"], zones["pct_bad_health"])[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(oas)

12 aggregations tested
  correlation range   0.409 to 0.671
  median              0.521
  sign consistent     True

sign consistent: True is the claim that survives. The direction of the relationship is robust to every aggregation tested down to 5 km cells; the magnitude is not, and quoting one magnitude would be quoting the grid. Extend the sweep to 12 km and the sign fails too (Example 1) — which is why the report skips grids with fewer than ten zones.

MAUP, the ecological fallacy and inflated correlation

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 invites both. The scale sweep is MAUP. Reading the 5 km figure of 0.67 as "people without cars are in much worse health" is the ecological fallacy: the census does not link a household's car to a resident's health, and a place with many car-free households and many ill residents need not contain a single car-free person in bad health. Even the output-area figure, 0.43, describes places, not people.

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 a car-free share that can only be 0%, 20%, 40%… 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, roughly by the square root of the group size — with a median of 40 output areas per cell at 5 km, by about a factor of six. So the correlation tends to rise, from 0.43 to 0.67 here.

"Tends to" is doing work. Averaging also mixes unlike places into one cell, which can dilute the shared signal as easily as it concentrates it: that is why 2 km and 3 km came out lower than 1 km. And once there are too few cells for the statistic to be estimated at all, it wanders — 0.36 at 8 km, −0.05 at 12 km.

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.

Natural units, denominators and small numbers

• 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.67 on a 5 km grid of Census 2021 output areas" is reproducible. "Correlation of 0.67" 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
• Census geographies explained — the zones most analyses inherit
• Small numbers explained — why rates in small zones jump around
• Exposure assessment explained — the reporting scale decides who counts as exposed

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.