Analisis cluster merupakan salah satu metode statistika multivariat yang digunakan untuk mengelompokkan objek ke dalam beberapa kelompok berdasarkan tingkat kemiripan karakteristik yang dimiliki. Pada penelitian ini, analisis cluster digunakan untuk mengelompokkan kecamatan di Kabupaten Pekalongan berdasarkan indikator pertanian dan peternakan tahun 2024.
Metode clustering yang digunakan adalah Partitioning Around Medoids (PAM). Metode ini dipilih karena lebih robust terhadap outlier dibandingkan metode k-means. Penentuan jumlah cluster terbaik dilakukan dengan menggunakan average silhouette width.
Data yang digunakan dalam analisis ini merupakan data sekunder tingkat kecamatan di Kabupaten Pekalongan tahun 2024 yang diperoleh dari website resmi BPS Kabupaten Pekalongan (https://pekalongankab.bps.go.id). Dataset ini memuat indikator sektor pertanian dan peternakan yang digunakan sebagai dasar dalam pengelompokan kecamatan.
Variabel yang digunakan dalam penelitian ini terdiri atas satu variabel identitas dan lima variabel numerik, yaitu:
Variabel yang digunakan dalam proses clustering adalah seluruh variabel numerik, yaitu padi_prod, sapi, kambing, ayam, dan itik. Variabel kecamatan digunakan sebagai penanda untuk masing-masing objek.
| kecamatan | padi_prod | sapi | kambing | ayam | itik |
|---|---|---|---|---|---|
| Bojong | 23359 | 97 | 1436 | 1162902 | 17781 |
| Buaran | 3104 | 20 | 1635 | 6905 | 1141 |
| Doro | 7630 | 243 | 3636 | 1345364 | 15728 |
| Kajen | 24933 | 252 | 3284 | 1777428 | 19482 |
| Kandangserang | 12540 | 1876 | 5305 | 80873 | 490 |
| Karanganyar | 10670 | 416 | 3038 | 1162902 | 17781 |
| Karangdadap | 8242 | 244 | 1937 | 64395 | 9956 |
| Kedungwuni | 7855 | 182 | 2098 | 72319 | 11019 |
| Kesesi | 36901 | 243 | 3299 | 2364226 | 29271 |
| Lebakbarang | 3537 | 1661 | 5395 | 24923 | 925 |
| Paninggaran | 9150 | 4305 | 5446 | 471650 | 931 |
| Petungkriyono | 2768 | 2873 | 4537 | 10687 | 230 |
| Siwalan | 10046 | 106 | 1080 | 303458 | 14797 |
| Sragi | 23778 | 43 | 1148 | 1342046 | 14525 |
| Talun | 7528 | 436 | 3106 | 1858363 | 4852 |
| Tirto | 3794 | 273 | 1458 | 102309 | 7870 |
| Wiradesa | 4098 | 259 | 1206 | 20732 | 2035 |
| Wonokerto | 1713 | 235 | 1294 | 4725 | 1042 |
| Wonopringgo | 4314 | 59 | 1717 | 1168898 | 12609 |
Data awal menunjukkan objek pengamatan berupa kecamatan di Kabupaten Pekalongan beserta beberapa indikator pertanian dan peternakan. Data ini selanjutnya digunakan sebagai dasar dalam proses standarisasi dan analisis cluster.
variabel_tbl <- data.frame(
Nama_Variabel = c("kecamatan", "padi_prod", "sapi", "kambing", "ayam", "itik"),
Jenis = c("Identitas", "Numerik", "Numerik", "Numerik", "Numerik", "Numerik"),
Keterangan = c(
"Nama kecamatan",
"Produksi padi per kecamatan",
"Jumlah ternak sapi",
"Jumlah ternak kambing",
"Jumlah ternak ayam",
"Jumlah ternak itik"
)
)
kable(variabel_tbl, caption = "Daftar Variabel Penelitian")| Nama_Variabel | Jenis | Keterangan |
|---|---|---|
| kecamatan | Identitas | Nama kecamatan |
| padi_prod | Numerik | Produksi padi per kecamatan |
| sapi | Numerik | Jumlah ternak sapi |
| kambing | Numerik | Jumlah ternak kambing |
| ayam | Numerik | Jumlah ternak ayam |
| itik | Numerik | Jumlah ternak itik |
Standarisasi dilakukan menggunakan metode z-score agar seluruh variabel numerik memiliki skala yang seragam.
| padi_prod | sapi | kambing | ayam | itik |
|---|---|---|---|---|
| 1.3097 | -0.5475 | -0.8513 | 0.5931 | 0.9746 |
| -0.8093 | -0.6144 | -0.7214 | -0.8956 | -1.0086 |
| -0.3358 | -0.4208 | 0.5852 | 0.8281 | 0.7299 |
| 1.4744 | -0.4129 | 0.3554 | 1.3845 | 1.1773 |
| 0.1779 | 0.9973 | 1.6750 | -0.8004 | -1.0862 |
| -0.0178 | -0.2705 | 0.1948 | 0.5931 | 0.9746 |
| -0.2718 | -0.4199 | -0.5242 | -0.8216 | 0.0420 |
| -0.3123 | -0.4737 | -0.4190 | -0.8114 | 0.1687 |
| 2.7265 | -0.4208 | 0.3652 | 2.1402 | 2.3440 |
| -0.7640 | 0.8106 | 1.7338 | -0.8724 | -1.0343 |
| -0.1768 | 3.1066 | 1.7671 | -0.2971 | -1.0336 |
| -0.8445 | 1.8631 | 1.1736 | -0.8908 | -1.1171 |
| -0.0831 | -0.5397 | -1.0838 | -0.5137 | 0.6190 |
| 1.3536 | -0.5944 | -1.0394 | 0.8238 | 0.5866 |
| -0.3465 | -0.2532 | 0.2392 | 1.4887 | -0.5663 |
| -0.7371 | -0.3947 | -0.8369 | -0.7728 | -0.2066 |
| -0.7053 | -0.4069 | -1.0015 | -0.8778 | -0.9020 |
| -0.9549 | -0.4277 | -0.9440 | -0.8984 | -1.0204 |
| -0.6827 | -0.5805 | -0.6678 | 0.6008 | 0.3582 |
Hasil standarisasi menunjukkan bahwa seluruh variabel numerik telah diubah ke dalam skala yang sama dengan rata-rata mendekati 0 dan simpangan baku mendekati 1. Proses ini diperlukan agar tidak ada variabel yang mendominasi perhitungan jarak dalam analisis cluster.
cek_z <- data.frame(
Variabel = num_cols,
Mean_Z = round(colMeans(Xz), 4),
SD_Z = round(apply(Xz, 2, sd), 4)
)
kable(cek_z, caption = "Pemeriksaan Mean dan Standar Deviasi Hasil Z-Score")| Variabel | Mean_Z | SD_Z | |
|---|---|---|---|
| padi_prod | padi_prod | 0 | 1 |
| sapi | sapi | 0 | 1 |
| kambing | kambing | 0 | 1 |
| ayam | ayam | 0 | 1 |
| itik | itik | 0 | 1 |
Berdasarkan hasil pemeriksaan, nilai rata-rata masing-masing variabel hasil standarisasi berada di sekitar 0 dan simpangan baku berada di sekitar 1. Hal ini menunjukkan bahwa proses standarisasi z-score telah dilakukan dengan benar.
Outlier multivariat dideteksi menggunakan jarak Mahalanobis, kemudian dibandingkan dengan cutoff distribusi chi-square.
D2 <- mahalanobis(Xz, center = colMeans(Xz), cov = cov(Xz))
p <- ncol(Xz)
chi_cut <- qchisq(1 - alpha_md, df = p)
out_md <- D2 > chi_cut
outlier_tbl <- df %>%
transmute(
!!id_col := .data[[id_col]],
D2 = D2,
outlier_md = out_md
) %>%
arrange(desc(D2))## Alpha = 0.05
## Derajat bebas = 5
## Cutoff Chi-square = 11.0705
## Jumlah outlier = 1
Deteksi outlier multivariat dilakukan menggunakan jarak Mahalanobis dengan pembanding distribusi chi-square. Jumlah observasi yang terdeteksi sebagai outlier menunjukkan adanya kecamatan yang memiliki karakteristik cukup berbeda dibandingkan kecamatan lainnya.
| kecamatan | D2 | outlier_md |
|---|---|---|
| Paninggaran | 12.233170 | TRUE |
| Talun | 10.637597 | FALSE |
| Kesesi | 7.873484 | FALSE |
| Kandangserang | 7.694313 | FALSE |
| Lebakbarang | 6.167878 | FALSE |
| Doro | 5.362307 | FALSE |
| Sragi | 4.997877 | FALSE |
| Wonopringgo | 4.187205 | FALSE |
| Petungkriyono | 4.006916 | FALSE |
| Siwalan | 3.473339 | FALSE |
Tabel di atas menampilkan kecamatan dengan nilai jarak Mahalanobis terbesar. Semakin besar nilai tersebut, semakin jauh karakteristik suatu kecamatan dari pola umum data secara multivariat.
| kecamatan | D2 | outlier_md |
|---|---|---|
| Paninggaran | 12.23317 | TRUE |
Berdasarkan hasil deteksi outlier multivariat, terdapat 1 kecamatan yang teridentifikasi sebagai outlier, yaitu Paninggaran dengan nilai jarak Mahalanobis sebesar 12.23317. Meskipun demikian, observasi tersebut tetap dipertahankan dalam analisis karena metode K-Medoids (PAM) bersifat lebih robust terhadap outlier, sehingga keberadaan outlier tidak terlalu memengaruhi hasil pengelompokan.
Uji KMO digunakan untuk melihat kecukupan sampel dalam analisis multivariat.
kmo_res <- KMO(df %>% select(all_of(num_cols)))
kmo_value <- data.frame(KMO_MSA = kmo_res$MSA)
kable(kmo_value, caption = "Nilai KMO")| KMO_MSA |
|---|
| 0.6493824 |
Nilai KMO sebesar 0.6494 menunjukkan bahwa data memiliki kecukupan sampel yang cukup untuk analisis multivariat. Oleh karena itu, variabel yang digunakan dapat dinilai layak untuk proses analisis cluster.
Multikolinearitas diuji menggunakan Variance Inflation Factor (VIF) secara manual.
vif_tbl <- lapply(num_cols, function(y_var) {
x_vars <- setdiff(num_cols, y_var)
fit <- lm(as.formula(paste(y_var, "~", paste(x_vars, collapse = " + "))), data = df)
r2 <- summary(fit)$r.squared
data.frame(
y_target = y_var,
R2 = r2,
VIF = 1 / (1 - r2)
)
}) %>%
bind_rows() %>%
arrange(desc(VIF)) %>%
mutate(
R2 = round(R2, 4),
VIF = round(VIF, 4)
)| y_target | R2 | VIF |
|---|---|---|
| itik | 0.7845 | 4.6410 |
| sapi | 0.7404 | 3.8525 |
| padi_prod | 0.6990 | 3.3223 |
| kambing | 0.6826 | 3.1508 |
| ayam | 0.6637 | 2.9733 |
Berdasarkan hasil uji multikolinearitas, seluruh variabel memiliki nilai VIF kurang dari 5, yaitu berkisar antara 2.9733 hingga 4.6410. Hal ini menunjukkan bahwa tidak terdapat masalah multikolinearitas yang kuat antar variabel, sehingga seluruh variabel masih layak dipertahankan dalam analisis cluster.
Data yang digunakan pada proses clustering adalah data hasil standarisasi z-score.
X_final <- Xz
D <- dist(X_final, method = dist_method)
Dm <- as.matrix(D)
n <- nrow(X_final)
k_max_eff <- min(k_max, n - 1)
ks <- k_min:k_max_eff
pam_models <- setNames(vector("list", length(ks)), as.character(ks))
metrics <- data.frame(k = ks, avg_silhouette = NA_real_)
for (i in seq_along(ks)) {
k <- ks[i]
fit <- pam(D, k = k, diss = TRUE)
pam_models[[as.character(k)]] <- fit
metrics$avg_silhouette[i] <- mean(silhouette(fit$clustering, D)[, 3])
}
best_k <- metrics$k[which.max(metrics$avg_silhouette)]
best_si <- max(metrics$avg_silhouette)
best_model <- pam_models[[as.character(best_k)]]
df_best <- df %>% mutate(cluster = factor(best_model$clustering))| k | avg_silhouette |
|---|---|
| 2 | 0.4650 |
| 3 | 0.4656 |
| 4 | 0.4430 |
| 5 | 0.3811 |
| 6 | 0.3935 |
| 7 | 0.3521 |
| 8 | 0.3481 |
| 9 | 0.3424 |
| 10 | 0.3518 |
## Jumlah cluster terbaik = 3
## Nilai average silhouette terbaik = 0.4656
Berdasarkan nilai average silhouette, jumlah cluster terbaik yang diperoleh adalah 3 cluster dengan nilai silhouette sebesar 0.4656. Nilai ini menunjukkan bahwa hasil pengelompokan sudah cukup baik, sehingga tiga cluster dianggap sebagai jumlah yang paling optimal untuk mewakili pola kemiripan antar kecamatan.
ggplot(metrics, aes(k, avg_silhouette)) +
geom_line() +
geom_point() +
geom_vline(xintercept = best_k, linetype = "dashed") +
geom_hline(yintercept = best_si, linetype = "dashed") +
annotate(
"text",
x = best_k,
y = best_si,
label = paste0("Best k = ", best_k, "\nSilhouette = ", round(best_si, 3)),
vjust = -1,
hjust = 0.5
) +
theme_minimal() +
labs(
title = "Average Silhouette Width untuk Berbagai Nilai k",
x = "Jumlah Cluster (k)",
y = "Average Silhouette Width"
)anggota_cluster <- df_best %>%
select(all_of(id_col), cluster) %>%
arrange(cluster)
kable(anggota_cluster, caption = "Anggota Cluster")| kecamatan | cluster |
|---|---|
| Bojong | 1 |
| Doro | 1 |
| Kajen | 1 |
| Karanganyar | 1 |
| Kesesi | 1 |
| Sragi | 1 |
| Talun | 1 |
| Buaran | 2 |
| Karangdadap | 2 |
| Kedungwuni | 2 |
| Siwalan | 2 |
| Tirto | 2 |
| Wiradesa | 2 |
| Wonokerto | 2 |
| Wonopringgo | 2 |
| Kandangserang | 3 |
| Lebakbarang | 3 |
| Paninggaran | 3 |
| Petungkriyono | 3 |
Berdasarkan hasil pengelompokan, diperoleh 3 cluster kecamatan di Kabupaten Pekalongan. Cluster 1 beranggotakan 7 kecamatan, cluster 2 beranggotakan 8 kecamatan, dan cluster 3 beranggotakan 4 kecamatan. Hasil ini menunjukkan adanya perbedaan karakteristik antar kelompok kecamatan berdasarkan variabel pertanian dan peternakan yang digunakan.
medoid_tbl <- data.frame(
cluster = seq_along(best_model$id.med),
medoid_row_index = best_model$id.med,
medoid_kecamatan = df_best[[id_col]][best_model$id.med]
)
kable(medoid_tbl, caption = "Medoid Tiap Cluster")| cluster | medoid_row_index | medoid_kecamatan |
|---|---|---|
| 1 | 4 | Kajen |
| 2 | 16 | Tirto |
| 3 | 10 | Lebakbarang |
Medoid merupakan objek yang paling mewakili masing-masing cluster. Berdasarkan hasil analisis, medoid untuk Cluster 1 adalah Kajen, medoid untuk Cluster 2 adalah Tirto, dan medoid untuk Cluster 3 adalah Lebakbarang. Dengan demikian, ketiga kecamatan tersebut dapat dianggap sebagai kecamatan yang paling representatif dalam cluster masing-masing.
medoid_idx <- best_model$id.med
dist_to_medoids <- Dm[, medoid_idx, drop = FALSE]
colnames(dist_to_medoids) <- paste0("Medoid_", seq_along(medoid_idx))
tab_jarak <- data.frame(
Kecamatan = df[[id_col]],
cluster = best_model$clustering,
dist_to_medoids
)
own_medoid <- medoid_idx[tab_jarak$cluster]
tab_jarak$dist_to_own_medoid <- Dm[cbind(1:nrow(df), own_medoid)]
tab_jarak <- tab_jarak[order(tab_jarak$cluster, tab_jarak$dist_to_own_medoid), ]tab_jarak_tampil <- tab_jarak
tab_jarak_tampil[ , sapply(tab_jarak_tampil, is.numeric)] <- round(tab_jarak_tampil[ , sapply(tab_jarak_tampil, is.numeric)], 4)
kable(tab_jarak_tampil, caption = "Jarak Setiap Anggota terhadap Medoid")| Kecamatan | cluster | Medoid_1 | Medoid_2 | Medoid_3 | dist_to_own_medoid | |
|---|---|---|---|---|---|---|
| 4 | Kajen | 1 | 0.0000 | 6.9633 | 9.3090 | 0.0000 |
| 1 | Bojong | 1 | 2.5001 | 4.7611 | 9.4915 | 2.5001 |
| 6 | Karanganyar | 1 | 2.7894 | 4.4223 | 6.8409 | 2.7894 |
| 14 | Sragi | 1 | 2.8486 | 4.8826 | 9.6129 | 2.8486 |
| 3 | Doro | 1 | 3.0517 | 4.3869 | 6.2729 | 3.0517 |
| 9 | Kesesi | 1 | 3.1920 | 10.1553 | 12.4814 | 3.1920 |
| 15 | Talun | 1 | 3.9448 | 4.2295 | 5.8051 | 3.9448 |
| 16 | Tirto | 2 | 6.9633 | 0.0000 | 4.7303 | 0.0000 |
| 17 | Wiradesa | 2 | 7.8844 | 1.0090 | 4.1492 | 1.0090 |
| 7 | Karangdadap | 2 | 5.9741 | 1.1007 | 5.1079 | 1.1007 |
| 18 | Wonokerto | 2 | 8.2241 | 1.2972 | 4.1469 | 1.2972 |
| 2 | Buaran | 2 | 8.0280 | 1.3323 | 3.9745 | 1.3323 |
| 8 | Kedungwuni | 2 | 5.8264 | 1.3357 | 5.1530 | 1.3357 |
| 13 | Siwalan | 2 | 5.5800 | 2.1305 | 6.8609 | 2.1305 |
| 19 | Wonopringgo | 2 | 4.9508 | 2.3477 | 6.7398 | 2.3477 |
| 10 | Lebakbarang | 3 | 9.3090 | 4.7303 | 0.0000 | 0.0000 |
| 5 | Kandangserang | 3 | 8.4748 | 5.7262 | 1.3113 | 1.3113 |
| 12 | Petungkriyono | 3 | 9.9828 | 5.4042 | 1.7944 | 1.7944 |
| 11 | Paninggaran | 3 | 10.4751 | 7.9684 | 3.4926 | 3.4926 |
Jarak anggota terhadap medoid menunjukkan tingkat kedekatan masing-masing kecamatan terhadap pusat cluster-nya. Semakin kecil nilai jarak ke medoid, semakin representatif kecamatan tersebut terhadap cluster yang diikutinya. Pada hasil ini, kecamatan medoid memiliki jarak 0, sedangkan kecamatan dengan jarak yang lebih besar menunjukkan karakteristik yang relatif lebih jauh dari pusat kelompoknya.
profil_mean <- df_best %>%
group_by(cluster) %>%
summarise(across(all_of(num_cols), mean, na.rm = TRUE), .groups = "drop") %>%
mutate(across(where(is.numeric), ~round(.x, 4)))
kable(profil_mean, caption = "Rata-rata Variabel Tiap Cluster dalam Skala Asli")| cluster | padi_prod | sapi | kambing | ayam | itik |
|---|---|---|---|---|---|
| 1 | 19257.00 | 247.1429 | 2706.714 | 1573318.7 | 17060.000 |
| 2 | 5395.75 | 172.2500 | 1553.125 | 217967.6 | 7558.625 |
| 3 | 6998.75 | 2678.7500 | 5170.750 | 147033.2 | 644.000 |
Berdasarkan rata-rata variabel pada tiap cluster, terlihat bahwa masing-masing cluster memiliki karakteristik yang berbeda. Cluster 1 cenderung memiliki nilai produksi padi dan jumlah ayam yang paling tinggi. Cluster 2 menunjukkan nilai yang relatif lebih rendah atau sedang pada sebagian besar variabel. Sementara itu, Cluster 3 memiliki ciri menonjol pada jumlah sapi dan kambing yang paling tinggi, tetapi memiliki nilai itik yang paling rendah. Perbedaan ini menunjukkan bahwa hasil clustering mampu memisahkan kecamatan berdasarkan pola karakteristik pertanian dan peternakan yang berbeda.
pca <- prcomp(X_final, center = TRUE, scale. = FALSE)
var_pct <- (pca$sdev^2) / sum(pca$sdev^2) * 100
plot_df <- data.frame(
pca1 = pca$x[, 1],
pca2 = pca$x[, 2],
cluster = factor(best_model$clustering),
label = df[[id_col]],
is_medoid = FALSE
)
plot_df$is_medoid[best_model$id.med] <- TRUE
hulls <- plot_df %>%
group_by(cluster) %>%
slice(chull(pca1, pca2))ggplot(plot_df, aes(pca1, pca2, color = cluster, shape = cluster)) +
geom_polygon(data = hulls, aes(fill = cluster), alpha = 0.18, color = NA) +
geom_point(size = 3) +
geom_text(aes(label = label), vjust = -0.7, show.legend = FALSE) +
geom_point(
data = subset(plot_df, is_medoid),
shape = 8,
size = 4,
stroke = 1.2,
show.legend = FALSE
) +
theme_minimal() +
labs(
title = paste0("PCA Cluster Plot Metode PAM dengan k = ", best_k),
x = paste0("PC1 (", round(var_pct[1], 1), "%)"),
y = paste0("PC2 (", round(var_pct[2], 1), "%)")
) +
guides(fill = "none")norm_kec <- function(x){
x %>%
str_to_lower() %>%
str_trim() %>%
str_replace_all("\\s+", " ") %>%
str_replace_all("^kec\\.?\\s*", "") %>%
str_replace_all("^kecamatan\\s*", "") %>%
str_replace_all("kab\\.?\\s*", "") %>%
str_replace_all("[^a-z\\s]", "") %>%
str_trim()
}
cluster_df <- df_best %>%
transmute(
kecamatan = .data[[id_col]],
cluster = as.factor(cluster)
) %>%
mutate(kec_key = norm_kec(kecamatan))
gadm_idn_lv3 <- geodata::gadm(country = "IDN", level = 3, path = tempdir()) |>
st_as_sf()
pekalongan_kab <- gadm_idn_lv3 %>%
filter(
str_to_lower(NAME_1) %in% c("jawa tengah", "central java"),
str_to_lower(NAME_2) == "pekalongan",
!str_detect(str_to_lower(NAME_2), "\\bkota\\b")
) %>%
mutate(kec_key = norm_kec(NAME_3))
map_df <- pekalongan_kab %>%
left_join(cluster_df %>% select(kec_key, cluster), by = "kec_key")
label_pts <- map_df %>%
st_point_on_surface() %>%
cbind(st_coordinates(.))
ggplot(map_df) +
geom_sf(aes(fill = cluster), color = "white", linewidth = 0.3) +
geom_text_repel(
data = label_pts,
aes(X, Y, label = NAME_3),
size = 3,
min.segment.length = 0,
box.padding = 0.25,
point.padding = 0.10,
max.overlaps = Inf,
seed = 123
) +
theme_minimal() +
theme(legend.position = "right") +
labs(
title = "Peta Hasil Clustering Kecamatan di Kabupaten Pekalongan",
fill = "Cluster"
)Berdasarkan analisis cluster menggunakan metode K Medoid Partitioning Around Medoids (PAM), diperoleh jumlah cluster terbaik sebanyak 3 cluster dengan nilai average silhouette sebesar 0.4656. Hasil ini menunjukkan bahwa pengelompokan kecamatan di Kabupaten Pekalongan berdasarkan indikator pertanian dan peternakan sudah cukup baik.
Deteksi outlier multivariat menunjukkan adanya 1 outlier, yaitu Paninggaran, namun observasi tersebut tetap dipertahankan karena metode K-Medoids bersifat robust terhadap outlier. Nilai KMO sebesar 0.6494 menunjukkan bahwa data cukup layak untuk analisis multivariat, sedangkan hasil VIF yang seluruhnya di bawah 5 menunjukkan tidak adanya multikolinearitas yang kuat antar variabel.
Secara keseluruhan, hasil clustering berhasil mengelompokkan kecamatan ke dalam beberapa kelompok yang memiliki karakteristik pertanian dan peternakan yang berbeda, sehingga dapat digunakan sebagai dasar untuk memahami pola kemiripan antar kecamatan di Kabupaten Pekalongan.