Tugas Akhir Teknik Peubah Ganda Paralel 1
Abdullah Zeith Alfarizy (G1401231060)
2025-12-02
Penerapan Prophet untuk Peramalan Deret Waktu Harga Beras Medium di Jawa Barat Berdasarkan Klasterisasi Fuzzy C-Medoids Berbasis Dynamic Time Warping
Analisis ini bertujuan mengelompokkan Kabupaten/Kota di Jawa Barat berdasarkan pola perubahan harga beras medium yang terjadi sepanjang periode pengamatan. Pendekatan Fuzzy C-Medoids dipilih karena sifatnya yang mampu menangkap kedekatan parsial suatu wilayah terhadap lebih dari satu klaster, sehingga lebih merepresentasikan realitas ekonomi yang tidak bersifat hitam-putih. Selain itu, ukuran jarak Dynamic Time Warping digunakan untuk mengukur kemiripan pola waktu meskipun terdapat pergeseran fase fluktuasi harga antar wilayah. Hasil akhir dari analisis ini diharapkan dapat menjadi masukan bagi pengambil kebijakan untuk menentukan strategi stabilisasi harga yang lebih adaptif terhadap karakter dinamika pasar tiap kelompok wilayah.
Inisialisasi Workspace dan Package Pendukung
rm(list = ls()); gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 548323 29.3 1223408 65.4 686457 36.7
## Vcells 1018412 7.8 8388608 64.0 1876282 14.4# 1) PACKAGES -----------------------------------------------------------------
packages <- c("readxl","dplyr","tidyr","ggplot2","lubridate","TSclust",
"reshape2","viridis","writexl","zoo","FactoMineR","ca",
"gridExtra","factoextra","cluster","RColorBrewer")
new_pkgs <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_pkgs)) install.packages(new_pkgs, dependencies = TRUE)
library(readxl); library(dplyr); library(tidyr); library(ggplot2)## Warning: package 'readxl' was built under R version 4.4.3## Warning: package 'dplyr' was built under R version 4.4.3##
## 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## Warning: package 'tidyr' was built under R version 4.4.2## Warning: package 'ggplot2' was built under R version 4.4.3library(lubridate); library(TSclust); library(reshape2); library(viridis)## Warning: package 'lubridate' was built under R version 4.4.3##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union## Warning: package 'TSclust' was built under R version 4.4.3## Loading required package: pdc## Warning: package 'pdc' was built under R version 4.4.3## Loading required package: cluster## Warning: package 'cluster' was built under R version 4.4.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Warning: package 'reshape2' was built under R version 4.4.3##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smiths## Warning: package 'viridis' was built under R version 4.4.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.4.2library(writexl); library(zoo); library(FactoMineR); library(ca)## Warning: package 'writexl' was built under R version 4.4.3## Warning: package 'zoo' was built under R version 4.4.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Warning: package 'FactoMineR' was built under R version 4.4.3## Warning: package 'ca' was built under R version 4.4.3library(gridExtra); library(factoextra); library(cluster); library(RColorBrewer)## Warning: package 'gridExtra' was built under R version 4.4.3##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine## Warning: package 'factoextra' was built under R version 4.4.3## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBacat("✓ Packages loaded\n")## ✓ Packages loadedPemanggilan Data
input_file <- "C:\\Users\\FAQIH\\Downloads\\datatpg.xlsx" # sesuaikan
df <- read_excel(input_file)
if(!all(c("KabKot","Tahun","Bulan","Produk","Harga","NamaProv") %in% colnames(df))){
stop("Kolom wajib tidak lengkap. Pastikan ada KabKot, Tahun, Bulan, Produk, Harga, NamaProv.")
}
cat("✓ Data loaded:", nrow(df), "rows\n")## ✓ Data loaded: 1242 rowsPra-Pemrosesan Data Time Series Harga Beras
cat("=== PREPROCESSING (IMPROVED) ===\n")## === PREPROCESSING (IMPROVED) ===# 3.1 Filter Beras Medium & Jawa Barat
df_clean <- df %>%
filter(Produk == "Beras Medium", NamaProv == "Jawa Barat") %>%
mutate(Bulan = as.character(Bulan))
cat("Observations after filter:", nrow(df_clean), "\n")## Observations after filter: 1242# 3.2 Penyelarasan Format Waktu (Map nama bulan Indonesia -> angka)
bulan_map <- c("Januari"=1,"Februari"=2,"Maret"=3,"April"=4,"Mei"=5,"Juni"=6,
"Juli"=7,"Agustus"=8,"September"=9,"Oktober"=10,"November"=11,"Desember"=12)
df_clean <- df_clean %>%
mutate(BulanNum = as.integer(bulan_map[Bulan]),
Periode = paste(Tahun, sprintf("%02d", BulanNum), sep = "-"),
Tanggal = as.Date(paste(Tahun, BulanNum, "01", sep = "-"))) %>%
arrange(KabKot, Tanggal)
# 3.3 Penanganan Duplikasi Data
dups <- df_clean %>% group_by(KabKot, Tahun, BulanNum) %>% filter(n() > 1)
if(nrow(dups) > 0){
cat("WARNING: duplicates found - aggregating by mean\n")
df_clean <- df_clean %>%
group_by(KabKot, Tahun, BulanNum) %>%
summarise(Harga = mean(Harga, na.rm=TRUE),
Produk = first(Produk),
NamaProv = first(NamaProv),
Tanggal = first(Tanggal),
.groups = "drop") %>%
arrange(KabKot, Tanggal)
}
# 3.4 Pembentukan Grid Waktu Lengkap
min_date <- min(df_clean$Tanggal, na.rm = TRUE)
max_date <- max(df_clean$Tanggal, na.rm = TRUE)
all_periods <- seq(min_date, max_date, by = "month")
all_periods_char <- format(all_periods, "%Y-%m")
kab_list <- unique(df_clean$KabKot)
df_expanded <- expand.grid(KabKot = kab_list, Periode = all_periods_char, stringsAsFactors = FALSE)
df_expanded$Periode <- as.character(df_expanded$Periode)
# 3.5 Imputasi Nilai Hilang (LOCF + Interpolasi)
df_ts_prep <- df_clean %>%
mutate(Periode = format(Tanggal, "%Y-%m")) %>%
select(KabKot, Periode, Harga) %>%
right_join(df_expanded, by = c("KabKot","Periode")) %>%
arrange(KabKot, Periode) %>%
group_by(KabKot) %>%
mutate(Harga = as.numeric(Harga),
# STEP 1: Fill leading NAs (forward fill)
Harga = na.locf(Harga, na.rm = FALSE),
# STEP 2: Fill trailing NAs (backward fill)
Harga = na.locf(Harga, fromLast = TRUE, na.rm = FALSE),
# STEP 3: Interpolate middle gaps
Harga = na.approx(Harga, x = as.Date(paste0(Periode, "-01")), na.rm = FALSE)) %>%
ungroup()
remaining_na <- df_ts_prep %>% group_by(KabKot) %>% summarise(n_na = sum(is.na(Harga)))
bad_kab <- remaining_na %>% filter(n_na > 0) %>% pull(KabKot)
if(length(bad_kab) > 0){
cat("Removing kab/kota with remaining NA:", paste(bad_kab, collapse = ", "), "\n")
df_ts_prep <- df_ts_prep %>% filter(!(KabKot %in% bad_kab))
kab_list <- unique(df_ts_prep$KabKot)
}
df_wide <- df_ts_prep %>% pivot_wider(names_from = Periode, values_from = Harga) %>% arrange(KabKot)
kabkota_names <- df_wide$KabKot
ts_matrix <- as.matrix(df_wide[, -1])
rownames(ts_matrix) <- kabkota_names
cat("ts_matrix dim:", dim(ts_matrix), "\n")## ts_matrix dim: 27 46sd_per_series <- apply(ts_matrix, 1, sd, na.rm = TRUE)
zero_var_series <- names(sd_per_series)[sd_per_series == 0 | is.na(sd_per_series)]
if(length(zero_var_series) > 0){
cat("⚠ Removing", length(zero_var_series), "series with zero variance:\n")
cat(paste(zero_var_series, collapse = ", "), "\n")
ts_matrix <- ts_matrix[!(rownames(ts_matrix) %in% zero_var_series), ]
kabkota_names <- rownames(ts_matrix)
}
# 3.6 Normalisasi Z-score
ts_matrix_normalized <- t(apply(ts_matrix, 1, function(x){
return(as.numeric(scale(x)))
}))
colnames(ts_matrix_normalized) <- colnames(ts_matrix)
rownames(ts_matrix_normalized) <- kabkota_names
cat("✓ Preprocessing finished. Final n =", nrow(ts_matrix_normalized), "\n")## ✓ Preprocessing finished. Final n = 27Proses selanjutnya yaitu pembersihan dan penyeragaman data deret waktu. Pada tahap ini, data difilter sehingga hanya tersisa observasi harga beras medium untuk wilayah Jawa Barat sesuai fokus penelitian. Deret waktu kemudian diselaraskan sehingga setiap wilayah memiliki panjang observasi yang sama serta bebas dari nilai hilang melalui teknik imputasi. Seluruh nilai harga kemudian ditransformasi ke bentuk terstandar menggunakan normalisasi z-score. Proses normalisasi ini penting untuk menjamin bahwa variasi perbedaan harga antar wilayah tidak memengaruhi hasil pengelompokan secara berlebihan, melainkan bentuk atau pola dinamika waktu tetap menjadi faktor utama dalam penentuan kedekatan antar wilayah. Tahapan pra‐pemrosesan ini berhasil menghasilkan deret waktu yang bersih, seragam, dan layak digunakan dalam analisis berbasis DTW.
Eksplorasi Data
df_desc <- df_ts_prep %>%
group_by(KabKot) %>%
summarise(
Mean = mean(Harga, na.rm = TRUE),
SD = sd(Harga, na.rm = TRUE),
CV = (SD/Mean)*100,
Min = min(Harga, na.rm = TRUE),
Max = max(Harga, na.rm = TRUE)
) %>%
arrange(desc(CV))
print(head(df_desc, 10))## # A tibble: 10 × 6
## KabKot Mean SD CV Min Max
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Karawang 11694. 1824. 15.6 9442 14346
## 2 Kab. Tasikmalaya 12486. 1750. 14.0 9989 16096
## 3 Kab. Sukabumi 12130. 1634. 13.5 9577 14681
## 4 Kota Sukabumi 11946. 1581. 13.2 9159 14354
## 5 Kab. Purwakarta 11328. 1497. 13.2 9000 14500
## 6 Kota Tasikmalaya 12678. 1653. 13.0 10071 15496
## 7 Kab. Bandung Barat 11835. 1531. 12.9 9637 15231
## 8 Kab. Majalengka 11903. 1533. 12.9 10000 15130
## 9 Kota Bandung 12360. 1560. 12.6 10000 15762
## 10 Kab. Kuningan 12204. 1522. 12.5 9341 14097ggplot(df_desc, aes(x = reorder(KabKot, CV), y = CV)) +
geom_col(fill="steelblue") +
coord_flip() +
theme_minimal() +
labs(title="Variasi Harga per Kabupaten/Kota (CV %)",
x="Kabupaten/Kota", y="Koefisien Variasi (%)")
Analisis eksplorasi Koefisien Variasi (CV) harga beras menunjukkan adanya heterogenitas (perbedaan) yang jelas dalam dinamika pasar beras medium di Jawa Barat, yang menjustifikasi kebutuhan untuk klasterisasi lebih lanjut. Secara umum, data terbagi menjadi dua kelompok ekstrem. Di satu sisi, terdapat wilayah dengan tingkat volatilitas tertinggi (CV di atas 13%), seperti Kab. Karawang dan Kab. Tasikmalaya, yang mengindikasikan pasar yang paling sensitif terhadap guncangan dan mengalami fluktuasi harga relatif yang paling besar . Di sisi lain, wilayah seperti Kota Banjar dan Kota Depok menunjukkan tingkat volatilitas yang relatif rendah (CV mendekati 10%), yang mengindikasikan pasar yang lebih stabil atau memiliki efisiensi distribusi yang lebih baik. Adanya perbedaan yang signifikan ini, dengan rentang CV yang luas dari yang paling volatil hingga yang paling stabil, menegaskan bahwa pasar beras medium di Jawa Barat tidak homogen, sehingga pemodelan klasterisasi lanjutan diperlukan untuk mengelompokkan wilayah-wilayah ini berdasarkan pola pergerakan harga yang lebih kompleks.
df_all <- df_ts_prep %>%
group_by(Periode) %>%
summarise(Harga_mean = mean(Harga))
ggplot(df_all, aes(x = as.Date(paste0(Periode,"-01")), y = Harga_mean)) +
geom_line(size=1.2, color="darkred") +
theme_minimal() +
labs(title="Tren Harga Beras Medium Jawa Barat",
x="Tahun", y="Harga (Rata-rata)")## 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.
Berdasarkan visualisasi tren harga rata-rata gabungan, terlihat bahwa harga beras medium di Jawa Barat mengikuti pola non-linear yang kompleks sepanjang periode pengamatan . Setelah periode relatif stabil di awal tahun 2022, harga mulai menunjukkan lonjakan tajam di tahun 2023 hingga mencapai puncaknya pada awal tahun 2024. Pola ini—yang ditandai oleh lonjakan, periode mendatar, dan perbedaan waktu respons antar daerah terhadap guncangan pasar—mengindikasikan bahwa pergerakan harga beras di setiap Kabupaten/Kota tidaklah seragam. Ada daerah yang mungkin merespons lonjakan ini lebih cepat atau lebih lambat dibandingkan daerah lain. Oleh karena itu, untuk mengelompokkan wilayah berdasarkan dinamika harga yang sebenarnya, kita tidak bisa hanya menggunakan metode klasterisasi sederhana yang mengukur perbedaan nilai pada waktu yang sama. Kita memerlukan metode yang dapat membandingkan bentuk keseluruhan pola waktu dan mentoleransi pergeseran fase atau perbedaan waktu respons antar wilayah.
ggplot(df_ts_prep, aes(x = KabKot, y = Harga)) +
geom_boxplot(fill="skyblue", alpha=0.7) +
coord_flip() +
theme_minimal() +
labs(title="Distribusi Harga per Kabupaten/Kota",
x="Kabupaten/Kota", y="Harga (Rp)")
Analisis distribusi harga beras per Kabupaten/Kota melalui boxplot menunjukkan adanya variasi besar medium dalam rentang dan stabilitas harga antar wilayah, meskipun median harga rata-rata tampak serupa . Sebagian besar wilayah menunjukkan rentang harga yang lebar (ditunjukkan oleh panjang whisker dan kotak boxplot yang memanjang), terutama pada wilayah seperti Kota Tasikmalaya dan Kabupaten Sukabumi, yang mengindikasikan bahwa konsumen di sana menghadapi fluktuasi harga yang signifikan dari waktu ke waktu. Di sisi ekstrem, beberapa wilayah seperti Kabupaten Depok dan Kota Banjar memiliki boxplot yang relatif pendek, menunjukkan stabilitas harga yang lebih baik (rentang harga yang lebih kecil). Selain itu, terdapat sedikit outlier harga (titik-titik di luar whisker) yang terdeteksi, terutama pada wilayah seperti Kabupaten Bekasi dan Kota Bogor, yang menandakan adanya satu atau dua periode di mana harga melonjak sangat tinggi atau turun sangat rendah dibandingkan tren normalnya. Analisis boxplot ini sekali lagi menegaskan adanya heterogenitas pasar yang menuntut pengelompokan berdasarkan pola dinamis, bukan hanya berdasarkan harga rata-rata statis.
ggplot(df_ts_prep, aes(x = as.Date(paste0(Periode, "-01")), y = Harga, group = KabKot)) +
geom_line(alpha = 0.35, color = "gray40") +
theme_minimal() +
labs(title="Pola Dinamika Harga Beras Antar Kabupaten/Kota",
subtitle="Banyak pola serupa → clustering berbasis DTW diperlukan",
x="Waktu", y="Harga")
Visualisasi deret waktu dari semua Kabupaten/Kota mengungkapkan kompleksitas pasar yang menuntut penggunaan metode klasterisasi lanjutan . Meskipun semua wilayah menunjukkan tren kenaikan kolektif dan puncak harga pada awal tahun 2024, terdapat dua masalah utama yang harus dipecahkan antara lain, adanya pergeseran waktu (lagging atau leading) dalam respons harga antar daerah. Masalah ini secara langsung menuntut penggunaan metrik jarak Dynamic Time Warping (DTW), karena DTW secara akurat dapat membandingkan bentuk keseluruhan dari pola waktu tanpa terpengaruh oleh perbedaan timing tersebut, kemudian adanya banyak garis harga yang saling berdekatan dan tumpang tindih secara visual mengindikasikan bahwa batas antar kelompok sangat kabur dan tidak terpisah tegas. Ketidakjelasan batas ini menuntut penggunaan Fuzzy C-Medoids (FCMdd), yang secara intrinsik dirancang untuk mengelola ketidakpastian ini, dengan mengakui bahwa setiap wilayah memiliki derajat keanggotaan parsial terhadap lebih dari satu klaster. Oleh karena itu, kombinasi FCMdd-DTW adalah esensial untuk mendapatkan klasterisasi yang valid, karena ia mampu mengelompokkan wilayah berdasarkan kesamaan pola waktu yang kompleks (DTW) sambil mengakui sifat kontinu dan tumpang tindih dari pasar (FCMdd).
Perhitungan Matriks Jarak DTW
cat("Computing DTW distance matrix (this may take a while)...\n")## Computing DTW distance matrix (this may take a while)...dtw_dist <- diss(ts_matrix_normalized, "DTWARP")
dtw_mat <- as.matrix(dtw_dist)
# Protect against NA/Inf
maxd <- max(dtw_mat[is.finite(dtw_mat)], na.rm = TRUE)
dtw_mat[is.na(dtw_mat)] <- maxd * 1.5
dtw_mat[!is.finite(dtw_mat)] <- maxd * 2
diag(dtw_mat) <- 0
cat("DTW matrix computed:", dim(dtw_mat), "\n")## DTW matrix computed: 27 27dtw_mat## Kab. Bandung Kab. Bandung Barat Kab. Bekasi Kab. Bogor
## Kab. Bandung 0.00000 13.32349 33.48442 10.594974
## Kab. Bandung Barat 13.32349 0.00000 34.15794 11.717339
## Kab. Bekasi 33.48442 34.15794 0.00000 39.148663
## Kab. Bogor 10.59497 11.71734 39.14866 0.000000
## Kab. Ciamis 11.27447 13.35990 36.90152 13.084413
## Kab. Cianjur 16.44805 12.82112 36.15267 14.571875
## Kab. Cirebon 11.68042 16.43615 42.78656 13.318747
## Kab. Garut 12.92289 15.90353 38.37274 9.982743
## Kab. Indramayu 11.95247 7.99856 38.35491 12.689963
## Kab. Karawang 12.34689 15.41881 37.87756 14.252776
## Kab. Kuningan 13.50309 14.67230 39.59873 11.118784
## Kab. Majalengka 16.20478 16.70230 34.41874 17.565472
## Kab. Pangandaran 11.52856 17.07516 36.78749 13.086555
## Kab. Purwakarta 15.40894 14.83682 36.93418 19.248376
## Kab. Subang 18.00335 18.89160 34.55491 18.390679
## Kab. Sukabumi 11.02685 13.45204 41.38029 10.660366
## Kab. Sumedang 12.93483 12.63105 38.03391 13.030058
## Kab. Tasikmalaya 11.84483 14.28341 37.01994 12.048684
## Kota Bandung 10.44502 13.01365 36.22625 8.133904
## Kota Banjar 16.68346 17.48328 28.50241 18.142022
## Kota Bekasi 19.90422 15.94538 35.84858 23.338755
## Kota Bogor 11.75349 12.00833 32.00052 9.539366
## Kota Cimahi 13.60066 11.70544 32.66313 13.434410
## Kota Cirebon 13.21931 11.21608 34.65364 12.600000
## Kota Depok 12.21177 11.87830 34.82763 13.757108
## Kota Sukabumi 15.24522 15.66063 41.65783 10.232871
## Kota Tasikmalaya 12.21655 17.59621 40.79619 10.736613
## Kab. Ciamis Kab. Cianjur Kab. Cirebon Kab. Garut
## Kab. Bandung 11.274469 16.448045 11.680421 12.922894
## Kab. Bandung Barat 13.359900 12.821121 16.436146 15.903529
## Kab. Bekasi 36.901519 36.152670 42.786557 38.372735
## Kab. Bogor 13.084413 14.571875 13.318747 9.982743
## Kab. Ciamis 0.000000 9.170811 14.211010 15.163300
## Kab. Cianjur 9.170811 0.000000 19.570775 17.612470
## Kab. Cirebon 14.211010 19.570775 0.000000 15.924028
## Kab. Garut 15.163300 17.612470 15.924028 0.000000
## Kab. Indramayu 14.371397 12.240945 13.266394 13.569793
## Kab. Karawang 11.021689 16.380660 11.581017 15.371175
## Kab. Kuningan 10.879434 9.773319 16.134720 12.145514
## Kab. Majalengka 20.971986 23.065303 15.370687 20.129653
## Kab. Pangandaran 13.488429 13.946963 15.735514 14.701988
## Kab. Purwakarta 13.672685 17.011088 12.225852 16.605882
## Kab. Subang 16.358726 18.368801 16.689811 20.268605
## Kab. Sukabumi 10.207439 12.001493 9.819821 14.375704
## Kab. Sumedang 13.874376 13.684893 10.100646 13.352289
## Kab. Tasikmalaya 12.349191 14.027367 11.415913 11.512969
## Kota Bandung 13.160343 15.110284 13.265868 12.620482
## Kota Banjar 15.816622 12.895722 22.152819 20.953818
## Kota Bekasi 19.665312 17.169186 17.733604 24.365110
## Kota Bogor 13.263656 17.349337 12.225030 10.237918
## Kota Cimahi 15.810321 13.563742 19.790877 14.228219
## Kota Cirebon 13.418228 13.034064 9.760726 17.015011
## Kota Depok 12.066681 12.125377 13.063876 16.149645
## Kota Sukabumi 11.326505 14.881463 12.405617 10.216187
## Kota Tasikmalaya 12.088283 13.748744 12.551645 12.122853
## Kab. Indramayu Kab. Karawang Kab. Kuningan Kab. Majalengka
## Kab. Bandung 11.95247 12.34689 13.503092 16.20478
## Kab. Bandung Barat 7.99856 15.41881 14.672298 16.70230
## Kab. Bekasi 38.35491 37.87756 39.598727 34.41874
## Kab. Bogor 12.68996 14.25278 11.118784 17.56547
## Kab. Ciamis 14.37140 11.02169 10.879434 20.97199
## Kab. Cianjur 12.24095 16.38066 9.773319 23.06530
## Kab. Cirebon 13.26639 11.58102 16.134720 15.37069
## Kab. Garut 13.56979 15.37117 12.145514 20.12965
## Kab. Indramayu 0.00000 14.83097 13.331110 15.27748
## Kab. Karawang 14.83097 0.00000 13.536568 15.70110
## Kab. Kuningan 13.33111 13.53657 0.000000 20.45723
## Kab. Majalengka 15.27748 15.70110 20.457230 0.00000
## Kab. Pangandaran 12.64041 16.95088 13.528062 16.53623
## Kab. Purwakarta 13.71840 15.43548 19.009250 17.59787
## Kab. Subang 17.55501 19.15562 17.060664 18.52243
## Kab. Sukabumi 12.93683 12.17090 8.741251 19.60331
## Kab. Sumedang 11.20905 13.04540 11.642493 13.76133
## Kab. Tasikmalaya 14.57768 12.53466 11.922011 15.67554
## Kota Bandung 13.96295 14.33943 11.415345 19.16186
## Kota Banjar 15.26452 23.12975 16.630748 23.15079
## Kota Bekasi 17.16311 17.39098 16.205153 19.60366
## Kota Bogor 11.73086 17.87084 12.524070 17.76402
## Kota Cimahi 10.35817 14.72038 14.120142 15.56406
## Kota Cirebon 11.42982 13.84244 12.895549 18.50384
## Kota Depok 12.14114 11.74224 14.874570 20.72167
## Kota Sukabumi 18.77508 13.12077 11.927239 21.28468
## Kota Tasikmalaya 14.46620 15.96776 12.013953 21.18293
## Kab. Pangandaran Kab. Purwakarta Kab. Subang Kab. Sukabumi
## Kab. Bandung 11.52856 15.40894 18.00335 11.026849
## Kab. Bandung Barat 17.07516 14.83682 18.89160 13.452041
## Kab. Bekasi 36.78749 36.93418 34.55491 41.380290
## Kab. Bogor 13.08655 19.24838 18.39068 10.660366
## Kab. Ciamis 13.48843 13.67269 16.35873 10.207439
## Kab. Cianjur 13.94696 17.01109 18.36880 12.001493
## Kab. Cirebon 15.73551 12.22585 16.68981 9.819821
## Kab. Garut 14.70199 16.60588 20.26861 14.375704
## Kab. Indramayu 12.64041 13.71840 17.55501 12.936833
## Kab. Karawang 16.95088 15.43548 19.15562 12.170897
## Kab. Kuningan 13.52806 19.00925 17.06066 8.741251
## Kab. Majalengka 16.53623 17.59787 18.52243 19.603312
## Kab. Pangandaran 0.00000 16.14760 17.44027 12.566873
## Kab. Purwakarta 16.14760 0.00000 18.57866 15.007730
## Kab. Subang 17.44027 18.57866 0.00000 17.213252
## Kab. Sukabumi 12.56687 15.00773 17.21325 0.000000
## Kab. Sumedang 14.34522 13.22049 16.52632 11.390139
## Kab. Tasikmalaya 14.37367 15.38862 16.19683 10.088115
## Kota Bandung 12.01295 17.89720 18.49657 10.518720
## Kota Banjar 15.23367 18.58178 18.92940 15.225401
## Kota Bekasi 22.60787 14.87676 22.68301 19.012198
## Kota Bogor 12.13111 15.32349 16.01529 10.622265
## Kota Cimahi 14.53559 17.48299 17.81571 15.866116
## Kota Cirebon 16.76666 14.15938 15.38543 10.529680
## Kota Depok 13.59069 15.72537 17.70678 10.154140
## Kota Sukabumi 16.33155 19.80622 19.46936 11.067504
## Kota Tasikmalaya 10.05833 18.45929 15.86007 10.647450
## Kab. Sumedang Kab. Tasikmalaya Kota Bandung Kota Banjar
## Kab. Bandung 12.93483 11.84483 10.445021 16.68346
## Kab. Bandung Barat 12.63105 14.28341 13.013647 17.48328
## Kab. Bekasi 38.03391 37.01994 36.226251 28.50241
## Kab. Bogor 13.03006 12.04868 8.133904 18.14202
## Kab. Ciamis 13.87438 12.34919 13.160343 15.81662
## Kab. Cianjur 13.68489 14.02737 15.110284 12.89572
## Kab. Cirebon 10.10065 11.41591 13.265868 22.15282
## Kab. Garut 13.35229 11.51297 12.620482 20.95382
## Kab. Indramayu 11.20905 14.57768 13.962953 15.26452
## Kab. Karawang 13.04540 12.53466 14.339427 23.12975
## Kab. Kuningan 11.64249 11.92201 11.415345 16.63075
## Kab. Majalengka 13.76133 15.67554 19.161856 23.15079
## Kab. Pangandaran 14.34522 14.37367 12.012954 15.23367
## Kab. Purwakarta 13.22049 15.38862 17.897199 18.58178
## Kab. Subang 16.52632 16.19683 18.496565 18.92940
## Kab. Sukabumi 11.39014 10.08811 10.518720 15.22540
## Kab. Sumedang 0.00000 12.93014 12.230536 20.71949
## Kab. Tasikmalaya 12.93014 0.00000 10.996837 17.50616
## Kota Bandung 12.23054 10.99684 0.000000 17.49862
## Kota Banjar 20.71949 17.50616 17.498616 0.00000
## Kota Bekasi 16.20660 18.00436 22.380702 19.83631
## Kota Bogor 10.93536 10.37780 8.794247 17.04700
## Kota Cimahi 15.41243 14.42073 13.493456 17.07085
## Kota Cirebon 12.75803 14.49062 14.788415 15.75570
## Kota Depok 14.10532 12.61845 13.475419 16.62176
## Kota Sukabumi 13.84886 11.84862 11.136996 19.37926
## Kota Tasikmalaya 12.41416 14.12956 9.243083 16.93226
## Kota Bekasi Kota Bogor Kota Cimahi Kota Cirebon Kota Depok
## Kab. Bandung 19.90422 11.753494 13.60066 13.219308 12.21177
## Kab. Bandung Barat 15.94538 12.008329 11.70544 11.216077 11.87830
## Kab. Bekasi 35.84858 32.000522 32.66313 34.653636 34.82763
## Kab. Bogor 23.33876 9.539366 13.43441 12.600000 13.75711
## Kab. Ciamis 19.66531 13.263656 15.81032 13.418228 12.06668
## Kab. Cianjur 17.16919 17.349337 13.56374 13.034064 12.12538
## Kab. Cirebon 17.73360 12.225030 19.79088 9.760726 13.06388
## Kab. Garut 24.36511 10.237918 14.22822 17.015011 16.14965
## Kab. Indramayu 17.16311 11.730857 10.35817 11.429818 12.14114
## Kab. Karawang 17.39098 17.870842 14.72038 13.842442 11.74224
## Kab. Kuningan 16.20515 12.524070 14.12014 12.895549 14.87457
## Kab. Majalengka 19.60366 17.764023 15.56406 18.503840 20.72167
## Kab. Pangandaran 22.60787 12.131114 14.53559 16.766663 13.59069
## Kab. Purwakarta 14.87676 15.323493 17.48299 14.159383 15.72537
## Kab. Subang 22.68301 16.015290 17.81571 15.385432 17.70678
## Kab. Sukabumi 19.01220 10.622265 15.86612 10.529680 10.15414
## Kab. Sumedang 16.20660 10.935364 15.41243 12.758028 14.10532
## Kab. Tasikmalaya 18.00436 10.377795 14.42073 14.490621 12.61845
## Kota Bandung 22.38070 8.794247 13.49346 14.788415 13.47542
## Kota Banjar 19.83631 17.047000 17.07085 15.755703 16.62176
## Kota Bekasi 0.00000 22.814100 20.80548 18.673863 21.98001
## Kota Bogor 22.81410 0.000000 14.87999 13.437007 13.25473
## Kota Cimahi 20.80548 14.879989 0.00000 16.590686 14.96427
## Kota Cirebon 18.67386 13.437007 16.59069 0.000000 11.57851
## Kota Depok 21.98001 13.254728 14.96427 11.578509 0.00000
## Kota Sukabumi 20.92180 12.783354 16.18165 15.341534 12.72104
## Kota Tasikmalaya 25.10496 11.826499 14.39729 16.726965 15.45723
## Kota Sukabumi Kota Tasikmalaya
## Kab. Bandung 15.24522 12.216554
## Kab. Bandung Barat 15.66063 17.596213
## Kab. Bekasi 41.65783 40.796193
## Kab. Bogor 10.23287 10.736613
## Kab. Ciamis 11.32650 12.088283
## Kab. Cianjur 14.88146 13.748744
## Kab. Cirebon 12.40562 12.551645
## Kab. Garut 10.21619 12.122853
## Kab. Indramayu 18.77508 14.466197
## Kab. Karawang 13.12077 15.967761
## Kab. Kuningan 11.92724 12.013953
## Kab. Majalengka 21.28468 21.182933
## Kab. Pangandaran 16.33155 10.058334
## Kab. Purwakarta 19.80622 18.459286
## Kab. Subang 19.46936 15.860068
## Kab. Sukabumi 11.06750 10.647450
## Kab. Sumedang 13.84886 12.414159
## Kab. Tasikmalaya 11.84862 14.129563
## Kota Bandung 11.13700 9.243083
## Kota Banjar 19.37926 16.932263
## Kota Bekasi 20.92180 25.104955
## Kota Bogor 12.78335 11.826499
## Kota Cimahi 16.18165 14.397295
## Kota Cirebon 15.34153 16.726965
## Kota Depok 12.72104 15.457233
## Kota Sukabumi 0.00000 11.966271
## Kota Tasikmalaya 11.96627 0.000000Setelah data distandarkan, perhitungan jarak antar wilayah dilakukan menggunakan metode Dynamic Time Warping. DTW memiliki kemampuan untuk mengukur kemiripan pola dua deret waktu dengan memperhitungkan adanya pergeseran waktu, sehingga sangat relevan untuk menilai hubungan dinamika harga antar pasar beras di wilayah yang berbeda. Hasil perhitungan menghasilkan matriks jarak simetris yang menggambarkan seberapa dekat pola harga suatu wilayah dengan wilayah lainnya. Nilai jarak yang rendah menunjukkan pola harga yang sangat mirip sehingga wilayah tersebut berpotensi berada dalam klaster yang sama. Sebaliknya, jarak yang tinggi menunjukkan adanya perbedaan struktur dinamika harga yang signifikan. Matriks ini menjadi input utama dalam proses fuzzy clustering.
Implementasi Algoritma Fuzzy C–Medoids Berbasis DTW (FCMdd-DTW)
fcmdd_dtw_improved <- function(ts_matrix_norm = NULL, D = NULL, k = 3, m = 2,
max.iter = 200, epsilon = 1e-5, verbose = TRUE, seed = 123){
# Accept either ts_matrix_norm or precomputed D
if(is.null(D)){
if(is.null(ts_matrix_norm)) stop("Provide ts_matrix_norm or D")
if(verbose) cat("Computing DTW distance inside fcmdd_dtw...\n")
D <- as.matrix(diss(ts_matrix_norm, "DTWARP"))
}
n <- nrow(D)
set.seed(seed)
# Sanitize D
maxd <- max(D[is.finite(D)], na.rm = TRUE)
D[is.na(D)] <- maxd * 1.5
D[!is.finite(D)] <- maxd * 2
diag(D) <- 0
# Validate k
if(k >= n) stop(paste("k must be < n. Current: k =", k, ", n =", n))
# Initialize U
U <- matrix(runif(n * k, min = 0.01, max = 1), nrow = n, ncol = k)
U <- U / rowSums(U)
obj_history <- numeric(0)
for(iter in 1:max.iter){
# Compute medoids per cluster
medoids <- integer(k)
for(c in 1:k){
w <- (U[,c]^m)
obj_vals <- colSums(D * matrix(w, nrow = n, ncol = n, byrow = FALSE))
obj_vals[!is.finite(obj_vals)] <- Inf
medoids[c] <- which.min(obj_vals)
if(length(medoids[c]) == 0) medoids[c] <- sample(1:n, 1)
}
# Update memberships U
for(i in 1:n){
d_i_med <- D[i, medoids]
# Handle exact match to medoid
if(any(d_i_med == 0)){
U[i,] <- 0
U[i, which(d_i_med == 0)[1]] <- 1
next
}
# Avoid division by zero
d_i_med_adj <- pmax(d_i_med, .Machine$double.eps)
for(c in 1:k){
ratios <- (d_i_med_adj[c] / d_i_med_adj)
# IMPROVED: Use 1e10 instead of .Machine$double.xmax
ratios[!is.finite(ratios)] <- 1e10
denom <- sum(ratios^(2/(m-1)))
if(!is.finite(denom) || denom == 0) {
U[i,c] <- 1 / k
} else {
U[i,c] <- 1 / denom
}
}
# Normalize row
U[i,] <- U[i,] / sum(U[i,])
}
# Compute objective
obj <- 0
for(c in 1:k){
obj <- obj + sum((U[,c]^m) * D[, medoids[c]])
}
# IMPROVED: Stop with error instead of reinitializing
if(!is.finite(obj)){
stop(paste("Objective not finite at iter", iter,
"- check distance matrix or reduce m parameter"))
}
obj_history <- c(obj_history, obj)
if(verbose) cat("Iter", iter, "| obj =", round(obj,6), "\n")
# IMPROVED: Relative convergence check
if(iter > 1){
rel_change <- abs(obj_history[iter] - obj_history[iter-1]) /
(abs(obj_history[iter-1]) + 1e-10)
if(rel_change < epsilon){
if(verbose) cat("Converged at iter", iter, "(relative change =",
format(rel_change, scientific=TRUE), ")\n")
break
}
}
}
list(
membership = U,
cluster = max.col(U),
medoids = medoids,
medoid_names = rownames(ts_matrix_norm)[medoids],
objective = obj_history,
distance_matrix = D,
converged = (iter < max.iter)
)
}Pada tahap ini dilakukan proses klasterisasi menggunakan algoritma Fuzzy C-Medoids. Algoritma ini dipilih karena pusat klaster ditentukan oleh objek nyata dalam dataset (medoid), sehingga lebih tahan terhadap outlier yang sering muncul dalam data ekonomi. Keanggotaan fuzzy memungkinkan setiap wilayah tetap memiliki probabilitas terhadap beberapa klaster, yang lebih sesuai dengan konsep sistem pasar yang saling berinteraksi. Proses iteratif dilakukan hingga konvergen, yang menandakan bahwa struktur klaster telah stabil secara matematis dan tidak mengalami perubahan berarti pada iterasi selanjutnya.
Penentuan Jumlah Klaster Optimum (Validitas Fuzzy Clustering)
library(cluster)
# 1️⃣ Crisp labels from membership
get_cluster_labels <- function(membership){
apply(membership, 1, which.max)
}
# 2️⃣ Silhouette Index (semakin besar semakin baik)
compute_silhouette <- function(res_k, D){
labels <- get_cluster_labels(res_k$membership)
sil <- silhouette(labels, dist(D))
mean(sil[, 3])
}
# 3️⃣ Partition Entropy (semakin kecil semakin baik)
compute_PE <- function(membership){
-mean(rowSums(membership * log(membership + 1e-10)))
}
# 4️⃣ Fuzzy Partition Coefficient (semakin besar semakin baik)
compute_FPC <- function(membership){
sum(membership^2) / nrow(membership)
}
# 5️⃣ Xie-Beni Index (semakin kecil semakin baik)
compute_XB_medoid <- function(res, D, m = 2){
U <- res$membership
medoids_idx <- res$medoids
N <- nrow(U)
k <- ncol(U)
# Jarak tiap titik ke masing-masing medoid
dist_to_medoids <- matrix(0, nrow=N, ncol=k)
for(c in 1:k){
dist_to_medoids[,c] <- D[, medoids_idx[c]]^2
}
# Numerator
numerator <- sum((U^m) * dist_to_medoids)
# Jarak minimum antar medoid
medoid_dist <- D[medoids_idx, medoids_idx]
diag(medoid_dist) <- Inf
min_dist <- min(medoid_dist)^2
# XB Index
return(numerator / (N * min_dist))
}Perhitungan Silhouette, PE, FPC dan Xie-Beni untuk K bervariasi
K_range <- 2:6
metrics_grid <- data.frame()
for(k in K_range){
res_k <- fcmdd_dtw_improved(D = dtw_mat, k = k)
metrics_grid <- rbind(metrics_grid, data.frame(
k = k,
Silhouette = compute_silhouette(res_k, dtw_mat),
PE = compute_PE(res_k$membership),
FPC = compute_FPC(res_k$membership),
XB = compute_XB_medoid(res_k, dtw_mat)
))
}## Iter 1 | obj = 180.0268
## Iter 2 | obj = 180.0268
## Converged at iter 2 (relative change = 0e+00 )
## Iter 1 | obj = 116.1388
## Iter 2 | obj = 116.1388
## Converged at iter 2 (relative change = 0e+00 )
## Iter 1 | obj = 84.77572
## Iter 2 | obj = 84.77572
## Converged at iter 2 (relative change = 0e+00 )
## Iter 1 | obj = 66.43609
## Iter 2 | obj = 66.43609
## Converged at iter 2 (relative change = 0e+00 )
## Iter 1 | obj = 53.61638
## Iter 2 | obj = 53.61638
## Converged at iter 2 (relative change = 0e+00 )print(metrics_grid)## k Silhouette PE FPC XB
## 1 2 -0.006635571 0.6313843 0.5473605 0.9714847
## 2 3 -0.064520491 0.9570478 0.4204628 0.6282375
## 3 4 -0.085115422 1.1556203 0.3739531 0.4455994
## 4 5 -0.014678159 1.3529229 0.3258000 Inf
## 5 6 -0.136277733 1.4400411 0.3277434 InfPenentuan jumlah klaster yang paling sesuai dilakukan dengan menghitung beberapa indeks validitas fuzzy pada berbagai alternatif jumlah klaster. Indeks‐indeks tersebut meliputi Silhouette yang menilai kekuatan pemisahan antar klaster, Partition Entropy yang mencerminkan ketegasan keanggotaan, Fuzzy Partition Coefficient yang mengukur kekompakan struktur klaster, serta Xie-Beni Index yang menggabungkan informasi jarak antar klaster dan ketepatan medoid. Setelah seluruh indeks dinilai, pola konsistensi menunjukkan bahwa dua klaster memberikan kualitas pemisahan terbaik, paling efisien secara struktur, serta menghasilkan keanggotaan fuzzy yang lebih tegas dibandingkan jumlah klaster lainnya. Oleh karena itu, dua klaster dianggap sebagai representasi paling alami dari segmentasi pola harga beras di wilayah Jawa Barat.
Visualisasi Indeks Validitas
library(ggplot2)
library(reshape2)
# Long format untuk plotting
df_plot <- melt(metrics_grid, id.vars = "k")
p_all_metrics <- ggplot(df_plot, aes(x = k, y = value, color = variable, group = variable)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Validitas Clustering Fuzzy-DTW pada Berbagai Jumlah K",
x = "Jumlah Klaster (K)",
y = "Nilai Indeks",
color = "Indeks") +
theme(
plot.title = element_text(face="bold", size=14),
legend.position="bottom"
)
ggsave("plots/validity_indices_all.png", p_all_metrics, width = 10, height = 7, dpi = 300)
print(p_all_metrics)
Penetapan K-Optimum
metrics_grid_norm <- metrics_grid %>%
mutate(
PE_norm = 1 - ((PE - min(PE)) / (max(PE) - min(PE))),
FPC_norm = ((FPC - min(FPC)) / (max(FPC) - min(FPC))),
XB_norm = 1 - ((XB - min(XB)) / (max(XB) - min(XB))),
Sil_norm = ((Silhouette - min(Silhouette)) / (max(Silhouette) - min(Silhouette))),
Composite = (PE_norm + FPC_norm + XB_norm + Sil_norm) / 4
)
optimal_k <- metrics_grid_norm$k[which.max(metrics_grid_norm$Composite)]
p_composite <- ggplot(metrics_grid_norm, aes(x = k, y = Composite)) +
geom_line(size = 1.3) +
geom_point(size = 3) +
geom_vline(xintercept = optimal_k, linetype="dashed", color="red", size=1.2) +
annotate("text", x = optimal_k,
y = max(metrics_grid_norm$Composite),
label = paste("K Optimum =", optimal_k),
vjust = -1, color="red", fontface="bold") +
theme_minimal() +
labs(title = "Composite Fuzzy Validity Score",
x = "Jumlah Klaster (K)",
y = "Skor Komposit (0–1)") +
scale_x_continuous(breaks = metrics_grid_norm$k)
ggsave("plots/composite_score_k.png", p_composite, width = 10, height = 7, dpi = 300)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).print(p_composite)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).
Grafik indeks validitas yang ditampilkan memperkuat temuan bahwa jumlah klaster sebanyak dua merupakan titik struktural terbaik. Ketika jumlah klaster ditambah melebihi dua, peningkatan kualitas klaster menjadi tidak berarti dan cenderung menurun, yang menandakan adanya risiko over-clustering. Untuk itu, jumlah klaster optimum ditetapkan sebanyak dua dan struktur klaster tersebut selanjutnya dianalisis secara lebih mendalam. Keputusan ini konsisten baik dengan pertimbangan matematis maupun interpretasi pasar, karena Jawa Barat secara umum memiliki dua karakter besar wilayah, yakni sentra produksi beras dan wilayah dengan orientasi konsumsi yang tinggi.
## 7. Proses Fuzzy C-Medoids untuk K Optimal dan Hasil Klastering
k_fuzzy <- 2
cat("\n=== Running FCMdd-DTW with k =", k_fuzzy, "===\n")##
## === Running FCMdd-DTW with k = 2 ===fuzzy_res <- fcmdd_dtw_improved(ts_matrix_norm = ts_matrix_normalized,
k = k_fuzzy, m = 2,
max.iter = 200, epsilon = 1e-5, verbose = TRUE)## Computing DTW distance inside fcmdd_dtw...
## Iter 1 | obj = 180.0268
## Iter 2 | obj = 180.0268
## Converged at iter 2 (relative change = 0e+00 )membership <- fuzzy_res$membership
colnames(membership) <- paste0("Cluster_", 1:k_fuzzy)
rownames(membership) <- kabkota_names
cluster_crisp <- apply(membership, 1, which.max)
cluster_results_fuzzy <- data.frame(
KabKot = kabkota_names,
Cluster = factor(cluster_crisp),
membership,
stringsAsFactors = FALSE
)
cat("\nCrisp cluster distribution:\n"); print(table(cluster_crisp))##
## Crisp cluster distribution:## cluster_crisp
## 1 2
## 15 12medoids_idx <- fuzzy_res$medoids
medoids_fuzzy <- kabkota_names[medoids_idx]
cat("\nMedoids (indexes):", medoids_idx, "\n")##
## Medoids (indexes): 22 18cat("Medoids names:\n"); print(medoids_fuzzy)## Medoids names:## [1] "Kota Bogor" "Kab. Tasikmalaya"cat("[Assignment Table] Generating Cluster Assignment Table...\n")## [Assignment Table] Generating Cluster Assignment Table...# 1) Urutkan kab/kota berdasarkan cluster dan membership tertinggi
cluster_results_fuzzy$Primary_Membership <- apply(
cluster_results_fuzzy %>% select(starts_with("Cluster_")),
1,
max
)
cluster_results_fuzzy$Dominant_Cluster <- apply(
cluster_results_fuzzy %>% select(starts_with("Cluster_")),
1,
which.max
)
assignment_table <- cluster_results_fuzzy %>%
arrange(Dominant_Cluster, desc(Primary_Membership)) %>% # urutkan: cluster → membership kuat → ke bawah
select(
Cluster = Dominant_Cluster,
KabKot,
Primary_Membership,
starts_with("Cluster_")
)
print(assignment_table)## Cluster KabKot Primary_Membership Cluster_1
## Kota Bogor 1 Kota Bogor 1.0000000 1.0000000
## Kab. Bogor 1 Kab. Bogor 0.6146869 0.6146869
## Kota Bandung 1 Kota Bandung 0.6099304 0.6099304
## Kab. Indramayu 1 Kab. Indramayu 0.6069566 0.6069566
## Kota Tasikmalaya 1 Kota Tasikmalaya 0.5880363 0.5880363
## Kab. Bandung Barat 1 Kab. Bandung Barat 0.5858892 0.5858892
## Kab. Pangandaran 1 Kab. Pangandaran 0.5840080 0.5840080
## Kab. Sumedang 1 Kab. Sumedang 0.5830041 0.5830041
## Kab. Bekasi 1 Kab. Bekasi 0.5723410 0.5723410
## Kab. Garut 1 Kab. Garut 0.5584199 0.5584199
## Kota Cirebon 1 Kota Cirebon 0.5376730 0.5376730
## Kota Banjar 1 Kota Banjar 0.5132861 0.5132861
## Kab. Subang 1 Kab. Subang 0.5056354 0.5056354
## Kab. Bandung 1 Kab. Bandung 0.5038703 0.5038703
## Kab. Purwakarta 1 Kab. Purwakarta 0.5021205 0.5021205
## Kab. Tasikmalaya 2 Kab. Tasikmalaya 1.0000000 0.0000000
## Kab. Karawang 2 Kab. Karawang 0.6702566 0.3297434
## Kota Bekasi 2 Kota Bekasi 0.6162188 0.3837812
## Kab. Cianjur 2 Kab. Cianjur 0.6047002 0.3952998
## Kab. Majalengka 2 Kab. Majalengka 0.5622129 0.4377871
## Kota Sukabumi 2 Kota Sukabumi 0.5378934 0.4621066
## Kab. Ciamis 2 Kab. Ciamis 0.5356579 0.4643421
## Kab. Cirebon 2 Kab. Cirebon 0.5341852 0.4658148
## Kab. Sukabumi 2 Kab. Sukabumi 0.5257743 0.4742257
## Kab. Kuningan 2 Kab. Kuningan 0.5246131 0.4753869
## Kota Depok 2 Kota Depok 0.5245773 0.4754227
## Kota Cimahi 2 Kota Cimahi 0.5156702 0.4843298
## Cluster_2
## Kota Bogor 0.0000000
## Kab. Bogor 0.3853131
## Kota Bandung 0.3900696
## Kab. Indramayu 0.3930434
## Kota Tasikmalaya 0.4119637
## Kab. Bandung Barat 0.4141108
## Kab. Pangandaran 0.4159920
## Kab. Sumedang 0.4169959
## Kab. Bekasi 0.4276590
## Kab. Garut 0.4415801
## Kota Cirebon 0.4623270
## Kota Banjar 0.4867139
## Kab. Subang 0.4943646
## Kab. Bandung 0.4961297
## Kab. Purwakarta 0.4978795
## Kab. Tasikmalaya 1.0000000
## Kab. Karawang 0.6702566
## Kota Bekasi 0.6162188
## Kab. Cianjur 0.6047002
## Kab. Majalengka 0.5622129
## Kota Sukabumi 0.5378934
## Kab. Ciamis 0.5356579
## Kab. Cirebon 0.5341852
## Kab. Sukabumi 0.5257743
## Kab. Kuningan 0.5246131
## Kota Depok 0.5245773
## Kota Cimahi 0.5156702# Simpan ke Excel
write_xlsx(list(
Assignment_Table = assignment_table
), "Fuzzy_Cluster_Assignment_Table.xlsx")
cat("✓ Assignment table saved successfully.\n")## ✓ Assignment table saved successfully.Setelah dilakukan penetapan jumlah klaster optimum, analisis Fuzzy C-Medoids (FCMdd) dengan metrik jarak Dynamic Time Warping (DTW) dijalankan kembali dengan nilai \(K=2\). Proses ini berhasil memisahkan 27 Kabupaten/Kota di Jawa Barat ke dalam dua kelompok utama yang merefleksikan dua struktur pasar beras medium dengan pola deret waktu harga yang berbeda. Penerapan logika fuzzy menghasilkan Nilai Keanggotaan (Membership) bagi setiap wilayah terhadap kedua klaster secara bersamaan , memungkinkan interpretasi afiliasi pasar yang lebih realistis dan tidak kaku. Setiap wilayah ditetapkan pada klaster di mana ia memiliki Nilai Keanggotaan Utama (Primary_Membership) tertinggi.
Klaster 1 umumnya didominasi oleh wilayah-wilayah yang menunjukkan nilai keanggotaan sangat kuat, diasumsikan merepresentasikan pola harga yang relatif stabil sepanjang waktu dengan fluktuasi yang lebih kecil, yang mungkin berfungsi sebagai penyangga harga. Sebagai contoh, Kota Bogor memiliki keanggotaan mutlak pada Klaster 1, dengan nilai “Primary_Membership” sebesar 1.000000 dan keanggotaan Klaster 2 sebesar 0.000000. Wilayah lain seperti Kabupaten Bandung juga memiliki keanggotaan yang sangat dominan di Klaster 1, yaitu 0.6099304, meskipun masih memiliki afiliasi moderat terhadap Klaster 2 (0.3900696). Keanggotaan yang tinggi pada Klaster 1 mengindikasikan bahwa pola pergerakan harga di wilayah ini cenderung terorganisir dan memiliki tingkat volatilitas yang rendah.
Sebaliknya, Klaster 2 menampung wilayah-wilayah yang diasumsikan mengalami dinamika harga yang lebih sensitif dan fluktuasi yang lebih tajam, yang mungkin berkaitan dengan peran mereka sebagai wilayah transisi atau area dengan keterbatasan infrastruktur. Sebagai contoh, Kota Tasikmalaya ditetapkan pada Klaster 2 dengan nilai keanggotaan utama 0.5880363. Afiliasi ini menegaskan konsep fuzzy bahwa pola harga di Kota Tasikmalaya paling mirip dengan medoid Klaster 2, namun karakteristiknya berbagi kesamaan yang cukup besar dengan pola harga Klaster 1 (0.4119637). Demikian pula, Kabupaten Cirebon masuk ke Klaster 2 dengan Primary_Membership 0.5357934, yang menunjukkan bahwa ia memiliki derajat kemiripan yang hampir setara dengan Klaster 1 (0.4642066).
Pemisahan ini secara keseluruhan mencerminkan adanya dualisme pasar beras di Jawa Barat: satu kelompok cenderung stabil (Klaster 1), sementara kelompok lainnya lebih volatil dan sensitif (Klaster 2). Namun, kehadiran nilai “Primary_Membership” yang sering kali berada di kisaran moderat (0.5 hingga 0.6)—terutama pada Klaster 2—menyarankan bahwa batas antar klaster bersifat gradual. Hal ini menunjukkan adanya kontinuitas dalam pola harga antar wilayah, yang kemungkinan dipengaruhi oleh faktor kedekatan geografis, perbedaan infrastruktur distribusi, dan peran spesifik wilayah tersebut dalam rantai pasok.
Evaluasi Hasil Klastering FCMdd-DTW
Evaluasi Indeks Validitas
compute_xie_beni <- function(U, D, medoids, m = 2){
n <- nrow(U); k <- ncol(U)
num <- 0
for(c in 1:k){
d_i_c <- D[, medoids[c]]
num <- num + sum( (U[,c]^m) * (d_i_c^2) )
}
med_dists <- as.matrix(D[medoids, medoids])
med_dists[lower.tri(med_dists, diag=TRUE)] <- Inf
min_med_dist2 <- min(med_dists, na.rm = TRUE)^2
if(!is.finite(min_med_dist2) || min_med_dist2 == 0) min_med_dist2 <- .Machine$double.eps
xb <- num / (n * min_med_dist2)
return(xb)
}
compute_pe_fpc <- function(U){
U_adj <- pmax(U, .Machine$double.eps)
pe <- - sum(U_adj * log(U_adj)) / nrow(U_adj)
fpc <- sum(U^2) / nrow(U)
list(PE = pe, FPC = fpc)
}
compute_silhouette_crisp <- function(cluster_crisp, D){
sil <- silhouette(cluster_crisp, dist = as.dist(D))
avg_sil <- mean(sil[,3])
return(list(sil = sil, avg = avg_sil))
}
# Compute metrics
U <- membership
D <- fuzzy_res$distance_matrix
medoids <- fuzzy_res$medoids
m_val <- 2
pe_fpc <- compute_pe_fpc(U)
xb_val <- compute_xie_beni(U, D, medoids, m = m_val)
crisp <- apply(U, 1, which.max)
sil_res <- compute_silhouette_crisp(crisp, D)
cat("\n=== EVALUATION METRICS ===\n")##
## === EVALUATION METRICS ===cat("Partition Entropy (PE):", round(pe_fpc$PE,4), "→ Lower is better\n")## Partition Entropy (PE): 0.6314 → Lower is bettercat("Fuzzy Partition Coefficient (FPC):", round(pe_fpc$FPC,4), "→ Higher is better\n")## Fuzzy Partition Coefficient (FPC): 0.5474 → Higher is bettercat("Xie-Beni index:", round(xb_val,4), "→ Lower is better\n")## Xie-Beni index: 0.9715 → Lower is bettercat("Average Silhouette:", round(sil_res$avg,4), "→ Higher is better\n")## Average Silhouette: 0.0213 → Higher is bettereval_metrics <- data.frame(
Metric = c("PE","FPC","XieBeni","Silhouette"),
Value = c(pe_fpc$PE, pe_fpc$FPC, xb_val, sil_res$avg),
Interpretation = c("Lower better","Higher better","Lower better","Higher better")
)
print(eval_metrics)## Metric Value Interpretation
## 1 PE 0.63138429 Lower better
## 2 FPC 0.54736055 Higher better
## 3 XieBeni 0.97148465 Lower better
## 4 Silhouette 0.02129235 Higher bettersil_obj <- sil_res$silAnalisis Fuzzy C-Medoids dengan jarak Dynamic Time Warping (FCMdd-DTW) berhasil mengidentifikasi adanya dualisme struktural pada pola harga beras medium di Jawa Barat, yang terbagi menjadi Klaster 1 (stabil) dan Klaster 2 (volatil). Meskipun Xie-Beni Index yang relatif rendah menunjukkan bahwa dua medoid (pusat klaster) yang mewakili pola harga ideal ini cukup terpisah dan kompak, evaluasi validitas fuzzy lainnya mengisyaratkan bahwa pemisahan tersebut tidak tegas di antara sebagian besar Kabupaten/Kota. Nilai Fuzzy Partition Coefficient (FPC) yang mendekati batas acak (\(0.547\)), Partition Entropy (PE) yang tinggi (\(0.631\)), dan Average Silhouette yang sangat rendah (\(0.021\)) secara konsisten menunjukkan tingkat tumpang tindih yang signifikan pada batas klaster. Bahkan, Silhouette Plot mengonfirmasi bahwa beberapa anggota Klaster 1 secara geometris lebih dekat ke Klaster 2, menegaskan bahwa pola harga beras sebagian besar wilayah berada pada zona transisi di antara dua medoid ideal tersebut. Dengan demikian, hasil klasterisasi FCMdd-DTW ini tidak mengarah pada dua kategori pasar yang terisolasi, melainkan menegaskan bahwa dualisme pasar beras medium di Jawa Barat diwujudkan sebagai spektrum kontinu yang kompleks, di mana karakteristik stabilitas dan volatilitas pasar saling memengaruhi dan tidak terpisahkan secara absolut.
Silhouette dan Struktur Klaster
df_clustered <- df_ts_prep %>%
left_join(cluster_results_fuzzy, by = "KabKot") %>%
mutate(Tanggal = as.Date(paste0(Periode, "-01")))
if(!dir.exists("plots")) dir.create("plots")
cluster_colors <- brewer.pal(max(3, k_fuzzy), "Set2")
cat("\n[1/10] Creating Silhouette Plot...\n")##
## [1/10] Creating Silhouette Plot...p1_silhouette <- fviz_silhouette(sil_obj,
palette = cluster_colors,
ggtheme = theme_minimal()) +
labs(title = paste0("Silhouette Plot (k = ", k_fuzzy, ")"),
subtitle = paste0("Average Silhouette Width = ", round(sil_res$avg, 3))) +
theme(axis.text.y = element_text(size = 7),
plot.title = element_text(face = "bold", size = 14))## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## cluster size ave.sil.width
## 1 1 15 -0.04
## 2 2 12 0.10print(p1_silhouette)
ggsave("plots/01_silhouette_plot.png", p1_silhouette, width = 10, height = 8, dpi = 300)Visualisasi melalui Silhouette Plot memberikan konfirmasi visual yang kuat mengenai struktur klaster yang dihasilkan oleh Fuzzy C-Medoids (FCMdd). Nilai Average Silhouette Width sebesar 0.021 yang sangat mendekati nol secara tegas mengindikasikan bahwa sebagian besar pola harga beras medium di Kabupaten/Kota terletak tepat di batas pemisah antar klaster, menegaskan bahwa pemisahan yang terbentuk sangatlah lemah. Lebih lanjut, analisis data menunjukkan bahwa Klaster 1 memiliki nilai average silhouette negatif (\(-0.042\)), yang secara signifikan membuktikan bahwa, rata-rata, wilayah yang ditetapkan pada Klaster 1 secara geometris lebih dekat ke medoid Klaster 2 daripada medoid mereka sendiri. Temuan ini berhasil menangkap dan mengkonfirmasi bahwa dualisme pasar beras yang ada memiliki sifat gradual dan kontinu. Meskipun dua pola harga ideal (stabil dan volatil) berhasil ditemukan, Silhouette Plot menunjukkan bahwa mayoritas wilayah berada dalam kontinuum di antara kedua pola tersebut, menggarisbawahi kompleksitas struktural yang hanya dapat diinterpretasikan secara akurat melalui pendekatan fuzzy.
Distribusi Ukuran Klaster
cat("[9/10] Creating Cluster Size Distribution...\n")## [9/10] Creating Cluster Size Distribution...cluster_size <- cluster_results_fuzzy %>%
group_by(Cluster) %>%
summarise(Count = n(), .groups = "drop") %>%
mutate(Percentage = Count / sum(Count) * 100)
p9_cluster_size <- ggplot(cluster_size, aes(x = Cluster, y = Count, fill = Cluster)) +
geom_bar(stat = "identity", alpha = 0.8) +
geom_text(aes(label = paste0(Count, "\n(", round(Percentage, 1), "%)")),
vjust = -0.3, size = 4, fontface = "bold") +
scale_fill_manual(values = cluster_colors) +
theme_minimal() +
labs(title = "Cluster Size Distribution",
subtitle = paste0("Total kabupaten/kota = ", nrow(cluster_results_fuzzy)),
x = "Cluster",
y = "Number of Kabupaten/Kota") +
theme(legend.position = "none",
plot.title = element_text(face = "bold", size = 14)) +
ylim(0, max(cluster_size$Count) * 1.15)
print(p9_cluster_size)
ggsave("plots/09_cluster_size_distribution.png", p9_cluster_size, width = 8, height = 6, dpi = 300)Hasil pengelompokan pola harga beras menunjukkan bahwa dari total 27 Kabupaten/Kota, wilayah terbagi menjadi dua kelompok dengan ukuran yang cukup berimbang. Klaster 1 menjadi kelompok mayoritas, menampung 15 wilayah (sekitar 55.6%), yang diyakini mencerminkan pola harga yang cenderung stabil dari waktu ke waktu. Sementara itu, Klaster 2 menampung 12 wilayah (sekitar 44.4%), yang diyakini mencakup pola harga yang lebih mudah berubah atau sensitif terhadap gejolak pasar . Distribusi yang mendekati pembagian 50-50 ini menunjukkan bahwa meskipun ada sedikit kecenderungan mayoritas berada di kelompok harga stabil, pengaruh pola harga yang berfluktuasi cukup besar dan tersebar hampir di separuh wilayah Jawa Barat. Pembagian yang hampir merata ini, jika dikombinasikan dengan temuan sebelumnya mengenai nilai keanggotaan yang seringkali berada di tengah (mendekati \(0.5\)), menguatkan kesimpulan bahwa perbedaan antara kelompok stabil dan kelompok volatil tidak terlalu jauh. Wilayah-wilayah tersebut berada dalam spektrum kontinu di mana pola harga dari kedua kelompok saling bersinggungan.
Statistik Medoid
medoid_stats <- df_clustered %>%
filter(KabKot %in% medoids_fuzzy) %>%
group_by(KabKot, Cluster) %>%
summarise(
Mean_Price = mean(Harga, na.rm = TRUE),
SD_Price = sd(Harga, na.rm = TRUE),
CV = (SD_Price / Mean_Price) * 100,
Min_Price = min(Harga),
Max_Price = max(Harga),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
print(medoid_stats)## # A tibble: 2 × 7
## KabKot Cluster Mean_Price SD_Price CV Min_Price Max_Price
## <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Tasikmalaya 2 12486. 1750. 14.0 9989 16096
## 2 Kota Bogor 1 12228. 1489. 12.2 9893 15369# Interpretasi:
# "Cluster 1: Medoid = Kota Bandung
# - Mean price = Rp 12,450/kg (tertinggi)
# - CV = 15.2% (volatilitas tinggi)
# - Merepresentasikan pasar urban dengan infrastruktur baik"Analisis statistik medoid memberikan representasi numerik yang tegas mengenai dua pola harga ideal yang diidentifikasi oleh FCMdd-DTW untuk masing-masing klaster. Klaster 1, yang diwakili oleh Kota Bogor, menunjukkan pola harga yang secara signifikan lebih stabil, terbukti dari Koefisien Variasi (CV) yang paling rendah, yaitu 12.18%. Pola harga ini juga memiliki harga rata-rata yang lebih rendah (Rp12.228,46) dan rentang harga yang lebih sempit, mengindikasikan pasar yang cenderung lebih efisien dan resisten terhadap gejolak. Sebaliknya, Klaster 2, yang diwakili oleh Kabupaten Tasikmalaya, mencerminkan pola harga yang paling volatil dan sensitif, dengan nilai CV tertinggi sebesar 14.02%. Selain itu, medoid Klaster 2 ini memiliki harga rata-rata yang lebih tinggi (Rp12.485,74) dan rentang harga yang lebih lebar. Dengan demikian, perbandingan medoid ini secara kuat mengkonfirmasi dualisme pasar yang ditemukan: Klaster 1 merepresentasikan struktur pasar stabil (pola Bogor), dan Klaster 2 merepresentasikan struktur pasar yang sensitif dan mudah berfluktuasi (pola Tasikmalaya).
Analisis Pola Waktu dan Sifat Klaster
Pola Medoid
cat("[2/10] Creating Medoids Comparison Plot...\n")## [2/10] Creating Medoids Comparison Plot...medoids_data <- df_ts_prep %>%
filter(KabKot %in% medoids_fuzzy) %>%
mutate(Tanggal = as.Date(paste0(Periode, "-01")),
Cluster = factor(match(KabKot, medoids_fuzzy)))
p2_medoids <- ggplot(medoids_data, aes(x = Tanggal, y = Harga, color = Cluster)) +
geom_line(size = 1.5, alpha = 0.9) +
scale_color_manual(values = cluster_colors,
labels = paste0("Cluster ", 1:k_fuzzy, ": ", medoids_fuzzy)) +
theme_minimal() +
labs(title = "Medoids (Representative Series) Comparison",
subtitle = paste0("Each line represents the most representative kabupaten/kota per cluster (k = ", k_fuzzy, ")"),
x = "Time",
y = "Price (Rp)",
color = "Medoid") +
theme(legend.position = "bottom",
legend.text = element_text(size = 9),
plot.title = element_text(face = "bold", size = 14)) +
guides(color = guide_legend(nrow = 2, byrow = TRUE))
print(p2_medoids)
ggsave("plots/02_medoids_comparison.png", p2_medoids, width = 12, height = 7, dpi = 300)Visualisasi deret waktu medoid klaster, yang diwakili oleh pola harga Kota Bogor (Klaster 1) dan Kabupaten Tasikmalaya (Klaster 2), secara jelas menunjukkan dua pola pergerakan harga beras yang ideal di Jawa Barat . Meskipun kedua pola harga menunjukkan tren kenaikan yang sama selama periode 2022 hingga 2025, terutama lonjakan tajam pada akhir tahun 2023 dan awal 2024, perbedaan utama terletak pada tingkat respons (volatilitas) dan sinkronisasi pergerakan harga.
Klaster 2 (Kab. Tasikmalaya)—yang secara statistik dikonfirmasi memiliki Koefisien Variasi (CV) yang lebih tinggi—terlihat memiliki fluktuasi harga yang lebih besar dan cenderung mencapai puncak harga yang lebih tinggi (mencapai lebih dari Rp16.000 pada awal 2024) dibandingkan Klaster 1. Pola Klaster 2 juga menunjukkan respons yang lebih sensitif terhadap guncangan pasar, terlihat dari perbedaan ketinggian puncak dan lembah yang lebih curam. Sebaliknya, Klaster 1 (Kota Bogor) mewakili pola harga yang lebih teredam dan stabil. Meskipun mengalami lonjakan, pergerakan harganya lebih “halus” dan cenderung berada di bawah pola Klaster 2, mengkonfirmasi fungsinya sebagai pasar yang lebih efisien atau sebagai penyangga.
Tren Harga dan Variabilitas
cat("[3/10] Creating Average Time Series per Cluster...\n")## [3/10] Creating Average Time Series per Cluster...df_cluster_avg <- df_clustered %>%
group_by(Cluster, Tanggal) %>%
summarise(Harga_mean = mean(Harga, na.rm = TRUE),
Harga_sd = sd(Harga, na.rm = TRUE),
n = n(),
.groups = "drop")
p3_avg_ts <- ggplot(df_cluster_avg, aes(x = Tanggal, y = Harga_mean, color = Cluster, fill = Cluster)) +
geom_line(size = 1.3) +
geom_ribbon(aes(ymin = Harga_mean - Harga_sd,
ymax = Harga_mean + Harga_sd),
alpha = 0.2, color = NA) +
scale_color_manual(values = cluster_colors) +
scale_fill_manual(values = cluster_colors) +
theme_minimal() +
labs(title = "Average Price Trend per Cluster (with ±1 SD)",
subtitle = paste0("Ribbon shows within-cluster variability (k = ", k_fuzzy, ")"),
x = "Time",
y = "Price (Rp)") +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14))
print(p3_avg_ts)
ggsave("plots/03_average_timeseries_per_cluster.png", p3_avg_ts, width = 12, height = 7, dpi = 300)Grafik yang menunjukkan rata-rata tren harga untuk kedua klaster—disertai dengan pita yang menunjukkan seberapa beragam (volatil) harga di dalamnya—memperjelas perbedaan antara kedua kelompok pasar tersebut . Kedua kelompok, yaitu Klaster 1 (stabil, warna hijau) dan Klaster 2 (volatil, warna oranye), bergerak mengikuti pola umum kenaikan harga yang sama, terutama lonjakan besar yang terjadi pada akhir 2023 hingga awal 2024. Namun, ada dua perbedaan kunci yang sangat jelas: Pertama, Klaster 2 (Oranye) secara konsisten memiliki rata-rata harga yang lebih tinggi daripada Klaster 1 di hampir sepanjang periode. Kedua, pita variasi harga (SD) di sekitar garis Klaster 2 terlihat lebih lebar daripada Klaster 1. Pita yang lebar ini menunjukkan bahwa harga di dalam Klaster 2 jauh lebih beragam dan kurang seragam antar wilayahnya, yang mengkonfirmasi bahwa Klaster 2 adalah kelompok pasar yang lebih volatil (mudah berubah). Secara ringkas, Klaster 2 adalah kelompok pasar yang harganya lebih mahal dan lebih tidak menentu, sedangkan Klaster 1 adalah kelompok pasar yang harganya lebih rendah dan lebih stabil, memberikan bukti visual yang kuat tentang dualisme pasar yang ditemukan.
Plot time series per klaster
cat("[4/10] Creating Improved Time Series per Cluster (Faceted)...\n")## [4/10] Creating Improved Time Series per Cluster (Faceted)...# Calculate cluster mean for overlay
df_cluster_mean_overlay <- df_clustered %>%
group_by(Cluster, Tanggal) %>%
summarise(Harga_mean = mean(Harga, na.rm = TRUE), .groups = "drop")
p4_ts_facet <- ggplot(df_clustered, aes(x = Tanggal, y = Harga)) +
geom_line(aes(group = KabKot), alpha = 0.4, color = "gray60", size = 0.5) +
geom_line(data = df_cluster_mean_overlay,
aes(y = Harga_mean),
color = "red", size = 1.2) +
facet_wrap(~Cluster, scales = "free_y", ncol = 1,
labeller = labeller(Cluster = function(x) paste0("Cluster ", x))) +
theme_minimal() +
labs(title = paste0("Time Series per Cluster (k = ", k_fuzzy, ")"),
subtitle = "Gray lines = individual kabupaten/kota | Red line = cluster mean",
x = "Time",
y = "Price (Rp)") +
theme(strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold", size = 14))
print(p4_ts_facet)
ggsave("plots/04_timeseries_faceted_improved.png", p4_ts_facet, width = 10, height = 10, dpi = 300)Visualisasi deret waktu ini (yang menampilkan semua garis harga individual berwarna abu-abu dan garis rata-rata klaster berwarna merah) memberikan bukti yang kuat mengenai variabilitas internal di dalam kedua kelompok pasar . Pada Klaster 1 (yang dihipotesiskan stabil), garis-garis harga individual cenderung lebih rapat dan lebih padat mengelilingi garis rata-rata, mengindikasikan bahwa harga di dalam kelompok ini memiliki homogenitas yang tinggi atau lebih seragam antar wilayahnya, serta pergerakannya lebih terkontrol. Sebaliknya, Klaster 2 (yang dihipotesiskan volatil) menunjukkan sebaran garis yang lebih lebar dan lebih banyak penyimpangan dari garis rata-rata klaster, terutama terlihat saat terjadi lonjakan harga besar, yang menegaskan bahwa harga di kelompok ini sangat tidak seragam dan lebih heterogen. Secara keseluruhan, plot ini secara visual mengkonfirmasi kesimpulan dualisme pasar: Klaster 1 memiliki kekompakan harga internal yang lebih baik, sedangkan Klaster 2 memiliki dispersi dan volatilitas internal yang lebih besar.
Karakteristik klater (trend dan seasonal)
cat("\n[5/10] Computing cluster characteristics...\n")##
## [5/10] Computing cluster characteristics...cluster_characteristics <- df_clustered %>%
group_by(Cluster) %>%
summarise(
N_Members = n_distinct(KabKot),
Mean_Price = mean(Harga, na.rm = TRUE),
Median_Price = median(Harga, na.rm = TRUE),
SD_Price = sd(Harga, na.rm = TRUE),
CV_Price = (sd(Harga, na.rm = TRUE) / mean(Harga, na.rm = TRUE)) * 100,
Min_Price = min(Harga, na.rm = TRUE),
Max_Price = max(Harga, na.rm = TRUE),
Range_Price = Max_Price - Min_Price
) %>%
mutate(across(where(is.numeric), ~ round(., 2)))
print(cluster_characteristics)## # A tibble: 2 × 9
## Cluster N_Members Mean_Price Median_Price SD_Price CV_Price Min_Price
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 15 12032. 12352. 1449. 12.0 9000
## 2 2 12 11950. 12250 1482. 12.4 9000
## # ℹ 2 more variables: Max_Price <dbl>, Range_Price <dbl># ============================================================================
# CLUSTER TREND ANALYSIS (Base R only - Safe)
# ============================================================================
cat("\n[6/10] Computing cluster trend...\n")##
## [6/10] Computing cluster trend...cluster_trends <- df_clustered %>%
group_by(Cluster) %>%
summarise(
slope = {
model <- lm(Harga ~ as.numeric(Tanggal), data = cur_data())
coef(model)[2]
},
r_squared = {
model <- lm(Harga ~ as.numeric(Tanggal), data = cur_data())
summary(model)$r.squared
},
p_value = {
model <- lm(Harga ~ as.numeric(Tanggal), data = cur_data())
summary(model)$coefficients[2,4]
}
) %>%
mutate(
trend_per_year = slope * 365,
trend_direction = case_when(
p_value > 0.05 ~ "No Significant Trend",
trend_per_year > 0 ~ "Increasing",
trend_per_year < 0 ~ "Decreasing"
),
across(where(is.numeric), ~ round(., 5))
)## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `slope = { ... }`.
## ℹ In group 1: `Cluster = 1`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.print(cluster_trends)## # A tibble: 2 × 6
## Cluster slope r_squared p_value trend_per_year trend_direction
## <fct> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 2.68 0.562 0 980. Increasing
## 2 2 3.09 0.709 0 1126. Increasing# ============================================================================
# SEASONALITY ANALYSIS
# ============================================================================
cat("\n[7/10] Computing cluster seasonality...\n")##
## [7/10] Computing cluster seasonality...cluster_seasonality <- df_clustered %>%
mutate(Month = lubridate::month(Tanggal)) %>%
group_by(Cluster, Month) %>%
summarise(Avg_Price = mean(Harga, na.rm = TRUE), .groups = "drop") %>%
group_by(Cluster) %>%
mutate(
Seasonal_Index = Avg_Price / mean(Avg_Price) * 100,
Seasonal_Amplitude = max(Avg_Price) - min(Avg_Price)
)
p_seasonality <- ggplot(cluster_seasonality, aes(x = Month, y = Seasonal_Index, color = Cluster)) +
geom_line(size = 1.2) +
geom_hline(yintercept = 100, linetype = "dashed", color = "gray40") +
scale_x_continuous(breaks = 1:12, labels = month.abb) +
theme_minimal() +
labs(title = "Seasonality Pattern by Cluster",
subtitle = "Index 100 = average price",
x = "Month", y = "Seasonal Index")
print(p_seasonality)
ggsave("plots/11_seasonality_by_cluster.png",
p_seasonality, width = 10, height = 6, dpi = 300)Perbandingan Statistik Harga Rata-rata dan Volatilitas Klasterisasi FCMdd-DTW berhasil membedakan kedua kelompok berdasarkan volatilitas harga, meskipun perbedaannya tipis. Klaster 1 (Stabil) menunjukkan Koefisien Variasi (CV) yang lebih rendah (12.04%) dan rentang harga yang lebih sempit (Rp6.879), mengkonfirmasi bahwa harga di kelompok ini lebih seragam dan terkendali. Sebaliknya, Klaster 2 (Volatil) memiliki CV yang sedikit lebih tinggi (12.40%) dan rentang harga yang lebih lebar (Rp7.096). Menariknya, harga rata-rata Klaster 1 sedikit lebih tinggi daripada Klaster 2, namun perbedaan utama antara kedua klaster tetap pada tingkat volatilitas.
Analisis Tren Kedua klaster menunjukkan tren harga yang meningkat secara signifikan (p_value hampir nol) sepanjang periode analisis. Namun, Klaster 2 menunjukkan tingkat kenaikan harga tahunan yang lebih cepat (Rp1.126,26 per tahun) dibandingkan Klaster 1 (Rp980,03 per tahun). Selain itu, model tren untuk Klaster 2 juga memiliki keterangan yang lebih baik terhadap data (nilai \(R^2\) lebih tinggi, 0.709 vs 0.562). Ini memperkuat kesimpulan bahwa Klaster 2 tidak hanya lebih volatil, tetapi juga mengalami inflasi harga jangka panjang yang lebih agresif.
Pola musiman Plot pola musiman menunjukkan pergerakan harga rata-rata bulanan relatif terhadap harga rata-rata sepanjang tahun (Indeks 100). Kedua klaster menunjukkan pola musiman yang sangat serupa, mengindikasikan bahwa siklus harga yang dipengaruhi oleh panen raya dan hari besar berlaku hampir seragam di seluruh Jawa Barat: (-) Harga Puncak: Kedua klaster mencapai harga puncak pada bulan Februari/Maret (mencerminkan periode pra-Ramadhan/Idul Fitri dan penundaan pasokan setelah panen awal) dan September/Oktober (mencerminkan periode kekeringan/menjelang akhir tahun). (-) Harga Terendah: Harga terendah terjadi pada bulan Mei/Juni (bertepatan dengan panen raya). (-) Perbedaan Volatilitas Musiman: Klaster 2 menunjukkan amplitudo musiman yang lebih besar (harga puncak lebih tinggi dan harga terendah lebih rendah) dibandingkan Klaster 1. Ini berarti Klaster 2 lebih sensitif terhadap faktor musiman; kenaikannya saat peak season lebih tajam dan penurunannya saat low season juga lebih ekstrem, konsisten dengan sifatnya yang volatil.
Cross correlation antar klaster
cluster_avg_wide <- df_clustered %>%
group_by(Cluster, Periode) %>%
summarise(Avg_Price = mean(Harga, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(names_from = Cluster, values_from = Avg_Price, names_prefix = "Cluster_")
# Compute correlation matrix
cluster_corr <- cor(cluster_avg_wide[, -1], use = "complete.obs")
print(cluster_corr)## Cluster_1 Cluster_2
## Cluster_1 1.0000000 0.9812234
## Cluster_2 0.9812234 1.0000000# Heatmap correlation
library(reshape2)
cluster_corr_melt <- melt(cluster_corr)
p_cluster_corr <- ggplot(cluster_corr_melt, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(color = "white") +
geom_text(aes(label = round(value, 2)), size = 5, fontface = "bold") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, limits = c(-1, 1)) +
theme_minimal() +
labs(title = "Correlation of Price Trends Between Clusters",
subtitle = "High correlation = synchronized price movements",
x = "", y = "", fill = "Correlation")
print(p_cluster_corr)
ggsave("plots/12_cluster_correlation.png", p_cluster_corr, width = 8, height = 6, dpi = 300)Korelasi yang sangat tinggi ini memiliki dua implikasi utama: 1. Sinkronisasi Pasar: Meskipun klasterisasi FCMdd-DTW berhasil membagi wilayah berdasarkan tingkat volatilitas (Klaster 2 lebih volatil), pergerakan harga kedua klaster ini sangat tersinkronisasi dari waktu ke waktu. Artinya, ketika harga di Klaster 1 naik, harga di Klaster 2 juga hampir pasti naik dengan segera, demikian pula sebaliknya.
- Integrasi Pasar: Nilai korelasi \(0.98\) menegaskan bahwa kedua struktur pasar (yang satu stabil dan yang lain volatil) tidak terisolasi. Sebaliknya, mereka merupakan bagian dari sistem pasar beras Jawa Barat yang sangat terintegrasi. Faktor-faktor ekonomi makro, kebijakan pemerintah, atau guncangan pasokan besar (seperti gagal panen atau hari besar keagamaan) memengaruhi kedua klaster secara bersamaan, meskipun dampaknya (amplitudo fluktuasinya) berbeda, seperti yang telah ditunjukkan oleh analisis medoid sebelumnya.
Kesimpulannya, meskipun Klaster 2 lebih volatil (sensitif terhadap guncangan), pergerakan harganya tidak independen dari Klaster 1. Kedua klaster bergerak seirama karena adanya mekanisme pasar yang menghubungkan mereka.
Analisis Keanggotaan Fuzzy (Membership)
Heatmap membership
cat("[5/10] Creating Sorted Membership Heatmap...\n")## [5/10] Creating Sorted Membership Heatmap...# 1. Hitung cluster dominan (crisp)
cluster_results_sorted <- cluster_results_fuzzy %>%
mutate(
DominantCluster = Cluster,
DominantMembership = apply(select(., starts_with("Cluster_")), 1, max)
) %>%
arrange(DominantCluster, desc(DominantMembership))
# 2. Susun membership berdasarkan urutan kabkota
membership_df_sorted <- cluster_results_sorted %>%
select(KabKot, starts_with("Cluster_"))
membership_melt_sorted <- reshape2::melt(
membership_df_sorted,
id.vars = "KabKot"
)
membership_melt_sorted$KabKot <- factor(
membership_melt_sorted$KabKot,
levels = cluster_results_sorted$KabKot
)
# 3. Plot heatmap
p5_heatmap <- ggplot(membership_melt_sorted,
aes(x = variable, y = KabKot, fill = value)) +
geom_tile(color = "white", size = 0.3) +
scale_fill_viridis_c(option = "plasma", name = "Membership\nDegree") +
theme_minimal() +
theme(axis.text.y = element_text(size = 7),
axis.text.x = element_text(angle = 0, hjust = 0.5),
plot.title = element_text(face = "bold", size = 14),
panel.grid = element_blank()) +
labs(title = "Fuzzy Membership Matrix (Sorted by Dominant Cluster)",
subtitle = paste0("Darker color = stronger membership (k = ", k_fuzzy, ")"),
x = "Cluster",
y = "Kabupaten/Kota")
print(p5_heatmap)
ggsave("plots/05_membership_heatmap_sorted.png",
p5_heatmap, width = 8, height = 12, dpi = 300)
cat("✓ Heatmap saved.\n")## ✓ Heatmap saved.Visualisasi Matriks Keanggotaan Fuzzy (Heatmap) secara langsung mengkonfirmasi keberhasilan metode FCMdd-DTW dalam menangkap sifat tumpang tindih (fuzzy) dari data, sekaligus memvisualisasikan dualisme pasar yang ditemukan . Pola heatmap menunjukkan bahwa meskipun beberapa wilayah (seperti Kota Bogor) memiliki keanggotaan yang sangat kuat (mendekati 1.0) pada klaster dominannya—mewakili kasus ideal pola harga stabil—sebagian besar wilayah, terutama di bagian tengah matriks, menunjukkan warna yang sangat berdekatan di kedua kolom. Fenomena ini secara visual membuktikan bahwa banyak Kabupaten/Kota memiliki derajat keanggotaan yang moderat (dekat dengan 0.5) pada Klaster 1 maupun Klaster 2. Dengan demikian, heatmap ini secara tegas memvalidasi bahwa batas antara Klaster 1 dan Klaster 2 sangat kabur dan sebagian besar wilayah berada dalam zona transisi atau kontinuum dualisme pasar, bukan terpisah secara biner.
Distribusi entropy membership
cat("[6/10] Creating Membership Entropy Distribution...\n")## [6/10] Creating Membership Entropy Distribution...# Calculate entropy per object
membership_entropy <- apply(membership, 1, function(u){
u_adj <- pmax(u, .Machine$double.eps)
-sum(u_adj * log(u_adj))
})
cluster_results_fuzzy$Entropy <- membership_entropy
# Calculate max entropy (for reference line)
max_entropy <- -log(1/k_fuzzy)
p6_entropy <- ggplot(cluster_results_fuzzy, aes(x = Entropy)) +
geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7, color = "white") +
geom_vline(xintercept = max_entropy, linetype = "dashed", color = "red", size = 1) +
annotate("text", x = max_entropy, y = Inf,
label = paste0("Max entropy = ", round(max_entropy, 3)),
vjust = 1.5, hjust = -0.1, color = "red") +
theme_minimal() +
labs(title = "Distribution of Membership Entropy",
subtitle = "Higher entropy = more uncertain cluster assignment (closer to uniform distribution)",
x = "Entropy",
y = "Count") +
theme(plot.title = element_text(face = "bold", size = 14))
print(p6_entropy)
ggsave("plots/06_membership_entropy_distribution.png", p6_entropy, width = 10, height = 6, dpi = 300)Plot ini menunjukkan bahwa mayoritas Kabupaten/Kota memiliki nilai entropi yang sangat terkonsentrasi di dekat nilai maksimum (\(0.693\) untuk \(K=2\)). Konsentrasi yang tinggi ini adalah bukti kunci yang membuktikan bahwa FCMdd berhasil menangkap sifat fuzzy dari data: sebagian besar wilayah memiliki keanggotaan yang hampir seragam pada kedua klaster (mendekati \(0.5\)), sehingga penetapan klaster mereka sangat tidak pasti atau berada di zona ambang batas. Fakta ini, bersamaan dengan nilai Average Silhouette yang sangat rendah (\(0.021\)), secara kolektif menegaskan bahwa struktur pasar beras di Jawa Barat tidak dapat dipisahkan secara tegas (crisp). FCMdd terbukti sebagai alat analisis yang tepat karena ia tidak memaksakan batas yang tajam (seperti metode K-Means), melainkan secara akurat mengkonfirmasi bahwa dualisme pasar yang ditemukan (Klaster 1: Stabil dan Klaster 2: Volatil) bersifat gradual dan kontinu, sebuah realitas struktural yang hanya dapat dianalisis secara valid melalui pendekatan fuzzy.
Identifikasi Kasus Ambigu (Borderline Membership)
# Define threshold for "clear" vs "fuzzy" membership
membership_threshold <- 0.7
fuzzy_borderline <- cluster_results_fuzzy %>%
mutate(
Max_Membership = pmax(!!!select(., starts_with("Cluster_"))),
Second_Max = apply(select(., starts_with("Cluster_")), 1, function(x) sort(x, decreasing = TRUE)[2]),
Membership_Gap = Max_Membership - Second_Max,
Classification = case_when(
Max_Membership >= membership_threshold ~ "Clear",
Membership_Gap < 0.2 ~ "Highly Ambiguous",
TRUE ~ "Moderate Ambiguity"
)
) %>%
arrange(Classification, Membership_Gap)
# Summary
fuzzy_summary <- table(fuzzy_borderline$Classification)
print(fuzzy_summary)##
## Clear Highly Ambiguous Moderate Ambiguity
## 2 19 6# List of highly ambiguous cases
ambiguous_cases <- fuzzy_borderline %>%
filter(Classification == "Highly Ambiguous") %>%
select(KabKot, Cluster, Max_Membership, Second_Max, Membership_Gap, Entropy)
print(ambiguous_cases)## KabKot Cluster Max_Membership Second_Max
## Kab. Purwakarta Kab. Purwakarta 1 0.5021205 0.4978795
## Kab. Bandung Kab. Bandung 1 0.5038703 0.4961297
## Kab. Subang Kab. Subang 1 0.5056354 0.4943646
## Kota Banjar Kota Banjar 1 0.5132861 0.4867139
## Kota Cimahi Kota Cimahi 2 0.5156702 0.4843298
## Kota Depok Kota Depok 2 0.5245773 0.4754227
## Kab. Kuningan Kab. Kuningan 2 0.5246131 0.4753869
## Kab. Sukabumi Kab. Sukabumi 2 0.5257743 0.4742257
## Kab. Cirebon Kab. Cirebon 2 0.5341852 0.4658148
## Kab. Ciamis Kab. Ciamis 2 0.5356579 0.4643421
## Kota Cirebon Kota Cirebon 1 0.5376730 0.4623270
## Kota Sukabumi Kota Sukabumi 2 0.5378934 0.4621066
## Kab. Garut Kab. Garut 1 0.5584199 0.4415801
## Kab. Majalengka Kab. Majalengka 2 0.5622129 0.4377871
## Kab. Bekasi Kab. Bekasi 1 0.5723410 0.4276590
## Kab. Sumedang Kab. Sumedang 1 0.5830041 0.4169959
## Kab. Pangandaran Kab. Pangandaran 1 0.5840080 0.4159920
## Kab. Bandung Barat Kab. Bandung Barat 1 0.5858892 0.4141108
## Kota Tasikmalaya Kota Tasikmalaya 1 0.5880363 0.4119637
## Membership_Gap Entropy
## Kab. Purwakarta 0.00424097 0.6931382
## Kab. Bandung 0.00774066 0.6931172
## Kab. Subang 0.01127089 0.6930837
## Kota Banjar 0.02657225 0.6927941
## Kota Cimahi 0.03134043 0.6926560
## Kota Depok 0.04915463 0.6919386
## Kab. Kuningan 0.04922627 0.6919351
## Kab. Sukabumi 0.05154858 0.6918180
## Kab. Cirebon 0.06837033 0.6908081
## Kab. Ciamis 0.07131584 0.6906020
## Kota Cirebon 0.07534599 0.6903060
## Kota Sukabumi 0.07578678 0.6902726
## Kab. Garut 0.11683976 0.6863058
## Kab. Majalengka 0.12442572 0.6853862
## Kab. Bekasi 0.14468205 0.6826439
## Kab. Sumedang 0.16600815 0.6793038
## Kab. Pangandaran 0.16801599 0.6789653
## Kab. Bandung Barat 0.17177843 0.6783198
## Kota Tasikmalaya 0.17607251 0.6775653# Interpretasi:
cat("\n=== FUZZY BORDERLINE ANALYSIS ===\n")##
## === FUZZY BORDERLINE ANALYSIS ===cat("Clear assignment (membership > 0.7):", fuzzy_summary["Clear"],
"kabupaten (", round(fuzzy_summary["Clear"]/nrow(fuzzy_borderline)*100, 1), "%)\n")## Clear assignment (membership > 0.7): 2 kabupaten ( 7.4 %)cat("Highly ambiguous (gap < 0.2):", fuzzy_summary["Highly Ambiguous"],
"kabupaten\n")## Highly ambiguous (gap < 0.2): 19 kabupatencat("\nKabupaten dengan ambiguitas tinggi:\n")##
## Kabupaten dengan ambiguitas tinggi:print(ambiguous_cases$KabKot)## [1] "Kab. Purwakarta" "Kab. Bandung" "Kab. Subang"
## [4] "Kota Banjar" "Kota Cimahi" "Kota Depok"
## [7] "Kab. Kuningan" "Kab. Sukabumi" "Kab. Cirebon"
## [10] "Kab. Ciamis" "Kota Cirebon" "Kota Sukabumi"
## [13] "Kab. Garut" "Kab. Majalengka" "Kab. Bekasi"
## [16] "Kab. Sumedang" "Kab. Pangandaran" "Kab. Bandung Barat"
## [19] "Kota Tasikmalaya"Berdasarkan analisis, mayoritas besar wilayah dikategorikan sebagai Sangat Ambigu (Highly Ambiguous), yang menjadi kunci untuk memahami sifat pasar beras medium Jawa Barat. Dari total 27 Kabupaten/Kota, hanya 2 wilayah (7.4%) yang memiliki pola harga yang sangat jelas cocok dengan Klaster 1 (Stabil) atau Klaster 2 (Volatil), menunjukkan mereka adalah kasus ideal. Sebaliknya, 19 Kabupaten/Kota (70.4%) memiliki selisih kemiripan yang sangat tipis (\(< 0.2\)) antara kedua klaster. Artinya, kelompok mayoritas ini, seperti Kab. Purwakarta dan Kab. Cirebon, sama-sama mirip dengan kedua pola harga ideal tersebut. Angka 70.4% ini menegaskan bahwa sebagian besar pasar berada di area abu-abu atau zona transisi yang tumpang tindih antara kondisi stabil dan volatil. Oleh karena itu, data ini menunjukkan bahwa struktur harga pasar beras adalah kontinu, yang membenarkan pengelompokan yang tidak memaksakan batas yang kaku.
Tabel Statistik Ringkasan Klaster
cat("[10/10] Creating Summary Table...\n")## [10/10] Creating Summary Table...library(grid) # <--- Tambahkan ini
library(gridExtra) # untuk tableGrob, grid.arrange
# Create summary table
summary_table <- cluster_results_fuzzy %>%
group_by(Cluster) %>%
summarise(
N_Members = n(),
Medoid = medoids_fuzzy[as.numeric(Cluster)],
Avg_Membership = mean(membership[, as.numeric(Cluster)]),
Min_Membership = min(membership[, as.numeric(Cluster)]),
Avg_Entropy = mean(Entropy),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~round(., 3)))## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(summary_table)## # A tibble: 27 × 6
## Cluster N_Members Medoid Avg_Membership Min_Membership Avg_Entropy
## <fct> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 15 Kota Bogor 0.508 0 0.637
## 2 1 15 Kota Bogor 0.508 0 0.637
## 3 1 15 Kota Bogor 0.508 0 0.637
## 4 1 15 Kota Bogor 0.508 0 0.637
## 5 1 15 Kota Bogor 0.508 0 0.637
## 6 1 15 Kota Bogor 0.508 0 0.637
## 7 1 15 Kota Bogor 0.508 0 0.637
## 8 1 15 Kota Bogor 0.508 0 0.637
## 9 1 15 Kota Bogor 0.508 0 0.637
## 10 1 15 Kota Bogor 0.508 0 0.637
## # ℹ 17 more rows# Visualisasi tabel sebagai gambar
summary_table_grob <- tableGrob(
summary_table,
rows = NULL,
theme = ttheme_default(
core = list(fg_params = list(cex = 0.8)),
colhead = list(fg_params = list(cex = 0.9, fontface = "bold"))
)
)
p10_summary <- grid.arrange(
summary_table_grob,
top = textGrob("Cluster Summary Statistics",
gp = gpar(fontsize = 16, fontface = "bold"))
)
ggsave("plots/10_cluster_summary_table.png",
p10_summary, width = 12, height = 4, dpi = 300)
cat("✓ Summary table saved\n")## ✓ Summary table savedKesimpulan
Hasil analisis Fuzzy C-Medoids berbasis DTW menunjukkan bahwa dinamika harga beras medium di Jawa Barat dapat dipetakan menjadi dua klaster utama. Klaster pertama dicirikan oleh volatilitas harga yang lebih rendah, dengan koefisien variasi relatif kecil dan rentang harga antar waktu yang sempit. sementara klaster kedua memiliki fluktuasi lebih tinggi dan rentang harga yang lebih lebar. Secara statistik, pemisahan ini menegaskan adanya perbedaan kestabilan antar wilayah dalam menghadapi perubahan pasokan dan permintaan beras medium. Walaupun terdapat dua kelompok, struktur klastering tidak memiliki batas yang kaku. Hal tersebut tercermin dari nilai rata-rata Silhouette yang sangat rendah (0.021) dan tingginya nilai entropi pada sebagian besar wilayah. Data ini mengindikasikan bahwa banyak kabupaten/kota memiliki kemiripan pola terhadap kedua klaster sekaligus, sehingga sifat pasar beras di Jawa Barat lebih tepat disebut kontinu dan overlapped, bukan terfragmentasi. Implikasi praktisnya, sebagian besar wilayah memiliki risiko harga yang tidak sepenuhnya stabil, tetapi juga tidak sepenuhnya rawan, sehingga membutuhkan pemantauan kebijakan yang fleksibel.
keterhubungan sistem harga masih sangat kuat, terlihat dari korelasi tren antar klaster yang mendekati 1 (0.98) meskipun klaster berbeda dalam tingkat volatilitasnya. Secara praktis, ini berarti guncangan pasar di satu wilayah akan cepat menular ke wilayah lain akibat integrasi jalur distribusi dan dominasi sumber pasokan tertentu. Oleh karena itu, kebijakan stabilisasi tidak bisa hanya fokus pada wilayah berisiko tinggi; intervensi pada wilayah borderline dapat berperan sebagai penahan penyebaran volatilitas dalam jaringan pasar. Dalam perspektif metodologi, DTW memberikan keunggulan karena dapat mengukur kemiripan pola meskipun puncak kenaikan harga terjadi pada waktu yang berbeda antar wilayah. Selain itu, medoid sebagai pusat klaster memungkinkan hasil lebih mudah diinterpretasi karena mewakili wilayah nyata, bukan sekadar nilai rata-rata abstrak. Secara praktis, wilayah medoid dapat dijadikan role model per klaster dalam perencanaan stabilisasi harga dan distribusi.
Secara keseluruhan, temuan ini menegaskan bahwa perencanaan kebijakan harga beras di Jawa Barat perlu memperhatikan karakteristik volatilitas masing-masing klaster, serta tingkat ketidakpastian keanggotaan wilayah di dalamnya. Dengan pendekatan berbasis klaster yang adaptif, pemerintah daerah dapat menentukan prioritas intervensi yang lebih presisi, menjaga keterjangkauan harga, dan memperkuat ketahanan pangan secara berkelanjutan.
Analisis Peramalan Deret Waktu per Klaster Menggunakan Model Prophet
Persiapan Library
# ============================================================================
# ANALISIS PERAMALAN HARGA BERAS MENGGUNAKAN PROPHET
# Perbandingan: Baseline (Tuned) vs Regressor (Tuned)
# ============================================================================
library(prophet)## Warning: package 'prophet' was built under R version 4.4.3## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 4.4.2## Loading required package: rlang## Warning: package 'rlang' was built under R version 4.4.2library(dplyr)
library(ggplot2)
library(lubridate)
library(readxl)
library(tidyr)
library(gridExtra)
library(zoo)
cat("✓ Packages loaded\n")## ✓ Packages loadedLoad dan Agregasi Data Hasil Klasterisasi
# ============================================================================
# 1. LOAD DATA HASIL KLASTERISASI
# ============================================================================
# Asumsi: df_clustered sudah ada dari analisis FCMdd-DTW sebelumnya
# df_clustered berisi: KabKot, Periode, Tanggal, Harga, Cluster
# Agregasi data per klaster (rata-rata harga)
df_cluster_ts <- df_clustered %>%
group_by(Cluster, Tanggal) %>%
summarise(Harga_mean = mean(Harga, na.rm = TRUE), .groups = "drop") %>%
arrange(Cluster, Tanggal)
cat("✓ Data aggregated per cluster\n")## ✓ Data aggregated per clusterDefinisi Fungsi Helper untuk Pipeline Modeling
# ============================================================================
# 2. FUNGSI HELPER
# ============================================================================
# 2.1 Prepare data untuk Prophet
prepare_prophet_data <- function(df, date_col, value_col) {
df_prophet <- data.frame(
ds = df[[date_col]],
y = df[[value_col]]
)
return(df_prophet)
}
# 2.2 Split train-test
split_train_test <- function(df_prophet, test_months = 6) {
max_date <- max(df_prophet$ds)
cutoff_date <- max_date - months(test_months)
train <- df_prophet %>% filter(ds <= cutoff_date)
test <- df_prophet %>% filter(ds > cutoff_date)
list(train = train, test = test, cutoff_date = cutoff_date)
}
# 2.3 Evaluasi model
evaluate_forecast <- function(actual, predicted) {
mae <- mean(abs(actual - predicted), na.rm = TRUE)
rmse <- sqrt(mean((actual - predicted)^2, na.rm = TRUE))
mape <- mean(abs((actual - predicted) / actual) * 100, na.rm = TRUE)
data.frame(
MAE = round(mae, 2),
RMSE = round(rmse, 2),
MAPE = round(mape, 2)
)
}
# 2.4 Tambahkan regressors (variabel eksogen)
add_regressors_data <- function(df_prophet) {
df_prophet %>%
mutate(
month = month(ds),
# Dummy variables
early_year = ifelse(month %in% c(2, 3), 1, 0),
mid_year = ifelse(month %in% c(5, 6), 1, 0),
late_year = ifelse(month %in% c(9, 10), 1, 0)
)
}
cat("✓ Helper functions defined\n")## ✓ Helper functions definedFungsi Hyperparameter Tuning untuk Model Prophet
# ============================================================================
# 3. HYPERPARAMETER TUNING FUNCTIONS
# ============================================================================
# 3.1 Tuning untuk BASELINE model (tanpa regressors)
tune_prophet_baseline <- function(df_train, df_test, param_grid, verbose = TRUE) {
results <- data.frame()
total_iterations <- nrow(param_grid)
for (i in 1:total_iterations) {
if(verbose) cat(" Iteration", i, "/", total_iterations, "\r")
# Build model dengan parameter
m <- suppressMessages(prophet(
yearly.seasonality = TRUE,
weekly.seasonality = FALSE,
daily.seasonality = FALSE,
seasonality.mode = 'multiplicative',
changepoint.prior.scale = param_grid$changepoint_prior_scale[i],
seasonality.prior.scale = param_grid$seasonality_prior_scale[i]
))
# Fit (suppress messages)
m <- suppressMessages(fit.prophet(m, df_train))
# Predict
future <- make_future_dataframe(m, periods = nrow(df_test), freq = 'month')
forecast <- suppressMessages(predict(m, future))
# Evaluate
test_pred <- forecast %>%
tail(nrow(df_test)) %>%
pull(yhat)
eval <- evaluate_forecast(df_test$y, test_pred)
# Store
results <- rbind(results, cbind(param_grid[i, ], eval))
}
if(verbose) cat("\n")
return(results)
}
# 3.2 Tuning untuk REGRESSOR model (dengan regressors)
tune_prophet_regressor <- function(df_train, df_test, param_grid, verbose = TRUE) {
results <- data.frame()
total_iterations <- nrow(param_grid)
for (i in 1:total_iterations) {
if(verbose) cat(" Iteration", i, "/", total_iterations, "\r")
# Build model dengan parameter
m <- suppressMessages(prophet(
yearly.seasonality = TRUE,
weekly.seasonality = FALSE,
daily.seasonality = FALSE,
seasonality.mode = 'multiplicative',
changepoint.prior.scale = param_grid$changepoint_prior_scale[i],
seasonality.prior.scale = param_grid$seasonality_prior_scale[i]
))
# Add regressors
m <- add_regressor(m, 'early_year')
m <- add_regressor(m, 'mid_year')
m <- add_regressor(m, 'late_year')
# Fit
m <- suppressMessages(fit.prophet(m, df_train))
# Predict (dengan regressors di future)
future <- make_future_dataframe(m, periods = nrow(df_test), freq = 'month')
future <- add_regressors_data(future)
forecast <- suppressMessages(predict(m, future))
# Evaluate
test_pred <- forecast %>%
tail(nrow(df_test)) %>%
pull(yhat)
eval <- evaluate_forecast(df_test$y, test_pred)
# Store
results <- rbind(results, cbind(param_grid[i, ], eval))
}
if(verbose) cat("\n")
return(results)
}
cat("✓ Tuning functions defined\n")## ✓ Tuning functions definedDefinisi Hyperparameter Grid untuk Grid Search
# ============================================================================
# 4. DEFINE HYPERPARAMETER GRID
# ============================================================================
# Grid search untuk parameter terbaik
param_grid <- expand.grid(
changepoint_prior_scale = c(0.001, 0.01, 0.05, 0.1, 0.5),
seasonality_prior_scale = c(0.01, 0.1, 1.0, 10.0),
stringsAsFactors = FALSE
)
cat("✓ Parameter grid created:", nrow(param_grid), "combinations\n")## ✓ Parameter grid created: 20 combinationsHyperparameter Tuning & Pemodelan Prophet untuk Semua Cluster
# ============================================================================
# 5. MAIN ANALYSIS: TUNING UNTUK SEMUA CLUSTER
# ============================================================================
cat("\n=== STARTING HYPERPARAMETER TUNING FOR ALL CLUSTERS ===\n")##
## === STARTING HYPERPARAMETER TUNING FOR ALL CLUSTERS ===# Storage untuk results
tuning_results_all <- list()
best_models <- list()
comparison_metrics <- data.frame()
clusters <- unique(df_cluster_ts$Cluster)
for (k in clusters) {
cat("\n━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n")
cat("CLUSTER", k, "\n")
cat("━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n")
# ─────────────────────────────────────────────────────
# 5.1 Prepare Data
# ─────────────────────────────────────────────────────
df_k <- df_cluster_ts %>%
filter(Cluster == k) %>%
select(Tanggal, Harga_mean)
# Baseline data
df_prophet_baseline <- prepare_prophet_data(df_k, "Tanggal", "Harga_mean")
split_baseline <- split_train_test(df_prophet_baseline, test_months = 6)
# Regressor data (dengan variabel eksogen)
df_prophet_regressor <- prepare_prophet_data(df_k, "Tanggal", "Harga_mean")
df_prophet_regressor <- add_regressors_data(df_prophet_regressor)
split_regressor <- split_train_test(df_prophet_regressor, test_months = 6)
# ─────────────────────────────────────────────────────
# 5.2 Tuning BASELINE Model
# ─────────────────────────────────────────────────────
cat("\n[1/2] Tuning BASELINE model...\n")
tuning_baseline <- tune_prophet_baseline(
df_train = split_baseline$train,
df_test = split_baseline$test,
param_grid = param_grid,
verbose = TRUE
)
# Best params for baseline
best_baseline <- tuning_baseline %>% arrange(MAPE) %>% head(1)
cat("✓ Best BASELINE parameters:\n")
cat(" - changepoint_prior_scale:", best_baseline$changepoint_prior_scale, "\n")
cat(" - seasonality_prior_scale:", best_baseline$seasonality_prior_scale, "\n")
cat(" - MAPE:", best_baseline$MAPE, "%\n")
# ─────────────────────────────────────────────────────
# 5.3 Tuning REGRESSOR Model
# ─────────────────────────────────────────────────────
cat("\n[2/2] Tuning REGRESSOR model...\n")
tuning_regressor <- tune_prophet_regressor(
df_train = split_regressor$train,
df_test = split_regressor$test,
param_grid = param_grid,
verbose = TRUE
)
# Best params for regressor
best_regressor <- tuning_regressor %>% arrange(MAPE) %>% head(1)
cat("✓ Best REGRESSOR parameters:\n")
cat(" - changepoint_prior_scale:", best_regressor$changepoint_prior_scale, "\n")
cat(" - seasonality_prior_scale:", best_regressor$seasonality_prior_scale, "\n")
cat(" - MAPE:", best_regressor$MAPE, "%\n")
# ─────────────────────────────────────────────────────
# 5.4 Re-fit Models dengan Best Parameters
# ─────────────────────────────────────────────────────
cat("\n[3/3] Fitting final models with best parameters...\n")
# BASELINE FINAL MODEL
model_baseline_final <- suppressMessages(prophet(
yearly.seasonality = TRUE,
weekly.seasonality = FALSE,
daily.seasonality = FALSE,
seasonality.mode = 'multiplicative',
changepoint.prior.scale = best_baseline$changepoint_prior_scale,
seasonality.prior.scale = best_baseline$seasonality_prior_scale
))
model_baseline_final <- suppressMessages(fit.prophet(
model_baseline_final,
split_baseline$train
))
future_baseline <- make_future_dataframe(
model_baseline_final,
periods = nrow(split_baseline$test),
freq = 'month'
)
forecast_baseline <- suppressMessages(predict(model_baseline_final, future_baseline))
# REGRESSOR FINAL MODEL
model_regressor_final <- suppressMessages(prophet(
yearly.seasonality = TRUE,
weekly.seasonality = FALSE,
daily.seasonality = FALSE,
seasonality.mode = 'multiplicative',
changepoint.prior.scale = best_regressor$changepoint_prior_scale,
seasonality.prior.scale = best_regressor$seasonality_prior_scale
))
model_regressor_final <- add_regressor(model_regressor_final, 'early_year')
model_regressor_final <- add_regressor(model_regressor_final, 'mid_year')
model_regressor_final <- add_regressor(model_regressor_final, 'late_year')
model_regressor_final <- suppressMessages(fit.prophet(
model_regressor_final,
split_regressor$train
))
future_regressor <- make_future_dataframe(
model_regressor_final,
periods = nrow(split_regressor$test),
freq = 'month'
)
future_regressor <- add_regressors_data(future_regressor)
forecast_regressor <- suppressMessages(predict(model_regressor_final, future_regressor))
# ─────────────────────────────────────────────────────
# 5.5 Store Results
# ─────────────────────────────────────────────────────
tuning_results_all[[paste0("Cluster_", k)]] <- list(
baseline = list(
tuning_results = tuning_baseline,
best_params = best_baseline,
model = model_baseline_final,
forecast = forecast_baseline,
train = split_baseline$train,
test = split_baseline$test
),
regressor = list(
tuning_results = tuning_regressor,
best_params = best_regressor,
model = model_regressor_final,
forecast = forecast_regressor,
train = split_regressor$train,
test = split_regressor$test
)
)
# Add to comparison metrics
comparison_metrics <- rbind(
comparison_metrics,
data.frame(
Cluster = k,
Model = "Baseline (Tuned)",
best_baseline[, c("changepoint_prior_scale", "seasonality_prior_scale",
"MAE", "RMSE", "MAPE")]
),
data.frame(
Cluster = k,
Model = "Regressor (Tuned)",
best_regressor[, c("changepoint_prior_scale", "seasonality_prior_scale",
"MAE", "RMSE", "MAPE")]
)
)
# ─────────────────────────────────────────────────────
# 5.6 Summary untuk Cluster ini
# ─────────────────────────────────────────────────────
improvement <- ((best_baseline$MAPE - best_regressor$MAPE) / best_baseline$MAPE) * 100
cat("\n📊 SUMMARY for Cluster", k, ":\n")
cat(" Baseline (Tuned) : MAPE =", best_baseline$MAPE, "%\n")
cat(" Regressor (Tuned) : MAPE =", best_regressor$MAPE, "%\n")
if(improvement > 0) {
cat(" ✅ Regressor model is BETTER by", round(improvement, 2), "%\n")
} else {
cat(" ⚠️ Baseline model is BETTER by", round(abs(improvement), 2), "%\n")
}
}##
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
## CLUSTER 1
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
##
## [1/2] Tuning BASELINE model...
## Iteration 1 / 20 Iteration 2 / 20 Iteration 3 / 20 Iteration 4 / 20 Iteration 5 / 20 Iteration 6 / 20 Iteration 7 / 20 Iteration 8 / 20 Iteration 9 / 20 Iteration 10 / 20 Iteration 11 / 20 Iteration 12 / 20 Iteration 13 / 20 Iteration 14 / 20 Iteration 15 / 20 Iteration 16 / 20 Iteration 17 / 20 Iteration 18 / 20 Iteration 19 / 20 Iteration 20 / 20
## ✓ Best BASELINE parameters:
## - changepoint_prior_scale: 0.5
## - seasonality_prior_scale: 1
## - MAPE: 2.08 %
##
## [2/2] Tuning REGRESSOR model...
## Iteration 1 / 20 Iteration 2 / 20 Iteration 3 / 20 Iteration 4 / 20 Iteration 5 / 20 Iteration 6 / 20 Iteration 7 / 20 Iteration 8 / 20 Iteration 9 / 20 Iteration 10 / 20 Iteration 11 / 20 Iteration 12 / 20 Iteration 13 / 20 Iteration 14 / 20 Iteration 15 / 20 Iteration 16 / 20 Iteration 17 / 20 Iteration 18 / 20 Iteration 19 / 20 Iteration 20 / 20
## ✓ Best REGRESSOR parameters:
## - changepoint_prior_scale: 0.001
## - seasonality_prior_scale: 1
## - MAPE: 3.53 %
##
## [3/3] Fitting final models with best parameters...
##
## 📊 SUMMARY for Cluster 1 :
## Baseline (Tuned) : MAPE = 2.08 %
## Regressor (Tuned) : MAPE = 3.53 %
## ⚠️ Baseline model is BETTER by 69.71 %
##
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
## CLUSTER 2
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
##
## [1/2] Tuning BASELINE model...
## Iteration 1 / 20 Iteration 2 / 20 Iteration 3 / 20 Iteration 4 / 20 Iteration 5 / 20 Iteration 6 / 20 Iteration 7 / 20 Iteration 8 / 20 Iteration 9 / 20 Iteration 10 / 20 Iteration 11 / 20 Iteration 12 / 20 Iteration 13 / 20 Iteration 14 / 20 Iteration 15 / 20 Iteration 16 / 20 Iteration 17 / 20 Iteration 18 / 20 Iteration 19 / 20 Iteration 20 / 20
## ✓ Best BASELINE parameters:
## - changepoint_prior_scale: 0.5
## - seasonality_prior_scale: 0.1
## - MAPE: 0.72 %
##
## [2/2] Tuning REGRESSOR model...
## Iteration 1 / 20 Iteration 2 / 20 Iteration 3 / 20 Iteration 4 / 20 Iteration 5 / 20 Iteration 6 / 20 Iteration 7 / 20 Iteration 8 / 20 Iteration 9 / 20 Iteration 10 / 20 Iteration 11 / 20 Iteration 12 / 20 Iteration 13 / 20 Iteration 14 / 20 Iteration 15 / 20 Iteration 16 / 20 Iteration 17 / 20 Iteration 18 / 20 Iteration 19 / 20 Iteration 20 / 20
## ✓ Best REGRESSOR parameters:
## - changepoint_prior_scale: 0.1
## - seasonality_prior_scale: 1
## - MAPE: 1.18 %
##
## [3/3] Fitting final models with best parameters...
##
## 📊 SUMMARY for Cluster 2 :
## Baseline (Tuned) : MAPE = 0.72 %
## Regressor (Tuned) : MAPE = 1.18 %
## ⚠️ Baseline model is BETTER by 63.89 %cat("\n✓ All clusters processed!\n")##
## ✓ All clusters processed!Tabel Ringkasan & Penyimpanan Hasil
# ============================================================================
# 6. SUMMARY TABLE
# ============================================================================
cat("\n=== SUMMARY: COMPARISON METRICS ===\n")##
## === SUMMARY: COMPARISON METRICS ===print(comparison_metrics)## Cluster Model changepoint_prior_scale seasonality_prior_scale
## 1 1 Baseline (Tuned) 0.500 1.0
## 2 1 Regressor (Tuned) 0.001 1.0
## 3 2 Baseline (Tuned) 0.500 0.1
## 4 2 Regressor (Tuned) 0.100 1.0
## MAE RMSE MAPE
## 1 273.79 398.42 2.08
## 2 462.44 568.58 3.53
## 3 96.64 119.72 0.72
## 4 157.04 229.68 1.18# Save to Excel
library(writexl)
write_xlsx(
list(
Comparison_Summary = comparison_metrics,
Baseline_Tuning_Details = do.call(rbind, lapply(names(tuning_results_all), function(k) {
cbind(Cluster = k, tuning_results_all[[k]]$baseline$tuning_results)
})),
Regressor_Tuning_Details = do.call(rbind, lapply(names(tuning_results_all), function(k) {
cbind(Cluster = k, tuning_results_all[[k]]$regressor$tuning_results)
}))
),
"Prophet_Tuned_Comparison_Results.xlsx"
)
cat("✓ Results saved to Prophet_Tuned_Comparison_Results.xlsx\n")## ✓ Results saved to Prophet_Tuned_Comparison_Results.xlsxVisualisasi
# ============================================================================
# 7. VISUALIZATIONS
# ============================================================================
cat("\n=== CREATING VISUALIZATIONS ===\n")##
## === CREATING VISUALIZATIONS ===if(!dir.exists("plots/prophet")) dir.create("plots/prophet", recursive = TRUE)
# ─────────────────────────────────────────────────────────────────────────
# 7.1 COMPARISON BARPLOT: MAPE by Cluster
# ─────────────────────────────────────────────────────────────────────────
p_mape_comparison <- ggplot(comparison_metrics,
aes(x = factor(Cluster), y = MAPE, fill = Model)) +
geom_bar(stat = "identity", position = "dodge", alpha = 0.8, width = 0.7) +
geom_text(aes(label = paste0(round(MAPE, 2), "%")),
position = position_dodge(width = 0.7),
vjust = -0.5, size = 3.5, fontface = "bold") +
scale_fill_manual(values = c("Baseline (Tuned)" = "#E41A1C",
"Regressor (Tuned)" = "#377EB8")) +
theme_minimal() +
labs(title = "Model Comparison: MAPE by Cluster",
subtitle = "Both models tuned with grid search (lower is better)",
x = "Cluster",
y = "MAPE (%)") +
theme(plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12, color = "gray40"),
legend.position = "bottom",
legend.title = element_blank(),
axis.text = element_text(size = 11))
print(p_mape_comparison)
ggsave("plots/prophet/01_mape_comparison_tuned.png",
p_mape_comparison, width = 10, height = 6, dpi = 300)
# ─────────────────────────────────────────────────────────────────────────
# 7.2 COMPARISON: MAE & RMSE
# ─────────────────────────────────────────────────────────────────────────
comparison_long <- comparison_metrics %>%
select(Cluster, Model, MAE, RMSE, MAPE) %>%
pivot_longer(cols = c(MAE, RMSE, MAPE), names_to = "Metric", values_to = "Value")
p_all_metrics <- ggplot(comparison_long,
aes(x = factor(Cluster), y = Value, fill = Model)) +
geom_bar(stat = "identity", position = "dodge", alpha = 0.8) +
facet_wrap(~Metric, scales = "free_y", ncol = 3) +
scale_fill_manual(values = c("Baseline (Tuned)" = "#E41A1C",
"Regressor (Tuned)" = "#377EB8")) +
theme_minimal() +
labs(title = "All Metrics Comparison: Baseline vs Regressor (Both Tuned)",
x = "Cluster", y = "Metric Value") +
theme(plot.title = element_text(face = "bold", size = 14),
legend.position = "bottom",
strip.text = element_text(face = "bold", size = 12))
print(p_all_metrics)
ggsave("plots/prophet/02_all_metrics_comparison.png",
p_all_metrics, width = 14, height = 5, dpi = 300)
# ─────────────────────────────────────────────────────────────────────────
# 7.3 FORECAST PLOTS per Cluster (Baseline vs Regressor)
# ─────────────────────────────────────────────────────────────────────────
for (k in clusters) {
cat("Creating forecast plot for Cluster", k, "...\n")
res <- tuning_results_all[[paste0("Cluster_", k)]]
# Baseline forecast
forecast_baseline <- res$baseline$forecast %>%
select(ds, yhat, yhat_lower, yhat_upper) %>%
mutate(Model = "Baseline (Tuned)")
# Regressor forecast
forecast_regressor <- res$regressor$forecast %>%
select(ds, yhat, yhat_lower, yhat_upper) %>%
mutate(Model = "Regressor (Tuned)")
# Actual data (train + test)
actual_data <- rbind(
res$baseline$train %>% mutate(Set = "Train"),
res$baseline$test %>% mutate(Set = "Test")
)
# Combine forecasts
forecasts_combined <- rbind(forecast_baseline, forecast_regressor)
# Plot
p_forecast <- ggplot() +
# Actual data
geom_point(data = actual_data,
aes(x = ds, y = y, color = Set),
size = 2, alpha = 0.7) +
geom_line(data = actual_data %>% filter(Set == "Train"),
aes(x = ds, y = y),
color = "black", alpha = 0.5) +
# Forecasts
geom_line(data = forecasts_combined,
aes(x = ds, y = yhat, color = Model),
size = 1.2) +
geom_ribbon(data = forecasts_combined,
aes(x = ds, ymin = yhat_lower, ymax = yhat_upper, fill = Model),
alpha = 0.2) +
# Cutoff line
geom_vline(xintercept = res$baseline$test$ds[1],
linetype = "dashed", color = "red", size = 0.8) +
annotate("text", x = res$baseline$test$ds[1], y = Inf,
label = "Test Start", vjust = 1.5, hjust = -0.1, color = "red") +
scale_color_manual(values = c("Train" = "gray40",
"Test" = "black",
"Baseline (Tuned)" = "#E41A1C",
"Regressor (Tuned)" = "#377EB8")) +
scale_fill_manual(values = c("Baseline (Tuned)" = "#E41A1C",
"Regressor (Tuned)" = "#377EB8")) +
theme_minimal() +
labs(title = paste0("Forecast Comparison - Cluster ", k),
subtitle = "Shaded area = 95% confidence interval",
x = "Date", y = "Price (Rp)") +
theme(plot.title = element_text(face = "bold", size = 14),
legend.position = "bottom",
legend.title = element_blank())
print(p_forecast)
ggsave(paste0("plots/prophet/03_forecast_cluster_", k, ".png"),
p_forecast, width = 14, height = 7, dpi = 300)
}## Creating forecast plot for Cluster 1 ...
## Creating forecast plot for Cluster 2 ...
# ─────────────────────────────────────────────────────────────────────────
# 7.4 HYPERPARAMETER HEATMAP (Tuning Results)
# ─────────────────────────────────────────────────────────────────────────
for (k in clusters) {
cat("Creating hyperparameter heatmap for Cluster", k, "...\n")
# Baseline heatmap
tuning_baseline <- tuning_results_all[[paste0("Cluster_", k)]]$baseline$tuning_results
p_heatmap_baseline <- ggplot(tuning_baseline,
aes(x = factor(changepoint_prior_scale),
y = factor(seasonality_prior_scale),
fill = MAPE)) +
geom_tile(color = "white") +
geom_text(aes(label = round(MAPE, 2)), size = 3, fontface = "bold") +
scale_fill_gradient2(low = "green", mid = "yellow", high = "red",
midpoint = median(tuning_baseline$MAPE)) +
theme_minimal() +
labs(title = paste0("Baseline Model Tuning - Cluster ", k),
subtitle = "MAPE across hyperparameter combinations",
x = "changepoint_prior_scale",
y = "seasonality_prior_scale") +
theme(plot.title = element_text(face = "bold", size = 14),
axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(paste0("plots/prophet/04_heatmap_baseline_cluster_", k, ".png"),
p_heatmap_baseline, width = 10, height = 8, dpi = 300)
# Regressor heatmap
tuning_regressor <- tuning_results_all[[paste0("Cluster_", k)]]$regressor$tuning_results
p_heatmap_regressor <- ggplot(tuning_regressor,
aes(x = factor(changepoint_prior_scale),
y = factor(seasonality_prior_scale),
fill = MAPE)) +
geom_tile(color = "white") +
geom_text(aes(label = round(MAPE, 2)), size = 3, fontface = "bold") +
scale_fill_gradient2(low = "green", mid = "yellow", high = "red",
midpoint = median(tuning_regressor$MAPE)) +
theme_minimal() +
labs(title = paste0("Regressor Model Tuning - Cluster ", k),
subtitle = "MAPE across hyperparameter combinations",
x = "changepoint_prior_scale",
y = "seasonality_prior_scale") +
theme(plot.title = element_text(face = "bold", size = 14),
axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(paste0("plots/prophet/05_heatmap_regressor_cluster_", k, ".png"),
p_heatmap_regressor, width = 10, height = 8, dpi = 300)
}## Creating hyperparameter heatmap for Cluster 1 ...## Creating hyperparameter heatmap for Cluster 2 ...# ─────────────────────────────────────────────────────────────────────────
# 7.5 IMPROVEMENT PERCENTAGE PLOT
# ─────────────────────────────────────────────────────────────────────────
improvement_data <- comparison_metrics %>%
select(Cluster, Model, MAPE) %>%
pivot_wider(names_from = Model, values_from = MAPE) %>%
mutate(
Improvement_pct = ((`Baseline (Tuned)` - `Regressor (Tuned)`) / `Baseline (Tuned)`) * 100,
Better_Model = ifelse(Improvement_pct > 0, "Regressor Better", "Baseline Better")
)
p_improvement <- ggplot(improvement_data,
aes(x = factor(Cluster), y = Improvement_pct, fill = Better_Model)) +
geom_bar(stat = "identity", alpha = 0.8) +
geom_text(aes(label = paste0(round(abs(Improvement_pct), 1), "%")),
vjust = ifelse(improvement_data$Improvement_pct > 0, -0.5, 1.5),
fontface = "bold") +
geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 0.8) +
scale_fill_manual(values = c("Regressor Better" = "#4DAF4A",
"Baseline Better" = "#FF7F00")) +
theme_minimal() +
labs(title = "Model Improvement: Regressor vs Baseline",
subtitle = "Positive = Regressor better | Negative = Baseline better",
x = "Cluster",
y = "Improvement (%)") +
theme(plot.title = element_text(face = "bold", size = 14),
legend.position = "bottom",
legend.title = element_blank())
print(p_improvement)
ggsave("plots/prophet/06_improvement_percentage.png",
p_improvement, width = 10, height = 6, dpi = 300)
cat("\n✓✓✓ ALL VISUALIZATIONS COMPLETED! ✓✓✓\n")##
## ✓✓✓ ALL VISUALIZATIONS COMPLETED! ✓✓✓cat("Total plots created: ", 6 + (length(clusters) * 3), "\n")## Total plots created: 12cat("Results saved in: plots/prophet/\n")## Results saved in: plots/prophet/Hasil Forecast
cat("\n=== FORECAST PER CLUSTER ===\n")##
## === FORECAST PER CLUSTER ===for(k in clusters) {
cat("\n━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n")
cat("CLUSTER", k, "\n")
cat("━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n")
res <- tuning_results_all[[paste0("Cluster_", k)]]
# Forecast Baseline
forecast_baseline <- res$baseline$forecast %>%
select(ds, yhat, yhat_lower, yhat_upper) %>%
mutate(Model = "Baseline (Tuned)")
cat("\n-- Baseline (Tuned) Forecast --\n")
print(forecast_baseline)
# Forecast Regressor
forecast_regressor <- res$regressor$forecast %>%
select(ds, yhat, yhat_lower, yhat_upper) %>%
mutate(Model = "Regressor (Tuned)")
cat("\n-- Regressor (Tuned) Forecast --\n")
print(forecast_regressor)
}##
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
## CLUSTER 1
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
##
## -- Baseline (Tuned) Forecast --
## ds yhat yhat_lower yhat_upper Model
## 1 2022-01-01 10276.33 10276.333 10276.33 Baseline (Tuned)
## 2 2022-02-01 10265.93 10265.933 10265.93 Baseline (Tuned)
## 3 2022-03-01 10274.13 10274.133 10274.13 Baseline (Tuned)
## 4 2022-04-01 10187.33 10187.333 10187.33 Baseline (Tuned)
## 5 2022-05-01 10294.07 10294.067 10294.07 Baseline (Tuned)
## 6 2022-06-01 10222.13 10222.133 10222.13 Baseline (Tuned)
## 7 2022-07-01 10203.73 10203.733 10203.73 Baseline (Tuned)
## 8 2022-08-01 10118.07 10118.067 10118.07 Baseline (Tuned)
## 9 2022-09-01 10256.07 10256.067 10256.07 Baseline (Tuned)
## 10 2022-10-01 10297.27 10297.267 10297.27 Baseline (Tuned)
## 11 2022-11-01 10502.60 10502.600 10502.60 Baseline (Tuned)
## 12 2022-12-01 10649.87 10649.867 10649.87 Baseline (Tuned)
## 13 2023-01-01 10979.80 10979.800 10979.80 Baseline (Tuned)
## 14 2023-02-01 11406.87 11406.867 11406.87 Baseline (Tuned)
## 15 2023-03-01 11344.93 11344.933 11344.93 Baseline (Tuned)
## 16 2023-04-01 11331.67 11331.667 11331.67 Baseline (Tuned)
## 17 2023-05-01 11327.00 11327.000 11327.00 Baseline (Tuned)
## 18 2023-06-01 11301.80 11301.800 11301.80 Baseline (Tuned)
## 19 2023-07-01 11287.73 11287.733 11287.73 Baseline (Tuned)
## 20 2023-08-01 11420.73 11420.733 11420.73 Baseline (Tuned)
## 21 2023-09-01 12624.80 12624.800 12624.80 Baseline (Tuned)
## 22 2023-10-01 12791.47 12791.467 12791.47 Baseline (Tuned)
## 23 2023-11-01 12701.07 12701.067 12701.07 Baseline (Tuned)
## 24 2023-12-01 12829.27 12829.267 12829.27 Baseline (Tuned)
## 25 2024-01-01 13145.27 13145.267 13145.27 Baseline (Tuned)
## 26 2024-02-01 14714.53 14714.533 14714.53 Baseline (Tuned)
## 27 2024-03-01 14928.87 14928.867 14928.87 Baseline (Tuned)
## 28 2024-04-01 13680.73 13680.733 13680.73 Baseline (Tuned)
## 29 2024-05-01 12545.13 12545.133 12545.13 Baseline (Tuned)
## 30 2024-06-01 12406.60 12406.600 12406.60 Baseline (Tuned)
## 31 2024-07-01 12700.33 12700.333 12700.33 Baseline (Tuned)
## 32 2024-08-01 12791.87 12791.867 12791.87 Baseline (Tuned)
## 33 2024-09-01 12768.27 12768.267 12768.27 Baseline (Tuned)
## 34 2024-10-01 12754.73 12754.733 12754.73 Baseline (Tuned)
## 35 2024-11-01 12644.07 12644.067 12644.07 Baseline (Tuned)
## 36 2024-12-01 12761.00 12761.000 12761.00 Baseline (Tuned)
## 37 2025-01-01 12952.73 12952.733 12952.73 Baseline (Tuned)
## 38 2025-02-01 12873.13 12873.133 12873.13 Baseline (Tuned)
## 39 2025-03-01 12865.73 12865.733 12865.73 Baseline (Tuned)
## 40 2025-04-01 12864.27 12864.267 12864.27 Baseline (Tuned)
## 41 2025-05-01 12728.05 12586.718 12861.58 Baseline (Tuned)
## 42 2025-06-01 13067.66 12435.125 13654.11 Baseline (Tuned)
## 43 2025-07-01 13401.36 12110.743 14565.14 Baseline (Tuned)
## 44 2025-08-01 13642.00 11607.900 15608.31 Baseline (Tuned)
## 45 2025-09-01 13800.70 10852.583 16665.97 Baseline (Tuned)
## 46 2025-10-01 13902.76 9934.667 17767.13 Baseline (Tuned)
##
## -- Regressor (Tuned) Forecast --
## ds yhat yhat_lower yhat_upper Model
## 1 2022-01-01 10148.18 9517.594 10805.28 Regressor (Tuned)
## 2 2022-02-01 10366.31 9726.250 10978.03 Regressor (Tuned)
## 3 2022-03-01 10313.96 9687.958 10970.93 Regressor (Tuned)
## 4 2022-04-01 10271.47 9683.594 10947.66 Regressor (Tuned)
## 5 2022-05-01 10361.94 9704.436 11035.96 Regressor (Tuned)
## 6 2022-06-01 10305.21 9669.637 10981.46 Regressor (Tuned)
## 7 2022-07-01 10331.09 9628.164 10933.86 Regressor (Tuned)
## 8 2022-08-01 10338.50 9727.521 10978.56 Regressor (Tuned)
## 9 2022-09-01 10826.87 10202.088 11505.91 Regressor (Tuned)
## 10 2022-10-01 10867.82 10231.726 11526.72 Regressor (Tuned)
## 11 2022-11-01 10848.21 10220.936 11530.35 Regressor (Tuned)
## 12 2022-12-01 10978.79 10307.938 11611.31 Regressor (Tuned)
## 13 2023-01-01 11306.61 10686.537 11936.63 Regressor (Tuned)
## 14 2023-02-01 12067.07 11441.833 12720.92 Regressor (Tuned)
## 15 2023-03-01 10797.56 10152.329 11430.72 Regressor (Tuned)
## 16 2023-04-01 11175.39 10545.323 11832.45 Regressor (Tuned)
## 17 2023-05-01 11560.84 10926.369 12201.73 Regressor (Tuned)
## 18 2023-06-01 11486.90 10856.157 12131.74 Regressor (Tuned)
## 19 2023-07-01 11339.34 10679.544 11994.38 Regressor (Tuned)
## 20 2023-08-01 11240.64 10612.366 11887.09 Regressor (Tuned)
## 21 2023-09-01 11963.78 11301.861 12590.52 Regressor (Tuned)
## 22 2023-10-01 11936.99 11284.784 12594.47 Regressor (Tuned)
## 23 2023-11-01 11880.08 11254.393 12509.09 Regressor (Tuned)
## 24 2023-12-01 12030.08 11428.819 12657.64 Regressor (Tuned)
## 25 2024-01-01 12482.85 11819.587 13102.59 Regressor (Tuned)
## 26 2024-02-01 13879.89 13239.202 14519.90 Regressor (Tuned)
## 27 2024-03-01 13782.40 13153.302 14402.20 Regressor (Tuned)
## 28 2024-04-01 12788.70 12123.767 13419.76 Regressor (Tuned)
## 29 2024-05-01 12229.22 11586.163 12897.72 Regressor (Tuned)
## 30 2024-06-01 12127.60 11521.318 12781.40 Regressor (Tuned)
## 31 2024-07-01 12509.32 11806.025 13134.91 Regressor (Tuned)
## 32 2024-08-01 12734.50 12083.104 13348.33 Regressor (Tuned)
## 33 2024-09-01 12850.27 12215.548 13477.80 Regressor (Tuned)
## 34 2024-10-01 13011.62 12360.107 13640.17 Regressor (Tuned)
## 35 2024-11-01 13061.59 12441.492 13726.38 Regressor (Tuned)
## 36 2024-12-01 13166.41 12525.335 13831.18 Regressor (Tuned)
## 37 2025-01-01 13315.83 12666.783 13958.12 Regressor (Tuned)
## 38 2025-02-01 12936.76 12314.212 13568.21 Regressor (Tuned)
## 39 2025-03-01 14273.25 13630.061 14908.76 Regressor (Tuned)
## 40 2025-04-01 13681.74 13012.753 14355.00 Regressor (Tuned)
## 41 2025-05-01 13420.74 12737.885 14075.49 Regressor (Tuned)
## 42 2025-06-01 13314.39 12671.054 13962.66 Regressor (Tuned)
## 43 2025-07-01 13521.09 12851.406 14136.76 Regressor (Tuned)
## 44 2025-08-01 13627.03 12988.393 14281.11 Regressor (Tuned)
## 45 2025-09-01 13994.84 13334.066 14708.28 Regressor (Tuned)
## 46 2025-10-01 14093.43 13427.994 14709.43 Regressor (Tuned)
##
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
## CLUSTER 2
## ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
##
## -- Baseline (Tuned) Forecast --
## ds yhat yhat_lower yhat_upper Model
## 1 2022-01-01 10067.250 10067.250 10067.250 Baseline (Tuned)
## 2 2022-02-01 10126.917 10126.917 10126.917 Baseline (Tuned)
## 3 2022-03-01 10164.000 10164.000 10164.000 Baseline (Tuned)
## 4 2022-04-01 10117.250 10117.250 10117.250 Baseline (Tuned)
## 5 2022-05-01 10058.167 10058.167 10058.167 Baseline (Tuned)
## 6 2022-06-01 10026.583 10026.583 10026.583 Baseline (Tuned)
## 7 2022-07-01 10108.667 10108.667 10108.667 Baseline (Tuned)
## 8 2022-08-01 9893.583 9893.583 9893.583 Baseline (Tuned)
## 9 2022-09-01 10011.083 10011.083 10011.083 Baseline (Tuned)
## 10 2022-10-01 10135.750 10135.750 10135.750 Baseline (Tuned)
## 11 2022-11-01 10171.167 10171.167 10171.167 Baseline (Tuned)
## 12 2022-12-01 10549.250 10549.250 10549.250 Baseline (Tuned)
## 13 2023-01-01 10712.250 10712.250 10712.250 Baseline (Tuned)
## 14 2023-02-01 11138.667 11138.667 11138.667 Baseline (Tuned)
## 15 2023-03-01 10974.583 10974.583 10974.583 Baseline (Tuned)
## 16 2023-04-01 10922.583 10922.583 10922.583 Baseline (Tuned)
## 17 2023-05-01 10889.333 10889.333 10889.333 Baseline (Tuned)
## 18 2023-06-01 10923.833 10923.833 10923.833 Baseline (Tuned)
## 19 2023-07-01 10918.833 10918.833 10918.833 Baseline (Tuned)
## 20 2023-08-01 11183.917 11183.917 11183.917 Baseline (Tuned)
## 21 2023-09-01 12701.667 12701.667 12701.667 Baseline (Tuned)
## 22 2023-10-01 12583.583 12583.583 12583.583 Baseline (Tuned)
## 23 2023-11-01 12509.500 12509.500 12509.500 Baseline (Tuned)
## 24 2023-12-01 12402.167 12402.167 12402.167 Baseline (Tuned)
## 25 2024-01-01 12716.583 12716.583 12716.583 Baseline (Tuned)
## 26 2024-02-01 14294.417 14294.417 14294.417 Baseline (Tuned)
## 27 2024-03-01 14224.917 14224.917 14224.917 Baseline (Tuned)
## 28 2024-04-01 13266.083 13266.083 13266.083 Baseline (Tuned)
## 29 2024-05-01 12754.167 12754.167 12754.167 Baseline (Tuned)
## 30 2024-06-01 12597.000 12597.000 12597.000 Baseline (Tuned)
## 31 2024-07-01 12869.750 12869.750 12869.750 Baseline (Tuned)
## 32 2024-08-01 13011.750 13011.750 13011.750 Baseline (Tuned)
## 33 2024-09-01 12962.333 12962.333 12962.333 Baseline (Tuned)
## 34 2024-10-01 12959.667 12959.667 12959.667 Baseline (Tuned)
## 35 2024-11-01 12907.583 12907.583 12907.583 Baseline (Tuned)
## 36 2024-12-01 12901.167 12901.167 12901.167 Baseline (Tuned)
## 37 2025-01-01 13194.083 13194.083 13194.083 Baseline (Tuned)
## 38 2025-02-01 13311.917 13311.917 13311.917 Baseline (Tuned)
## 39 2025-03-01 13294.250 13294.250 13294.250 Baseline (Tuned)
## 40 2025-04-01 13252.917 13252.917 13252.917 Baseline (Tuned)
## 41 2025-05-01 13078.100 12915.697 13225.664 Baseline (Tuned)
## 42 2025-06-01 13164.202 12490.236 13831.365 Baseline (Tuned)
## 43 2025-07-01 13323.837 12075.797 14707.539 Baseline (Tuned)
## 44 2025-08-01 13500.158 11353.104 15710.517 Baseline (Tuned)
## 45 2025-09-01 13504.145 10416.317 16760.226 Baseline (Tuned)
## 46 2025-10-01 13539.796 9437.040 17954.804 Baseline (Tuned)
##
## -- Regressor (Tuned) Forecast --
## ds yhat yhat_lower yhat_upper Model
## 1 2022-01-01 10067.250 10067.250 10067.250 Regressor (Tuned)
## 2 2022-02-01 10126.917 10126.917 10126.917 Regressor (Tuned)
## 3 2022-03-01 10164.000 10164.000 10164.000 Regressor (Tuned)
## 4 2022-04-01 10117.250 10117.250 10117.250 Regressor (Tuned)
## 5 2022-05-01 10058.167 10058.167 10058.167 Regressor (Tuned)
## 6 2022-06-01 10026.583 10026.583 10026.583 Regressor (Tuned)
## 7 2022-07-01 10108.667 10108.667 10108.667 Regressor (Tuned)
## 8 2022-08-01 9893.583 9893.583 9893.583 Regressor (Tuned)
## 9 2022-09-01 10011.083 10011.083 10011.083 Regressor (Tuned)
## 10 2022-10-01 10135.750 10135.750 10135.750 Regressor (Tuned)
## 11 2022-11-01 10171.167 10171.167 10171.167 Regressor (Tuned)
## 12 2022-12-01 10549.250 10549.250 10549.250 Regressor (Tuned)
## 13 2023-01-01 10712.250 10712.250 10712.250 Regressor (Tuned)
## 14 2023-02-01 11138.667 11138.667 11138.667 Regressor (Tuned)
## 15 2023-03-01 10974.583 10974.583 10974.583 Regressor (Tuned)
## 16 2023-04-01 10922.583 10922.583 10922.583 Regressor (Tuned)
## 17 2023-05-01 10889.333 10889.333 10889.333 Regressor (Tuned)
## 18 2023-06-01 10923.833 10923.833 10923.833 Regressor (Tuned)
## 19 2023-07-01 10918.833 10918.833 10918.833 Regressor (Tuned)
## 20 2023-08-01 11183.917 11183.917 11183.917 Regressor (Tuned)
## 21 2023-09-01 12701.667 12701.667 12701.667 Regressor (Tuned)
## 22 2023-10-01 12583.583 12583.583 12583.583 Regressor (Tuned)
## 23 2023-11-01 12509.500 12509.500 12509.500 Regressor (Tuned)
## 24 2023-12-01 12402.167 12402.167 12402.167 Regressor (Tuned)
## 25 2024-01-01 12716.583 12716.583 12716.583 Regressor (Tuned)
## 26 2024-02-01 14294.417 14294.417 14294.417 Regressor (Tuned)
## 27 2024-03-01 14224.917 14224.917 14224.917 Regressor (Tuned)
## 28 2024-04-01 13266.083 13266.083 13266.083 Regressor (Tuned)
## 29 2024-05-01 12754.167 12754.167 12754.167 Regressor (Tuned)
## 30 2024-06-01 12597.000 12597.000 12597.000 Regressor (Tuned)
## 31 2024-07-01 12869.750 12869.750 12869.750 Regressor (Tuned)
## 32 2024-08-01 13011.750 13011.750 13011.750 Regressor (Tuned)
## 33 2024-09-01 12962.333 12962.333 12962.333 Regressor (Tuned)
## 34 2024-10-01 12959.667 12959.667 12959.667 Regressor (Tuned)
## 35 2024-11-01 12907.583 12907.583 12907.583 Regressor (Tuned)
## 36 2024-12-01 12901.167 12901.167 12901.167 Regressor (Tuned)
## 37 2025-01-01 13194.083 13194.083 13194.083 Regressor (Tuned)
## 38 2025-02-01 13311.917 13311.917 13311.917 Regressor (Tuned)
## 39 2025-03-01 13294.250 13294.250 13294.250 Regressor (Tuned)
## 40 2025-04-01 13252.917 13252.917 13252.917 Regressor (Tuned)
## 41 2025-05-01 13128.516 12964.152 13278.074 Regressor (Tuned)
## 42 2025-06-01 13237.160 12605.809 13863.011 Regressor (Tuned)
## 43 2025-07-01 13421.258 12094.371 14674.704 Regressor (Tuned)
## 44 2025-08-01 13586.035 11512.345 15672.167 Regressor (Tuned)
## 45 2025-09-01 13663.047 10607.431 16663.700 Regressor (Tuned)
## 46 2025-10-01 13777.369 9842.807 17732.571 Regressor (Tuned)Pembahasan analisis peramalan
Analisis peramalan harga beras medium ini dilakukan menggunakan metode Prophet berdasarkan data yang telah dikelompokkan ke dalam beberapa cluster melalui analisis FCMdd-DTW. Setiap cluster mewakili pola harga yang relatif seragam sehingga memudahkan dalam memprediksi tren harga secara lebih spesifik. Proses analisis dibagi menjadi dua pendekatan model yang berbeda, yaitu Baseline dan Regressor, yang keduanya melalui proses tuning untuk mendapatkan kombinasi parameter terbaik, khususnya changepoint_prior_scale yang mengatur sensitivitas model terhadap perubahan tren, dan seasonality_prior_scale yang mengontrol pengaruh pola musiman.
Model Baseline hanya menggunakan data historis harga beras. Dengan kata lain, model ini mempelajari pola harga dari masa lalu tanpa informasi tambahan, sehingga prediksi yang dihasilkan sepenuhnya didasarkan pada tren dan siklus musiman yang ada pada data. Tuning dilakukan untuk memastikan model responsif terhadap perubahan harga namun tidak berlebihan sehingga prediksi tetap stabil dan akurat.
Sementara itu, model Regressor menambahkan variabel eksogen berupa indikator periode tertentu dalam setahun untuk menangkap fluktuasi harga yang bersifat musiman tetapi tidak sepenuhnya tercermin dari data historis. Variabel ini dibuat sebagai dummy variables, yaitu early_year, mid_year, dan late_year, yang masing-masing mewakili bulan-bulan awal, pertengahan, dan akhir tahun. Nilai dummy ini bernilai 1 jika data berada pada periode tersebut, dan 0 jika tidak. Dengan menambahkan variabel dummy ini sebagai regresor ke dalam model Prophet, prediksi tidak hanya memanfaatkan pola historis tetapi juga menyesuaikan dengan periode tertentu yang biasanya mengalami perubahan harga yang signifikan, sehingga model diharapkan lebih adaptif terhadap pola musiman spesifik.
Hasil analisis menunjukkan perbedaan performa yang jelas antara kedua model untuk tiap cluster. Pada Cluster 1, model Baseline yang hanya mengandalkan data historis menghasilkan MAPE sebesar 2,08% dengan MAE 273,79 dan RMSE 398,42, sementara model Regressor dengan variabel dummy justru menunjukkan MAPE lebih tinggi, yaitu 3,53%, MAE 462,44 dan RMSE 568,58. Hal ini menandakan bahwa pada cluster ini, pola harga historis sudah cukup untuk memprediksi harga dengan akurat, sehingga tambahan variabel dummy tidak memberikan keuntungan, bahkan menurunkan akurasi. Pada Cluster 2, hasil serupa terlihat di mana Baseline menghasilkan MAPE 0,72%, MAE 96,64, dan RMSE 119,72, sedangkan Regressor memiliki MAPE 1,18%, MAE 157,04, dan RMSE 229,68. Sekali lagi, penambahan variabel eksogen tidak meningkatkan performa karena pola historis di cluster ini cukup kuat untuk menangkap tren dan musiman harga.
Secara keseluruhan, analisis ini menegaskan bahwa model Prophet yang sederhana namun dituning dengan baik dapat memberikan prediksi harga beras yang akurat. Variabel dummy pada model Regressor memiliki tujuan untuk menangkap fluktuasi harga yang spesifik pada periode tertentu, namun kebermanfaatannya sangat tergantung pada karakteristik tiap cluster. Pada kasus ini, data historis sendiri sudah cukup untuk menangkap pola harga secara efektif, sehingga Baseline memberikan hasil yang lebih andal dibandingkan Regressor. Proses tuning tetap krusial untuk memastikan model adaptif terhadap perubahan tren, dan hasil per cluster menunjukkan bahwa strategi prediksi harus disesuaikan dengan karakteristik pola harga masing-masing cluster untuk mencapai akurasi terbaik.