Reversible Jump PDMP Samplers

Provides an implementation of the reversible jump piecewise deterministic Markov processes (PDMPs) methods developed in the paper Reversible Jump PDMP Samplers for Variable Selection (Chevallier, Fearnhead, Sutton 2020, ). It also contains an implementation of a Gibbs sampler for variable selection in Logistic regression based on Polya-Gamma augmentation.

Reversible Jump PDMP methods R package

rjpdmp is an R package that provides an implementation of the reversible jump PDMP methods developed in the paper Reversible Jump PDMP Samplers for Variable Selection (Chevallier, Fearnhead, Sutton 2020, https://arxiv.org/abs/2010.11771). It also contains an implementation of a Gibbs sampler for variable selection in Logistic regression based on Polya-Gamma augmentation (code for simulating PG variables adapted from https://github.com/jgscott/helloPG).

Installation

The preferred install is from CRAN:

install.packages("rjpdmp")

Using devtools install with:

# library(devtools)
# install_github("matt-sutton/rjpdmp")

Example Usage (Logistic regression)

generate.logistic.data <- function(beta, n.obs, Sig) {
p <- length(beta)
dataX <- MASS::mvrnorm(n=n.obs,mu=rep(0,p),Sigma=Sig)
vals <- dataX %*% as.vector(beta)
generateY <- function(p) { rbinom(1, 1, p)}
dataY <- sapply(1/(1 + exp(-vals)), generateY)
return(list(dataX = dataX, dataY = dataY))
}

n <- 15
p <- 25
beta <- c(1, rep(0, p-1))
Siginv <- diag(1,p,p)
Siginv[1,2] <- Siginv[2,1] <- 0.9
set.seed(1)
data <- generate.logistic.data(beta, n, solve(Siginv))
ppi <- 2/p

zigzag_fit <- zigzag_logit(maxTime = 5, dataX = data$dataX, datay = data$dataY,
                           prior_sigma2 = 10,theta0 = rep(0, p), x0 = rep(0, p), 
                           rj_val = 0.6, ppi = ppi)

gibbs_fit <- gibbs_logit(maxTime = 5, dataX = data$dataX, datay = data$dataY,
                         prior_sigma2 = 10,beta = rep(0,p), gamma = rep(0,p),
                         ppi = ppi)

# Plotting
plot_pdmp(zigzag_fit, coords = 1:2, inds = 1:10^4,burn = .1, nsamples = 5*1e4, 
          mcmc_samples = t(gibbs_fit$beta*gibbs_fit$gamma))

# Check how frequently models were visited
models_seen <- models_visited(zigzag_fit$theta)

# Look at the top ten visited models
models_seen[1:10,]

# Calculate model probabilities:
analysis <- model_probabilities(zigzag_fit$times, zigzag_fit$theta,
                    models = models_seen[1:10,1:p], marginals=1:p)

# Probability of specific models 
analysis$prob_mod

# Marginal probability of inclusion
analysis$marginal_prob