第09回演習マニュアル

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

データの準備

  • 最初に使うデータは,ノースカロライナ州における乳幼児突然死症候群(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))

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

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

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

GAM: Geographical Analysis Machine (Openshaw's clustering)

  • opgamというコマンドを使います.
  • opgamコマンドに,データセット,検索半径(radius),グリッドの大きさ(step),有意水準(alpha)を設定します.
  • GAMの検定では,有意水準(alpha)をかなり小さめの値(例えば,p<0.002)に設定することが多い.
# 上で作成したsidsにXY座標を追加
sids<-cbind(sids, x=nc.sids$x, y=nc.sids$y)
# opgam(data, radius, step, alpha)
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)
    plot(pp)
    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)

カーネルスムージング1

  • splancsパッケージのkernel2d()関数を使って、カーネル密度関数をグリッドにスムージングし、表示しましょう。
  • データは、上述の青少年犯罪データcardiffを用います。
    plot(cardiff$poly, asp=1,type="n")
    image(kernel2d(as.points(cardiff),cardiff$poly,h0=2,nx=100,ny=100),
    add=TRUE,col=terrain.colors(20))
    pointmap(as.points(cardiff),add=TRUE)
  • 次に、カーネル関数のバンド幅bandwidthを変えて、密度をプロットしてみましょう。
    plot(cardiff$poly, asp=1,type="n")
    image(kernel2d(as.points(cardiff),cardiff$poly,h0=3,nx=100,ny=100),
    add=TRUE,col=terrain.colors(20))
    pointmap(as.points(cardiff),add=TRUE)

カーネルスムージング2

  • カーネル関数を用いた平滑化は,splancsパッケージのkernel2d()関数やspkernel2d()関数,又は,spatstatパッケージのdensity.ppp関数を使っても可能です.
  • 花崗岩を含む地点のサンプリングデータbodminを使い、kernel2d関数を用いる方法では表示してみましょう.
    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))
    pointmap(as.points(bodmin), add=TRUE)
    polymap(bodmin$poly, add=TRUE)
  • spkernel2d()関数を使う場合は、あらかじめ表示するグリッドを定義します。
    bodmin.xy <- coordinates(bodmin[1:2])
    apply(bodmin$poly, 2, range)
    grd <- GridTopology(cellcentre.offset=c(-5.2, -11.5), cellsize=c(0.2, 0.2),  
    cells.dim=c(75,100))
    k10 <- spkernel2d(bodmin.xy, bodmin$poly, h0=1, grd)
    k15 <- spkernel2d(bodmin.xy, bodmin$poly, h0=1.5, grd)
    k20 <- spkernel2d(bodmin.xy, bodmin$poly, h0=2, grd)
    k25 <- spkernel2d(bodmin.xy, bodmin$poly, h0=2.5, grd)
    if (.sp_lt_0.9()) {
     df <- AttributeList(list(k100=k100, k150=k150, k200=k200, k250=k250))
    } else {
     df <- data.frame(k100=k100, k150=k150, k200=k200, k250=k250)
    }
    kernels <- SpatialGridDataFrame(grd1, data=df)
    spplot(kernels, checkEmptyRC=FALSE, col.regions=terrain.colors(16), cuts=15)
  • density.ppp()関数を用いる方法では,cellsというサンプルデータを使います.
    data(cells)
    Z <- density.ppp(cells, 0.05)
    plot(Z)

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