45  Machine Learning xgboost

As usual, we start by uploading the libraries we’re going to use.

We add here 2 new libraries : caret and timetk.

# ETL
library(tidyverse)
library(readxl)
library(reshape2)
library(scales)
library(data.table)
library(zoo)

# Forecasting
library(xgboost)
library(caret)
library(timetk)

# Charts
library(highcharter)

The caret package is for building machine learning models, including random forest models. It integrates various activities such as data preprocessing, model training, hyperparameter tuning, and evaluation.

The timetk is designed to make time series analysis easier, providing tools for visualization, data wrangling, and feature engineering of time series data.

46 Set Up Raw Data

Upload Raw Data

We’re going to upload a dataset from the file “actual_sales_data.csv”. This data frame contains historical sales quantities per product and per country.

We will then calculate some statistical forecasts, using the Machine Learning approach xgboost.

# Upload dataset

# Define the URL of the raw CSV file
url <- "https://raw.githubusercontent.com/nguyennico/sales_forecasting_practice/main/actual_sales_data.csv"

# Read the CSV file from the URL
df1 <- read.csv(url)

# pivot
df1 <- df1 |> gather(key = "period", 
                     value = "sales", 
                     3:length(df1))

# remove the "X" in front of the period
df1$period <- gsub("X", "",df1$period)

# Format Date
df1$period <- as.Date(df1$period, format = '%m.%d.%Y')

# keep results
initial_data <- df1

head(initial_data)

  country product_description     period sales
1   Spain            ProductA 2016-01-01  1370
2  Mexico            ProductB 2016-01-01     0
3  Brazil            ProductB 2016-01-01    12
4  Brazil            ProductC 2016-01-01   369
5 Germany            ProductD 2016-01-01     0
6  France            ProductA 2016-01-01   263

Transform

# set a working df
df1 <- initial_data

# create a field DFU
df1$DFU <- paste(df1$country, df1$product_description, sep="_")

# filter on period
df1 <- df1 |> filter(period <= "2019-06-01")

# keep only needed variables
df1 <- df1 |> select(DFU, period, sales)

# rename
df1 <- df1 |> rename(product = DFU)

# count number of SKUs
list_skus <- unique(df1$product)

# Ensure product is a factor
df1$product <- factor(df1$product)

# replace missing values by zero
df1$sales <- df1$sales |> replace_na(0)

# keep results
sales_data <- df1

Check levels :

# Check the levels of the 'product' factor
levels(df1$product)

[1] "Australia_ProductC" "Brazil_ProductB"    "Brazil_ProductC"   
[4] "France_ProductA"    "France_ProductB"    "Germany_ProductD"  
[7] "Mexico_ProductB"    "Spain_ProductA"    

47 Prepare the Train & Test Data

We’re going to create 2 datasets, splitting the initial data frame according to a “split date” :

• for training : before the split date.

• for testing : after the split date.

47.1 define split date

# Split the data into training and test sets based on time
split_date <- as.Date("2018-12-01")

split_date

[1] "2018-12-01"

47.2 create 2 datasets

# set a working df
df1 <- sales_data

# Sort the data by period
df1 <- df1 |> arrange(period)


# create training dataset
train_data <- df1 |> filter(period < split_date)

# create testing dataset
test_data <- df1 |> filter(period >= split_date)

48 Feature Engineering

Create dates features such as year, month, day, etc., from the period variable.

# training dataset
train_data$year <- year(train_data$period)
train_data$month <- month(train_data$period)
train_data$day <- day(train_data$period)

# testing dataset
test_data$year <- year(test_data$period)
test_data$month <- month(test_data$period)
test_data$day <- day(test_data$period)

49 Prepare Data for XGBoost

Now, let’s convert the data into a matrix format, which is suitable to apply the xgboost algorithm.

#---------------------------
# Create model matrices
#---------------------------

# training dataset
train_matrix <- model.matrix(sales ~ product + year + month + day, 
                             data = train_data)[, -1]

# testing dataset
test_matrix <- model.matrix(sales ~ product + year + month + day, 
                            data = test_data)[, -1]


# Convert the target variable to a numeric vector
train_labels <- train_data$sales
test_labels <- test_data$sales

# Create DMatrix objects for XGBoost
dtrain <- xgb.DMatrix(data = train_matrix, label = train_labels)
dtest <- xgb.DMatrix(data = test_matrix, label = test_labels)

50 Train the XGBoost Model

50.1 Parameters

Set the parameters and train the model. There are 5 parameters :

• objective.

  • specifies the learning task and the corresponding learning objective.

  • the value “reg:squarederror” is used for regression tasks and it minimizes the squared error between the predicted and actual values.

• eta.

  • this is the learning rate.

  • it controls the step size at each iteration while moving toward a minimum of the loss function. A lower value makes the model more robust but requires more rounds to converge. Typical values range from 0.01 to 0.3.

• max_depth.

  • defines the maximum depth of each tree.

  • increasing this value makes the model more complex and more likely to overfit. Typical values range from 3 to 10.

• subsample.

  • specifies the fraction of the training data to be used for growing each tree.

  • setting it to 0.8 means that 80% of the training data is used for each tree. This can help prevent overfitting.

• colsample_bytree.

  • specifies the fraction of features to be randomly sampled for each tree

  • setting it to 0.8 means that 80% of the features are used for each tree. This can also help prevent overfitting.

Let’s create a list called “params” with the values for those parameters.

params <- list(
  objective = "reg:squarederror", # for regression
  eta = 0.1,
  max_depth = 6,
  subsample = 0.8,
  colsample_bytree = 0.8
)

50.2 Model Training

We are going to apply the function xgb.train() from the package xgboost, with 6 inputs :

• params: the list of parameters defined above.

• data: the training data in the form of a DMatrix object (dtrain).

• nrounds.

  • the number of boosting rounds, i.e., the number of trees to be built.

  • here, it is set to 100.

• watchlist.

  • a list of DMatrix objects to be evaluated during training.

  • in this case, it includes both the training data (dtrain) and the test data (dtest) created previously. This allows the model to be evaluated on both datasets at each round.

• early_stopping_rounds.

  • specifies the number of rounds with no improvement after which training will be stopped.

  • if the performance on the test set does not improve for 10 consecutive rounds, training will stop early. This helps to prevent overfitting.

• print_every_n.

  • specifies how often (in terms of boosting rounds) the progress will be printed.

  • here, it is set to print every 10 rounds.

xgb_model <- xgb.train(
  params = params,
  data = dtrain,
  nrounds = 100,
  watchlist = list(train = dtrain, test = dtest),
  early_stopping_rounds = 10,
  print_every_n = 10
)

Multiple eval metrics are present. Will use test_rmse for early stopping.
Will train until test_rmse hasn't improved in 10 rounds.

[1] train-rmse:329.421991   test-rmse:218.915814 
[11]    train-rmse:150.737861   test-rmse:93.670311 
[21]    train-rmse:82.924302    test-rmse:81.936550 
Stopping. Best iteration:
[29]    train-rmse:60.951007    test-rmse:85.789493

[29]    train-rmse:60.951007    test-rmse:85.789493 

In summary :

• the list params: defines the parameters for the XGBoost model.

• the function xgb.train() :

  • trains the XGBoost model using the specified parameters and training data.

  • evaluates the model on the test data.

  • includes early stopping to prevent overfitting.

    • by tuning these parameters and using early stopping, we can improve the performance of the XGBoost model and prevent overfitting.

51 Evaluate the Model

Let’s evaluate the model’s performance on the test set :

• we calculate first the predicted values preds, using the function predict().

• then we calculate the RMSE (Root Mean Square Error) to evaluate the quality of the prediction.

Note :

• low RMSE means good model accuracy, which suggests the model’s predictions are reliable.

• high RMSE means poor performance.

However, whether RMSE is high or low depends on the context and the problem domain. Different fields have different scales and tolerances for errors.

# Make predictions
preds <- predict(xgb_model, dtest)

# Calculate RMSE
rmse <- sqrt(mean((preds - test_labels)^2))
print(paste("RMSE: ", rmse))

[1] "RMSE:  80.6042260130753"

In this case, is a RMSE of 84 good?

We can think about the Relative Error:

• we can consider the RMSE in relation to the average sales value.

• for example :

  • if our average sales per period are 1,000, an RMSE of 84 represents an error of 8.4%, which is pretty good.

  • if our average sales per period are 100, an RMSE of 84 represents an error of 84%, which is not really good.

52 Make Future Predictions

We will create a matrix with all the future periods for all the items, and then apply to it the model xgb_model that we previously calculated.

52.1 Create Future Periods matrix

# get the kist of SKUs as a dataframe
list_skus_data <- as.data.frame(list_skus)

# create a dataframe with the time horizon for forecasts
period_data <- seq.Date(as.Date("2019-01-01"), by = "month", length.out = 12)
period_data <- as.data.frame(period_data)

# cross joint
df1 <- tidyr::crossing(list_skus_data, period_data)

# rename
df1 <- df1 |> rename(product = list_skus,
                     period = period_data)

# Ensure product is a factor with the same levels as the original data
df1$product <- factor(df1$product)

# Create features for future periods
df1$year <- year(df1$period)
df1$month <- month(df1$period)
df1$day <- day(df1$period)

# keep results
future_periods <- df1

glimpse(df1)

Rows: 96
Columns: 5
$ product <fct> Australia_ProductC, Australia_ProductC, Australia_ProductC, Au…
$ period  <date> 2019-01-01, 2019-02-01, 2019-03-01, 2019-04-01, 2019-05-01, 2…
$ year    <int> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 20…
$ month   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8,…
$ day     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

52.2 Convert to matrix

# Convert to matrix
future_matrix <- model.matrix(~ product + year + month + day, 
                              data = future_periods)[, -1]

glimpse(future_matrix)

 num [1:96, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:96] "1" "2" "3" "4" ...
  ..$ : chr [1:10] "productBrazil_ProductB" "productBrazil_ProductC" "productFrance_ProductA" "productFrance_ProductB" ...

52.3 Predict Sales Forecasts

We :

• create an object xgb.DMatrix using the function xgb.DMatrix() on the object future_matrix that we just created.

• apply on this object the model xgb_model previously trained, using the function predict().

# Create DMatrix
dfuture <- xgb.DMatrix(data = future_matrix)

# Predict future sales
future_preds <- predict(xgb_model, dfuture)

# Add predictions to the future_periods dataframe
future_periods$predicted_sales <- future_preds

print(future_periods)

# A tibble: 96 × 6
   product            period      year month   day predicted_sales
   <fct>              <date>     <int> <int> <int>           <dbl>
 1 Australia_ProductC 2019-01-01  2019     1     1            254.
 2 Australia_ProductC 2019-02-01  2019     2     1            254.
 3 Australia_ProductC 2019-03-01  2019     3     1            254.
 4 Australia_ProductC 2019-04-01  2019     4     1            262.
 5 Australia_ProductC 2019-05-01  2019     5     1            261.
 6 Australia_ProductC 2019-06-01  2019     6     1            261.
 7 Australia_ProductC 2019-07-01  2019     7     1            261.
 8 Australia_ProductC 2019-08-01  2019     8     1            261.
 9 Australia_ProductC 2019-09-01  2019     9     1            261.
10 Australia_ProductC 2019-10-01  2019    10     1            261.
# ℹ 86 more rows

We have calculated the forecasts in the variable predicted_sales .

Now, let’s :

• keep only needed variables : product, period, predicted_sales.

• add a “dummy” variable actuals.

We get a data frame that we will call forecasts_data . As we did previously, we will then assemble it with the historical data in the coming parts.

# set a working df
df1 <- future_periods

# keep only needed variables
df1 <- df1 |> select(product, period, predicted_sales)

# rename
df1 <- df1 |> rename(forecasts = predicted_sales)

# add a dummy actuals variable
df1$actuals <- 0

# reorder
df1  <- df1 |> select(product, period, actuals, forecasts)

# keep results
forecasts_data <- df1

glimpse(df1)

Rows: 96
Columns: 4
$ product   <fct> Australia_ProductC, Australia_ProductC, Australia_ProductC, …
$ period    <date> 2019-01-01, 2019-02-01, 2019-03-01, 2019-04-01, 2019-05-01,…
$ actuals   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ forecasts <dbl> 253.63519, 253.63519, 253.63519, 262.28351, 260.90933, 260.9…

53 Assemble

Now we’re going to stack the historical (actual) values with the forecasts we just calculated.

We aim to have a single data frame with especially 2 variables : actuals and forecasts.

We first prepare the actuals, and then stack with the calculated forecasts.

Prepare actuals

# set a working df
df1 <- sales_data

# rename
df1 <- df1 |> rename(actuals = sales)

# add a dummy forecasts variable
df1$forecasts <- 0

# reorder
df1  <- df1 |> select(product, period, actuals, forecasts)

# keep results
actuals_data <- df1

glimpse(df1)

Rows: 336
Columns: 4
$ product   <fct> Spain_ProductA, Mexico_ProductB, Brazil_ProductB, Brazil_Pro…
$ period    <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01, 2016-01-01,…
$ actuals   <int> 1370, 0, 12, 369, 0, 263, 12, 213, 1528, 0, 47, 368, 0, 298,…
$ forecasts <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

Stack

# stack
df1 <- rbind(actuals_data, forecasts_data)

# replace missing values by zero
df1$actuals <- df1$actuals |> replace_na(0)
df1$forecasts <- df1$forecasts |> replace_na(0)

# aggregate
df1 <- df1 |> group_by(product, period) |>
  summarise(actuals = sum(actuals),
            forecasts = sum(forecasts)
            )

# rename
df1 <- df1 |> rename(DFU = product)


# replace any values at zero by NA for a better display
df1$actuals <- if_else(df1$actuals == 0, NA, df1$actuals)
df1$forecasts <- if_else(df1$forecasts == 0, NA, df1$forecasts)


# keep results
Set_Up_ML_boost_Forecasts_data <- df1

glimpse(df1)

Rows: 384
Columns: 4
Groups: DFU [8]
$ DFU       <fct> Australia_ProductC, Australia_ProductC, Australia_ProductC, …
$ period    <date> 2016-01-01, 2016-02-01, 2016-03-01, 2016-04-01, 2016-05-01,…
$ actuals   <dbl> 213, 232, 105, 225, 240, 228, 194, 178, 233, 167, 164, 208, …
$ forecasts <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

54 Add Dates Features

Most of the time we will create some charts to visualize the results.

For this purpose it’s good to add some classic dates features :

• calendar year and month.

• abbreviation of calendar month name, factorised.

# set a working df
df1 <- Set_Up_ML_boost_Forecasts_data

#-----------------------
# Calendar
#-----------------------

# Get Calendar.Month
df1$calendar_month <- month(df1$period)

# Get Calendar.Year
df1$calendar_year <- year(df1$period)


# Add Calendar.Month.abb
df1$calendar_month_abb <- month.abb[df1$calendar_month]

# create a factor
df1$calendar_month_abb <- factor(df1$calendar_month_abb,
                              levels = c("Jan","Feb","Mar","Apr","May","Jun",
                                        "Jul","Aug","Sep", "Oct", "Nov", "Dec"))



# keep results
Set_Up_ML_boost_Forecasts_data <- df1

glimpse(df1)

Rows: 384
Columns: 7
Groups: DFU [8]
$ DFU                <fct> Australia_ProductC, Australia_ProductC, Australia_P…
$ period             <date> 2016-01-01, 2016-02-01, 2016-03-01, 2016-04-01, 20…
$ actuals            <dbl> 213, 232, 105, 225, 240, 228, 194, 178, 233, 167, 1…
$ forecasts          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ calendar_month     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, …
$ calendar_year      <int> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 201…
$ calendar_month_abb <fct> Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, N…

55 Train | Test | Validate

Now let’s check the accuracy of our model :

• with a chart displaying the actuals and the calculated forecasts.

• with the calculation of the MAPE.

55.1 Chart

We can create a chart and compare our calculated forecasts vs the actuals, to validate the results.

Below is a line chart using the library highcharter .

# select product
df1 <- Set_Up_ML_boost_Forecasts_data  |> filter(DFU == "Spain_ProductA")


# chart
highchart() |>
  
  hc_add_series(name = "Actuals", 
                color = "steelblue", 
                data = df1$actuals) |>
  
  hc_add_series(name = "Forecasts", 
                color = "limegreen", 
                data = df1$forecasts) |>
  
  hc_title(text = "ML XGboost Sales Forecasting") |>
  hc_subtitle(text = "in units") |>
  hc_xAxis(categories = df1$period) |> 
  hc_add_theme(hc_theme_google())

Compared to the Time Series calculations in the previous chapters, we can notice that :

• there are some months where the forecasts fit pretty well the actuals : Feb, May and Jun 2019.

• and some months with a significant difference : March and April 2019.

The MAPE should probably be less good.

55.2 Calculation of the MAPE

We get a MAPE of 20%.

Though still a good result, in this example, the classic Time Series calculations provided a better performance.

# select product
df1 <- Set_Up_ML_boost_Forecasts_data  |> filter(DFU == "Spain_ProductA")

# calculate monthly difference
df1$delta <- df1$forecasts - df1$actuals
df1$delta_pc <- df1$delta / df1$actuals


#----------------------------
# calculate the Mean absolute percentage error (MAPE)
#----------------------------

# get the absolute percentage of error
df1$abs_delta_pc <- abs(df1$delta_pc)

# define timeframe for the measurement of the MAPE

# start date
Start.Date<- '2019-01-01'

# end date
End.Date<- '2019-06-01'

# select Period (historical window)
df2 <- df1 |> filter(period >= Start.Date & df1$period <= End.Date)

# calculate the mean
MAPE <- mean(df2$abs_delta_pc)
MAPE <- percent(MAPE)

MAPE

[1] "18%"