R/plot_Scatterplots.R

Defines functions ScatterSamples plot_Scatterplots

Documented in plot_Scatterplots ScatterSamples

#' @title Display Scatter Plot Matrix of the Bayesian Age Results
#'
#' @description Create a hexbin plot matrix ([hexbin::hexplom]) of age results returned by the Bayesian age calculation.
#'
#' @details
#'
#' **Additional supported plot arguments**\cr
#'
#' The following table lists additional arguments supported by the function in order to fine tune the
#' graphical output. Such arguments, can just be added in the function call. Example, for disabling
#' the [graphics::rug] in the plot mode `smoothScatter` you can type `plot_Scatterplots(..., rug = FALSE)`
#' Please note that not all arguments are supported by all plot types.
#'
#' \tabular{lll}{
#' **ARGUMENT** \tab **SUPPORTED BY PLOT TYPE** \tab **DESCRIPTION** \cr
#' `colramp` \tab `hexbin` and `smoothScatter` \tab Option to define an own colour ramp, by defining an own function, e.g., `function(n) heat.colors(n, alpha = 1)`. \cr
#' `pscales` \tab `hexbin` and `smoothScatter` \tab Controls the number of ticks shown on the plot axes, please note that the number works proportionally. \cr
#' `bw_smoothScatter` \tab `smoothScatter` \tab Controls the bandwidth of the smooth scatter, cf. [graphics::smoothScatter] \cr
#' `rug` \tab `smoothScatter` \tab enables/disables rugs \cr
#' `nlevels` \tab `smoothScatter` \tab controls the number of isolines shown (cf. [graphics::contour]) \cr
#' `nrpoints` \tab `smoothScatter` \tab defines the number of `nrpoints` to be plotted [graphics::smoothScatter] \cr
#' `col_contour` \tab `smoothScatter` \tab defines the colour of the contour lines \cr
#' `col_nrpoints` \tab `smoothScatter` \tab sets colour of the `nrpoints` in the scatter plot
#'
#' }
#'
#'
#' @param object [coda::mcmc.list] or a [data.frame] (**required**): mcmc list objects generated by [rjags::jags.model] in [AgeS_Computation],
#' [AgeC14_Computation] or [Age_OSLC14]. If a [data.frame] is provided, only the first two columns are taken and `NA` values are automatically
#' removed. The function can also handle `BayLum.list` objects directly for certain functions
#'
#' @param variables [character] (*with default*): variable to be selected for the scatter plot, e.g., `"A"`. Please note
#' that you can only select one variable at the time
#'
#' @param sample_names [character] (*optional*): sample names shown in the plot matrix
#'
#' @param sample_selection [numeric] (*with default*): vector of samples to be plotted in the scatter matrix, e.g.,
#' `c(1,2)` will plot the first two samples, `c(1,3)` will plot samples 1 and 3 and `c(1:3)` will plot the first
#' three samples
#'
#' @param n.chains [integer] (*with default*): allows to limit the number of chains shown,
#' by default the results of all chains are plotted.
#'
#' @param plot_type [character] (*with default*): switch between different plot types, `"hexbin"` (the default), based on
#' the function [hexbin::hexplom] and `smoothScatter` (the alternative) based on a highly customised plot function using the
#' function [graphics::smoothScatter]
#'
#' @param plot_mode [character] (*with default*): switch between a `matrix` plot mode and a `single` plot mode. The plot mode `single`
#' only works for `plot_type = smoothScatter` and creates a single plot panel for each sample. Please note that this cannot be further
#' combined with other par settings.
#'
#' @param ... further arguments to control the plot output, standard plot arguments supported are `main`, `xlab`, `ylab`, `xlim`, `ylim`, `cex`. For additional
#' arguments supporting a fine tuning of the plot, see details.
#'
#' @return
#' A scatter plot based on [hexbin::hexplom]
#'
#' @section Function version: 0.3.2
#'
#' @author Sebastian Kreutzer, Institute of Geography, Ruprecht-Karl University of Heidelberg (Germany) ,
#' based on the function 'ScatterSamples()' by Claire Christophe, Anne Philippe, Guillaume Guérin
#'
#' @seealso [Age_Computation], [AgeS_Computation], [AgeC14_Computation],
#' and [rjags] packages.
#'
#' @examples
#' data(AgeS,envir = environment())
#'
#' ##hexbin
#' plot_Scatterplots(
#'    object = AgeS$Sampling,
#'    sample_names = c("GDB5", "GDB3"),
#'    sample_selection = c(1,2)
#'  )
#'
#' ##scatter smooth (matrix)
#' plot_Scatterplots(
#'    object = AgeS$Sampling,
#'    sample_names = c("GDB5", "GDB3"),
#'    sample_selection = c(1,2),
#'    plot_type = "smoothScatter")
#'
#'
#' ##scatter smooth (single)
#' plot_Scatterplots(
#'    object = AgeS$Sampling,
#'    sample_names = c("GDB5", "GDB3"),
#'    sample_selection = c(1,2),
#'    plot_type = "smoothScatter",
#'    plot_mode = "single")
#'
#' @md
#' @export
plot_Scatterplots <- function(
  object,
  variables = c("A"),
  sample_names = NULL,
  sample_selection = NULL,
  n.chains = NULL,
  plot_type = "hexbin",
  plot_mode= "matrix",
  ...
){


  # Verify input --------------------------------------------------------------------------------
  if(inherits(object, "BayLum.list")) {
    if(!is.null(attributes(object)$originator)){
      ##select what to do for different functions
      switch(attributes(object)$originator,
             AgeS_Computation = object <- object$Sampling,
             Age_OSLC14 = object <- object$Sampling
      )
    }
  }

  if (is.null(attributes(object)$class) || attributes(object)$class != "mcmc.list")
    if(!inherits(object, "data.frame")){
      stop("[plot_Scatterplots()] Wrong input, only objects of type 'mcmc.list' or single 'data.frame' are allowed. Please check the manual!",
           call. = FALSE
      )

    }else{
      ##check whether the data.frame has two columns
      if(ncol(object) < 2)
        stop("[plot_Scatterplots()] Your 'data.frame' needs at least two columns!",
             call. = FALSE
        )

      ##select only the fist two, no experiments
      object <- object[,1:2]

      ##check for NA
      object <- na.exclude(object)

      ##check whether this is numeric
      if(any(!is.numeric(unlist(object)))){
        stop("[plot_Scatterplots()] Only numeric values are allowed!",
             call. = FALSE
        )
      }

      ##reset column names
      colnames(object) <- paste0("A[",1:2,"]")

      ##transform data.frame to mcmc.list
      object <- coda::as.mcmc.list(coda::as.mcmc(as.matrix(object)))

    }


  # Extract wanted parameters -------------------------------------------------------------------
  if(length(variables) > 1)
    stop("[plot_Scatterplots()] You can only select one variable at the time!", call. = FALSE)


  if(!all(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = "") %in% variables)){
    sel <- which(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = "") %in% variables)
    if(length(sel) == 0){
      allowed <- unique(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = ""))
      stop(paste0("[plot_Scatterplots()] Invalid 'variables', they did not match your dataset. Variable names of your dataset: ",
                  paste(allowed, collapse = ", "), "."), call. = FALSE)
    }

    ##create new object
    object <- as.mcmc.list(lapply(object, function(x){
      x[,sel, drop = FALSE]

    }))

  }

  # Preset values  ------------------------------------------------------------------------------
  ##get number of samples
  n.samples <- length(coda::varnames(object))

  ##make a sample selection
  if(is.null(sample_selection)){
    sample_selection <- 1:n.samples

  }else{
    if(length(sample_selection) > n.samples || length(sample_selection) < 2 || max(sample_selection) > n.samples || min(sample_selection) < 1){
      warning(paste0("[plot_Scatterplots()] You have only ", n.samples, " samples. 'sample_selection' wrong, reset to default!"), call. = FALSE)
      sample_selection <- 1:n.samples
    }

    ##reset n.samples
    n.samples <- length(sample_selection)

  }

  ##set sample names
  if(is.null(sample_names)){
    sample_names <- paste0("Sample ", 1:n.samples)

  }else{

    ##remove line breaks and trailing and leading space
    sample_names <- trimws(gsub(pattern = "\n", replacement = "", x = sample_names, fixed = TRUE))

    if(length(sample_names) < n.samples){
      warning("[plot_Scatterplots()] length of 'sample_names' shorter than the number of samples; default values used!", call. = FALSE)
      sample_names <- paste0("Sample ", 1:n.samples)

    }

    ##limit samples to selection, otherwise we risk to create a subscript out of bounds error
    sample_names <- sample_names[sample_selection]

  }

  ##get number of chains
  if(is.null(n.chains)){
    n.chains <- coda::nchain(object)

  }else{
    if(n.chains > coda::nchain(object) || n.chains < coda::nchain(object)){
      n.chains <- coda::nchain(object)
      warning(paste0("[plot_Scatterplots()] 'n.chains' setting wrong. You have ", n.chains, " chains, reset to default."), call. = FALSE)

    }

  }


  # Combine to matrix ---------------------------------------------------------------------------
  m_data <- do.call(rbind,object[1:n.chains])

    ##remove all columns we don't need, this makes the plot simple
    m_data <- m_data[, sample_selection]

    ##reset sample names
    colnames(m_data) <- sample_names[sample_selection]


  # # Plot output ---------------------------------------------------------------------------------
  ##make sure we do not screw up the par settings
  par.default <- par(no.readonly = TRUE)
    on.exit(do.call(
      what = par,
      args = list(
        mfrow = par.default$mfrow,
        cex = par.default$cex,
        fig = par.default$fig,
        mar = par.default$mar,
        omi = par.default$omi
      )
    ))

  ##set plot settings
  plot_settings <- list(
    xlab = "Age (ka)",
    ylab = "Age (ka)",
    colramp = if(plot_type == "hexbin") {
      function(n) terrain.colors(n, alpha = 1)
      } else {
        colorRampPalette(c("white", blues9))
      },
    pscales = 3,
    main = "Scatter Plots",
    bw_smoothScatter = 0.8,
    nlevels = 10,
    rug = TRUE,
    cex = 1.0,
    nrpoints = 100,
    col_nrpoints = "darkgray",
    col_contour = "black",
    xlim = NULL,
    ylim = NULL

  )

  ##reset list on demand
  plot_settings <- modifyList(x = plot_settings, val = list(...))

  ###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ###hexbin scatter plot
  ###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ##create plot
  if(plot_type == 'hexbin'){
    hexbin::hexplom(
      x = m_data,
      upper.panel = NULL,
      xlab = plot_settings$xlab,
      ylab = plot_settings$ylab,
      pscales = plot_settings$pscales,
      colramp = plot_settings$colramp,
      main = plot_settings$main
    )


  }else{
    ###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ###BayLum scatter plot function
    ###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    ###(0) extract datasets
    ##get possible combinations of the dataset
    cmb <- combn(colnames(m_data), m = 2)
    cmb <- cmb[2:1,ncol(cmb):1, drop = FALSE]

    ##create list of matricies with that combinations, and in the correct order, this is the
    ##tricky part
    data_list <- rev(unlist(lapply(colnames(m_data), function(x){
       temp <- cmb[,rev(which(cmb[2,] == x)), drop = FALSE]

       if(length(temp) > 0){
         temp <- lapply(1:ncol(temp), function(y){
            m_data[,temp[,y]]

         })
       }else{
         temp <- NULL
       }

       return(rev(temp))
    }), recursive = FALSE))

    ##rest list names
    names(data_list) <- vapply(data_list, function(x){
      paste(colnames(x), collapse = " <> ")
    }, character(1))

    ###(1) create plot matrix+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ###CREATE GRID MATRIX
      n <- ncol(m_data)
      name <- rep(sample_names,n)

      ##set prototype grid matrix with index
      m <- matrix(1:n^2, ncol = n, nrow = n, byrow = TRUE)

      ##reshuffle matrix, otherwise we get always the opposite
      m_re <- m[1:nrow(m),ncol(m):1]

      ##identify all diagonal members
      m_diag <- diag(m_re)

      ##identify upper triangle
      m_upper <- sort(m_re[upper.tri(m_re)])

      ##identify lower triangle
      m_lower <- sort(m_re[lower.tri(m_re)])

      ##get subdiagonal members, only this will have axis labelling
      m_diag_sub <- m_diag[2:length(m_diag)] + 1


    ###(2) - SET PROTOYPE FUNCTIONS+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ##PLOT 1 - EMPTY PLOT
    empty_plot <- function(){
        plot(NA,NA,xlim = c(0,1), ylim = c(0,1), xlab = "", ylab = "", xaxt = "n", yaxt = "n")

    }

    ##PLOT 2 - NAME PLOT +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    name_plot <-
      function(data,
               name,
               xdens = TRUE,
               ydens = TRUE,
               frame.plot = TRUE,
               ydens_opposite = FALSE,
               cex_scale = 1,
               xlim = plot_settings$xlim,
               ylim = plot_settings$ylim,
               main = if(plot_mode == "matrix") "" else name
               ) {

      ##define transfer function
      transfer <- function(x, min_scale, max_scale){
        m <- (max(x) - min(x)) / (max_scale - min_scale)
        n <- min(x)  - m * min_scale
        return((x - n)/m)
      }

      ##rset axes to make sure that the density plot is shown with the correct scale
      if(is.null(xlim) || is.null(ylim)){
        con <- KernSmooth::bkde2D(x = data[,1:2], bandwidth = plot_settings$bw_smoothScatter)

        ##xlim
        if(is.null(xlim))
          xlim <- range(con$x1)

        ##ylim
        if(is.null(ylim))
          ylim <- range(con$x2)

      }

      ##create plot area
      plot(
        NA,
        NA,
        xlim = xlim,
        ylim = ylim,
        xlab = "",
        ylab = "",
        xaxt = "n",
        yaxt = "n",
        main = main,
        frame.plot = frame.plot
      )

        ##calculate density
        density1 <- density(data[,1])
        density2 <- density(data[,2])

        ##draw density polygones
        ##marginal density on the xaxis
        if(xdens){
          x <- density1$x
          y <- transfer(density1$y, min_scale = ylim[1], max_scale = ylim[2])
          y <- ((y - ylim[1]) * 0.25 * cex_scale )+ ylim[1]
          polygon(
            x = c(x, rev(x)),
            y = c(y, rep(ylim[1], length(x))) - (ylim[1] - par()$usr[3]),
            col = rgb(0,0,0,0.1),
            border = "darkgray"
            )

        }

        ##marginal density y-axis
        if(ydens){
          if(ydens_opposite == FALSE){
             x <- transfer(-density2$y, min_scale = xlim[1], max_scale = xlim[2])
             x <- -(((-x + max(x)) * 0.25 * cex_scale) - max(x))
             y <- density2$x

             polygon(
               x = c(x, rep(xlim[2], length(x))) + (par()$usr[2] - xlim[2]),
               y = c(y, rev(y)),
               col = rgb(0,0,0,0.1),
               border = "darkgray"
             )

          }else{
              x <- transfer(density2$y, min_scale = xlim[1], max_scale = xlim[2])
              x <- ((x - xlim[1]) * 0.25 * cex_scale) + xlim[1]
              y <- density2$x

              polygon(
                x = c(x, rep(xlim[1], length(x))) - (xlim[1] - par()$usr[1]),
                y = c(y, rev(y)),
                col = rgb(0,0,0,0.1),
                border = "darkgray"
            )

          }

        }

        ##add sample name
        if(!is.null(names) && (is.null(main) || main == ""))
          text(x = mean(xlim), y = mean(ylim), labels = name, cex = 1.2 * cex_scale)

      }

      ##PLOT 3 - SCATTER PLOT
      scatter_plot <- function(
        data,
        xaxt = TRUE,
        yaxt = TRUE,
        xrug = TRUE,
        yrug = TRUE,
        axes_sides = c(x = 3, y = 2),
        xlab = plot_settings$xlab,
        ylab = plot_settings$ylab,
        xlim = plot_settings$xlim,
        ylim = plot_settings$ylim
        ){

        ##get information from the kernel smoother to crate
        con <- KernSmooth::bkde2D(x = data[,1:2], bandwidth = plot_settings$bw_smoothScatter)

        ##reset axes to make sure that the scaling is good enough
        ##xlim
        if(is.null(xlim))
          xlim <- range(con$x1)

        ##ylim
        if(is.null(ylim))
          ylim <- range(con$x2)

        ##add scatter
        smoothScatter(
          x = data[, 1],
          y = data[, 2],
          bandwidth = plot_settings$bw_smoothScatter,
          xlab = if(axes_sides[["x"]] == 1) xlab else "",
          ylab = if(axes_sides[["x"]] == 1) ylab else "",
          xaxt = "n",
          yaxt = "n",
          nrpoints = plot_settings$nrpoints,
          col = plot_settings$col_nrpoints,
          colramp = plot_settings$colramp,
          xlim = xlim,
          ylim = ylim
        )


        ##add contour lines
        contour(
          y = con$x2,
          x = con$x1,
          z = con$fhat,
          add = TRUE,
          labcex = 0.75 * plot_settings$cex,
          col = plot_settings$col_contour,
          nlevels = plot_settings$nlevels,
          xlim = xlim,
          ylim = ylim
        )

        ##add rug
        if(xrug)
          rug(x = data[,1], side = axes_sides[["x"]], col = rgb(0,0,0,0.3))

        if(yrug)
          rug(x = data[,2], side = axes_sides[["y"]], col = rgb(0,0,0,0.3))


        ##add x-axis above
        if(xaxt){
          at <- pretty(data[,1], n = plot_settings$pscales)
          at <- at[-c(1,length(at))]
          axis(side = axes_sides[["x"]], at = at, labels = NULL)

        }

        ##add y-axis above
        if(yaxt){
          at <- pretty(data[,2], n = plot_settings$pscales)
          at <- at[-c(1,length(at))]
          axis(side = axes_sides[["y"]], at = at, labels = NULL)

        }

      }

      ###PLOT
      ##set par
      if(plot_mode == "matrix")
        par(mfrow = c(n,n), mar = c(0,0,0,0), oma = c(4,4,3,3), cex = plot_settings$cex)

      ##now we have to create all the plots and then decide what
      ##needs to be plotted when
      ##set plot counter
      p <- 1

      ##set number of windows to plot, these differs for the single mode
      if(plot_mode == "matrix"){
        windows <- 1:n^2

      }else{
        windows <- m_lower

      }

      ##START LOOP
      for(i in windows){
        ##(A) plot empty plots
        if(i %in% m_upper && plot_mode == "matrix")
          empty_plot()

        ##(B) plot diagonale plots
        if(i %in% m_diag && plot_mode == "matrix"){

          ##do not plot density for the first and the last, otherwise it looks odd
          if(i %in% m_diag[1]){
            name_plot(data = data_list[[p]][,c(1,1)], name = name[i], ydens = FALSE)

          }else if(i %in% m_diag[length(m_diag)]){
            name_plot(data = data_list[[p]][,c(2,2)], name = name[i], xdens = FALSE)

          }else{
           name_plot(
              data = matrix(rep(m_data[,name[i]],2), ncol = 2), name = name[i])

          }


        }

        ##(C) plot scatter plot
        ## (here we have to decide whether the axis is plotted or not)
        if(i %in% m_lower){
          if(plot_mode != "matrix" || i %in% m_diag_sub){
           if(plot_settings$rug){
            xaxt = TRUE
            yaxt = TRUE
            xrug = TRUE
            yrug = TRUE

           }else{
             xaxt = FALSE
             yaxt = FALSE
             xrug = FALSE
             yrug = FALSE

           }

          }else{
            xaxt = FALSE
            yaxt = FALSE
            xrug = FALSE
            yrug = FALSE

          }

          ##allow two plot modes
          if(plot_mode == "matrix"){
            scatter_plot(
              data = data_list[[p]],
              xaxt = xaxt,
              yaxt = yaxt,
              xrug = xrug,
              yrug = yrug
            )

          }else{
            par(fig = c(0,0.9,0,0.9), mar = c(4,4,0,0), omi = c(0.42,.42,.42,.42), cex = plot_settings$cex)
            scatter_plot(
              data = data_list[[p]],
              xaxt = xaxt,
              yaxt = yaxt,
              xrug = xrug,
              yrug = yrug,
              axes_side = c(x = 1, y = 2)
            )

            par(fig = c(0.9,1,0,.9),  mar = c(4,0,0,0.7), new = TRUE)
            name_plot(
              data = data_list[[p]],
              name = NULL,
              xdens = FALSE,
              frame.plot = FALSE,
              ydens_opposite = TRUE,
              cex_scale = 3 * plot_settings$cex
            )

            par(fig = c(0,0.9,0.9,1), mar = c(0,4,0.7,0), new = TRUE)
            name_plot(
              data = data_list[[p]],
              name = names(data_list)[p],
              xdens = TRUE,
              ydens = FALSE,
              frame.plot = FALSE,
              cex_scale = 3 * plot_settings$cex)

          }


          ##update counter
          p <- p + 1

        }

      } #end loop for plotting

      ##ADD MTEXT
      ##add axis labelling
      if(plot_mode == "matrix"){
      mtext(plot_settings$xlab, side = 1, outer = TRUE, line = 1,cex = 0.95 * plot_settings$cex)
      mtext(plot_settings$ylab, side = 2, outer = TRUE, line = 1,cex = 0.95 * plot_settings$cex)

      ##add main
      mtext(plot_settings$main, side = 3, outer = TRUE, line = 1, cex = 1.1 * plot_settings$cex)
      }

  }

}

#TODO
#'Old function ScatterPlots()
#'@rdname plot_Scatterplots
#'@md
#'@export
ScatterSamples <- function(...){
  .Defunct(
    new = "plot_Scatterplots()",
    package = "BayLum"
  )

}
R-Lum/BayLum documentation built on March 19, 2024, 12:42 p.m.