#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.