R/Age_OSLC14.R
Age_OSLC14.RdThis 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.
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,
...
)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.
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.
numeric: corresponding to the error of C-14 age estimates.
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).
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).
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.
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.
logical (with default): if TRUE save graphs in pdf file
named OutputFileName in folder OutputFilePath.
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.
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 "/".
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.
character (with default): name of the table that will
be generated by the function if SaveEstimates=TRUE.
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 "/".
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.
character (with default): if StratiConstraints is character,
indicate column separator in StratiConstraints CSV-file.
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.
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.
character (with default): if THETA is character,
indicate column separator in THETA CSV-file.
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.
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.
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.
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.
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.
integer (with default): the number of iterations to run and who will be used to assess convergence and ages (see runjags::run.jags]).
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]).
integer (with default): the number of iterations used in the adaptive phase of the simulation (see runjags::run.jags]).
numeric (with default): 1 every t iterations of the MCMC is
considered for sampling the posterior distribution (for more information see [rjags::jags.model.
numeric (with default): number of independent chains for the model (for more information see [rjags::jags.model).
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])
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).
logical (with default): enables/disables rjags messages
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
NUMERICAL OUTPUT
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.
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.
Credible intervals and Bayes estimates: prints the Bayes' estimates, the credible intervals at 95% and 68% for the age parameters for each sample.
JAGS model output: returns the JAGS model with class "runjags".
PLOT OUTPUT
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.
Age estimate and HPD at 95% of 14C samples on calibration curve: plot age estimate and HPD on calibration plot.
Summary of sample age estimates: plot credible intervals and Bayes estimate of each sample age on a same graph.
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
| ARGUMENT | INPUT | METHOD | DEFAULT | DESCRIPTION |
max.time | character | rjparallel | Inf | maximum allowed processing time, e.g.,
10m for 10 minutes (cf. runjags::autorun.jags) |
interactive | logical | rjparallel | FALSE | enable/disable interactive mode (cf. runjags::autorun.jags) |
startburnin | integer | rjparallel | 4000 | number of burn-in iterations (cf. runjags::autorun.jags) |
startsample | integer | rjparallel | 10000 | total number of samples to assess convergence (cf. runjags::autorun.jags) |
inits | named list | rjparallel | NA | fine 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.
Size of the matrix: row number of StratiConstraints matrix is equal to Nb_sample+1,
and column number is equal to Nb_sample.
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.
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.
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.
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.2. https://CRAN.r-project.org/package=BayLum
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
## 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) {
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)
}