In BayesianPlayground/bayesplay: The Bayes Factor Playground

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

To calculate a Bayes factor you specify three things:

1. A likelihood that provides a description of the data

2. A prior that specifies the predictions of the first model to be compared

3. And a prior that specifies the predictions of the second model to be compared.

By convention, the two models to be compared are usually called the null and the alternative models.

The Bayes factor is defined as the ratio of two (weighted) average likelihoods where the prior provides the weights. Mathematically, the weighted average likelihood is given by the following integral

$$\mathcal{M}H = \int{\theta\in\Theta_H}\mathcal{l}_H(\theta|\mathbf{y})p(\theta)d\theta$$

Where $\mathcal{l}_H(\theta|\mathbf{y})$ represents the likelihood function, $p(\theta)$ represents the prior on the parameter, with the integral defined over the parameter space of the hypothesis ($\Theta_H$).

To demonstrate how to calculate Bayes factors using bayesplay, we can reproduce an example from @dienesmclatchie2018 and

library(bayesplay)

Examples

Example 1

The first example from @dienesmclatchie2018 that we'll reproduce describes a study from @brandt2014. In the study by @brandt2014, they obtained a mean difference for 5.5 Watts (t statistic = 0.17, SE = r round(5.5 / 0.17,2)).

We can describe this observation using a normal likelihood using the likelihood() function. We first specify the distribution, and then the mean and se parameters.

data_mod <- likelihood(distribution = "normal", mean = 5.5, sd = 32.35)

Following this, @dienesmclatchie2018 describe the two models they intend to compare. First, the null model is described as a point prior centred at 0. We can specify this with the prior() function, setting distribution to point and setting point as 0 (the default value).

h0_mod <- prior(distribution = "point", point = 0)

Next, @dienesmclatchie2018 describe the alternative model. For this they use a half-normal distribution with a mean of 0 and a standard deviation of 13.3. This can again be specified using the prior() function setting distribution to normal and setting the mean and sd parameters as required. Additionally, because they specify a half-normal distribution, an additional range value is needed to restrict the parameter range to 0 (the mean) to positive infinity.

h1_mod <- prior(distribution = "normal", mean = 0, sd = 13.3, range = c(0, Inf))

With the three parts specified we can compute the Bayes factor. Following the equation above, the first step is to calculate $\mathcal{M}_H$ for each model. To do this, we multiply the likelihood by the prior and integrate.

m1 <- integral(data_mod * h1_mod)
m0 <- integral(data_mod * h0_mod)

With $\mathcal{M}{H_1}$ and $\mathcal{M}{H_0}$ calculated, we now simply divide the two values to obtain the Bayes factor.

bf <- m1 / m0
bf

This gives a Bayes factor of ~r round(bf,2), the same value reported by @dienesmclatchie2018.

Example 2

The second example, from @dienes2014, we'll reproduce relates to an experiment where a mean difference of 5% was observed with a standard error of 10%. We can describe this observation using a normal likelihood.

data_mod <- likelihood(distribution = "normal", mean = 5, sd = 10)

Next, we specify a prior which described the alternative hypothesis. In this case, @dienes2014 uses a uniform prior that ranges from 0 to 20.

h1_mod <- prior(distribution = "uniform", min = 0, max = 20)

Following this, we specify a prior that describes the null hypothesis. Here, @dienes2014 again uses a point null centered at 0.

h0_mod <- prior(distribution = "point", point = 0)

This only thing left is to calculate the Bayes factor.

bf <- integral(data_mod * h1_mod) / integral(data_mod * h0_mod)
bf

This gives a Bayes factor of ~r round(bf,2), the same value reported by @dienes2014.

References

BayesianPlayground/bayesplay documentation built on Jan. 8, 2021, 12:42 a.m.