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)

## 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
Rmodel <- nimbleModel(code=code, constants=constants, data=data, inits=inits, check = FALSE)
## defining model...
## building model...
## setting data and initial values...
## model building finished

Default MCMC Algorithm

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

Rmcmc <- buildMCMC(Rmodel)

Compile the model and MCMC algorithm

Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

Execute MCMC algorithm and extract samples

Cmcmc$run(10000)
## NULL
samples <- as.matrix(Cmcmc$mvSamples)

Customize an MCMC Algorithm

First we make a new instance of the model:

Rmodel <- nimbleModel(code=code, constants=constants, data=data, inits=inits)
## defining model...
## building model...
## setting data and initial values...
## checking model...   (use nimbleModel(..., check = FALSE) to skip model check)
## NAs were detected in model variables: beta2, logProb_beta2, p, logProb_r.
## model building finished

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

spec <- configureMCMC(Rmodel, nodes=NULL)
spec$addSampler(type = 'slice', target ='beta0')
## [1] slice sampler: beta0,  adaptive: TRUE,  adaptInterval: 200,  sliceWidth: 1,  sliceMaxSteps: 100
spec$addSampler(type = 'slice', target ='beta1')
## [2] slice sampler: beta1,  adaptive: TRUE,  adaptInterval: 200,  sliceWidth: 1,  sliceMaxSteps: 100
spec$addSampler(type = 'RW', target ='sigma_RE')
## [3] RW sampler: sigma_RE,  adaptive: TRUE,  adaptInterval: 200,  scale: 1
spec$addSampler(type = 'RW_block', target ='beta2[1:10]')
## [4] RW_block sampler: beta2[1:10],  adaptive: TRUE,  adaptScaleOnly: FALSE,  adaptInterval: 200,  scale: 1,  propCov: identity

Then we build the MCMC from the configuration:

Rmcmc <- buildMCMC(spec)

Compile model and custom MCMC algorithm

Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

Execute custom MCMC algorithm and extract samples

Cmcmc$run(10000)
## NULL
samples <- as.matrix(Cmcmc$mvSamples)