* ESTRELA連載 第76回 [#ra9c50f8]

** 今回使うデータ [#a59f19a0]
- [[データ>http://web.sfc.keio.ac.jp/~maunz/ESTRELA/76.zip]]

** 地図表示 [#efdd6723]
- パッケージの読み込み
 library(spdep)
 library(DCluster)
- 地図データの読み込み
 jpn_pref <- readShapePoly("jpn_pref.shp",IDvar="PREF_CODE")
 pref_pnt <- readShapePoints("pref_gov.shp")
 hd06 <- read.table("76dat.csv", sep=",", header=T)
 summary(hd06)
- 都道府県境界ポリゴンデータへの属性データhd06のマッチング
 ID.match <- match(jpn_pref$PREF_CODE, hd06$PREF_CODE)
 jpn_hd06 <- hd06[ID.match,]
 jpn_pref_hd06 <- spCbind(jpn_pref, jpn_hd06)
 summary(jpn_pref_hd06)

** 経験ベイズ法(つづき) [#u7b512db]
- ポアソンーガンマモデル
 jpn_hd06_pm <- probmap(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100)
 summary(jpn_hd06_pm)
 jpn_hd06_sm <- empbaysmooth(jpn_hd06$HD06, jpn_hd06_pm$expCount)
 jpn_hd06_sm
- 対数正規モデル
 jpn_hd06_ln <- lognormalEB(jpn_hd06$HD06, jpn_hd06_pm$expCount)
 jpn_hd06_ln

** 空間集積性 [#sa7e492c]
*** 分析用データ [#ea11dbce]
- データの読み込みと作成
 # 配布したデータ名が「data76_test2.csv」となっている場合は、ファイル名を「data76.csv」に変更してください。
 data76 <- read.table("data76.csv", sep=",", header=T, 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)

*** 空間が均質であるかどうかの検定 [#zfce6d5b]
- χ2検定
 achisq.stat(data76_OE, lambda=1)
 achisq.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100)
 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)
 pottwhitt.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100)
 data76_pw_pb <- boot(data76_OE, statistic=pottwhitt.pboot, sim="parametric", ran.gen=poisson.sim, R=100)
 plot(data76_pw_pb)

- Tangoの一般的なクラスタリング検定
 data76_OE <- cbind(data76_OE, x=data76_spdf$Easting, y=data76_spdf$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)
 tango.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, list=col.W.tango, zero.policy=TRUE)
 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の統計量
 col.W.whitt <- col.W.tango
 whittermore.stat(data76_OE, col.W.whitt, zero.policy=TRUE)
 whittermore.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, list=col.W.whitt, zero.policy=TRUE)
 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)

*** クラスター位置の検出 [#ec2e79d6]
- Openshawの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")
 points(data76_opg$x, data76_opg$y, col="red", pch="*")

- Besag-Newell's test
 data76_bn_perboot <- boot(data76_OE, statistic=besagnewell.boot, R=100, k=20)
 plot(data76_bn_perboot)
 # besagnewell.stat(data76_OE, k=30)

- Kulldorffの空間スキャン統計量
 # Permutation model
 data76_kn_perboot <- boot(data76_OE, statistic=kullnagar.boot, R=100, fractpop=.5)
 plot(data76_kn_perboot)
 # Multinomial model
 data76_m_pboot <- boot(data76_OE, statistic=kullnagar.pboot, sim="parametric", ran.gen=multinom.sim, R=100, fractpop=0.5)
 plot(data76_kn_mboot)
 # Poisson model
 data76_kn_pboot <- boot(data76_OE, statistic=kullnagar.pboot, sim="parametric", ran.gen=poisson.sim, R=100, fractpop=0.5)
 plot(data76_kn_pboot)
 # Poisson-Gamma model
 data76_kn_pgboot <- boot(data76_OE, statistic=kullnagar.pboot, sim="parametric",ran.gen=negbin.sim, R=100, fractpop=0.5)
 plot(data76_kn_pgboot)
 ##
 dist <- (data76_OE$x-data76_OE$x[1])^2+(data76_OE$y-data76_OE$y[1])^2
 index<-order(dist)
 kullnagar.stat(sids[index,], fractpop=.5)

- Stone's test
 stone.stat(data76_OE, region=which(row.names(data76_OE)=="20"), lambda=1)
 stone.test(Observed~offset(log(Expected)), data76_OE, model="poisson", R=100, region=which(row.names(data76_OE)=="20"), lambda=1)
 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)

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS