A collection of various R functions for Bayesian analysis of luminescence data and C14 age estimates. This includes, amongst others, data import, export, application of age and palaeodose models.
This package is based on the functions: Generate_DataFile and Generate_DataFile_MG to import luminescence data. These functions create a list containing all informations to compute age of singlegrain OSL measurements for the first function and multigrain OSL measurements for the second.
The functions: Age_Computation and AgeS_Computation use Bayesian analysis for OSL age estimation for one or various samples according to difference models (e.g. different doseresponse curves and different equivalent dose distributions around the palaeodose).
It is possible to consider various BIN/BINXfiles per sample, to compute ages of samples in stratigraphic constraints and to integrate systematic errors.
It is possible to calibrate C14 age with the function AgeC14_Computation. We can also estimate chronology containing 14C age and OSL samples with the function Age_OSLC14.
Philippe, A., Guérin, G., Kreutzer, S., 2019. BayLum  An R package for Bayesian analysis of OSL ages: An introduction. Quaternary Geochronology 49, 16–24. doi:10.1016/j.quageo.2018.05.009
This function computes the age (in ka) of a sample according to the model developed in Combes and Philippe (2017),
based on an output of Generate_DataFile or Generate_DataFile_MG.
A sample, for which data is available in several BIN files, can be analysed.
Age_Computation(
DATA,
SampleName,
PriorAge = c(0.01, 100),
BinPerSample = c(1),
SavePdf = FALSE,
OutputFileName = c("MCMCplot"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(""),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
I = 1,
Iter = 50000,
t = 5,
n.chains = 3,
quiet = FALSE,
roundingOfValue = 3
)
DATA 
list of objects:

SampleName 
character: name of the sample. 
PriorAge 
numeric (with default): lower and upper bounds for the sample age parameter (in ka).
Note that, 
BinPerSample 
integer (with default): vector with the number of BIN files per sample.
If in 
SavePdf 
logical (with default): if TRUE save graph in pdf file named 
OutputFileName 
character (with default): name of the pdf file that will be generated by the function if 
OutputFilePath 
character (with default): path to the pdf file that will be generated by the function if 
SaveEstimates 
logical (with default): if TRUE save Bayes estimates and credible interval at level 68% and 95% and
the result of the gelman en Rubin test of convergency, in a csv table named 
OutputTableName 
character (with default): name of the table that will be generated by the function if 
OutputTablePath 
character (with default): path to the table that will be generated by the function if 
LIN_fit 
logical (with default): if 
Origin_fit 
logical (with default): if 
distribution 
character (with default): type of distribution that defines
how individual equivalent dose values are distributed around the palaeodose.
Allowed inputs are 
I 
integer (with default): if 
Iter 
integer (with default): number of iterations for the MCMC computation (for more information see jags.model). 
t 
integer (with default): 1 every 
n.chains 
integer (with default): number of independent chains for the model (for more information see jags.model). 
quiet 

roundingOfValue 
integer (with default): Integer indicating the number of decimal places to be used, default = 3. 
Option on growth curves
As for AgeS_Computation and Palaeodose_Computation, the user can choose from 4 dose response curves:
Saturating exponential plus linear growth (AgeMultiBF_EXPLIN
):
for all x
in IR+, $f(x)=a(1exp(x/b))+cx+d$
; select
LIN_fit=TRUE
Origin_fit=FALSE
Saturating exponential growth (AgeMultiBF_EXP
):
for all x
in IR+, $f(x)=a(1exp(x/b))+d$
; select
LIN_fit = FALSE
Origin_fit = FALSE
Saturating exponential plus linear growth and fitting through the origin (AgeMultiBF_EXPLINZO
):
for all x
in IR+, $f(x)=a(1exp(x/b))+cx$
; select
LIN_fit=TRUE
Origin_fit=TRUE
Saturating exponential growth and fitting through the origin (AgeMultiBF_EXPZO
):
for all x
in IR+, $f(x)=a(1exp(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 lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
: a lognormal distribution with Median equal to the palaeodose of the sample
NUMERICAL OUTPUT
A list containing the following objects:
Sampling that corresponds to a sample of the posterior distributions of the age (in ka), palaeodose (in Gy) and equivalent dose dispersion (in Gy) parameters.
Model_GrowthCurve, stating which dose response fitting option was chosen;
Distribution, stating which distribution was chosen to model the dispersion of individual equivalent dose values around the palaeodose of the sample;
PriorAge, stating the priors used for the age parameter (in ka).
The Gelman and Rubin test of convergency: prints the result of the Gelman and Rubin test of convergency for the age, palaeodose and equivalent dose dispersion parameters.
A result close to one is expected.
In addition, the user must visually assess the convergency of the trajectories by looking at the graph
generated by the function (see PLOT OUTPUT for more informations).
If both convergencies (Gelman and Rubin test and plot checking) are satisfactory,
the user can consider the printed estimates as valid. Otherwise, the user may try increasing the number of MCMC interations
(Iter
), or being more precise on the PriorAge
parameter (for example specify if it is a young sample c(0.01,10)
an old sample c(10,100)
),
or changing the parameter distribution
or the growth curve, to reach convergency.to reach convergency.
Credible intervals and Bayes estimates: prints the Bayes esitmates, the credible intervals at 95% and 68% for the age, palaeodose and equivalent dose dispersion parameters of the sample.
PLOT OUTPUT
A graph with the MCMC trajectories and posterior distributions of the age, palaeodose and equivalent dose dispersion parameters is displayed.
The first line of the figure correponds to the age parameter, the second to the palaeodose parameter and the third to the equivalent dose dispersion parameter.
On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
To give the results in a publication, we recommend to give the Bayes estimate of the parameter as well as the credible interval at 95% or 68%.
Christophe, C., Kreutzer, S., Philippe, A., Guérin, G., 2023. Age_Computation(): Bayesian analysis for the OSL age estimation of one sample. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
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.
Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guérin
Combes, Benoit and Philippe, Anne, 2017. Bayesian analysis of multiplicative Gaussian error for multiple ages estimation in optically stimulated luminescence dating. Quaternary Geochronology (39, 2434)
Combes, B., Philippe, A., Lanos, P., Mercier, N., Tribolo, C., Guerin, G., Guibert, P., Lahaye, C., 2015. A Bayesian central equivalent dose model for optically stimulated luminescence dating. Quaternary Geochronology 28, 6270. doi:10.1016/j.quageo.2015.04.001
Generate_DataFile, Generate_DataFile_MG, rjags, plot_MCMC, AgeS_Computation, Palaeodose_Computation
## load data file generated by the function Generate_DataFile
data(DATA1,envir = environment())
priorage < c(10,60) # GDB3 is an old sample
Age < Age_Computation(
DATA = DATA1,
SampleName = "GDB3",
PriorAge = priorage,
Iter = 100,
quiet = TRUE)
This function computes an age of OSL data consisting of at least two samples and calibrate
C14 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).
Singlegrain or Multigrain 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 C14 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,
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,
...
)
DATA 
Two types of inputs are supported.
(1): a list of objects: 
Data_C14Cal 
numeric vector: corresponding to C14 age estimate
(in years, conversion in ka is automatically done in the function).
If there is stratigraphic relations between samples, 
Data_SigmaC14Cal 
numeric: corresponding to the error of C14 age estimates. 
Nb_sample 
numeric: number of samples (OSL data and C14 age),
( 
SampleNames 
character: sample names for both OSL data and C14 data.
The length of this vector is equal to 
SampleNature 
matrix: states the nature of the sample.
Row number of 
PriorAge 
numeric (with default): lower and upper bounds for age parameter of each sample (in ka).
Note that, 
SavePdf 
logical (with default): if 
OutputFileName 
character (with default): name of the pdf file that
will be generated by the function if 
OutputFilePath 
character (with default): path to the pdf file that will
be generated by the function if 
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 CSVtable named 
OutputTableName 
character (with default): name of the table that will
be generated by the function if 
OutputTablePath 
character (with default): path to the table that will
be generated by the function if 
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 
sepSC 
character (with default): if 
BinPerSample 
numeric (with default): vector with the number of BINfiles
per OSL sample. The length of this vector is equal to the number of OSL samples.

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 
sepTHETA 
character (with default): if 
LIN_fit 
logical (with default): if 
Origin_fit 
plogical (with default): if 
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 
Model_C14 
character (with default): if 
CalibrationCurve 
character (with default): calibration curve chosen, for C14 samples. Allowed inputs are

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 
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 
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 

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 
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 burnin 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[2i1]
)
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 j1
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(1exp(x/b))+cx+d
; select
LIN_fit=TRUE
Origin_fit=FALSE
Saturating exponential growth (AgesMultiCS2_EXP
):
for all x
in IR+, f(x)=a(1exp(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(1exp(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(1exp(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 lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
: a lognormal 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.
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.
Christophe, C., Philippe, A., Kreutzer, S., Baumgarten, F.H., 2023. Age_OSLC14(): Bayesian analysis for age estimation of OSL measurements and C14 ages of various samples. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
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.
Claire Christophe, Anne Philippe, Guillaume Guerin, Sebastian Kreutzer, Frederik Harly Baumgarten
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
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, 050000 years cal BP. Radiocarbon 55(4):18891903
runjags, plot_MCMC, SCMatrix, plot_Ages
## 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
## Not run:
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)
## End(Not run)
This function calibrates the C14 age of samples to get an age (in ka). The user can choose one of the following radiocarbon calibration curve: Northern or Sourthen Hemisphere or marine atmospheric. It must be the same curve for all samples.
AgeC14_Computation(
Data_C14Cal,
Data_SigmaC14Cal,
SampleNames,
Nb_sample,
PriorAge = rep(c(10, 50), Nb_sample),
SavePdf = FALSE,
OutputFileName = c("MCMCplot", "HPD_CalC14Curve", "summary"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(""),
StratiConstraints = c(),
sepSC = c(","),
Model = c("full"),
CalibrationCurve = c("IntCal20"),
Iter = 50000,
t = 5,
n.chains = 3,
quiet = FALSE,
roundingOfValue = 3
)
Data_C14Cal 
numeric (required): corresponding to C14 age estimate. 
Data_SigmaC14Cal 
numeric (required): correponding to the error of C14 age estimates. 
SampleNames 
character (required): names of sample. The length of this vector is equal to 
Nb_sample 
integer: number of samples. 
PriorAge 
numeric (with default): lower and upper bounds for age parameter of each sample in years (not in ka).
Note that, 
SavePdf 
logical (with default): if TRUE save graphs in pdf file named 
OutputFileName 
character (with default): name of the pdf file that will be generated by the function if 
OutputFilePath 
character (with default): path to the pdf file that will be generated by the function if 
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 convergency, in a csv table named 
OutputTableName 
logical (with default): name of the table that will be generated by the function if 
OutputTablePath 
character (with default): path to the table that will be generated by the function if 
StratiConstraints 
numeric matrix or character(with default): input object for the statigraphic relation between samples.
If there is stratigraphic relation between samples see the details section for instructions regarding how to correctly fill 
sepSC 
character (with default): if 
Model 
character (with default): if "full", error on estimate calibration curve is taken account. If "naive" this error is not taken account in the age estimate. 
CalibrationCurve 
character (with default): calibration curve choosen. Allowed inputs are

Iter 
integer (with default): number of iterations for the MCMC computation (for more information see rjags::jags.model). 
t 
integer (with default): 1 every 
n.chains 
integer (with default): number of independent chains for the model (for more information see rjags::jags.model. 
quiet 

roundingOfValue 
integer (with default): Integer indicating the number of decimal places to be used, default = 3. 
** How to fill StratiConstraints? **
If there is stratigraphic relations between samples, C14 age in Data_C14Cal
must be ordered by order of increasing ages.
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[2i1]
)
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 j1
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) functions
to construct the StratiConstraints
matrix.
The user can also refer to a .csv file that containts the relation between samples as defined above.
The user must take care about the separator used in the csv file using the argument sepSC
.
** 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 developped by Christen, JA (1994), assume multiplicative effect of errors to address the problem of outliers. In addition, to not penalyse 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.
NUMERICAL OUTPUT
A list containing the following objects:
Sampling: that corresponds to a sample of the posterior distributions of the age parameters;
Outlier: stating the names of samples that are considered as outliers;
Model: stating which model was chosen ("full"
or "naive"
);
CalibrationCurve: stating which radiocarbon calibration curve was chosen;
PriorAge: stating the priors used for the age parameter;
StratiConstraints: stating the stratigraphic relations between samples considered in the model.
The Gelman and Rubin test of convergency: print the result of the Gelman and Rubin test of convergency for the age estimate for each sample.
A result close to one is expected.
In addition, the user must visually assess the convergency of the trajectories by looking at the graph
generated by the function (see PLOT OUTPUT for more informations).
If both convergencies (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 interations (Iter
)
or being more precise if it is possible on the PriorAge
parameter to reach convergency.
Credible intervals and Bayes estimates: prints the Bayes estimates, the credible intervals at 95% and 68% for the age parameters for each sample.
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.
Summary of sample age estimates: plot credible intervals and Bayes estimate of each sample age on one graph.
To give the results in a publication, we recommend to give the Bayes estimate of the parameters as well as the credible interval at 95% or 68%.
Christophe, C., Philippe, A., Guérin, G., Kreutzer, S., 2023. AgeC14_Computation(): Bayesian analysis for C14 age estimations of various samples. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
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.
Claire Christophe, Anne Philippe, Guillaume Guérin, Sebastian Kreutzer
Christen, JA (1994). Summarizing a set of radiocarbon determinations: a robust approach. Applied Statistics, 489503.
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
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, 050000 years cal BP. Radiocarbon 55(4):18891903
rjags, plot_MCMC, SCMatrix, plot_Ages
## Load data
data(DATA_C14,envir = environment())
C14Cal < DATA_C14$C14[,1]
SigmaC14Cal < DATA_C14$C14[,2]
Names < DATA_C14$Names
nb_sample < length(Names)
## Age computation of samples without stratigraphic relations
Age < AgeC14_Computation(
Data_C14Cal = C14Cal,
Data_SigmaC14Cal = SigmaC14Cal,
SampleNames = Names,
Nb_sample = nb_sample,
PriorAge = rep(c(20,60),nb_sample),
Iter = 500,
quiet = TRUE,
roundingOfValue=3)
AgeS_Computation
function for the samples: "GDB5" and "GDB3"
Output of AgeS_Computation
function for the samples: "GDB5" and "GDB3", there is no stratigraphic relation neither systematic errors.
data("AgeS")
A list containing
Sampling
MCMC.list that corresponds to a sample of the posterior distributions of the ages (in ka), palaeodoses (in Gy) and equivalent dose dispersions (in Gy) parameters of samples "GDB5" and "GDB3";
Model_GrowthCurve
stating which dose response fitting option was chosen to run the function
Distribution
stating which distribution was chosen to model the dispersion of individual equivalent dose values around the palaeodose of the sample;
PriorAge
stating the priors used for the age parameter (in ka);
StratiConstraints
stating the matrix of stratigraphic relations between samples considered in the model;
CovarianceMatrix
stating the covariance matrix of error used in the model, highlighting not common errors between samples in our cases (diagonal matrix).
Tribolo, C., Asrat, A., Bahain, J. J., Chapon, C., Douville, E., Fragnol, C., Hernandez, M., Hovers, E., Leplongeon, A., Martin, L, Pleurdeau, D, Pearson, O , Puaud, S, Assefa, Z. (2017). Across the Gap: Geochronological and Sedimentological Analyses from the Late PleistoceneHolocene Sequence of Goda Buticha, Southeastern Ethiopia. PloS one, 12(1), e0169418.
data(AgeS)
str(AgeS)
This function computes the age (in ka) of at least two samples
according to the model developed in Combès and Philippe (2017),
based on outputs 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.
Singlegrain or Multigrain OSL measurements can be analysed simultaneously.
AgeS_Computation(
DATA,
SampleNames,
Nb_sample,
PriorAge = rep(c(0.01, 100), Nb_sample),
BinPerSample = rep(1, Nb_sample),
SavePdf = FALSE,
OutputFileName = c("MCMCplot", "summary"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(""),
THETA = c(),
sepTHETA = c(","),
StratiConstraints = c(),
sepSC = c(","),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
model = NULL,
Iter = 10000,
burnin = 4000,
adapt = 1000,
t = 5,
n.chains = 3,
jags_method = "rjags",
autorun = FALSE,
quiet = FALSE,
roundingOfValue = 3,
...
)
DATA 
(required) Two inputs are possible:
(1): DATA list of objects: 
SampleNames 
character vector: names of samples. The length of this vector is equal to 
Nb_sample 
integer: number of samples, 
PriorAge 
numeric vector (with default): lower and upper bounds for age parameter of each sample (in ka).
Note that, 
BinPerSample 
integer vector (with default): vector with the number of BIN files per sample.
The length of this vector is equal to 
SavePdf 
logical (with default): if TRUE save graphs in pdf file named 
OutputFileName 
character (with default): name of the pdf file that will be generated by the function if 
OutputFilePath 
character (with default): path to the pdf file that will be generated by the function if 
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 
OutputTableName 
character (with default): name of the table that will be generated by the function if 
OutputTablePath 
character (with default): path to the table that will be generated by the function if 
THETA 
numeric matrix or character (with default): input object for systematic and individual error.
If systematic errors are considered, see the details section for instructions regarding how to correctly fill 
sepTHETA 
character (with default): if 
StratiConstraints 
numeric 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 
sepSC 
character (with default): if 
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. See details section for more informations on the proposed dose response curves. 
Origin_fit 
logical (with default): if TRUE, forces the dose response curves to pass through the origin. See details section for more informations on the proposed growth curves. 
distribution 
character (with default): type of distribution that defines
how individual equivalent dose values are distributed around the palaeodose.
Allowed inputs are 
model 
character (optional): allows to provide a custom model to the function as text string. Please note that if this option is chosen the parameter 
Iter 
integer (with default): the number of iterations to run which 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 
integer (with default): 1 every 
n.chains 
integer (with default): number of independent chains for the model (for more information see runjags::run.jags). 
jags_method 
character (with default): select which method to use in order to call JAGS. jags_methods 
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 
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. 1) When creating a new JAGS model, see details for supported arguments. 2) If extending a JAGS model see arguments of runjags::extend.JAGS. 
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 burnin 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 is stratigraphic relations between samples,
informations in DATA must be ordered by order of increasing ages.
To do this the user can either fill right Names
in Generate_DataFile
or in Generate_DataFile_MG (as it is indicated in Details section of these function),
or ordered by order of increasing ages outputs of Generate_DataFile
or Generate_DataFile_MG in combine_DataFiles
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[2i1]
)
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 j1
is lower than sample age whose number ID is equal to i
.
Otherwise, StratiConstraints[j,i]=0
.
Note that StratiConstraints_{2:Nb_sample+A,1:Nb_sample}
is a upper triangular matrix.
The user can also use SCMatrix
or SC_Ordered
(if all samples are ordered) functions
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 take care about the separator 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 follows.
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 symetric matrix.
The user can also refer to a .csv file that contains the errors as defined above.
Alternatively you can use the function create_ThetaMatrix.
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(1exp(x/b))+cx+d
; select
LIN_fit=TRUE
Origin_fit=FALSE
Saturating exponential growth (AgesMultiCS2_EXP
):
for all x
in IR+, f(x)=a(1exp(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(1exp(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(1exp(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 lognormal distribution with mean or Average equal to the palaeodose of the sample:
lognormal_M
: a lognormal distribution with Median equal to the palaeodose of the sample.
NUMERICAL OUTPUT
A list of type BayLum.list
containing the following objects:
Sampling: that corresponds to a sample of the posterior distributions of the age (in ka), palaeodose (in Gy) and equivalent dose dispersion (in Gy) parameters for each sample;
Model_GrowthCurve: stating which dose response fitting option was chosen;
Distribution: stating which distribution was chosen to model the dispersion of individual equivalent dose values around the palaeodose of the sample;
PriorAge: stating the priors used for the age parameter (in ka);
StratiConstraints: stating the stratigraphic relations between samples considered in the model;
CovarianceMatrix: stating the covariance matrix of error used in the model, highlighting common errors between samples or not.
model: returns the model that was used for the Bayesian modelling as a character
JAGS model output: returns the JAGS model with class "runjags".
The Gelman and Rubin test of convergency: prints the result of the Gelman and Rubin test of convergence for
the age, palaeodose and equivalent dose dispersion parameters 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 being more precise on the PriorAge
parameter (for example specify if it is a young sample c(0.01,10)
an old sample c(10,100)
),
or changing the parameter distribution
or the growth curve, to reach convergence.
Credible intervals and Bayes estimates: prints the Bayes estimates, the credible intervals at 95% and 68% for the age, palaeodose and equivalent dose dispersion parameters for each sample.
PLOT OUTPUT
MCMC trajectories: A graph with the MCMC trajectories and posterior distributions of the age, palaeodose and equivalent dose dispersion parameters
is displayed, there is one page per sample.
The first line of the figure corresponds to the age parameter, the second to the palaeodose parameter
and the third to the equivalent dose dispersion parameter.
On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
Summary of sample age estimates: plot credible intervals and Bayes estimate of each sample age on a same graph.
To give the results in a publication, we recommend to give the Bayes' estimate of the parameters as well as the credible interval at 95% or 68%.
Christophe, C., Philippe, A., Guérin, G., Kreutzer, S., 2023. AgeS_Computation(): Bayesian analysis for OSL age estimation of various samples. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
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.
Claire Christophe, Anne Philippe, Guillaume Guérin, Sebastian Kreutzer
Combes, Benoit and Philippe, Anne, 2017. Bayesian analysis of multiplicative Gaussian error for multiple ages estimation in optically stimulated luminescence dating. Quaternary Geochronology (39, 2434)
Combes, B., Philippe, A., Lanos, P., Mercier, N., Tribolo, C., Guerin, G., Guibert, P., Lahaye, C., 2015. A Bayesian central equivalent dose model for optically stimulated luminescence dating. Quaternary Geochronology 28, 6270. doi:10.1016/j.quageo.2015.04.001
Generate_DataFile, Generate_DataFile_MG, rjags, plot_MCMC, SCMatrix, Age_Computation, Palaeodose_Computation, plot_Ages, create_ThetaMatrix, runjags
## Age computation of samples GDB5 and GDB3,
## load data
data(DATA2) # GD85
data(DATA1) # GD83
## produce DataFile
Data < combine_DataFiles(DATA2, DATA1)
## without common error, assuming GDB5 age younger than GDB3 age
Nb_sample < 2
SC < matrix(
data = c(1,1,0,1,0,0),
ncol = 2,
nrow = (Nb_sample+1),
byrow = TRUE)
## Not run:
## run standard
Age < AgeS_Computation(
DATA = Data,
Nb_sample = Nb_sample,
SampleNames = c("GDB5","GDB3"),
PriorAge = rep(c(1,100), 2),
StratiConstraints = SC,
Iter = 100,
quiet = FALSE,
jags_method = "rjags"
)
## extend model
Age_extended < AgeS_Computation(
DATA = Age,
Nb_sample = Nb_sample,
SampleNames = c("GDB5","GDB3"),
PriorAge = rep(c(1,100), 2),
StratiConstraints = SC,
adapt = 0,
burnin = 500,
Iter = 1000,
quiet = FALSE,
jags_method = "rjags"
)
## autorun mode
Age < AgeS_Computation(
DATA = Data,
Nb_sample = Nb_sample,
SampleNames = c("GDB5","GDB3"),
PriorAge = rep(c(1,100), 2),
StratiConstraints = SC,
Iter = 10000,
quiet = FALSE,
jags_method = "rjags",
autorun = TRUE)
## parallel mode
Age < AgeS_Computation(
DATA = Data,
Nb_sample = Nb_sample,
SampleNames = c("GDB5","GDB3"),
PriorAge = rep(c(1,100), 2),
StratiConstraints = SC,
Iter = 10000,
quiet = FALSE,
jags_method = "rjparallel")
## End(Not run)
Combine objects generated by Generate_DataFile and Generate_DataFile_MG
combine_DataFiles(...)
Concat_DataFile(...)
... 
list objects generated by Generate_DataFile or Generate_DataFile_MG 
The function allows to combine data already generated by Generate_DataFile or
Generate_DataFile_MG. The number of input objects is not limited and the function
works similar to the standard base R function c()
, but preserves the particular structure of the
objects imported and generated by 'BayLum'. The elements are combined by list
element names.
Combining such data is rather useful in two scenarios:
The data have been already imported and treated and then stored in RDatafiles. Using the function
combine_DataFiles()
will significantly speed up the processing time,
simultaneous analysis of single and multigrain OSL measurements.
A nested list combining the input objects.
0.1.1
Kreutzer, S., Christophe, C., 2023. combine_DataFiles(): Combine objects. Function version 0.1.1. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, IRAMATCRP2A, UMR 5060, CNRS  Université Bordeaux Montaigne (France), adapting the idea from the function 'Concat_DataFile()' by Claire Christophe.
Generate_DataFile, Generate_DataFile_MG
# load data files
data(DATA1,envir = environment())
data(DATA2,envir = environment())
#combine objects
DATA3 < combine_DataFiles(DATA1, DATA2)
str(DATA3)
Create file and folder structure templates on the user hard drive
as expected by Generate_DataFile and Generate_DataFile_MG. Files and data in the
folders must then be overwritten manually with user data. The function intends to
minimise the errors going along with the creation of these folder structures.
The function uses the example data of BayLum
to create the templates.
create_FolderTemplates(
path,
mode = "SG",
n_folders = 1,
names = paste("Sample_", 1:n_folders),
verbose = TRUE
)
path 
character (required): path to the folder where the templates should be created 
mode 
character (with default): depending on the dataset you can
create templates or single grain ( 
n_folders 
numeric (with default): number of template folders to be created 
names 
character (optional): allows give own names to the subfolders. 
verbose 
logical (with default): enables/disables verbose mode 
If the templates were created successfully on the hard drive, the function returns nothing.
0.1.0
Kreutzer, S., 2023. create_FolderTemplates(): Create Folder Templates. Function version 0.1.0. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, Geography & Earth Sciences, Aberystwyth University (United Kingdom)
Generate_DataFile, Generate_DataFile_MG
create_FolderTemplates(tempdir())
Create the $\Theta$
matrix with the shared uncertainties
that can used as input in, e.g., AgeS_Computation and Age_OSLC14 which is used for the
covariance matrix $\Sigma$
(Combès & Philippe, 2017)
create_ThetaMatrix(input, output_file = NULL, sigma_s, ...)
input 
character or data.frame (optional): input data frame or file connection
to import a CSVfile with the needed information. If nothing is provided the function returns
an input template. The argument 
output_file 
character (optional): file path for the output CSVfile, the field separator
is hard set to 
sigma_s 
numeric (required): named character with values for systematic uncertainties. Those values
are labspecific. Can be set to 
... 
further arguments that can be passed to utils::read.table (for the CSVfile import) 
The function intends to ease the creation of the $Theta$
matrix, which cannot be
created straight forward, e.g., base R functions such as stats::cov
.
The relationship between the covariance matrix $Sigma$
and $Theta$
is given with
$\Sigma_ij = A_i * A_j * \Theta_ij$
For details see Combès & Philippe, 2017 and Guérin et al. (2021).
Input modes
The function supports two different operation modes:
input
is left empty: the function returns a data.frame template that can be used as input (the option output_file
works as well)
input
is fed with a data.frame or a character
(file path), the $\Theta$
matrix is returned
Input format
The function expects either a CSVfile or a data.frame as input. To create template you can
run the function leaving the argument input
empty (see example). Please note the format
of the input table (data.frame) needs to kept as specified in the template.
The following table lists the meaning of the columns:
COLUMN  DESCRIPTION  UNIT 
SAMPLE_ID 
sample name   
DR_BETA_K 
standard error betadose rate K  Gy/ka 
DR_BETA_U 
standard error betadose rate U  Gy/ka 
DR_BETA_Th 
standard error betadose rate Th  Gy/ka 
DR_GAMMA_K 
standard error gammadose rate K  Gy/ka 
DR_GAMMA_U 
standard error gammadose rate U  Gy/ka 
DR_GAMMA_Th 
standard error gammadose rate Th  Gy/ka 
DR_GAMMA_TOTAL 
standard error total gammadose rate  Gy/ka 
DR_TOTAL 
total dose rate  Gy/ka 
DR_TOTAL_X 
standard error total dose rate  Gy/ka 
Note: All columns can be set to 0 or NA
but no column must be left empty! If a value > 0 is provided
for DR_GAMMA_TOTAL
this value is taken and values in, e.g., DR_GAMMA_K
are discarded (set to 0)!
Systematic uncertainties
The following table provides information on the named argument
that can be provided via the argument sigma_s
. Missing values are not allowed, all
values must be set.
ARGUMENT  DESCRIPTION  UNIT 
s_betaK 
relative uncertainty K concentration   
s_betaU 
relative uncertainty U concentration   
s_betaTh 
relative uncertainty Th concentration   
s_gammaK 
relative uncertainty K concentration   
s_gammaU 
relative uncertainty U concentration   
s_gammaTh 
relative uncertainty Th concentration   
s_gammaDR 
relative uncertainty gammadose rate   
s_CAL 
relative uncertainty betasource calibration   
s_intDR 
absolute uncertainty internal dose rate  Gy/ka 
A symmetric $Theta$
matrix or if input
is missing, a data.frame with an input
template
0.1.0
Kreutzer, S., Guérin, G., 2023. create_ThetaMatrix(): Create Theta Matrix. Function version 0.1.0. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, IRAMATCRP2A, UMR 5060, CNRSUniversité Bordeaux Montaigne (France), based on an 'MS Excel' sheet by Guillaume Guérin, IRAMATCRP2A, UMR 5060, CNRSUniversité Bordeaux Montaigne (France)
Combès, B., Philippe, A., 2017. Bayesian analysis of individual and systematic multiplicative errors for estimating ages with stratigraphic constraints in optically stimulated luminescence dating. Quaternary Geochronology 39, 24–34. doi:10.1016/j.quageo.2017.02.003
Guérin, G., Lahaye, C., Heydari, M., Autzen, M., Buylaert, J.P., Guibert, P., Jain, M., Kreutzer, S., Lebrun, B., Murray, A.S., Thomsen, K.J., Urbanova, P., Philippe, A., 2021. Towards an improvement of optically stimulated luminescence (OSL) age uncertainties: modelling OSL ages with systematic errors, stratigraphic constraints and radiocarbon ages using the R package BayLum. Geochronology 3, 229—245. doi:10.5194/gchron32292021
AgeS_Computation, Age_OSLC14, utils::read.table, utils::write.table
##(1) return template data.frame (no file output)
create_ThetaMatrix()
## Not run:
##(2) return template as data.frame + file
file_path < tempfile(fileext = ".csv")
create_ThetaMatrix(output_file = file_path )
##NOT RUNNING EXAMPLE for sigma_s
calc_ThetaMatrix(...,
sigma_s = c(
s_betaK = 0.010,
s_betaU = 0.007,
s_betaTh = 0.006,
s_gammaK = 0.010,
s_gammaU = 0.007,
s_gammaTh = 0.006,
s_gammaDR = 0.05,
s_CAL = 0.020,
s_intDR = 0.030))
## End(Not run)
C14 cal age estiamtes and theirs error of samples SEVA26510, SEVA26506, SEVA26507, SEVA26508.
data("DATA_C14")
A list containing:
Names:
character vector of the sample names;
C14:
numeric matrix, in the first column the 14C Cal age of the samples, and in the second column theirs errors.
For more informations on this sample we refer to the following publication:
Guerin, G., Frouin, M., Talamo, S., Aldeias, V., Bruxelles, L., Chiotti, L., Goldberg, P., Hublin, J.J., Jain, M., Lahaye, C., Madelaine, S., Maureille, B., McPherron, S., Mercier, N., Murray, A., Sandgathe, D., Steele, T., Thomsen, K., Turq, A. (2015). A multimethod luminescence dating of the Palaeolithic sequence of La Ferrassie based on new excavations adjacent to the La Ferrassie 1 and 2 skeletons. Journal of Archaeological Science, 58, 147166.
data(DATA_C14)
(DATA_C14)
list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J,K,Nb_measurement obtained using Generate_DataFile
function with singlegrain OSL measurementsl of the sample GDB3.
data("DATA1")
A list containing:
LT:
(one list per sample): each list contains all L/T values for the corresponding sample;
sLT:
(one list per sample): each list contains all uncertainties on L/T values for the corresponding sample;
ITimes:
(one list per sample): each list contains irradiation time values for the corresponding sample;
dLab=
a matrix containing in line i
, the laboratory dose rate and its variance for sample i
;
ddot_env:
a matrix containing in line i
, the environmental dose rate and its variance (excluding the common error terms) for sample i
;
regDose:
(one list per sample): each list contains all regenerated doses;
J:
a vector giving, for each BIN file, the number of aliquots selected for the analysis;
K:
a vector giving, for each BIN file, the number of regenerative doses in the SAR protocol;
Nb_measurement:
a vector giving, for each BIN file, the number of measurements;
For more informations on this sample we refer to the following publication:
Tribolo, C., Asrat, A., Bahain, J. J., Chapon, C., Douville, E., Fragnol, C., Hernandez, M., Hovers, E., Leplongeon, A., Martin, L, Pleurdeau, D, Pearson, O , Puaud, S, Assefa, Z. (2017). Across the Gap: Geochronological and Sedimentological Analyses from the Late PleistoceneHolocene Sequence of Goda Buticha, Southeastern Ethiopia. PloS one, 12(1), e0169418.
data(DATA1)
str(DATA1)
list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J,K,Nb_measurement obtained using Generate_DataFile
function with singlegrain OSL measurementsl of the sample GDB5.
data("DATA2")
A data frame containing:
LT:
(one list per sample): each list contains all L/T values for the corresponding sample;
sLT:
(one list per sample): each list contains all uncertainties on L/T values for the corresponding sample;
ITimes:
(one list per sample): each list contains irradiation time values for the corresponding sample;
dLab:
a matrix containing in line i
, the laboratory dose rate and its variance for sample i
;
ddot_env:
a matrix containing in line i
, the environmental dose rate and its variance (excluding the common error terms) for sample i
;
regDose:
(one list per sample): each list contains all regenerated doses;
J:
a vector giving, for each BIN file, the number of aliquots selected for the analysis;
K:
a vector giving, for each BIN file, the number of regenerative doses in the SAR protocol;
Nb_measurement:
, a vector giving, for each BIN file, the number of measurements;
For more informations on this sample we refer to the following publication:
Tribolo, C., Asrat, A., Bahain, J. J., Chapon, C., Douville, E., Fragnol, C., Hernandez, M., Hovers, E., Leplongeon, A., Martin, L, Pleurdeau, D, Pearson, O , Puaud, S, Assefa, Z. (2017). Across the Gap: Geochronological and Sedimentological Analyses from the Late PleistoceneHolocene Sequence of Goda Buticha, Southeastern Ethiopia. PloS one, 12(1), e0169418.
data(DATA2)
str(DATA2)
list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J,K,Nb_measurement obtained using Generate_DataFile
function with multigrain OSL measurementsl of the sample FER1.
data("DATA3")
A list containing:
LT:
(one list per sample): each list contains all L/T values for the corresponding sample;
sLT:
(one list per sample): each list contains all uncertainties on L/T values for the corresponding sample;
ITimes:
(one list per sample): each list contains irradiation time values for the corresponding sample;
dLab=
a matrix containing in line i
, the laboratory dose rate and its variance for sample i
;
ddot_env:
a matrix containing in line i
, the environmental dose rate and its variance (excluding the common error terms) for sample i
;
regDose:
(one list per sample): each list contains all regenerated doses;
J:
a vector giving, for each BIN file, the number of aliquots selected for the analysis;
K:
a vector giving, for each BIN file, the number of regenerative doses in the SAR protocol;
Nb_measurement:
a vector giving, for each BIN file, the number of measurements;
For more informations on this sample we refer to the following publication:
Guerin, G., Frouin, M., Talamo, S., Aldeias, V., Bruxelles, L., Chiotti, L., Goldberg, P., Hublin, J.J., Jain, M., Lahaye, C., Madelaine, S., Maureille, B., McPherron, S., Mercier, N., Murray, A., Sandgathe, D., Steele, T., Thomsen, K., Turq, A. (2015). A multimethod luminescence dating of the Palaeolithic sequence of La Ferrassie based on new excavations adjacent to the La Ferrassie 1 and 2 skeletons. Journal of Archaeological Science, 58, 147166.
data(DATA3)
str(DATA3)
This function is used to generate, from the BIN file(s), a list of values of:
Singlegrain OSL intensities and associated uncertainties, regenerative doses, etc., which will be the input of the Bayesian models.
To be easytouse, this function requires a rigorous organisation  all needed files should be arranged in one folder 
of informations concerning each BIN file.
It is possible to process data for various samples simultaneously and to consider more than one BIN file per sample.
Generate_DataFile(
Path,
FolderNames,
Nb_sample,
Nb_binfile = length(FolderNames),
BinPerSample = rep(1, Nb_sample),
sepDP = c(","),
sepDE = c(","),
sepDS = c(","),
sepR = c("="),
verbose = TRUE,
...
)
Path 
character (required): the path to the project folder, containing one or more sub folders in which the BIN files are located. If it is not equal to "", it must be terminated by "/". 
FolderNames 
character (required): list of names of the subfolders containing the BIN files

Nb_sample 
integer (required): number of samples. 
Nb_binfile 
integer (with default): number of BIN files. It must be equal to, or greater than 
BinPerSample 
integer vector (with default): vector with the number of BIN files per sample.
The length of this vector must be equal to 
sepDP 
character (with default): column separator in the 
sepDE 
character (with default): column separator in the 
sepDS 
character (with default): column separator in the 'DoseLab.csv“ files. 
sepR 
character (with default): column separator in the Rule.csv files. 
verbose 
logical (with default): enable/disable verbose mode 
... 
further arguments that can be passed to Luminescence::read_BIN2R. 
With Path
and FolderNames
, this function goes to the sub folders containing the BIN files and associated information to compute
the luminescence data.
** What are the required files in each subfolder? **
Each sub folder can be named, for example, as the sample name followed by a number; it must contain:
bin.BIN: the bin file renamed as bin.BIN (note: the name of all files matters);
DiscPos.csv: a two columns csv file containing the list of disc and grain position number of the previously selected grains (typically this list will include the position of grains based on their sensitivity, recycling or other properties);
DoseEnv.csv: a two columns file containing the observation of the natural (or environmental), dose rate, and its nonshared variance (i.e. after removing all shared errors), both in Gy. Note: the user shall provide the squared value of the error associated with the dose rate experienced by the sample grains in nature;
DoseSourve.csv: a two columns file containing the observation of the laboratory dose rate, and its variance (squared error) both in Gy;
rule.csv: a csv file containing information on
beginSignal=
the first channel for summing the natural or regenerative OSL signal (typically 1 or 6);
endSignal=
the last channel for summing the natural or regenerative OSL signal (typically 5 or 10);
beginBackground=
the first channel for background estimation of the natural or regenerative OSL signal (typically 76 or 81);
endBackground=
the last channel for background estimation of the natural or regenerative OSL signal (typically 95 or 100);
beginTest=
,
endTest=
,
beginTestBackground=
,
endTestBackground=
same values as above, for the test dose response (typically the same values should be used);
inflatePercent=
uncertainty arising from the instrument reproducibility (typically 0.02, i.e. 2\
nbOfLastCycleToRemove=
number of cycles at the end of the SAR protocol which should not be included in the dose response curve fitting
(typically 1 if only a recycling test is performed, or 2 if both recycling and IR depletion are tested).
** How to fill the FolderNames
vector? **
FolderNames
is a vector of length Nb_binfile
. FolderNames[i]
is the name (e.g., Sample1File1, or successive names separated by "/" signs,
if BIN files are in subfolders, e.g. Sample1/File1) of the subfolder containing all informations on the BIN file of ID number i
.
The names in FolderNames
are ordered following two rules:
The names in the FolderNames
vector must be ordered following the sample order
(the names of subfolders containing BIN files for the same sample should follow each other in the FolderNames
vector, e.g.
Sample1, Sample2File1, Sample2File2, etc.).
If stratigraphic constraints apply to samples, and so a Bayesian model with stratigraphic constraints is implemented,
then the names in the FolderNames
vector must be ordered by order of increasing ages.
For example, FolderNames=c(noun1,noun2)
, in which case noun1
(respectively, noun2
) corresponds to the subfolder
name containing the BIN file of sample 1 (respectively of sample 2).
In addition, if we know that sample 1 is younger than sample 2, then FolderNames
vector is correctly filled.
If conversely, FolderNames=c(noun2,noun1)
, the analysis performed by AgeS_Computation
would not be consistent.
** How to fill the BinPerSample
vector? **
BinPerSample[i]
correponds to the number of BIN files for the sample whose number ID is equal to i
.
For example, let us consider a case with two samples (Sample1 and Sample2), with 2 BIN files for Sample1 and 1 for Sample2.
In this case, Nb_binfile
=3 and Nb_sample
=2.
The user may then set FolderNames=c("Sample1File1", "Sample1File2", "Sample2File1")
, in which case "Sample1File1"
is the name of the subfolder
containing the first BIN file for Sample1, "Sample1File2"
the name of the subfolder for the second BIN file of Sample1; eventually,
"Sample2File1"
is the name of the subfolder containing the BIN file for the second sample. In this case, BinPerSample=c(2,1)
.
For the general BINfile structure, the reader is referred to the following website:
http://www.nutech.dtu.dk/
The function read_BIN2R
developped in Luminescence
package is used to read the BIN files.
A list containing the following objects:
LT (one list per sample); each list contains all L/T values for the corresponding sample;
sLT (one list per sample); each list contains all uncertainties on L/T values for the corresponding sample;
ITimes (one list per sample); each list contains irradiation time values for the corresponding sample;
dLab, a matrix containing in line i
, the laboratory dose rate and its variance for sample i
;
ddot_env, a matrix containing in line i
, the environmental dose rate and its variance (excluding the common error terms) for sample i
;
regDose (one list per sample); each list contains all regenerated doses;
J, a vector giving, for each BIN file, the number of aliquots selected for the analysis;
K, a vector giving, for each BIN file, the number of regenerative doses in the SAR protocol;
Nb_measurement, a vector giving, for each BIN file, the number of measurements.
** How to save this list **
You can save this list in a .RData object. To do this, you can use the fonction save
.
Then, to load this list you can use the function load
(see example section fore more details).
Christophe, C., Kreutzer, S., Philippe, A., 2023. Generate_DataFile(): Generates, from one (or several) BINfile(s) of Singlegrain OSL measurements, a list of luminescence data and information before statistical analysis. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guerin
read_BIN2R
, combine_DataFiles
, Generate_DataFile_MG
, LT_RegenDose
Age_Computation
, AgeS_Computation
, Palaeodose_Computation
## Example for one sample with one Bin File
path< system.file("extdata/samp1", "", package="BayLum")
folder=""
nbsample=1 # give the number of sample
Data < Generate_DataFile(
Path = path,
FolderNames = folder,
Nb_sample = nbsample,
verbose = FALSE)
str(Data)
## to save information in RData object in folder containing bin file
# save(Data,file=c(paste(path,folder,'Data.RData',sep="")))
## to load information containing Data.RData object
# load(file=c(paste(path,folder,"Data.RData",sep="")))
This function is used to generate, from the BIN file(s), a list of values of:
Multigrain OSL intensities and associated uncertainties, regenerative doses, etc., which will be the input of the Bayesian models.
To be easytouse, this function requires a rigorous organisation  all needed files should be arranged in one folder 
of informations concerning each BIN file.
It is possible to process data for various samples simultaneously and to consider more
than one BINfile per sample.
Generate_DataFile_MG(
Path,
FolderNames,
Nb_sample,
Nb_binfile = length(FolderNames),
BinPerSample = rep(1, Nb_sample),
sepD = c(","),
sepDE = c(","),
sepDS = c(","),
sepR = c("="),
verbose = TRUE,
force_run1_at_a_time = FALSE,
...
)
Path 
character (required): the path to the project folder, containing one or more sub folders in which the BIN files
are located. If it is not equal to 
FolderNames 
character (required) vector: list of names of the subfolders containing the BIN files

Nb_sample 
integer (required): number of samples 
Nb_binfile 
integer (with default): number of BIN files. It must be equal to, or greater than 
BinPerSample 
integer vector (with default): vector with the number of BIN files per sample.
The length of this vector must be equal to 
sepD 
character (with default): column separator in the 
sepDE 
character (with default): column separator in the 
sepDS 
character (with default): column separator in the 
sepR 
character (with default): column separator in the 
verbose 
logical (with default): enable/disable verbose mode 
force_run1_at_a_time 
logical (with default): if set to 
... 
further arguments that can be passed to Luminescence::read_BIN2R. 
With Path
and FolderNames
, this function goes to the sub folders containing the BIN files and associated information to compute
the luminescence data.
** What are the required files in each subfolder? **
Each subfolder can be named, for example, as the sample name followed by a number; it must contain:
bin.BIN, the bin file renamed as bin.BIN (note: the name of all files matters);
Disc.csv, a one columns csv file containing the list of disc number of the previously selected grains (typically this list will include the position of grains based on their sensitivity, recycling or other properties);
DoseEnv.csv, a two columns file containing the observation of the natural (or environmental), dose rate, and its nonshared variance (i.e. after removing all shared errors), both in Gy. Note: the user shall provide the squared value of the error associated with the dose rate experienced by the sample grains in nature;
DoseSourve.csv, a two columns file containing the observation of the laboratory dose rate, and its variance (squared error), both in Gy;
rule.csv, a csv file containing information on
beginSignal=
the first channel for summing the natural or regenerative OSL signal (typically 1 or 6);
endSignal=
the last channel for summing the natural or regenerative OSL signal (typically 5 or 10);
beginBackground=
the first channel for background estimation of the natural or regenerative OSL signal (typically 76 or 81);
endBackground=
the last channel for background estimation of the natural or regenerative OSL signal (typically 95 or 100);
beginTest
,
endTest
,
beginTestBackground
,
endTestBackground=
same values as above, for the test dose response (typically the same values should be used);
inflatePercent=
uncertainty arising from the instrument reproducibility (typically 0.02, i.e. 2\
nbOfLastCycleToRemove=
number of cycles at the end of the SAR protocol which should not be included in the dose response curve fitting
(typically 1 if only a recycling test is performed, or 2 if both recycling and IR depletion are tested).
** How to fill the FolderNames
vector? **
FolderNames
is a vector of length Nb_binfile
. FolderNames[i]
is the name (e.g., Sample1File1, or successive names separated by "/" signs,
if BIN files are in subfolders, e.g. Sample1/File1) of the subfolder containing all informations on the BIN file of ID number i
.
The names in FolderNames
are ordered following two rules:
The names in the FolderNames
vector must be ordered following the sample order
(the names of subfolders containing BIN files for the same sample should follow each other in the FolderNames
vector, e.g.
Sample1, Sample2File1, Sample2File2, etc.).
If stratigraphic constraints apply to samples, and so a Bayesian model with stratigraphic constraints is implemented,
then the names in the FolderNames
vector must be ordered by order of increasing ages.
For example, FolderNames=c(noun1,noun2)
, in which case noun1
(respectively, noun2
) corresponds to the subfolder
name containing the BIN file of sample 1 (respectively of sample 2).
In addition, if we know that sample 1 is younger than sample 2, then FolderNames
vector is correctly filled.
If conversely, FolderNames=c(noun2,noun1)
, the analysis performed by AgeS_Computation
would not be consistent.
** How to fill the BinPerSample
vector? **
BinPerSample[i]
corresponds to the number of BIN files for the sample whose number ID is equal to i
.
For example, let us consider a case with two samples (Sample1 and Sample2), with 2 BIN files for Sample1 and 1 for Sample2.
In this case, Nb_binfile
=3 and Nb_sample
=2.
The user may then set FolderNames=c("Sample1File1", "Sample1File2", "Sample2File1")
, in which case "Sample11"
is the name of the subfolder containing the first BIN file
for Sample1, "Sample1File2"
the name of the subfolder for the second BIN file of Sample1; eventually, "Sample21"
is the name of the subfolder containing the BIN file
for the second sample. In this case, BinPerSample=c(2,1)
.
For the general BINfile structure, the reader is referred to the following website: http://www.nutech.dtu.dk/
The function Luminescence::read_BIN2R is used to read the BIN files.
A list containing the following objects:
LT (one list per sample); each list contains all L/T values for the corresponding sample;
sLT (one list per sample); each list contains all uncertainties on L/T values for the corresponding sample;
ITimes (one list per sample); each list contains irradiation time values for the corresponding sample;
dLab, a matrix containing in line i
, the laboratory dose rate and its variance for sample i
;
ddot_env, a matrix containing in line i
, the environmental dose rate and its variance (excluding the common error terms) for sample i
;
regDose (one list per sample); each list contains all regenerated doses;
J, a vector giving, for each BIN file, the number of aliquots selected for the analysis;
K, a vector giving, for each BIN file, the number of regenerative doses in the SAR protocol;
Nb_measurement, a vector giving, for each BIN file, the number of measurements;
** How to save this list **
You can save this list in a .RData object. To do this, you can use the function save
.
Then, to load this list you can use the function load
(see example section fore more details).
Christophe, C., Kreutzer, S., Philippe, A., Guérin, G., 2023. Generate_DataFile_MG(): Generates, from one (or several) BIN file(s) of Multigrain OSL measurements a list of luminescence data and information before statistical analysis. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
The function imports only BIN/BINXfile records which have been previously selected.
Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guérin
read_BIN2R
, combine_DataFiles
, LT_RegenDose
Age_Computation
, AgeS_Computation
, Palaeodose_Computation
path < system.file("extdata/FER1", "", package="BayLum")
folder < ""
# give the number of sample
nbsample < 1
DATA < Generate_DataFile_MG(
Path = path,
FolderNames = folder,
Nb_sample = nbsample)
str(DATA)
# to save information in RData object in folder containing bin file
#save(DATA,file=c(paste(path,folder,'DATA.RData',sep="")))
# to load information containing DATA.RData object
#load(file=c(paste(path,folder,"DATA.RData",sep="")))
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in AtmosphericNorth_CalC14 allows a calibration for midlatitude Northern Hemisphere atmospher reservoir.
data("IntCal13")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
data(IntCal13)
## maybe str(IntCal13) ; head(IntCal13) ...
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in AtmosphericNorth_CalC14 allows a calibration for midlatitude Northern Hemisphere atmospher reservoir.
data("IntCal20")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
Reimer, P., Austin, W., Bard, E., Bayliss, A., Blackwell, P., Bronk Ramsey, C., . . . Talamo, S. (2020). The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0–55 cal kBP). Radiocarbon, 62(4), 725757. doi:10.1017/RDC.2020.41
data(IntCal20)
## maybe str(IntCal20) ; head(IntCal20) ...
This function plots Lx/Tx
values as a function of regenerative dose,
for every selected aliquot and for each sample.
LT_RegenDose(
DATA,
Path,
FolderNames,
SampleNames = FolderNames,
Nb_sample,
BinPerSample = rep(1, Nb_sample),
SG = rep(TRUE, Nb_sample),
sepDP = c(","),
nrow = 3L,
ncol = nrow
)
DATA 
list (required): list of objects 
Path 
character (required): path to the project folder
(the same as the one used in Generate_DataFile or Generate_DataFile_MG to provide 
FolderNames 
character (required): vector of names of the subfolders
containing the BINfiles, which were used by Generate_DataFile or Generate_DataFile_MG
to generate the 
SampleNames 
character (with default): Names of samples. To use if there is more than one bin file per sample. 
Nb_sample 
integer (required): ID number (in 
BinPerSample 
integer (with default): integer vector (with default):
vector with the number of BIN files per sample, which was used in Generate_DataFile or
Generate_DataFile_MG to generate the 
SG 
logical (with default): vector to set the type of measurement for
each sample 
sepDP 
character (with default): column separator in the 
nrow 
integer (with default): controls the arrangement of the plots,
here the number of rows. Can be set to 
ncol 
integer (with default): controls the arrangement of the plots,
here the number of columns. Can be set to 
To fill FolderNames
and BinPerSample
, we refer to the Detail section from the
Generate_DataFile or Generate_DataFile_MG function.
As well for a precise description of input DATA
.
Lx/Tx plots; there are as many plots as selected aliquots in the DiscPos.csv
file.
There are 9 plots per page.
There is not interpolation.
Christophe, C., Kreutzer, S., Philippe, A., Guérin, G., 2023. LT_RegenDose(): Plots Lx/Tx as a function of the regenerative dose. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guérin
Generate_DataFile, Generate_DataFile_MG
## load data file generated by the function Generate_DataFile
data(DATA3,envir = environment())
path< system.file("extdata/FER1", "", package="BayLum")
folder=""
samplename < "FER1"
LT_RegenDose(
DATA = DATA3,
Path = path,
FolderNames = folder,
SampleNames = samplename,
Nb_sample = 1,
SG = FALSE)
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in marine_CalC14 allows a calibration for hypothetical "global" marine reservoir.
data("Marine13")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
data(Marine13)
## maybe str(Marine13) ; head(Marine13) ...
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in marine_CalC14 allows a calibration for hypothetical "global" marine reservoir.
data("Marine20")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
Heaton, T., Köhler, P., Butzin, M., Bard, E., Reimer, R., Austin, W., . . . Skinner, L. (2020). Marine20—The Marine Radiocarbon Age Calibration Curve (0–55,000 cal BP). Radiocarbon, 62(4), 779820. doi:10.1017/RDC.2020.68
data(Marine20)
## maybe str(Marine20) ; head(Marine20) ...
MCMC samples from the posterior distribution of "A" for age, "D" for palaeodose and "sD" for dispersion of equivalent doses around "D", of the data set GDB5.
data("MCMCsample")
It is a matric with 6000 row and tree column.
A
The first column of the matrice are sampled from the posterior distribution of the paramete A
D
The first column of the matrice are sampled from the posterior distribution of the paramete D
sD
The first column of the matrice are sampled from the posterior distribution of the paramete sD
Tribolo, C., Asrat, A., Bahain, J. J., Chapon, C., Douville, E., Fragnol, C., Hernandez, M., Hovers, E., Leplongeon, A., Martin, L, Pleurdeau, D, Pearson, O , Puaud, S, Assefa, Z. (2017). Across the Gap: Geochronological and Sedimentological Analyses from the Late PleistoceneHolocene Sequence of Goda Buticha, Southeastern Ethiopia. PloS one, 12(1), e0169418.
data(MCMCsample)
## maybe str(MCMCsample) ; plot(MCMCsample[,1],type="l") ...
Age_Computation
A list of JAGS models use to a Bayesian analysis of OSL age of one sample. There are models for various growth curves and various distrubution to describe equivalent dose distribution around the palaeodose.
data("Model_Age")
This list contains:
AgeMultiBF_EXPLIN
a list of 4 models that all consider a saturating exponential plus linear growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgeMultiBF_EXP
a list of 4 models that all consider a saturating exponential growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgeMultiBF_EXPZO
a list of 4 models that all consider a saturating exponential plus linear growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgeMultiBF_EXPLINZO
a list of 4 models that all consider a saturating exponential growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
The different distibutions to describe equivalent dose values around the palaeodose are:
cauchy
a Cauchy distribution with postition parameter equal to the palaeodose of the sample
gaussian
a Gaussian distribution with mean equal to the palaeodose of the sample
lognormal_A
a lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
a lognormal distribution with Median equal to the palaeodose of the sample
For more information we refer to the function Age_Computation
, section Details.
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 125. Technische Universit at Wien, Austria.
Plummer, M. (2015). JAGS Version 4.0. 0 user manual.
data(Model_Age)
## Terminal print
## The JAGS model for a saturating exponential plus linear growth
## (a function of the type \code{f(x)=a(1exp(x/b))+cx+d})
## and a gaussian distribution of equivalent doses around the palaeodose:
writeLines(Model_Age$AgeMultiBF_EXPLIN$cauchy)
AgeC14_Computation
A list of JAGS models use to a Bayesian analysis of C14 calibration age of various sample. Stratigraphic relations can be taken in count to calibrate C14 ages. This ages take into account that some data can be an outlier.
data("Model_AgeC14")
This list contains:
full
a model considering error on calibration curve.
naive
a model not considering error on calibration curve.
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
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, 050000 years cal BP. Radiocarbon 55(4):18891903
data(Model_AgeC14)
writeLines(Model_AgeC14$full)
AgeS_Computation
A list of JAGS models use to a Bayesian analysis of OSL age of various samples. There are models for various growth curves and various distrubution to describe equivalent dose distribution around the palaeodose.
data("Model_AgeS")
This list contains:
AgesMultiCS2_EXPLIN
a list of 4 models that all consider a saturating exponential plus linear growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXP
a list of 4 models that all consider a saturating exponential growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXPZO
a list of 4 models that all consider a saturating exponential plus linear growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXPLINZO
a list of 4 models that all consider a saturating exponential growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
The different distibutions to describe equivalent dose values around the palaeodose are:
cauchy
a Cauchy distribution with postition parameter equal to the palaeodose of the sample
gaussian
a Gaussian distribution with mean equal to the palaeodose of the sample
lognormal_A
a lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
a lognormal distribution with Median equal to the palaeodose of the sample
For more information we refer to the function AgeS_Computation
, section Details.
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 125. Technische Universit at Wien, Austria.
Plummer, M. (2015). JAGS Version 4.0. 0 user manual.
data(Model_AgeS)
## The JAGS model for a saturating exponential plus linear growth
## (a function of the type \code{f(x)=a(1exp(x/b))+cx+d})
## and a gaussian distribution of equivalent doses around the palaeodose:
writeLines(Model_AgeS$AgesMultiCS2_EXP$gaussian)
Palaeodose_Computation
A list of JAGS models use to a Bayesian analysis of OSL palaeodose of one or various samples. There are models for various growth curves and various distrubution to describe equivalent dose distribution around the palaeodose.
data("Model_Palaeodose")
This list contains:
PalaeodosesMultiBF_EXPLIN
a list of 4 models that all consider a saturating exponential plus linear growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
PalaeodosesMultiBF_EXP
a list of 4 models that all consider a saturating exponential growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
PalaeodosesMultiBF_EXPZO
a list of 4 models that all consider a saturating exponential plus linear growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
PalaeodosesMultiBF_EXPLINZO
a list of 4 models that all consider a saturating exponential growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
The different distibutions to describe equivalent dose values around the palaeodose are:
cauchy
a Cauchy distribution with postition parameter equal to the palaeodose of the sample
gaussian
a Gaussian distribution with mean equal to the palaeodose of the sample
lognormal_A
a lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
a lognormal distribution with Median equal to the palaeodose of the sample
For more information we refer to the function Palaeodose_Computation
, section Details.
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 125. Technische Universit at Wien, Austria.
Plummer, M. (2015). JAGS Version 4.0. 0 user manual.
data(Model_Palaeodose)
writeLines(Model_Palaeodose$PalaeodosesMultiBF_EXPLIN$gaussian)
Age_OSLC14
A list of models for C14 data to define likelyhood in JAGS models.
data("ModelC14")
This list contains:
full
a model considering error on calibration curve.
naive
a model not considering error on calibration curve.
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
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, 050000 years cal BP. Radiocarbon 55(4):18891903
data(Model_AgeC14)
writeLines(Model_AgeC14$full)
Age_OSLC14
A list of models for OSL data to define likelyhood in JAGS models.
data("ModelOSL")
This list contains:
AgesMultiCS2_EXPLIN
a list of 4 models that all consider a saturating exponential plus linear growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXP
a list of 4 models that all consider a saturating exponential growth. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXPZO
a list of 4 models that all consider a saturating exponential plus linear growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
AgesMultiCS2_EXPLINZO
a list of 4 models that all consider a saturating exponential growth and fitting through the origin. These 4 models have different distribution to describe equivalent dose values around the palaeodose.
The different distibutions to describe equivalent dose values around the palaeodose are:
cauchy
a Cauchy distribution with postition parameter equal to the palaeodose of the sample
gaussian
a Gaussian distribution with mean equal to the palaeodose of the sample
lognormal_A
a lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
a lognormal distribution with Median equal to the palaeodose of the sample
For more information we refer to the function AgeS_Computation
, section Details.
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 125. Technische Universit at Wien, Austria.
Plummer, M. (2015). JAGS Version 4.0. 0 user manual.
data(ModelOSL)
## The JAGS model of the likelyhood for a saturating exponential plus linear growth
## (a function of the type \code{f(x)=a(1exp(x/b))+cx+d})
## and a gaussian distribution of equivalent doses around the palaeodose:
writeLines(ModelOSL$AgesMultiOSL_EXPLIN$gaussian)
Age_OSLC14
A list to define prior in JAGS models, taking acount OSL data and C14 data in stratigraphic constraint. The difficulty is in the fact that each cases is different. The youngest sample can be a C14 as well as a OSL sample. To resolve this problem we consider diferent cases thanks to this list.
data("ModelPrior")
This list contains:
Sample1_C14
model considering that the youngest sample is a C14 sample
Sample1_OSL
model considering that the youngest sample is a OSL sample
C14_OSL
model considering that the second sample is a C14 sample
OSL_C14
model considering that the second sample is a OSL sample
C14
model considering that the last sample is a C14 sample
OSL
model considering that the last sample is a OSL sample
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing, volume 124, page 125. Technische Universit at Wien, Austria.
Plummer, M. (2015). JAGS Version 4.0. 0 user manual.
data(ModelPrior)
## ModelPrior[[OSL]]
writeLines(ModelPrior$OSL)
This function computes the palaeodose (in Gy) of one or various samples according to the model developed in Combes et al (2015),
based on an output of Generate_DataFile
or Generate_DataFile_MG
or both of them using combine_DataFiles
.
Samples, for which data is avalilable in several BIN files, can be analysed.
Singlegrain or Multigrain OSL measurements can be analysed simultaneouly.
Palaeodose_Computation(
DATA,
SampleNames,
Nb_sample,
BinPerSample = rep(1, Nb_sample),
SavePdf = FALSE,
OutputFileName = c("MCMCplot"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(""),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
Iter = 50000,
t = 5,
n.chains = 3
)
DATA 
list of objects: LT, sLT, ITimes, dLab, ddot_env, regDose, J, K, Nb_measurement,
provided by 
SampleNames 
character vector: names of sample. The length of this vector is equal to 
Nb_sample 
integer: number of samples. 
BinPerSample 
integer vector (with default): vector with the number of BIN files per sample.
The length of this vector is equal to 
SavePdf 
boolean (with default): if TRUE save graph in pdf file named 
OutputFileName 
character (with default): name of the pdf files that will be generated by the function. 
OutputFilePath 
character (with default): path to the pdf files that will be generated by the function. 
SaveEstimates 
boolean (with default): if TRUE save Bayes estimates and credible interval at level 68
in a csv table named 
OutputTableName 
character (with default): name of the table that will be generated by the function if 
OutputTablePath 
character (with default): path to the table that will be generated by the function if 
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. Please see details for more informations on the proposed dose response curves. 
Origin_fit 
logical (with default): if TRUE, forces the dose response curves to pass through the origin. Please see details for more informations on the proposed growth curves. 
distribution 
character (with default): type of distribution that defines how individual equivalent dose values are distributed around the palaeodose. Allowed inputs are "cauchy", "gaussian", "lognormal_A" and "lognormal_M". 
Iter 
integer (with default): number of iterations for the MCMC computation (for more information see 
t 
integer (with default): 1 every 
n.chains 
integer (with default): number of independent chains for the model (for more information see 
** Option on growth curves **
As for Age_Computation
and AgeS_Computation
, the user can choose from 4 dose response curves:
Saturating exponential plus linear growth (PalaeodosesMultiBF_EXPLIN
):
for all x
in IR+, f(x)=a(1exp(x/b))+cx+d
; select
LIN_fit=TRUE
Origin_fit=FALSE
Saturating exponential growth (PalaeodosesMultiBF_EXP
):
for all x
in IR+, f(x)=a(1exp(x/b))+d
; select
LIN_fit=FALSE
Origin_fit=FALSE
Saturating exponential plus linear growth and fitting through the origin (PalaeodosesMultiBF_EXPLINZO
):
for all x
in IR+, f(x)=a(1exp(x/b))+cx
; select
LIN_fit=TRUE
Origin_fit=TRUE
Saturating exponential growth and fitting through the origin (PalaeodosesMultiBF_EXPZO
):
for all x
in IR+, f(x)=a(1exp(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 lognormal distribution with mean or Average equal to the palaeodose of the sample
lognormal_M
: a lognormal distribution with Median equal to the palaeodose of the sample
NUMERICAL OUTPUT
A list containing the following objects:
Sampling that corresponds to a sample of the posterior distributions of palaeodose and equivalent dose dispersion parameters (both in Gy).
Model_GrowthCurve, stating which dose response fitting option was chosen;
Distribution, stating which distribution was chosen to model the dispersion of individual equivalent dose values around the palaeodose of the sample.
The Gelman and Rubin test of convergency: prints the result of the Gelman and Rubin test of convergency for palaeodose and equivalent dose dispersion parameters for each sample.
A result close to one is expected.
In addition, the user must visually assess the convergency of the trajectories by looking at the pdf file
generated by the function (see PLOT OUTPUT for more informations).
If both convergencies (Gelman and Rubin test and plot checking) are satisfactory,
the user can consider the printed estimates as valid. Otherwise, the user may try increasing the number of MCMC interations
(Iter
) to reach convergency.
Credible intervals and Bayes estimates: prints the Bayes esitmates, the credible intervals at 95 the palaeodose and equivalent dose dispersion parameters for each sample.
PLOT OUTPUT
MCMC trajectories
A graph with the MCMC trajectories and posterior distributions of the palaeodose and equivalent dose dispersion parameters
is displayed, there is one page per sample.
The first line of the figure correponds to the palaeodose parameter and the second to the equivalent dose dispersion parameter.
On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
Summary of palaeodose estimates: plot credible intervals and Bayes estimate of each sample palaeodose on a same graph.
To give result in a publication, we recommend to give the Bayes estimate of the parameters as well as the credible interval at 95
Christophe, C., Kreutzer, S., Philippe, A., Guérin, G., 2023. Palaeodose_Computation(): Bayesian analysis for the palaeodose estimation of various samples. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
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.
Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guérin
Combes, B., Philippe, A., Lanos, P., Mercier, N., Tribolo, C., Guerin, G., Guibert, P., Lahaye, C., 2015. A Bayesian central equivalent dose model for optically stimulated luminescence dating. Quaternary Geochronology 28, 6270. doi:10.1016/j.quageo.2015.04.001
Generate_DataFile
, Generate_DataFile_MG
, combine_DataFiles
,
rjags
, plot_MCMC
,
Age_Computation
, AgeS_Computation
## Load data
data(DATA1,envir = environment())
## Palaeodose computation of samples GDB3
P=Palaeodose_Computation(DATA=DATA1,Nb_sample=1,SampleNames=c("GDB5"),Iter=100)
Create Age Plot
plot_Ages(
object,
sample_names = NULL,
sample_order = NULL,
plot_mode = "ages",
...
)
object 
list or data.frame (required): Output as created by functions like AgeC14_Computation, which
is a list of class 
sample_names 
character (optional): alternative sample names used for the plotting. If the length of the provided character vector is shorter than the real number of samples, the names are recycled. 
sample_order 
numeric (optional): argument to rearrange the sample order, e.g., 
plot_mode 
character (with default): allows to switch from displaying ages as points with lines ( 
... 
further arguments to control the plot output,
standard arguments are: 
This function creates an age plot showing the mean ages along with the credible intervals. The function provides various arguments to modify the plot output, however, for an ultimate control the function returns the data.frame extracted from the input object for own plots.
The function returns a plot and the data.frame used to display the data
0.1.5
Kreutzer, S., Christophe, C., 2023. plot_Ages(): Create Age Plot. Function version 0.1.5. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, Institute of Geography, RuprechtKarlUniversity of Heidelberg (Germany), based on code written by Claire Christophe
AgeC14_Computation, AgeS_Computation
## load data
data(DATA_C14,envir = environment())
C14Cal < DATA_C14$C14[,1]
SigmaC14Cal < DATA_C14$C14[,2]
Names < DATA_C14$Names
nb_sample < length(Names)
## Age computation
Age < AgeC14_Computation(
Data_C14Cal = C14Cal,
Data_SigmaC14Cal = SigmaC14Cal,
SampleNames = Names,
Nb_sample = nb_sample,
PriorAge = rep(c(20,60),nb_sample),
Iter = 500,
quiet = TRUE)
## plot output
plot_Ages(Age)
## plot output
plot_Ages(Age, plot_mode = "density")
This function uses the output of rjags::jags.model to visualise the traces of the MCMC and the
corresponding densities. In particular it displays the posterior distributions of the age, if it is calculated,
palaeodose and the equivalent dose dispersion parameters of the sample. The function output is very
similar to plot output produced with the 'coda' package, but tailored to meet the needs in
the context of the 'BayLum'
package.
plot_MCMC(
object,
sample_names = NULL,
variables = c("A", "D", "sD"),
axes_labels = c(A = "Age (ka)", D = "D (Gy)", sD = "sD (Gy)"),
n.chains = NULL,
n.iter = 1000,
smooth = FALSE,
rug = TRUE,
plot_single = FALSE,
...
)
object 
coda::mcmc.list or coda::mcmc (required): Output generated by rjags::jags.model,
e.g., in Age_Computation. Alternatively, limited support is provided for 
sample_names 
character (optional): Names of the used samples. This argument overrides the optional
argument 
variables 
character (with default): Variables in your coda::mcmc object to be plotted. 
axes_labels 
character (with default): Axes labels used for the trace and density plots. The labels should be provided as named character vector with the parameter names as the names used to assign the axes labelling. The labelling for the xaxis (trace plots) and yaxis (density plot) cannot be modified. 
n.chains 
numeric (optional): Set the number of chains to visualise, if nothing is provided the number of chains is determined from the input object 
n.iter 
integer (with default): Set the number of iterations to be visualised in the trace plots, regardless
of the size of the input dataset as long as the real number of iterations is > 
smooth 
logical (with default): Enable/disables smooth of trace plots using stats::smooth 
rug 
logical (with default): Enable/disables rug under density plots 
plot_single 
logical (with default): Enables/disables the single plot mode of the function, i.e.
if set to 
... 
further arguments that can be passed to modify the plot output. Supported arguments are

The function is used in the function Age_Computation, AgeS_Computation and Palaeodose_Computation, but can be used also as standalone plot function.
Two plots: Traces of the MCMC chains and the corresponding density plots. This plots are similar to coda::traceplot and coda::densplot.
0.1.5
Kreutzer, S., Christophe, C., 2023. plot_MCMC(): Plot MCMC trajectories and posterior distributions. Function version 0.1.5. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, Institute of Geography, RuprechtKarl University of Heidelberg (Germany). This function is a rewritten version of the function 'MCMC_plot()' by Claire Christophe
Age_Computation, AgeS_Computation, Palaeodose_Computation, rjags::coda.samples and rjags packages.
data(MCMCsample,envir = environment())
object < coda::as.mcmc(MCMCsample)
plot_MCMC(object)
Create a hexbin plot matrix (hexbin::hexplom) of age results returned by the Bayesian age calculation.
plot_Scatterplots(
object,
variables = c("A"),
sample_names = NULL,
sample_selection = NULL,
n.chains = NULL,
plot_type = "hexbin",
plot_mode = "matrix",
...
)
ScatterSamples(...)
object 
coda::mcmc.list or a data.frame (required): mcmc list objects generated by rjags::jags.model in AgeS_Computation,
AgeC14_Computation or Age_OSLC14. If a data.frame is provided, only the first two columns are taken and 
variables 
character (with default): variable to be selected for the scatter plot, e.g., 
sample_names 
character (optional): sample names shown in the plot matrix 
sample_selection 
numeric (with default): vector of samples to be plotted in the scatter matrix, e.g.,

n.chains 
integer (with default): allows to limit the number of chains shown, by default the results of all chains are plotted. 
plot_type 
character (with default): switch between different plot types, 
plot_mode 
character (with default): switch between a 
... 
further arguments to control the plot output, standard plot arguments supported are 
Additional supported plot arguments
The following table lists additional arguments supported by the function in order to fine tune the
graphical output. Such arguments, can just be added in the function call. Example, for disabling
the graphics::rug in the plot mode smoothScatter
you can type plot_Scatterplots(..., rug = FALSE)
Please note that not all arguments are supported by all plot types.
ARGUMENT  SUPPORTED BY PLOT TYPE  DESCRIPTION 
colramp 
hexbin and smoothScatter 
Option to define an own colour ramp, by defining an own function, e.g., function(n) heat.colors(n, alpha = 1) . 
pscales 
hexbin and smoothScatter 
Controls the number of ticks shown on the plot axes, please note that the number works proportionally. 
bw_smoothScatter 
smoothScatter 
Controls the bandwidth of the smooth scatter, cf. graphics::smoothScatter 
rug 
smoothScatter 
enables/disables rugs 
nlevels 
smoothScatter 
controls the number of isolines shown (cf. graphics::contour) 
nrpoints 
smoothScatter 
defines the number of nrpoints to be plotted graphics::smoothScatter 
col_contour 
smoothScatter 
defines the colour of the contour lines 
col_nrpoints 
smoothScatter 
sets colour of the nrpoints in the scatter plot

A scatter plot based on hexbin::hexplom
0.3.2
Kreutzer, S., Christophe, C., Philippe, A., Guérin, G., 2023. plot_Scatterplots(): Display Scatter Plot Matrix of the Bayesian Age Results. Function version 0.3.2. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Sebastian Kreutzer, Institute of Geography, RuprechtKarl University of Heidelberg (Germany) , based on the function 'ScatterSamples()' by Claire Christophe, Anne Philippe, Guillaume Guérin
Age_Computation, AgeS_Computation, AgeC14_Computation, and rjags packages.
data(AgeS,envir = environment())
##hexbin
plot_Scatterplots(
object = AgeS$Sampling,
sample_names = c("GDB5", "GDB3"),
sample_selection = c(1,2)
)
##scatter smooth (matrix)
plot_Scatterplots(
object = AgeS$Sampling,
sample_names = c("GDB5", "GDB3"),
sample_selection = c(1,2),
plot_type = "smoothScatter")
##scatter smooth (single)
plot_Scatterplots(
object = AgeS$Sampling,
sample_names = c("GDB5", "GDB3"),
sample_selection = c(1,2),
plot_type = "smoothScatter",
plot_mode = "single")
Construct the stratigraphic matrix used in the functions AgeS_Computation and AgeC14_Computation for samples that are all ordered by increasing age.
SC_Ordered(Nb_sample)
Nb_sample 
integer (required): the number of samples 
Stratigraphic matrix where each sample are ordered by increasing order. This matrix can be intergrated in the function AgeS_Computation. Please see AgeS_Computation for more information on this matrix.
Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., 2023. SC_Ordered(): Create stratigraphically ordered sample matrix. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Claire Christophe, Anne Philippe, Sebastian Kreutzer, Guillaume Guérin
# compute the stratigraphic matrix for 3 samples such that the first sample is younger
# than the second, and the second is younger than the third
SC < SC_Ordered(Nb_sample = 3)
This function helps to define the stratigraphic relation between samples, with questions.
The output of this function can be used in function AgeS_Computation
.
SCMatrix(Nb_sample, SampleNames)
Nb_sample 
interger: the sample number. 
SampleNames 
charcater vector: sample names. 
Ask if sample i
is younger than sample j
to construc the stratigraphic constrain matrix.
A Matrix that summarise the ordered relation between samples.
This matrix can be intergrate in AgeS_Computation
function.
We refer to detail on AgeS_Computation
for more information concerning this matrix.
Christophe, C., Philippe, A., 2023. SCMatrix(): Definition of the stratigraphic constraint matrix. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Claire Christophe, Anne Philippe, Guillaume Guerin
AgeS_Computation
## Assume that "sample1" is younger than "sample2"
## That means the expected value is 1.
## It is an interactive function.
## Not run:
SCMatrix(Nb_sample=2,SampleNames=c("sample1","sample2"))
## Enter the value 1
## End(Not run)
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in AtmosphericSouth_CalC14 allows a calibration for midlatitude Southern Hemisphere atmospher reservoir.
data("SHCal13")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years (in ky) befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
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 050000 years cal BP. Radiocarbon 55(4)=18691887.
data(SHCal13)
## maybe str(SHCal13) ; head(SHCal13) ...
As 14C years is not equal to calendar years because atmospheric 14C concentration varies through time. Hence, data in AtmosphericSouth_CalC14 allows a calibration for midlatitude Southern Hemisphere atmospher reservoir.
data("SHCal20")
A data frame with 3 variables.
CAL.BP
a numeric vector correpondig to calendar years (in ky) befor present
X14C.age
a numeric vector correponding to 14C age
Error
a numeric vector correponding to error arround 14C age measurement
Hogg, A., Heaton, T., Hua, Q., Palmer, J., Turney, C., Southon, J., . . . Wacker, L. (2020). SHCal20 Southern Hemisphere Calibration, 0–55,000 Years cal BP. Radiocarbon, 62(4), 759778. doi:10.1017/RDC.2020.59
data(SHCal20)
## maybe str(SHCal20) ; head(SHCal20) ...
This function allows the user to write all .csv files expected by Generate_DataFile and Generate_DataFile_MG. Unlike create_FolderTemplates, this function makes it possible to write .csv files with all information directly from R. No further modification of .csv files are required. The purpose of this function is (i) to reduce tedious manual editing of .csvfiles and the errors that result (ii) to introduce an easy way to review information inside .csvfiles (by revisiting code rather than opening individual .csvfiles) and (iii) to streamline folder and file creation when preparing data to run BayLum's modelling functions. Note: the user will still need to move the appropriate .binfiles into all the sample folders.
write_BayLumFiles(
folder,
SampleNames = "Sample_1",
BinPerSample = rep(1, length(SampleNames)),
SubSampleNames = NULL,
DiscPos = NULL,
DRenv = 1,
DRenv.error = 0.04,
DRsource = 0.1,
DRsource.error = 0.002,
signal.integral.min = 6,
signal.integral.max = 10,
background.integral.min = 346,
background.integral.max = 395,
inflatePercent = 0.025,
nbOfLastCycleToRemove = 2
)
folder 
character (required*): The name of the main folder in which all subsequent BayLum files and folders will be located. This could be a path to an already existing folder, or the path/name of a folder to be created. 
SampleNames 
character (required): Vector of sample names. 
BinPerSample 
numeric (with default): Vector of numbers indicating the number of .binfiles per sample. 
SubSampleNames 
character (optional): Vector of names to give each subfolder within a sample when the number of .binfiles in a sample counts more than one. If omitted or NULL, the subfolders are named by the subfolder count number. 
DiscPos 
numeric (with default): List of data frames with each data frame having one or two columns to identify aliquots/grains to be included in the analysis. The first column corresponds to the position number, and the second column corresponds to the grain number. If the data frame has only one column, a Disc.csv will be written. If the data.frame has two columns, a DiscPos.csv will be written. The length of the list should be the number of .binfiles included. 
DRenv 
numeric (with default): Vector where 
DRenv.error 
numeric (with default): Vector where 
DRsource 
numeric (with default): Vector where 
DRsource.error 
numeric (with default): Vector where 
signal.integral.min 
numeric (with default): Vector where 
signal.integral.max 
numeric (with default): Vector where 
background.integral.min 
numeric (with default): Vector where 
background.integral.max 
numeric (with default): Vector where 
inflatePercent 
numeric (with default): Vector where 
nbOfLastCycleToRemove 
numeric (with default): Vector where 
The function returns nothing, but writes the folder structure.
0.1.0
Baumgarten, F.H., 2023. write_BayLumFiles(): Write BayLum .csvfiles. Function version 0.1.0. In: Christophe, C., Philippe, A., Kreutzer, S., Guérin, G., Baumgarten, F.H., 2023. BayLum: Chronological Bayesian Models Integrating Optically Stimulated. R package version 0.3.1.900017. https://CRAN.rproject.org/package=BayLum
Frederik Baumgarten, RadPhys, DTU Physics, Technical University of Denmark (Denmark)
Generate_DataFile, Generate_DataFile_MG
# example samples
SampleNames < c("OSL1MG","OSL2SG")
# number of .binfiles for each sample
BinPerSample < c(1,3)
# List of data.frames of accepted aliquot/grain to be included
# in the analysis for each .binfile.
DiscPos < list(
data.frame("position" = 1:48),
data.frame("position" = c(1,1,1,1), "grain" = c(4,67,92,99)),
data.frame("position" = c(2,2,2,2), "grain" = c(7,13,41,72)),
data.frame("position" = c(3,3,3,3), "grain" = c(7,52,67,88)))
# example 1: write to disk (all together)
write_BayLumFiles(
folder = paste(tempdir(),"new_BayLum_folder",sep = "/"),
SampleNames = SampleNames,
BinPerSample = BinPerSample,
DiscPos = DiscPos,
DRenv = c(1.75, 1.52, 1.52, 1.52),
DRenv.error = c(0.04, 0.03, 0.03, 0.03),
DRsource = c(0.2075, 0.1501, 0.1501, 0.1501),
DRsource.error = c(0.0010, 0.0008, 0.0008, 0.0008))
# example 2: write to disk (by sample)
write_BayLumFiles(
folder = paste(tempdir(),"new_BayLum_folder",sep = "/"),
SampleNames = "OSL1MG",
BinPerSample = 1,
DiscPos = DiscPos[[1]],
DRenv = 1.75,
DRenv.error = 0.04,
DRsource = 0.2075,
DRsource.error = 0.0010)
write_BayLumFiles(
folder = paste(tempdir(),"new_BayLum_folder",sep = "/"),
SampleNames = "OSL2SG",
BinPerSample = 3,
DiscPos = DiscPos[2:4],
DRenv = 1.75,
DRenv.error = 0.04,
DRsource = 0.2075,
DRsource.error = 0.0010)