session 2b

This page contains the R code that will be ran in the second session of the GSA Short Course on age-modeling - part 2b on the practicalities of running rbacon.

The presentation slides of this session can be downloaded from these links: .odp or .pptx.

The first slides which we could reproduce here were done using the R package clam. Assuming you have it installed already (see here), then you simply need to load the code:

The basic, default model is linear interpolation through the dated levels of the standard core that comes with clam:

Bacon

Now we move on to Bacon. If you haven’t installed it yet (or if an updated version might be available for download), do that now:

install.packages('rbacon')

Then load the code:

require(rbacon)

Now run the example core, using default settings (type y or press Enter to accept the defaults):

Bacon()

Output graph of Bacon age-depth model

Let’s run Bacon with thicker, fewer sections (30 cm instead of the 5 cm default; not a good idea, but anyway, in this case do NOT accept the suggested alternative section thickness):

Bacon(, 30)

Output graph of Bacon age-depth model with 30cm thick sections

Note that bad mixing of the MCMC run in the upper left panel. We should run this for longer (using more MCMC iterations, using the ssize parameter), or find different settings that work better.

We could also impose a hiatus at 23 cm depth:

Bacon(, hiatus.depth=23)

Output graph of Bacon age-depth model with a hiatus

Or, imagine that there was a thick tephra layer at 23-20 cm depth:

Bacon(, slump=c(23,20))

Output graph of Bacon age-depth model with a slump

It is also important to play around with the settings for the prior information - not least to see how robust your age-depth model is to different parameter settings:

Bacon(, acc.mean=30)

Output graph of Bacon age-depth model with a non-default prior for accumulation rate

Or set a boundary, and a different prior for accumulation rate below it than the one above it:

Bacon(, boundary=26, acc.mean=c(10, 50))

Output graph of Bacon age-depth model with different prior settings for accumulation rate below and above the boundary

To check if there are any problems with the MCMC run, you can also run your core a few times and compare the MCMC variability between the runs. For example, run the default core 5 times using very few iterations each time:

Baconvergence("MSB2K", 5, ssize=100)

Did 5 Bacon runs.
Gelman and Rubin Reduction Factor 1.11150186669561 (smaller and closer to 1 is better).
Probably not a robust MCMC run! Too much difference between runs, above the 1.05 threshold. Increase sample size?

Output graph of Baconvergence()

The used ssize of just 100 iterations is clearly too small for reliable, stable MCMC runs. In the end you’d want around 3,000+ stored iterations and a stable run, which looks like ‘white noise’ without features such as local optima and long horizontal stretches.

If you have a very long MCMC output and want to thin or cut it to reach a more manageable size, check out the functions scissors() and thinner().

RLGH3 is another core that comes with Bacon. It has fewer dates than MSB2K, and some of them seem to disagree with the others. Accept the suggestions by pressing y and Enter:

Bacon("RLGH3")

Output graph of Bacon age-depth model with core RLGH3

By default, age estimates are calculated for every cm from the top to the bottom core depth. The interval, d.by, can be adapted, and alternatively depths can be provided as a variable, or a file with depths read in (the file should live in the core’s folder, start with the core’s name and end in _depths.txt).

Bacon(d.by=2)
my.depths <- seq(0, 100, by=3)
Bacon(depths=my.depths)
Bacon(depths.file=TRUE)

To see all Bacon’s options, type:

?Bacon

Post-run analyses

You can load a previous, existing run without having to run it again (assuming here that minimum and maximum core depths haven’t changed). First load it and then draw the age-model plot to get all the data in memory:

Bacon("RLGH3", run=FALSE)
agedepth()

Output graph of Bacon age-depth model with core RLGH3

Plot accumulation rates against depth (first setting the layout back to one plot per panel):

layout(1)
accrate.depth.ghost()

A ghost graph of accumulation rate against core depth

or against age:

accrate.age.ghost()

A ghost graph of accumulation rate over time

Additional Bacon tutorials are available here:

You should also be able access the above by typing browseVignettes('rbacon').

prev: session 2a
next: session 3