Analisis Principal Component Analysis dan Clustering pada Data Kualitas Kopi Arabika

1 Pendahuluan dan Deskripsi Dataset

Laporan ini menerapkan Principal Component Analysis (PCA) dan clustering K-Means pada data penilaian sensori kopi arabika.

Sumber data. Coffee Quality Database dari Coffee Quality Institute (CQI), dikumpulkan oleh James LeDoux dan tersedia di GitHub (https://github.com/jldbc/coffee-quality-database) serta Kaggle (https://www.kaggle.com/datasets/volpatto/coffee-quality-database-from-cqi). Berkas yang dipakai: arabica_data_cleaned.csv.

Ukuran data. Berkas asli berisi 1311 observasi dan 44 kolom. Setelah membuang 1 baris yang seluruh skornya bernilai 0 (data cacat), tersisa 1310 observasi dengan 10 variabel numerik.

Variabel yang dianalisis adalah 10 skor sensori hasil cupping (skala 0-10):

Variabel	Keterangan
Aroma, Flavor, Aftertaste	Penilaian aroma, rasa, dan kesan rasa setelah ditelan
Acidity, Body, Balance	Keasaman, kekentalan, dan keseimbangan
Uniformity, Clean.Cup, Sweetness	Konsistensi antar cangkir, kebersihan cangkir, dan kemanisan
Cupper.Points	Penilaian keseluruhan oleh penguji

Variabel Total.Cup.Points tidak dipakai karena hampir sama dengan jumlah kesepuluh skor di atas sehingga akan menduplikasi informasi.

2 Tahapan Analisis

1. Eksplorasi data (statistik deskriptif, boxplot, histogram, korelasi).
2. Preprocessing: standarisasi (mean 0, simpangan baku 1).
3. PCA dan penentuan jumlah komponen utama.
4. Penentuan jumlah cluster optimal (Elbow Method dan Silhouette Score) pada skor PCA.
5. K-Means sebagai metode utama dan Hierarchical Clustering (Ward.D2) sebagai pembanding.
6. Visualisasi cluster pada ruang PCA dan interpretasi profil cluster.

3 Hasil

3.1 Eksplorasi Data Awal

statdesk <- data.frame(
  Mean   = sapply(kopi_num, mean),
  SD     = sapply(kopi_num, sd),
  Min    = sapply(kopi_num, min),
  Median = sapply(kopi_num, median),
  Max    = sapply(kopi_num, max)
)
kable(round(statdesk, 3), caption = "Statistik deskriptif skor sensori kopi arabika")

Statistik deskriptif skor sensori kopi arabika

	Mean	SD	Min	Median	Max
Aroma	7.570	0.316	5.08	7.58	8.75
Flavor	7.524	0.342	6.08	7.58	8.83
Aftertaste	7.403	0.350	6.17	7.42	8.67
Acidity	7.539	0.320	5.25	7.50	8.75
Body	7.523	0.293	5.25	7.50	8.58
Balance	7.523	0.349	6.08	7.50	8.75
Uniformity	9.841	0.489	6.00	10.00	10.00
Clean.Cup	9.841	0.722	0.00	10.00	10.00
Sweetness	9.911	0.455	1.33	10.00	10.00
Cupper.Points	7.504	0.427	5.17	7.50	10.00

boxplot(kopi_num, las = 2, col = "lightblue",
        main = "Boxplot Skor Sensori Kopi Arabika (Data Asli)")

kopi_num %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Skor") %>%
  ggplot(aes(x = Skor)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  facet_wrap(~Variabel, scales = "free") +
  labs(title = "Distribusi Skor Sensori Kopi Arabika",
       x = "Skor", y = "Frekuensi") +
  theme_minimal()

kor <- cor(kopi_num)
corrplot(kor, method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.6,
         tl.col = "black", tl.srt = 45)

Interpretasi eksplorasi. (Isi setelah melihat hasil: bentuk sebaran tiap variabel, variabel yang hampir selalu bernilai 10, keberadaan pencilan, dan kelompok variabel yang berkorelasi tinggi sebagai dasar dilakukannya PCA.)

3.2 Preprocessing: Standarisasi

Skor memiliki ragam yang berbeda antar variabel, sehingga dilakukan standarisasi agar setiap variabel berkontribusi setara dan PCA berbasis matriks korelasi.

kopi_scale <- scale(kopi_num)
boxplot(kopi_scale, las = 2, col = "lightgreen",
        main = "Boxplot Skor Sensori Setelah Standarisasi")

3.3 Principal Component Analysis

hasil_pca <- prcomp(kopi_num, scale. = TRUE)
eig <- hasil_pca$sdev^2
var_exp <- round(eig / sum(eig) * 100, 1)

tabel_eig <- data.frame(
  PC               = paste0("PC", seq_along(eig)),
  Eigenvalue       = round(eig, 4),
  Varians_Persen   = round(eig / sum(eig) * 100, 2),
  Kumulatif_Persen = round(cumsum(eig / sum(eig) * 100), 2)
)
kable(tabel_eig, caption = "Eigenvalue dan proporsi varians tiap komponen utama")

Eigenvalue dan proporsi varians tiap komponen utama

PC	Eigenvalue	Varians_Persen	Kumulatif_Persen
PC1	5.2303	52.30	52.30
PC2	1.7027	17.03	69.33
PC3	0.6153	6.15	75.48
PC4	0.5638	5.64	81.12
PC5	0.4698	4.70	85.82
PC6	0.3922	3.92	89.74
PC7	0.3707	3.71	93.45
PC8	0.2916	2.92	96.36
PC9	0.2313	2.31	98.68
PC10	0.1324	1.32	100.00

df_scree <- data.frame(PC = seq_along(eig), Eigenvalue = eig,
                       Persen = eig / sum(eig) * 100)

ggplot(df_scree, aes(x = PC, y = Persen)) +
  geom_col(fill = "steelblue", alpha = 0.8) +
  geom_line() + geom_point() +
  geom_text(aes(label = paste0(round(Persen, 1), "%")), vjust = -0.6, size = 3) +
  scale_x_continuous(breaks = df_scree$PC) +
  labs(title = "Scree Plot (Persen Varians)",
       x = "Komponen Utama", y = "Persen Varians") +
  theme_minimal()

ggplot(df_scree, aes(x = PC, y = Eigenvalue)) +
  geom_line() + geom_point(size = 2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = df_scree$PC) +
  labs(title = "Scree Plot (Eigenvalue) dengan Garis Kaiser",
       x = "Komponen Utama", y = "Eigenvalue") +
  theme_minimal()

# Kaiser (eigenvalue > 1) -> 2 PC; varians kumulatif 70-80% -> 3 PC.
# Dipakai 3 PC (ubah ke 2 bila ingin lebih ringkas; plot 3D butuh n_pc = 3).
n_pc <- 3
skor_pca <- as.data.frame(hasil_pca$x[, 1:n_pc])

Berdasarkan kriteria Kaiser, komponen dengan eigenvalue lebih dari 1 adalah 2 komponen, sedangkan 3 komponen utama pertama menjelaskan 75.48% keragaman data. Karena itu dipertahankan 3 komponen utama untuk analisis selanjutnya. PC1 menjelaskan 52.3% dan PC2 menjelaskan 17.03% keragaman.

loading <- round(hasil_pca$rotation[, 1:n_pc], 3)
kable(loading, caption = "Loading (koefisien eigenvector) komponen utama")

Loading (koefisien eigenvector) komponen utama

	PC1	PC2	PC3
Aroma	-0.348	-0.098	-0.001
Flavor	-0.404	-0.045	-0.008
Aftertaste	-0.399	-0.052	0.009
Acidity	-0.357	-0.116	-0.020
Body	-0.338	-0.159	-0.143
Balance	-0.374	-0.042	0.019
Uniformity	-0.126	0.547	0.692
Clean.Cup	-0.148	0.555	0.015
Sweetness	-0.081	0.580	-0.705
Cupper.Points	-0.362	-0.039	0.045

rot <- hasil_pca$rotation
kontribusi <- data.frame(Variabel = rownames(rot),
                         PC1 = rot[, 1]^2 * 100,
                         PC2 = rot[, 2]^2 * 100)

buat_kontribusi <- function(kolom, judul) {
  ggplot(kontribusi, aes(x = reorder(Variabel, .data[[kolom]]),
                         y = .data[[kolom]])) +
    geom_col(fill = "steelblue") +
    coord_flip() +
    geom_hline(yintercept = 100 / nrow(kontribusi),
               linetype = "dashed", color = "red") +
    labs(title = judul, x = NULL, y = "Kontribusi (%)") +
    theme_minimal()
}
buat_kontribusi("PC1", "Kontribusi Variabel terhadap PC1")

buat_kontribusi("PC2", "Kontribusi Variabel terhadap PC2")

koord <- data.frame(Variabel = rownames(rot),
                    PC1 = rot[, 1] * hasil_pca$sdev[1],
                    PC2 = rot[, 2] * hasil_pca$sdev[2])
lingkaran <- data.frame(x = cos(seq(0, 2 * pi, length.out = 200)),
                        y = sin(seq(0, 2 * pi, length.out = 200)))

ggplot() +
  geom_path(data = lingkaran, aes(x = x, y = y), color = "grey60") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey70") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey70") +
  geom_segment(data = koord, aes(x = 0, y = 0, xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.2, "cm")), color = "red") +
  geom_text(data = koord, aes(x = PC1 * 1.1, y = PC2 * 1.1, label = Variabel),
            size = 3) +
  coord_fixed(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
  labs(title = "Variabel pada Ruang PC1-PC2",
       x = paste0("PC1 (", var_exp[1], "%)"),
       y = paste0("PC2 (", var_exp[2], "%)")) +
  theme_minimal()

faktor <- 0.8 * max(abs(skor_pca[, 1:2])) / max(abs(koord[, c("PC1", "PC2")]))

ggplot() +
  geom_point(data = skor_pca, aes(x = PC1, y = PC2),
             color = "grey50", alpha = 0.3, size = 1.2) +
  geom_segment(data = koord,
               aes(x = 0, y = 0, xend = PC1 * faktor, yend = PC2 * faktor),
               arrow = arrow(length = unit(0.2, "cm")), color = "red") +
  geom_text(data = koord,
            aes(x = PC1 * faktor * 1.1, y = PC2 * faktor * 1.1, label = Variabel),
            color = "red", size = 3) +
  labs(title = "Biplot PCA",
       x = paste0("PC1 (", var_exp[1], "%)"),
       y = paste0("PC2 (", var_exp[2], "%)")) +
  theme_minimal()

Interpretasi PCA. (Isi setelah melihat tabel loading dan grafik: variabel apa yang memiliki loading terbesar pada PC1, PC2, dan PC3, lalu beri nama tiap komponen sesuai maknanya. Catatan: tanda positif atau negatif pada loading bersifat bebas dan tidak mengubah interpretasi.)

3.4 Penentuan Jumlah Cluster Optimal

Clustering dilakukan pada skor 3 komponen utama pertama.

jarak <- dist(skor_pca)

tabel_k <- data.frame(k = 1:10, WSS = NA_real_, Silhouette = NA_real_)
for (i in 1:10) {
  set.seed(123)
  km <- kmeans(skor_pca, centers = i, nstart = 25)
  tabel_k$WSS[i] <- round(km$tot.withinss, 2)
  if (i > 1) {
    tabel_k$Silhouette[i] <- round(
      mean(silhouette(km$cluster, jarak)[, "sil_width"]), 4)
  }
}
kable(tabel_k, caption = "Nilai WSS dan rata-rata Silhouette untuk tiap jumlah cluster (k)")

Nilai WSS dan rata-rata Silhouette untuk tiap jumlah cluster (k)

k	WSS	Silhouette
1	9880.79	NA
2	6014.42	0.4422
3	4523.83	0.4415
4	3453.38	0.4317
5	2923.99	0.4106
6	2472.61	0.4247
7	2174.52	0.4270
8	1927.22	0.4081
9	1765.77	0.4048
10	1563.29	0.4058

ggplot(tabel_k, aes(x = k, y = WSS)) +
  geom_line() + geom_point(size = 2) +
  scale_x_continuous(breaks = 1:10) +
  labs(title = "Elbow Method", x = "Jumlah Cluster (k)",
       y = "Total Within Sum of Squares") +
  theme_minimal()

sil_data <- tabel_k[tabel_k$k >= 2, ]
ggplot(sil_data, aes(x = k, y = Silhouette)) +
  geom_line() + geom_point(size = 2) +
  geom_vline(xintercept = sil_data$k[which.max(sil_data$Silhouette)],
             linetype = "dashed", color = "blue") +
  scale_x_continuous(breaks = 2:10) +
  labs(title = "Silhouette Score", x = "Jumlah Cluster (k)",
       y = "Rata-rata Silhouette") +
  theme_minimal()

# Tetapkan k berdasarkan Elbow dan Silhouette (lihat tabel di atas).
k_opt <- 3

Penentuan k. (Isi setelah melihat hasil: k dengan Silhouette tertinggi, ada atau tidaknya siku pada Elbow, dan alasan memilih k = 3.)

3.5 Hasil K-Means

set.seed(123)
km_final <- kmeans(skor_pca, centers = k_opt, nstart = 50)
kopi_ara$Cluster <- factor(km_final$cluster)

kable(as.data.frame(table(Cluster = kopi_ara$Cluster)),
      col.names = c("Cluster", "Jumlah Observasi"),
      caption = "Ukuran tiap cluster (K-Means)")

Ukuran tiap cluster (K-Means)

Cluster	Jumlah Observasi
1	651
2	600
3	59

Proporsi keragaman antar-cluster (between SS / total SS) sebesar 54.22%.

sil_km <- silhouette(km_final$cluster, jarak)
plot(sil_km, col = 2:(k_opt + 1), border = NA,
     main = "Silhouette Plot K-Means")

skor_plot <- cbind(skor_pca, Cluster = kopi_ara$Cluster)

ggplot(skor_plot, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(alpha = 0.6, size = 1.8) +
  stat_ellipse(level = 0.90) +
  labs(title = "Hasil K-Means pada Ruang PCA (2D)",
       x = paste0("PC1 (", var_exp[1], "%)"),
       y = paste0("PC2 (", var_exp[2], "%)")) +
  theme_minimal()

# Tampilan statis PC1 vs PC3 (selalu bisa tampil di HTML maupun PDF)
ggplot(skor_plot, aes(x = PC1, y = PC3, color = Cluster)) +
  geom_point(alpha = 0.6, size = 1.8) +
  labs(title = "Hasil K-Means pada Ruang PCA (PC1 vs PC3)",
       x = paste0("PC1 (", var_exp[1], "%)"),
       y = paste0("PC3 (", var_exp[3], "%)")) +
  theme_minimal()

# Plot 3D interaktif hanya bisa tampil pada keluaran HTML.
# Untuk PDF: jalankan chunk ini di RStudio, ambil screenshot jendela Viewer,
# simpan sebagai output/cluster_3d.png, lalu knit ke PDF.
p_3d <- plot_ly(skor_plot, x = ~PC1, y = ~PC2, z = ~PC3,
                color = ~Cluster, type = "scatter3d", mode = "markers",
                marker = list(size = 3, opacity = 0.7)) %>%
  layout(title = "Hasil K-Means pada Ruang PCA (3D)",
         scene = list(
           xaxis = list(title = paste0("PC1 (", var_exp[1], "%)")),
           yaxis = list(title = paste0("PC2 (", var_exp[2], "%)")),
           zaxis = list(title = paste0("PC3 (", var_exp[3], "%)"))))

if (knitr::is_html_output()) {
  print(p_3d)
} else if (file.exists("output/cluster_3d.png")) {
  knitr::include_graphics("output/cluster_3d.png")
} else {
  cat("[Tambahkan screenshot plot 3D sebagai output/cluster_3d.png]")
}

3.6 Profil Cluster

profil <- kopi_ara %>%
  group_by(Cluster) %>%
  summarise(n = n(), across(all_of(vars), ~ round(mean(.x), 2)))
kable(profil, caption = "Rata-rata skor sensori (data asli) tiap cluster")

Rata-rata skor sensori (data asli) tiap cluster

Cluster	n	Aroma	Flavor	Aftertaste	Acidity	Body	Balance	Uniformity	Clean.Cup	Sweetness	Cupper.Points
1	651	7.76	7.76	7.64	7.74	7.71	7.75	9.90	9.94	9.94	7.77
2	600	7.40	7.33	7.21	7.36	7.35	7.34	9.91	9.96	9.97	7.28
3	59	7.17	6.94	6.82	7.14	7.25	6.96	8.52	7.55	8.98	6.83

kopi_ara %>%
  select(Cluster, all_of(vars)) %>%
  pivot_longer(-Cluster, names_to = "Variabel", values_to = "Skor") %>%
  ggplot(aes(x = Cluster, y = Skor, fill = Cluster)) +
  geom_boxplot(outlier.size = 0.8) +
  facet_wrap(~Variabel, scales = "free_y") +
  labs(title = "Profil Skor Sensori per Cluster") +
  theme_minimal() +
  theme(legend.position = "none")

3.7 Pembanding: Hierarchical Clustering (Ward.D2)

hc <- hclust(jarak, method = "ward.D2")
plot(hc, labels = FALSE, hang = -1,
     main = "Dendrogram Hierarchical (Ward.D2)", xlab = "", sub = "")
rect.hclust(hc, k = k_opt, border = 2:(k_opt + 1))

kopi_ara$Cluster_HC <- factor(cutree(hc, k = k_opt))

kable(as.data.frame.matrix(table(KMeans = kopi_ara$Cluster,
                                 Hierarki = kopi_ara$Cluster_HC)),
      caption = "Kesesuaian keanggotaan cluster: K-Means (baris) vs Hierarchical (kolom)")

Kesesuaian keanggotaan cluster: K-Means (baris) vs Hierarchical (kolom)

1	2	3
249	389	13
0	428	172
0	0	59

sil_hc <- silhouette(cutree(hc, k = k_opt), jarak)
perbandingan <- data.frame(
  Metode     = c("K-Means", "Hierarchical (Ward.D2)"),
  Silhouette = round(c(mean(sil_km[, "sil_width"]),
                       mean(sil_hc[, "sil_width"])), 4)
)
kable(perbandingan, caption = "Perbandingan rata-rata Silhouette kedua metode")

Perbandingan rata-rata Silhouette kedua metode

Metode	Silhouette
K-Means	0.4411
Hierarchical (Ward.D2)	0.4165

ggplot(cbind(skor_pca, Cluster = kopi_ara$Cluster_HC),
       aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point(alpha = 0.6, size = 1.8) +
  labs(title = "Hierarchical (Ward.D2) pada Ruang PCA",
       x = paste0("PC1 (", var_exp[1], "%)"),
       y = paste0("PC2 (", var_exp[2], "%)")) +
  theme_minimal()

4 Diskusi

Interpretasi PCA. (Jelaskan makna tiap komponen utama berdasarkan loading, besar varians yang dijelaskan, dan alasan memilih jumlah komponen.)

Interpretasi clustering. (Jelaskan karakteristik tiap cluster berdasarkan tabel profil dan boxplot, kualitas pemisahan cluster berdasarkan Silhouette, dan perbandingan dengan Hierarchical Clustering.)

Keterbatasan. (Contoh: banyak skor mendekati 10 pada Uniformity, Clean.Cup, dan Sweetness sehingga sebarannya miring; pengaruh pencilan terhadap K-Means; hasil K-Means bergantung pada inisialisasi sehingga dipakai set.seed dan nstart.)

5 Kesimpulan

(Tuliskan 3-4 kalimat: jumlah komponen dan varians yang dijelaskan, makna komponen utama, jumlah cluster optimal beserta ciri tiap cluster, dan kesesuaian K-Means dengan Hierarchical.)

6 Lampiran: Kode Analisis

Kode analisis lengkap (R) dapat diakses melalui tautan Google Drive berikut: [tempel tautan Google Drive di sini]

sessionInfo()

## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Shanghai
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.52      plotly_4.12.1   cluster_2.1.8.2 corrplot_0.95  
## [5] ggplot2_4.0.3   tidyr_1.3.2     dplyr_1.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6        jsonlite_2.0.0      compiler_4.6.1     
##  [4] tidyselect_1.2.1    jquerylib_0.1.4     scales_1.4.0       
##  [7] yaml_2.3.12         fastmap_1.2.0       R6_2.6.1           
## [10] labeling_0.4.3      generics_0.1.4      MASS_7.3-65        
## [13] htmlwidgets_1.6.4   tibble_3.3.1        bslib_0.12.0       
## [16] pillar_1.11.1       RColorBrewer_1.1-3  rlang_1.3.0        
## [19] cachem_1.1.0        xfun_0.61           sass_0.4.10        
## [22] S7_0.2.2            otel_0.2.0          viridisLite_0.4.3  
## [25] cli_3.6.6           withr_3.0.3         magrittr_2.0.5     
## [28] crosstalk_1.2.2     digest_0.6.39       grid_4.6.1         
## [31] lifecycle_1.0.5     vctrs_0.7.3         data.table_1.18.6.1
## [34] evaluate_1.0.5      glue_1.8.1          farver_2.1.2       
## [37] httr_1.4.9          rmarkdown_2.32      purrr_1.2.2        
## [40] tools_4.6.1         pkgconfig_2.0.3     htmltools_0.5.9