Package 'mvgam'

Calculate randomized quantile residuals for mvgam objects

Description

Calculate randomized quantile residuals for mvgam objects

Usage

add_residuals(object, ...)

## S3 method for class 'mvgam'
add_residuals(object, ...)

Arguments

object	list object of class mvgam. See mvgam()
...	unused

Details

For each series, randomized quantile (i.e. Dunn-Smyth) residuals are calculated for inspecting model diagnostics If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent draws from the model's posterior distribution

Value

A list object of class mvgam with residuals included in the 'resids' slot

---

NEON Amblyomma and Ixodes tick abundance survey data

Description

A dataset containing timeseries of Amblyomma americanum and Ixodes scapularis nymph abundances at NEON sites

Usage

all_neon_tick_data

Format

A tibble/dataframe containing covariate information alongside the main fields of:

Year

Year of sampling

epiWeek

Epidemiological week of sampling

plot_ID

NEON plot ID for survey location

siteID

NEON site ID for survey location

amblyomma_americanum

Counts of A. americanum nymphs

ixodes_scapularis

Counts of I. scapularis nymphs

Source

https://www.neonscience.org/data

---

Augment an mvgam object's data

Description

Add fits and residuals to the data, implementing the generic augment from the package broom.

Usage

## S3 method for class 'mvgam'
augment(x, robust = FALSE, probs = c(0.025, 0.975), ...)

Arguments

x	An object of class mvgam.
robust	If FALSE (the default) the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and the median absolute deviation (MAD) are applied instead.
probs	The percentiles to be computed by the quantile function.
...	Unused, included for generic consistency only.

Details

A list is returned if class(x$obs_data) == 'list', otherwise a tibble is returned, but the contents of either object is the same.

The arguments robust and probs are applied to both the fit and residuals calls (see fitted.mvgam() and residuals.mvgam() for details).

Value

A list or tibble (see details) combining:

• The data supplied to mvgam().

• The outcome variable, named as .observed.

• The fitted backcasts, along with their variability and credible bounds.

• The residuals, along with their variability and credible bounds.

See Also

residuals.mvgam, fitted.mvgam

Other tidiers: tidy.mvgam()

Examples

## Not run: 
set.seed(0)
dat <- sim_mvgam(
  T = 80,
  n_series = 3,
  mu = 2,
  trend_model = AR(p = 1),
  prop_missing = 0.1,
  prop_trend = 0.6
)

mod1 <- mvgam(
  formula = y ~ s(season, bs = 'cc', k = 6),
  data = dat$data_train,
  trend_model = AR(),
  family = poisson(),
  noncentred = TRUE,
  chains = 2,
  silent = 2
)

augment(mod1, robust = TRUE, probs = c(0.25, 0.75))

## End(Not run)

---

Stan code and data objects for mvgam models

Description

Generate Stan code and data objects for mvgam models

Usage

code(object)

## S3 method for class 'mvgam_prefit'
stancode(object, ...)

## S3 method for class 'mvgam'
stancode(object, ...)

## S3 method for class 'mvgam_prefit'
standata(object, ...)

Arguments

object	An object of class mvgam or mvgam_prefit, returned from a call to mvgam
...	ignored

Value

Either a character string containing the fully commented Stan code to fit a mvgam model or a named list containing the data objects needed to fit the model in Stan.

Examples

## Not run: 
simdat <- sim_mvgam()
mod <- mvgam(y ~ s(season) +
               s(time, by = series),
             family = poisson(),
             data = simdat$data_train,
             run_model = FALSE)

# View Stan model code
stancode(mod)

# View Stan model data
sdata <- standata(mod)
str(sdata)

## End(Not run)

---

Display conditional effects of predictors for mvgam models

Description

Display conditional effects of one or more numeric and/or categorical predictors in models of class mvgam and jsdgam, including two-way interaction effects.

Usage

## S3 method for class 'mvgam'
conditional_effects(
  x,
  effects = NULL,
  type = "expected",
  points = FALSE,
  rug = FALSE,
  ...
)

## S3 method for class 'mvgam_conditional_effects'
plot(x, plot = TRUE, ask = FALSE, ...)

## S3 method for class 'mvgam_conditional_effects'
print(x, ...)

Arguments

x	Object of class mvgam, jsdgam or mvgam_conditional_effects
effects	An optional character vector naming effects (main effects or interactions) for which to compute conditional plots. Interactions are specified by a : between variable names. If NULL (the default), plots are generated for all main effects and two-way interactions estimated in the model. When specifying effects manually, all two-way interactions (including grouping variables) may be plotted even if not originally modeled.
type	character specifying the scale of predictions. When this has the value link the linear predictor is calculated on the link scale. If expected is used (the default), predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. When response is used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. Two special cases are also allowed: type latent_N will return the estimated latent abundances from an N-mixture distribution, while type detection will return the estimated detection probability from an N-mixture distribution.
points	Logical. Indicates if the original data points should be added, but only if type == 'response'. Default is TRUE.
rug	Logical. Indicates if displays tick marks should be plotted on the axes to mark the distribution of raw data, but only if type == 'response'. Default is TRUE.
...	other arguments to pass to plot_predictions
plot	Logical; indicates if plots should be plotted directly in the active graphic device. Defaults to TRUE.
ask	Logical. Indicates if the user is prompted before a new page is plotted. Only used if plot is TRUE. Default is FALSE.

Details

This function acts as a wrapper to the more flexible plot_predictions. When creating conditional_effects for a particular predictor (or interaction of two predictors), one has to choose the values of all other predictors to condition on. By default, the mean is used for continuous variables and the reference category is used for factors. Use plot_predictions to change these and create more bespoke conditional effects plots.

Value

conditional_effects returns an object of class mvgam_conditional_effects which is a named list with one slot per effect containing a ggplot object, which can be further customized using the ggplot2 package. The corresponding plot method will draw these plots in the active graphic device.

Author(s)

Nicholas J Clark

See Also

plot_predictions, plot_slopes

Examples

## Not run: 
# Simulate some data
simdat <- sim_mvgam(
  family = poisson(),
  seasonality = 'hierarchical'
)

# Fit a model
mod <- mvgam(
  y ~ s(season, by = series, k = 5) + year:series,
  family = poisson(),
  data = simdat$data_train,
  chains = 2,
  silent = 2
)

# Plot all main effects on the response scale
conditional_effects(mod)

# Change the prediction interval to 70% using plot_predictions() argument
# 'conf_level'
conditional_effects(mod, conf_level = 0.7)

# Plot all main effects on the link scale
conditional_effects(mod, type = 'link')

# Works the same for smooth terms, including smooth interactions
set.seed(0)
dat <- mgcv::gamSim(1, n = 200, scale = 2)
mod <- mvgam(
  y ~ te(x0, x1, k = 5) + s(x2, k = 6) + s(x3, k = 6),
  data = dat,
  family = gaussian(),
  chains = 2,
  silent = 2
)
conditional_effects(mod)
conditional_effects(mod, conf_level = 0.5, type = 'link')

# ggplot objects can be modified and combined with the help of many
# additional packages. Here is an example using the patchwork package

# Simulate some nonlinear data
dat <- mgcv::gamSim(1, n = 200, scale = 2)
mod <- mvgam(
  y ~ s(x1, bs = 'moi') + te(x0, x2),
  data = dat,
  family = gaussian(),
  chains = 2,
  silent = 2
)

# Extract the list of ggplot conditional_effect plots
m <- plot(conditional_effects(mod), plot = FALSE)

# Add custom labels and arrange plots together using patchwork::wrap_plots()
library(patchwork)
library(ggplot2)
wrap_plots(
  m[[1]] + labs(title = 's(x1, bs = "moi")'),
  m[[2]] + labs(title = 'te(x0, x2)')
)

## End(Not run)

---

Defining dynamic coefficients in mvgam formulae

Description

Set up time-varying (dynamic) coefficients for use in mvgam models. Currently, only low-rank Gaussian Process smooths are available for estimating the dynamics of the time-varying coefficient.

Usage

dynamic(variable, k, rho = 5, stationary = TRUE, scale = TRUE)

Arguments

variable	The variable that the dynamic smooth will be a function of
k	Optional number of basis functions for computing approximate GPs. If missing, k will be set as large as possible to accurately estimate the nonlinear function.
rho	Either a positive numeric stating the length scale to be used for approximating the squared exponential Gaussian Process smooth (see gp.smooth for details) or missing, in which case the length scale will be estimated by setting up a Hilbert space approximate GP.
stationary	Logical. If TRUE (the default) and rho is supplied, the latent Gaussian Process smooth will not have a linear trend component. If FALSE, a linear trend in the covariate is added to the Gaussian Process smooth. Leave at TRUE if you do not believe the coefficient is evolving with much trend, as the linear component of the basis functions can be hard to penalize to zero. This sometimes causes divergence issues in Stan. See gp.smooth for details. Ignored if rho is missing (in which case a Hilbert space approximate GP is used).
scale	Logical; If TRUE (the default) and rho is missing, predictors are scaled so that the maximum Euclidean distance between two points is 1. This often improves sampling speed and convergence. Scaling also affects the estimated length-scale parameters in that they resemble those of scaled predictors (not of the original predictors) if scale is TRUE.

Details

mvgam currently sets up dynamic coefficients as low-rank squared exponential Gaussian Process smooths via the call s(time, by = variable, bs = "gp", m = c(2, rho, 2)). These smooths, if specified with reasonable values for the length scale parameter, will give more realistic out of sample forecasts than standard splines such as thin plate or cubic. But the user must set the value for rho, as there is currently no support for estimating this value in mgcv. This may not be too big of a problem, as estimating latent length scales is often difficult anyway. The rho parameter should be thought of as a prior on the smoothness of the latent dynamic coefficient function (where higher values of rho lead to smoother functions with more temporal covariance structure). Values of k are set automatically to ensure enough basis functions are used to approximate the expected wiggliness of the underlying dynamic function (k will increase as rho decreases).

Value

a list object for internal usage in 'mvgam'

Author(s)

Nicholas J Clark

Examples

## Not run: 
# Simulate a time-varying coefficient
# (as a Gaussian Process with length scale = 10)
set.seed(1111)
N <- 200

# A function to simulate from a squared exponential Gaussian Process
sim_gp <- function(N, c, alpha, rho) {
  Sigma <- alpha ^ 2 *
    exp(-0.5 * ((outer(1:N, 1:N, "-") / rho) ^ 2)) +
    diag(1e-9, N)
  c + mgcv::rmvn(1, mu = rep(0, N), V = Sigma)
}

beta <- sim_gp(alpha = 0.75, rho = 10, c = 0.5, N = N)
plot(
  beta, type = 'l', lwd = 3, bty = 'l',
  xlab = 'Time', ylab = 'Coefficient', col = 'darkred'
)

# Simulate the predictor as a standard normal
predictor <- rnorm(N, sd = 1)

# Simulate a Gaussian outcome variable
out <- rnorm(N, mean = 4 + beta * predictor, sd = 0.25)
time <- seq_along(predictor)
plot(
  out, type = 'l', lwd = 3, bty = 'l',
  xlab = 'Time', ylab = 'Outcome', col = 'darkred'
)

# Gather into a data.frame and fit a dynamic coefficient model
data <- data.frame(out, predictor, time)

# Split into training and testing
data_train <- data[1:190, ]
data_test <- data[191:200, ]

# Fit a model using the dynamic function
mod <- mvgam(
  out ~
    # mis-specify the length scale slightly as this
    # won't be known in practice
    dynamic(predictor, rho = 8, stationary = TRUE),
  family = gaussian(),
  data = data_train,
  chains = 2,
  silent = 2
)

# Inspect the summary
summary(mod)

# Plot the time-varying coefficient estimates
plot(mod, type = 'smooths')

# Extrapolate the coefficient forward in time
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = 'dashed', lwd = 2)

# Overlay the true simulated time-varying coefficient
lines(beta, lwd = 2.5, col = 'white')
lines(beta, lwd = 2)

## End(Not run)

---

Combine forecasts from mvgam models into evenly weighted ensembles

Description

Generate evenly weighted ensemble forecast distributions from mvgam_forecast objects.

Usage

ensemble(object, ...)

## S3 method for class 'mvgam_forecast'
ensemble(object, ..., ndraws = 5000)

Arguments

object	list object of class mvgam_forecast. See forecast.mvgam()
...	More mvgam_forecast objects.
ndraws	Positive integer specifying the number of draws to use from each forecast distribution for creating the ensemble. If some of the ensemble members have fewer draws than ndraws, their forecast distributions will be resampled with replacement to achieve the correct number of draws

Details

It is widely recognised in the forecasting literature that combining forecasts from different models often results in improved forecast accuracy. The simplest way to create an ensemble is to use evenly weighted combinations of forecasts from the different models. This is straightforward to do in a Bayesian setting with mvgam as the posterior MCMC draws contained in each mvgam_forecast object will already implicitly capture correlations among the temporal posterior predictions.

Value

An object of class mvgam_forecast containing the ensemble predictions. This object can be readily used with the supplied S3 functions plot and score.

Author(s)

Nicholas J Clark

See Also

plot.mvgam_forecast, score.mvgam_forecast

Examples

## Not run: 
# Simulate some series and fit a few competing dynamic models
set.seed(1)
simdat <- sim_mvgam(
  n_series = 1,
  prop_trend = 0.6,
  mu = 1
)

plot_mvgam_series(
  data = simdat$data_train,
  newdata = simdat$data_test
)

m1 <- mvgam(
  y ~ 1,
  trend_formula = ~ time +
    s(season, bs = 'cc', k = 9),
  trend_model = AR(p = 1),
  noncentred = TRUE,
  data = simdat$data_train,
  newdata = simdat$data_test,
  chains = 2,
  silent = 2
)

m2 <- mvgam(
  y ~ time,
  trend_model = RW(),
  noncentred = TRUE,
  data = simdat$data_train,
  newdata = simdat$data_test,
  chains = 2,
  silent = 2
)

# Calculate forecast distributions for each model
fc1 <- forecast(m1)
fc2 <- forecast(m2)

# Generate the ensemble forecast
ensemble_fc <- ensemble(fc1, fc2)

# Plot forecasts
plot(fc1)
plot(fc2)
plot(ensemble_fc)

# Score forecasts
score(fc1)
score(fc2)
score(ensemble_fc)

## End(Not run)

---

Evaluate forecasts from fitted mvgam objects

Description

Evaluate forecasts from fitted mvgam objects

Usage

eval_mvgam(
  object,
  n_samples = 5000,
  eval_timepoint = 3,
  fc_horizon = 3,
  n_cores = 1,
  score = "drps",
  log = FALSE,
  weights
)

roll_eval_mvgam(
  object,
  n_evaluations = 5,
  evaluation_seq,
  n_samples = 5000,
  fc_horizon = 3,
  n_cores = 1,
  score = "drps",
  log = FALSE,
  weights
)

compare_mvgams(
  model1,
  model2,
  n_samples = 1000,
  fc_horizon = 3,
  n_evaluations = 10,
  n_cores = 1,
  score = "drps",
  log = FALSE,
  weights
)

Arguments

object	list object returned from mvgam
n_samples	integer specifying the number of samples to generate from the model's posterior distribution
eval_timepoint	integer indexing the timepoint that represents our last 'observed' set of outcome data
fc_horizon	integer specifying the length of the forecast horizon for evaluating forecasts
n_cores	Deprecated. Parallel processing is no longer supported
score	character specifying the type of ranked probability score to use for evaluation. Options are: variogram, drps or crps
log	logical. Should the forecasts and truths be logged prior to scoring? This is often appropriate for comparing performance of models when series vary in their observation ranges
weights	optional vector of weights (where length(weights) == n_series) for weighting pairwise correlations when evaluating the variogram score for multivariate forecasts. Useful for down-weighting series that have larger magnitude observations or that are of less interest when forecasting. Ignored if score != 'variogram'
n_evaluations	integer specifying the total number of evaluations to perform
evaluation_seq	Optional integer sequence specifying the exact set of timepoints for evaluating the model's forecasts. This sequence cannot have values <3 or > max(training timepoints) - fc_horizon
model1	list object returned from mvgam representing the first model to be evaluated
model2	list object returned from mvgam representing the second model to be evaluated

Details

eval_mvgam may be useful when both repeated fitting of a model using update.mvgam for exact leave-future-out cross-validation and approximate leave-future-out cross-validation using lfo_cv are impractical. The function generates a set of samples representing fixed parameters estimated from the full mvgam model and latent trend states at a given point in time. The trends are rolled forward a total of fc_horizon timesteps according to their estimated state space dynamics to generate an 'out-of-sample' forecast that is evaluated against the true observations in the horizon window. This function therefore simulates a situation where the model's parameters had already been estimated but we have only observed data up to the evaluation timepoint and would like to generate forecasts from the latent trends that have been observed up to that timepoint. Evaluation involves calculating an appropriate Rank Probability Score and a binary indicator for whether or not the true value lies within the forecast's 90% prediction interval

roll_eval_mvgam sets up a sequence of evaluation timepoints along a rolling window and iteratively calls eval_mvgam to evaluate 'out-of-sample' forecasts. Evaluation involves calculating the Rank Probability Scores and a binary indicator for whether or not the true value lies within the forecast's 90% prediction interval

compare_mvgams automates the evaluation to compare two fitted models using rolling window forecast evaluation and provides a series of summary plots to facilitate model selection. It is essentially a wrapper for roll_eval_mvgam

Value

For eval_mvgam, a list object containing information on specific evaluations for each series (if using drps or crps as the score) or a vector of scores when using variogram.

For roll_eval_mvgam, a list object containing information on specific evaluations for each series as well as a total evaluation summary (taken by summing the forecast score for each series at each evaluation and averaging the coverages at each evaluation)

For compare_mvgams, a series of plots comparing forecast Rank Probability Scores for each competing model. A lower score is preferred. Note however that it is possible to select a model that ultimately would perform poorly in true out-of-sample forecasting. For example if a wiggly smooth function of 'year' is included in the model then this function will be learned prior to evaluating rolling window forecasts, and the model could generate very tight predictions as a result. But when forecasting ahead to timepoints that the model has not seen (i.e. next year), the smooth function will end up extrapolating, sometimes in very strange and unexpected ways. It is therefore recommended to only use smooth functions for covariates that are adequately measured in the data (i.e. 'seasonality', for example) to reduce possible extrapolation of smooths and let the latent trends in the mvgam model capture any temporal dependencies in the data. These trends are time series models and so will provide much more stable forecasts

See Also

forecast, score, lfo_cv

Examples

## Not run: 
# Simulate from a Poisson-AR2 model with a seasonal smooth
set.seed(1)
dat <- sim_mvgam(
  T = 75,
  n_series = 1,
  prop_trend = 0.75,
  trend_model = AR(p = 2),
  family = poisson()
)

# Fit an appropriate model
mod_ar2 <- mvgam(
  formula = y ~ s(season, bs = 'cc'),
  trend_model = AR(p = 2),
  family = poisson(),
  data = dat$data_train,
  newdata = dat$data_test,
  chains = 2,
  silent = 2
)

# Fit a less appropriate model
mod_rw <- mvgam(
  formula = y ~ 1,
  trend_model = RW(),
  family = poisson(),
  data = dat$data_train,
  newdata = dat$data_test,
  chains = 2,
  silent = 2
)

# Compare Discrete Ranked Probability Scores for the testing period
fc_ar2 <- forecast(mod_ar2)
fc_rw <- forecast(mod_rw)
score_ar2 <- score(
  object = fc_ar2,
  score = 'drps'
)
score_rw <- score(
  object = fc_rw,
  score = 'drps'
)
sum(score_ar2$series_1$score)
sum(score_rw$series_1$score)

# Use rolling evaluation for approximate comparisons of 3-step ahead
# forecasts across the training period
compare_mvgams(
  model1 = mod_ar2,
  model2 = mod_rw,
  fc_horizon = 3,
  n_samples = 1000,
  n_evaluations = 5
)

# A more appropriate comparison would be to use approximate
# leave-future-out CV to compare forecasts (see ?mvgam::lfo_cv())

## End(Not run)

---

Calculate latent VAR forecast error variance decompositions

Description

Compute forecast error variance decompositions from mvgam models with Vector Autoregressive dynamics

Usage

fevd(object, ...)

## S3 method for class 'mvgam'
fevd(object, h = 10, ...)

Arguments

object	list object of class mvgam resulting from a call to mvgam() that used a Vector Autoregressive latent process model (either as VAR(cor = FALSE) or VAR(cor = TRUE); see VAR() for details)
...	ignored
h	Positive integer specifying the forecast horizon over which to calculate the IRF

Value

See mvgam_fevd-class for a full description of the quantities that are computed and returned by this function, along with key references.

Author(s)

Nicholas J Clark

References

Lütkepohl, H. (2007). New Introduction to Multiple Time Series Analysis. 2nd ed. Springer-Verlag Berlin Heidelberg.

See Also

VAR(), irf(), stability(), mvgam_fevd-class

Examples

## Not run: 
# Simulate some time series that follow a latent VAR(1) process
simdat <- sim_mvgam(
  family = gaussian(),
  n_series = 4,
  trend_model = VAR(cor = TRUE),
  prop_trend = 1
)
plot_mvgam_series(data = simdat$data_train, series = "all")

# Fit a model that uses a latent VAR(1)
mod <- mvgam(
  formula = y ~ -1,
  trend_formula = ~ 1,
  trend_model = VAR(cor = TRUE),
  family = gaussian(),
  data = simdat$data_train,
  chains = 2,
  silent = 2
)

# Plot the autoregressive coefficient distributions;
# use 'dir = "v"' to arrange the order of facets
# correctly
mcmc_plot(
  mod,
  variable = 'A',
  regex = TRUE,
  type = 'hist',
  facet_args = list(dir = 'v')
)

# Calulate forecast error variance decompositions for each series
fevds <- fevd(mod, h = 12)

# Plot median contributions to forecast error variance
plot(fevds)

# View a summary of the error variance decompositions
summary(fevds)

## End(Not run)

---

Expected values of the posterior predictive distribution for mvgam objects

Description

This method extracts posterior estimates of the fitted values (i.e. the actual predictions, including estimates for any trend states, that were obtained when fitting the model). It also includes an option for obtaining summaries of the computed draws.

Usage

## S3 method for class 'mvgam'
fitted(
  object,
  process_error = TRUE,
  scale = c("response", "linear"),
  summary = TRUE,
  robust = FALSE,
  probs = c(0.025, 0.975),
  ...
)

Arguments

object	An object of class mvgam
process_error	Logical. If TRUE and a dynamic trend model was fit, expected uncertainty in the process model is accounted for by using draws from a stationary, zero-centred multivariate Normal distribution using any estimated process variance-covariance parameters. If FALSE, uncertainty in the latent trend component is ignored when calculating predictions
scale	Either "response" or "linear". If "response", results are returned on the scale of the response variable. If "linear", results are returned on the scale of the linear predictor term, that is without applying the inverse link function or other transformations.
summary	Should summary statistics be returned instead of the raw values? Default is TRUE..
robust	If FALSE (the default) the mean is used as the measure of central tendency and the standard deviation as the measure of variability. If TRUE, the median and the median absolute deviation (MAD) are applied instead. Only used if summary is TRUE.
probs	The percentiles to be computed by the quantile function. Only used if summary is TRUE.
...	Further arguments passed to prepare_predictions that control several aspects of data validation and prediction.

Details

This method gives the actual fitted values from the model (i.e. what you will see if you generate hindcasts from the fitted model using hindcast.mvgam with type = 'expected'). These predictions can be overly precise if a flexible dynamic trend component was included in the model. This is in contrast to the set of predict functions (i.e. posterior_epred.mvgam or predict.mvgam), which will assume any dynamic trend component has reached stationarity when returning hypothetical predictions.

Value

An array of predicted mean response values.

If summary = FALSE the output resembles those of posterior_epred.mvgam and predict.mvgam.

If summary = TRUE the output is an n_observations x E matrix. The number of summary statistics E is equal to 2 + length(probs): The Estimate column contains point estimates (either mean or median depending on argument robust), while the Est.Error column contains uncertainty estimates (either standard deviation or median absolute deviation depending on argument robust). The remaining columns starting with Q contain quantile estimates as specified via argument probs.

Author(s)

Nicholas J Clark

See Also

hindcast.mvgam

Examples

## Not run: 
# Simulate some data and fit a model
simdat <- sim_mvgam(n_series = 1, trend_model = AR())

mod <- mvgam(
  y ~ s(season, bs = 'cc'),
  trend_model = AR(),
  data = simdat$data_train,
  chains = 2,
  silent = 2
)

# Extract fitted values (posterior expectations)
expectations <- fitted(mod)
str(expectations)

## End(Not run)

---

Extract or compute hindcasts and forecasts for a fitted mvgam object

Description

Extract or compute hindcasts and forecasts for a fitted mvgam object

Usage

## S3 method for class 'mvgam'
forecast(object, newdata, data_test, n_cores = 1, type = "response", ...)

Arguments

object	list object of class mvgam or jsdgam. See mvgam()
newdata	Optional dataframe or list of test data containing the same variables that were included in the original data used to fit the model. If included, the covariate information in newdata will be used to generate forecasts from the fitted model equations. If this same newdata was originally included in the call to mvgam, then forecasts have already been produced by the generative model and these will simply be extracted and plotted. However if no newdata was supplied to the original model call, an assumption is made that the newdata supplied here comes sequentially after the data supplied in the original model (i.e. we assume there is no time gap between the last observation of series 1 in the original data and the first observation for series 1 in newdata)
data_test	Deprecated. Still works in place of newdata but users are recommended to use newdata instead for more seamless integration into R workflows
n_cores	Deprecated. Parallel processing is no longer supported
type	When this has the value link (default) the linear predictor is calculated on the link scale. If expected is used, predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. When response is used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. When variance is used, the variance of the response with respect to the mean (mean-variance relationship) is returned. When type = "terms", each component of the linear predictor is returned separately in the form of a list (possibly with standard errors, if summary = TRUE): this includes parametric model components, followed by each smooth component, but excludes any offset and any intercept. Two special cases are also allowed: type latent_N will return the estimated latent abundances from an N-mixture distribution, while type detection will return the estimated detection probability from an N-mixture distribution
...	Ignored

Details

Posterior predictions are drawn from the fitted mvgam and used to simulate a forecast distribution

Value

An object of class mvgam_forecast containing hindcast and forecast distributions. See mvgam_forecast-class for details.

See Also

hindcast.mvgam(), plot.mvgam_forecast(), summary.mvgam_forecast(), score.mvgam_forecast() ensemble.mvgam_forecast()

Examples

## Not run: 
  # Simulate data with 3 series and AR trend model
  simdat <- sim_mvgam(n_series = 3, trend_model = AR())

  # Fit mvgam model
  mod <- mvgam(
    y ~ s(season, bs = 'cc', k = 6),
    trend_model = AR(),
    noncentred = TRUE,
    data = simdat$data_train,
    chains = 2,
    silent = 2
  )

  # Hindcasts on response scale
  hc <- hindcast(mod)
  str(hc)

  # Use summary() to extract hindcasts / forecasts for custom plotting
  head(summary(hc), 12)

  # Or just use the plot() function for quick plots
  plot(hc, series = 1)
  plot(hc, series = 2)
  plot(hc, series = 3)

  # Forecasts on response scale
  fc <- forecast(
    mod,
    newdata = simdat$data_test
  )
  str(fc)
  head(summary(fc), 12)
  plot(fc, series = 1)
  plot(fc, series = 2)
  plot(fc, series = 3)

  # Forecasts as expectations
  fc <- forecast(
    mod,
    newdata = simdat$data_test,
    type = 'expected'
  )
  head(summary(fc), 12)
  plot(fc, series = 1)
  plot(fc, series = 2)
  plot(fc, series = 3)

  # Dynamic trend extrapolations
  fc <- forecast(
    mod,
    newdata = simdat$data_test,
    type = 'trend'
  )
  head(summary(fc), 12)
  plot(fc, series = 1)
  plot(fc, series = 2)
  plot(fc, series = 3)

## End(Not run)

---

Extract formulae from mvgam objects

Description

Extract formulae from mvgam objects

Usage

## S3 method for class 'mvgam'
formula(x, trend_effects = FALSE, ...)

## S3 method for class 'mvgam_prefit'
formula(x, trend_effects = FALSE, ...)

Arguments

x	mvgam, jsdgam or mvgam_prefit object
trend_effects	logical, return the formula from the observation model (if FALSE) or from the underlying process model (ifTRUE)
...	Ignored

Value

A formula object

Author(s)

Nicholas J Clark

---

Extract information on default prior distributions for an mvgam model

Description

This function lists the parameters that can have their prior distributions changed for a given model, as well listing their default distributions

Usage

get_mvgam_priors(
  formula,
  trend_formula,
  factor_formula,
  knots,
  trend_knots,
  trend_model = "None",
  family = poisson(),
  data,
  unit = time,
  species = series,
  use_lv = FALSE,
  n_lv,
  trend_map,
  ...
)

Arguments

formula	A formula object specifying the GAM observation model formula. These are exactly like the formula for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, nonparametric gp() terms and offsets using offset(), can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these). In nmix() family models, the formula is used to set up a linear predictor for the detection probability. Details of the formula syntax used by mvgam can be found in mvgam_formulae
trend_formula	An optional formula object specifying the GAM process model formula. If supplied, a linear predictor will be modelled for the latent trends to capture process model evolution separately from the observation model. Important notes: Should not have a response variable specified on the left-hand side (e.g., ~ season + s(year)) Use trend instead of series for effects that vary across time series Only available for RW(), AR() and VAR() trend models In nmix() family models, sets up linear predictor for latent abundance Consider dropping one intercept using - 1 convention to avoid estimation challenges
factor_formula	Can be supplied instead trend_formula to match syntax from jsdgam
knots	An optional list containing user specified knot values for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied. Different terms can use different numbers of knots, unless they share a covariate.
trend_knots	As for knots above, this is an optional list of knot values for smooth functions within the trend_formula.
trend_model	character or function specifying the time series dynamics for the latent trend. Available options: None: No latent trend component (GAM component only, like gam) ZMVN or ZMVN(): Zero-Mean Multivariate Normal (Stan only) 'RW' or RW(): Random Walk 'AR1', 'AR2', 'AR3' or AR(p = 1, 2, 3): Autoregressive models 'CAR1' or CAR(p = 1): Continuous-time AR (Ornstein–Uhlenbeck process) 'VAR1' or VAR(): Vector Autoregressive (Stan only) 'PWlogistic', 'PWlinear' or PW(): Piecewise trends (Stan only) 'GP' or GP(): Gaussian Process with squared exponential kernel (Stan only) Additional features: Moving average and/or correlated process error terms available for most types (e.g., RW(cor = TRUE) for multivariate Random Walk) Hierarchical correlations possible for structured data See mvgam_trends for details and ZMVN() for examples
family	family specifying the exponential observation family for the series. Supported families: gaussian(): Real-valued data betar(): Proportional data on ⁠(0,1)⁠ lognormal(): Non-negative real-valued data student_t(): Real-valued data Gamma(): Non-negative real-valued data bernoulli(): Binary data poisson(): Count data (default) nb(): Overdispersed count data binomial(): Count data with imperfect detection when number of trials is known (use cbind() to bind observations and trials) beta_binomial(): As binomial() but allows for overdispersion nmix(): Count data with imperfect detection when number of trials is unknown (State-Space N-Mixture model with Poisson latent states and Binomial observations) See mvgam_families for more details.
data	A dataframe or list containing the model response variable and covariates required by the GAM formula and optional trend_formula. Required columns for most models: series: A factor index of the series IDs (number of levels should equal number of unique series labels) time: numeric or integer index of time points. For most dynamic trend types, time should be measured in discrete, regularly spaced intervals (i.e., c(1, 2, 3, ...)). Irregular spacing is allowed for trend_model = CAR(1), but zero intervals are adjusted to 1e-12 to prevent sampling errors. Special cases: Models with hierarchical temporal correlation (e.g., AR(gr = region, subgr = species)) should NOT include a series identifier Models without temporal dynamics (trend_model = 'None' or trend_model = ZMVN()) don't require a time variable
unit	The unquoted name of the variable that represents the unit of analysis in data over which latent residuals should be correlated. This variable should be either a numeric or integer variable in the supplied data. Defaults to time to be consistent with other functionalities in mvgam, though note that the data need not be time series in this case. See examples below for further details and explanations
species	The unquoted name of the factor variable that indexes the different response units in data (usually 'species' in a JSDM). Defaults to series to be consistent with other mvgam models
use_lv	logical. If TRUE, use dynamic factors to estimate series' latent trends in a reduced dimension format. Only available for RW(), AR() and GP() trend models. Default is FALSE. See lv_correlations for examples.
n_lv	integer specifying the number of latent dynamic factors to use if use_lv == TRUE. Cannot exceed n_series. Default is min(2, floor(n_series / 2)).
trend_map	Optional data.frame specifying which series should depend on which latent trends. Enables multiple series to depend on the same latent trend process with different observation processes. Required structure: Column series: Single unique entry for each series (matching factor levels in data) Column trend: Integer values indicating which trend each series depends on Notes: Sets up latent factor model by enabling use_lv = TRUE Process model intercept is NOT automatically suppressed Not yet supported for continuous time models (CAR())
...	Not currently used

Details

Users can supply a model formula, prior to fitting the model, so that default priors can be inspected and altered. To make alterations, change the contents of the prior column and supplying this data.frame to the mvgam or jsdgam functions using the argument priors. If using Stan as the backend, users can also modify the parameter bounds by modifying the new_lowerbound and/or new_upperbound columns. This will be necessary if using restrictive distributions on some parameters, such as a Beta distribution for the trend sd parameters for example (Beta only has support on (0,1)), so the upperbound cannot be above 1. Another option is to make use of the prior modification functions in brms (i.e. prior) to change prior distributions and bounds (just use the name of the parameter that you'd like to change as the class argument; see examples below)

Value

either a data.frame containing the prior definitions (if any suitable priors can be altered by the user) or NULL, indicating that no priors in the model can be modified

Note

Only the prior, new_lowerbound and/or new_upperbound columns of the output should be altered when defining the user-defined priors for the model. Use only if you are familiar with the underlying probabilistic programming language. There are no sanity checks done to ensure that the code is legal (i.e. to check that lower bounds are smaller than upper bounds, for example)

Author(s)

Nicholas J Clark

See Also

mvgam, mvgam_formulae, prior

Examples

## Not run: 
# ========================================================================
# Example 1: Simulate data and inspect default priors
# ========================================================================

dat <- sim_mvgam(trend_rel = 0.5)

# Get a model file that uses default mvgam priors for inspection (not
# always necessary, but this can be useful for testing whether your
# updated priors are written correctly)
mod_default <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  run_model = FALSE
)

# Inspect the model file with default mvgam priors
stancode(mod_default)

# Look at which priors can be updated in mvgam
test_priors <- get_mvgam_priors(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2)
)
test_priors

# ========================================================================
# Example 2: Modify priors manually
# ========================================================================

# Make a few changes; first, change the population mean for the
# series-level random intercepts
test_priors$prior[2] <- "mu_raw ~ normal(0.2, 0.5);"

# Now use stronger regularisation for the series-level AR2 coefficients
test_priors$prior[5] <- "ar2 ~ normal(0, 0.25);"

# Check that the changes are made to the model file without any warnings
# by setting 'run_model = FALSE'
mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = test_priors,
  run_model = FALSE
)
stancode(mod)

# No warnings, the model is ready for fitting now in the usual way with
# the addition of the 'priors' argument

# ========================================================================
# Example 3: Use brms syntax for prior modification
# ========================================================================

# The same can be done using 'brms' functions; here we will also change
# the ar1 prior and put some bounds on the ar coefficients to enforce
# stationarity; we set the prior using the 'class' argument in all brms
# prior functions
brmsprior <- c(
  prior(normal(0.2, 0.5), class = mu_raw),
  prior(normal(0, 0.25), class = ar1, lb = -1, ub = 1),
  prior(normal(0, 0.25), class = ar2, lb = -1, ub = 1)
)
brmsprior

mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = brmsprior,
  run_model = FALSE
)
stancode(mod)

# ========================================================================
# Example 4: Error handling example
# ========================================================================

# Look at what is returned when an incorrect spelling is used
test_priors$prior[5] <- "ar2_bananas ~ normal(0, 0.25);"
mod <- mvgam(
  y ~ s(series, bs = "re") + s(season, bs = "cc") - 1,
  family = nb(),
  data = dat$data_train,
  trend_model = AR(p = 2),
  priors = test_priors,
  run_model = FALSE
)
stancode(mod)

# ========================================================================
# Example 5: Parametric (fixed effect) priors
# ========================================================================

simdat <- sim_mvgam()

# Add a fake covariate
simdat$data_train$cov <- rnorm(NROW(simdat$data_train))

priors <- get_mvgam_priors(
  y ~ cov + s(season),
  data = simdat$data_train,
  family = poisson(),
  trend_model = AR()
)

# Change priors for the intercept and fake covariate effects
priors$prior[1] <- "(Intercept) ~ normal(0, 1);"
priors$prior[2] <- "cov ~ normal(0, 0.1);"

mod2 <- mvgam(
  y ~ cov + s(season),
  data = simdat$data_train,
  trend_model = AR(),
  family = poisson(),
  priors = priors,
  run_model = FALSE
)
stancode(mod2)

# ========================================================================
# Example 6: Alternative brms syntax for fixed effects
# ========================================================================

# Likewise using 'brms' utilities (note that you can use Intercept rather
# than `(Intercept)`) to change priors on the intercept
brmsprior <- c(
  prior(normal(0.2, 0.5), class = cov),
  prior(normal(0, 0.25), class = Intercept)
)
brmsprior

mod2 <- mvgam(
  y ~ cov + s(season),
  data = simdat$data_train,
  trend_model = AR(),
  family = poisson(),
  priors = brmsprior,
  run_model = FALSE
)
stancode(mod2)

# ========================================================================
# Example 7: Bulk prior assignment
# ========================================================================

# The "class = 'b'" shortcut can be used to put the same prior on all
# 'fixed' effect coefficients (apart from any intercepts)
set.seed(0)
dat <- mgcv::gamSim(1, n = 200, scale = 2)
dat$time <- 1:NROW(dat)
mod <- mvgam(
  y ~ x0 + x1 + s(x2) + s(x3),
  priors = prior(normal(0, 0.75), class = "b"),
  data = dat,
  family = gaussian(),
  run_model = FALSE
)
stancode(mod)

## End(Not run)

---

Specify dynamic Gaussian process trends in mvgam models

Description

Set up low-rank approximate Gaussian Process trend models using Hilbert basis expansions in mvgam. This function does not evaluate its arguments – it exists purely to help set up a model with particular GP trend models.

Usage

GP(...)

Arguments

...

unused

Details

A GP trend is estimated for each series using Hilbert space approximate Gaussian Processes. In mvgam, latent squared exponential GP trends are approximated using by default 20 basis functions and using a multiplicative factor of c = 5/4, which saves computational costs compared to fitting full GPs while adequately estimating GP alpha and rho parameters.

Value

An object of class mvgam_trend, which contains a list of arguments to be interpreted by the parsing functions in mvgam.

Author(s)

Nicholas J Clark

References

Riutort-Mayol G, Burkner PC, Andersen MR, Solin A and Vehtari A (2023). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing 33, 1. https://doi.org/10.1007/s11222-022-10167-2

See Also

gp

---

Enhance post-processing of mvgam models using gratia functionality

Description

These evaluation and plotting functions exist to allow some popular gratia methods to work with mvgam or jsdgam models

Usage

drawDotmvgam(
  object,
  trend_effects = FALSE,
  data = NULL,
  select = NULL,
  parametric = FALSE,
  terms = NULL,
  residuals = FALSE,
  scales = c("free", "fixed"),
  ci_level = 0.95,
  n = 100,
  n_3d = 16,
  n_4d = 4,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  constant = NULL,
  fun = NULL,
  dist = 0.1,
  rug = TRUE,
  contour = TRUE,
  grouped_by = FALSE,
  ci_alpha = 0.2,
  ci_col = "black",
  smooth_col = "black",
  resid_col = "steelblue3",
  contour_col = "black",
  n_contour = NULL,
  partial_match = FALSE,
  discrete_colour = NULL,
  discrete_fill = NULL,
  continuous_colour = NULL,
  continuous_fill = NULL,
  position = "identity",
  angle = NULL,
  ncol = NULL,
  nrow = NULL,
  guides = "keep",
  widths = NULL,
  heights = NULL,
  crs = NULL,
  default_crs = NULL,
  lims_method = "cross",
  wrap = TRUE,
  envir = environment(formula(object)),
  ...
)

eval_smoothDothilbertDotsmooth(
  smooth,
  model,
  n = 100,
  n_3d = NULL,
  n_4d = NULL,
  data = NULL,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  dist = NULL,
  ...
)

eval_smoothDotmodDotsmooth(
  smooth,
  model,
  n = 100,
  n_3d = NULL,
  n_4d = NULL,
  data = NULL,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  dist = NULL,
  ...
)

eval_smoothDotmoiDotsmooth(
  smooth,
  model,
  n = 100,
  n_3d = NULL,
  n_4d = NULL,
  data = NULL,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  dist = NULL,
  ...
)

Arguments

object	a fitted mvgam, the result of a call to mvgam()
trend_effects	logical specifying whether smooth terms from the trend_formula should be drawn. If FALSE, only terms from the observation formula are drawn. If TRUE, only terms from the trend_formula are drawn
data	a data frame of covariate values at which to evaluate the model's smooth functions
select	character, logical, or numeric; which smooths to plot. If NULL, the default, then all model smooths are drawn. Character select matches the labels for smooths as shown for example in the output from summary(object). Logical select operates as per numeric select in the order that smooths are stored
parametric	logical; plot parametric terms also? Note that select is used for selecting which smooths to plot. The terms argument is used to select which parametric effects are plotted. The default, as with mgcv::plot.gam(), is to not draw parametric effects
terms	character; which model parametric terms should be drawn? The Default of NULL will plot all parametric terms that can be drawn.
residuals	currently ignored for mvgam models
scales	character; should all univariate smooths be plotted with the same y-axis scale? If scales = "free", the default, each univariate smooth has its own y-axis scale. If scales = "fixed", a common y axis scale is used for all univariate smooths. Currently does not affect the y-axis scale of plots of the parametric terms
ci_level	numeric between 0 and 1; the coverage of credible interval.
n	numeric; the number of points over the range of the covariate at which to evaluate the smooth
n_3d, n_4d	numeric; the number of points over the range of last covariate in a 3D or 4D smooth. The default is NULL which achieves the standard behaviour of using n points over the range of all covariate, resulting in n^d evaluation points, where d is the dimension of the smooth. For d > 2 this can result in very many evaluation points and slow performance. For smooths of d > 4, the value of n_4d will be used for all dimensions ⁠> 4⁠, unless this is NULL, in which case the default behaviour (using n for all dimensions) will be observed
unconditional	ignored for mvgam models as all appropriate uncertainties are already included in the posterior estimates
overall_uncertainty	ignored for mvgam models as all appropriate uncertainties are already included in the posterior estimates
constant	numeric; a constant to add to the estimated values of the smooth. constant, if supplied, will be added to the estimated value before the confidence band is computed
fun	function; a function that will be applied to the estimated values and confidence interval before plotting. Can be a function or the name of a function. Function fun will be applied after adding any constant, if provided
dist	numeric; if greater than 0, this is used to determine when a location is too far from data to be plotted when plotting 2-D smooths. The data are scaled into the unit square before deciding what to exclude, and dist is a distance within the unit square. See mgcv::exclude.too.far() for further details
rug	logical; draw a rug plot at the bottom of each plot for 1-D smooths or plot locations of data for higher dimensions.
contour	logical; should contours be draw on the plot using ggplot2::geom_contour()
grouped_by	logical; should factor by smooths be drawn as one panel per level of the factor (FALSE, the default), or should the individual smooths be combined into a single panel containing all levels (TRUE)?
ci_alpha	numeric; alpha transparency for confidence or simultaneous interval
ci_col	colour specification for the confidence/credible intervals band. Affects the fill of the interval
smooth_col	colour specification for the smooth line
resid_col	colour specification for residual points. Ignored
contour_col	colour specification for contour lines
n_contour	numeric; the number of contour bins. Will result in n_contour - 1 contour lines being drawn. See ggplot2::geom_contour()
partial_match	logical; should smooths be selected by partial matches with select? If TRUE, select can only be a single string to match against
discrete_colour	a suitable colour scale to be used when plotting discrete variables
discrete_fill	a suitable fill scale to be used when plotting discrete variables.
continuous_colour	a suitable colour scale to be used when plotting continuous variables
continuous_fill	a suitable fill scale to be used when plotting continuous variables
position	Position adjustment, either as a string, or the result of a call to a position adjustment function
angle	numeric; the angle at which the x axis tick labels are to be drawn passed to the angle argument of ggplot2::guide_axis()
ncol, nrow	numeric; the numbers of rows and columns over which to spread the plots
guides	character; one of "keep" (the default), "collect", or "auto". Passed to patchwork::plot_layout()
widths, heights	The relative widths and heights of each column and row in the grid. Will get repeated to match the dimensions of the grid. If there is more than 1 plot and widths = NULL, the value of widths will be set internally to widths = 1 to accommodate plots of smooths that use a fixed aspect ratio.=
crs	the coordinate reference system (CRS) to use for the plot. All data will be projected into this CRS. See ggplot2::coord_sf() for details
default_crs	the coordinate reference system (CRS) to use for the non-sf layers in the plot. If left at the default NULL, the CRS used is 4326 (WGS84), which is appropriate for spline-on-the-sphere smooths, which are parameterized in terms of latitude and longitude as coordinates. See ggplot2::coord_sf() for more details
lims_method	character; affects how the axis limits are determined. See ggplot2::coord_sf(). Be careful; in testing of some examples, changing this to "orthogonal" for example with the chlorophyll-a example from Simon Wood's GAM book quickly used up all the RAM in my test system and the OS killed R. This could be incorrect usage on my part; right now the grid of points at which SOS smooths are evaluated (if not supplied by the user) can produce invalid coordinates for the corners of tiles as the grid is generated for tile centres without respect to the spacing of those tiles
wrap	logical; wrap plots as a patchwork? If FALSE, a list of ggplot objects is returned, 1 per term plotted
envir	an environment to look up the data within
...	additional arguments passed to other methods
smooth	a smooth object of class "gp.smooth" (returned from a model using either the dynamic() function or the gp() function) or of class "moi.smooth" or "mod.smooth" (returned from a model using the 'moi' or 'mod' basis)
model	a fitted mgcv model of clas gam or bam

Details

These methods allow mvgam models to be Enhanced if users have the gratia package installed, making available the popular draw() function to plot partial effects of mvgam smooth functions using ggplot2::ggplot() utilities

Author(s)

Nicholas J Clark

Examples

## Not run: 
# Fit a simple GAM and draw partial effects of smooths using 'gratia'
set.seed(0)
dat <- mgcv::gamSim(
  eg = 1,
  n = 200,
  scale = 2
)

mod <- mvgam(
  formula = y ~ s(x1, bs = 'moi') +
    te(x0, x2),
  data = dat,
  family = gaussian(),
  chains = 2,
  silent = 2
)

if (require("gratia")) {
  gratia::draw(mod)
}

## End(Not run)

---

Extract hindcasts for a fitted mvgam object

Description

Extract hindcasts for a fitted mvgam object

Usage

hindcast(object, ...)

## S3 method for class 'mvgam'
hindcast(object, type = "response", ...)

Arguments

object	list object of class mvgam or jsdgam. See mvgam()
...	Ignored
type	When this has the value link (default) the linear predictor is calculated on the link scale. If expected is used, predictions reflect the expectation of the response (the mean) but ignore uncertainty in the observation process. When response is used, the predictions take uncertainty in the observation process into account to return predictions on the outcome scale. When variance is used, the variance of the response with respect to the mean (mean-variance relationship) is returned. When type = "terms", each component of the linear predictor is returned separately in the form of a list (possibly with standard errors, if summary = TRUE): this includes parametric model components, followed by each smooth component, but excludes any offset and any intercept. Two special cases are also allowed: type latent_N will return the estimated latent abundances from an N-mixture distribution, while type detection will return the estimated detection probability from an N-mixture distribution

Details

Posterior hindcasts (i.e. retrodictions) are drawn from the fitted mvgam and organized into a convenient format for plotting

Value

An object of class mvgam_forecast containing hindcast distributions. See mvgam_forecast-class for details.

See Also

plot.mvgam_forecast(), summary.mvgam_forecast(), forecast.mvgam(), fitted.mvgam(), predict.mvgam()

Examples

## Not run: 
simdat <- sim_mvgam(n_series = 3, trend_model = AR())
mod <- mvgam(y ~ s(season, bs = 'cc'),
             trend_model = AR(),
             noncentred = TRUE,
             data = simdat$data_train,
             chains = 2,
             silent = 2)

# Hindcasts on response scale
hc <- hindcast(mod)
str(hc)
head(summary(hc), 12)
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)

# Hindcasts as expectations
hc <- hindcast(mod, type = 'expected')
head(summary(hc), 12)
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)

# Estimated latent trends
hc <- hindcast(mod, type = 'trend')
head(summary(hc), 12)
plot(hc, series = 1)
plot(hc, series = 2)
plot(hc, series = 3)

## End(Not run)

---

Generate a methods description for mvgam models

Description

Create a brief but fully referenced methods description, along with a useful list of references, for fitted mvgam and jsdgam models.

Usage

how_to_cite(object, ...)

## S3 method for class 'mvgam'
how_to_cite(object, ...)

Arguments

object	list object of class mvgam resulting from a call to mvgam() or jsdgam()
...	ignored

Details

This function uses the model's structure to come up with a very basic but hopefully useful methods description that can help users to appropriately acknowledge the hard work of developers and champion open science. Please do not consider the text returned by this function to be a completely adequate methods section; it is only meant to get you started.

Value

An object of class how_to_cite containing a text description of the methods as well as lists of both primary and additional references.

Author(s)

Nicholas J Clark

See Also

citation, mvgam, jsdgam

Examples

## Not run: 
#--------------------------------------------------
# Simulate 4 time series with hierarchical seasonality
# and a VAR(1) dynamic process
#--------------------------------------------------
set.seed(0)

simdat <- sim_mvgam(
  seasonality = 'hierarchical',
  trend_model = VAR(cor = TRUE),
  family = gaussian()
)

# Fit an appropriate model
mod1 <- mvgam(
  y ~ s(season, bs = 'cc', k = 6),
  data = simdat$data_train,
  family = gaussian(),
  trend_model = VAR(cor = TRUE),
  chains = 2,
  silent = 2
)

how_to_cite(mod1)

#--------------------------------------------------
# For a GP example, simulate data using the mgcv package
#--------------------------------------------------
dat <- mgcv::gamSim(1, n = 30, scale = 2)

# Fit a model that uses an approximate GP from brms
mod2 <- mvgam(
  y ~ gp(x2, k = 12),
  data = dat,
  family = gaussian(),
  chains = 2,
  silent = 2
)

how_to_cite(mod2)

## End(Not run)

---

Index mvgam objects

Description

Index mvgam objects

Usage

## S3 method for class 'mvgam'
variables(x, ...)

Arguments

x	list object returned from mvgam. See mvgam()
...	Arguments passed to individual methods (if applicable).

Value

a list object of the variables that can be extracted, along with their aliases

Author(s)

Nicholas J Clark

Examples

## Not run: 
# Simulate data and fit a model
simdat <- sim_mvgam(
  n_series = 1,
  trend_model = AR()
)

mod <- mvgam(
  y ~ s(season, bs = 'cc', k = 6),
  trend_model = AR(),
  data = simdat$data_train,
  chains = 2,
  silent = 2
)

# Extract model variables
variables(mod)

## End(Not run)

---

Calculate latent VAR impulse response functions

Description

Compute Generalized or Orthogonalized Impulse Response Functions (IRFs) from mvgam models with Vector Autoregressive dynamics

Usage

irf(object, ...)

## S3 method for class 'mvgam'
irf(object, h = 10, cumulative = FALSE, orthogonal = FALSE, ...)

Arguments

object	list object of class mvgam resulting from a call to mvgam() that used a Vector Autoregressive latent process model (either as VAR(cor = FALSE) or VAR(cor = TRUE); see VAR() for details)
...	ignored
h	Positive integer specifying the forecast horizon over which to calculate the IRF
cumulative	Logical flag indicating whether the IRF should be cumulative
orthogonal	Logical flag indicating whether orthogonalized IRFs should be calculated. Note that the order of the variables matters when calculating these

Details

See mvgam_irf-class for a full description of the quantities that are computed and returned by this function, along with key references.

Value

An object of mvgam_irf-class containing the posterior IRFs. This object can be used with the supplied S3 functions plot.mvgam_irf() and summary.mvgam_irf()

Author(s)

Nicholas J Clark

See Also

mvgam_irf-class, VAR(), plot.mvgam_irf(), stability(), fevd()

Examples

## Not run: 
# Fit a model to the portal time series that uses a latent VAR(1)
mod <- mvgam(
  formula = captures ~ -1,
  trend_formula = ~ trend,
  trend_model = VAR(cor = TRUE),
  family = poisson(),
  data = portal_data,
  chains = 2,
  silent = 2
)

# Plot the autoregressive coefficient distributions;
# use 'dir = "v"' to arrange the order of facets
# correctly
mcmc_plot(
  mod,
  variable = 'A',
  regex = TRUE,
  type = 'hist',
  facet_args = list(dir = 'v')
)

# Calulate Generalized IRFs for each series
irfs <- irf(
  mod,
  h = 12,
  cumulative = FALSE
)

# Plot them
plot(irfs, series = 1)
plot(irfs, series = 2)
plot(irfs, series = 3)
plot(irfs, series = 4)

# Calculate posterior median, upper and lower 95th quantiles
# of the impulse responses
summary(irfs)

## End(Not run)

---

Fit Joint Species Distribution Models in mvgam

Description

This function sets up a Joint Species Distribution Model whereby the residual associations among species can be modelled in a reduced-rank format using a set of latent factors. The factor specification is extremely flexible, allowing users to include spatial, temporal or any other type of predictor effects to more efficiently capture unmodelled residual associations, while the observation model can also be highly flexible (including all smooth, GP and other effects that mvgam can handle)

Usage

jsdgam(
  formula,
  factor_formula = ~-1,
  knots,
  factor_knots,
  data,
  newdata,
  family = poisson(),
  unit = time,
  species = series,
  share_obs_params = FALSE,
  priors,
  n_lv = 2,
  backend = getOption("brms.backend", "cmdstanr"),
  algorithm = getOption("brms.algorithm", "sampling"),
  control = list(max_treedepth = 10, adapt_delta = 0.8),
  chains = 4,
  burnin = 500,
  samples = 500,
  thin = 1,
  parallel = TRUE,
  threads = 1,
  silent = 1,
  run_model = TRUE,
  return_model_data = FALSE,
  residuals = TRUE,
  ...
)

Arguments

formula	A formula object specifying the GAM observation model formula. These are exactly like the formula for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, nonparametric gp() terms and offsets using offset(), can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these). Details of the formula syntax used by mvgam can be found in mvgam_formulae
factor_formula	A formula object specifying the linear predictor effects for the latent factors. Use by = trend within calls to functional terms (i.e. s(), te(), ti(), t2(), dynamic(), or gp()) to ensure that each factor captures a different axis of variation. See the example below as an illustration
knots	An optional list containing user specified knot values for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied. Different terms can use different numbers of knots, unless they share a covariate.
factor_knots	An optional list containing user specified knot values to be used for basis construction of any smooth terms in factor_formula. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied (note that the number of knots is not always just k). Different terms can use different numbers of knots, unless they share a covariate
data	A dataframe or list containing the model response variable and covariates required by the GAM formula and factor_formula objects
newdata	Optional dataframe or list of test data containing the same variables as in data. If included, observations in variable y will be set to NA when fitting the model so that posterior simulations can be obtained.
family	family specifying the observation family for the outcomes. Currently supported families are: gaussian() for real-valued data betar() for proportional data on ⁠(0,1)⁠ lognormal() for non-negative real-valued data student_t() for real-valued data Gamma() for non-negative real-valued data bernoulli() for binary data poisson() for count data nb() for overdispersed count data binomial() for count data with imperfect detection when the number of trials is known; note that the cbind() function must be used to bind the discrete observations and the discrete number of trials beta_binomial() as for binomial() but allows for overdispersion Default is poisson(). See mvgam_families for more details
unit	The unquoted name of the variable that represents the unit of analysis in data over which latent residuals should be correlated. This variable should be either a numeric or integer variable in the supplied data. Defaults to time to be consistent with other functionalities in mvgam, though note that the data need not be time series in this case. See examples below for further details and explanations
species	The unquoted name of the factor variable that indexes the different response units in data (usually 'species' in a JSDM). Defaults to series to be consistent with other mvgam models
share_obs_params	logical. If TRUE and the family has additional family-specific observation parameters (e.g., variance components, dispersion parameters), these will be shared across all outcome variables. Useful when multiple outcomes share properties. Default is FALSE.
priors	An optional data.frame with prior definitions (in Stan syntax) or, preferentially, a vector containing objects of class brmsprior (see. prior for details). See get_mvgam_priors and for more information on changing default prior distributions
n_lv	integer the number of latent factors to use for modelling residual associations. Cannot be ⁠> n_species⁠. Defaults arbitrarily to 2
backend	Character string naming the package for Stan model fitting. Options are "cmdstanr" (default) or "rstan". Can be set globally via "brms.backend" option. See https://mc-stan.org/rstan/ and https://mc-stan.org/cmdstanr/ for details.
algorithm	Character string naming the estimation approach: "sampling": MCMC (default) "meanfield": Variational inference with factorized normal distributions "fullrank": Variational inference with multivariate normal distribution "laplace": Laplace approximation (cmdstanr only) "pathfinder": Pathfinder algorithm (cmdstanr only) Can be set globally via "brms.algorithm" option. Limited testing suggests "meanfield" performs best among non-MCMC approximations for dynamic GAMs.
control	Named list for controlling sampler behaviour. Valid elements include max_treedepth, adapt_delta and init.
chains	integer specifying the number of parallel chains for the model. Ignored for variational inference algorithms.
burnin	integer specifying the number of warmup iterations to tune sampling algorithms. Ignored for variational inference algorithms.
samples	integer specifying the number of post-warmup iterations for sampling the posterior distribution.
thin	Thinning interval for monitors. Ignored for variational inference algorithms.
parallel	logical specifying whether to use multiple cores for parallel MCMC simulation. If TRUE, uses min(c(chains, parallel::detectCores() - 1)) cores.
threads	integer Experimental option to use multithreading for within-chain parallelisation in Stan. We recommend its use only if you are experienced with Stan's reduce_sum function and have a slow running model that cannot be sped up by any other means. Currently works for all families when using cmdstanr as the backend
silent	Verbosity level between 0 and 2. If 1 (default), most informational messages are suppressed. If 2, even more messages are suppressed. Sampling progress is still printed - set refresh = 0 to disable. For backend = "rstan", also set open_progress = FALSE to prevent additional progress bars.
run_model	logical. If FALSE, the model is not fitted but instead the function returns the model file and the data/initial values needed to fit the model outside of mvgam.
return_model_data	logical. If TRUE, the list of data needed to fit the model is returned, along with initial values for smooth and AR parameters, once the model is fitted. Helpful for users who wish to modify the model file to add other stochastic elements. Default is FALSE unless run_model == FALSE.
residuals	logical. Whether to compute series-level randomized quantile residuals. Default is TRUE. Set to FALSE to save time and reduce object size (can add later using add_residuals).
...	Other arguments to pass to mvgam

Details

Joint Species Distribution Models allow for responses of multiple species to be learned hierarchically, whereby responses to environmental variables in formula can be partially pooled and any latent, unmodelled residual associations can also be learned. In mvgam, both of these effects can be modelled with the full power of latent factor Hierarchical GAMs, providing unmatched flexibility to model full communities of species. When calling jsdgam, an initial State-Space model using trend = 'None' is set up and then modified to include the latent factors and their linear predictors. Consequently, you can inspect priors for these models using get_mvgam_priors by supplying the relevant formula, factor_formula, data and family arguments and keeping the default trend = 'None'.

In a JSDGAM, the expectation of response $Y_{ij}$ is modelled with

$g(\mu_{ij}) = X_i\beta + u_i\theta_j,$

where $g(.)$ is a known link function, $X$ is a design matrix of linear predictors (with associated $\beta$ coefficients), $u$ are $n_{lv}$-variate latent factors ($n_{lv}$<<$n_{species}$) and $\theta_j$ are species-specific loadings on the latent factors, respectively. The design matrix $X$ and $\beta$ coefficients are constructed and modelled using formula and can contain any of mvgam's predictor effects, including random intercepts and slopes, multidimensional penalized smooths, GP effects etc... The factor loadings $\theta_j$ are constrained for identifiability but can be used to reconstruct an estimate of the species' residual variance-covariance matrix using $\Theta \Theta'$ (see the example below and residual_cor() for details). The latent factors are further modelled using:

$u_i \sim \text{Normal}(Q_i\beta_{factor}, 1)$

where the second design matrix $Q$ and associated $\beta_{factor}$ coefficients are constructed and modelled using factor_formula. Again, the effects that make up this linear predictor can contain any of mvgam's allowed predictor effects, providing enormous flexibility for modelling species' communities.

Value

A list object of class mvgam containing model output, the text representation of the model file, the mgcv model output (for easily generating simulations at unsampled covariate values), Dunn-Smyth residuals for each species and key information needed for other functions in the package. See mvgam-class for details. Use methods(class = "mvgam") for an overview on available methods

Author(s)

Nicholas J Clark

References

Nicholas J Clark & Konstans Wells (2023). Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series. Methods in Ecology and Evolution. 14:3, 771-784.

David I Warton, F Guillaume Blanchet, Robert B O'Hara, Otso Ovaskainen, Sara Taskinen, Steven C Walker & Francis KC Hui (2015). So many variables: joint modeling in community ecology. Trends in Ecology & Evolution 30:12, 766-779.

See Also

mvgam(), residual_cor()

Examples

## Not run: 
# ========================================================================
# Example 1: Basic JSDGAM with Portal Data
# ========================================================================

# Fit a JSDGAM to the portal_data captures
mod <- jsdgam(
  formula = captures ~
    # Fixed effects of NDVI and mintemp, row effect as a GP of time
    ndvi_ma12:series + mintemp:series + gp(time, k = 15),
  factor_formula = ~ -1,
  data = portal_data,
  unit = time,
  species = series,
  family = poisson(),
  n_lv = 2,
  silent = 2,
  chains = 2
)

# Plot covariate effects
library(ggplot2); theme_set(theme_bw())
plot_predictions(
  mod,
  condition = c('ndvi_ma12', 'series', 'series')
)

plot_predictions(
  mod,
  condition = c('mintemp', 'series', 'series')
)

# A residual correlation plot
plot(residual_cor(mod))

# An ordination biplot can also be constructed
# from the factor scores and their loadings
if(requireNamespace('ggrepel', quietly = TRUE)){
  ordinate(mod, alpha = 0.7)
}

# ========================================================================
# Example 2: Advanced JSDGAM with Spatial Predictors
# ========================================================================

# Simulate latent count data for 500 spatial locations and 10 species
set.seed(0)
N_points <- 500
N_species <- 10

# Species-level intercepts (on the log scale)
alphas <- runif(N_species, 2, 2.25)

# Simulate a covariate and species-level responses to it
temperature <- rnorm(N_points)
betas <- runif(N_species, -0.5, 0.5)

# Simulate points uniformly over a space
lon <- runif(N_points, min = 150, max = 155)
lat <- runif(N_points, min = -20, max = -19)

# Set up spatial basis functions as a tensor product of lat and lon
sm <- mgcv::smoothCon(
  mgcv::te(lon, lat, k = 5),
  data = data.frame(lon, lat),
  knots = NULL
)[[1]]

# The design matrix for this smooth is in the 'X' slot
des_mat <- sm$X
dim(des_mat)

# Function to generate a random covariance matrix where all variables
# have unit variance (i.e. diagonals are all 1)
random_Sigma = function(N){
  L_Omega <- matrix(0, N, N);
  L_Omega[1, 1] <- 1;
  for (i in 2 : N) {
    bound <- 1;
    for (j in 1 : (i - 1)) {
      L_Omega[i, j] <- runif(1, -sqrt(bound), sqrt(bound));
      bound <- bound - L_Omega[i, j] ^ 2;
    }
    L_Omega[i, i] <- sqrt(bound);
  }
  Sigma <- L_Omega %*% t(L_Omega);
  return(Sigma)
}

# Simulate a variance-covariance matrix for the correlations among
# basis coefficients
Sigma <- random_Sigma(N = NCOL(des_mat))

# Now simulate the species-level basis coefficients hierarchically, where
# spatial basis function correlations are a convex sum of a base correlation
# matrix and a species-level correlation matrix
basis_coefs <- matrix(NA, nrow = N_species, ncol = NCOL(Sigma))
base_field <- mgcv::rmvn(1, mu = rep(0, NCOL(Sigma)), V = Sigma)
for(t in 1:N_species){
  corOmega <- (cov2cor(Sigma) * 0.7) +
    (0.3 * cov2cor(random_Sigma(N = NCOL(des_mat))))
  basis_coefs[t, ] <- mgcv::rmvn(1, mu = rep(0, NCOL(Sigma)), V = corOmega)
}

# Simulate the latent spatial processes
st_process <- do.call(rbind, lapply(seq_len(N_species), function(t){
  data.frame(
    lat = lat,
    lon = lon,
    species = paste0('species_', t),
    temperature = temperature,
    process = alphas[t] +
      betas[t] * temperature +
      des_mat %*% basis_coefs[t,]
  )
}))

# Now take noisy observations at some of the points (60)
obs_points <- sample(1:N_points, size = 60, replace = FALSE)
obs_points <- data.frame(
  lat = lat[obs_points],
  lon = lon[obs_points],
  site = 1:60
)

# Keep only the process data at these points
st_process %>%
  dplyr::inner_join(obs_points, by = c('lat', 'lon')) %>%
  # now take noisy Poisson observations of the process
  dplyr::mutate(count = rpois(NROW(.), lambda = exp(process))) %>%
  dplyr::mutate(species = factor(
    species,
    levels = paste0('species_', 1:N_species)
  )) %>%
  dplyr::group_by(lat, lon) -> dat

# View the count distributions for each species
ggplot(dat, aes(x = count)) +
  geom_histogram() +
  facet_wrap(~ species, scales = 'free')

ggplot(dat, aes(x = lon, y = lat, col = log(count + 1))) +
  geom_point(size = 2.25) +
  facet_wrap(~ species, scales = 'free') +
  scale_color_viridis_c()

# ------------------------------------------------------------------------
# Model Fitting with Custom Priors
# ------------------------------------------------------------------------

# Inspect default priors for a joint species model with three spatial factors
priors <- get_mvgam_priors(
  formula = count ~
    # Environmental model includes random slopes for
    # a linear effect of temperature
    s(species, bs = 're', by = temperature),

  # Each factor estimates a different nonlinear spatial process, using
  # 'by = trend' as in other mvgam State-Space models
  factor_formula = ~ gp(lon, lat, k = 6, by = trend) - 1,
  n_lv = 3,

  # The data and grouping variables
  data = dat,
  unit = site,
  species = species,

  # Poisson observations
  family = poisson()
)
head(priors)

# Fit a JSDM that estimates hierarchical temperature responses
# and that uses three latent spatial factors
mod <- jsdgam(
  formula = count ~
    # Environmental model includes random slopes for a
    # linear effect of temperature
    s(species, bs = 're', by = temperature),

  # Each factor estimates a different nonlinear spatial process, using
  # 'by = trend' as in other mvgam State-Space models
  factor_formula = ~ gp(lon, lat, k = 6, by = trend) - 1,
  n_lv = 3,

  # Change default priors for fixed random effect variances and
  # factor GP marginal deviations to standard normal
  priors = c(
    prior(std_normal(), class = sigma_raw),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend1`),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend2`),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend3`)
  ),

  # The data and the grouping variables
  data = dat,
  unit = site,
  species = species,

  # Poisson observations
  family = poisson(),
  chains = 2,
  silent = 2
)

# ------------------------------------------------------------------------
# Model Visualization and Diagnostics
# ------------------------------------------------------------------------

# Plot the implicit species-level intercept estimates
plot_predictions(mod, condition = 'species', type = 'link')

# Plot species' hierarchical responses to temperature
plot_predictions(
  mod,
  condition = c('temperature', 'species', 'species'),
  type = 'link'
)

# Plot posterior median estimates of the latent spatial factors
plot(mod, type = 'smooths', trend_effects = TRUE)

# Or using gratia, if you have it installed
if(requireNamespace('gratia', quietly = TRUE)){
  gratia::draw(mod, trend_effects = TRUE, dist = 0)
}

# Plot species' randomized quantile residual distributions
# as a function of latitude
pp_check(
  mod,
  type = 'resid_ribbon_grouped',
  group = 'species',
  x = 'lat',
  ndraws = 200
)

# ------------------------------------------------------------------------
# Residual Correlation Analysis
# ------------------------------------------------------------------------

# Calculate residual spatial correlations
post_cors <- residual_cor(mod)
names(post_cors)

# Look at lower and upper credible interval estimates for
# some of the estimated correlations
post_cors$cor[1:5, 1:5]
post_cors$cor_upper[1:5, 1:5]
post_cors$cor_lower[1:5, 1:5]

# Plot of the posterior median correlations for those estimated
# to be non-zero
plot(post_cors, cluster = TRUE)

# An ordination biplot can also be constructed
# from the factor scores and their loadings
if(requireNamespace('ggrepel', quietly = TRUE)){
  ordinate(mod)
}

# ------------------------------------------------------------------------
# Model Validation and Prediction
# ------------------------------------------------------------------------

# Posterior predictive checks and ELPD-LOO can ascertain model fit
pp_check(
  mod,
  type = "pit_ecdf_grouped",
  group = "species",
  ndraws = 200
)
loo(mod)

# Forecast log(counts) for entire region (site value doesn't matter as long
# as each spatial location has a different and unique site identifier);
# note this calculation takes a few minutes because of the need to calculate
# draws from the stochastic latent factors
newdata <- st_process %>%
  dplyr::mutate(species = factor(
    species,
    levels = paste0('species_', 1:N_species)
  )) %>%
  dplyr::group_by(lat, lon) %>%
  dplyr::mutate(site = dplyr::cur_group_id()) %>%
  dplyr::ungroup()
preds <- predict(mod, newdata = newdata)

# Plot the median log(count) predictions on a grid
newdata$log_count <- preds[,1]
ggplot(newdata, aes(x = lon, y = lat, col = log_count)) +
  geom_point(size = 1.5) +
  facet_wrap(~ species, scales = 'free') +
  scale_color_viridis_c() +
  theme_classic()

# Not needed for general use; cleans up connections for automated testing
closeAllConnections()

## End(Not run)

---

Approximate leave-future-out cross-validation of fitted mvgam objects

Description

Approximate leave-future-out cross-validation of fitted mvgam objects

Usage

lfo_cv(object, ...)

## S3 method for class 'mvgam'
lfo_cv(
  object,
  data,
  min_t,
  fc_horizon = 1,
  pareto_k_threshold = 0.7,
  silent = 1,
  ...
)

Arguments

object	list object of class mvgam. See mvgam()
...	Ignored
data	A dataframe or list containing the model response variable and covariates required by the GAM formula. Should include columns: 'series' (character or factor index of the series IDs) 'time' (numeric index of the time point for each observation). Any other variables to be included in the linear predictor of formula must also be present
min_t	Integer specifying the minimum training time required before making predictions from the data. Default is either the 30th timepoint in the observational data, or whatever training time allows for at least 10 lfo-cv calculations, if possible. This value is essentially arbitrary so it is highly recommended to change it to something that is more suitable to the data and models being evaluated.
fc_horizon	Integer specifying the number of time steps ahead for evaluating forecasts
pareto_k_threshold	Proportion specifying the threshold over which the Pareto shape parameter is considered unstable, triggering a model refit. Default is 0.7
silent	Verbosity level between 0 and 2. If 1 (the default), most of the informational messages of compiler and sampler are suppressed. If 2, even more messages are suppressed. The actual sampling progress is still printed. Set refresh = 0 to turn this off as well. If using backend = "rstan" you can also set open_progress = FALSE to prevent opening additional progress bars.

Details

Approximate leave-future-out cross-validation uses an expanding training window scheme to evaluate a model on its forecasting ability. The steps used in this function mirror those laid out in the lfo vignette from the loo package, written by Paul Bürkner, Jonah Gabry, Aki Vehtari. First, we refit the model using the first min_t observations to perform a single exact fc_horizon-ahead forecast step. This forecast is evaluated against the min_t + fc_horizon out of sample observations using the Expected Log Predictive Density (ELPD). Next, we approximate each successive round of expanding window forecasts by moving forward one step at a time ⁠for i in 1:N_evaluations⁠ and re-weighting draws from the model's posterior predictive distribution using Pareto Smoothed Importance Sampling (PSIS). In each iteration i, PSIS weights are obtained for the next observation that would have been included in the model if we had re-fit (i.e. the last observation that would have been in the training data, or min_t + i). If these importance ratios are stable, we consider the approximation adequate and use the re-weighted posterior's forecast for evaluating the next holdout set of testing observations ((min_t + i + 1):(min_t + i + fc_horizon)). At some point the importance ratio variability will become too large and importance sampling will fail. This is indicated by the estimated shape parameter k of the generalized Pareto distribution crossing a certain threshold pareto_k_threshold. Only then do we refit the model using all of the observations up to the time of the failure. We then restart the process and iterate forward until the next refit is triggered (Bürkner et al. 2020).

Value

A list of class mvgam_lfo containing the approximate ELPD scores, the Pareto-k shape values and 'the specified pareto_k_threshold

Author(s)

Nicholas J Clark

References

Paul-Christian Bürkner, Jonah Gabry & Aki Vehtari (2020). Approximate leave-future-out cross-validation for Bayesian time series models Journal of Statistical Computation and Simulation. 90:14, 2499-2523.

See Also

forecast, score, compare_mvgams

Examples

## Not run: 
# Simulate from a Poisson-AR2 model with a seasonal smooth
set.seed(100)
dat <- sim_mvgam(T = 75,
                n_series = 1,
                prop_trend = 0.75,
                trend_model = 'AR2',
                family = poisson())

# Plot the time series
plot_mvgam_series(data = dat$data_train,
                 newdata = dat$data_test,
                 series = 1)

# Fit an appropriate model
mod_ar2 <- mvgam(y ~ s(season, bs = 'cc', k = 6),
               trend_model = AR(p = 2),
               family = poisson(),
               data = dat$data_train,
               newdata = dat$data_test,
               chains = 2,
               silent = 2)

# Fit a less appropriate model
mod_rw <- mvgam(y ~ s(season, bs = 'cc', k = 6),
              trend_model = RW(),
              family = poisson(),
              data = dat$data_train,
              newdata = dat$data_test,
              chains = 2,
              silent = 2)

# Compare Discrete Ranked Probability Scores for the testing period
fc_ar2 <- forecast(mod_ar2)
fc_rw <- forecast(mod_rw)
score_ar2 <- score(fc_ar2, score = 'drps')
score_rw <- score(fc_rw, score = 'drps')
sum(score_ar2$series_1$score)
sum(score_rw$series_1$score)

# Now use approximate leave-future-out CV to compare
# rolling forecasts; start at time point 40 to reduce
# computational time and to ensure enough data is available
# for estimating model parameters
lfo_ar2 <- lfo_cv(mod_ar2,
                 min_t = 40,
                 fc_horizon = 3,
                 silent = 2)
lfo_rw <- lfo_cv(mod_rw,
                min_t = 40,
                fc_horizon = 3,
                silent = 2)

# Plot Pareto-K values and ELPD estimates
plot(lfo_ar2)
plot(lfo_rw)

# Proportion of timepoints in which AR2 model gives better forecasts
length(which((lfo_ar2$elpds - lfo_rw$elpds) > 0)) /
      length(lfo_ar2$elpds)

# A higher total ELPD is preferred
lfo_ar2$sum_ELPD
lfo_rw$sum_ELPD

## End(Not run)

---

Compute pointwise Log-Likelihoods from fitted mvgam objects

Description

Compute pointwise Log-Likelihoods from fitted mvgam objects

Usage

## S3 method for class 'mvgam'
logLik(object, linpreds, newdata, family_pars, include_forecast = TRUE, ...)

Arguments

object	list object of class mvgam or jsdgam
linpreds	Optional matrix of linear predictor draws to use for calculating pointwise log-likelihoods.
newdata	Optional data.frame or list object specifying which series each column in linpreds belongs to. If linpreds is supplied, then newdata must also be supplied.
family_pars	Optional list containing posterior draws of family-specific parameters (i.e. shape, scale or overdispersion parameters). Required if linpreds and newdata are supplied.
include_forecast	Logical. If newdata were fed to the model to compute forecasts, should the log-likelihood draws for these observations also be returned. Defaults to TRUE.
...	Ignored

Value

A matrix of dimension ⁠n_samples x n_observations⁠ containing the pointwise log-likelihood draws for all observations in newdata. If no newdata is supplied, log-likelihood draws are returned for all observations that were originally fed to the model (training observations and, if supplied to the original model via the newdata argument in mvgam, testing observations).

Author(s)

Nicholas J Clark

Examples

## Not run: 
# Simulate some data and fit a model
simdat <- sim_mvgam(
  n_series = 1,
  trend_model = AR()
)

mod <- mvgam(
  y ~ s(season, bs = 'cc', k = 6),
  trend_model = AR(),
  data = simdat$data_train,
  chains = 2,
  silent = 2
)

# Extract log-likelihood values
lls <- logLik(mod)
str(lls)

## End(Not run)

---

LOO information criteria for mvgam models

Description

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

Usage

## S3 method for class 'mvgam'
loo(x, incl_dynamics = FALSE, ...)

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

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

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).

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).

Author(s)

Nicholas J Clark

Examples

## Not run: 
#--------------------------------------------------
# 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')

## End(Not run)

---

Calculate trend correlations based on latent factor loadings for mvgam models

Description

This function uses factor loadings from a fitted dynamic factor mvgam model to calculate temporal correlations among series' trends.

Usage

lv_correlations(object)

Arguments

object

list object of class mvgam that used latent factors, either with use_lv = TRUE or by supplying a trend_map. See mvgam() for details and for an example.

Details

Although this function will still work, it is now recommended to use residual_cor() to obtain residual correlation information in a more user-friendly format that allows for a deeper investigation of relationships among the time series.

Value

A list object containing the mean posterior correlations and the full array of posterior correlations.

See Also

residual_cor(), plot.mvgam_residcor()

Examples

## Not run: 
#--------------------------------------------------
# Fit a model that uses two AR(1) dynamic factors to model
# the temporal dynamics of the four rodent species in the portal_data
#--------------------------------------------------
mod <- mvgam(
  captures ~ series,
  trend_model = AR(),
  use_lv = TRUE,
  n_lv = 2,
  data = portal_data,
  chains = 2,
  silent = 2
)

# Plot the two dynamic factors
plot(mod, type = 'factors')

# Calculate correlations among the series
lvcors <- lv_correlations(mod)
names(lvcors)
lapply(lvcors, class)

# Recommended: use residual_cor() instead
lvcors <- residual_cor(mod)
names(lvcors)
lvcors$cor

# Plot credible correlations as a matrix
plot(lvcors, cluster = TRUE)

# Not needed for general use; cleans up connections for automated testing
closeAllConnections()

## End(Not run)

---

MCMC plots of mvgam parameters, as implemented in bayesplot

Description

Convenient way to call MCMC plotting functions implemented in the bayesplot package for mvgam models

Usage

## S3 method for class 'mvgam'
mcmc_plot(
  object,
  type = "intervals",
  variable = NULL,
  regex = FALSE,
  use_alias = TRUE,
  ...
)

Arguments

object	An R object typically of class brmsfit
type	The type of the plot. Supported types are (as names) hist, dens, hist_by_chain, dens_overlay, violin, intervals, areas, areas_ridges, combo, acf, acf_bar, trace, trace_highlight, scatter, hex, pairs, violin, rhat, rhat_hist, neff, neff_hist and nuts_energy. For an overview on the various plot types see MCMC-overview.
variable	Names of the variables (parameters) to plot, as given by a character vector or a regular expression (if regex = TRUE). By default, a hopefully not too large selection of variables is plotted.
regex	Logical; Indicates whether variable should be treated as regular expressions. Defaults to FALSE.
use_alias	Logical. If more informative names for parameters are available (i.e. for beta coefficients b or for smoothing parameters rho), replace the uninformative names with the more informative alias. Defaults to TRUE.
...	Additional arguments passed to the plotting functions. See MCMC-overview for more details.

Value

A ggplot object that can be further customized using the ggplot2 package.

See Also

mvgam_draws for an overview of some of the shortcut strings that can be used for argument variable

Examples

## Not run: 
simdat <- sim_mvgam(n_series = 1, trend_model = AR())
mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
             trend_model = AR(),
             noncentred = TRUE,
             data = simdat$data_train,
             chains = 2,
             silent = 2)
mcmc_plot(mod)
mcmc_plot(mod, type = 'neff_hist')
mcmc_plot(mod, variable = 'betas', type = 'areas')
mcmc_plot(mod, variable = 'trend_params', type = 'combo')

## End(Not run)

---

Extract model.frame from a fitted mvgam object

Description

Extract model.frame from a fitted mvgam object

Usage

## S3 method for class 'mvgam'
model.frame(formula, trend_effects = FALSE, ...)

## S3 method for class 'mvgam_prefit'
model.frame(formula, trend_effects = FALSE, ...)

Arguments

formula	a model formula or terms object or an R object.
trend_effects	logical, return the model.frame from the observation model (if FALSE) or from the underlying process model (if TRUE)
...	Ignored

Value

A matrix containing the fitted model frame

Author(s)

Nicholas J Clark

---

Monotonic splines in mvgam models

Description

Uses constructors from package splines2 to build monotonically increasing or decreasing splines. Details also in Wang & Yan (2021).

Usage

## S3 method for class 'moi.smooth.spec'
smooth.construct(object, data, knots)

## S3 method for class 'mod.smooth.spec'
smooth.construct(object, data, knots)

## S3 method for class 'moi.smooth'
Predict.matrix(object, data)

## S3 method for class 'mod.smooth'
Predict.matrix(object, data)

Arguments

object	A smooth specification object, usually generated by a term s(x, bs = "moi", ...) or s(x, bs = "mod", ...)
data	a list containing just the data (including any by variable) required by this term, with names corresponding to object$term (and object$by). The by variable is the last element.
knots	a list containing any knots supplied for basis setup — in same order and with same names as data. Can be NULL. See details for further information.

Details

The constructor is not normally called directly, but is rather used internally by mvgam. If they are not supplied then the knots of the spline are placed evenly throughout the covariate values to which the term refers: For example, if fitting 101 data with an 11 knot spline of x then there would be a knot at every 10th (ordered) x value. The spline is an implementation of the closed-form I-spline basis based on the recursion formula given by Ramsay (1988), in which the basis coefficients must be constrained to either be non-negative (for monotonically increasing functions) or non-positive (monotonically decreasing)

Take note that when using either monotonic basis, the number of basis functions k must be supplied as an even integer due to the manner in which monotonic basis functions are constructed

Value

An object of class "moi.smooth" or "mod.smooth". In addition to the usual elements of a smooth class documented under smooth.construct, this object will contain a slot called boundary that defines the endpoints beyond which the spline will begin extrapolating (extrapolation is flat due to the first order penalty placed on the smooth function)

Note

This constructor will result in a valid smooth if using a call to gam or bam, however the resulting functions will not be guaranteed to be monotonic because constraints on basis coefficients will not be enforced

Author(s)

Nicholas J Clark

References

Wang, Wenjie, and Jun Yan. "Shape-Restricted Regression Splines with R Package splines2." Journal of Data Science 19.3 (2021).

Ramsay, J. O. (1988). Monotone regression splines in action. Statistical Science, 3(4), 425–441.

Examples

## Not run: 
# Simulate data from a monotonically increasing function
set.seed(123123)

x <- runif(80) * 4 - 1
x <- sort(x)
f <- exp(4 * x) / (1 + exp(4 * x))
y <- f + rnorm(80) * 0.1
plot(x, y)

# A standard TRPS smooth doesn't capture monotonicity
library(mgcv)

mod_data <- data.frame(y = y, x = x)
mod <- gam(
  y ~ s(x, k = 16),
  data = mod_data,
  family = gaussian()
)

library(marginaleffects)
plot_predictions(
  mod,
  by = 'x',
  newdata = data.frame(
    x = seq(min(x) - 0.5, max(x) + 0.5, length.out = 100)
  ),
  points = 0.5
)

# Using the 'moi' basis in mvgam rectifies this
mod_data$time <- 1:NROW(mod_data)
mod2 <- mvgam(
  y ~ s(x, bs = 'moi', k = 18),
  data = mod_data,
  family = gaussian(),
  chains = 2,
  silent = 2
)

plot_predictions(
  mod2,
  by = 'x',
  newdata = data.frame(
    x = seq(min(x) - 0.5, max(x) + 0.5, length.out = 100)
  ),
  points = 0.5
)

plot(mod2, type = 'smooth', realisations = TRUE)

# 'by' terms that produce a different smooth for each level of the 'by'
# factor are also allowed

x <- runif(80) * 4 - 1
x <- sort(x)

# Two different monotonic smooths, one for each factor level
f <- exp(4 * x) / (1 + exp(4 * x))
f2 <- exp(3.5 * x) / (1 + exp(3 * x))
fac <- c(rep('a', 80), rep('b', 80))
y <- c(
  f + rnorm(80) * 0.1,
  f2 + rnorm(80) * 0.2
)

plot(x, y[1:80])
plot(x, y[81:160])

# Gather all data into a data.frame, including the factor 'by' variable
mod_data <- data.frame(y, x, fac = as.factor(fac))
mod_data$time <- 1:NROW(mod_data)

# Fit a model with different smooths per factor level
mod <- mvgam(
  y ~ s(x, bs = 'moi', by = fac, k = 8),
  data = mod_data,
  family = gaussian(),
  chains = 2,
  silent = 2
)

# Visualise the different monotonic functions
plot_predictions(
  mod,
  condition = c('x', 'fac', 'fac'),
  points = 0.5
)

plot(mod, type = 'smooth', realisations = TRUE)

# First derivatives (on the link scale) should never be
# negative for either factor level
(derivs <- slopes(
  mod,
  variables = 'x',
  by = c('x', 'fac'),
  type = 'link'
))

all(derivs$estimate > 0)

## End(Not run)

---

Fit a Bayesian Dynamic GAM to Univariate or Multivariate Time Series

Description

This function estimates the posterior distribution for Generalised Additive Models (GAMs) that can include smooth spline functions, specified in the GAM formula, as well as latent temporal processes, specified by trend_model.

Further modelling options include State-Space representations to allow covariates and dynamic processes to occur on the latent 'State' level while also capturing observation-level effects. Prior specifications are flexible and explicitly encourage users to apply prior distributions that actually reflect their beliefs.

In addition, model fits can easily be assessed and compared with posterior predictive checks, forecast comparisons and leave-one-out / leave-future-out cross-validation.

Usage

mvgam(
  formula,
  trend_formula,
  knots,
  trend_knots,
  trend_model = "None",
  noncentred = FALSE,
  family = poisson(),
  share_obs_params = FALSE,
  data,
  newdata,
  use_lv = FALSE,
  n_lv,
  trend_map,
  priors,
  run_model = TRUE,
  prior_simulation = FALSE,
  residuals = TRUE,
  return_model_data = FALSE,
  backend = getOption("brms.backend", "cmdstanr"),
  algorithm = getOption("brms.algorithm", "sampling"),
  control = list(max_treedepth = 10, adapt_delta = 0.8),
  chains = 4,
  burnin = 500,
  samples = 500,
  thin = 1,
  parallel = TRUE,
  threads = 1,
  save_all_pars = FALSE,
  silent = 1,
  autoformat = TRUE,
  refit = FALSE,
  lfo = FALSE,
  ...
)

Arguments

formula	A formula object specifying the GAM observation model formula. These are exactly like the formula for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, nonparametric gp() terms and offsets using offset(), can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these). In nmix() family models, the formula is used to set up a linear predictor for the detection probability. Details of the formula syntax used by mvgam can be found in mvgam_formulae
trend_formula	An optional formula object specifying the GAM process model formula. If supplied, a linear predictor will be modelled for the latent trends to capture process model evolution separately from the observation model. Important notes: Should not have a response variable specified on the left-hand side (e.g., ~ season + s(year)) Use trend instead of series for effects that vary across time series Only available for RW(), AR() and VAR() trend models In nmix() family models, sets up linear predictor for latent abundance Consider dropping one intercept using - 1 convention to avoid estimation challenges
knots	An optional list containing user specified knot values for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied. Different terms can use different numbers of knots, unless they share a covariate.
trend_knots	As for knots above, this is an optional list of knot values for smooth functions within the trend_formula.
trend_model	character or function specifying the time series dynamics for the latent trend. Available options: None: No latent trend component (GAM component only, like gam) ZMVN or ZMVN(): Zero-Mean Multivariate Normal (Stan only) 'RW' or RW(): Random Walk 'AR1', 'AR2', 'AR3' or AR(p = 1, 2, 3): Autoregressive models 'CAR1' or CAR(p = 1): Continuous-time AR (Ornstein–Uhlenbeck process) 'VAR1' or VAR(): Vector Autoregressive (Stan only) 'PWlogistic', 'PWlinear' or PW(): Piecewise trends (Stan only) 'GP' or GP(): Gaussian Process with squared exponential kernel (Stan only) Additional features: Moving average and/or correlated process error terms available for most types (e.g., RW(cor = TRUE) for multivariate Random Walk) Hierarchical correlations possible for structured data See mvgam_trends for details and ZMVN() for examples
noncentred	logical. Use non-centred parameterisation for autoregressive trend models? Can improve efficiency by avoiding degeneracies in latent dynamic random effects estimation. Benefits vary by model - highly informative data may perform worse with this option. Available for RW(), AR(), CAR(), or trend = 'None' with trend_formula. Not available for moving average or correlated error models.
family	family specifying the exponential observation family for the series. Supported families: gaussian(): Real-valued data betar(): Proportional data on ⁠(0,1)⁠ lognormal(): Non-negative real-valued data student_t(): Real-valued data Gamma(): Non-negative real-valued data bernoulli(): Binary data poisson(): Count data (default) nb(): Overdispersed count data binomial(): Count data with imperfect detection when number of trials is known (use cbind() to bind observations and trials) beta_binomial(): As binomial() but allows for overdispersion nmix(): Count data with imperfect detection when number of trials is unknown (State-Space N-Mixture model with Poisson latent states and Binomial observations) See mvgam_families for more details.
share_obs_params	logical. If TRUE and the family has additional family-specific observation parameters (e.g., variance components, dispersion parameters), these will be shared across all outcome variables. Useful when multiple outcomes share properties. Default is FALSE.
data	A dataframe or list containing the model response variable and covariates required by the GAM formula and optional trend_formula. Required columns for most models: series: A factor index of the series IDs (number of levels should equal number of unique series labels) time: numeric or integer index of time points. For most dynamic trend types, time should be measured in discrete, regularly spaced intervals (i.e., c(1, 2, 3, ...)). Irregular spacing is allowed for trend_model = CAR(1), but zero intervals are adjusted to 1e-12 to prevent sampling errors. Special cases: Models with hierarchical temporal correlation (e.g., AR(gr = region, subgr = species)) should NOT include a series identifier Models without temporal dynamics (trend_model = 'None' or trend_model = ZMVN()) don't require a time variable
newdata	Optional dataframe or list of test data containing the same variables as in data. If included, observations in variable y will be set to NA when fitting the model so that posterior simulations can be obtained.
use_lv	logical. If TRUE, use dynamic factors to estimate series' latent trends in a reduced dimension format. Only available for RW(), AR() and GP() trend models. Default is FALSE. See lv_correlations for examples.
n_lv	integer specifying the number of latent dynamic factors to use if use_lv == TRUE. Cannot exceed n_series. Default is min(2, floor(n_series / 2)).
trend_map	Optional data.frame specifying which series should depend on which latent trends. Enables multiple series to depend on the same latent trend process with different observation processes. Required structure: Column series: Single unique entry for each series (matching factor levels in data) Column trend: Integer values indicating which trend each series depends on Notes: Sets up latent factor model by enabling use_lv = TRUE Process model intercept is NOT automatically suppressed Not yet supported for continuous time models (CAR())
priors	An optional data.frame with prior definitions or, preferably, a vector of brmsprior objects (see prior()). See get_mvgam_priors() and Details for more information.
run_model	logical. If FALSE, the model is not fitted but instead the function returns the model file and the data/initial values needed to fit the model outside of mvgam.
prior_simulation	logical. If TRUE, no observations are fed to the model, and instead simulations from prior distributions are returned.
residuals	logical. Whether to compute series-level randomized quantile residuals. Default is TRUE. Set to FALSE to save time and reduce object size (can add later using add_residuals).
return_model_data	logical. If TRUE, the list of data needed to fit the model is returned, along with initial values for smooth and AR parameters, once the model is fitted. Helpful for users who wish to modify the model file to add other stochastic elements. Default is FALSE unless run_model == FALSE.
backend	Character string naming the package for Stan model fitting. Options are "cmdstanr" (default) or "rstan". Can be set globally via "brms.backend" option. See https://mc-stan.org/rstan/ and https://mc-stan.org/cmdstanr/ for details.
algorithm	Character string naming the estimation approach: "sampling": MCMC (default) "meanfield": Variational inference with factorized normal distributions "fullrank": Variational inference with multivariate normal distribution "laplace": Laplace approximation (cmdstanr only) "pathfinder": Pathfinder algorithm (cmdstanr only) Can be set globally via "brms.algorithm" option. Limited testing suggests "meanfield" performs best among non-MCMC approximations for dynamic GAMs.
control	Named list for controlling sampler behaviour. Valid elements include max_treedepth, adapt_delta and init.
chains	integer specifying the number of parallel chains for the model. Ignored for variational inference algorithms.
burnin	integer specifying the number of warmup iterations to tune sampling algorithms. Ignored for variational inference algorithms.
samples	integer specifying the number of post-warmup iterations for sampling the posterior distribution.
thin	Thinning interval for monitors. Ignored for variational inference algorithms.
parallel	logical specifying whether to use multiple cores for parallel MCMC simulation. If TRUE, uses min(c(chains, parallel::detectCores() - 1)) cores.
threads	integer. Experimental option for within-chain parallelisation in Stan using reduce_sum. Recommended only for experienced Stan users with slow models. Currently works for all families except nmix() and when using Cmdstan backend.
save_all_pars	logical. Save draws from all variables defined in Stan's parameters block. Default is FALSE.
silent	Verbosity level between 0 and 2. If 1 (default), most informational messages are suppressed. If 2, even more messages are suppressed. Sampling progress is still printed - set refresh = 0 to disable. For backend = "rstan", also set open_progress = FALSE to prevent additional progress bars.
autoformat	logical. Use stanc parser to automatically format Stan code and check for deprecations. For development purposes - leave as TRUE.
refit	logical. Indicates whether this is a refit called using update.mvgam(). Users should leave as FALSE.
lfo	logical. Indicates whether this is part of lfo_cv.mvgam call. Returns lighter model version for speed. Users should leave as FALSE.
...	Further arguments passed to Stan: For backend = "rstan": passed to sampling() or vb() For backend = "cmdstanr": passed to cmdstanr::sample, cmdstanr::variational, cmdstanr::laplace or cmdstanr::pathfinder methods

Details

Dynamic GAMs are useful when we wish to predict future values from time series that show temporal dependence but we do not want to rely on extrapolating from a smooth term (which can sometimes lead to unpredictable and unrealistic behaviours). In addition, smooths can often try to wiggle excessively to capture any autocorrelation that is present in a time series, which exacerbates the problem of forecasting ahead.

As GAMs are very naturally viewed through a Bayesian lens, and we often must model time series that show complex distributional features and missing data, parameters for mvgam models are estimated in a Bayesian framework using Markov Chain Monte Carlo by default.

Getting Started Resources:

• General overview: vignette("mvgam_overview") and vignette("data_in_mvgam")

• Full list of vignettes: vignette(package = "mvgam")

• Real-world examples: mvgam_use_cases

• Quick reference: mvgam cheatsheet

Value

A list object of class mvgam containing model output, the text representation of the model file, the mgcv model output (for easily generating simulations at unsampled covariate values), Dunn-Smyth residuals for each series and key information needed for other functions in the package. See mvgam-class for details. Use methods(class = "mvgam") for an overview on available methods.

Model Specification Details

Formula Syntax: Details of the formula syntax used by mvgam can be found in mvgam_formulae. Note that it is possible to supply an empty formula where there are no predictors or intercepts in the observation model (i.e. y ~ 0 or y ~ -1). In this case, an intercept-only observation model will be set up but the intercept coefficient will be fixed at zero. This can be handy if you wish to fit pure State-Space models where the variation in the dynamic trend controls the average expectation, and/or where intercepts are non-identifiable (as in piecewise trends).

Families and Link Functions: Details of families supported by mvgam can be found in mvgam_families.

Trend Models: Details of latent error process models supported by mvgam can be found in mvgam_trends.

Prior Specifications

Default priors for intercepts and any variance parameters are chosen to be vaguely informative, but these should always be checked by the user. Prior distributions for most important model parameters can be altered (see get_mvgam_priors() for details). Note that latent trends are estimated on the link scale so choose priors accordingly.

However more control over the model specification can be accomplished by setting run_model = FALSE and then editing the model code (found in the model_file slot in the returned object) before running the model using either rstan or cmdstanr. This is encouraged for complex modelling tasks.

Important: No priors are formally checked to ensure they are in the right syntax so it is up to the user to ensure these are correct.

Model Components

Random Effects: For any smooth terms using the random effect basis (smooth.construct.re.smooth.spec), a non-centred parameterisation is automatically employed to avoid degeneracies that are common in hierarchical models. Note however that centred versions may perform better for series that are particularly informative, so as with any foray into Bayesian modelling, it is worth building an understanding of the model's assumptions and limitations by following a principled workflow. Also note that models are parameterised using drop.unused.levels = FALSE in jagam to ensure predictions can be made for all levels of the supplied factor variable.

Observation Level Parameters: When more than one series is included in data and an observation family that contains more than one parameter is used, additional observation family parameters (i.e. phi for nb() or sigma for gaussian()) are by default estimated independently for each series. But if you wish for the series to share the same observation parameters, set share_obs_params = TRUE.

Model Diagnostics

Residuals: For each series, randomized quantile (i.e. Dunn-Smyth) residuals are calculated for inspecting model diagnostics. If the fitted model is appropriate then Dunn-Smyth residuals will be standard normal in distribution and no autocorrelation will be evident. When a particular observation is missing, the residual is calculated by comparing independent draws from the model's posterior distribution.

Computational Backend

Using Stan: mvgam is primarily designed to use Hamiltonian Monte Carlo for parameter estimation via the software Stan (using either the cmdstanr or rstan interface). There are great advantages when using Stan over Gibbs / Metropolis Hastings samplers, which includes the option to estimate nonlinear effects via Hilbert space approximate Gaussian Processes, the availability of a variety of inference algorithms (i.e. variational inference, laplacian inference etc...) and capabilities to enforce stationarity for complex Vector Autoregressions.

Because of the many advantages of Stan over JAGS, further development of the package will only be applied to Stan. This includes the planned addition of more response distributions, plans to handle zero-inflation, and plans to incorporate a greater variety of trend models. Users are strongly encouraged to opt for Stan over JAGS in any proceeding workflows.

Recommended Workflow

How to Start: The mvgam cheatsheet is a good starting place if you are just learning to use the package. It gives an overview of the package's key functions and objects, as well as providing a reasonable workflow that new users can follow.

Recommended Steps:

1. Data Preparation: Check that your data are in a suitable tidy format for mvgam modeling (see the data formatting vignette for guidance)

2. Data Exploration: Inspect features of the data using plot_mvgam_series. Now is also a good time to familiarise yourself with the package's example workflows that are detailed in the vignettes:

   • Getting started vignette

   • Shared latent states vignette

   • Time-varying effects vignette

   • State-Space models vignette

   • "Fitting N-mixture models in mvgam"

   • "Joint Species Distribution Models in mvgam"

   • "Incorporating time-varying seasonality in forecast models"

   • "Temporal autocorrelation in GAMs and the mvgam package"

3. Model Structure: Carefully think about how to structure linear predictor effects (i.e. smooth terms using s(), te() or ti(), GPs using gp(), dynamic time-varying effects using dynamic(), and parametric terms), latent temporal trend components (see mvgam_trends) and the appropriate observation family (see mvgam_families). Use get_mvgam_priors() to see default prior distributions for stochastic parameters.

4. Prior Specification: Change default priors using appropriate prior knowledge (see prior()). When using State-Space models with a trend_formula, pay particular attention to priors for any variance parameters such as process errors and observation errors. Default priors on these parameters are chosen to be vaguely informative and to avoid zero (using Inverse Gamma priors), but more informative priors will often help with model efficiency and convergence.

5. Model Fitting: Fit the model using either Hamiltonian Monte Carlo or an approximation algorithm (i.e. change the backend argument) and use summary.mvgam(), conditional_effects.mvgam(), mcmc_plot.mvgam(), pp_check.mvgam(), pairs.mvgam() and plot.mvgam() to inspect / interrogate the model.

6. Model Comparison: Update the model as needed and use loo_compare.mvgam() for in-sample model comparisons, or alternatively use forecast.mvgam(), lfo_cv.mvgam() and score.mvgam_forecast() to compare models based on out-of-sample forecasts (see the forecast evaluation vignette for guidance).

7. Inference and Prediction: When satisfied with the model structure, use predict.mvgam(), plot_predictions() and/or plot_slopes() for more targeted simulation-based inferences (see "How to interpret and report nonlinear effects from Generalized Additive Models" for some guidance on interpreting GAMs). For time series models, use hindcast.mvgam(), fitted.mvgam(), augment.mvgam() and forecast.mvgam() to inspect posterior hindcast / forecast distributions.

8. Documentation: Use how_to_cite() to obtain a scaffold methods section (with full references) to begin describing this model in scientific publications.

Author(s)

Nicholas J Clark

References

Nicholas J Clark & Konstans Wells (2023). Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series. Methods in Ecology and Evolution. 14:3, 771-784.

Nicholas J Clark, SK Morgan Ernest, Henry Senyondo, Juniper Simonis, Ethan P White, Glenda M Yenni, KANK Karunarathna (2025). Beyond single-species models: leveraging multispecies forecasts to navigate the dynamics of ecological predictability. PeerJ. 13:e18929 https://doi.org/10.7717/peerj.18929

See Also

jagam(), gam(), gam.models, get_mvgam_priors(), jsdgam(), hindcast.mvgam(), forecast.mvgam(), predict.mvgam()

Examples

## Not run: 
# =============================================================================
# Basic Multi-Series Time Series Modeling
# =============================================================================

# Simulate three time series that have shared seasonal dynamics,
# independent AR(1) trends, and Poisson observations
set.seed(0)
dat <- sim_mvgam(
  T = 80,
  n_series = 3,
  mu = 2,
  trend_model = AR(p = 1),
  prop_missing = 0.1,
  prop_trend = 0.6
)

# Plot key summary statistics for a single series
plot_mvgam_series(data = dat$data_train, series = 1)

# Plot all series together
plot_mvgam_series(data = dat$data_train, series = "all")

# Formulate a model using Stan where series share a cyclic smooth for
# seasonality and each series has an independent AR1 temporal process.
# Note that 'noncentred = TRUE' will likely give performance gains.
# Set run_model = FALSE to inspect the returned objects
mod1 <- mvgam(
  formula = y ~ s(season, bs = "cc", k = 6),
  data = dat$data_train,
  trend_model = AR(),
  family = poisson(),
  noncentred = TRUE,
  run_model = FALSE
)

# View the model code in Stan language
stancode(mod1)

# View the data objects needed to fit the model in Stan
sdata1 <- standata(mod1)
str(sdata1)

# Now fit the model
mod1 <- mvgam(
  formula = y ~ s(season, bs = "cc", k = 6),
  data = dat$data_train,
  trend_model = AR(),
  family = poisson(),
  noncentred = TRUE,
  chains = 2,
  silent = 2
)

# Extract the model summary
summary(mod1)

# Plot the historical trend and hindcast distributions for one series
hc_trend <- hindcast(mod1, type = "trend")
plot(hc_trend)

hc_predicted <- hindcast(mod1, type = "response")
plot(hc_predicted)

# Residual diagnostics
plot(mod1, type = "residuals", series = 1)
resids <- residuals(mod1)
str(resids)

# Fitted values and residuals can be added directly to the training data
augment(mod1)

# Compute the forecast using covariate information in data_test
fc <- forecast(mod1, newdata = dat$data_test)
str(fc)
fc_summary <- summary(fc)
head(fc_summary, 12)
plot(fc)

# Plot the estimated seasonal smooth function
plot(mod1, type = "smooths")

# Plot estimated first derivatives of the smooth
plot(mod1, type = "smooths", derivatives = TRUE)

# Plot partial residuals of the smooth
plot(mod1, type = "smooths", residuals = TRUE)

# Plot posterior realisations for the smooth
plot(mod1, type = "smooths", realisations = TRUE)

# Plot conditional response predictions using marginaleffects
conditional_effects(mod1)
plot_predictions(mod1, condition = "season", points = 0.5)

# Generate posterior predictive checks using bayesplot
pp_check(mod1)

# Extract observation model beta coefficient draws as a data.frame
beta_draws_df <- as.data.frame(mod1, variable = "betas")
head(beta_draws_df)
str(beta_draws_df)

# Investigate model fit
mc.cores.def <- getOption("mc.cores")
options(mc.cores = 1)
loo(mod1)
options(mc.cores = mc.cores.def)


# =============================================================================
# Vector Autoregressive (VAR) Models
# =============================================================================

# Fit a model to the portal time series that uses a latent
# Vector Autoregression of order 1
mod <- mvgam(
  formula = captures ~ -1,
  trend_formula = ~ trend,
  trend_model = VAR(cor = TRUE),
  family = poisson(),
  data = portal_data,
  chains = 2,
  silent = 2
)

# Plot the autoregressive coefficient distributions;
# use 'dir = "v"' to arrange the order of facets correctly
mcmc_plot(
  mod,
  variable = 'A',
  regex = TRUE,
  type = 'hist',
  facet_args = list(dir = 'v')
)

# Plot the process error variance-covariance matrix in the same way
mcmc_plot(
  mod,
  variable = 'Sigma',
  regex = TRUE,
  type = 'hist',
  facet_args = list(dir = 'v')
)

# Calculate Generalized Impulse Response Functions for each series
irfs <- irf(
  mod,
  h = 12,
  cumulative = FALSE
)

# Plot some of them
plot(irfs, series = 1)
plot(irfs, series = 2)

# Calculate forecast error variance decompositions for each series
fevds <- fevd(mod, h = 12)

# Plot median contributions to forecast error variance
plot(fevds)


# =============================================================================
# Dynamic Factor Models
# =============================================================================

# Now fit a model that uses two RW dynamic factors to model
# the temporal dynamics of the four rodent species
mod <- mvgam(
  captures ~ series,
  trend_model = RW(),
  use_lv = TRUE,
  n_lv = 2,
  data = portal_data,
  chains = 2,
  silent = 2
)

# Plot the factors
plot(mod, type = 'factors')

# Plot the hindcast distributions
hcs <- hindcast(mod)
plot(hcs, series = 1)
plot(hcs, series = 2)
plot(hcs, series = 3)
plot(hcs, series = 4)

# Use residual_cor() to calculate temporal correlations among the series
# based on the factor loadings
lvcors <- residual_cor(mod)
names(lvcors)
lvcors$cor

# For those correlations whose credible intervals did not include
# zero, plot them as a correlation matrix (all other correlations
# are shown as zero on this plot)
plot(lvcors, cluster = TRUE)


# =============================================================================
# Shared Latent Trends with Custom Trend Mapping
# =============================================================================

# Example of supplying a trend_map so that some series can share
# latent trend processes
sim <- sim_mvgam(n_series = 3)
mod_data <- sim$data_train

# Here, we specify only two latent trends; series 1 and 2 share a trend,
# while series 3 has its own unique latent trend
trend_map <- data.frame(
  series = unique(mod_data$series),
  trend = c(1, 1, 2)
)

# Fit the model using AR1 trends
mod <- mvgam(
  formula = y ~ s(season, bs = "cc", k = 6),
  trend_map = trend_map,
  trend_model = AR(),
  data = mod_data,
  return_model_data = TRUE,
  chains = 2,
  silent = 2
)

# The mapping matrix is now supplied as data to the model in the 'Z' element
mod$model_data$Z

# The first two series share an identical latent trend; the third is different
plot(residual_cor(mod))
plot(mod, type = "trend", series = 1)
plot(mod, type = "trend", series = 2)
plot(mod, type = "trend", series = 3)


# =============================================================================
# Time-Varying (Dynamic) Coefficients
# =============================================================================

# Example of how to use dynamic coefficients
# Simulate a time-varying coefficient for the effect of temperature
set.seed(123)
N <- 200
beta_temp <- vector(length = N)
beta_temp[1] <- 0.4
for (i in 2:N) {
  beta_temp[i] <- rnorm(1, mean = beta_temp[i - 1] - 0.0025, sd = 0.05)
}
plot(beta_temp)

# Simulate a covariate called 'temp'
temp <- rnorm(N, sd = 1)

# Simulate some noisy Gaussian observations
out <- rnorm(N,
  mean = 4 + beta_temp * temp,
  sd = 0.5
)

# Gather necessary data into a data.frame; split into training / testing
data <- data.frame(out, temp, time = seq_along(temp))
data_train <- data[1:180, ]
data_test <- data[181:200, ]

# Fit the model using the dynamic() function
mod <- mvgam(
  formula = out ~ dynamic(
    temp,
    scale = FALSE,
    k = 40
  ),
  family = gaussian(),
  data = data_train,
  newdata = data_test,
  chains = 2,
  silent = 2
)

# Inspect the model summary, forecast and time-varying coefficient distribution
summary(mod)
plot(mod, type = "smooths")
fc <- forecast(mod, newdata = data_test)
plot(fc)

# Propagating the smooth term shows how the coefficient is expected to evolve
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 180, lty = "dashed", lwd = 2)
points(beta_temp, pch = 16)


# =============================================================================
# Working with Offset Terms
# =============================================================================

# Example showing how to incorporate an offset; simulate some count data
# with different means per series
set.seed(100)
dat <- sim_mvgam(
  prop_trend = 0,
  mu = c(0, 2, 2),
  seasonality = "hierarchical"
)

# Add offset terms to the training and testing data
dat$data_train$offset <- 0.5 * as.numeric(dat$data_train$series)
dat$data_test$offset <- 0.5 * as.numeric(dat$data_test$series)

# Fit a model that includes the offset in the linear predictor as well as
# hierarchical seasonal smooths
mod <- mvgam(
  formula = y ~ offset(offset) +
    s(series, bs = "re") +
    s(season, bs = "cc") +
    s(season, by = series, m = 1, k = 5),
  data = dat$data_train,
  chains = 2,
  silent = 2
)

# Inspect the model file to see the modification to the linear predictor (eta)
stancode(mod)

# Forecasts for the first two series will differ in magnitude
fc <- forecast(mod, newdata = dat$data_test)
plot(fc, series = 1, ylim = c(0, 75))
plot(fc, series = 2, ylim = c(0, 75))

# Changing the offset for the testing data should lead to changes in
# the forecast
dat$data_test$offset <- dat$data_test$offset - 2
fc <- forecast(mod, newdata = dat$data_test)
plot(fc)

# Relative Risks can be computed by fixing the offset to the same value
# for each series
dat$data_test$offset <- rep(1, NROW(dat$data_test))
preds_rr <- predict(mod,
  type = "link",
  newdata = dat$data_test,
  summary = FALSE
)
series1_inds <- which(dat$data_test$series == "series_1")
series2_inds <- which(dat$data_test$series == "series_2")

# Relative Risks are now more comparable among series
layout(matrix(1:2, ncol = 2))
plot(preds_rr[1, series1_inds],
  type = "l", col = "grey75",
  ylim = range(preds_rr),
  ylab = "Series1 Relative Risk", xlab = "Time"
)
for (i in 2:50) {
  lines(preds_rr[i, series1_inds], col = "grey75")
}

plot(preds_rr[1, series2_inds],
  type = "l", col = "darkred",
  ylim = range(preds_rr),
  ylab = "Series2 Relative Risk", xlab = "Time"
)
for (i in 2:50) {
  lines(preds_rr[i, series2_inds], col = "darkred")
}
layout(1)


# =============================================================================
# Binomial Family Models
# =============================================================================

# Example showcasing how cbind() is needed for Binomial observations
# Simulate two time series of Binomial trials
trials <- sample(c(20:25), 50, replace = TRUE)
x <- rnorm(50)
detprob1 <- plogis(-0.5 + 0.9 * x)
detprob2 <- plogis(-0.1 - 0.7 * x)
dat <- rbind(
  data.frame(
    y = rbinom(n = 50, size = trials, prob = detprob1),
    time = 1:50,
    series = "series1",
    x = x,
    ntrials = trials
  ),
  data.frame(
    y = rbinom(n = 50, size = trials, prob = detprob2),
    time = 1:50,
    series = "series2",
    x = x,
    ntrials = trials
  )
)
dat <- dplyr::mutate(dat, series = as.factor(series))
dat <- dplyr::arrange(dat, time, series)
plot_mvgam_series(data = dat, series = "all")

# Fit a model using the binomial() family; must specify observations
# and number of trials in the cbind() wrapper
mod <- mvgam(
  formula = cbind(y, ntrials) ~ series + s(x, by = series),
  family = binomial(),
  data = dat,
  chains = 2,
  silent = 2
)
summary(mod)
pp_check(mod,
  type = "bars_grouped",
  group = "series", ndraws = 50
)
pp_check(mod,
  type = "ecdf_overlay_grouped",
  group = "series", ndraws = 50
)
conditional_effects(mod, type = "link")

# To view predictions on the probability scale,
# use ntrials = 1 in datagrid()
plot_predictions(
  mod,
  by = c('x', 'series'),
  newdata = datagrid(
    x = runif(100, -2, 2),
    series = unique,
    ntrials = 1
  ),
  type = 'expected'
)

# Not needed for general use; cleans up connections for automated testing
closeAllConnections()

## End(Not run)

---