# 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)
  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[vars] <- as.data.frame(scale(df[vars]))
  df[vars][is.nan(as.matrix(df[vars]))] <- 0

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

  ctrl <- glmmTMBControl(optCtrl = list(iter.max = 400, eval.max = 600))
  m    <- try(glmmTMB(form, data = df, 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 {
    df$y_glmm <- plogis(predict(m, type = "link", re.form = NULL)) * 100

    # Jackknife — seluruh loop dibungkus tryCatch
    # Jika salah satu iterasi fatal crash, fallback ke MSE empiris
    mse_v <- tryCatch({
      n_d      <- nrow(df)
      jack_mat <- matrix(NA_real_, nrow = n_d, ncol = n_d)

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

      vapply(seq_len(n_d), function(j) {
        ev <- jack_mat[, j][!is.na(jack_mat[, j])]
        if (length(ev) >= 2) {
          tb <- mean(ev)
          ((length(ev) - 1) / length(ev)) * mean((ev - tb)^2)
        } else {
          (df$RSE[j] / 100 * df$y_glmm[j])^2   # fallback per-domain
        }
      }, numeric(1))

    }, error = function(e) {
      # Jackknife crash total → fallback ke MSE empiris semua domain
      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
}

run_hbbeta <- function(df, vars,
                       has_jags = get("has_jags", envir = .GlobalEnv)) {
  if (!has_jags) {
    df$y_hb_pct <- df$Estimasi; df$sd_hb <- NA_real_; df$RSE_hb <- df$RSE
    return(df)
  }
  eps       <- 1e-4
  df$y_prop <- pmax(pmin(df$Estimasi / 100, 1 - eps), eps)
  df[vars]  <- as.data.frame(scale(df[vars]))
  df[vars][is.nan(as.matrix(df[vars]))] <- 0
  form_hb   <- as.formula(paste("y_prop ~", paste(vars, collapse = "+")))
  hb <- try(saeHB::Beta(
    formula = form_hb, data = df,
    iter.update = 15, iter.mcmc = 4000, thin = 10, burn.in = 1000
  ), silent = TRUE)
  if (inherits(hb, "try-error") || is.null(hb$Est)) {
    message("  HB Beta gagal → fallback direct")
    df$y_hb_pct <- df$Estimasi; df$sd_hb <- NA_real_; df$RSE_hb <- df$RSE
  } else {
    df$y_hb_pct <- hb$Est$MEAN * 100
    df$sd_hb    <- hb$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")