Spatial Epi Public Health GIS

Bootstrapping Exposure Surfaces for Confidence Intervals

This guide, part of Exposure Measurement Error & Uncertainty, provides the fallback for exposure models that return no prediction variance — land-use regression, inverse distance weighting, machine-learned surfaces — by generating the uncertainty from resampling instead of from a covariance model.

Problem Context & Constraints

Kriging hands you a variance because it assumes a covariance model. Land-use regression does not, and neither do the ensemble methods increasingly used for exposure surfaces. That leaves a study with a point prediction per location and nothing to say about how well determined it is, which blocks every correction in the parent guide.

The bootstrap replaces the missing model with resampling. Draw a resample of the monitoring sites, rebuild the entire surface from that resample, assign exposures, fit the health model, and repeat. The spread of the resulting effect estimates includes both the health model’s sampling variation and the uncertainty in the exposure surface, because the surface was rebuilt every time.

Three constraints govern whether the result means anything. The whole pipeline must be inside the loop. Resampling monitors but reusing the variable selection, the fitted variogram or the chosen exponent from the full data leaks information and understates the spread — the same leakage that makes flat cross-validation optimistic. The resample must respect spatial structure. Independent resampling of sites breaks the spatial dependence between them; a block bootstrap that resamples spatial clusters of sites preserves it and gives a wider, more honest interval. The monitor network is the population being resampled, so the bootstrap describes uncertainty due to which sites were measured, not uncertainty due to the model form being wrong.

Resample Blocks, Not Sites Two resampling schemes over the same monitor network. Site-level resampling draws individual monitors with replacement, so a duplicated site sits beside its own copy and the spatial dependence between neighbouring sites is broken, producing replicates that are more alike than they should be. Block resampling divides the study area into spatial blocks and draws whole blocks with replacement, so an entire region is either over-represented or absent from a replicate. The block scheme produces a wider spread of surfaces and therefore a wider, more honest interval. Independent site resampling produces replicates that are too similar Site-level resample every region still represented in every replicate Block resample block absent a whole region can be missing, as in reality The block bootstrap asks: what if this region had never been monitored? which is the uncertainty the exposure surface actually carries

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4, statsmodels 0.14.2
  • A surface-building function that takes monitor sites and returns predictions at targets, performing any selection or tuning internally
  • Monitor sites with coordinates in a projected metric system
  • A block size, chosen from the spatial scale of dependence in the network

Step-by-Step Solution

Resample blocks of sites, rebuild the surface from scratch inside each replicate, and pool the effect estimates into a percentile interval.

# Block bootstrap of an exposure surface, propagated into a health effect estimate.
# Pinned: numpy==1.26.4, pandas==2.2.2, geopandas==0.14.4, statsmodels==0.14.2
import logging

import numpy as np
import pandas as pd
import geopandas as gpd
import statsmodels.api as sm

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

SEED = 20260918


def spatial_blocks(sites: gpd.GeoDataFrame, block_m: float) -> np.ndarray:
    """Assign each site to a square block of the given size."""
    x = sites.geometry.x.to_numpy()
    y = sites.geometry.y.to_numpy()
    bx = np.floor((x - x.min()) / block_m).astype(int)
    by = np.floor((y - y.min()) / block_m).astype(int)
    return bx * 100_000 + by


def block_bootstrap(sites: gpd.GeoDataFrame, value_col: str,
                    build_surface, health: pd.DataFrame, outcome: str,
                    covariates: list[str], target_geom: gpd.GeoDataFrame,
                    block_m: float, n_boot: int = 500, seed: int = SEED) -> dict:
    """Rebuild the whole surface in every replicate and collect effect estimates.

    build_surface(train_sites, value_col, targets) -> array of predictions.
    It MUST perform its own variable selection / tuning internally, or the
    interval will be too narrow.
    """
    rng = np.random.default_rng(seed)
    sites = sites.sort_values("site_id").reset_index(drop=True)
    blocks = spatial_blocks(sites, block_m)
    uniq = np.unique(blocks)
    log.info("bootstrap over %d sites in %d blocks of %.0f m", len(sites), len(uniq), block_m)

    betas, failures = [], 0
    for b in range(n_boot):
        pick = rng.choice(uniq, size=len(uniq), replace=True)
        idx = np.concatenate([np.flatnonzero(blocks == p) for p in pick])
        train = sites.iloc[idx].reset_index(drop=True)
        if train[value_col].nunique() < 5:
            failures += 1
            continue
        try:
            z = build_surface(train, value_col, target_geom)
        except Exception as exc:                # a replicate can be unfittable
            failures += 1
            log.debug("replicate %d failed: %s", b, exc)
            continue
        X = sm.add_constant(pd.DataFrame({"exposure": np.asarray(z, dtype=float),
                                          **{c: health[c].to_numpy() for c in covariates}}))
        betas.append(float(sm.OLS(health[outcome].to_numpy(float), X).fit().params["exposure"]))

    betas = np.asarray(betas)
    if len(betas) < 0.8 * n_boot:
        log.warning("%d/%d replicates failed — the interval rests on a biased subset "
                    "of resamples", failures, n_boot)
    out = {"beta_mean": float(betas.mean()), "beta_sd": float(betas.std(ddof=1)),
           "ci_low": float(np.percentile(betas, 2.5)),
           "ci_high": float(np.percentile(betas, 97.5)),
           "n_ok": int(len(betas)), "n_failed": failures, "block_m": block_m}
    log.info("block bootstrap: beta=%.5f 95%% CI [%.5f, %.5f] from %d replicates",
             out["beta_mean"], out["ci_low"], out["ci_high"], out["n_ok"])
    return out

The build_surface contract is the whole design. It must be a function that goes from a set of monitors to predictions, performing every data-dependent choice inside itself: the variable selection for a land-use regression, the exponent sweep for inverse distance weighting, the variogram fit for kriging. Passing a pre-selected model into the loop produces an interval that measures only the noise in the coefficients of a fixed specification, which is the smallest and least interesting part of the uncertainty.

Counting failed replicates matters because failures are not random. A resample missing the only sites in a region will often be unfittable, and those are exactly the replicates that would have produced the most extreme estimates. An interval computed from the survivors is narrower than the truth, and the failure count is what lets a reader judge by how much.

What Each Bootstrap Design Includes Four intervals for the same effect estimate. The naive interval from a fixed exposure column is narrowest and includes only the health model's sampling variation. A site-level bootstrap that reuses the fitted model is slightly wider, adding coefficient noise. A site-level bootstrap that refits the whole pipeline is wider again, adding selection and tuning uncertainty. A block bootstrap that refits the whole pipeline is widest, additionally allowing whole regions to be unmonitored. Each interval is annotated with the source of uncertainty it newly includes. Each design adds one more source of uncertainty that was previously hidden naive outcome sampling only site boot, fixed model plus coefficient noise site boot, full refit plus selection and tuning block boot, full refit plus unmonitored regions Report which design produced the interval — the label changes what it means None of them covers the possibility that the model form itself is wrong

Validation & Edge Cases

Block size chosen arbitrarily. Blocks smaller than the spatial dependence in the network behave like site resampling. Size them from the fitted range where one exists, or from the typical distance over which monitor values correlate.

Too few blocks. A network in six blocks gives at most a coarse resampling distribution and a lumpy interval. Below about fifteen blocks, report the interval as approximate and say how many there were.

Duplicated sites inside a replicate. Resampling with replacement puts the same site in twice, which makes a kriging system singular. The surface builder must handle duplicates — by jittering, averaging or using a least-squares solve — rather than crashing every other replicate.

Reusing the health model’s own data. The bootstrap here resamples monitors, not people. Uncertainty in the health outcome is captured by each replicate’s own standard errors, and the percentile interval across replicates captures exposure uncertainty. Resampling both requires a nested design and considerably more computation.

Computational cost that forces shortcuts. Five hundred replicates of a full land-use regression refit plus a health model can run for hours. The tempting economy is to rebuild the surface less often than the health model, reusing one surface across several replicates. That halves the runtime and removes most of what the bootstrap was measuring. If the budget is genuinely binding, reduce the replicate count and report it rather than reducing what each replicate does.

Intervals that are asymmetric. Percentile intervals from a bootstrap of a nonlinear pipeline are frequently asymmetric around the point estimate, and that asymmetry is information rather than an artefact to be tidied away. Report the percentile bounds as they come out; forcing a symmetric interval by doubling a standard deviation discards the shape the resampling revealed.

Confusing this interval with a model-form interval. The bootstrap describes uncertainty given the chosen model family. It cannot tell you that land-use regression was the wrong choice. That question needs a comparison across families, reported separately.

INFO exposure.boot bootstrap over 58 sites in 22 blocks of 60000 m
WARNING exposure.boot 71/500 replicates failed — the interval rests on a biased subset of resamples
INFO exposure.boot block bootstrap: beta=0.01042 95% CI [0.00188, 0.01903] from 429 replicates

The failure count earns its place in the report because failed replicates are not a random subset of the resamples:

Failed Replicates Come From the Tails A distribution of effect estimates across bootstrap replicates. The successfully fitted replicates form the central body. Replicates that could not be fitted are shown as a separate shaded region concentrated at both extremes, because a resample missing every site in a region is both the most likely to fail and the most likely to produce an extreme estimate. An interval computed from the successful replicates alone is therefore narrower than the true bootstrap interval. Dropping failures quietly trims the tails of the interval fitted replicates failed failed interval from survivors only effect estimate across replicates report the failure count with the interval

Compliance Notes

Record the resampling unit and block size with the reasoning behind it, the replicate count and seed, the surface-building function and confirmation that all data-dependent steps ran inside the loop, the number of failed replicates with their cause, and the resulting interval alongside the naive one.

State the design in the methods with enough specificity that a reader can tell which sources of uncertainty the interval includes. “Bootstrapped confidence interval” is not sufficient — a site-level bootstrap around a fixed model and a block bootstrap around a full refit differ by a factor that routinely exceeds two, and they answer different questions. Where replicates failed, report the count; an interval computed from the resamples that happened to be fittable is not the interval the method promises.