R/Age_OSLC14.R
In R-Lum/BayLum: Chronological Bayesian Models Integrating Optically Stimulated Luminescence and Radiocarbon Age Dating
Defines functions
Age_OSLC14
Documented in
Age_OSLC14
#' @title Bayesian analysis for age estimation of OSL measurements and C-14 ages of various samples
#'
#' @description This function computes an age of OSL data consisting of at least two samples and calibrate
#' C-14 ages of samples to get an age (in ka).\cr
#' Ages of OSL data are computed according to the model given in Combes and Philippe (2017).
#' Single-grain or Multi-grain OSL measurements can be analysed simultaneously
#' (with output of [Generate_DataFile] or [Generate_DataFile_MG] or both of them using [combine_DataFiles]).
#' Samples, for which data is available in several BIN files, can be analysed.\cr
#' For C-14 data, the user can choose one of the following radiocarbon calibration curve:
#' Northern or Southern Hemisphere or marine atmospheric.
#'
#' @param DATA Two types of inputs are supported.
#' (1): a list of objects: `LT`, `sLT`, `ITimes`, `dLab`, `ddot_env`, `regDose`, `J`,
#' `K`, `Nb_measurement`, provided by the function [Generate_DataFile], [Generate_DataFile_MG] or [combine_DataFiles].
#' \code{DATA} contains information for more than one sample.
#' If there is stratigraphic relations between samples, informations in DATA must be ordered by order of
#' increasing ages. See the details section to for more informations.
#' (2): an object of class `BayLum.list` which is provided by the output of [Age_OSLC14]. When input of class `BayLum.list` is identified, no new JAGS model is created. Instead, the JAGS model specified within the `BayLum.list` is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.
#'
#' @param Data_C14Cal [numeric] vector: corresponding to C-14 age estimate
#' (in years, conversion in ka is automatically done in the function).
#' If there is stratigraphic relations between samples, `Data_C14Cal` must be ordered
#' by order of increasing ages.
#'
#' @param Data_SigmaC14Cal [numeric]: corresponding to the error of C-14 age estimates.
#'
#' @param Nb_sample [numeric]: number of samples (OSL data and C-14 age),
#' (`Nb_sample>3`, at least to sample of OSL data and one sample of C-14 age).
#'
#' @param SampleNames [character]: sample names for both OSL data and C14 data.
#' The length of this vector is equal to `Nb_sample`.
#' If there is stratigraphic relation, this vector must be ordered by increasing
#' order (to mix OSL samples and C-14 ages if it is needed).
#'
#' @param SampleNature [matrix]: states the nature of the sample.
#' Row number of `SampleNature` matrix is equal to `2` and column number is equal
#' to `Nb_sample`. First line of the matrix: `SampleNature[1,i]` states if sample
#' whose number ID is equal to `i`, is an OSL sample `1` or not `0`.
#' Second line of the matrix: `SampleNature[2,i]` states if sample whose number ID
#' is equal to `i`, is an C-14 sample `1` or not `0.
#'
#' @param PriorAge [numeric] (with default): lower and upper bounds for age parameter of each sample **(in ka)**.
#' Note that, `length(PriorAge) = 2*Nb_sample` sand `PriorAge[2i-1,2i]` corresponds
#' to the lower and upper bounds of sample whose number ID is equal to `i`.
#'
#' @param SavePdf [logical] (with default): if `TRUE` save graphs in pdf file
#' named `OutputFileName` in folder `OutputFilePath`.
#'
#' @param OutputFileName [character] (with default): name of the pdf file that
#' will be generated by the function if `SavePdf=TRUE`, `length(OutputFileName)=3`,
#' see **PLOT OUTPUT** in **Value** section for more informations.
#'
#' @param 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 "/".
#'
#' @param 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`.
#'
#' @param OutputTableName [character] (with default): name of the table that will
#' be generated by the function if `SaveEstimates=TRUE`.
#'
#' @param 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 "/".
#'
#' @param StratiConstraints [matrix] or [character] (with default): input object for the
#' stratigraphic relation between samples. If there is stratigraphic relation between
#' samples see the details section for instructions regarding how to correctly fill `StratiConstraints`,
#' the user can refer to a [matrix] or to a CSV-file [character]. Otherwise, default value is suitable.
#'
#' @param sepSC [character] (with default): if `StratiConstraints` is character,
#' indicate column separator in `StratiConstraints` CSV-file.
#'
#' @param BinPerSample [numeric] (with default): vector with the number of BIN-files
#' per OSL sample. The length of this vector is equal to the number of OSL samples.
#' `BinPerSample[i]` corresponds to the number of BIN files for the sample whose
#' number ID is equal to `i`. For more information to fill this vector, we refer
#' to details in [Generate_DataFile] or in [Generate_DataFile_MG].
#'
#' @param THETA numeric [matrix] or [character] (with default): input object for
#' systematic and individual error for OSL samples. If systematic errors are considered,
#' see the details section for instructions regarding how to correctly fill `THETA`;
#' the user can refer to a matrix (numeric matrix) or to a csv file (character).
#' Otherwise, default value is suitable, and only individual error is considered.
#'
#' @param sepTHETA [character] (with default): if `THETA` is character,
#' indicate column separator in `THETA` CSV-file.
#'
#' @param LIN_fit [logical] (with default): if `TRUE` (default) allows a linear component,
#' on top of the (default) saturating exponential curve, for the fitting of dose response curves,
#' for OSL samples. See details for more informations on the proposed dose response curves.
#'
#' @param Origin_fit p[logical] (with default): if `TRUE`, forces the dose response
#' curves to pass through the origin. See details for more informations on the proposed
#' growth curves, for OSL samples.
#'
#' @param distribution [character] (with default): type of distribution that defines
#' how individual equivalent dose values are distributed around the palaeodose, for OSL samples.
#' Allowed inputs are `"cauchy"`, `"gaussian"`, `"lognormal_A"` and `"lognormal_M"`,
#' see details for more informations.
#'
#' @param Model_C14 [character] (with default): if `"full"`, error on estimate calibration
#' curve is taken account, for C-14 samples. If `"naive"` this error is not taken account in the age estimate.
#'
#' @param CalibrationCurve [character] (with default): calibration curve chosen, for C-14 samples.
#' Allowed inputs are
#' \itemize{
#' \item \bold{"Intcal13"} or \bold{"Intcal13"} for Northern Hemisphere atmospheric radiocarbon calibration curve,
#' \item \bold{"Marine13"} or \bold{"Marine13"} for Marine radiocarbon calibration curve,
#' \item \bold{"SHCal13"} or \bold{"SHCal20"} for Southern Hemisphere atmospheric radiocarbon calibration curve
#' \item \bold{a csv file, with tree columns, the first column is dedicated to "Cal.BP", the second to "XC-14.age", the third to "Error".
#' The decimal of this file must be a dot, and the separator must be a comma. }
#' }
#'
#' @param adapt [integer] (with default): the number of iterations used in the adaptive phase of the simulation (see [runjags::run.jags]]).
#' @param 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]]).
#' @param Iter [integer] (with default): the number of iterations to run and who will be used to assess convergence and ages (see [runjags::run.jags]]).
#'
#' @param t [numeric] (with default): 1 every `t` iterations of the MCMC is
#' considered for sampling the posterior distribution (for more information see [[rjags::jags.model].
#'
#' @param n.chains [numeric] (with default): number of independent chains for the model
#' (for more information see [[rjags::jags.model]).
#'
#' @param jags_method [character] (with default): select which method to use in order to call JAGS, supported are `"rjags"` (the default), `rjparallel`, `simple`, `interruptible`, `parallel`, and `snow` (for more information about each of these possibilities, see [runjags::run.jags]])
#'
#' @param 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]).
#'
#' @param quiet [logical] (with default): enables/disables [rjags] messages
#'
#' @param roundingOfValue [integer] (with default): Integer indicating the number of decimal places to be used, default = 3.
#
#' @param ... further arguments that can be passed to control the Bayesian process, see details
#' for supported arguments
#'
#' @details
#'
#' Note that there are three types of arguments in the previous list.
#' There are arguments for information concerning only OSL samples: \code{DATA}, \code{BinPerSample}, \code{THETA},
#' \code{sepTHETA}, \code{LIN_fit}, \code{Origin_fit}, \code{distribution}.
#'
#' There are arguments for information concerning only C14 samples: \code{Data_C14Cal}, \code{Data_SigmaC14Cal},
#' \code{Model_C14}, \code{CalibrationCurve}.
#'
#' There are arguments for information concerning all the samples: \code{Nb_sample}, \code{SampleNames}, \code{SampleNature},
#' \code{PriorAge}, \code{SavePdf}, \code{OutputFileName}, \code{OutputFilePath}, \code{SaveEstimates}, \code{OutputTableName},
#' \code{OutputTablePath}, \code{StratiConstraints}, \code{sepSC}.\cr
#'
#' **Supported `...` arguments**
#' \tabular{lllll}{
#' ARGUMENT \tab INPUT \tab METHOD \tab DEFAULT \tab DESCRIPTION\cr
#' `max.time` \tab [character] \tab `rjparallel` \tab `Inf` \tab maximum allowed processing time, e.g.,
#' `10m` for 10 minutes (cf. [runjags::autorun.jags])\cr
#' `interactive` \tab [logical] \tab `rjparallel` \tab `FALSE` \tab enable/disable interactive mode (cf. [runjags::autorun.jags])\cr
#' `startburnin` \tab [integer] \tab `rjparallel` \tab `4000` \tab number of burn-in iterations (cf. [runjags::autorun.jags]) \cr
#' `startsample` \tab [integer] \tab `rjparallel` \tab `10000` \tab total number of samples to assess convergence
#' (cf. [runjags::autorun.jags]) \cr
#' `inits` \tab named [list] \tab `rjparallel` \tab `NA` \tab fine control over random seeds and random number generator [runjags::autorun.jags]
#' }
#'
#' **How to fill `StratiConstraints?**\cr
#'
#' If there are stratigraphic relations between samples, \bold{14C estimate age in \code{Data_C14Cal} must be ordered by order of increasing ages,
#' as informations in \code{DATA}}. Names in \code{SampleNames} must be ordered and corresponds to the order in \code{Data_C14Cal} and in \code{DATA},
#' also if it is needed to mix names of OSL samples and 14C samples.
#'
#' The user can fill the \code{StratiConstraints} matrix as follow.
#' \enumerate{
#' \item \bold{Size of the matrix}: row number of \code{StratiConstraints} matrix is equal to \code{Nb_sample+1},
#' and column number is equal to \code{Nb_sample}.
#' \item \bold{First line of the matrix}:
#' for all \code{i in {1,...,Nb_Sample}}, \code{StratiConstraints[1,i]=1} that means the lower bound of the sample age (given in \code{PriorAge[2i-1]})
#' for the sample whose number ID is equal to \code{i}, is taken into account.
#' \item \bold{Sample relations}: for all \code{j in {2,...,Nb_Sample+1}} and all \code{i in {j,...,Nb_Sample}},
#' \code{StratiConstraints[j,i]=1} if sample age whose number ID is equal to \code{j-1} is lower than sample age whose number ID is equal to \code{i}.
#' Otherwise, \code{StratiConstraints[j,i]=0}.
#' }
#' Note that \code{StratiConstraints_{2:Nb_sample+1,1:Nb_sample}} is a upper triangular matrix.
#'
#' The user can also use \code{\link{SCMatrix}} or \code{\link{SC_Ordered}} (if all samples are ordered) function to construct
#' the \code{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 \code{sepSC}.\cr
#'
#' **How to fill `THETA` covariance matrix concerning common and individual error?**\cr
#'
#' If systematic errors are considered, the user can fill the \code{THETA} matrix as follow.
#' \itemize{
#' \item row number of \code{THETA} is equal the column number, equal to \code{Nb_sample}.
#' \item For all \code{i in {1,...,Nb_sample}}, \code{THETA[i,i]} contains individual error
#' plus systematic error of the sample whose number ID is equal to \code{i}.
#' \item For all \code{i,j in {1,...,Nb_sample}} and \code{i} different from \code{j} ,
#' \code{THETA[i,j]} contains common error between samples whose number ID are equal to \code{i} and \code{j}.
#' }
#' Note that \code{THETA[i,j]} is a symmetric matrix.
#'
#' The user can also refer to a .csv file that contains the errors as defined above.\cr
#'
#' **Option on growth curves**\cr
#'
#' As for \code{\link{Age_Computation}} and \code{\link{Palaeodose_Computation}}, the user can choose from 4 dose response curves:
#' \itemize{
#' \item \bold{Saturating exponential plus linear growth} (\code{AgesMultiCS2_EXPLIN}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx+d}; select
#' \itemize{
#' \item \code{LIN_fit=TRUE}
#' \item \code{Origin_fit=FALSE}
#' }
#' \item \bold{Saturating exponential growth} (\code{AgesMultiCS2_EXP}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+d}; select
#' \itemize{
#' \item \code{LIN_fit=FALSE}
#' \item \code{Origin_fit=FALSE}
#' }
#' \item \bold{Saturating exponential plus linear growth and fitting through the origin} (\code{AgesMultiCS2_EXPLINZO}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx}; select
#' \itemize{
#' \item \code{LIN_fit=TRUE}
#' \item \code{Origin_fit=TRUE}
#' }
#' \item \bold{Saturating exponential growth and fitting through the origin} (\code{AgesMultiCS2_EXPZO}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))}; select
#' \itemize{
#' \item \code{LIN_fit=FALSE}
#' \item \code{Origin_fit=TRUE}
#' }
#' }
#'
#' **Option on equivalent dose distribution around the palaeodose**\cr
#'
#' The use can choose between :
#' \itemize{
#' \item \code{cauchy}: a Cauchy distribution with location parameter equal to the palaeodose of the sample
#' \item \code{gaussian}: a Gaussian distribution with mean equal to the palaeodose of the sample
#' \item \code{lognormal_A}: a log-normal distribution with mean or \bold{A}verage equal to the palaeodose of the sample
#' \item \code{lognormal_M}: a log-normal distribution with \bold{M}edian equal to the palaeodose of the sample
#' }
#'
#' **More precision on `Model`**\cr
#'
#' We propose two models "full" or "naive". If \code{Model='full'} that means measurement error and error on calibration curve are taken account in
#' the Bayesian model; if \code{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.
#'
#' @return
#' \bold{NUMERICAL OUTPUT}\cr
#'
#' \enumerate{
#' \item \bold{A list containing the following objects:}
#' \itemize{
#' \item \bold{Sampling}: that corresponds to a sample of the posterior distributions of the age parameters (in ka for both C14 samples and OSL samples);
#' \item \bold{PriorAge}: stating the priors used for the age parameter;
#' \item \bold{StratiConstraints}: stating the stratigraphic relations between samples considered in the model;
#' \item \bold{Model_OSL_GrowthCurve}: stating which dose response fitting option was chosen;
#' \item \bold{Model_OSL_Distribution}: stating which distribution was chosen to model the dispersion of
#' individual equivalent dose values around the palaeodose of the sample;
#' \item \bold{Model_C14}: stating which model was chosen (\code{"full"} or \code{"naive"});
#' \item \bold{CalibrationCurve}: stating which radiocarbon calibration curve was chosen;
#' \item \bold{Outlier}: stating the names of samples that must be outliers.
#' }
#'
#' \item \bold{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.\cr
#' In addition, the user must visually assess the convergence of the trajectories by looking at the graph
#' generated by the function (see \bold{PLOT OUTPUT} for more informations).\cr
#' 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 (\code{Iter})
#' or be more precise on the \code{PriorAge} parameter to reach convergence.
#' \item \bold{Credible intervals and Bayes estimates}: prints the Bayes' estimates, the credible intervals at 95% and 68% for
#' the age parameters for each sample.
#' \item \bold{JAGS model output}: returns the JAGS model with class "runjags".
#' }
#'
#' \bold{PLOT OUTPUT}
#'
#' \enumerate{
#' \item\bold{MCMC trajectories}: A graph with the MCMC trajectories and posterior distributions of the age parameter is displayed. \cr
#' On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
#' \item \bold{Age estimate and HPD at 95% of 14C samples on calibration curve}: plot age estimate and HPD on calibration plot.
#' \item \bold{Summary of sample age estimates}: plot credible intervals and Bayes estimate of each sample age on a same graph.
#' }
#'
#' @author Claire Christophe, Anne Philippe, Guillaume Guerin, Sebastian Kreutzer, Frederik Harly Baumgarten
#'
#' @note Please note that the initial values for all MCMC are currently all the same for all chains since we rely on the automatic
#' initial value generation of JAGS. This is not optimal and will be changed in future. However, it does not affect the quality
#' of the age estimates if the chains have converged.
#'
#' @seealso
#' [runjags], [plot_MCMC], [SCMatrix], [plot_Ages]
#'
#' @references
#' Reimer PJ, Bard E, Bayliss A, Beck JW, Blackwell PC, Bronl Ramsey C, Buck CE, Cheng H, Edwards RL, Friedrich M,
#' Grootes PM, Guilderson TP, Haflidason H, Hajdas I, Hatte C, Heaton TJ, Hoffmann DL, Hogg AG, Hughen KA, Kaiser KF, Kromer B,
#' Manning SW, Niu M, Reimer RW, Richards DA, Scott EM, Southon JR, Staff RA, Turney CSM, van der Plicht J. 2013.
#' IntCal13 ans Marine13 radiocarbon age calibration curves 0-50000 years cal BP. Radiocarbon 55(4)=1869-1887.
#'
#' Hogg AG, Hua Q, Blackwell PG, Niu M, Buck CE, Guilderson TP, Heaton TJ, Palmer JG, Reimer PJ, Reimer RW, Turney CSM, Zimmerman SRH.
#' 2013. SHCal13 Southern Hemisphere calibration, 0-50000 years cal BP. Radiocarbon 55(4):1889-1903
#'
#' @examples
#' ## Load data
#' # OSL data
#' data(DATA1,envir = environment())
#' data(DATA2,envir = environment())
#' Data <- combine_DataFiles(DATA2,DATA1)
#'
#' # 14C data
# data(DATA_C14,envir = environment())
#' 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
#' \dontrun{
#' Age <- Age_OSLC14(
#' DATA = Data,
#' Data_C14Cal = C14Cal,
#' Data_SigmaC14Cal = SigmaC14Cal,
#' SampleNames = c("GDB5",Names,"GDB3"),
#' Nb_sample = 3,
#' SampleNature = samplenature,
#' PriorAge = prior,
#' StratiConstraints = SC,
#' Iter = 20,
#' burnin = 20,
#' adapt = 20,
#' n.chains = 2)
#' }
#'
#' @md
#' @export
Age_OSLC14 <- function(
DATA,
Data_C14Cal,
Data_SigmaC14Cal,
Nb_sample = DATA$Nb_sample,
SampleNames = DATA$SampleNames,
SampleNature,
PriorAge = rep(c(10, 60), Nb_sample),
SavePdf = FALSE,
OutputFileName = c('MCMCplot', 'HPD_Cal14CCurve', "summary"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(''),
StratiConstraints = c(),
sepSC = c(','),
BinPerSample = rep(1, sum(SampleNature[1,])),
THETA = c(),
sepTHETA = c(','),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
Model_C14 = c("full"),
CalibrationCurve = c("IntCal20"),
Iter = 10000,
burnin = 4000,
adapt = 1000,
t = 5,
n.chains = 3,
jags_method = "rjags",
autorun = FALSE,
quiet = FALSE,
roundingOfValue = 3,
...
) {
if (inherits(DATA, "BayLum.list")) {
## reattach mcmc-list to runjags_object
DATA$runjags_object$mcmc <- DATA$Sampling
ind_OSL <- which(DATA$runjags_object$args$SampleNature[1,] == 1)
CS_OSL <- cumsum(DATA$runjags_object$args$SampleNature[1,])
ind_C14 <- which(DATA$runjags_object$args$SampleNature[2,] == 1)
CS_C14 <- cumsum(DATA$runjags_object$args$SampleNature[2,])
AgeBP = rev(DATA$runjags_object$args$TableauCalib[, 1])
CalC14 = rev(DATA$runjags_object$args$TableauCalib[, 2])
SigmaCalC14 = rev(DATA$runjags_object$args$TableauCalib[, 3])
results_runjags <-
runjags::extend.JAGS(
runjags.object = DATA$runjags_object,
adapt = adapt,
burnin = burnin,
sample = Iter,
thin = t,
method = jags_method,
silent.jags = quiet,
combine = FALSE,
...
)
# storing the arguments used for the orignal BayLum run (as to not lose them when results are processed)
results_runjags$args <- list(
"Model_OSL_GrowthCurve" = DATA$runjags_object$args$Model_OSL_GrowthCurve,
"Model_OSL_Distribution" = DATA$runjags_object$args$Model_OSL_Distribution,
"PriorAge" = DATA$runjags_object$args$PriorAge,
"StratiConstraints" = DATA$runjags_object$args$StratiConstraints,
"CovarianceMatrix" = DATA$runjags_object$args$CovarianceMatrix,
"Model_C14" = DATA$runjags_object$args$Model_C14,
"TableauCalib" = DATA$runjags_object$args$TableauCalib,
"Outlier" = DATA$runjags_object$args$Outlier,
"SampleNature" = DATA$runjags_object$args$SampleNature,
"Data_C14Cal" = DATA$runjags_object$args$Data_C14Cal,
"Nb_sample" = DATA$runjags_object$args$Nb_sample
)
}
if(!inherits(DATA, "BayLum.list")) {
#--- StratiConstraints matrix ####
if (length(StratiConstraints) == 0) {
StratiConstraints = matrix(
data = c(rep(1, Nb_sample), rep(0, Nb_sample * Nb_sample)),
ncol = Nb_sample,
nrow = (Nb_sample + 1),
byrow = T
)
} else{
if (is(StratiConstraints)[1] == "character") {
SCMatrix = read.csv(StratiConstraints, sep = sepSC)
StratiConstraints = as.matrix(SCMatrix)
}
}
#--- Calibration curve ####
TableauCalib = c()
if (CalibrationCurve %in% c(
"IntCal13", "IntCal20", "Marine13", "Marine20", "SHCal13" , "SHCal20")) {
TableauCalib = get(data(list = CalibrationCurve, envir = environment()))
} else {
TableauCalib = read.csv(file = CalibrationCurve,
sep = ",",
dec = ".")
}
AgeBP = rev(TableauCalib[, 1])
CalC14 = rev(TableauCalib[, 2])
SigmaCalC14 = rev(TableauCalib[, 3])
# #--- C14 preparation: Calibration curve
# TableauCalib=read.csv(file=paste("inst/extdata/",CalibrationCurve,"_CalC14.csv",sep=""),sep=",",dec=".")
# AgeBP=rev(TableauCalib[,1])/1000
# CalC14=rev(TableauCalib[,2])
# SigmaCalC14=rev(TableauCalib[,3])
#--- OSL preparation
#- Index preparation ####
CSBinPerSample = cumsum(BinPerSample)
LengthSample = c()
for (ns in 1:sum(SampleNature[1, ])) {
LengthSample = c(LengthSample, length(DATA$LT[[ns]][, 1]))
}
CSLengthSample = c()
CSLengthSample = c(0, cumsum(LengthSample))
index2 = c(0, cumsum(DATA$J))
#- File preparation ####
LT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K) + 1))
sLT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K) + 1))
IrrT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K)))
for (ns in 1:sum(SampleNature[1, ])) {
LT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$LT[[ns]][1, ])] <-
DATA$LT[[ns]]
sLT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$sLT[[ns]][1, ])] <-
DATA$sLT[[ns]]
IrrT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$ITimes[[ns]][1, ])] <-
DATA$ITimes[[ns]]
}
#- THETA matrix ####
if (length(THETA[, 1]) == 0) {
THETA = diag(DATA$ddot_env[2, CSBinPerSample] + (DATA$ddot_env[1, CSBinPerSample]) ^
2 * DATA$dLab[2, CSBinPerSample])
} else{
if (is(THETA)[1] == "character") {
errorMatrix = read.csv(THETA, sep = sepTHETA)
THETA = as.matrix(errorMatrix)
}
}
#--- Index preparation ####
ind_OSL <- which(SampleNature[1,] == 1)
CS_OSL <- cumsum(SampleNature[1,])
ind_C14 <- which(SampleNature[2,] == 1)
CS_C14 <- cumsum(SampleNature[2,])
ind_change <- c(1)
for (i in 2:(Nb_sample - 1)) {
if (SampleNature[1, i] != SampleNature[1, i + 1]) {
ind_change <- c(ind_change, i)
}
}
ind_change <- c(ind_change, Nb_sample)
q <- length(ind_change) %/% 2
r <- length(ind_change) %% 2
##--- description du model BUG ####
BUGModel <- c()
#- Prior
ModelPrior <- 0
data(ModelPrior, envir = environment())
BUGPrior <- c()
if (r == 1) {
if (SampleNature[1, 1] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_OSL)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_C14)
}
if (SampleNature[1, 2] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$OSL_C14)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$C14_OSL)
}
} else{
q <- q - 1
if (q == 0) {
stop(
"[Age_OSLC14()] If you see this message, you are probably trying to run the model with a small number of samples.
You can still use the function, but the C-14 sample cannot be the first sample.",
call. = FALSE
)
}
if (SampleNature[1, 1] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_OSL)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_C14)
}
if (SampleNature[1, 2] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$OSL_C14)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$C14_OSL)
}
if (SampleNature[1, Nb_sample] == 1) {
BUGPrior = paste(BUGPrior, ModelPrior$OSL)
} else{
BUGPrior = paste(BUGPrior, ModelPrior$C14)
}
}
#- partie C14
ModelC14 <- 0
data(ModelC14, envir = environment())
if (Model_C14 == "full") {
BUGModel = paste(ModelC14$full, BUGPrior)
} else{
BUGModel = paste(ModelC14$naive, BUGPrior)
}
#- partie OSL
ModelOSL <- 0
data(ModelOSL, envir = environment())
if (LIN_fit == TRUE) {
cLIN = c('LIN')
} else{
cLIN = c()
}
if (Origin_fit == TRUE) {
cO = c("ZO")
} else{
cO = c()
}
Model_GrowthCurve = c(paste("AgesMultiOSL_EXP", cLIN, cO, sep = ""))
BUGModel = c(paste("model{", ModelOSL[[Model_GrowthCurve]][[distribution]], BUGModel, "}"))
if (Model_C14 == "full") {
dataList = list(
'q' = q,
"ind_change" = ind_change,
"ind_OSL" = ind_OSL,
"ind_C14" = ind_C14,
"CS_OSL" = CS_OSL,
"CS_C14" = CS_C14,
'X' = Data_C14Cal,
"sigma" = Data_SigmaC14Cal,
"xTableauCalib" = AgeBP,
"yTableauCalib" = CalC14,
"zTableauCalib" = SigmaCalC14,
'N' = LT,
'sN' = sLT,
"IT" = IrrT,
"sDlab" = DATA$dLab[1, ],
'J' = DATA$J,
'K' = DATA$K,
"ddot" = DATA$ddot_env[1, CSBinPerSample],
"Gamma" = THETA,
"index" = index2,
"BinPerSample" = BinPerSample,
"CSBinPerSample" = CSBinPerSample,
"xbound" = PriorAge,
"StratiConstraints" = StratiConstraints
)
} else{
dataList = list(
'q' = q,
"ind_change" = ind_change,
"ind_OSL" = ind_OSL,
"ind_C14" = ind_C14,
"CS_OSL" = CS_OSL,
"CS_C14" = CS_C14,
'X' = Data_C14Cal,
"sigma" = Data_SigmaC14Cal,
"xTableauCalib" = AgeBP,
"yTableauCalib" = CalC14,
'N' = LT,
'sN' = sLT,
"IT" = IrrT,
"sDlab" = DATA$dLab[1, ],
'J' = DATA$J,
'K' = DATA$K,
"ddot" = DATA$ddot_env[1, CSBinPerSample],
"Gamma" = THETA,
"index" = index2,
"BinPerSample" = BinPerSample,
"CSBinPerSample" = CSBinPerSample,
"xbound" = PriorAge,
"StratiConstraints" = StratiConstraints
)
}
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# JAGS RUN --------------------- START
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
##further settings provided eventually
process_settings <- modifyList(x = list(
max.time = Inf,
interactive = FALSE,
startburnin = 4000,
startsample = 10000,
inits = NA
), val = list(...))
if (autorun == FALSE) {
##as input a text file is wanted, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(BUGModel, con = temp_file)
##run JAGS
results_runjags <-
runjags::run.JAGS(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "Z"),
adapt = adapt,
burnin = burnin,
sample = Iter,
silent.jags = quiet,
method = jags_method,
thin = t
)
}
if (autorun == TRUE) {
##as input a text file is wanted, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(BUGModel, con = temp_file)
##run the auto processing
results_runjags <-
runjags::autorun.jags(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "Z"),
adapt = adapt,
silent.jags = quiet,
method = jags_method,
inits = process_settings$inits,
max.time = process_settings$max.time,
interactive = process_settings$interactive,
startburnin = process_settings$startburnin,
startsample = process_settings$startsample
)
}
results_runjags$args <- list(
"Model_OSL_GrowthCurve" = Model_GrowthCurve,
"Model_OSL_Distribution" = distribution,
"PriorAge" = PriorAge,
"StratiConstraints" = StratiConstraints,
"CovarianceMatrix" = THETA,
"Model_C14" = Model_C14,
"TableauCalib" = TableauCalib,
"SampleNature" = SampleNature,
"Data_C14Cal" = Data_C14Cal,
"Nb_sample" = Nb_sample
)
}
##extract mcmc list
echantillon <- results_runjags$mcmc
U <- summary(echantillon)
##remove mcmc-list from runjags output to reduce output object size
results_runjags$mcmc <- list("MCMC-list is not here. Go to first level -> object named *Sampling*")
##combine chains into one data.frame
Sample <- as.data.frame(runjags::combine.mcmc(echantillon))
## attach "A_" to sample names for the A-parameter
nom = c()
for (i in 1:results_runjags$args$Nb_sample) {
nom = c(nom, paste("A_", SampleNames[i], sep = ""))
}
##plot MCMC
if (SavePdf) {
pdf(file = paste(OutputFilePath, OutputFileName[1], '.pdf', sep = ""))
}
plot_MCMC(echantillon, sample_names = SampleNames)
if (SavePdf) {
dev.off()
}
Outlier <-
SampleNames[ind_C14[which(U$statistics[(results_runjags$args$Nb_sample + 1):(results_runjags$args$Nb_sample + sum(results_runjags$args$SampleNature[2, ])), 1] <
1.5)]]
##- Gelman and Rubin test of convergency of the MCMC
CV = gelman.diag(echantillon, multivariate = FALSE)
cat(paste("\n\n>> Convergence of MCMC for the age parameters <<\n"))
cat("----------------------------------------------\n")
cat(paste("Sample name ", " Bayes estimate ", " Uppers credible interval\n"))
for (i in 1:results_runjags$args$Nb_sample) {
#cat(paste(" Sample name: ", SampleNames[i],"\n"))
#cat("---------------------\n")
cat(paste(
paste("A_", SampleNames[i], sep = ""),
"\t",
round(CV$psrf[i, 1], roundingOfValue),
"\t\t",
round(CV$psrf[i, 2], roundingOfValue),
"\n"
))
}
cat(
"\n\n________________________________________________________________________________\n"
)
cat(" *** WARNING: following informations are only valid if MCMC chains converged ***\n")
cat(
"________________________________________________________________________________\n"
)
# Matrix of results
rnames = c()
for (i in 1:results_runjags$args$Nb_sample) {
rnames = c(rnames, paste("A_", SampleNames[i], sep = ""))
}
R = matrix(
data = NA,
ncol = 8,
nrow = results_runjags$args$Nb_sample,
dimnames = list(
rnames,
c(
"lower bound at 95%",
"lower bound at 68%",
"Bayes estimate",
"upper bound at 68%",
"upper bound at 95%",
"",
"Convergencies: Bayes estimate",
"Convergencies: uppers credible interval"
)
)
)
##- Bayes estimate and credible interval
cat(paste(
"\n\n>> Bayes estimates of Age for each sample and credible interval <<\n"
))
AgePlot95 = matrix(data = NA, nrow = results_runjags$args$Nb_sample, ncol = 3)
AgePlot68 = matrix(data = NA, nrow = results_runjags$args$Nb_sample, ncol = 3)
AgePlotMoy = rep(0, results_runjags$args$Nb_sample)
for (i in 1:results_runjags$args$Nb_sample) {
cat("------------------------------------------------------\n")
#cat(paste(" Sample name: ", SampleNames[i],"\n"))
#cat("---------------------\n")
cat(paste("Sample name", "\t", "Bayes estimate", " Credible interval: \n"))
cat(paste(
paste("A_", SampleNames[i], sep = ""),
"\t",
round(mean(Sample[, i]), roundingOfValue),
'\n'
))
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
HPD_95 = CredibleInterval(Sample[, i], 0.95, roundingOfValue =
roundingOfValue)
HPD_68 = CredibleInterval(Sample[, i], 0.68, roundingOfValue =
roundingOfValue)
cat(
"\t\t\t\t at level 95% \t",
round(c(HPD_95[2]), roundingOfValue),
"\t\t",
round(c(HPD_95[3]), roundingOfValue),
"\n"
)
cat(
"\t\t\t\t at level 68% \t",
round(c(HPD_68[2]), roundingOfValue),
"\t\t",
round(c(HPD_68[3]), roundingOfValue),
"\n"
)
AgePlot95[i, ] = HPD_95
AgePlot68[i, ] = HPD_68
AgePlotMoy[i] = round(mean(Sample[, i]), roundingOfValue)
R[i, 3] = round(mean(Sample[, i]), roundingOfValue)
R[i, c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
R[i, c(2, 4)] = round(HPD_68[2:3], roundingOfValue)
R[i, 6] = c('')
R[i, 7] = round(CV$psrf[i, 1], roundingOfValue)
R[i, 8] = round(CV$psrf[i, 2], roundingOfValue)
}
cat("\n------------------------------------------------------\n")
R[, c(7, 8)] <- round(CV$psrf[1:results_runjags$args$Nb_sample, ], roundingOfValue)
# Representation graphique des resultats
# des HPD sur la courbe de calibration
if (sum(results_runjags$args$SampleNature[2, ]) > 1) {
couleur = rainbow(results_runjags$args$Nb_sample)
par(mfrow = c(1, 1),
las = 0,
mar = c(5, 5, 2, 2))
xl = c(min(results_runjags$args$PriorAge[seq(1, (2 * results_runjags$args$Nb_sample - 1), 2)]), max(results_runjags$args$PriorAge[seq(2, (2 *
results_runjags$args$Nb_sample), 2)]))
plot(
xl,
xl,
col = "white",
xlab = c("Age (in ka)"),
ylab = c("cal C14"),
xaxt = "n",
yaxt = "n",
cex.lab = 1.8
)
axis(2, cex.axis = 2)
axis(1, cex.axis = 2)
polygon(c(AgeBP, rev(AgeBP)),
c(CalC14 + 2 * SigmaCalC14, rev(CalC14 - 2 * SigmaCalC14)),
col = "gray",
border = "black")
for (i in ind_C14) {
lines(c(AgePlot95[i, 2:3]),
rep(results_runjags$args$Data_C14Cal[CS_C14[i]], 2),
col = couleur[i],
lwd = 4)
lines(
AgePlotMoy[i],
results_runjags$args$Data_C14Cal[CS_C14[i]],
col = "black",
lwd = 2,
type = 'p'
)
}
legend(
"topleft",
SampleNames[ind_C14],
lty = rep(1, results_runjags$args$Nb_sample),
lwd = rep(2, results_runjags$args$Nb_sample),
cex = 1,
col = couleur[ind_C14]
)
if (SavePdf == TRUE) {
dev.print(
pdf,
file = paste(OutputFilePath, OutputFileName[2], '.pdf', sep = ""),
width = 8,
height = 10
)
}
}
##CSV output
if (SaveEstimates == TRUE) {
write.csv(R, file = c(
paste(OutputTablePath, "Estimates", OutputTableName, ".csv", sep = "")
))
}
.list_BayLum <- function(..., originator = sys.call(which = -1)[[1]]){
## set list
l <- list(...)
## update originators
attr(l, "class") <- "BayLum.list"
attr(l, "originator") <- as.character(originator)
return(l)
}
# Create return object -------------------------------------------------------------------------
output <- .list_BayLum(
"Ages" = data.frame(
SAMPLE = SampleNames,
AGE = AgePlotMoy,
HPD68.MIN = AgePlot68[, 2],
HPD68.MAX = AgePlot68[, 3],
HPD95.MIN = AgePlot95[, 2],
HPD95.MAX = AgePlot95[, 3],
stringsAsFactors = FALSE
),
"Sampling" = echantillon,
"PriorAge" = results_runjags$args$PriorAge,
"StratiConstraints" = results_runjags$args$StratiConstraints,
"Model_OSL_GrowthCurve" = results_runjags$args$Model_OSL_GrowthCurve,
"Model_OSL_Distribution" = results_runjags$args$Model_OSL_Distribution,
"CovarianceMatrix" = results_runjags$args$CovarianceMatrix,
"Model_C14" = results_runjags$args$Model_C14,
"CalibrationCurve" = results_runjags$args$CalibrationCurve,
"Outlier" = Outlier,
"runjags_object" = results_runjags
)
# Plot ages -----------------------------------------------------------------------------------
plot_Ages(object = output, legend.pos = "bottomleft")
##TODO: get rid of this
if (SavePdf) {
dev.print(
pdf,
file = paste(OutputFilePath, OutputFileName[3], '.pdf', sep = ""),
width = 8,
height = 10
)
}
# Return output -------------------------------------------------------------------------------
return(output)
}
R-Lum/BayLum documentation built on April 19, 2024, 9:33 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Age_OSLC14
Documented in Age_OSLC14
#' @title Bayesian analysis for age estimation of OSL measurements and C-14 ages of various samples
#'
#' @description This function computes an age of OSL data consisting of at least two samples and calibrate
#' C-14 ages of samples to get an age (in ka).\cr
#' Ages of OSL data are computed according to the model given in Combes and Philippe (2017).
#' Single-grain or Multi-grain OSL measurements can be analysed simultaneously
#' (with output of [Generate_DataFile] or [Generate_DataFile_MG] or both of them using [combine_DataFiles]).
#' Samples, for which data is available in several BIN files, can be analysed.\cr
#' For C-14 data, the user can choose one of the following radiocarbon calibration curve:
#' Northern or Southern Hemisphere or marine atmospheric.
#'
#' @param DATA Two types of inputs are supported.
#' (1): a list of objects: `LT`, `sLT`, `ITimes`, `dLab`, `ddot_env`, `regDose`, `J`,
#' `K`, `Nb_measurement`, provided by the function [Generate_DataFile], [Generate_DataFile_MG] or [combine_DataFiles].
#' \code{DATA} contains information for more than one sample.
#' If there is stratigraphic relations between samples, informations in DATA must be ordered by order of
#' increasing ages. See the details section to for more informations.
#' (2): an object of class `BayLum.list` which is provided by the output of [Age_OSLC14]. When input of class `BayLum.list` is identified, no new JAGS model is created. Instead, the JAGS model specified within the `BayLum.list` is extended. Useful for when convergence was not originally achieved and a complete restart is not desirable.
#'
#' @param Data_C14Cal [numeric] vector: corresponding to C-14 age estimate
#' (in years, conversion in ka is automatically done in the function).
#' If there is stratigraphic relations between samples, `Data_C14Cal` must be ordered
#' by order of increasing ages.
#'
#' @param Data_SigmaC14Cal [numeric]: corresponding to the error of C-14 age estimates.
#'
#' @param Nb_sample [numeric]: number of samples (OSL data and C-14 age),
#' (`Nb_sample>3`, at least to sample of OSL data and one sample of C-14 age).
#'
#' @param SampleNames [character]: sample names for both OSL data and C14 data.
#' The length of this vector is equal to `Nb_sample`.
#' If there is stratigraphic relation, this vector must be ordered by increasing
#' order (to mix OSL samples and C-14 ages if it is needed).
#'
#' @param SampleNature [matrix]: states the nature of the sample.
#' Row number of `SampleNature` matrix is equal to `2` and column number is equal
#' to `Nb_sample`. First line of the matrix: `SampleNature[1,i]` states if sample
#' whose number ID is equal to `i`, is an OSL sample `1` or not `0`.
#' Second line of the matrix: `SampleNature[2,i]` states if sample whose number ID
#' is equal to `i`, is an C-14 sample `1` or not `0.
#'
#' @param PriorAge [numeric] (with default): lower and upper bounds for age parameter of each sample **(in ka)**.
#' Note that, `length(PriorAge) = 2*Nb_sample` sand `PriorAge[2i-1,2i]` corresponds
#' to the lower and upper bounds of sample whose number ID is equal to `i`.
#'
#' @param SavePdf [logical] (with default): if `TRUE` save graphs in pdf file
#' named `OutputFileName` in folder `OutputFilePath`.
#'
#' @param OutputFileName [character] (with default): name of the pdf file that
#' will be generated by the function if `SavePdf=TRUE`, `length(OutputFileName)=3`,
#' see **PLOT OUTPUT** in **Value** section for more informations.
#'
#' @param 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 "/".
#'
#' @param 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`.
#'
#' @param OutputTableName [character] (with default): name of the table that will
#' be generated by the function if `SaveEstimates=TRUE`.
#'
#' @param 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 "/".
#'
#' @param StratiConstraints [matrix] or [character] (with default): input object for the
#' stratigraphic relation between samples. If there is stratigraphic relation between
#' samples see the details section for instructions regarding how to correctly fill `StratiConstraints`,
#' the user can refer to a [matrix] or to a CSV-file [character]. Otherwise, default value is suitable.
#'
#' @param sepSC [character] (with default): if `StratiConstraints` is character,
#' indicate column separator in `StratiConstraints` CSV-file.
#'
#' @param BinPerSample [numeric] (with default): vector with the number of BIN-files
#' per OSL sample. The length of this vector is equal to the number of OSL samples.
#' `BinPerSample[i]` corresponds to the number of BIN files for the sample whose
#' number ID is equal to `i`. For more information to fill this vector, we refer
#' to details in [Generate_DataFile] or in [Generate_DataFile_MG].
#'
#' @param THETA numeric [matrix] or [character] (with default): input object for
#' systematic and individual error for OSL samples. If systematic errors are considered,
#' see the details section for instructions regarding how to correctly fill `THETA`;
#' the user can refer to a matrix (numeric matrix) or to a csv file (character).
#' Otherwise, default value is suitable, and only individual error is considered.
#'
#' @param sepTHETA [character] (with default): if `THETA` is character,
#' indicate column separator in `THETA` CSV-file.
#'
#' @param LIN_fit [logical] (with default): if `TRUE` (default) allows a linear component,
#' on top of the (default) saturating exponential curve, for the fitting of dose response curves,
#' for OSL samples. See details for more informations on the proposed dose response curves.
#'
#' @param Origin_fit p[logical] (with default): if `TRUE`, forces the dose response
#' curves to pass through the origin. See details for more informations on the proposed
#' growth curves, for OSL samples.
#'
#' @param distribution [character] (with default): type of distribution that defines
#' how individual equivalent dose values are distributed around the palaeodose, for OSL samples.
#' Allowed inputs are `"cauchy"`, `"gaussian"`, `"lognormal_A"` and `"lognormal_M"`,
#' see details for more informations.
#'
#' @param Model_C14 [character] (with default): if `"full"`, error on estimate calibration
#' curve is taken account, for C-14 samples. If `"naive"` this error is not taken account in the age estimate.
#'
#' @param CalibrationCurve [character] (with default): calibration curve chosen, for C-14 samples.
#' Allowed inputs are
#' \itemize{
#' \item \bold{"Intcal13"} or \bold{"Intcal13"} for Northern Hemisphere atmospheric radiocarbon calibration curve,
#' \item \bold{"Marine13"} or \bold{"Marine13"} for Marine radiocarbon calibration curve,
#' \item \bold{"SHCal13"} or \bold{"SHCal20"} for Southern Hemisphere atmospheric radiocarbon calibration curve
#' \item \bold{a csv file, with tree columns, the first column is dedicated to "Cal.BP", the second to "XC-14.age", the third to "Error".
#' The decimal of this file must be a dot, and the separator must be a comma. }
#' }
#'
#' @param adapt [integer] (with default): the number of iterations used in the adaptive phase of the simulation (see [runjags::run.jags]]).
#' @param 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]]).
#' @param Iter [integer] (with default): the number of iterations to run and who will be used to assess convergence and ages (see [runjags::run.jags]]).
#'
#' @param t [numeric] (with default): 1 every `t` iterations of the MCMC is
#' considered for sampling the posterior distribution (for more information see [[rjags::jags.model].
#'
#' @param n.chains [numeric] (with default): number of independent chains for the model
#' (for more information see [[rjags::jags.model]).
#'
#' @param jags_method [character] (with default): select which method to use in order to call JAGS, supported are `"rjags"` (the default), `rjparallel`, `simple`, `interruptible`, `parallel`, and `snow` (for more information about each of these possibilities, see [runjags::run.jags]])
#'
#' @param 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]).
#'
#' @param quiet [logical] (with default): enables/disables [rjags] messages
#'
#' @param roundingOfValue [integer] (with default): Integer indicating the number of decimal places to be used, default = 3.
#
#' @param ... further arguments that can be passed to control the Bayesian process, see details
#' for supported arguments
#'
#' @details
#'
#' Note that there are three types of arguments in the previous list.
#' There are arguments for information concerning only OSL samples: \code{DATA}, \code{BinPerSample}, \code{THETA},
#' \code{sepTHETA}, \code{LIN_fit}, \code{Origin_fit}, \code{distribution}.
#'
#' There are arguments for information concerning only C14 samples: \code{Data_C14Cal}, \code{Data_SigmaC14Cal},
#' \code{Model_C14}, \code{CalibrationCurve}.
#'
#' There are arguments for information concerning all the samples: \code{Nb_sample}, \code{SampleNames}, \code{SampleNature},
#' \code{PriorAge}, \code{SavePdf}, \code{OutputFileName}, \code{OutputFilePath}, \code{SaveEstimates}, \code{OutputTableName},
#' \code{OutputTablePath}, \code{StratiConstraints}, \code{sepSC}.\cr
#'
#' **Supported `...` arguments**
#' \tabular{lllll}{
#' ARGUMENT \tab INPUT \tab METHOD \tab DEFAULT \tab DESCRIPTION\cr
#' `max.time` \tab [character] \tab `rjparallel` \tab `Inf` \tab maximum allowed processing time, e.g.,
#' `10m` for 10 minutes (cf. [runjags::autorun.jags])\cr
#' `interactive` \tab [logical] \tab `rjparallel` \tab `FALSE` \tab enable/disable interactive mode (cf. [runjags::autorun.jags])\cr
#' `startburnin` \tab [integer] \tab `rjparallel` \tab `4000` \tab number of burn-in iterations (cf. [runjags::autorun.jags]) \cr
#' `startsample` \tab [integer] \tab `rjparallel` \tab `10000` \tab total number of samples to assess convergence
#' (cf. [runjags::autorun.jags]) \cr
#' `inits` \tab named [list] \tab `rjparallel` \tab `NA` \tab fine control over random seeds and random number generator [runjags::autorun.jags]
#' }
#'
#' **How to fill `StratiConstraints?**\cr
#'
#' If there are stratigraphic relations between samples, \bold{14C estimate age in \code{Data_C14Cal} must be ordered by order of increasing ages,
#' as informations in \code{DATA}}. Names in \code{SampleNames} must be ordered and corresponds to the order in \code{Data_C14Cal} and in \code{DATA},
#' also if it is needed to mix names of OSL samples and 14C samples.
#'
#' The user can fill the \code{StratiConstraints} matrix as follow.
#' \enumerate{
#' \item \bold{Size of the matrix}: row number of \code{StratiConstraints} matrix is equal to \code{Nb_sample+1},
#' and column number is equal to \code{Nb_sample}.
#' \item \bold{First line of the matrix}:
#' for all \code{i in {1,...,Nb_Sample}}, \code{StratiConstraints[1,i]=1} that means the lower bound of the sample age (given in \code{PriorAge[2i-1]})
#' for the sample whose number ID is equal to \code{i}, is taken into account.
#' \item \bold{Sample relations}: for all \code{j in {2,...,Nb_Sample+1}} and all \code{i in {j,...,Nb_Sample}},
#' \code{StratiConstraints[j,i]=1} if sample age whose number ID is equal to \code{j-1} is lower than sample age whose number ID is equal to \code{i}.
#' Otherwise, \code{StratiConstraints[j,i]=0}.
#' }
#' Note that \code{StratiConstraints_{2:Nb_sample+1,1:Nb_sample}} is a upper triangular matrix.
#'
#' The user can also use \code{\link{SCMatrix}} or \code{\link{SC_Ordered}} (if all samples are ordered) function to construct
#' the \code{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 \code{sepSC}.\cr
#'
#' **How to fill `THETA` covariance matrix concerning common and individual error?**\cr
#'
#' If systematic errors are considered, the user can fill the \code{THETA} matrix as follow.
#' \itemize{
#' \item row number of \code{THETA} is equal the column number, equal to \code{Nb_sample}.
#' \item For all \code{i in {1,...,Nb_sample}}, \code{THETA[i,i]} contains individual error
#' plus systematic error of the sample whose number ID is equal to \code{i}.
#' \item For all \code{i,j in {1,...,Nb_sample}} and \code{i} different from \code{j} ,
#' \code{THETA[i,j]} contains common error between samples whose number ID are equal to \code{i} and \code{j}.
#' }
#' Note that \code{THETA[i,j]} is a symmetric matrix.
#'
#' The user can also refer to a .csv file that contains the errors as defined above.\cr
#'
#' **Option on growth curves**\cr
#'
#' As for \code{\link{Age_Computation}} and \code{\link{Palaeodose_Computation}}, the user can choose from 4 dose response curves:
#' \itemize{
#' \item \bold{Saturating exponential plus linear growth} (\code{AgesMultiCS2_EXPLIN}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx+d}; select
#' \itemize{
#' \item \code{LIN_fit=TRUE}
#' \item \code{Origin_fit=FALSE}
#' }
#' \item \bold{Saturating exponential growth} (\code{AgesMultiCS2_EXP}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+d}; select
#' \itemize{
#' \item \code{LIN_fit=FALSE}
#' \item \code{Origin_fit=FALSE}
#' }
#' \item \bold{Saturating exponential plus linear growth and fitting through the origin} (\code{AgesMultiCS2_EXPLINZO}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))+cx}; select
#' \itemize{
#' \item \code{LIN_fit=TRUE}
#' \item \code{Origin_fit=TRUE}
#' }
#' \item \bold{Saturating exponential growth and fitting through the origin} (\code{AgesMultiCS2_EXPZO}):
#'
#' for all \code{x} in IR+, \code{f(x)=a(1-exp(-x/b))}; select
#' \itemize{
#' \item \code{LIN_fit=FALSE}
#' \item \code{Origin_fit=TRUE}
#' }
#' }
#'
#' **Option on equivalent dose distribution around the palaeodose**\cr
#'
#' The use can choose between :
#' \itemize{
#' \item \code{cauchy}: a Cauchy distribution with location parameter equal to the palaeodose of the sample
#' \item \code{gaussian}: a Gaussian distribution with mean equal to the palaeodose of the sample
#' \item \code{lognormal_A}: a log-normal distribution with mean or \bold{A}verage equal to the palaeodose of the sample
#' \item \code{lognormal_M}: a log-normal distribution with \bold{M}edian equal to the palaeodose of the sample
#' }
#'
#' **More precision on `Model`**\cr
#'
#' We propose two models "full" or "naive". If \code{Model='full'} that means measurement error and error on calibration curve are taken account in
#' the Bayesian model; if \code{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.
#'
#' @return
#' \bold{NUMERICAL OUTPUT}\cr
#'
#' \enumerate{
#' \item \bold{A list containing the following objects:}
#' \itemize{
#' \item \bold{Sampling}: that corresponds to a sample of the posterior distributions of the age parameters (in ka for both C14 samples and OSL samples);
#' \item \bold{PriorAge}: stating the priors used for the age parameter;
#' \item \bold{StratiConstraints}: stating the stratigraphic relations between samples considered in the model;
#' \item \bold{Model_OSL_GrowthCurve}: stating which dose response fitting option was chosen;
#' \item \bold{Model_OSL_Distribution}: stating which distribution was chosen to model the dispersion of
#' individual equivalent dose values around the palaeodose of the sample;
#' \item \bold{Model_C14}: stating which model was chosen (\code{"full"} or \code{"naive"});
#' \item \bold{CalibrationCurve}: stating which radiocarbon calibration curve was chosen;
#' \item \bold{Outlier}: stating the names of samples that must be outliers.
#' }
#'
#' \item \bold{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.\cr
#' In addition, the user must visually assess the convergence of the trajectories by looking at the graph
#' generated by the function (see \bold{PLOT OUTPUT} for more informations).\cr
#' 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 (\code{Iter})
#' or be more precise on the \code{PriorAge} parameter to reach convergence.
#' \item \bold{Credible intervals and Bayes estimates}: prints the Bayes' estimates, the credible intervals at 95% and 68% for
#' the age parameters for each sample.
#' \item \bold{JAGS model output}: returns the JAGS model with class "runjags".
#' }
#'
#' \bold{PLOT OUTPUT}
#'
#' \enumerate{
#' \item\bold{MCMC trajectories}: A graph with the MCMC trajectories and posterior distributions of the age parameter is displayed. \cr
#' On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
#' \item \bold{Age estimate and HPD at 95% of 14C samples on calibration curve}: plot age estimate and HPD on calibration plot.
#' \item \bold{Summary of sample age estimates}: plot credible intervals and Bayes estimate of each sample age on a same graph.
#' }
#'
#' @author Claire Christophe, Anne Philippe, Guillaume Guerin, Sebastian Kreutzer, Frederik Harly Baumgarten
#'
#' @note Please note that the initial values for all MCMC are currently all the same for all chains since we rely on the automatic
#' initial value generation of JAGS. This is not optimal and will be changed in future. However, it does not affect the quality
#' of the age estimates if the chains have converged.
#'
#' @seealso
#' [runjags], [plot_MCMC], [SCMatrix], [plot_Ages]
#'
#' @references
#' Reimer PJ, Bard E, Bayliss A, Beck JW, Blackwell PC, Bronl Ramsey C, Buck CE, Cheng H, Edwards RL, Friedrich M,
#' Grootes PM, Guilderson TP, Haflidason H, Hajdas I, Hatte C, Heaton TJ, Hoffmann DL, Hogg AG, Hughen KA, Kaiser KF, Kromer B,
#' Manning SW, Niu M, Reimer RW, Richards DA, Scott EM, Southon JR, Staff RA, Turney CSM, van der Plicht J. 2013.
#' IntCal13 ans Marine13 radiocarbon age calibration curves 0-50000 years cal BP. Radiocarbon 55(4)=1869-1887.
#'
#' Hogg AG, Hua Q, Blackwell PG, Niu M, Buck CE, Guilderson TP, Heaton TJ, Palmer JG, Reimer PJ, Reimer RW, Turney CSM, Zimmerman SRH.
#' 2013. SHCal13 Southern Hemisphere calibration, 0-50000 years cal BP. Radiocarbon 55(4):1889-1903
#'
#' @examples
#' ## Load data
#' # OSL data
#' data(DATA1,envir = environment())
#' data(DATA2,envir = environment())
#' Data <- combine_DataFiles(DATA2,DATA1)
#'
#' # 14C data
# data(DATA_C14,envir = environment())
#' 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
#' \dontrun{
#' Age <- Age_OSLC14(
#' DATA = Data,
#' Data_C14Cal = C14Cal,
#' Data_SigmaC14Cal = SigmaC14Cal,
#' SampleNames = c("GDB5",Names,"GDB3"),
#' Nb_sample = 3,
#' SampleNature = samplenature,
#' PriorAge = prior,
#' StratiConstraints = SC,
#' Iter = 20,
#' burnin = 20,
#' adapt = 20,
#' n.chains = 2)
#' }
#'
#' @md
#' @export
Age_OSLC14 <- function(
DATA,
Data_C14Cal,
Data_SigmaC14Cal,
Nb_sample = DATA$Nb_sample,
SampleNames = DATA$SampleNames,
SampleNature,
PriorAge = rep(c(10, 60), Nb_sample),
SavePdf = FALSE,
OutputFileName = c('MCMCplot', 'HPD_Cal14CCurve', "summary"),
OutputFilePath = c(""),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(''),
StratiConstraints = c(),
sepSC = c(','),
BinPerSample = rep(1, sum(SampleNature[1,])),
THETA = c(),
sepTHETA = c(','),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
Model_C14 = c("full"),
CalibrationCurve = c("IntCal20"),
Iter = 10000,
burnin = 4000,
adapt = 1000,
t = 5,
n.chains = 3,
jags_method = "rjags",
autorun = FALSE,
quiet = FALSE,
roundingOfValue = 3,
...
) {
if (inherits(DATA, "BayLum.list")) {
## reattach mcmc-list to runjags_object
DATA$runjags_object$mcmc <- DATA$Sampling
ind_OSL <- which(DATA$runjags_object$args$SampleNature[1,] == 1)
CS_OSL <- cumsum(DATA$runjags_object$args$SampleNature[1,])
ind_C14 <- which(DATA$runjags_object$args$SampleNature[2,] == 1)
CS_C14 <- cumsum(DATA$runjags_object$args$SampleNature[2,])
AgeBP = rev(DATA$runjags_object$args$TableauCalib[, 1])
CalC14 = rev(DATA$runjags_object$args$TableauCalib[, 2])
SigmaCalC14 = rev(DATA$runjags_object$args$TableauCalib[, 3])
results_runjags <-
runjags::extend.JAGS(
runjags.object = DATA$runjags_object,
adapt = adapt,
burnin = burnin,
sample = Iter,
thin = t,
method = jags_method,
silent.jags = quiet,
combine = FALSE,
...
)
# storing the arguments used for the orignal BayLum run (as to not lose them when results are processed)
results_runjags$args <- list(
"Model_OSL_GrowthCurve" = DATA$runjags_object$args$Model_OSL_GrowthCurve,
"Model_OSL_Distribution" = DATA$runjags_object$args$Model_OSL_Distribution,
"PriorAge" = DATA$runjags_object$args$PriorAge,
"StratiConstraints" = DATA$runjags_object$args$StratiConstraints,
"CovarianceMatrix" = DATA$runjags_object$args$CovarianceMatrix,
"Model_C14" = DATA$runjags_object$args$Model_C14,
"TableauCalib" = DATA$runjags_object$args$TableauCalib,
"Outlier" = DATA$runjags_object$args$Outlier,
"SampleNature" = DATA$runjags_object$args$SampleNature,
"Data_C14Cal" = DATA$runjags_object$args$Data_C14Cal,
"Nb_sample" = DATA$runjags_object$args$Nb_sample
)
}
if(!inherits(DATA, "BayLum.list")) {
#--- StratiConstraints matrix ####
if (length(StratiConstraints) == 0) {
StratiConstraints = matrix(
data = c(rep(1, Nb_sample), rep(0, Nb_sample * Nb_sample)),
ncol = Nb_sample,
nrow = (Nb_sample + 1),
byrow = T
)
} else{
if (is(StratiConstraints)[1] == "character") {
SCMatrix = read.csv(StratiConstraints, sep = sepSC)
StratiConstraints = as.matrix(SCMatrix)
}
}
#--- Calibration curve ####
TableauCalib = c()
if (CalibrationCurve %in% c(
"IntCal13", "IntCal20", "Marine13", "Marine20", "SHCal13" , "SHCal20")) {
TableauCalib = get(data(list = CalibrationCurve, envir = environment()))
} else {
TableauCalib = read.csv(file = CalibrationCurve,
sep = ",",
dec = ".")
}
AgeBP = rev(TableauCalib[, 1])
CalC14 = rev(TableauCalib[, 2])
SigmaCalC14 = rev(TableauCalib[, 3])
# #--- C14 preparation: Calibration curve
# TableauCalib=read.csv(file=paste("inst/extdata/",CalibrationCurve,"_CalC14.csv",sep=""),sep=",",dec=".")
# AgeBP=rev(TableauCalib[,1])/1000
# CalC14=rev(TableauCalib[,2])
# SigmaCalC14=rev(TableauCalib[,3])
#--- OSL preparation
#- Index preparation ####
CSBinPerSample = cumsum(BinPerSample)
LengthSample = c()
for (ns in 1:sum(SampleNature[1, ])) {
LengthSample = c(LengthSample, length(DATA$LT[[ns]][, 1]))
}
CSLengthSample = c()
CSLengthSample = c(0, cumsum(LengthSample))
index2 = c(0, cumsum(DATA$J))
#- File preparation ####
LT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K) + 1))
sLT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K) + 1))
IrrT = matrix(data = 0,
nrow = sum(DATA$J),
ncol = (max(DATA$K)))
for (ns in 1:sum(SampleNature[1, ])) {
LT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$LT[[ns]][1, ])] <-
DATA$LT[[ns]]
sLT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$sLT[[ns]][1, ])] <-
DATA$sLT[[ns]]
IrrT[seq(CSLengthSample[ns] + 1, CSLengthSample[ns + 1], 1), 1:length(DATA$ITimes[[ns]][1, ])] <-
DATA$ITimes[[ns]]
}
#- THETA matrix ####
if (length(THETA[, 1]) == 0) {
THETA = diag(DATA$ddot_env[2, CSBinPerSample] + (DATA$ddot_env[1, CSBinPerSample]) ^
2 * DATA$dLab[2, CSBinPerSample])
} else{
if (is(THETA)[1] == "character") {
errorMatrix = read.csv(THETA, sep = sepTHETA)
THETA = as.matrix(errorMatrix)
}
}
#--- Index preparation ####
ind_OSL <- which(SampleNature[1,] == 1)
CS_OSL <- cumsum(SampleNature[1,])
ind_C14 <- which(SampleNature[2,] == 1)
CS_C14 <- cumsum(SampleNature[2,])
ind_change <- c(1)
for (i in 2:(Nb_sample - 1)) {
if (SampleNature[1, i] != SampleNature[1, i + 1]) {
ind_change <- c(ind_change, i)
}
}
ind_change <- c(ind_change, Nb_sample)
q <- length(ind_change) %/% 2
r <- length(ind_change) %% 2
##--- description du model BUG ####
BUGModel <- c()
#- Prior
ModelPrior <- 0
data(ModelPrior, envir = environment())
BUGPrior <- c()
if (r == 1) {
if (SampleNature[1, 1] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_OSL)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_C14)
}
if (SampleNature[1, 2] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$OSL_C14)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$C14_OSL)
}
} else{
q <- q - 1
if (q == 0) {
stop(
"[Age_OSLC14()] If you see this message, you are probably trying to run the model with a small number of samples.
You can still use the function, but the C-14 sample cannot be the first sample.",
call. = FALSE
)
}
if (SampleNature[1, 1] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_OSL)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$Sample1_C14)
}
if (SampleNature[1, 2] == 1) {
BUGPrior <- paste(BUGPrior, ModelPrior$OSL_C14)
} else{
BUGPrior <- paste(BUGPrior, ModelPrior$C14_OSL)
}
if (SampleNature[1, Nb_sample] == 1) {
BUGPrior = paste(BUGPrior, ModelPrior$OSL)
} else{
BUGPrior = paste(BUGPrior, ModelPrior$C14)
}
}
#- partie C14
ModelC14 <- 0
data(ModelC14, envir = environment())
if (Model_C14 == "full") {
BUGModel = paste(ModelC14$full, BUGPrior)
} else{
BUGModel = paste(ModelC14$naive, BUGPrior)
}
#- partie OSL
ModelOSL <- 0
data(ModelOSL, envir = environment())
if (LIN_fit == TRUE) {
cLIN = c('LIN')
} else{
cLIN = c()
}
if (Origin_fit == TRUE) {
cO = c("ZO")
} else{
cO = c()
}
Model_GrowthCurve = c(paste("AgesMultiOSL_EXP", cLIN, cO, sep = ""))
BUGModel = c(paste("model{", ModelOSL[[Model_GrowthCurve]][[distribution]], BUGModel, "}"))
if (Model_C14 == "full") {
dataList = list(
'q' = q,
"ind_change" = ind_change,
"ind_OSL" = ind_OSL,
"ind_C14" = ind_C14,
"CS_OSL" = CS_OSL,
"CS_C14" = CS_C14,
'X' = Data_C14Cal,
"sigma" = Data_SigmaC14Cal,
"xTableauCalib" = AgeBP,
"yTableauCalib" = CalC14,
"zTableauCalib" = SigmaCalC14,
'N' = LT,
'sN' = sLT,
"IT" = IrrT,
"sDlab" = DATA$dLab[1, ],
'J' = DATA$J,
'K' = DATA$K,
"ddot" = DATA$ddot_env[1, CSBinPerSample],
"Gamma" = THETA,
"index" = index2,
"BinPerSample" = BinPerSample,
"CSBinPerSample" = CSBinPerSample,
"xbound" = PriorAge,
"StratiConstraints" = StratiConstraints
)
} else{
dataList = list(
'q' = q,
"ind_change" = ind_change,
"ind_OSL" = ind_OSL,
"ind_C14" = ind_C14,
"CS_OSL" = CS_OSL,
"CS_C14" = CS_C14,
'X' = Data_C14Cal,
"sigma" = Data_SigmaC14Cal,
"xTableauCalib" = AgeBP,
"yTableauCalib" = CalC14,
'N' = LT,
'sN' = sLT,
"IT" = IrrT,
"sDlab" = DATA$dLab[1, ],
'J' = DATA$J,
'K' = DATA$K,
"ddot" = DATA$ddot_env[1, CSBinPerSample],
"Gamma" = THETA,
"index" = index2,
"BinPerSample" = BinPerSample,
"CSBinPerSample" = CSBinPerSample,
"xbound" = PriorAge,
"StratiConstraints" = StratiConstraints
)
}
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# JAGS RUN --------------------- START
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
##further settings provided eventually
process_settings <- modifyList(x = list(
max.time = Inf,
interactive = FALSE,
startburnin = 4000,
startsample = 10000,
inits = NA
), val = list(...))
if (autorun == FALSE) {
##as input a text file is wanted, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(BUGModel, con = temp_file)
##run JAGS
results_runjags <-
runjags::run.JAGS(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "Z"),
adapt = adapt,
burnin = burnin,
sample = Iter,
silent.jags = quiet,
method = jags_method,
thin = t
)
}
if (autorun == TRUE) {
##as input a text file is wanted, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(BUGModel, con = temp_file)
##run the auto processing
results_runjags <-
runjags::autorun.jags(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "Z"),
adapt = adapt,
silent.jags = quiet,
method = jags_method,
inits = process_settings$inits,
max.time = process_settings$max.time,
interactive = process_settings$interactive,
startburnin = process_settings$startburnin,
startsample = process_settings$startsample
)
}
results_runjags$args <- list(
"Model_OSL_GrowthCurve" = Model_GrowthCurve,
"Model_OSL_Distribution" = distribution,
"PriorAge" = PriorAge,
"StratiConstraints" = StratiConstraints,
"CovarianceMatrix" = THETA,
"Model_C14" = Model_C14,
"TableauCalib" = TableauCalib,
"SampleNature" = SampleNature,
"Data_C14Cal" = Data_C14Cal,
"Nb_sample" = Nb_sample
)
}
##extract mcmc list
echantillon <- results_runjags$mcmc
U <- summary(echantillon)
##remove mcmc-list from runjags output to reduce output object size
results_runjags$mcmc <- list("MCMC-list is not here. Go to first level -> object named *Sampling*")
##combine chains into one data.frame
Sample <- as.data.frame(runjags::combine.mcmc(echantillon))
## attach "A_" to sample names for the A-parameter
nom = c()
for (i in 1:results_runjags$args$Nb_sample) {
nom = c(nom, paste("A_", SampleNames[i], sep = ""))
}
##plot MCMC
if (SavePdf) {
pdf(file = paste(OutputFilePath, OutputFileName[1], '.pdf', sep = ""))
}
plot_MCMC(echantillon, sample_names = SampleNames)
if (SavePdf) {
dev.off()
}
Outlier <-
SampleNames[ind_C14[which(U$statistics[(results_runjags$args$Nb_sample + 1):(results_runjags$args$Nb_sample + sum(results_runjags$args$SampleNature[2, ])), 1] <
1.5)]]
##- Gelman and Rubin test of convergency of the MCMC
CV = gelman.diag(echantillon, multivariate = FALSE)
cat(paste("\n\n>> Convergence of MCMC for the age parameters <<\n"))
cat("----------------------------------------------\n")
cat(paste("Sample name ", " Bayes estimate ", " Uppers credible interval\n"))
for (i in 1:results_runjags$args$Nb_sample) {
#cat(paste(" Sample name: ", SampleNames[i],"\n"))
#cat("---------------------\n")
cat(paste(
paste("A_", SampleNames[i], sep = ""),
"\t",
round(CV$psrf[i, 1], roundingOfValue),
"\t\t",
round(CV$psrf[i, 2], roundingOfValue),
"\n"
))
}
cat(
"\n\n________________________________________________________________________________\n"
)
cat(" *** WARNING: following informations are only valid if MCMC chains converged ***\n")
cat(
"________________________________________________________________________________\n"
)
# Matrix of results
rnames = c()
for (i in 1:results_runjags$args$Nb_sample) {
rnames = c(rnames, paste("A_", SampleNames[i], sep = ""))
}
R = matrix(
data = NA,
ncol = 8,
nrow = results_runjags$args$Nb_sample,
dimnames = list(
rnames,
c(
"lower bound at 95%",
"lower bound at 68%",
"Bayes estimate",
"upper bound at 68%",
"upper bound at 95%",
"",
"Convergencies: Bayes estimate",
"Convergencies: uppers credible interval"
)
)
)
##- Bayes estimate and credible interval
cat(paste(
"\n\n>> Bayes estimates of Age for each sample and credible interval <<\n"
))
AgePlot95 = matrix(data = NA, nrow = results_runjags$args$Nb_sample, ncol = 3)
AgePlot68 = matrix(data = NA, nrow = results_runjags$args$Nb_sample, ncol = 3)
AgePlotMoy = rep(0, results_runjags$args$Nb_sample)
for (i in 1:results_runjags$args$Nb_sample) {
cat("------------------------------------------------------\n")
#cat(paste(" Sample name: ", SampleNames[i],"\n"))
#cat("---------------------\n")
cat(paste("Sample name", "\t", "Bayes estimate", " Credible interval: \n"))
cat(paste(
paste("A_", SampleNames[i], sep = ""),
"\t",
round(mean(Sample[, i]), roundingOfValue),
'\n'
))
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
HPD_95 = CredibleInterval(Sample[, i], 0.95, roundingOfValue =
roundingOfValue)
HPD_68 = CredibleInterval(Sample[, i], 0.68, roundingOfValue =
roundingOfValue)
cat(
"\t\t\t\t at level 95% \t",
round(c(HPD_95[2]), roundingOfValue),
"\t\t",
round(c(HPD_95[3]), roundingOfValue),
"\n"
)
cat(
"\t\t\t\t at level 68% \t",
round(c(HPD_68[2]), roundingOfValue),
"\t\t",
round(c(HPD_68[3]), roundingOfValue),
"\n"
)
AgePlot95[i, ] = HPD_95
AgePlot68[i, ] = HPD_68
AgePlotMoy[i] = round(mean(Sample[, i]), roundingOfValue)
R[i, 3] = round(mean(Sample[, i]), roundingOfValue)
R[i, c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
R[i, c(2, 4)] = round(HPD_68[2:3], roundingOfValue)
R[i, 6] = c('')
R[i, 7] = round(CV$psrf[i, 1], roundingOfValue)
R[i, 8] = round(CV$psrf[i, 2], roundingOfValue)
}
cat("\n------------------------------------------------------\n")
R[, c(7, 8)] <- round(CV$psrf[1:results_runjags$args$Nb_sample, ], roundingOfValue)
# Representation graphique des resultats
# des HPD sur la courbe de calibration
if (sum(results_runjags$args$SampleNature[2, ]) > 1) {
couleur = rainbow(results_runjags$args$Nb_sample)
par(mfrow = c(1, 1),
las = 0,
mar = c(5, 5, 2, 2))
xl = c(min(results_runjags$args$PriorAge[seq(1, (2 * results_runjags$args$Nb_sample - 1), 2)]), max(results_runjags$args$PriorAge[seq(2, (2 *
results_runjags$args$Nb_sample), 2)]))
plot(
xl,
xl,
col = "white",
xlab = c("Age (in ka)"),
ylab = c("cal C14"),
xaxt = "n",
yaxt = "n",
cex.lab = 1.8
)
axis(2, cex.axis = 2)
axis(1, cex.axis = 2)
polygon(c(AgeBP, rev(AgeBP)),
c(CalC14 + 2 * SigmaCalC14, rev(CalC14 - 2 * SigmaCalC14)),
col = "gray",
border = "black")
for (i in ind_C14) {
lines(c(AgePlot95[i, 2:3]),
rep(results_runjags$args$Data_C14Cal[CS_C14[i]], 2),
col = couleur[i],
lwd = 4)
lines(
AgePlotMoy[i],
results_runjags$args$Data_C14Cal[CS_C14[i]],
col = "black",
lwd = 2,
type = 'p'
)
}
legend(
"topleft",
SampleNames[ind_C14],
lty = rep(1, results_runjags$args$Nb_sample),
lwd = rep(2, results_runjags$args$Nb_sample),
cex = 1,
col = couleur[ind_C14]
)
if (SavePdf == TRUE) {
dev.print(
pdf,
file = paste(OutputFilePath, OutputFileName[2], '.pdf', sep = ""),
width = 8,
height = 10
)
}
}
##CSV output
if (SaveEstimates == TRUE) {
write.csv(R, file = c(
paste(OutputTablePath, "Estimates", OutputTableName, ".csv", sep = "")
))
}
.list_BayLum <- function(..., originator = sys.call(which = -1)[[1]]){
## set list
l <- list(...)
## update originators
attr(l, "class") <- "BayLum.list"
attr(l, "originator") <- as.character(originator)
return(l)
}
# Create return object -------------------------------------------------------------------------
output <- .list_BayLum(
"Ages" = data.frame(
SAMPLE = SampleNames,
AGE = AgePlotMoy,
HPD68.MIN = AgePlot68[, 2],
HPD68.MAX = AgePlot68[, 3],
HPD95.MIN = AgePlot95[, 2],
HPD95.MAX = AgePlot95[, 3],
stringsAsFactors = FALSE
),
"Sampling" = echantillon,
"PriorAge" = results_runjags$args$PriorAge,
"StratiConstraints" = results_runjags$args$StratiConstraints,
"Model_OSL_GrowthCurve" = results_runjags$args$Model_OSL_GrowthCurve,
"Model_OSL_Distribution" = results_runjags$args$Model_OSL_Distribution,
"CovarianceMatrix" = results_runjags$args$CovarianceMatrix,
"Model_C14" = results_runjags$args$Model_C14,
"CalibrationCurve" = results_runjags$args$CalibrationCurve,
"Outlier" = Outlier,
"runjags_object" = results_runjags
)
# Plot ages -----------------------------------------------------------------------------------
plot_Ages(object = output, legend.pos = "bottomleft")
##TODO: get rid of this
if (SavePdf) {
dev.print(
pdf,
file = paste(OutputFilePath, OutputFileName[3], '.pdf', sep = ""),
width = 8,
height = 10
)
}
# Return output -------------------------------------------------------------------------------
return(output)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.