PENDAHULUAN DAN PENJELASAN DATASET

Indeks Pembangunan Manusia (IPM) atau Human Development Index (HDI) merupakan indikator penting untuk menilai tingkat kesejahteraan suatu negara. Menurut United Nations Development Programme (UNDP), HDI adalah indeks komposit yang mencakup tiga dimensi utama, yaitu kesehatan (angka harapan hidup), pendidikan (rata-rata lama sekolah dan harapan lama sekolah), serta standar hidup layak (penghasilan).

Dataset yang digunakan dalam penelitian ini diperoleh dari platform Kaggle, yaitu Human Development Index Dataset (1990–2022) yang dapat diakses melalui: https://www.kaggle.com/datasets/lucasyukioimafuko/human-development-index-hdr-dataset-1990-2022 . Dataset ini memuat data pembangunan manusia di berbagai negara dalam rentang waktu 1990 hingga 2022.

Analisis dalam penelitian ini difokuskan pada hubungan antara HDI dengan berbagai indikator pendukung, seperti faktor kesehatan, pendidikan, ekonomi, kesetaraan gender, populasi, dan lingkungan, yang dianggap merepresentasikan dimensi penting dalam pembangunan manusia secara komprehensif.

Load Library

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.3
## Warning: package 'lubridate' was built under R version 4.5.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
## Warning: package 'zoo' was built under R version 4.5.3
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(ggplot2)
library(tidyr)
library(dplyr)
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.5.3
library(ggcorrplot)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:DescTools':
## 
##     Recode
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(factoextra)
## Warning: package 'factoextra' was built under R version 4.5.3
## Welcome to factoextra!
## Want to learn more? See two factoextra-related books at https://www.datanovia.com/en/product/practical-guide-to-principal-component-methods-in-r/
library(dbscan)
## Warning: package 'dbscan' was built under R version 4.5.3
## 
## Attaching package: 'dbscan'
## 
## The following object is masked from 'package:stats':
## 
##     as.dendrogram
library(factoextra)
library(ppclust)
## Warning: package 'ppclust' was built under R version 4.5.3
library(cluster)
## Warning: package 'cluster' was built under R version 4.5.3
library(e1071)
## 
## Attaching package: 'e1071'
## 
## The following object is masked from 'package:ggplot2':
## 
##     element
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(flexclust)
## Warning: package 'flexclust' was built under R version 4.5.3
## 
## Attaching package: 'flexclust'
## 
## The following object is masked from 'package:e1071':
## 
##     bclust
library(meanShiftR)
library(fpc)
## Warning: package 'fpc' was built under R version 4.5.3
## 
## Attaching package: 'fpc'
## 
## The following object is masked from 'package:ppclust':
## 
##     plotcluster
## 
## The following object is masked from 'package:dbscan':
## 
##     dbscan

Load Data

data <- read.csv("hdr_general.csv")
target_fitur <- c("life_expectancy", "pop_millions", "expec_yr_school", 
              "mean_yr_school", "gross_inc_percap", "gender_development", 
              "gender_inequality", "co2_emission_tons")
fitur_hdi <- c("hdi", target_fitur)
data <- data %>% select(country, year, all_of(fitur_hdi))

Missing Values: Pengecekan dan Handling

# Cek Missing Values
data_kotor <- data %>%
  summarise(across(where(is.numeric), ~sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "Fitur", values_to = "Jumlah_Data_Kosong") %>%
  filter(Jumlah_Data_Kosong > 0)

print("Total Missing Value Sebelum Imputasi")
## [1] "Total Missing Value Sebelum Imputasi"
print(data_kotor)
## # A tibble: 7 × 2
##   Fitur              Jumlah_Data_Kosong
##   <chr>                           <int>
## 1 hdi                               627
## 2 expec_yr_school                   248
## 3 mean_yr_school                    544
## 4 gross_inc_percap                  139
## 5 gender_development               1784
## 6 gender_inequality                2087
## 7 co2_emission_tons                  85
# Menghapus HDI dengan baris kosong selama 1990-2022
data <- data %>% 
  group_by(country) %>%
  filter(sum(!is.na(hdi))>0) %>%
  ungroup()

data <- data %>%
  group_by(country) %>% # agar tidak tertukar antar negara
  mutate(across(where(is.numeric), ~na.approx(., na.rm = FALSE, rule = 2))) %>%
  ungroup()

sum(is.na(data))
## [1] 1447
jumlah_negara_bersih <- nrow(data)
jumlah_kolom <- ncol(data)
total_missing_setelah <- sum(is.na(data))
cat("Jumlah Negara:", jumlah_negara_bersih, "\n")
## Jumlah Negara: 6732
cat("Jumlah Fitur:", jumlah_kolom, "\n")
## Jumlah Fitur: 11
cat("Total Missing Value Setelah Imputasi:", total_missing_setelah, "\n")
## Total Missing Value Setelah Imputasi: 1447
if(total_missing_setelah ==0){
  cat("Status: Semua Missing Value sudah berhasil ditangani")
} else{
  cat("Status: Missing Value masih ditemukan", total_missing_setelah, "NA tersisa. \n")
}
## Status: Missing Value masih ditemukan 1447 NA tersisa.
data <- data %>%
  group_by(country) %>%
  fill(everything(), .direction = "updown") %>% # mengisi dengan data terdekat
  ungroup() %>%
  na.omit() # Menghapus negara yang benar-benar tidak ada datanya

cat("Total Missing Value Akhir:", sum(is.na(data)))
## Total Missing Value Akhir: 0

Feature Engineering: Menghitung Rata-rata Per Negara

data_bersih <- data %>%
  group_by(country) %>%
  summarise(across(all_of(fitur_hdi), \(x) mean(x, na.rm=TRUE))) %>%
  ungroup()

print(paste("Jumlah baris sekarang:", nrow(data_bersih)))
## [1] "Jumlah baris sekarang: 177"

Pada Perhitung rata-rata ini sama saja dilakukan kompresi ukuran data dengan mengurangi baris dari dataset sehingga dataset hanya berisi *177 baris**. Namun, pada baris-baris tersebut ditemukan beberapa baris yang kurang relevan dengan tujuan analisis yakni “Very high human development”, “High human development”, “Medium human development”, “Low human development”, “World”, yang akan dihapus.

list_hapus <- c("Very high human development", 
                "High human development", 
                "Medium human development", 
                "Low human development",
                "World")
data_bersih <- data_bersih %>%
  filter(!(country %in% list_hapus))

Outlier: Pengecekan & Handling

# Mengecek Outlier: Visualisasi
boxplot(data_bersih %>%
          select(all_of(target_fitur)),
        col="skyblue", las=2)

# Mengecek Outlier: List Kolom
count_outliers <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR_val <- Q3 - Q1
  lower_bound <- Q1 - (1.5 * IQR_val)
  upper_bound <- Q3 + (1.5 * IQR_val)
  return(sum(x < lower_bound | x > upper_bound, na.rm = TRUE))
  
  x <- ifelse(x < lower_bound, lower_bound, x)
  x <- ifelse(x > upper_bound, upper_bound, x)
  return(x)
}

data_outlier <- data_bersih %>%
  summarise(across(all_of(target_fitur), count_outliers)) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Jumlah_Outlier") %>%
  filter(Jumlah_Outlier > 0)
print("Kolom dengan outlier:")
## [1] "Kolom dengan outlier:"
print(data_outlier)
## # A tibble: 5 × 2
##   Variabel           Jumlah_Outlier
##   <chr>                       <int>
## 1 pop_millions                   23
## 2 expec_yr_school                 1
## 3 gross_inc_percap               10
## 4 gender_development              5
## 5 co2_emission_tons              11
# Penanganan Outlier: Winsorize Method
winsorize_manual <- function(x, probs = c(0.20, 0.80)) {
  q_bounds <- quantile(x, probs = probs, na.rm = TRUE)
  x[x < q_bounds[1]] <- q_bounds[1]
  x[x > q_bounds[2]] <- q_bounds[2]
  return(x)
}
data_bersih <- data_bersih %>%
  mutate(across(all_of(target_fitur), ~ winsorize_manual(.)))

boxplot(data_bersih %>%
          select(all_of(target_fitur)),
        col="tomato", las=2)

Scaling: Perubahan Rentang Data

list_scale <- data_bersih %>%
  select(all_of(target_fitur))
data_scaled <- scale(list_scale)
data_scaled <- as.data.frame(data_scaled)

Pengujian Data

# Uji Multikolinearitas
corr_matrix <- cor(data_scaled)

corrplot(corr_matrix, method = "color",
         type = "lower", addCoef.col = "black",
         number.cex = 0.8, tl.col = "black",
         mar = c(0,0,1,0))

Berdasarkan hasil pengujian multikolinearitas tersebut, ditemukan 2 variabel yang memiliki hubungan tinggi yakni variabel co2_emission_tons dan gross_income_percapita, yang menyebabkan salah satu variabel sebaiknya dihapus untuk mencegah redundasi data. Karena analisis ini berfokus pada sosial-ekonomi, maka variabel yang dihapus adalah variabel co2_emission_tons.

# Menghapus Varibel co2_emission_tons
data_scaled <- data_scaled %>%
  select(-co2_emission_tons)
target_fitur <- setdiff(target_fitur, "co2_emission_tons")
# Uji Hopkins
hopkins_test <- get_clust_tendency(data_scaled,
                                   n = nrow(data_scaled)-1,
                                   graph = FALSE)
## Warning: Hopkins statistic uses the corrected formula (Wright 2022); results
## differ from legacy factoextra. Set options(factoextra.warn_hopkins = FALSE) to
## silence this warning.
print(paste("Hopkins Statistic:", round(hopkins_test$hopkins_stat, 4)))
## [1] "Hopkins Statistic: 0.9775"

Berdasarkan hasil pengujian, diperoleh nilai Hopkins sebesar 0.9863 yang menunjukkan bahwa dataset memiliki kecenderungan yang kuat untuk membentuk Cluster. Oleh karena itu, data dinyatakan layak untuk dianalisis lebih lanjut.

# Uji Densitas
data_ujidensi <- data_scaled %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")
ggplot(data_ujidensi, aes(x = Nilai, fill = Variabel)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~Variabel, scales = "free") +
  theme_minimal() +
  labs(title = "Distribusi Densitas Fitur", x = "Z-Score", y = "Density")

Berdasarkan hasil visualisasi distribusi densitas, dapat diketahui bahwa sebagian besar variabel menunjukkan pola distribusi yang cenderung skewed (miring). Beberapa variabel menunjukkan kemiringan ke kiri dan sebagian data menunjukkan kemiringan ke kanan, yang artinya data tidak berdistribusi normal. Kondisi ini menunjukkan bahwa data memiliki kompleksitas yang cukup tinggi sehingga penggunaan metode clustering masih tetap relevan, khususnya metode clustering yang tidak bergantung pada asumsi distribusi normal.

# Uji Bandwidth
kNNdistplot(data_scaled, k = 16)
abline(h = 1.2, col = "red", lty = 2)

Berdasarkan hasil visualisasi, terlihat bahwa kurva mulai mengalami perubahan lengkungan yang signifikan pada nilai jarak sekitar 1.2. Hal ini menunjukkan bahwa, titik-titik yang memiliki jarak ke tetangga ke-16 dibawah 1.2 berada dalam area yang cukup padat untuk membentuk sebuah cluster.

# Uji Elbow & Uji Silhouette
elbow_test <- fviz_nbclust(data_scaled,
                           kmeans,
                           method = "wss") +
  labs(title = "Elbow Method")
silhouette_test <- fviz_nbclust(data_scaled,
                                kmeans,
                                method = "silhouette") +
  labs(title = "Silhouette Method")

grid.arrange(elbow_test, 
             silhouette_test,
             ncol=2)

Berdasarkan hasil visualisasi, kedua Uji menunjukkan hasil yang konsisten yaitu jumlah cluster paling optimal berada pada k = 2. Pada Grafik Silhouette menunjukkan puncak nilai rata-rata tertinggi berada pada titik K = 2, hal ini juga didukung oleh grafik Elbow yang menunjukkan penurunan variansi paling drastis terjadi di k = 2. Meskipun k = 2 merupakan jumlah cluster yang paling optimal, pada penelitian ini jumlah cluster k = 3 dipilih sebagai jumlah cluster final untuk mendapatkan segmentasi yang lebih mendalam serta spesifik pada penelitian.

CLUSTERING: 5 Metode Utama

K-Means

set.seed(123)
jumlah_k <- 3

# Pembuatan Cluster
kmeans_clus <- kmeans(data_scaled,
                      centers = jumlah_k,
                      nstart = 25)
# Visualisasi
fviz_cluster(kmeans_clus, data = data_scaled, geom = "point", 
             ellipse.type = "convex", main = "K-Means Clustering")

# Dimasukan dalam Data Bersih
data_bersih$cluster_kmeans <- as.factor(kmeans_clus$cluster)

Berdasarkan hasil pengujian elbow dan silhouette yang telah dilakukan maka untuk clusteting K-Means dan K-Medians nantinya akan menggunakan 3 cluster untuk memperluas pandangan analisis. Dalam sebelum pembuatan Cluster menggunakan metode K-means ini, terdapat kode set.seed(123). Kode tersebut digunakan untuk menjaga konsistensi angka pusat (centroid) agar setiap kode dijalankan dan cluster k-means dibuat hasilnya tidak akan berbeda-beda.

Lalu, pada pembuatan cluster dengan fungsi kmeans() digunakan parameter nstart=25 yang digunakan untuk mencoba 25 konfigurasi awal yang berbeda dan memilih yang terbaik. Setelah 3 cluster berhasil dibentuk, hasil cluster akan dimasukan pada data_bersih (data yang belum discale) dengan mengubah tipe data dari hasil clustering menjadi bentuk factor menggunakan fungsi as.factor()

Untuk pembuatan visualisasi digunakan metode PCA untuk mereduksi dimensi data sehingga visualisasi dapat dilakukan menggunakan fungsi fviz_cluster dari library factorextra. Parameter yang digunakan untuk menyusun visualisai PCA K-Means ini adalah geom = "points" untuk menampilkan data sebagai titik dan ellipse.type = "convex" yang membuat garis pembatas antar cluster. Hasil visualisasi terlihat bahwa cluster terbagi menjadi 3 kelompok sesuai dengan jumlah_k yang di-set. Cluster 3 merupakan kelompok yang condong ke kiri, cluster 1 merupakan kelompok menengah, dan cluster 2 merupakan kelompok yang condong ke kanan.

K- Median

# Pembentuk Cluster
kmedian_clus <- kcca(data_scaled,
                     k = jumlah_k,
                     family = kccaFamily("kmedians"))
## Found more than one class "kcca" in cache; using the first, from namespace 'flexclust'
## Also defined by 'kernlab'
## Found more than one class "kcca" in cache; using the first, from namespace 'flexclust'
## Also defined by 'kernlab'
# Memasukan ke Dat Bersih
data_bersih$cluster_kmedians <- as.factor(predict(kmedian_clus))

# Visualisasi
cluster_label <- predict(kmedian_clus)

fviz_cluster(list(data = data_scaled, cluster = cluster_label), 
             geom = "point", 
             ellipse.type = "convex",
             palette = "jco",
             main = "Visualisasi K-Medians (PCA Plot)",
             ggtheme = theme_minimal())

Dalam pembentukan cluster menggunakan metode K-Median menggunakan fungsi kcca dari library flexclust terdapat penggunaan parameter family = kccaFamily("kmedians"). Penggunaan parameter family tersebut dikarenakan metode perhitungan jarak yang digunakan untuk pembuatan cluster k-medians ini menggunakan metode Manhattan.

Untuk menggabungkan hasil clustering pada tabel data_bersih dilakukan perubahan tipe data menjasi faktor menggunakan fungsi as.factor namun ditambah dengan fungsi predict() . Tujuan ditambahkan fungsi tersebut adalah struktur hasil cluster yang dihasilkan oleh metode k-medians berbeda dengan k-means, fungsi predict() akan mendapatkan vektor angka cluster yang akan dimasukan pada data bersih.

Untuk kode pembuatan visualisasi tidak jauh berbeda dengan kode pada k-means, namun terdapat parameter tambahan palette = "jco" yang bertugas untuk memilih warna cluster dengan warna yang bersih, lalu parameter ggtheme = theme_minimal() digunakan untuk menfokuskan tampilan pada sebaran data. Hasil dari clustering yang dilakukan menggunakan metode k-median ini terbagi menjadi tiga cluster utama dengan masing-masing karakteristik yang berbeda dengan cluster metode k-means. Cluster 1 merupakan kelompok yang condong ke kanan, cluster 2 kelompok menengah, dan cluster 3 kelompok yang condong ke kiri.

Meanshift

# Pembentukan Cluster
meanshift_test <- meanShift(as.matrix(data_scaled),
                            bandwidth = rep(1.2, ncol(data_scaled)))

# Memasukan ke Data Bersih
data_bersih$cluster_meanshift <- as.factor(meanshift_test$assignment)

# Visualisasi
means_clus <- as.factor(meanshift_test$assignment)
fviz_cluster(list(data = data_scaled, cluster = means_clus),
             geom = "point",
             ellipse.type = "norm",
             palette = "Set2",
             main = "Visualisasi Mean Shift Clustering",
             ggtheme = theme_minimal())
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning: Removed 46 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'
## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

Clustering dengan metode meanshift bekerja dengan menggeser titik ke arah daerah yang padat data paling tinggi secara iterasi. Dalam metode meanshift ini digunakan fungsi meanShift dengan parameter as.matrix() untuk data yang digunakan, tujuannya adalah untuk mengubah data_scaled ke bentuk matriks sehingga dapat diproses oleh fungsi meanShift. Parameter lain yang digunakan adalah bandwidth= rep(1.2, ncol(data_scaled)), parameter tersebut bertujuan untuk menentukan jumlah cluster yang dapat dibuat oleh fungsi meanShift, disini digunakan radius 1.2 untuk setiap variabel.

Setelah cluster berhasil dibentuk, terlihat divisualisasi PCA hasil cluster yang dibentuk oelh funsgi meanShift berjumlah 25 cluster. Hasil tersebut akan dimasukan pada data data_bersih dengan bantuan fungsi as.factor dan $assigment yang mengambil hasil clustering dari meanshift_test.

DBSCAN (Density-Based Spatial Clustering of Applications with Noise)

# Pembentukan Cluster
dbscan_test <- dbscan(data_scaled,
                      eps=2.2, MinPts = 12)

# Memasukan ke Data Bersih
data_bersih$cluster_dbscan <- as.factor(dbscan_test$cluster)

# Visualisasi
fviz_cluster(dbscan_test, 
             data = data_scaled, 
             geom = "point", 
             ellipse.type = "convex", 
             main = "Visualisasi DBSCAN Clustering",
             ggtheme = theme_minimal())

DBSCAN adalah metode yang mencari area dengan kepadatan data dan mencari data dengan jumlah tertentu pada jarak tertentu, jika ditemukan akan menjadi sebuah cluster. Dalam metode DBSCAN ini digunakan fungsi dbscan() dengan parameter custom eps = 2.2 yang digunakan sebagai jarak minimum antar dua titik agar dapat dianggap sebagai tetangga dan MinPts = 12 yang digunakan sebagai jumlah minimum titik dalam radius eps untuk menjadi sebuah cluster. DBSCAN memiliki keunikan dimana dirinya termasuk robust pada data noise dan dapat mendeteksi data noise.

Hasil dari visualisasi PCA yang dilakukan pada hasil clustering metode DBSCAN ini terlihat bahwa DBSCAN dengan parameter eps=2.2 dan MinPts = 12 gagal memisahkan dataset yang digunakan, sehingga cluster yang terbentuk hanyalah 1 cluster dengan semua data di dalamnya. Selain itu, DBSCAN tidak mendeteksi adanya noise dalam data yang digunakan untuk analisis ini.

Dikarenakan hasil clustering yang dilakukan dengan metode DBSCAN gagal membentuk kelompok maka untuk perbandingan metrik performa di akhir analisis, DBSCAN *tidak diuji kembali** karena performanya yang kurang pada parameter yang digunakan.

Fuzzy C-Means

# Uji Fuzzy Partition Coefficient
matrix_data <- as.matrix(data_scaled[,
                                     sapply(data_scaled,
                                            is.numeric)])

hitung_fpc <- function(fcm_obj) {
  u <- fcm_obj$membership
  n <- nrow(u)
  fpc <- sum(u^2) / n
  return(fpc)
}
k_range <- 2:5
hasil_fpc <- numeric(length(k_range))

for (i in seq_along(k_range)) {
  res <- cmeans(matrix_data, centers = k_range[i], m = 2)
  hasil_fpc[i] <- hitung_fpc(res)
  cat("K =", k_range[i], "| FPC manual:", hasil_fpc[i], "\n")
}
## K = 2 | FPC manual: 0.7579984 
## K = 3 | FPC manual: 0.614591 
## K = 4 | FPC manual: 0.5409891 
## K = 5 | FPC manual: 0.4734282
evaluasi_e1071 <- data.frame(k = k_range, Skor_FPC = hasil_fpc)
print(evaluasi_e1071)
##   k  Skor_FPC
## 1 2 0.7579984
## 2 3 0.6145910
## 3 4 0.5409891
## 4 5 0.4734282

Pada metode Fuzzy Cmeans terdapat uji yang harus dilakukan untuk menentukan jumlah cluster terbaik yang sebaiknya dibuat untuk analisis, uji tersebut adalah uji Fuzzy Partition Coefficient (FPC).

Metode ini dihitung dengan data yang berbentuk matriks, oleh karena itu sebelum dilakukan perhitungan data perlu diubah menjadi matriks menggunakan fungsi as.matrix(). Dalam menghitung FPC akan dilakukan iterasi untuk menentukan jumlah k terbaik antara 2-5 cluster. Berdasarkan hasil perhitungan FPC, ditemukan hasil k=2 memilik skor FPC paling tinggi yakni 0.7579984, disusul dengan k=3 dengan skor 0.6145910. Seperti halnya hasil uji elbow dan silhouette yang menunjukan hasil sama, pada analisis Fuzzy Cmeans ini juga digunakan k=3 untuk pandangan analisis yang lebih luas.

# Pembentukan Cluster
fcm_test <- cmeans(matrix_data, centers = 3, m=2)

# Memasukan ke Data Bersih
data_bersih$cluster_fcm <- as.factor(fcm_test$cluster)

# Visualisasi
fviz_cluster(list(data = matrix_data, cluster = fcm_test$cluster),
             geom = "point", 
             ellipse.type = "norm",
             palette = "jco",
             main = "Visualisasi Fuzzy C-Means Clustering",
             ggtheme = theme_minimal())

Fuzzy Cmeans memiliki sifat Soft Clustering dimana program tidak memaksa suatu data untuk masuk ke satu cluster saja, hal ini dikarenakan Fuzzy Cmeans menganggap setiap data memiliki bobot untuk setiap cluster. Namun untuk memperjelas clusternya, digunakan cluster dengan bobot tertinggi untuk pelabelan.

Dalam pembentukan cluster dengan metode Fuzzy ini digunakan fungsi cmeans dengan parameter center = 3 yang berarti pusat cluster terdapat 3 sesuai dengan hasil uji FPC yang telah dilakukan dan parameter m = 2 yang mengatur fleksibilitas data dalam pengelompokan, sehingga data dapat berada di beberappa cluster.

Hasil dari clustering yang dilakukan dengan metode Fuzzy Cmeans ini menghasilkan 3 cluster dengan karakteristik yang berbeda-beda. Cluster 1 kelompok yang condong ke kiri, cluster 2 kelompok yang condong ke kanan, dan cluster 3 kelompok menengah. Diantara tiga cluster tersebut terdapat tumpang tindih yang jelas pada cluster 2 dan 3, hal ini mengartikan bahwa bobot data tersebut pada cluster 2 dan 3 cukup mirip. Meski begitu untuk pelabelan akhir nanti, cluster yang digunakan adalah cluster dengan bobot paling tinggi pada data tersebut.

Uji Metode Terbaik

Pengujian hanya dilakukan untuk metode yang dapat membuat pengelompokan dengan baik, yakni K-Means, K-Medians, Meanshift, dan Fuzzy C-Means.

# Pembuatan Fungsi Uji Metrik
hitung_metrik <- function(data, labels){
  dist_matrik <- dist(data)
  stats <- cluster.stats(dist_matrik, labels)
  
  return(c(
    Silhouette = stats$avg.silwidth,
    Davies_Bouldin = stats$db,
    Dunn_index = stats$dunn
  ))
}

# Perbandingan Metode
perbandingan_metode <- rbind(
  K_Means    = hitung_metrik(data_scaled, kmeans_clus$cluster),
  K_Medians  = hitung_metrik(data_scaled, cluster_label),
  Mean_Shift = hitung_metrik(data_scaled, meanshift_test$assignment),
  Fuzzy_CM   = hitung_metrik(data_scaled, fcm_test$cluster)
)

print(round(perbandingan_metode, 4))
##            Silhouette Dunn_index
## K_Means        0.3635     0.1849
## K_Medians      0.3607     0.1824
## Mean_Shift     0.0304     0.0842
## Fuzzy_CM       0.3629     0.1841

Untuk membandingkan dan menentukan metode yang dapat melakukan clustering terbaik dilakukan perhitungan metrik untuk silhouette, davies-bouldin, dan dunn index. Silhouette digunakan untuk mengukur kedekatan titik dengan titik lain di satu cluster yang sama, lalu dibandingkan dengan cluster lain pada rentang -1 ke 1. Davies-Bouldin digunakan untuk mengukur rasio jarak dalam cluster antar cluster. Dunn Index digunakan untuk mengukur rasio minimum antar cluster terhadap diameter cluster maksimum.

Dalam kode perbandingan metode tersebut dibuat fungsi hitung_metrik yang menghitung semua metrik yang akan diujikan dalam satu waktu dengan data yang sudah terstandarisasi dan labels yang berupa vektor identitas cluster. Untuk menjalakan perhitungan dibuat fungsi perbandingan_metode yang memanggil fungsi hitung_metrik dan diimplementasikan pada setiap metode clustering yang dilakukan.

Hasil dari perhitungan matriks yang dilakukan pada setiap metode clustering yang digunakan menunjukan bahwa metode terbaik adalah K-Means dengan skor Silhouette 0.3635 yang membuatnya dinilai mampu membuat struktur cluster dengan cukup kuat, selain itu skor Dunn Index untuk metode ini juga paling tinggi di antara lainnya yang berarti metode k-means berhasil memisahkan tiap cluster dengan baik dan padat untuk tiap clusternya. Meski begitu, terdapat dua metode lain yang memiliki skor Silhouette dan Dunn Index dengan selisih tipis dengan metode K-Means yakni Fuzzy Cmeans dan K-Medians, hal ini menjelaskan bahwa kedua metode tersebut juga mampu melakukan clustering dengan baik dan hampir mirip dengan K-Means.

Berdasarkan hasil perbandingan metode yang telah dilakukan, maka analisis data eksploratif yang dilakukan pada tahapan selanjutnya akan menggunakan hasil clustering dengan metode K-Means. # Analisis Data Eksploratif ## Profiling Cluster

profil_cluster <- data_bersih %>%
  group_by(cluster_kmeans) %>%
  summarise(
    n = n(),
    Mean_HDI = mean(hdi),
    Mean_LifeExp = mean(life_expectancy),
    Mean_Popmillions = mean(pop_millions),
    Mean_EYS = mean(expec_yr_school),
    Mean_MYS = mean(mean_yr_school),
    Mean_GD = mean(gender_development),
    Mean_GI = mean(gender_inequality),
    Mean_GNI = mean(gross_inc_percap),
    Mean_CO2 = mean(co2_emission_tons)
  ) %>%
  arrange(desc(Mean_HDI))

print(profil_cluster)
## # A tibble: 3 × 11
##   cluster_kmeans     n Mean_HDI Mean_LifeExp Mean_Popmillions Mean_EYS Mean_MYS
##   <fct>          <int>    <dbl>        <dbl>            <dbl>    <dbl>    <dbl>
## 1 2                 57    0.834         75.1             13.9    14.0     10.1 
## 2 1                 59    0.685         70.6             17.0    12.3      8.11
## 3 3                 56    0.472         61.2             20.1     9.72     4.88
## # ℹ 4 more variables: Mean_GD <dbl>, Mean_GI <dbl>, Mean_GNI <dbl>,
## #   Mean_CO2 <dbl>

Tahapan pertama dalam analisis data eksploratif ini adalah melakukan profiling tiap cluster atau mengenal baik karakteristik tiap cluster. Untuk melakukan profiling tersebut dilakukan agregasi berdasarkan hasil cluster k-means yang nantinya tiap cluster akan dihitung rata-rata per variabelnya.

Untuk membuat perhitungan rata-rata variabel untuk tiap cluster digunakan fungsi group_by(cluster_kmeans) yang mengelompokan data berdasarkan cluster mereka. Lalu, untuk mempermudah menginterpretasikan hasil perhitungan digunakan fungsi summarise() untuk meringkas. Selain menghitung rata-rata setiap variabel per cluster menggunakan fungsi mean(), dilakuakn perhitungan jumlah data tiap cluster dengan fungsi n = n(). Setelah perhitungan dilakukan, data akan diurutkan dengan rata-rata HDI (Human Development Index) paling tinggi menggunakan fungsi arrange(desc(Mean_HDI)) sehingga dapat dimudahkan untuk interpretasi nantinya.

Berdasarkan hasil perhitungan rata-rata variabel untuk tiap cluster, berikut merupakan penjelasan detail untuk masing-masing cluster yang dibangun menggunakan metode K-Means:

Cluster 2: HDI Tinggi Cluster kedua ini memiliki rata-rata peningkatan pembangunan manusia paling tinggi dengan nilai rata-rata HDI sebesar 0.833 diikuti dengan nilai pendidikan, ekspektasi kehidupan dan ekonomi yang tinggi. Hal ini cukup menjelaskan bahwa cluster 2 merupakan kelompok negara-negara maju. Negara maju berdasarkan hasil clustering ini memiliki karakteristik rata-rata HDI yang tinggi dengan aspek ekonomi, pendidikan, usia kehidupan penduduknya yang tinggi. Namun, negara-negara pada cluster tersebut memiliki emisi CO2 paling tinggi diantara cluster lainnya, yang berarti negara-negara maju ini menjadi penyumbang gas emisi paling tinggi.

Cluster 1: HDI Menengah Cluster pertama memiliki rata-rata HDI tertinggi kedua dengan rata-rata HDI sebesar 0.684 diikuti dengan nilai rata-rata pendidikan, ekspektasi kehidupan dan ekonomi yang cukup tinggi, namun tetap di bawah cluster 2 ini. Hal ini cukup untuk menyimpulkan bahwa cluster 1 merupakan kelompok negara-negara berkembang. Berdasarkan hasil rata-rata yang telah dilakukan cluster 1 ini memiliki karakteristik perkembangan manusia yang cukup baik meski kalah dengan negara maju, berbeda dengan negara maju yang minim isu gender, pada negara berkembang isu gender cukup banyak didapati bahkan naik sekiranya dua kali lipat dibandingkan dengan isu gender pada negara maju. Namun, penyumbangan gas emisi dari negara berkembang ini turun drastis dibandingkan dengan negara maju.

Cluster 3: HDI Rendah Cluster ketiga ini memiliki nilai rata-rata HDI paling rendah dibandingkan dengan dua cluster lainnya yakni 0.472 diikuti dengan nilai ekonomi, pendidikan, dan ekspektasi hidup yang sama rendahnya. Berdasarkan rata-rata tersebut dapat disimpulkan bahwa cluster 3 merupakan kelompok negara tertinggal. Sebagai negara yang tertinggal ekspektasi hidup, nilai ekonomi, dan jalan pendidikan penduduk mereka tergolong rendah bahkan sangat rendah dikarenakan kurangnya fasilitas yang memadai. Lalu, untuk isu gender pada negara tertinggal ini sangat tinggi yang menandakan bahwa penduduk disana sering kali melanggar isu antar gender. Meski begitu, sebagai negara tertinggal yang minim sektor industrial dan kendaran berbasis gas, kelompok ini menjadi kelompok dengan penyumbang gas emisi paling rendah di antara tiga cluster lainnya.

Visualisasi PCA: Struktur Parsial

fviz_cluster(list(data = data_scaled, cluster = data_bersih$cluster_kmeans),
             geom = "point",
             ellipse.type = "convex",
             palette = "jco",
             main = "Peta Klaster K-Means (PCA Plot)")