第05回演習マニュアル

演習の目的

  • R言語を使って,空間重み付け行列やMoran’s I,Geary’s Cなどの空間的自己相関指標を計算する.

演習課題1

  • R言語の空間統計パッケージspdepで利用可能なコロンバス市の犯罪発生データを使った演習.

Rの起動とパッケージの読み込み

  • Rを起動し,「Spdep」パッケージを読み込みます.
  • パッケージを読み込むには,「パッケージ」→「パッケージの読み込み」→「spdep」を選んでください.
  • パッケージが読み込まれると,R consoleに以下のように表示されます.
> local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
要求されたパッケージ tripack をロード中です
要求されたパッケージ maptools をロード中です
要求されたパッケージ foreign をロード中です
要求されたパッケージ sp をロード中です
要求されたパッケージ SparseM をロード中です
[1] "SparseM library loaded"
要求されたパッケージ boot をロード中です

地図データの表示

  • まずは,地図データを表示してみましょう.ここでは,空間統計データとしてよく知られるコロンバス市の犯罪発生データを用います.
  • 以下の手順で,R concoleに入力していきましょう(>以下がコマンド入力部分です).
> x <- read.shape(system.file("shapes/columbus.shp", package="maptools")[1])
Shapefile type: Polygon, (5), # of Shapes: 49
> plot(x)
Warning message:
'plot.Map' は廃止予定です
'plot.Spatial' を代わりに使って下さい
help("Deprecated") と help("maptools-deprecated") を見て下さい 
> n <- attr(x$Shapes,'nshps')
> nParts <- integer(n)
> for (i in 1:n) nParts[i] <- attr(x$Shapes[[i]], "nParts")
> table(nParts)
nParts
 1 
49 
> cols <- c("azure", "blue", "orange")
> fgs <- cols[nParts]
> plot(x, fg=fgs)
Warning message:
'plot.Map' は廃止予定です
'plot.Spatial' を代わりに使って下さい
help("Deprecated") と help("maptools-deprecated") を見て下さい 
> res <- plot(x, auxvar=x$att.data$CRIME)
Warning message:
'plot.Map' は廃止予定です
'plot.Spatial' を代わりに使って下さい
help("Deprecated") と help("maptools-deprecated") を見て下さい 
> str(res)
List of 3
 $ ramp     : chr [1:5] "#FFCCCC" "#FF9999" "#FF6666" "#FF3333" ...
 $ col.class: int [1:49] 1 1 3 3 4 2 1 3 3 3 ...
 $ breaks   : num [1:6]  0.178 19.023 29.326 39.025 53.161 ...
> 
  • すると,以下のような地図が表示されます.この地図では,1000世帯あたりの家宅侵入(窃盗)犯罪発生件数及び自動車盗難件数の合計が多いほど濃い赤色,少ないほど薄い赤色で示しています.

 left

隣接行列の作成

  • 隣接行列の作成方法としてspdepでは,
  1. ドロネー三角網を用いる方法,
  2. 近隣k地点に隣接関係があると定義する方法,
  3. 周辺の一定距離l以内の地点に隣接関係があると定義する方法,及び
  4. 隣り合うポリゴンに隣接関係があると定義する方法,
    がサポートされている.
  • colunbus.shpの属性データとして,columbusとcoordsがある.
    • 前者の属性の内容は,spdep.pdfを参照のこと.ただし,columbusに含まれるXとYはshpファイルの座標とは一致しない(デジタイズした際のものであるらしい)ため,そのXY座標はcoordに含まれている.
    • 関連して,columbus.shpの位相関係はpolysというファイルに予め定義されている.
    • Rのmaptoolsでは,shpファイルの位相関係を独自の形式で保管している.ポリゴンデータの場合,以下の手順でデータ形式を変換するが,これを既に行ったものがpolysである(ので,以下のコマンドはcolumbus.shpでは実行しなくてよい).

ドロネー三角網を用いる方法 (tri2nb)

  • ドロネー三角網とは,最近隣の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")

 left

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

  • まず,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")

 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")

 left

周辺の一定距離l以内の地点に隣接関係があると定義する方法

  • この方法は,ある点から半径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=""))

 left

隣り合うポリゴンに隣接関係があると定義する方法(polu2nb)

  • この方法は,ポリゴン同士が隣り合う場合に隣接関係があるとする方法です.
  • 下図のように,ドロネー三角網図を用いた場合,ポリンゴンが隣り合わない場合にも隣接関係があるとみなされることがありますので,ポリゴンデータを利用する場合にこの方法は有益です.
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: drauney triangulation weights \n Black: and polygon generated queen weights")

 left

Moran’s IとGeary’s Cの計算

  • さて,いよいよ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:横浜市の公示地価(住宅地地価)の空間的自己相関指標を計算

データのダウンロード

  • 横浜市のデータをここからダウンロードしてください.
  • 解凍後,データが「C:\Program Files\R\R-2.2.1\yoko」にあるかどうか確認し,ない場合はフォルダを移動させてください.

データの読み込みと表示

  • まず,シェープファイルを読み込みます.
    > h14lp <- read.shape("yoko/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")
  • すると,以下のような図が表示されます.
     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年住宅地地価は強い正の空間的自己相関を持つことが分かります.
  • 隣接行列の定義の仕方でどのように異なるのか,また誤差項の空間的自己相関についてはどうかを計算してみましょう.

演習課題レポートの提出

  • 以上の二つの課題についてレポートを作成し,提出してください.
  • 計算結果だけでなく,地価や誤差項の空間分布を具体的に図示し,空間的自己相関との関係を考察すること.
  • 分量:A4 ページ数は特に指定しないが,読みやすくなるよう心がけること.
  • 提出期限・方法:11月27日(火)の授業時間中にレポートを提出すること.授業時間内の提出が困難な場合,提出期限である授業開始前までに担当者にメールで送付すること.

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


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