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 for the entire dataset, it estimates a separate coefficient vector at each location by re-fitting a weighted least squares model in which nearby observations count more than distant ones. The estimator at location is:
where is the design matrix of covariates, is the disease-rate response, and is a diagonal matrix of geographic weights that depends on the location being fitted. Every unit gets its own , 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,
and the bisquare kernel, which is compact — it assigns zero weight beyond the bandwidth and so bounds each local fit to a finite neighborhood,
Here is the distance from fit location to observation and 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 and an adaptive nearest-neighbor 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:
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 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 -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 -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) andresults.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 , 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 -value surface — a strongly-coloured but non-significant region is the single most common way GWR output is over-read.
Production Implementation Checklist
Related Topics
- Disease Clustering & Spatial Statistical Modeling — the parent section placing local regression alongside the autocorrelation and scan-statistic families.
- Fitting GWR for Disease-Rate Covariates in Python — a focused end-to-end mgwr walkthrough with provenance output.
- Choosing Adaptive vs Fixed Bandwidth for GWR — matching the kernel support to uneven population density.
- Spatial Lag & Error Regression — the global alternative that models dependence rather than varying the coefficients.
- Coordinate Reference Systems for Public Health — the projection discipline every distance-weighted method depends on.