R/crs-cic.R

In paleolimbot/pb210: Lead-210 dating utilities

#' Apply the Constant Rate of Supply or Constant Initial Concentration model
#'
#' The CRS model was first published by Robbins (1978) and Appleby and Oldfield (1978),
#' and is behind nearly every lead-210-based age-depth model. These functions
#' compute age-depth models based on excess lead-210 activities (perhaps
#' calculated by [pb210_excess()]) and propagate error in quadrature using the
#' [errors::errors()] package. For a more robust estimation of error, consider
#' using [pb210_crs_monte_carlo()] or [pb210_cic_monte_carlo()].
#'
#' @param cumulative_dry_mass The cumulative dry mass of the core (in kg), starting at the
#'   surface sample and including all samples in the core.
#'   These must be greater than 0 and in increasing order.
#' @param excess An excess (non-erosional) lead-210 specific activity (in Bq/kg)
#'   for samples where this was measured, and NA where lead-210 was not measured. Use
#'   [errors::set_errors()] to use quadrature error propogation.
#' @param inventory The cumulative excess lead-210 activity (in Bq), starting at the bottom
#'   of the core. By default, this is estimated by the default [pb210_inventory_calculator()].
#'   If specifying a vector of values, ensure  that the surface (0 cumulative mass) value is
#'   specified.
#' @param model_top A fit object, such as one generated by [pb210_fit_exponential()] or a
#'   constant specifying the surface `excess`. The choice of this value has
#'   considerable impact on young dates.
#' @param core_area The internal area of the corer (in m^2^). This can be calculated
#'   from an internal diameter using [pb210_core_area()].
#' @param decay_constant The decay contstant for lead-210 (in 1/years). This is an argument
#'   rather than a constant because we have found that different spreadsheets in the wild
#'   use different decay constants. See [pb210_decay_constant()].
#' @param object A fit object generated by [pb210_crs()] or [pb210_cic()].
#' @param ... Unused.
#'
#' @references
#' Appleby, P.G., and Oldfield, F. 1983. The assessment of ^210^Pb data from sites with
#' varying sediment accumulation rates.
#' Hydrobiologia, 103: 29–35. <https://doi.org/10.1007/BF00028424>
#'
#' Appleby, P.G., and Oldfield, F. 1978. The calculation of lead-210 dates assuming a
#' constant rate of supply of unsupported ^210^Pb to the sediment.
#' CATENA, 5: 1–8. <https://doi.org/10.1016/S0341-8162(78)80002-2>
#'
#' Robbins, J.A. 1978. Geochemical and geophysical applications of radioactive
#' lead isotopes. In The Biogeochemistry of lead in the environment.
#' Edited by J.O. Nriagu. Elsevier/North-Holland Biomedical Press, Amsterdam.
#' pp. 285–393. <https://books.google.com/books?id=N4wMAQAAMAAJ>
#'
#' @return `predict()` methods return a tibble with (at least)
#' components `age` and `age_sd` (both in years).
#'   CRS model `predict()` function output also contains `inventory`, `inventory_sd`,
#'   `mar` and `mar_sd` (in kg / m^2^ / year).
#' @export
#'
#' @examples
#' # simulate a core
#' core <- pb210_simulate_core() %>%
#'   pb210_simulate_counting()
#'
#' # calculate ages using the CRS model
#' crs <- pb210_crs(
#'   pb210_cumulative_mass(core$slice_mass),
#'   set_errors(
#'     core$activity_estimate,
#'     core$activity_se
#'   )
#' )
#'
#' predict(crs)
#'
pb210_cic <- function(cumulative_dry_mass, excess,
                      model_top = ~pb210_fit_exponential(..1, ..2),
                      decay_constant = pb210_decay_constant()) {
  check_mass_and_activity(cumulative_dry_mass, without_errors(excess))
  stopifnot(
    is.numeric(decay_constant), length(decay_constant) == 1,
    without_errors(decay_constant) > 0
  )

  # it's about 2x faster to do these calculations if we don't have to consider
  # error. this makes a differece during MC simulations
  use_errors <- any(is.finite(extract_errors(excess)))

  data <- tibble::tibble(
    cumulative_dry_mass = cumulative_dry_mass,
    excess = excess
  )

  # resolve the top model  and surface excess
  model_top <- pb210_as_fit(pb210_as_fit(model_top), data = data)
  surface_excess <- stats::predict(model_top, tibble::tibble(x = 0))

  # capture relevant info
  structure(
    list(
      data = data,
      surface_excess = surface_excess,
      use_errors = use_errors,
      decay_constant = decay_constant
    ),
    class = c("pb210_fit_cic", "pb210_fit")
  )
}

#' @rdname pb210_cic
#' @export
pb210_crs <- function(cumulative_dry_mass, excess,
                      inventory = pb210_inventory_calculator(),
                      core_area = pb210_core_area(),
                      decay_constant = pb210_decay_constant()) {

  check_mass_and_activity(cumulative_dry_mass, without_errors(excess))
  stopifnot(
    inherits(inventory, "inventory_calculator") || is.numeric(inventory),
    length(core_area) == 1, is.numeric(core_area),
    length(decay_constant) == 1, is.numeric(decay_constant)
  )

  # it's about 2x faster to do these calculations if we don't have to consider
  # error. this makes a differece during MC simulations
  use_errors <- any(is.finite(extract_errors(excess)))

  if (use_errors) {
    excess <- with_errors(excess)
    na <- with_errors(NA_real_)
    if (!inherits(inventory, "inventory_calculator")) {
      inventory <- with_errors(inventory)
      max_inventory <- with_errors(max_finite(inventory))
    }
  } else {
    excess <- without_errors(excess)
    na <- NA_real_
    if (!inherits(inventory, "inventory_calculator")) {
      inventory <- without_errors(inventory)
      max_inventory <- max(inventory, na.rm = TRUE)
    }
  }

  # need a surface value in cumulative dry mass for calculation purposes
  if (0 %in% cumulative_dry_mass) {
    cumulative_dry_mass_calc <- cumulative_dry_mass
    excess_calc <- excess
  } else {
    cumulative_dry_mass_calc <- c(0, cumulative_dry_mass)
    excess_calc <- c(na, excess)
    if (!inherits(inventory, "inventory_calculator")) {
      # this is an issue with manually specified inventories...not possible to specify
      # surface inventory currently (assuming max() for now)
      inventory <- c(max_inventory, inventory)
    }
  }

  if (inherits(inventory, "inventory_calculator")) {
    inventory <- predict.inventory_calculator(
      inventory,
      cumulative_dry_mass = cumulative_dry_mass_calc,
      excess = excess_calc
    )

    # inventory errors may be introduced by the inventory calculator
    use_errors <- use_errors || any(is.finite(extract_errors(inventory)))
  }

  # check inventory
  stopifnot(
    is.numeric(inventory),
    length(inventory) == length(cumulative_dry_mass_calc),
    # inventory must be decreasing everywhere
    all(diff(without_errors(inventory[is.finite(inventory)])) <= 0)
  )

  # the CRS model is the CIC model using inventory instead of concentration
  fit_cic <- pb210_cic(
    cumulative_dry_mass_calc,
    inventory,
    model_top = without_errors(inventory[1]),
    decay_constant = decay_constant
  )

  # capture relevant info
  structure(
    list(
      data = tibble::tibble(
        cumulative_dry_mass = cumulative_dry_mass,
        excess = excess
      ),
      data_calc = tibble::tibble(
        cumulative_dry_mass = cumulative_dry_mass_calc,
        excess = excess_calc,
        inventory = inventory
      ),
      fit_cic = fit_cic,
      use_errors = use_errors,
      core_area = core_area
    ),
    class = c("pb210_fit_crs", "pb210_fit")
  )
}

#' @rdname pb210_cic
#' @export
predict.pb210_fit_cic <- function(object, cumulative_dry_mass = NULL, ...) {
  if (is.null(cumulative_dry_mass)) {
    cumulative_dry_mass <- object$data$cumulative_dry_mass
  }

  if (object$use_errors) {
    # interpolate excess for all mass values
    excess <- approx_error(
      object$data$cumulative_dry_mass,
      with_errors(object$data$excess),
      xout = cumulative_dry_mass
    )

    one <- with_errors(1)
    decay_constant <- with_errors(object$decay_constant)
    surface_excess <- with_errors(object$surface_excess)
    excess <- with_errors(excess)
  } else {
    # interpolate excess for all mass values
    excess <- approx_no_error(
      object$data$cumulative_dry_mass,
      object$data$excess,
      xout = cumulative_dry_mass
    )

    one <- 1
    surface_excess <- without_errors(object$surface_excess)
    decay_constant <- without_errors(object$decay_constant)
  }

  # calculate ages
  age <- one / decay_constant * log(surface_excess / excess)

  if (object$use_errors) {
    age_sd <- extract_errors(age)
    has_error <- is.finite(errors(excess)) & (errors(excess) > 0)
    age_sd[!has_error] <- NA_real_
    age <- without_errors(age)
  } else {
    age_sd <- NA_real_
  }

  tibble::tibble(age = age, age_sd = age_sd)
}

#' @rdname pb210_cic
#' @export
predict.pb210_fit_crs <- function(object, cumulative_dry_mass = NULL, ...) {
  if (is.null(cumulative_dry_mass)) {
    cumulative_dry_mass <- object$data$cumulative_dry_mass
  }

  # calculate ages using the CIC fit
  ages <- predict.pb210_fit_cic(object$fit_cic, cumulative_dry_mass)

  # the CRS model lets us estimate the mass accumulation rate directly
  # this also calculates the inventory specific to each point
  if (object$use_errors) {
    excess <- approx_error(
      object$data_calc$cumulative_dry_mass,
      object$data_calc$excess,
      cumulative_dry_mass
    )
    inventory <- approx_error(
      object$data_calc$cumulative_dry_mass,
      object$data_calc$inventory,
      cumulative_dry_mass
    )
    has_error <- is.finite(extract_errors(excess)) & (extract_errors(excess) > 0)

    mar <- with_errors(object$fit_cic$decay_constant) *
      inventory / excess / with_errors(object$core_area)

    ages$mar <- without_errors(mar)
    ages$mar_sd <- extract_errors(mar)
    ages$mar_sd[!has_error] <- NA_real_
  } else {
    excess <- approx_no_error(
      object$data_calc$cumulative_dry_mass,
      object$data_calc$excess,
      cumulative_dry_mass
    )
    inventory <- approx_no_error(
      object$data_calc$cumulative_dry_mass,
      object$data_calc$inventory,
      cumulative_dry_mass
    )

    ages$mar <- without_errors(object$fit_cic$decay_constant) *
      without_errors(inventory) / without_errors(excess) / object$core_area
    ages$mar_sd <- NA_real_
  }

  ages$inventory <- without_errors(inventory)
  ages$inventory_sd <- extract_errors(inventory)

  ages
}