Spatial Epi Public Health GIS

Fitting GWR for Disease-Rate Covariates in Python

When a global regression reports one coefficient for poverty on disease incidence, it collapses every regional regime into a single average — and that average can hide that the effect is strong in one part of the study area and reversed in another. This guide, part of Geographically Weighted Regression for Disease-Rate Modeling, walks a focused Python pipeline that maps exactly where a covariate’s effect is strong, weak, or flipped, and writes the result as an auditable coefficient layer.

Problem Context & Constraints

Suppose you regress county-level respiratory-disease incidence on a poverty index. Ordinary least squares returns a coefficient of, say, +0.31+0.31 with a tidy p-value, and a report concludes “poverty raises incidence.” The trouble is that this one number is the population-weighted blend of several distinct local relationships. In a dense industrial belt the true local slope might be +0.7+0.7; in an affluent commuter ring it might be near zero; in a high-altitude rural interior, where a different exposure dominates, it might be −0.2-0.2. The global fit averages these into a coefficient that describes none of the three regions and would change if the study boundary were redrawn.

Geographically weighted regression (GWR) fixes this by re-fitting a distance-weighted regression at every unit, so each unit gets its own coefficient. The output is not a number but a surface — one local slope per unit that you can map, letting an analyst see the industrial belt light up positive and the rural interior turn negative. Three constraints shape a correct implementation:

  • Distances must be metric. GWR weights neighbors by geographic distance, so the whole model is invalid unless coordinates are in a projected CRS. Degrees will not do.
  • Coefficients must be comparable. Local slopes are only interpretable side by side if the covariates are on a common scale, which means standardizing them before the fit.
  • Local significance is a multiple-comparison problem. Testing one hypothesis per unit across hundreds of units manufactures false “significant” cells unless the local t-values are corrected.

The pipeline below is deliberately linear so each of those constraints maps to one step:

Focused GWR Fitting Pipeline A six-step pipeline in two rows. Row one: build coordinates from projected centroids, then standardize covariates to z-scores, then search a bandwidth with Sel_BW by minimizing AICc. A wrapping arrow drops to row two: fit GWR at the selected bandwidth, then extract local t-values with a multiple-comparison correction, then write a GeoParquet of local coefficients with embedded provenance. 1 · Projected centroids metric coords 2 · Standardize covariates z-scores 3 · Sel_BW search minimize AICc selected bandwidth 4 · GWR fit local betas 5 · Local t-vals + correction filter_tvals 6 · GeoParquet + provenance audit-ready Each step enforces one constraint — metric distance, comparable scale, corrected local significance

Prerequisites

Pin the stack so the bandwidth search and the output hash reproduce exactly across reviewers:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • mgwr 2.2.1 — the GWR/MGWR estimator and bandwidth search
  • libpysal 4.12.1 — spatial data structures mgwr builds on
  • numpy 1.26.4, pyarrow 16.1.0 for GeoParquet output

Input state assumed by the code:

  • Units as a polygon GeoDataFrame with a stable geoid column, one variance-stabilized disease-rate column, and one or more numeric covariate columns.
  • CRS: a single projected/metric CRS for the layer. If you are unsure the reprojection is correct, resolve it with Coordinate Reference Systems for Public Health before running anything here — GWR on geographic degrees is silently wrong, not merely imprecise.
  • Response already smoothed upstream (empirical-Bayes or a log/arcsine transform) so a handful of tiny-denominator units do not dominate the local fits.

Step-by-Step Solution

The function below builds coordinates from projected centroids, standardizes the covariates, searches an adaptive bandwidth by AICc, fits GWR, and extracts local coefficients alongside multiple-comparison-filtered local t-values. It logs a structured audit line and returns a tidy GeoDataFrame.

# python 3.11
# mgwr==2.2.1  libpysal==4.12.1  geopandas==0.14.4  numpy==1.26.4
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.fit")


def fit_disease_gwr(gdf, y_col, x_cols, *, kernel="bisquare", fixed=False):
    # --- guard: metric CRS or the kernel bandwidth is meaningless -------------
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("Reproject to a projected CRS before GWR (degrees != metres).")

    # --- deterministic order => reproducible search and byte-stable output -----
    gdf = gdf.sort_values("geoid").reset_index(drop=True)

    # --- step 1: coordinates from projected representative points --------------
    pts = gdf.geometry.representative_point()      # guaranteed inside the polygon
    coords = np.column_stack([pts.x.to_numpy(), pts.y.to_numpy()])

    # --- step 2: response + standardized design matrix -------------------------
    y = gdf[[y_col]].to_numpy(float)
    X_raw = gdf[x_cols].to_numpy(float)
    mu, sd = X_raw.mean(axis=0), X_raw.std(axis=0)
    if np.any(sd == 0):
        raise ValueError(f"Zero-variance covariate(s): {np.array(x_cols)[sd == 0]}")
    X = (X_raw - mu) / sd                          # z-scores => comparable betas

    # --- step 3: bandwidth search (adaptive kNN by default) --------------------
    bw = Sel_BW(coords, y, X, kernel=kernel, fixed=fixed).search(criterion="AICc")
    log.info("selected bandwidth: %s (kernel=%s fixed=%s)", bw, kernel, fixed)

    # --- step 4: fit -----------------------------------------------------------
    res = GWR(coords, y, X, bw=bw, kernel=kernel, fixed=fixed).fit()
    log.info("GWR fit: AICc=%.2f R2=%.3f ENP=%.1f", res.aicc, res.R2, res.ENP)

    # --- step 5: local coefficients + corrected local t-values -----------------
    # filter_tvals() applies a dependent multiple-comparison correction and
    # zeroes local t-values that fail it, so only trustworthy locals survive.
    betas = res.params
    t_adj = res.filter_tvals()

    out = gdf[["geoid", "geometry"]].copy()
    names = ["intercept"] + list(x_cols)
    for k, name in enumerate(names):
        out[f"beta_{name}"] = betas[:, k]
        out[f"t_{name}"] = t_adj[:, k]
        out[f"sig_{name}"] = t_adj[:, k] != 0       # locally significant flag

    log.info("audit=%s", json.dumps({
        "mgwr_version": "2.2.1", "kernel": kernel, "fixed": fixed,
        "bandwidth": float(bw), "criterion": "AICc",
        "aicc": round(float(res.aicc), 2), "enp": round(float(res.ENP), 1),
        "n_units": int(len(gdf)), "covariates": list(x_cols),
    }))
    return out, res, float(bw)

Two design choices are load-bearing. First, representative_point() rather than centroid guarantees the coordinate falls inside its polygon even for concave counties, so no unit is anchored in a lake or across a bay. Second, filter_tvals() is what turns a wall of raw per-unit t-tests into defensible local significance — without it, a few dozen units will read “significant” purely by chance and light up the map with phantom local effects.

To compare against the global baseline you were trying to improve on, fit an OLS with the same standardized inputs and read the AICc difference:

# mgwr bundles a plain OLS with matching diagnostics for a like-for-like AICc.
from mgwr.gwr import GWR
import numpy as np

def gwr_vs_ols_aicc(coords, y, X, bw, kernel="bisquare", fixed=False):
    gwr = GWR(coords, y, X, bw=bw, kernel=kernel, fixed=fixed).fit()
    # a "global" GWR with an enormous bandwidth reproduces the OLS fit
    n = len(y)
    ols = GWR(coords, y, X, bw=n - 1, kernel=kernel, fixed=fixed).fit()
    return {"gwr_aicc": float(gwr.aicc), "ols_aicc": float(ols.aicc),
            "improvement": float(ols.aicc - gwr.aicc)}

A meaningfully positive improvement (a rule of thumb is an AICc drop greater than 2–3) is the evidence that letting the coefficients vary was worth the extra parameters — that non-stationarity is real and not an artefact.

Finally, persist the surface with provenance so the layer is auditable:

# geopandas==0.14.4  pyarrow==16.1.0
import datetime as dt, hashlib, json

def write_local_coeffs(out, *, bandwidth, kernel, fixed, covariates, path):
    cfg = {"kernel": kernel, "fixed": fixed, "bandwidth": bandwidth,
           "criterion": "AICc", "covariates": list(covariates)}
    cfg_hash = hashlib.sha256(json.dumps(cfg, sort_keys=True).encode()).hexdigest()
    out.to_parquet(path, index=True, schema_metadata={
        "title": "GWR local coefficient surface",
        "mgwr_version": "2.2.1",
        "crs_authority": out.crs.to_authority()[1],
        "config_sha256": cfg_hash,
        "generated_utc": dt.datetime.now(dt.timezone.utc).isoformat(),
        "lineage": "GWR on projected representative points; standardized covariates; "
                   "AICc bandwidth; local t-values multiple-comparison filtered",
    })
    return cfg_hash

Validation & Edge Cases

1. Local multicollinearity inside the kernel window. Two covariates can be uncorrelated across the whole study area yet nearly collinear inside a small local neighborhood, which makes those local coefficients swing wildly. Check it before you trust any regional pattern:

lc = res.local_collinearity()      # returns local condition numbers (and VIFs)
cond = lc[0].ravel()
flagged = int((cond > 30).sum())
log.warning(json.dumps({"event": "local_collinearity",
                        "units_cond_gt_30": flagged, "max_cond": float(cond.max())}))
WARNING {"event":"local_collinearity","units_cond_gt_30":14,"max_cond":58.7}

Condition numbers above ~30 mark regions where the local betas are not separable — report those areas as indeterminate rather than interpreting the coefficient.

2. A bandwidth so small the surface is noise. If Sel_BW returns an adaptive neighbor count barely above the covariate count, the effective number of parameters approaches the unit count and GWR is memorizing, not modeling:

if res.ENP > 0.8 * len(out):
    log.warning(json.dumps({"event": "overfit_risk", "enp": round(float(res.ENP), 1),
                            "n_units": int(len(out))}))

The fix is not to trust the pretty surface — widen the bandwidth (or drop a covariate) and re-search. Tuning the bandwidth criterion is covered in Choosing Adaptive vs Fixed Bandwidth for GWR.

3. Residual spatial structure the local model failed to absorb. GWR should leave residuals that look spatially random. Map the significant-coefficient surface (sig_<covariate>) and inspect the residuals; if clustering remains, a global spatial-dependence model such as Spatial Lag & Error Regression may be the better tool for this dataset, or the response needs further stabilization.

Compliance Notes

For regulatory defensibility, three artifacts must be logged for every run. Persist the resolved bandwidth, kernel, and fixed flag together with the pinned mgwr version — because the bandwidth is searched rather than set, this triplet is what lets a reviewer re-fit deterministically instead of re-running the search. Record the SHA-256 configuration hash and the CRS authority code embedded by write_local_coeffs, so the exact input parameters are bound to the output layer. And never publish a coefficient map without its companion significance flag: a strongly-coloured but non-significant region is the most common way a GWR surface is over-read, and masking those units is a disclosure-control step, not a cosmetic one.