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 1.3.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
## 
## Note for advanced users who have written their own MCMC samplers:
##   As of version 0.13.0, NIMBLE's protocol for handling posterior
##   predictive nodes has changed in a way that could affect user-defined
##   samplers in some situations. Please see Section 15.5.1 of the User Manual.

## 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
##   [Note] Any error reports that follow may simply reflect missing values in model variables.

## Checking model sizes and dimensions

##   [Note] 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).

Default MCMC Algorithm

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

glmmMCMC <- buildMCMC(glmmModel)

## ===== Monitors =====
## thin = 1: beta0, beta1, sigma_RE
## ===== Samplers =====
## RW sampler (13)
##   - beta0
##   - beta1
##   - sigma_RE
##   - beta2[]  (10 elements)

Compile the model and MCMC algorithm

CglmmModel <- compileNimble(glmmModel)

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

CglmmMCMC <- compileNimble(glmmMCMC, project = glmmModel)

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

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
##   [Note] Any error reports that follow may simply reflect missing values in model variables.

## Checking model sizes and dimensions

##   [Note] 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).

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

spec <- configureMCMC(glmmModel, nodes = NULL)

## ===== Monitors =====
## thin = 1: beta0, beta1, sigma_RE
## ===== Samplers =====
## (no samplers assigned)
## ===== Comments =====

##   [Warning] No samplers assigned for 13 nodes, use conf$getUnsampledNodes() for node names.

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 "barker" sampler instead.

Then we build the MCMC from the configuration:

customGlmmMCMC <- buildMCMC(spec)

Compile model and custom MCMC algorithm

CglmmModel <- compileNimble(glmmModel)

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

CglmmMCMC <- compileNimble(customGlmmMCMC, project = glmmModel)

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

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.