【地理空間データ・R】同時自己回帰モデルの実装方法


# はじめに

地理空間データの勉強を始めようとして、村上大輔さんの<a href = "https://bookclub.kodansha.co.jp/title?code=1000042158">「Ｒではじめる地理空間データの統計解析入門」</a>(2022年4月6日発売)を読み初めました。**まだ半分しか読んでないですが、間違い無く良書**なので、本の中で紹介されていた手法を使い手を動かしてみました。なお、書いているコードの9割は本書に記載されている物を使います。前回は同時自己回帰モデル(Simultaneous Autoregressive; SAR)と呼ばれる空間相関をモデルを紹介すると同時に、そのモデルの取捨選択方法に触れました。今回はモデルを実装し、回帰係数の解釈までを行いたいと思います。

**前回**

https://qiita.com/Shu25/items/41ea4eafbf8109483034


# 空間ラグモデルの推定

*spaticalreg*パッケージを用い空間計量健在モデルを推定してみましょう。 空間ラグモデルは*lagsarlm*関数を使用します。
*listw*で空間重み行列を指定し、*tol.solve*で逆行列計算における特異性許容度合いを調整するパラメータを指定します。（この値を小さくする事でエラーが減りますが、計算時間が長くなります） 

```r
slm = spatialreg::lagsarlm(formula, data = boston.c, listw = w, tol.solve = 1.0e-20)
summary(slm)
```

![Screen Shot 2022-05-06 at 22.46.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/028b4c44-2514-235a-c86b-5bf4ff2a9fea.png)


空間ラグパラメーター$\rho$は0.54581、p~valueは2.22e-16 で、正の空間相関が証明されました。AICも2945.2と、lmモデルの3159と比べ、小さくなっています。

# 空間エラーモデル  

続いては空間エラーモデルです。*errorsarlm*関数を使用します。

```r
sem = spatialreg::errorsarlm(formula, data = boston.c, listw = w, tol.solve = 1.0e-20)
summary(sem)
```
![Screen Shot 2022-05-06 at 22.47.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/20b9d7b9-a23d-8b5d-493b-6804ba00ff15.png)

空間ラグパラメーター$\lambda$は0.66646、p~valueは2.22e-16 で、正の空間相関が証明されました。AICも2918.9と、lmモデルの3159や空間ラグモデル2945.2と比べ、小さくなっています。

# 空間ダービンモデル　

続いては空間ダービンモデルです。*lagsarlm*関数のパラメータ*type*を'mixed'にする事で使用できます。

```r
sdm = spatialreg::lagsarlm(formula, data = boston.c, listw = w, tol.solve = 1.0e-20, type = 'mixed')
summary(sdm)
```
![Screen Shot 2022-05-06 at 22.47.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/cfb97323-edb5-0557-f696-24c68d237b50.png)


パラメーター$\rho$は0.66147、p~valueは2.22e-16 で、正の空間相関が証明されました。AICも2909.7と、空間ラグモデル2945.2や空間エラーモデルの2918.9と比べ、小さくなっています。
係数のうち、*lag*がついているものが周辺への影響を表し、*LSTAT*と*RM*は周辺への影響があるものとされます。


### モデル毎の関数表　

|モデル|関数|パラメータ|
|---|:---:|:---:|
|空間ラグモデル(SLM)|lagsarlm|-|
|SLXモデル|lmSLX|-|
|空間エラーモデル(SEM)|errorsarlm|-|
|空間ダービンモデル(SDM)|lagsarlm|type = 'mixed'|
|空間ダービンエラーモデル(SDEM)|errorsarlm|etype = 'emixed'|
|SARARモデル|sacsarlm|-|
|Manski(GNM)モデル|sacsarlm|Durbin=TRUE| 


# Direct・Indirect Effects  

空間ダービンエラーモデルを例に直接・間接効果を評価していきます。

```r
sdem = spatialreg::errorsarlm(formula, data = boston.c, listw = w, tol.solve = 1.0e-20, etype = 'emixed')
summary(sdem)
```

![Screen Shot 2022-05-06 at 22.48.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/913c9a04-d826-3400-2035-2cbf628d1dbd.png)


### モンテカルロシミュレーションを使い直接・間接効果を評価 

```r
ires_sim = spatialreg::impacts(sdem, R = 1000)
summary(ires_sim)

# 対角成分の和を算出し、計算速度を早くする↓
# w2 = as(w, 'CasparseMatrix') 対角成分の和を算出し、計算速度を早くする
# trMat = spatial::trW(w2, type = 'mult')
# ires_sim = spatialreg::impacts(sdem, tr = reMat, R = 1000)
```
![Screen Shot 2022-05-06 at 22.48.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/a430e7b7-d29a-bcc6-626e-df8a3a8503d9.png)



*Impact measures* の欄に着目すると、  

- LSTAT（社会的弱者数）が増えると、同エリアの住宅価値がさがる  
- LSTAT（社会的弱者数）が増えると、周辺エリアの住宅価値が上がる  
- RM（部屋数）が増えると、同エリアの住宅価値が上がる  
- RM（部屋数）が増えると、周辺エリアの住宅価値が上がる  
- CRIM（犯罪数）が増えると、同エリアの住宅価値がさがる  
- CRIM（犯罪数）が増えると、周辺エリアの住宅価値がさがる  

このうち、5%の水準で統計的優位なものは、以下の4つになります。

- **LSTAT（社会的弱者数）が増えると、同エリアの住宅価値がさがる** 
- **RM（部屋数）が増えると、同エリアの住宅価値が上がる**
- **RM（部屋数）が増えると、周辺エリアの住宅価値が上がる**  
- **CRIM（犯罪数）が増えると、同エリアの住宅価値がさがる**  


# おわりに
村上大輔さんの<a href = "https://bookclub.kodansha.co.jp/title?code=1000042158">「Ｒではじめる地理空間データの統計解析入門」</a>で紹介されている、空間相関モデルの紹介とモデルの実装の仕方についてまとめました。また、モンテカルロシミュレーションを用い、直接・間接効果の評価も行いました。



### 参考文献

https://bookclub.kodansha.co.jp/title?code=1000042158