R/ppd_plot.r

In GmAMisc: 'Gianmarco Alberti' Miscellaneous

#' R function for plotting Posterior Probability Densities for Bayesian modeled 14C dates/parameters
#'
#' The function allows plot Posterior Probability Densities with a nice outlook thanks to
#' 'ggplot2'.\cr It takes as input a dataframe that must be organized as follows (it is rather easy
#' to do that once the data have been exported from OxCal):\cr -calendar dates (first column to the
#' left);\cr -posterior probabilities (second column);\cr -grouping variables (third column), which
#' could contain the names of the events of interest (e.g., phase 1 start, phase 1 end, phase 2
#' start, phase 2 end, etc).
#' @param data Dataframe containing the data as returned by the OxCal program.
#' @param lower Lower limit of the calendar date axis.
#' @param upper Upper limit of the calendar date axis; if the lower and upper parameters are not
#'   provided, the default values will be the earliest and latest calendar dates.
#' @param type Type of plot the user wishes to plot (a: curves outlined by a line; b: curves plotted
#'   as solid areas; c: combination of a and b).
#'
#' @keywords ppdPlot
#'
#'
#' @export
#'
#' @examples
#' #load a toy dataset
#' data(radioc_data)
#'
#' #plot the Posterior Probability Densities for the phases' parameters
#' ppdPlot(radioc_data, type="a")
#'
#' #plot the Posterior Probability Densities for the phases' parameters,
#' # setting different boundaries for the x-axis and using filled curves instead of simple outlines
#' ppdPlot(radioc_data, -1000, 100, type="b")
#'
ppdPlot <- function(data,lower=min(data[,1]),upper=max(data[,1]),type) {
  if (type == "a") {
    a <- ggplot2::ggplot(data=data) + geom_line(aes(x=data[,1],y=data[,2],color=data[,3])) + scale_x_continuous(limits =c(lower,upper)) + xlab("Calibrated date") + ylab("Probability density") + labs(color="Events")
    print(a)
  } else {
    if (type == "b") {
      b <- ggplot2::ggplot(data=data) + geom_area(position="identity", aes(x=data[,1],y=data[,2],fill=data[,3]), alpha=0.5) + scale_x_continuous(limits =c(lower,upper)) + xlab("Calibrated date") + ylab("Probability density") + labs(fill="Events")
      print(b)
    } else {
      if (type == "c") {
        a <- ggplot2::ggplot(data=data) + geom_line(aes(x=data[,1],y=data[,2],color=data[,3])) + scale_x_continuous(limits =c(lower,upper)) + xlab("Calibrated date")
        c <- a + geom_area(position="identity",aes(x=data[,1],y=data[,2],fill=data[,3]), alpha=0.5) + xlab("Calibrated date") + ylab("Probability density") + labs(fill="Events") + guides(color=FALSE)
        print(c)
      }
    }
  }
}