---
title: "Basic usage"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Basic usage}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
references:
- id: dienesmclatchie2018
title: Four reasons to prefer Bayesian analyses over significance testing
author:
- family: Dienes
given: Zoltan
- family: Mclatchie
given: Neil
container-title: Psychonomic Bulletin & Review
volume: 25
URL: https://link.springer.com/article/10.3758/s13423-017-1266-z
DOI: 10.3758/s13423-017-1266-z
page: 207–218
type: article-journal
issued:
year: 2018
- id: brandt2014
title: Does recalling moral behavior change the perception of brightness?
author:
- family: Brandt
given: Mark J
- family: IJzerman
given: Hans
- family: Blaken
given: Irene
container-title: Social Psychology
volume: 45
URL: https://econtent.hogrefe.com/doi/10.1027/1864-9335/a000191
DOI: 10.1027/1864-9335/a000191
page: 246–252
type: article-journal
issued:
year: 2014
- id: dienes2014
title: Using Bayes to get the most out of non-significant results
author:
- family: Dienes
given: Zoltan
container-title: Frontiers in Psychology
volume: 5
number: 781
URL: https://www.frontiersin.org/articles/10.3389/fpsyg.2014.00781/full
DOI: 10.3389/fpsyg.2014.00781
type: article-journal
issued:
year: 2014
- id: rouderbf
title: Bayesian t tests for accepting and rejecting the null hypothesis
author:
- family: Rouder
given: Jeffery N
- family: Speckman
given: Paul L
- family: Sun
given: Dongchu
- family: Morey
given: Richard D
container-title: Psychonomic Bulletin \& Review
volume: 20
number: 2
URL: http://dx.doi.org/10.3758/PBR.16.2.225
DOI: 10.3758/PBR.16.2.225
page: 225-237
type: article-journal
issued:
year: 2009
csl: apa.csl
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
To calculate a Bayes factor you need to 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 examples from @dienesmclatchie2018, @dienes2014, and from @rouderbf.
```{r setup}
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 `family`, and then the `mean` and
`se` parameters.
```{r}
data_mod <- likelihood(family = "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 `family` to `point`
and setting `point` as 0 (the default value).
```{r}
h0_mod <- prior(family = "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 `family` 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.
```{r}
h1_mod <- prior(family = "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.
```{r}
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.
```{r}
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.
```{r}
data_mod <- likelihood(family = "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.
```{r}
h1_mod <- prior(family = "uniform", min = 0, max = 20)
```
Following this, we specify a prior that describes the null hypothesis. Here,
@dienes2014 again uses a point null centred at 0.
```{r}
h0_mod <- prior(family = "point", point = 0)
```
This only thing left is to calculate the Bayes factor.
```{r}
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.
### Example 3
In Example three we'll reproduce an example from @rouderbf. @rouderbf specify
their models in terms of effect size units (*d*) rather than raw units as in
the example above. In this example by @rouderbf, they report a finding of a *t*
value of 2.03, with *n* of 80. To compute the Bayes factor, we first convert
this *t* value to a standardized effect size *d*. This *t* value equates to a
*d* of `r round(2.03 / sqrt(80),5)`. To model the data we use the
`noncentral_d` likelihood function, which is a rescaled noncentral *t*
distribution, with is parametrised in terms of *d* and *n*. We specify a null
model using a point prior at 0, and we specify the alternative model using a
Cauchy distribution centred at 0 (location parameter) with a scale parameter of
1.
```{r}
d <- 2.03 / sqrt(80) # convert t to d
data_model <- likelihood(family = "noncentral_d", d, 80)
h0_mod <- prior(family = "point", point = 0)
h1_mod <- prior(family = "cauchy", scale = 1)
bf <- integral(data_model * h0_mod) / integral(data_model * h1_mod)
bf
```
Performing the calculation as a above yields Bayes factor of ~`r round(bf,2)`,
the same value reported by @rouderbf.
To demonstrate the sensitivity of Bayes factor to prior specification,
@rouderbf recompute the Bayes factor for this example using a unit-information
(a standard normal) prior for the alternative model.
```{r}
d <- 2.03 / sqrt(80) # convert t to d
data_model <- likelihood(family = "noncentral_d", d, 80)
h0_mod <- prior(family = "point", point = 0)
h1_mod <- prior(family = "normal", mean = 0, sd = 1)
bf <- integral(data_model * h0_mod) / integral(data_model * h1_mod)
bf
```
Similarly recomputing our Bayes factor yields a value of ~`r round(bf,2)`, the
same value reported @rouderbf.
Although the Bayes factor outlined above is parametrised in terms of the effect
size *d*, it's also possible to performed the calculation directly on the *t*
statistic. To do this, however, we can't use the same Cauchy prior as above.
Instead, the Cauchy prior needs to be rescaled according to the same size. This
is because *t* values scale with sample size in a way that *d* values do not.
That is, for a given *d* the corresponding *t* value will be different
depending on the sample size. We can employ this alternative parametrisation in
the `Bayesplay` package by using the `noncentral_t` likelihood distribution.
The scale value for the Cauchy prior is now just multiplied by $\sqrt{n}$
```{r}
data_model <- likelihood(family = "noncentral_t", 2.03, 79)
h0_mod <- prior(family = "point", point = 0)
h1_mod <- prior(family = "cauchy", location = 0, scale = 1 * sqrt(80))
bf <- integral(data_model * h0_mod) / integral(data_model * h1_mod)
bf
```
Recomputing our Bayes factor now yields a value of ~`r round(bf,2)`, the same
value reported @rouderbf, and the same value reported above.
## References