生態学のデータ解析 - R2WinBUGS

• ここは WinBUGS を R から使う package R2WinBUGS (CRAN 内) についてのあれこれ書いてます．
  • JAGS を R から使う方法については R で JAGS (rjags 編) 参照
  • ラッパー関数については R2WBwrapper 参照

    

• 参照: ベイズ統計 & MCMC, WinBUGS のインストール, WinBUGS 雑, bugs() 関数の返り値, coda 雑, 例/car.normal()

[もくじ]

• R2WinBUGS とは
• R2WinBUGS をとりあえず動かす (Windows user むけ)
• R2WinBUGS をとりあえず動かす (Linux user むけ)
• R2WinBUGS をとりあえず動かす (MacOS X user むけ)
• WinBUGS の解除 key
  • bugs() 関数の返り値 (bugs オブジェクト)
• 計算結果の処理: R の library(coda) との連携
  • 計算結果の処理: chain ごとに結果を取りだす

R2WinBUGS とは 

• ベイズモデルの MCMC 計算をする Gibbs sampler ソフトウェアのひとつである WinBUGS を R から使う package R2WinBUGS です
  • library(help = R2WinBUGS)
  • WinBUGS のインストール
  • R ユーザーにとっては WinBUGS を直接つかうより，R から R2WinBUGS を介在して WinBUGS を使ったほうが便利だろう，という理由で作られた package です
  • というのも，WinBUGS は Gibbs sampling 専用ソフトウェアで，他の統計学的データ解析ツールはほとんど入っていないので
  • R で WinBUGS にわたすデータを整えたり，あるいは WinBUGS の計算結果 (bugs() から bugs オブジェクトとして返されます) を作図・解析するの，という使いかたをします
  • しかしながら library(R2WinBUGS) もかなり 使いにくい ので，私 (久保) はR2WBwrapper.R というラッパー関数を定義して使っています
    • 参照: R2WBwrapper, 例/car.normal()
• RjpWiki 内の解説
  • Rでベイズ統計学
• lindley さんの解説 blog: R と WinBUGS http://d.hatena.ne.jp/lindley/
• WinBUGS については 本/環境科学の階層モデル なども参考になります．
• library(R2WinBUGS) のインストール方法: R を起動して
  • update.packages()
  • install.packages("R2WinBUGS")
  • WinBUGS のインストールは WinBUGS のインストール を参照

R2WinBUGS をとりあえず動かす (Windows user むけ) 

• Windows Vista, 7 以降の注意
  • R を起動するときに，右クリックして 「管理者」 として起動する
• help(bugs) に掲載されてる example が実行するには……

library(R2WinBUGS) # R2WinBUGS のよみこみ
example(bugs) # モデルファイルやパラメーター・データが準備される
schools.sim <- bugs(
	data, inits, parameters, model.file,
	n.chains = 3, n.iter = 5000,
	bugs.directory = "c:/Program Files/WinBUGS14/",
	working.directory = NULL,
	clearWD = TRUE,
	debug = TRUE # これは FALSE でも動作する
)
print(schools.sim) # bugs オブジェクトのまとめ
plot(schools.sim) # bugs オブジェクトの作図

R2WinBUGS をとりあえず動かす (Linux user むけ) 

• まず Wine をインストール (参照: ぎょーむ日誌 2006-10-26)
• winecfg で $HOME/.wine/ を作ってから wine WinuBUGS14.exe で WinBUGS をインストール
• WinBUGS 1.4→1.4.3 の patch をあてる
• help(bugs) に掲載されてる example が実行するには (wine が /usr/local/bin/ にインストールされてる場合) ……

library(R2WinBUGS) # R2WinBUGS のよみこみ
example(bugs) # モデルファイルやパラメーター・データが準備される
WINEPATH <- "/usr/local/bin/winepath" # ここで定義する必要あり
schools.sim <- bugs(
	data, inits, parameters, model.file,
	n.chains = 3,
	n.iter = 5000,
	bugs.directory = paste( # そこでこのように変更
		Sys.getenv("HOME"),
		".wine/drive_c/Program\ Files/WinBUGS14/",
		sep = "/"
	),
	working.directory = NULL,
	clearWD = TRUE,
	program = "WinBugs",
	useWINE = TRUE,
	newWINE = TRUE,
	WINE = "/usr/local/bin/wine",
	WINEPATH = WINEPATH,
	debug = TRUE # これは FALSE でも動作する
)
print(schools.sim) # bugs オブジェクトのまとめ
plot(schools.sim) # bugs オブジェクトの作図

R2WinBUGS をとりあえず動かす (MacOS X user むけ) 

• 伊東さんの調査結果 http://blog.so-net.ne.jp/ito-hi/2006-11-23 などを参考にされてください

WinBUGS の解除 key 

• educational version cannot do this model というエラーが出たら
  • WinBUGS のユーザー登録 (無料)
  • 解除キーがメイルで送られてくるので WinBUGS で open して decode (で解除される)

bugs() 関数の返り値 (bugs オブジェクト) 

• DIC などが含まれている
• 一覧

	 [1] "n.chains"        "n.iter"          "n.burnin"        "n.thin"         
	 [5] "n.keep"          "n.sims"          "sims.array"      "sims.list"      
	 [9] "sims.matrix"     "summary"         "mean"            "sd"             
	[13] "median"          "root.short"      "long.short"      "dimension.short"
	[17] "indexes.short"   "last.values"     "pD"              "DIC"            
	[21] "DICbyR"          "model.file"      "is.DIC"          "program"        

• 特に重要なのは:
  • sims.matrix: MCMC 計算で得られたサンプルが格納されている (下記の coda との連携を参照)
  • sims.array: MCMC 計算で得られたサンプルが chain ごとに格納されている
  • last.values: いちばん最後のサンプルが格納されている
• 詳しくは bugs() 関数の返り値 参照
• print(bugs オブジェクト) するときに digits.summary=3 といったオプションを追加すると表示されるケタ数が増える

計算結果の処理: R の library(coda) との連携 

• R2WinBUGS の bugs オブジェクトを coda の mcmc オブジェクトに変換する:

	schools.mcmc <- as.mcmc(schools.sim$sims.matrix)

• MCMC 計算の結果処理に便利な library(coda) については coda 雑 を参照

計算結果の処理: chain ごとに結果を取りだす 

http://rgonzui.nakama.ne.jp/CRAN/markup/R2WinBUGS/R/bugs.plot.inferences.R

などをみればわかるように，bugs オブジェクトの中で， chain ごとの結果は sims.array に格納されている．

example(bugs) で生成される例題で考えてみよう． ss が bugs() で得られた bugs オブジェクトだとしよう． 格納されている step 数が 358，n.chain = 3, パラメーター数 11 なので ss$sims.array は

	> dim(ss$sims.array)
	[1] 358   3  11

となっており [step,chain,parameter] という三次元配列であることがわかる． したがって一番目の chain の第二パラメーターをとりだしたければ

	ss$sims.array[, 1, 2]

とすればよく，またこれを library(coda) の gelman.diag() に渡すには

library(coda)
ss.list <-  mcmc.list(
	lapply(
		1:ss$n.chain,
		function(chain) as.mcmc(ss$sims.array[, chain,])
	)
)

なる mcmc オブジェクトの list (つまり mcmc.list) を作って gelman.diag(ss.list) とすれば良い．

• この mcmc.list オブジェクトを plot() すると chain ごとに経時変化が作図される

(coda 雑 参照)