Hatena::ブログ(Diary)

盆栽日記

2012-05-22

緑本メモ

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が算出される

はてなユーザーのみコメントできます。はてなへログインもしくは新規登録をおこなってください。

トラックバック - http://d.hatena.ne.jp/dichika/20120522/1337681689
リンク元