Analysis of Archaeological Phases using the CRAN Package ArchaeoPhases

We propose new statistical tools to analyse and to estimate the characteristics of time periods based on the posterior distribution of the associated collection of dates. These tools are implemented in the R package ArchaeoPhases . The required inputs are simulated samples from the posterior distribution of a collection of dates. Such Markov Chains Monte Carlo samples are provided, for instance, by BCal , ChronoModel or OxCal , freely available software applications build for the chronological modeling of estimated dates. We give a practical introduction to the package ArchaeoPhases using published data and comment the statistical results.


Introduction
Statistical modelling within the Bayesian framework is widely used for constructing chronologies based on relative and absolute chronological information. This consists in the estimation of a collection of dates where • the observations are measurements (along with their error) coming from possibly different dating methods (radiocarbon, archaeomagnetism, thermoluminescence...), • the prior information about these dates is induced by archaeological, geological, historical, environmental or any other considerations.
Bayesian inference is a probabilistic estimation method, in the sense that the information on the estimated parameters is given by a probability distribution called the posterior distribution. Usually the analytic expression of the posterior distribution is not easily computable. However, inferences may be drawn using a sample generated from the posterior distribution. More precisely, it is possible to compute a Markov chain whose target distribution is the posterior distribution. Indeed, Markov chain Monte Carlo (MCMC) algorithms provide a way of drawing samples from the joint posterior distribution in high-dimensional Bayesian models.
Software applications available to estimate chronologies include: BCal (see Buck, Boden, Christen, James, and Sonnenwald (2018); Buck, Christen, and James (1999)), Chronomodel (see ; Vibet, Philippe, Lanos, and Dufresne (2016); Philippe (2017, 2018)) and OxCal (see Bronk Ramsey (2016Ramsey ( , 2009), but also JAGS via the package rjags or via the package ArchaeoChron (see Philippe and Vibet (2017)), and STAN via the package rstan. These applications return the marginal posterior distribution of all the estimated dates with elementary statistics (mean, standard deviation, HPD region,. . . ). It is also possible to export the output of MCMC simulations from the joint posterior distribution of the dates. These simulated samples can be imported into R. Hence, we can explore the joint distribution of a collection of dates using the ArchaeoPhases package.
The ArchaeoPhases package provides new statistical tools for the analysis of chronologies based on the joint distribution of estimated dates. Two first tools explore, for a collection of dates, the rate of occurrence of these dates (TempoPlot() et TempoActivityPlot()). Three other tools summarise periods of time, or what we call phases, defined by some archaeological, paleo-environmental, geographical... criterion. A phase is associated with a collection of calendar dates estimated from the chronological model.
In this article, we will explain and illustrate the use of the following new statistical tools: • Tempo plot introduced by Thomas Dye (see Dye (2016)), it is a statistical graphic designed for the archaeological study of rhythms of the long term that embodies a theory of archaeological evidence for the occurrence of events" and its associated Activity plot.
• Time range interval : a time interval that characterizes the period of time during which a phase happened.
• Transition range interval : a time interval that characterizes the period of time during which a transition between two successive phases happened.
• Testing procedure to check the presence of a gap between two successive phases. A gap range interval is estimated, if we accept its existence, and characterizes the hiatus between two successive phases.
The CRAN package ArchaeoPhases presents a list of functions for calculating these statistics. The inputs of these R functions are samples simulated from the posterior distribution of dates (e.g. the Markov Chains simulated by BCal, Chronomodel or OxCal). For non-R users, a web application has been developed in order to take advantages of these functions without having to know R. This application is freely available on https://archaeology-bayesian-modelling. shinyapps.io/ArchaeoPhases and within the ArchaeoPhases package (app_ArchaeoPhases()).
In Section 2, we describe the statistical aspects of our new tools developed for the analysis of phases (or their associated collection of dates) in a Bayesian framework. Section 3 provides a short description of the data published by Bosch, Mannino, Prendergast, O'Connell, Demarchi, Taylor, Niven, van der Plicht, and Hublin (2015) used to illustrate the use of Ar-chaeoPhases and explanations about the problematic. Then we construct a Bayesian model with Chronomodel in order to estimate four Palaeolithic periods based the information of the stratigraphic sequence and on radiocarbon dates found at Ksâr 'Akil (Lebanon). In Section 4, we give a practical introduction to the first steps with ArchaeoPhases (version 1.3). Then in Section 5, we analyse the MCMC output of the modelling done with Chronomodel and estimate the four Palaeolithic periods thanks to the information found at Ksâr 'Akil.

Statistical aspects
We define a chronology as a collection of calendar dates τ 1 , . . . , τ n . Usually, the dates τ 1 , . . . , τ n are not observed, they are estimated from a Bayesian chronological model that includes relative and absolute chronological information (e.g BCal, ChronoModel, OxCal). The relative information relates to the temporal ordering of dates whereas the absolute information usually arises from possibly different dating methods.

Notation:
We denote by M the set of measurements coming from dating methods. We assume that MCMC samples from the joint posterior distribution p(τ 1 , . . . , τ n | M) of all the dates τ 1 , . . . , τ n are available.
Two types of questions may arise when considering a chronology of dates : the rate of occurrence of these dates and the periods of time during which these collections of dates happened. We propose new statistical tools build as post-processing steps of the MCMC sample of the chronology in order to give answers to these questions.

Rate of occurrence of dates
Dye (2016) defined a graphical tool, called tempo plot, to evaluate the rate of occurrence of the dates. We propose a Bayesian interpretation of this estimate. The quantity of interest cannot be viewed as a counting process because the date of the events (τ 1 , . . . , τ n ) are not observed. For each date t, the aim is to estimate the number of dates N (t) which occurs before the date t, we have As N (t) is a function of the parameters (τ 1 , . . . , τ n ), we can easily estimate N (t) from the joint posterior distribution of (τ 1 , . . . , τ n ). The Bayes estimate of N (t) (under quadratic loss) is the posterior mean of N (t) i.e., Using the output of the MCMC algorithm (τ ..,N M CM C , we can approximate the Bayes estimateN (t) by taking where F (k) is the empirical cumulative distribution of the sample (τ The N M CM C functions nF (k) (.) provide a sample from the posterior distribution of N . Therefore we can easily build a credible confidence region for the function N . Indeed for each t we take the smaller posterior interval approximated from the sample nF (k) (t), k = 1, . . . , N M CM C . An alternative is to use Gaussian approximation to get the confidence intervals. To visualize the rhythms, an alternative to the Tempo Plot consists in representing the first derivative ofN defined in (1). This graphical tool is called Activity Plot.

Periods of time
The question about a chronology is not only a matter of estimation of these calendar dates τ 1 , . . . , τ n but also the characterization of historical, or geological, or archaeological... periods of time, or what we call phases. A phase is a collection of dates estimated from the chronology.
We define a phase (denoted P ) by a collection of dates gathered on some criterion, τ i for i ∈ I where I is a subset of {1, . . . , n} and n the total number of calendar dates.
In this section, we present the classical way of summarising a periods of time, or a phase, and we introduce our new concepts whose aim are to further interpret the MCMC sample from the joint distribution of the chronology of dates.

Classical summary
A phase is commonly summarized by two parameters: the start and the end dates of the corresponding period of time. Two different approaches are considered in order to estimate these two dates.
1. The start and end dates are modelled, and so additional parameters are incorporated in the model. This approach is implemented in OxCal application where the parameters are denoted t a and t b . This requires to do prior assumption on the distribution of the dates belonging to the phase on the estimated period [t a , t b ], such as a uniform distribution on [t a , t b ]. See Naylor and Smith (1988); Buck, Litton, and Smith (1992); Christen (1994); Buck, Litton, and Cavanagh (1996); Nicholls and Jones (2002) for different choice of prior distributions.
2. The start (resp. the end) date is estimated by the date of the minimum (resp. the maximum) of the associated collection of dates. This approach does not model any additional parameter to define the phase. It is only a reflection of the time period covering the collection of dates, a post-processing step based on the MCMC sample of the collection of dates. These statistics are implemented in ChronoModel and OxCal applications.
With both approaches, a phase is characterised by two dates (start and end) for which a credible interval may be given. However, these two estimated parameters do not provide a localisation of the phase in time. For instance, one could use the Bayes estimate of the start and end dates of a phase in order to build an interval. However, no one knows what is the probability that all the dates of the collection belong to that interval.
Our new tools are build in order to characterise a phase by a time interval with fixed probability that the dates of the collection belong to it.

Time Range interval
We propose to characterize a phase by a time period that contains all the dates of the associated collection with a given probability. The so-called time range interval gives an idea of the start, the end and the duration of the phase. It locates a phase in time.
The 100(1 − γ)% time range is the shortest interval (a, b) such that This means that the time interval [a, b] contains all the collection of dates {τ i i ∈ I} with a fixed posterior probability 1 − γ (e.g., 95 %, 68 % . . . ).
We denote α = start(τ i , ∀i ∈ I) and β = end(τ i , ∀i ∈ I). Equation (2) can also be rewritten in the form Therefore the construction of the time range interval depends only on the joint posterior distribution of (α, β).
The intervals [a, b] satisfying Equation (3) are those that satisfy: is the quantile function of the posterior distribution of α, is the quantile function of the conditional distribution of β given α ≥ a( ) and M.

So we have
The last equality comes from the definition of a( ) and b( ). Therefore the shortest interval is obtained by taking In practice, the values of a( ) and b( ) are estimated by the empirical quantiles calculated on the MCMC outputs. Using the output of the MCMC algorithm (τ we produce (α (k) , β (k) ) t=1,...,N M CM C the MCMC sample corresponding to the posterior distribution of (α, β),

Transition interval between two successive phases
We propose to characterize the transition of time between two successive phases, for phases in temporal order relationship, by a transition interval. This tool gives an estimation the period of time during which the older phase has not yet ended and the younger phase has just begun. For phases in temporal order relationship, the transition range interval between two successive phases is the shortest interval that covers the end date of the older phase and the start date of the younger phase. From a computational point of view this is equivalent to the time range interval calculated between the end of the older phase and the start of the younger phase.
Definition 2. Consider a succession of two phases, P 1 and P 2 . Assume P 1 is older than P 2 . We denote by α i (respectively β i ) the start date (respectively the end date) of the collection of dates The construction of this interval is obtained the same way as the time range interval by replacing (α, β) by (β 1 , α 2 ).

Testing procedure for gap between two successive phases
Phases in temporal order relationship may be separated in time by a gap, or a hiatus. We propose a testing procedure in order to check whether a gap of time exists between two successive phases with fixed probability. For phases in temporal order relationship, if a gap exists, it is the longest interval in between the end date of the older phase and the start date of the youngest one with fixed posterior probability.
Definition 3. Consider a succession of two phases, P 1 and P 2 . Assume P 1 is older than P 2 . We denote by α i (respectively β i ) the start date (respectively the end date) of the collection of dates P i (i = 1, 2).
The first step consists in the construction of all the couples a( ), b( ) such that is the quantile function of the posterior distribution of β 1 , is the quantile function of the conditional distribution of α 2 given β 1 < a( ) and M. Indeed is empty, then we conclude that no gap exists between these successive phases, with probability 1 − γ. Note that the probability 1 − γ cannot be interpreted as a confidence level. But it means that we cannot construct an interval satisfying (5) with posterior probabilty 1 − γ. Otherwise, the gap interval is the

Bayesian modeling for the stratigraphy of Ksâr 'Akil
We illustrate the use of the ArchaeoPhases package and mainly the functions described in Section 2 with the construction and the estimation of a Palaeolithic chronology based on the stratigraphic information and measurements made on the site of Ksâr 'Akil (Lebanon).
We estimate the Palaeolithic chronology thanks to a Bayesian model implemented in ChronoModel (see Vibet et al. (2016)). This software application has a user-friendly graphical interface easily manipulated in order to define the Bayesian hierarchical model developed in Philippe (2017, 2018). As this application is not connected with R, we should first construct the model with ChronoModel, then extract from this application the MCMC sample of the joint posterior distribution, and finally import the MCMC sample into R application and apply the ArchaeoPhases functions.

The stratigraphy of Ksâr 'Akil
At Ksâr 'Akil (Lebanon), a deep Palaeolithic stratigraphic sequence was investigated in order to established the chronology of the site (see Bosch et al. (2015)). This stratigraphic sequence is constituted of a succession of thirty-six layers. Based on the lithic assemblages found throughout the sequence, twenty-five layers were associated with four Palaeolithic deposits : Initial Upper Palaeolithic (IUP), from layer XXV (the bottom of the sequence) to layer XXI; Ahmarian, from layer XX to layer XIV; Upper Palaeolithic (UP), from layer XV to layer VI; and Epi-Palaeolithic (EPI), from layer V to layer I. The aim of this modelling is to establish the chronology of these four successive Palaeolithic periods (IUP, Ahmarian, UP and EPI) thanks to the chronology of the stratigraphy. Sixteen consumed shellfish were found throughout the stratigraphy and dated by AMS radiocarbon technique in order to estimate the date of their death. Dating the death of these consumed shellfish, whose shells were directly discarded after consumption, gives an indication of human activity on the site. Hence it gives also an indication of the different Palaeolithic periods during associated with this human activity. Thus, the chronology of the site was established on these evidences of human activity and on the stratigraphic information.

Modelling with Chronomodel
Based on the information given by Bosch et al. (2015), we propose a Bayesian model for the stratigraphy of Ksâr 'Akil implemented in ChronoModel (version 1.5).
In this model, the study period is set from -50 000 to -10 000 (dates are in calendar years Before Common Era) out of historical information on the Upper Palaeolithic.
As the aim of this modelling is to date the human occupations of the site, sixteen molluskan species collected for consumption and whose shells were directly discarded after consumption were dated by radiocarbon. This technique estimates the date of the death of these mollusks. Hence, dated events (death of mollusk) and target events (evidence of human activity) are nearly identical (see Dean (1978)). We use these dated events in order to date human activity in the layer in which these food remains were found. Figure 1 presents the modelling realised with ChronoModel application. Dated events (death of mollusk) are called according to the name given by the laboratory (for instance : GrA-53081). Each dated event is a proxy of the human occupation in the layer in which the mollusk has been found. Target events (evidence of human activity) are named according to the associated layer. In the best cases, several mollusks were found in the same layer, and so all of these dates give an estimation of the period of time corresponding to human activity in that layer. Stratigraphic constraints, that give a temporal order to the occupation of the different layers, are also included in the model, they are symbolised by black arrows in Figure 1. Then, layers and the 16 related events, are gathered in Palaeolithic periods in order to give an estimation of these Palaeolithic periods. In ChronoModel, four phases were constructed based on the evidences of human activity found in several layers. Figure 1 shows the stratigraphy of the layers in which sellfish were found, and the four Palaeolithic deposits.

Extracting the MCMC sample of the joint distribution
For this Bayesian modelling, 3 Markov chains are run in parallel. For each chain, 1 000 iterations are generated during the Burn-in period, 20 batches of 500 iterations are used in the Adapt period, 100 000 iterations are drawn in the Acquire period by only 1 out of 10 are kept in order to break the correlation structure. Now, the Markov chains simulated by ChronoModel have to be extracted and then imported into R in order to use the ArchaeoPhases package. To do that, use on the right hand side of the window in the Results tab (see Figure 2). Several CSV files are created: a file called "events.csv" containing the MCMC samples of all events (target events), a file called "phases.csv" containing the date of the minimum (called alpha) and the date of the maximum (called beta) of each phase, if at least one phase is modelled, and a file per phase (if any) containing the MCMC samples of the parameters of the phase (the minimum and the maximum date) and all events included in it. Note that the minimum (resp. maximum) of a group of dates is an estimation of the start (resp. end) of the group of dates. Figure 1: ChronoModel interface showing the modeling step of the chronology. Rectangular nodes represent the evidences of human activity (target events) associated with the radiocarbon date of the death of a mollusk (dated events). An arrow between two rectangular nodes symbolises a time order relationship. Gray lines symbolise the four Palaeolithic deposits and the layers associated with these depositis.
All those files may be imported into R in order to use the ArchaeoPhases package. Then the first step is to load the library: R> library("ArchaeoPhases")

Importing data into R software
To import a data file into R, you may use ImportCSV() function. For CSV files extracted from ChronoModel software, there is no need to specify any other parameters than the name of the file (and the path to it). Otherwise, you may change the specification after sep= and dec=. The parameter comment.char= is used to define how comments are written in the file to be imported. Comments of all csv files generated by ChronoModel are specified by '#'. For example, R> KADatesChronoModel = ImportCSV("events.csv") In the package ArchaeoPhases, a dataset containing the MCMC samples of the posterior distribution of the 16 events identified in Ksâr 'Akil , and generated by ChronoModel application as explained in the previous section, is already included.

R> data(KADatesChronoModel)
This dataset contains 3 Markov chains run in parallel. The first column corresponds to the iteration number, the following ones correspond to the posterior distribution of the 16 target events.

Diagnostic tools
Before any analysis, the convergence of the MCMC sample should always be checked (for more details on the diagnostic of Markov chain, see Robert and Casella (2009)). In this section, we show how this can be checked using R and the coda package. Note that ChronoModel application has also its own diagnostic tools.
In order to use any function from the coda package, data has to be transformed in a mcmc.list using the function coda.mcmc(). This function takes two parameters: the name of the dataset and the number of parallel chains generated by the application. R> Dates_mcmc = coda.mcmc(data = KADatesChronoModel, numberChains = 3) Now, we can trace the history plot of each parameter and see whether the chains have reach their equilibrium before the Acquire period. The column 1 containing the iteration number should be withdrawn.

R> plot(Dates_mcmc[,-1])
A part of the results is given in Figure 3. We can see that the 3 chains give similar results. For each parameter, traces visit the same time space of values, traces are correctly overlaid, and they seem to have a stationary behaviour. Indeed, within each subinterval, the same time space is visited. Hence, we can assume that the 3 chains have properly reach the same equilibrium, and finally that the overall MCMC sample (the concatenation of the 3 chains) has converged.  The values of the Gelman-Rubin criterion are equal to 1 meaning that the different Markov chains run in parallel have converged.
In summary, the Markov chains have properly reached their equilibrium. We can concatenate the 3 chains and go further in the analysis.

Examining a series of dates
First, we analyse the individual marginal posterior density of the 16 events that correspond to evidence of human activity on the site. In a Bayesian framework, an estimation of a date can by summarised by its confidence interval given at a 95% probability thanks to the marginal posterior density of the related date. This gives an interval of time that has 95% to contain the true date of interest. This interval is not assumed to have equal tail probabilities on the left and right but it is optimised in order to be the shortest. However, this statistical tool is not adapted if the posterior density is bimodal, Highest Posterior Density (HPD) regions would be preferred. As expected, the evidence of human activity in Layer XXII, the deepest dated layer, is associated with the oldest 95% confidence interval. Layer XVII and layer XVI are associated with four dates as four shellfish were found in them. The 95% confidence interval of the four evidences of human activity in Layer XVII are more or less similar.
A Highest Posterior Density (HPD) region (see Hyndman (2018) Figure 4. R> MultiDatesPlot(KADatesChronoModel, position = c(2:17), level = 0.95, intervals="CI", title=" 95% CI of the dates", order = "increasing") R> MultiDatesPlot(KADatesChronoModel, position = c(2:17), level = 0.95, intervals="HPD", title=" 95% HPD regions of the dates", order = "increasing") We can see that dates are properly sorted according to the stratigraphic order. Confidence intervals are close to each other from layer to layer, except between layer XI and layer VI. Indeed no shellfish were found in the layers VII to X leading to a gap of time in this series of dates.
We can examine the rate of occurrence of the human activity throughout the stratigraphy.
To do that two functions may be used : TempoPlot and TempoActivityPlot. at Ksâr 'Akil. From these graphs, we can see that the highest part of the human activity is dated between -43 000 to -35 000 and two evidences are younger, at about -31 000 and -28 000. The activity plot shows that there is an increasing amount of activity after -43 000 with a peak at -40 000, then the activity decreases until -35 000 where it is null. Then two other activities appear later in time. There seem to be a gap of time between the first 14 dates and the two remaining.

First conclusions
From these three last graphs, we can see that there is a continuity in the occurrences of human activity found at Ksâr 'Akil between -43 000 and -35 000. This correspond to the activity found in layers XXII to XI. Then there is a hiatus until -29 000. Then two other occurrences of human activity appears. These occurrences were found in layer VI and layer V. The hiatus highlighted there may correspond to the hiatus linked to the time elapsed between the deposit of layer XI and of layer VI. In this context, the Tempo Plot and the activity plot reflect more the sampling of the evidences of human activity, that is the sampling of a peculiar kind of consumed shellfish throughout the stratigraphy. In addition, only 16 of these shellfish could be dated by radiocarbon. This is probably due to a random process. In some layers, several shells were found. and dated. Hence human activity related to these layers were more precisely dated than in layer with only one shell. The gap of time is probably due to the absence of dated shells rather than to In all, human activity found at Ksâr 'Akil is dated at a period of time in coherence with other results, such as Bosch et al. (2015); Douka, Bergman, Hedges, Wesselingh, and Higham (2013). However, the analysis can be further refined by the inspection of the Palaeolithic periods of time, called phases hereafter, and by the use of our new statistical tools.

Palaeolithic chronology
The lithic assemblages found throughout the stratigraphy of Ksâr 'Akil allowed us to associate the first XXV layers to four Palaeolithic periods: Initial Upper Palaeolithic (IUP), from layer XXV (the bottom of the sequence) to layer XXI; Ahmarian, from layer XX to layer XIV; Upper Palaeolithic (UP), from layer XV to layer VI; and Epi-Palaeolithic (EPI), from layer V to layer I. he aim of this modelling is to establish the chronology of these four successive Palaeolithic periods thanks to the chronology of the stratigraphy.
In ChronoModel, we defined 4 phases in order to estimate the characteristics of these four Palaeolithic periods. As said in Section 3, several CSV files are created by ChronoModel and in particular, a file called "phases.csv" containing an estimation of the minimum date (called alpha) and maximum date (called beta) of each phase, as long as at least one phase is modelled. In ChronoModel the start date (resp. end date) of a phase is estimated by the minimum date (resp. the maximum date) of the collection of dates associated to this phase.
In the ArchaeoPhases package, a dataset containing the MCMC samples of the start date (called alpha) and the end date (called beta) of each Palaeolithic phase, and generated by ChronoModel application, is already included.

R> data(KAPhasesChronoModel) R> attach(KAPhasesChronoModel)
Let us first analyse a unique period of time, the period called Ahmarian, then we will examine the succession of the four Palaeolithic periods.

Examining the Ahmarian period
Let us analyse the Ahmarian period that gathers 11 dates of the 16 from the chronology.

Classical summary
A period of time is classically summarised by its start and end dates. From these two dates, we can estimate the duration of the period. The function PhaseStatistics() gives the classical summary statistics of the characteristics of a phase: mainly mean, median, credible interval and HPD region. This function takes a minimum of two parameters : a MCMC sample of the marginal posterior distribution of the start date and a MCMC sample of the marginal posterior distribution of the end date. R> PhaseStatistics(Ahmarian.alpha, Ahmarian.beta, level = 0.95) The output of the R console is then :

Time Range interval
We can estimate locate in time the Ahmarian period by its time range interval. This gives an interval of time that resumes the start, the end and the duration of the period at a given confidence level.
The endpoints of the time range of Ahmarian are close to the inferior endpoint of the 95% credible interval of the start date and the superior endpoint of the 95% credible interval of the end date. However, only the time range gives an interval that covers the collection of dates of the phase with a 95% probability.

Graphical illustration
Function PhasePlot() may be used to draw a plot of the marginal posterior density of the start and end dates of a phase on the same graph. In addition, the time interval of the phase at the desired level (default = 95 %) is also presented by a solid line above the curves. The result is shown on Figure 6 and displays the characteristics of a phase.

Characteristics of the Palaeolithic period called Ahmarian
Calendar

Examining the four Palaeolithic periods
Now, let us analyse the succession of the four Palaeolithic periods: IUP, Ahmarian, UP and EPI.

Time range intervals
The function MultiPhaseTimeRange() gives the time range intervals of several phases. It takes at least two parameters, data the name of the dataset, and position_minimum is a vector of positions that identify the start dates of the corresponding phases in the data matrix. By default, the vector of positions that identify the end dates is assumed to be position_maximum = position_minimum + 1. Let's estimate the time range of phases IUP, Ahmarian, UP and EPI, whose start date outputs are respectively in column number 8, 6, 4 and 2 in the data matrix. Phase "EPI" whose start date output is in column 2 has a time range = [-29 071, -27 102] at a confidence level of 95% and the phase "IUP" whose start date output is in column 8 has a time range = [-43 217, -41 106].

Transition intervals
Now, let us estimate the transition intervals between two successive phases. R> MultiPhasesTransition (KAPhasesChronoModel,position_minimum=c(8,6,4,2)) For this functions the order of the succession is important. The vector of positions of the start date (referring to the column number in the data matrix) should be ordered from the oldest phase to the youngest one.

Gap intervals
We can test if a the gap between these successive phases exists. This may be done using the following code : R> MultiPhasesGap (KAPhasesChronoModel,position_minimum=c(8,6,4,2)) Again, for this function, the order of the succession is important. The vector of the start position (referring to the column number in the data matrix) should be ordered from the oldest phase to the youngest one. With probability 95% there is no gap between the succession of phases IUP, Ahmarian and UP. Indeed we cannot construct an interval satisfying (5) with posterior probability 95% . On the other hand there exists a gap of 203 years between phase UP and phase EPI. This gap interval, with posterior probabability 95% , starts at -29 180 and ends at -28 977.

Graphical illustration
The function MultiSuccessionPlot() draws the characteristics of a succession of groups. Figure 7 presents the resulting plot.
• The density curves represent the densities of the minimum (solid line) or the maximum (dashed line) of each phase. When a group is defined by one only date, the minimum equals the maximum date of the group. Hence, in that case, only one curve is presented (see the phase IUP in Figure 7).
• Segments above the density curves correspond to time range of the phase.
• Two-coloured segments correspond to transition interval and to the gap range. When there is no gap between the succession of phases, the two-coloured segment is replaced by a a cross (see Figure 7 between phases IUP and Ahmarian, and Ahmarian and UP). Two-coloured segments correspond to transition interval or to the gap range. As there are no gaps between phases IUP and Ahmarian, and Ahmarian and UP, a cross is drawn instead.

Summary
We can summary the succession of the four Palaeolithic periods estimated from the site of Ksâr 'Akil, using the data published by Bosch et al. (2015) and the modeling done with ChronoModel described in section 3.2, by the following tables : Results from Table 1 are similar as those given by Bosch et al. (2015). We can see that for phases containing only one events, such as EPI and IUP, the 95% credible interval of the start  (or end) date is equal (up to numeric approximations) to the 95% time range interval. For phases containing several events such as Ahmarian and UP the time range interval provide rigourous alternative to the usual approach consisting to take the interval whose endpoints are the inferior value of the credible interval of the start date and the superior value of the credible interval of the end date. However, only the time range interval gives a 95% interval for the period covering all the dates included in the phase. Table 2 gives also the endpoints of the transition range intervals and the gap range intervals if existing. Again, the 95% transition intervals are similar to the interval build with the inferior value of the 95% credible interval of the end date of the older phase and the superior value of the 95% credible interval of the start date of the younger phase. However, the only interval that gives a 95% probability to include all the dates is the transition interval. The gap range intervals are similar to an interval build with the superior value of the 95% credible interval of the end date of the older phase and the inferior value of the 95% credible interval of the start date of the younger phase.
Our new statistical tools aim at summarizing periods of time by intervals that cover all the dates with a fixed probability. Instead of summarizing a phase by two credible intervals, one for the start date, one for the end date, the time range interval is a unique interval given with fixed probability.
In conclusion, we can summarize the succession of the Paleolithic phases by saying that from the data found at from the site of Ksâr 'Akil, using the data published by Bosch et al. (2015) and the modeling done with ChronoModel described in section 3.2, with a 95% probability, Phase IUP is estimated from -43 217 BCE to -41 106 BCE, Phase Ahmarian is estimated from -42 189 to -37 461, Phase UP is estimated from -38 559 to -29 335 and Phase EPI is estimated from -29 071 to -27 102 (years in BCE). The transition periods between successive phases are estimated with a 95% probability.

Conclusion
The ArchaeoPhases package provides new statistical tools that help interpret the MCMC sample generated from the posterior distribution of a chronology of dates. Among these tools, two of them explore the rate of occurrence of dates (TempoPlot() et TempoActivityPlot()). Three other tools summarise periods of time, or what we call phases, defined by some archaeological, paleo-environmental, geographical... criterion. A phase is associated with a collection of dates estimated from the chronology. PhaseTimeRange() (or MultiPhaseTimeRange()) gives an interval of time that characterizes the period corresponding to that phase with a fixed probability. PhasesTransition() (or MultiPhasesTransition()) gives an interval of time that characterizes the transition between two successive phases with a fixed probability. PhasesGap() (or MultiPhasesGap()) tests the existence of a hiatus between two successive phases, and if so, gives an interval of time that characterizes the hiatus between these two successive phases with a fixed probability. See Guérin, Antoine, Schmidt, Goval, Hérisson, Jamet, Reyss, Shao, Philippe, Vibet, and Bahain (2017) and  for an application of these tools. We chose to illustrate the use of these statistical tools, in comparison with the classical ones, by the analysis of four Palaeolithic periods and the stratigraphy of Ksâr 'Akil. The Bayesian modelling made with ChronoModel, a software application adapted for the construction of chronologies of estimated dates, resulted in a chronology similar to the one given by Bosch et al. (2015). Our new tools add complementary information to classical ones. The time range interval summarises a phase by a unique interval with a fixed probability. Transition range and gap range intervals result in a unique interval that defines a transition between two phases or a gap between two phases.
The ArchaeoPhases package has been applied to archaeology (see ; Duvivier and Florent (2018); Philippe, Vibet, and Bahain (2018)) but it could also be used in other fields such as paleo-environnement (see Guérin et al. (2017)), forensic physics and so on as long as dates are estimated.