1. 从全局回归到多尺度地理加权回归为什么MGWR值得学先聊一个我自己的经历。前几年做城市房价影响因素分析用普通最小二乘回归跑出来的结果很规整R方也不低但画完残差空间分布图之后后背发凉残差在空间上明显聚集老城区一片红新区一片蓝。这说明一个根本问题——全局回归假设影响因素对因变量的作用在空间各处是均质的但地理现象几乎从来不是这样。学区对房价的影响在老城区和新区完全不同地铁站的作用在核心区和边缘区也不同。这种空间非平稳性才是地理数据的常态。这时候就该上地理加权回归GWR了。GWR的基本思想不复杂它对每个样本点单独拟合一套系数样本位置不同系数也不同。它解决的核心问题是让模型自己看到空间上下文。系数在哪个区大、在哪个区小就能告诉你这个因素的作用在空间上是如何变化的。但GWR有个隐含假设所有自变量都共享同一种空间尺度也就是同一个带宽。所谓带宽可以理解成回归时使用的局部邻域有多大。GWR用单一带宽去拟合所有变量等于假设房价、地铁距离、学校密度这些因素在影响房价时作用尺度全都一样。这显然不符合真实世界——地铁对房价的影响可能涉及几公里范围而社区绿化率的有效作用半径可能只有几百米。多尺度地理加权回归MGWR就是为了打破这个限制。它允许每个解释变量拥有自己的带宽由数据本身去确定哪个变量在较大尺度上起作用哪个变量在较小尺度上变化。这既保留了GWR的局部拟合能力又能进一步回答变量的空间作用尺度这一层问题。Python生态里用mgwr包实现这条路已经相当成熟而且配合geopandas、matplotlib能做出很漂亮的分析链条。这篇文章我会按实际使用的流程来写原理讲到够用的深度代码直接可用结合作者踩过的坑提醒你关键节点。适合刚接触GWR/MGWR的GIS和数据方向从业者也适合那些已经在ArcGIS里点过GWR工具、想转向Python做更灵活分析的读者。后面涉及的内容较多建议按自己不熟的部分跳读。2. 理解MGWR的精髓带宽、局部回归与后拟合标定2.1 从GWR到MGWR模型形式的差异先看GWR的模型形式这样对比容易理解y_i β0(u_i, v_i) Σ_k βk(u_i, v_i) x_ik ε_i其中(u_i, v_i)是第i个样本点的坐标βk是第k个变量在位置i处的系数。GWR通过在每个样本位置做加权最小二乘给离i越近的样本赋予越高的权重。权重是带宽的函数常见的核函数有bi-square和gaussian两种bi-square在带宽之外权重直接为0gaussian则随距离递减但不会截断。MGWR的公式可以写成y_i β0(u_i, v_i) Σ_k βk(u_i, v_i) x_ik ε_i表面上形式一样区别在背后的标定方式GWR的所有βk共享同一个带宽b而MGWR允许每个βk拥有自己的带宽b_k。也就是说模型会估计出一组带宽可能地铁距离变量的带宽是2.5公里绿化率变量的带宽是800米。这个结果本身就包含信息带宽越大说明该变量在更大的空间范围内保持稳定接近全局性带宽越小说明它只在很小的邻域内变化局部性很强。2.2 多尺度里面的尺度到底是尺度还是带宽很多初学者会把多尺度理解成多分辨率数据或不同图层粒度其实MGWR里的尺度特指空间作用范围用带宽作为实际表现。带宽是核函数中控制权重衰减快慢的参数可以简单理解为一个有效半径。模型在估计某一点的系数时主要参考该半径内样本的信息。这里有个容易误解的点MGWR并没有把空间划分成几个尺度区间而是用连续带宽来表达每个变量的尺度特征。输出结果中每个变量有一个最优带宽数值分析者可以根据这些数值判断变量的空间过程是偏全局还是偏局部。通常如果某个变量的最优带宽接近样本数量或者说接近整个研究区域范围可以理解为该变量在全局尺度上作用系数几乎不因位置变化。在Python的mgwr输出中参数表会列出每个变量的带宽以及标定的尺度比例之类的统计量。利用这个结果你可以构造带宽柱状图直接展示不同影响因素的作用尺度差异这在撰写论文或分析报告时很加分。2.3 后拟合算法MGWR是怎么估计参数的MGWR的估计不像GWR那样可以通过矩阵运算一步得到它的核心是后拟合back-fitting算法。你可以把后拟合理解成逐个变量轮流优化的过程先固定其他变量的系数和带宽只更新第一个变量的系数然后固定其余变量更新第二个变量的系数重复多轮直到所有系数不再明显变化。这种算法源自加性模型的拟合思想好处是能处理多个变量带宽不一致的问题代价是计算量大幅上升而且对初始带宽和初始系数比较敏感。实际操作中MGWR的带宽初值通常由全局GWR的带宽或者用户指定的数值范围来初始化然后通过迭代优化。这也意味着同样的数据不同的初始设置可能得到有细微差异的结果不是bug而是这类算法的内在特点。因此我的建议是不要把MGWR当作一个黑盒跑完直接抄结果。至少要做两轮不同初值的估计对比带宽和系数的稳定性再下结论。3. Python环境搭建与数据准备80%的问题出在这一步3.1 安装mgwr及配套库MGWR在Python中的核心库是mgwr它隶属于pysal地理空间分析库家族依赖libpysal、spglm等模块。我用conda创建独立环境的方式避免污染其他项目conda create -n mgwr python3.9 conda activate mgwr pip install mgwr geopandas shapely matplotlib scikit-learn如果你之前已经装过pysal全家桶那基本只缺mgwrpip install mgwr安装过程中最容易出坑的是libpysal和其它空间库的版本兼容问题。如果import报错提示找不到某个模块优先尝试升级libpysal和spglmpip install --upgrade libpysal spglm esda必须提醒的一点mgwr包目前不算活跃在一些新Python版本上可能存在依赖编译问题用3.9或3.10版本最稳。用conda管理环境不只是为了干净更重要的是好回滚——依赖冲突时直接删掉环境重建比自己折腾系统环境省时得多。3.2 数据格式什么样的数据才能跑MGWRMGWR的输入数据结构说简单也简单说麻烦也麻烦。你需要准备以下几样东西一个二维数组y形状为(n_samples, 1)表示因变量一个二维数组X形状为(n_samples, n_features)表示自变量通常不包含截距项mgwr会自动加截距一个二维坐标数组coords形状为(n_samples, 2)表示每个样本的空间坐标一般为投影坐标系的x、y不要用经纬度除非你的研究范围很小且你清楚投影误差的影响。你需要的是样本级别的数据每一行一个样本每一列一个变量。它不像栅格那样需要邻域结构因为空间权重完全由样本间的距离计算。所以本质上你可以用任何来源的数据只要带上坐标并且做一遍清洗即可。比较省心的做法是直接从GeoDataFrame开始比如读一个Shapefile或GeoJSONimport geopandas as gpd gdf gpd.read_file(your_data.shp) gdf gdf.dropna(subset[price, subway_dist, school_density]) # 计算投影坐标这里选用适合你研究区的投影 gdf gdf.to_crs(EPSG:3857) coords list(zip(gdf.geometry.x, gdf.geometry.y)) y gdf[[price]].values.astype(float) X gdf[[subway_dist, school_density, green_cover]].values.astype(float)注意三个细节一定要剔除含有空值的行。MGWR对缺失值没有自动处理能力NaN进入模型一般会直接报错。使用投影坐标而不是经纬度。因为MGWR的距离计算基于欧氏距离经纬度带来的距离变形在投影坐标系下会小很多。量纲问题最好先做标准化。MGWR的back-fitting过程对变量尺度有一定敏感性虽然理论上模型能自己调节但实践中标准化后的收敛更稳、结果更易解释。用sklearn的StandardScaler处理X即可后面解释系数时再按缩放比例还原。3.3 变量选择与共线性前置检查MGWR比GWR更敏感于解释变量之间的多重共线性。原因在于当带宽不同时某些局部区域可能出现严重的局部共线性也就是变量在该区域内的局部变异高度重叠导致系数出现反常的大正值或负值。这不是MGWR的缺陷而是空间回归中普遍存在的数据病。因此在跑MGWR之前至少做两件事计算变量间的方差膨胀因子VIF把VIF大于10的变量剔除或合并用GWR跑一次基线模型看局部系数的取值范围是否异常宽泛。如果GWR系数在空间上出现剧烈正负交替往往说明共线性或带宽不合理这时候直接跑MGWR只会让问题更隐蔽。我习惯先用相关性热图做快速筛查再上VIF最后用GWR做基线。整个过程半小时能完成但能避免后面很多解释上的尴尬。4. 完整实战用Python跑一个MGWR模型4.1 构造或加载案例数据为了演示完整流程我这边就不贴某个具体城市的真实数据分析结果了而是采用一个结构化的模拟数据集。真实项目里你只需要把数据读取部分替换成自己的业务数据即可。模拟场景是某城市约500个小区的二手房交易数据目标变量是成交单价解释变量包括到地铁站距离、周边学校密度、绿化率和容积率。数据生成时人为设置了不同的作用尺度地铁距离在大尺度上影响房价学校密度在中尺度绿化率在小尺度容积率接近全局。这样可以体现MGWR的优势。import numpy as np import pandas as pd np.random.seed(42) n 500 x_coord np.random.uniform(0, 100, n) y_coord np.random.uniform(0, 100, n) subway_dist np.abs(np.random.normal(50, 20, n)) school_density np.random.gamma(2, 2, n) green_cover np.random.uniform(0, 1, n) floor_rate np.random.normal(2.5, 0.8, n) # 不同变量在不同空间位置上的真实系数 b_subway -0.002 0.0003 * x_coord b_school 0.5 0.2 * np.sin(0.1 * y_coord) b_green 2.0 * (1 - np.abs(x_coord - 50) / 50) b_floor -0.8 price (50 b_subway * subway_dist b_school * school_density b_green * green_cover b_floor * floor_rate np.random.normal(0, 2, n)) df pd.DataFrame({ x: x_coord, y: y_coord, price: price, subway_dist: subway_dist, school_density: school_density, green_cover: green_cover, floor_rate: floor_rate })这段代码只是为了演示真实业务数据至少千级别以上且需要仔细解决采样偏差问题。MGWR对样本覆盖密度较敏感如果某些区域样本过少局部系数估计会变得很不稳定这是后面检查结果时重点要看的。4.2 推荐先从GWR开始同一套数据的基线任何MGWR分析都应该先跑GWR。这句话值得强调。GWR的结果有两个用途一是作为模型对比的基线二是为MGWR提供合理的初始带宽。GWR的运行时间远低于MGWR先跑它成本很低。from mgwr.gwr import GWR from mgwr.sel_bw import Sel_BW y df[[price]].values X df[[subway_dist, school_density, green_cover, floor_rate]].values coords list(zip(df[x], df[y])) # 标准化一下减少量纲影响 from sklearn.preprocessing import StandardScaler X_scaled StandardScaler().fit_transform(X) # 选择GWR带宽 gwr_selector Sel_BW(coords, y, X_scaled, kernelbisquare) gwr_bw gwr_selector.search() print(GWR最优带宽:, gwr_bw) # 拟合GWR模型 gwr_model GWR(coords, y, X_scaled, gwr_bw, kernelbisquare) gwr_results gwr_model.fit() print(gwr_results.summary())这里的关键是带宽选择。sel_bw默认用AICc准则搜索带宽你也可以用交叉验证CV。实际使用中AICc的结果更平滑CV的结果更偏向预测导向两者在数据量较大时差异不大。带宽搜索是一个较耗时的过程尤其样本量上几千之后建议先跑一次小范围搜索比如使用bw参数指定候选范围避免全量网格搜索带来的等待。我还建议你把kernel参数固定在bisquare。因为bi-square核函数在带宽外权重直接归零计算效率和可解释性都更好gaussian核虽然平滑但每个点的计算都要考虑全部样点慢很多而且带宽膨胀时结果难解释。4.3 跑通MGWR模型拟合与常见报错排查接着用GWR的带宽作为MGWR初始带宽的参照跑MGWRfrom mgwr.mgwr import MGWR from mgwr.sel_bw import MGWRSelector # MGWR带宽优化器multi_bw参数用来指定每个变量的初始带宽范围 mgwr_selector MGWRSelector(coords, y, X_scaled, kernelbisquare) # 先用GWR带宽作为所有变量的公共初始带宽 mgwr_selector.multi_bw [gwr_bw] * X_scaled.shape[1] mgwr_bws mgwr_selector.search() print(MGWR各变量带宽:, mgwr_bws) # 拟合MGWR mgwr_model MGWR(coords, y, X_scaled, mgwr_bws, kernelbisquare) mgwr_results mgwr_model.fit() print(mgwr_results.summary())有一个API细节需要注意不同版本的mgwr在MGWRSelector的初始化参数上略有变化有些版本要求传入multi_bw作为初始值有些版本接受默认值再手动赋值。跑之前建议先看一眼函数签名help(MGWRSelector)。我遇到过不少读者按网上旧教程跑结果在MGWRSelector(...)这一步就报错多半就是版本差异。MGWR拟合过程中的常见报错有这么几类我按出现频率排一下输入数组维度不对。y必须是二维形状为(n,1)而不是一维(n,)。X不要包含全为常数的列否则会跟截距项冲突。带宽选择器搜索过程不收敛。通常是因为变量标准化不到位或样本太少。把数据标准化加大迭代上限。内存不足。MGWR在后拟合过程中需要反复计算空间权重矩阵样本超过5000后内存压力陡增。这时可以尝试减少自变量数量或者改用gaussian核配合稀疏矩阵方式处理但mgwr官方支持有限。拟合完成后summary()会输出一大段结果包括每个变量的带宽、局部系数的均值/标准差/最小值/最大值、模型的AICc、R2等。下面围绕怎么读这段输出专门写一节。5. 读懂MGWR结果带宽对比、系数解读与显著性判断5.1 先看带宽再看系数秩序很重要拿到MGWR结果后我建议不要第一时间看系数均值先把带宽表梳理出来。带宽告诉你的是尺度故事哪些变量在多大空间范围内起作用。你可以手工提取每个变量的带宽bw_df pd.DataFrame({ variable: [subway_dist, school_density, green_cover, floor_rate], bandwidth: mgwr_results.bw }) print(bw_df)在这个模拟案例里预期的结果是subway_dist的带宽接近200以上大尺度green_cover带宽较小floor_rate带宽很大甚至接近样本数。这种差异正是MGWR相对GWR的核心增量信息。如果某个变量的带宽小到只有几个样本点你要警惕这个变量在很局部的范围内剧烈变化局部系数估计可能不稳定。如果某个变量的带宽大到覆盖所有样本这个变量基本等价于回归中的全局变量它的系数在空间上变化很小。5.2 系数的空间分布提取并绘制系数地图用MGWR模型可以拿到所有样本点的局部系数形状为(n, k1)第一列是截距项。提取并写回GeoDataFrameparams mgwr_results.params # 形状 (n, 4) tvalues mgwr_results.tvalues gdf[b_subway] params[:, 1] gdf[b_school] params[:, 2] gdf[b_green] params[:, 3] gdf[b_floor] params[:, 4] # 显著性标记一般 |t| 1.96 近似为显著 sig_subway np.abs(tvalues[:, 1]) 1.96 gdf[sig_subway] sig_subway画系数地图时我建议用quantiles分类而不是等间距分类因为局部系数的分布往往偏态等间距会掩盖空间模式。做法很简单import matplotlib.pyplot as plt fig, axes plt.subplots(2, 2, figsize(12, 10)) for ax, col in zip(axes.flatten(), [b_subway, b_school, b_green, b_floor]): gdf.plot(columncol, schemequantiles, k5, cmapRdBu_r, legendTrue, axax) ax.set_title(col) ax.set_axis_off() plt.tight_layout() plt.savefig(mgwr_coef_maps.png, dpi200)画图时还有一个细节如果系数有正有负使用RdBu这类发散色带并设置对称的色标范围不然颜色会误导读者。具体做法是计算全部系数绝对值的最大范围然后设定vmin和vmax为对称值。不同变量使用相同图例范围对比更公平但有时变量系数本身数量级差异很大硬套统一范围会损失细节。我的习惯是分两种情况做展示图时每个变量用各自合适的范围但要在图题注明做严谨对比时使用统一range并配箱线图说明分布。5.3 显著性判断不显著的区域别硬解释MGWR结果里一个常被忽略的点是局部系数的显著性。很多分析者在论文里画出系数变化图就完事却没有区分哪些区域的系数是统计显著的。结论的可靠性会打折扣。mgwr的results对象提供了tvalues属性可以直接用来判断局部系数的显著性。使用方式见上面的代码。绘制时可以将不显著的样本点做空心化处理只填充显著样本点可以更清楚地看出显著性空间结构。另外要注意局部检验存在多重检验问题严格的研究会做FDR校正。虽然这不是必须操作但如果你的结论高度依赖某些区域的显著性建议用p值校正的方法复核一下。简单起见可以直接用statsmodels里的multipletests函数。from statsmodels.stats.multitest import multipletests pvalues 2 * (1 - stats.norm.cdf(np.abs(tvalues[:, 1]))) _, p_corrected, _, _ multipletests(pvalues, methodfdr_bh) gdf[sig_subway_fdr] p_corrected 0.05这样处理后的显著性地图会稳健很多审稿人或业务方的质询也能扛得住。6. 带宽搜索的效率问题样本规模、计算时间与优化策略MGWR最让人烦躁的问题就是慢。数据量到达几千条以后一次完整的MGWR带宽搜索加模型拟合可能从几分钟到几十分钟不等。为了不被漫长的等待消磨耐心有几个实际可用的优化技巧。第一把带宽搜索的候选范围缩小。如果你已经跑了GWR得到全局最优带宽为gwr_bw那么MGWR各个变量的带宽大概率不会偏离GWR带宽太多。在MGWRSelector里传入搜索范围时可以限制在gwr_bw乘以0.5到1.5倍区间。比如希望带宽从20搜索到400可以用range参数或相应的初始化方法限定候选列表。第二先用一小部分样本做一次预演。比如随机抽取20%的样本先跑通流程并得到一个粗略的带宽搜索范围再用全量数据在该范围内精细搜索。这个方法不能直接用于最终结果因为样本数量对带宽有影响但可以有效排查代码问题、估计耗时。第三合理选择核函数。bi-square核让距离超过带宽的样本权重变为0这样在局部拟合时只需要计算带宽内的样本整体开销远低于gaussian核。大量实测下来bi-square的时间优势非常明显且结果差异几乎不影响解释。我在一个约3000条记录的项目里做过粗略测试gaussian核的MGWR带宽搜索耗时接近40分钟而换成bisquare核后直接降到12分钟。如果你的数据还在不断增长强烈建议选择bi-square。如果你确实需要在大规模数据上跑MGWR还可以考虑并行化思路。mgwr自带的多进程支持目前比较有限但带宽搜索本身是串行依赖的不好直接并行。一个可行的变通办法是把研究区域按子区切分每个子区独立建模然后对比子区之间的参数稳定性。这不算标准的MGWR多尺度全局拟合但在业务场景中常常是实用性优先的选择。7. 结果解释中的常见陷阱多尺度不等于多分辨率MGWR给了每个变量一个带宽很容易让人陷入一个解释误区把带宽大小直接等同于影响因素重要性。这两个概念不能混为一谈。带宽描述的只是空间尺度也就是作用范围的大与小重要性描述的是强度比如系数绝对值的大小。一个变量的带宽可以很大但它对因变量的作用可能非常微弱反之一个变量带宽很小局部作用强却只影响一小撮样本。以房价为例到地铁站距离可能拥有很大的带宽表示它的影响范围跨越全市小区物业费可能带宽很小只影响周边一两公里。但二者的系数强度可能完全不同。解释时要分开报告该变量在xx尺度上空间变化系数变化范围为...其中xx区域作用最强。不要笼统说带宽越大越重要。第二个常见坑是对带宽数值的直接外推。MGWR得到的带宽是针对当前样本点的邻域尺度它依赖于样本分布密度。样本点密集的区域局部信息充足带宽可能偏小样本点稀疏的区域为了获得足够的局部估计信息带宽会偏大。也就是说带宽数值里混杂了数据采样密度的影响。跨区域对比不同研究的MGWR带宽时务必谨慎。第三个坑是过度解释系数的局部符号变化。如果一个变量的局部系数在某个区域显著为正在另一个区域显著为负这确实能说明该变量的作用方向发生变化但必须先验证是否由局部共线性导致。一个简单检查方式是对符号变化剧烈的变量绘制与其它变量的相关系数局部地图观察是否在符号变化区域出现相关性突变。我在审阅相关分析报告时还习惯做一个敏感性检查删去最相关的另一个变量后重新跑MGWR看目标变量的带宽和系数图是否发生明显改变。如果明显改变说明结果受到变量组构成的影响结论需要更保守地表达。8. 把MGWR写进报告需要配齐哪些图表和数据到了实际汇报或写论文这一步不少人的模型跑得很好但不知道怎么把结果呈现得让人信服。按照我的经验一份完成度高的MGWR分析报告至少要包含以下内容。首先是GWR和MGWR的对比表。表中列出两者的带宽数、AICc、R2、调整R2、残差平方和以及各自的最优带宽范围。这个对比的核心价值是说明我们用了多尺度模型而且它确实比单尺度GWR更好。其次是带宽条形图展示每个变量的带宽并按从小到大排序。带宽条形图能让读者一眼看出哪些变量是全局性的哪些是局部性的。如果有的变量带宽趋近样本总量建议在图中特别标注近似全局变量。第三是系数空间分布图。每个解释变量一幅图按分位数划分色级叠加上显著性标记。图名不要直接用变量代码要写业务含义比如到地铁站距离的影响系数。最后建议附一段结果解读文字。避免罗列统计量而是讲故事模型发现哪些变量在局地尺度上变化、哪些变量相对稳定、显著区域的空间格局意味着什么。一个好的MGWR分析并不止于得出系数在空间上变化更要回答这种变化背后的城市规划/交通/人口结构因素是什么。我习惯把报告的结论部分限定在从数据中看到的现象和需要结合专业知识解释的现象两层不强行给出因果判断。空间回归模型说明的是空间相关模式不是因果机制。这一点在公开写作时特别重要。9. 写在后面先用GWR打底再让MGWR说话跑过几十次MGWR之后我最深的体会是MGWR不是替代GWR的神器而是GWR的升级补充。它更适合用来回答除了空间非平稳性是否存在多尺度效应这个问题。如果你只是想得到一套平滑的局部系数地图GWR已经够用如果想把不同变量的作用尺度差异也纳入结论才轮到MGWR登场。实际操作中我的标准流程始终是OLS打底、GWR优化、MGWR深化。先看OLS残差有没有空间自相关再跑GWR确定是否存在空间非平稳性最后用MGWR解释这种非平稳性是否来源于多尺度过程。每一层模型都在回答不同层次的问题这样出来的结论逻辑完整也经得起业务或学术上的质疑。一个小技巧作为收尾跑MGWR时把每次带宽搜索和模型拟合的中间结果都保留下来特别是GWR的带宽和AICc。因为不同项目之间需要反复对比这些中间数据能帮你在模型结果异常时快速定位到底是数据问题还是初始参数问题。毕竟地理加权回归这类模型数据的质量和模型初值的选择往往比算法本身的差异更能影响最终结论。