Ringkasan analisis: 32 skenario SAE — 2 metode seleksi aux (Backward vs Top-n/10) × 4 model (EBLUP, GLMM-Logit, EB Beta-Binomial, HB Beta) × 4 partisi (All, RSE Natural Break, RSE Equal Size, Cluster-Aux). Baseline: Direct Estimator.
1. Library & Cek Paket
Library
# 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::recodeCek JAGS & saeHB
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")## Status JAGS/saeHB: Tersedia2. Load Data & Persiapan
Load Data
Estimasi_531 <- read_excel("C:/Users/User/Downloads/Hasil_Sulawesi_Kab2025.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")## Jumlah domain: 81## Jumlah variabel: 118Variabel Rasio
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
)Pre-filter Kandidat Global
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")## Variabel kandidat ( 125 ):
## X1, X2, X8, X9, X13, X14, X15, X16, X17, X18, X19, X21, X22, X23, X24, X25, X26, X32, X37, X38, X39, X40, X41, X43, X44, X45, X1_jumlah, X1_mean, X2_jumlah, X2_mean, X8_jumlah, X8_mean, X9_jumlah, X9_mean, X13_jumlah, X13_mean, X14_jumlah, X14_mean, X15_jumlah, X15_mean, X16_jumlah, X16_mean, X17_jumlah, X17_mean, X18_jumlah, X18_mean, X19_jumlah, X19_mean, X21_jumlah, X21_mean, X22_jumlah, X22_mean, X23_jumlah, X23_mean, X24_jumlah, X24_mean, X25_jumlah, X25_mean, X26_jumlah, X26_mean, X32_jumlah, X32_mean, X34_jumlah, X34_mean, X37_jumlah, X37_mean, X38_jumlah, X38_mean, X39_jumlah, X39_mean, X40_jumlah, X40_mean, X41_jumlah, X41_mean, X43_jumlah, X43_mean, X44_jumlah, X44_mean, X45_jumlah, X45_mean, X6_jumlah, X6_mean, X7_jumlah, X7_mean, X10_jumlah, X10_mean, X11_jumlah, X11_mean, X12_jumlah, X12_mean, X27_jumlah, X27_mean, X28_jumlah, X28_mean, X29_jumlah, X29_mean, X30_jumlah, X30_mean, X31_jumlah, X31_mean, X33_jumlah, X33_mean, X35_jumlah, X35_mean, X36_jumlah, X36_mean, X46, X3, X4, X5, X47, X48, X6, X7, X10, X33, X11, X12, X27, X28, X29, X30, X31, X35, X36##
## Total dikecualikan: 9 kolom3. Fungsi Pembantu
Seleksi A1: Backward
# 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 A2: Top-n/10
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)
}Model: EBLUP
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
}Model: GLMM (glmmTMB)
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
}Model: EB Beta-Binomial
#' 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
}Model: HB Beta
# ── 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
}Fungsi Plot & Evaluasi
# ── 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
)
}
# ── Korelasi aux terpilih vs Direct dan Logit(Direct) ────────────────────
# Dipanggil oleh run_on_segments setelah seleksi variabel selesai,
# sebelum model dijalankan — agar transparan hubungan tiap aux dengan respons.
show_var_corr <- function(df, vars, estimasi_col = "Estimasi") {
if (length(vars) == 0 || !estimasi_col %in% names(df)) return(invisible(NULL))
eps <- 1e-4
y_dir <- df[[estimasi_col]]
p <- pmax(pmin(y_dir / 100, 1 - eps), eps)
y_lgit <- log(p / (1 - p))
cat(sprintf(" %-16s %12s %12s\n", "Variabel", "cor(Direct)", "cor(Logit)"))
cat(" ", strrep("-", 42), "\n", sep = "")
for (v in vars) {
x <- df[[v]]
r1 <- tryCatch(cor(x, y_dir, use = "complete.obs"), error = function(e) NA_real_)
r2 <- tryCatch(cor(x, y_lgit, use = "complete.obs"), error = function(e) NA_real_)
cat(sprintf(" %-16s %12.4f %12.4f\n", v, r1, r2))
}
cat("\n")
}
# ── 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")
show_var_corr(df_s, vars)
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)
}4. Pembagian Data
C2 — RSE Natural Break (Jenks)
# 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")## Cut-off Jenks: 58.89 %##
## G1_RSE_Rendah G2_RSE_Tinggi
## 60 21C3 — RSE Equal Size (Median)
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")## Median RSE: 48.37 %##
## G1_RSE_Bawah G2_RSE_Atas
## 41 40C4 — Cluster-Aux (k=2)
# 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")## Cluster-A1-Normal vars: X40, X29_meandf_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")## Cluster-A2-Normal vars: X40, X8, X45, X36_mean, X29_mean, X27, X31_meandf_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")## Cluster-A1-Logit vars: X36_mean, X40df_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")## Cluster-A2-Logit vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45df_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"))| Cluster_Tipe | Variabel_Global |
|---|---|
| A1-Normal | X40, X29_mean |
| A2-Normal | X40, X8, X45, X36_mean, X29_mean, X27, X31_mean |
| A1-Logit | X36_mean, X40 |
| A2-Logit | X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45 |
Distribusi Semua Partisi
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)
5. EBLUP — Normal Link
Skenario 1–8 — Model EBLUP Fay-Herriot pada skala normal. Skenario 1–4: seleksi Backward (A1). Skenario 5–8: seleksi Top-n/10 (A2). Masing-masing 4 partisi: All, RSE-NB, RSE-ES, Cluster.
S1 — Backward · All
## === Skenario 1: EBLUP Backward All ===df_s1 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_backward,
model_fn = run_eblup,
y_col = "Estimasi"
)## [ All ] Estimasi | n = 81 | Vars: X40, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3690 -0.3705
## X29_mean -0.3220 -0.3921## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.572 Min. :26.1 Min. :20.7
## 1st Qu.: 4.10 1st Qu.: 4.531 1st Qu.:39.3 1st Qu.:26.2
## Median : 8.24 Median : 6.379 Median :48.4 Median :33.1
## Mean : 8.65 Mean : 6.384 Mean :52.3 Mean :34.8
## 3rd Qu.:11.12 3rd Qu.: 7.923 3rd Qu.:61.0 3rd Qu.:40.7
## Max. :26.70 Max. :14.102 Max. :98.8 Max. :87.2
S1: RSE Direct vs EBLUP Backward (All)
S2 — Backward · RSE-NB
## === Skenario 2: EBLUP Backward RSE-NB ===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"
)## [ G1_RSE_Rendah ] Estimasi | n = 60 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828## [ G2_RSE_Tinggi ] Estimasi | n = 21 | Vars: X41_jumlah, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.4969 0.5170
## X28_mean -0.4481 -0.6093## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.519 Min. :26.1 Min. : 19.6
## 1st Qu.: 4.10 1st Qu.:3.982 1st Qu.:39.3 1st Qu.: 55.4
## Median : 8.24 Median :8.250 Median :48.4 Median : 55.4
## Mean : 8.65 Mean :6.680 Mean :52.3 Mean : 50.8
## 3rd Qu.:11.12 3rd Qu.:8.250 3rd Qu.:61.0 3rd Qu.: 55.4
## Max. :26.70 Max. :8.250 Max. :98.8 Max. :184.6plot_rse(df_s2,"RSE_direct","RSE_eblup","Direct","EBLUP-BK RSE-NB","S2: EBLUP Backward — RSE Natural Break")
S2: EBLUP Backward RSE Natural Break
S3 — Backward · RSE-ES
## === Skenario 3: EBLUP Backward RSE-ES ===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"
)## [ G1_RSE_Bawah ] Estimasi | n = 41 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4689 0.4949
## X31_jumlah -0.4520 -0.5349## [ G2_RSE_Atas ] Estimasi | n = 40 | Vars: X40_mean, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_mean -0.4352 -0.3638
## X28_mean -0.4137 -0.5322## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.331 Min. :26.1 Min. : 15.1
## 1st Qu.: 4.10 1st Qu.:2.654 1st Qu.:39.3 1st Qu.: 20.6
## Median : 8.24 Median :7.940 Median :48.4 Median : 42.0
## Mean : 8.65 Mean :5.407 Mean :52.3 Mean : 34.1
## 3rd Qu.:11.12 3rd Qu.:7.940 3rd Qu.:61.0 3rd Qu.: 42.0
## Max. :26.70 Max. :7.940 Max. :98.8 Max. :195.8plot_rse(df_s3,"RSE_direct","RSE_eblup","Direct","EBLUP-BK RSE-ES","S3: EBLUP Backward — RSE Equal Size")
S3: EBLUP Backward RSE Equal Size
S4 — Backward · Cluster
## === Skenario 4: EBLUP Backward Cluster ===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"
)## [ Cluster_1 ] Estimasi | n = 60 | Vars: X47, X8
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.3911 0.3619
## X8 0.3176 0.3266
##
## [ Cluster_2 ] Estimasi | n = 21 | Vars: X35, X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X35 -0.3773 -0.2798
## X40_jumlah -0.3419 -0.3562## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.636 Min. :26.1 Min. :20.6
## 1st Qu.: 4.10 1st Qu.: 3.973 1st Qu.:39.3 1st Qu.:28.3
## Median : 8.24 Median : 6.432 Median :48.4 Median :34.4
## Mean : 8.65 Mean : 6.525 Mean :52.3 Mean :36.2
## 3rd Qu.:11.12 3rd Qu.: 8.427 3rd Qu.:61.0 3rd Qu.:42.5
## Max. :26.70 Max. :12.879 Max. :98.8 Max. :78.3plot_rse(df_s4,"RSE_direct","RSE_eblup","Direct","EBLUP-BK Cluster","S4: EBLUP Backward — Cluster Aux")
S4: EBLUP Backward Cluster-Aux
S5 — Top-n/10 · All
## === Skenario 5: EBLUP TopN All ===df_s5 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_topn, model_fn = run_eblup, y_col = "Estimasi"
)## [ All ] Estimasi | n = 81 | Vars: X40, X8, X45, X36_mean, X29_mean, X27, X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3690 -0.3705
## X8 0.3439 0.3797
## X45 -0.3224 -0.3604
## X36_mean -0.3220 -0.4311
## X29_mean -0.3220 -0.3921
## X27 0.3197 0.3124
## X31_mean -0.3188 -0.3953## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.545 Min. :26.1 Min. :22.6
## 1st Qu.: 4.10 1st Qu.: 4.469 1st Qu.:39.3 1st Qu.:28.0
## Median : 8.24 Median : 6.658 Median :48.4 Median :35.7
## Mean : 8.65 Mean : 6.520 Mean :52.3 Mean :37.3
## 3rd Qu.:11.12 3rd Qu.: 8.162 3rd Qu.:61.0 3rd Qu.:43.5
## Max. :26.70 Max. :13.824 Max. :98.8 Max. :92.2
S5: EBLUP Top-n/10 All
S6 — Top-n/10 · RSE-NB
## === Skenario 6: EBLUP TopN RSE-NB ===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"
)## [ G1_RSE_Rendah ] Estimasi | n = 60 | Vars: X3, X4, X27, X31_jumlah, X29_mean, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X4 0.4425 0.4608
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X39_jumlah -0.3610 -0.3605## [ G2_RSE_Tinggi ] Estimasi | n = 21 | Vars: X41_jumlah, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.4969 0.5170
## X28_mean -0.4481 -0.6093## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.519 Min. :26.1 Min. : 19.6
## 1st Qu.: 4.10 1st Qu.:3.982 1st Qu.:39.3 1st Qu.: 55.4
## Median : 8.24 Median :8.250 Median :48.4 Median : 55.4
## Mean : 8.65 Mean :6.680 Mean :52.3 Mean : 50.8
## 3rd Qu.:11.12 3rd Qu.:8.250 3rd Qu.:61.0 3rd Qu.: 55.4
## Max. :26.70 Max. :8.250 Max. :98.8 Max. :184.6plot_rse(df_s6,"RSE_direct","RSE_eblup","Direct","EBLUP-TN RSE-NB","S6: EBLUP Top-n/10 — RSE Natural Break")
S6: EBLUP Top-n/10 RSE Natural Break
S7 — Top-n/10 · RSE-ES
## === Skenario 7: EBLUP TopN RSE-ES ===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"
)## [ G1_RSE_Bawah ] Estimasi | n = 41 | Vars: X4, X3, X27, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X4 0.4782 0.5323
## X3 0.4689 0.4949
## X27 0.4608 0.5500
## X31_jumlah -0.4520 -0.5349## [ G2_RSE_Atas ] Estimasi | n = 40 | Vars: X25_mean, X47, X40_mean, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X25_mean -0.4527 -0.5200
## X47 0.4523 0.4163
## X40_mean -0.4352 -0.3638
## X28_mean -0.4137 -0.5322## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.468 Min. :26.1 Min. : 13.2
## 1st Qu.: 4.10 1st Qu.:2.857 1st Qu.:39.3 1st Qu.: 24.4
## Median : 8.24 Median :7.940 Median :48.4 Median : 42.0
## Mean : 8.65 Mean :5.451 Mean :52.3 Mean : 36.8
## 3rd Qu.:11.12 3rd Qu.:7.940 3rd Qu.:61.0 3rd Qu.: 42.0
## Max. :26.70 Max. :7.940 Max. :98.8 Max. :145.7plot_rse(df_s7,"RSE_direct","RSE_eblup","Direct","EBLUP-TN RSE-ES","S7: EBLUP Top-n/10 — RSE Equal Size")
S7: EBLUP Top-n/10 RSE Equal Size
S8 — Top-n/10 · Cluster
## === Skenario 8: EBLUP TopN Cluster ===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"
)## [ Cluster_1 ] Estimasi | n = 11 | Vars: X47
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.7038 0.4875
##
## [ Cluster_2 ] Estimasi | n = 70 | Vars: X40, X48, X8, X35_jumlah, X45_jumlah, X44_jumlah, X3
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3196 -0.3335
## X48 -0.2726 -0.2762
## X8 0.2662 0.3094
## X35_jumlah -0.2568 -0.2167
## X45_jumlah -0.2369 -0.1636
## X44_jumlah -0.2345 -0.2020
## X3 0.2332 0.1541## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.718 Min. :26.1 Min. :22.3
## 1st Qu.: 4.10 1st Qu.: 4.221 1st Qu.:39.3 1st Qu.:29.1
## Median : 8.24 Median : 6.966 Median :48.4 Median :38.2
## Mean : 8.65 Mean : 6.720 Mean :52.3 Mean :39.0
## 3rd Qu.:11.12 3rd Qu.: 8.955 3rd Qu.:61.0 3rd Qu.:45.8
## Max. :26.70 Max. :14.196 Max. :98.8 Max. :76.2plot_rse(df_s8,"RSE_direct","RSE_eblup","Direct","EBLUP-TN Cluster","S8: EBLUP Top-n/10 — Cluster Aux")
S8: EBLUP Top-n/10 Cluster-Aux
6. GLMM — Logit Link glmmTMB
Skenario 9–16 — GLMM binomial via
glmmTMB, random intercept per domain (1|Kako).
Proporsi dikonversi ke count efektif n_eff = p(1-p)/SE².
MSE via Jackknife leave-one-domain-out. Seleksi variabel pada skala
logit.
S9 — Backward · All
## === Skenario 9: GLMM Backward All ===df_s9 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_backward, model_fn = run_glmm, y_col = "y_logit"
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 1.07 Min. :26.1 Min. : 1.18
## 1st Qu.: 4.10 1st Qu.: 5.15 1st Qu.:39.3 1st Qu.: 9.07
## Median : 8.24 Median : 7.44 Median :48.4 Median : 18.04
## Mean : 8.65 Mean : 7.54 Mean :52.3 Mean : 25.79
## 3rd Qu.:11.12 3rd Qu.: 8.90 3rd Qu.:61.0 3rd Qu.: 33.13
## Max. :26.70 Max. :17.64 Max. :98.8 Max. :196.33
S9: GLMM Backward All
S10 — Backward · RSE-NB
## === Skenario 10: GLMM Backward RSE-NB ===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"
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.567 Min. :26.1 Min. : 2.35
## 1st Qu.: 4.10 1st Qu.: 3.519 1st Qu.:39.3 1st Qu.: 2.91
## Median : 8.24 Median : 7.865 Median :48.4 Median : 5.74
## Mean : 8.65 Mean : 7.608 Mean :52.3 Mean : 10.20
## 3rd Qu.:11.12 3rd Qu.:10.369 3rd Qu.:61.0 3rd Qu.: 12.69
## Max. :26.70 Max. :16.145 Max. :98.8 Max. :122.66plot_rse(df_s10,"RSE_direct","RSE_glmm","Direct","GLMM-BK RSE-NB","S10: GLMM Backward — RSE Natural Break")
S10: GLMM Backward RSE Natural Break
S11 — Backward · RSE-ES
## === Skenario 11: GLMM Backward RSE-ES ===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"
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.417 Min. :26.1 Min. : 0.843
## 1st Qu.: 4.10 1st Qu.: 4.191 1st Qu.:39.3 1st Qu.: 1.036
## Median : 8.24 Median : 6.490 Median :48.4 Median : 3.031
## Mean : 8.65 Mean : 7.338 Mean :52.3 Mean : 8.194
## 3rd Qu.:11.12 3rd Qu.:10.685 3rd Qu.:61.0 3rd Qu.:14.320
## Max. :26.70 Max. :13.162 Max. :98.8 Max. :65.603plot_rse(df_s11,"RSE_direct","RSE_glmm","Direct","GLMM-BK RSE-ES","S11: GLMM Backward — RSE Equal Size")
S11: GLMM Backward RSE Equal Size
S12 — Backward · Cluster
## === Skenario 12: GLMM Backward Cluster ===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"
)## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.944 Min. :26.1 Min. : 1.83
## 1st Qu.: 4.10 1st Qu.: 4.887 1st Qu.:39.3 1st Qu.: 3.62
## Median : 8.24 Median : 7.455 Median :48.4 Median : 16.36
## Mean : 8.65 Mean : 7.623 Mean :52.3 Mean : 24.40
## 3rd Qu.:11.12 3rd Qu.:10.790 3rd Qu.:61.0 3rd Qu.: 30.45
## Max. :26.70 Max. :14.525 Max. :98.8 Max. :220.65plot_rse(df_s12,"RSE_direct","RSE_glmm","Direct","GLMM-BK Cluster","S12: GLMM Backward — Cluster Aux")
S12: GLMM Backward Cluster-Aux
S13 — Top-n/10 · All
## === Skenario 13: GLMM TopN All ===df_s13 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_topn, model_fn = run_glmm, y_col = "y_logit"
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 1.08 Min. :26.1 Min. : 1.77
## 1st Qu.: 4.10 1st Qu.: 5.09 1st Qu.:39.3 1st Qu.: 9.33
## Median : 8.24 Median : 7.50 Median :48.4 Median : 18.37
## Mean : 8.65 Mean : 7.55 Mean :52.3 Mean : 25.10
## 3rd Qu.:11.12 3rd Qu.: 9.36 3rd Qu.:61.0 3rd Qu.: 33.50
## Max. :26.70 Max. :18.16 Max. :98.8 Max. :161.04
S13: GLMM Top-n/10 All
S14 — Top-n/10 · RSE-NB
## === Skenario 14: GLMM TopN RSE-NB ===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"
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.605 Min. :26.1 Min. : 2.00
## 1st Qu.: 4.10 1st Qu.: 3.382 1st Qu.:39.3 1st Qu.: 2.96
## Median : 8.24 Median : 7.698 Median :48.4 Median : 5.86
## Mean : 8.65 Mean : 7.582 Mean :52.3 Mean : 9.90
## 3rd Qu.:11.12 3rd Qu.:10.504 3rd Qu.:61.0 3rd Qu.: 10.82
## Max. :26.70 Max. :15.383 Max. :98.8 Max. :118.59plot_rse(df_s14,"RSE_direct","RSE_glmm","Direct","GLMM-TN RSE-NB","S14: GLMM Top-n/10 — RSE Natural Break")
S14: GLMM Top-n/10 RSE Natural Break
S15 — Top-n/10 · RSE-ES
## === Skenario 15: GLMM TopN RSE-ES ===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"
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.65 Min. :26.1 Min. : 1.04
## 1st Qu.: 4.10 1st Qu.: 4.51 1st Qu.:39.3 1st Qu.: 1.55
## Median : 8.24 Median : 6.58 Median :48.4 Median : 3.71
## Mean : 8.65 Mean : 7.36 Mean :52.3 Mean : 8.28
## 3rd Qu.:11.12 3rd Qu.:10.90 3rd Qu.:61.0 3rd Qu.:13.87
## Max. :26.70 Max. :14.75 Max. :98.8 Max. :58.20plot_rse(df_s15,"RSE_direct","RSE_glmm","Direct","GLMM-TN RSE-ES","S15: GLMM Top-n/10 — RSE Equal Size")
S15: GLMM Top-n/10 RSE Equal Size
S16 — Top-n/10 · Cluster
## === Skenario 16: GLMM TopN Cluster ===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"
)## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 1.17 Min. :26.1 Min. : 2.05
## 1st Qu.: 4.10 1st Qu.: 5.05 1st Qu.:39.3 1st Qu.: 7.17
## Median : 8.24 Median : 7.57 Median :48.4 Median : 18.12
## Mean : 8.65 Mean : 7.67 Mean :52.3 Mean : 23.18
## 3rd Qu.:11.12 3rd Qu.: 9.31 3rd Qu.:61.0 3rd Qu.: 30.26
## Max. :26.70 Max. :19.18 Max. :98.8 Max. :166.25plot_rse(df_s16,"RSE_direct","RSE_glmm","Direct","GLMM-TN Cluster","S16: GLMM Top-n/10 — Cluster Aux")
S16: GLMM Top-n/10 Cluster-Aux
7. EB Beta-Binomial Empirical Bayes
Skenario 17–24 — EB Beta-Binomial. Regresi
logit(p) ~ vars menghasilkan μᵢ (prior mean
domain-spesifik). Presisi φ diestimasi via Method of
Moments. Estimasi EB: p̂ᵢ = (aᵢ + yᵢ) / (aᵢ + bᵢ + nᵢ). MSE
dari variansi posterior Beta analytik.
S17 — Backward · All
## === Skenario 17: EB Backward All ===df_s17 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_backward, model_fn = run_eb_beta, y_col = "y_logit"
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.585 Min. :26.1 Min. :24.6
## 1st Qu.: 4.10 1st Qu.: 4.513 1st Qu.:39.3 1st Qu.:34.5
## Median : 8.24 Median : 7.550 Median :48.4 Median :41.8
## Mean : 8.65 Mean : 7.824 Mean :52.3 Mean :41.8
## 3rd Qu.:11.12 3rd Qu.: 9.579 3rd Qu.:61.0 3rd Qu.:46.3
## Max. :26.70 Max. :18.326 Max. :98.8 Max. :86.9
S17: EB Beta Backward All
S18 — Backward · RSE-NB
## === Skenario 18: EB Backward RSE-NB ===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"
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.497 Min. :26.1 Min. :24.0
## 1st Qu.: 4.10 1st Qu.: 4.043 1st Qu.:39.3 1st Qu.:32.2
## Median : 8.24 Median : 8.093 Median :48.4 Median :37.1
## Mean : 8.65 Mean : 8.023 Mean :52.3 Mean :40.1
## 3rd Qu.:11.12 3rd Qu.:10.637 3rd Qu.:61.0 3rd Qu.:46.8
## Max. :26.70 Max. :17.925 Max. :98.8 Max. :91.2plot_rse(df_s18,"RSE_direct","RSE_eb","Direct","EB-BK RSE-NB","S18: EB Beta-Binomial Backward — RSE NB")
S18: EB Beta Backward RSE Natural Break
S19 — Backward · RSE-ES
## === Skenario 19: EB Backward RSE-ES ===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"
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.498 Min. :26.1 Min. :22.1
## 1st Qu.: 4.10 1st Qu.: 4.314 1st Qu.:39.3 1st Qu.:29.8
## Median : 8.24 Median : 8.221 Median :48.4 Median :39.4
## Mean : 8.65 Mean : 7.903 Mean :52.3 Mean :41.8
## 3rd Qu.:11.12 3rd Qu.:10.416 3rd Qu.:61.0 3rd Qu.:51.1
## Max. :26.70 Max. :17.019 Max. :98.8 Max. :94.9plot_rse(df_s19,"RSE_direct","RSE_eb","Direct","EB-BK RSE-ES","S19: EB Beta-Binomial Backward — RSE ES")
S19: EB Beta Backward RSE Equal Size
S20 — Backward · Cluster
## === Skenario 20: EB Backward Cluster ===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"
)## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.574 Min. :26.1 Min. :22.1
## 1st Qu.: 4.10 1st Qu.: 4.646 1st Qu.:39.3 1st Qu.:33.1
## Median : 8.24 Median : 7.571 Median :48.4 Median :41.7
## Mean : 8.65 Mean : 7.831 Mean :52.3 Mean :41.0
## 3rd Qu.:11.12 3rd Qu.:10.308 3rd Qu.:61.0 3rd Qu.:46.4
## Max. :26.70 Max. :17.053 Max. :98.8 Max. :87.6plot_rse(df_s20,"RSE_direct","RSE_eb","Direct","EB-BK Cluster","S20: EB Beta-Binomial Backward — Cluster")
S20: EB Beta Backward Cluster-Aux
S21 — Top-n/10 · All
## === Skenario 21: EB TopN All ===df_s21 <- run_on_segments(
seg_list = list("All" = df_base),
aux_fn = seleksi_topn, model_fn = run_eb_beta, y_col = "y_logit"
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.579 Min. :26.1 Min. :24.5
## 1st Qu.: 4.10 1st Qu.: 4.418 1st Qu.:39.3 1st Qu.:34.7
## Median : 8.24 Median : 7.676 Median :48.4 Median :41.9
## Mean : 8.65 Mean : 7.833 Mean :52.3 Mean :41.8
## 3rd Qu.:11.12 3rd Qu.: 9.804 3rd Qu.:61.0 3rd Qu.:47.3
## Max. :26.70 Max. :18.702 Max. :98.8 Max. :87.3
S21: EB Beta Top-n/10 All
S22 — Top-n/10 · RSE-NB
## === Skenario 22: EB TopN RSE-NB ===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"
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.503 Min. :26.1 Min. :23.7
## 1st Qu.: 4.10 1st Qu.: 4.180 1st Qu.:39.3 1st Qu.:32.6
## Median : 8.24 Median : 8.107 Median :48.4 Median :37.0
## Mean : 8.65 Mean : 8.016 Mean :52.3 Mean :40.1
## 3rd Qu.:11.12 3rd Qu.:10.670 3rd Qu.:61.0 3rd Qu.:47.1
## Max. :26.70 Max. :17.506 Max. :98.8 Max. :90.6plot_rse(df_s22,"RSE_direct","RSE_eb","Direct","EB-TN RSE-NB","S22: EB Beta-Binomial Top-n/10 — RSE NB")
S22: EB Beta Top-n/10 RSE Natural Break
S23 — Top-n/10 · RSE-ES
## === Skenario 23: EB TopN RSE-ES ===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"
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.512 Min. :26.1 Min. :22.4
## 1st Qu.: 4.10 1st Qu.: 4.423 1st Qu.:39.3 1st Qu.:29.8
## Median : 8.24 Median : 8.259 Median :48.4 Median :39.6
## Mean : 8.65 Mean : 7.929 Mean :52.3 Mean :41.7
## 3rd Qu.:11.12 3rd Qu.:10.285 3rd Qu.:61.0 3rd Qu.:51.6
## Max. :26.70 Max. :17.124 Max. :98.8 Max. :93.6plot_rse(df_s23,"RSE_direct","RSE_eb","Direct","EB-TN RSE-ES","S23: EB Beta-Binomial Top-n/10 — RSE ES")
S23: EB Beta Top-n/10 RSE Equal Size
S24 — Top-n/10 · Cluster
## === Skenario 24: EB TopN Cluster ===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"
)## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.604 Min. :26.1 Min. :24.1
## 1st Qu.: 4.10 1st Qu.: 4.898 1st Qu.:39.3 1st Qu.:34.4
## Median : 8.24 Median : 8.065 Median :48.4 Median :39.8
## Mean : 8.65 Mean : 7.878 Mean :52.3 Mean :40.6
## 3rd Qu.:11.12 3rd Qu.: 9.915 3rd Qu.:61.0 3rd Qu.:45.6
## Max. :26.70 Max. :19.615 Max. :98.8 Max. :80.5plot_rse(df_s24,"RSE_direct","RSE_eb","Direct","EB-TN Cluster","S24: EB Beta-Binomial Top-n/10 — Cluster")
S24: EB Beta Top-n/10 Cluster-Aux
8. HB Beta-Binomial saeHB
Prasyarat JAGS: Skenario 25–32 memerlukan JAGS
terinstal di OS dan paket saeHB. Download: https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/
lalu restart R. Jika JAGS tidak tersedia, otomatis fallback ke direct
estimator. Setiap skenario HB menyertakan plot
posterior dengan credible interval 95%.
S25 — Backward · All
## === Skenario 25: HB Backward All ===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
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705
##
## [HB Iter 1] Vars aktif (2): X36_mean, X40
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.5303 0.03631 -2.6012 -2.5546 -2.5306 -2.5057 -2.4578
## X36_mean -0.2699 0.03937 -0.3460 -0.2966 -0.2703 -0.2429 -0.1930
## X40 -0.2094 0.03385 -0.2747 -0.2330 -0.2101 -0.1871 -0.1419
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean, X40## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 1.35 Min. :26.1 Min. :16.1
## 1st Qu.: 4.10 1st Qu.: 4.89 1st Qu.:39.3 1st Qu.:23.7
## Median : 8.24 Median : 8.31 Median :48.4 Median :26.5
## Mean : 8.65 Mean : 8.65 Mean :52.3 Mean :28.2
## 3rd Qu.:11.12 3rd Qu.:10.71 3rd Qu.:61.0 3rd Qu.:32.2
## Max. :26.70 Max. :24.97 Max. :98.8 Max. :46.1
S25: HB Beta Backward RSE
S25: HB Beta Posterior — All Domain
S26 — Backward · RSE-NB
## === Skenario 26: HB Backward RSE-NB ===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
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [HB Iter 1] Vars aktif (3): X3, X31_jumlah, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.2648 0.02958 -2.3236 -2.2849 -2.2646 -2.24501 -2.20726
## X3 0.1724 0.03141 0.1098 0.1515 0.1729 0.19380 0.23320
## X31_jumlah -0.1420 0.03403 -0.2087 -0.1648 -0.1422 -0.11894 -0.07459
## X29_mean -0.1208 0.03317 -0.1848 -0.1436 -0.1204 -0.09859 -0.05673
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X3, X31_jumlah, X29_mean
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.5093 0.06017 -3.6275 -3.5508 -3.5084 -3.4675 -3.3900
## X36_mean -0.5082 0.06340 -0.6311 -0.5514 -0.5084 -0.4647 -0.3852
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.514 Min. :26.1 Min. : 5.61
## 1st Qu.: 4.10 1st Qu.: 4.289 1st Qu.:39.3 1st Qu.:12.59
## Median : 8.24 Median : 8.282 Median :48.4 Median :14.91
## Mean : 8.65 Mean : 8.663 Mean :52.3 Mean :15.53
## 3rd Qu.:11.12 3rd Qu.:10.971 3rd Qu.:61.0 3rd Qu.:16.50
## Max. :26.70 Max. :25.013 Max. :98.8 Max. :49.39
S26: HB Beta Backward RSE Natural Break
S26: HB Beta Posterior — RSE Natural Break
S27 — Backward · RSE-ES
## === Skenario 27: HB Backward RSE-ES ===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
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [HB Iter 1] Vars aktif (2): X8, X39_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.1532 0.03127 -2.2171 -2.1735 -2.1527 -2.1322 -2.0942
## X8 0.2000 0.03036 0.1403 0.1793 0.2005 0.2210 0.2581
## X39_jumlah -0.1875 0.03148 -0.2501 -0.2088 -0.1878 -0.1661 -0.1263
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X8, X39_jumlah
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547
##
## [HB Iter 1] Vars aktif (1): X31_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0428 0.05493 -3.152 -3.0807 -3.0417 -3.0051 -2.9385
## X31_mean -0.5112 0.05829 -0.625 -0.5502 -0.5122 -0.4717 -0.3969
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X31_mean## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.517 Min. :26.1 Min. : 2.75
## 1st Qu.: 4.10 1st Qu.: 4.122 1st Qu.:39.3 1st Qu.: 3.88
## Median : 8.24 Median : 8.220 Median :48.4 Median : 6.87
## Mean : 8.65 Mean : 8.647 Mean :52.3 Mean : 9.07
## 3rd Qu.:11.12 3rd Qu.:11.130 3rd Qu.:61.0 3rd Qu.:12.69
## Max. :26.70 Max. :26.415 Max. :98.8 Max. :31.97
S27: HB Beta Backward RSE Equal Size
S27: HB Beta Posterior — RSE Equal Size
S28 — Backward · Cluster
## === Skenario 28: HB Backward Cluster ===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
)## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.7517 0.04742 -2.8456 -2.7831 -2.7515 -2.7196 -2.6611
## X36_mean -0.3425 0.04706 -0.4359 -0.3743 -0.3429 -0.3113 -0.2509
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562
##
## [HB Iter 1] Vars aktif (1): X40_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.0344 0.04982 -2.1328 -2.0679 -2.0338 -2.0010 -1.93689
## X40_jumlah -0.1766 0.05146 -0.2768 -0.2113 -0.1762 -0.1419 -0.07777
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X40_jumlah## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.919 Min. :26.1 Min. : 8.82
## 1st Qu.: 4.10 1st Qu.: 4.422 1st Qu.:39.3 1st Qu.:16.44
## Median : 8.24 Median : 8.011 Median :48.4 Median :20.64
## Mean : 8.65 Mean : 8.650 Mean :52.3 Mean :22.47
## 3rd Qu.:11.12 3rd Qu.:10.909 3rd Qu.:61.0 3rd Qu.:27.32
## Max. :26.70 Max. :25.845 Max. :98.8 Max. :52.09
S28: HB Beta Backward Cluster-Aux
S28: HB Beta Posterior — Cluster Aux (A1-Logit)
S29 — Top-n/10 · All
## === Skenario 29: HB TopN All ===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
)## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604
##
## [HB Iter 1] Vars aktif (7): X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.52086 0.03603 -2.59070 -2.54531 -2.52057 -2.49726 -2.44832
## X36_mean -0.03394 0.07773 -0.18760 -0.08680 -0.03411 0.01994 0.11924
## X30_mean -0.05760 0.05945 -0.16997 -0.09853 -0.05787 -0.01751 0.05904
## X31_mean -0.03315 0.07428 -0.18293 -0.08231 -0.03226 0.01765 0.10868
## X29_mean -0.12430 0.06712 -0.25771 -0.16863 -0.12321 -0.08024 0.00570
## X8 0.07834 0.05785 -0.03760 0.03964 0.07843 0.11776 0.19076
## X40 -0.20224 0.03530 -0.27020 -0.22615 -0.20230 -0.17785 -0.13389
## X45 0.02432 0.05725 -0.09039 -0.01427 0.02504 0.06315 0.13705
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X36_mean, X30_mean, X31_mean, X29_mean, X8, X45]
## [HB] → Keluarkan 'X36_mean' (SD = 0.07773, terbesar)
##
## [HB Iter 2] Vars aktif (6): X30_mean, X31_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.53143 0.03617 -2.60391 -2.55549 -2.53129 -2.50684 -2.46202
## X30_mean -0.07204 0.05442 -0.17843 -0.10885 -0.07245 -0.03517 0.03557
## X31_mean -0.03954 0.06611 -0.16885 -0.08340 -0.03988 0.00438 0.08871
## X29_mean -0.13435 0.05978 -0.24685 -0.17782 -0.13402 -0.09364 -0.01707
## X8 0.08597 0.05348 -0.01885 0.04925 0.08734 0.12333 0.18914
## X40 -0.20617 0.03555 -0.27727 -0.23032 -0.20605 -0.18171 -0.13658
## X45 0.02640 0.05597 -0.08047 -0.01201 0.02533 0.06467 0.13688
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X30_mean, X31_mean, X8, X45]
## [HB] → Keluarkan 'X31_mean' (SD = 0.06611, terbesar)
##
## [HB Iter 3] Vars aktif (5): X30_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51803 0.03611 -2.58811 -2.54255 -2.51794 -2.49341 -2.44678
## X30_mean -0.07316 0.04992 -0.17337 -0.10703 -0.07211 -0.04005 0.02331
## X29_mean -0.15256 0.05242 -0.25587 -0.18727 -0.15218 -0.11748 -0.05003
## X8 0.08885 0.05209 -0.01348 0.05357 0.08844 0.12489 0.19002
## X40 -0.20321 0.03535 -0.27205 -0.22742 -0.20321 -0.17945 -0.13342
## X45 0.01497 0.05125 -0.08697 -0.01967 0.01557 0.04955 0.11314
##
## [HB] ✗ Iter 3 — CI menyeberang nol: [X30_mean, X8, X45]
## [HB] → Keluarkan 'X8' (SD = 0.05209, terbesar)
##
## [HB Iter 4] Vars aktif (4): X30_mean, X29_mean, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51738 0.03538 -2.5869 -2.5407 -2.51701 -2.49333 -2.44780
## X30_mean -0.07768 0.04877 -0.1744 -0.1104 -0.07726 -0.04562 0.01741
## X29_mean -0.18611 0.04859 -0.2817 -0.2185 -0.18623 -0.15352 -0.09174
## X40 -0.21140 0.03441 -0.2760 -0.2349 -0.21185 -0.18849 -0.14141
## X45 -0.02148 0.04616 -0.1118 -0.0528 -0.02188 0.00971 0.06744
##
## [HB] ✗ Iter 4 — CI menyeberang nol: [X30_mean, X45]
## [HB] → Keluarkan 'X30_mean' (SD = 0.04877, terbesar)
##
## [HB Iter 5] Vars aktif (3): X29_mean, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51625 0.03513 -2.5865 -2.53987 -2.51600 -2.49271 -2.44743
## X29_mean -0.22527 0.04243 -0.3066 -0.25417 -0.22572 -0.19672 -0.13970
## X40 -0.21375 0.03495 -0.2825 -0.23731 -0.21345 -0.18945 -0.14647
## X45 -0.04969 0.04181 -0.1306 -0.07771 -0.05074 -0.02117 0.03277
##
## [HB] ✗ Iter 5 — CI menyeberang nol: [X45]
## [HB] → Keluarkan 'X45' (SD = 0.04181, terbesar)
##
## [HB Iter 6] Vars aktif (2): X29_mean, X40
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.5115 0.03618 -2.5804 -2.5356 -2.5112 -2.4873 -2.4397
## X29_mean -0.2569 0.03707 -0.3286 -0.2818 -0.2571 -0.2317 -0.1835
## X40 -0.2255 0.03274 -0.2901 -0.2474 -0.2254 -0.2035 -0.1612
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 6.
## [HB] ✓ Variabel final: X29_mean, X40## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 1.86 Min. :26.1 Min. :17.4
## 1st Qu.: 4.10 1st Qu.: 5.24 1st Qu.:39.3 1st Qu.:25.8
## Median : 8.24 Median : 8.25 Median :48.4 Median :28.0
## Mean : 8.65 Mean : 8.64 Mean :52.3 Mean :29.4
## 3rd Qu.:11.12 3rd Qu.:10.67 3rd Qu.:61.0 3rd Qu.:32.4
## Max. :26.70 Max. :24.76 Max. :98.8 Max. :42.6
S29: HB Beta Top-n/10 All
S29: HB Beta Posterior — All Domain
S30 — Top-n/10 · RSE-NB
## === Skenario 30: HB TopN RSE-NB ===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
)## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [HB Iter 1] Vars aktif (6): X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.26329 0.02909 -2.31999 -2.28290 -2.26349 -2.24370 -2.20663
## X3 0.17086 0.04404 0.08330 0.14146 0.17107 0.20021 0.25660
## X27 0.04713 0.04366 -0.03898 0.01750 0.04757 0.07592 0.13427
## X31_jumlah -0.16663 0.03700 -0.23780 -0.19141 -0.16729 -0.14199 -0.09229
## X29_mean -0.24692 0.05407 -0.35267 -0.28277 -0.24727 -0.21019 -0.14036
## X4 -0.04988 0.04328 -0.13374 -0.07891 -0.05023 -0.02031 0.03326
## X28_mean 0.14472 0.05204 0.04072 0.10944 0.14529 0.18050 0.24258
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X27, X4]
## [HB] → Keluarkan 'X27' (SD = 0.04366, terbesar)
##
## [HB Iter 2] Vars aktif (5): X3, X31_jumlah, X29_mean, X4, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.26048 0.02960 -2.31708 -2.28021 -2.26082 -2.24022 -2.20175
## X3 0.20406 0.03868 0.12916 0.17835 0.20437 0.22969 0.27910
## X31_jumlah -0.16525 0.03819 -0.24046 -0.19030 -0.16517 -0.13956 -0.09130
## X29_mean -0.27945 0.05407 -0.38528 -0.31584 -0.27947 -0.24247 -0.17510
## X4 -0.05243 0.04084 -0.13225 -0.07994 -0.05271 -0.02509 0.02757
## X28_mean 0.16013 0.05290 0.05604 0.12553 0.15983 0.19527 0.26198
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X4]
## [HB] → Keluarkan 'X4' (SD = 0.04084, terbesar)
##
## [HB Iter 3] Vars aktif (4): X3, X31_jumlah, X29_mean, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.2616 0.02819 -2.31622 -2.2805 -2.2614 -2.2428 -2.20612
## X3 0.1696 0.03147 0.10795 0.1485 0.1700 0.1915 0.22980
## X31_jumlah -0.1527 0.03621 -0.22271 -0.1772 -0.1532 -0.1281 -0.07935
## X29_mean -0.2497 0.05220 -0.35234 -0.2851 -0.2494 -0.2138 -0.14802
## X28_mean 0.1375 0.04874 0.04298 0.1046 0.1370 0.1706 0.23247
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 3.
## [HB] ✓ Variabel final: X3, X31_jumlah, X29_mean, X28_mean
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324
##
## [HB Iter 1] Vars aktif (1): X31_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.509 0.06046 -3.6248 -3.5493 -3.5102 -3.4698 -3.3874
## X31_mean -0.485 0.06442 -0.6101 -0.5289 -0.4853 -0.4407 -0.3604
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X31_mean## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.519 Min. :26.1 Min. : 5.44
## 1st Qu.: 4.10 1st Qu.: 4.354 1st Qu.:39.3 1st Qu.:14.46
## Median : 8.24 Median : 8.433 Median :48.4 Median :16.89
## Mean : 8.65 Mean : 8.658 Mean :52.3 Mean :16.74
## 3rd Qu.:11.12 3rd Qu.:11.102 3rd Qu.:61.0 3rd Qu.:18.86
## Max. :26.70 Max. :24.208 Max. :98.8 Max. :42.31
S30: HB Beta Top-n/10 RSE Natural Break
S30: HB Beta Posterior — RSE Natural Break
S31 — Top-n/10 · RSE-ES
## === Skenario 31: HB TopN RSE-ES ===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
)## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [HB Iter 1] Vars aktif (4): X44_mean, X8, X39_jumlah, X27
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.15909 0.02871 -2.21533 -2.17833 -2.15924 -2.13991 -2.1026
## X44_mean -0.02522 0.04869 -0.12133 -0.05893 -0.02485 0.01072 0.0646
## X8 0.12573 0.04327 0.04454 0.09573 0.12558 0.15592 0.2090
## X39_jumlah -0.16478 0.03217 -0.22997 -0.18602 -0.16342 -0.14259 -0.1034
## X27 0.09056 0.03520 0.02294 0.06687 0.08997 0.11399 0.1617
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X44_mean]
## [HB] → Keluarkan 'X44_mean' (SD = 0.04869, terbesar)
##
## [HB Iter 2] Vars aktif (3): X8, X39_jumlah, X27
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.15593 0.02769 -2.20940 -2.17496 -2.15688 -2.1374 -2.0992
## X8 0.13896 0.03518 0.06892 0.11451 0.13973 0.1636 0.2049
## X39_jumlah -0.17079 0.03334 -0.23518 -0.19445 -0.17096 -0.1483 -0.1050
## X27 0.09676 0.03565 0.02828 0.07256 0.09625 0.1211 0.1673
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 2.
## [HB] ✓ Variabel final: X8, X39_jumlah, X27
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315
##
## [HB Iter 1] Vars aktif (3): X27_mean, X28_mean, X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0697 0.05422 -3.1745 -3.1077 -3.0703 -3.03274 -2.96241
## X27_mean -0.2138 0.08502 -0.3804 -0.2700 -0.2166 -0.15841 -0.04124
## X28_mean -0.1123 0.08213 -0.2756 -0.1668 -0.1109 -0.05766 0.04929
## X36_mean -0.2259 0.08458 -0.3956 -0.2819 -0.2281 -0.17053 -0.05721
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X28_mean]
## [HB] → Keluarkan 'X28_mean' (SD = 0.08213, terbesar)
##
## [HB Iter 2] Vars aktif (2): X27_mean, X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0475 0.05381 -3.1519 -3.0844 -3.0467 -3.0112 -2.9428
## X27_mean -0.2905 0.06954 -0.4281 -0.3372 -0.2904 -0.2429 -0.1565
## X36_mean -0.2610 0.06894 -0.3972 -0.3054 -0.2612 -0.2153 -0.1280
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 2.
## [HB] ✓ Variabel final: X27_mean, X36_mean## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.544 Min. :26.1 Min. : 5.01
## 1st Qu.: 4.10 1st Qu.: 4.141 1st Qu.:39.3 1st Qu.: 6.84
## Median : 8.24 Median : 8.209 Median :48.4 Median : 8.61
## Mean : 8.65 Mean : 8.647 Mean :52.3 Mean :10.54
## 3rd Qu.:11.12 3rd Qu.:11.143 3rd Qu.:61.0 3rd Qu.:12.49
## Max. :26.70 Max. :26.443 Max. :98.8 Max. :31.40
S31: HB Beta Top-n/10 RSE Equal Size
S31: HB Beta Posterior — RSE Equal Size
S32 — Top-n/10 · Cluster
## === Skenario 32: HB TopN Cluster ===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
)## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [HB Iter 1] Vars aktif (7): X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41550 0.03837 -2.48810 -2.44245 -2.41547 -2.38931 -2.34054
## X40 -0.16230 0.03813 -0.23613 -0.18873 -0.16262 -0.13554 -0.08937
## X25_jumlah -0.08516 0.04158 -0.16658 -0.11356 -0.08566 -0.05612 -0.00425
## X33_mean -0.10986 0.04003 -0.18889 -0.13669 -0.10974 -0.08361 -0.02971
## X31_mean 0.05369 0.04969 -0.04501 0.02022 0.05369 0.08810 0.14863
## X48 -0.03763 0.04390 -0.12448 -0.06676 -0.03808 -0.00845 0.04829
## X35_jumlah -0.09226 0.03992 -0.17282 -0.11871 -0.09282 -0.06529 -0.01337
## X12_mean -0.08932 0.04610 -0.17796 -0.11983 -0.08986 -0.05828 0.00197
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X31_mean, X48, X12_mean]
## [HB] → Keluarkan 'X31_mean' (SD = 0.04969, terbesar)
##
## [HB Iter 2] Vars aktif (6): X40, X25_jumlah, X33_mean, X48, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41980 0.03830 -2.49566 -2.44568 -2.41982 -2.39380 -2.34585
## X40 -0.16942 0.03822 -0.24415 -0.19496 -0.16947 -0.14481 -0.09331
## X25_jumlah -0.07603 0.03895 -0.15120 -0.10224 -0.07684 -0.04978 0.00258
## X33_mean -0.10647 0.03913 -0.18291 -0.13295 -0.10579 -0.07924 -0.03172
## X48 -0.01870 0.04086 -0.09921 -0.04611 -0.01855 0.00827 0.06272
## X35_jumlah -0.08286 0.03941 -0.16018 -0.10872 -0.08318 -0.05705 -0.00545
## X12_mean -0.06330 0.04025 -0.14064 -0.09092 -0.06321 -0.03674 0.01569
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X25_jumlah, X48, X12_mean]
## [HB] → Keluarkan 'X48' (SD = 0.04086, terbesar)
##
## [HB Iter 3] Vars aktif (5): X40, X25_jumlah, X33_mean, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41679 0.03798 -2.4892 -2.44331 -2.41681 -2.39115 -2.34257
## X40 -0.17695 0.03595 -0.2471 -0.20124 -0.17709 -0.15291 -0.10543
## X25_jumlah -0.07801 0.03923 -0.1532 -0.10495 -0.07789 -0.05132 -0.00122
## X33_mean -0.10967 0.03921 -0.1871 -0.13638 -0.10943 -0.08343 -0.03333
## X35_jumlah -0.08244 0.03962 -0.1578 -0.10930 -0.08215 -0.05598 -0.00453
## X12_mean -0.06770 0.03795 -0.1435 -0.09293 -0.06803 -0.04256 0.00830
##
## [HB] ✗ Iter 3 — CI menyeberang nol: [X12_mean]
## [HB] → Keluarkan 'X12_mean' (SD = 0.03795, terbesar)
##
## [HB Iter 4] Vars aktif (4): X40, X25_jumlah, X33_mean, X35_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41544 0.03879 -2.4920 -2.4418 -2.41503 -2.38945 -2.33925
## X40 -0.17219 0.03627 -0.2450 -0.1967 -0.17222 -0.14747 -0.10130
## X25_jumlah -0.09541 0.03865 -0.1716 -0.1221 -0.09494 -0.06954 -0.02069
## X33_mean -0.12129 0.03856 -0.1960 -0.1480 -0.12224 -0.09581 -0.04437
## X35_jumlah -0.09018 0.03909 -0.1685 -0.1162 -0.08994 -0.06373 -0.01535
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 4.
## [HB] ✓ Variabel final: X40, X25_jumlah, X33_mean, X35_jumlah
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.2833 0.1056 -3.4891 -3.3540 -3.283 -3.213 -3.0710
## X36_mean -0.4585 0.1098 -0.6738 -0.5316 -0.460 -0.385 -0.2457
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.782 Min. :26.1 Min. :17.1
## 1st Qu.: 4.10 1st Qu.: 4.924 1st Qu.:39.3 1st Qu.:24.1
## Median : 8.24 Median : 8.409 Median :48.4 Median :26.5
## Mean : 8.65 Mean : 8.659 Mean :52.3 Mean :28.0
## 3rd Qu.:11.12 3rd Qu.:10.563 3rd Qu.:61.0 3rd Qu.:30.4
## Max. :26.70 Max. :24.867 Max. :98.8 Max. :46.9
S32: HB Beta Top-n/10 Cluster-Aux
S32: HB Beta Posterior — Cluster Aux (A2-Logit)
9. Skenario Refined S33–S64 · RSE > 25% Re-modelled
Skenario 33–64 adalah versi refined dari S1–S32. Setelah run model utama, domain dengan RSE > 25% di-pool per grup asalnya dan dimodelkan ulang dengan seleksi aux baru (aux bisa berbeda dari model utama). Estimasi dan RSE akhir = gabungan domain ≤25% (dari model utama) + domain >25% (dari re-run). Jika domain >25% < 4 dalam satu grup, fallback ke hasil model utama.
# ══════════════════════════════════════════════════════════════════════════════
# run_with_refine(): wrapper run_on_segments + 1x refinement domain RSE > 25%
#
# Alur:
# [1] run_on_segments() normal → df_result
# [2] Per segmen dalam seg_list: cari domain RSE > rse_threshold
# [3] Jika cukup (≥ min_n_refine): ambil baris ASLI dari df_orig, re-run
# dengan aux_fn BARU (seleksi ulang pada subset high-RSE)
# [4] Update df_result: ganti kolom model-output domain bad dengan hasil re-run
# [5] Return df_result final
#
# Catatan:
# - "tetap split per grup asal": tiap segmen diproses mandiri
# - model_cols dideteksi otomatis = kolom di df_result yg tidak ada di df_orig
# (menghindari overwrite kolom aux yang sudah di-scale)
# ══════════════════════════════════════════════════════════════════════════════
run_with_refine <- function(seg_list, aux_fn, model_fn, y_col,
rse_col, df_orig,
rse_threshold = 25,
min_n_refine = 4,
min_n = 4,
...) {
# ── [1] Run model utama ────────────────────────────────────────────────────
cat(" >>> Tahap 1: model utama\n")
df_result <- run_on_segments(seg_list, aux_fn, model_fn, y_col,
min_n = min_n, ...)
if (is.null(df_result)) return(NULL)
# Kolom model = semua kolom baru yang ditambahkan model (bukan dari df_orig)
model_cols <- setdiff(names(df_result), names(df_orig))
# ── [2-4] Per segmen: identifikasi & refine domain high-RSE ───────────────
cat(" >>> Tahap 2: refinement domain RSE >", rse_threshold, "%\n")
for (nm in names(seg_list)) {
kako_seg <- seg_list[[nm]]$Kako
idx_seg <- which(df_result$Kako %in% kako_seg)
if (!rse_col %in% names(df_result)) {
cat(sprintf(" [Refine|%s] Kolom '%s' tidak ada → skip\n", nm, rse_col))
next
}
rse_vals <- df_result[[rse_col]][idx_seg]
bad_mask <- !is.na(rse_vals) & rse_vals > rse_threshold
bad_kako <- df_result$Kako[idx_seg[bad_mask]]
cat(sprintf(" [Refine|%s] Domain RSE > %g%%: %d / %d\n",
nm, rse_threshold, length(bad_kako), length(idx_seg)))
if (length(bad_kako) < min_n_refine) {
cat(sprintf(" [Refine|%s] < %d domain → skip, pakai hasil awal\n",
nm, min_n_refine))
next
}
# [3] Subset ASLI dari df_orig (belum di-scale)
df_bad <- df_orig[df_orig$Kako %in% bad_kako, ]
cat(sprintf(" [Refine|%s] Re-run pada %d domain (aux dipilih ulang) ...\n",
nm, nrow(df_bad)))
df_ref <- tryCatch(
run_on_segments(
seg_list = setNames(list(df_bad), paste0(nm, "_ref")),
aux_fn = aux_fn,
model_fn = model_fn,
y_col = y_col,
min_n = min_n,
...
),
error = function(e) {
cat(sprintf(" [Refine|%s] Error: %s → pakai hasil awal\n", nm, e$message))
NULL
}
)
if (is.null(df_ref)) next
# [4] Update df_result: ganti model_cols untuk domain bad
upd_kako <- intersect(df_ref$Kako, bad_kako)
upd_res <- match(upd_kako, df_result$Kako) # posisi di df_result
upd_ref <- match(upd_kako, df_ref$Kako) # posisi di df_ref
for (col in model_cols) {
if (col %in% names(df_ref)) {
df_result[[col]][upd_res] <- df_ref[[col]][upd_ref]
}
}
cat(sprintf(" [Refine|%s] ✓ %d domain diperbarui\n", nm, length(upd_kako)))
}
df_result
}EBLUP Refined
S33 — Backward · All
## === S33: EBLUP Backward All + Refine ===df_s33 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_backward,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] Estimasi | n = 81 | Vars: X40, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3690 -0.3705
## X29_mean -0.3220 -0.3921
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 67 / 81
## [Refine|All] Re-run pada 67 domain (aux dipilih ulang) ...
## [ All_ref ] Estimasi | n = 67 | Vars: X31_jumlah, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_jumlah -0.3016 -0.2975
## X45 -0.2999 -0.3190
##
## [Refine|All] ✓ 67 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.594 Min. :26.1 Min. :20.7
## 1st Qu.: 4.10 1st Qu.: 4.677 1st Qu.:39.3 1st Qu.:30.2
## Median : 8.24 Median : 6.321 Median :48.4 Median :34.8
## Mean : 8.65 Mean : 6.260 Mean :52.3 Mean :35.8
## 3rd Qu.:11.12 3rd Qu.: 7.382 3rd Qu.:61.0 3rd Qu.:39.4
## Max. :26.70 Max. :14.102 Max. :98.8 Max. :83.9plot_rse(df_s33,"RSE_direct","RSE_eblup","Direct","EBLUP-BK-R All","S33: EBLUP Backward All — Refined")
S33: EBLUP Backward All + Refine
S34 — Backward · RSE-NB
## === S34: EBLUP Backward RSE-NB + Refine ===df_s34 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_backward,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] Estimasi | n = 60 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828## [ G2_RSE_Tinggi ] Estimasi | n = 21 | Vars: X41_jumlah, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.4969 0.5170
## X28_mean -0.4481 -0.6093
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 60 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 60 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] Estimasi | n = 60 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828## [Refine|G1_RSE_Rendah] ✓ 60 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 16 / 21
## [Refine|G2_RSE_Tinggi] Re-run pada 16 domain (aux dipilih ulang) ...
## [ G2_RSE_Tinggi_ref ] Estimasi | n = 16 | Vars: X41_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.6334 0.6140
##
## [Refine|G2_RSE_Tinggi] ✓ 16 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.745 Min. :26.1 Min. :19.6
## 1st Qu.: 4.10 1st Qu.:4.885 1st Qu.:39.3 1st Qu.:55.4
## Median : 8.24 Median :8.250 Median :48.4 Median :55.4
## Mean : 8.65 Mean :6.675 Mean :52.3 Mean :50.0
## 3rd Qu.:11.12 3rd Qu.:8.250 3rd Qu.:61.0 3rd Qu.:55.4
## Max. :26.70 Max. :8.250 Max. :98.8 Max. :98.9plot_rse(df_s34,"RSE_direct","RSE_eblup","Direct","EBLUP-BK-R RSE-NB","S34: EBLUP Backward RSE-NB — Refined")
S34: EBLUP Backward RSE-NB + Refine
S35 — Backward · RSE-ES
## === S35: EBLUP Backward RSE-ES + Refine ===df_s35 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_backward,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] Estimasi | n = 41 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4689 0.4949
## X31_jumlah -0.4520 -0.5349## [ G2_RSE_Atas ] Estimasi | n = 40 | Vars: X40_mean, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_mean -0.4352 -0.3638
## X28_mean -0.4137 -0.5322
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 41 / 41
## [Refine|G1_RSE_Bawah] Re-run pada 41 domain (aux dipilih ulang) ...
## [ G1_RSE_Bawah_ref ] Estimasi | n = 41 | Vars: X3, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4689 0.4949
## X31_jumlah -0.4520 -0.5349## [Refine|G1_RSE_Bawah] ✓ 41 domain diperbarui
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 9 / 40
## [Refine|G2_RSE_Atas] Re-run pada 9 domain (aux dipilih ulang) ...
## [ G2_RSE_Atas_ref ] Estimasi | n = 9 | Vars: X1_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X1_mean -0.8859 -0.8874
##
## [Refine|G2_RSE_Atas] ✓ 9 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.589 Min. :26.1 Min. : 15.1
## 1st Qu.: 4.10 1st Qu.:2.604 1st Qu.:39.3 1st Qu.: 20.6
## Median : 8.24 Median :7.940 Median :48.4 Median : 42.0
## Mean : 8.65 Mean :5.270 Mean :52.3 Mean : 35.1
## 3rd Qu.:11.12 3rd Qu.:7.940 3rd Qu.:61.0 3rd Qu.: 42.0
## Max. :26.70 Max. :7.940 Max. :98.8 Max. :104.2plot_rse(df_s35,"RSE_direct","RSE_eblup","Direct","EBLUP-BK-R RSE-ES","S35: EBLUP Backward RSE-ES — Refined")
S35: EBLUP Backward RSE-ES + Refine
S36 — Backward · Cluster
## === S36: EBLUP Backward Cluster + Refine ===df_s36 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_bk_normal), aux_fn = seleksi_backward,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] Estimasi | n = 60 | Vars: X47, X8
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.3911 0.3619
## X8 0.3176 0.3266
##
## [ Cluster_2 ] Estimasi | n = 21 | Vars: X35, X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X35 -0.3773 -0.2798
## X40_jumlah -0.3419 -0.3562
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 60 / 60
## [Refine|Cluster_1] Re-run pada 60 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] Estimasi | n = 60 | Vars: X47, X8
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.3911 0.3619
## X8 0.3176 0.3266
##
## [Refine|Cluster_1] ✓ 60 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 8 / 21
## [Refine|Cluster_2] Re-run pada 8 domain (aux dipilih ulang) ...
## [ Cluster_2_ref ] Estimasi | n = 8 | Vars: X6_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X6_jumlah 0.8317 0.7331
##
## [Refine|Cluster_2] ✓ 8 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.636 Min. :26.1 Min. :20.6
## 1st Qu.: 4.10 1st Qu.: 3.973 1st Qu.:39.3 1st Qu.:29.5
## Median : 8.24 Median : 6.411 Median :48.4 Median :35.1
## Mean : 8.65 Mean : 6.407 Mean :52.3 Mean :37.0
## 3rd Qu.:11.12 3rd Qu.: 7.739 3rd Qu.:61.0 3rd Qu.:43.1
## Max. :26.70 Max. :15.818 Max. :98.8 Max. :78.3plot_rse(df_s36,"RSE_direct","RSE_eblup","Direct","EBLUP-BK-R Cluster","S36: EBLUP Backward Cluster — Refined")
S36: EBLUP Backward Cluster + Refine
S37 — Top-n/10 · All
## === S37: EBLUP TopN All + Refine ===df_s37 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_topn,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] Estimasi | n = 81 | Vars: X40, X8, X45, X36_mean, X29_mean, X27, X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3690 -0.3705
## X8 0.3439 0.3797
## X45 -0.3224 -0.3604
## X36_mean -0.3220 -0.4311
## X29_mean -0.3220 -0.3921
## X27 0.3197 0.3124
## X31_mean -0.3188 -0.3953
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 73 / 81
## [Refine|All] Re-run pada 73 domain (aux dipilih ulang) ...
## [ All_ref ] Estimasi | n = 73 | Vars: X8_mean, X44_mean, X31_jumlah, X29_mean, X36_mean, X4, X3
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8_mean 0.3481 0.3659
## X44_mean -0.3225 -0.3743
## X31_jumlah -0.3167 -0.3129
## X29_mean -0.3130 -0.3883
## X36_mean -0.3101 -0.4189
## X4 0.2954 0.2928
## X3 0.2935 0.2333
##
## [Refine|All] ✓ 73 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.559 Min. :26.1 Min. :22.6
## 1st Qu.: 4.10 1st Qu.: 4.628 1st Qu.:39.3 1st Qu.:30.4
## Median : 8.24 Median : 6.886 Median :48.4 Median :37.4
## Mean : 8.65 Mean : 6.602 Mean :52.3 Mean :38.3
## 3rd Qu.:11.12 3rd Qu.: 8.198 3rd Qu.:61.0 3rd Qu.:42.2
## Max. :26.70 Max. :13.824 Max. :98.8 Max. :90.0plot_rse(df_s37,"RSE_direct","RSE_eblup","Direct","EBLUP-TN-R All","S37: EBLUP Top-n/10 All — Refined")
S37: EBLUP Top-n/10 All + Refine
S38 — Top-n/10 · RSE-NB
## === S38: EBLUP TopN RSE-NB + Refine ===df_s38 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_topn,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] Estimasi | n = 60 | Vars: X3, X4, X27, X31_jumlah, X29_mean, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X4 0.4425 0.4608
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X39_jumlah -0.3610 -0.3605## [ G2_RSE_Tinggi ] Estimasi | n = 21 | Vars: X41_jumlah, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.4969 0.5170
## X28_mean -0.4481 -0.6093
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 60 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 60 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] Estimasi | n = 60 | Vars: X3, X4, X27, X31_jumlah, X29_mean, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X4 0.4425 0.4608
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X39_jumlah -0.3610 -0.3605## [Refine|G1_RSE_Rendah] ✓ 60 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 16 / 21
## [Refine|G2_RSE_Tinggi] Re-run pada 16 domain (aux dipilih ulang) ...
## [ G2_RSE_Tinggi_ref ] Estimasi | n = 16 | Vars: X41_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X41_jumlah 0.6334 0.6140
##
## [Refine|G2_RSE_Tinggi] ✓ 16 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.745 Min. :26.1 Min. :19.6
## 1st Qu.: 4.10 1st Qu.:4.885 1st Qu.:39.3 1st Qu.:55.4
## Median : 8.24 Median :8.250 Median :48.4 Median :55.4
## Mean : 8.65 Mean :6.675 Mean :52.3 Mean :50.0
## 3rd Qu.:11.12 3rd Qu.:8.250 3rd Qu.:61.0 3rd Qu.:55.4
## Max. :26.70 Max. :8.250 Max. :98.8 Max. :98.9plot_rse(df_s38,"RSE_direct","RSE_eblup","Direct","EBLUP-TN-R RSE-NB","S38: EBLUP Top-n/10 RSE-NB — Refined")
S38: EBLUP Top-n/10 RSE-NB + Refine
S39 — Top-n/10 · RSE-ES
## === S39: EBLUP TopN RSE-ES + Refine ===df_s39 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_topn,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] Estimasi | n = 41 | Vars: X4, X3, X27, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X4 0.4782 0.5323
## X3 0.4689 0.4949
## X27 0.4608 0.5500
## X31_jumlah -0.4520 -0.5349## [ G2_RSE_Atas ] Estimasi | n = 40 | Vars: X25_mean, X47, X40_mean, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X25_mean -0.4527 -0.5200
## X47 0.4523 0.4163
## X40_mean -0.4352 -0.3638
## X28_mean -0.4137 -0.5322
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 41 / 41
## [Refine|G1_RSE_Bawah] Re-run pada 41 domain (aux dipilih ulang) ...
## [ G1_RSE_Bawah_ref ] Estimasi | n = 41 | Vars: X4, X3, X27, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X4 0.4782 0.5323
## X3 0.4689 0.4949
## X27 0.4608 0.5500
## X31_jumlah -0.4520 -0.5349## [Refine|G1_RSE_Bawah] ✓ 41 domain diperbarui
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 18 / 40
## [Refine|G2_RSE_Atas] Re-run pada 18 domain (aux dipilih ulang) ...
## [ G2_RSE_Atas_ref ] Estimasi | n = 18 | Vars: X1
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X1 -0.7871 -0.8524
##
## [Refine|G2_RSE_Atas] ✓ 18 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. :0.46 Min. :26.1 Min. : 13.2
## 1st Qu.: 4.10 1st Qu.:3.23 1st Qu.:39.3 1st Qu.: 24.1
## Median : 8.24 Median :7.94 Median :48.4 Median : 42.0
## Mean : 8.65 Mean :5.91 Mean :52.3 Mean : 36.6
## 3rd Qu.:11.12 3rd Qu.:7.94 3rd Qu.:61.0 3rd Qu.: 42.0
## Max. :26.70 Max. :9.33 Max. :98.8 Max. :135.7plot_rse(df_s39,"RSE_direct","RSE_eblup","Direct","EBLUP-TN-R RSE-ES","S39: EBLUP Top-n/10 RSE-ES — Refined")
S39: EBLUP Top-n/10 RSE-ES + Refine
S40 — Top-n/10 · Cluster
## === S40: EBLUP TopN Cluster + Refine ===df_s40 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_tn_normal), aux_fn = seleksi_topn,
model_fn = run_eblup, y_col = "Estimasi",
rse_col = "RSE_eblup", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] Estimasi | n = 11 | Vars: X47
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.7038 0.4875
##
## [ Cluster_2 ] Estimasi | n = 70 | Vars: X40, X48, X8, X35_jumlah, X45_jumlah, X44_jumlah, X3
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3196 -0.3335
## X48 -0.2726 -0.2762
## X8 0.2662 0.3094
## X35_jumlah -0.2568 -0.2167
## X45_jumlah -0.2369 -0.1636
## X44_jumlah -0.2345 -0.2020
## X3 0.2332 0.1541
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 11 / 11
## [Refine|Cluster_1] Re-run pada 11 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] Estimasi | n = 11 | Vars: X47
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X47 0.7038 0.4875
##
## [Refine|Cluster_1] ✓ 11 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 63 / 70
## [Refine|Cluster_2] Re-run pada 63 domain (aux dipilih ulang) ...
## [ Cluster_2_ref ] Estimasi | n = 63 | Vars: X40_mean, X8_mean, X3, X48, X31_mean, X45_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_mean -0.3335 -0.3164
## X8_mean 0.2564 0.2951
## X3 0.2497 0.1638
## X48 -0.2386 -0.2380
## X31_mean -0.2339 -0.2743
## X45_jumlah -0.2299 -0.1366
##
## [Refine|Cluster_2] ✓ 63 domain diperbarui## Estimasi y_eblup RSE_direct RSE_eblup
## Min. : 0.51 Min. : 0.718 Min. :26.1 Min. :22.3
## 1st Qu.: 4.10 1st Qu.: 4.280 1st Qu.:39.3 1st Qu.:28.7
## Median : 8.24 Median : 6.759 Median :48.4 Median :38.6
## Mean : 8.65 Mean : 6.626 Mean :52.3 Mean :38.6
## 3rd Qu.:11.12 3rd Qu.: 8.843 3rd Qu.:61.0 3rd Qu.:45.8
## Max. :26.70 Max. :15.059 Max. :98.8 Max. :72.6plot_rse(df_s40,"RSE_direct","RSE_eblup","Direct","EBLUP-TN-R Cluster","S40: EBLUP Top-n/10 Cluster — Refined")
S40: EBLUP Top-n/10 Cluster + Refine
GLMM Refined
S41 — Backward · All
## === S41: GLMM Backward All + Refine ===df_s41 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_backward,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 30 / 81
## [Refine|All] Re-run pada 30 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 30 | Vars: X31_jumlah, X33_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_jumlah -0.4510 -0.5430
## X33_mean -0.4403 -0.5255
##
## [Refine|All] ✓ 30 domain diperbarui## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.606 Min. :26.1 Min. : 1.18
## 1st Qu.: 4.10 1st Qu.: 4.763 1st Qu.:39.3 1st Qu.: 8.10
## Median : 8.24 Median : 7.491 Median :48.4 Median : 15.83
## Mean : 8.65 Mean : 7.643 Mean :52.3 Mean : 23.21
## 3rd Qu.:11.12 3rd Qu.: 8.904 3rd Qu.:61.0 3rd Qu.: 24.12
## Max. :26.70 Max. :18.196 Max. :98.8 Max. :121.33plot_rse(df_s41,"RSE_direct","RSE_glmm","Direct","GLMM-BK-R All","S41: GLMM Backward All — Refined")
S41: GLMM Backward All + Refine
S42 — Backward · RSE-NB
## === S42: GLMM Backward RSE-NB + Refine ===df_s42 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_backward,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 4 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 4 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] y_logit | n = 4 | Vars: X48
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X48 -0.9524 -0.9788
##
## [Refine|G1_RSE_Rendah] ✓ 4 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 0 / 21
## [Refine|G2_RSE_Tinggi] < 4 domain → skip, pakai hasil awal## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.567 Min. :26.1 Min. : 2.35
## 1st Qu.: 4.10 1st Qu.: 3.355 1st Qu.:39.3 1st Qu.: 2.91
## Median : 8.24 Median : 7.865 Median :48.4 Median : 5.74
## Mean : 8.65 Mean : 7.611 Mean :52.3 Mean : 7.95
## 3rd Qu.:11.12 3rd Qu.:10.369 3rd Qu.:61.0 3rd Qu.:11.81
## Max. :26.70 Max. :19.359 Max. :98.8 Max. :24.93plot_rse(df_s42,"RSE_direct","RSE_glmm","Direct","GLMM-BK-R RSE-NB","S42: GLMM Backward RSE-NB — Refined")
S42: GLMM Backward RSE-NB + Refine
S43 — Backward · RSE-ES
## === S43: GLMM Backward RSE-ES + Refine ===df_s43 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_backward,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 0 / 41
## [Refine|G1_RSE_Bawah] < 4 domain → skip, pakai hasil awal
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 4 / 40
## [Refine|G2_RSE_Atas] Re-run pada 4 domain (aux dipilih ulang) ...
## [ G2_RSE_Atas_ref ] y_logit | n = 4 | Vars: X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X28_mean -0.9506 -0.9817
##
## [Refine|G2_RSE_Atas] ✓ 4 domain diperbarui## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.417 Min. :26.1 Min. : 0.843
## 1st Qu.: 4.10 1st Qu.: 4.191 1st Qu.:39.3 1st Qu.: 1.036
## Median : 8.24 Median : 6.490 Median :48.4 Median : 3.031
## Mean : 8.65 Mean : 7.266 Mean :52.3 Mean : 6.319
## 3rd Qu.:11.12 3rd Qu.:10.685 3rd Qu.:61.0 3rd Qu.:10.378
## Max. :26.70 Max. :13.162 Max. :98.8 Max. :22.858plot_rse(df_s43,"RSE_direct","RSE_glmm","Direct","GLMM-BK-R RSE-ES","S43: GLMM Backward RSE-ES — Refined")
S43: GLMM Backward RSE-ES + Refine
S44 — Backward · Cluster
## === S44: GLMM Backward Cluster + Refine ===df_s44 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_bk_logit), aux_fn = seleksi_backward,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 29 / 60
## [Refine|Cluster_1] Re-run pada 29 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 29 | Vars: X33_mean, X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X33_mean -0.3775 -0.4178
## X31_mean -0.2626 -0.4096
##
## [Refine|Cluster_1] ✓ 29 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 0 / 21
## [Refine|Cluster_2] < 4 domain → skip, pakai hasil awal## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.609 Min. :26.1 Min. : 1.83
## 1st Qu.: 4.10 1st Qu.: 4.569 1st Qu.:39.3 1st Qu.: 3.62
## Median : 8.24 Median : 7.536 Median :48.4 Median : 15.85
## Mean : 8.65 Mean : 7.713 Mean :52.3 Mean : 26.13
## 3rd Qu.:11.12 3rd Qu.:10.896 3rd Qu.:61.0 3rd Qu.: 43.64
## Max. :26.70 Max. :18.077 Max. :98.8 Max. :211.26plot_rse(df_s44,"RSE_direct","RSE_glmm","Direct","GLMM-BK-R Cluster","S44: GLMM Backward Cluster — Refined")
S44: GLMM Backward Cluster + Refine
S45 — Top-n/10 · All
## === S45: GLMM TopN All + Refine ===df_s45 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_topn,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 30 / 81
## [Refine|All] Re-run pada 30 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 30 | Vars: X33_mean, X36_jumlah, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X33_mean -0.4148 -0.5075
## X36_jumlah -0.3662 -0.4857
## X31_jumlah -0.4011 -0.4710
##
## [Refine|All] ✓ 30 domain diperbarui## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.77 Min. :26.1 Min. : 1.77
## 1st Qu.: 4.10 1st Qu.: 4.87 1st Qu.:39.3 1st Qu.: 9.33
## Median : 8.24 Median : 7.38 Median :48.4 Median : 17.18
## Mean : 8.65 Mean : 7.72 Mean :52.3 Mean : 26.94
## 3rd Qu.:11.12 3rd Qu.: 9.09 3rd Qu.:61.0 3rd Qu.: 33.27
## Max. :26.70 Max. :19.35 Max. :98.8 Max. :136.25plot_rse(df_s45,"RSE_direct","RSE_glmm","Direct","GLMM-TN-R All","S45: GLMM Top-n/10 All — Refined")
S45: GLMM Top-n/10 All + Refine
S46 — Top-n/10 · RSE-NB
## === S46: GLMM TopN RSE-NB + Refine ===df_s46 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_topn,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 4 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 4 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] y_logit | n = 4 | Vars: X48
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X48 -0.9524 -0.9788
##
## [Refine|G1_RSE_Rendah] ✓ 4 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 0 / 21
## [Refine|G2_RSE_Tinggi] < 4 domain → skip, pakai hasil awal## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.605 Min. :26.1 Min. : 2.00
## 1st Qu.: 4.10 1st Qu.: 3.307 1st Qu.:39.3 1st Qu.: 2.96
## Median : 8.24 Median : 7.698 Median :48.4 Median : 5.86
## Mean : 8.65 Mean : 7.574 Mean :52.3 Mean : 7.77
## 3rd Qu.:11.12 3rd Qu.:10.504 3rd Qu.:61.0 3rd Qu.:10.54
## Max. :26.70 Max. :19.359 Max. :98.8 Max. :24.82plot_rse(df_s46,"RSE_direct","RSE_glmm","Direct","GLMM-TN-R RSE-NB","S46: GLMM Top-n/10 RSE-NB — Refined")
S46: GLMM Top-n/10 RSE-NB + Refine
S47 — Top-n/10 · RSE-ES
## === S47: GLMM TopN RSE-ES + Refine ===df_s47 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_topn,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 0 / 41
## [Refine|G1_RSE_Bawah] < 4 domain → skip, pakai hasil awal
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 3 / 40
## [Refine|G2_RSE_Atas] < 4 domain → skip, pakai hasil awal## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 0.65 Min. :26.1 Min. : 1.04
## 1st Qu.: 4.10 1st Qu.: 4.51 1st Qu.:39.3 1st Qu.: 1.55
## Median : 8.24 Median : 6.58 Median :48.4 Median : 3.71
## Mean : 8.65 Mean : 7.36 Mean :52.3 Mean : 8.28
## 3rd Qu.:11.12 3rd Qu.:10.90 3rd Qu.:61.0 3rd Qu.:13.87
## Max. :26.70 Max. :14.75 Max. :98.8 Max. :58.20plot_rse(df_s47,"RSE_direct","RSE_glmm","Direct","GLMM-TN-R RSE-ES","S47: GLMM Top-n/10 RSE-ES — Refined")
S47: GLMM Top-n/10 RSE-ES + Refine
S48 — Top-n/10 · Cluster
## === S48: GLMM TopN Cluster + Refine ===df_s48 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_tn_logit), aux_fn = seleksi_topn,
model_fn = run_glmm, y_col = "y_logit",
rse_col = "RSE_glmm", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 23 / 70
## [Refine|Cluster_1] Re-run pada 23 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 23 | Vars: X33_mean, X44_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X33_mean -0.3129 -0.4100
## X44_jumlah -0.3119 -0.3675
##
## [Refine|Cluster_1] ✓ 23 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 1 / 11
## [Refine|Cluster_2] < 4 domain → skip, pakai hasil awal## Estimasi y_glmm RSE_direct RSE_glmm
## Min. : 0.51 Min. : 1.12 Min. :26.1 Min. : 2.05
## 1st Qu.: 4.10 1st Qu.: 4.40 1st Qu.:39.3 1st Qu.: 7.17
## Median : 8.24 Median : 7.48 Median :48.4 Median : 16.43
## Mean : 8.65 Mean : 7.73 Mean :52.3 Mean : 28.24
## 3rd Qu.:11.12 3rd Qu.: 9.27 3rd Qu.:61.0 3rd Qu.: 27.86
## Max. :26.70 Max. :20.13 Max. :98.8 Max. :280.17plot_rse(df_s48,"RSE_direct","RSE_glmm","Direct","GLMM-TN-R Cluster","S48: GLMM Top-n/10 Cluster — Refined")
S48: GLMM Top-n/10 Cluster + Refine
EB Beta Refined
S49 — Backward · All
## === S49: EB Backward All + Refine ===df_s49 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_backward,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 80 / 81
## [Refine|All] Re-run pada 80 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 80 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3404 -0.4446
##
## [Refine|All] ✓ 80 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.564 Min. :26.1 Min. :24.6
## 1st Qu.: 4.10 1st Qu.: 5.054 1st Qu.:39.3 1st Qu.:35.0
## Median : 8.24 Median : 7.692 Median :48.4 Median :41.8
## Mean : 8.65 Mean : 7.702 Mean :52.3 Mean :41.7
## 3rd Qu.:11.12 3rd Qu.: 9.445 3rd Qu.:61.0 3rd Qu.:46.4
## Max. :26.70 Max. :16.532 Max. :98.8 Max. :88.4
S49: EB Beta Backward All + Refine
S50 — Backward · RSE-NB
## === S50: EB Backward RSE-NB + Refine ===df_s50 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_backward,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 58 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 58 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] y_logit | n = 58 | Vars: X3, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4986 0.4995
## X29_mean -0.3950 -0.4783
##
## [Refine|G1_RSE_Rendah] ✓ 58 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 21 / 21
## [Refine|G2_RSE_Tinggi] Re-run pada 21 domain (aux dipilih ulang) ...
## [ G2_RSE_Tinggi_ref ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631
##
## [Refine|G2_RSE_Tinggi] ✓ 21 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.497 Min. :26.1 Min. :24.0
## 1st Qu.: 4.10 1st Qu.: 4.043 1st Qu.:39.3 1st Qu.:32.3
## Median : 8.24 Median : 8.029 Median :48.4 Median :37.1
## Mean : 8.65 Mean : 7.964 Mean :52.3 Mean :40.1
## 3rd Qu.:11.12 3rd Qu.:10.645 3rd Qu.:61.0 3rd Qu.:46.8
## Max. :26.70 Max. :17.764 Max. :98.8 Max. :91.2plot_rse(df_s50,"RSE_direct","RSE_eb","Direct","EB-BK-R RSE-NB","S50: EB Beta Backward RSE-NB — Refined")
S50: EB Beta Backward RSE-NB + Refine
S51 — Backward · RSE-ES
## === S51: EB Backward RSE-ES + Refine ===df_s51 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_backward,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 39 / 41
## [Refine|G1_RSE_Bawah] Re-run pada 39 domain (aux dipilih ulang) ...
## [ G1_RSE_Bawah_ref ] y_logit | n = 39 | Vars: X44, X4
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44 -0.4590 -0.6038
## X4 0.5305 0.5801
##
## [Refine|G1_RSE_Bawah] ✓ 39 domain diperbarui
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 40 / 40
## [Refine|G2_RSE_Atas] Re-run pada 40 domain (aux dipilih ulang) ...
## [ G2_RSE_Atas_ref ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547
##
## [Refine|G2_RSE_Atas] ✓ 40 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.498 Min. :26.1 Min. :22.1
## 1st Qu.: 4.10 1st Qu.: 4.314 1st Qu.:39.3 1st Qu.:29.4
## Median : 8.24 Median : 8.169 Median :48.4 Median :40.0
## Mean : 8.65 Mean : 7.866 Mean :52.3 Mean :41.8
## 3rd Qu.:11.12 3rd Qu.:10.491 3rd Qu.:61.0 3rd Qu.:51.1
## Max. :26.70 Max. :17.641 Max. :98.8 Max. :94.9plot_rse(df_s51,"RSE_direct","RSE_eb","Direct","EB-BK-R RSE-ES","S51: EB Beta Backward RSE-ES — Refined")
S51: EB Beta Backward RSE-ES + Refine
S52 — Backward · Cluster
## === S52: EB Backward Cluster + Refine ===df_s52 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_bk_logit), aux_fn = seleksi_backward,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 60 / 60
## [Refine|Cluster_1] Re-run pada 60 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [Refine|Cluster_1] ✓ 60 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 20 / 21
## [Refine|Cluster_2] Re-run pada 20 domain (aux dipilih ulang) ...
## [ Cluster_2_ref ] y_logit | n = 20 | Vars: X9, X45_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X9 0.2909 0.3747
## X45_jumlah -0.3457 -0.3734
##
## [Refine|Cluster_2] ✓ 20 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.574 Min. :26.1 Min. :22.1
## 1st Qu.: 4.10 1st Qu.: 4.646 1st Qu.:39.3 1st Qu.:32.6
## Median : 8.24 Median : 7.447 Median :48.4 Median :41.7
## Mean : 8.65 Mean : 7.831 Mean :52.3 Mean :41.1
## 3rd Qu.:11.12 3rd Qu.:10.301 3rd Qu.:61.0 3rd Qu.:46.4
## Max. :26.70 Max. :17.512 Max. :98.8 Max. :87.6plot_rse(df_s52,"RSE_direct","RSE_eb","Direct","EB-BK-R Cluster","S52: EB Beta Backward Cluster — Refined")
S52: EB Beta Backward Cluster + Refine
S53 — Top-n/10 · All
## === S53: EB TopN All + Refine ===df_s53 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_topn,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 80 / 81
## [Refine|All] Re-run pada 80 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 80 | Vars: X36_mean, X30_mean, X29_mean, X31_mean, X45_mean, X8, X27_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3404 -0.4446
## X30_mean -0.3213 -0.4098
## X29_mean -0.3455 -0.4089
## X31_mean -0.3307 -0.4035
## X45_mean -0.3683 -0.3942
## X8 0.3418 0.3772
## X27_mean -0.3016 -0.3754
##
## [Refine|All] ✓ 80 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.591 Min. :26.1 Min. :24.5
## 1st Qu.: 4.10 1st Qu.: 4.950 1st Qu.:39.3 1st Qu.:35.0
## Median : 8.24 Median : 7.814 Median :48.4 Median :41.4
## Mean : 8.65 Mean : 7.727 Mean :52.3 Mean :41.7
## 3rd Qu.:11.12 3rd Qu.: 9.650 3rd Qu.:61.0 3rd Qu.:45.8
## Max. :26.70 Max. :16.619 Max. :98.8 Max. :86.3
S53: EB Beta Top-n/10 All + Refine
S54 — Top-n/10 · RSE-NB
## === S54: EB TopN RSE-NB + Refine ===df_s54 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_topn,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 58 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 58 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] y_logit | n = 58 | Vars: X3, X27, X29_mean, X4, X31_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4986 0.4995
## X27 0.4386 0.4867
## X29_mean -0.3950 -0.4783
## X4 0.4568 0.4726
## X31_jumlah -0.4120 -0.4711
##
## [Refine|G1_RSE_Rendah] ✓ 58 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 21 / 21
## [Refine|G2_RSE_Tinggi] Re-run pada 21 domain (aux dipilih ulang) ...
## [ G2_RSE_Tinggi_ref ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324
##
## [Refine|G2_RSE_Tinggi] ✓ 21 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.503 Min. :26.1 Min. :23.7
## 1st Qu.: 4.10 1st Qu.: 4.180 1st Qu.:39.3 1st Qu.:32.2
## Median : 8.24 Median : 8.051 Median :48.4 Median :37.2
## Mean : 8.65 Mean : 7.979 Mean :52.3 Mean :40.1
## 3rd Qu.:11.12 3rd Qu.:10.558 3rd Qu.:61.0 3rd Qu.:47.1
## Max. :26.70 Max. :17.697 Max. :98.8 Max. :90.6plot_rse(df_s54,"RSE_direct","RSE_eb","Direct","EB-TN-R RSE-NB","S54: EB Beta Top-n/10 RSE-NB — Refined")
S54: EB Beta Top-n/10 RSE-NB + Refine
S55 — Top-n/10 · RSE-ES
## === S55: EB TopN RSE-ES + Refine ===df_s55 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_topn,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 38 / 41
## [Refine|G1_RSE_Bawah] Re-run pada 38 domain (aux dipilih ulang) ...
## [ G1_RSE_Bawah_ref ] y_logit | n = 38 | Vars: X44, X27, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44 -0.4490 -0.5974
## X27 0.4771 0.5723
## X39_jumlah -0.4366 -0.5647
##
## [Refine|G1_RSE_Bawah] ✓ 38 domain diperbarui
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 40 / 40
## [Refine|G2_RSE_Atas] Re-run pada 40 domain (aux dipilih ulang) ...
## [ G2_RSE_Atas_ref ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315
##
## [Refine|G2_RSE_Atas] ✓ 40 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.512 Min. :26.1 Min. :22.4
## 1st Qu.: 4.10 1st Qu.: 4.423 1st Qu.:39.3 1st Qu.:29.7
## Median : 8.24 Median : 8.218 Median :48.4 Median :39.7
## Mean : 8.65 Mean : 7.857 Mean :52.3 Mean :41.8
## 3rd Qu.:11.12 3rd Qu.:10.185 3rd Qu.:61.0 3rd Qu.:51.6
## Max. :26.70 Max. :16.840 Max. :98.8 Max. :93.6plot_rse(df_s55,"RSE_direct","RSE_eb","Direct","EB-TN-R RSE-ES","S55: EB Beta Top-n/10 RSE-ES — Refined")
S55: EB Beta Top-n/10 RSE-ES + Refine
S56 — Top-n/10 · Cluster
## === S56: EB TopN Cluster + Refine ===df_s56 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_tn_logit), aux_fn = seleksi_topn,
model_fn = run_eb_beta, y_col = "y_logit",
rse_col = "RSE_eb", df_orig = df_base
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313
##
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 69 / 70
## [Refine|Cluster_1] Re-run pada 69 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 69 | Vars: X40_mean, X25_mean, X31_mean, X12_mean, X44_mean, X11_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_mean -0.2979 -0.3105
## X25_mean -0.1733 -0.2886
## X31_mean -0.2304 -0.2558
## X12_mean -0.2209 -0.2472
## X44_mean -0.1930 -0.2361
## X11_mean -0.1959 -0.2295
##
## [Refine|Cluster_1] ✓ 69 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 11 / 11
## [Refine|Cluster_2] Re-run pada 11 domain (aux dipilih ulang) ...
## [ Cluster_2_ref ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313
##
## [Refine|Cluster_2] ✓ 11 domain diperbarui## Estimasi y_eb_pct RSE_direct RSE_eb
## Min. : 0.51 Min. : 0.604 Min. :26.1 Min. :24.1
## 1st Qu.: 4.10 1st Qu.: 4.834 1st Qu.:39.3 1st Qu.:34.4
## Median : 8.24 Median : 7.913 Median :48.4 Median :40.0
## Mean : 8.65 Mean : 7.794 Mean :52.3 Mean :40.6
## 3rd Qu.:11.12 3rd Qu.: 9.605 3rd Qu.:61.0 3rd Qu.:45.2
## Max. :26.70 Max. :18.089 Max. :98.8 Max. :80.5plot_rse(df_s56,"RSE_direct","RSE_eb","Direct","EB-TN-R Cluster","S56: EB Beta Top-n/10 Cluster — Refined")
S56: EB Beta Top-n/10 Cluster + Refine
HB Beta Refined
S57–S64 memerlukan JAGS. Setiap run HB melakukan juga CI-check loop (hapus var dengan 2.5%–97.5% menyeberang nol). Re-run domain high-RSE akan memilih variabel baru dari subset tersebut.
S57 — Backward · All
## === S57: HB Backward All + Refine ===df_s57 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_backward,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X40
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X40 -0.3690 -0.3705
##
## [HB Iter 1] Vars aktif (2): X36_mean, X40
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.5303 0.03631 -2.6012 -2.5546 -2.5306 -2.5057 -2.4578
## X36_mean -0.2699 0.03937 -0.3460 -0.2966 -0.2703 -0.2429 -0.1930
## X40 -0.2094 0.03385 -0.2747 -0.2330 -0.2101 -0.1871 -0.1419
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean, X40
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 52 / 81
## [Refine|All] Re-run pada 52 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 52 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.1824 -0.3268
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.9644 0.04523 -3.0532 -2.9945 -2.9636 -2.9336 -2.8770
## X36_mean -0.1956 0.04580 -0.2848 -0.2264 -0.1963 -0.1639 -0.1037
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean
## [Refine|All] ✓ 52 domain diperbarui## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 1.27 Min. :26.1 Min. :16.1
## 1st Qu.: 4.10 1st Qu.: 4.42 1st Qu.:39.3 1st Qu.:22.1
## Median : 8.24 Median : 7.66 Median :48.4 Median :23.6
## Mean : 8.65 Mean : 8.30 Mean :52.3 Mean :25.7
## 3rd Qu.:11.12 3rd Qu.:10.71 3rd Qu.:61.0 3rd Qu.:26.7
## Max. :26.70 Max. :24.97 Max. :98.8 Max. :48.1
S57: HB Backward All + Refine RSE
S57: HB Beta Posterior All Refined
S58 — Backward · RSE-NB
## === S58: HB Backward RSE-NB + Refine ===df_s58 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_backward,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X31_jumlah, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
##
## [HB Iter 1] Vars aktif (3): X3, X31_jumlah, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.2648 0.02958 -2.3236 -2.2849 -2.2646 -2.24501 -2.20726
## X3 0.1724 0.03141 0.1098 0.1515 0.1729 0.19380 0.23320
## X31_jumlah -0.1420 0.03403 -0.2087 -0.1648 -0.1422 -0.11894 -0.07459
## X29_mean -0.1208 0.03317 -0.1848 -0.1436 -0.1204 -0.09859 -0.05673
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X3, X31_jumlah, X29_mean
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.4287 -0.6631
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.5093 0.06017 -3.6275 -3.5508 -3.5084 -3.4675 -3.3900
## X36_mean -0.5082 0.06340 -0.6311 -0.5514 -0.5084 -0.4647 -0.3852
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 4 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 4 domain (aux dipilih ulang) ...## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 0 / 21
## [Refine|G2_RSE_Tinggi] < 4 domain → skip, pakai hasil awal## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.514 Min. :26.1 Min. : 5.61
## 1st Qu.: 4.10 1st Qu.: 4.289 1st Qu.:39.3 1st Qu.:12.59
## Median : 8.24 Median : 8.282 Median :48.4 Median :14.91
## Mean : 8.65 Mean : 8.663 Mean :52.3 Mean :15.53
## 3rd Qu.:11.12 3rd Qu.:10.971 3rd Qu.:61.0 3rd Qu.:16.50
## Max. :26.70 Max. :25.013 Max. :98.8 Max. :49.39
S58: HB Backward RSE-NB + Refine RSE
S58: HB Beta Posterior RSE-NB Refined
S59 — Backward · RSE-ES
## === S59: HB Backward RSE-ES + Refine ===df_s59 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_backward,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X8, X39_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
##
## [HB Iter 1] Vars aktif (2): X8, X39_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.1532 0.03127 -2.2171 -2.1735 -2.1527 -2.1322 -2.0942
## X8 0.2000 0.03036 0.1403 0.1793 0.2005 0.2210 0.2581
## X39_jumlah -0.1875 0.03148 -0.2501 -0.2088 -0.1878 -0.1661 -0.1263
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X8, X39_jumlah
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3763 -0.5547
##
## [HB Iter 1] Vars aktif (1): X31_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0428 0.05493 -3.152 -3.0807 -3.0417 -3.0051 -2.9385
## X31_mean -0.5112 0.05829 -0.625 -0.5502 -0.5122 -0.4717 -0.3969
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X31_mean
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 0 / 41
## [Refine|G1_RSE_Bawah] < 4 domain → skip, pakai hasil awal
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 4 / 40
## [Refine|G2_RSE_Atas] Re-run pada 4 domain (aux dipilih ulang) ...## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.517 Min. :26.1 Min. : 2.75
## 1st Qu.: 4.10 1st Qu.: 4.122 1st Qu.:39.3 1st Qu.: 3.88
## Median : 8.24 Median : 8.220 Median :48.4 Median : 6.87
## Mean : 8.65 Mean : 8.647 Mean :52.3 Mean : 9.07
## 3rd Qu.:11.12 3rd Qu.:11.130 3rd Qu.:61.0 3rd Qu.:12.69
## Max. :26.70 Max. :26.415 Max. :98.8 Max. :31.97
S59: HB Backward RSE-ES + Refine RSE
S59: HB Beta Posterior RSE-ES Refined
S60 — Backward · Cluster
## === S60: HB Backward Cluster + Refine ===df_s60 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_bk_logit), aux_fn = seleksi_backward,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 60 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3236 -0.4251
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.7517 0.04742 -2.8456 -2.7831 -2.7515 -2.7196 -2.6611
## X36_mean -0.3425 0.04706 -0.4359 -0.3743 -0.3429 -0.3113 -0.2509
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean
## [ Cluster_2 ] y_logit | n = 21 | Vars: X40_jumlah
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40_jumlah -0.3419 -0.3562
##
## [HB Iter 1] Vars aktif (1): X40_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.0344 0.04982 -2.1328 -2.0679 -2.0338 -2.0010 -1.93689
## X40_jumlah -0.1766 0.05146 -0.2768 -0.2113 -0.1762 -0.1419 -0.07777
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X40_jumlah
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 24 / 60
## [Refine|Cluster_1] Re-run pada 24 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 24 | Vars: X43_mean, X33_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X43_mean 0.4972 0.4229
## X33_mean -0.3596 -0.3738
##
## [HB Iter 1] Vars aktif (2): X43_mean, X33_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.6766 0.04812 -3.7728 -3.7093 -3.6750 -3.6448 -3.5847
## X43_mean 0.2115 0.04697 0.1181 0.1799 0.2118 0.2425 0.3028
## X33_mean -0.1732 0.04967 -0.2715 -0.2059 -0.1735 -0.1404 -0.0737
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X43_mean, X33_mean
## [Refine|Cluster_1] ✓ 24 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 0 / 21
## [Refine|Cluster_2] < 4 domain → skip, pakai hasil awal## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.941 Min. :26.1 Min. : 8.82
## 1st Qu.: 4.10 1st Qu.: 4.004 1st Qu.:39.3 1st Qu.:14.24
## Median : 8.24 Median : 8.011 Median :48.4 Median :17.64
## Mean : 8.65 Mean : 8.483 Mean :52.3 Mean :17.68
## 3rd Qu.:11.12 3rd Qu.:10.909 3rd Qu.:61.0 3rd Qu.:20.88
## Max. :26.70 Max. :25.845 Max. :98.8 Max. :40.59plot_rse(df_s60,"RSE","RSE_hb","Direct","HB-BK-R Cluster","S60: HB Beta Backward Cluster — Refined")
S60: HB Backward Cluster + Refine RSE
plot_hb_posterior(df_s60, "S60: HB Beta Backward Cluster — Refined", group_col = "cluster_bk_logit")
S60: HB Beta Posterior Cluster Refined
S61 — Top-n/10 · All
## === S61: HB TopN All + Refine ===df_s61 <- run_with_refine(
seg_list = list("All" = df_base), aux_fn = seleksi_topn,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ All ] y_logit | n = 81 | Vars: X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3220 -0.4311
## X30_mean -0.3041 -0.3975
## X31_mean -0.3188 -0.3953
## X29_mean -0.3220 -0.3921
## X8 0.3439 0.3797
## X40 -0.3690 -0.3705
## X45 -0.3224 -0.3604
##
## [HB Iter 1] Vars aktif (7): X36_mean, X30_mean, X31_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.52086 0.03603 -2.59070 -2.54531 -2.52057 -2.49726 -2.44832
## X36_mean -0.03394 0.07773 -0.18760 -0.08680 -0.03411 0.01994 0.11924
## X30_mean -0.05760 0.05945 -0.16997 -0.09853 -0.05787 -0.01751 0.05904
## X31_mean -0.03315 0.07428 -0.18293 -0.08231 -0.03226 0.01765 0.10868
## X29_mean -0.12430 0.06712 -0.25771 -0.16863 -0.12321 -0.08024 0.00570
## X8 0.07834 0.05785 -0.03760 0.03964 0.07843 0.11776 0.19076
## X40 -0.20224 0.03530 -0.27020 -0.22615 -0.20230 -0.17785 -0.13389
## X45 0.02432 0.05725 -0.09039 -0.01427 0.02504 0.06315 0.13705
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X36_mean, X30_mean, X31_mean, X29_mean, X8, X45]
## [HB] → Keluarkan 'X36_mean' (SD = 0.07773, terbesar)
##
## [HB Iter 2] Vars aktif (6): X30_mean, X31_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.53143 0.03617 -2.60391 -2.55549 -2.53129 -2.50684 -2.46202
## X30_mean -0.07204 0.05442 -0.17843 -0.10885 -0.07245 -0.03517 0.03557
## X31_mean -0.03954 0.06611 -0.16885 -0.08340 -0.03988 0.00438 0.08871
## X29_mean -0.13435 0.05978 -0.24685 -0.17782 -0.13402 -0.09364 -0.01707
## X8 0.08597 0.05348 -0.01885 0.04925 0.08734 0.12333 0.18914
## X40 -0.20617 0.03555 -0.27727 -0.23032 -0.20605 -0.18171 -0.13658
## X45 0.02640 0.05597 -0.08047 -0.01201 0.02533 0.06467 0.13688
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X30_mean, X31_mean, X8, X45]
## [HB] → Keluarkan 'X31_mean' (SD = 0.06611, terbesar)
##
## [HB Iter 3] Vars aktif (5): X30_mean, X29_mean, X8, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51803 0.03611 -2.58811 -2.54255 -2.51794 -2.49341 -2.44678
## X30_mean -0.07316 0.04992 -0.17337 -0.10703 -0.07211 -0.04005 0.02331
## X29_mean -0.15256 0.05242 -0.25587 -0.18727 -0.15218 -0.11748 -0.05003
## X8 0.08885 0.05209 -0.01348 0.05357 0.08844 0.12489 0.19002
## X40 -0.20321 0.03535 -0.27205 -0.22742 -0.20321 -0.17945 -0.13342
## X45 0.01497 0.05125 -0.08697 -0.01967 0.01557 0.04955 0.11314
##
## [HB] ✗ Iter 3 — CI menyeberang nol: [X30_mean, X8, X45]
## [HB] → Keluarkan 'X8' (SD = 0.05209, terbesar)
##
## [HB Iter 4] Vars aktif (4): X30_mean, X29_mean, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51738 0.03538 -2.5869 -2.5407 -2.51701 -2.49333 -2.44780
## X30_mean -0.07768 0.04877 -0.1744 -0.1104 -0.07726 -0.04562 0.01741
## X29_mean -0.18611 0.04859 -0.2817 -0.2185 -0.18623 -0.15352 -0.09174
## X40 -0.21140 0.03441 -0.2760 -0.2349 -0.21185 -0.18849 -0.14141
## X45 -0.02148 0.04616 -0.1118 -0.0528 -0.02188 0.00971 0.06744
##
## [HB] ✗ Iter 4 — CI menyeberang nol: [X30_mean, X45]
## [HB] → Keluarkan 'X30_mean' (SD = 0.04877, terbesar)
##
## [HB Iter 5] Vars aktif (3): X29_mean, X40, X45
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.51625 0.03513 -2.5865 -2.53987 -2.51600 -2.49271 -2.44743
## X29_mean -0.22527 0.04243 -0.3066 -0.25417 -0.22572 -0.19672 -0.13970
## X40 -0.21375 0.03495 -0.2825 -0.23731 -0.21345 -0.18945 -0.14647
## X45 -0.04969 0.04181 -0.1306 -0.07771 -0.05074 -0.02117 0.03277
##
## [HB] ✗ Iter 5 — CI menyeberang nol: [X45]
## [HB] → Keluarkan 'X45' (SD = 0.04181, terbesar)
##
## [HB Iter 6] Vars aktif (2): X29_mean, X40
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.5115 0.03618 -2.5804 -2.5356 -2.5112 -2.4873 -2.4397
## X29_mean -0.2569 0.03707 -0.3286 -0.2818 -0.2571 -0.2317 -0.1835
## X40 -0.2255 0.03274 -0.2901 -0.2474 -0.2254 -0.2035 -0.1612
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 6.
## [HB] ✓ Variabel final: X29_mean, X40
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|All] Domain RSE > 25%: 69 / 81
## [Refine|All] Re-run pada 69 domain (aux dipilih ulang) ...
## [ All_ref ] y_logit | n = 69 | Vars: X36_mean, X44_mean, X30_mean, X31_mean, X8, X29_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.2706 -0.3904
## X44_mean -0.2870 -0.3471
## X30_mean -0.2315 -0.3445
## X31_mean -0.2567 -0.3421
## X8 0.3062 0.3326
## X29_mean -0.2500 -0.3311
##
## [HB Iter 1] Vars aktif (6): X36_mean, X44_mean, X30_mean, X31_mean, X8, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.69187 0.04223 -2.77430 -2.72087 -2.69139 -2.66298 -2.61097
## X36_mean -0.12483 0.09189 -0.30310 -0.18647 -0.12602 -0.06241 0.05793
## X44_mean 0.00522 0.08988 -0.16886 -0.05620 0.00526 0.06732 0.17948
## X30_mean -0.05638 0.06365 -0.17678 -0.10029 -0.05574 -0.01375 0.06841
## X31_mean -0.04347 0.07652 -0.19490 -0.09361 -0.04360 0.00750 0.10679
## X8 0.11661 0.07153 -0.02102 0.06949 0.11627 0.16445 0.25792
## X29_mean 0.02021 0.07904 -0.13543 -0.03230 0.02056 0.07274 0.17752
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X36_mean, X44_mean, X30_mean, X31_mean, X8, X29_mean]
## [HB] → Keluarkan 'X36_mean' (SD = 0.09189, terbesar)
##
## [HB Iter 2] Vars aktif (5): X44_mean, X30_mean, X31_mean, X8, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.67882 0.04177 -2.75987 -2.70716 -2.67808 -2.65146 -2.59782
## X44_mean -0.04118 0.07610 -0.18928 -0.09401 -0.04075 0.01044 0.10705
## X30_mean -0.07626 0.06103 -0.19548 -0.11777 -0.07595 -0.03575 0.04344
## X31_mean -0.06086 0.07404 -0.20638 -0.10978 -0.06026 -0.01066 0.08506
## X8 0.12778 0.06614 -0.00527 0.08309 0.12835 0.17246 0.25676
## X29_mean -0.00748 0.07209 -0.15108 -0.05427 -0.00858 0.03913 0.13522
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X44_mean, X30_mean, X31_mean, X8, X29_mean]
## [HB] → Keluarkan 'X44_mean' (SD = 0.07610, terbesar)
##
## [HB Iter 3] Vars aktif (4): X30_mean, X31_mean, X8, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.67046 0.03926 -2.75005 -2.69660 -2.67012 -2.64430 -2.59486
## X30_mean -0.07102 0.05676 -0.18209 -0.10890 -0.07088 -0.03288 0.03912
## X31_mean -0.06517 0.06620 -0.19705 -0.11041 -0.06384 -0.01983 0.06169
## X8 0.15230 0.05277 0.04915 0.11707 0.15269 0.18824 0.25326
## X29_mean -0.02202 0.06350 -0.14466 -0.06357 -0.02237 0.01962 0.10144
##
## [HB] ✗ Iter 3 — CI menyeberang nol: [X30_mean, X31_mean, X29_mean]
## [HB] → Keluarkan 'X31_mean' (SD = 0.06620, terbesar)
##
## [HB Iter 4] Vars aktif (3): X30_mean, X8, X29_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.67061 0.04131 -2.75169 -2.6990 -2.67030 -2.64329 -2.59000
## X30_mean -0.08535 0.05299 -0.19049 -0.1205 -0.08528 -0.04977 0.01727
## X8 0.16353 0.05070 0.06532 0.1294 0.16293 0.19813 0.26344
## X29_mean -0.05569 0.05377 -0.15935 -0.0918 -0.05595 -0.01896 0.04842
##
## [HB] ✗ Iter 4 — CI menyeberang nol: [X30_mean, X29_mean]
## [HB] → Keluarkan 'X29_mean' (SD = 0.05377, terbesar)
##
## [HB Iter 5] Vars aktif (2): X30_mean, X8
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.6723 0.04091 -2.75468 -2.699 -2.6723 -2.64531 -2.59323
## X30_mean -0.1214 0.04664 -0.21293 -0.153 -0.1211 -0.09005 -0.03099
## X8 0.1795 0.04563 0.09183 0.148 0.1790 0.21095 0.26992
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 5.
## [HB] ✓ Variabel final: X30_mean, X8
## [Refine|All] ✓ 69 domain diperbarui## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 1.66 Min. :26.1 Min. :17.4
## 1st Qu.: 4.10 1st Qu.: 4.57 1st Qu.:39.3 1st Qu.:24.5
## Median : 8.24 Median : 8.14 Median :48.4 Median :27.1
## Mean : 8.65 Mean : 8.52 Mean :52.3 Mean :29.0
## 3rd Qu.:11.12 3rd Qu.:10.87 3rd Qu.:61.0 3rd Qu.:32.3
## Max. :26.70 Max. :24.76 Max. :98.8 Max. :48.9
S61: HB Top-n/10 All + Refine RSE
S61: HB Beta Posterior All Refined
S62 — Top-n/10 · RSE-NB
## === S62: HB TopN RSE-NB + Refine ===df_s62 <- run_with_refine(
seg_list = split(df_base, df_base$grup_jenks), aux_fn = seleksi_topn,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ G1_RSE_Rendah ] y_logit | n = 60 | Vars: X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X3 0.4821 0.4882
## X27 0.4405 0.4880
## X31_jumlah -0.4258 -0.4828
## X29_mean -0.3784 -0.4642
## X4 0.4425 0.4608
## X28_mean -0.3437 -0.4251
##
## [HB Iter 1] Vars aktif (6): X3, X27, X31_jumlah, X29_mean, X4, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.26329 0.02909 -2.31999 -2.28290 -2.26349 -2.24370 -2.20663
## X3 0.17086 0.04404 0.08330 0.14146 0.17107 0.20021 0.25660
## X27 0.04713 0.04366 -0.03898 0.01750 0.04757 0.07592 0.13427
## X31_jumlah -0.16663 0.03700 -0.23780 -0.19141 -0.16729 -0.14199 -0.09229
## X29_mean -0.24692 0.05407 -0.35267 -0.28277 -0.24727 -0.21019 -0.14036
## X4 -0.04988 0.04328 -0.13374 -0.07891 -0.05023 -0.02031 0.03326
## X28_mean 0.14472 0.05204 0.04072 0.10944 0.14529 0.18050 0.24258
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X27, X4]
## [HB] → Keluarkan 'X27' (SD = 0.04366, terbesar)
##
## [HB Iter 2] Vars aktif (5): X3, X31_jumlah, X29_mean, X4, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.26048 0.02960 -2.31708 -2.28021 -2.26082 -2.24022 -2.20175
## X3 0.20406 0.03868 0.12916 0.17835 0.20437 0.22969 0.27910
## X31_jumlah -0.16525 0.03819 -0.24046 -0.19030 -0.16517 -0.13956 -0.09130
## X29_mean -0.27945 0.05407 -0.38528 -0.31584 -0.27947 -0.24247 -0.17510
## X4 -0.05243 0.04084 -0.13225 -0.07994 -0.05271 -0.02509 0.02757
## X28_mean 0.16013 0.05290 0.05604 0.12553 0.15983 0.19527 0.26198
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X4]
## [HB] → Keluarkan 'X4' (SD = 0.04084, terbesar)
##
## [HB Iter 3] Vars aktif (4): X3, X31_jumlah, X29_mean, X28_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.2616 0.02819 -2.31622 -2.2805 -2.2614 -2.2428 -2.20612
## X3 0.1696 0.03147 0.10795 0.1485 0.1700 0.1915 0.22980
## X31_jumlah -0.1527 0.03621 -0.22271 -0.1772 -0.1532 -0.1281 -0.07935
## X29_mean -0.2497 0.05220 -0.35234 -0.2851 -0.2494 -0.2138 -0.14802
## X28_mean 0.1375 0.04874 0.04298 0.1046 0.1370 0.1706 0.23247
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 3.
## [HB] ✓ Variabel final: X3, X31_jumlah, X29_mean, X28_mean
## [ G2_RSE_Tinggi ] y_logit | n = 21 | Vars: X31_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X31_mean -0.3973 -0.6324
##
## [HB Iter 1] Vars aktif (1): X31_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.509 0.06046 -3.6248 -3.5493 -3.5102 -3.4698 -3.3874
## X31_mean -0.485 0.06442 -0.6101 -0.5289 -0.4853 -0.4407 -0.3604
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X31_mean
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Rendah] Domain RSE > 25%: 5 / 60
## [Refine|G1_RSE_Rendah] Re-run pada 5 domain (aux dipilih ulang) ...
## [ G1_RSE_Rendah_ref ] y_logit | n = 5 | Vars: X43
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X43 0.8217 0.8607
##
## [HB Iter 1] Vars aktif (1): X43
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.6361 0.1456 -3.9103 -3.7381 -3.6379 -3.5403 -3.3368
## X43 0.5263 0.1519 0.2166 0.4272 0.5327 0.6298 0.8112
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X43
## [Refine|G1_RSE_Rendah] ✓ 5 domain diperbarui
## [Refine|G2_RSE_Tinggi] Domain RSE > 25%: 0 / 21
## [Refine|G2_RSE_Tinggi] < 4 domain → skip, pakai hasil awal## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.519 Min. :26.1 Min. : 5.44
## 1st Qu.: 4.10 1st Qu.: 4.119 1st Qu.:39.3 1st Qu.:13.36
## Median : 8.24 Median : 8.433 Median :48.4 Median :16.48
## Mean : 8.65 Mean : 8.573 Mean :52.3 Mean :15.79
## 3rd Qu.:11.12 3rd Qu.:11.102 3rd Qu.:61.0 3rd Qu.:18.55
## Max. :26.70 Max. :24.208 Max. :98.8 Max. :25.67
S62: HB Top-n/10 RSE-NB + Refine RSE
S62: HB Beta Posterior RSE-NB Refined
S63 — Top-n/10 · RSE-ES
## === S63: HB TopN RSE-ES + Refine ===df_s63 <- run_with_refine(
seg_list = split(df_base, df_base$grup_equal), aux_fn = seleksi_topn,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ G1_RSE_Bawah ] y_logit | n = 41 | Vars: X44_mean, X8, X39_jumlah, X27
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X44_mean -0.4399 -0.5825
## X8 0.4454 0.5695
## X39_jumlah -0.4373 -0.5570
## X27 0.4608 0.5500
##
## [HB Iter 1] Vars aktif (4): X44_mean, X8, X39_jumlah, X27
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.15909 0.02871 -2.21533 -2.17833 -2.15924 -2.13991 -2.1026
## X44_mean -0.02522 0.04869 -0.12133 -0.05893 -0.02485 0.01072 0.0646
## X8 0.12573 0.04327 0.04454 0.09573 0.12558 0.15592 0.2090
## X39_jumlah -0.16478 0.03217 -0.22997 -0.18602 -0.16342 -0.14259 -0.1034
## X27 0.09056 0.03520 0.02294 0.06687 0.08997 0.11399 0.1617
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X44_mean]
## [HB] → Keluarkan 'X44_mean' (SD = 0.04869, terbesar)
##
## [HB Iter 2] Vars aktif (3): X8, X39_jumlah, X27
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.15593 0.02769 -2.20940 -2.17496 -2.15688 -2.1374 -2.0992
## X8 0.13896 0.03518 0.06892 0.11451 0.13973 0.1636 0.2049
## X39_jumlah -0.17079 0.03334 -0.23518 -0.19445 -0.17096 -0.1483 -0.1050
## X27 0.09676 0.03565 0.02828 0.07256 0.09625 0.1211 0.1673
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 2.
## [HB] ✓ Variabel final: X8, X39_jumlah, X27
## [ G2_RSE_Atas ] y_logit | n = 40 | Vars: X27_mean, X28_mean, X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X27_mean -0.4069 -0.5399
## X28_mean -0.4137 -0.5322
## X36_mean -0.3461 -0.5315
##
## [HB Iter 1] Vars aktif (3): X27_mean, X28_mean, X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0697 0.05422 -3.1745 -3.1077 -3.0703 -3.03274 -2.96241
## X27_mean -0.2138 0.08502 -0.3804 -0.2700 -0.2166 -0.15841 -0.04124
## X28_mean -0.1123 0.08213 -0.2756 -0.1668 -0.1109 -0.05766 0.04929
## X36_mean -0.2259 0.08458 -0.3956 -0.2819 -0.2281 -0.17053 -0.05721
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X28_mean]
## [HB] → Keluarkan 'X28_mean' (SD = 0.08213, terbesar)
##
## [HB Iter 2] Vars aktif (2): X27_mean, X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.0475 0.05381 -3.1519 -3.0844 -3.0467 -3.0112 -2.9428
## X27_mean -0.2905 0.06954 -0.4281 -0.3372 -0.2904 -0.2429 -0.1565
## X36_mean -0.2610 0.06894 -0.3972 -0.3054 -0.2612 -0.2153 -0.1280
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 2.
## [HB] ✓ Variabel final: X27_mean, X36_mean
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|G1_RSE_Bawah] Domain RSE > 25%: 0 / 41
## [Refine|G1_RSE_Bawah] < 4 domain → skip, pakai hasil awal
## [Refine|G2_RSE_Atas] Domain RSE > 25%: 3 / 40
## [Refine|G2_RSE_Atas] < 4 domain → skip, pakai hasil awal## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.544 Min. :26.1 Min. : 5.01
## 1st Qu.: 4.10 1st Qu.: 4.141 1st Qu.:39.3 1st Qu.: 6.84
## Median : 8.24 Median : 8.209 Median :48.4 Median : 8.61
## Mean : 8.65 Mean : 8.647 Mean :52.3 Mean :10.54
## 3rd Qu.:11.12 3rd Qu.:11.143 3rd Qu.:61.0 3rd Qu.:12.49
## Max. :26.70 Max. :26.443 Max. :98.8 Max. :31.40
S63: HB Top-n/10 RSE-ES + Refine RSE
S63: HB Beta Posterior RSE-ES Refined
S64 — Top-n/10 · Cluster
## === S64: HB TopN Cluster + Refine ===df_s64 <- run_with_refine(
seg_list = split(df_base, df_base$cluster_tn_logit), aux_fn = seleksi_topn,
model_fn = run_hbbeta, y_col = "y_logit",
rse_col = "RSE_hb", df_orig = df_base, min_n = 5
)## >>> Tahap 1: model utama
## [ Cluster_1 ] y_logit | n = 70 | Vars: X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X40 -0.3091 -0.3196
## X25_jumlah -0.2038 -0.2746
## X33_mean -0.2294 -0.2340
## X31_mean -0.1967 -0.2304
## X48 -0.2223 -0.2146
## X35_jumlah -0.2500 -0.2117
## X12_mean -0.1739 -0.2109
##
## [HB Iter 1] Vars aktif (7): X40, X25_jumlah, X33_mean, X31_mean, X48, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41550 0.03837 -2.48810 -2.44245 -2.41547 -2.38931 -2.34054
## X40 -0.16230 0.03813 -0.23613 -0.18873 -0.16262 -0.13554 -0.08937
## X25_jumlah -0.08516 0.04158 -0.16658 -0.11356 -0.08566 -0.05612 -0.00425
## X33_mean -0.10986 0.04003 -0.18889 -0.13669 -0.10974 -0.08361 -0.02971
## X31_mean 0.05369 0.04969 -0.04501 0.02022 0.05369 0.08810 0.14863
## X48 -0.03763 0.04390 -0.12448 -0.06676 -0.03808 -0.00845 0.04829
## X35_jumlah -0.09226 0.03992 -0.17282 -0.11871 -0.09282 -0.06529 -0.01337
## X12_mean -0.08932 0.04610 -0.17796 -0.11983 -0.08986 -0.05828 0.00197
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X31_mean, X48, X12_mean]
## [HB] → Keluarkan 'X31_mean' (SD = 0.04969, terbesar)
##
## [HB Iter 2] Vars aktif (6): X40, X25_jumlah, X33_mean, X48, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41980 0.03830 -2.49566 -2.44568 -2.41982 -2.39380 -2.34585
## X40 -0.16942 0.03822 -0.24415 -0.19496 -0.16947 -0.14481 -0.09331
## X25_jumlah -0.07603 0.03895 -0.15120 -0.10224 -0.07684 -0.04978 0.00258
## X33_mean -0.10647 0.03913 -0.18291 -0.13295 -0.10579 -0.07924 -0.03172
## X48 -0.01870 0.04086 -0.09921 -0.04611 -0.01855 0.00827 0.06272
## X35_jumlah -0.08286 0.03941 -0.16018 -0.10872 -0.08318 -0.05705 -0.00545
## X12_mean -0.06330 0.04025 -0.14064 -0.09092 -0.06321 -0.03674 0.01569
##
## [HB] ✗ Iter 2 — CI menyeberang nol: [X25_jumlah, X48, X12_mean]
## [HB] → Keluarkan 'X48' (SD = 0.04086, terbesar)
##
## [HB Iter 3] Vars aktif (5): X40, X25_jumlah, X33_mean, X35_jumlah, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41679 0.03798 -2.4892 -2.44331 -2.41681 -2.39115 -2.34257
## X40 -0.17695 0.03595 -0.2471 -0.20124 -0.17709 -0.15291 -0.10543
## X25_jumlah -0.07801 0.03923 -0.1532 -0.10495 -0.07789 -0.05132 -0.00122
## X33_mean -0.10967 0.03921 -0.1871 -0.13638 -0.10943 -0.08343 -0.03333
## X35_jumlah -0.08244 0.03962 -0.1578 -0.10930 -0.08215 -0.05598 -0.00453
## X12_mean -0.06770 0.03795 -0.1435 -0.09293 -0.06803 -0.04256 0.00830
##
## [HB] ✗ Iter 3 — CI menyeberang nol: [X12_mean]
## [HB] → Keluarkan 'X12_mean' (SD = 0.03795, terbesar)
##
## [HB Iter 4] Vars aktif (4): X40, X25_jumlah, X33_mean, X35_jumlah
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.41544 0.03879 -2.4920 -2.4418 -2.41503 -2.38945 -2.33925
## X40 -0.17219 0.03627 -0.2450 -0.1967 -0.17222 -0.14747 -0.10130
## X25_jumlah -0.09541 0.03865 -0.1716 -0.1221 -0.09494 -0.06954 -0.02069
## X33_mean -0.12129 0.03856 -0.1960 -0.1480 -0.12224 -0.09581 -0.04437
## X35_jumlah -0.09018 0.03909 -0.1685 -0.1162 -0.08994 -0.06373 -0.01535
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 4.
## [HB] ✓ Variabel final: X40, X25_jumlah, X33_mean, X35_jumlah
## [ Cluster_2 ] y_logit | n = 11 | Vars: X36_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X36_mean -0.3819 -0.6313
##
## [HB Iter 1] Vars aktif (1): X36_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -3.2833 0.1056 -3.4891 -3.3540 -3.283 -3.213 -3.0710
## X36_mean -0.4585 0.1098 -0.6738 -0.5316 -0.460 -0.385 -0.2457
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 1.
## [HB] ✓ Variabel final: X36_mean
## >>> Tahap 2: refinement domain RSE > 25 %
## [Refine|Cluster_1] Domain RSE > 25%: 50 / 70
## [Refine|Cluster_1] Re-run pada 50 domain (aux dipilih ulang) ...
## [ Cluster_1_ref ] y_logit | n = 50 | Vars: X21, X26_jumlah, X2_mean, X25_mean, X12_mean
## Variabel cor(Direct) cor(Logit)
## ------------------------------------------
## X21 0.3539 0.3351
## X26_jumlah -0.2369 -0.2743
## X2_mean 0.3109 0.2653
## X25_mean -0.0658 -0.2323
## X12_mean -0.1729 -0.2267
##
## [HB Iter 1] Vars aktif (5): X21, X26_jumlah, X2_mean, X25_mean, X12_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.74122 0.03960 -2.81937 -2.7677 -2.74179 -2.71544 -2.66158
## X21 0.15643 0.03750 0.08251 0.1314 0.15604 0.18235 0.22927
## X26_jumlah -0.20215 0.04213 -0.28732 -0.2303 -0.20177 -0.17351 -0.12178
## X2_mean 0.22466 0.03987 0.14717 0.1977 0.22521 0.25117 0.30252
## X25_mean -0.11480 0.04326 -0.19917 -0.1442 -0.11442 -0.08557 -0.03035
## X12_mean -0.05013 0.04318 -0.13522 -0.0796 -0.05077 -0.02109 0.03453
##
## [HB] ✗ Iter 1 — CI menyeberang nol: [X12_mean]
## [HB] → Keluarkan 'X12_mean' (SD = 0.04318, terbesar)
##
## [HB Iter 2] Vars aktif (4): X21, X26_jumlah, X2_mean, X25_mean
## Mean SD 2.5% 25% 50% 75% 97.5%
## intercept -2.7484 0.03943 -2.82794 -2.7742 -2.7479 -2.7214 -2.67191
## X21 0.1649 0.03695 0.09445 0.1405 0.1645 0.1895 0.23965
## X26_jumlah -0.2118 0.04198 -0.29682 -0.2399 -0.2120 -0.1836 -0.12807
## X2_mean 0.2210 0.03930 0.14343 0.1943 0.2210 0.2475 0.29605
## X25_mean -0.1432 0.04137 -0.22611 -0.1707 -0.1426 -0.1156 -0.06284
##
## [HB] ✓ Semua CI aman — selesai pada iterasi 2.
## [HB] ✓ Variabel final: X21, X26_jumlah, X2_mean, X25_mean
## [Refine|Cluster_1] ✓ 50 domain diperbarui
## [Refine|Cluster_2] Domain RSE > 25%: 3 / 11
## [Refine|Cluster_2] < 4 domain → skip, pakai hasil awal## Estimasi y_hb_pct RSE RSE_hb
## Min. : 0.51 Min. : 0.782 Min. :26.1 Min. :17.1
## 1st Qu.: 4.10 1st Qu.: 4.456 1st Qu.:39.3 1st Qu.:21.0
## Median : 8.24 Median : 7.755 Median :48.4 Median :23.5
## Mean : 8.65 Mean : 8.329 Mean :52.3 Mean :24.6
## 3rd Qu.:11.12 3rd Qu.:11.062 3rd Qu.:61.0 3rd Qu.:26.0
## Max. :26.70 Max. :24.867 Max. :98.8 Max. :46.9plot_rse(df_s64,"RSE","RSE_hb","Direct","HB-TN-R Cluster","S64: HB Beta Top-n/10 Cluster — Refined")
S64: HB Top-n/10 Cluster + Refine RSE
plot_hb_posterior(df_s64, "S64: HB Beta Top-n/10 Cluster — Refined", group_col = "cluster_tn_logit")
S64: HB Beta Posterior Cluster Refined
10. Evaluasi Komparatif — 32 Skenario
Kumpulkan Hasil
# Standarisasi: tiap df hasil ditambah y_model & RSE_model
standardize_result <- function(df, y_col, rse_col, skenario_id) {
if (is.null(df)) return(NULL)
df %>%
mutate(
y_model = .data[[y_col]],
RSE_model = .data[[rse_col]],
Skenario = skenario_id
) %>%
dplyr::select(Kako, Estimasi, RSE, y_model, RSE_model, Skenario)
}
# Baseline direct
df_direct_eval <- df_base %>%
mutate(RSE_direct_eval = RSE_direct_col) %>%
dplyr::select(Kako, Estimasi, RSE, RSE_direct_eval)
# ── Tabel distribusi RSE ─────────────────────────────────────────────────
make_tabel_row <- function(df, rse_col, label) {
if (is.null(df) || !rse_col %in% names(df)) {
return(data.frame(Skenario = label, `RSE>=25` = NA, `RSE<25` = NA,
`%<25` = NA, `%<15` = NA, `CV Mean` = NA,
check.names = FALSE))
}
tabel_rse(df, rse_col, label)
}
eval_rows <- list(
make_tabel_row(df_base, "RSE_direct_col", "0. Direct"),
# EBLUP
make_tabel_row(df_s1, "RSE_eblup", "S01 EBLUP-BK All"),
make_tabel_row(df_s2, "RSE_eblup", "S02 EBLUP-BK RSE-NB"),
make_tabel_row(df_s3, "RSE_eblup", "S03 EBLUP-BK RSE-ES"),
make_tabel_row(df_s4, "RSE_eblup", "S04 EBLUP-BK Cluster"),
make_tabel_row(df_s5, "RSE_eblup", "S05 EBLUP-TN All"),
make_tabel_row(df_s6, "RSE_eblup", "S06 EBLUP-TN RSE-NB"),
make_tabel_row(df_s7, "RSE_eblup", "S07 EBLUP-TN RSE-ES"),
make_tabel_row(df_s8, "RSE_eblup", "S08 EBLUP-TN Cluster"),
# GLMM
make_tabel_row(df_s9, "RSE_glmm", "S09 GLMM-BK All"),
make_tabel_row(df_s10, "RSE_glmm", "S10 GLMM-BK RSE-NB"),
make_tabel_row(df_s11, "RSE_glmm", "S11 GLMM-BK RSE-ES"),
make_tabel_row(df_s12, "RSE_glmm", "S12 GLMM-BK Cluster"),
make_tabel_row(df_s13, "RSE_glmm", "S13 GLMM-TN All"),
make_tabel_row(df_s14, "RSE_glmm", "S14 GLMM-TN RSE-NB"),
make_tabel_row(df_s15, "RSE_glmm", "S15 GLMM-TN RSE-ES"),
make_tabel_row(df_s16, "RSE_glmm", "S16 GLMM-TN Cluster"),
# EB Beta
make_tabel_row(df_s17, "RSE_eb", "S17 EB-BK All"),
make_tabel_row(df_s18, "RSE_eb", "S18 EB-BK RSE-NB"),
make_tabel_row(df_s19, "RSE_eb", "S19 EB-BK RSE-ES"),
make_tabel_row(df_s20, "RSE_eb", "S20 EB-BK Cluster"),
make_tabel_row(df_s21, "RSE_eb", "S21 EB-TN All"),
make_tabel_row(df_s22, "RSE_eb", "S22 EB-TN RSE-NB"),
make_tabel_row(df_s23, "RSE_eb", "S23 EB-TN RSE-ES"),
make_tabel_row(df_s24, "RSE_eb", "S24 EB-TN Cluster"),
# HB Beta
make_tabel_row(df_s25, "RSE_hb", "S25 HB-BK All"),
make_tabel_row(df_s26, "RSE_hb", "S26 HB-BK RSE-NB"),
make_tabel_row(df_s27, "RSE_hb", "S27 HB-BK RSE-ES"),
make_tabel_row(df_s28, "RSE_hb", "S28 HB-BK Cluster"),
make_tabel_row(df_s29, "RSE_hb", "S29 HB-TN All"),
make_tabel_row(df_s30, "RSE_hb", "S30 HB-TN RSE-NB"),
make_tabel_row(df_s31, "RSE_hb", "S31 HB-TN RSE-ES"),
make_tabel_row(df_s32, "RSE_hb", "S32 HB-TN Cluster")
)
tabel_eval <- bind_rows(eval_rows)Tabel Distribusi RSE
kable(tabel_eval,
caption = "Distribusi RSE — 32 Skenario + Direct",
col.names = c("Skenario", "RSE≥25", "RSE<25", "% <25%", "% <15%", "CV Mean")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"),
full_width = FALSE, font_size = 12) %>%
row_spec(1, background = "#fff3cd") %>%
row_spec(
which(tabel_eval$`%<25` == max(tabel_eval$`%<25`, na.rm = TRUE))[-1],
bold = TRUE, background = "#d4edda"
) %>%
pack_rows("EBLUP", 2, 9) %>%
pack_rows("GLMM", 10, 17) %>%
pack_rows("EB Beta-Binomial", 18, 25) %>%
pack_rows("HB Beta-Binomial", 26, 33)| Skenario | RSE≥25 | RSE<25 | % <25% | % <15% | CV Mean |
|---|---|---|---|---|---|
|
81 | 0 | 0.0 | 0.0 | 52.28 |
| EBLUP | |||||
| S01 EBLUP-BK All | 67 | 14 | 17.3 | 0.0 | 34.84 |
| S02 EBLUP-BK RSE-NB | 76 | 5 | 6.2 | 0.0 | 50.84 |
| S03 EBLUP-BK RSE-ES | 50 | 31 | 38.3 | 0.0 | 34.14 |
| S04 EBLUP-BK Cluster | 68 | 13 | 16.0 | 0.0 | 36.15 |
| S05 EBLUP-TN All | 73 | 8 | 9.9 | 0.0 | 37.32 |
| S06 EBLUP-TN RSE-NB | 76 | 5 | 6.2 | 0.0 | 50.84 |
| S07 EBLUP-TN RSE-ES | 59 | 22 | 27.2 | 2.5 | 36.80 |
| S08 EBLUP-TN Cluster | 74 | 7 | 8.6 | 0.0 | 39.02 |
| GLMM | |||||
| S09 GLMM-BK All | 30 | 51 | 63.0 | 39.5 | 25.79 |
| S10 GLMM-BK RSE-NB | 4 | 77 | 95.1 | 82.7 | 10.20 |
| S11 GLMM-BK RSE-ES | 4 | 77 | 95.1 | 76.5 | 8.19 |
| S12 GLMM-BK Cluster | 29 | 52 | 64.2 | 46.9 | 24.40 |
| S13 GLMM-TN All | 30 | 51 | 63.0 | 39.5 | 25.10 |
| S14 GLMM-TN RSE-NB | 4 | 77 | 95.1 | 81.5 | 9.90 |
| S15 GLMM-TN RSE-ES | 3 | 78 | 96.3 | 77.8 | 8.28 |
| S16 GLMM-TN Cluster | 24 | 57 | 70.4 | 44.4 | 23.18 |
| EB Beta-Binomial | |||||
| S17 EB-BK All | 80 | 1 | 1.2 | 0.0 | 41.84 |
| S18 EB-BK RSE-NB | 79 | 2 | 2.5 | 0.0 | 40.08 |
| S19 EB-BK RSE-ES | 79 | 2 | 2.5 | 0.0 | 41.76 |
| S20 EB-BK Cluster | 80 | 1 | 1.2 | 0.0 | 41.04 |
| S21 EB-TN All | 80 | 1 | 1.2 | 0.0 | 41.85 |
| S22 EB-TN RSE-NB | 79 | 2 | 2.5 | 0.0 | 40.09 |
| S23 EB-TN RSE-ES | 78 | 3 | 3.7 | 0.0 | 41.71 |
| S24 EB-TN Cluster | 80 | 1 | 1.2 | 0.0 | 40.59 |
| HB Beta-Binomial | |||||
| S25 HB-BK All | 52 | 29 | 35.8 | 0.0 | 28.22 |
| S26 HB-BK RSE-NB | 4 | 77 | 95.1 | 53.1 | 15.53 |
| S27 HB-BK RSE-ES | 4 | 77 | 95.1 | 81.5 | 9.07 |
| S28 HB-BK Cluster | 24 | 57 | 70.4 | 24.7 | 22.47 |
| S29 HB-TN All | 69 | 12 | 14.8 | 0.0 | 29.40 |
| S30 HB-TN RSE-NB | 5 | 76 | 93.8 | 34.6 | 16.74 |
| S31 HB-TN RSE-ES | 3 | 78 | 96.3 | 81.5 | 10.54 |
| S32 HB-TN Cluster | 53 | 28 | 34.6 | 0.0 | 28.00 |
Metrik Akurasi (RB & RMSE)
get_est <- function(df, col) if (!is.null(df) && col %in% names(df)) df[[col]] else NULL
metrik_rows <- list(
metrik(get_est(df_s1, "y_eblup"), df_s1$Estimasi, "S01 EBLUP-BK All"),
metrik(get_est(df_s2, "y_eblup"), df_s2$Estimasi, "S02 EBLUP-BK RSE-NB"),
metrik(get_est(df_s3, "y_eblup"), df_s3$Estimasi, "S03 EBLUP-BK RSE-ES"),
metrik(get_est(df_s4, "y_eblup"), df_s4$Estimasi, "S04 EBLUP-BK Cluster"),
metrik(get_est(df_s5, "y_eblup"), df_s5$Estimasi, "S05 EBLUP-TN All"),
metrik(get_est(df_s6, "y_eblup"), df_s6$Estimasi, "S06 EBLUP-TN RSE-NB"),
metrik(get_est(df_s7, "y_eblup"), df_s7$Estimasi, "S07 EBLUP-TN RSE-ES"),
metrik(get_est(df_s8, "y_eblup"), df_s8$Estimasi, "S08 EBLUP-TN Cluster"),
metrik(get_est(df_s9, "y_glmm"), df_s9$Estimasi, "S09 GLMM-BK All"),
metrik(get_est(df_s10, "y_glmm"), df_s10$Estimasi, "S10 GLMM-BK RSE-NB"),
metrik(get_est(df_s11, "y_glmm"), df_s11$Estimasi, "S11 GLMM-BK RSE-ES"),
metrik(get_est(df_s12, "y_glmm"), df_s12$Estimasi, "S12 GLMM-BK Cluster"),
metrik(get_est(df_s13, "y_glmm"), df_s13$Estimasi, "S13 GLMM-TN All"),
metrik(get_est(df_s14, "y_glmm"), df_s14$Estimasi, "S14 GLMM-TN RSE-NB"),
metrik(get_est(df_s15, "y_glmm"), df_s15$Estimasi, "S15 GLMM-TN RSE-ES"),
metrik(get_est(df_s16, "y_glmm"), df_s16$Estimasi, "S16 GLMM-TN Cluster"),
metrik(get_est(df_s17, "y_eb_pct"), df_s17$Estimasi, "S17 EB-BK All"),
metrik(get_est(df_s18, "y_eb_pct"), df_s18$Estimasi, "S18 EB-BK RSE-NB"),
metrik(get_est(df_s19, "y_eb_pct"), df_s19$Estimasi, "S19 EB-BK RSE-ES"),
metrik(get_est(df_s20, "y_eb_pct"), df_s20$Estimasi, "S20 EB-BK Cluster"),
metrik(get_est(df_s21, "y_eb_pct"), df_s21$Estimasi, "S21 EB-TN All"),
metrik(get_est(df_s22, "y_eb_pct"), df_s22$Estimasi, "S22 EB-TN RSE-NB"),
metrik(get_est(df_s23, "y_eb_pct"), df_s23$Estimasi, "S23 EB-TN RSE-ES"),
metrik(get_est(df_s24, "y_eb_pct"), df_s24$Estimasi, "S24 EB-TN Cluster"),
metrik(get_est(df_s25, "y_hb_pct"), df_s25$Estimasi, "S25 HB-BK All"),
metrik(get_est(df_s26, "y_hb_pct"), df_s26$Estimasi, "S26 HB-BK RSE-NB"),
metrik(get_est(df_s27, "y_hb_pct"), df_s27$Estimasi, "S27 HB-BK RSE-ES"),
metrik(get_est(df_s28, "y_hb_pct"), df_s28$Estimasi, "S28 HB-BK Cluster"),
metrik(get_est(df_s29, "y_hb_pct"), df_s29$Estimasi, "S29 HB-TN All"),
metrik(get_est(df_s30, "y_hb_pct"), df_s30$Estimasi, "S30 HB-TN RSE-NB"),
metrik(get_est(df_s31, "y_hb_pct"), df_s31$Estimasi, "S31 HB-TN RSE-ES"),
metrik(get_est(df_s32, "y_hb_pct"), df_s32$Estimasi, "S32 HB-TN Cluster")
)
tabel_metrik <- bind_rows(Filter(Negate(is.null), metrik_rows))
kable(tabel_metrik,
caption = "Relative Bias (%) & RMSE terhadap Direct Estimator") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"),
full_width = FALSE, font_size = 12) %>%
row_spec(which(abs(tabel_metrik$`RB (%)`) == min(abs(tabel_metrik$`RB (%)`), na.rm=TRUE)),
bold = TRUE, background = "#d4edda")| Skenario | RB (%) | RMSE |
|---|---|---|
| S01 EBLUP-BK All | -8.422 | 4.6325 |
| S02 EBLUP-BK RSE-NB | 6.319 | 5.0908 |
| S03 EBLUP-BK RSE-ES | -20.218 | 5.6677 |
| S04 EBLUP-BK Cluster | -8.550 | 4.4194 |
| S05 EBLUP-TN All | -7.718 | 4.4533 |
| S06 EBLUP-TN RSE-NB | 6.319 | 5.0908 |
| S07 EBLUP-TN RSE-ES | -19.383 | 5.6638 |
| S08 EBLUP-TN Cluster | -7.079 | 4.1715 |
| S09 GLMM-BK All | 10.195 | 3.1854 |
| S10 GLMM-BK RSE-NB | 2.445 | 3.2588 |
| S11 GLMM-BK RSE-ES | 3.032 | 4.3255 |
| S12 GLMM-BK Cluster | 9.899 | 3.5810 |
| S13 GLMM-TN All | 10.797 | 3.2160 |
| S14 GLMM-TN RSE-NB | 2.619 | 3.3256 |
| S15 GLMM-TN RSE-ES | 3.334 | 4.2808 |
| S16 GLMM-TN Cluster | 12.907 | 3.1497 |
| S17 EB-BK All | 2.883 | 2.2960 |
| S18 EB-BK RSE-NB | 0.505 | 2.2617 |
| S19 EB-BK RSE-ES | 0.488 | 2.5714 |
| S20 EB-BK Cluster | 2.853 | 2.4682 |
| S21 EB-TN All | 2.888 | 2.2704 |
| S22 EB-TN RSE-NB | 0.546 | 2.2771 |
| S23 EB-TN RSE-ES | 0.531 | 2.5200 |
| S24 EB-TN Cluster | 3.727 | 2.2786 |
| S25 HB-BK All | 19.361 | 1.3446 |
| S26 HB-BK RSE-NB | 5.948 | 0.6386 |
| S27 HB-BK RSE-ES | 1.097 | 0.0818 |
| S28 HB-BK Cluster | 8.893 | 0.5607 |
| S29 HB-TN All | 25.637 | 1.8293 |
| S30 HB-TN RSE-NB | 6.637 | 0.8475 |
| S31 HB-TN RSE-ES | 1.262 | 0.1314 |
| S32 HB-TN Cluster | 18.857 | 1.4800 |
Barplot % Domain RSE < 25%
direct_pct <- pct_ok(df_base$RSE_direct_col)
plot_df <- tabel_eval %>%
filter(Skenario != "0. Direct") %>%
mutate(
pct_25 = `% < 25%`, # << nama persis dari debug
Model = case_when(
grepl("EBLUP", Skenario) ~ "EBLUP",
grepl("GLMM", Skenario) ~ "GLMM",
grepl("EB-BK|EB-TN", Skenario) ~ "EB Beta",
grepl("HB-BK|HB-TN", Skenario) ~ "HB Beta",
TRUE ~ "Lain"
),
AuxMetode = ifelse(grepl("-BK", Skenario), "Backward", "Top-n/10"),
Skenario = factor(Skenario, levels = rev(unique(Skenario)))
)
ggplot(plot_df, aes(x = Skenario, y = pct_25, fill = Model, alpha = AuxMetode)) +
geom_col(width = 0.72) +
geom_text(aes(label = paste0(pct_25, "%")),
hjust = -0.12, size = 3, fontface = "bold") +
geom_vline(xintercept = direct_pct,
linetype = "dashed", color = "#c0392b", linewidth = 0.9) +
scale_fill_manual(values = c(
"EBLUP" = "#2980b9", "GLMM" = "#8e44ad",
"EB Beta" = "#e67e22", "HB Beta" = "#27ae60"
)) +
scale_alpha_manual(values = c("Backward" = 1, "Top-n/10" = 0.65)) +
scale_y_continuous(limits = c(0, 115), expand = c(0, 0)) +
coord_flip() +
labs(title = "% Domain RSE < 25% — 32 Skenario SAE",
subtitle = "Garis merah = baseline Direct | Gelap = Backward, Terang = Top-n/10",
x = NULL, y = "% Domain RSE < 25%", fill = "Model", alpha = "Seleksi Aux") +
theme_minimal(base_size = 11) +
theme(legend.position = "top", plot.title = element_text(face = "bold", size = 13),
panel.grid.major.y = element_blank())
Scatter RSE Direct vs Model (Facet 32)
make_scatter_df <- function(df, y_col, rse_col, label) {
if (is.null(df) || !all(c(y_col, rse_col, "RSE_direct") %in% names(df))) return(NULL)
df %>%
mutate(RSE_model = .data[[rse_col]], Skenario = label) %>%
dplyr::select(Kako, RSE_direct, RSE_model, Skenario) # << RSE_direct bukan RSE_direct_col
}
scatter_list <- list(
make_scatter_df(df_s1, "y_eblup", "RSE_eblup", "S01 EBLUP-BK All"),
make_scatter_df(df_s2, "y_eblup", "RSE_eblup", "S02 EBLUP-BK RSE-NB"),
make_scatter_df(df_s3, "y_eblup", "RSE_eblup", "S03 EBLUP-BK RSE-ES"),
make_scatter_df(df_s4, "y_eblup", "RSE_eblup", "S04 EBLUP-BK Clust"),
make_scatter_df(df_s5, "y_eblup", "RSE_eblup", "S05 EBLUP-TN All"),
make_scatter_df(df_s6, "y_eblup", "RSE_eblup", "S06 EBLUP-TN RSE-NB"),
make_scatter_df(df_s7, "y_eblup", "RSE_eblup", "S07 EBLUP-TN RSE-ES"),
make_scatter_df(df_s8, "y_eblup", "RSE_eblup", "S08 EBLUP-TN Clust"),
make_scatter_df(df_s9, "y_glmm", "RSE_glmm", "S09 GLMM-BK All"),
make_scatter_df(df_s10, "y_glmm", "RSE_glmm", "S10 GLMM-BK RSE-NB"),
make_scatter_df(df_s11, "y_glmm", "RSE_glmm", "S11 GLMM-BK RSE-ES"),
make_scatter_df(df_s12, "y_glmm", "RSE_glmm", "S12 GLMM-BK Clust"),
make_scatter_df(df_s13, "y_glmm", "RSE_glmm", "S13 GLMM-TN All"),
make_scatter_df(df_s14, "y_glmm", "RSE_glmm", "S14 GLMM-TN RSE-NB"),
make_scatter_df(df_s15, "y_glmm", "RSE_glmm", "S15 GLMM-TN RSE-ES"),
make_scatter_df(df_s16, "y_glmm", "RSE_glmm", "S16 GLMM-TN Clust"),
make_scatter_df(df_s17, "y_eb_pct", "RSE_eb", "S17 EB-BK All"),
make_scatter_df(df_s18, "y_eb_pct", "RSE_eb", "S18 EB-BK RSE-NB"),
make_scatter_df(df_s19, "y_eb_pct", "RSE_eb", "S19 EB-BK RSE-ES"),
make_scatter_df(df_s20, "y_eb_pct", "RSE_eb", "S20 EB-BK Clust"),
make_scatter_df(df_s21, "y_eb_pct", "RSE_eb", "S21 EB-TN All"),
make_scatter_df(df_s22, "y_eb_pct", "RSE_eb", "S22 EB-TN RSE-NB"),
make_scatter_df(df_s23, "y_eb_pct", "RSE_eb", "S23 EB-TN RSE-ES"),
make_scatter_df(df_s24, "y_eb_pct", "RSE_eb", "S24 EB-TN Clust"),
make_scatter_df(df_s25, "y_hb_pct", "RSE_hb", "S25 HB-BK All"),
make_scatter_df(df_s26, "y_hb_pct", "RSE_hb", "S26 HB-BK RSE-NB"),
make_scatter_df(df_s27, "y_hb_pct", "RSE_hb", "S27 HB-BK RSE-ES"),
make_scatter_df(df_s28, "y_hb_pct", "RSE_hb", "S28 HB-BK Clust"),
make_scatter_df(df_s29, "y_hb_pct", "RSE_hb", "S29 HB-TN All"),
make_scatter_df(df_s30, "y_hb_pct", "RSE_hb", "S30 HB-TN RSE-NB"),
make_scatter_df(df_s31, "y_hb_pct", "RSE_hb", "S31 HB-TN RSE-ES"),
make_scatter_df(df_s32, "y_hb_pct", "RSE_hb", "S32 HB-TN Clust")
)
df_scatter <- bind_rows(Filter(Negate(is.null), scatter_list)) %>%
filter(!is.na(RSE_direct) & !is.na(RSE_model))
lim_max <- quantile(c(df_scatter$RSE_direct, df_scatter$RSE_model), 0.99, na.rm = TRUE)
ggplot(df_scatter, aes(x = RSE_direct, y = RSE_model)) +
geom_point(alpha = 0.5, size = 1.3, color = "#2980b9") +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "#c0392b") +
geom_hline(yintercept = 25, linetype = "dotted", color = "#e67e22") +
geom_vline(xintercept = 25, linetype = "dotted", color = "#e67e22") +
facet_wrap(~ Skenario, ncol = 4, scales = "free") +
coord_cartesian(xlim = c(0, lim_max), ylim = c(0, lim_max)) +
labs(title = "RSE Direct vs RSE Model — 32 Skenario SAE",
subtitle = "Di bawah diagonal merah = lebih presisi dari direct",
x = "RSE Direct (%)", y = "RSE Model (%)") +
theme_minimal(base_size = 9) +
theme(strip.text = element_text(face = "bold", size = 7),
plot.title = element_text(face = "bold"))
Heatmap Ringkasan
heat_df <- tabel_eval %>%
filter(Skenario != "0. Direct") %>%
mutate(
pct_25 = `% < 25%`, # << nama persis
Model = case_when(
grepl("EBLUP", Skenario) ~ "EBLUP",
grepl("GLMM", Skenario) ~ "GLMM",
grepl("EB-BK|EB-TN", Skenario) ~ "EB-Beta",
TRUE ~ "HB-Beta"
),
AuxMetode = ifelse(grepl("-BK", Skenario), "Backward", "Top-n/10"),
Partisi = case_when(
grepl("All", Skenario) ~ "All",
grepl("RSE-NB", Skenario) ~ "RSE-NB",
grepl("RSE-ES", Skenario) ~ "RSE-ES",
TRUE ~ "Cluster"
),
Partisi = factor(Partisi, levels = c("All","RSE-NB","RSE-ES","Cluster")),
Model = factor(Model, levels = c("EBLUP","GLMM","EB-Beta","HB-Beta")),
AuxMetode = factor(AuxMetode, levels = c("Backward","Top-n/10"))
)
ggplot(heat_df, aes(x = Partisi, y = Model, fill = pct_25)) +
geom_tile(color = "white", linewidth = 0.8) +
geom_text(aes(label = paste0(pct_25, "%")), size = 4, fontface = "bold") +
facet_wrap(~ AuxMetode, ncol = 2) +
scale_fill_gradient2(low = "#e74c3c", mid = "#f39c12", high = "#27ae60",
midpoint = 60, name = "% Domain\nRSE < 25%") +
labs(title = "Heatmap Performa — % Domain RSE < 25%",
subtitle = "Kiri = Backward | Kanan = Top-n/10",
x = "Partisi", y = "Model") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold", size = 12),
axis.text = element_text(size = 11))
10. Ringkasan & Kesimpulan
best_rse <- tabel_eval %>%
filter(Skenario != "0. Direct") %>%
mutate(pct_25 = `% < 25%`) %>% # << tarik dulu
slice_max(pct_25, n = 3, with_ties = FALSE)
best_rb <- tabel_metrik %>%
slice_min(abs(`RB (%)`), n = 3, with_ties = FALSE)
cat("=== Top 3 Model (% domain RSE < 25%) ===\n")## === Top 3 Model (% domain RSE < 25%) ===## Skenario pct_25 % < 15% CV Mean
## 1 S15 GLMM-TN RSE-ES 96.3 77.8 8.28
## 2 S31 HB-TN RSE-ES 96.3 81.5 10.54
## 3 S10 GLMM-BK RSE-NB 95.1 82.7 10.20##
## === Top 3 Model (Relative Bias terkecil) ===## Skenario RB (%) RMSE
## 1 S19 EB-BK RSE-ES 0.488 2.571
## 2 S18 EB-BK RSE-NB 0.505 2.262
## 3 S23 EB-TN RSE-ES 0.531 2.520Panduan interpretasi 32 skenario:
Seleksi Aux (Backward vs Top-n/10): Backward lebih selektif (stepAIC membuang variabel non-signifikan); Top-n/10 mempertahankan variabel dengan korelasi tertinggi. Keduanya menggunakan
floor(n_partisi/10)variabel awal dari pre-filter korelasi.EBLUP (Normal): Baseline kuat untuk domain tidak terlalu ekstrim proporsinya. Rentan estimasi negatif jika proporsi mendekati 0.
GLMM (Logit, glmmTMB): Menangani bounded response [0,1] lebih proper. Random intercept per domain memberi shrinkage empiris. MSE dari Jackknife leave-one-domain-out.
EB Beta-Binomial: Analytical — cepat dan tidak perlu MCMC. Prior mean domain-spesifik dari regresi logit. Posterior variance analytik dari distribusi Beta konjugat.
HB Beta (saeHB/JAGS): Paling proper secara Bayesian. Memanfaatkan full posterior distribution. Plot posterior dengan 95% CI memperlihatkan ketidakpastian per domain.
Partisi RSE-NB vs RSE-ES: Natural Break (Jenks) memisahkan berdasarkan gap alami distribusi RSE; Equal Size memastikan kedua grup seimbang. RSE-NB lebih heterogen antar grup, RSE-ES lebih stabil dalam estimasi.
Cluster-Aux: Variabel clustering berbeda per kombinasi (Aux method, Link function), membuat cluster lebih relevan secara statistik untuk masing-masing model.
Dianalisis dengan R. Seluruh 32 skenario merupakan kombinasi dari 2 metode seleksi aux × 4 model × 4 partisi.