The methods paper for ErrorTracer contains a one-page
Box 1 that walks through every stage of the pipeline
using a single predictor and a single response. The box gives the
intuition in a compressed form. This vignette is the long-form
companion: same example, same numbers, but with the code, the math, and
the diagnostic plots laid out one stage at a time so a reader new to
Bayesian modelling can follow along, run it, and modify it.
The whole example fits one model: phenology shift (days) as a function of a single standardised temperature anomaly. We know the data-generating coefficients because we simulate the data ourselves — so you can see how well each stage recovers what was put in.
Shelf life is a temporal concept (“for how many years is my
forecast still useful?”), so the toy data is generated as a year-indexed
time series. The temperature predictor x drifts upward at a
fixed warming rate, the response y follows it through a
linear coefficient \(\beta\), and the
training window is the past while the forecast window is the future.
set.seed(20260523L)
n_train <- 12L # 12 training years
true_beta <- 4.0 # days of phenology shift per SD of temperature
true_sigma <- 1.5 # residual SD (days)
true_alpha <- 0.0 # intercept
warming_slope <- 0.15 # SD of standardised temperature per year
year0 <- 2000L # year at which the trend is referenced to zero
years_train <- year0 + seq_len(n_train) - 1L # 2000..2011
train <- data.frame(
year = years_train,
x = warming_slope * (years_train - year0) +
rnorm(n_train, 0, 0.5) # noisy observed temperature
)
train$y <- true_alpha + true_beta * train$x +
rnorm(n_train, 0, true_sigma)We have 12 noisy training observations along the warming trajectory. The data-generating slope is \(\beta = 4\) days per SD of temperature; everything below should recover something close to this.
extract_priors() is the bridge from a “fast,
regularized” world to a “slow, Bayesian” world. The recipe is the same
for lm, glm, cv.glmnet, and
ranger: take the point estimate (or coefficient) as the
prior mean and scale the prior SD by the standard error (or by the
regularised coefficient magnitude, for glmnet).
lm_fit <- lm(y ~ x, data = train)
summary(lm_fit)$coefficients
priors <- extract_priors(lm_fit, multiplier = 2.0, min_sd = 0.1)
priorsConcretely, with \(m = 2\) (the
multiplier default) and \(\sigma_{\min} = 0.1\), the prior on \(\beta\) is
\[ \beta \sim \mathcal{N}\!\left(\hat\beta_{\text{lm}}, \bigl[\max(2 \cdot \mathrm{SE},\, \sigma_{\min})\bigr]^2\right). \]
For our simulated data, \(\hat\beta_{\text{lm}} \approx 4.02\) with \(\mathrm{SE} \approx 0.39\), so the prior is roughly \(\mathcal{N}(4.02, 0.79^2)\).
bfit <- et_fit(
formula = y ~ x,
data = train,
priors = priors,
chains = 4L, iter = 2000L, warmup = 1000L,
cores = 4L, seed = 1L, silent = 2L, refresh = 0
)
print(bfit)Under the hood et_fit() calls brms::brm()
with the prior we supplied for \(\beta\) plus weakly-informative defaults
for the intercept and the residual scale. The output is a full
posterior:
beta_draws <- as.matrix(bfit$fit, variable = "b_x")[, 1]
c(mean = mean(beta_draws), sd = sd(beta_draws))For our seed, the posterior collapses to \(\beta \approx 4.01\) with SD \(\approx 0.38\). The data didn’t drag the prior very far because the prior was already centred on the lm estimate — but the posterior SD is sharper than the prior SD, which is the regular Bayesian shrinkage we wanted.
A picture is worth a thousand draws. The package ships a helper:
The forecast window is years 2012–2045. At each future year, the
projected temperature is the trend-line value \(x(t) = \mathrm{warming\_slope}\cdot(t -
2000)\). We pass env_noise = list(x = 0.3) to tell
et_predict() that the projected temperature is itself
uncertain by \(\pm 0.3\) SD (a
reasonable error for a downscaled climate projection).
year_min <- year0
year_max <- year0 + 45L
years_grid <- seq(year_min, year_max, by = 1L)
forecast_df <- data.frame(
year = years_grid,
x = warming_slope * (years_grid - year0)
)
pred <- et_predict(
model = bfit,
newdata = forecast_df,
env_noise = list(x = 0.3),
n_draws = 2000L,
n_perturb = 500L,
ci_levels = 0.90,
include_env_in_ci = TRUE
)Setting include_env_in_ci = TRUE folds the environmental
noise into the reported credible intervals. With it FALSE
(the default), the CI captures parameter + residual uncertainty only;
environmental variance is still reported separately in
decompose_uncertainty(), but the user-facing CI is “what
would the response be if we knew the future temperature exactly”. The
choice depends on whether you’re forecasting under projection
uncertainty or given a known driver.
decompose_uncertainty() returns one row per forecast
year with the three (or four, if you have autocorrelation) variance
components:
Three things to notice as the forecast horizon extends into the future (year increases, so \(x\) on the warming trajectory increases too):
param_var grows roughly as \(\hat\sigma_\beta^2 \cdot x(t)^2\). That’s
the extrapolation cost of being uncertain about \(\beta\): the further past the training
window you forecast, the more it matters.env_var is roughly constant at \(\hat\beta^2 \cdot \sigma_x^2\). It scales
with how steep the dose–response is and how poorly we know the
future dose, but it does not grow with horizon.residual_var is the biological-noise floor. It does not
change with the forecast year because it’s a property of the response,
not of when you forecast.A quick numerical look at a representative mid-forecast year:
A forecast is informative when the predictive CI is narrower than the plausible range of the response. The plausible range here is “any phenology shift between \(-4\) and \(+4\) days” — anything wider than that and we are basically saying “we have no idea”.
sl <- shelf_life(
predictions = pred,
response_scale = c(-4, 4),
ci_level = 0.90,
threshold = 1.0,
time_col = "year"
)
attr(sl, "horizon")shelf_life() walks along the forecast time axis and
reports the first year at which the 90 % CI width first exceeds
the response range. For this example it lands at \(t^{\star} \approx 2027\) — i.e. the model
still produces a useful forecast for roughly 16 years past the end of
training; beyond that, the credible interval is wider than the plausible
response and should not be over-interpreted. That answer — a calendar
year, not a predictor value — is the quantity ecologists and
conservation practitioners actually need.
The four diagnostic panels — prior vs posterior, the nested CI fan,
the stacked variance components, and the shelf-life curve — are produced
together by the script that generates Box 1 of the paper. The exact same
code is what appears in the box’s figure; see the package source under
src/box1_worked_example.R in the proof-of-concept
repository (or run the chunk below in your own session):
Two limits of this minimal walk-through deserve flagging:
extract_priors(). The fundamental
decomposition is the same.et_fit(..., grouping = "g") runs the same pipeline
independently per level of a grouping factor (per SNP cluster, per
population, per site). The toy example fits a single global model; the
multi-group form returns a named list of et_model objects,
and decompose_uncertainty() / shelf_life()
preserve the grouping.sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=es_CU.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=es_CU.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=es_CU.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=es_CU.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/Chicago
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0 xfun_0.57
#> [5] cachem_1.1.0 knitr_1.51 htmltools_0.5.9 rmarkdown_2.31
#> [9] lifecycle_1.0.5 cli_3.6.6 sass_0.4.10 jquerylib_0.1.4
#> [13] compiler_4.6.0 tools_4.6.0 evaluate_1.0.5 bslib_0.10.0
#> [17] yaml_2.3.12 otel_0.2.0 rlang_1.2.0 jsonlite_2.0.0