R/calc_SedimentationRate.R

In tzerk/leadR: Sediment Deposition and Accumulation Rates Determined by 210Pb Isotope Methods

#' CF:CS model to estimate sedimentation rate from 210Pb activity
#'
#' This function fits a "Constant Flux Constant Sedimentation Rate" (CF:CS)
#' model to 210Pb activity data from different soil depths
#' to estimate the sedimentation rate of the profile.
#'
#' **Model**
#'
#' \deqn{
#'   y = a0 * exp((l/s) * x)
#' }
#' - `y` = 210Pb activity at depth `x` (cm)
#' - `a0` = 210Pb activity at depth `x=0`
#' - `l` = 210Pb decay constant (= 0.03114 a-1)
#' - `s` = sedimentation rate (cm/a)
#' - `x` = depth below surface (cm)
#'
#' @param x [`data.frame`] (**required**):
#' Measured 210Pb activity at different soil depths with the following
#' structure:
#'
#'
#'```
#'      | depth (cm )| activity (Bq/kg) |
#'      |     [ ,1]  |       [ ,2]      |
#'      |------------|------------------|
#' [1, ]|    ~~~~    |       ~~~~       |
#' [2, ]|    ~~~~    |       ~~~~       |
#'  ... |     ...    |        ...       |
#' [x, ]|    ~~~~    |       ~~~~       |
#' ```
#'
#' @param reverse [`logical`] (**optional**):
#' Plot the y-axis in ascending order so that the soil surface is on top.
#' Defaults to `TRUE`.
#'
#' @param fix_a0 [`logical`] (**optional**):
#' In the fitted model parameter `a0` is the 210Pb activity at depth `z=0`.
#' If the argument `fix_a0` is `TRUE`, this function will take the first
#' activity value (i.e., the value closest to the surface) as a fixed value
#' in the equation. For `fix_a0=FALSE` parameter `a0` is estimated by the model.
#' Defaults to `FALSE`.
#'
#' @param verbose [`logical`] (**optional**):
#'
#' @param ... Futher arguments passed to [`plot`].
#'
#' @return
#'
#' A [`list`] with the following structure:
#'
#' \tabular{lll}{
#'  Element \tab Data type \tab Description \cr
#' `$summary` \tab [`data.frame`] \tab summary table of important model results \cr
#' `$data`    \tab [`data.frame`] \tab original input data \cr
#' `$model`    \tab [`nls`] \tab model output generated by [`minpack.lm::nlsLM`]
#' }
#'
#' @author
#' Christoph Burow, University of Cologne
#'
#' @references
#' Szmytkiewicz, A., Zalewska, T., 2014. Sediment deposition and accumulation
#' rates determined by sediment trap and 210Pb isotope methods in the
#' Outer Puck Bay (Baltic Sea). Oceanologia 56(1), 85-106.
#'
#' @examples
#'
#' # Load example data (synthetic)
#' data(Pb)
#'
#' # Apply the model
#' results <- calc_SedimentationRate(x = Pb,
#'                                   reverse = TRUE,
#'                                   fix_a0 = FALSE,
#'                                   verbose = TRUE)
#'
#'
#' @export
#' @md
calc_SedimentationRate <- function(x, reverse = TRUE, fix_a0 = FALSE, verbose = TRUE, ...) {

  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## PRE-CHECKS
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##

  # Basic checks for data.frames
  if (!is.data.frame(x))
    stop("'x' must be a data.frame.", call. = FALSE)
  if (ncol(x) != 2)
    stop("'x' must have exactly two columns (depth, activity).", call. = FALSE)

  # Set column names for more descriptive code
  colnames(x) <- c("depth", "activity")

  # Make sure that the data.frame is sorted by depth in ascending order
  # This is especially important when fixing a0, as we naively assume that the
  # first activity value is from the surface (or closest to)
  x <- x[order(x$depth), ]


  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## DEFAULT SETTINGS
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  settings <- list(
    main = "",
    ylim = if (reverse) rev(range(pretty(x$depth))) else range(pretty(x$depth)),
    xlim = range(pretty(x$activity)),
    xlab = "210Pb activity (Bq/kg)",
    ylab = "Depth (cm)",
    pch = 16,
    col = "black",
    cex = 1.1
  )

  settings <- modifyList(settings, list(...))


  ##
  if (fix_a0)
    a0 <- x$activity[1]
  # if ("a0" %in% names(list(...)))
  #   a0 <- list(...)$a0

  ## 210Pb decay constant
  l <- 0.03114

  ## Fitted model equation
  equation <- formula(activity ~ a0 * exp( (-l/s) * depth ))


  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## CALCULATIONS
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##

  # Start values
  start <- list(
    s = 1
  )
  if (!fix_a0)
    start <- modifyList(start,
                        list(a0 = 100))

  # (un)constrained fitting
  fit <- tryCatch({
    minpack.lm::nlsLM(formula = equation,
                      data = x,
                      start = start,
                      control = list(maxiter = 1024, maxfev = 5000))
  }, error = function(e) { e }
  )


  ## Get standard errors
  if (!inherits(fit, "error")) {

    # estimates
    s <- coef(fit)["s"]
    a0 <- ifelse(fix_a0, x$activity[1], coef(fit)["a0"])

    s_se <- summary(fit)$coefficients["s", "Std. Error"]

    # standard errors
    if (!fix_a0)
      a0_se <- summary(fit)$coefficients["a0", "Std. Error"]
    else
      a0_se <- NA

    # fit failed
  } else {
    a0 <- ifelse(fix_a0, x$activity[1], NA)
    a0_se <- NA
    s <- NA
    s_se <- NA
  }

  summary <- data.frame(
    a0_fixed = fix_a0,
    a0 = a0,
    a0_se = a0_se,
    s = s,
    s_se = s_se
  )

  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## PLOT
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##

  ## Save & recover plot parameters
  par.old <- par(no.readonly = TRUE)
  on.exit(par(par.old))

  ## Set margins
  par(mar = c(2, 6, 6, 5) + 0.1)
  par(cex = settings$cex)

  ## Main plot
  plot(
    x = x$activity,
    y = x$depth,
    ylab = settings$ylab,
    xlab = "",
    xlim = settings$xlim,
    ylim = settings$ylim,
    pch = settings$pch,
    col = settings$col,
    xaxt = "n"
  )

  ## Add x-axis
  axis(3)
  mtext(settings$xlab, side = 3, line = 3, cex = settings$cex)

  # Add fitted model
  if (!inherits(fit, "error")) {
    newx <- seq(0, range(x[ ,1])[2], length.out = 50)
    newy <- suppressWarnings(predict(fit, newdata = list(depth = newx)))
    points(newy, newx, type = "l", col = "red")
  }

  ## Add model parameters
  legend("bottomright",
         legend = c(
           paste0("Activity at z=0 (Bq/kg): ",
                  round(summary$a0, 2), " \u00b1 ", round(summary$a0_se, 2)
           ),
           paste0("Sedimentation rate (cm/a): ",
                  round(summary$s, 2), " \u00b1 ", round(summary$s_se, 2))
         )
  )

  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## CONSOLE OUTPUT
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  if (verbose) {
    f_out <- substitute(activity ~ a0 * exp((l/s) * depth ),
                        list(a0 = round(summary$a0, 3),
                             s = round(summary$s, 3),
                             l = l)
    )

    cat("\n")
    message("Equation: ", capture.output(f_out))
    message("Sedimentation rate: ", round(summary$s, 3), " \u00b1 ", round(summary$s_se, 3), " cm/a")
    cat("\n")
  }

  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  ## RETURN VALUE
  ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##
  output <- list(
    summary = summary,
    data = x,
    model = fit
  )

  return(output)
}