【地理空間データ】新型コロナ感染者の空間相関性について

# はじめに

地理空間データの勉強を始めようとして、村上大輔さんの<a href = "https://bookclub.kodansha.co.jp/title?code=1000042158">「Ｒではじめる地理空間データの統計解析入門」</a>(2022年4月6日発売)を読み初めました。まだ半分も読んでないのですが、間違い無く良書（英語でも名称が翻訳されていてさらに好き）なので、本の中で紹介されていた手法を使い手を動かしてみました。分析テーマは「新型コロナウィルスの感染者の割合は、近隣の都道府県に影響するのか」です。本書で紹介されているモランI統計量(Moran's I Statistics)を用い、分析しています。なお、扱うデータは違えど、書いているコードの9割は本書に記載されている物を使います。

:::note warn
警告
以下分析はあくまで試験的な物です。このデータ分析結果に対して一切の責任を負いかねます。
ネットの情報に惑わされず正確な情報を信じましょう。
:::


# モラン I 統計量(Moran's I Statistics) 
モラン I 統計量(Moran's I Statistics)とは空間相関（地理的に近いところと遠いところのデータの相関）を評価するものです。「エリアAの数値と全体の平均の差が、エリアAの近隣のエリアの差と近しいか」を比べています。以下の式で表されます。

$$
I = \dfrac{N}{\sum_i \sum_j w_{ij}}\dfrac{\sum_i \sum_jw_{ij}(y_j-\bar{y})(y_i-\bar{y})}{\sum_i(y_i-\bar{y})^2}
$$

$w_{i,j}$は近接行列を表し、ゾーン$i$と$j$の近さを表す重みです。

※詳しくは本書を読んでください

モランI統計量は相関係数のように「1に近づくほど正の相関が強い」ということになります。この場合の「正の相関」とは、空間相関が強い、つまり、周りのゾーンに数値が依存するという事です。

例）大阪府のたこ焼きの消費量が高ければ、大阪府近隣の県も高い


# 実装

### データ
データは以下のオープンソースと、*NipponMap*ライブラリの日本地図データを使用します。

https://www3.nhk.or.jp/news/special/coronavirus/data/

### コード

```r:setup  
library(rgdal)
library(spdep)
library(RColorBrewer)
library(NipponMap)
library(dplyr)
df = read.csv('nhk_news_covid19_prefectures_daily_data.csv')
df = df %>% 
  filter(日付 == "2022/4/20") %>% 
  filter(都道府県名 != "沖縄県")

shp = system.file('shapes/jpn.shp', package = 'NipponMap')[1] #日本地図読込
pref = read_sf(shp)　#日本地図読込
pref = pref[pref$name != "Okinawa",]　#可視化の都合上沖縄県を除外（ごめんなさい！！）
st_crs(pref) = 4326　#WGS84のEPSGコードを指定(地域のコードの形式を指定)
pref = st_transform(pref, crs = 6677)　#軽度緯度情報をメートルに変換
pref$jiscode = as.integer(pref$jiscode)
df = merge(df, pref, by.x = c('都道府県コード'), by.y = c('jiscode'))
```
>EPSGコードについては<a href = "https://lemulus.me/column/epsg-list-gis#WGS84">こちら</a>に詳しい解説が記載されています。

### 可視化
昨日(2022/4/20)の各地の直近1週間の人口10万人あたりの感染者数を表示しています。

```r: 可視化
Plot = data.frame(region = tolower(df$name), value = df$各地の直近1週間の人口10万人あたりの感染者数)
non_oki = c(Plot$region)
admin1_choropleth(country.name = "japan",
                  df           = Plot,
                  title        = paste("mean:", mean(Plot$value)),
                  legend       = "",
                  num_colors   = 4,
                  zoom = non_oki)
```


![Screen Shot 2022-04-21 at 22.18.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/54d3c940-6fc1-a93b-c2a0-fa53fe5915cc.png)
　
むむ、可視化だと何とも言えない...

### 検定
可視化だと何とも言えないので検定をしてみましょう。

```r
coords = st_coordinates(st_centroid(pref)) #都道府県の中心を算出
knn = knearneigh(coords, 4)　#近隣4ゾーンを指定
nb2 = knn2nb(knn)#nb形式に変換(近接情報を格納する型)
w = nb2listw(nb2) #近接情報をマトリクス形式に変換（距離に対しての相関係数行列的なもの）
```
> 近隣ゾーンについては<a href = "https://pro.arcgis.com/ja/pro-app/2.7/tool-reference/spatial-statistics/average-nearest-neighbor.htm">こちら</a>を参照してみてください。要はKNNです。

モランI統計量検定は*moran.test*で行えます。
*moran.test*はデフォルトで「H0(帰無仮説):$I \le E(I)$」「Ha(対立仮説):$I>E(I)$」が、設定されています。
```r:mran test
moran = moran.test(df$各地の直近1週間の人口10万人あたりの感染者数, list = w)
moran
```
### 結論
![Screen Shot 2022-04-21 at 22.59.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/2dd40ed1-f788-f3c6-4ef4-720cedcb5d86.png)

p-valueを見てみると0.05より小さい事がわかります。
よって、「H0(帰無仮説):$I \le E(I)$」は棄却されます(Reject Null) 
つまり、「Ha(対立仮説):$I>E(I)$を採用し、「5%の有意水準に基づき、ある県の新型コロナウィルスの感染者の割合（各地の直近1週間の人口10万人あたりの感染者数）は、近隣の都道府県に影響する」

要は、**「新型コロナウィルスの感染者の割合は、近隣の都道府県に影響する」**


# おわりに
村上大輔さんの<a href = "https://bookclub.kodansha.co.jp/title?code=1000042158">「Ｒではじめる地理空間データの統計解析入門」</a>で紹介されている、モラン I 統計量(Moran's I Statistics)を用い、「新型コロナウィルスの感染者の割合は、近隣の都道府県に影響するのか」という検定をしたところ、**「新型コロナウィルスの感染者の割合は、近隣の都道府県に影響する」** という分析結果がでました！（あくまでモランI統計量検定を用いた私個人の見解です。ネットの情報に惑わされず正確な情報を信じましょう）




### 参考文献

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

- https://lemulus.me/column/epsg-list-gis#WGS84
- https://www3.nhk.or.jp/news/special/coronavirus/data/