% 17th July 2019 Version

\documentclass[a4paper,10pt]{article}
\usepackage{psfrag}
\usepackage{graphicx}
\usepackage{afterpage}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{setspace}
\usepackage{subfigure}
\usepackage{enumerate}
\usepackage{natbib}
\usepackage{rotating}
\usepackage{setspace}
\usepackage{lineno}
\usepackage{url}
\usepackage{xcolor}
\usepackage[margin=3cm]{geometry}
\usepackage{color}


\usepackage{bbm}
\usepackage[labelfont=it,textfont={it}]{caption}


\newcommand{\dC}{$\Delta ^{14}C$ }
\newcommand{\rC}{$^{14}C$ }
\newcommand{\rCyrs}{$^{14}C$ yrs }
\newcommand{\fC}{$F^{14}C$ }
\DeclareMathOperator{\EX}{\mathbb{E}}
\providecommand{\keywords}[1]{\textbf{\textit{Keywords:}} #1}

%%%% To mark notes or added text by Maarten (mb)
\newcommand{\mb}{\color{blue}}  %% begin
\newcommand{\bm}{\color{black}} %% end
\newcommand{\diag}{\textrm{diag}}



%% for Paul B
\usepackage{changes}
% and then \added, \deleted, \replaced, \textcolor{red}

\renewcommand{\thepage}{S\arabic{page}} 
\renewcommand{\thesection}{S\arabic{section}}  
\renewcommand{\thetable}{S\arabic{table}}  
\renewcommand{\thefigure}{S\arabic{figure}} 

%opening
\title{Supplementary Information \\
The IntCal20 approach to radiocarbon calibration
curve construction: \\ A new methodology
using Bayesian splines and errors-in-variables}
\author{Heaton TJ \and Blaauw M \and Blackwell PG \and Bronk Ramsey C \and
Reimer PJ \and Scott EM}

\begin{document}

\maketitle
     

\section{Mixing and Convergence}
For the mainly dendrochronologically-dated tine period from approximately 0--14 cal kBP, the curve is predominantly obtained via Gibbs sampling and so convergence is not a significant concern. We therefore concentrate on assessing convergence for the older period of the curve from 14--55 cal kBP whereby the MCMC updating is Metropolis-within-Gibbs. In this older period, all the data have uncertain calendar ages, i.e.\ have errors-in-variables, which are updated during curve construction within the MCMC according to a Metropolis-Hastings (MH) step. It is this calendar age updating step where we might be most concerned about MCMC mixing and convergence, both in terms of ensuring the chain has reached equilibrium and that it has fully explored the potential parameter space to overcome potential multi-modality in the posterior. Addressing these concerns is a significant factor in our change to Bayesian splines from the previous random walk approach. 

In using splines, the MCMC updating is less constrained and so hopefully enables the chain to move around and explore the space more freely. With splines, given the current set of calendar ages, we can update the entire curve simultaneously in a single step. This is a key improvement on the previous random walk approach whereby the curve value was updated a single year at a time conditional upon its value at all other times, with the consequence it was very slow and cumbersome to move. Further, the faster updating possible with splines enable us to perform diagnostics and run multiple chains. This was not computationally feasible with the previous random walk approach. The use of tempering also aims to address convergence concerns, specifically potential multi-modality. 

The large number of parameters in our MCMC mean that it is difficult to assess convergence using standard statistical tests. However, in order to assess whether the target MCMC chain has reached its equilibrium distribution, and whether it fully explores the potential parameter space to overcome potential multi-modality, we performed the following visual checks. 

\subsection{Log-likelihood}
\begin{figure}
  \centering
      \includegraphics[width=0.95\textwidth, keepaspectratio]{MHlikelihood.png}
  \caption{A trace-plot of the log-likelihood of the target MCMC chain, thinned to every $5^{th}$ realisation i.e.\ iteration 10,000 equates to running the MCMC for 50,000 updates. The first 25,000 thinned iterations (i.e. 125,000 MCMC updates pre-thinning) of our sampler are discarded as burn-in. Our curve estimate is based upon the iterations after this burn-in period.  \label{F:MHlikelihood}}
\end{figure}

For the older part of the curve, where we have errors-in-variables, we provide in Figure \ref{F:MHlikelihood} a trace plot of the log-likelihood of the target (non-tempered) MCMC chain. This indicates that the likelihood has reached a settled value by the point we begin to sample from it and hence provides support that hopefully we have reached equilibrium in the chain itself. 

\subsection{Tempering \label{S:Tempering}}


\begin{figure}
  \centering
      \includegraphics[width=0.495\textwidth, keepaspectratio]{TemperedCurveComparison14to16calkBP.png}
      \includegraphics[width=0.495\textwidth, keepaspectratio]{TemperedCurveComparison18to26calkBP.png}
  \caption{Effect of tempering. The left hand plot shows the curves in the region of the floating B\o lling-Aller\o d tree ring sequences from 14--16 cal kBP, the right hand plot shows the curve over a longer period from 18--26 cal kBP.  Note that the data are shown at their prior/observed calendar ages. Consequently assessing their fit to the estimated curve is not appropriate as they are relocated during our Bayesian errors-in-variables curve fitting approach. This is perhaps most apparent in the case of the Cariaco varved sediment (shown in green on the left hand plot) which demonstrates the same shape as the fitted curve yet has a prior calendar age estimate somewhat more recent. The calendar ages of these data will be jointly updated during curve construction to fit with the rest of the data. \label{F:TemperedCurve}}
\end{figure}

\begin{figure}
  \centering
      \includegraphics[width=0.32\textwidth, keepaspectratio]{TemperedTraceLocationA15.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{TemperedTraceLocationP305u.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{TemperedTraceLocationP317.png}
  \caption{Trace plots of estimated location of floating trees for each tempered chain. The first 25,000 plotted values (corresponding to 125,000 MCMC updates pre-thinning) are considered burn-in and not used in our final curve estimate.  \label{F:TemperedFloatingCalEstimates}}
\end{figure} 

Perhaps the most significant concern regarding mixing lies in the potential for calendar age posteriors of the constituent data to be multi-modal, i.e.\ have multiple potential values, with the consequence of a multi-modal posterior for the curve itself. Our MCMC sampler alternates between updating the curve given the current set of calendar age estimates via Gibbs; and then proposing new calendar ages, given the current curve estimate, in a Metropolis-Hastings step. Enabling mixing of calendar age estimates is a particular challenge since, within the Metropolis-Hastings step, the curve which helps determine whether a proposed new set of calendar ages should be accepted or not has been created based upon a fit to the previous set of calendar ages. This means the probability of large changes to calendar age estimates being accepted in any single update is low. Hence, we have the potential for calendar age estimates to get stuck in a single mode of the posterior as switching to another mode would require a significant change to the curve. 

The potential for such an issue is most evident when considering the floating tree-ring sequences. Here, proposed calendar age changes result in moving the complete tree-ring sequence as a whole while rigidly maintaining its internal structure. In the case such a floating sequence has multiple potential posterior calendar age modes, a change from one mode to another in a single update step is challenging since it would be highly probable any such proposed move would be rejected by the MH --- placing the complete tree-ring sequence at the alternate location would be unlikely to fit with the current curve estimate which had been obtained based upon the sequence's previous location.  

Through tempering we hope to reduce this risk. Specifically, we aim to create higher temperature chains in which the calendar ages of the data cam move much more freely within the update steps and cover the full range of potential different modes. In proposing exchanges between these freer chains at high temperatures, which explore more of the space, and those that are less flexible as we reduce the temperatures, we can hopefully enable our sampler to switch between potential modes in our target untempered chain. 

To demonstrate the effect, and benefit, of tempering we present in Figure \ref{F:TemperedCurve} the posterior curve estimates obtained for the target chain (i.e. not tempered) and the chain run at the highest temperature (i.e. increasing up the observational noise). In these plots, the blue curve is the estimate obtained at the highest temperature while the red curve is the target. For interest, we also show, in green, the naive spline estimate obtained were we to ignore the errors-in-variables. As expected, increasing the temperature provides a smoother curve compared with the target. This is particularly evident in the left-hand plot where we might be most concerned about the mixing of the three floating tree-ring sequences which are present in this time period. The higher temperature curve is much smoother and does not show the same significant inversion around 14.75 cal kBP indicating that, at these higher temperatures, the fitted locations of the floating tree-ring sequences are more easily able to move in calendar age. Hence, by proposing state swaps between different temperature chains, we provide the target chain with the potential to switch between multiple modes if it is appropriate. While this is perhaps most important for our tree-ring sequences, it also enables improved mixing for all our uncertain calendar age data.

Relatedly, Figure \ref{F:TemperedFloatingCalEstimates} provides the trace of the calendar age estimates for the three floating B\o lling-Aller\o d tree ring sequences \citep{AdolphiQSR2017} obtained with each tempered chain. Since these trees shown the most significant internal \rC variation, they have the potential to have the most difficult mixing. As described in the main manuscript we increase the amount of tempering on tree P305u and P317 further, in comparison to the other uncertain age data used for IntCal20, to address this. At the higher temperatures (e.g.\ the yellow, most tempered chain) the calendar ages can move relatively freely, seen most clearly in the case of P305u and P317. In the case of P305u, the highest temperature chain moves evenly from approximately 14.35 to 14.45 cal kBP; however as seen in Figure \ref{F:PosteriorFloatCalEstimates}, the untempered target has a bimodal posterior with peaks around 14.36 cal kBP and 14.41 cal kBP. If we were not able to propose swaps from a chain which moved more freely and covered both these modes we might not be able to satisfactorily explore them both in the target but would rather remain solely in one or the other. Conversely, in the case of tree P317, the highest temperature chain appears to alternate between two distinctly different modes, one around 14.5 cal kBP and another around 14.8 cal kBP. As the temperature reduces towards the target, the 14.5 cal kBP mode becomes less likely with probability concentrating around a single mode located around 14.8 cal kBP for the target posterior. The rejection of proposed swaps to the more recent 14.5 cal kBP mode suggests that, for this tree and its quoted \rC uncertainties the older mode around 14.8 cal kBP is the appropriate fitted age.    

\section{Posterior information}

\subsection{Fitting locations of floating tree-ring sequences}
\begin{figure}
  \centering
      \includegraphics[width=0.32\textwidth, keepaspectratio]{PosteriorAgeA15.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{PosteriorAgeP305u.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{PosteriorAgeP317.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{PosteriorAgeFinlayson.png}
      \includegraphics[width=0.32\textwidth, keepaspectratio]{PosteriorAgeMangawhai.png}
  \caption{Posterior calendar age estimates for three B\o lling-Aller\o d floating tree ring sequences \citep{AdolphiQSR2017}; and two older Southern Hemisphere kauri trees \citep{TurneyMangawhai2010, TurneyFinalyson2017}. \label{F:PosteriorFloatCalEstimates} }
\end{figure}


We highlight some of the posterior information provided by the MCMC. Firstly we present, in Figure \ref{F:PosteriorFloatCalEstimates}, the posterior calendar ages of the five floating tree ring sequences --- three from the  B\o lling-Aller\o d \citep{AdolphiQSR2017} and the two older SH kauri trees \citep{TurneyMangawhai2010, TurneyFinalyson2017}. An uninformative prior was placed on the calendar ages of these five sequences and hence the posteriors are determined according to their fit with the rest of the constituent IntCal20 data on which some calendar age information is available. In using such uninformative priors, as opposed to for example the ice-core based calendar age estimate provided by \citet{AdolphiQSR2017}, we maintain independence between the radiocarbon and ice-core timescales. This allows consideration of potential transfer functions between the radiocarbon and ice-core timescales without concerns over circularity. Note in particular the bi-modality for tree P305u indicating there are two potential fitting calendar age locations for this sequence, as already highlighted in Section \ref{S:Tempering}. Since this ring sequence exhibits significant internal \rC variability \citep[see][]{AdolphiQSR2017}, this will also lead to some bi-modality in the posterior for the curve itself. Approximately half the curve realisations will correspond to \rC history where this tree-ring sequence begins around 14.36 cal kBP and other half of the curve realisations will correspond to \rC history where this feature begins around 14.41 cal kBP.  

\subsection{Posterior on smoothing parameters}
\begin{figure}
  \centering
      \includegraphics[width=0.495\textwidth, keepaspectratio]{PosteriorLambdaRecentPeriod.png}
      \includegraphics[width=0.495\textwidth, keepaspectratio]{PosteriorLambdaOldPeriod.png}
  \caption{Posterior estimate of the spline smoothing parameter $\lambda$ in a) the predominantly dendrochronologically-dated based section of the curve from 0--13.9 cal kBP; and b) the older time period from 13.9--55 cal kBP where the curve is based on a wider range of material. Note that due to the differing density of knots in each section, we would not expect the posterior estimates to be the same in each region. \label{F:PosteriorLambda}}
\end{figure}

In Figure \ref{F:PosteriorLambda} we present the posterior estimate of $\lambda$, the spline smoothing parameter, for both sections of the curve separately --- firstly the predominantly dendrochronologically-dated section extending from 0--13,910 cal BP; and secondly the older section of curve from 13,910--55,000 cal BP based upon a wider variety of material including speleothems, lake and marine sediment cores, corals and floating tree-ring sequences. Since the number, and density, of knots in the two separate sections of curve are not equivalent the two posteriors are not directly comparable with one another.  



\section{Choice of additive error model}

\begin{figure}
  \centering
      \includegraphics[width=0.93\textwidth, keepaspectratio]{F14Cboxplots.png}
      \includegraphics[width=0.93\textwidth, keepaspectratio]{ConstanttauSiri.png}
      \includegraphics[width=0.93\textwidth, keepaspectratio]{ScaledtauSiri.png}
      \includegraphics[width=0.93\textwidth, keepaspectratio]{SquareRootScaledtauSiri.png}
  \caption{Comparing different additive error models. An appropriate model would provide the same/similar estimates of $\tau$ when fitted to both the more recent and the older trees taken from SIRI. We model the $F^{14}C$ on tree $i$ in laboratory $j$ as $F_{ij} = f(\theta_i) + \epsilon_{ij} + \zeta_{ij}$. Here $\epsilon_{ij} \sim N(0, \sigma_{ij}^2)$ are independent with $\sigma_{ij}$ the laboratory-reported uncertainty; and $\zeta_{ij}$ denotes the over-dispersion. Panel a) shows box plots of the measured $F^{14}C$ laboratory values for each of the three recent SIRI tree-rings and the two older SIRI tree-rings. Panel b) shows the posterior estimates obtained for $\tau$, on both the recent and older sets of trees separately, if we model the over-dispersion seen within any tree-ring as independent of the value of $f(\theta)$ i.e. $\zeta_{ij} \sim N(0, \tau^2)$. Panel c) if we model the over-dispersion seen within any tree-ring as scaling linearly with the value of $f(\theta)$ i.e. $\zeta_{ij} \sim N(0, \tau^2 f(\theta_i)^2)$. Panel d) if we model it as scaling linearly with $\sqrt{f(\theta)}$ i.e. $\zeta_{ij} \sim N(0, \tau^2 f(\theta_i))$. \label{F:Additive error}}
\end{figure}

We introduce an additive error term to recognise and adapt to potential additional variability in \rC determinations beyond the lab-reported uncertainty. Specifically, in the $F^{14}C$ domain, we model the observed values as
$$
F_i = f(\theta_i) + \epsilon_i + \zeta_i \quad \textrm{for } i = 1, \ldots, n,
$$
where $\theta_i$ is the calendar year of determination $i$, and the function $f(\theta_i)$ is our estimate of the atmospheric level of $F^{14}C$ in that calendar year. Under this model, the observational uncertainty is decomposed into two terms: $\epsilon_i \sim N(0, \sigma_i^2)$, with $\sigma_i$ known, being the independent uncertainty reported by the laboratory; and $\zeta_i \sim N(0, \tau_i^2)$ the additive error which models the potential over-dispersion present in our \rC observations and takes the form of independent random effects with unknown variance $\tau_i^2$. We considered three potential models for the form of $\tau_i$ (the standard deviation of the additive error):
\begin{enumerate}
 \item \textit{Constant additive error model} --- A constant standard deviation model i.e. 
 $$
 \zeta_i \sim N(0, \tau^2).
 $$ Under this model, the level of over-dispersion on $F_i$ is not dependent upon the level of \rC in the object. 
 \item \textit{Scaled (linear) additive error model} --- \label{Scale} A model where the standard deviation of the additive error scales linearly with $f(\theta)$, i.e.\
 $$
 \zeta_{i} \sim N(0, \tau^2 f(\theta_i)^2)
 $$
 \item \textit{Scaled (square root) additive error model} --- \label{SquareRootScale} A model where the standard deviation of the additive error scales linearly with $\sqrt{f(\theta)}$, i.e.\
 $$
 \zeta_{i} \sim N(0, \tau^2 f(\theta_i))
 $$
\end{enumerate}

In order to assess which of these models was most suitable we considered data from the SIRI intercomparison exercise  \citep{SIRI2017} consisting of repeat measurements of trees from two distinct time periods --- 3 trees with radiocarbon determinations of around 300 $^{14}C$yrs BP; and 2 trees with radiocarbon determinations of about 10,000 $^{14}C$yrs BP. Approximately 80 repeat measurements were available on each of the younger group of trees; and 76 repeat measurements on each of the older trees. We implemented separate Bayesian analyses fitting the three above models on each distinct time period independently. Uninformative priors were placed on the values of $f(\theta)$ for each tree and the value of $\tau$. MCMC was then performed for 100,000 iterations, with the first 50,000 discarded as burn-in.      

For a given additive error model, this provided posterior estimates of $\tau$ for both the recent and older sets of trees independently. To assess model suitability, the consistency of these two posterior estimates over the differing time periods were considered. An appropriate model would be expected to provide similar posterior estimates for $\tau$ when fitted to the older trees as when fitted to the more recent trees. The relative estimates are shown in Fig \ref{F:Additive error}. As can be seen, both the model where the standard deviation of the additive error is constant (option 1) and the model it scales linearly with $f(\theta)$ (option 2), provide posterior $\tau$ estimates that are considerably different when fitted to the more recent trees as compared to the older trees. This suggests these additive error models are not appropriate. Conversely, the model whereby the additive error scales linearly with $\sqrt{f(\theta)}$ provides similar estimates for $\tau$ in both the recent and older time periods. It was therefore decided to implement this square-root additive error model (option 3) to adapt to over-dispersion within the IntCal20 data. To provide a prior on the value of $\tau$ to use for curve construction, the SIRI data were combined into a single dataset and another Bayesian analysis performed. The mean and variance of the resultant joint estimate for $\tau$ were calculated and used to determine our prior for IntCal20 as described in the main manuscript.                      

\bibliographystyle{apalike}
%\bibliographystyle{elsarticle-harv}
\bibliography{IntCal20Refs}

\end{document}

