【解説】とても便利なPubmedのメールアラート機能
Makoto KYOUGOKU
京極真のWebマガジン
[備忘録]glmmstanとbrmsでwaicの比較
Makoto KYOUGOKU
記事内にプロモーションを含む場合があります
本記事では「glmmstanとbrmsでwaicの比較」についてサクッと解説しています。
Contents
[備忘録]glmmstanとbrmsでwaicの比較
glmmstanとbrmsで研究データを解析しています。
waicを比較する方法を調べたので、忘れないうちにメモです。
以下すべて太字がwaicの比較に必要なコードです。
間違っていたり、もっと良いやり方があったら、ぜひご指摘ください。
glmmstanでwaicの比較
glmmstanでwaicの比較は、どうしたらよいかわからなかったので、ちょっと調べたら以下のように記述すると出力できました。
looパッケージを呼び出して比較しています。
library(glmmstan)
library(loo)
fit0 <- Pglmmstan(a ~ b + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)
fit1 <- Pglmmstan(a ~ b + c + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)
log_lik1 <- extract_log_lik(fit0)
waic1 <- waic(log_lik1)
print(waic1)
log_lik2 <- extract_log_lik(fit1)
waic2 <- waic(log_lik2)
print(waic2)
compare(waic1, waic2)
結果は、以下のようにelpd_waicの差で示されます。
> log_lik1 <- extract_log_lik(fit0)
> waic1 <- waic(log_lik1)
> print(waic1)
Computed from 10000 by 86 log-likelihood matrix
Estimate SE
elpd_waic -158.5 3.6
p_waic 4.8 0.5
waic 317.0 7.3
> log_lik2 <- extract_log_lik(fit1)
> waic2 <- waic(log_lik2)
> print(waic2)
Computed from 10000 by 86 log-likelihood matrix
Estimate SE
elpd_waic -139.4 4.4
p_waic 53.8 3.4
waic 278.9 8.8
> compare(waic1, waic2)elpd_diff se weight1 weight2 19.1 2.0 0.0 1.0
ちなみにglmmstanの並列化は、Pglmmstanと書くだけでOKです。
brmsでwaicの比較
brmsはパッケージに書いてある通り、以下のように記述するとwaicの比較ができます。
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit0 <- brm(a ~ b + (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)
fit1 <- brm(a ~ b + c+ (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)
WAIC(fit0)
WAIC(fit1)
WAIC(fit0,fit1)
結果は、以下のようにwaicの差で示されます。
> WAIC(fit0)
WAIC SE
316.48 7.3
> WAIC(fit1) WAIC SE
301.9 8.91
> WAIC(fit0,fit1) WAIC SE Weights
fit0 316.48 7.30 0
fit1 301.90 8.91 1
fit0 – fit1 14.57 3.89
ちなみに、library(brms)の直下にある以下のコードは、並列化するために書いています。
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
その他、参考になる資料は以下です。
brmsパッケージ
glmmstanパッケージを作ってみた from Hiroshi Shimizu
rstanで簡単にGLMMができるglmmstan()を作ってみた from Hiroshi Shimizu
階層ベイズとWAIC from Hiroshi Shimizu
まとめ：【備忘録】glmmstanとbrmsでwaicの比較
本記事では「glmmstanとbrmsでwaicの比較」についてサクッと解説しました。
役に立てば嬉しいです。
著者紹介
1976年大阪府生まれ。Ph.D、OT。Thriver Project代表。首都大学東京大学院人間健康科学研究科博士後期課程・終了。吉備国際大学ならびに同大学大学院・教授（役職：人間科学部長、保健科学研究科長）。『この一冊でわかる！セラピストのための研究論文の書き方ガイド』（三輪書店）、『医療関係者のための信念対立解明アプローチ』（誠信書房）、『OCP・OFP・OBPで学ぶ作業療法実践の教科書』（メジカルビュー）、『作業で創るエビデンス』（医学書院）など著書・論文多数。
