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

# ── Visualisasi ───────────────────────────────────────────────────
library(highcharter)       # interactive charts
library(htmltools)         # tagList, HTML

# ── Pemodelan ─────────────────────────────────────────────────────
library(ResourceSelection) # Hosmer-Lemeshow
library(lmtest)            # lrtest
library(logistf)           # Firth logistic regression
library(smotefamily)       # SMOTE oversampling
library(pROC)              # ROC-AUC
library(car)               # VIF

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

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

allowed_vals <- c(1, 2, 3, 5, 6, 7)

dataku_filtered <- data_susenas %>%
  mutate(
    across(R1002:R1009, ~ as.numeric(.x)),

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

    prov = as.numeric(R101),
    kabu = as.numeric(R101) * 100 + as.numeric(R102),

    # Y: Partisipasi Sekolah (1 = Bersekolah, 0 = Tidak Bersekolah)
    status = ifelse(
      R407 >= 13 & R407 <= 15 & JUMLAH_JENIS_DISABILITAS > 0,
      ifelse(R610 == 2, 1, 0),
      NA_real_
    )
  ) %>%
  filter(
    R407 >= 13 & R407 <= 15,
    JUMLAH_JENIS_DISABILITAS > 0,
    as.numeric(R101) %in% c(52:53, 71:76, 81:82, 91:97)
  )

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

id_rt <- c("URUT")

data_krt <- data_susenas %>%
  filter(as.numeric(R403) == 1) %>%
  mutate(
    R405 = as.numeric(R405),
    R614 = as.numeric(R614),
    R705 = as.numeric(R705),
    R706 = as.numeric(R706),

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

    # X5: Tingkat Pendidikan KRT
    X5 = case_when(
      R614 >= 11 & R614 <= 24   ~ 1L,
      R614 %in% c(0:10, 25)     ~ 0L,
      TRUE                      ~ NA_integer_
    ),

    # 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(id_rt), X4, X5, X6)

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 = as.numeric(R101),
    R102 = as.numeric(R102),
    URUT = as.numeric(URUT)
  )

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 = as.numeric(R101))

data_kp_x9 <- data_kp_43 %>%
  left_join(GK, by = c("R101" = "R101", "R105" = "R105")) %>%
  mutate(
    # X9: Status Kemiskinan RT (1 = Miskin, 0 = Tidak Miskin)
    X9 = if_else(KAPITA < GK, 1L, 0L, missing = NA_integer_)
  ) %>%
  select(all_of(id_rt), X9)

dataku_final <- dataku_filtered %>%
  mutate(
    across(R1002:R1009, ~ as.numeric(.x)),
    R105 = as.numeric(R105),
    R405 = as.numeric(R405),
    R615 = as.numeric(R615),
    R616 = as.numeric(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 (Tunggal / Ganda)
    X2 = case_when(
      JUMLAH_JENIS_DISABILITAS > 1  ~ 1L,
      JUMLAH_JENIS_DISABILITAS == 1 ~ 0L,
      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
    )
  ) %>%
  left_join(data_krt,   by = id_rt) %>%
  left_join(data_kp_x9, by = id_rt) %>%
  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 tunggal", "Disabilitas ganda")),
    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("Tidak miskin", "Miskin"))
  )

data.frame(
  Kode     = c("Y", paste0("X", 1:9)),
  Nama     = c(
    "Partisipasi Sekolah",
    "Jenis Kelamin Anak",
    "Jenis Disabilitas",
    "Tingkat Kesulitan Disabilitas",
    "Jenis Kelamin KRT",
    "Pendidikan KRT",
    "Pekerjaan KRT",
    "Daerah Tempat Tinggal",
    "Kepemilikan KIP/PIP",
    "Status Kemiskinan RT"
  ),
  Kategori = c(
    "Dependen",
    rep("Karakteristik Anak", 3),
    rep("Karakteristik KRT",  3),
    rep("Sosial-Ekonomi RT",  3)
  ),
  Koding = c(
    "1 = Bersekolah; 0 = Tidak bersekolah",
    "1 = Laki-laki; 0 = Perempuan",
    "1 = Ganda (≥2 jenis); 0 = Tunggal",
    "0 = Sangat kesulitan; 1 = Kesulitan; 2 = Sedikit kesulitan",
    "1 = Laki-laki; 0 = Perempuan",
    "1 = ≥ SMA; 0 = ≤ SMP",
    "0 = Tidak bekerja; 1 = Pertanian; 2 = Nonpertanian",
    "1 = Perkotaan; 0 = Perdesaan",
    "1 = Memiliki; 0 = Tidak memiliki",
    "1 = Miskin; 0 = Tidak miskin"
  )
) %>%
  kable(caption = "Referensi Variabel Penelitian (Y dan X1–X9)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE) %>%
  column_spec(1, bold = TRUE, width = "4em") %>%
  column_spec(3, italic = TRUE) %>%
  row_spec(1, background = "#fff3e0")

cat("Jumlah observasi :", nrow(dataku_final), "\n")

cat("Jumlah variabel  :", ncol(dataku_final), "\n\n")

glimpse(dataku_final)

summary(dataku_final)

missing_tbl <- dataku_final %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(everything(),
               names_to  = "Variabel",
               values_to = "Jumlah_NA") %>%
  mutate(Persen_NA = round(Jumlah_NA / nrow(dataku_final) * 100, 2)) %>%
  arrange(desc(Jumlah_NA))

kable(missing_tbl,
      caption = "Jumlah dan Persentase Missing Value per Variabel",
      col.names = c("Variabel", "Jumlah NA", "% NA")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

tbl_y <- dataku_final %>%
  count(Y) %>%
  mutate(Persen = round(n / sum(n) * 100, 2))

kable(tbl_y,
      caption = "Distribusi Partisipasi Sekolah Anak Disabilitas Usia 13–15 Tahun",
      col.names = c("Status Sekolah", "Frekuensi", "%")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)

# ── Pie Chart Distribusi Y ──────────────────────────────────────
pd_y <- tbl_y %>%
  mutate(name = as.character(Y), y = Persen) %>%
  select(name, y)

highchart() %>%
  hc_chart(type = "pie", backgroundColor = "transparent") %>%
  hc_title(text  = "Distribusi Partisipasi Sekolah",
           style = list(fontSize = "15px", fontWeight = "700", color = "#333")) %>%
  hc_subtitle(text = "Anak Penyandang Disabilitas Usia 13\u201315 Tahun (SUSENAS 2023)",
              style = list(fontSize = "12px", color = "#666")) %>%
  hc_series(list(name = "Persentase", colorByPoint = TRUE,
                 data = list_parse(pd_y))) %>%
  hc_colors(c("#D73027","#4575B4")) %>%
  hc_plotOptions(pie = list(
    allowPointSelect = TRUE, cursor = "pointer",
    dataLabels = list(enabled = TRUE,
                      format  = "<b>{point.name}</b><br>n = {point.y:.0f}%",
                      style   = list(fontSize = "13px", textOutline = "none")),
    showInLegend = TRUE
  )) %>%
  hc_tooltip(pointFormat = "<b>{point.name}</b>: {point.y:.2f}%") %>%
  hc_exporting(enabled = TRUE, filename = "distribusi_y",
               buttons = list(contextButton = list(
                 menuItems = c("downloadPNG","downloadSVG","downloadPDF"),
                 text = "Unduh HD")))

tbl_dis <- dataku_filtered %>%
  count(JUMLAH_JENIS_DISABILITAS) %>%
  mutate(Persen = round(n / sum(n) * 100, 2))

kable(tbl_dis,
      caption  = "Distribusi Jumlah Jenis Disabilitas",
      col.names = c("Jumlah Jenis Disabilitas", "Frekuensi", "%")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)

make_col(
  cats  = tbl_dis$JUMLAH_JENIS_DISABILITAS,
  vals  = tbl_dis$n,
  title = "Distribusi Jumlah Jenis Disabilitas",
  color = "#4575B4",
  xlab  = "Jumlah Jenis Disabilitas",
  ylab  = "Frekuensi"
)

tbl_prov <- dataku_final %>%
  count(prov) %>%
  mutate(Persen = round(n / sum(n) * 100, 2)) %>%
  arrange(desc(n))

kable(tbl_prov,
      caption   = "Distribusi Observasi per Provinsi",
      col.names = c("Kode Provinsi", "Frekuensi", "%")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)

tbl_prov_sorted <- tbl_prov %>% arrange(n)

make_col(
  cats       = tbl_prov_sorted$prov,
  vals       = tbl_prov_sorted$n,
  title      = "Jumlah Observasi per Provinsi",
  color      = "#74ADD1",
  xlab       = "Kode Provinsi",
  ylab       = "Frekuensi",
  horizontal = TRUE
)

vars_for_tab <- c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9")
label_vars   <- c(
  "Jenis Kelamin Anak (X1)",
  "Jenis Disabilitas (X2)",
  "Tingkat Kesulitan (X3)",
  "Jenis Kelamin KRT (X4)",
  "Pendidikan KRT (X5)",
  "Pekerjaan KRT (X6)",
  "Daerah Tinggal (X7)",
  "Kepemilikan KIP/PIP (X8)",
  "Status Kemiskinan RT (X9)"
)

for (i in seq_along(vars_for_tab)) {
  var <- vars_for_tab[i]
  cat("\n###", label_vars[i], "\n\n")

  df_tmp <- dataku_final %>%
    filter(!is.na(.data[[var]]), !is.na(Y))

  tab     <- table(df_tmp[[var]], df_tmp$Y)
  tab_pct <- prop.table(tab, margin = 1) * 100

  tbl_out <- as.data.frame.matrix(tab) %>%
    mutate(
      Total          = rowSums(.),
      `% Bersekolah` = round(tab_pct[, "Bersekolah"], 1)
    )

  chi_res <- tryCatch(chisq.test(tab), error = function(e) NULL)
  p_val   <- if (!is.null(chi_res)) round(chi_res$p.value, 4) else NA

  print(
    kable(tbl_out,
          caption = paste0(label_vars[i], " — Chi-Square p-value: ", p_val)) %>%
      kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                    full_width = FALSE)
  )
}

ringkasan <- lapply(seq_along(vars_for_tab), function(i) {
  var <- vars_for_tab[i]
  dataku_final %>%
    filter(!is.na(.data[[var]]), !is.na(Y)) %>%
    group_by(Variabel = label_vars[i], Kategori = .data[[var]]) %>%
    summarise(
      n              = n(),
      n_bersekolah   = sum(Y == "Bersekolah"),
      pct_bersekolah = round(n_bersekolah / n * 100, 1),
      .groups        = "drop"
    )
}) %>% bind_rows()

kable(ringkasan,
      caption   = "Ringkasan % Bersekolah per Kategori Prediktor",
      col.names = c("Variabel", "Kategori", "n", "n Bersekolah", "% Bersekolah")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  collapse_rows(columns = 1, valign = "top")

library(survey)

allowed_vals_ind <- c(1, 2, 3, 5, 6, 7)

# ── Bentuk data indikator dari data_susenas mentah ────────────────────────
indicator_disab <- data_susenas %>%
  mutate(
    across(R1002:R1009, ~ as.numeric(.x)),
    R101 = as.numeric(R101),
    R102 = as.numeric(R102),
    R105 = as.numeric(R105),
    R407 = as.numeric(R407),
    R610 = as.numeric(R610),
    FWT  = as.numeric(FWT),

    # Jumlah jenis disabilitas (nilai 1,2,3,5,6,7 = ada kesulitan)
    JUMLAH_DISAB = rowSums(
      across(R1002:R1009, ~ .x %in% allowed_vals_ind),
      na.rm = TRUE
    ),

    # Strata: kombinasi kode provinsi × daerah (perkotaan/perdesaan)
    STRATA_IND = paste0(sprintf("%02d", R101), R105),

    jml_pddk = 100,

    # Variabel respons: 100 = bersekolah (R610==2), 0 = tidak
    pa = ifelse(R610 == 2, 100, 0)
  ) %>%
  filter(
    R407 >= 13 & R407 <= 15,          # Usia 13–15 tahun
    JUMLAH_DISAB > 0,                 # Penyandang disabilitas
    R101 %in% c(52:53, 71:76, 81:82, 91:97)
    
  )

cat("Jumlah observasi domain indikator:", nrow(indicator_disab), "\n")

# ── Desain survei ─────────────────────────────────────────────────────────
options(survey.adjust.domain.lonely = TRUE)
options(survey.lonely.psu           = "adjust")

susenas_ind_design <- svydesign(
  id      = ~PSU + SSU,
  strata  = ~STRATA_IND,
  data    = indicator_disab,
  weights = ~FWT,
  nest    = TRUE
)

summary(susenas_ind_design$prob)

# ── Estimasi langsung per provinsi ────────────────────────────────────────
resultP_ind <- svyby(
  formula = ~pa,
  denom   = ~jml_pddk,
  by      = ~R101,
  design  = susenas_ind_design,
  deff    = TRUE,
  svyratio,
  vartype = c("se", "ci", "cv", "cvpct", "var")
)
resultP_ind[is.na(resultP_ind)] <- 0

# ── Hitung turunan statistik ──────────────────────────────────────────────
resultP_ind <- resultP_ind %>%
  mutate(
    theta    = round(`pa/jml_pddk` * 100, 2),
    SE       = round(`se.pa/jml_pddk` * 100, 2),
    VAR      = round(SE^2, 2),
    CI_LOWER = round(ci_l * 100, 2),
    CI_UPPER = round(ci_u * 100, 2),
    RSE      = round(`cv%`, 2),
    DEFF     = round(DEff, 2)
  )
resultP_ind$CI_LOWER[resultP_ind$CI_LOWER < 0] <- 0

# ── Label Nama Provinsi Lengkap (Kode BPS) ───────────────────────────────────
prov_label_ind <- data.frame(
  R101 = c(
    11, 12, 13, 14, 15, 16, 17, 18, 19, 21, 
    31, 32, 33, 34, 35, 36, 
    51, 52, 53, 
    61, 62, 63, 64, 65, 
    71, 72, 73, 74, 75, 76, 
    81, 82, 
    91, 94
  ),
  Nama_Prov = c(
    "Aceh", "Sumatera Utara", "Sumatera Barat", "Riau", "Jambi", 
    "Sumatera Selatan", "Bengkulu", "Lampung", "Kepulauan Bangka Belitung", "Kepulauan Riau",
    "DKI Jakarta", "Jawa Barat", "Jawa Tengah", "DI Yogyakarta", "Jawa Timur", 
    "Banten",
    "Bali", "Nusa Tenggara Barat", "Nusa Tenggara Timur",
    "Kalimantan Barat", "Kalimantan Tengah", "Kalimantan Selatan", "Kalimantan Timur", "Kalimantan Utara",
    "Sulawesi Utara", "Sulawesi Tengah", "Sulawesi Selatan", "Sulawesi Tenggara", "Gorontalo", 
    "Sulawesi Barat",
    "Maluku", "Maluku Utara",
    "Papua Barat", "Papua"
  ),
  Wilayah = c(
    rep("Sumatera", 10),
    rep("Jawa", 6),
    "Bali", rep("Nusa Tenggara", 2),
    rep("Kalimantan", 5),
    rep("Sulawesi", 6),
    rep("Maluku", 2),
    rep("Papua", 2)
  )
)

# ── Tabel output bersih ───────────────────────────────────────────────────
outputP_ind <- resultP_ind %>%
  select(R101, theta, SE, VAR, CI_LOWER, CI_UPPER, RSE, DEFF) %>%
  rename(
    Prov      = R101,
    Estimasi  = theta,
    CI_LOWER  = CI_LOWER,
    CI_UPPER  = CI_UPPER
  ) %>%
  left_join(prov_label_ind, by = c("Prov" = "R101")) %>%
  arrange(Prov)

kable(
  outputP_ind[, c("Nama_Prov", "Wilayah", "Estimasi", "SE",
                  "CI_LOWER", "CI_UPPER", "RSE", "DEFF")],
  caption   = "Estimasi Langsung: % Penyandang Disabilitas 13–15 Tahun yang Bersekolah per Provinsi",
  col.names = c("Provinsi", "Wilayah", "Estimasi (%)", "SE",
                "CI Lower", "CI Upper", "RSE (%)", "DEFF"),
  digits    = 2
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE) %>%
  column_spec(3, bold = TRUE, color = "#1565c0") %>%
  column_spec(7, color = ifelse(outputP_ind$RSE > 25, "red", "darkgreen"),
              bold = ifelse(outputP_ind$RSE > 25, TRUE, FALSE)) %>%
  row_spec(which(outputP_ind$RSE > 25), background = "#fff3e0")

xvars_global <- paste0("X", 1:9)

# ── Konversi kable ke string HTML ─────────────────────────────────
kbl_h <- function(df, caption = NULL) {
  knitr::kable(df, format = "html", caption = caption, escape = FALSE) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed"),
      full_width = FALSE, font_size = 13, position = "left"
    ) %>%
    as.character()
}

# ── Pra-pemrosesan data standar ───────────────────────────────────
# Y (factor "Tidak bersekolah"/"Bersekolah") dikonversi ke status_miskin 0/1
# agar semua fungsi helper berjalan langsung tanpa modifikasi tambahan
prep_data <- function(data) {
  df <- data %>%
    select(all_of(c("Y", xvars_global))) %>%
    mutate(
      status_miskin = as.integer(Y) - 1L,   # 0 = Tidak bersekolah, 1 = Bersekolah
      across(all_of(xvars_global), as.numeric)
    ) %>%
    select(-Y) %>%
    na.omit()

  konstan     <- xvars_global[sapply(xvars_global, function(v) length(unique(df[[v]])) < 2)]
  xvars_model <- setdiff(xvars_global, konstan)
  formula_mod <- as.formula(paste("status_miskin ~", paste(xvars_model, collapse = " + ")))
  list(df = df, xvars_model = xvars_model, formula_mod = formula_mod, konstan = konstan)
}

# ── 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(df$status_miskin)
  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 <- length(xvars_model)
  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) {
  koef <- summary(model)$coefficients
  kbl_h(data.frame(
    Variabel   = rownames(koef),
    B          = round(koef[, 1], 4),
    SE         = round(koef[, 2], 4),
    Wald       = round((koef[, 1] / koef[, 2])^2, 4),
    `p-value`  = round(koef[, 4], 4),
    Kesimpulan = ifelse(koef[, 4] < 0.05, "✅ Signifikan", "❌ Tdk Signifikan"),
    check.names = FALSE
  ), "Uji Parsial (Wald Test)")
}

# ── 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)")
}

# ── Odds Ratio (GLM-based) ───────────────────────────────────────
or_html <- function(model) {
  or <- exp(coef(model))
  ci <- exp(confint.default(model))
  kbl_h(data.frame(
    Variabel     = names(or),
    `exp(B)`     = round(or, 4),
    `CI 2.5%`    = round(ci[, 1], 4),
    `CI 97.5%`   = round(ci[, 2], 4),
    Interpretasi = ifelse(names(or) == "(Intercept)", "Baseline odds",
                   ifelse(or > 1,
                          "↑ Lebih berisiko tidak bersekolah",
                          "↓ Lebih kecil risiko tidak bersekolah")),
    check.names = FALSE
  ), "Interpretasi Parameter (Odds Ratio)")
}

# ── VIF (Multikolinearitas) ──────────────────────────────────────
vif_html <- function(model, df, xm, method_label = "GLM Biasa") {
  vif_val    <- tryCatch(car::vif(model), error = function(e) NULL)
  note_extra <- ""
  if (is.null(vif_val)) {
    fm_vif  <- as.formula(paste("status_miskin ~", paste(xm, collapse = " + ")))
    mdl_vif <- tryCatch(glm(fm_vif, data = df, family = binomial), error = function(e) NULL)
    vif_val <- if (!is.null(mdl_vif)) tryCatch(car::vif(mdl_vif), error = function(e) NULL) else NULL
    note_extra <- paste0(
      '<p><em><small>&#9432; VIF dihitung ulang via GLM biasa (binomial) karena model utama (',
      method_label, ') tidak kompatibel langsung dengan car::vif(). ',
      'VIF adalah properti matriks desain X — hasil tetap valid.</small></em></p>')
  }
  if (is.null(vif_val))
    return('<p><em>⚠️ VIF tidak dapat dihitung.</em></p>')

  if (is.matrix(vif_val)) {
    vif_vec  <- setNames(vif_val[, "GVIF^(1/(2*Df))"]^2, rownames(vif_val))
    col_note <- " (GVIF disesuaikan untuk faktor multi-level)"
  } else {
    vif_vec  <- vif_val
    col_note <- ""
  }

  vif_df <- data.frame(
    Variabel = names(vif_vec),
    VIF      = round(vif_vec, 4),
    Kategori = ifelse(vif_vec > 10, "❌ Parah (>10)",
               ifelse(vif_vec > 5,  "⚠️ Sedang (5–10)",
                                    "✅ Tidak Ada (<5)")),
    check.names = FALSE
  )

  status <- if (any(vif_vec > 10, na.rm = TRUE)) {
    '<p style="color:#c62828;font-weight:600;">⚠️ Terdeteksi multikolinearitas parah.</p>'
  } else if (any(vif_vec > 5, na.rm = TRUE)) {
    '<p style="color:#e65100;font-weight:600;">⚠️ Terdeteksi multikolinearitas sedang. Perlu diwaspadai.</p>'
  } else {
    '<p style="color:#2e7d32;font-weight:600;">✅ Tidak terdeteksi masalah multikolinearitas (semua VIF < 5).</p>'
  }

  paste0(status, note_extra, kbl_h(vif_df, paste0("Variance Inflation Factor (VIF)", col_note)))
}

# ── 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

  if (method == "quasi") {
    phi_est  <- tryCatch(summary(model)$dispersion, error = function(e) 1)
    dev_val  <- model$deviance
    df_r     <- model$df.residual
    qaic_val <- round(dev_val / phi_est + 2 * k, 3)
    rows <- paste0(
      '<tr><td><strong>QAIC</strong></td><td>', qaic_val,
      '</td><td>Quasi-AIC (φ̂ = ', round(phi_est, 3), ')</td></tr>',
      '<tr><td>Deviance Residual</td><td>', round(dev_val, 3),
      '</td><td>df = ', df_r, '</td></tr>',
      '<tr><td>Dispersi (φ̂)</td><td>', round(phi_est, 3),
      '</td><td>Estimasi dispersi quasibinomial</td></tr>',
      '<tr><td style="color:#888;font-style:italic">AIC / BIC</td>',
      '<td style="color:#888">N/A</td>',
      '<td style="color:#888;font-style:italic">Tidak terdefinisi untuk quasibinomial</td></tr>'
    )
  } else {
    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></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  <- 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 ────────────────────────────────
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) 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, quiet = TRUE), error = function(e) NULL)
  cut_you <- if (!is.null(roc_obj)) {
    co <- pROC::coords(roc_obj, "best", best.method = "youden", ret = "threshold")
    as.numeric(co[1, 1])
  } 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, 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 (GLM-based) ─────────────────────────
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)))
  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>'
  )
}

# ── D2. Evaluasi SMOTE (2 data: asli + SMOTE) ────────────────────
evaluasi_smote_html <- function(model,
                                 y_true_ori, y_prob_ori,
                                 y_true_sm,  y_prob_sm,
                                 plot_id) {
  coefs <- coef(model)
  ses   <- tryCatch(sqrt(diag(vcov(model))), error = function(e) rep(NA, length(coefs)))
  pid   <- gsub("[^a-zA-Z0-9]", "_", plot_id)

  panel_ori <- paste0(
    '<div class="smote-panel active" id="panel_ori_', pid, '">',
    cm_html(y_true_ori, y_prob_ori),
    roc_html(y_true_ori, y_prob_ori, paste0(plot_id, "_ori")),
    '</div>'
  )
  panel_sm <- paste0(
    '<div class="smote-panel" id="panel_sm_', pid, '">',
    cm_html(y_true_sm, y_prob_sm),
    roc_html(y_true_sm, y_prob_sm, paste0(plot_id, "_smote")),
    '</div>'
  )
  tabs <- paste0(
    '<div class="smote-tabs">',
    '<button class="smote-tab-btn active" ',
    paste0('onclick="switchSmote(\'', pid, '\',\'ori\')">&#128202; Data Asli (Fair)</button>'),
    '<button class="smote-tab-btn" ',
    paste0('onclick="switchSmote(\'', pid, '\',\'sm\')">&#129512; Data SMOTE (In-Sample)</button>'),
    '</div>'
  )
  paste0(
    '<div class="eval-wrap">',
    separation_html(coefs, ses, FALSE),
    '<div class="eval-sub"><h6>Confusion Matrix, Metrik &amp; ROC — Pilih Data</h6>',
    tabs, panel_ori, panel_sm,
    '</div></div>'
  )
}

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

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

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

  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 '',
    '</p>')

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

  paste0(info,
    sec("📊", "EDA — Distribusi Y",               eda_html(df)),
    sec("📐", "Uji Simultan (G²)",                simultan_html(mdl, mdl0, xm)),
    sec("📋", "Uji Parsial (Wald Test)",           wald_html(mdl)),
    sec("🔗", "Multikolinearitas (VIF)",           vif_html(mdl, df, xm, "Reglog Biasa")),
    sec("✅", "Goodness of Fit (H-L Test)",        hl_html(mdl, df$status_miskin)),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC)", ic_html(mdl, n_biasa, k_biasa, "glm")),
    sec("📈", "Interpretasi Parameter (OR)",       or_html(mdl)),
    sec("🔍", "Evaluasi Model",
        evaluasi_html(mdl, df$status_miskin, fitted(mdl), pid))
  )
}

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

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

  mdl <- tryCatch(logistf(fm, data = df), error = function(e) NULL)
  if (is.null(mdl)) return("<p>⚠️ Model Firth gagal di-fit.</p>")

  # Uji Simultan — penalized LRT
  G2     <- 2 * abs(diff(mdl$loglik))
  df_G2  <- length(xm)
  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 & uji parsial
  coefs   <- mdl$coefficients
  se_raw  <- sqrt(diag(mdl$var))
  se      <- setNames(se_raw[seq_along(coefs)], names(coefs))
  pvals   <- mdl$prob[names(coefs)]
  chi_sq  <- (coefs / se)^2

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

  # Odds Ratio dengan Firth CI
  or_tbl <- kbl_h(data.frame(
    Variabel     = names(coefs),
    `exp(B)`     = round(exp(coefs), 4),
    `CI 2.5%`    = round(exp(mdl$ci.lower), 4),
    `CI 97.5%`   = round(exp(mdl$ci.upper), 4),
    Interpretasi = ifelse(names(coefs) == "(Intercept)", "Baseline odds",
                   ifelse(exp(coefs) > 1,
                          "↑ Lebih berisiko tidak bersekolah",
                          "↓ Lebih kecil risiko tidak bersekolah")),
    check.names = FALSE
  ), "Interpretasi Parameter (Odds Ratio — Firth CI)")

  info <- paste0('<p><strong>n valid:</strong> ', nrow(df),
    ' obs. | <em>Firth menangani quasi-complete separation & sampel kecil</em></p>')

  # GoF Hosmer-Lemeshow
  hl_firth <- tryCatch(
    hoslem.test(df$status_miskin, mdl$predict, 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)")
  }

  k_firth <- length(coefs)
  n_firth <- nrow(df)
  pid     <- paste0("firth-", gsub("[^a-zA-Z0-9]", "_", nama_dataset))

  paste0(info,
    sec("📊", "EDA — Distribusi Y",                     eda_html(df)),
    sec("📐", "Uji Simultan (Penalized LRT)",            sim_tbl),
    sec("📋", "Uji Per Variabel",                        wald_tbl),
    sec("🔗", "Multikolinearitas (VIF)",                 vif_html(mdl, df, xm, "Firth")),
    sec("✅", "Goodness of Fit (H-L Test)",              hl_firth_tbl),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC)",       ic_html(mdl, n_firth, k_firth, "firth")),
    sec("📈", "Interpretasi Parameter (OR — Firth CI)",  or_tbl),
    sec("🔍", "Evaluasi Model",
        evaluasi_html(mdl, df$status_miskin, mdl$predict, pid, is_firth = TRUE))
  )
}

run_smote <- function(data, nama_dataset) {
  p  <- prep_data(data)
  df <- p$df; xm <- p$xvars_model; fm <- p$formula_mod

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

  df_sm <- tryCatch({
    res <- smotefamily::SMOTE(
      X        = df[, xm, drop = FALSE],
      target   = as.character(df$status_miskin),
      K        = 5,
      dup_size = 0
    )$data
    res$status_miskin <- as.numeric(as.character(res$class))
    res$class         <- NULL
    res[, c("status_miskin", xm)]
  }, error = function(e) NULL)

  if (is.null(df_sm))
    return("<p>⚠️ SMOTE gagal — kemungkinan n terlalu kecil atau kelas sangat tidak seimbang.</p>")

  mdl  <- tryCatch(glm(fm, data = df_sm, family = binomial), error = function(e) NULL)
  mdl0 <- tryCatch(glm(status_miskin ~ 1, data = df_sm, family = binomial), error = function(e) NULL)
  if (is.null(mdl)) return("<p>⚠️ Model GLM setelah SMOTE gagal.</p>")

  # EDA before & after
  t_ori   <- table(df$status_miskin)
  t_sm    <- table(df_sm$status_miskin)
  eda_tbl <- kbl_h(data.frame(
    Kategori  = c("Tidak Bersekolah (0)", "Bersekolah (1)"),
    `N Asal`  = as.integer(c(t_ori["0"], t_ori["1"])),
    `N SMOTE` = as.integer(c(t_sm["0"], t_sm["1"])),
    check.names = FALSE
  ), "Distribusi Y Sebelum & Sesudah SMOTE")

  info <- paste0('<p><strong>n asli:</strong> ', nrow(df),
    ' | <strong>n setelah SMOTE:</strong> ', nrow(df_sm), '</p>')

  k_smote <- length(coef(mdl))
  n_smote <- nrow(df_sm)
  pid     <- paste0("smote-", gsub("[^a-zA-Z0-9]", "_", nama_dataset))

  paste0(info,
    sec("📊", "EDA — Distribusi Y (Sebelum & Sesudah SMOTE)", eda_tbl),
    sec("📐", "Uji Simultan (G²)",       simultan_html(mdl, mdl0, xm)),
    sec("📋", "Uji Parsial (Wald Test)",  wald_html(mdl)),
    sec("🔗", "Multikolinearitas (VIF)",
        paste0(vif_html(mdl, df_sm, xm, "SMOTE"),
               '<p><em><small>&#9432; VIF dihitung dari data SMOTE (in-sample). ',
               'VIF pada data asli mungkin sedikit berbeda.</small></em></p>')),
    sec("✅", "Goodness of Fit (H-L Test)", hl_html(mdl, df_sm$status_miskin)),
    sec("📊", "Kriteria Informasi (AIC/AICc/BIC)",
        paste0(ic_html(mdl, n_smote, k_smote, "glm"),
               '<p><em><small>&#9432; AIC/AICc/BIC dihitung dari data SMOTE (n = ',
               n_smote, '). Gunakan sebagai perbandingan; evaluasi fair tetap di data asli.</small></em></p>')),
    sec("📈", "Interpretasi Parameter (OR)", or_html(mdl)),
    sec("🔍", "Evaluasi Model",
        evaluasi_smote_html(
          mdl,
          df$status_miskin,    predict(mdl, newdata = df,    type = "response"),
          df_sm$status_miskin, predict(mdl, newdata = df_sm, type = "response"),
          pid
        ))
  )
}