--- title: "Getting Started with ErrorTracer" subtitle: "Bayesian error propagation and forecast uncertainty decomposition" author: "Luis Javier Madrigal-Roca, John Kelly" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true vignette: > %\VignetteIndexEntry{Getting Started with ErrorTracer} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE, # Stan/brms chunks are too slow for automated builds fig.width = 7, fig.height = 4.5 ) ``` > **How to use this vignette.** > All interactive code chunks have `eval = FALSE` so the vignette can build > without running Stan. Copy and run each chunk in sequence in your own R > session. > Before the first use, generate the bundled dataset by running once from the > package root: > > ```r > source("data-raw/generate_et_sim.R") > ``` > > For real analyses use `chains = 4` and `iter` ≥ 2000. The demo below uses > 2 chains × 1 000 iterations to keep runtime short. --- # Introduction **ErrorTracer** addresses a gap in the ecological and genomic forecasting toolkit: existing Bayesian packages (`brms`, `rstanarm`) do not decompose forecast uncertainty into its contributing sources, and no R package quantifies the *forecast shelf life* — the time horizon at which a credible interval becomes too wide to be informative. The package provides a seven-step pipeline: | Step | Function | What it does | |---|---|---| | 1 | `extract_priors()` / `et_priors_split()` | Convert a regularized or standard model to `brms` priors | | 2 | `et_fit()` / `as_et_model()` | Fit the informed Bayesian model (or wrap an existing `brms` fit) | | 3 | `et_diagnose()` | Convergence, ESS, LOO-CV | | 4 | `et_predict()` | Posterior prediction + environmental noise / ensemble propagation | | 5 | `decompose_uncertainty()` / `et_sobol()` | Variance decomposition (param / group / env / residual / temporal) & Sobol indices | | 6 | `shelf_life()` / `et_shelf_life_pool()` | Forecast horizon / shelf-life analysis & rolling-origin pooling | | 7 | `et_calibrate()`, `et_pit()`, `et_skill_score()` | Coverage, PIT calibration, and CRPS skill against nulls | Visualisation functions (`et_plot_*`) return `ggplot2` objects at every step. ## 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. | Predictor | Description | Units (raw) | |---|---|---| | `Tmean` | Mean growing-season temperature | °C (standardised) | | `PPT` | Total growing-season precipitation | mm (standardised) | | `SWE` | Peak snow water equivalent | mm (standardised) | | `d1`..`d10` | Ten independent nuisance predictors (zero true effect) | standardised | 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) $$ | Parameter | Cluster A | Cluster B | |---|---|---| | $\alpha$ (intercept) | 0.00 | 0.10 | | $\beta_1$ (Tmean) | **+0.50** | **+0.30** | | $\beta_2$ (PPT) | **−0.30** | **−0.20** | | $\beta_3$ (SWE) | **+0.20** | **−0.25** | | $\sigma$ | 0.30 | 0.35 | 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. --- # Setup ```{r libraries} library(ErrorTracer) library(glmnet) # for elastic net (Suggests) library(ggplot2) # Load the bundled simulated dataset data(et_sim) ``` ```{r demo-settings} # Use small MCMC settings for this tutorial. # In a real analysis: CHAINS = 4, ITER = 2000 (or more). CHAINS <- 2L ITER <- 1000L SEED <- 42L ``` --- ```{r data-setup, include = FALSE, eval = TRUE} library(ErrorTracer) data(et_sim) ``` # Exploring the dataset ```{r data-overview, eval = TRUE} str(et_sim, max.level = 2) ``` ```{r head-train, eval = TRUE} head(et_sim$train[, 1:6]) ``` ```{r true-params, eval = TRUE} # True coefficients — known because the data are simulated et_sim$true_params ``` ```{r forecast-predictors, eval = TRUE} # Forecast predictors sit above the training range in Tmean (warming trend) head(et_sim$forecast[et_sim$forecast$cluster_id == "A", 1:5]) ``` ## Visualise the training data ```{r plot-training-data, fig.cap="Simulated allele-frequency change over the training period."} 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() ``` ```{r plot-predictors, fig.cap="Standardised climate predictors over the training and forecast periods."} # 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() ``` --- # Single-cluster workflow: Cluster A We walk through the complete pipeline for cluster A, then repeat for both clusters simultaneously using the `grouping` argument. ## 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. ```{r enet-fit} 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). ```{r enet-plot, fig.cap="Cross-validated elastic net for cluster A."} plot(cv_fit_A) ``` ## 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). ```{r extract-priors} 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))`). ## Step 3 — Fit the Bayesian model `et_fit()` wraps `brms::brm()`, attaches the prior specification, and returns an `et_model` object. ```{r et-fit} 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 ``` ```{r et-summary} 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. ## Step 4 — Diagnose convergence and model fit ```{r et-diagnose} 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 ``` ```{r loo-results} 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. ## 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. ```{r et-predict} 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. ```{r ci-table} # 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 ``` ## Step 6 — Uncertainty decomposition ```{r decompose} 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)`. ```{r plot-decomp, fig.cap="Variance decomposition for cluster A over the 2020-2034 forecast window."} 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. ## 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`). ```{r shelf-life} # 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) ``` ```{r shelf-life-table} 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. ```{r plot-shelf-life, fig.cap="Shelf life for cluster A: 90% CI width divided by the training range, on the arcsin-square-root scale."} et_plot_shelf_life(sl_A, show_ratio = TRUE) ``` ## Step 8 — Calibration assessment We compare the posterior predictive intervals against the held-out true values in `et_sim$validation`. ```{r calibrate} 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. ```{r plot-calibration, fig.cap="Posterior predictive calibration for cluster A."} et_plot_calibration(cal_A) ``` --- # Visualisations ## Forecast fan chart ```{r plot-forecast, fig.cap="Forecast for cluster A. Shaded ribbons are the 50, 80, 90 and 95% credible intervals."} 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. ## Prior versus posterior ```{r plot-prior-posterior, fig.cap="How 20 observations update the elastic-net prior, cluster A."} 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. ## Forest plot: Bayesian vs. elastic net ```{r plot-coefficients, fig.cap="Bayesian 95% credible intervals (blue) against the elastic-net coefficients (red crosses), cluster A."} 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. --- # Parameter recovery validation A key advantage of simulated data is that we can check whether the 95% posterior credible intervals contain the known true values. ```{r parameter-recovery} post <- brms::fixef(fit_A$fit) true_A <- et_sim$true_params$A post_df <- data.frame( parameter = c("Intercept", "Tmean", "PPT", "SWE"), true_value = c(true_A["intercept"], true_A["Tmean"], true_A["PPT"], true_A["SWE"]), post_mean = post[, "Estimate"], lower_95 = post[, "Q2.5"], upper_95 = post[, "Q97.5"] ) post_df$covered <- with(post_df, true_value >= lower_95 & true_value <= upper_95) post_df[, c("parameter", "true_value", "post_mean", "lower_95", "upper_95", "covered")] # parameter true_value post_mean lower_95 upper_95 covered # 1 Intercept 0.00 0.041 -0.162 0.246 TRUE # 2 Tmean 0.50 0.666 0.438 0.888 TRUE # 3 PPT -0.30 -0.487 -0.909 -0.123 TRUE # 4 SWE 0.20 0.406 0.023 0.827 TRUE ``` All four true parameter values fall within their 95% posterior credible intervals. The posterior means for Tmean, PPT, and SWE are consistently above the true values in absolute magnitude — this is typical over-estimation bias with n = 20 and informative priors near the OLS estimate. The bias is small relative to the interval widths, and would shrink with more training data. --- # Multi-cluster grouped workflow In practice you often need to fit the same model to many units (SNP clusters, species, sites) simultaneously. `et_fit()` handles this with the `grouping` argument, returning an `et_model_list`. ## Fitting ```{r grouped-priors} # 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") ``` ```{r grouped-fit} # 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 ``` ## Diagnostics ```{r grouped-diagnose} 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. ## Prediction and decomposition ```{r grouped-predict} 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 ``` ```{r grouped-decompose} 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). ## Shelf life comparison ```{r grouped-shelf-life} # 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. ```{r plot-grouped-shelf-life, fig.cap="Shelf life for both clusters: 90% CI width divided by the plausible range."} et_plot_shelf_life(sl_all, show_ratio = TRUE) + labs( colour = "Cluster" ) ``` ## Calibration and forecast plots ```{r grouped-calibrate, fig.cap="Posterior predictive calibration for both clusters."} 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) ``` ```{r grouped-coef-plot, fig.cap="Bayesian 95% credible intervals (blue) against the elastic-net coefficients (red crosses), both clusters."} # Forest plot for all groups in a single faceted figure et_plot_coefficients(fit_all) ``` --- # Long-horizon forecasting with GCM projections The shelf-life analysis becomes most powerful when combined with climate projections from Global Circulation Models (GCMs). The idea is: 1. **Validate** that the model works on real held-out years (e.g.\ 2022–2025). 2. **Project** onto GCM-simulated climate for 2025–2075 to ask: *over what horizon does the forecast remain informative?* A critical refinement is that GCM uncertainty is **not constant** — ensemble spread (the disagreement among model runs) typically grows with the projection horizon. The `env_noise` argument now accepts **per-row vectors** so you can pass growing noise SDs directly. ```{r gcm-workflow, eval = FALSE} # Suppose you have GCM ensemble-mean predictors for 2025-2075, # standardised using your training-period statistics. gcm_years <- 2025:2075 gcm_df <- data.frame( year = gcm_years, Tmean = ..., # standardised GCM ensemble-mean temperature PPT = ..., # standardised GCM ensemble-mean precipitation SWE = ... # standardised GCM ensemble-mean snowpack ) # GCM ensemble spread (SD across runs) grows roughly linearly with time. # Here we model spread as starting at the validation-period noise level # and increasing by ~0.01 SD per year for temperature. base_year <- 2025 tmean_noise <- 0.30 + 0.01 * (gcm_years - base_year) # 0.30 → 0.80 ppt_noise <- 0.20 + 0.005 * (gcm_years - base_year) # 0.20 → 0.45 swe_noise <- 0.15 + 0.005 * (gcm_years - base_year) # 0.15 → 0.40 # Predict over the 50-year GCM window with time-varying noise pred_gcm <- et_predict( model = fit_A, newdata = gcm_df, env_noise = list( Tmean = tmean_noise, # per-row vector, length = nrow(gcm_df) PPT = ppt_noise, SWE = swe_noise ), ci_levels = c(0.90, 0.95), n_draws = 1000L ) # Shelf life over the 50-year horizon sl_gcm <- shelf_life( pred_gcm, response_scale = range(train_A$z_diff), time_col = "year", max_extrapolation_factor = Inf # allow projection beyond the 50-yr window ) print(sl_gcm) # If the 90% CI stays below the plausible range throughout 2025-2075: # Shelf life : > 2075 (lower bound — all periods informative, no trend) # If it crosses (say, around 2058): # Shelf life : ~2058 (observed — threshold first exceeded) ``` The per-row noise means that `decompose_uncertainty()` now shows `env_var` **increasing** across forecast years — directly capturing the epistemic uncertainty from the GCM ensemble spread. ```{r gcm-decompose, eval = FALSE} decomp_gcm <- decompose_uncertainty(pred_gcm) # env_var should be small in 2025 and growing toward 2075 as GCM spread # accumulates; param_var and residual_var remain roughly constant. plot(gcm_years, decomp_gcm$env_var, type = "l", xlab = "Year", ylab = "Environmental variance", main = "Growing GCM ensemble uncertainty propagated into forecast") ``` > **Interpreting the result for conservation.** > If the shelf life exceeds your planning horizon (e.g.\ "> 2075"), you can > state: *"Under the [GCM scenario], the 90% credible interval for allele- > frequency change remains within the plausible biological range throughout > the 50-year planning horizon, supporting its use in long-term conservation > genomic forecasting."* > If a crossing year is observed (e.g.\ 2058), you can report the horizon > directly and note which uncertainty source (parameter, environmental, or > residual) drives the loss of informativeness. --- # Using `lm` instead of `glmnet` `extract_priors()` is not limited to elastic net. With the default `shrinkage = "zero"`, the `lm` method scales prior SDs by coefficient magnitude (`multiplier * |coef|`, floored at `min_sd`); with `shrinkage = "estimate"` it uses OLS standard errors (`multiplier * SE`) around the OLS point estimates. ```{r lm-priors} lm_fit_A <- lm(z_diff ~ Tmean + PPT + SWE, data = train_A) prior_lm_A <- extract_priors( model = lm_fit_A, multiplier = 2.0, min_sd = 0.10 ) print(prior_lm_A) # ErrorTracer prior specification # Method : lm # Predictors : 3 # Multiplier : 2 # Min SD : 0.1 # Shrinkage : zero (prior means at 0; regularizing) # Coefficients: # Tmean prior mean = 0.0000 sd = 1.0023 (coef = 0.5012) # PPT prior mean = 0.0000 sd = 0.5490 (coef = -0.2745) # SWE prior mean = 0.0000 sd = 0.3416 (coef = 0.1708) ``` In this case the OLS and elastic-net estimates on the three true predictors are very close, so the two prior specifications lead to nearly identical posteriors. The key advantage of the elastic-net pathway arises when candidate nuisance predictors (`d1`..`d10`) are included and regularisation is needed to screen them out. --- # Standardising and back-transforming The `standardization` slot stores the training-period mean and SD for each predictor, enabling back-transformation to original units with `unstandardize()`. ```{r back-transform} # Standardisation constants used when generating et_sim et_sim$standardization # $Tmean # mean sd # 15.74798 1.24258 # # $PPT # mean sd # 108.60625 15.94841 # # $SWE # mean sd # 61.23573 20.64495 # A Tmean of +1 SD corresponds to: unstandardize(z = 1.0, mu = et_sim$standardization$Tmean["mean"], s = et_sim$standardization$Tmean["sd"]) # [1] 16.99 ← 1 SD above training mean = 17.0 °C ``` --- # Summary This vignette has demonstrated the complete ErrorTracer pipeline on simulated data where the ground truth is known: 1. **`extract_priors()`** converted elastic-net coefficients into `brms`-compatible priors (centred at zero by default with `shrinkage = "zero"` to avoid prior/likelihood double-use, or at the point estimates with `shrinkage = "estimate"`). 2. **`et_fit()`** fitted the Bayesian model for a single cluster and for both clusters simultaneously via `grouping = "cluster_id"` (or an existing `brmsfit` can be passed directly via `as_et_model()`). 3. **`et_diagnose()`** confirmed convergence (Rhat < 1.05, zero divergences) and adequate model fit (no high Pareto-k observations for cluster A). 4. **`et_predict()`** generated reproducible posterior predictive draws and propagated environmental measurement uncertainty into the credible intervals and total predictive variance. 5. **`decompose_uncertainty()`** partitioned measured forecast variance (`total_var`) into parameter, environmental, and residual components, identifying residual process noise as the dominant bottleneck in both clusters. 6. **`shelf_life()`** quantified that cluster A forecasts are more informative than cluster B, consistent with cluster B's higher residual variance. Both clusters remain well below the uninformative threshold (1.0) across the forecast window — reported as lower bounds. The three-mode horizon logic (`observed` / `projected` / `lower_bound`), persistence requirement (`min_run = 2`), trend significance test (`projection_alpha = 0.05`), and `max_extrapolation_factor` cap prevent spurious horizons or far-future projections when the CI/range trend is flat. 7. **`et_calibrate()`** confirmed near-nominal posterior predictive coverage on held-out validation data. 8. **Parameter recovery** confirmed that all four true coefficients fell within their 95% posterior credible intervals, validating the pipeline's statistical integrity. --- # Session information ```{r session-info, eval = TRUE, echo = FALSE} sessionInfo() ```