空间滞后–混合地理加权回归模型中的数据分区及参数估计

Data Partition and Parameter Estimation in Spatial Lag-Mixed Geographical Weighted Regression Model

摘要: 空间滞后模型和地理加权回归模型均为经典的地统计学模型，分别用于处理带有空间自相关性或异质性的数据，但是在处理同时带有空间自相关性和异质性的数据时拟合效果较差。为了同时考虑数据的自相关性和异质性，提升模型的拟合效果，本文在空间滞后模型和地理加权回归模型的基础上做出改进。首先针对空间数据的异质性，使用改进的k均值聚类方法对空间数据进行分区处理。其次，在分区内部引入空间的自相关性，给出空间滞后–混合地理加权回归模型，并提出了基于莫兰指数与权重矩阵的关系进行估计的莫兰指数优化法。通过在真实数据集上的实验研究，证明了本文方法相比传统方法具有更好的拟合效果。

Abstract: Both the spatial lag model and the geographically weighted regression model are classically geostatistical models, which are used to deal with data with spatial autocorrelation or heterogeneity respectively, but the fitting effect is poor when dealing with data with both spatial autocorrelation and heterogeneity. In order to consider the autocorrelation and heterogeneity of the data at the same time and improve the fitting effect of the model, this paper makes improvements on the basis of the spatial lag model and the geographically weighted regression model. Firstly, according to the heterogeneity of spatial data, the improved k-means clustering method is used to partition the spatial data. Secondly, the spatial autocorrelation is introduced into the interior of the zone, and the spatial lag-mixed geographically weighted regression model is given, and the Moran’s I optimization method based on the relationship between Moran’s I and weighted matrix is proposed. Through experimental research on real data sets, it is proved that this method has better fitting effect than traditional methods.

1. 引言

普通线性回归模型用于处理独立、正态和方差齐性的数据，但是由于空间数据存在自相关性和异质性，很难满足独、正态、方差齐性这些假设条件，因此需要构建适用于空间数据的空间回归模型。

空间回归模型包括空间滞后模型、地理加权回归模型等，主要用于处理空间数据，空间数据的两个重要特征分别是空间自相关性和空间异质性。针对空间自相关性的研究最初是在线性回归模型中加入含有空间权重矩阵的项，使得研究区域的因变量不仅与本区域的解释变量有关，还与邻近区域的因变量相关。1988年Anselin [1] 提出了空间线性回归(Spatial Linear Regression, SLR)模型，当模型中的参数取值不同时，SLR模型可派生出普通线性回归模型、一阶段空间自回归模型、空间滞后模型(Spatial Lag Model, SLM)、空间误差模型以及空间杜宾模型。其中空间滞后模型可以用来分析许多具有空间自相关性的问题，如房价预测，犯罪率分析，经济增长等。

针对空间数据的异质性，1996年Brunsdon等人 [2] 首先将异质性的思想融入到SLM的参数求解中，提出局部最大似然估计法，将SLM中所有的空间系数看做随空间位置变化的变量。1998年Fotheringham等人 [3] 将SLM中的常系数变为随空间位置变化的变系数来处理空间数据的异质性，提出著名的地理加权回归(Geographically Weighted Regression, GWR)模型。GWR模型将所有解释变量作为局部解释变量，由于局部解释变量既可能与空间位置有关，也可能与空间位置无关，因此2011年Paez等人 [4] 将GWR模型中部分解释变量的系数变为常系数，提出混合地理加权回归(Mixed Geographically Weighted Regression, MGWR)模型，使得模型拟合效果相比于GWR模型更精确。目前，混合地理加权回归模型在空间统计学领域已经得到了广泛应用，可以用于城市规划、环境保护和资源管理等方面，通过精确预测，为决策者提供更准确的决策依据。

MGWR模型以及SLM有利有弊：SLM考虑了因变量间的空间自相关性，却忽略了空间异质性；MGWR模型虽然基于局部光滑的思想减弱了空间异质性，但是却没有考虑空间自相关性。目前，同时考虑空间自相关性和空间异质性模型的研究寥寥无几，相应的参数估计以及模型的拟合效果检验更是凤毛麟角。1996年Brunsdon等人 [5] 在考虑空间自相关的情况下，使用极大似然估计法对GWR模型进行参数估计，在每个观测点处都进行线性回归计算。2005年魏传华等人 [6] 通过模拟实验验证了具有空间自相关的GWR模型中常系数估计的精确性和稳健性。2011年Geniaux等人 [7] 证明了在MGWR模型中，空间部分回归系数是非平稳的，并采用局部两步最小二乘法对MGWR模型进行参数估计。2013年乔宁宁 [8] 提出了混合地理加权空间滞后回归模型，给出相应的参数估计方法，即先由局部最小二乘法得到变系数的表达式，再使用极大似然估计的方法对参数进行求解。本文在2017年乔宁宁 [8] 的研究基础上，通过在MGWR模型上增加空间滞后项来降低空间自相关性对结果的影响，提出空间滞后–混合地理加权回归(Spatial Lag-Geographical Weighted Regression, SL-MGWR)模型，SL-MGWR模型的预测能力和拟合效果均优于原始模型。针对SL-MGWR模型的参数估计问题，将混合地理加权回归模型的两步估计法 [7] [8] 扩展到空间滞后–混合地理加权回归模型中，提出了莫兰指数优化法。莫兰指数优化法能够有效降低计算复杂度，提高模型预测能力并降低了自相关性和异质性对模型计算结果的影响程度。SL-MGWR模型可以处理原始模型的所有问题，既可以应用在房价预测，犯罪率分析，经济增长，又可以用于城市规划、环境保护和资源管理等方面，且其解释能力及预测效果都优于原始模型。

本文第一部分是引言，介绍本文研究背景。第二部分是模型介绍，介绍了地理加权回归模型、空间滞后模型以及本文提出的空间滞后–混合地理加权回归(SL-MGWR)模型。第三部分是参数估计方法介绍，介绍了两步估计法和莫兰指数优化法的推导过程。第四部分是实验研究，分为三小节，第一小节为使用两步估计法 [7] [8] 比较普通地理加权回归模型以及分区后的地理加权回归模型，来说明数据分区之后模型的模型效果更好，第二小节为使用两步估计法将GWR以及SL-MGWR模型进行对比，说明SL-MGWR模型的拟合效果较好。第三小节是SL-MGWR模型传统及莫兰指数优化法的拟合效果对比，得到莫兰指数优化法的效果优于传统两步估计法。第五部分是结论，得到最后的结果。

2. 模型介绍

2.1. 地理加权回归模型

由于空间数据存在空间异质性，空间地区内的空间单元本身不是均质的，在形状及面积上有较大差别，而在进行空间数据分析时，变量的观测值一般都是以给定的空间单元为抽样单位得到的，随着地理位置的变化，变量间的关系或者结构也会引起变化。在全局模型中，我们假设不存在空间异质性，所得结果有一定误差，所以需要对传统方法进行改进。

1996年Brunsdon [2] 及1998年Fotheringham等人 [3] 基于局部平滑的思想，提出了地理加权回归模型，将数据的空间位置嵌入到回归参数中，利用局部加权最小二乘法进行逐点参数估计。即，

$y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\beta }_k\left(u_i,v_i\right)x_{ik}}+{\epsilon }_i$

其中 $\left(u_i,v_i\right)$ 为第i个采样点的坐标， ${\beta }_k\left(u_i,v_i\right)$ 为第i个采样点上的第k个回归参数，是地理位置的函数； ${\epsilon }_i$ 是第i个区域的随机误差，满足零均值、同方差、相互独立等基本假定。

2.2. 空间滞后模型

传统经典的线性回归模型具有严格的前提假设条件——独立、正态、齐方差性。但是，由于空间数据存在空间自相关性和异质性，使得这些条件很难满足，在使用传统线性回归模型解决空间问题时，会造成模型参数错误并降低模型的有效性。1988年Anselin [1] 考虑空间数据的自相关性，给出了如下形式，使用空间滞后模型来处理数据，

$Y=\rho WY+X\beta +\epsilon ,\epsilon ~N\left[O,{\sigma }^{2}I\right].$

其中Y为 $n\times 1$ 维响应变量 $Y=\left(Y_1,Y_2,\cdots ,Y_n\right)$ ， $X=\left(X_1,X_2,\cdots ,X_K\right)$ 为包含K个解释变量的 $n\times K$ 维矩阵。 $W=\left[\begin{array}{cccc}0& b_{12}& \cdots & b_{1n}\\ b_{21}& 0& \cdots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ b_{n1}& b_{n2}& \cdots & 0\end{array}\right]$ 为空间权重矩阵， $b_{ij}=\frac{1}{d_{ij}^2}$ ， $d_{ij}$ 为位置i和位置j之间的距离，且 $d_{ij}=\sqrt{{\left(u_i-u_j\right)}^2+{\left(v_i-v_j\right)}^2}$ ，其中 $u\left(i\right),v\left(i\right)$ 分别为位置i的经度和纬度。WY为空间滞后因子，空间滞后

考虑一个特定的观测地区会受到相邻地区的影响，使用空间滞后因子可以得出实际上相邻的空间地区观测值依距离加权后的平均值， $\rho \in \left[0,1\right]$ 为空间滞后项系数，取值越接近1，说明相邻地区的因变量取值越相似。 $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_K\right)}^{\text{T}}$ 为参数向量， $\epsilon ={\left[{\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_n\right]}^{\text{T}}$ 为服从正态分布的n维随机误差向量且 $Cov\left({\epsilon }_i,{\epsilon }_j\right)=0,i\ne j$ ，I为 $n\times n$ 的单位矩阵。

2.3. 空间滞后–混合地理加权回归模型

为了同时考虑空间数据的异质性及自相关性，本文将空间滞后模型中的空间滞后因子加入地理加权回归模型中，将空间数据的自相关性和异质性结合起来，得到具有空间自相关性的地理加权回归模型：

$y_i=\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{b}_{ij}}{y}_j+{\beta }_0\left(u_i,v_i\right)+{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\beta }_k\left(u_i,v_i\right)x_{ik}}+{\displaystyle \underset{k=p+1}{\overset{m}{\sum }}{\beta }_k{x}_{ik}}+{\epsilon }_i$

为避免混淆，假设总共有n个样本点，模型中具有m个解释变量，令 ${\beta }_0\left(u_i,v_i\right)={\beta }_{i0},{\beta }_k\left(u_i,v_i\right)={\beta }_{ik}$ 。可得：

$\begin{array}{c}{y}_i=\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{b}_{ij}}{y}_i+{\beta }_0\left(u_i,v_i\right)+{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\beta }_k\left(u_i,v_i\right)x_{ik}}+{\displaystyle \underset{k=p+1}{\overset{m}{\sum }}{\beta }_k{x}_{ik}}+{\epsilon }_i\\ =\rho {\displaystyle \underset{j=1}{\overset{n}{\sum }}{b}_{ij}}{y}_i+{\beta }_{i0}+{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\beta }_{ik}}{x}_{ik}+{\displaystyle \underset{k=p+1}{\overset{m}{\sum }}{\beta }_k{x}_{ik}}+{\epsilon }_i\end{array}$

可以简写为：

$Y=\rho WY+X\beta +M+\epsilon .$ (1)

其中 $M={\left[\begin{array}{cccc}{\displaystyle \underset{k=0}{\overset{p}{\sum }}{\beta }_{1k}{x}_{1k}}& {\displaystyle \underset{k=0}{\overset{p}{\sum }}{\beta }_{2k}{x}_{2k}}& \cdots & {\displaystyle \underset{k=0}{\overset{p}{\sum }}{\beta }_{nk}{x}_{nk}}\end{array}\right]}^{\text{T}}$ ， $\beta ={\left({\beta }_{p+1},{\beta }_{p+2},\cdots ,{\beta }_m\right)}^{\text{T}}$ 为常参数向量， $x_{10}=x_{20}=\cdots =x_{n0}=1$ ， $X=\left[\begin{array}{cccc}{x}_{1\left(p+1\right)}& x_{1\left(p+2\right)}& \cdots & x_{1m}\\ x_{2\left(p+1\right)}& x_{21}& \cdots & x_{2m}\\ ⋮& ⋮& & ⋮\\ x_{n\left(p+1\right)}& x_{n\left(p+2\right)}& \cdots & x_{nm}\end{array}\right]$ 为 $n\times \left(m-p\right)$ 维矩阵。W为空间权重矩阵， $W=\left[\begin{array}{cccc}0& b_{12}& \cdots & b_{1n}\\ b_{21}& 0& \cdots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ b_{n1}& b_{n2}& \cdots & 0\end{array}\right]$ ， $b_{ij}=\frac{1}{d_{ij}^2}$ ， $d_{ij}$ 为位置i和位置j之间的距离，且 $d_{ij}=\sqrt{{\left(u_i-u_j\right)}^2+{\left(v_i-v_j\right)}^2}$ ，其中 $u\left(i\right),v\left(i\right)$ 分别为位置i的经度和纬度。WY为空间滞后因子， $\rho \in \left[0,1\right]$ 为空间滞后项系数，取值越接近1，说明相邻地区的因变量取值越相似。 $\epsilon ={\left[{\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_n\right]}^{\text{T}}$ 为服从正态分布的n维随机误差向量且 $Cov\left({\epsilon }_i,{\epsilon }_j\right)=0,i\ne j$ ，I为 $n\times n$ 的单位矩阵。

若M中的系数项变为常系数，则模型变为空间滞后模型(SLM)。

3. 参数估计法

3.1. 两步估计法

根据模型特点，可知模型由空间滞后部分 $\rho WY$ ，变系数部分M以及常系数部分 $X\beta$ 组成，因此可以优化两步估计法 [7] [8] 来估计模型参数 $\rho ,\beta ,{\sigma }^2$ ，具体步骤如下：

第一步，假设空间滞后部分 $\rho WY$ 中的 $\rho$ 一常系数部分 $X\beta$ 中的 $\beta$ 已知，则此模型 $Y=\rho WY+X\beta +M+\epsilon$ 可以化为

$\left(I_n-\rho W\right)Y-X\beta =M+\epsilon .$ (2)

令 $A=I_n-\rho W$ ， $Y^{*}=\left(I_n-\rho W\right)Y-X\beta =AY-X\beta$ ，则式(2)变为 $Y^{*}=M+\epsilon$ 。对于该模型，利用空间局部加权最小二乘法可知：

$\hat{M}=S{Y}^{*}$ (3)

其中令 $S_i=X_i^{\text{T}}{\left(X^{\text{T}}{W}_{i}X\right)}^{-1}{X}^{\text{T}}{W}_i$ ，将 $S_i$ 称为i点的帽子向量 [9] ， $X_i^{\text{T}}$ 为矩阵X的第i行。 $S=\left[\begin{array}{c}{S}_1\\ S_2\\ ⋮\\ S_n\end{array}\right]=\left[\begin{array}{c}{X}_1^{\text{T}}{\left(X^{\text{T}}{W}_{1}X\right)}^{-1}{X}^{\text{T}}{W}_1\\ X_2^{\text{T}}{\left(X^{\text{T}}{W}_{2}X\right)}^{-1}{X}^{\text{T}}{W}_2\\ ⋮\\ X_n^{\text{T}}{\left(X^{\text{T}}{W}_{n}X\right)}^{-1}{X}^{\text{T}}{W}_n\end{array}\right]$ 为帽子矩阵。在帽子矩阵S中， $W_i=diag\left(w_{i1},w_{i2},\cdots ,w_{in}\right)$ ，此时 $w_{ij}$ 为高斯空间权函数，带宽b由 $CV={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-{\hat{y}}_{\ne i}\left(b\right)\right]}^2}$ [10] 确定。

将 $\hat{M}=S{Y}^{*}$ 代入 $Y^{*}=M+\epsilon$ 中可得

$Y^{*}=S{Y}^{*}+\epsilon$ (4)

第二步，利用最大似然估计的方法对 $Y^{*}=S{Y}^{*}+\epsilon$ 进行求解，由于 $\epsilon =\left(I_n-S\right)Y^{*}$ ，则由变换定理可得Y的似然函数 $p\left(Y\right)=p\left(\epsilon \right)\left|\frac{\partial \epsilon }{\partial Y}\right|=p\left(\epsilon \right)\left|\frac{\partial \left[\left(I_n-S\right)\left(AY-X\beta \right)\right]}{\partial Y}\right|=p\left(\epsilon \right)\left|I_n-S\right|\left|A\right|$ ，可得Y的对数似然函数为：

$\begin{array}{c}\ln L=-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln {\sigma }^2+\ln \left|I_n-S\right|+\ln \left|A\right|-\frac{1}{2{\sigma }^2}{\left(Y^{*}-S{Y}^{*}\right)}^{\text{T}}\left(Y^{*}-S{Y}^{*}\right)\\ =-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln {\sigma }^2+\ln \left|I_n-S\right|+\ln \left|I_n-\rho W\right|-\frac{1}{2{\sigma }^2}{\left[\left(I_n-S\right)Y^{*}\right]}^{\text{T}}\left[\left(I_n-S\right)Y^{*}\right].\end{array}$

求解 $\frac{\partial \ln L}{\partial \beta }=0$ ， $\frac{\partial \ln L}{\partial {\sigma }^2}=0$ ，可得：

$\begin{array}{l}\hat{\beta }\left(\rho \right)={\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY,\\ {\hat{\sigma }}^2\left(\rho \right)=\frac{1}{n}{\epsilon }^{\text{T}}\epsilon =\frac{1}{n}{\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right]}^{\text{T}}Q\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right].\end{array}$

其中 $Q={\left(I_n-S\right)}^{\text{T}}\left(I_n-S\right)$ ，将 $\hat{\beta }\left(\rho \right),{\hat{\sigma }}^2\left(\rho \right)$ 代入上述对数似然函数 $\ln L$ ，可得：

$\ln L=C-\frac{n}{2}\ln \left({\left|A\right|}^{-\frac{2}{n}}{\epsilon }^{\text{T}}\epsilon \right)$

其中C为常数。

假设W的特征值为 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ ，则W可对角化，因此存在不可逆矩阵P，使得W变为对角阵 $\Lambda$ ，即 $\Lambda =P^{-1}WP$ ，可得：

$P^{-1}\left(\lambda I_n-W\right)P=P^{-1}\lambda P-P^{-1}WP=\lambda I_n-\Lambda ,$

$\left|\lambda I_n-W\right|={\displaystyle \underset{i=1}{\overset{n}{\prod }}\left(\lambda -{\lambda }_i\right)}$

可得 [11]

$\left|A\right|=\left|I_n-\rho W\right|={\displaystyle \underset{i=1}{\overset{n}{\prod }}\left(1-\rho \cdot {\lambda }_i\right)}$

因此求下式极小值，可得 $\rho$ 的估计值 $\hat{\rho }$

${\left|A\right|}^{-\frac{2}{n}}{\epsilon }^{\text{T}}\epsilon ={\displaystyle \underset{i=1}{\overset{n}{\prod }}{\left(1-\rho \cdot {\lambda }_i\right)}^{-\frac{2}{n}}}{\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right]}^{\text{T}}\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right]$ (5)

将 $\rho$ 的估计值 $\hat{\rho }$ 代入 $\hat{\beta }\left(\rho \right),{\hat{\sigma }}^2\left(\rho \right)$ ，即可得到 $\beta \left(\rho \right),{\sigma }^2\left(\rho \right)$ 的最终估计值 $\hat{\beta }\left(\hat{\rho }\right),{\hat{\sigma }}^2\left(\hat{\rho }\right)$ 。由此可得，模型中变系数部分 $\hat{M}=S{Y}^{*}=S\left(AY-X\beta \right)$ ，此时模型中 $\rho ,W,\beta ,M,{\sigma }^2$ 均已知，因此可以使用SLM-MGWR模型进行预测。

模型 $Y=\rho WY+X\beta +M+\epsilon$ 的参数 $\hat{\beta }\left(\hat{\rho }\right),{\hat{\sigma }}^2\left(\hat{\rho }\right)$ 求得分别为

$\begin{array}{l}\hat{\beta }\left(\hat{\rho }\right)={\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY,\\ {\hat{\sigma }}^2\left(\hat{\rho }\right)=\frac{1}{n}{\epsilon }^{\text{T}}\epsilon =\frac{1}{n}{\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right]}^{\text{T}}Q\left[AY-X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}QAY\right].\end{array}$

其中 $Q={\left(I_n-S\right)}^{\text{T}}\left(I_n-S\right)$ ， $A=I_n-\hat{\rho }W$ 。

3.2. 莫兰指数优化法

3.2.1. $\rho$ 和莫兰指数的关系

莫兰指数是用来度量空间相关性的重要指标，在SL-MGWR模型中，由于空间滞后因子的存在，所以我们考虑将空间滞后项转化为莫兰指数，找到莫兰指数和空间滞后项的关系，由文献 [12] 可知莫兰指数具有如下形式：

$\begin{array}{l}{M}_c=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}\left(y_i-\bar{y}\right)\left(y_j-\bar{y}\right)}}{\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-\bar{y}\right)}^2}}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}}=\frac{n}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-\bar{y}\right)}^2}}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}\left(y_i-\bar{y}\right)\left(y_j-\bar{y}\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}}\\ =\frac{n}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(y_i-\bar{y}\right)}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}\left(y_j-\bar{y}\right)}{\left(n-1\right)s^2}=\frac{n}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{z}_i{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}{z}_j}}{n-1}.\end{array}$ (6)

其中 $W=\left[\begin{array}{cccc}0& w_{12}& \cdots & w_{1n}\\ w_{21}& 0& \cdots & w_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ w_{n1}& w_{n2}& \cdots & 0\end{array}\right]$ 为空间权重矩阵， $Y^{\text{T}}=\left[y_1,y_2,\cdots ,y_n\right]$ 为随机变量。上式(6)右侧强调了莫兰

指数MC对双变量的回归，由文献 [12] 可知MC为回归斜率的系数。对随机变量Y值进行标准差标准化处理之后得到响应变量Z，可知WZ关于Z的回归系数为

${\left(Z^{\text{T}}Z\right)}^{-1}{Z}^{\text{T}}WZ=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Z}_i}\left({\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}{Z}_j\right)}{n-1}.$

则向量W1关于 $1_{n\times 1}={\left(1,1,\cdots ,1\right)}^{\text{T}}$ 的回归系数为

${\left(1^{\text{T}}1\right)}^{-1}{1}^{\text{T}}W1=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}}{n}.$

如果空间权重矩阵为行标准化矩阵，则有 ${\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}=1$ ， ${\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}}=n$ ，MC变为 $M_c=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{z}_i{\displaystyle \underset{j=1}{\overset{n}{\sum }}{w}_{ij}}{z}_j}}{n-1}$ ，此时 $M_c={\left(Z^{\text{T}}Z\right)}^{-1}{Z}^{\text{T}}WZ$ 。由此可知：

$WZ=M_c\cdot Z+\epsilon .$ (7)

其中 $M_c$ 为莫兰指数， $\epsilon ~N\left(O,{\sigma }^2{I}_n\right)$ ，且误差之间相互独立。

3.2.2. 模型参数估计

根据模型特点，可知模型由空间滞后部分 $\rho WY$ ，变系数部分M，以及常系数部分 $X\beta$ 组成，因此可以借鉴两步估计法 [7] [8] 以及文献 [12] 中 $\rho$ 和莫兰指数的关系，来估计模型未知参数 $\rho ,\beta ,{\sigma }^2$ ，具体步骤如下：

第一步，假设SL-MGWR中莫兰指数 $M_c$ 、滞后项系数 $\rho$ 和常系数部分 $\beta$ 已知，同时若空间权重矩阵W为行标准化矩阵，则由式(7)可知：

$WY=M_c\cdot Y+{\epsilon }_1,{\epsilon }_1~N\left[O,{\sigma }_1^2{I}_n\right]$

在原模型(1)的基础上，可知：

$Y=\rho WY+M+X\beta +{\epsilon }_2,{\epsilon }_2~N\left[O,{\sigma }_2^2{I}_n\right]$

则此模型可以变为

$\begin{array}{c}Y=\rho WY+M+X\beta +{\epsilon }_2=\rho \left(M_c\cdot Y+{\epsilon }_1\right)+M+X\beta +{\epsilon }_2\\ =\rho M_c\cdot Y+M+X\beta +\left(\rho {\epsilon }_1+{\epsilon }_2\right)=\rho M_c\cdot Y+M+X\beta +\epsilon \end{array}$

其中 $\epsilon =\rho {\epsilon }_1+{\epsilon }_2$ ， $\epsilon ~N\left(O,\left({\rho }^2{\sigma }_1^2+{\sigma }_2^2\right)I_n\right)$ 。

即：

$\left(1-\rho M_c\right)Y-X\beta =M+\epsilon$ (8)

令 $Y^{*}=\left(1-\rho \cdot M_c\right)Y-X\beta$ ，则式(2)变为 $Y^{*}=M+\epsilon$ 。对于该模型，利用空间局部加权最小二乘法可知：

$\hat{M}=S{Y}^{*}$ (9)

其中令 $S_i=X_i^{\text{T}}{\left(X^{\text{T}}{W}_{i}X\right)}^{-1}{X}^{\text{T}}{W}_i$ ，将称 $S_i$ 为i点的帽子向量 [9] ， $X_i^{\text{T}}$ 为矩阵X的第i行。

$S=\left[\begin{array}{c}{S}_1\\ S_2\\ ⋮\\ S_n\end{array}\right]=\left[\begin{array}{c}{X}_1^{\text{T}}{\left(X^{\text{T}}{W}_{1}X\right)}^{-1}{X}^{\text{T}}{W}_1\\ X_2^{\text{T}}{\left(X^{\text{T}}{W}_{2}X\right)}^{-1}{X}^{\text{T}}{W}_2\\ ⋮\\ X_n^{\text{T}}{\left(X^{\text{T}}{W}_{n}X\right)}^{-1}{X}^{\text{T}}{W}_n\end{array}\right]$

在帽子矩阵S中， $W_i=diag\left(w_{i1},w_{i2},\cdots ,w_{in}\right)$ ，此时 $w_{ij}$ 为高斯空间权函数，带宽b由

$CV={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_i-{\hat{y}}_{\ne i}\left(b\right)\right]}^2}$ [10] 确定。

将(9)式中M的估计值代入 $Y^{*}=M+\epsilon$ 中可得

$Y^{*}=S{Y}^{*}+\epsilon$ (10)

第二步，利用最大似然估计的方法对式(10)进行求解，由于 $\epsilon =\left(I_n-S\right)Y^{*}$ ，则由变换定理可得Y的似然函数 $p\left(Y\right)=p\left(\epsilon \right)\left|\frac{\partial \epsilon }{\partial Y}\right|=p\left(\epsilon \right)\left|\frac{\partial \left[\left(I_n-S\right)\left(\left(1-\rho M_c\right)Y-X\beta \right)\right]}{\partial Y}\right|=p\left(\epsilon \right)\left|I_n-S\right|\left|1-\rho M_c\right|$ ，可得Y的对数似然函数为：

$\begin{array}{c}\ln L=-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln {\sigma }^2+\ln \left|I_n-S\right|+\ln \left(1-\rho M_c\right)-\frac{1}{2{\sigma }^2}{\left(Y^{*}-S{Y}^{*}\right)}^{\text{T}}\left(Y^{*}-S{Y}^{*}\right)\\ =-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln {\sigma }^2+\ln \left|I_n-S\right|+\ln \left(1-\rho M_c\right)-\frac{1}{2{\sigma }^2}{\left[\left(I_n-S\right)Y^{*}\right]}^{\text{T}}\left[\left(I_n-S\right)Y^{*}\right]\end{array}$

求解 $\frac{\partial \ln L}{\partial \beta }=0$ ， $\frac{\partial \ln L}{\partial {\sigma }^2}=0$ ，可得：

$\begin{array}{l}\hat{\beta }\left(\rho \right)={\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}Q\left(1-\rho \cdot M_c\right)Y,\\ {\hat{\sigma }}^2\left(\rho \right)=\frac{1}{n}{\epsilon }^{\text{T}}\epsilon =\frac{1}{n}{\left[\left(1-\rho M_c\right)Y-H\right]}^{\text{T}}Q\left[\left(1-\rho M_c\right)Y-H\right]\end{array}$

其中 $Q={\left(I_n-S\right)}^{\text{T}}\left(I_n-S\right)$ ， $H=X\hat{\beta }=X{\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}Q\left(1-\rho \cdot M_c\right)Y$ 。

将 $\hat{\beta }\left(\rho \right),{\hat{\sigma }}^2\left(\rho \right)$ 代入 $\ln L$ 可得关于 $\rho$ 的集中对数似然函数，

$\begin{array}{c}\ln L=-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln \frac{1}{n}-\frac{n}{2}+\ln \left(1-\rho M_c\right)+\ln \left|I_n-S\right|\\ -\frac{n}{2}\ln {\left[\left(1-\rho M_c\right)Y-H\right]}^{\text{T}}Q\left[\left(1-\rho M_c\right)Y-H\right]\end{array}$ (11)

上述式子为 $\rho$ 的非线性函数，因此我们采用优化算法0.618法求其极大值，得到 $\rho$ 的估计值 $\hat{\rho }$ ，再将 $\hat{\rho }$ 代入到 $\hat{\beta }\left(\rho \right),{\hat{\sigma }}^2\left(\rho \right)$ 即可得到其的最终估计值 $\hat{\beta }\left(\hat{\rho }\right),{\hat{\sigma }}^2\left(\hat{\rho }\right)$ 。

模型 $Y=\rho M_c\cdot Y+M+X\beta +\epsilon$ 的参数 $\hat{\beta }\left(\hat{\rho }\right),{\hat{\sigma }}^2\left(\hat{\rho }\right)$ 求得分别为

$\begin{array}{l}\hat{\beta }\left(\hat{\rho }\right)={\left(X^{\text{T}}QX\right)}^{-1}{X}^{\text{T}}Q\left(1-\rho \cdot M_c\right)Y,\\ {\hat{\sigma }}^2\left(\hat{\rho }\right)=\frac{1}{n}{\epsilon }^{\text{T}}\epsilon =\frac{1}{n}{\left[\left(1-\rho M_c\right)Y-H\right]}^{\text{T}}Q\left[\left(1-\rho M_c\right)Y-H\right]\end{array}$

其中 $Q={\left(I_n-S\right)}^{\text{T}}\left(I_n-S\right)$ 。

4. 实验研究

4.1. 数据分层实验

由于空间数据具有异质性，即在描述空间关系的参数在研究区域的不同地方是不同的，但在局部的变化是一致的。因此，空间异质性的存在导致在空间数据分析过程中，需要强调对局部的识别和分析，否则很难保证结果的可靠性，甚至会得到错误的结论。本文为证明分层的有效性，在GWR模型和十折交叉验证的基础上，通过添加数据分层来比较拟合效果。实验通过对空间数据进行聚类，得到分层的数据，在每一层中使用十折交叉验证的方法，选择出训练集和验证集，使用加权最小二乘法对GWR参数进行估计，并使用差异函数均值 $E\left(r\right)$ 给出最后的拟合效果。

本文所有实验设置的对比指标均为分层数(group_index)，验证集数目(validation_num)，训练集数目(train_num)，总数(ALL_num)，第k层差异函数值 $r_k$ ，第k层权重 $w_k$ 以及差异函数均值 $E\left(r\right)$ ，实验结果如表1和表2所示。

Table 1. Ten fold cross validation GWR

表1. 十折交叉验证GWR

Table 2. Stratified tenfold cross-validation GWR

表2. 分层十折交叉验证GWR

对比两表可得对数据进行分层处理，能够得到较好的预测效果，数据分区后相比分区前提高了0.33，由此发现数据分区可以提高模型的拟合效果，证明了数据分层的有效性。因此本节之后的所有实验均采用分层十折交叉验证处理数据。

4.2. 空间滞后–混合地理加权回归模型对比实验

为同时考虑数据的异质性及自相关性，降低因为地理位置远近而造成数据差异的影响及空间相邻单元相互作用的影响，我们提出SL-MGWR模型。本节使用两步估计法将SL-MGWR模型的拟合效果与GWR的拟合效果进行比较，来证明模型SL-MGWR的有效性。首先使用聚类分析将数据总共分为8层，每层使用十折交叉验证并选择合适带宽。根据十折交叉验证的方法可知，每折对应一个相应的滞后因子ρ。确定每层带宽和训练集、验证集划分后，使用两步估计法可得到模型参数对数据进行预测和拟合。如果某分层内数据量过小，我们对本层整体采用SL-MGWR模型进行处理。

在本文中，我们使用优化算法0.618法可以求得式(5)的极小值，在极小值处确定对应的滞后因子ρ，每折滞后因子及每层的差异函数均值如表3所示。

Table 3. Two-step estimation method of SL-MGWR model

表3. SL-MGWR模型两步估计法

由表3可知，在第2层和第8层内数据不存在自相关性，每折滞后因子ρ均取0，因此第2层和第8层数据均使用混合地理加权模型进行处理。第3层只有第6折处的数据存在自相关性，其他折不存在自相关性。SL-MGWR模型两步估计法进一步说明模型进行十折交叉验证的有效性，能够较全面地处理数据。

由表2和表3可得 $E\left(r\right)$ 的大小，可知SL-MGWR模型的预测误差比GWR模型的预测误差低0.49，SL-MGWR模型的拟合效果优于GWR的拟合效果，证明SL-MGWR模型的有效性。

4.3. 莫兰指数估计法参数估计对比实验

本节分别使用两步估计法及莫兰指数优化法对SL-MGWR模型进行预测和拟合，通过对比模型的拟合效果，证明莫兰指数优化法的优越性。在莫兰指数优化法中，代入权重矩阵和莫兰指数I的关系式对模型(1)进行处理，使得最后的参数估计更加准确，得到更为精确的结果，此时同时输出每折的莫兰指数I以及滞后因子ρ，具体数值如表4所示。

比较表3、表4可知，同一组数据聚类分析结果一致，且每层数据的训练集和验证集数目相同。在表4中，也存在滞后因子为0的情况，说明对应数据不存在空间滞后效果。对比表3、表4中 $E\left(r\right)$ 的大小，莫兰指数优化法比两步估计法得到预测值误差低0.6，因此莫兰指数优化法拟合效果更好，证明莫兰指数优化法的优越性。

Table 4. Moran index optimization method of SL-MGWR model

表4. SL-MGWR模型莫兰指数优化法

5. 结论

本文在地理加权回归模型和空间滞后模型基础上做出改进，提出空间滞后–混合地理加权回归模型，空间滞后–混合地理加权回归模型优于原有模型，具有有效性。在空间数据处理过程中，使用k均值聚类方法对空间数据进行分区处理，可以减少空间异质性对数据的影响，降低预测数据的误差。在模型参数求解过程中，莫兰指数优化法给出空间滞后项与莫兰指数的数量关系，进一步提高了拟合效果。通过实验证明分区和莫兰指数优化法对提升模型拟合效果具有积极意义。通过在真实数据集上的实验证明，本文所提方法可以有效提高模型拟合效果。