オンライン教育プラットフォーム|Thriver Project
生活・仕事
PR

[備忘録]glmmstanとbrmsでwaicの比較

Makoto KYOUGOKU
記事内にプロモーションを含む場合があります
京極真
京極真

本記事では「glmmstanとbrmsでwaicの比較」についてサクッと解説しています。

[備忘録]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代表。吉備国際大学ならびに同大学大学院・教授。作業療法学科長、保健科学研究科長、(通信制)保健科学研究科長。首都大学東京大学院人間健康科学研究科博士後期課程・終了。『医療関係者のための信念対立解明アプローチ』『作業療法リーズニングの教科書』『作業で創るエビデンス』など著書・論文多数。
記事URLをコピーしました