Output graph of Bacon age-depth model
This page contains the R code that was used in the third and final session of the GSA Short Course on age-modeling - part 3 on running Bacon, Plum and your own cores.
The presentation slides of this session can be downloaded from these links: .odp or .pptx.
The .Rmd script is here: session 3.
The streamed session is available here:
video on youtube
chats
transcript
The first slides which we could reproduce here were done using the R package rbacon. Assuming you have it installed already (see here), then you simply need to load the code:
require(rbacon)First we run the default core, using all default settings, in order to look at the prior (green) and posterior (grey) of the accumulation rates:
Bacon()Output graph of Bacon age-depth model
Compare the upper panels with those of the following run, NOT accepting the suggested change in accumulation rate prior (so keeping it at the default of 20 yr/cm):
Bacon("RLGH3", suggest=FALSE)
Output graph of Bacon age-depth model
Now we run the above core again, but accepting the suggested change in acc.rate prior:
Bacon("RLGH3", accept.suggestions=TRUE)
Output graph of Bacon age-depth model
You can also set different acc.rate priors below and above specific core depths, e.g. at 140 cm depth, with a prior of 200 yr/cm below it and 5 yr/cm above it:
Bacon("RLGH3", boundary=140, acc.mean=c(5, 200))
Output graph of Bacon age-depth model
Or model a hiatus instead of a boundary (which is a 0-yr hiatus):
Bacon("RLGH3", hiatus.depths=140, acc.mean=c(5, 200))
Output graph of Bacon age-depth model
Now we move on to the MCMC panel. Let’s do a bad run, with far too few sections, and which will hopefully result in an awful MCMC process:
Bacon(Bacon, 30)
Output graph of Bacon age-depth model
Note that the upper-left panel does not show a nice ‘white-noise’ structure. Instead there is lots of correlation between the iterations, so they are not independent and likely not a good approximation of the true distribution of the accumulation rate parameters.
If you encounter such difficult runs, it is probably a good idea to check for the robustness of the MCMC runs of your core with the applied settings. You can use Baconvergence for this (here we run MSB2K, 5 times, using a very small sample size of only 100):
Baconvergence("MSB2K", 5, ssize=100)
Output graph of Bacon age-depth model
Let’s run some other cores. Amy kindly provided us with two cores, which I ran using the default settings:
Bacon("Swiftcurrent41223")
Output graph of Bacon age-depth model
And also Steel98:
Bacon("Steel98")
Output graph of Bacon age-depth model
You should also be able access the above by typing browseVignettes('rbacon').
Now we move on to Plum, for Pb-210 age-depth modelling. Again we assume you have it installed already (see here), so we simply load it:
require(rplum)Now we run the default core. This takes a while to run, so we might not run the core during the session. rplum runs tend to take longer than rbacon runs, because generally more parameters have to be estimated, not least because the default section thickness for rplum is 1 cm, not rbacon’s 5 (why would you think we made that choice?).
Let’s run the default core (accepting all defaults):
Plum()
Output graph of rplum age-depth model
As with rbacon, we have top panels (but more), and the main panel. All panels are useful and should be presented in papers! The first three are as in rbacon, and then there are two more for the Pb influx and the supported Pb. The main panel also has a difference with rbacon in that it shows the Pb measurements (blue squares) and their modeled values (bluescale).
If your core also has radium data to estimate supported Pb, you can run the core either assuming that the radium values are constant throughout the core (ra.case 1; this requires fewer parameters and thus runs faster):
Plum(“LL14”, ra.case=1)
Output graph of rplum age-depth model
Or assume that the values fluctuate throughout the core (takes longer to run):
Plum(“LL14”, ra.case=2)
Output graph of rplum age-depth model
If you have additional dating information such as historical markers or radiocarbon dates, these can be added as well (using the data format of rbacon’s .csv files):
Plum(“LL14”, ra.case=1, otherdates=”LL14_14C.csv”)
Output graph of rplum age-depth model