Welcome to bayclumpr! Before we get started with this
tutorial, we would like to remind you that there is an associated shiny
app that accompanies this R package. You can access
BayClump directly from your browser using by clicking here. Now, let’s
go ahead and discuss some of the basic functions in
bayclumpr.
Performing calibrations using bayclumpr
First, we will need some data to work with. We can use
bayclumpr to generate simulated datasets with uncertainty
values described in Roman-Palacios et al. (2022). For this example, we
will simulate 50 observations under a low-eror scenario. Note
that the functions in bayclumpr expect users to provide
uncertainty in terms of standard deviation. The resulting
dataset will be stored in the ds object.
ds <- cal.dataset(error = "S1", nobs = 50)
head(ds)
#> x_TRUE Temperature TempError y_TRUE D47error D47 Material
#> 1 10.076133 10.089943 0.013809939 0.6475262 -0.0032283733 0.6442978 1
#> 2 11.414949 11.399569 -0.015379377 0.6841481 0.0065876134 0.6907357 1
#> 3 12.517576 12.522651 0.005074617 0.7430624 0.0012176807 0.7442800 1
#> 4 9.695736 9.662728 -0.033008011 0.6333012 0.0021347308 0.6354360 1
#> 5 12.391566 12.364749 -0.026817078 0.7379670 0.0027211068 0.7406881 1
#> 6 12.060248 12.051630 -0.008617473 0.7206252 0.0005650349 0.7211903 1Now, let’s start by fitting different models in the simulated
dataset. For instance, let’s fit a Deming regression model using the
cal.deming function in bayclumpr:
cal.deming(data = ds, replicates = 10)
#> alpha beta
#> 1 0.2440497 0.03906877
#> 2 0.2614754 0.03670698
#> 3 0.2760968 0.03574009
#> 4 0.2419661 0.03857797
#> 5 0.2752911 0.03585697
#> 6 0.2526154 0.03764163
#> 7 0.2497403 0.03784485
#> 8 0.2505127 0.03785190
#> 9 0.2590315 0.03708834
#> 10 0.2057244 0.04158558Alternatively, you can fit an unweighted or weighted OLS regression
using cal.ols and cal.wols functions,
respectively:
cal.ols(data = ds, replicates = 10)
#> alpha beta
#> 1 0.2839352 0.03554307
#> 2 0.2872762 0.03527608
#> 3 0.2733101 0.03642586
#> 4 0.2640842 0.03709674
#> 5 0.2652946 0.03701428
#> 6 0.2668595 0.03679151
#> 7 0.2660233 0.03687870
#> 8 0.2826984 0.03542715
#> 9 0.2693398 0.03663859
#> 10 0.2850873 0.03507551
cal.wols(data = ds, replicates = 10)
#> alpha beta
#> 1 0.2557905 0.03789321
#> 2 0.2734807 0.03632590
#> 3 0.2716440 0.03643221
#> 4 0.2711929 0.03621119
#> 5 0.2743809 0.03625556
#> 6 0.2716777 0.03629290
#> 7 0.2765445 0.03604140
#> 8 0.2747726 0.03619440
#> 9 0.2710448 0.03650991
#> 10 0.2699076 0.03685035York regression models are also implemented in
bayclumpr:
cal.york(data = ds, replicates = 10)
#> alpha beta
#> 1 0.2502012 0.03822585
#> 2 0.2417208 0.03889872
#> 3 0.2430982 0.03835195
#> 4 0.2689329 0.03625636
#> 5 0.2714125 0.03611063
#> 6 0.2638400 0.03685362
#> 7 0.2897728 0.03520313
#> 8 0.2755461 0.03600768
#> 9 0.2630668 0.03715828
#> 10 0.2678943 0.03654801Finally, bayclumpr implements three types of Bayesian
linear models that are used for calibrations and temperature
reconstructions. Let’s fit all three models using the
cal.bayesian function:
BayesCal <- cal.bayesian(calibrationData = ds)The results are here stored in the BayesCal object and
corresponds to stan objects summarizing posterior
distributions of the parameters:
BayesCal
#> $BLM1_fit
#> Inference for Stan model: anon_model.
#> 2 chains, each with iter=2500; warmup=1000; thin=1;
#> post-warmup draws per chain=1500, total post-warmup draws=3000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> alpha 0.27 0.00 0.01 0.24 0.26 0.27 0.27 0.29 2566 1
#> beta 0.04 0.00 0.00 0.03 0.04 0.04 0.04 0.04 2593 1
#> sigma 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.02 3138 1
#> log_lik[1] 3.29 0.00 0.13 3.01 3.21 3.30 3.37 3.51 2588 1
#> log_lik[2] 3.34 0.00 0.11 3.12 3.27 3.34 3.42 3.54 2970 1
#> log_lik[3] 2.71 0.00 0.18 2.31 2.59 2.72 2.83 3.02 2638 1
#> log_lik[4] 2.95 0.00 0.21 2.47 2.82 2.97 3.09 3.30 2742 1
#> log_lik[5] 2.50 0.00 0.21 2.03 2.37 2.52 2.65 2.87 2838 1
#> log_lik[6] 3.11 0.00 0.12 2.85 3.03 3.11 3.19 3.33 2636 1
#> log_lik[7] 3.31 0.00 0.11 3.09 3.24 3.31 3.39 3.50 2921 1
#> log_lik[8] 1.50 0.01 0.51 0.38 1.19 1.55 1.87 2.39 3203 1
#> log_lik[9] 3.01 0.00 0.19 2.59 2.89 3.03 3.15 3.33 2844 1
#> log_lik[10] 2.81 0.00 0.26 2.20 2.66 2.84 2.99 3.23 2904 1
#> log_lik[11] 3.25 0.00 0.12 3.00 3.18 3.25 3.33 3.46 2985 1
#> log_lik[12] 3.36 0.00 0.11 3.12 3.29 3.36 3.43 3.56 2855 1
#> log_lik[13] 3.36 0.00 0.11 3.14 3.29 3.36 3.43 3.56 3039 1
#> log_lik[14] 1.60 0.01 0.39 0.76 1.35 1.63 1.87 2.26 3173 1
#> log_lik[15] 3.33 0.00 0.11 3.11 3.26 3.33 3.40 3.53 3042 1
#> log_lik[16] 3.36 0.00 0.11 3.14 3.29 3.36 3.43 3.55 3091 1
#> log_lik[17] 3.30 0.00 0.12 3.03 3.22 3.31 3.38 3.52 2626 1
#> log_lik[18] 3.37 0.00 0.11 3.14 3.30 3.37 3.44 3.57 2986 1
#> log_lik[19] 3.25 0.00 0.14 2.91 3.16 3.26 3.35 3.48 3035 1
#> log_lik[20] 3.34 0.00 0.11 3.12 3.27 3.34 3.41 3.53 3005 1
#> log_lik[21] 2.58 0.00 0.22 2.09 2.45 2.60 2.74 2.96 2731 1
#> log_lik[22] 3.31 0.00 0.12 3.04 3.23 3.32 3.39 3.53 2631 1
#> log_lik[23] 3.21 0.00 0.11 2.97 3.14 3.22 3.29 3.42 2728 1
#> log_lik[24] 3.23 0.00 0.16 2.85 3.15 3.25 3.34 3.48 2822 1
#> log_lik[25] 3.18 0.00 0.12 2.91 3.11 3.19 3.26 3.41 2655 1
#> log_lik[26] 3.27 0.00 0.11 3.04 3.20 3.28 3.35 3.48 2986 1
#> log_lik[27] 3.04 0.00 0.19 2.58 2.93 3.06 3.18 3.36 2752 1
#> log_lik[28] 2.99 0.00 0.19 2.56 2.88 3.01 3.13 3.31 2794 1
#> log_lik[29] 1.83 0.01 0.34 1.08 1.62 1.85 2.07 2.41 3005 1
#> log_lik[30] 2.87 0.00 0.22 2.38 2.74 2.90 3.03 3.24 2779 1
#> log_lik[31] 3.11 0.00 0.16 2.75 3.01 3.12 3.22 3.38 2801 1
#> log_lik[32] 3.36 0.00 0.11 3.13 3.29 3.36 3.43 3.56 2948 1
#> log_lik[33] 2.20 0.01 0.32 1.48 2.01 2.23 2.44 2.74 2928 1
#> log_lik[34] 2.32 0.01 0.29 1.67 2.14 2.35 2.53 2.83 3076 1
#> log_lik[35] 3.37 0.00 0.11 3.14 3.30 3.37 3.44 3.57 3027 1
#> log_lik[36] 0.85 0.01 0.58 -0.40 0.49 0.90 1.26 1.83 3258 1
#> log_lik[37] 2.65 0.01 0.31 1.93 2.48 2.69 2.87 3.15 2913 1
#> log_lik[38] 3.22 0.00 0.11 2.99 3.15 3.23 3.30 3.42 2977 1
#> log_lik[39] 3.09 0.00 0.16 2.74 2.99 3.10 3.20 3.35 2906 1
#> log_lik[40] 3.08 0.00 0.16 2.73 2.99 3.10 3.19 3.36 2779 1
#> log_lik[41] 2.65 0.00 0.24 2.10 2.51 2.67 2.81 3.04 3034 1
#> log_lik[42] 1.71 0.01 0.37 0.89 1.47 1.74 1.97 2.35 2971 1
#> log_lik[43] 3.11 0.00 0.22 2.58 3.00 3.15 3.27 3.45 2603 1
#> log_lik[44] 3.37 0.00 0.11 3.13 3.30 3.37 3.44 3.57 2949 1
#> log_lik[45] 2.30 0.00 0.26 1.71 2.15 2.32 2.48 2.75 3088 1
#> log_lik[46] 3.29 0.00 0.16 2.88 3.22 3.32 3.40 3.55 1936 1
#> log_lik[47] 2.84 0.00 0.15 2.50 2.75 2.85 2.95 3.12 2728 1
#> log_lik[48] 3.20 0.00 0.16 2.84 3.10 3.21 3.31 3.46 2639 1
#> log_lik[49] 2.08 0.01 0.30 1.43 1.90 2.10 2.29 2.58 3158 1
#> log_lik[50] 3.30 0.00 0.12 3.02 3.22 3.30 3.38 3.52 2666 1
#> lp__ 136.49 0.14 5.13 125.45 133.17 136.81 140.20 145.80 1332 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jul 22 08:17:58 2025.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> $BLM1_fit_NoErrors
#> Inference for Stan model: anon_model.
#> 2 chains, each with iter=2500; warmup=1000; thin=1;
#> post-warmup draws per chain=1500, total post-warmup draws=3000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> alpha 0.27 0.00 0.01 0.24 0.26 0.27 0.28 0.29 1104 1
#> beta 0.04 0.00 0.00 0.03 0.04 0.04 0.04 0.04 1113 1
#> sigma 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.02 923 1
#> log_lik[1] 3.29 0.00 0.12 3.04 3.21 3.29 3.37 3.51 1229 1
#> log_lik[2] 3.34 0.00 0.10 3.12 3.27 3.34 3.41 3.54 1088 1
#> log_lik[3] 2.70 0.00 0.17 2.32 2.59 2.72 2.82 3.01 2779 1
#> log_lik[4] 2.95 0.01 0.20 2.53 2.82 2.97 3.09 3.28 1435 1
#> log_lik[5] 2.50 0.00 0.21 2.04 2.37 2.51 2.64 2.85 2387 1
#> log_lik[6] 3.11 0.00 0.11 2.87 3.03 3.11 3.18 3.32 2451 1
#> log_lik[7] 3.31 0.00 0.10 3.09 3.24 3.31 3.38 3.51 1220 1
#> log_lik[8] 1.51 0.01 0.51 0.40 1.19 1.56 1.88 2.37 1498 1
#> log_lik[9] 3.01 0.01 0.19 2.56 2.90 3.03 3.15 3.31 1312 1
#> log_lik[10] 2.80 0.01 0.26 2.22 2.63 2.83 2.99 3.22 1692 1
#> log_lik[11] 3.25 0.00 0.12 3.00 3.17 3.26 3.33 3.46 1065 1
#> log_lik[12] 3.36 0.00 0.11 3.14 3.29 3.36 3.43 3.56 893 1
#> log_lik[13] 3.36 0.00 0.11 3.14 3.29 3.37 3.43 3.56 901 1
#> log_lik[14] 1.59 0.01 0.38 0.72 1.35 1.62 1.86 2.24 1581 1
#> log_lik[15] 3.33 0.00 0.10 3.11 3.26 3.33 3.40 3.52 1029 1
#> log_lik[16] 3.36 0.00 0.10 3.14 3.29 3.36 3.43 3.55 924 1
#> log_lik[17] 3.30 0.00 0.12 3.06 3.22 3.30 3.38 3.52 1191 1
#> log_lik[18] 3.37 0.00 0.10 3.15 3.30 3.37 3.44 3.56 938 1
#> log_lik[19] 3.24 0.00 0.14 2.92 3.16 3.26 3.34 3.49 1197 1
#> log_lik[20] 3.34 0.00 0.10 3.12 3.27 3.34 3.41 3.53 1099 1
#> log_lik[21] 2.59 0.00 0.21 2.10 2.46 2.60 2.73 2.94 1926 1
#> log_lik[22] 3.31 0.00 0.11 3.07 3.23 3.31 3.39 3.52 1168 1
#> log_lik[23] 3.21 0.00 0.11 2.99 3.14 3.21 3.28 3.42 1725 1
#> log_lik[24] 3.23 0.01 0.16 2.84 3.14 3.25 3.34 3.49 969 1
#> log_lik[25] 3.18 0.00 0.11 2.94 3.10 3.18 3.26 3.40 1737 1
#> log_lik[26] 3.27 0.00 0.11 3.04 3.20 3.28 3.35 3.48 1018 1
#> log_lik[27] 3.03 0.00 0.19 2.60 2.91 3.06 3.17 3.34 1723 1
#> log_lik[28] 2.99 0.00 0.19 2.58 2.88 3.01 3.13 3.30 1755 1
#> log_lik[29] 1.83 0.01 0.33 1.08 1.65 1.87 2.06 2.38 2022 1
#> log_lik[30] 2.87 0.01 0.21 2.42 2.74 2.89 3.03 3.23 1448 1
#> log_lik[31] 3.10 0.00 0.16 2.75 3.00 3.12 3.21 3.37 1719 1
#> log_lik[32] 3.36 0.00 0.11 3.14 3.29 3.36 3.43 3.56 956 1
#> log_lik[33] 2.19 0.01 0.32 1.49 2.00 2.22 2.42 2.73 1815 1
#> log_lik[34] 2.33 0.01 0.29 1.72 2.15 2.35 2.54 2.81 1997 1
#> log_lik[35] 3.37 0.00 0.10 3.15 3.30 3.37 3.44 3.56 967 1
#> log_lik[36] 0.85 0.02 0.57 -0.42 0.49 0.91 1.27 1.83 1467 1
#> log_lik[37] 2.64 0.01 0.30 1.97 2.45 2.68 2.86 3.14 1653 1
#> log_lik[38] 3.22 0.00 0.11 3.01 3.16 3.23 3.30 3.42 1412 1
#> log_lik[39] 3.09 0.00 0.15 2.73 2.99 3.10 3.20 3.35 1308 1
#> log_lik[40] 3.09 0.00 0.15 2.76 2.99 3.10 3.19 3.34 1791 1
#> log_lik[41] 2.65 0.00 0.23 2.16 2.51 2.68 2.82 3.04 2122 1
#> log_lik[42] 1.71 0.01 0.36 0.89 1.51 1.74 1.96 2.31 1819 1
#> log_lik[43] 3.12 0.01 0.22 2.60 3.00 3.15 3.27 3.44 1365 1
#> log_lik[44] 3.37 0.00 0.11 3.15 3.29 3.37 3.44 3.56 937 1
#> log_lik[45] 2.29 0.01 0.25 1.72 2.14 2.31 2.47 2.72 2051 1
#> log_lik[46] 3.29 0.00 0.15 2.94 3.21 3.32 3.40 3.54 998 1
#> log_lik[47] 2.85 0.00 0.14 2.53 2.76 2.86 2.94 3.10 2862 1
#> log_lik[48] 3.20 0.00 0.15 2.86 3.11 3.22 3.30 3.45 1401 1
#> log_lik[49] 2.07 0.01 0.29 1.46 1.89 2.09 2.28 2.57 1560 1
#> log_lik[50] 3.30 0.00 0.11 3.06 3.22 3.30 3.38 3.52 1201 1
#> lp__ 185.91 0.04 1.20 182.72 185.35 186.20 186.80 187.31 888 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jul 22 08:18:22 2025.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> $BLM3_fit
#> Inference for Stan model: anon_model.
#> 2 chains, each with iter=2500; warmup=1000; thin=1;
#> post-warmup draws per chain=1500, total post-warmup draws=3000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> alpha[1] 0.27 0.00 0.01 0.24 0.26 0.27 0.27 0.29 979 1
#> beta[1] 0.04 0.00 0.00 0.03 0.04 0.04 0.04 0.04 981 1
#> sigma 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.02 926 1
#> log_lik[1] 3.29 0.00 0.12 3.02 3.22 3.30 3.37 3.51 1140 1
#> log_lik[2] 3.34 0.00 0.10 3.12 3.28 3.35 3.41 3.54 1021 1
#> log_lik[3] 2.71 0.00 0.19 2.30 2.59 2.72 2.84 3.02 2898 1
#> log_lik[4] 2.95 0.01 0.20 2.51 2.82 2.97 3.09 3.29 1434 1
#> log_lik[5] 2.50 0.00 0.22 2.03 2.37 2.52 2.66 2.87 2655 1
#> log_lik[6] 3.11 0.00 0.12 2.85 3.04 3.12 3.19 3.32 2390 1
#> log_lik[7] 3.31 0.00 0.11 3.09 3.24 3.32 3.38 3.51 1127 1
#> log_lik[8] 1.48 0.02 0.53 0.34 1.15 1.53 1.86 2.38 1212 1
#> log_lik[9] 3.01 0.01 0.19 2.58 2.89 3.03 3.15 3.32 1279 1
#> log_lik[10] 2.81 0.01 0.27 2.19 2.65 2.84 3.00 3.23 1385 1
#> log_lik[11] 3.25 0.00 0.12 3.00 3.17 3.25 3.33 3.46 1123 1
#> log_lik[12] 3.36 0.00 0.11 3.13 3.29 3.37 3.44 3.56 839 1
#> log_lik[13] 3.36 0.00 0.10 3.15 3.29 3.37 3.44 3.55 893 1
#> log_lik[14] 1.60 0.01 0.39 0.74 1.34 1.63 1.88 2.26 1673 1
#> log_lik[15] 3.33 0.00 0.11 3.12 3.26 3.33 3.40 3.52 1057 1
#> log_lik[16] 3.36 0.00 0.10 3.15 3.29 3.36 3.43 3.55 936 1
#> log_lik[17] 3.30 0.00 0.12 3.04 3.23 3.31 3.38 3.52 1104 1
#> log_lik[18] 3.37 0.00 0.10 3.15 3.30 3.37 3.44 3.56 907 1
#> log_lik[19] 3.25 0.00 0.14 2.91 3.17 3.26 3.35 3.48 1124 1
#> log_lik[20] 3.34 0.00 0.10 3.12 3.27 3.34 3.41 3.53 1027 1
#> log_lik[21] 2.58 0.00 0.22 2.10 2.44 2.60 2.73 2.95 1999 1
#> log_lik[22] 3.31 0.00 0.12 3.05 3.24 3.31 3.39 3.52 1084 1
#> log_lik[23] 3.21 0.00 0.11 2.98 3.15 3.22 3.29 3.42 1617 1
#> log_lik[24] 3.23 0.01 0.17 2.85 3.14 3.26 3.35 3.49 946 1
#> log_lik[25] 3.18 0.00 0.12 2.93 3.11 3.19 3.26 3.40 1666 1
#> log_lik[26] 3.27 0.00 0.11 3.03 3.20 3.28 3.35 3.48 1061 1
#> log_lik[27] 3.04 0.01 0.20 2.58 2.93 3.07 3.19 3.36 1350 1
#> log_lik[28] 2.98 0.01 0.20 2.55 2.86 3.00 3.12 3.32 1421 1
#> log_lik[29] 1.82 0.01 0.33 1.06 1.61 1.85 2.04 2.39 1766 1
#> log_lik[30] 2.87 0.01 0.21 2.41 2.74 2.89 3.03 3.23 1464 1
#> log_lik[31] 3.11 0.00 0.16 2.75 3.02 3.13 3.23 3.37 1430 1
#> log_lik[32] 3.36 0.00 0.11 3.14 3.29 3.36 3.43 3.56 914 1
#> log_lik[33] 2.21 0.01 0.33 1.48 2.00 2.24 2.44 2.75 1688 1
#> log_lik[34] 2.31 0.01 0.30 1.66 2.12 2.34 2.53 2.81 1524 1
#> log_lik[35] 3.37 0.00 0.10 3.15 3.30 3.37 3.44 3.56 934 1
#> log_lik[36] 0.82 0.02 0.59 -0.42 0.46 0.87 1.24 1.84 1238 1
#> log_lik[37] 2.65 0.01 0.31 1.93 2.47 2.69 2.88 3.15 1380 1
#> log_lik[38] 3.22 0.00 0.11 3.00 3.15 3.22 3.29 3.43 1490 1
#> log_lik[39] 3.09 0.00 0.16 2.75 2.99 3.10 3.20 3.35 1361 1
#> log_lik[40] 3.08 0.00 0.16 2.74 2.98 3.09 3.19 3.35 1539 1
#> log_lik[41] 2.64 0.01 0.24 2.11 2.49 2.66 2.81 3.05 1603 1
#> log_lik[42] 1.69 0.01 0.36 0.88 1.47 1.73 1.95 2.32 1645 1
#> log_lik[43] 3.11 0.01 0.22 2.59 2.98 3.14 3.26 3.45 1154 1
#> log_lik[44] 3.37 0.00 0.11 3.15 3.30 3.37 3.44 3.56 900 1
#> log_lik[45] 2.30 0.01 0.27 1.71 2.14 2.32 2.49 2.74 2187 1
#> log_lik[46] 3.30 0.01 0.15 2.91 3.22 3.31 3.40 3.54 729 1
#> log_lik[47] 2.84 0.00 0.15 2.51 2.74 2.86 2.94 3.11 3094 1
#> log_lik[48] 3.19 0.00 0.16 2.84 3.10 3.21 3.30 3.46 1252 1
#> log_lik[49] 2.08 0.01 0.30 1.44 1.89 2.10 2.29 2.58 1852 1
#> log_lik[50] 3.30 0.00 0.12 3.04 3.23 3.30 3.38 3.51 1117 1
#> lp__ 185.86 0.05 1.28 182.51 185.34 186.19 186.76 187.29 696 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jul 22 08:18:47 2025.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> attr(,"loo")
#> elpd_diff se_diff
#> BLM1_NE 0.0 0.0
#> BLM1_E 0.0 0.0
#> BLM3 -0.1 0.1Reconstructing temperatures in bayclumpr
bayclumpr implements two functions to perform
temperature reconstructions under frequentist (rec.clumped)
and Bayesian frameworks (rec.bayesian). Let’s review how
each of these functions work by generating a synthetic dataset for two
samples.
recData <- data.frame(Sample = paste("Sample", 1:9),
D47 = rep(c(0.6, 0.7, 0.8), 3),
D47error = c(rep(0.005,3), rep(0.01,3), rep(0.02,3)),
N = rep(2, 9),
Material = rep(1, 9))As for the calibration step, bayclumpr expects
uncertainty (D47error) to be expressed in terms of standard
deviation. Note that the recData object generated above
includes the smallest number of columns that are needed to perform
reconstructions in bayclumpr.
From this point, we will need to either specify the distribution of
parameter estimates from the calibration step. For instance, let’s
assume that we were interested in reconstructing temperatures for our
recData under an OLS model. First, we would have to run our
calibration analyses:
paramdist <- cal.ols(data = ds, replicates = 10)From this point, we can use the rec.clumped to
reconstruct temperatures based on the reconstruction dataset
(recData argument) and the observed calibration object
(obCal argument):
rec.clumped(recData = recData, obCal = paramdist)
#> Sample D47 D47error meanTemp error
#> 1 Sample 1 0.6 0.005 60.03159 2.524450
#> 2 Sample 2 0.7 0.005 18.33601 1.694860
#> 3 Sample 3 0.8 0.005 -10.81881 1.237513
#> 4 Sample 4 0.6 0.010 60.03159 4.992374
#> 5 Sample 5 0.7 0.010 18.33601 3.360496
#> 6 Sample 6 0.8 0.010 -10.81881 2.457676
#> 7 Sample 7 0.6 0.020 60.03159 9.766853
#> 8 Sample 8 0.7 0.020 18.33601 6.607381
#> 9 Sample 9 0.8 0.020 -10.81881 4.847547The resulting object includes information from the template
reconstruction dataset but also information on the reconstructed
temperature and associated uncertainty (1 SD). Let’s now
perform reconstructions but under a Bayesian framework. For this, we
will again need parameter estimates derived from the calibration step
(see the BayesCal created above). We will perform
reconstructions under only a single of the Bayesian models equivalent to
the OLS but fit under a Bayesian framework
(BayesCal$BLM1_fit_NoErrors).
PredsBay <- rec.bayesian(calModel = BayesCal$BLM1_fit_NoErrors, recData = recData, iter = 1000, postcalsamples = 100, MC = FALSE)
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 4.9e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.459 seconds (Warm-up)
#> Chain 1: 0.291 seconds (Sampling)
#> Chain 1: 0.75 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 2.7e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.384 seconds (Warm-up)
#> Chain 2: 0.29 seconds (Sampling)
#> Chain 2: 0.674 seconds (Total)
#> Chain 2:The associated reconstructions to this Bayesian model are shown below:
PredsBay
#> Sample D47 D47error meanTemp error
#> 1 Sample 1 0.6 0.005 59.57497 0.7403076
#> 2 Sample 2 0.7 0.005 18.66416 0.5285738
#> 3 Sample 3 0.8 0.005 -10.12629 0.3530472
#> 4 Sample 4 0.6 0.010 59.55347 0.8024151
#> 5 Sample 5 0.7 0.010 18.63890 0.5782041
#> 6 Sample 6 0.8 0.010 -10.12668 0.4293646
#> 7 Sample 7 0.6 0.020 59.51917 1.2027646
#> 8 Sample 8 0.7 0.020 18.65469 0.8807962
#> 9 Sample 9 0.8 0.020 -10.12319 0.5602955Including a co-variate (with error in the covariate)
Now, instead of having a single predictor, we will use a second column (i.e., Ion) and it’s error (i.e., IonError) as predictors. Let’s first create a calibration dataset.
ds <- cal.dataset(error = "S1", nobs = 50)Note that we will be using not super informative values for Ion and its error. We just need to include those two new columns in the dataset:
At this point, we can fit the regression model that accounts for uncertainty in each of the two predictors
ionmodel <- cal.ion.bayesian(calibrationData = ds, useIonError = TRUE)
ionmodel
#> $BLM1_fit
#> Inference for Stan model: anon_model.
#> 2 chains, each with iter=2500; warmup=1000; thin=1;
#> post-warmup draws per chain=1500, total post-warmup draws=3000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> alpha 0.27 0.00 0.01 0.24 0.26 0.27 0.28 0.30 6450 1
#> beta 0.04 0.00 0.00 0.03 0.04 0.04 0.04 0.04 6339 1
#> gamma 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 1793 1
#> sigma 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.02 4431 1
#> lp__ 99.57 0.26 7.16 84.74 94.80 99.82 104.67 112.48 776 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jul 22 08:19:38 2025.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).From here, we can simply run a reconstruction based on a synthetic reconstruction dataset with the following structure:
recData <- data.frame(Sample = paste("Sample", 1:9),
D47 = rep(c(0.6, 0.7, 0.8), 3),
D47error = c(rep(0.005,3), rep(0.01,3), rep(0.02,3)),
Ion = rep(c(0.6, 0.7, 0.8), 3),
IonError = c(rep(0.005,3), rep(0.01,3), rep(0.02,3)),
N = rep(2, 9),
Material = rep(1, 9))The reconstruction can be performed as follows:
PredsBay <- rec.ion.bayesian(calModel = ionmodel[[1]],
recData = recData,
iter = 1000,
postcalsamples = 100, MC = FALSE)
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 5.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.575 seconds (Warm-up)
#> Chain 1: 0.323 seconds (Sampling)
#> Chain 1: 0.898 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.534 seconds (Warm-up)
#> Chain 2: 0.634 seconds (Sampling)
#> Chain 2: 1.168 seconds (Total)
#> Chain 2:And the resulting reconstructions are as follow:
PredsBay
#> Sample D47 D47error meanTemp error
#> 1 Sample 1 0.6 0.005 59.83046 7.665063
#> 2 Sample 2 0.7 0.005 18.44528 5.025854
#> 3 Sample 3 0.8 0.005 -10.52767 3.740377
#> 4 Sample 4 0.6 0.010 59.90502 8.854410
#> 5 Sample 5 0.7 0.010 18.50280 5.825532
#> 6 Sample 6 0.8 0.010 -10.49236 4.318791
#> 7 Sample 7 0.6 0.020 60.25275 12.588328
#> 8 Sample 8 0.7 0.020 18.72431 8.353161
#> 9 Sample 9 0.8 0.020 -10.40502 6.101526Including a co-variate (without error in the covariate)
We will now do the same for the models outlined above but ignoring uncertainty in the covariate (Ion):
ds <- cal.dataset(error = "S1", nobs = 50)
ds$Ion <- rnorm(nrow(ds))
ionmodel <- cal.ion.bayesian(calibrationData = ds, useIonError = FALSE)
ionmodel
#> $BLM1_fit
#> Inference for Stan model: anon_model.
#> 2 chains, each with iter=2500; warmup=1000; thin=1;
#> post-warmup draws per chain=1500, total post-warmup draws=3000.
#>
#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> alpha 0.27 0.00 0.01 0.25 0.26 0.27 0.28 0.30 2908 1
#> beta 0.04 0.00 0.00 0.03 0.04 0.04 0.04 0.04 2931 1
#> gamma 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5064 1
#> sigma 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.02 2766 1
#> lp__ 136.05 0.14 5.24 124.71 132.76 136.44 139.86 145.28 1500 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Jul 22 08:20:27 2025.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).From here, we can just run the reconstructions
recData <- data.frame(Sample = paste("Sample", 1:9),
D47 = rep(c(0.6, 0.7, 0.8), 3),
D47error = c(rep(0.005,3), rep(0.01,3), rep(0.02,3)),
Ion = rep(c(0.6, 0.7, 0.8), 3),
N = rep(2, 9),
Material = rep(1, 9))
PredsBay <- rec.ion.bayesian(calModel = ionmodel[[1]],
recData = recData,
iter = 1000,
postcalsamples = 100, MC = FALSE,
useIonError = FALSE)
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 6.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.64 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.48 seconds (Warm-up)
#> Chain 1: 0.322 seconds (Sampling)
#> Chain 1: 0.802 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.46 seconds (Warm-up)
#> Chain 2: 0.594 seconds (Sampling)
#> Chain 2: 1.054 seconds (Total)
#> Chain 2:
PredsBay
#> Sample D47 D47error meanTemp error
#> 1 Sample 1 0.6 0.005 60.50896 7.835126
#> 2 Sample 2 0.7 0.005 18.71759 5.071897
#> 3 Sample 3 0.8 0.005 -10.46470 3.745609
#> 4 Sample 4 0.6 0.010 60.56685 9.001496
#> 5 Sample 5 0.7 0.010 18.76788 5.894212
#> 6 Sample 6 0.8 0.010 -10.43113 4.367098
#> 7 Sample 7 0.6 0.020 60.91822 12.720127
#> 8 Sample 8 0.7 0.020 18.95777 8.401952
#> 9 Sample 9 0.8 0.020 -10.33175 6.151552