This repository contains sample code for reproducing results in Temporal deconvolution of vascular plant-derived fatty acids exported from terrestial watersheds by Vonk et.al., to appear in Geochimica et Cosmochimica Acta.
The model considered in the paper can be described by the equation ([5] in the paper)
for 
, where 
and 
is a normally distributed (Gaussian) measurement error with known standard deviation. Prior to fitting, the integral is approximated by the sum
where 
and 
. Here, the second approximation uses that 
is essentially zero prior to the bomb-spike.
We set vague prior distributions as follows; 
is uniform on 
, 
is 
distributed, 
is 
and finally 
is apriori uniform on 
(corresponding to an exponential prior for 
with mean 
).
We take a Bayesian approach and fit the model using JAGS called from the R computing environment. Model code can be found in file deconvolve_f.jags and an MCMC sample from the posterior distribution of parameters is drawn using the helper function deconvolve_f. This is illustrated below for the Mackenzie data loaded together with the atmospheric data by load.R:
source("load.R") # Load data and libraries (rjags requires installation of JAGS)
source("functions.R") # Load helper functions
mcmc_sample <- deconvolve_f(D14_Mackenzie, D14_atm) # Draw sample from posteriorThe function draw_tau_hist draws a sample from the posterior predictive distribution of residence times in the younger segment (0-50 years) and plots a histogram:
draw_tau_hist(mcmc_sample, 0:50) + ggtitle("Mackenzie delta, posterior predictive residence times (younger)")
In order to check model fit, plot_fit plots the posterior mean curve together with a shaded 90% credibility interval and data with error bars (four standard deviations wide):
plot_fit(mcmc_sample, D14_Mackenzie, D14_atm) + ggtitle("Mackenzie delta, fitted curve (posterior mean) with data")
Finally table_pars creates a table of posterior parameter summaries:
table_pars(mcmc_sample)| parameter | mean (sd) |
|---|---|
| b | 1.33 (0.13) |
| f_Old | 0.52 (0.02) |
| tau_Old | 28386.23 (9424.77) |
sessionInfo()## R version 3.4.4 (2018-03-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Swedish_Sweden.1252 LC_CTYPE=Swedish_Sweden.1252
## [3] LC_MONETARY=Swedish_Sweden.1252 LC_NUMERIC=C
## [5] LC_TIME=Swedish_Sweden.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 rjags_4-6 coda_0.19-1
## [4] forcats_0.3.0 stringr_1.3.1 dplyr_0.7.6
## [7] purrr_0.2.5 readr_1.1.1 tidyr_0.8.1
## [10] tibble_1.4.2 ggplot2_3.0.0.9000 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.4 haven_1.1.2 lattice_0.20-35 colorspace_1.3-2
## [5] htmltools_0.3.6 yaml_2.2.0 rlang_0.2.2 pillar_1.3.0
## [9] glue_1.3.0 withr_2.1.2 modelr_0.1.2 readxl_1.1.0
## [13] bindr_0.1.1 plyr_1.8.4 munsell_0.5.0 gtable_0.2.0
## [17] cellranger_1.1.0 rvest_0.3.2 codetools_0.2-15 evaluate_0.11
## [21] labeling_0.3 knitr_1.20 highr_0.7 broom_0.5.0
## [25] Rcpp_0.12.19 scales_1.0.0 backports_1.1.2 jsonlite_1.5
## [29] hms_0.4.2 digest_0.6.17 stringi_1.1.7 grid_3.4.4
## [33] rprojroot_1.3-2 cli_1.0.1 tools_3.4.4 magrittr_1.5
## [37] lazyeval_0.2.1 crayon_1.3.4 pkgconfig_2.0.2 xml2_1.2.0
## [41] lubridate_1.7.4 assertthat_0.2.0 rmarkdown_1.10 httr_1.3.1
## [45] rstudioapi_0.8 R6_2.2.2 nlme_3.1-131.1 compiler_3.4.4