diff --git a/DESCRIPTION b/DESCRIPTION index a4cb41de..5cabb8ff 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mvgam Title: Multivariate Bayesian Generalised Additive Models for discrete time series -Version: 0.0.0.9000 +Version: 1.0.0 Authors@R: person(given = "Nicholas", family = "Clark", @@ -9,8 +9,8 @@ Authors@R: Description: This package is for fitting Bayesian Generalised Additive Models to sets of discrete time series. The primary purpose of the package is to build a version of dynamic factor models that incorporates the flexibility of GAMs in the linear predictor and latent dynamic trend components that are useful - for time series analysis and forecasting. Estimation is performed using Markov Chain Monte Carlo (MCMC) by - the Gibbs sampling software JAGS. The package also includes + for time series analysis and forecasting. Estimation is performed using Markov Chain Monte Carlo (MCMC) with + either the Gibbs sampling software JAGS or with Hamiltonian Monte Carlo in the software Stan. The package also includes utilities for online updating of forecast distributions with a recursive particle filter that uses sequential importance sampling to assimilate new observations as they become available. Maintainer: Nicholas J Clark diff --git a/NAMESPACE b/NAMESPACE index 624835f9..ab37485c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,5 @@ # Generated by roxygen2: do not edit by hand -S3method(dic,mvgam) S3method(plot,mvgam) S3method(ppc,mvgam) S3method(predict,mvgam) @@ -15,11 +14,10 @@ export(add_stan_data) export(add_trend_lines) export(add_tweedie_lines) export(compare_mvgams) -export(dic) export(eval_mvgam) export(get_monitor_pars) export(lv_correlations) -export(mvjagam) +export(mvgam) export(pfilter_mvgam_fc) export(pfilter_mvgam_init) export(pfilter_mvgam_online) diff --git a/NEON_manuscript/Case studies/mvgam_case_study1.Rmd b/NEON_manuscript/Case studies/mvgam_case_study1.Rmd index 4ded03d4..bd26685f 100644 --- a/NEON_manuscript/Case studies/mvgam_case_study1.Rmd +++ b/NEON_manuscript/Case studies/mvgam_case_study1.Rmd @@ -93,9 +93,9 @@ head(fake_data$data_train) head(fake_data$data_test) ``` -As a first pass at modelling this series, we will fit a GAM that includes a seasonal smooth and a yearly smooth, as well as their tensor product interaction. As the seasonality is cyclic, we will use a cyclic cubic regression spline for `season`. The knots are set so that the boundaries of the cyclic smooth match up between December 31 and January 1. We will stick with the default thin plate regression spline for `year`. This is similar to what we might do when fitting a model in `mgcv` to try and forecast ahead, except here we also have an explicit model for the residual component. `mvgam` uses the `jagam` function from the `mgcv` package to generate a skeleton `JAGS` file and updates that file to incorporate any dynamic trend components (so far limited to no trend, random walks or AR models up to order 3). This is advantageous as any GAM that is allowed in `mgcv` can in principle be used in `mvgam` to fit dynamic linear models with smooth functions for nonlinear covariate effects. For multivariate series, `mvgam` also includes an option for dimension reduction by inducing trend correlations via a dynamic factor process. Here we use the Negative Binomial family and a burnin length of `10000` iterations (in practice we would probably run these models for a bit longer, but they tend to converge rather quickly for univariate series thanks to the useful starting values supplied by `jagam`). Note that feeding the `data_test` object does not mean that these values are used in training of the model. Rather, they are included as `NA` so that we can automatically create a forecast from the posterior predictions for these observations. This is useful for plotting forecasts without needing to run new data through the model's equations later on +As a first pass at modelling this series, we will fit a GAM that includes a seasonal smooth and a yearly smooth, as well as their tensor product interaction. As the seasonality is cyclic, we will use a cyclic cubic regression spline for `season`. The knots are set so that the boundaries of the cyclic smooth match up between December 31 and January 1. We will stick with the default thin plate regression spline for `year`. This is similar to what we might do when fitting a model in `mgcv` to try and forecast ahead, except here we also have an explicit model for the residual component. `mvgam` uses the `jagam` function from the `mgcv` package to generate a skeleton `JAGS` file and updates that file to incorporate any dynamic trend components (so far limited to no trend, random walks or AR models up to order 3 for `JAGS` models). This is advantageous as any GAM that is allowed in `mgcv` can in principle be used in `mvgam` to fit dynamic linear models with smooth functions for nonlinear covariate effects. For multivariate series, `mvgam` also includes an option for dimension reduction by inducing trend correlations via a dynamic factor process. Here we use the Negative Binomial family and a burnin length of `10000` iterations (in practice we would probably run these models for a bit longer, but they tend to converge rather quickly for univariate series thanks to the useful starting values supplied by `jagam`). Note that feeding the `data_test` object does not mean that these values are used in training of the model. Rather, they are included as `NA` so that we can automatically create a forecast from the posterior predictions for these observations. This is useful for plotting forecasts without needing to run new data through the model's equations later on ```{r, message=FALSE, warning=FALSE} -mod1 <- mvjagam(data_train = fake_data$data_train, +mod1 <- mvgam(data_train = fake_data$data_train, data_test = fake_data$data_test, formula = y ~ s(season, bs = c('cc'), k = 12) + s(year, k = 5) + @@ -107,7 +107,7 @@ mod1 <- mvjagam(data_train = fake_data$data_train, chains = 4) ``` -We can view the `JAGS` model file that has been generated to see the additions that have been made to the base `gam` model. If the user selects `return_model_data = TRUE` when calling `mvjagam`, this file can be modified and the resulting `model_data` object can also be modified to fit more complex Bayesian models. Note that here the AR and phi terms have been set to zero as this model does not include a dynamic trend component +We can view the `JAGS` model file that has been generated to see the additions that have been made to the base `gam` model. If the user selects `return_model_data = TRUE` when calling `mvgam`, this file can be modified and the resulting `model_data` object can also be modified to fit more complex Bayesian models. Note that here the AR and phi terms have been set to zero as this model does not include a dynamic trend component ```{r} mod1$model_file ``` @@ -201,7 +201,7 @@ ppc(mod1, series = 1, type = 'pit', data_test = fake_data$data_test) There are some problems with the way this model is generating future predictions. Perhaps a different smooth function for `year` can help? Here we fit our original model again but use a different smooth term for `year` to try and capture the long-term trend (using `B` splines with multiple penalties, following the excellent example by Gavin Simpson about [extrapolating with smooth terms](https://fromthebottomoftheheap.net/2020/06/03/extrapolating-with-gams/)). This is similar to what we might do when trying to forecast ahead from a more wiggly function, as `B` splines have useful properties by allowing the penalty to be extended into the range of values we wish to predict (in this case, the years in `data_test`). Note that we drop the `year*season` tensor product interaction to help emphasize the extrapolation behaviour of the `B` spline ```{r, message=FALSE, warning=FALSE} -mod2 <- mvjagam(data_train = fake_data$data_train, +mod2 <- mvgam(data_train = fake_data$data_train, data_test = fake_data$data_test, formula = y ~ s(season, bs = c('cc'), k = 12) + s(year, bs = "bs", m = c(2, 1)), @@ -259,7 +259,7 @@ ppc(mod2, series = 1, type = 'pit', data_test = fake_data$data_test) Can we improve the forecasts by removing our reliance on extrapolation? Now we will fit a model in which the GAM component of the linear predictor captures the repeated seasonality (again with a cyclic smooth) and a dynamic latent trend captures the temporal residual process using AR parameters (up to order 3). This model is a mildly modified version of the base `mgcv` model where the linear predictor is augmented with the latent trend component. Slightly longer burnin is used here due to the added complexity of the time series component, but the model still fits in ~ 30 seconds on most machines ```{r, message=FALSE, warning=FALSE} -mod3 <- mvjagam(data_train = fake_data$data_train, +mod3 <- mvgam(data_train = fake_data$data_train, data_test = fake_data$data_test, formula = y ~ s(season, bs = c('cc'), k = 12), knots = list(season = c(0.5, 12.5)), @@ -271,7 +271,7 @@ mod3 <- mvjagam(data_train = fake_data$data_train, chains = 4) ``` -In this case the fitted model is more different to the base `mgcv` model that was used in `jagam` to produce the skeleton `JAGS` file, so the summary of that base model is less accurate. But we can still check the model summary for the `mvjagam` mdoel to examine convergence for key parameters +In this case the fitted model is more different to the base `mgcv` model that was used in `jagam` to produce the skeleton `JAGS` file, so the summary of that base model is less accurate. But we can still check the model summary for the `mvgam` mdoel to examine convergence for key parameters ```{r} summary(mod3) ``` @@ -316,7 +316,7 @@ Benchmarking against "null" models is a very important part of evaluating a prop ```{r, message=FALSE, warning=FALSE} fake_data$data_train$fake_cov <- rnorm(NROW(fake_data$data_train)) fake_data$data_test$fake_cov <- rnorm(NROW(fake_data$data_test)) -mod4 <- mvjagam(data_train = fake_data$data_train, +mod4 <- mvgam(data_train = fake_data$data_train, data_test = fake_data$data_test, formula = y ~ s(fake_cov, k = 3), family = 'poisson', @@ -607,13 +607,13 @@ lynx_mgcv2_drps <- drps_mcmc_object(truth = lynx_test$population[1:10], This forecast is highly overconfident, with very unrealistic uncertainty intervals due to the interpolation behaviours of the lag smooths. You can certainly keep trying different formulations (our experience is that the B spline variant above produces the best forecasts from any tested `mgcv` model, but we did not test an exhaustive set), but hopefully it is clear that forecasting using splines is tricky business and it is likely that each time you do it you'll end up honing in on different combinations of penalties, knot selections etc.... Now we will fit an `mvgam` model for comparison. This model fits a similar model to the `mgcv` model directly above but with a full time series model for the errors (in this case an `AR1` process), rather than smoothing splines that do not incorporate a concept of the future. We do not use a `year` term to reduce any possible extrapolation and because the latent dynamic component should capture this temporal variation. ```{r, message=FALSE, warning=FALSE} -lynx_mvgam <- mvjagam(data_train = lynx_train, +lynx_mvgam <- mvgam(data_train = lynx_train, data_test = lynx_test, formula = y ~ s(season, bs = 'cc', k = 19), knots = list(season = c(0.5, 19.5)), family = 'poisson', trend_model = 'AR1', - burnin = 5000, + burnin = 8000, chains = 4) ``` diff --git a/NEON_manuscript/Case studies/mvgam_case_study1.html b/NEON_manuscript/Case studies/mvgam_case_study1.html index d71602ad..cea419a2 100644 --- a/NEON_manuscript/Case studies/mvgam_case_study1.html +++ b/NEON_manuscript/Case studies/mvgam_case_study1.html @@ -14,12 +14,49 @@ mvgam case study 1: model comparison and data assimilation - + + - - - - + + + + - - - - - + + + + + + + + - - - - - + + + + + + + + - - - - - + + + + +