BoomSpikeSlab (version 0.5.2)

plot.marginal.inclusion.probabilities: 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.

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)