R/MarginalPlot.R

Defines functions marginal_plot MarginalPlot

Documented in marginal_plot MarginalPlot

#####################################################
#          Marginal posterior Density               #
#       NEW version in ArchaeoPhases 1.4            #
#####################################################
#' Plot a marginal posterior density
#'
#' Draws a plot of the estimated marginal posterior density for the one-parameter and adds the mean
#' and the credible interval at the desired level
#'
#' @param a_chain Numeric vector containing the output of the MCMC algorithm
#' for the parameter.
#' @param level Probability corresponding to the level of confidence.
#' @param GridLength Length of the grid used to estimate the density.
#' @param title Title of the graph.
#' @param subtitle Subtitle of the graph.
#' @param caption Caption of the graph.
#' @param x.label Label of the x-axis.
#' @param y.label Label of the y-axis.
#' @param width Plot width in \code{units}.
#' @param height Plot height in \code{units}.
#' @param units String recognized by the \code{ggsave()} function,
#' one of "in", "cm", "mm".
#' @param x.min Minimum x axis value.
#' @param x.max Maximum x axis value.
#' @param x.scale One of "calendar" for calendar years,
#' "BP" for years before present,
#' or "elapsed" for time elapsed from a specified origin.
#' @param elapsed.origin.position Position of the column to use as the
#' origin for elapsed time calculations.
#' @param y.grid Switch for horizontal grid lines.
#' @param file Name of the file that will be saved if chosen, default = \code{NULL}.
#' @param newWindow Whether or not the plot is drawn within a new window.
#'
#' @return \code{NULL}, called for its side effects
#'
#' @details The density is estimated using \code{density()} function with \code{n = GridLength}.
#'
#' @author Anne Philippe, \email{Anne.Philippe@@univ-nantes.fr} and
#'
#' @author  Marie-Anne Vibet, \email{Marie-Anne.Vibet@@univ-nantes.fr}
#'
#' @examples
#'   data(Events);
#'   MarginalPlot(a_chain = Events$Event.1, level = 0.95)
#'
#' @importFrom stats density
#' @importFrom grDevices dev.new
#'
#' @export
#'
MarginalPlot <- function(a_chain, level=0.95, GridLength=1024,
                  title="Characteristics of a date", subtitle = NULL,caption = "ArchaeoPhases",
                  x.label = "Calendar year",y.label = NULL,y.grid = TRUE,
                  x.scale = "calendar", elapsed.origin.position = NULL,x.min = NULL, x.max = NULL,
                  height = 7, width = 7, units = "in",file = NULL, newWindow=TRUE){
  # x.scale can either be "calendar", "BP" or "elapsed" if any other origin that 0 and 1950
  if (x.scale == "BP") {
    a_chain <- 1950-a_chain
  }
  if (x.scale == "elapsed") {
    a_chain <- elapsed.origin.position-a_chain
  }
  # credible interval
  CR <- CredibleInterval(a_chain, level=level, roundingOfValue=4)
  CR <- c(CR,"y"=0,"yend"=0)
  # Mean
  Mean = mean(a_chain)
  # new dataframe
  data = data.frame("a_chain" = as.vector(a_chain))
  dataCR = data.frame("Inf" = CR[2], "Sup" = CR[3], "y"=0,"yend"=0, "Mean"=Mean)

  h <- ggplot2::ggplot(data = data, ggplot2::aes(x=a_chain))
  h <- h + ggplot2::geom_density(n = GridLength)
  h <- h + ggplot2::geom_segment(data = dataCR, ggplot2::aes(x=dataCR[1,1], xend=dataCR[1,2], y = dataCR[1,3], yend = dataCR[1,4], colour="steelblue"), size=3.5, show.legend=F)
  h <- h + ggplot2::geom_point(data = dataCR, ggplot2::aes(x=Mean, y = dataCR[1,3]), size = 2)
  #h <- h + ggplot2::scale_color_manual(values =c('#56B4E9', '#FC4EO7'))
  h <- h + ggplot2::labs(x = x.label, y = y.label, title = title, subtitle = subtitle, caption = caption)

  if (y.grid==FALSE) {
    h <- h + ggplot2::theme(panel.grid.major.y=ggplot2::element_blank())
  }
  # x abscisses
  if (is.null(x.min) ) {
    x.min <- min(density(a_chain, n=GridLength)$x)
  }
  if (is.null(x.max)) {
    x.max <- max(density(a_chain, n=GridLength)$x)
  }
  h <- h + ggplot2::xlim(x.min, x.max)

  # export file
  if (!is.null(file)) {
    ggplot2::ggsave(filename = file, plot = h, height = height,width = width, units = units)
  }
  if(newWindow == TRUE) {
    dev.new(height = height, width = width)
  }
  print(h)

}

#' Plot a marginal posterior density
#'
#' Draws a plot of the marginal posterior density for a single parameter, with
#' an option to add the mean and the credible interval at the desired level
#'
#' @details
#' The plot is drawn with the current theme and color scales; the function
#' does not alter or override theme elements.
#'
#' @param data Data frame containing the output of the MCMC algorithm.
#' @param position Index of the column corresponding to the MCMC chain of
#' interest, or a column name.
#' @param level Probability corresponding to the level of confidence.
#' @param grid_length Length of the grid used to estimate the density.
#' @param title Title of the graph.  The default uses the \code{data}
#' column name.
#' @param subtitle Subtitle of the graph.  The default is
#' "Marginal posterior density".
#' @param caption Caption of the graph.  The default describes the
#' confidence of the credible interval.
#' @param x_label Label of the x-axis.
#' @param y_label Label of the y-axis.
#' @param width Plot width in \code{units}.
#' @param height Plot height in \code{units}.
#' @param units String recognized by the \code{ggsave()} function,
#' one of "in", "cm", "mm".  This parameter has no effect on the
#' display plot.
#' @param x_min Minimum x axis value.
#' @param x_max Maximum x axis value.
#' @param x_scale One of "calendar" for calendar years,
#' "BP" for years before present,
#' or "elapsed" for time elapsed from a specified origin.
#' @param elapsed_origin_position Position of the column to use as the
#' origin for elapsed time calculations.
#' @param y_grid Switch for horizontal grid lines.
#' @param file Name of the file that will be saved if chosen,
#' default = \code{NULL}.
## ' @param new_window Whether or not the plot is drawn within a new window.
#' @param plot_result If \code{TRUE}, then draw a plot on the display,
#' else suppress drawing.
#' @param mean_linetype The \code{linetype} used to indicate the mean density.
#' @param mean_color The color of the line used to indicate mean density.
#' @param mean_size The width of the line used to indicate the mean density.
#' @param ci_linetype The \code{linetype} used to indicate the credible
#' intervals.
#' @param ci_color The color of the lines used to indicate the credible
#' intervals.
#' @param ci_size The width of the lines used to indicate the credible
#' intervals.
#' @param line_linetype The \code{linetype} used to indicate the density.
#' @param line_color The color of the line used to indicate the density.
#' @param line_size The width of the line used to indicate the density.
#' @param density_color Color to use if fill_palette is not specified.
#' @param fill_palette Palette to use for fills.
#'
#' @return An \code{archaeophases_plot} object with the data and metadata
#' needed to reproduce the plot.
#'
#' @author Anne Philippe, \email{Anne.Philippe@@univ-nantes.fr};
#' @author Marie-Anne Vibet, \email{Marie-Anne.Vibet@@univ-nantes.fr}; and
#' @author Thomas S. Dye, \email{tsd@tsdye.online}
#'
#' @examples
#'   data(Events)
#'   mp <- marginal_plot(data = Events, position = 2, level = 0.95)
#'   ## View data and metadata
#'   str(mp)
#'
#' @importFrom stats density approx
#' @importFrom grDevices dev.new
#' @importFrom ggplot2 ggplot aes geom_ribbon geom_segment geom_line labs theme xlim ggsave
#' @importFrom magrittr %$%
#' @export
#'
marginal_plot <- function(data,
                          position = 1,
                          level = 0.95,
                          grid_length = 1024,
                          title = if (is.numeric(position)) names(data)[position] else position,
                          subtitle = "Marginal posterior density",
                          caption = paste(level * 100, "% credible interval", sep=""),
                          x_label = "Calendar year",
                          y_label = "Density",
                          y_grid = TRUE,
                          x_scale = "calendar",
                          elapsed_origin_position = NULL,
                          x_min = NULL, x_max = NULL,
                          height = 7, width = 7,
                          units = "in",
                          file = NULL,
                          ## new_window = TRUE,
                          plot_result = TRUE,
                          mean_linetype = "dashed",
                          mean_color = "white",
                          mean_size = 0.5,
                          ci_linetype = "dotted",
                          ci_color = mean_color,
                          ci_size = mean_size,
                          line_linetype = "solid",
                          line_color = "black",
                          line_size = 1,
                          density_color = "gray30",
                          fill_palette = NULL)
{

    if (!is.data.frame(data)) stop("Data format not recognized.")

    if (length(position) > 1) stop("Expected one marginal posterior.")

    ## a_chain variable is a vector
    if (is.element("archaeophases_plot", class(data))) {
        a_chain <- as.vector(unlist(data[, 1]))
    }
    else {
        a_chain <- as.vector(unlist(data[, position]))

        if (x_scale == "BP") {
            a_chain <- 1950 - a_chain
        }

        if (x_scale == "elapsed") {
            elapsed_origin <- data[, elapsed_origin_position]
            a_chain <- a_chain - elapsed_origin
        }
    }

    ## credible interval
    c_i <- credible_interval(a_chain, level = level, round_to = 4)$ci

    ## mean
    chain_mean = mean(a_chain)

    chain_density <- density(a_chain, n = grid_length) %$%
        data.frame(x = x, y = y) %>%
        mutate(mid = (x > c_i["inf"] & x < c_i["sup"]))

    dens_mean <- approx(chain_density, xout = chain_mean)
    dens_low <- approx(chain_density, xout = c_i["inf"])
    dens_high <- approx(chain_density, xout = c_i["sup"])

    if (!is.null(fill_palette)) {
        fill_color <- unname(fill_palette[density_color])
    }
    else {
        fill_color <- density_color
    }

    a_chain_df <- as.data.frame(a_chain)
    names(a_chain_df) <- "x"

    h <- ggplot(data = a_chain_df,
                mapping = aes(x = x))

    h <- h + geom_density(mapping = aes(y = ..density..),
                          color = line_color,
                          linetype = line_linetype,
                          size = line_size,
                          fill = fill_color,
                          n = grid_length)

    h <- h + geom_segment(mapping = aes(x = c_i[["inf"]],
                                        xend = c_i[["inf"]],
                                        y = 0,
                                        yend = dens_low[["y"]]),
                          linetype = ci_linetype,
                          color = ci_color,
                          size = ci_size)

    h <- h + geom_segment(mapping = aes(x = c_i[["sup"]],
                                        xend = c_i[["sup"]],
                                        y = 0,
                                        yend = dens_high[["y"]]),
                          linetype = ci_linetype,
                          color = ci_color,
                          size = ci_size)

    h <- h + geom_segment(mapping = aes(x = chain_mean,
                                        xend = chain_mean,
                                        y = 0,
                                        yend = dens_mean[["y"]]),
                          linetype = mean_linetype,
                          color = mean_color,
                          size = mean_size)

    h <- h + labs(x = x_label, y = y_label,
                  title = title, subtitle = subtitle,
                  caption = caption)

    if (y_grid==FALSE) {
        h <- h + theme(panel.grid.major.y = ggplot2::element_blank())
    }

    ## x abscisses
    if (is.null(x_min) ) {
        x_min <- min(density(unlist(a_chain), n = grid_length)$x)
    }

    if (is.null(x_max)) {
        x_max <- max(density(unlist(a_chain), n = grid_length)$x)
    }

    h <- h + xlim(x_min, x_max)

    ## export file
    if (!is.null(file)) {
        ggsave(filename = file, plot = h, height = height,
               width = width, units = units)
    }

    if (plot_result == TRUE) {
        ## if(new_window == TRUE) {
        ##     dev.new(height = height, width = width)
        ## }
        print(h)
    }
    new_archaeophases_plot(x = a_chain_df,
                           mcmc = data,
                           call = match.call())
}

Try the ArchaeoPhases package in your browser

Any scripts or data that you put into this service are public.

ArchaeoPhases documentation built on June 22, 2022, 1:05 a.m.