環境ガバナンスのデータサイエンス/空間モデリング特論　R演習8

演習の概要

• 空間計量経済モデル

注意事項

• 慶應義塾大学SFCで開講している「環境ガバナンスのデータサイエンス」「空間モデリング特論」の授業履修者向けの演習用ページです。
• 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.5.3を使っています。
• Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

データ分析の準備

データのダウンロード

• 演習で用いるデータはここからダウンロードしてください。
• ここでは./直下にgisdataフォルダを作成し、setwd("gisdata")とディレクトリを指定しています

パッケージのインストール

• install.packages()でパッケージをインストールし、library()でパッケージを呼び出す
• ここでは以下のパッケージを用います

install.packages("spdep")
install.packages("spatstat")
install.packages("spgwr")
install.packages("nlme")
install.packages("lme4")

library(spdep)
library(spatstat)
library(spgwr)
library(nlme)
library(lme4)

以下の方法でパイプ %>% の優先順位を高める

needs::prioritize(magrittr)

データの作成

関東圏の市区町村境界データkanto_area.shpに、市区町村別住宅地平均地価データlph.csvを結合します。

kanto <- sf::st_read('kanto_area.shp')
lph <- readr::read_csv("lph.csv")
kanto_lph <- dplyr::inner_join(kanto, lph, by=c("JCODE"="JCODE"))

住宅地地価の空間分布を描画します。

ggplot2::ggplot(data = kanto, aes(fill=LPH)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues", trans = "reverse") + 
  theme_bw() +
  labs(title = "Land Price of House in Kanto Region")

空間計量経済分析

線形回帰モデル

まず、地価を被説明変数、夜間人口密度と第三次産業従業人口密度を説明変数とする重回帰モデルを推定しよう。 #### 重回帰モデルの推定

lph.lm <- lm(LPH~POPD+EMP3D,data=kanto)
summary(lph.lm)

Call:
lm(formula = LPH ~ POPD + EMP3D, data = kanto)

Residuals:
    Min      1Q  Median      3Q     Max 
-85.110  -1.873  -0.581   1.230  67.060 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.55112    0.56647   4.504 9.05e-06 ***
POPD         1.68162    0.10611  15.848  < 2e-16 ***
EMP3D        2.24666    0.07841  28.654  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.188 on 359 degrees of freedom
Multiple R-squared:  0.8289,    Adjusted R-squared:  0.8279 
F-statistic: 869.6 on 2 and 359 DF,  p-value: < 2.2e-16

誤差項の空間分布

kanto$lm.resid <- resid(lph.lm)
ggplot2::ggplot(data = kanto, aes(fill=lm.resid)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues") + 
  theme_bw() +
  labs(title = "Residuals of LPH model")

空間的自己相関

各変数に空間的自己相関があるかどうかを確認しよう。

空間隣接行列の定義

ここでは、ドロネー三角網を使って空間隣接行列を定義する

lph.tri.nb <- sf::st_coordinates(sf::st_centroid(kanto)) %>%
                    spdep::tri2nb()

Moran’s Iの算出

moran.test()関数を使ってMoran’s Iを計算します。

空間重み付け行列のスタイルをstyle=で指定する。ここでは、style="W"とした場合の例を示します。

住宅地地価のMoran’s I

spdep::moran.test(kanto$LPH,nb2listw(kanto.tri.nb, style="W"))

    Moran I test under randomisation

data:  kanto$LPH  
weights: nb2listw(kanto.tri.nb, style = "W")    

Moran I statistic standard deviate = 26.269, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.7455363769     -0.0027700831      0.0008114547 

夜間人口密度のMoran’sI

spdep::moran.test(kanto$POPD,nb2listw(kanto.tri.nb, style="W"))

    Moran I test under randomisation

data:  kanto$POPD  
weights: nb2listw(kanto.tri.nb, style = "W")    

Moran I statistic standard deviate = 28.329, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.8600771990     -0.0027700831      0.0009277053 

第三次産業従業人口密度のMoran’sI

spdep::moran.test(kanto$EMP3D,nb2listw(kanto.tri.nb, style="W"))

    Moran I test under randomisation

data:  kanto$EMP3D  
weights: nb2listw(kanto.tri.nb, style = "W")    

Moran I statistic standard deviate = 23.07, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.6230867560     -0.0027700831      0.0007359486 

一般化線形回帰モデル

空間データを用いたモデリングには、例えば面積によるオフセット項を考慮するモデル、緯度経度の影響を考慮するモデル、空間フィルタリングによる空間ラグを表現したモデル、空間ラグを考慮したモデルなどが考案されている。

オフセット項を考慮したモデル

面積によるオフセット項を考慮したモデルは、以下のように推定する。

lph.glm <- glm(LPH~POPD+EMP3D+offset(log(S)), data=lph)
summary(lph.glm)

Call:
glm(formula = LPH ~ POPD + EMP3D + offset(log(S)), data = lph)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-84.993   -2.106   -0.703    1.460   67.424  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.29070    0.56909  -2.268   0.0239 *  
POPD         1.75745    0.10660  16.487   <2e-16 ***
EMP3D        2.25733    0.07877  28.658   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 67.66322)

    Null deviance: 144897  on 361  degrees of freedom
Residual deviance:  24291  on 359  degrees of freedom
AIC: 2558

Number of Fisher Scoring iterations: 2

一般化加法モデル

mgcvパッケージのgam()関数を用いると、一般化加法モデルを以下のように推定できる。

lph.gam1 <- mgcv::gam(LPH~POPD+EMP3D+s(Easting, Northing), data=lph)
summary(lph.gam1)

Family: gaussian 
Link function: identity 

Formula:
LPH ~ POPD + EMP3D + s(Easting, Northing)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.56805    1.00229   8.548 4.55e-16 ***
POPD         0.21133    0.23183   0.912    0.363    
EMP3D        2.01558    0.08435  23.897  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F  p-value    
s(Easting,Northing) 25.27  28.18 2.213 0.000552 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.849   Deviance explained =   86%
GCV =   63.8  Scale est. = 58.818    n = 362

lph.gam2 <- mgcv::gam(LPH~POPD+EMP3D+offset(log(S))+s(Easting, Northing), data=lph)
summary(lph.gam2)

Family: gaussian 
Link function: identity 

Formula:
LPH ~ POPD + EMP3D + offset(log(S)) + s(Easting, Northing)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.93844    1.00134   4.932 1.29e-06 ***
POPD         0.23006    0.23165   0.993    0.321    
EMP3D        2.02890    0.08417  24.105  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F  p-value    
s(Easting,Northing) 25.39  28.23 2.445 9.01e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.85   Deviance explained = 86.6%
GCV = 63.383  Scale est. = 58.412    n = 362

AIC(lph.lm, lph.glm, lph.gam1, lph.gam2)

空間フィルタリングによる空間ラグモデル

Moran固有ベクトルによる空間ラグを表現したモデルは、以下のように推定できる。

lph.SF <- SpatialFiltering(LPH~POPD+EMP3D, data=lph, nb=lph.tri.nb, style="W")
lph.lm.SF <- lm(LPH~POPD+EMP3D+fitted(lph.SF), data=lph)
summary(lph.lm.SF)

Call:
lm(formula = LPH ~ POPD + EMP3D + fitted(lph.SF), data = lph)

Residuals:
    Min      1Q  Median      3Q     Max 
-85.101  -1.940  -0.434   1.265  66.787 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     2.55112    0.56457   4.519 8.46e-06 ***
POPD            1.68162    0.10575  15.902  < 2e-16 ***
EMP3D           2.24666    0.07814  28.751  < 2e-16 ***
fitted(lph.SF) 15.09787    8.16049   1.850   0.0651 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.16 on 358 degrees of freedom
Multiple R-squared:  0.8305,    Adjusted R-squared:  0.8291 
F-statistic: 584.8 on 3 and 358 DF,  p-value: < 2.2e-16

空間隣接行列による空間ラグを表現したモデル

空間隣接行列による空間ラグを表現したモデルは、以下のように推定できる。

lph.ME <- ME(LPH~POPD+EMP3D, data=lph, offset=log(S), 
             listw=nb2listw(lph.tri.nb, style="W"), alpha=0.5)
lph.glm.ME <- glm(LPH~POPD+EMP3D+offset(log(S))+fitted(lph.ME), data=lph)
summary(lph.glm.ME)

Call:
glm(formula = LPH ~ POPD + EMP3D + offset(log(S)) + fitted(lph.ME), 
    data = lph)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-85.184   -1.955   -0.375    1.368   65.987  

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -0.66441    0.59771  -1.112  0.26706    
POPD                   1.58652    0.11864  13.373  < 2e-16 ***
EMP3D                  2.27008    0.07829  28.994  < 2e-16 ***
fitted(lph.ME)vec1   -29.22642    9.26982  -3.153  0.00175 ** 
fitted(lph.ME)vec13   -1.34616    8.16682  -0.165  0.86917    
fitted(lph.ME)vec244  -1.06420    8.14866  -0.131  0.89617    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 66.37098)

    Null deviance: 144897  on 361  degrees of freedom
Residual deviance:  23628  on 356  degrees of freedom
AIC: 2553.9

Number of Fisher Scoring iterations: 2

自己回帰モデル

空間隣接行列に基づく空間重み付け行列を考慮した自己回帰モデルには、同時自己回帰モデルと条件付き自己回帰モデルがある。すでに空間隣接行列が定義されているものとして、これらのモデルを推定してみよう。

同時自己回帰モデル

同時自己回帰モデルはspdepパッケージのspautolm()関数を用いて、以下のように推定できる。family="SAR"と指定する。

lph.sar <- spautolm(LPH~POPD+EMP3D, data=kanto, 
nb2listw(lph.tri.nb, style="W"), family="SAR", method="eigen",  
verbose=TRUE)

Jacobian calculated using neighbourhood matrix eigenvalues
Computing eigenvalues ...

lambda: -0.8485824 function: -1309.331 Jacobian -20.23452 SSE 26261.24 
lambda: -0.1424868 function: -1274.356 Jacobian -0.6002126 SSE 24127.13 
lambda: 0.2939044 function: -1277.192 Jacobian -2.85387 SSE 24204.9 
lambda: 0.009451693 function: -1273.316 Jacobian -0.002721559 SSE 24068.19 
lambda: 0.006444693 function: -1273.315 Jacobian -0.001264476 SSE 24068.24 
lambda: 0.004091694 function: -1273.314 Jacobian -0.0005094297 SSE 24068.3 
lambda: 0.004051692 function: -1273.314 Jacobian -0.0004995132 SSE 24068.3 
lambda: 0.004052254 function: -1273.314 Jacobian -0.0004996521 SSE 24068.3 
lambda: 0.004052259 function: -1273.314 Jacobian -0.0004996532 SSE 24068.3 
lambda: 0.004052257 function: -1273.314 Jacobian -0.0004996526 SSE 24068.3 
lambda: 0.004052256 function: -1273.314 Jacobian -0.0004996524 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996522 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0 function: -1273.315 Jacobian 0 SSE 24068.47 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.00405228 function: -1273.314 Jacobian -0.0004996583 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052231 function: -1273.314 Jacobian -0.0004996462 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.00405228 function: -1273.314 Jacobian -0.0004996583 SSE 24068.3 
lambda: 0.00405228 function: -1273.314 Jacobian -0.0004996583 SSE 24068.3 
lambda: 0.00405228 function: -1273.314 Jacobian -0.0004996583 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 
lambda: 0.004052255 function: -1273.314 Jacobian -0.0004996523 SSE 24068.3 

summary(lph.sar)

Call: 
spautolm(formula = LPH ~ POPD + EMP3D, data = kanto, listw = nb2listw(lph.tri.nb, 
    style = "W"), family = "SAR", method = "eigen", verbose = TRUE)

Residuals:
      Min        1Q    Median        3Q       Max 
-85.12895  -1.87580  -0.58569   1.22201  67.15587 

Coefficients: 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 2.556939   0.566271  4.5154  6.32e-06
POPD        1.680409   0.105972 15.8570 < 2.2e-16
EMP3D       2.245907   0.078231 28.7088 < 2.2e-16

Lambda: 0.0040523 LR test value: 0.0015774 p-value: 0.96832 
Numerical Hessian standard error of lambda: 0.037053 

Log likelihood: -1273.314 
ML residual variance (sigma squared): 66.487, (sigma: 8.154)
Number of observations: 362 
Number of parameters estimated: 5 
AIC: 2556.6

条件付き自己回帰モデル

条件付き自己回帰モデルは、family="CAR"と指定することにより、以下のように推定できる。

lph.car <- spautolm(LPH~POPD+EMP3D, data=kanto,
nb2listw(lph.tri.nb, style="W"), family="CAR", method="eigen",  
verbose=TRUE)

Non-symmetric spatial weights in CAR model

Jacobian calculated using neighbourhood matrix eigenvalues
Computing eigenvalues ...

lambda: -0.8485824 function: -1282.453 Jacobian -20.23452 SSE 23938.56 
lambda: -0.1424868 function: -1273.593 Jacobian -0.6002126 SSE 24065.55 
lambda: 0.2939044 function: -1274.364 Jacobian -2.85387 SSE 24018.27 
lambda: 0.005206338 function: -1273.315 Jacobian -0.0008249936 SSE 24068.33 
lambda: 0.005921649 function: -1273.315 Jacobian -0.001067433 SSE 24068.31 
lambda: 0.007836143 function: -1273.314 Jacobian -0.001870019 SSE 24068.25 
lambda: 0.007813614 function: -1273.314 Jacobian -0.001859272 SSE 24068.25 
lambda: 0.007813422 function: -1273.314 Jacobian -0.001859181 SSE 24068.25 
lambda: 0.007813423 function: -1273.314 Jacobian -0.001859181 SSE 24068.25 
lambda: 0.007090829 function: -1273.315 Jacobian -0.001530956 SSE 24068.27 
lambda: 0.007537416 function: -1273.314 Jacobian -0.001730044 SSE 24068.26 
lambda: 0.007707997 function: -1273.314 Jacobian -0.001809305 SSE 24068.25 
lambda: 0.007773153 function: -1273.314 Jacobian -0.00184005 SSE 24068.25 
lambda: 0.007798041 function: -1273.314 Jacobian -0.001851862 SSE 24068.25 
lambda: 0.007807547 function: -1273.314 Jacobian -0.001856383 SSE 24068.25 
lambda: 0.007811178 function: -1273.314 Jacobian -0.001858112 SSE 24068.25 
lambda: 0.007812565 function: -1273.314 Jacobian -0.001858772 SSE 24068.25 
lambda: 0.007813095 function: -1273.314 Jacobian -0.001859025 SSE 24068.25 
lambda: 0.007813297 function: -1273.314 Jacobian -0.001859121 SSE 24068.25 
lambda: 0.007813361 function: -1273.314 Jacobian -0.001859151 SSE 24068.25 
lambda: 0.007813391 function: -1273.314 Jacobian -0.001859166 SSE 24068.25 
lambda: 0.00781338 function: -1273.314 Jacobian -0.00185916 SSE 24068.25 
lambda: 0.007813372 function: -1273.314 Jacobian -0.001859157 SSE 24068.25 
lambda: 0.007813384 function: -1273.314 Jacobian -0.001859162 SSE 24068.25 
lambda: 0.007813387 function: -1273.314 Jacobian -0.001859164 SSE 24068.25 
lambda: 0.007813382 function: -1273.314 Jacobian -0.001859162 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813382 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0 function: -1273.315 Jacobian 0 SSE 24068.47 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813429 function: -1273.314 Jacobian -0.001859184 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813334 function: -1273.314 Jacobian -0.001859139 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813429 function: -1273.314 Jacobian -0.001859184 SSE 24068.25 
lambda: 0.007813429 function: -1273.314 Jacobian -0.001859184 SSE 24068.25 
lambda: 0.007813429 function: -1273.314 Jacobian -0.001859184 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 
lambda: 0.007813381 function: -1273.314 Jacobian -0.001859161 SSE 24068.25 

summary(lph.car)

Call: 
spautolm(formula = LPH ~ POPD + EMP3D, data = kanto, listw = nb2listw(lph.tri.nb, 
    style = "W"), family = "CAR", method = "eigen", verbose = TRUE)

Residuals:
      Min        1Q    Median        3Q       Max 
-85.18887  -1.86743  -0.58564   1.21226  67.20937 

Coefficients: 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 2.556808   0.566197  4.5158 6.309e-06
POPD        1.680729   0.105962 15.8616 < 2.2e-16
EMP3D       2.245737   0.078225 28.7086 < 2.2e-16

Lambda: 0.0078134 LR test value: 0.0015214 p-value: 0.96889 
Numerical Hessian standard error of lambda: 0.05786 

Log likelihood: -1273.314 
ML residual variance (sigma squared): 66.487, (sigma: 8.1539)
Number of observations: 362 
Number of parameters estimated: 5 
AIC: 2556.6

空間的自己相関モデル

空間的自己回帰モデル（ 空間同時自己回帰モデル、spatial auto-regression model）は、被説明変数に空間的従属性を表現したモデルである。

空間従属性

被説明変数と、空間的重み付け（空間的従属性）を考慮した被説明変数との間に相関があるのかどうかをみてみよう。

LPH.lag <- lag.listw(nb2listw(lph.tri.nb, style="W"), kanto$LPH)
plot(LPH.lag, kanto$LPH, ylab="y", xlab="Wy", cex=1.5, lwd=2,
cex.axis=1.3, cex.lab=1.2)

非説明変数の経験累積分布は以下のようになる。

plot(ecdf(kanto$LPH), main="", cex=1.5, cex.axis=1.3, cex.lab=1.2)

空間重み付け行列を考慮した被説明変数の経験累積分布は、以下のようになる。

plot(ecdf(LPH.lag), main="", cex=1.5, cex.axis=1.3, cex.lab=1.2)

空間的自己回帰モデル（空間同時自己回帰ラグモデル）

空間的自己回帰モデルは、lagsarlm()関数を用いて推定する。

lph.lag <- lagsarlm(LPH~POPD+EMP3D, data=kanto,
nb2listw(lph.tri.nb, style="W"))
summary(lph.lag)

Call:
lagsarlm(formula = LPH ~ POPD + EMP3D, data = kanto, listw = nb2listw(lph.tri.nb, 
    style = "W"))

Residuals:
      Min        1Q    Median        3Q       Max 
-77.54050  -1.50614  -0.48205   1.32737  74.22407 

Type: lag 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  1.71625    0.53941  3.1817  0.001464
POPD         0.73909    0.14002  5.2783 1.304e-07
EMP3D        1.64289    0.10459 15.7079 < 2.2e-16

Rho: 0.42498, LR test value: 47.93, p-value: 4.4176e-12
Asymptotic standard error: 0.051562
    z-value: 8.2421, p-value: 2.2204e-16
Wald statistic: 67.931, p-value: 2.2204e-16

Log likelihood: -1249.35 for lag model
ML residual variance (sigma squared): 56.255, (sigma: 7.5003)
Number of observations: 362 
Number of parameters estimated: 5 
AIC: 2508.7, (AIC for lm: 2554.6)
LM test for residual autocorrelation
test value: 27.939, p-value: 1.252e-07

誤差項の空間的自己回帰モデル（空間同時自己回帰誤差モデル）

このモデルは、誤差項に空間的自己相関を明示的に用いたモデルである。

errorsarlm()関数を用いて、以下のように推定できる。

lph.err <- errorsarlm(LPH~POPD+EMP3D, data=kanto,
                      nb2listw(lph.tri.nb, style="W"))
summary(lph.err)

Call:
errorsarlm(formula = LPH ~ POPD + EMP3D, data = kanto, listw = nb2listw(lph.tri.nb, 
    style = "W"))

Residuals:
      Min        1Q    Median        3Q       Max 
-85.12895  -1.87580  -0.58569   1.22201  67.15587 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 2.556939   0.566271  4.5154  6.32e-06
POPD        1.680409   0.105972 15.8570 < 2.2e-16
EMP3D       2.245907   0.078231 28.7088 < 2.2e-16

Lambda: 0.0040523, LR test value: 0.0015774, p-value: 0.96832
Asymptotic standard error: 0.089641
    z-value: 0.045205, p-value: 0.96394
Wald statistic: 0.0020435, p-value: 0.96394

Log likelihood: -1273.314 for error model
ML residual variance (sigma squared): 66.487, (sigma: 8.154)
Number of observations: 362 
Number of parameters estimated: 5 
AIC: 2556.6, (AIC for lm: 2554.6)

Spatial Durbinモデル

空間ダービンモデルは、説明変数と被説明変数の両方に空間的従属性を取り入れたモデルである。

lagsarlm()関数を用いて、以下のように推定できる。

lph.durbin <- lagsarlm(LPH~POPD+EMP3D, data=kanto,
                       nb2listw(lph.tri.nb, style="W"), type="mixed")
summary(lph.durbin)

Call:
lagsarlm(formula = LPH ~ POPD + EMP3D, data = kanto, listw = nb2listw(lph.tri.nb, 
    style = "W"), type = "mixed")

Residuals:
      Min        1Q    Median        3Q       Max 
-64.43702  -1.76364  -0.36154   1.76424  62.27820 

Type: mixed 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  1.56274    0.54384  2.8735  0.004059
POPD        -0.75261    0.26640 -2.8251  0.004727
EMP3D        1.32107    0.12098 10.9193 < 2.2e-16
lag.POPD     2.42189    0.32409  7.4729 7.838e-14
lag.EMP3D    1.17284    0.25608  4.5799 4.652e-06

Rho: 0.052348, LR test value: 0.36311, p-value: 0.54679
Asymptotic standard error: 0.084488
    z-value: 0.61959, p-value: 0.53553
Wald statistic: 0.38389, p-value: 0.53553

Log likelihood: -1220.279 for mixed model
ML residual variance (sigma squared): 49.577, (sigma: 7.0411)
Number of observations: 362 
Number of parameters estimated: 7 
AIC: 2454.6, (AIC for lm: 2452.9)
LM test for residual autocorrelation
test value: 4.6991, p-value: 0.030178

空間的自己相関モデルを用いた予測

predict()関数を用いると、推定した各モデルの予測値を得ることができる。

predict(lph.err)
predict(lph.lag)
predict(lph.durbin)

空間従属性の検定

観測データや誤差項に空間的従属性を取り入れるかどうかについては、空間的従属性を取り入れなかった場合、すなわち通常の線形回帰モデルと比較して、空間的従属性を考慮することの統計的有意性を判断する方法が提案されている。これを空間従属性に関するラグランジュ乗数検定という。

lm.LMtests()関数を用いて、各モデルのラグランジュ乗数とp値を得ることができる。

lm.LMtests(lph.lm, nb2listw(lph.tri.nb), test=c("LMerr", "LMlag", "RLMerr",
"RLMlag", "SARMA"))

    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = LPH ~ POPD + EMP3D, data = kanto)
weights: nb2listw(lph.tri.nb)

LMerr = 0.0012212, df = 1, p-value = 0.9721


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = LPH ~ POPD + EMP3D, data = kanto)
weights: nb2listw(lph.tri.nb)

LMlag = 40.563, df = 1, p-value = 1.903e-10


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = LPH ~ POPD + EMP3D, data = kanto)
weights: nb2listw(lph.tri.nb)

RLMerr = 35.098, df = 1, p-value = 3.136e-09


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = LPH ~ POPD + EMP3D, data = kanto)
weights: nb2listw(lph.tri.nb)

RLMlag = 75.66, df = 1, p-value < 2.2e-16


    Lagrange multiplier diagnostics for spatial dependence

data:  
model: lm(formula = LPH ~ POPD + EMP3D, data = kanto)
weights: nb2listw(lph.tri.nb)

SARMA = 75.661, df = 2, p-value < 2.2e-16

地理的加重回帰モデル

空間的異質性と空間的従属性の両方を考慮した空間モデルとして、地理的加重回帰モデル（ geographically weighted regression model: GWR）がある。地域\(i\)ごとに異なる回帰係数\(\beta_i\)と空間重み付け関数\(W_i\)とで構成される。

空間重み付け関数\(W_i\)は、地域\(i\)と地域間距離地域\(d_i\)及びバンド幅数\(\theta\)を用いた距離低減関数を用いて表現される。

バンド幅\(\theta\)は、spgwrパッケージのgwr.sel()関数を用いて計算する。

coords <- matrix(0, nrow=length(kanto$LPH), ncol=2) 
coords[,1] <- kanto$Easting 
coords[,2] <- kanto$Northing 
lph.bw <- spgwr::gwr.sel(LPH~POPD+EMP3D, data=lph, coords=coords)

計算したバンド幅\(\theta\)を用いて、spgwrパッケージのgwr()関数により地理的加重回帰モデルを推定する。

lph.gwr <- spgwr::gwr(LPH~POPD+EMP3D, data=lph, coords=coords, 
               bandwidth=lph.bw,  hatmatrix=TRUE)
summary(lph.gwr$SDF)

Object of class SpatialPointsDataFrame
Coordinates:
              min       max
coord.x -120012.6  85196.34
coord.y -114091.4 115714.52
Is projected: NA 
proj4string : [NA]
Number of points: 362
Data attributes:
     sum.w         (Intercept)         POPD           EMP3D      
 Min.   : 14.85   Min.   :1.114   Min.   :1.378   Min.   :2.191  
 1st Qu.: 71.48   1st Qu.:2.265   1st Qu.:1.556   1st Qu.:2.250  
 Median :114.71   Median :3.110   Median :1.610   Median :2.264  
 Mean   :109.29   Mean   :3.158   Mean   :1.613   Mean   :2.260  
 3rd Qu.:150.33   3rd Qu.:4.006   3rd Qu.:1.654   3rd Qu.:2.274  
 Max.   :170.85   Max.   :5.659   Max.   :1.941   Max.   :2.302  
 (Intercept)_se      POPD_se          EMP3D_se           gwr.e         
 Min.   :0.7196   Min.   :0.1163   Min.   :0.07789   Min.   :-85.9098  
 1st Qu.:0.7505   1st Qu.:0.1186   1st Qu.:0.07817   1st Qu.: -1.5152  
 Median :0.8049   Median :0.1241   Median :0.07871   Median : -0.5827  
 Mean   :0.8794   Mean   :0.1508   Mean   :0.08415   Mean   : -0.1439  
 3rd Qu.:0.9506   3rd Qu.:0.1393   3rd Qu.:0.07992   3rd Qu.:  0.8215  
 Max.   :1.6788   Max.   :1.3407   Max.   :0.37839   Max.   : 65.3154  
      pred            pred.se          localR2       (Intercept)_se_EDF
 Min.   :  1.353   Min.   :0.5168   Min.   :0.7689   Min.   :0.7222    
 1st Qu.:  3.742   1st Qu.:0.6535   1st Qu.:0.8023   1st Qu.:0.7532    
 Median :  6.240   Median :0.7391   Median :0.8235   Median :0.8078    
 Mean   : 13.241   Mean   :0.8492   Mean   :0.8237   Mean   :0.8826    
 3rd Qu.: 18.461   3rd Qu.:0.9491   3rd Qu.:0.8467   3rd Qu.:0.9540    
 Max.   :164.160   Max.   :5.1645   Max.   :0.8778   Max.   :1.6849    
  POPD_se_EDF      EMP3D_se_EDF        pred.se      
 Min.   :0.1167   Min.   :0.07817   Min.   :0.5187  
 1st Qu.:0.1190   1st Qu.:0.07845   1st Qu.:0.6558  
 Median :0.1245   Median :0.07899   Median :0.7418  
 Mean   :0.1514   Mean   :0.08446   Mean   :0.8523  
 3rd Qu.:0.1398   3rd Qu.:0.08021   3rd Qu.:0.9526  
 Max.   :1.3455   Max.   :0.37976   Max.   :5.1832  

kanto$gwr.popd <- lph.gwr$SDF$POPD
ggplot2::ggplot(data = kanto, aes(fill=gwr.popd)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues") + 
  theme_bw() +
  labs(title = "Estiamted parameters of `POPD`")

kanto$gwr.emp3d <- lph.gwr$SDF$EMP3D
ggplot2::ggplot(data = kanto, aes(fill=gwr.emp3d)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues") + 
  theme_bw() +
  labs(title = "Estiamted parameters of `EMP3D`")

lph.gwr.resid <- resid(lph.gwr)
kanto$gwr.pred.se <- lph.gwr$SDF$pred.se
ggplot2::ggplot(data = kanto, aes(fill=gwr.pred.se)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues") + 
  theme_bw() +
  labs(title = "Predicted standard errors")

kanto$gwr.localR2 <- lph.gwr$SDF$localR2
ggplot2::ggplot(data = kanto, aes(fill=gwr.localR2)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Blues", trans = "reverse") + 
  theme_bw() +
  labs(title = "Local R^2")