一、空间非平稳性
社会科学家长期以来,面对一个困难的问题和一个潜在的困境:是否有任何“定理”支配着社会活动,如果没有,定量的方法还有任何有效性吗?这个问题从下面两个子问题来看更清楚。第一个问题是模拟社会科学的模型不是完全正确的,第二个问题是从一个系统得到的结果几乎很少能复制到另一个系统。物理活动常常是趋于平稳的(如著名的能量和质量公式E=MC2,不管采用什么度量方法,也不管你在哪个国家或者城市,这个公式都是成立的),而社会活动经常不是平稳的,一个关系的度量部分地依赖于所采用的度量方法。在空间活动下,我们称这种非平稳性为空间非平稳性。本质上,我们试图调查的空间活动在空间中可能不是保持不变的。
在建模时,参数可能会随着空间而变化。通常有以下三种可能原因导致空间非平稳性:第一是,随机抽样变异,如假定我们从一个数据集中取几个空间子数据集,然后分别刻画同一个模型,则我们不会期望所获得的参数估计值是完全相同的:由于使用不同的抽样数据,差异是存在的。第二是,不管什么原因,一些关系本质上在不同的空间是不同的。如由于人们的态度或者偏好的差异,或者由于不同管理的、政治的观点,在空间中对相同的刺激会产生不同的反应。第三是,用来估计关系的模型是一个对现实有误的描述,一个或更多相关的变量要么被忽略了,要么被一个不正确函数来表示了。
假定空间非平稳性是存在的,那么用像OLS这样的全局分析方法就无法反映非平稳性,因此有必要用局部的技术手段分析空间非平稳性。FOTHERINGHAM ET AL。(2002)指出分析空间数据的局部分析方法分为三类:单变量空间数据局部统计学、多变量空间数据局部统计学、运动的空间模式。
单变量空间数据局部分析方法有4种类型:点模式分析的局部、局部地理分析、局部过滤和空间依赖的局部度量。点模式分析的局部形式首先由奥彭肖等人(OPENSHAW ET AL。,1987)发展,然后由福瑟林厄姆和张(FOTHERINGHAM AND ZHANG,1996)修正的地理分析机器(GEOGRAPGICAL ANALYSIS MACHINE,GAM),GAM的基本思想是识别数据集的局部感兴趣部分,如识别疾病集聚。局部地理分析是随着可视化数据的技术发展而产生的,主要有克雷西(CRESS-IE,1984)的空间滞后散点图,阿斯莱等人(HASLETT ET AL。,1991)的变量云图(VARIOGRAM CLOUD PLOT),安瑟兰(ANSELIN,1996)的MORAN散点图。空间局部过滤是把图像处理中高通滤波器和低通滤波器的技术应用于空间分析,可从如下的学者的研究中找出局部过滤例子:施密德和麦卡坎内尔(SCHMID&MACCANNELL,1955),昂温(UNWIN,1981),郑(CHENG,1996)。空间依赖的局部度量是从空间依赖的全局度量基础上发展而来的,主要是研究空间集聚和空间自相关。许多学者对此进行了广泛的讨论:格蒂什和奥德(GETIS&ORD,1992),奥德和格蒂什(ORD&GETIS,1995,2001),安瑟兰(ANSELIN,1995,1998),罗杰森(ROGERSON,1999),罗森堡(ROSENBERG,2000),布伦斯坦(BRUNSDON,1998)。
多变量空间数据局部分析方法有:空间扩展方法(THE SPATIAL EXPANSION METHOD,CASETTI,1972,1997;JONES&CASETTI,1992;FOTHERINGHAM&PITTS,1995)、空间自适应过滤法(SPATIALLY ADAPTIVE FILTERING,WIDROW&HOFF,1960,TRIGG&LEACH,1968)、多层模型(MULTILEVEL MODELLING,GOLDSTEIN,1987)、随机系数模型(RANDOM COEFFICIENT MODELS,AIKIN,1997)、空间回归模型、地理加权回归(BRUNSDON ET AL。,1996,1999;FOTHERINGHAM ET AL。,1997B,1998,1999,2002)。
空间扩展方法通过设定参数为其他属性(如,地理坐标)的函数使得参数估计得以局部地变化,但是也有一些局限。该方法显示变量在空间上的变化趋势,却受限制于扩展方程的复杂性以及扩展方程的选择需要先验经验。空间自适应过滤法也被用于空间非平稳性,但是该方法采用一种相当特殊的方式,而且得到的参数估计值无法进行统计检验。多层模型(MULTILEVEL MODEL)和随机系数模型都假定所估计参数是随机变量,前者假定为高斯分布,而后者假定为有限混合分布。两种方法都是采用贝叶斯定理来估计参数,没有考虑到空间依赖性。空间回归模型是混合的模型,尽管认识到了数据间的局部关系,但是局部关系差不多是通过一个全局自相关统计量来度量的,而且模型的结果也是一个全局参数估计的集合,因此可以把它理解为“半局部”。地理加权回归模型把这些缺点都克服了。
二、地理加权回归方法的机制
地理加权回归(GWR)方法是1996年以后发展起来的,广泛用于处理空间非稳定性的统计方法。它吸引了下列学者的广泛注意:布伦斯坦等人(BRUNSDON ET AL。,1996,1999);福瑟林厄姆等人(FOTHERINGHAM ET AL。,1997B,1998,1999,2002);梁怡等人(LEUNG ET AL。,2000A,2000B);黄和梁(HUANG AND LEUNG,2002);佩茨等人(PAEZ ET AL。,2002A,2002B);余丹林等人(YU ET AL。,2004)。简而言之,GWR发展了卡塞蒂(CASETTI,1972,1992)扩展回归方法的思想。特别地,根据托布勒(TOBLER,1970)的地理学第一定理,GWR容许回归系数在空间中变化。
(一)基本方法
考虑如下的全局回归模型:
地理加权回归(GWR)扩展了传统的回归框架,容许局部而不是全局的参数估计,扩展后的模型如下:
其中,(UI,VI)是第I个样本点的空间坐标,βK(UI,VI)是连续函数βK(U,V)在I点的值。如果βK(U,V)在空间上保持不变,则模型(5.2)就变为全局模型(5.1)。因此,GWR方程(5.2)认可空间变化关系可能是存在的,并且提供了一种可度量的方法。
按照现在的情况,由于未知的变量比已观察到的变量还要多,在校准方程(5.2)时将会出现许多问题,然而这类模型在统计文献和讨论(ROSENBERG,1973,HASTIE&TIBSHIRANI,1990,LOADER,1999)中可以发现。福瑟林厄姆,布伦斯坦,查尔顿(FOTHERING-HAM,BRUNSDON AND CHARLTON)借鉴黑斯蒂和蒂希兰尼以及洛亚德尔的经验,假定系数不是随机的,而是其他一些变量(如,空间位置)的决定函数。用通常的方法处理像这样的模型时应该注意到,尽管局部系数的无偏估计是不可能的,但是只有小的偏差的参数估计还是可以找到的(此处,偏差产生于根据除了位置I以外所收集的数据,来推断在位置I处的非稳定过程的结果)。假设参数表现出一定程度的空间一致性,那么接近被估计位置的估计值有相似的大小和符号。因此,当对给定位置I的参数值进行估计时,可以用全局方程(5.1)来近似代替位置I的方程(5.2),使用和位置I相近的位置数据集来进行回归。对于其他位置的参数估计,也采用类似的方法,以此类推。
如上所述,方程(5.2)的刻画过程含蓄地假设:接近位置I的观察数据比那些离位置I远一些的数据对βK(UI,VI)的估计有更多的影响。故加权最小二乘法为理解 GWR提供了一个基础,在GWR中,一个观察值是通过与位置I的邻近来加权的,因此,一个观察值的加权在刻画过程中不再保持不变,而是随着I而变化。
就是说:β^(UI,VI)=(XTW(UI,VI)X)-1XTW(UI,VI)Y(5.3)
其中:
β^是β的估计值,N是空间样本数,K是自变量的个数,WIN是对位置I刻画模型时赋予数据点N的权重。
(二)空间权重函数的选择
首先,考虑全局模型(5.1)所隐含的权重函数:
其中,J代表空间中数据可观测到的特定点,I代表空间中参数被是1
另一个选择是排除一个给定距离D以外的所有点,这些点的权重被赋予0.如下所示:
这个权重函数是一个移动窗口,会造成不连续问题。福瑟林厄姆(FOTHERINGHAM,1996),查尔顿(CHARLTON,1997)用此函数作为权重函数进行了相关的研究。
一个克服权重不连续问题的方法是把WIJ指定为距离DIJ的连续且单调递减的函数,如采用如下的高斯函数:
其中,B是带宽,如果点I的数据被观测,则其他点的权重将根据高斯曲线随着距离DIJ的增加而减少。给定带宽B,距离DIJ越大,位置J所赋予的权重越小;另一方面,给定B,离点I足够远的点的权重将会趋于0.在现实中,与此对应的情形是如果两个地区的距离充分远,则一地区的区域经济发展对另一地区的经济发展没有影响。或者,采用如下的双重平方函数:
该函数特别有用,因为它提供了一个连续、近似高斯权重函数。
(三)校准(CALIBRATE)空间加权函数
从GWR的角度来看,参数估计部分地依赖于加权函数或核函数的选择,如在方程(5.5)中,如果D变得越大,则局部模型的解越趋于全局模型的解;如果D等于所研究空间任意两点间的最大距离,则两个模型将是相等的。又如在方程(5.6)中,若B趋于无穷大,任意两点的权重将趋于1,则被估计的参数变成一致的,此时GWR也等于OLS;反之,当带宽变得更小时,参数估计将愈加依赖于接近I的观测值。因此,问题的关键是如何选择一个适当的带宽或衰减函数。有许多标准适用于带宽的选择。
考虑方程(5.6)中B的选择,一个可能的方法是根据“最小平方”标准来选择B,使如下的数值最小:
其中,^YI(B)是使用带宽B来计算得到的YI的拟合值。为了得到YI的拟合值,需要在每个数据点估计βK(UI,VI),然后结合X-值计算YI的拟合值。然而正如FOTHERINGHAM ET AL。(1997A)提到的一样,在最小化方程(5.8)残差平方和的过程中会有个问题。假设选择很小的带宽,除了I点外所有其他点的权重都变得可以忽略,则样本点的拟合值将趋于其真实值,结果导致方程(5.8)的值为0,显然带宽B趋于0对分析是无益的。
为了解决这个问题,克利夫兰(CLEVELAND,1979),鲍曼(BOWMAN,1984)建议采用一种称作交叉确认(CROSS-VALIDATION,CV)的方法,此处CV如下:
其中^Y≠I(B)是YI的拟合值,在刻画过程中省略了点I的观测值。这样,当B变得很小时,模型仅仅刻画点I附近的样本而没有包括I本身。
为了取得最优的带宽,一个普遍采用的方法是使GWR模型的AIC最小(FOTHERINGHAM ET AL。,2002)。GWR模型的 AIC是根据HURVICH ET AL。(1998)的研究结果定义的。
其中,下标C表示“修正后的”AIC估计值,N是样本的大小,^σ是误差项估计的标准离差,TR(S)是GWR的S矩阵的迹,它是带宽的函数。S定义如下:
其中,Y和^Y是因变量和其估计值的向量。把不同的自由度考虑进两个模型后,这个AIC有利于评价GWR是否比OLS更好地模拟了数据。此外,式(5.10)还有更简单的形式:
(四)空间变异的显著性检验
根据上面讨论的方法得到GWR模型的参数估计值后,还有两个关键问题需要解决:一个是GWR模型是否比OLS模型更好地、显著地描述变量间的关系,另一个是每个参数估计集合是否在所研究的区域展示了空间变异。LEUNG ET AL。(2000 A),布伦斯坦等人(BRUNSDON ET AL。,1999)分别给出了解决上述两个问题的方法。在这里,简单介绍梁怡等人的工作:
解决第一个问题的方法之一是:
(1)给出零假设H0:GWR与OLS在描述变量间的关系上没有显著的差异;
(2)构造统计量F1:
GWR的残差平方和是:
其中,I是N阶单位矩阵,L记为:
L是N阶方阵,X TI是矩阵X的第I行,W I是W(I)。OLS的残差平方和是:
其中,Q=X(XTX)-1XT,δ1 =TR((I-L)T(I-L)),δ2 =TR[(I-L)T(I-L)]2
(3)检验假设。由于近似地有F1~F(δ21/δ2,N-K-1),给定一个显著水平 α,F1-α(δ21/δ2,N -K -1)表示上100(1 -α)分位数。如果F1 <F1-α(δ21/δ2,N-K-1),就拒绝零假设,由此推断GWR模型比OLS模型更好地、显著地描述变量间的关系;反之,就可以说和OLS模型相比,GWR模型不能显著提高拟合效果。
解决第一个问题的方法之二是FOTHERINGHAM ET AL。(2002)提出的利用式(5.10)和式(5.12)来判别。
解决第二个问题的方法是:
(1)给出零假设:
H0:βI(U1,V1)=βI(U2,V2)=…=βI(UN,VN),I=0,1,…,K
H1:不是所有的βI(UJ,VJ)(J=1,2,…,N)都相等。(2)构造统计量F3(I):
这里,F3的分布近似于F-分布。其中:
这里,β^(I)=β^I(U1,V1)=β^I(U2,V2)=…=β^I(UN,VN),J是所有
元素都为1的N阶方阵:
EI是第(I+1)个元素为1,其他元素都为0的列向量。
(3)检验假设。由于近似地有F3(I)~F(γ21/γ2,δ21/δ2),给定一个显著水平α,Fα(γ21/γ2,δ21/δ2)表示上100α分位数。如果F3(I)≥Fα(γ21/γ2,δ21/δ2),就拒绝零假设,否则,接受零假设。