Geographically Weighted Regression in Python

Geographically weighted regression (GWR) answers a question that global models cannot: does the relationship between a predictor and an outcome stay the same everywhere, or does it drift across space? A single global slope for elevation on temperature, or for income on house price, averages away real geographic structure. GWR fits a local regression at every location, weighting nearby observations more heavily through a distance-decay kernel, and returns a coefficient surface you can map. It is the standard tool for diagnosing and describing spatial non-stationarity in environmental analysis, urban economics, epidemiology, and ecology, and it sits inside the broader Python Workflows for Spatial Modeling & Regression pipeline, downstream of the same GeoDataFrame preparation and weight construction that feed the global spatial regression models.

Prerequisites #

• Python 3.9+
• mgwr>=2.2, libpysal>=4.8, spreg>=1.4, geopandas>=0.13, numpy>=1.22, scipy>=1.9
• A GeoDataFrame reprojected to a metric CRS (e.g., EPSG:32633). Kernel distances are computed directly from the coordinate array, so geographic degrees produce meaningless bandwidths.
• Coordinate array coords of shape (n, 2), response y of shape (n, 1), and design matrix X of shape (n, k) — all row-aligned after reset_index(drop=True). mgwr adds the intercept column itself.
• Predictors that are not globally collinear; local multicollinearity is amplified inside small kernel windows and destabilises local coefficients.

Mathematical Core #

Local Weighted Least Squares #

GWR replaces the single global estimate $\hat{\boldsymbol{\beta}}$ with a location-specific estimate. At each regression point $i$ with coordinates $(u_i, v_i)$, the coefficients are the weighted least squares solution:

$$\hat{\boldsymbol{\beta}}(u_i, v_i) = \left(\mathbf{X}^\top \mathbf{W}(u_i, v_i)\, \mathbf{X}\right)^{-1} \mathbf{X}^\top \mathbf{W}(u_i, v_i)\, \mathbf{y}$$

$\mathbf{W}(u_i, v_i)$ is an $n \times n$ diagonal matrix whose $j$-th diagonal entry $w_{ij}$ is the geographic weight of observation $j$ when estimating the model at location $i$. Observations close to $i$ receive weights near 1; distant ones receive weights near 0. The estimator is solved $n$ times — once per location — so GWR produces an $n \times (k+1)$ matrix of local coefficients rather than a single vector.

The fitted value at $i$ uses only that location’s coefficients: $\hat{y}_i = \mathbf{x}_i \hat{\boldsymbol{\beta}}(u_i, v_i)$. Because each fit reuses overlapping neighbours, GWR is a form of local smoothing, and its effective number of parameters lies between $k+1$ (fully global) and $n$ (fully local).

Kernel Functions #

The weight $w_{ij}$ is a decreasing function of the distance $d_{ij}$ between locations $i$ and $j$, controlled by a bandwidth $b$. Three kernels are used in practice.

The bisquare kernel is compact — it assigns exactly zero weight beyond the bandwidth, so each local fit uses a finite neighbourhood:

$$w_{ij} = \begin{cases} \left(1 - \left(\dfrac{d_{ij}}{b}\right)^2\right)^2 & d_{ij} \le b \\ 0 & d_{ij} > b \end{cases}$$

The gaussian kernel never reaches zero; every observation contributes, but with exponentially diminishing weight:

$$w_{ij} = \exp\!\left(-\tfrac{1}{2}\left(\dfrac{d_{ij}}{b}\right)^2\right)$$

The exponential kernel also has infinite support but decays faster in the near field:

$$w_{ij} = \exp\!\left(-\dfrac{d_{ij}}{b}\right)$$

Bisquare is the usual default because its compact support keeps each local regression genuinely local and computationally sparse. The gaussian kernel is smoother and preferred when you want every observation to contribute a nonzero weight.

Fixed vs Adaptive Bandwidth #

A fixed bandwidth holds $b$ constant in map units — the kernel is the same physical radius everywhere. This is appropriate only when sampling density is uniform, because in sparse regions a fixed radius may enclose too few points to estimate a stable local model.

An adaptive bandwidth fixes the number of nearest neighbours $b = k$ and lets the physical radius expand in sparse areas and contract in dense ones. For $d_{ij}$ the kernel then uses the distance to the $k$-th nearest neighbour as the local scaling. Adaptive bandwidths are the safer default for irregularly sampled data, which is nearly all real-world spatial data.

Bandwidth Selection #

The bandwidth is the single most consequential parameter: too small overfits noise into jagged local surfaces, too large collapses GWR back toward global OLS. It is chosen by optimising a data-driven criterion. The corrected Akaike Information Criterion penalises the effective degrees of freedom $\text{tr}(\mathbf{S})$ (the trace of the hat matrix that maps $\mathbf{y}$ to $\hat{\mathbf{y}}$):

$$\text{AICc} = 2n\ln(\hat{\sigma}) + n\ln(2\pi) + n\,\frac{n + \text{tr}(\mathbf{S})}{n - 2 - \text{tr}(\mathbf{S})}$$

The alternative is leave-one-out cross-validation, which minimises the sum of squared prediction errors where each point is excluded from its own local fit. AICc is generally preferred because it accounts for model complexity directly and is less prone to selecting overly small bandwidths. mgwr searches the bandwidth space with a golden-section search, covered in depth in choosing GWR bandwidth with golden-section search.

MGWR: Per-Covariate Bandwidths #

Standard GWR forces every covariate to vary at the same spatial scale — one bandwidth governs all of them. That is rarely true: a demographic driver may operate near-globally while a local amenity effect is highly localised. Multiscale GWR (MGWR) estimates a separate bandwidth $b_m$ for each covariate $m$ by fitting the model as a generalised additive back-fitting procedure:

$$y_i = \sum_{m=0}^{k} \hat{\beta}_m(u_i, v_i)\, x_{im} + \varepsilon_i, \qquad \text{each } \hat{\beta}_m \text{ smoothed at its own bandwidth } b_m$$

Each covariate’s optimal bandwidth is itself an interpretable result: a near-global bandwidth signals a stationary process, a small bandwidth a strongly local one.

Annotated Implementation #

Step 1 — Build Coordinate and Design Arrays #

import geopandas as gpd
import numpy as np

# Load and project to a metric CRS so kernel distances are in metres
gdf = gpd.read_file("data/tracts.gpkg").to_crs("EPSG:32633")
gdf = gdf.reset_index(drop=True)

# mgwr consumes plain NumPy arrays, not a GeoDataFrame
coords = np.array(list(zip(gdf.geometry.centroid.x, gdf.geometry.centroid.y)))
y = gdf["median_price"].values.reshape(-1, 1)          # shape (n, 1)
X = gdf[["income", "rooms", "dist_cbd"]].values         # shape (n, 3)

# Standardising X aids numerical stability and comparability of local betas
X = (X - X.mean(axis=0)) / X.std(axis=0)

Standardising the predictors is strongly recommended: local coefficients become comparable across covariates and the golden-section search behaves better numerically.

Step 2 — Select the Bandwidth #

from mgwr.sel_bw import Sel_BW

# kernel: 'bisquare' (default, compact) | 'gaussian' | 'exponential'
# fixed=False -> adaptive (nearest-neighbour count); fixed=True -> map distance
selector = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)

# Golden-section search minimising the corrected AIC
bw = selector.search(criterion="AICc", search_method="golden_section")
print(f"Optimal adaptive bandwidth: {bw:.0f} nearest neighbours")

Sel_BW returns a single scalar for GWR: either a neighbour count (adaptive) or a distance (fixed). The criterion may be "AICc", "AIC", "BIC", or "CV".

Step 3 — Fit the GWR Model #

from mgwr.gwr import GWR

gwr_model = GWR(coords, y, X, bw, kernel="bisquare", fixed=False)
gwr_results = gwr_model.fit()

# Local coefficient surface: (n, k+1); column 0 is the intercept
local_betas = gwr_results.params
print("Local slope range for 'income':",
      local_betas[:, 1].min().round(3), "to", local_betas[:, 1].max().round(3))

gwr_results.summary()   # global diagnostics: R2, AICc, effective params

gwr_results.params is the headline output — the mappable local coefficient surface. gwr_results.summary() prints global fit statistics including the pseudo $R^2$, AICc, and the effective number of parameters $\text{tr}(\mathbf{S})$.

Step 4 — Fit MGWR (Per-Covariate Bandwidths) #

from mgwr.sel_bw import Sel_BW
from mgwr.gwr import MGWR

# multi=True triggers the multiscale back-fitting search
mgwr_selector = Sel_BW(coords, y, X, multi=True, kernel="bisquare", fixed=False)
mgwr_bws = mgwr_selector.search(criterion="AICc")
print("Per-covariate bandwidths (intercept, income, rooms, dist_cbd):", mgwr_bws)

mgwr_model = MGWR(coords, y, X, mgwr_selector, kernel="bisquare", fixed=False)
mgwr_results = mgwr_model.fit()
mgwr_results.summary()

The returned mgwr_bws array has one entry per column (including the intercept). A large bandwidth for income and a small one for dist_cbd, for example, tells you the income relationship is nearly stationary while distance-to-centre is a strongly local effect.

Diagnostic Configuration #

Two configuration choices dominate GWR quality: the kernel and the bandwidth type. Set them explicitly rather than relying on defaults.

# Compare adaptive vs fixed under the same kernel
sel_adaptive = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)
sel_fixed    = Sel_BW(coords, y, X, kernel="bisquare", fixed=True)

bw_adaptive = sel_adaptive.search(criterion="AICc")   # neighbour count
bw_fixed    = sel_fixed.search(criterion="AICc")      # distance in metres

# Fit both and compare AICc from the summary
res_a = GWR(coords, y, X, bw_adaptive, fixed=False).fit()
res_f = GWR(coords, y, X, bw_fixed,    fixed=True ).fit()
print(f"Adaptive AICc: {res_a.aicc:.1f}   Fixed AICc: {res_f.aicc:.1f}")

Lower AICc wins. For irregularly sampled data the adaptive kernel almost always prevails because it avoids empty windows in sparse regions.

Output Interpretation #

Local Coefficient Surfaces #

results.params is an $(n, k+1)$ array; map each column to see where a covariate’s effect strengthens, weakens, or reverses sign. Never map raw coefficients without significance filtering — many local estimates are not distinguishable from zero. Use the corrected critical $t$-value that adjusts for multiple dependent local tests:

import matplotlib.pyplot as plt

# Boolean mask of locally significant estimates (adjusted for dependency)
tvals_filtered = gwr_results.filter_tvals()   # (n, k+1); 0 where non-significant
gdf["income_beta"] = gwr_results.params[:, 1]
gdf["income_sig"]  = tvals_filtered[:, 1] != 0

fig, ax = plt.subplots(figsize=(9, 8))
gdf.plot(column="income_beta", cmap="RdBu_r", legend=True, ax=ax)
# Hatch or grey-out non-significant tracts
gdf[~gdf["income_sig"]].plot(color="none", edgecolor="grey", hatch="///", ax=ax)
ax.set_title("Local coefficient for income (non-significant hatched)")
plt.tight_layout()

Local R-squared #

results.localR2 is an $(n, 1)$ array giving the goodness of fit of each local model. Mapping it reveals where the specified covariates explain the outcome well and where they fail — low local $R^2$ pockets suggest a missing local driver.

Condition Number and Local Multicollinearity #

Local windows can be near-collinear even when predictors are globally independent, which inflates local coefficient variance and produces the classic “checkerboard” of wildly signed betas. Check it explicitly:

# Local collinearity diagnostics for GWR
local_CN, VDP = gwr_results.local_collinearity()
n_flagged = int((local_CN > 30).sum())
print(f"{n_flagged} of {len(local_CN)} locations have condition number > 30")

A local condition number above roughly 30 flags a window where coefficients cannot be reliably separated. Treat local estimates at those locations with caution: increase the bandwidth, drop a correlated covariate, or fall back to a global model where non-stationarity is weak.

Production Considerations #

Computational cost. GWR solves $n$ local regressions, each $O(k^3)$ in the covariate dimension plus $O(n)$ for distance weighting, so a naive fit is roughly $O(n^2)$ per bandwidth evaluation, and bandwidth selection multiplies that by the number of golden-section iterations. Above about 10,000 observations, cache the distance computations, use the compact bisquare kernel so most weights are zero, and consider tiling the study area. MGWR is heavier still because back-fitting iterates over covariates until convergence.

Memory. The hat matrix $\mathbf{S}$ is $n \times n$; do not materialise it for large $n$. mgwr computes the trace incrementally, but requesting full local covariance surfaces reintroduces the $O(n^2)$ cost.

Reproducibility. Store the kernel, fixed flag, and selected bandwidth alongside the model — GWR is not portable without them, and re-selecting the bandwidth on new data will silently change the smoothing scale.

When GWR is the wrong tool. If the selected bandwidth approaches the maximum (near-global), the data are effectively stationary and a global spatial lag or error specification from the spatial regression models family is simpler and more efficient. GWR describes non-stationarity; it does not manufacture it.

Troubleshooting #

Next Steps #

For an end-to-end runnable walkthrough — arrays, search, fit, and reading results.params and localR2 — see how to run geographically weighted regression in mgwr. To master the parameter that governs everything, the bandwidth, work through choosing GWR bandwidth with golden-section search. When a global model is the right choice instead, compare specifications in the spatial regression models guide, and validate any fitted surface with the geographically honest folds in cross-validation strategies.

Frequently Asked Questions #

What is the difference between GWR and a global spatial regression model?

A global spatial lag or spatial error model estimates one coefficient vector for the whole study area, assuming the relationship between predictors and outcome is stationary in space. Geographically weighted regression relaxes that assumption: it fits a separate weighted least squares regression at every location using a distance-decay kernel, producing a surface of local coefficients that reveals where and how relationships change across the map. GWR is a diagnostic and descriptive tool for spatial non-stationarity, not a replacement for a global model when relationships are genuinely constant.

Should I use an adaptive or a fixed bandwidth in GWR?

Use an adaptive bandwidth, expressed as a nearest-neighbour count, when sample density varies across the study area, which is the common case for point data and administrative units. The kernel then expands in sparse regions and contracts in dense ones, keeping the effective sample per local regression roughly constant. Use a fixed bandwidth, expressed in map distance units, only when observations are evenly and regularly spaced, such as a gridded raster sample. mgwr defaults to adaptive when you pass fixed=False.

What is MGWR and when do I need it?

Multiscale GWR (MGWR) estimates a separate bandwidth for each covariate instead of forcing all relationships to vary at the same spatial scale. Some processes are near-global (a large bandwidth) while others are highly local (a small bandwidth); standard GWR compromises with a single average bandwidth that oversmooths some effects and overfits others. Use MGWR when you expect covariates to operate at different scales, and read each covariate’s bandwidth as an estimate of the spatial extent of its process.

Why are my local GWR coefficients unstable or wildly signed?

Unstable local coefficients usually indicate local multicollinearity: within a small kernel window the predictors are nearly collinear even if they are not globally correlated. Check results.local_collinearity(); local condition numbers above about 30 flag windows where coefficients cannot be separated reliably. Remedies include increasing the bandwidth, dropping or combining correlated covariates, or switching to a global spatial regression model where the non-stationarity is weak.

Related

• How to Run Geographically Weighted Regression in mgwr — the runnable end-to-end fit from arrays to local surfaces
• Choosing GWR Bandwidth with Golden-Section Search — deep dive on Sel_BW, AICc vs CV, and sensitivity
• Spatial Regression Models — the global lag, error, and Durbin models GWR is compared against
• Stationarity & Trend Analysis — testing whether the stationarity assumption GWR relaxes actually holds

← Back to Python Workflows for Spatial Modeling & Regression