緑本メモ
MCMCで混合効果モデルのパラメータ推定をやりたい。
R(+外部アプリ)でやろうとするととりあえず以下の方法があるようで。
(Task-viewを読む限り他にもいろいろあるけど、日本語情報の量等考慮して以下を選んだ。)
- MCMC+パラメータ推定を一括でRでやってしまう
- lmmパッケージを使う
- MCMCglmmパッケージを使う
- lme4パッケージを使う
- MCMCを外部アプリに任せる
- WinBUGS(JAGS/OpenBUGS)をR2WinBUGS(rjags or R2jags/BRugs)から使う
一括でRでやる方が簡単そうだけど、今読んでいる本(http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html 以降、緑本と呼ぶ)が最後の方法を用いているので勉強がてらこれでいく。
ただしWinBUGSはWindows7にインストールできないのでJAGSを使う。
JAGSを用いる場合、本に掲載されているWinBUGSのコードとの相違点からいくつか詰まりそうなので本を読みながら適宜メモしていく。
JAGSをインストールする
JAGSをインストールしていると、問題が出て動作を停止しましたとかなるがインストールはされている。
PATHを通してやればOK。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html#toc5
RからJAGSを使うためのパッケージをインストールする
rjagsとR2jagsの2つのパッケージがある。
後者はR2WinBUGSを真似たものとのこと。
http://yusung.blogspot.jp/2008/05/r2jags-package-for-running-jags-from-r.html
R2WinBUGSは色々と使いにくい(以下参照)という話があるので、rjagsでとりあえずいく。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoWinBUGS.html#toc1
RからJAGSを使うチュートリアルを眺める
とりあえず日本語で読む。
(日本語)http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html#toc4
(日本語)http://rmecab.jp/wiki/index.php?cmd=read&page=R_Bayes_jags&word=rjags
http://streaming.stat.iastate.edu/~stat444x/R%20and%20BUGS%20Info/using-jags-inR.pdf
http://www.rochester.edu/college/psc/thestarlab/help/JAGS.pdf
WinBUGSとの相違点に注意しながら本の実例を試す
JAGSの場合、行末にセミコロンが必要。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html#toc4
とりいそぎ本のコードについては行末にセミコロンを加えるだけで動いた。
むしろR2WinBUGsからrjagsへの書き換えの方が面倒。
(コードは一番下)
収束判断するためのR-hatについて
収束判断の際のRhatの具体的な値についての記述はこちらに。
http://d.hatena.ne.jp/takehiko-i-hayashi/20100512/1273643029
緑本にも書いてあるけど。
なお、rjagsからだとR-hatの値が見つからない(暫定)のでどうしても欲しければR2jagsで取得する。
とりあえずのまとめ
自由にモデリングしたい!って場合じゃなければlmmとかMCMCglmmを使った方が楽かもしれないので、そっちも調べる。
その他眺めた関連資料
北大の久保先生(今読んでいる本の著者)のまとめが有用。
http://hosho.ees.hokudai.ac.jp/~kubo/ce/LinksGlmm.html
http://hosho.ees.hokudai.ac.jp/~kubo/ce/BayesianMcmc.html
http://hosho.ees.hokudai.ac.jp/~kubo/stat/2010/ism/kuboISM2011feb.pdf
CRAN task-viewのベイズ推定の項
http://cran.r-project.org/web/views/Bayesian.html
GLM他のメモ集
http://www.lowtem.hokudai.ac.jp/plantecol/akihiro/obenkyou/R-study_%20Linear_Models_101203.pdf
混合効果モデルを説明する際のモデルデータ
http://ito-hi.blog.so-net.ne.jp/2011-07-23-1
ベイズ統計学のMCMCとの出会い(製薬)
http://www.sas.com/offices/asiapacific/japan/usergroups/wg/archive/041015mati.pdf
Rでマルコフ連鎖モンテカルロ法
http://d.hatena.ne.jp/isseing333/20100611/1276235951
混合効果モデル(変量効果モデル、mixed effect model)について
http://d.hatena.ne.jp/isseing333/20110413/1302695785
GLMMの実装比較(6年前の時点)
http://ito-hi.blog.so-net.ne.jp/2006-10-30-1
マルコフ連鎖モンテカルロ法入門
http://www.slideshare.net/teramonagi/ss-5190440
http://www.slideshare.net/teramonagi/ss-5344006
近似ベイズ計算によるカジュアルなベイズ推定
http://www.slideshare.net/kos59125/ss-9290131
FMEパッケージを用いたMCMCでのモデルのパラメータ推定
http://shimana7.seesaa.net/article/267132804.html
コード
library(rjags) library(R2WinBUGS) # BUGSファイルの書き出しにのみ利用 # 緑本サンプルデータのダウンロード tmpd <- tempdir() setwd(tmpd) download.file(destfile = "ch9.RData", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData") load("ch9.RData") # BUGSファイルの書き出し:write.model関数はR2WinBUGSパッケージから samplemodel <- function(){ for(i in 1:N){ Y[i] ~ dpois(lambda[i]); log(lambda[i]) <- beta1 + beta2 * (X[i] - Mean.X); } beta1 ~ dnorm(0, 1.0E-4); beta2 ~ dnorm(0, 1.0E-4); } filename <- file.path("sample.bug") R2WinBUGS::write.model(samplemodel, filename) # データ list の準備 list.data <- list( Y = d$y, X = d$x, Mean.X = mean(d$x), N = nrow(d) ) # 初期値 list の準備 inits <- list( beta = c(0, 0) ) # R 内でのモデルの定義 m <- jags.model( file = filename, data = list.data, inits = list(inits, inits, inits), n.chain = 3 ) # burn-in のためのカラまわし MCMC 計算 update(m, 100) # MCMC 計算で事後分布からサンプリング,その結果をうけとる post.list <- coda.samples( m, c("beta1", "beta2"), thin = 3, n.iter = 1500 ) summary(post.list) # 図9.5 plot(post.list) # 図9.6(A) beta1_smp <- c(post.list[,"beta1"][[1]], post.list[,"beta1"][[2]], post.list[,"beta1"][[3]]) beta2_smp <- c(post.list[,"beta2"][[1]], post.list[,"beta2"][[2]], post.list[,"beta2"][[3]]) beta1_median <- summary(post.list)$quantile[,"50%"][[1]] beta2_median <- summary(post.list)$quantile[,"50%"][[2]] plot(d, xlab = "植物の体サイズ", ylab = "種子数") par(new= TRUE) for(i in 1:length(beta1_smp)){ curve(exp(beta1_smp[i] + (x - mean(d$x)) * beta2_smp[i]), from = 3, to = 7, ylim = c(3, 12), col=rgb(0.5, 0.5, 0.5, alpha=0.1), xlab = "", ylab = "") par(new= TRUE) } curve(exp(beta1_median + (x - mean(d$x)) * beta2_median), from = 3, to = 7, ylim = c(3, 12), xlab = "", ylab = "") # 図9.6(B) plot(as.matrix(post.list)[,c("beta1", "beta2")]) ################################################################# # R2jags を使ってR-hatを算出する library(R2jags) init <- list(beta = c(0,0)) post.jags <- jags(data = list.data, inits = c(init, init, init), parameters.to.save = c("beta1", "beta2"), n.iter = 1600, model.file = filename, n.chains = 3, n.thin = 3, n.burnin = 100) print(post.jags) # R-hatが算出される