R/PlotCalendarAgeDensityIndividualSample.R

In carbondate: Calibration and Summarisation of Radiocarbon Dates

#' Plot Posterior Calendar Age Estimate for an Individual Determination after Joint Calibration
#'
#' @description
#' Once a joint calibration function (any of [carbondate::PolyaUrnBivarDirichlet], [carbondate::WalkerBivarDirichlet]
#' or [carbondate::PPcalibrate]) has been run to calibrate a set of related radiocarbon
#' determinations, this function plots the posterior calendar age estimate for a given single
#' determination. Shown are a (direct) histogram of the posterior calendar ages generated by the MCMC
#' chain and also a (smoothed) kernel density estimate obtained using a Gaussian kernel. The highest
#' posterior density (HPD) interval is also shown for the interval width specified (default 2\eqn{\sigma}).
#'
#' For more information read the vignettes: \cr
#' \code{vignette("Non-parametric-summed-density", package = "carbondate")} \cr
#' \code{vignette("Poisson-process-modelling", package = "carbondate")}
#'
#' \strong{Note:} The output of this function will provide different results from independent calibration
#' of the determination. By jointly, and simultaneously, calibrating all the related \eqn{{}^{14}}C determinations
#' using the library functions we are able to share the available calendar information between the samples. This should result in improved
#' individual calibration.
#'
#' @inheritParams PlotPredictiveCalendarAgeDensity
#' @param ident The individual determination for which you want to plot the posterior
#' density estimate of the calendar age.
#' @param output_data The return value either from one of the Bayesian non-parametric DPMM
#' functions ([carbondate::PolyaUrnBivarDirichlet] or [carbondate::WalkerBivarDirichlet]); or
#' from the Poisson process modelling function ([carbondate::PPcalibrate]).
#' @param hist_resolution The distance between histogram breaks when plotting the individual
#' posterior calendar age density. Default is 5.
#' @param density_resolution The distance between calendar ages for the returned smoothed calendar age
#' probability. Default is 1.
#' @param interval_width The confidence intervals to show for the calibration curve and
#' for the highest posterior density ranges.
#' Choose from one of "1sigma" (68.3%), "2sigma" (95.4%) and "bespoke". Default is "2sigma".
#' @param bespoke_probability The probability to use for the confidence interval
#' if `"bespoke"` is chosen above. E.g., if 0.95 is chosen, then the 95% confidence
#' interval is calculated. Ignored if `"bespoke"` is not chosen.
#' @param show_hpd_ranges Set to `TRUE` to also show the highest posterior density (HPD) range
#' on the plot.
#' @param show_unmodelled_density Set to `TRUE` to also show the unmodelled density (i.e., the
#' result of independent calibration using [carbondate::CalibrateSingleDetermination]) on the plot.
#' Default is `FALSE`.
#' @param plot_pretty logical, defaulting to `TRUE`. If set `TRUE` then will select pretty plotting
#' margins (that create sufficient space for axis titles and rotates y-axis labels). If `FALSE` will
#' implement current user values.
#' @param n_end The iteration number of the last sample in `output_data` to use in the calculations.
#' Assumed to be the total number of (thinned) realisations stored if not given.
#'
#' @export
#'
#' @return A data frame with one column `calendar_age_BP` containing the calendar ages, and the other
#' column `probability` containing the (smoothed) kernel density estimate of the probability at that
#' calendar age.
#'
#' @seealso [carbondate::CalibrateSingleDetermination] for independent calibration of a sample against
#' a calibration curve.
#'
#' @examples
#' # NOTE 1: These examples are shown with a small n_iter to speed up execution.
#' # When you run ensure n_iter gives convergence (try function default).
#'
#' # NOTE 2: The examples only show application to  PolyaUrnBivarDirichlet output.
#' # The function can also be used with WalkerBivarDirichlet and PPcalibrate output.
#'
#' polya_urn_output <- PolyaUrnBivarDirichlet(
#'     two_normals$c14_age,
#'     two_normals$c14_sig,
#'     intcal20,
#'     n_iter = 500,
#'     n_thin = 2,
#'     show_progress = FALSE)
#'
#' # Result for 15th determination
#' PlotCalendarAgeDensityIndividualSample(15, polya_urn_output)
#'
#' # Now change to show 1 sigma interval for HPD range and calibration curve
#' # and plot in yr AD
#' PlotCalendarAgeDensityIndividualSample(
#'     15,
#'     polya_urn_output,
#'     plot_cal_age_scale = "AD",
#'     interval_width = "1sigma",
#'     show_hpd_ranges = TRUE,
#'     show_unmodelled_density = TRUE)
#'
#' # Plot and then assign the returned probability
#' posterior_dens <- PlotCalendarAgeDensityIndividualSample(15, polya_urn_output)
#' # Use this to find the mean posterior calendar age
#' weighted.mean(posterior_dens$calendar_age_BP, posterior_dens$probability)
PlotCalendarAgeDensityIndividualSample <- function(
    ident,
    output_data,
    calibration_curve = NULL,
    plot_14C_age = TRUE,
    plot_cal_age_scale = "BP",
    hist_resolution = 5,
    density_resolution = 1,
    interval_width = "2sigma",
    bespoke_probability = NA,
    n_burn = NA,
    n_end = NA,
    show_hpd_ranges = FALSE,
    show_unmodelled_density = FALSE,
    plot_pretty = TRUE) {

  arg_check <- .InitializeErrorList()
  .CheckInteger(arg_check, ident)
  .CheckOutputData(arg_check, output_data,  c("Polya Urn", "Walker", "RJPP"))
  n_iter <- output_data$input_parameters$n_iter
  n_thin <- output_data$input_parameters$n_thin

  .CheckCalibrationCurveFromOutput(arg_check, output_data, calibration_curve)
  .CheckNumber(arg_check, hist_resolution, lower = 0.01)
  .CheckNumber(arg_check, density_resolution, lower = 0.01)
  .CheckChoice(arg_check, plot_cal_age_scale, c("BP", "AD", "BC"))
  .CheckNBurnAndNEnd(arg_check, n_burn, n_end, n_iter, n_thin)
  .ReportErrors(arg_check)

  # Ensure revert to main environment par on exit of function
  oldpar <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(oldpar))

  # Set nice plotting parameters
  if(plot_pretty) {
    graphics::par(
      mgp = c(3, 0.7, 0),
      xaxs = "i",
      yaxs = "i",
      mar = c(5, 4.5, 4, 2) + 0.1,
      las = 1)
  }

  n_burn <- .SetNBurn(n_burn, n_iter, n_thin)
  n_end <- .SetNEnd(n_end, n_iter, n_thin)

  if (is.null(calibration_curve)) {
    calibration_curve <- get(output_data$input_data$calibration_curve_name)
  }
  rc_determinations <- output_data$input_data$rc_determinations
  rc_sigmas <- output_data$input_data$rc_sigmas
  F14C_inputs <-output_data$input_data$F14C_inputs

  if (plot_14C_age == TRUE) {
    calibration_curve <- .AddC14ageColumns(calibration_curve)
    if (F14C_inputs == TRUE) {
      converted <- .ConvertF14cTo14Cage(rc_determinations, rc_sigmas)
      rc_determinations <- converted$c14_age
      rc_sigmas <- converted$c14_sig
    }
  } else {
    calibration_curve <- .AddF14cColumns(calibration_curve)
    if (F14C_inputs == FALSE) {
      converted <- .Convert14CageToF14c(rc_determinations, rc_sigmas)
      rc_determinations <- converted$f14c
      rc_sigmas <- converted$f14c_sig
    }
  }

  calendar_age_BP <- output_data$calendar_ages[(n_burn+1):n_end, ident]
  rc_age <- rc_determinations[ident]
  rc_sig <- rc_sigmas[ident]

  if(output_data$update_type == "RJPP") {
    bandwidth_selector <- "nrd0" # As all calendar ages are integers
  } else {
    bandwidth_selector <- "SJ" # As continuous
  }

  smoothed_density <- stats::density(calendar_age_BP, bw = bandwidth_selector)
  xrange_BP <- range(calendar_age_BP)

  # Interpolate for the smoothed density to return
  returned_calendar_age_BP <- seq(
    from = ceiling(xrange_BP[1] / density_resolution) * density_resolution,
    to = floor(xrange_BP[2] / density_resolution) * density_resolution,
    by = density_resolution)
  returned_density <- data.frame(
    calendar_age_BP = returned_calendar_age_BP,
    probability = stats::approx(smoothed_density$x, smoothed_density$y, returned_calendar_age_BP)$y)

  # Expand the calendar age range for plotting
  xrange_BP <- xrange_BP + 0.1 * c(-1, 1) * diff(xrange_BP)
  xrange_BP[1] <- floor(xrange_BP[1] / hist_resolution) * hist_resolution
  xrange_BP[2] <- ceiling(xrange_BP[2] / hist_resolution) * hist_resolution

  if (plot_14C_age == FALSE) {
    title <- substitute(
      paste(
        "Posterior of ",
        i^th,
        " determination ",
        f14c_age,
        "\u00B1",
        f14c_sig,
        "F"^14,
        "C"),
      list(i = ident, f14c_age = signif(rc_age, 2), f14c_sig = signif(rc_sig, 2)))
  } else {
    title <- substitute(
      paste(
        "Posterior of ",
        i^th,
        " determination ",
        c14_age,
        "\u00B1",
        c14_sig,
        ""^14,
        "C yr BP"),
      list(i = ident, c14_age = round(rc_age), c14_sig = round(rc_sig, 1)))
  }

  .PlotCalibrationCurve(
    plot_cal_age_scale,
    xlim = rev(xrange_BP),
    plot_14C_age = plot_14C_age,
    calibration_curve = calibration_curve,
    calibration_curve_colour = "blue",
    calibration_curve_bg = grDevices::rgb(0, 0, 1, .3),
    interval_width = interval_width,
    bespoke_probability = bespoke_probability,
    title = title)

  calendar_age <- .ConvertCalendarAge(plot_cal_age_scale, calendar_age_BP)
  xrange <- .ConvertCalendarAge(plot_cal_age_scale, xrange_BP)
  smoothed_density$x <- .ConvertCalendarAge(plot_cal_age_scale, smoothed_density$x)

  # Plot the 14C determination on the y-axis
  yfromto <- seq(rc_age - 4 * rc_sig, rc_age + 4 * rc_sig, length.out = 100)
  radpol <- cbind(
    c(0, stats::dnorm(yfromto, mean =rc_age, sd = rc_sig), 0),
    c(min(yfromto), yfromto, max(yfromto))
  )
  relative_height <- 0.1
  radpol[, 1] <- radpol[, 1] * (xrange[2] - xrange[1]) / max(radpol[, 1])
  radpol[, 1] <- radpol[, 1] * relative_height
  radpol[, 1] <- xrange[2] - radpol[, 1]
  graphics::polygon(radpol, col = grDevices::rgb(1, 0, 0, .5))

  # Plot the posterior cal age on the x-axis
  graphics::par(new = TRUE)
  # Create hist but do not plot - works out sensible ylim
  if (plot_cal_age_scale == "AD") {
    breaks <-seq(xrange[2], xrange[1], by=hist_resolution)
  } else {
    breaks <-seq(xrange[1], xrange[2], by=hist_resolution)
  }
  temphist <- graphics::hist(calendar_age, breaks = breaks, plot = FALSE)

  graphics::hist(
    calendar_age,
    prob = TRUE,
    breaks = breaks,
    xlim = rev(xrange),
    axes = FALSE,
    xlab = NA,
    ylab = NA,
    main = "",
    xaxs = "i",
    ylim = c(0, 3 * max(temphist$density)),
    col = "lavender",
    border = "purple")
  graphics::lines(smoothed_density, lwd=2, col='purple', lty=2)

  if (show_unmodelled_density) {
    unmodelled <- CalibrateSingleDetermination(
      rc_age, rc_sig, calibration_curve, resolution = density_resolution)
    unmodelled$calendar_age <- .ConvertCalendarAge(plot_cal_age_scale, unmodelled$calendar_age_BP)

    graphics::polygon(
      c(unmodelled$calendar_age, rev(unmodelled$calendar_age)),
      c(unmodelled$probability, rep(0, length(unmodelled$probability))),
      border = NA,
      col = grDevices::grey(0.1, alpha = 0.25))
  }
  if (show_hpd_ranges) {
    hpd_probability <- switch(
      interval_width,
      "1sigma" = 0.682,
      "2sigma" = 0.954,
      "bespoke" = bespoke_probability)
    hpd <- FindHPD(smoothed_density$x, smoothed_density$y, hpd_probability)
    graphics::segments(hpd$start_ages, hpd$height, hpd$end_ages, hpd$height, lwd=4, col='black', lend='butt')
  }

  .AddLegendToDensitySamplePlot(
    interval_width,
    bespoke_probability,
    output_data$input_data$calibration_curve_name,
    show_hpd_ranges,
    show_unmodelled_density,
    hpd)

  invisible(returned_density)
}

.AddLegendToDensitySamplePlot <- function(
    interval_width,
    bespoke_probability,
    calcurve_name,
    show_hpd_ranges,
    show_unmodelled_density,
    hpd){
  ci_label <- switch(
    interval_width,
    "1sigma" = expression(paste("1", sigma, " interval")),
    "2sigma"  = expression(paste("2", sigma, " interval")),
    "bespoke" = paste0(round(100 * bespoke_probability), "% interval"))

  legend_labels <- c(
    substitute(paste(""^14, "C determination ")),
    gsub("intcal", "IntCal", gsub("shcal", "SHCal", calcurve_name)), # Both IntCal and SHCal
    ci_label)
  lty <- c(-1, 1, 2)
  lwd <- c(-1, 1, 1)
  pch <- c(15, NA, NA)
  col <- c(grDevices::rgb(1, 0, 0, .5), "blue", "blue")

  if (show_unmodelled_density) {
    legend_labels <- c(legend_labels, "Unmodelled density")
    lty <- c(lty, -1)
    lwd <- c(lwd, -1)
    pch <- c(pch, 15)
    col <- c(col, grDevices::grey(0.1, alpha = 0.25))
  }

  legend_labels <- c(legend_labels, "Posterior density")
  lty <- c(lty, 1)
  lwd <- c(lwd, 1)
  pch <- c(pch, NA)
  col <- c(col, "purple")

  if (show_hpd_ranges) {
    hpd_label <- switch(
      interval_width,
      "1sigma" = "68.3% HPD",
      "2sigma"  = "95.4% HPD",
      "bespoke" = paste0(round(100 * bespoke_probability, 3), "% HPD"))
    legend_labels <- c(legend_labels, hpd_label)
    lty <- c(lty, 1)
    lwd <- c(lwd, 2)
    pch <- c(pch, NA)
    col <- c(col, "black")

    for (i in rev(seq_along(hpd$start_ages))) {
      auc_string <- paste0("(", round(hpd$area_under_curve[i] * 100, digits = 1), "%)")
      legend_labels <- c(
        legend_labels,
        paste("   ", round(hpd$end_ages[i]), "-", round(hpd$start_ages[i]), auc_string))
      lty <- c(lty, NA)
      lwd <- c(lwd, NA)
      pch <- c(pch, NA)
      col <- c(col, NA)
    }
  }

  graphics::legend(
    "topright", legend = legend_labels, lty = lty, lwd=lwd, pch = pch, col = col)

}