1 Library Data

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)
library(cluster)      # pam() for k-medoids
library(factoextra)   # visualization

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

library(readr)
library(readxl)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.4     ✔ tibble    3.3.0
## ✔ purrr     1.2.0

## ── 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

2 Input Data

df <- read_excel("C:/Users/user/Downloads/0_df_long.xlsx")
df

## # A tibble: 340,620 × 9
##    KabKot          Tahun Bulan   Produk Harga Kategori KodeBPS KodeProv NamaProv
##    <chr>           <dbl> <chr>   <chr>  <chr> <chr>      <dbl>    <dbl> <chr>   
##  1 Kab. Aceh Barat  2022 Januari Beras… 11429 Beras       11.0       11 Aceh    
##  2 Kab. Aceh Barat  2022 Januari Beras… 9979  Beras       11.0       11 Aceh    
##  3 Kab. Aceh Barat  2022 Januari Bawan… 31000 Bawang      11.0       11 Aceh    
##  4 Kab. Aceh Barat  2022 Januari Bawan… 28636 Bawang      11.0       11 Aceh    
##  5 Kab. Aceh Barat  2022 Januari Cabai… 19409 Cabai       11.0       11 Aceh    
##  6 Kab. Aceh Barat  2022 Januari Cabai… 45706 Cabai       11.0       11 Aceh    
##  7 Kab. Aceh Barat  2022 Januari Dagin… 1438… Protein     11.0       11 Aceh    
##  8 Kab. Aceh Barat  2022 Januari Dagin… 32955 Protein     11.0       11 Aceh    
##  9 Kab. Aceh Barat  2022 Januari Telur… 26136 Protein     11.0       11 Aceh    
## 10 Kab. Aceh Barat  2022 Januari Gula … 14545 Gula        11.0       11 Aceh    
## # ℹ 340,610 more rows

3 Eksplorasi Data

summary(df)

##     KabKot              Tahun         Bulan              Produk         
##  Length:340620      Min.   :2022   Length:340620      Length:340620     
##  Class :character   1st Qu.:2022   Class :character   Class :character  
##  Mode  :character   Median :2023   Mode  :character   Mode  :character  
##                     Mean   :2023                                        
##                     3rd Qu.:2024                                        
##                     Max.   :2025                                        
##     Harga             Kategori            KodeBPS         KodeProv   
##  Length:340620      Length:340620      Min.   :11.01   Min.   :11.0  
##  Class :character   Class :character   1st Qu.:17.09   1st Qu.:17.0  
##  Mode  :character   Mode  :character   Median :35.21   Median :35.0  
##                                        Mean   :43.81   Mean   :43.6  
##                                        3rd Qu.:71.01   3rd Qu.:71.0  
##                                        Max.   :92.71   Max.   :92.0  
##    NamaProv        
##  Length:340620     
##  Class :character  
##  Mode  :character  
##                    
##                    
##

glimpse(df)

## Rows: 340,620
## Columns: 9
## $ KabKot   <chr> "Kab. Aceh Barat", "Kab. Aceh Barat", "Kab. Aceh Barat", "Kab…
## $ Tahun    <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2…
## $ Bulan    <chr> "Januari", "Januari", "Januari", "Januari", "Januari", "Janua…
## $ Produk   <chr> "Beras Premium", "Beras Medium", "Bawang Merah", "Bawang Puti…
## $ Harga    <chr> "11429", "9979", "31000", "28636", "19409", "45706", "143864"…
## $ Kategori <chr> "Beras", "Beras", "Bawang", "Bawang", "Cabai", "Cabai", "Prot…
## $ KodeBPS  <dbl> 11.05, 11.05, 11.05, 11.05, 11.05, 11.05, 11.05, 11.05, 11.05…
## $ KodeProv <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 1…
## $ NamaProv <chr> "Aceh", "Aceh", "Aceh", "Aceh", "Aceh", "Aceh", "Aceh", "Aceh…

3.1 Mengubah Tipe Data pada Harga Menjadi Numerik

df <- df %>%
  mutate(Harga = as.character(Harga),
         Harga = na_if(Harga, "NA"),
         Harga = as.numeric(Harga))

3.2 Mengurutkan Bulan dari Januari ke Desember

df$Bulan <- factor(
  df$Bulan,
  levels = c("Januari", "Februari", "Maret", "April", "Mei", "Juni",
             "Juli", "Agustus", "September", "Oktober", "November", "Desember")
)

3.3 Cek Missing Value pada Kolom Harga secara keseluruhan

sum(is.na(df$Harga))  # ini NA asli

## [1] 61513

4 Filter Data Beras di Jawa Tengah Tahun 2023

df_jateng_beras_2023 <- df %>%
  filter(
    NamaProv == "Jawa Tengah",
    Tahun == 2023,
    Kategori == "Beras"
  )
df_jateng_beras_2023

## # A tibble: 840 × 9
##    KabKot            Tahun Bulan Produk Harga Kategori KodeBPS KodeProv NamaProv
##    <chr>             <dbl> <fct> <chr>  <dbl> <chr>      <dbl>    <dbl> <chr>   
##  1 Kab. Banjarnegara  2023 Janu… Beras… 12000 Beras       33.0       33 Jawa Te…
##  2 Kab. Banjarnegara  2023 Janu… Beras… 11067 Beras       33.0       33 Jawa Te…
##  3 Kab. Banjarnegara  2023 Febr… Beras… 12667 Beras       33.0       33 Jawa Te…
##  4 Kab. Banjarnegara  2023 Febr… Beras… 11870 Beras       33.0       33 Jawa Te…
##  5 Kab. Banjarnegara  2023 Maret Beras… 12617 Beras       33.0       33 Jawa Te…
##  6 Kab. Banjarnegara  2023 Maret Beras… 11833 Beras       33.0       33 Jawa Te…
##  7 Kab. Banjarnegara  2023 April Beras… 12926 Beras       33.0       33 Jawa Te…
##  8 Kab. Banjarnegara  2023 April Beras… 12130 Beras       33.0       33 Jawa Te…
##  9 Kab. Banjarnegara  2023 Mei   Beras… 12000 Beras       33.0       33 Jawa Te…
## 10 Kab. Banjarnegara  2023 Mei   Beras… 11000 Beras       33.0       33 Jawa Te…
## # ℹ 830 more rows

4.1 Eksplorasi Data Beras Jawa Tengah 2023

4.1.1 Cek Missing Value pada Data Beras Jawa Tengah 2023

sum(is.na(df_jateng_beras_2023$Harga))

## [1] 0

Tidak terdapat missing value pada kategori beras di Jawa Tengah Tahun 2023.

4.1.2 Statistik Deskriptif Harga Beras Jawa Tengah 2023

stats_beras <- df_jateng_beras_2023 %>%
  group_by(Produk) %>%
  summarise(
    mean_price = mean(Harga, na.rm = TRUE),
    median_price = median(Harga, na.rm = TRUE),
    sd_price = sd(Harga, na.rm = TRUE),
    min_price = min(Harga, na.rm = TRUE),
    max_price = max(Harga, na.rm = TRUE),
    n = n()
  )

stats_beras

## # A tibble: 2 × 7
##   Produk        mean_price median_price sd_price min_price max_price     n
##   <chr>              <dbl>        <dbl>    <dbl>     <dbl>     <dbl> <int>
## 1 Beras Medium      11907.       11728.     943.      9917     14033   420
## 2 Beras Premium     13239.       13018.     854.     11067     15033   420

4.1.3 Boxplot Perbandingan Harga Premium vs Medium

ggplot(df_jateng_beras_2023, aes(x = Produk, y = Harga, fill = Produk)) +
  geom_boxplot(alpha = 0.6) +
  scale_fill_manual(values = c("#FFD27F", "#FF8C00")) +  # warna orange
  theme_minimal(base_size = 13) +
  labs(
    title = "Perbandingan Harga Beras Premium vs Medium di Jawa Tengah (2023)",
    x = "Jenis Beras",
    y = "Harga (Rp)"
  )

5 Filter Data Beras Medium Jawa Tengah 2023

beras__medium_jateng_2023 <- df %>%
  filter(
    NamaProv == "Jawa Tengah",
    Tahun == 2023,
    Produk == "Beras Medium"
  )
beras__medium_jateng_2023

## # A tibble: 420 × 9
##    KabKot            Tahun Bulan Produk Harga Kategori KodeBPS KodeProv NamaProv
##    <chr>             <dbl> <fct> <chr>  <dbl> <chr>      <dbl>    <dbl> <chr>   
##  1 Kab. Banjarnegara  2023 Janu… Beras… 11067 Beras       33.0       33 Jawa Te…
##  2 Kab. Banjarnegara  2023 Febr… Beras… 11870 Beras       33.0       33 Jawa Te…
##  3 Kab. Banjarnegara  2023 Maret Beras… 11833 Beras       33.0       33 Jawa Te…
##  4 Kab. Banjarnegara  2023 April Beras… 12130 Beras       33.0       33 Jawa Te…
##  5 Kab. Banjarnegara  2023 Mei   Beras… 11000 Beras       33.0       33 Jawa Te…
##  6 Kab. Banjarnegara  2023 Juni  Beras… 11667 Beras       33.0       33 Jawa Te…
##  7 Kab. Banjarnegara  2023 Juli  Beras… 12000 Beras       33.0       33 Jawa Te…
##  8 Kab. Banjarnegara  2023 Agus… Beras… 12000 Beras       33.0       33 Jawa Te…
##  9 Kab. Banjarnegara  2023 Sept… Beras… 12385 Beras       33.0       33 Jawa Te…
## 10 Kab. Banjarnegara  2023 Okto… Beras… 13482 Beras       33.0       33 Jawa Te…
## # ℹ 410 more rows

6 Eksplorasi Data Beras Medium Jawa Tengah 2023

glimpse(beras__medium_jateng_2023)

## Rows: 420
## Columns: 9
## $ KabKot   <chr> "Kab. Banjarnegara", "Kab. Banjarnegara", "Kab. Banjarnegara"…
## $ Tahun    <dbl> 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2…
## $ Bulan    <fct> Januari, Februari, Maret, April, Mei, Juni, Juli, Agustus, Se…
## $ Produk   <chr> "Beras Medium", "Beras Medium", "Beras Medium", "Beras Medium…
## $ Harga    <dbl> 11067, 11870, 11833, 12130, 11000, 11667, 12000, 12000, 12385…
## $ Kategori <chr> "Beras", "Beras", "Beras", "Beras", "Beras", "Beras", "Beras"…
## $ KodeBPS  <dbl> 33.04, 33.04, 33.04, 33.04, 33.04, 33.04, 33.04, 33.04, 33.04…
## $ KodeProv <dbl> 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, 3…
## $ NamaProv <chr> "Jawa Tengah", "Jawa Tengah", "Jawa Tengah", "Jawa Tengah", "…

6.1 Statistik Deskriptif per Kabupaten Kota

stat_kabkot <- beras__medium_jateng_2023 %>%
  group_by(KabKot) %>%
  summarise(
    mean_price = mean(Harga, na.rm = TRUE),
    sd_price   = sd(Harga, na.rm = TRUE),
    min_price  = min(Harga, na.rm = TRUE),
    max_price  = max(Harga, na.rm = TRUE),
    n_obs      = n(),
    na_obs     = sum(is.na(Harga))
  )

stat_kabkot

## # A tibble: 35 × 7
##    KabKot            mean_price sd_price min_price max_price n_obs na_obs
##    <chr>                  <dbl>    <dbl>     <dbl>     <dbl> <int>  <int>
##  1 Kab. Banjarnegara     12180      839.     11000     13482    12      0
##  2 Kab. Banyumas         11733.     512.     11161     12545    12      0
##  3 Kab. Batang           12074.     689.     11007     13159    12      0
##  4 Kab. Blora            11187.     771.     10480     12508    12      0
##  5 Kab. Boyolali         11637.    1077.     10297     13452    12      0
##  6 Kab. Brebes           12284.     837.     11500     13821    12      0
##  7 Kab. Cilacap          11384.     908.      9917     12933    12      0
##  8 Kab. Demak            12230.     810.     11370     13517    12      0
##  9 Kab. Grobogan         11372.    1159.      9964     13000    12      0
## 10 Kab. Jepara           11894.     799.     11000     13162    12      0
## # ℹ 25 more rows

6.2 Grafik Tren Harga per Bulan

beras__medium_jateng_2023 %>%
  group_by(Bulan) %>%
  summarise(mean_harga = mean(Harga, na.rm = TRUE)) %>%
  ggplot(aes(x = Bulan, y = mean_harga, group = 1)) +
  geom_line(size = 1.2, color = "orange") +
  geom_point(size = 3, color = "darkred") +
  labs(title = "Tren Harga Rata-Rata Beras Medium di Jawa Tengah Tahun 2023",
       x = "Bulan",
       y = "Harga") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

7 Tahapan Analisis Cluster

7.1 Data Scaling

# Scaling variabel mean
scaled_data <- stat_kabkot %>%
  select(mean_price) %>%
  scale()
scaled_data <- as.data.frame(scaled_data)
scaled_data <- mutate_all(scaled_data, as.numeric)
# Menambahkan nama kabupaten sebagai rowname
rownames(scaled_data) <- stat_kabkot$KabKot
scaled_data

##                    mean_price
## Kab. Banjarnegara  0.67680712
## Kab. Banyumas     -0.43213025
## Kab. Batang        0.41506650
## Kab. Blora        -1.78462893
## Kab. Boyolali     -0.66825458
## Kab. Brebes        0.93338316
## Kab. Cilacap      -1.29729971
## Kab. Demak         0.79993056
## Kab. Grobogan     -1.32518842
## Kab. Jepara       -0.03135929
## Kab. Karanganyar   0.04073816
## Kab. Kebumen      -1.76025214
## Kab. Kendal        0.66833721
## Kab. Klaten        0.10271305
## Kab. Kudus         0.36982483
## Kab. Magelang     -1.42558774
## Kab. Pati          0.16696035
## Kab. Pekalongan   -0.34598515
## Kab. Pemalang      0.79084091
## Kab. Purbalingga  -0.06461915
## Kab. Purworejo    -1.75488099
## Kab. Rembang      -0.13712977
## Kab. Semarang      1.13294231
## Kab. Sragen       -0.97978169
## Kab. Sukoharjo    -0.22781970
## Kab. Tegal         0.02007986
## Kab. Temanggung    0.18142116
## Kab. Wonogiri      0.72638703
## Kab. Wonosobo     -0.40692713
## Kota Magelang      0.87285435
## Kota Pekalongan    1.89048205
## Kota Salatiga      1.26102375
## Kota Semarang      1.57688910
## Kota Surakarta    -1.32725425
## Kota Tegal         1.34241743

7.2 Pemilihan K (Elbow Method)

fviz_nbclust(scaled_data, kmeans, method = "wss") + ggtitle("Elbow method")

Elbow Method menunjukkan bahwa penurunan WSS mulai melandai pada k = 3, sehingga jumlah cluster optimal yang digunakan adalah 3 cluster.

7.3 Menentukan K pada K-Means

# range k yang ingin diuji
k_range <- 2:5

# wadah hasil silhouette
sil_scores <- numeric(length(k_range))

# looping
for (i in seq_along(k_range)) {
  k <- k_range[i]
  set.seed(123)
  
  km <- kmeans(scaled_data, centers = k, nstart = 25)
  sil <- silhouette(km$cluster, dist(scaled_data))
  
  sil_scores[i] <- mean(sil[, 3])
}

# tampilkan hasil
data.frame(k = k_range, silhouette_score = sil_scores)

##   k silhouette_score
## 1 2        0.5657356
## 2 3        0.6331908
## 3 4        0.6116777
## 4 5        0.5715312

Di antara hasil k yang lain, k=3 menunjukkan score silhouette yang paling besar.

7.3.1 Menambahkan Hasil Cluster pada Data

set.seed(123)
km3 <- kmeans(scaled_data[, "mean_price"], centers = 3)
stat_kabkot$cluster <- as.factor(km3$cluster)

7.3.2 Visualisasi Scatter Plot K-Means Clustering k=3

# Buat data frame untuk visualisasi
viz_data <- data.frame(
  mean_price = scaled_data[, "mean_price"],   # ganti jika nama kolom berbeda
  cluster = factor(km3$cluster)               # pakai hasil k = 3
)

# Plot visualisasi
ggplot(viz_data, aes(x = mean_price, y = cluster, color = cluster)) +
  geom_point(size = 3, alpha = 0.8) +
  labs(
    title = "K-Means Clustering Berdasarkan Harga Beras 2023 (k = 3)",
    x = "Rata-rata Harga Beras (Scaled)",
    y = "Cluster"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
  )

Interpretasi:

Cluster 1 (merah) → daerah dengan mean price tertinggi (paling mahal).
Cluster 2 (hijau) → daerah dengan mean price kategori menengah.
Cluster 3 (biru) → daerah dengan mean price terendah (paling murah).

Pemisahan titik-titik cukup jelas hanya berdasarkan satu variabel mean price, sehingga:

Cluster 1 berada di sisi paling kanan (harga tinggi),
Cluster 2 di tengah,
Cluster 3 di sisi kiri (harga rendah).

Kesimpulannya, data dapat dipisah dengan baik menjadi 3 kelompok harga, yaitu tinggi, sedang, dan rendah.

7.3.3 Visualisasi Bar Plot K-Means Clustering k=3

ggplot(stat_kabkot, aes(x = reorder(KabKot, mean_price), 
                        y = mean_price, 
                        fill = cluster)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Bar Plot Mean Price per Kabupaten/Kota Berdasarkan Cluster (K=3)",
    x = "Kabupaten/Kota",
    y = "Mean Price (Scaled)"
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(size = 6),   # biar nama kab/kota muat
    legend.position = "bottom"
  )

Interpretasi:

Cluster 1 (harga tinggi) Kabupaten/kota dalam cluster ini memiliki rata-rata harga beras yang paling tinggi dibandingkan wilayah lain. Pada scatter plot terlihat titik-titik berada di sisi kanan (nilai scaled positif), dan pada bar plot cluster ini menempati posisi nilai mean price paling besar.
Cluster 2 (harga sedang) Cluster ini berisi wilayah dengan harga beras kategori sedang, yaitu harga tidak setinggi cluster 1 namun juga tidak serendah cluster 3. Pada bar plot, kelompok ini mendominasi bagian tengah diagram.
Cluster 3 (harga rendah) Wilayah dalam cluster ini memiliki harga beras relatif lebih rendah. Pada scatter plot terlihat nilai scaled berada di sekitar posisi negatif, dan pada bar plot posisi cluster 3 berada di bagian bawah dengan mean price paling rendah.

Secara keseluruhan, k=3 mampu memisahkan kabupaten/kota berdasarkan tingkat harga beras menjadi tiga kelompok yang jelas, yakni tinggi, sedang, dan rendah, sehingga hasil ini dapat digunakan untuk analisis perbandingan harga antar wilayah atau dasar kebijakan terkait stabilisasi harga pangan.

7.4 Pemilihan K pada K-Medoids (PAM)

fviz_nbclust(scaled_data, kmeans, method = "silhouette") + ggtitle("Silhouette")

Nilai Silhouette tertinggi pada k=3, sekitar 0.68. Ini menandakan bahwa 3 cluster adalah jumlah paling optimal karena memberikan pemisahan antar cluster yang paling baik.

7.5 Menentukan K pada K-Medoids (PAM)

# range k yang ingin diuji
k_range <- 2:5

# wadah hasil silhouette
sil_scores_pam <- numeric(length(k_range))

# looping
for (i in seq_along(k_range)) {
  k <- k_range[i]
  set.seed(123)
  
  # K-Medoids (PAM)
  pam_model <- pam(scaled_data, k = k)
  
  # hitung silhouette
  sil <- silhouette(pam_model$clustering, dist(scaled_data))
  
  sil_scores_pam[i] <- mean(sil[, 3])
}

# tampilkan hasil
data.frame(
  k = k_range,
  silhouette_score = sil_scores_pam
)

##   k silhouette_score
## 1 2        0.6035651
## 2 3        0.6331908
## 3 4        0.6116777
## 4 5        0.5799488

Di antara hasil k yang lain, k=3 menunjukkan score silhouette yang paling besar.

7.5.1 Visualisasi Scatter Plot K-Medoids (PAM) Clustering k=3

ggplot(viz_data, aes(x = mean_price, y = cluster, color = cluster)) +
  geom_point(size = 3, alpha = 0.8) +
  labs(
    title = "K-Medoids (PAM) Clustering Berdasarkan Mean Price (k = 3)",
    x = "Rata-rata Harga (Scaled)",
    y = "Cluster"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
  )

Interpretasi:

Cluster 1 (merah) → daerah dengan mean price tertinggi (paling mahal).
Cluster 2 (hijau) → daerah dengan mean price kategori menengah.
Cluster 3 (biru) → daerah dengan mean price terendah (paling murah).

Pemisahan titik-titik cukup jelas hanya berdasarkan satu variabel mean price, sehingga:

Cluster 1 berada di sisi paling kanan (harga tinggi),
Cluster 2 di tengah,
Cluster 3 di sisi kiri (harga rendah).

Kesimpulannya, data dapat dipisah dengan baik menjadi 3 kelompok harga, yaitu tinggi, sedang, dan rendah.

7.5.2 Visualisasi Bar Plot K-Medoids (PAM) Clustering k=3

ggplot(stat_kabkot, aes(x = reorder(KabKot, mean_price),
                        y = mean_price,
                        fill = cluster)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Bar Plot Mean Price per Kabupaten/Kota Berdasarkan Cluster PAM (K = 3)",
    x = "Kabupaten/Kota",
    y = "Mean Price (Scaled)"
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(size = 6),
    legend.position = "bottom"
  )

Interpretasi:

Cluster 1 (harga tinggi) Kabupaten/kota dalam cluster ini memiliki rata-rata harga beras yang paling tinggi dibandingkan wilayah lain. Pada scatter plot terlihat titik-titik berada di sisi kanan (nilai scaled positif), dan pada bar plot cluster ini menempati posisi nilai mean price paling besar.
Cluster 2 (harga sedang) Cluster ini berisi wilayah dengan harga beras kategori sedang, yaitu harga tidak setinggi cluster 1 namun juga tidak serendah cluster 3. Pada bar plot, kelompok ini mendominasi bagian tengah diagram.
Cluster 3 (harga rendah) Wilayah dalam cluster ini memiliki harga beras relatif lebih rendah. Pada scatter plot terlihat nilai scaled berada di sekitar posisi negatif, dan pada bar plot posisi cluster 3 berada di bagian bawah dengan mean price paling rendah.

8 Perbandingan K-Means dan K-Medoids

k_range <- 2:5

results <- data.frame(
  k = integer(),
  method = character(),
  silhouette = numeric()
)

for (k in k_range) {
  
  # --- K-MEANS ---
  set.seed(123)
  km <- kmeans(scaled_data, centers = k, nstart = 25)
  sil_km <- silhouette(km$cluster, dist(scaled_data))
  results <- rbind(
    results,
    data.frame(k = k, method = "K-Means", silhouette = mean(sil_km[, 3]))
  )
  
  # --- K-MEDOIDS (PAM) ---
  set.seed(123)
  pam_model <- pam(scaled_data, k = k)
  sil_pam <- silhouette(pam_model$clustering, dist(scaled_data))
  results <- rbind(
    results,
    data.frame(k = k, method = "K-Medoids", silhouette = mean(sil_pam[, 3]))
  )
}

results

##   k    method silhouette
## 1 2   K-Means  0.5657356
## 2 2 K-Medoids  0.6035651
## 3 3   K-Means  0.6331908
## 4 3 K-Medoids  0.6331908
## 5 4   K-Means  0.6116777
## 6 4 K-Medoids  0.6116777
## 7 5   K-Means  0.5715312
## 8 5 K-Medoids  0.5799488

9 Kesimpulan

k = 3 adalah jumlah cluster terbaik untuk data ini.
K-Medoids sedikit lebih baik dalam stabilitas, tetapi kedua metode memberikan hasil yang hampir identik.

Dengan demikian, baik K-Means maupun K-Medoids dapat digunakan, namun k = 3 tetap menjadi pilihan optimal untuk clustering.

Tugas Akhir Individu TPG

2025-11-27

1 Library Data

2 Input Data

3 Eksplorasi Data

3.1 Mengubah Tipe Data pada Harga Menjadi Numerik

3.2 Mengurutkan Bulan dari Januari ke Desember

3.3 Cek Missing Value pada Kolom Harga secara keseluruhan

4 Filter Data Beras di Jawa Tengah Tahun 2023

4.1 Eksplorasi Data Beras Jawa Tengah 2023

4.1.1 Cek Missing Value pada Data Beras Jawa Tengah 2023

4.1.2 Statistik Deskriptif Harga Beras Jawa Tengah 2023

4.1.3 Boxplot Perbandingan Harga Premium vs Medium

5 Filter Data Beras Medium Jawa Tengah 2023

6 Eksplorasi Data Beras Medium Jawa Tengah 2023

6.1 Statistik Deskriptif per Kabupaten Kota

6.2 Grafik Tren Harga per Bulan

7 Tahapan Analisis Cluster

7.1 Data Scaling

7.2 Pemilihan K (Elbow Method)

7.3 Menentukan K pada K-Means

7.3.1 Menambahkan Hasil Cluster pada Data

7.3.2 Visualisasi Scatter Plot K-Means Clustering k=3

7.3.3 Visualisasi Bar Plot K-Means Clustering k=3

7.4 Pemilihan K pada K-Medoids (PAM)

7.5 Menentukan K pada K-Medoids (PAM)

7.5.1 Visualisasi Scatter Plot K-Medoids (PAM) Clustering k=3

7.5.2 Visualisasi Bar Plot K-Medoids (PAM) Clustering k=3

8 Perbandingan K-Means dan K-Medoids

9 Kesimpulan