# 点パターン分析・・・地点データの集積性の分析

# 完全空間ランダム分布(CSR)

# 地点分布と完全空間ランダム分布を比較することで地点分布の集積性を検定

# CSR分布に比べて地点分布がどの程度偏っているかを測る手法として

# 区画法や最近距離法などの方法が提案されている

# cardiff(カーディフ市に住む168人の少年犯罪者の居住地データ)の例

# 少年犯罪者の居住地の地図

library(splancs)

data(cardiff)

polymap(cardiff$poly,xlab="Easting",ylab="Northing")

points(cardiff$x,cardiff$y,pch=19)

title(main="少年犯罪者の居住地")

# CSR分布での地図

polymap(cardiff$poly,xlab="Easting",ylab="Northing")

points(csr(cardiff$poly,length(cardiff$x)),pch=1)

title(main="CSR分布")

# 面域データの空間自己相関分析

# 空間重み行列・・・隣接または近接する区画の関係からその影響の大きさを定義する行列

# 2進数的重み係数:要素が0,1の空間重み行列, 必要に応じて標準化(区域ごと or 全区域)

# 一般化重み係数:連続値で表される重み係数

# 例)区域の共通な境界線の長さと区域間距離を用いた空間的重み係数

# ただし,R単独では一般化重み係数は計算でず,他のGISソフトとの連携が必要

# Rでは,区域同士が接続しているか否かの情報を含む近隣リストを作成し,空間重み行列を作成

# 米国ノースカロライナ州SIDS死亡データの例

# まず,SpatialPolygonsDataFrameクラスのncを作成

library(maptools)

pj <- CRS("+proj=longlat +sllps=clrk66") # 投影法=地理座標系,準拠楕円体=Clarke1866

f <- system.file("shapes/sids.shp",package="maptools") # シェープファイルのフルパス

nc <- readShapePoly(f,IDvar="FIPSNO",proj4string=pj) # シェープファイル読込

# poly2nb()でSpatialPolygonsDataFrameなどからnbクラスを出力

library(spdep)

nc.nb <- poly2nb(nc)

# 近隣リストのクラスとしてnbクラスがあり,nb2listw()でnbクラスから空間重み行列を生成

# nb2listw()の出力は行列ではない,行列で出力したい場合はnb2mat()がある

# B:標準化なし,W:区域ごとの標準化,C:全区域で標準化

# U:全区域で標準化し区域数で割る,S:Tiefelsdorfら(1999)による分散安定化方式

nb2listw(nc.nb,style="B") # 標準化なしの2進数的重み

# ドロネー三角網を利用した隣接リスト

# 最近隣の2点を結んで作成する三角形のネットワーク

nc.c <- coordinates(nc) # 面域の幾何中心を出力

nc.tri.nb <- tri2nb(nc.c) # 幾何中心に基づくドロネー三角網に基づく近隣リストを出力

# 隣接区域を1次隣接とし,隣接区域に隣接する区域を2次隣接を考える(一般にはk次隣接区域)

# knearneigh()でk次隣接区域の近隣関係を計算する

nc.knn.nb <- knn2nb(knearneigh(nc.c,k=3))

# 一定距離以内の地点に隣接関係があると定義するには

# dnearneigh()を使う,最小距離(d1)と最大距離(d2)

nc.dn.nb <- dnearneigh(nc.c,d1=0,d2=40,longlat=T) # 半径40km以内の円に含む区域を近隣区域とする

# 作成したnbクラスの視覚化

par(mfrow=c(2,2))

plot(nc,border="grey"); plot(nc.nb,nc.c,col="red",add=T); title(main="隣接の有無")

plot(nc,border="grey"); plot(nc.tri.nb,nc.c,col="red",add=T); title(main="ドロネー三角網")

plot(nc,border="grey"); plot(nc.knn.nb,nc.c,col="red",add=T); title(main="k次隣接(k=3)")

plot(nc,border="grey"); plot(nc.dn.nb,nc.c,col="red",add=T); title(main="一定距離以内(40km)")

# diffnb()を使ってnbクラスの差分がわかる

diffs <- diffnb(nc.nb,nc.knn.nb)

plot(nc,border="grey"); plot(diffs,nc.c,col="red",add=T); title(main="隣接の有無と3次隣接の差分")

# 空間自己相関指標

# MoranのI統計量・・・空間自己共分散を標準化した統計量,相関係数と同様な解釈

# Gearyのc統計量・・・小さいc統計量は高い空間相関を意味する

# ローカル検定・・・局所的に空間自己相関を検討する検定,区別して上記2つをグローバル検定と呼ぶ

# sids死亡データの例(空間自己相関指標)

s <- nc@data$SID74/nc@data$BIR74 * 1e+5 # 群別の10万人対SIDS死亡率

moran.test(s,nb2listw(nc.knn.nb)) # MoranのI統計量

geary.test(s,nb2listw(nc.nb,style="B")) # Gearyのc統計量

lmo <- localmoran(s,nb2listw(nc.nb,style="B")) # ローカルMoranのI統計量

# spdepパッケージのmoran.plot()を使って,

# 地区における観測値と近隣地区の重みつき観測値の散布図を描く

moran.plot(s,nb2listw(nc.nb,style="B"))

# 地域集積の検出と検定

# 空間スキャン統計量・・・移動窓(moving window)をくまなく走査(scan)することで集積性を分析

# フォーカスド検定・・・ある特定の地点に近づくと罹患が増えるかどうかを検定

# 空間スキャン統計量

# 対象地域にグリッドを作成し,グリッド点を中心として円を設定して計算

#(ポリゴンとその中心点といった不規則な配置で利用される場合が多い)

# 円の半径を変更しながら,円の内外での症例数・人口を比較する

# DClusterパッケージのkullnagar.stat()はKulldorff-Nagarwallaの空間スキャン統計量を計算する

# sids死亡データの例(Kulldorff-Nagarwallaの空間スキャン統計量)

library(DCluster)

Observed <- nc@data$SID74

Expected <- nc@data$BIR74*sum(nc@data$SID74)/sum(nc@data$BIR74)

Population <- nc@data$BIR74

dat <- data.frame(Observed,Expected,Population)

summary(dat)

xy <- coordinates(nc) # 各群の中心点

d <- spDistsN1(xy,xy[1,],longlat=T) # 1番目の群からの距離を計算

i <- order(d) # 距離の近い順に並べる

kullnagar.stat(dat[i,],fractpop=0.5) # 空間スキャン統計量を計算