poisson.spike: Spike and slab Poisson regression

Description

MCMC algorithm for Poisson regression models with a 'spike-and-slab' prior that places some amount of posterior probability at zero for a subset of the coefficients.

Usage

poisson.spike(formula,
              exposure = 1,
              niter,
              data,
              subset,
              prior = NULL,
              na.action = options("na.action"),
              contrasts = NULL,
              drop.unused.levels = TRUE,
              initial.value = NULL,
              ping = niter / 10,
              nthreads = 4,
              seed = NULL,
              ...)

Value

Returns an object of class poisson.spike. The returned object is a list with the following elements.

beta: A niter by ncol(x) matrix of regression coefficients, many of which may be zero. Each row corresponds to an MCMC iteration.
prior: The prior used to fit the model. If a prior was supplied as an argument it will be returned. Otherwise this will be the automatically generated prior based on the other function arguments.

Arguments

formula: A model formula, as would be passed to glm, specifying the maximal model (i.e. the model with all predictors included).
exposure: A vector of exposure durations matching the length of the response vector. If exposure is of length 1 it will be recycled.
niter: The number of MCMC iterations to run.
data: An optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which poisson.spike is called.
subset: An optional vector specifying a subset of observations to be used in the fitting process.
prior: A list such as that returned by SpikeSlabPrior. If prior is supplied it will be used. Otherwise a prior distribution will be built using the remaining arguments. See SpikeSlabPrior.
na.action: A function which indicates what should happen when the data contain NAs. The default is set by the na.action setting of options, and is na.fail if that is unset. The factory-fresh default is na.omit. Another possible value is NULL, no action. Value na.exclude can be useful.
contrasts: An optional list. See the contrasts.arg of model.matrix.default.
drop.unused.levels: A logical value indicating whether factor levels that are unobserved should be dropped from the model.
initial.value: Initial value for the MCMC algorithm. Can either be a numeric vector, a glm object (from which the coefficients will be used), or a poisson.spike object. If a poisson.spike object is supplied, it is assumed to be from a previous MCMC run for which niter additional draws are desired. If a glm object is supplied then its coefficients will be used as the initial values for the simulation.
ping: If positive, then print a status update to the console every ping MCMC iterations.
nthreads: The number of CPU-threads to use for data augmentation.
seed: Seed to use for the C++ random number generator. It should be NULL or an int. If NULL the seed value will be taken from the global .Random.seed object.
...: Extra arguments to be passed to SpikeSlabPrior.

Author

Steven L. Scott

Details

The MCMC algorithm used here is based on the auxiliary mixture sampling algorithm published by Fruhwirth-Schnatter, Fruhwirth, Held, and Rue (2009).

References

Sylvia Fruhwirth-Schnatter, Rudolf Fruhwirth, Leonhard Held, and Havard Rue. Statistics and Computing, Volume 19 Issue 4, Pages 479-492. December 2009

Examples

Run this code

simulate.poisson.spike <- function(n = 100, p = 10, ngood = 3, niter=1000){
  x <- cbind(1, matrix(rnorm(n * (p-1)), nrow=n))
  beta <- c(rnorm(ngood), rep(0, p - ngood))
  lambda <- exp(x %*% beta)
  y <- rpois(n, lambda)
  x <- x[,-1]
  model <- poisson.spike(y ~ x, niter=niter)
  return(invisible(model))
}
model <- simulate.poisson.spike()
plot(model)
summary(model)

Run the code above in your browser using DataLab