Basic Usage
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.
Calibration (R_Date)
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 stored successful at /tmp/Rtmp8ISDC3!## 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)## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## ℹ The deprecated feature was likely used in the oxcAAR package.
## Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
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)
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 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"##
## =============================
## 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)## 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"
)
Simulation (R_Simulate)
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 = 4685, std = 25
##
##
## unmodelled: posterior:
##
## one sigma
## 3514 BC - 3492 BC (15.21%)
## 3457 BC - 3376 BC (53.05%)
##
## two sigma
## 3524 BC - 3370 BC (95.45%)
##
## three sigma
## 3623 BC - 3576 BC (2.13%)
## 3568 BC - 3560 BC (0.11%)
## 3532 BC - 3364 BC (97.5%)
##
## Calibrated with:
## Atmospheric data from Reimer et al (2020)
Simulate Sum Calibration
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': 297 obs. of 2 variables:
## ..$ dates : num [1:297] -4350 -4344 -4340 -4334 -4330 ...
## ..$ probabilities: num [1:297] 0.00 0.00 0.00 0.00 4.01e-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"## Warning in max(this_bp_distribution$x): no non-missing arguments to max;
## returning -Inf
Execute custom OxCal code
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 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)"
...## List of 3 calibrated dates:
##
## =============================
## R_Date: Lab-1
## =============================
##
##
## BP = 5000, std = 25
...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.4.4 Bronk Ramsey (2021); r:5"
## .. ..$ comment[1]: chr "Atmospheric data from Reimer et al (2020)"
## .. ..$ comment[2]: chr "( Phase Phase1"
...