Example of NIMBLE’s MCMC Suite

Daniel Turek

20 August 2015

Back to NIMBLE Vignettes

 

The MCMCsuite() function in NIMBLE can be used to compare MCMC algorithms of NIMBLE, WinBUGS, OpenBUGS, and JAGS, in addition to custom-defined NIMBLE MCMC algorithms. In this example I’ll only include NIMBLE’s default MCMC and JAGS, since I’m running on OSX.

load NIMBLE library

library(nimble)

Define model, constants, data, and initial values

Here we present a simple logistic random effects model. Any BUGS- or JAGS-compatible model would work. In fact, NIMBLE extends the allowable model syntax beyond that of BUGS and JAGS substantially.

code <- nimbleCode({
    beta0 ~ dnorm(0, 0.001)
    beta1 ~ dnorm(0, 0.001)
    tau_RE ~ dgamma(1, 0.001)
    for(i in 1:N) {
        beta2[i] ~ dnorm(0, tau_RE)
        logit(p[i]) <- beta0 + beta1 * x[i] + beta2[i]
        r[i] ~ dbin(p[i], n[i])
    }
})

constants <- list(
    N = 10,
    x = c(0,  0,  0,  0,  0,  1, 1,  1,  1,  1),
    n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79)
)

data <- list(
    r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46)
)

inits <- list(beta0 = 0, beta1 = 0, tau_RE = 1, beta2 = rep(0, 10))

NIMBLE’s MCMC Suite

Use ?MCMCsuite for more details of NIMBLE’s MCMC Suite options and functionality.

For now, we’ll just run NIMBLE’s default and JAGS MCMCs.

out <- MCMCsuite(
    code = code,
    constants = constants,
    data = data,
    inits = inits,
    MCMCs = c('nimble', 'jags'),   ## could also include 'WinBUGS' or 'OpenBUGS'
    niter = 10000,
    burnin = 2000,
    makePlot = FALSE   ## optionally creates posterior density and trace plots for each parameter
)

Summary statistics and timing

By default, MCMC Suite calculates the following summary statistics. They can be changed, or custom statistics can be added using the summaryStats argument to MCMCsuite().

out$summary

## , , beta0
##  
##              mean     median        sd   CI95_low   CI95_upp
## nimble -0.5626499 -0.5591506 0.1320016 -0.8356139 -0.3111702
## jags   -0.5656397 -0.5630073 0.1431934 -0.8450105 -0.3015801
##  
## , , beta1
##  
##            mean   median        sd  CI95_low CI95_upp
## nimble 1.311850 1.310690 0.1926143 0.9429042 1.696294
## jags   1.311689 1.305311 0.1937822 0.9426383 1.703211
##  
## , , tau_RE
##  
##            mean   median       sd CI95_low CI95_upp
## nimble 805.1271 489.2565 922.1882 7.659215 3271.556
## jags   819.6766 517.6022 928.4239 9.218182 3367.238

Runtime of each algorithm, in seconds, and also the time required for NIMBLE compilation to C++.

out$timing

##        nimble           jags nimble_compile 
##         0.154          0.788          8.369 

All samples are also avaialble in out$samples, which is a 3-D array indexed by (1) MCMC algorithm, (2) model parameter, and (3) MCMC iteration. So in this case, a 2 x 3 x 8000 array.

dim(out$samples)

## [1]    2    3 8000

dimnames(out$samples)

## [[1]]
## [1] "nimble" "jags"  
##  
## [[2]]
## [1] "beta0"  "beta1"  "tau_RE"
##  
## [[3]]
## NULL

apply(out$samples, c(1,2), mean)

##             beta0    beta1   tau_RE
## nimble -0.5626499 1.311850 805.1271
## jags   -0.5656397 1.311689 819.6766

Closing

Yes, you can easily compare MCMC algorithms using NIMBLE. And do a whole lot more, too!

I encourage you to check-out our User Manual, or join the NIMBLE User’s Google Group (nimble-users) for help or further questions.

Cheers!

Daniel