The chronoverse

Towards a common language for chronological modelling in R

Joe Roe
<https://joeroe.io>

‘Analysis of Dates and Chronological Patterns’ in R

Marwick et al., CTV Archaeology, https://github.com/benmarwick/ctv-archaeology

Radiometric Dating
rcarbon, rintcal, Bchron, clam, nimbleCarbon, rbacon, coffee, oxcAAR, ArchaeoPhases, spDates, IsoplotR, c14, c14bazAAR, UThwigl, ADMUR, ArchaeoChron
Dendrochronology
dendroNetwork, fellingdater
Luminescence Dating
Luminescence, numOSL, BayLum
Paleoenvironmental Proxies
tidypaleo, shoredate
Archaeological Time Series
era, aion, aoristic, kairos, aoristAAR, SPARTAAS, datplot, archSeries, eratosthenes
Stratigraphic Analysis
archeofrag, stratigraphr
Radiocarbon Datasets and APIs
BSDA, c14bazAAR, p3k14c, xronos, rintchron, ArchaeoData

Idiomatic R

  • Vectorisation (everything is a vector)
  • Rectangular data frames
  • Functional programming
  • Function-based OOP (S3, S7)
  • Metaprogramming
  • (Tidy) data as first argument
  • Pipes for function composition (%>%, |>)

📦 c14

📦 c14 <https://c14.joeroe.io>

Tidy radiocarbon data:

  • Converting and calibrating radiocarbon measurements
  • S3 vector for calendar probability distributions (cal, e.g. calibrated radiocarbon dates)
  • Summarising and aggregating calendar probability distributions
  • Cleaning radiocarbon datasets (e.g. c14_control_lab_id())
  • Conversion between types in common radiocarbon packages

Tidy radiocarbon calibration

library(c14)
library(dplyr)
ppnd |>
  mutate(cal = c14_calibrate(cra, error)) |>
  filter(site == "Ganj Dareh") |>
  select(lab_id, cal)
# A tibble: 9 × 2
  lab_id               cal
  <chr>              <cal>
1 GaK 807  c. 12200 cal BP
2 OxA 2099  c. 9900 cal BP
3 OxA 2100 c. 10200 cal BP
# ℹ 6 more rows

Tidy radiocarbon aggregation

ppnd |>
  mutate(cal = c14_calibrate(cra, error)) |>
  summarise(cal_sum = cal_sum(cal), .by = site) # or cal_density()
# A tibble: 83 × 2
  site                cal_sum
  <chr>                 <cal>
1 Tell 'Abr 3 c. 11185 cal BP
2 Abu Gosh    c. 10140 cal BP
3 Abu Hureyra c. 12745 cal BP
# ℹ 80 more rows

cal vectors

vctrs-based S3 class for calendar probability distributions (e.g. calibrated radiocarbon dates)

ganj_dareh <- ppnd[ppnd$site == "Ganj Dareh",]
asiab <- ppnd[ppnd$site == "Asiab",]

ganj_dareh_cal <- c14_calibrate(ganj_dareh$cra, ganj_dareh$error)
asiab_cal <- c14_calibrate(asiab$cra, asiab$error)

c(ganj_dareh_cal, asiab_cal)
<c14_cal[9]>
[1] c. 12200 cal BP  c. 9900 cal BP c. 10200 cal BP c. 10100 cal BP
[5]  c. 9560 cal BP  c. 9540 cal BP  c. 9545 cal BP  c. 9550 cal BP
[9]  c. 9440 cal BP

📦 stratigraphr

📦 stratigraphr <https://stratigraphr.joeroe.io>

Tidy framework archaeological stratigraphy and chronology:

  • Read & write stratigraphic data
  • Graph representation of stratigraphic sequences (Dye & Buck 2015, https://doi.org/10.1016/j.jas.2015.08.008)
  • Validate and clean stratigraphs
  • R interface to OxCal’s Chronological Query Language (CQL)

Import stratigraphic data

library(readr)
read_csv("example_data/harris12.csv")
# A tibble: 10 × 4
  context above below equal
  <chr>   <dbl> <chr> <dbl>
1 1          NA 2,3,4    NA
2 2           1 5        NA
3 3           1 5        NA
# ℹ 7 more rows
library(stratigraphr)
read_lst("example_data/bonn.lst")
# A tibble: 19 × 5
  name  above     contemporary_with equal_to below    
  <chr> <list>    <chr>             <chr>    <list>   
1 +     <chr [1]> <NA>              <NA>     <chr [4]>
2 -     <chr [5]> <NA>              <NA>     <chr [1]>
3 1     <chr [1]> <NA>              <NA>     <chr [1]>
# ℹ 16 more rows

Stratigraphs

h12 <- stratigraph(harris12, "context", "above")
strg_prune(h12) # Harris' 'Law of Succession'
# A tbl_graph: 10 nodes and 12 edges
#
# A directed acyclic simple graph with 1 component
#
# Node Data: 10 × 4 (active)
  context above     below     equal
  <chr>   <list>    <list>    <chr>
1 1       <chr [1]> <chr [3]> <NA> 
2 2       <chr [1]> <chr [1]> <NA> 
3 3       <chr [1]> <chr [1]> <NA> 
# ℹ 7 more rows
#
# Edge Data: 12 × 2
   from    to
  <int> <int>
1     1     2
2     1     3
3     1     4
# ℹ 9 more rows

Plotting ‘Harris matrixes’

E.g. with ggraph <https://ggraph.data-imaginist.com/>:

library(ggraph)
ggraph(h12, layout = "sugiyama") +
  geom_edge_elbow() +
  geom_node_label(aes(label = context), label.r = unit(0, "mm")) +
  theme_graph()

Validating stratigraphs

strat_is_valid(h12)
[1] TRUE
bad_h12 <- tidygraph::bind_edges(h12, data.frame(from = 6, to = 5))
strat_is_valid(bad_h12)
Warning in strat_is_valid(bad_h12): Invalid stratigraphic graph: contains
cycles
[1] FALSE

Bayesian stratigraphic models with stratigraphr + c14

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
 // CQL2 generated by stratigraphr v0.3.0.9000
 Sequence("Shubayqa 1")
 {
  Boundary("");
  Phase("Phase 7")
  {
   R_Date("RTD-7951", 12166, 55);
   R_Date("Beta-112146", 12310, 60);
   R_Date("RTD-7317", 12289, 46);
   R_Date("RTD-7318", 12332, 46);
   R_Date("RTD-7948", 12478, 38);
  };
  Boundary("");
  Phase("Phase 6")
  {
   R_Date("RTD-7947", 12322, 38);
   R_Date("RTD-7313", 12346, 46);
   R_Date("RTD-7311", 12367, 65);
   R_Date("RTD-7312", 12405, 50);
   R_Date("RTD-7314", 12273, 48);
   R_Date("RTD-7316", 12337, 46);
   R_Date("RTD-7315", 12445, 70);
  };
  Boundary("");
  Phase("Phase 5")
  {
   R_Date("RTK-6818", 12477, 76);
   R_Date("RTK-6820", 12385, 75);
   R_Date("RTK-6821", 12385, 78);
   R_Date("RTK-6822", 12412, 79);
   R_Date("RTK-6823", 12321, 78);
  };
  Boundary("");
  Phase("Phase 4")
  {
   R_Date("RTK-6813", 12344, 85);
   R_Date("RTK-6816", 12389, 78);
  };
  Boundary("");
  Phase("Phase 3")
  {
   R_Date("RTK-6819", 11325, 74);
  };
  Boundary("");
  Phase("Phase 2")
  {
   R_Date("RTK-6812", 11365, 72);
   R_Date("RTK-6817", 11322, 75);
  };
  Boundary("");
  Phase("Phase 1")
  {
   R_Date("RTD-8904", 10317, 38);
   R_Date("RTK-6814", 10229, 70);
   R_Date("RTD-8902", 10107, 53);
   R_Date("RTD-8903", 10095, 52);
  };
  Boundary("");
 };

📦 oxcAAR <https://github.com/ISAAKiel/oxcAAR>

# Not executed (too slow)
library(oxcAAR)
quickSetupOxcal(path = fs::path_temp())
shub1_oxcal <- executeOxcalScript(shub1_cql)
readOxcalOutput(shub1_oxcal) |> 
  parseOxcalOutput() |> 
  plot()

📦 era

📦 era <https://era.joeroe.io>

  • S3 vector class for year-based real time scales (yr)
  • Standard for defining eras (e.g. Common Era, Before Present, Hijri)
  • General function for conversion between eras
  • Principled concatenation, coercion, and arithmetic rules
  • Chronology-aware methods (e.g. yr_earliest(), yr_sort())

Years with era

library(era)
years <- yr(c(1100, 1200, 1300), "cal BP")
years
# cal BP years <yr[3]>:
[1] 1100 1200 1300
# Era: Before Present (cal BP): Gregorian years (365.2425 days), counted backwards from 1950
yr_transform(years, "AH")
# AH years <yr[3]>:
[1] 235.4748 132.4059  29.3369
# Era: Anno Hegirae (AH): Islamic lunar years (354.36708 days), counted forwards from 621.53662977337

Era definition

  • label
  • epoch
  • name
  • unit
  • scale
  • direction
era(
    label = "T.A.", 
    epoch = -9021, 
    name = "Third Age", 
    unit = era_year("Gregorian"),
    scale = 1,
    direction = 1
)
<era[1]>
[1] Third Age (T.A.): Gregorian years (365.2425 days), counted forwards from -9021

When it’s working, you don’t notice it

x <- yr(1:3, "BP")
y <- yr(4:6, "cal BP")
z <- yr(7:9, "uncal BP")

y - x
# cal BP years <yr[3]>:
[1] 3 3 3
# Era: Before Present (cal BP): Gregorian years (365.2425 days), counted backwards from 1950
c(y, z)
Error in `vec_ptype2.era_yr.era_yr()`:
! Can't combine `..1` <yr (cal BP)> and `..2` <yr (uncal BP)>.
Reconcile eras with yr_transform() first.
yr_transform(z, yr_era(x)) # No shortcut to calibration!
Error in `yr_transform()`:
! Cannot transform uncalibrated Before Present to Before Present years:
✖ Calendar length of a radiocarbon year is undefined.

📦 tempo

https://github.com/joeroe/tempo

  • Temporal intervals (periods)
  • Temporal relations (Allen + Levy’s extended typology)
library("tempo")
#> Loading required package: rlang

period1 <- c(1500, 1900)
period2 <- c(1800, 1950)

contemporary_with(period1, period2)
[1] TRUE

A ‘chronoverse’?

An ‘ecosystem’ of packages

Packages that work well together or alone:

  • Build up on base R (vector) types and generics
  • Use explicit classes (S3, S4, or S7)
  • Use understandable, predictable and unambiguous function names
  • Have a common and predictable function signature (e.g. .data-first)
  • Avoid complex/idiosyncratic data structures

R package ecosystems

  • tidyverse
  • tidymodels
  • easystats
  • r-spatial / rspatial
  • tesselle

Rudiments of a ‘chronoverse’

c14, stratigraphr, era:
(And tempo, rintchron, xronos, chronologr)

  • vctrs-based formalisms to represent different types of chronological vectors
  • Functional paradigm for composing programs
  • Assume/encourage the use of a tidy data

https://github.com/r-chrono

The chronoverse – towards a common language for chronological modelling in R

Joe Roe <https://joeroe.io>
University of Copenhagen

📦 c14
https://c14.joeroe.io
📦 stratigraphr
https://stratigraphr.joeroe.io
📦 era
https://era.joeroe.io
📦 tempo (v0.0.9000!)
https://github.com/joeroe/tempo
these slides:
https://github.com/joeroe/caa25_chronoverse