We propose a robust event date model to estimate the date of a target event by a combination of individual dates obtained from archaeological artifacts assumed to be contemporaneous. These dates are affected by errors of different types: laboratory and calibration curve errors, irreducible errors related to contaminations, and taphonomic disturbances, hence the possible presence of outliers. Modeling based on a hierarchical Bayesian statistical approach provides a simple way to automatically penalize outlying data without having to remove them from the dataset. Prior information on individual irreducible errors is introduced using a uniform shrinkage density with minimal assumptions about Bayesian parameters. We show that the event date model is more robust than models implemented in BCal or OxCal, although it generally yields less precise credibility intervals. The model is extended in the case of stratigraphic sequences that involve several events with temporal order constraints (relative dating), or with duration, hiatus constraints. Calculations are based on Markov chain Monte Carlo (MCMC) numerical techniques and can be performed using ChronoModel software which is freeware, open source and cross-platform. Features of the software are presented in Vibet
Bayesian chronological modeling appears as an important issue in archeology and palaeo-environ mental science. This methodology has been developed since the 1990s (Bayliss, 2009, 2015) and is now the method of choice for the interpretation of radiocarbon dates. Most applications are undertaken using the flexible software packages, BCal (Buck
All of these models provide an estimation of a chronology of dated events (DE), the estimated dates correspond to the dates of events that are actually dated by any chronometric technique. However calibrated radiocarbon dates and date estimates from other chronometric dating methods such as thermoluminescence, archaeomagnetism, and dendrochronology. can be combined with various prior archeological information to produce a combined chronology that should be more reliable than its individual components.
We propose a new chronological Bayesian model based on the concept of a target event. This is related to the concepts of the dated event and the target event proposed by Dean (1978). The target event (TE) is the event to which the date is to be applied by the chronometrician. Usually, the target events are not directly related to the dated events. According to Dean (1978), the dated event is contemporaneous with the target event when there is convergence (the two dates are coeval: DE = TE) and when the date is relevant to the target event. Relevance refers to the degree to which the date is applicable to the TE : it must be demonstrated or argued based on archaeological or other evidence. It is generally not easy or possible to assure the relevance of the dated event to the target event date of interest; therefore, it is recommended to get many dates of the dated event, and if possible from different dating techniques.
These dates can be outliers without having a means to determine if they are dating anomalies. In other words, we do not have any convincing archaeological arguments for rejecting them before modeling. This motivates the development of a robust statistical model for combining theses dates in such a way that it is very little sensitive to outliers (Lanos and Philippe, 2017).
The target event model is a statistical model introduced in (Lanos and Philippe, 2017) for estimating the date of an event called a target event. This model allows to combine in a robust way the dates of artifacts, which are assumed to be contemporary to this target event. To validate the robustness of this model to outliers, we provide a comparison of our model with the
The target event model is then integrated to a global model for constructing chronologies of target events. Prior information is brought upon the dates of the target events. We show how the prior archeological information based on relative dating between the target events (in a stratigraphic sequence for instance) or based on duration, hiatus or terminus post quem (TPQ) or terminus ante quem (TAQ), assessment can be included in the model. The simulations and applications illustrates the improvement brought in term of robustness. The proposed model is implemented in ChronoModel software, whose features are described in Vibet
Many issues in archaeology raise the problem of phasing or how to characterize the beginning, the end and the duration of a given period. This question can be viewed as a post processing of the chronological model as in Philippe and Vibet (2017a), Guérin
The paper is organized as: In Section 2 we recall the construction of the target event model and we provide a theoretical comparison with the
We propose a statistical Bayesian approach for combining dates, which is based on a robust statistical model. This combination of dates is aiming to date a target event defined on the basis of archaeological/historical arguments.
We propose to use a hierarchical Bayesian model to estimate the date
The event dates
At this step, the random effect model can be written as:
where (
We assume that the experimental error provided by the laboratory are independent. Dependence structure could be added between the (
If we assume that all the measurements can be calibrated with a common calibration curve (i.e.,
The target event date
Consequently, we model the over-dispersion by an individual error
where (
Finally, the joint distribution of the probabilistic model can be written according to a Bayesian hierarchical structure:
where the conditional distributions that appear in the decomposition are given by:
The parameter of interest
This interval, called a study period, is fixed by the user based on historical or archeological evidence.
This appears as an important a priori temporal information. Note that we do not assume that the same information is imparted to the dates
The uniform shrinkage distribution for the variance
where 1
The motivation for this choice of prior is described in detail in (Lanos and Philippe, 2017). Parameter
An individual calibration step is done for each measurement
using the same notation as (
(a) Sample from the posterior distribution of
where
(b) Approximate the posterior variance var(
Take as shrinkage parameter
The calibration curve
We suggest to extend the calibration curve by an arbitrary constant value with a very large variance with respect to the known reference curve: for instance
and
In the case of TL/OSL or Gauss measurements, there is no need for an extension because
This statistical approach does not model the outliers (i.e., we do not estimate the posterior probability that a date is an outlier). Outlier modeling (Section 2.2) can provide more accurate results, but it often requires two (or more) estimations of the model: the outliers are identified after a first estimation and thus discarded from the dataset. Then the final model is estimated again from the new dataset (however, the question remains to know when to stop, because new outliers can appear during the subsequent runs). The event model which is based on the choice of robustness, avoids this two-step procedure.
Several outlier models have been implemented since the 90’s, in the framework of radiocarbon dating.
So we can distinguish between mainly three outlier modelings:
Outlier model with respect to the measurement parameter (Christen, 1994; Buck
Outlier model with respect to the laboratory variance parameter (Christen and Pérez, 2009).
Outlier model with respect to the time parameter. This corresponds to the
The
The model with random effect becomes :
where
(
the prior on
When parameter
where the conditional distributions that appear in the decomposition are given by:
This modeling can be compared to the event date model if we consider several dates nested in the Oxcal function Combine with
If we consider linear calibration curves with constant errors, by setting:
and knowing that
it is possible to analytically integrate (
where the conditional distribution of
Note that the outlier model operates well as a model averaging thanks to the structure of mixture distribution in (
However, posterior density of event date
where
Figure 1 shows a a graphical representation of the densities defined in (
We now consider a group of dates
The concept of phase is currently implemented in chronological modeling by using a specific parametrization as discussed in Section 4. Hereafter, we extend the target event date model in a simple way to the case of a group of dates related to each other by order or duration relationships. Hence we propose to estimate the beginning, the end and the duration of a phase directly from the group of target dates, without adding any supplementary parametrization.
We consider a collection of
where
When there is no supplementary information, we assumed that the dates
Figure 2 describes the overall Bayesian model using a directed acyclic graph. Such a graph describes the dependencies in the joint distribution of the probabilistic model. Each random variable of the model (that is an observation or a parameter) appears as a node in the graph. Any node is conditionally independent of its non-descendants, given its parents. The circles correspond to all random variables of the model. With the color of the circles, we distinguish between observations (red), parameters (blue) and exogenous variables (green).
We often do not know how a group of event dates can a priori be distributed in a phase. It means that, in our approach, a phase does not respond to a statistical model. From a chronological point of view, all the information is carried by the target event dates themselves. Consequently, we estimate the beginning and the end of a phase as a measurable application of the parameters
In the same way, the end of a phase is estimated by the maximum
By the plug-in principle, we estimate the duration of the phase by
Considering two phases
where
These parameters are estimated knowing the data from the target event dates which are available in the phase of interest. In particular, a good estimation of the beginning (resp. the end) of a phase requires that the archaeologist has sampled artifacts belonging to target events which are very near to this beginning (resp. this end).
These estimates are valid whatever the prior on the event dates are. Precision on the estimation of these parameters can be gained if it is possible to add supplementary information on the target event dates within the study period
Different prior on dates
Relative dating based on stratigraphy as defined in Harris (1989) and Desachy (2005, 2008), can imply antero-posteriority relationships between target dates
We can have some prior information about the maximal duration of a group of target events.
We can also have some prior information about the minimal temporal hiatus between two groups of target events.
It has been demonstrated that these types of supplementary prior information can significantly improve chronometric dates (see initial works of Naylor and Smith (1988); Buck
Target events or phases (groups) of target events can have to check order relationships. This order can be defined in different ways: by the stratigraphic relationship (physical relationship observed in the field) or by stylistic, technical, architectural,
If we consider a stratigraphic sequence composed of target events, the prior on vector
with
We can also consider two groups of target events,
or
or
The
Estimating the date of a target event needs to incorporate several dates (Ed), otherwise the event date modeling will not yield a better posterior information. However, it is possible to nest only one date per target event provided that the group of events is constrained by temporal order.
It is not rare to encounter dating results which contradict the stratigraphic order. This situation of stratigraphic inversion often occurs when some artifact movements are provoked for example by bioturbations or establishment of backfill soils. The target event date model makes it possible to manage such situations thanks to the individual variances
Prior information can be included on the duration of a phase, that is on a group of target events. We can impose a maximal duration
The
Prior information about a hiatus
The event dates will have to check total or partial stratigraphic constraints and also to satisfy the inequality (
Bounds, such as historical dates, TPQ or TAQ, may also be introduced in order to constrain one or several event dates
This condition must be included in the set of constraints that define the support of the prior distribution of the event dates.
The bound can be also defined with an uncertainty, i.e.,
where
and
It is important to note that the introduction of bounds
In this section we discuss the question of the specification of the prior on the event dates
A span of 2Δ is favored over a span of Δ by a factor of about Δ
Naylor and Smith (1988), Buck
In the absence of supplementary information, a non informative prior density is assigned to
with
From these hypotheses, it is possible to calculate the prior joint probability density of the event dates (
Finally, we have:
The joint probability density for
where
Starting from the prior densities (
These distributions are shown in Figure 4. When
Instead of looking at span Δ, we can look at the variance of the event dates
obtained with
Note that we simulate a sample from the prior distribution of (
Generate
Take
Generate
The same result is obtained whatever the partial or total order of the event dates.
The concentration effect is particularly visible when considering dates with different uncertainties. In this case, the posterior results for
Our result shows that the prior on event dates
We consider a phase with 5 non ordered Gaussian dates: 3 dates 0±30 AD, 1 date −500±200 AD and 1 date +500±200 AD. We can see in Figure 6(a) that the posterior densities for
In conclusion, the specification of an appropriate prior on target event dates remains an open question. We showed that the NSBC phase model favors, in an underlying way, the concentration of the event dates in a phase. We also showed that the uniform prior favors spreading dates, especially when temporal order constraints act.
In this section, Monte Carlo experiments are done to illustrate the performance of the proposed model on simulated samples. All the simulations are done using the R package ArchaeoChron (Philippe and Vibet, 2017c).
To illustrate the properties of the event date modeling, especially its robustness, we simulate sample of measurements
The date of the target event is 0, and
the r-combine model (Remark 1).
the
We assume that all the measurements have the same calibration curve, this condition is required for the r-combine model. We take
The lengths of the boxplots also indicate that the loss of accuracy if less important for event model than the
In this part, we illustrate the behavior of our chronological model in presence of stratigraphic inversion. We want to date target events that satisfy stratigraphy constraints
We simulate sample of 10 measurements such that the true dates satisfy the condition
Figure 9 gives the results for the target event model. The value of
Figure 10 gives the results of the
The model is implemented in the cross-platform application ChronoModel which is free and open source software (Lanos
The aim is to date the last firing of a medieval potter’s kiln recovered at the Maison-de-Retraite-Publique site (Mennessier-Jouannet
The observations are composed of 3 TL dates (CLER 202a, 202b, 203), 2 AM dates (inclination, declination), and 1 radiocarbon date (Ly-5212) As the kiln belongs to the historical period, the prior time range
In some favorable archaeological contexts, it is possible to get information about the duration of a phase. It is the case here with a household cluster excavated on the Malpaís Prieto site (Michoacán State, Mexico), in the framework of the archaeological project Uacusecha (Pereira
The simple modeling does not allow to significantly improve the prior archaeological information. To do this, various archaeological evidences can better constrain the occupation duration: stratigraphic evidence, accumulation processes of the occupation remains, and durability of the partly perishable architecture. In this example, these information indicate that the occupation has hardly exceeded one century, the most plausible estimation being between 60 and 90 years. As a trial, we consider that the phase duration cannot exceed
Posterior densities of the target event dates
The Bayesian event date model aims to estimate the date
This project is supported by the grant ANR-11-MONU-007 ChronoModel. The authors expressed their gratitude to Grégory Pereira (CNRS, UMR ArchAm, University of Paris 1 Panthéon-Sorbonne) for the exchange of information about the example of Malpaís Prieto site in Mexico.
We also thank Marie Anne Vibet for her very helpful discussions during the preparation of this paper.
The posterior distributions of the parameter of interest
Here we give a more detailed insight on the way to estimate the date
where
To estimate the posterior density of
MH-1: A Metropolis-Hastings algorithm where the proposal is the prior distribution of the parameter. This method is recommended when no calibration is needed, namely for TL/OSL, Gaussian measurements or typo-chronological references.
MH-2: A Metropolis-Hastings algorithm where the proposal is an adaptive Gaussian random walk. This method is adapted when the density to be approximated is unimodal. The variance of this proposal density is adapted during the process.
MH-3: Metropolis-Hastings algorithm where the proposal mimics the individual calibration density. This method is adapted for multimodal densities, such as calibrated measurements.
We are frequently confronted with multimodal target distributions as, for instance, in archaeomagnetic dating (Example 6.1). In this case, algorithms MH-1 and MH-2 are not well adapted to ensure a good mixing of the Markov chain. Consequently, an alternative is to choose the MH-3 algorithm with a proposal distribution that mimics the individual calibration density defined in (
where
where (
An alternative is to choose the distribution with density
where (
Using the packages
Figures A.1 and A.2 provide graphical tools for the diagnostic of MCMC sampler. We consider the example of Lezoux and only give the result for the parameter of interest
Within the Gibbs sampler, the
The Gelman diagnostic (evaluated from 5 parallel chains) is equal to 1, and confirms the convergence of MCMC samplers.
NSBC = Naylor-Smith-Buck-Christen.
95% highest posterior density region for estimating the phase (begin/end/duration)
Phase | |
---|---|
Begin | [1209; 1355] |
End | [1372; 1450] |
Duration | [48; 211] |
95% highest posterior density regions for estimating the phase (begin/end/duration) (
Phase | |
---|---|
Begin | [1286; 1397] |
End | [1335; 1443] |
Duration | [28; 60] |