利用Python构建空间滞后模型分析网约车出行量影响因素

在毕业论文写作之余，为大家带来一篇干货文章。通常我们在做空间计量模型的时候，最常用的应该是Geoda（然而我并不会用），但又想学习一下怎么办？Pysal帮我们解决了这个大问题，Pysal集成了众多空间分析与空间计量的功能，十分强大，其中的子模块spreg即可实现众多的空间计量模型，对于空间计量经济学方面的理论我一无所知，因此本篇文章是从GIS及数据分析角度出发的，如果有错误的话，还请各位老师同学多多指出，欢迎私信公众号后台与我交流。

---

正常针对数据建模的话，我们最常用的便是基于最小二乘OLS的多重线性回归。多重线性回归假设样本之间都是相互独立的。而对于空间数据，这一假设常常是不成立的。举个例子，如果我想研究网约车出行量与POI数量之间的关系，把交通小区内网约车出行量设置为y，把各类POI个数设置为x，然后做OLS回归，其实这时我们默认了交通小区之间都是相互独立的。但是，实际上，根据地理学第一定律，相邻的交通小区之间，必然存在相关关系，比如，网约车出行量大的小区，通常是扎堆的，一个小区的出行量与其相邻小区的出行量通常是相关的，这就是空间依赖(spatial dependence)或空间自相关(spatial autocorrection)。OLS忽略了这种空间依赖，会导致模型设定错误，造成错误或有偏的结果。通过空间计量模型，我们可以捕捉空间依赖的结构，将其放入方程中，使模型更为准确。常用的空间计量模型包括空间滞后模型(SLM)和空间误差模型(SEM)，前者假设空间依赖存在于因变量中，后者假设空间依赖存在于误差项（不可观测的因素）中。下面上代码，给大家介绍一下应用空间计量模型的基本流程。用到的数据还是GWR那篇文章中的数据，感谢滴滴盖亚计划。

首先，我们对数据进行处理，具体的方法在前几期都讲过了，这里不展开。本次分析是以成都市的街道作为分析单位的，共有90个街道，也就是说有90个样本。处理过后的数据长这样：

这张表就是每个街道包括的各类POI数量和网约车出行量，顺便说一下，滴滴盖亚提供的数据1天成都大概有20w的订单，我觉得如果是全样本的话，似乎有点少。不过我看东南有人拿这个成都的盖亚发了Journal of Transport Geography，用1天20w的数据研究整个成都市域网约车出行量，具体怎么回事咱也不清楚，就先假设这20w是全部的吧。

首先我们算一下moran指数，用的也是GWR那篇文章里的方法。

自相关还是挺明显的，都显著。

然后画一下散点图：

from libpysal.weights import lag_spatial
lag_ridership = lag_spatial(w_queen,stat['ridership'])
plt.figure(figsize=(8,6))
plt.plot(stat['ridership'],lag_ridership,'.',color='red')
plt.vlines(stat['ridership'].mean(),lag_ridership.min(),lag_ridership.max(),linestyle='--',alpha=0.8)
plt.hlines(lag_ridership.mean(),stat['ridership'].min(),stat['ridership'].max(),linestyle='--',alpha=0.8)
plt.title('Moran Scatterplot',fontsize=15)
plt.xlabel('Ridership',fontsize=13)
plt.ylabel('Spatial lag of ridership',fontsize=13）

出行量的HH LL集聚现象还是挺明显的。

在建模前，先看一下变量之间的共线性：

var.corr()

把相关系数大于0.7的去掉避免共线性，最后保留以下三个：

然后开始OLS建模：

ols = spreg.OLS(y,x,w_queen,spat_diag=True,moran=True,name_y='ridership',name_x=list(var_3.columns),name_ds='DIDI_chengdu',name_w='queen')
print(ols.summary)

结果的话，我们首先看一下空间依赖性的诊断，残差空间分布的moran指数显著为正，说明这里使用OLS是有很大问题的（模型如果设置正确的话，残差空间分布应该是随机的），LMlag、LMerror和RLMlag检验都显著，RLMerror不显著，说明存在因变量的空间依赖，因此我们需要使用空间滞后模型(SLM)。

针对OLS残差的分析的话，正常到这里就结束了，不过我们可以玩的更酷一点，在这里我模仿一下杨老师新发的Applied Geography中的残差可视化方法，画一个残差空间分布图：

杨老师原版（高端大气上档次，第一看的时候内心os（这也太好看太高级了吧））：

我模仿的：

base = chengdu.plot(figsize=(12,12),facecolor='none',edgecolor='grey')
for i in stat.index:
    if stat.loc[i,'ols_res'] > 0:
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='royalblue',alpha=0.6,s=abs(stat.loc[i,'ols_res'])/2)
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='',edgecolor='black',s=abs(stat.loc[i,'ols_res'])/2)
    else:
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='gold',alpha=0.6,s=abs(stat.loc[i,'ols_res'])/2)
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='',edgecolor='black',s=abs(stat.loc[i,'ols_res'])/2)
plt.title('Residuals of OLS',fontsize=20)
plt.axis('off')
plt.show()

哈哈哈，精髓只学到一半（这里说一下，配色分别是gold和royalblue，看这个名字就特别高级）。看这个残差，确实是扎堆了，存在明显的空间自相关。

前文的各种检验结果告诉我们，空间滞后模型更为恰当，下面我们开始建模：

splag = spreg.ML_Lag(y,x,w_queen,name_y='ridership',name_x=list(var_3.columns),name_ds='DIDI_chengdu',name_w='queen')
print(splag.summary)
lr = spreg.diagnostics.likratiotest(ols,splag)

这里我用的是Queen矩阵，结果如下：

可以看到，在考虑了空间依赖后，只有residential的poi是显著的，此外空间滞后项的系数显著为正，说明出行量高的小区周围的小区出行量也高。似然比检验的结果说明，SLM拟合优度比OLS有了显著提升。

然后让我们看一下，用了SLM之后，模型的残差的空间自相关情况与空间分布：

可以看到，moran指数并不显著，说明SLM的残差并没有明显的空间自相关存在，该模型设定是正确的。

base = chengdu.plot(figsize=(12,12),facecolor='none',edgecolor='grey')
for i in stat.index:
    if stat.loc[i,'lag_res'] > 0:
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='royalblue',alpha=0.6,s=abs(stat.loc[i,'lag_res'])/2)
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='',edgecolor='black',s=abs(stat.loc[i,'lag_res'])/2)
    else:
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='gold',alpha=0.6,s=abs(stat.loc[i,'lag_res'])/2)
        plt.scatter(stat.loc[i,'x'],stat.loc[i,'y'],marker='o',color='',edgecolor='black',s=abs(stat.loc[i,'lag_res'])/2)
plt.title('Residuals of SLM',fontsize=20)
plt.axis('off')
plt.show()

空间分布长这样，可以看到，比OLS残差的空间分布好了很多，扎堆的现象不再那么明显。

至此，本篇文章就结束啦，如果需要全部代码的话可以后台私信我。欢迎各位老板猛烈点赞猛烈转发，如果再能打赏的话就再好不过了，我们下期见！

往期内容