MCMC for logistic regression with random effects

This example shows how to build and run MCMC for a Generalized Linear Mixed Model (GLMM), specifically a logistic regression model with random effects.

Model Creation

## load the NIMBLE library
library(nimble, warn.conflicts = FALSE)
## nimble version 0.9.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.
## define the model
code <- nimbleCode({
    beta0 ~ dnorm(0, sd = 10000)
    beta1 ~ dnorm(0, sd = 10000)
    sigma_RE ~ dunif(0, 1000)
    for (i in 1:N) {
        beta2[i] ~ dnorm(0, sd = sigma_RE)
        logit(p[i]) <- beta0 + beta1 * x[i] + beta2[i]
        r[i] ~ dbin(p[i], n[i])
    }
})

## constants, data, and initial values
constants <- list(N = 10)

data <- list(
    r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46),
    n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79),
    x = c(0,  0,  0,  0,  0,  1, 1,  1,  1,  1)
)

inits <- list(beta0 = 0, beta1 = 0, sigma_RE = 1)

## create the model object
glmmModel <- nimbleModel(code = code, constants = constants, data = data, 
                         inits = inits, check = FALSE)
## defining model...
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
## model building finished.

Default MCMC Algorithm

Now we are ready to create the default MCMC algorithm from the model object.

glmmMCMC <- buildMCMC(glmmModel)

Compile the model and MCMC algorithm

CglmmModel <- compileNimble(glmmModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
CglmmMCMC <- compileNimble(glmmMCMC, project = glmmModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.

Execute MCMC algorithm and extract samples

samples <- runMCMC(CglmmMCMC, niter = 10000)
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Customize an MCMC Algorithm

We may want to customize the MCMC with specific samplers. To do so, we first make a new instance of the model:

glmmModel <- nimbleModel(code = code, constants = constants, 
                         data = data, inits = inits)
## defining model...
## building model...
## setting data and initial values...
## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
## model building finished.

Then we make an empty MCMC configuration and add some samplers of our choice to it:

spec <- configureMCMC(glmmModel, nodes = NULL)
spec$addSampler(type = 'slice', target = 'beta0')
spec$addSampler(type = 'slice', target = 'beta1')
spec$addSampler(type = 'RW', target = 'sigma_RE')
spec$addSampler(type = 'RW_block', target = 'beta2[1:10]')
## Note: Assigning an RW_block sampler to nodes with very different scales can result in low MCMC efficiency.  If all nodes assigned to RW_block are not on a similar scale, we recommend providing an informed value for the "propCov" control list argument, or using the AFSS sampler instead.

Then we build the MCMC from the configuration:

customGlmmMCMC <- buildMCMC(spec)

Compile model and custom MCMC algorithm

CglmmModel <- compileNimble(glmmModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.
CglmmMCMC <- compileNimble(customGlmmMCMC, project = glmmModel)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## compilation finished.

Execute custom MCMC algorithm and extract samples

samples <- runMCMC(CglmmMCMC, niter = 10000)
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Note that NIMBLE does not provide its own tools for working with posterior samples.