...Ini bukan cerita tentang wilkerst*t, tapi tentang data yang kita sayangi dan cintai



Kalau ada alat yang bisa digunakan untuk melihat-lihat dan mendapatkan insight dari data yang kita punya, exploratory data analysis (EDA) adalah jawabannya. Ia adalah tool yang powerfull dan sangat bermanfaat bagi pengguna data. Lalu, apa itu EDA?

Dalam universe statistik dan analisis data, pemahaman terhadap data merupakan langkah fundamental sebelum dilakukan pemodelan atau pengambilan keputusan. Salah satu pendekatan penting dalam proses ini adalah Exploratory Data Analysis (EDA). Konsep EDA pertama kali diperkenalkan secara formal oleh John Tukey, seorang ahli statistik terkemuka, dalam bukunya “Exploratory Data Analysis” pada tahun 1977.

Pak Tukey bilang:

“Exploratory data analysis is an attitude, a flexibility, and a reliance on display, not a bundle of techniques.”

Pernyataan ini memberi penekanan bahwa EDA bukan sekadar kumpulan metode statistik, tetapi merupakan pendekatan yang bersifat fleksibel dan eksploratif, dengan fokus yang kuat pada visualisasi data. EDA bertujuan untuk menggali informasi tersembunyi, mengidentifikasi pola dan anomali, serta menghasilkan pemahaman intuitif tentang struktur data, semua dilakukan sebelum masuk ke tahap pemodelan matematis.

Karakteristik Unik Data Time Series

Data time series berbeda dari data tabular biasa karena adanya dependensi waktu antar observasi. Hal ini mengharuskan pendekatan EDA yang mempertimbangkan aspek-aspek seperti:

  1. Tren (Trend): Pola umum yang menunjukkan peningkatan atau penurunan jangka panjang.
  2. Musiman (Seasonality): Pola berulang secara periodik dalam interval waktu tertentu, misalnya bulanan atau tahunan.
  3. Siklus (Cyclic): Fluktuasi jangka panjang yang tidak bersifat periodik tetap seperti musiman.
  4. Kejadian tak terduga (Outlier): Nilai ekstrem yang bisa disebabkan oleh faktor eksternal.

Yang Akan dilakukan..

Ada beberapa hal yang akan dilakukan dalam proses EDA ini, antara lain:

  1. Membuat statistik deskriptif
  2. Membuat plot data time series
  3. Melakukan dekomposisi data
  4. Uji stasioner

Nah, data yang akan kita EDA-in adalah data produksi gabah kering giling (GKG) bulanan periode 2018 sampai 2025. Demi menyebarluaskan FOSS, saya akan menggunakan bahasa pemrograman R (free) yang dijalankan di mesin ber-OS Linux Zorin (juga free).

1. Descriptive Statistics

Langkah pertama, kita perlu mengimport dulu file ke RStudio. Panggil dulu library yang akan digunakan, setelah itu import datanya.

# Load library
library(readxl)
library(dplyr)
library(ggplot2)
library(lubridate)
library(tidyr)
library(scales)

# Baca data dari Excel
data <- read_excel("C:/Users/acer/Documents/data.xlsx")

Kita inspect datanya dengan menggunakan str().

str(data)
## tibble [88 × 2] (S3: tbl_df/tbl/data.frame)
##  $ bulan            : POSIXct[1:88], format: "2018-01-01" "2018-02-01" ...
##  $ produksi_padi_ton: num [1:88] 114130 148112 323383 554272 260567 ...

Dari hasil inspeksi data, terdapat dua fitur, yaitu bulan dan produksi_padi_ton. FYI, fitur bulan merupakan indeks waktu bulanan dan fitur produksi_padi_ton merupakan produksi GKG dalam ton. Periode data sebanyak 88 bulan.

Kita buat dulu statistik deskriptifnya.

# Ubah kolom waktu menjadi format tanggal
data <- data %>%
  mutate(bulan = as.Date(bulan))

# Statistik deskriptif produksi padi
summary_stats <- data %>%
  summarise(
    mean = mean(produksi_padi_ton, na.rm = TRUE),
    median = median(produksi_padi_ton, na.rm = TRUE),
    sd = sd(produksi_padi_ton, na.rm = TRUE),
    min = min(produksi_padi_ton, na.rm = TRUE),
    max = max(produksi_padi_ton, na.rm = TRUE),
    q1 = quantile(produksi_padi_ton, 0.25, na.rm = TRUE),
    q3 = quantile(produksi_padi_ton, 0.75, na.rm = TRUE)
  )

print(summary_stats)
## # A tibble: 1 × 7
##      mean  median      sd   min    max     q1      q3
##     <dbl>   <dbl>   <dbl> <dbl>  <dbl>  <dbl>   <dbl>
## 1 219456. 144412. 188801. 19177 797741 71144. 326386.

Statistik lima serangkainya sudah ditampilkan di atas. Sila dibaca dan interpretasi sendiri, ya.

2. Plot

Setelah statistik deskriptif, berikutnya kita membuat visualisasi dalam bentuk line plot.

# Membuat line plot
library(ggplot2)

ggplot(data, aes(x = bulan, y = produksi_padi_ton)) +
  geom_line(color = "steelblue", size = 1) +
  geom_point(color = "darkred", size = 2) +
  scale_y_continuous(labels = comma)+
  labs(
    title = "Produksi GKG Bulanan",
    x = "Bulan",
    y = "Produksi (Ton)"
  ) +
  theme_minimal()

Plot menunjukkan pola yang teratur dari tahun ke tahun dengan periode puncak produksi GKG antara Bulan April atau Mei untuk musim tanam I (MT 1) kemudian menurun dan melandai. Di MT 2, puncak produksi sebagian besar terjadi pada Bulan September - Oktober, setelah itu panen menurun dan kembali melandai sampai akhir tahun. Pola 2 puncak periode panen ini menunjukkan kebiasaan tanam para petani di wilayah tersebut, umumnya petani melakukan 2 kali penanaman dalam 1 tahun.

Kita coba lihat bagaimana jika plot dikelompokkan berdasarkan tahun.

# Tambahkan kolom tahun dan bulan numerik
data <- data %>%
  mutate(
    bulan = as.Date(bulan),
    tahun = year(bulan),
    bulan_num = month(bulan)
  )

# Plot garis berdasarkan tahun
ggplot(data, aes(x = bulan_num, y = produksi_padi_ton, color = factor(tahun))) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  scale_y_continuous(labels = comma)+
  labs(
    title = "Produksi GKG Bulanan per Tahun",
    x = "Bulan",
    y = "Produksi (Ton)",
    color = "Tahun"
  ) +
  theme_minimal()

Plot berdasarkan tahun menguatkan informasi dari sebelumnya, di mana puncak panen pertama terjadi di Bulan April, akan tetapi terjadi pergeseran panen di tahun 2024, di mana puncak panen MT 1 terjadi pada Bulan Mei.

Selanjutnya, boxplot dan violinplot.

# Plot boxplot + violin seluruh data
ggplot(data, aes(x = "", y = produksi_padi_ton)) +
  geom_violin(fill = "lightblue", alpha = 0.5, color = NA) +
  geom_boxplot(width = 0.15, fill = "white", outlier.color = "red") +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Distribusi Produksi GKG",
    x = NULL,
    y = "Produksi (Ton)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_blank())

Bentuk violin plot menunjukkan distribusi data yang cukup menyebar, dengan sebagian besar data berada di kisaran bawah (sekitar 0–300 ribu ton). Ada kepadatan data tinggi di area bawah (violin lebih lebar) dan semakin menyempit ke arah atas, menandakan jumlah observasi berkurang untuk nilai produksi yang besar.

Garis horizontal di dalam boxplot menunjukkan median produksi padi sekitar 150 ribu ton. Artinya, setengah bulan dalam data memiliki produksi ≤ 150 ribu ton, dan setengahnya lagi ≥ 150 ribu ton. Kotak (box) mencakup rentang kuartil 1 (Q1) hingga kuartil 3 (Q3), kira-kira dari 75 ribu ton hingga 320 ribu ton. Ini menunjukkan variasi yang cukup besar antar bulan.

Titik merah di bagian atas adalah outlier, yaitu bulan dengan produksi sangat tinggi dibanding bulan-bulan lain, sekitar ≥ 780 ribu ton. Outlier ini mungkin disebabkan oleh faktor musiman, panen raya, atau kesalahan pencatatan. Who knows?

Whisker (garis vertikal) menunjukkan produksi terendah mendekati nol ton, dan tertinggi (di luar outlier) sekitar 650 ribu ton.

Bagaimana jika diplot per tahun?

# Tambahkan kolom tahun & bulan
# Baca data
data <- read_excel("C:/Users/acer/Documents/data.xlsx") %>%
  mutate(
    bulan = as.Date(bulan),
    tahun = year(bulan),
    bulan_label = month(bulan, label = TRUE, abbr = TRUE)
  )

# Fungsi untuk mendeteksi outlier per tahun
find_outliers <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  x > (Q3 + 1.5 * IQR) | x < (Q1 - 1.5 * IQR)
}

# Tandai outlier
data_outliers <- data %>%
  group_by(tahun) %>%
  filter(find_outliers(produksi_padi_ton)) %>%
  ungroup()

# Plot violin + boxplot + label outlier
ggplot(data, aes(x = factor(tahun), y = produksi_padi_ton)) +
  geom_violin(fill = "orange", alpha = 0.5, color = NA) +
  geom_boxplot(width = 0.15, fill = "white", outlier.color = "red") +
  scale_y_continuous(labels = comma) +
  geom_text(data = data_outliers,
            aes(label = bulan_label),
            vjust = -0.5, color = "blue", size = 3) +
  labs(
    title = "Distribusi Produksi GKG per Tahun",
    x = "Tahun",
    y = "Produksi (Ton)"
  ) +
  theme_minimal()

Distribusi data setiap tahun menunjukkan pola yang relatif serupa dengan data secara keseluruhan, mirip-mirip je. Sedikit catatan untuk data tahun 2025 baru ada data sampai dengan bulan April. Jadi violinnya nampak kurus begitu.

3. Decomposition

Dari asal katanya, decomposition berarti proses menguraikan data menjadi komponen-komponen penyusunnya, yaitu komponen seasonal, trend, dan error atau remainder atau sisaan. Dekomposisi ini sendiri ada yang sifatnya aditif, ada yang multiplikatif, tinggal dilihat saja dari datanya. Kalau datanya terus membesar nilainya, biasanya dia multiplikatif. Kalau datanya relatif datar kayak hidup kau, kemungkinan besar aditif.

# Baca data
data <- read_excel("C:/Users/acer/Documents/data.xlsx") %>%
  mutate(bulan = as.Date(bulan)) %>%
  arrange(bulan)

# Ubah jadi time series
ts_data <- ts(data$produksi_padi_ton,
              start = c(year(min(data$bulan)), month(min(data$bulan))),
              frequency = 12)

# STL decomposition
stl_result <- stl(ts_data, s.window = "periodic")

# Konversi hasil STL ke data frame
stl_df <- data.frame(
  date = data$bulan,
  data = as.numeric(stl_result$time.series[, "remainder"] + 
                      stl_result$time.series[, "trend"] +
                      stl_result$time.series[, "seasonal"]),
  seasonal = as.numeric(stl_result$time.series[, "seasonal"]),
  trend = as.numeric(stl_result$time.series[, "trend"]),
  remainder = as.numeric(stl_result$time.series[, "remainder"])
)

# Ubah ke format long untuk ggplot
stl_long <- stl_df %>%
  pivot_longer(cols = -date, names_to = "component", values_to = "value")

# Plot dengan ggplot
ggplot(stl_long, aes(x = date, y = value)) +
  geom_line(color = "steelblue", size = 1) +
  facet_wrap(~ component, scales = "free_y", ncol = 1) +
  labs(
    title = "STL Decomposition Produksi GKG Bulanan",
    x = "Tahun",
    y = "Ton"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold", size = 12),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5)
  )

Grafik menunjukkan hasil dekomposisi berdasarkan komponen-komponen yang sudah disebutkan di atas. Komponen seasonal menunjukkan bahwa pola tahunan selalu berulang dengan puncak panen MT 1 konsisten di Bulan April - Mei, MT 2 konsisten di Bulan September - Oktober. Bulan Desember - Januari memiliki panen terendah karena pada saat itu biasanya petani baru memulai musim tanamnya. Amplitudo musiman cukup besar, artinya faktor musim sangat berpengaruh terhadap produksi GKG.

Komponen trend menunjukkan pola data dalam jangka panjang. Secara umum, polanya relatif stabil dengan sedikit fluktuatif khas dari pola pertanaman padi, di mana setiap bulan padi mengalami fase mulai dari pengolahan lahan sampai panen. Dari 2018 hingga sekitar 2020, tren sedikit menurun, kemudian naik lagi mulai 2021 hingga 2024. Kenaikan tren akhir-akhir ini menunjukkan kemungkinan peningkatan produksi GKG secara struktural, mungkin karena perbaikan teknologi, irigasi, atau perluasan lahan panen.

Nilai remainder relatif kecil dibandingkan komponen musiman dan tren, tapi ada lonjakan besar pada titik-titik tertentu, kemungkinan akibat kejadian khusus seperti cuaca ekstrem, bencana, atau perubahan kebijakan pertanian. Variasi remainder yang kecil menunjukkan model STL mampu menangkap pola utama data.

4. Stationery Testing

Stasioneritas itu penting di analisis time series karena banyak metode statistik dan machine learning klasik seperti ARIMA, SARIMA, dan VAR secara matematis mengasumsikan bahwa data yang diolah bersifat stasioner. Kalau data tidak stasioner, model bisa salah membaca pola dan hasil prediksi jadi bias.

Data stasioner berarti sifat statistiknya tidak berubah seiring waktu: rata-rata tetap, variansi tetap, dan autokorelasi hanya tergantung pada jarak lag, bukan waktu.

Kita akan coba melihat apakah datanya stasioner atau tidak menggunakan dua uji, yaitu uji ADF (Augmented Dickey-Fuller test) dan uji KPSS (The Kwiatkowski-Phillips-Schmidt-Shin test).

# load package
library(tseries)

# Baca data
df <- read_excel("C:/Users/acer/Documents/data.xlsx")

# Ubah jadi time series bulanan mulai Januari 2018
ts_data <- ts(df$produksi_padi_ton, start = c(2018, 1), frequency = 12)

# Uji ADF
adf_result <- adf.test(ts_data)

print(adf_result)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_data
## Dickey-Fuller = -4.7922, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Pada uji ADF, hipotesis awal atau \(h_0\) adalah data mengandung akar unit (tidak stasioner) dan p-value hasil uji adalah 0.01. Karena p-value lebih kecil dari \(\alpha\) = 0.05, maka \(h_0\) tertolak, artinya data stasioner.

# Uji KPSS
kpss_result <- kpss.test(ts_data)
print(kpss_result)
## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_data
## KPSS Level = 0.096712, Truncation lag parameter = 3, p-value = 0.1

Uji KPSS sedikit berbeda dengan uji ADF, tetapi interpretasinya memberikan hasil yang sama. \(H_0\) pada uji KPSS yaitu data stasioner. P-value hasil uji 0.1, lebih besar dari \(\alpha\) = 0.05. Kesimpulannya? Tidak tolak \(h_0\), kau terima saja \(h_0\), yang artinya data bersifat stasioner.

Sekian ya untuk materi EDA for time series data. Tulisan ini hanya gambaran kicik je apa yang kita lakukan pada tahap pra-pemodelan. Gambaran kicik, karena masih banyak yang bisa dilakukan selain tahapan di atas. Kenapa EDA? Karena perlu dan penting. Kenapa story tally?

Hehehehe, ini hanya guyonan saja karena di kantor kami sedang ada kegiatan perlombaan untuk menyemarakkan momen perayaan HUT kemerdekaan. Salah sebut, harusnya story telling tapi jadi story tally. Mungkin gara-gara bulan ini sedang wilkerst*t, ya?

Next tulisan, kita akan bahas model yang umum digunakan dalam time series. Rencananya akan dilengkapi juga dengan konsep Matematikanya. Stay tuned, ya. Kata pak Mikael Dewabrata, simak dan bookmark!