carbondate

An R package to analyse, and summarise, multiple radiocarbon (¹⁴C) determinations. The package provides two linked (but distinct) Bayesian approaches that can both be used to obtain rigorous and robust alternatives to summed probability distributions (SPDs):

Bayesian Non-Parametric Density Estimation (Modelling the joint, summary, calendar distribution as an unknown mixture of calendar age clusters) — see “Non-parametric calibration of multiple related radiocarbon determinations and their calendar age summarisation” (Heaton, 2022).
Variable-Rate Poisson Process Estimation (Modelling changes in the occurrence rate of the ¹⁴C samples over calendar time) — see “A New Approach to Radiocarbon Summarisation: Rigorous Identification of Variations/Changepoints in the Occurrence Rate of Radiocarbon Samples using a Poisson Process” (Heaton, Al-assam & Bard, 2025)

Both methods are implemented using Markov Chain Monte Carlo (MCMC). In terms of which method, non-parametric density estimation or poisson process estimation, is more appropriate for a given problem we would suggest that users might apply both approaches and compare outputs as part of the research process.

Installation

The easiest way to install the latest release is via CRAN, by typing the following into your R console:

install.packages("carbondate")

You can alternatively install the development version of carbondate from GitHub with:

# devtools::install_github("TJHeaton/carbondate")
install.packages("carbondate", repos = c('https://tjheaton.r-universe.dev', 'https://cloud.r-project.org'))

Once you have installed the library with either of the above methods, you need to load it using:

library(carbondate)

Data

There are a few example datasets of radiocarbon determinations (e.g., two_normals, kerr, pp_uniform_phase, buchanan, alces, equus, human, …) provided, which can be used to try out the calibration functions. two_normals is a small simulated dataset for which the underlying calendar ages were drawn from a known mixture of two normals. It is included simply to give a quick-to-run example for the Bayesian Non-Parametric Density calibration functions. pp_uniform_phase is another small simulated dataset, for which the underlying calendar ages were drawn uniformly at random from a short calendar interval (equivalent to a single uniform phase). This is included to give a quick-to-run example for the Poisson Process modelling functions. There are analogous simulated datasets from marine environments two_normals_marine and pp_uniform_phase_marine; and also a dataset containing simulated data from SH, NH, and Marine environments pp_uniform_phase_mixed. The remaining datasets are from real-life data. The Northern Hemisphere IntCalXX calibration curves, Southern Hemisphere SHCalXX calibration curves, and surface-ocean MarineXX calibration curves are also provided (the most recent curves, published in 2020, are intcal20,shcal20, and marine20).

Bayesian Non-Parametric Summarisation: Quick-Start Example

The below example implements the Bayesian non-parametric calibration and summarisation approach on the simulated two_normals data using the IntCal20 curve via two (slightly) different MCMC methods.

polya_urn_output <- PolyaUrnBivarDirichlet(
  rc_determinations = two_normals$c14_age,
  rc_sigmas = two_normals$c14_sig,
  calibration_curve = intcal20)

walker_output <- WalkerBivarDirichlet(
  rc_determinations = two_normals$c14_age,
  rc_sigmas = two_normals$c14_sig,
  calibration_curve = intcal20)

Once the calibration has been run, the resultant estimate for the joint (summary) calendar age density can be plotted.

PlotPredictiveCalendarAgeDensity(
  output_data = list(polya_urn_output, walker_output),
  show_SPD = TRUE)

The summary estimates generated by the Bayesian non-parametric approach (shown in purple and green) are called predictive densities. They aim to predict the calendar age of a future (new) sample based upon the set of ¹⁴C samples that you have summarised.

Poisson Process Modelling: Quick-Start Example

The below example implements the Poisson process modelling approach using the simulated pp_uniform_phase data and the IntCal20 curve.

pp_output <- PPcalibrate(
  rc_determinations = pp_uniform_phase$c14_age,
  rc_sigmas = pp_uniform_phase$c14_sig,
  calibration_curve = intcal20)

Once the calibration has been run, the posterior estimate for the occurrence rate of the samples can be plotted.

PlotPosteriorMeanRate(output_data = pp_output)

This estimated sample occurrence rate, and changes within it, can be interpreted equivalently to methods of calendar age summarisation.