Model averaging for tulpaRatio fits

Compute model-averaged predictions using stacking or pseudo-BMA weights. Model weights are derived from LOO-CV or WAIC, and predictions are combined accordingly.

Usage

ratiod_average(
  ...,
  weights = c("loo", "waic", "pbma", "pbma+"),
  newdata = NULL,
  type = c("ratio", "numerator", "denominator"),
  summary = TRUE
)

Arguments

...

Multiple ratiod_fit objects to average

weights

Method for computing weights: "loo" (stacking, default), "waic", "pbma" (pseudo-BMA), or "pbma+" (pseudo-BMA+ with Bayesian bootstrap).

newdata

Optional new data for predictions. If NULL, uses fitted values.

type

Type of prediction: "ratio" (default), "numerator", or "denominator".

summary

Logical; if TRUE (default), return summary statistics. If FALSE, return full posterior draws.

Value

A ratiod_average object containing:

• weights: Model weights

• predictions: Model-averaged predictions

• models: List of model names

Details

Model averaging accounts for model uncertainty by combining predictions from multiple candidate models. The weighting methods are:

• stacking (weights = "loo"): Optimal linear combination minimizing leave-one-out prediction error. Recommended default.

• pseudo-BMA (weights = "pbma"): Weights proportional to exp(-ELPD). Can be unstable with similar models.

• pseudo-BMA+ (weights = "pbma+"): Bayesian bootstrap variant that accounts for estimation uncertainty.

For ratio predictions, averaging is performed on the log scale to respect the multiplicative nature of ratios.

See also

ratiod_compare() for model comparison, loo::stacking_weights()

Examples

# \donttest{
# Fit candidate models (slow, not run on CRAN)
set.seed(123)
n <- 60
df <- data.frame(
  count = rpois(n, lambda = 8),
  effort = rgamma(n, shape = 4, rate = 1),
  x = rnorm(n),
  z = rnorm(n),
  site = factor(rep(1:10, each = 6))
)
fit1 <- tratio(count | effort ~ x + (1 | site), data = df,
              family = ratiod_poisson_gamma(),
              control = list(iter = 200, warmup = 100, chains = 1))
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#>   Reason: default (full MCMC)
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 60
#>   Iterations: 200 (warmup: 100)
#> Running NUTS sampler...
#>   Parameters: 16
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 1)
#> Warning: 9 divergent transition(s) after warmup. Increase max_treedepth or reparameterize.
fit2 <- tratio(count | effort ~ x + z + (1 | site), data = df,
              family = ratiod_poisson_gamma(),
              control = list(iter = 200, warmup = 100, chains = 1))
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#>   Reason: default (full MCMC)
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 60
#>   Iterations: 200 (warmup: 100)
#> Running NUTS sampler...
#>   Parameters: 18
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 1)
# Model averaging (requires loo package)
# avg <- ratiod_average(fit1, fit2)
# print(avg)
# }