# chunk options
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE
)

1 Pendahuluan

Artikel ini dibuat untuk memenuhi tugas Learning by Building (LBB) pada materi TS - Time Series pada Machine Learning Course. Membuat model time series dari data penjualan pizza dalam setahun untuk membuat model time series yang nanti performanya akan diukur dengan dataset itu sendiri dan dianalisa apakah model layak atau tidak untuk membuat prediksi pada periode waktu selanjutnya.

Plato’s Pizza is a Greek-inspired pizza place located in New Jersey, USA. The restaurant is looking for operations improvement, and they think a year worth of data would be helpful for this.

2 Import Data dan Pre-Processing

Pada bagian ini, kita akan coba menjabarkan proses pengambilan data dan proses awal apa saja yang harus dilakukan agar sebuah dataset siap untuk dianalisa.

2.1 Import Library

Kita akan menggunakan beberapa library pada Artikel ini. Berhubung dataset dan artikel ini kemungkinan akan dikembangkan untuk Tugas LBB yang lainnya, kita akan mengupload library yang akan digunakan pada materi selanjutnya sehingga dokumen lebih futureproof.

# load library
options(scipen = 99) #me-non-aktifkan scientific annotation
library(tidyverse) #koleksi beberapa library R
library(dplyr) #grammar of data manipulation
library(readr) #membaca data
library(ggplot2) #plot statis
library(gridExtra) #menambah plot dalam 1 grid
library(lubridate) #treatment date column
library(scales) # mengatur skala pada plot
library(readxl) #import data excel

# load library for time series
library(padr) # untuk padding data
library(forecast) # time series library
library(tseries) # time series library
library(MLmetrics) # kalkulasi error

2.2 Import dan Cek Data

Data ditaruh pada folder data_input pada working directory dimana file R Markdown kita tempatkan.

#import Dataset
pizza_sales <- read_excel("datasets/Data Model - Pizza Sales.xlsx", sheet = "pizza_sales")

#glimpse Dataset
glimpse(pizza_sales)
## Rows: 48,620
## Columns: 12
## $ order_details_id  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ order_id          <dbl> 1, 2, 2, 2, 2, 2, 3, 3, 4, 5, 6, 6, 7, 8, 9, 9, 9, 9…
## $ pizza_id          <chr> "hawaiian_m", "classic_dlx_m", "five_cheese_l", "ita…
## $ quantity          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ order_date        <dttm> 2015-01-01, 2015-01-01, 2015-01-01, 2015-01-01, 201…
## $ order_time        <dttm> 1899-12-31 11:38:36, 1899-12-31 11:57:40, 1899-12-3…
## $ unit_price        <dbl> 13.25, 16.00, 18.50, 20.75, 16.00, 20.75, 16.50, 20.…
## $ total_price       <dbl> 13.25, 16.00, 18.50, 20.75, 16.00, 20.75, 16.50, 20.…
## $ pizza_size        <chr> "M", "M", "L", "L", "M", "L", "M", "L", "M", "M", "S…
## $ pizza_category    <chr> "Classic", "Classic", "Veggie", "Supreme", "Veggie",…
## $ pizza_ingredients <chr> "Sliced Ham, Pineapple, Mozzarella Cheese", "Peppero…
## $ pizza_name        <chr> "The Hawaiian Pizza", "The Classic Deluxe Pizza", "T…

2.3 Contents

This pizza sales dataset make up 12 relevant features:

  • order_id: Unique identifier for each order placed by a table

  • order_details_id: Unique identifier for each pizza placed within each order (pizzas of the same type and size are kept in the same row, and the quantity increases)

  • pizza_id: Unique key identifier that ties the pizza ordered to its details, like size and price

  • quantity: Quantity ordered for each pizza of the same type and size

  • order_date: Date the order was placed (entered into the system prior to cooking & serving)

  • order_time: Time the order was placed (entered into the system prior to cooking & serving)

  • unit_price: Price of the pizza in USD

  • total_price: unit_price * quantity

  • pizza_size: Size of the pizza (Small, Medium, Large, X Large, or XX Large)

  • pizza_type: Unique key identifier that ties the pizza ordered to its details, like size and price

  • pizza_ingredients: ingredients used in the pizza as shown in the menu (they all include Mozzarella Cheese, even if not specified; and they all include Tomato Sauce, unless another sauce is specified)

  • pizza_name: Name of the pizza as shown in the menu

2.4 Cleansing Data

Pada bagian ini, kita akan melakukan beberapa pengecekan yang umum dilakukan pada sebuah dataset sebelum memulai proses data analisis seperti cek ukuran data awal, cek missing value, cek tipe data pada kolom. Lalu setelah itu, kita melakukan treatment untuk berbagai kondisi yang ditemui pada dataset yang dicek.

#cek ukuran data sebelum cleansing
nrow(pizza_sales)
## [1] 48620
#cek missing value
colSums(is.na(pizza_sales))
##  order_details_id          order_id          pizza_id          quantity 
##                 0                 0                 0                 0 
##        order_date        order_time        unit_price       total_price 
##                 0                 0                 0                 0 
##        pizza_size    pizza_category pizza_ingredients        pizza_name 
##                 0                 0                 0                 0
#cek jenis class untuk kolom pizza_ingredients
class(pizza_sales$pizza_ingredients[5])
## [1] "character"

Kita merasa bahwa value di Column “pizza_ingredients” memiliki potensi untuk dianalisa, sehingga akan lebih mudah apabila kita mengkonversikan value tersebut menjadi data list agar memudahkan data processing

#converting value di column 'pizza_ingredients' menjadi list

pizza_sales$pizza_ingredients <- lapply(strsplit(as.character(pizza_sales$pizza_ingredients),
                                                 "[][']|,\\s*"), function(x) x[nzchar(x)])
class(pizza_sales$pizza_ingredients[5])
## [1] "list"
pizza_sales$pizza_ingredients[5]
## [[1]]
## [1] "Tomatoes"         "Red Peppers"      "Jalapeno Peppers" "Red Onions"      
## [5] "Cilantro"         "Corn"             "Chipotle Sauce"   "Garlic"
#mengganti jenis kolom
pizza_sales_clean <-
  pizza_sales %>%
  mutate(
    unit_price = as_factor(unit_price),
    pizza_size = as_factor(pizza_size),
    pizza_category = as_factor(pizza_category),
    pizza_name = as_factor(pizza_name)
  )
glimpse(pizza_sales_clean)
## Rows: 48,620
## Columns: 12
## $ order_details_id  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ order_id          <dbl> 1, 2, 2, 2, 2, 2, 3, 3, 4, 5, 6, 6, 7, 8, 9, 9, 9, 9…
## $ pizza_id          <chr> "hawaiian_m", "classic_dlx_m", "five_cheese_l", "ita…
## $ quantity          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ order_date        <dttm> 2015-01-01, 2015-01-01, 2015-01-01, 2015-01-01, 201…
## $ order_time        <dttm> 1899-12-31 11:38:36, 1899-12-31 11:57:40, 1899-12-3…
## $ unit_price        <fct> 13.25, 16, 18.5, 20.75, 16, 20.75, 16.5, 20.75, 16.5…
## $ total_price       <dbl> 13.25, 16.00, 18.50, 20.75, 16.00, 20.75, 16.50, 20.…
## $ pizza_size        <fct> M, M, L, L, M, L, M, L, M, M, S, S, S, S, S, S, L, L…
## $ pizza_category    <fct> Classic, Classic, Veggie, Supreme, Veggie, Chicken, …
## $ pizza_ingredients <list> <"Sliced Ham", "Pineapple", "Mozzarella Cheese">, <…
## $ pizza_name        <fct> "The Hawaiian Pizza", "The Classic Deluxe Pizza", "T…
#cek ukuran data
nrow(pizza_sales_clean)
## [1] 48620

Data sudah dianggap ‘clean’ dan bisa digunakan untuk proses selanjutnya yaitu data analisis.

3 Analysis Data Time Series

3.1 EDA

Kita membuat eksplorasi Data Time Series pada Dataset ini agar lebih mengenal jenis data time series yang terjadi.

head(pizza_sales_clean)

Dari data di atas, dapat terlihat bahwa kolom order_date adalah kolom date yang sudah sesuai. kita akan membuat dataframe baru yang berisi data penjualan untuk setiap hari dan dataframe ini selanjutnya akan digunakan untuk analisis time series. Untuk volume penjualan kita akan mengambil jumlah pizza yang terjual, bukan jumlah transaksi. Hal ini karena menurut kita akan lebih bermanfaat untuk analisa berbagai kepentingan bisnis.

pizza_day_sales <- pizza_sales_clean %>%
  group_by(order_date) %>%
  summarise(total_pizza = sum(quantity)) %>%
  ungroup() %>%
  arrange(order_date) %>%
  padr::pad() #memastikan tidak ada tanggal yang terlewat
anyNA(pizza_day_sales)
## [1] TRUE
pizza_day_sales[is.na(pizza_day_sales$total_pizza),]

Beberapa hari di atas merupakan tanggal dimana tidak ada data penjualan, kita asumsikan bahwa di tanggal tersebut toko tutup sehingga penjualannya = 0

#imputasi data ke hari libur
pizza_day_sales[is.na(pizza_day_sales$total_pizza), "total_pizza"] <- 0
anyNA(pizza_day_sales)
## [1] FALSE

3.2 Visualisasi Data Time Series

Selanjutnya adalah pembuatan data time series. Berdasarkan data di atas, kita ingin melihat pola pergerakan data dengan rentang waktu bulanan. Karena kita ingin mendeteksi data seasonal mingguan lebih baik, kita menggunakan (7 x 4) untuk frequensinya

#pembuatan data time series
pizza_ts <- ts(
  data = pizza_day_sales$total_pizza, # yang kita masukkan adalah colomn dari target kita 
  frequency = (7*4) # menggunakan frequensi bulanan 7*4
)

# visualisasi dan inspect pola data
autoplot(pizza_ts)

Insight :

  • Data tidak memiliki trend naik atau turun, dari plot diatas dapat dilihat bahwa data memiliki sifat stationer atau bergerak di sekitar nilai rata-rata.

  • Terdapat banyak data yang bernilai 0 di bulan 10 hinggal bulan 13

Analisa lebih lanjut menggunakan plot hasil decompose untuk menentukan sifat-sifat pada data tersebut.

#decompose time series
pizza_dc <- decompose(pizza_ts) # yang dimasukkan object TS

# plot decompose
pizza_dc %>% autoplot()

Dari plot terlihat bahwa masih terdapat seasonal pada trend. Maka dari itu dicoba decompose lebih lanjut dengan menggunakan metode multi seasonal

# making 2nd ts object
pizza_msts <- pizza_day_sales$total_pizza %>% msts(seasonal.periods = c(7,7*4))  # multiseasonal ts (mingguan, dan bulanan)
pizza_msts %>% mstl() %>% autoplot()

adf.test(pizza_msts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pizza_msts
## Dickey-Fuller = -5.347, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

Insight :

  • Data tidak memiliki trend naik atau trend turun pada rentang waktu bulan 1 hingga 10, namun terjadi trend turun di bulan 10 hingga 11 dan di bulan 13. Sedangkan rentang waktu bulan 11 - 12 memiliki trend naik

  • Terdapat 2 Data Seasonal dengan pola penjualan mingguan dan pola penjualan bulanan. Sehingga data ini merupakan multi-seasonal data time series.

  • Pola pergerakan data pada trend adalah additive, dimana peningkatannya adalah peningkatan yang linear bukan eksponensial.

  • Data Time Series sudah stationer karena p-value < 0.05, yaitu data bergerak di sekitar rata-rata.

Hipotesa Awal : Apabila pembuatan model time series menggunakan rentang waktu bulan 1-10 sebagai data train dan rentang waktu bulan 10-13 sebagai data test, kita merasa bahwa model tidak akan mampu mendeteksi penurunan dan kenaikan ekstrim yang terjadi di rentang bulan 10-13. Sehingga untuk dijadikan kasus forecasting, data ini tidak cukup lengkap karena akan ada data seasonal tahunan yang terjadi namun tidak terwakili di dalam data.

Mungkin harus diambil data minimal rentang waktu 2 tahun untuk mampu membuat forecasting yang akurat atau memperkecil rentang waktu observasi untuk lebih mempersempit cakupan dari trend data time series yang terjadi.

4 Modelling Data Time Series

Berdasarkan eksplorasi data time series sebelumnya, dapat diketahui bahwa data time series ini merupakan jenis multi-seasonal data sehingga untuk pembuatan modelnya mungkin memerlukan pendekatan berbeda dengan data time series yang hanya memiliki 1 seasonal. Untuk itu kita membuat dua pendekatan berbeda yaitu :

  1. Modelling biasa (dengan asumsi 1 seasonal) dengan Holts Winters dan SARIMA, karena data memiliki trend additif dan terdapat seasonal

  2. Modelling TBATS, TBATS adalah model yang umum digunakan untuk membuat model time series yang memiliki seasonal lebih dari 1.

4.1 Model Holts Winter dan SARIMA

Membagi dataset menjadi data train dan data test untuk validasi dan measuring performance model menggunakan fungsi head. Kenapa fungsi head() ? karena dataset yang diambil harus berurutan sesuai dengan kaidah time series dan tidak random

kita mengambil data satu bulan terakhir sebagai dataset untuk test

pizza_msts_train <- head(pizza_msts, length(pizza_msts) - 7*4)
pizza_msts_test <- tail(pizza_msts,  7*4)

autoplot(pizza_msts_train)

#Pembuatan Model Holt Winters
pizza_hw_model <- pizza_msts_train %>%
  stlm(method = "ets", 
       lambda = "auto") #auto, 0 apabila tidak ada log transform
#Pembuatan Model Seasonal ARIMA dengan Auto
pizza_sarima_model <- pizza_msts_train %>%
  auto.arima(D=0, seasonal = T)

kita tidak membuat model SARIMA dengan tuning manual karena dataset merupakan multi seasonal sehingga menurut kita model SARIMA pada dasarnya tidak akan menghasilkan model yang optimal, sehingga kita hanya menggunakan model Auto SARIMA

4.2 Model TBATS

TBATS adalah sebuah model yang dapat digunakan untuk memprediksi data time series. TBATS merupakan singkatan dari “Troncoso-Brockwell-Akima-Tsay”. Model ini dikembangkan oleh Peter J. Brockwell dan Rafael D.S. Troncoso, dan dapat digunakan untuk memprediksi data time series yang memiliki pola seperti musiman, trend, dan siklikal. Model TBATS menggabungkan beberapa metode analisis time series, termasuk model ARIMA dan Fourier, dan dapat memberikan hasil prediksi yang lebih akurat dibandingkan dengan model-model terpisah.

#Pembuatan Model TBATS
pizza_tbats_model <- pizza_msts_train %>% 
  tbats(use.box.cox = NULL, 
                  use.trend = NULL, 
                  use.damped.trend = NULL)

5 Model Forecasting and Evaluation

5.1 Forecasting Each Model

Membuat Forecast sesuai dengan rentang waktu yang ada di data test agar perbandingan benar-benar sesuai.

#forecast Holt Winters
hw_forecast <- forecast(pizza_hw_model, h=length(pizza_msts_test))

#forecast Auto SARIMA
sarima_forecast <- forecast(pizza_sarima_model, h=length(pizza_msts_test))

#forecast TBATS
tbats_forecast <- forecast(pizza_tbats_model, h=length(pizza_msts_test))

5.2 Visualisasi Hasil Forecast

autoplot(hw_forecast$mean, series = "Holt Winters", aes(x = hw_forecast$mean)) + 
  autolayer(pizza_msts_train, color = "black") + 
  autolayer(hw_forecast, color = "red") + 
  autolayer(pizza_msts_test, series = "Actual", color = "black") +
  labs(subtitle = "Pizza Sold for a Year",
       y = "Pizza Sold", 
       x = "Month Number")

Insight :

  • Berdasarkan data di atas, model tidak terlalu baik dalam memprediksi dan range untuk forecast masih terlalu lebar.
autoplot(sarima_forecast$mean, series = "Auto SARIMA", aes(x = sarima_forecast$mean)) + 
  autolayer(pizza_msts_train, color = "black") + 
  autolayer(sarima_forecast, color = "red") + 
  autolayer(pizza_msts_test, series = "Actual", color = "black") +
  labs(subtitle = "Pizza Sold for a Year",
       y = "Pizza Sold", 
       x = "Month Number")

Insight :

  • Dibandingkan hasil forecast Holt Winters, Auto Sarima mampu menghasilkan range forecast yang lebih sempit, namun untuk hasil prediksi data rata-rata nya masih meleset di akhir-akhir data.
autoplot(tbats_forecast$mean, series = "TBATS", aes(x = tbats_forecast$mean)) + 
  autolayer(pizza_msts_train, color = "black") + 
  autolayer(tbats_forecast, color = "red") + 
  autolayer(pizza_msts_test, series = "Actual", color = "black") +
  labs(subtitle = "Pizza Sold for a Year",
       y = "Pizza Sold", 
       x = "Month Number")

Insight :

  • Model mampu memprediksi penurunan dan kenaikan data lebih baik dari Model Holt Winters dan Model Auto SARIMA

5.3 Perbandingan Antar Model

kita mengambil data rata-rata untuk tiap forecast dari tiap model untuk membandingkan performa model terhadap aktual

#membuat Dataframe untuk mempermudah plotting Data
accuracyData <- data.frame(datetime= tail(pizza_day_sales$order_date,7*4),
  hw_value = as.vector(hw_forecast$mean), 
  sarima_value = as.vector(sarima_forecast$mean), 
  tbats_value = as.vector(tbats_forecast$mean),
  actual = as.vector(pizza_msts_test)
)
#plotting perbandingan hasil prediksi dan aktual
accuracyData %>% 
 ggplot() +
  geom_line(aes(y = hw_value,
                x = accuracyData$datetime,
                colour = "Holt Winters"),
            size = 1,
            show.legend = TRUE
            ) + 
  geom_line(aes(y = sarima_value,
                x = accuracyData$datetime,
                colour = "Auto Sarima"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = tbats_value,
                x = accuracyData$datetime,
                colour = "TBATS"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = tbats_value,
                x = accuracyData$datetime,
                colour = "TBATS"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = actual,
                x = accuracyData$datetime,
                colour = "Actual"),
            size = 1,
            show.legend = TRUE
            ) +
  labs(
    title = "Daily Pizza Sold",
    x = "Date",
    y = "Pizza Sold"
  )

Subsetting Data dimana penjualan pizza = 0 untuk mengetahui hari libur.

pizza_day_sales[pizza_day_sales$total_pizza==0,]

Insight :

  • Semua model tidak mampu memprediksi penurunan di tanggal 25 Desember karena memang pada tanggal tersebut Toko Pizza Libur
#Imputasi Data 0 pada hari libur di Value Prediksi

accuracyData[22,2] <- 0
accuracyData[22,3] <- 0
accuracyData[22,4] <- 0
accuracyData %>% 
 ggplot() +
  geom_line(aes(y = hw_value,
                x = accuracyData$datetime,
                colour = "Holt Winters"),
            size = 1,
            show.legend = TRUE
            ) + 
  geom_line(aes(y = sarima_value,
                x = accuracyData$datetime,
                colour = "Auto Sarima"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = tbats_value,
                x = accuracyData$datetime,
                colour = "TBATS"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = tbats_value,
                x = accuracyData$datetime,
                colour = "TBATS"),
            size = 1,
            show.legend = TRUE
            ) +
  geom_line(aes(y = actual,
                x = accuracyData$datetime,
                colour = "Actual"),
            size = 1,
            show.legend = TRUE
            ) +
  labs(
    title = "Daily Pizza Sold",
    x = "Date",
    y = "Pizza Sold"
  )

Insight Plot :

  • Berdasarkan plotting yang baru, terlihat bahwa nilai rata-rata forecast dari ketiga model cukup baik untuk memprediksi kenaikan dan penurunan dari penjualan pizza

  • Tanggal 25 Desember merupakan tanggal libur Natal, kemungkinan toko pizza tidak buka di hari tersebut

Insight Bisnis :

  • Tanggal 25 hingga 30 Desember hasil prediksi sangat jauh dari data aktual, hal ini mungkin disebabkan karena setelah Natal orang-orang lebih memilih untuk memasak di rumah menikmati hidangan yang biasa dihidangkan pada hari libur natal. Ini bisa jadi adalah data seasonal tahunan yang tidak tertangkap oleh modelling karena data kurang.

  • Tanggal 31 Desember penjualan meningkat kembali, kemungkinan karena tanggal tersebut orang-orang mengadakan pesta tahun baru yang identik dengan **standing party* *sehingga Pizza merupakan menu yang sangat ideal.

5.4 Model Evaluation

kita menggunakan metode RMSE atau Root Mean Squared Error untuk mengukur performa model. Untuk pengukuran hasil model, kita menggunakan data yang sudah diimputasi nilai 0 untuk hari libur tanggal 25 Desember.

# model evaluation: root mean squared error (RMSE)
data.frame(HoltWinters = RMSE(accuracyData$hw_value, accuracyData$actual),
           SARIMA = RMSE(accuracyData$sarima_value, accuracyData$actual),
           TBATS = RMSE(accuracyData$tbats_value, accuracyData$actual))

Berdasarkan hasil RMSE, model TBATS memiliki nilai yang paling baik dibandingkan model yang lainnya. Error yang terjadi yaitu 23,6.

5.5 Conclusion

Berdasarkan hasil analisis di artikel ini, dapat terlihat bahwa untuk menganalisa data time series yang memiliki multi seasonal perlu ditreatment dengan modelling yang berbeda, salah satu contohnya yaitu menggunakan TBATS.

Dari hasil analisis model TBATS menghasilkan nilai RMSE yaitu 23,6 yang bisa diinterpretasikan bahwa ketika memprediksi nilai penjualan pizza ke depan, terdapat potensi error 24 pizza dari hasil rata-rata forecast dengan nilai keyakinan 95%.

Terdapat kemungkinan seasonal tahunan yang tidak terdeteksi karena dataset yang tersedia hanya 1 tahun.

6 Assumption Check

kita menggunakan model TBATS untuk cek asumsi, karena model tersebut yang performanya paling baik.

6.1 Normality : Shapiro Test

  • H0 : residuals terdistribusi normal

  • H1 : residuals tidak terdistribusi normal

  • p-value < 0.05; tolak H0; terima H1

shapiro.test(tbats_forecast$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  tbats_forecast$residuals
## W = 0.81873, p-value < 0.00000000000000022
hist(tbats_forecast$residuals)

plot(tbats_forecast$residuals)

Insight :

  • H0, Data residual model tidak terdistribusi Normal

Data mungkin menjadi tidak terdistribusi normal karena imputasi nilai = 0 pada hari libur (data real), sedangkan pada aktualnya apabila tidak libur bisa saja Toko Pizza menjual pizza tidak jauh dari angka prediksi. Hal ini bisa diatasi dengan menginput nilai mean atau median pada data penjualan di hari libur.

6.2 Autocorrelation : Box.test - Ljng-Box

  • H0 : Tidak ada Autocorrelation pada Residuals

  • H1 : Ada Autocorrelation pada Residual

  • p-value < 0.05; tolak H0; terima H1

Box.test(tbats_forecast$residuals, type = "Ljung-Box") 
## 
##  Box-Ljung test
## 
## data:  tbats_forecast$residuals
## X-squared = 2.6284, df = 1, p-value = 0.105

Insight :

  • Residual tidak ada autocorrelation terhadap data penjualan pizza.