mvgam (version 1.1.594)

loo.mvgam: LOO information criteria for mvgam models

Description

Extract the LOOIC (leave-one-out information criterion) using loo::loo().

Usage

# S3 method for mvgam
loo(x, incl_dynamics = FALSE, ...)
# S3 method for mvgam
loo_compare(x, ..., model_names = NULL, incl_dynamics = FALSE)

Value

For loo.mvgam, an object of class psis_loo (see loo::loo()

for details). For loo_compare.mvgam, an object of class compare.loo

(see loo::loo_compare() for details).

Arguments

x

Object of class mvgam

incl_dynamics

Deprecated and currently ignored

...

More mvgam objects

model_names

If NULL (the default) will use model names derived from deparsing the call. Otherwise will use the passed values as model names

Author

Nicholas J Clark

Details

When comparing two (or more) fitted mvgam models, we can estimate the difference in their in-sample predictive accuracies using the Expected Log Predictive Density (ELPD). This metric can be approximated using Pareto Smoothed Importance Sampling (PSIS), which re-weights posterior draws to approximate predictions for a datapoint had it not been included in the original model fit (i.e. leave-one-out cross-validation).

See loo::loo() and loo::loo_compare() for further details on how this importance sampling works.

Note: In-sample predictive metrics such as PSIS-LOO can sometimes be overly optimistic for models that include process error components (e.g. those with trend_model, trend_formula, or factor_formula). Consider using out-of-sample evaluations for further scrutiny (see forecast.mvgam, score.mvgam_forecast, lfo_cv).

Examples

if (FALSE) {
#--------------------------------------------------
# Simulate 4 time series with hierarchical seasonality
# and independent AR1 dynamic processes
#--------------------------------------------------
set.seed(111)

simdat <- sim_mvgam(
  seasonality = 'hierarchical',
  trend_model = AR(),
  family = gaussian()
)

# Fit a model with shared seasonality
mod1 <- mvgam(
  y ~ s(season, bs = 'cc', k = 6),
  data = rbind(simdat$data_train, simdat$data_test),
  family = gaussian(),
  chains = 2,
  silent = 2
)

conditional_effects(mod1)

mc.cores.def <- getOption('mc.cores')
options(mc.cores = 1)
loo(mod1)

# Fit a model with hierarchical seasonality
mod2 <- update(
  mod1,
  formula = y ~ s(season, bs = 'cc', k = 6) +
    s(season, series, bs = 'fs', xt = list(bs = 'cc'), k = 4),
  chains = 2,
  silent = 2
)

conditional_effects(mod2)
loo(mod2)

# Add AR1 dynamic errors to mod2
mod3 <- update(
  mod2,
  trend_model = AR(),
  chains = 2,
  silent = 2
)

conditional_effects(mod3)
plot(mod3, type = 'trend')
loo(mod3)

#--------------------------------------------------
# Compare models using LOO
#--------------------------------------------------
loo_compare(mod1, mod2, mod3)
options(mc.cores = mc.cores.def)

#--------------------------------------------------
# Compare forecast abilities using LFO-CV
#--------------------------------------------------

lfo_mod2 <- lfo_cv(mod2, min_t = 92)
lfo_mod3 <- lfo_cv(mod3, min_t = 92)

# Plot forecast ELPD differences
plot(
  y = lfo_mod2$elpds - lfo_mod3$elpds,
  x = lfo_mod2$eval_timepoints,
  pch = 16,
  ylab = 'ELPD_mod2 - ELPD_mod3',
  xlab = 'Evaluation timepoint'
)

abline(h = 0, lty = 'dashed')
}