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

# ── Visualisasi ───────────────────────────────────────────────────
library(ggplot2)
library(scales)
library(gridExtra)

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

data_susenas <- read.dbf("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) >= 10 & as.numeric(R101) < 100
  )

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),
    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
    X6 = case_when(
      R706 >= 1 & R706 <= 6     ~ 0L,
      R706 > 6  | R706 == 0     ~ 1L,
      TRUE                      ~ NA_integer_
    )
  ) %>%
  select(all_of(id_rt), X4, X5, X6)

files_43 <- c(
  "46654/Susenas KP Maret 2023/bv2_2026_02_22_08_55_02_ssn202303kp_blok43_51_94.dbf",
  "46654/Susenas KP Maret 2023/bv2_2026_02_22_08_49_01_ssn202303kp_blok43_11_31.dbf",
  "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(
  "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),    labels = c("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",
    "1 = Non-pertanian; 0 = Pertanian",
    "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)

ggplot(tbl_y, aes(x = Y, y = n, fill = Y)) +
  geom_bar(stat = "identity", width = 0.5, show.legend = FALSE) +
  geom_text(aes(label = paste0(n, "\n(", Persen, "%)")),
            vjust = -0.4, size = 4) +
  scale_fill_manual(values = c("#D73027", "#4575B4")) +
  labs(title    = "Distribusi Partisipasi Sekolah",
       subtitle = "Anak Penyandang Disabilitas Usia 13–15 Tahun\nSulawesi & Kalimantan (SUSENAS 2023)",
       x = NULL, y = "Frekuensi") +
  theme_minimal(base_size = 13)

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

plot_list <- lapply(seq_along(vars_cat), function(i) {
  var    <- vars_cat[i]
  df_tmp <- dataku_final %>%
    count(.data[[var]]) %>%
    filter(!is.na(.data[[var]])) %>%
    mutate(pct = round(n / sum(n) * 100, 1))

  ggplot(df_tmp, aes(x = .data[[var]], y = n, fill = .data[[var]])) +
    geom_bar(stat = "identity", show.legend = FALSE) +
    geom_text(aes(label = paste0(pct, "%")), vjust = -0.4, size = 3) +
    labs(title = labels_cat[i], x = NULL, y = "n") +
    theme_minimal(base_size = 10) +
    theme(plot.title = element_text(size = 10, face = "bold"))
})

do.call(grid.arrange, c(plot_list, ncol = 3))

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)

ggplot(tbl_dis, aes(x = factor(JUMLAH_JENIS_DISABILITAS), y = n)) +
  geom_bar(stat = "identity", fill = "#4575B4", width = 0.6) +
  geom_text(aes(label = paste0(n, "\n(", Persen, "%)")),
            vjust = -0.3, size = 3.5) +
  labs(title = "Distribusi Jumlah Jenis Disabilitas",
       x = "Jumlah Jenis Disabilitas", y = "Frekuensi") +
  theme_minimal(base_size = 12)

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)

ggplot(tbl_prov, aes(x = reorder(factor(prov), n), y = n)) +
  geom_bar(stat = "identity", fill = "#74ADD1") +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title = "Jumlah Observasi per Provinsi",
       x = "Kode Provinsi", y = "Frekuensi") +
  theme_minimal(base_size = 12)

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

plot_stacked <- function(var, label) {
  df_tmp <- dataku_final %>%
    filter(!is.na(.data[[var]]), !is.na(Y)) %>%
    count(.data[[var]], Y) %>%
    group_by(.data[[var]]) %>%
    mutate(pct = round(n / sum(n) * 100, 1)) %>%
    ungroup()

  ggplot(df_tmp, aes(x = .data[[var]], y = pct, fill = Y)) +
    geom_bar(stat = "identity", position = "stack") +
    geom_text(aes(label = paste0(pct, "%")),
              position = position_stack(vjust = 0.5),
              size = 3, color = "white", fontface = "bold") +
    scale_fill_manual(values = c("Tidak bersekolah" = "#D73027",
                                 "Bersekolah"        = "#4575B4")) +
    labs(title = label, x = NULL, y = "%", fill = NULL) +
    theme_minimal(base_size = 10) +
    theme(legend.position = "bottom",
          plot.title      = element_text(size = 10, face = "bold"))
}

p_list <- lapply(seq_along(vars_for_tab), function(i) {
  plot_stacked(vars_for_tab[i], label_vars[i])
})

do.call(grid.arrange, c(p_list, ncol = 3))

df_prov_y <- dataku_final %>%
  filter(!is.na(Y)) %>%
  count(prov, Y) %>%
  group_by(prov) %>%
  mutate(pct = round(n / sum(n) * 100, 1)) %>%
  ungroup()

ggplot(df_prov_y, aes(x = factor(prov), y = pct, fill = Y)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(label = paste0(pct, "%")),
            position = position_stack(vjust = 0.5),
            size = 2.8, color = "white", fontface = "bold") +
  scale_fill_manual(values = c("Tidak bersekolah" = "#D73027",
                               "Bersekolah"        = "#4575B4")) +
  labs(title    = "Partisipasi Sekolah Anak Disabilitas Usia 13–15 Tahun\nper Provinsi (Sulawesi & Kalimantan)",
       x = "Kode Provinsi", y = "%", fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

df_pulau_y <- dataku_final %>%
  filter(!is.na(Y)) %>%
  mutate(
    pulau = case_when(
      floor(prov / 10) == 6 ~ "Kalimantan",
      floor(prov / 10) == 7 ~ "Sulawesi",
      TRUE                   ~ NA_character_
    ),
    pulau = factor(pulau, levels = c("Kalimantan", "Sulawesi"))
  ) %>%
  filter(!is.na(pulau)) %>%
  count(pulau, Y) %>%
  group_by(pulau) %>%
  mutate(
    total = sum(n),
    pct   = round(n / total * 100, 1)
  ) %>%
  ungroup()

ggplot(df_pulau_y, aes(x = pulau, y = pct, fill = Y)) +
  geom_bar(stat = "identity", position = "stack", width = 0.65) +
  geom_text(aes(label = paste0(pct, "%")),
            position = position_stack(vjust = 0.5),
            size = 3.5, color = "white", fontface = "bold") +
  geom_text(
    data = df_pulau_y %>% distinct(pulau, total),
    aes(x = pulau, y = 102, label = paste0("n=", total)),
    inherit.aes = FALSE, size = 3.2, color = "grey40"
  ) +
  scale_fill_manual(values = c("Tidak bersekolah" = "#D73027",
                               "Bersekolah"        = "#4575B4")) +
  scale_y_continuous(limits = c(0, 110), labels = function(x) paste0(x, "%")) +
  labs(
    title    = "Partisipasi Sekolah Anak Penyandang Disabilitas Usia 13–15 Tahun",
    subtitle = "per Pulau — SUSENAS 2023 (Sulawesi & Kalimantan)",
    x        = NULL,
    y        = "Persentase (%)",
    fill     = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "bottom",
    plot.title      = element_text(face = "bold"),
    axis.text.x     = element_text(size = 11)
  )

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
    
  )

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),
    rep("Nusa Tenggara & Bali", 3),
    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")

library(dplyr)
library(ggplot2)
library(tidyr)

# ── 1. Restrukturisasi Data untuk Faceting Vertikal ───────────────────────
plot_data_ind <- outputP_ind %>%
  select(Nama_Prov, Wilayah, Estimasi, RSE, CI_LOWER, CI_UPPER) %>%
  mutate(
    # Urutkan provinsi berdasarkan nilai Estimasi dari terbesar ke terkecil
    Nama_Prov = factor(Nama_Prov, levels = Nama_Prov[order(-Estimasi)]),
    RSE_flag  = factor(
      ifelse(RSE > 25, "RSE > 25% (Tidak Reliabel)", "RSE ≤ 25% (Reliabel)"),
      levels = c("RSE ≤ 25% (Reliabel)", "RSE > 25% (Tidak Reliabel)")
    )
  ) %>%
  # Mengubah data menjadi format long agar bisa di-facet secara vertikal
  pivot_longer(cols = c(Estimasi, RSE), names_to = "Indikator", values_to = "Nilai")

# ── 2. Pembuatan Single Plot dengan Facet ─────────────────────────────────
ggplot(plot_data_ind, aes(x = Nama_Prov, y = Nilai, fill = Wilayah)) +
  # Batang utama grafik
  geom_bar(stat = "identity", width = 0.65, alpha = 0.88, color = "white", linewidth = 0.2) +
  
  # Error Bar Khusus untuk panel Estimasi saja
  geom_errorbar(
    data = filter(plot_data_ind, Indikator == "Estimasi"),
    aes(ymin = CI_LOWER, ymax = CI_UPPER),
    width = 0.2, color = "gray20", linewidth = 0.6
  ) +
  
  # Garis Batas RSE 25% Khusus untuk panel RSE saja
  geom_hline(
    data = filter(plot_data_ind, Indikator == "RSE"),
    aes(yintercept = 25), linetype = "dashed", color = "#C62828", linewidth = 0.7
  ) +
  
  # Label angka di atas bar (fleksibel mengikuti posisi CI_UPPER atau Nilai RSE)
  geom_text(
    aes(
      label = paste0(round(Nilai, 1), "%"),
      y = ifelse(Indikator == "Estimasi", CI_UPPER, Nilai)
    ),
    vjust = -0.5, size = 3.3, fontface = "bold", color = "gray10"
  ) +
  
  # Pewarnaan khusus untuk wilayah Kalimantan & Sulawesi
  scale_fill_manual(values = c("Kalimantan" = "#1976D2", "Sulawesi" = "#E64A19")) +
  
  # Pengaturan skala Y otomatis melonggar di atas agar text tidak terpotong
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  
  # Pembagian Panel Vertikal (Atas: Estimasi, Bawah: RSE) dengan skala Y bebas
  facet_wrap(~ Indikator, ncol = 1, scales = "free_y", 
             labeller = as_labeller(c("Estimasi" = "Estimasi Angka Sekolah (95% CI)", 
                                     "RSE" = "Relative Standard Error / RSE (Batas 25%)"))) +
  
  # Labeling Grafis
  labs(
    title    = "Analisis Indikator Disabilitas Usia 13–15 Tahun yang Bersekolah",
    subtitle = "Sulawesi & Kalimantan — SUSENAS 2023 | Estimasi Langsung",
    x        = NULL,
    y        = "Persentase (%)",
    fill     = "Wilayah:",
    caption  = "Garis putus-putus merah menyatakan batas maksimal reliabilitas RSE BPS (25%)"
  ) +
  
  # Tema visual minimalis & profesional
  theme_minimal(base_size = 12) +
  theme(
    plot.title         = element_text(face = "bold", size = 14, color = "gray10"),
    plot.subtitle      = element_text(color = "gray40", size = 10, margin = margin(b = 12)),
    
    # Merapikan text nama provinsi di sumbu X agar vertikal/miring dan tidak bertumpuk
    axis.text.x        = element_text(angle = 45, hjust = 1, vjust = 1, face = "bold", color = "gray20"),
    axis.title.y       = element_text(size = 10, color = "gray30"),
    
    # Desain strip panel facet yang clean
    strip.background   = element_rect(fill = "gray95", color = NA),
    strip.text         = element_text(face = "bold", size = 11, color = "gray20", margin = margin(t = 6, b = 6)),
    
    # Menghapus grid line yang tidak perlu
    panel.grid.major.x = element_blank(),
    panel.grid.minor   = element_blank(),
    
    legend.position    = "bottom",
    legend.margin      = margin(t = -5)
  )

## 1. AGREGASI DATA PER PULAU (Menggabungkan Nilai Varians & SE)
# Menggabungkan data provinsi ke tingkat pulau dengan pendekatan survey
pulau_data <- aggregate(cbind(Estimasi, VAR) ~ Wilayah, data = outputP_ind, FUN = mean)

# Menghitung ulang SE, CI, dan RSE tingkat Pulau agar akurat secara statistik
pulau_data$SE       <- sqrt(pulau_data$VAR)
pulau_data$CI_LOWER <- pulau_data$Estimasi - (1.96 * pulau_data$SE)
pulau_data$CI_UPPER <- pulau_data$Estimasi + (1.96 * pulau_data$SE)
pulau_data$RSE      <- (pulau_data$SE / pulau_data$Estimasi) * 100

# Urutkan pulau berdasarkan nilai Estimasi terbesar ke terkecil agar plot rapi
pulau_data <- pulau_data[order(-pulau_data$Estimasi), ]


## 2. SETTING LAYOUT GRAFIK (Base R)
# Membagi layar menjadi 2 kolom (Kiri: Estimasi, Kanan: RSE)
old_par <- par(mfrow = c(1, 2), mar = c(6, 5, 4, 2))

# Palette warna estetik modern untuk masing-masing pulau
warna_pulau <- c("#1976D2", "#E64A19", "#388E3C", "#0097A7", "#7B1FA2", "#FBC02D", "#E91E63")[1:nrow(pulau_data)]


## ── PLOT 1: ESTIMASI LENGKAP DENGAN CI PER PULAU ──────────────────────────
y_max_est <- max(pulau_data$CI_UPPER, na.rm = TRUE) * 1.15

bp1 <- barplot(pulau_data$Estimasi,
               names.arg = pulau_data$Wilayah,
               las = 2,                  # Nama pulau tegak lurus/vertikal
               col = warna_pulau,
               border = NA,
               ylim = c(0, y_max_est),
               ylab = "Estimasi (%)",
               main = "Estimasi Angka Sekolah Per Pulau",
               font.main = 2, cex.main = 1.1, cex.names = 0.9)

# Tambahkan Grid & Gambar Ulang Bar agar bersih
grid(nx = NA, ny = NULL, col = "gray93", lty = "solid")
barplot(pulau_data$Estimasi, col = warna_pulau, border = NA, add = TRUE, las = 2, names.arg = NULL)

# Tambahkan Error Bar (Confidence Interval)
arrows(x0 = bp1, y0 = pulau_data$CI_LOWER, 
       x1 = bp1, y1 = pulau_data$CI_UPPER, 
       code = 3, angle = 90, length = 0.05, col = "gray20", lwd = 1.5)

# Tampilkan angka persentase di atas bar
text(x = bp1, y = pulau_data$Estimasi + (y_max_est * 0.02), 
     labels = paste0(round(pulau_data$Estimasi, 1), "%"), 
     cex = 0.85, font = 2)


## ── PLOT 2: RELATIVE STANDARD ERROR (RSE) PER PULAU ──────────────────────
y_max_rse <- max(c(pulau_data$RSE, 30), na.rm = TRUE) * 1.15

# Warna RSE otomatis: Hijau jika <= 25% (Reliabel), Merah jika > 25%
warna_rse <- ifelse(pulau_data$RSE > 25, "#EF5350", "#66BB6A")

bp2 <- barplot(pulau_data$RSE,
               names.arg = pulau_data$Wilayah,
               las = 2,
               col = warna_rse,
               border = NA,
               ylim = c(0, y_max_rse),
               ylab = "RSE (%)",
               main = "Nilai RSE (Reliabilitas Data) Per Pulau",
               font.main = 2, cex.main = 1.1, cex.names = 0.9)

# Tambahkan Grid & Gambar Ulang Bar RSE
grid(nx = NA, ny = NULL, col = "gray93", lty = "solid")
barplot(pulau_data$RSE, col = warna_rse, border = NA, add = TRUE, las = 2, names.arg = NULL)

# Tambahkan Garis Batas RSE 25% (Standard BPS untuk data reliabel)
abline(h = 25, col = "#C62828", linetype = "dashed", lwd = 1.5)
text(x = max(bp2)*0.3, y = 26.5, labels = "Batas RSE 25%", col = "#C62828", cex = 0.8, font = 3)

# Tampilkan nilai RSE di atas bar
text(x = bp2, y = pulau_data$RSE + (y_max_rse * 0.02), 
     labels = paste0(round(pulau_data$RSE, 1), "%"), 
     cex = 0.85, font = 2)

# Reset Layout Parameter R
par(old_par)

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

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

  n  <- nrow(df)
  n1 <- sum(df$status_miskin == 1)
  n0 <- sum(df$status_miskin == 0)

  # Simpan bobot sebagai kolom di df agar glm selalu menemukannya
  df$.w <- ifelse(df$status_miskin == 1, n / (2 * n1), n / (2 * n0))

  ctrl <- glm.control(maxit = 200, epsilon = 1e-8)

  mdl <- suppressWarnings(tryCatch(
    glm(fm, data = df, family = quasibinomial(link = "logit"),
        weights = .w, control = ctrl),
    error = function(e) NULL
  ))

  mdl0 <- suppressWarnings(tryCatch(
    glm(status_miskin ~ 1, data = df,
        family  = quasibinomial(link = "logit"),
        weights = .w, control = ctrl),
    error = function(e) NULL
  ))

  if (is.null(mdl)) return("<p>⚠️ Model Class Weight gagal di-fit.</p>")

  # EDA + weight info
  t   <- table(df$status_miskin)
  pct <- round(prop.table(t) * 100, 2)
  eda_tbl <- kbl_h(data.frame(
    Kategori       = c("Tidak Bersekolah (0)", "Bersekolah (1)"),
    N              = as.integer(c(t["0"], t["1"])),
    Persen         = paste0(c(pct["0"], pct["1"]), " %"),
    `Class Weight` = c(round(n / (2 * n0), 4), round(n / (2 * n1), 4)),
    check.names    = FALSE
  ), "Distribusi Y + Class Weight")

  info <- paste0(
    '<p><strong>n valid:</strong> ', nrow(df), ' obs.',
    ' | <strong>w bersekolah:</strong> ', round(n / (2 * n1), 3),
    ' | <strong>w tidak bersekolah:</strong> ', round(n / (2 * n0), 3),
    ' | <em>family: quasibinomial</em></p>'
  )

  k_weight <- length(coef(mdl))
  n_weight <- nrow(df)
  pid      <- paste0("weight-", gsub("[^a-zA-Z0-9]", "_", nama_dataset))

  paste0(info,
    sec("📊", "EDA — Distribusi Y + Class Weight",    eda_tbl),
    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, "Class Weight")),
    sec("✅", "Goodness of Fit (H-L Test)",            hl_html(mdl, df$status_miskin)),
    sec("📊", "Kriteria Informasi (QAIC/Deviance)",    ic_html(mdl, n_weight, k_weight, "quasi")),
    sec("📈", "Interpretasi Parameter (OR)",           or_html(mdl)),
    sec("🔍", "Evaluasi Model",
        evaluasi_html(mdl, df$status_miskin, fitted(mdl), pid))
  )
}