| mvgam {mvgam} | R Documentation |
Fit a Bayesian dynamic GAM to a univariate or multivariate set of 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,
data,
data_train,
newdata,
data_test,
run_model = TRUE,
prior_simulation = FALSE,
return_model_data = FALSE,
family = "poisson",
share_obs_params = FALSE,
use_lv = FALSE,
n_lv,
trend_map,
trend_model = "None",
drift = FALSE,
noncentred = FALSE,
chains = 4,
burnin = 500,
samples = 500,
thin = 1,
parallel = TRUE,
threads = 1,
priors,
refit = FALSE,
lfo = FALSE,
residuals = TRUE,
use_stan = TRUE,
backend = getOption("brms.backend", "cmdstanr"),
algorithm = getOption("brms.algorithm", "sampling"),
autoformat = TRUE,
save_all_pars = FALSE,
max_treedepth = 12,
adapt_delta = 0.85,
silent = 1,
jags_path,
...
)
Arguments
formula |
A |
trend_formula |
An optional |
knots |
An optional |
trend_knots |
As for |
data |
A
Should also include any other variables to be included in the linear predictor of |
data_train |
Deprecated. Still works in place of |
newdata |
Optional |
data_test |
Deprecated. Still works in place of |
run_model |
|
prior_simulation |
|
return_model_data |
|
family |
Note that only |
share_obs_params |
|
use_lv |
|
n_lv |
|
trend_map |
Optional |
trend_model |
For all trend types apart from |
drift |
|
noncentred |
|
chains |
|
burnin |
|
samples |
|
thin |
Thinning interval for monitors. Ignored
if |
parallel |
|
threads |
|
priors |
An optional |
refit |
Logical indicating whether this is a refit, called using update.mvgam. Users should leave
as |
lfo |
Logical indicating whether this is part of a call to lfo_cv.mvgam. Returns a
lighter version of the model with no residuals and fewer monitored parameters to speed up
post-processing. But other downstream functions will not work properly, so users should always
leave this set as |
residuals |
Logical indicating whether to compute series-level randomized quantile residuals and include
them as part of the returned object. Defaults to |
use_stan |
Logical. If |
backend |
Character string naming the package to use as the backend for fitting
the Stan model (if |
algorithm |
Character string naming the estimation approach to use.
Options are |
autoformat |
|
save_all_pars |
|
max_treedepth |
positive integer placing a cap on the number of simulation steps evaluated during each iteration when
|
adapt_delta |
positive numeric between |
silent |
Verbosity level between |
jags_path |
Optional character vector specifying the path to the location of the |
... |
Further arguments passed to Stan.
For |
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. A general overview is provided
in the primary vignettes: vignette("mvgam_overview") and vignette("data_in_mvgam").
For a full list of available vignettes see vignette(package = "mvgam")
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, see examples below)
Families and link functions: Details of families supported by mvgam can be found in
mvgam_families.
Trend models: Details of latent trend dynamic models supported by mvgam can be found in
mvgam_trends.
Priors: A jagam model file is generated from formula and
modified to include any latent
temporal processes. Default priors for intercepts and any scale parameters are generated
using the same practice as brms. Prior distributions for most important model parameters can be
altered by the user to inspect model
sensitivities to given priors (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 first using mvgam as a
baseline, then editing the returned model accordingly. The model file can be edited and run outside
of mvgam by setting run_model = FALSE and this is encouraged for complex
modelling tasks. Note, no priors are
formally checked to ensure they are in the right syntax for the respective
probabilistic modelling framework, so it is
up to the user to ensure these are correct (i.e. use dnorm for normal
densities in JAGS, with the mean and precision
parameterisation; but use normal for normal densities in Stan, with the mean
and standard deviation parameterisation)
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
Factor regularisation: When using a dynamic factor model for the trends with JAGS factor precisions are given
regularized penalty priors to theoretically allow some factors to be dropped from the model by squeezing increasing
factors' variances to zero. This is done to help protect against selecting too many latent factors than are needed to
capture dependencies in the data, so it can often be advantageous to set n_lv to a slightly larger number. However
larger numbers of factors do come with additional computational costs so these should be balanced as well. When using
Stan, all factors are parameterised with fixed variance parameters
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
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 smooth latent trends via Hilbert space approximate Gaussian Processes.
This often makes sense for ecological series, which we expect to change smoothly. In mvgam, latent squared
exponential GP trends are approximated using by default 20 basis functions, which saves computational
costs compared to fitting full GPs while adequately estimating
GP alpha and rho parameters. 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
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.
Author(s)
Nicholas J Clark
References
Nicholas J Clark & Konstans Wells (2020). Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series. Methods in Ecology and Evolution. 14:3, 771-784.
See Also
jagam, gam, gam.models,
Examples
# Simulate a collection of three time series that have shared seasonal dynamics
# and independent AR1 trends, with a Poisson observation process
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;
# 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,
use_stan = TRUE,
run_model = FALSE)
# View the model code in Stan language
code(mod1)
# Now fit the model, noting that 'noncentred = TRUE' will likely give performance gains
mod1 <- mvgam(formula = y ~ s(season, bs = 'cc', k = 6),
data = dat$data_train,
trend_model = AR(),
family = poisson(),
noncentred = TRUE,
chains = 2)
# Extract the model summary
summary(mod1)
# Plot the estimated historical trend and forecast for one series
plot(mod1, type = 'trend', series = 1)
plot(mod1, type = 'forecast', series = 1)
# Residual diagnostics
plot(mod1, type = 'residuals', series = 1)
resids <- residuals(mod1)
str(resids)
# Compute the forecast using covariate information in data_test
fc <- forecast(mod1, newdata = dat$data_test)
str(fc)
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
library(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)
# 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 it's 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(y ~ s(season, bs = 'cc', k = 6),
trend_map = trend_map,
trend_model = AR(),
data = mod_data,
return_model_data = TRUE,
chains = 2)
# The mapping matrix is now supplied as data to the model in the 'Z' element
mod$model_data$Z
code(mod)
# The first two series share an identical latent trend; the third is different
plot(mod, type = 'trend', series = 1)
plot(mod, type = 'trend', series = 2)
plot(mod, type = 'trend', series = 3)
# 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 the Gaussian observation process
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() formula helper
mod <- mvgam(out ~
dynamic(temp,
scale = FALSE,
k = 40),
family = gaussian(),
data = data_train,
newdata = data_test,
chains = 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)
# 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)
# Inspect the model file to see the modification to the linear predictor
# (eta)
code(mod)
# Forecasts for the first two series will differ in magnitude
fc <- forecast(mod, newdata = dat$data_test)
layout(matrix(1:2, ncol = 2))
plot(fc, series = 1, ylim = c(0, 75))
plot(fc, series = 2, ylim = c(0, 75))
layout(1)
# 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)
# 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(cbind(y, ntrials) ~ series + s(x, by = series),
family = binomial(),
data = dat,
chains = 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')