Dates

Samples can be dated using historical information (e.g., when fossils are found of a species that was introduced to a region at a known historical time), or using radiometric dating where the ratio of certain radioactive isotopes is measured, with this ratio changing in a predictable manner over time. If the calendar age ${\theta _i}$ (in cal BP, calendar years before AD 1950) of a sample ${x_i}$ is estimated by a measurement ${y_i}$ with a degree of uncertainty ${\sigma _i}$ , it is often assumed that this measurement is distributed normally around the (unknown) true age,

(1) $${y_i}\sim N\left( {{\theta _i},{\sigma _i}} \right).$$

If the age of a sample is known exactly, then ${y_i} = {\theta _i}$ , while in other cases we know that an event took place before a specific date, ${\theta _i} \gt {\theta _{i - 1}}$ , or after it, ${\theta _i} \lt {\theta _{i + 1}}$ .

For radiocarbon dates, we have a good understanding of how ${{\rm{\;}}^{14}}{\rm{C}}$ ages relate to cal BP years through the different calibration curves such as IntCal20 (Reimer et al., Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards, Friedrich, Grootes, Guilderson, Hajdas, Heaton, Hogg, Hughen, Kromer, Manning, Muscheler, Palmer, Pearson, van der Plicht, Reimer, Richards, Scott, Southon, Turney, Wacker, Adolphi, Büntgen, Capano, Fahrni, Fogtmann-Schulz, Friedrich, Köhler, Kudsk, Miyake, Olsen, Reinig, Sakamoto, Sookdeo and Talamo2020), Marine20 (Heaton et al., Reference Heaton, Köhler, Butzin, Bard, Reimer, Austin, Bronk Ramsey, Grootes, Hughen, Kromer, Reimer, Adkins, Burke, Cook, Olsen and Skinner2020), SHCal20 (Hogg et al., Reference Hogg, Heaton, Hua, Palmer, Turney, Southon, Bayliss, Blackwell, Boswijk, Bronk Ramsey, Petchey, Reimer, Reimer and Wacker2020) or postbomb curves (Hua et al., Reference Hua, Turnbull, Santos, Rakowski, Ancapichún, De Pol–Holz, Hammer, Lehman, Levin, Miller, Palmer and Turney2022). Calibration curves essentially consist of three columns with the first a range of cal BP ages $\theta $ , followed by their associated ${{\rm{\;}}^{14}}{\rm{C}}$ ages $\mu \left( \cdot \right)$ and finally by their uncertainties $\sigma \left( \cdot \right)$ . Then, radiocarbon ages are assumed to be normally distributed given the calibration curve’s ${{\rm{\;}}^{14}}{\rm{C}}$ age belonging to the calendar year of interest:

(2) $${y_i}\sim N\left( {\mu \left( {{\theta _i}} \right),\omega \left( {{\theta _i}} \right)} \right),$$

where $\omega \left( {{\theta _i}} \right) = \sqrt {{\sigma ^2}\left( {{\theta _i}} \right) + \sigma _i^2} $ combines the uncertainties of both the calibration curve and of the date (e.g., Christen and Litton, Reference Christen and Litton1995).

Outliers

Analysis of material dated through multiple radiocarbon measurements often shows a scatter beyond what the individual lab errors would suggest (Christen and Pérez, Reference Christen and Pérez2009; Scott et al., Reference Scott, Naysmith and Cook2018; Scott, Reference Scott2023). This overdispersion can be modelled in several ways (Christen, 1994a,b; Bronk Ramsey, Reference Bronk Ramsey2009b). The most common approach to modelling outliers is to enable a shift of the measured mean:

(3) $$y{{\rm{'}}_{\!\!i}}\sim N\left( {\mu \left( {{\theta _i}} \right) + {\psi _i}{\delta _i},\omega \left( {{\theta _i}} \right)} \right),$$

where ${\psi _i} = \left( {0,1} \right)$ is a ‘flag’ and ${\delta _i} = N\left( {0,\beta \sigma _i^2} \right)$ is a shift on the radiocarbon scale (often, $\beta $ is set to 2; shifts on the cal BP scale can also be modelled). If a date is not outlying, ${\psi _i}$ is set to 0 and the date is treated as usual. If the date is an outlier, ${\psi _i}$ is set to 1 and a shift of its age is modelled. This prior outlier probability is often set to 5% (e.g., Blaauw and Christen, Reference Blaauw and Christen2005; Haslett and Parnell, Reference Haslett and Parnell2008; Bronk Ramsey, Reference Bronk Ramsey2009b) although Haslett and Parnell (Reference Haslett and Parnell2008) model two types of outliers (moderate outliers with a 5% prior probability and more extreme outliers with a 1% prior probability). This ‘flag and shift’ approach requires several additional parameters to be modelled for each date in a site. Posterior outlier probabilities are then calculated as the proportion of times during modelling that the outlier flag had to be set to 1. It is thus usually not necessary or recommendable to manually remove outliers - they are downweighted automatically based on the other dates and the modelling constraints.

Another approach to overdispersion is to treat the uncertainty of each date not as a known single value but as a distribution (Christen and Pérez, Reference Christen and Pérez2009). For example, if a lab reports a date to have an uncertainty of, say, 20 years, the true uncertainty could be lower or higher than the reported value. By multiplying the lab error by a gamma distribution, in effect we are replacing the normal distribution of Eq. (2) with a Student-t distribution, with parameters $t.a$ and $t.b$ . Depending on the parameter values used, the latter distribution will appear very similar to the normal distribution, but will have wider tails. $Bacon$ (Blaauw and Christen, Reference Blaauw and Christen2011) uses the Student-t distribution to model dates, as it makes the age-model robust to dating scatter without the need to model additional parameters (see also Blaauw et al., Reference Blaauw, Christen, Bennett and Reimer2018). We will be using the Student-t distribution in the approach developed below. Also with this approach, usually there is no need to manually remove outlying dates.

Stratigraphy

The dates can be integrated with any information regarding the time spans between (dated or undated) sample positions. For example, if we have radiocarbon samples from a tree with $n$ reliable yearly rings but no master dendrochronology to align it with, we can still precisely count the exact gaps, or time spans, between specific (sections of) tree-rings, and combine these with radiocarbon–dated rings to produce a chronology. If say 10 annual tree-rings were deposited between the older sample 1 and the younger sample 2, then we know that whatever the known or unknown age is of sample 1, it must be 10 years older than that of sample 2. In statistical terms (e.g., Christen, 1994b; Christen and Litton, Reference Christen and Litton1995), this can be written as there being $m$ dated positions in a tree, with each position ${x_i}$ having an associated calendar age ${\theta _i}$ separated by known integer time-spans from other rings within the tree (here we start counting from the youngest, outermost ring):

(4) $${\delta _i} = {\theta _i} - {\theta _{i - 1}} \in \mathbb{N}$$

Some rings will have one or more radiocarbon dates and some none. We will leave this implicit in the definition of the likelihood in Eq. (7) to avoid a more cumbersome notation. Other positions of interest, while not radiocarbon dated, may still require an age estimate. Given ${\theta _0}$ (the age of the outermost ring), the age of any ring can then be calculated by adding up the gaps:

(5) $${\theta _i} = {\theta _0} + \sum\limits_{j = 1}^i {} {\delta _j}.$$

In sites with less exact information on the time spans between levels, but where we can still safely assume a chronological ordering of the different levels based on the site’s stratigraphy (e.g., Harris, Reference Harris1989), often all we know is that between two stratigraphically separated levels ${x_i}$ and ${x_{i - 1}}$ , the time span ${\theta _i} - {\theta _{i - 1}} = {\delta _i} \gt 0$ . Some layers will not have such information and then dates within said layer cannot safely be assumed to be in chronological order (we will call these layers ‘blocks’). In some cases however, there could be additional information, say, if two samples are separated by a soil which is known from other sources to take a certain amount of time to form. These stratigraphical gaps could be modelled using a range of distributions, e.g., uniform: ${\delta _i}\sim U\left( {{a_i},{b_i}} \right)$ (see Bronk Ramsey, Reference Bronk Ramsey2000; Steier and Rom, Reference Steier and Rom2000), normal: ${\delta _i}\sim N\left( {{\mu _i},{\sigma _i}} \right)$ , gamma: ${\delta _i}\sim {\rm{\Gamma }}\left( {{\mu _i},{\sigma _i}} \right)$ , or even an exact gap: ${\delta _i} = {\mu _i}$ . (For sites such as bogs, lakes or marine deposits which can reliably be assumed to have accumulated approximately continuously over time, age–depth models such as BChron (Haslett and Parnell, Reference Haslett and Parnell2008), Bacon (Blaauw and Christen, Reference Blaauw and Christen2011) or OxCal’s P_Sequence (Bronk Ramsey, Reference Bronk Ramsey2009b) would be more appropriate.)

Combining dates and stratigraphy

Trees with their exactly known time-gaps between ${{\rm{\;}}^{14}}{\rm{C}}$ -dated rings are relatively easy to model. Combining Equations (2) and (5), we can calculate the likelihood of a calendar age for each ${{\rm{}}^{14}}{\rm{C}}$ -dated tree-ring ${y_i}\sim N\left( {\mu \left( {{\theta _i}} \right),\omega _i^2} \right)$ where, as before, ${\theta _i} = {\theta _0} + \sum\nolimits_{j = 1}^i {} {\delta _j}$ and $\omega _i^2 = \sqrt {{\sigma ^2}\left( {{\theta _i}} \right) + \sigma _i^2} $ . For any chosen calendar year ${\theta _0}$ , the likelihood of each ${{\rm{\;}}^{14}}{\rm{C}}$ -dated ring is calculated as $p({x_i}|{\theta _i})$ , and we can use the formula for the normal distribution:

(6) $$p({x_i}|{\theta _i}) = {1 \over {{\omega _i}\left( {{\theta _i}} \right)\sqrt {2\pi } }}\,{\rm{exp}}\left\{ { - {{{{({x_i} - \mu \left( {{\theta _i}} \right))}^2}} \over {2\omega _i^2\left( {{\theta _i}} \right)}}} \right\}.$$

Here, the vertical bar means ‘given’, so, we calculate the probability of a chosen age ${x_i}$ given the (unknown) true age ${\theta _i}$ (e.g., Christen and Litton, Reference Christen and Litton1995). As outlined above, instead of the normal distribution, the Student-t distribution could also be used. Now, we calculate the product of all likelihoods,

(7) $$p({\boldsymbol x}|{\boldsymbol \theta} ) = \prod _{i = 1}^mp({x_i}|{\theta _i}),$$

where ${\boldsymbol x} = \left( {{x_1},{x_2}, \ldots, {x_n}} \right)$ and ${\boldsymbol \theta} = \left( {{\theta _1},{\theta _2}, \ldots, {\theta _n}} \right)$ .

Bayes’ Theorem (Buck et al., Reference Buck, Cavanagh and Litton1996; Bronk Ramsey, Reference Bronk Ramsey2009a) states that

(8) $$p(\theta |x) \propto p(x|\theta )p\left( \theta \right).$$

In words, the posterior probability (our distribution of interest, here the age of a sample given the data) is proportional to the probability of the data given the age (see above), times the prior. This brings us to the prior information.

Prior information is any knowledge we have that could help interpreting the data — information that existed prior to obtaining or seeing the dating results. In the case of producing chronologies, this is the relative information between the dated levels, such as Eq. (5) for tree-rings, or any knowledge about the stratigraphical ordering of sites (e.g., Harris, Reference Harris1989; Buck et al., Reference Buck, Kenworthy, Litton and Smith1991; Christen and Litton, Reference Christen and Litton1995; Bayliss, Reference Bayliss2009). A very basic prior could be to constrain the ages to be within a certain timeframe, e.g., if we know that a site is of Holocene age, then we do not need to consider pre-Holocene ages, and radiocarbon-derived ages should be within the ${{\rm{\;}}^{14}}{\rm{C}}$ dating limit (< 55,000 cal BP).

For trees with year-rings, the prior information could take either of two values (Christen and Litton, Reference Christen and Litton1995):

(9) $$p({\boldsymbol \theta} ) = \left \{\matrix{ 1 & {0 \lt {\theta _0} \lt \infty ;{\theta _i} = {\theta _0} + \sum\nolimits_{j = 1}^i {} {\delta _j}\forall i = 1,2,...,m} \cr 0 & {\rm otherwise}. \hfill\cr } \right. $$

So, if a proposed modelling solution lies within a timeframe and has the correct exact time-gaps between all ${{\rm{\;}}^{14}}{\rm{C}}$ -dated rings, then all prior information is obeyed, and the prior becomes 1. If however any of these constraints are violated, the prior becomes 0.

For stratigraphical deposits, as outlined above, we usually assume chronological ordering (unless other information is available). For this, we measure the timespans ${\delta _i}$ between all neighboring ages ${\theta _i}$ and ${\theta _{i - 1}}$ , check that they are positive (i.e., in chronological order or ${\theta _{i - 1}} \lt {\theta _i}$ ), and calculate the corresponding probability ${p_i}\left( {{\theta _i} - {\theta _{i - 1}}} \right)$ , where ${p_i}$ is some prior distribution for the (positive) length of time-span ${\delta _i}$ . The joint prior is calculated by multiplying all individual probabilities (i.e. a priori independent):

(10) $$p({\boldsymbol \theta} ) = \left \{\matrix{ {\prod _{i = 1}^n} & {p_i({\theta _i} - {\theta _{i - 1}})\;{\rm{for}}\;{\theta _0} \lt {\theta _1} \lt ... \lt {\theta _n}} \cr 0 & {\rm otherwise}. \hfill\cr } \right. $$

The prior gap probabilities ${p_i}$ depend on both the gap length and on the chosen distributions. Most commonly used in stratigraphies are uniform priors where any gap length between certain limits has the same prior probability. Since Jones and Nicholls (Reference Jones and Nicholls2002) show that uniform prior densities favour a span of $2\delta $ over that of $\delta $ by a factor of c. ${2^{M - 1}}$ , the uniform distributions are weighted by the amount of events within the time-span (see Bronk Ramsey, Reference Bronk Ramsey2000; Steier and Rom, Reference Steier and Rom2000; Jones and Nicholls, Reference Jones and Nicholls2002):

(11) $${p_i}({\delta _i}) = \left \{\matrix{ {{1 \over {(R - {\delta _i})}}{1 \over {{{({\delta _i})}^{M - 1}}}},\;{\rm if}\;0 \lt {\delta _i} \lt R} \cr 0 \;{\rm otherwise,}\hfill \cr } \right.$$

where $R = P- A$ is the total possible time span based on the upper and lower bounds, and $M$ is the number of phases within the gap. If the gap length is assumed to be normally distributed as ${\mu _i} \pm {\sigma _i}$ , we use

(12) $${p_i}({\delta _i}|{\mu _i},{\sigma _i}) = {1 \over {{\sigma _i}\sqrt {2\pi } }}{\rm{\;exp}}\left\{ { - {{{{({\delta _i} - {\mu _i})}^2}} \over {2\sigma _i^2\left( {{\delta _i}} \right)}}} \right\}.$$

A Gamma distribution ensures that ${\delta _i} \gt 0$ :

(13) $$p_i({\delta _i}|{\alpha _i},{\beta _i}) = {{\beta _i^{{\alpha _i}}} \over {\Gamma \left( {{\alpha _i}} \right)}}\delta _i^{{\alpha _i} - 1}{e^{ - {\beta _i}{\delta _i}}}$$

with mean ${\mu _i} = {{{\alpha _i}} \over {{\beta _i}}}$ and variance $\sigma _i^2 = {{{\alpha _i}} \over {\beta _i^2}}$ . For exact gaps, we use a Dirac’s delta:

(14) $$p_i({\delta _i}) = \left \{\matrix{ 1 & {{\delta _i} = {\mu _i}} \cr 0 & {\rm otherwise} \cr } \right.$$

Implementation

Let ${\theta _{0,j}}$ be a calendar age assigned to the outermost, youngest ring ${x_0}$ of a ${{\rm{\;}}^{14}}{\rm{C}}$ -dated tree. Since this youngest ring now has an assigned age, all of the ${{\rm{\;}}^{14}}{\rm{C}}$ -dated levels ${x_1},{x_2}, \ldots, {x_m}$ too have their corresponding ages ${\theta _1},{\theta _2}, \ldots {\theta _m}$ (after Eq. 5, which ensures the correct exact time gaps), and these ages are used to calculate the ‘fit’ of the ${{\rm{\;}}^{14}}{\rm{C}}$ ages (Eq. 6) given the age ${\theta _{0,j}}$ . The product (Eq. 7) gives the ‘fit’ of the entire ${{\rm{\;}}^{14}}{\rm{C}}$ -dated tree. Now the same calculations are repeated for a sufficiently wide and dense range of calendar ages ${\theta _{0,j}}$ , providing a ‘wiggle-matched’ age distribution for the ${{\rm{\;}}^{14}}{\rm{C}}$ -dated tree (Christen and Litton, Reference Christen and Litton1995; Bronk Ramsey et al., Reference Bronk Ramsey, van der Plicht and Weninger2001).

For sedimentary sites, we resort to Markov Chain Monte Carlo (MCMC) (Robert and Casella, Reference Robert and Casella2004; Christen and Fox, Reference Christen and Fox2010). The process starts with a simulated set of ages through Monte Carlo simulation of parameters from the prior distribution. Namely, initially an age ${\theta _{0,j}}$ is assigned to the topmost, youngest level ${x_0}$ . This age is sampled from a uniform prior distribution with very wide boundaries (e.g., across the entire range of ${{\rm{\;}}^{14}}{\rm{C}}$ dating). Then, an age gap ${\delta _{i,j}}$ is simulated between levels ${x_0}$ and ${x_1}$ (this will provide a simulated age ${\theta _{1,j}} = {\theta _{0,j}} + {\delta _{1,j}}$ for ${x_1}$ ). This process is repeated for all levels in the site, providing gap and age estimates for all levels. For most levels, the gaps will be uniform with ${\delta _{i,j}} \gt 0$ to ensure chronological ordering, but levels within ‘blocks’ will allow for ${\delta _{i,j}} \le 0$ , and other gaps will need sampling from a normal or gamma distribution or require an exact gap. This series of ages and gaps then provides the first of $j = 1:k$ iterations of modelled years ${\theta _{0:n,k}}$ — but many more will be needed.

From this initial set, randomly one of the ages will be chosen and changed somehow. This becomes a ‘proposal’. For this proposal, the ‘energy’ will be calculated as the product of the likelihoods of the dates given the assigned ages (according to Eq. 6), times the product of all priors (the gaps, according to Eq. 11-13 as required). The proposal will also be checked to ensure that it still obeys all constraints (e.g., no reversals). Now, this proposed iteration may either be accepted as the next iteration, or it may be rejected (in which case, the previous iteration is copied to be the next iteration) using a Metropolis-Hastings acceptance probability (Robert and Casella, Reference Robert and Casella2004).

Over time, the MCMC process provably approximates the underlying distribution of all parameters involved (Robert and Casella, Reference Robert and Casella2004) — that is, the ages and gaps. Since each iteration depends on the previous iteration, there is a large degree of ‘memory’ between neighbouring iterations. To remove this memory, after running a few hundreds of thousands to millions of iterations, ‘thinning’ is carried out, storing only 1 every say few dozen iterations (depending on the number of parameters that are estimated). Additionally, since the first iteration is simulated from the prior, it represents sub-optimal ‘ballpark’ values, and commonly the first dozens of MCMC iterations often fit poorly – therefore, this ‘burn-in’ will have to be removed as well (by default the first 100 iterations). It is important to check that after removing the burn-in and thinning, the time-series of iterations (e.g., any of the parameters, or, the energy), does not exhibit an initial clear jump (too short burn-in) or clear auto-correlation. Several implementations of the MCMC approach are available, and here we use the t-walk (Christen and Fox, Reference Christen and Fox2010) as it runs natively within $R$ and is flexible enough to work with many parameters without the need for expert tuning.

To evaluate the effectiveness of our model in enhancing the accuracy of dating, either in dendrochronology or stratigraphy, a ‘fit’ measure is calculated for each date, by finding the percentage of modelled ages that fall within any of the highest posterior density (hpd) ranges of the unmodelled date. By default, this is done on the 95% hpd ranges. For the tree, a goodness-of-fit is also calculated as the standardized offset (z-score) between the dates and the calibration curve’s ${{\rm{\;}}^{14}}{\rm{C}}$ ages (at the assigned ‘best’ calendar age for each date).