R/plots.R

In dlebauer/pecan-priors: PEcAn module to support prior elicitation

##' Plot univariate response of model output to a trait parameter.
##'
##' @title Sensitivity plot 
##' @param sa.sample trait quantiles used in sensitivity analysis 
##' @param sa.spline spline function estimated from sensitivity analysis
##' @param trait trait name for title
##' @param y.range 
##' @param median.i index of median value in sa.sample; \code{median.i == which(as.numeric(rownames(sa.sample)) == 50) }
##' @param prior.sa.sample similar to sa.sample, but for prior distribution. If given, plots sensitivity for prior run
##' @param prior.sa.spline similar to sa.spline, but for prior trait distribution. 
##' @param fontsize (optional) list with three arguments that can be set to vary the fontsize of the title, axis labels, and axis title in the sensitivity plots
##' @return object of class ggplot
plot.sensitivity <- function(sa.sample, sa.spline, trait,
                             y.range = c(0,50), median.i = 4,
                             prior.sa.sample = NULL, prior.sa.spline = NULL,
                             fontsize = list(title = 24, axis = 18),
                             linesize = 1,
                             dotsize = 2) {
  LENGTH_OUT <- 1000
  
  units <- trait.dictionary(trait)$units
  saplot <- ggplot()

  post.x <- seq(from = min(sa.sample),
                to = max(sa.sample),
                length.out = LENGTH_OUT)
  
  saplot <- saplot + 
    ## plot spline function 
    geom_line(aes(x,y), data = data.frame(x = post.x, y = sa.spline(post.x)), size = linesize) + 
      ## plot points used to evaluate spline
      geom_point(aes(x,y), data = data.frame(x = sa.sample, y = sa.spline(sa.sample)), size = dotsize) +
                                        #indicate median with larger point
        geom_point(aes(x,y), data = data.frame(x = sa.sample[median.i], y = sa.spline(sa.sample[median.i])), size = dotsize * 1.3) + 
          scale_y_continuous(limits = range(pretty(y.range)), breaks = pretty(y.range, n = 3)[1:3]) +
            theme_bw() +
              opts(title= trait.dictionary(trait)$figid, 
                   axis.text.x = theme_text(size = fontsize$axis*1.5),
                   axis.text.y = theme_text(size = fontsize$axis*1.5),
                   axis.title.x = theme_text(size = fontsize$axis*1.5),
                   axis.title.y = theme_blank(),
                   plot.title = theme_text(size = fontsize$title*1.5),
                   panel.border = theme_blank(),
                   panel.grid.minor = theme_blank(),
                   panel.grid.major = theme_blank())
  ## Following conditional can be removed to only plot posterior sa
  prior.x <- post.x
  if(!is.null(prior.sa.sample) & !is.null(prior.sa.spline)){
    prior.x <- seq(from = min(prior.sa.sample), to = max(prior.sa.sample), length.out = LENGTH_OUT)
    saplot <- saplot +
      ## plot spline
      geom_line(aes(x,y), data = data.frame(x = prior.x, y = prior.sa.spline(prior.x)),
                size = linesize, color = 'grey') +
                  ## plot points used to evaluate spline 
                  geom_point(aes(x,y), data = data.frame(x = prior.sa.sample, y = prior.sa.spline(prior.sa.sample)),
                             size = dotsize, color = 'grey') +
                               ## indicate location of medians
                               geom_point(aes(x,y), data = data.frame(x = prior.sa.sample[median.i], 
                                                      y = prior.sa.spline(prior.sa.sample[median.i])),
                                          size = dotsize * 1.5, color = 'grey') 
  }
  max.x <- max(prior.x)
  min.x <- 0
  x.breaks <- pretty(c(min.x, max.x), 4)

  saplot <- saplot + scale_x_continuous(units, limits = range(x.breaks),
                                        breaks = x.breaks)
                                        #  print(saplot)
  return(saplot)
}

plot.variance.decomposition <- function(plot.inputs,
                                        prior.plot.inputs = NULL,
                                        fontsize = list(title = 24, axis = 18)){
  traits    <- names(plot.inputs$variances)
  units     <- trait.dictionary(traits)$units
  trait.labels <- merge(data.frame(id = traits),
                        trait.dictionary(traits),
                        by = 'id', sort = FALSE)$figid
  .plot.data <- data.frame(trait.labels  = trait.labels,
                           units         = units,
                           coef.vars     = plot.inputs$coef.vars * 100,
                           elasticities  = plot.inputs$elasticities,
                           variances     = sqrt(plot.inputs$variances))
                                        #  recover()
  if(!is.null(prior.plot.inputs)) {
    prior.plot.data <- data.frame(trait.labels       = trait.labels,
                                  units              = units,
                                  prior.coef.vars    = prior.plot.inputs$coef.vars * 100,
                                  prior.elasticities = prior.plot.inputs$elasticities,
                                  prior.variances    = sqrt(prior.plot.inputs$variances))
    .plot.data <- merge(.plot.data, prior.plot.data, by = 'trait.labels')
  }
  pv.order <- order(.plot.data$variances, decreasing = FALSE)

  ## location of words and lollipops set by 'points'
  ##    these points can be moved up or down by adjusting the offset X in 1:length(traits) - X
  plot.data <- data.frame(.plot.data[pv.order, ], points = 1:length(traits) - 0.5)
  cv.xticks <<- pretty(as.matrix(plot.data[,grep('coef.var', colnames(plot.data))]), 4)
  pv.xticks <<- pretty(as.matrix(plot.data[,grep('variances', colnames(plot.data))]), 4)  
  el.xticks <<- pretty(as.matrix(plot.data[,grep('elasticities', colnames(plot.data))]), 3)
  el.xrange <<- range(pretty(as.matrix(plot.data[,grep('elasticities', colnames(plot.data))]), 4))

  
  ## Notes on fine-tuning plots below
  ## axis lines and ticks drawn for each plot using geom_segment  
  ## size of x axis tick set by xend = ...
  ## vertical location of axis numbers set in base.plot using vjust
  
  base.plot <- ggplot(plot.data) +
    coord_flip() +
      theme_bw() +
        opts(axis.line.y = theme_blank(),
             axis.text.x = theme_text(size=fontsize$axis, vjust = -1),
             axis.text.y = theme_blank(),
             axis.title.x = theme_blank(), 
             axis.title.y = theme_blank(),
             axis.ticks = theme_blank(),
             panel.grid.major = theme_blank(),
             panel.grid.minor = theme_blank(),
             panel.border = theme_blank())

  if(!is.null(prior.plot.inputs)) {
    .cv.plot <-  base.plot +
      geom_pointrange(aes(x = points, y = prior.coef.vars,
                          ymin = 0, ymax = prior.coef.vars),
                      size = 1.25, color = 'grey')
    
    .el.plot <- base.plot +
      geom_pointrange(aes(x = points, prior.elasticities,
                          ymin = 0, ymax = prior.elasticities),
                      size = 1.25, color = 'grey') 

    .pv.plot <- base.plot +
      geom_pointrange(aes(x = points, y = prior.variances,
                          ymin = 0, ymax = prior.variances),
                      size = 1.25, color = 'grey') 
  } else {
    .cv.plot <- base.plot + scale_y_continuous(breaks = cv.xticks)
    .el.plot <- base.plot + scale_y_continuous(breaks = el.xrange)
    .pv.plot <- base.plot + scale_y_continuous(breaks = pv.xticks)
  }


  trait.plot <- base.plot + 
    opts(title = 'Parameter',
         plot.title = theme_text(hjust = 0.96, size = fontsize$title),
         axis.text.x = theme_text(colour='white'),
         axis.line.x = theme_blank()) +
           geom_text(aes(y = 1, x = points,
                         label=trait.labels, hjust = 1),
                     size = fontsize$axis/3) +
                       scale_y_continuous( breaks = c(0,0), limits = c(0,1)) +
                         ##  Add Invisible Axes to resize like other plots
                         geom_segment(aes(x = c(0,0), y = c(0,0),
                                          yend = c(0, max(cv.xticks)),
                                          xend = c(length(traits), 0)), colour = 'white')  + 
                                            ## Add invisible ticks
                                            geom_segment(aes(x = 0,
                                                             y = cv.xticks,
                                                             xend = -0.1,
                                                             yend = cv.xticks), colour = 'white') 

  cv.plot <- .cv.plot +
    opts(title = 'CV (%)', plot.title = theme_text(size = fontsize$title)) +
      scale_y_continuous(breaks = cv.xticks, limits = range(cv.xticks)) +
        geom_pointrange(aes(x = points, y = coef.vars, ymin = 0, ymax = coef.vars),
                        size = 1.25) + 
                          ##  Add Axes
                          geom_segment(aes(x = c(0,0), y = c(0,0),
                                           yend = c(0, max(cv.xticks)),
                                           xend = c(length(traits), 0)))  + 
                                             ## Add Ticks
                                             geom_segment(aes(x = 0,
                                                              y = cv.xticks,
                                                              xend = -0.1,
                                                              yend = cv.xticks))


  if (diff(range(el.xticks)) < 4) el.xticks <- c(-1,0,1)
  el.plot <- .el.plot + 
    opts(title = 'Elasticity', plot.title = theme_text(size = fontsize$title)) +
      scale_y_continuous(breaks = el.xticks, limits = range(el.xrange)) +
        geom_pointrange(aes(x = points, y = elasticities, ymin = 0, ymax = elasticities),
                        size = 1.25) +
                          ##  Add Axes
                          geom_segment(aes(x = c(0,0), y = c(0, min(el.xrange)),
                                           yend = c(0, max(el.xrange)),
                                           xend = c(length(traits), 0)))  +
                                             ## Add Ticks
                                             geom_segment(aes(x = 0,
                                                              y = el.xticks,
                                                              xend = -0.1,
                                                              yend = el.xticks)) 

  pv.plot <- .pv.plot + 
    opts(title = 'Root Variance (SD)',
         plot.title = theme_text(size = fontsize$title)) +
           scale_y_continuous(breaks = pv.xticks, limits = range(pv.xticks)) +
             geom_pointrange(aes(x = points, variances,
                                 ymin = 0, ymax = variances), size = 1.25) +
                                   ##  Add Axes
                                   geom_segment(aes(x = c(0,0), y = c(0,0),
                                                    yend = c(0, max(pv.xticks)),
                                                    xend = c(length(traits), 0)))  + 
                                                      ## Add Ticks
                                                      geom_segment(aes(x = 0,
                                                                       y = pv.xticks,
                                                                       xend = -0.1,
                                                                       yend = pv.xticks)) 
  
  
  return(list(trait.plot = trait.plot, cv.plot = cv.plot, el.plot = el.plot, pv.plot = pv.plot))


}

##' Plot functions and quantiles used in sensitivity analysis
##'
##' 
##' @title Plot Sensitivities
##' @param sensitivity.plot.inputs 
##' @param prior.sensitivity.plot.inputs 
##' @param ... arguments passed to \code{\link{plot.sensitivity}}
##' @param sensitivity.results list containing sa.samples and sa.splines 
##' @return list of plots, one per trait
plot.sensitivities <- function(sensitivity.plot.inputs, prior.sensitivity.plot.inputs = NULL, ...){
  sa.samples <- sensitivity.plot.inputs$sa.samples
  sa.splines <- sensitivity.plot.inputs$sa.splines
  if(!is.null(prior.sensitivity.plot.inputs)) {
    prior.sa.samples <- prior.sensitivity.plot.inputs$sa.samples
    prior.sa.splines <- prior.sensitivity.plot.inputs$sa.splines
  }
  traits <- names(sa.samples)

  y.range <- c(0, max(mapply(do.call, sa.splines, lapply(sa.samples, list)), na.rm = TRUE))

  sensitivity.plots <- list()
  for(trait in traits) {
    if(!is.null(prior.sensitivity.plot.inputs)) {
      prior.sa.sample <- prior.sa.samples[,trait]
      prior.sa.spline <- prior.sa.splines[[trait]]
    } else {
      prior.sa.sample <- NULL
      prior.sa.spline <- NULL
    }
    sensitivity.plots[[trait]] <-plot.sensitivity(sa.sample =  sa.samples[,trait],
                                                  sa.spline = sa.splines[[trait]],
                                                  trait <- trait,
                                                  y.range = y.range,
                                                  median.i =  4,#which(as.numeric(rownames(sa.samples)) == 50),
                                                  prior.sa.sample = prior.sa.sample,
                                                  prior.sa.spline = prior.sa.spline,
                                                  ...)
  }
  return(sensitivity.plots)
}

##' Variable-width (dagonally cut) histogram
##'
##' 
##' When constructing a histogram, it is common to make all bars the same width.
##' One could also choose to make them all have the same area.
##' These two options have complementary strengths and weaknesses; the equal-width histogram oversmooths in regions of high density, and is poor at identifying sharp peaks; the equal-area histogram oversmooths in regions of low density, and so does not identify outliers.
##' We describe a compromise approach which avoids both of these defects. We regard the histogram as an exploratory device, rather than as an estimate of a density. 
##' @title Diagonally Cut Histogram 
##' @param x is a numeric vector (the data)
##' @param a is the scaling factor, default is 5 * IQR
##' @param nbins is the number of bins, default is assigned by the Stuges method
##' @param rx  is the range used for the left of the left-most bin to the right of the right-most bin  
##' @param eps used to set artificial bound on min width / max height of bins as described in Denby and Mallows (2009) on page 24.
##' @param xlab is label for the x axis 
##' @param plot = TRUE produces the plot, FALSE returns the heights, breaks and counts
##' @param lab.spikes = TRUE labels the \% of data in the spikes
##' @return list with two elements, heights of length n and breaks of length n+1 indicating the heights and break points of the histogram bars. 
##' @author Lorraine Denby, Colin Mallows
##' @references Lorraine Denby, Colin Mallows. Journal of Computational and Graphical Statistics. March 1, 2009, 18(1): 21-31. doi:10.1198/jcgs.2009.0002.
dhist<-function(x, a=5*iqr(x),
                nbins=nclass.Sturges(x), rx = range(x,na.rm = TRUE),
                eps=.15, xlab = "x", plot = TRUE,lab.spikes = TRUE)
{

  if(is.character(nbins))
    nbins <- switch(casefold(nbins),
                    sturges = nclass.Sturges(x),
                    fd = nclass.FD(x),
                    scott = nclass.scott(x),
                    stop("Nclass method not recognized"))
  else if(is.function(nbins))
    nbins <- nbins(x)

  x <- sort(x[!is.na(x)])
  if(a == 0)
    a <- diff(range(x))/100000000
  if(a != 0 & a != Inf) {
    n <- length(x)
    h <- (rx[2] + a - rx[1])/nbins
    ybr <- rx[1] + h * (0:nbins)
    yupper <- x + (a * (1:n))/n
                                        # upper and lower corners in the ecdf
    ylower <- yupper - a/n
                                        #
    cmtx <- cbind(cut(yupper, breaks = ybr), cut(yupper, breaks = 
                                ybr, left.include = TRUE), cut(ylower, breaks = ybr),
                  cut(ylower, breaks = ybr, left.include = TRUE))
    cmtx[1, 3] <- cmtx[1, 4] <- 1
                                        # to replace NAs when default r is used
    cmtx[n, 1] <- cmtx[n, 2] <- nbins
                                        #
                                        #checksum <- apply(cmtx, 1, sum) %% 4
    checksum <- (cmtx[, 1] + cmtx[, 2] + cmtx[, 3] + cmtx[, 4]) %%
    4
                                        # will be 2 for obs. that straddle two bins
    straddlers <- (1:n)[checksum == 2]
                                        # to allow for zero counts
    if(length(straddlers) > 0) {
      counts <- table(c(1:nbins, cmtx[ - straddlers, 1]))
    } else {
      counts <- table(c(1:nbins, cmtx[, 1]))
    }
    counts <- counts - 1
                                        #
    if(length(straddlers) > 0) {
      for(i in straddlers) {
        binno <- cmtx[i, 1]
        theta <- ((yupper[i] - ybr[binno]) * n)/a
        counts[binno - 1] <- counts[binno - 1] + (
                                                  1 - theta)
        counts[binno] <- counts[binno] + theta
      }
    }
    xbr <- ybr
    xbr[-1] <- ybr[-1] - (a * cumsum(counts))/n
    spike<-eps*diff(rx)/nbins
    flag.vec<-c(diff(xbr)<spike,F)
    if ( sum(abs(diff(xbr))<=spike) >1) {
      xbr.new<-xbr
      counts.new<-counts
      diff.xbr<-abs(diff(xbr))
      amt.spike<-diff.xbr[length(diff.xbr)]
      for (i in rev(2:length(diff.xbr))) {
        if (diff.xbr[i-1]<=spike&diff.xbr[i]<=spike&
            !is.na(diff.xbr[i])) {
          amt.spike<-amt.spike+diff.xbr[i-1]
          counts.new[i-1]<-counts.new[i-1]+counts.new[i]
          xbr.new[i]<-NA
          counts.new[i]<-NA
          flag.vec[i-1]<-T
        }
        else amt.spike<-diff.xbr[i-1]
      }
      flag.vec<-flag.vec[!is.na(xbr.new)]
      flag.vec<-flag.vec[-length(flag.vec)]
      counts<-counts.new[!is.na(counts.new)]
      xbr<-xbr.new[!is.na(xbr.new)]

    }
    else flag.vec<-flag.vec[-length(flag.vec)]
    widths <- abs(diff(xbr))
    ## N.B. argument "widths" in barplot must be xbr
    heights <- counts/widths
  }
  bin.size <- length(x)/nbins
  cut.pt <- unique(c(min(x) - abs(min(x))/1000,
                     approx(seq(length(x)), x, (1:(nbins - 1)) * bin.size, rule = 2)$y, max(x)))
  aa <- hist(x, breaks = cut.pt, plot = FALSE, probability = TRUE)
  if(a == Inf) {
    heights <- aa$counts
    xbr <- aa$breaks
  }
  amt.height<-3
  q75<-quantile(heights,.75)
  if (sum(flag.vec)!=0) {
    amt<-max(heights[!flag.vec])
    ylim.height<-amt*amt.height
    ind.h<-flag.vec&heights> ylim.height
    flag.vec[heights<ylim.height*(amt.height-1)/amt.height]<-F
    heights[ind.h] <- ylim.height
  }
  amt.txt<-0
  end.y<-(-10000)
  if(plot) {
    barplot(heights, abs(diff(xbr)), space = 0, density = -1, xlab = 
            xlab, plot = TRUE, xaxt = "n",yaxt='n')
    at <- pretty(xbr)
    axis(1, at = at - xbr[1], labels = as.character(at))
    if (lab.spikes) {
      if (sum(flag.vec)>=1) {
        usr<-par('usr')
        for ( i in seq(length(xbr)-1)) {
          if (!flag.vec[i]) {
            amt.txt<-0
            if (xbr[i]-xbr[1]<end.y) amt.txt<-1
          }
          else {
            amt.txt<-amt.txt+1
            end.y<-xbr[i]-xbr[1]+3*par('cxy')[1]
          }
          if (flag.vec[i]) {
            txt<-paste(' ',format(round(counts[i]/
                                        sum(counts)*100)),'%',sep='')
            par(xpd = TRUE)
            text(xbr[i+1]-xbr[1],ylim.height-par('cxy')[2]*(amt.txt-1),txt, adj=0)
          }}
      }
      else print('no spikes or more than one spike')
    }
    invisible(list(heights = heights, xbr = xbr))
  }
  else {
    return(list(heights = heights, xbr = xbr,counts=counts))
  }
}

##' Calculate interquartile range
##'
##' Calculates the 25th and 75th quantiles given a vector x; used in function \link{dhist}.
##' @title Interquartile range
##' @param x vector
##' @return numeric vector of length 2, with the 25th and 75th quantiles of input vector x. 
iqr<-function(x){
  return(diff(quantile(x, c(0.25, 0.75), na.rm = TRUE)))
}

create.base.plot <- function() {
  base.plot <- ggplot()
  return(base.plot)
}
##' Plots a prior density from a parameterized probability distribution  
##'
##' @title Add Prior Density
##' @param prior.density 
##' @param base.plot a ggplot object (grob), created by \code{\link{create.base.plot}} if none provided
##' @param prior.color color of line to be plotted
##' @export
##' @return plot with prior density added
##' @seealso \code{\link{pr.dens}}
##' @examples
##' add.prior.density(c('norm', 0, 1))
add.prior.density <- function(prior.density, base.plot = NULL, prior.color = 'black' ) {
  if(is.null(base.plot)) base.plot <- create.base.plot()
  new.plot <- base.plot +  geom_line(data = prior.density,
                                     aes(x = x, y = y),
                                     color = prior.color)
  return(new.plot)
}

##' Returns a data frame from \link{stats::density} function 
##'
##' @title Create Density Data Frame from Sample
##' @param sample 
##' @param ... additional arguments to density 
##' @return data frame with x and y
create.density.df <- function(samps = NULL,
                              zero.bounded = FALSE,
                              distribution = NULL) {
  samp.exists <- !is.null(samps)
  dist.exists <- !is.null(distribution)
  if(identical(samp.exists, dist.exists)){
    stop('create.density.df requires one and only one of:
             samps: a vector of samples, e.g. MCMC chain,
               OR
             distribution: a named distribution supported by R')
  }
  if(samp.exists){
    if(zero.bounded) {
      new.density <- zero.bounded.density(samps)
    } else {    
      new.density <- density(samps)
    }
    density.df <- with(new.density,
                       data.frame(x = x,
                                  y = y))
  }
  
  if(dist.exists) {
    density.df <- do.call(pr.dens, c(distribution[1:3], n = 1000))
  }
  return(density.df)
}

##'  Add posterior density to a plot
##'
##' @title Add posterior density. 
##' @param posterior.density 
##' @param base.plot a ggplot object (grob), created by \code{\link{create.base.plot}} if none provided
##' @return plot with posterior density line added
add.posterior.density <- function(posterior.density, base.plot = NULL) {
  if(is.null(base.plot)) base.plot <- create.base.plot()
  new.plot <- base.plot +  geom_line(data = posterior.density,
                                     aes(x = x, y = y))
  return(new.plot)  
}

##' Add data to an existing plot or create a new one from \code{\link{create.base.plot}}
##'
##' Used to add raw data or summary statistics to the plot of a distribution.
##' The height of Y is arbitrary, and can be set to optimize visualization.
##' If SE estimates are available, tehse wil be plotted
##' @title Add data to plot 
##' @param trait.data data to be plotted
##' @param base.plot a ggplot object (grob), created by \code{\link{create.base.plot}} if none provided
##' @param ymax maximum height of y
##' @seealso \code{\link{create.base.plot}}
##' @return updated plot object
##' @examples
##' add.data(data.frame(Y = c(1, 2), se = c(1,2)), base.plot = NULL, ymax = 10)
add.data <- function(trait.data, base.plot = NULL, ymax, color = 'black') {
  if(is.null(base.plot)) base.plot <- create.base.plot()
  n.pts <- nrow(trait.data)
  ymax <- ymax * sqrt(n.pts + 1)/sqrt((n.pts + 1) * 64)
  y.pts <- seq(0, ymax, length.out = 1 + n.pts)[-1]
  plot.data <- data.frame(x = trait.data$Y,
                          y = y.pts,
                          se = trait.data$se)
  new.plot <- base.plot +
    geom_point(data = plot.data,
               aes(x = x, y = y),
               color = color) +
                 geom_segment(data = plot.data,
                              aes(x = x - se, y = y, xend = x + se, yend = y),
                              color = color)
  return(new.plot)
}

##' Plot trait density and data
##'
##' @title Plot trait density
##' @param trait character, name of trait to be plotted
##' @param prior named distribution with parameters
##' @param posterior.sample 
##' @param trait.data 
##' @param fontsize 
##' @return plot (grob) object
##' @examples
##' prior1 <- data.frame(distn = 'norm',
##'                      parama = 20,
##'                      paramb = 5)
##' data1  <- data.frame(Y = c(19, 21), se = c(1,1))
##' plot.trait(trait = 'Vcmax',
##'           prior = prior1,
##'           trait.data = data1)
plot.trait <- function(trait,
                       prior = NULL,
                       posterior.sample = NULL,
                       trait.df = NULL,
                       fontsize = list(title = 24, axis = 18),
                       x.lim = NULL,
                       y.lim = NULL,
                       logx = FALSE) {
  ## Determine plot components
  plot.posterior <- !is.null(posterior.sample)
  plot.prior     <- !is.null(prior)
  plot.data      <- !is.null(trait.df)
  
  ## If no posterior, set prior density to black

  ## get units for plot title
  units <- trait.dictionary(trait)$units

  prior.density <- posterior.density <- data.frame(x = NA, y = NA)
  if(plot.prior){
    prior.color = 'black'
    prior.density   <- create.density.df(distribution = prior)
    prior.density <- prior.density[prior.density$x > 0, ]
  }
  if(plot.posterior){
    posterior.density <- create.density.df(samps = posterior.sample)
    posterior.density <- posterior.density[posterior.density$x > 0, ]
    prior.color <-  'grey'
  }

  if(is.null(x.lim)){
    x.lim <- range(posterior.density$x, prior.density$x, na.rm = TRUE)
  }
  if(is.null(y.lim)){
    y.lim <- range(posterior.density$y, prior.density$y, na.rm = TRUE)
  }
  print(trait)
  print(x.lim)
  print(y.lim)
  ticks <<- lapply(data.frame(x = x.lim, y = y.lim), pretty, 4)

  print(ticks)
  x.lim <<- range(ticks$x)
  y.lim <<- range(ticks$y) 


  base.plot <- create.base.plot() + theme_bw()
  if(logx == TRUE){
    ticks$x <- ticks$x[ticks$x > 0]
    base.plot <- base.plot + scale_x_log10(units,
                                           breaks = ticks$x,
                                           labels = ticks$x,
                                           limits = range(ticks$x))
  } else {
    base.plot <- base.plot + scale_x_continuous(units,
                                                breaks = ticks$x,
                                                labels = ticks$x,
                                                limits = range(ticks$x))
  }
  
  if(plot.prior){
    keep <- (prior.density$x > x.lim[1] &
             prior.density$x < x.lim[2] &
             prior.density$y > y.lim[1] &
             prior.density$y < y.lim[2])
    
    prior.density <- prior.density[keep,]
    base.plot <- add.prior.density(prior.density,
                                   base.plot = base.plot,
                                   prior.color = prior.color)
  }
  if(plot.posterior){
    keep <- (posterior.density$x > x.lim[1] &
             posterior.density$x < x.lim[2] &
             posterior.density$y > y.lim[1] &
             posterior.density$y < y.lim[2])
    posterior.density <- posterior.density[keep, ]
    base.plot <- add.posterior.density(posterior.density,
                                       base.plot = base.plot)
  }
  if(plot.data){
    base.plot <- add.data(trait.df,
                          base.plot = base.plot,
                          ymax = y.lim[2])
  }

  trait.plot <- base.plot +
    opts(title= trait.dictionary(trait)$figid, 
         axis.text.x = theme_text(size = fontsize$axis),
         axis.text.y = theme_blank(),
         axis.title.x = theme_text(size = fontsize$axis),
         axis.title.y = theme_blank(),
         axis.ticks = theme_blank(),
         axis.line = theme_blank(),
         plot.title = theme_text(size = fontsize$title),
         panel.grid.major = theme_blank(),
         panel.grid.minor = theme_blank(),
         panel.border = theme_blank()) +
           geom_segment(aes(x = 0, y = 0, yend = 0, xend = x.lim[2])) +
             geom_segment(aes(x = pretty(x.lim, 4),
                              y = 0,
                              xend = pretty(x.lim, 4),
                              yend = -y.lim[2]/50)) 
  return(trait.plot)
}


##' Plot probability density and data
##'
##' @title Plot Trait Probability Densities
##' @param sensitivity.results list containing sa.samples and sa.splines 
##' @param outdir 
##' @return outputs plots in outdir/sensitivity.analysis.pdf file 
plot.densities <- function(density.plot.inputs, outdir, ...){
  trait.samples          <- density.plot.inputs$trait.samples
  trait.df               <- density.plot.inputs$trait.df
  prior.trait.samples    <- density.plot.inputs$trait.df
  
  traits <- names(trait.samples)
  pdf(paste(outdir, 'trait.densities.pdf', sep=''), height = 12, width = 20)

  for(trait in traits) {
    density.plot <- plot.density(trait.sample       =  trait.samples[,trait],
                                 trait.df         = trait.df[[trait]],
                                 ...)
    print(sensitivity.plot)
  }
  dev.off()
}


##' Calculate the density of a distribution for use in plotting
##'
##' @title Prior Density 
##' @param distribution one of R's supported distributions (character)
##' @param a first parameter of distribution (numeric)
##' @param b second parameter of distribution (numeric)
##' @return data frame with values of x, the density at each x and the probability at each x
##' @author David LeBauer
prior.density <- function(distribution = 'norm', a = 0, b = 1, xlim = NA){
  distribution <- gsub('lognormal', 'lnorm', distribution)
  if(distribution != 'beta'){
    if(isTRUE(is.na(xlim))){
      range.x <- range(pretty(do.call(paste('q', distribution, sep = ''), list(c(0.005, 0.995),a,b))))
    } else if (isTRUE(!is.na(xlim))) {
      range.x <- xlim
    }
    prior.x <- seq(from=range.x[1], to = range.x[2], length = 1000)
  } else {
    range.x <- c(0,1)
    prior.x <- seq(from=0, to = 1, length = 1000)  
  }
  dens.x  <- do.call(paste('d', distribution, sep=''),list(prior.x, a, b))
  prob.x  <- do.call(paste('p', distribution, sep=''),list(prior.x, a, b))
  return(data.frame(prior.x, dens.x, prob.x))
}

##' Plot prior density and data
##'
##' @title Prior Figure 
##' @param priordata observations to be plotted as points
##' @param priordensity density of prior distribution, calculated by \code{\link{prior.density}}
##' @param trait name of trait
##' @param xlim limits for x axis
##' @return plot / grob of prior distribution with data used to inform the distribution 
priorfig <- function(priordata = NA, priordensity = NA, trait = '', xlim = 'auto', fontsize = 24){
  if(is.data.frame(priordata)){
    colnames(priordata) <- 'x'
  }

  if(isTRUE(xlim == 'auto')) {
    x.breaks <- pretty(c(signif(priordensity$prior.x, 2)), 4)
    xlim <- range(x.breaks)
  } else {
    x.breaks <- pretty(signif(xlim, 2), 4)
    xlim <- range(c(x.breaks, xlim))
  }

  priorfigure <- ggplot() + theme_bw() + 
    scale_x_continuous(limits = xlim, breaks = x.breaks, trait.dictionary(trait)$units) +
      scale_y_continuous(breaks=NA)+
        opts(title = trait.dictionary(trait)$figid,
             panel.grid.major = theme_blank(),    
             panel.grid.minor = theme_blank(),
             axis.text.y = theme_blank(),
             axis.text.x = theme_text(size= fontsize),
             axis.title.y = theme_blank(), ## hide y axis label
             axis.title.x = theme_text(size = fontsize * 0.9), ## hide y axis label
             plot.title = theme_text(size = fontsize*1.1),
             ## plot.margin = unit(c(0.4, 0.1, 0.1, 0.1), 'lines'), 
             panel.border = theme_border(c('bottom'))
             ) 
  
  
  if(is.data.frame(priordata)){
    priordata <- subset(priordata, subset = !is.na(x))
    dx <- with(priordata,
               min(abs(diff(x)[diff(x)!=0])))
    priordata <- transform(priordata,
                           x = x + runif(length(x), -dx/2, dx/2))# this adds jitter to separate equal values 
    rug <- geom_rug(data = priordata, aes(x))

                                        #rug <- geom_point(data = priordata, aes(x=x, y = seq(0,max(priordensity$dens.x)/30, length = length(x))), size = 3, alpha = 3/sqrt(nrow(priordata)))
                                        #hist <-  geom_histogram(data = priordata, aes(x=x, y = ..density..),  alpha = 0.5, binwidth = diff(range(priordata))/sqrt(nrow(priordata)))
    priorfigure <- priorfigure + rug
  } 
  if(is.data.frame(priordensity[1])){
    dens.line <- geom_line(data=priordensity, aes(x=prior.x, y=dens.x))
    qpts <- get.quantiles.from.density(priordensity)
    dens.ci <- geom_point(data = qpts, aes(x,y))
    priorfigure <- priorfigure + dens.line + dens.ci
  }
  return(priorfigure)
} 


get.quantiles.from.density <- function(priordensity){
  qi <- c(which.min(abs(priordensity$prob.x - 0.025)),
          which.min(abs(priordensity$prob.x - 0.5)),
          which.min(abs(priordensity$prob.x - 0.975)))
  qs <- priordensity[qi,c('prior.x', 'dens.x')]
  colnames(qs) <- c('x', 'y')
  return(qs)
}
##' .. content for \description{} (no empty lines) ..
##'
##' .. content for \details{} ..
##' @title panel.border
##' @param type 
##' @param colour 
##' @param size 
##' @param linetype 
##' @return a plot border
##' @author Rudolf Cardinal http://goo.gl/YThRT
theme_border <- function(type = c("left", "right", "bottom", "top", 
                           "none"), colour = "black", size = 1, linetype = 1) {
  type <- match.arg(type, several.ok=TRUE) 
  structure( 
            function(x = 0, y = 0, width = 1, height = 1, ...) { 
              xlist <- c() 
              ylist <- c() 
              idlist <- c() 
              if ("bottom" %in% type) { # bottom 
                xlist <- append(xlist, c(x, x+width)) 
                ylist <- append(ylist, c(y, y)) 
                idlist <- append(idlist, c(1,1)) 
              } 
              if ("top" %in% type) { # top 
                xlist <- append(xlist, c(x, x+width)) 
                ylist <- append(ylist, c(y+height, y+height)) 
                idlist <- append(idlist, c(2,2)) 
              } 
              if ("left" %in% type) { # left 
                xlist <- append(xlist, c(x, x)) 
                ylist <- append(ylist, c(y, y+height)) 
                idlist <- append(idlist, c(3,3)) 
              } 
              if ("right" %in% type) { # right 
                xlist <- append(xlist, c(x+width, x+width)) 
                ylist <- append(ylist, c(y, y+height)) 
                idlist <- append(idlist, c(4,4)) 
              } 
              polylineGrob( 
                           x=xlist, y=ylist, id=idlist, ..., default.units = "npc", 
                           gp=gpar(lwd=size, col=colour, lty=linetype), 
                           ) 
            }, 
            class = "theme", 
            type = "box", 
            call = match.call() 
            ) 
}