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

glmmstanとbrmsで研究データを解析しています。

waicを比較する方法を調べたので、忘れないうちにメモです。

以下すべて太字がwaicの比較に必要なコードです。

間違っていたり、もっと良いやり方があったら、ぜひご指摘ください。

◯

glmmstanでwaicの比較は、どうしたらよいかわからなかったので、ちょっと調べたら以下のように記述すると出力できました。

looパッケージを呼び出して比較しています。

library(glmmstan) <strong>library(loo)</strong>

1 2	library(glmmstan) <strong>library(loo)</strong>

fit0 &lt;- Pglmmstan(a ~ b + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)

1	fit0 <- Pglmmstan(a ~ b + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)

fit1 &lt;- Pglmmstan(a ~ b + c + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)

1	fit1 <- Pglmmstan(a ~ b + c + (1|d), data = dat, chains = 4, tier = 5000, family = “gaussian”)

<strong>log_lik1 &lt;- extract_log_lik(fit0) waic1 &lt;- waic(log_lik1) print(waic1) log_lik2 &lt;- extract_log_lik(fit1) waic2 &lt;- waic(log_lik2) print(waic2) compare(waic1, waic2)</strong>

1 2 3 4 5 6 7

<strong>log_lik1 &lt;- extract_log_lik(fit0) waic1 &lt;- waic(log_lik1) print(waic1) log_lik2 &lt;- extract_log_lik(fit1) waic2 &lt;- waic(log_lik2) print(waic2) compare(waic1, waic2)</strong>

結果は、以下のようにelpd_waicの差で示されます。

&gt; log_lik1 &lt;- extract_log_lik(fit0) &gt; waic1 &lt;- waic(log_lik1) &gt; print(waic1) Computed from 10000 by 86 log-likelihood matrix

1 2 3 4

> 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 &gt; log_lik2 &lt;- extract_log_lik(fit1) &gt; waic2 &lt;- waic(log_lik2) &gt; print(waic2) Computed from 10000 by 86 log-likelihood matrix

1 2 3 4 5 6 7 8

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 &gt; compare(waic1, waic2)elpd_diff se weight1 weight2 19.1 2.0 0.0 1.0

1 2 3 4 5 6

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の比較ができます。

library(brms)

1	library(brms)

rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())

1 2	rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())

fit0 &lt;- brm(a ~ b + (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)

1	fit0 <- brm(a ~ b + (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)

fit1 &lt;- brm(a ~ b + c+ (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)

1	fit1 <- brm(a ~ b + c+ (1|d), data = dat, n.chains = 4, n.iter = 2000, n.warmup = 1000, n.cluster = 2, family = “gaussian”)

<strong>WAIC(fit0) WAIC(fit1) WAIC(fit0,fit1)</strong>

1 2 3

<strong>WAIC(fit0) WAIC(fit1) WAIC(fit0,fit1)</strong>

結果は、以下のようにwaicの差で示されます。

&gt; WAIC(fit0) WAIC SE 316.48 7.3 &gt; WAIC(fit1) WAIC SE 301.9 8.91 &gt; WAIC(fit0,fit1) WAIC SE Weights fit0 316.48 7.30 0 fit1 301.90 8.91 1 fit0 - fit1 14.57 3.89

1 2 3 4 5 6 7 8 9

> 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())

1 2	rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())

◯

その他、参考になる資料は以下です。

brmsパッケージ