A. Philippe, M.-A. Vibet
chronos provides a list of functions for the statistical analysis of archaeological dates and groups of dates. It is based on the post-processing of the Markov Chains whose stationary distribution is the posterior distribution of a series of dates. Such MCMC output can be simulated by different applications as for instance ChronoModel (Lanos et al. 2020), Oxcal (Bronk Ramsey 2009) or BCal (Buck, Christen, and James 1999). The only requirement is to have a CSV file containing a sample from the posterior distribution.
## Install the latest version install.packages("fasti", repos = "https://tesselle.r-universe.dev")
Let’s use the data of Ksar Akil generated by ChronoModel (Bosch et al. 2015).
Two different files are generated by ChronoModel:
Chain_all_Events.csv that contains the MCMC samples of each event created in the modeling, and
Chain_all_Phases.csv that contains all the MCMC samples of the minimum and the maximum of each group of dates if at least one group is created.
chrono_path <- "chronomodel/ksarakil" ## Read events from ChronoModel output_events <- system.file(chrono_path, "Chain_all_Events.csv", package = "fasti") chrono_events <- read_chronomodel_events(output_events) ## Read phases from ChronoModel output_phases <- system.file(chrono_path, "Chain_all_Phases.csv", package = "fasti") chrono_phases <- read_chronomodel_phases(output_phases)
vignette("import") for more details on how to import MCMC samples.
For more details on the diagnostic of Markov chain, see Robert and Casella (2010).
To assess the agreement between the posterior distributions and the numerical approximations, three Markov chains were run in parallel by ChronoModel. For each chain, 1 000 iterations were used during the Burn-in period, 20 batches of 500 iterations were used in the Adapt period, 100 000 iterations were drawn in the Acquire period by only 1 out of 10 were kept in order to break the correlation structure.
From the analysis of the history plot, all Markov chains reach their equilibrium before the Acquire period. The autocorrelations of the three Markov chains are not significant, meaning the rate of subsample (1 over 10) is enough.
Now, using the package chronos and the package coda, we can verify whether the MCMC samples are correctly generated by the software.
## Create an mcmc.list coda_events <- as_coda(chrono_events, chains = 3)
Indeed, the MCMC samples should have no autocorrelation and should have reached their equilibrium (that is the posterior density of the parameter under investigation).
The autocorrelation plots show that each of these three chains are not significant. That means that we actually generated a non correlated sample, which was the aim the MCMC process.
We can also check whether the chains reached equilibrium. For example, let’s consider the first date of the dataset.
plot(coda_events[, 1, ], trace = TRUE, density = FALSE)
The plot shows that the three chains corresponding to the first date reached the same stationary process.
We can test the Gelman-Rubin criterion. The expected value to confirm that all of the Markov chains reached equilibrium is 1.
coda::gelman.diag(coda_events) #> Potential scale reduction factors: #> #> Point est. Upper C.I. #> V 1 1 #> VI 1 1 #> XI 1 1 #> XII 1 1 #> XVI_1 1 1 #> XVI_2 1 1 #> XVI_3 1 1 #> XVI_4 1 1 #> XVII_1 1 1 #> XVII_2 1 1 #> XVII_3 1 1 #> XVII_4 1 1 #> XVIII 1 1 #> XIX 1 1 #> XX 1 1 #> XXII 1 1 #> #> Multivariate psrf #> #> 1
The Gelman-Rubin criterion confirms that all of the Markov chains reached equilibrium.
We can also test the Geweke criterion. The expected value to confirm that all of the Markov chains reached equilibrium is strickly less than 1.
coda::geweke.diag(coda_events[, 1, ], frac1 = 0.1, frac2 = 0.5) #> [] #> #> Fraction in 1st window = 0.1 #> Fraction in 2nd window = 0.5 #> #> var1 #> 0.4978 #> #> #> [] #> #> Fraction in 1st window = 0.1 #> Fraction in 2nd window = 0.5 #> #> var1 #> 0.7736 #> #> #> [] #> #> Fraction in 1st window = 0.1 #> Fraction in 2nd window = 0.5 #> #> var1 #> 0.7666
The Geweke criterion criterion confirms that all of the Markov chains reached equilibrium.
ChronoModel generated correct samples of the posterior distribution. Now gathering the three chains, a total of 30 000 iterations was collected in order to give estimations of the posterior distribution of each parameter.
plot(chrono_events) #> Picking joint bandwidth of 49.2
The tempo plot has been introduced by Dye (2016). See Philippe and Vibet (2020) for more statistical details.
The tempo plot is one way to measure change over time: it estimates the cumulative occurrence of archaeological events in a Bayesian calibration. The tempo plot yields a graphic where the slope of the plot directly reflects the pace of change: a period of rapid change yields a steep slope and a period of slow change yields a gentle slope. When there is no change, the plot is horizontal. When change is instantaneous, the plot is vertical.
From these graphs, we can see that the highest part of the sampled activity is dated between -45 000 to -35 000 but two dates are younger, at about -32 000 and -28 000.
The activity plot displays the derivative of the Bayes estimate of the Tempo plot. It is an other way to see changes over time.
The Occurrence plot calculates the calendar date \(t\) corresponding to the smallest date such that the number of events observed before \(t\) is equal to \(k\), for \(k = [(1, 16)]\). The Occurrence plot draws the credible intervals or the highest posterior density (HPD) region of those dates associated to a desired level of confidence.
A group of dates (phase) is defined by the date of the minimum and the date of the maximum of the group. In this part, we will use the data containing these values for each group of dates.
We can estimate the time range of a group of dates as the shortest interval that contains all the dates of the group at a given confidence level (Philippe and Vibet 2020).
The following code gives the endpoints of the time range of all groups of dates of Ksar Akil data at a given confidence level.
boundaries(chrono_groups, level = 0.95) #> lower upper #> EPI -28978.53 -26969.82 #> UP -38570.37 -29368.75 #> Ahmarian -42168.47 -37433.31 #> IUP -43240.37 -41161.00
The time range interval of the group of dates is a way to summarize the estimation of its minimum, the estimation of its maximum and their uncertainties at the same time.
plot(chrono_groups, level = 0.95)
We may also be interested in a succession of phases. This is actually the case of the succession of IUP, Ahmarian, UP and EPI that are in stratigraphic order. Hence, we can estimate the transition interval and, if it exists, the gap between these successive phases.
The transition interval between two successive phases is the shortest interval that covers the end of the oldest group of dates and the start of the youngest group of dates. The start and the end are estimated by the minimum and the maximum of the dates included in the group of dates. It gives an idea of the transition period between two successive group of dates. From a computational point of view this is equivalent to the time range calculated between the end of the oldest group of dates and the start of the youngest group of dates.
transition(chrono_groups, level = 0.95) #> An object of class "TimeRange" #> Slot "lower": #> EPI UP Ahmarian IUP #> EPI NA NA NA NA #> UP -31479.80 NA NA NA #> Ahmarian -39138.82 -39118.07 NA NA #> IUP -43487.53 -43395.89 -43212.31 NA #> #> Slot "upper": #> EPI UP Ahmarian IUP #> EPI NA NA NA NA #> UP -26905.04 NA NA NA #> Ahmarian -27122.05 -36741.08 NA NA #> IUP -26866.99 -36480.26 -40733.77 NA #> #> Slot "names": #> EPI UP Ahmarian IUP #> EPI "EPI-EPI" "EPI-UP" "EPI-Ahmarian" "EPI-IUP" #> UP "UP-EPI" "UP-UP" "UP-Ahmarian" "UP-IUP" #> Ahmarian "Ahmarian-EPI" "Ahmarian-UP" "Ahmarian-Ahmarian" "Ahmarian-IUP" #> IUP "IUP-EPI" "IUP-UP" "IUP-Ahmarian" "IUP-IUP" #> #> Slot "calendar": #>  "CE" #> #> Slot "hash": #>  "496feac16967682860342e99fbd5189d"
Successive phases may also be separated in time. Indeed there may exist a gap between them. This testing procedure check whether a gap exists between two successive groups of dates with fixed probability. If a gap exists, it is an interval that covers the end of one group of dates and the start of the successive one with fixed posterior probability.
hiatus(chrono_groups, level = 0.95) #> An object of class "TimeRange" #> Slot "lower": #> EPI UP Ahmarian IUP #> EPI NA NA NA NA #> UP -29188.56 NA NA NA #> Ahmarian -37368.33 NA NA NA #> IUP -41220.64 -41282.64 NA NA #> #> Slot "upper": #> EPI UP Ahmarian IUP #> EPI NA NA NA NA #> UP -28961.79 NA NA NA #> Ahmarian -28884.37 NA NA NA #> IUP -28814.70 -38421.49 NA NA #> #> Slot "names": #> EPI UP Ahmarian IUP #> EPI "EPI-EPI" "EPI-UP" "EPI-Ahmarian" "EPI-IUP" #> UP "UP-EPI" "UP-UP" "UP-Ahmarian" "UP-IUP" #> Ahmarian "Ahmarian-EPI" "Ahmarian-UP" "Ahmarian-Ahmarian" "Ahmarian-IUP" #> IUP "IUP-EPI" "IUP-UP" "IUP-Ahmarian" "IUP-IUP" #> #> Slot "calendar": #>  "CE" #> #> Slot "hash": #>  "496feac16967682860342e99fbd5189d"
At a confidence level of 95%, there is no gap between the succession of phases IUP, Ahmarian and UP, but there exists one of 203 years between phase UP and phase EPI.
Bosch, M. D., M. A. Mannino, A. L. Prendergast, T. C. O’Connell, B. Demarchi, S. M. Taylor, L. Niven, J. van der Plicht, and J.-J. Hublin. 2015. “New Chronology for Ksâr ’Akil (Lebanon) Supports Levantine Route of Modern Human Dispersal into Europe.” Proceedings of the National Academy of Sciences 112 (25): 7683–8. https://doi.org/10.1073/pnas.1501529112.
Bronk Ramsey, C. 2009. “Bayesian Analysis of Radiocarbon Dates.” Radiocarbon 51 (1): 337–60. https://doi.org/10.1017/S0033822200033865.
Buck, C. E., J. A. Christen, and G. E. James. 1999. “BCal: An on-Line Bayesian Radiocarbon Calibration Tool.” Internet Archaeology 7. https://doi.org/10.11141/ia.7.1.
Dye, T. S. 2016. “Long-Term Rhythms in the Development of Hawaiian Social Stratification.” Journal of Archaeological Science 71: 1–9. https://doi.org/10.1016/j.jas.2016.05.006.
Lanos, Ph., A. Philippe, H. Lanos, and Ph. Dufresne. 2020. Chronomodel: Chronological Modeling of Archaeological Data Using Bayesian Statistics. https://www.chronomodel.fr.
Philippe, A., and M.-A. Vibet. 2020. “Analysis of Archaeological Phases Using the R Package Archaeophases.” Journal of Statistical Software, Code Snippets 93 (1): 1–25. https://doi.org/10.18637/jss.v093.c01.
Robert, C. P., and G. Casella. 2010. Introducing Monte Carlo Methods with R. Use R! New York: Springer.