update ppc docs

Author: nicholasjclark

R/forecast.mvgam.R

@@ -574,7 +574,7 @@ forecast.mvgam = function(object,
                                use_lv = object$use_lv,
                                fit_engine = object$fit_engine,
                                type = type,
-                               series_names = factor(unique(data_train$series),
+                               series_names = factor(levels(data_train$series),
                                                      levels = levels(data_train$series)),
                                train_observations = series_obs,
                                train_times = unique(data_train$index..time..index),
@@ -654,32 +654,35 @@ forecast_draws = function(object,
     Xp_trend <- trend_Xp_matrix(newdata = sort_data(data_test),
                                 trend_map = object$trend_map,
                                 series = series,
-                                mgcv_model = object$trend_mgcv_model)
+                                mgcv_model = object$trend_mgcv_model,
+                                forecast = TRUE)
 
     # For trend_formula models with autoregressive processes,
     # the process model operates as: AR * (process[t - 1] - mu[t-1]])
     # We therefore need the values of mu at the end of the training set
     # to correctly propagate the process model forward
     if(use_lv & attr(object$model_data, 'trend_model') != 'GP'){
       # Get the observed trend predictor matrix
-      Xp_trend_last <- trend_Xp_matrix(newdata = object$obs_data,
-                                       trend_map = object$trend_map,
-                                       series = series,
-                                       mgcv_model = object$trend_mgcv_model)
+      newdata <- trend_map_data_prep(object$obs_data,
+                                     object$trend_map,
+                                     forecast = TRUE)
+      Xp_trend_last <- predict(object$trend_mgcv_model,
+                               newdata = newdata,
+                               type = 'lpmatrix')
 
       # Ensure the last three values are used, in case the obs_data
       # was not supplied in order
-      data.frame(time = object$obs_data$index..time..index,
-                 series = object$obs_data$series,
-                 row_id = 1:length(object$obs_data$index..time..index)) %>%
+      data.frame(time = newdata$index..time..index,
+                 series = newdata$series,
+                 row_id = 1:length(newdata$index..time..index)) %>%
         dplyr::arrange(time, series) %>%
         dplyr::pull(row_id) -> sorted_inds
-
-      linpred_order <- vector(length = 3 * n_series)
-      last_rows <- tail(sort(sorted_inds), 3 * n_series)
-      for(i in seq_along(last_rows)){
-        linpred_order[i] <- which(sorted_inds == last_rows[i])
-      }
+      n_processes <- length(unique(object$trend_map$trend))
+      linpred_order <- tail(sorted_inds, 3 * n_processes)
+      # last_rows <- tail(sorted_inds, 3 * n_processes)
+      # for(i in seq_along(last_rows)){
+      #   linpred_order[i] <- which(sorted_inds == last_rows[i])
+      # }
 
       # Deal with any offsets
       if(!all(attr(Xp_trend_last, 'model.offset') == 0)){

R/get_linear_predictors.R

@@ -73,13 +73,10 @@ obs_Xp_matrix = function(newdata, mgcv_model){
   return(Xp)
 }
 
-
-#' Function to prepare trend linear predictor matrix, ensuring ordering and
-#' indexing is correct with respect to the model structure
+#' Function to prepare trend linear predictor matrix in the presence of a
+#' trend_map
 #' @noRd
-trend_Xp_matrix = function(newdata, trend_map, series = 'all',
-                           mgcv_model){
-
+trend_map_data_prep = function(newdata, trend_map, forecast = FALSE){
   trend_test <- newdata
   trend_indicators <- vector(length = length(trend_test$series))
   for(i in 1:length(trend_test$series)){
@@ -92,6 +89,44 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all',
   trend_test$series <- trend_indicators
   trend_test$y <- NULL
 
+  # Only keep one time observation per trend, in case this is a reduced dimensionality
+  # State-Space model (with a trend_map) and we are forecasting ahead
+  if(forecast){
+    data.frame(series = trend_test$series,
+               time = trend_test$index..time..index,
+               row_num = 1:length(trend_test$index..time..index)) %>%
+      dplyr::group_by(series, time) %>%
+      dplyr::slice_head(n = 1) %>%
+      dplyr::pull(row_num) -> inds_keep
+    inds_keep <- sort(inds_keep)
+
+    if(inherits(trend_test, 'list')){
+      trend_test <- lapply(trend_test, function(x){
+        if(is.matrix(x)){
+          matrix(x[inds_keep,], ncol = NCOL(x))
+        } else {
+          x[inds_keep]
+        }
+
+      })
+    } else {
+      trend_test <- trend_test[inds_keep, ]
+    }
+  }
+
+  return(trend_test)
+}
+
+#' Function to prepare trend linear predictor matrix, ensuring ordering and
+#' indexing is correct with respect to the model structure
+#' @noRd
+trend_Xp_matrix = function(newdata, trend_map, series = 'all',
+                           mgcv_model, forecast = FALSE){
+
+  trend_test <- trend_map_data_prep(newdata,
+                                    trend_map,
+                                    forecast = forecast)
+
   suppressWarnings(Xp_trend  <- try(predict(mgcv_model,
                                             newdata = trend_test,
                                             type = 'lpmatrix'),

R/ppc.mvgam.R

@@ -1,4 +1,4 @@
-#'@title Plot mvgam posterior predictive checks for a specified series
+#'@title Plot mvgam conditional posterior predictive checks for a specified series
 #'@importFrom stats quantile density ecdf formula terms
 #'@importFrom graphics hist abline box rect lines polygon par
 #'@importFrom grDevices rgb
@@ -26,16 +26,17 @@
 #'@param xlab label for x axis.
 #'@param ylab label for y axis.
 #'@param ... further \code{\link[graphics]{par}} graphical parameters.
-#'@details Posterior predictions are drawn from the fitted \code{mvgam} and compared against
+#'@details Conditional posterior predictions are drawn from the fitted \code{mvgam} and compared against
 #'the empirical distribution of the observed data for a specified series to help evaluate the model's
 #'ability to generate unbiased predictions. For all plots apart from `type = 'rootogram'`, posterior predictions
 #'can also be compared to out of sample observations as long as these observations were included as
 #''data_test' in the original model fit and supplied here. Rootograms are currently only plotted using the
 #''hanging' style.
 #'\cr
-#'Note that the predictions used for these plots are those that have been generated directly within
-#'the `mvgam()` model, so they can be misleading if the model included flexible dynamic trend components. For
-#'a broader range of posterior checks that are created using "new data" predictions, see
+#'Note that the predictions used for these plots are *conditional on the observed data*, i.e. they
+#'are those predictions that have been generated directly within
+#'the `mvgam()` model. They can be misleading if the model included flexible dynamic trend components. For
+#'a broader range of posterior checks that are created using *unconditional* "new data" predictions, see
 #'\code{\link{pp_check.mvgam}}
 #'@return A base \code{R} graphics plot showing either a posterior rootogram (for \code{type == 'rootogram'}),
 #'the predicted vs observed mean for the
@@ -759,7 +760,7 @@ ppc.mvgam = function(object, newdata, data_test, series = 1, type = 'hist',
 
 #' Posterior Predictive Checks for \code{mvgam} Objects
 #'
-#' Perform posterior predictive checks with the help
+#' Perform unconditional posterior predictive checks with the help
 #' of the \pkg{bayesplot} package.
 #'
 #' @aliases pp_check
@@ -777,11 +778,14 @@ ppc.mvgam = function(object, newdata, data_test, series = 1, type = 'hist',
 #' @return A ggplot object that can be further
 #'  customized using the \pkg{ggplot2} package.
 #'
-#' @details For a detailed explanation of each of the ppc functions,
+#' @details Unlike the conditional posterior checks provided by \code{\link{ppc}},
+#' This function computes *unconditional* posterior predictive checks (i.e. it generates
+#' predictions for fake data without considering the true observations associated with those
+#' fake data). For a detailed explanation of each of the ppc functions,
 #' see the \code{\link[bayesplot:PPC-overview]{PPC}}
 #' documentation of the \pkg{\link[bayesplot:bayesplot-package]{bayesplot}}
 #' package.
-#' @seealso \code{\link{ppc}} \code{\link{predict.mvgam}}
+#' @seealso \code{\link{ppc}}, \code{\link{predict.mvgam}}
 #' @examples
 #' \dontrun{
 #'simdat <- sim_mvgam(seasonality = 'hierarchical')

R/sysdata.rda

@@ -152,10 +152,13 @@ mvgam_example1 <- mvgam(y ~ s(season, k = 5),
                         samples = 30,
                         chains = 1)
 
-# Univariate process with trend_formula and correlated process errors
+# Univariate process with trend_formula, trend_map and correlated process errors
+trend_map <- data.frame(series = unique(mvgam_examp_dat$data_train$series),
+                        trend = c(1,1,2))
 mvgam_example2 <- mvgam(y ~ 1,
                         trend_formula = ~ s(season, k = 5),
                         trend_model = RW(cor = TRUE),
+                        trend_map = trend_map,
                         family = gaussian(),
                         data = mvgam_examp_dat$data_train,
                         burnin = 300,

man/pp_check.mvgam.Rd

@@ -42,3 +42,4 @@ gaus_data <- sim_mvgam(family = gaussian(),
                        mu = c(-1, 0, 1),
                        trend_rel = 0.5,
                        prop_missing = 0.2)
+

man/ppc.mvgam.Rd

@@ -107,7 +107,6 @@ test_that("binomial() post-processing works", {
   expect_no_error(ppc(mod, type = 'pit'))
   expect_no_error(ppc(mod, type = 'cdf'))
   expect_no_error(ppc(mod, type = 'rootogram'))
-
   expect_no_error(plot(mod, type = 'residuals'))
 
   expect_no_error(plot_mvgam_series(object = mod))