Age-depth models and what they do to a proxy

palaeoecology

time series

uncertainty

ecology tutorial

A sediment core with a perfectly constant accumulation rate acquires a fluctuating influx record once dating error is turned into a chronology. Measured in R.

Author

Tidy Ecology

Published

2026-08-07

Modified

2026-09-30

Corrected 30 September 2026: the post said the ensemble interval undercovers for a structural reason common to all ensembles; it now says the shortfall comes from the simulated dating error not being centred on zero, which the ensemble inherits and doubles, and it no longer says that the ensemble is not what its label suggests or that more dates would improve its centring. It also said a gradual transition reverses the balance between proxy and chronology; it now says that reversal comes from a few runs in which the simple detector jumped to an end of the core.

Every record that comes out of a core is measured against depth and reported against time, and the thing that converts one into the other is an age-depth model fitted to a handful of dated levels. That model is an estimate, built from six or eight ages with errors of a century each, and everything downstream inherits it: accumulation rates, influx, the date of a transition, the rate of change either side of it.

The inheritance is easy to miss because the chronology is usually built once, early, by someone else, and then used as though it were a measurement. This post takes a core where the truth is known, puts realistic dating error on the dated levels, and measures what comes out the other end.

A core with nothing happening in it

The cleanest test is a core with a constant accumulation rate. Sediment arrives at the same speed from bottom to top, a proxy arrives at a constant concentration per unit volume, and therefore the influx per unit time is constant too. Anything the analysis reports as variation is manufactured.

library(ggplot2)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          plot.title       = element_text(colour = te_ink, face = "bold"),
          axis.text        = element_text(colour = te_body))
}

dated  <- c(15, 70, 130, 190, 245, 295)     # cm, the dated levels
sig_lo <- 60                                # dating error, years
sig_hi <- 110                               # longer tail towards older ages
conc   <- 1200                              # grains per cubic cm, constant
rate_a <- 0.08                              # cm per year, constant
age_a  <- function(d) d / rate_a

# calibrated ages are skewed, so use a two-piece normal, and keep draws in order
two_piece <- function(n, s_lo, s_hi)
  ifelse(runif(n) < s_lo / (s_lo + s_hi),
         -abs(rnorm(n, 0, s_lo)), abs(rnorm(n, 0, s_hi)))

draw_ages <- function(true_age) {
  repeat {
    a <- true_age + two_piece(length(true_age), sig_lo, sig_hi)
    if (all(diff(a) > 0) && a[1] > 0) return(a)
  }
}

set.seed(44)
obs <- draw_ages(age_a(dated))
round(rbind(depth = dated, true_age = age_a(dated), dated_age = obs))

          [,1] [,2] [,3] [,4] [,5] [,6]
depth       15   70  130  190  245  295
true_age   188  875 1625 2375 3062 3688
dated_age  279  803 1661 2455 3053 3591

Six levels, each out by up to about a century, and every one of them within what a laboratory would report. The chronology is then whatever function you choose to run through those points. Linear interpolation and a monotone spline are the two that dominate the literature, and both are one line.

slices <- seq(15, 295, by = 10)             # the sampled slices

lin <- approxfun(c(0, dated), c(0, obs), rule = 2)
spl <- splinefun(c(0, dated), c(0, obs), method = "monoH.FC")

influx <- function(f) conc * diff(slices) / diff(f(slices))
i_lin  <- influx(lin)
i_spl  <- influx(spl)
i_true <- conc * rate_a

round(c(true_influx = i_true,
        linear_cv   = sd(i_lin) / mean(i_lin),
        spline_cv   = sd(i_spl) / mean(i_spl),
        linear_range = max(i_lin) / min(i_lin),
        spline_range = max(i_spl) / min(i_spl)), 3)

 true_influx    linear_cv    spline_cv linear_range spline_range 
      96.000        0.145        0.182        1.501        1.898

The true influx is 96 grains per square centimetre per year at every depth in the core. The linear chronology reports a series with a coefficient of variation of 0.145 and a 1.50-fold range from its lowest slice to its highest. The spline reports 0.182 and 1.90.

mid <- slices[-length(slices)] + diff(slices) / 2
dd  <- rbind(data.frame(depth = mid, influx = i_lin, model = "linear interpolation"),
             data.frame(depth = mid, influx = i_spl, model = "monotone spline"))

ggplot(dd, aes(depth, influx, colour = model, shape = model)) +
  geom_hline(yintercept = i_true, colour = te_ink, linewidth = 0.7) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 2.2) +
  geom_rug(data = data.frame(depth = dated), aes(x = depth), inherit.aes = FALSE,
           sides = "t", colour = te_ink, linewidth = 0.8,
           length = unit(0.035, "npc")) +
  annotate("text", x = 292, y = i_true, hjust = 1, vjust = -0.7, size = 3.4,
           colour = te_ink, label = "the truth") +
  scale_colour_manual(values = c("linear interpolation" = te_rust,
                                 "monotone spline" = te_forest)) +
  scale_shape_manual(values = c(16, 17)) +
  labs(x = "depth, cm", y = "influx, grains per square cm per year",
       colour = NULL, shape = NULL,
       title = "A flat record with structure in it") +
  theme_datasheet() +
  theme(legend.position = "top")

Two curves of influx against depth wandering above and below a flat horizontal truth line. The linear curve holds constant values between dated levels and jumps at them; the spline curve bends smoothly and overshoots sharply near the top of the core. Six tick marks along the top mark the dated levels. — Figure 1: Influx along a core whose true accumulation rate is constant. The horizontal line is the truth. Both chronologies were fitted to the same six dated levels, marked along the top.

The shape of the artefact is worth reading off the picture. Linear interpolation holds the accumulation rate constant between dated levels and changes it abruptly at each one, so the influx record is a staircase whose steps sit exactly where the dates are. The spline spreads the same error out into smooth undulations, which look far more like an environmental signal and are for that reason more dangerous, and it overshoots near the top of the core where it has an anchor on one side only.

One core is one draw of the dating error, so repeat it.

cv_and_range <- function(f) {
  v <- conc * diff(slices) / diff(f(slices))
  c(cv = sd(v) / mean(v), range = max(v) / min(v))
}

set.seed(31)
core_a <- t(replicate(500, {
  o <- draw_ages(age_a(dated))
  c(lin = cv_and_range(approxfun(c(0, dated), c(0, o), rule = 2)),
    spl = cv_and_range(splinefun(c(0, dated), c(0, o), method = "monoH.FC")))
}))
round(rbind(median = apply(core_a, 2, median),
            upper_decile = apply(core_a, 2, quantile, 0.9)), 3)

             lin.cv lin.range spl.cv spl.range
median        0.154     1.543  0.182     1.856
upper_decile  0.253     2.004  0.323     2.844

Across 500 cores with a genuinely constant accumulation rate, the median manufactured coefficient of variation is 0.154 for linear interpolation and 0.182 for the spline, against a truth of zero. One core in ten reports a range of 2.84-fold or more from the spline. A manufactured swing of 18 per cent in influx is the size of result that gets a paragraph in a results section.

What it does to a change that is really there

Manufacturing variation from nothing is one failure. The other is getting the size of a real change wrong, and that one is harder to notice because the change is genuinely present.

rate_b <- function(d) ifelse(d > 120 & d < 200, 0.02, 0.12)   # a slow section
dg <- seq(0, 300, by = 0.25)
age_b <- approxfun(dg, cumsum(c(0, diff(dg) / rate_b(dg[-1]))))

true_b   <- conc * diff(slices) / diff(age_b(slices))
true_amp <- max(true_b) / min(true_b)

set.seed(77)
core_b <- t(replicate(500, {
  o <- draw_ages(age_b(dated))
  c(lin = cv_and_range(approxfun(c(0, dated), c(0, o), rule = 2))[["range"]],
    spl = cv_and_range(splinefun(c(0, dated), c(0, o), method = "monoH.FC"))[["range"]])
}))
round(c(true_amplitude = true_amp,
        linear_median  = median(core_b[, "lin"]),
        linear_p90     = quantile(core_b[, "lin"], 0.9),
        spline_median  = median(core_b[, "spl"]),
        spline_p90     = quantile(core_b[, "spl"], 0.9)), 2)

true_amplitude  linear_median linear_p90.90%  spline_median spline_p90.90% 
          6.00           7.22          11.35          11.81          33.50

The core now has a real slow section, and the true ratio between the highest and lowest influx across the sampled slices is 6.00. Linear interpolation reports a median of 7.22. The monotone spline reports 11.81, roughly double the truth, and in the upper decile of cores it reports 33.5.

The direction is not an accident. A spline forced through dated levels that are individually displaced has to bend to reach them, and the bending shows up as very low apparent accumulation somewhere, which turns into a very high apparent influx. The estimator that produces the smoother, better-looking chronology is the one that exaggerates the amplitude most.

The date of a transition is a chronological quantity

Ask when something happened and the answer has two sources of error: where in the core the transition is, and what date the chronology assigns to that depth. It is worth knowing which one dominates, because effort spent sampling the core more finely is wasted if the answer is the other one.

step_depth <- 160
samp <- seq(20, 290, by = 1)

detect <- function(sharpness, noise) {
  y <- 1 + 2 * plogis((samp - step_depth) / sharpness) + rnorm(length(samp), 0, noise)
  k <- 6:(length(samp) - 5)
  gap <- vapply(k, function(i) mean(y[i:length(y)]) - mean(y[1:(i - 1)]), 1)
  samp[k[which.max(gap)]]
}

set.seed(88)
n_rep      <- 600
chron_only <- replicate(n_rep, approxfun(c(0, dated), c(0, draw_ages(age_a(dated))),
                                         rule = 2)(step_depth))
sharp_only <- replicate(n_rep, age_a(detect(0.4, 0.35)))
soft_only  <- replicate(n_rep, age_a(detect(3.5, 0.90)))

round(c(true_date = age_a(step_depth),
        sd_chronology = sd(chron_only),
        sd_sharp_transition = sd(sharp_only),
        sd_gradual_transition = sd(soft_only)), 1)

            true_date         sd_chronology   sd_sharp_transition 
               2000.0                  57.1                   6.5 
sd_gradual_transition 
                250.5

For a sharp transition sampled at one centimetre the proxy side contributes a standard deviation of 6.5 years and the chronology 57.1. The date of that event is a chronological quantity almost entirely, and no amount of extra sampling will improve it.

Make the transition gradual, spread over several centimetres and measured with more noise, and the proxy side’s standard deviation rises to 251 years against the same 57.1 from the chronology. Nearly all of that comes from a few runs in which this simple detector, misled by the noise, jumped to a depth near one end of the core; in the other runs the proxy side stays below the chronology. Which term dominates depends on the transition and on how well it is located, not on the core, so it has to be worked out rather than assumed.

An ensemble is the honest output, and it inherits any dating bias

The usual advice is to stop treating the chronology as fixed: draw many chronologies compatible with the dated levels, push the proxy through all of them, and report the spread. That is right, and it is worth knowing how well the resulting interval performs.

set.seed(404)
ens_check <- replicate(300, {
  o  <- draw_ages(age_a(dated))               # the chronology this core actually has
  dd <- detect(0.4, 0.35)
  draws <- replicate(400, approxfun(c(0, dated), c(0, draw_ages(o)), rule = 2)(dd))
  ci <- quantile(draws, c(0.025, 0.975))
  c(covers = age_a(step_depth) >= ci[1] && age_a(step_depth) <= ci[2],
    width  = unname(diff(ci)))
})
round(c(coverage = 100 * mean(ens_check["covers", ]),
        median_width = median(ens_check["width", ])), 1)

    coverage median_width 
        68.7        236.5

The ensemble interval has a median width of 237 years and contains the true date in 69 per cent of cores rather than 95. The shortfall comes from how the dating errors are simulated here rather than from ensembles as such: the two-piece normal used for the dated levels is not centred on zero, so the reported ages are biased old on average, and drawing the ensemble around those reported ages adds the same bias a second time. With a zero-mean dating error the same ensemble comes close to its nominal coverage; an ensemble can only be as well centred as the ages it is built on, and it inherits any bias in them. It is a large improvement on a single chronology, which offers no interval at all.

What to report

Report the age-depth model as a choice, with the software, the interpolation family and the calibration curve named, because a reader cannot reconstruct any of the numbers above without them.

Push the proxy through an ensemble of chronologies rather than one, and report the spread of the quantity you actually care about: the date of the transition, the rate of change, the influx maximum. Propagating uncertainty into the final quantity is not the same as quoting the age errors of the dated levels and leaving the reader to imagine what they do.

Show accumulation rate against depth somewhere. It is the derivative that everything else runs on, it is where the chronology’s structure is visible, and a reader who sees a rate series that jumps at every dated level knows exactly how much of the influx record to believe.

Be careful with influx as a currency. Concentration multiplied by an estimated accumulation rate carries the chronology’s error into what looks like a measurement, and the concentration series itself is often the more honest object.

Honest limits

The dating errors here are a two-piece normal, which gives the right skew and none of the real structure. Calibrated radiocarbon dates have multimodal probability densities, sometimes with two well-separated peaks a century apart, and no smooth interpolation through point estimates represents that. A proper treatment samples from the calibrated densities themselves.

Monotonicity is imposed throughout, both in the truth and in every fitted chronology. Cores with hiatuses, slumps and reworked material break that assumption, and when they do the failure is not a wobble in the influx record but a section of core in the wrong order.

The accumulation model has no autocorrelation. Real sedimentation is persistent, which is what the accumulation rate priors in the Bayesian age-depth packages encode, and a prior of that kind removes some of the manufactured variation measured above. It cannot remove all of it, because the dating error is still there, and adding a prior trades one set of assumptions for another.

Two chronologies were compared and there are many more, including smoothing splines that do not pass through the dated points at all and the Bayesian models that dominate current practice. The comparison here is about what interpolation does to error, not a recommendation between named packages.

Finally, the coverage figure is specific to this design. Six dated levels over three metres is a modest chronology, but adding dates would not cure the shortfall here, because every added date carries the same off-centre error and the ensemble doubles it.

References

Telford RJ, Heegaard E, Birks HJB 2004 Quaternary Science Reviews 23(1-2):1-5 (10.1016/j.quascirev.2003.11.003)

Blaauw M 2010 Quaternary Geochronology 5(5):512-518 (10.1016/j.quageo.2010.01.002)

Blaauw M, Christen JA 2011 Bayesian Analysis 6(3):457-474 (10.1214/11-BA618)

Trachsel M, Telford RJ 2016 The Holocene 27(6):860-869 (10.1177/0959683616675939)

--- title: "Age-depth models and what they do to a proxy" description: "A sediment core with a perfectly constant accumulation rate acquires a fluctuating influx record once dating error is turned into a chronology. Measured in R." date: "2026-08-07 10:00" date-modified: "2026-09-30" categories: [R, palaeoecology, time series, uncertainty, ecology tutorial] image: thumbnail.png image-alt: "Inferred influx along a core whose true accumulation rate is constant, wobbling by a factor of one and a half because the age-depth model bends at each dated level." --- *Corrected 30 September 2026: the post said the ensemble interval undercovers for a structural reason common to all ensembles; it now says the shortfall comes from the simulated dating error not being centred on zero, which the ensemble inherits and doubles, and it no longer says that the ensemble is not what its label suggests or that more dates would improve its centring. It also said a gradual transition reverses the balance between proxy and chronology; it now says that reversal comes from a few runs in which the simple detector jumped to an end of the core.* Every record that comes out of a core is measured against depth and reported against time, and the thing that converts one into the other is an age-depth model fitted to a handful of dated levels. That model is an estimate, built from six or eight ages with errors of a century each, and everything downstream inherits it: accumulation rates, influx, the date of a transition, the rate of change either side of it. The inheritance is easy to miss because the chronology is usually built once, early, by someone else, and then used as though it were a measurement. This post takes a core where the truth is known, puts realistic dating error on the dated levels, and measures what comes out the other end. ## A core with nothing happening in it The cleanest test is a core with a constant accumulation rate. Sediment arrives at the same speed from bottom to top, a proxy arrives at a constant concentration per unit volume, and therefore the influx per unit time is constant too. Anything the analysis reports as variation is manufactured. ```{r setup} #| message: false #| warning: false library(ggplot2) te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), plot.title = element_text(colour = te_ink, face = "bold"), axis.text = element_text(colour = te_body)) } dated <- c(15, 70, 130, 190, 245, 295) # cm, the dated levels sig_lo <- 60 # dating error, years sig_hi <- 110 # longer tail towards older ages conc <- 1200 # grains per cubic cm, constant rate_a <- 0.08 # cm per year, constant age_a <- function(d) d / rate_a # calibrated ages are skewed, so use a two-piece normal, and keep draws in order two_piece <- function(n, s_lo, s_hi) ifelse(runif(n) < s_lo / (s_lo + s_hi), -abs(rnorm(n, 0, s_lo)), abs(rnorm(n, 0, s_hi))) draw_ages <- function(true_age) { repeat { a <- true_age + two_piece(length(true_age), sig_lo, sig_hi) if (all(diff(a) > 0) && a[1] > 0) return(a) } } set.seed(44) obs <- draw_ages(age_a(dated)) round(rbind(depth = dated, true_age = age_a(dated), dated_age = obs)) ``` Six levels, each out by up to about a century, and every one of them within what a laboratory would report. The chronology is then whatever function you choose to run through those points. Linear interpolation and a monotone spline are the two that dominate the literature, and both are one line. ```{r chronologies} slices <- seq(15, 295, by = 10) # the sampled slices lin <- approxfun(c(0, dated), c(0, obs), rule = 2) spl <- splinefun(c(0, dated), c(0, obs), method = "monoH.FC") influx <- function(f) conc * diff(slices) / diff(f(slices)) i_lin <- influx(lin) i_spl <- influx(spl) i_true <- conc * rate_a round(c(true_influx = i_true, linear_cv = sd(i_lin) / mean(i_lin), spline_cv = sd(i_spl) / mean(i_spl), linear_range = max(i_lin) / min(i_lin), spline_range = max(i_spl) / min(i_spl)), 3) ``` The true influx is `r sprintf("%.0f", i_true)` grains per square centimetre per year at every depth in the core. The linear chronology reports a series with a coefficient of variation of `r sprintf("%.3f", sd(i_lin)/mean(i_lin))` and a `r sprintf("%.2f", max(i_lin)/min(i_lin))`-fold range from its lowest slice to its highest. The spline reports `r sprintf("%.3f", sd(i_spl)/mean(i_spl))` and `r sprintf("%.2f", max(i_spl)/min(i_spl))`. ```{r fig-influx} #| fig-cap: "Influx along a core whose true accumulation rate is constant. The horizontal line is the truth. Both chronologies were fitted to the same six dated levels, marked along the top." #| fig-alt: "Two curves of influx against depth wandering above and below a flat horizontal truth line. The linear curve holds constant values between dated levels and jumps at them; the spline curve bends smoothly and overshoots sharply near the top of the core. Six tick marks along the top mark the dated levels." mid <- slices[-length(slices)] + diff(slices) / 2 dd <- rbind(data.frame(depth = mid, influx = i_lin, model = "linear interpolation"), data.frame(depth = mid, influx = i_spl, model = "monotone spline")) ggplot(dd, aes(depth, influx, colour = model, shape = model)) + geom_hline(yintercept = i_true, colour = te_ink, linewidth = 0.7) + geom_line(linewidth = 0.9) + geom_point(size = 2.2) + geom_rug(data = data.frame(depth = dated), aes(x = depth), inherit.aes = FALSE, sides = "t", colour = te_ink, linewidth = 0.8, length = unit(0.035, "npc")) + annotate("text", x = 292, y = i_true, hjust = 1, vjust = -0.7, size = 3.4, colour = te_ink, label = "the truth") + scale_colour_manual(values = c("linear interpolation" = te_rust, "monotone spline" = te_forest)) + scale_shape_manual(values = c(16, 17)) + labs(x = "depth, cm", y = "influx, grains per square cm per year", colour = NULL, shape = NULL, title = "A flat record with structure in it") + theme_datasheet() + theme(legend.position = "top") ``` The shape of the artefact is worth reading off the picture. Linear interpolation holds the accumulation rate constant between dated levels and changes it abruptly at each one, so the influx record is a staircase whose steps sit exactly where the dates are. The spline spreads the same error out into smooth undulations, which look far more like an environmental signal and are for that reason more dangerous, and it overshoots near the top of the core where it has an anchor on one side only. One core is one draw of the dating error, so repeat it. ```{r core-a} cv_and_range <- function(f) { v <- conc * diff(slices) / diff(f(slices)) c(cv = sd(v) / mean(v), range = max(v) / min(v)) } set.seed(31) core_a <- t(replicate(500, { o <- draw_ages(age_a(dated)) c(lin = cv_and_range(approxfun(c(0, dated), c(0, o), rule = 2)), spl = cv_and_range(splinefun(c(0, dated), c(0, o), method = "monoH.FC"))) })) round(rbind(median = apply(core_a, 2, median), upper_decile = apply(core_a, 2, quantile, 0.9)), 3) ``` Across `r sprintf("%d", nrow(core_a))` cores with a genuinely constant accumulation rate, the median manufactured coefficient of variation is `r sprintf("%.3f", median(core_a[, "lin.cv"]))` for linear interpolation and `r sprintf("%.3f", median(core_a[, "spl.cv"]))` for the spline, against a truth of zero. One core in ten reports a range of `r sprintf("%.2f", quantile(core_a[, "spl.range"], 0.9))`-fold or more from the spline. A manufactured swing of `r sprintf("%.0f", 100 * median(core_a[, "spl.cv"]))` per cent in influx is the size of result that gets a paragraph in a results section. ## What it does to a change that is really there Manufacturing variation from nothing is one failure. The other is getting the size of a real change wrong, and that one is harder to notice because the change is genuinely present. ```{r core-b} rate_b <- function(d) ifelse(d > 120 & d < 200, 0.02, 0.12) # a slow section dg <- seq(0, 300, by = 0.25) age_b <- approxfun(dg, cumsum(c(0, diff(dg) / rate_b(dg[-1])))) true_b <- conc * diff(slices) / diff(age_b(slices)) true_amp <- max(true_b) / min(true_b) set.seed(77) core_b <- t(replicate(500, { o <- draw_ages(age_b(dated)) c(lin = cv_and_range(approxfun(c(0, dated), c(0, o), rule = 2))[["range"]], spl = cv_and_range(splinefun(c(0, dated), c(0, o), method = "monoH.FC"))[["range"]]) })) round(c(true_amplitude = true_amp, linear_median = median(core_b[, "lin"]), linear_p90 = quantile(core_b[, "lin"], 0.9), spline_median = median(core_b[, "spl"]), spline_p90 = quantile(core_b[, "spl"], 0.9)), 2) ``` The core now has a real slow section, and the true ratio between the highest and lowest influx across the sampled slices is `r sprintf("%.2f", true_amp)`. Linear interpolation reports a median of `r sprintf("%.2f", median(core_b[, "lin"]))`. The monotone spline reports `r sprintf("%.2f", median(core_b[, "spl"]))`, roughly double the truth, and in the upper decile of cores it reports `r sprintf("%.1f", quantile(core_b[, "spl"], 0.9))`. The direction is not an accident. A spline forced through dated levels that are individually displaced has to bend to reach them, and the bending shows up as very low apparent accumulation somewhere, which turns into a very high apparent influx. The estimator that produces the smoother, better-looking chronology is the one that exaggerates the amplitude most. ## The date of a transition is a chronological quantity Ask when something happened and the answer has two sources of error: where in the core the transition is, and what date the chronology assigns to that depth. It is worth knowing which one dominates, because effort spent sampling the core more finely is wasted if the answer is the other one. ```{r step-date} step_depth <- 160 samp <- seq(20, 290, by = 1) detect <- function(sharpness, noise) { y <- 1 + 2 * plogis((samp - step_depth) / sharpness) + rnorm(length(samp), 0, noise) k <- 6:(length(samp) - 5) gap <- vapply(k, function(i) mean(y[i:length(y)]) - mean(y[1:(i - 1)]), 1) samp[k[which.max(gap)]] } set.seed(88) n_rep <- 600 chron_only <- replicate(n_rep, approxfun(c(0, dated), c(0, draw_ages(age_a(dated))), rule = 2)(step_depth)) sharp_only <- replicate(n_rep, age_a(detect(0.4, 0.35))) soft_only <- replicate(n_rep, age_a(detect(3.5, 0.90))) round(c(true_date = age_a(step_depth), sd_chronology = sd(chron_only), sd_sharp_transition = sd(sharp_only), sd_gradual_transition = sd(soft_only)), 1) ``` For a sharp transition sampled at one centimetre the proxy side contributes a standard deviation of `r sprintf("%.1f", sd(sharp_only))` years and the chronology `r sprintf("%.1f", sd(chron_only))`. The date of that event is a chronological quantity almost entirely, and no amount of extra sampling will improve it. Make the transition gradual, spread over several centimetres and measured with more noise, and the proxy side's standard deviation rises to `r sprintf("%.0f", sd(soft_only))` years against the same `r sprintf("%.1f", sd(chron_only))` from the chronology. Nearly all of that comes from a few runs in which this simple detector, misled by the noise, jumped to a depth near one end of the core; in the other runs the proxy side stays below the chronology. Which term dominates depends on the transition and on how well it is located, not on the core, so it has to be worked out rather than assumed. ## An ensemble is the honest output, and it inherits any dating bias The usual advice is to stop treating the chronology as fixed: draw many chronologies compatible with the dated levels, push the proxy through all of them, and report the spread. That is right, and it is worth knowing how well the resulting interval performs. ```{r ensemble} set.seed(404) ens_check <- replicate(300, { o <- draw_ages(age_a(dated)) # the chronology this core actually has dd <- detect(0.4, 0.35) draws <- replicate(400, approxfun(c(0, dated), c(0, draw_ages(o)), rule = 2)(dd)) ci <- quantile(draws, c(0.025, 0.975)) c(covers = age_a(step_depth) >= ci[1] && age_a(step_depth) <= ci[2], width = unname(diff(ci))) }) round(c(coverage = 100 * mean(ens_check["covers", ]), median_width = median(ens_check["width", ])), 1) ``` The ensemble interval has a median width of `r sprintf("%.0f", median(ens_check["width", ]))` years and contains the true date in `r sprintf("%.0f", 100 * mean(ens_check["covers", ]))` per cent of cores rather than 95. The shortfall comes from how the dating errors are simulated here rather than from ensembles as such: the two-piece normal used for the dated levels is not centred on zero, so the reported ages are biased old on average, and drawing the ensemble around those reported ages adds the same bias a second time. With a zero-mean dating error the same ensemble comes close to its nominal coverage; an ensemble can only be as well centred as the ages it is built on, and it inherits any bias in them. It is a large improvement on a single chronology, which offers no interval at all. ## What to report Report the age-depth model as a choice, with the software, the interpolation family and the calibration curve named, because a reader cannot reconstruct any of the numbers above without them. Push the proxy through an ensemble of chronologies rather than one, and report the spread of the quantity you actually care about: the date of the transition, the rate of change, the influx maximum. Propagating uncertainty into the final quantity is not the same as quoting the age errors of the dated levels and leaving the reader to imagine what they do. Show accumulation rate against depth somewhere. It is the derivative that everything else runs on, it is where the chronology's structure is visible, and a reader who sees a rate series that jumps at every dated level knows exactly how much of the influx record to believe. Be careful with influx as a currency. Concentration multiplied by an estimated accumulation rate carries the chronology's error into what looks like a measurement, and the concentration series itself is often the more honest object. ## Honest limits The dating errors here are a two-piece normal, which gives the right skew and none of the real structure. Calibrated radiocarbon dates have multimodal probability densities, sometimes with two well-separated peaks a century apart, and no smooth interpolation through point estimates represents that. A proper treatment samples from the calibrated densities themselves. Monotonicity is imposed throughout, both in the truth and in every fitted chronology. Cores with hiatuses, slumps and reworked material break that assumption, and when they do the failure is not a wobble in the influx record but a section of core in the wrong order. The accumulation model has no autocorrelation. Real sedimentation is persistent, which is what the accumulation rate priors in the Bayesian age-depth packages encode, and a prior of that kind removes some of the manufactured variation measured above. It cannot remove all of it, because the dating error is still there, and adding a prior trades one set of assumptions for another. Two chronologies were compared and there are many more, including smoothing splines that do not pass through the dated points at all and the Bayesian models that dominate current practice. The comparison here is about what interpolation does to error, not a recommendation between named packages. Finally, the coverage figure is specific to this design. Six dated levels over three metres is a modest chronology, but adding dates would not cure the shortfall here, because every added date carries the same off-centre error and the ensemble doubles it. ## References Telford RJ, Heegaard E, Birks HJB 2004 Quaternary Science Reviews 23(1-2):1-5 (10.1016/j.quascirev.2003.11.003) Blaauw M 2010 Quaternary Geochronology 5(5):512-518 (10.1016/j.quageo.2010.01.002) Blaauw M, Christen JA 2011 Bayesian Analysis 6(3):457-474 (10.1214/11-BA618) Trachsel M, Telford RJ 2016 The Holocene 27(6):860-869 (10.1177/0959683616675939) ## Related tutorials - [Weighted averaging transfer functions in R](../weighted-averaging-transfer-functions/) - [Summed radiocarbon dates: an honest, nearly blind test](../summed-radiocarbon-probabilities/) - [How long a calibration overlap do you need?](../how-long-a-calibration-overlap/) - [Interpolated cores and false slowing down](../interpolated-cores-and-false-slowing-down/)