---
title: "Inspecting Convergence of a Model"
author: "Jacob Carey, Steven Cristiano, and Robert Scharpf"
date: "`r doc_date()`"
output: BiocStyle::pdf_document
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Overview of CNPBayes package}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc} 
---

# Introduction

A Markov Chain Monte Carlo posterior simulation should be visually inspected to assess convergence.

```{r packages, message=FALSE}
# load CNPBayes
suppressMessages(library(CNPBayes))
# load packages for manipulating and visualizing data
suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(ggplot2))
```

# Workflow

The number of starting values, burnin MCMC simulatioms, and MCMC simulations
after burnin are controlled by a `McmcParams` object. Here, we specify a small
number so that the example runs quickly. As discussed in [overview](Overview.Rmd), four types of models are possible--SB, MB, SBP, and MBP.  Below we fit only the SBP model with 2 components.

```{r post-sim}
set.seed(1)
N <- 200
n <- 81
lrr <- c(rnorm(100, -0.5, sd=0.1), rnorm(100, 0, sd=0.1))
mp <- McmcParams(iter=600, burnin=10, nStarts=4)
mlist <- gibbs(model="SBP", mp=mp, dat=lrr, k_range=c(2, 2))
```

Had more than 1 model been fit, the models in `mlist` would be sorted by decreasing value of the marginal likelihood.  Note, the marginal likelihood is only meaningful if the models have converged.  By default, the function `gibbs` requires that the effective sample size (number of independent MCMC draws) for all parameter chains is at least 500 and the multivariate Gelman Rubin diagnostic to be less than 1.2.  Otherwise, starting values for 4 (`nStarts`) new models will be independently initialized and the burnin and number of `thin` iterations will be increased by a factor of 2.  This process is repeated until the maximum number of burnin iterations has been reached (see `max_burnin`) or the effective size and Gelman Rubin criteria are both satisfied.  For each model returned, the `nStarts` chains will be combined.  In addition to the above automated diagnostics, we suggest plotting the chains as illustrated below.

```{r chains}
model <- mlist[[1]]
figs <- ggChains(model)
figs[[1]]
figs[[2]]
```

The marginal likelihood of the models can be retrieved using the `marginal_lik` accessor.

```{r marginal-lik, eval=FALSE}
m.lik <- round(sapply(mlist, marginal_lik), 1)
m.lik
```

# References