Spatial Epi Public Health GIS

Separating Berkson and Classical Error in Area Exposure

This guide, part of Exposure Measurement Error & Uncertainty, takes the most common real situation — an exposure averaged over a census tract and assigned to everyone in it — and decomposes its error into the component that biases the estimate and the component that only widens the interval.

Problem Context & Constraints

An area-level exposure is built in two steps, and each introduces a different error. First, a surface is predicted over the area; that surface is wrong by its prediction error, and the area mean inherits a version of it. Second, the area mean is assigned to every resident, whose individual exposures vary around it.

The second step is Berkson: the truths scatter around the assigned value, which leaves the effect estimate approximately unbiased and widens its interval. The first step is classical: the assigned area mean scatters around the true area mean, which attenuates. Both are present in every area-level exposure, and treating the whole thing as one structure gets the correction wrong in one direction or the other.

Getting the split right matters because the two components respond to different investments. A large classical component means the surface is poorly determined over that area, and more monitors or a better model would fix it. A large Berkson component means the area is internally heterogeneous, and finer geography — block groups instead of tracts, or address-level assignment — would fix it. An agency deciding where to spend has to know which one dominates.

The constraint is that the decomposition needs two quantities that must be estimated separately: the prediction variance of the area mean, which the surface supplies, and the within-area variance of true exposures, which requires the surface at finer resolution than the area.

Two Errors Stacked in One Assigned Value A single census tract shown with a fine exposure surface inside it. The true area mean sits at one level and the predicted area mean sits slightly away from it, and that gap is labelled the classical component arising from surface prediction error. Individual residents are shown scattered above and below the predicted mean according to where they live within the tract, and that spread is labelled the Berkson component. The two are drawn as separate vertical segments stacked on one axis so their relative sizes can be compared, and a note records that they call for different remedies. One assigned number, two different errors, two different fixes one tract, fine surface inside assigned area mean true area mean classical: surface error fix: more monitors, better model Berkson: within-area spread fix: finer geography Only the classical part attenuates the estimate; only the Berkson part responds to finer areas

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4
  • The exposure surface at a resolution meaningfully finer than the analysis areas
  • The surface’s prediction variance, or its cross-validated error where no variance band exists
  • Area polygons with population weights, so within-area spread is weighted by where people live

Step-by-Step Solution

Compute the population-weighted within-area variance from the fine surface, take the prediction variance of the area mean from the surface, and report the two components with their implications.

# Decompose area-level exposure error into classical and Berkson components.
# Pinned: numpy==1.26.4, pandas==2.2.2, geopandas==0.14.4
import logging

import numpy as np
import pandas as pd

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


def area_error_components(cells: pd.DataFrame, area_col: str = "geoid",
                          value_col: str = "exposure",
                          var_col: str = "exposure_var",
                          weight_col: str = "population") -> pd.DataFrame:
    """Per-area classical and Berkson variance components.

    cells: fine-resolution surface cells with an area identifier, predicted value,
           prediction variance and a population weight.
    """
    rows = []
    for geoid, g in cells.groupby(area_col):
        w = g[weight_col].to_numpy(float).clip(min=0.0)
        if w.sum() <= 0:
            w = np.ones(len(g))
        w = w / w.sum()
        z = g[value_col].to_numpy(float)

        # Berkson: how much true exposure varies AROUND the assigned area mean.
        mean_z = float(w @ z)
        berkson = float(w @ (z - mean_z) ** 2)

        # Classical: prediction error of the area mean itself. Cells are strongly
        # correlated, so the weighted mean of cell variances is the defensible
        # conservative figure, not that mean divided by the cell count.
        classical = float(w @ g[var_col].to_numpy(float))

        rows.append({area_col: geoid, "exposure": mean_z,
                     "var_classical": classical, "var_berkson": berkson,
                     "n_cells": int(len(g))})

    df = pd.DataFrame(rows)
    df["var_total"] = df["var_classical"] + df["var_berkson"]
    df["berkson_share"] = df["var_berkson"] / df["var_total"].replace(0, np.nan)
    log.info("decomposition over %d areas: median classical=%.4f median berkson=%.4f "
             "median berkson share=%.2f",
             len(df), df["var_classical"].median(), df["var_berkson"].median(),
             df["berkson_share"].median())
    return df


def attenuation_from_classical(area_exposure: np.ndarray,
                               var_classical: np.ndarray) -> dict:
    """Reliability using ONLY the classical component — Berkson does not attenuate."""
    var_assigned = float(np.nanvar(area_exposure, ddof=1))
    var_u = float(np.nanmean(var_classical))
    var_signal = max(var_assigned - var_u, 1e-12)
    lam = var_signal / (var_signal + var_u)
    log.info("classical-only reliability=%.3f (assigned var=%.4f classical var=%.4f)",
             lam, var_assigned, var_u)
    return {"reliability": lam, "var_assigned": var_assigned,
            "var_classical_mean": var_u}


def diagnose(df: pd.DataFrame) -> str:
    """Which investment would improve this exposure most?"""
    share = float(df["berkson_share"].median())
    if share > 0.65:
        msg = ("Berkson dominates: the areas are internally heterogeneous. Finer "
               "geography or address-level assignment would help; more monitors "
               "would not.")
    elif share < 0.35:
        msg = ("Classical dominates: the surface is poorly determined over these "
               "areas. More monitors or a better model would help; finer geography "
               "would not.")
    else:
        msg = ("Mixed: both the surface accuracy and the area size are limiting, "
               "and neither investment alone resolves the exposure error.")
    log.info("diagnosis: %s", msg)
    return msg

The key line is attenuation_from_classical, which uses only the classical variance in the reliability ratio. Including the Berkson component would over-correct, because Berkson error does not attenuate — it is precisely the mistake that follows from lumping the two together and applying a single correction. Studies that report implausibly large corrected effects for area-level exposures have frequently made exactly this error.

The diagnosis function exists because the decomposition is actionable. A public health agency asked to improve its exposure data wants to know whether to fund monitors or to move to finer geography, and the median Berkson share answers that directly — something neither a validation statistic nor a corrected effect estimate can do.

The Diagnosis Varies Across the Study Area Berkson share of total exposure error plotted against area size for the tracts in a study. Small dense urban tracts sit low on the Berkson share because they are internally homogeneous and lie near monitors, so their error is mostly classical. Large rural tracts sit high because they contain very different environments and the surface within them varies a lot. Two horizontal bands mark the thresholds where the diagnosis changes from classical-dominated through mixed to Berkson-dominated, and the pattern shows that one study can need different remedies in different places. One study can be surface-limited downtown and geography-limited in the country Berkson dominates: use finer areas classical dominates: add monitors mixed small urban large rural area size Berkson share

Validation & Edge Cases

No fine surface available. Without a resolution finer than the analysis areas there is no way to estimate within-area spread, and the Berkson component cannot be computed. Do not assume it is zero; that assumption converts the whole error into classical and over-corrects.

Population weights that are uniform. An unweighted within-area variance treats an uninhabited hillside as equal to a dense block. The Berkson component should reflect the spread across people, not across area.

Cells strongly correlated within an area. The classical component computed as the weighted mean of cell variances is conservative. Dividing by the cell count would assume independence and understate it severely, since neighbouring kriged cells are nearly perfectly correlated.

Areas with one cell. The Berkson component is zero by construction, which is an artefact of resolution rather than a homogeneous area. Flag those areas; their apparent homogeneity is not measured.

Differential error across areas. If the classical component is much larger in rural areas — which it usually is, because they are further from monitors — then attenuation differs across the study area. A single reliability ratio averages over that, and a sensitivity analysis restricted to well-monitored areas is worth reporting.

INFO exposure.decomp decomposition over 1842 areas: median classical=0.7421 median berkson=1.9068 median berkson share=0.72
INFO exposure.decomp classical-only reliability=0.881 (assigned var=6.2140 classical var=0.7421)
INFO exposure.decomp diagnosis: Berkson dominates: the areas are internally heterogeneous. Finer geography or address-level assignment would help; more monitors would not.

The over-correction that follows from lumping both components into one reliability ratio is large enough to change a study’s headline:

Including the Berkson Part Over-Corrects Three effect estimates from the same area-level analysis. The naive estimate is slightly attenuated by the classical component of the error. Correcting with the classical variance alone restores it modestly and lands near the truth. Correcting with the total variance, classical plus Berkson, treats the within-area spread as if it attenuated the slope, which it does not, and inflates the estimate far beyond the truth. The dashed reference line marks the true effect used in the simulation. Berkson error does not attenuate, so it must not be corrected for true effect naive classical only total variance inflated simulated area-level study with a known true effect

Compliance Notes

Record the fine surface resolution and its prediction variance source, the population weighting used for the within-area spread, the per-area classical and Berkson variances with their distributions, the classical-only reliability ratio, and the diagnosis with the thresholds that produced it.

Where a correction is applied, state explicitly that it used the classical component only and give both variances. A reviewer’s first question about a corrected area-level effect is whether the Berkson component was mistakenly included, and answering it in the methods is faster than answering it in review. Where the diagnosis is reported to programme staff rather than to a journal, the actionable sentence — monitors or finer geography — is the part that matters, and it should lead.