直接说结论吧:多尺度地理加权回归(MGWR)这个东西,光看名字就能劝退一拨人,但用过之后再回头去看全局回归或者普通GWR,总觉得哪里差了点火候。我自己的体会是,GWR把"关系随空间变化"这件事做了出来,但它在带宽上搞的是"一刀切"——所有解释变量共用同一个最优带宽。而MGWR打破了这条限制,允许每个变量拥有自己的带宽,有的变量在几百公里尺度上平稳变化,有的变量只在几公里范围内剧烈波动,这种"多尺度"视角才是真实空间过程的写照。
这篇指南的目标是带你把一整条MGWR实证流程走通:从环境准备、数据预处理、权重矩阵构建,到模型拟合、带宽搜索、结果解读与可视化出图,每一步我都会给出可直接复现的代码、参数说明和最容易翻车的坑。我用的是Python生态下的mgwr库,配合geopandas、libpysal和matplotlib,完全避开ArcGIS里那个按钮化流程的黑盒问题——至少你得知道每个参数背后在算什么。
适合阅读这篇文章的人有三类:一是在ArcGIS里点过MGWR工具但完全不理解输出结果的人;二是正在做空间计量毕业论文、需要可解释稳健结果的研究生;三是想把手里的截面空间数据挖得更深,但不满足于全局OLS和普通GWR的进阶用户。如果你只是需要一个"快速跑通"的脚本,后面第三、四节直接抄;如果你想弄懂为什么这么做,第一、二、五节值得细看。
1. 为什么是MGWR:单带宽GWR的局限与多尺度逻辑
1.1 普通GWR的"一刀切"病根在哪
空间异质性(spatial heterogeneity)是空间数据的常识性特征:同一解释变量在不同地理位置上对因变量的作用强度可能完全不同。GWR通过在每个回归点进行局部加权回归来刻画这种变化,权重由该点与周围样本之间的空间距离决定,距离越近权重越大。这思路没毛病,但问题出在带宽的选择方式上。
标准的GWR流程会搜寻一个全局最优带宽——比如通过校正AICc(小样本Akaike信息准则)最小化来确定——然后用这个唯一带宽去计算所有解释变量的局部系数。这意味着代码里所有的变量都被强制假定为"在相同的空间范围内产生效应"。可现实里,一个解释变量可能是区域级的慢变过程(比如海拔对气温的影响,可能在几百公里尺度上都很稳定),另一个则是局地性的快变过程(比如商铺密度对房价的影响,可能只在一两公里内起作用)。如果你让这两类变量共用同一个带宽,结果一定是谁都被扭曲:慢变变量被过度局部化,快变变量被过度平滑化。
我从实操角度补一个例子:在研究城市房价时,"到CBD的距离"这个变量通常的效应尺度较大,跨片区时系数变化平缓;而"周边学校质量"的效应尺度很小,换个学区可能系数就从正转负。普通GWR一旦带宽选小了,"到CBD"的系数会变得异常跳脱,出现毫无经济含义的正负震荡;带宽选大了,"学校质量"的异质性又被抹平。MGWR的出现就是为了解决这种尺度混杂。
1.2 MGWR的计算思路:逐个变量去寻找自己的最优带宽
多说一句,这里的"多尺度"和"多分辨率"是两码事。MGWR的多尺度体现在回归系数估计过程中的带宽定制化,而不是对数据空间网格做多重采样。其核心思路可以这样理解:在每次迭代中,只允许一个解释变量拥有独立的带宽,其余变量暂时采用当前全局最优带宽或上一次迭代结果,再通过后向拟合(back-fitting)反复修正各变量带宽和系数,直到收敛。
它本质上是一种广义可加模型(GAM)在空间扩展上的拟合策略,每个解释变量对应一个光滑项,只是这里的"光滑"受到空间权重的约束,且各光滑项的光滑程度(带宽)自由、独立地被数据驱动出来。正因为拟合过程是迭代的,MGWR远比GWR费时——如果你跑几千个样本还嫌慢,那是正常现象,放到下一节环境部分说。
你不需要掌握后向拟合的全部数学推导,但需要记住三个关键输出:
- 每个解释变量最终对应的最优带宽(连续带宽的核函数带宽值,或者自适应带宽的邻居数),这直接表示该变量效应的空间尺度;
- 每个样本点上每个变量的局部系数估计值;
- 每个变量的局部显著性检验(基于伪t值),这告诉你在哪些地区这个变量显著、在哪些地区不显著。
由此带来的解读价值比普通GWR上了一个台阶。带宽小的变量说明它是局地力量,带宽大的变量更接近全局性因素,甚至接近全局OLS的系数——这个逻辑链条是从结果解读的关键前提。
1.3 MGWR的适用边界:不是所有数据都值得用"多尺度"
我从实际项目里总结出的适用范围判断标准,供你参考:
| 场景 | 适合用MGWR? | 说明 |
|---|---|---|
| 样本量超过500、空间单元为区县/街道级别 | 是 | 局部拟合才有足够样本支撑 |
| 解释变量个数在3~8个之间 | 是 | 变量过多导致后向拟合收敛慢且解释困难 |
| 变量之间的空间作用尺度差异明显 | 强烈推荐 | 这正是MGWR存在的价值 |
| 怀疑存在全局性过程(如宏观政策变量) | 是 | 你可以观察该变量带宽是否极大化 |
| 样本量低于200 | 建议谨慎 | 局部回归邻域内样本太少,伪t值可靠性下降 |
| 只需要总体平均效应 | 否 | 直接跑OLS或者空间误差模型更快、更稳 |
我自己见过不少把MGWR套在小样本乡镇数据的案例,结果每个变量的带宽都跑到样本量上限,局部系数和全局回归几乎一致,只有标准误略有变化,等于白算。所以在开跑之前,用上面的表先卡一遍自己的数据。
2. 环境准备:MGWR依赖链与最容易翻车的地理坐标坑
2.1 安装全套依赖库的具体步骤
MGWR在Python生态中最常用的实现是mgwr库(基于spreg和libpysal),你要装的不止这一个包。完整的依赖链是:geopandas负责矢量数据读写与坐标转换,libpysal提供空间权重矩阵构建,mgwr负责模型拟合,matplotlib用于结果可视化,numpy/pandas做常规数据操作。
从清华镜像源一次性装齐最省事,直接在终端里逐行执行:
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple geopandas libpysal mgwr matplotlib numpy pandas如果你之前装过,怕版本冲突,可以先创建一个独立环境(强烈建议):
conda create -n mgwr_env python=3.9 -y conda activate mgwr_env pip install -i https://pypi.tuna.tsinghua.edu.cn/simple geopandas libpysal mgwr matplotlib numpy pandas这里我特意强调Python版本:mgwr库对Python 3.10及以上版本的兼容性此前一直有些磨合问题,3.8到3.9是最稳的区间。到2025年写这篇文章时,新版mgwr已经基本兼容3.10,但为了减少不必要的包依赖编译报错,新创建环境时选3.9依然是最稳妥的策略。
2.2 地理坐标投影:比任何依赖库都重要的一步
这一节我放到环境部分而不是数据部分,是因为投影错误导致的报错几乎会发生在任何一个后续环节,而且错误信息极其迷惑。
mgwr在构建空间权重矩阵时,需要传入坐标数据来计算样本点之间的距离。如果你直接把经纬度(WGS84)坐标扔进去,会出现两个问题:第一,距离单位是十进制度数,会导致空间核函数的带宽参数含义混乱;第二,在较高纬度地区,经度方向的实际距离会被严重压缩,不同方向的空间关系被扭曲,模型结果变形。
所以你需要先把地理坐标转为投影坐标。最常用的是UTM投影,但不同地区要用对应分带。判断方法很简单:如果你的研究区跨越多个UTM分带,建议改用Albers等积圆锥投影或兰伯特等角圆锥投影,以研究区中心为中央经线自定义投影参数。
用geopandas做投影转换时,代码长这样:
import geopandas as gpd # 读入带几何信息的矢量面数据(比如区县边界加属性表) gdf = gpd.read_file("your_data.shp") # 先确保数据本身是WGS84地理坐标 gdf = gdf.to_crs(epsg=4326) # 转为投影坐标系,这里以UTM 50N为例(中央经线117°E) gdf_proj = gdf.to_crs(epsg=32650) # 提取每个样本点的质心坐标,作为回归坐标 gdf_proj["X"] = gdf_proj.geometry.centroid.x gdf_proj["Y"] = gdf_proj.geometry.centroid.y如果样本数据是点数据,则直接用点坐标即可;如果是面数据,质心坐标是默认选择。需要注意,mgwr拟合时使用的是样本点坐标来构建距离矩阵,而不是共享边界——但空间权重矩阵可以选择基于质心距离或基于邻居关系,这两者的区别在下一节说。
2.3 常见的依赖加载报错与解决
安装顺利不代表导入顺利。我遇到过太多人卡在第一步就心态崩了,这里把几个高频报错列出来:
ImportError: cannot import name 'weights' from 'libpysal':这通常是因为libpysal版本过旧,weights模块的导入路径变了。升级到最新版即可:pip install --upgrade libpysal。AttributeError: module 'spreg' has no attribute 'GMW':说明你单独装了一个很老的spreg,而mgwr依赖新版接口。安装顺序会影响依赖解析,建议统一从镜像源安装最新版本,不要手动指定旧版本。MemoryError在构建距离矩阵时出现:当样本量达到几千时,完整距离矩阵会占用大量内存。建议用libpysal的稀疏权重矩阵,并在MGWR拟合时使用sparse=True参数(后面代码里会给)。
这些报错基本都是环境层面的,和数据质量无关,解决掉之后整个流程就能顺畅往下走了。
3. 数据准备与权重矩阵:别让坐标系和数据格式毁掉全部努力
3.1 数据表结构:从原始数据到回归数据框的整理流程
MGWR的输入数据本质上就是一个结构化数据框,每一行是一个空间样本点(或面单元),每一列是一个变量。变量分为两类:一类是因变量(y),另一类是解释变量(X)。但这里有几个必须提前处理干净的脏数据问题:
- 缺失值处理:
mgwr不会像statsmodels那样自动删掉缺失行,缺失值会直接造成矩阵运算失败。所以用dropna()之前,要想清楚缺失机制——如果是随机缺失,直接删掉问题不大;如果是系统缺失(比如只在偏远地区缺失),删除会损伤样本代表性,最好用插补。 - 多重共线性:MGWR每个变量都拥有独立带宽,但这不等于它能化解共线性问题。如果两个解释变量的相关系数超过0.7,它们会在局部邻域内争夺解释力,导致系数符号极端不稳定。建议在拟合前专门跑一个OLS的VIF检验,把VIF大于10的变量剔除。
- 变量标准化不等于必须做,但强烈推荐:MGWR的带宽搜索基于距离矩阵,变量的量纲差异不会直接影响距离计算,因为距离矩阵只依赖坐标。但是,在后向拟合的迭代算法中,变量的量纲会影响系数初始值和收敛速度,甚至影响伪t值的稳定性。我在实践中会对所有解释变量做z-score标准化(均值0,标准差1),因变量保留原始量纲以便解读。
这里给出一个完成上述过程的参考代码:
import pandas as pd df = pd.read_csv("your_data.csv") # 合并坐标(来自geopandas转换后的字段) df = df.merge(gdf_proj[["id", "X", "Y"]], on="id") # 剔除缺失 df = df.dropna(subset=["y", "x1", "x2", "x3", "X", "Y"]) # z-score标准化(排除坐标) feature_cols = ["x1", "x2", "x3"] df[feature_cols] = (df[feature_cols] - df[feature_cols].mean()) / df[feature_cols].std()3.2 空间权重矩阵的两种选择:连续核函数与自适应核函数
空间权重矩阵是距离加权的基础数据结构。mgwr里最常见的是基于核函数(Kernel)的距离权重矩阵,分为固定带宽型(fixed bandwidth)和自适应带宽型(adaptive bandwidth)。两者区别用一个比喻说清楚:
固定带宽是"以每家店为中心画一个固定半径的圆",圆内顾客都参与加权,圆外不参与;自适应带宽是"不管店铺开在哪里,都招募距离最近的K个人作为顾客",这样人口稀疏地区的圆会自动变大,人口密集地区的圆会自动缩小。
实际数据中,样本点经常分布不均——市中心密集,郊区稀疏。固定带宽下,郊区样本的邻域可能只有两三个邻居,局部回归完全无法估计;自适应带宽能保证每个回归点都有足够的邻居数来拟合。默认并且绝大多数情况下推荐使用自适应带宽(kernel='bisquare'配合邻居数)。只有样本点分布非常均匀时,固定带宽才更好。
用代码构建权重矩阵:
import numpy as np import libpysal from mgwr.gwr import MGWR from mgwr.sel import search_mgwr_paras coords = list(zip(df["X"], df["Y"])) y = df["y"].values.reshape(-1, 1) X = df[feature_cols].values # 构建自适应均匀核权重矩阵的搜索对象 # 其中的gwr_type参数后面会用到 gwr_selector = sel.MGWRSelector(coords, y, X, kernel='bisquare', fixed=False, sparse=True)这里有一个经常被忽略的细节:sparse参数。构建权重矩阵时,libpysal默认生成的是稀疏矩阵,但这只影响存储方式和计算效率,不会改变计算结果。样本量超过1000时务必保持sparse=True,否则内存占用会快速膨胀。
3.3 更精细的空间关系处理:基于邻居权重还是距离权重
除了距离核权重外,有时研究者还想把样本单元之间的邻接关系(如区县共享边界)纳入权重。mgwr本身在拟合时使用的是核函数权重,但你可以通过设置weights参数传入一个二值邻接权重矩阵来限制邻域范围。比如只允许"共享边界的区县"互相进入对方的局部邻域,距离核再在这个子集上计算权重。
这个做法适合行政区单元数据,能避免"隔座山也能被认作近邻"的荒谬情况(比如两个区县质心距离很近但被山脉完全隔开)。实现方式不复杂:
from libpysal.weights import Queen w = Queen.from_dataframe(gdf_proj) gwr_selector = sel.MGWRSelector(coords, y, X, kernel='bisquare', fixed=False, sparse=True, w=w)但我要提醒一句:传入二值邻接权重后,自适应带宽的搜索范围会被限制在该权重矩阵的非零连接内。如果你的数据存在孤岛(极少数区县没有任何邻接单元),这些样本会直接失去邻居,导致拟合崩溃。处理方法是先把邻接矩阵中度为0的样本排除,或者改用纯距离核。
4. 核心实操:带宽搜索、模型拟合与收敛诊断
4.1 带宽搜索:从原理到参数上限设置的实战策略
带宽搜索是MGWR最耗时也最关键的环节。search_mgwr_paras函数会在指定的带宽候选范围内,通过校正AICc完成每个变量的带宽筛选,搜索策略是网格化穷举加黄金分割优化的组合方式。
两个最需要关注的参数:
bp_tol:搜索收敛容差,默认值通常为0.05,含义是两次迭代之间带宽变化小于5%时认为收敛。追求精度可调低到0.01,但耗时显著增加;search_multi:是否启用并行搜索,多核CPU上设为True可显著提速;max_iter:后向拟合的最大迭代次数,默认200,若模型不收敛需增大。
下面是一段标准调用:
search = search_mgwr_paras( gwr_selector, selector_type='aicc', search_multi=True, max_iter=200, scale=True, )其中scale=True表示在拟合时对解释变量做标准化处理,这会提高数值稳定性,但我前面已经在数据准备阶段手动做过标准化了,这里设True也不会有冲突。最终搜索结果会给出每个变量的最优带宽值,你需要把它打印出来检查合理性:
opt_paras = search[0] print(opt_paras)如果某个变量的最优带宽恰好等于你设置的带宽搜索上限(比如默认等于样本量附近),说明该变量实际上接近全局尺度,在自适应带宽语境下就是"邻居数趋于样本总数"。这是合理结果,不要惊慌,后面我会讲这种带宽在结果解读中的含义。
4.2 模型拟合与关键输出对象解析
带宽搜索完成之后,拟合就很快了。用搜索到的最优带宽数组传入MGWR类即可:
model = MGWR(gwr_selector, sel_opt=opt_paras) results = model.fit()results对象内部包含大量信息,最常见的提取方式如下:
# 每个变量的局部系数,shape为 (n_samples, n_vars) params = results.params # 局部系数的伪t值 t_vals = results.tvalues # 拟合优度信息 aicc = results.aicc r2 = results.r2 adj_r2 = results.adj_r2 # 每个变量的带宽汇总 bw = results.bw这几个字段的用途得讲清楚:
params:每一行是一个样本点,每一列是一个变量的局部回归系数。注意这些系数的量纲对应你传入的因变量量纲和标准化后的解释变量量纲;t_vals:伪t值,绝对值大于1.96(对应95%置信水平约等于正态分布临界值)可认为该样本点上该变量显著。但严格说,MGWR的伪t值并不完全服从标准正态分布,所以更适合作为"探索性显著"的判断依据,而非严格假设检验;bw:每个变量的最终带宽列表,哪怕拟合成功也要重新检查一次。
我还会额外算一个influence信息,用来排查异常样本点。results对象里有get_beta_se()方法可以获取系数标准误,不过实操上更常见的是直接输出各变量的系数分位数,用于后续地图分级配色。
4.3 收敛诊断:不看这个就下结论,迟早出事
MGWR的收敛特性比GWR脆弱。由于是后向拟合迭代算法,变量的初始带宽可能陷入局部最优,特别是变量间相关性较强时。results对象中有个字段iterations(或类似命名,不同版本略有差异),记录了实际迭代次数。我在实践中的判断标准有两条:
- 实际迭代次数明显小于最大迭代上限,说明达到了收敛容忍度;
- 用不同初始带宽重跑一次,最终带宽和系数应该基本一致,如果差异很大,说明陷入了局部最优,需要调整搜索容差或带宽搜索范围。
如果你发现AICc在某次运行中不稳定,或者某个变量带宽反复横跳,优先检查两件事:一是样本点分布是否严重不均造成大量样本邻居数不足;二是是否有多重共线性在局部邻域内爆表。这两个问题的修复都发生在数据准备环节,而不是模型参数环节。
4.4 完整可复现脚本:从数据读入到结果保存
我把前面所有步骤串成一个可直接套用的流水线脚本,注释写得比较详细:
import pandas as pd import geopandas as gpd import numpy as np from mgwr.sel import MGWRSelector, search_mgwr_paras from mgwr.gwr import MGWR # ========== 1. 数据读入与坐标转换 ========== gdf = gpd.read_file("your_data.shp").to_crs(epsg=32650) gdf["X"] = gdf.geometry.centroid.x gdf["Y"] = gdf.geometry.centroid.y df = gdf.drop(columns=["geometry"]) # ========== 2. 变量整理 ========== feature_cols = ["x1", "x2", "x3"] df = df.dropna(subset=["y"] + feature_cols + ["X", "Y"]) for col in feature_cols: df[col] = (df[col] - df[col].mean()) / df[col].std() # ========== 3. 构建选择器 ========== coords = list(zip(df["X"], df["Y"])) y = df["y"].values.reshape(-1, 1) X = df[feature_cols].values selector = MGWRSelector( coords, y, X, kernel="bisquare", fixed=False, sparse=True, ) # ========== 4. 带宽搜索(耗时环节) ========== search_result = search_mgwr_paras( selector, selector_type="aicc", search_multi=True, max_iter=200, scale=True, ) bw_opt = search_result[0] # ========== 5. 拟合与结果输出 ========== model = MGWR(selector, sel_opt=bw_opt) results = model.fit() np.save("mgwr_params.npy", results.params) np.save("mgwr_tvalues.npy", results.tvalues) print("R2:", results.r2, "adj R2:", results.adj_r2) print("bandwidths:", results.bw)这个脚本基本能应对大多数截面空间数据。如果你的数据样本量超过3000,建议search_multi=True加多进程并行,同时在带宽搜索前关掉其他占用内存的程序。
5. 结果解读与可视化:从系数表到空间格局叙事
5.1 带宽结果的第一层解读:每个变量的空间作用尺度
拿到results.bw后,我会先不看系数,先看带宽。因为带宽决定了后续说故事的方向。
在自适应带宽且使用二平方核(bisquare)的前提下,带宽的含义是"每个局部回归实际使用的最近邻样本数量"。举例:某研究有1000个区县,变量A的最优带宽是980,变量B的最优带宽是120。这意味着变量A在估计任意区县的系数时,几乎用到了全部样本,它的空间异质性极弱,系数在全域内接近常数;变量B则在每个回归点只使用最近的120个区县参与加权,空间异质性极强,系数变化剧烈。
把这个转成研究结论:变量A可能是一个结构性、制度性的缓慢空间过程,比如地区人均受教育年限对经济产出的影响,在省域尺度上相对一致;变量B则可能是充满地方色彩的快速过程,比如人口密度对房价的影响,越靠近大城市核心区,影响强度越大。
注意:如果你用的是固定带宽,带宽单位是距离(米或千米),此时解读为"每个局部回归使用方圆多少范围内的样本"。
很多初学者拿到带宽后喜欢用"带宽越小越重要"来解读,这是错误的。带宽大小只反映空间过程尺度,不反映变量重要性。变量重要性要看系数大小、显著样本比例和模型贡献。
5.2 局部系数地图:一张图里怎么放四个变量
结果可视化的核心任务是把每个变量的局部系数空间分布画出来。为了让地图表达可读,我一般会分四步走:
- 把
results.params按样本顺序拼回原始GeoDataFrame; - 对每个系数列做分级(等间隔分5类,或分位数分5类);
- 叠加研究区边界底图;
- 用
matplotlib的subplots多子图排版。
参考代码:
import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap gdf_coef = gdf_proj.copy() for i, col in enumerate(feature_cols): gdf_coef[f"coef_{col}"] = results.params[:, i] fig, axes = plt.subplots(1, len(feature_cols), figsize=(15, 5)) for j, col in enumerate(feature_cols): gdf_coef.plot( column=f"coef_{col}", cmap="RdBu_r", legend=True, ax=axes[j], edgecolor="white", linewidth=0.2, ) axes[j].set_title(f"Coef of {col}", fontsize=11) axes[j].axis("off") plt.tight_layout() plt.savefig("mgwr_coef_maps.png", dpi=300)出图时一个容易破坏专业感的问题是分级配色。RdBu这类发散色带适合系数有正有负的情况,色带中心对应0值;如果系数全部为正或全部为负,用发散色带反而会误导视觉。因此出图前先看一眼系数范围再决定色带类型。
地图表达之外,必须搭配一个统计表:每个变量系数的最小值、最大值、均值、标准差以及显著样本占比。这个表才是审稿人或导师真正看的干货,地图只是直观补充。我通常这样生成:
summary_list = [] for i, col in enumerate(feature_cols): coef = results.params[:, i] tv = results.tvalues[:, i] n_sig = np.sum(np.abs(tv) > 1.96) summary_list.append({ "变量": col, "带宽": results.bw[i], "系数均值": coef.mean(), "系数最小值": coef.min(), "系数最大值": coef.max(), "显著性样本占比": n_sig / len(coef), }) summary_df = pd.DataFrame(summary_list)这个汇总表里的"显著性样本占比"是MGWR结果解读里最有信息量的指标之一。占比高说明该变量的空间异质性交互明显,占比低说明变量在大部分区域都不起作用,只在少数热点地区凸显——这也是多尺度视角才能揭示的结论。
5.3 伪t值与局部显著性:找出真正起作用的子区域
伪t值的可视化方式比系数地图更需谨慎。因为成千上万个样本点的t值地图,如果没有经过统计校正(如FDR校正),会存在大量伪显著点,直接上图会显得整片区域都是"显著红色",信息量反而降低。
mgwr库本身提供了results.tvalues,但校正方法需要自己实现。我常用的做法是先用Benjamini-Hochberg(BH)步骤做FDR校正,再进行显著性二值化(显著=1,不显著=0),然后画"棉花地图"——只有显著区域显示颜色深浅,不显著区域用浅灰填充。这种地图的叙事效果远好于单纯展示t值连续分布。
具体校正逻辑不复杂:把某个变量在全样本上的所有伪t值对应的p值排序,然后按BH公式计算每个p值的显著性阈值,逐一比较判断。对上千个样本,这个过程计算量可以忽略。
5.4 一篇实证论文里怎么系统性呈现MGWR结果
根据我写论文和修改稿子的经验,一套完整且不会挨批的MGWR结果呈现结构大概是这样的三件套:
第一步,带宽结果表:展示每个变量的最优带宽、带宽的绝对值/相对范围、解释变量显著性样本占比。这个表回答的问题是"哪些过程是局部性的,哪些过程是全局性的"。
第二步,系数地图组:把各变量局部系数空间分布展示出来,可以按4到6个等级配色,搭配简洁的图例和气象图形。地图下方用文字标注关键区域——比如"变量A在东部城市群表现出正效应,而在西部农牧区转为负效应"——把读图方向给出来,不能让读者自己瞎看。
第三步,与GWR对比表:把普通GWR的AICc、adj R2和MGWR做对比,展示MGWR的拟合优势。这不是可有可无的环节,而是多尺度方法存在的意义证明。如果MGWR的adj R2只比GWR高0.01,那这篇文章的价值会被质疑。但不是每份数据都能得到大幅提升,遇到这种情况,就重点强调带宽的差异揭示了哪些GWR掩盖掉的尺度结构,这本身就是发现。
6. 我踩过的那些坑和后续可以玩的花活
6.1 几个印象深刻的报错和奇怪现象
用MGWR这几年,我积累了几个频率最高、且教程里不太写的坑,逐个说。
第一个坑:带宽搜索搜索上限紧贴样本量。默认的自适应带宽搜索范围会拉到一个很大的值附近(接近样本总量),如果你的变量确实接近全局尺度,理想带宽就会很大。但如果所有变量的最优带宽都顶到上限,说明这个模型按MGWR应该退化成OLS或者GWR,而不是真的存在多尺度结构。这种情况多半是数据空间异质性很弱,或者变量共线性太强,导致后向拟合无法有效分离各变量尺度。别硬撑着说"多尺度发现了新规律",老老实实把GWR结果也跑出来做对比。
第二个坑:伪t值大面积显著。当样本量超过2000时,伪t值的分布会比标准正态分布更肥尾,直接卡1.96会得到过多显著样本,地图上一片红。你应该改用更保守的阈值,比如|t| > 2.58,或者做FDR校正后再判断显著性。我个人通常以FDR校正后的结果为准,在论文里注明校正方法。
第三个坑:变量符号在大区域内整体反转。这种情况未必是模型算错了,可能是带宽搜索时搜索空间太大导致局部最优,也可能是某个解释变量与其他变量在局部邻域内的相关性结构翻转了。处理方案是先检查局部Pearson相关矩阵,再尝试把初始带宽的范围收紧(通过手动传入搜索候选值来缩小搜索域)。
第四个坑:坐标的投影单位错了还在自我感动。我见过有人用WGS84经纬度直接跑MGWR,带宽结果单位是"度",画系数图时还发现南北部差异巨大——那是因为纬度跨度造成的投影变形。检查办法很简单:打印两个距离较远样本点之间的欧氏距离,看是不是米级。UTM投影下,同一城市内不同街道样本点的距离应该在几百到几万米之间,如果数值是0.000X,那你的坐标没转对。
6.2 进阶方向:时空扩展、异方差调整与可解释性集成
如果你跑通了截面MGWR并觉得不过瘾,有几个自然的后续方向值得我在这里点一下。
第一个方向是时空多尺度地理加权回归(GTWR/MGTWR)。把时间维度纳入带宽搜索,每个变量同时拥有空间带宽和时间带宽,这在面板型空间数据上会派上用场。Python里没有像mgwr一样成熟的MGTWR库,但你可以基于mgwr的框架自己加时间距离矩阵,实现思路不算太远。
第二个方向是异方差和空间自相关同时存在的处理。MGWR假定局部回归残差独立同分布,但真实空间数据往往有残差空间自相关。此时可以考虑将MGWR与空间误差模型(SEM)结合,或者使用稳健标准误。spreg库里提供了一些工具,但是组合使用起来代码要自己写。
第三个方向也最实用:把MGWR结果接入SHAP型可解释性框架。我试过把每个样本点的局部系数当作特征的局部作用值,然后和SHAP值做对比——两者在空间维度上的分布形态往往高度一致,但MGWR提供的带宽尺度信息是SHAP给不了的。将MGWR的"尺度"和机器学习模型的"预测精度"结合,是一个写作空间很大的研究方向。
6.3 最后的实操建议
一定不要把MGWR当成一个按钮工具来用。它适合在论文里做实证深挖,但它的计算过程敏感性强,信任它之前你得先理解带宽搜索的对象和空间权重矩阵的含义。如果你只是项目里临时用一次,那就严格照着我第三、四节的一体化脚本走,数据预处理做好,大概率能跑通。如果是为了学术发表,建议多跑几个不同的核函数(bisquare、gaussian)、不同带宽搜索范围作稳健性检验,并在附录里报告这些检验的对比结果——审稿人看到这一步,基本就不会再质疑你的结果稳健性了。
我自己现在跑一份区县级截面数据,习惯是先用OLS看共线性,再跑GWR看异质性,最后用MGWR看尺度结构,三个模型层层递进,每一步都有明确目的。这套思路走下来,虽然工作量比只跑一个模型多出不少,但写文章时每个结论都有据可依,不再心虚。希望这篇文章能让你的MGWR之路少走几个我最开始走过的弯路。