library(readr); library(tidyverse); library(knitr); library(kableExtra)

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

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) %>%
  row_spec(which(tabel_target$status_stunting == "Normal"),
           background = "#E3F2FD", color = "#0D47A1") %>%
  row_spec(0, background = "#1A237E", color = "white", bold = TRUE)

pct_min <- tabel_target %>% filter(status_stunting == "Stunting") %>% pull(persentase)
cat(sprintf(
  "Kelas minoritas (Stunting) : %s%%\nKelas mayoritas (Normal)   : %s%%\nRasio imbalance            : %.2f : 1\nKategori imbalance         : Mild Imbalance (20–40%%)\n",
  pct_min, 100 - pct_min, (100 - pct_min) / pct_min))

ggplot(data_stunting_raw, aes(x = status_stunting, fill = status_stunting)) +
  geom_bar(width = 0.5, colour = "white", linewidth = 0.5) +
  geom_text(stat = "count",
            aes(label = paste0(..count.., "\n(",
                               round(..count.. / nrow(data_stunting_raw) * 100, 1), "%)")),
            vjust = -0.3, size = 4.5, fontface = "bold") +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18)), limits = c(0, NA)) +
  labs(title    = "Distribusi Variabel Target: Status Stunting",
       subtitle = "Dataset data_stunting_raw.csv — 500 Balita, 50 Kab/Kota Prioritas",
       x = "Status Stunting", y = "Jumlah Balita", fill = "Status",
       caption = "Sumber: data_stunting_raw.csv | Rasio Imbalance 2,88 : 1 (Mild Imbalance)") +
  annotate("label", x = 1.5, y = 420,
           label = "Rasio Imbalance\n2,88 : 1\n(Mild Imbalance)",
           fill = "#FFF9C4", colour = "#E65100", size = 3.2,
           label.size = 0.4, fontface = "italic")

str(data_stunting_raw)

head(data_stunting_raw, 5)

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, digits = 2,
      caption = "Statistik Deskriptif Variabel Numerik",
      align   = c("l", rep("r", 8))) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE) %>%
  column_spec(1, bold = TRUE, background = "#E8EAF6") %>%
  column_spec(9, background = "#FFF9C4",
              color = ifelse(abs(as.numeric(stat_desk$Skew)) >= 1, COL_STUNTING, "black"),
              bold  = ifelse(abs(as.numeric(stat_desk$Skew)) >= 1, TRUE, FALSE)) %>%
  add_header_above(c(" " = 1, "Nilai Ekstrem" = 2, "Pemusatan" = 2,
                     "Sebaran" = 3, "Bentuk" = 1)) %>%
  row_spec(0, background = "#1A237E", color = "white", bold = TRUE)

ggplot(data_stunting_raw, aes(x = zscore_tbu, fill = status_stunting)) +
  geom_histogram(bins = 42, alpha = 0.75, position = "identity",
                 colour = "white", linewidth = 0.15) +
  geom_vline(xintercept = -2, linetype = "dashed", colour = COL_ORANGE, linewidth = 1) +
  annotate("rect", xmin = -Inf, xmax = -2, ymin = -Inf, ymax = Inf,
           fill = "#FFCCBC", alpha = 0.12) +
  annotate("label", x = -4.2, y = 25, label = "Zona Stunting\n(< −2 SD)",
           fill = "#FFF3E0", colour = COL_ORANGE, size = 3, label.size = 0.35) +
  scale_fill_manual(values = c("Normal" = COL_NORMAL, "Stunting" = COL_STUNTING)) +
  scale_x_continuous(breaks = seq(-6, 4, 1)) +
  coord_cartesian(xlim = c(-8, 6)) +
  labs(title    = "Distribusi Z-Score TB/U berdasarkan Status Stunting",
       subtitle = "Outlier ekstrem (di luar rentang ±6 SD) terpotong untuk keterbacaan",
       x = "Z-Score TB/U", y = "Frekuensi", fill = "Status Stunting",
       caption = "Garis putus-putus oranye: cut-off WHO −2 SD")

library(patchwork)

make_box <- function(var, label, unit) {
  ggplot(data_stunting_raw, aes(x = status_stunting, y = .data[[var]],
                                 fill = status_stunting)) +
    geom_violin(alpha = 0.25, trim = TRUE, colour = NA) +
    geom_boxplot(width = 0.3, outlier.colour = COL_ORANGE, outlier.shape = 21,
                 outlier.size = 1.8, outlier.fill = "#FFB300", linewidth = 0.45) +
    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 = label, x = "Status Stunting", y = unit) +
    theme(legend.position = "none")
}

p1 <- make_box("tb_sekarang", "Tinggi Badan (cm)", "Tinggi Badan (cm)")
p2 <- make_box("bb_sekarang", "Berat Badan (kg)",  "Berat Badan (kg)")
p3 <- make_box("zscore_tbu",  "Z-Score TB/U",       "Z-Score TB/U")

(p1 | p2 | p3) +
  plot_annotation(
    title    = "Distribusi Variabel Antropometri per Status Stunting",
    subtitle = "Violin + Boxplot | Berlian putih = rata-rata | Titik oranye = outlier",
    theme    = theme(plot.title    = element_text(face = "bold", colour = "#1A237E", size = 13),
                     plot.subtitle = element_text(colour = "#546E7A", size = 9.5))
  )

ggplot(data_stunting_raw %>% filter(zscore_tbu > -8, zscore_tbu < 6),
       aes(x = tb_sekarang, y = zscore_tbu, colour = status_stunting)) +
  geom_point(alpha = 0.5, size = 1.8, shape = 16) +
  geom_smooth(aes(group = status_stunting), method = "lm", se = TRUE,
              linewidth = 0.9) +
  geom_hline(yintercept = -2, linetype = "dashed", colour = COL_ORANGE, linewidth = 0.9) +
  annotate("label", x = 115, y = -1.55, label = "Cut-off −2 SD",
           fill = "#FFF3E0", colour = COL_ORANGE, 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 = "Korelasi positif sangat kuat (r = 0,880) | Outlier ekstrem dipotong",
       x = "Tinggi Badan Sekarang (cm)", y = "Z-Score TB/U",
       colour = "Status Stunting",
       caption = "Sumber: data_stunting_raw.csv")

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",
      col.names = c("Variabel","Nilai Skewness","Kategori")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(skew_df$Kategori == "Skewed Berat"),
           background = "#FFEBEE", color = COL_STUNTING, bold = TRUE) %>%
  row_spec(which(skew_df$Kategori == "Simetris"),
           background = "#E8F5E9") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE)

buat_tabel_proporsi <- function(data, var) {
  data %>%
    group_by(across(all_of(c(var, "status_stunting")))) %>%
    summarise(n = n(), .groups = "drop") %>%
    group_by(across(all_of(var))) %>%
    mutate(proporsi = round(n / sum(n) * 100, 1)) %>%
    ungroup()
}

for (v in c("pend_ibu","akses_air","sanitasi","asi_eksklusif","imunisasi")) {
  tbl <- buat_tabel_proporsi(data_stunting_raw, v)
  cat("\n\n**Tabel Proporsi:", v, "vs status_stunting**\n\n")
  print(kable(tbl, caption = paste("Proporsi:", v, "vs Status Stunting")) %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed"),
                  full_width = FALSE, position = "center"))
}

data_stunting_raw %>%
  filter(!is.na(pend_ibu)) %>%
  mutate(pend_ibu = factor(pend_ibu,
         levels = c("Tidak Sekolah","SD","SMP","SMA","PT"))) %>%
  count(pend_ibu, status_stunting) %>%
  group_by(pend_ibu) %>%
  mutate(pct = round(n / sum(n) * 100, 1)) %>%
  ggplot(aes(x = pend_ibu, y = pct, fill = status_stunting)) +
  geom_col(position = "fill", colour = "white", linewidth = 0.4, width = 0.65) +
  geom_text(aes(label = paste0(pct, "%")),
            position = position_fill(vjust = 0.5),
            colour = "white", size = 3.5, fontface = "bold") +
  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 Tingkat Pendidikan Ibu",
       subtitle = "Semakin tinggi pendidikan ibu → prevalensi stunting cenderung menurun",
       x = "Tingkat Pendidikan Ibu", y = "Proporsi (%)", fill = "Status Stunting",
       caption = "Prevalensi stunting tertinggi: SD (31,8%) | terendah: PT (15,9%)")

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")) %>%
  count(variabel, kategori, status_stunting) %>%
  group_by(variabel, kategori) %>%
  mutate(pct = round(n / sum(n) * 100, 1)) %>%
  ggplot(aes(x = kategori, y = pct, fill = status_stunting)) +
  geom_col(position = "fill", colour = "white", linewidth = 0.35, width = 0.6) +
  geom_text(aes(label = paste0(pct, "%")),
            position = position_fill(vjust = 0.5),
            colour = "white", size = 3.2, fontface = "bold") +
  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 kesehatan dan perilaku kesehatan",
       x = NULL, y = "Proporsi (%)", fill = "Status Stunting") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

outlier_tb <- data_stunting_raw %>%
  filter(tb_sekarang < 61.2 | tb_sekarang > 124.0) %>%
  select(id_balita, usia_bulan, tb_sekarang, status_stunting)

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

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

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

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

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

print(head(outlier_tb, 10))

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

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

print(outlier_bb)

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

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

print(outlier_bbl)

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

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

print(outlier_zscore)

rekap_outlier <- data.frame(
  Variabel    = c("tb_sekarang","bb_sekarang","bb_lahir","zscore_tbu","TOTAL"),
  Jumlah      = c(nrow(outlier_tb), nrow(outlier_bb), nrow(outlier_bbl),
                  nrow(outlier_zscore),
                  nrow(outlier_tb)+nrow(outlier_bb)+nrow(outlier_bbl)+nrow(outlier_zscore)),
  Persen      = round(c(nrow(outlier_tb), nrow(outlier_bb), nrow(outlier_bbl),
                         nrow(outlier_zscore),
                         nrow(outlier_tb)+nrow(outlier_bb)+nrow(outlier_bbl)+nrow(outlier_zscore)
                         ) / 500 * 100, 1),
  Batas_WHO   = c("61,2–124,0 cm","5,5–24,0 kg","2,0–4,5 kg","−6 s/d +6","—")
)

kable(rekap_outlier, caption = "Rekap Total Nilai Tidak Valid Berdasarkan Batas WHO",
      col.names = c("Variabel","Jumlah Tidak Valid","Persentase (%)","Batas WHO")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE,
                position = "center") %>%
  row_spec(nrow(rekap_outlier), background = "#B71C1C", color = "white", bold = TRUE) %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE)

cat("\n=== Missing Values per Variabel ===\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()

n_obs <- nrow(data_stunting_raw); n_pred <- 14; min_obs <- n_pred * 10
cat(sprintf(
  "\n=== Kecukupan Data (Rule of Thumb: 10 obs/prediktor) ===\nObservasi       : %d\nPrediktor aktif : %d\nMinimum (10x)   : %d\nRasio           : %.2f\nKesimpulan      : %s\n",
  n_obs, n_pred, min_obs, n_obs/n_pred,
  ifelse(n_obs >= min_obs, "CUKUP ✓ (rasio 35,71 ≥ 10)", "TIDAK CUKUP ✗")))

cat("=== Dimensi Dataset ===\n")

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

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

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

cat("\n=== Simulasi Stratified Split 80:20 ===\n")

cat("Training (80%): Normal ≈", round(371*0.8), "| Stunting ≈", round(129*0.8), "\n")

cat("Test     (20%): Normal ≈", round(371*0.2), "| Stunting ≈", round(129*0.2), "\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()

# ── Kondisi 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()

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(sprintf("\n--- Outlier (Batas Soal) ---\ntb_sekarang  (30-130 cm): %d\nbb_sekarang  (2-20 kg)  : %d\nbb_lahir     (0.5-5.5 kg): %d\nDuplikat: %d\n",
    n_tb, n_bb, n_bbl, sum(duplicated(data_stunting_raw))))

# ── Langkah 1: Hapus outlier ──────────────────────────────────────────────────
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:", nrow(data_step1), "baris (dihapus:",
    nrow(data_stunting_raw) - nrow(data_step1), ")\n")

# ── Langkah 2: Imputasi ───────────────────────────────────────────────────────
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"))

median_bb  <- data_step2 %>% group_by(kelompok_usia) %>%
              summarise(med_bb = median(bb_lahir, na.rm=TRUE))
median_kun <- data_step2 %>% group_by(kelompok_usia) %>%
              summarise(med_kun = median(kunjungan_posyandu, na.rm=TRUE))
modus_pend <- data_step2 %>% filter(!is.na(pend_ibu)) %>%
              count(pend_ibu, sort=TRUE) %>% slice(1) %>% pull(pend_ibu)

cat("\nModus pend_ibu:", modus_pend, "\n")

data_step2 <- data_step2 %>%
  left_join(median_bb,  by = "kelompok_usia") %>%
  left_join(median_kun, by = "kelompok_usia") %>%
  mutate(
    bb_lahir           = ifelse(is.na(bb_lahir),           med_bb,  bb_lahir),
    kunjungan_posyandu = ifelse(is.na(kunjungan_posyandu), med_kun, kunjungan_posyandu),
    pend_ibu           = ifelse(is.na(pend_ibu),           modus_pend, pend_ibu)
  ) %>% select(-med_bb, -med_kun)

cat("Missing setelah imputasi:\n")

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

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

# ── Langkah 3: Hapus duplikat ─────────────────────────────────────────────────
data_step3 <- data_step2 %>% distinct()
cat("\nSetelah hapus duplikat:", nrow(data_step3), "baris\n")

cat("\n=== SUMMARY SEBELUM vs SESUDAH ===\n")

cat("--- Sebelum ---\n")

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

cat("--- Sesudah ---\n")

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

# ── Simpan ────────────────────────────────────────────────────────────────────
data_stunting_clean <- data_step3 %>% select(-kelompok_usia)
write_csv(data_stunting_clean, "C:/Users/User/Downloads/data_stunting_clean.csv")
cat("\n✓ data_stunting_clean.csv tersimpan |",
    nrow(data_stunting_clean), "baris ×", ncol(data_stunting_clean), "kolom\n")

make_box_clean <- function(df, var, unit, label, col_fill) {
  ggplot(df, aes(y = .data[[var]])) +
    geom_boxplot(fill = col_fill, colour = darken_col(col_fill, 0.3),
                 outlier.colour = COL_ORANGE, outlier.shape = 21,
                 outlier.size = 1.5, width = 0.45) +
    labs(title = label, y = unit) +
    theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
}

darken_col <- function(col, factor = 0.3) {
  rgb_col <- col2rgb(col) / 255
  rgb(rgb_col[1]*factor, rgb_col[2]*factor, rgb_col[3]*factor)
}

p1a <- ggplot(data_stunting_raw, aes(y = tb_sekarang)) +
  geom_boxplot(fill = "#FFCDD2", colour = COL_STUNTING,
               outlier.colour = COL_ORANGE, outlier.shape = 21, width = 0.45) +
  labs(title = "Tinggi Badan\n(Sebelum)", 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", colour = COL_GREEN,
               outlier.colour = COL_ORANGE, outlier.shape = 21, width = 0.45) +
  labs(title = "Tinggi Badan\n(Sesudah)", 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", colour = COL_STUNTING,
               outlier.colour = COL_ORANGE, outlier.shape = 21, width = 0.45) +
  labs(title = "Berat Badan\n(Sebelum)", 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", colour = COL_GREEN,
               outlier.colour = COL_ORANGE, outlier.shape = 21, width = 0.45) +
  labs(title = "Berat Badan\n(Sesudah)", y = "kg") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

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

prev_sblm <- data_stunting_raw %>%
  count(status_stunting) %>%
  mutate(persen = round(n/sum(n)*100, 2), dataset = "Sebelum Cleaning")

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

perbandingan <- bind_rows(prev_sblm, prev_ssdh) %>%
  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") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE)

pct_s <- prev_sblm %>% filter(status_stunting == "Stunting") %>% pull(persen)
pct_a <- prev_ssdh %>% filter(status_stunting == "Stunting") %>% pull(persen)
cat(sprintf("Prevalensi Stunting sebelum : %s%%\nPrevalensi Stunting sesudah : %s%%\nSelisih                     : %s poin persentase\n",
    pct_s, pct_a, round(pct_a - pct_s, 2)))

bind_rows(prev_sblm, prev_ssdh) %>%
  ggplot(aes(x = dataset, y = persen, fill = status_stunting)) +
  geom_col(position = "fill", colour = "white", linewidth = 0.5, 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")

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

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("=== Hasil Encoding pend_ibu (Ordinal) ===\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)

data_encoded <- data_encoded %>%
  mutate(
    bmi      = bb_sekarang / (tb_sekarang / 100)^2,
    rasio_bb = bb_sekarang / bb_lahir
  )

cat("=== Ringkasan 2 Fitur Baru per Status Stunting ===\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)

make_violin <- function(var, label, unit) {
  ggplot(data_encoded, aes(x = status_stunting, y = .data[[var]],
                            fill = status_stunting)) +
    geom_violin(alpha = 0.3, trim = FALSE, colour = NA) +
    geom_boxplot(width = 0.22, outlier.colour = COL_ORANGE, outlier.shape = 21,
                 outlier.size = 1.8, outlier.fill = "#FFB300", linewidth = 0.45) +
    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 = label, x = "Status Stunting", y = unit) +
    theme(legend.position = "none")
}

(make_violin("bmi",      "Indeks Massa Tubuh (BMI)",          "BMI (kg/m²)") |
 make_violin("rasio_bb", "Rasio BB Sekarang / BB Lahir", "Rasio BB")) +
  plot_annotation(
    title    = "Distribusi Fitur Rekayasa per Status Stunting",
    subtitle = "Violin = distribusi | Box = kuartil | Berlian putih = rata-rata",
    theme    = theme(plot.title    = element_text(face = "bold", colour = "#1A237E", size = 12),
                     plot.subtitle = element_text(colour = "#546E7A", size = 9))
  )

minmax <- function(x) (x - min(x, na.rm=TRUE)) / (max(x, na.rm=TRUE) - min(x, na.rm=TRUE))

vars_norm <- 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_norm), minmax, .names = "{.col}_norm"))

cat("=== 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)

# ── Library ───────────────────────────────────────────────────────────────────
library(car)
library(dplyr)
library(kableExtra)

# ── Persiapan data & helper ───────────────────────────────────────────────────
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()

semua_prediktor <- c("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")

kandidat_drop <- c("tb_sekarang", "bb_sekarang", "bmi", "rasio_bb")

hitung_vif <- function(data, prediktor) {
  formula_str <- paste("status_stunting_num ~", paste(prediktor, collapse = " + "))
  vif(lm(as.formula(formula_str), data = data))
}

render_vif_table <- function(vif_vals, caption_text) {
  df <- data.frame(Fitur = names(vif_vals), VIF = round(vif_vals, 3)) %>%
    arrange(desc(VIF))
  kable(df, caption = caption_text,
        col.names = c("Fitur", "VIF"), align = "lr", booktabs = TRUE) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE, position = "center") %>%
    row_spec(0, background = "#37474F", color = "white", bold = TRUE) %>%
    row_spec(which(df$VIF > 10), background = "#FFEBEE",
             color = COL_STUNTING, bold = TRUE) %>%
    row_spec(which(df$VIF <= 5), background = "#E8F5E9")
}

# ══════════════════════════════════════════════════════════════════════════════
# [0] BASELINE — Semua variabel
# ══════════════════════════════════════════════════════════════════════════════
vif_base    <- hitung_vif(fitur_model, semua_prediktor)
vif_base_df <- data.frame(Fitur = names(vif_base), VIF = round(vif_base, 3)) %>%
  arrange(desc(VIF))

render_vif_table(vif_base, "Baseline — Nilai VIF Semua Prediktor")

# ══════════════════════════════════════════════════════════════════════════════
# [1] ITERASI 1 — Keluarkan 1 variabel
# ══════════════════════════════════════════════════════════════════════════════
hasil_iter1 <- lapply(kandidat_drop, function(var) {
  v <- hitung_vif(fitur_model, setdiff(semua_prediktor, var))
  data.frame(Variabel_Dikeluarkan = var,
             Max_VIF              = round(max(v), 3),
             N_VIF_gt10           = sum(v > 10))
}) %>% bind_rows() %>% arrange(N_VIF_gt10, Max_VIF)

kable(hasil_iter1,
      caption = "Iterasi 1 — Efek Pengeluaran 1 Variabel terhadap VIF",
      align = "lrr", booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE) %>%
  row_spec(1, background = "#E3F2FD", bold = TRUE)

# ══════════════════════════════════════════════════════════════════════════════
# [2] ITERASI 2 — Keluarkan 2 variabel
# ══════════════════════════════════════════════════════════════════════════════
hasil_iter2 <- combn(kandidat_drop, 2, simplify = FALSE) %>%
  lapply(function(combo) {
    v <- hitung_vif(fitur_model, setdiff(semua_prediktor, combo))
    data.frame(Variabel_Dikeluarkan = paste(combo, collapse = " + "),
               Max_VIF              = round(max(v), 3),
               N_VIF_gt10           = sum(v > 10))
  }) %>% bind_rows() %>% arrange(N_VIF_gt10, Max_VIF)

kable(hasil_iter2,
      caption = "Iterasi 2 — Efek Pengeluaran 2 Variabel terhadap VIF",
      align = "lrr", booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE) %>%
  row_spec(1, background = "#E3F2FD", bold = TRUE)

# ══════════════════════════════════════════════════════════════════════════════
# [3] ITERASI 3 — Keluarkan 3 variabel
# ══════════════════════════════════════════════════════════════════════════════
hasil_iter3 <- combn(kandidat_drop, 3, simplify = FALSE) %>%
  lapply(function(combo) {
    v <- hitung_vif(fitur_model, setdiff(semua_prediktor, combo))
    data.frame(Variabel_Dikeluarkan = paste(combo, collapse = " + "),
               Max_VIF              = round(max(v), 3),
               N_VIF_gt10           = sum(v > 10))
  }) %>% bind_rows() %>% arrange(N_VIF_gt10, Max_VIF)

kable(hasil_iter3,
      caption = "Iterasi 3 — Efek Pengeluaran 3 Variabel terhadap VIF",
      align = "lrr", booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE) %>%
  row_spec(1, background = "#E3F2FD", bold = TRUE)

# ══════════════════════════════════════════════════════════════════════════════
# [4] RINGKASAN KOMPARATIF & KEPUTUSAN OPTIMAL
# ══════════════════════════════════════════════════════════════════════════════
ringkasan <- data.frame(
  Skenario      = c("Baseline (0 dikeluarkan)",
                    paste0("Iterasi 1: keluarkan ", hasil_iter1$Variabel_Dikeluarkan[1]),
                    paste0("Iterasi 2: keluarkan ", hasil_iter2$Variabel_Dikeluarkan[1]),
                    paste0("Iterasi 3: keluarkan ", hasil_iter3$Variabel_Dikeluarkan[1])),
  Fitur_Tersisa = c(length(semua_prediktor),
                    length(semua_prediktor) - 1,
                    length(semua_prediktor) - 2,
                    length(semua_prediktor) - 3),
  Max_VIF       = c(round(max(vif_base), 3),
                    hasil_iter1$Max_VIF[1],
                    hasil_iter2$Max_VIF[1],
                    hasil_iter3$Max_VIF[1]),
  N_VIF_gt10    = c(sum(vif_base > 10),
                    hasil_iter1$N_VIF_gt10[1],
                    hasil_iter2$N_VIF_gt10[1],
                    hasil_iter3$N_VIF_gt10[1]),
  Status        = c("❌ Multikolinear berat",
                    ifelse(hasil_iter1$N_VIF_gt10[1] == 0, "✅ Bersih", "⚠️ Masih ada"),
                    ifelse(hasil_iter2$N_VIF_gt10[1] == 0, "✅ Bersih", "⚠️ Masih ada"),
                    ifelse(hasil_iter3$N_VIF_gt10[1] == 0, "✅ Bersih", "⚠️ Masih ada")),
  stringsAsFactors = FALSE
)

optimal_idx <- which(ringkasan$N_VIF_gt10 == 0 &
                     ringkasan$Fitur_Tersisa ==
                       max(ringkasan$Fitur_Tersisa[ringkasan$N_VIF_gt10 == 0]))
if (length(optimal_idx) == 0) optimal_idx <- which.min(ringkasan$Max_VIF)

kable(ringkasan,
      caption = "Ringkasan Komparatif Skenario Pengeluaran Variabel",
      align = "lrrrc", booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, background = "#37474F", color = "white", bold = TRUE) %>%
  row_spec(1, background = "#FFEBEE") %>%
  row_spec(optimal_idx, background = "#C8E6C9", bold = TRUE)

# ══════════════════════════════════════════════════════════════════════════════
# [5] VIF FINAL & KORELASI
# ══════════════════════════════════════════════════════════════════════════════
n_drop_optimal <- length(semua_prediktor) - ringkasan$Fitur_Tersisa[optimal_idx]

vars_optimal_drop <- switch(as.character(n_drop_optimal),
  "1" = hasil_iter1$Variabel_Dikeluarkan[1],
  "2" = unlist(strsplit(hasil_iter2$Variabel_Dikeluarkan[1], " \\+ ")),
  "3" = unlist(strsplit(hasil_iter3$Variabel_Dikeluarkan[1], " \\+ "))
)

pred_final <- setdiff(semua_prediktor, vars_optimal_drop)
vif_final  <- hitung_vif(fitur_model, pred_final)

render_vif_table(vif_final,
  paste("VIF Final — Variabel dikeluarkan:",
        paste(vars_optimal_drop, collapse = ", ")))

# ── Matriks korelasi (pasangan r > 0.9) ──────────────────────────────────────
cor_mat  <- cor(fitur_model %>%
                  select(usia_bulan, bb_lahir, tb_sekarang, bb_sekarang,
                         pendapatan, diare_3bln, kunjungan_posyandu,
                         bmi, rasio_bb),
                use = "complete.obs")

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

if (nrow(cor_tinggi) == 0) {
  cat("\n✅ Tidak ada pasangan variabel dengan korelasi r > 0.9.",
      "\n   Multikolinearitas bersifat multivariat, bukan berpasangan.\n")
} else {
  kable(cor_tinggi,
        caption = "Pasangan Variabel dengan Korelasi r > 0.9",
        col.names = c("Var 1", "Var 2", "Korelasi (r)"),
        align = "llr", digits = 3, booktabs = TRUE) %>%
    kable_styling(bootstrap_options = c("striped", "hover"),
                  full_width = FALSE, position = "center") %>%
    row_spec(0, background = "#37474F", color = "white", bold = TRUE)
}

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("=== Dataset Final untuk Modeling ===\n")

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

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

cat("Missing  :", sum(is.na(data_stunting_model)), "\n")

write_csv(data_stunting_model, "C:/Users/User/Downloads/data_stunting_model.csv")
cat("\n✓ data_stunting_model.csv tersimpan!\n")

library(caret)

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

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

cat("Distribusi target:\n")

print(prop.table(table(data_stunting_model$status_stunting)) %>% round(3))

# ── Skenario 1: 80:20 Split ───────────────────────────────────────────────────
set.seed(42)
split_idx <- createDataPartition(data_stunting_model$status_stunting, p = 0.8, list = FALSE)
train_set <- data_stunting_model[split_idx,]
test_set  <- data_stunting_model[-split_idx,]

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

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

cat("Distribusi Training:\n")

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

cat("Distribusi Test:\n")

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

# ── Skenario 2: 5-Fold CV ────────────────────────────────────────────────────
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")

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

ctrl_cv <- trainControl(method = "cv", number = 5, classProbs = TRUE,
                        summaryFunction = twoClassSummary,
                        savePredictions = "final", sampling = NULL)
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("\n✓ 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 ~ . - bmi - rasio_bb,
                data   = train_set,
                family = binomial(link = "logit"))

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

summary(model_lr)

# ── Cek Perfect Separation ────────────────────────────────────────────────────
cat("\n=== CEK PERFECT SEPARATION ===\n")

cat(sprintf("  Null deviance    : %.4f (df = %d)\n",
            model_lr$null.deviance, model_lr$df.null))

cat(sprintf("  Residual deviance: %.6e (df = %d)\n",
            model_lr$deviance, model_lr$df.residual))

cat(sprintf("  Iterasi Fisher   : %d\n", model_lr$iter))

if (model_lr$deviance < 1e-3) {
  cat("  ⚠️  Perfect/quasi-complete separation terdeteksi!\n")
  cat("     SE sangat besar → p-value tidak reliable.\n")
  cat("     CI dihitung via Wald (±1.96×SE) sebagai aproksimasi.\n")
}

# ── Koefisien + p-value, diurutkan dari Pr(>|z|) terkecil ────────────────────
cat("\n=== KOEFISIEN MODEL (Diurutkan: p-value terkecil → terbesar) ===\n")

smry      <- summary(model_lr)$coefficients
coef_df   <- as.data.frame(smry) %>%
  rownames_to_column("Variabel") %>%
  rename(Estimate   = Estimate,
         Std_Error  = `Std. Error`,
         z_value    = `z value`,
         p_value    = `Pr(>|z|)`) %>%
  arrange(p_value) %>%                      # ← urutkan terkecil
  mutate(
    Signifikan = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01  ~ "**",
      p_value < 0.05  ~ "*",
      p_value < 0.1   ~ ".",
      TRUE            ~ ""
    )
  )

# Tampilkan semua, highlight top-3
kable(coef_df %>%
        mutate(across(c(Estimate, Std_Error, z_value), ~round(., 4)),
               p_value = formatC(p_value, format = "e", digits = 4)),
      caption = "Koefisien LR — Diurutkan Berdasarkan p-value (Terkecil ke Terbesar)",
      align   = "lrrrrc") %>%
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = FALSE, position = "center") %>%
  row_spec(1:3, background = "#C8E6C9", bold = TRUE) %>%       # top-3 hijau
  row_spec(which(coef_df$p_value < 0.05 & !((1:nrow(coef_df)) %in% 1:3)),
           background = "#FFF9C4") %>%                         # sig lain kuning
  add_header_above(c(" " = 1, "Estimasi" = 3, "Uji" = 2))

# ── Cetak Top-3 Paling Signifikan ─────────────────────────────────────────────
cat("\n=== TOP 3 VARIABEL PALING SIGNIFIKAN (p-value terkecil) ===\n")

top3 <- coef_df %>% slice(1:3)
for (i in 1:3) {
  cat(sprintf("  [%d] %-25s | p = %s | z = %.4f | Estimate = %.4f\n",
              i,
              top3$Variabel[i],
              formatC(top3$p_value[i], format = "e", digits = 4),
              top3$z_value[i],
              top3$Estimate[i]))
}

# ── Odds Ratio dengan CI Wald (robust terhadap separation) ───────────────────
cat("\n=== ODDS RATIO + 95% CI Wald ===\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(coef(model_lr) - 1.96 * sqrt(diag(vcov(model_lr)))), 4),
  CI_Upper  = round(exp(coef(model_lr) + 1.96 * sqrt(diag(vcov(model_lr)))), 4)
) %>%
  left_join(coef_df %>% select(Variabel, p_value), by = "Variabel") %>%
  arrange(p_value)                                   # ← ikut urutan signifikansi

kable(or_df %>% select(-p_value),
      caption = "Odds Ratio + 95% CI Wald (diurutkan: paling signifikan)",
      align   = "lrrrr") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(1:3, background = "#C8E6C9", bold = TRUE)

library(rpart); library(rpart.plot)

set.seed(42)
ctrl_dt   <- trainControl(method = "cv", number = 5, classProbs = TRUE,
                           summaryFunction = twoClassSummary)
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("=== HASIL TUNING cp ===\n")

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

cat("\ncp Optimal:", model_dt_cv$bestTune$cp, "\n")

res_dt    <- model_dt_cv$results
best_roc  <- max(res_dt$ROC)
se_roc    <- sd(res_dt$ROC) / sqrt(5)
cp_1se    <- res_dt$cp[res_dt$ROC >= (best_roc - se_roc)][1]
cat(sprintf("Best ROC: %.4f | 1-SE threshold: %.4f | cp (1-SE): %s\n",
            best_roc, best_roc - se_roc, cp_1se))

set.seed(42)
model_dt_final <- rpart(status_stunting ~ ., data = train_set, method = "class",
                         control = rpart.control(cp = cp_1se))
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 = "gray75",
           col = "#1A237E", main = paste0("Decision Tree Prediksi Stunting  (cp = ", cp_1se, ")"),
           cex = 0.78, box.palette = c(COL_NORMAL, COL_STUNTING))

library(randomForest)

set.seed(42)
ctrl_rf <- trainControl(method = "cv", number = 5, classProbs = TRUE,
                         summaryFunction = twoClassSummary)
mtry_grid   <- expand.grid(mtry = 2:8)
model_rf_cv <- train(status_stunting ~ ., data = train_set, method = "rf",
                     trControl = ctrl_rf, tuneGrid = mtry_grid,
                     metric = "ROC", ntree = 200)

cat("=== HASIL TUNING mtry ===\n")

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

cat("\nmtry Optimal:", model_rf_cv$bestTune$mtry, "\n")

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

imp_df <- varImp(model_rf_cv)$importance %>%
  as.data.frame() %>%
  tibble::rownames_to_column("Fitur") %>%
  rename(Importance = Overall) %>%
  arrange(Importance)

ggplot(imp_df, aes(x = reorder(Fitur, Importance), y = Importance, fill = Importance)) +
  geom_col(colour = "white", width = 0.72) +
  geom_text(aes(label = round(Importance, 1)), hjust = -0.1,
            size = 3.2, colour = "#37474F") +
  scale_fill_gradient(low = "#BBDEFB", high = COL_NORMAL) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  labs(title    = "Feature Importance — Random Forest",
       subtitle = paste0("mtry optimal = ", model_rf_cv$bestTune$mtry, " | ntree = 200"),
       x = NULL, y = "Importance Score (ROC)",
       caption  = "Semakin tinggi skor → semakin penting fitur dalam prediksi") +
  theme(legend.position = "none")

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))

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

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=== RF TANPA ROSE ===\n"); print(cm_no)