---
title: "Inspecting Convergence of a Model"
author: "Jacob Carey, Steven Cristiano, and Robert Scharpf"
date: \today
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 lib} 
# load CNPBayes
suppressMessages(library(CNPBayes))

# load packages for manipulating and visualizing data
suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(ggplot2))
```

# Workflow

```{r post-sim}
set.seed(1)

N <- 7524
n <- 81
lrr <- replicate(N, mean(rnorm(n)))

mp <- McmcParams(iter=1000, burnin=5000, thin=1, nStarts=1)

model <- MarginalModel(data=lrr, mcmc.params=mp) 
m.list <- posteriorSimulation(model, k=1:4)
m.lik <- marginalLikelihood(m.list)
m.lik
```

```{r plot1}
data1 <- as.data.frame(theta(chains(m.list[[1]]))) %>%
    mutate(iter=1:1000) %>%
    gather(component, theta, V1) %>%
    mutate(model=1)

data2 <- as.data.frame(theta(chains(m.list[[2]]))) %>%
    mutate(iter=1:1000) %>%
    gather(component, theta, V1:V2) %>%
    mutate(model=2)

data3 <- as.data.frame(theta(chains(m.list[[3]]))) %>%
    mutate(iter=1:1000) %>%
    gather(component, theta, V1:V3) %>%
    mutate(model=3)

data4 <- as.data.frame(theta(chains(m.list[[4]]))) %>%
    mutate(iter=1:1000) %>%
    gather(component, theta, V1:V4) %>%
    mutate(model=4)

data <- bind_rows(data1, data2, data3, data4) %>%
    mutate(component=gsub("V", "", component))

ggplot(data, aes(x=iter, y=theta)) +
    geom_line(aes(colour=component, linetype=component)) +
    facet_wrap(~model, nrow=2, ncol=2) +
    theme_classic() +
    xlab("")
```

```{r plot2}
data_2.4 <- data %>%
    filter(model > 1, (iter <= 200 | iter >= 800)) %>%
    mutate(iter.cat=cut(iter, c(0, 200, 799, 1000),
                        c("bottom", "middle", "top"))) %>%
    group_by(component, iter.cat, model) %>%
    mutate(cat.model.median=median(theta)) %>%
    ungroup() %>%
    mutate(box.cat=paste(iter.cat, component))

ggplot(data_2.4, aes(x=box.cat, y=theta)) +
    geom_boxplot(aes(colour=box.cat)) +
    geom_hline(yintercept=0.0, linetype="dashed", colour="gray") +
    facet_wrap(~model, scales="free_x") +
    guides(colour=FALSE) +
    theme_classic() +
    xlab("")
```

# References
