spatialreg | 空間滯後模型、空間誤差模型和空間杜賓模型簡單形式的R語言實現…

舊文重發，原文連結：https://mp.weixin.qq.com/s/pHRS9BfkCMe1uQOSkHmqAw

關於空間計量模型，小編是透過閱讀勒沙傑（James LeSage）和佩斯（R.Kelley Pace）合著的《空間計量經濟學導論》（Introduction of Spatial Econometrics）入門的，但是當時著重的是理解這些模型，並沒有用程式碼去實現。

作者也提供了書中模型的Matlab程式碼，可以在網站http://www.spatial-econometrics.com/（Econometrics Toolbox: by James P. LeSage）上檢視。另外，作者之一的佩斯還有一個網站http://www.spatial-statistics.com/（Spatial Statistics Software and Articles）,上面有一些空間統計領域的Matlab工具箱和論文。

由於小編目前也在學習階段，因此本篇只會使用R語言去實現這些模型的一些比較簡單的形式，使用的工具包是spatialreg，後續可能會分享一些更複雜的內容。

1 模型形式

1.1 空間滯後模型

自回歸模型（Autoregressive Model，AR）在時間序列分析中很易理解，即因變數與它的時間滯後值（Lag）存在相關性，這也意味著自回歸模型放棄了因變數獨立性的假設。

在空間計量模型中，空間滯後值被認為是鄰近空間單元的屬性（加權）值，因此下面是一個形式比較簡單的空間自回歸模型（最簡單的形式應該是不包含引數），也就是空間滯後模型（Spatial Lagged Model，SLM）：

  • 空間自回歸模型一般就是指空間滯後模型，但它也有一個更廣義的概念，即所有包含因變數的空間滯後項的模型；

  • 模型估計時會首先對引數進行估計，再使用廣義最小二乘法估計和其他引數。

1.2 空間誤差模型

空間誤差模型（Spatial Error Model，SEM）可以分解成如下兩步：

上述兩式合併得，

402 Payment Required

  • 第一個式子並非線性模型，因為不需要服從正態分布；

  • 模型估計時會首先對引數進行估計，再使用廣義最小二乘法估計和其他引數。

1.3 空間杜賓模型

空間杜賓模型（Spatial Durbin Model，SDM）假定因變數取值除受本地引數的影響外，還會受到鄰近地區的引數影響，即在模型中加入引數的空間滯後值：

402 Payment Required

1.4 複合模型

上述三個模型各自有針對性的假設，但這些假設相互之間並不排斥，可以在同一個模型中存在。

1.4.1 空間自相關模型

  • 該模型綜合了空間滯後模型和空間誤差模型，稱作空間自相關模型（Spatial Autocorrelation Model，SAC）；

  • 當或其中一個為0時，SAC就退化成了SEM或SLM。

1.4.2 空間杜賓（滯後）模型

402 Payment Required

  • 該模型綜合了空間滯後模型和空間杜賓模型，但習慣上仍稱作是空間杜賓模型。

1.4.3 空間杜賓誤差模型

  • 該模型綜合了空間誤差模型和空間杜賓模型，稱作空間杜賓誤差模型（Spatial Durbin Error Model，SDEM）。

1.4.4 空間自相關杜賓模型

  • 該模型是最一般的形式，同時綜合了三種基本模型，也可以認為是綜合了空間自相關模型和空間杜賓模型。

2 R語言程式碼

2.1 函式概述

上述7個模型形式可以透過spatialreg工具包中的3個函式來實現：

  • lagsarlm()：空間滯後模型、空間杜賓（滯後）模型

  • errorsarlm()：空間誤差模型、空間杜賓誤差模型

  • sacsarlm()：空間自相關模型、空間自相關杜賓模型

2.2 線性模型

範例資料如下：

1
2
3
4
5
library(sf)
library(tidyverse)
usa <- albersusa::counties_sf(proj = "laea") %>%
  mutate(fips = as.character(fips)) %>%
  left_join(socviz::county_data, by = c("fips" = "id"))

因變數取收入，引數取黑人比例。在執行空間計量模型前，先使用線性模型進行建模：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
data <- st_drop_geometry(usa)
model <- lm(hh_income ~ black, data = data)

summary(model)
## Call:
## lm(formula = hh_income ~ black, data = data)
##
## Residuals:
##    Min     1Q Median     3Q    Max
## -27680  -7341  -2171   4653  76006
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 47765.53     244.55  195.32   <2e-16 ***
## black        -199.18      14.29  -13.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11570 on 3141 degrees of freedom
## Multiple R-squared:  0.05825,    Adjusted R-squared:  0.05795
## F-statistic: 194.3 on 1 and 3141 DF,  p-value: < 2.2e-16

2.3 莫蘭指數

使用莫蘭指數檢驗因變數的空間自相關性：

1
2
3
4
5
6
7
8
9
10
library(spdep)
nb <- poly2nb(usa)
listW <- nb2listw(nb, zero.policy = T)

jpeg("1.jpeg", width = 8, height = 6, res = 600, units = "in")
moran.plot(usa$hh_income, listW,
           zero.policy = T,
           labels = F,
           pch = 20, cex = 0.1)
dev.off()

  • 部分縣沒有鄰接單元，設定zero.policy = T可以允許空間權重矩陣（實際空間資料結構是list）存在空元素。

a5fe18ff7ece3123f81adf06f75d29f7.png

2.4 空間計量模型

以lagsarlm()函式為例，它的完整語法結構如下：

1
2
3
4
5
6
7
lagsarlm(formula, data = list(),
         listw, na.action,
         Durbin, type,
         method="eigen", quiet=NULL,
         zero.policy=NULL, interval=NULL,
         tol.solve=.Machine$double.eps,
         trs=NULL, control=list())

本篇僅涉及以下幾個引數，其餘引數使用...代替：

1
2
3
4
lagsarlm(formula, data = list(),
         listw, Durbin,
         zero.policy = NULL,
         ...,)

  • formula：與對應的線性模型的運算式一致；

  • data：變數所在的資料框；

  • listw：空間權重矩陣；

  • Durbin：是否在模型中加入引數的空間滯後值；

  • zero.policy：針對權重矩陣存在空元素的應對方案，TRUE表示對應的空間滯後的權重為0。

2.4.1 空間滯後模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
library(spatialreg)
sl_model <- lagsarlm(hh_income ~ black, data = data,
                     listw = listW, zero.policy = T)

summary(sl_model)
## Call:lagsarlm(formula = hh_income ~ black, data = data, listw = listW,
##     zero.policy = T)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -26914.4  -4772.3  -1298.9   2996.4  72040.0
##
## Type: lag
## Regions with no neighbours included:
##  2788 2836 2995 3135 3140 3141 3143
## Coefficients: (numerical Hessian approximate standard errors)
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 14401.438    628.923  22.899 < 2.2e-16
## black        -121.955     10.236 -11.915 < 2.2e-16
##
## Rho: 0.71234, LR test value: 1783.5, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.012906
##     z-value: 55.195, p-value: < 2.22e-16
## Wald statistic: 3046.5, p-value: < 2.22e-16
##
## Log likelihood: -32973.68 for lag model
## ML residual variance (sigma squared): 67335000, (sigma: 8205.8)
## Number of observations: 3143
## Number of parameters estimated: 4
## AIC: 65955, (AIC for lm: 67737)

2.4.2 空間誤差模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
library(spatialreg)
se_model <- errorsarlm(hh_income ~ black, data = data,
                       listw = listW, zero.policy = T)

summary(se_model)
## Call:errorsarlm(formula = hh_income ~ black, data = data, listw = listW,
##     zero.policy = T)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -30053.7  -4590.2  -1184.3   2843.3  56012.9
##
## Type: error
## Regions with no neighbours included:
##  2788 2836 2995 3135 3140 3141 3143
## Coefficients: (asymptotic standard errors)
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 50696.822    725.463  69.882 < 2.2e-16
## black        -400.074     18.035 -22.183 < 2.2e-16
##
## Lambda: 0.82333, LR test value: 2384.3, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.011128
##     z-value: 73.99, p-value: < 2.22e-16
## Wald statistic: 5474.5, p-value: < 2.22e-16
##
## Log likelihood: -32673.29 for error model
## ML residual variance (sigma squared): 52510000, (sigma: 7246.4)
## Number of observations: 3143
## Number of parameters estimated: 4
## AIC: 65355, (AIC for lm: 67737)

2.4.3 空間杜賓誤差模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
library(spatialreg)
sd_model <- errorsarlm(hh_income ~ black, data = data,
                       listw = listW, zero.policy = T,
                       Durbin = T)

summary(sd_model)
## Call:errorsarlm(formula = hh_income ~ black, data = data, listw = listW,
##     Durbin = T, zero.policy = T)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -29657.6  -4639.4  -1123.4   2817.7  56324.1
##
## Type: error
## Regions with no neighbours included:
##  2788 2836 2995 3135 3140 3141 3143
## Coefficients: (asymptotic standard errors)
##              Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) 48577.976    753.720  64.4509 < 2.2e-16
## black        -392.654     17.888 -21.9512 < 2.2e-16
## lag.black     223.080     36.404   6.1278 8.909e-10
##
## Lambda: 0.81169, LR test value: 2384.6, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.011396
##     z-value: 71.228, p-value: < 2.22e-16
## Wald statistic: 5073.5, p-value: < 2.22e-16
##
## Log likelihood: -32655.17 for error model
## ML residual variance (sigma squared): 52280000, (sigma: 7230.5)
## Number of observations: 3143
## Number of parameters estimated: 5
## AIC: 65320, (AIC for lm: 67703)

6d5e6a242f4a3a4b060d618287e99d5c.jpeg