-
Notifications
You must be signed in to change notification settings - Fork 206
/
demo2_3.Rmd
60 lines (48 loc) · 1.48 KB
/
demo2_3.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
---
title: "Bayesian data analysis demo 2.3"
author: "Aki Vehtari, Markus Paasiniemi"
date: "`r format(Sys.Date())`"
output:
html_document:
theme: readable
code_download: true
---
## Probability of a girl birth given placenta previa (BDA3 p. 37).
Simulate samples from Beta(438,544), draw a histogram with
quantiles, and do the same for a transformed variable.
ggplot2 is used for plotting, tidyr for manipulating data frames
```{r setup, message=FALSE, error=FALSE, warning=FALSE}
library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
```
Sample from posterior Beta(438,544).
Obtain all draws at once and store them in vector 'theta'
```{r }
a <- 438
b <- 544
theta <- rbeta(10000, a, b)
```
Compute odds ratio for all draws
```{r }
phi <- (1 - theta) / theta
```
Compute 2.5% and 97.5% quantile approximation using samples
```{r }
quantiles <- c(0.025, 0.975)
thetaq <- quantile(theta, quantiles)
phiq <- quantile(phi, quantiles)
```
Histogram plots with 30 bins for theta and phi
```{r }
# merge the data into one data frame for plotting
df1 <- data.frame(phi,theta) %>% pivot_longer(everything())
# merge quantiles into one data frame for plotting
df2 <- data.frame(phi = phiq, theta = thetaq) %>% pivot_longer(everything())
ggplot() +
geom_histogram(data = df1, aes(value), bins = 30) +
geom_vline(data = df2, aes(xintercept = value), linetype = 'dotted') +
facet_wrap(~name, ncol = 1, scales = 'free_x') +
labs(x = '', y = '') +
scale_y_continuous(breaks = NULL)
```