---
title: "Post-hoc MCMC with glmmTMB"
author: "Ben Bolker"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Post-hoc MCMC with glmmTMB}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

One commonly requested feature is to be able to run a *post hoc* Markov chain Monte Carlo analysis based on the results of a frequentist fit. This is often a reasonable shortcut for computing confidence intervals and p-values that allow for finite-sized samples rather than relying on asymptotic sampling distributions. This vignette shows an example of such an analysis. Some caveats:

- when using such a "pseudo-Bayesian" approach, be aware that using a scaled likelihood (implicit, improper priors) can often cause problems, especially when the model is poorly constrained by the data
- in particular, models with poorly constrained random effects (singular or nearly singular) are likely to give bad results
- as shown below, even models that are well-behaved for frequentist fitting may need stronger priors to give well-behaved MCMC results
- as with all MCMC analysis, it is the *user's responsibility to check for proper mixing and convergence of the chains* before drawing conclusions
- the first MCMC sampler illustrated below is simple (Metropolis with a multivariate normal candidate distribution). Users who want to do MCMC sampling on a regular basis should consider the [tmbstan package](https://CRAN.R-project.org/package=tmbstan), which does much more efficient hybrid/Hamiltonian Monte Carlo sampling.

```{r knitr_setup, include=FALSE, message=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE,
               eval = identical(Sys.getenv("NOT_CRAN"), "true"))
rc <- knitr::read_chunk
rc(system.file("vignette_data","mcmc.R",package="glmmTMB"))
```

Load packages:

```{r libs,message=FALSE}
library(glmmTMB)
library(coda)     ## MCMC utilities
library(reshape2) ## for melt()
## graphics
library(lattice)
library(ggplot2); theme_set(theme_bw())
```

Fit basic model:
```{r fit1}
```

Set up for MCMC: define scaled log-posterior function (in this case the log-likelihood function); extract coefficients and variance-covariance matrices as starting points.
```{r setup}
```

This is a basic block Metropolis sampler, based on Florian Hartig's code [here](https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/).

```{r run_MCMC}
```

Run the chain:

```{r do_run_MCMC,eval=FALSE}
```

```{r load_MCMC, echo=FALSE, eval=TRUE}
L <- load(system.file("vignette_data", "mcmc.rda", package="glmmTMB"))
```
(running this chain takes `r round(t1["elapsed"],1)` seconds)

Add more informative names and transform correlation parameter (see vignette on covariance structures and parameters):

```{r add_names}
colnames(m1) <- c(names(fixef(fm1)[[1]]),
                  "log(sigma)",
                  c("log(sd_Intercept)","log(sd_Days)","cor"))
m1[,"cor"] <- sapply(m1[,"cor"], get_cor)
```

```{r traceplot,fig.width=7}
xyplot(m1,layout=c(2,3),asp="fill")
```

The trace plots are poor, especially for the correlation; the effective sample size backs this up, as would any other diagnostics we did.

```{r effsize}
print(effectiveSize(m1),digits=3)
```

**In a real analysis we would stop and fix the mixing/convergence problems before proceeding**; for this simple sampler, some of our choices would be (1) simply run the chain for longer; (2) tune the candidate distribution (e.g. by using `tune` to scale some parameters, or perhaps by switching to a multivariate Student t distribution [see the `mvtnorm` package]); (3) add regularizing priors.

Ignoring the problems and proceeding, we can compute column-wise quantiles or highest posterior density intervals (`coda::HPDinterval`) to get confidence intervals.
Plotting posterior distributions, omitting
the intercept because it's on a very different scale.

```{r violins,echo=FALSE}
ggplot(reshape2::melt(as.matrix(m1[,-1])),aes(x=Var2,y=value))+
         geom_violin(fill="gray")+coord_flip()+labs(x="")
```

## tmbstan

The `tmbstan` package allows direct, simple access to a hybrid/Hamiltonian Monte Carlo algorithm for sampling from a TMB object; the `$obj` component of a `glmmTMB` fit is such an object. (To run this example you'll need to install the `tmbstan` package and its dependencies.)

```{r do_tmbstan,eval=FALSE}
```
(running this command, which creates 4 chains, takes `r round(t2["elapsed"],1)` seconds)

However, there are many indications (warning messages; trace plots) that the correlation parameter needs to be given a more informative prior. (In the plot below, the correlation parameter is shown on its unconstrained scale; the actual correlation would be $\theta_3/\sqrt{1+\theta_3^2}$.)

```{r show_traceplot,echo=FALSE,fig.width=8,fig.height=5}
library(png)
library(grid)
img <- readPNG(system.file("vignette_data","tmbstan_traceplot.png",package="glmmTMB"))
grid.raster(img)
```

## To do

- solve mixing for cor parameter
- more complex example - e.g. `Owls`