-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
112 lines (83 loc) · 5.21 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
fig.width = 8,
fig.height = 4,
dev = "png",
out.width = "100%"
)
```
# bremla
<!-- badges: start -->
[](https://github.com/eirikmn/bremla/actions/workflows/R-CMD-check.yaml)
[](https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html)
<!-- badges: end -->
This repository contains the bremla package (**B**ayesian **Re**gression **M**odeling of **L**ayer-counted **a**rchives) for uncertainty quantification and synchronization for layer-counted proxy archives, including 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. (202x). Tie-point paper (title TBD)
## Installation
The simplest way to install the package is to install the remotes package and run
```{r install, eval=FALSE}
#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](https://www.iceandclimate.nbi.ku.dk) at the Niels Bohr Institute of the University of Copenhagen, Denmark, and the table of stadial-interstadial events presented by [Rasmussen et al. (2014)](https://www.sciencedirect.com/science/article/pii/S0277379114003485). Also includes the tie-point presented in [Adolphi et al. (2018)](https://cp.copernicus.org/articles/14/1755/2018/) and [Muscheler et al. (2020)](https://www.cambridge.org/core/journals/radiocarbon/article/testing-and-improving-the-intcal20-calibration-curve-with-independent-records/D72D9214C47FE9441B5E730D33DCCE3D). 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.
```{r example, message=FALSE, warning=FALSE, include=TRUE}
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)
```
```{r plot, message=FALSE, eval=TRUE, echo=FALSE}
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.
## License
The code in this repository is made available under the terms of the GNU (version >=2) License. For details, see LICENSE.md file.