In eirikmn/bremla: Bayesian Regression Modeling of Layer-Counted Archives

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  fig.width = 8,
  fig.height = 4,
  dev = "png",
  out.width = "100%"
)

bremla

This repository contains the bremla package (Bayesian Regression Modeling of Layer-counted archives) for uncertainty quantification and synchronization for layer-counted proxy archives. This repository also includes, in the 'inst/' folder, the data and code to reproduce the results for

Myrvoll-Nilsen, E., Riechers, K., Rypdal, M. & Boers, N. (2022). Comprehensive uncertainty estimation of the timing of Greenland warmings of the Greenland Ice core records. Climate of the Past, 18, 1275-1294. doi.org/10.5194/cp-18-1275-2022

and

Myrvoll-Nilsen, E., Riechers, K. & Boers, N. (2025). Synchronization of layer-counted paleoclimatic proxy archives using a Bayesian regression modeling framework. Bayesian Analysis. In review.

Installation

The simplest way to install the package is to install the remotes package and run

#install.packages("remotes")
remotes::install_github("eirikmn/bremla")

To get started, see the documentation '?bremla' and the associated example.

The package includes the NGRIP/GICC05 data set 'NGRIP_d18O_and_dust_5cm.xls' downloaded from Centre for Ice and Climate at the Niels Bohr Institute of the University of Copenhagen, Denmark, and the table of stadial-interstadial events presented by Rasmussen et al. (2014). Also includes the tie-point presented in Adolphi et al. (2018) and Muscheler et al. (2020). Other data and tie-points should work as well.

Warning: The provided real data example generates several thousand simulated chronologies with individual lengths exceeding 18,000. More if both synchronized and unsynchronized chronologies computed. This will require a significant amount of memory. If needed, reduce the number of samples generated or free up memory before use.

Example

This is a real data example for the GICC05 chronology of the NGRIP ice core record which produces 5000 simulations from a time scale synchronized using the Adolphi tie-points.

require(INLA)
library(bremla)
data("event_intervals")
data("events_rasmussen")
data("NGRIP_5cm")

age = NGRIP_5cm$age
depth = NGRIP_5cm$depth
depth2 = depth^2/depth[1]^2 #normalize for stability
d18O = NGRIP_5cm$d18O

data = data.frame(age=age,dy=c(NA,diff(age)),depth=depth,depth2=depth2,d18O=d18O)
formula = dy~-1+depth2+d18O

events=list(locations = events_rasmussen$depth,
            locations_unit="depth",degree=1)

results = bremla(formula,data,reference.label="GICC05",
                nsims=5000,
                events=events,
                synchronization=list(method="adolphi"),
                control.fit=list(method="inla"),
                control.sim=list(synchronized=TRUE) )
summary(results)

require(ggplot2)
free_indexes = results$tie_points$free_indexes
tie_indexes = results$tie_points$tie_indexes
ageref_free = results$data$age[free_indexes]
ageref_tie = results$data$age[tie_indexes]

fullpd = data.frame(depth = results$data$depth[free_indexes],
                medians=results$simulation$summary_sync$median[free_indexes]-ageref_free,
                lower=results$simulation$summary_sync$lower[free_indexes]-ageref_free,
                upper=results$simulation$summary_sync$upper[free_indexes]-ageref_free)

gg2 = ggplot(data=fullpd,aes(x=.data$depth)) +
geom_line(aes(y=0),color="blue",linetype="dotted",size=0.2)+
theme_bw()+ylab("Estimated age - reference (years)")+
xlab("NGRIP depth (m)")


gg2 = gg2+geom_line(aes(y=.data$medians))+
geom_ribbon(aes(ymin=.data$lower,ymax=.data$upper),color="red",fill="red",alpha=0.3)


gg2 = gg2 + geom_segment(data=data.frame(depth=results$data$depth[tie_indexes],                                   median=results$simulation$summary_sync$median[tie_indexes]-ageref_tie,                         lower=results$simulation$summary_sync$lower[tie_indexes]-ageref_tie,
                   upper=results$simulation$summary_sync$upper[tie_indexes]-ageref_tie),
                   aes(x=.data$depth,y=.data$lower,xend=.data$depth,yend=.data$upper),
                   col="magenta")
print(gg2)

Attribution

This code is associated and written for the papers Myrvoll-Nilsen et al. (2022) and Myrvoll-Nilsen et al. (202x) mentioned above. Feel free to use the code, but please cite the accompanying papers.