3_First_Stan_examples.Rmd

---
title: "3. First Stan examples"
output: html_notebook
---


## Introduction

Stan is a platform for running statistical computations, especially Bayesian data analysis. The mechanism by which it works (Hamiltonian Markov Chain Monte Carlo) is sophisticated. Here we'll describe how to set up a model and compute a posterior distribution.

The first time you run the code chunks in this file, it will take a while. Be patient. I've included options and code that will allow Stan to bypass a lot of the work R has to do the second time around (and every time thereafter, as long as you're in the same R session).


## Preliminaries

Load necessary packages:

```{r}
library(tidyverse)
library(rstan)
library(tidybayes)
library(bayesplot)
```

Stan will take a fair amount of time the first time it's run on a model. The following bit of code tells Stan not to recompile code that has already been compiled once so that it won't take so long after the first time.

```{r}
rstan_options(auto_write = TRUE)
```

Finally, the following code will detect if there are multiple cores available for parallel processing, and if so, use those cores to speed up the sampling process.

```{r}
options(mc.cores = parallel::detectCores())
```


## Stan's structure

Stan code looks quite different from R code. (Stan is actually coded in C.) We'll see later, for example, that each line must terminate with a semicolon, and comments use two forward slashes.

There are three basic code blocks in a simple Stan program. The first is called `data` and it reads in the raw data. The second is `parameters` in which you tell Stan the parameters for which you're conducting statistical inference. The third is `model` where you describe the statistical model.

There are other, more advanced possibilities for code blocks, but those three suffice for a complete model.


## Storing the data

We re-create the analysis of 12 successes in 18 trials from the `2_Continuous_Bayes.Rmd` notes.

Stan is a little particular about the form of the data. It requires a `list` from R. Here's how we do that. First, just store the values we need.

```{r}
N1 <- 18  # Define the sample size
y1 <- c(rep(1, 12), rep(0, 6))  # 12 successes and 6 failures
```

The line defining `y` is just using an R trick to get 12 ones and 6 zeros. Observe:

```{r}
y1
```

Now we bundle `N` and `y` together into a list.

```{r}
stan_data <- list(N = N1, y = y1)
stan_data
```

This is the format required by Stan for accessing the data.


## Uniform prior

We'll start with assuming a uniform prior.

### The Stan code

R Markdown is awesome, so it gives us a way to create a Stan model within a code chunk. The resulting model is stored in the variable name specified by the `output.var` chunk option (in this case, `bin_unif`).

We also use the option `cache = TRUE` in the header of the code chunk. This chunk takes a while the first time it's run; however, once the result is cached, it won't run again even when you "Run All" code chunks. (You may also ignore any warnings that may appear below the chunk once it's finished running.)


```{stan, output.var = "bin_unif", cache = TRUE}
data {
    int<lower = 0> N;
    array[N] int<lower = 0, upper = 1> y;
}
parameters {
    real<lower = 0, upper = 1> theta;
}
model {
    theta ~ uniform(0, 1);  // prior
    y ~ bernoulli(theta);   // likelihood
}
```


Let's take a closer look at various pieces of the Stan code. Here is the `data` block:

```
data {
    int<lower = 0> N;
    array[N] int<lower = 0, upper = 1> y;
}
```

Notice that the two variables declared here are precisely the variables in the list `stan_data` (namely, `N`, and `y`). The `lower` and/or `upper` bounds in the definitions just serve as a check to make sure we are passing in data of the right form. (If invalid data is passed to Stan, the program will stop with an error.) The notation `array[N]` means that the variable `y` should be an array of `N` integers.

The `parameters` block is

```
parameters {
    real<lower = 0, upper = 1> theta;
}
```

This simply declares that we are interested in doing inference for a value called `theta`. The `lower` and `upper` bounds here serve a slightly different purpose. These tell Stan that it shouldn't consider values of `theta` unless they lie between 0 and 1.

The `model` block is

```
model {
    theta ~ uniform(0, 1);  // prior
    y ~ bernoulli(theta);   // likelihood
}
```

The first line is saying that we start with a uniform prior on `theta`. Technically, this line is not required: the parameter declaration `real<lower = 0, upper = 1> theta` defines a uniform $(0, 1)$ prior on `theta` by default. But we include it for completeness because in most situations, we won't use a uniform prior.

The second line is the likelihood function. An event that consists of a success/failure trial is called a "Bernoulli" trial, so `bernoulli(theta)` describes a probability model that considers each of the `y` values as being generated by a Bernoulli trial with probability `theta`. It is, of course, the true value of `theta` that we are trying to determine with our data.

### Sampling from the model

Now we sample from the model using our data. (Since there are random processes in play, we'll use the `seed` argument to make our work reproducible.)

```{r, cache = TRUE}
fit_bin_unif <- sampling(bin_unif,
                         data = stan_data,
                         seed = 54321)
```

There is a lot of output here, showing the progress of running each chain.

### Summarizing the model

Here are some summary statistics for the posterior distribution.

```{r}
fit_bin_unif
```

Ignore the line that starts with `lp__`; we only care about the values of `theta`. With a uniform prior, remember that the posterior is identical to the (scaled) likelihood function. Since 12 out of 18 is $2/3$, the mean of 0.65 makes sense. Of course, there is a range of values for `theta` that is consistent with obtaining 12 successes in 18 trials, and that's represented in the table using the standard deviation `sd` and the percentiles to its right. (Don't worry about `se_mean`, `n_eff`, or `Rhat` for now.)

### Visualizing the model

We can also make plots of the posterior distribution.

```{r}
mcmc_hist(fit_bin_unif, pars = "theta")
```

The following plot is called a "density plot". It is a smooth approximation of the histogram above. In other words, it's our best guess as to the shape of the function representing the posterior probability density.

```{r}
mcmc_dens(fit_bin_unif, pars = "theta")
```

Although it's a little rough in shape due to the fact that the posterior is simulated, you can see that it's close to the theoretically correct posterior shown in the `2_Continuous_Bayes.Rmd` notes. It will be easier to see if we put it on the same x-axis scale as it was there:

```{r}
mcmc_dens(fit_bin_unif, pars = "theta") +
    xlim(0,1)
```

Let's try some other plotting options.

The default option for `mcmc_intervals` gives 50% and 90% credible intervals for our parameters. We would likely need a much higher "effective sample size" to get good estimates of these, especially at the 95% level, but at least we can see the code we need.

```{r}
mcmc_intervals(fit_bin_unif, pars = "theta")
```

The following is a "traceplot" showing the chains. You use this to diagnose chain mixing.

```{r}
mcmc_trace(fit_bin_unif, pars = "theta")
```


## Normal prior (far from data)

The next example in `2_Continuous_Bayes.Rmd` is about selecting the normal prior

$$
\theta \sim N(0.3, 0.1).
$$

Here is the Stan code used to define this model. The only change is the line
```
theta ~ normal(0.3, 0.1);  // prior
```
in the `model` block.

```{stan, output.var = "bin_norm1", cache = TRUE}
data {
    int<lower = 0> N;
    array[N] int<lower = 0, upper = 1> y;
}
parameters {
    real<lower = 0, upper = 1> theta;
}
model {
    theta ~ normal(0.3, 0.1);   // prior
    y ~ bernoulli(theta);       // likelihood
}
```

Once again, we sample from the model using our data. We'll also tell the `sampling` function not to print out all the information about processing the chains by using the argument `refresh = 0`.

```{r, cache = TRUE}
fit_bin_norm1 <- sampling(bin_norm1,
                          data = stan_data,
                          seed = 54321,
                          refresh = 0)
```

```{r}
fit_bin_norm1
```

```{r}
mcmc_dens(fit_bin_norm1, pars = "theta") +
    xlim(0, 1)
```

Again, compare this to the theoretical posterior in the `2_Continuous_Bayes.Rmd` notes: the mode is a little to the left of 0.5 and most of the distribution is concentrated between about 0.25 to 0.7 or so.


## Normal prior (close to data)

Suppose the prior is

$$\theta \sim N(0.7, 0.1).$$

In the code chunk below, we make the necessary change.

```{stan, output.var = "bin_norm2", cache = TRUE}
data {
    int<lower = 0> N;
    array[N] int<lower = 0, upper = 1> y;
}
parameters {
    real<lower = 0, upper = 1> theta;
}
model {
    theta ~ normal(0.7, 0.1);   // prior
    y ~ bernoulli(theta);       // likelihood
}
```

```{r, cache = TRUE}
set.seed(54321)
fit_bin_norm2 <- sampling(bin_norm2,
                          data = stan_data,
                          refresh = 0)
```

```{r}
fit_bin_norm2
```

```{r}
mcmc_dens(fit_bin_norm2, pars = "theta") +
    xlim(0,1)
```


## Triangular prior

Stan does not have a built-in expression for a triangular distribution, so we'll leave that example out. (There is a hacky way to make it work, but it's not really worth it. You don't really want to use a triangular distribution for real-world problems anyway.)


## Change the data

The `2_Continuous_Bayes.Rmd` notes also show the effect of the prior in situations of sparse data. What if we have only 2 successes in 3 trials?

```{r}
N2 <- 3  # Define the sample size
y2 <- c(1, 1, 0)  # 2 successes and 1 failure
stan_data_small <- list(N = N2, y = y2)
```

Let's apply the three models already developed to this new data. Note that we do not need to run any more Stan code. This is just sampling from existing models using new data.

```{r, cache = TRUE}
fit_bin_unif_small <- sampling(bin_unif,
                               data = stan_data_small,
                               seed = 54321,
                               refresh = 0)
fit_bin_norm1_small <- sampling(bin_norm1,
                                data = stan_data_small,
                                seed = 54321,
                                refresh = 0)
fit_bin_norm2_small <- sampling(bin_norm2,
                                data = stan_data_small,
                                seed = 54321,
                                refresh = 0)
```

Summaries and plots to follow. Compare these to the ones shows in the `2_Continuous_Bayes.Rmd` notes.

```{r}
fit_bin_unif_small
```

```{r}
mcmc_dens(fit_bin_unif_small, pars = "theta") +
    xlim(0,1)
```

```{r}
fit_bin_norm1_small
```

```{r}
mcmc_dens(fit_bin_norm1_small, pars = "theta") +
    xlim(0,1)
```

```{r}
fit_bin_norm2_small
```

```{r}
mcmc_dens(fit_bin_norm2_small, pars = "theta") +
    xlim(0,1)
```