Indikator IPM jawa barat

library(readxl)
df_jabar<-read_excel("ipm_jabar.xlsx")
head(df_jabar)

# A tibble: 6 × 6
  Provinsi   KOTA            AHH   RLS   PPK   IPM
  <chr>      <chr>         <dbl> <dbl> <dbl> <dbl>
1 Jawa Barat BANDUNG        74.9  9.1  11018  73.7
2 Jawa Barat BANDUNG BARAT  74.8  8.23  9392  69.6
3 Jawa Barat BEKASI         75.1  9.57 12123  75.8
4 Jawa Barat BOGOR          74.7  8.37 11153  71.8
5 Jawa Barat CIAMIS         75.0  8.09  9750  72.0
6 Jawa Barat CIANJUR        74.6  7.22  8626  66.6

library(tidyverse)
library(dplyr)
  print(colnames(df_jabar))

[1] "Provinsi" "KOTA"     "AHH"      "RLS"      "PPK"      "IPM"     

  df_jabar2 <- df_jabar %>%
  dplyr::select(-one_of("Provinsi")) %>%
  column_to_rownames(var = "KOTA")
  head(df_jabar2,5)

                AHH  RLS   PPK   IPM
BANDUNG       74.92 9.10 11018 73.74
BANDUNG BARAT 74.78 8.23  9392 69.61
BEKASI        75.09 9.57 12123 75.76
BOGOR         74.67 8.37 11153 71.78
CIAMIS        74.96 8.09  9750 72.05

Uji normalitas

library(tidyverse)
normalitas <- df_jabar2 %>%
  gather(key = "Indicator", value = "Value") %>%
  group_by(Indicator) %>%
  summarise(
    Shapiro_Wilk_Statistic = round(shapiro.test(Value)$statistic,3),
    p_value = shapiro.test(Value)$p.value
  )
print(normalitas)

# A tibble: 4 × 3
  Indicator Shapiro_Wilk_Statistic p_value
  <chr>                      <dbl>   <dbl>
1 AHH                        0.923 0.0475 
2 IPM                        0.919 0.0369 
3 PPK                        0.875 0.00370
4 RLS                        0.893 0.00930

library(ggplot2)
library(gridExtra)
library(grid)

# Function to Create Blue Histogram with Gridlines
plot_histogram <- function(data, variable_name) {
  ggplot(data, aes_string(x = variable_name)) +
    geom_histogram(fill = "skyblue", color = "white", bins = 30) +  # Adjust bins as needed
    labs(x = variable_name, y = "Frequency") +
    theme_minimal() +  # Clean theme
    theme(panel.grid.major.y = element_line(color = "lightgray"))  # Add y-axis gridlines
}

# Create Histograms for each variable
a <- plot_histogram(df_jabar, "AHH")
b <- plot_histogram(df_jabar, "RLS")
d <- plot_histogram(df_jabar, "PPK")
e <- plot_histogram(df_jabar, "IPM")

# Combine the plots into one frame with a title
grid.arrange(
  arrangeGrob(a, b, d, e, ncol = 2),
  top = textGrob("", gp = gpar(fontsize = 16, fontface = "bold"))
)

library(tidyverse)
library(dplyr)

# Membuat kolom status IPM berdasarkan rentang IPM yang ditentukan
df_jabar <- df_jabar %>%
  mutate(rentang_IPM = case_when(
    IPM >= 80 ~ "Sangat Tinggi",
    IPM >= 70 & IPM < 80 ~ "Tinggi",
    IPM >= 60 & IPM < 70 ~ "Sedang",
    IPM < 60 ~ "Rendah"
  ))

if ("IPM_Kategori" %in% names(df_jabar)) {
  df_jabar <- df_jabar[, !names(df_jabar) %in% "IPM_Kategori"]
}

df_jabar <- df_jabar %>%
  arrange(desc(IPM))

head(df_jabar,5)

# A tibble: 5 × 7
  Provinsi   KOTA           AHH   RLS   PPK   IPM rentang_IPM  
  <chr>      <chr>        <dbl> <dbl> <dbl> <dbl> <chr>        
1 Jawa Barat KOTA BANDUNG  75.5  11.1 18236  83.0 Sangat Tinggi
2 Jawa Barat KOTA BEKASI   75.9  11.7 16479  83.0 Sangat Tinggi
3 Jawa Barat KOTA DEPOK    75.5  11.6 16279  82.4 Sangat Tinggi
4 Jawa Barat KOTA CIMAHI   75.3  11.4 12883  79.5 Tinggi       
5 Jawa Barat KOTA BOGOR    75.5  10.6 12656  77.8 Tinggi       

library(tidyverse)
library(dplyr)

df_jabar <- df_jabar %>%
  mutate(rentang_IPM = case_when(
    IPM >= 80 ~ "Sangat Tinggi",
    IPM >= 70 & IPM < 80 ~ "Tinggi",
    IPM >= 60 & IPM < 70 ~ "Sedang",
    IPM < 60 ~ "Rendah"
  ))

df_jabar <- df_jabar %>%
  arrange(desc(IPM))
ggplot(df_jabar, aes(x = rentang_IPM, fill = rentang_IPM)) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  scale_fill_manual(values = c("Sangat Tinggi" = "green", "Tinggi" = "yellow", "Sedang" = "orange", "Rendah" = "red")) +
  labs(title = "Distribusi Rentang IPM", x = "Rentang IPM", y = "Jumlah") +
  theme_minimal() +
  theme(legend.position = "none")

library(tidyverse)
library(dplyr)
library(corrplot)

numeric_df_jabar <- df_jabar %>%
  dplyr::select(-IPM) %>%
  dplyr::select_if(is.numeric)

# Menghitung matriks korelasi
cor_matrix <- cor(numeric_df_jabar, use = "complete.obs")

corrplot(cor_matrix, method = "color", col = colorRampPalette(c("skyblue", "white", "orange"))(200), 
         type = "upper", order = "hclust", tl.col = "black", tl.srt = 45, 
         addCoef.col = "black", number.cex = 0.7)

head(df_jabar)

# A tibble: 6 × 7
  Provinsi   KOTA           AHH   RLS   PPK   IPM rentang_IPM  
  <chr>      <chr>        <dbl> <dbl> <dbl> <dbl> <chr>        
1 Jawa Barat KOTA BANDUNG  75.5  11.1 18236  83.0 Sangat Tinggi
2 Jawa Barat KOTA BEKASI   75.9  11.7 16479  83.0 Sangat Tinggi
3 Jawa Barat KOTA DEPOK    75.5  11.6 16279  82.4 Sangat Tinggi
4 Jawa Barat KOTA CIMAHI   75.3  11.4 12883  79.5 Tinggi       
5 Jawa Barat KOTA BOGOR    75.5  10.6 12656  77.8 Tinggi       
6 Jawa Barat KOTA CIREBON  75.2  10.4 12506  76.5 Tinggi       

summary(df_jabar2)

      AHH             RLS              PPK             IPM       
 Min.   :73.87   Min.   : 6.940   Min.   : 8562   Min.   :66.55  
 1st Qu.:74.69   1st Qu.: 7.865   1st Qu.: 9880   1st Qu.:69.50  
 Median :74.89   Median : 8.230   Median :11136   Median :72.09  
 Mean   :74.93   Mean   : 8.870   Mean   :11526   Mean   :73.25  
 3rd Qu.:75.08   3rd Qu.: 9.970   3rd Qu.:12449   3rd Qu.:76.04  
 Max.   :75.86   Max.   :11.660   Max.   :18236   Max.   :83.04  

df_jabar2<-df_jabar2[,1:3]
scaled_jabar<-scale(df_jabar2)
scaled_jabar<-data.frame(scaled_jabar)

cek Pencilan

library(ggplot2)
library(tidyverse)
# Mengubah dataframe ke format long untuk ggplot
df_long <- scaled_jabar %>%
  pivot_longer(cols = c(AHH, RLS, PPK), 
               names_to = "Variable", 
               values_to = "Value")

# Membuat boxplot dengan ggplot
boxplot <- ggplot(df_long, aes(x = Variable, y = Value, fill = Variable)) +
  geom_boxplot(outlier.color = "red", outlier.shape = 16, outlier.size = 2) +
  labs(title = "Boxplot for AHH, RLS, PPK", 
       x = "Variable", 
       y = "Value") +
  theme_minimal() +
  scale_fill_brewer(palette = "Pastel1") +
  theme(legend.position = "none")

# Menampilkan boxplot
print(boxplot)

# Fungsi untuk menghitung outliers menggunakan IQR
calculate_outliers <- function(df, variable) {
  Q1 <- quantile(df[[variable]], 0.25, na.rm = TRUE)
  Q3 <- quantile(df[[variable]], 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower_bound <- Q1 - 1.5 * IQR
  upper_bound <- Q3 + 1.5 * IQR
  outliers <- df %>%
    filter(df[[variable]] < lower_bound | df[[variable]] > upper_bound) %>%
    mutate(Row = row_number()) %>%
  return(outliers)
}

variables <- c("AHH", "RLS", "PPK")

# Membuat dataframe untuk menyimpan hasil outliers
outliers_summary <- data.frame(
  Variable = character(),
  Jumlah_Pencilan = integer(),
  Baris = character(),
  stringsAsFactors = FALSE
)

for (var in variables) {
  outliers <- calculate_outliers(df_jabar, var)
  jumlah_pencilan <- nrow(outliers)
  baris <- if (jumlah_pencilan > 0) paste(outliers$Row, collapse = ", ") else "Tidak ada pencilan"
  outliers_summary <- rbind(outliers_summary, data.frame(
    Variable = var,
    Jumlah_Pencilan = jumlah_pencilan,
    Baris = baris,
    stringsAsFactors = FALSE
  ))
}
print(outliers_summary)

  Variable Jumlah_Pencilan              Baris
1      AHH               2               1, 2
2      RLS               0 Tidak ada pencilan
3      PPK               2               1, 2

Gerombol hirarki

head(scaled_jabar)

                      AHH        RLS        PPK
BANDUNG       -0.02954018  0.1579374 -0.2144463
BANDUNG BARAT -0.38973985 -0.4394781 -0.9010936
BEKASI         0.40784512  0.4806791  0.2521867
BOGOR         -0.67275387 -0.3433422 -0.1574368
CIAMIS         0.07337401 -0.5356139 -0.7499129
CIANJUR       -0.82712515 -1.1330294 -1.2245695

library(factoextra)

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

# Plot cluster results
p1 <- fviz_nbclust(scaled_jabar, FUN = hcut, method = "wss", 
                   k.max = 10) +
  ggtitle("(A) Elbow method")
p2 <- fviz_nbclust(scaled_jabar, FUN = hcut, method = "silhouette", 
                   k.max = 10) +
  ggtitle("(B) Silhouette method")
#p3 <- fviz_nbclust(scaled_jabar, FUN = hcut, method = "gap_stat", 
                   #k.max = 16) +
  #ggtitle("(C) Gap statistic")

# Display plots side by side
gridExtra::grid.arrange(p1, p2, nrow = 1)

dist_jabar <- dist(scaled_jabar, method = "euclidean")  
distance_jabar<-data.frame(dist_jabar)

Complete linkage

set.seed(123)
library(cluster)
hclust_complete<-hclust(dist_jabar, method = "complete")
#dend<-as.dendrogram(hclust_complete)

k = 3

k3_hclustcomplete<-cutree(hclust_complete, k=3)
table(k3_hclustcomplete)

k3_hclustcomplete
 1  2  3 
 8 16  3 

fviz_dend( hclust_complete, k = 3, horiz = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco", k_colors = "jco", cex = 0.1)

k = 4

k4_hclustcomplete<-cutree(hclust_complete, k=4)
table(k4_hclustcomplete)

k4_hclustcomplete
 1  2  3  4 
 8 15  3  1 

Average linkage

set.seed(123)
library(cluster)
hclust_average<-hclust(dist_jabar, method = "average")
#dend<-as.dendrogram(hclust_average)

k = 3

k3_hclustaverage<-cutree(hclust_average, k=3)
table(k3_hclustaverage)

k3_hclustaverage
 1  2  3 
23  3  1 

fviz_dend( hclust_average, k = 3, horiz = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco", k_colors = "jco", cex = 0.1)

k = 4

k4_hclustaverage<-cutree(hclust_average, k=4)
table(k4_hclustaverage)

k4_hclustaverage
 1  2  3  4 
 8 15  3  1 

Single linkage

set.seed(123)
library(cluster)
hclust_single<-hclust(dist_jabar, method = "single")
#dend<-as.dendrogram(hclust_single)

k = 3

k3_hclustsingle<-cutree(hclust_single, k=3)
table(k3_hclustsingle)

k3_hclustsingle
 1  2  3 
23  3  1 

fviz_dend( hclust_single, k = 4, horiz = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco", k_colors = "jco", cex = 0.1)

k = 4

k4_hclustsingle<-cutree(hclust_single, k=4)
table(k4_hclustsingle)

k4_hclustsingle
 1  2  3  4 
23  1  2  1 

Ward linkage

set.seed(123)
library(cluster)
hclust_ward<-hclust(dist_jabar, method = "ward")
#dend<-as.dendrogram(hclust_ward)

k = 3

k3_hclustward<-cutree(hclust_ward, k=3)
table(k3_hclustward)

k3_hclustward
 1  2  3 
10 11  6 

fviz_dend( hclust_ward, k = 3, horiz = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco", k_colors = "jco", cex = 0.1)

k = 4

k4_hclustward<-cutree(hclust_ward, k=4)
table(k4_hclustward)

k4_hclustward
 1  2  3  4 
10 11  3  3 

Koefisien aglomeratif

# Menghitung koefisien aglomeratif untuk metode 'complete'
agnes_complete <- agnes(dist_jabar, method = "complete")
coef_complete <- agnes_complete$ac

# Menghitung koefisien aglomeratif untuk metode 'single'
agnes_single <- agnes(dist_jabar, method = "single")
coef_single <- agnes_single$ac

# Menghitung koefisien aglomeratif untuk metode 'average'
agnes_average <- agnes(dist_jabar, method = "average")
coef_average <- agnes_average$ac

# Menghitung koefisien aglomeratif untuk metode 'ward'
agnes_ward <- agnes(dist_jabar, method = "ward")
coef_ward <- agnes_ward$ac

# Membuat dataframe dengan hasil koefisien aglomeratif
ac_jabar <- data.frame(
  Method = c("Complete", "Single", "Average", "Ward"),
  Agglomerative_Coefficient = c(round(coef_complete,3), round(coef_single,3), round(coef_average,3), round(coef_ward,3))
)

print(ac_jabar)

    Method Agglomerative_Coefficient
1 Complete                     0.901
2   Single                     0.717
3  Average                     0.841
4     Ward                     0.929

library(cluster)
library(fpc)
library(clusterSim)

# Fungsi untuk menghitung S_W dan S_B 
calculate_sw_sb_ratio_hierarchical <- function(data, cluster_labels) {
  # Total Sum of Squares (TSS)
  overall_mean <- colMeans(data)
  tss <- sum((data - overall_mean)^2)
  
  # Within-cluster Sum of Squares (S_W)
  sw <- 0
  unique_clusters <- unique(cluster_labels)
  cluster_centers <- sapply(unique_clusters, function(cl) {
    colMeans(data[cluster_labels == cl, ])
  })
  cluster_centers <- t(cluster_centers) # Transpose to match dimensions
  
  for (i in 1:nrow(data)) {
    cluster <- cluster_labels[i]
    center <- cluster_centers[cluster, ]
    sw <- sw + sum((data[i, ] - center)^2)
  }
  
  # Between-cluster Sum of Squares (S_B)
  sb <- 0
  for (k in unique_clusters) {
    center <- cluster_centers[k, ]
    num_points <- sum(cluster_labels == k)
    sb <- sb + num_points * sum((center - overall_mean)^2)
  }
  
  # Ratio S_W / S_B
  ratio <- sw / sb
  
  return(list(S_W = sw, S_B = sb, ratio = ratio))
}

calculate_internal_metrics <- function(data, cluster_labels) {
  dist_matrix <- dist(data)
  cluster_internal <- cluster.stats(dist_matrix, cluster_labels, silhouette = TRUE)
  cluster_db <- index.DB(data, cluster_labels, dist_matrix, centrotypes = "medoids", p = 2, q = 2)
  sw_sb <- calculate_sw_sb_ratio_hierarchical(data, cluster_labels)
  dunn_index <- cluster_internal$dunn
  silhouette_score <- cluster_internal$avg.silwidth
  calinski_harabasz_index <- cluster_internal$ch
  db_index <- cluster_db$DB
  swsb_ratio <- sw_sb$ratio
  
  return(list(Dunn_Index = dunn_index, Silhouette = silhouette_score, Calinski_Harabasz = calinski_harabasz_index, Davies_Bouldin = db_index, sw_sb=swsb_ratio))
}

run_hierarchical_clustering <- function(data, method) {
  # Melakukan clustering hierarkis
  dist_jabar <- dist(data)
  hierarchical_clustering <- hclust(dist_jabar, method = method)
  
  # Membuat dataframe untuk menyimpan hasil metrik internal
  internal_metrics_all <- data.frame(
    k = integer(),
    Dunn_Index = numeric(),
    Silhouette = numeric(),
    Calinski_Harabasz = numeric(),
    Davies_Bouldin = numeric(),
    sw_sb = numeric()
  )
  
  for (k in 3:4) {
    cluster_labels <- cutree(hierarchical_clustering, k = k)
    
    # Menghitung metrik validasi internal
    internal_metrics <- calculate_internal_metrics(data, cluster_labels)
    
    # Tambahkan hasil ke dataframe
    internal_metrics_all <- rbind(internal_metrics_all, data.frame(
      k = k,
      Dunn_Index = round(internal_metrics$Dunn_Index, 3),
      Silhouette = round(internal_metrics$Silhouette, 3),
      Calinski_Harabasz = round(internal_metrics$Calinski_Harabasz, 3),
      Davies_Bouldin = round(internal_metrics$Davies_Bouldin, 3),
      sw_sb = round(internal_metrics$sw_sb, 3)
    ))
  }
  
  return(internal_metrics_all)
}

# Data yang akan digunakan
data <- scaled_jabar

# Menjalankan clustering hierarkis dengan berbagai metode linkage
methods <- c("ward")
results <- list()

for (method in methods) {
  cat("Running hierarchical clustering with method:", method, "\n")
  result <- run_hierarchical_clustering(data, method)
  results[[method]] <- result
  print(result)
}

Running hierarchical clustering with method: ward 
  k Dunn_Index Silhouette Calinski_Harabasz Davies_Bouldin sw_sb
1 3      0.163      0.322            33.079          1.146 0.363
2 4      0.163      0.328            33.633          1.124 0.228

# Pilih k optimal berdasarkan Dunn Index (nilai tertinggi)
best_dunn_index <- result[which.max(result$Dunn_Index), ]

# Pilih k optimal berdasarkan Silhouette (nilai tertinggi)
best_silhouette <- result[which.max(result$Silhouette), ]

# Pilih k optimal berdasarkan Calinski-Harabasz (nilai tertinggi)
best_calinski_harabasz <- result[which.max(result$Calinski_Harabasz), ]

# Pilih k optimal berdasarkan Davies-Bouldin (nilai terendah)
best_davies_bouldin <- result[which.min(result$Davies_Bouldin), ]

# Pilih k optimal berdasarkan sw_sb (nilai terendah)
best_sw_sb <- result[which.min(result$sw_sb), ]

# Gabungkan hasil-hasil tersebut dalam satu dataframe
best_clusters <- data.frame(
  Metric = c("Dunn_Index", "Silhouette", "Calinski_Harabasz", "Davies_Bouldin", "sw_sb"),
  k = c(best_dunn_index$k, best_silhouette$k, best_calinski_harabasz$k, best_davies_bouldin$k, best_sw_sb$k),
  Score = c(best_dunn_index$Dunn_Index, best_silhouette$Silhouette, best_calinski_harabasz$Calinski_Harabasz, best_davies_bouldin$Davies_Bouldin, best_sw_sb$sw_sb)
)

print(best_clusters)

             Metric k  Score
1        Dunn_Index 3  0.163
2        Silhouette 4  0.328
3 Calinski_Harabasz 4 33.633
4    Davies_Bouldin 4  1.124
5             sw_sb 4  0.228

dari empat indeks validitas internal terpilih metode ward dengan empat gerombol optimum

dendogram

# Memuat perpustakaan yang diperlukan
library(dplyr)
library(heatmaply)
library(RColorBrewer)

# Menyiapkan data matriks untuk heatmap
data_matrix <- as.matrix(scaled_jabar)

# Membuat heatmap dan dendrogram dengan metode linkage ward.D2
heatmaply(
  data_matrix,
  colors = colorRampPalette(brewer.pal(3, "RdBu"))(256),
  k_col = 2, 
  k_row = 4, # Menggunakan 4 cluster pada baris seperti yang ditentukan
  scale = "none",
  main = "Heatmap dan Dendrogram Berdasarkan Cluster",
  xlab = "Indikator",
  ylab = "Kota/Kabupaten",
  dendrogram = "both",
  fontsize_row = 8,
  fontsize_col = 10,
  labRow = rownames(data_matrix),
  labCol = colnames(data_matrix),
  hclust_method = "ward.D2" # Menentukan metode linkage ward.D2
)

# Menyimpan heatmap ke file HTML dengan metode linkage ward.D2
heatmaply(
  data_matrix,
  colors = colorRampPalette(brewer.pal(3, "RdBu"))(256),
  k_col = 2, 
  k_row = 4, # Menggunakan 4 cluster pada baris seperti yang ditentukan
  scale = "none",
  main = "Heatmap dan Dendrogram Berdasarkan Cluster",
  xlab = "Indikator",
  ylab = "Kota/Kabupaten",
  dendrogram = "both",
  fontsize_row = 8,
  fontsize_col = 10,
  labRow = rownames(data_matrix),
  labCol = colnames(data_matrix),
  hclust_method = "ward.D2", # Menentukan metode linkage ward.D2
  #file = "heatmap_single_linkage_k4.html"
)

dendro_data_k <- function(hc, k) {
  
  hcdata    <-  ggdendro::dendro_data(hc, type = "rectangle")
  seg       <-  hcdata$segments
  labclust  <-  cutree(hc, k)[hc$order]
  segclust  <-  rep(0L, nrow(seg))
  heights   <-  sort(hc$height, decreasing = TRUE)
  height    <-  mean(c(heights[k], heights[k - 1L]), na.rm = TRUE)
  
  for (i in 1:k) {
    xi      <-  hcdata$labels$x[labclust == i]
    idx1    <-  seg$x    >= min(xi) & seg$x    <= max(xi)
    idx2    <-  seg$xend >= min(xi) & seg$xend <= max(xi)
    idx3    <-  seg$yend < height
    idx     <-  idx1 & idx2 & idx3
    segclust[idx] <- i
  }
  
  idx                    <-  which(segclust == 0L)
  segclust[idx]          <-  segclust[idx + 1L]
  hcdata$segments$clust  <-  segclust
  hcdata$segments$line   <-  as.integer(segclust < 1L)
  hcdata$labels$clust    <-  labclust
  
  hcdata
}

set_labels_params <- function(nbLabels,
                              direction = c("tb", "bt", "lr", "rl"),
                              fan       = FALSE) {
  if (fan) {
    angle       <-  360 / nbLabels * 1:nbLabels + 90
    idx         <-  angle >= 90 & angle <= 270
    angle[idx]  <-  angle[idx] + 180
    hjust       <-  rep(0, nbLabels)
    hjust[idx]  <-  1
  } else {
    angle       <-  rep(0, nbLabels)
    hjust       <-  0
    if (direction %in% c("tb", "bt")) { angle <- angle + 45 }
    if (direction %in% c("tb", "rl")) { hjust <- 1 }
  }
  list(angle = angle, hjust = hjust, vjust = 0.5)
}
plot_ggdendro <- function(hcdata,
                          direction   = c("lr", "rl", "tb", "bt"),
                          fan         = FALSE,
                          scale.color = NULL,
                          branch.size = 1,
                          label.size  = 3,
                          nudge.label = 0.01,
                          expand.y    = 0.1) {
  
  direction <- match.arg(direction) # if fan = FALSE
  ybreaks   <- pretty(segment(hcdata)$y, n = 5)
  ymax      <- max(segment(hcdata)$y)
  
  ## branches
  p <- ggplot() +
    geom_segment(data         =  segment(hcdata),
                 aes(x        =  x,
                     y        =  y,
                     xend     =  xend,
                     yend     =  yend,
                     linetype =  factor(line),
                     colour   =  factor(clust)),
                 lineend      =  "round",
                 show.legend  =  FALSE,
                 size         =  branch.size)
  
  ## orientation
  if (fan) {
    p <- p +
      coord_polar(direction = -1) +
      scale_x_continuous(breaks = NULL,
                         limits = c(0, nrow(label(hcdata)))) +
      scale_y_reverse(breaks = ybreaks)
  } else {
    p <- p + scale_x_continuous(breaks = NULL)
    if (direction %in% c("rl", "lr")) {
      p <- p + coord_flip()
    }
    if (direction %in% c("bt", "lr")) {
      p <- p + scale_y_reverse(breaks = ybreaks)
    } else {
      p <- p + scale_y_continuous(breaks = ybreaks)
      nudge.label <- -(nudge.label)
    }
  }
  
  # labels
  labelParams <- set_labels_params(nrow(hcdata$labels), direction, fan)
  hcdata$labels$angle <- labelParams$angle
  
  p <- p +
    geom_text(data        =  label(hcdata),
              aes(x       =  x,
                  y       =  y,
                  label   =  label,
                  colour  =  factor(clust),
                  angle   =  angle),
              vjust       =  labelParams$vjust,
              hjust       =  labelParams$hjust,
              nudge_y     =  ymax * nudge.label,
              size        =  label.size,
              show.legend =  FALSE)
  
  # colors and limits
  if (!is.null(scale.color)) {
    p <- p + scale_color_manual(values = scale.color)
  }
  
  ylim <- -round(ymax * expand.y, 1)
  p    <- p + expand_limits(y = ylim)
  
  p
}

library(ggdendro)
hcdata <- dendro_data_k(hclust_ward, 4)

p <- plot_ggdendro(hcdata,
                   direction   = "lr",
                   expand.y    = 0.2)
p

cols <- c("#a9a9a9","#bdd7e7", "#6baed6", "#3182bd",  "#08519c")

p <- plot_ggdendro(hcdata,
                   direction   = "tb",
                   scale.color = cols,
                   label.size  = 2.5,
                   branch.size = 0.5,
                   expand.y    = 0.2)

p <- p + theme_void()
p

Stabilitas gerombol

complete linkage

library(clValid)
stability_complete <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "complete", validation="stability")
summary(stability_complete)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                       3      4
                               
hierarchical APN  0.1821 0.2082
             AD   1.3268 1.0182
             ADM  0.7664 0.4581
             FOM  0.6187 0.5956

Optimal Scores:

    Score  Method       Clusters
APN 0.1821 hierarchical 3       
AD  1.0182 hierarchical 4       
ADM 0.4581 hierarchical 4       
FOM 0.5956 hierarchical 4       

internal_complete <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "complete", validation="internal")
summary(internal_complete)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                                 3       4
                                          
hierarchical Connectivity  11.6302 14.5591
             Dunn           0.1492  0.2252
             Silhouette     0.4476  0.4423

Optimal Scores:

             Score   Method       Clusters
Connectivity 11.6302 hierarchical 3       
Dunn          0.2252 hierarchical 4       
Silhouette    0.4476 hierarchical 3       

single linkage

stability_single <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "single", validation="stability")
summary(stability_single)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                       3      4
                               
hierarchical APN  0.0329 0.1514
             AD   1.4251 1.3090
             ADM  0.2596 0.4074
             FOM  0.8319 0.7066

Optimal Scores:

    Score  Method       Clusters
APN 0.0329 hierarchical 3       
AD  1.3090 hierarchical 4       
ADM 0.2596 hierarchical 3       
FOM 0.7066 hierarchical 4       

internal_single <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "single", validation="internal")
summary(internal_single)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                                 3       4
                                          
hierarchical Connectivity   7.9187 10.4187
             Dunn           0.4197  0.2379
             Silhouette     0.4380  0.3596

Optimal Scores:

             Score  Method       Clusters
Connectivity 7.9187 hierarchical 3       
Dunn         0.4197 hierarchical 3       
Silhouette   0.4380 hierarchical 3       

average linkage

stability_average <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "average", validation="stability")
summary(stability_average)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                       3      4
                               
hierarchical APN  0.1267 0.2082
             AD   1.3380 1.0182
             ADM  0.3899 0.4581
             FOM  0.6967 0.5956

Optimal Scores:

    Score  Method       Clusters
APN 0.1267 hierarchical 3       
AD  1.0182 hierarchical 4       
ADM 0.3899 hierarchical 3       
FOM 0.5956 hierarchical 4       

internal_average <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "average", validation="internal")
summary(internal_average)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                                 3       4
                                          
hierarchical Connectivity   7.9187 14.5591
             Dunn           0.4197  0.2252
             Silhouette     0.4380  0.4423

Optimal Scores:

             Score  Method       Clusters
Connectivity 7.9187 hierarchical 3       
Dunn         0.4197 hierarchical 3       
Silhouette   0.4423 hierarchical 4       

ward linkage

stability_ward <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "ward", validation="stability")

The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"
The "ward" method has been renamed to "ward.D"; note new "ward.D2"

summary(stability_ward)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                       3      4
                               
hierarchical APN  0.2358 0.2394
             AD   1.2341 1.0223
             ADM  0.5263 0.4766
             FOM  0.6768 0.5900

Optimal Scores:

    Score  Method       Clusters
APN 0.2358 hierarchical 3       
AD  1.0223 hierarchical 4       
ADM 0.4766 hierarchical 4       
FOM 0.5900 hierarchical 4       

internal_ward <-clValid(scaled_jabar, 3:4, clMethods = c("hierarchical"), method = "ward", validation="internal")

The "ward" method has been renamed to "ward.D"; note new "ward.D2"

summary(internal_ward)

Clustering Methods:
 hierarchical 

Cluster sizes:
 3 4 

Validation Measures:
                                 3       4
                                          
hierarchical Connectivity  17.8556 20.9083
             Dunn           0.1632  0.1632
             Silhouette     0.3217  0.3285

Optimal Scores:

             Score   Method       Clusters
Connectivity 17.8556 hierarchical 3       
Dunn          0.1632 hierarchical 3       
Silhouette    0.3285 hierarchical 4       

# Membuat data frame untuk setiap metode linkage
complete <- data.frame(
  Score = c(0.1821, 1.0182, 0.4581, 0.5956),
  Method = rep("complete", 4),
  Clusters = c(3, 4, 4, 4),
  Metric = c("APN", "AD", "ADM", "FOM")
)

single <- data.frame(
  Score = c(0.0329, 1.3090, 0.2596, 0.7066),
  Method = rep("single", 4),
  Clusters = c(3, 4, 3, 4),
  Metric = c("APN", "AD", "ADM", "FOM")
)

average <- data.frame(
  Score = c(0.1267, 1.0182, 0.3899, 0.5956),
  Method = rep("average", 4),
  Clusters = c(3, 4, 3, 4),
  Metric = c("APN", "AD", "ADM", "FOM")
)

ward <- data.frame(
  Score = c(0.2358, 1.0223, 0.4766, 0.5900),
  Method = rep("ward", 4),
  Clusters = c(3, 4, 4, 4),
  Metric = c("APN", "AD", "ADM", "FOM")
)

# Menggabungkan semua data frame menjadi satu
all_data <- rbind(complete, single, average, ward)

# Mengidentifikasi baris dengan skor terkecil untuk setiap metrik
best_scores <- all_data %>% 
  group_by(Metric) %>% 
  filter(Score == min(Score)) %>% 
  distinct()

# Menyusun ulang kolom sesuai dengan permintaan
best_scores <- best_scores %>% 
  dplyr::select(Metric, Score, Method, Clusters)

print(best_scores)

# A tibble: 5 × 4
# Groups:   Metric [4]
  Metric  Score Method   Clusters
  <chr>   <dbl> <chr>       <dbl>
1 AD     1.02   complete        4
2 APN    0.0329 single          3
3 ADM    0.260  single          3
4 AD     1.02   average         4
5 FOM    0.59   ward            4

Gerombol non-hirarki

library(factoextra)
# Plot cluster results
p1 <- fviz_nbclust(scaled_jabar, FUN = kmeans, method = "wss", 
                   k.max = 10) +
  ggtitle("(A) Elbow method")
p2 <- fviz_nbclust(scaled_jabar, FUN = kmeans, method = "silhouette", 
                   k.max = 10) +
  ggtitle("(B) Silhouette method")
#p3 <- fviz_nbclust(scaled_jabar, FUN = hcut, method = "gap_stat", 
                   #k.max = 16) +
  #ggtitle("(C) Gap statistic")

# Display plots side by side
gridExtra::grid.arrange(p1, p2, nrow = 1)

set.seed(123)
k3_means<-kmeans(scaled_jabar, centers = 3, iter.max = 100, nstart = 25)
k3_means

K-means clustering with 3 clusters of sizes 6, 18, 3

Cluster means:
         AHH        RLS        PPK
1  0.6136735  0.9911146  0.2727383
2 -0.5083770 -0.6237384 -0.4760564
3  1.8229152  1.7602012  2.3108616

Clustering vector:
         BANDUNG    BANDUNG BARAT           BEKASI            BOGOR 
               2                2                1                2 
          CIAMIS          CIANJUR          CIREBON            GARUT 
               2                2                2                2 
       INDRAMAYU         KARAWANG     KOTA BANDUNG      KOTA BANJAR 
               2                2                3                2 
     KOTA BEKASI       KOTA BOGOR      KOTA CIMAHI     KOTA CIREBON 
               3                1                1                1 
      KOTA DEPOK    KOTA SUKABUMI KOTA TASIKMALAYA         KUNINGAN 
               3                1                1                2 
      MAJALENGKA      PANGANDARAN       PURWAKARTA           SUBANG 
               2                2                2                2 
        SUKABUMI         SUMEDANG      TASIKMALAYA 
               2                2                2 

Within cluster sum of squares by cluster:
[1]  3.0309365 14.3566051  0.9942159
 (between_SS / total_SS =  76.4 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

head(scaled_jabar)

                      AHH        RLS        PPK
BANDUNG       -0.02954018  0.1579374 -0.2144463
BANDUNG BARAT -0.38973985 -0.4394781 -0.9010936
BEKASI         0.40784512  0.4806791  0.2521867
BOGOR         -0.67275387 -0.3433422 -0.1574368
CIAMIS         0.07337401 -0.5356139 -0.7499129
CIANJUR       -0.82712515 -1.1330294 -1.2245695

set.seed(123)
k4_means<-kmeans(scaled_jabar, centers = 4, iter.max = 100, nstart = 25)
k4_means

K-means clustering with 4 clusters of sizes 1, 6, 17, 3

Cluster means:
         AHH        RLS        PPK
1 -2.7310377 -0.6248829 -1.2515962
2  0.6136735  0.9911146  0.2727383
3 -0.3776323 -0.6236711 -0.4304364
4  1.8229152  1.7602012  2.3108616

Clustering vector:
         BANDUNG    BANDUNG BARAT           BEKASI            BOGOR 
               3                3                2                3 
          CIAMIS          CIANJUR          CIREBON            GARUT 
               3                3                3                3 
       INDRAMAYU         KARAWANG     KOTA BANDUNG      KOTA BANJAR 
               3                3                4                3 
     KOTA BEKASI       KOTA BOGOR      KOTA CIMAHI     KOTA CIREBON 
               4                2                2                2 
      KOTA DEPOK    KOTA SUKABUMI KOTA TASIKMALAYA         KUNINGAN 
               4                2                2                3 
      MAJALENGKA      PANGANDARAN       PURWAKARTA           SUBANG 
               3                3                3                3 
        SUKABUMI         SUMEDANG      TASIKMALAYA 
               3                3                1 

Within cluster sum of squares by cluster:
[1] 0.0000000 3.0309365 8.4889402 0.9942159
 (between_SS / total_SS =  84.0 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

k4_means$cluster

         BANDUNG    BANDUNG BARAT           BEKASI            BOGOR 
               3                3                2                3 
          CIAMIS          CIANJUR          CIREBON            GARUT 
               3                3                3                3 
       INDRAMAYU         KARAWANG     KOTA BANDUNG      KOTA BANJAR 
               3                3                4                3 
     KOTA BEKASI       KOTA BOGOR      KOTA CIMAHI     KOTA CIREBON 
               4                2                2                2 
      KOTA DEPOK    KOTA SUKABUMI KOTA TASIKMALAYA         KUNINGAN 
               4                2                2                3 
      MAJALENGKA      PANGANDARAN       PURWAKARTA           SUBANG 
               3                3                3                3 
        SUKABUMI         SUMEDANG      TASIKMALAYA 
               3                3                1 

library(cluster)
library(fpc)
library(clusterSim)

# Function to calculate S_W and S_B
calculate_sw_sb_ratio <- function(data, kmeans_result) {
  # Total Sum of Squares
  tss <- sum((data - colMeans(data))^2)
  
  # Within-cluster Sum of Squares (S_W)
  sw <- sum(kmeans_result$withinss)
  
  # Between-cluster Sum of Squares (S_B)
  sb <- tss - sw
  
  # Ratio S_W / S_B
  ratio <- sw / sb
  
  return(list(S_W = sw, S_B = sb, ratio = ratio))
}

# Initialize an empty dataframe to store results
internal_metrics <- data.frame(
  K = integer(),
  Dunn_Index = numeric(),
  Silhouette = numeric(),
  Calinski_Harabasz = numeric(),
  Davies_Bouldin = numeric(),
  sw_sb = numeric()
)

# List of k-means results for k = 3 and k = 4
kmeans_list <- list(k3_means, k4_means)

# Loop through the list of k-means results
for (k in 3:4) {
  k_means <- kmeans_list[[k - 2]]
  
  # Calculate internal validation metrics
  distance_matrix <- dist(scaled_jabar)
  
  # Ensure that the length of the clustering vector matches the number of rows in the distance matrix
  if(length(k_means$cluster) != nrow(as.matrix(distance_matrix))) {
    stop("The length of the clustering vector does not match the number of rows in the distance matrix for k = ", k)
  }
  
  kmeans_internal <- cluster.stats(distance_matrix, k_means$cluster, silhouette = TRUE)
  kmeans_db <- index.DB(scaled_jabar, k_means$cluster, distance_matrix, centrotypes = "centroids", p = 2, q = 2)
  sw_sb_k <- calculate_sw_sb_ratio(scaled_jabar, k_means)
  
  # Extract metrics
  dunn_index <- round(kmeans_internal$dunn, 3)
  silhouette_score <- round(kmeans_internal$avg.silwidth, 3)
  calinski_harabasz_index <- round(kmeans_internal$ch, 3)
  db_index <- round(kmeans_db$DB, 3)
  swsb_ratio <- round(sw_sb_k$ratio, 3)
  
  # Add to the dataframe
  internal_metrics <- rbind(internal_metrics, data.frame(
    K = k,
    Dunn_Index = dunn_index,
    Silhouette = silhouette_score,
    Calinski_Harabasz = calinski_harabasz_index,
    Davies_Bouldin = db_index,
    sw_sb = swsb_ratio
  ))
}

# Print the dataframe
print(internal_metrics)

  K Dunn_Index Silhouette Calinski_Harabasz Davies_Bouldin sw_sb
1 3      0.158      0.477            38.920          0.680 0.308
2 4      0.276      0.469            40.119          0.551 0.191

internal_metrics_kmeans <- internal_metrics

# Determine the optimal k based on different criteria
optimal_k_silhouette <- internal_metrics_kmeans$K[which.max(internal_metrics_kmeans$Silhouette)]
optimal_k_dunn <- internal_metrics_kmeans$K[which.max(internal_metrics_kmeans$Dunn_Index)]
optimal_k_ch <- internal_metrics_kmeans$K[which.max(internal_metrics_kmeans$Calinski_Harabasz)]
optimal_k_db <- internal_metrics_kmeans$K[which.min(internal_metrics_kmeans$Davies_Bouldin)]
optimal_k_sw_sb <- internal_metrics_kmeans$K[which.min(internal_metrics_kmeans$sw_sb)]

# Print the optimal k for each metric
cat("Optimal k based on Silhouette Score:", optimal_k_silhouette, "\n")

Optimal k based on Silhouette Score: 3 

cat("Optimal k based on Dunn Index:", optimal_k_dunn, "\n")

Optimal k based on Dunn Index: 4 

cat("Optimal k based on Calinski-Harabasz Index:", optimal_k_ch, "\n")

Optimal k based on Calinski-Harabasz Index: 4 

cat("Optimal k based on Davies-Bouldin Index:", optimal_k_db, "\n")

Optimal k based on Davies-Bouldin Index: 4 

cat("Optimal k based on sw/sb Ratio:", optimal_k_sw_sb, "\n")

Optimal k based on sw/sb Ratio: 4 

stab_mean <-clValid(scaled_jabar, 3:4, clMethods = c("kmeans"), validation="stability")
summary(stab_mean)

Clustering Methods:
 kmeans 

Cluster sizes:
 3 4 

Validation Measures:
                 3      4
                         
kmeans APN  0.1835 0.2489
       AD   1.2835 1.0043
       ADM  0.5207 0.4374
       FOM  0.6597 0.5857

Optimal Scores:

    Score  Method Clusters
APN 0.1835 kmeans 3       
AD  1.0043 kmeans 4       
ADM 0.4374 kmeans 4       
FOM 0.5857 kmeans 4       

Best algorithm

# Menentukan metode terbaik untuk setiap metrik
best_dunn_index <- score_data$method[which.max(score_data$Dunn_Index)]
best_silhouette <- score_data$method[which.max(score_data$Silhouette)]
best_calinski_harabasz <- score_data$method[which.max(score_data$Calinski_Harabasz)]
best_davies_bouldin <- score_data$method[which.min(score_data$Davies_Bouldin)]
best_sw_sb <- score_data$method[which.min(score_data$sw_sb)]

# Membuat data frame hasil terbaik
best_algorithms <- data.frame(
  Metric = c("Dunn_Index", "Silhouette", "Calinski_Harabasz", "Davies_Bouldin", "sw_sb"),
  Best_Method = c(best_dunn_index, best_silhouette, best_calinski_harabasz, best_davies_bouldin, best_sw_sb)
)

# Menampilkan data frame hasil terbaik
print(best_algorithms)

             Metric               Best_Method
1        Dunn_Index kmeans (non-hierarchical)
2        Silhouette kmeans (non-hierarchical)
3 Calinski_Harabasz kmeans (non-hierarchical)
4    Davies_Bouldin kmeans (non-hierarchical)
5             sw_sb kmeans (non-hierarchical)

Plot kota

library(dplyr)
merged_jabar <- gdf %>%
  right_join(df_jabar2, by = c("KAB_KOTA" = "KOTA"))
head(merged_jabar,10)

Simple feature collection with 10 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 106.4011 ymin: -7.74015 xmax: 108.8338 ymax: -5.91377
Geodetic CRS:  WGS 84
        KAB_KOTA   Provinsi   AHH  RLS   PPK   IPM rentang_IPM g4_kmeans
1        BANDUNG Jawa Barat 74.92 9.10 11018 73.74      Tinggi         3
2  BANDUNG BARAT Jawa Barat 74.78 8.23  9392 69.61      Sedang         3
3         BEKASI Jawa Barat 75.09 9.57 12123 75.76      Tinggi         2
4          BOGOR Jawa Barat 74.67 8.37 11153 71.78      Tinggi         3
5         CIAMIS Jawa Barat 74.96 8.09  9750 72.05      Tinggi         3
6        CIANJUR Jawa Barat 74.61 7.22  8626 66.55      Sedang         3
7        CIREBON Jawa Barat 74.71 7.64 11128 70.95      Tinggi         3
8          GARUT Jawa Barat 74.66 7.84  8685 68.11      Sedang         3
9      INDRAMAYU Jawa Barat 74.61 6.94 10580 69.25      Sedang         3
10      KARAWANG Jawa Barat 74.90 8.04 12392 72.35      Tinggi         3
                         geometry
1  MULTIPOLYGON (((107.7327 -6...
2  MULTIPOLYGON (((107.4393 -6...
3  MULTIPOLYGON (((107.0314 -5...
4  MULTIPOLYGON (((106.9709 -6...
5  MULTIPOLYGON (((108.5653 -7...
6  MULTIPOLYGON (((107.23 -6.6...
7  MULTIPOLYGON (((108.6733 -6...
8  MULTIPOLYGON (((107.918 -6....
9  MULTIPOLYGON (((108.3674 -6...
10 MULTIPOLYGON (((107.1123 -5...

boxplot

# Membuat boxplot untuk variabel AHH
p1 <- ggplot(sorted_jabar, aes(x = g4_kmeans, y = AHH, fill = g4_kmeans)) +
  geom_boxplot() +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "",
       x = "",
       y = "AHH") +
  theme(legend.position = "none")

# Membuat boxplot untuk variabel RLS
p2 <- ggplot(sorted_jabar, aes(x = g4_kmeans, y = RLS, fill = g4_kmeans)) +
  geom_boxplot() +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "",
       x = "",
       y = "RLS") +
  theme(legend.position = "none")

# Membuat boxplot untuk variabel PPK
p3 <- ggplot(sorted_jabar, aes(x = g4_kmeans, y = PPK, fill = g4_kmeans)) +
  geom_boxplot() +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "",
       x = "",
       y = "PPK") +
  theme(legend.position = "none")

# Membuat boxplot untuk variabel IPM
p4 <- ggplot(sorted_jabar, aes(x = g4_kmeans, y = IPM, fill = g4_kmeans)) +
  geom_boxplot() +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "",
       x = "Cluster",
       y = "IPM") +
  theme(legend.position = "none")
grid.arrange(p1, p2, p3, nrow = 2)

barplot

# Menghitung rataan IPM per cluster
df_ipm_means <- sorted_jabar %>%
  group_by(g4_kmeans) %>%
  summarise(IPM = mean(IPM, na.rm = TRUE))

# Menampilkan dataframe hasil
print(df_ipm_means)

# A tibble: 4 × 2
  g4_kmeans   IPM
  <chr>     <dbl>
1 G1         82.8
2 G2         76.7
3 G3         70.7
4 G4         67.8

# Membuat plot bar untuk variabel IPM
ggplot(df_ipm_means, aes(x =g4_kmeans, y = IPM, fill = g4_kmeans)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "Rataan IPM Berdasarkan Cluster",
       x = "Cluster",
       y = "Rataan IPM") +
  theme(legend.position = "none") +
  geom_text(aes(label = round(IPM, 2)), vjust = -0.5)

# Memuat perpustakaan yang diperlukan
library(dplyr)
library(tidyr)
library(ggplot2)

# Menghitung rataan tiap indikator per cluster
df_means <- sorted_jabar %>%
  group_by(g4_kmeans) %>%
  summarise(
    AHH = mean(AHH, na.rm = TRUE),
    RLS = mean(RLS, na.rm = TRUE),
    PPK = mean(PPK, na.rm = TRUE)
    #IPM = mean(IPM, na.rm = TRUE)
  )

# Mengubah dataframe dari format lebar ke panjang untuk ggplot
df_means_long <- df_means %>%
  pivot_longer(cols = AHH:PPK, names_to = "Indikator", values_to = "Rata-rata")

# Membuat plot bar untuk rataan tiap indikator
ggplot(df_means_long, aes(x = g4_kmeans, y = `Rata-rata`, fill = g4_kmeans)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ Indikator, scales = "free_y") +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c")) +
  theme_minimal() +
  labs(title = "",
       x = "Cluster",
       y = "Rata-rata") +
  theme(legend.position = "none") +
  geom_text(aes(label = round(`Rata-rata`, 2)), position = position_dodge(width = 0.9), vjust = -0.5)

# Memuat perpustakaan yang diperlukan
library(sf)
library(ggplot2)
library(gridExtra)

# Membuat plot untuk rentang_ipm
p1 <- ggplot(merged_jabar2) +
  geom_sf(aes(fill = rentang_IPM)) +
  scale_fill_manual(values = c("Sangat Tinggi" = "#bdd7e7", "Tinggi" = "#6baed6", "Sedang" = "#3182bd", "Rendah" = "#08519c"),
                    name = "Rentang IPM") +
  theme_minimal() +
  labs(title = "Sebaran Kota Berdasarkan status IPM",
       x = "Longitude",
       y = "Latitude") +
  theme(panel.background = element_blank())

# Membuat plot untuk cluster
p2 <- ggplot(merged_jabar2) +
  geom_sf(aes(fill = g4_kmeans)) +
  scale_fill_manual(values = c("G1" = "#bdd7e7", "G2" = "#6baed6", "G3" = "#3182bd", "G4" = "#08519c"),
                    name = "Gerombol") +
  theme_minimal() +
  labs(title = "Sebaran Kota Berdasarkan k-means",
       x = "Longitude",
       y = "Latitude") +
  theme(panel.background = element_blank())


# Menampilkan kedua plot dalam satu grid
grid.arrange(p1, p2, ncol = 1)