読者です 読者をやめる 読者になる 読者になる

ill-identified diary

所属組織の業務や見解とは一切無関係なアフィリエイト付きメモ帳。所属とは関係ないけどここを見て所属先にも興味を持っていただけると喜びます。

[QGIS] [R] QGIS と空間統計モデル (CARモデル)

GIS ggplot2 R 計量経済学

概要・前置き

  • 以前も何度か R で地図を作る方法を紹介していたが, 自分のプログラミングテクが雑なこともあり, 冗長なコードの掲載であまり便利でないのではという印象を持たれる恐れもあった. そこで, GUI で操作のできるわりに高機能な QGIS (Quantum GIS; QGISプロジェクトへようこそ!) の使い方についても紹介したいと考えていた.

  • すると先日, 岩波DSが発売され, CAR モデルが取り上げられていたので, これと絡めて書くことにした.

  • ただ, 実際は QGIS で直接 CAR の推定はできないので, 今回はあまり QGIS の出番はなく, 空間統計モデルの説明と R へのつなげ方がメインになってしまった…… 次回以降にデータハンドリングの場面での QGIS の強みを紹介したい……

  • 分量は PDF 換算 6ページ.

空間計量経済モデル

主題である CAR ( Conditional AutoRegressive; 条件付き自己回帰) モデルに触れる前に, まず一般的な空間自己相関モデルについて話しておく. 株価などの金融時系列データならば, 過去の値を説明変数にして現在の値に当てはめ, 予測する時系列的な自己回帰モデルがよく使われる. つまり,  \{y_{1},\,y_{2},\,\cdots,\,y_{T}\} という時系列データがあって,

\begin{align}
y_{t}= & f(y_{t-j}),\,t=0,\,1,\,2,\,3,\, \cdots ,\, T \end{align}

というような関係を意味する. 一方で, 自身の過去の値ではなく, 自身の周囲の周りの値によって, 自身の値を決定する「空間的な」自己回帰モデルというのもある. 時間以外のインデックス  i のついた  y_{1},\,\cdots,\,y_{i},\,\cdots,\,y_{N} があって, やはり

\begin{align}
y_{i}= & f(y_{j}),\,j\neq i\end{align}

のように表される. 地理経済学の他にも, 伊東 (2015) の解説のように生態の調査でも用いられているようで, Lesage (1999) で挙げられている例も, ミシガン州マイマイガの個体の分布である. マイマイガの幼虫は森林害虫であるため, 生態について昔から研究がなされているようだ. マイマイガの幼虫は, 飛ぶことはできないが風に吹かれることで移動する. つまり, まったくのランダムに自然に出現するのではなく, 最初にいた地点からわずかに近くに多く分布しやすい, という空間的な自己相関があると考えられる. ではこれを, 古典的な線形回帰モデルで表現できるだろうか?

空間統計モデリングでは空間的な関係を表す隣接行列 (contiguity matrices) を得ることが重要になる. そこで, まずは簡単な例として, ミシガン州を9の区画に分けてマイマイガの個体数を調査したとする. 区画は碁盤の目のように分割され, 区画ごとに調査がなされたとしよう. 区画  i の個体数を  y_{i} とする. 空間的な自己相関は, 隣接する区画にのみ作用し, 互いに隣接しない区画どうしでは相関がないものと仮定する*1. このとき, 全部で9区画の場合の中央の区画5は左右と上下の区画と隣接しているとみなせる. 他の区画も同様に上下左右と隣接しているので,以下のように各区画に順に番号を割り当てていくことで,

1 2 3
4 5 6
7 8 9


この各観測点ごとの隣接の情報を行列で表すことができる: 区画  i が区画  j に隣接しているとき,  (i,\,j) 成分が  1 , 隣接していないならゼロとなる ( i=j なら隣接ではないのでゼロ) ようにすると,

\begin{align}
\mathbf{W}= & \begin{bmatrix}0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\\
0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 0\\
 &  &  &  &  & \ddots\\
\\
\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0
\end{bmatrix}\end{align}

と表せる. 簡単のために碁盤の目状と言ったが, 都道府県や市区町村のようなより複雑な区域でも隣接行列を決定することはできる. また, 隣接の考え方にもいくつか種類があり, 上の例はチェスの駒にたとえて「ルーク型」である. 斜め方向を隣接というのなら「ビショップ型」, 上下左右と斜めを隣接していると定めるのなら「クイーン型」である. ただしどのタイプであれ, 定義から,  \mathbf{W} の対角成分は必ずゼロでなければならない. さらに, 重み付けに使う必要から, この  \mathbf{W} の各行の和が1になるように「標準化」する必要がある. よって, 上記のルーク型隣接行列は,

\begin{align}
\mathbf{W}= & \begin{bmatrix}0 & 1/2 & 0 & 1/2 & 0 & 0 & 0 & 0 & 0\\
1/3 & 0 & 1/3 & 0 & 1/3 & 0 & 0 & 0 & 0\\
0 & 1/3 & 0 & 1/3 & 0 & 1/3 & 0 & 0 & 0\\
1/3 & 0 & 0 & 0 & 1/3 & 0 & 0 & 1/3 & 0\\
0 & 1/4 & 0 & 1/4 & 0 & 1/4 & 0 & 1/4 & 0\\
 &  &  &  &  & \ddots\\
\\
\\
\\
\end{bmatrix}\end{align}

となる. この行列を掛けて,

\begin{align}
\begin{bmatrix}y_{1}\\
\vdots\\
y_{N}
\end{bmatrix}= & \rho\mathbf{W}\begin{bmatrix}y_{1}\\
\vdots\\
y_{N}
\end{bmatrix}+\begin{bmatrix}\varepsilon_{1}\\
\vdots\\
\varepsilon_{N}
\end{bmatrix}\end{align}

と表す. よって,  y_{i}^{\ast}=(\mathbf{W}\boldsymbol{y})_{i} とみなせば,  y_{i}^{\ast}  y_{i} に隣接する区画での観測値の加重平均とみなせるので,

\begin{align}
y_{i}= & \rho y_{i}^{\ast}+\varepsilon_{i}\label{eq:far}\end{align}

という通常の線形回帰モデルのようにあつかうことで, 隣接区画との自己相関パラメータ  \rho を推定できる1階の空間自己回帰 (FAR) モデル, あるいは時系列の類推から空間ラグモデルと呼ばれるものになる*2. これは定数項や他の説明変数を加えり, 誤差項にも自己相関を仮定したモデルへも拡張できて, それらは空間ダービン (Spatial Durbin) モデルとか, 空間自己回帰 (SAR) モデルとか呼ばれるものになる*3.

空間統計では, 単に自己相関といっても, 観測点どうしの距離が開くほど自己相関が減衰するという場合が多い. これは横軸に距離を, 縦軸を自己相関の大きさとしてプロットするとよくわかる. なお, このような 距離と自己共分散の関係はバリオグラムという指標で表し, 空間統計でよく使われる*4. ここから, 隣接行列の代わりに より一般に, 成分に観測点どうしの距離に依存する関数を与えたウエイト行列を用いる場合があり, 地理的加重回帰 (GWR) モデルと呼ばれる.

ここまで紹介したモデルは, GWR モデルを除き, 古典的な計量経済学の枠組みでパラメータを適切に推定できる. また, R の spdep パッケージで実際の計算ができる. なお, Lesage (1999) には空間計量経済学のモデル全般について説明が書かれているが, CAR については言及されておらず, またサンプルコードが多く掲載されているものの全て matlab である.

CAR モデル

SAR モデル全般は,  \boldsymbol{\varepsilon}\sim\mathcal{N}(0,\,\sigma^{2}\mathbf{I}_{N}) , というふうな正規分布か, さらに自己回帰項を合成して表されるという仮定だった. 一方で, CAR モデルの場合, さらに地点ごとに個別の誤差項を持つ分布であるとする. 言い換えると, 空間モデルに混合効果モデルを持ち込んだものとも言える. つまり, たとえば正規分布ならば, 分散共分散行列の対角成分以外にも非ゼロの成分が存在するため, 一般的な行列  \mathbf{V} を用いて表される  \boldsymbol{\varepsilon}\sim\mathcal{N}(\mathbf{0},\,\sigma^{2}\mathbf{V}) のような多変量正規分布となることを意味する. 伊東 (2015) などでは, CAR モデル, 特に intrinsic CAR, または Gaussian CAR と呼ばれるものの誤差項は

\begin{align}
\varepsilon_{i}|\varepsilon_{j\neq i}= & \mathcal{N}\left(\sum_{j}\frac{w_{i,j}\varepsilon_{j}}{\sum_{j}w_{i,j}},\,\frac{\sigma_{i}^{2}}{\sum_{j}w_{i,j}}\right)\end{align}

というふうに表される (前後の記述と整合をとるため, 記号を変えている) が, そもそもこれはどういう意味かというのも書いておく. Besag, York, and Mollié (1991) では, 誤差項の事前分布を \boldsymbol{\varepsilon} の同時分布として以下のように表している.

\begin{align}
p(\boldsymbol{\varepsilon}) & \propto\exp\left\{ -\sum_{i < j}w_{ij}\phi(\varepsilon_{i}-\varepsilon_{j})\right\} \end{align}

 \phi(z) は任意の |z| について増加する関数である.  w_{ij} は隣接行列の成分であるが, 先に紹介したものとは異なり, 隣接していれば 1, そうでなければゼロとなり, 列の和で除すことで標準化しない. 標準化しないのは最初の分布で表された完成形を見手の通りである. これは全ての  \{\varepsilon_{i}\} の同時分布だが,  \varepsilon_{i} 以外の残りの誤差項  \varepsilon_{-i} に対する \varepsilon_{i} の条件付き密度関数に書き換えると,

\begin{align}
p(\varepsilon_{i}|\varepsilon_{-i})\propto & \exp\left\{ -\sum_{j\in\delta_{i}}w_{ij}\phi(\varepsilon_{i}-\varepsilon_{j})\right\} \end{align}

となる. ただし, このとき,  \delta_{i}  i に隣接しているすべての地点という意味である.

さらに,  \phi を正数のパラメータ  \kappa を伴う  \phi(z)=z^{2}/2\kappa とおくと, 当初の正規分布の密度関数となる. これが intrinsic CAR である. 一方で, "intrinsic" でない CAR モデルとして, リモートセンシングなどでは  \phi(z)=|z|/\kappa がよく使われるらしい.

QGIS での操作

結論から言うと, 空間自己相関モデル全般を QGIS 単体で計算するのは難しい. しかし, 地理データのハンドリングに関しては QGIS は便利なので, モデルパラメータの推定の全段階までは QGIS でデータハンドリングを行い, それを R で読み込んで残りの工程を行う, という方針にする.

一番面倒なのは隣接行列の計算だろう. そこで, maps2WinBUGS というプラグインで WinBUGS フォーマットで隣接行列を出力し, R の CARBayes パッケージで CAR モデルを推定する, という予定だったのだが……, WinBUGS のグラフ機能用に地理データを変換する機能など, 今はあまり重要でない機能が多いので, 今回はこれを使わないことにした.


f:id:ill-identified:20151120002817p:plain
隣接行列の出力はできるが, R でやるよりも遅い……

QGIS で加工した地理データを R で読み込み, spdepCARBayes で CAR モデルを推定することにする. WinBUGSWindows 以外の OS での導入が面倒なので, より OS に依存しない CARBayes を利用することにした. もしも WinBUGS を利用したいというのなら, 伊東 (2015) を参考にすれば簡単にできるだろう. spdep は隣接行列を計算するのに必要である. 隣接行列を求めるだけなら rgeos パッケージでもいいが, spdep ならば SAR など他の空間自己相関モデルの推定ができるパッケージなのでこちらを採用した.

今回は例として, 国土数値情報ダウンロードサービス から東京都内の最新の公示地価データをダウンロードし, CAR モデルを適用してみた. 点データなので, ボロノイ分割してポリゴンに置き換える. 今回は, ボロノイポリゴンの隣接から隣接行列を求めることにした.

  1. 点データのレイヤを選択.

  2. [ベクタ] -> [ジオメトリツール] -> [ボロノイポリゴン] を選択.

  3. 入力レイヤと出力ファイル名を選択し, バッファに 5% を指定して実行.

この処理は重いのでノート PC など, スペックが低いマシンを使っている場合は注意.


f:id:ill-identified:20151120001618p:plain
ポリゴン分割し単価で色分けした様子 (端は東京都の境界線でクリップしている)

出力されたポリゴンには, 変換前の点に対応した属性情報が紐付けられているので, これを R に読み込ませる. 以前 [R] Rで学ぶ都知事選のデータ可視化【地理データ編】 - ill-identified diary で書いたように, maptools でR に読み込ませたら, まず隣接行列を計算させる. すぐに隣接関係を出力できる. ポリゴンの隣接関係は spdep::poly2nb() で計算できる. 点と点の距離で判定したいならば, nb2WB() などを使う (rgeos なら gTouches() を使う. ). poly2nb() は隣接情報をリストで返すので, これを matrix 型に変換して, CARBayes::S.CARia() に渡す. CAR モデルの場合, 標準化は必要ないので, 隣接行列の成分  w_{ij} 必ずゼロか1にする*5.

隣接行列を作成したら, あとは CARBayes::S.CARiar() 関数でフィッティングを行う. 空間自己相関モデルで都心の地価を推定した研究に, 唐渡 (2007) というのがあり, こちらはもっといろいろな説明変数を用意しているが, 今回はとりあえず実際のデータを持ってきて動くところを見せられれば, という考えなので, 説明変数は適当である (size は面積である.).

require(CARBayes) # 4.3
require(maptools) # 0.8-30
require(spdep)    # 0.5-88
require(ggplot2)  # 1.0.1
require(ggthemes) # 2.2.1
setwd("~/Documents/blog/20151027_GIS/") # 任意のディレクトリに変更
voronoi <- readShapePoly("out.shp")
# list 型から matrix へ
nbl2mat <- function(x, is.dev=F){
  w <- matrix(0, ncol=length(x), nrow=length(x))
  for (i in 1:length(x)){
    w[i,x[[i]]] <- 1
    if ( is.dev == T) w[i,] <- w[i,]/sum(w[i,])
  }
  return(w)
}
W <- nbl2mat(poly2nb(pl = voronoi, queen = F), is.dev = F)
fit <- S.CARiar(price ~ price + size, data=voronoi@data, family = "gaussian", W=W, burnin = 2000, n.sample = 32000, thin=10)
voronoi$predict <- fit$fitted.values
voronoi$id <- rownames(voronoi@data)
df <- fortify(voronoi)
df <- merge(x=df, y=voronoi[,c("id", "price","predict","region","region_2")], by="id")
df$region_2 <- as.character(df$region_2)
ggplot(df) + geom_polygon(aes(x=long, y=lat, group=group, fill=predict) ) + scale_fill_continuous_tableau() + coord_equal() 
ggplot(df) + geom_polygon(aes(x=long, y=lat, group=group, fill=price) ) + scale_fill_continuous_tableau() + coord_equal()
ggplot(df[grepl("区$",df$region_2),]) + geom_polygon(aes(x=long, y=lat, group=group, fill=predict ) ) + scale_fill_continuous_tableau() + coord_equal() 
ggplot(df[grepl("区$",df$region_2),]) + geom_polygon(aes(x=long, y=lat, group=group, fill=price ) ) + scale_fill_continuous_tableau() + coord_equal()

で, 適当な変数にしたからか全くフィットしていないようだ……


f:id:ill-identified:20151120002047p:plain
f:id:ill-identified:20151120002103p:plain
上:推定結果, 下:実際の値


f:id:ill-identified:20151120002145p:plain
f:id:ill-identified:20151120002153p:plain
上:推定結果 (23区),下:実際の値 (23区)

今回は QGIS の機能をほとんど紹介できなかったが, 地理情報のデータハンドリングに関しては便利な機能が多く存在するので, また別の機会にまとめておきたい.

参考文献


*1:ミシガン州は正方形ではないが, 話の簡単のためそう想定する.

*2:ガウシアン AR モデルと同様に, 最小二乗法ではなく最尤法で推定しなければならない.

*3:Lesage (1999) などを参照. ただし, 略称が同じである同時自己回帰モデル (SAR) はこれらの空間自己回帰モデル全般を指すようだ(Kissling and Carl 2008).

*4:どちらかというと地球統計学? 間瀬 (2010) などを参考

*5:先に書いた CAR モデルの分布の段階で標準化がなされるからである.