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

演習の概要

• 空間集積性の検定、ピアソンの\(\chi^2\)検定、Potthoff & Whittinghillの検定、Stone’s test、Tango’s test、Whittermoreの統計量、Besag-Newell’s test、OpenshawのGAM、Kulldorffの空間スキャン統計量

注意事項

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

データ分析の準備

データのダウンロード

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

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

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

install.packages("DCluster")
install.packages("spdep")

library(DCluster)
library(spdep)

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

needs::prioritize(magrittr)

データの読み込み

都道府県別出生率データを読み込み、出生率の期待値を計算します。

data76 <- read.table("data76.csv", sep=",", header=TRUE, row.names=1)
data76_OE <- data.frame(Observed=data76$n.birth)
data76_OE <- cbind(data76_OE,  
Expected=data76$pop*sum(data76$n.birth)/sum(data76$pop), 
x=data76$Easting, y=data76$Northing)

空間集積性の検定

ピアソンの\(\chi^2\)検定

DClusterパッケージのachisq.stat()関数で\(\chi^2\)統計量を計算し、achisq.test()関数で\(\chi^2\)検定の結果を示す。これにより、ランダム性に関する検定を行う。

achisq.stat(data76_OE, lambda=1)

$T
[1] 1616.169

$df
[1] 79

$pvalue
[1] 3.154174e-285

achisq.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100)

Chi-square test for overdispersion 

    Type of boots.: parametric 
    Model used when sampling: Poisson 
    Number of simulations: 100 
    Statistic:  1616.169 
    p-value :  0.00990099 

さらに、boot()関数を用いてブートストラップ法により標本を生成し、\(\chi^2\)検定を行うこともできる。

data76_achb_pb <- boot(data76_OE, statistic=achisq.pboot, sim="parametric", ran.gen=poisson.sim, R=100)
plot(data76_achb_pb)

Potthoff & Whittinghillの検定（過分散の有無に関する検定）

すべての地区において相対リスクが互いに類似しているかどうかを検証する方法です。

pottwhitt.stat(data76_OE)

$T
[1] 3840008393

$asintmean
[1] 3745501200

$asintvat
[1] 591789189600

$pvalue
[1] 0

pottwhitt.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100)

Potthoff-Whittinghill's test of overdispersion 

    Type of boots.: parametric 
    Model used when sampling: Poisson 
    Number of simulations: 100 
    Statistic:  3840008393 
    p-value :  0.00990099 

boot()関数を用いてブートストラップ法により標本を生成し、過分散の有無に関する検定を行う。

data76_pw_pb <- boot(data76_OE, statistic=pottwhitt.pboot, sim="parametric", ran.gen=poisson.sim, R=100)
plot(data76_pw_pb)

Stone’s test

各地区の相対リスクが等しいという帰無仮説\(H_0\)に対して、特定の地域からの距離が増加する毎に相対リスク\(\theta\)が減少するという対立仮説\(H_1\)を立て、仮説検定する方法です。

stone.stat(data76_OE, region=which(row.names(data76_OE)=="20"), lambda=1)

             region 
 1.000948 79.000000 

stone.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, region=which(row.names(data76_OE)=="20"), lambda=1)

Stone's Test for raised incidence around locations 

    Type of boots.: parametric 
    Model used when sampling: Poisson 
    Number of simulations: 100 
    Statistic:  1.000948 
    p-value :  0.8910891 

boot()関数を用いてブートストラップ法により標本を生成し、検定を行う。

data76_st_pb <- boot(data76_OE, statistic=stone.pboot, sim="parametric", ran.gen=poisson.sim, R=100, region=which(row.names(data76_OE)=="20"))
plot(data76_st_pb)

Tango’s test

空間オブジェクト同士の隣接性を考慮して空間属性の集積性を仮説検定する方法。

spdepパッケージを使って空間隣接行列と空間重み付け行列を計算する。さらにDClusterパッケージのtango.stat()関数により、Tango統計量\(T\)を用いた集積性の有無を\(\chi^2\)検定する。

data76_OE <- cbind(data76_OE, x=data76$Easting, y=data76$Northing)
coords <- as.matrix(data76_OE[,c("x","y")])
dlist <- dnearneigh(coords, 0, Inf)
dlist <- include.self(dlist)
dlist.d <- nbdists(dlist, coords)
col.W.tango <- nb2listw(dlist, glist=lapply(dlist.d, function(x){exp(-x)}), style="C")
tango.stat(data76_OE, col.W.tango, zero.policy=TRUE)

             [,1]
[1,] 0.0003874411

tango.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, list=col.W.tango, zero.policy=TRUE)

Tango's test of global clustering 

    Type of boots.: parametric 
    Model used when sampling: Poisson 
    Number of simulations: 100 
    Statistic:  0.0003874411 
    p-value :  0.00990099 

boot()関数を用いてブートストラップ法により標本を生成し、検定を行う。

data76_tn_pb <- boot(data76_OE, statistic=tango.pboot, sim="parametric", ran.gen=poisson.sim, R=100, listw=col.W.tango, zero.policy=TRUE)
plot(data76_tn_pb)

Whittermoreの統計量

空間オブジェクト間の距離\(d_ij\)を要素とする距離行列\(D\)を用いて、空間属性の集積性を仮説検定する。

col.W.whitt <- col.W.tango
whittermore.stat(data76_OE, col.W.whitt, zero.policy=TRUE)

           [,1]
[1,] 0.02664723

whittermore.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, list=col.W.whitt, zero.policy=TRUE)

Whittermore's test of global clustering 

    Type of boots.: parametric 
    Model used when sampling: Poisson 
    Number of simulations: 100 
    Statistic:  0.02664723 
    p-value :  0.00990099 

boot()関数を用いてブートストラップ法により標本を生成し、検定を行う。

data76_wt_pb <- boot(data76_OE, statistic=whittermore.pboot, sim="parametric", ran.gen=poisson.sim, R=100, listw=col.W.whitt, zero.policy=TRUE)
plot(data76_wt_pb)

Besag-Newell’s test

空間データの観測値\(O_i\)の発生が非常にまれであり、負の二項分布などに従うような場合に用いられる。

boot()関数を用いてブートストラップ法により標本を生成し、検定を行う。

data76_bn_perboot <- boot(data76_OE, statistic=besagnewell.boot, R=100, k=20)
plot(data76_bn_perboot)

OpenshawのGAM

opgam()関数を用いてOpenshaw のGeographical Analysis Machine(GAM)を計算する。

格子データの交点ではないが、市区町村代表点データ上に、市区町村別出生率を用いてGAM を適用してみる。

thegrid <- as(data76_OE, "data.frame")[,c("x","y")]
data76_opg <- opgam(data=as(data76_OE, "data.frame"), thegrid=thegrid,radius=20000, step=1000, alpha=0.05)
data76_opg

plot(data76_OE$x,data76_OE$y, xlab="Easting", ylab="Northing", cex=2,cex.axis=1.2, cex.lab=1.2)
points(data76_opg$x, data76_opg$y, col="red", pch="*", cex=4)

Kulldorffの空間スキャン統計量

最大尤度比となるウィンドゥをクラスターの候補と考える方法。最大尤度比を求めるために、モンテカルロシミュレーションで求めた値を用いて、統計的有意性を判断する。

kullnagar.stat()関数を用いて以下のように計算できる。

data76_OE <- data.frame(Observed=data76$n.birth)
data76_OE <- cbind(data76_OE,  
Expected=data76$pop*sum(data76$n.birth)/sum(data76$pop), 
x=data76$Easting, y=data76$Northing)
dist <- (data76_OE$x-data76_OE$x[1])^2+(data76_OE$y-data76_OE$y[1])^2
index<-order(dist)
kullnagar.stat(data76_OE[index,], fractpop=.5)

      value        size 
5.33617e+75 3.00000e+01