Extract or recompute uncertainty decomposition

Description

Returns a data.frame with the three-way uncertainty decomposition stored inside an et_prediction object:

param_var

Variance of the posterior linear predictor — captures uncertainty in fitted regression coefficients.

env_var

Additional variance arising from measurement or prediction uncertainty in the predictor values (estimated via perturbation in et_predict). Zero when env_noise = NULL.

residual_var

Posterior mean of \sigma^2 — biological process noise, unmeasured drivers, and drift. Using the mean (not median) ensures param_var + residual_var \approx total_var by the law of total variance.

total_var

Variance of the full posterior predictive draws (parameter + residual; note that env_var is an additive component measured separately from the perturbation step).

Usage

decompose_uncertainty(predictions, ...)

Arguments

Details

All variance components are guaranteed non-negative.

Value

A data.frame with columns obs_id, param_var, env_var, residual_var, total_var. For grouped predictions, a group column is prepended.

Examples

set.seed(1)
df  <- data.frame(y = rnorm(20), x1 = rnorm(20))
fit <- et_fit(y ~ x1, data = df,
              chains = 1, iter = 500, warmup = 250,
              cores = 1, refresh = 0)
new_df <- data.frame(x1 = rnorm(5))
pred   <- et_predict(fit, newdata = new_df,
                     env_noise = list(x1 = 0.2),
                     n_draws = 200, n_perturb = 50)
decomp <- decompose_uncertainty(pred)
head(decomp)

Assess calibration of posterior predictive intervals

Description

Computes observed coverage probability at multiple nominal CI levels. A well-calibrated Bayesian model should produce 90% CIs that contain the true value 90% of the time, etc.

Usage

et_calibrate(predictions, observed, response_col = NULL, ci_levels = NULL, ...)

Arguments

Value

A data.frame with columns:

ci_level

Nominal CI level.

nominal

Same as ci_level.

observed_coverage

Fraction of true values falling inside the CI.

n_obs

Number of observations used.

calibration_error

Signed difference: observed - nominal. Positive = over-coverage (CIs too wide / conservative). Negative = under-coverage (CIs too narrow / overconfident).

sharpness

Mean CI width across observations. Sharpness and calibration are complementary: a model can be calibrated but useless if sharpness is poor (very wide CIs).

For grouped predictions, a group column is prepended.

Examples

set.seed(1)
df  <- data.frame(y = rnorm(20), x1 = rnorm(20))
fit <- et_fit(y ~ x1, data = df,
              chains = 1, iter = 500, warmup = 250,
              cores = 1, refresh = 0)
valid_df <- data.frame(y = rnorm(5), x1 = rnorm(5))
pred <- et_predict(fit, newdata = valid_df,
                   n_draws = 200, n_perturb = 50)
cal <- et_calibrate(pred, observed = valid_df, response_col = "y")
print(cal)

Diagnose a fitted ErrorTracer model

Description

Computes Rhat, effective sample size ratios, divergent transitions, and leave-one-out cross-validation (LOO-CV) for a fitted et_model.

Usage

et_diagnose(model, loo = TRUE, ...)

Arguments

Value

A list with elements:

convergence

List: rhat_max, rhat_all_ok, neff_min, neff_all_ok, n_divergences.

loo

List or NULL: elpd_loo, p_loo, looic, n_bad_pareto_k, loo_object.

summary

Printed summary from brms::summary().

For et_model_list, a named list of per-group diagnostic lists plus an aggregated summary data.frame.

Fit a Bayesian regression model with informed priors

Description

Wraps brms::brm() and attaches the prior specification, training data reference, and configuration for downstream uncertainty decomposition. Pass priors from extract_priors to use regularized-model coefficients as prior means; omit it for default (weakly informative) priors.

Usage

et_fit(
  formula,
  data,
  priors = NULL,
  chains = 4L,
  iter = 2000L,
  warmup = floor(iter/2),
  cores = min(chains, parallel::detectCores()),
  seed = 42L,
  adapt_delta = 0.95,
  max_treedepth = 12L,
  grouping = NULL,
  eiv = NULL,
  silent = 2L,
  ...
)

Arguments

Value

An et_model object (or an et_model_list if grouping is specified).

Examples

set.seed(1)
df  <- data.frame(y = rnorm(20), x1 = rnorm(20), x2 = rnorm(20))
ps  <- extract_priors(lm(y ~ x1 + x2, data = df))
fit <- et_fit(y ~ x1 + x2, data = df, priors = ps,
              chains = 1, iter = 500, warmup = 250,
              cores = 1, refresh = 0)
print(fit)

Plot calibration: observed vs nominal coverage

Description

A well-calibrated model produces points along the 1:1 diagonal. Points above the diagonal indicate over-coverage (conservative); below indicates under-coverage (anti-conservative).

Usage

et_plot_calibration(cal)

Arguments

Value

A ggplot2 object.

Forest plot of regression coefficients

Description

Compares Bayesian posterior estimates (95% CI) with the regularized coefficient values used as prior means (shown as crosses).

Usage

et_plot_coefficients(model)

Arguments

Value

A ggplot2 object.

Plot three-way uncertainty decomposition

Description

Produces a stacked bar chart showing the relative contributions of parameter uncertainty, environmental uncertainty, and residual variance for each observation.

Usage

et_plot_decomposition(decomp, proportional = TRUE, group_col = NULL)

Arguments

Value

A ggplot2 object.

Plot posterior predictive fan chart

Description

Overlays nested credible interval ribbons on the median forecast, with optional observed values for calibration assessment.

Usage

et_plot_forecast(
  predictions,
  observed = NULL,
  response_col = NULL,
  time_col = NULL
)

Arguments

Value

A ggplot2 object.

Plot prior vs posterior distributions for model coefficients

Description

Overlays prior and posterior density for each predictor coefficient, visualising how much the data update the priors.

Usage

et_plot_prior_posterior(model, max_preds = 8L, n_prior_draws = 4000L)

Arguments

Value

A ggplot2 object.

Plot a sensitivity profile

Description

Visualises the output of et_sensitivity_profile: for each noise grid point, shows how the environmental share of total variance and the forecast horizon respond. Observed horizons are drawn as solid points; projected horizons are hollow points; lower-bound rows are drawn as upward arrows at the last informative time.

Usage

et_plot_sensitivity(sens, show = c("horizon", "env_share", "ratio"))

Arguments

Details

The x-axis is the noise fraction (when the grid was built from fraction_grid) or the grid step label otherwise. When both a numeric fraction and a descriptive label exist, the fraction is preferred for continuous x-positioning.

Value

A ggplot2 object.

See Also

et_sensitivity_profile

Plot forecast shelf life

Description

Shows how credible interval width grows over the forecast horizon and marks the threshold beyond which the forecast is uninformative.

Usage

et_plot_shelf_life(sl, show_ratio = TRUE)

Arguments

Value

A ggplot2 object.

Posterior prediction with uncertainty decomposition

Description

Generates posterior predictive draws for new observations, propagates environmental measurement uncertainty through the model, and computes credible intervals. The resulting et_prediction object is the input to decompose_uncertainty, shelf_life, and the plotting functions.

Usage

et_predict(
  model,
  newdata,
  env_noise = NULL,
  env_cov = NULL,
  env_dist = NULL,
  n_draws = 2000L,
  ci_levels = c(0.5, 0.8, 0.9, 0.95),
  n_perturb = NULL,
  interval_type = c("predictive", "linpred"),
  include_env_in_ci = FALSE,
  ...
)

Arguments

Value

An et_prediction object (list) containing:

posterior_predict

Matrix [n_draws x n_obs]: full posterior predictive draws (parameter + residual uncertainty).

posterior_linpred

Matrix [n_draws x n_obs]: linear predictor draws (parameter uncertainty only).

lp_perturbed

Matrix [n_perturb x n_obs]: linear predictor computed on environmentally perturbed inputs.

sigma_draws

Numeric vector: posterior draws of sigma.

credible_intervals

data.frame with columns row_id, ci_level, lower, median, upper, width.

decomposition

data.frame from decompose_uncertainty(): param_var, env_var, residual_var, total_var.

newdata

The input newdata.

model

Reference to the et_model used.

env_cov

The p \times p correlation matrix actually used for perturbation (identity for env_cov = NULL).

env_dist

Named character vector mapping each predictor to the distribution actually used for its perturbation.

See Also

decompose_uncertainty, shelf_life, et_calibrate

Sensitivity profile of the environmental uncertainty component

Description

Sweeps a grid of noise magnitudes, re-runs et_predict for each, and summarises how the environmental variance component and the forecast shelf life respond. This turns the subjective choice of env_noise into an auditable curve: instead of reporting a single shelf life at one (often hand-picked) measurement-error level, the user sees how the horizon degrades as assumed noise grows.

Usage

et_sensitivity_profile(
  model,
  newdata,
  plausible_range,
  fraction_grid = NULL,
  absolute_grid = NULL,
  env_cov = NULL,
  env_dist = NULL,
  include_env_in_ci = TRUE,
  ci_level = 0.9,
  ci_levels = c(0.5, 0.8, 0.9, 0.95),
  threshold = 1,
  time_col = NULL,
  n_draws = 2000L,
  n_perturb = NULL,
  max_extrapolation_factor = 10,
  verbose = TRUE
)

Arguments

Details

Three argument styles are supported for the grid:

• fraction_grid: scalar fractions of each predictor's SD in newdata. These are passed straight to the scalar form of et_predict's env_noise. Use this for an apples- to-apples sweep that scales all predictors together.

• absolute_grid: a list of named numeric lists / vectors; each list element becomes a single env_noise argument (so you can sweep absolute SDs for one predictor while holding others fixed).

• Neither supplied (NULL): a default log-spaced fraction grid c(0, 0.01, 0.025, 0.05, 0.1, 0.2, 0.4, 0.8) is used.

Each iteration calls et_predict(..., env_noise = g) then shelf_life and records:

• Mean parameter / environmental / residual / total variance across forecast observations.

• The horizon description ("observed", "projected", or "lower_bound") and numeric horizon if any.

• Mean / max CI width / plausible-range ratio.

Value

An et_sensitivity object: a data.frame with one row per grid point and columns

grid_id

1-indexed step.

label

Short descriptor of the env_noise used (e.g.\ "fraction = 0.05").

fraction

Scalar fraction when applicable; NA for absolute-grid rows.

param_var, env_var, residual_var, total_var

Mean across forecast observations.

env_share

env_var / (env_var + total_var) — the fraction of combined predictive variance attributable to predictor noise (bounded in [0, 1]).

ci_width_mean, ci_width_max

At ci_level.

ratio_mean, ratio_max

CI width / plausible range.

horizon_type

One of "observed", "projected", "lower_bound".

horizon

Numeric horizon (observed or projected) or NA for lower-bound rows.

horizon_description

The long description returned by shelf_life.

The returned object carries the original call and plausible_range as attributes for use by et_plot_sensitivity.

See Also

et_predict, shelf_life, et_plot_sensitivity

Simulated allele-frequency-change dataset for ErrorTracer tutorials

Description

A named list containing training data, forecast predictors, validation responses, and ground-truth parameters for two synthetic SNP clusters (A and B) at a hypothetical mountain plant site. The dataset is designed to exercise every step of the ErrorTracer pipeline and to support parameter-recovery validation (true coefficients are known).

Usage

et_sim

Format

A named list with seven elements:

train

data.frame with 40 rows and 6 columns (2 clusters \times 20 training years, 1995–2014):

year

Integer. Calendar year.

cluster_id

Character. SNP cluster identifier: "A" or "B".

Tmean

Numeric. Standardised mean growing-season temperature (original unit: °C; standardised using training-period mean and SD).

PPT

Numeric. Standardised total growing-season precipitation (original unit: mm).

SWE

Numeric. Standardised peak snow water equivalent (original unit: mm).

z_diff

Numeric. Simulated allele-frequency change on the arcsin-sqrt scale (\arcsin(\sqrt{f}) transformation). This transformation is unbounded and has no fixed [-1, 1] constraint; the plausible range should be derived from the training data or from the theoretical bounds ([0,\, \pi/2]\ on the arcsin scale).

forecast

data.frame with 10 rows and 5 columns (2 clusters \times 5 forecast years, 2015–2019). Same columns as train except z_diff is absent — these are the prediction targets. Predictors are standardised using the training-period statistics stored in standardization.

validation

data.frame with 10 rows and 3 columns (year, cluster_id, z_diff). True response values for the forecast period; used with et_calibrate to assess posterior predictive coverage.

true_params

Named list with one element per cluster. Each element is a named numeric vector of the true data-generating parameters: intercept, Tmean, PPT, SWE, and sigma (residual SD).

Cluster A

intercept = 0.00, Tmean = 0.50, PPT = -0.30, SWE = 0.20, sigma = 0.40

Cluster B

intercept = 0.10, Tmean = 0.30, PPT = -0.20, SWE = -0.10, sigma = 0.50

env_noise

Named list. Suggested per-predictor environmental noise SDs (on the standardised scale) for use with the env_noise argument of et_predict: Tmean = 0.30, PPT = 0.20, SWE = 0.15.

standardization

Named list with one element per predictor. Each element is a named numeric vector c(mean = ..., sd = ...) giving the training-period statistics used to standardise the data. Needed if you wish to back-transform predictions to the original (unstandardised) scale with unstandardize.

description

Character string. Human-readable description of the dataset and how it was generated.

Details

The climate time series are simulated as AR(1) processes with a warming trend in Tmean (+0.05 °C yr^{-1}). SWE is derived from Tmean and PPT to introduce physically motivated covariance (colder winters produce more snowpack; wetter winters also increase snowpack). All predictors are standardised using training-period statistics only, which is the correct procedure when predictors in the forecast period may fall outside the training range.

Responses are generated from a linear model with Gaussian noise:

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

The true parameters are stored in et_sim$true_params for parameter-recovery validation.

The dataset was generated with set.seed(111). The full generation script is at data-raw/generate_et_sim.R.

Source

Generated by data-raw/generate_et_sim.R with set.seed(111).

Examples

data(et_sim)

# Inspect structure
str(et_sim, max.level = 2)

# Training data
head(et_sim$train)

# True parameters for cluster A
et_sim$true_params$A

# Suggested noise SDs for et_predict()
et_sim$env_noise

Minimal ggplot2 theme for ErrorTracer plots

Description

Minimal ggplot2 theme for ErrorTracer plots

Usage

et_theme(base_size = 12)

Arguments

Value

A ggplot2 theme object.

Extract brms prior specification from a regularized or standard regression model

Description

Converts a fitted model object into a brms prior specification suitable for et_fit. The coefficient estimates (or importance weights for ranger) from the regularized or standard fit are used as informative prior means, so the Bayesian model starts close to the regularized or standard solution while remaining open to data-driven revision.

Usage

extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  ...
)

## S3 method for class 'lm'
extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  ...
)

## S3 method for class 'glm'
extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  ...
)

## S3 method for class 'cv.glmnet'
extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  lambda = "lambda.min",
  ...
)

## S3 method for class 'glmnet'
extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  s = 1L,
  ...
)

## S3 method for class 'ranger'
extract_priors(
  model,
  multiplier = 2,
  min_sd = 0.1,
  intercept_prior_sd = NULL,
  sigma_prior_scale = 1,
  ...
)

Arguments

Details

For ranger models, signed coefficients are not available. Priors are centred at zero (direction unknown) and the prior SD for each predictor is set to multiplier * importance_normalised, where importance is normalised to the [min_sd, 1] interval. Only variables with positive permutation importance are included.

Value

An et_prior_spec list containing:

prior

A brmsprior object for use in brms::brm().

pred_names

Character vector of included predictor names.

coefs

Named numeric vector of regularized coefficients (NULL for ranger, which uses importance instead).

method

Character: the dispatch method used.

multiplier, min_sd

Settings echo.

Examples

fit_lm <- lm(mpg ~ wt + hp + cyl, data = mtcars)
ps <- extract_priors(fit_lm, multiplier = 2, min_sd = 0.1)
print(ps)

Compute the forecast shelf life

Description

Quantifies when a forecast becomes uninformative by comparing the width of credible intervals to a plausible response range. A forecast is uninformative when its CI width exceeds threshold * plausible_range.

Usage

shelf_life(
  predictions,
  plausible_range,
  ci_level = 0.9,
  threshold = 1,
  time_col = NULL,
  min_slope_for_projection = 1e-04,
  max_extrapolation_factor = 10,
  ...
)

Arguments

Details

The function operates in three modes depending on the available data and whether the uninformative threshold is crossed within the forecast window:

Observed

The threshold is crossed within the forecast/validation window. The shelf life is the first time point at which ratio >= threshold.

Projected

All forecast periods remain informative but the CI/range ratio is trending upward. A linear trend is fitted to the ratios and extrapolated to estimate when the threshold would be reached.

Lower bound

All forecast periods are informative with no upward trend in the ratio. The shelf life is reported as a lower bound: > last observed time.

The intended workflow is:

1. Fit the model (et_fit).

2. Predict on held-out data or a future time window (et_predict).

3. Call shelf_life() on those predictions.

4. If the threshold is not crossed in the held-out window, the projected mode extrapolates the horizon automatically.

Value

An et_shelf_life object (a data.frame) with columns:

obs_id

Observation index.

time

Time axis value.

ci_width

Width of the credible interval at ci_level.

plausible_range

Effective plausible range (scalar).

ratio

CI width / plausible range.

informative

Logical; TRUE when ratio < threshold.

For grouped predictions a group column is prepended. The object carries three attributes that drive print():

horizon

Named list with elements value, type ("observed", "projected", or "lower_bound"), last_informative, and description.

horizon_by_group

For grouped objects: named list of per-group horizon lists.

threshold

The threshold value used.

Examples

set.seed(1)
df  <- data.frame(y = rnorm(20), year = 2001:2020, x1 = rnorm(20))
fit <- et_fit(y ~ x1, data = df,
              chains = 1, iter = 500, warmup = 250,
              cores = 1, refresh = 0)
new_df <- data.frame(x1 = rnorm(10), year = 2021:2030)
pred   <- et_predict(fit, newdata = new_df,
                     n_draws = 200, n_perturb = 50)
sl <- shelf_life(pred,
                 plausible_range          = range(df$y),
                 ci_level                 = 0.90,
                 threshold                = 1.0,
                 time_col                 = "year",
                 max_extrapolation_factor = 10)
print(sl)

Standardize a numeric vector (mean-centre, unit-variance)

Description

Standardize a numeric vector (mean-centre, unit-variance)

Usage

standardize(x)

Arguments

Value

Standardized numeric vector. Returns zeros if variance is zero.

Reverse standardization

Description

Reverse standardization

Usage

unstandardize(z, mu, s)

Arguments

Value

Values on the original scale.