/
robustness.Rmd
118 lines (91 loc) · 2.94 KB
/
robustness.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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
---
title: "Robustness regions"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Robustness regions}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
```{r setup}
library(bayesplay)
```
In typical Bayes factor analyses, a prior is specified for the alternative and
null models. For example, an analyst might specify a normal distribution with a
mean of 0 and standard deviation of 15 as the prior for the alternative model.
In `bayesplay` syntax, this would be done as follows:
```{r, eval=FALSE}
prior(family = "normal", mean = 0, sd = 15)
```
In specifying a prior in this way, the analyst chooses a set value for the two
parameters (mean and standard deviation). It might be of interest to
The `bayesplay` package includes functionality for systematically varying
parameters of the alternative prior for computing *robustness regions* and for
visualizing these robustness regions.
```{r}
data_model <- likelihood(family = "student_t", mean = 5, sd = 11, df = 12)
alternative_prior <- prior(family = "normal", mean = 0, sd = 15)
null_prior <- prior(family = "point", point = 0)
m1 <- data_model * alternative_prior
m0 <- data_model * null_prior
bf <- integral(m1) / integral(m0)
summary(bf)
```
One might then wonder how robust this conclusion is to variations in the
alternative prior. Would the conclusions change if a higher value or a lower
value was selected for the standard deviation of the prior? With the `bfrr`
function, it's possible to set an upper and lower bound of the standard deviation
and to compute and visualize the resulting Bayes factor for values within this
range.
In the example below, the standard deviation is varied across the range from 1
to 20 (in steps of 1).
```{r}
rr <- bfrr(
likelihood = data_model, # reuse the definition above
alternative_prior = alternative_prior, # reuse the definition above
null_prior = null_prior, # reuse the definition above
parameters = list(sd = c(1, 50)) # vary sd from 1 to 20
)
```
```{r}
plot(rr)
```
```{r}
data_model <- likelihood(family = "noncentral_d", d = 0.2, n = 50)
alternative_prior <- prior(family = "cauchy", location = 0, scale = 1)
null_prior <- prior(family = "point", point = 0)
m1 <- data_model * alternative_prior
m0 <- data_model * null_prior
bf <- integral(m1) / integral(m0)
summary(bf)
```
```{r}
rr <- bfrr(
likelihood = data_model, # reuse the definition above
alternative_prior = alternative_prior, # reuse the definition above
null_prior = null_prior, # reuse the definition above
parameters = list(scale = c(0.1, 2))
)
rr
```
```{r}
plot(rr)
```
```{r}
rr <- bfrr(
likelihood = data_model, # reuse the definition above
alternative_prior = alternative_prior, # reuse the definition above
null_prior = null_prior, # reuse the definition above
parameters = list(scale = c(0, 2), location = c(-1, 1)),
steps = 32
)
rr
```
```{r}
plot(rr)
```