#' @title Bayesian analysis for OSL age estimation of various samples
#'
#' @description 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].\cr
#' Samples, for which data is available in several BIN files, can be analysed.\cr
#' Single-grain or Multi-grain OSL measurements can be analysed simultaneously.
#'
#' @param 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].
#' \code{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.
#'
#' @param SampleNames [character] vector: names of samples. The length of this vector is equal to `Nb_sample`.
#'
#' @param Nb_sample [integer]: number of samples, `Nb_sample>1`.
#'
#' @param PriorAge [numeric] vector (with default): lower and upper bounds for age parameter of each sample (in ka).
#' Note that, \code{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`.
#'
#' @param 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].
#'
#' @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`;
#' \code{length(OutputFileName)=2}, see \bold{PLOT OUTPUT} in \bold{Value} section for more informations.
#'
#' @param OutputFilePath [character] (with default): path to the pdf file that will be generated by the function if \code{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 \code{OutputFileName} in folder \code{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 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 \code{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.
#'
#' @param sepTHETA [character] (with default): if `THETA` is character,
#' indicate column separator in `THETA` CSV-file.
#'
#' @param 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.
#'
#' @param sepSC [character] (with default): if \code{StratiConstraints} is character,
#' indicate column separator in \code{StratiConstraints} .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.
#' See details section for more informations on the proposed dose response curves.
#'
#' @param 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.
#'
#' @param 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.
#'
#' @param 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.
#'
#' @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 which will be used to assess convergence and ages (see [runjags::run.jags]).
#'
#' @param t [integer] (with default): 1 every \code{t} iterations of the MCMC is considered for sampling the posterior distribution.
#' (for more information see [runjags::run.jags]).
#'
#' @param n.chains [integer] (with default): number of independent chains for the model (for more information see [runjags::run.jags]).
#'
#' @param 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])
#'
#' @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. 1) When creating a new JAGS model, see details
#' for supported arguments. 2) If extending a JAGS model see arguments of [runjags::extend.JAGS].
#'
#'
#' @details **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 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 \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+A,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) functions
#' 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 take care about the separator 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 follows.
#' \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 symetric matrix.
#'
#' The user can also refer to a .csv file that contains the errors as defined above.\cr
#'
#' Alternatively you can use the function [create_ThetaMatrix].
#'
#' **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.
#' }
#'
#' @return
#' **NUMERICAL OUTPUT**
#'
#' \enumerate{
#' \item **A list of type `BayLum.list` containing the following objects:**
#' \itemize{
#' \item \bold{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;
#' \item \bold{Model_GrowthCurve}: stating which dose response fitting option was chosen;
#' \item \bold{Distribution}: stating which distribution was chosen to model the dispersion of
#' individual equivalent dose values around the palaeodose of the sample;
#' \item \bold{PriorAge}: stating the priors used for the age parameter (in ka);
#' \item \bold{StratiConstraints}: stating the stratigraphic relations between samples considered in the model;
#' \item \bold{CovarianceMatrix}: stating the covariance matrix of error used in the model, highlighting common errors between samples or not.
#' \item \bold{model}: returns the model that was used for the Bayesian modelling as a [character]
#' \item \bold{JAGS model output}: returns the JAGS model with class "runjags".
#' }
#' \item\bold{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.\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 being more precise on the \code{PriorAge} parameter (for example specify if it is a young sample \code{c(0.01,10)} an old sample \code{c(10,100)}),
#' or changing the parameter \code{distribution} or the growth curve, to reach convergence.
#' \item \bold{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**
#'
#' \enumerate{
#' \item\bold{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. \cr
#' 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.
#' \item \bold{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%.
#'
#' @author Claire Christophe, Anne Philippe, Guillaume Guérin, Sebastian Kreutzer
#'
#' @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 [Generate_DataFile], [Generate_DataFile_MG], [rjags::rjags-package], [plot_MCMC], [SCMatrix], [Age_Computation], [Palaeodose_Computation], [plot_Ages], [create_ThetaMatrix], [runjags::autorun.jags]
#'
#' @references
#' 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
#'
#' @examples
#' ## 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)
#'
#'\dontrun{
#' ## 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")
#' }
#'
#' @md
#' @export
AgeS_Computation <- function(
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,
...
) {
#---check to see if DATA input is a runjags-object and extend if so ####
if (inherits(DATA, "BayLum.list")) {
# reattach mcmc-list which was removed to reduce object size
DATA$runjags_object$mcmc <- DATA$Sampling
# extend via runjags
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 original BayLum run (as to not lose them when results are processed)
results_runjags$args <- list(
"Model_GrowthCurve" = DATA$runjags_object$args$Model_GrowthCurve,
"Distribution" = DATA$runjags_object$args$Distribution,
"PriorAge" = DATA$runjags_object$args$PriorAge,
"StratiConstraints" = DATA$runjags_object$args$StratiConstraints,
"CovarianceMatrix" = DATA$runjags_object$args$CovarianceMatrix,
"model" = DATA$runjags_object$args$model
)
}
#---check to see if DATA input is a DataFile and run JAGS ####
if (!inherits(DATA, "BayLum.list")) {
##---Index preparation ####
CSBinPerSample <- cumsum(BinPerSample)
LengthSample <- c()
for (ns in 1:Nb_sample) {
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:Nb_sample) {
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 ####
## cover character case
if (is(THETA)[1] == "character")
THETA <- as.matrix(read.csv(THETA, sep = sepTHETA, header = FALSE))
## cover data.frame case (users should not do it but they do)
if (inherits(THETA, "data.frame"))
THETA <- as.matrix(THETA)
if (length(THETA[,1]) == 0) {
THETA <- diag(DATA$ddot_env[2, CSBinPerSample] +
(DATA$ddot_env[1, CSBinPerSample]) ^ 2 * DATA$dLab[2, CSBinPerSample])
}
##JAGS will crash with a runtime error if the dimension of the theta matrix does not match the
##number of samples
if (sum(dim(THETA)) %% Nb_sample != 0)
stop(
"[AgeS_Computation()] The number of samples does not match the dimension of the THETA-matrix!",
call. = FALSE)
##---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)
}
}
##---BUG file selection ####
Model_AgeS <- 0
data(Model_AgeS, 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("AgesMultiCS2_EXP", cLIN, cO, sep = ""))
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
##---JAGS RUN --------------------- START ####
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
##set variables
dataList <- list(
'N' = LT,
'sN' = sLT,
"IT" = IrrT,
"sDlab" = DATA$dLab[1,],
'J' = DATA$J,
'K' = DATA$K,
"I" = Nb_sample,
"ddot" = DATA$ddot_env[1, CSBinPerSample],
"Gamma" = THETA,
"xbound" = PriorAge,
"StratiConstraints" = StratiConstraints,
"index" = index2,
"BinPerSample" = BinPerSample,
"CSBinPerSample" = CSBinPerSample
)
##select model
if (is.null(model))
model <- Model_AgeS[[Model_GrowthCurve]][[distribution]]
##there are two main ways of running JAGS: single-run vs auto-run
##(1) if user selects to do single-run:
if (!autorun) {
##a text file is wanted as input, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(model, con = temp_file)
##run JAGS
results_runjags <-
runjags::run.JAGS(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "D", "sD"),
adapt = adapt,
burnin = burnin,
sample = Iter,
silent.jags = quiet,
method = jags_method,
thin = t
)
}
##(2) if user selects to do auto-run:
if (autorun) {
##further settings provided eventually
process_settings <- modifyList(x = list(
max.time = Inf,
interactive = FALSE,
startburnin = 4000,
startsample = 10000,
inits = NA
), val = list(...))
##a text file is wanted as input, so we have to cheat a little bit
temp_file <- tempfile(fileext = ".txt")
writeLines(model, con = temp_file)
##run the autoprocessing
results_runjags <-
runjags::autorun.JAGS(
model = temp_file,
data = dataList,
n.chains = n.chains,
monitor = c("A", "D", "sD"),
adapt = adapt,
startburnin = process_settings$startburnin,
startsample = process_settings$startsample,
silent.jags = quiet,
method = jags_method,
thin = t,
inits = process_settings$inits,
max.time = process_settings$max.time,
interactive = process_settings$interactive
)
}
# storing the arguments used for the BayLum-run this way,
# because it allows us an easy way to code the storage of arguments when extending a JAGS-model.
results_runjags$args <- list(
"Model_GrowthCurve" = Model_GrowthCurve,
"Distribution" = distribution,
"PriorAge" = PriorAge,
"StratiConstraints" = StratiConstraints,
"CovarianceMatrix" = THETA,
"model" = model
)
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# JAGS RUN --------------------- END
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
}
#---processing of JAGS results
##extract mcmc list from runjags object
echantillon <- results_runjags$mcmc
##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))
#---plot MCMC ####
if (SavePdf) {
pdf(file = paste(OutputFilePath, OutputFileName[1], '.pdf', sep = ""))
}
##try makes sure that the function runs
try(plot_MCMC(echantillon, sample_names = SampleNames))
if (SavePdf) {
dev.off()
}
#---Gelman and Rubin test of convergence of the MCMC ####
CV <- gelman.diag(echantillon, multivariate = FALSE)
cat(paste(
"\n\n>> Results of the Gelman and Rubin criterion of convergence <<\n"
))
for (i in 1:Nb_sample) {
cat("----------------------------------------------\n")
cat(paste(" Sample name: ", SampleNames[i], "\n"))
cat("---------------------\n")
cat(paste("\t\t", "Point estimate", "Uppers confidence interval\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(paste(
paste("D_", SampleNames[i], sep = ""),
"\t",
round(CV$psrf[(Nb_sample + i), 1], roundingOfValue),
"\t\t",
round(CV$psrf[(Nb_sample + i), 2], roundingOfValue),
"\n"
))
cat(paste(
paste("sD_", SampleNames[i], sep = ""),
"\t",
round(CV$psrf[(2 * Nb_sample + i), 1], roundingOfValue),
"\t\t",
round(CV$psrf[(2 * Nb_sample + i), 2], roundingOfValue),
"\n"
))
}
cat(
"\n\n---------------------------------------------------------------------------------------------------\n"
)
cat(
" *** WARNING: The following information are only valid if the MCMC chains have converged ***\n"
)
cat(
"---------------------------------------------------------------------------------------------------\n\n"
)
#---print results ####
##Matrix of results
rnames <- rep(NA, 3 * Nb_sample)
for (i in 1:Nb_sample) {
rnames[i] = c(paste("A_", SampleNames[i], sep = ""))
rnames[Nb_sample + i] = c(paste("D_", SampleNames[i], sep = ""))
rnames[2 * Nb_sample + i] = c(paste("sD_", SampleNames[i], sep = ""))
}
R = matrix(
data = NA,
ncol = 8,
nrow = 3 * 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: Point estimate",
"Convergencies: uppers confidence interval"
)
)
)
##Bayes estimate and credible interval
cat(
paste(
"\n\n>> Bayes estimates of Age, Palaeodose and its dispersion for each sample and credible interval <<\n"
)
)
AgePlot95 = matrix(data = NA,
nrow = Nb_sample,
ncol = 3)
AgePlot68 = matrix(data = NA,
nrow = Nb_sample,
ncol = 3)
AgePlotMoy = rep(0, Nb_sample)
for (i in 1:Nb_sample) {
cat("----------------------------------------------\n")
cat(paste(" Sample name: ", SampleNames[i], "\n"))
cat("---------------------\n")
cat(paste(
"Parameter",
"\t",
"Bayes estimate",
"\t",
" 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)
cat(paste(
"\nParameter",
"\t",
"Bayes estimate",
"\t",
" Credible interval \n"
))
cat(paste(
paste(" D_", SampleNames[i], sep = ""),
"\t",
round(mean(sample[, (Nb_sample + i)]), roundingOfValue),
'\n'
))
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
HPD_95 = CredibleInterval(sample[, (Nb_sample + i)], 0.95, roundingOfValue =
roundingOfValue)
HPD_68 = CredibleInterval(sample[, (Nb_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"
)
R[(Nb_sample + i), 3] = round(mean(sample[, (Nb_sample + i)]), roundingOfValue)
R[(Nb_sample + i), c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
R[(Nb_sample + i), c(2, 4)] = round(HPD_68[2:3], roundingOfValue)
cat(paste(
"\nParameter",
"\t",
"Bayes estimate",
"\t",
" Credible interval \n"
))
cat(paste(
paste("sD_", SampleNames[i], sep = ""),
"\t",
round(mean(sample[, (2 * Nb_sample + i)]), roundingOfValue),
'\n'
))
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
HPD_95 = CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
i)], 0.95, roundingOfValue = roundingOfValue)
HPD_68 = CredibleInterval(echantillon[[1]][, (2 * Nb_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"
)
R[(2 * Nb_sample + i), 3] = round(mean(sample[, (2 * Nb_sample + i)]), roundingOfValue)
R[(2 * Nb_sample + i), c(1, 5)] = round(HPD_95[2:3], roundingOfValue)
R[(2 * Nb_sample + i), c(2, 4)] = round(HPD_68[2:3], roundingOfValue)
# R[(3*(i-1)+1):(3*i),6]=c('','','')
# R[(3*(i-1)+1):(3*i),7]=round(CV$psrf[c(i,i+Nb_sample,i+2*Nb_sample),1],2)
# R[(3*(i-1)+1):(3*i),8]=round(CV$psrf[c(i,i+Nb_sample,i+2*Nb_sample),2],2)
}
cat("\n----------------------------------------------\n")
R[, c(7, 8)] <- round(CV$psrf, roundingOfValue)
##check if there is stratigraphic order between sample.
##if(sum(StratiConstraints[2:Nb_sample,1:Nb_sample]>0)){
##}
#---print csv table ####
if (SaveEstimates == TRUE) {
write.csv(R, file = c(
paste(
OutputTablePath,
"Estimates",
OutputTableName,
".csv",
sep = ""
)
))
}
#---Create return object ####
.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)
}
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,
"Model_GrowthCurve" = results_runjags$args$Model_GrowthCurve,
"Distribution" = results_runjags$args$Distribution,
"PriorAge" = results_runjags$args$PriorAge,
"StratiConstraints" = results_runjags$args$StratiConstraints,
"CovarianceMatrix" = results_runjags$args$CovarianceMatrix,
"model" = results_runjags$model,
"runjags_object" = results_runjags
)
#---Plot ages ####
BayLum::plot_Ages(object = output, legend.pos = "bottomleft")
##TODO: get rid of this ... at some point
if (SavePdf) {
dev.print(
pdf,
file = paste(OutputFilePath, OutputFileName[3], '.pdf', sep = ""),
width = 8,
height = 10
)
}
#---Return output ####
return(output)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.