# ── 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(survey)            # estimasi tertimbang & confidence interval

# ── 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 (0 = Miskin/ref; 1 = Tidak Miskin)
    X9 = if_else(KAPITA < GK, 0L, 1L, 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 (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
    )
  ) %>%
  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 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", 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",
    "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"
  )
) %>%
  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") %>%
  footnote(
    general = "* = kategori referensi.",
    general_title = "Keterangan: ",
    footnote_as_chunk = TRUE,
    escape = FALSE
  )

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",
           height = PIE_CHART_HEIGHT,
           style = list(fontFamily = "inherit")) %>%
  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(
    size = PIE_SIZE,
    center = PIE_CENTER,
    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",
               sourceWidth = 800, sourceHeight = 500,
               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(),
      pct_tidak_bersekolah = round(sum(Y == "Tidak bersekolah") / n * 100, 1),
      pct_bersekolah       = round(sum(Y == "Bersekolah") / n * 100, 1),
      .groups              = "drop"
    )
}) %>% bind_rows()

kable(ringkasan,
      caption   = "Ringkasan % Status Sekolah per Kategori Prediktor",
      col.names = c("Variabel", "Kategori", "n", "% Tidak 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 dikonversi ke 0/1 (0 = tidak bersekolah; 1 = bersekolah).
# Seluruh X dipertahankan sebagai faktor agar kategori referensi mengikuti
# level pertama/faktor berkode 0 yang ditandai (*) pada tabel referensi.
prep_data <- function(data) {
  df <- data %>%
    select(all_of(c("Y", xvars_global))) %>%
    mutate(
      status_sekolah = if_else(Y == "Bersekolah", 1L, 0L),
      across(all_of(xvars_global), ~ droplevels(as.factor(.x)))
    ) %>%
    select(-Y) %>%
    na.omit()

  konstan <- xvars_global[sapply(xvars_global, function(v) dplyr::n_distinct(df[[v]]) < 2)]
  xvars_model <- setdiff(xvars_global, konstan)
  df <- df %>% mutate(across(all_of(xvars_model), ~ droplevels(as.factor(.x))))
  formula_mod <- as.formula(paste("status_sekolah ~", 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_sekolah)
  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) {
  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))
  kbl_h(data.frame(
    Variabel     = names(or),
    `exp(B)`     = round(or, 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)")
}

# ── 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_sekolah ~ 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("✅", "Goodness of Fit (H-L Test)",        hl_html(mdl, df$status_sekolah)),
    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_sekolah, fitted(mdl), pid))
  )
}

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

  # SMOTE memerlukan input numerik. Faktor X diubah menjadi dummy dengan
  # model.matrix sehingga referensi tetap kategori berkode 0/*.
  X_dummy <- model.matrix(fm, data = df)[, -1, drop = FALSE]
  X_dummy <- as.data.frame(X_dummy)
  names(X_dummy) <- make.names(names(X_dummy), unique = TRUE)
  xm_sm <- names(X_dummy)
  df_eval <- cbind(status_sekolah = df$status_sekolah, X_dummy)
  fm_sm <- as.formula(paste("status_sekolah ~", paste(xm_sm, collapse = " + ")))

  df_sm <- tryCatch({
    res <- smotefamily::SMOTE(
      X        = X_dummy,
      target   = as.character(df$status_sekolah),
      K        = 5,
      dup_size = 0
    )$data
    res$status_sekolah <- as.numeric(as.character(res$class))
    res$class         <- NULL
    res[, c("status_sekolah", xm_sm)]
  }, 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_sm, 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)) return("<p>⚠️ Model GLM setelah SMOTE gagal.</p>")

  # EDA before & after
  t_ori   <- table(df$status_sekolah)
  t_sm    <- table(df_sm$status_sekolah)
  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),
    ' | <em>X kategorik dikonversi menjadi dummy; referensi tetap kategori bertanda (*).</em></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_sm)),
    sec("📋", "Uji Parsial (Wald Test)",  wald_html(mdl)),
    sec("✅", "Goodness of Fit (H-L Test)", hl_html(mdl, df_sm$status_sekolah)),
    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_eval$status_sekolah, predict(mdl, newdata = df_eval, type = "response"),
          df_sm$status_sekolah,   predict(mdl, newdata = df_sm,   type = "response"),
          pid
        ))
  )
}