Geostatistical Modeling Spatial Statistics with Python

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

bash
pip install "mgwr>=2.2" "libpysal>=4.8" "spreg>=1.4" \
            "geopandas>=0.13" "numpy>=1.22" "scipy>=1.9"
python
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

python
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.

The three input arrays and the params surface they produce Four small tables are drawn side by side. coords holds two columns of projected metres, y holds one column of median price, X holds three standardised covariate columns, and results.params holds four columns because mgwr prepends a local intercept. Row 0 of every table is the same tract. Two boxes underneath contrast a flat one-dimensional y, which raises a shape error, with y.reshape(-1, 1), the two-dimensional column vector GWR expects. Three arrays in, one coefficient surface out — and the row order must match coords y X results.params column map row 0 row 1 row 2 row 3 x_m y_m 412380 5812100 413905 5810740 411220 5813560 414770 5811090 (n, 2) metres price 318000 402500 275900 351200 (n, 1) income rooms dist_cbd 0.42 −0.31 1.05 1.18 0.62 −0.44 −0.87 −1.02 0.36 0.05 0.77 −0.90 (n, 3) standardised Sel_BW(…).search() → bw = 48 NN GWR(…, bw).fit() b0 b1 b2 b3 5.62 0.41 0.12 −0.28 5.58 0.47 0.09 −0.31 5.71 0.33 0.15 −0.22 5.49 0.52 0.11 −0.35 (n, k+1) = (n, 4) b0 local intercept b1 income b2 rooms b3 dist_cbd Row i is the same tract in every array — reset_index(drop=True) first, then assert the row counts agree. ✗ y = gdf["median_price"].values shape (n,) — a flat 1-D array fit() raises a shape error deep inside mgwr ✓ y = gdf["median_price"].values.reshape(-1, 1) shape (n, 1) — a 2-D column vector what GWR expects; row order is untouched

Step 2 — Select the bandwidth with Sel_BW

python
# 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.

What bw = 48 nearest neighbours actually weights On the left a curve plots bisquare kernel weight against neighbour rank: it starts at 1.00 at the calibration point, falls through 0.56 at rank 12 and 0.25 at rank 24, reaches exactly zero at rank 48 and stays flat at zero for every further neighbour. On the right the same weights are drawn from above as four nested shaded rings around the calibration point, with a table giving the radius as a fraction of the bandwidth, the neighbour rank sitting at that radius, and the weight there. Three notes below cover compact support, why fixed=False lets the disc breathe with tract density, and why the fit must reuse the search settings. An adaptive bandwidth is a neighbour count — and bisquare weights stop dead at the 48th weight given to the jth nearest neighbour 1.00 0.75 0.50 0.25 0 kernel weight w 0 12 24 36 48 60 neighbour rank j (nearest first) bw = 48 w = 0 w ≈ (1 − j/48)² (locally even density: d(j) ∝ √j) the same weights, seen from above bw radius = bw = the distance out to the 48th nearest neighbour of point i r / bw rank w 0.25 3 0.88 0.50 12 0.56 0.75 27 0.19 1.00 48 0.00 compact support Neighbour 49 and beyond get weight exactly 0 — bisquare is truncated. fixed=False is adaptive bw is a count, so the disc shrinks where tracts crowd and grows where they thin. search and fit must agree Refitting with fixed=True reads 48 as 48 metres — no error, nonsense output.

Step 3 — Fit the GWR model

python
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

python
# 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

python
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)(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)(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 R2R^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 tt-value corrected for the many dependent local tests. Always significance-filter before mapping; raw local coefficients over-interpret noise.

What filter_tvals() removes before you map the surface Left, a five by four lattice of twenty locations draws the raw local income slope from results.params column 1 as graduated circles, ranging from 0.61 in the top left corner down to 0.05 in the bottom right. Centre, each location's absolute t-value is plotted on an axis from 0 to 5; a dashed line marks the naive threshold of 1.96 and a solid line the adjusted critical value of 2.61, and eleven of the twenty dots lie to the right of it. Right, the same lattice after filter_tvals: the eleven significant slopes keep their value and the nine that fall short are drawn as small dashed rings holding a zero. filter_tvals() keeps the local slopes the data can support and zeroes the rest results.params[:, 1] — raw local slope |t| at each of the 20 locations filter_tvals()[:, 1] — safe to map 0.61 0.58 0.52 0.47 0.39 0.55 0.49 0.44 0.36 0.28 0.46 0.38 0.31 0.22 0.14 0.33 0.25 0.17 0.09 0.05 every tract keeps an estimate, however thin its evidence naive 1.96 adjusted critical t = 2.61 0 1 2 3 4 5 |t| of the local income slope 11 clear the adjusted threshold 9 fall short and are set to 0 0.61 0.58 0.52 0.47 0.39 0.55 0.49 0.44 0.36 0 0.46 0.38 0 0 0 0 0 0 0 0 gdf["income_signif"] is False wherever a ring is drawn A zeroed tract is not one where income stops mattering — it is one where this bandwidth leaves too few weighted neighbours to separate the local slope from zero. Widen the bandwidth and some of them return.

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

Symptom Likely cause Fix
ValueError about array shape on fit y is 1-D y = y.reshape(-1, 1)
LinAlgError: singular matrix Duplicate intercept or collinear X in a window Remove your constant column; standardise X; raise bandwidth
Bandwidth search runs very slowly Large n with full distance recomputation Use adaptive bisquare; cap bw_max; subsample to prototype
All coefficients look non-significant Oversmoothing or weak signal Check global model first; try a smaller bandwidth
Kernel weights collapse to near-uniform Coordinates in degrees Reproject to a metric CRS before building coords
params has one more column than expected Forgot mgwr adds the intercept Column 0 is the intercept; covariates start at index 1

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 tt-value for multiple dependent tests. Map only the surviving coefficients.

Related

← Back to Geographically Weighted Regression in Python

Related

Choosing GWR Bandwidth with Golden-Section Search Multiscale GWR vs Classic GWR