How to Run Geographically Weighted Regression in mgwr

TL;DR — Build three arrays (coords shape (n,2), y shape (n,1), X shape (n,k)), find the bandwidth with bw = Sel_BW(coords, y, X).search(), fit with results = GWR(coords, y, X, bw).fit(), then read the local coefficient surface results.params, per-location results.localR2, and results.summary(). mgwr adds the intercept for you; keep everything projected in metres.

Why this matters #

Geographically weighted regression turns a single global slope into a mappable surface of local slopes, exposing where a relationship strengthens, weakens, or flips sign across your study area. Running it correctly in mgwr is mostly about getting the array shapes and the bandwidth right; the fit itself is two lines. This how-to is the hands-on companion to the Geographically Weighted Regression in Python guide, and it plugs into the wider Python workflows for spatial modeling and regression pipeline that prepares the GeoDataFrame you start from.

Environment #

pip install "mgwr>=2.2" "libpysal>=4.8" "spreg>=1.4" \
            "geopandas>=0.13" "numpy>=1.22" "scipy>=1.9"

import numpy as np                 # >= 1.22
import geopandas as gpd            # >= 0.13
from mgwr.sel_bw import Sel_BW     # mgwr >= 2.2
from mgwr.gwr import GWR

mgwr consumes NumPy arrays, not GeoDataFrames, so the whole job is: prepare arrays, select a bandwidth, fit, read. All coordinates must be projected (metres) because the kernel measures distance directly from the coordinate array.

Step-by-step implementation #

Step 1 — Prepare projected coordinates, X, and y #

gdf = gpd.read_file("data/tracts.gpkg").to_crs("EPSG:32633")   # metric CRS
gdf = gdf.reset_index(drop=True)                                # 0-based index

# (n, 2) coordinate array from polygon centroids (or point geometry directly)
coords = np.array(list(zip(gdf.geometry.centroid.x, gdf.geometry.centroid.y)))

# Response as a COLUMN vector (n, 1) — a 1-D array raises a shape error
y = gdf["median_price"].values.reshape(-1, 1)

# Design matrix (n, k). Do NOT add a constant; mgwr adds the intercept.
X = gdf[["income", "rooms", "dist_cbd"]].values

# Standardise predictors: comparable local betas + a better-behaved search
X = (X - X.mean(axis=0)) / X.std(axis=0)

assert coords.shape[0] == y.shape[0] == X.shape[0], "row counts must match"

Row alignment across coords, y, and X is the single most common failure point — always reset_index(drop=True) before extracting arrays.

Step 2 — Select the bandwidth with Sel_BW #

# fixed=False -> adaptive bandwidth (a nearest-neighbour count)
# kernel='bisquare' -> compact support, the standard default
selector = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)

bw = selector.search(criterion="AICc", search_method="golden_section")
print(f"Selected bandwidth: {bw:.0f} nearest neighbours")

search() returns one scalar: for an adaptive kernel it is the optimal number of neighbours, for a fixed kernel it is a distance in metres. The golden-section search minimises AICc by default.

Step 3 — Fit the GWR model #

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

The bw, kernel, and fixed arguments must match what you searched with — otherwise the fitted smoothing scale silently diverges from the selected one.

Step 4 — Read params, localR2, and the summary #

# Local coefficient surface: (n, k+1); column 0 is the intercept
betas = results.params
print("income local slope: "
      f"{betas[:,1].min():.3f} to {betas[:,1].max():.3f} (global-ish mean {betas[:,1].mean():.3f})")

# Per-location goodness of fit: (n, 1)
print(f"local R2 range: {results.localR2.min():.2f} to {results.localR2.max():.2f}")

# Global diagnostics: R2, AICc, effective parameters, sigma
results.summary()
print(f"AICc: {results.aicc:.1f}   effective params tr(S): {results.tr_S:.1f}")

# Which local estimates are significant (adjusted critical t-value)?
sig_mask = results.filter_tvals()      # (n, k+1); 0 where not significant
n_sig_income = int((sig_mask[:, 1] != 0).sum())
print(f"income locally significant at {n_sig_income} of {len(betas)} locations")

Step 5 — Attach the surface back to the GeoDataFrame #

gdf["income_beta"]   = results.params[:, 1]
gdf["local_r2"]      = results.localR2.flatten()
gdf["income_signif"] = results.filter_tvals()[:, 1] != 0

# Now gdf.plot(column="income_beta", ...) maps the local coefficient surface

Interpreting the output #

results.params is the headline result: an $(n, k+1)$ matrix where each row is one location’s local regression and each column a covariate (column 0 is the intercept). Mapping a column shows the spatial pattern of that covariate’s effect.

results.localR2 is an $(n, 1)$ array of per-location fit. High values mark places where your covariates explain the outcome well; low pockets flag a missing local driver.

results.summary() prints the global picture — pseudo $R^2$, AICc, residual sum of squares, and tr(S), the effective number of parameters. Use AICc to compare kernels and bandwidth types; the model with the lowest AICc wins.

results.filter_tvals() returns the coefficient array with non-significant entries zeroed out, using a critical $t$-value corrected for the many dependent local tests. Always significance-filter before mapping; raw local coefficients over-interpret noise.

Critical best practices #

Reshape y to a column vector #

mgwr expects y with shape (n, 1). A flat (n,) array raises a broadcasting or shape error deep inside the fit. Always y = y.reshape(-1, 1) immediately after extraction.

Never add your own intercept #

mgwr prepends the constant column internally. If you add one, you get a duplicated intercept, a rank-deficient design, and a LinAlgError. Pass X with covariates only, and remember params therefore has k+1 columns.

Keep coordinates projected in metres #

The kernel computes Euclidean distances straight from coords. Latitude/longitude degrees make the bandwidth meaningless and distort every local weight. Reproject to a UTM zone or other metric CRS before extracting coordinates, exactly as covered in GeoPandas data preparation.

Match fit settings to search settings #

The kernel and fixed flags you pass to GWR must equal those you passed to Sel_BW. Searching an adaptive bandwidth and then fitting with fixed=True applies the neighbour count as a distance, producing nonsense with no error.

Standardise predictors #

Standardising X makes local coefficients comparable across covariates and keeps the bandwidth search numerically stable, especially when predictors span very different units.

Troubleshooting #

Next steps #

The bandwidth you plugged in here controls everything — learn to choose and stress-test it in choosing GWR bandwidth with golden-section search. Return to the Geographically Weighted Regression in Python guide for MGWR and local multicollinearity diagnostics.

Frequently Asked Questions #

Does mgwr add the intercept column for me?

Yes. Pass X without a constant column. mgwr prepends the intercept automatically, so results.params has k+1 columns where column 0 is the local intercept and columns 1…k follow the order of your X columns.

What shape must y and coords be for GWR?

coords must be an (n, 2) array of projected x, y coordinates, y must be a two-dimensional (n, 1) column vector, and X must be (n, k). A one-dimensional y raises a shape error; reshape with y.reshape(-1, 1).

How do I know which local coefficients are significant?

Call results.filter_tvals(), which returns the params array with non-significant local estimates set to zero after correcting the critical $t$-value for multiple dependent tests. Map only the surviving coefficients.

Related

• Geographically Weighted Regression in Python — kernels, MGWR, and full interpretation
• Choosing GWR Bandwidth with Golden-Section Search — how Sel_BW picks the bandwidth
• GeoPandas Data Preparation — projecting and cleaning geometry before extracting arrays

← Back to Geographically Weighted Regression in Python