plot.marginal.inclusion.probabilities {BoomSpikeSlab}

Plot marginal inclusion probabilities.

Description

Produces a barplot of the marginal inclusion probabilities for a set of model coefficients sampled under a spike and slab prior. The coefficients are sorted by the marginal inclusion probability, and shaded by the conditional probability that a coefficient is positive, given that it is nonzero.

Usage

PlotMarginalInclusionProbabilities(
     beta,
     burn = 0,
     inclusion.threshold = 0,
     unit.scale = TRUE,
     number.of.variables = NULL,
     ...)

Arguments

beta

A matrix of model coefficients. Each row represents an MCMC draw. Each column represents a coefficient for a variable.

burn

The number of MCMC iterations in the ojbect to be discarded as burn-in.

inclusion.threshold

Only plot coefficients with posterior inclusion probabilities exceeding this value.

unit.scale

A logical value indicating whether the scale of the plot should be from 0 to 1. Otherwise the scale is determined by the maximum inclusion probability.

number.of.variables

If non-NULL this specifies the number of coefficients to plot, taking precedence over inclusion.threshold.

...

Additional arguments to be passed to barplot.

Value

Invisibly returns a list with the following elements.

barplot

The midpoints of each bar, which is useful for adding to the plot.

inclusion.prob

The marginal inclusion probabilities of each variable, ordered smallest to largest (the same order as the plot).

positive.prob

The probability that each variable has a positive coefficient, in the same order as inclusion.prob.

permutation

The permutation of beta that puts the coefficients in the same order as positive.prob and inclusion.prob. That is: beta[, permutation] will have the most significant coefficients in the right hand columns.

Author(s)

Steven L. Scott

See Also

lm.spike SpikeSlabPrior summary.lm.spike predict.lm.spike

Examples

simulate.lm.spike <- function(n = 100, p = 10, ngood = 3, niter=1000, sigma = 8){
  x <- cbind(matrix(rnorm(n * (p-1)), nrow=n))
  beta <- c(rnorm(ngood), rep(0, p - ngood))
  y <- rnorm(n, beta[1] + x %*% beta[-1], sigma)
  draws <- lm.spike(y ~ x, niter=niter)
  return(invisible(draws))
}
model <- simulate.lm.spike(n = 1000, p = 50, sigma = .3)
plot(model, inclusion.threshold = .01)