In joeroe/stratigraphr: Archaeological Stratigraphy and Chronological Sequences

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

The Chronological Query Language (CQL) is a tool for formally describing chronological models [@Bronk_Ramsey1998-zw]. It is most commonly used to input data for Bayesian radiocarbon calibration in OxCal [@Bronk_Ramsey2009-yk]. stratigraphr includes an R interface for CQL2, the version used in OxCal v3+.

This vignette describes how to use this interface to generate CQL models in R.

Overview

• Individual CQL commands (see https://c14.arch.ox.ac.uk/oxcalhelp/hlp_commands.html) are implemented as cql_* functions.
• Use cql() to group together CQL functions or include arbitrary CQL code.
• Execute CQL scripts in OxCal using write_oxcal() or the oxcAAR package.

Writing CQL models

library("stratigraphr")

Generating CQL models from other types of data

Used in the simple way above, stratigraphr's CQL interface offers little benefit over writing CQL directly. Its real power is in combining cql() with other R tools to build models based on other data.

For the following examples, we will use stratigraphic and radiocarbon data from Shubayqa 1 [@Richter2017-xy]:

data("shub1")
data("shub1_radiocarbon")

Tabular data

With radiocarbon and stratigraphic data in a tabular format, you can take advantage of dplyr's powerful tools for data manipulation to build CQL models programmatically.

For example, use dplyr::mutate() to concisely express a table of dates as CQL R_Date commands:

library("purrr")
library("dplyr")

shub1_radiocarbon %>% 
  mutate(r_date = cql_r_date(lab_id, cra, error)) %>% 
  pluck("r_date") %>% 
  cql()

Or use dplyr::group_by() and dplyr::summarise() to build phase models. This is a three stage process, illustrated below with the phase model from Shubayqa 1:

1. Preprocess the dates, e.g. to remove outliers and order phases.
2. Group dates by phase and summarise them with cql_phase().
3. Group phases into a sequence, using boundaries to automatically add boundary constraints between them. If we had multiple independent sequences, we could include an additional grouping variable here, before concatenating them with cql().

shub1_radiocarbon %>% 
  filter(!outlier) %>% 
  group_by(phase) %>% 
  summarise(cql = cql_phase(phase, cql_r_date(lab_id, cra, error))) %>% 
  arrange(desc(phase)) %>% 
  summarise(cql = cql_sequence("Shubayqa 1", cql, boundaries = TRUE)) %>%
  pluck("cql") %>%
  cql() ->
  shub1_cql

shub1_cql

Stratigraphs

...

Running CQL models

You can run models generated by cql() using the desktop or online versions of OxCal by simply copying the output into the program. Alternatively, use write_oxcal() to create a .oxcal file:

oxcal_cql <- cql(
  cql_r_date("ABC-001", 9100, 30),
  cql_r_date("ABC-002", 9200, 30),
  cql_r_date("ABC-003", 9300, 30)
)

write_oxcal(oxcal_cql, "cql.oxcal")

file.remove("cql.oxcal")

You can also run OxCal directly through R using the oxcAAR package. This depends on a local installation of OxCal. If you already have one installed, you can set the path to the executable using oxcAAR::setOxcalExecutablePath(). Otherwise, use oxcAAR::quickSetupOxcal() to download one, for example to a temporary directory:

library("oxcAAR")
quickSetupOxcal(path = tempdir())

You can then use oxcAAR::executeOxcalScript() to run the CQL script and oxcAAR::readOxcalOutput() to read the output back into R.

executeOxcalScript(oxcal_cql) %>% 
  readOxcalOutput() ->
  oxcal_output

You can parse the output with oxcAAR::parseOxcalOutput() and visualise it using oxcAAR's built-in plotting functions:

oxcal_parsed <- oxcAAR::parseOxcalOutput(oxcal_output)
plot(oxcal_parsed)
calcurve_plot(oxcal_parsed)

The current CRAN version of oxcAAR (v. 1.0.0) does not read the posterior probabilities produced by a model with Bayesian calibration, so to work with these you need to install the latest development version (devtools::install_github("ISAAKiel/oxcAAR")). With this, oxcAAR::parseOxcalOutput() also contains the modelled results in $posterior_sigma_ranges and $posterior_probabilities. Again, you can quickly visualise these with the built-in plotting functions:

# Not run: slow
# shub1_oxcal <- executeOxcalScript(shub1_cql)
# readOxcalOutput(shub1_oxcal) %>% 
#   parseOxcalOutput() %>% 
#   plot()

References

joeroe/stratigraphr documentation built on May 17, 2023, 9:52 p.m.