Modelling Ideas

Suppose we have a series of \(n\) events, which occur at calendar times \(\{\theta_i, i= 1,\ldots,n \}\) cal yr BP. Each event is assumed to result in the generation of a radiocarbon sample, \(\{ X_i, i= 1,\ldots,n \}\). Examples of events might be:

• death of an animal (which leaves a bone)
• creation of a cave painting (which leaves a painting)
• construction of a settlement (which leaves a farmstead)
• a palaeofire (which leaves behind charcoal)

We model the occurrence of these events as random and aim to investigate whether, and how, the rate at which they occur varies over calendar time.

Specifically, we assume that the events occur according to a Poisson Process with a variable rate \(\lambda(t)\). This rate \(\lambda(t)\) is considered piecewise-constant, but changes an unknown number of times and at unknown calendar ages. We aim to estimate \(\lambda(t)\) and hence infer how the rate at which the events occur changes over calendar time. Time periods where the rate \(\lambda(t)\) is high will typically be expected to generate a greater number of events (perhaps since there is more activity, or a greater population, at the site/sites); while periods when the rate \(\lambda(t)\) is lower (perhaps due to lower activity at the site/sites) will be expected to have fewer events.

When the calendar ages of the events \(\{\theta_i , i=1,\ldots,n\}\) are known exactly, estimation of the variable rate \(\lambda(t)\) can be tackled directly via reversible jump Markov Chain Monte Carlo (RJ-MCMC) (Green 1995). However, when we only observe the radiocarbon ages of the samples (\(\{ X_i, i= 1,\ldots,n \}\)) the underlying calendar ages \(\theta_i\) are uncertain and must be estimated via radiocarbon calibration. This complicates the process of estimating \(\lambda(t)\) especially as the calibration must be done simultaneously to its estimation (as the modelling will influence the set of \(\theta_i\) and vice versa). In this library, we provide a rigorous approach to the estimation of \(\lambda(t)\).

Modelling sample occurence using a Poisson process with a latent occurence/activity rate. Left Panel: An illustration of a Poisson process with two changepoints in the occurrence rate. The samples (shown as a rug in purple) occur at random calendar times, but proportional to the underlying occurrence rate. There are therefore more samples between (2700, 2300) cal yr BP than other times. Right Panel: Each sample has a ¹⁴C age (shown as black ticks on the y-axis). We wish to reconstruct the underlying Poisson process rate given only these ¹⁴C determination.

An Inhomogeneous Poisson process

A Poisson process is a common statistical process used to model the occurrence of random events. It is described by specifying the random variables \(N(s,t)\) for \(0 \leq s \leq t\), where \(N(s,t)\) represents the number of events that have occurred in the time interval \((s,t]\). Any Poisson process has a parameter \(\lambda(t)\), known as the rate of the process. This controls the expected number of events in a particular time interval.

In the case of an inhomogeneous Poisson process, the rate \(\lambda(t)\) varies over time so that some time periods are expected to have a greater number of event than others. An (inhomogeneous) Poisson process is then defined by the following two properties:

• For any \(0 \leq s \leq t\), the distribution of \(N(s,t)\) is Poisson with parameter \(\int_{\theta = s}^t \lambda(\theta) d\theta\).
• If \((s_1,t_1] , (s_2,t_2 ], \ldots ,(s_n,t_n]\) are disjoint intervals (i.e., non-overlapping) then \(N(s_1,t_1),N(s_2,t_2 ),\ldots ,N(s_n,t_n)\) are independent random variables.

In particular, the first of these properties implies that the expected number of events in a specific time interval \((s,t]\) is \(\int_{\theta = s}^t \lambda(\theta) d\theta\). Hence time periods with higher rates \(\lambda(t)\) will be expected to have greater number of events. The second implies that, conditional on the rate \(\lambda(t)\), the number of events that have occurred in one time period is not affected by the number of events in a different, and non-overlapping time period. It also implies that events are not clustered beyond the effect of the variable rate, e.g., the occurrence of one event does not increase/decrease the probability that another event occurs shortly afterwards. Note that this last assumption is quite strong as clustering might occur in many instances, for example if the animal bones you date all result from the hunting of a single (short-lived) group of individuals, or you date multiple buildings in a settlement.

Final Radiocarbon Model

Given a set of multiple \(^{14}\)C samples, \[X_1, \ldots, X_n\] the function PPcalibrate will simultaneously:

• calibrate them to provide posterior calendar ages estimates \(\theta_1, \ldots, \theta_n\)
• estimate the rate \(\lambda(t)\) of the underlying process. This estimate \(\hat{\lambda}(t)\) may provide a useful proxy to make archaeological and/or environmental inference, e.g., to model changes in culture, or useage, or population size over time.

The (calendar) rate at which the samples occur is modelled as a Poisson process with a variable rate \(\lambda(t)\). Specifically, \(\lambda(t)\) is modelled as a piecewise constant function with an unknown number of changepoints: \[ \lambda(t) = \left\{ \begin{array}{ll} h_i \quad & \textrm{if } s_i \leq t < s_{i+1} \\ 0 & \textrm{otherwise} \end{array} \right. \] where \(T_A < s_1 < s_2 < \ldots < s_k < T_B\) are the calendar ages at which there are changes in the Poisson process rate, and \(T_A\) and \(T_B\) are bounding calendar ages. We use reversible jump MCMC (RJMCMC) to alternate between calibrating the samples and updating the estimate of the underlying rate lambda.

Example

This is a basic example which shows you how to solve a common problem. In this simulated example, we will assume that the events (samples) can must lie within the calendar period from [700, 300] cal yr BP. Within this period, we will create 40 artificial samples with \(^{14}\)C determinations corresponding to the (narrower) calendar period from [550, 500] cal yr BP. Effectively, the underlying model is a Poisson process with a rate \[ \lambda(t) \left\{ \begin{array}{ll} h_1 = 0 & \textrm{if } 300 \leq t < 500 \textrm{ cal yr BP} \\ h_2 \qquad & \textrm{if } 500 \leq t < 550 \textrm{ cal yr BP}\\ h_3 = 0 & \textrm{if } 550 \leq t < 700 \textrm{ cal yr BP} \end{array} \right. \] We will then investigate if we can reconstruct this underlying information, from just the \(^{14}\)C values, if we only knew that the calendar ages of the samples lay between the bounding [700, 300] cal yr BP.

set.seed(15)

# Set initial values 
n_observed <- 40
rc_sigmas <- rep(15, n_observed)

# Create artificial rc_determinations
calendar_age_range <- c(300, 700) 
observed_age_range <- c(500, 550) 
true_theta <- seq(
      from = observed_age_range[1],
      to = observed_age_range[2],
      length = n_observed)
intcal_mean <- approx(
      x = intcal20$calendar_age_BP,
      y = intcal20$c14_age,
      xout = true_theta)$y
intcal_sd <- approx(
      x = intcal20$calendar_age_BP,
      y = intcal20$c14_sig,
      xout = true_theta)$y
rc_determinations <- rnorm(
      n = n_observed,
      mean = intcal_mean,
      sd = sqrt(rc_sigmas^2 + intcal_sd^2))

# Fit the model
PP_fit_output <- PPcalibrate(
    rc_determinations = rc_determinations,
    rc_sigmas = rc_sigmas,
    calibration_curve = intcal20,
    calendar_age_range = calendar_age_range,
    calendar_grid_resolution = 1,
    n_iter = 1e5,
    n_thin = 10,
    show_progress = FALSE)

This creates an object PP_fit_output which you can access directly or use the in-built plotting functions.

PlotPosteriorMeanRate(PP_fit_output)

PlotNumberOfInternalChanges(PP_fit_output)

PlotPosteriorChangePoints(PP_fit_output)
#> Warning in PlotPosteriorChangePoints(PP_fit_output): No posterior samples with
#> 1 internal changes

PlotPosteriorHeights(PP_fit_output)
#> Warning in PlotPosteriorHeights(PP_fit_output): No posterior samples with 1
#> internal changes

PlotCalendarAgeDensityIndividualSample(
  9, PP_fit_output, show_hpd_ranges = TRUE, show_unmodelled_density = TRUE)

Changing the calendar age plotting scale

As elsewhere in the library, the calendar age scale (cal yr BP, cal AD, or cal BC) shown when plotting can be selected via the plot_cal_age_scale variable within the relevant functions. In the case of the Poisson process modelling, the calendar age scale can be set in PlotPosteriorMeanRate and PlotPosteriorChangePoints, e.g.,

PlotPosteriorChangePoints(PP_fit_output, 
                          plot_cal_age_scale = "AD")
#> Warning in PlotPosteriorChangePoints(PP_fit_output, plot_cal_age_scale = "AD"):
#> No posterior samples with 1 internal changes