Pendahuluan

Jumlah kunjungan pasien rawat jalan merupakan salah satu indikator penting dalam pengelolaan Balai Kesehatan Masyarakat (BKM). Informasi mengenai pola kunjungan pasien dibutuhkan untuk mendukung perencanaan sumber daya, penjadwalan tenaga kesehatan, serta peningkatan kualitas (Rustam et al., 2022) pelayanan. Oleh karena itu, peramalan jumlah kunjungan pasien menjadi langkah strategis dalam pengambilan keputusan manajerial berbasis data (Basri et al., n.d.)

Data kunjungan pasien yang dicatat secara bulanan membentuk deret waktu yang umumnya mengandung pola tren dan fluktuasi musiman. Untuk menganalisis karakteristik tersebut, diperlukan metode yang mampu menangkap unsur musiman secara efektif. Model Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan salah satu metode deret waktu yang sesuai karena dapat memodelkan pola tren dan musiman secara bersamaan (Denashurya, 2024).

Penelitian ini bertujuan untuk memodelkan dan meramalkan jumlah kunjungan pasien rawat jalan di Balai Kesehatan Masyarakat menggunakan model SARIMA. Hasil penelitian diharapkan dapat memberikan informasi yang bermanfaat bagi pihak manajemen dalam merencanakan pelayanan kesehatan yang lebih optimal dan efisien.

Deskripsi dan Sumber Data

Data yang digunakan dalam penelitian ini merupakan data sekunder berupa jumlah kunjungan pasien rawat jalan per bulan pada sebuah Balai Kesehatan Masyarakat (BKM). Data mencakup periode Januari 2014 hingga Desember 2018, sehingga diperoleh sebanyak 60 pengamatan. Variabel utama yang dianalisis adalah jumlah kunjungan pasien rawat jalan setiap bulan.

Data bersumber dari jurnal berjudul “Peramalan Jumlah Kunjungan Pasien Rawat Jalan Menggunakan Metode ARIMA, SES, dan Holt-Winters di Balai Kesehatan” yang ditulis oleh Ilham Basri K. dari Universitas Komputer Indonesia, Bandung. Dalam penelitian ini, identitas dan lokasi spesifik Balai Kesehatan Masyarakat disamarkan untuk menjaga kerahasiaan dan privasi institusi yang bersangkutan (Basri et al., n.d.)

Pengolahan dan analisis data dilakukan menggunakan perangkat lunak R dengan bantuan beberapa paket pendukung, antara lain readxl, psych, dplyr, forecast, dan tseries.

#Masukkan data
library(readxl)
data <- read_excel(file.choose())
str(data)
## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Bulan : chr [1:60] "JANUARI" "FEBRUARI" "MARET" "APRIL" ...
##  $ Tahun : num [1:60] 2014 2014 2014 2014 2014 ...
##  $ Jumlah: num [1:60] 2236 2285 2399 2734 2576 ...
head(data)

Exploratory Data Analysis (EDA)

Analisis data eksploratif (Exploratory Data Analysis/EDA) dilakukan sebagai tahap awal untuk memperoleh gambaran umum mengenai karakteristik data jumlah kunjungan pasien rawat jalan bulanan di Balai Kesehatan Masyarakat selama periode Januari 2014 hingga Desember 2018. Pada tahap ini dilakukan analisis statistik deskriptif serta visualisasi data dalam bentuk grafik deret waktu untuk mengidentifikasi pola perubahan, kecenderungan (tren), dan variasi musiman jumlah kunjungan pasien dari bulan ke bulan.

Sebelum analisis dilakukan, data terlebih dahulu diurutkan berdasarkan tahun dan bulan sehingga dapat dikenali dan diproses sebagai data deret waktu oleh perangkat lunak R.

#statistika deskriptif
library(psych)
describe(data$Jumlah)

Berdasarkan statistik deskriptif, data jumlah kunjungan pasien rawat jalan terdiri dari 60 pengamatan dengan nilai rata-rata sebesar 2623,42 dan simpangan baku 358,15, yang menunjukkan adanya variasi antar bulan. Nilai minimum dan maksimum masing-masing sebesar 1846 dan 3474, sehingga rentang data mencapai 1.628. Nilai skewness sebesar 0,27 mengindikasikan distribusi data sedikit menceng ke kanan, sedangkan nilai kurtosis −0,47 menunjukkan distribusi yang relatif mendatar dan mendekati normal.

#Urutkan data berdasarkan Tahun dan Bulan supaya R dapat mengenali deret waktunya
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
data <- data %>%
  mutate(
    Bulan_urut = match(Bulan, 
     c("JANUARI","FEBRUARI","MARET","APRIL","MEI","JUNI",
      "JULI","AGUSTUS","SEPTEMBER","OKTOBER","NOVEMBER","DESEMBER"))
  )
data <- data %>%
  arrange(Tahun, Bulan_urut)

head(data, 12)
tail(data, 12)
# Definisikan data menjadi time series
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
timeseries_data <- ts(
  as.numeric(data$Jumlah),
  start = c(2014,1),
  frequency = 12
)

timeseries_data
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2014 2236 2285 2399 2734 2576 2609 2091 2358 2513 2441 2190 2207
## 2015 2196 2365 2502 2556 2237 2190 2120 2157 2217 2538 2424 2565
## 2016 2457 2690 2857 2849 2789 2688 2161 3025 2730 2846 2953 2858
## 2017 3064 2710 3097 2903 3170 2233 3257 3281 3056 3474 3445 2793
## 2018 2932 2540 2677 2896 2748 1846 2744 2617 2264 2877 2550 2622
plot(timeseries_data)

Berdasarkan grafik runtun waktu, jumlah kunjungan pasien rawat jalan selama periode Januari 2014 hingga Desember 2018 menunjukkan fluktuasi dari bulan ke bulan dengan kecenderungan meningkat hingga sekitar tahun 2017. Setelah itu, terlihat penurunan pada beberapa periode di tahun 2018. Selain itu, grafik juga memperlihatkan adanya pola musiman yang ditandai dengan naik turunnya jumlah kunjungan pada bulan-bulan tertentu setiap tahunnya. Kondisi ini mengindikasikan bahwa data belum bersifat stasioner pada level dan memerlukan pengujian stasioneritas lebih lanjut sebelum dilakukan pemodelan runtun waktu.

Metode Analisis

Metode analisis yang digunakan dalam penelitian ini adalah analisis runtun waktu (time series) dengan pendekatan Seasonal Autoregressive Integrated Moving Average (SARIMA) (Hariadi, n.d.). Penelitian ini menerapkan pendekatan analisis data runtun waktu dengan menggunakan model Seasonal Autoregressive Integrated Moving Average (SARIMA). Pemilihan model SARIMA didasarkan pada kemampuannya dalam memodelkan serta memprediksi data univariat yang menunjukkan adanya komponen musiman (Junaedi et al., 2025), sebagaimana karakteristik data jumlah kunjungan pasien rawat jalan yang diamati secara bulanan. Tahapan analisis dilakukan secara bertahap sebagai berikut.

#stationeritas data
adf.test(timeseries_data)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  timeseries_data
## Dickey-Fuller = -2.1237, Lag order = 3, p-value = 0.5248
## alternative hypothesis: stationary

Berdasarkan hasil uji Augmented Dickey–Fuller (ADF), diperoleh nilai statistik Dickey–Fuller sebesar −2,1237 dengan p-value sebesar 0,5248. Nilai p-value yang lebih besar dari taraf signifikansi 5% menunjukkan bahwa hipotesis nol tidak dapat ditolak, sehingga data deret waktu pada level belum bersifat stasioner. Oleh karena itu, diperlukan proses differencing untuk mencapai kondisi stasioner sebelum dilakukan pemodelan runtun waktu lebih lanjut.

#differencing
ts_diff1 <- diff(timeseries_data, differences = 1)
adf.test(ts_diff1)
## Warning in adf.test(ts_diff1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_diff1
## Dickey-Fuller = -5.0128, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil uji Augmented Dickey–Fuller (ADF) terhadap data setelah dilakukan differencing orde pertama, diperoleh nilai statistik Dickey–Fuller sebesar −5,0128 dengan p-value sebesar 0,01. Nilai p-value yang lebih kecil dari taraf signifikansi 5% menunjukkan bahwa hipotesis nol ditolak, sehingga data deret waktu telah bersifat stasioner dan memenuhi asumsi untuk dilakukan pemodelan runtun waktu lebih lanjut.

#Plot ACF-PAcf
ts_diff1 <- diff(timeseries_data)
acf(ts_diff1,
    main = "ACF Data Setelah Differencing Orde 1")

pacf(ts_diff1,
     main = "PACF Data Setelah Differencing Orde 1")

Berdasarkan plot ACF setelah dilakukan differencing orde pertama, terlihat adanya satu spike yang signifikan pada lag awal, sedangkan lag-lag selanjutnya cenderung berada dalam batas kepercayaan. Pola ini mengindikasikan adanya komponen Moving Average (MA) pada bagian non-musiman.

Sementara itu, plot PACF menunjukkan adanya spike signifikan pada lag awal dan beberapa lag tertentu, serta indikasi spike pada lag musiman. Hal ini mengarah pada keberadaan komponen Autoregressive (AR), baik pada bagian non-musiman maupun musiman.

Dengan demikian, berdasarkan pola ACF dan PACF setelah differencing orde pertama, data diduga mengikuti model SARIMA dengan kombinasi komponen AR dan MA, yang selanjutnya dikonfirmasi melalui pemilihan model terbaik menggunakan prosedur auto.arima.

#Penentuan Best Model menggunakan Auto arima
model_arima <- auto.arima(
  timeseries_data,
  seasonal = TRUE,
  stepwise = FALSE,
  approximation = FALSE
)
summary(model_arima)
## Series: timeseries_data 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.6416  0.3737
## s.e.   0.0950  0.1336
## 
## sigma^2 = 76257:  log likelihood = -415.5
## AIC=837   AICc=837.44   BIC=843.24
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 2.005867 269.1554 202.5969 -0.743098 7.906012 0.5720383
##                     ACF1
## Training set -0.09371856

Model terbaik yang diperoleh dari proses pemodelan adalah SARIMA(0,1,1)(1,0,0)[12], yang mencakup komponen moving average non-musiman orde satu serta komponen autoregressive musiman orde satu dengan periode 12 bulan. Nilai koefisien yang dihasilkan menunjukkan adanya pengaruh komponen acak dan pola musiman dalam data. Selain itu, nilai AIC, AICc, dan BIC yang relatif rendah mengindikasikan bahwa model memiliki tingkat kecocokan yang baik. Berdasarkan ukuran kesalahan peramalan, nilai MAPE sebesar 7,91% menunjukkan bahwa model SARIMA yang digunakan memiliki tingkat akurasi peramalan yang baik.

#uji diagnostik
checkresiduals(model_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,0,0)[12]
## Q* = 7.4994, df = 10, p-value = 0.6776
## 
## Model df: 2.   Total lags used: 12

Berdasarkan uji Ljung–Box pada residual model ARIMA(0,1,1)(1,0,0)[12] 12, diperoleh p-value sebesar 0,6776 (> 0,05), sehingga tidak terdapat autokorelasi yang signifikan dan residual dapat dianggap sebagai white noise. Hal ini didukung oleh plot diagnostik residual yang menunjukkan fluktuasi di sekitar nol tanpa pola tertentu, nilai ACF residual berada dalam batas kepercayaan, serta distribusi residual yang relatif mendekati normal. Dengan demikian, model ARIMA(0,1,1)(1,0,0)[12] dinyatakan layak digunakan untuk peramalan jumlah kunjungan pasien rawat jalan.

Hasil Dan Pembahasan

Berdasarkan statistik deskriptif, jumlah kunjungan pasien rawat jalan bulanan memiliki rata-rata 2623,42 kunjungan dengan standar deviasi 358,15, yang menunjukkan variasi data relatif moderat. Jumlah kunjungan tertinggi tercatat sebesar 3474 kunjungan pada Oktober 2017, sedangkan terendah sebesar 1.846 kunjungan pada Juni 2018.

Hasil uji Augmented Dickey-Fuller (ADF) menunjukkan data pada level belum stasioner (p-value = 0,5248), sehingga dilakukan differencing satu kali. Setelah differencing orde pertama, data menjadi stasioner secara signifikan (p-value = 0,01).

Model terbaik yang diperoleh adalah SARIMA (0,1,1)(1,0,0)[12]. Uji diagnostik residual menunjukkan residual bersifat white noise dengan p-value uji Ljung–Box sebesar 0,6776, sehingga model dinyatakan layak digunakan untuk peramalan jumlah kunjungan pasien rawat jalan.

#prediksi 12 bulan atau setahun ke depan
forecast_arima <- forecast(model_arima, h = 12)
plot(forecast_arima,
     main = "Prediksi Jumlah Kunjungan Pasien Rawat Jalan",
     xlab = "Tahun",
     ylab = "Jumlah Kunjungan Pasien Rawat Jalan")

Hasil peramalan menunjukkan bahwa jumlah kunjungan pasien rawat jalan selama 12 bulan ke depan diperkirakan berfluktuasi secara musiman di sekitar nilai rata-rata tanpa perubahan tren yang signifikan. Interval prediksi yang semakin melebar menunjukkan meningkatnya ketidakpastian seiring bertambahnya periode peramalan, namun secara umum jumlah kunjungan diperkirakan relatif stabil.

Kesimpulan

Berdasarkan analisis runtun waktu terhadap data jumlah kunjungan pasien rawat jalan bulanan periode Januari 2014 hingga Desember 2018, data menunjukkan adanya pola musiman dan belum stasioner pada level. Setelah dilakukan differencing orde pertama, data menjadi stasioner dan dapat dimodelkan menggunakan metode SARIMA.

Model terbaik yang diperoleh adalah SARIMA (0,1,1)(1,0,0)[12], yang telah memenuhi uji diagnostik dengan residual bersifat white noise dan tingkat akurasi peramalan yang baik. Hasil peramalan menunjukkan bahwa jumlah kunjungan pasien rawat jalan pada 12 bulan ke depan diperkirakan berfluktuasi secara musiman di sekitar nilai rata-rata tanpa perubahan tren yang signifikan. Model ini dapat dimanfaatkan sebagai dasar perencanaan pelayanan dan pengelolaan sumber daya di Balai Kesehatan Masyarakat.

Referensi

Basri, I., Peramalan, M., Kunjungan, J., Rawat, P., Menggunakan, J., Arima, M., Holt-Winters, D., Kesehatan, B., & Xyz, M. (n.d.). Forecasting the Number of Outpatient Patient Visits Using the ARIMA, SES And Holt-Winters Methods at XYZ Community Health Center. Denashurya, N. I. (2024). 46 Pengembangan Model Sarima Untuk Meramalkan Produksi Tanaman Obat di Indonesia. In Jurnal Agribisnis Unisi (Vol. 13, Issue 2). Hariadi, W. (n.d.). Pemodelan Kasus Pasien Terkonfirmasi Positif Covid-19 Per-Hari Di Indonesia dengan Metode SARIMA. Jurnal UJMC, 7(2), 19–30. Junaedi, L., Damastuti, N., & Widodo, A. (2025). Penerapan Metode Seasonal ARIMA (SARIMA) untuk Peramalan Penjualan Barang dengan Pola Musiman Tahunan. JISEM Jurnal Program Studi Informatika Universitas Katolik Widya Mandala Surabaya, 01, 38–48. https://doi.org/10.33508/jisem.v1i01.7403 Rustam, M. Z. A., Amalia, N., & Riestiyowati, M. A. (2022). Analisis Prediksi Kunjungan Pasien Dengan Metode Autoregresiive Integrated Moving Average di Rumah Sakit Ibu dan Anak Putri Surabaya. Jurnal Manajemen Informasi Kesehatan Indonesia, 10(2), 135. https://doi.org/10.33560/jmiki.v10i2.441