library(tidyverse)
library(ggplot2)
library(dplyr)
library(forecast)
library(purrr)
library(yardstick)
library(lubridate)
library(recipes)
library(magrittr)
library(plotly)
library(MLmetrics)

1 Preparasi Data

# Memasukkan data
scot <- read_csv("data_input/scotty-ts-auto/data/data-train.csv")

## Parsed with column specification:
## 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()
## )

summary(scot)

##       id              trip_id           driver_id        
##  Length:90113       Length:90113       Length:90113      
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##    rider_id           start_time                     src_lat     
##  Length:90113       Min.   :2017-10-01 00:00:17   Min.   :41.00  
##  Class :character   1st Qu.:2017-10-19 19:25:25   1st Qu.:41.04  
##  Mode  :character   Median :2017-11-07 14:39:07   Median :41.06  
##                     Mean   :2017-11-04 16:37:48   Mean   :41.06  
##                     3rd Qu.:2017-11-18 22:40:43   3rd Qu.:41.07  
##                     Max.   :2017-12-02 23:59:20   Max.   :41.09  
##     src_lon        src_area         src_sub_area          dest_lat    
##  Min.   :28.96   Length:90113       Length:90113       Min.   :11.13  
##  1st Qu.:28.98   Class :character   Class :character   1st Qu.:41.04  
##  Median :28.99   Mode  :character   Mode  :character   Median :41.05  
##  Mean   :29.00                                         Mean   :41.05  
##  3rd Qu.:29.01                                         3rd Qu.:41.07  
##  Max.   :29.05                                         Max.   :52.37  
##     dest_lon       dest_area         dest_sub_area         distance       
##  Min.   : 2.407   Length:90113       Length:90113       Min.   :   0.000  
##  1st Qu.:28.974   Class :character   Class :character   1st Qu.:   2.313  
##  Median :28.994   Mode  :character   Mode  :character   Median :   3.839  
##  Mean   :28.993                                         Mean   :   5.695  
##  3rd Qu.:29.017                                         3rd Qu.:   6.453  
##  Max.   :80.192                                         Max.   :5860.002  
##     status          confirmed_time_sec
##  Length:90113       Min.   :  0.00    
##  Class :character   1st Qu.: 12.00    
##  Mode  :character   Median : 27.00    
##                     Mean   : 47.04    
##                     3rd Qu.: 67.00    
##                     Max.   :495.00

Melihat summary data di atas, kolom yang akan digunakan dalam analisis adalah kolom start_time (waktu order) dan count dari row yang akan dihitung sebagai data permintaan. Data permintaan tersebut akan di grup untuk setiap jamnya menjadi data permintaan setiap jam.

Selanjutnya, melengkapi data jam agar tersedia untuk setiap jamnya walaupun tidak ada permintaan pada jam tersebut. Hal ini dilakukan agar saat modeling, model akan akurat dalam memprediksi kondisi waktu tanpa permintaan tersebut.

# Membulatkan satuan data waktu ke jam
scot %<>% 
  mutate(datetime = floor_date(start_time,"hour")) %>% 
  group_by(src_sub_area,datetime) %>% 
  summarise(demand = length(status))

# Melengkapi data observasi jam yang kosong
scot %<>% 
  complete(datetime = seq(min(scot$datetime), max(scot$datetime), by = "hour")) %>% 
  mutate(demand = ifelse(is.na(demand), 0, demand))

head(scot)

## # A tibble: 6 x 3
## # Groups:   src_sub_area [1]
##   src_sub_area datetime            demand
##   <chr>        <dttm>               <dbl>
## 1 sxk97        2017-10-01 00:00:00      6
## 2 sxk97        2017-10-01 01:00:00      4
## 3 sxk97        2017-10-01 02:00:00      9
## 4 sxk97        2017-10-01 03:00:00      2
## 5 sxk97        2017-10-01 04:00:00      5
## 6 sxk97        2017-10-01 05:00:00      4

Visualisasi data permintaan per sub area dapat dilihat pada plot di bawah.

ggplotly(ggplot(scot,aes(x = datetime, y = demand)) +
           geom_line(aes(col = src_sub_area)) +
           labs(x = "", y = "Permintaan (order)",title = "PERMINTAAN UNTUK SETIAP SUB-AREA") +
           facet_wrap(~ src_sub_area, scale = "free_y", ncol = 1) + 
           scale_color_brewer(palette = "Pastel2") +
           theme_test(base_family = "Bebas Neue") +
           theme(legend.position = "none"))

Membuat visualisasi data dari salah satu sub-area untuk mengidentifikasi pola musiman yang muncul. Saya menggunakan skala kuadrat agar dapat memperjelas pola yang terjadi.

sxk97 <- scot %>% filter(src_sub_area == "sxk97") %>% .$demand

ggseasonplot(ts(sxk97,frequency = 24),polar = T) +
  theme_test(base_family = "Bebas Neue") +
  theme(legend.position = "none") +
  labs(title = "Pola Harian",x = "Jam") +
  scale_y_sqrt()

ggseasonplot(ts(sxk97,frequency = 24*7),polar = T) +
  theme_test(base_family = "Bebas Neue") +
  theme(legend.position = "none") +
  labs(title = "Pola Mingguan",x =NULL) +
  scale_y_sqrt() 

  # scale_x_time(limits = c(0, 360),breaks = seq(1, 359.99, by = 51.42), labels=c("Minggu","Senin","Selasa","Rabu","Kamis","Jumat","Sabtu")) 
  
ggseasonplot(ts(sxk97,frequency = 24*7*4),polar = T) +
  theme_test(base_family = "Bebas Neue") +
  theme(legend.position = "none") +
  labs(title = "Pola Bulanan",x = "") +
  scale_y_sqrt()

  # scale_x_time(limits = c(0, 360),breaks = seq(1, 359.99, by = 90), labels=c("Minggu 1","Minggu 2","Minggu 3","Minggu 4")) 

Pada plot pola harian, terlihat permintaan tertinggi ada pada jam 17-19 malam. Sedangkan permintaan terendah ada di jam 5-6 pagi. Pada plot pola mingguan (dimulai dari hari minggu pagi di arah jam 12), dapat kita lihat bahwa permintaan tertinggi muncul di hari jumat atau sabtu. Pada plot pola bulanan (dimulai tanggal 1 Oktober di arah jam 12), dapat dilihat bahwa permintaan tinggi di minggu awal bulan dan akhir bulan. Setelah melihat plot dari pola harian, mingguan, dan bulanan, saya tentukan untuk menggunakan pola harian dan mingguan saja dalam pemodelan ini. Hal ini dikarenakan pola bulanan tidak terlalu terlihat karena observasi pada dataset hanya selama 2 bulan.

Setelah menentukan pola yang ingin digunakan, saya membuat visualisasi data setelah di dekomposisi oleh pola harian dan mingguan.

dcm <- scot %>% filter(src_sub_area == "sxk97") %>% .$demand %>% msts(.,seasonal.periods = c(24,24*7))
autoplot(mstl(dcm)) + theme_test(base_family = "Bebas Neue") + labs(title = "Dekomposisi data oleh pola mingguan dan pola harian")

Trend sudah terlihat halus, hal ini mengindikasikan pola musiman yang digunakan sudah cukup tepat.

Selanjutnya, melakukan Cross Validation dengan cara membagi data menjadi data Train (untuk melatih model) dan data Test (untuk mengevaluasi model). Data Test diperoleh dari observasi selama 1 minggu terakhir dan sisanya dimasukkan ke data Train.

# Besaran waktu data test
test_size <- 24 * 7

# Menentukan awal dan akhir dari data train dan test
test_end <- max(scot$datetime)
test_start <- test_end - hours(test_size) + hours(1)

train_end <- test_start - hours(1)
train_start <- min(scot$datetime)

intrain <- interval(train_start, train_end)
intest <- interval(test_start, test_end)

# Aplikasikan interval ke dataset
scot %<>%
  mutate(sample = case_when(
    datetime %within% intrain ~ "train",
    datetime %within% intest ~ "test"
  )) %>% 
  drop_na()

head(scot)

## # A tibble: 6 x 4
## # Groups:   src_sub_area [1]
##   src_sub_area datetime            demand sample
##   <chr>        <dttm>               <dbl> <chr> 
## 1 sxk97        2017-10-01 00:00:00      6 train 
## 2 sxk97        2017-10-01 01:00:00      4 train 
## 3 sxk97        2017-10-01 02:00:00      9 train 
## 4 sxk97        2017-10-01 03:00:00      2 train 
## 5 sxk97        2017-10-01 04:00:00      5 train 
## 6 sxk97        2017-10-01 05:00:00      4 train

Melakukan scaling data sebelum memulai pemodelan. Hal ini dilakukan agar pemodelan menjadi lebih robust dan lebih tidak sensitif pada data outlier.

# Mengubah data menjadi format lebar 
scot %<>%
  spread(src_sub_area, demand)

# Jenis scaling yang dilakukan adalah square root, center, scale
rec <- recipe(~ ., filter(scot)) %>%
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  prep()

# Menjalankan fungsi scaling
scot <- bake(rec, scot)

# Mengubah kembali data menjadi format tinggi
scot %<>%
  gather(src_sub_area, demand, -datetime,-sample)

head(scot)

## # A tibble: 6 x 4
##   datetime            sample src_sub_area demand
##   <dttm>              <fct>  <chr>         <dbl>
## 1 2017-10-01 00:00:00 train  sxk97        -0.612
## 2 2017-10-01 01:00:00 train  sxk97        -0.855
## 3 2017-10-01 02:00:00 train  sxk97        -0.313
## 4 2017-10-01 03:00:00 train  sxk97        -1.17 
## 5 2017-10-01 04:00:00 train  sxk97        -0.727
## 6 2017-10-01 05:00:00 train  sxk97        -0.855

Setelah melakukan scaling, tidak lupa juga membuat fungsi untuk mengembalikan data menjadi nilai dengan skala sebenarnya.

# Fungsi untuk mengembalikan data menjadi ukuran skala sebenarnya
rec_revert <- function(vector, rec, varname) {
  rec_center <- rec$steps[[2]]$means[varname]
  rec_scale <- rec$steps[[3]]$sds[varname]
  results <- (vector * rec_scale + rec_center) ^ 2
  results <- round(results)
  results
}

Selanjutnya, melakukan nest kepada data Training dan data Test agar mempermudah proses pemilihan model.

scot %<>%
  group_by(src_sub_area, sample) %>%
  nest(.key = "data") %>%
  spread(sample, data)

scot

## # A tibble: 3 x 3
##   src_sub_area test               train               
##   <chr>        <list>             <list>              
## 1 sxk97        <tibble [168 × 2]> <tibble [1,344 × 2]>
## 2 sxk9e        <tibble [168 × 2]> <tibble [1,344 × 2]>
## 3 sxk9s        <tibble [168 × 2]> <tibble [1,344 × 2]>

Membuat object time series dan multiple seasonality time series. Saya akan mencoba kedua jenis objek time series pada pemodelan dan melihat objek mana yang memberikan error terkecil saat model melakukan prediksi.

# Membuat list dari fungsi objek
data_funs <- list(
  ts = function(x) ts(x$demand, frequency = 24),
  msts = function(x) msts(x$demand, seasonal.periods = c(24, 24 * 7))
)

# Mengubah bentuk fungsi kedalam data frame yang siap disatukan dengan dataset
data_funs %<>%
  rep(length(unique(scot$src_sub_area))) %>%
  enframe("data_fun_name", "data_fun") %>%
  mutate(src_sub_area =
    sort(rep(unique(scot$src_sub_area), length(unique(.$data_fun_name))))
  )

# Penggabungan dengan dataset
scot %<>% left_join(data_funs)

## Joining, by = "src_sub_area"

head(scot)

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

2 Model

Selanjutnya, membuat daftar algoritma model yang ingin diaplikasikan pada dataset. Saya menggunakan model:
1. Exponential Smoothing State Space Model (ets)
2. Seasonal and Trend decomposition using Loess (stlm)
3. Trigonometric seasonality, Box-Cox transformation, ARMA errors, Trend and Seasonal components (tbats)
4. Autoregressive integrated moving average (arima)
5. Holt Winter

Saya akan mencoba semua model tersebut, lalu melihat model mana yang memberikan error terkecil saat model melakukan prediksi. Khusus untuk model ets dan arima, tidak diaplikasikan pada objek multiple seasonal time series karena model tersebut tidak kompatibel dengan objek tersebut.

# Membuat list dari model 
models <- list(
  ets = function(x) ets(x),
  stlm = function(x) stlm(x),
  tbats = function(x) tbats(x),
  auto.arima = function(x) auto.arima(x),
  holt.winter = function(x) HoltWinters(x,seasonal = "additive")
)

# Mengubah bentuk fungsi kedalam data frame yang siap disatukan dengan dataset
models %<>%
  rep(length(unique(scot$src_sub_area))) %>%
  enframe("model_name", "model") %>%
  mutate(src_sub_area =
    sort(rep(unique(scot$src_sub_area), length(unique(.$model_name))))
  )
models

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

# Penggabungan dengan dataset dan seleksi model ets dan arima pada objek multiple seasonal ts
scot %<>% 
  left_join(models) %>% 
  filter(!(model_name == "ets" & data_fun_name == "msts"),
         !(model_name == "auto.arima" & data_fun_name == "msts"))

## Joining, by = "src_sub_area"

head(scot)

## # A tibble: 6 x 7
##   src_sub_area test      train      data_fun_name data_fun model_name model
##   <chr>        <list>    <list>     <chr>         <list>   <chr>      <lis>
## 1 sxk97        <tibble … <tibble [… ts            <fn>     ets        <fn> 
## 2 sxk97        <tibble … <tibble [… ts            <fn>     stlm       <fn> 
## 3 sxk97        <tibble … <tibble [… ts            <fn>     tbats      <fn> 
## 4 sxk97        <tibble … <tibble [… ts            <fn>     auto.arima <fn> 
## 5 sxk97        <tibble … <tibble [… ts            <fn>     holt.wint… <fn> 
## 6 sxk97        <tibble … <tibble [… msts          <fn>     stlm       <fn>

Mengaplikasikan fungsi pembuatan objek time series dan algoritma pemodelan pada dataset.

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

## Warning in HoltWinters(x, seasonal = "additive"): optimization
## difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

scot

## # A tibble: 24 x 8
##    src_sub_area test   train data_fun_name data_fun model_name model fitted
##    <chr>        <list> <lis> <chr>         <list>   <chr>      <lis> <list>
##  1 sxk97        <tibb… <tib… ts            <fn>     ets        <fn>  <ets> 
##  2 sxk97        <tibb… <tib… ts            <fn>     stlm       <fn>  <stlm>
##  3 sxk97        <tibb… <tib… ts            <fn>     tbats      <fn>  <tbat…
##  4 sxk97        <tibb… <tib… ts            <fn>     auto.arima <fn>  <ARIM…
##  5 sxk97        <tibb… <tib… ts            <fn>     holt.wint… <fn>  <HltW…
##  6 sxk97        <tibb… <tib… msts          <fn>     stlm       <fn>  <stlm>
##  7 sxk97        <tibb… <tib… msts          <fn>     tbats      <fn>  <tbat…
##  8 sxk97        <tibb… <tib… msts          <fn>     holt.wint… <fn>  <HltW…
##  9 sxk9e        <tibb… <tib… ts            <fn>     ets        <fn>  <ets> 
## 10 sxk9e        <tibb… <tib… ts            <fn>     stlm       <fn>  <stlm>
## # … with 14 more rows

3 Evaluasi

Menghitung error pada prediksi yang dilakukan setiap model terhadap data Test.

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

scot %>%
  select(src_sub_area, ends_with("_name"), error)

## # A tibble: 24 x 4
##    src_sub_area data_fun_name model_name  error
##    <chr>        <chr>         <chr>       <dbl>
##  1 sxk97        msts          tbats       10.1 
##  2 sxk97        msts          stlm        11.2 
##  3 sxk97        msts          holt.winter 11.3 
##  4 sxk97        ts            ets         11.8 
##  5 sxk97        ts            tbats       11.9 
##  6 sxk97        ts            holt.winter 12.1 
##  7 sxk97        ts            stlm        12.4 
##  8 sxk97        ts            auto.arima  15.2 
##  9 sxk9e        msts          tbats        7.57
## 10 sxk9e        msts          holt.winter  7.98
## # … with 14 more rows

Selanjutnya, membuat visualisasi yang memperlihatkan bagaimana perbedaan hasil prediksi dari setiap model (warna hijau) jika dibandingan dengan data permintaan riil (warna hitam).

# Melakukan forecast pada data test
scot_test <- scot %>%
  mutate(
    forecast =
      map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
      map2(test, ~ tibble(
        datetime = .y$datetime,
        demand = as.vector(.x$mean)
      )),
    key = paste(data_fun_name, model_name, sep = "-")
  )

# Membentuk data untuk visualisasi 
scot_test %<>%
  select(src_sub_area, key, actual = test, forecast) %>%
  spread(key, forecast) %>%
  gather(key, value, -src_sub_area) %>% 
  unnest(value) %>% 
  mutate(demand = rec_revert(demand,rec,src_sub_area))

# Visualisasi
color <- "forestgreen"

ggplotly(ggplot(scot_test,aes(x = datetime, y = demand)) +
           geom_line(data = scot_test %>% filter(key == "actual"),aes(y = demand),alpha = 0.2,size = 0.8)+
           geom_line(data = scot_test %>% filter(key != "actual"),aes(frame = key,col = key)) +
           labs(x = "", y = "Permintaan (order)",title = "PERBANDINGAN HASIL PREDIKSI MODEL", frame = "Models") +
           facet_wrap(~ src_sub_area, scale = "free_y", ncol = 1) +
           scale_color_manual(values = c(color,color,color,color,color,color,color,color,color)) +
           theme_test(base_family = "Bebas Neue") +
           theme(legend.position = "none"))

## Warning: Ignoring unknown aesthetics: frame

## Warning in p$x$data[firstFrame] <- p$x$frames[[1]]$data: number of items to
## replace is not a multiple of replacement length

Memilih model yang memiliki error prediksi terkecil untuk setiap sub-area.

# Seleksi error terkecil
scot %<>%
  select(-fitted) %>% # remove unused
  group_by(src_sub_area) %>%
  filter(error == min(error)) %>%
  ungroup()

# Menggabungkan data test dan train
scot %<>%
  mutate(fulldata = map2(train, test, ~ bind_rows(.x, .y))) %>%
  select(src_sub_area, fulldata, everything(), -train, -test)

scot

## # A tibble: 3 x 7
##   src_sub_area fulldata       data_fun_name data_fun model_name model error
##   <chr>        <list>         <chr>         <list>   <chr>      <lis> <dbl>
## 1 sxk97        <tibble [1,51… msts          <fn>     tbats      <fn>  10.1 
## 2 sxk9e        <tibble [1,51… msts          <fn>     tbats      <fn>   7.57
## 3 sxk9s        <tibble [1,51… msts          <fn>     tbats      <fn>   8.28

Model tbats memiliki error prediksi terkecil untuk setiap sub-area. Selanjutnya menggabungkan data Test dan data Train untuk membuat model final.

Melakukan pemodelan dengan model tbats terhadap dataset yang sudah di gabung, lalu melakukan prediksi permintaan selama 7 hari kedepan.

# Menjalankan model
scot %<>%
  mutate(
    params = map(fulldata, ~ list(x = .x)),
    data = invoke_map(data_fun, params),
    params = map(data, ~ list(x = .x)),
    fitted = invoke_map(model, params)
  ) %>%
  select(-data, -params)

# Melakukan prediksi
scot %<>%
  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)
    ))
  )

scot

## # A tibble: 3 x 9
##   src_sub_area fulldata data_fun_name data_fun model_name model error
##   <chr>        <list>   <chr>         <list>   <chr>      <lis> <dbl>
## 1 sxk97        <tibble… msts          <fn>     tbats      <fn>  10.1 
## 2 sxk9e        <tibble… msts          <fn>     tbats      <fn>   7.57
## 3 sxk9s        <tibble… msts          <fn>     tbats      <fn>   8.28
## # … with 2 more variables: fitted <list>, forecast <list>

Membuka nest data dan membuat visualisasi dari hasil prediksi.

# Membuka nest data
scot %<>%
  select(src_sub_area, actual = fulldata, forecast) %>%
  gather(key, value, -src_sub_area) %>%
  unnest(value) %>% 
  mutate(demand = rec_revert(demand,rec,src_sub_area))

# Visualisasi
ggplotly(ggplot(scot,aes(x = datetime, y = demand, colour = key)) +
           geom_line() +
           labs(y = "Permintaan (Order)", x = NULL, title = "HASIL PREDIKSI MODEL") +
           facet_wrap(~ src_sub_area, scale = "free_y", ncol = 1) +
           scale_color_brewer(palette = "Pastel2") +
           theme_test(base_family = "Bebas Neue") +
           theme(legend.position = "none"))

# predict <- scot %>% 
#   filter(key == "forecast") %>% 
#   select(-key)
# 
# head(predict)
# 
# write_csv(predict,"data_input/scotty-ts-auto/data/submission-0001HB.csv")