Spatial Epi Public Health GIS

Simulation Extrapolation for Exposure Misclassification

This guide, part of Exposure Measurement Error & Uncertainty, covers the correction to reach for when you know roughly how much error your exposure carries but have no gold-standard subset to calibrate against — which is the usual situation in environmental epidemiology.

Problem Context & Constraints

Regression calibration needs a validation sample where the true exposure is observed. Environmental studies rarely have one: personal monitoring campaigns are small, expensive and unrepresentative. Simulation extrapolation, usually called SIMEX, needs no such sample. It needs only an estimate of the error variance, which the exposure surface’s own cross-validation supplies.

The idea is a trick with a clear logic. If classical error attenuates the estimate, then adding more error should attenuate it further, in a way you can observe. So add a known multiple λ\lambda of the error variance, refit, and record the estimate. Do this for several values of λ\lambda and you trace out how the estimate decays as error grows. Then fit a smooth curve to that trajectory and extrapolate it back to λ=−1\lambda = -1, which corresponds to removing the error that was already there. The value at λ=−1\lambda = -1 is the SIMEX estimate.

Two constraints determine whether it works. The error variance must be roughly right. SIMEX corrects the amount of error you tell it about; supplying half the true variance under-corrects proportionally. The extrapolant must be appropriate. Extrapolating outside the observed range is inherently model-dependent, and the quadratic and nonlinear extrapolants in common use can disagree materially when the decay is strong.

A third consideration is specific to spatial exposure: the error is correlated in space. Adding independent noise in the simulation step assumes it is not, which makes SIMEX under-correct for the regional component. Where a covariance model exists, the conditional-simulation route in the sibling guide is preferable; SIMEX is the right tool when no covariance model is available.

The SIMEX Curve and Where the Answer Comes From A plot of the estimated effect against lambda, the multiple of extra error variance added. At lambda zero the naive estimate is plotted. As lambda increases to one and two the estimate decays steadily toward zero, and those points are computed by simulation. A fitted curve passes through them, and it is extrapolated leftwards past zero to lambda equal to minus one, where the value is read off as the corrected estimate. The extrapolated region is drawn with a dashed line and shaded to show that it is model-dependent, being outside the range where anything was observed. The answer is read where nothing was simulated — so the curve shape matters extrapolated naive estimate SIMEX estimate λ = −1 0 1 2 λ: multiple of the error variance added effect estimate Report the extrapolant used; a quadratic and a nonlinear fit can differ by a fifth

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, scipy 1.13.1, statsmodels 0.14.2
  • An estimate of the exposure error variance, per record or as a single value
  • A health model that refits quickly, since SIMEX fits it many times
  • A declared lambda grid and extrapolant, fixed before the run

Step-by-Step Solution

Add scaled noise at each lambda, average the estimate over repeated simulations at that lambda, fit the extrapolant and evaluate at minus one.

# SIMEX correction for a spatially assigned exposure with known error variance.
# 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
from scipy.optimize import curve_fit

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

LAMBDAS = (0.0, 0.5, 1.0, 1.5, 2.0)
SEED = 20260918


def _quadratic(lam, a, b, c):
    return a + b * lam + c * lam ** 2


def _nonlinear(lam, g0, g1, g2):
    """Rational extrapolant; usually the better-behaved choice for strong decay."""
    return g0 + g1 / (g2 + lam)


def simex(df: pd.DataFrame, outcome: str, exposure: str, error_var,
          covariates: list[str], lambdas=LAMBDAS, n_sim: int = 100,
          extrapolant: str = "nonlinear", seed: int = SEED) -> dict:
    """Simulation-extrapolation correction for classical measurement error."""
    rng = np.random.default_rng(seed)
    d = df.dropna(subset=[outcome, exposure] + covariates).reset_index(drop=True)
    z = d[exposure].to_numpy(float)
    sd = np.sqrt(np.clip(np.asarray(error_var, dtype=float) * np.ones(len(d)), 0, None))
    y = d[outcome].to_numpy(float)
    C = {c: d[c].to_numpy() for c in covariates}

    curve = []
    for lam in lambdas:
        reps = []
        n = 1 if lam == 0.0 else n_sim      # lambda 0 is the naive fit, no noise
        for _ in range(n):
            zz = z + (np.sqrt(lam) * sd * rng.standard_normal(len(d)) if lam > 0 else 0.0)
            X = sm.add_constant(pd.DataFrame({exposure: zz, **C}))
            reps.append(float(sm.OLS(y, X).fit().params[exposure]))
        curve.append({"lambda": float(lam), "beta": float(np.mean(reps)),
                      "sd_across_sims": float(np.std(reps, ddof=1)) if len(reps) > 1 else 0.0})
        log.info("lambda=%.2f mean_beta=%.5f (%d simulations)",
                 lam, curve[-1]["beta"], n)

    cv = pd.DataFrame(curve)
    lam_a = cv["lambda"].to_numpy()
    beta_a = cv["beta"].to_numpy()
    if extrapolant == "quadratic":
        popt, _ = curve_fit(_quadratic, lam_a, beta_a, maxfev=20000)
        corrected = float(_quadratic(-1.0, *popt))
    else:
        p0 = [beta_a[-1], (beta_a[0] - beta_a[-1]) * 2.0, 1.0]
        popt, _ = curve_fit(_nonlinear, lam_a, beta_a, p0=p0, maxfev=40000)
        corrected = float(_nonlinear(-1.0, *popt))

    naive = float(beta_a[0])
    log.info("SIMEX (%s): naive=%.5f corrected=%.5f ratio=%.3f",
             extrapolant, naive, corrected, corrected / naive if naive else np.nan)
    if not np.isfinite(corrected) or abs(corrected) > 10 * abs(naive):
        log.error("extrapolation is unstable — the decay curve does not support "
                  "reading a value at lambda = -1")
    return {"naive": naive, "corrected": corrected, "extrapolant": extrapolant,
            "curve": cv, "n": int(len(d))}

The instability guard matters because SIMEX’s failure mode is silent and dramatic. A nonlinear extrapolant fitted to a nearly flat decay curve can produce a pole near λ=−1\lambda = -1 and return a corrected estimate several times the naive one. That is not a large correction; it is a numerical artefact, and the honest response is to report that the correction could not be computed rather than to publish the number.

Running both extrapolants and reporting both is good practice. When they agree to within a few per cent, the extrapolation is well determined. When they differ by a fifth, the shape of the curve outside the observed range is doing the work, and the corrected estimate should be presented as a range rather than a point.

The Two Extrapolants Agree Where You Looked and Not Where You Read The same four simulated decay points fitted with a quadratic and with a rational extrapolant. Across the simulated range from lambda zero to two the two fitted curves are almost indistinguishable. Extrapolated back to lambda minus one they separate visibly, with the rational fit rising more steeply than the quadratic. The two corrected estimates are marked and their difference is annotated as the extrapolation uncertainty, which is not reflected in either fit's own standard error. Indistinguishable where simulated, divergent where it counts extrapolation uncertainty quadratic rational λ = −1 0 2 Report both corrected values; their spread is part of the result

Validation & Edge Cases

A decay curve that is flat. If the estimate barely moves as lambda rises, either the error variance supplied is far too small, or the exposure carries little error relative to its variance. Check the reliability ratio; a flat curve with a reliability near one is a correct result meaning no correction is needed.

A non-monotone curve. The estimate should decay steadily. Wobble means too few simulations per lambda; increase the count before fitting anything.

Error variance from the wrong source. Using the surface’s cross-validated error omits the residential-proxy component and under-corrects. Using an inflated guess over-corrects. Run SIMEX across a range of assumed variances and report the sensitivity.

Correlated error. SIMEX adds independent noise, so it addresses the independent component only. Where the surface error is regionally correlated — as kriging error always is — SIMEX under-corrects, and the conditional-simulation route is stronger.

Binary or categorical exposures. SIMEX extends to misclassification of a categorical exposure, but the simulation step is a reclassification driven by a misclassification matrix rather than additive noise, and that matrix has to come from somewhere. Applying the continuous version to a dichotomised exposure — adding Gaussian noise to a zero-one variable and rounding — is not the same procedure and has no supporting theory. Dichotomise after correcting, not before.

Covariates measured with error too. The correction addresses one variable. If a confounder is also mismeasured, correcting the exposure alone can move the estimate further from the truth than leaving both uncorrected, because residual confounding is now scaled up along with the exposure effect. Where a key confounder is a modelled quantity rather than a measured one, say so and treat the corrected estimate as an upper bound on what the correction achieves.

Standard errors. The naive standard errors from the final fit do not reflect the extrapolation. Use a bootstrap over the whole SIMEX procedure, or report the estimate with an explicit note that its interval is approximate.

INFO exposure.simex lambda=0.00 mean_beta=0.01120 (1 simulations)
INFO exposure.simex lambda=1.00 mean_beta=0.00842 (100 simulations)
INFO exposure.simex lambda=2.00 mean_beta=0.00692 (100 simulations)
INFO exposure.simex SIMEX (nonlinear): naive=0.01120 corrected=0.01684 ratio=1.504

Because the correction depends directly on the error variance supplied, the sensitivity of the answer to that input should be shown rather than assumed:

The Correction Scales With the Error You Assume The SIMEX-corrected effect estimate plotted against the error variance supplied to it, expressed as a multiple of the cross-validated surface error. At zero assumed error the correction returns the naive estimate. As the assumed error grows the corrected estimate rises, slowly at first and then steeply. A shaded band marks the plausible range between the surface's own cross-validated error and that error plus an allowance for the residential proxy, and the corrected estimate across that band is the range to report. Report the corrected estimate across the plausible error range plausible error range naive estimate report this span 0 surface CV error plus proxy allowance assumed exposure error variance corrected estimate

Compliance Notes

Record the error variance used and its derivation, the lambda grid, the number of simulations per lambda, the seed, the extrapolant and the corrected estimate under each extrapolant tried, and the decay curve itself so a reviewer can see whether the extrapolation was well determined.

State in the methods that SIMEX assumes independent classical error, and say what is known about the spatial correlation of the exposure error. That sentence is what separates a defensible use of the method from an application of a correction whose assumptions were never checked — and where the correlated component is large, it is also the sentence that tells a reader the correction is a lower bound.