-
Notifications
You must be signed in to change notification settings - Fork 0
/
bayes_factor.brms.Rmd
64 lines (52 loc) · 1.89 KB
/
bayes_factor.brms.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
---
title: "``brms::bayes_factor.brmsfit ``"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(brms)
```
Bayes Factors from Marginal Likelihoods
Description
Compute Bayes factors from marginal likelihoods.
Usage
<pre><code>
## S3 method for class 'brmsfit'
bayes_factor(x1, x2, log = FALSE, ...)
</code></pre>
#### Arguments
* x1 A brmsfit object
* x2 Another brmsfit object based on the same responses.
* log Report Bayes factors on the log-scale?
* ... Additional arguments passed to bridge_sampler.
#### Details
Computing the marginal likelihood requires samples of all variables defined in Stan’s parameters block to be saved. Otherwise bayes_factor cannot be computed. Thus, please set ``save_all_pars = TRUE``
in the call to brm, if you are planning to apply bayes_factor to your models.
The computation of Bayes factors based on bridge sampling requires a lot more posterior samples than usual. A good conservative rule of thump is perhaps 10-fold more samples (read: the default of 4000 samples may not be enough in many cases). If not enough posterior samples are provided, the bridge sampling algoritm tends to be unstable leading to considerably different results each time
it is run. We thus recommend running bayes_factor multiple times to check the stability of the results.
More details are provided under ``bridgesampling:bayes_factor``.
See Also
bridge_sampler, post_prob
#### Examples
```{r}
## Not run:
# model with the treatment effect
fit1 <- brm(
count ~ zAge + zBase + Trt,
data = epilepsy, family = negbinomial(),
prior = prior(normal(0, 1), class = b),
save_all_pars = TRUE
)
summary(fit1)
# model without the treatment effect
fit2 <- brm(
count ~ zAge + zBase,
data = epilepsy, family = negbinomial(),
prior = prior(normal(0, 1), class = b),
save_all_pars = TRUE
)
summary(fit2)
# compute the bayes factor
bayes_factor(fit1, fit2)
## End(Not run)
```