Spatial Epi Public Health GIS

Land-Use Regression for Air Pollution Exposure

This guide is part of the Environmental Exposure & Risk Surface Modeling section. Land-use regression is the method to reach for when the exposure varies over hundreds of metres rather than tens of kilometres, and a smooth interpolated surface would average away exactly the contrast the study is about.

Concept & Epidemiological Alignment

Land-use regression abandons the geostatistical premise. Instead of assuming nearby places are similar and borrowing strength from neighbouring measurements, it assumes concentration is a function of the built environment surrounding a location — road length within 100 m, industrial land within 500 m, population density within 1 km, elevation, distance to the nearest major road. Those predictors are computed at each monitoring site, a regression is fitted against the measured concentrations, and the fitted model is then evaluated at any location where the predictors can be computed, which is everywhere.

The epidemiological appeal is direct. Traffic-related pollutants — nitrogen dioxide, ultrafine particles, black carbon — fall by half within one to two hundred metres of a busy road. A cohort study of asthma incidence lives or dies on whether the exposure model can distinguish a residence on the corridor from one three streets back. Kriging a regulatory network cannot: the monitors are deliberately sited away from roadsides to measure background, and the fitted covariance smooths over the gradient. Land-use regression can, because the gradient is encoded in the predictors rather than inferred from the measurements.

The trade is severe and must be stated up front. The model is fitted on a small number of sites — 40 to 80 in a typical metropolitan campaign — against dozens or hundreds of candidate predictors generated by crossing a handful of land-use variables with a ladder of buffer radii. That ratio is a recipe for a model that fits its own monitors superbly and predicts poorly a kilometre away. Every serious design decision in this topic exists to manage that one risk.

A second assumption deserves naming: the model is a statistical description of the sampling campaign’s period and place, not a dispersion model. It contains no physics. It will not extrapolate to a source that did not exist during the campaign, it will not respond to a highway rebuilt after the fit, and its coefficients have no meaning outside the range of predictor values that were sampled.

Predictor Families and What Each One Captures

Predictor family Typical buffers Captures Expected sign Common failure
Road length by class 25–1000 m Direct traffic emissions Positive Road classes inconsistent between jurisdictions
Traffic volume × length 50–1000 m Emission intensity, not just presence Positive Volume counts missing on minor roads
Distance to major road Point measure Near-source decay Negative Undefined where no major road is within range
Land cover: industrial, port 100–5000 m Point-source influence Positive Class definitions change between raster vintages
Land cover: green, water 100–5000 m Dilution, absence of sources Negative Correlated with rurality, so it proxies everything
Population or dwelling density 100–5000 m Domestic and area sources Positive Near-collinear with road density
Elevation, terrain Point or buffer Drainage, inversion layers Context-dependent Sign varies by basin; do not constrain blindly

Crossing seven families with six buffer radii yields more than forty candidates before interactions, and a campaign with 50 sites can support perhaps four or five predictors in the final model.

How One Variable Becomes Six Correlated Candidates A monitoring site at the centre of five concentric buffers at twenty-five, one hundred, three hundred, one thousand and three thousand metres. Road segments cross the buffers. Road length is computed within each buffer, producing five candidate predictors from a single underlying variable. A panel on the right shows that these five candidates are strongly correlated with one another because each larger buffer contains the smaller ones, so a selection procedure that treats them as independent will pick whichever one best fits the noise. The panel notes that only one member of each ladder should enter the model. A buffer ladder is one variable wearing five hats monitoring site with 25, 100, 300, 1000, 3000 m buffers Why this breaks naive selection Each larger buffer contains the smaller ones, so the five road-length candidates correlate at 0.8 and above with each other. A stepwise procedure treats them as five independent chances to fit the residual noise, and picks the one that happens to fit best. Rule: at most one buffer per variable family enters the final model.

Spatial Data Prerequisites

  • Geometry: monitoring sites as points with a stable identifier; a road network with a usable functional class attribute; land-cover raster or polygons; population at block or grid resolution.
  • Coordinate reference system: one projected metric system for every layer. Buffer radii are the model, so a buffer computed in degrees invalidates every predictor at once.
  • Site count: a usable campaign has at least 40 sites; below 30, the model cannot be validated in a way anyone should believe.
  • Site siting: the campaign must sample the range of the predictors, not just the geography. Twenty background sites and two roadside sites cannot estimate a roadside gradient regardless of the total count.
  • Vintage alignment: road network, land cover and population must date from the measurement campaign, not from today. A model fitted with current land cover against a 2015 campaign attributes concentrations to development that had not happened.
  • Predictor completeness: every predictor must be computable at every prediction location. A distance-to-port predictor is undefined for most of a state and cannot enter a state-wide model.

Production Implementation

The routine below builds the buffer predictors, then runs supervised forward selection with two constraints that separate a defensible land-use regression from a data-dredging exercise: an expected sign for each predictor family, and a rule that only one buffer from each family may enter.

# Land-use regression: buffer predictors + supervised selection with sign constraints.
# Pinned: geopandas==0.14.4, numpy==1.26.4, pandas==2.2.2, scikit-learn==1.5.0
import logging
import re

import numpy as np
import pandas as pd
import geopandas as gpd
from sklearn.linear_model import LinearRegression

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

BUFFERS_M = (25, 50, 100, 300, 1000, 3000)
# Direction each family must take; a candidate with the wrong sign is rejected
# outright rather than accepted because it improved the fit.
EXPECTED_SIGN = {"road": +1, "traffic": +1, "industrial": +1, "port": +1,
                 "popdens": +1, "green": -1, "water": -1, "dist_major": -1}
MAX_VIF = 3.0
MIN_R2_GAIN = 0.01


def build_buffer_predictors(sites: gpd.GeoDataFrame, roads: gpd.GeoDataFrame,
                            landcover: gpd.GeoDataFrame,
                            buffers=BUFFERS_M) -> pd.DataFrame:
    """Road length and land-cover area within each buffer, per site."""
    if sites.crs.is_geographic:
        raise ValueError("project to metres — buffer radii are the model")
    sites = sites.sort_values("site_id").reset_index(drop=True)
    cols = {"site_id": sites["site_id"].to_numpy()}
    for r in buffers:
        ring = sites.copy()
        ring["geometry"] = ring.geometry.buffer(r)
        seg = gpd.overlay(roads[["fclass", "geometry"]], ring[["site_id", "geometry"]],
                          how="intersection", keep_geom_type=False)
        seg["m"] = seg.geometry.length
        major = seg[seg["fclass"].isin(("motorway", "trunk", "primary"))]
        cols[f"road_major_{r}"] = (major.groupby("site_id")["m"].sum()
                                   .reindex(sites["site_id"]).fillna(0.0).to_numpy())
        cols[f"road_all_{r}"] = (seg.groupby("site_id")["m"].sum()
                                 .reindex(sites["site_id"]).fillna(0.0).to_numpy())
        lc = gpd.overlay(landcover[["class", "geometry"]], ring[["site_id", "geometry"]],
                         how="intersection", keep_geom_type=False)
        lc["a"] = lc.geometry.area
        for cls, name in (("industrial", "industrial"), ("forest", "green")):
            cols[f"{name}_{r}"] = (lc[lc["class"] == cls].groupby("site_id")["a"].sum()
                                   .reindex(sites["site_id"]).fillna(0.0).to_numpy())
    log.info("built %d candidate predictors for %d sites", len(cols) - 1, len(sites))
    return pd.DataFrame(cols)


def _family(name: str) -> str:
    return re.sub(r"_\d+$", "", name)


def _vif(X: np.ndarray, j: int) -> float:
    others = np.delete(X, j, axis=1)
    if others.shape[1] == 0:
        return 1.0
    m = LinearRegression().fit(others, X[:, j])
    r2 = m.score(others, X[:, j])
    return 1.0 / max(1e-9, 1.0 - r2)


def supervised_forward_selection(X: pd.DataFrame, y: np.ndarray,
                                 max_terms: int | None = None) -> list[str]:
    """Add predictors one at a time, honouring sign, family and collinearity rules."""
    n = len(y)
    cap = max_terms or max(1, n // 10)      # one term per ten sites, hard ceiling
    chosen, used_families, r2 = [], set(), 0.0
    while len(chosen) < cap:
        best = None
        for col in X.columns:
            fam = _family(col)
            if col in chosen or fam in used_families:
                continue
            trial = chosen + [col]
            M = X[trial].to_numpy(float)
            fit = LinearRegression().fit(M, y)
            want = EXPECTED_SIGN.get(fam)
            if want is not None and np.sign(fit.coef_[-1]) != want:
                continue                     # wrong direction: reject, do not keep
            if _vif(M, M.shape[1] - 1) > MAX_VIF:
                continue
            gain = fit.score(M, y) - r2
            if gain > MIN_R2_GAIN and (best is None or gain > best[1]):
                best = (col, gain, fit.score(M, y))
        if best is None:
            break
        chosen.append(best[0])
        used_families.add(_family(best[0]))
        r2 = best[2]
        log.info("selected %-22s cumulative_r2=%.3f terms=%d", best[0], r2, len(chosen))
    log.info("final model: %s", " + ".join(chosen) or "(intercept only)")
    return chosen

The three guards in that loop do most of the work. The sign constraint refuses a predictor whose fitted direction contradicts physical reasoning — a negative coefficient on major road length within 100 m is not a discovery, it is collinearity producing a nonsense partial effect, and accepting it because the fit improved is how untransferable models are built. The one-buffer-per-family rule stops a correlated ladder from entering multiple times. The term cap of one predictor per ten sites is blunt and effective; a 50-site campaign gets at most five terms no matter how good the candidates look.

Parameter Selection & Tuning

The buffer ladder should span the physical scales that matter: tens of metres for near-road, hundreds for neighbourhood sources, kilometres for regional background. Adding more rungs does not add information — adjacent buffers correlate above 0.9 — and does inflate the selection search.

The term cap is the single most effective defence against overfitting. One term per ten sites is a common convention. Treat it as a hard ceiling, not a target: a three-term model that validates is better than a five-term model that does not.

The minimum gain threshold stops the procedure from adding terms that improve the fit by rounding error. One per cent of explained variance is a reasonable floor; the exact value matters less than declaring it before the run.

Collinearity tolerance. A variance inflation factor above about 3 means the new predictor is largely a linear combination of those already in the model, and its coefficient is unstable. Some published models tolerate 5 or 10; the higher the tolerance, the less transferable the result.

Transformations. Traffic-related concentrations are often closer to linear in the square root or logarithm of road length than in the raw metre count. Fix the transformation from the physical reasoning before selection, and apply it uniformly across the family — choosing a transformation per predictor by fit is another route to overfitting.

Edge Cases & Failure Modes

A campaign that samples geography rather than predictor range. If the sites span the metropolitan area evenly but only two sit within 50 m of a major road, the near-road coefficient is estimated from two observations, and its confidence interval will span zero even as the model reports a good overall fit. Design or subset the campaign around the predictor distribution.

Predictors undefined at prediction locations. A distance-to-port term fitted on a coastal campaign has no value inland. Every predictor must be computable everywhere the surface will be evaluated, or the surface has holes that get filled with an implicit extrapolation.

Extrapolation past the sampled range. A model fitted where road length within 100 m never exceeded 800 m will happily predict for a location with 3000 m, and the prediction is pure linear extrapolation. Flag prediction locations that fall outside the training hull of every predictor and treat them as unmodelled.

Vintage drift. Applying a model fitted on a 2016 campaign to 2024 land cover changes the predictor values without changing the coefficients, silently assuming the emission-per-unit-road relationship has not moved. Fleet turnover alone has cut roadside nitrogen dioxide substantially over such a period.

Spatial autocorrelation in the residuals. If the model residuals cluster in space, a systematic source is missing and the standard errors are too small. Test them with the diagnostics in Diagnosing Spatial Autocorrelation in Regression Residuals, and consider adding a kriged residual layer rather than chasing another predictor.

Where the Term Cap Comes From Two curves plotted against the number of predictors in a land-use regression model. The training fit rises steadily and never falls, because adding a predictor can only improve in-sample fit. The held-out fit rises at first, peaks at a small number of terms, and then declines as additional predictors fit campaign-specific noise. The gap between the curves widens continuously and is labelled as the overfitting gap. A marker at one term per ten sites shows the cap, sitting close to the peak of the held-out curve. In-sample fit only ever improves — which is why it decides nothing training fit held-out fit cap: one term per ten sites overfitting gap 0 predictors in the model explained variance

The training hull is the control that stops those extrapolation failures from reaching a published surface, and it is simple to picture:

Inside and Outside the Training Hull A plot of two selected predictors, road length within one hundred metres and industrial area within one kilometre. The campaign's monitoring sites form a cloud occupying the lower-left and centre of the plot, and a box drawn around their minimum and maximum values marks the training hull. Most prediction locations fall inside it. A group of locations beside a new motorway interchange falls far to the right, beyond any sampled road length, and those predictions are linear extrapolation with no data behind them and must be flagged. Outside the hull the model is extrapolating a straight line training hull new interchange flag, do not publish as ordinary road length within 100 m industrial area within 1 km filled points are campaign sites; open points are prediction locations beyond the hull

Compliance & Audit Controls

A land-use regression is a published claim about a place at a time, so its record must let someone rebuild it. Store the campaign dataset and measurement period; the road, land-cover and population sources with vintages; the full candidate list and the buffer ladder; the selection rules including expected signs, family rule, term cap, gain threshold and collinearity tolerance; the selected terms with coefficients and standard errors; the held-out validation statistics; and the training hull of each predictor so extrapolation can be detected later.

Determinism requires a stable site sort and a fixed candidate column order, because forward selection with near-ties is sensitive to ordering. Where the selection is bootstrapped, seed it. The model itself contains no protected data and should be published; the exposure values it generates at residences do, and enter the controls described in the parent section.

Production Implementation Checklist

Frequently Asked Questions

How many monitoring sites do I need? Forty is a working minimum for a metropolitan model and sixty is comfortable. What matters more than the raw count is coverage of the predictor range: a campaign of eighty sites that all sit in background locations cannot estimate a near-road gradient, while forty sites deliberately spread across the road-density distribution can.

Can I use a published model from another city? Only with local validation. Land-use regression coefficients encode a specific fleet, road classification and land-cover taxonomy, and they transfer poorly across those boundaries. Measure at a handful of local sites and check the transferred model’s error before adopting it; the procedure is set out in the transfer guide below.

Should I combine land-use regression with kriging? Frequently, yes. Fitting the regression first and then kriging its residuals captures both the street-scale structure the predictors explain and the regional structure they miss. The combination usually validates better than either alone, and it restores a prediction variance the regression does not provide.

Why did my model select an implausible negative coefficient on road length? Collinearity. When two correlated predictors enter together, the partial coefficients can take large offsetting values with meaningless signs. The sign constraint in the selection loop exists precisely to reject those, and a family that keeps producing them is a sign that the candidate set is too correlated for the site count.

How do I know the model will not be used outside its valid range? Store the minimum and maximum of every selected predictor across the training sites, and evaluate that hull at every prediction location. Locations outside it should be flagged in the output rather than shipped as ordinary predictions, because for those the model is extrapolating linearly with nothing to support it.