Pendahuluan

Permintaan darah merupakan bagian penting dalam sistem pelayanan kesehatan karena berkaitan langsung dengan keberlangsungan layanan medis dan tindakan darurat. Pengelolaan permintaan yang kurang tepat dapat menyebabkan ketidakseimbangan antara ketersediaan dan kebutuhan darah. Oleh karena itu, diperlukan pendekatan berbasis data untuk memahami pola permintaan darah dari waktu ke waktu.[1]

Analisis deret waktu (time series) merupakan metode statistik yang digunakan untuk menganalisis data yang diamati secara berurutan berdasarkan waktu. Salah satu model yang umum digunakan dalam analisis deret waktu adalah ARIMA (Autoregressive Integrated Moving Average), yang dikembangkan untuk memodelkan ketergantungan temporal dan pola tren pada data historis. [2], [3] Dalam penelitian ini, model ARIMA digunakan untuk menganalisis dan memprediksi permintaan darah berdasarkan data historis.

#library
library(readxl)
data <- read_excel("C:\\Users\\HP\\OneDrive\\Documents\\Kuliah\\Smt 5\\Komstat 2\\16\\Data UAS.xlsx")
head(data)
## # A tibble: 6 × 2
##   Periode             `Permintaan Darah`
##   <dttm>                           <dbl>
## 1 2012-01-01 00:00:00               1257
## 2 2012-02-01 00:00:00               1198
## 3 2012-03-01 00:00:00               1125
## 4 2012-04-01 00:00:00               1005
## 5 2012-05-01 00:00:00               1057
## 6 2012-06-01 00:00:00                866

Exploratory Data Analysis (EDA)

Exploratory Data Analysis (EDA) dilakukan untuk memperoleh gambaran awal mengenai karakteristik data permintaan darah sebelum dilakukan pemodelan. Tahapan EDA meliputi penyajian statistik deskriptif dan visualisasi deret waktu untuk mengidentifikasi pola umum, fluktuasi, serta kemungkinan adanya tren pada data.

Visualisasi deret waktu digunakan untuk mengamati dinamika perubahan permintaan darah dari waktu ke waktu. Selain itu, histogram digunakan untuk melihat sebaran data serta mendeteksi adanya pencilan. Analisis awal terhadap struktur ketergantungan data juga dilakukan melalui plot Autocorrelation Function (ACF) dan Partial Autocorrelation Function (PACF), yang umum digunakan dalam tahap awal identifikasi model deret waktu.[3]

library(ggplot2)
library(tseries)
## Warning: package 'tseries' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
Permintaan_Darah.ts <- ts(data$`Permintaan Darah`)

ggplot(data.frame(
  Periode = time(Permintaan_Darah.ts),
  Nilai = as.numeric(Permintaan_Darah.ts)
), aes(x = Periode, y = Nilai)) +
  geom_line(color = "red", linewidth = 1) +
  labs(title = "Time Series Permintaan Darah",
       x = "Periode",
       y = "Permintaan Darah") +
  theme_minimal()
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Berdasarkan hasil statistik deskriptif data permintaan darah diperoleh nilai minimum sebesar 866 dan nilai maksimum sebesar 1656. Nilai median sebesar 1281 menunjukkan bahwa setengah dari pengamatan berada di bawah nilai tersebut. Rata-rata (mean) permintaan darah sebesar 1262, yang relatif dekat dengan nilai median, mengindikasikan bahwa data tidak memiliki kemencengan (skewness) yang ekstrem. Plot deret waktu permintaan darah menunjukkan adanya fluktuasi nilai permintaan dari periode ke periode. Meskipun terdapat variasi yang cukup jelas, secara umum terlihat adanya kecenderungan perubahan permintaan sepanjang waktu. Pola fluktuasi ini mengindikasikan bahwa data tidak sepenuhnya bersifat acak dan kemungkinan memiliki ketergantungan antarwaktu.

Selain itu, tidak terlihat pola musiman yang berulang secara konsisten pada setiap periode, namun terdapat beberapa lonjakan dan penurunan tajam pada periode tertentu yang mencerminkan dinamika permintaan darah. Karakteristik ini menunjukkan bahwa pendekatan analisis deret waktu, seperti model ARIMA, relevan untuk digunakan dalam memodelkan data tersebut.

# Histogram
hist(Permintaan_Darah.ts,
     main = "Histogram Permintaan Darah",
     xlab = "Permintaan Darah",
     col = "red")

Histogram permintaan darah menunjukkan bahwa sebaran data cenderung terkonsentrasi pada kisaran nilai tengah, dengan frekuensi tertinggi berada di sekitar nilai rata-rata. Distribusi data tampak relatif simetris dan tidak menunjukkan adanya pencilan (outlier) yang ekstrem.

Pola distribusi seperti ini mengindikasikan bahwa variasi data masih dalam batas yang wajar, sehingga data layak untuk dianalisis lebih lanjut menggunakan metode statistik parametrik, termasuk model deret waktu.

Metode Analisis

Metode analisis yang digunakan dalam penelitian ini adalah model ARIMA (Autoregressive Integrated Moving Average). Sebelum pemodelan dilakukan, data diuji kestasionerannya menggunakan uji Augmented Dickey-Fuller (ADF). Uji ini digunakan untuk menentukan keberadaan akar unit pada data deret waktu.[4]

Apabila data belum stasioner, dilakukan proses differencing hingga diperoleh data yang stasioner. Penentuan orde model ARIMA dilakukan berdasarkan pola pada plot ACF dan PACF dari data hasil differencing sebagaimana disarankan dalam pendekatan Box–Jenkins.[2]

Beberapa kandidat model kemudian dibangun menggunakan data training dan dievaluasi berdasarkan nilai Akaike Information Criterion (AIC) serta hasil uji diagnostik residual. Selanjutnya, performa prediksi model dievaluasi menggunakan Mean Absolute Percentage Error (MAPE) pada data testing sebagai ukuran akurasi peramalan.[5]

Hasil dan Pembahasan

Hasil analisis menunjukkan bahwa beberapa model ARIMA dapat dibangun berdasarkan data training. Evaluasi model dilakukan melalui perbandingan nilai AIC dan pengujian diagnostik residual untuk memastikan bahwa model telah memenuhi asumsi yang diperlukan.

Model terbaik kemudian digunakan untuk menghasilkan prediksi permintaan darah pada periode data testing. Hasil prediksi dibandingkan dengan data aktual untuk menilai akurasi model. Tingkat akurasi peramalan dievaluasi menggunakan Mean Absolute Percentage Error (MAPE) sebagai ukuran kesalahan prediksi.

#Pembagian data
data_training <- Permintaan_Darah.ts[1:44]
data_testing<- Permintaan_Darah.ts[45:48]

Membagi data menjadi dua, training untuk pemodelan dan testing untuk menilai kelayakan model.

#Uji kestasioneran
adf_test_permintaan <- adf.test(data_training)
adf_test_permintaan
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_training
## Dickey-Fuller = -2.6181, Lag order = 3, p-value = 0.3291
## alternative hypothesis: stationary
#Differencing
diff1 <- diff(data_training)
adf_diff1 <- adf.test(diff1)
adf_diff1
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1
## Dickey-Fuller = -3.6847, Lag order = 3, p-value = 0.03792
## alternative hypothesis: stationary
#Plot Data Differensing
plot.ts(diff1, 
        col = "red",
        xlab="Periode", ylab="Permintaan Darah",
        main = "Plot Differensing Permintaan Darah 2012 - 2015")

#plot ACF dan PACF untuk identifikasi model
par(mfrow = c(1, 2))
acf(diff1, main = "Plot ACF Data Bangkitan")
pacf(diff1, main = "Plot PACF Data Bangkitan")

par(mfrow = c(1, 1))

Dilakukan uji kestasioneran pertama kali dan didapat hasil tidak stasioner, lalu dilakukan differencing sebanyak satu kali kemudian diuji kestasioneran lagi dan didapati data sudah stasioner. Selanjutnya, dengan membuat plot acf dan pacf dapat digunakan untuk membuat ordo ARIMA p,d,q dengan p = dari plot pacf d = dari n differencing q = dari plot acf dari output yang dihasilkan, plot acf menunjukkan dropdown pada 1-2, dan plot pacf menunjukkan dropdown pada lag-1

#ACF Max 2 lag dan Pacf max 1 lag
#Model yang didapatkan (1,1,1); (1,1,0); (0,1,1)

#Membangun model ARIMA
#Model 1
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.4.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
model1=arima(data_training, order= c (1,1,2))
summary(model1)
##           Length Class  Mode     
## coef       3     -none- numeric  
## sigma2     1     -none- numeric  
## var.coef   9     -none- numeric  
## mask       3     -none- logical  
## loglik     1     -none- numeric  
## aic        1     -none- numeric  
## arma       7     -none- numeric  
## residuals 44     ts     numeric  
## call       3     -none- call     
## series     1     -none- character
## code       1     -none- numeric  
## n.cond     1     -none- numeric  
## nobs       1     -none- numeric  
## model     10     -none- list
coeftest(model1)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.86792    0.14103 -6.1544 7.537e-10 ***
## ma1  0.33052    0.21763  1.5187    0.1288    
## ma2 -0.23004    0.20577 -1.1179    0.2636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model 2
model2=arima(data_training, order= c (1,1,1))
summary(model2)
##           Length Class  Mode     
## coef       2     -none- numeric  
## sigma2     1     -none- numeric  
## var.coef   4     -none- numeric  
## mask       2     -none- logical  
## loglik     1     -none- numeric  
## aic        1     -none- numeric  
## arma       7     -none- numeric  
## residuals 44     ts     numeric  
## call       3     -none- call     
## series     1     -none- character
## code       1     -none- numeric  
## n.cond     1     -none- numeric  
## nobs       1     -none- numeric  
## model     10     -none- list
coeftest(model2)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)   
## ar1 -0.70075    0.25555 -2.7421 0.006105 **
## ma1  0.21674    0.38228  0.5670 0.570747   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model 3
model3=arima(data_training, order = c (1,1,0))
summary(model3)
##           Length Class  Mode     
## coef       1     -none- numeric  
## sigma2     1     -none- numeric  
## var.coef   1     -none- numeric  
## mask       1     -none- logical  
## loglik     1     -none- numeric  
## aic        1     -none- numeric  
## arma       7     -none- numeric  
## residuals 44     ts     numeric  
## call       3     -none- call     
## series     1     -none- character
## code       1     -none- numeric  
## n.cond     1     -none- numeric  
## nobs       1     -none- numeric  
## model     10     -none- list
coeftest(model3)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.55584    0.12538 -4.4331 9.287e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model 4
model4=arima(data_training, order = c (0,1,1))

#AIC ARIMA dan Signifikansi Parameter
model_accuracy<- data.frame(
  "Model"=c("ARIMA(1,1,2)","ARIMA(1,1,1)", "ARIMA(1,1,0)", "ARIMA(0,1,1)"),
  "AIC"=c(model1$aic, model2$aic, model3$aic, model4$aic))
model_accuracy
##          Model      AIC
## 1 ARIMA(1,1,2) 566.1742
## 2 ARIMA(1,1,1) 565.3265
## 3 ARIMA(1,1,0) 563.6321
## 4 ARIMA(0,1,1) 565.9082
library(forecast)
## Warning: package 'forecast' was built under R version 4.4.3
auto_model <- auto.arima(data_training, stepwise = FALSE, trace = TRUE)
## 
##  ARIMA(0,1,0)                    : 577.5995
##  ARIMA(0,1,0) with drift         : 579.801
##  ARIMA(0,1,1)                    : 566.2082
##  ARIMA(0,1,1) with drift         : 568.4703
##  ARIMA(0,1,2)                    : 568.3746
##  ARIMA(0,1,2) with drift         : 570.8096
##  ARIMA(0,1,3)                    : 565.9397
##  ARIMA(0,1,3) with drift         : Inf
##  ARIMA(0,1,4)                    : 568.1699
##  ARIMA(0,1,4) with drift         : Inf
##  ARIMA(0,1,5)                    : 569.7859
##  ARIMA(0,1,5) with drift         : Inf
##  ARIMA(1,1,0)                    : 563.9321
##  ARIMA(1,1,0) with drift         : 566.2471
##  ARIMA(1,1,1)                    : 565.9419
##  ARIMA(1,1,1) with drift         : 568.378
##  ARIMA(1,1,2)                    : 567.2268
##  ARIMA(1,1,2) with drift         : 569.7956
##  ARIMA(1,1,3)                    : 567.9525
##  ARIMA(1,1,3) with drift         : Inf
##  ARIMA(1,1,4)                    : 570.6496
##  ARIMA(1,1,4) with drift         : Inf
##  ARIMA(2,1,0)                    : 566.1048
##  ARIMA(2,1,0) with drift         : 568.5411
##  ARIMA(2,1,1)                    : Inf
##  ARIMA(2,1,1) with drift         : Inf
##  ARIMA(2,1,2)                    : 566.8645
##  ARIMA(2,1,2) with drift         : Inf
##  ARIMA(2,1,3)                    : 569.5671
##  ARIMA(2,1,3) with drift         : Inf
##  ARIMA(3,1,0)                    : 567.4761
##  ARIMA(3,1,0) with drift         : 570.0446
##  ARIMA(3,1,1)                    : 567.0515
##  ARIMA(3,1,1) with drift         : Inf
##  ARIMA(3,1,2)                    : 569.5617
##  ARIMA(3,1,2) with drift         : Inf
##  ARIMA(4,1,0)                    : 569.9848
##  ARIMA(4,1,0) with drift         : 572.6947
##  ARIMA(4,1,1)                    : 569.756
##  ARIMA(4,1,1) with drift         : Inf
##  ARIMA(5,1,0)                    : 569.9028
##  ARIMA(5,1,0) with drift         : 572.7001
## 
## 
## 
##  Best model: ARIMA(1,1,0)
model_auto <- arima(data_training, order = c(1, 1, 0))
summary(model_auto)
## 
## Call:
## arima(x = data_training, order = c(1, 1, 0))
## 
## Coefficients:
##           ar1
##       -0.5558
## s.e.   0.1254
## 
## sigma^2 estimated as 26062:  log likelihood = -279.82,  aic = 563.63
## 
## Training set error measures:
##                     ME    RMSE      MAE       MPE     MAPE     MASE       ACF1
## Training set -1.003855 159.592 119.9782 -1.399346 10.36392 0.779079 0.02067679

Setelah ordo didapatkan, dilakukan pemodelan dari berbagai kemungkinan ordo yang terbentuk secara manual, lalu keakuratan model dibandingkan menggunakan nilai AIC dan didapati model ARIMA(1,1,0) yang terbaik. Dilakukan juga pemodelan auto dan didapati hasil yang sama, dengan koefisien ar = -0.558.

#Uji non-autokorelasi
Box.test(model3$residuals)
## 
##  Box-Pierce test
## 
## data:  model3$residuals
## X-squared = 0.018811, df = 1, p-value = 0.8909
#Uji Normalitas
ks.test(model3$residuals,"pnorm", mean(model3$residuals), sd(model3$residuals))
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  model3$residuals
## D = 0.098076, p-value = 0.7544
## alternative hypothesis: two-sided

Selanjutnya, dilakukan uji asumsi white noise dan normalitas. Dari output yang dihasilkan, didapati kedua asumsi terpenuhi yang mengindikasikan model ARIMA(1,1,0) layak digunakan.

#Forecasting ARIMA
arima_pred_model3 <- predict(model3, n.ahead = 4)
arima_pred_model3
## $pred
## Time Series:
## Start = 45 
## End = 48 
## Frequency = 1 
## [1] 1180.497 1244.142 1208.765 1228.429
## 
## $se
## Time Series:
## Start = 45 
## End = 48 
## Frequency = 1 
## [1] 161.4369 176.6445 214.4417 234.0820
#Perbandingan hasil
start <- 45
end <- 48
periode <- start:end

Hasil_Perbandingan_Arima <- data.frame(
  Periode = periode,
  Aktual = data_testing,
  Prediksi_ARIMA = arima_pred_model3$pred)
Hasil_Perbandingan_Arima
##   Periode Aktual Prediksi_ARIMA
## 1      45   1216       1180.497
## 2      46   1656       1244.142
## 3      47   1416       1208.765
## 4      48   1620       1228.429
#PLOT ARIMA
library(graphics)
library(ggplot2)
ggplot(arima_pred_model3$pred, aes(x = periode)) +
  geom_line(aes(y = data_testing, color = "Aktual"), size = 1) +
  geom_line(aes(y = arima_pred_model3$pred, color = "Prediksi ARIMA"), 
            linetype = "dashed", size = 1) +
  labs(title = "Perbandingan Prediksi ARIMA vs Aktual",
       y = "Permintaan Darah",
       x = "Periode") +
  scale_color_manual(name = "Legenda", 
                     values = c("Aktual" = "red",
                                "Prediksi ARIMA" = "blue")) +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Setelah diperoleh model terbaik dan diperiksa kelayakannya, forecasting (peramalan) dapat dilakukan. Dengan menggunakan ARIMA(1,1,0) dilakukan forecasting dari periode 45 hingga 48, lalu hasilnya dibandingkan menggunakan tabel dan plot. Secara visual melalui plot, grafik forecasting terlihat cukup jauh dari data aktual, namun untuk memastikannya kita perlu mengecek menggunakan MAPE.

#Perbandingan MAPE
library(Metrics)
## Warning: package 'Metrics' was built under R version 4.4.3
## 
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
## 
##     accuracy
arima_pred_model1 <- predict(model1, n.ahead = 4)
arima_pred_model2 <- predict(model3, n.ahead = 4)
arima_pred_model4 <- predict(model4, n.ahead = 4)

mape_model1 <- mape(data_testing, arima_pred_model1$pred)
mape_model2 <- mape(data_testing, arima_pred_model2$pred)
mape_model3 <- mape(data_testing, arima_pred_model3$pred)
mape_model4 <- mape(data_testing, arima_pred_model4$pred)

Perbandingan_MAPE <- data.frame(
  "Model"=c("ARIMA(1,1,2)", "ARIMA(1,1,1)", "ARIMA(1,1,0)", "ARIMA(0,1,1)"),
  "MAPE"=c(mape_model1, mape_model2, mape_model3, mape_model4))
Perbandingan_MAPE
##          Model      MAPE
## 1 ARIMA(1,1,2) 0.1636994
## 2 ARIMA(1,1,1) 0.1664914
## 3 ARIMA(1,1,0) 0.1664914
## 4 ARIMA(0,1,1) 0.1510546

Melalui MAPE, kita dapat mengetahui keakuratan model. Dari output dapat terlihat bahwa MAPE ARIMA(1,1,0) adalah 0.1664914, dan MAPE antar model bisa dikatakan tidak berbeda secara signifikan. Meskipun model ARIMA(0,1,1) menghasilkan nilai MAPE yang sedikit lebih kecil pada data testing, perbedaannya relatif kecil dan tidak signifikan secara praktis. Oleh karena itu, model ARIMA(1,1,0) tetap dipertahankan sebagai model yang digunakan dalam peramalan

Kesimpulan

Berdasarkan keseluruhan hasil yang didapatkan,

Referensi

[1]
T. A. Rini, “Sistem ketahanan darah (blood security system) di kota yogyakarta menggunakan metode arima box- jenkins dan pendekatan system dynamics (studi kasus: Pmi utd kota jogjakarta),” Universitas Islam Indonesia, 2016.
[2]
G. T. Wilson, “Time series analysis: Forecasting and control, 5th edition, by george e. P. Box, gwilym m. Jenkins, gregory c. Reinsel and greta m. Ljung, 2015. Published by john wiley and sons inc., hoboken, new jersey, pp. 712. ISBN: 978‐1‐118‐67502‐1,” J. Time Ser. Anal., vol. 37, no. 5, pp. 709–711, Sep. 2016.
[3]
R. Hyndman and G. Athanasopoulos, Forecasting: Principles and practice, 3rd ed. Australia: OTexts, 2021.
[4]
D. A. Dickey and W. A. Fuller, “Distribution of the estimators for autoregressive time series with a unit root,” J. Am. Stat. Assoc., vol. 74, no. 366a, pp. 427–431, Jun. 1979.
[5]
C. Beaumont, S. Makridakis, S. C. Wheelwright, and V. E. McGee, “Forecasting: Methods and applications,” J. Oper. Res. Soc., vol. 35, no. 1, p. 79, Jan. 1984.