R/plot_Histogram.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Plot a histogram with separate error plot
#'
#' @description
#' Function plots a predefined histogram with an accompanying error plot as
#' suggested by Rex Galbraith at the UK LED in Oxford 2010.
#'
#' If the normal curve is added, the y-axis in the histogram will show the
#' probability density.
#'
#' A statistic summary, i.e. a collection of statistic measures of
#' centrality and dispersion (and further measures) can be added by specifying
#' one or more of the following keywords:
#' - `"n"` (number of samples),
#' - `"mean"` (mean De value),
#' - `"mean.weighted"` (error-weighted mean),
#' - `"median"` (median of the De values),
#' - `"median.weighted"` (error-weighted median),
#' - `"sd.rel"` (relative standard deviation in percent),
#' - `"sd.rel.weighted"` (error-weighted relative standard deviation in percent),
#' - `"sd.abs"` (absolute standard deviation),
#' - `"sd.abs.weighted"` (error-weighted absolute standard deviation),
#' - `"se.rel"` (relative standard error in percent),
#' - `"se.rel.weighted"` (error-weighted relative standard error in percent),
#' - `"se.abs"` (absolute standard error),
#' - `"se.abs.weighted"` (error-weighted absolute standard error),
#' - `"kurtosis"` (kurtosis) and
#' - `"skewness"` (skewness).
#'
#' @param data [data.frame] or [Luminescence::RLum.Results-class] object (**required**):
#' for `data.frame`: two columns: De (`data[,1]`) and De error (`data[,2]`).
#' If the error column is missing or only contains `NA` values, then `se` is
#' set to `FALSE`, and the weighted statistics will match the unweighted ones.
#' In any other case, rows with `NA` values will be removed prior to plotting.
#'
#' @param mtext [character] (*optional*):
#' further sample information ([mtext]).
#'
#' @param cex.global [numeric] (*with default*):
#' global scaling factor.
#'
#' @param se [logical] (*with default*):
#' plots standard error points over the histogram, default is `TRUE`. It is
#' silently reset to `FALSE` if the error column in `data` is missing or only
#' contains `NA` values.
#'
#' @param rug [logical] (*with default*):
#' adds rugs to the histogram, default is `TRUE`.
#'
#' @param normal_curve [logical] (*with default*):
#' adds a normal curve to the histogram. Mean and standard deviation are
#' calculated from the input data. If `TRUE`, the y-axis in the histogram
#' will show the probability density.
#'
#' @param summary [character] (*with default*):
#' add statistic measures of centrality and dispersion to the plot.
#' Can be one or more of several keywords. See details for available keywords.
#'
#' @param summary.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`)
#' for the statistical summary. Alternatively, the keyword `"sub"` may be
#' specified to place the summary below the plot header. However, this latter
#' option in only possible if `mtext` is not used. In case of coordinate
#' specification, y-coordinate refers to the right y-axis.
#'
#' @param colour [numeric] or [character] (*with default*):
#' optional vector of length 4 which specifies the colours of the following
#' plot items in exactly this order: histogram bars, rug lines and summary
#' text, normal distribution curve, standard error points
#' (e.g., `c("grey", "black", "red", "grey")`).
#'
#' @param interactive [logical] (*with default*):
#' create an interactive histogram plot (requires the 'plotly' package)
#'
#' @param ... further arguments and graphical parameters to control the plot
#' output (see [plot] and [hist]). Supported are: `main`, `xlab`, `ylab`,
#' `xlim`, `ylim`, `breaks`, `pch`, `pt.cex` (point size), and `fun`. If y-axis labels are
#' provided, these must be specified as a vector of length 2 since the plot
#' features two axes (e.g. `ylab = c("axis label 1", "axis label 2")`). Y-axes
#' limits (`ylim`) must be provided as vector of length four, with the first two
#' elements specifying the left axes limits and the latter two elements giving
#' the right axis limits.
#'
#' @section Function version: 0.4.6
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)
#'
#' @seealso [hist], [plot]
#'
#' @examples
#'
#' ## load data
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <- convert_Second2Gray(ExampleData.DeValues$BT998,
#'                                             dose.rate = c(0.0438,0.0019))
#'
#' ## plot histogram the easiest way
#' plot_Histogram(ExampleData.DeValues)
#'
#' ## plot histogram with some more modifications
#' plot_Histogram(ExampleData.DeValues,
#'                rug = TRUE,
#'                normal_curve = TRUE,
#'                cex.global = 0.9,
#'                pch = 2,
#'                colour = c("grey", "black", "blue", "green"),
#'                summary = c("n", "mean", "sd.rel"),
#'                summary.pos = "topleft",
#'                main = "Histogram of De-values",
#'                mtext = "Example data set",
#'                ylab = c(expression(paste(D[e], " distribution")),
#'                         "Standard error"),
#'                xlim = c(100, 250),
#'                ylim = c(0, 0.1, 5, 20))
#'
#' @md
#' @export
plot_Histogram <- function(
  data,
  mtext = "",
  cex.global = 1,
  se = TRUE,
  rug = TRUE,
  normal_curve = FALSE,
  summary = "",
  summary.pos = "sub",
  colour = c("white", "black", "red", "black"),
  interactive = FALSE,
  ...
) {
  .set_function_name("plot_Histogram")
  on.exit(.unset_function_name(), add = TRUE)

  ## Integrity checks -------------------------------------------------------

  .validate_class(data, c("data.frame", "RLum.Results"))
  .validate_not_empty(data)
  if (inherits(data, "RLum.Results")) {
    data <- get_RLum(data)
  }

  ## we don't check the second column, as that can be NA
  .validate_class(data[, 1], c("numeric", "integer"),
                  name = "All columns of 'data'")

  .validate_class(mtext, "character", length = 1)
  .validate_class(cex.global, "numeric")
  .validate_class(summary, "character")
  summary.pos <- .validate_position(summary.pos, sub = TRUE)
  .validate_logical_scalar(se)
  .validate_logical_scalar(rug)
  .validate_logical_scalar(normal_curve)
  .validate_class(colour, "character", length = 4)

  ## handle error-free data sets
  if (length(data) < 2 || all(is.na(data[, 2]))) {
    data[, 2] <- 1
    se <- FALSE
  }

  ## Set general parameters ---------------------------------------------------

  ## read out additional arguments list
  extraArgs <- list(...)

  ## count and remove NA values
  n.NA <- sum(is.na(data[,1]))
  if (n.NA > 0) {
      .throw_message(sprintf("%d NA value%s excluded\n",
                             n.NA, ifelse(n.NA > 1, "s", "")),
                     error = FALSE)
      data <- data[!is.na(data[, 1]), ]
  }

  if("breaks" %in% names(extraArgs)) {
    breaks.plot <- extraArgs$breaks
    breaks_calc <- hist(x = data[,1],
                        breaks = breaks.plot,
                        plot = FALSE)$breaks
  } else {
    breaks.plot <- hist(x = data[,1],
                        plot = FALSE)$breaks
    breaks_calc <- breaks.plot
  }

  if("ylim" %in% names(extraArgs)) {
    ylim.plot <- extraArgs$ylim
    .validate_class(ylim.plot, "numeric", length = 4, name = "'ylim'")
  } else {
    H.lim <- hist(data[,1],
                  breaks = breaks.plot,
                  plot = FALSE)
    left.ylim <- c(0, if (normal_curve) max(H.lim$density) else max(H.lim$counts))
    range.error <- try(expr = range(data[,2], na.rm = TRUE),
                       silent = TRUE)
    range.error[1] <- ifelse(is.infinite(range.error[1]), 0, range.error[1])
    range.error[2] <- ifelse(is.infinite(range.error[2]), 0, range.error[2])
    ylim.plot <- c(left.ylim, range.error)
  }

  settings <- modifyList(list(
      main = "Histogram",
      xlab = expression(paste(D[e], " [Gy]")),
      ylab = c("Frequency", "Standard error"),
      xlim = range(breaks_calc),
      pch = 1,
      pt.cex = 1,
      fun = FALSE
  ), extraArgs)

  ## Set plot area format
  par.default <- .par_defaults()
  on.exit(par(par.default), add = TRUE)
  par(mar = c(4.5, 4.5, 4.5, if (se) 4.5 else 1),
      cex = cex.global)

  ## Plot histogram -----------------------------------------------------------
  HIST <- hist(data[,1],
               main = settings$main,
               xlab = settings$xlab,
               ylab = settings$ylab[1],
               xlim = settings$xlim,
               ylim = settings$ylim[1:2],
               breaks = breaks.plot,
               freq = !normal_curve,
               col = colour[1]
  )

  ## Optionally, add rug ------------------------------------------------------
  if (rug) {
    graphics::rug(data[, 1], col = colour[2], quiet = TRUE)
  }

  ## Optionally, add a normal curve based on the data -------------------------
  if (normal_curve) {
    ## silence notes raised by R CMD check
    x <- NULL

    ## add normal distribution curve
    curve(dnorm(x,
                mean = mean(na.exclude(data[,1])),
                sd = sd(na.exclude(data[,1]))),
          col = colour[3],
          add = TRUE,
          lwd = 1.2 * cex.global)
  }

  ## calculate and paste statistical summary
  data.stats <- list(data = data)

  ## calculate and paste statistical summary
  De.stats <- matrix(nrow = length(data), ncol = 18)
  colnames(De.stats) <- c("n",
                          "mean",
                          "mean.weighted",
                          "median",
                          "median.weighted",
                          "kde.max",
                          "sd.abs",
                          "sd.rel",
                          "se.abs",
                          "se.rel",
                          "q25",
                          "q75",
                          "skewness",
                          "kurtosis",
                          "sd.abs.weighted",
                          "sd.rel.weighted",
                          "se.abs.weighted",
                          "se.rel.weighted")

  ## placeholder for the summary label text
  label.text <- list()
  is.sub <- summary.pos[1] == "sub"

  for(i in 1:length(data)) {
    statistics <- calc_Statistics(data)
    De.stats[i,1] <- statistics$weighted$n
    De.stats[i,2] <- statistics$unweighted$mean
    De.stats[i,3] <- statistics$weighted$mean
    De.stats[i,4] <- statistics$unweighted$median
    De.stats[i,5] <- statistics$weighted$median
    De.stats[i,7] <- statistics$unweighted$sd.abs
    De.stats[i,8] <- statistics$unweighted$sd.rel
    De.stats[i,9] <- statistics$unweighted$se.abs
    De.stats[i,10] <- statistics$weighted$se.rel
    De.stats[i,11] <- quantile(data[,1], 0.25)
    De.stats[i,12] <- quantile(data[,1], 0.75)
    De.stats[i,13] <- statistics$unweighted$skewness
    De.stats[i,14] <- statistics$unweighted$kurtosis
    De.stats[i,15] <- statistics$weighted$sd.abs
    De.stats[i,16] <- statistics$weighted$sd.rel
    De.stats[i,17] <- statistics$weighted$se.abs
    De.stats[i,18] <- statistics$weighted$se.rel

    ##kdemax - here a little doubled as it appears below again
    De.density <- NA
    De.stats[i,6] <- NA
    if(nrow(data) >= 2){
      De.density <-density(x = data[,1],
                           kernel = "gaussian",
                           from = settings$xlim[1],
                           to = settings$xlim[2])

      De.stats[i,6] <- De.density$x[which.max(De.density$y)]
    }

    ## convert to a list of lists, like the object produced by calc_Statistics
    De.stats.list <- list(as.list(De.stats[i, ]))
    names(De.stats.list) <- "unweighted" # dummy placeholder

    ## generate the summary label text
    label.text[[i]] <- .create_StatisticalSummaryText(
        De.stats.list,
        keywords = summary,
        sep = ifelse(is.sub, " | ", "\n")
    )
  }

  ## convert keywords into summary placement coordinates
  coords <- .get_keyword_coordinates(summary.pos, settings$xlim, ylim.plot[1:2])
  summary.pos <- coords$pos
  summary.adj <- coords$adj

  ## add summary content
  for(i in 1:length(data.stats)) {
    if(summary.pos[1] != "sub") {
      text(x = summary.pos[1],
           y = summary.pos[2],
           adj = summary.adj,
           labels = label.text[[i]],
           col = colour[2],
           cex = 0.8)
    } else if (mtext == "") {
        mtext(side = 3,
              line = 1 - i,
              text = label.text[[i]],
              col = colour[2],
              cex = cex.global * 0.8)
    }
  }

  ## Optionally, add standard error plot --------------------------------------

  if (se) {
    par(new = TRUE)
    plot.data <- data[!is.na(data[,2]),]

    plot(x = plot.data[,1],
         y = plot.data[,2],
         xlim = settings$xlim,
         ylim = ylim.plot[3:4],
         pch = settings$pch,
         col = colour[4],
         cex = settings$pt.cex,
         main = "",
         xlab = "",
         ylab = "",
         axes = FALSE,
         frame.plot = FALSE
    )
    axis(side = 4,
         labels = TRUE,
         cex = cex.global
    )
    mtext(settings$ylab[2],
          side = 4,
          line = 3,
          cex = cex.global)
  }

  ## Optionally add user-defined mtext
  mtext(side = 3,
        line = 0.5,
        text = mtext,
        cex = 0.8 * cex.global)

  ## FUN by R Luminescence Team
  if (settings$fun && !interactive) sTeve() # nocov

  ## Optionally: Interactive Plot ----------------------------------------------
  if (interactive) {
    .require_suggested_package("plotly", "The interactive histogram")

    ## tidy data ----
    data <- as.data.frame(data)
    colnames(data) <- c("x", "y")
    x <- y <- NULL # suffice CRAN check for no visible binding
    if (any(grepl("paste", as.character(settings$xlab), fixed = TRUE)))
      settings$xlab <- "Equivalent dose [Gy]"

    ## create plots ----

    # histogram
    hist <- plotly::plot_ly(data = data, x = ~x,
                            type = "histogram",
                            showlegend = FALSE,
                            name = "Bin", opacity = 0.75,
                            marker = list(color = "#428BCA",
                                          line = list(width = 1.0,
                                                      color = "white")),
                            histnorm = ifelse(normal_curve, "probability density", ""),
                            yaxis = "y"
    )

    # normal curve ----
    if (normal_curve) {
      density.curve <- density(data$x)
      normal.curve <- data.frame(x = density.curve$x, y = density.curve$y)
      hist <- plotly::add_trace(hist, data = normal.curve, x = ~x, y = ~y,
                                inherit = FALSE,
                                type = "scatter", mode = "lines",
                                line = list(color = "red"),
                                name = "Normal curve",
                                yaxis = "y")
    }

    # scatter plot of individual errors
    if (se) {
      yaxis2 <- list(overlaying = "y", side = "right",
                     showgrid = FALSE, title = settings$ylab[2],
                     ticks = "", showline = FALSE)

      se.text <- paste0("Measured value:</br>",
                        data$x, " &plusmn; ", data$y,"</br>")

      hist <- plotly::add_trace(hist, data = data, x = ~x, y = ~y,
                                inherit = FALSE,
                                type = "scatter", mode = "markers",
                                name = "Error", hoverinfo = "text",
                                text = se.text,
                                marker = list(color = "black"),
                                yaxis = "y2")

      hist <- plotly::layout(hist, yaxis2 = yaxis2)
    }

    # set layout ----
    hist <- plotly::layout(hist, hovermode = "closest",
                           title = paste("<b>", settings$main, "</b>"),
                           margin = list(r = 90),
                           xaxis = list(title = settings$xlab,
                                        ticks = ""),
                           yaxis = list(title = settings$ylab[1],
                                        ticks = "",
                                        showline = FALSE,
                                        showgrid = FALSE)
    )

    ## show and return plot ----
    print(hist)
    return(hist)
  }
}