Analisis Perbandingan Clustering Harga Cabai Rawit Merah di Jawa Barat Tahun 2022–2025 Menggunakan DTW+K-Medoids dan Pendekatan Tsfeatures

Pengenalan Kasus

Harga cabe rawit merah merupakan salah satu komoditas hortikultura yang memiliki volatilitas tinggi dan memberikan kontribusi signifikan terhadap inflasi pangan di Indonesia, termasuk di Provinsi Jawa Barat. Fluktuasi harga yang terjadi dipengaruhi oleh berbagai faktor seperti kondisi pasokan, pola musim tanam dan panen, distribusi antarwilayah, cuaca ekstrem, serta dinamika permintaan rumah tangga. Ketidakstabilan harga ini tidak hanya berdampak pada kesejahteraan petani dan pedagang, tetapi juga memengaruhi stabilitas ekonomi daerah karena cabe rawit termasuk dalam komoditas volatile food yang sensitif terhadap guncangan pasokan. Oleh karena itu, pemahaman mengenai pola harga cabe rawit antar kabupaten/kota menjadi penting untuk mendukung kebijakan pengendalian harga dan perencanaan logistik pangan.

Tujuan

Mengelompokkan kabupaten/kota di Jawa Barat berdasarkan kemiripan pola harga cabe rawit merah selama tahun 2022–2025.

Data dan Peubah yang digunakan

Data yang digunakan dalam penelitian ini merupakan data harga bulanan cabe rawit merah di seluruh kabupaten/kota di Provinsi Jawa Barat. Rentang waktu data yang tersedia adalah dari tahun 2022 hingga 2025.

Input Data

library(readxl)
data <- read.csv("C:\\Users\\ASUS\\Downloads\\0_df_long.csv")
head(data)
##         KabKot Tahun    Bulan            Produk Harga
## 1 Kab. Bandung  2022  Januari Cabai Rawit Merah 51667
## 2 Kab. Bandung  2022 Februari Cabai Rawit Merah 38259
## 3 Kab. Bandung  2022    Maret Cabai Rawit Merah 53621
## 4 Kab. Bandung  2022    April Cabai Rawit Merah 40967
## 5 Kab. Bandung  2022      Mei Cabai Rawit Merah 40667
## 6 Kab. Bandung  2022     Juni Cabai Rawit Merah 81111
str(data)
## 'data.frame':    1242 obs. of  5 variables:
##  $ KabKot: chr  "Kab. Bandung" "Kab. Bandung" "Kab. Bandung" "Kab. Bandung" ...
##  $ Tahun : int  2022 2022 2022 2022 2022 2022 2022 2022 2022 2022 ...
##  $ Bulan : chr  "Januari" "Februari" "Maret" "April" ...
##  $ Produk: chr  "Cabai Rawit Merah" "Cabai Rawit Merah" "Cabai Rawit Merah" "Cabai Rawit Merah" ...
##  $ Harga : num  51667 38259 53621 40967 40667 ...

Tahap diatas dilakukan untuk imput data dan mengecek strukturnya.

Praprosesing data

#Mengubah bulan ke angka dan membuat kolom waktu
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
data2 <- data %>%
  mutate(
    Bulan = match(Bulan, c("Januari","Februari","Maret","April","Mei","Juni",
                           "Juli","Agustus","September","Oktober","November","Desember")),
    Waktu = paste(Tahun, sprintf("%02d", Bulan), sep = "-")
  )
head(data2)
##         KabKot Tahun Bulan            Produk Harga   Waktu
## 1 Kab. Bandung  2022     1 Cabai Rawit Merah 51667 2022-01
## 2 Kab. Bandung  2022     2 Cabai Rawit Merah 38259 2022-02
## 3 Kab. Bandung  2022     3 Cabai Rawit Merah 53621 2022-03
## 4 Kab. Bandung  2022     4 Cabai Rawit Merah 40967 2022-04
## 5 Kab. Bandung  2022     5 Cabai Rawit Merah 40667 2022-05
## 6 Kab. Bandung  2022     6 Cabai Rawit Merah 81111 2022-06
library(tidyr)

# Mengambil rata-rata per  Kab/Kota per bulan bila ada banyak baris dalam 1 bulan
data2 <- data2 %>%
  group_by(KabKot, Waktu) %>%
  summarise(Harga = mean(Harga, na.rm = TRUE), .groups = "drop")

# Mengubah ke format wide agar K-means dapar membaca tiap baris=1 objek tiap kolom 1 fitur

data_wide <- data2 %>%
  pivot_wider(names_from = Waktu, values_from = Harga)
data_wide
## # A tibble: 27 × 47
##    KabKot  `2022-01` `2022-02` `2022-03` `2022-04` `2022-05` `2022-06` `2022-07`
##    <chr>       <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 Kab. B…     51667     38259     53621     40967     40667     81111     81379
##  2 Kab. B…     60473     39583     58785     40478     41544     84845     81699
##  3 Kab. B…     65000     48929     67944     49286     49653     99792     91667
##  4 Kab. B…     60103     44036     67655     40233     43571     97640     91607
##  5 Kab. C…     72334     46852     59167     53500     53833     84846    100968
##  6 Kab. C…     56484     45190     63672     40655     42957     90630     85323
##  7 Kab. C…     67194     48786     67241     55000     53467     93483    101467
##  8 Kab. G…     49634     35095     51378     38967     39444     78286     77862
##  9 Kab. I…     53600     38714     56806     40833     45871     80345     83871
## 10 Kab. K…     77679     43464     57933     47167     50429     86793     91833
## # ℹ 17 more rows
## # ℹ 39 more variables: `2022-08` <dbl>, `2022-09` <dbl>, `2022-10` <dbl>,
## #   `2022-11` <dbl>, `2022-12` <dbl>, `2023-01` <dbl>, `2023-02` <dbl>,
## #   `2023-03` <dbl>, `2023-04` <dbl>, `2023-05` <dbl>, `2023-06` <dbl>,
## #   `2023-07` <dbl>, `2023-08` <dbl>, `2023-09` <dbl>, `2023-10` <dbl>,
## #   `2023-11` <dbl>, `2023-12` <dbl>, `2024-01` <dbl>, `2024-02` <dbl>,
## #   `2024-03` <dbl>, `2024-04` <dbl>, `2024-05` <dbl>, `2024-06` <dbl>, …
# Cek ukuran data
dim(data_wide)
## [1] 27 47
# Cek missing value
sum(is.na(data_wide))
## [1] 0

Dimensi data terdiri dari 27 baris (Kab/KOta) dan 47 kolom (banyaknya bulan Januari 2022-Oktober 2025).  Tidak ditemukan missing value pada data sehingga tak perlu penanganan khusus.

Eksplorasi

Ringkasan Statistik Harga

summary(as.vector(as.matrix(data_wide[, -1])))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25367   41370   51715   55639   65407  105000

Berdasarkan ringkasan ini dapat diketahui harga minimum cabai rawit merah sepanjang periode pengamatan sebesar Rp25.376 sedangkan maksimumnya adalah Rp105.000, hal ini menunjukkan adanya lonjakan harga yang extrim. Median harga berada pada RP51.715, relatif dekat dengan nilai rata-rata sebesar Rp55.639, yang mengindikasikan distribusi harga cenderung sedikit menceng ke kanan akibat beberapa periode dengan harga sangat tinggi. Kuartil bawah (Q1) sebesar 41.370 dan kuartil atas (Q3) sebesar 65.407 menunjukkan bahwa sebagian besar harga berada pada rentang Rp41.000–65.000 per kilogram.

Plot harga cabai perbulan diJawa Barat

rata_bulanan <- colMeans(data_wide[, -1], na.rm = TRUE)

plot(rata_bulanan, type = "l", main = "Rata-rata Harga Cabe Rawit Merah Jawa Barat",
     xlab = "Waktu", ylab = "Harga")

Berdasarkan plot diatas terlihat bahwa rata-rata harga cabai rawit merah di Jawa Barat mengalami fluktuasi kuat tanpa tren naik atau turun yang jelas. Polanya didominasi oleh lonjakan harga yang terjadi berulang, sehingga mengindikasikan adanya pola musiman atau siklikal yang dipengaruhi kondisi pasokan dan permintaan pada periode tertentu. Variabilitas harga juga cukup tinggi, terlihat dari pergerakan harga yang naik dan turun secara tajam dari waktu ke waktu. Secara keseluruhan, data mencerminkan pasar yang sangat volatil, dengan perubahan harga yang sensitif terhadap perubahan musim maupun gangguan pasokan.

Boxplot sebaran harga antar wilayah

# Boxplot harga bulanan perkabupaten 
boxplot(t(data_wide[, -1]), outline = TRUE,
        main = "Sebaran Harga Bulanan Antar Kabupaten",
        ylab = "Harga")

Urutan wilayah Kab/Kota Kab. Bandung
Kab. Bandung Barat
Kab. Bekasi
Kab. Bogor
Kab. Ciamis
Kab. Cianjur
Kab. Cirebon
Kab. Garut
Kab. Indramayu
Kab. Karawang
Kab. Kuningan
Kab. Majalengka
Kab. Pangandaran
Kab. Purwakarta
Kab. Subang
Kab. Sukabumi
Kab. Sumedang
Kab. Tasikmalaya
Kota Bandung
Kota Banjar
Kota Bekasi
Kota Bogor
Kota Cimahi
Kota Cirebon
Kota Depok
Kota Sukabumi
Kota Tasikmalaya

Terlihat hampir semua kabupaten memiliki rentang harga yang cukup lebar, menunjukkan fluktuasi harga antar bulan yang cukup besar di dalam satu daerah. Median antar kabupaten relatif berbeda beberapa kabupaten memiliki median harga lebih tinggi (misalnya titik sekitar 7–11 dan 18–23), sedangkan yang lain lebih rendah. Banyaknya outlier ke arah atas menunjukkan bahwa di sejumlah bulan terjadi lonjakan harga yang tidak umum di beberapa wilayah. Secara keseluruhan, pola ini memperlihatkan bahwa perbedaan tingkat harga antar kabupaten cukup nyata, dan stabilitas harga dalam satu kabupaten tidak sama: ada yang fluktuasinya tinggi, ada yang relatif sempit.

Heatmap Pola Harga cabai

mat <- as.matrix(data_wide[, -1])
rownames(mat) <- data_wide$KabKot
mat_scaled <- scale(mat)

library(pheatmap)
## Warning: package 'pheatmap' was built under R version 4.5.2
pheatmap(
  mat_scaled,
  scale = "none",
  cluster_rows = FALSE,   # dendogram baris hilang
  cluster_cols = FALSE,   # dendogram kolom hilang
  show_rownames = FALSE,
  show_colnames = FALSE,
  main = "Heatmap Pola Harga Cabe Rawit Merah Jawa Barat (2022–2025)"
)

Secara visual, terlihat adanya blok-blok warna yang cenderung berkumpul, yang mengisyaratkan bahwa beberapa kabupaten memiliki perilaku harga yang mirip pada periode tertentu. Misalnya, beberapa baris memperlihatkan pola kemunculan warna merah atau biru yang berdekatan dalam rentang waktu yang sama, menandakan bahwa kabupaten-kabupaten tersebut kemungkinan berada dalam klaster wilayah dengan dinamika harga yang serupa, entah karena berada di sentra produksi, memiliki kondisi distribusi yang mirip, atau sama-sama rentan terhadap periode kelangkaan.

Pada sisi waktu, pola warna juga tidak acak sepenuhnya. Terdapat bagian memanjang secara vertikal yang memperlihatkan warna yang lebih homogen (misalnya deretan kuning atau pecahan merah dalam beberapa kolom berdekatan), yang menunjukkan adanya klaster waktu atau bulan di mana pergerakan harga relatif konsisten antar wilayah. Ini bisa merupakan periode musiman seperti panen raya, musim hujan, atau gelombang permintaan tertentu, yang menyebabkan harga di banyak kabupaten bergerak ke arah yang sama pada bulan-bulan tersebut.

Merah = Harga lebih tinggi dari rata-rata periode tersebut
Biru = Harga lebih rendah dari rata-rata periode tersebut

DTW+K-medoids

Penyesuaian data dan standarisasi

#Melakukan standarisasi data
# Ambil hanya kolom harga
data_ts <- data_wide[, -1]
# Standarisasi per wilayah (row-wise)
kep_scaled<- t(scale(t(data_ts)))
# Mengubah data menjadi list time series yang merupakan format dari DTW
ts_list <- lapply(1:nrow(kep_scaled), function(i) {
  as.numeric(kep_scaled[i, ])
})

# Beri nama wilayah
names(ts_list) <- data_wide$KabKot

Sintaks diatas mengubah matrix standar menjadi list time series untuk DTW. Setiap kabupaten menjadi satu elemen list, tiap elemen berisi vektor harga bulanan yang telah distandarisasi. Kemudian nama kabupaten disimpan untuk memudahkan interpretasi hasil clustering.

Penentuan jumlah k optimal

#Menentukan jumlah k untuk  dengan metide elbow dan silhoutte
library(cluster)
## Warning: package 'cluster' was built under R version 4.5.2
## Elbow
diss <- c()

for (k in 1:10) {
  pam_k <- pam(kep_scaled, k = k)
  diss[k] <- pam_k$objective
}
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
## Warning in diss[k] <- pam_k$objective: number of items to replace is not a
## multiple of replacement length
plot(1:10, diss, type = "b", pch = 19,
     xlab = "Jumlah Cluster (k)",
     ylab = "Total Dissimilarity",
     main = "Elbow Method (PAM)")

## silhoutte
library(cluster)

avg_sil <- c()

for (k in 2:10) {
  pam_k <- pam(kep_scaled, k = k)
  avg_sil[k] <- pam_k$silinfo$avg.width
}

plot(2:10, avg_sil[2:10], type = "b", pch = 19,
     xlab = "Jumlah Cluster (k)",
     ylab = "Rata-rata Silhouette Width",
     main = "Silhouette untuk K-Medoids (PAM)")

Penentuan jumlah cluster optimal dilakukan menggunakan metode Silhouette dan Elbow. Hasil Silhouette menunjukkan nilai optimal pada k = 2, namun nilai silhouette untuk k = 3 masih berada pada rentang yang dapat diterima. Sementara itu, grafik Elbow menunjukkan bahwa penurunan total dissimilarity masih signifikan hingga k = 3 dan cenderung melandai setelahnya. Oleh karena itu, jumlah cluster k = 3 dipilih sebagai kompromi antara kualitas struktur cluster dan kebutuhan interpretasi yang lebih informatif. Meskipun Silhouette terbaik berada pada k = 2, grafik Elbow menunjukkan bahwa struktur data masih membaik hingga k = 3. Karena silhouette k = 3 masih positif dan stabil, maka k = 3 dipilih sebagai jumlah cluster yang representatif dan lebih interpretatif.

Proses analisis DTW + K-medoids

library(dtwclust)
## Warning: package 'dtwclust' was built under R version 4.5.2
## Loading required package: proxy
## 
## Attaching package: 'proxy'
## The following objects are masked from 'package:stats':
## 
##     as.dist, dist
## The following object is masked from 'package:base':
## 
##     as.matrix
## Loading required package: dtw
## Warning: package 'dtw' was built under R version 4.5.2
## Loaded dtw v1.23-1. See ?dtw for help, citation("dtw") for use in publication.
## dtwclust:
## Setting random number generator to L'Ecuyer-CMRG (see RNGkind()).
## To read the included vignettes type: browseVignettes("dtwclust").
## See news(package = "dtwclust") after package updates.
set.seed(123)

# DTW + K-medoids
dtw_pam <- tsclust(ts_list, #list time series untuk DTW yang sebelumnya sudah dibuat
                   type = "partitional", # melakukan clustering secara langsung berdasarkan k yang sudah ditentukan
                   k = 3,                # jumlah cluster optimal
                   distance = "dtw_basic", 
                   centroid = "pam",      # menggunakan medoid
                   seed = 123)

# Simpan hasil cluster
data_wide$Cluster_DTW_PAM <- dtw_pam@cluster

tsclust(ts_list, …) → menjalankan clustering partitional (K-medoids) dengan:
distance = “dtw_basic” → jarak antar time series dihitung dengan DTW
centroid = “pam” → centroid cluster adalah medoid nyata
k = 3 → jumlah cluster yang ingin dibentuk
set.seed(123) → memastikan hasil clustering reproducible
→ menyimpan cluster assignment tiap kabupaten

Jadi data harga cabe tiap kabupaten diubah terlebih dahulu menjadi time series yang distandarisasi agar setiap kabupaten memiliki skala yang sama, sehingga fokus analisis berada pada pola fluktuasi harga. Selanjutnya, digunakan metode Dynamic Time Warping (DTW) untuk mengukur kesamaan pola antara setiap pasangan kabupaten. Keunggulan DTW adalah kemampuannya mengenali kemiripan pola meskipun puncak harga atau fluktuasi terjadi pada waktu yang berbeda. Setelah menghitung jarak antar time series dengan DTW, dilakukan clustering menggunakan K-medoids (Partitioning Around Medoids / PAM). Berbeda dengan K-means yang menggunakan rata-rata, K-medoids memilih satu kabupaten nyata sebagai medoid untuk setiap cluster, sehingga hasil clustering lebih stabil dan robust terhadap kabupaten dengan pola harga unik atau ekstrem. Algoritma secara iteratif menyesuaikan medoid agar total jarak antar kabupaten ke medoid cluster masing-masing minimum, hingga pembagian cluster stabil.

Hasil clustering

# Lihat pembagian cluster
table(data_wide$Cluster_DTW_PAM) #membuat tabel frekuensi untuk masing-masing cluster.
## 
##  1  2  3 
##  7 19  1
split(data_wide$KabKot, data_wide$Cluster_DTW_PAM)#memisahkan nama kabupaten berdasarkan cluster sehingga dapat dilihat kabupaten mana saja yang tergabung dalam tiap cluster.
## $`1`
## [1] "Kab. Bekasi"      "Kab. Ciamis"      "Kab. Karawang"    "Kab. Subang"     
## [5] "Kab. Tasikmalaya" "Kota Cirebon"     "Kota Sukabumi"   
## 
## $`2`
##  [1] "Kab. Bandung"       "Kab. Bandung Barat" "Kab. Bogor"        
##  [4] "Kab. Cianjur"       "Kab. Cirebon"       "Kab. Garut"        
##  [7] "Kab. Indramayu"     "Kab. Kuningan"      "Kab. Pangandaran"  
## [10] "Kab. Purwakarta"    "Kab. Sukabumi"      "Kab. Sumedang"     
## [13] "Kota Bandung"       "Kota Banjar"        "Kota Bekasi"       
## [16] "Kota Bogor"         "Kota Cimahi"        "Kota Depok"        
## [19] "Kota Tasikmalaya"  
## 
## $`3`
## [1] "Kab. Majalengka"
# Plot centroid (medoid) dan series
plot(dtw_pam, type = "centroids", main = "Centroid / Medoid per Cluster")
## Warning in ggplot2::geom_line(data = dfcm[dfcm$cl %in% clus, ], ...): Ignoring
## unknown parameters: `main`

plot(dtw_pam, type = "series", main = "Time Series per Cluster")

Jumlah cluster 1 sebanyak cluster 2 sebanyak 15, dan cluster 3 terdiri dari1 daerah saja. Cluster 1 dan 2 terlihat mirip namun tetap terdapat perbedaannya.

Visualisasi

library(ggplot2)
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 4.5.2
# Ambil data yang sudah distandarisasi
ts_matrix <- kep_scaled  # baris = kabupaten, kolom = bulan

# PCA (data sudah distandarisasi)
pca_res <- prcomp(ts_matrix, scale. = FALSE)

# Buat dataframe untuk plotting
pca_df <- data.frame(
  PC1 = pca_res$x[,1],
  PC2 = pca_res$x[,2],
  KabKot = data_wide$KabKot,
  Cluster = factor(data_wide$Cluster_DTW_PAM)
)

# Plot PCA dengan label rapi
ggplot(pca_df, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(size = 3) +
  geom_text_repel(aes(label = KabKot), size = 3) +
  labs(title = "PCA Plot DTW + K-medoids (k=3)",
       x = "PC1", y = "PC2") +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

Alur kerja PCA:
Input data: Data standarisasi tiap kabupaten (kep_scaled), di mana baris = kabupaten, kolom = bulan. Standarisasi memastikan semua kabupaten memiliki skala yang sama, sehingga PCA fokus pada pola fluktuasi harga.
Reduksi dimensi: PCA mengubah data berdimensi tinggi (47 bulan) menjadi beberapa komponen utama (PC) yang menangkap sebagian besar variasi pola harga.
PC1 menangkap variasi terbesar, PC2 variasi terbesar berikutnya, dan seterusnya.
Proyeksi kabupaten: Setiap kabupaten dipetakan ke ruang 2D (PC1 vs PC2). Kabupaten dengan pola harga mirip akan berdekatan di plot.
Visualisasi hasil clustering: Dengan menambahkan warna berdasarkan cluster DTW + K-medoids, dapat dilihat apakah kabupaten dalam cluster yang sama memang memiliki pola harga yang serupa.

Jadi PCA ini membantu kita untuk melihat visualisasi hasil DTW+K-medoids tadi.

interpretasi hasil clustering:
Cluster 1 : Kabupaten/kota ini memiliki pola harga bulanan yang relatif mirip satu sama lain, sehingga digabung menjadi satu cluster
Cluster 2 : Memiliki pola fluktuasi harga mirip, meskipun level absolut mungkin berbeda. Cluster ini lebih besar karena ada lebih banyak kabupaten dengan pola serupa.
Cluster 3: Kab. Majalengka menjadi 1 cluster tersendiri hal ini bisa disebabkan oleh puncak harga yang terjadi di bulan berbeda. Fluktuasi lebih tajam atau lebih datar dibanding kabupaten lain.
DTW menilai shape / pola, sehingga kabupaten dengan pola unik akan membentuk cluster tersendiri. Hal ini didukung oleh peta heatmap pada eksplorasi awal yang menunjukkan Kab.Majalengka memang memiliki pola yang berbeda dibandingkan dengan Kab/Kota lainnya, terlihat pada Kab.Majalengka memiliki fluktuasi harga yang cukup tinggi ditandai warna yang merah dan orange.

DTW + K-medoids menekankan kemiripan pola (shape), bukan level harga mutlak. Jika sebuah kabupaten memiliki fluktuasi harga unik atau pergeseran waktu puncak berbeda dari kabupaten lain, jarak DTW ke cluster manapun cukup besar. Karena K-medoids memilih medoid nyata untuk meminimalkan jarak total dalam cluster, algoritma mungkin menempatkan Majalengka di cluster sendiri agar total jarak antar anggota cluster lain tetap minimum.

Karakteristik tiap cluster

# Hitung rata-rata harga per bulan untuk tiap cluster dan membandingkan dengan rata-rata secara keseluruhan
library(dplyr)

# Hitung rata-rata keseluruhan semua harga
overall_mean <- mean(as.numeric(unlist(data_wide[,-1])), na.rm = TRUE)

# Hitung rata-rata tiap cluster
data_cluster <- data_wide %>%
  mutate(Cluster = factor(Cluster_DTW_PAM)) %>%
  group_by(Cluster) %>%
  summarise(across(-c(KabKot, Cluster_DTW_PAM), mean, na.rm = TRUE)) %>%
  rowwise() %>%
  mutate(
    MeanPriceCluster = mean(c_across(-Cluster)),   # rata-rata semua bulan per cluster
    Category = case_when(
      MeanPriceCluster < overall_mean * 0.95 ~ "Murah",      
      MeanPriceCluster > overall_mean * 1.05 ~ "Mahal",       
      TRUE ~ "Sedang"                                        
    )
  ) %>%
  ungroup() %>%
  dplyr::select(Cluster, MeanPriceCluster, Category)
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(-c(KabKot, Cluster_DTW_PAM), mean, na.rm = TRUE)`.
## ℹ In group 1: `Cluster = 1`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))
data_cluster
## # A tibble: 3 × 3
##   Cluster MeanPriceCluster Category
##   <fct>              <dbl> <chr>   
## 1 1                 56230. Sedang  
## 2 2                 54844. Sedang  
## 3 3                 66599. Mahal

Terlihat bahwa cluster 1 dan 2 memiliki nilai harga cabai yang masuk dalam kategori sedang, hal ini menunjukkan harga pada cluster 1 dan 2 cenderung mirip, namun perlu diingat DTW+K-medoids memisahkan keduanya karena algoritma mempertimbangkan pola fluktuasi harga per bulan, sehingga perbedaan temporal kecil tapi konsisten cukup untuk membedakan cluster. Ini wajar dan menunjukkan bahwa clustering berbasis time series menangkap dinamika musiman dan perubahan harga, bukan hanya level rata-rata. Untuk cluster 3 yang terdiri dari Kab.Majalengka harga cabe secara konsisten lebih tinggi dibanding kabupaten lain, sehingga algoritma DTW + K-medoids menempatkannya dalam cluster tersendiri. Tidak hanya harga rata-rata, pola fluktuasi bulanan Majalengka juga berbeda: puncak harga lebih ekstrem dan timing puncaknya berbeda dari cluster lain yang mana hal ini dipertimbangkan oleh DTW.

#plot visual
library(ggplot2)
library(dplyr)
library(tidyr)

# Pastikan semua kolom harga numeric
data_numeric <- data_wide %>%
  mutate(across(-c(KabKot, Cluster_DTW_PAM), ~as.numeric(.)))

# Data long format
data_long <- data_numeric %>%
  pivot_longer(-c(KabKot, Cluster_DTW_PAM), names_to="Month", values_to="Price") %>%
  group_by(KabKot) %>%
  mutate(MonthNum = row_number()) %>%
  ungroup() %>%
  mutate(Cluster = factor(Cluster_DTW_PAM))

# Hitung rata-rata harga per cluster
cluster_avg <- data_long %>%
  group_by(Cluster, MonthNum) %>%
  summarise(AvgPrice = mean(Price, na.rm = TRUE), .groups = "drop")

# Plot semua cluster termasuk Majalengka (Cluster 3) dengan warna berbeda
ggplot(cluster_avg, aes(x = MonthNum, y = AvgPrice, color = Cluster, group = Cluster)) +
  geom_line(size = 1.5) +
  labs(
    title = "Rata-rata Harga Bulanan per Cluster ",
    x = "Bulan ke-",
    y = "Harga",
    color = "Cluster"
  ) +
  scale_color_manual(values = c("1" = "orange", "2" = "purple", "3" = "red")) +
  scale_x_continuous(breaks = seq(1, max(cluster_avg$MonthNum), by = 2)) +
  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.

dapat dilihat berdasarkan plot tersebut terlihat bahwa ke-3 cluster memiliki bentuk pola yang mirip namum berbeda dan DTW menangkap perbedaan tersebut sehingga menjadikannya cluster tersendiri.

library(dplyr)

data_cluster_stats <- data_wide %>%
  group_by(Cluster = factor(Cluster_DTW_PAM)) %>%
  summarise(
    MeanPrice = mean(as.numeric(unlist(select(cur_data(), -KabKot, -Cluster_DTW_PAM))), na.rm = TRUE),
    SDPrice = sd(as.numeric(unlist(select(cur_data(), -KabKot, -Cluster_DTW_PAM))), na.rm = TRUE)
  )
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `MeanPrice = mean(...)`.
## ℹ In group 1: `Cluster = 1`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
data_cluster_stats 
## # A tibble: 3 × 3
##   Cluster MeanPrice SDPrice
##   <fct>       <dbl>   <dbl>
## 1 1          56230.  17067.
## 2 2          54844.  17838.
## 3 3          66599.  17469.

Simpulan:
Untuk Cluster 1 & 2, rata-rata harga per bulan mirip, dan fluktuasi (standar deviasi) juga relatif sama. Namun, meskipun level harga mirip, pola bulanan tiap kabupaten dalam cluster sedikit berbeda — misalnya timing puncak harga atau tren naik-turunnya berbeda antar kabupaten. Inilah yang diperhatikan oleh DTW, sehingga kedua cluster tetap terpisah.
Untuk Cluster 3 (Majalengka), rata-rata harga per bulan tertinggi dibanding cluster lain. Karena cluster ini hanya berisi satu kabupaten, pola bulanan yang diamati adalah rata-rata harga per bulan di Majalengka itu sendiri, sehingga membentuk cluster tersendiri.

Peta sebaran

library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(ggplot2)
library(dplyr)

# Load shapefile kabupaten/kota Jawa Barat
jabar_kab <- st_read("C:\\Users\\ASUS\\Downloads\\Peta.json")
## Reading layer `Kab_Kota' from data source `C:\Users\ASUS\Downloads\Peta.json' using driver `TopoJSON'
## Simple feature collection with 27 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3704 ymin: -7.820982 xmax: 108.8468 ymax: -5.806538
## CRS:           NA
head(jabar_kab)
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.8977 ymin: -7.450505 xmax: 108.8468 ymax: -5.918429
## CRS:           NA
##      id KODE_KK KODE_PROV         KAB_KOTA   PROVINSI
## 1 32.16   32.16        32           Bekasi Jawa Barat
## 2 32.75   32.75        32      Kota Bekasi Jawa Barat
## 3 32.74   32.74        32     Kota Cirebon Jawa Barat
## 4 32.78   32.78        32 Kota Tasikmalaya Jawa Barat
## 5 32.09   32.09        32          Cirebon Jawa Barat
## 6 32.04   32.04        32          Bandung Jawa Barat
##                         geometry
## 1 MULTIPOLYGON (((106.991 -6....
## 2 MULTIPOLYGON (((106.9186 -6...
## 3 MULTIPOLYGON (((108.5621 -6...
## 4 MULTIPOLYGON (((108.1933 -7...
## 5 MULTIPOLYGON (((108.398 -6....
## 6 MULTIPOLYGON (((107.4768 -7...
library(dplyr)
#menyesuaikan nama daerah anatar peta json dan di data
data_wides <- data_wide %>%
  mutate(KabKot_clean = ifelse(grepl("^Kab\\.\\s*", KabKot),
                               gsub("^Kab\\.\\s*", "", KabKot),  # hapus "Kab. "
                               KabKot))  

data_wides$KabKot_clean <- as.character(data_wides$KabKot_clean)
jabar_kab$KAB_KOTA <- as.character(jabar_kab$KAB_KOTA)

# melakukan joining dengan peta json
jabar_cluster <- left_join(jabar_kab, data_wides, by = c("KAB_KOTA" = "KabKot_clean"))
jabar_cluster
## Simple feature collection with 27 features and 53 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3704 ymin: -7.820982 xmax: 108.8468 ymax: -5.806538
## CRS:           NA
## First 10 features:
##       id KODE_KK KODE_PROV         KAB_KOTA   PROVINSI           KabKot 2022-01
## 1  32.16   32.16        32           Bekasi Jawa Barat      Kab. Bekasi   65000
## 2  32.75   32.75        32      Kota Bekasi Jawa Barat      Kota Bekasi   59833
## 3  32.74   32.74        32     Kota Cirebon Jawa Barat     Kota Cirebon   66677
## 4  32.78   32.78        32 Kota Tasikmalaya Jawa Barat Kota Tasikmalaya   52034
## 5  32.09   32.09        32          Cirebon Jawa Barat     Kab. Cirebon   67194
## 6  32.04   32.04        32          Bandung Jawa Barat     Kab. Bandung   51667
## 7  32.71   32.71        32       Kota Bogor Jawa Barat       Kota Bogor   57741
## 8  32.18   32.18        32      Pangandaran Jawa Barat Kab. Pangandaran   56290
## 9  32.06   32.06        32      Tasikmalaya Jawa Barat Kab. Tasikmalaya   59067
## 10 32.14   32.14        32       Purwakarta Jawa Barat  Kab. Purwakarta   50786
##    2022-02 2022-03 2022-04 2022-05 2022-06 2022-07 2022-08 2022-09 2022-10
## 1    48929   67944   49286   49653   99792   91667   66296   65385   48951
## 2    41429   63500   47100   40290   84036   92679   60172   60828   54484
## 3    38857   58966   53833   48419   86897   97097   58167   62833   61129
## 4    32358   50882   37878   42367   81840   80208   53538   51988   47284
## 5    48786   67241   55000   53467   93483  101467   55793   62567   54645
## 6    38259   53621   40967   40667   81111   81379   48517   51500   53000
## 7    46074   65000   39800   52433   99496   91759   58931   65767   54516
## 8    41607   57581   38267   41613   86923   83226   49333   57233   50345
## 9    35893   55533   45000   44333   81250   92957   57931   53500   49500
## 10   46250   58300   37400   45935   93704   85926   53214   59241   46533
##    2022-11 2022-12 2023-01 2023-02 2023-03 2023-04 2023-05 2023-06 2023-07
## 1    32045   44844   59032   66042   69762   36996   29444   37308   31607
## 2    37321   49828   61750   67343   71683   40074   33397   39967   34500
## 3    43833   45719   61970   69828   78548   50000   41452   43700   37032
## 4    37797   43640   56074   54926   68705   41750   31602   37617   31372
## 5    37567   46938   64281   69071   80935   45517   40097   47552   34645
## 6    39000   41250   61235   55345   63500   38333   33871   37552   32387
## 7    40967   53419   67576   70172   77667   43071   37290   46367   33567
## 8    36333   43226   51875   56786   69097   38448   32903   35333   29033
## 9    49310   50926   58636   58036   68387   49483   40000   47000   45167
## 10   37750   48406   62727   65793   75871   44250   37581   40828   35742
##    2023-08 2023-09 2023-10 2023-11 2023-12 2024-01 2024-02 2024-03 2024-04
## 1    38400   33792   52222   85714   91667   60714   62857   55517   39750
## 2    45862   35429   52917   89241   91103   55651   62393   68865   42292
## 3    43387   36667   56935   89000   86129   60429   60000   54016   41140
## 4    39518   31655   45900   80231   75184   45108   51149   63045   45550
## 5    43677   34893   56774   90833   86833   71742   69310   58656   47593
## 6    37000   33467   42774   79207   78935   55833   53551   71290   43604
## 7    46387   34800   60258   91333   97167   55930   66105   59672   45944
## 8    38867   26633   49500   81000   83548   50968   50500   46563   35179
## 9    44516   38000   37167   84500   79167   45269   48500   55109   48426
## 10   41871   37334   53903   85414   87233   55968   60024   53406   44840
##    2024-05 2024-06 2024-07 2024-08 2024-09 2024-10 2024-11 2024-12 2025-01
## 1    31852   33478   63542   63036   40167   40962   36923   42321   97143
## 2    34707   39309   59941   68042   41481   47976   43432   45339   93051
## 3    34839   39818   48525   63145   38333   40984   37593   35862   85641
## 4    38947   41021   57633   68349   37872   47512   40878   42448   84113
## 5    39226   42500   63968   70871   46133   48323   40833   52643   95217
## 6    35926   37313   56596   58055   36833   45269   37870   45991   86073
## 7    33262   36625   63919   65565   44596   47869   40233   48661   99810
## 8    30667   32679   55484   56613   35667   42581   34138   37931   80068
## 9    43532   42037   57414   65968   38862   42500   33667   40690   78829
## 10   33100   35607   60833   61107   39724   43733   38111   46519   71214
##    2025-02 2025-03 2025-04 2025-05 2025-06 2025-07 2025-08 2025-09 2025-10
## 1    72037   94194   73103   39167   51213   60556   37619   41548   39074
## 2    73632  102243   85358   38687   52129   58202   36784   39084   40426
## 3    71545   89692   64097   40027   50661   60770   41253   41322   41629
## 4    71039   83594   81925   40522   46603   56358   37956   40000   35898
## 5    79452   98373   84345   41280   52888   60822   39129   39596   42419
## 6    67721   87939   76259   37163   48104   54405   32867   34549   34313
## 7    77750  102781   82789   38471   56485   66402   45654   41648   43115
## 8    69939   88274   81190   42613   48874   57849   38849   40500   40624
## 9    68107   87821   71494   33742   48567   57199   36581   36747   34430
## 10   72824   85885   76595   36811   50148   69366   36516   36738   36730
##    Cluster_DTW_PAM                       geometry
## 1                1 MULTIPOLYGON (((106.991 -6....
## 2                2 MULTIPOLYGON (((106.9186 -6...
## 3                1 MULTIPOLYGON (((108.5621 -6...
## 4                2 MULTIPOLYGON (((108.1933 -7...
## 5                2 MULTIPOLYGON (((108.398 -6....
## 6                2 MULTIPOLYGON (((107.4768 -7...
## 7                2 MULTIPOLYGON (((106.7952 -6...
## 8                2 MULTIPOLYGON (((108.3408 -7...
## 9                1 MULTIPOLYGON (((107.9101 -7...
## 10               2 MULTIPOLYGON (((107.3022 -6...
# Cek hasil join
table(is.na(jabar_cluster$Cluster_DTW_PAM))
## 
## FALSE 
##    27
#petaggplot(jabar_cluster) +
 library(ggplot2)

ggplot(jabar_cluster) +
  geom_sf(aes(fill = factor(Cluster_DTW_PAM)), color = "black") +
  scale_fill_manual(values = c("1" = "orange", "2" = "yellow", "3" = "red")) +
  labs(title = "Sebaran Cluster DTW+K-medoids (k=3)",
       fill = "Cluster") +
  theme_minimal()

Peta diatas merupakan peta cluster. Berdasarkan hasil klasterisasi DTW dan pemetaan geografis, diperoleh bahwa wilayah dalam Cluster 1 memiliki tingkat harga yang lebih tinggi dibandingkan Cluster 2, tetapi masih berada di bawah Cluster 3 yang didominasi oleh Kabupaten Majalengka. Secara spasial, wilayah Cluster 2 terkonsentrasi di bagian barat Jawa Barat yang relatif dekat dengan sentra produksi dan pusat distribusi, sehingga harga cenderung lebih rendah. Sebaliknya, wilayah Cluster 1 berada lebih jauh dari pusat produksi utama dan lebih bergantung pada pasokan dari luar, yang mendorong kenaikan harga. Majalengka yang terpisah dalam Cluster 3 menunjukkan karakteristik wilayah dengan akses distribusi yang relatif terbatas serta kondisi topografi yang bervariasi, sehingga membentuk harga tertinggi dibandingkan wilayah lain.

Mengapa Majalengka Memiliki Harga Cabai Yang Mahal??

Kabupaten majalengka memiliki harga cabai yang sangat mahal berbeda dengan daerah lainnya bisa saja diaebabkan oleh beberapa faktor berikut:
1. Majalengka bukanlah sentra cabai. Hal ini membuat Kabupaten Majalengka cenderung harus mendatangkan cabai dari daerah lain, sehingga menaikan biaya operasional dan otomatis meningkatkan harga cabai.
2. Cuaca yang tidak mendukung sehingga membuat cabai tidak tahan lama.
3. Petani lokal banyak mengalami gagal panen akibat hama seperti diaserang oleh kera (Koran Pikiran Rakyat 1 Desember 2025), kemudian diserang oleh hama patek yang menyebabkan gagal panen cabai (Media Indonesia, 2022).
4. Kurangnya perhatian pemerintah dalam melakukan distribusi cabai dan penanganan hama yang terus menerus terjadi.
Faktor faktor ini memungkinkan harga cabai di Majalengka melambung tinggi.

Tsfeature+K-medoids+DTW

Metode tsfeatures → K-Medoids → DTW refine digunakan karena deret waktu memiliki karakter yang kompleks sehingga tidak cukup dianalisis hanya dari satu sudut pandang. Tsfeatures dipakai terlebih dahulu karena mampu mengekstraksi ciri-ciri statistik penting dari setiap deret waktu—seperti tren, musiman, autokorelasi, entropi, hingga tingkat ketidakstabilan—sehingga memberikan representasi global yang ringkas namun informatif.Pada ruang fitur yang lebih sederhana ini, proses clustering menjadi lebih cepat, lebih stabil, dan tidak sensitif terhadap noise dibanding langsung memakai DTW. Setelah fitur terbentuk, digunakan K-Medoids sebagai metode clustering awal karena algoritma ini lebih robust terhadap outlier dibanding KMeans dan memilih pusat cluster dari data asli, sehingga lebih cocok untuk karakter harga/deret waktu yang fluktuatif. Tahap ini menghasilkan struktur cluster kasar berdasarkan sifat statistik global.
Namun, dua deret waktu bisa saja memiliki fitur global mirip tetapi pola bentuk (shape) yang berbeda. Karena itu dilakukan tahap DTW refine, yaitu menghitung jarak DTW antar deret waktu dan menggunakan hasil clustering awal sebagai inisialisasi untuk menyempurnakan pembagian cluster berdasarkan kemiripan bentuk sebenarnya. DTW mampu menangkap pola yang bergeser dalam waktu, memiliki panjang berbeda, atau tidak linier, sehingga pada tahap ini cluster menjadi lebih akurat dan lebih sesuai dengan bentuk dinamika deret waktu. Kombinasi dua pendekatan—fitur statistik dan jarak bentuk—membuat metode ini menghasilkan pengelompokan yang lebih cerdas, stabil, dan interpretatif dibanding memakai satu metode saja.

# Paket yang diperlukan
library(tsfeatures)
## Warning: package 'tsfeatures' was built under R version 4.5.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(dtw)
library(cluster)
library(pbapply)
## Warning: package 'pbapply' was built under R version 4.5.2
library(dplyr)
library(tidyr)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(factoextra)
## Warning: package 'factoextra' was built under R version 4.5.2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(ggrepel)

#jumlah cluster
k <- 3

# Deteksi kolom bulan (format YYYY-MM)
bulan_cols <- grep("^\\d{4}-\\d{2}$", names(data_wides), value = TRUE)
if(length(bulan_cols) < 6) stop("Sepertinya kolom bulan tidak terdeteksi. Pastikan nama kolom format 'YYYY-MM'.")

# Buat list of ts (satu elemen = satu wilayah/kab)
name_col <- if("KabKot_clean" %in% names(data_wides)) "KabKot_clean" else "KabKot"

ts_list <- lapply(seq_len(nrow(data_wides)), function(i){
  vals <- as.numeric(data_wides[i, bulan_cols])
  # buat ts; start diambil dari bulan pertama kolom
  start_year  <- as.integer(substr(bulan_cols[1], 1, 4))
  start_month <- as.integer(substr(bulan_cols[1], 6, 7))
  ts(vals, start = c(start_year, start_month), frequency = 12)
})
names(ts_list) <- data_wides[[name_col]]

Tahap diatas adalah meload library yang dipakai dan mempersiapkan data nya untuk analisis berikutnya.

# Ekstraksi tsfeatures
feat <- tsfeatures(ts_list, features = c("mean", "var", "entropy", "lumpiness", 
                                         "stability", "acf_features", "pacf_features"))

feat_scaled <- scale(feat)

Mengubah deret waktu mentah menjadi ciri numerik, lalu menormalkannya agar bisa dibandingkan dan dikelompokkan secara seimbang.
Fungsi tiap fitur:
mean → rata-rata harga
var → tingkat variasi / volatilitas harga
entropy → tingkat ketidakpastian / kompleksitas pola
lumpiness → seberapa tidak stabil varians antar waktu
stability → kestabilan mean dari waktu ke waktu
acf_features → pola hubungan lag (autokorelasi)
pacf_features → ketergantungan langsung antar periode

# Scaling fitur dan K-Medoids (pam) sebagai clustering awal
set.seed(123)
init_kmed <- pam(feat_scaled, k = k)
init_clusters <- init_kmed$clustering

cat("Cluster awal (fitur):\n")
## Cluster awal (fitur):
print(table(init_clusters)) 
## init_clusters
##  1  2  3 
## 11 12  4

Tahap ini mengerjakan:
1. Standardisasi fitur (agar fair)
2. Clustering awal memakai K-Medoids berbasis fitur tsfeatures
3. Menyimpan output cluster & medoid untuk dipakai sebagai “starting point”
4. Menampilkan distribusi cluster sebelum masuk ke tahap DTW Refinement

#menghitung medoids indeks
medoid_indices <- tapply(1:length(ts_list), init_clusters, function(x) x[1])
medoid_indices <- unlist(medoid_indices)

cat("\nMedoid awal dipakai untuk DTW (index valid):\n")
## 
## Medoid awal dipakai untuk DTW (index valid):
print(medoid_indices)
##  1  2  3 
##  1  3 10

Menentukan titik pusat awal (medoid) untuk algoritma K-medoids berbasis DTW, dengan memilih satu deret waktu per cluster sebagai wakil awal.
Tahap ini penting karena:
Medoid awal mempengaruhi hasil akhir clustering,
Pemilihan dari cluster awal membuat proses lebih stabil dibanding acak murni.

#Hitung matriks jarak DTW
n <- length(ts_list)
dtw_mat <- matrix(0, nrow = n, ncol = n)
rownames(dtw_mat) <- names(ts_list)
colnames(dtw_mat) <- names(ts_list)

pairs <- which(upper.tri(dtw_mat), arr.ind = TRUE)

pboptions(type = "timer")

compute_one <- function(idx_pair){
  i <- idx_pair[1]; j <- idx_pair[2]
  a <- ts_list[[i]]; b <- ts_list[[j]]
  d <- dtw(a, b, distance.only = TRUE)$distance
  c(i, j, d)
}

res <- pbapply(pairs, 1, compute_one)

for(col in seq_len(ncol(res))){
  i <- res[1,col]; j <- res[2,col]; d <- res[3,col]
  dtw_mat[i,j] <- dtw_mat[j,i] <- d
}

diag(dtw_mat) <- 0
dtw_dist <- as.dist(dtw_mat)
# Ambil medoids awal dari clustering fitur
init_medoids <- init_kmed$medoids   # karakter atau numeric tergantung output
init_clusters <- init_kmed$clustering

#Konversi medoids menjadi indeks numerik valid
if (is.character(init_medoids)) {
  medoid_indices <- match(init_medoids, rownames(dtw_mat))
} else {
  medoid_indices <- as.integer(init_medoids)
}


# Jika ada NA atau jumlah tidak sesuai atau tidak unik → invalid
if (any(is.na(medoid_indices)) ||
    length(medoid_indices) != k ||
    length(unique(medoid_indices)) != k ||
    any(medoid_indices < 1 | medoid_indices > length(ts_list))) {

  valid <- FALSE
}

# Jika tidak valid → ambil medoid aman dari cluster awal
if (!valid) {

  medoid_indices <- tapply(1:length(ts_list), init_clusters, function(x) x[1])
  medoid_indices <- as.integer(unlist(medoid_indices))
}

Tahap diatas memastikan medoid awal yang dipakai dalam DTW valid, unik, dan sesuai jumlah cluster; jika bermasalah, otomatis diganti dengan pilihan aman dari cluster awal.

#FINAL PAM dengan DTW distance
set.seed(123)
final_pam <- pam(
  x = dtw_dist,
  k = k,
  diss = TRUE,
  medoids = medoid_indices
)

final_clusters <- final_pam$clustering

cat("\nRefined K-Medoids (DTW) - jumlah anggota tiap cluster:\n")
## 
## Refined K-Medoids (DTW) - jumlah anggota tiap cluster:
print(table(final_clusters))
## final_clusters
##  1  2  3 
## 15 10  2
#menghitung jumlah anggota tiap cluster
data_wides$cluster_tsfeature_kmedoids_dtw <- final_clusters[
  match(data_wides[[name_col]], names(ts_list))
]

# Ringkasan jumlah anggota tiap cluster
cat("\nJumlah anggota cluster berdasarkan DTW-refined K-Medoids:\n")
## 
## Jumlah anggota cluster berdasarkan DTW-refined K-Medoids:
print(table(data_wides$cluster_tsfeature_kmedoids_dtw, useNA = "ifany"))
## 
##  1  2  3 
## 15 10  2
#Melihat anggota tiap cluster
cluster_members <- split(data_wides[[name_col]],
                         data_wides$cluster_tsfeature_kmedoids_dtw)

cat("\nDaftar wilayah per cluster:\n")
## 
## Daftar wilayah per cluster:
print(cluster_members)
## $`1`
##  [1] "Bandung"          "Bandung Barat"    "Cianjur"          "Garut"           
##  [5] "Indramayu"        "Karawang"         "Pangandaran"      "Purwakarta"      
##  [9] "Subang"           "Sukabumi"         "Sumedang"         "Tasikmalaya"     
## [13] "Kota Banjar"      "Kota Cirebon"     "Kota Tasikmalaya"
## 
## $`2`
##  [1] "Bekasi"        "Bogor"         "Cirebon"       "Kuningan"     
##  [5] "Kota Bandung"  "Kota Bekasi"   "Kota Bogor"    "Kota Cimahi"  
##  [9] "Kota Depok"    "Kota Sukabumi"
## 
## $`3`
## [1] "Ciamis"     "Majalengka"

Terdapat perbedaan yang menarik pada cluster 3 untuk metode ini ciamis bergabung dengan Majalengka.hal ini mungkin saja terjadi sebab:
1. tsfeatures clustering (fitur statistik)
Berdasarkan rangkumannya → mean, var, entropi, acf, pacf
Tidak melihat bentuk kurva asli
Lebih mirip “summary” heuristic
Ciamis & Majalengka terlihat ekstrem → dikelompokkan terpisah

  1. DTW K-Medoids clustering
    Berdasarkan bentuk kurva waktu bulanan
    Sensitif terhadap urutan naik–turun, bukan angka absolut
    Lebih akurat untuk time-series yang musiman
    Majalengka → bentuk terlalu unik → terpisah sendirian
    Ciamis menunjukkan kemiripan pola temporal dengan wilayah seperti Karawang, Subang, dan Tasikmalaya, sehingga dalam pendekatan DTW digabungkan dalam cluster yang sama.

Perbandingan

Hasil clustering menggunakan dua pendekatan—tsfeatures dan DTW K-Medoids—menunjukkan perbedaan yang menarik karena prinsip kerjanya berbeda. Clustering berbasis tsfeatures mengelompokkan kabupaten/kota berdasarkan ciri statistik ringkas dari deret waktunya, seperti rata-rata, varians, entropi, stabilitas, dan pola autokorelasi. Dalam pendekatan ini, wilayah dengan pola statistik ekstrem, seperti Ciamis dan Majalengka, ditempatkan pada cluster tersendiri meskipun bentuk kurva bulanan mereka mungkin mirip dengan wilayah lain. Sementara itu, clustering menggunakan DTW K-Medoids mempertimbangkan bentuk kurva harga bulanan secara langsung, sehingga wilayah yang memiliki pola naik-turun yang sinkron akan berada dalam cluster yang sama meskipun level harga absolutnya berbeda. Akibatnya, Majalengka tetap berdiri sendiri karena bentuk fluktuasinya unik, sedangkan Ciamis bergabung dengan cluster lain yang kurvanya mirip, meskipun secara statistik sebelumnya dianggap ekstrem. Perbedaan ini menegaskan bahwa tsfeatures menekankan ringkasan statistik, sedangkan DTW menekankan kesamaan bentuk kurva, sehingga interpretasi cluster bisa berbeda.