# Load urutan: MASS & car DULU, baru dplyr (cegah masking select/recode)
library(MASS)
library(car)
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(sae)        # mseFH
library(glmmTMB)    # GLMM binomial
library(classInt)   # Jenks Natural Breaks
library(factoextra)
library(knitr)
library(kableExtra)

# Ikat namespace kritis
select <- dplyr::select
recode <- dplyr::recode

has_jags <- tryCatch({
  library(rjags); library(saeHB); TRUE
}, error = function(e) {
  message("saeHB/JAGS tidak tersedia — skenario HB Beta akan fallback ke direct.")
  FALSE
})
cat("Status JAGS/saeHB:", ifelse(has_jags, "Tersedia", "Tidak tersedia"), "\n")

Estimasi_531  <- read_excel("C:/Users/User/Downloads/Hasil_sulawesi_KabKota.xlsx")
Podes_kab2024 <- read_excel("C:/Users/User/Downloads/Podes_kab_lengkap.xlsx")
Data_Penduduk <- read_excel("C:/Users/User/Downloads/Data_Penduduk.xlsx")

Estimasi_531  <- Estimasi_531  %>% mutate(Kako = sprintf("%04d", Kako))
Podes_kab2024 <- Podes_kab2024 %>% mutate(kode_kab = sprintf("%04s", kode_kab))

df_gabungan <- Estimasi_531 %>%
  left_join(Podes_kab2024, by = c("Kako" = "kode_kab")) %>%
  dplyr::select(Kako, Estimasi, RSE, starts_with("X"))

df_gabungan$Kako <- as.numeric(df_gabungan$Kako)
df_final <- df_gabungan %>% left_join(Data_Penduduk, by = "Kako")

cat("Jumlah domain:", nrow(df_final), "\n")

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

df_final <- df_final %>%
  mutate(
    X6  = X6_jumlah  / Jumlah_Penduduk * 1000,
    X7  = X7_jumlah  / Jumlah_Penduduk * 1000,
    X10 = X10_jumlah / Jumlah_Penduduk * 1000,
    X33 = X33_jumlah / Jumlah_Penduduk * 1000,
    X11 = X11_jumlah / Jumlah_P * 1000,
    X12 = X12_jumlah / Jumlah_P * 1000,
    X27 = X27_jumlah / Jumlah_P * 1000,
    X28 = X28_jumlah / Jumlah_P * 1000,
    X29 = X29_jumlah / Jumlah_P * 1000,
    X30 = X30_jumlah / Jumlah_P * 1000,
    X31 = X31_jumlah / Jumlah_P * 1000,
    X35 = X35_jumlah / Jumlah_P * 1000,
    X36 = X36_jumlah / Jumlah_P * 1000
  )

eps_logit <- 1e-4

df_base <- df_final %>%
  mutate(
    vardir  = (RSE / 100 * Estimasi)^2,
    y_prop  = pmax(pmin(Estimasi / 100, 1 - eps_logit), eps_logit),
    y_logit = log(y_prop / (1 - y_prop))
  ) %>%
  as.data.frame()

# ── Kandidat aux: HANYA kolom yang diawali "X" ───────────────────────────
# RSE_direct_col, grup_jenks, cluster_*, dll yang ditambahkan kemudian
# otomatis TIDAK masuk karena tidak berawalan X
kandidat_vars <- names(df_base)[startsWith(names(df_base), "X")]

# ── excl_always: otomatis semua kolom SELAIN X* ──────────────────────────
# Diupdate di setiap fungsi seleksi agar mencakup kolom baru yang
# ditambahkan setelah chunk ini (RSE_direct_col, grup_*, cluster_*)
excl_always <- setdiff(names(df_base), kandidat_vars)

cat("Variabel kandidat (", length(kandidat_vars), "):\n",
    paste(kandidat_vars, collapse = ", "), "\n")

cat("\nTotal dikecualikan:", length(excl_always), "kolom\n")

# Helper: pilih 1 representatif terbaik per grup X (X19, X19_mean, X19_jumlah → 1)
# Pemilihan berdasarkan |cor| tertinggi ke y_col pada data segmen berjalan
pilih_per_grup <- function(df, kandidat, y_col) {
  get_base <- function(x) sub("(_jumlah|_mean)$", "", x)
  bases    <- unique(sapply(kandidat, get_base))
  terpilih <- c()
  for (b in bases) {
    anggota <- kandidat[get_base(kandidat) == b]
    if (length(anggota) == 1) { terpilih <- c(terpilih, anggota); next }
    cor_v   <- sapply(anggota, function(v)
      abs(cor(df[[v]], df[[y_col]], use = "complete.obs")))
    cor_v   <- cor_v[!is.na(cor_v)]
    if (length(cor_v) == 0) { terpilih <- c(terpilih, anggota[1]); next }
    terpilih <- c(terpilih, names(which.max(cor_v)))
  }
  unique(terpilih)
}

seleksi_backward <- function(df, y_col, excl_cols = NULL, r2_cap = 0.8) {
  if (is.null(excl_cols)) excl_cols <- names(df)[!startsWith(names(df), "X")]
  n_domain <- nrow(df)
  max_top  <- max(1L, floor(n_domain / 10L))

  keep <- sapply(df, function(x) length(unique(na.omit(x))) > 1)
  df   <- df[, keep]

  kandidat <- setdiff(names(df), excl_cols)
  kandidat <- kandidat[startsWith(kandidat, "X")]
  if (length(kandidat) == 0) return(character(0))

  # Pilih 1 representatif per grup X (sesuai y_col & data segmen ini)
  kandidat <- pilih_per_grup(df, kandidat, y_col)
  if (length(kandidat) == 0) return(character(0))

  cor_v <- sapply(df[kandidat], function(x)
    abs(cor(x, df[[y_col]], use = "complete.obs")))
  cor_v <- cor_v[!is.na(cor_v)]
  if (length(cor_v) == 0) return(character(0))
  vars  <- names(sort(cor_v, decreasing = TRUE))[1:min(max_top, length(cor_v))]
  vars  <- vars[!is.na(vars)]

  step_m <- try(stepAIC(
    lm(as.formula(paste(y_col, "~", paste(vars, collapse = "+"))), data = df),
    direction = "backward", trace = FALSE), silent = TRUE)
  if (!inherits(step_m, "try-error")) {
    sv <- names(coef(step_m))[-1]
    if (length(sv) > 0) vars <- sv
  }

  repeat {
    if (length(vars) <= 1) break
    m <- try(lm(as.formula(paste(y_col,"~",paste(vars,collapse="+"))),data=df),silent=TRUE)
    if (inherits(m,"try-error")) { vars <- vars[1]; break }
    v <- try(vif(m), silent=TRUE)
    if (inherits(v,"try-error")) break
    v_ok <- v[!is.na(v)]
    if (length(v_ok)==0 || all(v_ok<=10)) break
    vars <- vars[vars != names(which.max(v_ok))]
  }

  if (length(vars) > 1) {
    m <- lm(as.formula(paste(y_col,"~",paste(vars,collapse="+"))),data=df)
    if (summary(m)$r.squared > r2_cap) vars <- vars[1]
  }
  return(vars)
}

seleksi_topn <- function(df, y_col, excl_cols = NULL, r2_cap = 0.8) {
  if (is.null(excl_cols)) excl_cols <- names(df)[!startsWith(names(df), "X")]
  n_domain <- nrow(df)
  max_top  <- min(max(1L, floor(n_domain / 10L)), 8L)

  keep <- sapply(df, function(x) length(unique(na.omit(x))) > 1)
  df   <- df[, keep]

  kandidat <- setdiff(names(df), excl_cols)
  kandidat <- kandidat[startsWith(kandidat, "X")]
  if (length(kandidat) == 0) return(character(0))

  # Pilih 1 representatif per grup X
  kandidat <- pilih_per_grup(df, kandidat, y_col)
  if (length(kandidat) == 0) return(character(0))

  cor_v <- sapply(df[kandidat], function(x)
    abs(cor(x, df[[y_col]], use = "complete.obs")))
  cor_v <- cor_v[!is.na(cor_v)]
  if (length(cor_v) == 0) return(character(0))
  vars  <- names(sort(cor_v, decreasing = TRUE))[1:min(max_top, length(cor_v))]
  vars  <- vars[!is.na(vars)]

  repeat {
    if (length(vars) <= 1) break
    m <- try(lm(as.formula(paste(y_col,"~",paste(vars,collapse="+"))),data=df),silent=TRUE)
    if (inherits(m,"try-error")) { vars <- vars[1]; break }
    v <- try(vif(m), silent=TRUE)
    if (inherits(v,"try-error")) break
    v_ok <- v[!is.na(v)]
    if (length(v_ok)==0 || all(v_ok<=10)) break
    vars <- vars[vars != names(which.max(v_ok))]
  }

  if (length(vars) > 1) {
    m <- lm(as.formula(paste(y_col,"~",paste(vars,collapse="+"))),data=df)
    if (summary(m)$r.squared > r2_cap) vars <- vars[1]
  }
  return(vars)
}

run_eblup <- function(df, vars) {
  df$vardir <- pmax(df$vardir, 1e-10)

  # Trim ke kolom yang dibutuhkan mseFH saja (seperti pola SAE_Lengkap_Fixed2)
  # Kolom ekstra (grup_jenks, cluster_*, RSE_direct_col, dll) dikeluarkan
  # agar mseFH tidak terganggu
  df_m        <- df[, c("Estimasi", "vardir", vars), drop = FALSE]

  # Scale per kolom — lebih aman dari scale() bulk untuk kolom near-constant
  for (v in vars) {
    mu <- mean(df_m[[v]], na.rm = TRUE)
    sd <- sd(df_m[[v]],   na.rm = TRUE)
    df_m[[v]] <- if (is.na(sd) || sd == 0) 0 else (df_m[[v]] - mu) / sd
    df_m[[v]][is.na(df_m[[v]]) | is.nan(df_m[[v]])] <- 0
  }

  form <- as.formula(paste("Estimasi ~", paste(vars, collapse = "+")))
  m    <- try(mseFH(form, vardir = vardir, data = df_m), silent = TRUE)

  if (inherits(m, "try-error") || is.null(m$est) || is.null(m$mse)) {
    df$y_eblup <- df$Estimasi
    df$mse     <- df$vardir
  } else {
    ev  <- as.numeric(m$est$eblup)
    mv  <- as.numeric(m$mse)
    # Per-domain fallback: domain NaN/NA → pakai direct, bukan fallback seluruh segmen
    bad <- is.na(ev) | is.nan(ev) | is.na(mv) | is.nan(mv)
    if (any(bad)) {
      message("  EBLUP: ", sum(bad), " domain NaN → fallback direct per domain")
      ev[bad] <- df$Estimasi[bad]
      mv[bad] <- df$vardir[bad]
    }
    df$y_eblup <- ev
    df$mse     <- mv
  }

  df$mse        <- pmax(df$mse, 0)
  df$RSE_direct <- sqrt(df$vardir) / pmax(df$Estimasi, 1e-6) * 100
  df$RMSE_eblup <- sqrt(df$mse)
  df$RSE_eblup  <- df$RMSE_eblup / pmax(abs(df$y_eblup), 1e-6) * 100
  df
}

run_glmm <- function(df, vars) {
  eps       <- 1e-4
  df$y_prop <- pmax(pmin(df$Estimasi / 100, 1 - eps), eps)
  se_prop   <- pmax(df$RSE / 100 * df$y_prop, eps)

  # n_eff minimum 5 agar binomial tidak quasi-saturated
  df$n_eff   <- pmax(as.integer(round(df$y_prop * (1 - df$y_prop) / se_prop^2)), 5L)
  df$y_count <- pmin(as.integer(round(df$y_prop * df$n_eff)), df$n_eff - 1L)
  df$y_count <- pmax(df$y_count, 1L)   # hindari 0 dan n_eff (separasi)

  df_m        <- df[, c("Kako", "y_count", "n_eff", vars)]
  df_m[vars]  <- as.data.frame(scale(df_m[vars]))
  df_m[vars][is.nan(as.matrix(df_m[vars]))] <- 0
  df_m$Kako   <- as.factor(df_m$Kako)

  form <- as.formula(paste(
    "cbind(y_count, n_eff - y_count) ~",
    paste(vars, collapse = " + "),
    "+ (1 | Kako)"
  ))

  ctrl <- glmmTMBControl(optCtrl = list(iter.max = 500, eval.max = 800))
  m    <- try(glmmTMB(form, data = df_m, family = binomial, control = ctrl),
              silent = TRUE)

  if (inherits(m, "try-error") || is.null(m)) {
    message("  GLMM fit gagal → fallback direct")
    df$y_glmm   <- df$Estimasi
    df$mse_glmm <- df$vardir
  } else {
    # Prediksi: fixed + random effect (sesuai slide: X_i*beta + v_i_hat)
    df$y_glmm <- plogis(predict(m, type = "link", re.form = NULL)) * 100

    # ── Jackknife MSE (M2i saja) — M1i diabaikan karena g1i GLMM tidak
    #    memiliki bentuk analytik seperti FH, jadi full MSE ≈ M2i
    n_d   <- nrow(df)
    theta <- df$y_glmm   # estimasi penuh

    mse_v <- tryCatch({
      jack_pred <- matrix(NA_real_, nrow = n_d, ncol = n_d)

      for (i in seq_len(n_d)) {
        df_j <- df_m[-i, ]
        m_j  <- try(glmmTMB(form, data = df_j, family = binomial,
                             control = glmmTMBControl(
                               optCtrl = list(iter.max = 300, eval.max = 500))),
                    silent = TRUE)
        if (!inherits(m_j, "try-error")) {
          pj <- try(plogis(predict(m_j, newdata = df_m, type = "link",
                                    re.form = NA,
                                    allow.new.levels = TRUE)) * 100,
                    silent = TRUE)
          if (!inherits(pj, "try-error") && length(pj) == n_d)
            jack_pred[i, ] <- pj
        }
      }

      # M2i = (m-1)/m * sum(theta_i,-l - theta_i)^2
      vapply(seq_len(n_d), function(j) {
        ev <- jack_pred[, j][!is.na(jack_pred[, j])]
        if (length(ev) >= 2) {
          ((length(ev) - 1) / length(ev)) * mean((ev - theta[j])^2)
        } else {
          (df$RSE[j] / 100 * theta[j])^2
        }
      }, numeric(1))

    }, error = function(e) {
      message("  Jackknife crash: ", e$message, " → fallback MSE empiris")
      (df$RSE / 100 * df$y_glmm)^2
    })

    df$mse_glmm <- pmax(mse_v, 0)
  }

  df$RSE_direct <- df$RSE
  df$RMSE_glmm  <- sqrt(df$mse_glmm)
  df$RSE_glmm   <- df$RMSE_glmm / pmax(df$y_glmm, 1e-6) * 100
  df
}

#' EB Beta-Binomial dengan covariate:
#' (1) Regresi logit(p) ~ vars → mu_i (domain-specific prior mean)
#' (2) Estimasi phi (precision) via MOM dari {p_i}
#' (3) EB: p_hat_i = (a_i + y_i) / (a_i + b_i + n_eff_i) x 100
#' (4) Var posterior analytik → RSE_eb
run_eb_beta <- function(df, vars) {
  eps       <- 1e-4
  df$y_prop <- pmax(pmin(df$Estimasi / 100, 1 - eps), eps)

  se_prop    <- pmax(df$RSE / 100 * df$y_prop, eps)
  df$n_eff   <- pmax(as.integer(round(df$y_prop * (1 - df$y_prop) / se_prop^2)), 2L)
  df$y_count <- pmin(as.integer(round(df$y_prop * df$n_eff)), df$n_eff)

  df$y_logit_eb <- log(df$y_prop / (1 - df$y_prop))

  df[vars] <- as.data.frame(scale(df[vars]))
  df[vars][is.nan(as.matrix(df[vars]))] <- 0

  # Regresi untuk mu_i (prior mean per domain)
  form_reg <- as.formula(paste("y_logit_eb ~", paste(vars, collapse = "+")))
  reg      <- try(lm(form_reg, data = df), silent = TRUE)
  mu_i     <- if (inherits(reg, "try-error")) {
    df$y_prop
  } else {
    pmax(pmin(plogis(fitted(reg)), 1 - eps), eps)
  }

  # MOM untuk phi (precision)
  p_bar  <- mean(df$y_prop)
  var_p  <- var(df$y_prop)
  phi    <- max((p_bar * (1 - p_bar) / var_p) - 1, 0.01)

  a_i <- mu_i * phi
  b_i <- (1 - mu_i) * phi
  n_i <- df$n_eff
  y_i <- df$y_count

  # EB estimator
  a_post   <- a_i + y_i
  b_post   <- b_i + (n_i - y_i)
  denom    <- a_i + b_i + n_i
  p_hat_eb <- a_post / denom

  df$y_eb_pct <- p_hat_eb * 100

  # Var posterior analytik (Beta conjugate)
  var_post      <- (a_post * b_post) / (denom^2 * (denom + 1))
  df$var_eb_pct <- var_post * 100^2

  df$RSE_direct <- df$RSE
  df$RMSE_eb    <- sqrt(df$var_eb_pct)
  df$RSE_eb     <- df$RMSE_eb / pmax(df$y_eb_pct, 1e-6) * 100
  df
}

# ── Helper: identifikasi kovariat yang CI-nya menyeberang nol ────────────────
# "Tidak aman" := 2.5% < 0 DAN 97.5% > 0  (intercept selalu dikecualikan)
#
# Desain: akses POSISIONAL, bukan berbasis nama kolom/baris.
# saeHB bisa menyimpan rownames sebagai "beta[1]","beta[2]"... (notasi JAGS)
# meski di-print tampak seperti nama variabel — intersect() gagal dalam kasus itu.
#
# Struktur baku saeHB::Beta()$coefficient:
#   Baris 1     = intercept  (di-skip)
#   Baris 2..k+1 = kovariat urutan sesuai formula = urutan var_names
#   Kolom 1=Mean | 2=SD | 3=2.5% | 4=25% | 5=50% | 6=75% | 7=97.5%
#
.vars_cross_ci <- function(coef_mat, var_names) {
  if (is.null(coef_mat) || nrow(coef_mat) < 2L || ncol(coef_mat) < 7L)
    return(character(0))

  # Baris kovariat: lewati baris-1 (intercept), petakan ke var_names
  p    <- min(length(var_names), nrow(coef_mat) - 1L)
  sub  <- coef_mat[1L + seq_len(p), , drop = FALSE]

  # Kolom 3 = 2.5%, kolom 7 = 97.5% — paksa numerik agar aman
  ci_lo <- as.numeric(sub[, 3L])
  ci_hi <- as.numeric(sub[, 7L])

  cross <- which(ci_lo < 0 & ci_hi > 0)
  if (length(cross) == 0L) return(character(0))

  var_names[cross]   # mapping posisional kembali ke nama variabel asli
}

# ── Fungsi utama HB Beta dengan seleksi CI otomatis ─────────────────────────
#
# Algoritma:
#  1. Skala semua vars sekali di awal (per kolom, tahan near-constant).
#  2. Jalankan saeHB::Beta() — output internalnya di-suppress agar tidak duplikat.
#  3. Cetak tabel koefisien kita sendiri (bersih, per iterasi).
#  4. Cek 2.5% & 97.5%: jika ada yang menyeberang nol → hapus SD terbesar.
#  5. Ulangi (2–4) sampai semua CI aman atau tidak ada vars tersisa.
#  6. Pakai hasil run terakhir yang valid untuk output domain.
#
# Catatan penting:
#  - saeHB::Beta() di-wrap capture.output() agar tidak mencetak tabel sendiri
#    (yang menyebabkan output duplikat). Kita cetak sendiri versi yang lebih bersih.
#  - Semua status memakai cat() bukan message() supaya tampil di stdout
#    (ikut copy-paste, tidak tenggelam di stderr/merah RStudio).
#
run_hbbeta <- function(df, vars,
                       has_jags    = get("has_jags", envir = .GlobalEnv),
                       iter.mcmc   = 50000,
                       burn.in     = 2000,
                       thin        = 10,
                       iter.update = 5) {

  # ── Fallback jika JAGS tidak tersedia ──────────────────────────────────────
  if (!has_jags) {
    df$y_hb_pct <- df$Estimasi
    df$sd_hb    <- NA_real_
    df$RSE_hb   <- df$RSE
    return(df)
  }

  # ── Persiapan data ──────────────────────────────────────────────────────────
  eps       <- 1e-4
  df$y_prop <- pmax(pmin(df$Estimasi / 100, 1 - eps), eps)

  # Scale per kolom — lebih aman dari scale() bulk untuk kolom near-constant
  for (v in vars) {
    mu_v <- mean(df[[v]], na.rm = TRUE)
    sd_v <- sd(df[[v]],   na.rm = TRUE)
    df[[v]] <- if (is.na(sd_v) || sd_v == 0) 0 else (df[[v]] - mu_v) / sd_v
    df[[v]][is.na(df[[v]]) | is.nan(df[[v]])] <- 0
  }

  # ── Loop seleksi berbasis CI ────────────────────────────────────────────────
  current_vars <- vars
  hb_last      <- NULL
  max_iter     <- length(vars) + 1L   # batas aman: maks 1 drop per iterasi

  for (iter in seq_len(max_iter)) {

    # Tidak ada vars → fallback langsung
    if (length(current_vars) == 0L) {
      cat("  [HB] Tidak ada variabel tersisa setelah seleksi CI → fallback direct\n")
      df$y_hb_pct <- df$Estimasi
      df$sd_hb    <- NA_real_
      df$RSE_hb   <- df$RSE
      return(df)
    }

    form_hb <- as.formula(paste("y_prop ~", paste(current_vars, collapse = "+")))

    # ── saeHB::Beta() di-wrap capture.output() agar tidak duplikat ─────────
    # Output internal saeHB (cat/print) ditangkap & dibuang;
    # kita cetak sendiri tabel koefisien di bawah.
    capture.output(
      hb_cur <- try(
        saeHB::Beta(
          formula     = form_hb,
          data        = df,
          iter.update = iter.update,
          iter.mcmc   = iter.mcmc,
          thin        = thin,
          burn.in     = burn.in
        ),
        silent = TRUE
      )
    )

    # JAGS crash / konvergensi gagal → fallback
    if (inherits(hb_cur, "try-error") || is.null(hb_cur$Est)) {
      cat("  [HB] Iterasi", iter, "gagal (JAGS error) → fallback direct\n")
      df$y_hb_pct <- df$Estimasi
      df$sd_hb    <- NA_real_
      df$RSE_hb   <- df$RSE
      return(df)
    }

    hb_last  <- hb_cur
    coef_mat <- hb_cur$coefficient

    # saeHB tidak mengembalikan tabel koefisien → lewati cek CI
    if (is.null(coef_mat)) {
      cat("  [HB] Tabel koefisien tidak tersedia → lewati cek CI\n")
      break
    }

    # ── Cetak tabel koefisien kita sendiri (bersih, berlabel iterasi) ───────
    cat(sprintf("\n  [HB Iter %d] Vars aktif (%d): %s\n",
                iter, length(current_vars),
                paste(current_vars, collapse = ", ")))
    print(round(coef_mat, 5))
    cat("\n")

    # ── Cek CI ───────────────────────────────────────────────────────────────
    cross_vars <- .vars_cross_ci(coef_mat, current_vars)

    if (length(cross_vars) == 0L) {
      # ✓ Semua aman
      cat(sprintf("  [HB] ✓ Semua CI aman — selesai pada iterasi %d.\n", iter))
      break
    }

    # ✗ Ada yang menyeberang nol → keluarkan yang SD-nya terbesar
    # SD lookup juga posisional (kolom 2 = SD, baris = 1 + posisi di current_vars)
    pos_cross <- match(cross_vars, current_vars)          # posisi di current_vars
    sd_cross  <- setNames(
      as.numeric(coef_mat[1L + pos_cross, 2L]),
      cross_vars
    )
    to_drop   <- names(which.max(sd_cross))

    cat(sprintf("  [HB] ✗ Iter %d — CI menyeberang nol: [%s]\n",
                iter, paste(cross_vars, collapse = ", ")))
    cat(sprintf("  [HB]   → Keluarkan '%s'  (SD = %.5f, terbesar)\n",
                to_drop, sd_cross[[to_drop]]))

    current_vars <- setdiff(current_vars, to_drop)
    # lanjut ke iterasi berikutnya dengan current_vars − 1
  }

  # ── Tulis hasil ke df ───────────────────────────────────────────────────────
  if (is.null(hb_last) || is.null(hb_last$Est)) {
    cat("  [HB] Tidak ada run valid → fallback direct\n")
    df$y_hb_pct <- df$Estimasi
    df$sd_hb    <- NA_real_
    df$RSE_hb   <- df$RSE
  } else {
    cat("  [HB] ✓ Variabel final:", paste(current_vars, collapse = ", "), "\n")
    df$y_hb_pct <- hb_last$Est$MEAN * 100
    df$sd_hb    <- hb_last$Est$SD
    df$RSE_hb   <- df$sd_hb / pmax(df$y_hb_pct / 100, 1e-6) * 100
  }

  df
}

# ── Plot RSE line (2 metode) ─────────────────────────────────────────────
plot_rse <- function(data, col1, col2, lab1, lab2, title) {
  miss <- setdiff(c(col1, col2), names(data))
  if (length(miss) > 0) { message("Kolom tidak ada: ", paste(miss)); return(invisible(NULL)) }
  df_plt <- data %>%
    dplyr::select(Kako, all_of(c(col1, col2))) %>%
    pivot_longer(-Kako, names_to = "Metode", values_to = "RSE") %>%
    mutate(Metode = case_when(Metode == col1 ~ lab1, Metode == col2 ~ lab2, TRUE ~ Metode),
           Kako = as.character(Kako)) %>%
    filter(!is.na(RSE))
  ggplot(df_plt, aes(x = Kako, y = RSE, group = Metode, color = Metode)) +
    geom_line(linewidth = 0.8, alpha = 0.85) +
    geom_point(size = 1.6) +
    geom_hline(yintercept = 25, linetype = "dashed", color = "#c0392b", linewidth = 0.9) +
    annotate("text", x = 1, y = 26.5, label = "Batas 25%", color = "#c0392b",
             hjust = 0, size = 3) +
    scale_color_manual(values = c("#2980b9", "#27ae60")) +
    labs(title = title, x = "Kabupaten/Kota", y = "RSE (%)", color = NULL) +
    theme_minimal(base_size = 11) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6),
          legend.position = "top", plot.title = element_text(face = "bold", size = 12))
}

# ── Plot HB Beta Posterior (dot + CI per domain) ─────────────────────────
plot_hb_posterior <- function(df, title, group_col = NULL) {
  eps   <- 1e-4
  df_p  <- df %>%
    filter(!is.na(y_hb_pct) & !is.na(sd_hb)) %>%
    mutate(
      ci_lo  = pmax((y_hb_pct / 100 - 1.96 * sd_hb) * 100, 0),
      ci_hi  = pmin((y_hb_pct / 100 + 1.96 * sd_hb) * 100, 100),
      Domain = reorder(as.character(Kako), y_hb_pct)
    )
  if (!is.null(group_col) && group_col %in% names(df_p)) {
    df_p$Grup <- df_p[[group_col]]
  } else {
    df_p$Grup <- "Semua Domain"
  }

  p <- ggplot(df_p, aes(x = Domain)) +
    geom_errorbar(aes(ymin = ci_lo, ymax = ci_hi, color = Grup),
                  width = 0.3, alpha = 0.55) +
    geom_point(aes(y = y_hb_pct, color = Grup), size = 2) +
    geom_point(aes(y = Estimasi), color = "#c0392b", size = 1.2, alpha = 0.6, shape = 4) +
    geom_hline(yintercept = mean(df_p$y_hb_pct, na.rm = TRUE),
               linetype = "dashed", color = "#2ecc71", linewidth = 0.8) +
    scale_color_manual(values = c("#2980b9","#e67e22","#8e44ad","#16a085")) +
    labs(title    = title,
         subtitle = "Titik biru = HB | Error bar = 95% CI | Silang merah = Direct | Garis hijau = Rata-rata HB",
         x = "Kabupaten/Kota (urut HB)", y = "Estimasi (%)", color = "Grup") +
    theme_minimal(base_size = 10) +
    theme(axis.text.x  = element_text(angle = 90, vjust = 0.5, size = 5),
          legend.position = "top",
          plot.title      = element_text(face = "bold", size = 12),
          plot.subtitle   = element_text(size = 9))
  p
}

# ── Utilitas evaluasi ────────────────────────────────────────────────────
pct_ok     <- function(x)     round(mean(x < 25, na.rm = TRUE) * 100, 1)
pct_ok_15  <- function(x)     round(mean(x < 15, na.rm = TRUE) * 100, 1)
cv_mean    <- function(x)     round(mean(x, na.rm = TRUE), 2)

tabel_rse <- function(data, var_rse, label) {
  data.frame(
    Skenario      = label,
    `RSE >= 25`   = sum(data[[var_rse]] >= 25, na.rm = TRUE),
    `RSE < 25`    = sum(data[[var_rse]] < 25,  na.rm = TRUE),
    `% < 25%`     = pct_ok(data[[var_rse]]),
    `% < 15%`     = pct_ok_15(data[[var_rse]]),
    `CV Mean`     = cv_mean(data[[var_rse]]),
    check.names   = FALSE
  )
}

metrik <- function(est_model, est_direct, label) {
  data.frame(
    Skenario  = label,
    `RB (%)`  = round(mean((est_model - est_direct) / est_direct, na.rm=TRUE)*100, 3),
    `RMSE`    = round(sqrt(mean((est_model - est_direct)^2, na.rm=TRUE)), 4),
    check.names = FALSE
  )
}

# ── Fungsi utama: jalankan model pada list segmen ────────────────────────
# aux_fn  : seleksi_backward atau seleksi_topn
# model_fn: run_eblup / run_glmm / run_eb_beta / run_hbbeta
# y_col   : "Estimasi" (EBLUP) atau "y_logit" (GLMM/EB/HB)
# Catatan: excl_cols tidak perlu dioper — kedua fungsi seleksi
#          sudah auto-detect kolom non-X secara dinamis
run_on_segments <- function(seg_list, aux_fn, model_fn, y_col,
                             min_n = 4, ...) {
  hasil <- list()
  for (nm in names(seg_list)) {
    df_s <- seg_list[[nm]]
    if (nrow(df_s) < min_n) { message("  Skip ", nm, " (n=", nrow(df_s), ")"); next }
    vars <- tryCatch(
      aux_fn(df_s, y_col = y_col),
      error = function(e) { message("  VarSel error [", nm, "]: ", e$message); character(0) }
    )
    if (length(vars) == 0) {
      message("  Tidak ada var terpilih untuk ", nm, " — skip")
      next
    }
    cat("  [", nm, "]", y_col, "| n =", nrow(df_s),
        "| Vars:", paste(vars, collapse=", "), "\n")
    res <- tryCatch(
      model_fn(df_s, vars, ...),
      error = function(e) {
        message("  Model error [", nm, "]: ", e$message)
        NULL
      }
    )
    if (!is.null(res)) hasil[[nm]] <- res
  }
  if (length(hasil) == 0) return(NULL)
  bind_rows(hasil)
}

# RSE_direct_col: hanya untuk pembagian partisi, BUKAN kandidat aux
# Karena nama tidak berawalan "X", otomatis tereksklusi oleh kedua fungsi seleksi
df_base$RSE_direct_col <- sqrt(df_base$vardir) / pmax(df_base$Estimasi, 1e-6) * 100

brk_jenks <- classIntervals(df_base$RSE_direct_col, n = 2, style = "jenks")
df_base$grup_jenks <- cut(df_base$RSE_direct_col,
                          breaks = brk_jenks$brks,
                          labels = c("G1_RSE_Rendah", "G2_RSE_Tinggi"),
                          include.lowest = TRUE)

cat("Cut-off Jenks:", round(brk_jenks$brks[2], 2), "%\n")

table(df_base$grup_jenks)

med_rse <- median(df_base$RSE_direct_col, na.rm = TRUE)
df_base$grup_equal <- ifelse(df_base$RSE_direct_col <= med_rse,
                              "G1_RSE_Bawah", "G2_RSE_Atas")
cat("Median RSE:", round(med_rse, 2), "%\n")

table(df_base$grup_equal)

# Fungsi helper cluster
make_cluster <- function(df_full, vars_cl, seed = 42) {
  mat <- scale(df_full[, vars_cl, drop = FALSE])
  mat[is.nan(mat)] <- 0
  set.seed(seed)
  km  <- kmeans(mat, centers = 2, nstart = 50, iter.max = 200)
  paste0("Cluster_", km$cluster)
}

# === Global var selection (ALL data) untuk menentukan vars clustering ===

# A1 (Backward) — Normal link
vars_gbk_normal <- seleksi_backward(df_base, y_col = "Estimasi")
cat("Cluster-A1-Normal vars:", paste(vars_gbk_normal, collapse=", "), "\n")

df_base$cluster_bk_normal <- make_cluster(df_base, vars_gbk_normal)

# A2 (TopN) — Normal link
vars_gtn_normal <- seleksi_topn(df_base, y_col = "Estimasi")
cat("Cluster-A2-Normal vars:", paste(vars_gtn_normal, collapse=", "), "\n")

df_base$cluster_tn_normal <- make_cluster(df_base, vars_gtn_normal)

# A1 (Backward) — Logit link
vars_gbk_logit <- seleksi_backward(df_base, y_col = "y_logit")
cat("Cluster-A1-Logit vars:", paste(vars_gbk_logit, collapse=", "), "\n")

df_base$cluster_bk_logit <- make_cluster(df_base, vars_gbk_logit)

# A2 (TopN) — Logit link
vars_gtn_logit <- seleksi_topn(df_base, y_col = "y_logit")
cat("Cluster-A2-Logit vars:", paste(vars_gtn_logit, collapse=", "), "\n")

df_base$cluster_tn_logit <- make_cluster(df_base, vars_gtn_logit)

# Ringkasan cluster
kable(
  data.frame(
    Cluster_Tipe    = c("A1-Normal","A2-Normal","A1-Logit","A2-Logit"),
    Variabel_Global = c(paste(vars_gbk_normal, collapse=", "),
                        paste(vars_gtn_normal, collapse=", "),
                        paste(vars_gbk_logit,  collapse=", "),
                        paste(vars_gtn_logit,  collapse=", "))
  ), caption = "Variabel Global untuk Pembentukan Cluster") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"))

p1 <- ggplot(df_base, aes(x = RSE_direct_col, fill = grup_jenks)) +
  geom_histogram(bins = 25, color = "white", alpha = 0.85) +
  geom_vline(xintercept = brk_jenks$brks[2], linetype="dashed", color="#c0392b") +
  scale_fill_manual(values = c("#2980b9","#e67e22")) +
  labs(title = "C2: RSE Natural Break (Jenks)", x = "RSE Direct (%)", y = "n", fill=NULL) +
  theme_minimal(base_size = 11) + theme(legend.position = "top")

p2 <- ggplot(df_base, aes(x = RSE_direct_col, fill = grup_equal)) +
  geom_histogram(bins = 25, color = "white", alpha = 0.85) +
  geom_vline(xintercept = med_rse, linetype="dashed", color="#c0392b") +
  scale_fill_manual(values = c("#2980b9","#e67e22")) +
  labs(title = "C3: RSE Equal Size (Median)", x = "RSE Direct (%)", y = "n", fill=NULL) +
  theme_minimal(base_size = 11) + theme(legend.position = "top")

p3 <- ggplot(df_base, aes(x = Estimasi, fill = cluster_bk_normal)) +
  geom_boxplot(alpha = 0.8) +
  scale_fill_manual(values = c("#2980b9","#27ae60")) +
  labs(title = "C4 Cluster (A1-Normal)", x = "Estimasi (%)", fill=NULL) +
  theme_minimal(base_size = 11) + theme(legend.position = "top")

p4 <- ggplot(df_base, aes(x = Estimasi, fill = cluster_bk_logit)) +
  geom_boxplot(alpha = 0.8) +
  scale_fill_manual(values = c("#2980b9","#27ae60")) +
  labs(title = "C4 Cluster (A1-Logit)", x = "Estimasi (%)", fill=NULL) +
  theme_minimal(base_size = 11) + theme(legend.position = "top")

(p1 | p2) / (p3 | p4)

cat("=== Skenario 1: EBLUP Backward All ===\n")

df_s1 <- run_on_segments(
  seg_list  = list("All" = df_base),
  aux_fn    = seleksi_backward,
  model_fn  = run_eblup,
  y_col     = "Estimasi"
)

summary(df_s1[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s1, "RSE_direct","RSE_eblup","Direct","EBLUP-BK All","S1: EBLUP Backward — All Domain")

cat("=== Skenario 2: EBLUP Backward RSE-NB ===\n")

df_s2 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_backward, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s2[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s2,"RSE_direct","RSE_eblup","Direct","EBLUP-BK RSE-NB","S2: EBLUP Backward — RSE Natural Break")

cat("=== Skenario 3: EBLUP Backward RSE-ES ===\n")

df_s3 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_backward, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s3[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s3,"RSE_direct","RSE_eblup","Direct","EBLUP-BK RSE-ES","S3: EBLUP Backward — RSE Equal Size")

cat("=== Skenario 4: EBLUP Backward Cluster ===\n")

df_s4 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_bk_normal),
  aux_fn   = seleksi_backward, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s4[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s4,"RSE_direct","RSE_eblup","Direct","EBLUP-BK Cluster","S4: EBLUP Backward — Cluster Aux")

cat("=== Skenario 5: EBLUP TopN All ===\n")

df_s5 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_topn, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s5[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s5,"RSE_direct","RSE_eblup","Direct","EBLUP-TN All","S5: EBLUP Top-n/10 — All Domain")

cat("=== Skenario 6: EBLUP TopN RSE-NB ===\n")

df_s6 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_topn, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s6[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s6,"RSE_direct","RSE_eblup","Direct","EBLUP-TN RSE-NB","S6: EBLUP Top-n/10 — RSE Natural Break")

cat("=== Skenario 7: EBLUP TopN RSE-ES ===\n")

df_s7 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_topn, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s7[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s7,"RSE_direct","RSE_eblup","Direct","EBLUP-TN RSE-ES","S7: EBLUP Top-n/10 — RSE Equal Size")

cat("=== Skenario 8: EBLUP TopN Cluster ===\n")

df_s8 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_tn_normal),
  aux_fn   = seleksi_topn, model_fn = run_eblup, y_col = "Estimasi"
)

summary(df_s8[, c("Estimasi","y_eblup","RSE_direct","RSE_eblup")])

plot_rse(df_s8,"RSE_direct","RSE_eblup","Direct","EBLUP-TN Cluster","S8: EBLUP Top-n/10 — Cluster Aux")

cat("=== Skenario 9: GLMM Backward All ===\n")

df_s9 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_backward, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s9[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s9,"RSE_direct","RSE_glmm","Direct","GLMM-BK All","S9: GLMM Backward — All Domain")

cat("=== Skenario 10: GLMM Backward RSE-NB ===\n")

df_s10 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_backward, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s10[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s10,"RSE_direct","RSE_glmm","Direct","GLMM-BK RSE-NB","S10: GLMM Backward — RSE Natural Break")

cat("=== Skenario 11: GLMM Backward RSE-ES ===\n")

df_s11 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_backward, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s11[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s11,"RSE_direct","RSE_glmm","Direct","GLMM-BK RSE-ES","S11: GLMM Backward — RSE Equal Size")

cat("=== Skenario 12: GLMM Backward Cluster ===\n")

df_s12 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_bk_logit),
  aux_fn   = seleksi_backward, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s12[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s12,"RSE_direct","RSE_glmm","Direct","GLMM-BK Cluster","S12: GLMM Backward — Cluster Aux")

cat("=== Skenario 13: GLMM TopN All ===\n")

df_s13 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_topn, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s13[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s13,"RSE_direct","RSE_glmm","Direct","GLMM-TN All","S13: GLMM Top-n/10 — All Domain")

cat("=== Skenario 14: GLMM TopN RSE-NB ===\n")

df_s14 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_topn, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s14[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s14,"RSE_direct","RSE_glmm","Direct","GLMM-TN RSE-NB","S14: GLMM Top-n/10 — RSE Natural Break")

cat("=== Skenario 15: GLMM TopN RSE-ES ===\n")

df_s15 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_topn, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s15[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s15,"RSE_direct","RSE_glmm","Direct","GLMM-TN RSE-ES","S15: GLMM Top-n/10 — RSE Equal Size")

cat("=== Skenario 16: GLMM TopN Cluster ===\n")

df_s16 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_tn_logit),
  aux_fn   = seleksi_topn, model_fn = run_glmm, y_col = "y_logit"
)

summary(df_s16[, c("Estimasi","y_glmm","RSE_direct","RSE_glmm")])

plot_rse(df_s16,"RSE_direct","RSE_glmm","Direct","GLMM-TN Cluster","S16: GLMM Top-n/10 — Cluster Aux")

cat("=== Skenario 17: EB Backward All ===\n")

df_s17 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_backward, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s17[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s17,"RSE_direct","RSE_eb","Direct","EB-BK All","S17: EB Beta-Binomial Backward — All")

cat("=== Skenario 18: EB Backward RSE-NB ===\n")

df_s18 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_backward, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s18[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s18,"RSE_direct","RSE_eb","Direct","EB-BK RSE-NB","S18: EB Beta-Binomial Backward — RSE NB")

cat("=== Skenario 19: EB Backward RSE-ES ===\n")

df_s19 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_backward, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s19[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s19,"RSE_direct","RSE_eb","Direct","EB-BK RSE-ES","S19: EB Beta-Binomial Backward — RSE ES")

cat("=== Skenario 20: EB Backward Cluster ===\n")

df_s20 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_bk_logit),
  aux_fn   = seleksi_backward, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s20[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s20,"RSE_direct","RSE_eb","Direct","EB-BK Cluster","S20: EB Beta-Binomial Backward — Cluster")

cat("=== Skenario 21: EB TopN All ===\n")

df_s21 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_topn, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s21[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s21,"RSE_direct","RSE_eb","Direct","EB-TN All","S21: EB Beta-Binomial Top-n/10 — All")

cat("=== Skenario 22: EB TopN RSE-NB ===\n")

df_s22 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_topn, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s22[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s22,"RSE_direct","RSE_eb","Direct","EB-TN RSE-NB","S22: EB Beta-Binomial Top-n/10 — RSE NB")

cat("=== Skenario 23: EB TopN RSE-ES ===\n")

df_s23 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_topn, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s23[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s23,"RSE_direct","RSE_eb","Direct","EB-TN RSE-ES","S23: EB Beta-Binomial Top-n/10 — RSE ES")

cat("=== Skenario 24: EB TopN Cluster ===\n")

df_s24 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_tn_logit),
  aux_fn   = seleksi_topn, model_fn = run_eb_beta, y_col = "y_logit"
)

summary(df_s24[, c("Estimasi","y_eb_pct","RSE_direct","RSE_eb")])

plot_rse(df_s24,"RSE_direct","RSE_eb","Direct","EB-TN Cluster","S24: EB Beta-Binomial Top-n/10 — Cluster")

cat("=== Skenario 25: HB Backward All ===\n")

df_s25 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_backward, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s25[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s25,"RSE","RSE_hb","Direct","HB-BK All","S25: HB Beta Backward — All Domain")

plot_hb_posterior(df_s25, "S25: HB Beta Backward — All Domain")

cat("=== Skenario 26: HB Backward RSE-NB ===\n")

df_s26 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_backward, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s26[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s26,"RSE","RSE_hb","Direct","HB-BK RSE-NB","S26: HB Beta Backward — RSE Natural Break")

plot_hb_posterior(df_s26, "S26: HB Beta Backward — RSE Natural Break", group_col = "grup_jenks")

cat("=== Skenario 27: HB Backward RSE-ES ===\n")

df_s27 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_backward, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s27[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s27,"RSE","RSE_hb","Direct","HB-BK RSE-ES","S27: HB Beta Backward — RSE Equal Size")

plot_hb_posterior(df_s27, "S27: HB Beta Backward — RSE Equal Size", group_col = "grup_equal")

cat("=== Skenario 28: HB Backward Cluster ===\n")

df_s28 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_bk_logit),
  aux_fn   = seleksi_backward, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s28[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s28,"RSE","RSE_hb","Direct","HB-BK Cluster","S28: HB Beta Backward — Cluster Aux")

plot_hb_posterior(df_s28, "S28: HB Beta Backward — Cluster Aux",
                  group_col = "cluster_bk_logit")

cat("=== Skenario 29: HB TopN All ===\n")

df_s29 <- run_on_segments(
  seg_list = list("All" = df_base),
  aux_fn   = seleksi_topn, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s29[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s29,"RSE","RSE_hb","Direct","HB-TN All","S29: HB Beta Top-n/10 — All Domain")

plot_hb_posterior(df_s29, "S29: HB Beta Top-n/10 — All Domain")

cat("=== Skenario 30: HB TopN RSE-NB ===\n")

df_s30 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_jenks),
  aux_fn   = seleksi_topn, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s30[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s30,"RSE","RSE_hb","Direct","HB-TN RSE-NB","S30: HB Beta Top-n/10 — RSE Natural Break")

plot_hb_posterior(df_s30, "S30: HB Beta Top-n/10 — RSE Natural Break", group_col = "grup_jenks")

cat("=== Skenario 31: HB TopN RSE-ES ===\n")

df_s31 <- run_on_segments(
  seg_list = split(df_base, df_base$grup_equal),
  aux_fn   = seleksi_topn, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s31[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s31,"RSE","RSE_hb","Direct","HB-TN RSE-ES","S31: HB Beta Top-n/10 — RSE Equal Size")

plot_hb_posterior(df_s31, "S31: HB Beta Top-n/10 — RSE Equal Size", group_col = "grup_equal")

cat("=== Skenario 32: HB TopN Cluster ===\n")

df_s32 <- run_on_segments(
  seg_list = split(df_base, df_base$cluster_tn_logit),
  aux_fn   = seleksi_topn, model_fn = run_hbbeta, y_col = "y_logit", min_n = 5
)

summary(df_s32[, c("Estimasi","y_hb_pct","RSE","RSE_hb")])

plot_rse(df_s32,"RSE","RSE_hb","Direct","HB-TN Cluster","S32: HB Beta Top-n/10 — Cluster Aux")

plot_hb_posterior(df_s32, "S32: HB Beta Top-n/10 — Cluster Aux",
                  group_col = "cluster_tn_logit")