ESTRELA連載 第75回

お詫び

  • 本文中、経験ベイズ推定値のMoran's I値を、経験ベイズ推定量と間違って記述してしまいました。

今回使うデータ

地図表示

  • パッケージの読み込み
    library(spdep)
  • 地図データの読み込み
    jpn_pref <- readShapePoly("jpn_pref.shp",IDvar="PREF_CODE")
    pref_pnt <- readShapePoints("pref_gov.shp")
    hd06 <- read.table("75dat.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)
  • 都道府県庁所在地ポイントデータへの属性データhd06のマッチング
    # 演習では使いません
    ID.match2 <- match(pref_pnt$KENCODE, hd06$PREF_CODE)
    jpn_hd06_2 <- hd06[ID.match2,]
    pref_pnt_hd06 <- spCbind(pref_pnt, jpn_hd06_2)
    summary(pref_pnt_hd06)

確率地図の作成

  • probmap()関数
    jpn_hd06_pm <- probmap(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100)
    summary(jpn_hd06_pm)
  • 主題図の作成
    library(classInt)
    brks1 <- c(0,1500, 2500, 3500, 4500, 20000)
    brks2 <- c(0,100,125,150,175,200)
    brks3 <- c(0, 0.2, 0.6, 0.9, 1.0)
    cols <- grey(6:2/6)
  • 観測値
    plot(jpn_pref, col=cols[findInterval(jpn_pref_hd06$HD06, brks1, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks1), bty="n")
    title(main="Observed number of heart disease death (2006)")
  • 期待値(jpn_hd06_pm$expCount)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_pm$expCount, brks1, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks1), bty="n")
    title(main="Expected number of heart disease death (2006)")
  • 粗率(jpn_hd06_pm$raw)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_pm$raw, brks2, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks2), bty="n")
    title(main="Raw rate of heart disease per 100,000 (2006)")
  • 相対危険度(jpn_hd06_pm$relRisk)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_pm$relRisk, brks2, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks2), bty="n")
    title(main="Relevant risk of heart disease per 100,000 (2006)")
  • ポアソン確率(jpn_hd06_pm$pmap)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_pm$pmap, brks3, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks3), bty="n")
    title(main="Poisson probability of heart disease (2006)")
    hist(jpn_hd06_pm$pmap)

経験ベイズ法

  • 経験ベイズ推定(Marshallのグローバルな経験ベイズ推定値)
    jpn_hd06_ebg <- EBest(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_ebg$estmm, brks2, all.inside=TRUE)])
    legend("topleft", fill=cols, legend=leglabs(brks2), bty="n")
    title(main="Empirical Bayes estimates of heart disease death rate per 100,000 (2006)", cex.main=0.9)
  • localな経験ベイズ推定(Marshallのローカルな経験ベイズ推定値)
    pref_gov <- read.table("pref_gov.txt",",",header=T, row.names=2)
    coords <- matrix(0,nrow(pref_gov),2)
    coords[,1] <- pref_gov$X
    coords[,2] <- pref_gov$Y
    pref.knn <- knearneigh(coords, k=4)
    pref.knn.nb <- knn2nb(pref.knn, row.names=rownames(pref_gov))
    jpn_hd06_ebl <- EBlocal(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100, pref.knn.nb)
    plot(jpn_pref, col=cols[findInterval(jpn_hd06_ebl$est, brks2, all.inside=TRUE)])
    plot(pref.knn.nb, coords, add=TRUE)
    legend("topleft", fill=cols, legend=leglabs(brks2), bty="n")
    title(main="Local empirical Bayes estimates of heart disease death rate per 100,000 (2006)", cex.main=0.8)
  • 経験ベイズ推定値によるMoran's Iの繰り返し検定
    EBImoran.mc(jpn_pref_hd06$HD06, jpn_pref_hd06$POPJ06/100, 
    nb2listw(pref.knn.nb, style="B", zero.policy=TRUE), 
    nsim=999, zero.policy=TRUE)
  • 通常のMoran's Iの繰り返し検定
    moran.mc(jpn_pref_hd06$HD06/(jpn_pref_hd06$POPJ06/100), 
    nb2listw(pref.knn.nb, style="B", zero.policy=TRUE), 
    nsim=999, zero.policy=TRUE)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-13 (月) 00:50:04 (3328d)