---
title: 'Package: Bayesplot'
author: "SongZhai"
date: ""
output:
  pdf_document: default
  html_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
```

### Description of 'bayesplot'

- Plotting functions for posterior analysis, model checking, and MCMC diagnostics. 

### Packages I will use for this tutorial

```{r, message=FALSE, warning=FALSE}
library("bayesplot")
library("ggplot2")
# install.packages("retanarm")
library("rstanarm") 
library("mcmcse")
```

### \textcolor{red}{0. Outline}

1. Simulate the data
- stan_glm
2. Visualize statistic information of the chain
- MCMC-intervals: mean, median, quantiles, $95\%$ confidence interval
- MCMC-recover: compare MCMC estimate with ``true''
- MCMC-distribution: histogram, density estimation
3. MCMC-diagnostics
- mcmc_trace: trace
- mcmc_acf: autocorrelation
- mcmc_rhat: $\hat R$ statistic
- mcmc_neff: effective sample size

### \textcolor{red}{1. Simulate the data}

\hrulefill

### stan_glm: Solve Bayesian generalized linear models via MCMC

\hrulefill

```{r}
data(mtcars)
mtcars <- mtcars[,1:4]
head(mtcars)
```

```{r}
fit <- stan_glm(mpg ~ .,data = mtcars,seed = 1)
```
```{r}
# Posterior MCMC draws with 4 chains
posterior <- as.array(fit)

# Take 1st chain as example
chain <- posterior[,1,]

# Gibbs Sampler will give us a 1000*5 matrix
# 1000 = # of iterations
# 5 = # of parameters
head(chain)
#
summary(fit)
# Parameters we are interested in
mypara <- chain[,c("cyl","sigma")]
#
head(mypara)
```
### As a conclusion
- fit: Bayesian generalized linear regression by MCMC along with Stan method
- posterior: 3-D matrix, [# of iterations, # of chains, # of parameteres]
- chain: 2-D matrix / 1st chain of posterior, [# of iterations, # of parameters]
- mypara: 2-D matrix / subset of chain, [# of iterations, # of parameters we are interested in]

### \textcolor{red}{2. Visualize statistic information of the chain}

\hrulefill

### MCMC-intervals: Plot interval estimates from MCMC draws.

\hrulefill

### Description
- mcmc_intervals: Plot central (quantile-based) posterior interval estimates from MCMC draws.
- mcmc_areas: To show the uncertainty intervals as shaded areas under the estimated posterior density curves.

### Usage
- mcmc_intervals(x,pars = character(),prob = 0.5,prob_outer = 0.9,point_est = c("median", "mean", "none"))

### Arguments
- x: A 2-D matrix / data frame of MCMC draws.
- par: Parameters you are interested in.
- prob: The probability mass to include in the inner interval (for mcmc_intervals) or in the shaded region (for mcmc_areas). The default is 0.5 (50% interval).
- prob_outer: The probability mass to include in the outer interval. The default is 0.9 for mcmc_intervals (90% interval) and 1 for mcmc_areas.
- point_est: The point estimate to show. Either "median" (the default), "mean", or "none".

```{r, message=FALSE, warning=FALSE}
color_scheme_set("brightblue")
#
mcmc_intervals(mypara, pars = c("cyl","sigma"))
#
mcmc_areas(
  mypara, 
  pars = c("cyl","sigma"),
  prob = 0.5,
  prob_outer = 0.9, 
  point_est = "mean"
)
```

\hrulefill

### MCMC-recover: Compare MCMC estimates to "true" parameter values

\hrulefill

### Description
- Plots comparing MCMC estimates to "true" parameter values. Before fitting a model to real data it is useful to simulate data according to the model using known (fixed) parameter values and to check that these "true" parameter values are (approximately) recovered by fitting the model to the
simulated data.

### Usage
- mcmc_recover_intervals(x,true,prob = 0.5,prob_outer = 0.9,point_est = c("median", "mean", "none"),size = 4, alpha = 1)
- mcmc_recover_scatter(x,true,point_est = c("median", "mean"), size = 3,alpha = 1)

### Arguments
- x: A 2-D matrix / data frame of MCMC draws.
- true: A numeric vector of "true" values of the parameters in x. There should be one value in true for each parameter included in x and the order of the parameters in true should be the same as the order of the parameters in x.
- prob: The probability mass to include in the inner interval for mcmc_recover_intervals. The default is 0.5 (50% interval).
- prob_outer: The probability mass to include in the outer interval. The default is 0.9 for mcmc_recover_intervals (90% interval).
- point_est: The point estimate to show. Either "median" (the default), "mean", or "none".
- size, alpha: Passed to 'geom_point' to control the appearance of plotted points.


```{r}
true <- apply(mypara,2,mean)
true
#
color_scheme_set("red")
mcmc_recover_intervals(mypara,true,
                       prob_outer = 0.95,
                       size = 4,alpha = 0.3)
mcmc_recover_scatter(mypara,true,
                       size = 4,alpha = 0.3)
```

\hrulefill

### MCMC-distributions: Histograms and kernel density plots of MCMC draws

\hrulefill

### Description
- Various types of histograms and kernel density plots of MCMC draws.

```{r, message=FALSE, warning=FALSE}
mcmc_hist(mypara)
mcmc_dens(mypara,pars = c("cyl","sigma"))
mcmc_combo(mypara,c("dens","hist"))
```

### \textcolor{red}{3. Whether or not my chain is good enough to use}

\hrulefill

### MCMC-diagnostics: General MCMC diagnostics

\hrulefill

\hrulefill

### mcmc_trace: Trace plot (time series plot) of MCMC draws

\hrulefill
```{r}
color_scheme_set("pink")

mcmc_trace(mypara,pars = c("cyl","sigma"),
           facet_args = list(ncol=1,strip.position="left"))
# compare results with function 'plot.ts' from package 'mcmcse'
plot.ts(mypara)
```

\hrulefill

### mcmc_acf: Plots of autocorrelation of MCMC draws

\hrulefill

```{r}
# bar plot
mcmc_acf_bar(mypara,lags = 30)

# line plot
mcmc_acf(mypara,pars = c("cyl","sigma"),
         # facet_args = list(ncol=1,strip.position="left"),
         lags = 30)
# compare results with function 'acf' from package 'mcmcse' 
par(mfrow=c(1,2))
acf(mypara[,"cyl"])
acf(mypara[,"sigma"])
par(mfrow=c(1,1))
```

\hrulefill

### mcmc_rhat: Plots of Rhat statistics of MCMC draws

### mcmc_neff: Plots of ratios of effective sample size to total sample size ofMCMC draws

\hrulefill

### Rhat: potential scale reduction statistic

One way to monitor whether a chain has converged to the equilibrium distribution is to compare its behavior to other randomly initialized chains. This is the motivation for the Gelman and Rubin (1992) potential scale reduction statistic, $\hat R$. The $\hat R$ statistic measures the ratio of the average variance of samples within each chain to the variance of the pooled samples across chains; if all chains are at equilibrium, these will be the same and $\hat R$ will be one. If the chains have not converged to a common distribution, the $\hat R$ statistic will be greater than one.

```{r}
# Go back to 'fit'
# fit <- stan_glm(mpg ~ .,data = mtcars,seed = 1)
#
# rhat
rhat <- rhat(fit)[c(2,5)]
rhat
# ratio of effective sample size to the total sample size by 'neff_ratio'
ratio <- neff_ratio(fit)[c(2,5)]
ratio
# # ratio of effective sample size to the total sample size by 'ess'
ess(mypara)/1000
```
```{r}
color_scheme_set("brightblue")
mcmc_rhat(rhat) + yaxis_text(hjust = 1)
```

```{r}
mcmc_neff(ratio) + yaxis_text(hjust = 1)
```