Introduction

What is bayfoxr?

bayfoxr is a suite of linear Bayesian calibration models for planktic core top foraminiferal δ18O (δ18Oc) and sea surface temperature (SST). These calibrations are especially useful because they capture the uncertainty in the relationship between modern SSTs and coretop δ18Oc. This package is a companion to a paper currently under preparation for the journal “Paleoceanography and Paleoclimatology”.

Citing bayfoxr in your research

Please cite our work if you use bayfoxr in your research. We have a paper currently in preparation and I’ll be sure to update this section with the citation as soon as the paper is out.

To cite the code repository directly use:

Malevich, Steven B., 2019. bayfoxr. <https://github.com/brews/bayfoxr >.

Alternatively, you can cite the package in R’s CRAN repository. You can see this information by running citation("bayfoxr") in an R session.

Installing and loading bayfoxr

You can install bayfoxr from the CRAN repository with:

install.packages("bayfoxr")

Development for bayfoxr is hosted at https://github.com/brews/bayfoxr. You can download and install directly from this github repository with the devtools package:

devtools::install_github("brews/bayfoxr")

This will give you the development versions of the package. It may be unstable. I only recommend this for advanced users.

Once bayfoxr has been installed, you can load it into your R session like any other R package:

library("bayfoxr")

Predicting sea surface temperature and foraminiferal δ18O

bayfoxr revolves around two main functions, predict_d18oc() and predict_seatemp(). Unsurprisingly, these functions help us to predict δ18Oc or SST. Let’s walk through a very simple case with some example sediment core data from Bass River. This data is bundled with bayfoxr, but feel free to use your own data – just be sure it doesn’t contain missing values.

data("bassriver")  # Load the "bassriver" dataframe.

head(bassriver)  # Take a quick look at what the data looks like.
#>    depth  d18o
#> 1 345.17 -3.49
#> 2 345.48 -2.42
#> 3 346.82 -3.03
#> 4 347.06 -3.04
#> 5 349.04 -4.24
#> 6 349.19 -3.37

The example data are marine core samples from John et al. (2008). The dataframe has two columns: “depth”, giving down-core depth in meters, and “d18o”, foraminiferal (Morozovella spp.) calcite δ18O samples (‰ VPDB). The core samples cover the Paleocene-Eocene thermal maximum (PETM).

Now run this information through the predict_seatemp() function:

sst <- predict_seatemp(bassriver$d18o, d18osw = 0.0, prior_mean = 30.0, 
                       prior_std = 20.0)

The predict function spits out a prediction object. Note that we need to specify δ18O for seawater in units ‰ VSMOW (d18osw), and a prior mean and standard deviation for our SST inference.

The sst variable contains an ensemble, or empirical distribution, rather than single prediction points because the calibration is a Bayesian regression model. This ensemble is in sst[['ensemble']]. Here we get median and 90% interval for the prediction:

quants <- quantile(sst, probs = c(0.05, 0.50, 0.95))
head(quants)  # Just see the top of the data...
#>            5%      50%      95%
#> [1,] 27.49985 31.41207 35.34253
#> [2,] 22.89976 26.79676 30.64125
#> [3,] 25.47354 29.38227 33.20536
#> [4,] 25.44748 29.41367 33.26830
#> [5,] 30.85150 34.60419 38.44847
#> [6,] 26.99319 30.84986 34.71672

We can also make a quick and dirty plot to visualize sst:

predictplot(x = bassriver$depth, y = sst, ylim = c(20, 40), 
            ylab = "SST (°C)", xlab = "Depth (m)")

You can see more options for predictplot() with help(predictplot). This applies to any of the functions I’ve mentioned in this tutorial.

Of course, predictions with the predict_d18oc() function are very similar to what we’ve already seen. Here we get a δ18Oc prediction from some made-up SST values. Note that we don’t need to specify a prior mean or standard deviation with predict_d18oc().

d18oc <- predict_d18oc(c(24, 25, 23), d18osw = 0.0)

We can use with this prediction object just like before:

print(quantile(d18oc))
#>             0%       25%       50%       75%       100%
#> [1,] -3.830263 -2.154400 -1.790512 -1.421858  0.6450355
#> [2,] -4.079929 -2.390524 -2.017354 -1.645911 -0.2286378
#> [3,] -3.489937 -1.914891 -1.556574 -1.196472  0.5486787

predictplot(y=d18oc, ylab="δ18Oc (‰ VPDB)", xlab='Sample', ylim = c(-3.5, 0))

The four calibration models

There are four calibration models available for the predict_d18oc() and predict_seatemp() prediction functions.

• Pooled annual
• Hierarchical annual
• Pooled seasonal
• Hierarchical seasonal

The “pooled annual” model is the simplest case, as it has all foraminiferal species pooled together and calibrated against annual SST. The “hierarchical annual” model is similar, but accounts for some species-specific differences in calibration parameters. There are also the “pooled seasonal” and “hierarchical seasonal” models. These are similar to those described above, but use SSTs that are averaged to reflect general temperature preferences of foraminifera. See our 2019 paper for further details on these calibration models and their implementation.

By default, predict_d18oc() and predict_seatemp() use the pooled annual calibration model. To use a seasonal model set the seasonal_seatemp argument to TRUE. To use a hierarchical model, pass a foraminiferal species name string to the foram argument. So, for example, our earlier

d18oc <- predict_d18oc(c(24, 25, 23), d18osw = 0.0)

uses the pooled annual model. To get other variations:

#d18oc <- predict_d18oc(c(24, 25, 23), d18osw = 0.0, seasonal_seatemp = TRUE)  # Pooled seasonal

d18oc <- predict_d18oc(c(24, 25, 23), d18osw = 0.0, foram = "G. bulloides")  # Hierarchical annual for bulloides

# Hierarchical seasonal for bulloides:
d18oc <- predict_d18oc(c(24, 25, 23), d18osw = 0.0, seasonal_seatemp = TRUE,
                       foram = "G. bulloides")

See help(predict_d18oc) or help(predict_seatemp) for more details including a full list of supported foram species.

Generally, the hierarchical annual, hierarchical seasonal, and pooled annual models are recommended. See our 2019 paper for further details