Tugas Akhir Teknik Pembelajaran Mesin

Email: zenr.prog@gmail.com

Github: https://github.com/Zen-Rofiqy

️Rpubs: https://rpubs.com/ZenR_Prog/

Kaggle: https://www.kaggle.com/zenrofiqy

Linkedin: https://www.linkedin.com/in/angga-fathan-rofiqy/

Perkenalan

Angga Fathan Rofiqy adalah mahasiswa semester 6 di IPB University yang tengah menempuh studi di bidang Statistika. Untuk memenuhi tugas akhir mata kuliah Teknik Pembelajaran Mesin, Angga dan timnya melakukan penelitian berjudul “Perbandingan Metode Clustering dalam Analisis Kesejahteraan di Indonesia tahun 2021: K-Means, Fuzzy C-Means, dan Gaussian Mixture Model (GMM).” Penelitian ini dipublikasikan di RPubs dengan tujuan memberikan bahan pembelajaran yang bermanfaat bagi semua orang.

Anggota Kelompok

Kelompok 8 Pararel 2

Nama	NIM
Angga Fathan Rofiqy	G1401211006
Natasha Muti Hafiza	G1401211019
Arfiah Kania Sektiaruni	G1401211023
Naswa Nabila Zahrani	G1401211056
Kamilah Nurul Azizah	G1401211073
Angel Martha Pradina .P.	G1401211105

Data

Input Data

Data untuk analais bisa di dapatkan di data.github dan data kode provinsi ada di kode.github

raw.data <- read.csv("https://raw.githubusercontent.com/Zen-Rofiqy/STA1382-TPM/main/Project/Tugas%20Kelompok/Data.csv")
data <- raw.data # Untuk menjaga data awal
colnames(data) <- c("Prov","X1","X2",'X3','X4','X5','X6','X7')

kode <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1382-TPM/main/Project/Tugas%20Kelompok/Kode%20Provinsi.csv")

datatable(raw.data,  options = list(pageLength = 5))

Check Kualitas Data

Cek tipe data

str(raw.data)

## 'data.frame':    34 obs. of  8 variables:
##  $ Provinsi   : chr  "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
##  $ uhh        : num  70 69.2 69.6 71.7 71.2 ...
##  $ puskesmas  : int  363 618 282 238 212 352 185 310 65 92 ...
##  $ imun       : num  21.6 42.8 42.8 44.6 48.3 ...
##  $ tenagamedis: int  2683 5611 2626 2634 1533 2464 787 2392 751 1188 ...
##  $ bpjs       : num  83.4 32 37.9 27.1 25 ...
##  $ sanitasi   : num  77.5 82 68.7 83.6 80.4 ...
##  $ sumberair  : num  88.8 90.9 83.4 89.8 79.7 ...

Cek data hilang

colSums(is.na(raw.data))

##    Provinsi         uhh   puskesmas        imun tenagamedis        bpjs 
##           0           0           0           0           0           0 
##    sanitasi   sumberair 
##           0           0

Tidak ada data hilang pada setiap kolom.

Eksplorasi data

Sebaran data

colors <- c("#0f2c3a","#357a87","#6ec1c1","#a9ded7","#9a3b1a","#ce572b","#e9865b")
titles <- c("\nX1 - UHH","\nX2 - Puskesmas","\nX3 - Imun","\nX4 - Tenaga Medis",
            "\nX5 - BPJS", "\nX6 - Sanitasi","\nX7 - Sumber Air")
bin <- c(.8, 85, 10, 700, 7, 5, 6)
all_plots <- list()

for (i in 1:(ncol(data)-1) ) {
  plots <- list()
  
  # Histogram with density plot
  p1 <- ggplot(data, aes(x = data[,1+i] )) +
    geom_histogram(aes(y = ..density.., weight=data[,1+i] ), 
                   binwidth = bin[i], 
                   fill = colors[i], color = "black", lwd=1) +
    geom_density(color = "black", lwd=2, fill = colors[i], alpha = .5) +
    ylab("Density") + theme_void() +
    xlab(colnames(data[,1+i])[i]) 
  
  # Boxplot
  p2 <- ggplot(data, aes(y = data[,1+i] )) +
    geom_boxplot(fill = colors[i], color = "black", 
                 outlier.size = 3, lwd=1, alpha=.5) +
    xlab("Nilai") + theme_set(theme_void() + theme(axis.text.x = element_text())) + 
    coord_flip()
  
  # Gabung plot
  combined_plot <- plot_grid(p1, p2, ncol = 1, align = 'v', 
                             rel_heights = c(1, 0.15))
  
  # Menambahkan judul utama
  title <- ggdraw() + 
    draw_label(titles[i], fontface = 'bold', x = 0.5, 
               hjust = 0.5, size = 30)
  
  # Menampilkan plot dengan judul utama
  plots <- 
  plot_grid(title, combined_plot, ncol = 1, rel_heights = c(0.1, 1))
  
  all_plots[[i]] <- plots
}
chart <-
plot_grid(plotlist = all_plots, ncol=4)
chart

#Export Chart
ggsave("Sebaran Data.png", chart, path = export.chart,
        dpi = 300, height = 15, width = 30)

Peubah X2, X4, dan X5 menyebar dengan menjulur ke kanan. Ini bisa dilihat dari density plot dan histogram nya yang menjulur ke kanan, serta dengan adanya outlier di sebelah kanan pada boxplotnya. Artinya mayoritas provinsi memiliki nilai yang jauh lebih sedikit daripada provinsi lainnya. Dalam kasus ini nilai yang dimaksud adalah banyaknya puskesmas, tenaga medis, dan persentase penerima BPJS.

Peubah X3, X6, dan X7 menyebar dengan menjulur ke kiri. Ini bisa dilihat dari density plot dan histogram nya yang menjulur ke kiri, serta dengan adanya outlier di sebelah kiri pada boxplotnya. Artinya masih ada beberapa provinsi yang nilai nya tertinggal jauh dari pada kebanyakan provinsi. Dalam kasus ini nilai yang dimaksud adalah Persentase balita sudah imunisasi dasar lengkap, Persentase rumah tangga dengan sanitasi layak, dan Persentase rumah tangga dengan pelayanan sumber air minum layak.

Sedangkan untuk umur harapan hidup nampak hampir simetris jika dilihat dari histogram, density plot dan boxplot nya. Artinya sebaran UHH di setiap provinsi terdistribusi secara merata di sekitar rata-rata, tanpa ada kemenjuluran yang signifikan ke arah tertentu.

Sebaran per provinsi

# List of variable names (excluding Provinsi and Cluster)
variables <- names(raw.data)[!names(raw.data) %in% c("Provinsi")]

plots <- list()

for (var in variables) {
  # Sort data based on the current variable
  sorted_data <- raw.data %>%
    arrange(!!sym(var)) %>%
    mutate(Provinsi = factor(Provinsi, levels = Provinsi))
  
  p <- ggplot(sorted_data, aes_string(x = var, y = "Provinsi")) +
    geom_bar(stat = "identity", fill = "#357a87") +
    labs(x = var, y = "Provinsi") +
    theme_minimal()
  
  plots[[var]] <- p
}

combined_plot <- plot_grid(plotlist = plots, ncol = 4)
combined_plot

#Export Chart
ggsave("Bar Chart.png", combined_plot, path = export.chart,
        dpi = 300, height = 15, width = 30)

Untuk sebaran lebih detail nya seperti ini. Bisa dilihat bahwa selaras dengan yang sudah diinterpretasikan di plot sebelumnya, bahwa UHH hampir menyebar simetris, lalu X2, X4, dan X5 menyebar dengan menjulur ke kanan yang ditandai banyaknya bar yang kecil, lalu terakhir X3, X6, dan X7 menyebar dengan menjulur ke kiri yang ditandai dengan ada beberapa bar yang ukurannya lebih kecil dari mayoritas.

Hubungan antar peubah

# Standarisasi Data
dt <- scale(data[,-1])
# Korelasi
pairs.panels(dt,
             method = "pearson", # correlation method
             hist.col = "#357a87", #Coloring histogram
             density = TRUE,  # show density plots
             ellipses = TRUE, # show correlation ellipses
             smooth = TRUE, #show loess smooths
             pch = 20, #Scatter = cirlce / dot
             rug = TRUE, #Rug under histogram
             stars = TRUE #Significance of corr
             )

Jika kita melihat dari scatterplotnya, maka kita akan melihat bahwa pola hubungan antar peubah itu beragam, ada yang memilki hubungan linear ada juga yang nampaknya tidak meiliki hubungan linear. Jika dihitung hubungan linearnya dengan korelasi pearson, maka kita bisa melihat bahwa mayoritas hubungan antarr peubahnya itu lemah. Hubungan lineaer yang paling kuat adalah hubungan antara peubah X2 dan X4 dengan nilai 0.78. Bisa dilihat juga dari scatterplitnya yang memilki pola lurus positif. Ini didukung dnegan garis LOESS Smooth nya yang paling lurus dan paling 45 derajat diantara yang lain. Juga dari bentuk elips korelasinya yang nampak lancip nan bersudut 45 derajat. Yang mana ini tentu make sense karena Jumlah Puskesmas selaras dengan juumlah tenaga medis yang ada.

Multikolinearitas

CekVIF <- function(data) {
  corr = as.matrix(cor(data))
  VIF = diag(solve(corr))
  
  return(VIF)
}

CekVIF(data[,-1])

##       X1       X2       X3       X4       X5       X6       X7 
## 2.027411 3.895705 1.356273 4.427064 1.131544 2.935106 1.892231

Pada setiap variabel tidak ada yang memiliki nilai VIF > 10 sehingga antar variabel tidak terjadi multikolinearitas.

Clustering

GMM

n Cluster Optimum

Pemilihan banyak culster

Pada dasarnya nilai k optimum dapat ditentukan sendiri oleh peneliti, selain itu nilai k dapat ditentukan juga menggunakan bantuan nilai dari setiap sebaran gaussian yang dihasilkan dan dilihat mana yang memiliki nilai BIC terbesar. Tak harus melihat nilai dari sebarannya, ini juga bisa dibuat plot agar mudah untuk menentukan mana jumlah cluster yang optimal.

# GMM clustering dengan mclust
gmm_model <- Mclust(dt)

gmm_model$BIC

## Bayesian Information Criterion (BIC): 
##         EII       VII       EEI       VEI       EVI       VVI       EEE
## 1 -696.5206 -696.5206 -717.6788 -717.6788 -717.6788 -717.6788 -688.7975
## 2 -678.3782 -691.6535 -651.8682 -709.7373 -659.3102 -719.7595 -695.5856
## 3 -684.8992 -672.6377 -677.8125 -669.9904 -654.0443 -696.2913 -685.1682
## 4 -685.5860 -669.0659 -689.8244 -646.4406 -688.4053 -696.4109 -695.3764
## 5 -695.5070 -669.6889 -693.1522 -664.2876 -694.5112 -683.8182 -708.9387
## 6 -705.1183 -667.5226 -708.6766 -679.6809 -667.3445 -712.4262 -723.7749
## 7 -658.4365        NA -670.1641        NA        NA        NA -697.2630
## 8 -648.0953        NA -646.6228        NA        NA        NA -659.2519
## 9 -661.2521        NA -650.3789        NA        NA        NA -682.3787
##         VEE       EVE       VVE       EEV       VEV       EVV       VVV
## 1 -688.7975 -688.7975 -688.7975 -688.7975 -688.7975 -688.7975 -688.7975
## 2 -683.4840 -633.3132 -683.4810 -664.2698 -717.8221 -742.4286 -715.2660
## 3        NA        NA        NA -729.0026 -739.9755        NA        NA
## 4        NA        NA        NA -834.7233 -764.4534        NA        NA
## 5        NA        NA        NA -858.6840 -797.5165        NA        NA
## 6        NA        NA        NA -801.6111 -740.7355        NA        NA
## 7        NA        NA        NA -889.6800        NA        NA        NA
## 8        NA        NA        NA        NA        NA        NA        NA
## 9        NA        NA        NA        NA        NA        NA        NA
## 
## Top 3 models based on the BIC criterion: 
##     EVE,2     VEI,4     EEI,8 
## -633.3132 -646.4406 -646.6228

# Menampilkan hasil BIC untuk berbagai jumlah cluster
plot(gmm_model, what = "BIC")

optimal_clusters <- gmm_model$G
cat("Optimal number of clusters: ", optimal_clusters, "\n")

## Optimal number of clusters:  2

Juumlah culster optimal adalah 2, ini bisa dilihat dengan mncari titik dengan nilai BIC terbesar. Dalam kasus ini niali BIC terbesar adalah model EVE dengan jumlah komponen 2.

Hasil Clustering

Summary Jumlah Cluster optimal

gmm_model2 = Mclust(dt, G = 2, modelNames = c("EVE"))
summary(gmm_model2, parameters = TRUE)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 2 components: 
## 
##  log-likelihood  n df       BIC       ICL
##       -230.2607 34 49 -633.3132 -633.6239
## 
## Clustering table:
##  1  2 
##  8 26 
## 
## Mixing probabilities:
##         1         2 
## 0.2371723 0.7628277 
## 
## Means:
##          [,1]        [,2]
## X1  0.1144555 -0.03558557
## X2  1.0713155 -0.33308480
## X3 -0.3672036  0.11416799
## X4  1.0983322 -0.34148459
## X5  0.6390915 -0.19870118
## X6 -0.3399775  0.10570305
## X7  0.1613344 -0.05016080
## 
## Variances:
## [,,1]
##           X1         X2         X3         X4         X5         X6        X7
## X1 1.1966799  0.5396505  0.3526419  0.7200114  0.2609832  1.0641955 1.0539213
## X2 0.5396505  1.6567542  0.2596593  0.6619644 -1.0051928 -0.5272227 0.1763555
## X3 0.3526419  0.2596593  2.4694787  0.3673152 -1.1227186  1.2063338 0.4237561
## X4 0.7200114  0.6619644  0.3673152  1.0466275 -0.2464518  0.3670720 0.6283520
## X5 0.2609832 -1.0051928 -1.1227186 -0.2464518  1.7777280  0.9367601 0.5730257
## X6 1.0641955 -0.5272227  1.2063338  0.3670720  0.9367601  2.7662602 1.6061668
## X7 1.0539213  0.1763555  0.4237561  0.6283520  0.5730257  1.6061668 1.2603627
## [,,2]
##             X1          X2           X3           X4         X5          X6
## X1  0.60395032 -0.02468338  0.074122367  0.083040279 0.04273550  0.17955786
## X2 -0.02468338  0.32662187 -0.053175085  0.182971337 0.11459858 -0.02550871
## X3  0.07412237 -0.05317508  0.549237732 -0.007011136 0.16691270  0.19222887
## X4  0.08304028  0.18297134 -0.007011136  0.201864312 0.07161786  0.19500003
## X5  0.04273550  0.11459858  0.166912696  0.071617863 0.44924615  0.01653687
## X6  0.17955786 -0.02550871  0.192228868  0.195000033 0.01653687  0.83416164
## X7  0.10548568 -0.02111411  0.129677631  0.082658961 0.02383007  0.17123913
##             X7
## X1  0.10548568
## X2 -0.02111411
## X3  0.12967763
## X4  0.08265896
## X5  0.02383007
## X6  0.17123913
## X7  0.71452709

Proporsi cluster

# Menghitung frekuensi observasi di setiap kluster
cluster_frequencies <- table(gmm_model2$classification)

# Mengurutkan kluster berdasarkan frekuensinya
sorted_clusters <- names(sort(cluster_frequencies, decreasing = TRUE))

# Membuat urutan kluster yang diinginkan (1, 2, 3, 4, 5)
new_order <- 1:length(sorted_clusters)

# Menukar isi kluster dengan urutan yang dihasilkan
gmm_model2$classification <- recode(gmm_model2$classification, 
                               !!!setNames(as.character(new_order), sorted_clusters))

table(gmm_model2$classification)

## 
##  1  2 
## 26  8

Proporsi jumlah cluster 1 dan 2 adalah 26 : 8.

Plot culster

colors <- c("1" = "#d3a83c", "2" = "#4fbaa9")

chart <-
fviz_cluster(gmm_model2, data = dt, repel = TRUE, labelsize =8) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)
chart

# Eksport
ggsave("GMM 2D.png", chart, path = export.chart,
        dpi = 300, height = 5, width = 7)

Terlihat bahwa GMM tidak membuat cluster berdararkan kedekatan titik, melainkan berdasarkan sebaran berbagai sebaran gaussian.

Cluster Profil

data.clust1 <- cbind(dt %>% as.data.frame() %>% mutate(across(everything(), as.numeric)), 
                     Cluster = gmm_model2[["classification"]]) 

# Calculate the mean of each variable for each cluster
cluster_profiles1 <- aggregate(. ~ Cluster, data.clust1, mean)

# Print the cluster profiles
print(cluster_profiles1)

##   Cluster          X1         X2         X3         X4         X5         X6
## 1       1 -0.04121201 -0.3346306  0.1121178 -0.3430521 -0.1977580  0.1031076
## 2       2  0.13393902  1.0875495 -0.3643828  1.1149194  0.6427135 -0.3350998
##            X7
## 1 -0.05254404
## 2  0.17076813

# Convert the data to long format for plotting
cluster_profiles_long1 <- tidyr::pivot_longer(cluster_profiles1, -Cluster, 
                                             names_to = "Variable", values_to = "Value")

# Create the bar plot
chart <-
ggplot(cluster_profiles_long1, aes(x = Cluster, y = Value, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "Cluster", y = "Mean Value", fill = "Variable") +
  theme_minimal() +
  ggtitle("Cluster Profiles") +
  scale_fill_brewer(palette = "GnBu")

chart

# Eksport
ggsave("GMM Prof.png", chart, path = export.chart,
        dpi = 300, height = 4, width = 7)

Terlihat dari nilai yang sudah distandarisasi bahwa profil dari cuslter 2 memilki nilai yang hampir unggul dalam segala peubah dibandingan dengan cluster 1. Notes: adanya nilai negatif karena hasil dari standarisasi memang mengandung nilai positif dan negatif. Agar mudah perbandingannya, mari kita lihat dengan diagram spider web dibawah ini.

data.akhir1 <- cbind(raw.data, Cluster =  gmm_model2[["classification"]]) %>% 
  relocate(Cluster, .before = 2)

# Radar Plot
chart <- 
ggRadar(
  data = data.akhir1,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)

chart

# Eksport
ggsave("GMM Radar.png", chart, path = export.chart,
        dpi = 300, height = 10, width = 10)

Selaras dengan yang sudah saya sampaikan diatas, bahwa cluster 2 hampir unggul dalam segala aspek. Cluster 1 hanya bisa dunggul dalam sanitasi dan imun saja, sedangkan yang lainnya tidak mampu mengungguli atau bahkan tidak ada harapan untuk mengungguli cluster 2.

Map

indo <- st_read(dsn= "C:/Users/Fathan/Documents/Obsidian Vault/2. Kuliah/Smt 5/8. Pengantar Sains Data/Tugas/Tugas Akhir/SHP Indonesia/prov.shp", 
                quiet = TRUE)

data.map <- cbind(data.clust1, KODE=kode$KODE)
  
data.indo <- indo %>%
  inner_join(data.map, by = c("KODE" = "KODE"))

chart <-
ggplot() +  
  geom_sf(data=data.indo, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "#e0c377", "2" = "#80cdc0"), 
                    name = "Keterangan") +
  labs(title = "GMM Clustering - Kesejahteraan \n pada Provinsi Indonesia 2021",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5)) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))
chart

# Eksport
ggsave("GMM Map.png", chart, path = export.chart,
        dpi = 300, height = 9, width = 16)

Ada yang aneh? ya memang. Kenapa bisa provinsi papua masuk kedalam cluster 2 yang sebelumnya disebut-sebut sebagai cluster yang unggul dalam segala aspek kesejahteraan. Mari kita telurusi lebih lanjut.

Anggota Cluster

Tampilan Data dengan cluster

datatable(data.akhir1, options = list(pageLength = 5))

#Export data
write.csv(data.akhir1, "data_clust1.csv", row.names = FALSE)

Jika anda ingin melihat bagaimana tampilan data dengan cluster nya. Atau lebih mudah dengan bar chart dibawah.

colors1 <- c("1" = "#e7dec5", "2" = "#7dcdc0")

# List of variable names (excluding Provinsi and Cluster)
variables <- names(data.akhir1)[!names(data.akhir1) %in% c("Provinsi", "Cluster")]

plots <- list()

for (var in variables) {
  # Sort data based on the current variable
  sorted_data <- data.akhir1 %>%
    arrange(!!sym(var)) %>%
    mutate(Provinsi = factor(Provinsi, levels = Provinsi))
  
  p <- ggplot(sorted_data, aes_string(x = var, y = "Provinsi", fill = "factor(Cluster)")) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = colors1) +
    labs(x = var, y = "Provinsi", fill = "Cluster") +
    theme_minimal()
  
  plots[[var]] <- p
}

combined_plot <- plot_grid(plotlist = plots, ncol = 4)
combined_plot

# Eksport
ggsave("GMM Dist.png", combined_plot, path = export.chart,
        dpi = 300, height = 15, width = 30)

Memang terlihat keanehannya, bahwa salah satu anggota cluster 2 tidak bisa meungguli semua aspek kesejahteraan yang ada. Bahwakan kalah jauh daripada yang lain. Sepertinya memang clustering dengan GMM bukanlah ide yang tepat.

K-Means

n Cluster Optimum

Menentukan jumlah klaster paling optimum untuk membagi 34 provinsi yang ada di Indonesia. Pada dasarnya nilai k optimum dapat ditentukan sendiri oleh peneliti, selain itu nilai k dapat ditentukan juga menggunakan bantuan grafik Elbow berikut. Dimana k minimum akan didapat ketika grafik sudah tidak bergerak menurun dengan tajam.

Silhouette

fviz_nbclust(x = dt, 
             FUNcluster = kmeans,
             method = 'silhouette',
             k.max = 10)

Dari hasil diatas didapatkan nilai K optimum adalah 2.

Gap Statistic

#calculate gap statistic for each number of clusters (up to 10 clusters)
gap_stat <- clusGap(dt, FUN = hcut, nstart = 25, K.max = 10, B = 50)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat)

Dari hasil diatas didapatkan nilai K optimum adalah 1.

WSS (Elbow)

fviz_nbclust(dt, kmeans, method = "wss")

Untuk wss nampaknya agak sulit untuk melihat berapa jumlah cluster yang optimumnya, namun nampaknya 6 atau 2 merupakan jumlah cluster optimumnya.

Hubert Statistic & Dindex

nb <- NbClust(data = dt, distance = "euclidean", method="complete")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 8 proposed 2 as the best number of clusters 
## * 5 proposed 3 as the best number of clusters 
## * 2 proposed 4 as the best number of clusters 
## * 3 proposed 5 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 4 proposed 14 as the best number of clusters 
## * 1 proposed 15 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  2 
##  
##  
## *******************************************************************

Untuk hubert statistic dan dindex, jumlah cluster optimal yang disarankan adalah 2.

Hasil Clustering 1

Beberapa metode untuk menentukan jumlah cluster memiliki jumlah cluster optimal yang cukup beragam. Namun 2 merupakan jumlah cluster yang nampaknya bisa diterima oleh semua metode yang ada. Sehingga jumlah cluster yang akan dicobakan adalah 2.

Mari kita tambahkan paramter nstart = 6 untuk menentukan jumlah awalisasi acak yang akan dilakukan oleh algoritma K-means. K-means adalah algoritma yang sensitif terhadap titik awal centroid yang dipilih secara acak. Untuk mengurangi kemungkinan menemukan solusi yang hanya merupakan minimum lokal (bukan global), sering kali beberapa awalizasi acak dilakukan, dan hasil dengan total inersia (atau sum of squared distances from points to their assigned cluster centroids) terkecil dipilih sebagai hasil akhir.

set.seed(100)

data_kmeans <- eclust(dt, "kmeans", k = 2,
                      nstart = 6, graph = FALSE)

chart <-
fviz_cluster(object = data_kmeans, data=dt)  + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)
chart

# Eksport
ggsave("Kmeans 2D 1.png", chart, path = export.chart,
        dpi = 300, height = 5, width = 7)

Proporsi cluster

table(data_kmeans$cluster)

## 
##  1  2 
## 30  4

Proporsi jumlah cluster 1 dan 2 adalah 30 : 4.

Validasi Cluster

Silhouette

chart <- 
fviz_silhouette(data_kmeans, palette = "jco",
                ggtheme = theme_minimal())

##   cluster size ave.sil.width
## 1       1   30          0.39
## 2       2    4          0.41

chart

# Eksport
ggsave("Kmeans Valid.png", chart, path = export.chart,
        dpi = 300, height = 3, width = 7)

Terlihat bahwa sudah tidak ada pengamatan dengan niali siluet negatif, sehingga seluruh pengamatan dikelompokkan secara tidak benar sudah kemunkinannya sudah kecil.

Cluster Profil

data.clust2 <- cbind(dt %>% as.data.frame() %>% mutate(across(everything(), as.numeric)), 
                     Cluster = data_kmeans$cluster) 

# Calculate the mean of each variable for each cluster
cluster_profiles2 <- aggregate(. ~ Cluster, data.clust2, mean)

# Print the cluster profiles
print(cluster_profiles2)

##   Cluster         X1         X2          X3         X4           X5          X6
## 1       1 -0.1524433 -0.2679779 -0.04791201 -0.3412296 -0.009305985 -0.02421519
## 2       2  1.1433248  2.0098344  0.35934005  2.5592220  0.069794888  0.18161395
##          X7
## 1 -0.138005
## 2  1.035037

# Convert the data to long format for plotting
cluster_profiles_long2 <- tidyr::pivot_longer(cluster_profiles2, -Cluster, 
                                             names_to = "Variable", values_to = "Value")

# Create the bar plot
chart <- 
ggplot(cluster_profiles_long2, aes(x = Cluster, y = Value, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "Cluster", y = "Mean Value", fill = "Variable") +
  theme_minimal() +
  ggtitle("Cluster Profiles") +
  scale_fill_brewer(palette = "GnBu")
chart

# Eksport
ggsave("Kmeans Prof.png", chart, path = export.chart,
        dpi = 300, height = 4, width = 7) 

Berlainan dengan hasil dari GMM, cluster 2 merupakan cluster dengan aspek kesejahteraan yang paling tinggi dibandingkan dengan cluster 1. Benar-benar unggul semua. Mari kita lihat lebih jelas dengan radar chart dibawah.

data.akhir2 <- cbind(raw.data, Cluster =  data_kmeans$cluster) %>% 
  relocate(Cluster, .before = 2)

colors <- c("1" = "#e0c377", "2" = "#7dcdc0")

# Radar Plot
chart <- 
ggRadar(
  data = data.akhir2,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)

chart

# Eksport
ggsave("Kmeans Radar.png", chart, path = export.chart,
        dpi = 300, height = 10, width = 10)

Sejalan dengan yang saya bilang diatas, bahwa cluster 2 benar-benar mengungguli cluster 1.

Map

data.map <- cbind(data.clust2, KODE=kode$KODE)
  
data.indo <- indo %>%
  inner_join(data.map, by = c("KODE" = "KODE"))

chart <- 
ggplot() +  
  geom_sf(data=data.indo, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "#e0c377", "2" = "#80cdc0"), 
                    name = "Keterangan") +
  labs(title = "K-means Clustering -  Kesejahteraan \n pada Provinsi Indonesia 2021",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5)) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))
chart

# Eksport
ggsave("Kmeans Map.png", chart, path = export.chart,
        dpi = 300, height = 9, width = 16)

Berbeda dengan GMM, hasil dari k-means itu lebih make sense, mari kita lihat keanggotaan clusternya.

Anggota Cluster

Tampilan Data dengan cluster

datatable(data.akhir2, options = list(pageLength = 5))

#Export data
write.csv(data.akhir2, "data_clust2.csv", row.names = FALSE)

# List of variable names (excluding Provinsi and Cluster)
variables <- names(data.akhir2)[!names(data.akhir2) %in% c("Provinsi", "Cluster")]

plots <- list()

for (var in variables) {
  # Sort data based on the current variable
  sorted_data <- data.akhir2 %>%
    arrange(!!sym(var)) %>%
    mutate(Provinsi = factor(Provinsi, levels = Provinsi))
  
  p <- ggplot(sorted_data, aes_string(x = var, y = "Provinsi", fill = "factor(Cluster)")) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = colors1) +
    labs(x = var, y = "Provinsi", fill = "Cluster") +
    theme_minimal()
  
  plots[[var]] <- p
}

combined_plot <- plot_grid(plotlist = plots, ncol = 4)
combined_plot

# Eksport
ggsave("Kmeans Dist.png", combined_plot, path = export.chart,
        dpi = 300, height = 15, width = 30)

Bisa dilihat bahwa memang K-Means menghasilkan cluster yang lebih make sense daripada GMM. Karena hampir di semua aspek memang cluster 1 lebih unggul daripada cluster 2.

Fuzzy C-Means

Untuk FCM pemilihan jumlah clusternya sama dengan K-Means, sehingga jumlah cluster yang akan digunakan adalah 2.

Hasil Clustering 1

# clustering menggunakan fanny()
data_fcm <- fanny(x = dt, k = 2, metric = "euclidean")

chart <-
fviz_cluster(object = data_fcm, data=dt) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)
chart

# Eksport
ggsave("FCM 2D 1.png", chart, path = export.chart,
        dpi = 300, height = 5, width = 7) 

Proporsi Cluster

table(data_fcm$cluster)

## 
##  1  2 
## 16 18

Berbeda dengan 2 metode sebelumnya, FCM menghasilkan proporsi jumlah cluster yang berbeda, yakni lebih proporsional antar clusternya yakni 16 : 18.

Validasi Cluster

Dunn Index

data_fcm$coeff

##   dunn_coeff   normalized 
## 5.000000e-01 1.998401e-15

Nilai Dunn index yang kita peroleh sebesar 0.5, sangat kecil tbh. Untuk mengetahui hasil cluster yang kita bentuk sudah baik berdasarkan nilai Dunn index yaitu dengan mendapatkan Dunn index yang semakin besar.

Kita perlu mengetahui berapa banyak cluster yang terbentuk dari nilai derajat keanggotaan yang ada.

Pengelompokkan

data_fcm$clustering %>% unique()

## [1] 1 2

Dari hasil pengelompokan menjadi 2 kelompok, ternyata 2 kelompok tersebut berhasil terisi. Artinya bahwa dengan menggunakan k=2, fuzzy c-means mampu menseparasikan data berdasarkan 2 kelompok yang berbeda. Sehingga juumlah cluster 2 merupakan jumlah cluster yang sudah tepat.

Silhouette

chart <- 
fviz_silhouette(data_fcm, palette = "jco",
                ggtheme = theme_minimal())

##   cluster size ave.sil.width
## 1       1   16          0.03
## 2       2   18          0.31

chart

# Eksport
ggsave("FCM Valid 1.png", chart, path = export.chart,
        dpi = 300, height = 3, width = 7)

Karena pengamatan dengan nilai siluet negatif yang dapat dilihat pada cluster 1 yang berwarna biru pada plot mengindikasikan bahwa pengamatan tersebut mungkin dikelompokkan secara tidak benar. Mari lihat angka indeks dari pengamatan tersebut.

sil <- data_fcm$silinfo$widths
neg_sil_index <- which(sil[, 'sil_width']<0)
data[neg_sil_index, 1]

## [1] "KEPULAUAN RIAU" "DKI JAKARTA"    "JAWA BARAT"     "JAWA TENGAH"   
## [5] "D I YOGYAKARTA" "JAWA TIMUR"     "BANTEN"

Terlihat bahwa 7 provinsi tersebut merupakan provinsi yang mungkin blm dapat dikelompokkan secara benar oleh FCM. Mari kita perbaiki.

Hasil Cluster 2

data_fcm <- fanny(x = dt, k = 2, metric = "manhattan", memb.exp = 3)

chart <-
fviz_cluster(object = data_fcm, data=dt) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)
chart

# Eksport
ggsave("FCM 2D 2.png", chart, path = export.chart,
        dpi = 300, height = 5, width = 7)  

Proporsi Cluster

table(data_fcm$cluster)

## 
##  1  2 
## 19 15

Berbeda dengan 2 metode sebelumnya, FCM menghasilkan proporsi jumlah cluster yang berbeda, yakni lebih proporsional antar clusternya yakni 19 : 15.

Validasi Cluster

Dunn Index

data_fcm$coeff

## dunn_coeff normalized 
##        0.5        0.0

Nilai Dunn index yang kita peroleh masih sebesar 0.5.

Dari hasil pengelompokan menjadi 2 kelompok, ternyata 2 kelompok tersebut berhasil terisi. Artinya bahwa dengan menggunakan k=2, fuzzy c-means mampu menseparasikan data berdasarkan 2 kelompok yang berbeda. Sehingga juumlah cluster 2 merupakan jumlah cluster yang sudah tepat.

Kita perlu mengetahui berapa banyak cluster yang terbentuk dari nilai derajat keanggotaan yang ada.

Pengelompokkan

data_fcm$clustering %>% unique()

## [1] 1 2

Dari hasil pengelompokan menjadi 2 kelompok, ternyata 2 kelompok tersebut berhasil terisi. Artinya bahwa dengan menggunakan k=2, fuzzy c-means mampu menseparasikan data berdasarkan 2 kelompok yang berbeda. Sehingga juumlah cluster 2 merupakan jumlah cluster yang sudah tepat.

chart <- 
fviz_silhouette(data_fcm, palette = "jco",
                ggtheme = theme_minimal())

##   cluster size ave.sil.width
## 1       1   19          0.27
## 2       2   15          0.14

chart

# Eksport
ggsave("FCM Valid 2.png", chart, path = export.chart,
        dpi = 300, height = 3, width = 7)

Masih ada pengamatan dengan nilai siluet negatif yang dapat dilihat pada cluster 1 yang berwarna kuning pada plot mengindikasikan bahwa pengamatan tersebut mungkin dikelompokkan secara tidak benar. Mari lihat angka indeks dari pengamatan tersebut.

sil <- data_fcm$silinfo$widths
neg_sil_index <- which(sil[, 'sil_width']<0)
data[neg_sil_index, 1]

## [1] "MALUKU UTARA" "PAPUA BARAT"  "PAPUA"

Terlihat bahwa 3 provinsi tersebut merupakan provinsi yang mungkin blm dapat dikelompokkan secara benar oleh FCM. Namun saya rasa, saya sudah mencoba berbagai hal untuk memperbaikinya, dan ini adalah hasil terbaik dari FCM. Jadi ya, kita sudahi.

Cluster Profil

data.clust3 <- cbind(dt %>% as.data.frame() %>% mutate(across(everything(), as.numeric)), 
                     Cluster = data_fcm$cluster) 

# Calculate the mean of each variable for each cluster
cluster_profiles3 <- aggregate(. ~ Cluster, data.clust3, mean)

# Print the cluster profiles
print(cluster_profiles3)

##   Cluster         X1         X2         X3         X4         X5         X6
## 1       1 -0.4986192 -0.1677914 -0.5238502 -0.3578339 -0.1385835 -0.4146049
## 2       2  0.6315843  0.2125358  0.6635436  0.4532563  0.1755392  0.5251663
##           X7
## 1 -0.5125450
## 2  0.6492236

# Convert the data to long format for plotting
cluster_profiles_long3 <- tidyr::pivot_longer(cluster_profiles3, -Cluster, 
                                             names_to = "Variable", values_to = "Value")

# Create the bar plot
chart <-
ggplot(cluster_profiles_long3, aes(x = Cluster, y = Value, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "Cluster", y = "Mean Value", fill = "Variable") +
  theme_minimal() +
  ggtitle("Cluster Profiles")  +
  scale_fill_brewer(palette = "GnBu")
chart

# Eksport
ggsave("FCM Prof.png", chart, path = export.chart,
        dpi = 300, height = 4, width = 7) 

Dari sini kita tahu bahwa cluster 1 dan cluster 2 cenderung memiliki bentuk bar yang mirip namun kebalik. Tidak tahu pasti apa artinya, namun cluster 2 merupakan cluster yang lebih unggul dalam segala aspek daripada cluster 1. Mari kita lihat radar chart nya.

data.akhir3 <- cbind(raw.data, Cluster =  data_fcm$cluster) %>% 
  relocate(Cluster, .before = 2)

# Radar Plot
chart <-
ggRadar(
  data = data.akhir3,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
) + 
  scale_colour_manual(values = colors) +
  scale_fill_manual(values = colors)
chart

# Eksport
ggsave("FCM Radar.png", chart, path = export.chart,
        dpi = 300, height = 10, width = 10)

Selaras dengan yang saya sampaikan di atas, bawahwa cluster 1 dan cluster 2 memang mirip, namun cluster 2 memiliki nilai yang lebih besar daripada cluster 1.

Map

data.map <- cbind(data.clust3, KODE=kode$KODE)
  
data.indo <- indo %>%
  inner_join(data.map, by = c("KODE" = "KODE"))

chart <-
ggplot() +  
  geom_sf(data=data.indo, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "#e0c377", "2" = "#80cdc0"), 
                    name = "Keterangan") +
  labs(title = "FCM Clustering -  Kesejahteraan \n pada Provinsi Indonesia 2021",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5)) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))
chart

# Eksport
ggsave("FCM Map.png", chart, path = export.chart,
        dpi = 300, height = 9, width = 16)

Beginilah penampakan dari hasil clustering FCM. Apakah lebih bagus daripada metode yang lain? Mari kita pastikan dengan melihat keanggotaan clusternya.

Anggota Cluster

Tampilan Data dengan cluster

datatable(data.akhir3, options = list(pageLength = 5))

#Export data
write.csv(data.akhir3, "data_clust1.csv", row.names = FALSE)

# List of variable names (excluding Provinsi and Cluster)
variables <- names(data.akhir3)[!names(data.akhir3) %in% c("Provinsi", "Cluster")]

plots <- list()

for (var in variables) {
  # Sort data based on the current variable
  sorted_data <- data.akhir3 %>%
    arrange(!!sym(var)) %>%
    mutate(Provinsi = factor(Provinsi, levels = Provinsi))
  
  p <- ggplot(sorted_data, aes_string(x = var, y = "Provinsi", fill = "factor(Cluster)")) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = colors1) +
    labs(x = var, y = "Provinsi", fill = "Cluster") +
    theme_minimal()
  
  plots[[var]] <- p
}

combined_plot <- plot_grid(plotlist = plots, ncol = 4)
combined_plot

# Eksport
ggsave("FCM Dist.png", combined_plot, path = export.chart,
        dpi = 300, height = 15, width = 30)

Walaupun dalam validasi cluster dengan metode siluet masih ada yang bernilai negatif, namun saya rasa metode FCM mampu membuat cluster yang jauh lebih baik daripada yang lainnya. Tidak seperti K-means yang pada profil nya seolah berkata bahwa cluster 2 itu jauh lebih baik daripada cluster 1 di segala aspek. Namun pada saat kita tampilkan sebaran tiap peubah pada semua provinsi nya, cluster 2 tidak selalu berada di puncak.

Sedangkan untuk FCM itu sesuai dengan apa yang dikatakan oleh profil nya. Yakni clsuter 1 dan cluster 2 itu mirip, namun cluster 2 lebih tinggi daripada cluster 1. Ini juga sesuai dengan sebaran yang kita lihat di atas.

Referensi

• Validasi Cluster

• Clustering

• https://rpubs.com/inayatus/fuzzy-c-means

• https://rpubs.com/Anoe/lbb-kmeans

• https://blog.datawrapper.de/