1.1 The simulated dataset

We work with et_sim, a list bundled with the package. It represents a hypothetical mountain plant site where allele-frequency change on the arcsin-sqrt scale (z_diff) in two SNP clusters (A and B) is linked to three growing-season climate predictors (plus ten independent nuisance predictors d1..d10 with zero true effect).

Note on the arcsin-sqrt transformation. Raw allele frequencies \(f \in [0,\,1]\) are transformed via \(z = \arcsin(\sqrt{f})\). The result is unbounded (theoretically \([0,\,\pi/2] \approx [0,\,1.57]\) for a single chromosome, but changes \(\Delta z\) are unbounded in principle). Do not use response_scale = c(-1, 1) for this response — instead derive the range from the training data (the default when response_scale = NULL) or use a biologically motivated interval.

The dataset has a training period (1920–2019, 100 years per cluster) and a forecast period (2020–2034, 15 years per cluster). The forecast climate is warmer than the training period on average (reflecting a mild warming trend). Because the data are simulated, we know the true regression coefficients — stored in et_sim$true_params — which lets us validate parameter recovery and calibration.

True data-generating model:

\[ z\_diff_t = \alpha + \beta_1 \cdot Tmean_t + \beta_2 \cdot PPT_t + \beta_3 \cdot SWE_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0,\, \sigma^2) \]

Both clusters have a positive response to warming (Tmean) and a negative response to precipitation (PPT). They differ in the direction of their response to snowpack (SWE) and in residual variability — cluster B is noisier (σ = 0.35 vs. 0.30), which will be reflected in its wider forecast intervals and higher shelf-life ratio.

3.1 Visualise the training data

ggplot(et_sim$train, aes(x = year, y = z_diff, colour = cluster_id)) +
  geom_line() +
  geom_point(size = 2) +
  labs(
    x       = "Year",
    y       = expression(italic(z)[diff]),
    colour  = "Cluster"
  ) +
  et_theme()

# Climate predictors over the full period — warming trend visible in Tmean
all_years <- rbind(
  transform(et_sim$train[et_sim$train$cluster_id == "A",
                          c("year", "Tmean", "PPT", "SWE")],
            period = "Training"),
  transform(et_sim$forecast[et_sim$forecast$cluster_id == "A",
                             c("year", "Tmean", "PPT", "SWE")],
            period = "Forecast")
)

all_long <- reshape(all_years, direction = "long",
                    varying = c("Tmean", "PPT", "SWE"),
                    v.names = "value", timevar = "predictor",
                    times   = c("Tmean", "PPT", "SWE"))

ggplot(all_long, aes(x = year, y = value,
                     colour = period, group = period)) +
  geom_line() +
  geom_point(size = 1.8) +
  facet_wrap(~ predictor, ncol = 3) +
  labs(
    x      = "Year", y = "Standardised value", colour = "Period"
  ) +
  et_theme()

4.1 Step 1 — Elastic net for prior extraction

We fit a cross-validated elastic net (alpha = 0.5) on the training data. The elastic net performs regularised estimation (and variable selection when nuisance predictors d1..d10 are included) and provides coefficient magnitudes that scale the priors in the subsequent Bayesian model.

train_A <- et_sim$train[et_sim$train$cluster_id == "A", ]

x_train <- as.matrix(train_A[, c("Tmean", "PPT", "SWE")])
y_train <- train_A$z_diff

set.seed(42)
cv_fit_A <- glmnet::cv.glmnet(
  x     = x_train,
  y     = y_train,
  alpha = 0.5          # 0 = ridge, 1 = lasso, 0.5 = elastic net
)

# Coefficient estimates at lambda.min
coef(cv_fit_A, s = "lambda.min")
# 4 x 1 sparse Matrix of class "dgCMatrix"
#             lambda.min
# (Intercept)  0.0453
# Tmean        0.4983
# PPT         -0.2722
# SWE          0.1683

The elastic net closely recovers the signs and magnitudes of all three predictors (true values: Tmean = 0.50, PPT = −0.30, SWE = 0.20).

plot(cv_fit_A)

4.2 Step 2 — Extract informed priors

extract_priors() converts the elastic net coefficients into a brmsprior object. With the default shrinkage = "zero", each coefficient’s prior is centred at zero with SD set to multiplier × |coef| (floored at min_sd) — a regularising prior whose scale is data-informed without reusing the point estimate as a prior mean (pass shrinkage = "estimate" to centre the priors on the fitted coefficients, or use et_priors_split() to fit the first stage on a disjoint subset of rows).

prior_spec_A <- extract_priors(
  model      = cv_fit_A,
  lambda     = "lambda.min",
  multiplier = 2.0,   # prior SD = 2 × |elastic-net coef|
  min_sd     = 0.10   # floor: avoids spike priors on near-zero coefficients
)

print(prior_spec_A)
# ErrorTracer prior specification
#   Method      : glmnet
#   Predictors  : 3
#   Multiplier  : 2
#   Min SD      : 0.1
#   Shrinkage   : zero (prior means at 0; regularizing)
#   Coefficients:
#     Tmean                 prior mean =   0.0000  sd =   0.9966  (coef =   0.4983)
#     PPT                   prior mean =   0.0000  sd =   0.5444  (coef =  -0.2722)
#     SWE                   prior mean =   0.0000  sd =   0.3365  (coef =   0.1683)

The prior SDs (1.00 for Tmean, 0.54 for PPT, 0.34 for SWE) reflect the multiplier = 2 setting — wide enough for the likelihood to dominate on strong effects while regularising smaller ones.

Design note: The intercept is excluded from the informed coefficient priors, and intercept_prior_sd = NULL (default) lets brms use its data-aware default (student_t(3, median(y), mad(y))).

4.3 Step 3 — Fit the Bayesian model

et_fit() wraps brms::brm(), attaches the prior specification, and returns an et_model object.

fit_A <- et_fit(
  formula = z_diff ~ Tmean + PPT + SWE,
  data    = train_A,
  priors  = prior_spec_A,
  chains  = CHAINS,
  iter    = ITER,
  warmup  = ITER %/% 2L,
  cores   = CHAINS,
  seed    = SEED
)

print(fit_A)
# ErrorTracer model (et_model)
#   Formula : z_diff ~ Tmean + PPT + SWE
#   n obs   : 100
#   Chains  : 2   Iter: 1000   Warmup: 500
#   Priors  : informed (glmnet, 3 predictors)
#   Rhat max: 1.014

summary(fit_A)
# === ErrorTracer model summary ===
#
# --- Fixed effects ---
#              Estimate  Est.Error     Q2.5    Q97.5
# Intercept    0.041      0.104    -0.162    0.246
# Tmean        0.666      0.113     0.438    0.888
# PPT         -0.487      0.200    -0.909   -0.123
# SWE          0.406      0.206     0.023    0.827

The posterior means are shrunk slightly relative to the elastic-net estimates, particularly for PPT and SWE — the likelihood is pulling the estimates toward the true values (Tmean = 0.50, PPT = −0.30, SWE = 0.20). All 95% credible intervals exclude zero for Tmean and PPT; SWE is positive but its lower bound just crosses zero, consistent with the true value being 0.20 with σ = 0.40 and only 20 training observations.

4.4 Step 4 — Diagnose convergence and model fit

diag_A <- et_diagnose(fit_A, loo = TRUE)

# Convergence
diag_A$convergence
# $rhat_max
# [1] 1.004         # all chains well below 1.05 threshold
#
# $rhat_all_ok
# [1] TRUE
#
# $neff_min
# [1] 0.275         # lowest effective-sample-size ratio
#
# $neff_all_ok
# [1] TRUE          # all Neff ratios > 0.10
#
# $n_divergences
# [1] 0             # no divergent transitions

diag_A$loo[c("elpd_loo", "p_loo", "looic", "n_bad_pareto_k")]
# $elpd_loo
# [1] -14.55
#
# $p_loo
# [1] 4.00          # effective parameters ≈ 4 (intercept, 3 betas); sigma
#                   # is implicit — consistent with a 4-parameter model
# $looic
# [1] 29.10
#
# $n_bad_pareto_k
# [1] 0             # no influential observations (Pareto k > 0.7)

All Rhat < 1.05, zero divergences, and no Pareto-k warnings — the sampler has converged and the posterior is not dominated by any single observation.

4.5 Step 5 — Predict with environmental noise propagation

et_predict() generates posterior predictive draws and propagates uncertainty about the predictor values themselves via the env_noise argument. Each perturbation draw adds Gaussian noise to the predictor values before computing the linear predictor, capturing sensitivity to measurement or forecast error in the climate inputs.

env_noise accepts three forms: - A named list of scalars — constant noise SD per predictor across all observations (suitable for short validation windows). - A named list of vectors of length nrow(newdata) — time-varying noise, where each observation gets its own SD. This is the recommended form for GCM projections, where ensemble spread grows with forecast horizon. - A single scalar — applied as a fraction of each predictor’s empirical SD.

fcast_A <- et_sim$forecast[et_sim$forecast$cluster_id == "A", ]

pred_A <- et_predict(
  model     = fit_A,
  newdata   = fcast_A,
  env_noise = et_sim$env_noise,   # Tmean = 0.30, PPT = 0.20, SWE = 0.15, d1..d10 = 0
  n_draws   = 1000L,
  ci_levels = c(0.50, 0.80, 0.90, 0.95)
)

print(pred_A)
# ErrorTracer prediction (et_prediction)
#   Observations  : 15
#   Draws         : 1000
#   CI levels     : 0.5, 0.8, 0.9, 0.95
#   Interval type : predictive
#   Mean var decomposition (across observations):
#     Parameter    : 0.0040
#     Environmental: 0.0264  (MC SE: 0.0014)
#     Residual     : 0.0870
#     Total        : 0.1174  (measured predictive variance)
#     Closure gap  : 0.0001  (Monte Carlo; compare MC SE above)

In this dataset the environmental uncertainty is several times larger than the parameter uncertainty — with 100 training observations the regression coefficients are well estimated, so measurement noise in the climate inputs contributes substantially more forecast variance than coefficient imprecision. Residual process noise remains the largest single component.

# 90% credible intervals by forecast year
ci90 <- pred_A$credible_intervals[pred_A$credible_intervals$ci_level == 0.90, ]
ci90$year <- fcast_A$year
head(ci90[, c("year", "lower", "median", "upper", "width")])
#    year      lower     median     upper    width
# 31 2020 -0.9124320 -0.3710932 0.2328418 1.145274
# 32 2021 -0.4276259  0.1268306 0.6656614 1.093287
# 33 2022 -0.7306874 -0.1616093 0.4136616 1.144349
# 34 2023 -0.2588986  0.2843572 0.8218421 1.080741
# 35 2024 -0.4706744  0.1134329 0.6489682 1.119643
# 36 2025 -0.3844870  0.1803172 0.7336616 1.118149

4.6 Step 6 — Uncertainty decomposition

decomp_A <- decompose_uncertainty(pred_A)
head(decomp_A)
#   obs_id total_var param_var env_var v_env_mcse residual_var closure_gap
# 1      1    0.1241    0.0061  0.0279     0.0054       0.0894      0.0007
# 2      2    0.1218    0.0049  0.0262     0.0053       0.0894      0.0013
# ...

residual_var is constant across observations (for a Gaussian identity model) — it equals the posterior mean of σ² and represents process noise that cannot be reduced by better predictors or more training data. total_var is the measured variance of the environmentally inclusive posterior predictive draws, and closure_gap records the Monte Carlo remainder total_var - (param_var + env_var + residual_var).

et_plot_decomposition(decomp_A, proportional = TRUE)

The proportional view makes it clear that process noise (residual) accounts for ~80% of total forecast variance, parameter uncertainty for ~12–20%, and environmental measurement noise for ~17–25%. To substantially improve forecast precision for this system, the primary target should be the residual variance — identifying additional drivers of allele-frequency change.

Autocorrelation note. If you fit the model with a residual autocorrelation term (e.g. et_fit(y ~ x + ar(time = t, p = 1), ...)), decompose_uncertainty() adds a fourth column temporal_var, the remainder total_var − (param_var + env_var + residual_var). Since 1.4.0 it is signed: it is a covariance term, and a negative value is meaningful, since autocorrelated errors can reduce predictive variance relative to an independent-innovation baseline. It captures the autocorrelation-induced spread that brms::posterior_predict() accumulates iteratively beyond a single innovation, and residual_var is then interpretable as the innovation variance rather than the stationary marginal variance. et_plot_decomposition() shows the fourth segment automatically when the column is present.

4.7 Step 7 — Forecast shelf life

The shelf life analysis asks: at which forecast year does the 90% CI width exceed the biologically plausible response scale? Because z_diff is on the arcsin-sqrt scale (unbounded), we derive the response scale from the training data (which is also the default when response_scale = NULL).

# Derive response scale from training data (arcsin-sqrt scale — unbounded)
pr_A <- range(train_A$z_diff)   # c(-1.73, 1.57) => width ~ 3.30

sl_A <- shelf_life(
  predictions              = pred_A,
  response_scale           = pr_A,
  ci_level                 = 0.90,
  threshold                = 1.0,        # uninformative when CI width >= response scale
  time_col                 = "year",
  min_run                  = 2L,         # require >= 2 consecutive exceedances
  projection_alpha         = 0.05,       # require significantly positive slope to project
  min_slope_for_projection = 1e-4,       # minimum slope magnitude to project
  max_extrapolation_factor = 10          # cap: project <= last_year + 10 * window_width
)

print(sl_A)

head(as.data.frame(sl_A))

All forecast years are strongly informative (ratio well below the threshold of 1.0). Because the CI/range ratio neither sustains a crossing (min_run = 2) nor exhibits a statistically significant upward slope (projection_alpha = 0.05), the result is reported as a lower bound (> 2034): the forecast remains informative across the entire observed forecast window.

Interpreting the three horizon modes:

Mode	When it applies	What is returned
observed	Threshold crossed and sustained for min_run periods within the window	First sustained crossing time
projected	No sustained crossing; significantly positive trend within cap	Extrapolated crossing time + se_t_star
lower_bound	No sustained crossing; no significant trend, or projection beyond cap	> last observed time

Increase max_extrapolation_factor (or set it to Inf) to allow longer projections if the biological context supports it.

et_plot_shelf_life(sl_A, show_ratio = TRUE)

4.8 Step 8 — Calibration assessment

We compare the posterior predictive intervals against the held-out true values in et_sim$validation.

valid_A <- et_sim$validation[et_sim$validation$cluster_id == "A", ]

cal_A <- et_calibrate(
  predictions  = pred_A,
  observed     = valid_A,
  response_col = "z_diff",
  ci_levels    = c(0.50, 0.80, 0.90, 0.95)
)

print(cal_A)
#   ci_level nominal observed_coverage n_obs calibration_error
# 1     0.50    0.50              0.60     5              0.10
# 2     0.80    0.80              0.80     5              0.00
# 3     0.90    0.90              1.00     5              0.10
# 4     0.95    0.95              1.00     5              0.05

The 80% CI achieves exactly nominal coverage. The 90% and 95% CIs are slightly over-conservative (all 5 held-out values fall inside), which is expected with only 5 validation observations and regularising priors. Coverage estimates based on n = 5 have high sampling variance; the important result is that no CI is systematically under-covering.

et_plot_calibration(cal_A)

5.1 Forecast fan chart

et_plot_forecast(
  predictions  = pred_A,
  observed     = valid_A,
  response_col = "z_diff",
  time_col     = "year"
)

Nested credible-interval ribbons (darker = narrower interval) overlay the posterior median. Red points are the held-out validation values. All five observed values fall within the 95% CI.

5.2 Prior versus posterior

et_plot_prior_posterior(fit_A, max_preds = 3L)

Each panel overlays the prior density (grey, wide) with the posterior density (blue, narrower). A tightly shifted posterior means the data are informative; a posterior that closely mirrors the prior means that the likelihood has little leverage on that coefficient. Tmean shows the largest update, consistent with its largest true effect size.

5.3 Forest plot: Bayesian vs. elastic net

et_plot_coefficients(fit_A)

Red crosses mark the elastic-net estimates used as prior means. The Bayesian posterior means (blue points) with their 95% CIs show how the likelihood shifted the estimates. For all three predictors the posterior is pulled slightly toward zero relative to the elastic-net estimate, reflecting regularisation by the wide but informative priors.

7.1 Fitting

# Build one et_prior_spec per cluster from individual elastic nets
prior_list <- lapply(c("A", "B"), function(cl) {
  df <- et_sim$train[et_sim$train$cluster_id == cl, ]
  x  <- as.matrix(df[, c("Tmean", "PPT", "SWE")])
  set.seed(42)
  suppressWarnings(cv <- glmnet::cv.glmnet(x, df$z_diff, alpha = 0.5))
  extract_priors(cv, multiplier = 2.0, min_sd = 0.10)
})
names(prior_list) <- c("A", "B")

# One brms model per cluster, returned as et_model_list
fit_all <- et_fit(
  formula  = z_diff ~ Tmean + PPT + SWE,
  data     = et_sim$train,
  priors   = prior_list,       # named list: one et_prior_spec per cluster
  grouping = "cluster_id",
  chains   = CHAINS,
  iter     = ITER,
  warmup   = ITER %/% 2L,
  cores    = CHAINS,
  seed     = SEED
)

print(fit_all)
# ErrorTracer grouped model list (et_model_list)
#   Grouping : cluster_id
#   Formula  : z_diff ~ Tmean + PPT + SWE
#   Groups   : 2
#   Fitted   : 2 / 2

summary(fit_all)
# --- Per-group Rhat max ---
#   A                     Rhat max = 1.004
#   B                     Rhat max = 1.007

7.2 Diagnostics

diag_all <- et_diagnose(fit_all, loo = TRUE)
diag_all$summary
#   group rhat_ok neff_ok n_divergences elpd_loo n_bad_pareto_k
# 1     A    TRUE    TRUE             0   -14.55              0
# 2     B    TRUE    TRUE             0   -17.15              1

Cluster B has one observation with Pareto-k > 0.7, which is a mild flag worth investigating (it suggests that observation has an outsized influence on the posterior). With only 20 training points this is not unusual.

7.3 Prediction and decomposition

pred_all <- et_predict(
  model     = fit_all,
  newdata   = et_sim$forecast,
  env_noise = et_sim$env_noise,
  n_draws   = 1000L,
  ci_levels = c(0.50, 0.80, 0.90, 0.95)
)

print(pred_all)
# ErrorTracer grouped predictions (et_prediction_list)
#   Grouping : cluster_id
#   Groups   : 2 / 2 predicted

decomp_all <- decompose_uncertainty(pred_all)
head(decomp_all)

Cluster B has a systematically larger residual_var, directly reflecting its higher σ (0.35 vs. 0.30).

7.4 Shelf life comparison

# Use the pooled training range across both clusters (arcsin-sqrt scale)
pr_all <- range(et_sim$train$z_diff)   # c(-1.73, 1.57) => width ~ 3.30

sl_all <- shelf_life(
  pred_all,
  response_scale           = pr_all,
  ci_level                 = 0.90,
  threshold                = 1.0,
  time_col                 = "year",
  max_extrapolation_factor = 10
)

print(sl_all)

Both clusters remain strongly informative across all 15 forecast years. Cluster B’s mean ratio is noticeably higher than cluster A’s, directly reflecting cluster B’s higher residual variance (σ = 0.35 vs. 0.30). Both are well below the uninformative threshold of 1.0, and neither crosses within the 2020–2034 window, so both are reported as lower bounds.

The difference in ratios has a practical implication: under a longer forecast horizon or with larger environmental noise inputs, cluster B would approach the uninformative threshold first — a direct consequence of its higher biological process noise.

et_plot_shelf_life(sl_all, show_ratio = TRUE) +
  labs(
    colour   = "Cluster"
  )

7.5 Calibration and forecast plots

cal_all <- et_calibrate(
  pred_all,
  observed     = et_sim$validation,
  response_col = "z_diff",
  ci_levels    = c(0.50, 0.80, 0.90, 0.95)
)

et_plot_calibration(cal_all)

# Forest plot for all groups in a single faceted figure
et_plot_coefficients(fit_all)