Introduction to the rice package
Introduction
Radiocarbon dating requires a range of calculations, for example
calibration, translations between pMC, F14C, C14 age
and D14C, and assessing the impacts of contamination. This package
provides functions to do so in R.
Installation
On first usage of the package, it has to be installed:
The companion data package ‘rintcal’ which has the radiocarbon
calibration curves will be installed if it isn’t already. New versions
of R packages appear regularly, so please re-issue the above command
regularly to remain up-to-date, or use:
To obtain access to the calibration curves and radiocarbon functions,
first the package has to be loaded:
## Loading required package: rintcal
Calibration curves
The calibration curves can be plotted:
Or, comparing two calibration curves, on the BC/AD scale:
draw.ccurve(1000, 2020, BCAD=TRUE, cc2='marine20', add.yaxis=TRUE)
Or zooming in to between AD 1600 and 2000:
draw.ccurve(1600, 1950, BCAD=TRUE)
Note the wiggles, or recurring C-14 ages? Calibration curves such as
IntCal20 come with 3 columns: cal BP (\(\theta\)), the corresponding C-14 ages
(\(\mu\)), and their errors (\(\sigma\)). Whereas each cal BP or BC/AD
year has a single \(\mu\) and \(\sigma\), owing to wiggles in the
calibration curves, any single \(\mu\)
will often have multiple \(\theta\)’s:
## [,1] [,2]
## [1,] 129 9
## [1] 1694.000 1725.500 1810.667 1873.000 1876.500 1917.000
Let’s add these values to the previous plot:
draw.ccurve(1600, 1950, BCAD=TRUE)
abline(h=BCADtoC14(1694)[1], lty=3)
abline(v=C14toBCAD(129), lty=3)
Atmospheric nuclear tests around AD 1950-1964 caused a huge
human-made C-14 peak (resulting in highly negative C-14 ages):
draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1')
Since the postbomb curve dwarfs the IntCal20 curve, we could also
plot both on separate vertical axes:
draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1', add.yaxis=TRUE)
The calibration curves can also be plotted in the ‘realms’ of F14C,
pMC or D14C (see below), e.g.:
draw.ccurve(50000, 35000, realm="D")
Realms
This package provides functions to translate values between the
radiocarbon-relevant ‘realms’ of cal BP (calendar years before AD 1950),
BC/AD (here we mean cal BC/AD, but shorten this to BC/AD), C14, F14C,
pMC and D14C. The following Table lists the available functions:
| calBP |
- |
calBPtoBCAD |
calBPtoC14 |
calBPtoF14C |
calBPtopMC |
calBPtoD14C |
| BCAD |
BCADtocalBP |
- |
BCADtoC14 |
BCADtoF14C |
BCADtopMC |
BCADtoD14C |
| C14 |
C14tocalBP |
C14toBCAD |
- |
C14toF14C |
C14topMC |
C14toD14C |
| F14C |
NA |
NA |
F14CtoC14 |
- |
F14CtopMC |
F14CtoD14C |
| pMC |
NA |
NA |
pMCtoC14 |
pMCtoF14C |
- |
pMCtoD14C |
| D14C |
NA |
NA |
D14CtoC14 |
D14CtoF14C |
D14CtopMC |
- |
As an example of the above functions, the IntCal20 C14 age and error
belonging to one or more cal BP or cal BCAD ages can be found
(interpolating linearly where necessary):
## [,1] [,2]
## [1,] 168 11
## [,1] [,2]
## [1,] 170 11
## [2,] 174 11
## [3,] 178 11
## [4,] 181 11
## [5,] 185 11
## [6,] 188 11
## [7,] 190 11
## [8,] 193 11
## [9,] 195 11
## [10,] 197 11
## [11,] 199 11
To translate between cal BC/AD and cal BP ages, we can use (the last
example avoids 0 BC/AD, since some calendars do not include zero):
## [1] -74
BCADtocalBP(-1, zero=TRUE)
## [1] 1951
BCADtocalBP(-1, zero=FALSE)
## [1] 1950
D14C (\(\Delta\)14C, a
proxy for atmospheric 14C concentration at t cal BP)
can be transferred to F14C, and the other way around:
## [1] 0.7100983
## [1] 151.8406
This can also be done with C14 ages:
## [,1]
## [1,] 591.9085
## [1] 151.51
These functions can be used to investigate \(\Delta^{14}C\) over time:
x <- seq(0, 55e3, length=1e3)
cc <- calBPtoC14(x)
Dcc <- C14toD14C(cc[,1], cc[,2], x)
par(mar=c(4,3,1,3), bty="l")
plot(x/1e3, Dcc[,1]+Dcc[,2], type="l", xlab="kcal BP", ylab="")
mtext(expression(paste(Delta, ""^{14}, "C")), 2, 1.7)
lines(x/1e3, Dcc[,1]-Dcc[,2])
par(new=TRUE)
plot(x/1e3, (cc[,1]-cc[,2])/1e3, type="l", xaxt="n", yaxt="n", col=4, xlab="", ylab="")
lines(x/1e3, (cc[,1]+cc[,2])/1e3, col=4)
mtext(expression(paste(""^{14}, "C kBP")), 4, 2, col=4)
axis(4, col=4, col.axis=4, col.ticks=4)
Pooling dates
Sometimes, the same material is measured using multiple radiocarbon
dates. If we can be sure that the dated material stems from one single
age in time (e.g., multiple dates on the same single bone, or perhaps
the cereal grains within one bowl which could be assumed to all stem
from the same season), then we can check to which degree the dates agree
using a Chi2-test (Ward & Wilson 1978). If they do agree,
then a pooled mean can be calculated. For example, take the Shroud of
Turin, which was dated multiple times in three different labs:
## ID y er
## 1 AA-3367.1 591 30
## 2 AA-3367.2 690 35
## 3 AA-3367.3 606 41
## 4 AA-3367.4 701 33
## 5 Ox-2575.1 795 65
## 6 Ox-2575.2 730 45
## 7 Ox-2575.3 745 55
## 8 ETH-3883.1 733 61
## 9 ETH-3883.2 722 56
## 10 ETH-3883.3 635 57
## 11 ETH-3883.4 639 45
## 12 ETH-3883.5 679 51
## ! Scatter too large to calculate the pooled mean
## Chisq value: 20.697, p-value 0.037 < 0.05
Zu <- grep("ETH", shroud$ID) # Zurich lab only
pool(shroud$y[Zu],shroud$er[Zu])
## pooled mean: 676.1 +- 23.7
## Chi2: 2.745, p-value 0.601 (>0.05, OK)
## [1] 676.1406 23.7443
Then the weighted mean can be calculated, with the error taken as the
largest of the weighted uncertainty and the standard deviation:
weighted_means(shroud$y[Zu],shroud$er[Zu])
## wmean 676.1, error is max of weighted uncertainty (23.7) & sdev (44.3): 44.3
## [1] 676.1 44.3
If it can indeed be safely assumed that all dates stem from the same
(unknown) calendar year, then the age distribution of that single year
can be plotted together with the individual calibrated ages:
as.one(shroud$y,shroud$er)
## point estimates (mean, median, mode and midpoint): 630, 652, 653 & 618 cal BP
## 95% hpd ranges: 583-570 (29.5%), 667-647 (65.5%)
It would however often be much safer to assume that the multiple
dates were deposited over not just one calendar year but rather over a
period, e.g., over 50 years. To do so, a moving bin is made, and for
each bin placement it is checked how much of the calibrated distribution
of each date fits within that bin. Here is an example, using a bin width
of 100 years, moving at 25 year steps (the top value indicates that
around 660 cal BP, a total of around 4-5 dates fit within the 100-year
bin):
as.bin(shroud$y,shroud$er, 100, move.by=25)
## point estimates (mean, median, mode and midpoint): 641, 627, 625 & 636 cal BP
## 95% hpd ranges: 757-516 (95%)
The spread of multiple calibrated dates can be visualised and
summarised:
spread(shroud$y,shroud$er)
## average spread: 65.6 calendar years (median 56.5)
## 95% range: 2.2 to 187.3
We can assess the overlap between the calibrated distributions of
multiple dates (as the minimum calibrated height for each calendar year
in a range). For example, imagine two dates from the same archaeological
layer, a twig (3820 ± 40 C-14 BP, requiring calibration with IntCal20)
and a marine shell (4430 ± 40 C-14 BP, requiring Marine20, with a
regional marine offset of 90 ± 25). To what degree do their calibrated
distributions overlap?
y <- c(3820, 4590-90)
er <- c(40, sqrt(40^2 + 25^2))
cc <- c(1,2)
overlap(y, er, cc=cc)
## Overlap: 20.1%
Contamination
To calculate the effect of contamination on radiocarbon ages,
e.g. what age would be observed if material with a “true” radiocarbon
age of 5000 ± 20 14C BP would be contaminated with 10% (± 0%)
of modern carbon (F14C=1)?
contaminate(5000, 20, 10, 0, 1)
## True age in F14C: 0.53664 +- 0.00134
## Observed F14C: (0.9*0.53664) + (0.1*1) = 0.58425 +- 0
## Observed C14 age: 4317.1 +- 0
Or imagine that you measured a dinosaur bone, dating to far beyond
the limit of radiocarbon dating, and the sample is very clean as it
contains only 0.5% ± 0.1% modern contamination:
contaminate(66e6, 1e6, 0.5, 0.1, 1)
## True age in F14C: 0 +- 0
## Observed F14C: (0.995*0) + (0.005*1) = 0.005 +- 0.00099
## Observed C14 age: 42719.7 +- 1606.7
The other way round, e.g., inferring what would happen to an observed
age if its assumed 10% modern contamination were to be removed:
## Observed age as F14C: 0.32616 +- 0.00409
## True F14C: (0.9*0.32616) - (0.1*1) = 0.24559 +- 0
## True C14 age: 11279 +- 0
We can also calculate the amount of contamination, or muck, required
to ‘explain away’ certain ages. For example, one of the measurements of
the Shroud of Turin dates to 591 ± 30 C14 BP. How much modern
contamination would have to be inferred for the material to really date
to, say, AD 40?
muck(591, 30, BCADtoC14(40)[,1], 0, 1)
## 1
## Observed age: 591+-30 C14 BP (0.929+-0.003 F14C)
## Target age: 1970 C14 BP (0.783 F14C)
## Calculation: (0.929-0.783)/(1 - 0.783) = 0.67399
## Contamination required: 67.4+-1.6%
So we’d require the sample to have been contaminated by 67% modern
carbon to still date to around AD 40. But what if the sample had been
repaired in, say, AD 1400, thus adding material of an
not-entirely-modern F14C value (i.e., taking into account both the
atmospheric C-14 concentrations in AD 1400 and the fact that some of the
C14 will have decayed since then)?
perFaith <- BCADtoC14(40)
repairF <- BCADtoF14C(1400)
muck(591, 30, perFaith[,1], perFaith[,2], repairF[,1], repairF[,2])
## 0.9317454
## Observed age: 591+-30 C14 BP (0.929+-0.003 F14C)
## Target age: 1970 C14 BP (0.783 F14C)
## Calculation: (0.929-0.783)/(0.932 - 0.783) = 0.98208
## Contamination required: 98.2+-2.4%
This means that the dated sample would have to consist almost
entirely of Medieval age material - which is exactly what was found by
Damon et al. (1989).
The effect of different levels of contamination can be
visualised:
real.14C <- seq(0, 50e3, length=200)
contam <- seq(0, 10, length=101) # 0 to 10% contamination
contam.col <- rainbow(length(contam))
plot(0, type="n", xlim=c(0, 55e3), xlab="real 14C age", ylim=range(real.14C), ylab="observed 14C age")
for (i in 1:length(contam)) {
observed <- contaminate(real.14C, 0, contam[i], 0, 1, decimals=5, talk=FALSE, MC=FALSE)
lines(real.14C, observed[,1], col = contam.col[i])
}
If that is too much code for you, try this function instead:
Fractions
Sometimes, one needs to estimate a missing radiocarbon age from a
sample which has C14 dates on both the entire sample and on fractions,
but where one of the samples was too small to be dated. This can be used
in for example soils separated into size fractions, or samples dated
using both bulk and humic/humin fractions, where one of the samples
turns out to be too small to be dated. This equation requires the bulk
age, the ages of the dated fractions, and the carbon contents and
weights of all fractions.
Cs <- c(.02, .05, .03, .04) # carbon contents of each fraction
wghts <- c(5, 4, 2, .5) # weights for all fractions, e.g., in mg
ages <- c(130, 130, 130, NA) # ages of all fractions. The unmeasured one is NA
errors <- c(10, 12, 10, NA) # errors, unmeasured is NA
fractions(150, 20, Cs, wghts, ages, errors) # assuming a bulk age of 150 +- 20 C14 BP
## estimated C14 age of fraction 4: 519.3 +- 436.7
Calibration
Now on to calibration of radiocarbon dates. We can obtain the
calibrated probability distributions from radiocarbon dates, e.g., one
of 130 ± 10 C14 BP:
calib.130 <- caldist(130, 10, BCAD=TRUE)
plot(calib.130, type="l")
It is also possible to find the likelihood of a single calendar year
for our radiocarbon age, e.g., 145 cal BP:
## [1] 0.0052052
To sample n=10 random values from a calibrated distribution (where
the probability of a year being sampled is proportional to its
calibrated height):
dice <- r.calib(100, 130, 10)
plot(density(dice))
rug(dice)
For reporting purposes, calibrated dates are often reduced to their
95% highest posterior density (hpd) ranges (please report all, not just
your favourite one!):
## from to perc
## [1,] 1710 1684 13.0
## [2,] 1733 1719 8.1
## [3,] 1758 1758 0.1
## [4,] 1824 1803 8.3
## [5,] 1892 1832 50.8
## [6,] 1928 1905 14.8
Additionally, calibrated dates are often reduced to single point
estimates. Note however how poor representations they are of the entire
calibrated distribution!
calib.2450 <- caldist(2450, 20)
plot(calib.2450, type="l")
points.2450 <- point.estimates(calib.2450)
points.2450
## weighted mean median mode midpoint
## 2539.9 2512.2 2666.0 2531.5
abline(v=points.2450, col=1:4, lty=2)
Want a plot of the radiocarbon and calibrated distributions, together
with their hpd ranges?
Sometimes one would want to smooth a calibration curve to take into
account the fact that material has accumulated over a certain time
(e.g., a tree, or peat). To do so, a calibration curve can be smoothed
to produce a tailor-made calibration curve, after which this one is used
to calibrate the date:
mycurve <- smooth.ccurve(smooth=50)
calibrate(2450, 40, thiscurve=mycurve)
Calibrating ‘young’ radiocarbon dates (close to 0 C14 BP) can cause
an error, because a bomb curve might be required to capture the youngest
ages. Do not worry, there is an option to avoid that error:
## Error in calibrate(130, 30) :
## This appears to be a postbomb age (or is close to being one). Please provide a postbomb curve
calibrate(130, 30, bombalert=FALSE)
## Date falls partly beyond calibration curve and will be truncated!
It is also possible to analyse the calibrated probability
distributions, e.g. what is the probability (between 0 and 1) that the
date stems from material that is of the age of 150 cal BP or younger? Or
that it is older than that age?
## [1] 0.7669044
## [1] 0.2330956
Or if you wish to calculate the probability covered between two age
points (e.g., here between 300 and 150 cal BP):
p.range(300, 150, 130, 10)
## [1] 0.2330956
The probabilities reported by p.range should be similar
to those reported by the hpd function mentioned earlier.
Any discrepancies can be explained by rounding errors in the
calculations (e.g., hpd’s settings of every,
age.round and prob.round).
At times there could be a lag in the date, for example if material
has accumulated over decades. We can calibrate a date and then add a
lag, by adding or subtracting either a normal distribution or a gamma
one:
push.normal(2450,40, 400,20)