**Load packages**

```
library(tidyr)
library(dplyr)
library(rstan)
library(rstanarm)
options(mc.cores = 1)
library(loo)
library(shinystan)
library(ggplot2)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(ggdist)
library(gridExtra)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()
SEED <- 48927 # set random seed for reproducability
```

This notebook contains several examples of how to use Stan in R with **rstanarm**. This notebook assumes basic knowledge of Bayesian inference and MCMC. The examples are related to Bayesian data analysis course.

Note that you can easily analyse Stan fit objects returned by `stan_glm()`

with a ShinyStan package by calling `launch_shinystan(fit)`

.

The models are not exactly equal to the models at rstan_demo.Rmd, but rather serve as examples of how to implement similar models with **rstanarm**.

Toy data with sequence of failures (0) and successes (1). We would like to learn about the unknown probability of success.

`data_bern <- data.frame(y = c(1, 1, 1, 0, 1, 1, 1, 0, 1, 0))`

Uniform prior (beta(1,1)) is achieved by setting the prior to NULL, which is not recommended in general. y ~ 1 means y depends only on the intercept term

```
fit_bern <- stan_glm(y ~ 1, family = binomial(), data = data_bern,
prior_intercet = NULL, seed = SEED, refresh = 0)
```

You can use ShinyStan examine and diagnose the fitted model is to call shinystan in R terminal as `launch_shinystan(fit_bern)`

Monitor provides summary statistics and diagnostics

`monitor(fit_bern$stanfit)`

```
Inference for the input samples (4 chains: each with iter = 2000; warmup = 0):
Q5 Q50 Q95 Mean SD Rhat Bulk_ESS Tail_ESS
(Intercept) -0.2 0.8 2.1 0.9 0.7 1.01 1449 1151
mean_PPD 0.3 0.7 1.0 0.7 0.2 1.00 2243 4000
log-posterior -10.0 -8.3 -8.0 -8.5 0.8 1.01 1373 1268
For each parameter, Bulk_ESS and Tail_ESS are crude measures of
effective sample size for bulk and tail quantities respectively (an ESS > 100
per chain is considered good), and Rhat is the potential scale reduction
factor on rank normalized split chains (at convergence, Rhat <= 1.05).
```

To see the parameter values on the ouput space, do the inverse logistic transformation (plogis in R) on the intercept

```
draws <- as.data.frame(fit_bern)
mean(draws$`(Intercept)`)
```

`[1] 0.8606056`

Probability of success

```
draws$theta <- plogis(draws$`(Intercept)`)
mean(draws$theta)
```

`[1] 0.683932`

Histogram of theta

`mcmc_hist(draws, pars='theta') + xlab('theta')`