Spatial Epi Public Health GIS

Propagating Kriging Variance Into Health Effect Estimates

This guide, part of Exposure Measurement Error & Uncertainty, takes the one input family that hands you an uncertainty for free — geostatistical prediction — and carries that uncertainty all the way into the health model instead of dropping it at the join.

Problem Context & Constraints

Kriging produces Z^(s)\hat{Z}(s) and σ2(s)\sigma^2(s) at every location. The second quantity is discarded in the overwhelming majority of applied studies, which then report a confidence interval that describes sampling variation in the health outcome and nothing at all about how well the exposure was known. In a study area where half the population lives beyond the monitor network, that omission is not a rounding detail.

The naive fix is to add independent noise to each person’s exposure, drawn with their own prediction standard deviation, and repeat. It is better than nothing and it understates the uncertainty, for a specific reason: kriging errors are spatially correlated. If the surface is too high in one county, it is too high for everyone in that county at once. Independent per-point draws let those errors cancel across the cohort, producing a between-draw variance far smaller than the truth. The effect estimate, which depends on contrasts across the study area, is sensitive to exactly the correlated component that independent draws destroy.

The correct instrument is conditional simulation: generate whole realisations of the exposure field that honour the observed monitor values and reproduce the fitted covariance. Each realisation is a plausible version of reality, complete with its spatial structure. Assign exposures from each realisation, fit the health model, and pool.

Two constraints follow. Conditional simulation is expensive on a fine grid, so it is usually run at the resolution the study actually uses rather than the surface’s native one. And it requires a fitted covariance model, which means it is available for kriging and not for methods that return no variance.

Independent Noise Destroys the Structure That Matters Two panels showing exposure error across a study area. In the independent-noise panel each location receives its own unrelated perturbation, so the errors form a speckled pattern that cancels when averaged across any group of people. In the conditional-simulation panel the error is spatially coherent, with whole regions shifted up or down together, because the realisation honours the fitted covariance. A note records that the health effect estimate depends on contrasts between regions, which the coherent error moves and the speckled error does not. Kriging error is regional, not speckled — and the estimate cares about regions Independent per-point noise errors cancel across any group One conditional realisation too high here too low whole regions move together Independent draws give a reassuringly narrow interval that is simply wrong Each realisation must honour the monitors and reproduce the fitted covariance

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, scipy 1.13.1, statsmodels 0.14.2
  • A fitted variogram with nugget, sill and range from Interpolating Environmental Exposure Surfaces
  • Monitor sites with values, and the analysis locations — usually area centroids or population-weighted points
  • A health model that can be refitted cheaply, since it will run once per realisation

Step-by-Step Solution

Generate unconditional realisations from the covariance model, condition them on the observed monitors by kriging the simulated residuals, assign, refit and pool.

# Conditional simulation of a kriged exposure field, propagated into a health model.
# Pinned: numpy==1.26.4, pandas==2.2.2, scipy==1.13.1, statsmodels==0.14.2
import logging

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

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

SEED = 20260918


def exponential_cov(h: np.ndarray, nugget: float, sill: float, rng_m: float) -> np.ndarray:
    """Covariance implied by the fitted exponential semivariogram."""
    partial = sill - nugget
    c = partial * np.exp(-h / rng_m)
    c = np.where(h == 0.0, sill, c)
    return c


def conditional_realisations(sites_xy: np.ndarray, sites_z: np.ndarray,
                             target_xy: np.ndarray, krige_fn,
                             nugget: float, sill: float, rng_m: float,
                             n_real: int = 100, seed: int = SEED) -> np.ndarray:
    """Return an (n_real, n_target) array of plausible exposure fields.

    Method: simulate an UNCONDITIONAL field with the right covariance at both the
    sites and the targets, krige the simulated site values, and add the difference
    between the simulated field and its own kriging to the real kriged surface.
    The result honours the observed monitors exactly and carries the fitted
    spatial structure — which independent per-point noise does not.
    """
    rng = np.random.default_rng(seed)
    all_xy = np.vstack([sites_xy, target_xy])
    n_s = len(sites_xy)
    d = np.sqrt(((all_xy[:, None, :] - all_xy[None, :, :]) ** 2).sum(-1))
    C = exponential_cov(d, nugget, sill, rng_m)
    C[np.diag_indices_from(C)] += 1e-8 * sill          # numerical jitter
    L = np.linalg.cholesky(C)

    base, _ = krige_fn(sites_xy, sites_z, target_xy)    # the kriged surface itself
    out = np.empty((n_real, len(target_xy)), dtype=float)
    for r in range(n_real):
        field = L @ rng.standard_normal(len(all_xy))
        sim_at_sites, sim_at_targets = field[:n_s], field[n_s:]
        sim_krige, _ = krige_fn(sites_xy, sim_at_sites, target_xy)
        out[r] = base + (sim_at_targets - sim_krige)
    log.info("conditional simulation: %d realisations, %d targets, "
             "mean between-realisation sd=%.4f",
             n_real, len(target_xy), float(out.std(axis=0).mean()))
    return out


def pool_over_realisations(realisations: np.ndarray, health: pd.DataFrame,
                           target_index: np.ndarray, outcome: str,
                           covariates: list[str]) -> dict:
    """Fit the health model once per realisation and pool with Rubin's rules."""
    betas, vars_ = [], []
    for r in range(realisations.shape[0]):
        z = realisations[r][target_index]
        X = sm.add_constant(pd.DataFrame({"exposure": z,
                                          **{c: health[c].to_numpy() for c in covariates}}))
        fit = sm.OLS(health[outcome].to_numpy(float), X).fit()
        betas.append(float(fit.params["exposure"]))
        vars_.append(float(fit.bse["exposure"]) ** 2)
    m = realisations.shape[0]
    b = float(np.mean(betas))
    within = float(np.mean(vars_))
    between = float(np.var(betas, ddof=1))
    total = within + (1.0 + 1.0 / m) * between
    frac_exposure = (1.0 + 1.0 / m) * between / total
    log.info("pooled beta=%.5f se=%.5f; %.0f%% of the variance comes from "
             "exposure uncertainty", b, np.sqrt(total), 100.0 * frac_exposure)
    return {"beta": b, "se": float(np.sqrt(total)), "within_var": within,
            "between_var": between, "exposure_variance_share": frac_exposure,
            "n_realisations": m}

The conditioning step — adding the difference between a simulated field and the kriging of its own site values — is what makes each realisation pass exactly through the observed monitors while varying freely in between. Without it, the realisations would be unconstrained fields that ignore the data entirely.

The reported share of variance coming from exposure uncertainty is the number to lead with. A study where 5% of the standard error is exposure uncertainty can reasonably present the naive interval with a footnote; one where it is 40% cannot, and the difference is visible only if the propagation is done.

How Much the Interval Actually Widens Two distributions of the effect estimate across draws. Under independent per-point noise the estimates cluster tightly around the naive value, because the added errors cancel across the cohort and the contrasts barely move. Under conditional simulation the distribution is much wider, because each realisation shifts whole regions together and therefore changes the exposure contrast the model is fitted to. The naive point estimate is marked and sits near the centre of both, showing that the correction is about the interval rather than the estimate. Both are centred correctly; only one reports the right uncertainty naive estimate independent noise conditional simulation effect estimate across draws The wide distribution is the honest one, and it is the one that decides whether the finding holds

Validation & Edge Cases

Realisations that do not reproduce the variogram. Compute the empirical variogram of a few realisations and compare it with the fitted model. A mismatch means the covariance matrix or the conditioning step is wrong, and everything downstream inherits it.

Cholesky failure on a large target set. The covariance matrix is dense and its factorisation is cubic. Simulate on a coarser grid and sample the analysis locations from it, or use a sequential or spectral simulation method; do not silently drop to independent noise.

Realisations at the wrong resolution. Simulating at 1 km and assigning to addresses assumes exposure is constant within the cell, which reintroduces a Berkson component that the simulation does not represent. State the simulation resolution.

Too few realisations. The between-realisation variance is estimated from the realisation count. Below about fifty it is itself noisy; a hundred is a comfortable default for a study that will be published.

A health model too slow to refit a hundred times. Reduce the realisation count and report it, or fit on a summarised design. Do not substitute a single “mean plus variance” adjustment and describe it as propagation.

INFO exposure.condsim conditional simulation: 100 realisations, 4218 targets, mean between-realisation sd=1.3402
INFO exposure.condsim pooled beta=0.01184 se=0.00412; 38% of the variance comes from exposure uncertainty

The verification step in that list compares each realisation’s own variogram with the one it was simulated from, and a correct run looks like this:

Realisations Should Reproduce the Model The fitted exponential variogram drawn as a solid curve, with the empirical variograms of five conditional realisations drawn as thin lines around it. The realisations scatter on either side of the model and follow its shape, rising from the nugget and levelling at the sill near the fitted range. A realisation whose empirical variogram sat consistently above or below the model, or failed to level off, would indicate an error in the covariance matrix or the conditioning step. Scatter around the model is correct; a systematic offset is a bug fitted model five realisations separation distance semivariance check a handful before trusting a hundred

Compliance Notes

Record the variogram parameters used for simulation, the simulation method, the resolution and target set, the realisation count and the seed, the health model specification, and the pooled estimate with its within- and between-realisation variance components and the exposure share of the total.

Report the exposure variance share wherever the effect estimate appears. It is the single most useful number this procedure produces: it tells a reader how much of the reported uncertainty comes from not knowing the exposure, which is information no conventionally reported interval contains. Where the share is large, it also tells the agency something actionable — that the cheapest route to a more precise finding is more monitors, not more cases.