inst/doc/mcmc.R

In glmmTMB: Generalized Linear Mixed Models using Template Model Builder

## ----knitr_setup, include=FALSE, message=FALSE--------------------------------
library(knitr)
opts_chunk$set(echo = TRUE)
## OK to evaluate on CRAN since we have stored all the slow stuff ...
##               eval = identical(Sys.getenv("NOT_CRAN"), "true"))
knitr::read_chunk(system.file("vignette_data", "mcmc.R", package="glmmTMB"))
L <- load(system.file("vignette_data", "mcmc.rda", package="glmmTMB"))
library(png)
library(grid)

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

## ----fit1---------------------------------------------------------------------
fm1 <- glmmTMB(count ~ mined + (1|site),
    zi=~mined,
    family=poisson, data=Salamanders)

## ----run_MCMC-----------------------------------------------------------------
##' @param start starting value
##' @param V variance-covariance matrix of MVN candidate distribution
##' @param iterations total iterations
##' @param nsamp number of samples to store
##' @param burnin number of initial samples to discard
##' @param thin thinning interval
##' @param tune tuning parameters; expand/contract V
##' @param seed random-number seed
run_MCMC <- function(start,
                     V,   
                     logpost_fun,
                     iterations = 10000,
                     nsamp = 1000,
                     burnin = iterations/2,
                     thin = floor((iterations-burnin)/nsamp),
                     tune = NULL,
                     seed=NULL
                     ) {
    ## initialize
    if (!is.null(seed)) set.seed(seed)
    if (!is.null(tune)) {
        tunesq <- if (length(tune)==1) tune^2 else outer(tune,tune)
        V <-  V*tunesq
    }
    chain <- matrix(NA, nsamp+1, length(start))
    chain[1,] <- cur_par <- start
    postval <- logpost_fun(cur_par)
    j <- 1
    for (i in 1:iterations) {
        proposal = MASS::mvrnorm(1,mu=cur_par, Sigma=V)
        newpostval <- logpost_fun(proposal)
        accept_prob <- exp(newpostval - postval)
        if (runif(1) < accept_prob) {
            cur_par <- proposal
            postval <- newpostval
        }
        if ((i>burnin) && (i %% thin == 1)) {
            chain[j+1,] <- cur_par
            j <- j + 1
        }
    }
    chain <- na.omit(chain)
    colnames(chain) <- names(start)
    chain <- coda::mcmc(chain)
    return(chain)
}

## ----setup--------------------------------------------------------------------
## FIXME: is there a better way for user to extract full coefs?
rawcoef <- with(fm1$obj$env,last.par[-random])
names(rawcoef) <- make.names(names(rawcoef), unique=TRUE)
## log-likelihood function 
## (run_MCMC wants *positive* log-lik)
logpost_fun <- function(x) -fm1$obj$fn(x)
## check definitions
stopifnot(all.equal(c(logpost_fun(rawcoef)),
                    c(logLik(fm1)),
          tolerance=1e-7))
V <- vcov(fm1,full=TRUE)

## ----do_run_MCMC,eval=FALSE---------------------------------------------------
# t1 <- system.time(m1 <- run_MCMC(start=rawcoef,
#                                  V=V, logpost_fun=logpost_fun,
#                                  seed=1001))

## ----add_names----------------------------------------------------------------
colnames(m1) <- colnames(vcov(fm1, full = TRUE))
colnames(m1)[ncol(m1)] <- "sd_site"

## ----traceplot,fig.width=10, fig.height = 7, eval = FALSE---------------------
# lattice::xyplot(m1,layout=c(2,3),asp="fill")

## ----effsize------------------------------------------------------------------
print(effectiveSize(m1),digits=3)

## ----violins,echo=FALSE, fig.width = 6, fig.height = 6------------------------
m_long <- reshape2::melt(as.matrix(m1[,-1]))
ggplot(m_long, aes(x=Var2, y=value))+
    geom_violin(fill="gray")+
    coord_flip()+labs(x="")

## ----do_tmbstan,eval=FALSE----------------------------------------------------
# library(tmbstan)
# t2 <- system.time(m2 <- tmbstan(fm1$obj, seed = 101))

## ----diagnostic_tab, echo = FALSE---------------------------------------------
knitr::kable(dp2, digits = c(0, 0, 3, 3))

## ----show_traceplot,echo=FALSE,fig.width=10,fig.height=5,eval = FALSE---------
# img <- readPNG(system.file("vignette_data","tmbstan_traceplot.png",package="glmmTMB"))
# grid.raster(img)

## ----show_pairsplot,echo=FALSE,fig.width=8,fig.height=8, eval=FALSE-----------
# img <- readPNG(system.file("vignette_data","tmbstan_pairsplot.png",package="glmmTMB"))
# grid.raster(img)

## ----sleepstudy_tmbstan, eval = FALSE-----------------------------------------
# data("sleepstudy", package = "lme4")
# fm2 <- glmmTMB(Reaction ~ Days + (Days | Subject), data = sleepstudy)
# t3 <- system.time(m3 <- tmbstan(fm2$obj, seed = 101))

## ----sleepstudy_diag, eval = FALSE--------------------------------------------
# dp3 <- bayestestR::diagnostic_posterior(m3)

## ----sleepstudy_diag_tab, echo = FALSE----------------------------------------
knitr::kable(dp3, digits = c(0, 0, 3, 3))

## ----sleepstudy_trace,fig.width=10,fig.height=5, echo = FALSE, eval = FALSE----
# img <- readPNG(system.file("vignette_data","sleepstudy_traceplot.png",package="glmmTMB"))
# grid.raster(img)

## ----sleepstudy_trace_theta3,fig.width=5,fig.height=5, echo = FALSE-----------
img <- readPNG(system.file("vignette_data","sleepstudy_traceplot_theta3.png",package="glmmTMB"))
grid.raster(img)

## ----sleepstudy_tmbstan_bounds, eval = FALSE----------------------------------
# sdrsum <- summary(fm2$sdr)
# par_est <- sdrsum[,"Estimate"]
# par_sd <- sdrsum[,"Std. Error"]
# t4 <- system.time(m4 <- tmbstan(fm2$obj,
#                                 lower = par_est - 5*par_sd,
#                                 upper = par_est + 5*par_sd,
#                                 seed = 101))

## ----sleepstudy_bounds_diag, eval = FALSE-------------------------------------
# dp4 <- bayestestR::diagnostic_posterior(m4)

## ----sleepstudy_bounds_diag_tab, echo = FALSE---------------------------------
knitr::kable(dp4, digits = c(0, 0, 3, 3))

## ----sleepstudy_trace_bounds,fig.width=10,fig.height=5, echo = FALSE, eval = FALSE----
# img <- readPNG(system.file("vignette_data","sleepstudy_traceplot_bounds.png",package="glmmTMB"))
# grid.raster(img)

## ----sleepstudy_trace_bounds_theta3,fig.width = 5, fig.height=5, echo = FALSE----
img <- readPNG(system.file("vignette_data","sleepstudy_traceplot_bounds_theta3.png",package="glmmTMB"))
grid.raster(img)

## ----trans_param, eval = FALSE------------------------------------------------
# samples4 <- as.data.frame(extract(m4, pars=c("beta","betadisp","theta")))
# colnames(samples4) <- c(names(fixef(fm2)$cond),
#                   "log(sigma)",
#                   c("log(sd_Intercept)", "log(sd_Days)", "cor"))
# samples4$cor <- sapply(samples4$cor, get_cor)

## ----sleepstudy_hist, fig.width = 10, fig.height = 5--------------------------
m4_long <- reshape2::melt(as.matrix(samples4))
ggplot(m4_long, aes(x = value)) + geom_histogram(bins = 50) + facet_wrap(~Var2, scale = "free")