JAGS と Stan と lm
ここ最近、GLM や GLMM の学び直しをしていました。みどりぼんも2周目になり、6～7割は理解できるようになったと思います。そこで、唐突ですが３ヶ月ほど前に Stan＆Rstan でやった、単回帰分析を JAGS＆rjags でやります。元ネタ（？）はこちら。 hikaru1122.hatenadiary.jp
JAGS を使うにあたって参考にした書籍はこちら（以下、DBDA2 と呼びます）。また、JAGSのユーザーマニュアルもざっと読みました。
Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan
- 作者: John Kruschke
- 出版社/メーカー: Academic Press
- 発売日: 2014/11/17
- メディア: ハードカバー
- この商品を含むブログを見る
DBDA2 の第8章でやっと JAGS を使った分析例が出てきます。DBDA2 による JAGS を使ったベイズ推定の手順は次のとおりです*1。勝手に3つの段階「準備する→MCMCする*2→結果を見る」にわけてみました。
＜準備する＞
１．分析するデータを準備する。
２．リスト形式でまとめる。…A
３．モデルを作る。
４．モデルを書き出す。…B
５．（初期値を設定する…C）
※５は必ずしも必要ではない。
＜MCMCする＞
１．JAGSにすべての情報（A,B,C）を渡す。
→ rjags::jags.model を使う。
２．バーンイン期間のために空回しする。
→ rjags::update を使う。
３．結果を記録する
→ rjags::coda.samples
＜結果を見る＞
・収束しているかどうか確認する。
→ 著者が用意してくれた diagMCMC を使う*3。
・事後分布を確認する。
→ mode, median, mean を確認。mode がいちばんいいけど、不安定なことももある。median がたいていいちばん安定している。
→ 著者が用意してくれた plotMCMC を使う*4。
では、やってみます。分析するデータは「アイスクリーム統計学」からお借りします。
http://kogolab.chillout.jp/elearn/icecream/chap4/sec1.html
＜準備する＞
icecream <- read.delim("icecream.txt") kyaku <- icecream$kyaku kion <- icecream$kion Ntotal <- length(kyaku) # JAGS に渡すデータをリスト形式で作る。 bunsekiData <- list( kyaku = kyaku, kion = kion, Ntotal = Ntotal ) # モデルを作る。 Model <-" model{ for(i in 1:Ntotal){ kyaku[i] ~ dnorm(mu[i], tau) mu[i] <- alpha + beta * kion[i] } alpha ~ dnorm(0, 1.0E-4) #事前分布はこれでいいの？ beta ~ dnorm(0, 1.0E-4) #事前分布はこれでいいの？ tau ~ dgamma(1.0E-3, 1.0E-3) #事前分布はこれでいいの？ } " #モデルを書き出す。 writeLines(Model, con = "icecreamModel.txt")
＜MCMCする＞
#rjags をロード。 library(rjags) # jagsオブジェクトを作る。初期値は設定せず、おまかせ。 haahaa <- jags.model(file = "icecreamModel.txt", data = bunsekiData, n.chains = 3, n.adapt = 500) # バーンイン期間のための空回し。 update(haahaa, n.iter = 5000) # MCMCサンプルを kekka に収める。知りたいパラメータは alpha（切片）と beta（偏回帰係数）。 kekka <- coda.samples(haahaa, variable.names = c("alpha", "beta"), n.iter = 100000)
＜結果を見る＞
source("DBDA2E-utilities.R") # これは DBDA2 の付録で付いてくる。 # DBDA2 著者作成の関数で収束診断。 diagMCMC(codaObject = kekka, parName = "alpha") diagMCMC(codaObject = kekka, parName = "beta") # DBDA2 著者作成の関数で事後分布を見る。mean, mode, median をそれぞれ見てみる。 par(mfrow=c(2,3)) plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "mean") plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "mode") plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "median") plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "mean") plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "mode") plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "median")
収束はしているようです（下図のそれぞれ左上と左下）。
切片と係数は次のとおりです。それぞれ3つ出してみます。
median を採用するとしたら、回帰式は
ということですよね。
同じデータを lm でやった結果は次のとおりです。
> kekka_lm <- lm(kyaku ~ kion, data=icecream) > summary(kekka_lm) Call: lm(formula = kyaku ~ kion, data = icecream) Residuals: Min 1Q Median 3Q Max -47.969 -17.709 -1.218 17.413 51.031 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -229.98 73.79 -3.117 0.00596 ** kion 17.25 2.30 7.499 6.08e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 29.54 on 18 degrees of freedom Multiple R-squared: 0.7575, Adjusted R-squared: 0.744 F-statistic: 56.23 on 1 and 18 DF, p-value: 6.082e-07
Stan での結果も再掲します。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha -230.75 0.57 80.23 -389.66 -282.82 -230.41 -178.89 -72.06 19789 1 beta 17.27 0.02 2.50 12.32 15.66 17.26 18.90 22.23 19771 1
Stan＆Rstan でやると lm とほぼ同じ結果が得られたんですが、今回の JAGS だとなんか違います。事前分布の設定が悪かったのかな…。それとも初期値の設定が必要だったのか。それとも分析するデータが少なすぎるとか。やっとモデルを作る、見るのに慣れてきたと思ったら、また次の難問が出てきました。厳しい…。
しかし、Stan よりは JAGS のほうがとっつきやすい気がします。しばらく慣れ親しんでみようと思います。
