Package {ErrorTracer}

Extract or recompute uncertainty decomposition

Description

Returns a data.frame with the 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 (or its family-specific analogue) — biological process noise, unmeasured drivers, and drift. For autocorrelation models this is the innovation variance, not the stationary marginal variance; the autocorrelated accumulation is reported separately in temporal_var.

temporal_var

(Only present when the model formula contains an autocorrelation term such as ar(), ma(), arma(), cosy(), unstr(), sar(), or car().) Variance attributable to residual temporal or spatial dependence beyond the iid param + residual sum, computed as pmax(0, total_var - (param_var + residual_var)). env_var is deliberately excluded from this gap because it is an additive perturbation-based augmentation measured outside of posterior_predict.

total_var

Variance of the full posterior predictive draws, including any autocorrelation structure modelled by brms.

Usage

decompose_uncertainty(predictions, ...)

Arguments

Details

All variance components are guaranteed non-negative. When temporal_var is present, param_var + residual_var + temporal_var reconstructs total_var (modulo Monte Carlo error); when it is absent, param_var + residual_var does. env_var is always additive on top, representing the extra variance that would be contributed by perturbing predictors with env_noise.

Value

A data.frame with columns obs_id, param_var, env_var, residual_var, total_var (plus temporal_var for autocorrelation models, and a leading group column for grouped predictions).

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, group_col = NULL)

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 uncertainty decomposition

Description

Produces a stacked bar chart showing the relative contributions of parameter, environmental, and residual variance for each observation, plus a fourth temporal-autocorrelation component when present in decomp (see decompose_uncertainty).

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,
  n_env_draws = 1L,
  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 on the link scale (parameter uncertainty only).

lp_perturbed

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

sigma_draws

Numeric vector: posterior draws of sigma (NA for families without a sigma parameter, e.g.\ Binomial).

credible_intervals

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

decomposition

data.frame with columns obs_id, total_var, param_var, env_var, v_env_mcse, residual_var. All components are on the response scale. v_env_mcse is the Monte Carlo SE of env_var. residual_var is per-observation for non-Gaussian families (e.g.\ Binomial: E_s[\mu^{(s)}(1-\mu^{(s)})]) and constant for Gaussian.

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,
  response_scale,
  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,
  plausible_range = NULL
)

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 response_scale 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 100 rows and 6 columns (2 clusters \times 50 training years, 1970–2019):

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 30 rows and 5 columns (2 clusters \times 15 forecast years, 2020–2034). 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 30 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.30

Cluster B

intercept = 0.10, Tmean = 0.30, PPT = -0.20, SWE = -0.25, sigma = 0.35

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.015 °C yr^{-1}, cumulating to ~0.30 SD above the training mean over the 15-year forecast window). SWE is generated with a weak dependence on Tmean (negative coupling) and PPT (positive coupling) plus a dominant independent noise component, so that the three predictors have only mild pairwise correlation (|R| < 0.2 on the training subset). This makes all regression coefficients reliably identifiable in a single replicate. All predictors are standardised using training-period statistics only.

In addition to the three real climate predictors, the dataset includes ten independent standardised nuisance predictors d1, ..., d10 with zero true effect on the response. They give the regularized-prior extraction step a real selection job: cv.glmnet (alpha = 0.5) is expected to identify the three real predictors and shrink the ten dummies toward zero. Without them, regularization on three well-identified predictors merely biases the true coefficients without selecting anything.

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 response scale. A forecast is uninformative when its CI width exceeds threshold * response_scale.

Usage

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

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. The projected crossing time t^* = (\tau - a) / b (where \tau is the threshold, a the fitted intercept, b the fitted slope) is reported together with a Monte Carlo standard error se_t_star derived via the delta method.

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 response scale (scalar diff).

ratio

CI width / response scale.

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, description, and (for projected) se_t_star.

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,
                 response_scale            = 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.