Training dataset - forecastML::create_lagged_df

• We have 3 datasets for training models that forecast 1, 1 to 7, and 1 to 30 days into the future. We’ll view the 1-day-ahead training data below.

outcome_col <- 1  # The column position of our 'wind_spd' outcome (after removing the 'date' column).

horizons <- c(1, 7, 30)  # Forecast 1, 1:7, and 1:30 days into the future.

lookback <- c(1:30, 360:370)  # Features from 1 to 30 days in the past and annually.

dates <- data$date  # Grouped time series forecasting requires dates.
data$date <- NULL  # Dates, however, don't need to be in the input data.

frequency <- "1 day"  # A string that works in base::seq(..., by = "frequency").

dynamic_features <- c("day", "year")  # Features that change through time but which will not be lagged.

groups <- "buoy_id"  # 1 forecast for each group or buoy.

static_features <- c("lat", "lon")  # Features that do not change through time.

type <- "train"  # Create a model-training dataset.

data_train <- forecastML::create_lagged_df(data, type = type, outcome_col = outcome_col,
                                           horizons = horizons, lookback = lookback,
                                           dates = dates, frequency = frequency,
                                           dynamic_features = dynamic_features,
                                           groups = groups, static_features = static_features, 
                                           use_future = FALSE)

DT::datatable(head(data_train$horizon_1), options = list(scrollX = TRUE))

• The plot below shows the feature map for any lagged features across forecast horizons. Here, we set all non-dynamic and non-static features to have the same lags (refer to the custom lags vignette to see how this could be modified). Notice that features that don’t support direct forecasting to the given horizon–e.g., lags of 1 to 29 days for the 30-day-horizon model–are silently dropped.

p <- plot(data_train)  # plot.lagged_df() returns a ggplot object.
p <- p + geom_tile(NULL)  # Remove the gray border for a cleaner plot.
p

CV setup - forecastML::create_windows

• We’ll model with 3 validation datasets. Given that our measurements are taken daily, we’ll set the skip = 730 argument to skip 2 years between validation datasets. Custom validation windows could be defined with vectors of start and stop dates given to window_start and window_stop.

windows <- forecastML::create_windows(data_train, window_length = 365, skip = 730,
                                      include_partial_window = FALSE)

p <- plot(windows, data_train) + theme(legend.position = "none")
p

• Now we’ll use the group_filter = "buoy_id == 1" argument to get a closer look at 1 of our 14 time series. The user-supplied filter is passed to dplyr::filter() internally.

p <- plot(windows, data_train, group_filter = "buoy_id == 1") + 
  theme(legend.position = "none")
p

User-defined modeling function

• A user-defined wrapper function for model training that takes the following arguments:
  • 1: A horizon-specific data.frame made with create_lagged_df(..., type = "train") (e.g., my_lagged_df$horizon_h),
  • 2: optionally, any number of additional named arguments which can be passed as ‘…’ in train_model()
  • and returns a model object or list containing a model that will be passed into the user-defined predict() function.

Any data transformations, hyperparameter tuning, or inner loop cross-validation procedures should take place within this function, with the limitation that it ultimately needs to return() a model suitable for the user-defined predict() function; a list can be returned to capture meta-data such as pre-processing pipelines or hyperparameter results.

• Notice that the xgboost-specific input datasets are created within this wrapper function.

# The value of outcome_col can also be set in train_model() with train_model(outcome_col = 1).
model_function <- function(data, outcome_col = 1) {
  
  # xgboost cannot handle missing outcomes data.
  data <- data[!is.na(data[, outcome_col]), ]

  indices <- 1:nrow(data)
  
  set.seed(224)
  train_indices <- sample(1:nrow(data), ceiling(nrow(data) * .8), replace = FALSE)
  test_indices <- indices[!(indices %in% train_indices)]

  data_train <- xgboost::xgb.DMatrix(data = as.matrix(data[train_indices, 
                                                           -(outcome_col), drop = FALSE]),
                                     label = as.matrix(data[train_indices, 
                                                            outcome_col, drop = FALSE]))

  data_test <- xgboost::xgb.DMatrix(data = as.matrix(data[test_indices, 
                                                          -(outcome_col), drop = FALSE]),
                                    label = as.matrix(data[test_indices, 
                                                           outcome_col, drop = FALSE]))

  params <- list("objective" = "reg:linear")
  watchlist <- list(train = data_train, test = data_test)
  
  set.seed(224)
  model <- xgboost::xgb.train(data = data_train, params = params, 
                              max.depth = 8, nthread = 2, nrounds = 30,
                              metrics = "rmse", verbose = 0, 
                              early_stopping_rounds = 5, 
                              watchlist = watchlist)

  return(model)
}

Model training - forecastML::train_model

• This should take ~1 minute to train our ‘3 forecast horizons’ * ‘3 validation datasets’ = 9 models.

• The user-defined modeling wrapper function could be much more elaborate, in which case many more models could potentially be trained here.

• These models could be trained in parallel on any OS with the very flexible future package by un-commenting the code below and setting use_future = TRUE. To avoid nested parallelization, models are either trained in parallel across forecast horizons or validation windows, whichever is longer (when equal, the default is parallel across forecast horizons).

#future::plan(future::multiprocess)  # Multi-core or multi-session parallel training.

model_results_cv <- forecastML::train_model(lagged_df = data_train,
                                            windows = windows,
                                            model_name = "xgboost",
                                            model_function = model_function, 
                                            use_future = FALSE)

## [04:13:01] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:01] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:06] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:06] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:11] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:11] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:16] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:16] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:21] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:21] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:25] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:25] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:29] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:29] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:30] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:30] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.
## 
## [04:13:32] WARNING: src/objective/regression_obj.cu:213: reg:linear is now deprecated in favor of reg:squarederror.
## [04:13:32] WARNING: src/learner.cc:767: 
## Parameters: { "metrics" } are not used.

• We can access the xgboost model for any horizon or validation window. Here, we show a summary() of the 1-step-ahead model for the first validation window which is 2012.

summary(model_results_cv$horizon_1$window_1$model)

##                 Length Class              Mode       
## handle               1 xgb.Booster.handle externalptr
## raw             334427 -none-             raw        
## best_iteration       1 -none-             numeric    
## best_ntreelimit      1 -none-             numeric    
## best_score           1 -none-             numeric    
## best_msg             1 -none-             character  
## niter                1 -none-             numeric    
## evaluation_log       3 data.table         list       
## call                10 -none-             call       
## params               5 -none-             list       
## callbacks            2 -none-             list       
## feature_names      128 -none-             character  
## nfeatures            1 -none-             numeric

User-defined prediction function

The following user-defined prediction function is needed for each model:

• A wrapper function that takes the following 2 positional arguments:
  • 1: The model returned from the user-defined modeling function (could be a list containing the model).
  • 2: A data.frame of the model features from forecastML::create_lagged_df(..., type = "train").
• and returns a data.frame of predictions with 1 or 3 columns. A 1-column data.frame will produce point forecasts, and a 3-column data.frame can be used to return point, lower, and upper forecasts (column names and order do not matter).

# If 'model' is passed as a named list, the prediction model would be accessed with model$model or model["model"].
prediction_function <- function(model, data_features) {
  x <- xgboost::xgb.DMatrix(data = as.matrix(data_features))
  data_pred <- data.frame("y_pred" = predict(model, x),
                          "y_pred_lower" = predict(model, x) - 2,  # Optional; in practice, forecast bounds are not hard coded.
                          "y_pred_upper" = predict(model, x) + 2)  # Optional; in practice, forecast bounds are not hard coded.
  return(data_pred)
}

Historical model fit

• Here, we’re predicting on our 3 validation datasets.

data_pred_cv <- predict(model_results_cv, prediction_function = list(prediction_function), data = data_train)

• We’ll filter this plot for closer inspection below.

plot(data_pred_cv) + theme(legend.position = "none")

• It’s somewhat difficult to see how we’ve done here, so we’ll use the group_filter and facet arguments to focus on specific buoys.

plot(data_pred_cv, facet = group ~ model, group_filter = "buoy_id %in% c(1, 2, 3)", windows = 1) 

• Here is another plot of the same historical predictions but with facet = group ~ horizon. Use different combinations of model, horizon, group, along with . (on the right hand side of ~) in a formula with ~ to quickly explore results.

plot(data_pred_cv, facet = group ~ horizon, group_filter = "buoy_id %in% c(1, 2, 3)", windows = 1) 

Historical prediction error - forecastML::return_error

• Let’s take a quick look at our historical forecast error for buoys 1:3.

• We’ll look at mean absolute error (a) for each validation window, (b) for each of the direct forecast horizons (collapsed across validation windows), and global error collapsed across validation windows and direct forecast horizons.

data_error <- forecastML::return_error(data_pred_cv)

plot(data_error, type = "window", group_filter = "buoy_id %in% c(1, 2, 3)", metric = "mae")

plot(data_error, type = "horizon", group_filter = "buoy_id %in% c(1, 2, 3)", metric = "mae")

plot(data_error, type = "global", group_filter = "buoy_id %in% c(1, 2, 3)", metric = "mae")

Forecasting with multiple models from nested CV

• We have 3 datasets that support forecasting 1, 1 to 7, and 1 to 30 days into the future. We’ll view the 1-day-ahead forecasting data below.

• Note that the index and horizon columns are removed internally when passed into the user-defined predict() function.

type <- "forecast"  # Create a forecasting dataset for our predict() function.

data_forecast <- forecastML::create_lagged_df(data, type = type, outcome_col = outcome_col,
                                              horizons = horizons, lookback = lookback,
                                              dates = dates, frequency = frequency,
                                              dynamic_features = dynamic_features,
                                              groups = groups, static_features = static_features, 
                                              use_future = FALSE)

DT::datatable(head(data_forecast$horizon_1), options = list(scrollX = TRUE))

Dynamic features and forecasting

• Our dynamic features day and year were not lagged in our modeling dataset. This was the right choice from a modeling perspective; however, in order to forecast ‘h’ steps ahead, we need to know their future values for each forecast horizon. At present, there’s no function in forecastML to autofill the future values of dynamic, non-lagged features so we’ll simply do it manually below.

for (i in seq_along(data_forecast)) {
  data_forecast[[i]]$day <- lubridate::mday(data_forecast[[i]]$index)  # When dates are given, the 'index` is date-based.
  data_forecast[[i]]$year <- lubridate::year(data_forecast[[i]]$index)
}

Forecast

• Now we’ll forecast 1, 1:7, and 1:30 days into the future with predict(..., data = data_forecast).

• The first time step into the future is max(dates) + 1 * frequency. Here, this is 12-31-2018 + 1 * ‘1 day’ or 1-1-2019.

data_forecasts <- predict(model_results_cv, prediction_function = list(prediction_function), data = data_forecast)

• Plots for each model–just xgboost here–and direct forecast horizon. Because there are 3 xgboost models for each direct forecast horizon there are 3 points/lines for each buoy within each direct forecast horizon.

plot(data_forecasts)

plot(data_forecasts, facet = group ~ ., group_filter = "buoy_id %in% 1:3")