Skip to contents

Bayesian analysis for age estimation of OSL measurements and C-14 ages of various samples

Source: R/Age_OSLC14.R
Age_OSLC14.Rd

This function computes an age of OSL data consisting of at least two samples and calibrate C-14 ages of samples to get an age (in ka).
Ages of OSL data are computed according to the model given in Combes and Philippe (2017). Single-grain or Multi-grain OSL measurements can be analysed simultaneously (with output of Generate_DataFile or Generate_DataFile_MG or both of them using combine_DataFiles). Samples, for which data is available in several BIN files, can be analysed.
For C-14 data, the user can choose one of the following radiocarbon calibration curve: Northern or Southern Hemisphere or marine atmospheric.

Usage

Age_OSLC14(
  DATA,
  Data_C14Cal,
  Data_SigmaC14Cal,
  Nb_sample = DATA$Nb_sample,
  SampleNames = DATA$SampleNames,
  SampleNature,
  PriorAge = rep(c(10, 60), Nb_sample),
  SavePdf = FALSE,
  OutputFileName = c("MCMCplot", "HPD_Cal14CCurve", "summary"),
  OutputFilePath = c(""),
  SaveEstimates = FALSE,
  OutputTableName = c("DATA"),
  OutputTablePath = c(""),
  StratiConstraints = c(),
  sepSC = c(","),
  BinPerSample = rep(1, sum(SampleNature[1, ])),
  THETA = c(),
  sepTHETA = c(","),
  LIN_fit = TRUE,
  Origin_fit = FALSE,
  distribution = c("cauchy"),
  Model_C14 = c("full"),
  CalibrationCurve = c("IntCal20"),
  Iter = 10000,
  burnin = 4000,
  adapt = 1000,
  t = 5,
  n.chains = 3,
  jags_method = "rjags",
  autorun = FALSE,
  quiet = FALSE,
  roundingOfValue = 3,
  ...
)

Arguments

DATA

Two types of inputs are supported. (1): a list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J, K, Nb_measurement, provided by the function Generate_DataFile, Generate_DataFile_MG or combine_DataFiles. DATA contains information for more than one sample. If there is stratigraphic relations between samples, informations in DATA must be ordered by order of increasing ages. See the details section to for more informations. (2): an object of class BayLum.list which is provided by the output of Age_OSLC14. When input of class BayLum.list is identified, no new JAGS model is created. Instead, the JAGS model specified within the BayLum.list is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.

Data_C14Cal

numeric vector: corresponding to C-14 age estimate (in years, conversion in ka is automatically done in the function). If there is stratigraphic relations between samples, Data_C14Cal must be ordered by order of increasing ages.

Data_SigmaC14Cal

numeric: corresponding to the error of C-14 age estimates.

Nb_sample

numeric: number of samples (OSL data and C-14 age), (Nb_sample>3, at least to sample of OSL data and one sample of C-14 age).

SampleNames

character: sample names for both OSL data and C14 data. The length of this vector is equal to Nb_sample. If there is stratigraphic relation, this vector must be ordered by increasing order (to mix OSL samples and C-14 ages if it is needed).

SampleNature

matrix: states the nature of the sample. Row number of SampleNature matrix is equal to 2 and column number is equal to Nb_sample. First line of the matrix: SampleNature[1,i] states if sample whose number ID is equal to i, is an OSL sample 1 or not 0. Second line of the matrix: SampleNature[2,i] states if sample whose number ID is equal to i, is an C-14 sample 1 or not `0.

PriorAge

numeric (with default): lower and upper bounds for age parameter of each sample (in ka). Note that, length(PriorAge) = 2*Nb_sample sand PriorAge[2i-1,2i] corresponds to the lower and upper bounds of sample whose number ID is equal to i.

SavePdf

logical (with default): if TRUE save graphs in pdf file named OutputFileName in folder OutputFilePath.

OutputFileName

character (with default): name of the pdf file that will be generated by the function if SavePdf=TRUE, length(OutputFileName)=3, see PLOT OUTPUT in Value section for more informations.

OutputFilePath

character (with default): path to the pdf file that will be generated by the function if SavePdf=TRUE. If it is not equal to "", it must be terminated by "/".

SaveEstimates

logical (with default): if TRUE save Bayes' estimates, credible interval at level 68% and 95% and the result of the Gelman en Rubin test of convergence, in a CSV-table named OutputFileName in folder OutputFilePath.

OutputTableName

character (with default): name of the table that will be generated by the function if SaveEstimates=TRUE.

OutputTablePath

character (with default): path to the table that will be generated by the function if SaveEstimates=TRUE. If it is not equal to "", it must be terminated by "/".

StratiConstraints

matrix or character (with default): input object for the stratigraphic relation between samples. If there is stratigraphic relation between samples see the details section for instructions regarding how to correctly fill StratiConstraints, the user can refer to a matrix or to a CSV-file character. Otherwise, default value is suitable.

sepSC

character (with default): if StratiConstraints is character, indicate column separator in StratiConstraints CSV-file.

BinPerSample

numeric (with default): vector with the number of BIN-files per OSL sample. The length of this vector is equal to the number of OSL samples. BinPerSample[i] corresponds to the number of BIN files for the sample whose number ID is equal to i. For more information to fill this vector, we refer to details in Generate_DataFile or in Generate_DataFile_MG.

THETA

numeric matrix or character (with default): input object for systematic and individual error for OSL samples. If systematic errors are considered, see the details section for instructions regarding how to correctly fill THETA; the user can refer to a matrix (numeric matrix) or to a csv file (character). Otherwise, default value is suitable, and only individual error is considered.

sepTHETA

character (with default): if THETA is character, indicate column separator in THETA CSV-file.

LIN_fit

logical (with default): if TRUE (default) allows a linear component, on top of the (default) saturating exponential curve, for the fitting of dose response curves, for OSL samples. See details for more informations on the proposed dose response curves.

Origin_fit

plogical (with default): if TRUE, forces the dose response curves to pass through the origin. See details for more informations on the proposed growth curves, for OSL samples.

distribution

character (with default): type of distribution that defines how individual equivalent dose values are distributed around the palaeodose, for OSL samples. Allowed inputs are "cauchy", "gaussian", "lognormal_A" and "lognormal_M", see details for more informations.

Model_C14

character (with default): if "full", error on estimate calibration curve is taken account, for C-14 samples. If "naive" this error is not taken account in the age estimate.

CalibrationCurve

character (with default): calibration curve chosen, for C-14 samples. Allowed inputs are

  • "Intcal13" or "Intcal13" for Northern Hemisphere atmospheric radiocarbon calibration curve,

  • "Marine13" or "Marine13" for Marine radiocarbon calibration curve,

  • "SHCal13" or "SHCal20" for Southern Hemisphere atmospheric radiocarbon calibration curve

  • a csv file, with tree columns, the first column is dedicated to "Cal.BP", the second to "XC-14.age", the third to "Error". The decimal of this file must be a dot, and the separator must be a comma.

Iter

integer (with default): the number of iterations to run and who will be used to assess convergence and ages (see runjags::run.jags]).

burnin

integer (with default): the number of iterations used to "home in" on the stationary posterior distribution. These are not used for assessing convergence (see runjags::run.jags]).

adapt

integer (with default): the number of iterations used in the adaptive phase of the simulation (see runjags::run.jags]).

t

numeric (with default): 1 every t iterations of the MCMC is considered for sampling the posterior distribution (for more information see rjags::jags.model.

n.chains

numeric (with default): number of independent chains for the model (for more information see [rjags::jags.model).

jags_method

character (with default): select which method to use in order to call JAGS, supported are "rjags" (the default), rjparallel, simple, interruptible, parallel, and snow (for more information about each of these possibilities, see runjags::run.jags])

autorun

logical (with default): choose to automate JAGS processing. JAGS model will be automatically extended until convergence is reached (for more information see runjags::autorun.jags).

quiet

logical (with default): enables/disables rjags::rjags-package messages

roundingOfValue

integer (with default): Integer indicating the number of decimal places to be used, default = 3.

...

further arguments that can be passed to control the Bayesian process, see details for supported arguments

Value

NUMERICAL OUTPUT

  1. A list containing the following objects:

    • Sampling: that corresponds to a sample of the posterior distributions of the age parameters (in ka for both C14 samples and OSL samples);

    • PriorAge: stating the priors used for the age parameter;

    • StratiConstraints: stating the stratigraphic relations between samples considered in the model;

    • Model_OSL_GrowthCurve: stating which dose response fitting option was chosen;

    • Model_OSL_Distribution: stating which distribution was chosen to model the dispersion of individual equivalent dose values around the palaeodose of the sample;

    • Model_C14: stating which model was chosen ("full" or "naive");

    • CalibrationCurve: stating which radiocarbon calibration curve was chosen;

    • Outlier: stating the names of samples that must be outliers.

  2. The Gelman and Rubin test of convergency: prints the result of the Gelman and Rubin test of convergence for the age estimate for each sample. A result close to one is expected.
    In addition, the user must visually assess the convergence of the trajectories by looking at the graph generated by the function (see PLOT OUTPUT for more informations).
    If both convergences (Gelman and Rubin test and plot checking) are satisfactory, the user can consider the estimates as valid. Otherwise, the user may try increasing the number of MCMC iterations (Iter) or be more precise on the PriorAge parameter to reach convergence.

  3. Credible intervals and Bayes estimates: prints the Bayes' estimates, the credible intervals at 95% and 68% for the age parameters for each sample.

  4. JAGS model output: returns the JAGS model with class "runjags".

PLOT OUTPUT

  1. MCMC trajectories: A graph with the MCMC trajectories and posterior distributions of the age parameter is displayed.
    On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.

  2. Age estimate and HPD at 95% of 14C samples on calibration curve: plot age estimate and HPD on calibration plot.

  3. Summary of sample age estimates: plot credible intervals and Bayes estimate of each sample age on a same graph.

Details

Note that there are three types of arguments in the previous list. There are arguments for information concerning only OSL samples: DATA, BinPerSample, THETA, sepTHETA, LIN_fit, Origin_fit, distribution.

There are arguments for information concerning only C14 samples: Data_C14Cal, Data_SigmaC14Cal, Model_C14, CalibrationCurve.

There are arguments for information concerning all the samples: Nb_sample, SampleNames, SampleNature, PriorAge, SavePdf, OutputFileName, OutputFilePath, SaveEstimates, OutputTableName, OutputTablePath, StratiConstraints, sepSC.

Supported ... arguments

ARGUMENTINPUTMETHODDEFAULTDESCRIPTION
max.timecharacterrjparallelInfmaximum allowed processing time, e.g., 10m for 10 minutes (cf. runjags::autorun.jags)
interactivelogicalrjparallelFALSEenable/disable interactive mode (cf. runjags::autorun.jags)
startburninintegerrjparallel4000number of burn-in iterations (cf. runjags::autorun.jags)
startsampleintegerrjparallel10000total number of samples to assess convergence (cf. runjags::autorun.jags)
initsnamed listrjparallelNAfine control over random seeds and random number generator runjags::autorun.jags

How to fill `StratiConstraints?

If there are stratigraphic relations between samples, 14C estimate age in Data_C14Cal must be ordered by order of increasing ages, as informations in DATA. Names in SampleNames must be ordered and corresponds to the order in Data_C14Cal and in DATA, also if it is needed to mix names of OSL samples and 14C samples.

The user can fill the StratiConstraints matrix as follow.

  1. Size of the matrix: row number of StratiConstraints matrix is equal to Nb_sample+1, and column number is equal to Nb_sample.

  2. First line of the matrix: for all i in {1,...,Nb_Sample}, StratiConstraints[1,i]=1 that means the lower bound of the sample age (given in PriorAge[2i-1]) for the sample whose number ID is equal to i, is taken into account.

  3. Sample relations: for all j in {2,...,Nb_Sample+1} and all i in {j,...,Nb_Sample}, StratiConstraints[j,i]=1 if sample age whose number ID is equal to j-1 is lower than sample age whose number ID is equal to i. Otherwise, StratiConstraints[j,i]=0.

Note that StratiConstraints_{2:Nb_sample+1,1:Nb_sample} is a upper triangular matrix.

The user can also use SCMatrix or SC_Ordered (if all samples are ordered) function to construct the StratiConstraints matrix.

The user can also refer to a csv file that contains the relation between samples as defined above. The user must be careful about which separator is used in the csv file using the argument sepSC.

How to fill THETA covariance matrix concerning common and individual error?

If systematic errors are considered, the user can fill the THETA matrix as follow.

  • row number of THETA is equal the column number, equal to Nb_sample.

  • For all i in {1,...,Nb_sample}, THETA[i,i] contains individual error plus systematic error of the sample whose number ID is equal to i.

  • For all i,j in {1,...,Nb_sample} and i different from j , THETA[i,j] contains common error between samples whose number ID are equal to i and j.

Note that THETA[i,j] is a symmetric matrix.

The user can also refer to a .csv file that contains the errors as defined above.

Option on growth curves

As for Age_Computation and Palaeodose_Computation, the user can choose from 4 dose response curves:

  • Saturating exponential plus linear growth (AgesMultiCS2_EXPLIN):

    for all x in IR+, f(x)=a(1-exp(-x/b))+cx+d; select

    • LIN_fit=TRUE

    • Origin_fit=FALSE

  • Saturating exponential growth (AgesMultiCS2_EXP):

    for all x in IR+, f(x)=a(1-exp(-x/b))+d; select

    • LIN_fit=FALSE

    • Origin_fit=FALSE

  • Saturating exponential plus linear growth and fitting through the origin (AgesMultiCS2_EXPLINZO):

    for all x in IR+, f(x)=a(1-exp(-x/b))+cx; select

    • LIN_fit=TRUE

    • Origin_fit=TRUE

  • Saturating exponential growth and fitting through the origin (AgesMultiCS2_EXPZO):

    for all x in IR+, f(x)=a(1-exp(-x/b)); select

    • LIN_fit=FALSE

    • Origin_fit=TRUE

Option on equivalent dose distribution around the palaeodose

The use can choose between :

  • cauchy: a Cauchy distribution with location parameter equal to the palaeodose of the sample

  • gaussian: a Gaussian distribution with mean equal to the palaeodose of the sample

  • lognormal_A: a log-normal distribution with mean or Average equal to the palaeodose of the sample

  • lognormal_M: a log-normal distribution with Median equal to the palaeodose of the sample

More precision on Model

We propose two models "full" or "naive". If Model='full' that means measurement error and error on calibration curve are taken account in the Bayesian model; if Model="naive" that means only error on measurement are taken account in the mode.

More precisely, the model considered here, as the one developed by Christen, JA (1994), assume multiplicative effect of errors to address the problem of outliers. In addition, to not penalise variables that are not outliers and damage theirs estimation, we introduce a structure of mixture, that means only variable that are considered as outlier have in addition a multiplicative error.

Note

Please note that the initial values for all MCMC are currently all the same for all chains since we rely on the automatic initial value generation of JAGS. This is not optimal and will be changed in future. However, it does not affect the quality of the age estimates if the chains have converged.

How to cite

Christophe, C., Philippe, A., Kreutzer, S., Baumgarten, F.H., 2024. Age_OSLC14(): Bayesian analysis for age estimation of OSL measurements and C-14 ages of various samples. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., Frerebeau, N., 2024. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.3.9000-13. https://CRAN.r-project.org/package=BayLum

References

Reimer PJ, Bard E, Bayliss A, Beck JW, Blackwell PC, Bronl Ramsey C, Buck CE, Cheng H, Edwards RL, Friedrich M, Grootes PM, Guilderson TP, Haflidason H, Hajdas I, Hatte C, Heaton TJ, Hoffmann DL, Hogg AG, Hughen KA, Kaiser KF, Kromer B, Manning SW, Niu M, Reimer RW, Richards DA, Scott EM, Southon JR, Staff RA, Turney CSM, van der Plicht J. 2013. IntCal13 ans Marine13 radiocarbon age calibration curves 0-50000 years cal BP. Radiocarbon 55(4)=1869-1887.

Hogg AG, Hua Q, Blackwell PG, Niu M, Buck CE, Guilderson TP, Heaton TJ, Palmer JG, Reimer PJ, Reimer RW, Turney CSM, Zimmerman SRH. 2013. SHCal13 Southern Hemisphere calibration, 0-50000 years cal BP. Radiocarbon 55(4):1889-1903

See also

runjags::run.jags, plot_MCMC, SCMatrix, plot_Ages

Author

Claire Christophe, Anne Philippe, Guillaume Guerin, Sebastian Kreutzer, Frederik Harly Baumgarten

Examples

## Load data
# OSL data
data(DATA1,envir = environment())
data(DATA2,envir = environment())
Data <- combine_DataFiles(DATA2,DATA1)

# 14C data
C14Cal <- DATA_C14$C14[1,1]
SigmaC14Cal <- DATA_C14$C14[1,2]
Names <- DATA_C14$Names[1]

# Prior Age
prior <- rep(c(1,60),3)
samplenature <- matrix(
 data = c(1,0,1,0,1,0),
 ncol = 3,
 nrow = 2,
 byrow = TRUE)

SC <- matrix(
 data = c(1,1,1,0,1,1,0,0,1,0,0,0),
 ncol = 3,
 nrow =4 ,
 byrow = TRUE)

## Age computation of samples
if (FALSE) { # \dontrun{
Age <- Age_OSLC14(
 DATA = Data,
 Data_C14Cal = C14Cal,
 Data_SigmaC14Cal = SigmaC14Cal,
 SampleNames = c("GDB5",Names,"GDB3"),
 Nb_sample = 3,
 SampleNature = samplenature,
 PriorAge = prior,
 StratiConstraints = SC,
 Iter = 20,
 burnin = 20,
 adapt = 20,
 n.chains = 2)
} # }