# ── Muat Dataset Utama ───────────────────────────────────────────────────────
library(readr)
library(tidyverse)

data_stunting_raw <- read_csv("C:/Users/User/Downloads/data_stunting_raw.csv")
str(data_stunting_raw)

# ── Distribusi Frekuensi & Persentase ────────────────────────────────────────
library(knitr)
library(kableExtra)

tabel_target <- data_stunting_raw %>%
  count(status_stunting) %>%
  mutate(persentase = round(n / sum(n) * 100, 2))

kable(tabel_target,
      caption = "Distribusi Frekuensi Variabel Target: status_stunting",
      col.names = c("Status Stunting", "Frekuensi (n)", "Persentase (%)"),
      align = "lrr") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(tabel_target$status_stunting == "Stunting"),
           background = "#FFEBEE", color = "#B71C1C", bold = TRUE)

# ── Rasio Imbalance ──────────────────────────────────────────────────────────
pct_minoritas <- tabel_target %>%
  filter(status_stunting == "Stunting") %>% pull(persentase)

cat("Persentase kelas minoritas (Stunting) :", pct_minoritas, "%\n")

cat("Persentase kelas mayoritas (Normal)   :", 100 - pct_minoritas, "%\n")

cat("Rasio imbalance (Normal : Stunting)   :",
    round((100 - pct_minoritas) / pct_minoritas, 2), ": 1\n")

# ── Bar Chart Distribusi Target ───────────────────────────────────────────────
ggplot(data_stunting_raw, aes(x = status_stunting, fill = status_stunting)) +
  geom_bar(width = 0.55, colour = "white", linewidth = 0.4) +
  geom_text(stat = "count",
            aes(label = paste0(..count.., "\n(",
                round(..count.. / nrow(data_stunting_raw) * 100, 1), "%)")),
            vjust = -0.4, size = 4.2, fontface = "bold") +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING),
                    guide = guide_legend(title.position = "top")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title    = "Distribusi Variabel Target: Status Stunting",
       subtitle = "Dataset data_stunting_raw.csv — 500 Balita, 50 Kab/Kota",
       x        = "Status Stunting",
       y        = "Jumlah Balita",
       fill     = "Status",
       caption  = "Sumber: data_stunting_raw.csv") +
  annotate("text", x = 1.5, y = 430,
           label = paste0("Rasio Imbalance\n",
                          round((100-pct_minoritas)/pct_minoritas, 1), " : 1"),
           colour = COL_ACCENT, size = 3.5, fontface = "italic")

# ── Struktur Data ─────────────────────────────────────────────────────────────
library(corrplot)

str(data_stunting_raw)

head(data_stunting_raw)

# ── Matriks Korelasi Variabel Numerik ─────────────────────────────────────────
numerik_vars <- data_stunting_raw %>%
  select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
         pendapatan, diare_3bln, kunjungan_posyandu, zscore_tbu)

cor_matrix <- cor(numerik_vars, use = "complete.obs")

corrplot(cor_matrix,
         method      = "color",
         type        = "upper",
         order       = "hclust",
         addCoef.col = "#37474F",
         tl.col      = "#1A237E",
         tl.srt      = 45,
         tl.cex      = 0.85,
         number.cex  = 0.72,
         col         = colorRampPalette(c(COL_STUNTING, "white", COL_NORMAL))(200),
         cl.cex      = 0.75,
         title       = "Matriks Korelasi Variabel Numerik",
         mar         = c(0, 0, 2, 0))

# ── Pasangan Korelasi Tertinggi ───────────────────────────────────────────────
cor_df <- as.data.frame(as.table(cor_matrix)) %>%
  filter(Var1 != Var2) %>%
  mutate(abs_cor = abs(Freq)) %>%
  arrange(desc(abs_cor)) %>%
  filter(!duplicated(abs_cor))

kable(head(cor_df, 5),
      caption = "5 Pasangan Variabel dengan Korelasi Tertinggi",
      col.names = c("Variabel 1", "Variabel 2", "Korelasi", "|Korelasi|"),
      digits = 3, align = "llrr") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  column_spec(3:4, bold = TRUE)

# ── Statistik Deskriptif Lengkap ──────────────────────────────────────────────
library(e1071)

stat_desk <- data_stunting_raw %>%
  select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
         pendapatan, diare_3bln, kunjungan_posyandu, zscore_tbu) %>%
  summarise(across(everything(), list(
    Min    = ~min(., na.rm = TRUE),
    Maks   = ~max(., na.rm = TRUE),
    Mean   = ~mean(., na.rm = TRUE),
    Median = ~median(., na.rm = TRUE),
    SD     = ~sd(., na.rm = TRUE),
    P25    = ~quantile(., 0.25, na.rm = TRUE),
    P75    = ~quantile(., 0.75, na.rm = TRUE),
    Skew   = ~skewness(., na.rm = TRUE)
  ), .names = "{.col}__{.fn}")) %>%
  pivot_longer(everything(),
               names_to = c("Variabel", "Statistik"),
               names_sep = "__") %>%
  pivot_wider(names_from = Statistik, values_from = value)

kable(stat_desk,
      caption = "Statistik Deskriptif Variabel Numerik",
      digits  = 2,
      align   = c("l", rep("r", 8))) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  column_spec(1, bold = TRUE, background = "#E8EAF6") %>%
  add_header_above(c(" " = 1, "Ukuran Pemusatan & Sebaran" = 7, "Bentuk" = 1))

# ── Visualisasi 1: Histogram Z-Score TB/U ─────────────────────────────────────
ggplot(data_stunting_raw, aes(x = zscore_tbu, fill = status_stunting)) +
  geom_histogram(bins = 40, alpha = 0.78, position = "identity",
                 colour = "white", linewidth = 0.2) +
  geom_vline(xintercept = -2, linetype = "dashed",
             colour = "#E65100", linewidth = 1) +
  annotate("rect", xmin = -Inf, xmax = -2,
           ymin = -Inf, ymax = Inf,
           fill = "#FFCCBC", alpha = 0.15) +
  annotate("label", x = -3.2, y = 30,
           label = "Zona\nStunting\n(< -2 SD)",
           fill = "#FFF3E0", colour = "#E65100",
           size = 3, label.size = 0.3) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_x_continuous(breaks = seq(-6, 4, 1)) +
  labs(title    = "Distribusi Z-Score TB/U berdasarkan Status Stunting",
       subtitle = "Garis putus-putus: cut-off WHO −2 SD",
       x        = "Z-Score TB/U",
       y        = "Frekuensi",
       fill     = "Status Stunting",
       caption  = "Sumber: data_stunting_raw.csv")

# ── Visualisasi 2: Boxplot Tinggi & Berat Badan ───────────────────────────────
library(patchwork)

p1 <- ggplot(data_stunting_raw,
             aes(x = status_stunting, y = tb_sekarang,
                 fill = status_stunting)) +
  geom_boxplot(outlier.colour = "#FF6F00", outlier.shape = 21,
               outlier.size = 2.5, outlier.fill = "#FFB300",
               notch = FALSE, width = 0.5, linewidth = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3, fill = "white", colour = COL_ACCENT) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  labs(title    = "Tinggi Badan (cm)",
       x        = "Status",
       y        = "Tinggi Badan (cm)") +
  theme(legend.position = "none")

p2 <- ggplot(data_stunting_raw,
             aes(x = status_stunting, y = bb_sekarang,
                 fill = status_stunting)) +
  geom_boxplot(outlier.colour = "#FF6F00", outlier.shape = 21,
               outlier.size = 2.5, outlier.fill = "#FFB300",
               notch = FALSE, width = 0.5, linewidth = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3, fill = "white", colour = COL_ACCENT) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  labs(title    = "Berat Badan (kg)",
       x        = "Status",
       y        = "Berat Badan (kg)") +
  theme(legend.position = "none")

(p1 | p2) +
  plot_annotation(
    title    = "Distribusi Antropometri per Status Stunting",
    subtitle = "Berlian putih = nilai rata-rata; titik oranye = outlier",
    theme    = theme(plot.title    = element_text(face = "bold", colour = "#1A237E"),
                     plot.subtitle = element_text(colour = "#546E7A"))
  )

# ── Visualisasi 3: Scatter Plot Z-Score vs Tinggi Badan ───────────────────────
ggplot(data_stunting_raw,
       aes(x = tb_sekarang, y = zscore_tbu, colour = status_stunting)) +
  geom_point(alpha = 0.55, size = 2, shape = 16) +
  geom_smooth(aes(group = status_stunting),
              method = "lm", se = TRUE, linewidth = 0.9, linetype = "solid") +
  geom_hline(yintercept = -2, linetype = "dashed",
             colour = "#E65100", linewidth = 0.9) +
  annotate("label", x = 110, y = -1.7,
           label = "Cut-off −2 SD", fill = "#FFF3E0",
           colour = "#E65100", size = 3, label.size = 0.3) +
  scale_colour_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  labs(title    = "Scatter Plot: Tinggi Badan vs Z-Score TB/U",
       subtitle = "Garis regresi per kelompok status stunting",
       x        = "Tinggi Badan Sekarang (cm)",
       y        = "Z-Score TB/U",
       colour   = "Status Stunting",
       caption  = "Sumber: data_stunting_raw.csv")

# ── Nilai Skewness ─────────────────────────────────────────────────────────────
skew_df <- data_stunting_raw %>%
  select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
         pendapatan, diare_3bln, kunjungan_posyandu, zscore_tbu) %>%
  summarise(across(everything(), ~round(skewness(., na.rm = TRUE), 3))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Skewness") %>%
  mutate(Kategori = case_when(
    abs(Skewness) < 0.5  ~ "Simetris",
    abs(Skewness) < 1    ~ "Skewed Sedang",
    TRUE                 ~ "Skewed Berat"
  ))

kable(skew_df,
      caption = "Nilai Skewness Variabel Numerik",
      align   = "lrc") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  column_spec(3, color = ifelse(abs(skew_df$Skewness) >= 1, COL_STUNTING, "black"),
              bold = ifelse(abs(skew_df$Skewness) >= 1, TRUE, FALSE))

# ── Tabel Proporsi per Variabel Kategorik ─────────────────────────────────────
buat_tabel_proporsi <- function(data, var_prediktor) {
  data %>%
    group_by(across(all_of(c(var_prediktor, "status_stunting")))) %>%
    summarise(n = n(), .groups = "drop") %>%
    group_by(across(all_of(var_prediktor))) %>%
    mutate(proporsi = round(n / sum(n) * 100, 1)) %>%
    ungroup()
}

var_kategorik <- c("pend_ibu", "akses_air", "sanitasi",
                   "asi_eksklusif", "imunisasi")

for (v in var_kategorik) {
  cat("\n=== Tabel Proporsi:", v, "vs status_stunting ===\n")
  print(buat_tabel_proporsi(data_stunting_raw, v))
}

# ── Visualisasi 1: Pendidikan Ibu ─────────────────────────────────────────────
data_stunting_raw %>%
  filter(!is.na(pend_ibu)) %>%
  mutate(pend_ibu = factor(pend_ibu,
         levels = c("Tidak Sekolah", "SD", "SMP", "SMA", "PT"))) %>%
  ggplot(aes(x = pend_ibu, fill = status_stunting)) +
  geom_bar(position = "fill", colour = "white", linewidth = 0.3) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title    = "Proporsi Status Stunting per Pendidikan Ibu",
       subtitle = "Semakin rendah pendidikan ibu → proporsi stunting cenderung lebih tinggi",
       x        = "Tingkat Pendidikan Ibu",
       y        = "Proporsi (%)",
       fill     = "Status Stunting",
       caption  = "Sumber: data_stunting_raw.csv") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

# ── Visualisasi 2: Panel 4 Variabel Kategorik ─────────────────────────────────
data_stunting_raw %>%
  select(akses_air, sanitasi, asi_eksklusif, imunisasi, status_stunting) %>%
  pivot_longer(-status_stunting, names_to = "variabel", values_to = "kategori") %>%
  mutate(variabel = recode(variabel,
    akses_air    = "Akses Air Bersih",
    sanitasi     = "Sanitasi",
    asi_eksklusif = "ASI Eksklusif",
    imunisasi    = "Imunisasi")) %>%
  ggplot(aes(x = kategori, fill = status_stunting)) +
  geom_bar(position = "fill", colour = "white", linewidth = 0.3) +
  facet_wrap(~variabel, scales = "free_x", nrow = 2) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title    = "Proporsi Status Stunting per Variabel Kategorik",
       subtitle = "Faktor akses layanan dan perilaku kesehatan",
       x        = NULL,
       y        = "Proporsi (%)",
       fill     = "Status Stunting",
       caption  = "Sumber: data_stunting_raw.csv") +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

# ── IK 07a: Identifikasi Outlier Berdasarkan Batas WHO ───────────────────────
batas_who <- data.frame(
  Variabel = c("tb_sekarang", "bb_sekarang", "bb_lahir"),
  Batas_Bawah = c(55, 4, 1.5),
  Batas_Atas  = c(120, 25, 5.0)
)

kable(batas_who,
      caption = "Batas Wajar WHO untuk Balita Usia 6–59 Bulan",
      col.names = c("Variabel", "Batas Bawah", "Batas Atas"),
      align = "lrr") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  column_spec(1, bold = TRUE)

# ── Identifikasi outlier per variabel ────────────────────────────────────────
cat("\n=== tb_sekarang: Nilai di luar batas WHO (55-120 cm) ===\n")

tb_outlier <- data_stunting_raw %>%
  filter(tb_sekarang < 55 | tb_sekarang > 120) %>%
  select(id_balita, usia_bulan, tb_sekarang, status_stunting)
print(tb_outlier)

cat("Jumlah outlier tb_sekarang:", nrow(tb_outlier), "\n")

cat("\n=== bb_sekarang: Nilai di luar batas WHO (4-25 kg) ===\n")

bb_outlier <- data_stunting_raw %>%
  filter(bb_sekarang < 4 | bb_sekarang > 25) %>%
  select(id_balita, usia_bulan, bb_sekarang, status_stunting)
print(bb_outlier)

cat("Jumlah outlier bb_sekarang:", nrow(bb_outlier), "\n")

cat("\n=== bb_lahir: Nilai di luar batas WHO (1.5-5.0 kg) ===\n")

bblahir_outlier <- data_stunting_raw %>%
  filter(bb_lahir < 1.5 | bb_lahir > 5.0) %>%
  select(id_balita, usia_bulan, bb_lahir, status_stunting)
print(bblahir_outlier)

cat("Jumlah outlier bb_lahir:", nrow(bblahir_outlier), "\n")

cat("\n=== RINGKASAN OUTLIER ===\n")

cat("tb_sekarang:", nrow(tb_outlier), "outlier dari 500 observasi (",
    round(nrow(tb_outlier)/500*100,1), "%)\n")

cat("bb_sekarang:", nrow(bb_outlier), "outlier dari 500 observasi (",
    round(nrow(bb_outlier)/500*100,1), "%)\n")

cat("bb_lahir   :", nrow(bblahir_outlier), "outlier dari 500 observasi (",
    round(nrow(bblahir_outlier)/500*100,1), "%)\n")

# Cek batas soal
cat("\n=== Batas outlier SOAL IK09 untuk perbandingan ===\n")

cat("tb_sekarang soal: 30-130 cm\n")

cat("bb_sekarang soal: 2-20 kg\n")

cat("bb_lahir    soal: 0.5-5.5 kg\n")

tb_soal  <- data_stunting_raw %>% filter(tb_sekarang < 30 | tb_sekarang > 130)
bb_soal  <- data_stunting_raw %>% filter(bb_sekarang < 2  | bb_sekarang > 20)
bbl_soal <- data_stunting_raw %>% filter(bb_lahir < 0.5   | bb_lahir > 5.5)

cat("Outlier tb_sekarang (batas soal):", nrow(tb_soal), "\n")

cat("Outlier bb_sekarang (batas soal):", nrow(bb_soal), "\n")

cat("Outlier bb_lahir    (batas soal):", nrow(bbl_soal), "\n")

# Kecukupan data
cat("\n=== IK 07b: Kecukupan Data ===\n")

n_obs       <- nrow(data_stunting_raw)
n_prediktor <- 15
min_obs     <- n_prediktor * 10
rasio       <- n_obs / n_prediktor

cat("Jumlah observasi         :", n_obs, "\n")

cat("Jumlah variabel prediktor:", n_prediktor, "\n")

cat("Minimum obs (10x pred)   :", min_obs, "\n")

cat("Rasio observasi/prediktor:", round(rasio, 2), "\n")

cat("Kecukupan data:", ifelse(n_obs >= min_obs, "CUKUP ✓", "TIDAK CUKUP ✗"), "\n")

cat("\n=== Missing Values per Variabel ===\n")

data_stunting_raw %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Missing") %>%
  filter(Missing > 0) %>%
  mutate(Persen = round(Missing/500*100, 1)) %>%
  print()

# ── IK 07a — Identifikasi Outlier per Variabel (Batas WHO ±3SD) ──────────────
outlier_tb <- data_stunting_raw %>%
  filter(tb_sekarang < 61.2 | tb_sekarang > 124.0) %>%
  select(id_balita, usia_bulan, tb_sekarang, status_stunting)

cat("=== Outlier tb_sekarang (WHO: 61.2–124.0 cm) ===\n")

cat("Jumlah outlier:", nrow(outlier_tb), "\n")

print(outlier_tb)

outlier_bb <- data_stunting_raw %>%
  filter(bb_sekarang < 5.5 | bb_sekarang > 24.0) %>%
  select(id_balita, usia_bulan, bb_sekarang, status_stunting)

cat("\n=== Outlier bb_sekarang (WHO: 5.5–24.0 kg) ===\n")

cat("Jumlah outlier:", nrow(outlier_bb), "\n")

print(outlier_bb)

outlier_bbl <- data_stunting_raw %>%
  filter(!is.na(bb_lahir)) %>%
  filter(bb_lahir < 2.0 | bb_lahir > 4.5) %>%
  select(id_balita, usia_bulan, bb_lahir, status_stunting)

cat("\n=== Outlier bb_lahir (WHO: 2.0–4.5 kg) ===\n")

cat("Jumlah outlier:", nrow(outlier_bbl), "\n")

print(outlier_bbl)

outlier_zscore <- data_stunting_raw %>%
  filter(zscore_tbu < -6 | zscore_tbu > 6) %>%
  select(id_balita, usia_bulan, zscore_tbu, status_stunting)

cat("\n=== Outlier zscore_tbu (WHO: -6 s/d +6) ===\n")

cat("Jumlah outlier:", nrow(outlier_zscore), "\n")

print(outlier_zscore)

cat("\n=== Ringkasan Statistik Variabel Antropometri ===\n")

data_stunting_raw %>%
  select(tb_sekarang, bb_sekarang, bb_lahir, zscore_tbu) %>%
  summary()

cat("\n=== REKAP TOTAL NILAI TIDAK VALID (Batas WHO) ===\n")

cat("tb_sekarang  :", nrow(outlier_tb),     "data tidak valid\n")

cat("bb_sekarang  :", nrow(outlier_bb),     "data tidak valid\n")

cat("bb_lahir     :", nrow(outlier_bbl),    "data tidak valid\n")

cat("zscore_tbu   :", nrow(outlier_zscore), "data tidak valid\n")

cat("TOTAL        :", nrow(outlier_tb) + nrow(outlier_bb) +
                      nrow(outlier_bbl) + nrow(outlier_zscore),
    "data tidak valid\n")

# ── IK 07b — Rule of Thumb Kecukupan Data ────────────────────────────────────
n_observasi  <- nrow(data_stunting_raw)
n_prediktor  <- 14
min_obs      <- n_prediktor * 10
rasio        <- round(n_observasi / n_prediktor, 2)

cat("\n=== Penilaian Kecukupan Data (Rule of Thumb: 10 obs/prediktor) ===\n")

cat("Jumlah observasi          :", n_observasi, "\n")

cat("Jumlah variabel prediktor :", n_prediktor, "\n")

cat("Minimum observasi (10x)   :", min_obs, "\n")

cat("Rasio obs/prediktor       :", rasio, "\n")

if (n_observasi >= min_obs) {
  cat("Kesimpulan : CUKUP ✓ (rasio", rasio, ">= 10)\n")
} else {
  cat("Kesimpulan : TIDAK CUKUP ✗ (rasio", rasio, "< 10)\n")
}

cat("\n=== Missing Values per Variabel ===\n")

data_stunting_raw %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(),
               names_to  = "Variabel",
               values_to = "Jumlah_NA") %>%
  filter(Jumlah_NA > 0) %>%
  print()

# ── IK 08a: Ringkasan variabel & keputusan penggunaan ─────────────────────────
cat("=== Dimensi Dataset ===\n")

cat("Jumlah baris (observasi):", nrow(data_stunting_raw), "\n")

cat("Jumlah kolom (variabel) :", ncol(data_stunting_raw), "\n")

cat("\n=== Nama Variabel ===\n")

print(names(data_stunting_raw))

cat("\n=== Tipe Data per Variabel ===\n")

glimpse(data_stunting_raw)

cat("\n=== Distribusi Target: status_stunting ===\n")

data_stunting_raw %>%
  count(status_stunting) %>%
  mutate(persen = round(n/sum(n)*100, 2)) %>%
  print()

cat("\n=== Proporsi Stunting per Jenis Kelamin ===\n")

data_stunting_raw %>%
  group_by(jenis_kelamin, status_stunting) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(jenis_kelamin) %>%
  mutate(proporsi = round(n/sum(n)*100, 2)) %>%
  print()

cat("\n=== Proporsi Stunting per Kelompok Usia ===\n")

data_stunting_raw %>%
  mutate(kelompok_usia = case_when(
    usia_bulan >= 6  & usia_bulan <= 11 ~ "6-11 bln",
    usia_bulan >= 12 & usia_bulan <= 23 ~ "12-23 bln",
    usia_bulan >= 24 & usia_bulan <= 35 ~ "24-35 bln",
    usia_bulan >= 36 & usia_bulan <= 47 ~ "36-47 bln",
    usia_bulan >= 48 & usia_bulan <= 59 ~ "48-59 bln"
  )) %>%
  group_by(kelompok_usia, status_stunting) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(kelompok_usia) %>%
  mutate(proporsi = round(n/sum(n)*100, 2)) %>%
  arrange(kelompok_usia) %>%
  print()

set.seed(42)
cat("\n=== Simulasi Stratified Split 80:20 ===\n")

cat("Training set (80%):\n")

cat("  - Target Normal  : ekspektasi", round(371*0.8), "observasi\n")

cat("  - Target Stunting: ekspektasi", round(129*0.8), "observasi\n")

cat("Test set (20%):\n")

cat("  - Target Normal  : ekspektasi", round(371*0.2), "observasi\n")

cat("  - Target Stunting: ekspektasi", round(129*0.2), "observasi\n")

cat("\n=== Missing Values per Variabel Prediktor ===\n")

data_stunting_raw %>%
  select(-id_balita, -kab_kota, -status_stunting) %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(),
               names_to  = "Variabel",
               values_to = "Jumlah_NA") %>%
  print()

# ── IK 09a: RINGKASAN MASALAH DATA SEBELUM CLEANING ──────────────────────────
cat("=== KONDISI DATA SEBELUM CLEANING ===\n")

cat("Jumlah baris:", nrow(data_stunting_raw), "\n")

cat("\n--- Missing Values ---\n")

data_stunting_raw %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "NA_Count") %>%
  filter(NA_Count > 0) %>%
  mutate(Persen = round(NA_Count / 500 * 100, 1)) %>%
  print()

cat("\n--- Outlier Berdasarkan Batas Soal ---\n")

n_tb  <- sum(data_stunting_raw$tb_sekarang < 30 | data_stunting_raw$tb_sekarang > 130, na.rm = TRUE)
n_bb  <- sum(data_stunting_raw$bb_sekarang < 2  | data_stunting_raw$bb_sekarang > 20,  na.rm = TRUE)
n_bbl <- sum(data_stunting_raw$bb_lahir    < 0.5| data_stunting_raw$bb_lahir    > 5.5,  na.rm = TRUE)
cat("Outlier tb_sekarang  (batas 30-130 cm) :", n_tb,  "\n")

cat("Outlier bb_sekarang  (batas 2-20 kg)   :", n_bb,  "\n")

cat("Outlier bb_lahir     (batas 0.5-5.5 kg):", n_bbl, "\n")

cat("\n--- Data Duplikat ---\n")

n_dup <- sum(duplicated(data_stunting_raw))
cat("Jumlah baris duplikat:", n_dup, "\n")

# ── LANGKAH 1 — Hapus Outlier Ekstrem ────────────────────────────────────────
data_step1 <- data_stunting_raw %>%
  filter(
    tb_sekarang >= 30  & tb_sekarang <= 130,
    bb_sekarang >= 2   & bb_sekarang <= 20,
    (bb_lahir   >= 0.5 & bb_lahir   <= 5.5) | is.na(bb_lahir)
  )

cat("\nSetelah hapus outlier antropometri:", nrow(data_step1), "baris tersisa\n")

cat("Baris yang dihapus:", nrow(data_stunting_raw) - nrow(data_step1), "\n")

# ── LANGKAH 2 — Imputasi Missing Values ──────────────────────────────────────
data_step2 <- data_step1 %>%
  mutate(
    kelompok_usia = case_when(
      usia_bulan <= 12 ~ "0-12",
      usia_bulan <= 24 ~ "13-24",
      usia_bulan <= 36 ~ "25-36",
      usia_bulan <= 48 ~ "37-48",
      TRUE             ~ "49-59"
    )
  )

cat("\n--- Median bb_lahir per Kelompok Usia ---\n")

median_bb_lahir <- data_step2 %>%
  group_by(kelompok_usia) %>%
  summarise(median_bb = median(bb_lahir, na.rm = TRUE))
print(median_bb_lahir)

cat("\n--- Median kunjungan_posyandu per Kelompok Usia ---\n")

median_kunjungan <- data_step2 %>%
  group_by(kelompok_usia) %>%
  summarise(median_kunjungan = median(kunjungan_posyandu, na.rm = TRUE))
print(median_kunjungan)

cat("\n--- Modus pend_ibu ---\n")

modus_pend_ibu <- data_step2 %>%
  filter(!is.na(pend_ibu)) %>%
  count(pend_ibu, sort = TRUE) %>%
  slice(1) %>%
  pull(pend_ibu)
cat("Modus pend_ibu:", modus_pend_ibu, "\n")

data_step2 <- data_step2 %>%
  left_join(median_bb_lahir,  by = "kelompok_usia") %>%
  left_join(median_kunjungan, by = "kelompok_usia") %>%
  mutate(
    bb_lahir           = ifelse(is.na(bb_lahir),           median_bb,       bb_lahir),
    kunjungan_posyandu = ifelse(is.na(kunjungan_posyandu), median_kunjungan, kunjungan_posyandu),
    pend_ibu           = ifelse(is.na(pend_ibu),           modus_pend_ibu,  pend_ibu)
  ) %>%
  select(-median_bb, -median_kunjungan)

cat("\n--- Missing Values Setelah Imputasi ---\n")

data_step2 %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "NA_Count") %>%
  filter(NA_Count > 0) %>%
  print()

cat("(Jika tidak ada output = tidak ada missing values tersisa)\n")

# ── LANGKAH 3 — Hapus Duplikat ────────────────────────────────────────────────
data_step3 <- data_step2 %>% distinct()

cat("\nSetelah hapus duplikat:", nrow(data_step3), "baris tersisa\n")

cat("Duplikat yang dihapus:", nrow(data_step2) - nrow(data_step3), "\n")

cat("\n=== SUMMARY SEBELUM CLEANING ===\n")

data_stunting_raw %>%
  select(tb_sekarang, bb_sekarang, bb_lahir, kunjungan_posyandu) %>%
  summary() %>% print()

cat("\n=== SUMMARY SESUDAH CLEANING ===\n")

data_step3 %>%
  select(tb_sekarang, bb_sekarang, bb_lahir, kunjungan_posyandu) %>%
  summary() %>% print()

# ── Boxplot Sebelum vs Sesudah Cleaning ──────────────────────────────────────
library(patchwork)

p1a <- ggplot(data_stunting_raw, aes(y = tb_sekarang)) +
  geom_boxplot(fill = "#FFCDD2", outlier.colour = COL_STUNTING,
               colour = COL_STUNTING, width = 0.45) +
  labs(title = "Tinggi Badan\n(Sebelum Cleaning)", y = "cm") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

p1b <- ggplot(data_step3, aes(y = tb_sekarang)) +
  geom_boxplot(fill = "#C8E6C9", outlier.colour = "#2E7D32",
               colour = "#2E7D32", width = 0.45) +
  labs(title = "Tinggi Badan\n(Sesudah Cleaning)", y = "cm") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

p2a <- ggplot(data_stunting_raw, aes(y = bb_sekarang)) +
  geom_boxplot(fill = "#FFCDD2", outlier.colour = COL_STUNTING,
               colour = COL_STUNTING, width = 0.45) +
  labs(title = "Berat Badan\n(Sebelum Cleaning)", y = "kg") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

p2b <- ggplot(data_step3, aes(y = bb_sekarang)) +
  geom_boxplot(fill = "#C8E6C9", outlier.colour = "#2E7D32",
               colour = "#2E7D32", width = 0.45) +
  labs(title = "Berat Badan\n(Sesudah Cleaning)", y = "kg") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

(p1a | p1b | p2a | p2b) +
  plot_annotation(
    title    = "Perbandingan Distribusi Sebelum dan Sesudah Cleaning",
    subtitle = "Merah = sebelum cleaning | Hijau = sesudah cleaning",
    theme    = theme(plot.title    = element_text(face = "bold", colour = "#1A237E"),
                     plot.subtitle = element_text(colour = "#546E7A"))
  )

# ── Simpan data bersih ────────────────────────────────────────────────────────
data_stunting_clean <- data_step3 %>% select(-kelompok_usia)
write_csv(data_stunting_clean, "C:/Users/User/Downloads/data_stunting_clean.csv")
cat("\nFile data_stunting_clean.csv berhasil disimpan!\n")

cat("Dimensi final:", nrow(data_stunting_clean), "baris x",
    ncol(data_stunting_clean), "kolom\n")

# ── IK 10a: Prevalensi Sebelum vs Sesudah Cleaning ───────────────────────────
prev_sebelum <- data_stunting_raw %>%
  count(status_stunting) %>%
  mutate(persen  = round(n / sum(n) * 100, 2),
         dataset = "Sebelum Cleaning")

prev_sesudah <- data_stunting_clean %>%
  count(status_stunting) %>%
  mutate(persen  = round(n / sum(n) * 100, 2),
         dataset = "Sesudah Cleaning")

cat("=== PREVALENSI STUNTING SEBELUM CLEANING ===\n"); print(prev_sebelum)

cat("\n=== PREVALENSI STUNTING SESUDAH CLEANING ===\n"); print(prev_sesudah)

perbandingan <- bind_rows(prev_sebelum, prev_sesudah) %>%
  select(dataset, status_stunting, n, persen)

kable(perbandingan,
      caption   = "Perbandingan Prevalensi Stunting Sebelum dan Sesudah Cleaning",
      col.names = c("Dataset", "Status", "n", "%"),
      align     = "llrr") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  collapse_rows(columns = 1, valign = "middle")

pct_sblm <- prev_sebelum %>% filter(status_stunting == "Stunting") %>% pull(persen)
pct_ssdh <- prev_sesudah %>% filter(status_stunting == "Stunting") %>% pull(persen)
cat("\nPrevalensi Stunting sebelum:", pct_sblm, "%\n")

cat("Prevalensi Stunting sesudah:", pct_ssdh, "%\n")

cat("Selisih                    :", round(pct_ssdh - pct_sblm, 2), "poin persentase\n")

bind_rows(prev_sebelum, prev_sesudah) %>%
  ggplot(aes(x = dataset, y = persen, fill = status_stunting)) +
  geom_bar(stat = "identity", position = "fill",
           colour = "white", linewidth = 0.4, width = 0.5) +
  geom_text(aes(label = paste0(persen, "%")),
            position = position_fill(vjust = 0.5),
            colour = "white", size = 4.2, fontface = "bold") +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title    = "Perbandingan Prevalensi Stunting",
       subtitle = "Sebelum vs Sesudah Pembersihan Data",
       x        = NULL,
       y        = "Proporsi (%)",
       fill     = "Status Stunting",
       caption  = "Sumber: data_stunting_raw.csv & data_stunting_clean.csv")

# Status stunting data outlier yang dihapus
cat("\n=== STATUS STUNTING PADA DATA YANG DIHAPUS (OUTLIER) ===\n")

data_dihapus <- data_stunting_raw %>%
  filter(
    tb_sekarang < 30  | tb_sekarang > 130 |
    bb_sekarang < 2   | bb_sekarang > 20  |
    (!is.na(bb_lahir) & (bb_lahir < 0.5 | bb_lahir > 5.5))
  )

cat("Jumlah data outlier yang dihapus:", nrow(data_dihapus), "\n")

data_dihapus %>%
  count(status_stunting) %>%
  mutate(persen = round(n / sum(n) * 100, 1)) %>%
  print()

# ── Muat data clean ───────────────────────────────────────────────────────────
data_stunting_clean <- read_csv("C:/Users/User/Downloads/data_stunting_clean.csv")
cat("Dimensi data clean:", dim(data_stunting_clean), "\n")

# ── Ordinal Encoding: pend_ibu & Binary Encoding ─────────────────────────────
data_encoded <- data_stunting_clean %>%
  mutate(
    pend_ibu_enc = case_when(
      pend_ibu == "Tidak Sekolah" ~ 0,
      pend_ibu == "SD"            ~ 1,
      pend_ibu == "SMP"           ~ 2,
      pend_ibu == "SMA"           ~ 3,
      pend_ibu == "PT"            ~ 4,
      TRUE ~ NA_real_
    ),
    jenis_kelamin_enc = ifelse(jenis_kelamin == "L", 1, 0),
    akses_air_enc     = ifelse(akses_air     == "Ya", 1, 0),
    sanitasi_enc      = ifelse(sanitasi      == "Baik", 1, 0),
    asi_eksklusif_enc = ifelse(asi_eksklusif == "Ya", 1, 0),
    imunisasi_enc     = ifelse(imunisasi     == "Lengkap", 1, 0)
  )

cat("\n=== Hasil Encoding pend_ibu ===\n")

data_encoded %>%
  count(pend_ibu, pend_ibu_enc) %>%
  arrange(pend_ibu_enc) %>%
  print()

cat("\n=== Hasil Encoding Variabel Binary ===\n")

data_encoded %>%
  select(jenis_kelamin, jenis_kelamin_enc,
         akses_air,     akses_air_enc,
         sanitasi,      sanitasi_enc,
         asi_eksklusif, asi_eksklusif_enc,
         imunisasi,     imunisasi_enc) %>%
  distinct() %>%
  arrange(jenis_kelamin_enc) %>%
  print(n = 20)

# ── Buat 2 Fitur Baru ─────────────────────────────────────────────────────────
data_encoded <- data_encoded %>%
  mutate(
    bmi      = bb_sekarang / (tb_sekarang / 100)^2,
    rasio_bb = bb_sekarang / bb_lahir
  )

cat("\n=== Ringkasan Fitur Baru ===\n")

data_encoded %>%
  select(bmi, rasio_bb, status_stunting) %>%
  group_by(status_stunting) %>%
  summarise(
    n              = n(),
    BMI_mean       = round(mean(bmi,      na.rm=TRUE), 3),
    BMI_median     = round(median(bmi,    na.rm=TRUE), 3),
    BMI_sd         = round(sd(bmi,        na.rm=TRUE), 3),
    RasioBB_mean   = round(mean(rasio_bb, na.rm=TRUE), 3),
    RasioBB_median = round(median(rasio_bb, na.rm=TRUE), 3),
    RasioBB_sd     = round(sd(rasio_bb,  na.rm=TRUE), 3)
  ) %>% print()

cat("\n=== Uji Mann-Whitney: BMI vs status_stunting ===\n")

wilcox.test(bmi ~ status_stunting, data = data_encoded) %>% print()

cat("\n=== Uji Mann-Whitney: rasio_bb vs status_stunting ===\n")

wilcox.test(rasio_bb ~ status_stunting, data = data_encoded) %>% print()

library(patchwork)

p_bmi <- ggplot(data_encoded,
                aes(x = status_stunting, y = bmi, fill = status_stunting)) +
  geom_violin(alpha = 0.35, trim = FALSE, colour = NA) +
  geom_boxplot(width = 0.25, outlier.colour = "#FF6F00",
               outlier.shape = 21, outlier.size = 2,
               outlier.fill = "#FFB300", linewidth = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3.5, fill = "white", colour = COL_ACCENT) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  labs(title = "Indeks Massa Tubuh (BMI)",
       x     = "Status Stunting",
       y     = "BMI (kg/m²)") +
  theme(legend.position = "none")

p_rasio <- ggplot(data_encoded,
                  aes(x = status_stunting, y = rasio_bb, fill = status_stunting)) +
  geom_violin(alpha = 0.35, trim = FALSE, colour = NA) +
  geom_boxplot(width = 0.25, outlier.colour = "#FF6F00",
               outlier.shape = 21, outlier.size = 2,
               outlier.fill = "#FFB300", linewidth = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 23,
               size = 3.5, fill = "white", colour = COL_ACCENT) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  labs(title = "Rasio BB Sekarang / BB Lahir",
       x     = "Status Stunting",
       y     = "Rasio BB") +
  theme(legend.position = "none")

(p_bmi | p_rasio) +
  plot_annotation(
    title    = "Distribusi Fitur Rekayasa per Status Stunting",
    subtitle = "Violin + Boxplot | Berlian putih = rata-rata",
    theme    = theme(plot.title    = element_text(face = "bold", colour = "#1A237E"),
                     plot.subtitle = element_text(colour = "#546E7A"))
  )

# ── Fungsi Min-Max & Normalisasi ──────────────────────────────────────────────
minmax <- function(x) {
  (x - min(x, na.rm=TRUE)) / (max(x, na.rm=TRUE) - min(x, na.rm=TRUE))
}

vars_normalisasi <- c("usia_bulan", "bb_lahir", "tb_sekarang", "bb_sekarang",
                      "pendapatan", "diare_3bln", "kunjungan_posyandu",
                      "zscore_tbu", "bmi", "rasio_bb", "pend_ibu_enc")

data_normalized <- data_encoded %>%
  mutate(across(all_of(vars_normalisasi), minmax, .names = "{.col}_norm"))

cat("\n=== Verifikasi Range Setelah Min-Max Normalisasi ===\n")

data_normalized %>%
  select(ends_with("_norm")) %>%
  summarise(across(everything(), list(
    Min = ~round(min(., na.rm=TRUE), 3),
    Max = ~round(max(., na.rm=TRUE), 3)
  ))) %>%
  pivot_longer(everything(),
               names_to  = c("Variabel","Stat"),
               names_sep = "_(?=[^_]+$)") %>%
  pivot_wider(names_from = Stat, values_from = value) %>%
  print(n = 15)

# ── Deteksi Multikolinearitas dengan VIF ──────────────────────────────────────
library(car)

fitur_model <- data_encoded %>%
  select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
         pendapatan, diare_3bln, kunjungan_posyandu,
         bmi, rasio_bb, pend_ibu_enc,
         jenis_kelamin_enc, akses_air_enc, sanitasi_enc,
         asi_eksklusif_enc, imunisasi_enc,
         status_stunting) %>%
  mutate(status_stunting_num = ifelse(status_stunting == "Stunting", 1, 0)) %>%
  drop_na()

model_vif <- lm(status_stunting_num ~ usia_bulan + bb_lahir + tb_sekarang +
                  bb_sekarang + pendapatan + diare_3bln + kunjungan_posyandu +
                  bmi + rasio_bb + pend_ibu_enc +
                  jenis_kelamin_enc + akses_air_enc + sanitasi_enc +
                  asi_eksklusif_enc + imunisasi_enc,
                data = fitur_model)

vif_result <- vif(model_vif)
vif_df <- data.frame(Fitur = names(vif_result),
                     VIF   = round(vif_result, 3)) %>%
  arrange(desc(VIF))

kable(vif_df,
      caption = "Nilai VIF per Fitur Prediktor",
      col.names = c("Fitur", "VIF"),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(vif_df$VIF > 10),
           background = "#FFEBEE", color = COL_STUNTING, bold = TRUE) %>%
  row_spec(which(vif_df$VIF <= 5 & vif_df$VIF > 0),
           background = "#E8F5E9")

cat("\n=== Pasangan korelasi r > 0.9 ===\n")

cor_fitur <- cor(fitur_model %>%
                   select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
                          pendapatan, diare_3bln, kunjungan_posyandu,
                          bmi, rasio_bb),
                 use = "complete.obs")

cor_melt <- as.data.frame(as.table(cor_fitur)) %>%
  filter(Var1 != Var2, abs(Freq) > 0.9) %>%
  arrange(desc(abs(Freq)))
print(cor_melt)

# ── Dataset Final untuk Modeling ─────────────────────────────────────────────
data_stunting_model <- data_encoded %>%
  select(
    usia_bulan, jenis_kelamin_enc,
    bb_lahir, tb_sekarang, bb_sekarang,
    pend_ibu_enc, pendapatan,
    akses_air_enc, sanitasi_enc, kunjungan_posyandu,
    asi_eksklusif_enc, imunisasi_enc, diare_3bln,
    bmi, rasio_bb,
    status_stunting
  ) %>%
  drop_na()

cat("\n=== Dataset Final untuk Modeling ===\n")

cat("Dimensi:", dim(data_stunting_model), "\n")

cat("Variabel:", names(data_stunting_model), "\n")

cat("\nMissing values:\n")

print(colSums(is.na(data_stunting_model)))

write_csv(data_stunting_model, "C:/Users/User/Downloads/data_stunting_model.csv")
cat("\nFile data_stunting_model.csv berhasil disimpan!\n")

cat("\nHead 5 baris:\n")

print(head(data_stunting_model, 5))

# ── Muat data model ──────────────────────────────────────────────────────────
library(caret)

data_stunting_model <- read_csv("C:/Users/User/Downloads/data_stunting_model.csv")

data_stunting_model <- data_stunting_model %>%
  mutate(status_stunting = factor(status_stunting,
                                   levels = c("Normal", "Stunting")))

cat("Dimensi data:", dim(data_stunting_model), "\n")

cat("Distribusi target:\n")

print(table(data_stunting_model$status_stunting))

print(prop.table(table(data_stunting_model$status_stunting)))

# ── SKENARIO 1: Percentage Split 80:20 ───────────────────────────────────────
set.seed(42)

split_index <- createDataPartition(
  data_stunting_model$status_stunting,
  p    = 0.8,
  list = FALSE
)

train_set <- data_stunting_model[split_index, ]
test_set  <- data_stunting_model[-split_index, ]

cat("\n=== SKENARIO 1: Percentage Split 80:20 ===\n")

cat("Training set:", nrow(train_set), "observasi\n")

cat("Test set    :", nrow(test_set),  "observasi\n")

cat("\nDistribusi kelas — Training set:\n")

print(table(train_set$status_stunting))

print(round(prop.table(table(train_set$status_stunting))*100, 1))

cat("\nDistribusi kelas — Test set:\n")

print(table(test_set$status_stunting))

print(round(prop.table(table(test_set$status_stunting))*100, 1))

# ── SKENARIO 2: Stratified 5-Fold Cross Validation ───────────────────────────
set.seed(42)

cv_folds <- createFolds(
  data_stunting_model$status_stunting,
  k           = 5,
  list        = TRUE,
  returnTrain = FALSE
)

cat("\n=== SKENARIO 2: Stratified 5-Fold Cross Validation ===\n")

cat("Jumlah fold:", length(cv_folds), "\n\n")

for (i in 1:5) {
  fold_test  <- data_stunting_model[cv_folds[[i]], ]
  fold_train <- data_stunting_model[-cv_folds[[i]], ]
  cat(sprintf("Fold %d | Train: %d | Test: %d | Stunting di test: %d (%.1f%%)\n",
              i,
              nrow(fold_train),
              nrow(fold_test),
              sum(fold_test$status_stunting == "Stunting"),
              mean(fold_test$status_stunting == "Stunting") * 100))
}

ctrl_cv <- trainControl(
  method          = "cv",
  number          = 5,
  classProbs      = TRUE,
  summaryFunction = twoClassSummary,
  savePredictions = "final",
  sampling        = NULL
)

cat("\n=== trainControl (untuk UK-9) sudah disiapkan ===\n")

cat("Method: 5-fold CV | ClassProbs: TRUE | Summary: ROC\n")

saveRDS(train_set, "C:/Users/User/Downloads/train_set.rds")
saveRDS(test_set,  "C:/Users/User/Downloads/test_set.rds")
saveRDS(ctrl_cv,   "C:/Users/User/Downloads/ctrl_cv.rds")
cat("\nFile train_set.rds, test_set.rds, ctrl_cv.rds tersimpan!\n")

# ── Muat data & set factor ────────────────────────────────────────────────────
data_stunting_model <- read_csv("C:/Users/User/Downloads/data_stunting_model.csv") %>%
  mutate(status_stunting = factor(status_stunting, levels = c("Normal", "Stunting")))

train_set <- readRDS("C:/Users/User/Downloads/train_set.rds") %>%
  mutate(status_stunting = factor(status_stunting, levels = c("Normal", "Stunting")))
test_set  <- readRDS("C:/Users/User/Downloads/test_set.rds") %>%
  mutate(status_stunting = factor(status_stunting, levels = c("Normal", "Stunting")))

cat("Train:", nrow(train_set), "| Test:", nrow(test_set), "\n")

# ── Bangun Logistic Regression ────────────────────────────────────────────────
set.seed(42)
model_lr <- glm(status_stunting ~ .,
                data   = train_set,
                family = binomial(link = "logit"))

cat("\n=== RINGKASAN MODEL LOGISTIC REGRESSION ===\n")

summary(model_lr)

# ── Odds Ratio ────────────────────────────────────────────────────────────────
cat("\n=== ODDS RATIO (exp(coef)) ===\n")

or_df <- data.frame(
  Variabel  = names(coef(model_lr)),
  Koefisien = round(coef(model_lr), 4),
  OddsRatio = round(exp(coef(model_lr)), 4),
  CI_Lower  = round(exp(confint(model_lr))[,1], 4),
  CI_Upper  = round(exp(confint(model_lr))[,2], 4)
)
print(or_df)

cat("\n=== VARIABEL SIGNIFIKAN (p < 0.05) ===\n")

coef_summary <- summary(model_lr)$coefficients
sig_vars     <- coef_summary[coef_summary[,4] < 0.05, ]
print(round(sig_vars, 4))

# ── Decision Tree dengan Tuning cp ────────────────────────────────────────────
library(rpart)
library(rpart.plot)

set.seed(42)

ctrl_dt <- trainControl(
  method          = "cv",
  number          = 5,
  classProbs      = TRUE,
  summaryFunction = twoClassSummary,
  savePredictions = "final"
)

cp_grid <- expand.grid(cp = seq(0.001, 0.05, by = 0.001))

model_dt_cv <- train(
  status_stunting ~ .,
  data      = train_set,
  method    = "rpart",
  trControl = ctrl_dt,
  tuneGrid  = cp_grid,
  metric    = "ROC"
)

cat("\n=== HASIL TUNING cp (5-Fold CV) ===\n")

print(model_dt_cv$results[, c("cp","ROC","Sens","Spec")])

cat("\n=== cp OPTIMAL (best) ===\n")

print(model_dt_cv$bestTune)

# 1-SE Rule
results_dt <- model_dt_cv$results
best_roc   <- max(results_dt$ROC)
se_roc     <- sd(results_dt$ROC) / sqrt(5)
threshold  <- best_roc - se_roc
cp_1se     <- results_dt$cp[results_dt$ROC >= threshold][1]

cat("\nBest ROC       :", round(best_roc,  4), "\n")

cat("1-SE threshold :", round(threshold, 4), "\n")

cat("cp optimal (1-SE rule):", cp_1se, "\n")

set.seed(42)
model_dt_final <- rpart(
  status_stunting ~ .,
  data    = train_set,
  method  = "class",
  control = rpart.control(cp = cp_1se)
)

cat("\n=== RINGKASAN DECISION TREE FINAL ===\n")

printcp(model_dt_final)

print(model_dt_final$variable.importance)

rpart.plot(model_dt_final,
           type          = 4,
           extra         = 104,
           under         = TRUE,
           fallen.leaves = TRUE,
           branch.lty    = 3,
           shadow.col    = "gray70",
           col           = COL_ACCENT,
           main          = paste0("Decision Tree Prediksi Stunting (cp = ", cp_1se, ")"),
           cex           = 0.78,
           box.palette   = c(COL_NORMAL, COL_STUNTING))

# ── Random Forest dengan Tuning mtry ─────────────────────────────────────────
library(randomForest)

set.seed(42)

ctrl_rf <- trainControl(
  method          = "cv",
  number          = 5,
  classProbs      = TRUE,
  summaryFunction = twoClassSummary,
  savePredictions = "final"
)

mtry_grid <- expand.grid(mtry = c(2, 3, 4, 5, 6, 7, 8))

model_rf_cv <- train(
  status_stunting ~ .,
  data      = train_set,
  method    = "rf",
  trControl = ctrl_rf,
  tuneGrid  = mtry_grid,
  metric    = "ROC",
  ntree     = 200
)

cat("\n=== HASIL TUNING mtry (5-Fold CV) ===\n")

print(model_rf_cv$results[, c("mtry","ROC","Sens","Spec")])

cat("\n=== mtry OPTIMAL ===\n")

print(model_rf_cv$bestTune)

plot(model_rf_cv, main = "Tuning mtry Random Forest (5-Fold CV)")

cat("\n=== FEATURE IMPORTANCE ===\n")

imp_rf <- varImp(model_rf_cv)
print(imp_rf)

# ── Feature Importance yang lebih cantik ─────────────────────────────────────
imp_df <- imp_rf$importance %>%
  as.data.frame() %>%
  rownames_to_column("Fitur") %>%
  rename(Importance = Overall) %>%
  arrange(Importance)

ggplot(imp_df, aes(x = reorder(Fitur, Importance), y = Importance)) +
  geom_col(aes(fill = Importance), colour = "white", width = 0.7) +
  geom_text(aes(label = round(Importance, 1)), hjust = -0.1, size = 3.2) +
  scale_fill_gradient(low = "#BBDEFB", high = COL_NORMAL) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title    = "Feature Importance — Random Forest",
       subtitle = paste0("mtry optimal = ", model_rf_cv$bestTune$mtry),
       x        = NULL,
       y        = "Importance Score (ROC)",
       fill     = "Importance") +
  theme(legend.position = "none")

# ── Penanganan Class Imbalance dengan ROSE ────────────────────────────────────
library(ROSE)

set.seed(42)

cat("Distribusi kelas training (sebelum ROSE):\n")

print(table(train_set$status_stunting))

train_rose <- ROSE(status_stunting ~ ., data = train_set, seed = 42)$data

cat("\nDistribusi kelas training (sesudah ROSE):\n")

print(table(train_rose$status_stunting))

# RF tanpa ROSE
pred_rf_no <- predict(model_rf_cv, newdata = test_set)
cm_no      <- confusionMatrix(pred_rf_no, test_set$status_stunting, positive = "Stunting")

# RF dengan ROSE
set.seed(42)
model_rf_rose <- train(
  status_stunting ~ .,
  data      = train_rose,
  method    = "rf",
  trControl = trainControl(method = "cv", number = 5,
                           classProbs = TRUE,
                           summaryFunction = twoClassSummary),
  tuneGrid  = expand.grid(mtry = model_rf_cv$bestTune$mtry),
  metric    = "ROC",
  ntree     = 200
)

pred_rf_rose <- predict(model_rf_rose, newdata = test_set)
cm_rose      <- confusionMatrix(pred_rf_rose, test_set$status_stunting, positive = "Stunting")

cat("\n=== CONFUSION MATRIX — RF TANPA ROSE ===\n"); print(cm_no)

cat("\n=== CONFUSION MATRIX — RF DENGAN ROSE ===\n"); print(cm_rose)

cat("\n=== PERBANDINGAN SENSITIVITY (Recall) ===\n")

cat("Tanpa ROSE:", round(cm_no$byClass["Sensitivity"],   4), "\n")

cat("Dengan ROSE:", round(cm_rose$byClass["Sensitivity"], 4), "\n")

cat("\n=== PERBANDINGAN SPECIFICITY ===\n")

cat("Tanpa ROSE:", round(cm_no$byClass["Specificity"],   4), "\n")

cat("Dengan ROSE:", round(cm_rose$byClass["Specificity"], 4), "\n")

# ── Load dan persiapkan data ──────────────────────────────────────────────────
library(pROC)

test_set  <- readRDS("C:/Users/User/Downloads/test_set.rds") %>%
  mutate(status_stunting = factor(status_stunting, levels = c("Normal","Stunting")))
train_set <- readRDS("C:/Users/User/Downloads/train_set.rds") %>%
  mutate(status_stunting = factor(status_stunting, levels = c("Normal","Stunting")))

# ── MODEL 1: Logistic Regression ──────────────────────────────────────────────
set.seed(42)
model_lr <- glm(status_stunting ~ ., data = train_set,
                family = binomial(link = "logit"))

pred_lr_class <- predict(model_lr, newdata = test_set, type = "response")
pred_lr_label <- factor(ifelse(pred_lr_class >= 0.5, "Stunting", "Normal"),
                        levels = c("Normal", "Stunting"))

cm_lr  <- confusionMatrix(pred_lr_label, test_set$status_stunting, positive = "Stunting")
roc_lr <- roc(test_set$status_stunting, pred_lr_class,
              levels = c("Normal","Stunting"), quiet = TRUE)

cat("=== CONFUSION MATRIX — Logistic Regression ===\n"); print(cm_lr)

cat("AUC-ROC LR:", round(auc(roc_lr), 4), "\n")

# ── MODEL 2: Decision Tree ────────────────────────────────────────────────────
set.seed(42)
model_dt_final <- rpart(status_stunting ~ ., data = train_set,
                        method  = "class",
                        control = rpart.control(cp = 0.001))

pred_dt_prob  <- predict(model_dt_final, newdata = test_set, type = "prob")[, "Stunting"]
pred_dt_label <- predict(model_dt_final, newdata = test_set, type = "class")

cm_dt  <- confusionMatrix(pred_dt_label, test_set$status_stunting, positive = "Stunting")
roc_dt <- roc(test_set$status_stunting, pred_dt_prob,
              levels = c("Normal","Stunting"), quiet = TRUE)

cat("\n=== CONFUSION MATRIX — Decision Tree ===\n"); print(cm_dt)

cat("AUC-ROC DT:", round(auc(roc_dt), 4), "\n")

# ── MODEL 3: Random Forest ────────────────────────────────────────────────────
set.seed(42)
model_rf_final <- randomForest(status_stunting ~ ., data = train_set,
                               mtry = 6, ntree = 200)

pred_rf_prob  <- predict(model_rf_final, newdata = test_set, type = "prob")[, "Stunting"]
pred_rf_label <- predict(model_rf_final, newdata = test_set, type = "class")

cm_rf  <- confusionMatrix(pred_rf_label, test_set$status_stunting, positive = "Stunting")
roc_rf <- roc(test_set$status_stunting, pred_rf_prob,
              levels = c("Normal","Stunting"), quiet = TRUE)

cat("\n=== CONFUSION MATRIX — Random Forest ===\n"); print(cm_rf)

cat("AUC-ROC RF:", round(auc(roc_rf), 4), "\n")

# ── Rekap Perbandingan Metrik ─────────────────────────────────────────────────
rekap <- data.frame(
  Model       = c("Logistic Regression", "Decision Tree", "Random Forest"),
  Accuracy    = c(round(cm_lr$overall["Accuracy"], 4),
                  round(cm_dt$overall["Accuracy"], 4),
                  round(cm_rf$overall["Accuracy"], 4)),
  Sensitivity = c(round(cm_lr$byClass["Sensitivity"], 4),
                  round(cm_dt$byClass["Sensitivity"], 4),
                  round(cm_rf$byClass["Sensitivity"], 4)),
  Specificity = c(round(cm_lr$byClass["Specificity"], 4),
                  round(cm_dt$byClass["Specificity"], 4),
                  round(cm_rf$byClass["Specificity"], 4)),
  Precision   = c(round(cm_lr$byClass["Precision"], 4),
                  round(cm_dt$byClass["Precision"], 4),
                  round(cm_rf$byClass["Precision"], 4)),
  F1_Score    = c(round(cm_lr$byClass["F1"], 4),
                  round(cm_dt$byClass["F1"], 4),
                  round(cm_rf$byClass["F1"], 4)),
  AUC_ROC     = c(round(auc(roc_lr), 4),
                  round(auc(roc_dt), 4),
                  round(auc(roc_rf), 4)),
  Kappa       = c(round(cm_lr$overall["Kappa"], 4),
                  round(cm_dt$overall["Kappa"], 4),
                  round(cm_rf$overall["Kappa"], 4))
)

kable(rekap,
      caption = "Rekap Perbandingan Metrik Evaluasi — Ketiga Model",
      align   = c("l", rep("r", 7))) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  column_spec(1, bold = TRUE, background = "#E8EAF6") %>%
  column_spec(3, background = "#E8F5E9",
              color = ifelse(rekap$Sensitivity == max(rekap$Sensitivity),
                             "#1B5E20", "black"),
              bold  = ifelse(rekap$Sensitivity == max(rekap$Sensitivity), TRUE, FALSE)) %>%
  column_spec(7, background = "#E3F2FD",
              color = ifelse(rekap$AUC_ROC == max(rekap$AUC_ROC),
                             "#0D47A1", "black"),
              bold  = ifelse(rekap$AUC_ROC == max(rekap$AUC_ROC), TRUE, FALSE)) %>%
  add_header_above(c(" " = 1, "Performa Model" = 7))