* 空間モデリング(空間的自己相関) [#jb69b99d]

*** 演習の目的 [#c22f8cf8]
- R言語を使って,空間重み付け行列やMoran’s I,Geary’s Cなどの空間的自己相関指標を計算する.

** 演習課題1 [#r01ff7af]
- R言語の空間統計パッケージspdepで利用可能なコロンバス市の犯罪発生データを使った演習.

** ライブラリの読み込み [#y2391dc2]
- Rで空間的自己相関指標を計算するために、spdepというパッケージを読み込みます。
- Windowsの人は、以下の手順でパッケージを読み込んでください。
-- RGui画面上のメニューバーから、「パッケージ」→「パッケージのインストール」→「CRAN mirror」から適当なミラーサイトを選ぶ→「Packages」リストから「spdep」を選んでOKボタンをクリックする
-- 再びRGui画面上のメニューバーから、「パッケージの読み込み」→選択画面から「spdepと」を読み込む
- Macの人は、以下の手順でパッケージを読み込んでください。
-- RGui画面上のメニューバーから、「パッケージとデータ」→(初めてパッケージをインストールする場合は)「パッケージインストーラ」→「CRAN mirror」から適当なミラーサイトを選ぶ→「spdep」を選び、「パッケージ間の依存関係を解決」をチェックし、「選択をインストール」をクリック
-- 再びRGui画面上のメニューバーから、「パッケージとデータ」→「パッケージマネージャ」→「spdep」にチェックを入れる。
- すでに一度、パッケージをインストールしている場合は、次回から以下のコマンドをコンソール上で入力すればよい。

 library(spdep)

- また、地図描画のために、パッケージclassIntを読み込む。

 library(classInt)

** 地図データの表示 [#s54eb939]
- まずは,地図データを表示してみましょう.ここでは,空間統計データとしてよく知られるコロンバス市の犯罪発生データを用います.
- 以下の手順で,R concoleに入力していきましょう(>以下がコマンド入力部分です).

 x <- readShapePoly(system.file("shapes/columbus.shp", package="maptools")[1])
 pal <- c("white","red")
 q_x <- classIntervals(x$CRIME, n=4,style="quantile")
 plot(q_x, pal=pal)
 q_x_Col <- findColours(q_x,pal)
 plot(x,col=q_x_Col)

- すると,以下のような地図が表示されます.この地図では,1000世帯あたりの家宅侵入(窃盗)犯罪発生件数及び自動車盗難件数の合計が多いほど濃い赤色,少ないほど薄い赤色で示しています.~

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image002.jpg, left);~
~
** 隣接行列の作成 [#fbe373b4]
- 隣接行列の作成方法としてspdepでは,

+ ドロネー三角網を用いる方法,
+ 近隣k地点に隣接関係があると定義する方法,
+ 周辺の一定距離l以内の地点に隣接関係があると定義する方法,及び
+ 隣り合うポリゴンに隣接関係があると定義する方法,~
がサポートされている.

- colunbus.shpの属性データとして,columbusとcoordsがある.
-- 前者の属性の内容は,[[spdep.pdf>http://cran.r-project.org/web/packages/spdep/spdep.pdf]]を参照のこと.ただし,columbusに含まれるXとYはshpファイルの座標とは一致しない(デジタイズした際のものであるらしい)ため,そのXY座標はcoordに含まれている.
-- 関連して,columbus.shpの位相関係はpolysというファイルに予め定義されている.
-- Rのmaptoolsでは,shpファイルの位相関係を独自の形式で保管している.ポリゴンデータの場合,以下の手順でデータ形式を変換するが,これを既に行ったものがpolysである(ので,以下のコマンドはcolumbus.shpでは実行しなくてよい).

*** ドロネー三角網を用いる方法 (tri2nb) [#vbc37f76]
- ドロネー三角網とは,最近隣の2点同士を結び,三角形のネットワークを形成するものでしたね.
- これはボロノイ図をつくる際にも使われます.tri2nbを実行するには,tripackが必要ですが,このパッケージには,実際にボロノイ図を作るためのコマンドも用意されています.
- 以下の内容を実行すると,下図の結果が得られます.これがドロネー三角網図です.

 data(columbus)
 col.tri.nb <- tri2nb(coords, row.names=rownames(columbus))
 plot(polys, border="grey")
 plot(col.tri.nb, coords, add=TRUE)
 title(main="Raw triangulation links")

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image004.jpg, left);~
~

*** 近隣k地点に隣接関係があると定義する方法 (knn2nb) [#a8b26399]

- まず,k=4の場合として以下のコマンドを実行しましょう.
- ここで,knearneighは,座標行列coordsをもちいて近隣k=4個の隣接関係を定義するコマンドです.結果は下図に示されています.

 col.knn <- knearneigh(coords, k=4)
 plot(polys, border="grey")
 plot(knn2nb(col.knn), coords, add=TRUE)
 title(main="K nearest neighbours, k = 4")

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image006.jpg, left);~
~

- ここで,kの違いによる隣接関係の決まり方の違いを視覚的に理解してみましょう.
- diffnbコマンドを使い,両者の差異を示します.(図示する前に,一度plot図を消しておくとよい)下図では,k=3が黒線,k=4が赤点線で示されています.そもそも,(1)のドロネー三角網図の結果とは大きく異なりますね.

 knn3 <- knearneigh(coords, 3)
 knn4 <- knearneigh(coords, 4)
 nb3 <- knn2nb(knn3, row.names=rownames(columbus))
 nb4 <- knn2nb(knn4, row.names=rownames(columbus))
 diffs <- diffnb(nb3, nb4)
 plot(polys, border="grey")
 plot(nb3, coords, add=TRUE)
 plot(diffs, coords, add=TRUE, col="red", lty=2)
 title(main="Plot of third (black) and fourth (red)\nnearest neighbours")

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image008.jpg, left);~
~

*** 周辺の一定距離l以内の地点に隣接関係があると定義する方法 [#f800fd64]
- この方法は,ある点から半径l以内の距離にある(半径lの円の中に含まれる)地点を隣接関係があるものとみなします.(ポリゴンならポリゴンの中心点が含まれているかどうかを判断する場合が多い)
- このとき,dnearneighを使って,半径規模を指定します.
- l=1.0とした場合,以下の内容を実行すると,下図の結果が得られます.

 col.nb1 <- dnearneigh(coords, 0, 1.0, row.names=rownames(columbus))
 plot(polys, border="grey")
 plot(col.nb1, coords, add=TRUE)
 title(main=paste("Distance based neighbours 0-1 distance units", sep=""))

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image010.jpg, left);~
~

*** 隣り合うポリゴンに隣接関係があると定義する方法(poly2nb) [#e24bf442]
- この方法は,ポリゴン同士が隣り合う場合に隣接関係があるとする方法です.
- 下図のように,ドロネー三角網図を用いた場合,ポリンゴンが隣り合わない場合にも隣接関係があるとみなされることがありますので,ポリゴンデータを利用する場合にこの方法は有益です.

 pnb <- poly2nb(polys)
 col.tri.nb <- tri2nb(coords, row.names=rownames(columbus))
 dpnb <- diffnb(pnb, col.tri.nb)
 plot(polys, border="grey")
 plot(dpnb, coords, add=TRUE, col="red")
 plot(polys, border="grey")
 plot(col.tri.nb, coords, add=TRUE)
 plot(dpnb, coords, add=TRUE, col="red")
 title(main=" Red: Delaunay Triangulation Weights \n Black: and polygon generated queen weights")

&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image012.jpg, left);~
~

*** Moran’s IとGeary’s Cの計算 [#kc4b39bb]
- さて,いよいよMoran’s IとGeary’s Cを計算してみましょう.
- ここでは,3)-(4)の方法を用いて隣接行列を定義します.
- ここで,spNamedVec("CRIME", columbus)は行列columbusの"CRIME"という名前の列をベクトルとして返すことを意味し,nb2listw(pnb, style="W")は近隣行列pnbを行和で標準化した重み付け行列を用いた場合のリストを返します.
- nb2listwでの重み付け行列のオプションについては,spdep.pdfのnb2listwの項目を参照してください.

 > moran.test(spNamedVec("CRIME", columbus), nb2listw(pnb, style="W"))
         Moran's I test under randomisation
 data:  spNamedVec("CRIME", columbus)  
 weights: nb2listw(pnb, style = "W")  
 Moran I statistic standard deviate = 5.5894, p-value = 1.139e-08
 alternative hypothesis: greater 
 sample estimates:
 Moran I statistic       Expectation          Variance 
        0.50018856       -0.02083333        0.00868929 

- 同様に,Geary’s C統計量は以下のように求めます.

 > geary.test(spNamedVec("CRIME", columbus), nb2listw(pnb, style="W"))
         Geary's C test under randomisation
 data:  spNamedVec("CRIME", columbus) 
 weights: nb2listw(pnb, style = "W") 
 Geary C statistic standard deviate = -4.7431, p-value = 1.053e-06
 alternative hypothesis: less 
 sample estimates:
 Geary C statistic       Expectation          Variance 
       0.540528203       1.000000000       0.009384264


*** 演習課題2:横浜市の公示地価(住宅地地価)の空間的自己相関指標を計算 [#gf62b0df]
** データのダウンロード [#h8d0629c]
- 横浜市のデータを[[ここ>http://web.sfc.keio.ac.jp/~maunz/SI07/yoko.zip]]からダウンロードしてください.
- 解凍後,データが「C:\Program Files\R\R-2.x.x\yoko」にあるかどうか確認し,ない場合はフォルダを移動させてください.

** データの読み込みと表示 [#tad453a6]
- まず,シェープファイルを読み込みます.
 > h14lp <- read.shape("h14lp.shp")
 Shapefile type: Point, (1), # of Shapes: 691
- 次に,ポイントデータに変換します.
 > h14lp.pnt <- Map2points(h14lp)
- XY座標のみを取り出したcoords行列を作成します.
 > coords <- matrix(0,nrow(h14lp$att.dat),2)
 > coords[,1] <- h14lp$att.dat$X
 > coords[,2] <- h14lp$att.dat$Y
- ドロネー三角網で隣接行列lp.tri.nbを計算します.
 > lp.tri.nb <- tri2nb(coords, row.names=rownames(h14lp))
 > plot(h14lp, col="grey")
 Warning message:
 'plot.Map' は廃止予定です
 'plot.Spatial' を代わりに使って下さい
 help("Deprecated") と help("maptools-deprecated") を見て下さい 
 > plot(lp.tri.nb, coords, add=TRUE)
 > title(main="Raw triangulation links")
- すると,以下のような図が表示されます.~
&ref(http://web.sfc.keio.ac.jp/~maunz/SII04/SDM04-07ex.files/image023.gif, left);~
~
- 最後に,Moran’s Iを計算します.
 > moran.test(spNamedVec("H14LP", h14lp$att.data), nb2listw (lp.tri.nb,style="W"))
 
         Moran's I test under randomisation 
 
 data:  spNamedVec("H14LP", h14lp$att.data)  
 weights: nb2listw(lp.tri.nb, style = "W")  
  
 Moran I statistic standard deviate = 33.1096, p-value < 2.2e-16
 alternative hypothesis: greater 
 sample estimates:
 Moran I statistic       Expectation          Variance 
       0.732513488      -0.001449275       0.000491407
- この結果から,横浜市内の平成14年住宅地地価は強い正の空間的自己相関を持つことが分かります.
- 隣接行列の定義の仕方でどのように異なるのか,また誤差項の空間的自己相関についてはどうかを計算してみましょう.


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


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