【问题导向】GWR与MGWR——以南京市中心城区住宅小区为例

背景

大二的时候和导师进行过GWR(地理加权回归)的学习，当时只是囫囵吞枣，对于GWR的含义和细节并没有很好的理解。这几天在微信公众号看了城市数据派，古恒宇博士的讲座。讲座链接.

抽时间实现了GWR和MGWR，温故而知新，本次学习中有以下几点新收获：
1，如何理解GWR，带宽、核函数是什么。MGWR是如何优化GWR的不足。
2，如何解读GWR和MGWR的结果。
3，如何使用python读取和绘制shp文件（geopandas库的使用）

实例：南京市中心城区住宅小区分析

数据预处理

研究区域为南京市的主城区（玄武区、秦淮区、建邺区、鼓楼区、栖霞区、雨花台区）

从高地地图爬取的POI数据，然后利用Arcgis进行预处理，将研究区域用1km*1km格网划分，统计每个格网内的POI数据量。为了避免多重共线性问题，去除掉小区数为0的格网。

之后的处理都使用python执行。本次处理是第一次使用一直心念念的geopandas库。值得纪念。
geopandas是用来处理地理空间数据的python第三方库，它是在pandas的基础上建立的，完美地融合了pandas的数据类型，并且提供了操作地理空间数据的高级接口，使得在python中进行GIS操作变成可能。
执行MGWR时spyder会报错（建模时有一段动画），因此要使用Jupyter Notebook运行代码。

首先使用geopandas读取库，并可视化。

import numpy as np
from mgwr.gwr import GWR, MGWR
from mgwr.sel_bw import Sel_BW
import geopandas as gp
import matplotlib
import matplotlib.pyplot as plt
import time
import contextily as ctx
from mgwr.utils import shift_colormap, truncate_colormap
#导入shp数据，可视化并加载在线底图
georgia_shp = gp.read_file('D:\data\data.shp').to_crs('EPSG:3857')#变换投影为web墨卡托即EPSG:3857
fig, ax = plt.subplots(figsize=(10,10))
ax = georgia_shp.plot(ax=ax, alpha=0.1, edgecolor='k')
ax.axis('off')
georgia_shp.plot(ax=ax, alpha=0.5,**{'edgecolor':'black', 'facecolor':'white'})#alpha设置透明度
ctx.add_basemap(ax, 
                source='https://{s}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png',
                zoom=11)

加载查看一下数据，可以看见geopandas在pandas最后加了一列geopandas对应点线面形状。

构建模型

读取自变量，因变量，坐标等准备拟合。

#准备拟合参数
g_y = georgia_shp['小区'].values.reshape((-1,1))
georgia_shp['政府机关']=georgia_shp['政府机']
georgia_shp['休闲娱乐']=georgia_shp['休闲娱']
g_X = georgia_shp[['政府机关','休闲娱乐','零售','医疗','公园','餐饮','酒店']].values
u = georgia_shp['pointx']
v = georgia_shp['pointy']
g_coords = list(zip(u,v))
g_X = (g_X - g_X.mean(axis=0)) / g_X.std(axis=0)
g_y = g_y.reshape((-1,1))
g_y = (g_y - g_y.mean(axis=0)) / g_y.std(axis=0)

执行GWR,带宽唯一为59.0,耗时为0.5271022319793701

#执行GWR
time_start=time.time()
gwr_selector = Sel_BW(g_coords, g_y, g_X)
gwr_bw = gwr_selector.search()
print(gwr_bw)
gwr_results = GWR(g_coords, g_y, g_X, gwr_bw).fit()
time_end=time.time()
print('GWR cost',time_end-time_start)

执行MGWR,带宽为43；184；271；51；271；63；271；74；耗时为112.7812876701355；
可见MGWR相对于GWR能够对每个变量构建带宽，但同时更耗费时间。

#执行MGWR
time_start=time.time()
mgwr_selector = Sel_BW(g_coords, g_y, g_X, multi=True)
mgwr_bw = mgwr_selector.search()
print(mgwr_bw)
mgwr_results = MGWR(g_coords, g_y, g_X, mgwr_selector).fit()
time_end=time.time()
print('MGWR cost',time_end-time_start)

将模型计算的系数添加到shp数据中

georgia_shp['gwr_intercept'] = gwr_results.params[:,0]
georgia_shp['gwr_政府机关'] = gwr_results.params[:,1]
georgia_shp['gwr_休闲娱乐'] = gwr_results.params[:,2]
georgia_shp['gwr_零售'] = gwr_results.params[:,3]
georgia_shp['gwr_医疗'] = gwr_results.params[:,4]
georgia_shp['gwr_公园'] = gwr_results.params[:,5]
georgia_shp['gwr_餐饮'] = gwr_results.params[:,6]
georgia_shp['gwr_酒店'] = gwr_results.params[:,7]
gwr_filtered_t = gwr_results.filter_tvals()

georgia_shp['mgwr_intercept'] = gwr_results.params[:,0]
georgia_shp['mgwr_政府机关'] = gwr_results.params[:,1]
georgia_shp['mgwr_休闲娱乐'] = gwr_results.params[:,2]
georgia_shp['mgwr_零售'] = gwr_results.params[:,3]
georgia_shp['mgwr_医疗'] = gwr_results.params[:,4]
georgia_shp['mgwr_公园'] = gwr_results.params[:,5]
georgia_shp['mgwr_餐饮'] = gwr_results.params[:,6]
georgia_shp['mgwr_酒店'] = gwr_results.params[:,7]
mgwr_filtered_t = mgwr_results.filter_tvals()

对结果可视化

def deal(a,num):
    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(45,20))
    myfont = matplotlib.font_manager.FontProperties(fname="simsun.ttc")
    ax0 = axes[0]
    ax0.set_title('GWR '+a+' (BW: ' + str(gwr_bw) +')', fontsize=40,fontproperties=myfont)
    ax1 = axes[1]
    ax1.set_title('MGWR '+a+' (BW: ' + str(mgwr_bw[num]) +')', fontsize=40,fontproperties=myfont)

    #Set color map
    cmap = plt.cm.seismic

    #Find min and max values of the two combined datasets
    gwr_min = georgia_shp['gwr_'+a].min()
    gwr_max = georgia_shp['gwr_'+a].max()
    mgwr_min = georgia_shp['mgwr_'+a].min()
    mgwr_max = georgia_shp['mgwr_'+a].max()
    vmin = np.min([gwr_min, mgwr_min])
    vmax = np.max([gwr_max, mgwr_max])

    #If all values are negative use the negative half of the colormap
    if (vmin < 0) & (vmax < 0):
        cmap = truncate_colormap(cmap, 0.0, 0.5)
    #If all values are positive use the positive half of the colormap
    elif (vmin > 0) & (vmax > 0):
        cmap = truncate_colormap(cmap, 0.5, 1.0)
    #Otherwise, there are positive and negative values so the colormap so zero is the midpoint
    else:
        cmap = shift_colormap(cmap, start=0.0, midpoint=1 - vmax/(vmax + abs(vmin)), stop=1.)

    #Create scalar mappable for colorbar and stretch colormap across range of data values
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))

    #Plot GWR parameters
    georgia_shp.plot('gwr_'+a, cmap=sm.cmap, ax=ax0, vmin=vmin, vmax=vmax, **{'edgecolor':'black', 'alpha':.65})
    #If there are insignificnt parameters plot gray polygons over them
    if (gwr_filtered_t[:,0] == 0).any():
        georgia_shp[gwr_filtered_t[:,0] == 0].plot(color='lightgrey', ax=ax0, **{'edgecolor':'black'})

    #Plot MGWR parameters
    georgia_shp.plot('mgwr_'+a, cmap=sm.cmap, ax=ax1, vmin=vmin, vmax=vmax, **{'edgecolor':'black', 'alpha':.65})
    #If there are insignificnt parameters plot gray polygons over them
    if (mgwr_filtered_t[:,0] == 0).any():
        georgia_shp[mgwr_filtered_t[:,0] == 0].plot(color='lightgrey', ax=ax1, **{'edgecolor':'black'})
 
    #Set figure options and plot 
    fig.tight_layout()    
    fig.subplots_adjust(right=0.9)
    cax = fig.add_axes([0.92, 0.14, 0.03, 0.75])
    sm._A = []
    cbar = fig.colorbar(sm, cax=cax)
    cbar.ax.tick_params(labelsize=50) 
    ax0.get_xaxis().set_visible(False)
    ax0.get_yaxis().set_visible(False)
    ax1.get_xaxis().set_visible(False)
    ax1.get_yaxis().set_visible(False)
    plt.show()

num=0
for i in ['政府机关','休闲娱乐','零售','医疗','公园','餐饮','酒店']:
    deal(i,num)
    num=num+1

模型对比

从多个指标分析，可以看到MGWR优于GWR

结果分析

可以看见MGWR不同变量的带宽不同。比如零售、公园的带宽小，说明影响距离小，服务范围小，这也符合我们生活常识，去零售店或者公园都是就近原则。而餐饮、医疗、娱乐设施的影响距离远，说明服务范围大，吃饭就医远一些也能忍受。

数据较多，只挑选几个典型的分析。

休闲娱乐设施显著正向影响小区数量，说明娱乐设施越多居住环境越好，小区数量越多。

零售设施正向影响小区数量，且具有圈层结构。说明零售施能显著提升居住质量。但是距离影响较少，稍微远一些影响就呈现负效果。

同样公园对居住的影响也是圈层，但是与零售完全相反，且影响力很小。可能是公园的布局与零售不同。

理解GWR和MGWR

空间统计与空间过程

ps：后面的内容是结合古恒宇博士的讲座总结而来。讲座链接.

传统统计学中，统计数据时没有考虑到位置信息。如我们中学就学过的OLS。随着人们将数据的地理位置信息纳入考虑，空间统计学也发展而来。

空间数据是一种结果，我们需要对形成这种结果的空间过程进行一定的假设和分析；主要的分析方向是空间自相关（空间依赖性，地理学第一定律）和空间异质性（地理学第二定律，地理数据在空间上的差异性）。而GWR主要是围绕空间异质性展开。

如何对空间过程描叙？传统上是用过感性的语言进行描述，比如我们中学时做地理分析题的那种术语。然而随着数学引入地理学中，现代地理学对过程的描叙需要用严谨的数学公式进行。（左边的随机性大一些，右边的随机性小一些）

GWR更关注的特性是2，过程的异质性，也就是对关系的异质性的描述（A、B关系的异质性，A、B的异质性是关系异质性的结果）

GWR的理解

GWR中的空间尺度（也可以理解为带宽），描述了各个样本的空间关系在多大空间范围内是平稳的。

GWR将位置引入为回归参数，用空间权重来表达，一般是样本回归点到其他各点的距离衰减函数。
GWR能够很直观的表现、处理空间异质性。允许全局空间的变参数。

GWR的理解：
中间的红点，代表进行回归的点。只有一定距离（带宽，即以任意一点为中心，邻域的范围）的样本点参与回归，距离较远的点影响就趋近于0了。

彩色的“罩子”理解为核函数，他的范围就是带宽。通过核函数确定邻域对该点影响的大小即权重。（根据地理学第一定律，越近的数据点的权重越高）

通过回归就能可视化这种影响的变化。

求带宽的方式很多，如固定带宽，自适应带宽；带宽越大纳入的点越多。纳入所有的点，就是全局回归。

自适应带宽

GWR的不足

GWR的带宽是固定的，比如回归中有多个参数（x1 x2 x3…）；在GWR中这些参数的带宽都一样。如下图。

但是，在实际中不同的参数的影响是不一样的。如对于房价而言面积和教育的影响程度是不同的（学区房面积虽然小，但是价格仍然很高）。因此，对不同的参数，带宽应该是不同的。

MGWR解决了这一点，MGWR（多尺度地理加权回归）中不同参数的带宽是不同的，这比GWR描述异质性上进了一步。如下图

MGWR

MGWR的优势

推荐的文献：
第一篇提及GWR的文章

讨论权重求解问题

应用

分析方法

一般是对同一组数据分别进行OLS GWR MGWR三次实验，观察效果。
衡量不同模型精度关注的是：$R^2$、Sigma、AICc（赤则信息量，越低越好）、BIC（贝叶斯信息量，越小越好）
对GWR和MGWR应该分析每个变量的系数。分析改变量是正向影响还是负向影响y

对MGWR分析应该包括不同变量带宽的差异。

链接:
1, Jiang, W.; Wang, Y.; Dou, M.; Liu, S.; Shao, S.; Liu, H. Solving Competitive Location Problems with Social Media Data Based on Customers’ Local Sensitivities. ISPRS Int. J. Geo-Inf. 2019, 8, 202. https://doi.org/10.3390/ijgi8050202.

2, 模型选择的一些基本思想和方法.

3, 白话空间统计.

4, 沈体雁,于瀚辰,周麟,古恒宇,何泓浩.北京市二手住宅价格影响机制——基于多尺度地理加权回归模型(MGWR)的研究[J].经济地理,2020,40(03):75-83..

新开通了本人的公众号，欢迎关注：燕南路GISer ，专注GIS干货分享，不定期更新。
主要兴趣：GIS、时空数据挖掘、python、机器学习深度学习
优快云的部分内容会重写再搬迁到公众号，欢迎关注！