Using derived quantities in NIMBLE’s MCMC system

NIMBLE is a system for building and sharing analysis methods for statistical models, especially for hierarchical models and computationally-intensive methods (such as MCMC, Laplace approximation, and sequential Monte Carlo).

As of version 1.4.0, NIMBLE allows one to specify “derived quantities” to calculate or record additional quantities of interest during an MCMC. Potential derived quantities include:

• posterior means,
• posterior variances,
• log-probabilities (log-density values) for single or (summed) groups of nodes,
• posterior predictive samples for predictive nodes, and
• arbitrary calculations by providing a user-defined nimbleFunction that follows the required specification.

By using derived quantities, users can avoid adding additional nodes to the model itself and can potentially reduce how much is stored during an MCMC run and saved as output from the run, thereby reducing memory and disk use.

We’ll give a brief example of saving log-probability (log-posterior-density) values summed over a set of nodes in a simple Poisson random effects model.

library(nimble)

pumpCode <- nimbleCode({ 
    for(i in 1:N) {
        theta[i] ~ dgamma(alpha, beta)
        lambda[i] <- theta[i]*t[i]
        x[i] ~ dpois(lambda[i])
    }
    alpha ~ dexp(1.0)
    beta ~ dgamma(0.1, 1.0)
})

pumpConsts <- list(N = 10,
                   t = c(94.3, 15.7, 62.9, 126, 5.24,
                         31.4, 1.05, 1.05, 2.1, 10.5))
pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
pumpInits <- list(alpha = 1, beta = 1,
                  theta = rep(0.1, pumpConsts$N))
pumpModel <- nimbleModel(pumpCode, pumpConsts, pumpData, pumpInits)

Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Checking model sizes and dimensions

Now we add the logProb derived quantity, specified directly as an argument to configureMCMC.

conf <- configureMCMC(pumpModel, logProb = c('alpha', 'beta'))

===== Monitors =====
thin = 1: alpha, beta
===== Samplers =====
RW sampler (1)
  - alpha
conjugate sampler (11)
  - beta
  - theta[]  (10 elements)
===== DerivedQ =====
- logProb (1)

conf$printDerivedQuantities()

[1] derived quantity: logProb,  execution interval: thin,  nodes: c(alpha, beta)

When we run the MCMC, we see the derived quantities are saved as an element of the output list, in addition to the usual posterior samples of the parameters.

mcmc <- buildMCMC(conf)

cPumpModel <- compileNimble(pumpModel)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

cmcmc <- compileNimble(mcmc)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

runMCMC(cmcmc, niter = 100, thin = 20)

running chain 1...

|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

$samples
         alpha      beta
[1,] 0.7434957 0.8437510
[2,] 0.5046057 0.8152505
[3,] 0.5201513 1.3765454
[4,] 1.1023278 2.2564913
[5,] 1.4316249 2.1598157

$derived
$derived$logProb
          alpha      beta
[1,] -0.7434957 -2.943556
[2,] -0.5046057 -2.884129
[3,] -0.5201513 -3.916877
[4,] -1.1023278 -5.241634
[5,] -1.4316249 -5.105549

We could also have added the logProb derived quantity to an existing MCMC configuration argument (and also illustrating that one can modify the interval at which recording is done):

conf <- configureMCMC(pumpModel)
conf$addDerivedQuantity('logProb', nodes = c('alpha', 'beta'), interval = 10)

The operation of every derived quantity function is governed by an interval parameter. Each derived quantity function is executed at the end of every interval MCMC iterations (and only for iterations after the burnin period). The default value for interval is given by the thin interval of the MCMC. For example, if interval = 1 (regardless of the thin value), then that derived quantity function will execute at the end of every MCMC sampling iteration.

Note that user-defined derived quantity nimbleFunctions provide a general way for a user to intervene at the end of an MCMC iteration and do something. This would often be recording a value (to be returned at the end of the MCMC) or updating some summary statistic, but the derived quantities system doesn’t require this, giving users flexibility.

More details can be found in the NIMBLE User Manual and via help(derived).