第11回演習マニュアル

  • Rを起動し,splancs,spdep,DCluster,clusterパッケージをロードしてください.パッケージリストにない場合は,適宜インストールしてください.

データの準備1

  • 最初に使うデータは,ノースカロライナ州における乳幼児突然死症候群(SIDS)のデータ(nc.sids)です.DClusterパッケージで利用可能です.
  • nc.sids$SID74は'74年の乳幼児突然死症候群にかかった新生児数,nc.sids$BIR74は同年の新生児数です.Expectedは各都市での期待死亡数を意味します.
data(nc.sids)
sids<-data.frame(Observed=nc.sids$SID74)
sids<-cbind(sids, Expected=nc.sids$BIR74*sum(nc.sids$SID74)/sum(nc.sids$BIR74))
sids<-cbind(sids, x=nc.sids$x, y=nc.sids$y)

χ2統計量による点分布の一様性(ランダム性)の検定

  • achisq.statを使います.引数lambda=1としています.つまり,区画総数に対するポイント総数の比を1としています,
achisq.stat(sids, lambda=1)
$T
[1] 225.5723
$df
[1] 100
$pvalue
[1] 1.087968e-11

ローカルな空間集積性:Local Moran's I

  • spdepパッケージのlocalmoran関数を用いて計算できます.
  • ここでは,アフリカにおける紛争発生状況(1966-78年)のデータafconを使います.
    • totconが紛争件数を意味します.国間の隣接関係はpaper.nbというオブジェクトがすでに用意されています.
data(afcon)
oid <- order(afcon$id)
resI <- localmoran(spNamedVec("totcon", afcon), nb2listw(paper.nb))
printCoefmat(data.frame(resI[oid,], row.names=afcon$name[oid]),  
check.names=FALSE)

点データのクラスタリング

GAM: Geographical Analysis Machine (Openshaw's clustering)

  • opgamというコマンドを使います.
  • opgamコマンドに,データセット,検索半径(radius),グリッドの大きさ(step),有意水準(alpha)を設定します.
    • GAMの検定では,有意水準(alpha)をかなり小さめの値(例えば,p<0.002)に設定することが多い.
sidsgam<-opgam(data=sids, radius=30, step=10, alpha=.002)
plot(sids$x, sids$y, xlab="Easting", ylab="Northing")
points(sidsgam$x, sidsgam$y, col="red", pch="*")
radii <- rep(30, length(sidsgam$x))
symbols(sidsgam$x, sidsgam$y, circles=radii,add=T)

Besag and Newellの方法

  • Besag and Newell(1991)が提案したGAMの改良版で,空間クラスターの中心位置と円の大きさを,GAMに比べて絞っている.
  • opgam関数のisclusterオプションをiscluster=bn.isclusterとし,セルサイズk=20,有意水準alpha=0.05とします.
  • 円内部の点属性と対象地域全体の点属性を,ポアソン分布を用いてモンテカルロ検定するため,model="poisson",ブートストラップ数R=100,負の二項分布を計算するためのパラメータmle=calculate.mle(sids)を設定します.
  • ブートストラップを呼び出すパッケージ(boot)を使います.
library(boot)
bnresults<-opgam(data=sids, thegrid=sids[,c("x","y")], 
alpha=.05,iscluster=bn.iscluster, k=20, R=100, model="poisson", 
mle=calculate.mle(sids))
plot(sids$x, sids$y, main="Besag & Newell's method")
points(bnresults$x, bnresults$y, col="red", pch=19)

空間スキャン統計量

  • 授業では扱いません.作業メモ
niter<-100
kn.perboot<-boot(sids, statistic=kullnagar.boot, R=niter, fractpop=.2)
plot(kn.perboot)#Display results
kn.mboot<-boot(sids, statistic=kullnagar.pboot, sim="parametric",
ran.gen=multinom.sim, R=niter, fractpop=.2)
plot(kn.mboot)#Display results
kn.pboot<-boot(sids, statistic=kullnagar.pboot, sim="parametric",
ran.gen=poisson.sim, R=niter, fractpop=.2)
plot(kn.pboot)#Display results
kn.pgboot<-boot(sids, statistic=kullnagar.pboot, sim="parametric",
ran.gen=negbin.sim, R=niter, fractpop=.2)
plot(kn.pgboot)#Display results

カーネル関数とカーネルスムージング

カーネル関数

  • K関数を推定するだけであれば,spatstatパッケージのKest関数を使うことができます.
  • ランダムで一様なn個のポイントを作成するのに,runinfpoint関数が使えます.
pp <- runifpoint(50)
K<-Kest(pp,correction="Ripley")
plot(K)
    lty col
iso    1   1
theo   2   2
  • ここで,colは曲線の色で,1〜4に順に,黒,赤,緑,青となっています.
  • splancsパッケージのKhat関数を使って計算することもできます.
  • ここでは,青少年犯罪者データcardiffを使います.
data(cardiff)
s <- seq(2,30,2)
plot(s, sqrt(khat(as.points(cardiff), cardiff$poly, s)/pi) - s,
type="l", xlab="Splancs - polygon boundary", ylab="Estimated L",
ylim=c(-1,1.5))
newstyle <- khat(as.points(cardiff), cardiff$poly, s, newstyle=TRUE)
str(newstyle)
newstyle
apply(newstyle$khats, 2, sum)
plot(newstyle)

カーネルスムージング

  • カーネル関数を用いた平滑化を行うには,splancsパッケージのkernel2d関数,又は,spatstatパッケージのdensity.ppp関数を使います.
  • kernel2d関数を用いる方法では,花崗岩を含む地点のサンプリングデータbodminを使います.
data(bodmin)
plot(bodmin$poly, asp=1, type="n")
image(kernel2d(as.points(bodmin), bodmin$poly, h0=2, nx=100, ny=100),
add=TRUE, col=terrain.colors(20))
 Xrange is  -5.2 9.5 
 Yrange is  -11.5 8.3 
 Doing quartic kernel
pointmap(as.points(bodmin), add=TRUE)
polymap(bodmin$poly, add=TRUE)
  • density.ppp関数を用いる方法では,cellsというサンプルデータを使います.
data(cells)
Z <- density.ppp(cells, 0.05)
plot(Z)

階層クラスタリング

  • clusterパッケージを使って,空間データに対して(凝集型)階層クラスタリングを適用します.
  • ここでは,sidsのObserved変数をつかって,階層クラスタリングを行います.
sids.clust <- agnes(sids[,1])
plot(sids.clust)
 *適宜,図をクリックしてください.最後にデンドログラムが表示されます.
sids.clust
Call:    agnes(x = sids[, 1]) 
Agglomerative coefficient:  0.9913407 
Order of objects:
 [1]   1   4  10  19  20  21  23  40  77  80  83   2   7   8  22  32  35  41  45  56  73  78  87  90
[25]  11  24  36  38  48  55  66  81  17  50  58  69  70  75   3  43  44  46  52  60  67  71  72  95
[49] 100  28  34  54  63   9  13  14  15  18  29  39  59  79  84  88  89  97   5  53  25  33  64  74
[73]   6  47  57  86  31  42  65  92  96  27  91  49  76  61  99  12  30  37  51  85  98  16  62  26
[97]  93  94  68  82
Height (summary):
  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.0000  0.0000  0.0000  0.9711  0.0000 27.7200 
Available components:
[1] "order"  "height" "ac"     "merge"  "diss"   "call"   "method" "data"  
sids.clust$height
[1]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[10]  0.000000  1.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[19]  0.000000  0.000000  0.000000  0.000000  0.000000  1.970238  0.000000  0.000000  0.000000
[28]  0.000000  0.000000  0.000000  0.000000  1.000000  0.000000  0.000000  0.000000  0.000000
[37]  0.000000  3.494361  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[46]  0.000000  0.000000  0.000000  1.000000  0.000000  0.000000  0.000000  1.266667  0.000000
[55]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000
[64]  0.000000  0.000000  8.574713  0.000000  1.000000  0.000000  0.000000  0.000000  2.111111
[73]  0.000000  0.000000  0.000000  1.000000  0.000000  0.000000  0.000000  0.000000  3.600000
[82]  0.000000  1.500000  0.000000  1.000000  0.000000  6.571429  0.000000  0.000000  1.333333
[91]  1.000000  0.000000  2.666667  0.000000 27.715789  7.000000  2.000000 13.333333  6.000000
  • クラスターのグループ数に応じてAICを計算してみましょう.

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