Bayesian age modelling with Bacon, an example

This report is an example of the application of a Bayesian age model implemented in R with the software Bacon.

All the files necessary to run this example are hosted here. Download these files into a working directory of your choice. This is what my working directory looks like:

dir()

##  [1] "Bacon.R"              "Bay.Rproj"            "bin"                  "Cores"               
##  [5] "cpp"                  "Curves"               "default_settings.txt" "do.html"             
##  [9] "do.Rmd"               "manualBacon_2.2.pdf"  "README.html"          "README.Rmd"

First we load the software like this:

source("./bacon.R")

## Hi there, welcome to Bacon for Bayesian age-depth modelling

Bacon includes 2 example cores (“MSB2K” and “RLGH3”), and I added one of my cores (“nk1WithOutliers”).

dir("./Cores")

## [1] "MSB2K"           "nk1WithOutliers" "RLGH3"

Let’s have a look at core MSB2K.

head(read.csv("./Cores/MSB2K/MSB2K.csv"))

##       labID  age error depth
## 1 GrA-19478 4128    65   1.5
## 2 GrA-19143 4106    60   4.5
## 3 GrA-19144 4046    59   8.5
## 4 GrA-19146 4184    58  12.5
## 5 GrA-19147 4076    62  14.5
## 6 GrA-19148 4107    61  14.5

Now lets have a look at my core nk1WithOutliers

# What does my core data look like?
head(read.csv("./Cores/nk1WithOutliers/nk1WithOutliers.csv"))

##    name age error depth
## 1 MLB08 -33 2.655    27
## 2 MLE28   0 1.185    45
## 3 MLF45 332 3.500    70
## 4 MLE26 345 3.140    80
## 5 MLF44 379 2.150   110
## 6 MLC04 410 1.750   151

Note that the .csv file containing the data you want to model MUST be named as the folder that contains it and placed into “./Cores”.
You can name your columns as you wish but the order of the columns MUST be .csv “labID”, “age”, “error”, “depth”. This means that you could name you columns: labName, UseriesAge, UseriesAgeError, CoreDepth. But note that the order is the same as in the example above. Also, ages are by default in years BP, errors are 1 sigma and depths are in cm. But all these defaults can be changed and information is in the manual.

You can quickly run Bacon (using all the defaults) on one of the example cores typing:

Bacon("MSB2K")

…and answering y to the prompts.

The defaults need to be changed depending on the nature of your data (and you can learn what is appropriate in the manual). Some of the important parameters that I changed from the defaults are these:

Bacon("nk1WithOutliers", # core
  cc = 0, # U-Th ages already in calendar years, no need calibration
  postbomb=5, # core with postbomb ages (negative BP)
  t.a = 33,  # appropriate for narrow errors of U-series ages
  t.b = 34,  # appropriate for narrow errors of U-series ages
  d.by = 5  # I want ages for every 5 cm section in my core
  )

See the comments next to the code above to understand why I decided to change defaults in a particular way.

Now, you will find new files in your the folder containing your core data

# What files are in my core directory
dir("./Cores/nk1WithOutliers")

## [1] "nk1WithOutliers.csv"          "nk1WithOutliers_71.bacon"     "nk1WithOutliers_71.out"      
## [4] "nk1WithOutliers_71.pdf"       "nk1WithOutliers_71_ages.txt"  "nk1WithOutliers_settings.txt"

In particular, I’m interested in nk1WithOutliers_71_ages.txt

nk1WithOutliers_71_ages <- read.table(
        "./Cores/nk1WithOutliers/nk1WithOutliers_71_ages.txt", header = TRUE)
head(nk1WithOutliers_71_ages)

##   depth   min   max median wmean
## 1    27 -39.1 -28.5  -33.8 -33.7
## 2    32 -35.0 -11.8  -26.0 -25.1
## 3    37 -29.3  -4.1  -16.7 -16.7
## 4    42 -18.2   0.5   -6.9  -7.3
## 5    47   0.1  11.2    4.8   5.1
## 6    52  21.7 183.3  103.4 102.4

nk1WithOutliers_71_ages has all you need!
- depth can be rounded to your needs to match your xrf data (and merege your xrf dataframe with this new age modelled dataframe)

library(plyr)  # install it if you don't have it. It is awesome!!
nk1WithOutliers_71_ages$depth <- round_any(nk1WithOutliers_71_ages$depth, 5,
  floor)

set.seed(1234)  # for reproducible random numbers
# Here I create a fake dataset containing faka Calcium counts (generated
# randomly)
fake_XRF <- data.frame(depth = seq(25, 375, by = 5), 
  fake_Ca = rnorm(71) + 2 * 2)

dat <- join(nk1WithOutliers_71_ages, fake_XRF, type = "full")

## Joining by: depth

head(dat)

##   depth   min   max median wmean  fake_Ca
## 1    25 -39.1 -28.5  -33.8 -33.7 2.792934
## 2    30 -35.0 -11.8  -26.0 -25.1 4.277429
## 3    35 -29.3  -4.1  -16.7 -16.7 5.084441
## 4    40 -18.2   0.5   -6.9  -7.3 1.654302
## 5    45   0.1  11.2    4.8   5.1 4.429125
## 6    50  21.7 183.3  103.4 102.4 4.506056

Now we can plot our fake Calcium varsus age and use the minimun and maximun modelled ages as a measure of uncertainty

library(ggplot2)  # install it if you don't have it. It's awesome!!
ggplot(data = dat,  # this is the dataframe that we just created
  aes(x = median,  # this is the median modelled age
      y = fake_Ca)) +  # this is our face Calcium variable
geom_line(size = 1) +  # line that connets each point
geom_point() +  # add data points to show data density
geom_errorbarh(colour = "grey",  # grey horizontal error bar
  aes(xmin = min,  # minimun modelled age
  xmax = max)) +  # maximun modelled age
scale_x_reverse() + # Reverse x axis so that time progresses from L-R
theme_bw() + # Make plot background white (just because I like it that way)
theme(text = element_text(size = 20))  # Bigger text (just because I like it)

References

Blaauw M, Andres Christen J (2011) Flexible Paleoclimate Age-Depth Models Using an Autoregressive Gamma Process. Bayesian Analysis, 6, 457-474.