Nothing
#' Bayesian analysis for the OSL age estimation of one sample
#'
#' This function computes the age (in ka) of a sample according to the model developed in Combes and Philippe (2017),
#' based on an output of [Generate_DataFile] or [Generate_DataFile_MG].\cr
#' A sample, for which data is available in several BIN files, can be analysed.
#'
#' @param DATA [list] of objects: `LT`, `sLT`, `ITimes`, `dLab`, `ddot_env`, `regDose`, `J`, `K`, `Nb_measurement`,
#' provided by the function [Generate_DataFile] or [Generate_DataFile_MG].
#'
#' `DATA` can contain information for more than one sample.
#'
#' @param SampleName [character]: name of the sample.
#'
#' @param PriorAge [numeric] (with default): lower and upper bounds for the sample age parameter (in ka).
#' Note that, `length(PriorAge)=2`.
#'
#' @param BinPerSample [integer] (with default): vector with the number of BIN files per sample.
#' If in `DATA` there is more than one sample,
#' the `BinPerSample` vector must be the same as that used to run the function
#' [Generate_DataFile] or in [Generate_DataFile_MG] for generating the `DATA` object.
#'
#' @param SavePdf [logical] (with default): if TRUE save graph 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 = 2`, 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 and credible interval at level 68% and 95% and
#' the result of the gelman en Rubin test of convergency, 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 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 I [integer] (with default): if `DATA` contains data from more than one sample,
#' I indicates the ID number of the sample to be analysed.
#'
#' @param Iter [integer] (with default): number of iterations for the MCMC computation (for more information see [jags.model]).
#'
#' @param t [integer] (with default): 1 every `t` iterations of the MCMC is considered for sampling the posterior distribution
#' (for more information see [jags.model]).
#'
#' @param n.chains [integer] (with default): number of independent chains for the model (for more information see [jags.model]).
#'
#' @param quiet [logical] (with default): enables/disables \link{rjags} messages
#'
#' @param roundingOfValue [integer] (with default): Integer indicating the number of decimal places to be used, default = 3.
#'
#' @details
#'
#' **Option on growth curves**\cr
#'
#' As for [AgeS_Computation] and [Palaeodose_Computation], the user can choose from 4 dose response curves:
#' \itemize{
#' \item **Saturating exponential plus linear growth** (`AgeMultiBF_EXPLIN`):
#'
#' for all `x` in IR+, \eqn{f(x)=a(1-exp(-x/b))+cx+d}; select
#' \itemize{
#' \item `LIN_fit=TRUE`
#' \item `Origin_fit=FALSE`
#' }
#' \item **Saturating exponential growth** (`AgeMultiBF_EXP`):
#'
#' for all `x` in IR+, \eqn{f(x)=a(1-exp(-x/b))+d}; select
#' \itemize{
#' \item `LIN_fit = FALSE`
#' \item `Origin_fit = FALSE`
#' }
#' \item **Saturating exponential plus linear growth and fitting through the origin** (`AgeMultiBF_EXPLINZO`):
#'
#' for all `x` in IR+, \eqn{f(x)=a(1-exp(-x/b))+cx}; select
#' \itemize{
#' \item `LIN_fit=TRUE`
#' \item `Origin_fit=TRUE`
#' }
#' \item **Saturating exponential growth and fitting through the origin** (`AgeMultiBF_EXPZO`):
#'
#' for all `x` in IR+, \eqn{f(x)=a(1-exp(-x/b))}; select
#' \itemize{
#' \item `LIN_fit=FALSE`
#' \item `Origin_fit=TRUE`
#' }
#' }
#'
#' **Option on equivalent dose distribution around the palaeodose**\cr
#'
#' The use can choose between :
#' \itemize{
#' \item `cauchy`: a Cauchy distribution with location parameter equal to the palaeodose of the sample
#' \item `gaussian`: a Gaussian distribution with mean equal to the palaeodose of the sample
#' \item `lognormal_A`: a log-normal distribution with mean or **A**verage equal to the palaeodose of the sample
#' \item `lognormal_M`: a log-normal distribution with **M**edian equal to the palaeodose of the sample
#' }
#'
#' @return
#' **NUMERICAL OUTPUT**
#'
#' \enumerate{
#' \item **A list containing the following objects**:
#' \itemize{
#' \item **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.
#' \item **Model_GrowthCurve**, stating which dose response fitting option was chosen;
#' \item **Distribution**, stating which distribution was chosen to model the dispersion of
#' individual equivalent dose values around the palaeodose of the sample;
#' \item **PriorAge**, stating the priors used for the age parameter (in ka).
#' }
#' \item **The Gelman and Rubin test of convergency**: prints the result of the Gelman and Rubin test of convergency for the age, palaeodose and equivalent dose dispersion parameters.
#' A result close to one is expected.\cr
#' In addition, the user must visually assess the convergency of the trajectories by looking at the graph
#' generated by the function (see **PLOT OUTPUT** for more informations).\cr
#' If both convergencies (Gelman and Rubin test and plot checking) are satisfactory,
#' the user can consider the printed estimates as valid. Otherwise, the user may try increasing the number of MCMC interations
#' (`Iter`), or being more precise on the `PriorAge` parameter (for example specify if it is a young sample `c(0.01,10)` an old sample `c(10,100)`),
#' or changing the parameter `distribution` or the growth curve, to reach convergency.to reach convergency.
#' \item **Credible intervals and Bayes estimates**: prints the Bayes esitmates, the credible intervals at 95% and 68% for
#' the age, palaeodose and equivalent dose dispersion parameters of the sample.
#' }
#'
#' **PLOT OUTPUT**
#'
#' A graph with the MCMC trajectories and posterior distributions of the age, palaeodose and equivalent dose dispersion parameters is displayed.\cr
#' The first line of the figure correponds to the age parameter, the second to the palaeodose parameter and the third to the equivalent dose dispersion parameter.
#' On each line, the plot on the left represents the MCMC trajectories, and the one on the right the posterior distribution of the parameter.
#'
#' To give the results in a publication, we recommend to give the Bayes estimate of the parameter as well as the credible interval at 95% or 68%.
#'
#' @author Claire Christophe, Sebastian Kreutzer, Anne Philippe, Guillaume Guérin
#'
#' @seealso [Generate_DataFile], [Generate_DataFile_MG], [rjags], [plot_MCMC], [AgeS_Computation], [Palaeodose_Computation]
#'
#' @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
#'
#' @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.
#'
#' @examples
#' ## load data file generated by the function Generate_DataFile
#' data(DATA1,envir = environment())
#' priorage <- c(10,60) # GDB3 is an old sample
#' Age <- Age_Computation(
#' DATA = DATA1,
#' SampleName = "GDB3",
#' PriorAge = priorage,
#' Iter = 100,
#' quiet = TRUE)
#'
#' @md
#' @export
Age_Computation <- function(
DATA,
SampleName = DATA$SampleNames[1],
PriorAge = c(0.01, 100),
BinPerSample = c(1),
SavePdf = FALSE,
OutputFileName = c("MCMCplot"),
OutputFilePath = c(''),
SaveEstimates = FALSE,
OutputTableName = c("DATA"),
OutputTablePath = c(''),
LIN_fit = TRUE,
Origin_fit = FALSE,
distribution = c("cauchy"),
I = 1,
Iter = 50000,
t = 5,
n.chains = 3,
quiet = FALSE,
roundingOfValue = 3
){
Model_Age<-0
data(Model_Age,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("AgeMultiBF_EXP",cLIN,cO,sep=""))
CSBinPerSample=cumsum(BinPerSample)
if(BinPerSample[I]==1){
index2=0
}else{
index2=rep(0,BinPerSample[I])
index2[2:BinPerSample[I]]=cumsum(DATA$J[(CSBinPerSample[I]-BinPerSample[I]+1):(CSBinPerSample[I]-1)])
}
dataList = list('N'= DATA$LT[[I]],'sN'=DATA$sLT[[I]],"IT"=DATA$ITimes[[I]],
"sDlab"=DATA$dLab[1,(CSBinPerSample[I]-BinPerSample[I]+1):(CSBinPerSample[I])],
'J'=DATA$J[(CSBinPerSample[I]-BinPerSample[I]+1):(CSBinPerSample[I])],
'K'=DATA$K[(CSBinPerSample[I]-BinPerSample[I]+1):(CSBinPerSample[I])],
"ddot"=DATA$ddot_env[1,(CSBinPerSample[I]-BinPerSample[I]+1)],
"Sigma"=DATA$ddot_env[2,(CSBinPerSample[I]-BinPerSample[I]+1)]+DATA$ddot_env[1,(CSBinPerSample[I]-BinPerSample[I]+1)]^2*DATA$dLab[2,(CSBinPerSample[I]-BinPerSample[I]+1)],
"xbound"=PriorAge,"index"=index2,"BinPerSample"=BinPerSample[I])
##open text connection
con <- textConnection(Model_Age[[Model_GrowthCurve]][[distribution]])
jags <- rjags::jags.model(
file = con,
data = dataList,
n.chains = n.chains,
n.adapt = Iter,
quiet = quiet
)
##close connection
close(con)
##set progress.bar
if(quiet) progress.bar <- 'none' else progress.bar <- 'text'
update(jags,Iter, progress.bar = progress.bar)
echantillon <-
rjags::coda.samples(jags,
c("A", "D", "sD"),
min(Iter, 10000),
thin = t,
progress.bar = progress.bar)
CV <- gelman.diag(echantillon)
sample <- echantillon[[1]]
for(i in 2:n.chains){
sample=rbind(sample,echantillon[[i]])
}
##MCMC plot output
plot_MCMC(echantillon, sample_names = SampleName)
if(SavePdf){
dev.print(pdf,file=paste(OutputFilePath,OutputFileName,'.pdf',sep=""),width=8,height=10)
}
cat(paste("\n\n>> Sample name <<\n"))
cat("----------------------------------------------\n")
cat(paste(SampleName))
cat(paste("\n\n>> Results of the Gelman and Rubin criterion of convergence <<\n"))
cat("----------------------------------------------\n")
cat(paste("\t", "Point estimate", "Uppers confidence interval\n"))
cat(paste("A\t",round(CV$psrf[1,1],roundingOfValue),"\t\t",round(CV$psrf[1,2],roundingOfValue),"\n"))
cat(paste("D\t",round(CV$psrf[2,1],roundingOfValue),"\t\t",round(CV$psrf[2,2],roundingOfValue),"\n"))
cat(paste("sD\t",round(CV$psrf[3,1],roundingOfValue),"\t\t",round(CV$psrf[3,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")
# Matrix of results
R=matrix(data=NA,ncol=8,nrow=3,dimnames=list(c("A","D","sD"),
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")))
cat(paste("parameter", "\t","Bayes estimate","\t"," Credible interval \n"))
cat("----------------------------------------------\n")
cat(paste("A","\t\t", round(mean(sample[,1]),3),'\n'))
HPD_95=CredibleInterval(echantillon[[1]][,1],0.95,roundingOfValue = roundingOfValue)
HPD_68=CredibleInterval(echantillon[[1]][,1],0.68,roundingOfValue = roundingOfValue)
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
cat("\t\t\t\t at level 95%\t",round(c(HPD_95[2]),2),"\t\t",round(c(HPD_95[3]),roundingOfValue),"\n")
cat("\t\t\t\t at level 68%\t",round(c(HPD_68[2]),2),"\t\t",round(c(HPD_68[3]),roundingOfValue),"\n")
cat("----------------------------------------------\n")
R[1,3]=round(mean(sample[,1]),roundingOfValue)
R[1,c(1,5)]=round(HPD_95[2:3],roundingOfValue)
R[1,c(2,4)]=round(HPD_68[2:3],roundingOfValue)
cat(paste("D","\t\t", round(mean(sample[,2]),roundingOfValue),'\n'))
HPD_95=CredibleInterval(echantillon[[1]][,2],0.95,roundingOfValue = roundingOfValue)
HPD_68=CredibleInterval(echantillon[[1]][,2],0.68,roundingOfValue = roundingOfValue)
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
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,3]=round(mean(sample[,2]),roundingOfValue)
R[2,c(1,5)]=round(HPD_95[2:3],roundingOfValue)
R[2,c(2,4)]=round(HPD_68[2:3],roundingOfValue)
cat("----------------------------------------------\n")
cat(paste("sD","\t\t", round(mean(sample[,3]),3),'\n'))
HPD_95=CredibleInterval(echantillon[[1]][,3],0.95,roundingOfValue = roundingOfValue)
HPD_68=CredibleInterval(echantillon[[1]][,3],0.68,roundingOfValue = roundingOfValue)
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
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[3,3]=round(mean(sample[,3]),roundingOfValue)
R[3,c(1,5)]=round(HPD_95[2:3],roundingOfValue)
R[3,c(2,4)]=round(HPD_68[2:3],roundingOfValue)
R[,6]=c('','','')
R[,7]=round(CV$psrf[,1],roundingOfValue)
R[,8]=round(CV$psrf[,2],roundingOfValue)
if(SaveEstimates==TRUE){
write.csv(R,file=c(paste(OutputTablePath,"Estimates",OutputTableName,".csv",sep="")))
}
Info=list("Sampling"=echantillon,"Model_GrowthCurve"=Model_GrowthCurve, "Distribution"=distribution,"PriorAge"=PriorAge)
return(Info)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.