add check to ensure test data are in correct order when forecasting

Author: Nicholas Clark

R/add_nmixture.R

@@ -719,27 +719,3 @@ add_nmix_posterior = function(model_output,
 
   return(model_output)
 }
-
-# Get list object into correct order in case it is not already
-sort_data = function(obs_data){
-  if(inherits(obs_data, 'list')){
-    obs_data_arranged <- obs_data
-    temp_dat = data.frame(time = obs_data$time,
-                          series = obs_data$series) %>%
-      dplyr::mutate(index = dplyr::row_number()) %>%
-      dplyr::arrange(series, time)
-    obs_data_arranged <- lapply(obs_data, function(x){
-      if(is.matrix(x)){
-        matrix(x[temp_dat$index,], ncol = NCOL(x))
-      } else {
-        x[temp_dat$index]
-      }
-    })
-    names(obs_data_arranged) <- names(obs_data)
-  } else {
-    obs_data_arranged <- obs_data %>%
-      dplyr::arrange(series, time)
-  }
-
-  return(obs_data_arranged)
-}

R/forecast.mvgam.R

@@ -144,7 +144,8 @@ forecast.mvgam = function(object, newdata, data_test,
                                type = type,
                                series = series,
                                data_test = data_test,
-                               n_cores = n_cores)
+                               n_cores = n_cores,
+                               ...)
 
     # Extract hindcasts and forecasts into the correct format
     if(series == 'all'){
@@ -757,7 +758,10 @@ forecast_draws = function(object,
                           data_test,
                           n_cores = 1,
                           n_samples,
-                          ending_time){
+                          ending_time,
+                          b_uncertainty = TRUE,
+                          trend_uncertainty = TRUE,
+                          obs_uncertainty = TRUE){
 
   # Check arguments
   validate_pos_integer(n_cores)
@@ -769,9 +773,10 @@ forecast_draws = function(object,
     s_name <- levels(data_test$series)[series]
   }
 
-  # Generate the observation model linear predictor matrix
+  # Generate the observation model linear predictor matrix,
+  # ensuring the test data is sorted correctly (by time and then series)
   if(inherits(data_test, 'list')){
-    Xp <- obs_Xp_matrix(newdata = data_test,
+    Xp <- obs_Xp_matrix(newdata = sort_data(data_test),
                         mgcv_model = object$mgcv_model)
 
     if(series != 'all'){
@@ -800,15 +805,16 @@ forecast_draws = function(object,
       Xp <- obs_Xp_matrix(newdata = series_test,
                           mgcv_model = object$mgcv_model)
     } else {
-      Xp <- obs_Xp_matrix(newdata = data_test,
+      Xp <- obs_Xp_matrix(newdata = sort_data(data_test),
                           mgcv_model = object$mgcv_model)
       series_test <- NULL
     }
   }
 
-  # Generate linear predictor matrix from trend mgcv model
+  # Generate linear predictor matrix from trend mgcv model, ensuring
+  # the test data is sorted correctly (by time and then series)
   if(!is.null(object$trend_call)){
-    Xp_trend <- trend_Xp_matrix(newdata = data_test,
+    Xp_trend <- trend_Xp_matrix(newdata = sort_data(data_test),
                                 trend_map = object$trend_map,
                                 series = series,
                                 mgcv_model = object$trend_mgcv_model)
@@ -984,7 +990,10 @@ forecast_draws = function(object,
                                   'series_test',
                                   'Xp',
                                   'Xp_trend',
-                                  'fc_horizon'),
+                                  'fc_horizon',
+                                  'b_uncertainty',
+                                  'trend_uncertainty',
+                                  'obs_uncertainty'),
                           envir = environment())
   parallel::clusterExport(cl = cl,
                           unclass(lsf.str(envir = asNamespace("mvgam"),
@@ -998,18 +1007,31 @@ forecast_draws = function(object,
     samp_index <- i
 
     # Sample beta coefs
-    betas <- betas[samp_index, ]
+    if(b_uncertainty){
+      betas <- betas[samp_index, ]
+    } else {
+      betas <- betas[1, ]
+    }
 
     if(!is.null(betas_trend)){
-      betas_trend <- betas_trend[samp_index, ]
+      if(b_uncertainty){
+        betas_trend <- betas_trend[samp_index, ]
+      } else {
+        betas_trend <- betas_trend[1, ]
+      }
     }
 
     # Return predictions
     if(series == 'all'){
 
       # Sample general trend-specific parameters
-      general_trend_pars <- extract_general_trend_pars(trend_pars = trend_pars,
-                                                       samp_index = samp_index)
+      if(trend_uncertainty){
+        general_trend_pars <- extract_general_trend_pars(trend_pars = trend_pars,
+                                                         samp_index = samp_index)
+      } else {
+        general_trend_pars <- extract_general_trend_pars(trend_pars = trend_pars,
+                                                         samp_index = 1)
+      }
 
       if(use_lv || trend_model %in% c('VAR1', 'PWlinear', 'PWlogistic')){
         if(trend_model == 'PWlogistic'){
@@ -1107,9 +1129,17 @@ forecast_draws = function(object,
           # Family-specific parameters
           family_extracts <- lapply(seq_along(family_pars), function(x){
             if(is.matrix(family_pars[[x]])){
-              family_pars[[x]][samp_index, series]
+              if(obs_uncertainty){
+                family_pars[[x]][samp_index, series]
+              } else {
+                family_pars[[x]][1, series]
+              }
             } else {
-              family_pars[[x]][samp_index]
+              if(obs_uncertainty){
+                family_pars[[x]][samp_index]
+              } else {
+                family_pars[[x]][1]
+              }
             }
           })
           names(family_extracts) <- names(family_pars)

R/sim_mvgam.R

@@ -106,6 +106,11 @@ sim_mvgam = function(T = 100,
     use_lv <- FALSE
   }
 
+  if(trend_model %in% c('RWcor', 'AR1cor', 'AR2cor', 'AR3cor')){
+    warning(paste0('Simulation of correlated AR or RW trends not yet supported.\n',
+                   'Reverting to uncorrelated trends'))
+  }
+
   if(missing(trend_rel)){
     trend_rel <- prop_trend
   }
@@ -219,31 +224,32 @@ sim_mvgam = function(T = 100,
   }
 
   # Set trend parameters
-  if(trend_model == 'RW'){
+  if(trend_model %in% c('RW', 'RWcor')){
     ar1s <- rep(1, n_lv)
     ar2s <- rep(0, n_lv)
     ar3s <- rep(0, n_lv)
   }
 
-  if(trend_model == 'AR1'){
+  if(trend_model %in% c('AR1', 'AR1cor')){
     ar1s <- rnorm(n_lv, sd = 0.5)
     ar2s <- rep(0, n_lv)
     ar3s <- rep(0, n_lv)
   }
 
-  if(trend_model == 'AR2'){
+  if(trend_model %in% c('AR2', 'AR2cor')){
     ar1s <- rnorm(n_lv, sd = 0.5)
     ar2s <- rnorm(n_lv, sd = 0.5)
     ar3s <- rep(0, n_lv)
   }
 
-  if(trend_model == 'AR3'){
+  if(trend_model %in% c('AR3', 'AR3cor')){
     ar1s <- rnorm(n_lv, sd = 0.5)
     ar2s <- rnorm(n_lv, sd = 0.5)
     ar3s <- rnorm(n_lv, sd = 0.5)
   }
 
-  if(trend_model %in% c('RW', 'AR1', 'AR2', 'AR3', 'VAR1', 'VAR1cor')){
+  if(trend_model %in% c('RW', 'AR1', 'AR2', 'AR3',
+                        'VAR1', 'VAR1cor')){
     # Sample trend drift terms so they are (hopefully) not too correlated
     if(drift){
       trend_alphas <- rnorm(n_lv, sd = 0.15)

R/validations.R

@@ -64,6 +64,44 @@ validate_series_time = function(data, name = 'data'){
   return(data)
 }
 
+#'@noRd
+#' Get data objects into correct order in case it is not already
+sort_data = function(data, series_time = FALSE){
+  if(inherits(data, 'list')){
+    data_arranged <- data
+    if(series_time){
+      temp_dat = data.frame(time = data$time,
+                            series = data$series) %>%
+        dplyr::mutate(index = dplyr::row_number()) %>%
+        dplyr::arrange(series, time)
+    } else {
+      temp_dat = data.frame(time = data$time,
+                            series = data$series) %>%
+        dplyr::mutate(index = dplyr::row_number()) %>%
+        dplyr::arrange(time, series)
+    }
+
+    data_arranged <- lapply(data, function(x){
+      if(is.matrix(x)){
+        matrix(x[temp_dat$index,], ncol = NCOL(x))
+      } else {
+        x[temp_dat$index]
+      }
+    })
+    names(data_arranged) <- names(data)
+  } else {
+    if(series_time){
+      data_arranged <- data %>%
+        dplyr::arrange(series, time)
+    } else {
+      data_arranged <- data %>%
+        dplyr::arrange(time, series)
+    }
+  }
+
+  return(data_arranged)
+}
+
 #'@importFrom rlang warn
 #'@noRd
 validate_family = function(family, use_stan = TRUE){