空间计量经济学主要研究的是经济学中的空间效应问题，包括空间相关性和空间异质性。空间相关性表现为空间单元变量在一定程度上受到其他位置的空间单元影响，并且这种影响往往与距离存在一定相关性，而这种距离既可以是地理空间意义的距离，也可以是经济和社会的距离。空间异质性主要指地理空间上的区域缺乏均质性，从不同的空间观测点来观测会产生不同的观测结果。空间异质性反映了空间单元之间经济行为关系普遍存在的非均质性。应该说，空间相关性和空间异质性是空间计量经济学模型构建的理论基石。

在实际经济社会发展中，由于空间相关性的客观存在，空间自相关成为空间计量模型建立的依据。变量的空间自相关性表现为变量的数值不仅与所在空间位置点有关，还会受到其他空间点上变量取值的影响。自相关性的判断主要用于进行探索性空间数据分析，包括全局相关性和局部相关性，涉及常用的2种指数主要是Moran’s I指数和Geary G指数^[1-2]。Fisher^[3]、Kelejian和Prucha^[4-7]、Lee^[8-12]对考虑空间相关性和空间异质性的截面线性回归模型进行研究，认为当空间系数在各空间点上不再恒定，此时采取变系数形式更合适。空间变系数的提出反映了空间局域性特征，但是增加了模型估计中待求解未知系数的个数，一定程度上增加了估计求解的难度。可以通过引入空间权重矩阵采用加权最小二乘法加以估计，而最合适的空间权重设置则是通过交叉验证所有空间点因变量的观测值与估计值最小平方和来实现。Brunsdon等^[13-16]和Fotheringham等^[17]正是基于上述思路提出地理加权回归(geographically weighted regression, GWR)模型来揭示空间非均质性。GWR模型强调空间变量系数有别于传统的回归模型，它是变化的，并且与所在空间点和其他空间点的距离相关。这样，空间权重矩阵的引入也使得传统假设相互独立的空间样本点变得相互依赖起来。

近年来，GWR模型在处理非均质空间数据得到了广泛应用。从空间变量来看，空间相关性是变量在空间分布特征最直接的体现，空间异质性则是空间结构的不稳定性表现。对于GWR而言，自变量并非全部是不稳定的空间结构，因此Brunsdon等^[16]提出混合地理加权回归(mixed geographically weighted regression, MGWR)，考虑了仅有部分自变量是稳定的空间结构。除自变量的空间异质性外，因变量本身还存在着空间滞后效应，至少从模型等式关系来看，自变量的空间效应也会影响到因变量。乔宁宁^[18]提出常系数空间滞后的混合地理加权回归(mixed geographically weighted regression with constant-coefficient spatial lag regression, MGWR-CSLR)模型。事实上，空间异质性结构不仅影响自变量，也影响着因变量，Brunsdon等^[15]、Paze等^[19]对于传统GWR进一步研究了因变量的空间滞后项回归，将滞后项的常系数替代为变系数，并且Brunsdon等^[15]认为各区域点因变量空间滞后的系数即使没有显著差异，也并不意味着常系数适用于每个地区。综上，从已有的GWR模型理论成果总结来看，不论是因变量还是自变量，所对应的变系数应该看作是常系数的一般性推广。常系数体现为变系数的一种特殊情形，譬如当表征变量的空间点由显著的空间集聚(扩散)分布变为随机分布，表征变量所对应每个空间点的变系数也就退化为对应所有点均值的常系数。由于模型变量需要考虑的空间特征越多，系数估计及检验也就越复杂，至今关于GWR考虑各种空间特征的复杂形式的理论与实践研究均有待于解决。

为此，在前人研究基础上，本文试图对GWR模型做一个更具广泛意义的推广形式，将可能同时存在的因变量空间滞后性与自变量空间异质性2种效应融合并进行拓展研究，主要工作分为3方面：1)将变系数空间滞后回归与MGWR进行融合，放到一个全模型下进行考虑，因变量的空间滞后性与自变量空间异质性随着空间点变化能同时兼顾到并能充分展现，与乔宁宁^[18]提出的MGWR-CSLR模型相对应，将这一类新的模型称为变系数空间滞后的混合地理加权回归(mixed geographically weighted regression with varying-coefficient spatial Lag regression, MGWR-VSLR)模型。2)采用Xiong^[20]提出的重构参数化方法对MGWR-VSLR模型进行系数估计，实现同步异带宽估计并给出模型与系数估计值的显著性检验，以及变系数选取的判别检验。3)基于蒙特卡罗模拟和实际应用检验了不同GWR形式的模型在处理空间相关性和空间异质性时所反映出的优良性质。基于本文的分析，有助于处理空间问题时选取适宜的GWR模型。

1 MGWR-VSLR模型及其系数估计 1.1 MGWR-VSLR模型

对于横截面空间数据来说，采用普通最小二乘回归(ordinary least square regression, OLSR)模型通常假定空间样本点相互独立，空间结构是稳定的，因此回归系数实质是样本均值回归的结果，该系数是不变的且适用于所有空间样本点，自变量对因变量的影响是全局性的。对所有空间点i而言，OLSR模型自变量所对应的系数是不变的：

$ y_i=\beta_0+\sum\limits_{j=1}^p \alpha_j x_{i j}+\varepsilon_i, \quad i=1, 2, \cdots, n . $

(1)

式中：y_i和x_i1, x_i2, …, x_ip都是对应空间点可观测的一般变量，其中β₀, α₁, α₂, …, α_p为回归系数，ε_i为相应的随机误差项，ε_i~N(0, σ²)。对于任意第j(1≤j≤p) 个自变量X_j中所有n个独立空间点观测值x_1j, x_2j, …, x_nj，所对应回归系数α_j是一样的，此时自变量也被看成是全局变量，对应系数即为常系数。当考虑到空间结构并不稳定，存在空间异质性时，则每个空间点的回归系数都是变化的，此时系数即为变系数，而自变量也被看成是局域变量。在GWR中每个空间样本点对应的回归系数是变化的，自变量对因变量的影响是局域性的，Brundon等^[15]引入空间点距离的高斯权重后，通过采用局部加权的最小二乘估计来求解变系数。

$ \begin{gathered} y_i=\beta_0\left(u_i, v_i\right)+\sum\limits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, \\ i=1, 2, \cdots, n . \end{gathered} $

(2)

式(2)即为GWR模型，第k个自变量每个空间点i的系数均与所在空间位置的坐标(u_i, v_i)相关，不同空间点i所对应的回归系数也有所差异。随后，当考虑到仅有部分自变量在空间点上存在空间异质性，而其余变量空间异质性并不明显时，也即全局变量与局域变量同时存在于一个地理问题中，Brunsdon等^[16]提出MGWR模型表示如下

$ \begin{gathered} y_i=\sum\nolimits_{j=1}^p \alpha_j x_{i j}+\beta_0\left(u_i, v_i\right)+ \\ \sum\nolimits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, \\ i=1,2, \cdots, n . \end{gathered} $

(3)

式中: 有p个自变量对应的系数α_j(j=1, 2, …, p)被认为是常数，不会随空间位置变化而变化，其余剩下的自变量对应系数β_k以及截距项β₀则随着空间位置变化而变化，可以采用两步估计法^[21]求解。

尽管上述MGWR较为完整地反映了全局变量与局域变量空间结构的稳定性和不稳定性，却没有考虑到因变量存在的空间滞后性。乔宁宁^[18]在MGWR模型的基础上引入因变量的空间滞后性，建立了一种具有因变量空间滞后项的MGWR模型，由于假设因变量空间滞后项对不同空间点因变量的影响是一样的，也可以看成是MGWR-CSLR模型，如

$ \begin{gathered} y_i=\rho \sum\nolimits_{j=1}^n w_{i j} y_j+\sum\nolimits_{j=1}^p \alpha_j x_{i j}+\beta_0\left(u_i, v_i\right)+ \\ \sum\nolimits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, \quad i=1, \cdots, n . \end{gathered} $

(4)

$ \begin{gathered} \text { 或 } y_i=\rho \sum\nolimits_{j=1}^n w_{i j} y_j+\alpha_0+\sum\nolimits_{j=1}^p \alpha_j x_{i j}+ \\ \sum\nolimits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, \quad i=1, 2, \cdots, n . \end{gathered} $

(5)

在式(4)和式(5)中增加了因变量的空间滞后性，w_ij为预先给定的n×n阶空间邻接矩阵元素，模型假设因变量空间滞后性对每个空间点因变量的影响均是相同的，因此，这里因变量的空间滞后项对应的系数ρ为待估计的常系数。

其中，式(4)的截距项β₀(u_i, v_i)为一个与空间位置(u_i, v_i)有关的变参数，式(5)的截距项α₀为一个与空间位置无关的常参数。乔宁宁^[18]主要针对式(5)的求解过程做了详细介绍。在式(4)和式(5)中由于因变量的空间滞后项为常数ρ，表示因变量y_i所有空间点的滞后项对因变量影响均是一样的。Brunsdon^[15]曾提出变系数的空间滞后回归(varying-coefficient spatial lag regression, VSLR)模型，认为即使ρ(u_i, v_i)在各个空间点差异不显著也并不意味着各个空间点适用于固定的常数ρ。从空间滞后的MGWR模型的形式来看，自变量对应的变系数同样可能会影响因变量，促使因变量表现出局域性特征，将空间滞后项对应的常系数设为变系数，更能反映因变量局部情况。因此，与式(4)和式(5)的滞后项所对应的常系数ρ类似，当滞后项前系数为变系数ρ(u_i, v_i)，分别得到：

$ \begin{aligned} & y_i= \rho\left(u_i, v_i\right) \sum\nolimits_{j=1}^n w_{i j} y_j+\sum\nolimits_{j=1}^p \alpha_j x_{i j}+\beta_0\left(u_i, v_i\right)+ \\ & \sum\nolimits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, i=1, 2, \cdots, n . \quad(6) \end{aligned} $

(6)

$ \begin{gathered} \text { 或 } y_i=\rho\left(u_i, v_i\right) \sum\nolimits_{j=1}^n w_{i j} y_j+\alpha_0+\sum\nolimits_{j=1}^p \alpha_j x_{i j}+ \\ \sum\nolimits_{k=1}^q \beta_k\left(u_i, v_i\right) z_{i k}+\varepsilon_i, i=1, 2, \cdots, n . \quad \text { (7) } \end{gathered} $

(7)

式(6)和式(7)的差异主要体现在截距项的区别，在一般情况下，式(6)可以看作是式(7)的更一般形式。因此为使本文研究的模型更具一般性，后面的讨论均以式(6)为例。MGWR-VSLR模型的矩阵形式可写为

$ \boldsymbol{Y}=\boldsymbol{\rho}_n \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\alpha}+\boldsymbol{M}+\boldsymbol{\varepsilon}. $

(8)

其中：

$ \begin{aligned} & \boldsymbol{Y}=\left[\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right], \quad \boldsymbol{X}=\left[\boldsymbol{X}_{1}, \boldsymbol{X}_{2}, \cdots, \boldsymbol{X}_{p}\right] ,\\ & \boldsymbol{\rho}_{n}=\left[\begin{array}{cccc} \rho\left(u_{1}, v_{1}\right) & 0 & \cdots & 0 \\ 0 & \rho\left(u_{2}, v_{2}\right) & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & \rho\left(u_{n}, v_{n}\right) \end{array}\right], \\ & \boldsymbol{X}_{i}=\left[\begin{array}{c} x_{1 i} \\ x_{2 i} \\ \vdots \\ x_{n i} \end{array}\right], \boldsymbol{\alpha}=\left[\begin{array}{c} \alpha_{1} \\ \alpha_{2} \\ \vdots \\ \alpha_{p} \end{array}\right], \\ & \boldsymbol{\beta}=\left[\begin{array}{cccc} \beta_{0}\left(u_{1}, v_{1}\right) & \beta_{0}\left(u_{2}, v_{2}\right) & \cdots & \beta_{0}\left(u_{n}, v_{n}\right) \\ \beta_{1}\left(u_{1}, v_{1}\right) & \beta_{1}\left(u_{2}, v_{2}\right) & \cdots & \beta_{1}\left(u_{n}, v_{n}\right) \\ \vdots & \vdots & \vdots & \vdots \\ \beta_{q}\left(u_{1}, v_{1}\right) & \beta_{q}\left(u_{2}, v_{2}\right) & \cdots & \beta_{q}\left(u_{n}, v_{n}\right) \end{array}\right]= \end{aligned}\\ \begin{gathered} {\left[\begin{array}{c} \boldsymbol{\beta}_{0}(u, v)^{\mathrm{T}} \\ \boldsymbol{\beta}_{1}(u, v)^{\mathrm{T}} \\ \vdots \\ \boldsymbol{\beta}_{q}(u, v)^{\mathrm{T}} \end{array}\right], \boldsymbol{Z}=\left[\boldsymbol{1}_{n}, \boldsymbol{Z}_{1}, \cdots, \boldsymbol{Z}_{q}\right], } \\ \boldsymbol{Z}_{i}=\left[\begin{array}{c} z_{1 i} \\ z_{2 i} \\ \vdots \\ z_{n i} \end{array}\right], \boldsymbol{\varepsilon}=\left[\begin{array}{c} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{n} \end{array}\right] \\ \boldsymbol{W}=\left[\begin{array}{cccc} w_{11} & w_{12} & \cdots & w_{1 n} \\ w_{21} & w_{22} & \cdots & w_{2 n} \\ \vdots & \vdots & \vdots & \vdots \\ w_{n 1} & w_{n 2} & \cdots & w_{n n} \end{array}\right] . \end{gathered} $

设定{d₁, …, d_n}为Zβ的对角线元素，则 M=diag{d₁, …, d_n}1_n，ρ_nWY表示变系数空间滞后项。通过式(8)，可以看出把MGWR-VSLR模型作为全模型时: 1)当不同变量对应的系数不再显著，全模型分别退化为：Y=ρ_n WY+ε为变系数空间自回归(varying-coefficient spatial auto-regression, VSAR)模型；Y=ρ_nWY+Xα+ε为变系数空间滞后回归VSLR模型；Y=ρ_nWY+M+ε为变系数空间滞后的地理加权回归(geographically weighted regression with varying-coefficient spatial lag regression, GWR-VSLR)模型; 2)而当ρ(u_i, v_i)恒等于同一常数ρ时，可以分别相应地得到传统常系数的空间自回归(spatial auto-regression, SAR)模型Y=ρWY+ε；常系数的空间滞后回归(spatial lag-regression, SLR)模型Y=ρWY+Xα+ε；常系数空间滞后的地理加权回归(geographically weighted regression with constant-coefficient spatial lag regression, GWR-CSLR)模型Y= ρWY+M+ε; 3)而当ρ(u_i, v_i)≡0(i=1, …, n)时，此时 Y=Xα+M+ε为MGWR模型。因此，从上述不同模型的形式总结可以看出，MGWR-VSLR模型具有空间回归模型更一般的形式，研究其变量对应系数的估计与统计量检验，有助于更深刻地理解空间回归模型。

1.2 MGWR-VSLR模型系数估计

在前人的研究中，各类GWR都是利用空间局部加权最小二乘法来估计变系数部分，然而应用这个方法难以实现一步估计含有变系数空间滞后项的模型系数。另一方面，传统的GWR方法在估计系数时对所有系数使用同一个带宽，这对不同空间点对应的变系数施加了同样的空间比例，而在现实情况下，不同的变系数可能会在不同的空间尺度上产生影响。因此，本文采用Xiong^[20]提出的重构参数化方法，将未知变系数关于空间位置的非参数回归转化为插值形式后，整个模型的极大似然估计便可以给出，估计过程可以对不同自变量项对应的系数以及因变量滞后项对应的系数同时分配不同的带宽并进行优化。另外，重构参数化方法将需要估计全部样本空间点的系数值转化为对少量节点对应的系数值进行估计，既减少了计算复杂度，又可以得到变系数关于空间位置的插值函数表达，有助于进一步分析变系数在空间上的变化。

参照Xiong^[20]所给出的重构参数化方法求解过程, 假设函数$f(\boldsymbol{x}), \boldsymbol{x} \in \mathbb{R}^{d}$为待估函数, 取线性插值算子$\boldsymbol{\gamma}^{\mathrm{T}} \boldsymbol{b}(\boldsymbol{x})$, 其中$\boldsymbol{\gamma}=\left(\gamma_{1}, \cdots, \gamma_{m}\right)^{\mathrm{T}} \in \mathbb{R}^{m}$为$m$个实数域节点$\mathcal{A}=\left\{\boldsymbol{a}_{1}, \cdots, \boldsymbol{a}_{m}\right\}$对应的真实值, 即$\boldsymbol{\gamma}^{\mathrm{T}} \boldsymbol{b}\left(\boldsymbol{a}_{i}\right)=\gamma_{i}, i=1, \cdots, m_{\text {。}}$因此$f$可以用$\boldsymbol{\gamma}^{\mathrm{T}} \boldsymbol{b}(\boldsymbol{x})$近似重构, 其中$\boldsymbol{\gamma}=\boldsymbol{f}_{\mathcal{A}}=\left(\boldsymbol{f}\left(\boldsymbol{a}_{1}\right), \cdots\right.$, $\left.f\left(\boldsymbol{a}_{m}\right)\right)^{\mathrm{T}}$, 由此函数$f$的估计为

$ \tilde{f}(\boldsymbol{x})=\tilde{\boldsymbol{\gamma}}^{\mathrm{T}} \boldsymbol{b}(\boldsymbol{x}) $

其中$\tilde{\boldsymbol{\gamma}}=\left(\tilde{\gamma}_{1}, \cdots, \tilde{\gamma}_{m}\right)^{\mathrm{T}}$为下面目标函数的最优解：

$ \min\limits_{\gamma \in \mathbb{R}^{m}} \frac{1}{n} \sum\limits_{i=1}^{n}\left[f\left(\boldsymbol{x}_{i}\right)-\boldsymbol{\gamma}^{\mathrm{T}} \boldsymbol{b}\left(\boldsymbol{x}_{i}\right)\right]^{2} . $

在本文中采用核插值，其中b(x) 的计算方法如下

$ \begin{equation*} \boldsymbol{b}(\boldsymbol{x})=P g(\boldsymbol{x})+Q \boldsymbol{r}_{\mathcal{A}}(\boldsymbol{x}) \end{equation*}. $

(9)

其中：

$ \begin{gathered} \boldsymbol{P}=\boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\left(\boldsymbol{G}_{\mathcal{A}}^{\mathrm{T}} \boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\right)^{-1} ,\\ \boldsymbol{Q}=\left[\boldsymbol{I}_{m}-\boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\left(\boldsymbol{G}_{\mathcal{A}}^{\mathrm{T}} \boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\right)^{-1} \boldsymbol{G}_{\mathcal{A}}^{\mathrm{T}}\right] \boldsymbol{R}_{\mathcal{A}}^{-1}, \\ \boldsymbol{r}_{\mathcal{A}}(\boldsymbol{x})=\left(R\left(\boldsymbol{x}-\boldsymbol{a}_{1}\right), \cdots, R\left(\boldsymbol{x}-\boldsymbol{a}_{m}\right)\right)^{\mathrm{T}} ,\\ \boldsymbol{R}_{\mathcal{A}}=\left(\boldsymbol{R}\left(\boldsymbol{a}_{i}-\boldsymbol{a}_{j}\right)\right)_{i=1, \cdots, m, j=1, \cdots, m}, \\ \boldsymbol{g}(\boldsymbol{x})=\left(1, \boldsymbol{x}^{\mathrm{T}}\right)^{\mathrm{T}}, \boldsymbol{G}_{\mathcal{A}}=\left(\boldsymbol{g}\left(\boldsymbol{a}_{1}\right), \cdots, \boldsymbol{g}\left(\boldsymbol{a}_{m}\right)\right)^{\mathrm{T}}. \end{gathered} $

为与传统GWR方法一致，本文中取高斯核函数$R(\boldsymbol{h})=\exp \left(-\sum_{j=1}^{2} h_{j}^{2} / \theta^{2}\right), \boldsymbol{h}=\left(h_{1}, h_{2}\right)^{\mathrm{T}}, \theta>0$为核参数，这里也称为带宽，表示局部影响效力的空间范围。

运用这一方法，给定n个空间点并从中选取m个节点，变系数向量$\boldsymbol{\beta}_{k}=\left(\beta_{k}\left(u_{1}, v_{1}\right), \cdots\right.$, $\left.\beta_{k}\left(u_{n}, v_{n}\right)\right)^{\mathrm{T}}(k=0, 1, \cdots, q)$和空间滞后系数向量$\boldsymbol{\rho}=\left(\rho\left(u_{1}, v_{1}\right), \cdots, \rho\left(u_{n}, v_{n}\right)\right)^{\mathrm{T}}$可以近似写作

$ \begin{equation*} \tilde{\boldsymbol{\beta}}_{k}\left(u_{i}, v_{i}\right)=\boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{b}\left(u_{i}, v_{i}\right), k=0, 1, \cdots, q, \\ \tilde{\rho}\left(u_{i}, v_{i}\right)=\boldsymbol{\phi}^{\mathrm{T}} \boldsymbol{b}\left(u_{i}, v_{i}\right) \text {. } \end{equation*} $

(10)

式中，γ_k=(γ_k1, …, γ_km)^T，$\boldsymbol{\phi}=\left(\phi_{1}, \cdots, \phi_{m}\right)^{\mathrm{T}}$，并且在式(10)中进一步记$\boldsymbol{A}(\boldsymbol{\phi})=\operatorname{diag}\left\{\tilde{\boldsymbol{\rho}}\left(u_{1}, v_{1}\right)\right.$, $\left.\cdots, \tilde{\rho}\left(u_{n}, v_{n}\right)\right\}, \tilde{\boldsymbol{Z}}=\left(\tilde{\boldsymbol{Z}}_{0}, \tilde{\boldsymbol{Z}}_{1}, \cdots, \tilde{\boldsymbol{Z}}_{q}\right), \tilde{\boldsymbol{Z}}_{0}=\left(\boldsymbol{b}\left(u_{1}\right.\right.$, $\left.\left.v_{1}\right), \cdots, \boldsymbol{b}\left(u_{n}, v_{n}\right)\right)^{\mathrm{T}}, \tilde{\boldsymbol{Z}}_{k}=\left(\boldsymbol{b}\left(u_{1}, v_{1}\right) z_{1 k}, \cdots, \boldsymbol{b}\left(u_{n}\right.\right.$, $\left.\left.v_{n}\right) z_{n k}\right)^{\mathrm{T}}, k=1, \cdots, q$。令 γ=(γ₀^T, γ₁^T, …,γ_q^T)^T，$\boldsymbol{T}(\boldsymbol{\phi})=\boldsymbol{I}_{n}-\boldsymbol{A}(\boldsymbol{\phi}) \boldsymbol{W}, \boldsymbol{B}=(\boldsymbol{X}, \tilde{\boldsymbol{Z}})$以及 η=(α^T, γ^T)^T，则 Bη=Xα+M，式(8)可以简化为

$ \begin{equation*} \boldsymbol{Y}=\boldsymbol{T}(\boldsymbol{\phi})^{-1}(\boldsymbol{B} \boldsymbol{\eta}+\boldsymbol{\varepsilon}) \end{equation*}. $

(11)

由于ε~N(0, σ²I_n)，式(11)的对数似然函数为

$ \begin{gather*} \ln L\left(\boldsymbol{\phi}, \boldsymbol{\eta}, \sigma^{2}\right)=-\frac{n}{2} \ln 2 \mathsf{π} \sigma^{2}+\ln |\boldsymbol{T}(\boldsymbol{\phi})|- \\ \frac{1}{2 \sigma^{2}}\left\{[\boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}]^{\mathrm{T}}[\boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}]\right\} . \end{gather*} $

(12)

假设ϕ已知，由式(12)容易得到η和σ²的极大似然函数估计分别为

$ \left\{\begin{array}{l} \tilde{\boldsymbol{\eta}}(\boldsymbol{\phi})=\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}, \\ \tilde{\boldsymbol{\sigma}}^{2}(\boldsymbol{\phi})=\frac{1}{n}(\boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y})^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}. \end{array}\right. $

(13)

其中S=I_n- B(B^T B)^-1 B^T，则ϕ的集中对数似然函数如下所示

$ \begin{gather*} \ln L(\boldsymbol{\phi})=\frac{n}{2}[\ln (2 \mathsf{π})+1]-\frac{n}{2} \ln \tilde{\boldsymbol{\sigma}}^{2}(\boldsymbol{\phi})+ \\ \ln |\boldsymbol{T}(\boldsymbol{\phi})| . \end{gather*} $

(14)

只要进一步通过最小化式(14)，即$\min\limits_{\phi \in \mathbb{R}^{m}} \ln L(\boldsymbol{\phi})$，即可得到 ϕ的估计值$\widetilde{\boldsymbol{\phi}}=\left(\widetilde{\phi}_{1}, \cdots, \widetilde{\phi}_{m}\right)^{\mathrm{T}}$。在已知$\tilde{\boldsymbol{\phi}}$的基础上，再返回式(13)即可求解出$\tilde{\boldsymbol{\eta}}$以及$\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})$，因此，可得到式(11)相应的估计等式

$ \begin{equation*} \tilde{\boldsymbol{Y}}=\boldsymbol{T}(\tilde{\boldsymbol{\phi}})^{-1} \boldsymbol{B} \tilde{\boldsymbol{\eta}} . \end{equation*} $

(15)

对于本文讨论的模型，将不同的带宽应用在不同系数对应的插值形式中可以更好地适应系数的光滑度差异。因此引入kernel参数向量θ=(θ₁, …, θ_q, θ_a)^T, 式(10)分别对应改写为

$ \begin{gather*} \widetilde{\boldsymbol{\beta}_{k}}\left(u_{i}, v_{i}\right)=\boldsymbol{\gamma}_{k}^{\mathrm{T}} \boldsymbol{b}\left(u_{i}, v_{i}, \theta_{k}\right), k=0, 1, \cdots, q \text { 和 } \\ \tilde{\rho}\left(u_{i}, v_{i}\right)=\boldsymbol{\phi}^{\mathrm{T}} \boldsymbol{b}\left(u_{i}, v_{i}, \theta_{a}\right) . \end{gather*} $

(16)

对带宽θ可以通过变异系数(coefficient of variation, CV)方法来估计，对于给定的γ和ϕ，寻找最小的CV值，如下式所示，其中$\tilde{y}_{-i}(\boldsymbol{\theta})$是固定θ已删掉第i组观测数据后的数据集得到的估计值

$ \begin{gather*} \widetilde{\boldsymbol{\theta}}=\operatorname{argmin}_{\boldsymbol{\theta} \subset \mathbb{R}^{q+1}} \operatorname{CV}(\boldsymbol{\theta})= \\ \operatorname{argmin}_{\boldsymbol{\theta} \subset \mathbb{R}^{q+1}} \frac{1}{n} \sum\nolimits_{i=1}^{n}\left(y_{i}-\tilde{y}_{-i}(\boldsymbol{\theta})\right)^{2}. \end{gather*} $

(17)

2 MGWR-VSLR模型的相关统计检验 2.1 模型整体的显著性检验

对于回归模型整体的显著性检验，主要是考察模型中所有参数存在的合理性，以及对整个模型的解释力度。若所有参数不能通过显著性检验，则模型本身失去存在的意义。因此，模型的显著性检验中存在零假设H₀：τ=0，其中 τ=(ρ, α₁, …, α_p, β₁^T, …,β_q^T)^T，备择假设H₁：其中至少有一个τ_i≠0。在式(13)中得到相应参数的极大似然估计值，全模型的对数似然函数式(12)对应的极值问题转换为

$ \begin{gather*} \sup\limits_{\tau, \sigma^{2}} \ln L\left(\boldsymbol{\tau}, \sigma^{2}\right)=-\frac{n}{2}\left[\ln \frac{2 \mathsf{π}}{n}-\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|+\right. \\ \left.\ln (\boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}){ }^{\mathrm{T}} \boldsymbol{S}^{\mathbf{T}} \boldsymbol{S} \boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}+1\right], \end{gather*} $

(18)

其中, $\boldsymbol{S}=\boldsymbol{I}_{n}-\boldsymbol{B}\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{\mathrm{T}}, \boldsymbol{B}=(\boldsymbol{X}, \tilde{\boldsymbol{Z}}), \tilde{\boldsymbol{\phi}}$为H₀不成立时ϕ的极大似然估计量。

当τ=0时，代入式(11)得到约简模型

$ \boldsymbol{Y}=\boldsymbol{\beta}_{0}+\boldsymbol{\varepsilon}, $

其中β₀=(β₀(u₁, v₁), …, β₀(u_n, v_n))^T, 由式(10)可知$\widetilde{\boldsymbol{\beta}_{0}}\left(u_{i}, v_{i}\right)=\boldsymbol{\gamma}_{0}^{\mathrm{T}} \boldsymbol{b}\left(u_{i}, v_{i}\right)$。记 B₀=(b(u₁, v₁), …, b(u_n, v_n))^T, 故零假设成立时的对数似然函数为

$ \begin{gather*} \ln L\left(\boldsymbol{\beta}_{0}, \sigma^{2}\right)=\ln L\left(\boldsymbol{\gamma}_{0}, \sigma^{2}\right)=-\frac{n}{2} \ln \left(2 \mathsf{π} \sigma^{2}\right)- \\ \frac{1}{2 \sigma^{2}}\left(\boldsymbol{Y}-\boldsymbol{B}_{0} \boldsymbol{\gamma}_{0}\right)^{\mathrm{T}}\left(\boldsymbol{Y}-\boldsymbol{B}_{0} \boldsymbol{\gamma}_{0}\right). \end{gather*} $

(19)

则

$ \sup\limits_{\boldsymbol{\beta}_{0}, \sigma^{2}} \ln L\left(\boldsymbol{\beta}_{0}, \sigma^{2}\right)=-\frac{n}{2}\left[\ln \frac{2 \mathsf{π}}{n}+\ln \boldsymbol{Y}^{\mathrm{T}} \boldsymbol{S}_{0}^{\mathrm{T}} \boldsymbol{S}_{0} \boldsymbol{Y}+1\right] $

(20)

其中S₀=I_n- B₀(B₀^T B₀)^-1B₀^T。根据式(18)和式(20)，定义似然比检验统计量

$ \begin{aligned} & \Lambda_{\mathrm{T}}(\boldsymbol{B})=2 \ln \frac{\sup\limits_{\tau, \sigma^{2}} L\left(\boldsymbol{\tau}, \sigma^{2}\right)}{\sup\limits_{\boldsymbol{\beta}_{0}, \sigma^{2}} L\left(\boldsymbol{\beta}_{0}, \sigma^{2}\right)}= \\ & -n\left[\ln (\boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y})^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}-\right. \\ & \left.\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|-\ln \boldsymbol{Y}^{\mathrm{T}} \boldsymbol{S}_{0}^{\mathrm{T}} \boldsymbol{S}_{0} \boldsymbol{Y}\right]. \end{aligned} $

当样本量n足够大时，Λ_T(B)→χ²(p+m(q+1))。对给定的显著性水平a，当Λ_T(B)>χ_a²(p+m(q+1))时，拒绝原假设H₀。

2.2 模型系数的显著性检验

对于上述整体性检验，如果检验的结果是拒绝原假设，这意味着因变量Y的存在依赖于自变量X和 Z以及受到因变量空间滞后项的影响，因此这里需要对引起 Y变化的每个自变量逐一作显著性检验。对空间滞后项对应的变系数ρ进行显著性检验是对因变量Y滞后项是否存在独立影响Y的证明，需要作单独检验。类似地，对α和β进行显著性检验也是为了验证自变量 X和Z对Y是否存在影响。

1) 空间滞后项对应变系数的显著性检验

考虑零假设H₀:ρ=0的检验问题，对于在式(10)中ρ的核插值估计，检验问题近似等价于

$ H_{0}^{*}: \phi_{1}=\cdots=\phi_{m}=0 \text { vs } H_{1}^{*}: $

至少有一个ϕ_i≠0, i=1, …, m.

假设H₀^*成立，即模型不存在空间滞后项，此时模型退化为混合地理加权回归模型，记 η=(α^T, γ^T)^T，相应的对数似然函数为

$ \ln L\left(\boldsymbol{\eta}, \sigma^{2}\right)=-\frac{n}{2} \ln 2 \mathsf{π} \sigma^{2}-\frac{[\boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}]^{\mathrm{T}}[\boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}]}{2 \sigma^{2}}. $

(21)

则

$ \begin{equation*} \sup\limits_{\boldsymbol{\eta}, \sigma^{2}} \ln L\left(\boldsymbol{\eta}, \sigma^{2}\right)=-\frac{n}{2}\left[\ln \frac{2 \mathsf{π}}{n}+\ln \boldsymbol{Y}^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{Y}+1\right] . \end{equation*} $

(22)

其中S=I_n- B(B^TB)^-1 B^T。该检验所对应的全模型存在时的极大似然估计解为式(13)，由式(18)及式(22)，定义检验统计量

$ \begin{gathered} \Lambda_\rho(\boldsymbol{B})=2 \ln \frac{\sup\limits_{\phi, \boldsymbol{\eta}, \sigma^2} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}, \sigma^2\right)}{\sup\limits_{\boldsymbol{\eta}, \sigma^{\mathrm{2}}} L\left(\boldsymbol{\eta}, \sigma^2\right)}= \\ -n\left[\ln (\boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y})^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{ST}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}-\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|-\right. \\ \left.\ln \boldsymbol{Y}^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{Y}\right] . \end{gathered} $

其中$\widetilde{\boldsymbol{\phi}}$为全模型存在时ϕ的极大似然估计值。当样本量n足够大时，Λ_ρ(B)→χ²(m)。对给定的显著性水平a，当Λ_ρ(B)>χ_a²(m)时，拒绝原假设H₀。

2) 自变量对应常系数的显著性检验

下面考虑对于不同自变量对应的系数进行逐个显著性检验。首先讨论常系数α_i, i=1, …, p的显著性，记 α_-i=(α₁, …, α_i-1, α_i+1, …, α_p)^T，$\boldsymbol{B}_{-i}=\left(\boldsymbol{X}_{-i}, \tilde{\boldsymbol{Z}}\right)$, 其中X_-i为 X删去系数α_i对应的自变量列后的n×(p-1) 维矩阵，i=1, …p。考虑原假设H_0i: α_i=0，此时原假设下模型对应参数为η_-i= (α_-i^T,γ^T)^T，ϕ及σ_-i²，其中η_-i和σ_-i²的约束极大似然估计值与式(13)相似，分别为

$ \left\{\begin{array}{l} \widetilde{\boldsymbol{\eta}}_{-i}=\left(\boldsymbol{B}_{-i}^{\mathrm{T}} \boldsymbol{B}_{-i}\right)^{-1} \boldsymbol{B}_{-i}^{\mathrm{T}} \boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-i}\right) \boldsymbol{Y} , \\ \widetilde{\sigma_{-i}^{2}}=\frac{1}{n}\left(\boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-i}\right) \boldsymbol{Y}\right)^{\mathrm{T}} \boldsymbol{S}_{-i}^{\mathrm{T}} \boldsymbol{S}_{-i} T\left(\widetilde{\boldsymbol{\phi}}_{-i}\right) \boldsymbol{Y}. \end{array}\right. $

(23)

其中$\boldsymbol{S}_{-i}=\boldsymbol{I}_{n}-\boldsymbol{B}_{-i}\left(\boldsymbol{B}_{-i}^{\mathrm{T}} \boldsymbol{B}_{-i}\right)^{-1} \boldsymbol{B}_{-i}^{\mathrm{T}}, \widetilde{\boldsymbol{\phi}}_{-i}$为H_0i成立时 ϕ的极大似然估计值，结合式(18)，此时定义似然比统计量为

$ \begin{gathered} \Lambda_{-i}(\boldsymbol{B})=2 \ln \frac{\sup\limits_{\phi, \boldsymbol{\eta}, \sigma^{2}} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}, \boldsymbol{\sigma}^{2}\right)}{\sup\limits_{\boldsymbol{\phi}, \boldsymbol{\eta}_{-i}, \sigma^{2}} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}_{-i}, \sigma^{2}\right)}= \\ -n\left[\ln (\boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}){ }^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{ST}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}-\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|-\right. \\ \left.\ln \left(\boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-i}\right) \boldsymbol{Y}\right)^{\mathrm{T}} \boldsymbol{S}_{-i}^{\mathrm{T}} \boldsymbol{S}_{-i} \boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-i}\right) \boldsymbol{Y}+\frac{2}{n} \ln \left|\boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-i}\right)\right|\right]. \end{gathered} $

$\widetilde{\boldsymbol{\phi}}$为全模型存在时ϕ的极大似然估计值。当样本量n足够大时，Λ_-i(B)→χ²(1)。对给定的显著性水平a，当Λ_-i(B)>χ_a²(1)时，拒绝原假设H_0i。

3) 自变量对应变系数的显著性检验

与自变量对应常系数的显著性检验类似，定义$\boldsymbol{\beta}_{-j}=\left(\boldsymbol{\beta}_{1}, \cdots, \boldsymbol{\beta}_{j-1}, \boldsymbol{\beta}_{j+1}, \cdots, \boldsymbol{\beta}_{q}\right)^{\mathrm{T}}, \boldsymbol{B}_{-j}=(\boldsymbol{X}$, $\left.\tilde{\boldsymbol{Z}}_{-j}\right)$, 其中$\tilde{\boldsymbol{Z}}_{-j}=\left(\tilde{\boldsymbol{Z}}_{0}, \tilde{\boldsymbol{Z}}_{1}, \cdots, \tilde{\boldsymbol{Z}}_{j-1}, \tilde{\boldsymbol{Z}}_{j+1}, \cdots, \tilde{\boldsymbol{Z}}_{q}\right)$为$\tilde{\boldsymbol{Z}}$删去系数β_j对应的自变量列$\tilde{\boldsymbol{Z}}_{j}$后的n×mq维矩阵，j=1, …, q。

令γ=(γ₀^T, γ₁^T, …, γ_q^T)^T，γ_-j= (γ₀^T,γ₁^T, …,γ_j-1^T, γ_j+1^T, …, γ_q^T)^T，由式(10)中 β的核插值估计, 原假设H_0j:β_j=0近似等价于H_0j^*:γ_j=0。此时原假设下模型对应参数为 η_-j=(α^T, γ_-j^T)^T，ϕ及σ_-j²，其中η_-j和σ_-j²的约束极大似然估计值分别为

$ \left\{\begin{array}{l} \widetilde{\boldsymbol{\eta}}_{-j}=\left(\boldsymbol{B}_{-j}^{\mathrm{T}} \boldsymbol{B}_{-j}\right)^{-1} \boldsymbol{B}_{-j}^{\mathrm{T}} T\left(\widetilde{\boldsymbol{\phi}}_{-j}\right) \boldsymbol{Y}, \\ \widetilde{\sigma_{-j}^{2}}=\frac{1}{n}\left(\boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-j}\right) \boldsymbol{Y}\right){ }^{\mathrm{T}} \boldsymbol{S}_{-j}^{\mathrm{T}} \boldsymbol{S}_{-j} \boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}_{-j}\right) \boldsymbol{Y}. \end{array}\right. $

(24)

其中$\boldsymbol{S}_{-j}=\boldsymbol{I}_{n}-\boldsymbol{B}_{-j}\left(\boldsymbol{B}_{-j}^{\mathrm{T}} \boldsymbol{B}_{-j}\right)^{-1} \boldsymbol{B}_{-j}^{\mathrm{T}}, \widetilde{\boldsymbol{\phi}}_{-j}$为H_0j成立时 ϕ的极大似然估计值，此时定义似然比统计量为

$ \begin{array}{*{20}{L}} {{\Lambda _{ - j}}({\boldsymbol{B}}) = 2\ln \frac{{\mathop {\sup }\limits_{\boldsymbol{\phi} ,{\boldsymbol{\eta }},{\sigma ^2}} L\left( {\boldsymbol{\phi} ,{\boldsymbol{\eta }},{\sigma ^2}} \right)}}{{\mathop {\sup }\limits_{\boldsymbol{\phi} ,{{\boldsymbol{\eta }}_{ - j}},{\sigma ^2}} L\left( {\boldsymbol{\phi} ,{{\boldsymbol{\eta }}_{ - j}},{\sigma ^2}} \right)}} = }\\ { - n\left[ {\ln {{({\boldsymbol{T}}(\widetilde{\boldsymbol{\phi}} ){\boldsymbol{Y}})}^{\rm{T}}}{{\boldsymbol{S}}^{\rm{T}}}{\boldsymbol{ST}}(\widetilde{\boldsymbol{\phi}} ){\boldsymbol{Y}} - \frac{2}{n}\ln |{\boldsymbol{T}}(\widetilde{\boldsymbol{\phi}} )| - } \right.}\\ {\ln {{\left( {{\boldsymbol{T}}\left( \widetilde{\boldsymbol{\phi}}_{-j} \right){\boldsymbol{Y}}} \right)}^{\rm{T}}}{\boldsymbol{S}}_{ - j}^{\rm{T}}{{\boldsymbol{S}}_{ - j}}{\boldsymbol{T}}\left( \widetilde{\boldsymbol{\phi}}_{-j} \right){\boldsymbol{Y}} + \frac{2}{n}\ln \left| {{\boldsymbol{T}}\left( \widetilde{\boldsymbol{\phi}}_{-j} \right)} \right|.} \end{array} $

$\tilde{\boldsymbol{\phi}}$为全模型成立时ϕ的极大似然估计值。当样本量n足够大时，Λ_-j(B)→χ²(m)。对给定的显著性水平a，当Λ_-j(B)>χ_a²(m)时，拒绝原假设H_0j。

2.3 模型系数的空间局域性判别检验

由于常系数都可以看作是变系数在每个空间点上的系数值相等，反映了空间全局性的特征，MGWR模型、MGWR-CSLR模型等都是本文中所给出MGWR-VSLR模型的特殊情形，因此，检验估计的系数是否具有局域性而选取变系数形式，究竟应该选取常系数还是变系数对于设定适宜的模型形式是非常必要的。

1) 空间滞后项对应系数的空间局域性检验

设置零假设H₀：ρ(u_i, v_i)≡ρ₀，$ \forall\left(u_{i}, v_{i}\right) \in \mathbb{R}^{2}$，当H₀成立时，相应系数为空间全局性的常系数，MGWR-VSLR模型退化为MGWR-CSLR模型, 如下所示

$ \begin{equation*} \boldsymbol{Y}=\rho_{0} \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\alpha}+\boldsymbol{M}+\boldsymbol{\varepsilon} . \end{equation*} $

(25)

记C(ρ₀)=I-ρ₀W，与前文相同，取$\boldsymbol{B}=(\boldsymbol{X}, \tilde{\boldsymbol{Z}})$以及 η=(α^T, γ^T)^T，上式对应的对数似然函数为

$ \begin{gathered} \ln L\left(\rho_{0}, \boldsymbol{\eta}, \sigma^{2}\right)=-\frac{n}{2} \ln 2 \mathsf{π} \sigma^{2}+\ln \left|\boldsymbol{C}\left(\rho_{0}\right)\right|- \\ \frac{1}{2 \sigma^{2}}\left\{\left[\boldsymbol{C}\left(\rho_{0}\right) \boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}\right]^{\mathrm{T}}\left[\boldsymbol{C}\left(\rho_{0}\right) \boldsymbol{Y}-\boldsymbol{B} \boldsymbol{\eta}\right]\right\}. \end{gathered} $

此时上式对应的极值可写为

$ \begin{gather*} \sup\limits_{\rho_{0}, \boldsymbol{\eta}, \sigma^{2}} \ln L\left(\rho_{0}, \boldsymbol{\eta}, \sigma^{2}\right)=-\frac{n}{2}\left[\ln \frac{2 \mathsf{π}}{n}-\ln \left|\boldsymbol{C}\left(\tilde{\rho}_{0}\right)\right|+\right. \\ \left.\ln \left(\boldsymbol{C}\left(\tilde{\rho}_{0}\right) \boldsymbol{Y}\right)^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{C}\left(\tilde{\rho}_{0}\right) \boldsymbol{Y}+1\right] . \end{gather*} $

(26)

其中$\tilde{\rho}_{0}$为ρ₀在式(25)下的极大似然估计值，S=I_n- B。因此可以定义似然比统计量

$ \begin{gathered} \Lambda_{c p}(\boldsymbol{B})=2 \ln \frac{\sup\limits_{\phi, \boldsymbol{\eta}, \sigma^{2}} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}, \sigma^{2}\right)}{\sup\limits_{\rho_{0}, \boldsymbol{\eta}, \sigma^{2}} L\left(\rho_{0}, \boldsymbol{\eta}, \sigma^{2}\right)}= \\ -n\left[\ln (\boldsymbol{T}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}){ }^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{ST}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}-\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|-\right. \\ \left.\ln \left(\boldsymbol{C}\left(\tilde{\boldsymbol{\rho}}_{0}\right) \boldsymbol{Y}\right){ }^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{C}\left(\tilde{\boldsymbol{\rho}}_{0}\right) \boldsymbol{Y}+\frac{2}{n} \ln \left|\boldsymbol{C}\left(\tilde{\boldsymbol{\rho}}_{0}\right)\right|\right]. \end{gathered} $

其中$\tilde{\boldsymbol{\phi}}$为全模型成立时ϕ的极大似然估计值。当样本量n足够大时，Λ_cρ(B)→χ²(m-1)。对给定的显著性水平a，当Λ_cρ(B)>χ_a²(m-1)时，拒绝原假设H₀，空间滞后项对应的系数形式选取变系数。

2) 自变量对应系数的空间局域性检验

考虑零假设$H_{0 k}: \beta_{k}(u, v) \equiv \beta_{k 0}, \forall(u, v) \in$ $\mathbb{R}^{2}, k=0, 1, \cdots, q$，当假设H_0k成立时，相应系数为空间全局性的常系数，模型变为

$ \begin{equation*} \boldsymbol{Y}=\boldsymbol{\rho}_{n} \boldsymbol{W} \boldsymbol{Y}+\boldsymbol{X} \boldsymbol{\alpha}+\beta_{k 0} \boldsymbol{Z}_{k}+\boldsymbol{M}_{-k}+\varepsilon. \end{equation*} $

(27)

令$\boldsymbol{Z}_{-k}=\left[\boldsymbol{1}, \boldsymbol{Z}_{1}, \cdots, \boldsymbol{Z}_{k}, \cdots, \boldsymbol{Z}_{k+1}, \cdots \boldsymbol{Z}_{q}\right], \boldsymbol{\beta}_{-k}=\left[\boldsymbol{\beta}_{0}\right.$, $\left.\cdots, \boldsymbol{\beta}_{k-1}, \boldsymbol{\beta}_{k+1}, \cdots, \boldsymbol{\beta}_{q}\right]^{\mathrm{T}}, \left\{d_{-k, 1}, \cdots, d_{-k, n}\right\}$为$\boldsymbol{Z}_{-k} \boldsymbol{\beta}_{-k}$的对角线元素, 则$\boldsymbol{M}_{-k}=\operatorname{diag}\left\{d_{-k, 1}, \cdots\right.$, $\left.d_{-k, n}\right\} \boldsymbol{1}_{n}$。记$\boldsymbol{B}^{*}=\left(\boldsymbol{X}, \boldsymbol{Z}_{k}, \tilde{\boldsymbol{Z}}_{-k}\right)$, 其中$\tilde{\boldsymbol{Z}}_{-k}=\left(\tilde{\boldsymbol{Z}}_{0}\right.$, $\left.\tilde{\boldsymbol{Z}}_{1}, \cdots, \tilde{\boldsymbol{Z}}_{k-1}, \tilde{\boldsymbol{Z}}_{k+1}, \cdots, \tilde{\boldsymbol{Z}}_{q}\right), \tilde{\boldsymbol{Z}}_{j}=\left(\boldsymbol{b}\left(u_{1}, v_{1}\right) z_{1 j}, \cdots\right.$, $\left.\boldsymbol{b}\left(u_{n}, v_{n}\right) z_{n j}\right)^{\mathrm{T}}, j=1, \cdots, q$。另记$\boldsymbol{\eta}^{*}=\left(\boldsymbol{\alpha}^{\mathrm{T}}, \boldsymbol{\beta}_{i 0}\right.$, $\left.\boldsymbol{\gamma}_{-k}^{\mathrm{T}}\right)^{\mathrm{T}}, \boldsymbol{\gamma}_{-k}=\left(\boldsymbol{\gamma}_{0}^{\mathrm{T}}, \boldsymbol{\gamma}_{1}^{\mathrm{T}}, \cdots, \boldsymbol{\gamma}_{k-1}^{\mathrm{T}}, \boldsymbol{\gamma}_{k+1}^{\mathrm{T}}, \cdots, \boldsymbol{\gamma}_{q}^{\mathrm{T}}\right)^{\mathrm{T}}$, 则$\boldsymbol{B}^{*} \boldsymbol{\eta}^{*}=\boldsymbol{X} \boldsymbol{\alpha}+\boldsymbol{\beta}_{k 0} \boldsymbol{Z}_{k}+\boldsymbol{M}_{-k}$, 则式(27) 可以简化为

$ \boldsymbol{Y}=\boldsymbol{T}(\boldsymbol{\phi})^{-1}\left(\boldsymbol{B}^{*} \boldsymbol{\eta}^{*}+\boldsymbol{\varepsilon}\right) . $

上式的对数似然函数为

$ \begin{gathered} \ln L\left(\boldsymbol{\phi}, \boldsymbol{\eta}^{*}, \sigma^{2}\right)=-\frac{n}{2} \ln 2 \mathsf{π} \sigma^{2}+\ln |\boldsymbol{T}(\boldsymbol{\phi})|- \\ \frac{1}{2 \sigma^{2}}\left\{\left[\boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}-\boldsymbol{B}^{*} \boldsymbol{\eta}^{*}\right]^{\mathrm{T}}\left[\boldsymbol{T}(\boldsymbol{\phi}) \boldsymbol{Y}-\boldsymbol{B}^{*} \boldsymbol{\eta}^{*}\right]\right\}. \end{gathered} $

此时对应的极值可写为

$ \begin{array}{r} \sup\limits_{\boldsymbol{\phi}, \boldsymbol{\eta}^{*}, \sigma^{2}} \ln L\left(\boldsymbol{\phi}, \boldsymbol{\eta}^{*}, \sigma^{2}\right)=-\frac{n}{2}\left[\ln \frac{2 \mathsf{π}}{n}-\ln \left|\boldsymbol{T}\left(\tilde{\boldsymbol{\phi}}^{*}\right)\right|+\right. \\ \left.\ln \left(\boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}^{*}\right) \boldsymbol{Y}\right){ }^{\mathrm{T}} \boldsymbol{S}^{* \mathrm{~T}} \boldsymbol{S}^{*} \boldsymbol{T}\left(\widetilde{\boldsymbol{\phi}}^{*}\right) \boldsymbol{Y}+1\right] . \end{array} $

(28)

其中$\tilde{\boldsymbol{\phi}}^{*}$为$\boldsymbol{\phi}$在模型(27) 下的极大似然估计值, $\boldsymbol{S}^{*}=\boldsymbol{I}_{n}-\boldsymbol{B}^{*}\left(\boldsymbol{B}^{* \mathrm{~T}} \boldsymbol{B}^{*}\right)^{-1} \boldsymbol{B}^{* \mathrm{~T}}$。根据式(18) 及式(28), 定义似然比统计量

$ \begin{gathered} \Lambda_{c k}(\boldsymbol{B})=2 \ln \frac{\sup\limits_{\boldsymbol{\phi, \eta}^{*}, \sigma^{2}} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}, \sigma^{2}\right)}{\sup\limits_{\boldsymbol{\phi, \eta}^{*}, \sigma^{2}} L\left(\boldsymbol{\phi}, \boldsymbol{\eta}^{*}, \sigma^{2}\right)}= \\ -n\left[\ln (\boldsymbol{T}(\tilde{\boldsymbol{\phi}}) \boldsymbol{Y})^{\mathrm{T}} \boldsymbol{S}^{\mathrm{T}} \boldsymbol{ST}(\widetilde{\boldsymbol{\phi}}) \boldsymbol{Y}-\frac{2}{n} \ln |\boldsymbol{T}(\widetilde{\boldsymbol{\phi}})|-\right. \\ \left.\ln \left(\boldsymbol{T}\left(\tilde{\boldsymbol{\phi}}^{*}\right) \boldsymbol{Y}\right)^{\mathrm{T}} \boldsymbol{S}^{* \mathrm{~T}} \boldsymbol{S}^{*} \boldsymbol{T}\left(\tilde{\boldsymbol{\phi}}^{*}\right) \boldsymbol{Y}+\frac{2}{n} \ln \left|\boldsymbol{T}\left(\tilde{\boldsymbol{\phi}}^{*}\right)\right|\right]. \end{gathered} $

其中$\widetilde{\boldsymbol{\phi}}$为ϕ在全模型下的极大似然估计值。当样本量n足够大时，Λ_ck(B)→χ²(m-1)。对给定显著性水平a，当Λ_ck(B)>χ_a²(m-1)时，拒绝原假设H_0k，自变量对应的系数形式选取变系数。

3 蒙特卡罗模拟与实际算例应用 3.1 蒙特卡罗模拟的方案设定

乔宁宁^[18]通过数值模拟实验认为MGWR-CSLR模型比传统常系数空间滞后模型在空间异质性处理上表现出更大的优越性。为考察MGWR-VSLR模型的参数估计效果和检验统计量，本文同样基于数值模拟实验在比较MGWR-VSLR模型与MGWR-CSLR模型的同时，也引入了VSLR和GWR-VSLR 2种变系数的空间滞后回归模型，通过比较系数估计值的均方根误差(root mean square error，RMSE)及因变量预测值的RMSE，用以对比不同样本容量下系数估计值偏误和因变量预测能力的差异，4种模型的具体形式如下：

a) VSLR模型

$ y_{i}=\rho\left(u_{i}, v_{i}\right) \sum\nolimits_{j=1}^{n} w_{i j} y_{j}+\alpha+\beta x_{i}+\kappa z_{i}+\varepsilon_{i}, $

b) GWR-VSLR模型

$ \begin{aligned} y_{i}= & \rho\left(u_{i}, v_{i}\right) \sum\nolimits_{j=1}^{n} w_{i j} y_{j}+\alpha\left(u_{i}, v_{i}\right)+ \\ & \beta\left(u_{i}, v_{i}\right) x_{i}+\kappa\left(u_{i}, v_{i}\right) z_{i}+\varepsilon_{i}, \end{aligned} $

c) MGWR-CSLR模型

$ \begin{gathered} y_{i}=\rho \sum\nolimits_{j=1}^{n} w_{i j} y_{j}+\alpha\left(u_{i}, v_{i}\right)+\beta x_{i}+ \\ \kappa\left(u_{i}, v_{i}\right) z_{i}+\varepsilon_{i}, \end{gathered} $

d) MGWR-VSLR模型

$ \begin{gathered} y_{i}=\rho\left(u_{i}, v_{i}\right) \sum\nolimits_{j=1}^{n} w_{i j} y_{j}+\alpha\left(u_{i}, v_{i}\right)+ \\ \beta x_{i}+\kappa\left(u_{i}, v_{i}\right) z_{i}+\varepsilon_{i} . \end{gathered} $

在数值模拟的方案中, 设定$\boldsymbol{Y}=\left(y_{1}, \cdots, y_{n}\right)^{\mathrm{T}}$为因变量, $\rho \sum_{j=1}^{n} w_{i j} y_{j}$表示常系数空间滞后项, 而$\rho\left(u_{i}, v_{i}\right) \sum_{j=1}^{n} w_{i j} y_{j}$表示变系数空间滞后项, $\beta$表示自变量$\boldsymbol{X}=\left(x_{1}, \cdots, x_{n}\right)^{\mathrm{T}}$对应的全局性系数, $\beta\left(u_{i}, v_{i}\right)$则为自变量$\boldsymbol{X}=\left(x_{1}, \cdots, x_{n}\right)^{\mathrm{T}}$对应的局域性系数, $\alpha$和$\kappa$分别表示截距项和自变量$\boldsymbol{Z}=$ $\left(z_{1}, \cdots, z_{n}\right)^{\mathrm{T}}$对应的全局性系数, 而$\alpha\left(u_{i}, v_{i}\right)$和$\kappa\left(u_{i}, v_{i}\right)$则分别为截距项和自变量$\boldsymbol{Z}=\left(z_{1}, \cdots\right.$, $\left.z_{n}\right)^{\mathrm{T}}$对应的局域性系数, $\varepsilon_{i} \sim N\left(0, \sigma^{2}\right)$。设定研究空间点$i$的变量二维坐标空间$\left(u_{i}, v_{i}\right)$且横坐标与纵坐标均在$100 \mathrm{~km}$范围内, 以0值为坐标原点且以$100 \mathrm{~km}$为坐标单位, 则$u_{i}, v_{i} \in[0, 1]$。在VSLR模型中$\alpha=0.5$作为固定的值, $\boldsymbol{W}$作为空间权重矩阵其权重元素通过不同空间位置的欧式距离倒数获得。考虑到空间存在的线性与非线性关系，这里假设不同空间点$i$的可变系数$\rho\left(u_{i}, v_{i}\right)$, $\beta\left(u_{i}, v_{i}\right)$和$\kappa\left(u_{i}, v_{i}\right)$分别满足3种关于$u_{i}, v_{i}$的空间关系, 即$\rho\left(u_{i}, v_{i}\right)=\left(u_{i}-0.5\right)^{2}+v_{i}^{2} ; \beta\left(u_{i}, v_{i}\right)=$ $0.5\left(u_{i}-2 v_{i}\right) ; \kappa\left(u_{i}, v_{i}\right)=\exp \left\{0.6 u_{i}\right\}+\sin \left(5 u_{i} v_{i}\right)$。在实验中$\boldsymbol{Y}, \boldsymbol{X}, \boldsymbol{Z}$变量的观测数值由均匀分布$U(0, 1)$随机生成, 分别独立生成$n=200, 500$, 1 000和3 000的样本量并重复模拟$N=50$次, 同样由均匀分布$U(0, 1)$随机生成的测试集样本量为$n_{\text {test }}=n / 10$, 在重构参数化求解的模拟中取插值的节点数$m=[4 \ln 2 n], []$为向下取整函数。通过模型a) 至d) 得到不同模型形式的系数估计值分别为$\tilde{\alpha}, \tilde{\alpha}\left(u_{i}, v_{i}\right), \tilde{\beta}, \tilde{\beta}\left(u_{i}, v_{i}\right), \tilde{\kappa}, \tilde{\kappa}\left(u_{i}, v_{i}\right), \rho$和$\tilde{\rho}\left(u_{i}\right.$, $\left.v_{i}\right)$。

3.2 不同样本容量下4种模型的模拟结果与分析

采用Xiong^[20]的重构参数化方法在n=200, 500, 1 000, 及3 000不同样本容量下的数值实验模拟中，对上述4种模型进行相关系数估计，其结果为表 1至表 4，对比发现：

1) 当截距项对应的是全局性系数(VSLR模型)时，系数估计的偏误明显高于截距项对应的是局域性系数(GWR-VSLR模型、MGWR-CSLR模型和MGWR-VSLR模型)，这一特征在不同样本容量下均表现出一致性。由此可以看出，在处理空间异质性问题时，将截距项对应的系数设为局域性系数比全局性系数有助于减小系数估计的偏误。

2) 在不同样本容量下的数值实验模拟中，4种模型按照因变量预测值的RMSE由大到小依次排序为VSLR模型>MGWR-CSLR模型>GWR-VSLR模型>MGWR-VSLR模型，即变系数的空间滞后回归模型对于因变量的预测能力明显弱于其他3种模型。可以看出在处理空间异质性问题上，除了考虑因变量自身存在的空间相关性，自变量的空间异质性也决定着因变量的预测能力，自变量也需要考虑是否存在局域性自变量，所对应的系数用以反映局域性特征，提高模型整体的预测能力。

3) 进一步对比存在局域自变量的GWR-VSLR模型、MGWR-CSLR模型和MGWR-VSLR模型发现，MGWR-CSLR模型的自变量X所对应系数估计值$\tilde{\beta}\left(u_{k}, v_{k}\right)$的偏误小于GWR-VSLR模型，但是MGWR-CSLR模型因变量空间滞后项对应的变系数估计值$\tilde{\rho}\left(u_{k}, v_{k}\right)$的偏误大于GWR-VSLR模型。由此可以看出自变量X与Z所对应系数并不都适合反映局域性特征，这也是MGWR相比GWR具有的优势。同时可以看出GWR-VSLR模型采用反映因变量空间滞后项的变系数对整个模型的拟合能力也起着重要作用。

4) 相较于VSLR模型、GWR-VSLR模型和MGWR-CSLR模型，本文提出的MGWR-VSLR模型在不同样本容量的数值模拟实验中，各项系数估计的偏误整体上表现出较小的特征。并且随着样本容量的增加，系数估计偏误有所下降，测试集上因变量预测值的RMSE也在下降，体现了估计的稳定性。其他模型的系数估计偏误与因变量预测值的RMSE也随着样本容量增大有所下降，说明本文中的随机试验不同模型的模拟结果对比具有一定的可参考性。

3.3 一个实际算例的应用

经济发展与人口总数和海拔高度具有一定相关性。一般来说，一个区域人口越多，存在的经济活动也就越多，而海拔越高，交通条件受限导致经济发展的自然环境相对越差。考察中国2015年人口数与海拔高度对国内生产总值(gross domestic product, GDP)的影响，这里随机抽取1 200个栅格样本数据，其中GDP作为因变量，人口数POP(population)和海拔高度利用数字高程模型(digital elevation model, DEM)数据作为自变量，每个栅格数据的地理坐标为(u_i, v_i)，它们均匀覆盖中国31省区市，代表 1 km²范围的平均状况。选取1 000个栅格数据用于模型拟合，其余200个栅格数据则作为检验模型的预测能力。数据来源于中国科学院地理科学与资源研究所资源环境科学与数据中心(https://www.resdc.cn/)。从栅格单位面积的人口数来看，POP即人口密度，相较于人口密度，交通状况与经济增长更直接相关，而交通状况直接受制于海拔地形等自然条件约束。因此，可以考虑在MGWR-VSLR模型中把DEM变量作为局域变量，POP作为全局变量，同时各地区的经济发展与周边地区经济相关性存在着差异，采用变系数空间滞后性也具有一定合理性。对于不同空间点的GDP，建立MGWR-VSLR模型的具体计算公式为

$ \begin{align*} \mathrm{GDP}_{i}= & \rho\left(u_{i}, v_{i}\right) \sum\nolimits_{i=1}^{n} w_{i j} \mathrm{GDP}_{j}+\beta_{0}\left(u_{i}, v_{i}\right)+ \\ & \beta_{1} \mathrm{POP}_{i}+\beta_{2}\left(u_{i}, v_{i}\right) \mathrm{DEM}_{i}+\varepsilon_{i}, \end{align*} $

(29)

其中空间权重矩阵W=(w_ij)_{i, j=1}ⁿ同样采用欧式距离倒数获得，ρ(u_i, v_i), β₂(u_i, v_i), β₁和β₀(u_i, v_i) 为相关变量的系数和截距项。上述回归分析除采用经典的OLSR模型和GWR模型外，也引入了MGWR-CSLR模型，进行4种模型整体性能的对比(表 5)。由表 5可以看出，MGWR-VSLR模型通过重构参数化的方法在GDP数据上的拟合优度高于其他3种模型，同时AIC值和均方残差最小，说明MGWR-VSLR模型估计方法更加精确。另外，在模型的预测方面，MGWR-VSLR模型的均方预测误差最小，说明该方法能够更好地预测不同空间单元人口和海拔高度对GDP的影响。

4 结论

空间计量模型在国民经济社会发展的今天得到了广泛应用，尤其是在面对处理空间相关性和空间异质性等实际的经济地理学问题时，备受学者们的关注。本研究在现有GWR各种模型形式的研究基础上，提出一类新的MGWR-VSLR模型，该模型不仅将因变量的空间相关性与自变量的空间异质性同时纳入到模型中考虑，同时刻画了因变量的空间局域性特征。研究结论可以概括如下：

1) MGWR-VSLR模型对于处理空间异质性问题具有更广泛的适用性。在模型形式的构成上，MGWR-VSLR模型涵盖了GWR的多种模型形式，不仅反映了自变量存在空间全局性与局域性的特征，并且进一步考虑了因变量空间局域性的特征。在蒙特卡洛模拟与实际算例应用的结果中显示出MGWR-VSLR模型较好的因变量拟合和预测能力，在样本量增大的情况下也呈现出可靠性与稳健性。

2) MGWR-VSLR模型能够更加广泛地涵盖GWR的多种模型形式，并将MGWR-CSLR模型做了进一步扩展。基于似然比检验给出MGWR-VSLR模型系数估计的显著性检验，蒙特卡罗模拟结果显示，显著性检验在样本量增大的情况下具有稳健性。同时也给出了空间全局性常系数与局域性变系数选取的统计检验方法，为空间计量模型的合理设置提供了依据。

概言之，本文提出的MGWR-VSLR模型，能够为处理空间效应(空间相关性和空间异质性)问题提供一种更加全面的、更具有一般性的定量分析方法，对于常系数与变系数的显著性检验乃至选取的统计检验，也有望推广到其他非参数空间模型的研究中。