# clear-up the environment
rm(list = ls())
# chunk options
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE,
fig.align = "center",
comment = "#>"
)
Artikel ini dibuat untuk memenuhi tugas Learning by Building (LBB) pada materi TS - Time Series pada Machine Learning Course. Materi berikut ditujukan untuk semua orang yang tertarik pada bidang Data Science dan R-Language secara umum. Materi ini dipersilahkan untuk direproduksi, didistribusikan, diterjemahkan, atau diadaptasikan sesuai dengan kebutuhan pembaca.
Artikel ini dibuat melanjutkan dari LBB - Data Visualization : Pizza Sales dimana penulis mencoba memaparkan visualisasi data penjualan pizza dalam setahun.
Pada artikel kali ini, penulis mencoba 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.
knitr::include_graphics("images/pizza-enak.jpg")
Setiap langkah dari metodologi ini akan diikuti oleh penjelasan yang dijabarkan oleh Penulis untuk memperjelas maksud dan tujuan dari setiap langkah yang dilakukan. Selain itu juga dapat bermanfaat untuk pembaca ketika ingin mencoba memulai belajar untuk mengolah sebuah dataset.
Pada bagian ini, Penulis akan coba menjabarkan proses pengambilan data dan proses awal apa saja yang harus dilakukan agar sebuah dataset siap untuk dianalisa.
Penulis akan menggunakan beberapa library pada Artikel ini. Berhubung dataset dan artikel ini kemungkinan akan dikembangkan untuk Tugas LBB yang lainnya, penulis 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
Data ditaruh pada folder data_input pada working directory dimana file R Markdown penulis tempatkan.
#import Dataset
pizza_sales <- read_excel("data_input/datamodel-pizzasales.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…
Deskripsi Kolom :
Berdasarkan glimpse di atas, ada beberapa kolom yang perlu diganti jenisnya untuk memudahkan processing data.
Pada bagian ini, Penulis 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, Penulis 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"
Penulis merasa bahwa value di Column “pizza_ingredients” memiliki potensi untuk dianalisa, sehingga akan lebih mudah apabila Penulis 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.
Penulis 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. Penulis akan membuat dataframe baru
yang berisi data penjualan untuk setiap hari dan dataframe ini
selanjutnya akan digunakan untuk analisis time series. Untuk volume
penjualan penulis akan mengambil jumlah pizza yang terjual, bukan jumlah
transaksi. Hal ini karena menurut penulis 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
Selanjutnya adalah pembuatan data time series. Berdasarkan data di atas, penulis ingin melihat pola pergerakan data dengan rentang waktu bulanan. Karena penulis ingin mendeteksi data seasonal mingguan lebih baik, penulis 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 :
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 :
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, penulis 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.
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 penulis membuat dua pendekatan berbeda yaitu :
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
Penulis 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(seasonal = T)
Penulis tidak membuat model SARIMA dengan tuning manual karena dataset merupakan multi seasonal sehingga menurut penulis model SARIMA pada dasarnya tidak akan menghasilkan model yang optimal, sehingga penulis hanya menggunakan model Auto SARIMA
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)
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))
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 :
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 :
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 :
Penulis 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 :
#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 :
Insight Bisnis :
Penulis menggunakan metode RMSE atau Root Mean Squared Error untuk mengukur performa model. Untuk pengukuran hasil model, penulis 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.
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.
Penulis menggunakan model TBATS untuk cek asumsi, karena model tersebut yang performanya paling baik.
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 :
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.
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 :