knitr::opts_chunk$set(fig.width=7, fig.height = 5, fig.align = 'center', warning=FALSE, message=FALSE)
Bchron is an R package that enables quick calibration of radiocarbon dates under various calibration curves (including user generated ones); age-depth modelling as per the algorithm of Haslett and Parnell (2008); Relative sea level rate estimation incorporating time uncertainty in polynomial regression models (Parnell and Gehrels 2015); and non-parametric phase modelling via Gaussian mixtures as a means to determine the activity of a site (and as an alternative to the Oxcal function SUM; currently unpublished).
If you find bugs or want to suggest new features please visit the Bchron GitHub issues page.
Bchron will run in Windows, Mac OS X or Linux. To install Bchron you first need to install R. I would also recommend installing Rstudio as a nice desktop environment for using R. Once in R you can type:
at the R command prompt to install Bchron. If you then type:
it will load in all the Bchron functions.
Bchron will calibrate single or multiple dates under multiple (even user defined) calibration curves. By default, the intcal13, shcal13 and marine13 calibration curves are included. You can calibrate a single radiocarbon date with, e.g.
ages1 = BchronCalibrate(ages=11553, ageSds=230, calCurves='intcal13', ids='Date-1') summary(ages1) plot(ages1)
This will calibrate the radiocarbon age of 11,553 14C years BP with standard error 230 14C years on the intcal13 calibration curve. The
id given is optional and only used for summarising and plotting. The summary command then gives the highest density regions of the calibrated date and the plot command produces a simple plot of the density, together with a shaded region for the 95% highest density region.
Bchron can calibrate multiple dates simultaneously by inputting the dates as vectors:
ages2 = BchronCalibrate(ages=c(3445,11553,7456), ageSds=c(50,230,110), calCurves=c('intcal13','intcal13','shcal13')) summary(ages2) plot(ages2)
This will calibrate three different 14C ages with the calibration curves as specified in the
calCurves argument. The
plot commands will produce individual highest density regions and density plots for the three dates.
Finally, if you provide position information (e.g. depths) to the BchronCalibrate function it will create a plot with position on the y-axis, e.g.:
ages3 = BchronCalibrate(ages=c(3445,11553), ageSds=c(50,230), positions=c(100,150), calCurves=c('intcal13','normal')) summary(ages3) plot(ages3,withPositions=TRUE)
If, alternatively, you want a credible interval (i.e. a single contiguous range) rather than an HDR, you could use the
# First create age samples for each date age_samples = SampleAges(ages3) # Now summarise them with quantile - this gives a 95% credible interval apply(age_samples, 2, quantile, prob=c(0.025,0.975))
SampleAges function introduces some uncertainty as it is a simple random sample. This means that repeatedly running it will give you slightly different results. If you are introduced in precision at a yearly level, you might want to increate the
n_samp argument though will slow the function down a litt.e
If you wanted to find each date's median age (though this is often an inappropriate summary and should be treated with caution), use:
apply(age_samples, 2, quantile, prob=c(0.5))
The calibration code is very fast. On a standard PC you should have no trouble calibrating thousands of dates simultaneously without a noticeably long wait.
The Bchronology function fits the age-depth model outlined by Haslett and Parnell (2008). An illustrative data set from Glendalough is provided with the package, containing 5 radiocarbon dates and a known age for the top of the core. It can be called in via:
The top date is from the present and has the calibration curve 'normal' as it is not a 14C date. This core can be run through Bchron via:
GlenOut = Bchronology(ages=Glendalough$ages, ageSds=Glendalough$ageSds, calCurves=Glendalough$calCurves, positions=Glendalough$position, positionThicknesses=Glendalough$thickness, ids=Glendalough$id, predictPositions=seq(0,1500,by=10))
There are other arguments you can supply to Bchronology, including the date the core was extracted, the outlier probabilities for each individual date, and the number of iterations for which to run the algorithm. For more details see:
Once run, the summary commands will show various output:
summary(GlenOut, type='convergence') summary(GlenOut, type='outliers')
The first summary command produces ages for each position supplied in the
predictPositions argument above (output not shown as it's too long). The second provides convergence diagnostics. The third gives outlier probabilities. The
plot command will produce an age-position plot:
plot(GlenOut, main="Glendalough", xlab='Age (cal years BP)', ylab='Depth (cm)', las=1)
predict command will produce predicted ages for a newly specified set of positions with optional thicknesses:
predictAges = predict(GlenOut, newPositions = c(150,725,1500), newPositionThicknesses=c(5,0,20)) predictAges = predict(GlenOut, newPositions = seq(0,1500,by=10))
To run this model on a data set of your own, you will need to load in your data set via, e.g.
mydata = read.table(file='path/to/file.txt',header=TRUE) run = Bchronology(ages=mydata[,1],ageSds=mydata[,2], ...
The summary function used above also allows for calculation of sedimentation and accumulation rates. Here is an example of calculating and plotting the accumulation rate for the
GlenOut object created above:
acc_rate = summary(GlenOut, type = 'acc_rate')
plot(acc_rate[,'age_grid'], acc_rate[,'50%'], type='l', ylab = 'cm per year', xlab = 'Age (k cal years BP)', ylim = range(acc_rate[,-1])) lines(acc_rate[,'age_grid'], acc_rate[,'2.5%'], lty='dotted') lines(acc_rate[,'age_grid'], acc_rate[,'97.5%'], lty='dotted')
Clearly these accumulation rates are very uncertain, which is unsurprising given the paucity of dates.
To calculate sedimentation rates, you need to either make sure that the
predictPositions argument in the
Bchronology call above is given with unit differences (e.g.
predictPositions=seq(0,1500,by=1)) or set the argument
useExisting=FALSE in the command below:
sed_rate = summary(GlenOut, type = 'sed_rate', useExisting = FALSE)
This will be quite slow as it needs to predict ages on a regular unit depth grid. However, the plots can be created in the same way:
plot(sed_rate[,'position_grid'], sed_rate[,'50%'], type='l', ylab = 'Years per cm', xlab = 'Depth (cm)', ylim = range(sed_rate[,-1])) lines(sed_rate[,'position_grid'], sed_rate[,'2.5%'], lty='dotted') lines(sed_rate[,'position_grid'], sed_rate[,'97.5%'], lty='dotted')
You could then for example save this to a csv file with, e.g.
write.csv(sed_rate, file = 'Glendalough_sed_rates.csv', quote=FALSE, row.names = FALSE)
BchronRSL will produce estimated relative sea level rates from a regression model taking into account the uncertainties in age provided by a
Bchronology run as above. Two example data sets are provided:
These can be run through Bchronology and BchronRSL via:
RSLrun = Bchronology(ages = TestChronData$ages, ageSds = TestChronData$ageSds, positions = TestChronData$position, positionThicknesses = TestChronData$thickness, ids = TestChronData$id, calCurves = TestChronData$calCurves, jitterPositions = TRUE, predictPositions = TestRSLData$Depth) RSLrun2 = BchronRSL(RSLrun, RSLmean = TestRSLData$RSL, RSLsd = TestRSLData$Sigma, degree = 3)
Bchronology run is as described in the section above. The
BChronRSL run takes this object, an estimate of the RSL means and standard deviations, and a value of degree (here 3 indicating cubic regression). They can then be summarised and plotted via:
summary(RSLrun2, type = 'RSL', age_grid = seq(0, 2000, by = 250)) plot(RSLrun2, type = 'RSL', main = 'Relative sea level plot') plot(RSLrun2, type = 'rate', main = 'Rate of RSL change')
See the help files for more options, including outputting parameter values, and plots of acceleration of RSL itself.
Bchron contains two functions for running non-parametric phase models for estimating activity level in a site/region. The first, BchronDensity fits a full Bayesian Gaussian mixture model to the radiocarbon dates whilst the second BchronDensityFast fits an approximate version which will run on much larger data sets. An example run is
data(Sluggan) SlugDens = BchronDensity(ages=Sluggan$ages, ageSds=Sluggan$ageSds, calCurves=Sluggan$calCurves)
You can then output the possible start/end dates of phases:
summary(SlugDens, prob = 0.95)
The probability argument will specify the sensitivity of the algorithm to finding phases. Lower values of
prob will lead to more discovered phases.
You can plot the density with:
plot(SlugDens,xlab='Age (cal years BP)', xaxp=c(0, 16000, 16))
Phases will be shown with a thick black red line. These can be turned off by setting the argument
plotPhase=FALSE whilst the sensitivity can be changed with
BchronDensityFast is identical except for the function call:
SlugDensFast = BchronDensityFast(ages=Sluggan$ages, ageSds=Sluggan$ageSds, calCurves=Sluggan$calCurves)
BchronDensityFast function outputs the age grid and density by default so requires no extra arguments to store it.
Bchron allows for user-defined calibration curves to be created through the function
CreateCalCurve. If inputs to this function are:
nameof the calibration curve (e.g. intcal13c)
Once run, the calibration curve can be used for all future runs until the package is deleted or re-installed. The calibration curve can be used by simply calling it in any other Bchron function using the supplied name. Note that no name checking is done, so it is possible to over-write the pre-defined calibration curves if
As an example, consider installing the intcal09 calibration curve, as available from (http://www.radiocarbon.org/IntCal09%20files/intcal09.14c)[here]. We can run the following code:
# Load in the calibration curve with: intcal09 = read.table('http://www.radiocarbon.org/IntCal09%20files/intcal09.14c',sep=',') # Run CreateCalCurve CreateCalCurve(name='intcal09',cal_ages=intcal09[,1],uncal_ages=intcal09[,2],one_sigma=intcal09[,3])
The new calibration curve can now be used via its name. For example, we can now compare the calibration of a date under the intcal09 curve with that of intcal13:
age_09 = BchronCalibrate(age=15500,ageSds=150,calCurves = 'intcal09',ids='My Date') age_13 = BchronCalibrate(age=15500,ageSds=150,calCurves = 'intcal13') plot(age_09) lines(age_13$Date1$ageGrid,age_13$Date1$densities,col='red') legend('topleft',legend=c('intcal09','intcal13'),col=c('black','red'),lty=1)
For a description of the kind of things that Bchron can do:
Parnell, A. C., Haslett, J., Allen, J. R. M., Buck, C. E., & Huntley, B. (2008). A flexible approach to assessing synchroneity of past events using Bayesian reconstructions of sedimentation history. Quaternary Science Reviews, 27(19-20), 1872–1885.
For the maths behind Bchron:
Haslett, J., & Parnell, A. (2008). A simple monotone process with application to radiocarbon-dated depth chronologies. Journal of the Royal Statistical Society: Series C (Applied Statistics), 57(4), 399–418.
For a review of chronology models:
Parnell, A. C., Buck, C. E., & Doan, T. K. (2011). A review of statistical chronology models for high-resolution, proxy-based Holocene palaeoenvironmental reconstruction. Quaternary Science Reviews, 30(21-22), 2948–2960.
For relative sea level rate estimation:
Parnell, A. C., & Gehrels, W. R. (2015). Using chronological models in late Holocene sea level reconstructions from salt marsh sediments. In: I. Shennan, B.P. Horton, and A.J. Long (eds). Handbook of Sea Level Research. Chichester: Wiley.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.