R/clam.R
In Maarten14C/clam: Classical Age-Depth Modelling of Cores from Deposits
Defines functions
.model.clam
clam
Documented in
clam
#' clam
#'
#' Classical (non-Bayesian) age-depth modelling.
#'
#' @docType package
#' @author Maarten Blaauw <maarten.blaauw@qub.ac.uk>
#' @importFrom grDevices dev.off grey pdf png rgb
#' @importFrom graphics abline image layout legend lines par points polygon rect plot
#' @importFrom stats approx density dnorm lm loess pnorm predict qnorm quantile rnorm runif smooth.spline spline weighted.mean
#' @importFrom utils packageName read.csv read.table write.table
#' @import rintcal
#' @name clam
NULL
# do:
# done:
# so that functions such as pMC.age etc. are available immediately
library(rintcal)
#' @name clam
#' @title The main age-depth modelling function
#' @description Produce age-depth models for cores with dated depths.
#' @details
#' Cores containing several 14C and/or other dates can be processed semi-automatically in order to obtain age-depth models.
#' In the process, any 14C dates are calibrated, and age-depth curves are repeatedly drawn through point estimates sampled from the dates.
#' Age-depth models can be based on linear interpolation, linear/polynomial regression, or cubic, smooth or locally weighted splines.
#' For each date, the probability of a calendar year being sampled is proportionate to its calibrated probability (see Blaauw, 2010).
#' Uncertainty ranges as well as a 'best' age-model are calculated.
#'
#' Additional cores should be put in a comma-separated file in a sub-folder of the directory where the cores are stored.
#' By default this parent folder is called \code{coredir="clam_runs"} (if no folder called \code{"Cores"} already exists). If your core is called MyCore1, save MyCore1.csv as \code{clam_runs/MyCore1/MyCore1.csv}.
#' Ensure that the names of the core's folder and filename's root (the part before .csv) match, e.g., using exactly similar upper- and lower case letters.
#'
#' Avoid the use of spaces or non-standard (non-ASCII) characters within the file or in folder or file names.
#' The plain text file should consist of 6 or 7 columns (also called fields), containing in the following exact order (see the example below):
#' \enumerate{
#' \item Identification labels (e.g. 14C lab codes)
#' \item 14C ages for 14C-dated depths; leave empty for non-14C dated depths
#' \item cal BP ages (for any non-14C dates such as the core surface; leave empty for levels with 14C dates)
#' \item errors (reported 1 standard deviation errors. This column should never be left empty. Errors should always be larger than 0)
#' \item age offsets if known (otherwise leave empty)
#' \item depths (depths in the sequence were the dated samples were taken, default unit depth="cm"; this column should never be left empty)
#' \item thicknesses of the sampled slices (optional column; leave empty for default of 1)
#' }
#' Add a final empty line to your core's .csv file by pressing 'Enter' after the file's last value.
#'
#' These files can be made in spreadsheet software such as MS-Excel, but it is always a good idea to check the file's formatting in a plain-text editor such as WordPad. Remove any lines which contain only commas, and it is also recommended to remove quotes ()\code{\" or \'}) in the headers or elsewhere.
#'
#' Age-models for the core can then be produced by typing, e.g., \code{clam("MyCore1")}.
#'
#' By default the northern hemisphere terrestrial calibration curve is used (\code{cc=1}, \code{cc1="3Col_intcal20.14C"}).
#' To use alternative curves, change \code{cc} to \code{cc=2 (cc2="3Col_marine20.14C")}, \code{cc=3 (cc3="3Col_shcal20.14C")}, \code{cc=4 (cc4="mixed.14C")}.
#' You can also provide custom-built calibration curves, indicating its location using \code{ccdir}.
#'
#' The provided example (default \code{core="Example"}) is core Quilichao-1 which was sampled from a Colombian lake (Berrio et al., 2002).
#' This core was chosen because it was dated at a rather high resolution, and appears to contain a hiatus (e.g., try \code{hiatus=450}
#' for a hiatus at 450 cm depth).
#'
#' Each clam run will produce a range of files within the core's folder. One, ending with \code{"_calibrated.txt"} contains the calibrated
#' age ranges of the 14C and other dates. The others will be named according to the core's name followed by the model type,
#' and contain the age estimates for all depths (files ending with \code{"_ages.txt"}), settings (files ending with \code{"_settings.txt"})
#' and graphs (files ending with \code{".pdf"} and \code{".png"}).
#' The file containing the age estimates has 5 columns; first the depths, then the minima and maxima of the confidence intervals,
#' then a "best" estimate, and finally the reconstructed accumulation rates. The reported values are rounded to 0 decimals by default
#' (\code{decimals=0}). Accumulation rates are in yr/cm ("deposition time") by default (\code{cmyr=FALSE}), but can be reported in cm/yr (\code{cmyr=TRUE}).
#'
#' see Blaauw 2010 (Quaternary Geochronology 5: 512-518).
#'
#' @param core Name of the core, given using quotes. Defaults to the core provided with clam, \code{core="Example"}.
#' @param type The type of age-depth model. Five different types are provided:
#' \enumerate{
#' \item linear interpolation between neighbouring levels (1, "int", "inter" or "interp")
#' \item linear or higher polynomial regression (2, "reg", "regr", "poly" or "polyn", default linear)
#' \item cubic spline (3, "spl" or "spline")
#' \item smooth spline (4, "sm" or "smooth", default smoothing 0.3)
#' \item locally weighted spline (5, "loess" or "lowess", default smoothing 0.75, cannot extrapolate)
#' }
#' @param smooth Degree of smoothing. Gives polynomial degree for model type 2. Not relevant for \code{type=1 or type=3}.
#' \itemize{
#' \item for type=2: \code{smooth=1} (linear), \code{smooth=2} second-order polynomial, \code{smooth=3} for third-order polynomial, etc.
#' \item for type=4: \code{smooth=0.3}
#' \item for type=5: \code{smooth=0.75}
#' }
#' @param prob Confidence intervals (between 0 and 1), default \code{prob=0.95} or 95\%.
#' @param its Amount of age-model iterations; defaults to \code{its=1000}.
#' @param coredir The directory where core runs are stored (each core in its own directory named after the core's name).
#' @param ask By default, and as per R rules, clam will ask if it is OK to make or write to a directory.
#' Defaults to \code{coredir="clam_runs"}, or to \code{coredir="Cores"} if this folder exists where R is working.
#' @param wghts Weights can be applied to dated depths as follows:
#' \itemize{
#' \item 0 no weighting
#' \item 1 weighted to calibrated probabilities of sampled calendar years (default, \code{wghts=1}).
#' \item 2 weighted to (inverse squared) errors of the dates.
#' }
#' @param cc calibration curve for C14 dates (1, 2 or 3).
#' @param cc1 For terrestrial, northern hemisphere C14 dates.
#' @param cc2 For marine C14 dates.
#' @param cc3 For southern hemisphere C14 dates.
#' @param cc4 For mixed terrestrial/marine C14 dates.
#' @param postbomb Use a postbomb curve for negative (i.e. postbomb) 14C ages. \code{0 = none, 1 = NH1, 2 = NH2, 3 = NH3, 4 = SH1-2, 5 = SH3}. See \url{http://calib.org/CALIBomb/}.
#' @param pb1 For Northern hemisphere region 1 postbomb C14 dates.
#' @param pb2 For Northern hemisphere region 2 postbomb C14 dates.
#' @param pb3 For Northern hemisphere region 3 postbomb C14 dates.
#' @param pb4 For Southern hemisphere regions 1-2 postbomb C14 dates.
#' @param pb5 For Southern hemisphere region 3 postbomb C14 dates.
#' @param ccdir Directory where the calibration curves for C14 dates \code{cc} are located. By default \code{ccdir=""}.
#' For example, use \code{ccdir="."} to choose current working directory, or \code{ccdir="Curves/"} to choose sub-folder \code{Curves/}.
#' @param outliers The number of any dates to be considered outlying, e.g. \code{c(5,6)} for the fifth and sixth dated depth counting from the top of a core.
#' @param ignore The number of any dates that should be ignored, e.g., \code{c(5,6)} for the fifth and sixth date counting from the top of a core.
#' @param youngest The age beyond which dates should be truncated (e.g., \code{youngest=-60} if the core was sampled in -60 cal BP or AD 2010).
#' @param extradates Depths of any additional dates with their files of ages and probabilities.
#' @param slump Upper and lower depths of sections of abrupt accumulation that should be excised, e.g., \code{c(600, 550, 120, 100)} for two sections of 600-550 and 120-100 cm depth.
#' @param est Which point estimate to use as 'best' age. It is highly recommended to not only use these 'best' point estimates, as chronological uncertainties are often considerable and should not be ignored.
#' \enumerate{
#' \item averages of age-depth model derived ages (default, \code{est=1})
#' \item midpoints of age-depth model derived age estimates
#' \item midpoints of calibrated ranges
#' \item weighted means of calibrated ranges
#' \item medians of calibrated distributions
#' \item maximum densities of calibrated distributions
#' \item midpoints of entire calibrated distributions (including years outside the calibrated ranges)
#' }
#' @param calibt Calibration based on the student-t distribution. By default, the Gaussian distribution is used (\code{calibt=FALSE}). To use the student-t distribution, provide two parameters such as \code{calibt=c(3,4)}.
#' @param mixed.effect Set to \code{TRUE} to activate mixed-effect modelling.
#' @param dmin Minimum depth of age-depth model (e.g., extrapolate).
#' @param dmax Maximum depth of age-depth model (e.g., extrapolate).
#' @param every Resolution at which (ages for) depths are calculated.
#' @param yrmin Minimum of calendar axis of age-depth plot (calculate automatically by default).
#' @param yrmax Maximum of calendar axis of age-depth plot (calculated automatically by default).
#' @param yrsteps Temporal resolution at which calibrated ages are calculated (in calendar years).
#' @param pbsteps Temporal resolution at which postbomb C14 ages are calibrated (in calendar years).
#' @param hpdsteps Temporal resolution at which highest posterior density ranges are calibrated (in calendar years).
#' @param BCAD Use BC/AD or cal BP scale.
#' @param decimals Amount of decimals for rounding.
#' @param cmyr Accumulation rates can be provided as yr/cm (default, \code{cmyr=TRUE}, more accurately named deposition times) or cm/yr (\code{cmyr=FALSE}).
#' @param ageofdepth Calculate age estimates of a specific depth.
#' @param depth Depth units.
#' @param depthseq Sequence of depths for which age estimates are to be calculated (default: from \code{dmin} to \code{dmax} with steps of size every)
#' @param depths.file Use a file with depths for depthseq.
#' @param thickness Thickness of the dated samples.
#' @param hiatus Depths of any hiatuses, e.g., \code{c(500, 300)}. Each sub-section must have at least 2 dates (\code{4} for smoothing spline; does not work with loess as it cannot extrapolate).
#' @param remove.reverse Proportion of age-models with reversals that can be removed before prompting a warning.
#' Set at \code{FALSE} to avoid removing models with reversals.
#' @param times Half-range of calibration curve used to calibrate dates (multiplication factor for the dates' errors).
#' @param sep Separator between the fields of the plain text file containing the dating information.
#' @param ext Extension of the file containing the dating information.
#' @param runname Text to add to the core name for specific runs, e.g., "MyCore_Test1"
#' @param storedat Store the dates and age-model within R after a \code{clam} run. Defaults to \code{storedat=TRUE}.
#' @param threshold Below which value should probabilities be excluded from calculations.
#' @param proxies Set to \code{TRUE} to plot proxies against age after the run.
#' @param revaxes Set to \code{TRUE} to plot ages on the vertical axis and depth on the horizontal axis.
#' @param revd Plot depth axis in reverse.
#' @param revyr Plot age axis in reverse.
#' @param calhght Heights of the calibrated distributions in the age-depth plot.
#' @param maxhght Maximum height of age probability distributions.
#' @param mirror Plot the age distributions in "mirror" style (above and below depth).
#' @param plotrange Plot the confidence ranges of the age-model.
#' @param mar Plot margins (amount of white space along edges of axes 1-4).
#' @param mgp Axis text margins (where should titles, labels and tick marks be plotted).
#' @param bty Type of box to be drawn around plots. Draw a box around the graph (\code{"n"} for none,
#' and \code{"l"}, \code{"7"}, \code{"c"}, \code{"u"}, "]" or \code{"o"} for correspondingly shaped boxes).
#' @param plotpdf Produce a pdf file of the age-depth plot.
#' @param plotpng Produce a png file of the age-depth plot.
#' @param greyscale Produce a grey-scale representation of all age-models (number gives resolution, e.g., 500 bins; will cancel plotting of the confidence intervals).
#' @param yrlab Label of the calendar axis. Defaults to either cal BP or BC/AD. Alternative names can be provided.
#' @param dlab Label of the depth axis. Defaults to \code{dlab="Depth (cm)"} (assuming \code{depth="cm"}), but alternative names can be provided.
#' @param calcol Colour of the calibrated distributions in the age-depth plot.
#' @param C14col Colour of the calibrated ranges of the dates.
#' @param outcol Colour of outlying dates.
#' @param outlsize Size of symbols outlying dates.
#' @param bestcol Colour of the "best" age-depth model (based on chosen value for est).
#' @param rangecol Colour of plotted confidence ranges.
#' @param slumpcol Colour of slump.
#' @param plotname Print the core name on the graph.
#' @param ash Plot all distributions at the same height.
#' @param rule How should R's approx function deal with extrapolation. If \code{rule=1}, the default, then NAs are returned for such points and if it is 2, the value at the closest data extreme is used.
#'
#' @author Maarten Blaauw
#' @return Age model construction together with a text output and files saved to a folder in the \code{coredir/core} directory.
#' @examples
#' clam(, coredir=tempdir()) # Create the example in Cores/Example folder
#' clam(, coredir=tempdir(), extradates=470)
#'
#' @references
#' Berrio, J.C., Hooghiemstra, H., Marchant, R., Rangel, O., 2002. Late-glacial and Holocene history of the dry forest area
#' in the south Colombian Cauca Valley. Journal of Quaternary Science 17, 667-682
#'
#' Blaauw, M., 2010. Methods and code for 'classical' age-modelling of radiocarbon sequences. Quaternary Geochronology 5, 512-518
#' \doi{10.1016/j.quageo.2010.01.002}
#' @export
clam <- function(core="Example", type=1, smooth=NULL, prob=0.95, its=1000, coredir=NULL, ask=TRUE, wghts=1, cc=1, cc1="3Col_intcal20.14C", cc2="3Col_marine20.14C", cc3="3Col_shcal20.14C", cc4="mixed.14C", postbomb=FALSE, pb1="postbomb_NH1.14C", pb2="postbomb_NH2.14C", pb3="postbomb_NH3.14C", pb4="postbomb_SH1-2.14C", pb5="postbomb_SH3.14C", ccdir="", outliers=NULL, ignore=NULL, youngest=NULL, extradates=NULL, slump=NULL, est=1, calibt=FALSE, mixed.effect=FALSE, dmin=NULL, dmax=NULL, every=1, yrmin=NULL, yrmax=NULL, yrsteps=1, pbsteps=0.01, hpdsteps=1, BCAD=FALSE, decimals=0, cmyr=FALSE, ageofdepth=NULL, depth="cm", depthseq=NULL, depths.file=FALSE, thickness=1, hiatus=NULL, remove.reverse=0.5, times=5, sep=",", ext=".csv", runname=NULL, storedat=TRUE, threshold=1e-6, proxies=FALSE, revaxes=FALSE, revd=TRUE, revyr=TRUE, calhght=0.3, maxhght=0.01, mirror=TRUE, plotrange=TRUE, bty="l", mar=c(3.5,3,2,1), mgp=c(2,1,0), plotpdf=TRUE, plotpng=TRUE, greyscale=NULL, yrlab=NULL, dlab=NULL, calcol=rgb(0,0.5,0.5,0.5), C14col=rgb(0,0,1,0.5), outcol="red", outlsize=1, bestcol="black", rangecol=rgb(0,0,0,0.3), slumpcol=grey(0.75), plotname=TRUE, ash=FALSE, rule=1) {
# If coredir is left empty, check for a folder named Cores in the current working directory, and if this doesn't exist, for a folder called clam_runs (make this folder if it doesn't exist yet).
# Check if we have write access. If not, tell the user to provide a different, writeable location for coredir.
coredir <- assign_coredir(coredir, core, ask)
if(core == "Example") {
dir.create(file.path(coredir, "Example/"), showWarnings = FALSE, recursive = TRUE)
fileCopy <- system.file("extdata/Example/", package="clam")
file.copy(fileCopy, coredir, recursive = TRUE, overwrite=FALSE)
}
# warn and stop if unexpected settings are provided
if(type > 5 || type < 1 || prob < 0 || prob > 1 || its < 100 || wghts < 0 || wghts > 1 || est < 1 || est > 7 || yrsteps <= 0 || hpdsteps <= 0 || every <= 0 || decimals < 0 || cmyr < 0 || cmyr > 1 || thickness < 0 || times < 1 || calhght < 0 || (type==5 && length(hiatus)>0))
stop("\n Warning, clam cannot run with these settings! Please check the manual.\n\n", call.=FALSE)
csvFile <- paste0(coredir, core, "/", core, ext)
if(!file.exists(csvFile))
stop(paste("\nInput data file", csvFile, "not found."), call.=FALSE)
dets <- read.csv(csvFile, sep=sep)
d <- dets[,6]
these <- 1:length(d) # to possibly kick out dates later on; Nov 2020
if(min(diff(d)) < 0)
message(" Warning, depths not in ascending order (top ones should come first).")
# avoid Windows/Mac habit of silently adding .txt extension to plain text files
bill <- list.files(file.path(coredir, core), pattern=".csv.txt")
if(length(bill) > 0) {
message("Removing unnecessary .txt extension from .csv file", bill[1])
file.rename(file.path(coredir, core, paste0(core, ".csv.txt")),
file.path(coredir, core, "/", paste0(core, ".csv")))
}
# set the calibration curve
ccdir <- .validateDirectoryName(ccdir)
if(ccdir == "") # so, if no alternative folder provided, use clam's calibration curves
ccdir = system.file("extdata", package="rintcal")
if(cc==1) calcurve <- read.table(file.path(ccdir, cc1)) else
if(cc==2) calcurve <- read.table(file.path(ccdir, cc2)) else
if(cc==3) calcurve <- read.table(file.path(ccdir, cc3)) else
if(cc==4) calcurve <- read.table(file.path(ccdir, cc4)) else
stop("I do not understand which calibration curve you mean, check the manual", call.=FALSE)
if(cc==1) ccname <- cc1 else
if(cc==2) ccname <- cc2 else
if(cc==3) ccname <- cc3 else
if(cc==4) ccname <- cc4
# negative C14 ages need a postbomb curve
pbnames <- c(pb1, pb2, pb3, pb4, pb5)
cdat <- dets[,2]
if(length(cdat[!is.na(cdat)]) > 0)
if(min(cdat[!is.na(cdat)]) < 0)
if(postbomb==FALSE)
message("Warning, negative 14C ages, should I use a postbomb curve?") else {
if(postbomb>5)
stop("I do not understand which postbomb curve you mean, check the manual", call.=FALSE)
yrsteps <- min(pbsteps, yrsteps)
pb <- read.table(system.file("extdata", pbnames[postbomb], package="rintcal")) # now points to rintcal package
pb.x <- seq(min(pb[,1]), max(pb[,1]), by=yrsteps)
pb.y <- approx(pb[,1], pb[,2], pb.x)$y
pb.sd <- approx(pb[,1], pb[,3], pb.x)$y
calcurve <- cbind(c(pb.x, calcurve[,1]), c(pb.y, calcurve[,2]), c(pb.sd, calcurve[,3]))
}
# work in BC/AD if needed, and prepare for calculations in f14C
if(BCAD) {
theta <- 1950-calcurve[,1]
border <- max(which(theta > 0))
theta <- c(theta[1:(border-1)], theta[border]:theta[border+2], theta[(border+3):length(theta)])
mu <- approx(1950-calcurve[,1], calcurve[,2], theta, rule=rule)$y
sigma <- approx(1950-calcurve[,1], calcurve[,3], theta, rule=rule)$y
theta[theta <=0] <- theta[theta <= 0]-1
calcurve <- cbind(theta, mu, sigma)
} else theta <- calcurve[,1]
if(length(yrlab) == 0)
yrlab <- ifelse(BCAD, "cal BC/AD", "cal BP")
f.mu <- exp(-calcurve[,2]/8033)
f.sigma <- exp(-(calcurve[,2]-calcurve[,3])/8033) - f.mu
# prepare for slumps and hiatuses
#cat(paste("Core name Test:", core))
if(length(greyscale) > 0) storedat <- TRUE
if(length(slump) > 0) {
if(length(slump) %% 2 == 1)
stop("\n Warning, slumps need both upper and lower depths. Please check the manual", call.=FALSE)
slump <- matrix(sort(slump), ncol=2, byrow=TRUE)
if(length(dmax)==0)
dmax <- max(dets[,6])
if(length(extradates) > 0)
dmax <- max(dmax, extradates)
for(i in 1:nrow(slump)) {
d[d > min(slump[i,])] <- d[d > min(slump[i,])] - (max(slump[i,]) - min(slump[i,]))
dmax <- dmax - (max(slump[i,])-min(slump[i,]))
}
if(length(hiatus) > 0)
for(i in 1:nrow(slump)) {
below.slump <- which(hiatus > max(slump[i,]))
above.slump <- which(hiatus < min(slump[i,]))
hiatus[below.slump] <- hiatus[below.slump] - (max(slump[i,])-min(slump[i,]))
hiatus <- hiatus[c(above.slump, below.slump)]
}
}
# read in the data
dat <- .read.clam(core, coredir, ext, hpdsteps, yrsteps, prob, times, sep, BCAD, storedat, ignore, thickness, youngest, slump, threshold, theta, f.mu, f.sigma, calibt, extradates, calcurve, postbomb, rule=rule)
message(" Calibrating dates... ")
# calculate the depths to be used, based on the ranges and resolution
if(length(dmin)==0)
dmin <- floor(min(dat$depth))
if(length(dmax)==0)
dmax <- ceiling(max(dat$depth))
if(depths.file)
if(file.exists(dd <- paste0(coredir, core, "/", core, "_depths.txt"))) {
if(length(depthseq) == 0)
depthseq <- seq(dmin, dmax, by=every)
depthseq <- sort(unique(c(depthseq, suppressWarnings(read.table(dd))[,1])))
# dmin <- min(depthseq)#, read.table(dd)[,1])
# dmax <- max(depthseq)#, read.table(dd)[,1])
} else
stop(paste0("\nCannot find file ", dat$core, "_depths.txt!\n"), call.=FALSE)
if(length(depthseq) == 0)
depthseq <- seq(dmin, dmax, by=every)
dmin <- min(depthseq)
dmax <- max(depthseq)
if(proxies) {
storedat <- TRUE
if(file.exists(dd <- paste0(coredir, core, "/", core, "_proxies.csv")))
dat$proxies <- suppressWarnings(read.csv(dd, sep=sep)) else
stop(paste0("\nCannot find file ", dat$core, " _proxies.csv!\n"), call.=FALSE)
dmin <- min(depthseq, dat$proxies[,1])
dmax <- max(depthseq, dat$proxies[,1])
depthseq <- sort(unique(c(depthseq, dat$proxies[,1])))
}
if(length(ageofdepth) > 0)
depthseq <- sort(unique(c(ageofdepth, depthseq)))
# decide which models and point estimates should be used
if(any(type==c(1, "int", "inter", "interp"))) type <- 1 else
if(any(type==c(2, "reg", "regr", "poly", "polyn"))) type <- 2 else
if(any(type==c(3, "spline", "spl"))) type <- 3 else
if(any(type==c(4, "smooth", "sm"))) type <- 4 else
if(any(type==c(5, "loess", "lowess"))) type <- 5
best <- cbind(dat$mid1, dat$mid1, dat$mid1, dat$wmn, dat$med, dat$mode, dat$mid2)
Est <- best[,est]
# remove (manually chosen) outliers from the age-depth modelling
if(length(outliers) > 0) {
depths <- dat$depth[-outliers]
errors <- dat$error[-outliers]
calibs <- dat$calib[-outliers]
Est <- Est[-outliers]
} else {
depths <- dat$depth
errors <- dat$error
calibs <- dat$calib
}
#calibs <<- calibs # tmp
# age-depth modelling with curves through sampled age estimates
# in sections if one or more hiatuses are present
if(length(hiatus) > 0) {
allrange <- c(0,0,0,0)
hiatusseq <- sort(c(range(depthseq), hiatus))
for(i in 2:length(hiatusseq)) {
message(paste0("\n section ", i-1, ","))
section <- depthseq[min(which(depthseq >= hiatusseq[i-1])) : max(which(depthseq <= hiatusseq[i]))]
if(i>2) section <- section[-1]
sel <- min(which(depths >= min(section))):max(which(depths <= max(section)))
if(mixed.effect)
if(length(outliers) > 0)
smp <- .mixed.effect(its, depths, dat$cal[-outliers], dat$cage[-outliers], errors, calibs, Est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .mixed.effect(its, depths, dat$cal, dat$cage, errors, calibs, est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .smpl(its, depths[sel], calibs[sel], Est[sel])
calrange <- .model.clam(type, smooth, its, wghts, depths[sel], errors[sel], section, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
allrange <- rbind(allrange, calrange)
}
calrange <- allrange[2:nrow(allrange),]
} else {
if(mixed.effect) {
if(length(outliers) > 0)
smp <- .mixed.effect(its, depths, dat$cal[-outliers], dat$cage[-outliers], errors, calibs, est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .mixed.effect(its, depths, dat$cal, dat$cage, errors, calibs, Est, theta, f.mu, f.sigma, yrsteps, calibt)
} else
smp <- .smpl(its, depths, calibs, Est)
calrange <- .model.clam(type, smooth, its, wghts, depths, errors, depthseq, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
}
dat$model <- approx(calrange[,1], (calrange[,2]+calrange[,3])/2, dat$depth)$y
if(est==2) calrange[,4] <- (calrange[,2]+calrange[,3])/2
if(!BCAD && any(diff(calrange[,4]) < 0) || BCAD && any(diff(calrange[,4]) > 0))
reversal <- TRUE else reversal <- FALSE
gfit <- round(.gfit(theta, f.mu, f.sigma, dat, calrange, outliers), 2)
# re-correct the depths if slumps were applied
if(length(slump) > 0) {
dat <- .read.clam(core, coredir, ext, hpdsteps, yrsteps, prob, times, sep, BCAD, storedat, ignore, thickness, youngest, slump=NULL, threshold, theta, f.mu, f.sigma, calibt, extradates, calcurve, postbomb) # read in the original dates again
calrange <- calrange[which(calrange[,1] <= dmax),]
d <- calrange[,1]
for(i in 1:nrow(slump)) {
d[d > min(slump[i,])] <- d[d > min(slump[i,])] + (max(slump[i,]) - min(slump[i,]))
dmax <- dmax + (max(slump[i,]) - min(slump[i,]))
calrange[,1] <- d
hiatus[hiatus > min(slump[i,])] <- hiatus[hiatus > min(slump[i,])] + (max(slump[i,]) - min(slump[i,]))
}
}
# produce the age-depth plot, and a pdf copy if desired
if(length(yrmin)==0)
yrmin <- min(dat$mid1, calrange[,2])
if(length(yrmax)==0)
yrmax <- max(dat$mid1, calrange[,3])
if(length(ageofdepth > 0))
layout(matrix(c(1,2,1,3), nrow=2), heights=c(.7,.3))
.ageplot(yrmin, yrmax, dmin, dmax, revaxes, revd, revyr, yrlab, dlab, hiatus, depthseq, outliers, plotrange, BCAD, greyscale, if(length(greyscale)>0) get('chron') else NULL, C14col, outcol, outlsize, bestcol, rangecol, dat, calrange, depth, calhght, maxhght, mirror, calcol, slump, slumpcol, plotname, core, bty, mar, mgp, ash)
# write files providing calibrated dates, age-model and settings
colnames(calrange) <- c("Depth", paste0("min.", 100*prob, "%range"), paste0("max.", 100*prob, "%range"), "point")
.write.clam(dat, coredir, runname, calrange, core, prob, type, remove.reverse, smooth, wghts, its, outliers, ignore, est, BCAD, yrsteps, every, decimals, cmyr, depth, depthseq, hiatus, gfit, reversal, plotpdf, plotpng, yrmin, yrmax, dmin, dmax, dlab, yrlab, plotrange, greyscale, if(length(greyscale)>0) get('chron') else NULL, C14col, outcol, outlsize, bestcol, rangecol, calhght, maxhght, mirror, calcol, slump, slumpcol, revaxes, revyr, revd, calibt, youngest, extradates, plotname, calcurve, ccname, postbomb, pbnames, depths.file, bty, mar, mgp, ash)
if(storedat) {
calrange <<- calrange
dat <<- dat
smp <<- smp
}
# plot the age distribution of a provided depth
if(length(ageofdepth) > 0) {
if(revaxes)
abline(v=ageofdepth, lty=2) else
abline(h=ageofdepth, lty=2)
xlim <- range(.ageofdepth)
if(!BCAD) xlim <- xlim[2:1]
.ageofdepth <- get(".ageofdepth")
hst <- density(.ageofdepth, n=max(1, max(xlim)-min(xlim)))
yr <- seq(min(xlim), max(xlim), by=yrsteps)
hst <- cbind(c(yr, max(xlim), min(xlim)), c(approx(hst$x, hst$y, yr)$y, 0, 0))
plot(hst, type="n", main="", xlim=xlim, xlab=yrlab, ylab="")
polygon(hst, col="grey")
legend("topleft", paste(ageofdepth, depth), bty="n")
layout(matrix(1))
rng <- round(calrange[max(which(calrange[,1] <= ageofdepth)),])
message(" Age range of ", ageofdepth, " ", depth, ": ", rng[3], " to ", rng[2], ifelse(BCAD, " cal BC/AD", " cal BP"), " (", rng[3]-rng[2], " yr, ", prob, " % range)", sep="")
}
# report the confidence ranges, the goodness-of-fit, and whether any age-reversals occurred
rng <- round(calrange[,3]-calrange[,2])
message("\n", core, "'s ", 100*prob, "% confidence ranges span from ", min(rng), " to ", max(rng), " yr (average ", round(mean(rng)), " yr)")
message(" Fit (-log, lower is better): ", gfit)
if(reversal)
message(" Age reversals occurred. Try other model?\n")
}
# calculate the age-depth model and its uncertainty
.model.clam <- function(type, smooth, its, wghts, depths, errors, depthseq, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
{
# warn for extrapolation, refuse to do so for loess
if(min(depthseq) < min(dat$depth) || max(depthseq) > max(dat$depth))
if(type==5)
stop(" cannot extrapolate using loess! Change settings.\n ", call.=FALSE) else
cat(" extrapolating beyond dated levels, dangerous!\n ")
# choose model: interpolation, (polynomial) regression, spline, smooth spline or loess
chron <- array(0, dim=c(length(depthseq), its))
if(type==1) chron <- .interp(depthseq, depths, its, chron, smp) else
if(type==2) chron <- .poly(depthseq, smooth, wghts, errors, depths, its, chron, smp) else
if(type==3) chron <- .spline(depthseq, smooth, depths, its, chron, smp) else
if(type==4) chron <- .smooth(depthseq, smooth, wghts, errors, depths, its, chron, smp) else
if(type==5) chron <- .loess(depthseq, smooth, wghts, errors, depths, its, chron, smp)
# test against age reversals
# warp <- c()
warp <- NULL
if(remove.reverse!=FALSE)
for(i in 1:ncol(chron))
if(!BCAD && min(diff(chron[,i])) <= 0 || BCAD && max(diff(chron[,i])) >= 0)
warp <- c(warp, i)
if(length(warp) > 0)
if(length(warp) > remove.reverse*its)
message("\n !!! Too many models with age reversals!!!\n") else
{
message("\n Removing ", length(warp), " models with age reversals, ", its-length(warp), " models left...")
chron <- chron[,-warp]
smp <- smp[,-warp,]
}
if(length(ageofdepth) > 0)
if(ageofdepth %in% depthseq)
.assign_to_global(".ageofdepth", chron[which(depthseq==ageofdepth),])
if(storedat)
{
chron <<- chron
smp <<- smp
}
# find uncertainty ranges of calendar age for each depth of the core
calrange <- array(0, dim=c(nrow(chron), 2))
# mn <- c()
mn <- numeric(nrow(chron))
for(i in 1:nrow(chron))
{
x <- chron[i,2:ncol(chron)]
qp <- (1-prob)/2
calrange[i,] <- quantile(x, c(qp, 1-qp))
if(est==1) mn[i] <- mean(x)
}
if(est==1)
cbind(depthseq, cbind(calrange, mn)) else
cbind(depthseq, cbind(calrange, chron[,1]))
}
Maarten14C/clam documentation built on Aug. 18, 2022, 6:33 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .model.clam clam
Documented in clam
#' clam
#'
#' Classical (non-Bayesian) age-depth modelling.
#'
#' @docType package
#' @author Maarten Blaauw <maarten.blaauw@qub.ac.uk>
#' @importFrom grDevices dev.off grey pdf png rgb
#' @importFrom graphics abline image layout legend lines par points polygon rect plot
#' @importFrom stats approx density dnorm lm loess pnorm predict qnorm quantile rnorm runif smooth.spline spline weighted.mean
#' @importFrom utils packageName read.csv read.table write.table
#' @import rintcal
#' @name clam
NULL
# do:
# done:
# so that functions such as pMC.age etc. are available immediately
library(rintcal)
#' @name clam
#' @title The main age-depth modelling function
#' @description Produce age-depth models for cores with dated depths.
#' @details
#' Cores containing several 14C and/or other dates can be processed semi-automatically in order to obtain age-depth models.
#' In the process, any 14C dates are calibrated, and age-depth curves are repeatedly drawn through point estimates sampled from the dates.
#' Age-depth models can be based on linear interpolation, linear/polynomial regression, or cubic, smooth or locally weighted splines.
#' For each date, the probability of a calendar year being sampled is proportionate to its calibrated probability (see Blaauw, 2010).
#' Uncertainty ranges as well as a 'best' age-model are calculated.
#'
#' Additional cores should be put in a comma-separated file in a sub-folder of the directory where the cores are stored.
#' By default this parent folder is called \code{coredir="clam_runs"} (if no folder called \code{"Cores"} already exists). If your core is called MyCore1, save MyCore1.csv as \code{clam_runs/MyCore1/MyCore1.csv}.
#' Ensure that the names of the core's folder and filename's root (the part before .csv) match, e.g., using exactly similar upper- and lower case letters.
#'
#' Avoid the use of spaces or non-standard (non-ASCII) characters within the file or in folder or file names.
#' The plain text file should consist of 6 or 7 columns (also called fields), containing in the following exact order (see the example below):
#' \enumerate{
#' \item Identification labels (e.g. 14C lab codes)
#' \item 14C ages for 14C-dated depths; leave empty for non-14C dated depths
#' \item cal BP ages (for any non-14C dates such as the core surface; leave empty for levels with 14C dates)
#' \item errors (reported 1 standard deviation errors. This column should never be left empty. Errors should always be larger than 0)
#' \item age offsets if known (otherwise leave empty)
#' \item depths (depths in the sequence were the dated samples were taken, default unit depth="cm"; this column should never be left empty)
#' \item thicknesses of the sampled slices (optional column; leave empty for default of 1)
#' }
#' Add a final empty line to your core's .csv file by pressing 'Enter' after the file's last value.
#'
#' These files can be made in spreadsheet software such as MS-Excel, but it is always a good idea to check the file's formatting in a plain-text editor such as WordPad. Remove any lines which contain only commas, and it is also recommended to remove quotes ()\code{\" or \'}) in the headers or elsewhere.
#'
#' Age-models for the core can then be produced by typing, e.g., \code{clam("MyCore1")}.
#'
#' By default the northern hemisphere terrestrial calibration curve is used (\code{cc=1}, \code{cc1="3Col_intcal20.14C"}).
#' To use alternative curves, change \code{cc} to \code{cc=2 (cc2="3Col_marine20.14C")}, \code{cc=3 (cc3="3Col_shcal20.14C")}, \code{cc=4 (cc4="mixed.14C")}.
#' You can also provide custom-built calibration curves, indicating its location using \code{ccdir}.
#'
#' The provided example (default \code{core="Example"}) is core Quilichao-1 which was sampled from a Colombian lake (Berrio et al., 2002).
#' This core was chosen because it was dated at a rather high resolution, and appears to contain a hiatus (e.g., try \code{hiatus=450}
#' for a hiatus at 450 cm depth).
#'
#' Each clam run will produce a range of files within the core's folder. One, ending with \code{"_calibrated.txt"} contains the calibrated
#' age ranges of the 14C and other dates. The others will be named according to the core's name followed by the model type,
#' and contain the age estimates for all depths (files ending with \code{"_ages.txt"}), settings (files ending with \code{"_settings.txt"})
#' and graphs (files ending with \code{".pdf"} and \code{".png"}).
#' The file containing the age estimates has 5 columns; first the depths, then the minima and maxima of the confidence intervals,
#' then a "best" estimate, and finally the reconstructed accumulation rates. The reported values are rounded to 0 decimals by default
#' (\code{decimals=0}). Accumulation rates are in yr/cm ("deposition time") by default (\code{cmyr=FALSE}), but can be reported in cm/yr (\code{cmyr=TRUE}).
#'
#' see Blaauw 2010 (Quaternary Geochronology 5: 512-518).
#'
#' @param core Name of the core, given using quotes. Defaults to the core provided with clam, \code{core="Example"}.
#' @param type The type of age-depth model. Five different types are provided:
#' \enumerate{
#' \item linear interpolation between neighbouring levels (1, "int", "inter" or "interp")
#' \item linear or higher polynomial regression (2, "reg", "regr", "poly" or "polyn", default linear)
#' \item cubic spline (3, "spl" or "spline")
#' \item smooth spline (4, "sm" or "smooth", default smoothing 0.3)
#' \item locally weighted spline (5, "loess" or "lowess", default smoothing 0.75, cannot extrapolate)
#' }
#' @param smooth Degree of smoothing. Gives polynomial degree for model type 2. Not relevant for \code{type=1 or type=3}.
#' \itemize{
#' \item for type=2: \code{smooth=1} (linear), \code{smooth=2} second-order polynomial, \code{smooth=3} for third-order polynomial, etc.
#' \item for type=4: \code{smooth=0.3}
#' \item for type=5: \code{smooth=0.75}
#' }
#' @param prob Confidence intervals (between 0 and 1), default \code{prob=0.95} or 95\%.
#' @param its Amount of age-model iterations; defaults to \code{its=1000}.
#' @param coredir The directory where core runs are stored (each core in its own directory named after the core's name).
#' @param ask By default, and as per R rules, clam will ask if it is OK to make or write to a directory.
#' Defaults to \code{coredir="clam_runs"}, or to \code{coredir="Cores"} if this folder exists where R is working.
#' @param wghts Weights can be applied to dated depths as follows:
#' \itemize{
#' \item 0 no weighting
#' \item 1 weighted to calibrated probabilities of sampled calendar years (default, \code{wghts=1}).
#' \item 2 weighted to (inverse squared) errors of the dates.
#' }
#' @param cc calibration curve for C14 dates (1, 2 or 3).
#' @param cc1 For terrestrial, northern hemisphere C14 dates.
#' @param cc2 For marine C14 dates.
#' @param cc3 For southern hemisphere C14 dates.
#' @param cc4 For mixed terrestrial/marine C14 dates.
#' @param postbomb Use a postbomb curve for negative (i.e. postbomb) 14C ages. \code{0 = none, 1 = NH1, 2 = NH2, 3 = NH3, 4 = SH1-2, 5 = SH3}. See \url{http://calib.org/CALIBomb/}.
#' @param pb1 For Northern hemisphere region 1 postbomb C14 dates.
#' @param pb2 For Northern hemisphere region 2 postbomb C14 dates.
#' @param pb3 For Northern hemisphere region 3 postbomb C14 dates.
#' @param pb4 For Southern hemisphere regions 1-2 postbomb C14 dates.
#' @param pb5 For Southern hemisphere region 3 postbomb C14 dates.
#' @param ccdir Directory where the calibration curves for C14 dates \code{cc} are located. By default \code{ccdir=""}.
#' For example, use \code{ccdir="."} to choose current working directory, or \code{ccdir="Curves/"} to choose sub-folder \code{Curves/}.
#' @param outliers The number of any dates to be considered outlying, e.g. \code{c(5,6)} for the fifth and sixth dated depth counting from the top of a core.
#' @param ignore The number of any dates that should be ignored, e.g., \code{c(5,6)} for the fifth and sixth date counting from the top of a core.
#' @param youngest The age beyond which dates should be truncated (e.g., \code{youngest=-60} if the core was sampled in -60 cal BP or AD 2010).
#' @param extradates Depths of any additional dates with their files of ages and probabilities.
#' @param slump Upper and lower depths of sections of abrupt accumulation that should be excised, e.g., \code{c(600, 550, 120, 100)} for two sections of 600-550 and 120-100 cm depth.
#' @param est Which point estimate to use as 'best' age. It is highly recommended to not only use these 'best' point estimates, as chronological uncertainties are often considerable and should not be ignored.
#' \enumerate{
#' \item averages of age-depth model derived ages (default, \code{est=1})
#' \item midpoints of age-depth model derived age estimates
#' \item midpoints of calibrated ranges
#' \item weighted means of calibrated ranges
#' \item medians of calibrated distributions
#' \item maximum densities of calibrated distributions
#' \item midpoints of entire calibrated distributions (including years outside the calibrated ranges)
#' }
#' @param calibt Calibration based on the student-t distribution. By default, the Gaussian distribution is used (\code{calibt=FALSE}). To use the student-t distribution, provide two parameters such as \code{calibt=c(3,4)}.
#' @param mixed.effect Set to \code{TRUE} to activate mixed-effect modelling.
#' @param dmin Minimum depth of age-depth model (e.g., extrapolate).
#' @param dmax Maximum depth of age-depth model (e.g., extrapolate).
#' @param every Resolution at which (ages for) depths are calculated.
#' @param yrmin Minimum of calendar axis of age-depth plot (calculate automatically by default).
#' @param yrmax Maximum of calendar axis of age-depth plot (calculated automatically by default).
#' @param yrsteps Temporal resolution at which calibrated ages are calculated (in calendar years).
#' @param pbsteps Temporal resolution at which postbomb C14 ages are calibrated (in calendar years).
#' @param hpdsteps Temporal resolution at which highest posterior density ranges are calibrated (in calendar years).
#' @param BCAD Use BC/AD or cal BP scale.
#' @param decimals Amount of decimals for rounding.
#' @param cmyr Accumulation rates can be provided as yr/cm (default, \code{cmyr=TRUE}, more accurately named deposition times) or cm/yr (\code{cmyr=FALSE}).
#' @param ageofdepth Calculate age estimates of a specific depth.
#' @param depth Depth units.
#' @param depthseq Sequence of depths for which age estimates are to be calculated (default: from \code{dmin} to \code{dmax} with steps of size every)
#' @param depths.file Use a file with depths for depthseq.
#' @param thickness Thickness of the dated samples.
#' @param hiatus Depths of any hiatuses, e.g., \code{c(500, 300)}. Each sub-section must have at least 2 dates (\code{4} for smoothing spline; does not work with loess as it cannot extrapolate).
#' @param remove.reverse Proportion of age-models with reversals that can be removed before prompting a warning.
#' Set at \code{FALSE} to avoid removing models with reversals.
#' @param times Half-range of calibration curve used to calibrate dates (multiplication factor for the dates' errors).
#' @param sep Separator between the fields of the plain text file containing the dating information.
#' @param ext Extension of the file containing the dating information.
#' @param runname Text to add to the core name for specific runs, e.g., "MyCore_Test1"
#' @param storedat Store the dates and age-model within R after a \code{clam} run. Defaults to \code{storedat=TRUE}.
#' @param threshold Below which value should probabilities be excluded from calculations.
#' @param proxies Set to \code{TRUE} to plot proxies against age after the run.
#' @param revaxes Set to \code{TRUE} to plot ages on the vertical axis and depth on the horizontal axis.
#' @param revd Plot depth axis in reverse.
#' @param revyr Plot age axis in reverse.
#' @param calhght Heights of the calibrated distributions in the age-depth plot.
#' @param maxhght Maximum height of age probability distributions.
#' @param mirror Plot the age distributions in "mirror" style (above and below depth).
#' @param plotrange Plot the confidence ranges of the age-model.
#' @param mar Plot margins (amount of white space along edges of axes 1-4).
#' @param mgp Axis text margins (where should titles, labels and tick marks be plotted).
#' @param bty Type of box to be drawn around plots. Draw a box around the graph (\code{"n"} for none,
#' and \code{"l"}, \code{"7"}, \code{"c"}, \code{"u"}, "]" or \code{"o"} for correspondingly shaped boxes).
#' @param plotpdf Produce a pdf file of the age-depth plot.
#' @param plotpng Produce a png file of the age-depth plot.
#' @param greyscale Produce a grey-scale representation of all age-models (number gives resolution, e.g., 500 bins; will cancel plotting of the confidence intervals).
#' @param yrlab Label of the calendar axis. Defaults to either cal BP or BC/AD. Alternative names can be provided.
#' @param dlab Label of the depth axis. Defaults to \code{dlab="Depth (cm)"} (assuming \code{depth="cm"}), but alternative names can be provided.
#' @param calcol Colour of the calibrated distributions in the age-depth plot.
#' @param C14col Colour of the calibrated ranges of the dates.
#' @param outcol Colour of outlying dates.
#' @param outlsize Size of symbols outlying dates.
#' @param bestcol Colour of the "best" age-depth model (based on chosen value for est).
#' @param rangecol Colour of plotted confidence ranges.
#' @param slumpcol Colour of slump.
#' @param plotname Print the core name on the graph.
#' @param ash Plot all distributions at the same height.
#' @param rule How should R's approx function deal with extrapolation. If \code{rule=1}, the default, then NAs are returned for such points and if it is 2, the value at the closest data extreme is used.
#'
#' @author Maarten Blaauw
#' @return Age model construction together with a text output and files saved to a folder in the \code{coredir/core} directory.
#' @examples
#' clam(, coredir=tempdir()) # Create the example in Cores/Example folder
#' clam(, coredir=tempdir(), extradates=470)
#'
#' @references
#' Berrio, J.C., Hooghiemstra, H., Marchant, R., Rangel, O., 2002. Late-glacial and Holocene history of the dry forest area
#' in the south Colombian Cauca Valley. Journal of Quaternary Science 17, 667-682
#'
#' Blaauw, M., 2010. Methods and code for 'classical' age-modelling of radiocarbon sequences. Quaternary Geochronology 5, 512-518
#' \doi{10.1016/j.quageo.2010.01.002}
#' @export
clam <- function(core="Example", type=1, smooth=NULL, prob=0.95, its=1000, coredir=NULL, ask=TRUE, wghts=1, cc=1, cc1="3Col_intcal20.14C", cc2="3Col_marine20.14C", cc3="3Col_shcal20.14C", cc4="mixed.14C", postbomb=FALSE, pb1="postbomb_NH1.14C", pb2="postbomb_NH2.14C", pb3="postbomb_NH3.14C", pb4="postbomb_SH1-2.14C", pb5="postbomb_SH3.14C", ccdir="", outliers=NULL, ignore=NULL, youngest=NULL, extradates=NULL, slump=NULL, est=1, calibt=FALSE, mixed.effect=FALSE, dmin=NULL, dmax=NULL, every=1, yrmin=NULL, yrmax=NULL, yrsteps=1, pbsteps=0.01, hpdsteps=1, BCAD=FALSE, decimals=0, cmyr=FALSE, ageofdepth=NULL, depth="cm", depthseq=NULL, depths.file=FALSE, thickness=1, hiatus=NULL, remove.reverse=0.5, times=5, sep=",", ext=".csv", runname=NULL, storedat=TRUE, threshold=1e-6, proxies=FALSE, revaxes=FALSE, revd=TRUE, revyr=TRUE, calhght=0.3, maxhght=0.01, mirror=TRUE, plotrange=TRUE, bty="l", mar=c(3.5,3,2,1), mgp=c(2,1,0), plotpdf=TRUE, plotpng=TRUE, greyscale=NULL, yrlab=NULL, dlab=NULL, calcol=rgb(0,0.5,0.5,0.5), C14col=rgb(0,0,1,0.5), outcol="red", outlsize=1, bestcol="black", rangecol=rgb(0,0,0,0.3), slumpcol=grey(0.75), plotname=TRUE, ash=FALSE, rule=1) {
# If coredir is left empty, check for a folder named Cores in the current working directory, and if this doesn't exist, for a folder called clam_runs (make this folder if it doesn't exist yet).
# Check if we have write access. If not, tell the user to provide a different, writeable location for coredir.
coredir <- assign_coredir(coredir, core, ask)
if(core == "Example") {
dir.create(file.path(coredir, "Example/"), showWarnings = FALSE, recursive = TRUE)
fileCopy <- system.file("extdata/Example/", package="clam")
file.copy(fileCopy, coredir, recursive = TRUE, overwrite=FALSE)
}
# warn and stop if unexpected settings are provided
if(type > 5 || type < 1 || prob < 0 || prob > 1 || its < 100 || wghts < 0 || wghts > 1 || est < 1 || est > 7 || yrsteps <= 0 || hpdsteps <= 0 || every <= 0 || decimals < 0 || cmyr < 0 || cmyr > 1 || thickness < 0 || times < 1 || calhght < 0 || (type==5 && length(hiatus)>0))
stop("\n Warning, clam cannot run with these settings! Please check the manual.\n\n", call.=FALSE)
csvFile <- paste0(coredir, core, "/", core, ext)
if(!file.exists(csvFile))
stop(paste("\nInput data file", csvFile, "not found."), call.=FALSE)
dets <- read.csv(csvFile, sep=sep)
d <- dets[,6]
these <- 1:length(d) # to possibly kick out dates later on; Nov 2020
if(min(diff(d)) < 0)
message(" Warning, depths not in ascending order (top ones should come first).")
# avoid Windows/Mac habit of silently adding .txt extension to plain text files
bill <- list.files(file.path(coredir, core), pattern=".csv.txt")
if(length(bill) > 0) {
message("Removing unnecessary .txt extension from .csv file", bill[1])
file.rename(file.path(coredir, core, paste0(core, ".csv.txt")),
file.path(coredir, core, "/", paste0(core, ".csv")))
}
# set the calibration curve
ccdir <- .validateDirectoryName(ccdir)
if(ccdir == "") # so, if no alternative folder provided, use clam's calibration curves
ccdir = system.file("extdata", package="rintcal")
if(cc==1) calcurve <- read.table(file.path(ccdir, cc1)) else
if(cc==2) calcurve <- read.table(file.path(ccdir, cc2)) else
if(cc==3) calcurve <- read.table(file.path(ccdir, cc3)) else
if(cc==4) calcurve <- read.table(file.path(ccdir, cc4)) else
stop("I do not understand which calibration curve you mean, check the manual", call.=FALSE)
if(cc==1) ccname <- cc1 else
if(cc==2) ccname <- cc2 else
if(cc==3) ccname <- cc3 else
if(cc==4) ccname <- cc4
# negative C14 ages need a postbomb curve
pbnames <- c(pb1, pb2, pb3, pb4, pb5)
cdat <- dets[,2]
if(length(cdat[!is.na(cdat)]) > 0)
if(min(cdat[!is.na(cdat)]) < 0)
if(postbomb==FALSE)
message("Warning, negative 14C ages, should I use a postbomb curve?") else {
if(postbomb>5)
stop("I do not understand which postbomb curve you mean, check the manual", call.=FALSE)
yrsteps <- min(pbsteps, yrsteps)
pb <- read.table(system.file("extdata", pbnames[postbomb], package="rintcal")) # now points to rintcal package
pb.x <- seq(min(pb[,1]), max(pb[,1]), by=yrsteps)
pb.y <- approx(pb[,1], pb[,2], pb.x)$y
pb.sd <- approx(pb[,1], pb[,3], pb.x)$y
calcurve <- cbind(c(pb.x, calcurve[,1]), c(pb.y, calcurve[,2]), c(pb.sd, calcurve[,3]))
}
# work in BC/AD if needed, and prepare for calculations in f14C
if(BCAD) {
theta <- 1950-calcurve[,1]
border <- max(which(theta > 0))
theta <- c(theta[1:(border-1)], theta[border]:theta[border+2], theta[(border+3):length(theta)])
mu <- approx(1950-calcurve[,1], calcurve[,2], theta, rule=rule)$y
sigma <- approx(1950-calcurve[,1], calcurve[,3], theta, rule=rule)$y
theta[theta <=0] <- theta[theta <= 0]-1
calcurve <- cbind(theta, mu, sigma)
} else theta <- calcurve[,1]
if(length(yrlab) == 0)
yrlab <- ifelse(BCAD, "cal BC/AD", "cal BP")
f.mu <- exp(-calcurve[,2]/8033)
f.sigma <- exp(-(calcurve[,2]-calcurve[,3])/8033) - f.mu
# prepare for slumps and hiatuses
#cat(paste("Core name Test:", core))
if(length(greyscale) > 0) storedat <- TRUE
if(length(slump) > 0) {
if(length(slump) %% 2 == 1)
stop("\n Warning, slumps need both upper and lower depths. Please check the manual", call.=FALSE)
slump <- matrix(sort(slump), ncol=2, byrow=TRUE)
if(length(dmax)==0)
dmax <- max(dets[,6])
if(length(extradates) > 0)
dmax <- max(dmax, extradates)
for(i in 1:nrow(slump)) {
d[d > min(slump[i,])] <- d[d > min(slump[i,])] - (max(slump[i,]) - min(slump[i,]))
dmax <- dmax - (max(slump[i,])-min(slump[i,]))
}
if(length(hiatus) > 0)
for(i in 1:nrow(slump)) {
below.slump <- which(hiatus > max(slump[i,]))
above.slump <- which(hiatus < min(slump[i,]))
hiatus[below.slump] <- hiatus[below.slump] - (max(slump[i,])-min(slump[i,]))
hiatus <- hiatus[c(above.slump, below.slump)]
}
}
# read in the data
dat <- .read.clam(core, coredir, ext, hpdsteps, yrsteps, prob, times, sep, BCAD, storedat, ignore, thickness, youngest, slump, threshold, theta, f.mu, f.sigma, calibt, extradates, calcurve, postbomb, rule=rule)
message(" Calibrating dates... ")
# calculate the depths to be used, based on the ranges and resolution
if(length(dmin)==0)
dmin <- floor(min(dat$depth))
if(length(dmax)==0)
dmax <- ceiling(max(dat$depth))
if(depths.file)
if(file.exists(dd <- paste0(coredir, core, "/", core, "_depths.txt"))) {
if(length(depthseq) == 0)
depthseq <- seq(dmin, dmax, by=every)
depthseq <- sort(unique(c(depthseq, suppressWarnings(read.table(dd))[,1])))
# dmin <- min(depthseq)#, read.table(dd)[,1])
# dmax <- max(depthseq)#, read.table(dd)[,1])
} else
stop(paste0("\nCannot find file ", dat$core, "_depths.txt!\n"), call.=FALSE)
if(length(depthseq) == 0)
depthseq <- seq(dmin, dmax, by=every)
dmin <- min(depthseq)
dmax <- max(depthseq)
if(proxies) {
storedat <- TRUE
if(file.exists(dd <- paste0(coredir, core, "/", core, "_proxies.csv")))
dat$proxies <- suppressWarnings(read.csv(dd, sep=sep)) else
stop(paste0("\nCannot find file ", dat$core, " _proxies.csv!\n"), call.=FALSE)
dmin <- min(depthseq, dat$proxies[,1])
dmax <- max(depthseq, dat$proxies[,1])
depthseq <- sort(unique(c(depthseq, dat$proxies[,1])))
}
if(length(ageofdepth) > 0)
depthseq <- sort(unique(c(ageofdepth, depthseq)))
# decide which models and point estimates should be used
if(any(type==c(1, "int", "inter", "interp"))) type <- 1 else
if(any(type==c(2, "reg", "regr", "poly", "polyn"))) type <- 2 else
if(any(type==c(3, "spline", "spl"))) type <- 3 else
if(any(type==c(4, "smooth", "sm"))) type <- 4 else
if(any(type==c(5, "loess", "lowess"))) type <- 5
best <- cbind(dat$mid1, dat$mid1, dat$mid1, dat$wmn, dat$med, dat$mode, dat$mid2)
Est <- best[,est]
# remove (manually chosen) outliers from the age-depth modelling
if(length(outliers) > 0) {
depths <- dat$depth[-outliers]
errors <- dat$error[-outliers]
calibs <- dat$calib[-outliers]
Est <- Est[-outliers]
} else {
depths <- dat$depth
errors <- dat$error
calibs <- dat$calib
}
#calibs <<- calibs # tmp
# age-depth modelling with curves through sampled age estimates
# in sections if one or more hiatuses are present
if(length(hiatus) > 0) {
allrange <- c(0,0,0,0)
hiatusseq <- sort(c(range(depthseq), hiatus))
for(i in 2:length(hiatusseq)) {
message(paste0("\n section ", i-1, ","))
section <- depthseq[min(which(depthseq >= hiatusseq[i-1])) : max(which(depthseq <= hiatusseq[i]))]
if(i>2) section <- section[-1]
sel <- min(which(depths >= min(section))):max(which(depths <= max(section)))
if(mixed.effect)
if(length(outliers) > 0)
smp <- .mixed.effect(its, depths, dat$cal[-outliers], dat$cage[-outliers], errors, calibs, Est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .mixed.effect(its, depths, dat$cal, dat$cage, errors, calibs, est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .smpl(its, depths[sel], calibs[sel], Est[sel])
calrange <- .model.clam(type, smooth, its, wghts, depths[sel], errors[sel], section, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
allrange <- rbind(allrange, calrange)
}
calrange <- allrange[2:nrow(allrange),]
} else {
if(mixed.effect) {
if(length(outliers) > 0)
smp <- .mixed.effect(its, depths, dat$cal[-outliers], dat$cage[-outliers], errors, calibs, est, theta, f.mu, f.sigma, yrsteps, calibt) else
smp <- .mixed.effect(its, depths, dat$cal, dat$cage, errors, calibs, Est, theta, f.mu, f.sigma, yrsteps, calibt)
} else
smp <- .smpl(its, depths, calibs, Est)
calrange <- .model.clam(type, smooth, its, wghts, depths, errors, depthseq, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
}
dat$model <- approx(calrange[,1], (calrange[,2]+calrange[,3])/2, dat$depth)$y
if(est==2) calrange[,4] <- (calrange[,2]+calrange[,3])/2
if(!BCAD && any(diff(calrange[,4]) < 0) || BCAD && any(diff(calrange[,4]) > 0))
reversal <- TRUE else reversal <- FALSE
gfit <- round(.gfit(theta, f.mu, f.sigma, dat, calrange, outliers), 2)
# re-correct the depths if slumps were applied
if(length(slump) > 0) {
dat <- .read.clam(core, coredir, ext, hpdsteps, yrsteps, prob, times, sep, BCAD, storedat, ignore, thickness, youngest, slump=NULL, threshold, theta, f.mu, f.sigma, calibt, extradates, calcurve, postbomb) # read in the original dates again
calrange <- calrange[which(calrange[,1] <= dmax),]
d <- calrange[,1]
for(i in 1:nrow(slump)) {
d[d > min(slump[i,])] <- d[d > min(slump[i,])] + (max(slump[i,]) - min(slump[i,]))
dmax <- dmax + (max(slump[i,]) - min(slump[i,]))
calrange[,1] <- d
hiatus[hiatus > min(slump[i,])] <- hiatus[hiatus > min(slump[i,])] + (max(slump[i,]) - min(slump[i,]))
}
}
# produce the age-depth plot, and a pdf copy if desired
if(length(yrmin)==0)
yrmin <- min(dat$mid1, calrange[,2])
if(length(yrmax)==0)
yrmax <- max(dat$mid1, calrange[,3])
if(length(ageofdepth > 0))
layout(matrix(c(1,2,1,3), nrow=2), heights=c(.7,.3))
.ageplot(yrmin, yrmax, dmin, dmax, revaxes, revd, revyr, yrlab, dlab, hiatus, depthseq, outliers, plotrange, BCAD, greyscale, if(length(greyscale)>0) get('chron') else NULL, C14col, outcol, outlsize, bestcol, rangecol, dat, calrange, depth, calhght, maxhght, mirror, calcol, slump, slumpcol, plotname, core, bty, mar, mgp, ash)
# write files providing calibrated dates, age-model and settings
colnames(calrange) <- c("Depth", paste0("min.", 100*prob, "%range"), paste0("max.", 100*prob, "%range"), "point")
.write.clam(dat, coredir, runname, calrange, core, prob, type, remove.reverse, smooth, wghts, its, outliers, ignore, est, BCAD, yrsteps, every, decimals, cmyr, depth, depthseq, hiatus, gfit, reversal, plotpdf, plotpng, yrmin, yrmax, dmin, dmax, dlab, yrlab, plotrange, greyscale, if(length(greyscale)>0) get('chron') else NULL, C14col, outcol, outlsize, bestcol, rangecol, calhght, maxhght, mirror, calcol, slump, slumpcol, revaxes, revyr, revd, calibt, youngest, extradates, plotname, calcurve, ccname, postbomb, pbnames, depths.file, bty, mar, mgp, ash)
if(storedat) {
calrange <<- calrange
dat <<- dat
smp <<- smp
}
# plot the age distribution of a provided depth
if(length(ageofdepth) > 0) {
if(revaxes)
abline(v=ageofdepth, lty=2) else
abline(h=ageofdepth, lty=2)
xlim <- range(.ageofdepth)
if(!BCAD) xlim <- xlim[2:1]
.ageofdepth <- get(".ageofdepth")
hst <- density(.ageofdepth, n=max(1, max(xlim)-min(xlim)))
yr <- seq(min(xlim), max(xlim), by=yrsteps)
hst <- cbind(c(yr, max(xlim), min(xlim)), c(approx(hst$x, hst$y, yr)$y, 0, 0))
plot(hst, type="n", main="", xlim=xlim, xlab=yrlab, ylab="")
polygon(hst, col="grey")
legend("topleft", paste(ageofdepth, depth), bty="n")
layout(matrix(1))
rng <- round(calrange[max(which(calrange[,1] <= ageofdepth)),])
message(" Age range of ", ageofdepth, " ", depth, ": ", rng[3], " to ", rng[2], ifelse(BCAD, " cal BC/AD", " cal BP"), " (", rng[3]-rng[2], " yr, ", prob, " % range)", sep="")
}
# report the confidence ranges, the goodness-of-fit, and whether any age-reversals occurred
rng <- round(calrange[,3]-calrange[,2])
message("\n", core, "'s ", 100*prob, "% confidence ranges span from ", min(rng), " to ", max(rng), " yr (average ", round(mean(rng)), " yr)")
message(" Fit (-log, lower is better): ", gfit)
if(reversal)
message(" Age reversals occurred. Try other model?\n")
}
# calculate the age-depth model and its uncertainty
.model.clam <- function(type, smooth, its, wghts, depths, errors, depthseq, prob, est, dat, smp, greyscale, remove.reverse, storedat, ageofdepth, BCAD)
{
# warn for extrapolation, refuse to do so for loess
if(min(depthseq) < min(dat$depth) || max(depthseq) > max(dat$depth))
if(type==5)
stop(" cannot extrapolate using loess! Change settings.\n ", call.=FALSE) else
cat(" extrapolating beyond dated levels, dangerous!\n ")
# choose model: interpolation, (polynomial) regression, spline, smooth spline or loess
chron <- array(0, dim=c(length(depthseq), its))
if(type==1) chron <- .interp(depthseq, depths, its, chron, smp) else
if(type==2) chron <- .poly(depthseq, smooth, wghts, errors, depths, its, chron, smp) else
if(type==3) chron <- .spline(depthseq, smooth, depths, its, chron, smp) else
if(type==4) chron <- .smooth(depthseq, smooth, wghts, errors, depths, its, chron, smp) else
if(type==5) chron <- .loess(depthseq, smooth, wghts, errors, depths, its, chron, smp)
# test against age reversals
# warp <- c()
warp <- NULL
if(remove.reverse!=FALSE)
for(i in 1:ncol(chron))
if(!BCAD && min(diff(chron[,i])) <= 0 || BCAD && max(diff(chron[,i])) >= 0)
warp <- c(warp, i)
if(length(warp) > 0)
if(length(warp) > remove.reverse*its)
message("\n !!! Too many models with age reversals!!!\n") else
{
message("\n Removing ", length(warp), " models with age reversals, ", its-length(warp), " models left...")
chron <- chron[,-warp]
smp <- smp[,-warp,]
}
if(length(ageofdepth) > 0)
if(ageofdepth %in% depthseq)
.assign_to_global(".ageofdepth", chron[which(depthseq==ageofdepth),])
if(storedat)
{
chron <<- chron
smp <<- smp
}
# find uncertainty ranges of calendar age for each depth of the core
calrange <- array(0, dim=c(nrow(chron), 2))
# mn <- c()
mn <- numeric(nrow(chron))
for(i in 1:nrow(chron))
{
x <- chron[i,2:ncol(chron)]
qp <- (1-prob)/2
calrange[i,] <- quantile(x, c(qp, 1-qp))
if(est==1) mn[i] <- mean(x)
}
if(est==1)
cbind(depthseq, cbind(calrange, mn)) else
cbind(depthseq, cbind(calrange, chron[,1]))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.