Time Series and Forecasting Drug Sales in Pharma Retail

by Delvia

1. Introduction

Pembentukan model time series dan forecasting di pharma retail sangat bermanfaat dalam memperkirakan penjualan obat maupun sediaan farmasi lain pada periode selanjutnya. Dengan meramalkan permintaan obat, pharma retail dapat memastikan bahwa mereka memiliki stok yang cukup untuk memenuhi kebutuhan pelanggan serta meminimalkan kelebihan stok. Ini dapat membantu mengoptimalkan inventory, mengurangi kehabisan stok serta meningkatkan kepuasan pelanggan.

Time series dan forecasting juga dapat membantu pharma retail dalam merencanakan strategi penjualan secara lebih efektif. Dengan memprediksi produk mana yang cenderung laku dan kapan, hal ini akan dapat membantu strategi pemasaran dan promosi sehingga, dapat membantu meningkatkan kinerja penjualan dan memaksimalkan pendapatan. sebagai tambahan, manfaat lain adalah meningkatkan efisiensi operasional dalam mengalokasikan ruang penyimpanan obat misalnya sehingga hal ini akan mengurangi biaya operasional.

2. Import Library

Berikut merupakan packages yang digunakan untuk analisis data serta pembuatan model “Time series dan Forecasting”

# load library
library(dplyr) #data manipulation
library(lubridate) # date manipulation
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
library(tseries) # adf.test
library(fpp) # usconsumtion
library(padr)
library(psych)
library(fpp)
library(ggplot2)

3. Read Data

drug_sales <- read.csv("data_input/salesmonthly.csv", na.strings = c("N-A", " ", "N/A", "0", "0.0", "0.00", "0.000", "0.0000"))
rmarkdown::paged_table(drug_sales)

Description:
- M01AB : Anti-inflammatory and antirheumatic products, non-steroids, Acetic acid derivatives and related substances
- M01AE : Anti-inflammatory and antirheumatic products, non-steroids, Propionic acid derivatives
- N02BA : Other analgesics and antipyretics, Salicylic acid and derivatives
- N02BE : Other analgesics and antipyretics, Pyrazolones and Anilides
- N05B : Psycholeptics drugs, Anxiolytic drugs
- N05C : Psycholeptics drugs, Hypnotics and sedatives drugs
- R03 : Drugs for obstructive airway diseases
- R06 : Antihistamines for systemic use

4. Data Wrangling

Adapun syarat agar data dapat dibentuk menjadi model time series dan forecasting adalah sebagai berikut:
Syarat 1: Data harus terurut berdasarkan periode waktunya
Syarat 2: Tidak boleh ada waktu atau periode yang terlewat
Syarat 3: Data tidak boleh ada yang missing
Syarat 4: Interval waktu harus sama

glimpse(drug_sales)
#> Rows: 70
#> Columns: 9
#> $ datum <chr> "2014-01-31", "2014-02-28", "2014-03-31", "2014-04-30", "2014-05…
#> $ M01AB <dbl> 127.69, 133.32, 137.44, 113.10, 101.79, 112.07, 117.06, 134.79, …
#> $ M01AE <dbl> 99.090, 126.050, 92.950, 89.475, 119.933, 94.710, 95.010, 99.780…
#> $ N02BA <dbl> 152.100, 177.000, 147.655, 130.900, 132.100, 122.900, 129.300, 1…
#> $ N02BE <dbl> 878.030, 1001.900, 779.275, 698.500, 628.780, 548.225, 491.900, …
#> $ N05B  <dbl> 354.0, 347.0, 232.0, 209.0, 270.0, 323.0, 348.0, 420.0, 399.0, 4…
#> $ N05C  <dbl> 50, 31, 20, 18, 23, 23, 21, 29, 14, 30, 19, 25, 24, 9, 13, 5, 10…
#> $ R03   <dbl> 112.00, 122.00, 112.00, 97.00, 107.00, 57.00, 61.00, 37.00, 115.…
#> $ R06   <dbl> 48.20, 36.20, 85.40, 73.70, 123.70, 109.30, 69.10, 70.80, 58.80,…

Pre-processing Data

drug_sales <- drug_sales %>% 
  mutate(datum=as.Date(datum)) %>% 
  rename(date=datum)

drug_sales <- drug_sales[-70,]
drug_sales$monthly_date <- floor_date(drug_sales$date, "month")
drug_sales$month <- strftime(drug_sales$date, "%m")
drug_sales$year <- strftime(drug_sales$date, "%Y")
drug_sales_clean <- drug_sales[,c(10,12,2:9)]

Cek missing value

is.na(drug_sales_clean) %>% colSums()
#> monthly_date         year        M01AB        M01AE        N02BA        N02BE 
#>            0            0            1            1            1            1 
#>         N05B         N05C          R03          R06 
#>            0            1            1            1
summary(drug_sales_clean)
#>   monthly_date            year               M01AB           M01AE       
#>  Min.   :2014-01-01   Length:69          Min.   :101.8   Min.   : 88.27  
#>  1st Qu.:2015-06-01   Class :character   1st Qu.:138.0   1st Qu.:103.54  
#>  Median :2016-11-01   Mode  :character   Median :155.2   Median :114.98  
#>  Mean   :2016-10-30                      Mean   :153.8   Mean   :119.39  
#>  3rd Qu.:2018-04-01                      3rd Qu.:169.4   3rd Qu.:128.65  
#>  Max.   :2019-09-01                      Max.   :211.1   Max.   :222.35  
#>                                          NA's   :1       NA's   :1       
#>      N02BA            N02BE             N05B            N05C      
#>  Min.   : 75.20   Min.   : 479.4   Min.   :  1.0   Min.   : 5.00  
#>  1st Qu.: 98.46   1st Qu.: 651.7   1st Qu.:226.0   1st Qu.:12.00  
#>  Median :118.67   Median : 873.0   Median :250.6   Median :18.00  
#>  Mean   :118.10   Mean   : 914.5   Mean   :264.7   Mean   :18.26  
#>  3rd Qu.:133.94   3rd Qu.:1068.2   3rd Qu.:295.2   3rd Qu.:23.00  
#>  Max.   :191.60   Max.   :1856.8   Max.   :492.0   Max.   :50.00  
#>  NA's   :1        NA's   :1                        NA's   :1      
#>       R03             R06        
#>  Min.   : 37.0   Min.   : 30.60  
#>  1st Qu.:114.2   1st Qu.: 51.15  
#>  Median :162.0   Median : 75.70  
#>  Mean   :172.1   Mean   : 89.05  
#>  3rd Qu.:219.2   3rd Qu.:121.15  
#>  Max.   :386.0   Max.   :213.04  
#>  NA's   :1       NA's   :1
# handle missing value
drug_sales_clean <- drug_sales_clean %>% 
  mutate_if(is.numeric, na.aggregate, FUN = median)

#recheck missing value
is.na(drug_sales_clean) %>% colSums()
#> monthly_date         year        M01AB        M01AE        N02BA        N02BE 
#>            0            0            0            0            0            0 
#>         N05B         N05C          R03          R06 
#>            0            0            0            0

Cek duplikat data

anyDuplicated(drug_sales_clean)
#> [1] 0
drug_sales_clean %>% pad()
#>    monthly_date year  M01AB    M01AE   N02BA     N02BE  N05B N05C    R03    R06
#> 1    2014-01-01 2014 127.69  99.0900 152.100  878.0300 354.0   50 112.00  48.20
#> 2    2014-02-01 2014 133.32 126.0500 177.000 1001.9000 347.0   31 122.00  36.20
#> 3    2014-03-01 2014 137.44  92.9500 147.655  779.2750 232.0   20 112.00  85.40
#> 4    2014-04-01 2014 113.10  89.4750 130.900  698.5000 209.0   18  97.00  73.70
#> 5    2014-05-01 2014 101.79 119.9330 132.100  628.7800 270.0   23 107.00 123.70
#> 6    2014-06-01 2014 112.07  94.7100 122.900  548.2250 323.0   23  57.00 109.30
#> 7    2014-07-01 2014 117.06  95.0100 129.300  491.9000 348.0   21  61.00  69.10
#> 8    2014-08-01 2014 134.79  99.7800 123.800  583.8500 420.0   29  37.00  70.80
#> 9    2014-09-01 2014 108.78 109.0940 122.100  887.8200 399.0   14 115.00  58.80
#> 10   2014-10-01 2014 154.75 185.2410 191.600 1856.8150 472.0   30 182.00  74.50
#> 11   2014-11-01 2014 138.08 100.8600 142.700  723.8000 489.0   19 112.00  45.20
#> 12   2014-12-01 2014 131.90 121.4010 111.124 1015.6600 492.0   25 163.00  33.40
#> 13   2015-01-01 2015 135.91 130.3490 141.000 1044.2400 463.0   24 177.25  42.00
#> 14   2015-02-01 2015 115.71 123.7400 131.830  953.2520 243.0    9 208.00  47.00
#> 15   2015-03-01 2015 156.04 129.3860 133.800 1084.8500 208.0   13 195.00  54.00
#> 16   2015-04-01 2015 154.50 101.1150 122.100  940.1700 192.0    5  97.00 112.00
#> 17   2015-05-01 2015 160.02 119.1170 136.040  765.9000 194.0   10 100.00 159.50
#> 18   2015-06-01 2015 151.87 113.6900 145.460  746.7880 217.0   12 193.00 125.80
#> 19   2015-07-01 2015 175.61 113.8100 125.500  708.8280 203.0    6  60.00 130.30
#> 20   2015-08-01 2015 181.69 144.5190 133.400  790.7880 265.5   15  45.00  83.70
#> 21   2015-09-01 2015 166.22 134.1220 110.400  852.1250 243.5   11  91.00  71.00
#> 22   2015-10-01 2015 195.81 127.2310 146.200 1574.3350 222.0    8 184.00  72.00
#> 23   2015-11-01 2015 152.78 128.2330 145.900 1277.7250 228.0   18 195.00  44.00
#> 24   2015-12-01 2015 159.46 131.2910 137.000 1258.3490 286.0   28 231.00  41.73
#> 25   2016-01-01 2016 171.65 128.4020 172.500 1476.3240 248.0   24 174.00  56.50
#> 26   2016-02-01 2016 173.81 137.5280 134.200 1224.8620 239.0   20 245.00  58.00
#> 27   2016-03-01 2016 156.64 180.5890 148.400 1150.7000 250.0   13 253.00  97.84
#> 28   2016-04-01 2016 166.61 146.5260 147.700  998.3370 318.0   18 216.00 162.40
#> 29   2016-05-01 2016 167.36 120.8610 130.550  997.1500 275.0   18 131.00 137.10
#> 30   2016-06-01 2016 169.67 114.9610 117.750  760.0500 311.0   20 127.00 134.80
#> 31   2016-07-01 2016 203.97 141.0190 137.900  652.3620 240.0    8 109.00 116.83
#> 32   2016-08-01 2016 211.13 114.3750 132.700  753.0500 275.5   12 116.00  85.30
#> 33   2016-09-01 2016 172.96 126.2180 116.700 1118.6990 307.0   18 121.00  69.30
#> 34   2016-10-01 2016 186.76 142.0560 160.150 1617.2750 312.0   11 220.00  60.90
#> 35   2016-11-01 2016 175.18 116.8500 133.850 1062.6860 246.0   27 150.00  51.20
#> 36   2016-12-01 2016 169.32 135.0560 132.400 1624.3350 257.0   18 275.00  34.90
#> 37   2017-01-01 2017 155.18 114.9765 118.675  872.9645   1.0   18 162.00  75.70
#> 38   2017-02-01 2017 139.69 103.5170  97.000  526.3500 144.0    7 117.00  30.60
#> 39   2017-03-01 2017 162.85 111.0550 107.350  612.5000 165.0    9 139.00 100.10
#> 40   2017-04-01 2017 155.61 101.2150 100.500  540.2000 132.0    9 209.00 122.40
#> 41   2017-05-01 2017 143.66 118.1250  98.950  547.9400 148.0   23 128.00 161.81
#> 42   2017-06-01 2017 122.33 103.0060 119.600  496.1000 163.0    8 163.00 151.90
#> 43   2017-07-01 2017 159.67 116.2060  75.200  479.3500 219.0   15 115.00  81.10
#> 44   2017-08-01 2017 170.15 112.4700  84.400  549.3000 239.0   12  75.00  60.10
#> 45   2017-09-01 2017 138.33 118.7110  88.150  863.7500 223.0   23 139.00  66.90
#> 46   2017-10-01 2017 137.64  88.7370 100.400 1184.3500 226.0   15 247.00  51.00
#> 47   2017-11-01 2017 163.85 119.7800 104.450  867.8990 192.0   15 196.00  46.60
#> 48   2017-12-01 2017 160.01 121.6630 115.150 1007.1800 226.0    6 204.00  47.10
#> 49   2018-01-01 2018 132.28 109.4460 101.150 1134.3250 229.0   11 219.00  49.50
#> 50   2018-02-01 2018 128.36 132.8040 114.650 1255.3740 268.0   12 253.00  39.06
#> 51   2018-03-01 2018 146.16 111.7640 122.300  999.1230 381.0   42 269.00  85.50
#> 52   2018-04-01 2018 170.02 107.7230  84.600  836.0370 289.0   21 229.00 197.10
#> 53   2018-05-01 2018 160.52 103.5220  89.400  644.6480 259.0   13 192.00 213.04
#> 54   2018-06-01 2018 141.18 114.2260  86.800  584.3430 248.0   18 101.00 120.80
#> 55   2018-07-01 2018 150.18 132.5490  87.200  679.3500 283.0   19  90.00 122.20
#> 56   2018-08-01 2018 140.00 114.7190  88.250  733.8380 253.0   20 159.00 103.10
#> 57   2018-09-01 2018 153.52 114.9920  86.500 1058.2620 263.0   12 205.00  88.10
#> 58   2018-10-01 2018 144.71 129.4000  76.050 1129.2750 287.0   25 353.00  76.90
#> 59   2018-11-01 2018 172.29 105.4870 102.150  995.1500 252.2   22 311.00  48.40
#> 60   2018-12-01 2018 147.71 113.0240  84.750 1213.9500 254.0   27 384.00  53.10
#> 61   2019-01-01 2019 179.70 222.3510  99.700 1660.6120 295.2   23 386.00  41.30
#> 62   2019-02-01 2019 133.73 142.1550 110.200 1001.2120 249.4   12 226.00  69.50
#> 63   2019-03-01 2019 154.52 113.1180  83.350  941.0500 301.4   19 257.00 169.50
#> 64   2019-04-01 2019 161.39 100.1650  88.100  647.6500 299.4   22 259.00 179.10
#> 65   2019-05-01 2019 168.04  97.2580 104.100  703.5620 265.8   26 322.00 135.40
#> 66   2019-06-01 2019 151.54 101.6270 103.200  610.0000 193.0   25 142.00 156.04
#> 67   2019-07-01 2019 181.00 103.5410  92.800  649.8000 250.6   20 115.00 105.20
#> 68   2019-08-01 2019 181.91  88.2690  84.200  518.1000 237.0   26 145.00  97.30
#> 69   2019-09-01 2019 161.07 111.4370  93.500  984.4800 227.8   16 161.00 109.10

5. Exploratory Data Analysis (EDA)

Mari kita cek persebaran data numerik pada dataset

pairs.panels(drug_sales_clean)

library(tidyr)
data_viz <- drug_sales_clean %>% 
            select(-monthly_date) %>% 
            pivot_longer(cols = c(M01AB,M01AE,N02BA,N02BE,N05B,N05C,R03,R06),
            names_to = "Var2")

data_viz <- data_viz %>% 
      group_by(year, Var2) %>% 
      summarise(sum(value))

data_viz <- as.data.frame(data_viz)
ggplot(data = data_viz, mapping = aes(x = as.integer(year), y = `sum(value)`)) +
  geom_line(mapping = aes(color = Var2)) +
  
  labs(x = "Year",
       y = "Amount of Drug Sales Based on ATC",
       title = "Drug Sales Based on ATC nYear of 2014 to 2019") +
         theme_light()

💡 Insight:
* Dari grafik penjualan obat berdasarkan ATC di atas, diketahui bahwa penjualan terbanyak adalah obat dengan ATC N02BE dan yang terendah adalah atc N05C. Sehingga, kali ini akan dilakukan time series dan forecasting untuk kedua atc tersebut.

6. DRUG SALES ATC N02BE

Pembentukan Data Time Series

N02BE_sales_ts <- ts(data=drug_sales_clean$N02BE,
                     start = c(2014, 1),
                     frequency = 12)

autoplot(N02BE_sales_ts, series="Actual")

## Pengecekan ada/tidaknya trend/seasonal

N02BE_decompose <- decompose(N02BE_sales_ts)
autoplot(N02BE_decompose)

autoplot(N02BE_decompose$x - N02BE_decompose$trend)

autoplot(N02BE_decompose$x - N02BE_decompose$seasonal)

💡 Insight:
* Dari kedua plot trend dan seasonal di atas, dapat disimpulkan bahwa data penjualan obat tidak memiliki trend dan seasonal. Hal ini dapat dilihat bahwa pada plot trend, tidak terdapat trend naik ataupun turun, serta pada plot seasonal grafik yang dihasilkan bersifat random (acak).

Cross validation

# Data Train
N02BE_train <- head(N02BE_sales_ts, -9)

# Data Test
N02BE_test <- tail(N02BE_sales_ts, 9)

Model 1: Exponential Smoothing

Simple Exponential Smoothing

A = additive M = multiplicative N = None Z = auto error = Additive trend = None seasonal = None

N02BE_SES <- ets(y=N02BE_train, model="ANN")
# Lihat hasil modelnya
N02BE_SES$x %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014  878.030 1001.900  779.275  698.500  628.780  548.225
N02BE_SES$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 889.0775 881.8733 960.1437 842.1978 748.4913 670.4266
autoplot(N02BE_SES$fitted, series = "ANN") +
  autolayer(N02BE_train, series = "Actual")

### Double Exponential Smoothing Holt (DES)

Double Exponential Smoothing Holt digunakan ketika data memiliki linear trend dan tidak memiliki seasonal. Model ini mempertimbangkan kemungkinan seri beberapa trend.

N02BE_DES <- ets(N02BE_train, model = "AAN")
# Lihat hasil modelnya
N02BE_DES$x %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014  878.030 1001.900  779.275  698.500  628.780  548.225
N02BE_DES$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 713.7429 825.7928 945.5118 843.6065 755.3648 679.1000
autoplot(N02BE_DES$fitted, series = "AAN") +
  autolayer(N02BE_train, series = "Actual")

### Tuning Model

N02BE_DES_tuned <- ets(y=N02BE_train, model="AAN", alpha=0.9)
# Lihat hasil modelnya
N02BE_DES_tuned$x %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014  878.030 1001.900  779.275  698.500  628.780  548.225
N02BE_DES_tuned$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 691.6334 867.7850 996.8966 809.4235 717.9676 646.0651
autoplot(N02BE_DES_tuned$fitted, series = "AAN") +
  autolayer(N02BE_train, series = "Actual")

### Holt’s Exponential

N02BE_holts <-  HoltWinters(x = N02BE_train,
                           alpha = 0.8,
                           beta = F,
                           gamma = F)
# Lihat hasil modelnya
N02BE_holts$x %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014  878.030 1001.900  779.275  698.500  628.780  548.225
N02BE_holts$fitted %>% head()
#>              xhat    level
#> Feb 2014 878.0300 878.0300
#> Mar 2014 977.1260 977.1260
#> Apr 2014 818.8452 818.8452
#> May 2014 722.5690 722.5690
#> Jun 2014 647.5378 647.5378
#> Jul 2014 568.0876 568.0876
N02BE_train %>% 
  autoplot() +
  autolayer(N02BE_holts$fitted[,1])

accuracy(N02BE_SES)
#>                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
#> Training set 6.565732 294.4787 207.4496 -4.941181 22.22737 0.8025587 0.00877902
accuracy(N02BE_DES)
#>                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
#> Training set 2.535124 295.3906 211.5711 -5.488945 22.83408 0.8185035 0.01664166
accuracy(N02BE_DES_tuned)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 0.1579444 303.7301 213.0842 -4.648966 22.49961 0.8243571
#>                    ACF1
#> Training set -0.2159441
accuracy(N02BE_holts$fitted[,1], N02BE_train)
#>                ME     RMSE      MAE       MPE     MAPE       ACF1 Theil's U
#> Test set 6.279993 299.8968 210.5718 -4.332708 22.18452 -0.1312138  1.010088

💡 Insight:
* Dari hasil RMSE diatas dapat dikatakan bahwa model SES adalah model yang memiliki nilai error paling kecil dibanding dengan model lainnya yaitu 294.4787.

Model 2: ARIMA (Autoregressive Integrated Moving Average)

Cek stasioner pada data

adf.test(N02BE_train)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  N02BE_train
#> Dickey-Fuller = -3.9891, Lag order = 3, p-value = 0.01617
#> alternative hypothesis: stationary

nilai p-value < 0.05, artinya model sudah stasioner.

Model Arima dengan tuning nilai p,d,q

N02BE_train %>% tsdisplay()

Hasil PACF menunjukkan bahwa lag 1 memiliki spike tertinggi dan satu-satunya yang melewati garis putus-putus biru. Sehingga, perlu dilakukan tuning model. Sementara itu, nilai ACF yang melewati garis putus-putus adalah 1,2,6,12,17,18,19,20 dengan spike tertinggi adalah 1 dan 18.

# Please type your code
N02BEarima_101 <- Arima(y = N02BE_train, order = c(1, 0, 1))
N02BEarima_101$residuals %>% pacf()

N02BEarima_102 <- Arima(y = N02BE_train, order = c(1, 0, 2))
N02BEarima_102$residuals %>% pacf()

N02BEarima_106 <- Arima(y = N02BE_train, order = c(1, 0, 6))
N02BEarima_106$residuals %>% pacf()

N02BEarima_1012 <- Arima(y = N02BE_train, order = c(1, 0, 12))
N02BEarima_1012$residuals %>% pacf()

N02BEarima_1017 <- Arima(y = N02BE_train, order = c(1, 0, 17))
N02BEarima_1017$residuals %>% pacf()

N02BEarima_1018 <- Arima(y = N02BE_train, order = c(1, 0, 18))
N02BEarima_1018$residuals %>% pacf()

N02BEarima_1019 <- Arima(y = N02BE_train, order = c(1, 0, 19))
N02BEarima_1019$residuals %>% pacf()

N02BEarima_1020 <- Arima(y = N02BE_train, order = c(1, 0, 20))
N02BEarima_1020$residuals %>% pacf()

accuracy(N02BEarima_101)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 0.4909456 268.8322 199.8872 -7.821485 22.74628 0.7733021
#>                      ACF1
#> Training set -0.001605361
accuracy(N02BEarima_102)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 0.1852868 265.7599 194.4228 -7.622916 22.09939 0.7521619
#>                    ACF1
#> Training set 0.02303351
accuracy(N02BEarima_106)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -2.143463 250.2235 182.0728 -7.094827 20.42269 0.7043836
#>                     ACF1
#> Training set 0.008371126
accuracy(N02BEarima_1012)
#>                     ME     RMSE     MAE       MPE     MAPE      MASE       ACF1
#> Training set -2.507705 215.8824 158.767 -4.995203 17.44314 0.6142207 0.03588007
accuracy(N02BEarima_1017)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -1.829446 199.8284 145.5962 -5.487701 16.44833 0.5632669
#>                    ACF1
#> Training set 0.01770212
accuracy(N02BEarima_1018)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -2.079411 201.3239 147.2707 -5.531854 16.68026 0.5697449
#>                    ACF1
#> Training set 0.01903365
accuracy(N02BEarima_1019)
#>                      ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -0.2622911 199.3973 146.4209 -5.223848 16.55729 0.5664572
#>                    ACF1
#> Training set 0.01672197
accuracy(N02BEarima_1020)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -0.269586 199.5495 146.4173 -5.230184 16.54837 0.5664435
#>                    ACF1
#> Training set 0.01633607

Dari hasil tuning model di atas, maka model ARIMA N02BEarima_1019 dapat digunakan karena tidak terdapat lag yang mengalami autocol, serta memiliki nilai RMSE terkecil yaitu 199.3973.

Model Auto Arima

N02BEauto_arima <- auto.arima(y = N02BE_train, seasonal = F)
accuracy(N02BEauto_arima)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 0.5201541 269.3049 200.9224 -7.859875 22.87685 0.7773068
#>                    ACF1
#> Training set -0.0288778

Jika dibandingkan semua model yang telah dibuat diperoleh bahwa model dengan RMSE yang kecil adalah sebagai berikut:
- Model SES : 294.4787
- Model N02BEarima_1019 : 199.3973
- Model auto_arima : 269.3049
Model dengan RMSE terkecil untuk produk obat dengan ATC N02BE adalah N02BEarima_1019. Namun, untuk dijadikan perbandingan maka ketiga model ini tetap akan digunakan pada forecasting.

FORECASTING

N02BEarima_1019_forecast <- forecast(object = N02BEarima_1019,
                                h = 12)
N02BEarima_1019_forecast %>% 
  autoplot() +
  autolayer(N02BE_test)

N02BE_forecast_autoarima <- forecast(object = N02BEauto_arima,
                                h = 12)
N02BE_forecast_autoarima %>% 
  autoplot() +
  autolayer(N02BE_test)

#### Evaluasi Hasil Forecast

# Please type your code
accuracy(N02BEarima_1019_forecast$mean, N02BE_test)
#>                 ME     RMSE      MAE      MPE     MAPE        ACF1 Theil's U
#> Test set -79.99382 289.8798 208.6642 -18.8256 26.57397 -0.05216906 0.8830537
# Please type your code
accuracy(N02BE_forecast_autoarima$mean, N02BE_test)
#>                 ME    RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -98.86939 306.204 245.4844 -23.53855 32.96355 0.2216767 0.9589186

Uji Asumsi

Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena dengan menggunakan residual data, kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model. Metode forecasting yang baik menghasilkan nilai residual berikut ini:
1. Residual yang tidak berkorelasi : Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast
2. Residual memiliki rata-rata 0 : Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya

  • Norma Distribution Check
hist(N02BEarima_1019_forecast$residuals)

hist(N02BE_forecast_autoarima$residuals)

* No-autocorrelation residual yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

# menggunakan Ljung-Box test
Box.test(x=N02BEarima_1019$residuals, type="Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  N02BEarima_1019$residuals
#> X-squared = 0.017631, df = 1, p-value = 0.8944

p-value>0.05 —> tidak ada korelasi pada residual

# menggunakan Ljung-Box test
Box.test(x=N02BEauto_arima$residuals, type="Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  N02BEauto_arima$residuals
#> X-squared = 0.05258, df = 1, p-value = 0.8186

p-value>0.05 —> tidak ada korelasi pada residual

7. DRUG SALES ATC N05C

Pembentukan Data Time series

N05C_sales_ts <- ts(data=drug_sales_clean$N05C,
                     start = c(2014, 1),
                     frequency = 12)

autoplot(N05C_sales_ts, series="Actual")

## Pengecekan ada/tidaknya trend/seasonal

N05C_decompose <- decompose(N05C_sales_ts)
autoplot(N05C_decompose)

autoplot(N05C_decompose$x - N05C_decompose$trend)

autoplot(N05C_decompose$x - N05C_decompose$seasonal)

💡 Insight:
Dari kedua plot trend dan seasonal di atas, dapat disimpulkan bahwa data penjualan obat tidak memiliki trend dan seasonal. Hal ini dapat dilihat bahwa pada plot trend, tidak terdapat trend naik ataupun turun, serta pada plot seasonal grafik yang dihasilkan bersifat random (acak).

Cross validation

# Data Train
N05C_train <- head(N05C_sales_ts, -9)

# Data Test
N05C_test <- tail(N05C_sales_ts, 9)

Model 1: Exponential Smoothing

Simple Exponential Smoothing

A = additive M = multiplicative N = None Z = auto error = Additive trend = None seasonal = None

N05C_SES <- ets(y=N05C_train, model="ANN")
# Lihat hasil modelnya
N05C_SES$x %>% head()
#>      Jan Feb Mar Apr May Jun
#> 2014  50  31  20  18  23  23
N05C_SES$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 30.46741 35.61175 34.39714 30.60534 27.28544 26.15678
autoplot(N05C_SES$fitted, series = "ANN") +
  autolayer(N05C_train, series = "Actual")

### Double Exponential Smoothing Holt (DES)

Double Exponential Smoothing Holt digunakan ketika data memiliki linear trend dan tidak memiliki seasonal. Model ini mempertimbangkan kemungkinan seri beberapa trend.

N05C_DES <- ets(N05C_train, model = "AAN")
# Lihat hasil modelnya
N05C_DES$x %>% head()
#>      Jan Feb Mar Apr May Jun
#> 2014  50  31  20  18  23  23
N05C_DES$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 32.27091 29.28780 26.84020 24.81972 23.16418 21.83696
autoplot(N05C_DES$fitted, series = "AAN") +
  autolayer(N05C_train, series = "Actual")

### Tuning Model

N05C_DES_tuned <- ets(y=N05C_train, model="AAN", alpha=0.9)
# Lihat hasil modelnya
N05C_DES_tuned$x %>% head()
#>      Jan Feb Mar Apr May Jun
#> 2014  50  31  20  18  23  23
N05C_DES_tuned$fitted %>% head()
#>           Jan      Feb      Mar      Apr      May      Jun
#> 2014 48.12082 49.44244 32.47276 20.87454 17.91443 22.11893
autoplot(N05C_DES_tuned$fitted, series = "AAN") +
  autolayer(N05C_train, series = "Actual")

### Holt’s Exponential

N05C_holts <-  HoltWinters(x = N05C_train,
                           alpha = 0.8,
                           beta = F,
                           gamma = F)
# Lihat hasil modelnya
N05C_holts$x %>% head()
#>      Jan Feb Mar Apr May Jun
#> 2014  50  31  20  18  23  23
N05C_holts$fitted %>% head()
#>              xhat    level
#> Feb 2014 50.00000 50.00000
#> Mar 2014 34.80000 34.80000
#> Apr 2014 22.96000 22.96000
#> May 2014 18.99200 18.99200
#> Jun 2014 22.19840 22.19840
#> Jul 2014 22.83968 22.83968
N05C_train %>% 
  autoplot() +
  autolayer(N05C_holts$fitted[,1])

accuracy(N05C_SES)
#>                     ME     RMSE      MAE      MPE    MAPE      MASE       ACF1
#> Training set -0.549915 8.089737 5.950332 -22.6665 42.3755 0.6332947 0.07272393
accuracy(N05C_DES)
#>                     ME    RMSE      MAE       MPE     MAPE      MASE      ACF1
#> Training set 0.3686911 7.37536 5.851266 -19.42458 42.62516 0.6227511 0.1584898
accuracy(N05C_DES_tuned)
#>                       ME     RMSE    MAE       MPE     MAPE      MASE
#> Training set 0.005636369 8.958351 6.8903 -14.41829 46.27911 0.7333357
#>                    ACF1
#> Training set -0.3251323
accuracy(N05C_holts$fitted[,1], N05C_train)
#>                  ME     RMSE      MAE       MPE     MAPE       ACF1 Theil's U
#> Test set -0.5078733 8.805838 6.664687 -18.17166 45.41749 -0.2342406  0.924305

💡 Insight:
* Dari hasil RMSE diatas dapat dikatakan bahwa model DES adalah model yang memiliki nilai error paling kecil dibanding dengan model lainnya yaitu 7.37536.

Model 2: ARIMA (Autoregressive Integrated Moving Average)

Cek stasioner data

adf.test(N05C_train)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  N05C_train
#> Dickey-Fuller = -2.9343, Lag order = 3, p-value = 0.1971
#> alternative hypothesis: stationary

nilai p-value > 0.05, artinya model belum stasioner. Mari kita coba implementasikan fungsi diff() pada data N05C_train.

N05C_diff1 <- diff(N05C_train,differences = 1)

N05C_diff1 %>% adf.test()
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -6.368, Lag order = 3, p-value = 0.01
#> alternative hypothesis: stationary

Jumlah differencing yang kita lakukan akan menjadi nilai untuk d

Mendapatkan nilai p

N05C_diff1 %>% tsdisplay()

- Lag yang melewati garis putus-putus biru pada visual PACF adalah 1 - Lag yang melewati garis putus-putus biru pada visual ACF adalah 1

Model Arima dengan tuning nilai p,d,q

  • Additional Tips Dari Pembuat ARIMA, selalu masukan p & q = 0 List Kombinasi ARIMA:
  • p: 0, 1
  • d: 1
  • q: 0, 1
# Please type your code
N05Carima_010 <- Arima(y = N05C_train, order = c(0, 1, 0))
N05Carima_011 <- Arima(y = N05C_train, order = c(0, 1, 1))
N05Carima_110 <- Arima(y = N05C_train, order = c(1, 1, 0))
N05Carima_111 <- Arima(y = N05C_train, order = c(1, 1, 1))
N05Carima_000 <- Arima(y = N05C_train, order = c(0, 0, 0))
N05Carima_001 <- Arima(y = N05C_train, order = c(0,0, 1))
N05Carima_100 <- Arima(y = N05C_train, order = c(1,0, 0))
N05Carima_101 <- Arima(y = N05C_train, order = c(1,0, 1))
# Please type your code
N05Carima_010 <- Arima(y = N05C_train, order = c(0, 1, 0))
N05Carima_010$residuals %>% pacf()

N05Carima_011 <- Arima(y = N05C_train, order = c(0, 1, 1))
N05Carima_011$residuals %>% pacf()

N05Carima_110 <- Arima(y = N05C_train, order = c(1, 1, 0))
N05Carima_110$residuals %>% pacf()

N05Carima_111 <- Arima(y = N05C_train, order = c(1, 1, 1))
N05Carima_111$residuals %>% pacf()

N05Carima_000 <- Arima(y = N05C_train, order = c(0, 0, 0))
N05Carima_000$residuals %>% pacf()

N05Carima_001 <- Arima(y = N05C_train, order = c(0, 0, 1))
N05Carima_001$residuals %>% pacf()

N05Carima_100 <- Arima(y = N05C_train, order = c(1, 0, 0))
N05Carima_100$residuals %>% pacf()

N05Carima_101 <- Arima(y = N05C_train, order = c(1, 0, 1))
N05Carima_101$residuals %>% pacf()

accuracy(N05Carima_010)
#>                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
#> Training set -0.3825 9.292472 7.150833 -17.08497 48.76925 0.7610643 -0.4001772
accuracy(N05Carima_011)
#>                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
#> Training set -1.088989 8.092255 5.846669 -23.75094 42.17936 0.6222619 0.1433572
accuracy(N05Carima_110)
#>                      ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -0.5496674 8.463652 6.226353 -17.89377 42.44678 0.6626717
#>                      ACF1
#> Training set -0.005616885
accuracy(N05Carima_111)
#>                     ME    RMSE      MAE       MPE     MAPE      MASE      ACF1
#> Training set -1.227376 8.03456 5.964767 -24.89806 43.41019 0.6348311 0.0438414
accuracy(N05Carima_000)
#>                                    ME     RMSE   MAE     MPE     MAPE      MASE
#> Training set -0.000000000000005032939 8.421847 6.365 -25.667 47.52854 0.6774279
#>                  ACF1
#> Training set 0.260002
accuracy(N05Carima_001)
#>                      ME     RMSE      MAE       MPE     MAPE     MASE
#> Training set -0.1014243 8.147305 6.007301 -24.72756 44.84804 0.639358
#>                     ACF1
#> Training set 0.002915856
accuracy(N05Carima_100)
#>                     ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
#> Training set -0.217103 8.023137 5.885973 -24.03731 43.54247 0.626445 -0.1169403
accuracy(N05Carima_101)
#>                      ME     RMSE      MAE      MPE     MAPE      MASE
#> Training set -0.3856623 7.966158 5.840999 -24.2865 42.97569 0.6216584
#>                     ACF1
#> Training set -0.07893196

Dari hasil tuning model di atas, maka model ARIMA N05Carima_101 dapat digunakan karena tidak terdapat lag yang mengalami autocol, serta memiliki nilai RMSE terkecil yaitu 7.966158.

Model auto arima

N05Cauto_arima <- auto.arima(y = N05C_train, seasonal = F)
accuracy(N05Cauto_arima)
#>                     ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
#> Training set -0.217103 8.023137 5.885973 -24.03731 43.54247 0.626445 -0.1169403

💡 Insight:
Jika dibandingkan semua model yang telah dibuat diperoleh bahwa model dengan RMSE yang kecil adalah sebagai berikut:
- Model DES : 7.37536
- Model N05Carima_101 : 7.966158
- Model auto_arima : 8.023137
Model dengan RMSE terkecil untuk produk obat dengan ATC N05C adalah model DES. Namun, untuk dijadikan perbandingan maka ketiga model ini tetap akan digunakan pada forecasting.

FORECASTING

N05C_DES_forecast <- forecast(object = N05C_DES,
                                h = 12)
N05C_DES_forecast %>% 
  autoplot() +
  autolayer(N05C_test)

N05Carima_101_forecast <- forecast(object = N05Carima_101,
                                h = 12)
N05Carima_101_forecast %>% 
  autoplot() +
  autolayer(N05C_test)

N05C_forecast_autoarima <- forecast(object = N05Cauto_arima,
                                h = 12)
N05C_forecast_autoarima %>% 
  autoplot() +
  autolayer(N05C_test)

#### Evaluasi Hasil Forecast

# Please type your code
accuracy(N05C_DES_forecast$mean, N05C_test)
#>               ME     RMSE      MAE      MPE     MAPE        ACF1 Theil's U
#> Test set 4.71092 6.502487 5.728228 17.79595 26.10542 -0.07187641 0.9014166
# Please type your code
accuracy(N05Carima_101_forecast$mean, N05C_test)
#>                ME    RMSE      MAE       MPE     MAPE     ACF1 Theil's U
#> Test set 1.326797 5.01313 4.124834 0.3208044 21.63511 0.130055 0.7102691
# Please type your code
accuracy(N05C_forecast_autoarima$mean, N05C_test)
#>                ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
#> Test set 2.269708 5.207683 4.376932 5.305743 21.84243 0.1156496 0.7383353

Uji Asumsi

Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena dengan menggunakan residual data, kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model. Metode forecasting yang baik menghasilkan nilai residual berikut ini:
1. Residual yang tidak berkorelasi : Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast
2. Residual memiliki rata-rata 0 : Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

  • Norma Distribution Check
hist(N05C_DES_forecast$residuals)

hist(N05Carima_101_forecast$residuals)

hist(N05C_forecast_autoarima$residuals)

* No-autocorrelation residual yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

# menggunakan Ljung-Box test
Box.test(x=N05C_DES$residuals, type="Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  N05C_DES$residuals
#> X-squared = 1.5838, df = 1, p-value = 0.2082

p-value>0.05 —> tidak ada korelasi pada residual

# menggunakan Ljung-Box test
Box.test(x=N05Carima_101$residuals, type="Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  N05Carima_101$residuals
#> X-squared = 0.39282, df = 1, p-value = 0.5308

p-value>0.05 —> tidak ada korelasi pada residual

# menggunakan Ljung-Box test
Box.test(x=N05Cauto_arima$residuals, type="Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  N05Cauto_arima$residuals
#> X-squared = 0.86222, df = 1, p-value = 0.3531

p-value>0.05 —> tidak ada korelasi pada residual

8. Kesimpulan

  • Untuk forecasting penjualan obat dengan ATC N02BE yaitu golongan analgesics dan antipyretics: Pyrazolones dan Anilides, model forecasting terbaik dengan RMSE terkecil adalah N02BEarima_1019
  • Sementara itu, forecasting penjualan obat dengan ATC N05C yaitu golongan Psycholeptics, Hypnotics dan sedatives, model forecasting terbaik dengan RMSE terkecil adalah N05Carima_101