生態学のデータ解析 - 例/car.normal()

# runbugs1.R は下のようにデータの準備，パラメーターの初期化など
# そして WinBUGS の呼び出しとその結果のうけとりを実行する
source("R2WBwrapper.R") # R2WBwrapper.R の読みこみ
clear.data.param() # for initialization

# set data (WinBUGS にわたすデータの準備)
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)) # 「近傍」ID の定義
set.data("Weights", rep(1, 2 * N.site - 2)) # 近傍の重み (全部 1)
set.data("Num", c(1, rep(2, N.site - 2), 1)) # 近傍の個数 (1, 2, 2, ... , 2, 1)

# set parameter (WinBUGS にわたすパラメーター初期値などの準備)
set.param("beta", 0)
set.param("re", rnorm(N.site, 0, 0.1))
set.param("tau", 1)

post.bugs1 <- call.bugs( # ここまで準備したデータとパラメーターをわたして
    file.bug = "model1.bug.txt", # model1.bug.txt で定義されているベイズモデルを
    n.iter = 200, n.burnin = 100, n.thin = 1 # このような MCMC step で
) # 推定計算して結果を post.bugs1 に格納しろ，という命令

• R2WBwrapper.R: runbugs1.R で使う RtoWinBUGS のラッパー関数定義ファイル
  • library(R2WinBUGS) は何とも使いづらいのでこういう「R2WinBUGS をラップする関数群」を定義した

R2WinBUGS で推定計算 

• 上のように準備ができれば，あとは runbugs1.R を R で実行させるだけです
  • たとえば下のように R のコマンドラインで source("runbugs1.R") と指示すると， runbugs1.R の読み込み，その内容を実行します
• WinBUGS が計算を終えると何やら奇妙なグラフをだして停止するので，WinBUGS を「終了」します
• すると WinBUGS の推定計算結果が R 内の post.bugs1 に格納されます

WinBUGS で得られた推定結果を R で調べる 

• plot(post.bugs1) とすると，下のような図が得られます
  • 左側は収束のよさを R-hat 値であらわし，右は事後分布の中央値と 80% 区間が示されています

• print(post.bugs1,digits.summary=3) とすると mcmc1.txt にあるような事後分布表が出力されます

R を使って事後分布の作図 

• plotcar.R を source("plotcar.R") してから plotcar(post.bugs1) とすると，下のような図が出力されます

補足 1: 観測データに欠測がある場合 

• 下の図のように欠測 (黒丸) がある場合の推定計算方法について説明します
  • 黒丸は欠測つまり観測されなかった観測値です
• じつはこれは (BUGS 言語で書かれた) 統計モデルを変更する必要はありません
  • 観測データ (ここでは Y) に欠側値 NA をまぜてやるだけで OK です
    • WinBUGS user manual 参照してください

必要なファイル 

• Y.RData: 上の例と同じ
• model1.bug.txt: 上の例と同じ (変更の必要なし!)

• runbugs2.R: 上の例の runbugs1.R をちょっとだけ変更したもの

# (runbugs1.R と runbugs2.R の相違)
# set data
load("Y.RData")

v.na <- c(6, 9, 12, 13, 26:30) # 欠測の生じた場所
YwithNA <- Y
YwithNA[v.na] <- NA # つまりデータ vector に NA があると
                    # WinBUGS はその地点を欠測とみなす

set.data("N.site", length(YwithNA))
set.data("Y", YwithNA)

(途中おなじ，中略)

post.bugs2 <- call.bugs( # post.bugs2 に結果を受け取らせている
	file.bug = "model1.bug.txt", # BUGS コードは上の例と同じ，変更の必要なし
	n.iter = 200, n.burnin = 100, n.thin = 1
)

> # つまりデータがこう変わるということです
> # (欠測なしの Y)
> Y
 [1]  0  3  2  5  6 16  8 14 11 10 17 19 14 19 19 18 15 13 13  9 11 15 18 12 11
[26] 17 14 16 15  9  6 15 10 11 14  7 14 14 13 17  8  7 10  4  5  5  7  4  3  1

> # (欠測ありの Y)
> Y
 [1]  0  3  2  5  6 NA  8 14 NA 10 17 NA NA 19 19 18 15 13 13  9 11 15 18 12 11
[26] NA NA NA NA NA  6 15 10 11 14  7 14 14 13 17  8  7 10  4  5  5  7  4  3  1

• R2WBwrapper.R: 上の例と同じ

R2WinBUGS で推定計算，結果の作図 

• 上のように準備ができれば，あとは runbugs2.R を R で実行させるだけです
• 今回は推定計算の結果が R 内の post.bugs2 に格納されます
• plotcar.R を source("plotcar.R")
  • 注意: この plotcar.R は 2008-03-25 に改訂 (欠測値表示をできるようにした) したものです (古い version のものをもっている人はいれかえてください)
  • 次に plotcar(post.bugs2,v.na=v.na) とすると，下のような図が出力されます
    • v.na は runbugs2.R の中で定義されている「欠測のあった位置」です

補足 2: 各地点独立とした場合 

• ついでながら，空間相関なし random effects，つまりごくふつーの階層ベイズモデルで「場所差 random effects」(各地点完全独立を仮定している) を表現する場合はどう計算するのかを説明します

必要なファイル 

• Y.RData: 上の例と同じ
• model3.bug.txt: 空間相関なし版モデル
  • 蛇足: 下の BUGS code はわかりやすさを重視した書きかたで，じつはもっと計算速度が速く収束もよい書きかたがあります (参照: ぎょーむ日誌 2008-09-08) -- ただしこれは dnorm() などでしか使えません (car.normal() では使えない)

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[i] ~ dnorm(0.0, tau) # 各地点独立な random effects
	}
	beta ~ dnorm(0, Tau.noninformative)
	tau ~ dgamma(P.gamma, P.gamma)
}

• runbugs3.R: 上の例の runbugs1.R をちょっとだけ変更したもの

source("R2WBwrapper.R")
clear.data.param() # for initialization

# set data
load("Y.RData")
set.data("N.site", length(Y))
set.data("Y", Y)
# Adj, Weights, Num の設定は不用

# set parameter
set.param("beta", 0)
set.param("re", rnorm(N.site, 0, 0.1))
set.param("tau", 1)

post.bugs3 <- call.bugs( # post.bugs3 に結果を受け取らせている
	file.bug = "model3.bug.txt",
	n.iter = 200, n.burnin = 100, n.thin = 1
)

• R2WBwrapper.R: 上の例と同じ

R2WinBUGS で推定計算，結果の作図 

• 上のように準備ができれば，あとは runbugs3.R を R で実行させるだけです
• 今回は推定計算の結果が R 内の post.bugs3 に格納されます
• plotcar.R を source("plotcar.R")

> # (R 内で)
> plotcar(post.bugs3, col1 = "#0000ff", col2 = "#ccccff")