---
title: "Plotting MCMC draws using the bayesplot package"
author: "Jonah Gabry"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Plotting MCMC draws}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, settings, include=FALSE}
library(bayesplot)
knitr::opts_chunk$set(
  dev = "pdf",
  fig.align = "center",
  fig.width = 6,
  fig.height = 4,
  comment = NA
)
```

This vignette focuses on plotting parameter estimates from MCMC draws. MCMC
diagnostic plots are covered in the separate vignette [*Visual MCMC diagnostics
using the bayesplot package*]("MCMC-diagnostics.html"), and graphical posterior
predictive checks are covered in [*Graphical posterior predictive checks using
the bayesplot package*]("PPC.rmd").

## Plots for MCMC draws

The **bayesplot** package provides various plotting functions for visualizing 
Markov chain Monte Carlo (MCMC) draws from the posterior distribution of the 
parameters of a Bayesian model.

In this vignette we'll use draws obtained using the `stan_glm` function in the 
**rstanarm** package (Gabry and Goodrich, 2016), but MCMC draws from using 
any package can be used with the functions in the **bayesplot** package. See,
for example, **brms** (which, like **rstanarm**, calls the **rstan** package
internally to use [Stan](http://mc-stan.org/)'s MCMC sampler).

```{r, eval=FALSE, results='hide'}
library("rstanarm")
fit <- stan_glm(
  mpg ~ ., # ~ . includes all other variables in dataset
  data = mtcars, 
  chains = 4, 
  iter = 2000,
  seed = 1111
)
print(fit)
```

```{r stan_glm, echo=FALSE, results='hide'}
suppressPackageStartupMessages(library("rstanarm"))
fit <- stan_glm(
  mpg ~ ., data = mtcars, 
  chains = 4, 
  iter = 2000,
  seed = 1111
)
```

```{r, print-fit, echo=FALSE}
print(fit)
```

To use the posterior draws with the functions in the **bayesplot** package
we'll extract them from the fitted model object: 

```{r, get-draws}
posterior <- as.array(fit)
dim(posterior)
dimnames(posterior)
```

### Interval estimates

Posterior intervals for the parameters can be plotted using the `mcmc_intervals`
function.

```{r, mcmc_intervals}
library("bayesplot")
color_scheme_set("red")
mcmc_intervals(posterior, pars = c("cyl", "drat", "am", "sigma"))
```

The default is to show 50% intervals (the thick lines) and 90% intervals 
(the thinner outer lines). These defaults can be changed using the `prob`
and `prob_outer` arguments, respectively. The points in the above plot 
are posterior medians. The `point_est` argument can be used to select posterior
means instead or to omit the point estimates.

To show the uncertainty intervals as shaded areas under the estimated posterior
density curves we can use the `mcmc_areas` function:

```{r, mcmc_areas}
mcmc_areas(
  posterior, 
  pars = c("cyl", "drat", "am", "sigma"),
  prob = 0.8, # 80% intervals
  prob_outer = 0.99, # 99%
  point_est = "mean"
)
```


### Histograms and density estimates

Histograms or kernel density estimates of posterior distributions of the various
model parameters can be visualized using the functions described on the 
`MCMC-distributions` page in the **bayesplot** documentation.

#### Histograms

The `mcmc_hist` and `mcmc_dens` functions plot posterior distributions (combining all chains):

```{r, mcmc_hist, message=FALSE, fig.height=3}
color_scheme_set("green")
mcmc_hist(posterior, pars = c("wt", "am"))
mcmc_dens(posterior, pars = c("wt", "am"))
```

To view the four Markov chain separately we can use  `mcmc_hist_by_chain`, `mcmc_dens_overlay`:

```{r, mcmc_hist_by_chain, message=FALSE, fig.height=3}
color_scheme_set("brightblue")
mcmc_hist_by_chain(posterior, pars = c("wt", "am"))
mcmc_dens_overlay(posterior, pars = c("wt", "am"))
```

The `mcmc_violin` function also plots the density estimates of 
each chain as violins with horizontal lines at user-specified quantiles:

```{r, mcmc_violin, fig.height=3}
mcmc_violin(posterior, pars = c("wt", "am"), probs = c(0.1, 0.5, 0.9))
```

### Scatterplots

The `mcmc_scatter` function creates a scatterplot with two parameters:

```{r, mcmc_scatter}
color_scheme_set("gray")
mcmc_scatter(posterior, pars = c("(Intercept)", "wt"), size = 1.5, alpha = 0.5)
```

In addition to `mcmc_scatter`, an `mcmc_pairs` functions for creating pairs
plots with more than two parameters will be available in one of the next
releases of the package.

### Traceplots 

The `mcmc_trace` function creates standard traceplots:

```{r, mcmc_trace}
color_scheme_set("blue")
mcmc_trace(posterior, pars = c("wt", "sigma"))
```

If it's hard to see the difference between the chains we can 
change to a mixed color scheme, for example:

```{r, change-scheme}
color_scheme_set("mix-blue-red")
mcmc_trace(posterior, pars = c("wt", "sigma"))
```

The `mcmc_trace_highlight` uses points instead of lines and reduces the opacity
of all but a single chain (which is specified using the `highlight` argument).

```{r, mcmc_trace_highlight}
mcmc_trace_highlight(posterior, pars = c("wt", "sigma"), highlight = 3)
```

## References

Gabry, J., and Goodrich, B. (2016). rstanarm: Bayesian Applied Regression 
Modeling via Stan. R package version 2.12.1. 
http://mc-stan.org/interfaces/rstanarm.html

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and
Rubin, D. B. (2013). *Bayesian Data Analysis*. Chapman & Hall/CRC Press, London,
third edition.

Stan Development Team. (2016). *Stan Modeling Language Users
Guide and Reference Manual*. http://mc-stan.org/documentation/
