R/Age_Computation.R

In BayLum: Chronological Bayesian Models Integrating Optically Stimulated Luminescence and Radiocarbon Age Dating

#' 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,
  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=ArchaeoPhases::CredibleInterval(echantillon[[1]][,1],0.95,roundingOfValue = roundingOfValue)
  HPD_68=ArchaeoPhases::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=ArchaeoPhases::CredibleInterval(echantillon[[1]][,2],0.95,roundingOfValue = roundingOfValue)
  HPD_68=ArchaeoPhases::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=ArchaeoPhases::CredibleInterval(echantillon[[1]][,3],0.95,roundingOfValue = roundingOfValue)
  HPD_68=ArchaeoPhases::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)

}