確率と情報の科学：データ解析のための統計モデリング入門：久保拓弥、岩波書店
第11章 空間構造のある階層ベイズモデル
???????????
場所差の空間相関(spatial correlation)を考慮する統計モデリング。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/CarNormalExample.html

11.1　例題:一次空間上の個体数分布
調査区画を50個設定し、それが一本の直線上に等間隔で配置されるとする。
調整区画の番号jは、左から1,2,….,50とする。

ここでは架空のデータ{yi}は場所jごとに平均の異なるポアソン乱数として発生させたもの。

http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/m1/model.bug.txt
model
{
for (j in 1:N.site) {
Y[j] ~ dpois(mean[j])
log(mean[j]) <- beta + r[j] } r[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau) beta ~ dnorm(0, 1.0E-4) tau <- 1 / (s * s) s ~ dunif(0, 1.0E+4) } > source(“R2WBwrapper.R”)
# reading “R2WBwrapper.R” (written by kubo@ees.hokudai.ac.jp)…
> clear.data.param()
>
> load(“Y.RData”)
> set.data(“N.site”, length(Y))
> set.data(“Y”, Y)
>
> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1))) > set.data(“Adj”, c(2, Adj, N.site – 1))
> set.data(“Weights”, rep(1, 2 * N.site – 2))
> set.data(“Num”, c(1, rep(2, N.site – 2), 1))
>
> set.param(“beta”, 0)
> set.param(“r”, rnorm(N.site, 0, 0.1))
> set.param(“s”, 1)
>
> post.bugs <- call.bugs( + n.iter = 10100, n.burnin = 100, n.thin = 10 + ) > post.list <- to.list(post.bugs) > post.mcmc <- to.mcmc(post.bugs) > plot(Y, xlim=c(0,50), ylim=c(0,25))

11.2　階層ベイズモデルに空間構造をくみこむ
まずとっかかりとして、個体数yiがすべての区間で共通する平均値λのポアソン分布に従うとする。
データは過分散であり、区画jごとに平均λjが異なっているとする。

平均値λjを線形予測子と対数リンク関数を用いて、
logλj = β + ri
とする。
この線形予測子は、切片β（大域パラメータ）と場所差rj（局所パラメータ）で構成されている。
βには、無情報事前分布を、{rj}には階層事前分布を与える。

11.2.1 空間構造のない階層事前分布

P(rj|s) = 1/√2πs^2exp(-rj^2/2s^2)

11.2.1 空間構造のある階層事前分布
以下の仮定を設ける

区間の場所差は「近傍区間の場所差にしか影響しない
区画jの近傍の個体数njは有限であり、その区画が近傍であるかはモデル設計者が指定する。
近傍の直接の影響はそれも等しく1/nj

さらに簡略化して、この一時空間モデルでは、ある区画はその隣接している区画とだけ相互作用を持つ。
したがって、近傍数njは２で、左右の両端(n1,n50)では1。

p(rj|μj,s) = √(nj/2πs^2exp(-(rj-μj)^2/2s^2/nj)

μj = (rj-1+rj+1)/2、ただし左右の両端j=1とJ=50では、 μ1 = r2, μ50 = r49

このような事前分布は、条件付き自己回帰(conditional auto regressive, CAR)モデルと呼ばれる。
このモデルのようにさまざまの制約をつけて簡単にしたものはintrinsic Gaussian CARモデルと分類される。

この場所差{rj} = {r1,r2,….r50}全体の事前分布である同時分布p){rj}|s)は以下のようになる。
p({rj}|s) ∝ exp(^1/2s^2Σ(rj-rj’)^2

11.3 空間統計モデルをデータにあてはめる

事後分布は、

p(β,s,{rj}|Y) ∝ p({rj}|s)p(s)p(β)Πp(yi|λi)

データyjが得られる確率p(yj|λj)は平均λj=exp(β+rj)のポアソン分布とする。
大域的なパラメータβには、無情報事前分布p(β)を指定する。
局所的なパラメータである場所差rjの事前分布は、空間相関を考慮した階層事前分布として、同時分布p({rj}|s)を使っている。MCMCサンプリングでは、個々のrjの条件付き事前分布p(rj|μ,s)を使用し、これは平均μjで標準偏差s/√njの正規分布、sの事前分布p(s)は無情報事前分布であるとする。

http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/m1/model.bug.txt
model
{
for (j in 1:N.site) {
Y[j] ~ dpois(mean[j])
log(mean[j]) <- beta + r[j] } r[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau) beta ~ dnorm(0, 1.0E-4) tau <- 1 / (s * s) s ~ dunif(0, 1.0E+4) } > load(“Y.RData”)
> set.data(“N.site”, length(Y))
> set.data(“Y”, Y)
>
> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1))) > set.data(“Adj”, c(2, Adj, N.site – 1))
> set.data(“Weights”, rep(1, 2 * N.site – 2))
> set.data(“Num”, c(1, rep(2, N.site – 2), 1))
>
> set.param(“beta”, 0)
> set.param(“r”, rnorm(N.site, 0, 0.1))
> set.param(“s”, 1)
>
> post.bugs <- call.bugs( + n.iter = 10100, n.burnin = 100, n.thin = 10 + ) > post.list <- to.list(post.bugs) > post.mcmc <- to.mcmc(post.bugs) > plot(post.bugs)

BUGSコードでの注意点は、
r[1:N.site]~car.normal(Adj[], Weights[], Num[], tau)
で、car.normal()関数は、intrinsic Gaussian CARモデルから[rj}のMCMCサンプルを発生させる関数。Adj[]は隣接する場所の番号j、Weights[]はAdj[]に対応する重み、Num[]はそれぞれの区画委隣接する場所の数nj、tauは分散の逆数、つまり1/s^2。

> library(R2WinBUGS)
> if (!exists(“Y”)) load(“Y.RData”)
> if (!exists(“post.bugs”)) load(“post.bugs.RData”)
> #post.mcmc <- to.mcmc(post.bugs) > post.mcmc <- as.mcmc(post.bugs$sims.matrix) > v <- 1:length(Y) > y.max <- 27 >
> plot(
+ v, Y,
+ xlab = “location”,
+ ylab = “y_i”,
+ ylim = c(0, y.max)
+ )
> lines(m, lwd = 2, lty = 2)
>
> mre <- sapply( + v, function(i) quantile( + post.mcmc[, sprintf("r[%i]", i)], + probs = c(0.5, 0.025, 0.975) + ) + ) > b <- median(post.mcmc[,"beta"]) > polygon(
+ c(v, rev(v)),
+ exp(b + c(mre[2,], rev(mre[3,]))),
+ border = NA,
+ col = “#00000030”
+ )
> lines(exp(b + mre[1,]), lwd = 2)

http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/m2/model.bug.txt

model
{
for (j in 1:N.site) {
Y[j] ~ dpois(mean[j])
log(mean[j]) <- beta + r[j] r[j] ~ dnorm(0.0, tau) } beta ~ dnorm(0, 1.0E-4) tau <- 1 / (s * s) s ~ dunif(0, 1.0E+4) } > source(“R2WBwrapper.R”)
# reading “R2WBwrapper.R” (written by kubo@ees.hokudai.ac.jp)…
> clear.data.param()
>
> load(“Y.RData”)
> set.data(“N.site”, length(Y))
> set.data(“Y”, Y)
>
> set.param(“beta”, 0)
> set.param(“r”, rnorm(N.site, 0, 0.1))
> set.param(“s”, 1)
>
> post.bugs <- call.bugs( + n.iter = 10100, n.burnin = 100, n.thin = 10 + ) > post.list <- to.list(post.bugs) > post.mcmc <- to.mcmc(post.bugs) > plot(post.bugs)

http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/model.bug.txt

model
{
# j = 1: GLMM, no NA data
# j = 2: GLMM, NA data
# j = 3: CAR, no NA data
# j = 4: CAR, NA data
Tau.noninformative <- 1.0E-4 S.max <- 1.0E+4 for (i in 1:N.site) { for (j in 1:4) { Y[i, j] ~ dpois(mean [i, j]) log(mean[i, j]) <- beta[j] + r[j, i] } for (j in 1:2) { r[j, i] ~ dnorm(0.0, tau[j]) } } r[3, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[3]) r[4, 1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau[4]) for (j in 1:4) { beta[j] ~ dnorm(0, Tau.noninformative) tau[j] <- 1 / (s[j] * s[j]) s[j] ~ dunif(0, S.max) } } > source(“R2WBwrapper.R”)
# reading “R2WBwrapper.R” (written by kubo@ees.hokudai.ac.jp)…
> clear.data.param()
>
> load(“Y.RData”)
> no <- c(6, 9, 12, 13, 26:30) > Yno <- Y > Yno[no] <- NA >
> set.data(“N.site”, length(Y))
> set.data(“Y”, matrix(c(Y, Yno, Y, Yno), N.site, 4))
>
> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1))) > set.data(“Adj”, c(2, Adj, N.site – 1))
> set.data(“Weights”, rep(1, 2 * N.site – 2))
> set.data(“Num”, c(1, rep(2, N.site – 2), 1))

11.4　空間統計モデルが作り出す確率場

隣同士の相互作用をするrjの一次元ならびである{r1,r2,…r50}の関係を決めるために、intrinsic Gaussian CARモデルを部品として用いた。
一般に相互作用をする確率変数でうめつくされた空間は確率場random fieldと呼ばれる。ばらつきパラメータsの大小が確率場{rj}にあたえる影響を見ると、大域的なパラメータβを適当に固定して、sを0.0316, 0.224, 10.0とした。ｓが小さいときは両隣の平均と似た、rj全体のばらつきは小さくなり、sが大きいと全体にぎざぎざしたばらつきの大きい確率場となる。

11.5 空間相関モデルと欠測のある観測データ
空間相関を組み込んだ階層ベイズモデルでは「欠測データ」に対して、より良い予測が得られる。

全部まとめて推定 (欠測あり・なし，CAR と GLMM)
> source(“set.data.R”)
# reading “R2WBwrapper.R” (written by kubo@ees.hokudai.ac.jp)…
>
> set.param(“beta”, rep(0, 4))
> set.param(“r”, matrix(rnorm(N.site * 4, 0, 0.1), 4, N.site))
> set.param(“s”, rep(1, 4))
>
> post.bugs <- call.bugs(n.iter = 10100, n.burnin = 100, n.thin = 10) > post.list <- to.list(post.bugs) > post.mcmc <- to.mcmc(post.bugs) > plot(post.bugs)

http://hosho.ees.hokudai.ac.jp/~kubo/ce/2008/car/model1.bug.txt
model
{
Tau.noninformative <- 1.0E-2 P.gamma <- 1.0E-2 for (i in 1:N.site) { Y[i] ~ dpois(mean[i]) log(mean[i]) <- beta + re[i] } re[1:N.site] ~ car.normal(Adj[], Weights[], Num[], tau) beta ~ dnorm(0, Tau.noninformative) tau ~ dgamma(P.gamma, P.gamma) } > source(“R2WBwrapper.R”)
# reading “R2WBwrapper.R” (written by kubo@ees.hokudai.ac.jp)…
> clear.data.param() # for initialization
>
> # set data
> load(“Y.RData”)
>
> v.na <- c(6, 9, 12, 13, 26:30) > YwithNA <- Y > YwithNA[v.na] <- NA >
> set.data(“N.site”, length(YwithNA))
> set.data(“Y”, YwithNA)
>
> Adj <- c(sapply(2:(N.site - 1), function(a) c(a - 1, a + 1))) > set.data(“Adj”, c(2, Adj, N.site – 1))
> set.data(“Weights”, rep(1, 2 * N.site – 2))
> set.data(“Num”, c(1, rep(2, N.site – 2), 1))
>
> # set parameter
> set.param(“beta”, 0)
> set.param(“re”, rnorm(N.site, 0, 0.1))
> set.param(“tau”, 1)
>
> post.bugs2 <- call.bugs( + file = "model1.bug.txt",　　#file.bug ->fileへ変更
+ n.iter = 300, n.burnin = 100, n.thin = 1
+ )

plotcar.R

> source(“plotcar.R”)
> plotcar(post.bugs2,v.na=v.na)

11.6 この章のまとめ
空間構造のあるデータを統計モデル化する場合、近傍とは似ているけれど遠方とは似ていないといった空間相関を考慮しなければならない。

空間相関のある場所差を生成するintrinsic Gaussian CARモデルはWinBUGSで簡単に扱える。

空間相関のある場所差は確率場を使って表現できる

空間相関を考慮した階層ベイズモデルでは、観測データの欠損部分をよそくするような用途にも使える。