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

演習の概要

• リスク分析：粗率、相対危険度、ポアソン確率、相対リスクのベイズ推定

注意事項

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

データ分析の準備

データのダウンロード

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

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

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

install.packages("spdep")
install.packages("sf")
install.packages("spatstat")
install.packages("tidyverse")
install.packages("kokudosuuchi")

library(spdep)
library(sf)
library(spatstat)
library(tidyverse)
library(kokudosuuchi)

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

needs::prioritize(magrittr)

データの読み込み

都道府県境界データ（ポリゴン）、都道府県庁所在地データ（ポイント）及び都道府県別属性データを読み込みます。

pref_pnt <- sf::st_read('pref_gov.shp')
jpn_pref <- sf::st_read('jpn_pref.shp')
hd06 <- readr::read_csv("hd06.csv")

都道府県境界ポリゴンデータに属性データhd06を結合します。

jpn_pref_hd06 <- dplyr::inner_join(jpn_pref, hd06, by=c("PREF_CODE"="PREF_CODE"))

リスク分析

粗率

確率地図を作成する

spdepのprobmap()で確率地図を作成する。粗率は以下のように計算する。

jpn_hd06_pm <- spdep::probmap(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100)
summary(jpn_hd06_pm)

      raw            expCount        relRisk            pmap       
 Min.   : 95.15   Min.   :  822   Min.   : 69.45   Min.   :0.0000  
 1st Qu.:137.58   1st Qu.: 1581   1st Qu.:100.42   1st Qu.:0.5122  
 Median :152.10   Median : 2381   Median :111.02   Median :1.0000  
 Mean   :151.09   Mean   : 3677   Mean   :110.28   Mean   :0.7424  
 3rd Qu.:169.73   3rd Qu.: 3729   3rd Qu.:123.89   3rd Qu.:1.0000  
 Max.   :199.93   Max.   :16995   Max.   :145.93   Max.   :1.0000  

主題図を作成する

まず、都道府県別心疾患死亡者数を地図上に可視化する

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=HD06)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Number of Death by Heart Disease")

粗率（ここでは、心疾患死亡者数/人口数x100）の計算結果はjpn_hd06_pmのrawに格納されている

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_pm$raw)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Raw Rate of Death by Heart Disease")

相対危険度

主題図を作成する

リスクの期待値と相対危険度の計算結果はjpn_hd06_pmのexpCountとrelRiskに格納されている。

期待値の主題図は以下のように描画できる。

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_pm$expCount)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Expected Count of Death by Heart Disease")

相対危険度の主題図は以下のように描画できる。

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_pm$relRisk)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Expected Count of Death by Heart Disease")

ポアソン確率

主題図を作成する

ポアソン確率の計算結果はjpn_hd06_pmのpmapに格納されている。

主題図は以下のようにして描画できる。

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_pm$pmap)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Poisson Probability of Death by Heart Disease")

相対リスクのベイズ推定

グローバルな経験ベイズ推定

Marshallのグローバルな経験ベイズ推定を行うには、spdepのEBest()関数を用いる。

jpn_hd06_ebg <- EBest(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100)
summary(jpn_hd06_ebg)

      raw             estmm       
 Min.   : 95.15   Min.   : 95.95  
 1st Qu.:137.58   1st Qu.:137.55  
 Median :152.10   Median :151.82  
 Mean   :151.09   Mean   :150.73  
 3rd Qu.:169.73   3rd Qu.:168.51  
 Max.   :199.93   Max.   :198.81  

以下のようにグローバルな経験ベイズ推定結果の主題図を作成する。

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_ebg$estmm)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Global EB Estimate of Death by Heart Disease Risk")

ローカルな経験ベイズ推定

Marshallのローカルな経験ベイズ推定を行うには、まず隣接行列を定義したのち、spdepのEBlocal()関数を用いる。

隣接行列を定義する

pref.tri.nb <- sf::st_coordinates(st_centroid(pref.pnt)) %>%
                    spdep::tri2nb()

ローカルな経験ベイズ推定を行う

jpn_hd06_ebl <- EBlocal(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100, pref.tri.nb)
summary(jpn_hd06_ebl)

      raw              est       
 Min.   : 95.15   Min.   : 95.7  
 1st Qu.:137.58   1st Qu.:138.0  
 Median :152.10   Median :151.3  
 Mean   :151.09   Mean   :150.8  
 3rd Qu.:169.73   3rd Qu.:169.0  
 Max.   :199.93   Max.   :199.1  

ローカルな経験ベイズ推定結果の主題図を作成する。

ggplot2::ggplot(data = jpn_pref_hd06, aes(fill=jpn_hd06_ebl$est)) + 
  geom_sf() + 
  scale_fill_distiller(palette="Reds", trans = "reverse") + 
  theme_bw() +
  labs(title = "Local EB Estimate of Death by Heart Disease Risk")

ローカルな経験ベイズ推定によるMoran’s Iの繰り返し検定

EBImoran.mc(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100, 
            nb2listw(pref.tri.nb, style="W", zero.policy=TRUE), 
            nsim=999, zero.policy=TRUE)

The default for subtract_mean_in_numerator set TRUE from February 2016

    Monte-Carlo simulation of Empirical Bayes Index (mean subtracted)

data:  cases: jpn_pref_hd06$HD06, risk population: jpn_pref_hd06$POPJ06/100
weights: nb2listw(pref.tri.nb, style = "W", zero.policy = TRUE)
number of simulations + 1: 1000

statistic = 0.015911, observed rank = 680, p-value = 0.32
alternative hypothesis: greater

通常のMoran’s Iの繰り返し検定

set.seed(1234)
moran.mc(jpn_pref_hd06$POPJ06/100, 
            listw=nb2listw(pref.tri.nb, style="W", zero.policy=TRUE), 
            nsim=999, zero.policy=TRUE)

    Monte-Carlo simulation of Moran I

data:  jpn_pref_hd06$POPJ06/100 
weights: nb2listw(pref.tri.nb, style = "W", zero.policy = TRUE)  
number of simulations + 1: 1000 

statistic = -0.1045, observed rank = 149, p-value = 0.851
alternative hypothesis: greater