`R/AgeS_Computation.R`

`AgeS_Computation.Rd`

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.

Single-grain or Multi-grain OSL measurements can be analysed simultaneously.

```
AgeS_Computation(
DATA,
SampleNames = DATA$SampleNames,
Nb_sample = DATA$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:`LT`

,`sLT`

,`ITimes`

,`dLab`

,`ddot_env`

,`regDose`

,`J`

,`K`

,`Nb_measurement`

, provided by the function Generate_DataFile, Generate_DataFile_MG or combine_DataFiles.`DATA`

contains informations for more than one sample. If there is stratigraphic relations between samples, informations in DATA must be ordered by order of increasing ages. See the details section to for more informations. (2): An object of class`BayLum.list`

which is provided by the output of AgeS_Computation. When input of class`BayLum.list`

is identified, no new JAGS model is created. Instead, the JAGS model specified by the AgeS_Computation output is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.- SampleNames
character vector: names of samples. The length of this vector is equal to

`Nb_sample`

.- Nb_sample
integer: number of samples,

`Nb_sample>1`

.- PriorAge
numeric vector (with default): lower and upper bounds for age parameter of each sample (in ka). Note that,

`length(PriorAge)=2*Nb_sample`

and`PriorAge[2i-1,2i]`

corresponds to the lower and upper bounds of sample whose number ID is equal to`i`

.- BinPerSample
integer vector (with default): vector with the number of BIN files per sample. The length of this vector is equal to

`Nb_sample`

.`BinPerSample[i]`

corresponds to the number of BIN files for the sample whose number ID is equal to`i`

. For more information to fill this vector, we refer to details in Generate_DataFile or in Generate_DataFile_MG.- SavePdf
logical (with default): if TRUE save graphs in pdf file named

`OutputFileName`

in folder`OutputFilePath`

.- OutputFileName
character (with default): name of the pdf file that will be generated by the function if

`SavePdf = TRUE`

;`length(OutputFileName)=2`

, see**PLOT OUTPUT**in**Value**section for more informations.- OutputFilePath
character (with default): path to the pdf file that will be generated by the function if

`SavePdf`

=TRUE. If it is not equal to "", it must be terminated by "/".- SaveEstimates
logical (with default): if TRUE save Bayes' estimates, credible interval at level 68% and 95% and the result of the Gelman en Rubin test of convergence, in a csv table named

`OutputFileName`

in folder`OutputFilePath`

.- OutputTableName
character (with default): name of the table that will be generated by the function if

`SaveEstimates = TRUE`

.- OutputTablePath
character (with default): path to the table that will be generated by the function if

`SaveEstimates = TRUE`

. If it is not equal to "", it must be terminated by "/".- 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

`THETA`

; the user can refer to a matrix (numeric matrix) or to a csv file (character). Otherwise, default value is suitable, and only individual errors are considered.- sepTHETA
character (with default): if

`THETA`

is character, indicate column separator in`THETA`

CSV-file.- 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

`StratiConstraints`

; the user can refer to a matrix (numeric matrix) or to a csv file (character). If there is no stratigraphic relation default value is suitable.- sepSC
character (with default): if

`StratiConstraints`

is character, indicate column separator in`StratiConstraints`

.csv file.- LIN_fit
logical (with default): if TRUE (default) allows a linear component, on top of the (default) saturating exponential curve, for the fitting of dose response curves. 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

`"cauchy"`

,`"gaussian"`

,`"lognormal_A"`

and`"lognormal_M"`

, see details section for more informations.- 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`distribution`

is ignored and no safety net is applied. If the function crashes it is up to the user.- 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

`t`

iterations of the MCMC is considered for sampling the posterior distribution. (for more information see runjags::run.jags).- 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

`"rjags"`

(the default) and`"rjparallel"`

have been tested. (for more information about these possibilities and others, see runjags::run.jags)- autorun
logical (with default): choose to automate JAGS processing. JAGS model will be automatically extended until convergence is reached (for more information see runjags::autorun.jags).

- quiet
logical (with default): enables/disables

`rjags`

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

- ...
further arguments that can be passed to control the Bayesian process. 1) When creating a new JAGS model, see details for supported arguments. 2) If extending a JAGS model see arguments of runjags::extend.JAGS.

**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%.

**Supported ... arguments**

ARGUMENT | INPUT | METHOD | DEFAULT | DESCRIPTION |

`max.time` | character | `rjparallel` | `Inf` | maximum allowed processing time, e.g.,
`10m` for 10 minutes (cf. runjags::autorun.jags) |

`interactive` | logical | `rjparallel` | `FALSE` | enable/disable interactive mode (cf. runjags::autorun.jags) |

`startburnin` | integer | `rjparallel` | `4000` | number of burn-in iterations (cf. runjags::autorun.jags) |

`startsample` | integer | `rjparallel` | `10000` | total number of samples to assess convergence (cf. runjags::autorun.jags) |

`inits` | named list | `rjparallel` | `NA` | fine control over random seeds and random number generator runjags::autorun.jags |

**How to fill StratiConstraints**

If there 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[2i-1]`

) for the sample whose number ID is equal to`i`

, is taken into account.**Sample relations**: for all`j in {2,...,Nb_Sample+1}`

and all`i in {j,...,Nb_Sample}`

,`StratiConstraints[j,i]=1`

if sample age whose number ID is equal to`j-1`

is lower than sample age whose number ID is equal to`i`

. Otherwise,`StratiConstraints[j,i]=0`

.

Note that `StratiConstraints_{2:Nb_sample+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(1-exp(-x/b))+cx+d`

; select`LIN_fit=TRUE`

`Origin_fit=FALSE`

**Saturating exponential growth**(`AgesMultiCS2_EXP`

):for all

`x`

in IR+,`f(x)=a(1-exp(-x/b))+d`

; select`LIN_fit=FALSE`

`Origin_fit=FALSE`

**Saturating exponential plus linear growth and fitting through the origin**(`AgesMultiCS2_EXPLINZO`

):for all

`x`

in IR+,`f(x)=a(1-exp(-x/b))+cx`

; select`LIN_fit=TRUE`

`Origin_fit=TRUE`

**Saturating exponential growth and fitting through the origin**(`AgesMultiCS2_EXPZO`

):for all

`x`

in IR+,`f(x)=a(1-exp(-x/b))`

; select`LIN_fit=FALSE`

`Origin_fit=TRUE`

**Option on equivalent dose distribution around the palaeodose**

The use can choose between :

`cauchy`

: a Cauchy distribution with location parameter equal to the palaeodose of the sample;`gaussian`

: a Gaussian distribution with mean equal to the palaeodose of the sample;`lognormal_A`

: a log-normal distribution with mean or**A**verage equal to the palaeodose of the sample:`lognormal_M`

: a log-normal distribution with**M**edian equal to the palaeodose of the sample.

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

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

Combes, Benoit and Philippe, Anne, 2017. Bayesian analysis of multiplicative Gaussian error for multiple ages estimation in optically stimulated luminescence dating. Quaternary Geochronology (39, 24-34)

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, 62-70. 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)
if (FALSE) {
## 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")
}
```