Geographically Weighted Regression for Disease-Rate Modeling

A single global regression coefficient assumes a covariate’s effect on disease rates is identical everywhere; geographically weighted regression relaxes that by fitting a local model at every spatial unit, and it extends the regression toolkit in Disease Clustering & Spatial Statistical Modeling. GWR maps where a driver such as poverty, age structure, or air quality actually moves incidence — and, critically for public health, where its sign flips between regions that a global ordinary least squares (OLS) model would average into a single misleading number.

The problem GWR solves is spatial non-stationarity: the assumption that one relationship holds across the whole study area is often false in epidemiology. Poverty may raise respiratory-disease rates sharply in a dense industrial corridor, have almost no measurable effect in an affluent suburban belt, and be confounded by altitude or humidity in a rural interior. A global model reports the population-weighted average of those distinct regimes — a coefficient that describes nowhere in particular and can be significant, null, or wrong-signed depending only on how the study area was drawn.

Concept & Epidemiological Alignment

GWR is a local regression. Instead of estimating one coefficient vector β\beta for the entire dataset, it estimates a separate coefficient vector at each location (u,v)(u, v) by re-fitting a weighted least squares model in which nearby observations count more than distant ones. The estimator at location (u,v)(u, v) is:

β^(u,v)=(XW(u,v)X)1XW(u,v)y\hat{\beta}(u, v) = \left(X^{\top} W(u, v)\, X\right)^{-1} X^{\top} W(u, v)\, y

where XX is the design matrix of covariates, yy is the disease-rate response, and W(u,v)W(u, v) is a diagonal matrix of geographic weights that depends on the location being fitted. Every unit gets its own W(u,v)W(u, v), so the model produces one coefficient surface per covariate rather than one number. Where a spatial-lag or spatial-error model (covered under Spatial Lag & Error Regression) keeps the coefficients global and instead absorbs spatial dependence through an autoregressive term, GWR keeps the errors simple and lets the coefficients themselves carry the spatial structure. The two are complementary lenses on the same residual problem, not substitutes.

The weights come from a distance-decay kernel centered on the fit point. The two standard choices are the Gaussian kernel, which is continuous and never reaches zero,

wij=exp ⁣(12(dijb)2),w_{ij} = \exp\!\left(-\tfrac{1}{2}\left(\frac{d_{ij}}{b}\right)^{2}\right),

and the bisquare kernel, which is compact — it assigns zero weight beyond the bandwidth bb and so bounds each local fit to a finite neighborhood,

wij={[1(dijb)2]2,dijb0,dij>bw_{ij} = \begin{cases}\left[1 - \left(\dfrac{d_{ij}}{b}\right)^{2}\right]^{2}, & d_{ij} \le b \\[2mm] 0, & d_{ij} > b\end{cases}

Here dijd_{ij} is the distance from fit location ii to observation jj and bb is the bandwidth. Bandwidth is the single most consequential parameter: too small and each local fit sees too few observations, producing noisy, overfit coefficients; too large and GWR collapses back toward global OLS. The choice between a fixed-distance bb and an adaptive nearest-neighbor bb is important enough to have its own guide — see Choosing Adaptive vs Fixed Bandwidth for GWR.

The contrast between what a global model reports and what GWR reveals is the whole reason to reach for it:

Global OLS Coefficient vs GWR Coefficient Surface Two nine-cell map panels for the same covariate. The left panel, global OLS, shows the identical coefficient value 0.42 in every cell — one relationship imposed everywhere. The right panel, GWR, shows the coefficient varying cell by cell from strongly positive in the top-left through near zero in the centre to negative in the lower-left, so the sign of the covariate effect reverses across the study area. A labelled arrow between the panels marks that both describe the same covariate. Global OLS one coefficient for every unit 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 GWR coefficient surface B(u,v) 0.61 0.44 0.12 0.30 0.05 0.28 -0.19 0.08 -0.31 same covariate shaded cells mark where the local effect weakens or reverses sign GWR replaces one imposed coefficient with a per-unit surface that exposes regional reversals

Reach for GWR when the scientific question is explicitly where does this relationship hold, and where does it break, and when you have enough units (roughly 150+) to support a stable local fit. Prefer a global model when the relationship is genuinely constant, when the unit count is small, or when a single interpretable coefficient is what a policy audience needs. GWR is descriptive and exploratory first; it is not a substitute for a causal design, and a coefficient surface is a hypothesis generator, not proof of a local mechanism.

Method-Selection Table

Method Coefficients Handles spatial dependence via Use when Main risk
OLS Global, one per covariate Nothing — assumes independent errors Relationship is stationary; small n; you need one headline number Averages away regional reversals; autocorrelated residuals invalidate inference
Spatial lag / error Global An autoregressive term on yy or the error Dependence is a nuisance you want to model out, coefficients still assumed constant Still reports one coefficient per covariate
GWR Local, one surface per covariate, single bandwidth Letting coefficients vary over space Non-stationarity is the question; all covariates vary at a similar scale One bandwidth forces every covariate to the same spatial scale; overfitting
MGWR Local, one surface per covariate, covariate-specific bandwidth Per-covariate spatial scale Some drivers act regionally and others locally Slower to fit; heavier configuration and audit surface

The progression down the table trades interpretability and speed for spatial flexibility. Multiscale GWR (MGWR) is the natural next step once you accept non-stationarity but reject the assumption — baked into plain GWR — that every covariate operates at one shared bandwidth; population age structure may vary smoothly at a regional scale while a point-source pollutant acts within a few kilometres.

Spatial Data Prerequisites

GWR is built on distances, so its correctness depends entirely on a metric coordinate system. Running it on unprojected WGS84 degrees makes the kernel bandwidth a nonsense quantity — a “kilometre” of decay means something different at every latitude. Reproject every layer to a single projected CRS before you compute a centroid or a distance, following Coordinate Reference Systems for Public Health; an equal-area projection such as EPSG:5070 for CONUS-scale work, or the local UTM/State Plane zone for a single region, keeps the decay geometry honest.

Minimum inputs before a fit:

  • Response: a disease rate, ideally variance-stabilized (an empirical-Bayes smoothed rate or a log/arcsine-transformed rate) so that a few tiny-denominator units do not dominate the local fits.
  • Covariates: numeric, standardized to zero mean and unit variance so local coefficients are comparable in magnitude across space and across variables.
  • Geometry: polygon units reduced to representative interior points, or event points aggregated to units; centroids must fall inside their polygon (use a point-on-surface fallback for concave shapes).
  • CRS: one documented projected CRS for all layers; no geographic-coordinate arithmetic anywhere.
  • Unit count: enough units that the smallest local sample (set by the bandwidth) still exceeds the number of covariates by a comfortable margin.

The neighborhood structure GWR imposes through its kernel is conceptually the same object formalized in Spatial Weights Matrix Construction — a continuous distance-decay weighting rather than a binary contiguity rule, but a spatial-weights specification nonetheless, and one that must be logged with the same discipline.

Production Implementation

The reference implementation uses the mgwr package with libpysal and GeoPandas. It searches for a bandwidth by minimizing the corrected Akaike Information Criterion (AICc), fits the model, and extracts local coefficients together with local tt-values so that non-significant local estimates can be masked out on the map.

# Pinned: mgwr==2.2.1, libpysal==4.12.1, geopandas==0.14.4, numpy==1.26.4, scipy==1.13.1
import hashlib
import json
import logging
import numpy as np
import geopandas as gpd
from mgwr.sel_bw import Sel_BW
from mgwr.gwr import GWR

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


def fit_gwr(gdf, y_col, x_cols, *, kernel="bisquare", fixed=False, criterion="AICc"):
    """Fit GWR on a projected GeoDataFrame and return coefficients + local t-values.

    gdf     : GeoDataFrame with a projected CRS and a stable 'geoid' column
    y_col   : disease-rate response column (variance-stabilized upstream)
    x_cols  : list of covariate columns
    kernel  : 'bisquare' (compact) or 'gaussian' (continuous)
    fixed   : False -> adaptive nearest-neighbour bandwidth; True -> fixed distance
    """
    # 1. CRS gate: distances are meaningless in degrees.
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("Reproject to a projected/metric CRS before GWR.")

    # 2. Deterministic ordering => reproducible bandwidth search and output.
    gdf = gdf.sort_values("geoid").reset_index(drop=True)

    # 3. Coordinates from representative interior points (safe for concave polygons).
    pts = gdf.geometry.representative_point()
    coords = np.column_stack([pts.x.to_numpy(), pts.y.to_numpy()])

    # 4. Response and standardized design matrix.
    y = gdf[[y_col]].to_numpy(float)
    X = gdf[x_cols].to_numpy(float)
    X = (X - X.mean(axis=0)) / X.std(axis=0)  # z-scores => comparable local betas

    # 5. Data-driven bandwidth search (adaptive kNN count or fixed distance).
    selector = Sel_BW(coords, y, X, kernel=kernel, fixed=fixed)
    bw = selector.search(criterion=criterion)
    log.info("bandwidth search complete: bw=%s kernel=%s fixed=%s crit=%s",
             bw, kernel, fixed, criterion)

    # 6. Fit and summarize.
    model = GWR(coords, y, X, bw=bw, kernel=kernel, fixed=fixed)
    results = model.fit()
    log.info("GWR fit: AICc=%.2f R2=%.3f ENP=%.1f",
             results.aicc, results.R2, results.ENP)

    # 7. Local coefficients (col 0 = intercept) and filtered local t-values.
    #    filter_tvals applies a dependent multiple-comparison correction.
    betas = results.params                      # (n_units, k+1)
    t_filtered = results.filter_tvals()         # non-sig local t-values -> 0

    out = gdf.copy()
    for k, name in enumerate(["intercept"] + list(x_cols)):
        out[f"beta_{name}"] = betas[:, k]
        out[f"t_{name}"] = t_filtered[:, k]

    # 8. Provenance for audit: the run is only defensible if it re-runs identically.
    provenance = {
        "mgwr_version": "2.2.1",
        "kernel": kernel,
        "fixed_bandwidth": fixed,
        "bandwidth": float(bw),
        "criterion": criterion,
        "n_units": int(len(gdf)),
        "covariates": list(x_cols),
        "config_sha256": hashlib.sha256(
            json.dumps({"kernel": kernel, "fixed": fixed, "bw": float(bw),
                        "crit": criterion, "x": list(x_cols)}, sort_keys=True).encode()
        ).hexdigest(),
    }
    log.info("provenance=%s", json.dumps(provenance))
    return out, results, provenance

The returned GeoDataFrame carries one beta_<covariate> column and one t_<covariate> column per predictor. Mapping the beta column produces the coefficient surface; masking cells where the filtered tt-value is zero prevents you from over-interpreting local estimates that are not distinguishable from no effect. Because filter_tvals() applies a correction that accounts for the many local hypotheses being tested at once, the surviving significant cells are far more trustworthy than a raw per-unit cutoff would leave.

# Map local coefficients, masking locally non-significant units.
# Pinned: geopandas==0.14.4, matplotlib==3.8.4
import matplotlib.pyplot as plt

def map_local_coefficient(out, covariate, ax=None):
    col, tcol = f"beta_{covariate}", f"t_{covariate}"
    ax = ax or plt.gca()
    signif = out[out[tcol] != 0]     # keep only locally significant units
    masked = out[out[tcol] == 0]
    masked.plot(ax=ax, color="none", edgecolor="grey", linewidth=0.3)  # hatch-out
    signif.plot(ax=ax, column=col, cmap="RdBu", legend=True,
                edgecolor="white", linewidth=0.2)
    ax.set_title(f"Local coefficient for {covariate} (non-sig units left blank)")
    ax.set_axis_off()
    return ax

Parameter Selection & Tuning

Four decisions govern a GWR fit, and each must be recorded, not defaulted:

  • Kernel shape. Bisquare is compact and usually preferred for surveillance data because each local fit is bounded to a genuine neighborhood; Gaussian is smoother but lets every observation influence every fit, which can leak signal across sharp regional boundaries. Report which you used — the coefficient surface is not comparable across kernels.
  • Fixed vs adaptive bandwidth. A fixed-distance kernel uses the same radius everywhere; an adaptive kernel uses the same number of neighbors everywhere, expanding in sparse areas and contracting in dense ones. For the heterogeneous population densities typical of health geography, adaptive is the safer default — the reasoning and diagnostics are worked through in the dedicated bandwidth guide.
  • Bandwidth criterion. Sel_BW.search() minimizes AICc by default, which penalizes the effective number of parameters and so guards against the overfitting a raw cross-validation score can miss. Cross-validation (criterion="CV") is an alternative but tends to select smaller, noisier bandwidths.
  • Effective degrees of freedom. After fitting, check results.ENP (effective number of parameters) and results.tr_S. A model whose ENP approaches the unit count is overfit — the bandwidth is too small and every local fit is memorizing its neighborhood.

When you suspect covariates act at genuinely different spatial scales, refit with MGWR and compare AICc; a materially lower AICc for MGWR is evidence that forcing all coefficients through one bandwidth was distorting the surfaces. Whatever you choose, hold the random state and library versions fixed so the bandwidth search is reproducible.

Edge Cases & Failure Modes

  • Local multicollinearity. Two covariates can be uncorrelated globally yet nearly collinear inside a small local window, which inflates and destabilizes local coefficients. Compute the local condition number and local variance-inflation factors (mgwr exposes local_collinearity() on the results object); condition numbers above ~30 in a region mean those local betas are not interpretable there.
  • Small local samples. An adaptive bandwidth that resolves to a neighbor count barely above the covariate count yields high-variance, near-singular local fits. Raise the minimum neighbor count or drop a covariate rather than trusting the surface.
  • Overfitting from a tiny bandwidth. If the search returns a bandwidth so small that ENP nears nn, the “surface” is noise. Inspect AICc across the search path; a criterion curve with no clear minimum is a warning that the data do not support local estimation at all.
  • Edge effects. Units at the study-area boundary borrow from fewer neighbors on one side, biasing their coefficients. Pad the study area with a buffer of units beyond the reporting boundary and clip the surface to the reporting area only at output.
  • Unstable rates in the response. A raw rate from a tiny denominator is dominated by Poisson noise; feed GWR a variance-stabilized rate, or the local fits will chase sampling artefacts. Diagnosing whether residual spatial structure remains is exactly the job of Global & Local Moran’s I Implementation — run Moran’s I on the GWR residuals; near-zero autocorrelation confirms the local model has absorbed the spatial structure a global model left in its residuals.

Compliance & Audit Controls

A GWR coefficient surface that informs where a health department targets an intervention must be reconstructable on demand. The non-negotiable record for one run is the pinned mgwr version, the kernel, the fixed-vs-adaptive flag, the searched bandwidth, the selection criterion, the covariate list and their standardization, the CRS authority code, and a SHA-256 hash binding that configuration to the output. Because the bandwidth is searched rather than set by hand, logging the resolved value is what lets a reviewer skip the search and re-fit deterministically.

# Persist the coefficient surface with embedded provenance for audit.
# Pinned: geopandas==0.14.4, pyarrow==16.1.0
import datetime as dt

def write_gwr_surface(out, provenance, path):
    provenance = dict(provenance)
    provenance["crs_authority"] = out.crs.to_authority()[1]
    provenance["generated_utc"] = dt.datetime.now(dt.timezone.utc).isoformat()
    # GeoParquet preserves CRS, column types, and bbox; embed lineage in metadata.
    out.to_parquet(path, index=True, schema_metadata={
        "title": "GWR local coefficient surface",
        "lineage": "GWR over projected units; standardized covariates; "
                   "AICc-selected bandwidth; local t-values multiple-comparison filtered",
        **{k: str(v) for k, v in provenance.items()},
    })
    return provenance["config_sha256"]

Suppress or hatch-out locally non-significant units in any published map, and never release a coefficient surface without the companion tt-value surface — a strongly-coloured but non-significant region is the single most common way GWR output is over-read.

Production Implementation Checklist