# ── Manipulasi Data ───────────────────────────────────────────────
library(foreign)           # read.dbf
library(dplyr)
library(readxl)            # read_excel
library(tidyr)
library(htmltools)         # HTML

# ── Pemodelan ─────────────────────────────────────────────────────
library(ResourceSelection) # Hosmer-Lemeshow
library(logistf)           # Firth logistic regression
library(pROC)              # ROC-AUC

# ── Output Tabel ──────────────────────────────────────────────────
library(knitr)
library(kableExtra)

# Paket yang dipanggil dengan prefix (themis::smoten, car::vif)
for (pkg in c("themis", "car")) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    cat("PERINGATAN: paket", pkg, "belum terpasang. Jalankan install.packages('", pkg, "').\n", sep = "")
  }
}

data_susenas <- read.dbf("D:/Altytan Semester 7/Skripsi Rampunggg/Data Pengolahan/46654/Susenas KOR maret 2023/bv2_2026_02_22_08_59_08_ssn202303kor_ind_1.dbf")

# Konversi aman ke numerik (berlaku untuk kolom factor maupun character)
num <- function(x) as.numeric(as.character(x))

UMUR_MIN  <- 7
UMUR_MAX  <- 12
PROV_JAWA <- 31:36
allowed_vals <- c(1, 2, 3, 5, 6, 7)

dataku_filtered <- data_susenas %>%
  mutate(
    across(R1002:R1009, num),

    # Kolom kunci dijadikan numerik agar aman dipakai untuk join
    R101  = num(R101),
    R102  = num(R102),
    URUT  = num(URUT),
    umur  = num(R407),
    R610n = num(R610),

    JUMLAH_JENIS_DISABILITAS = rowSums(
      across(R1002:R1009, ~ .x %in% allowed_vals),
      na.rm = TRUE
    ),

    prov = R101,
    kabu = R101 * 100 + R102,

    # Y: Partisipasi Sekolah (1 = Bersekolah, 0 = Tidak Bersekolah)
    status = ifelse(
      !is.na(umur) & umur >= UMUR_MIN & umur <= UMUR_MAX & JUMLAH_JENIS_DISABILITAS > 0,
      ifelse(R610n == 2, 1, 0),
      NA_real_
    )
  ) %>%
  filter(
    !is.na(umur), umur >= UMUR_MIN, umur <= UMUR_MAX,
    JUMLAH_JENIS_DISABILITAS > 0,
    R101 %in% PROV_JAWA
  )

cat("Jumlah observasi setelah filter (SD, Jawa):", nrow(dataku_filtered), "\n")

kunci_rt <- c("R101", "R102", "URUT")

cek_kunci_unik <- function(df, kunci, nama) {
  n_baris <- nrow(df)
  n_unik  <- nrow(dplyr::distinct(df, dplyr::across(dplyr::all_of(kunci))))
  cat(sprintf("[Cek kunci] %-12s : %d baris, %d kunci unik (%s)\n",
              nama, n_baris, n_unik, paste(kunci, collapse = " + ")))
  if (n_baris != n_unik) {
    stop(sprintf(
      "Kunci %s tidak unik pada %s (%d baris vs %d kunci unik). Join akan menggandakan baris. Periksa data sumber.",
      paste(kunci, collapse = " + "), nama, n_baris, n_unik
    ))
  }
  invisible(TRUE)
}

join_aman <- function(x, y, by, nama) {
  n_sebelum <- nrow(x)
  out <- dplyr::left_join(x, y, by = by)
  cat(sprintf("[Cek join]  %-22s : %d baris sebelum, %d baris sesudah\n",
              nama, n_sebelum, nrow(out)))
  if (nrow(out) != n_sebelum) {
    stop(sprintf("Join '%s' mengubah jumlah baris (%d -> %d). Kemungkinan many-to-many; periksa kunci join.",
                 nama, n_sebelum, nrow(out)))
  }
  out
}

# ==================================================================
# FUNGSI REKODE IJAZAH (R614) — satu sumber kebenaran (sama dengan skripsi)
# ==================================================================
# Biner : 0 = ≤ SMP/Sederajat (ref) ; 1 = Minimal SMA/Sederajat
ijazah_biner <- function(x) {
  x <- as.numeric(x)
  dplyr::case_when(
    x >= 11 & x <= 24   ~ 1L,
    x %in% c(0:10, 25)  ~ 0L,
    TRUE                ~ NA_integer_
  )
}

data_krt <- data_susenas %>%
  filter(num(R403) == 1) %>%
  mutate(
    R101 = num(R101),
    R102 = num(R102),
    URUT = num(URUT),
    R405 = num(R405),
    R614 = num(R614),
    R705 = num(R705),
    R706 = num(R706),

    # X4: Jenis Kelamin KRT
    X4 = case_when(
      R405 == 1 ~ 1L,
      R405 == 2 ~ 0L,
      TRUE      ~ NA_integer_
    ),

    # X5: Tingkat Pendidikan KRT (biner) — dari ijazah R614
    X5 = ijazah_biner(R614),

    # X6: Pekerjaan KRT
    # 0 = Tidak bekerja; 1 = Pertanian; 2 = Nonpertanian
    X6 = case_when(
      R706 >= 1 & R706 <= 6   ~ 1L,
      R706 >= 7 & R706 <= 26  ~ 2L,
      R705 == 5               ~ 0L,
      TRUE                    ~ NA_integer_
    )
  ) %>%
  select(all_of(kunci_rt), X4, X5, X6)

cek_kunci_unik(data_krt, kunci_rt, "data_krt")

files_43 <- c(
  "D:/Altytan Semester 7/Skripsi Rampunggg/Data Pengolahan/46654/Susenas KP Maret 2023/bv2_2026_02_22_08_55_02_ssn202303kp_blok43_51_94.dbf",
  "D:/Altytan Semester 7/Skripsi Rampunggg/Data Pengolahan/46654/Susenas KP Maret 2023/bv2_2026_02_22_08_49_01_ssn202303kp_blok43_11_31.dbf",
  "D:/Altytan Semester 7/Skripsi Rampunggg/Data Pengolahan/46654/Susenas KP Maret 2023/bv2_2026_02_22_08_52_06_ssn202303kp_blok43_32_36.dbf"
)

data_kp_43 <- bind_rows(lapply(files_43, read.dbf, as.is = TRUE)) %>%
  mutate(
    R101 = num(R101),
    R102 = num(R102),
    R105 = num(R105),
    URUT = num(URUT)
  )

prov_hilang <- setdiff(PROV_JAWA, unique(data_kp_43$R101))
if (length(prov_hilang) > 0) {
  cat("PERINGATAN: provinsi berikut tidak ditemukan di data KP blok 4.3:",
      paste(prov_hilang, collapse = ", "), "\n")
} else {
  cat("Seluruh provinsi Jawa (31–36) ditemukan di data KP blok 4.3.\n")
}

GK <- read_excel(
  "D:/Altytan Semester 7/Skripsi Rampunggg/Data Pengolahan/Garis Kemiskinan (Rupiah_Kapita_Bulan) Menurut Provinsi dan Daerah , 2023.xlsx",
  sheet = "Sheet2"
) %>%
  mutate(R101 = num(R101), R105 = num(R105))

data_kp_x9 <- join_aman(data_kp_43, GK, by = c("R101", "R105"), nama = "KP blok 4.3 x GK") %>%
  mutate(
    # X9: Status Kemiskinan RT (0 = Miskin/ref; 1 = Tidak Miskin)
    X9 = if_else(KAPITA < GK, 0L, 1L, missing = NA_integer_)
  ) %>%
  select(all_of(kunci_rt), X9)

cek_kunci_unik(data_kp_x9, kunci_rt, "data_kp_x9")

dataku_final <- dataku_filtered %>%
  mutate(
    R105 = num(R105),
    R405 = num(R405),
    R615 = num(R615),
    R616 = num(R616),

    # Y: Partisipasi Sekolah
    Y  = as.integer(status),

    # X1: Jenis Kelamin Anak
    X1 = case_when(
      R405 == 1 ~ 1L,
      R405 == 2 ~ 0L,
      TRUE      ~ NA_integer_
    ),

    # X2: Jenis Disabilitas (0 = Ganda/ref; 1 = Tunggal)
    X2 = case_when(
      JUMLAH_JENIS_DISABILITAS > 1  ~ 0L,
      JUMLAH_JENIS_DISABILITAS == 1 ~ 1L,
      TRUE                          ~ NA_integer_
    ),

    # X3: Tingkat Kesulitan Disabilitas (prioritas terberat)
    X3 = case_when(
      rowSums(across(R1002:R1009, ~ .x %in% c(1, 5)), na.rm = TRUE) > 0 ~ 0L,
      rowSums(across(R1002:R1009, ~ .x %in% c(2, 6)), na.rm = TRUE) > 0 ~ 1L,
      rowSums(across(R1002:R1009, ~ .x %in% c(3, 7)), na.rm = TRUE) > 0 ~ 2L,
      TRUE ~ NA_integer_
    ),

    # X7: Daerah Tempat Tinggal
    X7 = case_when(
      R105 == 1 ~ 1L,
      R105 == 2 ~ 0L,
      TRUE      ~ NA_integer_
    ),

    # X8: Kepemilikan KIP/PIP
    X8 = if_else(
      R615 %in% c(1, 2) | R616 %in% c(1, 2),
      1L, 0L
    )
  ) %>%
  join_aman(data_krt,   by = kunci_rt, nama = "x data_krt (X4-X6)") %>%
  join_aman(data_kp_x9, by = kunci_rt, nama = "x data_kp_x9 (X9)") %>%
  select(prov, kabu, Y, X1, X2, X3, X4, X5, X6, X7, X8, X9)

dataku_final <- dataku_final %>%
  mutate(
    Y       = factor(Y,       levels = c(0, 1),    labels = c("Tidak bersekolah", "Bersekolah")),
    X1      = factor(X1,      levels = c(0, 1),    labels = c("Perempuan", "Laki-laki")),
    X2      = factor(X2,      levels = c(0, 1),    labels = c("Disabilitas ganda", "Disabilitas tunggal")),
    X3      = factor(X3,      levels = c(0, 1, 2), labels = c("Sangat kesulitan", "Kesulitan", "Sedikit kesulitan")),
    X4      = factor(X4,      levels = c(0, 1),    labels = c("Perempuan", "Laki-laki")),
    X5      = factor(X5,      levels = c(0, 1),    labels = c("≤ SMP", "≥ SMA")),
    X6      = factor(X6,      levels = c(0, 1, 2), labels = c("Tidak bekerja", "Pertanian", "Nonpertanian")),
    X7      = factor(X7,      levels = c(0, 1),    labels = c("Perdesaan", "Perkotaan")),
    X8      = factor(X8,      levels = c(0, 1),    labels = c("Tidak memiliki", "Memiliki")),
    X9      = factor(X9,      levels = c(0, 1),    labels = c("Miskin", "Tidak miskin"))
  )

data.frame(
  Kode     = c("Y", "X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9"),
  Nama     = c(
    "Partisipasi Sekolah",
    "Jenis Kelamin Anak",
    "Jenis Disabilitas",
    "Tingkat Kesulitan Disabilitas",
    "Jenis Kelamin KRT",
    "Pendidikan KRT (biner)",
    "Pekerjaan KRT",
    "Daerah Tempat Tinggal",
    "Kepemilikan KIP/PIP",
    "Status Kemiskinan RT"
  ),
  Sumber = c(
    "R610",
    "R405", "R1002–R1009", "R1002–R1009",
    "R403 == 1 & R405",
    "R403 == 1 & R614",
    "R403 == 1 & R705/R706",
    "R105", "R615/R616", "KAPITA vs GK"
  ),
  Kategori = c(
    "Dependen",
    rep("Karakteristik Anak", 3),
    rep("Karakteristik KRT",  3),
    rep("Sosial-Ekonomi RT",  3)
  ),
  Koding = c(
    "1 = Bersekolah; 0 = Tidak bersekolah",
    "0 = Perempuan*; 1 = Laki-laki",
    "0 = Disabilitas ganda*; 1 = Disabilitas tunggal",
    "0 = Sangat kesulitan*; 1 = Kesulitan; 2 = Sedikit kesulitan",
    "0 = Perempuan*; 1 = Laki-laki",
    "0 = ≤ SMP/Sederajat*; 1 = Minimal SMA/Sederajat",
    "0 = Tidak bekerja*; 1 = Pertanian; 2 = Nonpertanian",
    "0 = Perdesaan*; 1 = Perkotaan",
    "0 = Tidak memiliki*; 1 = Memiliki",
    "0 = Miskin*; 1 = Tidak miskin"
  ),
  check.names = FALSE
) %>%
  kable(caption = "Referensi Variabel Penelitian (Y, X1–X9)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE) %>%
  column_spec(1, bold = TRUE, width = "7em") %>%
  column_spec(4, italic = TRUE) %>%
  row_spec(1, background = "#fff3e0") %>%
  footnote(
    general = "* = kategori referensi.",
    general_title = "Keterangan: ",
    footnote_as_chunk = TRUE,
    escape = FALSE
  )

prov_lab <- data.frame(
  prov = 31:36,
  Provinsi = c("DKI Jakarta", "Jawa Barat", "Jawa Tengah", "DI Yogyakarta", "Jawa Timur", "Banten")
)

tab_prov <- dataku_final %>%
  filter(!is.na(Y)) %>%
  count(prov, Y) %>%
  pivot_wider(names_from = Y, values_from = n, values_fill = 0) %>%
  left_join(prov_lab, by = "prov") %>%
  mutate(
    Total = `Tidak bersekolah` + Bersekolah,
    `% Tidak bersekolah` = round(`Tidak bersekolah` / Total * 100, 2)
  ) %>%
  select(Kode = prov, Provinsi, `Tidak bersekolah`, Bersekolah, Total, `% Tidak bersekolah`) %>%
  arrange(Kode)

cat(kbl_h(tab_prov, "Partisipasi Sekolah per Provinsi di Jawa (sampel, tanpa pembobot)"))

sel_jarang <- bind_rows(lapply(xvars_global, function(v) {
  tb <- table(df_all[[v]], df_all$status_sekolah)
  data.frame(
    Variabel = LABEL_X[[v]],
    Kategori = rownames(tb),
    `Tidak bersekolah` = as.integer(tb[, "0"]),
    Bersekolah         = as.integer(tb[, "1"]),
    check.names = FALSE
  )
})) %>%
  mutate(Tanda = ifelse(pmin(`Tidak bersekolah`, Bersekolah) < 5, "⚠ < 5", ""))

tbl_sj <- knitr::kable(sel_jarang, format = "html", escape = FALSE,
                       caption = "Frekuensi kategori prediktor menurut Y") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13, position = "left")
if (any(sel_jarang$Tanda != "")) {
  tbl_sj <- row_spec(tbl_sj, which(sel_jarang$Tanda != ""), background = "#ffebee")
}
cat(as.character(tbl_sj))

cat('<p><em>', sum(sel_jarang$Tanda != ""), ' dari ', nrow(sel_jarang),
    ' kategori memiliki kurang dari 5 observasi pada salah satu kelas Y.</em></p>', sep = "")

# ── Wadah rekap odds ratio lintas kombinasi ───────────────────────
OR_RECAP <- data.frame()

catat_or <- function(kombinasi, metode, variabel, or, p, n, ci_l = NA_real_, ci_u = NA_real_) {
  OR_RECAP <<- bind_rows(OR_RECAP, data.frame(
    Kombinasi    = kombinasi,
    Metode       = metode,
    Variabel     = variabel,
    `Odds Ratio` = round(or, 3),
    `p-value`    = round(p, 4),
    `CI 2.5%`    = round(ci_l, 3),
    `CI 97.5%`   = round(ci_u, 3),
    n            = n,
    row.names    = NULL,
    check.names  = FALSE
  ))
}

# ── Wrapper HTML untuk tiap sub-bagian ───────────────────────────
sec <- function(icon, judul, isi) {
  paste0(
    '<div style="margin-bottom:18px">',
    '<h5 style="color:#1565c0;border-bottom:2px solid #bbdefb;',
    'padding-bottom:5px;margin-top:18px;font-size:13.5px;">',
    icon, ' ', judul, '</h5>', isi, '</div>'
  )
}

# ── EDA: distribusi Y ────────────────────────────────────────────
eda_html <- function(df) {
  tbl <- table(factor(df$status_sekolah, levels = c(0, 1)))
  pct <- round(prop.table(tbl) * 100, 2)
  kbl_h(data.frame(
    Kategori = c("Tidak Bersekolah (0)", "Bersekolah (1)"),
    N        = as.integer(c(tbl["0"], tbl["1"])),
    Persen   = paste0(c(pct["0"], pct["1"]), " %")
  ), "Distribusi Y — Partisipasi Sekolah")
}

# ── Uji Simultan (GLM-based) ─────────────────────────────────────
simultan_html <- function(model, model_null, xvars_model) {
  G2    <- model_null$deviance - model$deviance
  df_G2 <- model_null$df.residual - model$df.residual
  chi_k <- qchisq(0.95, df_G2)
  pval  <- pchisq(G2, df_G2, lower.tail = FALSE)
  kbl_h(data.frame(
    Statistik = c("-2LL Null", "-2LL Model", "G² Hitung", "df",
                  "χ² Tabel (α=5%)", "p-value", "Keputusan"),
    Nilai = c(
      round(model_null$deviance, 3), round(model$deviance, 3),
      round(G2, 3), df_G2, round(chi_k, 3), round(pval, 4),
      ifelse(G2 > chi_k, "✅ Tolak H0 — Model Signifikan", "❌ Gagal Tolak H0")
    )
  ), "Uji Simultan (Likelihood Ratio Test / G²)")
}

# ── Uji Parsial Wald (GLM-based) ────────────────────────────────
wald_html <- function(model, alpha = 0.05) {
  koef      <- summary(model)$coefficients
  wald_stat <- (koef[, 1] / koef[, 2])^2
  chi_tabel <- qchisq(1 - alpha, df = 1)

  kbl_h(data.frame(
    Variabel        = rownames(koef),
    B               = round(koef[, 1], 4),
    SE              = round(koef[, 2], 4),
    `Wald (Chi-Sq Hitung)` = round(wald_stat, 4),
    df              = 1,
    `Chi-Sq Tabel (df=1, α=5%)` = round(chi_tabel, 4),
    `p-value`       = round(koef[, 4], 4),
    Kesimpulan      = ifelse(wald_stat > chi_tabel & koef[, 4] < alpha,
                              "✅ Signifikan", "❌ Tdk Signifikan"),
    check.names = FALSE
  ), "Uji Parsial (Wald Test — Statistik Hitung dibandingkan Chi-Square Tabel)")
}

# ── Hosmer-Lemeshow ──────────────────────────────────────────────
hl_html <- function(model, y) {
  hl <- tryCatch(hoslem.test(y, fitted(model), g = 10), error = function(e) NULL)
  if (is.null(hl)) return("<p><em>⚠️ Hosmer-Lemeshow tidak dapat dihitung.</em></p>")
  kbl_h(data.frame(
    Statistik = c("χ² HL", "df", "p-value", "Keputusan"),
    Nilai = c(
      round(hl$statistic, 3), hl$parameter, round(hl$p.value, 4),
      ifelse(hl$p.value > 0.05,
             "✅ Gagal Tolak H0 — Model FIT",
             "❌ Tolak H0 — TIDAK FIT")
    )
  ), "Goodness of Fit (Hosmer-Lemeshow)")
}

hl_p <- function(y, p) {
  hl <- tryCatch(hoslem.test(y, p, g = 10), error = function(e) NULL)
  if (is.null(hl)) NA_real_ else as.numeric(hl$p.value)
}

g2_p_glm <- function(model) {
  tryCatch(
    pchisq(model$null.deviance - model$deviance,
           model$df.null - model$df.residual, lower.tail = FALSE),
    error = function(e) NA_real_
  )
}

# ── Odds Ratio (GLM-based) dengan CI 95% ─────────────────────────
or_html <- function(model) {
  or <- exp(coef(model))
  ci <- tryCatch(exp(confint.default(model)), error = function(e) NULL)
  ci_l <- if (!is.null(ci)) ci[, 1] else rep(NA_real_, length(or))
  ci_u <- if (!is.null(ci)) ci[, 2] else rep(NA_real_, length(or))
  kbl_h(data.frame(
    Variabel     = names(or),
    `exp(B)`     = round(or, 4),
    `CI 2.5%`    = round(ci_l, 4),
    `CI 97.5%`   = round(ci_u, 4),
    Interpretasi = ifelse(names(or) == "(Intercept)", "Baseline odds",
                   ifelse(or > 1,
                          "↑ Peluang bersekolah lebih besar dibanding referensi",
                          "↓ Peluang bersekolah lebih kecil dibanding referensi")),
    check.names = FALSE
  ), "Interpretasi Parameter (Odds Ratio, CI 95%)")
}

# ── Keputusan eksplisit variabel OR≈1 / tidak signifikan ─────────
keputusan_html_generic <- function(nama, or, pval) {
  pv <- ifelse(is.na(pval), 1, pval)
  df_dec <- data.frame(
    Variabel  = nama,
    `exp(B)`  = round(or, 3),
    `p-value` = round(pv, 4),
    check.names = FALSE
  ) %>% filter(Variabel != "(Intercept)")

  df_dec$Kandidat_Ditinjau <- ifelse(
    (df_dec$`exp(B)` >= 0.8 & df_dec$`exp(B)` <= 1.25) | df_dec$`p-value` >= 0.05,
    "⚠️ OR mendekati 1 dan/atau tidak signifikan", "—")

  df_dec$Keputusan <- ifelse(
    df_dec$Kandidat_Ditinjau != "—",
    "Dipertahankan dalam model (justifikasi teoretis, bukan dibuang otomatis)",
    "Dipertahankan (signifikan/berkontribusi)"
  )

  paste0(
    kbl_h(df_dec, "Keputusan Eksplisit: Variabel dengan OR \u2248 1 atau Tidak Signifikan"),
    '<p><em>Catatan: variabel yang ditandai tidak otomatis dikeluarkan dari model. ',
    'Keputusan mempertahankan atau mengeluarkan suatu prediktor perlu disertai justifikasi ',
    'teoretis (relevansi variabel terhadap partisipasi sekolah anak disabilitas), bukan semata-mata ',
    'berdasarkan signifikansi statistik, karena model regresi logistik di sini dipakai untuk tujuan ',
    'eksplanatori, bukan seleksi prediktif murni. Peneliti perlu menuliskan alasan spesifik untuk ',
    'tiap variabel yang ditandai pada baris di atas sebelum finalisasi model.</em></p>'
  )
}

keputusan_variabel_html <- function(model) {
  koef <- summary(model)$coefficients
  keputusan_html_generic(rownames(koef), exp(coef(model)[rownames(koef)]), koef[, 4])
}

# ── Kriteria Informasi: AIC / AICc / BIC ─────────────────────────
ic_html <- function(model, n, k, method = "glm") {
  if (method == "glm") {
    aic_val <- tryCatch(AIC(model),                error = function(e) NA)
    bic_val <- tryCatch(BIC(model),                error = function(e) NA)
    ll_val  <- tryCatch(as.numeric(logLik(model)), error = function(e) NA)
  } else if (method == "firth") {
    ll_val  <- max(model$loglik)
    aic_val <- -2 * ll_val + 2 * k
    bic_val <- -2 * ll_val + log(n) * k
  } else {
    aic_val <- NA; bic_val <- NA; ll_val <- NA
  }

  aicc_val <- if (!is.na(aic_val) && (n - k - 1) > 0) {
    aic_val + (2 * k^2 + 2 * k) / (n - k - 1)
  } else NA

  rows <- paste0(
    '<tr><td><strong>AIC</strong></td><td>',  round(aic_val,  3),
    '</td><td>Semakin kecil semakin baik</td></tr>',
    '<tr><td><strong>AICc</strong></td><td>', round(aicc_val, 3),
    '</td><td>Koreksi utk sampel kecil</td></tr>',
    '<tr><td><strong>BIC</strong></td><td>',  round(bic_val,  3),
    '</td><td>Penalti lebih besar pada parameter</td></tr>',
    '<tr><td>Log-Likelihood</td><td>',         round(ll_val,  3),
    '</td><td>', if (method == "firth") "Penalized log-likelihood" else "", '</td></tr>'
  )

  paste0(
    '<table class="table table-striped table-condensed" style="font-size:13px;width:auto;margin-bottom:0">',
    '<thead><tr>',
    '<th style="min-width:130px">Kriteria</th>',
    '<th style="min-width:90px">Nilai</th>',
    '<th>Keterangan</th></tr></thead><tbody>',
    rows,
    '<tr><td>k (parameter)</td><td>', k, '</td><td>Termasuk intercept</td></tr>',
    '<tr><td>n (observasi)</td><td>', n, '</td><td></td></tr>',
    '</tbody></table>'
  )
}

# ── A. Perfect Separation Check ──────────────────────────────────
separation_html <- function(coefs, ses, is_firth = FALSE) {
  if (is_firth) {
    return(paste0(
      '<div class="eval-sub">',
      '<h6>Perfect Separation Check</h6>',
      '<p class="sep-ok">&#10003; Firth (Penalized Likelihood) dirancang khusus ',
      'menangani complete/quasi-complete separation — hasil tetap valid tanpa bias.</p>',
      '</div>'
    ))
  }
  flag_b  <- !is.na(coefs) & abs(coefs) > 10
  flag_se <- !is.na(ses) & ses > 5
  flag    <- flag_b | flag_se
  status  <- if (any(flag, na.rm = TRUE)) {
    '<p class="sep-flag">&#9888; Terdeteksi potensi separation pada variabel berikut:</p>'
  } else {
    '<p class="sep-ok">&#10003; Tidak terdeteksi perfect/quasi-complete separation.</p>'
  }
  tbl <- kbl_h(data.frame(
    Variabel   = names(coefs),
    `|B|>10`   = ifelse(flag_b,  "&#9888; Ya", "&#10003; Tidak"),
    `SE>5`     = ifelse(flag_se, "&#9888; Ya", "&#10003; Tidak"),
    Status     = ifelse(flag,
                        '<span class="sep-flag">&#9888; Potensi Separation</span>',
                        '<span class="sep-ok">&#10003; Aman</span>'),
    check.names = FALSE
  ), "Perfect Separation Check")
  paste0('<div class="eval-sub"><h6>Perfect Separation Check</h6>', status, tbl, '</div>')
}

# ── B. Confusion Matrix + Metrics (blok seluruh data) ────────────
conf_metrics <- function(y_true, y_prob, cutoff) {
  y_pred <- as.integer(y_prob >= cutoff)
  tp <- sum(y_true == 1 & y_pred == 1)
  tn <- sum(y_true == 0 & y_pred == 0)
  fp <- sum(y_true == 0 & y_pred == 1)
  fn <- sum(y_true == 1 & y_pred == 0)
  acc  <- ifelse((tp + tn + fp + fn) > 0, (tp + tn) / (tp + tn + fp + fn), NA)
  sens <- ifelse((tp + fn) > 0, tp / (tp + fn), NA)
  spec <- ifelse((tn + fp) > 0, tn / (tn + fp), NA)
  prec <- ifelse((tp + fp) > 0, tp / (tp + fp), NA)
  f1   <- ifelse(!is.na(prec) & !is.na(sens) & (prec + sens) > 0,
                 2 * prec * sens / (prec + sens), NA)
  list(tp = tp, tn = tn, fp = fp, fn = fn,
       acc = acc, sens = sens, spec = spec, prec = prec, f1 = f1)
}

conf_block_html <- function(y_true, y_prob, cutoff, label_class, label_text) {
  m   <- conf_metrics(y_true, y_prob, cutoff)
  pct <- function(x) ifelse(is.na(x), "NA", paste0(round(x * 100, 2), " %"))
  cm_tbl <- kbl_h(data.frame(
    ` `          = c("Pred: Tidak Bersekolah (0)", "Pred: Bersekolah (1)"),
    `Aktual: 0`  = c(m$tn, m$fp),
    `Aktual: 1`  = c(m$fn, m$tp),
    check.names  = FALSE
  ), paste0("Confusion Matrix — ", label_text))
  chips <- paste0(
    '<div class="metric-row">',
    '<span class="metric-chip">Accuracy: ',    pct(m$acc),  '</span>',
    '<span class="metric-chip">Sensitivity: ', pct(m$sens), '</span>',
    '<span class="metric-chip">Specificity: ', pct(m$spec), '</span>',
    '<span class="metric-chip">Precision: ',   pct(m$prec), '</span>',
    '<span class="metric-chip">F1-Score: ',    pct(m$f1),   '</span>',
    '</div>'
  )
  paste0(
    '<div class="cm-box">',
    '<span class="cm-label ', label_class, '">', label_text, '</span>',
    cm_tbl, chips,
    '</div>'
  )
}

cm_html <- function(y_true, y_prob) {
  roc_obj <- tryCatch(
    pROC::roc(y_true, y_prob, levels = c(0, 1), direction = "<", quiet = TRUE),
    error = function(e) NULL
  )
  cut_you <- if (!is.null(roc_obj)) {
    co <- tryCatch(
      pROC::coords(roc_obj, "best", best.method = "youden", ret = "threshold"),
      error = function(e) NULL
    )
    v <- if (!is.null(co)) suppressWarnings(as.numeric(unlist(co)[1])) else NA_real_
    if (is.finite(v)) v else 0.5
  } else 0.5

  b50 <- conf_block_html(y_true, y_prob, 0.5,     "cut50",  "Cutoff 0.5")
  byo <- conf_block_html(y_true, y_prob, cut_you, "cutyou",
                         paste0("Youden Optimal (", round(cut_you, 3), ")"))
  paste0('<div class="eval-sub"><h6>Confusion Matrix &amp; Metrik Klasifikasi</h6>',
         '<div class="cm-wrap">', b50, byo, '</div></div>')
}

# ── C. ROC-AUC + Toggle Plot ─────────────────────────────────────
roc_html <- function(y_true, y_prob, plot_id) {
  roc_obj <- tryCatch(
    pROC::roc(y_true, y_prob, levels = c(0, 1), direction = "<", quiet = TRUE),
    error = function(e) NULL
  )
  if (is.null(roc_obj))
    return('<div class="eval-sub"><h6>ROC-AUC</h6><p><em>ROC tidak dapat dihitung.</em></p></div>')

  auc_val <- as.numeric(pROC::auc(roc_obj))
  auc_cat <- if (auc_val >= 0.9) list(cls = "auc-great", txt = "Sangat Baik")
        else if (auc_val >= 0.8) list(cls = "auc-good",  txt = "Baik")
        else if (auc_val >= 0.7) list(cls = "auc-fair",  txt = "Cukup")
        else                     list(cls = "auc-poor",  txt = "Kurang")

  tmp <- tempfile(fileext = ".png")
  png(tmp, width = 400, height = 400, res = 90)
  plot(roc_obj, col = "#1565c0", lwd = 2.5,
       main = paste0("ROC Curve\nAUC = ", round(auc_val, 4)),
       cex.main = 0.9, cex.axis = 0.8, cex.lab = 0.85)
  abline(a = 0, b = 1, lty = 2, col = "gray60")
  dev.off()
  b64 <- knitr::image_uri(tmp)
  unlink(tmp)

  pid <- gsub("[^a-zA-Z0-9]", "_", plot_id)
  paste0(
    '<div class="eval-sub"><h6>ROC-AUC</h6>',
    '<span class="auc-badge ', auc_cat$cls, '">AUC = ', round(auc_val, 4), '</span>',
    '<strong>', auc_cat$txt, '</strong>',
    ' &nbsp;|&nbsp; Skala: 0.9–1.0 Sangat Baik | 0.8–0.9 Baik | 0.7–0.8 Cukup | &lt;0.7 Kurang<br>',
    '<button class="btn-roc" onclick="toggleRoc(\'', pid, '\')">',
    '&#128200; Tampilkan / Sembunyikan Plot ROC</button>',
    '<div class="roc-plot-wrap" id="roc_', pid, '">',
    '<img src="', b64, '" alt="ROC Curve">',
    '</div></div>'
  )
}

# ── D. Full Evaluasi Section ─────────────────────────────────────
evaluasi_html <- function(model, y_true, y_prob, plot_id, is_firth = FALSE) {
  coefs <- coef(model)
  ses   <- tryCatch(sqrt(diag(vcov(model))), error = function(e) rep(NA, length(coefs)))
  if (length(ses) != length(coefs)) ses <- rep(NA_real_, length(coefs))
  paste0(
    '<div class="eval-wrap">',
    separation_html(coefs, ses, is_firth),
    cm_html(y_true, y_prob),
    roc_html(y_true, y_prob, plot_id),
    '</div>'
  )
}

# ==================================================================
# METRIK EVALUASI — event model tetap Y = 1 (Bersekolah)
# Kelas minoritas = Tidak Bersekolah (Y = 0), sehingga recall, precision
# dan F1 kelas 0 dilaporkan eksplisit bersama Balanced Accuracy dan G-mean.
# ==================================================================
smote_metrics_school <- function(y_true, prob_school, cutoff_school = 0.50) {
  y_true <- as.integer(y_true)
  y_pred <- as.integer(prob_school >= cutoff_school)

  tp <- sum(y_true == 1 & y_pred == 1)
  tn <- sum(y_true == 0 & y_pred == 0)
  fp <- sum(y_true == 0 & y_pred == 1)
  fn <- sum(y_true == 1 & y_pred == 0)

  acc <- ifelse((tp + tn + fp + fn) > 0,
                (tp + tn) / (tp + tn + fp + fn), NA_real_)

  # Kinerja kelas 1 = Bersekolah
  sens_school <- ifelse((tp + fn) > 0, tp / (tp + fn), NA_real_)
  prec_school <- ifelse((tp + fp) > 0, tp / (tp + fp), NA_real_)
  f1_school   <- ifelse(!is.na(prec_school) && !is.na(sens_school) &&
                        (prec_school + sens_school) > 0,
                        2 * prec_school * sens_school /
                          (prec_school + sens_school), NA_real_)

  # Kinerja kelas 0 = Tidak Bersekolah (minoritas)
  recall_noschool <- ifelse((tn + fp) > 0, tn / (tn + fp), NA_real_)
  prec_noschool   <- ifelse((tn + fn) > 0, tn / (tn + fn), NA_real_)
  f1_noschool     <- ifelse(!is.na(prec_noschool) && !is.na(recall_noschool) &&
                            (prec_noschool + recall_noschool) > 0,
                            2 * prec_noschool * recall_noschool /
                              (prec_noschool + recall_noschool), NA_real_)

  bal <- ifelse(!is.na(sens_school) && !is.na(recall_noschool),
                (sens_school + recall_noschool) / 2, NA_real_)
  gmn <- ifelse(!is.na(sens_school) && !is.na(recall_noschool),
                sqrt(sens_school * recall_noschool), NA_real_)

  roc_obj <- tryCatch(
    pROC::roc(response = y_true, predictor = prob_school,
              levels = c(0, 1), direction = "<", quiet = TRUE),
    error = function(e) NULL
  )
  auc_val <- if (!is.null(roc_obj)) as.numeric(pROC::auc(roc_obj)) else NA_real_

  list(
    tp = tp, tn = tn, fp = fp, fn = fn,
    acc = acc,
    sens_school = sens_school,
    spec_school = recall_noschool,
    prec_school = prec_school,
    f1_school = f1_school,
    recall_noschool = recall_noschool,
    prec_noschool = prec_noschool,
    f1_noschool = f1_noschool,
    bal = bal,
    gmean = gmn,
    auc = auc_val,
    cutoff = cutoff_school
  )
}

# ---- Cutoff Youden ditentukan HANYA dari training asli ----
smote_youden_from_train <- function(y_true, prob_school) {
  y_true <- as.integer(y_true)

  roc_obj <- tryCatch(
    pROC::roc(response = y_true, predictor = prob_school,
              levels = c(0, 1), direction = "<", quiet = TRUE),
    error = function(e) NULL
  )
  if (is.null(roc_obj)) return(0.50)

  co <- tryCatch(
    pROC::coords(roc_obj, "best", best.method = "youden",
                 ret = "threshold", transpose = FALSE),
    error = function(e) NULL
  )
  if (is.null(co)) return(0.50)

  cut <- suppressWarnings(as.numeric(unlist(co)[1]))
  if (!is.finite(cut)) 0.50 else cut
}

# ---- Baris ringkasan lintas hasil ----
mk_ringkas <- function(kelompok, metode, skema, n_tr, n_sint, pct_min,
                       m, hlp, g2p, pemenang = NA_character_) {
  data.frame(
    Kelompok            = kelompok,
    Metode              = metode,
    Skema               = skema,
    n_training          = n_tr,
    n_sintetis          = n_sint,
    pct_minoritas_akhir = pct_min,
    BalAcc              = m$bal * 100,
    Recall_Tdk          = m$recall_noschool * 100,
    Gmean               = m$gmean * 100,
    AUC                 = m$auc,
    HL_p                = hlp,
    G2_p                = g2p,
    Pemenang            = pemenang,
    stringsAsFactors    = FALSE
  )
}

# ---- Pemenang di antara metode (Balanced Accuracy > recall Tdk Bersekolah > AUC) ----
pilih_pemenang <- function(mlist) {
  nilai <- function(x, ganti = -1) ifelse(is.na(x), ganti, round(x, 8))
  tb <- data.frame(
    metode = names(mlist),
    bal  = unname(sapply(mlist, function(m) nilai(m$bal))),
    rec0 = unname(sapply(mlist, function(m) nilai(m$recall_noschool))),
    auc  = unname(sapply(mlist, function(m) nilai(m$auc))),
    stringsAsFactors = FALSE
  )
  tb  <- tb[order(-tb$bal, -tb$rec0, -tb$auc), ]
  top <- tb[tb$bal == tb$bal[1] & tb$rec0 == tb$rec0[1] & tb$auc == tb$auc[1], ]
  if (nrow(top) > 1) paste0("Setara: ", paste(top$metode, collapse = " = ")) else top$metode[1]
}

run_biasa <- function(data, nama_dataset, xvars = xvars_global) {
  p   <- prep_data(data, xvars)
  df  <- p$df; xm <- p$xvars_model; fm <- p$formula_mod

  if (length(xm) == 0) return(list(html = "<p>⚠️ Tidak ada variabel valid.</p>", ringkas = NULL))

  mdl  <- tryCatch(glm(fm, data = df, family = binomial), error = function(e) NULL)
  mdl0 <- tryCatch(glm(status_sekolah ~ 1, data = df, family = binomial), error = function(e) NULL)
  if (is.null(mdl)) return(list(html = "<p>⚠️ Model gagal di-fit.</p>", ringkas = NULL))

  info <- paste0('<p><strong>n valid:</strong> ', nrow(df), ' obs.',
    if (length(p$konstan) > 0) paste0(' | <em>Konstan (dikeluarkan): ',
                                       paste(p$konstan, collapse = ", "), '</em>') else '',
    ' | <em>Model dilatih pada seluruh data, tanpa pembagian latih-uji.</em></p>')

  k_biasa <- length(coef(mdl))
  n_biasa <- nrow(df)
  pid     <- paste0("biasa-", gsub("[^a-zA-Z0-9]", "_", nama_dataset))

  # ── Rekam OR (dengan CI) untuk rekap Lampiran ──
  koef_b <- summary(mdl)$coefficients
  ci_b   <- tryCatch(exp(confint.default(mdl)), error = function(e) NULL)
  nm_b   <- rownames(koef_b)
  catat_or(nama_dataset, "Reglog Biasa", nm_b, exp(coef(mdl)[nm_b]), koef_b[, 4], n_biasa,
           if (!is.null(ci_b)) ci_b[nm_b, 1] else NA_real_,
           if (!is.null(ci_b)) ci_b[nm_b, 2] else NA_real_)

  m_ins <- smote_metrics_school(df$status_sekolah, fitted(mdl), 0.50)
  ringkas <- mk_ringkas(
    "Seluruh data (in-sample)", "Reglog Biasa", "Seluruh data",
    n_biasa, 0, mean(df$status_sekolah == 0) * 100,
    m_ins, hl_p(df$status_sekolah, fitted(mdl)), g2_p_glm(mdl), "—"
  )

  html <- paste0(info,
    sec("📊", "EDA — Distribusi Y",               eda_html(df)),
    sec("✅", "Goodness of Fit (H-L Test)",        hl_html(mdl, df$status_sekolah)),
    sec("📐", "Uji Simultan (G²)",                if (!is.null(mdl0)) simultan_html(mdl, mdl0, xm) else "<p>⚠️ Model null gagal.</p>"),
    sec("📋", "Uji Parsial (Wald Test)",           wald_html(mdl)),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC)", ic_html(mdl, n_biasa, k_biasa, "glm")),
    sec("📈", "Interpretasi Parameter (OR, CI 95%)", or_html(mdl)),
    sec("🧭", "Keputusan Variabel OR≈1 / Tdk Signifikan", keputusan_variabel_html(mdl)),
    sec("🔍", "Evaluasi Model (in-sample, seluruh data)",
        evaluasi_html(mdl, df$status_sekolah, fitted(mdl), pid))
  )
  list(html = html, ringkas = ringkas)
}

# ---- Probabilitas Firth untuk data baru (dengan cadangan manual) ----
prob_firth <- function(mdl, newdata, xm) {
  out <- tryCatch(as.numeric(predict(mdl, newdata = newdata, type = "response")),
                  error = function(e) NULL)
  if (!is.null(out) && length(out) == nrow(newdata) && all(is.finite(out))) return(out)

  # Cadangan: hitung manual dari matriks desain
  b <- mdl$coefficients
  X <- model.matrix(reformulate(xm), data = newdata)
  X <- X[, names(b), drop = FALSE]
  as.numeric(plogis(X %*% b))
}

run_firth <- function(data, nama_dataset, xvars = xvars_global) {
  p  <- prep_data(data, xvars)
  df <- p$df; xm <- p$xvars_model; fm <- p$formula_mod

  if (length(xm) == 0) return(list(html = "<p>⚠️ Tidak ada variabel valid.</p>", ringkas = NULL))

  mdl <- tryCatch(logistf(fm, data = df), error = function(e) e)
  if (inherits(mdl, "error")) {
    return(list(html = paste0("<p>⚠️ Model Firth gagal di-fit: ", mdl$message, "</p>"), ringkas = NULL))
  }

  coefs <- mdl$coefficients
  k_f   <- length(coefs)
  n_f   <- nrow(df)

  # Uji Simultan — penalized LRT
  G2    <- 2 * abs(diff(as.numeric(mdl$loglik)))
  df_G2 <- k_f - 1
  chi_k <- qchisq(0.95, df_G2)
  pval  <- pchisq(G2, df_G2, lower.tail = FALSE)
  sim_tbl <- kbl_h(data.frame(
    Statistik = c("LogLik Null", "LogLik Model", "G² Hitung", "df",
                  "χ² Tabel (α=5%)", "p-value", "Keputusan"),
    Nilai = c(
      round(min(mdl$loglik), 3), round(max(mdl$loglik), 3),
      round(G2, 3), df_G2, round(chi_k, 3), round(pval, 4),
      ifelse(G2 > chi_k, "✅ Tolak H0 — Signifikan", "❌ Gagal Tolak H0")
    )
  ), "Uji Simultan (Penalized LRT)")

  # Koefisien, SE, p-value — disejajarkan secara posisional dengan koefisien
  se_raw <- sqrt(diag(mdl$var))
  se     <- setNames(se_raw[seq_along(coefs)], names(coefs))
  pvals  <- as.numeric(mdl$prob)
  if (length(pvals) != k_f) pvals <- rep(NA_real_, k_f)
  names(pvals) <- names(coefs)
  chi_sq <- (coefs / se)^2
  ci_l   <- as.numeric(mdl$ci.lower)
  ci_u   <- as.numeric(mdl$ci.upper)

  wald_tbl <- kbl_h(data.frame(
    Variabel   = names(coefs),
    B          = round(coefs, 4),
    SE         = round(se, 4),
    `Wald (B/SE)²` = round(chi_sq, 4),
    `p-value (Penalized LRT)` = round(pvals, 4),
    Kesimpulan = ifelse(!is.na(pvals) & pvals < 0.05, "✅ Signifikan", "❌ Tdk Signifikan"),
    check.names = FALSE
  ), "Uji Per Variabel (Penalized LRT)")

  # Odds Ratio dengan CI profile penalized likelihood
  or_tbl <- kbl_h(data.frame(
    Variabel     = names(coefs),
    `exp(B)`     = round(exp(coefs), 4),
    `CI 2.5%`    = round(exp(ci_l), 4),
    `CI 97.5%`   = round(exp(ci_u), 4),
    Interpretasi = ifelse(names(coefs) == "(Intercept)", "Baseline odds",
                   ifelse(exp(coefs) > 1,
                          "↑ Peluang bersekolah lebih besar dibanding referensi",
                          "↓ Peluang bersekolah lebih kecil dibanding referensi")),
    check.names = FALSE
  ), "Interpretasi Parameter (Odds Ratio — CI Firth/Profile Likelihood)")

  # Goodness of Fit Hosmer-Lemeshow — memakai mdl$predict
  prob_fit <- as.numeric(mdl$predict)
  hl_firth <- tryCatch(hoslem.test(df$status_sekolah, prob_fit, g = 10), error = function(e) NULL)
  hl_firth_tbl <- if (is.null(hl_firth)) {
    "<p><em>⚠️ Hosmer-Lemeshow tidak dapat dihitung.</em></p>"
  } else {
    kbl_h(data.frame(
      Statistik = c("χ² HL", "df", "p-value", "Keputusan"),
      Nilai = c(
        round(hl_firth$statistic, 3), hl_firth$parameter, round(hl_firth$p.value, 4),
        ifelse(hl_firth$p.value > 0.05,
               "✅ Gagal Tolak H0 — Model FIT",
               "❌ Tolak H0 — TIDAK FIT")
      )
    ), "Goodness of Fit (Hosmer-Lemeshow)")
  }

  info <- paste0('<p><strong>n valid:</strong> ', n_f, ' obs.',
    if (length(p$konstan) > 0) paste0(' | <em>Konstan (dikeluarkan): ',
                                       paste(p$konstan, collapse = ", "), '</em>') else '',
    ' | <em>Firth menangani quasi-complete separation &amp; kelas minoritas kecil. ',
    'Model dilatih pada seluruh data, tanpa pembagian latih-uji.</em></p>')

  pid <- paste0("firth-", gsub("[^a-zA-Z0-9]", "_", nama_dataset))

  # ── Rekam OR (dengan CI) untuk rekap Lampiran ──
  catat_or(nama_dataset, "Firth", names(coefs), exp(coefs), pvals, n_f, exp(ci_l), exp(ci_u))

  m_ins <- smote_metrics_school(df$status_sekolah, prob_fit, 0.50)
  ringkas <- mk_ringkas(
    "Seluruh data (in-sample)", "Firth", "Seluruh data",
    n_f, 0, mean(df$status_sekolah == 0) * 100,
    m_ins, hl_p(df$status_sekolah, prob_fit), pval, "—"
  )

  html <- paste0(info,
    sec("📊", "EDA — Distribusi Y",                    eda_html(df)),
    sec("✅", "Goodness of Fit (H-L Test)",            hl_firth_tbl),
    sec("📐", "Uji Simultan (Penalized LRT)",          sim_tbl),
    sec("📋", "Uji Per Variabel (Penalized LRT)",      wald_tbl),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC)",     ic_html(mdl, n_f, k_f, "firth")),
    sec("📈", "Interpretasi Parameter (OR — CI Firth)", or_tbl),
    sec("🧭", "Keputusan Variabel OR≈1 / Tdk Signifikan",
        keputusan_html_generic(names(coefs), exp(coefs), pvals)),
    sec("🔍", "Evaluasi Model (in-sample, seluruh data)",
        evaluasi_html(mdl, df$status_sekolah, prob_fit, pid, is_firth = TRUE))
  )
  list(html = html, ringkas = ringkas)
}

# ==================================================================
# SMOTEN — ALUR SAMA DENGAN SKRIPSI KTI (v15)
# 1) Split stratified training-testing dilakukan SEBELUM resampling.
# 2) SMOTEN hanya pada training set (themis::smoten, jarak VDM,
#    karena seluruh prediktor kategorik).
# 3) Testing set tidak pernah di-SMOTE; dipakai hanya untuk evaluasi final.
# 4) Reglog biasa, Firth, dan SMOTEN dibandingkan pada testing set yang SAMA.
# 5) Coding outcome tetap: 1 = Bersekolah, 0 = Tidak Bersekolah (minoritas).
# 6) Cutoff Youden ditentukan dari training asli, dikunci, lalu dipakai pada testing.
# ==================================================================

# ---- Split stratified yang menjaga kedua kelas dan level prediktor ----
smote_stratified_split <- function(df, xvars_model, train_prop = 0.80,
                                   seed = 2026, max_try = 200) {
  y <- df$status_sekolah
  if (length(unique(y)) != 2) stop("Y harus memiliki tepat dua kelas.")

  for (attempt in 0:(max_try - 1)) {
    set.seed(seed + attempt)

    idx_by_class <- split(seq_len(nrow(df)), y)
    idx_train <- unlist(lapply(idx_by_class, function(idx) {
      n_i <- length(idx)
      if (n_i < 2) stop("Salah satu kelas memiliki kurang dari 2 observasi.")
      n_tr <- floor(train_prop * n_i)
      n_tr <- max(1L, min(n_i - 1L, n_tr))
      sample(idx, size = n_tr, replace = FALSE)
    }), use.names = FALSE)

    idx_train <- sort(unique(idx_train))
    idx_test  <- setdiff(seq_len(nrow(df)), idx_train)

    tr <- df[idx_train, , drop = FALSE]
    te <- df[idx_test,  , drop = FALSE]

    kelas_ok <- length(unique(tr$status_sekolah)) == 2 &&
                length(unique(te$status_sekolah)) == 2

    pred_ok <- all(vapply(xvars_model, function(v) {
      lev_tr <- unique(as.character(tr[[v]]))
      lev_te <- unique(as.character(te[[v]]))
      length(lev_tr) >= 2 && all(lev_te %in% lev_tr)
    }, logical(1)))

    if (kelas_ok && pred_ok) {
      for (v in xvars_model) {
        tr[[v]] <- droplevels(as.factor(tr[[v]]))
        te[[v]] <- factor(as.character(te[[v]]), levels = levels(tr[[v]]))
      }
      return(list(train = tr, test = te, seed_used = seed + attempt))
    }
  }

  stop(paste0("Tidak berhasil membentuk split stratified yang valid dalam ", max_try,
              " percobaan (kemungkinan ada kategori prediktor yang sangat jarang). ",
              "Coba ubah seed atau train_prop."))
}

# ---- Konteks split: split + model pembanding (dilatih sekali per split) ----
siapkan_split <- function(data, label, train_prop, xvars = xvars_global, seed = 2026) {
  p  <- prep_data(data, xvars)
  df <- p$df; xm <- p$xvars_model; fm <- p$formula_mod
  if (length(xm) == 0) stop("Tidak ada variabel valid.")

  spl   <- smote_stratified_split(df, xm, train_prop = train_prop, seed = seed)
  train <- spl$train
  test  <- spl$test

  # Reglog biasa pada training asli
  mdl_b  <- glm(fm, data = train, family = binomial)
  mdl_b0 <- glm(status_sekolah ~ 1, data = train, family = binomial)

  # Firth pada training asli (pl = FALSE: CI/p profile tidak diperlukan di sini)
  mdl_f <- tryCatch(logistf(fm, data = train, pl = FALSE), error = function(e) NULL)

  pr_te <- list(); pr_tr <- list()
  pr_te[["Reglog Biasa"]] <- as.numeric(predict(mdl_b, newdata = test,  type = "response"))
  pr_tr[["Reglog Biasa"]] <- as.numeric(predict(mdl_b, newdata = train, type = "response"))
  if (!is.null(mdl_f)) {
    pr_te[["Firth"]] <- prob_firth(mdl_f, test,  xm)
    pr_tr[["Firth"]] <- prob_firth(mdl_f, train, xm)
  }

  ok <- vapply(names(pr_te), function(nm)
    all(is.finite(pr_te[[nm]])) && all(is.finite(pr_tr[[nm]])), logical(1))
  if (!ok[["Reglog Biasa"]])
    stop("Prediksi Reglog Biasa mengandung NA/Inf (kategori sangat jarang atau rank deficiency).")
  pr_te <- pr_te[ok]; pr_tr <- pr_tr[ok]

  cuts <- lapply(pr_tr, function(pp) smote_youden_from_train(train$status_sekolah, pp))

  n_min_train <- sum(train$status_sekolah == 0)
  n_maj_train <- sum(train$status_sekolah == 1)
  pct_min_tr  <- n_min_train / nrow(train) * 100

  # Baris ringkasan baseline (training asli -> testing)
  kel <- paste("Split", label)
  rb  <- list()
  rb[[1]] <- mk_ringkas(kel, "Reglog Biasa", "Training asli → testing",
                        nrow(train), 0, pct_min_tr,
                        smote_metrics_school(test$status_sekolah, pr_te[["Reglog Biasa"]], 0.50),
                        hl_p(train$status_sekolah, pr_tr[["Reglog Biasa"]]),
                        g2_p_glm(mdl_b), "—")
  if ("Firth" %in% names(pr_te)) {
    g2_f <- tryCatch(pchisq(2 * abs(diff(as.numeric(mdl_f$loglik))),
                            length(mdl_f$coefficients) - 1, lower.tail = FALSE),
                     error = function(e) NA_real_)
    rb[[2]] <- mk_ringkas(kel, "Firth", "Training asli → testing",
                          nrow(train), 0, pct_min_tr,
                          smote_metrics_school(test$status_sekolah, pr_te[["Firth"]], 0.50),
                          hl_p(train$status_sekolah, pr_tr[["Firth"]]),
                          g2_f, "—")
  }

  list(
    label = label, train_prop = train_prop, seed_used = spl$seed_used,
    df = df, xm = xm, fm = fm, train = train, test = test,
    pr_te = pr_te, pr_tr = pr_tr, cuts = cuts,
    n_min_train = n_min_train, n_maj_train = n_maj_train,
    ringkas_base = bind_rows(rb),
    firth_ada = "Firth" %in% names(pr_te)
  )
}

# ---- HTML confusion matrix pada testing asli ----
smote_test_block_html <- function(y_true, prob_school, cutoff, label_text) {
  m <- smote_metrics_school(y_true, prob_school, cutoff)
  pct <- function(x) ifelse(is.na(x), "NA", paste0(round(x * 100, 2), " %"))

  cm_tbl <- kbl_h(data.frame(
    ` ` = c("Pred: Tidak Bersekolah (0)", "Pred: Bersekolah (1)"),
    `Aktual: Tidak Bersekolah (0)` = c(m$tn, m$fp),
    `Aktual: Bersekolah (1)`       = c(m$fn, m$tp),
    check.names = FALSE
  ), paste0("Confusion Matrix Testing - ", label_text))

  paste0(
    '<div class="cm-box">',
    '<span class="cm-label cut50">', label_text, '</span>',
    cm_tbl,
    '<div class="metric-row">',
    '<span class="metric-chip">Accuracy: ', pct(m$acc), '</span>',
    '<span class="metric-chip">Sensitivity / Recall Bersekolah: ', pct(m$sens_school), '</span>',
    '<span class="metric-chip">Specificity / Recall Tidak Bersekolah: ', pct(m$recall_noschool), '</span>',
    '<span class="metric-chip">Balanced Accuracy: ', pct(m$bal), '</span>',
    '<span class="metric-chip">Precision Bersekolah: ', pct(m$prec_school), '</span>',
    '<span class="metric-chip">F1 Bersekolah: ', pct(m$f1_school), '</span>',
    '<span class="metric-chip">Precision Tidak Bersekolah: ', pct(m$prec_noschool), '</span>',
    '<span class="metric-chip">F1 Tidak Bersekolah: ', pct(m$f1_noschool), '</span>',
    '<span class="metric-chip">G-mean: ', pct(m$gmean), '</span>',
    '<span class="metric-chip">AUC (Y=1 Bersekolah): ',
      ifelse(is.na(m$auc), "NA", round(m$auc, 4)), '</span>',
    '</div></div>'
  )
}

# ---- Perbandingan tiga metode (Reglog Biasa, Firth, SMOTEN) pada testing yang sama ----
banding_tiga_html <- function(y_test, probs, cuts) {
  mkrow <- function(metode, skema, m) {
    data.frame(
      Metode = metode,
      Cutoff = skema,
      `Nilai cutoff P(Bersekolah)` = round(m$cutoff, 3),
      `Accuracy (%)` = round(m$acc * 100, 2),
      `Sensitivity / Recall Bersekolah (%)` = round(m$sens_school * 100, 2),
      `Specificity / Recall Tidak Bersekolah (%)` = round(m$recall_noschool * 100, 2),
      `Balanced Accuracy (%)` = round(m$bal * 100, 2),
      `F1 Bersekolah (%)` = round(m$f1_school * 100, 2),
      `F1 Tidak Bersekolah (%)` = round(m$f1_noschool * 100, 2),
      `G-mean (%)` = round(m$gmean * 100, 2),
      `AUC Y=1 Bersekolah` = round(m$auc, 4),
      check.names = FALSE
    )
  }

  rows50 <- list(); rowsY <- list(); mets50 <- list(); metsY <- list()
  for (nm in names(probs)) {
    mets50[[nm]] <- smote_metrics_school(y_test, probs[[nm]], 0.50)
    metsY[[nm]]  <- smote_metrics_school(y_test, probs[[nm]], cuts[[nm]])
    rows50[[nm]] <- mkrow(nm, "Fixed 0.50", mets50[[nm]])
    rowsY[[nm]]  <- mkrow(nm, "Youden dari training", metsY[[nm]])
  }
  tbl <- bind_rows(c(unname(rows50), unname(rowsY)))

  # Keputusan utama: testing asli, cutoff tetap 0.50
  pemenang <- pilih_pemenang(mets50)

  sm <- mets50[["SMOTEN + Reglog"]]
  kal <- paste0(
    '<p><strong>Metode terbaik (testing asli, cutoff 0.50):</strong> ', pemenang, '. ',
    'Urutan penentu: Balanced Accuracy, lalu recall kelas Tidak Bersekolah, lalu AUC.</p>'
  )
  for (pembanding in setdiff(names(mets50), "SMOTEN + Reglog")) {
    mb <- mets50[[pembanding]]
    kal <- paste0(kal,
      '<p>SMOTEN terhadap ', pembanding, ': selisih Balanced Accuracy = ',
      round((sm$bal - mb$bal) * 100, 2), ' poin persentase; selisih recall Tidak Bersekolah = ',
      round((sm$recall_noschool - mb$recall_noschool) * 100, 2),
      ' poin persentase; selisih sensitivity Bersekolah = ',
      round((sm$sens_school - mb$sens_school) * 100, 2), ' poin persentase.</p>')
  }
  kal <- paste0(kal,
    '<p><em>Outcome regresi tetap dikodekan 1 = Bersekolah dan 0 = Tidak Bersekolah. ',
    'Karena kelas Tidak Bersekolah merupakan kelas minoritas, keberhasilan penanganan imbalance ',
    'dinilai terutama dari Balanced Accuracy dan recall kelas Tidak Bersekolah, tanpa membalik ',
    'event model logistik.</em></p>')

  list(
    html = paste0(kbl_h(tbl, "Perbandingan Reglog Biasa, Firth, dan SMOTEN pada Testing Set Asli"), kal),
    mets50 = mets50, pemenang = pemenang
  )
}

# ---- Satu skenario SMOTEN (satu split x satu persen Chawla) ----
run_smote <- function(ctx, nama_dataset, pct, label_skema, K = 5) {
  xm <- ctx$xm; fm <- ctx$fm; train <- ctx$train; test <- ctx$test
  n_min <- ctx$n_min_train; n_maj <- ctx$n_maj_train
  gagal <- function(msg) list(html = paste0("<p>⚠️ ", msg, "</p>"), ringkas = NULL)

  if (!requireNamespace("themis", quietly = TRUE))
    return(gagal("Paket <code>themis</code> diperlukan untuk SMOTEN. Instal dengan <code>install.packages('themis')</code>."))

  pred_semua_kategorik <- all(vapply(train[xm], function(z) is.factor(z) || is.character(z), logical(1)))
  if (!pred_semua_kategorik)
    return(gagal("SMOTEN hanya tepat bila seluruh prediktor kategorik. Gunakan SMOTENC untuk prediktor numerik."))

  # Konversi persen Chawla -> over_ratio (themis::smoten)
  target_akhir <- round(n_min * (1 + pct / 100))
  over_ratio   <- min(1, target_akhir / n_maj)

  # Pengaman K: kelas minoritas pada training harus >= K + 1 observasi
  K_pakai <- K; cat_K <- ""
  if (n_min < K + 1) {
    K_pakai <- n_min - 1
    cat_K <- paste0(' | <em>K diturunkan otomatis dari ', K, ' menjadi ', K_pakai,
                    ' karena kelas minoritas pada training hanya ', n_min, ' observasi.</em>')
  }
  if (K_pakai < 1)
    return(gagal(paste0("SMOTEN tidak dapat dijalankan: kelas minoritas pada training hanya ", n_min, " observasi.")))

  # ================================================================
  # 1. SMOTEN HANYA PADA TRAINING SET
  # ================================================================
  train_sm_in <- train[, c("status_sekolah", xm), drop = FALSE]
  train_sm_in$status_sekolah <- factor(
    train_sm_in$status_sekolah,
    levels = c(0, 1),
    labels = c("Tidak Bersekolah", "Bersekolah")
  )
  for (v in xm) {
    train_sm_in[[v]] <- factor(as.character(train_sm_in[[v]]), levels = levels(train[[v]]))
  }

  set.seed(ctx$seed_used + 1000L)
  df_sm <- tryCatch(
    themis::smoten(df = train_sm_in, var = "status_sekolah", k = K_pakai, over_ratio = over_ratio),
    error = function(e) e
  )
  if (inherits(df_sm, "error"))
    return(gagal(paste0("SMOTEN gagal: ", df_sm$message)))

  df_sm <- as.data.frame(df_sm)
  df_sm$status_sekolah <- ifelse(as.character(df_sm$status_sekolah) == "Bersekolah", 1L, 0L)
  for (v in xm) {
    df_sm[[v]] <- factor(as.character(df_sm[[v]]), levels = levels(train[[v]]))
  }

  # ================================================================
  # 2. FIT REGRESI LOGISTIK SETELAH SMOTEN
  # ================================================================
  mdl  <- tryCatch(glm(fm, data = df_sm, family = binomial), error = function(e) NULL)
  mdl0 <- tryCatch(glm(status_sekolah ~ 1, data = df_sm, family = binomial), error = function(e) NULL)
  if (is.null(mdl) || is.null(mdl0))
    return(gagal("Model GLM setelah SMOTEN gagal di-fit."))

  # ================================================================
  # 3. EVALUASI FINAL HANYA PADA TESTING SET ASLI
  # ================================================================
  prob_te_s <- tryCatch(as.numeric(predict(mdl, newdata = test, type = "response")),
                        error = function(e) rep(NA_real_, nrow(test)))
  if (any(!is.finite(prob_te_s)))
    return(gagal("Prediksi testing SMOTEN mengandung NA/Inf. Periksa kategori sangat jarang atau rank deficiency."))

  prob_tr_s <- as.numeric(predict(mdl, newdata = train, type = "response"))
  cut_s     <- smote_youden_from_train(train$status_sekolah, prob_tr_s)

  probs <- c(ctx$pr_te, list("SMOTEN + Reglog" = prob_te_s))
  cuts  <- c(ctx$cuts,  list("SMOTEN + Reglog" = cut_s))

  bandingan <- banding_tiga_html(test$status_sekolah, probs, cuts)

  # ================================================================
  # 4. EDA DISTRIBUSI TOTAL, TRAINING, SMOTEN, TESTING
  # ================================================================
  tab2 <- function(x) as.integer(table(factor(x, levels = c(0, 1))))
  eda_tbl <- kbl_h(data.frame(
    Kategori = c("Tidak Bersekolah (0) - minoritas", "Bersekolah (1) - event model"),
    `N Total Asli`      = tab2(ctx$df$status_sekolah),
    `N Training Asli`   = tab2(train$status_sekolah),
    `N Training SMOTEN` = tab2(df_sm$status_sekolah),
    `N Testing Asli`    = tab2(test$status_sekolah),
    check.names = FALSE
  ), "Distribusi Y: Total, Training, Setelah SMOTEN, dan Testing")

  info <- paste0(
    '<p><strong>Skema:</strong> ', label_skema,
    ' | <strong>target minoritas pada training:</strong> ', target_akhir,
    ' (dari ', n_min, ' asli) | <strong>over_ratio:</strong> ', round(over_ratio, 4),
    ' | <strong>K:</strong> ', K_pakai, cat_K, '</p>',
    '<p><strong>n total valid:</strong> ', nrow(ctx$df),
    ' | <strong>training asli:</strong> ', nrow(train),
    ' | <strong>training setelah SMOTEN:</strong> ', nrow(df_sm),
    ' | <strong>testing asli:</strong> ', nrow(test),
    ' | <strong>split:</strong> ', round(ctx$train_prop * 100), ':', round((1 - ctx$train_prop) * 100),
    ' stratified | <strong>seed:</strong> ', ctx$seed_used, '.</p>',
    '<p><em>SMOTEN diterapkan hanya pada training set. Testing set tetap pada distribusi asli ',
    'dan tidak ikut dalam pembentukan observasi sintetis. Prediktor kategorik dipertahankan ',
    'sebagai faktor sehingga kategori sintetis dibentuk dengan metode nominal (VDM), bukan interpolasi dummy. ',
    'Coding outcome tidak dibalik: 1 = Bersekolah dan 0 = Tidak Bersekolah.</em></p>'
  )

  # Blok evaluasi testing untuk tiap metode, dua skema cutoff
  blok50 <- paste0(vapply(names(probs), function(nm)
    smote_test_block_html(test$status_sekolah, probs[[nm]], 0.50,
                          paste0(nm, " | cutoff P(Bersekolah) 0.50")), character(1)), collapse = "")
  blokY  <- paste0(vapply(names(probs), function(nm)
    smote_test_block_html(test$status_sekolah, probs[[nm]], cuts[[nm]],
                          paste0(nm, " | Youden training = ", round(cuts[[nm]], 3))), character(1)), collapse = "")

  ses_s <- tryCatch(sqrt(diag(vcov(mdl))), error = function(e) rep(NA, length(coef(mdl))))
  if (length(ses_s) != length(coef(mdl))) ses_s <- rep(NA_real_, length(coef(mdl)))

  eval_html <- paste0(
    '<div class="eval-wrap">',
    separation_html(coef(mdl), ses_s, FALSE),
    '<div class="eval-sub"><h6>Testing Set Asli - Cutoff Tetap 0.50</h6>',
    '<p><em>Event model tetap Y=1 = Bersekolah. Recall kelas Tidak Bersekolah (Y=0) dilaporkan eksplisit untuk menilai kelas minoritas.</em></p>',
    '<div class="cm-wrap">', blok50, '</div></div>',
    '<div class="eval-sub"><h6>Testing Set Asli - Cutoff Youden dari Training</h6>',
    '<p><em>Cutoff dipilih dari training asli dan kemudian dikunci sebelum diterapkan pada testing.</em></p>',
    '<div class="cm-wrap">', blokY, '</div></div>',
    '</div>'
  )

  # ================================================================
  # 5. INFERENSI MODEL SMOTEN PADA TRAINING RESAMPLED
  # ================================================================
  k_smote <- length(coef(mdl))
  n_smote <- nrow(df_sm)

  koef_s <- summary(mdl)$coefficients
  ci_s   <- tryCatch(exp(confint.default(mdl)), error = function(e) NULL)
  nm_s   <- rownames(koef_s)
  catat_or(nama_dataset, "SMOTE", nm_s, exp(coef(mdl)[nm_s]), koef_s[, 4], n_smote,
           if (!is.null(ci_s)) ci_s[nm_s, 1] else NA_real_,
           if (!is.null(ci_s)) ci_s[nm_s, 2] else NA_real_)

  ringkas <- mk_ringkas(
    paste("Split", ctx$label), "SMOTEN + Reglog", label_skema,
    nrow(train), nrow(df_sm) - nrow(train),
    mean(df_sm$status_sekolah == 0) * 100,
    bandingan$mets50[["SMOTEN + Reglog"]],
    hl_p(df_sm$status_sekolah, fitted(mdl)), g2_p_glm(mdl),
    bandingan$pemenang
  )

  html <- paste0(
    info,
    sec("📊", "EDA - Distribusi Y dan Pembagian Data", eda_tbl),
    sec("⚖️", "Perbandingan Reglog Biasa, Firth, dan SMOTEN pada Testing Asli", bandingan$html),
    sec("✅", "Goodness of Fit (H-L Test) - Model SMOTEN", hl_html(mdl, df_sm$status_sekolah)),
    sec("📐", "Uji Simultan (G²) - Model SMOTEN", simultan_html(mdl, mdl0, xm)),
    sec("📋", "Uji Parsial (Wald Test) - Model SMOTEN", wald_html(mdl)),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC) - Model SMOTEN",
        paste0(
          ic_html(mdl, n_smote, k_smote, "glm"),
          '<p><em><small>&#9432; AIC/AICc/BIC dihitung pada training setelah SMOTEN. ',
          'Perbandingan performa prediktif utama menggunakan testing set asli yang tidak di-resampling.</small></em></p>'
        )),
    sec("📈", "Interpretasi Parameter (OR, CI 95%) - Model SMOTEN", or_html(mdl)),
    sec("🧭", "Keputusan Variabel OR≈1 / Tdk Signifikan - Model SMOTEN",
        keputusan_variabel_html(mdl)),
    sec("🔍", "Evaluasi Model pada Testing Set Asli", eval_html)
  )
  list(html = html, ringkas = ringkas)
}

for (sp_id in names(split_list)) {
  sp <- split_list[[sp_id]]
  cat("\n\n### Split ", sp$label, "\n\n", sep = "")

  if (!is.null(ctx_err[[sp_id]])) {
    cat('<p style="color:#c62828;font-weight:600;">&#9888; Pembagian data gagal: ',
        ctx_err[[sp_id]], '</p>\n', sep = "")
    next
  }

  cx <- ctx_list[[sp_id]]
  g  <- grid_info[[sp_id]]
  n_min <- cx$n_min_train; n_maj <- cx$n_maj_train

  cat('<p><strong>Training asli:</strong> ', nrow(cx$train),
      ' observasi (Tidak Bersekolah = ', n_min, '; Bersekolah = ', n_maj, '; minoritas = ',
      round(n_min / nrow(cx$train) * 100, 2), '%). <strong>Testing asli:</strong> ', nrow(cx$test),
      ' observasi. <strong>Titik seimbang penuh:</strong> ', g$pct_full, '%.</p>\n', sep = "")

  if (length(g$grid) == 0) {
    cat('<p style="color:#c62828;font-weight:600;">&#9888; Kelas Tidak Bersekolah bukan kelas minoritas ',
        'pada training (atau titik seimbang tidak terdefinisi). SMOTEN tidak dijalankan pada split ini.</p>\n', sep = "")
    next
  }
  if (g$pct_full < PCT_MULAI) {
    cat('<p style="color:#e65100;font-weight:600;">&#9888; Titik seimbang penuh (', g$pct_full,
        '%) lebih kecil dari PCT_MULAI (', PCT_MULAI,
        '%). Hanya titik seimbang penuh yang dijalankan pada split ini.</p>\n', sep = "")
  }
  if (length(g$grid) > 15) {
    cat('<p style="color:#e65100;font-weight:600;">&#9888; Jumlah titik grid pada split ini = ', length(g$grid),
        ' (lebih dari 15). Waktu knit bisa lama; PCT_LANGKAH dapat diperbesar.</p>\n', sep = "")
  }

  ta <- round(n_min * (1 + g$grid / 100))
  tab_grid <- data.frame(
    `Persen Chawla` = paste0(g$grid, "%"),
    Jenis           = ifelse(g$grid == g$pct_full, "Seimbang penuh", "Chawla"),
    `Target minoritas (training)` = ta,
    over_ratio      = round(pmin(1, ta / n_maj), 4),
    `Proporsi minoritas akhir (%)` = round(ta / (ta + n_maj) * 100, 2),
    check.names = FALSE
  )
  cat(kbl_h(tab_grid, paste0("Grid SMOTEN — Split ", sp$label, " (", length(g$grid), " titik)")))
}