Workshop R Studio
BPS Babel
2025-10-24
1. ⚙️ Setup Awal dan Pengambilan Data
# === 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, unionlibrary(purrr) # Iterasi dan pemetaan fungsi##
## Attaching package: 'purrr'## The following object is masked from 'package:jsonlite':
##
## flattenlibrary(tidyr) # Ubah bentuk data
library(stringr) # Olah teks dan regex
library(ggplot2) # Visualisasi data
library(viridis) # Palet warna modern## Loading required package: viridisLitelibrary(factoextra) # Visualisasi dan analisis clustering## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(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 TRUElibrary(patchwork) # Gabungkan beberapa plot
library(grid) # Atur tata letak grafis2. 🛰️ Pengambilan Data API BPS
# === 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/e8629b07af3f4067624b0515d06a1ad0",
'Belitung' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1902/var/35/th/110/key/e8629b07af3f4067624b0515d06a1ad0",
'Bangka Barat' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1903/var/30/th/110/key/e8629b07af3f4067624b0515d06a1ad0",
'Bangka Tengah' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1904/var/46/th/110/key/e8629b07af3f4067624b0515d06a1ad0",
'Bangka Selatan' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1905/var/8/th/110/key/e8629b07af3f4067624b0515d06a1ad0",
'Belitung Timur' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1906/var/58/th/110/key/e8629b07af3f4067624b0515d06a1ad0",
'Kota Pangkal Pinang' = "https://webapi.bps.go.id/v1/api/list/model/data/lang/ind/domain/1971/var/48/th/110/key/e8629b07af3f4067624b0515d06a1ad0"
)
# === 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 Pinang3. 🧹 Parsing dan Pembersihan Data
# === 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))4. 🔢 Penentuan K Optimal dan Visualisasi K Optimal
# === 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): 4cat("📊 K Optimal (Silhouette):", k_silhouette, "\n")## 📊 K Optimal (Silhouette): 2cat("📊 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"
)
5. 🧩 Klasterisasi K-Means
# === 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]>, …6. Visualisasi Hasil K-Means
# === 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 offmap <- st_read(dsn = "C:/Users/BPS/Downloads/babel_shp/babel_shp/", quiet = TRUE) %>%
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))
💡 5. Kesimpulan Umum
Kota Pangkalpinang membentuk klaster tersendiri karena peran ekonominya sangat berbeda — berfokus pada jasa, administrasi, dan sektor pendidikan.
Enam kabupaten lain memiliki struktur ekonomi yang relatif homogen dan masih bergantung pada sektor ekstraktif dan pertanian.
Implikasi kebijakan:
Diversifikasi ekonomi kabupaten perlu ditingkatkan agar tidak bergantung pada pertambangan.
Penguatan sektor jasa di daerah lain bisa menjadi strategi mengurangi kesenjangan ekonomi antarwilayah.
Pangkalpinang dapat berperan sebagai pusat pertumbuhan regional (growth pole) untuk mendorong transformasi ekonomi kabupaten sekitarnya.