このページは、R上からWinBUGSを操作することを可能にするパッケージ、R2WinBUGSの使い方を解説するページです。
モデル式の部分を加筆中
累積訪問者数（from 2008/12/17）： -
今日来てくれた人： -
計算機の進歩により、Bayes統計（というかMCMC計算）が可能になり、その利用に関する関心が急速に高まっています。しかし、実行方法の難しさが多くのユーザーへの拡大を妨げているといえるでしょう。そこでMCMC（Gibbsサンプリングのみですが.......）計算を行うためのWinBUGSが開発されました。WinBUGSは単体で結果の表示、作図などを行ってくれます。しかし、それら以外の結果が取り出したいケースもありますし、何よりRで使えたら論文などへの二次利用も容易となるでしょう。
そこで開発されたのがR2WinBUGSです。名前は本当にそのまんまと言う感じです（2はtoをもじったものでしょう）。R2WinBUGSを使うと、データはRで準備すればよく、WinBUGSでのMCMC計算の過程などをRに格納することが可能となり、利便性が格段に向上します。
といっても、R2WinBUGSはR2WinBUGSでいろいろお作法というか指定しなければならないことがあります。ここでは初心者が陥りがちな点を記述するつもりです。
なお、armというパッケージは、R2WinBUGSや作図で使うcodaなど、複数のパッケージの集まりです。個別にパッケージを導入するのが面倒くさい方は（私もその一人＾＾；）、armを導入するといいと思います。
各部分の詳細について以下で述べますが、せっかくですから実データを使って解析をしてみましょう。使うデータはこちら。
R.csv
ここがある意味最も難しいところです。自分が考えているモデルを表現する必要がありますが、なれないうちは中々書くことができません。こればかりは慣れですので、頑張って習得してください。
今回使うデータの解析の目的は、
植物に窒素を与える区と与えない区で、葉の窒素濃度に差が出るか？
という問題です。窒素を与える区と与えない区でそれぞれ同種の植物が数個体あり、各個体から5枚ほど葉がサンプリングされています。
実際にファイルを開いてみると、データは以下のようになっていることがわかります。
BUGSコードを記述するファイルは、以下のように作る必要があります。
model {
......（モデルの中身）
}
このように記述したファイルを、.bugという拡張子をつけて保存します。
で、問題となるのはモデルの中身の書き方ですね。それを以下で解説します。
要するに代入する部分です。
<-
で指定すると決定論的指定になります。決定論的指定にしては誤差などを挟む余地はないので、主に自分が考えているモデルを記述する部分になります。
今回の仮説は「植物に窒素を与える区と与えない区で、葉の窒素濃度に差が出るか？」でした。なので単純に考えられるモデルを表現すると、以下のようになりそうです。
Ncon[i] <- a + coef * Ntreat[i] #aは切片（窒素が与えられないときの葉の窒素濃度）、coefは係数、iは各行を示す
aとcoefだけ決まれば葉の窒素濃度がわかりそうです。でも、こんなに現象は単純でしょうか？例えば以下のようにして図を描いてみると、
> plot(Ncon ~ Ntreat, col=No, d) #Noは各個体
Ntreatだけではっきりと区別できているわけではないですね。また、色の違いは個体の違いですが、個体によってずいぶんと窒素濃度が違うように見えます。
そうすると、まずモデルの予測値とデータの間にはある程度の誤差が含まれ、かつ個体によっても葉の窒素濃度に違いがありそうです。それを記述しようとすると、
Ncon[i] <- mu[i] + error[i] mu[i] <- a + coef * Ntreat[i] + catID[No[i]] #mu[i]はモデルの予測値、catIDは個体によって異なるカテゴリカル変数
となりそうです。
しかし、<-は代入ですから、未知の切片（a）や窒素処理の有無の係数（coef）、誤差（error）、個体間差（catID）は表現することができません。そこで、以下の確率論的指定が必要になります。
確率分布を使う箇所は、確率論的指定をする必要があります。
~
で指定すると確率論的指定になります。具体的には、
さて、これを踏まえて先ほどのモデルを考えると、aやcoef、error、catIDは確率論的指定を使えば表現できそうです。ただし、確率論的指定という名称からもわかるように、「どの確率分布」を使うかを決める必要があります。
特に、事前情報が無くどのような確率分布を使うかわからないときは、できるだけ無情報な確率分布を当てはめることが多いです。今回は、「平均が0、非常に分散が大きい正規分布」を用いることにします。
Ncon[i] ~ dnorm(mu[i], 0.001) #確率論的指定 mu[i] <- a + coef * Ntreat[i] + catID[ID[i]] #決定論的指定
a ~ dnorm(0.0. 0.001) #確率論的指定 coef ~ dnorm(0.0, 0.001) #確率論的指定 catID[1] ~ dnorm(0.0, 0.001) #確率論的指定 ...... catID[21] ~ dnorm(0.0, 0.001) #確率論的指定
さてさて、いきなりレベルが上がったように感じると思いますが.......
これらがわかったところで、もう一度上記のモデルを検討しましょう。気になるのは個体差の係数catIDですね。このままだと個体ごとに係数を独立に推定することになり、パラメータを21個も消費することになります。
しかし実際には個体差は興味のある対象ではないので、21個もパラメータを消費して個別に推定したいわけではありません。そこで、「個体差の係数は全個体でみると平均値0、分散σの正規分布に従う」と考えます。こうすると、正規分布に従うという制約が加わりますので、パラメータ数の消費を抑えることができます。
for (id in 1:21) {
catID[id] ~ dnorm(0.0, tau)
}
またまたわけのわからない形に見えますが.......
さて、この21個の係数達の関係（似ているか似ていないか）は、係数の分散が小さいか大きいかという問題と言えます（分散が大きければ
この分散については特に事前の情報がないので、やはりある確率分布に従う値として推定させます。ここでどのような確率分布を使うかですが、使われる分布は正規分布との積になるので、共役事前分布を使うと計算が楽になります。
共役事前分布とは、別の確率分布との組み合わせによっては、積の前後で同じ確率分布になる物です。よく知られているものとしては、
for (id 1:21) {
catID[id] ~ dnorm(0.0, tau)
}
tau ~ dgamma(0.001, 0.001)
完成品がこちら。bayes.bug
R上ではデータをデータフレーム形式で用意していることがほとんどかと思いますが、bugs()に渡すデータは以下の点に気をつけて準備してください。
後述しますが、WinBUGSが途中でこける場合の原因の多くが、初期値設定にあります。stochastic node（モデル式中で ~ で表されている部分）には初期値をきちんと与えましょう。
また、繰り返し回数やパラメータ数が多い時に全てのパラメータを結果としてみようとすると、bugs()計算が終了してRに戻るときに、「メモリーがいっぱいです」などと言われて結果が見れないことがあります。初期値はstochastic node全てに必ず与えないといけませんが、見るパラメータは絞りましょう。
基本的には以下のbugs()関数を使ってMCMC計算を行います。bugs()はR2WinBUGSに含まれる関数です。
> res <- bugs(d, inits, para, "***.bug",
n.chains=3, n.iter=3000, n.thin=10, n.burnin=1000,
bugs.directory="c:/自分がインストールしたディレクトリ",
working.directory=getwd(),
debug=TRUE)
それぞれの項目の意味ですが、
各種回数を設定するところが最も悩ましいところかと思います。どのぐらいで収束するかはデータとモデルのよしあしから決定されます。n.chainsは複数本とするのが普通で、データが大きいなど特別な理由がない限り1にはしないで下さい。普通は3ぐらいだと思います。
繰り返し数を増やせば収束しやすくなりますが、その分計算時間が増大します。今回のように少ないデータであれば問題となりませんが、1000を超えるようなデータではよく考えないと、計算時間が数時間に及びます。
n.thinはサンプリングの自己相関を除くために設定しますが、最終的に得るデータの数を減らす（きちんとシミュレーションした上で）という意味あいもあります。
慣れない内は上手くまわせないことが多いですから、debug=TRUEとしておきましょう。WinBUGSにサンプリングの図が出てくれば計算はきちんと終了していますので、WinBUGSを閉じるとRで作業できるようになります。一方、失敗した場合も、debug=TRUEとしておけばWinBUGSは立ち上がったままになりますから、原因を特定するのに役立ちます（見てもわからない場合も多々ありますが＾＾；）。
ちゃんと回ってよかったよかった、というわけにはいきません。収束していなければ意味がありません。収束を確認する最もよい方法は、MCMC計算の過程や事後分布を作図によって図示することです。
たとえばresというオブジェクトにbugs()の結果を付与している場合、
> res
と打つと要約された推定値などが出力します。表示される桁数を変えたい場合は、
> print(res, digits.summary=*)
と打ちます。*の数字が出力桁数になります。
これを図で見る場合、
> plot(res)
とします。これを行うと、各パラメータの80%信頼区間（何故80%かはわかりませんが）、収束判定に使うR-hat値（古谷 2008を見てください）などが表示されます。
> names(res)
とすると、bugs()の結果にどのような情報が含まれているのかわかります。このうち、
こちらはchainの平均の値が使われていますから、厳密な収束診断に用いるべきではありませんが、簡単に描画することができます。
> plot(as.mcmc(res$sims.matrix))
とすると、パラメータごとのサンプリング過程と事後分布図が表示されます。
as.mcmc()というのは、データをmcmcオブジェクトという形式に変換するための関数です。mcmcオブジェクトをplot()すると、自動的にサンプリング過程や事後分布図が表示される、という仕組みです。
こちらはchainごとのサンプリング値が取り出せるので、収束診断として最も適当な図です。ただし、図にするには少々工夫が要ります。
chainごとのデータは配列（array。行列matrixは二次元だが、それを多次元に拡張したもの）で保存されています。これをplot()で読める形に変換します。
res$sims.arrayを表示させると、3次元のデータになっていることがわかります。ではそれぞれの次元は何を表しているかというと、
res$sims.array[サンプリング数, chain数, パラメータ数]
となっています。
先ほどのmcmcオブジェクトに変換するための関数as.mcmc()は、一度に1chainのデータしか変換することができません。そこで、指定した列に関数を適用できるlapply()を使い、各chainのデータをmcmcオブジェクトに変換します。
lapply()は、
> lapply(使いたいリスト形式のデータ, 関数)
という書き方が基本です。ただし、今回のように、as.mcmc()を全部に適用するのではなく、各chainごとに、というように、自分で動作を指定する場合は、
> lapply(関数を適用する列番号, function(x) .....(動作の中身))
という書き方をします。
lapply()は結果をリスト形式で返しますが、そのままでは描画することができません。そこで、mcmc.list()という関数を使ってmcmcリストというこれまた変な形式のリストに変換して、一つのリストとします。これをplot()に放り込めば望むような図が出力されます。
以上を踏まえ、
> plot(mcmc.list(lapply(1:res$n.chains, #1からres$n.chainsまで（要は指定したchain数）
function(x) as.mcmc(res$sims.array[, x, ])
)
)
)
とすると、chainごとに自動で色分けされたサンプリング過程と事後分布図が表示されます。
