Bayesian Error Propagation and Forecast Uncertainty Decomposition

Provides a full pipeline from regularized or standard regression models (elastic net, linear models, generalized linear models, random forests) to informed Bayesian priors, structured forecast uncertainty decomposition (parameter / environmental / residual, plus a temporal component when the model carries an autocorrelation term), and forecast shelf life analysis (the quantification of when a forecast becomes uninformative). Designed for ecological and genomic forecasting with climate or environmental covariates. Methods build on Bürkner (2017) for Bayesian regression via 'Stan', Friedman, Hastie, and Tibshirani (2010) for elastic net regularization, Wright and Ziegler (2017) for random forests, and Vehtari, Gelman, and Gabry (2017) for leave-one-out cross-validation.

ErrorTracer: Bayesian Error Propagation and Forecast Uncertainty Decomposition

Overview

ErrorTracer provides a complete pipeline for ecological and genomic forecasting with climate or environmental covariates. It bridges the gap between regularized regression and fully Bayesian workflows in three steps:

1. Extract informed priors from a fitted regularized or standard model (glmnet, lm, glm, ranger) and carry those estimates forward as prior means into a Bayesian model — rather than discarding them.
2. Decompose forecast uncertainty into three independent sources — parameter uncertainty, environmental/covariate measurement uncertainty, and residual process variance — plus a fourth temporal-autocorrelation component whenever the model formula carries an ar(), ma(), arma(), cosy(), unstr(), sar(), or car() term.
3. Quantify forecast shelf life — the time horizon at which a forecast's credible interval becomes wider than the biologically plausible response range, rendering it uninformative. This concept has no equivalent in any existing R package.

The package is designed for ecological and genomic time-series forecasting in which a regularized (or plain) regression model is refit as a Bayesian model and predictions must be accompanied by a principled uncertainty budget. The API is fully general — response, predictors, and grouping are formula-driven.

Key Features

Installation

Install the development version from GitHub using remotes:

# install.packages("remotes")
remotes::install_github("madrigalrocalj/ErrorTracer")

ErrorTracer requires a Stan backend. Install either cmdstanr (recommended) or rstan separately:

# Option A: cmdstanr (recommended)
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
cmdstanr::install_cmdstan()

# Option B: rstan
install.packages("rstan")

Core Workflow

Step 1 — Extract informed priors from a regularized model

extract_priors() is an S3 generic that dispatches on the class of the fitted model. Coefficient estimates (or importance weights for ranger) are used as prior means; prior SDs are proportional to coefficient magnitude or standard error.

library(ErrorTracer)
library(glmnet)

# --- Elastic net (cv.glmnet) ---
cv_fit <- cv.glmnet(x_train, y_train, alpha = 0.5)

prior_spec <- extract_priors(
  model              = cv_fit,
  lambda             = "lambda.min",  # or "lambda.1se"
  multiplier         = 2.0,           # prior SD = multiplier * |coef|
  min_sd             = 0.1,           # floor for prior SD
  intercept_prior_sd = 2.5,           # intercept prior SD
  sigma_prior_scale  = 1.0            # half-Cauchy scale for sigma
)

print(prior_spec)
# ErrorTracer prior specification
#   Method      : glmnet
#   Predictors  : 5
#   Multiplier  : 2
#   Min SD      : 0.1
#   Coefficients:
#     Tmean                 mean =   0.3241  sd =   0.6482
#     PPT                   mean =  -0.1087  sd =   0.2174
#     ...

Supported model classes and their prior construction:

For ranger, only variables with positive permutation importance are included. Because random forests do not produce signed coefficients, priors are centered at zero and scaled by importance — encoding "this predictor matters" without a directional claim.

Note on centering: The glmnet intercept is intentionally excluded from the prior — it reflects the centering of training data, which may not apply to new data. Center predictors upstream and document this in your analysis.

Step 2 — Fit the Bayesian model

et_fit() wraps brms::brm() and attaches the prior specification and training data for downstream use.

fit <- et_fit(
  formula       = z_diff ~ Tmean + PPT + SWE,
  data          = train_df,
  priors        = prior_spec,       # from extract_priors(); or NULL for brms defaults
  chains        = 4L,
  iter          = 2000L,
  warmup        = 1000L,
  cores         = 4L,
  seed          = 42L,
  adapt_delta   = 0.95,
  max_treedepth = 12L
)

print(fit)
# ErrorTracer model (et_model)
#   Formula : z_diff ~ Tmean + PPT + SWE
#   n obs   : 48
#   Chains  : 4  Iter: 2000  Warmup: 1000
#   Priors  : informed (glmnet, 3 predictors)
#   Rhat max: 1.002

summary(fit)

Group-level fitting

Pass a grouping column name to fit one model per group and receive an et_model_list:

# Fits one brms model per SNP cluster
fit_list <- et_fit(
  formula  = z_diff ~ Tmean + PPT + SWE,
  data     = all_clusters_df,
  priors   = prior_spec,        # single spec applied to all groups
  grouping = "cluster_id",      # or pass a named list of et_prior_spec for per-group priors
  chains   = 4L,
  iter     = 2000L
)

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

Step 3 — Predict with uncertainty propagation

et_predict() draws from the posterior predictive distribution, propagates environmental measurement noise through the model, and computes credible intervals.

pred <- et_predict(
  model     = fit,
  newdata   = forecast_df,
  env_noise = list(Tmean = 0.5, PPT = 10.0, SWE = 5.0),
  # env_noise can also be a single fraction of each predictor's SD:
  # env_noise = 0.10   →  10% noise applied uniformly
  n_draws   = 2000L,
  ci_levels = c(0.5, 0.8, 0.90, 0.95),
  n_perturb = 500L    # draws used in the perturbation step; reduce for speed
)

print(pred)
# ErrorTracer prediction (et_prediction)
#   Observations : 24
#   Draws        : 2000
#   CI levels    : 0.5, 0.8, 0.9, 0.95
#   Mean var decomposition (across observations):
#     Parameter    : 0.0412
#     Environmental: 0.0083
#     Residual     : 0.1204
#     Total        : 0.1619

The returned et_prediction object contains:

Step 4 — Decompose uncertainty

decompose_uncertainty() extracts the variance decomposition table from a prediction object:

decomp <- decompose_uncertainty(pred)
head(decomp)
#   obs_id param_var  env_var residual_var total_var
# 1      1    0.0389  0.0071       0.1204    0.1601
# 2      2    0.0441  0.0098       0.1204    0.1637
# ...

Decomposition method:

total_var is defined as the component sum rather than read off Var[posterior_predict] separately, so the budget reconciles to 100% exactly and every channel is a genuine sub-share of a total that contains it. (Before 1.3.0 the total was computed independently and env_var sat outside it — a defect corrected in that release.)

Autocorrelation. For models whose formula carries ar(), ma(), arma(), cosy(), unstr(), sar() or car(), a temporal_var column appears. It captures the autocorrelation-induced spread that brms::posterior_predict() accumulates iteratively beyond a single innovation, and residual_var is then read as the innovation variance rather than the stationary marginal variance. The pmax(., 0) guard on this term is instrumented: it warns when the floored shortfall exceeds twice its Monte-Carlo standard error, rather than silently absorbing a systematic negative.

Hierarchical models. Predicting for a group the fit has seen adds no channel — that group's effect is an estimated parameter and its uncertainty is already inside param_var. Predicting for a new level integrates over the group-level distribution and adds group_var, carrying the between-group variance. et_predict() detects new levels automatically and reports which budget it computed; the two are genuinely different objects and should not be compared directly.

Step 5 — Forecast shelf life

shelf_life() computes the ratio of credible interval width to the plausible response range at each forecast time point, and flags when the forecast becomes uninformative:

# For arcsin-sqrt transformed responses, derive range from training data:
#   plausible_range = range(train_df$z_diff)
# For bounded responses, supply explicit bounds:
#   plausible_range = c(lower, upper)
sl <- shelf_life(
  predictions              = pred,
  response_scale           = range(train_df$z_diff),  # defaults to the training range
  ci_level                 = 0.90,   # must be present in the et_prediction object
  threshold                = 1.0,    # CI/range ratio above which forecast is uninformative
  time_col                 = "year", # column in newdata to use as time axis
  min_run                  = 2L,     # consecutive periods required before a crossing counts
  projection_alpha         = 0.05,   # a projected horizon needs a significantly positive slope
  max_extrapolation_factor = 10      # cap on linear projection beyond observed window
)

print(sl)
# ErrorTracer shelf life analysis
#   Observations    : 24
#   Plausible range : 2
#   Informative     : 18 / 24
#   Mean CI/range   : 0.847
#   Max CI/range    : 1.231

as.data.frame(sl)
#   obs_id  time ci_width plausible_range  ratio informative
# 1      1  2024    0.841               2  0.421        TRUE
# 2      2  2025    1.103               2  0.552        TRUE
# 3      3  2026    2.461               2  1.231       FALSE

A threshold of 1.0 (default) means the forecast is considered uninformative when the CI spans the entire plausible response range. Lower thresholds (e.g., 0.8) impose stricter requirements.

Two gates guard against reporting noise as a horizon (both added in 1.3.1, both changing the default output):

• min_run (default 2) requires the crossing to persist. A ratio hovering near the threshold without trending will exceed it once by chance, and the previous first-exceedance rule promoted that single excursion to a horizon — on a real 16-period series the reported horizon moved between four different answers under sampler settings alone.
• projection_alpha (default 0.05) requires a projected horizon to rest on a slope that is significantly positive. min_slope_for_projection is a magnitude gate only, so a flat, noisy trend could be extrapolated to a confident-looking crossing time whose interval spanned centuries.

Every horizon carries n_exceedances, frac_exceedance, first_exceedance and min_run (plus slope/slope_p in projection mode) so an isolated excursion is visible rather than silent. Set min_run = 1 and projection_alpha = 1 to reproduce pre-1.3.1 behaviour. Any horizon computed with an earlier release should be recomputed, not trusted.

Step 6 — Calibrate against validation data

et_calibrate() computes observed coverage probability at each nominal CI level. A well-calibrated model produces coverage that matches the nominal level:

cal <- et_calibrate(
  predictions  = pred,
  observed     = validation_df,  # same n_rows as newdata; matched positionally
  response_col = "z_diff",       # inferred from formula if omitted
  ci_levels    = c(0.5, 0.8, 0.90, 0.95)
)

print(cal)
#   ci_level nominal observed_coverage n_obs calibration_error
# 1     0.50    0.50             0.542    24             0.042
# 2     0.80    0.80             0.792    24             0.008
# 3     0.90    0.90             0.875    24             0.025
# 4     0.95    0.95             0.958    24             0.008

Step 7 — Diagnose convergence

et_diagnose() reports Rhat, effective sample size ratios, divergent transitions, and LOO-CV in a single call:

diag <- et_diagnose(fit, loo = TRUE)

# Convergence
diag$convergence$rhat_max     # should be < 1.05
diag$convergence$neff_all_ok  # all Neff ratios > 0.1
diag$convergence$n_divergences

# LOO-CV
diag$loo$elpd_loo
diag$loo$n_bad_pareto_k  # Pareto k > 0.7 indicate influential observations

For grouped models, et_diagnose() returns a per_group list and a summary data.frame with one row per group.

Forecast evaluation

Shelf life is a statement about precision on the response scale. A forecast can be precise and still wrong, so the package pairs it with an accuracy criterion and a full-distribution calibration check.

# Null-relative forecast limit (CRPS). Two nulls, hard in opposite regimes:
#   "climatology"  (default) -- the unconditional TRAINING-response distribution:
#                   what a forecaster with no covariates would have issued. Hard to
#                   beat for a mean-reverting series. Needs no external data.
#   "random_walk"            -- persistence: the last training value carried forward
#                   with SD growing as sd(diff(train)) * sqrt(lead). Hard to beat for
#                   a near-unit-root series.
# Report BOTH: which null is demanding is a property of the series, not the model.
sk_clim <- et_skill_score(pred, observed = valid_df, time_col = "year",
                          null = "climatology")
sk_rw   <- et_skill_score(pred, observed = valid_df, time_col = "year",
                          null = "random_walk")
attr(sk_clim, "forecast_limit")   # $value, $type, $n_below, $first_below, $min_run

# Gate the shelf life on skill, so a precise-but-biased forecast cannot pass:
sl_gated <- shelf_life(pred, response_scale = range(train_df$z_diff),
                       ci_level = 0.90, skill = sk_clim)

# Probability integral transform -- the continuous-response analogue of the rank
# histogram. Uniform = calibrated; U-shape = over-confident; tilt = biased.
pit <- et_pit(pred, observed = valid_df)
et_plot_pit(pit)

# Order-independent (Sobol) decomposition, as a cross-check on the additive budget's
# ordering: reports first-order parameter and environmental shares plus their interaction.
sob <- et_sobol(pred, seed = 1L)

et_skill_score() takes the same min_run as shelf_life() (default 2) for the same reason: skill hovering near zero crosses by chance, and a bare first-crossing rule turns one unlucky lead into a forecast limit.

Pooling across forecast origins

A horizon read off a single forecast origin inherits whatever was idiosyncratic about the final training years. et_shelf_life_pool() pools horizons from repeated origins, treating origins whose data ran out before the forecast degraded as right-censored rather than dropping them — dropping them biases the pooled horizon downward, because those are exactly the long horizons.

# horizons: a data.frame with one row per forecast origin and columns
#   lead     -- the horizon in lead time
#   censored -- TRUE if the forecast was still informative when the data ran out
#   group    -- optional; pools separately per group (e.g. per species)
pooled <- et_shelf_life_pool(horizons, conf_level = 0.95)   # Kaplan-Meier via {survival}

Avoiding the post-selection problem

When the first-stage model is an elastic net, variable selection and the likelihood normally see the same rows, so reported coverage is conditional on the selected set being correct. et_priors_split() removes the issue by construction:

sp <- et_priors_split(train_df, prior_fun = my_enet_fun, prop = 0.5, seed = 1L)
fit <- et_fit(z_diff ~ Tmean + PPT + SWE, data = sp$fit_data, priors = sp$priors)

The cost is sample size: the likelihood sees fewer rows and the posterior is correspondingly wider. Split when selection stability is in doubt; don't when the predictor set is fixed a priori and there is no selection to protect against.

Visualization

All plotting functions return ggplot2 objects that can be further customized with standard ggplot2 additions (+ theme(...), + labs(...), etc.).

# Stacked bar: proportional contribution of each variance component
et_plot_decomposition(decomp, proportional = TRUE)

# Line chart: CI width / plausible range over time, with uninformative threshold
et_plot_shelf_life(sl, show_ratio = TRUE)

# Calibration: observed vs. nominal coverage (points on the 1:1 diagonal = well-calibrated)
et_plot_calibration(cal)

# Fan chart: nested CI ribbons with optional observed values overlaid
et_plot_forecast(pred, observed = validation_df,
                 response_col = "z_diff", time_col = "year")

# Density overlay: prior vs. posterior for each regression coefficient
et_plot_prior_posterior(fit, max_preds = 8L)

# Forest plot: Bayesian 95% CI (blue) vs. regularized coefficient (red ×)
et_plot_coefficients(fit)

S3 Class System

Full Example: Allele Frequency Forecasting

library(ErrorTracer)
library(glmnet)

# --- Prepare data ---
# train_df: historical years with columns z_diff, Tmean, PPT, SWE
# forecast_df: future years with predictor values only

# 1. Fit elastic net
x_mat <- as.matrix(train_df[, c("Tmean", "PPT", "SWE")])
cv_fit <- cv.glmnet(x_mat, train_df$z_diff, alpha = 0.5)

# 2. Extract informed priors
prior_spec <- extract_priors(cv_fit, multiplier = 2.0, min_sd = 0.1)

# 3. Fit Bayesian model
fit <- et_fit(
  formula = z_diff ~ Tmean + PPT + SWE,
  data    = train_df,
  priors  = prior_spec,
  chains  = 4L, iter = 2000L, cores = 4L
)

# 4. Predict with environmental noise propagation
pred <- et_predict(
  model     = fit,
  newdata   = forecast_df,
  env_noise = list(Tmean = 0.5, PPT = 10.0, SWE = 5.0),
  ci_levels = c(0.5, 0.80, 0.90, 0.95)
)

# 5. Decompose uncertainty
decomp <- decompose_uncertainty(pred)

# 6. Assess forecast shelf life
# For arcsin-sqrt responses use the observed training range; for bounded
# responses supply explicit bounds.
sl <- shelf_life(
  pred,
  response_scale           = range(train_df$z_diff),
  ci_level                 = 0.90,
  threshold                = 1.0,
  time_col                 = "year",
  max_extrapolation_factor = 10
)

# 7. Calibrate (if validation data available)
cal <- et_calibrate(pred, observed = valid_df, response_col = "z_diff")

# 8. Diagnose
et_diagnose(fit)

# 9. Visualize
library(ggplot2)
et_plot_decomposition(decomp)
et_plot_shelf_life(sl)
et_plot_calibration(cal)
et_plot_forecast(pred, observed = valid_df,
                 response_col = "z_diff", time_col = "year")
et_plot_prior_posterior(fit)
et_plot_coefficients(fit)

Group-Level Analysis Example

# Fit one model per SNP cluster, then predict and summarize shelf life by cluster
fit_list <- et_fit(
  formula  = z_diff ~ Tmean + PPT + SWE,
  data     = all_df,
  priors   = prior_spec,
  grouping = "cluster_id",
  chains = 4L, iter = 2000L, cores = 4L
)

pred_list <- et_predict(
  fit_list,
  newdata   = forecast_df,
  env_noise = list(Tmean = 0.5, PPT = 10.0),
  ci_levels = c(0.50, 0.90, 0.95)
)

# Per-group decomposition and shelf life
decomp_all <- decompose_uncertainty(pred_list)
sl_all     <- shelf_life(pred_list,
                         plausible_range          = range(train_df$z_diff),
                         ci_level                 = 0.90,
                         time_col                 = "year",
                         max_extrapolation_factor = 10)

# Per-group diagnostics summary
diag_list  <- et_diagnose(fit_list)
diag_list$summary  # data.frame: rhat_ok, neff_ok, n_divergences, elpd_loo per cluster

# Faceted shelf life plot
et_plot_shelf_life(sl_all) +
  ggplot2::facet_wrap(~ group, ncol = 4)

Comparison with Existing R Packages

ErrorTracer's niche: the full pipeline from regularized feature selection → informed Bayesian refitting → structured uncertainty decomposition → forecast shelf life, designed for ecology and genomics.

Package Structure

ErrorTracer/
├── R/
│   ├── priors.R      # extract_priors() — dispatch for glmnet, lm, glm, ranger
│   ├── fit.R         # et_fit() — Bayesian model fitting via brms
│   ├── predict.R     # et_predict() — posterior prediction with decomposition
│   ├── decompose.R   # decompose_uncertainty() — param / env / residual split
│   ├── shelf_life.R  # shelf_life() — forecast horizon analysis
│   ├── calibrate.R   # et_calibrate(), et_diagnose()
│   ├── plot.R        # et_plot_*() — six ggplot2 visualizations
│   └── utils.R       # internal helpers (standardize, logging, etc.)
├── tests/
│   └── testthat/
├── vignettes/
│   ├── genomics.Rmd  # Drosophila allele frequency example
│   └── ecology.Rmd   # Phenology / species abundance example
├── DESCRIPTION
├── NAMESPACE

Dependencies

Imports (installed automatically):

Suggests (install separately as needed):

Backend (required; install separately):

• cmdstanr (recommended) or rstan

Caveats and Limitations

• Small-n regime: When n < 10 observations, informed priors from extract_priors() dominate the posterior. This is the setting where the enet→prior pipeline provides the most value, but is also most sensitive to prior misspecification. Document prior choices carefully.
• Intercept handling: The glmnet intercept is excluded from the prior by default, as it reflects upstream centering of training data. If your data is not pre-centered, verify that the intercept prior (default Normal(0, 2.5)) is appropriate.
• Environmental uncertainty is inside the total: env_var is estimated by perturbation and is a component of total_var, which is defined as the sum of the channels. Releases before 1.3.0 reported it outside the total; decomposition tables produced by those versions should be recomputed rather than compared.
• ranger priors are undirected: Because random forests provide no signed coefficients, ranger-derived priors are centered at zero. This is conservative but avoids false directional assumptions.

License

MIT License. See LICENSE for details.

Use of AI assistance

Parts of this package were developed with the assistance of Anthropic's Claude Opus, used as a coding and review tool. Its contributions include unit and regression tests, package and function documentation, implementation work on several functions, and code audits that identified and corrected defects in the variance decomposition and in the horizon-reporting rules.

All AI-assisted output was specified, supervised, reviewed and accepted by the author, who takes full responsibility for the correctness of the package. No AI system is an author or a contributor to this software.

Authors

• Luis Javier Madrigal-Roca (Creator, Maintainer) — [email protected]
• John Kelly (Author) — [email protected]