Scotty: Bring Me The Crystal Ball!

Pengantar
- Tujuan
Memuat R-Packages yang Dibutuhkan
Persiapan Data
- Memuat Data
Analisis Data
Membentuk Data dalam Format Time Series
- Cross Validation
- Data Time Series
Pembuatan Model Menggunakan Data Train
Perbaikan Model
Melakukan Forecast Final untuk Tanggal 3 Desember 2017 sampai 9 Desember 2017
Pembuatan Model Menggunakan Data Train Digabungkan Dengan Data Test
- Pembuatan Model Berdasarkan Data Train dan Data Test
- Melakukan Forecast Final untuk Tanggal 3 Desember 2017 sampai 9 Desember 2017
Summary

Pengantar

Scotty adalah sebuah perusahaan ride-sharing yang beroperasi di beberapa kota besar di Turki sejak tahun 2017. Perusahaan tersebut menyediakan layanan ride-sharing untuk sepeda motor bagi penduduk Turki dan menjadi solusi yang efisien untuk berpergian di tengah kemacetan yang ada.

Scotty akan menyediakan dataset dari transaksinya dan dengan data tersebut, kita diharapkan bisa membantu mereka untuk melakukan forecast untuk meningkatkan bisnis mereka.

Tujuan

Di penghujung tahun 2017 ini, kita diharapkan bisa membantu Scotty agar Scotty siap untuk menangani permintaan di akhir tahun. Sayangnya, karena Scotty baru berdiri di tahun 2017, kita tidak memiliki data untuk Desember 2016 sehingga kita tidak bisa melihat permintaan di tahun sebelumnya. Untuk itu, kita akan menggunakan transaksi beberapa bulan terakhir untuk melakukan forecast permintaan untuk Desember 2017 per jam nya berdasarkan sub-area permintaan sehingga Scotty dapat memastikan pengemudi di sub-area tersebut dapat mencukupi permintaan yang ada. Selain itu kita juga perlu membuat automated forecasting framework sehingga kita tidak perlu mengulang tahapan-tahapan forecast yang sama di masa yang akan mendatang.

Memuat R-Packages yang Dibutuhkan

Untuk mulai melakukan membuat automated forecasting framework, pertama-tama kita harus memuat package/library R yang dibutuhkan:

library(readr)        #untuk membaca data di CSV, contoh fungsi read_csv()
library(DT)           #untuk membuat tampilan data berupa table, contoh fungsi datatables()
library(ggplot2)      #untuk membuat plot
library(dplyr)        #untuk proses praprocess data
library(lubridate)    #untuk melakukan manipulasi data yang berupa tanggal
library(padr)         #untuk melengkapi data yang tidak lengkap, contoh fungsi padding()
library(zoo)          #untuk mengisi data yang nilainya N/A
library(forecast)     #untuk fungsi msts(), ets()
library(tidyr)        #untuk fungsi spread(), nest()
library(recipes)      #untuk fungsi recipe(), bake()
library(tibble)       #untuk fungsi enframe()
library(purrr)        #untuk fungsi map()
library(tidymodels)   #untuk fungsi rmse_vec()
library(png)          #untuk fungsi readPNG()
library(grid)         #untuk fungsi grid.raster()

Persiapan Data

Memuat Data

Berikutnya yang harus kita lakukan adalah memuat data yang akan kita gunakan untuk pemodelan kita nantinya, data tersebut telah disiapkan dalam file CSV yang bernama data-train.csv dan setelah dimuat ke R kita, data tsb akan diberi nama scotty_data:

scotty_data <- read_csv("data/data-train.csv")

Selanjutnya kita akan melakukan pengecekan tipe data atas scotty_data:

str(scotty_data)

## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 90113 obs. of  16 variables:
##  $ id                : chr  "59d005e1ffcfa261708ce9cd" "59d0066affcfa261708ceb11" "59d006a1ffcfa261708ceba4" "59d006d8cb564761a8fe5efd" ...
##  $ trip_id           : chr  "59d005e9cb564761a8fe5d3e" NA "59d006c131e39c618969343d" "59d007193d32b861760d4a77" ...
##  $ driver_id         : chr  "59a892c5568be44b2734f276" NA "599dc0dfa5b4fd5471ad8453" "59a5856573d56e7103b5d51d" ...
##  $ rider_id          : chr  "59ad2d6efba75a581666b506" "59cd704bcf482f6ce2fadfdb" "59bd62cc7e3c3b663924ba86" "59c17cd1f6da2d274e16d0a7" ...
##  $ start_time        : POSIXct, format: "2017-10-01 00:00:17" "2017-10-01 00:02:34" ...
##  $ src_lat           : num  41.1 41.1 41 41.1 41.1 ...
##  $ src_lon           : num  29 29 29 29 29 ...
##  $ src_area          : chr  "sxk9" "sxk9" "sxk9" "sxk9" ...
##  $ src_sub_area      : chr  "sxk9s" "sxk9e" "sxk9s" "sxk9e" ...
##  $ dest_lat          : num  41.1 41.1 41 41.1 41 ...
##  $ dest_lon          : num  29 29 29 29 29.1 ...
##  $ dest_area         : chr  "sxk9" "sxk9" "sxk9" "sxk9" ...
##  $ dest_sub_area     : chr  "sxk9u" "sxk9e" "sxk9e" "sxk9e" ...
##  $ distance          : num  5.38 1.13 4.17 3.36 11.69 ...
##  $ status            : chr  "confirmed" "nodrivers" "confirmed" "confirmed" ...
##  $ confirmed_time_sec: num  8 0 32 65 0 27 32 30 0 24 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   id = col_character(),
##   ..   trip_id = col_character(),
##   ..   driver_id = col_character(),
##   ..   rider_id = col_character(),
##   ..   start_time = col_datetime(format = ""),
##   ..   src_lat = col_double(),
##   ..   src_lon = col_double(),
##   ..   src_area = col_character(),
##   ..   src_sub_area = col_character(),
##   ..   dest_lat = col_double(),
##   ..   dest_lon = col_double(),
##   ..   dest_area = col_character(),
##   ..   dest_sub_area = col_character(),
##   ..   distance = col_double(),
##   ..   status = col_character(),
##   ..   confirmed_time_sec = col_double()
##   .. )

Dari data di atas, kita dapat melihat bahwa scotty_data terdiri dari 90,113 baris dan 16 kolom / variabel.

Analisis Data

Penjelasan Variabel

Berikut merupakan penjelasan dari masing-masing variabel yang ada di scotty_data:

id : merupakan kode unik per baris data kita.
trip_id : merupakan ID perjalanan.
driver_id : merupakan ID pengemudi.
rider_id : merupakan ID penumpang.
start_time : merupakan waktu ketika seorang penumpang melakukan permintaan kendaraan.
src_lat : merupakan Latitude ketika penumpang melakukan permintaan kendaraan.
src_lon : merupakan Longitude ketika penumpang melakukan permintaan kendaraan.
src_area : src_lat dan src_lon dikelompokkan menjadi area permintaan.
src_sub_area : src_lat dan src_lon dikelompokkan menjadi sub-area permintaan.
dest_lat : merupakan Latitude dari lokasi tujuan penumpang.
dest_lon : merupakan Longitude dari lokasi tujuan penumpang.
dest_area : dest_lat dan dest_lon dikelompokkan menjadi area tujuan.
dest_sub_area : dest_lat dan dest_lon dikelompokkan menjadi sub-area tujuan.
distance : merupakan jarak perjalanan (dalam Kilometer).
status : merupakan status perjalanan, apakah permintaan perjalanan tsb berhasil mendapatkan pengemudi (confirmed) atau tidak (nodrivers). Tapi semua status tersebut tetap dianggap sebagai permintaan.
confirmed_time_sec : merupakan jarak waktu sejak permintaan dilakukan sampai permintaan tersebut dikonfirmasi oleh pengemudi (dalam Detik).

Pengolahan Data (Data Preprocessing)

Karena kita akan melakukan forecast berdasarkan seri waktu, maka dari 16 kolom / variabel yang ada di scotty_data, kita hanya akan menggunakan 2 variabel yaitu start_time dan src_sub_area, maka kita akan membuat semua variabel lainnya:

scotty_data <- scotty_data %>% 
  select(start_time, src_sub_area)

Variabel start_time merupakan waktu ketika seorang penumpang melakukan permintaan kendaraan, maka kita akan mengkategorikan jam permintaan dengan melakukan pembulatan ke jam terdekat dan kategori jam permintaan tsb akan ditampung di variable datetime. Misalkan permintaan dilakukan pukul 10.03, maka kita kategorikan permintaan tsb ke jam 10. Kita akan menggunakan fungsi floor_date() dari library lubridate.

scotty_data$datetime <- floor_date(scotty_data$start_time, unit = "hour")

Karena scotty_data berupa raw data, maka tahapan berikutnya adalah menghitung jumlah permintaan per area dan per jam :

scotty_data_sum <- scotty_data %>% 
  group_by(datetime, src_sub_area) %>% 
  summarise(
    demand = n()
  )

Karakteristik dari Data Time Series

1. Tidak diperkenankan ada period waktu yang tidak lengkap.

Kita akan cek range data dari variabel datetime terlebih dulu:

range(scotty_data_sum$datetime)

## [1] "2017-10-01 00:00:00 UTC" "2017-12-02 23:00:00 UTC"

Dari data di atas, diketahui range data kita adalah dari 1 Oktober 2017 pukul 00:00 sampai 2 Desember 2017 pukul 23:00.

Tahap berikutnya adalah kita memastikan jumlah data per src_sub_area:

scotty_data_sum %>% 
  group_by(src_sub_area) %>% 
  summarise(
    demand = n()
  )

## # A tibble: 3 x 2
##   src_sub_area demand
##   <chr>         <int>
## 1 sxk97          1439
## 2 sxk9e          1419
## 3 sxk9s          1367

Berdasarkan data di atas, diketahui bahwa jumlah data per src_sub_area tidak sama, artinya di data kita ada period waktu yang tidak lengkap. Untuk itu kita akan melengkapi period waktu yang tidak lengkap pada variabel datetime menggunakan fungsi pad() dari library padr:

scotty_data_sum <- scotty_data_sum %>% 
  pad('hour', 
      start_val = as.POSIXct(min(scotty_data_sum$datetime)), 
      end_val = as.POSIXct(max(scotty_data_sum$datetime)), 
      group = "src_sub_area") %>%
  mutate(
    src_sub_area = as.factor(src_sub_area)
  )

scotty_data_sum %>% 
  group_by(src_sub_area) %>% 
  summarise(
    demand = n()
  )

## # A tibble: 3 x 2
##   src_sub_area demand
##   <fct>         <int>
## 1 sxk97          1512
## 2 sxk9e          1512
## 3 sxk9s          1512

Setelah dilakukan fungsi pad(), jumlah data per src_sub_area sudah sama yaitu sebanyak 1,512 data per sub area.

2. Data harus urut berdasarkan periode waktu

Kita akan mengurutkan data kita berdasarkan sub area dan waktu:

scotty_data_sum <- scotty_data_sum %>% 
  arrange(src_sub_area, datetime)

3. Tidak diperkenankan ada nilai yang hilang (N/A)

Cek apakah di data scotty_data_sum ada nilai yang hilang (N/A):

scotty_data_sum %>% 
  filter(is.na(demand) == 1) %>% 
  group_by() %>% 
  summarise(
    na_rows = n()
  )

## # A tibble: 1 x 1
##   na_rows
##     <int>
## 1     311

Ternyata terdapat 311 data yang tidak memiliki nilai, maka kita akan mengisi nilai yang hilang dengan nilai 0 menggunakan fungsi na.fill0() dari library zoo:

scotty_data_sum <- na.fill0(scotty_data_sum, fill = 0)

scotty_data_sum %>% 
  filter(is.na(demand) == 1) %>% 
  group_by() %>% 
  summarise(
    na_rows = n()
  )

## # A tibble: 1 x 1
##   na_rows
##     <int>
## 1       0

Sudah tidak ada lagi data yang tidak memiliki nilai.

Membentuk Data dalam Format Time Series

Menampilkan contoh data dari data 1 bulan terakhir ke dalam plot berdasarkan sub area nya:

scotty_data_sum %>%
  filter(datetime > (max(scotty_data_sum$datetime) - hours(24 * 7 * 4))) %>%
  ggplot(aes(x = datetime, y = demand)) +
    geom_line() +
    labs(x = NULL, y = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_light()

Cross Validation

Sebelum kita membuat model, kita harus membagi data kita menjadi data train dan data test terlebih dulu. Data train akan kita gunakan melakukan training terhadap model kita sedangkan data test akan digunakan menjadi pembanding dan melihat apakah model kita overfit / underfit untuk melakukan prediksi terhadap data baru di luar data training nya. Kita akan menggunakan 7 hari terakhir dari data kita untuk menjadi data test dan sisanya akan kita gunakan menjadi data train.

Karena data kita berupa tanggal, maka kita akan mencari range awal dan range akhir baik untuk data test maupun data train. Cara yang akan kita lakukan adalah mencari jumlah data untuk data train dan data test lalu kemudian mencari range awal dan range akhir:

# jumlah data test dan data train
test_size <- 24 * 7
train_size <- nrow(scotty_data_sum)/3 - test_size

# range awal dan range akhir untuk masing2 data
test_end <- max(scotty_data_sum$datetime)
test_start <- test_end - hours(test_size) + hours(1)

train_end <- test_start - hours(1)
train_start <- train_end - hours(train_size) + hours(1)

# interval untuk masing2 data
interval_train <- interval(train_start, train_end)
interval_test <- interval(test_start, test_end)

Setelah dibagi menjadi data train dan data test, berikut merupakan visualisasinya:

# plot data train dan data test
scotty_data_sum %>%
  mutate(sample = case_when(
    datetime %within% interval_train ~ "train",
    datetime %within% interval_test ~ "test"
  )) %>%
  mutate(sample = factor(sample, levels = c("train", "test"))) %>%
  ggplot(aes(x = datetime, y = demand, colour = sample)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_light() +
    tidyquant::scale_colour_tq()

Data Time Series

Sebelum kita membentuk data Time Series, kita akan melakukan pengolahan data terlebih dulu menggunakan library recipes.

Karena library recipes bekerja dalam bentuk kolom, maka kita akan mengubah data kita menjadi tabular berdasarkan src_sub_area menggunakan fungsi spread() dari library tidyr:

scotty_data_sum <- scotty_data_sum  %>% 
  spread(src_sub_area, demand)

Berikutnya kita akan melakukan manipulasi data yaitu scale yang akan dilakukan menggunakan fungsi recipe() library recipes:

scotty_recipe <- recipe(formula = ~.,  
                        data = scotty_data_sum) %>% 
  #step
  step_scale(all_numeric()) %>% #melakukan scale
  #prep
  prep()

Selanjutnya kita akan mengimplementasikan scotty_recipe ke data kita menggunakan fungsi bake():

scotty_data_recipe <- bake(scotty_recipe, scotty_data_sum)

Karena kita telah melakukan manipulasi data, kita harus membuat fungsi untuk mengembalikan data kita ke nilai awal yang akan kita namakan fungsi scotty_recipe_revert:

# revert back function
scotty_recipe_revert <- function(vector, recipe, varname) {

  # store recipe values
  results <- recipe$steps[[1]]$sds[varname] * vector

  # add additional adjustment if necessary
  results <- round(results)

  # return the results
  results

}

Setelah kita selesai melakukan manipulasi data, kita akan mengembalikan data kita menjadi bentuk baris lagi menggunakan fungsi gather():

scotty_data_recipe <- scotty_data_recipe %>% 
  gather(src_sub_area, demand, -datetime)

Menambahkan variabel baru yaitu variabel sample yang berfungsi untuk membedakan data train dan data test:

scotty_data_recipe <- scotty_data_recipe %>%
  mutate(sample = case_when(
    datetime %within% interval_train ~ "train",
    datetime %within% interval_test ~ "test"
  ))

Kita akan membuat kumpulan data berupa list berdasarkan sub area nya. Selanjutnya list tsb akan dibagi menjadi train dan test:

scotty_data_nest <- scotty_data_recipe %>%
  group_by(src_sub_area, sample) %>%
  nest(.key = "data") %>%
  pivot_wider(names_from = sample, values_from = data)

Untuk membuat data kita menjadi bentuk time series kita membutuhkan fungsi ts() atau msts(). Jadi kita akan membuat list data fungsi bernama ts_function_list yang berisi kedua fungsi tsb, dimana untuk time series akan menggunakan seasonality sebagai berikut:

ts yang dipergunakan untuk single seasonality, akan menggunakan seasonality harian (frekuensi = 24)
msts yyang dipergunakan untuk multiple / complex seasonality, akan menggunakan seasonality harian (frekuensi = 24) dan mingguan (frekuensi = 24 * 7)

# list fungsi yang akan digunakan untuk membentuk time series
ts_function_list <- list(
  ts = function(x) ts(x$demand, frequency = 24),
  msts = function(x) msts(x$demand, seasonal.periods = c(24, 24 * 7))
)

# kombinasikan ts_function_list dengan masing-masing sub area
ts_function_list <- ts_function_list %>% 
  rep(length(unique(scotty_data_nest$src_sub_area))) %>%
  enframe("func_name", "func") %>%
  mutate(
    src_sub_area = sort(rep(unique(scotty_data_nest$src_sub_area), length(unique(.$func_name))))
  )

ts_function_list

## # A tibble: 6 x 3
##   func_name func   src_sub_area
##   <chr>     <list> <chr>       
## 1 ts        <fn>   sxk97       
## 2 msts      <fn>   sxk97       
## 3 ts        <fn>   sxk9e       
## 4 msts      <fn>   sxk9e       
## 5 ts        <fn>   sxk9s       
## 6 msts      <fn>   sxk9s

Tahap terakhir adalah menggabungkan ts_function_list ke data scotty_data_nest:

scotty_data_nest <- scotty_data_nest %>% 
  left_join(ts_function_list)

## # A tibble: 6 x 5
## # Groups:   src_sub_area [3]
##   src_sub_area          train           test func_name func  
##   <chr>        <list<df[,2]>> <list<df[,2]>> <chr>     <list>
## 1 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>  
## 2 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>  
## 3 sxk9e           [1,344 × 2]      [168 × 2] ts        <fn>  
## 4 sxk9e           [1,344 × 2]      [168 × 2] msts      <fn>  
## 5 sxk9s           [1,344 × 2]      [168 × 2] ts        <fn>  
## 6 sxk9s           [1,344 × 2]      [168 × 2] msts      <fn>

Pembuatan Model Menggunakan Data Train

Membuat list model time series

Untuk membuat model bagi data yang bentuknya time series, kita akan membuat list data model bernama model_function_list yang berisi fungsi-fungsi untuk melakukan pemodelan data time series:

# list model yang akan digunakan untuk pemodelan time series
model_function_list <- list(
  ets = function(x) ets(x),
  auto.arima = function(x) auto.arima(x),
  stlm = function(x) stlm(x),
  tbats = function(x) tbats(x, use.box.cox = FALSE)
)

# kombinasikan model_function_list dengan masing-masing sub area
model_function_list <- model_function_list %>% 
  rep(length(unique(scotty_data_nest$src_sub_area))) %>%
  enframe("model_name", "model") %>%
  mutate(src_sub_area =
    sort(rep(unique(scotty_data_nest$src_sub_area), length(unique(.$model_name))))
  )

model_function_list

## # A tibble: 12 x 3
##    model_name model  src_sub_area
##    <chr>      <list> <chr>       
##  1 ets        <fn>   sxk97       
##  2 auto.arima <fn>   sxk97       
##  3 stlm       <fn>   sxk97       
##  4 tbats      <fn>   sxk97       
##  5 ets        <fn>   sxk9e       
##  6 auto.arima <fn>   sxk9e       
##  7 stlm       <fn>   sxk9e       
##  8 tbats      <fn>   sxk9e       
##  9 ets        <fn>   sxk9s       
## 10 auto.arima <fn>   sxk9s       
## 11 stlm       <fn>   sxk9s       
## 12 tbats      <fn>   sxk9s

Tahap berikutnya adalah menggabungkan model_function_list ke data scotty_data_nest. Dikarenakan model ets dan auto.arima tidak bisa digunakan untuk tipe time series yang multiseasonal, maka kita tidak akan menggabungkan kedua model tersebut dengan tipe time series yang msts :

scotty_data_nest <- scotty_data_nest %>% 
  left_join(model_function_list) %>%
  filter(
    !(model_name == "ets" & func_name == "msts"),
    !(model_name == "auto.arima" & func_name == "msts")
  )

## # A tibble: 6 x 7
## # Groups:   src_sub_area [1]
##   src_sub_area          train           test func_name func   model_name model 
##   <chr>        <list<df[,2]>> <list<df[,2]>> <chr>     <list> <chr>      <list>
## 1 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   ets        <fn>  
## 2 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   auto.arima <fn>  
## 3 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   stlm       <fn>  
## 4 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   tbats      <fn>  
## 5 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>   stlm       <fn>  
## 6 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>   tbats      <fn>

Pembuatan Model

Menjalankan seluruh data yang ada di scotty_data_nest terhadap fungsi dan model yang ada di data tsb. Tahap pertama adalah menjalankan fungsi untuk membentuk data time series berdasarkan variabel func lalu membentuk model berdasarkan variabel model dan hasilnya akan kita tampung ke variabel fitted:

scotty_data_nest <- scotty_data_nest %>% 
  mutate(
    params = map(train, ~ list(x = .x)),
    data = invoke_map(func, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  ) %>%
  select(-data, -params)

## # A tibble: 6 x 8
## # Groups:   src_sub_area [1]
##   src_sub_area        train        test func_name func  model_name model fitted 
##   <chr>        <list<df[,2> <list<df[,> <chr>     <lis> <chr>      <lis> <list> 
## 1 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  ets        <fn>  <ets>  
## 2 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  auto.arima <fn>  <fr_AR…
## 3 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  stlm       <fn>  <stlm> 
## 4 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  tbats      <fn>  <tbats>
## 5 sxk97         [1,344 × 2]   [168 × 2] msts      <fn>  stlm       <fn>  <stlm> 
## 6 sxk97         [1,344 × 2]   [168 × 2] msts      <fn>  tbats      <fn>  <tbats>

Membandingkan dengan Data Test

Setelah kita selesai melakukan pemodelan, berikutnya kita akan melihat nilai error kita jika model tersebut diaplikasikan ke data test kita:

scotty_data_nest <- scotty_data_nest %>% 
  mutate(
    error = map(fitted, ~ forecast(.x, h = 24 * 7)) %>% 
            map2_dbl(test, ~ rmse_vec(truth = .y$demand, estimate = .x$mean))
  ) %>%
  arrange(src_sub_area, error)

Berikutnya kita akan menampilkan data forecast kita ke plot untuk membandingkan data hasil forecast dan data test kita.

scotty_test <- scotty_data_nest %>%
  mutate(
    forecast =  map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
                map2(test, ~ tibble(datetime = .y$datetime, demand = as.vector(.x$mean))),
    key = paste(func_name, model_name, sep = "-")
  )

## # A tibble: 6 x 11
## # Groups:   src_sub_area [1]
##   src_sub_area       train      test func_name func  model_name model fitted
##   <chr>        <list<df[,> <list<df> <chr>     <lis> <chr>      <lis> <list>
## 1 sxk97        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## 2 sxk97        [1,344 × 2] [168 × 2] msts      <fn>  stlm       <fn>  <stlm>
## 3 sxk97        [1,344 × 2] [168 × 2] ts        <fn>  tbats      <fn>  <tbat…
## 4 sxk97        [1,344 × 2] [168 × 2] ts        <fn>  stlm       <fn>  <stlm>
## 5 sxk97        [1,344 × 2] [168 × 2] ts        <fn>  ets        <fn>  <ets> 
## 6 sxk97        [1,344 × 2] [168 × 2] ts        <fn>  auto.arima <fn>  <fr_A…
## # … with 3 more variables: error <dbl>, forecast <list>, key <chr>

Selanjutnya kita akan mengembalikan data kita yang tadinya berupa kolom menjadi baris lagi:

scotty_test <- scotty_test %>% 
  select(src_sub_area, key, actual = test, forecast) %>%
  spread(key, forecast) %>%
  gather(key, value, -src_sub_area)

## # A tibble: 6 x 3
## # Groups:   src_sub_area [3]
##   src_sub_area key       value             
##   <chr>        <chr>     <list>            
## 1 sxk97        actual    <tibble [168 × 2]>
## 2 sxk9e        actual    <tibble [168 × 2]>
## 3 sxk9s        actual    <tibble [168 × 2]>
## 4 sxk97        msts-stlm <tibble [168 × 2]>
## 5 sxk9e        msts-stlm <tibble [168 × 2]>
## 6 sxk9s        msts-stlm <tibble [168 × 2]>

Dikarenakan data hasil forecast kita masih dalam hasil fungsi scotty_recipe, maka kita akan mengembalikan nilai asli nya menggunakan fungsi scotty_recipe_revert:

scotty_test <- scotty_test %>% 
  unnest(value) %>%
  mutate(demand = scotty_recipe_revert(demand, scotty_recipe, src_sub_area))

Menampilkan data hasil forecast dan data test dalam bentuk plot:

scotty_test %>%
  ggplot(aes(x = datetime, y = demand, colour = key)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_light() +
    tidyquant::scale_colour_tq()

Pemilihan Model Otomatis

Dari model-model yang telah dibuat sebelumnya, kita akan mengambil model yang nilai errornya paling kecil berdasarkan sub area dengan langkah sebagai berikut:

scotty_model_selected <- scotty_data_nest %>% 
  group_by(src_sub_area) %>%
  filter(error == min(error)) %>%
  ungroup()

## # A tibble: 3 x 9
##   src_sub_area       train      test func_name func  model_name model fitted
##   <chr>        <list<df[,> <list<df> <chr>     <lis> <chr>      <lis> <list>
## 1 sxk97        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## 2 sxk9e        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## 3 sxk9s        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## # … with 1 more variable: error <dbl>

Evaluasi model menggunakan MAE untuk data test:

scotty_mae <- scotty_test %>% 
  filter(
    paste(key, src_sub_area, sep = "-") %in% 
      paste(scotty_model_selected$func_name, 
            scotty_model_selected$model_name, 
            scotty_model_selected$src_sub_area, 
            sep = "-") 
    | key == "actual"
  ) %>% 
  mutate(
    demand = as.numeric(demand),
    key = if_else(key == "actual", "actual", "forecast")
  ) %>% 
  spread(key, demand) 

# MAE per sub area
scotty_mae_by_src_sub_area <- scotty_mae %>%
  group_by(src_sub_area) %>%
  nest() %>% 
  mutate(
    tmp = map(data, ~invoke_map_dfc(list(mae = mae_vec),
                                    truth = .x$actual, estimate = .x$forecast))
  ) %>%
  select(-data) %>%
  unnest() %>% 
  ungroup() %>% 
  add_row(src_sub_area = "all-sub-area", 
          mae = MLmetrics::MAE(scotty_mae$forecast, scotty_mae$actual)) %>% 
  mutate(
    mae = round(mae,4)
  )

Perbaikan Model

Kita akan mencoba untuk melakukan perbaikan model yaitu dengan melakukan improvement pada bagian pre-processing data. Jika sebelumnya kita akan melakukan pengolahan data menggunakan fungsi scotty_recipe dimana kita hanya melakukan scaling data, maka sekarang kita akan mencoba membuat fungsi baru yang bernama scotty_recipe_improve dengan step pre-processing :

melakukan akar
melakukan pengurangan dengan nilai tengah
melakukan scaling data

scotty_recipe_improve <- recipe(formula = ~.,  
                        data = scotty_data_sum) %>% 
  #step
  step_sqrt(all_numeric()) %>% #melakukan akar
  step_center(all_numeric()) %>% #melakukan pengurangan dengan nilai tengah
  step_scale(all_numeric()) %>% #melakukan scaling data
  #prep
  prep()

Selanjutnya kita akan mengimplementasikan scotty_recipe_improve ke data scotty_data_sum kita menggunakan fungsi bake():

scotty_data_recipe_improve <- bake(scotty_recipe_improve, scotty_data_sum)

Karena kita telah melakukan manipulasi data, kita harus membuat fungsi untuk mengembalikan data kita ke nilai awal yang akan kita namakan fungsi scotty_recipe_improve_revert:

# revert back function
scotty_recipe_improve_revert <- function(vector, recipe, varname) {

  # store recipe values
  recipe_center <- recipe$steps[[2]]$means[varname]
  recipe_scale <- recipe$steps[[3]]$sds[varname]

  # convert back based on the recipe
  results <- (vector * recipe_scale + recipe_center) ^ 2

  # add additional adjustment if necessary
  results <- round(results)

  # return the results
  results

}

Setelah kita itu kita akan mengembalikan data kita menjadi bentuk baris lagi menggunakan fungsi gather() dan menambahkan variabel baru yaitu variabel sample yang berfungsi untuk membedakan data train dan data test:

scotty_data_recipe_improve <- scotty_data_recipe_improve %>% 
  gather(src_sub_area, demand, -datetime) %>%
  mutate(sample = case_when(
    datetime %within% interval_train ~ "train",
    datetime %within% interval_test ~ "test"
  ))

Kita akan membuat kumpulan data berupa list berdasarkan sub area nya. Selanjutnya list tsb akan dibagi menjadi train dan test:

scotty_data_nest_improve <- scotty_data_recipe_improve %>%
  group_by(src_sub_area, sample) %>%
  nest(.key = "data") %>%
  pivot_wider(names_from = sample, values_from = data)

Berikutnya kita akan mengkombinasikan fungsi ts_function_list (berisi list fungsi untuk membentuk time series) yang telah kita buat sebelumnya dengan data scotty_data_nest_improve kita :

scotty_data_nest_improve <- scotty_data_nest_improve %>% 
  left_join(ts_function_list)

## # A tibble: 6 x 5
## # Groups:   src_sub_area [3]
##   src_sub_area          train           test func_name func  
##   <chr>        <list<df[,2]>> <list<df[,2]>> <chr>     <list>
## 1 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>  
## 2 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>  
## 3 sxk9e           [1,344 × 2]      [168 × 2] ts        <fn>  
## 4 sxk9e           [1,344 × 2]      [168 × 2] msts      <fn>  
## 5 sxk9s           [1,344 × 2]      [168 × 2] ts        <fn>  
## 6 sxk9s           [1,344 × 2]      [168 × 2] msts      <fn>

Tahap berikutnya adalah menggabungkan model_function_list (berisi list model) ke data scotty_data_nest_improve:

scotty_data_nest_improve <- scotty_data_nest_improve %>% 
  left_join(model_function_list) %>%
  filter(
    !(model_name == "ets" & func_name == "msts"),
    !(model_name == "auto.arima" & func_name == "msts")
  )

## # A tibble: 6 x 7
## # Groups:   src_sub_area [1]
##   src_sub_area          train           test func_name func   model_name model 
##   <chr>        <list<df[,2]>> <list<df[,2]>> <chr>     <list> <chr>      <list>
## 1 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   ets        <fn>  
## 2 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   auto.arima <fn>  
## 3 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   stlm       <fn>  
## 4 sxk97           [1,344 × 2]      [168 × 2] ts        <fn>   tbats      <fn>  
## 5 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>   stlm       <fn>  
## 6 sxk97           [1,344 × 2]      [168 × 2] msts      <fn>   tbats      <fn>

Selanjutnya kita akan menjalankan seluruh data yang ada di scotty_data_nest_improve terhadap fungsi dan model yang ada di data tsb dan hasilnya akan kita tampung ke variabel fitted:

scotty_data_nest_improve <- scotty_data_nest_improve %>% 
  mutate(
    params = map(train, ~ list(x = .x)),
    data = invoke_map(func, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  ) %>%
  select(-data, -params)

## # A tibble: 6 x 8
## # Groups:   src_sub_area [1]
##   src_sub_area        train        test func_name func  model_name model fitted 
##   <chr>        <list<df[,2> <list<df[,> <chr>     <lis> <chr>      <lis> <list> 
## 1 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  ets        <fn>  <ets>  
## 2 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  auto.arima <fn>  <fr_AR…
## 3 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  stlm       <fn>  <stlm> 
## 4 sxk97         [1,344 × 2]   [168 × 2] ts        <fn>  tbats      <fn>  <tbats>
## 5 sxk97         [1,344 × 2]   [168 × 2] msts      <fn>  stlm       <fn>  <stlm> 
## 6 sxk97         [1,344 × 2]   [168 × 2] msts      <fn>  tbats      <fn>  <tbats>

Dari model-model yang telah dibuat sebelumnya, kita akan mengambil model yang nilai errornya paling kecil berdasarkan sub area dengan langkah sebagai berikut:

scotty_model_improve_selected <- scotty_data_nest_improve %>%
  # check error ketika diaplikasikan ke data test
  mutate(
    error = map(fitted, ~ forecast(.x, h = 24 * 7)) %>% 
            map2_dbl(test, ~ rmse_vec(truth = .y$demand, estimate = .x$mean))
  ) %>%
  arrange(src_sub_area, error) %>%
  # ambil nilai eror terkecil
  group_by(src_sub_area) %>%
  filter(error == min(error)) %>%
  ungroup()

## # A tibble: 3 x 9
##   src_sub_area       train      test func_name func  model_name model fitted
##   <chr>        <list<df[,> <list<df> <chr>     <lis> <chr>      <lis> <list>
## 1 sxk97        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## 2 sxk9e        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## 3 sxk9s        [1,344 × 2] [168 × 2] msts      <fn>  tbats      <fn>  <tbat…
## # … with 1 more variable: error <dbl>

Jika diimplementasikan ke data test kita, maka nilai MAEnya adalah :

scotty_mae_improve <- scotty_model_improve_selected %>% 
  # hitung nilai forecast berdasarkan data test
  mutate(
    forecast =  map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
                map2(test, ~ tibble(datetime = .y$datetime, demand = as.vector(.x$mean))),
    key = paste(func_name, model_name, sep = "-")
  ) %>% 
  # mengembalikan data kolom menjadi baris
  select(src_sub_area, key, actual = test, forecast) %>%
  spread(key, forecast) %>%
  gather(key, value, -src_sub_area) %>% 
  # mengembalikan nilai asli menggunakan `scotty_recipe_improve_revert`
  unnest(value) %>%
  mutate(demand = scotty_recipe_improve_revert(demand, scotty_recipe_improve, src_sub_area)) %>% 
  # membandingkan dengan model yang nilai error nya paling kecil
  filter(
    paste(key, src_sub_area, sep = "-") %in% 
      paste(scotty_model_improve_selected$func_name, 
            scotty_model_improve_selected$model_name, 
            scotty_model_improve_selected$src_sub_area, 
            sep = "-") 
    | key == "actual"
  ) %>% 
  mutate(
    demand = as.numeric(demand),
    key = if_else(key == "actual", "actual", "forecast")
  ) %>% 
  spread(key, demand)

scotty_mae_improve_by_sub_area <- scotty_mae_improve %>% 
  # hitung MAE per sub area
  group_by(src_sub_area) %>%
  nest() %>% 
  mutate(
    tmp = map(data, ~invoke_map_dfc(list(mae_improve = mae_vec),
                                    truth = .x$actual, estimate = .x$forecast))
  ) %>%
  select(-data) %>%
  unnest() %>% 
  ungroup() %>% 
  add_row(src_sub_area = "all-sub-area", 
          mae_improve = MLmetrics::MAE(scotty_mae_improve$forecast, scotty_mae_improve$actual)) %>% 
  mutate(
    mae_improve = round(mae_improve,4)
  )

Kesimpulan : Ketika kita melakukan preprocessing data menggunakan fungsi step_sqrt(), step_center(), dan step_scale() ternyata kita menghasilkan model yang lebih baik. Hal ini ditunjukkan dengan kita menghasilkan Mean Average Error (MAE) terhadap data test yang lebih kecil.

Melakukan Forecast Final untuk Tanggal 3 Desember 2017 sampai 9 Desember 2017

Dikarenakan data train kita hanya sampai tanggal 25 November 2017, maka kita harus menambahkan 14 hari ke depan untuk melakukan forecast data:

scotty_forecast <- scotty_model_improve_selected %>%
  select(src_sub_area, train, everything(), -test) %>% 
  mutate(forecast = map(fitted, ~ forecast(.x, h = 24 * 7 * 2)) %>%
                    map2(train, ~ tibble(datetime = timetk::tk_make_future_timeseries(.y$datetime, 24 * 7 * 2),
                                         demand = as.vector(.x$mean)))
  )

## # A tibble: 3 x 9
##   src_sub_area       train func_name func  model_name model fitted error
##   <chr>        <list<df[,> <chr>     <lis> <chr>      <lis> <list> <dbl>
## 1 sxk97        [1,344 × 2] msts      <fn>  tbats      <fn>  <tbat… 0.616
## 2 sxk9e        [1,344 × 2] msts      <fn>  tbats      <fn>  <tbat… 0.450
## 3 sxk9s        [1,344 × 2] msts      <fn>  tbats      <fn>  <tbat… 0.511
## # … with 1 more variable: forecast <list>

Mengembalikan nilai hasil forecast ke nilai asli nya menggunakan fungsi scotty_recipe_improve_revert:

scotty_forecast <- scotty_forecast %>% 
  select(src_sub_area, actual = train, forecast) %>%
  gather(key, value, -src_sub_area) %>%
  unnest(value) %>%
  mutate(demand = scotty_recipe_improve_revert(demand, scotty_recipe_improve, src_sub_area))

Melihat hasil forecast dalam bentuk plot:

scotty_forecast %>%
  ggplot(aes(x = datetime, y = demand, colour = key)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_light() +
    tidyquant::scale_colour_tq()

Export data hasil forecast kita ke dalam bentuk file CSV hanya untuk data di range tanggal 3 Desember 2017 sampai 9 Desember 2017 :

# jumlah data forecast
forecast_size <- 24 * 7

# range awal dan range akhir untuk data forecast
forecast_end <- max(scotty_forecast$datetime)
forecast_start <- forecast_end - hours(forecast_size) + hours(1)

# interval untuk data forecast
inforecast <- interval(forecast_start, forecast_end)

# filter data yang ada di interval data forecast saja
scotty_forecast <- scotty_forecast %>%
  mutate(type = case_when(
    datetime %within% inforecast ~ "forecast"
  )) %>%
  filter(
    type == 'forecast'
  ) %>% 
  select(-c(key, type))

# export ke file CSV
write.csv(scotty_forecast,file="result/submission-train-data.csv", row.names = FALSE)

Pembuatan Model Menggunakan Data Train Digabungkan Dengan Data Test

Menggabungkan data train dan data test:

scotty_combine <- scotty_model_improve_selected %>% 
  mutate(fulldata = map2(train, test, ~ bind_rows(.x, .y))) %>%
  select(src_sub_area, fulldata, everything(), -train, -test, -error) %>% 
  select(-fitted)

## # A tibble: 3 x 6
##   src_sub_area fulldata             func_name func   model_name model 
##   <chr>        <list>               <chr>     <list> <chr>      <list>
## 1 sxk97        <tibble [1,512 × 2]> msts      <fn>   tbats      <fn>  
## 2 sxk9e        <tibble [1,512 × 2]> msts      <fn>   tbats      <fn>  
## 3 sxk9s        <tibble [1,512 × 2]> msts      <fn>   tbats      <fn>

Pembuatan Model Berdasarkan Data Train dan Data Test

Menjalankan seluruh data yang ada di scotty_combine terhadap fungsi dan model yang ada di data tsb. Tahap pertama adalah menjalankan fungsi untuk membentuk data time series berdasarkan variabel func lalu membentuk model berdasarkan variabel model dan hasilnya akan kita tampung ke variabel fitted:

scotty_combine <- scotty_combine %>% 
  mutate(
    params = map(fulldata, ~ list(x = .x)),
    data = invoke_map(func, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  )

Melakukan Forecast Final untuk Tanggal 3 Desember 2017 sampai 9 Desember 2017

Dikarenakan data train dan data test kita sudah digabungkan, maka total data kita adalah sampai tanggal 2 Desember 2017, maka kita harus menambahkan 7 hari ke depan untuk melakukan forecast data:

scotty_forecast_full <- scotty_combine %>% 
  mutate(forecast =  map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
                     map2(fulldata, ~ tibble(datetime = timetk::tk_make_future_timeseries(.y$datetime, 24 * 7),
                                             demand = as.vector(.x$mean)))
  )

Mengembalikan nilai hasil forecast ke nilai asli nya menggunakan fungsi scotty_recipe_improve_revert:

scotty_forecast_full <- scotty_forecast_full %>% 
  select(src_sub_area, actual = fulldata, forecast) %>%
  gather(key, value, -src_sub_area) %>%
  unnest(value) %>%
  mutate(demand = scotty_recipe_improve_revert(demand, scotty_recipe_improve, src_sub_area))

Melihat hasil forecast dalam bentuk plot:

scotty_forecast_full %>%
  ggplot(aes(x = datetime, y = demand, colour = key)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    theme_light() +
    tidyquant::scale_colour_tq()

Export data hasil forecast kita ke dalam bentuk file CSV hanya untuk data di range tanggal 3 Desember 2017 sampai 9 Desember 2017 :

# filter data yang ada tipe forecast saja
scotty_forecast_full <- scotty_forecast_full %>% 
  filter(key=="forecast") %>% 
  select(-key)

# export ke file CSV
write.csv(scotty_forecast_full,file="result/submission-combine-data.csv", row.names = FALSE)

Summary

Pre-processing data menggunakan fungsi step_sqrt(), step_center(), dan step_scale() pada recipe scotty_recipe_improve ternyata menghasilkan nilai yang lebih baik daripada pre-processing data yang hanya menggunakan step_scale() pada recipe scotty_recipe dengan perbandingan MAE sebagai berikut:
Berdasarkan pemilihan model otomatis, data time series yang menggunakan pola Complex Seasonality (harian dan mingguan) menghasilkan nilai error yang lebih kecil dibandingkan dengan data time series yang menggunakan pola Single Seasonality (harian).
Setelah kedua data submission dikirimkan melalui Scoring Dashboard, hasil yang dicapai oleh hasil forecast untuk model yang dibuat dari data train saja adalah :

Sedangkan hasil yang dicapai oleh hasil forecast untuk model yang dibuat dari data train dan data test adalah :

Jadi dapat disimpulkan, bahwa model yang dibuat dari data train saja menghasilkan nilai forecast yang lebih baik.