nimble (version 0.8.0)

configureRJ: Configure Reversible Jump Sampler

Description

Modifies the MCMCconf object of a specific model to include a reversible jump MCMC sampler for variable selection using an univariate normal proposal distribution. Users can control the mean and scale of the proposal. This function supports two different types of model specification: with and without indicator variables.

Usage

configureRJ(mcmcConf, targetNodes, indicatorNodes = NULL, priorProb = NULL,
  control = list(mean = NULL, scale = NULL, fixedValue = NULL))

Arguments

mcmcConf

An MCMCconf object.

targetNodes

A character vector, specifying the nodes and/or variables for which variable selection is to be performed. Nodes may be specified in their indexed form, 'y[1, 3]'. Alternatively, nodes specified without indexing will be expanded fully, e.g., 'x' will be expanded to 'x[1]', 'x[2]', etc.

indicatorNodes

An optional character vector, specifying the indicator nodes and/or variables paired with targetNodes. Nodes may be specified in their indexed form, 'y[1, 3]'. Alternatively, nodes specified without indexing will be expanded fully, e.g., 'x' will be expanded to 'x[1]', 'x[2]', etc. Nodes must be provided consistently with targetNodes vector. See details.

priorProb

An optional value or vector of prior probabilities for each node to be in the model. See details.

control

An optional list of control arguments:

• mean. The mean of the normal proposal distribution (default = 0).

• scale. The standard deviation of the normal proposal distribution (default = 1).

• fixedValue. (Optional) Value for the variable when it is out of the model, which can be used only when priorProb is provided (default = 0). If specified when indicatorNodes is passed, a warning is given and fixedValue is ignored.

Value

NULL configureRJ modifies the input configuration MCMCconf in place.

Details

This function modifies the samplers in MCMCconf for each of the nodes provided in the targetNodes argument. To these elements two samplers are assigned: a specialized reversible jump sampler and a modified version of the existing sampler that is used only when the target node is already in the model. configureRJ can handle two different ways of writing a NIMBLE model, either using indicator variables or not as shown in the examples below and discussed further in the MCMC Chapter in the NIMBLE User Manual. When using indicator variables, the user should provide the indicatorNodes argument and no fixedValue argument should be used. When not using indicator variables, the user should provide the priorProb argument and no indicatorNodes argument should be provided. In the latter case, the user can provide a non-zero value for fixedValue if desired.

Note that this functionality is intended for variable selection in regression-style models but may be useful for other situations as well. At the moment, setting a variance component to zero and thereby removing a set of random effects that are explicitly part of a model will not work because MCMC sampling in that case would need to propose values for multiple parameters (the random effects), whereas this functionality only proposes to add a single parameter.

References

Peter J. Green. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4), 711-732.

See Also

sampler configureRJ

Examples

# NOT RUN {
# }
# NOT RUN {
## Linear regression with intercept and two covariates, using indicator variables

code <- nimbleCode({
  beta0 ~ dnorm(0, sd = 100)
  beta1 ~ dnorm(0, sd = 100)
  beta2 ~ dnorm(0, sd = 100)
  sigma ~ dunif(0, 100) 

  z1 ~ dbern(psi)  ## indicator variable associated with beta1
  z2 ~ dbern(psi)  ## indicator variable associated with beta2
  psi ~ dunif(0, 1) ## hyperprior on inclusion probability
  for(i in 1:N) {
    Ypred[i] <- beta0 + beta1 * z1 * x1[i] + beta2 * z2 * x2[i]
    Y[i] ~ dnorm(Ypred[i], sd = sigma)
  }
})

## simulate some data
set.seed(1)
N <- 100
x1 <- runif(N, -1, 1)
x2 <- runif(N, -1, 1) ## this covariate is not included
Y <- rnorm(N, 1 + 2.5 * x1, sd = 1)

## build the model
rIndicatorModel <- nimbleModel(code, constants = list(N = N),
                               data = list(Y = Y, x1 = x1, x2 = x2), 
                               inits=  list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y),
                               z1 = 1, z2 = 1, psi = 0.5))

indicatorModelConf <- configureMCMC(rIndicatorModel)

## Add reversible jump  
configureRJ(mcmcConf = indicatorModelConf,     ## model configuration
            targetNodes = c("beta1", "beta2"), ## coefficients for selection
            indicatorNodes = c("z1", "z2"),    ## indicators paired with coefficients
            control = list(mean = 0, scale = 2))

indicatorModelConf$addMonitors("beta1", "beta2", "z1", "z2")

rIndicatorMCMC <- buildMCMC(indicatorModelConf)
cIndicatorModel <- compileNimble(rIndicatorModel)
cIndicatorMCMC <- compileNimble(rIndicatorMCMC, project = rIndicatorModel)

set.seed(1)
samples <- runMCMC(cIndicatorMCMC, 10000, nburnin = 6000)

## posterior probability to be included in the mode
mean(samples[ , "z1"])
mean(samples[ , "z2"])

## posterior means when in the model
mean(samples[ , "beta1"][samples[ , "z1"] != 0])
mean(samples[ , "beta2"][samples[ , "z2"] != 0])


## Linear regression with intercept and two covariates, without indicator variables

code <- nimbleCode({
  beta0 ~ dnorm(0, sd = 100)
  beta1 ~ dnorm(0, sd = 100)
  beta2 ~ dnorm(0, sd = 100)
  sigma ~ dunif(0, 100)
  for(i in 1:N) {
    Ypred[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
    Y[i] ~ dnorm(Ypred[i], sd = sigma)
  }
})

rNoIndicatorModel <- nimbleModel(code, constants = list(N = N),
                                 data = list(Y = Y, x1 = x1, x2 = x2), 
                                 inits=  list(beta0 = 0, beta1 = 0, beta2 = 0, sigma = sd(Y)))

noIndicatorModelConf <- configureMCMC(rNoIndicatorModel)

## Add reversible jump  
configureRJ(mcmcConf = noIndicatorModelConf,   ## model configuration
            targetNodes = c("beta1", "beta2"), ## coefficients for selection   
            priorProb = 0.5,                   ## prior probability of inclusion
            control = list(mean = 0, scale = 2))

## add monitors
noIndicatorModelConf$addMonitors("beta1", "beta2")
rNoIndicatorMCMC <- buildMCMC(noIndicatorModelConf) 

cNoIndicatorModel <- compileNimble(rNoIndicatorModel)
cNoIndicatorMCMC <- compileNimble(rNoIndicatorMCMC, project = rNoIndicatorModel)

set.seed(1)
samples <- runMCMC(cNoIndicatorMCMC, 10000, nburnin = 6000 )

## posterior probability to be included in the mode
mean(samples[ , "beta1"] != 0)
mean(samples[ , "beta2"] != 0)

## posterior means when in the model
mean(samples[ , "beta1"][samples[ , "beta1"] != 0])
mean(samples[ , "beta2"][samples[ , "beta2"] != 0])
# }
# NOT RUN {
 
# }