add ensemble and standata methods; more tests

Author: nicholasjclark

NAMESPACE

@@ -14,6 +14,7 @@ S3method(as_draws_matrix,mvgam)
 S3method(as_draws_rvars,mvgam)
 S3method(coef,mvgam)
 S3method(conditional_effects,mvgam)
+S3method(ensemble,mvgam_forecast)
 S3method(find_predictors,mvgam)
 S3method(find_predictors,mvgam_prefit)
 S3method(fitted,mvgam)
@@ -59,6 +60,7 @@ S3method(smooth.construct,mod.smooth.spec)
 S3method(smooth.construct,moi.smooth.spec)
 S3method(stancode,mvgam)
 S3method(stancode,mvgam_prefit)
+S3method(standata,mvgam_prefit)
 S3method(summary,mvgam)
 S3method(summary,mvgam_prefit)
 S3method(update,mvgam)
@@ -76,6 +78,7 @@ export(code)
 export(compare_mvgams)
 export(drawDotmvgam)
 export(dynamic)
+export(ensemble)
 export(eval_mvgam)
 export(eval_smoothDothilbertDotsmooth)
 export(eval_smoothDotmodDotsmooth)
@@ -135,6 +138,7 @@ importFrom(brms,qstudent_t)
 importFrom(brms,rbeta_binomial)
 importFrom(brms,rstudent_t)
 importFrom(brms,stancode)
+importFrom(brms,standata)
 importFrom(brms,student)
 importFrom(ggplot2,scale_colour_discrete)
 importFrom(ggplot2,scale_fill_discrete)

NEWS.md

@@ -1,8 +1,9 @@
 # mvgam 1.1.3 (development version; not yet on CRAN)
 ## New functionalities
-* Allow intercepts to be included in process models when `trend_formula` is supplied. This breaks the assumption that the process has to be zero-centred, adding flexibility but also potentially inducing nonidentifiabilities with respect to any observation model intercepts. Thoughtful priors are a must for these models
-* Added `stancode.mvgam` and `stancode.mvgam_prefit` methods
+* Allow intercepts to be included in process models when `trend_formula` is supplied. This breaks the assumption that the process has to be zero-centred, adding more modelling flexibility but also potentially inducing nonidentifiabilities with respect to any observation model intercepts. Thoughtful priors are a must for these models
+* Added `standata.mvgam_prefit`, `stancode.mvgam` and `stancode.mvgam_prefit` methods for better alignment with 'brms' workflows
 * Added 'gratia' to *Enhancements* to allow popular methods such as `draw()` to be used for 'mvgam' models if 'gratia' is already installed
+* Added an `ensemble.mvgam_forecast` method to generate evenly weighted combinations of probabilistic forecast distributions
 
 ## Deprecations
 * The `drift` argument has been deprecated. It is now recommended for users to include parametric fixed effects of "time" in their respective GAM formulae to capture any expected drift effects

R/ensemble.R

@@ -0,0 +1,175 @@
+#' Combine mvgam forecasts into evenly weighted ensembles
+#'
+#' Generate evenly weighted ensemble forecast distributions from \code{mvgam_forecast} objects
+#'
+#'@name ensemble.mvgam_forecast
+#'@param object \code{list} object of class \code{mvgam_forecast}. See [forecast.mvgam()]
+#'@param ... More \code{mvgam_forecast} objects.
+#'@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 \code{mvgam_forecast} object will already implicitly capture correlations among
+#' the temporal posterior predictions.
+#'@return An object of class \code{mvgam_forecast} containing the ensemble predictions. This
+#'object can be readily used with the supplied S3 functions \code{plot} and \code{score}
+#'@author Nicholas J Clark
+#'@seealso \code{\link{plot.mvgam_forecast}}, \code{\link{score.mvgam_forecast}}
+#' @examples
+#' \donttest{
+#' # 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)
+#'
+#' m2 <- mvgam(y ~ time,
+#'             trend_model = RW(),
+#'             noncentred = TRUE,
+#'             data = simdat$data_train,
+#'             newdata = simdat$data_test)
+#'
+#' # 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)
+#' }
+#'@export
+ensemble <- function(object, ...){
+  UseMethod("ensemble", object)
+}
+
+#'@rdname ensemble.mvgam_forecast
+#'@method ensemble mvgam_forecast
+#'@param 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
+#'@export
+ensemble.mvgam_forecast <- function(object, ..., ndraws = 5000){
+  models <- split_fc_dots(object, ..., model_names = NULL)
+  n_models <- length(models)
+
+  # Check that series names and key dimensions match for all forecasts
+  allsame <- function(x) length(unique(x)) == 1
+  if(!allsame(purrr::map(models, 'series_names'))){
+    stop('Names of series must match for all forecast objects.',
+         call. = FALSE)
+  }
+
+  if(!allsame(lapply(models, function(x) length(x$forecasts)))){
+    stop('The number of forecast distributions must match for all forecast objects.',
+         call. = FALSE)
+  }
+
+  if(!allsame(lapply(models, function(x) length(x$test_observations)))){
+    stop('Validation data must match for all forecast objects.',
+         call. = FALSE)
+  }
+
+  if(!allsame(lapply(models, function(x) {
+    unlist(lapply(x$forecasts, function(y) dim(y)[2]),
+           use.names = FALSE) }))){
+    stop('Forecast horizons must match for all forecast objects.',
+         call. = FALSE)
+  }
+
+  validate_pos_integer(ndraws)
+
+  # End of checks; now proceed with ensembling
+  n_series <- length(models[[1]]$series_names)
+
+  # Function to random sample rows of a matrix with
+  # replacement (in case some forecasts contain fewer draws than others)
+  subsamp <- function(x, nsamps){
+    if(NROW(x) < nsamps){
+      sampinds <- sample(1:NROW(x), nsamps, replace = TRUE)
+    } else {
+      sampinds <- 1:nsamps
+    }
+
+    x[sampinds, ]
+  }
+
+  # Create evenly weighted ensemble forecasts
+  ens_fcs <- lapply(seq_len(n_series), function(series){
+    all_fcs <- do.call(rbind,
+                       lapply(models,
+                              function(x) x$forecasts[[series]]))
+    subsamp(all_fcs, ndraws)
+  })
+
+  # Initiate the ensemble forecast object
+  ens_fc <- models[[1]]
+
+  # Add in forecasts
+  ens_fc$forecasts <- ens_fcs
+  names(ens_fc$forecasts) <- names(models[[1]]$forecasts)
+
+  # Ensure hindcasts have same number of samples
+  ens_hcs <- lapply(seq_len(n_series), function(series){
+    subsamp(ens_fc$hindcasts[[series]], ndraws)
+  })
+  ens_fc$hindcasts <- ens_hcs
+  names(ens_fc$hindcasts) <- names(models[[1]]$hindcasts)
+
+  # Return
+  return(ens_fc)
+}
+
+
+#'@noRd
+split_fc_dots = function (x, ..., model_names = NULL, other = TRUE) {
+
+  dots <- list(x, ...)
+  names <- substitute(list(x, ...), env = parent.frame())[-1]
+  names <- ulapply(names, deparse)
+
+  if(!is.null(model_names)){
+    names <- model_names
+  }
+
+  if (length(names)) {
+    if (!length(names(dots))) {
+      names(dots) <- names
+    }
+    else {
+      has_no_name <- !nzchar(names(dots))
+      names(dots)[has_no_name] <- names[has_no_name]
+    }
+  }
+  is_mvgam_fc <- unlist(lapply(dots, function(y) inherits(y, 'mvgam_forecast')))
+  models <- dots[is_mvgam_fc]
+  out <- dots[!is_mvgam_fc]
+
+  if (length(out)) {
+    stop("Only mvgam_forecast objects can be passed to '...' for this method.",
+         call. = FALSE)
+  }
+  models
+}
+
+

R/families.R

@@ -454,7 +454,7 @@ mvgam_predict = function(Xp,
       mu <- ((matrix(Xp, ncol = NCOL(Xp)) %*%
                 betas)) + attr(Xp, 'model.offset')
       sd <- as.vector(family_pars$sigma_obs)
-      out <- (exp((sd) ^ 2) - 1) * exp((2 * mu + sd ^ 2))
+      out <- as.vector((exp((sd) ^ 2) - 1) * exp((2 * mu + sd ^ 2)))
 
     } else {
       mu <- as.vector((matrix(Xp, ncol = NCOL(Xp)) %*%
@@ -610,10 +610,10 @@ mvgam_predict = function(Xp,
       out <- ((n * p) * (1 - p)) * ((alpha + beta + n) / (alpha + beta + 1))
 
     } else {
-      out <- plogis(((matrix(Xp, ncol = NCOL(Xp)) %*%
+      out <- as.vector(plogis(((matrix(Xp, ncol = NCOL(Xp)) %*%
                         betas)) +
                       attr(Xp, 'model.offset')) *
-        as.vector(family_pars$trials)
+        as.vector(family_pars$trials))
     }
   }
 
@@ -640,9 +640,9 @@ mvgam_predict = function(Xp,
                              betas) + attr(Xp, 'model.offset')))
       out <- mu + mu^2 / as.vector(family_pars$phi)
     } else {
-      out <- exp(((matrix(Xp, ncol = NCOL(Xp)) %*%
+      out <- as.vector(exp(((matrix(Xp, ncol = NCOL(Xp)) %*%
                      betas)) +
-                   attr(Xp, 'model.offset'))
+                   attr(Xp, 'model.offset')))
     }
   }
 
@@ -735,9 +735,9 @@ mvgam_predict = function(Xp,
                    attr(Xp, 'model.offset')) ^ 1.5) *
         as.vector(family_pars$phi)
       } else {
-      out <- exp(((matrix(Xp, ncol = NCOL(Xp)) %*%
+      out <- as.vector(exp(((matrix(Xp, ncol = NCOL(Xp)) %*%
                      betas)) +
-                   attr(Xp, 'model.offset'))
+                   attr(Xp, 'model.offset')))
     }
   }
   return(out)

R/forecast.mvgam.R

@@ -19,7 +19,7 @@
 #'@details Posterior predictions are drawn from the fitted \code{mvgam} and used to simulate a forecast distribution
 #'@return An object of class \code{mvgam_forecast} containing hindcast and forecast distributions.
 #'See \code{\link{mvgam_forecast-class}} for details.
-#'@seealso \code{\link{hindcast}}, \code{\link{score}}
+#'@seealso \code{\link{hindcast}}, \code{\link{score}}, \code{\link{ensemble}}
 #'@export
 forecast <- function(object, ...){
   UseMethod("forecast", object)

R/mvgam.R

@@ -366,7 +366,8 @@
 #'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;
+#'# 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,
@@ -377,9 +378,13 @@
 #'              run_model = FALSE)
 #'
 #' # View the model code in Stan language
-#' code(mod1)
+#' stancode(mod1)
 #'
-#' # Now fit the model, noting that 'noncentred = TRUE' will likely give performance gains
+#' # 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(),

R/score.mvgam_forecast.R

@@ -49,7 +49,7 @@
 #'str(fc_scores)
 #'}
 #'@method score mvgam_forecast
-#'@seealso \code{\link{forecast.mvgam}}
+#'@seealso \code{\link{forecast.mvgam}}, \code{\link{ensemble}}
 #'@export
 score.mvgam_forecast = function(object, score = 'crps',
                                 log = FALSE, weights,

R/stan_utils.R

@@ -1,12 +1,12 @@
-#' Stan code for mvgam models
+#' Stan code and data objects for mvgam models
 #'
-#' Generate Stan code for \pkg{mvgam} models
+#' Generate Stan code and data objects for \pkg{mvgam} models
 #'
 #' @param object An object of class `mvgam` or `mvgam_prefit`,
 #' returned from a call to \code{mvgam}
-#' @return A character string containing the fully commented \pkg{Stan} code
-#'   to fit a \pkg{mvgam} model. It is of class \code{c("character", "brmsmodel")}
-#'   to facilitate pretty printing.
+#' @return Either a character string containing the fully commented \pkg{Stan} code
+#'   to fit a \pkg{mvgam} model or a named list containing the data objects needed
+#'   to fit the model in Stan.
 #' @export
 #' @examples
 #' simdat <- sim_mvgam()
@@ -15,8 +15,14 @@
 #'              family = poisson(),
 #'              data = simdat$data_train,
 #'              run_model = FALSE)
+#'
+#' # View Stan model code
 #' stancode(mod)
 #'
+#' # View Stan model data
+#' sdata <- standata(mod)
+#' str(sdata)
+#'
 code = function(object){
   if(!class(object) %in% c('mvgam', 'mvgam_prefit')){
     stop('argument "object" must be of class "mvgam" or "mvgam_prefit"')
@@ -49,6 +55,15 @@ stancode.mvgam = function(object, ...){
   code(object)
 }
 
+#' @export
+#' @importFrom brms standata
+#' @param ... ignored
+#' @rdname code
+standata.mvgam_prefit = function(object, ...){
+
+  object$model_data
+}
+
 #' @noRd
 remove_likelihood = function(model_file){
   like_line <- grep('// likelihood functions',