Overview of the mvgam package

Dynamic GAMs

mvgam is designed to propagate unobserved temporal processes to capture latent dynamics in the observed time series. This works in a state-space format, with the temporal trend evolving independently of the observation process. An introduction to the package and some worked examples are also shown in this seminar: Ecological Forecasting with Dynamic Generalized Additive Models. Briefly, assume $\tilde{\boldsymbol{y}}_{i,t}$ is the conditional expectation of response variable i at time t. Assuming y_i is drawn from an exponential distribution with an invertible link function, the linear predictor for a multivariate Dynamic GAM can be written as:

for i in 1 : N_series ... for t in 1 : N_timepoints ...

$$g^{-1}(\tilde{\boldsymbol{y}}_{i,t})=\alpha_{i}+\sum\limits_{j=1}^J\boldsymbol{s}_{i,j,t}\boldsymbol{x}_{j,t}+\boldsymbol{Z}\boldsymbol{z}_{k,t}\,,$$ Here α are the unknown intercepts, the s’s are unknown smooth functions of covariates (x’s), which can potentially vary among the response series, and z are dynamic latent processes. Each smooth function s_j is composed of basis expansions whose coefficients, which must be estimated, control the functional relationship between x_j and $g^{-1}(\tilde{\boldsymbol{y}})$. The size of the basis expansion limits the smooth’s potential complexity. A larger set of basis functions allows greater flexibility. For more information on GAMs and how they can smooth through data, see this blogpost on how to interpret nonlinear effects from Generalized Additive Models. Latent processes are captured with Zz_i, t, where Z is an i by k matrix of loading coefficients (which can be fixed or a combination of fixed and freely estimated parameters) and z_k, t are a set of K latent factors that can also include their own GAM linear predictors (see the State-Space models vignette), the N-mixtures vignette and the example in jsdgam to get an idea of how flexible these processes can be.

Several advantages of GAMs are that they can model a diversity of response families, including discrete distributions (i.e. Poisson, Negative Binomial, Gamma) that accommodate common ecological features such as zero-inflation or overdispersion, and that they can be formulated to include hierarchical smoothing for multivariate responses. mvgam supports a number of different observation families, which are summarized below:

Supported observation families

For all supported observation families, any extra parameters that need to be estimated (i.e. the σ in a Gaussian model or the ϕ in a Negative Binomial model) are by default estimated independently for each series. However, users can opt to force all series to share extra observation parameters using share_obs_params = TRUE in mvgam(). Note that default link functions cannot currently be changed.

Supported temporal dynamic processes

As stated above, the latent processes can take a wide variety of forms, some of which can be multivariate to allow the different observational variables to interact or be correlated. When using the mvgam() function, the user chooses between different process models with the trend_model argument. Available process models are described in detail below.

Regression formulae

mvgam supports an observation model regression formula, built off the mgcv package, as well as an optional process model regression formula. The formulae supplied to are exactly like those supplied to glm() except that smooth terms, s(), te(), ti() and t2(), time-varying effects using dynamic(), monotonically increasing (using s(x, bs = 'moi')) or decreasing splines (using s(x, bs = 'mod'); see ?smooth.construct.moi.smooth.spec for details), as well as Gaussian Process functions using gp(), can be added to the right hand side (and . is not supported in mvgam formulae). See ?mvgam_formulae for more guidance.

For setting up State-Space models, the optional process model formula can be used (see the State-Space model vignette and the shared latent states vignette for guidance on using trend formulae).

Example time series data

The ‘portal_data’ object contains time series of rodent captures from the Portal Project, a long-term monitoring study based near the town of Portal, Arizona. Researchers have been operating a standardized set of baited traps within 24 experimental plots at this site since the 1970’s. Sampling follows the lunar monthly cycle, with observations occurring on average about 28 days apart. However, missing observations do occur due to difficulties accessing the site (weather events, COVID disruptions etc…). You can read about the full sampling protocol in this preprint by Ernest et al on the Biorxiv.

data("portal_data")

As the data come pre-loaded with the mvgam package, you can read a little about it in the help page using ?portal_data. Before working with data, it is important to inspect how the data are structured, first using head:

head(portal_data)

But the glimpse function in dplyr is also useful for understanding how variables are structured

dplyr::glimpse(portal_data)

We will focus analyses on the time series of captures for one specific rodent species, the Desert Pocket Mouse Chaetodipus penicillatus. This species is interesting in that it goes into a kind of “hibernation” during the colder months, leading to very low captures during the winter period

Manipulating data for modelling

Manipulating the data into a ‘long’ format is necessary for modelling in mvgam. By ‘long’ format, we mean that each series x time observation needs to have its own entry in the dataframe or list object that we wish to use as data for modelling. A simple example can be viewed by simulating data using the sim_mvgam function. See ?sim_mvgam for more details

data <- sim_mvgam(n_series = 4, T = 24)
head(data$data_train, 12)

Notice how we have four different time series in these simulated data, but we do not spread the outcome values into different columns. Rather, there is only a single column for the outcome variable, labelled y in these simulated data. We also must supply a variable labelled time to ensure the modelling software knows how to arrange the time series when building models. This setup still allows us to formulate multivariate time series models, as you can see in the State-Space vignette. Below are the steps needed to shape our portal_data object into the correct form. First, we create a time variable, select the column representing counts of our target species (PP), and select appropriate variables that we can use as predictors

portal_data %>%
  
  # mvgam requires a 'time' variable be present in the data to index
  # the temporal observations. This is especially important when tracking 
  # multiple time series. In the Portal data, the 'moon' variable indexes the
  # lunar monthly timestep of the trapping sessions
  dplyr::mutate(time = moon - (min(moon)) + 1) %>%
  
  # We can also provide a more informative name for the outcome variable, which 
  # is counts of the 'PP' species (Chaetodipus penicillatus) across all control
  # plots
  dplyr::mutate(count = PP) %>%
  
  # The other requirement for mvgam is a 'series' variable, which needs to be a
  # factor variable to index which time series each row in the data belongs to.
  # Again, this is more useful when you have multiple time series in the data
  dplyr::mutate(series = as.factor('PP')) %>%
  
  # Select the variables of interest to keep in the model_data
  dplyr::select(series, year, time, count, mintemp, ndvi) -> model_data

The data now contain six variables:
series, a factor indexing which time series each observation belongs to
year, the year of sampling
time, the indicator of which time step each observation belongs to
count, the response variable representing the number of captures of the species PP in each sampling observation
mintemp, the monthly average minimum temperature at each time step
ndvi, the monthly average Normalized Difference Vegetation Index at each time step

Now check the data structure again

head(model_data)

dplyr::glimpse(model_data)

You can also summarize multiple variables, which is helpful to search for data ranges and identify missing values

summary(model_data)

We have some NAs in our response variable count. These observations will generally be thrown out by most modelling packages in . But as you will see when we work through the tutorials, mvgam keeps these in the data so that predictions can be automatically returned for the full dataset. The time series and some of its descriptive features can be plotted using plot_mvgam_series():

plot_mvgam_series(data = model_data, series = 1, y = 'count')

GLMs with temporal random effects

Our first task will be to fit a Generalized Linear Model (GLM) that can adequately capture the features of our count observations (integer data, lower bound at zero, missing values) while also attempting to model temporal variation. We are almost ready to fit our first model, which will be a GLM with Poisson observations, a log link function and random (hierarchical) intercepts for year. This will allow us to capture our prior belief that, although each year is unique, having been sampled from the same population of effects, all years are connected and thus might contain valuable information about one another. This will be done by capitalizing on the partial pooling properties of hierarchical models. Hierarchical (also known as random) effects offer many advantages when modelling data with grouping structures (i.e. multiple species, locations, years etc…). The ability to incorporate these in time series models is a huge advantage over traditional models such as ARIMA or Exponential Smoothing. But before we fit the model, we will need to convert year to a factor so that we can use a random effect basis in mvgam. See ?smooth.terms and ?smooth.construct.re.smooth.spec for details about the re basis construction that is used by both mvgam and mgcv

model_data %>%
  
  # Create a 'year_fac' factor version of 'year'
  dplyr::mutate(year_fac = factor(year)) -> model_data

Preview the dataset to ensure year is now a factor with a unique factor level for each year in the data

dplyr::glimpse(model_data)
levels(model_data$year_fac)

We are now ready for our first mvgam model. The syntax will be familiar to users who have previously built models with mgcv. But for a refresher, see ?formula.gam and the examples in ?gam. Random effects can be specified using the s wrapper with the re basis. Note that we can also suppress the primary intercept using the usual R formula syntax - 1. mvgam has a number of possible observation families that can be used, see ?mvgam_families for more information. We will use Stan as the fitting engine, which deploys Hamiltonian Monte Carlo (HMC) for full Bayesian inference. By default, 4 HMC chains will be run using a warmup of 500 iterations and collecting 500 posterior samples from each chain. The package will also aim to use the Cmdstan backend when possible, so it is recommended that users have an up-to-date installation of Cmdstan and the associated cmdstanr interface on their machines (note that you can set the backend yourself using the backend argument: see ?mvgam for details). Interested users should consult the Stan user’s guide for more information about the software and the enormous variety of models that can be tackled with HMC.

model1 <- mvgam(count ~ s(year_fac, bs = 're') - 1,
                family = poisson(),
                data = model_data)

The model can be described mathematically for each timepoint t as follows:

Where the β_year effects are drawn from a population distribution that is parameterized by a common mean (μ_year) and variance (σ_year). Priors on most of the model parameters can be interrogated and changed using similar functionality to the options available in brms. For example, the default priors on (μ_year) and (σ_year) can be viewed using the following code:

get_mvgam_priors(count ~ s(year_fac, bs = 're') - 1,
                 family = poisson(),
                 data = model_data)

See examples in ?get_mvgam_priors to find out different ways that priors can be altered. Once the model has finished, the first step is to inspect the summary to ensure no major diagnostic warnings have been produced and to quickly summarise posterior distributions for key parameters

summary(model1)

The diagnostic messages at the bottom of the summary show that the HMC sampler did not encounter any problems or difficult posterior spaces. This is a good sign. Posterior distributions for model parameters can be extracted in any way that an object of class brmsfit can (see ?mvgam::mvgam_draws for details). For example, we can extract the coefficients related to the GAM linear predictor (i.e. the β’s) into a data.frame using:

beta_post <- as.data.frame(model1, variable = 'betas')
dplyr::glimpse(beta_post)

With any model fitted in mvgam, the underlying Stan code can be viewed using the code function:

code(model1)

plot(model1, type = 're')

mcmc_plot(object = model1,
          variable = 'betas',
          type = 'areas')

pp_check(object = model1)

plot(model1, type = 'forecast')

hc <- hindcast(model1)
str(hc)

hc <- hindcast(model1, type = 'link')
range(hc$hindcasts$PP)

plot(model1, type = 'residuals')

Automatic forecasting for new data

These temporal random effects do not have a sense of “time”. Because of this, each yearly random intercept is not restricted in some way to be similar to the previous yearly intercept. This drawback becomes evident when we predict for a new year. To do this, we can repeat the exercise above but this time will split the data into training and testing sets before re-running the model. We can then supply the test set as newdata. For splitting, we will make use of the filter function from dplyr

model_data %>% 
  dplyr::filter(time <= 160) -> data_train 
model_data %>% 
  dplyr::filter(time > 160) -> data_test

model1b <- mvgam(count ~ s(year_fac, bs = 're') - 1,
                 family = poisson(),
                 data = data_train,
                 newdata = data_test)

We can view the test data in the forecast plot to see that the predictions do not capture the temporal variation in the test set

plot(model1b, type = 'forecast', newdata = data_test)

As with the hindcast function, we can use the forecast function to automatically extract the posterior distributions for these predictions. This also returns an object of class mvgam_forecast, but now it will contain both the hindcasts and forecasts for each series in the data:

fc <- forecast(model1b)
str(fc)

Adding predictors as “fixed” effects

Any users familiar with GLMs will know that we nearly always wish to include predictor variables that may explain some of the variation in our observations. Predictors are easily incorporated into GLMs / GAMs. Here, we will update the model from above by including a parametric (fixed) effect of ndvi as a linear predictor:

model2 <- mvgam(count ~ s(year_fac, bs = 're') + 
                  ndvi - 1,
                family = poisson(),
                data = data_train,
                newdata = data_test)

The model can be described mathematically as follows:

Where the β_year effects are the same as before but we now have another predictor (β_ndvi) that applies to the ndvi value at each timepoint t. Inspect the summary of this model

summary(model2)

Rather than printing the summary each time, we can also quickly look at the posterior empirical quantiles for the fixed effect of ndvi (and other linear predictor coefficients) using coef:

coef(model2)

Look at the estimated effect of ndvi using using a histogram. This can be done by first extracting the posterior coefficients:

beta_post <- as.data.frame(model2, variable = 'betas')
dplyr::glimpse(beta_post)

The posterior distribution for the effect of ndvi is stored in the ndvi column. A quick histogram confirms our inference that log(counts) respond positively to increases in ndvi:

hist(beta_post$ndvi,
     xlim = c(-1 * max(abs(beta_post$ndvi)),
              max(abs(beta_post$ndvi))),
     col = 'darkred',
     border = 'white',
     xlab = expression(beta[NDVI]),
     ylab = '',
     yaxt = 'n',
     main = '',
     lwd = 2)
abline(v = 0, lwd = 2.5)

conditional_effects(model2)

Adding predictors as smooths

Smooth functions, using penalized splines, are a major feature of mvgam. Nonlinear splines are commonly viewed as variations of random effects in which the coefficients that control the shape of the spline are drawn from a joint, penalized distribution. This strategy is very often used in ecological time series analysis to capture smooth temporal variation in the processes we seek to study. When we construct smoothing splines, the workhorse package mgcv will calculate a set of basis functions that will collectively control the shape and complexity of the resulting spline. It is often helpful to visualize these basis functions to get a better sense of how splines work. We’ll create a set of 6 basis functions to represent possible variation in the effect of time on our outcome.In addition to constructing the basis functions, mgcv also creates a penalty matrix S, which contains known coefficients that work to constrain the wiggliness of the resulting smooth function. When fitting a GAM to data, we must estimate the smoothing parameters (λ) that will penalize these matrices, resulting in constrained basis coefficients and smoother functions that are less likely to overfit the data. This is the key to fitting GAMs in a Bayesian framework, as we can jointly estimate the λ’s using informative priors to prevent overfitting and expand the complexity of models we can tackle. To see this in practice, we can now fit a model that replaces the yearly random effects with a smooth function of time. We will need a reasonably complex function (large k) to try and accommodate the temporal variation in our observations. Following some useful advice by Gavin Simpson, we will use a b-spline basis for the temporal smooth. Because we no longer have intercepts for each year, we also retain the primary intercept term in this model (there is no -1 in the formula now):

model3 <- mvgam(count ~ s(time, bs = 'bs', k = 15) + 
                  ndvi,
                family = poisson(),
                data = data_train,
                newdata = data_test)

The model can be described mathematically as follows:

Where the smooth function f_time is built by summing across a set of weighted basis functions. The basis functions (b) are constructed using a thin plate regression basis in mgcv. The weights (β_smooth) are drawn from a penalized multivariate normal distribution where the precision matrix (Ω) is multiplied by a smoothing penalty (λ). If λ becomes large, this acts to squeeze the covariances among the weights (β_smooth), leading to a less wiggly spline. Note that sometimes there are multiple smoothing penalties that contribute to the covariance matrix, but I am only showing one here for simplicity. View the summary as before

summary(model3)

The summary above now contains posterior estimates for the smoothing parameters as well as the basis coefficients for the nonlinear effect of time. We can visualize conditional_effects as before:

conditional_effects(model3, type = 'link')

Inspect the underlying Stan code to gain some idea of how the spline is being penalized:

code(model3)

The line below // prior for s(time)... shows how the spline basis coefficients are drawn from a zero-centred multivariate normal distribution. The precision matrix S is penalized by two different smoothing parameters (the λ’s) to enforce smoothness and reduce overfitting

Latent dynamics in mvgam

Forecasts from the above model are not ideal:

plot(model3, type = 'forecast', newdata = data_test)

Why is this happening? The forecasts are driven almost entirely by variation in the temporal spline, which is extrapolating linearly forever beyond the edge of the training data. Any slight wiggles near the end of the training set will result in wildly different forecasts. To visualize this, we can plot the extrapolated temporal functions into the out-of-sample test set for the two models. Here are the extrapolated functions for the first model, with 15 basis functions:

plot_mvgam_smooth(model3, smooth = 's(time)',
                  # feed newdata to the plot function to generate
                  # predictions of the temporal smooth to the end of the 
                  # testing period
                  newdata = data.frame(time = 1:max(data_test$time),
                                       ndvi = 0))
abline(v = max(data_train$time), lty = 'dashed', lwd = 2)

This model is not doing well. Clearly we need to somehow account for the strong temporal autocorrelation when modelling these data without using a smooth function of time. Now onto another prominent feature of mvgam: the ability to include (possibly latent) autocorrelated residuals in regression models. To do so, we use the trend_model argument (see ?mvgam_trends for details of different dynamic trend models that are supported). This model will use a separate sub-model for latent residuals that evolve as an AR1 process (i.e. the error in the current time point is a function of the error in the previous time point, plus some stochastic noise). We also include a smooth function of ndvi in this model, rather than the parametric term that was used above, to showcase that mvgam can include combinations of smooths and dynamic components:

model4 <- mvgam(count ~ s(ndvi, k = 6),
                family = poisson(),
                data = data_train,
                newdata = data_test,
                trend_model = 'AR1')

The model can be described mathematically as follows:

Here the term z_t captures autocorrelated latent residuals, which are modelled using an AR1 process. You can also notice that this model is estimating autocorrelated errors for the full time period, even though some of these time points have missing observations. This is useful for getting more realistic estimates of the residual autocorrelation parameters. Summarise the model to see how it now returns posterior summaries for the latent AR1 process:

summary(model4)

View posterior hindcasts / forecasts and compare against the out of sample test data

plot(model4, type = 'forecast', newdata = data_test)

The trend is evolving as an AR1 process, which we can also view:

plot(model4, type = 'trend', newdata = data_test)

In-sample model performance can be interrogated using leave-one-out cross-validation utilities from the loo package (a higher value is preferred for this metric):

loo_compare(model3, model4)

The higher estimated log predictive density (ELPD) value for the dynamic model suggests it provides a better fit to the in-sample data.

Though it should be obvious that this model provides better forecasts, we can quantify forecast performance for models 3 and 4 using the forecast and score functions. Here we will compare models based on their Discrete Ranked Probability Scores (a lower value is preferred for this metric)

fc_mod3 <- forecast(model3)
fc_mod4 <- forecast(model4)
score_mod3 <- score(fc_mod3, score = 'drps')
score_mod4 <- score(fc_mod4, score = 'drps')
sum(score_mod4$PP$score, na.rm = TRUE) - sum(score_mod3$PP$score, na.rm = TRUE)

A strongly negative value here suggests the score for the dynamic model (model 4) is much smaller than the score for the model with a smooth function of time (model 3)

Interested in contributing?

I’m actively seeking PhD students and other researchers to work in the areas of ecological forecasting, multivariate model evaluation and development of mvgam. Please reach out if you are interested (n.clark’at’uq.edu.au)