Bulk Baconizing

Introduction

Over the course of working with the Neotoma Paleoecology Database (Williams et al., 2018) and the age modeling software Bacon (Blaauw and Christen, 2011) it has often become neccessary to run a large number of Bacon models in an itterative process (for example the Stepps-Baconizing GitHub repository). The process of building and re-building a large number of age models is often complicated and can result in a large amount of code. This code might involve fine-tuning results for single sites, using information from multiple runs, the generation of multiple data files and a set of R code that can span hundreds of lines. In addition, the complexities of the underlying data (for example, issues with existing chronologies and data structures within Neotoma) may mean that users might unintentionally introduce errors into their workflows. Data repositories serve a unique position within the research community (Goring et al., 2018), and one service they can provide is the distillation and distribution of best-practices documents such as this resource for generating a large number of Bacon age models from Neotoma data.

This repository makes use of experience obtained in working through several projects (Dawson et al., 2016; Kujawa et al., 2016), and uses data from several projects as the basis for establishing certain prior values. We make an effort in this document to clearly explain the steps neccessary to undertake such a project in a fairly straightforward manner.

library(rmarkdown)
render('bulk_baconizing.Rmd')

Rscript -e "rmarkdown::render('bulk_baconizing.Rmd')"

Setting up the System

This document is intended to be used as a template for users, and not as a solution in and of itself. The process for generating chronologies is itterative, as such, the use of this Rmd script is intended to be an itterative process, whereby you select sites, run Bacon, revise parameters and run the script again. Each itteration will involve modifying the parameters file (created in the data/params folder during your first run), and also the settings.yaml file. Please be sure to check carefully as you do this. Errors may result in long wait times, or runs that provide no new useful information.

#  Should the program remove any existing files in the Cores directory and start again?
 clean_run: true
 bacon: true
 bchron: true
 
 #  What version of the run is this?  This will create a new parameters file if none exists.
 version: 1
 
 # If a parameters file exists with the current version number should that file overwritten?
 reset_parameters: true
 
 # Where will the Bacon core files go?
 core_path: 'Cores'
 
 #  For any records with assigned dates, what date should be used (either 'YYYY-MM-DD' or 'today')
 date: 'today'
 
 # Core top error (if unspecified)
 core_top_err: 2
 
 # Use modern hiatus break:
 modern: true
 
 # Prior thickness
 thickness: data/paleon_thick.csv
 
 # Settlement indicators
 settlement: data/expert_assessment.csv
 
 # In some work we have assignments for accumulation rates set as defauly.  Others may not.
 accumulation: false
 
 # How many processing cores should be used.
 parallel: 3

source('R/setup_runs.R', echo=FALSE, verbose=FALSE)

Getting Core Data from Neotoma

Here we download the sedimentary data. The key component here is the get_dataset() function. The newest version of neotoma (version >= 1.7.3) allows a user to pass in a vector of dataset IDs, so you can either get datasets using a set of parameters in the get_dataset() function (see help for the get_dataset() function), or you can pass in a vector of dataset IDs (for example, if you’ve already picked a set of sites to be examined). In this example I am using the state, Michigan, as my search parameter:

# Define using state names here.  Could also use altitudes, locations,
#  dataset types. . .

dataset_list <- get_dataset(datasettype='pollen',
                            gpid = 'Utah',
                            ageyoung = 0)

#  The load_downloads() function will download the datasets if either:
#  there is not currently a downloaded object saved with the same dataset version
#  or the setup parameter is TRUE, which calls for the whole data file to be re-written.

downloads <- suppressWarnings(load_downloads(dataset_list,
                            version = settings$version,
                            setup = settings$clean_run))

It might be the case that you want to do further validation on samples based on parameters in the downloads object (e.g., only sites with a certain taxon). The downloads are required for the generation of the core depths file used by Bacon. If you want to further adjust the records, do so with the dataset object, not the downloads. For example, you can remove elements of a list by setting them to NULL, so if you know you want to remove the dataset in the 5th position, you could assign dataset_list[[5]] <- NULL. You would then need to re-run load_downloads() and ensure that setup was TRUE in load_downloads():

# This code is not run:

dataset_list <- get_dataset(gpid = 'Michigan')

# Only accept sites with a site description:

for (i in length(dataset_list):1) {
  if (is.na(dataset_list[[i]]$site$description)) {
    dataset_list[[i]] <- NULL
  }
}

# We would go from approximately 810 datasets to 502.

downloads <- load_downloads(dataset_list,
                            version = settings$version,
                            setup = TRUE)

Using the neotoma R package’s plot_leaflet() function we can generate an interactive map of the dowloaded sites that can be used to explore the data further.

plot_leaflet(dataset_list)

Click on any point to obtain more information about the site. Note that this rendering is only interactive when the document is knit to an HTML file.

Default Bacon Parameters

The Baconizing tool records parameters in a parameter data.frame called params. This data.frame is generated and the saved to file to ensure that errors in the Bacon run do not affect any information stored in the parameters. The data.frame includes a number of columns.

One unique element in this workflow is the use of two accumulation rates, one identified as mod and one identified as old. This comes from work by Goring et al. (2012) and Dawson et al (2016) who identifiy a break in accumulation rates at or around a European settlement horizon. This treatment of modern sediments is the default, however, the modern flag in the settings.yaml file can be set to false if you do not wish to use these modern/historical settings.

Any of these parameters can be changed by altering the values in the code chunk below (in the Rmd file) before running the document. If you have already run the document using a set of parameters and decide to change defaults, you can either change the version in settings.yaml or choose to set reset_parameters to true. Once you have run the full script and you have produced Bacon files you can manually edit individual sites by opening the data/params/bacon_params_version.csv file.

existing_params <- file.exists(paste0('data/params/bacon_params_v', settings$version, '.csv'))

if ((!existing_params) | settings$reset_parameters == TRUE) {

  message("Writing a new parameters file.")

  if (existing_params) {
    file.copy(paste0('data/params/bacon_params_v', settings$version, '.csv'),
              paste0('data/params/bacon_params_v', settings$version, '_', lubridate::round_date(lubridate::now("UTC"), unit="day"),'.csv'))
  }
  ds_handles <- sapply(dataset_list, function(x) { x$dataset.meta$collection.handle })
  ds_ids <- as.integer(sapply(dataset_list, function(x) { x$dataset.meta$dataset.id }))

  params <- data.frame(handle = ds_handles,
                       datasetid = ds_ids,
                       acc.mean.mod = 3.02,
                       acc.mean.old = 15.,
                       acc.shape.mod = 0.53,
                       acc.shape.old = 0.9,
                       mem.strength = 2.,
                       mem.mean = 0.5,
                       hiatus = as.numeric(NA),
                       thick = 5.,
                       age.type = as.character(NA),
                       run = FALSE,
                       suitable = NA,
                       ndates = as.integer(NA),
                       success = NA,
                       reliableold = NA,
                       reliableyoung = NA,
                       notes = ".",
                       stringsAsFactors = FALSE)

  default <- params

  readr::write_csv(x = params,
                   path = paste0('data/params/bacon_params_v', settings$version, '.csv'))
} else {
  params <- suppressMessages(readr::read_csv(paste0('data/params/bacon_params_v', settings$version, '.csv'),
                            col_types = paste0(c('c','i', rep('n',8),'c', 'l','l','i','l','c', 'c','c'), collapse='')))
}

## Writing a new parameters file.

  if(is.null(settings$thickness) | settings$thickness == FALSE | "character" %in% class(settings$thickness)) {

    if ("character" %in% class(settings$thickness)) {
      if (file.exists(settings$thickness)) {
        params <- add_thickness(file = settings$thickness,
                                id_col = 1,
                                thick_col = 4,
                                parameters = params)
          readr::write_csv(x = params,
                     path = paste0('data/params/bacon_params_v', settings$version, '.csv'))

      } else {
        stop("The file defined in `settings.yaml` does not exist.")
      }
    } else {
      message("No thickness file is defined.")
    }
  } else {
    stop("The setting for thickness in `settings.yaml` is not set correctly.")
  }

## None of the dataset ids in the thickness file are in the parameter file.

## Modifying 0 records to update thicknesses.

if(is.null(settings$accumulation) | settings$accumulation == FALSE | "character" %in% class(settings$accumulation)) {

  if ("character" %in% class(settings$accumulation)) {
    if (file.exists(settings$accumulation)) {
      params <- add_accumulation(file = settings$accumulation,
                              id_col = 1,
                              accum_col = 4,
                              parameters = params)
    } else {
      stop("The file defined in `settings.yaml` does not exist.")
    }
  } else {
    message("No accumulation file is defined.")
  }
} else {
  stop("The setting for accumulations in `settings.yaml` is not set correctly.")
}

## No accumulation file is defined.

Generate Core Age and Depth Files

Here we begin to generate the csv files for the cores to allow Bacon to run. This requires calling Neotoma’s API and then resolving the return so that it is turned into a csv. There are several decisions that go into the construction of the CSV file. These decisions are documented both using the notes column of the parameter file, and also commented in the file R/build_agefiles.R. The main decisions are:

1. Given the age model precedence for chronologies in Neotoma (Calendar years BP > Calibrated radiocarbon years BP > Varve years BP > Calendar years BP (Ma) > Radiocarbon years BP), always choose the chronological controls used for the “best” age type. This is the default and no message is raised in the notes field.
2. Given that each age type has a default model, choose the chronologies used in the default model for the best age type. This is the default and no message is raised in the notes field.
3. In cases where a chronology exists for the best age type, but there is no assigned default model, choose the chonology that is most recent, based either on creation date or chronology id. Raises a message in notes
4. In cases where a chronology exists for the best age type, but there are multiple assigned default model, choose the chonology that is most recent, based either on creation date or chronology id. Raises a message in notes

ageorder <- get_table('agetypes')

for (i in 1:nrow(params)) {

  params[i, ] <- build_agefiles(param = params[i, ],
                                ageorder = ageorder,
                                datasets = dataset_list,
                                downloads = downloads,
                                settings = settings)
  readr::write_csv(x = params,
                 path = paste0("data/params/bacon_params_v", settings$version, ".csv"))

}

Possible warnings raised

The script is set up so that places where there is an objective choice, a clear notation of that choice is made and added to the parameter file’s notes column as a semi-colon separated list. In addition to notes regarding the choice of age file, notes are also raised if:

• There are significantly high accumulation rates (more than 100yr/cm)
• There are adjusted indicators for Euro-american settlement (defined in a file reported in settings.yaml.)
• If there is a core top with missing age information (age is assigned -60, to be adjusted by the user).
• If there are no uncertainty constraints on an age that is not associated with geochronological information (uncertainties are assigned to 0).
• If the chronological control is not associated with a geochronological object and comes from an age estimate such as Guess, or Deglaciation (age is dropped)
• If a geochronological object does not have a proper lab ID associated with it.
• If the geochronological object has a NULL age.
• If a 210Pb age has an assigned age greater than 500 (may have been assigned years AD/BC)
• If there was no uncertainty associated with a 210Pb age.
• If there was no uncertainty associated with a geochronological age other than 210Pb.
• If there is only one chronological constraint.
• If core files already existed in the Cores directory associated with the record.
• If for one reason or another something prevented the core file from being written.

Running Bacon

The final step, once all the files have been written, is to run Bacon. Using the rbacon package we can simplify everything:

params <- run_batch(params, settings = settings)

readr::write_csv(x = params,
              path = paste0("data/params/bacon_params_v", settings$version, ".csv"))

The ultimate product here is a Cores directory with a folder for each dataset that generates a suitable chronology for Bacon, a parameters file (in data/params/) that records the parameters used in the Bacon model, along with any notes about elements that may have been changed or were of note in the Bacon runs.

A file has been added to the standard Bacon output. In each directory the file handle_XXX_posteriorout.csv represents draws from the posterior provided by Bacon. You will see NAs in this result (but not in the Bacon ages file)

Run Summary