Analisis Profil Gizi Makanan Berdasarkan Jenis dan Kandungan Nutrisi untuk Pembuatan Kategori Diet Seimbang pada Remaja

```Program Studi S1 Sains Data Fakultas Matematika dan Ilmu Pengetahuan Alam Universitas Negeri Surabaya —

Abstract

Remaja merupakan kelompok usia yang rentan terhadap ketidakseimbangan asupan gizi akibat pola makan yang kurang teratur dan preferensi terhadap makanan cepat saji. Penelitian ini bertujuan untuk mengelompokkan makanan berdasarkan jenis dan kandungan nutrisinya guna membentuk kategori diet seimbang yang sesuai untuk remaja. Tiga metode clustering digunakan dalam analisis ini adalah DBSCAN, Fuzzy C-Means, dan Hierarchical Clustering. Dataset yang digunakan mencakup informasi nutrisi dari berbagai jenis makanan yang umum dikonsumsi oleh remaja Indonesia, yang diolah menggunakan pendekatan analisis multivariat. Penelitian ini diharapkan dapat membantu dalam merancang rekomendasi diet yang lebih personal dan sesuai dengan kebutuhan gizi remaja, sekaligus memberikan dasar bagi program gizi yang lebih efektif.

Index Terms: clustering, DBSCAN, Fuzzy C-Means, Hierarchical, statistical analysis, nutrition, adolescent diet, food composition.

I. Introduction

A. Importance of Balanced Nutrition for Adolescents

Masa remaja merupakan periode krusial yang ditandai oleh pertumbuhan fisik dan perkembangan kognitif yang pesat. Asupan nutrisi yang adekuat sangat penting untuk mendukung proses-proses ini, meningkatkan daya tahan tubuh, dan mencegah timbulnya masalah kesehatan di masa depan. Namun, kelompok usia ini seringkali rentan terhadap ketidakseimbangan gizi. Pola makan yang tidak teratur, ditambah dengan kecenderungan mengonsumsi makanan cepat saji, dapat menyebabkan defisiensi atau kelebihan nutrisi tertentu, yang pada gilirannya dapat menghambat pertumbuhan optimal dan memengaruhi kesehatan jangka panjang.

B. Challenges of Multidimensional Nutritional Data

Profil gizi makanan memiliki karakteristik yang kompleks dan multidimensional. Setiap makanan mengandung berbagai komponen nutrien yang saling terkait, mulai dari makronutrien seperti energi, protein, lemak, dan karbohidrat, hingga mikronutrien seperti vitamin dan mineral. Keterkaitan dan banyaknya dimensi data nutrisi ini menciptakan tantangan dalam mengidentifikasi pola konsumsi gizi secara efektif. Mengurai kompleksitas ini memerlukan pendekatan analisis data yang canggih untuk mengungkap struktur tersembunyi dalam data nutrisi.

C. Role of Clustering and Multivariate Analysis in Nutritional Epidemiology

Untuk mengatasi kompleksitas data nutrisi, teknik analisis data multivariat menjadi sangat diperlukan. Metode clustering seperti DBSCAN, Fuzzy C-Means, dan Hierarchical Clustering dapat mengelompokkan makanan berdasarkan kemiripan kandungan nutrisinya. Pendekatan ini memudahkan pemahaman terhadap pola diet yang tersebar dalam populasi, karena makanan dengan profil gizi serupa akan dikelompokkan bersama. Selain itu, penggunaan analisis multivariat juga membantu mengurangi dimensi data nutrisi yang banyak dan memudahkan interpretasi hasil. Dengan mengkombinasikan teknik-teknik ini, peneliti dapat mengidentifikasi kelompok makanan yang memiliki profil gizi serupa dan mengaitkannya dengan rekomendasi diet yang tepat, khususnya bagi remaja yang membutuhkan asupan nutrisi optimal untuk mendukung pertumbuhan dan aktivitas sehari-hari.

D. Objectives of the Research

Penelitian ini bertujuan untuk menerapkan metode clustering DBSCAN, Fuzzy C-Means, dan Hierarchical Clustering untuk menganalisis data profil gizi makanan Indonesia. Tujuan utamanya adalah mengelompokkan makanan berdasarkan jenis dan kandungan nutrisinya guna membentuk kategori diet seimbang yang sesuai untuk remaja. Pendekatan ini diharapkan dapat memfasilitasi pengembangan sistem rekomendasi diet yang lebih personal dan adaptif, berdasarkan pola konsumsi sebenarnya di masyarakat. Dengan demikian, program gizi dapat menjadi lebih efektif dan tepat sasaran, mendukung kebutuhan nutrisi remaja secara lebih mendalam dan komprehensif.

II. Data and Methods

A. Dataset Description and Preprocessing

Dataset yang digunakan dalam analisis ini berasal dari Tabel Komposisi Pangan Indonesia (TKPI) tahun 2020, yang diterbitkan oleh Kementerian Kesehatan Republik Indonesia. Dataset ini memuat informasi komposisi nutrisi dari berbagai jenis makanan yang umum dikonsumsi di Indonesia, mencakup makronutrien (energi, protein, lemak, karbohidrat) serta mikronutrien (vitamin dan mineral).

Pada saat pemuatan awal, dataset memiliki dimensi 1147 baris (item makanan) dan 22 kolom (variabel nutrisi). Namun, ditemukan adanya missing value yang cukup signifikan. Oleh karena itu, variabel dengan lebih dari 25% missing value dikeluarkan dari analisis. Variabel yang dipertahankan adalah: NATRIUM, SERAT, THIAMIN, BESI, KALSIUM, FOSFOR, LEMAK, PROTEIN, ABU, AIR, ENERGI, dan KH. Imputasi dilakukan untuk mengisi missing value yang tersisa dan scaling dilakukan terhadap data numerik agar rentangnya seragam.

Table 1. Summary of Dataset Characteristics and Preprocessing Steps

Characteristic / Step	Description / Value
Initial Data Dimensions	1147 Rows x 22 Columns
Variables with >25% Missing Values (Excluded)	RETINOL, KAR_TOTAL, B_KAR, TEMBAGA, RIBOFLAVIN, SENG, NIASIN, VIT_C
Variables Retained for Analysis	NATRIUM, SERAT, THIAMIN, BESI, KALSIUM, FOSFOR, LEMAK, PROTEIN, ABU, AIR, ENERGI, KH
Missing Values After Imputation	0
Final Data Dimensions	1147 Rows x 12 Columns

B. Clustering Algorithms

Dalam penelitian ini, tiga algoritma clustering yang berbeda diterapkan untuk mengelompokkan data profil gizi makanan yang telah dipreprocessing menggunakan Hierarchical Clustering (HC), Fuzzy C-Means (FCM), dan DBSCAN. Pemilihan ketiga algoritma ini merefleksikan pendekatan komprehensif untuk mengeksplorasi struktur data nutrisi yang mendasari. Ini mengakui bahwa karakteristik data yang berbeda mungkin paling baik ditangkap oleh filosofi algoritmik yang berbeda, sehingga memungkinkan perbandingan kinerja yang mendalam.

1. Hierarchical Clustering (HC): Algoritma ini membangun hirarki cluster baik secara aglomeratif (bottom-up, dimulai dengan setiap titik data sebagai cluster individu dan menggabungkannya) atau divisif (top-down, dimulai dengan satu cluster besar dan membaginya). Untuk studi ini, pendekatan aglomeratif umumnya digunakan. Penentuan jumlah cluster optimal (k) untuk HC adalah 5.[1] Kriteria linkage (misalnya, metode Ward, complete, average) menentukan bagaimana jarak antar cluster dihitung, dan metode spesifik akan dijelaskan dalam implementasi teknis.
2. Fuzzy C-Means (FCM): Berbeda dengan clustering tradisional yang menempatkan setiap titik data secara eksklusif ke dalam satu cluster, FCM adalah algoritma soft clustering di mana setiap titik data dapat memiliki derajat keanggotaan tertentu pada beberapa cluster. Ini sangat berguna ketika titik data mungkin secara alami berada di antara kategori yang berbeda. Parameter fuzziness (m) dalam FCM mengontrol tingkat keburaman; nilai m yang lebih tinggi menghasilkan cluster yang lebih fuzzy.
3. DBSCAN (Density-Based Spatial Clustering of Applications with Noise): DBSCAN adalah algoritma clustering berbasis kepadatan yang mengelompokkan titik-titik yang berdekatan satu sama lain, dan menandai titik-titik yang berada sendiri di daerah berdensitas rendah sebagai outlier (noise). Keunggulan DBSCAN adalah tidak memerlukan penentuan jumlah cluster di awal. Namun, algoritma ini membutuhkan dua parameter: eps (epsilon), yang mendefinisikan radius maksimum lingkungan, dan MinPts (minimum points), yang mendefinisikan jumlah minimum titik yang diperlukan untuk membentuk daerah padat.

C. Evaluation Metrics

Untuk mengevaluasi kualitas hasil clustering dan representasi data, dua metrik utama digunakan: 1. Silhouette Score: Metrik ini digunakan untuk mengukur seberapa mirip suatu objek dengan _cluster_nya sendiri (kohesi) dibandingkan dengan cluster lain (pemisahan). Skor berkisar dari -1 hingga +1, di mana nilai tinggi menunjukkan bahwa objek sangat cocok dengan _cluster_nya sendiri dan tidak cocok dengan cluster tetangga. Metrik ini digunakan untuk membandingkan kinerja Hierarchical Clustering dan Fuzzy C-Means. 2. Multidimensional Scaling (MDS): MDS adalah teknik visualisasi yang digunakan untuk merepresentasikan kemiripan kasus individu dari dataset. Tujuannya adalah untuk menggambarkan jarak antar titik data dalam ruang berdimensi tinggi seakurat mungkin dalam ruang berdimensi lebih rendah (biasanya 2D atau 3D). Nilai “stress” dalam MDS mengukur perbedaan antara jarak dalam ruang berdimensi tinggi asli dan jarak dalam representasi berdimensi rendah. Nilai stress yang lebih rendah menunjukkan kecocokan yang lebih baik.

Penggunaan kombinasi metrik kuantitatif kualitas clustering (Silhouette Score) dan teknik visualisasi dengan ukuran kecocokan (MDS Stress) memberikan evaluasi yang kuat dan multi-aspek terhadap hasil clustering. Ini memastikan bahwa cluster yang dipilih tidak hanya valid secara statistik tetapi juga dapat diinterpretasikan dalam ruang berdimensi rendah, yang penting untuk memahami dan mengomunikasikan pola nutrisi yang kompleks.

III. Implementasi Kode R

# ========================================
# R SCRIPT ANALISIS PROFIL GIZI MAKANAN REMAJA
# ========================================
# Script ini melakukan analisis multivariate pada data gizi makanan
# berdasarkan Rencana Analisis Ideal yang diberikan.

# ========================================
# 1. INSTALL & LOAD LIBRARIES
# ========================================
# Pastikan semua library yang dibutuhkan terinstall
packages <- c("cluster", "factoextra", "dbscan", "psych",
              "MASS","tidyverse", "corrplot", "naniar", "vegan", "e1071", "reshape2",
              "ggpubr", "clusterSim", "fpc", "RColorBrewer", "ggrepel")

# Cek dan install paket yang belum ada
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Load semua library
lapply(packages, library, character.only = TRUE)

## [[1]]
## [1] "cluster"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
##  [1] "factoextra" "ggplot2"    "cluster"    "stats"      "graphics"  
##  [6] "grDevices"  "utils"      "datasets"   "methods"    "base"      
## 
## [[3]]
##  [1] "dbscan"     "factoextra" "ggplot2"    "cluster"    "stats"     
##  [6] "graphics"   "grDevices"  "utils"      "datasets"   "methods"   
## [11] "base"      
## 
## [[4]]
##  [1] "psych"      "dbscan"     "factoextra" "ggplot2"    "cluster"   
##  [6] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
## [11] "methods"    "base"      
## 
## [[5]]
##  [1] "MASS"       "psych"      "dbscan"     "factoextra" "ggplot2"   
##  [6] "cluster"    "stats"      "graphics"   "grDevices"  "utils"     
## [11] "datasets"   "methods"    "base"      
## 
## [[6]]
##  [1] "lubridate"  "forcats"    "stringr"    "dplyr"      "purrr"     
##  [6] "readr"      "tidyr"      "tibble"     "tidyverse"  "MASS"      
## [11] "psych"      "dbscan"     "factoextra" "ggplot2"    "cluster"   
## [16] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
## [21] "methods"    "base"      
## 
## [[7]]
##  [1] "corrplot"   "lubridate"  "forcats"    "stringr"    "dplyr"     
##  [6] "purrr"      "readr"      "tidyr"      "tibble"     "tidyverse" 
## [11] "MASS"       "psych"      "dbscan"     "factoextra" "ggplot2"   
## [16] "cluster"    "stats"      "graphics"   "grDevices"  "utils"     
## [21] "datasets"   "methods"    "base"      
## 
## [[8]]
##  [1] "naniar"     "corrplot"   "lubridate"  "forcats"    "stringr"   
##  [6] "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"    
## [11] "tidyverse"  "MASS"       "psych"      "dbscan"     "factoextra"
## [16] "ggplot2"    "cluster"    "stats"      "graphics"   "grDevices" 
## [21] "utils"      "datasets"   "methods"    "base"      
## 
## [[9]]
##  [1] "vegan"      "lattice"    "permute"    "naniar"     "corrplot"  
##  [6] "lubridate"  "forcats"    "stringr"    "dplyr"      "purrr"     
## [11] "readr"      "tidyr"      "tibble"     "tidyverse"  "MASS"      
## [16] "psych"      "dbscan"     "factoextra" "ggplot2"    "cluster"   
## [21] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
## [26] "methods"    "base"      
## 
## [[10]]
##  [1] "e1071"      "vegan"      "lattice"    "permute"    "naniar"    
##  [6] "corrplot"   "lubridate"  "forcats"    "stringr"    "dplyr"     
## [11] "purrr"      "readr"      "tidyr"      "tibble"     "tidyverse" 
## [16] "MASS"       "psych"      "dbscan"     "factoextra" "ggplot2"   
## [21] "cluster"    "stats"      "graphics"   "grDevices"  "utils"     
## [26] "datasets"   "methods"    "base"      
## 
## [[11]]
##  [1] "reshape2"   "e1071"      "vegan"      "lattice"    "permute"   
##  [6] "naniar"     "corrplot"   "lubridate"  "forcats"    "stringr"   
## [11] "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"    
## [16] "tidyverse"  "MASS"       "psych"      "dbscan"     "factoextra"
## [21] "ggplot2"    "cluster"    "stats"      "graphics"   "grDevices" 
## [26] "utils"      "datasets"   "methods"    "base"      
## 
## [[12]]
##  [1] "ggpubr"     "reshape2"   "e1071"      "vegan"      "lattice"   
##  [6] "permute"    "naniar"     "corrplot"   "lubridate"  "forcats"   
## [11] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"     
## [16] "tibble"     "tidyverse"  "MASS"       "psych"      "dbscan"    
## [21] "factoextra" "ggplot2"    "cluster"    "stats"      "graphics"  
## [26] "grDevices"  "utils"      "datasets"   "methods"    "base"      
## 
## [[13]]
##  [1] "clusterSim" "ggpubr"     "reshape2"   "e1071"      "vegan"     
##  [6] "lattice"    "permute"    "naniar"     "corrplot"   "lubridate" 
## [11] "forcats"    "stringr"    "dplyr"      "purrr"      "readr"     
## [16] "tidyr"      "tibble"     "tidyverse"  "MASS"       "psych"     
## [21] "dbscan"     "factoextra" "ggplot2"    "cluster"    "stats"     
## [26] "graphics"   "grDevices"  "utils"      "datasets"   "methods"   
## [31] "base"      
## 
## [[14]]
##  [1] "fpc"        "clusterSim" "ggpubr"     "reshape2"   "e1071"     
##  [6] "vegan"      "lattice"    "permute"    "naniar"     "corrplot"  
## [11] "lubridate"  "forcats"    "stringr"    "dplyr"      "purrr"     
## [16] "readr"      "tidyr"      "tibble"     "tidyverse"  "MASS"      
## [21] "psych"      "dbscan"     "factoextra" "ggplot2"    "cluster"   
## [26] "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
## [31] "methods"    "base"      
## 
## [[15]]
##  [1] "RColorBrewer" "fpc"          "clusterSim"   "ggpubr"       "reshape2"    
##  [6] "e1071"        "vegan"        "lattice"      "permute"      "naniar"      
## [11] "corrplot"     "lubridate"    "forcats"      "stringr"      "dplyr"       
## [16] "purrr"        "readr"        "tidyr"        "tibble"       "tidyverse"   
## [21] "MASS"         "psych"        "dbscan"       "factoextra"   "ggplot2"     
## [26] "cluster"      "stats"        "graphics"     "grDevices"    "utils"       
## [31] "datasets"     "methods"      "base"        
## 
## [[16]]
##  [1] "ggrepel"      "RColorBrewer" "fpc"          "clusterSim"   "ggpubr"      
##  [6] "reshape2"     "e1071"        "vegan"        "lattice"      "permute"     
## [11] "naniar"       "corrplot"     "lubridate"    "forcats"      "stringr"     
## [16] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [21] "tidyverse"    "MASS"         "psych"        "dbscan"       "factoextra"  
## [26] "ggplot2"      "cluster"      "stats"        "graphics"     "grDevices"   
## [31] "utils"        "datasets"     "methods"      "base"

II. Load Data

# ========================================
# 2. LOAD DATA                                    
# ========================================

# Pastikan file dataset berada di direktori kerja atau gunakan path lengkap
file_path <- "ANMUL - DATASET PANGAN.csv"

# Periksa apakah file ada
if (!file.exists(file_path)) {
  stop("File dataset tidak ditemukan. Pastikan path file sudah benar.")
}

# Membaca data, melewati 2 baris pertama (jika ada header tambahan)
data_raw <- read.csv(file_path, skip = 2, stringsAsFactors = FALSE, na.strings = c("NA", "", "-"))

# Rename kolom sesuai dengan struktur data
colnames(data_raw) <- c("KODE", "NAMA_BAHAN", "SUMBER", "AIR", "ENERGI", "PROTEIN",
                        "LEMAK", "KH", "SERAT", "ABU", "KALSIUM", "FOSFOR", "BESI",
                        "NATRIUM", "KALIUM", "TEMBAGA", "ZINC", "RETINOL", "B_KAR",
                        "KAR_TOTAL", "THIAMIN", "RIBOFLAVIN", "NIASIN", "VIT_C", "BDD")

# Memproses kolom nutrien menjadi numerik
nutrien_cols_indices <- 4:24
nutrien_cols_names <- colnames(data_raw)[nutrien_cols_indices]

data_processed <- data_raw %>%
  mutate(across(all_of(nutrien_cols_names), ~ suppressWarnings(as.numeric(as.character(.)))))

# Menyiapkan data numerik untuk analisis
data_numeric <- data_processed %>%
  dplyr::select(NAMA_BAHAN, all_of(nutrien_cols_names))

# Tampilkan dimensi data sebagai validasi
cat("Dimensi data setelah load:", dim(data_numeric), "\n")

## Dimensi data setelah load: 1147 22

# ========================================
# 3. PREPROCESSING DATA
# ========================================
data_numeric_nutrients_only <- data_numeric %>% dplyr::select(-NAMA_BAHAN)
missing_summary <- miss_var_summary(data_numeric_nutrients_only)
print("Persentase Missing Value per Variabel (Hanya Nutrisi):")

## [1] "Persentase Missing Value per Variabel (Hanya Nutrisi):"

print(missing_summary)

## # A tibble: 21 × 3
##    variable   n_miss pct_miss
##    <chr>       <int>    <num>
##  1 RETINOL       607     52.9
##  2 KAR_TOTAL     542     47.3
##  3 B_KAR         472     41.2
##  4 TEMBAGA       356     31.0
##  5 RIBOFLAVIN    353     30.8
##  6 ZINC          350     30.5
##  7 NIASIN        344     30.0
##  8 KALIUM        330     28.8
##  9 VIT_C         308     26.9
## 10 NATRIUM       286     24.9
## # ℹ 11 more rows

vars_to_keep_nutrients <- missing_summary %>%
  filter(pct_miss < 25) %>%
  pull(variable)
data_filtered <- data_numeric %>%
  dplyr::select(NAMA_BAHAN, all_of(vars_to_keep_nutrients))
cat("\nVariabel yang digunakan setelah filtering missing value (>25%):")

## 
## Variabel yang digunakan setelah filtering missing value (>25%):

kept_vars_after_missing <- colnames(data_filtered %>% dplyr::select(-NAMA_BAHAN))
print(kept_vars_after_missing)

##  [1] "NATRIUM" "SERAT"   "THIAMIN" "BESI"    "KALSIUM" "FOSFOR"  "LEMAK"  
##  [8] "PROTEIN" "ABU"     "AIR"     "ENERGI"  "KH"

cat("Dimensi data setelah filtering variabel:", dim(data_filtered), "\n")

## Dimensi data setelah filtering variabel: 1147 13

# **Penting:** Pastikan nama_bahan selalu diambil dari data_filtered yang sudah bersih
data_for_analysis <- data_filtered %>% dplyr::select(-NAMA_BAHAN)
nama_bahan <- data_filtered$NAMA_BAHAN # Dipastikan diambil dari data_filtered


# Melihat distribusi data numerik setelah filtering variabel
cat("\nMembuat plot distribusi untuk variabel numerik...\n")

## 
## Membuat plot distribusi untuk variabel numerik...

plot_distribusi <- data_for_analysis %>%
  dplyr::select(where(is.numeric)) %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai") %>%
  ggplot(aes(x = Nilai)) +
  geom_histogram(bins = 30, fill = "skyblue", color = "black") +
  facet_wrap(~Variabel, scales = "free") +
  theme_minimal() +
  labs(title = "Distribusi Data Numerik Setelah Filtering Variabel", x = "Nilai", y = "Frekuensi")
print(plot_distribusi)

## Warning: Removed 748 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Menyimpan plot distribusi
ggsave("distribusi_data_numerik.png", plot_distribusi, width = 12, height = 8)

## Warning: Removed 748 rows containing non-finite outside the scale range
## (`stat_bin()`).

cat("Plot distribusi data numerik disimpan sebagai 'distribusi_data_numerik.png'.\n")

## Plot distribusi data numerik disimpan sebagai 'distribusi_data_numerik.png'.

# --- Perbaikan: Tangani nilai Inf, -Inf, dan NaN sebelum imputasi ---
# Mengubah nilai Inf, -Inf, dan NaN menjadi NA agar bisa diimputasi median
data_for_analysis <- data_for_analysis %>%
  mutate(across(where(is.numeric), ~na_if(., Inf))) %>%
  mutate(across(where(is.numeric), ~na_if(., -Inf))) %>%
  mutate(across(where(is.numeric), ~ifelse(is.nan(.), NA, .)))
# --- Akhir Perbaikan ---

if (any(is.na(data_for_analysis))) {
  cat("\nMelakukan imputasi median untuk missing values...\n")
  data_imputed <- data_for_analysis %>%
    mutate(across(everything(), ~ ifelse(is.na(.), median(., na.rm = TRUE), .)))
  cat("Jumlah NA setelah imputasi:", sum(is.na(data_imputed)), "\n")
} else {
  cat("\nTidak ada missing value yang perlu diimputasi setelah filtering variabel.\n")
  data_imputed <- data_for_analysis
}

## 
## Melakukan imputasi median untuk missing values...
## Jumlah NA setelah imputasi: 0

cat("\nMelakukan scaling data (Z-score)...\n")

## 
## Melakukan scaling data (Z-score)...

data_scaled_original <- as.data.frame(scale(data_imputed))
# data_scaled_original ini yang akan langsung kita gunakan untuk clustering
cat("Dimensi data setelah scaling awal:", dim(data_scaled_original), "\n")

## Dimensi data setelah scaling awal: 1147 12

# ========================================
# 4. CLUSTERING (Menggunakan Data Asli yang Sudah Skala)
# ========================================
# Inisialisasi variabel hasil clustering
# Penting: Sesuaikan panjang rep() dengan jumlah baris data_scaled_original
clusters_hc <- rep(NA, nrow(data_scaled_original))
clusters_fcm <- rep(NA, nrow(data_scaled_original))
clusters_db <- rep(NA, nrow(data_scaled_original))
k_optimal_hc <- NA_integer_
k_fcm <- NA_integer_
eps_chosen <- NA_real_
min_pts_dbscan_val <- NA_integer_
dist_mat <- NULL # Mengganti nama dari dist_mat_pca agar lebih jelas

# Pastikan data_for_clustering menggunakan data yang sudah diskalakan
if (nrow(data_scaled_original) > 0 && ncol(data_scaled_original) > 0) {
  cat("\nMemulai Analisis Clustering pada data asli yang sudah di-scaled...\n")
  data_for_clustering <- data_scaled_original # Perubahan penting di sini!
  
# 4.1. Hierarchical Clustering (HC)
  cat("\n--- Hierarchical Clustering (HC) ---\n")
  dist_mat <- dist(data_for_clustering, method = "euclidean") # dihitung dari data_for_clustering yang baru
  hc_ward <- hclust(dist_mat, method = "ward.D2")
  print("Membuat Dendrogram HC...")
  
# Simpan Dendrogram HC
  png("dendrogram_hc.png", width = 1000, height = 600, res = 100)
  plot(hc_ward, main = "Dendrogram HC (Ward.D2)", xlab = "Observasi", ylab = "Height", sub = "", cex = 0.6, labels = FALSE)
  dev.off()
  cat("Dendrogram HC disimpan sebagai 'dendrogram_hc.png'.\n")
  
  cat("\nMengevaluasi jumlah cluster optimal (k=2 to 5) untuk HC...\n")
  sil_scores_hc <- numeric(4)
  k_options_hc <- 2:5
  valid_k_found_hc <- FALSE
  if (nrow(data_for_clustering) >= max(k_options_hc) +1) { # +1 agar cukup data untuk cluster terakhir
    for (i_hc in seq_along(k_options_hc)) {
      k_hc <- k_options_hc[i_hc]
      if (nrow(data_for_clustering) < k_hc) { # Cek jika data kurang dari k
        cat("Data (",nrow(data_for_clustering),") < k (",k_hc,") untuk Silhouette. Lewati.\n")
        sil_scores_hc[i_hc] <- NA
        next
      }
      clusters_temp_hc <- cutree(hc_ward, k = k_hc)
      if (length(unique(clusters_temp_hc)) > 1 && min(table(clusters_temp_hc)) > 1) {
        sil_hc <- silhouette(clusters_temp_hc, dist_mat) # Menggunakan dist_mat yang baru
        sil_scores_hc[i_hc] <- mean(sil_hc[, 3])
        cat("Rata-rata Silhouette Score untuk k =", k_hc, ":", sil_scores_hc[i_hc], "\n")
        # print(fviz_silhouette(sil_hc) + ggtitle(paste("Silhouette Plot HC, k=", k_hc)))
        valid_k_found_hc <- TRUE
      } else {
        cat("Tidak dapat menghitung Silhouette untuk k =", k_hc, "(singleton cluster atau hanya 1 cluster).\n")
        sil_scores_hc[i_hc] <- NA
      }
    }
    if(valid_k_found_hc && any(!is.na(sil_scores_hc))){
      k_optimal_hc <- k_options_hc[which.max(sil_scores_hc[!is.na(sil_scores_hc)])] # Ambil max dari yg tidak NA
    } else {
      k_optimal_hc <- 3 # Default jika tidak ada yang valid
      cat("Tidak ada k optimal valid dari Silhouette, menggunakan default k=", k_optimal_hc, "\n")
    }
  } else {
    k_optimal_hc <- min(c(2, if(nrow(data_for_clustering)>1) nrow(data_for_clustering) -1 else 2), na.rm=TRUE)
    if(is.na(k_optimal_hc) || k_optimal_hc < 2) k_optimal_hc = 2
    cat("Data terlalu sedikit untuk evaluasi k=2-5, menggunakan k=", k_optimal_hc, "\n")
  }
  cat("\nJumlah cluster optimal (k) yang dipilih untuk HC:", k_optimal_hc, "\n")
  if (!is.na(k_optimal_hc) && k_optimal_hc >=2 && nrow(data_for_clustering) >= k_optimal_hc) {
    clusters_hc <- cutree(hc_ward, k = k_optimal_hc)
  } else {
    clusters_hc <- rep(NA, nrow(data_for_clustering)) # Sesuaikan panjangnya
  }
  
# 4.2. Fuzzy C-Means (FCM)
  cat("\n--- Fuzzy C-Means (FCM) ---\n")
  k_fcm <- k_optimal_hc # Menggunakan k optimal dari HC
  cat("Menggunakan k =", k_fcm, "untuk FCM (berdasarkan hasil HC).\n")
  if (!is.na(k_fcm) && k_fcm >= 2 && nrow(data_for_clustering) >= k_fcm) {
    fcm_result <- cmeans(data_for_clustering, centers = k_fcm, iter.max = 100, m = 2, method = "cmeans", dist = "euclidean")
    clusters_fcm <- fcm_result$cluster
    membership_matrix <- fcm_result$membership
    # (Evaluasi membership bisa ditambahkan kembali jika diperlukan)
  } else {
    cat("Tidak dapat melakukan FCM, k_fcm tidak valid atau data kurang.\n")
    clusters_fcm <- rep(NA, nrow(data_for_clustering))
  }
  
# 4.3. DBSCAN
  cat("\n--- DBSCAN ---\n")
  data_dbscan_input <- data_for_clustering # DBSCAN juga pakai data_for_clustering yang baru
  
  potential_min_pts_db <- ncol(data_dbscan_input) + 1
  min_pts_dbscan_val <- max(2, potential_min_pts_db)
  cat("Menggunakan MinPts =", min_pts_dbscan_val, "\n")
  
  k_for_knnplot_db <- min_pts_dbscan_val - 1
  eps_chosen <- 1.0 # Default fallback
  
  if (k_for_knnplot_db <= 0) {
    cat("Nilai k untuk kNNdistplot (MinPts - 1) adalah <= 0. Tidak dapat membuat plot.\n")
  } else if (nrow(data_dbscan_input) <= k_for_knnplot_db) {
    cat("Jumlah observasi (", nrow(data_dbscan_input), ") tidak cukup untuk k =", k_for_knnplot_db, "pada kNNdistplot.\n")
  } else {
    # Simpan kNNdistplot
    png("kNNdistplot_dbscan.png", width = 800, height = 600, res = 100)
    kNNdistplot(data_dbscan_input, k = k_for_knnplot_db)
    kth_nn_distances_db <- kNNdist(data_dbscan_input, k = k_for_knnplot_db)
    if (length(kth_nn_distances_db) > 0 && !all(is.na(kth_nn_distances_db))) {
      median_kth_dist_db <- median(kth_nn_distances_db, na.rm = TRUE)
      abline(h = median_kth_dist_db, col = "red", lty = 2)
      cat("Perhatikan 'elbow' pada plot kNNdistplot.\n")
      eps_chosen <- median_kth_dist_db
      cat("\nMenggunakan eps =", round(eps_chosen,2) ,", estimasi dari median k-NN distance. Sesuaikan berdasarkan plot.\n")
    } else {
      cat("Tidak dapat menghitung jarak kNN. Menggunakan eps default=", eps_chosen, "\n")
    }
    dev.off()
    cat("kNNdistplot disimpan sebagai 'kNNdistplot_dbscan.png'.\n")
  }
  
  if (!is.na(eps_chosen) && !is.na(min_pts_dbscan_val) && eps_chosen > 0 && min_pts_dbscan_val > 0) {
    if(min_pts_dbscan_val > nrow(data_dbscan_input)){
      cat("MinPts (", min_pts_dbscan_val, ") > data (", nrow(data_dbscan_input), "). DBSCAN tidak dijalankan.\n")
      clusters_db <- rep(NA, nrow(data_dbscan_input))
    } else {
      dbscan_res <- dbscan::dbscan(data_dbscan_input, eps = eps_chosen, minPts = min_pts_dbscan_val)
      clusters_db <- dbscan_res$cluster
    }
  } else {
    cat("Parameter eps atau MinPts tidak valid. DBSCAN tidak dijalankan.\n")
    clusters_db <- rep(NA, nrow(data_dbscan_input))
  }
} else {
  cat("Data tidak tersedia atau kosong. Clustering tidak dilakukan.\n")
  # Pastikan panjang clusters_xx cocok dengan nama_bahan yang tidak di PCA
  clusters_hc <- rep(NA, length(nama_bahan))
  clusters_fcm <- rep(NA, length(nama_bahan))
  clusters_db <- rep(NA, length(nama_bahan))
}

## 
## Memulai Analisis Clustering pada data asli yang sudah di-scaled...
## 
## --- Hierarchical Clustering (HC) ---
## [1] "Membuat Dendrogram HC..."

## Dendrogram HC disimpan sebagai 'dendrogram_hc.png'.
## 
## Mengevaluasi jumlah cluster optimal (k=2 to 5) untuk HC...
## Rata-rata Silhouette Score untuk k = 2 : 0.3530866 
## Rata-rata Silhouette Score untuk k = 3 : 0.369881 
## Rata-rata Silhouette Score untuk k = 4 : 0.3837861 
## Rata-rata Silhouette Score untuk k = 5 : 0.3863501 
## 
## Jumlah cluster optimal (k) yang dipilih untuk HC: 5 
## 
## --- Fuzzy C-Means (FCM) ---
## Menggunakan k = 5 untuk FCM (berdasarkan hasil HC).
## 
## --- DBSCAN ---
## Menggunakan MinPts = 13

## Perhatikan 'elbow' pada plot kNNdistplot.
## 
## Menggunakan eps = 0.94 , estimasi dari median k-NN distance. Sesuaikan berdasarkan plot.
## kNNdistplot disimpan sebagai 'kNNdistplot_dbscan.png'.

# Cek jumlah cluster dan noise DBSCAN
cat("\nHasil DBSCAN Clustering (0 = Noise):")

## 
## Hasil DBSCAN Clustering (0 = Noise):

print(table(clusters_db, useNA = "ifany"))

## clusters_db
##   0   1 
## 425 722

num_noise_dbscan <- sum(clusters_db == 0, na.rm = TRUE)
num_clusters_dbscan <- length(unique(clusters_db[!is.na(clusters_db) & clusters_db != 0]))
cat("Jumlah Noise Points:", num_noise_dbscan, "\n")

## Jumlah Noise Points: 425

cat("Jumlah Cluster Ditemukan (tidak termasuk noise):", num_clusters_dbscan, "\n")

## Jumlah Cluster Ditemukan (tidak termasuk noise): 1

if(!is.na(eps_chosen) && !is.na(min_pts_dbscan_val)) {
  cat("Saran DBSCAN: Coba eps sekitar", round(eps_chosen,2), "dan MinPts", min_pts_dbscan_val, ".\n")
}

## Saran DBSCAN: Coba eps sekitar 0.94 dan MinPts 13 .

# ========================================
# 5. EVALUASI CLUSTERING
# ========================================
cat("\nMemulai Evaluasi Hasil Clustering...\n")

## 
## Memulai Evaluasi Hasil Clustering...

# 5.1. Silhouette Score
# HC
cat("\nMenghitung Silhouette Score untuk HC (k=", k_optimal_hc, ")...\n")

## 
## Menghitung Silhouette Score untuk HC (k= 5 )...

mean_sil_hc <- NA
if (!is.na(k_optimal_hc) && length(clusters_hc) == nrow(data_for_clustering) && !all(is.na(clusters_hc)) && !is.null(dist_mat)) { # Menggunakan dist_mat
  valid_hc_indices_sil <- !is.na(clusters_hc)
  clusters_hc_valid_sil <- clusters_hc[valid_hc_indices_sil]
  if (length(unique(clusters_hc_valid_sil)) > 1 && min(table(clusters_hc_valid_sil)) > 1 && length(clusters_hc_valid_sil) > 1) {
    data_hc_sil_valid <- data_for_clustering[valid_hc_indices_sil, , drop = FALSE]
    if(nrow(data_hc_sil_valid) > 1){
      dist_mat_hc_sil_valid <- dist(data_hc_sil_valid)
      sil_hc_final <- silhouette(clusters_hc_valid_sil, dist_mat_hc_sil_valid)
      mean_sil_hc <- mean(sil_hc_final[, 3])
      cat("Rata-rata Silhouette Score HC:", mean_sil_hc, "\n")
    } else {cat("Data HC valid < 2 untuk Silhouette.\n")}
  } else { cat("Tidak bisa hitung Silhouette HC (singleton/1 cluster setelah filter NA).\n")}
} else { cat("Tidak bisa memulai Silhouette HC (k, clusters_hc, data_for_clustering, atau dist_mat tidak valid).\n")}

## Rata-rata Silhouette Score HC: 0.3863501

# FCM
cat("\nMenghitung Silhouette Score untuk FCM (k=", k_fcm, ")...\n")

## 
## Menghitung Silhouette Score untuk FCM (k= 5 )...

mean_sil_fcm <- NA
if (!is.na(k_fcm) && length(clusters_fcm) == nrow(data_for_clustering) && !all(is.na(clusters_fcm)) && !is.null(dist_mat)) { # Menggunakan dist_mat
  valid_fcm_indices_sil <- !is.na(clusters_fcm)
  clusters_fcm_valid_sil <- clusters_fcm[valid_fcm_indices_sil]
  if (length(unique(clusters_fcm_valid_sil)) > 1 && min(table(clusters_fcm_valid_sil)) > 1 && length(clusters_fcm_valid_sil) > 1) {
    data_fcm_sil_valid <- data_for_clustering[valid_fcm_indices_sil, , drop = FALSE]
    if(nrow(data_fcm_sil_valid) > 1){
      dist_mat_fcm_sil_valid <- dist(data_fcm_sil_valid)
      sil_fcm <- silhouette(clusters_fcm_valid_sil, dist_mat_fcm_sil_valid)
      mean_sil_fcm <- mean(sil_fcm[, 3])
      cat("Rata-rata Silhouette Score FCM:", mean_sil_fcm, "\n")
    } else {cat("Data FCM valid < 2 untuk Silhouette.\n")}
  } else { cat("Tidak bisa hitung Silhouette FCM (singleton/1 cluster setelah filter NA).\n")}
} else { cat("Tidak bisa memulai Silhouette FCM (k, clusters_fcm, data_for_clustering, atau dist_mat tidak valid).\n")}

## Rata-rata Silhouette Score FCM: 0.1911777

# DBSCAN
cat("\nMenghitung Silhouette Score untuk DBSCAN (hanya non-noise points)...\n")

## 
## Menghitung Silhouette Score untuk DBSCAN (hanya non-noise points)...

mean_sil_dbscan <- NA
# data_dbscan_input adalah data_for_clustering yang baru
if (exists("clusters_db") && length(clusters_db) == nrow(data_dbscan_input) && !all(is.na(clusters_db))) {
  non_noise_indices_db_sil <- which(clusters_db != 0 & !is.na(clusters_db))
  if (length(non_noise_indices_db_sil) > 1) {
    clusters_db_non_noise_sil <- clusters_db[non_noise_indices_db_sil]
    if (length(unique(clusters_db_non_noise_sil)) > 1 && min(table(clusters_db_non_noise_sil)) > 1) {
      data_dbscan_non_noise_sil <- data_dbscan_input[non_noise_indices_db_sil, , drop = FALSE]
      if (nrow(data_dbscan_non_noise_sil) > 1) {
        dist_mat_dbscan_non_noise_sil <- dist(data_dbscan_non_noise_sil, method = "euclidean")
        sil_dbscan <- silhouette(clusters_db_non_noise_sil, dist_mat_dbscan_non_noise_sil)
        mean_sil_dbscan <- mean(sil_dbscan[, 3])
        cat("Rata-rata Silhouette Score DBSCAN (non-noise):", mean_sil_dbscan, "\n")
      } else {cat("Data DBSCAN non-noise valid < 2 untuk Silhouette.\n")}
    } else { cat("Tidak bisa hitung Silhouette DBSCAN (kurang dari 2 cluster non-noise atau singleton).\n")}
  } else { cat("Tidak cukup titik non-noise untuk Silhouette DBSCAN.\n")}
} else { cat("Tidak bisa memulai Silhouette DBSCAN (clusters_db atau data_dbscan_input tidak valid).\n")}

## Tidak bisa hitung Silhouette DBSCAN (kurang dari 2 cluster non-noise atau singleton).

# 5.2. Visualisasi Cluster pada Plot
# Karena tidak ada PCA, kita tidak bisa lagi menggunakan fviz_pca_ind.
# Kita bisa membuat plot MDS dengan cluster, atau plot 2D dari nutrisi pilihan.
cat("\nMembuat visualisasi cluster pada plot...\n")

## 
## Membuat visualisasi cluster pada plot...

# ========================================
# 6. PERCEPTUAL MAPPING (MDS)
# ========================================
# Menggunakan data_scaled_original untuk MDS agar tidak terpengaruh KMO variable removal
cat("\nMelakukan Perceptual Mapping (Metric MDS) pada data_scaled_original...\n")

## 
## Melakukan Perceptual Mapping (Metric MDS) pada data_scaled_original...

mds_stress <- NA
mds_coords <- data.frame(Dim1=numeric(0), Dim2=numeric(0))
mds_df <- data.frame()

if (nrow(data_scaled_original) > 1 && ncol(data_scaled_original) > 0) {
  cat("Mengecek duplikat dalam data_scaled_original sebelum MDS...\n")
  duplicated_rows_mds <- duplicated(data_scaled_original)
  data_scaled_mds_unique <- data_scaled_original[!duplicated_rows_mds, ]
  nama_bahan_mds_unique <- nama_bahan[!duplicated_rows_mds] 
  
  # Sesuaikan juga clusters_hc dan clusters_db untuk plot MDS
  clusters_hc_mds_unique <- if(length(clusters_hc) == length(duplicated_rows_mds)) clusters_hc[!duplicated_rows_mds] else rep(NA, nrow(data_scaled_mds_unique))
  clusters_db_mds_unique <- if(length(clusters_db) == length(duplicated_rows_mds)) clusters_db[!duplicated_rows_mds] else rep(NA, nrow(data_scaled_mds_unique))
  
  if (nrow(data_scaled_mds_unique) < 2) {
    cat("Setelah menghapus duplikat, data untuk MDS < 2. MDS tidak dilakukan.\n")
  } else {
    dist_mat_mds_unique <- dist(data_scaled_mds_unique, method = "euclidean")
    
    if (any(dist_mat_mds_unique <= 1e-9)) {
      cat("PERINGATAN MDS: Jarak sangat kecil ditemukan. Mencoba cmdscale...\n")
      mds_result_cmd <- tryCatch({ cmdscale(dist_mat_mds_unique, k = 2, eig = TRUE) }, error = function(e) { cat("cmdscale juga gagal:", e$message, "\n"); return(NULL) })
      if (!is.null(mds_result_cmd) && !is.null(mds_result_cmd$points)) {
        mds_coords <- as.data.frame(mds_result_cmd$points); colnames(mds_coords) <- c("Dim1", "Dim2"); mds_stress <- NA; cat("cmdscale berhasil (Stress N/A).\n")
      } else { cat("cmdscale gagal. MDS tidak dilanjutkan.\n") }
    } else {
      mds_result_iso <- tryCatch({ MASS::isoMDS(dist_mat_mds_unique, k = 2, trace = FALSE) },
                                 error = function(e) { cat("isoMDS gagal:", e$message, ". Mencoba cmdscale...\n");
                                   tryCatch({ cmdscale(dist_mat_mds_unique, k = 2, eig = TRUE) }, error = function(e_cmd) { cat("cmdscale (fallback) juga gagal:", e_cmd$message, "\n"); return(NULL) })
                                 })
      if (!is.null(mds_result_iso) && !is.null(mds_result_iso$points)) {
        if("stress" %in% names(mds_result_iso)){ mds_stress <- mds_result_iso$stress; mds_coords <- as.data.frame(mds_result_iso$points); colnames(mds_coords) <- c("Dim1", "Dim2"); cat("isoMDS berhasil. Stress:", mds_stress, "\n")
        } else if ("points" %in% names(mds_result_iso)) { mds_coords <- as.data.frame(mds_result_iso$points); colnames(mds_coords) <- c("Dim1", "Dim2"); mds_stress <- NA; cat("isoMDS gagal, cmdscale (fallback) berhasil (Stress N/A).\n")
        } else { cat("Struktur hasil MDS tidak dikenali.\n") }
      } else { cat("Baik isoMDS maupun cmdscale (fallback) gagal.\n") }
    }
    if (nrow(mds_coords) > 0) {
      mds_df <- mds_coords %>%
        mutate(NAMA_BAHAN = nama_bahan_mds_unique,
               Cluster_HC = as.factor(clusters_hc_mds_unique), # Pastikan ini faktor
               Cluster_DBSCAN = as.factor(clusters_db_mds_unique)) # Pastikan ini faktor
    }
  }
} else { cat("MDS tidak dilakukan (data_scaled_original tidak valid).\n")}

## Mengecek duplikat dalam data_scaled_original sebelum MDS...
## isoMDS berhasil. Stress: 9.99271

if (nrow(mds_df) > 0) {
  cat("Membuat plot MDS...\n")
  # Plot MDS HC
  if (any(!is.na(mds_df$Cluster_HC))) { # Cek apakah ada cluster HC yang valid
    plot_mds_hc <- ggplot(mds_df, aes(x = Dim1, y = Dim2, color = Cluster_HC)) +
      geom_point(size = 2, alpha = 0.7) +
      ggrepel::geom_text_repel(aes(label = NAMA_BAHAN), size = 2.5, max.overlaps = 15, min.segment.length = 0) +
      scale_color_brewer(palette = "jco", na.value="grey80") + # Tambah na.value
      theme_minimal() +
      labs(title = paste("MDS with HC Clusters (k=", if(!is.na(k_optimal_hc)) k_optimal_hc else "NA", ")"),
           subtitle = paste("Stress:", if(!is.na(mds_stress)) round(mds_stress, 3) else "N/A"),
           x = "MDS Dim 1", y = "MDS Dim 2", color = "HC Cluster")
    print(plot_mds_hc)
    
    # Menyimpan MDS Plot HC
    ggsave("mds_plot_hc.png", plot_mds_hc, width = 10, height = 8)
    cat("MDS Plot HC disimpan sebagai 'mds_plot_hc.png'.\n")
    
  } else {cat("Tidak dapat membuat plot MDS untuk HC (tidak ada cluster HC valid di mds_df).\n")}
  # Plot MDS DBSCAN
  if (any(!is.na(mds_df$Cluster_DBSCAN))) {
    unique_db_mds_plot <- sort(unique(mds_df$Cluster_DBSCAN[!is.na(mds_df$Cluster_DBSCAN)]))
    num_actual_db_mds_plot <- length(unique_db_mds_plot[unique_db_mds_plot != "0"])
    palette_db_mds_values_plot <- RColorBrewer::brewer.pal(n = max(3, num_actual_db_mds_plot +1 ), name = "Dark2")
    final_palette_db_mds_plot <- setNames(rep("grey80", length(unique_db_mds_plot)), as.character(unique_db_mds_plot)) # default grey
    idx_color_mds_plot = 1
    if ("0" %in% names(final_palette_db_mds_plot)) { final_palette_db_mds_plot["0"] <- "grey50" }
    non_noise_mds_names_plot <- as.character(unique_db_mds_plot[unique_db_mds_plot != "0"])
    if (length(non_noise_mds_names_plot) > 0) {
      for(cl_name_mds_plot in non_noise_mds_names_plot){
        final_palette_db_mds_plot[cl_name_mds_plot] <- palette_db_mds_values_plot[idx_color_mds_plot]
        idx_color_mds_plot <- idx_color_mds_plot + 1
        if(idx_color_mds_plot > length(palette_db_mds_values_plot)) idx_color_mds_plot = 1
      }
    }
    plot_mds_dbscan <- ggplot(mds_df, aes(x = Dim1, y = Dim2, color = Cluster_DBSCAN)) +
      geom_point(size = 2, alpha = 0.7) +
      ggrepel::geom_text_repel(aes(label = NAMA_BAHAN), size = 2.5, max.overlaps = 15, min.segment.length = 0) +
      scale_color_manual(values = final_palette_db_mds_plot, na.value="grey80") +
      theme_minimal() +
      labs(title = paste("MDS with DBSCAN (eps=", if(!is.na(eps_chosen)) round(eps_chosen,2) else "NA", ", MinPts=", if(!is.na(min_pts_dbscan_val)) min_pts_dbscan_val else "NA", ")"),
           subtitle = paste("Stress:", if(!is.na(mds_stress)) round(mds_stress, 3) else "N/A"),
           x = "MDS Dim 1", y = "MDS Dim 2", color = "DBSCAN Cluster (0=Noise)")
    print(plot_mds_dbscan)
    
    # Menyimpan MDS Plot DBSCAN
    ggsave("mds_plot_dbscan.png", plot_mds_dbscan, width = 10, height = 8)
    cat("MDS Plot DBSCAN disimpan sebagai 'mds_plot_dbscan.png'.\n")
    
  } else {cat("Tidak dapat membuat plot MDS untuk DBSCAN (tidak ada cluster DBSCAN valid di mds_df).\n")}
} else { cat("Tidak ada data MDS yang valid untuk diplot.\n")}

## Membuat plot MDS...

## Warning: Unknown palette: "jco"

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

## Warning: ggrepel: 1137 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

## Warning: ggrepel: 1127 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## MDS Plot HC disimpan sebagai 'mds_plot_hc.png'.

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

## Warning: ggrepel: 1138 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

## Warning: ggrepel: 1128 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## MDS Plot DBSCAN disimpan sebagai 'mds_plot_dbscan.png'.

# ========================================
# 7. INTERPRETASI DAN INSIGHT
# ========================================
cat("\nMenyiapkan Tabel Karakteristik Cluster...\n")

## 
## Menyiapkan Tabel Karakteristik Cluster...

# `data_imputed` sekarang sudah disesuaikan jika ada variabel yang hilang karena KMO.
# Pastikan panjang `nama_bahan` masih sesuai dengan `data_imputed` (jumlah barisnya).
# Jika KMO menghilangkan variabel (kolom), jumlah baris `data_imputed` seharusnya tidak berubah.
# `clusters_xx` panjangnya harus sama dengan `data_imputed` dan `nama_bahan`.
# Panjang `clusters_xx` akan sama dengan `data_for_clustering` yang baru.
# Jika PCA gagal, `clusters_xx` diinisialisasi dengan panjang `nama_bahan`.

# Cek panjang sebelum membuat data_final
if (nrow(data_imputed) == length(nama_bahan) &&
    length(clusters_hc) == length(nama_bahan) &&
    length(clusters_fcm) == length(nama_bahan) &&
    length(clusters_db) == length(nama_bahan)
) {
  data_final <- data_imputed %>%
    mutate(NAMA_BAHAN = nama_bahan,
           Cluster_HC = as.factor(clusters_hc),
           Cluster_FCM = as.factor(clusters_fcm),
           Cluster_DBSCAN = as.factor(clusters_db)) %>%
    dplyr::select(NAMA_BAHAN, Cluster_HC, Cluster_FCM, Cluster_DBSCAN, everything())
  
  # Menyimpan data_final ke CSV
  write.csv(data_final, "hasil_clustering_makanan1.csv", row.names = FALSE)
  cat("Hasil clustering makanan disimpan sebagai 'hasil_clustering_makanan1.csv'.\n")
  
} else {
  cat("PERINGATAN: Ada ketidakcocokan panjang antara data_imputed, nama_bahan, atau hasil cluster.\n")
  cat("Dim data_imputed:", dim(data_imputed), "\n")
  cat("Length nama_bahan:", length(nama_bahan), "\n")
  cat("Length clusters_hc:", length(clusters_hc), "\n")
  cat("Membuat data_final dengan potensi data yang tidak sinkron atau NA.\n")
  # Upaya membuat data_final dengan baris minimum yang cocok
  min_rows_final <- min(nrow(data_imputed), length(nama_bahan), length(clusters_hc), length(clusters_fcm), length(clusters_db))
  if (min_rows_final > 0) {
    data_final <- data_imputed[1:min_rows_final, , drop=FALSE] %>%
      mutate(NAMA_BAHAN = nama_bahan[1:min_rows_final],
             Cluster_HC = as.factor(clusters_hc[1:min_rows_final]),
             Cluster_FCM = as.factor(clusters_fcm[1:min_rows_final]),
             Cluster_DBSCAN = as.factor(clusters_db[1:min_rows_final])) %>%
      dplyr::select(NAMA_BAHAN, Cluster_HC, Cluster_FCM, Cluster_DBSCAN, everything())
    
    # Menyimpan data_final yang disesuaikan ke CSV
    write.csv(data_final, "hasil_clustering_makanan1.csv", row.names = FALSE)
    cat("Hasil clustering makanan (disesuaikan) disimpan sebagai 'hasil_clustering_makanan.csv'.\n")
    
  } else {
    data_final <- data.frame() # Kosong jika tidak bisa
    cat("Tidak ada data_final yang dapat dibuat atau disimpan.\n")
  }
}

## Hasil clustering makanan disimpan sebagai 'hasil_clustering_makanan1.csv'.

# Fungsi create_summary_table (Anda yang sudah diperbaiki untuk syntax across)
create_summary_table <- function(df_final, cluster_col_name, df_imputed_source) {
  imputed_nutrient_cols_summary <- colnames(df_imputed_source)
  if (!cluster_col_name %in% names(df_final)) {
    message(paste("Kolom cluster", cluster_col_name, "tidak ada di df_final."))
    return(NULL)
  }
  if (all(is.na(df_final[[cluster_col_name]]))) {
    message(paste("Semua nilai di kolom cluster", cluster_col_name, "adalah NA."))
    return(NULL)
  }
  df_final[[cluster_col_name]] <- as.factor(df_final[[cluster_col_name]])
  # --- Perbaikan: Gunakan kept_vars_after_missing untuk memastikan kolom numerik yang benar ---
  # Pastikan kept_vars_after_missing tersedia di global environment atau diteruskan sebagai argumen
  # Untuk kasus ini, kita asumsikan kept_vars_after_missing tersedia dari global scope
  numeric_cols_summary <- intersect(kept_vars_after_missing, names(df_final))
  # --- Akhir Perbaikan ---
  
  if(length(numeric_cols_summary) == 0){
    message("Tidak ada kolom numerik cocok untuk summary.")
    return(NULL)
  }
  summary_t <- df_final %>%
    filter(!is.na(.data[[cluster_col_name]])) %>%
    group_by(!!sym(cluster_col_name)) %>%
    summarise(
      across(all_of(numeric_cols_summary), ~mean(., na.rm = TRUE)),
      .groups = 'drop'
    ) %>%
    left_join(
      df_final %>%
        filter(!is.na(.data[[cluster_col_name]])) %>%
        count(!!sym(cluster_col_name), name = "Jumlah_Anggota"),
      by = cluster_col_name
    )
  return(summary_t)
}


if(nrow(data_final) > 0) {
  cat("\n--- Tabel Karakteristik Cluster HC (Rata-rata Variabel Asli) ---\n")
  summary_hc <- create_summary_table(data_final, "Cluster_HC", data_imputed) # data_imputed sudah disesuaikan
  if(!is.null(summary_hc)) {
    print(as.data.frame(summary_hc))
    write.csv(summary_hc, "summary_table_hc.csv", row.names = FALSE)
    cat("Tabel ringkasan HC disimpan sebagai 'summary_table_hc.csv'.\n")
  }
  
  cat("\n--- Tabel Karakteristik Cluster FCM (Rata-rata Variabel Asli) ---\n")
  summary_fcm <- create_summary_table(data_final, "Cluster_FCM", data_imputed)
  if(!is.null(summary_fcm)) {
    print(as.data.frame(summary_fcm))
    write.csv(summary_fcm, "summary_table_fcm.csv", row.names = FALSE)
    cat("Tabel ringkasan FCM disimpan sebagai 'summary_table_fcm.csv'.\n")
  }
  
  cat("\n--- Tabel Karakteristik Cluster DBSCAN (Rata-rata Variabel Asli, Termasuk Noise=0) ---\n")
  summary_dbscan <- create_summary_table(data_final, "Cluster_DBSCAN", data_imputed)
  if(!is.null(summary_dbscan)) {
    print(as.data.frame(summary_dbscan))
    write.csv(summary_dbscan, "summary_table_dbscan.csv", row.names = FALSE)
    cat("Tabel ringkasan DBSCAN disimpan sebagai 'summary_table_dbscan.csv'.\n")
  }
  
  cat("\nContoh Anggota Cluster 1 HC:\n")
  if("Cluster_HC" %in% names(data_final) && any(data_final$Cluster_HC == 1, na.rm = TRUE)){
    cluster_1_hc_examples <- data_final %>% filter(Cluster_HC == 1) %>% dplyr::select(NAMA_BAHAN)
    print(head(cluster_1_hc_examples))
    write.csv(cluster_1_hc_examples, "cluster_1_hc_examples.csv", row.names = FALSE)
    cat("Contoh anggota Cluster 1 HC disimpan sebagai 'cluster_1_hc_examples.csv'.\n")
  } else {
    cat("Tidak dapat menampilkan contoh anggota cluster HC 1.\n")
  }
} else {
  cat("data_final kosong, tabel summary tidak dapat dibuat.\n")
}

## 
## --- Tabel Karakteristik Cluster HC (Rata-rata Variabel Asli) ---
##   Cluster_HC    NATRIUM     SERAT   THIAMIN      BESI    KALSIUM    FOSFOR
## 1          1   86.04751 2.0407670 0.1616761  2.128125   98.17756  99.47727
## 2          2   65.86122 4.9436735 0.2220000  3.787347  117.51429 172.96327
## 3          3  404.83571 0.9878571 0.3230714  3.722143  111.04286 229.20000
## 4          4  290.78431 1.2431373 0.3713725 16.056863 1112.50980 818.62745
## 5          5 4972.14286 0.9000000 0.0700000  3.100000  124.71429 231.00000
##       LEMAK   PROTEIN       ABU      AIR    ENERGI        KH Jumlah_Anggota
## 1  2.247869  5.818324  1.359517 77.34020  95.08381 13.193750            704
## 2  9.202041  8.583265  2.275510 14.92122 374.11020 64.285306            245
## 3 28.222857 21.143571  3.866429 37.11643 370.87857  8.483571            140
## 4 11.772549 34.558824 10.566667 28.01961 312.96078 18.376471             51
## 5  2.457143 24.185714 33.742857 46.65714 172.71429 12.885714              7
## Tabel ringkasan HC disimpan sebagai 'summary_table_hc.csv'.
## 
## --- Tabel Karakteristik Cluster FCM (Rata-rata Variabel Asli) ---
##   Cluster_FCM   NATRIUM    SERAT   THIAMIN     BESI   KALSIUM    FOSFOR
## 1           1 678.17333 0.848000 0.3958667 5.204000 196.46667 270.93333
## 2           2  56.52767 2.524000 0.1220889 1.729778  77.54444  60.74444
## 3           3 177.37037 1.406667 0.2348889 2.789259 127.98519 164.79630
## 4           4  73.69091 4.204091 0.1826818 3.264545  93.27727 137.19545
## 5           5 322.90152 2.665152 0.3381061 8.695455 502.71212 496.43939
##       LEMAK   PROTEIN      ABU      AIR    ENERGI        KH Jumlah_Anggota
## 1 14.313333 22.930667 5.057333 47.68933 258.76000  7.630667             75
## 2  1.144667  2.457333 1.069333 82.46444  69.19333 12.734000            450
## 3  5.256296 12.237407 2.080370 69.01630 142.15926 11.392222            270
## 4  7.505455  5.929545 1.801818 16.95136 361.68182 67.345000            220
## 5 28.384848 27.702273 7.669697 15.19470 448.42424 23.174242            132
## Tabel ringkasan FCM disimpan sebagai 'summary_table_fcm.csv'.
## 
## --- Tabel Karakteristik Cluster DBSCAN (Rata-rata Variabel Asli, Termasuk Noise=0) ---
##   Cluster_DBSCAN   NATRIUM    SERAT   THIAMIN     BESI   KALSIUM    FOSFOR
## 1              0 315.18824 3.650588 0.3310824 5.940235 278.28471 294.82118
## 2              1  67.96461 1.806648 0.1276454 1.749584  73.12188  86.65374
##       LEMAK   PROTEIN      ABU      AIR   ENERGI       KH Jumlah_Anggota
## 1 14.720000 17.535529 4.744941 36.53224 307.8212 26.41129            425
## 2  2.977562  5.039197 1.127978 68.59958 134.1620 22.20028            722
## Tabel ringkasan DBSCAN disimpan sebagai 'summary_table_dbscan.csv'.
## 
## Contoh Anggota Cluster 1 HC:
##                    NAMA_BAHAN
## 1                        <NA>
## 2 Jagung muda, kuning, mentah
## 3                        Nasi
## 4                    Nasi tim
## 5                Beras, tapai
## 6           Beras merah, Nasi
## Contoh anggota Cluster 1 HC disimpan sebagai 'cluster_1_hc_examples.csv'.

# Menyimpan log hasil penting ke file teks
sink("analysis_log.txt")
cat("--- Ringkasan Analisis Profil Gizi Makanan Remaja ---\n\n")
cat("Dimensi data setelah load:", dim(data_numeric), "\n\n")
cat("Persentase Missing Value per Variabel (Hanya Nutrisi):\n")
print(missing_summary)
cat("\nVariabel yang digunakan setelah filtering missing value (>25%):", paste(kept_vars_after_missing, collapse=", "), "\n")
cat("Dimensi data setelah filtering variabel:", dim(data_filtered), "\n\n")
cat("Jumlah NA setelah imputasi:", sum(is.na(data_imputed)), "\n\n")
cat("Dimensi data setelah scaling awal:", dim(data_scaled_original), "\n\n")

cat("--- Hasil Clustering --- \n")
cat("Jumlah cluster optimal (k) yang dipilih untuk HC:", k_optimal_hc, "\n")
cat("Rata-rata Silhouette Score HC:", mean_sil_hc, "\n")
cat("Rata-rata Silhouette Score FCM:", mean_sil_fcm, "\n")
cat("Rata-rata Silhouette Score DBSCAN (non-noise):", mean_sil_dbscan, "\n")
cat("DBSCAN - Jumlah Noise Points:", num_noise_dbscan, "\n")
cat("DBSCAN - Jumlah Cluster Ditemukan (tidak termasuk noise):", num_clusters_dbscan, "\n\n")

# Bagian MDS tetap karena bisa dilakukan tanpa PCA
cat("--- Hasil Perceptual Mapping (MDS) ---\n")
if(exists("mds_stress") && !is.na(mds_stress)) {
  cat("MDS Stress:", mds_stress, "\n\n")
} else {
  cat("MDS Stress: N/A (MDS mungkin tidak dijalankan atau gagal)\n\n")
}

sink()
cat("Log analisis lengkap disimpan sebagai 'analysis_log.txt'.\n")

## Log analisis lengkap disimpan sebagai 'analysis_log.txt'.

# ========================================
# AKHIR SCRIPT
# ========================================
cat("\nAnalisis Selesai.\n")

## 
## Analisis Selesai.