Post-hoc MCMC with glmmTMB

Overview

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 examples 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 (e.g. with functions from the coda package) before drawing conclusions. tmbstan does some checking automatically, so it produces warnings where the DIY sampler below does not.

The first MCMC sampler illustrated below is conceptually simple (Metropolis with a multivariate normal candidate distribution). Users who want to do MCMC sampling with less code or on a regular basis should consider the tmbstan package, which does more efficient hybrid/Hamiltonian Monte Carlo sampling, and can take full advantage of glmmTMB's ability to provide gradients of the log-posterior with respect to the parameters. More info on tmbstan and the methods that are available to use on the samples can be found on Github (e.g. methods(class="stanfit")). It is safe to skip the Metropolis-Hastings (DIY) section of the salamander example below if the user prefers to use tmbstan.

In general, you can plug the objective function from glmmTMB (i.e., fitted_model$obj$fn) into any general-purpose MCMC package (e.g. adaptMCMC, MCMCpack (using MCMCmetrop1R()), ensemblemcmc, ...). However, most of these algorithms do not take advantage of glmmTMB's automatic computation of gradients; for that, we will demonstrate the use of tmbstan, which uses glmmTMB's objective and gradient functions in conjunction with the sampling algorithms from Stan.

Setup

Load packages:

library(glmmTMB)
library(coda)     ## MCMC utilities
library(reshape2) ## for melt()
## graphics
library(lattice)
library(ggplot2); theme_set(theme_bw())

Salamander example

Fit basic model:

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

##' @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)
}

## 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)

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

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

lattice::xyplot(m1,layout=c(2,3),asp="fill")

print(effectiveSize(m1),digits=3)

##    (Intercept)        minedno zi~(Intercept)     zi~minedno        sd_site 
##            262            243            258            259            249

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

Sleep study example

Now using the (now) classic sleepstudy data set included with the lme4 package:

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

Running this chain produces many warnings about divergent transitions, low effective sample size, etc.; the diagnostic table confirms this.

Diagnostics:

dp3 <- bayestestR::diagnostic_posterior(m3)

Most of the trace plot (redacted) looks OK. However, one sub-panel shows that on of the parameters (theta[3], controlling the intercept/slope correlation) is problematic for one of the chains:

One way to address these problems is to add bounds that prevent the chains from going outside a range of \(\pm 5\) standard errors from the MLE (we originally used \(\pm 3 \sigma\), but it looked like these bounds were being hit too frequently - if possible we want bounds that will prevent crazy behaviour but not otherwise constrain the chains ...)

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))

The diagnostics and the trace plot look much better now (only theta[3] shown) ...

dp4 <- bayestestR::diagnostic_posterior(m4)

The last step (for now) is to extract the sampled values, give them interpretable names, and transform the correlation parameter back to a more interpretable scale (see the "Covariance structures" vignette):

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)

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

A more principled, properly Bayesian way to handle this problem would be to add sensible regularizing priors. (To be done/tested; as of this writing, priors can be set on fixed effects and on random effect standard deviations, but not on random effect correlations ...)