Martin Hinz 2021-07-05
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.
library(ggplot2)
library(oxcAAR)
quickSetupOxcal()## Oxcal doesn't seem to be installed. Downloading it now:
## Oxcal stored successful at /tmp/Rtmp0I5Ks0!
## Oxcal path set!
## NULL#setOxcalExecutablePath("~/OxCal/bin/OxCalLinux")
my_date <- oxcalCalibrate(5000,25,"KIA-12345")
my_date##
## =============================
## R_Date: KIA-12345
## =============================
##
##
## BP = 5000, std = 25
##
##
## unmodelled: posterior:
##
## one sigma
## 3890 BC - 3882 BC (4.84%)
## 3796 BC - 3754 BC (36.14%)
## 3745 BC - 3710 BC (27.29%)
##
## two sigma
## 3936 BC - 3872 BC (21.24%)
## 3805 BC - 3702 BC (70.36%)
## 3672 BC - 3656 BC (3.85%)
##
## three sigma
## 3946 BC - 3650 BC (99.73%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)plot(my_date)
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:
##
## =============================
## R_Date: Date 1
## =============================
##
##
## BP = 5000, std = 45
##
##
## unmodelled: posterior:
##
## one sigma
## 3926 BC - 3922 BC (1.25%)
## 3913 BC - 3875 BC (15.11%)
## 3802 BC - 3705 BC (48.04%)
## 3669 BC - 3658 BC (3.87%)
##
## two sigma
## 3944 BC - 3852 BC (28.55%)
## 3846 BC - 3830 BC (2.15%)
## 3818 BC - 3652 BC (64.75%)
##
## three sigma
## 3954 BC - 3644 BC (99.73%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)
##
## =============================
## R_Date: Date 2
## =============================
##
##
## BP = 4500, std = 35
##
##
## unmodelled: posterior:
##
## one sigma
## 3336 BC - 3308 BC (10.7%)
## 3298 BC - 3282 BC (5.55%)
## 3272 BC - 3266 BC (2.14%)
## 3240 BC - 3206 BC (13.46%)
## 3194 BC - 3102 BC (36.41%)
##
## two sigma
## 3355 BC - 3090 BC (94.3%)
## 3049 BC - 3039 BC (1.15%)
##
## three sigma
## 3368 BC - 3012 BC (99.73%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)
##
## =============================
## R_Date: Date 3
## =============================
##
##
## BP = 4000, std = 25
##
##
## unmodelled: posterior:
##
## one sigma
## 2565 BC - 2528 BC (42.42%)
## 2494 BC - 2472 BC (25.85%)
##
## two sigma
## 2572 BC - 2466 BC (95.45%)
##
## three sigma
## 2624 BC - 2454 BC (99.73%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)plot(my_cal_dates)
And you might like to plot them on the calibration curve:
calcurve_plot(my_cal_dates)
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:
str(my_cal_dates, max.level = 1)## List of 3
## $ Date 1:List of 9
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## $ Date 2:List of 9
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## $ Date 3:List of 9
## ..- attr(*, "class")= chr "oxcAARCalibratedDate"
## - attr(*, "class")= chr [1:2] "list" "oxcAARCalibratedDatesList"my_cal_dates[[1]] # equivalent to my_cal_dates[["Date 1"]] or my_cal_dates$`Date 1`##
## =============================
## R_Date: Date 1
## =============================
##
##
## BP = 5000, std = 45
##
##
## unmodelled: posterior:
##
## one sigma
## 3926 BC - 3922 BC (1.25%)
## 3913 BC - 3875 BC (15.11%)
## 3802 BC - 3705 BC (48.04%)
## 3669 BC - 3658 BC (3.87%)
##
## two sigma
## 3944 BC - 3852 BC (28.55%)
## 3846 BC - 3830 BC (2.15%)
## 3818 BC - 3652 BC (64.75%)
##
## three sigma
## 3954 BC - 3644 BC (99.73%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)str(my_cal_dates$`Date 1`)## List of 9
## $ name : chr "Date 1"
## $ type : chr "R_Date"
## $ bp : int 5000
## $ std : int 45
## $ cal_curve :List of 5
## ..$ name : chr "Atmospheric data from Reimer et al (2020)"
## ..$ resolution: num 5
## ..$ bp : num [1:11000] 50100 50095 50090 50086 50081 ...
## ..$ bc : num [1:11000] -53050 -53044 -53040 -53034 -53030 ...
## ..$ sigma : num [1:11000] 1024 1022 1021 1020 1018 ...
## $ sigma_ranges :List of 3
## ..$ one_sigma :'data.frame': 4 obs. of 3 variables:
## .. ..$ start : num [1:4] -3926 -3913 -3802 -3669
## .. ..$ end : num [1:4] -3922 -3875 -3705 -3658
## .. ..$ probability: num [1:4] 1.25 15.11 48.04 3.87
## ..$ two_sigma :'data.frame': 3 obs. of 3 variables:
## .. ..$ start : num [1:3] -3944 -3846 -3818
## .. ..$ end : num [1:3] -3852 -3830 -3652
## .. ..$ probability: num [1:3] 28.55 2.15 64.75
## ..$ three_sigma:'data.frame': 1 obs. of 3 variables:
## .. ..$ start : num -3954
## .. ..$ end : num -3644
## .. ..$ probability: num 99.7
## $ raw_probabilities :'data.frame': 110 obs. of 2 variables:
## ..$ dates : num [1:110] -4060 -4054 -4050 -4044 -4040 ...
## ..$ probabilities: num [1:110] 0.00 0.00 0.00 1.15e-08 8.65e-08 ...
## $ posterior_sigma_ranges :List of 3
## ..$ one_sigma : logi NA
## ..$ two_sigma : logi NA
## ..$ three_sigma: logi NA
## $ posterior_probabilities: logi NA
## - 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##
## =============================
## R_Simulate: SimDate_1
## =============================
##
##
## BP = 4691, std = 25
##
##
## unmodelled: posterior:
##
## one sigma
## 3516 BC - 3494 BC (17.9%)
## 3437 BC - 3377 BC (50.37%)
##
## two sigma
## 3598 BC - 3589 BC (1.04%)
## 3527 BC - 3481 BC (25.47%)
## 3474 BC - 3370 BC (68.95%)
##
## three sigma
## 3624 BC - 3559 BC (3.75%)
## 3532 BC - 3366 BC (95.98%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)plot(my_simulated_dates)
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 calendar 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 9
## $ name : chr " Sum "
## $ type : chr "Sum"
## $ bp : logi NA
## $ std : logi NA
## $ cal_curve :List of 5
## ..$ name : chr "Atmospheric data from Reimer et al (2020)"
## ..$ resolution: num 5
## ..$ bp : num [1:11000] 50100 50095 50090 50086 50081 ...
## ..$ bc : num [1:11000] -53050 -53044 -53040 -53034 -53030 ...
## ..$ sigma : num [1:11000] 1024 1022 1021 1020 1018 ...
## $ sigma_ranges :List of 3
## ..$ one_sigma : logi NA
## ..$ two_sigma : logi NA
## ..$ three_sigma: logi NA
## $ raw_probabilities :'data.frame': 333 obs. of 2 variables:
## ..$ dates : num [1:333] -4344 -4340 -4334 -4330 -4324 ...
## ..$ probabilities: num [1:333] 0.00 0.00 0.00 0.00 3.22e-09 ...
## $ posterior_sigma_ranges :List of 3
## ..$ one_sigma : logi NA
## ..$ two_sigma : logi NA
## ..$ three_sigma: logi NA
## $ posterior_probabilities: logi NA
## - attr(*, "class")= chr "oxcAARCalibratedDate"plot(my_sum_sim)## Warning in max(this_bp_distribution$x): kein nicht-fehlendes Argument für max;
## gebe -Inf zurück
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:
R_Simulate(-4000, 25, "MySimDate")## [1] "R_Simulate(\"MySimDate\",\n -4000, 25);"my_dates <- R_Date(c("Lab-12345","Lab-54321"), c(5000,4500), 25)
cat(my_dates)## R_Date("Lab-12345", 5000, 25);
## R_Date("Lab-54321", 4500, 25);my_sum <- oxcal_Sum(my_dates)
cat(my_sum)## 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:
my_result_data <- parseOxcalOutput(my_result_text)
str(my_result_data)## List of 3
## $ Lab-1:List of 9
## ..$ name : chr "Lab-1"
## ..$ type : chr "R_Date"
## ..$ bp : int 5000
## ..$ std : int 25
## ..$ cal_curve :List of 5
## .. ..$ name : chr "Atmospheric data from Reimer et al (2020)"
...print(my_result_data)## List of 3 calibrated dates:
##
## =============================
## R_Date: Lab-1
## =============================
##
##
## BP = 5000, std = 25
...Or you get the whole output of Oxcal as object:
my_result_data <- parseFullOxcalOutput(my_result_text)
str(my_result_data)## List of 12
## $ ocd[0] :List of 4
## ..$ likelihood:List of 5
## .. ..$ comment :List of 1
## .. .. ..$ : list()
## .. ..$ comment[0]: chr "OxCal v4.4.4 Bronk Ramsey (2021); r:5"
## .. ..$ comment[1]: chr "Atmospheric data from Reimer et al (2020)"
## .. ..$ comment[2]: chr "( Phase Phase1"
...