INTRODUCTION

    Semua makhluk hidup di dunia ini pasti membutuhkan makanan. Dan makan itu sendiri banyak ragam dan variasinya. Dan salah satu makanan yang biasa dikonsumsi adalah roti. Roti bisa di jadikan sebagai sarapan, jamuan untuk tamu, atau mungkin sebagai cemilan sebagai pengganjal lapar.
    Pada data kali ini akan mengolah data dari sebuah toko roti "The Bread Basket" yang berlokasi di Edinburgh. Kumpulan data tersebut memberikan informasi mengenai rincian transaksi pembelian dari toko roti "The Bread Basket" yang memesan item berbeda dari toko roti ini secara online selama periode waktu 2016-01-11 hingga 2017-12-03. Adapun tujuan dari analisis ini yaitu mengetahui pola pembelian dari customer atau pelanggan pada sebuah toko roti.

link : https://www.kaggle.com/datasets/akashdeepkuila/bakery

IMPORT LIBRARY

library(dplyr)
library(lubridate) 
library(tidyr)
library(ggplot2)
library(TSstudio) 
library(forecast)
library(readr)
library(padr)
library(zoo)
library(ggsci)

READ DATA

Hal pertama yang harus kita lakukan yaitu import data kita terlebih dahulu read.csv.

data_bakerys <- read.csv("bakerys.csv")

DATA INSPECTIONS

Kita gunakan head() untuk melihat urutan data teratas dan gunakan tail() untuk melihat urutan data paling bawah.

head(data_bakerys)
tail(data_bakerys)

Cek struktur tipe data yang kita miliki menggunakan glimpse()

glimpse(data_bakerys)
#> Rows: 20,507
#> Columns: 5
#> $ TransactionNo <int> 1, 2, 2, 3, 3, 3, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8,…
#> $ Items         <chr> "Bread", "Scandinavian", "Scandinavian", "Hot chocolate"…
#> $ DateTime      <chr> "2016-10-30 09:58:11", "2016-10-30 10:05:34", "2016-10-3…
#> $ Daypart       <chr> "Morning", "Morning", "Morning", "Morning", "Morning", "…
#> $ DayType       <chr> "Weekend", "Weekend", "Weekend", "Weekend", "Weekend", "…

Dari fungsi glimps() yang telah kita jalankan diatas, kita mendapatkan informasi bahwa terdapat 20.507 rows dan 5 columns. Memuat informasi columns sebagai berikut :

`TransactionNo`    : Identifikasi nilai unik untuk setiap transaksi.
`Items`            : Barang yang dibeli.
`DateTime`         : Stempel tanggal dan waktu transaksi.
`Daypart`          : Bagian dari hari ketika transaksi dilakukan.
`DayType`          : Mengklasifikasikan apakah suatu transaksi telah dilakukan pada akhir pekan atau hari kerja.

CLEANSING DATA

Kita ubah tipe data yang kita miliki dengan tipe data yang sesuai menggunakan fungsi mmutate(). Dan memilih columns yang kita ingin kan menggunakan fungsi select().

bakerys_Z <- data_bakerys %>% 
  mutate(DateTime = ymd_hms(DateTime),
         ds = floor_date(DateTime, unit = "hours")) %>%
  select(DateTime, TransactionNo, ds)
# Aggregate/summarise
bakerys_Z <- bakerys_Z %>%
  group_by(ds) %>%
  summarise(customer_number = n_distinct(TransactionNo))

Time Series Padding, untuk mengisi rangkaian data yang hilang (dalam satuan jam), sehingga setiap hari data rangkaian waktu akan terdiri dari 24 record - 24 jam / per hari.

bakerys_Z <- bakerys_Z %>% pad()
# Replace NA Value
bakerys_Z <- bakerys_Z %>%
  mutate(customer_number = na.fill(customer_number, 0))

Dengan menggunakan fungsi anyNA() atau colsums(is.na()) yaitu untuk mengecheck apakah terdapat missing value atau tidak.

anyNA(bakerys_Z)
#> [1] FALSE
colSums(is.na(data_bakerys))
#> TransactionNo         Items      DateTime       Daypart       DayType 
#>             0             0             0             0             0

Filter the datetime dari 07.00 AM ke 16.00 PM.

bakerys_Z <- bakerys_Z %>% filter(hour(ds) >= 7 & hour(ds) <= 16)
bakerys_Z

Periksa mulai & akhir interval waktu setelah pemfilteran menggunakan fungsi range() dan min() & max().

#range
range(bakerys_Z$ds)
#> [1] "2016-01-11 07:00:00 UTC" "2017-12-03 16:00:00 UTC"
#min & max
min_date <- min(bakerys_Z$ds)
max_date <- max(bakerys_Z$ds)

min_date
#> [1] "2016-01-11 07:00:00 UTC"
max_date
#> [1] "2017-12-03 16:00:00 UTC"

SEASONALITY ANALYSIS

Observation & Decomposition

Buat Objek Time Series menggunakan fungsi ts() dengan frekuensi 10, dimana 10 adalah rentang waktu buka outlet setiap hari, mulai pukul 07.00 hingga 16.00.

bakerys_ts <- ts(data = bakerys_Z$customer_number,
                 frequency = 10)

Decompose and Plotting untuk memeriksa Seasonal, Trend and Irregular component of time series.

bakerys_ts_dc <- bakerys_ts %>%
  decompose(type = "additive")

bakerys_ts %>% tail(10*7*4) %>%
  decompose() %>% 
  autoplot()

Dari plot di atas, komponen Estimasi Tren menunjukkan sebuah pola, yang mungkin berisi ekstra musiman yang tidak tertangkap dan ini dapat dianggap sebagai data multimusim. Untuk mengatasi multi-musiman, kita perlu mengkonversi data menjadi “Multi-seasonality Time Series” yang menerima beberapa pengaturan frekuensi.

Kemudian Buat Objek “Multi-seasonality Time Series” menggunakan fungsi msts(), dengan Frequency : Harian (10) dan Mingguan (10*7), ini akan menangkap musiman dalam Harian dan Mingguan. Kemudian dekomposisi dan plot.

bakerys_msts <- msts(data = bakerys_Z$customer_number, seasonal.periods = c(10,10*7))

bakerys_msts_dc <- mstl(bakerys_msts)

bakerys_msts %>% tail(10*10*7) %>% stl(s.window = "periodic") %>% autoplot()

Berdasarkan Plot di atas, kini Estimasi Tren dengan “Multiple Seasonality Time Series” menjadi lebih halus dan jelas. Dan juga lebih jelas mengenai Musiman Harian dan Mingguan, lebih mudah dijelaskan untuk analisa lebih lanjut.

Membuat seasonality by hour.

bakerys_Z %>% 
  mutate(Hour = hour(ds), Seasonal = bakerys_ts_dc$seasonal) %>% 
  distinct(Hour, Seasonal) %>% 
  ggplot() +
  geom_bar(aes(x = Hour, y = Seasonal, fill = Seasonal), stat ="identity", width = 0.7)+
  scale_fill_gradient(low = "yellow", high = "red") +
  scale_x_continuous(breaks = seq(0,22,1)) +
  labs(title = "Seasonality Analysis  - Hourly")

Berdasarkan plot di atas, rentang Peak Time Outlet adalah antara pukul 07.00 - 16.00, dan pada pukul 11.00 merupakan waktu dimana Pengunjung paling banyak pengunjung datang. Hal tersebut terjadi dikarenakan pada pukul 11.00 sudah memasuki waktu makan siang. Dan pukul 07.00 merupakan waktu yang paling sedikit pengunjung yang datang ke Outlet karena Outlet baru dibuka.

Membuat Seasonality by Week & Hour

bakerys_msts_dc_fr <- data.frame(bakerys_msts_dc)

bakerys_msts_dc_fr %>%
  mutate(ds = bakerys_Z$ds) %>% 
  mutate(Day_of_Week  = wday(ds, label = TRUE, abbr = FALSE), Hour = (hour(ds))) %>% 
  group_by(Day_of_Week, Hour) %>%
  summarise(Seasonal = sum(Seasonal10 + Seasonal70)) %>%
  ggplot() +
  geom_bar(aes(x = Hour, y = Seasonal, fill = Day_of_Week),stat ="identity", position = "stack", width = 0.7)+
  scale_x_continuous(breaks = seq(0,22,1)) +
  scale_fill_locuszoom()+
  labs(title = "Multi-Seasonality Analysis  - Weekly & Hourly")

Berdasarkan plot di atas, rentang Peak Time Outlet adalah antara pukul 10.00 - 14.00 setiap hari dalam seminggu, dan pada hari Sabtu pukul 11.00 adalah waktu dimana Pengunjung paling banyak datang. Dan pukul 07.00 - setiap hari dalam seminggu, merupakan waktu dimana pengunjung paling sedikit datang karena Outlet baru saja dibuka.

Model Fitting and Evaluation

#Cross Validation
test_bakerys_msts  <- tail(bakerys_msts, 10*7)
train_bakerys_msts <- head(bakerys_msts, length(bakerys_msts)-length(test_bakerys_msts))

Modelling

#Holt’s Winters Exponential
model_hw_msts <- HoltWinters(train_bakerys_msts)
#ARIMA Model
model_arima_msts <- stlm(train_bakerys_msts, method = "arima")

Forecasting

Memperkirakan kedua model dalam 1 minggu ke depan.

hw_forecast_msts <- forecast(model_hw_msts, h = 10*7)
arima_forecast_msts <- forecast(model_arima_msts, h=10*7)

Vizualize Each Models

Visualisasi Holt Winter Model - Actual vs Estimated

bakerys_msts %>% autoplot(series = "Actual") +
  autolayer(hw_forecast_msts$mean, series = "Predict HW")

Visualisasi Arima Model - Actual vs Estimated

bakerys_msts %>% autoplot(series = "Actual") +
  autolayer(arima_forecast_msts$mean, series = "Predict Arima")

Summary Models

Berikut beberapa model yang dibuat dalam Multi-Seasonal Time Series

#Accuracy
acc_hw_msts <- accuracy(object = hw_forecast_msts$mean, test_bakerys_msts)
acc_arima_msts <- accuracy(object = arima_forecast_msts$mean, test_bakerys_msts)

result <- rbind(acc_hw_msts, acc_arima_msts)
rownames(result) <- c("Holt Winter Model", "Arima Model")
result
#>                         ME     RMSE      MAE MPE MAPE      ACF1 Theil's U
#> Holt Winter Model 2.184262 4.498085 2.751223 NaN  Inf 0.7372316       NaN
#> Arima Model       2.000002 4.376451 2.279497 NaN  Inf 0.7200583       NaN

Berdasarkan Rangkuman Akurasi di atas, Akurasi “ARIMA Model” lebih baik, MAE 2.279497 (lebih rendah dari 6), model ini akan dipilih untuk digunakan dalam Prediksi ==> Model Terbaik “ARIMA Model”.

Visualization Actual vs Estimated

#presiapan pembuatan plot
accuracyData <- data.frame(ds= bakerys_Z$ds %>% tail(10*7),
  actual = as.vector(test_bakerys_msts) ,
  HWForecastmsts = as.vector(hw_forecast_msts$mean) ,
  ArimaForecastmsts = as.vector(arima_forecast_msts$mean))
#Visualization of Actual vs Estimated number of visitors (Best Model)
accuracyData %>% 
 ggplot() +
  geom_line(aes(x = ds, y = actual, colour = "Actual"),size=1)+
  geom_line(aes(x = ds, y = ArimaForecastmsts, colour = "Arima Model"),size=1)+
  labs(title = "Hourly Visitor - Actual Vs Arima Model",x = "Date",y = "Customer Number",colour = "")

#Visualization of Actual vs Estimated number of visitors (All Models)
accuracyData %>% 
 ggplot() +
  geom_line(aes(x = ds, y = actual, colour = "Actual"),size=0.5)+
  geom_line(aes(x = ds, y = HWForecastmsts, colour = "Holt Winter Model"),size=0.1)+
  geom_line(aes(x = ds, y = ArimaForecastmsts, colour = "Arima Model"),size=0.5)+
  labs(title = "Hourly Visitor - Actual Vs All Models",x = "Date",y = "Visitor",colour = "")

Prediction Performance

MAE < 6 in validation dataset

accuracy(object = arima_forecast_msts$mean, test_bakerys_msts)
#>                ME     RMSE      MAE MPE MAPE      ACF1 Theil's U
#> Test set 2.000002 4.376451 2.279497 NaN  Inf 0.7200583       NaN

Model terbaik adalah Model Arima karena memiliki akurasi MAE sebesar 2.279497.

MAE < 6 in test dataset

Prediksi menggunakan Multi-Seasonal Data

#membaca test Test Data
test_data <- read.csv("bakerys.csv")

#Membuat ARIMA model dengan Multi-Seasonal Time Series Data 
model_arima_test <- stlm(bakerys_msts, method = "arima")

#Forecast model 1 minggu kedepan
forecast_arima_test <- forecast(model_arima_test, h=10*7)

Conclusion

Autocorrelation (Asumsi Autocorrelation menggunakan Ljung-box)

Box.test(model_arima_msts$residuals)
#> 
#>  Box-Pierce test
#> 
#> data:  model_arima_msts$residuals
#> X-squared = 10.154, df = 1, p-value = 0.00144

p-value = 0.00144 < 0.05 (alpha), itu artinya residu terdapat autocorrelation.

Normality of residuals (Normality of residuals diuji menngunakan saphiro.test)

train_bakerys_msts_shapiro <- head(bakerys_msts, 5000)

model_arima_msts_shapiro <- stlm(train_bakerys_msts_shapiro, method = "arima")

shapiro.test(model_arima_msts_shapiro$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_arima_msts_shapiro$residuals
#> W = 0.84063, p-value < 0.00000000000000022

p-value = 0.00000000000000022 < 0.05 (alpha), yang dapat disimpulkan bahwa residu terdapat autocorrelation.

SUMMARY

Dari hasil diatas didapatkan informasi bahwa Peak Time Outlet adalah antara pukul 10.00 - 14.00 setiap hari dalam seminggu, dan pada hari Sabtu pukul 11.00 adalah waktu dimana Pengunjung paling banyak datang. Dan pukul 07.00 - setiap hari dalam seminggu, merupakan waktu dimana pengunjung paling sedikit datang karena Outlet baru saja dibuka.