LBB: Wikipedia Forecasting

1. Pendahuluan

library(pageviews)
library(dplyr)
library(forecast)
library(plotly)
library(TTR)
library(MLmetrics)
library(tseries)

options(scipen = 100, max.print = 1e+06)

Pada LBB (Learning by Building) ini akan dilakukan forecasting dari data wikipedia tahun 2020 sampai dengan 3 Juli 2022. Data wikipedia ini menginformasikan jumlah viewers yang mengakses wikipedia dengan artikel mengenai Indonesia. Saya melakukan perubahan data input dari Rpubs milik saudara Wahyudi Tri (https://rpubs.com/wahyuditr/ts).

wiki <- article_pageviews(project = "en.wikipedia", 
                  article = "Indonesia", 
                  start =as.Date('2020-01-01'), end = as.Date("2022-07-03"))
str(wiki)
#> 'data.frame':    915 obs. of  8 variables:
#>  $ project    : chr  "wikipedia" "wikipedia" "wikipedia" "wikipedia" ...
#>  $ language   : chr  "en" "en" "en" "en" ...
#>  $ article    : chr  "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...
#>  $ access     : chr  "all-access" "all-access" "all-access" "all-access" ...
#>  $ agent      : chr  "all-agents" "all-agents" "all-agents" "all-agents" ...
#>  $ granularity: chr  "daily" "daily" "daily" "daily" ...
#>  $ date       : POSIXct, format: "2020-01-01" "2020-01-02" ...
#>  $ views      : num  10174 12289 11914 12534 12074 ...
head(wiki,3)

2. Menginspeksi Data

Karena dalam Time Series variabel yang digunakan adalah time dan y, data yang digunakan yaitu date dan views

wiki1 <-  wiki %>% select(date, views)
head(wiki1,3)
dim (wiki1)
#> [1] 915   2

Terlihat bahwa data wiki1 terdiri dari 915 baris dan 2 kolom. Kemudian, cek missing value apakah ada atau tidak

anyNA(wiki1)
#> [1] FALSE
colSums(is.na(wiki1))
#>  date views 
#>     0     0

Terlihat bahwa tidak ada missing value, sehingga data siap diolah ke tahap selanjutnya.

3. Mengolah Data

Data Time Series

Untuk membuat object time series, digunakan fungsi ‘ts()’, yang didalamnya terdapat beberapa parameter:

  1. data = nilai y yang ingin diamati adalah views atau jumlah viewers di data wikipedia
  2. start = waktu mulai dari time series yang digunakan adalah 2020
  3. frequency = banyaknya data dalam 1 pola musiman, pola yang digunakan adalah pola tahunan yaitu 365.25
wiki_ts <- ts(wiki1$views, start = c(2020,1), frequency = 365.25)
autoplot(wiki_ts, series = "Actual")

Decompose Data Time Series

wiki_dec <- decompose(wiki_ts)
autoplot(wiki_dec)

Untuk Melihat Pola/Seasonal:

autoplot(wiki_dec$x - wiki_dec$trend)

Kesimpulan: Seasonal datanya adalah acak atau random (tidak memiliki seasonal).

Untuk Melihat Trend Naik atau Turun: menggunakan adjusted seasonal (buang efek seasonal data)

autoplot(wiki_dec$x - wiki_dec$seasonal)

Kesimpulan: Trend tidak naik ataupun turun.

4. Membangun Model

Bagi data menjadi data train dan data test:

wiki_train <- wiki_ts %>% 
  head(-365.25)

wiki_test <- wiki_ts %>% 
  tail(365.25)

Simple Exponential Smoothing

wiki_ses <- ets(wiki_train, model = "ANN")
summary(wiki_ses)
#> ETS(A,N,N) 
#> 
#> Call:
#>  ets(y = wiki_train, model = "ANN") 
#> 
#>   Smoothing parameters:
#>     alpha = 0.2349 
#> 
#>   Initial states:
#>     l = 11956.7147 
#> 
#>   sigma:  2317.196
#> 
#>      AIC     AICc      BIC 
#> 11974.57 11974.61 11987.49 
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
#> Training set -3.496463 2312.971 1016.011 -1.641239 7.645228 0.6124478 0.1666354
autoplot(wiki_ses$fitted, series = "ANN") +
  autolayer(wiki_train, series = "Actual")

Double Exponential Smoothing (Holt)

wiki_holt <- ets(wiki_train, model = "AAN")
autoplot(wiki_holt$fitted, series = "AAN") +
  autolayer(wiki_train, series = "Actual")

summary(wiki_holt)
#> ETS(A,Ad,N) 
#> 
#> Call:
#>  ets(y = wiki_train, model = "AAN") 
#> 
#>   Smoothing parameters:
#>     alpha = 0.2339 
#>     beta  = 0.0001 
#>     phi   = 0.8 
#> 
#>   Initial states:
#>     l = 11495.0792 
#>     b = 431.5964 
#> 
#>   sigma:  2322.702
#> 
#>      AIC     AICc      BIC 
#> 11980.16 11980.31 12006.01 
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
#> Training set -13.33763 2312.101 1014.708 -1.72055 7.638419 0.6116624 0.1669997

Triple Exponential Smoothing (Holtwinter)

wiki_hw <- HoltWinters(wiki_train, beta = F, gamma = F)
summary(wiki_hw)
#>              Length Class  Mode     
#> fitted       1096   mts    numeric  
#> x             549   ts     numeric  
#> alpha           1   -none- numeric  
#> beta            1   -none- logical  
#> gamma           1   -none- logical  
#> coefficients    1   -none- numeric  
#> seasonal        1   -none- character
#> SSE             1   -none- numeric  
#> call            4   -none- call

AUTO ARIMA

Model ARIMA juga bisa dibuat dengan cara automatis, yaitu dengan fungsi auto.arima()

wiki_auto <- auto.arima(wiki_train, seasonal = F)
summary(wiki_auto)
#> Series: wiki_train 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>          ar1     ar2      ma1
#>       0.3295  0.0692  -0.9716
#> s.e.  0.0455  0.0452   0.0159
#> 
#> sigma^2 = 4961046:  log likelihood = -5001.39
#> AIC=10010.79   AICc=10010.86   BIC=10028.01
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE      MPE     MAPE     MASE
#> Training set -4.813909 2219.212 940.9296 -1.74357 7.058432 0.567189
#>                       ACF1
#> Training set -0.0005259233

BANDINGKAN ERRORNYA MASING - MASING

sqrt(wiki_ses$mse)
#> [1] 2312.971
sqrt(wiki_holt$mse)
#> [1] 2312.101
sqrt(MSE(wiki_hw$fitted[,1], wiki_train))
#> [1] 2318
sqrt(MSE(wiki_auto$fitted, wiki_train)) 
#> [1] 2219.212

KESIMPULAN: Dari hasil RMSE diatas dapat dikatakan bahwa MODEL AUTO ARIMA adalah model yang memiliki nilai error paling kecil dibanding dengan model lainnya. Sehingga saya akan menggunakan ini untuk untuk memprediksi data.

5. Forecasting

Menggunakan Model Auto Arima

forecasting_auto <- forecast(wiki_auto, h = 365.25)
forecasting_auto %>% 
  autoplot() +
  autolayer(wiki_test)

6. Evaluasi

Cek Distribusi Normal

shapiro.test(wiki_auto$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  wiki_auto$residuals
#> W = 0.44407, p-value < 0.00000000000000022
hist(wiki_auto$residuals, breaks = 100)

plot(wiki_auto$residuals)

Kesimpulan: p-value < 0.05, forecast’s residuals tidak berdistribusi normal.

Autocorrelation: Box.test - Ljng-Box

Box.test(wiki_auto$residuals, type = "Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  wiki_auto$residuals
#> X-squared = 0.00015268, df = 1, p-value = 0.9901

Kesimpulan: p-value > 0.05, No Autocorrelation.

7. Kesimpulan

Dari 4 model yang digunakan yaitu Simple Exponential Smoothing (SES), Holt, Holtwinter, dan Auto Arima, Model Auto Arima merupakan model yang terbaik dengan nilai error paling kecil, dan Box.test - Ljng-Box nya dengan p-value > 0.05, No Autocorrelation.