# === Library utama ===
library(jsonlite) # Baca data JSON dari API
library(dplyr) # Manipulasi data
##
## 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(purrr) # Iterasi dan pemetaan fungsi
##
## Attaching package: 'purrr'
## The following object is masked from 'package:jsonlite':
##
## flatten
library(tidyr) # Ubah bentuk data
library(stringr) # Olah teks dan regex
library(ggplot2) # Visualisasi data
library(viridis) # Palet warna modern
## Loading required package: viridisLite
library(factoextra) # Visualisasi dan analisis clustering
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(cluster) # Analisis cluster lanjut
library(sf) # Data spasial (shapefile)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(patchwork) # Gabungkan beberapa plot
library(grid) # Atur tata letak grafis
# === Daftar API PDRB per daerah (contoh tahun 2023/2024) ===
# Setiap elemen berisi URL API BPS untuk masing-masing kabupaten/kota di Provinsi Babel.
urls <- list(
'Bangka' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1901/var/45/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Belitung' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1902/var/35/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Bangka Barat' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1903/var/30/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Bangka Tengah' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1904/var/46/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Bangka Selatan' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1905/var/8/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Belitung Timur' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1906/var/58/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d",
'Kota Pangkal Pinang' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1971/var/48/th/110/key/dec6ad5937bd46f1da5da3dbbf15cd4d"
)
# === Fungsi untuk mengambil data JSON dari API BPS ===
# Fungsi ini:
# - Menampilkan nama daerah yang sedang diambil
# - Mengambil data dari URL menggunakan fromJSON()
# - Menangani error jika koneksi gagal atau format salah
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)
}
# === Mengambil semua data dari daftar URL ===
# Menggunakan imap() dari purrr agar bisa iterasi dengan nama daerah
# Hasil akhirnya berupa list berisi JSON mentah untuk tiap kabupaten/kota
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
# === Fungsi untuk parsing data PDRB per kabupaten/kota ===
# Tujuan:
# - Mengambil dan menyusun data PDRB dari hasil JSON API BPS
# - Menghubungkan kode variabel dengan label lapangan usaha
# - Menghasilkan data frame: daerah, lapangan_usaha, nilai
parse_pdrb_daerah <- function(res, daerah) {
# Abaikan jika data kosong
if (is.null(res) || is.null(res$datacontent)) return(NULL)
# Ambil komponen data dari struktur JSON
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)
# Buat kombinasi unik key berdasarkan struktur BPS API
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))
# Gabungkan nilai data dengan label variabel
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 untuk semua kabupaten/kota ===
# Fungsi imap_dfr() akan menjalankan parse_pdrb_daerah() untuk setiap elemen raw_list
# dan langsung menggabungkannya menjadi satu data frame besar.
pdrb_kabkota <- imap_dfr(raw_list, parse_pdrb_daerah)
# === Fungsi pembersihan label lapangan usaha ===
# Tujuan:
# - Membersihkan teks label dari karakter HTML
# - Menstandarkan format teks (huruf kecil, tanpa spasi berlebih)
# - Mengelompokkan ke dalam kategori lapangan usaha resmi (A–U)
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_
)
}
# === Membersihkan dan mengelompokkan data PDRB ===
# Langkah-langkah:
# 1. Bersihkan label menjadi kategori resmi
# 2. Hilangkan data tanpa kategori (NA)
# 3. Agregasi nilai per kategori per daerah
# 4. Ubah format ke bentuk wide (kolom = kategori)
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)
# === Penyesuaian satuan data tertentu ===
# Contoh: untuk Bangka Tengah, satuan nilainya 1000 kali lebih besar, jadi dibagi 1000
pdrb_clean <- pdrb_clean %>%
mutate(across(
where(is.numeric),
~ if_else(daerah == "Bangka Tengah", .x / 1000, .x)
))
# === Standardisasi (scaling) data untuk analisis cluster ===
# Tujuan: menormalkan nilai numerik agar tidak bias terhadap skala besar/kecil
# Setiap variabel diubah menjadi skala (mean = 0, sd = 1)
pdrb_scaled <- pdrb_clean %>%
group_by(daerah) %>%
ungroup() %>%
mutate(across(where(is.numeric), scale))
# === Standarisasi data ===
# Skala semua variabel numerik agar setara (mean = 0, sd = 1)
data_scaled <- scale(select(pdrb_clean, -daerah))
# === Tentukan jumlah cluster maksimum ===
# Batasi sampai 10 atau sesuai banyaknya observasi - 1
max_k <- min(10, nrow(data_scaled) - 1)
# --- (a) Elbow Method ---
# Cari titik "tekukan" grafik WSS (Within Sum of Squares)
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 # fallback
# --- (b) Silhouette Method ---
# Mengukur seberapa baik tiap observasi cocok dengan clusternya
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 ---
# Bandingkan variasi antar cluster dengan data acak
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 dari ketiga metode ===
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 jumlah cluster final ===
# Urutan prioritas: Silhouette → Gap → Elbow
k_optimal <- if (!is.na(k_silhouette)) k_silhouette else if (!is.na(k_gap)) k_gap else k_elbow
# === Visualisasi 3 Metode Penentuan Jumlah Cluster ===
# --- (1) Elbow Method ---
# Menampilkan grafik WSS dengan garis vertikal di titik 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)
# --- (2) Silhouette Method ---
# Menampilkan nilai rata-rata silhouette untuk tiap jumlah cluster
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)
# --- (3) Gap Statistic Method ---
# Menampilkan perbandingan variasi antar cluster dengan data acak
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)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# --- Gabungkan ketiga plot dalam satu tampilan ---
(p1 | p2 | p3) +
plot_annotation(
title = "Penentuan Jumlah Cluster Optimal",
subtitle = "Metode: Elbow, Silhouette, dan Gap Statistic"
)
# === Langkah Klasterisasi K-Means ===
# Menetapkan nilai acak awal agar hasil dapat direplikasi
set.seed(123)
# --- Proses K-Means ---
# data_scaled : data yang sudah distandarisasi
# centers : jumlah cluster optimal (hasil dari tahap sebelumnya)
# nstart : jumlah percobaan awal acak untuk mencari hasil terbaik
km <- kmeans(data_scaled, centers = k_optimal, nstart = 25)
# --- Menambahkan hasil cluster ke dataset utama ---
# Menyimpan nomor cluster tiap daerah ke dalam kolom baru "cluster"
pdrb_scaled$cluster <- km$cluster
# --- Menampilkan hasil ---
# Data sudah berisi kolom tambahan 'cluster' yang menunjukkan kelompok tiap daerah
pdrb_scaled
## # A tibble: 7 × 20
## daerah A. Pertanian, Kehuta…¹ B. Pertambangan dan …² C. Industri Pengolah…³
## <chr> <dbl> <dbl> <dbl>
## 1 Bangka 0.758 0.775 0.373
## 2 Bangka B… -0.164 0.376 1.99
## 3 Bangka S… 1.38 1.16 -1.05
## 4 Bangka T… -1.03 0.366 -0.189
## 5 Belitung 0.568 -0.754 -0.664
## 6 Belitung… -0.0522 -0.133 -0.549
## 7 Kota Pan… -1.46 -1.79 0.0936
## # ℹ abbreviated names: ¹`A. Pertanian, Kehutanan, dan Perikanan`[,1],
## # ²`B. Pertambangan dan Penggalian`[,1], ³`C. Industri Pengolahan`[,1]
## # ℹ 16 more variables: `D. Pengadaan Listrik dan Gas` <dbl[,1]>,
## # `E. Pengadaan Air, Pengelolaan Sampah, Limbah, dan Daur Ulang` <dbl[,1]>,
## # `F. Konstruksi` <dbl[,1]>,
## # `G. Perdagangan Besar dan Eceran; Reparasi Mobil dan Sepeda Motor` <dbl[,1]>,
## # `H. Transportasi dan Pergudangan` <dbl[,1]>, …
# === Menambahkan hasil cluster ke data asli ===
pdrb_clean$cluster <- km$cluster
# === Peta Klasterisasi Ekonomi Bangka Belitung ===
sf_use_s2(FALSE) # Matikan validasi geometri agar peta bisa dibaca tanpa error
## Spherical geometry (s2) switched off
map <- st_read(dsn = "C:/Users/LOQ/Downloads/babel_shp/babel_shp", quiet = TRUE) %>% # params$shp_path -> sesuaikan dengan tempat folder peta
rename(daerah = WADMKK) %>%
mutate(daerah = str_to_title(daerah))
# Gabungkan data peta dengan hasil klasterisasi
map_clustered <- map %>%
left_join(pdrb_clean %>% mutate(cluster = factor(cluster)), by = "daerah")
# Visualisasi peta
ggplot(map_clustered) +
geom_sf(aes(fill = cluster), color = "white", size = 0.3) +
labs(title = "Peta Klasterisasi Ekonomi Bangka Belitung") +
theme_void()
# === Visualisasi PCA (proyeksi 2D hasil klasterisasi) ===
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
# === Identifikasi Faktor Dominan Tiap Cluster ===
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) %>%
ungroup()
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(where(is.numeric), mean, na.rm = TRUE)`.
## ℹ In group 1: `cluster = 1`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
# === Visualisasi sektor dominan ===
# Konfigurasi warna pada plot
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))
