Spatial Epi Public Health GIS

Choosing Distance-Decay Functions for 2SFCA

The enhanced two-step floating catchment area score computed in the spatial equity index calculation section is acutely sensitive to the distance-decay function and its parameters, so the wrong decay can manufacture or hide access deserts. This guide implements Gaussian, gravity, and piecewise decay in Python, sweeps the accessibility index across them, and quantifies how stably the most-underserved units keep their rank.

Problem Context & Constraints

The E2SFCA accessibility score weights each reachable facility by a decay function of travel cost: a clinic across the street counts more than one at the edge of the catchment. That weighting is not a cosmetic refinement — it is the mechanism that turns a binary “reachable within the threshold” into a graded accessibility surface. Because the same decay kernel discounts competitors in the supply step and facilities in the demand step, its shape propagates into every score twice. A decay that falls off too slowly flattens the surface and washes out genuine disparities; one that falls off too sharply concentrates all the weight on the nearest facility and can invent a desert wherever the closest clinic is just past a threshold.

The step-by-step mechanics of the two floating catchments — how supply-to-demand ratios are built and summed back onto residents — are derived in calculating the two-step floating catchment area index; this page holds those mechanics fixed and varies only the decay function inside them. The problem is that practitioners routinely pick a decay family by convention — Gaussian because a paper used it, a step function because it is simple — and then report a single accessibility ranking as if it were the decay-independent truth. It is not. Two defensible decay choices can disagree about which tract is the most underserved, and if an allocation decision rides on that ranking, the choice of function is a policy decision disguised as a technical default.

Three constraints frame a disciplined approach.

  • The decay family must match care-seeking behavior for the service. Emergency care is near-binary inside a response-time guarantee (a piecewise cliff); routine primary care declines smoothly with travel (Gaussian); specialty referral tolerates long trips with a heavy tail (gravity/power). Picking a family that contradicts observed behavior biases every score.
  • Parameters must be justified, not defaulted. A Gaussian bandwidth set to half the catchment threshold is a convenience, not a calibration. The parameter that controls how fast access falls off should track the travel cost at which care-seeking is actually observed to drop.
  • Sensitivity must be reported, not assumed away. Because the score has no analytic standard error, defensibility comes from showing how the ranking of underserved units holds up across plausible decay choices — not from asserting one is correct.

The three decay families weight travel cost very differently across the catchment, which is exactly why the resulting accessibility ranking can shift:

Distance-Decay Weight Curves by Function Family A line chart. The horizontal axis is travel time from zero to the catchment threshold d-zero; the vertical axis is the decay weight W from zero to one. Three curves start near weight one at zero travel time. The piecewise step holds near one then drops in stages to zero at the threshold. The Gaussian falls slowly at first then steeply through the middle distances. The gravity power-law drops sharply near the origin then keeps a long thin tail toward the threshold. The three curves cross, showing they rank mid-distance facilities differently. W=1 W=0 decay weight travel time d₀ piecewise Gaussian gravity The families weight mid-distance facilities differently, so the underserved ranking can depend on which is chosen

Prerequisites

Pin the stack so the sweep is reproducible:

  • python 3.11
  • numpy 1.26.4
  • geopandas 0.14.4 for the population and score layers
  • scipy 1.13.1 for the rank-correlation diagnostic

Input state assumed by the code below:

  • A travel-cost vector or matrix in consistent units (minutes) between population units and facilities, from an upstream routed matrix or drive-time isochrone generation.
  • Supply and demand attributes already harmonized to one projected CRS, as the parent section requires.
  • A catchment threshold d0d_0 in the same units, beyond which the weight is zero.

Step-by-Step Solution

Implement the three decay families as pure functions of travel cost dd so they are interchangeable inside the E2SFCA machinery. The Gaussian weight uses a bandwidth β\beta that controls how fast access falls off:

Wgauss(d)=exp⁡ ⁣(−d22β2),d≤d0W_{\text{gauss}}(d) = \exp\!\left(-\frac{d^{2}}{2\beta^{2}}\right), \qquad d \le d_0

The gravity (power-law) weight uses a friction exponent γ\gamma and keeps a heavy tail for long trips:

Wgrav(d)=d−γ,d≥dmin⁡,  d≤d0W_{\text{grav}}(d) = d^{-\gamma}, \qquad d \ge d_{\min}, \; d \le d_0

where dmin⁡d_{\min} is a small floor that prevents the weight from diverging as d→0d \to 0. The piecewise weight assigns a constant within each of several travel-time bands, encoding a stepped guarantee.

# python 3.11
# numpy==1.26.4  scipy==1.13.1
import logging
import numpy as np

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

def gaussian_decay(d, d0, beta):
    """Smooth decay for routine care; beta sets how fast access falls off (minutes)."""
    w = np.exp(-(d ** 2) / (2.0 * beta ** 2))
    return np.where(d <= d0, w, 0.0)          # hard zero beyond the catchment

def gravity_decay(d, d0, gamma, d_min=1.0):
    """Heavy-tailed decay for specialty referral; d_min floors the singularity at d->0."""
    dd = np.clip(d, d_min, None)
    w = dd ** (-gamma)
    return np.where(d <= d0, w / w.max(), 0.0)   # normalize peak to 1 for comparability

def piecewise_decay(d, bands=((10, 1.0), (20, 0.68), (30, 0.22))):
    """Stepped guarantee: constant weight inside each travel-time band, else 0."""
    w = np.zeros_like(np.asarray(d, dtype=float))
    prev = 0.0
    for upper, weight in bands:
        w = np.where((d > prev) & (d <= upper), weight, w)
        prev = upper
    return w

With interchangeable decay functions, the sensitivity sweep runs the accessibility index once per decay configuration and stacks the results. The E2SFCA scoring itself is unchanged from the parent section; only the decay_fn argument varies.

def sweep_decay(compute_index, pop_gdf, fac_gdf, configs):
    """
    compute_index(pop_gdf, fac_gdf, decay_fn) -> np.ndarray of per-unit accessibility.
    configs: list of (label, decay_fn). Returns (labels, score_stack) with stack
    shape (n_configs, n_units), rows in the order of `configs`.
    """
    labels, rows = [], []
    for label, decay_fn in configs:
        scores = compute_index(pop_gdf, fac_gdf, decay_fn=decay_fn)
        labels.append(label)
        rows.append(scores)
        log.info("swept decay=%s min=%.4f median=%.4f", label,
                 float(np.nanmin(scores)), float(np.nanmedian(scores)))
    return labels, np.vstack(rows)

Validation & Edge Cases

A decay choice is only defensible if you can show what it does and does not change. Three diagnostics do that.

1. Rank correlation across decay choices. The operational question is rarely the absolute score — it is which units are worst off. Compute Spearman’s rank correlation between every pair of decay configurations; a high correlation means the ranking is robust to the decay choice, a low one means the ranking is decay-driven and any allocation must acknowledge that:

from scipy.stats import spearmanr

def rank_stability(labels, score_stack):
    """Pairwise Spearman rho between decay configurations. Low rho = fragile ranking."""
    n = len(labels)
    for a in range(n):
        for b in range(a + 1, n):
            rho, _ = spearmanr(score_stack[a], score_stack[b])
            log.info("rank_corr %s vs %s rho=%.3f", labels[a], labels[b], rho)
INFO rank_corr gaussian_b8 vs gravity_g1.5 rho=0.71
INFO rank_corr gaussian_b8 vs piecewise    rho=0.94

A ρ\rho of 0.71 between Gaussian and gravity is a warning: roughly a third of the ranking reorders, so the most-underserved list is partly an artifact of the decay family.

2. Rank stability of the most-underserved units specifically. Aggregate correlation can hide churn at the tail that is exactly what matters. Track how many of the bottom-kk units stay in the bottom-kk across configurations:

def bottom_k_overlap(score_stack, k=20):
    """Jaccard overlap of the k lowest-scoring units across all decay configs."""
    bottoms = [set(np.argsort(row)[:k]) for row in score_stack]
    common = set.intersection(*bottoms)
    union = set.union(*bottoms)
    log.info("bottom_%d stable=%d/%d jaccard=%.2f", k, len(common), len(union),
             len(common) / len(union))
    return common

Only the units in common are unambiguously underserved regardless of decay; those are the defensible priorities.

3. Boundary behavior at d→0d \to 0 and d→d0d \to d_0. The gravity weight diverges as travel cost approaches zero unless floored, and any family must go cleanly to zero at the threshold rather than leaving a discontinuous crumb of weight just inside it. Assert both:

d = np.array([0.0, 0.5, 15.0, 29.9, 30.0, 30.1])
assert np.all(np.isfinite(gravity_decay(d, d0=30, gamma=1.5))), "gravity must not diverge at d=0"
assert gaussian_decay(np.array([30.1]), d0=30, beta=8)[0] == 0.0, "weight must vanish beyond d0"

Compliance Notes

Three records make a decay-based accessibility surface defensible:

  • Log the decay family and every parameter. The chosen family, the Gaussian bandwidth or gravity exponent, the piecewise band edges, and the catchment threshold must all be recorded with the output. Because the decay is effectively a policy choice, an unlogged decay makes the resulting shortage designation impossible to reproduce or contest.
  • Publish the sensitivity result, not just the point estimate. Persist the pairwise rank correlations and the bottom-kk overlap alongside the score. Reporting a single underserved ranking without its rank stability overstates certainty; the overlap set is the honest, decay-robust finding.
  • Deterministic, version-pinned sweep. Fix the decay configuration list, sort units by stable ID, and record the numpy and scipy versions, so a reviewer re-running the sweep reproduces the same scores and the same rank-stability numbers for peer review and interagency handoff.