R/AgeC14_Computation.R

In R-Lum/BayLum: Chronological Bayesian Models Integrating Optically Stimulated Luminescence and Radiocarbon Age Dating

#' @title Bayesian analysis for C-14 age estimations of various samples
#'
#' @description This function calibrates the C-14 age of samples to get an age (in ka).
#' The user can choose one of the following radiocarbon calibration curve:
#' Northern or Southern Hemisphere or marine atmospheric. It must be the same curve for all samples.
#'
#' @param Data_C14Cal [numeric] (**required**): corresponding to C-14 age estimate.
#'
#' @param Data_SigmaC14Cal [numeric] (**required**): corresponding to the error of C-14 age estimates.
#'
#' @param SampleNames [character] (**required**): names of sample. The length of this vector is equal to `Nb_sample`.
#'
#' @param Nb_sample [integer]: number of samples.
#'
#' @param PriorAge [numeric] (with default): lower and upper bounds for age parameter of each sample in years (not in ka).
#' Note that, `length(PriorAge) == 2 * Nb_sample`
#' and `PriorAge[2i-1,2i]` corresponds to the lower and upper bounds of sample whose number ID is equal to `i`.
#'
#' @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 `SavePd=TRUE`,
#' `length(OutputFileName) = =3`, 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 `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 and Rubin test of convergence, in a csv table named `OutputFileName` in folder `OutputFilePath`.
#'
#' @param OutputTableName [logical] (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 [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 `StratiConstraints` is character,
#' indicate column separator in `StratiConstraints` csv file.
#'
#' @param Model [character] (with default): if \bold{"full"}, error on estimate calibration curve is taken account.
#' If \bold{"naive"} this error is not taken account in the age estimate.
#'
#' @param CalibrationCurve [character] (with default): calibration curve chosen. 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 Iter [integer] (with default): number of iterations for the MCMC computation (for more information see [rjags::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 [rjags::jags.model].
#'
#' @param n.chains [integer] (with default): number of independent chains for the model (for more information see [rjags::jags.model].
#'
#' @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 set to 3.
#'
#' @details
#'
#' **How to fill `StratiConstraints`?**
#'
#' If there is stratigraphic relations between samples, *C-14 age in `Data_C14Cal` must be ordered by order of increasing ages*.
#'
#' The user can fill the `StratiConstraints` matrix as follow.
#'\enumerate{
#'  \item **Size of the matrix**: row number of `StratiConstraints` matrix is equal to `Nb_sample+1`,
#'  and column number is equal to `Nb_sample`.
#'  \item **First line of the matrix**:
#' for all `i in {1,...,Nb_Sample}`, `StratiConstraints[1,i]=1` that means the lower bound
#' of the sample age (given in `PriorAge[2i-1]`)
#' for the sample whose number ID is equal to `i`, is taken into account.
#'  \item **Sample relations**: for all  `j in {2,...,Nb_Sample+1}` and all `i in {j,...,Nb_Sample}`,
#' `StratiConstraints[j,i]=1` if sample age whose number ID is equal to `j-1` is lower than
#' sample age whose number ID is equal to `i`. Otherwise, `StratiConstraints[j,i]=0`.
#' }
#'
#' Note that `StratiConstraints_{2:Nb_sample+1,1:Nb_sample}` is a upper triangular matrix.
#'
#' The user can also use `SCMatrix` or [SC_Ordered] (if all samples are ordered) functions
#' to construct the `StratiConstraints` matrix.
#'
#' The user can also refer to a .csv file that contains the relation between samples as defined above.
#' The user must take care about the separator used in the csv file using the argument `sepSC`.\cr
#'
#' ** More precision on `Model` **\cr
#'
#' We propose two models "full" or "naive". If `Model = 'full'` that means
#' measurement error and error on calibration curve are taken account in
#' the Bayesian model; if `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}
#'
#' \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;
#'   \item \bold{Outlier}: stating the names of samples that are considered as outliers;
#'   \item \bold{Model}: stating which model was chosen (`"full"` or `"naive"`);
#'   \item \bold{CalibrationCurve}: stating which radiocarbon calibration curve was chosen;
#'   \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{The Gelman and Rubin test of convergency}: print 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 (`Iter`)
#' or being more precise if it is possible on the `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.
#' }
#'
#' \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{Summary of sample age estimates}: plot credible intervals and Bayes' estimate of each sample age on one 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
#' [rjags], [plot_MCMC], [SCMatrix], [plot_Ages]
#'
#' @references
#' Christen, JA (1994). Summarizing a set of radiocarbon determinations: a robust approach. Applied Statistics, 489-503.
#'
#' 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
#' data(DATA_C14,envir = environment())
#' C14Cal <- DATA_C14$C14[,1]
#' SigmaC14Cal <- DATA_C14$C14[,2]
#' Names <- DATA_C14$Names
#' nb_sample <- length(Names)
#'
#' ## Age computation of samples without stratigraphic relations
#' Age <- AgeC14_Computation(
#'  Data_C14Cal = C14Cal,
#'  Data_SigmaC14Cal = SigmaC14Cal,
#'  SampleNames = Names,
#'  Nb_sample = nb_sample,
#'  PriorAge = rep(c(20,60),nb_sample),
#'  Iter = 500,
#'  quiet = TRUE,
#'  roundingOfValue = 3)
#'
#' @md
#' @export
AgeC14_Computation <- function(Data_C14Cal,
 Data_SigmaC14Cal,
 SampleNames,
 Nb_sample,
 PriorAge = rep(c(10, 50), Nb_sample),
 SavePdf = FALSE,
 OutputFileName = c('MCMCplot', "HPD_CalC-14Curve", "summary"),
 OutputFilePath = c(""),
 SaveEstimates = FALSE,
 OutputTableName = c("DATA"),
 OutputTablePath = c(''),
 StratiConstraints = c(),
 sepSC = c(','),
 Model = c("full"),
 CalibrationCurve = c("IntCal20"),
 Iter = 50000,
 t = 5,
 n.chains = 3,
 quiet = FALSE,
 roundingOfValue = 3
){


   #--- BUG file selection
   Model_AgeC14<-0
   data(Model_AgeC14,envir = environment())

   #--- 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])

   # #--- Calibration curve
   # TableauCalib=read.csv(file=paste("inst/extdata/",CalibrationCurve,"_CalC14.csv",sep=""),sep=",",dec=".")
   # AgeBP=rev(TableauCalib[,1])
   # CalC14=rev(TableauCalib[,2])
   # SigmaCalC14=rev(TableauCalib[,3])

   #--- 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)
     }
   }

   if(Model=="full"){
    dataList = list('X'=Data_C14Cal,"sigma"=Data_SigmaC14Cal,'N'= Nb_sample,
                   "xTableauCalib"=AgeBP,"yTableauCalib"=CalC14,"zTableauCalib"=SigmaCalC14,
                   "xbound"=PriorAge,"StratiConstraints"=StratiConstraints)
   }else{
     dataList = list('X'=Data_C14Cal,"sigma"=Data_SigmaC14Cal,'N'= Nb_sample,
                     "xTableauCalib"=AgeBP,"yTableauCalib"=CalC14,
                     "xbound"=PriorAge,"StratiConstraints"=StratiConstraints)
   }

   ##set connection to model
   con <- textConnection(Model_AgeC14[[Model]])

   jags <-
     rjags::jags.model(
       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('Age','Z'),min(Iter,10000),thin=t, progress.bar = progress.bar)
   U <- summary(echantillon)

   Sample=echantillon[[1]]
   for(i in 2:n.chains){
     Sample=rbind(Sample,echantillon[[i]])
   }

   nom=c()
   for(i in 1: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, variables = "Age")

   if(SavePdf){
    dev.off()
   }

   Outlier <- SampleNames[which(U$statistics[(Nb_sample+1):(2*Nb_sample),1]<1.5)]

   ##if Outliers are detect, return warning
    if(length(Outlier) > 0){
      warning(paste("[Age14_Computation()] Outliers detected in sample:", paste(Outlier, collapse = ", ")), call. = FALSE)

    }


   ##- Gelman and Rubin test of convergency of the MCMC
   CV=gelman.diag(echantillon,multivariate=FALSE)
   cat(paste("\n\n>> MCMC Convergence of Age parameters <<\n"))
   cat("----------------------------------------------\n")
   cat(paste("Sample name ", " Bayes estimate ", " Uppers credible interval\n"))
   for(i in 1: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: The following information are only valid if the MCMC chains have converged  ***\n")
   cat("---------------------------------------------------------------------------------------------------\n\n")

   # Matrix of results
   rnames=c()
   for(i in 1:Nb_sample){
     rnames=c(rnames,paste("A_",SampleNames[i],sep=""))
   }
   R=matrix(data=NA,ncol=8,nrow=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=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("Sample name", "\t","Bayes estimate"," Credible interval: \n"))
     cat(paste(paste("A_",SampleNames[i],sep=""),"\t",(mean(Sample[,i])),'\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",(c(HPD_95[2])),"\t\t",(c(HPD_95[3])),"\n")
     cat("\t\t\t\t at level 68% \t",(c(HPD_68[2])),"\t\t",(c(HPD_68[3])),"\n")
     AgePlot95[i,] <- HPD_95
     AgePlot68[i,] <- HPD_68
     AgePlotMoy[i]=(mean(Sample[,i]))

     R[i,3]=(mean(Sample[,i]))
     R[i,c(1,5)]=(HPD_95[2:3])
     R[i,c(2,4)]=(HPD_68[2:3])
     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")

   # Representation graphique des resultats
   #        des HPD sur la courbe de calibration
   couleur=rainbow(Nb_sample)
   par(mfrow=c(1,1),las = 0,mar=c(5,5,2,2))
   xl <- c(min(PriorAge[seq(1,(2*Nb_sample-1),2)]),max(PriorAge[seq(2,(2*Nb_sample),2)]))
   plot(
      xl,
      xl,
      col = "white",
      xlab = "C-14 Age BP (ka)",
      ylab = "C-14 cal. BP (ka)",
      xaxt = "n",
      yaxt = "n",
      cex.lab = 1.8
   )
   axis(2,cex.axis=2)
   axis(1,cex.axis=2)
   polygon(c(TableauCalib[, 1], rev(TableauCalib[, 1])),
           c(TableauCalib[, 2] + 2 * TableauCalib[, 3],
             rev(TableauCalib[, 2] - 2 * TableauCalib[, 3]))/1000,
           col = "gray", border = "black")
  for (i in 1:Nb_sample) {
    lines(c(AgePlot95[i, 2:3]), rep(Data_C14Cal[i], 2)/1000, col = couleur[i],
          lwd = 4)
    lines(AgePlotMoy[i], Data_C14Cal[i]/1000, col = "black", lwd = 2,
          type = "p")
  }
   legend("topleft",SampleNames,lty=rep(1,Nb_sample),lwd=rep(2,Nb_sample),cex=1,col=couleur)

   if(SavePdf){
     dev.print(pdf,file=paste(OutputFilePath,OutputFileName[2],'.pdf',sep=""),width=8,height=10)
   }


   # Create return objecty -------------------------------------------------------------------------
   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,
      "Outlier" = Outlier,
      "Model" = Model,
      "CalibrationCurve" = CalibrationCurve,
      "PriorAge" = PriorAge,
      "StratiConstraints" = StratiConstraints
   )

   # 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)
}