notebook.community

Edit and run

第1講 空間的相関の検定

0 概要

maptoolsパッケージに含まれる米国ノースカロライナ州の乳幼児突然死症候群(SIDS)のデータを用いて分析を行う。

必要なライブラリーとデータの読み込み。





In [1]:



    


library(sp)
library(maptools)
library(spdep)
pj <- CRS("+proj=longlat +ellps=clrk66")
f   <- system.file("shapes/sids.shp", package = "maptools")
nc <- readShapePoly(f, IDvar ="FIPSNO", proj4string = pj)



    


















    






Checking rgeos availability: TRUE

Attaching package: ‘maptools’

The following object is masked from ‘package:sp’:

    nowrapSpatialLines

Loading required package: Matrix

0.1 空間データの扱い方

複数のファイルから構成されるシェープファイルを取り扱う。Rに取り込む方法はいくつかあるが、ここでは主にmaptoolsパッケージのreadShapeSpatial()を使う。
Rのパッケージにはサンプルデータが付属している場合が多く、今回もそのサンプルデータを使う。
しかし、自分の研究の際には自分で作らなければならず、その際にはGISソフトが役に立つ（特にArcGISは空間統計のコマンドも豊富であり、そもそもRに取り込む必要すらない）。

0.2 投影法と座標系

地図には投影法と座標系が指定されている。これらを正しく設定しないと正確な結果が得られない。
詳細は各自で学んでもらいたい（社会基盤学科の空間情報学Ⅰを履修するとよい）。
例えば、緯度経度で表現されている地図とメートルで定義されている地図を混同したり、日本なのにアメリカの座標系を指定したりする間違いは致命的である。

0.3 座標系の変換　　

CRSクラスで座標系を指定する。
例えば、

pj <- CRS("+proj=longlat +ellps=clrk66")

の場合、投影法を投影法を地球座標系に、準拠楕円体をClark1866にしている　　

1 データの準備と空間データ

データはSpatialPolygonsDataFrameクラスで与えられている。





In [2]:



    


plot(nc)

coordinates()でポリゴンの幾何中心を求めてプロットする。

asp = 1は縦横比を1に揃える引数である。
xlab = "Longitude", ylab = "Latitude"はそれぞれx軸、y軸のラベルを意味している。





In [37]:



    


plot(coordinates(nc), asp = 1, xlab = "Longitude", ylab = "Latitude")

2 面域データの自己相関分析

2.1 空間重み行列

2.2 空間自己相関指標

空間的自己相関を検定する。
a 大域的自己相関
・MoranのI統計量
・Gearyのc統計量
b　局所的自己相関
・ローカルMoran





In [10]:



    


nc.nb <- poly2nb(nc)



    











In [11]:



    


nc.nb



    


















    
Out[11]:








Neighbour list object:
Number of regions: 100 
Number of nonzero links: 490 
Percentage nonzero weights: 4.9 
Average number of links: 4.9 

















In [12]:



    


nb2listw(nc.nb, style = "B")



    


















    
Out[12]:








Characteristics of weights list object:
Neighbour list object:
Number of regions: 100 
Number of nonzero links: 490 
Percentage nonzero weights: 4.9 
Average number of links: 4.9 

Weights style: B 
Weights constants summary:
    n    nn  S0  S1    S2
B 100 10000 490 980 10696

















In [14]:



    


nc.c <- coordinates(nc)
nc.tri.nb <- tri2nb(nc.c)



    











In [16]:



    


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

MoranのI統計量

$$ I = \frac{n}{\sum^n_{i=1}\sum^n_{j=1}\omega_{ij}}\frac{\sum^n_{i=1}\sum^n_{j=1}(x_i-\bar{x})}{\sum^n_{i=1}(x_i-\bar{x})^2} $$


・正の相関がある場合は正の値を、負の相関がある場合は負の値をとる。
・相関がない場合は$\frac{-1}{(n-1}$に近づく。





In [17]:



    


library(spdep)
s <- nc@data$SID74/nc@data$BIR74 * 1e+05
moran.test(s, nb2listw(nc.knn.nb))



    


















    
Out[17]:








	Moran I test under randomisation

data:  s  
weights: nb2listw(nc.knn.nb)  

Moran I statistic standard deviate = 3.0285, p-value = 0.001229
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.217763140      -0.010101010       0.005660983

Gearyのc統計量

$$ c = \frac{(n-1)\sum_{i}\sum_{j}\omega_{ij}(z_i-z_j)^2}{2(\sum + i\sum_{j})\sum_{k}(z_i-z_j)^2} $$


・正の相関がある場合は正の値を、負の相関がある場合は負の値をとる。
・相関がない場合は$\frac{-1}{(n-1}$に近づく。





In [19]:



    


geary.test(s, nb2listw(nc.nb, style = "B"))



    


















    
Out[19]:








	Geary C test under randomisation

data:  s 
weights: nb2listw(nc.nb, style = "B") 

Geary C statistic standard deviate = 3.0989, p-value = 0.0009711
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic       Expectation          Variance 
       0.67796679        1.00000000        0.01079878 


















In [35]:



    


Imo <- localmoran(s, nb2listw(nc.nb, style = "B"))
head(Imo)



    


















    
Out[35]:









IiE.IiVar.IiZ.IiPr(z > 0)


	
37001-0.38296680-0.06060606 5.66778683-0.13540522 0.55385425

	
37003 2.84892433-0.04040404 3.77639636 1.48681905 0.06853130

	
37005 1.98692859-0.03030303 2.83039504 1.19903605 0.11525696

	
37007-3.68378017-0.04040404 3.77639636-1.87484437 0.96959293

	
37009 1.89322430-0.03030303 2.83039504 1.14333852 0.12644903

	
37011 3.20602844-0.05050505 4.72219363 1.49859284 0.06698965





















In [27]:



    


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

参考：nb2listw()

nb2listw

Details

Starting from a binary neighbours list, in which regions are either listed as neighbours or are absent (thus not in the set of neighbours for some definition), the function adds a weights list with values given by the coding scheme style chosen. B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme proposed by Tiefelsdorf et al. 1999, p. 167-168 (sums over all links to n).

補題:ボロノイ図、ドロネー図





In [ ]:

Content source: NlGG/SpatialStatistics

Similar notebooks:

notebook.community | gallery | about