-
Notifications
You must be signed in to change notification settings - Fork 205
/
Copy pathdemo5_2.R
117 lines (102 loc) · 4.42 KB
/
demo5_2.R
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
#' ---
#' title: "Bayesian data analysis demo 5.2"
#' author: "Aki Vehtari, Markus Paasiniemi"
#' date: "`r format(Sys.Date())`"
#' output:
#' html_document:
#' theme: readable
#' code_download: true
#' ---
#' ## Hierarchical model for SAT-example data (BDA3, p. 120)
#'
#' ggplot2, grid, and gridExtra are used for plotting, tidyr for
#' manipulating data frames
#+ setup, message=FALSE, error=FALSE, warning=FALSE
library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)
library(latex2exp)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()
#' Data
y <- c(28,8,-3,7,-1,1,18,12)
s <- c(15,10,16,11,9,11,10,18)
#' Plot data, use normpdf to plot both the y_j and sigma_j
x <- seq(-40, 60, length.out = 500)
df_sep <- mapply(function(y, s, x) dnorm(x, y, s), y, s, MoreArgs = list(x = x)) %>%
as.data.frame() %>%
setNames(LETTERS[1:8]) %>%
cbind(x) %>%
pivot_longer(cols = !x, names_to = "school", values_to = "p")
labs1 <- c('Other Schools', 'School A')
plot_sep <- ggplot(data = df_sep) +
geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
labs(x = 'Treatment effect', y = '', title = 'Separate model', color = '') +
scale_y_continuous(breaks = NULL) +
scale_color_manual(values = c('blue','red'), labels = labs1) +
theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
#+ fig.height=4
plot_sep
#' A pooled model
df_pool <- data.frame(x = x, p = dnorm(x, sum(y/s^2)/sum(1/s^2), sqrt(1/sum(1/s^2))))
#' Create plot for the pooled model
plot_pool <- ggplot(data = df_pool) +
geom_line(aes(x = x, y = p, color = '1')) +
labs(x = 'Treatment effect', y = '', title = 'Pooled model', color = '') +
scale_y_continuous(breaks = NULL) +
scale_color_manual(values = 'red', labels = 'All schools') +
theme(legend.background = element_blank(), legend.position = c(0.7,0.9))
#' Load the pre-computed results for the hierarchical model
#' Replace this with your own code in the related exercise
load(root("demos_ch5","demo5_2.RData"))
#' Create plot for the hierarchical model
df_hier <- as.data.frame(t(pxm)) %>%
setNames(LETTERS[1:8]) %>%
cbind(x) %>%
pivot_longer(cols = !x, names_to = "school", values_to = "p")
plot_hier <- ggplot(data = df_hier) +
geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
labs(x = 'Treatment effect', y = '', title = 'Hierarchical model', color = '') +
scale_y_continuous(breaks = NULL) +
scale_color_manual(values = c('blue','red'), labels = labs1) +
theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
#' Plot separate, pooled, and hierarchical model
#+ fig.height=6
grid.arrange(plot_sep, plot_pool, plot_hier)
#' Various marginal and conditional posterior summaries
df_margpost = data.frame(x = t(tt), p = t(tp))
title1 <- TeX('Marginal posterior $p(\\tau | y)$')
plot_margpost <-
ggplot(data = df_margpost) +
geom_line(aes(x = x, y = p), color='forestgreen') +
labs(x = expression(tau), y = TeX('$p(\\tau | y)$'), title = title1) +
scale_y_continuous(breaks = 0) +
xlim(c(0, 35))
df_condmeans <- as.data.frame(t(tm)) %>%
setNames(LETTERS[1:8]) %>%
cbind(x = t(tt)) %>%
pivot_longer(cols = !x, names_to = "school", values_to = "p")
yl <- c(-5, 40)
title2 <- TeX('Conditional means E\\[$\\theta_j | \\tau , y $\\]')
plot_condmeans <- ggplot(data = df_condmeans) +
geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
coord_cartesian(ylim = yl, xlim = c(0, 35)) +
labs(x = expression(tau), y = TeX('E\\[$\\theta_j | \\tau , y $\\]'), title = title2, color = '') +
scale_color_manual(values = c('blue','red'), labels = labs1) +
theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
df_condsds <- as.data.frame(t(tsd)) %>%
setNames(LETTERS[1:8]) %>%
cbind(x = t(tt)) %>%
pivot_longer(cols = !x, names_to = "school", values_to = "p")
yl <- c(0, 25)
title3 <- TeX('Conditional standard deviations sd\\[$\\theta_j | \\tau , y $\\]')
plot_condsds <- ggplot(data = df_condsds) +
geom_line(aes(x = x, y = p, color = (school=='A'), group = school)) +
coord_cartesian(ylim = yl, xlim = c(0,35)) +
labs(x = expression(tau), y = TeX('sd\\[$\\theta_j | \\tau , y $\\]'), title = title3, color = '') +
scale_color_manual(values = c('blue','red'), labels = labs1) +
theme(legend.background = element_blank(), legend.position = c(0.8,0.9))
#' Plot the posterior summaries
#+ fig.height=8
grid.arrange(plot_margpost, plot_condmeans, plot_condsds)