R/bmaImage.R

In DataScienceSalon/Bayesian-Regression: Bayesian Linear Regression and Default Priors

#' Images of models used in Bayesian model averaging
#'
#' Creates an image of the models selected using \code{\link{bas}}.
#'
#' Creates an image of the model space sampled using \code{\link{bas}}.  If a
#' subset of the top models are plotted, then probabilities are renormalized
#' over the subset.
#'
#'
#' @param x A BMA object of type 'bas' created by BAS
#' @param top.models Number of the top ranked models to plot
#' @param intensity Logical variable, when TRUE image intensity is proportional
#' to the probability or log(probability) of the model, when FALSE, intensity
#' is binary indicating just presence (light) or absence (dark) of a variable.
#' @param prob Logical variable for whether the area in the image for each
#' model should be proportional to the posterior probability (or log
#' probability) of the model (TRUE) or with equal area (FALSE).
#' @param log Logical variable indicating whether the intensities should be
#' based on log posterior odds (TRUE) or posterior probabilities (FALSE).  The
#' log of the posterior odds is for comparing the each model to the worst model
#' in the top.models.
#' @param rotate Should the image of models be rotated so that models are on
#' the y-axis and variables are on the x-axis (TRUE)
#' @param color The color scheme for image intensities. The value "rainbow"
#' uses the rainbow palette. The value "blackandwhite" produces a black and
#' white image (greyscale image)
#' @param subset indices of variables to include in plot; 1 is the intercept
#' @param offset numeric value to add to intensity
#' @param digits number of digits in posterior probabilities to keep
#' @param vlas las parameter for placing variable names; see par
#' @param plas las parameter for posterior probability axis
#' @param rlas las parameter for model ranks
#' @param ... Other parameters to be passed to the \code{image} and \code{axis}
#' functions.
#' @note Suggestion to allow area of models be proportional to posterior
#' probability due to Thomas Lumley
#' @author Merlise Clyde \email{clyde@@stat.duke.edu}
#' @seealso \code{\link{bas}}
#' @references Clyde, M. (1999) Bayesian Model Averaging and Model Search
#' Strategies (with discussion). In Bayesian Statistics 6. J.M. Bernardo, A.P.
#' Dawid, J.O. Berger, and A.F.M. Smith eds. Oxford University Press, pages
#' 157-185.
#' @keywords regression
#' @examples
#'
#' require(graphics)
#' data("Hald")
#' hald.ZSprior =  bas.lm(Y~ ., data=Hald,  prior="ZS-null")
#' image(hald.ZSprior, subset=-1)
#'
bmaImage <- function (x, top.models=20, intensity=TRUE, prob=TRUE, log=TRUE, rotate=TRUE, color="rainbow", subset=NULL, offset=.75, digits=3, vlas=2,plas=0,rlas=0, ...)
{
  postprob = x$postprobs
  top.models = min(top.models, x$n.models)
  best = order(-x$postprobs)[1:top.models]
  postprob=postprob[best]/sum(postprob[best])
  which.mat <-  list2matrix.which(x, best)
  nvar <- ncol(which.mat)
  
  
  if (is.null(subset)) subset=1:nvar
  
  which.mat =  which.mat[,subset, drop=FALSE]
  nvar = ncol(which.mat)
  namesx = x$namesx[subset]
  
  scale = postprob
  prob.lab= "Posterior Probability"
  
  if (log)   {
    scale = log(postprob) - min(log(postprob))
    prob.lab="Log Posterior Odds"
  }
  
  if (intensity)  which.mat = sweep(which.mat, 1, scale+offset,"*")
  
  if (rotate) scale = rev(scale)
  
  if (prob) m.scale = cumsum(c(0,scale))
  else  m.scale = seq(0, top.models)
  
  mat = (m.scale[-1] +  m.scale[-(top.models+1)])/2
  
  colors = switch(color,
                  "rainbow" = c("black", rainbow(top.models+1, start=.75, end=.05)),
                  "blackandwhite" =  gray(seq(0, 1, length=top.models))
  )
  
  par.old = par()$mar
  
  if (rotate) {
    par(mar = c(6,10,6,2) + .1, cex.main = 2)
    image(0:nvar, mat, t(which.mat[top.models:1,]),
          xaxt="n", yaxt="n",
          ylab="",
          xlab="",
          zlim=c(0, max(which.mat)),
          col=colors, ...)
    
    axis(2,at=mat,labels=round(scale, digits=digits), las=plas, ...)
    axis(4,at=mat,labels=top.models:1, las=rlas, ...)
    #mtext("Model Rank", side=4, line=3, las=0)
    mtext(prob.lab, side=2, line=4, las=0)
    axis(1,at=(1:nvar -.5), labels=namesx, las=vlas, ...)
  }
  else{
    par(mar = c(6,10,6,2) + .1, cex.main = 2)
    image(mat, 0:nvar, which.mat[ , nvar:1],
          xaxt="n", yaxt="n",
          xlab="",
          ylab="",
          zlim=c(0, max(which.mat)),
          col=colors, ...)
    
    axis(1,at=mat,labels=round(scale, digits=digits), las=plas,...)
    axis(3,at=mat,labels=1:top.models, las=rlas, ...)
    #mtext("Model Rank", side=3, line=3)
    mtext(prob.lab, side=1, line=4)
    axis(2,at=(1:nvar -.5), labels=rev(namesx), las=vlas, ...)
  }
  
  box()
  
  par(par.old)
  invisible()
}