ESTRELA連載 第76回

今回使うデータ

地図表示

  • パッケージの読み込み
    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)

経験ベイズ法(つづき)

  • ポアソンーガンマモデル
    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

空間集積性

分析用データ

  • データの読み込みと作成
    # 配布したデータ名が「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)

空間が均質であるかどうかの検定

  • χ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)

クラスター位置の検出

  • 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
Last-modified: 2010-09-16 (木) 11:38:58 (3230d)