Multiscale GWR vs Classic GWR

TL;DR: Classic GWR fits one bandwidth for the whole model. Multiscale GWR fits one per covariate: build Sel_BW(coords, y, X, multi=True), call .search(multi_bw_min=[20]) to get an array of bandwidths, then pass that selector into MGWR(coords, y, X, selector).fit(). Standardise y and X first or back-fitting will not converge. Read the bandwidth array as a result, not as a setting.

Why This Matters #

A geographically weighted regression with one bandwidth makes a strong and rarely examined assumption: that every relationship in the model changes over the same distance. Road density might influence measured particulate concentration over a couple of kilometres, while winter temperature acts over tens of kilometres and an industry index does not vary meaningfully at all across a single county. One bandwidth cannot serve all three. The golden-section search will find the value that minimises AICc for the model as a whole, which means it over-smooths the sharp term and under-smooths the flat one, and it does both silently — the summary output looks perfectly healthy. That single compromise number is the main limitation of the workflow described in How to Run Geographically Weighted Regression in mgwr, and it is why Choosing GWR Bandwidth with Golden Search spends so much effort on a quantity that is, in the end, an average of several different scales.

Multiscale GWR relaxes exactly that. It treats the model as an additive sum of smooth terms, one per covariate plus the intercept, and gives each its own bandwidth found by back-fitting. The bandwidths that come back are the headline output. They are an estimate of the spatial scale at which each relationship operates, which is a substantive finding about the process rather than a tuning artefact — the empirical counterpart to the questions raised in Spatial Scale & the Modifiable Areal Unit Problem. The whole method sits inside the broader family covered under Spatial Regression Models, and the price is real: back-fitting costs roughly $(k+1)$ bandwidth searches per sweep, and it needs several sweeps.

Classic GWR estimates

$$y_i = \beta_0(\mathbf{u}_i) + \sum_{j=1}^{k} \beta_j(\mathbf{u}_i)\, x_{ij} + \varepsilon_i$$

with every $\beta_j(\cdot)$ smoothed by the same kernel of bandwidth $b$. MGWR keeps the same expression but attaches a separate bandwidth $b_j$ to each term, so the smoother for covariate $j$ is $S_j$ built from $b_j$ alone. The fitted surface is then a sum of terms with different roughness, which is why the model is estimated by back-fitting rather than in one weighted least-squares pass.

Environment and Version Pinning #

MGWR lives in the same package as GWR, so nothing extra is required beyond a recent mgwr. The multiscale code paths were substantially reworked at 2.1, and n_jobs on the selector arrived at 2.2, so pin at or above that.

pip install "mgwr>=2.2.1" "libpysal>=4.9.0" "spreg>=1.4.0" \
            "geopandas>=1.0" "numpy>=1.23,<2.1" "scipy>=1.9" "pandas>=2.0"

import numpy as np                     # >= 1.23
import pandas as pd                    # >= 2.0
from mgwr.sel_bw import Sel_BW         # mgwr >= 2.2.1
from mgwr.gwr import GWR, MGWR

Step-by-Step Implementation #

Step 1 — Build a dataset whose scales are known #

Comparing the two models on real data tells you which fits better but not whether the bandwidths are right. Constructing the coefficient surfaces yourself makes the bandwidth table testable: you know what the answer should be.

rng = np.random.default_rng(11)
SIDE_KM, STEP_KM = 40.0, 2.0

gx, gy = np.meshgrid(np.arange(1.0, SIDE_KM, STEP_KM),
                     np.arange(1.0, SIDE_KM, STEP_KM))
gx, gy = gx.ravel(), gy.ravel()          # 20 x 20 = 400 sites
ux, uy = gx / SIDE_KM, gy / SIDE_KM

# Coefficient surfaces at three deliberately different scales.
b_road = 1.2 * np.sin(2 * np.pi * 2 * ux)      # two cycles across x -> local
b_temp = 0.9 * (ux + uy) / 2                   # smooth ramp        -> regional
b_ind  = np.full_like(ux, 0.7)                 # constant           -> global
b_0    = 0.5 * np.sin(np.pi * ux)

x_road = rng.normal(size=400)
x_temp = rng.normal(size=400)
x_ind  = rng.normal(size=400)

y_raw = (b_0 + b_road * x_road + b_temp * x_temp + b_ind * x_ind
         + rng.normal(0, 0.25, size=400))

coords = np.column_stack([gx * 1000.0, gy * 1000.0])   # metres, projected
sd_y = y_raw.std()
print(f"n = {len(y_raw)}, sd(y) = {sd_y:.3f}")
for name, b in [("road_density", b_road), ("winter_temp", b_temp),
                ("industry_index", b_ind)]:
    print(f"true {name:<15} slope on standardised y: "
          f"{(b / sd_y).min():+.2f} to {(b / sd_y).max():+.2f}")

n = 400, sd(y) = 1.238
true road_density    slope on standardised y: -0.97 to +0.97
true winter_temp     slope on standardised y: +0.02 to +0.71
true industry_index  slope on standardised y: +0.57 to +0.57

Step 2 — Standardise both sides #

MGWR back-fits one term at a time and judges convergence with a single score over the whole additive model. Columns on different numeric scales make that score dominated by whichever term happens to be largest, and the sweep stalls. Centre and scale everything.

X = np.column_stack([x_road, x_temp, x_ind])
X = (X - X.mean(axis=0)) / X.std(axis=0)

y = y_raw.reshape(-1, 1)
y = (y - y.mean()) / y.std()          # (400, 1) column vector

names = ["intercept", "road_density", "winter_temp", "industry_index"]
assert coords.shape == (400, 2) and y.shape == (400, 1) and X.shape == (400, 3)

Do not add a constant column: both GWR and MGWR take constant=True by default and prepend one for you, which is why names has four entries for three covariates.

Step 3 — Fit classic GWR as the baseline #

Always run the single-bandwidth model first. It is fast, and if it cannot beat OLS there is no local variation for MGWR to decompose into scales.

gwr_sel = Sel_BW(coords, y, X, kernel="bisquare", fixed=False)
gwr_bw = gwr_sel.search(criterion="AICc", search_method="golden_section")

gwr_res = GWR(coords, y, X, gwr_bw, kernel="bisquare", fixed=False).fit()
print(f"GWR bandwidth : {gwr_bw:.0f} nearest neighbours")
print(f"GWR AICc      : {gwr_res.aicc:.1f}   R2 = {gwr_res.R2:.3f}   "
      f"tr(S) = {gwr_res.tr_S:.1f}")

GWR bandwidth : 97 nearest neighbours
GWR AICc      : 301.2   R2 = 0.925   tr(S) = 58.2

For reference, ordinary least squares on the same arrays returns an AICc of 936.9 and an $R^2$ of 0.412, so local variation plainly exists. The single bandwidth of 97 is the compromise: it lies between the three true scales and equals none of them.

Step 4 — Search per-covariate bandwidths #

multi=True switches Sel_BW from a one-dimensional golden-section search to the back-fitting routine. search() then returns an array of length $k+1$ instead of a scalar.

mgwr_sel = Sel_BW(coords, y, X, kernel="bisquare", fixed=False,
                  multi=True, constant=True)

bws = mgwr_sel.search(criterion="AICc",
                      multi_bw_min=[20],        # floor, broadcast to all columns
                      tol_multi=1.0e-5,         # convergence on the SOC-RSS score
                      max_iter_multi=200,
                      verbose=False)

print("bandwidths:", bws.astype(int))

bandwidths: [118  43 186 399]

The floor of 20 stops a noisy column from collapsing onto a handful of neighbours; the ceiling defaults to $n-1 = 399$, and industry_index has been pushed all the way to it.

Step 5 — Fit MGWR from the fitted selector #

MGWR takes the selector object itself, not the bandwidth array, because it needs the intermediate state the back-fitting search produced.

mgwr_res = MGWR(coords, y, X, mgwr_sel, kernel="bisquare",
                fixed=False, constant=True, sigma2_v1=True).fit()

print(f"MGWR AICc : {mgwr_res.aicc:.1f}   R2 = {mgwr_res.R2:.3f}   "
      f"tr(S) = {mgwr_res.tr_S:.1f}")
print(f"AICc improvement over GWR: {gwr_res.aicc - mgwr_res.aicc:.1f}")

MGWR AICc : 198.0   R2 = 0.933   tr(S) = 42.8
AICc improvement over GWR: 103.2

Step 6 — Assemble the bandwidth table #

This table is the output people actually read. Convert each bandwidth into an approximate kernel radius using the sampling density, and pair it with ENP_j, the effective number of parameters that covariate’s surface consumes.

density = len(coords) / (SIDE_KM ** 2)              # 0.25 sites per km^2
radius_km = np.sqrt(bws / (np.pi * density))

table = pd.DataFrame({
    "covariate": names,
    "bw_nn": bws.astype(int),
    "radius_km": radius_km.round(1),
    "ENP_j": mgwr_res.ENP_j.round(1),
    "beta_min": mgwr_res.params.min(axis=0).round(2),
    "beta_max": mgwr_res.params.max(axis=0).round(2),
})
print(table.to_string(index=False))

     covariate  bw_nn  radius_km  ENP_j  beta_min  beta_max
     intercept    118       12.3    9.2     -0.21      0.44
  road_density     43        7.4   27.4     -0.92      0.95
   winter_temp    186       15.4    4.9      0.04      0.69
industry_index    399       22.5    1.3      0.56      0.58

Step 7 — Compare the coefficient surfaces #

The bandwidths explain the surfaces, so put the two models side by side against the truth you constructed.

comparison = pd.DataFrame({
    "covariate": names[1:],
    "true_min": [(b / sd_y).min() for b in (b_road, b_temp, b_ind)],
    "true_max": [(b / sd_y).max() for b in (b_road, b_temp, b_ind)],
    "mgwr_min": mgwr_res.params[:, 1:].min(axis=0),
    "mgwr_max": mgwr_res.params[:, 1:].max(axis=0),
    "gwr_min":  gwr_res.params[:, 1:].min(axis=0),
    "gwr_max":  gwr_res.params[:, 1:].max(axis=0),
}).round(2)
print(comparison.to_string(index=False))

     covariate  true_min  true_max  mgwr_min  mgwr_max  gwr_min  gwr_max
  road_density     -0.97      0.97     -0.92      0.95    -0.61     0.64
   winter_temp      0.02      0.71      0.04      0.69    -0.05     0.82
industry_index      0.57      0.57      0.56      0.58     0.42     0.72

Interpreting the Output #

Read the bandwidth column first and the coefficients second. A bandwidth near the search floor means the relationship changes over a short distance and the term is consuming many degrees of freedom — road_density at 43 neighbours takes 27.4 of the model’s 42.8 effective parameters, more than the other three combined. A bandwidth at or near the ceiling of $n-1$ means the kernel includes nearly every observation at every point, so the coefficient surface is flat and the term is global: industry_index spans 0.56 to 0.58, which is a constant with rounding noise. That is a finding worth writing down, because it says the data contain no evidence that this relationship varies at all, and it justifies reporting a single slope for that covariate.

The comparison table shows both ways a single bandwidth fails, in the same fit. With $b = 97$, classic GWR flattens road_density to a range of $-0.61$ to $0.64$ against a true amplitude of $0.97$, losing more than a third of the signal — the kernel is too wide to follow a surface that reverses sign twice across the square. At the same time it grants industry_index a spurious range of $0.42$ to $0.72$, because a bandwidth narrow enough to keep some of the local term also lets a genuinely constant coefficient wander with the noise. Anyone mapping that GWR surface would describe a spatial pattern in the industry effect that does not exist.

AICc adjudicates. The drop from 301.2 to 198.0 is decisive here, though this dataset was built with scales that differ by nearly a factor of ten and real data rarely separates so cleanly; improvements of five to forty are more typical. Note that $R^2$ barely moved, from 0.925 to 0.933, while tr(S) fell from 58.2 to 42.8. MGWR did not fit the data much better — it fit it about as well with substantially fewer effective parameters, spending them where the variation actually is. That is the shape of a real multiscale improvement, and if you see a large $R^2$ gain instead, check for over-fitting at a collapsed bandwidth.

Significance filtering works as it does for GWR, but the critical $t$-value is computed per covariate because each term has its own effective parameter count:

sig = mgwr_res.filter_tvals()                 # (400, 4), zero where not significant
for j, name in enumerate(names):
    print(f"{name:<15} adj alpha = {mgwr_res.adj_alpha_j[j, 1]:.4f}   "
          f"significant at {int((sig[:, j] != 0).sum())} of 400 sites")

intercept       adj alpha = 0.0163   significant at 271 of 400 sites
road_density    adj alpha = 0.0055   significant at 388 of 400 sites
winter_temp     adj alpha = 0.0306   significant at 344 of 400 sites
industry_index  adj alpha = 0.1154   significant at 400 of 400 sites

The adjusted alpha rises as the bandwidth rises, because a globally smoothed term involves far fewer independent tests than a local one. Using an unadjusted 0.05 everywhere would over-declare significance for road_density and under-declare it for industry_index.

Critical Best Practices #

Establish that local variation exists before reaching for MGWR #

MGWR decomposes local variation into scales. It cannot manufacture it. Fit classic GWR first and compare its AICc against OLS; if the improvement is negligible, every MGWR bandwidth will come back near the ceiling after a long wait, and the right answer was a global model all along. Running the cheap model first also gives back-fitting a sensible starting point, which is exactly what Sel_BW uses.

Standardise, and standardise both sides #

Back-fitting sweeps one column at a time and stops when a single score over the additive model settles. On raw units that score is dominated by the term with the largest numeric scale, so the other bandwidths drift for many iterations before the tolerance is met — or max_iter_multi is reached and the search silently returns whatever it had. Scale X column-wise and standardise y as well. The side benefit is that the coefficients in the table are directly comparable in magnitude.

Treat the bandwidth array as a result with uncertainty #

A bandwidth of 43 and a bandwidth of 51 are not meaningfully different; the AICc surface near its optimum is nearly flat, as anyone who has watched a golden-section search knows. Report bandwidths in bands — local, regional, global — rather than as point estimates, and be sceptical of narrative built on a difference of a few neighbours. Where the distinction matters, bootstrap the sites and look at the spread of each bandwidth across resamples.

Budget the runtime honestly #

Classic GWR needs one bandwidth search and one fit. MGWR needs a search per column per sweep, so a four-column model converging in six sweeps performs roughly 24 univariate searches, each involving many weighted least-squares passes over all $n$ points. On this 400-site problem the GWR fit takes well under a second and MGWR takes about half a minute. At $n = 5000$ the gap becomes hours against seconds, because each pass is $O(n^2)$. Use n_jobs=-1 on both the selector and the model, keep hat_matrix=False unless you need the full inference machinery, and never put an MGWR fit inside a cross-validation loop without measuring one fit first.

Watch for local collinearity before trusting a small bandwidth #

Two covariates that are only weakly correlated globally can be strongly correlated inside a small kernel. When that happens back-fitting can push one bandwidth down and the other up, producing a pair of surfaces that trade signal back and forth rather than a genuine scale separation. The symptom is a bandwidth at or near multi_bw_min combined with a wildly swinging coefficient surface. Compute local condition numbers, or refit with the suspect column dropped and check that the remaining bandwidths hold still.

Troubleshooting #

Next Steps #

Once the bandwidths are in hand, map the surfaces and filter them — mgwr.utils.compare_surfaces plots a GWR and an MGWR coefficient map on a shared colour scale, which is the fastest way to show a reviewer what the compromise bandwidth was hiding. If MGWR reports every covariate as effectively global, step back to Spatial Regression Models, where a spatial lag or error specification handles residual dependence without a coefficient surface.

Frequently Asked Questions #

When is MGWR worth the extra cost over classic GWR?

When you have a substantive reason to think your covariates operate at different scales, and classic GWR has already shown that some local variation exists. A model of air quality mixing traffic density, land cover and regional meteorology is the obvious case. If every covariate plausibly acts over the same distance, or if GWR’s AICc barely improves on OLS, MGWR will spend many minutes of back-fitting to return bandwidths that all sit near the classic one, and the simpler model is the honest report.

What does a bandwidth equal to the sample size mean?

It means the kernel for that covariate includes essentially every observation at every calibration point, so its coefficient surface is flat and the term is effectively global. In mgwr the adaptive search ceiling is $n-1$, so a bandwidth reported at that value has been pushed as far towards global as the search allows. Read it as a finding: the data contain no evidence that this relationship varies across the study area, which is a claim worth stating.

Do I have to standardise the variables for MGWR?

In practice yes. Back-fitting updates one covariate at a time against the partial residuals of the others, and the convergence criterion is a single score over the whole additive model. Columns on wildly different numeric scales make that score dominated by one term, so the sweep either converges slowly or stops at max_iter_multi without settling. Centring and scaling y and every column of X also makes the returned coefficients directly comparable in magnitude, which is how the bandwidth table is usually read.

Can I compare GWR and MGWR with AICc directly?

Yes, provided both were fitted to identical y and X arrays with the same kernel and the same adaptive or fixed setting. Both models compute AICc from the trace of their hat matrix, so the effective parameter count is on a comparable footing. A difference of less than about three is not evidence for either model; differences in the tens mean the multiscale decomposition is capturing structure the single bandwidth cannot. Never compare AICc across different standardisations of y.

Related

• How to Run Geographically Weighted Regression in mgwr — the single-bandwidth fit this page compares against
• Choosing GWR Bandwidth with Golden Search — the one-dimensional search that back-fitting runs once per covariate per sweep
• Spatial Scale & the Modifiable Areal Unit Problem — why the scale a relationship operates at is a substantive question

← Back to Geographically Weighted Regression