Implementasi Metode Clustering pada Data Tingkat Kemiskinan Kabupaten/Kota di Indonesia

1. Pendahuluan

1.1 Latar Belakang

Kemiskinan merupakan permasalahan multidimensi yang dipengaruhi oleh berbagai faktor sosial, ekonomi, dan demografis. Pemahaman yang mendalam terhadap pola distribusi kemiskinan antarwilayah di Indonesia menjadi krusial untuk mendukung perencanaan pembangunan yang tepat sasaran. Analisis clustering memungkinkan pengelompokan kabupaten/kota berdasarkan kemiripan profil indikator kemiskinan, sehingga kebijakan dapat dirancang secara lebih spesifik sesuai karakteristik tiap klaster.

1.2 Deskripsi Dataset

Dataset ini merupakan kumpulan indikator sosial ekonomi seluruh kabupaten/kota di Indonesia tahun 2024 yang bersumber dari Badan Pusat Statistik (BPS). Data mencakup berbagai dimensi kemiskinan yang meliputi aspek pendidikan, kesehatan, ekonomi, dan pola konsumsi.

Variabel yang tersedia beserta penjelasannya:

Rata-rata Lama Sekolah (RLS): Jumlah tahun pendidikan yang telah ditempuh penduduk usia 15 tahun ke atas
Indeks Pembangunan Gender (IPG): Tingkat kesetaraan pembangunan antara laki-laki dan perempuan
Usia Harapan Hidup (UHH): Kondisi kesehatan dan harapan hidup penduduk di suatu wilayah
Pengeluaran Per Kapita: Kemampuan ekonomi rumah tangga dalam memenuhi kebutuhan
Harapan Lama Sekolah (HLS): Potensi lama pendidikan yang akan ditempuh oleh anak usia sekolah
Produk Domestik Regional Bruto (PDRB): Aktivitas ekonomi dan produktivitas wilayah
Indeks Kemahalan Konstruksi (IKK): Tingkat kemahalan konstruksi yang terkait dengan kondisi sosial ekonomi daerah
Pengeluaran Rokok Per Kapita: Pola konsumsi yang kurang produktif pada rumah tangga
Kabupaten Kota: Satuan wilayah administrasi tingkat II yang menjadi unit analisis dalam data, mencakup daerah kabupaten dan kota
Tingkat Penduduk Miskin: Persentase penduduk yang berada di bawah garis kemiskinan di suatu wilayah, yang mencerminkan tingkat kesejahteraan masyarakat

1.3 Tujuan Analisis

Mengelompokkan kabupaten/kota di Indonesia ke dalam klaster berdasarkan kesamaan profil kemiskinan menggunakan lima metode clustering yang berbeda.
Membandingkan performa masing-masing metode menggunakan indeks validitas internal (Silhouette, Dunn, Calinski-Harabasz).
Mengidentifikasi metode terbaik dan menginterpretasikan profil tiap klaster untuk mendukung rekomendasi kebijakan.

2. Persiapan

2.1 Memuat Package

suppressPackageStartupMessages({
  library(readxl)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(cluster)
  library(factoextra)
  library(dbscan)
  library(e1071)
  library(clusterCrit)
  library(gridExtra)
  library(RColorBrewer)
  library(fpc)
  library(knitr)
  library(kableExtra)
})

2.2 Memuat Data

file_path <- "Data_Tingkat_Kemiskinan.xlsx"
if (!file.exists(file_path)) stop("File tidak ditemukan di working directory!")

df_raw   <- read_excel(file_path)
kabkota  <- df_raw$Kabupaten_Kota
num_cols <- names(df_raw)[names(df_raw) != "Kabupaten_Kota"]
df_num   <- as.data.frame(lapply(df_raw[, num_cols], as.numeric))

Dataset berhasil dimuat dengan 514 kabupaten/kota dan 9 variabel numerik.

3. Eksplorasi Data

3.1 Statistika Deskriptif

desc_stats <- data.frame(
  Variabel = num_cols,
  Min      = round(sapply(df_num, min,              na.rm = TRUE), 2),
  Q1       = round(sapply(df_num, quantile, 0.25,   na.rm = TRUE), 2),
  Median   = round(sapply(df_num, median,            na.rm = TRUE), 2),
  Mean     = round(sapply(df_num, mean,              na.rm = TRUE), 2),
  Q3       = round(sapply(df_num, quantile, 0.75,   na.rm = TRUE), 2),
  Max      = round(sapply(df_num, max,               na.rm = TRUE), 2),
  Skewness = round(sapply(df_num, function(x) {
    mean((x - mean(x, na.rm=TRUE))^3) / sd(x, na.rm=TRUE)^3
  }), 3)
)

kable(desc_stats, row.names = FALSE,
      caption = "Tabel 1. Statistika Deskriptif Variabel Indikator Kemiskinan") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE, font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 1. Statistika Deskriptif Variabel Indikator Kemiskinan
Variabel	Min	Q1	Median	Mean	Q3	Max	Skewness
Tingkat_Penduduk_Miskin	2.23	6.38	9.32	11.17	13.77	41.42	1.637
RLS	1.92	7.82	8.61	8.74	9.61	13.10	-0.382
IPG	58.04	87.80	91.43	90.68	94.78	99.56	-1.507
UHH	55.74	68.26	70.59	70.41	72.78	78.26	-0.506
PengeluaranPerKapita	4597.00	9648.00	11305.00	11435.44	12896.75	25573.00	0.687
HLS	4.45	12.60	13.13	13.20	13.74	17.94	-1.209
PDRB	304.95	7322.40	17918.20	42958.83	39678.26	922861.10	5.906
IKK	77.03	93.84	99.66	105.52	106.62	379.81	5.687
PengeluaranPerkapita_Rokok	1190.00	17492.75	21883.50	22062.99	26127.25	47983.00	0.288

Interpretasi: Terdapat variasi skala antar variabel yang cukup besar, misalnya PDRB dan Pengeluaran Per Kapita memiliki rentang nilai jauh lebih lebar dibandingkan RLS atau UHH. Nilai skewness* yang positif pada beberapa variabel (terutama PDRB dan IKK) mengindikasikan distribusi yang right-skewed, menandakan adanya kabupaten/kota dengan nilai ekstrem tinggi yang perlu ditangani dalam tahap preprocessing. Perbedaan skala ini mengonfirmasi kebutuhan normalisasi sebelum analisis clustering.

3.2 Visualisasi Boxplot Seluruh Variabel

df_long <- df_num %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")

p_box <- ggplot(df_long, aes(x = Variabel, y = Nilai)) +
  geom_boxplot(fill = "steelblue", alpha = 0.7, outlier.color = "red", outlier.size = 1.2) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Gambar 1. Boxplot Seluruh Variabel",
       x = "Variabel", y = "Nilai")

print(p_box)

3.3 Visualisasi Histogram dengan Kurva Densitas Seluruh Variabel

# Reshape data ke format panjang
df_all_long <- df_num %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")

# Buat plot dengan facet_wrap (4 kolom x 3 baris)
p_all_hist <- ggplot(df_all_long, aes(x = Nilai)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, 
                 fill = "steelblue", alpha = 0.6, color = "black", linewidth = 0.2) +
  geom_density(color = "darkred", linewidth = 1.2, bw = "nrd0") +
  facet_wrap(~ Variabel, scales = "free", ncol = 3) +
  theme_bw(base_size = 11) +
  theme(strip.background = element_rect(fill = "lightgray"),
        strip.text = element_text(face = "bold")) +
  labs(title = "Distribusi Seluruh Variabel Penelitian",
       subtitle = "Histogram (batang) dengan kurva densitas (garis merah)",
       x = "Nilai", y = "Densitas")

print(p_all_hist)

Gambar 2. Histogram dan Density Plot untuk Seluruh Variabel (Data Asli)

4. Preprocessing Data

Strategi preprocessing yang diterapkan dirancang untuk memaksimalkan kualitas clustering pada data yang bersifat heterogen dan mengandung outlier:

Winsorizing P5–P95: memangkas nilai ekstrem di luar persentil 5 dan 95 untuk mereduksi pengaruh outlier berat.
Transformasi Log1p: menormalisasi distribusi heavy-tailed secara konsisten pada seluruh variabel.
RobustScaler (Median/IQR): menstandarkan skala menggunakan median dan IQR agar tidak terpengaruh outlier residual pascawinsorisasi.

# 1. Winsorizing P5-P95
winsorize_col <- function(x, low = 0.05, high = 0.95) {
  pmax(pmin(x, quantile(x, high, na.rm = TRUE)),
       quantile(x, low, na.rm = TRUE))
}
df_win <- as.data.frame(lapply(df_num, winsorize_col))

# 2. Log1p transform semua variabel
df_log1p <- as.data.frame(mapply(function(col) {
  log1p(col - min(col) + 1)
}, df_win))
names(df_log1p) <- num_cols

# 3. RobustScaler: x_scaled = (x - median) / IQR
robust_scale <- function(x) {
  med <- median(x, na.rm = TRUE)
  iqr <- IQR(x, na.rm = TRUE)
  if (iqr == 0) return(x - med)
  (x - med) / iqr
}
df_scaled    <- as.data.frame(lapply(df_log1p, robust_scale))
df_scaled_mx <- as.matrix(df_scaled)

# Matriks jarak
dist_euc <- dist(df_scaled, method = "euclidean")
dist_man <- dist(df_scaled, method = "manhattan")

Interpretasi: Setelah preprocessing, seluruh variabel berada pada skala yang sebanding dan distribusi yang lebih simetris. Penggunaan RobustScaler (berbasis median/IQR) dipilih karena lebih tahan terhadap outlier residual dibandingkan standardisasi Z-score biasa, sehingga hasil clustering lebih stabil dan representatif.

5. Fungsi Validitas Clustering

Tiga indeks validitas internal digunakan untuk mengevaluasi kualitas clustering:

Silhouette Score: mengukur seberapa mirip sebuah observasi dengan klasternya sendiri dibandingkan klaster lain. Nilai berkisar \([-1, 1]\); semakin tinggi semakin baik (>0.5 = struktur klaster kuat).
Dunn Index: rasio jarak minimum antar-klaster terhadap diameter klaster maksimum. Nilai lebih tinggi menunjukkan klaster yang lebih kompak dan terpisah.
Calinski-Harabasz (CH) Index: rasio dispersi antar-klaster terhadap dispersi dalam-klaster. Nilai lebih tinggi menunjukkan pemisahan klaster yang lebih baik.

get_silhouette <- function(labels, dist_mat) {
  valid <- labels > 0
  if (sum(valid) < 10 || length(unique(labels[valid])) < 2) return(NA_real_)
  sil <- silhouette(as.integer(labels[valid]),
                    as.dist(as.matrix(dist_mat)[valid, valid]))
  round(mean(sil[, 3]), 6)
}

get_dunn <- function(labels, dist_mat) {
  valid <- labels > 0
  if (sum(valid) < 10 || length(unique(labels[valid])) < 2) return(NA_real_)
  tryCatch(
    intCriteria(df_scaled_mx[valid, ], as.integer(labels[valid]), "Dunn")$dunn,
    error = function(e) NA_real_
  )
}

get_ch <- function(labels) {
  valid <- labels > 0
  if (sum(valid) < 10 || length(unique(labels[valid])) < 2) return(NA_real_)
  tryCatch(
    intCriteria(df_scaled_mx[valid, ], as.integer(labels[valid]),
                "Calinski_Harabasz")$calinski_harabasz,
    error = function(e) NA_real_
  )
}

6. Metode Clustering

6.1 K-Means

K-Means adalah algoritma partitional clustering yang meminimalkan jumlah kuadrat jarak (within-cluster sum of squares/WSS) antara setiap observasi dengan pusat klasternya (centroid). Pemilihan \(k\) optimal dilakukan melalui dua pendekatan komplementer: Elbow Method (melihat titik belok pada kurva WSS) dan Silhouette Score (memilih \(k\) dengan nilai rata-rata silhouette tertinggi).

Penentuan k Optimal

set.seed(123)
k_range      <- 2:10
wss_km       <- numeric(length(k_range))
sil_km_range <- numeric(length(k_range))

for (i in seq_along(k_range)) {
  km_i            <- kmeans(df_scaled, centers = k_range[i],
                            nstart = 100, iter.max = 500,
                            algorithm = "Hartigan-Wong")
  wss_km[i]       <- km_i$tot.withinss
  sil_km_range[i] <- mean(silhouette(km_i$cluster, dist_euc)[, 3])
}

k_opt_km <- k_range[which.max(sil_km_range)]

df_elbow <- data.frame(k = k_range, wss = wss_km, sil = sil_km_range)

p_elbow <- ggplot(df_elbow, aes(x = k, y = wss)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(size = 3, color = "steelblue") +
  geom_vline(xintercept = k_opt_km, linetype = "dashed", color = "red") +
  labs(title = "Elbow Method", x = "Jumlah Klaster (k)",
       y = "Total Within-Cluster SS") +
  theme_minimal(base_size = 12)

p_sil_km_plot <- ggplot(df_elbow, aes(x = k, y = sil)) +
  geom_line(color = "darkgreen", linewidth = 1) +
  geom_point(aes(color = k == k_opt_km), size = 3.5) +
  scale_color_manual(values = c("FALSE" = "darkgreen", "TRUE" = "red"),
                     guide = "none") +
  labs(title = "Silhouette Score vs k", x = "Jumlah Klaster (k)",
       y = "Rata-rata Silhouette") +
  theme_minimal(base_size = 12)

grid.arrange(p_elbow, p_sil_km_plot, ncol = 2,
             top = "Gambar 1. Penentuan k Optimal — K-Means")

Gambar 1. Penentuan k Optimal — K-Means

Interpretasi: Panel kiri (Elbow Method) menunjukkan penurunan WSS yang melambat secara signifikan pada \(k\) tertentu, titik “siku” ini mengindikasikan jumlah klaster yang efisien. Panel kanan mempertegas pilihan tersebut: \(k\) = 3 menghasilkan Silhouette Score tertinggi (ditandai titik merah), menandakan klaster yang paling kohesif dan terpisah satu sama lain pada algoritma K-Means.

Model Final K-Means

set.seed(42)
km_final <- kmeans(df_scaled, centers = k_opt_km,
                   nstart = 200, iter.max = 1000,
                   algorithm = "Hartigan-Wong")

sil_km  <- get_silhouette(km_final$cluster, dist_euc)
dunn_km <- get_dunn(km_final$cluster, dist_euc)
ch_km   <- get_ch(km_final$cluster)

km_summary <- data.frame(
  Parameter   = c("Jumlah Klaster (k)", "Silhouette Score", "Dunn Index",
                  "Calinski-Harabasz Index", "Distribusi Klaster"),
  Nilai       = c(k_opt_km,
                  round(sil_km, 4),
                  round(dunn_km, 4),
                  round(ch_km, 2),
                  paste(table(km_final$cluster), collapse = " / "))
)

kable(km_summary, col.names = c("Parameter", "Nilai"),
      caption = "Tabel 2. Ringkasan Hasil K-Means Final") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 2. Ringkasan Hasil K-Means Final
Parameter	Nilai
Jumlah Klaster (k)	3
Silhouette Score	0.6876
Dunn Index	0.2293
Calinski-Harabasz Index	378.63
Distribusi Klaster	26 / 466 / 22

Interpretasi: Model final K-Means dijalankan dengan nstart = 200 untuk memastikan konvergensi ke solusi global yang optimal. Silhouette Score sebesar 0.6876 menunjukkan kualitas pemisahan klaster; nilai Dunn dan CH melengkapi evaluasi dari sisi kompaktness dan separabilitas klaster.

6.2 K-Median

K-Median merupakan varian dari K-Means yang menggunakan median (bukan mean) sebagai pusat klaster dan jarak Manhattan sebagai ukuran kedekatan. Pendekatan ini lebih robust terhadap outlier karena median tidak sensitif terhadap nilai ekstrem, menjadikannya pilihan yang tepat untuk data dengan distribusi miring seperti indikator kemiskinan.

Implementasi dan Penentuan k Optimal

# Implementasi K-Median manual (median sebagai centroid, jarak Manhattan)
kmedian <- function(data, k, max_iter = 150, n_init = 30) {
  X <- as.matrix(data); n <- nrow(X); p <- ncol(X)
  best_obj <- Inf; best_out <- NULL

  for (run in seq_len(n_init)) {
    set.seed(run * 13 + 7)
    centers <- X[sample(n, k), , drop = FALSE]
    labels  <- integer(n)

    for (it in seq_len(max_iter)) {
      dmat <- matrix(0, n, k)
      for (j in seq_len(k))
        dmat[, j] <- rowSums(abs(sweep(X, 2, centers[j, ], "-")))
      new_labels <- max.col(-dmat)
      new_centers <- matrix(0, k, p)
      for (j in seq_len(k)) {
        mem <- X[new_labels == j, , drop = FALSE]
        new_centers[j, ] <- if (nrow(mem) > 0) apply(mem, 2, median) else centers[j, ]
      }
      if (all(new_labels == labels) &&
          max(abs(new_centers - centers)) < 1e-8) break
      labels <- new_labels; centers <- new_centers
    }

    obj <- sum(vapply(seq_len(k), function(j) {
      mem <- X[labels == j, , drop = FALSE]
      if (nrow(mem) == 0) return(0)
      sum(rowSums(abs(sweep(mem, 2, centers[j, ], "-"))))
    }, numeric(1)))

    if (obj < best_obj) {
      best_obj <- obj
      best_out <- list(cluster = labels, centers = centers)
    }
  }
  best_out
}

# Pencarian k optimal
sil_kmed_range <- numeric(length(k_range))
for (i in seq_along(k_range)) {
  res <- kmedian(df_scaled, k_range[i], n_init = 15)
  if (is.null(res) || length(unique(res$cluster)) < 2) {
    sil_kmed_range[i] <- NA; next
  }
  sil_kmed_range[i] <- mean(silhouette(res$cluster, dist_man)[, 3])
}

k_opt_kmed <- k_range[which.max(sil_kmed_range)]
kmed_final <- kmedian(df_scaled, k_opt_kmed, n_init = 50)
sil_kmed   <- get_silhouette(kmed_final$cluster, dist_man)
dunn_kmed  <- get_dunn(kmed_final$cluster, dist_man)
ch_kmed    <- get_ch(kmed_final$cluster)

kmed_summary <- data.frame(
  Parameter = c("Jumlah Klaster (k)", "Jarak yang Digunakan", "Silhouette Score",
                "Dunn Index", "Calinski-Harabasz Index", "Distribusi Klaster"),
  Nilai     = c(k_opt_kmed, "Manhattan",
                round(sil_kmed, 4), round(dunn_kmed, 4), round(ch_kmed, 2),
                paste(table(kmed_final$cluster), collapse = " / "))
)

kable(kmed_summary, col.names = c("Parameter", "Nilai"),
      caption = "Tabel 3. Ringkasan Hasil K-Median Final") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 3. Ringkasan Hasil K-Median Final
Parameter	Nilai
Jumlah Klaster (k)	2
Jarak yang Digunakan	Manhattan
Silhouette Score	0.6449
Dunn Index	0.1766
Calinski-Harabasz Index	280.57
Distribusi Klaster	31 / 483

Interpretasi: K-Median dengan jarak Manhattan menghasilkan klaster yang lebih tahan terhadap pengaruh kabupaten/kota dengan nilai ekstrem (misalnya daerah dengan PDRB sangat tinggi). Nilai \(k\) optimal = 2 dipilih berdasarkan Silhouette Score tertinggi. Perbandingan antara K-Means dan K-Median akan memberikan gambaran tentang seberapa besar outlier memengaruhi struktur pengelompokan.

6.3 DBSCAN

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) adalah algoritma berbasis kepadatan yang mengelompokkan observasi berdasarkan kerapatan titik di sekitarnya. Keunggulan utama DBSCAN adalah kemampuannya mendeteksi klaster berbentuk arbitrer dan mengidentifikasi noise/outlier secara otomatis tanpa perlu menentukan jumlah klaster terlebih dahulu. Parameter utamanya adalah eps (radius lingkungan) dan minPts (jumlah minimum titik dalam radius).

Visualisasi k-Distance Plot

# Pastikan data dalam bentuk matriks numerik
df_scaled_mat <- as.matrix(df_scaled)

# Hitung jarak ke tetangga ke-5
knn_dist <- dbscan::kNNdist(df_scaled_mat, k = 5)

# Karena k=5, hasil adalah matriks 514x5. Ambil kolom ke-5.
if (is.matrix(knn_dist) && ncol(knn_dist) >= 5) {
  k_dist_vec <- knn_dist[, 5]
} else {
  # Fallback jika k=1 (tidak mungkin di sini)
  k_dist_vec <- knn_dist
}

knn_sorted <- sort(k_dist_vec)

p_kdist <- ggplot(
  data.frame(Index = seq_along(knn_sorted), KDist = knn_sorted),
  aes(x = Index, y = KDist)
) +
  geom_line(color = "steelblue", linewidth = 0.8) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "red", linewidth = 0.8) +
  labs(title = "Gambar 2. k-Distance Plot (k = 5) — Panduan Pemilihan eps DBSCAN",
       x = "Titik (diurutkan)", y = "Jarak ke Tetangga ke-5") +
  theme_minimal(base_size = 12)

print(p_kdist)

Gambar 2. k-Distance Plot untuk Penentuan eps Optimal (DBSCAN)

Interpretasi: k-Distance Plot menampilkan jarak setiap titik ke tetangga ke-5 terdekatnya, diurutkan secara menaik. Titik “siku” pada kurva ini menjadi panduan pemilihan nilai eps yang optimal, di bawah titik siku, titik-titik berdekatan membentuk klaster; di atasnya, titik dianggap sebagai noise. Garis putus-putus merah menandai area referensi untuk eps.

Tuning Parameter DBSCAN

df_scaled_mat <- as.matrix(df_scaled)
dimnames(df_scaled_mat) <- NULL
dist_man_mat  <- as.matrix(dist_man)

eps_candidates   <- seq(0.3, 2.2, by = 0.1)
minPts_candidates <- c(3, 4, 5, 6, 7, 8, 9, 10, 12, 15)

best_sil_db <- -Inf; best_eps_db <- NA; best_mp_db <- NA
best_db_out <- NULL; best_metric_db <- "euclidean"

# Pencarian dengan Euclidean distance
for (mp in minPts_candidates) {
  for (ep in eps_candidates) {
    db  <- dbscan::dbscan(df_scaled_mat, eps = ep, minPts = mp)
    nc  <- length(unique(db$cluster[db$cluster > 0]))
    if (nc < 2 || nc > 20) next
    if (mean(db$cluster == 0) > 0.4) next
    sv  <- get_silhouette(db$cluster, dist_euc)
    if (!is.na(sv) && sv > best_sil_db) {
      best_sil_db <- sv; best_eps_db <- ep; best_mp_db <- mp
      best_db_out <- db; best_metric_db <- "euclidean"
    }
  }
}

# Pencarian dengan Manhattan distance
for (mp in minPts_candidates) {
  for (ep in eps_candidates) {
    db  <- fpc::dbscan(dist_man_mat, eps = ep, MinPts = mp, method = "dist")
    nc  <- length(unique(db$cluster[db$cluster > 0]))
    if (nc < 2 || nc > 20) next
    if (mean(db$cluster == 0) > 0.4) next
    sv  <- get_silhouette(db$cluster, dist_man)
    if (!is.na(sv) && sv > best_sil_db) {
      best_sil_db <- sv; best_eps_db <- ep; best_mp_db <- mp
      best_db_out <- db; best_metric_db <- "manhattan"
    }
  }
}

# Fallback jika tidak ada kombinasi memenuhi syarat
if (is.null(best_db_out)) {
  db <- dbscan::dbscan(df_scaled_mat, eps = 0.8, minPts = 5)
  best_db_out <- db; best_sil_db <- get_silhouette(db$cluster, dist_euc)
  best_eps_db <- 0.8; best_mp_db <- 5; best_metric_db <- "euclidean (fallback)"
}

nc_db     <- length(unique(best_db_out$cluster[best_db_out$cluster > 0]))
nnoise_db <- sum(best_db_out$cluster == 0)
dunn_db   <- get_dunn(best_db_out$cluster,
                       if (best_metric_db == "euclidean") dist_euc else dist_man)
ch_db     <- get_ch(best_db_out$cluster)
best_sil_db <- ifelse(is.na(best_sil_db), NA, best_sil_db)

db_summary <- data.frame(
  Parameter = c("eps Optimal", "minPts Optimal", "Metrik Jarak",
                "Jumlah Klaster", "Jumlah Noise",
                "Silhouette Score", "Dunn Index", "Calinski-Harabasz Index"),
  Nilai     = c(round(best_eps_db, 2), best_mp_db, best_metric_db,
                nc_db, sprintf("%d (%.1f%%)", nnoise_db, nnoise_db / nrow(df_scaled) * 100),
                round(best_sil_db, 4), round(dunn_db, 4), round(ch_db, 2))
)

kable(db_summary, col.names = c("Parameter", "Nilai"),
      caption = "Tabel 4. Ringkasan Hasil DBSCAN Terbaik") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 4. Ringkasan Hasil DBSCAN Terbaik
Parameter	Nilai
eps Optimal	1.5
minPts Optimal	4
Metrik Jarak	euclidean
Jumlah Klaster	2
Jumlah Noise	100 (19.5%)
Silhouette Score	0.7714
Dunn Index	1.5227
Calinski-Harabasz Index	168.86

Interpretasi: Tuning DBSCAN dilakukan secara ekstensif melalui grid search pada kombinasi eps (0.3–2.2) dan minPts (3–15) untuk kedua metrik jarak. Parameter optimal eps = 1.5 dan minPts = 4 menghasilkan 2 klaster dengan 100 titik noise (19.5%). Kehadiran noise mencerminkan kabupaten/kota dengan profil kemiskinan yang sangat berbeda dari pola umum, hal ini merupakan informasi berharga untuk identifikasi wilayah yang membutuhkan penanganan khusus.

6.4 Mean Shift

Mean Shift adalah algoritma berbasis kepadatan yang secara iteratif menggeser setiap titik menuju daerah dengan kepadatan tertinggi (mode) di sekitarnya. Tidak seperti K-Means, Mean Shift tidak memerlukan spesifikasi jumlah klaster, algoritma secara otomatis menemukan jumlah dan bentuk klaster berdasarkan parameter bandwidth (h) yang menentukan lebar “jendela” pencarian.

Tuning Bandwidth

# Grid bandwidth berbasis Silverman dan rentang fixed
bw_base <- 1.06 * mean(apply(df_scaled, 2, sd)) *
           nrow(df_scaled)^(-1 / (4 + ncol(df_scaled)))

bw_grid <- sort(unique(c(
  bw_base * c(0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0),
  seq(1.0, 8.0, by = 0.5)
)))
bw_grid <- bw_grid[bw_grid > 0]

best_sil_ms <- -Inf; best_bw_ms <- NA
best_ms_out <- NULL; best_lab_ms <- NULL

for (bw in bw_grid) {
  ms <- tryCatch(
    meanshift(df_scaled, h = bw, iter = 200, tol = 1e-6),
    error = function(e) NULL
  )
  if (is.null(ms)) next
  labs <- ms$labels + 1L
  nc   <- length(unique(labs))
  if (nc < 2 || nc > 20) next
  sv <- tryCatch(mean(silhouette(labs, dist_euc)[, 3]), error = function(e) NA)
  if (!is.na(sv) && sv > best_sil_ms) {
    best_sil_ms <- sv; best_bw_ms <- bw
    best_ms_out <- ms; best_lab_ms <- labs
  }
}

if (!is.null(best_ms_out)) {
  nc_ms   <- length(unique(best_lab_ms))
  dunn_ms <- get_dunn(best_lab_ms, dist_euc)
  ch_ms   <- get_ch(best_lab_ms)
} else {
  best_sil_ms <- NA; dunn_ms <- NA; ch_ms <- NA; nc_ms <- NA
}

ms_summary <- data.frame(
  Parameter = c("Bandwidth (h) Optimal", "Jumlah Klaster",
                "Silhouette Score", "Dunn Index", "Calinski-Harabasz Index"),
  Nilai     = c(ifelse(is.na(best_bw_ms), "Tidak valid", round(best_bw_ms, 4)),
                ifelse(is.na(nc_ms), "Tidak valid", nc_ms),
                ifelse(is.na(best_sil_ms), "NA", round(best_sil_ms, 4)),
                ifelse(is.na(dunn_ms),    "NA", round(dunn_ms, 4)),
                ifelse(is.na(ch_ms),      "NA", round(ch_ms, 2)))
)

kable(ms_summary, col.names = c("Parameter", "Nilai"),
      caption = "Tabel 5. Ringkasan Hasil Mean Shift Terbaik") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 5. Ringkasan Hasil Mean Shift Terbaik
Parameter	Nilai
Bandwidth (h) Optimal	Tidak valid
Jumlah Klaster	Tidak valid
Silhouette Score	NA
Dunn Index	NA
Calinski-Harabasz Index	NA

Interpretasi: Pencarian bandwidth dilakukan pada 25 nilai berbasis aturan Silverman dan rentang fixed (1.0–8.0). Bandwidth optimal tidak ditemukan menghasilkan NA klaster. Bandwidth yang lebih besar cenderung menghasilkan lebih sedikit klaster (penggabungan mode yang berdekatan), sementara bandwidth kecil menghasilkan terlalu banyak klaster kecil. Nilai Silhouette yang diperoleh mencerminkan sejauh mana struktur densitas data mendukung pemisahan klaster yang jelas.

6.5 Fuzzy C-Means

Fuzzy C-Means (FCM) merupakan metode soft clustering di mana setiap observasi tidak ditugaskan secara eksklusif ke satu klaster, melainkan memiliki derajat keanggotaan di semua klaster. Parameter fuzziness \(m\) (\(m > 1\)) mengontrol tingkat “keaburaman”, semakin tinggi \(m\), semakin merata distribusi keanggotaan; nilai \(m = 2\) adalah standar yang dikemukakan Bezdek (1981). FCM ideal untuk data di mana batas antarklaster tidak tajam, seperti pada indikator kemiskinan yang cenderung kontinu.

Tuning Parameter c dan m

m_grid <- c(1.5, 2.0, 2.5)
c_grid <- 2:7

best_sil_fcm <- -Inf; best_c_fcm <- NA; best_m_fcm <- NA; best_fcm_out <- NULL

for (c_i in c_grid) {
  for (m_i in m_grid) {
    set.seed(2024)
    fcm <- tryCatch(
      cmeans(df_scaled, centers = c_i, m = m_i, iter.max = 500,
             dist = "euclidean", method = "cmeans"),
      error = function(e) NULL
    )
    if (is.null(fcm) || length(unique(fcm$cluster)) < 2) next
    sv <- tryCatch(mean(silhouette(fcm$cluster, dist_euc)[, 3]), error = function(e) NA)
    if (!is.na(sv) && sv > best_sil_fcm) {
      best_sil_fcm <- sv; best_c_fcm <- c_i; best_m_fcm <- m_i; best_fcm_out <- fcm
    }
  }
}

if (!is.null(best_fcm_out)) {
  dunn_fcm <- get_dunn(best_fcm_out$cluster, dist_euc)
  ch_fcm   <- get_ch(best_fcm_out$cluster)
} else {
  best_sil_fcm <- NA; dunn_fcm <- NA; ch_fcm <- NA
}

fcm_summary <- data.frame(
  Parameter = c("Jumlah Klaster (c) Optimal", "Fuzziness (m) Optimal",
                "Silhouette Score", "Dunn Index", "Calinski-Harabasz Index",
                "Distribusi Klaster (Hard Assignment)"),
  Nilai     = c(ifelse(is.na(best_c_fcm), "NA", best_c_fcm),
                ifelse(is.na(best_m_fcm), "NA", best_m_fcm),
                ifelse(is.na(best_sil_fcm), "NA", round(best_sil_fcm, 4)),
                ifelse(is.na(dunn_fcm),    "NA", round(dunn_fcm, 4)),
                ifelse(is.na(ch_fcm),      "NA", round(ch_fcm, 2)),
                ifelse(is.null(best_fcm_out), "NA",
                       paste(table(best_fcm_out$cluster), collapse = " / ")))
)

kable(fcm_summary, col.names = c("Parameter", "Nilai"),
      caption = "Tabel 6. Ringkasan Hasil Fuzzy C-Means Terbaik") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13) %>%
  column_spec(1, bold = TRUE)

Tabel 6. Ringkasan Hasil Fuzzy C-Means Terbaik
Parameter	Nilai
Jumlah Klaster (c) Optimal	2
Fuzziness (m) Optimal	1.5
Silhouette Score	0.6764
Dunn Index	0.1253
Calinski-Harabasz Index	371.52
Distribusi Klaster (Hard Assignment)	465 / 49

Interpretasi: Grid search pada kombinasi \(c \in \{2, ..., 7\}\) dan \(m \in \{1.5, 2.0, 2.5\}\) mengidentifikasi parameter optimal \(c = 2\) dan \(m = 1.5\). Nilai \(m\) optimal memberikan gambaran tentang ambiguitas batas klaster dalam data, nilai \(m\) yang lebih kecil menandakan klaster yang lebih tegas (mendekati hard clustering), cocok jika ada pemisahan yang nyata antarkelompok kabupaten/kota.

7. Perbandingan Metode

7.1 Tabel Perbandingan Indeks Validitas

comparison_tbl <- data.frame(
  Metode    = c("K-Means", "K-Median", "DBSCAN", "Mean Shift", "Fuzzy C-Means"),
  Parameter = c(
    sprintf("k = %d", k_opt_km),
    sprintf("k = %d (Manhattan)", k_opt_kmed),
    if (!is.null(best_db_out))
      sprintf("eps = %.2f, minPts = %d (%s)", best_eps_db, best_mp_db, best_metric_db)
    else "Tidak valid",
    if (!is.null(best_ms_out)) sprintf("bw = %.4f", best_bw_ms) else "Tidak valid",
    if (!is.null(best_fcm_out))
      sprintf("c = %d, m = %.2f", best_c_fcm, best_m_fcm) else "Tidak valid"
  ),
  n_Cluster = c(
    k_opt_km, k_opt_kmed,
    ifelse(is.null(best_db_out),  NA, nc_db),
    ifelse(is.null(best_ms_out),  NA, nc_ms),
    ifelse(is.null(best_fcm_out), NA, best_c_fcm)
  ),
  Silhouette = round(c(
    sil_km, sil_kmed,
    ifelse(is.null(best_db_out),  NA, best_sil_db),
    ifelse(is.null(best_ms_out),  NA, best_sil_ms),
    ifelse(is.null(best_fcm_out), NA, best_sil_fcm)
  ), 4),
  Dunn = round(c(
    dunn_km, dunn_kmed,
    ifelse(is.null(best_db_out),  NA, dunn_db),
    ifelse(is.null(best_ms_out),  NA, dunn_ms),
    ifelse(is.null(best_fcm_out), NA, dunn_fcm)
  ), 4),
  CH_Index = round(c(
    ch_km, ch_kmed,
    ifelse(is.null(best_db_out),  NA, ch_db),
    ifelse(is.null(best_ms_out),  NA, ch_ms),
    ifelse(is.null(best_fcm_out), NA, ch_fcm)
  ), 2)
)

best_idx     <- which.max(comparison_tbl$Silhouette)
best_method  <- comparison_tbl$Metode[best_idx]
best_sil_val <- comparison_tbl$Silhouette[best_idx]

kable(comparison_tbl,
      col.names = c("Metode", "Parameter Terbaik", "Jumlah Klaster",
                    "Silhouette", "Dunn", "CH Index"),
      caption = "Tabel 7. Perbandingan Performa Semua Metode Clustering",
      align = c("l", "l", "c", "c", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE, font_size = 13) %>%
  row_spec(best_idx, bold = TRUE, background = "#fff3cd") %>%
  column_spec(1, bold = TRUE)

Tabel 7. Perbandingan Performa Semua Metode Clustering
Metode	Parameter Terbaik	Jumlah Klaster	Silhouette	Dunn	CH Index
K-Means	k = 3	3	0.6876	0.2293	378.63
K-Median	k = 2 (Manhattan)	2	0.6449	0.1766	280.57
DBSCAN	eps = 1.50, minPts = 4 (euclidean)	2	0.7714	1.5227	168.86
Mean Shift	Tidak valid	NA	NA	NA	NA
Fuzzy C-Means	c = 2, m = 1.50	2	0.6764	0.1253	371.52

Metode Terbaik: DBSCAN dengan Silhouette Score = 0.7714 .

7.2 Visualisasi Perbandingan Silhouette

comp_plot_df <- comparison_tbl %>%
  filter(!is.na(Silhouette)) %>%
  mutate(Metode  = reorder(Metode, Silhouette),
         Terbaik = (as.character(Metode) == best_method))

p5 <- ggplot(comp_plot_df, aes(x = Metode, y = Silhouette, fill = Terbaik)) +
  geom_bar(stat = "identity", alpha = 0.85, color = "white") +
  geom_text(aes(label = sprintf("%.4f", Silhouette)),
            hjust = -0.1, size = 4, fontface = "bold") +
  geom_hline(yintercept = 0.5, linetype = "dashed",
             color = "gray50", linewidth = 0.8) +
  coord_flip() +
  scale_fill_manual(values = c("FALSE" = "steelblue", "TRUE" = "gold3"),
                    guide = "none") +
  ylim(0, min(1.0, max(comp_plot_df$Silhouette, na.rm = TRUE) * 1.25)) +
  labs(title  = "Gambar 5. Perbandingan Silhouette Score Semua Metode",
       subtitle = "Garis putus-putus = ambang batas 0.5 (struktur klaster kuat)",
       x = "Metode", y = "Silhouette Score") +
  theme_minimal(base_size = 12)

print(p5)

Gambar 5. Perbandingan Silhouette Score Semua Metode Clustering

Interpretasi: Bar berwarna emas menandai metode terbaik (DBSCAN). Garis putus-putus pada nilai 0.5 adalah ambang batas umum: di atas 0.5 menandakan struktur klaster yang kuat dan dapat diandalkan untuk interpretasi. Metode yang jatuh di bawah ambang ini masih menghasilkan klaster yang valid, namun dengan pemisahan yang kurang tegas, bisa diakibatkan oleh asumsi bentuk klaster yang tidak sesuai dengan distribusi data kemiskinan.

8. Hasil Clustering Terbaik

8.1 Penetapan Metode Terbaik

best_labels <- switch(best_method,
  "K-Means"       = km_final$cluster,
  "K-Median"      = kmed_final$cluster,
  "DBSCAN"        = best_db_out$cluster,
  "Mean Shift"    = best_lab_ms,
  "Fuzzy C-Means" = best_fcm_out$cluster
)

df_result         <- df_raw
df_result$Cluster <- best_labels

Metode terpilih adalah DBSCAN dengan Silhouette Score = 0.7714.

8.2 Distribusi Kabupaten/Kota per Cluster

dist_cluster <- as.data.frame(table(Cluster = best_labels)) %>%
  rename(Jumlah_KabKota = Freq) %>%
  mutate(Persentase = sprintf("%.1f%%", Jumlah_KabKota / sum(Jumlah_KabKota) * 100))

kable(dist_cluster,
      caption = "Tabel 8. Distribusi Kabupaten/Kota per Cluster",
      align = c("c", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, position = "left", font_size = 13)

Tabel 8. Distribusi Kabupaten/Kota per Cluster
Cluster	Jumlah_KabKota	Persentase
0	100	19.5%
1	409	79.6%
2	5	1.0%

Interpretasi: Distribusi anggota klaster menunjukkan apakah pengelompokan bersifat seimbang atau ada klaster dominan. Klaster dengan anggota sangat sedikit (< 5) perlu ditelaah lebih lanjut, karena bisa merupakan kelompok “ekstrem” yang memang berbeda secara signifikan, atau artefak dari parameter yang kurang optimal.

8.3 Daftar Kabupaten/Kota per Cluster

valid_cls <- sort(unique(best_labels[best_labels > 0]))
for (cl in valid_cls) {
  kab_cl <- df_result$Kabupaten_Kota[df_result$Cluster == cl]
  cat(sprintf("\n**Cluster %d** (%d Kab/Kota):\n\n", cl, length(kab_cl)))
  cat(paste(kab_cl, collapse = ", "), "\n\n")
}

Cluster 1 (409 Kab/Kota):

Aceh Singkil, Aceh Selatan, Aceh Tenggara, Aceh Timur, Aceh Tengah, Aceh Barat, Aceh Besar, Pidie, Bireuen, Aceh Utara, Aceh Barat Daya, Gayo Lues, Nagan Raya, Aceh Jaya, Bener Meriah, Pidie Jaya, Kota Banda Aceh, Kota Langsa, Kota Lhokseumawe, Mandailing Natal, Tapanuli Selatan, Tapanuli Tengah, Tapanuli Utara, Toba Samosir, Labuhan Batu, Asahan, Simalungun, Dairi, Karo, Deli Serdang, Langkat, Humbang Hasundutan, Samosir, Serdang Bedagai, Batu Bara, Padang Lawas Utara, Padang Lawas, Labuhan Batu Selatan, Labuhan Batu Utara, Kota Sibolga, Kota Tanjung Balai, Kota Pematang Siantar, Kota Tebing Tinggi, Kota Medan, Kota Binjai, Kota Padangsidimpuan, Pesisir Selatan, Solok, Sijunjung, Tanah Datar, Padang Pariaman, Agam, Lima Puluh Kota, Pasaman, Solok Selatan, Dharmasraya, Pasaman Barat, Kota Padang, Kota Solok, Kota Sawah Lunto, Kota Padang Panjang, Kota Bukittinggi, Kota Payakumbuh, Kota Pariaman, Kuantan Singingi, Indragiri Hulu, Indragiri Hilir, Pelalawan, Siak, Kampar, Rokan Hulu, Bengkalis, Rokan Hilir, Kepulauan Meranti, Kota Pekanbaru, Kota Dumai, Kerinci, Merangin, Sarolangun, Batang Hari, Muaro Jambi, Tanjung Jabung Timur, Tanjung Jabung Barat, Tebo, Bungo, Kota Jambi, Kota Sungai Penuh, Ogan Komering Ulu, Ogan Komering Ilir, Muara Enim, Lahat, Musi Rawas, Musi Banyuasin, Banyu Asin, Ogan Komering Ulu Selatan, Ogan Komering Ulu Timur, Ogan Ilir, Empat Lawang, Penukal Abab Lematang Ilir, Musi Rawas Utara, Kota Palembang, Kota Prabumulih, Kota Pagar Alam, Kota Lubuklinggau, Bengkulu Selatan, Rejang Lebong, Bengkulu Utara, Kaur, Seluma, Mukomuko, Lebong, Kepahiang, Bengkulu Tengah, Kota Bengkulu, Lampung Barat, Tanggamus, Lampung Selatan, Lampung Timur, Lampung Tengah, Lampung Utara, Way Kanan, Tulangbawang, Pesawaran, Pringsewu, Mesuji, Tulang Bawang Barat, Pesisir Barat, Kota Bandar Lampung, Kota Metro, Bangka, Belitung, Bangka Barat, Bangka Tengah, Bangka Selatan, Belitung Timur, Kota Pangkal Pinang, Karimun, Bintan, Natuna, Lingga, Kepulauan Anambas, Kota Batam, Kota Tanjung Pinang, Kepulauan Seribu, Kota Jakarta Selatan, Kota Jakarta Timur, Kota Jakarta Pusat, Kota Jakarta Barat, Kota Jakarta Utara, Bogor, Sukabumi, Cianjur, Bandung, Garut, Tasikmalaya, Ciamis, Kuningan, Cirebon, Majalengka, Sumedang, Indramayu, Subang, Purwakarta, Karawang, Bekasi, Bandung Barat, Pangandaran, Kota Bogor, Kota Sukabumi, Kota Bandung, Kota Cirebon, Kota Bekasi, Kota Depok, Kota Cimahi, Kota Tasikmalaya, Kota Banjar, Cilacap, Banyumas, Purbalingga, Banjarnegara, Kebumen, Purworejo, Wonosobo, Magelang, Boyolali, Klaten, Sukoharjo, Wonogiri, Karanganyar, Sragen, Grobogan, Blora, Rembang, Pati, Kudus, Jepara, Demak, Semarang, Temanggung, Kendal, Batang, Pekalongan, Pemalang, Tegal, Brebes, Kota Magelang, Kota Surakarta, Kota Salatiga, Kota Semarang, Kota Pekalongan, Kota Tegal, Kulon Progo, Bantul, Gunung Kidul, Sleman, Kota Yogyakarta, Ponorogo, Tulungagung, Blitar, Kediri, Malang, Lumajang, Jember, Banyuwangi, Bondowoso, Situbondo, Probolinggo, Pasuruan, Sidoarjo, Mojokerto, Jombang, Nganjuk, Madiun, Magetan, Ngawi, Bojonegoro, Tuban, Lamongan, Gresik, Bangkalan, Sampang, Pamekasan, Kota Kediri, Kota Blitar, Kota Malang, Kota Probolinggo, Kota Pasuruan, Kota Mojokerto, Kota Madiun, Kota Surabaya, Kota Batu, Tangerang, Kota Tangerang, Kota Cilegon, Kota Serang, Kota Tangerang Selatan, Jembrana, Tabanan, Badung, Gianyar, Bangli, Buleleng, Kota Denpasar, Lombok Barat, Lombok Tengah, Lombok Timur, Sumbawa, Dompu, Bima, Sumbawa Barat, Kota Mataram, Kota Bima, Sumba Timur, Manggarai, Manggarai Barat, Sambas, Bengkayang, Landak, Pontianak, Sanggau, Ketapang, Sintang, Kapuas Hulu, Sekadau, Melawi, Kubu Raya, Kota Pontianak, Kota Singkawang, Kotawaringin Barat, Kotawaringin Timur, Kapuas, Barito Selatan, Barito Utara, Sukamara, Lamandau, Seruyan, Katingan, Pulang Pisau, Gunung Mas, Barito Timur, Murung Raya, Kota Palangka Raya, Tanah Laut, Kota Baru, Banjar, Barito Kuala, Tapin, Hulu Sungai Selatan, Hulu Sungai Tengah, Tabalong, Tanah Bumbu, Balangan, Kota Banjarmasin, Kota Banjar Baru, Paser, Kutai Barat, Kutai Kartanegara, Kutai Timur, Berau, Penajam Paser Utara, Kota Balikpapan, Kota Samarinda, Kota Bontang, Malinau, Bulungan, Nunukan, Kota Tarakan, Bolaang Mongondow, Minahasa, Kepulauan Sangihe, Minahasa Selatan, Minahasa Utara, Bolaang Mongondow Utara, Minahasa Tenggara, Bolaang Mongondow Timur, Kota Manado, Kota Bitung, Kota Tomohon, Kota Kotamobagu, Banggai Kepulauan, Banggai, Morowali, Poso, Donggala, Toli-Toli, Buol, Parigi Moutong, Tojo Una-Una, Sigi, Banggai Laut, Morowali Utara, Kota Palu, Kepulauan Selayar, Bulukumba, Bantaeng, Jeneponto, Takalar, Gowa, Sinjai, Maros, Pangkajene Dan Kepulauan, Barru, Bone, Soppeng, Wajo, Sidenreng Rappang, Pinrang, Enrekang, Luwu, Tana Toraja, Luwu Utara, Luwu Timur, Toraja Utara, Kota Makassar, Kota Parepare, Kota Palopo, Konawe, Kolaka, Konawe Selatan, Bombana, Wakatobi, Kolaka Utara, Buton Utara, Konawe Utara, Kolaka Timur, Muna Barat, Kota Kendari, Kota Baubau, Gorontalo, Pohuwato, Bone Bolango, Kota Gorontalo, Polewali Mandar, Mamuju, Mamuju Utara, Mamuju Tengah, Maluku Tengah, Buru, Kota Ambon, Halmahera Barat, Halmahera Utara, Kota Ternate, Kota Tidore Kepulauan, Manokwari, Kota Sorong, Jayapura, Biak Numfor, Merauke, Mimika

Cluster 2 (5 Kab/Kota):

Manokwari Selatan, Pegunungan Arfak, Tambrauw, Intan Jaya, Puncak Jaya

if (any(best_labels == 0)) {
  noise_kab <- df_result$Kabupaten_Kota[df_result$Cluster == 0]
  cat(sprintf("\n**Noise/Outlier** (%d titik):\n\n", length(noise_kab)))
  cat(paste(noise_kab, collapse = ", "), "\n\n")
}

Noise/Outlier (100 titik):

Simeulue, Aceh Tamiang, Kota Sabang, Kota Subulussalam, Nias, Nias Selatan, Pakpak Bharat, Nias Utara, Nias Barat, Kota Gunungsitoli, Kepulauan Mentawai, Pacitan, Trenggalek, Sumenep, Pandeglang, Lebak, Serang, Klungkung, Karang Asem, Lombok Utara, Sumba Barat, Kupang, Timor Tengah Selatan, Timor Tengah Utara, Belu, Alor, Lembata, Flores Timur, Sikka, Ende, Ngada, Rote Ndao, Sumba Tengah, Sumba Barat Daya, Nagekeo, Manggarai Timur, Sabu Raijua, Malaka, Kota Kupang, Kayong Utara, Hulu Sungai Utara, Mahakam Hulu, Tana Tidung, Kepulauan Talaud, Siau Tagulandang Biaro, Bolaang Mongondow Selatan, Buton, Muna, Konawe Kepulauan, Buton Tengah, Buton Selatan, Boalemo, Gorontalo Utara, Majene, Mamasa, Maluku Tenggara Barat, Maluku Tenggara, Kepulauan Aru, Seram Bagian Barat, Seram Bagian Timur, Maluku Barat Daya, Buru Selatan, Kota Tual, Halmahera Tengah, Kepulauan Sula, Halmahera Selatan, Halmahera Timur, Pulau Morotai, Pulau Taliabu, Fakfak, Kaimana, Teluk Wondama, Teluk Bintuni, Raja Ampat, Sorong, Sorong Selatan, Maybrat, Kepulauan Yapen, Sarmi, Keerom, Waropen, Supiori, Mamberamo Raya, Kota Jayapura, Boven Digoel, Mappi, Asmat, Dogiyai, Deiyai, Nabire, Paniai, Puncak, Nduga, Jayawijaya, Lanny Jaya, Tolikara, Mamberamo Tengah, Yalimo, Yahukimo, Pegunungan Bintang

8.4 Visualisasi Cluster (PCA 2D)

pca_res <- prcomp(df_scaled, scale. = FALSE)
pca_df  <- as.data.frame(pca_res$x[, 1:2])
pca_df$Cluster <- as.factor(best_labels)
var_exp <- round(summary(pca_res)$importance[2, 1:2] * 100, 1)

p3 <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(alpha = 0.6, size = 1.8) +
  stat_ellipse(data = subset(pca_df, Cluster != "0"),
               aes(group = Cluster), type = "norm", level = 0.75,
               linewidth = 0.8, linetype = "dashed") +
  scale_color_brewer(palette = "Set1") +
  labs(title    = sprintf("Gambar 3. Visualisasi Cluster — %s", best_method),
       subtitle = sprintf("Sil = %.4f  |  PC1 = %.1f%%  |  PC2 = %.1f%%",
                          best_sil_val, var_exp[1], var_exp[2]),
       x = sprintf("PC1 (%.1f%% variansi)", var_exp[1]),
       y = sprintf("PC2 (%.1f%% variansi)", var_exp[2]),
       color = "Cluster") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "right")

print(p3)

Gambar 3. Visualisasi Cluster dalam Ruang PCA 2D

Interpretasi: PCA mereduksi dimensi data ke dua komponen utama yang menjelaskan 51.8% (PC1) dan 26% (PC2) variansi total. Titik-titik yang berwarna berbeda mewakili klaster berbeda; elips putus-putus menunjukkan batas 75% confidence region tiap klaster. Klaster yang terpisah dengan jelas dalam plot ini mengonfirmasi bahwa pemisahan oleh DBSCAN memiliki dasar geometrik yang nyata dalam ruang fitur aslinya. Tumpang-tindih antar-elips yang minimal adalah tanda struktur klaster yang kuat.

9. Profil dan Interpretasi Cluster

9.1 Profil Rata-Rata per Cluster

df_profile <- df_result %>%
  filter(Cluster > 0) %>%
  group_by(Cluster) %>%
  summarise(
    N                    = n(),
    Kemiskinan_pct       = round(mean(Tingkat_Penduduk_Miskin), 2),
    RLS_thn              = round(mean(RLS), 2),
    UHH_thn              = round(mean(UHH), 2),
    HLS_thn              = round(mean(HLS), 2),
    IPG_pct              = round(mean(IPG), 2),
    PengeluaranPerKapita = round(mean(PengeluaranPerKapita)),
    PDRB_juta            = round(mean(PDRB)),
    IKK                  = round(mean(IKK), 2),
    Pengeluaran_Rokok    = round(mean(PengeluaranPerkapita_Rokok)),
    .groups = "drop"
  )

kable(as.data.frame(df_profile), row.names = FALSE,
      col.names = c("Cluster", "N", "Kemiskinan (%)", "RLS (thn)",
                    "UHH (thn)", "HLS (thn)", "IPG (%)",
                    "Pengeluaran/Kap", "PDRB (juta)", "IKK",
                    "Rokok/Kap"),
      caption = "Tabel 9. Profil Rata-Rata Indikator per Cluster",
      align = c("c", rep("r", 10))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE, font_size = 12) %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(3, bold = TRUE, color = "white",
              background = spec_color(df_profile$Kemiskinan_pct,
                                      option = "D", direction = -1))

Tabel 9. Profil Rata-Rata Indikator per Cluster
Cluster	N	Kemiskinan (%)	RLS (thn)	UHH (thn)	HLS (thn)	IPG (%)	Pengeluaran/Kap	PDRB (juta)	IKK	Rokok/Kap
1	409	8.96	9.01	71.28	13.39	91.74	12234	51868	99.51	22911
2	5	33.17	5.30	65.93	10.45	72.48	5981	962	216.72	29335

Interpretasi: Kolom “Kemiskinan (%)” disorot dengan gradasi warna, semakin gelap menandakan tingkat kemiskinan yang lebih tinggi. Profil ini memungkinkan perbandingan multidimensi antar-klaster: klaster dengan tingkat kemiskinan tinggi umumnya disertai RLS rendah, UHH rendah, dan pengeluaran per kapita yang lebih rendah, konsisten dengan teori lingkaran kemiskinan.

9.2 Kategori Tingkat Kemiskinan per Cluster

interp <- df_profile %>%
  arrange(desc(Kemiskinan_pct)) %>%
  mutate(Kategori = case_when(
    Kemiskinan_pct >= 20 ~ "Kemiskinan Sangat Tinggi",
    Kemiskinan_pct >= 12 ~ "Kemiskinan Tinggi",
    Kemiskinan_pct >= 7  ~ "Kemiskinan Sedang",
    TRUE                 ~ "Kemiskinan Rendah"
  ))

kable(interp[, c("Cluster", "N", "Kemiskinan_pct", "Kategori")],
      col.names   = c("Cluster", "Jumlah Kab/Kota", "Rata-rata Kemiskinan (%)", "Kategori"),
      caption = "Tabel 10. Kategorisasi Tingkat Kemiskinan per Cluster",
      row.names = FALSE,
      align = c("c", "c", "c", "l")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13) %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(4, bold = TRUE, color = case_when(
    interp$Kategori == "Kemiskinan Sangat Tinggi" ~ "red",
    interp$Kategori == "Kemiskinan Tinggi"        ~ "darkorange",
    interp$Kategori == "Kemiskinan Sedang"        ~ "goldenrod",
    TRUE                                          ~ "darkgreen"
  ))

Tabel 10. Kategorisasi Tingkat Kemiskinan per Cluster
Cluster	Jumlah Kab/Kota	Rata-rata Kemiskinan (%)	Kategori
2	5	33.17	Kemiskinan Sangat Tinggi
1	409	8.96	Kemiskinan Sedang

Interpretasi: Kategorisasi ini memberikan label deskriptif yang mudah dipahami oleh pemangku kebijakan:

Kemiskinan Sangat Tinggi (≥ 20%): Wilayah prioritas utama intervensi, membutuhkan program perlindungan sosial, peningkatan akses layanan dasar, dan stimulus ekonomi yang masif.
Kemiskinan Tinggi (12–20%): Wilayah yang masih memerlukan perhatian serius, namun sudah memiliki kapasitas dasar yang dapat dikembangkan melalui program pemberdayaan.
Kemiskinan Sedang (7–12%): Wilayah transisi yang perlu didorong untuk melampaui batas kemiskinan melalui peningkatan kualitas SDM dan produktivitas ekonomi.
Kemiskinan Rendah (< 7%): Wilayah relatif maju yang dapat menjadi model dan pusat pertumbuhan regional.

9.3 Distribusi Kemiskinan per Cluster (Boxplot)

p4 <- ggplot(
  df_result %>%
    filter(Cluster > 0) %>%
    mutate(Cluster = as.factor(Cluster)),
  aes(x = Cluster, y = Tingkat_Penduduk_Miskin, fill = Cluster)
) +
  geom_boxplot(alpha = 0.75, outlier.shape = 21) +
  scale_fill_brewer(palette = "Set1") +
  labs(title    = "Gambar 4. Distribusi Tingkat Kemiskinan per Cluster",
       subtitle = sprintf("Metode: %s | Silhouette = %.4f", best_method, best_sil_val),
       x = "Cluster", y = "Tingkat Penduduk Miskin (%)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

print(p4)

Gambar 4. Distribusi Tingkat Kemiskinan per Cluster

Interpretasi: Boxplot menyajikan distribusi tingkat kemiskinan dalam tiap klaster secara lebih rinci. Lebar kotak (IQR) menggambarkan heterogenitas internal klaster, klaster dengan kotak sempit berarti anggotanya homogen; kotak yang lebar mengindikasikan keragaman yang tinggi di dalam klaster tersebut. Keberadaan outlier (titik di luar whisker) dalam suatu klaster menunjukkan kabupaten/kota dengan profil kemiskinan yang menyimpang dari rata-rata kelompoknya, layak mendapat perhatian analisis lebih lanjut.

10. Kesimpulan

final_tbl <- comparison_tbl[, c("Metode", "n_Cluster", "Silhouette", "Dunn", "CH_Index")]

kable(final_tbl,
      col.names = c("Metode", "Jumlah Klaster", "Silhouette", "Dunn", "CH Index"),
      caption   = "Tabel 11. Ringkasan Akhir Perbandingan Metode Clustering",
      align     = c("l", "c", "c", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"),
                full_width = TRUE, font_size = 13) %>%
  row_spec(best_idx, bold = TRUE, background = "#d4edda") %>%
  column_spec(1, bold = TRUE)

Tabel 11. Ringkasan Akhir Perbandingan Metode Clustering
Metode	Jumlah Klaster	Silhouette	Dunn	CH Index
K-Means	3	0.6876	0.2293	378.63
K-Median	2	0.6449	0.1766	280.57
DBSCAN	2	0.7714	1.5227	168.86
Mean Shift	NA	NA	NA	NA
Fuzzy C-Means	2	0.6764	0.1253	371.52

Berdasarkan hasil analisis komparatif lima metode clustering terhadap data indikator kemiskinan 514 kabupaten/kota di Indonesia tahun 2024, diperoleh kesimpulan sebagai berikut:

Metode terbaik adalah DBSCAN dengan Silhouette Score tertinggi sebesar 0.7714, yang menghasilkan 2 klaster dengan pemisahan paling jelas di antara semua metode yang diuji.
Preprocessing tiga tahap (Winsorize P5–P95 –> Log1p –> RobustScaler) terbukti efektif dalam meningkatkan kualitas clustering pada data yang heterogen dan mengandung outlier seperti indikator kemiskinan antarwilayah.
Profil klaster mengungkap pola yang konsisten: wilayah dengan tingkat kemiskinan tinggi umumnya dicirikan oleh rendahnya RLS, UHH, pengeluaran per kapita, dan PDRB, hal ini mengindikasikan bahwa kemiskinan bersifat multidimensi dan saling terkait antar-indikator.
Implikasi kebijakan: pengelompokan berbasis klaster ini dapat menjadi dasar desain kebijakan yang lebih tepat sasaran. Intervensi tidak perlu seragam, namun setiap klaster membutuhkan pendekatan yang disesuaikan dengan karakteristik dominannya, misalnya fokus pada peningkatan pendidikan untuk klaster dengan RLS rendah, atau stimulus ekonomi lokal untuk klaster dengan PDRB dan pengeluaran per kapita yang rendah.