-
Notifications
You must be signed in to change notification settings - Fork 206
/
rstan_demo.Rmd
483 lines (341 loc) · 14.7 KB
/
rstan_demo.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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
---
title: "Bayesian data analysis - RStan demos"
author: "Aki Vehtari, Markus Paasiniemi"
date: "First version 2017-07-17. Last modified `r format(Sys.Date())`."
output:
html_document:
fig_caption: yes
toc: TRUE
toc_depth: 2
number_sections: TRUE
toc_float:
smooth_scroll: FALSE
theme: readable
code_download: true
---
# Setup {.unnumbered}
```{r setup, include=FALSE}
knitr::opts_chunk$set(cache=FALSE, message=FALSE, error=FALSE, warning=TRUE, comment=NA, out.width='95%')
```
**Load packages**
```{r }
library(tidyr)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = 1)
library(loo)
library(ggplot2)
library(gridExtra)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(shinystan)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()
SEED <- 48927 # set random seed for reproducability
```
# Introduction
This notebook contains several examples of how to use [Stan](https://mc-stan.org) in R with __rstan__. This notebook assumes basic knowledge of Bayesian inference and MCMC. The Stan models are stored in separate .stan-files. The examples are related to [Bayesian data analysis course](https://avehtari.github.io/BDA_course_Aalto/).
Note that you can easily analyse Stan fit objects returned by `stan()` with a ShinyStan package by calling `launch_shinystan(fit)`.
# Bernoulli model
Toy data with sequence of failures (0) and successes (1). We would like to learn about the unknown probability of success.
```{r }
data_bern <- list(N = 10, y = c(1, 1, 1, 0, 1, 1, 1, 0, 1, 0))
```
Bernoulli model with a proper Beta(1,1) (uniform) prior
```{r }
code_bern <- root("demos_rstan", "bern.stan")
writeLines(readLines(code_bern))
```
Sample form the posterior and show the summary
```{r results='hide'}
fit_bern <- stan(file = code_bern, data = data_bern, seed = SEED)
monitor(fit_bern)
```
Plot the histogram of the posterior draws
```{r }
draws <- as.data.frame(fit_bern)
mcmc_hist(draws, pars='theta')
# or with base R
# hist(draws[,'theta'])
```
# Binomial model
Instead of sequence of 0's and 1's, we can summarize the data with the number of experiments and the number successes:
```{r }
data_bin <- list(N = 10, y = 7)
```
And then we use Binomial model with Beta(1,1) prior for the probability of success.
```{r }
code_binom <- root("demos_rstan","binom.stan")
writeLines(readLines(code_binom))
```
Sample from the posterior and plot the posterior. The histogram should look similar as in the Bernoulli case.
```{r results='hide'}
fit_bin <- stan(file = code_binom, data = data_bin, seed = SEED)
monitor(fit_bin)
draws <- as.data.frame(fit_bin)
mcmc_hist(draws, pars = 'theta')
```
Re-run the model with a new data. The compiled Stan program is re-used making the re-use faster.
```{r results='hide'}
data_bin <- list(N = 100, y = 70)
fit_bin <- stan(file = code_binom, data = data_bin, seed = SEED)
monitor(fit_bin)
draws <- as.data.frame(fit_bin)
mcmc_hist(draws, pars = 'theta')
```
## Explicit transformation of variables
In the above examples the probability of success $\theta$ was declared as
`real<lower=0,upper=1> theta;`
Stan makes automatic transformation of the variable to the unconstrained space using logit transofrmation for interval constrained and log transformation for half constraints.
The following example shows how we can also make an explicit transformation and use binomial_logit function which takes the unconstrained parameter as an argument and uses logit transformation internally. This form can be useful for better numerical stability.
```{r }
code_binomb <- root("demos_rstan", "binomb.stan")
writeLines(readLines(code_binomb))
```
Here we have used Gaussian prior in the unconstrained space, which produces close to uniform prior for theta.
Sample from the posterior and plot the posterior. The histogram should look similar as with the previous models.
```{r results='hide'}
data_bin <- list(N = 100, y = 70)
fit_bin <- stan(file = code_binomb, data = data_bin, seed = SEED)
monitor(fit_bin)
draws <- as.data.frame(fit_bin)
mcmc_hist(draws, pars = 'theta')
```
# Comparison of two groups with Binomial
An experiment was performed to estimate the effect of beta-blockers on mortality of cardiac patients. A group of patients were randomly assigned to treatment and control groups:
- out of 674 patients receiving the control, 39 died
- out of 680 receiving the treatment, 22 died
Data:
```{r }
data_bin2 <- list(N1 = 674, y1 = 39, N2 = 680, y2 = 22)
```
To analyse whether the treatment is useful, we can use Binomial model for both groups and compute odds-ratio:
```{r }
code_binom2 <- root("demos_rstan", "binom2.stan")
writeLines(readLines(code_binom2))
```
Sample from the posterior and plot the posterior
```{r results='hide'}
fit_bin2 <- stan(file = code_binom2, data = data_bin2, seed = SEED)
monitor(fit_bin2)
```
```{r warning=FALSE}
draws <- as.data.frame(fit_bin2)
mcmc_hist(draws, pars = 'oddsratio') +
geom_vline(xintercept = 1) +
scale_x_continuous(breaks = c(seq(0.25,1.5,by=0.25)))
```
# Linear Gaussian model
The following file has Kilpisjärvi summer month temperatures 1952-2013:
```{r }
data_kilpis <- read.delim(root("demos_rstan","kilpisjarvi-summer-temp.csv"), sep = ";")
data_lin <-list(N = nrow(data_kilpis),
x = data_kilpis$year,
xpred = 2016,
y = data_kilpis[,5])
```
Plot the data
```{r }
ggplot() +
geom_point(aes(x, y), data = data.frame(data_lin), size = 1) +
labs(y = 'Summer temp. @Kilpisjärvi', x= "Year") +
guides(linetype = "none")
```
To analyse whether the average summer month temperature is rising, we use a linear model with Gaussian model for the unexplained variation.
## Gaussian linear model with adjustable priors
The folloing Stan code allows also setting hyperparameter values as data allowing easier way to use different priors in different analyses:
```{r }
code_lin <- root("demos_rstan", "lin.stan")
writeLines(readLines(code_lin))
```
Create another list with data and priors
```{r }
data_lin_priors <- c(list(
pmualpha = mean(unlist(data_kilpis[,5])), # centered
psalpha = 100, # weakly informative
pmubeta = 0, # a priori incr. and decr. as likely
psbeta = (.1--.1)/6, # avg temp prob does does not incr. more than a degree per 10 years
pssigma = 1), # total variation in summer average temperatures is less +-3 degrees
data_lin)
```
Run Stan
```{r results='hide'}
fit_lin <- stan(file = code_lin, data = data_lin_priors, seed = SEED)
```
Stan gives a warning: There were {{get_num_max_treedepth(fit_lin)}} transitions after warmup that exceeded the maximum treedepth. You can use ShinyStan (`launch_shinystan(fit_lin)`) to look at the treedepth info and joint posterior of alpha and beta, to get a hint for the reason. ShinyStan helps also checking divergences, energy diagnostic, ESS's and Rhats.
Instead of interactive ShinyStan, we can also check the diagnostics as follows
```{r }
monitor(fit_lin)
```
The following diagnostics are explained in [CmdStan User's Guide](https://mc-stan.org/docs/cmdstan-guide/diagnose.html).
```{r message=TRUE}
check_hmc_diagnostics(fit_lin)
```
Compute the probability that the summer temperature is increasing.
```{r }
draws_lin <- rstan::extract(fit_lin, permuted = T)
mean(draws_lin$beta>0) # probability that beta > 0
```
Plot the data, the model fit and prediction for year 2016.
```{r }
mu <- apply(draws_lin$mu, 2, quantile, c(0.05, 0.5, 0.95)) %>%
t() %>% data.frame(x = data_lin$x, .) %>% gather(pct, y, -x)
pfit <- ggplot() +
geom_point(aes(x, y), data = data.frame(data_lin), size = 1) +
geom_line(aes(x, y, linetype = pct), data = mu, color = 'red') +
scale_linetype_manual(values = c(2,1,2)) +
labs(y = 'Summer temp. @Kilpisjärvi', x= "Year") +
guides(linetype = "none")
pars <- intersect(names(draws_lin), c('beta','sigma','ypred'))
draws <- as.data.frame(fit_lin)
phist <- mcmc_hist(draws, pars = pars)
grid.arrange(pfit, phist, nrow = 2)
```
## Gaussian linear model with standardized data
In the above we used the unnormalized data and as x values are far away from zero, this will cause very strong posterior dependency between alpha and beta (did you use ShinyStan for the above model?). The strong posterior dependency can be removed by normalizing the data to have zero mean. The following Stan code makes it in Stan. In generated quantities we do correspnding transformation back to the original scale.
```{r }
code_lin_std <- root("demos_rstan", "lin_std.stan")
writeLines(readLines(code_lin_std))
```
```{r results='hide'}
fit_lin_std <- stan(file = code_lin_std, data = data_lin, seed = SEED)
```
Now there were no warnings. You can use ShinyStan (`launch_shinystan(fit_lin)`) to look at the posterior and diagnostics and compare to the previous model results. We can also check diagnostics with the following commands.
```{r message=TRUE}
monitor(fit_lin_std)
check_hmc_diagnostics(fit_lin_std)
```
We see that there are no warnings by diagnostics and ESS's are higher than with the previous case with non-standardized data.
Next we check that we get similar probability for beta>0.
```{r }
draws_lin_std <- rstan::extract(fit_lin_std, permuted = T)
mean(draws_lin_std$beta>0) # probability that beta > 0
```
# Linear Student's $t$ model.
The temperatures used in the above analyses are averages over three months, which makes it more likely that they are normally distributed, but there can be extreme events in the feather and we can check whether more robust Student's $t$ observation model woul give different results.
```{r }
code_lin_std_t <- root("demos_rstan", "lin_std_t.stan")
writeLines(readLines(code_lin_std_t))
```
```{r results='hide'}
fit_lin_std_t <- stan(file = code_lin_std_t, data = data_lin, seed = SEED)
```
We get some warnings, but these specific warnings are not critical if counts are small as here.
Let's examine further diagnostics.
```{r message=TRUE}
monitor(fit_lin_std_t)
check_hmc_diagnostics(fit_lin_std_t)
```
We get similar diagnostics as for the linear Gaussian model with non-standardised data.
Compute the probability that the summer temperature is increasing.
```{r }
draws_lin_std_t <- extract(fit_lin_std_t, permuted = T)
mean(draws_lin_std_t$beta>0) # probability that beta > 0
```
We get similar probability as with Gaussian obervation model.
Plot data and the model fit
```{r }
mu <- apply(draws_lin_std_t$mu, 2, quantile, c(0.05, 0.5, 0.95)) %>%
t() %>% data.frame(x = data_lin$x, .) %>% gather(pct, y, -x)
pfit <- ggplot() +
geom_point(aes(x, y), data = data.frame(data_lin), size = 1) +
geom_line(aes(x, y, linetype = pct), data = mu, color = 'red') +
scale_linetype_manual(values = c(2,1,2)) +
labs(y = 'Summer temp. @Kilpisjärvi', x= "Year") +
guides(linetype = "none")
pars <- intersect(names(draws_lin_std_t), c('beta','sigma','nu','ypred'))
draws <- as.data.frame(fit_lin_std_t)
phist <- mcmc_hist(draws, pars = pars)
grid.arrange(pfit, phist, nrow = 2)
```
We see also that the marginal posterior of nu is wide with lot of mass for values producing distrbution really close to Gaussian.
# Pareto-smoothed importance-sampling leave-one-out cross-validation (PSIS-LOO)
We can use leave-one-out cross-validation to compare the expected predictive performance. For the following lines to work, the log-likelihood needs to be evaluated in the stan code. For an example, see lin.stan and [Computing approximate leave-one-out cross-validation usig PSIS-LOO](http://mc-stan.org/loo/articles/loo2-with-rstan.html).
```{r }
loo_lin_std <- loo(fit_lin_std)
loo_lin_std_t <- loo(fit_lin_std_t)
loo_compare(loo_lin_std, loo_lin_std_t)
```
There is no practical difference between Gaussian and Student's $t$ observation model for this data.
# Comparison of $k$ groups with hierarchical models
Let's compare the temperatures in three summer months.
```{r }
data_kilpis <- read.delim(root("demos_rstan","kilpisjarvi-summer-temp.csv"), sep = ";")
data_grp <-list(N = 3*nrow(data_kilpis),
K = 3,
x = rep(1:3, nrow(data_kilpis)),
y = c(t(data_kilpis[,2:4])))
```
## Common variance (ANOVA) model
```{r }
code_grp_aov <- root("demos_rstan", "grp_aov.stan")
writeLines(readLines(code_grp_aov))
```
Fit the model
```{r results='hide'}
fit_grp <- stan(file = code_grp_aov, data = data_grp, seed = SEED)
monitor(fit_grp)
check_hmc_diagnostics(fit_grp)
```
## Common variance and hierarchical prior for mean.
Results do not differ much from the previous, because there is only
few groups and quite much data per group, but this works as an example of a hierarchical model.
```{r }
code_grp_prior_mean <- root("demos_rstan", "grp_prior_mean.stan")
writeLines(readLines(code_grp_prior_mean))
```
Fit the model
```{r results='hide'}
fit_grp <- stan(file = code_grp_prior_mean, data = data_grp, seed = SEED)
monitor(fit_grp)
check_hmc_diagnostics(fit_grp)
```
We got a small number of divergences, so we increase
`adapt_delta=0.95`. Note that thisis useful only in case of small
number of divergences, and increasing `adapt_delta>0.99` doesn't
usually make sense, and in such cases there is need to investigate the
identifiability or parameter transformations.
```{r results='hide'}
fit_grp <- stan(file = code_grp_prior_mean, data = data_grp, seed = SEED, control=list(adapt_delta=0.95));
monitor(fit_grp)
check_hmc_diagnostics(fit_grp)
```
## Unequal variance and hierarchical prior for mean and variance
```{r }
code_grp_prior_mean_var <- root("demos_rstan", "grp_prior_mean_var.stan")
writeLines(readLines(code_grp_prior_mean_var))
```
Fit the model
```{r results='hide'}
fit_grp <- stan(file = code_grp_prior_mean_var, data = data_grp, seed = SEED)
monitor(fit_grp)
check_hmc_diagnostics(fit_grp)
```
We got a small number of divergences, so we increase
`adapt_delta=0.95`. Note that thisis useful only in case of small
number of divergences, and increasing `adapt_delta>0.99` doesn't
usually make sense, and in such cases there is need to investigate the
identifiability or parameter transformations.
```{r results='hide'}
fit_grp <- stan(file = code_grp_prior_mean_var, data = data_grp, seed = SEED, control=list(adapt_delta=0.95));
monitor(fit_grp)
check_hmc_diagnostics(fit_grp)
```
Plot the results
```{r }
draws_grp <- extract(fit_grp, permuted = T)
temps <- data.frame(draws_grp$mu) %>%
setNames(c('June','July','August'))
mcmc_areas(temps) + xlab('Temperature')
```
Probabilities that June is hotter than July, June is hotter than August
and July is hotter than August:
```{r }
paste('p(TempJun > TempJul) = ', mean(temps$June > temps$July))
paste('p(TempJun > TempAug) = ', mean(temps$June > temps$August))
paste('p(TempJul > TempAug) = ', mean(temps$July > temps$August))
```
<br />
# Licenses {.unnumbered}
* Code © 2017-2020, Aki Vehtari, 2017 Markus Paasiniemi, licensed under BSD-3.
* Text © 2017-2020, Aki Vehtari, licensed under CC-BY-NC 4.0.