空間データと空間計量経済モデル

特殊講義「インターネットを活用した経済データの分析」講義資料

2/16/2017

1 Rにおける空間データと空間計量経済モデル

講義スライドで説明した通り,空間データは新たな情報を提供する一方で,解析上配慮を必要とする問題の増加にもつながる.空間計量経済学は,特に空間データの解析を用いた解析やデータの空間構造のモデル化を中心課題として発展してきた計量経済学の一分野である.ここでは,こうした空間データの解析に固有の問題を確認した上で,基礎的な空間計量経済 (spatial econometrics) モデルをRを用いて実行する.まずは,パッケージとサンプルデータを読み込む.

library(arm)
library(classInt)
library(maptools)
library(mgcv)
library(tidyverse)
library(spdep)
library(stargazer)
boston_tr <- readShapePoly(system.file("etc/shapes/boston_tracts.shp",
    package="spdep")[1], ID="poltract",
    proj4string=CRS(paste("+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66",
    "+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")))

2 データの確認と線形回帰の推定

ここでは,空間計量経済/統計モデルを扱う定番のパッケージであるspdepパッケージのサンプルデータbostonを用いて,空間データと空間計量経済モデルの初歩的な取り扱いを学ぶ. まずは,上で読み込んだboston_trオブジェクトからデータフレームを取り出す.

class(boston_tr)    # check obj class
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
boston_dat = boston_tr@data %>% # pull out data.frame obj
    mutate_at(  # convert some variables to log scale
        vars(DIS, RAD, LSTAT),
        funs(log)   
        )

1行目の出力からわかる通り,boston_trはSpatialPolygonsDataFrameクラスのオブジェクトである.SpatialPolygonsDataFrameクラスはRにおける空間データクラスの1つであり,行列やデータフレームオブジェクトよりも複雑な構造をしている.おおまかにいえば,次のコードで描ける「地図 (ここでは行政区画のポリゴン)」の部分と,それぞれの行政区画に対応する「データ」を保持するデータフレームを主要な構成要素とするオブジェクトといえる.その他,空間データの処理のための細々とした情報も保持しているが,ここでは省略する.SpatialPolygonsDataFrameクラスのオブジェクトを引数にしてplot()関数を使えば,簡単に地図が描ける.

par(mar=c(0,0,1,0),oma=c(0,0,0,0))
plot(boston_tr, lwd = .2, col = "gray90", border="gray20")

さて,boston_datには地図の行政区画1つ1つに対応するデータフレームが格納されている.このデータフレームを用いて,線形回帰を実行しよう.

## Standard linear regression
mod_0 = as.formula(MEDV ~ CRIM + ZN + INDUS + CHAS + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT + NOX + I(NOX^2) + RM + I(RM^2))
res_0 = lm(mod_0, data=boston_dat)
stargazer(res_0, single.row = TRUE, type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                MEDV            
## -----------------------------------------------
## CRIM                     -0.156*** (0.026)     
## ZN                         0.018 (0.011)       
## INDUS                      0.001 (0.050)       
## CHAS1                    2.253*** (0.703)      
## AGE                       -0.002 (0.012)       
## DIS                      -5.198*** (0.776)     
## RAD                      1.916*** (0.408)      
## TAX                      -0.008*** (0.003)     
## PTRATIO                  -0.686*** (0.107)     
## B                        0.006*** (0.002)      
## LSTAT                    -7.970*** (0.529)     
## NOX                       11.739 (21.943)      
## I(NOX2)                  -23.788 (16.302)      
## RM                      -19.428*** (2.497)     
## I(RM2)                   1.757*** (0.196)      
## Constant                111.123*** (10.226)    
## -----------------------------------------------
## Observations                    506            
## R2                             0.831           
## Adjusted R2                    0.826           
## Residual Std. Error      3.839 (df = 490)      
## F Statistic          160.572*** (df = 15; 490) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

解析結果が面白いかどうかはともかくとして,ここではこの回帰分析の結果の残差に何らかの空間構造・空間的な偏りがあるか否かを検討しよう.このような

  1. 非空間的な回帰モデルとのその残差の空間的分布を検討し,
  2. 必要であれば空間構造を明示的に取り込んだ回帰モデル (e.g., 空間計量経済モデル) を用いて解析する,

という手順は,空間データを解析する上での基本的な手順である. まずは,デフォルトのresiduals()関数を用いて残差を抽出しよう.

## Pull out and store residuals
boston_dat$Rsd_BM = residuals(res_0)
boston_tr@data = boston_dat

上記で得た残差は,空間的にどのように分布しているのだろうか.この点を確かめるため,残差を地図上にプロットする.

nn = 5
br.palette = colorRampPalette(c("red3", "gray85","royalblue"), space = "rgb")
classes_km = classIntervals(boston_tr$Rsd_BM, n=nn, style = "fisher", rtimes = 1)
pal = br.palette(nn)
cols = findColours(classes_km, pal) ## Set cols
par(mar = c(0,0,1,0), oma = c(0,0,0,0))
plot(boston_tr, col=cols, border="black", lwd=.2)
legend(
    "topleft", inset=c(0,0), title="OLS residuals", bty="n", 
    cex=.85, fill=attr(cols,"palette"), 
    legend=names(attr(cols,"table")), 
    ncol=1)

地図を眺める限り,残差が空間的に偏っているように見える. 残差に有意な空間的偏りがあるか否かを厳密に検討してみよう. 空間計量経済モデルでは,観察/分析単位 (ここでは行政区画) 間の空間的配置や近接関係を記述するために空間重み行列 (spatial weight matrix, SWM)と呼ばれる行列を用いることが多い.SWMは\(N\times N\)の行列\(\boldsymbol{W}\) (\(N\)は分析単位の数) を用いて記述され,分析単位\(i\)\(j\)が相互に影響を与え合う (e.g., 「隣接」している) 場合\(ij\)要素\(w_{ij}\)\(>0\)に,与えない場合は\(0\)となる.ここでは,もっとも単純なSWMの定義の1つである隣接行列 (contiguity based SWM; \(i\)\(j\)が隣接している場合は\(1\), していなければ\(0\)) を用いる.

(ポリゴンデータから) Rで隣接行列型のSWMを生成するには,spdep::poly2nb()関数とspdep::nb2listw()関数を用いる.なお,ここでは\(\boldsymbol{W}\)は行和\(w_{i+}=1\)に行基準化してある (空間自己相関指標の理解がりやすくなるため).

## Generate SWM
poly_cods = coordinates(boston_tr) # 行政区画ポリゴンの幾何中心点 (作図用)
poly_cont = poly2nb(boston_tr, queen = TRUE) 
poly_SWM = nb2listw(poly_cont, style="W", zero.policy=TRUE)

以上で,SWMを定義できた.SWMは,次の地図のようにある種のネットワークとして定義される (点は各行政区画の幾何中心点,実線は隣接関係を示す).

par(mar=c(0,0,1,0),oma=c(0,0,0,0))
plot(boston_tr, lwd = .2, col = "gray90", border="gray60")
plot(poly_SWM, coords = poly_cods, cex=0.35, lwd=.5, pch=19, col = "royalblue", add=TRUE)

空間自己相関 (spatial autocorrelation) の統計的指標はいくつか提案されているが,古典的かつ頻繁に用いられる指標に (Global) Moran’s Iがある.この指標や下で用いる空間計量経済 (回帰) モデルの詳細は,

を参照のこと.

詳細は省略するが,Moran’s Iは「相関係数の空間番」とざっくりと理解しておいてここでは差し支えない.SWMが\(\boldsymbol{W}\)は行和\(w_{i+}=1\)に行基準化されているとき,Moran’s Iは-1から1の値をとり,ピアソンの相関係数と同様に解釈できる.

線形回帰の残差についてMoran’s Iの検定を行なうには,spdep::lm.morantest()関数を用いる.

## Moran's I
lm.morantest(model=res_0, listw=poly_SWM)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = mod_0, data = boston_dat)
## weights: poly_SWM
## 
## Moran I statistic standard deviate = 10.449, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.2553506767    -0.0181115321     0.0006849124

出力から分かる通り,上記の線形回帰の残差には有意な (正の) 空間的自己相関があるという結果が得られた.

3 データの空間構造に配慮した回帰モデル

空間構造を明示的にモデル化する回帰モデルは,一般に空間計量経済モデルと呼ばれる.近年,様々なモデルが提案され,同一のモデルの呼称さえも共有されていない節があるが,ここではもっとも基本的なモデルである空間ラグモデル (spatial lag model, SLM) と (空間自己回帰的な誤差構造をもつ) 空間誤差モデル (spatial autoregressive error model, SEM) を実行してみよう.

SLMは従属変数の空間自己相関構造を導入した回帰モデルで,次のように定義される \[ y_i = \rho\sum^{n}_{j=1}w_{ij}y_{j} + \beta_{1} X_{i1} + \cdots + \beta_{k} X_{ik} + \epsilon_{i},\,\,\,\mathrm{for}\, i = 1,\ldots, n. \] ただし,\(\epsilon_{i} \sim N(0,\sigma^2)\)は誤差項,また\(X_{i1} = 1\)であり,\(\rho\)は従属変数の空間自己相関を捉える空間パラメータである.

他方,SEMは誤差項の空間自己相関構造を導入した回帰モデルで,次のように定義される \[ y_i = \beta_{1} X_{i1} + \cdots + \beta_{k} X_{ik} + \mu_{i},\\ \mu_{i} = \lambda \sum^{n}_{j=1}w_{ij}\mu_{j} + \epsilon_{i}. \]

ここで,\(\epsilon_{i} \sim N(0,\sigma^2)\), また\(\lambda\)は誤差項の空間自己相関を捉える空間パラメータである.\(\mu_{i}\)の定式化から分かる通り,SEMでは誤差項の空間自己相関構造がポイントになる.ただし,誤差こうの空間構造を定式化する手法は様々なものが提案されているため,同じ「空間誤差モデル (SEM)」という言葉を使っていても,論文・書籍等によって異なるモデルを意味する場合があることに注意してほしい.

いずれのモデルも,spdepパッケージの関数で実行できる.

mat_meth = "Chebyshev"
## SLM
res_SLM_ml = lagsarlm(mod_0, listw = poly_SWM, method=mat_meth, data=boston_tr@data, tol.solve = 1e-15)
## SEM
res_SAEM_ml = errorsarlm(mod_0, listw = poly_SWM, method=mat_meth, data=boston_tr@data)

これで,空間構造を除く回帰式は上述の線形回帰と等しいSLMとSEMの推定ができた. 結果を確認してみる.

summary(res_SLM_ml)
## 
## Call:lagsarlm(formula = mod_0, data = boston_tr@data, listw = poly_SWM, 
##     method = mat_meth, tol.solve = 1e-15)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -24.08709  -1.71892  -0.14888   1.51619  23.84639 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept)  79.5566503   9.6898754   8.2103 2.220e-16
## CRIM         -0.1139147   0.0240612  -4.7344 2.197e-06
## ZN            0.0193806   0.0101769   1.9044 0.0568611
## INDUS         0.0162607   0.0450196   0.3612 0.7179557
## CHAS1         1.3138230   0.6400697   2.0526 0.0401090
## AGE          -0.0039243   0.0104993  -0.3738 0.7085804
## DIS          -4.3793973   0.7078814  -6.1866 6.147e-10
## RAD           1.4392839   0.3714392   3.8749 0.0001067
## TAX          -0.0075470   0.0023345  -3.2328 0.0012257
## PTRATIO      -0.3661131   0.1020828  -3.5864 0.0003352
## B             0.0043609   0.0019701   2.2135 0.0268656
## LSTAT        -6.2752092   0.5085571 -12.3392 < 2.2e-16
## NOX          16.7528676  19.7953845   0.8463 0.3973845
## I(NOX^2)    -22.2944600  14.6900581  -1.5177 0.1291011
## RM          -15.8687801   2.2631501  -7.0118 2.353e-12
## I(RM^2)       1.4811586   0.1781206   8.3155 < 2.2e-16
## 
## Rho: 0.31651, LR test value: 79.454, p-value: < 2.22e-16
## Asymptotic standard error: 0.034314
##     z-value: 9.2239, p-value: < 2.22e-16
## Wald statistic: 85.08, p-value: < 2.22e-16
## 
## Log likelihood: -1350.787 for lag model
## ML residual variance (sigma squared): 11.964, (sigma: 3.4589)
## Number of observations: 506 
## Number of parameters estimated: 18 
## AIC: 2737.6, (AIC for lm: 2815)
## LM test for residual autocorrelation
## test value: 17.435, p-value: 2.9723e-05
summary(res_SAEM_ml)
## 
## Call:errorsarlm(formula = mod_0, data = boston_tr@data, listw = poly_SWM, 
##     method = mat_meth)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -23.90525  -1.79213  -0.25984   1.58774  23.70230 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept)  94.1620818  12.3996281   7.5939 3.109e-14
## CRIM         -0.1251484   0.0254162  -4.9240 8.481e-07
## ZN            0.0138644   0.0122983   1.1273 0.2595957
## INDUS         0.0046164   0.0606767   0.0761 0.9393542
## CHAS1         0.7965288   0.8072466   0.9867 0.3237785
## AGE          -0.0132518   0.0121578  -1.0900 0.2757201
## DIS          -5.5789118   1.0673388  -5.2269 1.723e-07
## RAD           1.6004820   0.4707197   3.4001 0.0006737
## TAX          -0.0087506   0.0028411  -3.0800 0.0020700
## PTRATIO      -0.5107428   0.1276614  -4.0008 6.314e-05
## B             0.0068842   0.0026766   2.5720 0.0101112
## LSTAT        -6.8391907   0.5511540 -12.4089 < 2.2e-16
## NOX          -9.8247343  29.0105710  -0.3387 0.7348655
## I(NOX^2)     -5.8411233  21.5827772  -0.2706 0.7866693
## RM          -14.0293287   2.4483858  -5.7300 1.004e-08
## I(RM^2)       1.3817962   0.1892958   7.2997 2.884e-13
## 
## Lambda: 0.56225, LR test value: 82.569, p-value: < 2.22e-16
## Asymptotic standard error: 0.049103
##     z-value: 11.45, p-value: < 2.22e-16
## Wald statistic: 131.11, p-value: < 2.22e-16
## 
## Log likelihood: -1349.23 for error model
## ML residual variance (sigma squared): 11.321, (sigma: 3.3647)
## Number of observations: 506 
## Number of parameters estimated: 18 
## AIC: 2734.5, (AIC for lm: 2815)

いずれの結果でも,lm()関数を用いた結果の出力にはなかったrhoあるいはlambdaの情報が出力されている.AICLog likelihoodという出力があるのは,これらの空間回帰モデルがOLSではなく最尤法 (maximum likelihood estimation) という推定方法 (本講義では扱わない) を用いている (OLSでは適切な推定ができない) ためであり,空間回帰モデルだからという訳ではない.rholambdaはそれぞれ,上述のSLMとSEMの回帰式の\(\rho\)\(\lambda\)に対応する.いずれのモデルでも,有意な正の空間自己相関が確認できる.

では,空間構造に配慮した回帰モデルから得た回帰係数と,非空間的な線形回帰から得た回帰係数の推定値には差異があるのだろうか.こうした比較には,rholambdaは無視されてしまうがstargazer()関数を用いると手っ取り早い.

stargazer(res_0, res_SLM_ml, res_SAEM_ml, single.row = TRUE, type = "text")
## 
## ===================================================================================
##                                           Dependent variable:                      
##                     ---------------------------------------------------------------
##                                                  MEDV                              
##                                OLS                 spatial            spatial      
##                                                 autoregressive         error       
##                                (1)                   (2)                (3)        
## -----------------------------------------------------------------------------------
## CRIM                    -0.156*** (0.026)     -0.114*** (0.024)  -0.125*** (0.025) 
## ZN                        0.018 (0.011)         0.019* (0.010)     0.014 (0.012)   
## INDUS                     0.001 (0.050)         0.016 (0.045)      0.005 (0.061)   
## CHAS1                   2.253*** (0.703)       1.314** (0.640)     0.797 (0.807)   
## AGE                      -0.002 (0.012)         -0.004 (0.010)     -0.013 (0.012)  
## DIS                     -5.198*** (0.776)     -4.379*** (0.708)  -5.579*** (1.067) 
## RAD                     1.916*** (0.408)       1.439*** (0.371)   1.600*** (0.471) 
## TAX                     -0.008*** (0.003)     -0.008*** (0.002)  -0.009*** (0.003) 
## PTRATIO                 -0.686*** (0.107)     -0.366*** (0.102)  -0.511*** (0.128) 
## B                       0.006*** (0.002)       0.004** (0.002)    0.007** (0.003)  
## LSTAT                   -7.970*** (0.529)     -6.275*** (0.509)  -6.839*** (0.551) 
## NOX                      11.739 (21.943)       16.753 (19.795)    -9.825 (29.011)  
## I(NOX2)                 -23.788 (16.302)       -22.294 (14.690)   -5.841 (21.583)  
## RM                     -19.428*** (2.497)     -15.869*** (2.263) -14.029*** (2.448)
## I(RM2)                  1.757*** (0.196)       1.481*** (0.178)   1.382*** (0.189) 
## Constant               111.123*** (10.226)    79.557*** (9.690)  94.162*** (12.400)
## -----------------------------------------------------------------------------------
## Observations                   506                   506                506        
## R2                            0.831                                                
## Adjusted R2                   0.826                                                
## Log Likelihood                                    -1,350.787         -1,349.230    
## sigma2                                              11.964             11.321      
## Akaike Inf. Crit.                                 2,737.574          2,734.460     
## Residual Std. Error     3.839 (df = 490)                                           
## F Statistic         160.572*** (df = 15; 490)                                      
## Wald Test (df = 1)                                85.080***          131.111***    
## LR Test (df = 1)                                  79.454***          82.569***     
## ===================================================================================
## Note:                                                   *p<0.1; **p<0.05; ***p<0.01

CHAS変数の回帰係数の推定値に,分かりやすいモデル間の違いが見られる.CHAS変数は非空間的な線形回帰では正の回帰係数をもち有意になっているものの,SLM (真ん中の列) では回帰係数が小さくなっている.SEMでは回帰係数がさらに小さくなり,有意ではなくなっている.

いずれのモデルがデータを最もよく説明しているかという疑問に答えるには別の作業が必要になるが,(非空間的な) 線形回帰から得られたCHAS変数についての結果は,データの空間構造を無視したことに影響されている可能性があることを示唆する. ただし,非空間的な回帰モデルの回帰係数とSLMの回帰係数は同様に解釈することができず,直接比較できない (空間構造を明示的にモデル化することに起因して,限界効果が異なる) ことを付言しておく (堤・瀬谷 2014: 第6章).