空間モデリング(09年度春学期第07回演習)

演習の目的

  • 空間データの分類方法のうち、local Moran's I、階層クラスター分析、非階層クラスター分析について、統計言語Rを使って演習を行う。
  • 今回は、日本の都道府県データを使って演習します。データは、SFC-SFSの授業ページからダウンロードしてください。
  • 都道府県データ(ポリゴン、都道府県庁所在地ポイント)の属性として、人口(密度)、社会生活統計指標 -都道府県の指標- の要因別死亡者数(2005年)のデータが含まれています。

パッケージの読み込み

  • 以下のパッケージをインストールし読み込みます。
    library(spdep)
    library(maptools)
    library(splancs)
    library(rgdal) # Macではうまく読み込めないかもしれません
    library(gpclib)
    library(classInt)
    library(RColorBrewer)

主題図の作成

  • 都道府県マップを読み込む。
    pref_poly <- readShapePoly("pref_poly_COD.shp",IDvar="PREF_CODE")
    plot(pref_poly)
  • 四分位区分により人口密度をプロットする。
    pal <- gray.colors(4,0.9,0.3,2.2)
    q_popd <- classIntervals(pref_poly$Pop_Dens, n=4, style="quantile")
    plot(q_popd, pal=pal)
    q_popd_Col <- findColours(q_popd,pal)
    # plot(pref_poly$Pop_Dens,col=q_popd_Col)
    plot(pref_poly,col=q_popd_Col)
    title("Population Density (quantile)")
    legend("topleft",fill=attr(q_popd_Col,"palette"),
    legend=names(attr(q_popd_Col,"table")),bty="n")
  • 標準偏差区分により人口密度をプロットする。
    sd_popd <- classIntervals(pref_poly$Pop_Dens, style="sd")
    plot(sd_popd, pal=pal)
    sd_popd_Col <- findColours(sd_popd,pal)
    plot(pref_poly,col=sd_popd_Col)
    title("Population Density (standard deviations)")
    legend("topleft",fill=attr(sd_popd_Col,"palette"),
    legend=names(attr(sd_popd_Col,"table")),bty="n")
  • 等間隔区分により人口密度をプロットする。
    eq_popd <- classIntervals(pref_poly$Pop_Dens, n=5, style="equal")
    plot(eq_popd, pal=pal)
    eq_popd_Col <- findColours(eq_popd,pal)
    plot(pref_poly,col=eq_popd_Col)
    title("Population Density (equql)")
    legend("topleft",fill=attr(eq_popd_Col,"palette"),
    legend=names(attr(eq_popd_Col,"table")),bty="n")
  • 階層クラスタリングによる区分に基づき人口密度をプロットする。
    hc_popd <- classIntervals(pref_poly$Pop_Dens, n=5, style="hclust")
    plot(hc_popd, pal=pal)
    hc_popd_Col <- findColours(hc_popd,pal)
    plot(pref_poly,col=hc_popd_Col)
    title("Population Density (hierarchical clustering)")
    legend("topleft",fill=attr(hc_popd_Col,"palette"),
    legend=names(attr(hc_popd_Col,"table")),bty="n")
  • 非階層クラスタリングによる区分に基づき人口密度をプロットする。
    km_popd <- classIntervals(pref_poly$Pop_Dens, n=5, style="kmeans")
    plot(km_popd, pal=pal)
    km_popd_Col <- findColours(km_popd,pal)
    plot(pref_poly,col=km_popd_Col)
    title("Population Density (K-means)")
    legend("topleft",fill=attr(km_popd_Col,"palette"),
    legend=names(attr(km_popd_Col,"table")),bty="n")

Local Moran's I

  • データを読み込み、隣接行列を定義する。
    pref_pnt_COD <- readShapePoints("pref_pnt_COD.shp")
    coords <- matrix(0,nrow(pref_pnt_COD),2)
    coords[,1] <- pref_pnt_COD$X
    coords[,2] <- pref_pnt_COD$Y
    pref.tri.nb <- tri2nb(coords, row.names=rownames(pref_pnt_COD$KENCODE))
    plot(pref_pnt_COD)
    plot(pref.tri.nb, coords, add=TRUE)
  • Moran'sI Geary's Cをそれぞれ計算する
    moran.test(pref_pnt_COD$Cancer,nb2listw(pref.tri.nb,style="W"))
    geary.test(pref_pnt_COD$Cancer,nb2listw(pref.tri.nb,style="W"))
  • Local Moran's Iを計算してみましょう。
    resI <- localmoran(pref_pnt_COD$Cancer, nb2listw(pref.tri.nb))
    resI
    COD.lm <- data.frame(cbind(resI[,1],(pref_pnt_COD$Cancer- 
    mean(pref_pnt_COD$Cancer))/sd(pref_pnt_COD$Cancer)),
    row.names=pref_pnt_COD$KENCODE)
    colnames(COD.lm) <- c("Ii","standardized totcon")
    COD.lm
    plot(COD.lm,xlab="Local Moran's I", ylab="Standardized Cancer Death")
    text(COD.lm,rownames(COD.lm))

階層クラスター分析

  • ライブラリの読み込み
    library(MASS)
  • データの読み込み
    pref_COD.dat <- read.table("pref_COD.csv",sep=",",header=T)
    summary(pref_COD.dat)
    pref_COD <- pref_COD.dat[,5:8]
    rownames(pref_COD) <- pref_COD.dat[,2]
  • hclust()関数による階層クラスター分析
  • ユークリッド距離を用いたクラスタリング
  • 最近隣距離法
    pref_COD.hc1 <- hclust(dist(pref_COD), "single")
    pref_COD.hc1
    plot(pref_COD.hc1, hang=-1)
  • 最遠距離法
    pref_COD.hc2 <- hclust(dist(pref_COD), "complete") 
    pref_COD.hc2
    plot(pref_COD.hc2, hang=-1)
  • 群平均法
    pref_COD.hc3 <- hclust(dist(pref_COD), "average")
    pref_COD.hc3
    plot(pref_COD.hc3, hang=-1)
  • 重心法
    pref_COD.hc4 <- hclust(dist(pref_COD), "centroid")
    pref_COD.hc4
    plot(pref_COD.hc4, hang=-1)
  • メディアン法
    pref_COD.hc5 <- hclust(dist(pref_COD), "median")
    pref_COD.hc5
    plot(pref_COD.hc5, hang=-1)
  • ウォード法
    pref_COD.hc6 <- hclust(dist(pref_COD), "ward")
    pref_COD.hc6
    plot(pref_COD.hc6, hang=-1)
  • 最近隣距離法
    pref_COD.hc7 <- hclust(dist(pref_COD, "can"), "single")
    pref_COD.hc7
    plot(pref_COD.hc7, hang=-1)

非階層クラスター分析

  • kmean()関数を用いた分析
    pref_COD.km <- kmeans(pref_COD,4)
    summary(pref_COD.km)
    pref_COD.km

以上で今回の演習は終了です.


トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2009-07-16 (木) 20:07:05 (3657d)