空間モデリング

演習の目的

  • 空間データの分類方法のうち、local Moran's I、多次元尺度法(MDS)、階層クラスター分析、非階層クラスター分析について、統計言語Rを使って演習を行う。
  • 今回は、Rのサンプルデータを使って演習します。

Local Moran's I

  • ライブラリspdepを読み込みます。(先週インストール済み)
    library(spdep)
  • Anselin (1995) のアフリカ紛争データを使ってLocal Moran's Iを計算してみます。データやコマンドの詳細は、CRAN spdepを参照してください。
    • afcon: 42カ国の紛争データ
    • africa.rook.nb: rookタイプの隣接行列
    • paper.nb: Anselin (1995) で使われた隣接行列
  • データを読み込み、隣接関係を表示する。
    > data(afcon) 
    > plot(africa.rook.nb, afxy) 
    > plot(diffnb(paper.nb, africa.rook.nb), afxy, col="red", add=TRUE) 
    Neighbour difference for region id: MO in relation to id: MR 
    Neighbour difference for region id: MR in relation to id: MO 
    Neighbour difference for region id: CG in relation to id: TZ 
    Neighbour difference for region id: TZ in relation to id: CG 
    Neighbour difference for region id: AO in relation to id: SF 
    Neighbour difference for region id: ZA in relation to id: BC SF 
    Neighbour difference for region id: BC in relation to id: ZA 
    Neighbour difference for region id: SF in relation to id: AO ZA 
    > text(afxy, labels=attr(africa.rook.nb, "region.id"), pos=4, offset=0.4) 
  • Moran'sI Geary's Cをそれぞれ計算する
    > moran.test(afcon$totcon, nb2listw(africa.rook.nb)) 
    	Moran's I test under randomisation
    data:  afcon$totcon  
    weights: nb2listw(africa.rook.nb)    
    Moran I statistic standard deviate = 4.1251, p-value = 1.853e-05
    alternative hypothesis: greater 
    sample estimates:
    Moran I statistic       Expectation          Variance 
           0.40981723       -0.02439024        0.01107950 
    > moran.test(afcon$totcon, nb2listw(paper.nb)) 
    	Moran's I test under randomisation
    data:  afcon$totcon  
    weights: nb2listw(paper.nb)   
    Moran I statistic standard deviate = 4.3485, p-value = 6.854e-06
    alternative hypothesis: greater 
    sample estimates:
    Moran I statistic       Expectation          Variance 
           0.41679563       -0.02439024        0.01029358 
    > geary.test(afcon$totcon, nb2listw(paper.nb))
    	Geary's C test under randomisation 
    data:  afcon$totcon 
    weights: nb2listw(paper.nb)  
    Geary C statistic standard deviate = -2.8988, p-value = 0.001873
    alternative hypothesis: less 
    sample estimates:
    Geary C statistic       Expectation          Variance 
           0.58395772        1.00000000        0.02059931 
    > plot(paper.nb, afxy)
  • Moran'I値よりGeary'C値が大きいことがわかります。局地的な自己相関が強いと考えられます。
  • そこで、Local Moran's Iを計算してみましょう。
    resI <- localmoran(afcon$totcon, nb2listw(paper.nb))
    resI
                   Ii        E.Ii     Var.Ii        Z.Ii    Pr(z > 0)
    [1,]  0.005387339 -0.02439024 0.42611302  0.04561703 4.818077e-01
    [2,] -0.010036897 -0.02439024 0.12966843  0.03985989 4.841024e-01
    [3,] -0.096960693 -0.02439024 0.42611302 -0.11117251 5.442602e-01
     :
  • 紛争数とLocal Moran's Iをプロットしましょう。ここで、紛争数は標準偏差と平均値を使って標準化しています。
    > afcon.lm <- data.frame(cbind(resI[,1],(afcon$totcon- 
    mean(afcon$totcon))/sd(afcon$totcon)),row.names=afcon$name)
    > colnames(afcon.lm) <- c("Ii","standardized totcon")
    > afcon.lm
                                       Ii standardized totcon
    TUNISIA                   0.005387339          0.01100754
    ALGERIA                  -0.010036897          0.06257364
    MOROCCO                  -0.096960693          0.45376473
     :
    > plot(afcon.lm,xlab="Local Moran's I", ylab="Standardized Total Conflicts")
    > text(afcon.lm,rownames(afcon.lm))
  • テキストを表示すると判読しにくい(できない)ので、エクセルか何かにデータをエクスポートしてグラフ表示するとよいでしょう。

多次元尺度法

  • ここでは、datasetsというパッケージにあるヨーロッパの都市間距離eurodistデータに、多次元尺度法を適用してみましょう。
    > library(datasets)
    > library(stats)
    > eurodist
                    Athens Barcelona Brussels Calais Cherbourg Cologne
    Barcelona         3313                                            
    Brussels          2963      1318       
     :
    > eur.cmd <- cmdscale(eurodist,k=2)
    > plot(eur.cmd)
    > text(eur.cmd,rownames(eur.cmd ),cex=0.6)

階層クラスター分析

  • ここで使うデータは,ノースカロライナ州における乳幼児突然死症候群(SIDS)のデータ(nc.sids)です.DClusterパッケージで利用可能です.
  • nc.sids$SID74は'74年の乳幼児突然死症候群にかかった新生児数,nc.sids$BIR74は同年の新生児数です.Expectedは各都市での期待死亡数を意味します.
    > library(DCluster)
    > 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)
  • clusterパッケージを使って,空間データに対して(凝集型)階層クラスタリングを適用します.
  • ここでは,sidsのObserved変数をつかって,階層クラスタリングを行います.
    >  library(cluster)
    >  sids.clust <- agnes(sids[,1])
    >  plot(sids.clust)
     次の図を見るためには<Return>キーを押して下さい:  
     次の図を見るためには<Return>キーを押して下さい:  
     *適宜,図をクリックしてください.最後にデンドログラムが表示されます.
    > 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

非階層クラスター分析

  • K-平均法は、kmeans()関数を使って計算できます。
  • kmeans(データ, クラスタ数)
    > sids.km<-kmeans(sids[,1],5)
    > plot(sids[,1],col=sids.km$cluster)
    > points(sids.km$centers,col=1:5,pch=8,cex=2)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-05-24 (月) 16:48:18 (3345d)