Fluktuasi harga bawang merah di Indonesia dipengaruhi oleh musim, pasokan, dan kondisi wilayah sehingga pola harganya berbeda-beda antar daerah. Untuk mengidentifikasi kesamaan pola tersebut, digunakan metode Fuzzy C-Means (FCM) yang mampu mengelompokkan kabupaten/kota secara fleksibel berdasarkan data harga bawang merah 2022–2024. Analisis ini membantu memetakan wilayah dengan dinamika harga serupa sebagai dasar perencanaan kebijakan pangan.

1 Import Library and Data

library(readr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(dplyr)
library(ggplot2)
library(tidyr)
library(ppclust)
library(cluster)
library(clValid)
library(sf)

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

library(tmap)
library(stringr)
library(patchwork)
library(clusterSim)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:patchwork':
## 
##     area

## The following object is masked from 'package:dplyr':
## 
##     select

df <- read_delim("~/Aliyah/STK/TPG/0_df_long.csv")

## Rows: 340620 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (5): KabKot, Bulan, Produk, Kategori, NamaProv
## dbl (4): Tahun, Harga, KodeBPS, KodeProv
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df)

## # A tibble: 6 × 9
##   KabKot          Tahun Bulan   Produk  Harga Kategori KodeBPS KodeProv NamaProv
##   <chr>           <dbl> <chr>   <chr>   <dbl> <chr>      <dbl>    <dbl> <chr>   
## 1 Kab. Aceh Barat  2022 Januari Beras … 11429 Beras       11.0       11 Aceh    
## 2 Kab. Aceh Barat  2022 Januari Beras …  9979 Beras       11.0       11 Aceh    
## 3 Kab. Aceh Barat  2022 Januari Bawang… 31000 Bawang      11.0       11 Aceh    
## 4 Kab. Aceh Barat  2022 Januari Bawang… 28636 Bawang      11.0       11 Aceh    
## 5 Kab. Aceh Barat  2022 Januari Cabai … 19409 Cabai       11.0       11 Aceh    
## 6 Kab. Aceh Barat  2022 Januari Cabai … 45706 Cabai       11.0       11 Aceh

colSums(is.na(df))

##   KabKot    Tahun    Bulan   Produk    Harga Kategori  KodeBPS KodeProv 
##        0        0        0        0    61513        0        0        0 
## NamaProv 
##     3450

Terdapat total 61513 missing value pada kolom Harga

table(df$Tahun)

## 
##  2022  2023  2024  2025 
## 88275 88815 89730 73800

table(df$Produk)

## 
##           Bawang Merah Bawang Putih (Bonggol)           Beras Medium 
##                  22708                  22708                  22708 
##          Beras Premium      Cabai Merah Besar   Cabai Merah Keriting 
##                  22708                  22708                  22708 
##      Cabai Rawit Merah        Daging Ayam Ras      Daging Sapi Murni 
##                  22708                  22708                  22708 
##          Gula Konsumsi    Jagung Tk. Peternak    Minyak Goreng Curah 
##                  22708                  22708                  22708 
##  Minyak Goreng Kemasan              Minyakita         Telur Ayam Ras 
##                  22708                  22708                  22708

Dataset ini memuat beragam komoditas pangan utama di Indonesia—mulai dari bawang, beras, cabai, daging, gula, jagung, minyak goreng, hingga telur ayam—dengan masing-masing memiliki jumlah pengamatan sebanyak 22.708 data.

1.1 Shapefile

Shapefile batas administrasi kabupaten/kota di Indonesia diimpor untuk keperluan visualisasi peta klaster.

indonesia_shp <- st_read("~/Aliyah/LOMBA/09 SMATIC 2025/datashapefile/idn_admbnda_adm2_bps_20200401.shp")

## Reading layer `idn_admbnda_adm2_bps_20200401' from data source 
##   `C:\Users\User\Documents\Aliyah\LOMBA\09 SMATIC 2025\datashapefile\idn_admbnda_adm2_bps_20200401.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 522 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84

2 Pre-Processing Data

2.1 Filter Bawang Merah

bulan_id <- 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_sub <- df %>%
  filter(Produk == "Bawang Merah") %>%
  mutate(
    BulanNum = bulan_id[Bulan])

head(df_sub)

## # A tibble: 6 × 10
##   KabKot    Tahun Bulan Produk Harga Kategori KodeBPS KodeProv NamaProv BulanNum
##   <chr>     <dbl> <chr> <chr>  <dbl> <chr>      <dbl>    <dbl> <chr>       <dbl>
## 1 Kab. Ace…  2022 Janu… Bawan… 31000 Bawang      11.0       11 Aceh            1
## 2 Kab. Ace…  2022 Febr… Bawan… 33107 Bawang      11.0       11 Aceh            2
## 3 Kab. Ace…  2022 Maret Bawan… 37600 Bawang      11.0       11 Aceh            3
## 4 Kab. Ace…  2022 April Bawan… 35800 Bawang      11.0       11 Aceh            4
## 5 Kab. Ace…  2022 Mei   Bawan… 41600 Bawang      11.0       11 Aceh            5
## 6 Kab. Ace…  2022 Juni  Bawan… 57481 Bawang      11.0       11 Aceh            6

2.2 Ubah ke Format Wide

df_wide_raw <- df_sub %>%
  mutate(
    Bulan = factor(
      BulanNum,
      levels = 1:12,
      labels = month.abb
    )
  ) %>%
  dplyr::select(KabKot, KodeProv, Tahun, Bulan, Harga) %>%
  pivot_wider(
    names_from = Bulan,
    values_from = Harga
  ) %>%
  arrange(KabKot, Tahun)

head(df_wide_raw)

## # A tibble: 6 × 15
##   KabKot    KodeProv Tahun   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
##   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Ace…       11  2022 31000 33107 37600 35800 41600 57481 57233 35567 32367
## 2 Kab. Ace…       11  2023 35500 38931 34167 31760 34065 35321 36759 31655 25071
## 3 Kab. Ace…       11  2024 38386 37140 38623 55768 56542 50893 38508 28000 29200
## 4 Kab. Ace…       11  2025 44640 40187 36688 47461 43941 42325 48045 59456 47663
## 5 Kab. Ace…       11  2022 31231 31913 37731 36962 32860 51313 53065 33800 29400
## 6 Kab. Ace…       11  2023 26273 26897 33871 33040 25167 25071 25032 25000 24900
## # ℹ 3 more variables: Oct <dbl>, Nov <dbl>, Dec <dbl>

df_long1 <- df_wide_raw %>%
  pivot_longer(
    cols = Jan:Dec,
    names_to = "Bulan",
    values_to = "Nilai"
  )

df_long1 <- df_long1 %>%
  mutate(BulanTahun = paste0(Bulan, "_", Tahun))

df_wide1 <- df_long1 %>%
  dplyr::select(KabKot, KodeProv, BulanTahun, Nilai) %>%
  pivot_wider(
    names_from = BulanTahun,
    values_from = Nilai
  ) %>%
  arrange(KodeProv, KabKot)

colSums(is.na(df_wide1))

##   KabKot KodeProv Jan_2022 Feb_2022 Mar_2022 Apr_2022 May_2022 Jun_2022 
##        0        0       17       17       16       17       18       16 
## Jul_2022 Aug_2022 Sep_2022 Oct_2022 Nov_2022 Dec_2022 Jan_2023 Feb_2023 
##       16       16       16       16       17       13       13       13 
## Mar_2023 Apr_2023 May_2023 Jun_2023 Jul_2023 Aug_2023 Sep_2023 Oct_2023 
##       14       14       14       14       14       13       13       13 
## Nov_2023 Dec_2023 Jan_2024 Feb_2024 Mar_2024 Apr_2024 May_2024 Jun_2024 
##       13       13        8        9       10       11        8        6 
## Jul_2024 Aug_2024 Sep_2024 Oct_2024 Nov_2024 Dec_2024 Jan_2025 Feb_2025 
##        6        6        7        9        9       10       25       19 
## Mar_2025 Apr_2025 May_2025 Jun_2025 Jul_2025 Aug_2025 Sep_2025 Oct_2025 
##       19       17       15       14       12       12        9       11 
## Nov_2025 Dec_2025 
##      506      506

Bulan November dan Desember menunjukkan jumlah missing value yang sangat tinggi karena data aktual hanya tersedia hingga Oktober 2025. Selanjutnya dilakukan filter hanya tahun 2022-2024 saja dan imputasi dengan menggunakan rata-rata tiap provinsi pada tahun dan bulan yang sama.

2.3 Imputasi

df_long_for_impute <- df_wide1 %>%
  dplyr::select(-Nov_2025, -Dec_2025) %>%
  pivot_longer(
    cols = Jan_2022:Oct_2025,
    names_to = "Bulan_Tahun",
    values_to = "Harga"
  ) %>%
  separate(
    col = Bulan_Tahun,
    into = c("Bulan", "Tahun"),
    sep = "_",
    remove = FALSE
  )  

head(df_long_for_impute)

## # A tibble: 6 × 6
##   KabKot          KodeProv Bulan_Tahun Bulan Tahun Harga
##   <chr>              <dbl> <chr>       <chr> <chr> <dbl>
## 1 Kab. Aceh Barat       11 Jan_2022    Jan   2022  31000
## 2 Kab. Aceh Barat       11 Feb_2022    Feb   2022  33107
## 3 Kab. Aceh Barat       11 Mar_2022    Mar   2022  37600
## 4 Kab. Aceh Barat       11 Apr_2022    Apr   2022  35800
## 5 Kab. Aceh Barat       11 May_2022    May   2022  41600
## 6 Kab. Aceh Barat       11 Jun_2022    Jun   2022  57481

colSums(is.na(df_long_for_impute))

##      KabKot    KodeProv Bulan_Tahun       Bulan       Tahun       Harga 
##           0           0           0           0           0         608

df_impute <- df_long_for_impute %>%
  filter(Tahun <= 2024) %>% 
  group_by(KodeProv, Tahun, Bulan) %>% 
  mutate(
    Harga = replace_na(Harga, mean(Harga, na.rm = TRUE))
  ) %>%
  ungroup()

head(df_impute)

## # A tibble: 6 × 6
##   KabKot          KodeProv Bulan_Tahun Bulan Tahun Harga
##   <chr>              <dbl> <chr>       <chr> <chr> <dbl>
## 1 Kab. Aceh Barat       11 Jan_2022    Jan   2022  31000
## 2 Kab. Aceh Barat       11 Feb_2022    Feb   2022  33107
## 3 Kab. Aceh Barat       11 Mar_2022    Mar   2022  37600
## 4 Kab. Aceh Barat       11 Apr_2022    Apr   2022  35800
## 5 Kab. Aceh Barat       11 May_2022    May   2022  41600
## 6 Kab. Aceh Barat       11 Jun_2022    Jun   2022  57481

colSums(is.na(df_impute))

##      KabKot    KodeProv Bulan_Tahun       Bulan       Tahun       Harga 
##           0           0           0           0           0           0

df_wide_final <- df_impute %>%
  dplyr::select(KabKot, KodeProv, Tahun, Bulan, Harga) %>%
  pivot_wider(
    names_from = Bulan,
    values_from = Harga
  ) %>%
  arrange(KabKot, Tahun)
head(df_wide_final)

## # A tibble: 6 × 15
##   KabKot    KodeProv Tahun   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
##   <chr>        <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Ace…       11 2022  31000 33107 37600 35800 41600 57481 57233 35567 32367
## 2 Kab. Ace…       11 2023  35500 38931 34167 31760 34065 35321 36759 31655 25071
## 3 Kab. Ace…       11 2024  38386 37140 38623 55768 56542 50893 38508 28000 29200
## 4 Kab. Ace…       11 2022  31231 31913 37731 36962 32860 51313 53065 33800 29400
## 5 Kab. Ace…       11 2023  26273 26897 33871 33040 25167 25071 25032 25000 24900
## 6 Kab. Ace…       11 2024  30419 30828 33484 50630 54419 50148 33067 24700 22357
## # ℹ 3 more variables: Oct <dbl>, Nov <dbl>, Dec <dbl>

colSums(is.na(df_wide_final))

##   KabKot KodeProv    Tahun      Jan      Feb      Mar      Apr      May 
##        0        0        0        0        0        0        0        0 
##      Jun      Jul      Aug      Sep      Oct      Nov      Dec 
##        0        0        0        0        0        0        0

Setelah proses imputasi dilakukan, seluruh variabel, termasuk bulan Januari hingga Desember, sudah tidak memiliki missing value, sehingga dataset kini sepenuhnya lengkap dan siap untuk dianalisis lebih lanjut.

3 Eksplorasi Data

colnames(df_wide_final)

##  [1] "KabKot"   "KodeProv" "Tahun"    "Jan"      "Feb"      "Mar"     
##  [7] "Apr"      "May"      "Jun"      "Jul"      "Aug"      "Sep"     
## [13] "Oct"      "Nov"      "Dec"

# 1. Filter hanya tahun 2022
data_2022 <- df_wide_final %>% 
  filter(Tahun == 2022)

# 2. Hitung rata-rata setiap bulan
mean_2022 <- data_2022 %>% 
  summarise(
    Jan = mean(Jan, na.rm = TRUE),
    Feb = mean(Feb, na.rm = TRUE),
    Mar = mean(Mar, na.rm = TRUE),
    Apr = mean(Apr, na.rm = TRUE),
    May = mean(May, na.rm = TRUE),
    Jun = mean(Jun, na.rm = TRUE),
    Jul = mean(Jul, na.rm = TRUE),
    Aug = mean(Aug, na.rm = TRUE),
    Sep = mean(Sep, na.rm = TRUE),
    Oct = mean(Oct, na.rm = TRUE),
    Nov = mean(Nov, na.rm = TRUE),
    Dec = mean(Dec, na.rm = TRUE)
  )

# 3. Ubah wide → long untuk plot
mean_long <- mean_2022 %>% 
  pivot_longer(cols = everything(),
               names_to = "Bulan",
               values_to = "Rata_rata")

# 4. Urutkan bulan
mean_long$Bulan <- factor(mean_long$Bulan,
                          levels = c("Jan","Feb","Mar","Apr","May",
                                     "Jun","Jul","Aug","Sep","Oct","Nov","Dec"))


# 4. Line chart
ggplot(mean_long, aes(x = Bulan, y = Rata_rata, group = 1)) +
  geom_line() +
  geom_point() +
  labs(title = "Rata-rata Harga Bawang Merah per Bulan (2022)",
       x = "Bulan",
       y = "Rata-rata Harga") +
  theme_minimal()

Harga bawang merah pada tahun 2022 meningkat pada bulan Juni-Juli, kemudian mulai menurun kembali di bulan Agustus-Oktober.

# 1. Filter hanya tahun 2022
data_2023 <- df_wide_final %>% 
  filter(Tahun == 2023)

# 2. Hitung rata-rata setiap bulan
mean_2023 <- data_2023 %>% 
  summarise(
    Jan = mean(Jan, na.rm = TRUE),
    Feb = mean(Feb, na.rm = TRUE),
    Mar = mean(Mar, na.rm = TRUE),
    Apr = mean(Apr, na.rm = TRUE),
    May = mean(May, na.rm = TRUE),
    Jun = mean(Jun, na.rm = TRUE),
    Jul = mean(Jul, na.rm = TRUE),
    Aug = mean(Aug, na.rm = TRUE),
    Sep = mean(Sep, na.rm = TRUE),
    Oct = mean(Oct, na.rm = TRUE),
    Nov = mean(Nov, na.rm = TRUE),
    Dec = mean(Dec, na.rm = TRUE)
  )

# 3. Ubah wide → long untuk plot
mean_long <- mean_2023 %>% 
  pivot_longer(cols = everything(),
               names_to = "Bulan",
               values_to = "Rata_rata")

# 4. Urutkan bulan
mean_long$Bulan <- factor(mean_long$Bulan,
                          levels = c("Jan","Feb","Mar","Apr","May",
                                     "Jun","Jul","Aug","Sep","Oct","Nov","Dec"))


# 4. Line chart
ggplot(mean_long, aes(x = Bulan, y = Rata_rata, group = 1)) +
  geom_line() +
  geom_point() +
  labs(title = "Rata-rata Harga Bawang Merah per Bulan (2023)",
       x = "Bulan",
       y = "Rata-rata Harga") +
  theme_minimal()

Grafik menunjukkan bahwa harga bawang merah tahun 2023 berfluktuasi dengan dua puncak harga pada Februari dan Juni. Setelah itu harga menurun tajam dan mencapai titik terendah pada September–Oktober, lalu kembali naik di akhir tahun.

# 1. Filter hanya tahun 2022
data_2024 <- df_wide_final %>% 
  filter(Tahun == 2024)

# 2. Hitung rata-rata setiap bulan
mean_2024 <- data_2024 %>% 
  summarise(
    Jan = mean(Jan, na.rm = TRUE),
    Feb = mean(Feb, na.rm = TRUE),
    Mar = mean(Mar, na.rm = TRUE),
    Apr = mean(Apr, na.rm = TRUE),
    May = mean(May, na.rm = TRUE),
    Jun = mean(Jun, na.rm = TRUE),
    Jul = mean(Jul, na.rm = TRUE),
    Aug = mean(Aug, na.rm = TRUE),
    Sep = mean(Sep, na.rm = TRUE),
    Oct = mean(Oct, na.rm = TRUE),
    Nov = mean(Nov, na.rm = TRUE),
    Dec = mean(Dec, na.rm = TRUE)
  )

# 3. Ubah wide → long untuk plot
mean_long <- mean_2024 %>% 
  pivot_longer(cols = everything(),
               names_to = "Bulan",
               values_to = "Rata_rata")

# 4. Urutkan bulan
mean_long$Bulan <- factor(mean_long$Bulan,
                          levels = c("Jan","Feb","Mar","Apr","May",
                                     "Jun","Jul","Aug","Sep","Oct","Nov","Dec"))


# 4. Line chart
ggplot(mean_long, aes(x = Bulan, y = Rata_rata, group = 1)) +
  geom_line() +
  geom_point() +
  labs(title = "Rata-rata Harga Bawang Merah per Bulan (2023)",
       x = "Bulan",
       y = "Rata-rata Harga") +
  theme_minimal()

Harga bawang merah pada tahun 2024 meningkat pada bulan April-Mei. Kemudian menurun kembali di bulan Juli-Agustus dan meningkat kembali di akhir tahun.

4 Pemodelan Fuzzy C-Means

Fuzzy C-Means (FCM) adalah algoritma pengelompokan (clustering) yang memungkinkan suatu data memiliki tingkat keanggotaan pada beberapa cluster sekaligus. Tingkat keanggotaan ini menunjukkan seberapa kuat hubungan data tersebut dengan setiap cluster, dengan nilai antara 0 dan 1. Dalam FCM, setiap data tidak hanya dimiliki oleh satu cluster, melainkan dapat menjadi bagian dari beberapa cluster dengan nilai keanggotaan yang berbeda. Jika jumlah cluster yang ditentukan adalah \(c\), maka setiap data akan memiliki \(c\) nilai keanggotaan.

Algoritma Fuzzy C-Means membagi suatu himpunan \(X = \{X_1, X_2, \ldots, X_n\}\) menjadi \(c\) cluster berdasarkan kriteria tertentu. Setiap elemen \(X_i\) merupakan vektor berdimensi \(d\). Cluster diwakili oleh pusat-pusatnya \(V = \{V_1, V_2, \ldots, V_c\}\). Matriks \(U\) merepresentasikan tingkat keanggotaan setiap elemen dalam setiap cluster. Matriks \(U\) memiliki beberapa karakteristik berikut:

\(U(i,k)\) adalah nilai keanggotaan elemen \(X_i\) dalam cluster dengan pusat \(V_k\).
\(0 \leq U(i,k) \leq 1\) untuk semua \(i\) dan \(k\), dan jumlah nilai keanggotaan untuk setiap elemen di semua cluster adalah \(\sum_{k=1}^{c} U(i,k) = 1\).
Semakin besar \(U(i,k)\), semakin tinggi kemungkinan \(X_i\) termasuk dalam cluster \(k\).

Fungsi objektif \(J(U,V)\) didefinisikan sebagai:

\[ J(U,V) = \sum_{k=1}^{c} \sum_{i=1}^{n} U(i,k)^m D(i,k)^2 \]

dengan:
- \(D(i,k)^2 = \|X_i - V_k\|^2\) adalah kuadrat jarak antara elemen \(X_i\) dan pusat cluster \(V_k\).
- \(m\) adalah koefisien fuzzifikasi (\(m > 1\)).

Langkah-langkah algoritma FCM:

Inisialisasi
- Tentukan pusat cluster awal \(V\).
- Tetapkan \(t = 0\) (penghitung iterasi), \(\varepsilon > 0\) (kriteria penghentian), dan \(m > 1\) (koefisien fuzzifikasi).
Perbarui Matriks Keanggotaan
- Untuk setiap elemen \(X_i\) dan setiap cluster \(k\), hitung:
  \[ U(i,k) = \left( \sum_{j=1}^{c} \left( \frac{D(i,k)}{D(i,j)} \right)^{\frac{2}{m-1}} \right)^{-1} \]
Perbarui Pusat Cluster \(V\)
- Hitung pusat cluster baru dengan:
  \[ V_k = \frac{\sum_{i=1}^{n} U(i,k)^m X_i}{\sum_{i=1}^{n} U(i,k)^m} \]
Pemeriksaan Konvergensi
- Jika \(\|V^{(t+1)} - V^{(t)}\| < \varepsilon\), hentikan iterasi.
- Jika tidak, set \(t = t+1\) dan kembali ke langkah 2.
Akhir
- Algoritma selesai dan menghasilkan cluster-cluster yang telah terbentuk.

Algoritma ini cocok digunakan untuk data yang memiliki batasan antar cluster yang tidak jelas, seperti dalam analisis pola harga komoditas yang dapat mengalami fluktuasi dan tumpang tindih antar kategori (Sumargo dan Handhayani 2025).

4.1 Pemodelan Fuzzy C-Means Tahunan

4.1.1 Data

Data terlebih dahulu dihitung rata-rata harga tahunannya dari bulan Januari hingga Desember. Hasilnya kemudian diubah ke format lebar agar setiap tahun menjadi kolom tersendiri. Setelah dicek dan dibersihkan dari nilai hilang, data non-numerik dihapus, dan seluruh variabel numerik distandardisasi agar siap digunakan untuk analisis clusterin

df_tahunan <- df_wide_final %>%
  rowwise() %>% 
  mutate(Rata_Tahunan = mean(c_across(Jan:Dec), na.rm = TRUE)) %>%
  ungroup() %>%
  dplyr::select(KabKot, KodeProv, Tahun, Rata_Tahunan)

df_wide_tahunan <- df_tahunan %>%
  mutate(Tahun = paste0("Rata_Tahunan_", Tahun)) %>% 
  pivot_wider(
    names_from = Tahun,
    values_from = Rata_Tahunan
  )

colSums(is.na(df_wide_tahunan))

##            KabKot          KodeProv Rata_Tahunan_2022 Rata_Tahunan_2023 
##                 0                 0                 0                 0 
## Rata_Tahunan_2024 
##                 0

Terlihat tidak ada missing value pada data. Selanjutnya data distandarisasi.

# Pisahkan data numerik
data_numeric_yearly <- df_wide_tahunan %>%
  dplyr::select(-KabKot, -KodeProv)

# Handle missing values jika ada
data_numeric_yearly <- na.omit(data_numeric_yearly)

# Standardisasi data
data_scaled_yearly <- scale(data_numeric_yearly)

# Cek dimensi data
dim(data_scaled_yearly)

## [1] 506   3

4.1.2 Pemodelan

Selanjutnya dilakukan penentuan jumlah klaster optimal dengan mengevaluasi beberapa pilihan jumlah klaster menggunakan Silhouette Score dan Davies–Bouldin Index, sehingga dapat dipilih jumlah klaster yang paling sesuai sebelum menjalankan Fuzzy C-Means.

# Inisialisasi list yearly
fcm_results_yearly <- list()
cluster_centers_yearly <- list()

# Looping untuk centers 2 sampai 5
for(k in 2:5) {
  cat("Running Fuzzy C-Means YEARLY dengan", k, "clusters...\n")
  
  set.seed(123)
  fcm_yearly <- fcm(data_scaled_yearly, centers = k, m = 2, nstart = 10)
  
  # Simpan
  fcm_results_yearly[[paste0("k", k)]] <- fcm_yearly
  
  centers_df_y <- as.data.frame(fcm_yearly$v)
  colnames(centers_df_y) <- colnames(data_scaled_yearly)
  centers_df_y$Cluster <- paste0("Cluster_", 1:k)
  centers_df_y$K <- k
  
  cluster_centers_yearly[[paste0("k", k)]] <- centers_df_y
}

## Running Fuzzy C-Means YEARLY dengan 2 clusters...
## Running Fuzzy C-Means YEARLY dengan 3 clusters...
## Running Fuzzy C-Means YEARLY dengan 4 clusters...
## Running Fuzzy C-Means YEARLY dengan 5 clusters...

# Data frame validitas untuk yearly
validitas_yearly <- data.frame(
  K = 2:5,
  Silhouette_Score = numeric(4),
  Davies_Bouldin = numeric(4)
)

# Hitung validitas untuk setiap K
for(k in 2:5) {
  cat("\n📊 Menghitung validitas YEARLY untuk K =", k, "...\n")
  
  fcm_result_yearly <- fcm_results_yearly[[paste0("k", k)]]
  
  # Hard clusters dari membership
  hard_clusters_yearly <- apply(fcm_result_yearly$u, 1, which.max)
  
  # 1. Silhouette Score
  sil_score_yearly <- silhouette(hard_clusters_yearly, dist(data_scaled_yearly))
  avg_sil_width_yearly <- mean(sil_score_yearly[, 3])
  
  # 2. Davies-Bouldin Index
  db_index_yearly <- index.DB(data_scaled_yearly, hard_clusters_yearly)$DB
  
  # Simpan hasil
  validitas_yearly[validitas_yearly$K == k, "Silhouette_Score"] <- avg_sil_width_yearly
  validitas_yearly[validitas_yearly$K == k, "Davies_Bouldin"] <- db_index_yearly
  
  cat("✅ K =", k, 
      "| Silhouette:", round(avg_sil_width_yearly, 4),
      "| Davies-Bouldin:", round(db_index_yearly, 4), "\n")
}

## 
## 📊 Menghitung validitas YEARLY untuk K = 2 ...
## ✅ K = 2 | Silhouette: 0.7004 | Davies-Bouldin: 0.6355 
## 
## 📊 Menghitung validitas YEARLY untuk K = 3 ...
## ✅ K = 3 | Silhouette: 0.4668 | Davies-Bouldin: 0.7908 
## 
## 📊 Menghitung validitas YEARLY untuk K = 4 ...
## ✅ K = 4 | Silhouette: 0.4404 | Davies-Bouldin: 0.843 
## 
## 📊 Menghitung validitas YEARLY untuk K = 5 ...
## ✅ K = 5 | Silhouette: 0.3839 | Davies-Bouldin: 0.8754

print(validitas_yearly)

##   K Silhouette_Score Davies_Bouldin
## 1 2        0.7003975      0.6355358
## 2 3        0.4667828      0.7907938
## 3 4        0.4403714      0.8429639
## 4 5        0.3838934      0.8753999

# Plot Silhouette Score
p1_yearly <- ggplot(validitas_yearly, aes(x = K, y = Silhouette_Score)) +
  geom_line(size = 1, color = "steelblue") +
  geom_point(size = 3, color = "steelblue") +
  labs(
    title = "Silhouette Score vs Jumlah Cluster (YEARLY)",
    subtitle = "Semakin tinggi semakin baik",
    x = "Jumlah Cluster (K)",
    y = "Silhouette Score"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = 2:5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red")

## 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.

# Plot Davies-Bouldin Index
p2_yearly <- ggplot(validitas_yearly, aes(x = K, y = Davies_Bouldin)) +
  geom_line(size = 1, color = "coral") +
  geom_point(size = 3, color = "coral") +
  labs(
    title = "Davies-Bouldin Index vs Jumlah Cluster (YEARLY)",
    subtitle = "Semakin rendah semakin baik",
    x = "Jumlah Cluster (K)",
    y = "Davies-Bouldin Index"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = 2:5)

# Tampilkan
if(require(patchwork)) {
  library(patchwork)
  p1_yearly + p2_yearly
} else {
  print(p1_yearly)
  print(p2_yearly)
}

Berdasarkan hasil evaluasi, jumlah klaster terbaik adalah 2 klaster. Hal ini ditunjukkan oleh nilai Silhouette Score tertinggi (0.7004) serta Davies–Bouldin Index terendah (0.6355) dibandingkan jumlah klaster lainnya. Kedua indikator tersebut mengonfirmasi bahwa pemisahan klaster pada K = 2 menghasilkan struktur klaster yang paling jelas dan paling baik kualitasnya.

4.1.3 Visualisasi Hasil Cluster

fcm_result_k2_yearly <- fcm_results_yearly[["k2"]]  # KOREKSI: fcm_results_yearly
centers_k2_yearly <- cluster_centers_yearly[["k2"]]

centers_original_yearly <- centers_k2_yearly
for(i in 1:(ncol(centers_original_yearly)-2)) {
  centers_original_yearly[, i] <- centers_k2_yearly[, i] * 
    attr(data_scaled_yearly, 'scaled:scale')[i] + 
    attr(data_scaled_yearly, 'scaled:center')[i]
}

# 1. TRANSFORM DATA
# Identifikasi kolom yang berisi tahun (bukan Cluster dan K)
numeric_cols <- setdiff(colnames(centers_original_yearly), c("Cluster", "K"))

# Ekstrak tahun dari nama kolom
extracted_years <- as.numeric(gsub("Rata_Tahunan_", "", numeric_cols))

# 2. TRANSFORM KE LONG FORMAT DENGAN HANDLING NA
centers_long_yearly <- centers_original_yearly %>%
  pivot_longer(
    cols = all_of(numeric_cols),
    names_to = "Tahun_Label",
    values_to = "Harga_Rata_Rata"
  ) %>%
  mutate(
    # Coba ekstrak tahun dari label
    Tahun = as.numeric(gsub("[^0-9]", "", Tahun_Label)),  # Ekstrak angka saja
    Tahun = ifelse(is.na(Tahun), Tahun_Label, Tahun),     # Jika gagal, pakai label
    Cluster = as.factor(Cluster)
  )

4.1.3.1 Profil

# Konversi Tahun ke numeric jika masih character
  centers_long_yearly$Tahun <- as.numeric(centers_long_yearly$Tahun)
  
  p_trend_yearly <- ggplot(centers_long_yearly, 
                           aes(x = Tahun, y = Harga_Rata_Rata, 
                               color = Cluster, group = Cluster)) +
    geom_line(size = 2, alpha = 0.8) +
    geom_point(size = 4, alpha = 0.9) +
    geom_text(aes(label = format(round(Harga_Rata_Rata, 0), big.mark = ",")), 
              size = 3, vjust = -0.8, show.legend = FALSE) +
    labs(
      title = "Trend Rata-rata Harga Bawang per Cluster (Tahunan)",
      subtitle = "Perkembangan harga dari tahun ke tahun untuk setiap cluster",
      x = "Tahun",
      y = "Harga Rata-rata (Rp)",
      color = "Cluster"
    ) +
    theme_minimal() +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks = unique(centers_long_yearly$Tahun)) +
    scale_color_manual(values = c("#E74C3C", "#3498DB")) +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      plot.subtitle = element_text(size = 10)
    )
  
  print(p_trend_yearly)

Grafik menunjukkan bahwa kedua cluster memiliki pola harga yang sama dari 2022–2024, yaitu penurunan pada 2023 lalu kembali naik pada 2024. Cluster 1 selalu berada pada kisaran harga lebih tinggi, turun dari Rp53.151 (2022) ke Rp50.137 (2023), lalu naik menjadi Rp54.287 (2024). Cluster 2 juga mengalami pola yang sama dengan level harga lebih rendah, yakni Rp35.877 (2022), turun ke Rp31.435 (2023), dan naik ke Rp34.651 (2024). Pola ini menunjukkan adanya titik harga terendah pada 2023 sebelum terjadi pemulihan di 2024.

4.1.3.2 Distribusi Probabilitas

# HITUNG MEMBERSHIP DAN HARD CLUSTERS
membership_k2_yearly <- as.data.frame(fcm_result_k2_yearly$u)
colnames(membership_k2_yearly) <- c("Cluster_1", "Cluster_2")
hard_clusters_yearly <- apply(membership_k2_yearly, 1, which.max)

# TAMBAHKAN CLUSTER KE DATA ASLI
data_with_clusters_yearly <- df_wide_tahunan
data_with_clusters_yearly$Cluster <- paste0("Cluster_", hard_clusters_yearly)
data_with_clusters_yearly$Prob_Max <- apply(membership_k2_yearly, 1, max)

membership_long_yearly <- membership_k2_yearly %>%
  mutate(KabKot = data_with_clusters_yearly$KabKot) %>%
  pivot_longer(cols = c(Cluster_1, Cluster_2), 
               names_to = "Cluster", values_to = "Probability")

p_membership_yearly <- ggplot(membership_long_yearly, aes(x = Probability, fill = Cluster)) +
  geom_histogram(alpha = 0.7, bins = 20, position = "identity") +
  labs(
    title = "Distribusi Probabilitas Keanggotaan Cluster",
    subtitle = "Semakin dekat ke 1, semakin pasti keanggotaannya",
    x = "Probabilitas Keanggotaan",
    y = "Frekuensi"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("#E74C3C", "#3498DB")) +
  facet_wrap(~Cluster, ncol = 2)

print(p_membership_yearly)

Histogram menunjukkan bahwa probabilitas keanggotaan kedua cluster sangat terpisah dengan jelas. Cluster 1 didominasi oleh observasi dengan probabilitas yang sangat rendah (mendekati 0), sedangkan Cluster 2 didominasi oleh probabilitas yang sangat tinggi (mendekati 1). Pola ini menandakan bahwa proses clustering berhasil memisahkan data dengan baik, karena sebagian besar observasi memiliki tingkat kepastian keanggotaan yang tinggi terhadap cluster masing-masing. Hanya sedikit titik yang memiliki probabilitas di area tengah, sehingga ambiguitas dalam penentuan cluster hampir tidak terjadi.

4.1.3.3 Visualisasi PCA

# 3. PCA VISUALIZATION 2D
pca_result_yearly <- prcomp(data_scaled_yearly, scale. = TRUE)
pca_df_yearly <- as.data.frame(pca_result_yearly$x)
pca_df_yearly$Cluster <- data_with_clusters_yearly$Cluster
pca_df_yearly$KabKot <- data_with_clusters_yearly$KabKot

# Hitung variance explained
variance_pc1 <- round(summary(pca_result_yearly)$importance[2,1] * 100, 1)
variance_pc2 <- round(summary(pca_result_yearly)$importance[2,2] * 100, 1)

# 4. PLOT PCA
p_pca_yearly <- ggplot(pca_df_yearly, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(alpha = 0.7, size = 2) +
  stat_ellipse(level = 0.7, size = 1) +
  labs(
    title = "Visualisasi Cluster dalam 2D PCA - Data Tahunan",
    subtitle = "Setiap titik mewakili satu kabupaten/kota",
    x = paste("PC1 (", variance_pc1, "%)"),
    y = paste("PC2 (", variance_pc2, "%)"),
    color = "Cluster"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("#E74C3C", "#3498DB")) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10),
    legend.position = "right"
  )

print(p_pca_yearly)

Plot PCA dua dimensi menunjukkan bahwa kedua cluster terpisah dengan jelas pada ruang PC1–PC2. Cluster 1 terkonsentrasi pada nilai PC1 yang lebih rendah, sedangkan Cluster 2 berada pada rentang PC1 yang lebih tinggi. Pemisahan yang tegas di sepanjang PC1, yang menjelaskan 94,3% variasi data, mengindikasikan bahwa dimensi utama ini sangat menentukan perbedaan karakteristik antar cluster. Sebaran titik dalam masing-masing cluster juga terlihat relatif kompak, menegaskan bahwa struktur clustering stabil dan mampu membedakan kelompok kabupaten/kota dengan baik. Plot ini mendukung hasil probabilitas keanggotaan, karena tidak banyak area tumpang-tindih antar cluster.

4.1.3.4 Distribusi Kab/Kota

# Plot distribusi cluster
p_distribusi <- ggplot(data_with_clusters_yearly, aes(x = Cluster, fill = Cluster)) +
  geom_bar(alpha = 0.8) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, fontface = "bold") +
  labs(
    title = "Distribusi Kabupaten/Kota per Cluster",
    subtitle = paste("Total", nrow(data_with_clusters_yearly), "kabupaten/kota"),
    x = "Cluster",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("#E74C3C", "#3498DB"))

print(p_distribusi)

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Berdasarkan barchart, terlihat bahwa Cluster 1 beranggotakan 68 kabupaten/kota, sedangkan Cluster 2 memiliki anggota jauh lebih banyak, yaitu 438 kabupaten/kota. Perbedaan jumlah anggota ini menunjukkan bahwa sebagian besar wilayah cenderung memiliki karakteristik yang lebih mirip dengan Cluster 2, sementara Cluster 1 merupakan kelompok yang lebih kecil dan lebih spesifik.

4.1.3.5 Choropleth Map

4.1.3.5.1 Merge Data

data_with_clusters_yearly$KabKot_clean <- data_with_clusters_yearly$KabKot |>
  str_replace("^Kab\\.\\s*", "") |>   # hapus "Kab. "
  str_replace("^Kota\\s*", "")        # hapus "Kota "

indonesia_shp$ADM2_EN_baru <- indonesia_shp$ADM2_EN |>
  str_replace("^Kab\\.\\s*", "") |>   # hapus "Kab. "
  str_replace("^Kota\\s*", "")        # hapus "Kota "

# Data Anda punya tapi shapefile tidak:
missing_in_shp <- setdiff(data_with_clusters_yearly$KabKot_clean, indonesia_shp$ADM2_EN_baru)
cat("Kabupaten di data Anda tapi tidak di shapefile (", length(missing_in_shp), "):\n")

## Kabupaten di data Anda tapi tidak di shapefile ( 9 ):

missing_in_shp

## [1] "Karangasem"         "Kepulauan Tanimbar" "Mahakam Ulu"       
## [4] "Pasangkayu"         "Toba"               "Tojo Una-una"      
## [7] "Toli-toli"          "Lubuk Linggau"      "Palangkaraya"

# Shapefile punya tapi data Anda tidak:
missing_in_data <- setdiff(indonesia_shp$ADM2_EN_baru, data_with_clusters_yearly$KabKot_clean)
cat("\nKabupaten di shapefile tapi tidak di data Anda (", length(missing_in_data), "):\n")

## 
## Kabupaten di shapefile tapi tidak di data Anda ( 22 ):

missing_in_data

##  [1] "Asmat"                 "Danau"                 "Danau Toba"           
##  [4] "Deiyai"                "Hutan"                 "Karang Asem"          
##  [7] "Lubuklinggau"          "Palangka Raya"         "Lanny Jaya"           
## [10] "Mahakam Hulu"          "Maluku Tenggara Barat" "Mamberamo Raya"       
## [13] "Mamberamo Tengah"      "Paniai"                "Pegunungan Arfak"     
## [16] "Supiori"               "Toba Samosir"          "Tojo Una-Una"         
## [19] "Toli-Toli"             "Tolikara"              "Waduk Cirata"         
## [22] "Wadung Kedungombo"

nama_fix <- data.frame(
  wrong = c(
    "Karangasem",
    "Kepulauan Tanimbar",
    "Mahakam Ulu",
    "Pasangkayu",
    "Toba",
    "Tojo Una-una",
    "Toli-toli",
    "Lubuk Linggau",
    "Palangkaraya"
  ),
  right = c(
    "Karang Asem",
    "Maluku Tenggara Barat",
    "Mahakam Hulu",
    "Pasangkayu",          # shapefile memang tidak punya
    "Toba Samosir",
    "Tojo Una-Una",
    "Toli-Toli",
    "Lubuklinggau",
    "Palangka Raya"
  ),
  stringsAsFactors = FALSE
)

for (i in seq_len(nrow(nama_fix))) {
  data_with_clusters_yearly$KabKot_clean[
    data_with_clusters_yearly$KabKot_clean == nama_fix$wrong[i]
  ] <- nama_fix$right[i]
}

missing_in_shp <- setdiff(data_with_clusters_yearly$KabKot_clean,
                          indonesia_shp$ADM2_EN_baru)

missing_in_data <- setdiff(indonesia_shp$ADM2_EN_baru,
                           data_with_clusters_yearly$KabKot_clean)

cat("Masih mismatch (data → shp):\n")

## Masih mismatch (data → shp):

missing_in_shp

## [1] "Pasangkayu"

cat("\nMasih mismatch (shp → data):\n")

## 
## Masih mismatch (shp → data):

missing_in_data

##  [1] "Asmat"             "Danau"             "Danau Toba"       
##  [4] "Deiyai"            "Hutan"             "Lanny Jaya"       
##  [7] "Mamberamo Raya"    "Mamberamo Tengah"  "Paniai"           
## [10] "Pegunungan Arfak"  "Supiori"           "Tolikara"         
## [13] "Waduk Cirata"      "Wadung Kedungombo"

shp_clusteryearly <- indonesia_shp %>%
  left_join(data_with_clusters_yearly, by = c("ADM2_EN_baru" = "KabKot_clean"))

## Warning in sf_column %in% names(g): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 18 of `x` matches multiple rows in `y`.
## ℹ Row 17 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

4.1.3.5.2 Peta Klaster Harga Bawang Merah

# Pastikan Cluster jadi factor tanpa NA levels
shp_clusteryearly$Cluster <- factor(shp_clusteryearly$Cluster)

ggplot(shp_clusteryearly) +
  geom_sf(aes(fill = Cluster), color = "black", size = 0.001) +  # stroke tipis
  scale_fill_manual(
    values = c(
      "Cluster_1" = "#F8766D",
      "Cluster_2" = "#00BFC4"
    ),
    drop = TRUE,        # Hilangkan NA dari legend
    na.value = NA,
    name = "Cluster"
  ) +
  labs(title = "Peta Klaster Fuzzy C-Means (2 Klaster)") +
  theme_minimal()

Hasil pemetaan menunjukkan bahwa wilayah timur Indonesia cenderung memiliki harga bawang merah yang lebih tinggi, sejalan dengan dominasi Cluster 2 di kawasan tersebut. Sebaliknya, daerah di bagian barat, yang lebih banyak masuk ke dalam Cluster 1, memiliki tingkat harga yang relatif lebih rendah. Pola ini mengindikasikan adanya perbedaan regional dalam struktur harga, yang kemungkinan dipengaruhi oleh faktor distribusi, akses pasokan, serta kondisi geografis antar wilayah.

4.2 Pemodelan Fuzzy C-Means Bulanan (Tahun 2024)

Pada tahap selanjutnya dilakukan pemodelan Fuzzy C-Means berdasarkan data bulanan tahun 2024 untuk melihat pola pergerakan harga bawang secara lebih rinci setiap bulan.

4.2.1 Data

Data difilter menjadi tahun 2024 saja, lalu distandardisasi.

tahun2024 <- df_wide_final %>%
  filter(Tahun == 2024)

head(tahun2024)

## # A tibble: 6 × 15
##   KabKot    KodeProv Tahun   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep
##   <chr>        <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kab. Ace…       11 2024  38386 37140 38623 55768 56542 50893 38508 28000 29200
## 2 Kab. Ace…       11 2024  30419 30828 33484 50630 54419 50148 33067 24700 22357
## 3 Kab. Ace…       11 2024  35097 32345 37031 53077 58710 45750 32700 25484 26467
## 4 Kab. Ace…       11 2024  40290 38345 39250 56786 59839 51407 39032 29032 29467
## 5 Kab. Ace…       11 2024  41858 38103 37156 52661 50871 44750 35419 27129 26233
## 6 Kab. Ace…       11 2024  38857 33929 40000 51154 55484 47885 30185 24310 26556
## # ℹ 3 more variables: Oct <dbl>, Nov <dbl>, Dec <dbl>

tahun2024 <- tahun2024 %>%
  rename_with(
    .cols = Jan:Dec,
    .fn = ~ paste0(., "_", tahun2024$Tahun[1])
  ) %>%
  dplyr::select(-Tahun)

head(tahun2024)

## # A tibble: 6 × 14
##   KabKot KodeProv Jan_2024 Feb_2024 Mar_2024 Apr_2024 May_2024 Jun_2024 Jul_2024
##   <chr>     <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Kab. …       11    38386    37140    38623    55768    56542    50893    38508
## 2 Kab. …       11    30419    30828    33484    50630    54419    50148    33067
## 3 Kab. …       11    35097    32345    37031    53077    58710    45750    32700
## 4 Kab. …       11    40290    38345    39250    56786    59839    51407    39032
## 5 Kab. …       11    41858    38103    37156    52661    50871    44750    35419
## 6 Kab. …       11    38857    33929    40000    51154    55484    47885    30185
## # ℹ 5 more variables: Aug_2024 <dbl>, Sep_2024 <dbl>, Oct_2024 <dbl>,
## #   Nov_2024 <dbl>, Dec_2024 <dbl>

# Pisahkan data numerik
data_numeric_2024 <- tahun2024 %>%
  dplyr::select(-KabKot, -KodeProv)

# Handle missing values jika ada
data_numeric_2024 <- na.omit(data_numeric_2024)

# Standardisasi data
data_scaled_2024 <- scale(data_numeric_2024)

# Cek dimensi data
dim(data_scaled_2024)

## [1] 506  12

4.2.2 Pemodelan

Selanjutnya dilakukan penentuan jumlah klaster optimal dengan mengevaluasi beberapa pilihan jumlah klaster menggunakan Silhouette Score dan Davies–Bouldin Index, sehingga dapat dipilih konfigurasi klaster yang paling sesuai sebelum menjalankan Fuzzy C-Means.

# Inisialisasi list untuk menyimpan hasil
fcm_results <- list()
cluster_centers_list <- list()

# Looping untuk centers 2 sampai 5
for(k in 2:5) {
  cat("Running Fuzzy C-Means dengan", k, "clusters...\n")
  
  set.seed(123) # untuk reproduktibilitas
  fcm_result <- fcm(data_scaled_2024, centers = k, m = 2, nstart = 10)
  
  # Simpan hasil
  fcm_results[[paste0("k", k)]] <- fcm_result
  
  # Simpan cluster centers
  centers_df <- as.data.frame(fcm_result$v)
  colnames(centers_df) <- colnames(data_scaled_2024)
  centers_df$Cluster <- paste0("Cluster_", 1:k)
  centers_df$K <- k
  
  cluster_centers_list[[paste0("k", k)]] <- centers_df
}

## Running Fuzzy C-Means dengan 2 clusters...
## Running Fuzzy C-Means dengan 3 clusters...
## Running Fuzzy C-Means dengan 4 clusters...
## Running Fuzzy C-Means dengan 5 clusters...

# Data frame untuk menyimpan hasil validitas
validitas_df <- data.frame(
  K = 2:5,
  Silhouette_Score = numeric(4),
  Davies_Bouldin = numeric(4)
)

# Gunakan package clusterSim untuk Davies-Bouldin
if(!require(clusterSim)) install.packages("clusterSim")
library(clusterSim)

# Hitung validitas untuk setiap K
for(k in 2:5) {
  cat("\n📊 Menghitung validitas untuk K =", k, "...\n")
  
  fcm_result <- fcm_results[[paste0("k", k)]]
  
  # Hitung hard clusters dari membership
  hard_clusters <- apply(fcm_result$u, 1, which.max)
  
  # 1. Silhouette Score
  sil_score <- silhouette(hard_clusters, dist(data_scaled_2024))
  avg_sil_width <- mean(sil_score[, 3])
  
  # 2. Davies-Bouldin Index dari package clusterSim
  db_index <- index.DB(data_scaled_2024, hard_clusters)$DB
  
  # Simpan hasil
  validitas_df[validitas_df$K == k, "Silhouette_Score"] <- avg_sil_width
  validitas_df[validitas_df$K == k, "Davies_Bouldin"] <- db_index
  
  cat("✅ K =", k, 
      "| Silhouette:", round(avg_sil_width, 4),
      "| Davies-Bouldin:", round(db_index, 4), "\n")
}

## 
## 📊 Menghitung validitas untuk K = 2 ...
## ✅ K = 2 | Silhouette: 0.6019 | Davies-Bouldin: 0.867 
## 
## 📊 Menghitung validitas untuk K = 3 ...
## ✅ K = 3 | Silhouette: 0.3186 | Davies-Bouldin: 1.1122 
## 
## 📊 Menghitung validitas untuk K = 4 ...
## ✅ K = 4 | Silhouette: 0.2359 | Davies-Bouldin: 1.2953 
## 
## 📊 Menghitung validitas untuk K = 5 ...
## ✅ K = 5 | Silhouette: 0.2211 | Davies-Bouldin: 1.2017

print(validitas_df)

##   K Silhouette_Score Davies_Bouldin
## 1 2        0.6019215      0.8670243
## 2 3        0.3186115      1.1121670
## 3 4        0.2359387      1.2952580
## 4 5        0.2210804      1.2017076

# Plot Silhouette Score
p1 <- ggplot(validitas_df, aes(x = K, y = Silhouette_Score)) +
  geom_line(size = 1, color = "steelblue") +
  geom_point(size = 3, color = "steelblue") +
  labs(
    title = "Silhouette Score vs Jumlah Cluster",
    subtitle = "Semakin tinggi semakin baik",
    x = "Jumlah Cluster (K)",
    y = "Silhouette Score"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = 2:5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red")

# Plot Davies-Bouldin Index
p2 <- ggplot(validitas_df, aes(x = K, y = Davies_Bouldin)) +
  geom_line(size = 1, color = "coral") +
  geom_point(size = 3, color = "coral") +
  labs(
    title = "Davies-Bouldin Index vs Jumlah Cluster",
    subtitle = "Semakin rendah semakin baik",
    x = "Jumlah Cluster (K)",
    y = "Davies-Bouldin Index"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = 2:5)

# Tampilkan kedua plot
if(require(patchwork)) {
  library(patchwork)
  p1 + p2
} else {
  print(p1)
  print(p2)
}

Berdasarkan hasil evaluasi, jumlah klaster terbaik adalah 2 klaster. Hal ini ditunjukkan oleh nilai Silhouette Score tertinggi (0.6019215) serta Davies–Bouldin Index terendah (0.8670243) dibandingkan jumlah klaster lainnya. Kedua indikator tersebut mengonfirmasi bahwa pemisahan klaster pada K = 2 menghasilkan struktur klaster yang paling jelas dan paling baik kualitasnya.

4.2.3 Visualisasi Hasil Cluster

4.2.3.1 Profil

# Ambil hasil untuk K=2
fcm_result_k2 <- fcm_results[["k2"]]
centers_k2 <- cluster_centers_list[["k2"]]

# 1. TRANSFORM CLUSTER CENTERS KE SKALA ASLI
centers_original <- centers_k2
for(i in 1:(ncol(centers_original)-2)) { # -2 untuk kolom Cluster dan K
  centers_original[, i] <- centers_k2[, i] * attr(data_scaled_2024, 'scaled:scale')[i] + 
                          attr(data_scaled_2024, 'scaled:center')[i]
}

# 2. VISUALISASI PROFIL CLUSTER - TREND BULANAN
centers_long <- centers_original %>%
  pivot_longer(
    cols = -c(Cluster, K),
    names_to = "Bulan",
    values_to = "Harga"
  ) %>%
  mutate(
    Bulan = str_remove(Bulan, "_2024"),
    Bulan = factor(Bulan, levels = month.abb),
    Tanggal = as.Date(paste("2024", match(Bulan, month.abb), "15", sep = "-"))
  )

# Plot trend cluster centers
p_trend <- ggplot(centers_long, aes(x = Tanggal, y = Harga, color = Cluster, group = Cluster)) +
  geom_line(size = 2) +
  geom_point(size = 3) +
  labs(
    title = "Profil Harga Bawang per Cluster - 2024",
    subtitle = "Cluster centers menunjukkan pola harga typical",
    x = "Bulan",
    y = "Harga (Rp)",
    color = "Cluster"
  ) +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_color_manual(values = c("#E74C3C", "#3498DB"))

print(p_trend)

Plot profil harga bulanan tahun 2024 menunjukkan bahwa Cluster 1 konsisten memiliki harga yang lebih tinggi sepanjang tahun. Puncak harga pada cluster ini terjadi pada bulan Juni dengan kecenderungan naik sejak April sebelum akhirnya turun kembali setelah Juli. Sementara itu, Cluster 2 memiliki harga yang lebih rendah dan menunjukkan pola kenaikan tajam pada April hingga Juni, kemudian penurunan signifikan pada Juli dan Agustus. Kedua cluster mengalami pola naik lagi menjelang akhir tahun, walaupun level harganya tetap berbeda cukup jauh. Pola ini menggambarkan bahwa perbedaan harga antar cluster tidak hanya muncul pada level tahunan tetapi juga pada dinamika bulanan.

4.2.3.2 Distribusi Probabilitas

# Hitung hard clusters
membership_k2 <- as.data.frame(fcm_result_k2$u)
colnames(membership_k2) <- c("Cluster_1", "Cluster_2")
hard_clusters <- apply(membership_k2, 1, which.max)

# Tambahkan ke data asli (pastikan data_2024_clean sudah ada)
data_with_clusters <- tahun2024
data_with_clusters$Cluster <- paste0("Cluster_", hard_clusters)
data_with_clusters$Prob_Max <- apply(membership_k2, 1, max)

# VISUALISASI MEMBERSHIP PROBABILITIES
membership_long <- membership_k2 %>%
  mutate(KabKot = data_with_clusters$KabKot) %>%
  pivot_longer(cols = c(Cluster_1, Cluster_2), 
               names_to = "Cluster", values_to = "Probability")

p_membership <- ggplot(membership_long, aes(x = Probability, fill = Cluster)) +
  geom_histogram(alpha = 0.7, bins = 20, position = "identity") +
  labs(
    title = "Distribusi Probabilitas Keanggotaan Cluster",
    subtitle = "Semakin dekat ke 1, semakin pasti keanggotaannya",
    x = "Probabilitas Keanggotaan",
    y = "Frekuensi"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("#E74C3C", "#3498DB")) +
  facet_wrap(~Cluster, ncol = 2)

print(p_membership)

4.2.3.3 Visualisasi PCA

# PCA VISUALIZATION 2D
pca_result <- prcomp(data_scaled_2024, scale. = TRUE)
pca_df <- as.data.frame(pca_result$x)
pca_df$Cluster <- data_with_clusters$Cluster
pca_df$KabKot <- data_with_clusters$KabKot

p_pca <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(alpha = 0.7, size = 2) +
  stat_ellipse(level = 0.7, size = 1) +
  labs(
    title = "Visualisasi Cluster dalam 2D PCA",
    subtitle = "Setiap titik mewakili satu kabupaten/kota",
    x = paste("PC1 (", round(summary(pca_result)$importance[2,1]*100, 1), "%)"),
    y = paste("PC2 (", round(summary(pca_result)$importance[2,2]*100, 1), "%)")
  ) +
  theme_minimal() +
  scale_color_manual(values = c("#E74C3C", "#3498DB"))

print(p_pca)

Plot PCA dua dimensi memperlihatkan bahwa kedua cluster terpisah dengan cukup tegas pada bidang PC1 dan PC2. Cluster 1 mengelompok pada bagian dengan nilai PC1 yang lebih rendah, sementara Cluster 2 berada pada rentang PC1 yang lebih tinggi. Pemisahan yang jelas sepanjang PC1, yang menyumbang porsi terbesar variasi data, menunjukkan bahwa komponen utama ini menjadi faktor paling berpengaruh dalam membedakan karakteristik antar cluster. Pola sebaran titik dalam tiap cluster tampak cukup rapat dan konsisten sehingga menguatkan bahwa hasil pengelompokan stabil dan mampu memisahkan kabupaten atau kota dengan baik. Visualisasi ini juga sejalan dengan pola probabilitas keanggotaan yang minim tumpang tindih antara kedua cluster.

4.2.3.4 Distribusi Kab/Kota

# VISUALISASI DISTRIBUSI ANGGOTA CLUSTER
# Plot distribusi cluster
p_distribusi <- ggplot(data_with_clusters, aes(x = Cluster, fill = Cluster)) +
  geom_bar(alpha = 0.8) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, fontface = "bold") +
  labs(
    title = "Distribusi Kabupaten/Kota per Cluster",
    subtitle = paste("Total", nrow(data_with_clusters), "kabupaten/kota"),
    x = "Cluster",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("#E74C3C", "#3498DB"))

print(p_distribusi)

Berdasarkan barchart, terlihat bahwa Cluster 1 beranggotakan 70 kabupaten/kota, sedangkan Cluster 2 memiliki anggota jauh lebih banyak, yaitu 436 kabupaten/kota. Perbedaan jumlah anggota ini menunjukkan bahwa sebagian besar wilayah cenderung memiliki karakteristik yang lebih mirip dengan Cluster 2, sementara Cluster 1 merupakan kelompok yang lebih kecil dan lebih spesifik.

4.2.3.5 Choropleth Map

4.2.3.5.1 Merge Data

data_with_clusters$KabKot_clean <- data_with_clusters$KabKot |>
  str_replace("^Kab\\.\\s*", "") |>   # hapus "Kab. "
  str_replace("^Kota\\s*", "")        # hapus "Kota "

data_with_clusters

## # A tibble: 506 × 17
##    KabKot         KodeProv Jan_2024 Feb_2024 Mar_2024 Apr_2024 May_2024 Jun_2024
##    <chr>             <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
##  1 Kab. Aceh Bar…       11    38386    37140    38623    55768    56542    50893
##  2 Kab. Aceh Bar…       11    30419    30828    33484    50630    54419    50148
##  3 Kab. Aceh Bes…       11    35097    32345    37031    53077    58710    45750
##  4 Kab. Aceh Jaya       11    40290    38345    39250    56786    59839    51407
##  5 Kab. Aceh Sel…       11    41858    38103    37156    52661    50871    44750
##  6 Kab. Aceh Sin…       11    38857    33929    40000    51154    55484    47885
##  7 Kab. Aceh Tam…       11    37452    36259    37952    54040    54000    47275
##  8 Kab. Aceh Ten…       11    38022    33696    36841    49259    56393    49444
##  9 Kab. Aceh Ten…       11    37571    28321    30474    41377    47307    46754
## 10 Kab. Aceh Tim…       11    36419    31964    34906    46107    53065    47714
## # ℹ 496 more rows
## # ℹ 9 more variables: Jul_2024 <dbl>, Aug_2024 <dbl>, Sep_2024 <dbl>,
## #   Oct_2024 <dbl>, Nov_2024 <dbl>, Dec_2024 <dbl>, Cluster <chr>,
## #   Prob_Max <dbl>, KabKot_clean <chr>

shp_cluster2024 <- indonesia_shp %>%
  left_join(data_with_clusters, by = c("ADM2_EN_baru" = "KabKot_clean"))

## Warning in sf_column %in% names(g): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 18 of `x` matches multiple rows in `y`.
## ℹ Row 17 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

4.2.3.5.2 Peta Klaster Harga Bawang

# Pastikan Cluster jadi factor tanpa NA levels
shp_cluster2024$Cluster <- factor(shp_cluster2024$Cluster)

ggplot(shp_cluster2024) +
  geom_sf(aes(fill = Cluster), color = "black", size = 0.001) +  # stroke tipis
  scale_fill_manual(
    values = c(
      "Cluster_1" = "#F8766D",
      "Cluster_2" = "#00BFC4"
    ),
    drop = TRUE,        # Hilangkan NA dari legend
    na.value = NA,
    name = "Cluster"
  ) +
  labs(title = "Peta Klaster Fuzzy C-Means (2 Klaster)") +
  theme_minimal()

Daftar Pustaka:

Sumargo EG, Handhayani T. 2025. ANALISIS CLUSTERING HARGA PANGAN DI PASAR TRADISIONAL WILAYAH KALIMANTAN MENGGUNAKAN ALGORITMA K-MEANS DAN FUZZY C-MEANS. Jurnal Mahasiswa Teknik Informatika. 9(5): 8703-8710

Analisis Pola Harga Bawang Merah di Indonesia (2022–2024) dengan Metode Fuzzy C-Means

Aliyah Nadhratunnaim NIM G1401231087

2025-12-01