統計学エンドユーザーのための JAGS 入門

岩波データサイエンス vol.1 のどこかのページに「初心者はJAGS（BUGS）」と書いてあったので*1、こつこつ学び始めています。世の中は Stan 大人気ですが、気にしない…。まず JAGS に慣れてから Stan をやります…。本エントリーは、前回のエントリーと多くがかぶっています。

　
今回は、次のページを参考にさせていただいています。

生態学のデータ解析 - R で JAGS (rjags 編)
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html

生態学のデータ解析 - JAGS 雑
http://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html

　
分析に使用したのは、JAGS 4.0.0 および rjags 4.4 です。

JAGS と rjags のインストール

JAGS と rjags はそれぞれ別にインストールします。JAGS 3.0 系をインストールしていて、rjags 4.4 にしていた場合は動きません。JAGS も４にする必要があります。

　
JAGS はこちらから .exe をダウンロードできます。
http://sourceforge.net/projects/mcmc-jags/files/

rjags はいつもの install.packages でOKです。

分析の手順

１．データの準備

今回は先ほどの「RでJAGS（rjags 編）」のデータを使用しました。上記ページにある csv ファイルから、ページのとおりにデータをセットします。JAGSに渡すデータはリスト形式である必要があります。

２．モデルを作る

「RでJAGS（rjags 編）」にあるモデルを自分なりにいじってみました。いまのJAGSマニュアルによると、行末にセミコロンは不要のようです。モデルを作る前に図解しました。これも岩波データサイエンスにそうしろ、と書いてあったように思います。分布の形はテキトーです*2。

f:id:hikaru1122:20151120002559j:plain

　
さて、モデルです。ロジスティック回帰でリンク関数はロジットリンクです。詳しくはみどりぼん第６章を参照してください。

Model <-"
model{
    for (i in 1:N) {
        Positive[i] ~ dbin(p[i], Total[i])
        logit(p[i]) <- z[i]
        z[i] <- beta[1] + beta[2] * Logdose[i]
    }
    beta[1] ~ dnorm(0, tau)
    beta[2] ~ dnorm(0, tau)
    tau ~ dgamma(1.0E-6, 1.0E-6)
}
"

追記2015年11月20日
そういえば、JAGS の dnorm の分散は逆数でした。
tau <- 1/(sd*sd)
sd ~ dgamma(1.0E-3, 1.0E-3)
がより正しい書き方なのかな…？
また、twitter にてアドバイスをいただきました。
一様分布や半コーシー分布のほうがよい、とのことです。
ありがとうございます。

　
いちばんの変更点は tau の事前分布です。一度、同じように dgamma(1.0E-3, 1.0E-3) でやってみたのですが、うまく行きませんでした。代わりに、dgamma(1.0E-6, 1.0E-6) を使いました。

　
その理由は2つあります。
（１）JAGS 4.0 マニュアルに載っていた例で dgamma(1.0E-6, 1.0E-6) が使われていたこと
（２）前回のエントリーで dgamma(1.0E-6, 1.0E-6) を使ってやり直したら、Stan とほぼ同じ結果が得られたこと
です。1.0E-3 と 1.0E-6 で結果が変わるのってシビア。事前分布って難しい。

　
モデルを書いたら、書き出します。

writeLines(Model, con ="testModel.txt") # 拡張子は特に不要ですが、いちおう。

３．MCMCする

モデルを jags オブジェクトにしていきます。

# rjags のロード。
library(rjags)
# jags モデルを作る。
testModel <- jags.model(file = "testModel.txt", data = list.data, n.chain = 3)

　
jags.model の file = は別ファイルである必要があります。なので、先ほど書きだした次第です。次にバーンイン、そして本格的なサンプリングを行います。n.chain は２以上でないと、後で Gelman & Rubin 収束診断指標のグラフが描けないようなので、３にしています。初期値は設定していません。まだ理解できていないので。

　
では、本格的なMCMCサンプリングを行います。

# バーンイン期間のための試運転。
update(testModel, n.iter = 1000)
# 本番。MCMC サンプルを kekka に収める。
kekka <- coda.samples(
            testModel,
            variable.names = c("beta[1]", "beta[2]", "tau"),
            n.iter = 20000)

　
coda.samples の variable.names では結果で見たいパラメータを文字列として指定します。これで、下ごしらえは終わりました。rjags パッケージで使うのは、jags.model, update, coda.samples の３つだけです。writeLines は普通に使えます。

結果を見る。

収束診断をします。codamenu() でメニューを選びながら、確認します。引数はいらないです。メニューで「2→kekka（MCMCサンプリングを格納したオブジェクト）→2→2」とするとGelman & Rubin の指標が見られます。1.1 より下なので、収束しているようです。他にもcodamenu()でいろいろ指標を選べます。

Potential scale reduction factors:
        Point est. Upper C.I.
beta[1]       1.01       1.02
beta[2]       1.01       1.02
tau           1.00       1.00

　
図でも確認します。plot(kekka)です。

f:id:hikaru1122:20151120011135p:plain

　
推定値は summary(kekka) で確認します。

> summary(kekka) 
            Mean      SD  Naive SE Time-series SE
beta[1] -5.99179 1.10916 0.0045281       0.053419
beta[2]  1.18754 0.19647 0.0008021       0.009398
tau      0.06042 0.07121 0.0002907       0.001165