局所空間統計量(LISA)でエリアごとに空間的自己相関があるかどうかを調べる

:::note note

この記事は、「Rではじめる地理空間データの統計解析入門」の第8章を読んで自分なりの理解を整理するために書いたものです。
間違った理解もあるかもしれませんので、正確な情報、詳しい説明は「[Rではじめる地理空間データの統計解析入門](https://www.kspub.co.jp/book/detail/5273036.html)」をお読みください。

:::




以前の記事では、空間相関の強さを評価する大域空間統計量(GISA)を調べる方法について学んだ。

https://qiita.com/nbayashi/items/c5a389afbecddc202bff


ここでは、局所空間統計量(`LISA: Local Indicator of Spatial Association`)ついて整理する。


# 局所空間統計量
各エリア間の局所的な空間相関を評価する統計量。
大域空間統計量が対象エリア全体で空間相関があるかどうかの評価であったのに対して、局所空間統計量は、エリアごとに空間相関があるかの評価が可能。

代表的な指標として、
	•	ローカルモランI統計量
	•	ローカルギアリC統計量
	•	ローカルG（ローカルG*）統計量
がある。



# ローカルモランI統計量
モランI統計量をエリアごとに分解した統計量らしい。
数値が正の時は自分と周囲が似た傾向を持つ正の空間相関がある、逆に負の値の時は自分と周囲が逆の傾向を持つ関係性が存在するということになる。
つまり以下のような判定

| 判定条件                  | 状態                   | 説明                            |
|:----------------------------|:-------------------------|:---------------------------------|
| `Ii > 0` ＆ `値が高い`        | High-High | 自分も周りも高い <br>（ホットスポット）                 |
| `Ii > 0` ＆ `値が低い`        | Low-Low | 自分も周りも低い  <br>（コールドスポット）              |
| `Ii < 0` ＆ `値が高い`        | High-Low    | 自分は高いが周りは低い<br>（一人勝ち）         |
| `Ii < 0` ＆ `値が低い`        | Low-High    | 自分は低いが周りは高い<br>（一人負け）          |

図で示すとこんなイメージだと思う。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61492/2b0b1b1a-2b4b-4dfd-8c68-bd2af315ae33.png)





---
## RでモランI統計量を求める

### データの読み込みと近接行列の用意

```R
library(sf);library(spdep); 

hokkaido <- read_sf("hokkaido.shp")
st_crs(hokkaido) <- 6668
hokkaido <- st_transform(hokkaido,crs=6680)

```
前回のGISAと同様北海道の市町村別人口のデータを使用する。
（データ[総務省統計局e-Stat](https://www.e-stat.go.jp/) 2020年国勢調査の境界データ、統計データを一部編集して使用）
自分でも試したいという人は、[書籍](https://www.kspub.co.jp/book/detail/5273036.html)の中ではNipponMapライブラリのサンプルデータを使って実践してますので、そちらをご参照ください。

SETAIというカラムに人口の値が入っている。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61492/c62f8900-d661-4af7-bbbb-7a7d3954f8a6.png)

```R
# 近接行列の作成
nb <- poly2nb(hokkaido,queen = TRUE)
w <- nb2listw(nb, zero.policy = TRUE)

```
近接行列は必要に応じて、最近隣kゾーン(knearneigh)や距離ベース（dnearneigh）で求めて良い。

### モランI統計量を求める

```R
# モランI統計量
lmoran <- localmoran(hokkaido$JINKO, list=w)
```

結果の最初5行分を出力すると以下のようになる。
```R
> lmoran[1:5,]
            Ii          E.Ii       Var.Ii       Z.Ii Pr(z != E(Ii))
1  0.013255149 -4.382847e-05 1.270557e-03  0.3730964      0.7090767
2 -0.036329373 -1.174023e-04 5.164100e-03 -0.5039122      0.6143231
3  0.000000000  0.000000e+00 0.000000e+00        NaN            NaN
4  0.022803548 -1.632128e-04 9.626757e-03  0.2340775      0.8149248
5 -0.001315811 -2.874308e-06 8.332763e-05 -0.1438299      0.8856348
```

- Ii：ローカルモランI統計量
- E.Ii：期待値。空間的にランダムだと仮定した時の理想的なローカルモランI統計量の平均。普通は0付近
- Var.Ii：分散：期待値に対するばらつきの大きさ。どれくらいIiがランダムにぶれるかの指標
- Z.Ii：Zスコア。E.Ii（期待されるランダム）と比べてIiがどれくらい離れているか
- Pr(z != E(Ii))：P値
がエリアごとに求められる。


## ローカルモランI統計量をプロットしてみる
```R
library(RColorBrewer)
hokkaido$lmoran <- lmoran[,"Ii"]
breaks <-c(-5,-2,-1,-0.5,0,0.5,1,2,5)
nc <-length(breaks)-1
pal <-rev(brewer.pal(n=nc,name="RdYlBu"))
plot(hokkaido[,"lmoran"],pal=pal,breaks=breaks)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61492/64baaa55-e261-467c-95b4-3a162d4437b6.png)


札幌付近ではローカルモランI統計量が高く、正の相関になっており、旭川や稚内、函館などいくつかの市町村では値がマイナスになっており、周囲は低いけど自分自身は高いという一人勝ちの状態になっていることがわかる。

## p値をプロットしてみる
```R
hokkaido$lmoran_p <- lmoran[,"Pr(z != E(Ii))"]
breaks <-c(0,0.01,0.05,0.10,1)
nc <-length(breaks)-1
pal <-rev(brewer.pal(n=nc,name="YlOrRd"))
plot(hokkaido[,"lmoran_p"],pal=pal,breaks=breaks)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61492/9ad05653-7556-4c85-9b1a-f32a8618cf35.png)

札幌市やその近郊の石狩市、小樽市、恵庭市などではp値が0.01未満で、有意な正の空間相関が見られる状態となった。
旭川市や稚内市、函館市などではp値>0.1となり、これらの市で見られた空間相関は統計的に有意差なしという結果となった。


# ローカルギアリC統計量
ローカルギアリC統計量は、グローバルはギアリC統計量と同様に、0以上の値を取り
- C<1の時に正の空間相関
- C=1の時にランダム
- C>1の時に負の空間相関を示す

その求め方は書籍には記載がなかった。どうやらローカルギアリC 統計量を直接求める関数はないらしい。
以前はelsaというパッケージでも求められたらしいが、このパッケージはすでにCRANからアーカイブ化されてる。一応手作業でも求めることはできるらしい。




# ローカルG統計量とローカルG*統計量

グローバルなG統計量と同じく、エリアごとのG統計量は人口や地価など非負の数値に対する統計量で、0以上の値を取り、0に近いほどランダム、数値が大きいほど大きな観測値が特定の地域に集中していることを指す。


G統計量、G*統計量の算出にはいずれも`localG`関数を用いる。
G*統計量では、G*統計量では自身のエリアを含んで近隣関係を考慮した統計量とするため、`nb2listw(include.self(nb1))`とする

## ローカルG/ローカルG*統計量を求める

```R
# ローカルG
w_b <- nb2listw(nb, zero.policy = TRUE)
lG <- localG(hokkaido$JINKO, list=w_b)

# ローカルG＊
w_b2 <- nb2listw(include.self(nb))
Gstar <- localG(hokkaido$JINKO, list=w_b2)
```

結果をプロット

```R

hokkaido$lG <- lG
breaks <-c(-5,-2.58,-1.96,-1.65,0,1.65,1.96,2.58,5)
nc <-length(breaks)-1
pal <-rev(brewer.pal(n=nc,name="RdYlBu"))
plot(hokkaido[,"lG"],pal=pal,breaks=breaks)

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61492/eb4ec56a-a7c2-4672-8c6e-6c59352f2107.png)

ここで、G統計量は0以上の値をとるはずなのに、図を見ると値がマイナスの地域もあるのがわかる。
localGで求められるのはG統計量そのものではなく、G統計量を標準化したZ値ということらしい。
つまり大きい値であるほど自分の値が高く、周囲も高い（ホットスポット）、小さい値であるほど自分の値が低く、周囲も低い（コールドスポット）ということになるらしい。



# まとめ
LISAによる統計量は色々あるけど、ローカルモランI統計量がわかりやすいきがすr

| 統計量                  | 値の幅     | 値の意味                                 |コード                                 |
|-------------------------|------------|------------------------------------------|-----------------------------------------|
| ローカルモランI統計量           | -1 ~ 1     | 正の値：正の空間相関<br>0：ランダム<br>負の値：負の空間相関 |`localmoran(data, list)`                                 |
| ローカルギアリC統計量           | 0以上の値  | 1未満：正の空間相関<br>1：ランダム<br>1以上：負の空間相関を示す |かつては`elsa`パッケージにあったが、今は使えない                                    |
| ローカルG統計量<br>（ローカルG*統計量）     | Zスコア  | Z値が大きいほど自分の値が高く、周囲も高い（ホットスポット）、Z値が小さいほど自分の値が低く、周囲も低い（コールドスポット）となる |`localG(data, list)`                                 |