ESTRELA連載 第74回

今回使うデータ

地図表示

  • パッケージの読み込み
    library(spdep)
  • 地図データの読み込みと表示
  • ESTRELA記事では、read.shape()とMap2poly()を組み合わせた方法を紹介していますが、R最新版を利用している場合は、readShapePoly?()を使う方法を適用して下さい。
    jpn.pref <- readShapePoly("jpn_pref.shp",IDvar="PREF_CODE")
    # または
    # jpn.shp <- read.shape("jpn_pref.shp",IDvar="PREF_CODE")
    # jpn.pref <- Map2poly(jpn.shp)
    pref.pnt <- readShapePoints("pref_gov.shp")
    # または
    # pref.shp <- read.shape("pref_gov.shp")
    # pref.pnt <- Map2points(pref.shp)
    plot(jpn.pref)
    points(pref.pnt)

空間隣接行列

  • 座標テーブルの作成
    pref_gov <- read.table("pref_gov.txt",",",header=T, row.names=2)
    coords <- as.matrix(cbind(pref_gov$X,pref_gov$Y))
    # coords <- matrix(0,nrow(pref_gov),2)
    # coords[,1] <- pref_gov$X
    # coords[,2] <- pref_gov$Y
  • ドロネー三角網図
    pref.tri.nb <- tri2nb(coords, row.names=rownames(pref_gov))
    plot(pref.tri.nb, coords, add=T)
  • 最近隣k地点
    pref.knn <- knearneigh(coords, k=4)
    pref.knn.nb <- knn2nb(pref.knn, row.names=rownames(pref_gov))
    plot(jpn.pref)
    plot(pref.knn.nb, coords, add=TRUE)
  • k=4とk=3の違い
    pref.knn3 <- knearneigh(coords, k=3)
    pref.knn3.nb <- knn2nb(pref.knn3)
    diffs <- diffnb(pref.knn3.nb, pref.knn.nb)
    plot(jpn.pref)
    plot(pref.knn3.nb, coords, add=TRUE)
    plot(diffs, coords, add=TRUE, col="red", lty=2)
  • 一定距離以内に隣接関係があると定義する場合(r=2)
    pref.r.nb <- dnearneigh(coords, 0, 2, row.names=rownames(pref_gov))
    plot(jpn.pref)
    plot(pref.r.nb, coords, add=TRUE)
  • ポリゴンの隣接関係で定義する場合
    honshu <- read.shape("honshu.shp")
    honshu <- Map2poly(honshu)
    honshu.poly.nb <- poly2nb(honshu)
    plot(honshu)
    plot(honshu.poly.nb, coords, add=TRUE)
    dpnb <- diffnb(pref.poly.nb, pref.tri.nb)
    plot(jpn.pref)
    plot(pref.tri.nb, coords, add=TRUE)
    plot(dpnb, coords, add=TRUE, col="red")

空間的自己相関

  • Moran's I
    moran.test(pref.pnt$diabetes, nb2listw (pref.tri.nb,style="W"))
    moran.test(pref.pnt$geriatric, nb2listw (pref.tri.nb,style="W"))
    moran.test(pref.pnt$malignant, nb2listw (pref.tri.nb,style="W"))
    moran.test(pref.pnt$hypertensi, nb2listw (pref.tri.nb,style="W"))
  • Geary's C
    geary.test(pref.pnt$diabetes, nb2listw(pref.tri.nb, style="W"))
    geary.test(pref.pnt$geriatric, nb2listw(pref.tri.nb, style="W"))
    geary.test(pref.pnt$malignant, nb2listw(pref.tri.nb, style="W"))
    geary.test(pref.pnt$hypertensi, nb2listw(pref.tri.nb, style="W"))
  • Local Moran's I
    LMI1 <- localmoran(pref.pnt$diabetes, nb2listw(pref.tri.nb, style="W"))
    pref.lm <- data.frame(cbind(LMI1[,1],(pref.pnt$diabetes- 
    mean(pref.pnt$diabetes))/sd(pref.pnt$diabetes)),row.names=pref.pnt$KENCODE)
    colnames(pref.lm) <- c("Ii","standardized diabetes")
    pref.lm
    plot(pref.lm,xlab="Local Moran's I", ylab="Standardized diabetes")
    text(pref.lm,rownames(pref.lm))

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-05-11 (火) 23:03:35 (3479d)