Forecasting Jumlah View dari Halaman Wikipedia Game of Thrones

Sekilas tentang Game of Thrones

Game of Thrones adalah serial televisi populer dari Amerika Serikat yang diadaptasi dari seri novel fantasi berjudul A Song of Ice and Fire karya George R. R. Martin. Ditayangkan pertama kali di HBO pada tanggal 17 April 2011, dan akhirnya akan segera tamat pada season ke-8 yang akan tayang perdana pada tanggal 14 April 2019 mendatang.

Pengenalan

Dalam Analisis Time Series ini, saya mencoba melakukan forecasting terhadap perkiraan jumlah view halaman wikipedia Game of Throne dalam kurun waktu satu tahun ke depan.

Manfaat untuk Bisnis

Analisis Time Series dapat digunakan untuk berbagai macam aplikasi bisnis untuk forecasting sebuah kuantitas di masa depan dan menjelaskan historical pattern-nya. Beberapa contoh use-case nya:

Menjelaskan seasonal pattern dalam penjualan suatu barang atau entitas tertentu
Memprediksi jumlah customer yang akan bergabung ataupun yang akan berhenti (churning)
Mengestimasi efek dari peluncuran sebuah produk baru terhadap jumlah penjualan unit
Mendeteksi kejadian yang tidak biasa dan memperkirakan besarnya efek yang ditimbulkan

Membaca Data

Data wiki_got_views.csv merupakan data jumlah view harian dari dari halaman wikipedia Game of Thrones dimulai dari tanggal 01 Juli 2015 hingga 16 Maret 2019. Proses pengunduhan data dilakukan pada file wiki_got.R, dimana saya menggunakan library R pageviews dan fungsi article_pageviews() untuk mendapatkan semua data viewnya.

got <- read.csv("data_input/wiki_got_views.csv")
str(got)

## 'data.frame':    1355 obs. of  8 variables:
##  $ project    : Factor w/ 1 level "wikipedia": 1 1 1 1 1 1 1 1 1 1 ...
##  $ language   : Factor w/ 1 level "en": 1 1 1 1 1 1 1 1 1 1 ...
##  $ article    : Factor w/ 1 level "Game_of_Thrones": 1 1 1 1 1 1 1 1 1 1 ...
##  $ access     : Factor w/ 1 level "all-access": 1 1 1 1 1 1 1 1 1 1 ...
##  $ agent      : Factor w/ 1 level "all-agents": 1 1 1 1 1 1 1 1 1 1 ...
##  $ granularity: Factor w/ 1 level "daily": 1 1 1 1 1 1 1 1 1 1 ...
##  $ date       : Factor w/ 1355 levels "2015-07-01","2015-07-02",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ views      : int  42619 40759 38163 38153 42783 44893 40231 38047 36339 34203 ...

Menggunakan library dplyr saya melakukan perubahan nama dan tipe variabel dengan mutate() dan kemudian dilakukan seleksi variabel yang dibutuhkan untuk proses forecasting selanjutnya.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

got <- got %>% 
  mutate(ds = as.Date(date), y = views) %>% 
  select(ds, y)

Saya menggunakan library prophet untuk proses forecasting dalam analisis ini. Dikembangkan secara open-source oleh tim Data Science milik Facebook, prophet sangat kuat dalam menangani permasalahan data yang hilang, pergeseran tren, dan juga nilai outlier

Untuk dapat menggunakan library ini, dibutuhkan sebuah dataframe dengan dua buah kolom bernama ds untuk data bertipe tanggal dan y sebagai nilai dari time series-nya

library(prophet)

## Loading required package: Rcpp

## Loading required package: rlang

got_prophet_model <- prophet(got)

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Selanjutnya digunakan make_future_dataframe() untuk membuat sebuah objek dataframe baru yang mempunyai tambahan baris data baru yang merupakan data 365 hari (1355 + 365 = 1720) ke depan dari data sebelumnya. Dataframe baru ini mempunyai 1720 baris dengan hanya satu variabel bernama ds.

got_future_df <- make_future_dataframe(got_prophet_model, periods = 365)
tail(got_future_df)

##              ds
## 1715 2020-03-10
## 1716 2020-03-11
## 1717 2020-03-12
## 1718 2020-03-13
## 1719 2020-03-14
## 1720 2020-03-15

Menggunakan fungsi predict() bawaan dari R, dilakukan sebuah forecast menggunakan model prophet dan dataframe baru yang sudah kita buat sebelumnya.

got_pred <- predict(got_prophet_model, got_future_df)
str(got_pred)

## 'data.frame':    1720 obs. of  19 variables:
##  $ ds                        : POSIXct, format: "2015-07-01" "2015-07-02" ...
##  $ trend                     : num  10134 10303 10471 10640 10809 ...
##  $ additive_terms            : num  20226 19223 18878 20976 29916 ...
##  $ additive_terms_lower      : num  20226 19223 18878 20976 29916 ...
##  $ additive_terms_upper      : num  20226 19223 18878 20976 29916 ...
##  $ weekly                    : num  -2711 -4794 -6163 -5007 3094 ...
##  $ weekly_lower              : num  -2711 -4794 -6163 -5007 3094 ...
##  $ weekly_upper              : num  -2711 -4794 -6163 -5007 3094 ...
##  $ yearly                    : num  22937 24017 25041 25983 26821 ...
##  $ yearly_lower              : num  22937 24017 25041 25983 26821 ...
##  $ yearly_upper              : num  22937 24017 25041 25983 26821 ...
##  $ multiplicative_terms      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ multiplicative_terms_lower: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ multiplicative_terms_upper: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ yhat_lower                : num  -14221 -17123 -15794 -13063 -4315 ...
##  $ yhat_upper                : num  75926 75253 75351 79015 83487 ...
##  $ trend_lower               : num  10134 10303 10471 10640 10809 ...
##  $ trend_upper               : num  10134 10303 10471 10640 10809 ...
##  $ yhat                      : num  30360 29525 29349 31616 40724 ...

actual_pred <- cbind(got, got_pred[1:nrow(got), c("yhat_lower", "yhat", "yhat_upper")])
head(actual_pred)

##           ds     y yhat_lower     yhat yhat_upper
## 1 2015-07-01 42619 -14220.623 30359.78   75926.33
## 2 2015-07-02 40759 -17122.713 29525.43   75252.54
## 3 2015-07-03 38163 -15793.764 29349.29   75350.77
## 4 2015-07-04 38153 -13063.304 31616.27   79015.14
## 5 2015-07-05 42783  -4314.695 40724.43   83486.78
## 6 2015-07-06 44893   6718.285 52291.73   98065.36

Dengan kemampuan plotting dari prophet, kita dapat mem-plot forecast berdasarkan observasi data sebelumnya

plot(got_prophet_model, got_pred)

Dari plot di atas, dapat kita simpulkan bahwa view dari halaman ini akan meningkat secara signifikan dimulai dari Q2 hingga Q3 setiap tahunnya, dengan alasan serial GoT tayang perdana setiap musimnya pada waktu tersebut (bulan April hingga September)

Ada pengecualian untuk tahun 2018 dikarenakan serial ini absen tayang pada tahun tersebut.

Terdapat juga satu buah nilai outlier di sekitar pertengahan tahun 2017, dimana saat itu sedang ditayangkan season 7 sebagai season yang paling baru dan ditunggu fans-nya.

Dengan prophet, dapat dilakukan proses decompose menggunakan fungsi prophet_plot_components()

prophet_plot_components(got_prophet_model, got_pred)

Melihat grafik trend, tahun 2018 terjadi penurunan secara signifikan jumlah view dikarenakan hiatus-nya serial ini. Namun menjelang akhir tahun, jumlah view akan semaking meningkat seiring hype dari seri semakin membumbung tinggi.

Hari Minggu dan Senin merupakan dua hari dengan tingkat kenaikan view yang tinggi, yang sangat dipengaruhi oleh waktu penayangan perdana serial ini di hari Minggu malam waktu Amerika Serikat.

Awal bulan April hingga akhir bulan Agustus merupakan saat dimana jumlah view cenderung meningkat, karena merupakan waktu penayangan perdana serial tersebut setiap musimnya.

Forecasting dengan ARIMA

Kali ini saya akan melakukan forecasting dengan ARIMA. ARIMA adalah singkatan dari Auto Regressive Integrated Moving Average dan mempunyai 3 parameter: p, d, q

Penjelasan p, d, q:

p is the order (number of time lags) of the AR model
d is the degree of differencing (number of times the data have had past values subtracted)
q is the order of the moving-average model

ARIMA merupakah model yang populer untuk forecasting dan mempunyai flexible class yang memanfaatkan informasi historical untuk membuat prediksi. Model ini menjadi dasar dari teknik forecasting yang selanjutnya dapat digunakan sebagai fondasi untuk model yang lebih kompleks.

Dalam analisis ini, saya akan mengolah time series dari jumlah view halaman wikipedia Game of Thrones, melakukan fitting ARIMA model, dan membuat forecast sederhana.

Library yang Dibutuhkan

library(ggplot2)
library(forecast)
library(tseries)
library(plotly)

str(got)

## 'data.frame':    1355 obs. of  2 variables:
##  $ ds: Date, format: "2015-07-01" "2015-07-02" ...
##  $ y : int  42619 40759 38163 38153 42783 44893 40231 38047 36339 34203 ...

Manipulasi Data

Di sini saya membatasi cakupan hanya untuk data dari satu tahun terakhir yaitu dari 16 Maret 2018 hingga 16 Maret 2019

got2 <- got %>% 
  mutate(Date = ds, views = y) %>% 
  filter(Date >= "2018-03-16" & Date <= "2019-03-17") %>% 
  select(Date, views)

str(got2)

## 'data.frame':    366 obs. of  2 variables:
##  $ Date : Date, format: "2018-03-16" "2018-03-17" ...
##  $ views: int  19053 17559 18803 17823 16991 16812 17678 16461 16553 17871 ...

ggplotly(ggplot(got2, aes(Date, views)) + geom_line() + ylab("Daily Views") + xlab("Date"))

Ada beberapa outlier dalam plot di atas yang dapat membuat model menjadi bias. R menyediakan sebuah metode untuk menghilangkan outlier dari time series yaitu tsclean() dari package forecast. Fungsi ini mengidentifikasi dan mengganti outlier dengan series smoothing and decomposition.

count_ts <- ts(got2[, c("views")])
got2$clean_cnt <- tsclean(count_ts)

ggplot() + 
  geom_line(data = got2, aes(x = Date, y = clean_cnt)) + ylab("Cleaned Page Views") +
  ylab("Daily Views") + xlab("Date")

Di sini saya mengambil weekly dan monthly moving average untuk proses smoothing sehingga data lebih stabil dan predictable

got2$cnt_ma <- ma(got2$clean_cnt, order = 7)
got2$cnt_ma30 <- ma(got2$clean_cnt, order = 30)

ggplot() +
  geom_line(data = got2, aes(x = Date, y = clean_cnt, colour = "Normal Views")) +
  geom_line(data = got2, aes(x = Date, y = cnt_ma, colour = "Moving Average (Weekly)")) +
  geom_line(data = got2, aes(x = Date, y = cnt_ma30, colour = "Moving Average (Monthly)")) +
  ylab("Daily Views") + xlab("Date")

Pada analisis ini, saya akan membuat model smoothing time series cukup dari weekly moving average (warna hijau)

Proses Decompose Data

Komponenn dari analisis time series adalah seasonality, trend, dan cycle.

Seasonal berarti fluktuasi data yang berhubungan dengan siklus kalender. Sebagai contoh, orang-orang lebih suka membuka halaman wikipedia pada saat musim baru serial tv ini akan segera tayang, dan akan jauh lebih sedikit setelah penayangan berakhir.

Trend berarti keseluruhan pattern dari time series. Contohnya, apakah jumlah view meningkat atau menurun dari waktu ke waktu?

Cycle berarti pattern penurunan atau peningkatan yang tidak ada hubungannya dengan seasonal. Biasanya, trend dan cycle dikelompokkan bersama. Trend-cycle bisa diestimasi menggunakan moving averages

Terakhir, bagian dari time series yang tidak dapat dikategorikan sebagai seasonal, cycle, atau trend akan dianggap sebagai residual atau error

Decomposition merupakan proses meng-esktrak komponen-komponen tersebut.

Langkah pertama, kita hitung seasonal menggunakan stl(). STL adalah fungsi untuk decomposing dan juga forecasting.

count_ma <- ts(na.omit(got2$cnt_ma), frequency = 30)
decomp <- stl(count_ma, s.window = "periodic")
deseasonal_cnt <- seasadj(decomp)
plot(decomp)

Fitting sebuah ARIMA model

Fungsi auto.arima() melakukan otomatisasi terhadap nilai optimal (p,d,q)

auto.arima(deseasonal_cnt, seasonal=FALSE)

## Series: deseasonal_cnt 
## ARIMA(3,1,4) 
## 
## Coefficients:
##          ar1     ar2      ar3     ma1      ma2     ma3     ma4
##       0.4088  0.3820  -0.4727  0.5033  -0.3644  0.4933  0.8958
## s.e.  0.0566  0.0574   0.0556  0.0306   0.0329  0.0262  0.0238
## 
## sigma^2 estimated as 82018:  log likelihood=-2540.11
## AIC=5096.21   AICc=5096.62   BIC=5127.28

Evaluasi dan Iterasi

Setelah kita selesai melakukan fitting pada ARIMA model dan membuat forecast, apakah itu masuk akal? Dapatkan kita mempercayai model ini? Kita dapat memerika ACF dan PACF plot untuk model residual. Jika model urutan parameter dan struktur sudah benar, kita dapat berharap tidak ada autocorellations yang signifikan.

tsdisplay() digunakan untuk memplot model diagnostic.

fit <- auto.arima(deseasonal_cnt, seasonal=FALSE)
tsdisplay(residuals(fit), lag.max=45, main='(3,1,4) Model Residuals')

Berdasarkan ACF/PACF, ada lag di sekitar nilai 7 sehingga kita buat fitting model selanjutnya fit2

fit2 = arima(deseasonal_cnt, order=c(1,1,7))
tsdisplay(residuals(fit2), lag.max=15, main='Seasonal Model Residuals')

Forecasting menggunakan fitted model yang terbaru fit2

fcast <- forecast(fit2, h=30)
plot(fcast)

Kesimpulan

Berdasarkan hasil analisa forecasting dengan library prophet dan juga ARIMA, saya bisa menyimpulkan adanya sebuah pola kenaikan dan penurunan jumlah view dari halaman Wikipedia Game of Thrones ini. Pola tersebut adalah:

Sekitar beberapa hari sebelum penayangan perdana musim terbaru dari Game of Thrones, akan terjadi peningkatan jumlah view secara signifikan
Hari Minggu dan Senin akan menjadi hari dengan trafik tertinggi setiap minggunya selama serial ini masih tayang.
Q1 dan Q4 tiap tahunnya adalah saat terendah jumlah view

References

Prophet Library from Facebook

Game of Thrones Wikipedia

Introduction to Forecasting with ARIMA in R