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:
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:
library(glmmTMB) library(coda) ## MCMC utilities library(reshape2) ## for melt() ## graphics library(lattice) library(ggplot2); theme_set(theme_bw())
Fit basic model:
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.
This is a basic block Metropolis sampler, based on Florian Hartig's code here.
Run the chain:
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):
colnames(m1) <- c(names(fixef(fm1)[[1]]), "log(sigma)", c("log(sd_Intercept)","log(sd_Days)","cor")) m1[,"cor"] <- sapply(m1[,"cor"], get_cor)
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.
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.
ggplot(reshape2::melt(as.matrix(m1[,-1])),aes(x=Var2,y=value))+ geom_violin(fill="gray")+coord_flip()+labs(x="")
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.)
(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 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}$.)
library(png) library(grid) img <- readPNG(system.file("vignette_data","tmbstan_traceplot.png",package="glmmTMB")) grid.raster(img)
Owls
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.