# 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.13.0 is loaded.
##
## 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 AFSS 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.