Martin Hinz 2018-02-09
The package oxcAAR is designed to represent an interface between R and Oxcal. It offers the possibility to parse Oxcal scripts from data within R, execute them and reread the results from the Oxcal output files. There are other packages that can also calibrate 14C data, like eg. Bchron, which will be sufficient and probably also faster than using oxcAAR. But this package is intended to use especially the algorithms of Oxcal, which is a quasi-standard for archaeological research these days.
Lets assume we want to calibrate the date 5000 BP +- 25 years. oxcAAR has the function oxcalCalibrate for doing so. But at first we have to load the package and tell it where to find the local path to the Oxcal distribution. Afterwards you can calibrate the date using bp, std and name.
## Oxcal doesn't seem to be installed. Downloading it now:
## Oxcal download to /tmp/RtmptWoiBO/Rbuild21626981abaf/oxcAAR/vignettes successful!
## Oxcal path set!
## NULL
#setOxcalExecutablePath("~/OxCal/bin/OxCalLinux")
my_date <- oxcalCalibrate(5000,25,"KIA-12345")
my_date##
## KIA-12345
## BP = 5000, std = 25
##
## one sigma: 3796 BC - 3712 BC
## two sigma: 3931 BC - 3704 BC
## three sigma: 3944 BC - 3692 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
You can also calibrate multiple dates at once:
my_uncal_dates <- data.frame(bp=c(5000,4500,4000),
std=c(45,35,25),
names=c("Date 1", "Date 2", "Date 3")
)
my_cal_dates <- oxcalCalibrate(my_uncal_dates$bp, my_uncal_dates$std, my_uncal_dates$names)
my_cal_dates## List of 3 calibrated dates:
##
## Date 1
## BP = 5000, std = 45
##
## one sigma: 3910 BC - 3706 BC
## two sigma: 3942 BC - 3668 BC
## three sigma: 3954 BC - 3652 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
##
## Date 2
## BP = 4500, std = 35
##
## one sigma: 3334 BC - 3104 BC
## two sigma: 3353 BC - 3089 BC
## three sigma: 3363 BC - 3020 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
##
## Date 3
## BP = 4000, std = 25
##
## one sigma: 2564 BC - 2478 BC
## two sigma: 2572 BC - 2470 BC
## three sigma: 2618 BC - 2460 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
And you might like to plot them on the calibration curve:
The resulting object from the calibration is a list of class oxcAARCalibratedDatesList, containing elements of class oxcAARCalibratedDate. Each of these dates is again a list of the essential informations of the calibrated date including the raw probabilities, which can be extracted for additional analysis:
## List of 3
## $ Date 1:List of 6
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## $ Date 2:List of 6
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## $ Date 3:List of 6
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## - attr(*, "class")= chr [1:2] "list" "oxcAARCalibratedDatesList"
##
## Date 1
## BP = 5000, std = 45
##
## one sigma: 3910 BC - 3706 BC
## two sigma: 3942 BC - 3668 BC
## three sigma: 3954 BC - 3652 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
## List of 6
## $ name : chr "Date 1"
## $ bp : int 5000
## $ std : int 45
## $ cal_curve :List of 5
## ..$ name : chr " IntCal13 atmospheric curve (Reimer et al 2013)"
## ..$ resolution: num 5
## ..$ bp : num [1:10001] 46401 46396 46391 46386 46381 ...
## ..$ bc : num [1:10001] -48050 -48044 -48040 -48034 -48030 ...
## ..$ sigma : num [1:10001] 274 274 274 273 273 ...
## $ sigma_ranges :List of 3
## ..$ one_sigma :'data.frame': 2 obs. of 3 variables:
## .. ..$ start : num [1:2] -3910 -3802
## .. ..$ end : num [1:2] -3876 -3706
## .. ..$ probability: num [1:2] 14.9 53.3
## ..$ two_sigma :'data.frame': 2 obs. of 3 variables:
## .. ..$ start : num [1:2] -3942 -3676
## .. ..$ end : num [1:2] -3692 -3668
## .. ..$ probability: num [1:2] 94.3 1.1
## ..$ three_sigma:'data.frame': 1 obs. of 3 variables:
## .. ..$ start : num -3954
## .. ..$ end : num -3652
## .. ..$ probability: num 99.7
## $ raw_probabilities:'data.frame': 108 obs. of 2 variables:
## ..$ dates : num [1:108] -4054 -4050 -4044 -4040 -4034 ...
## ..$ probabilities: num [1:108] 0.00 0.00 0.00 1.94e-08 1.49e-07 ...
## - attr(*, "class")= chr "oxcAARCalibratedDate"
plot(
my_cal_dates$`Date 1`$raw_probabilities$dates,
my_cal_dates$`Date 1`$raw_probabilities$probabilities,
type = "l",
xlab = "years",
ylab = "probs"
)
You can also use oxcAAR to simulate 14C dates the same way as the OxCal R_Simulate function works. You enter a calibrated year (1000 for 1000 AD, -1000 for 1000 BC) and OxCal will simulate a BP value using a bit of randomisation. Resulting in the fact that each run will have a slightly different BP value.
my_cal_date <- data.frame(bp=c(-3400),
std=c(25),
names=c("SimDate_1")
)
my_simulated_dates <- oxcalSimulate(my_cal_date$bp,
my_cal_date$std,
my_cal_date$names
)
# equivalent to
my_simulated_dates <- oxcalSimulate(-3400, 25, "SimDate_1")
my_simulated_dates##
## SimDate_1
## BP = 4730, std = 25
##
## one sigma: 3627 BC - 3382 BC
## two sigma: 3632 BC - 3376 BC
## three sigma: 3634 BC - 3372 BC
##
## Calibrated after:
## IntCal13 atmospheric curve (Reimer et al 2013)
This package was originally intended to support a series of articles dealing with the investigation of sum calibration. This is why a function is implemented to simulate sum calibration. You can use it to simulate a series of 14C dates and explore the sum calibrated results. You can specify the beginning and end of the time span that should be used for the simulation (in calender years), the number of 14C dates that should be simulated, their standard deviation, either as vector of length n or as one number for all dates, and the type of distribution that should be used (either equally spaced in time, or random uniform). The result is again of class oxcAARCalibratedDate so you can access the raw probabilities for further analysis.
my_sum_sim<-oxcalSumSim(
timeframe_begin = -4000,
timeframe_end = -3000,
n = 50,
stds = 35,
date_distribution = "uniform"
)
str(my_sum_sim)## List of 6
## $ name : chr " Sum "
## $ bp : int(0)
## $ std : int(0)
## $ cal_curve :List of 5
## ..$ name : chr " IntCal13 atmospheric curve (Reimer et al 2013)"
## ..$ resolution: num 5
## ..$ bp : num [1:10001] 46401 46396 46391 46386 46381 ...
## ..$ bc : num [1:10001] -48050 -48044 -48040 -48034 -48030 ...
## ..$ sigma : num [1:10001] 274 274 274 273 273 ...
## $ sigma_ranges :List of 3
## ..$ one_sigma :'data.frame': 0 obs. of 3 variables:
## .. ..$ start : num(0)
## .. ..$ end : num(0)
## .. ..$ probability: num(0)
## ..$ two_sigma :'data.frame': 0 obs. of 3 variables:
## .. ..$ start : num(0)
## .. ..$ end : num(0)
## .. ..$ probability: num(0)
## ..$ three_sigma:'data.frame': 0 obs. of 3 variables:
## .. ..$ start : num(0)
## .. ..$ end : num(0)
## .. ..$ probability: num(0)
## $ raw_probabilities:'data.frame': 295 obs. of 2 variables:
## ..$ dates : num [1:295] -4334 -4330 -4324 -4320 -4314 ...
## ..$ probabilities: num [1:295] 0.00 0.00 0.00 1.76e-09 3.52e-09 ...
## - attr(*, "class")= chr "oxcAARCalibratedDate"
You can also use the package to execute your own OxCal code from within R, and import the results back into the workspace. You can use R_Date, R_Simulate and oxcal_Sum to produce that OxCal code:
## [1] "R_Simulate(\"MySimDate\",\n -4000, 25);"
## R_Date("Lab-12345", 5000, 25);
## R_Date("Lab-54321", 4500, 25);
## Sum(" Sum "){
## R_Date("Lab-12345", 5000, 25);
## R_Date("Lab-54321", 4500, 25);
## };
or use your own script as string variable.
knitr::opts_chunk$set(cache=TRUE)
my_oxcal_code <- ' Plot()
{
Sequence("Sequence1")
{
Boundary("Beginn");
Phase("Phase1")
{
R_Date("Lab-1",5000,25);
R_Date("Lab-2",4900,37);
};
Boundary("Between");
Phase("Phase2")
{
R_Date("Lab-3",4800,43);
};
Boundary("End");
};
};'
my_result_file <- executeOxcalScript(my_oxcal_code)
my_result_text <- readOxcalOutput(my_result_file)You can either parse the result to a ‘standard’ oxcAAR object:
## List of 3
## $ Lab-1:List of 6
## ..$ name : chr "Lab-1"
## ..$ bp : int 5000
## ..$ std : int 25
## ..$ cal_curve :List of 5
## .. ..$ name : chr " IntCal13 atmospheric curve (Reimer et al 2013)"
## .. ..$ resolution: num 5
...
## List of 3 calibrated dates:
##
## Lab-1
## BP = 5000, std = 25
##
## one sigma: 3796 BC - 3712 BC
## two sigma: 3931 BC - 3704 BC
## three sigma: 3944 BC - 3692 BC
...
Or you get the whole output of Oxcal as object:
## List of 12
## $ ocd[0] :List of 4
## ..$ likelihood:List of 5
## .. ..$ comment :List of 1
## .. .. ..$ : list()
## .. ..$ comment[0]: chr "OxCal v4.3.2 Bronk Ramsey (2017); r:5"
## .. ..$ comment[1]: chr " IntCal13 atmospheric curve (Reimer et al 2013)"
## .. ..$ comment[2]: chr "( Phase Phase1"
...