📍 Pendahuluan

Analisis ini bertujuan untuk: - Mengelompokkan kabupaten/kota di Provinsi Kepulauan Bangka Belitung berdasarkan struktur PDRB menurut lapangan usaha. - Menentukan jumlah klaster optimal dengan Elbow, Silhouette, dan Gap Statistic. - Menafsirkan hasil klasterisasi untuk memahami struktur ekonomi wilayah.

⚙️ Setup Awal dan Pengambilan Data

# === Library utama ===
library(jsonlite)   # untuk membaca JSON dari API BPS
library(dplyr)      # manipulasi data (filter, mutate, summarise, dll)
library(purrr)      # iterasi fungsional (map/imap)
library(tidyr)      # reshaping data (pivot_wider, pivot_longer)
library(stringr)    # pembersihan teks (regex)
library(ggplot2)    # visualisasi data
library(viridis)    # palet warna modern
library(factoextra) # fungsi K-Means dan visualisasi cluster
library(cluster)    # Gap Statistic
library(sf)         # data spasial shapefile
library(patchwork)  # menggabungkan plot
library(grid)

🛰️ Pengambilan Data API BPS

# Daftar API PDRB per daerah (contoh tahun 2023/2024)
urls <- list(
  'Bangka'              = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1901/var/45/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Belitung'            = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1902/var/35/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Bangka Barat'        = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1903/var/30/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Bangka Tengah'       = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1904/var/46/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Bangka Selatan'      = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1905/var/8/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Belitung Timur'      = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1906/var/58/th/124/key/e8629b07af3f4067624b0515d06a1ad0",
  'Kota Pangkal Pinang' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1971/var/48/th/124/key/e8629b07af3f4067624b0515d06a1ad0"
)

# Fungsi ambil JSON API
get_raw_bps <- function(url, daerah) {
  cat("📦 Mengambil data:", daerah, "\n")
  res <- tryCatch(fromJSON(url), error = function(e) NULL)
  if (is.null(res)) {
    warning(paste("❌ Gagal mengambil data untuk", daerah))
    return(NULL)
  }
  return(res)
}

# Ambil semua data
raw_list <- imap(urls, get_raw_bps)
## 📦 Mengambil data: Bangka 
## 📦 Mengambil data: Belitung 
## 📦 Mengambil data: Bangka Barat 
## 📦 Mengambil data: Bangka Tengah 
## 📦 Mengambil data: Bangka Selatan 
## 📦 Mengambil data: Belitung Timur 
## 📦 Mengambil data: Kota Pangkal Pinang

Penjelasan:

fromJSON() membaca API dan mengonversinya menjadi list di R.

imap() dari purrr digunakan agar bisa mengambil data sekaligus menyimpan nama daerah.

🧹 Parsing dan Pembersihan Data

# Fungsi parsing tiap kabupaten/kota
parse_pdrb_daerah <- function(res, daerah) {
  if (is.null(res) || is.null(res$datacontent)) return(NULL)
  vervar_df <- res$vervar %>% mutate(val = as.character(val))
  var_val <- as.character(res$var$val)
  turvar_val <- as.character(res$turvar$val)
  tahun_val <- as.character(res$tahun$val)
  turtahun_val <- as.character(res$turtahun$val)

  key_combos <- expand.grid(vervar = vervar_df$val, var = var_val,
                            turvar = turvar_val, tahun = tahun_val,
                            turtahun = turtahun_val, stringsAsFactors = FALSE) %>%
    mutate(key = paste0(vervar, var, turvar, tahun, turtahun))
  
  df_raw <- tibble(key = names(res$datacontent), nilai = unlist(res$datacontent))
  df_final <- key_combos %>%
    left_join(df_raw, by = "key") %>%
    left_join(vervar_df, by = c("vervar" = "val")) %>%
    transmute(daerah, lapangan_usaha = label, nilai = as.numeric(nilai))
  return(df_final)
}

# Eksekusi parsing
pdrb_kabkota <- imap_dfr(raw_list, parse_pdrb_daerah)

🧭 Klasifikasi Lapangan Usaha & Pivot ke Format Wide

# Fungsi bersihkan label lapangan usaha
clean_lapangan_usaha <- function(label) {
  label_clean <- label %>%
    str_remove_all("<.*?>") %>%
    str_trim() %>% str_squish() %>% str_to_lower()
  
  case_when(
    str_detect(label_clean, "pertanian|kehutanan|perikanan") ~ "A. Pertanian, Kehutanan, dan Perikanan",
    str_detect(label_clean, "pertambangan") ~ "B. Pertambangan dan Penggalian",
    str_detect(label_clean, "industri pengolahan") ~ "C. Industri Pengolahan",
    str_detect(label_clean, "listrik") ~ "D. Pengadaan Listrik dan Gas",
    str_detect(label_clean, "air|sampah|limbah") ~ "E. Pengadaan Air, Pengelolaan Sampah, Limbah, dan Daur Ulang",
    str_detect(label_clean, "konstruksi") ~ "F. Konstruksi",
    str_detect(label_clean, "perdagangan") ~ "G. Perdagangan Besar dan Eceran; Reparasi Mobil dan Sepeda Motor",
    str_detect(label_clean, "transportasi") ~ "H. Transportasi dan Pergudangan",
    str_detect(label_clean, "akomodasi|makan minum") ~ "I. Penyediaan Akomodasi dan Makan Minum",
    str_detect(label_clean, "informasi|komunikasi") ~ "J. Informasi dan Komunikasi",
    str_detect(label_clean, "keuangan|asuransi") ~ "K. Jasa Keuangan dan Asuransi",
    str_detect(label_clean, "real estat") ~ "L. Real Estat",
    str_detect(label_clean, "perusahaan") ~ "M,N. Jasa Perusahaan",
    str_detect(label_clean, "administrasi|pemerintahan|jaminan sosial") ~ "O. Administrasi Pemerintahan, Pertahanan, dan Jaminan Sosial Wajib",
    str_detect(label_clean, "pendidikan") ~ "P. Jasa Pendidikan",
    str_detect(label_clean, "kesehatan") ~ "Q. Jasa Kesehatan dan Kegiatan Sosial",
    str_detect(label_clean, "lainnya") ~ "R,S,T,U. Jasa Lainnya",
    str_detect(label_clean, "pdrb|produk domestik regional bruto") ~ "Total PDRB",
    TRUE ~ NA_character_
  )
}

pdrb_clean <- pdrb_kabkota %>%
  mutate(kategori = clean_lapangan_usaha(lapangan_usaha)) %>%
  filter(!is.na(kategori)) %>%
  group_by(daerah, kategori) %>%
  summarise(nilai = sum(nilai, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = kategori, values_from = nilai, values_fill = 0)

pdrb_clean <- pdrb_clean %>%
  mutate(across(
    where(is.numeric),
    ~ if_else(daerah == "Bangka Tengah", .x / 1000, .x)
  ))

pdrb_scaled <- pdrb_clean %>%
  group_by(daerah) %>%
  ungroup() %>%
  mutate(across(where(is.numeric), scale))

🔢 Klasterisasi K-Means

# Standarisasi agar tiap variabel setara
data_scaled <- scale(select(pdrb_clean, -daerah))

# Menentukan jumlah cluster optimal
max_k <- min(10, nrow(data_scaled) - 1) # aman untuk data kecil

# --- (a) Elbow Method ---
set.seed(123)
fviz_elbow <- fviz_nbclust(data_scaled, kmeans, method = "wss", k.max = max_k)
elbow_data <- fviz_elbow$data
k_elbow <- which(diff(diff(elbow_data$y)) == min(diff(diff(elbow_data$y)))) + 1
if (length(k_elbow) == 0 || is.na(k_elbow)) k_elbow <- 2

# --- (b) Silhouette Method ---
set.seed(123)
sil_result <- fviz_nbclust(data_scaled, kmeans, method = "silhouette", k.max = max_k)
k_silhouette <- sil_result$data$clusters[which.max(sil_result$data$y)] %>% as.numeric()

# --- (c) Gap Statistic ---
set.seed(123)
gap_stat <- clusGap(data_scaled, FUN = kmeans, nstart = 25, K.max = max_k, B = 100)
k_gap <- maxSE(gap_stat$Tab[, "gap"], gap_stat$Tab[, "SE.sim"], method = "firstSEmax")

# === Ringkasan Hasil ===
cat("📊 K Optimal (Elbow):", k_elbow, "\n")
## 📊 K Optimal (Elbow): 4
cat("📊 K Optimal (Silhouette):", k_silhouette, "\n")
## 📊 K Optimal (Silhouette): 2
cat("📊 K Optimal (Gap Statistic):", k_gap, "\n")
## 📊 K Optimal (Gap Statistic): 1
# Pilih K Final (prioritas: Silhouette → Gap → Elbow)
k_optimal <- if (!is.na(k_silhouette)) k_silhouette else if (!is.na(k_gap)) k_gap else k_elbow
cat("\n✅ Jumlah Cluster Final:", k_optimal, "\n")
## 
## ✅ Jumlah Cluster Final: 2

Penjelasan Fungsi:

fviz_nbclust(…, method=“wss”) → mencari titik “tekukan” (Elbow) antara jumlah cluster & total variasi.

fviz_nbclust(…, method=“silhouette”) → mengukur seberapa baik tiap observasi cocok dalam cluster-nya.

clusGap() → menghitung Gap Statistic (pembeda signifikan antar cluster).

Plot Penentuan K Optimal

p1 <- fviz_nbclust(data_scaled, kmeans, method = "wss", k.max = max_k) +
  geom_vline(xintercept = k_elbow, linetype = "dashed", color = "#E4572E") +
  annotate("text", x = k_elbow, y = max(elbow_data$y), label = paste("K =", k_elbow),
           vjust = -1, color = "#E4572E", size = 3.5, fontface = "bold") +
  ggtitle("Elbow Method (WSS)") +
  theme_minimal(base_size = 11)

p2 <- fviz_nbclust(data_scaled, kmeans, method = "silhouette", k.max = max_k) +
  geom_vline(xintercept = k_silhouette, linetype = "dashed", color = "#4CAF50") +
  annotate("text", x = k_silhouette, y = 0, label = paste("K =", k_silhouette),
           vjust = -1, color = "#4CAF50", size = 3.5, fontface = "bold") +
  ggtitle("Silhouette Method") +
  theme_minimal(base_size = 11)

p3 <- fviz_gap_stat(gap_stat) +
  geom_vline(xintercept = k_gap, linetype = "dashed", color = "yellow") +
  annotate("text", x = k_gap, y = max(gap_stat$Tab[,"gap"]), label = paste("K =", k_gap),
           vjust = -1, color = "yellow", size = 3.5, fontface = "bold") +
  ggtitle("Gap Statistic Method") +
  theme_minimal(base_size = 11)

(p1 | p2 | p3) +
  plot_annotation(
    title = "Penentuan Jumlah Cluster Optimal",
    subtitle = "Metode: Elbow, Silhouette, dan Gap Statistic"
  )

🧩 Hasil Akhir K-Means

set.seed(123)
km <- kmeans(data_scaled, centers = k_optimal, nstart = 25)
pdrb_scaled$cluster <- km$cluster

pdrb_clean$cluster <- km$cluster
fviz_cluster(km, data = data_scaled, geom = "point", ellipse.type = "norm") +
  labs(title = paste("K-Means Clustering (K =", k_optimal, ")"))
## Too few points to calculate an ellipse

pdrb_scaled %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>%
  pivot_longer(cols = -cluster, names_to = "Sektor", values_to = "Rata2") %>%
  group_by(cluster) %>%
  slice_max(Rata2, n = 1, with_ties = FALSE) %>%
  ungroup()
faktor_dominan <- pdrb_scaled %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>%
  pivot_longer(cols = -cluster, names_to = "Sektor", values_to = "Rata2") %>%
  group_by(cluster) %>%
  slice_max(Rata2, n = 1, with_ties = FALSE) %>%
  ungroup()

pastel_colors <- c("#FFB3BA", "#BAFFC9", "#FFF5BA", "#D5BAFF", "#FFE4BA", "#B4E1FF", "#D1BAFF", "#BAFFD5", "#FFD1BA", "#F5FFBA")

ggplot(faktor_dominan, aes(x = factor(cluster), y = Rata2, fill = Sektor)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  geom_text(aes(label = Sektor), vjust = -0.8, size = 3.5, fontface = "bold") +
  scale_fill_manual(values = pastel_colors) +
  labs(
    title = "Faktor Dominan Tiap Cluster",
    x = "Cluster",
    y = "Rata-rata Z-score"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 14)
  )

🗺️ Peta Klasterisasi

sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
map <- st_read(dsn = params$shp_path, quiet = TRUE) %>%
  rename(daerah = WADMKK) %>%
  mutate(daerah = str_to_title(daerah))

map_clustered <- map %>%
  left_join(pdrb_clean %>% mutate(cluster = factor(cluster)), by = "daerah")

ggplot(map_clustered) +
  geom_sf(aes(fill = cluster), color = "white", size = 0.3) +
  labs(title = "Peta Klasterisasi Ekonomi Bangka Belitung") +
  theme_void()

# Buat data frame PCA
pca_res <- prcomp(data_scaled, scale. = FALSE)
pca_df <- data.frame(pca_res$x[, 1:2]) %>%
  mutate(daerah = pdrb_scaled$daerah, cluster = factor(km$cluster))

ggplot(pca_df, aes(x = PC1, y = PC2, color = cluster)) +
  geom_point(size = 3) +
  geom_text(aes(label = daerah), size = 3, vjust = -0.6, check_overlap = TRUE) +
  stat_ellipse(aes(fill = cluster), geom = "polygon", alpha = 0.12) +
  labs(title = paste("PCA - Proyeksi 2D (K =", k_optimal, ")")) +
  theme_minimal()
## Too few points to calculate an ellipse

🧠 Interpretasi Hasil ⚙️ Ringkasan Hasil

Jumlah cluster: 2

Ukuran cluster:

Cluster 1: 1 kabupaten/kota

Cluster 2: 6 kabupaten/kota

Proporsi variasi antar cluster: r round(km\(betweenss / km\)totss * 100, 1)%

→ Artinya, sekitar r round(km\(betweenss / km\)totss * 100, 1)% variasi data dijelaskan oleh dua cluster ini — pemisahan cukup moderat dan jelas terlihat.

🧩 Interpretasi Tiap Cluster 🌆 Cluster 1 — Ekonomi Jasa dan Perkotaan

Ciri:

Z-score sangat tinggi di sektor tersier: perdagangan, jasa keuangan, real estat, pendidikan, kesehatan.

Nilai negatif di sektor primer: pertanian, pertambangan.

Makna:

Mewakili daerah service-oriented dan modern, pusat aktivitas ekonomi jasa dan administrasi.

Kemungkinan besar: Kota Pangkalpinang.

🌾 Cluster 2 — Ekonomi Primer dan Tradisional

Ciri:

Z-score sedikit positif di sektor pertanian dan pertambangan.

Negatif pada jasa dan industri tersier.

Makna:

Menggambarkan kabupaten dengan basis ekonomi primer dan sekunder (agraris dan tambang).

Kemungkinan besar: Kabupaten Bangka, Bangka Tengah, Bangka Barat, Bangka Selatan, Belitung, Belitung Timur.

🧭 Kesimpulan Umum Aspek Cluster 1 Cluster 2 Jumlah daerah 1 6 Tipe ekonomi Jasa dan perkotaan Primer & sekunder (agraris-tambang) Contoh daerah Kota Pangkalpinang Kabupaten Bangka, Belitung, dkk Ciri ekonomi Terdiversifikasi, modern, tersier dominan Tradisional, berbasis SDA Arah kebijakan Penguatan sektor jasa, digitalisasi Diversifikasi ke industri & jasa

💡 Kesimpulan ringkas: Hasil K-Means menunjukkan dualisme struktur ekonomi di Provinsi Kep. Bangka Belitung — kota utama (Pangkalpinang) berorientasi jasa, sedangkan kabupaten lainnya masih berbasis sumber daya alam.

💾 Ekspor Hasil write.csv(pdrb_clean, file.path(params$out_dir, “hasil_cluster.csv”), row.names = FALSE)

✅ Analisis selesai. Silakan buka output HTML untuk eksplorasi visualisasi dan hasil interpretasi.