## ----setup, include = FALSE, echo=FALSE---------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = NULL
)

## ----packages-----------------------------------------------------------------
# Load packages
library(tabula)
library(magrittr)

## ----mcd-model, fig.width=7, fig.height=3.5, fig.align="center"---------------
## Coerce the zuni dataset to an abundance (count) matrix
zuni_count <- as_count(zuni)

## Set the start and end dates for each ceramic type
zuni_dates <- list(
  LINO = c(600, 875), KIAT = c(850, 950), RED = c(900, 1050), 
  GALL = c(1025, 1125), ESC = c(1050, 1150), PUBW = c(1050, 1150),
  RES = c(1000, 1200), TULA = c(1175, 1300), PINE = c(1275, 1350),
  PUBR = c(1000, 1200), WING = c(1100, 1200), WIPO = c(1125, 1225),
  SJ = c(1200, 1300), LSJ = c(1250, 1300), SPR = c(1250, 1300),
  PINER = c(1275, 1325), HESH = c(1275, 1450), KWAK = c(1275, 1450)
)

## Calculate date midpoint
zuni_mid <- vapply(X = zuni_dates, FUN = mean, FUN.VALUE = numeric(1))
zuni_error <- vapply(X = zuni_dates, FUN = diff, FUN.VALUE = numeric(1)) / 2

## Calculate MCD
## (we use a bootstrapping procedure to estimate the confidence interval)
zuni_mcd <- date_mcd(zuni_count, dates = zuni_mid, errors = zuni_error)
head(zuni_mcd)

## Plot dates
set_dates(zuni_count) <- list(value = zuni_mcd$date, error = zuni_mcd$error)
plot_date(zuni_count, select = 1:15, sort = "asc") +
  ggplot2::labs(title = "Mean Ceramic Date") +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none")

## ----event-model--------------------------------------------------------------
## Bellanger et al. did not publish the data supporting their
## demonstration: no replication of their results is possible. 
## Here is a pseudo-replication using the zuni dataset and results of the 
## previous example.

## Coerce the zuni dataset to an abundance (count) matrix
zuni_count <- as_count(zuni)

## Randomly picks 100 assemblages and assume that they are reliably dated
## (this is NOT a real example)
set.seed(12345)
train_index <- sample(seq_len(nrow(zuni_count)), size = 100, replace = FALSE)
train_dates <- zuni_mcd[train_index, 2]
names(train_dates) <- zuni_mcd[train_index, 1]
set_dates(zuni_count) <- train_dates

## Model the event and accumulation date for each assemblage
model <- date_event(zuni_count, cutoff = 90)
summary(model[["model"]])

## Estimated event dates
head(model[["rows"]])

## ----event-plot, fig.show="hold"----------------------------------------------
## Activity plot
plot_date(model, type = "activity", select = "LZ1105") +
  ggplot2::theme_bw()
## Tempo plot
plot_date(model, type = "tempo", select = "LZ1105") +
  ggplot2::theme_bw()

## ----event-refine-------------------------------------------------------------
## Check model variability
## Jackknife fabrics
refined_jack <- refine_dates(model, method = "jackknife", n = 1000)
head(refined_jack)

## Bootstrap of assemblages
refined_boot <- refine_dates(model, method = "bootstrap", n = 1000)
head(refined_boot)

