Dataset incomedata

Dataset incomedata merupakan data pendapatan rumah tangga yang tersedia sebagai built-in di paket sae. Dataset ini digunakan untuk mempraktikkan berbagai metode estimasi area kecil.

# ============================================================
# SLIDE 9 — Estimator Horvitz-Thompson (HT) dengan incomedata
# ============================================================

library(sae)

# Load dataset bawaan package sae
data(incomedata)

# Lihat struktur data
knitr::kable(head(incomedata), caption = "6 Baris Pertama Dataset `incomedata`",
             digits = 2)

str(incomedata)

'data.frame':   17199 obs. of  21 variables:
 $ provlab: Factor w/ 52 levels "Alava","Albacete",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ prov   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ ac     : int  16 16 16 16 16 16 16 16 16 16 ...
 $ gen    : int  2 1 1 1 1 2 1 2 1 2 ...
 $ age    : int  5 3 5 5 3 5 4 4 3 3 ...
 $ nat    : int  2 1 1 1 1 1 1 1 1 1 ...
 $ educ   : int  2 2 2 0 1 1 1 2 0 3 ...
 $ labor  : int  3 1 1 0 1 3 3 3 0 1 ...
 $ age2   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ age3   : int  0 1 0 0 1 0 0 0 1 1 ...
 $ age4   : int  0 0 0 0 0 0 1 1 0 0 ...
 $ age5   : int  1 0 1 1 0 1 0 0 0 0 ...
 $ educ1  : int  0 0 0 0 1 1 1 0 0 0 ...
 $ educ2  : int  1 1 1 0 0 0 0 1 0 0 ...
 $ educ3  : int  0 0 0 0 0 0 0 0 0 1 ...
 $ nat1   : num  0 1 1 1 1 1 1 1 1 1 ...
 $ labor1 : int  0 1 1 0 1 0 0 0 0 1 ...
 $ labor2 : int  0 0 0 0 0 0 0 0 0 0 ...
 $ labor3 : int  1 0 0 0 0 1 1 1 0 0 ...
 $ income : num  12022 8610 3384 2343 9024 ...
 $ weight : num  2804 1348 2568 1168 953 ...

Dataset sizeprov

Dataset sizeprov berisi ukuran populasi per provinsi, digunakan sebagai argumen domsize pada fungsi direct().

data("sizeprov")
knitr::kable(sizeprov, caption = "Dataset `sizeprov` — Ukuran Populasi per Provinsi")

Result 1 — Domain Numerik, domsize kolom 2:3

Pada bagian ini dilakukan estimasi rata-rata pendapatan per provinsi menggunakan Horvitz–Thompson direct estimator dengan desain sampling tanpa pengembalian (WOR) di dalam setiap domain.

Variabel prov digunakan sebagai domain numerik, sedangkan bobot sampling diberikan melalui sweight = weight. Ukuran populasi tiap domain (domsize) diambil dari kolom ke-2 dan ke-3 pada sizeprov, yang berisi kode domain dan ukuran populasinya.

# ============================================================
# RESULT 1
# domain = prov (numerik), dengan bobot & domsize kolom 2:3
# ============================================================
result1 <- direct(
  y        = income,
  dom      = prov,
  sweight  = weight,
  domsize  = sizeprov[, c(2, 3)],
  data     = incomedata
)

knitr::kable(result1, caption = "Result 1 — domain numerik `prov`, domsize kolom 2:3",
             digits = 4)

Result 2 — Variabel Eksplisit $, domsize kolom 1 & 3

Estimasi yang sama dilakukan seperti pada Result 1, namun semua variabel dipanggil secara eksplisit menggunakan operator $ tanpa menyertakan argumen data.

Pada kasus ini, domain menggunakan provlab (label provinsi), dan domsize diambil dari kolom ke-1 dan ke-3 sizeprov, yang menyesuaikan dengan format label domain.

# ============================================================
# RESULT 2
# Sama seperti result1, tapi variabel dipanggil dengan $ eksplisit
# ============================================================
result2 <- direct(
  y        = incomedata$income,
  dom      = incomedata$provlab,
  sweight  = incomedata$weight,
  domsize  = sizeprov[, c(1, 3)]
)

knitr::kable(result2, caption = "Result 2 — variabel eksplisit `$`, domsize kolom 1 & 3",
             digits = 4)

Result 3 — Domain Label provlab, domsize kolom 1 & 3

Pada bagian ini digunakan asumsi Simple Random Sampling Without Replacement (SRSWOR) dalam setiap domain.

Karena desain sampling dianggap SRS, maka bobot sampling (sweight) tidak perlu disertakan. Variabel domain tetap menggunakan provlab, dan ukuran domain (domsize) tetap diperlukan untuk menghitung estimator dan varians berbasis SRSWOR.

# ============================================================
# RESULT 3
# domain = provlab (nama), domsize kolom c(1,3) tanpa bobot eksplisit
# ============================================================
result3 <- direct(
  y        = income,
  dom      = provlab,
  domsize  = sizeprov[, c(1, 3)],
  data     = incomedata
)

knitr::kable(result3, caption = "Result 3 — domain label `provlab`, domsize kolom 1 & 3",
             digits = 4)

Result 4 — Tanpa domsize, Sampling dengan Pengembalian

Pada bagian ini digunakan asumsi Simple Random Sampling With Replacement (SRSWR).

Dalam kondisi ini, ukuran populasi (domsize) tidak diperlukan karena varians dihitung berdasarkan asumsi sampling dengan pengembalian. Argumen `replace = TRUE digunakan untuk mengaktifkan pendekatan SRSWR.

# ============================================================
# RESULT 4
# Tanpa domsize, menggunakan replace = TRUE (SRSWR)
# ============================================================
result4 <- direct(
  y       = income,
  dom     = provlab,
  data    = incomedata,
  replace = TRUE
)

knitr::kable(result4, caption = "Result 4 — tanpa `domsize`, `replace = TRUE` (SRSWR)",
             digits = 4)

Persiapan Data dan Helper

# ================================================================
# SLIDE 12 — ESTIMASI TPT: 4 METODE LENGKAP
# ================================================================

data_tpt <- data.frame(
  j    = 1:15,
  y_ij = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0),
  w_ij = c(665, 808, 1234, 605, 480, 586, 582, 507, 605,
            684, 742, 605, 586, 808, 507)
)

n      <- nrow(data_tpt)
sum_w  <- sum(data_tpt$w_ij)
sum_wy <- sum(data_tpt$w_ij * data_tpt$y_ij)
p_srs  <- sum(data_tpt$y_ij) / n
p_hat  <- sum_wy / sum_w
data_tpt$pi_i <- 1 / data_tpt$w_ij

LINE <- paste0(strrep("=", 68), "\n")
line <- paste0(strrep("-", 68), "\n")
sec  <- function(x) cat(sprintf("\n  [%s]\n", x))

cat(LINE); cat("  DATA INPUT\n"); cat(LINE)
cat(sprintf("  n (ukuran sampel)         = %d\n", n))
cat(sprintf("  Pengangguran (y=1)        = %d  (unit 7, 10, 14)\n", sum(data_tpt$y_ij)))
cat(sprintf("  Σ wi   (pop AK)           = %d\n", sum_w))
cat(sprintf("  Σ wi×yi (pengangguran)    = %d\n\n", sum_wy))

cat(sprintf("  %-4s  %-8s  %-6s  %-10s\n", "j", "Status AK", "w_ij", "(wij)×(yij)"))
cat(line)
for (i in 1:n) {
  status <- ifelse(data_tpt$y_ij[i] == 1, "Nganggur", "Bekerja")
  cat(sprintf("  %-4d  %-8s  %-6d  %-10d\n",
      i, status, data_tpt$w_ij[i], data_tpt$w_ij[i] * data_tpt$y_ij[i]))
}
cat(line)
cat(sprintf("  %-4s  %-8s  %-6d  %-10d\n\n", "TOT", "", sum_w, sum_wy))

====================================================================
  DATA INPUT
====================================================================
  n (ukuran sampel)         = 15
  Pengangguran (y=1)        = 3  (unit 7, 10, 14)
  Σ wi   (pop AK)           = 10004
  Σ wi×yi (pengangguran)    = 2074

  j     Status AK  w_ij    (wij)×(yij)
--------------------------------------------------------------------
  1     Bekerja   665     0         
  2     Bekerja   808     0         
  3     Bekerja   1234    0         
  4     Bekerja   605     0         
  5     Bekerja   480     0         
  6     Bekerja   586     0         
  7     Nganggur  582     582       
  8     Bekerja   507     0         
  9     Bekerja   605     0         
  10    Nganggur  684     684       
  11    Bekerja   742     0         
  12    Bekerja   605     0         
  13    Bekerja   586     0         
  14    Nganggur  808     808       
  15    Bekerja   507     0         
--------------------------------------------------------------------
  TOT             10004   2074      

Metode 1 — SRS Versi PPT

cat(LINE); cat("  METODE 1: SRS — VERSI PPT\n"); cat(LINE)

sec("1A  Estimasi TPT — persis seperti di slide")
cat(line)
cat(sprintf("  Populasi AK = %d\n\n", sum_w))
cat("  Rumus slide:\n")
cat("            Σ (wij × yij)\n")
cat("   TPT  =  ────────────── × 100%\n")
cat("              Σ (wij)\n\n")
cat(sprintf("           %d\n", sum_wy))
cat("   TPT  =  ───────  × 100%\n")
cat(sprintf("           %d\n\n", sum_w))
cat(sprintf("        =  %.1f%%\n", p_hat * 100))

sec("1B  Varians — persis seperti di slide")
cat(line)
vp_m1 <- p_hat * (1 - p_hat) / (n - 1)
cat(sprintf("  v(p) = p(1-p) / (n-1)\n"))
cat(sprintf("       = %.4f × (1 - %.4f) / (%d - 1)\n", p_hat, p_hat, n))
cat(sprintf("       = %.4f × %.4f / %d\n", p_hat, 1 - p_hat, n - 1))
cat(sprintf("       = %.6f\n", vp_m1))

sec("1C  SE & RSE — persis seperti di slide")
cat(line)
SE_m1  <- 0.00587
RSE_m1 <- (SE_m1 / p_hat) * 100
cat(sprintf("  SE  (slide) = 0.00587\n"))
cat(sprintf("  RSE (slide) = (0.00587 / %.4f) × 100 = %.2f%%\n\n", p_hat, RSE_m1))
cat(sprintf("  ► Metode 1: TPT = %.1f%% | SE = %.5f | RSE = %.2f%%\n\n",
            p_hat * 100, SE_m1, RSE_m1))

====================================================================
  METODE 1: SRS — VERSI PPT
====================================================================

  [1A  Estimasi TPT — persis seperti di slide]
--------------------------------------------------------------------
  Populasi AK = 10004

  Rumus slide:
            Σ (wij × yij)
   TPT  =  ────────────── × 100%
              Σ (wij)

           2074
   TPT  =  ───────  × 100%
           10004

        =  20.7%

  [1B  Varians — persis seperti di slide]
--------------------------------------------------------------------
  v(p) = p(1-p) / (n-1)
       = 0.2073 × (1 - 0.2073) / (15 - 1)
       = 0.2073 × 0.7927 / 14
       = 0.011738

  [1C  SE & RSE — persis seperti di slide]
--------------------------------------------------------------------
  SE  (slide) = 0.00587
  RSE (slide) = (0.00587 / 0.2073) × 100 = 2.83%

  ► Metode 1: TPT = 20.7% | SE = 0.00587 | RSE = 2.83%

Metode 2 — SRS yang Benar

cat(LINE); cat("  METODE 2: SRS — SEHARUSNYA (HITUNG BENAR)\n"); cat(LINE)

sec("2A  Estimasi TPT — rumus SRS murni tanpa bobot")
cat(line)
cat(sprintf("  TPT = Σyi / n = %d / %d × 100%%\n", sum(data_tpt$y_ij), n))
cat(sprintf("      = %.4f%%\n", p_srs * 100))

sec("2B  Varians")
cat(line)
vp_m2 <- p_srs * (1 - p_srs) / (n - 1)
cat(sprintf("  v(p) = p(1-p) / (n-1)\n"))
cat(sprintf("       = %.4f × (1 - %.4f) / (%d - 1)\n", p_srs, p_srs, n))
cat(sprintf("       = %.4f × %.4f / %d\n", p_srs, 1 - p_srs, n - 1))
cat(sprintf("       = %.6f\n", vp_m2))

sec("2C  SE — koreksi kesalahan slide")
cat(line)
SE_m2  <- sqrt(vp_m2)
RSE_m2 <- (SE_m2 / p_srs) * 100
cat(sprintf("  SE  = √v(p) = √%.6f = %.6f\n\n", vp_m2, SE_m2))
cat("  ⚠  Kesalahan di slide:\n")
cat(sprintf("     √0.012 yang benar = %.6f\n", sqrt(vp_m2)))
cat(sprintf("     Slide tulis 0.00587, padahal 0.00587² = %.7f ≠ %.4f\n",
            0.00587^2, vp_m2))

sec("2D  RSE")
cat(line)
cat(sprintf("  RSE = (SE / p) × 100\n"))
cat(sprintf("      = (%.6f / %.4f) × 100\n", SE_m2, p_srs))
cat(sprintf("      = %.4f%%\n\n", RSE_m2))
cat(sprintf("  ► Metode 2: TPT = %.4f%% | SE = %.6f | RSE = %.4f%%\n\n",
            p_srs * 100, SE_m2, RSE_m2))

====================================================================
  METODE 2: SRS — SEHARUSNYA (HITUNG BENAR)
====================================================================

  [2A  Estimasi TPT — rumus SRS murni tanpa bobot]
--------------------------------------------------------------------
  TPT = Σyi / n = 3 / 15 × 100%
      = 20.0000%

  [2B  Varians]
--------------------------------------------------------------------
  v(p) = p(1-p) / (n-1)
       = 0.2000 × (1 - 0.2000) / (15 - 1)
       = 0.2000 × 0.8000 / 14
       = 0.011429

  [2C  SE — koreksi kesalahan slide]
--------------------------------------------------------------------
  SE  = √v(p) = √0.011429 = 0.106904

  ⚠  Kesalahan di slide:
     √0.012 yang benar = 0.106904
     Slide tulis 0.00587, padahal 0.00587² = 0.0000345 ≠ 0.0114

  [2D  RSE]
--------------------------------------------------------------------
  RSE = (SE / p) × 100
      = (0.106904 / 0.2000) × 100
      = 53.4522%

  ► Metode 2: TPT = 20.0000% | SE = 0.106904 | RSE = 53.4522%

Metode 3 — HT Estimasi Total (Dibedah Lengkap)

cat(LINE); cat("  METODE 3: HT — ESTIMASI TOTAL (DIBEDAH LENGKAP)\n"); cat(LINE)

sec("3A  Estimasi TPT — Ratio HT Estimator")
cat(line)
cat(sprintf("  Ŷ_HT = Σ(wi × yi) = %d\n", sum_wy))
cat(sprintf("  X̂_HT = Σ(wi)      = %d\n", sum_w))
cat(sprintf("  TPT  = %d / %d × 100%% = %.4f%%\n\n", sum_wy, sum_w, p_hat * 100))

sec("3B  Varians — SUKU-1: Σ(wi² - wi) × yi²")
cat(line)
data_tpt$koef1   <- data_tpt$w_ij^2 - data_tpt$w_ij
data_tpt$term1_i <- data_tpt$koef1 * data_tpt$y_ij^2

cat(sprintf("  %-4s  %-5s  %-6s  %-12s  %-12s  %-5s  %-12s\n",
            "j","y_ij","wi","wi²","wi²-wi","yi²","kontribusi"))
cat(sprintf("  %s\n", strrep("-", 62)))
for (i in 1:n) {
  cat(sprintf("  %-4d  %-5d  %-6d  %-12s  %-12s  %-5d  %-12.2f\n",
    data_tpt$j[i], data_tpt$y_ij[i], data_tpt$w_ij[i],
    format(data_tpt$w_ij[i]^2, big.mark=","),
    format(data_tpt$koef1[i],  big.mark=","),
    data_tpt$y_ij[i]^2, data_tpt$term1_i[i]))
}
SUKU1 <- sum(data_tpt$term1_i)
cat(sprintf("  %s\n", strrep("-", 62)))
cat(sprintf("  %-38s  %-12s\n", "SUKU-1 = Σ kontribusi",
            format(round(SUKU1), big.mark=",")))

cat("\n  Detail unit y=1 (satu-satunya yang berkontribusi):\n")
for (i in which(data_tpt$y_ij == 1)) {
  w <- data_tpt$w_ij[i]
  cat(sprintf("    Unit %-2d (w=%d):  %s - %d  =  %s\n",
      data_tpt$j[i], w,
      format(w^2, big.mark=","), w,
      format(w^2-w, big.mark=",")))
}

sec("3C  Varians — SUKU-2: 2 × ΣΣ(wi×wj - 1/πij) × yi × yj")
cat(line)
cat("  Hanya pasangan (i,j) di mana yi=1 DAN yj=1 yang aktif\n")
cat("  → Pasangan aktif: (7,10), (7,14), (10,14)\n\n")

N_est  <- sum_w
pi_ij  <- (n * (n-1)) / (N_est * (N_est-1))
cat(sprintf("  πij (SRSWOR) = n(n-1) / N(N-1)\n"))
cat(sprintf("               = %d×%d / %d×%d\n", n, n-1, N_est, N_est-1))
cat(sprintf("               = %d / %s\n", n*(n-1), format(N_est*(N_est-1), big.mark=",")))
cat(sprintf("               = %.10f\n\n", pi_ij))

cat(sprintf("  %-12s  %-6s  %-6s  %-12s  %-14s  %-14s  %-10s\n",
            "Pasangan","wi","wj","wi×wj","1/πij","wi×wj-1/πij","2×kontrib"))
cat(sprintf("  %s\n", strrep("-", 78)))

unit_y1   <- which(data_tpt$y_ij == 1)
SUKU2_raw <- 0
for (a in 1:(length(unit_y1)-1)) {
  for (b in (a+1):length(unit_y1)) {
    i  <- unit_y1[a]; j <- unit_y1[b]
    wi <- data_tpt$w_ij[i]; wj <- data_tpt$w_ij[j]
    wiwj      <- wi * wj
    inv_pi_ij <- 1 / pi_ij
    koef      <- wiwj - inv_pi_ij
    kontrib   <- koef * data_tpt$y_ij[i] * data_tpt$y_ij[j]
    SUKU2_raw <- SUKU2_raw + kontrib
    cat(sprintf("  (%-2d, %-2d)      %-6d  %-6d  %-12s  %-14s  %-14s  %-10.2f\n",
      data_tpt$j[i], data_tpt$j[j], wi, wj,
      format(round(wiwj),      big.mark=","),
      format(round(inv_pi_ij), big.mark=","),
      format(round(koef),      big.mark=","),
      2 * kontrib))
  }
}
SUKU2 <- 2 * SUKU2_raw
cat(sprintf("  %s\n", strrep("-", 78)))
cat(sprintf("  SUKU-2 = 2 × %.2f = %.2f\n", SUKU2_raw, SUKU2))

sec("3D  Total v(Ŷ_HT) → v(TPT) → SE → RSE")
cat(line)
v_HT_total <- SUKU1 + SUKU2
var_tpt_3  <- v_HT_total / sum_w^2
SE_m3      <- sqrt(var_tpt_3) * 100
RSE_m3     <- (SE_m3 / (p_hat * 100)) * 100

cat(sprintf("  v(Ŷ_HT) = SUKU-1 + SUKU-2\n"))
cat(sprintf("          = %s + (%s)\n",
            format(round(SUKU1), big.mark=","),
            format(round(SUKU2), big.mark=",")))
cat(sprintf("          = %s\n\n", format(round(v_HT_total), big.mark=",")))
cat(sprintf("  v(TPT)  = v(Ŷ_HT) / (Σwi)²\n"))
cat(sprintf("          = %s / %s\n",
            format(round(v_HT_total), big.mark=","),
            format(sum_w^2, big.mark=",")))
cat(sprintf("          = %.6f\n\n", var_tpt_3))
cat(sprintf("  SE      = √%.6f × 100 = %.4f%%\n", var_tpt_3, SE_m3))
cat(sprintf("  RSE     = (%.4f / %.4f) × 100 = %.4f%%\n\n", SE_m3, p_hat*100, RSE_m3))
cat(sprintf("  ► Metode 3: TPT = %.4f%% | SE = %.4f%% | RSE = %.4f%%\n\n",
            p_hat*100, SE_m3, RSE_m3))

====================================================================
  METODE 3: HT — ESTIMASI TOTAL (DIBEDAH LENGKAP)
====================================================================

  [3A  Estimasi TPT — Ratio HT Estimator]
--------------------------------------------------------------------
  Ŷ_HT = Σ(wi × yi) = 2074
  X̂_HT = Σ(wi)      = 10004
  TPT  = 2074 / 10004 × 100% = 20.7317%


  [3B  Varians — SUKU-1: Σ(wi² - wi) × yi²]
--------------------------------------------------------------------
  j     y_ij   wi      wi²          wi²-wi       yi²   kontribusi  
  --------------------------------------------------------------
  1     0      665     442,225       441,560       0      0.00        
  2     0      808     652,864       652,056       0      0.00        
  3     0      1234    1,522,756     1,521,522     0      0.00        
  4     0      605     366,025       365,420       0      0.00        
  5     0      480     230,400       229,920       0      0.00        
  6     0      586     343,396       342,810       0      0.00        
  7     1      582     338,724       338,142       1      338142.00   
  8     0      507     257,049       256,542       0      0.00        
  9     0      605     366,025       365,420       0      0.00        
  10    1      684     467,856       467,172       1      467172.00   
  11    0      742     550,564       549,822       0      0.00        
  12    0      605     366,025       365,420       0      0.00        
  13    0      586     343,396       342,810       0      0.00        
  14    1      808     652,864       652,056       1      652056.00   
  15    0      507     257,049       256,542       0      0.00        
  --------------------------------------------------------------
  SUKU-1 = Σ kontribusi                  1,457,370   

  Detail unit y=1 (satu-satunya yang berkontribusi):
    Unit 7  (w=582):  338,724 - 582  =  338,142
    Unit 10 (w=684):  467,856 - 684  =  467,172
    Unit 14 (w=808):  652,864 - 808  =  652,056

  [3C  Varians — SUKU-2: 2 × ΣΣ(wi×wj - 1/πij) × yi × yj]
--------------------------------------------------------------------
  Hanya pasangan (i,j) di mana yi=1 DAN yj=1 yang aktif
  → Pasangan aktif: (7,10), (7,14), (10,14)

  πij (SRSWOR) = n(n-1) / N(N-1)
               = 15×14 / 10004×10003
               = 210 / 100,070,012
               = 0.0000020985

  Pasangan      wi      wj      wi×wj        1/πij          wi×wj-1/πij   2×kontrib
  ------------------------------------------------------------------------------
  (7 , 10)      582     684     398,088       476,524         -78,436         -156871.73
  (7 , 14)      582     808     470,256       476,524         -6,268          -12535.73 
  (10, 14)      684     808     552,672       476,524         76,148          152296.27 
  ------------------------------------------------------------------------------
  SUKU-2 = 2 × -8555.60 = -17111.20

  [3D  Total v(Ŷ_HT) → v(TPT) → SE → RSE]
--------------------------------------------------------------------
  v(Ŷ_HT) = SUKU-1 + SUKU-2
          = 1,457,370 + (-17,111)
          = 1,440,259

  v(TPT)  = v(Ŷ_HT) / (Σwi)²
          = 1,440,259 / 100,080,016
          = 0.014391

  SE      = √0.014391 × 100 = 11.9963%
  RSE     = (11.9963 / 20.7317) × 100 = 57.8644%

  ► Metode 3: TPT = 20.7317% | SE = 11.9963% | RSE = 57.8644%

Metode 4 — HT Delta Method (Ratio Estimator)

cat(LINE); cat("  METODE 4: HT — DELTA METHOD (RATIO ESTIMATOR)\n"); cat(LINE)

sec("4A  Estimasi TPT")
cat(line)
cat(sprintf("  TPT = %.4f%%  (sama dengan Metode 3)\n", p_hat * 100))

sec("4B  Kenapa Metode 3 over-estimasi varians?")
cat(line)
cat("  Metode 3 → v(Ŷ_HT) untuk estimasi TOTAL → yi² pakai nilai asli\n")
cat("  TPT adalah RASIO: TPT = Ŷ/X̂  →  harus pakai Delta Method\n\n")
cat("  Kunci: ganti  yi  →  ei = yi - p_hat  sebelum hitung varians\n")
cat(sprintf("  Karena semua unit adalah AK (xi=1): ei = yi - %.4f\n", p_hat))

sec("4C  Varians SUKU-1: Σ(wi² - wi) × ei²")
cat(line)
data_tpt$e_i    <- data_tpt$y_ij - p_hat
data_tpt$c_delt <- (data_tpt$w_ij^2 - data_tpt$w_ij) * data_tpt$e_i^2

cat(sprintf("  e_i jika y=1:  1 - %.4f =  %.4f\n",  p_hat,  1 - p_hat))
cat(sprintf("  e_i jika y=0:  0 - %.4f = %.4f\n\n", p_hat, -p_hat))

cat(sprintf("  %-4s  %-5s  %-6s  %-8s  %-12s  %-8s  %-12s\n",
            "j","y_ij","wi","e_i","wi²-wi","ei²","kontribusi"))
cat(sprintf("  %s\n", strrep("-", 62)))
for (i in 1:n) {
  cat(sprintf("  %-4d  %-5d  %-6d  %-8.4f  %-12s  %-8.4f  %-12.2f\n",
    data_tpt$j[i], data_tpt$y_ij[i], data_tpt$w_ij[i],
    data_tpt$e_i[i],
    format(data_tpt$w_ij[i]^2 - data_tpt$w_ij[i], big.mark=","),
    data_tpt$e_i[i]^2, data_tpt$c_delt[i]))
}
v_d1 <- sum(data_tpt$c_delt)
cat(sprintf("  %s\n", strrep("-", 62)))
cat(sprintf("  %-38s  %-12s\n", "SUKU-1(delta) = Σ kontribusi",
            format(round(v_d1), big.mark=",")))

sec("4D  Varians SUKU-2: 2 × ΣΣ(wi×wj - 1/πij) × ei × ej")
cat(line)
cat("  Semua pasangan aktif karena ei ≠ 0 untuk semua unit\n\n")

cat(sprintf("  %-10s  %-6s  %-6s  %-8s  %-8s  %-14s  %-8s  %-10s\n",
            "Pasangan","wi","wj","ei","ej","wi×wj-1/πij","ei×ej","2×kontrib"))
cat(sprintf("  %s\n", strrep("-", 82)))

d2_raw <- 0; shown <- 0
for (i in 1:(n-1)) {
  for (j in (i+1):n) {
    wi <- data_tpt$w_ij[i]; wj <- data_tpt$w_ij[j]
    ei <- data_tpt$e_i[i];  ej <- data_tpt$e_i[j]
    koef    <- wi*wj - 1/pi_ij
    kontrib <- koef * ei * ej
    d2_raw  <- d2_raw + kontrib
    shown   <- shown + 1
    if (shown <= 10 || data_tpt$y_ij[i]==1 || data_tpt$y_ij[j]==1) {
      cat(sprintf("  (%-2d,%-2d)    %-6d  %-6d  %-8.4f  %-8.4f  %-14s  %-8.4f  %-10.2f\n",
        data_tpt$j[i], data_tpt$j[j], wi, wj, ei, ej,
        format(round(koef), big.mark=","), ei*ej, 2*kontrib))
    } else if (shown == 11) {
      cat(sprintf("  ... (%d pasangan lainnya dihitung tapi tidak ditampilkan)\n",
                  n*(n-1)/2 - 10))
    }
  }
}
d2 <- 2 * d2_raw
cat(sprintf("  %s\n", strrep("-", 82)))
cat(sprintf("  Total %d pasangan  |  SUKU-2(delta) = 2 × %.2f = %.2f\n",
            n*(n-1)/2, d2_raw, d2))

sec("4E  Total v(TPT) → SE → RSE")
cat(line)
v_d_tot  <- v_d1 + d2
var_tpt4 <- v_d_tot / sum_w^2
SE_m4    <- sqrt(var_tpt4) * 100
RSE_m4   <- (SE_m4 / (p_hat*100)) * 100

cat(sprintf("  v_HT(ê) = SUKU-1 + SUKU-2\n"))
cat(sprintf("          = %s + (%s)\n",
            format(round(v_d1), big.mark=","),
            format(round(d2),   big.mark=",")))
cat(sprintf("          = %s\n\n", format(round(v_d_tot), big.mark=",")))
cat(sprintf("  v(TPT)  = %s / %s\n",
            format(round(v_d_tot), big.mark=","),
            format(sum_w^2, big.mark=",")))
cat(sprintf("          = %.6f\n\n", var_tpt4))
cat(sprintf("  SE      = √%.6f × 100 = %.4f%%\n", var_tpt4, SE_m4))
cat(sprintf("  RSE     = (%.4f / %.4f) × 100 = %.4f%%\n\n", SE_m4, p_hat*100, RSE_m4))
cat(sprintf("  ► Metode 4: TPT = %.4f%% | SE = %.4f%% | RSE = %.4f%%\n\n",
            p_hat*100, SE_m4, RSE_m4))

====================================================================
  METODE 4: HT — DELTA METHOD (RATIO ESTIMATOR)
====================================================================

  [4A  Estimasi TPT]
--------------------------------------------------------------------
  TPT = 20.7317%  (sama dengan Metode 3)

  [4B  Kenapa Metode 3 over-estimasi varians?]
--------------------------------------------------------------------
  Metode 3 → v(Ŷ_HT) untuk estimasi TOTAL → yi² pakai nilai asli
  TPT adalah RASIO: TPT = Ŷ/X̂  →  harus pakai Delta Method

  Kunci: ganti  yi  →  ei = yi - p_hat  sebelum hitung varians
  Karena semua unit adalah AK (xi=1): ei = yi - 0.2073

  [4C  Varians SUKU-1: Σ(wi² - wi) × ei²]
--------------------------------------------------------------------
  e_i jika y=1:  1 - 0.2073 =  0.7927
  e_i jika y=0:  0 - 0.2073 = -0.2073

  j     y_ij   wi      e_i       wi²-wi       ei²      kontribusi  
  --------------------------------------------------------------
  1     0      665     -0.2073   441,560       0.0430    18978.41    
  2     0      808     -0.2073   652,056       0.0430    28025.61    
  3     0      1234    -0.2073   1,521,522     0.0430    65395.58    
  4     0      605     -0.2073   365,420       0.0430    15705.89    
  5     0      480     -0.2073   229,920       0.0430    9882.05     
  6     0      586     -0.2073   342,810       0.0430    14734.10    
  7     1      582     0.7927    338,142       0.6283    212470.25   
  8     0      507     -0.2073   256,542       0.0430    11026.27    
  9     0      605     -0.2073   365,420       0.0430    15705.89    
  10    1      684     0.7927    467,172       0.6283    293545.76   
  11    0      742     -0.2073   549,822       0.0430    23631.55    
  12    0      605     -0.2073   365,420       0.0430    15705.89    
  13    0      586     -0.2073   342,810       0.0430    14734.10    
  14    1      808     0.7927    652,056       0.6283    409716.92   
  15    0      507     -0.2073   256,542       0.0430    11026.27    
  --------------------------------------------------------------
  SUKU-1(delta) = Σ kontribusi           1,160,285   

  [4D  Varians SUKU-2: 2 × ΣΣ(wi×wj - 1/πij) × ei × ej]
--------------------------------------------------------------------
  Semua pasangan aktif karena ei ≠ 0 untuk semua unit

  Pasangan    wi      wj      ei        ej        wi×wj-1/πij   ei×ej    2×kontrib
  ----------------------------------------------------------------------------------
  (1 ,2 )    665     808     -0.2073   -0.2073   60,796          0.0430    5226.08   
  (1 ,3 )    665     1234    -0.2073   -0.2073   344,086         0.0430    29577.90  
  (1 ,4 )    665     605     -0.2073   -0.2073   -74,199         0.0430    -6378.19  
  (1 ,5 )    665     480     -0.2073   -0.2073   -157,324        0.0430    -13523.68 
  (1 ,6 )    665     586     -0.2073   -0.2073   -86,834         0.0430    -7464.30  
  (1 ,7 )    665     582     -0.2073   0.7927    -89,494         -0.1643   29414.25  
  (1 ,8 )    665     507     -0.2073   -0.2073   -139,369        0.0430    -11980.25 
  (1 ,9 )    665     605     -0.2073   -0.2073   -74,199         0.0430    -6378.19  
  (1 ,10)    665     684     -0.2073   0.7927    -21,664         -0.1643   7120.34   
  (1 ,11)    665     742     -0.2073   -0.2073   16,906          0.0430    1453.26   
  ... (95 pasangan lainnya dihitung tapi tidak ditampilkan)
  (1 ,14)    665     808     -0.2073   0.7927    60,796          -0.1643   -19982.07 
  (2 ,7 )    808     582     -0.2073   0.7927    -6,268          -0.1643   2060.08   
  (2 ,10)    808     684     -0.2073   0.7927    76,148          -0.1643   -25027.87 
  (2 ,14)    808     808     -0.2073   0.7927    176,340         -0.1643   -57958.31 
  (3 ,7 )    1234    582     -0.2073   0.7927    241,664         -0.1643   -79428.57 
  (3 ,10)    1234    684     -0.2073   0.7927    367,532         -0.1643   -120798.04
  (3 ,14)    1234    808     -0.2073   0.7927    520,548         -0.1643   -171090.33
  (4 ,7 )    605     582     -0.2073   0.7927    -124,414        -0.1643   40891.53  
  (4 ,10)    605     684     -0.2073   0.7927    -62,704         -0.1643   20609.09  
  (4 ,14)    605     808     -0.2073   0.7927    12,316          -0.1643   -4047.99  
  (5 ,7 )    480     582     -0.2073   0.7927    -197,164        -0.1643   64802.52  
  (5 ,10)    480     684     -0.2073   0.7927    -148,204        -0.1643   48710.67  
  (5 ,14)    480     808     -0.2073   0.7927    -88,684         -0.1643   29148.03  
  (6 ,7 )    586     582     -0.2073   0.7927    -135,472        -0.1643   44526.00  
  (6 ,10)    586     684     -0.2073   0.7927    -75,700         -0.1643   24880.53  
  (6 ,14)    586     808     -0.2073   0.7927    -3,036          -0.1643   997.81    
  (7 ,8 )    582     507     0.7927    -0.2073   -181,450        -0.1643   59637.75  
  (7 ,9 )    582     605     0.7927    -0.2073   -124,414        -0.1643   40891.53  
  (7 ,10)    582     684     0.7927    0.7927    -78,436         0.6283    -98569.76 
  (7 ,11)    582     742     0.7927    -0.2073   -44,680         -0.1643   14685.08  
  (7 ,12)    582     605     0.7927    -0.2073   -124,414        -0.1643   40891.53  
  (7 ,13)    582     586     0.7927    -0.2073   -135,472        -0.1643   44526.00  
  (7 ,14)    582     808     0.7927    0.7927    -6,268          0.6283    -7876.78  
  (7 ,15)    582     507     0.7927    -0.2073   -181,450        -0.1643   59637.75  
  (8 ,10)    507     684     -0.2073   0.7927    -129,736        -0.1643   42640.73  
  (8 ,14)    507     808     -0.2073   0.7927    -66,868         -0.1643   21977.69  
  (9 ,10)    605     684     -0.2073   0.7927    -62,704         -0.1643   20609.09  
  (9 ,14)    605     808     -0.2073   0.7927    12,316          -0.1643   -4047.99  
  (10,11)    684     742     0.7927    -0.2073   31,004          -0.1643   -10190.23 
  (10,12)    684     605     0.7927    -0.2073   -62,704         -0.1643   20609.09  
  (10,13)    684     586     0.7927    -0.2073   -75,700         -0.1643   24880.53  
  (10,14)    684     808     0.7927    0.7927    76,148          0.6283    95694.78  
  (10,15)    684     507     0.7927    -0.2073   -129,736        -0.1643   42640.73  
  (11,14)    742     808     -0.2073   0.7927    123,012         -0.1643   -40430.82 
  (12,14)    605     808     -0.2073   0.7927    12,316          -0.1643   -4047.99  
  (13,14)    586     808     -0.2073   0.7927    -3,036          -0.1643   997.81    
  (14,15)    808     507     0.7927    -0.2073   -66,868         -0.1643   21977.69  
  ----------------------------------------------------------------------------------
  Total 105 pasangan  |  SUKU-2(delta) = 2 × -11814.49 = -23628.98

  [4E  Total v(TPT) → SE → RSE]
--------------------------------------------------------------------
  v_HT(ê) = SUKU-1 + SUKU-2
          = 1,160,285 + (-23,629)
          = 1,136,656

  v(TPT)  = 1,136,656 / 100,080,016
          = 0.011357

  SE      = √0.011357 × 100 = 10.6571%
  RSE     = (10.6571 / 20.7317) × 100 = 51.4050%

  ► Metode 4: TPT = 20.7317% | SE = 10.6571% | RSE = 51.4050%

Tabel Perbandingan 4 Metode

# Tabel ringkasan menggunakan knitr::kable
tabel_metode <- data.frame(
  Metode = c("1. SRS versi PPT (berbobot)",
             "2. SRS benar (tanpa bobot)",
             "3. HT Estimasi Total",
             "4. HT Delta Method"),
  TPT_pct  = round(c(p_hat*100, p_srs*100, p_hat*100, p_hat*100), 4),
  SE_pct   = round(c(SE_m1, SE_m2, SE_m3, SE_m4), 6),
  RSE_pct  = round(c(RSE_m1, RSE_m2, RSE_m3, RSE_m4), 4),
  Keterangan = c("TPT berbobot; SE slide salah",
                 "Murni SRS tanpa bobot",
                 "Over-estimasi varians",
                 "✓ Benar untuk ratio estimator")
)

knitr::kable(
  tabel_metode,
  caption = "Ringkasan Perbandingan 4 Metode Estimasi TPT",
  col.names = c("Metode", "TPT (%)", "SE (%)", "RSE (%)", "Keterangan"),
  digits = 4,
  align = c("l","r","r","r","l")
)

Catatan Metodologis

cat("  CATATAN:\n")
cat("  [M1] Slide berjudul 'SRS' tapi pakai rumus berbobot (HT).\n")
cat("       v(p) dihitung dari p_hat=0.2073. SE ditulis 0.00587\n")
cat(sprintf("       padahal √%.4f = %.5f → RSE slide tidak valid.\n\n",
            vp_m1, sqrt(vp_m1)))
cat("  [M2] SRS murni: TPT = Σyi/n = 20.00%% (tanpa bobot).\n")
cat(sprintf("       SE yang benar = √%.4f = %.5f.\n\n", vp_m2, SE_m2))
cat("  [M3] Rumus HT lengkap (Suku-1+Suku-2) untuk estimasi TOTAL.\n")
cat("       Jika dipakai untuk rasio TPT → varians over-estimasi.\n\n")
cat("  [M4] Delta Method: ganti yi → ei=yi-p_hat sebelum hitung\n")
cat("       varians. Ini metode yang benar untuk ratio estimator.\n\n")
cat("  RSE ~50%% pada M2-M4 adalah WAJAR: n=15 sangat kecil,\n")
cat("  hanya 3 dari 15 unit yang pengangguran.\n")

  CATATAN:
  [M1] Slide berjudul 'SRS' tapi pakai rumus berbobot (HT).
       v(p) dihitung dari p_hat=0.2073. SE ditulis 0.00587
       padahal √0.0117 = 0.10834 → RSE slide tidak valid.

  [M2] SRS murni: TPT = Σyi/n = 20.00%% (tanpa bobot).
       SE yang benar = √0.0114 = 0.10690.

  [M3] Rumus HT lengkap (Suku-1+Suku-2) untuk estimasi TOTAL.
       Jika dipakai untuk rasio TPT → varians over-estimasi.

  [M4] Delta Method: ganti yi → ei=yi-p_hat sebelum hitung
       varians. Ini metode yang benar untuk ratio estimator.

  RSE ~50%% pada M2-M4 adalah WAJAR: n=15 sangat kecil,
  hanya 3 dari 15 unit yang pengangguran.

Package dan Opsi Global

# ================================================================
# SLIDE 18–19 — Estimasi Langsung TPT dengan Taylor Linearization
# ================================================================

library(readxl)
library(survey)
library(haven)
library(dplyr)
library(writexl)

options(scipen = 99)
options(survey.lonely.psu = "adjust")

# Pilihan lonely.psu:
#   "adjust"    → SE = 0 di stratum itu (REKOMENDASI BPS)
#   "average"   → ambil rata-rata SE stratum lain
#   "remove"    → hapus stratum dari perhitungan
#   "certainty" → anggap PSU pasti terpilih

Load Data SAKERNAS

cat("================================================================\n")
cat("  LOAD DATA\n")
cat("================================================================\n")

data_coba <- read_excel("C:/Users/User/Downloads/DATASAKERNASSIAPOLAH.xlsx")

cat(sprintf("  Jumlah baris   = %d\n",   nrow(data_coba)))
cat(sprintf("  Jumlah kolom   = %d\n",   ncol(data_coba)))
cat(sprintf("  Kolom tersedia : %s\n\n", paste(names(data_coba), collapse = ", ")))

================================================================
  LOAD DATA
================================================================
  Jumlah baris   = 12129
  Jumlah kolom   = 20
  Kolom tersedia : tahun, urutan, SUPASfinal_weightR, KODE_PROV, NAMA_PROV, KODE_KAB2, NAMA_KAB, KODE_KEC2, NAMA_KEC, KLASIFIKASI2, no_urutNKS, NURT2, B1_R11, B4_K8, B5_R5A1, B5_R6, B5_R15A, B5_R15B, B5_R20A, STATUS

Variabel Turunan

Dari kolom STATUS dibuat dua variabel biner: pengangguran dan angkatankerja.

cat("================================================================\n")
cat("  VARIABEL TURUNAN\n")
cat("================================================================\n")

data_coba$pengangguran <- as.numeric(
  ifelse(data_coba$STATUS == "PENGANGGURAN", "1", "0")
)

data_coba$angkatankerja <- as.numeric(
  ifelse(data_coba$STATUS == "PENGANGGURAN", "1",
    ifelse(data_coba$STATUS == "BEKERJA",    "1", "0"))
)

data_coba$pengangguran[is.na(data_coba$pengangguran)] <- 0

cat("  Distribusi STATUS:\n")
print(table(data_coba$STATUS))
cat(sprintf("\n  Total pengangguran (n)   = %d\n",  sum(data_coba$pengangguran,  na.rm = TRUE)))
cat(sprintf("  Total angkatan kerja (n) = %d\n\n", sum(data_coba$angkatankerja, na.rm = TRUE)))

================================================================
  VARIABEL TURUNAN
================================================================
  Distribusi STATUS:

             BEKERJA BUKAN ANGKATAN KERJA         PENGANGGURAN 
                6856                 4633                  640 

  Total pengangguran (n)   = 640
  Total angkatan kerja (n) = 7496

Definisi Strata, PSU, dan SSU

cat("================================================================\n")
cat("  DESAIN SAMPLING SAKERNAS\n")
cat("================================================================\n")
cat("  Strata : KODE_PROV + KLASIFIKASI2 (perkotaan/perdesaan)\n")
cat("  PSU    : Blok Sensus  (no_urutNKS)\n")
cat("  SSU    : Rumahtangga  (NURT2)\n")
cat("  Bobot  : SUPASfinal_weightR\n\n")

data_coba$STAT   <- paste(data_coba$KODE_PROV, data_coba$KODE_KAB2, data_coba$KODE_KEC2)
data_coba$PSU    <- paste(data_coba$STAT, data_coba$no_urutNKS)
data_coba$SSU    <- paste(data_coba$PSU,  data_coba$NURT2)
data_coba$STRATA <- paste(data_coba$KODE_PROV, data_coba$KLASIFIKASI2)

cat(sprintf("  Jumlah kecamatan (STAT)  = %d\n", length(unique(data_coba$STAT))))
cat(sprintf("  Jumlah PSU (blok sensus) = %d\n", length(unique(data_coba$PSU))))
cat(sprintf("  Jumlah SSU (rumahtangga) = %d\n", length(unique(data_coba$SSU))))
cat(sprintf("  Jumlah strata            = %d\n\n", length(unique(data_coba$STRATA))))

================================================================
  DESAIN SAMPLING SAKERNAS
================================================================
  Strata : KODE_PROV + KLASIFIKASI2 (perkotaan/perdesaan)
  PSU    : Blok Sensus  (no_urutNKS)
  SSU    : Rumahtangga  (NURT2)
  Bobot  : SUPASfinal_weightR

  Jumlah kecamatan (STAT)  = 144
  Jumlah PSU (blok sensus) = 448
  Jumlah SSU (rumahtangga) = 12129
  Jumlah strata            = 2

Penanganan Lonely PSU

psu_per_strata <- data_coba %>%
  group_by(STRATA) %>%
  summarise(n_PSU = n_distinct(PSU), .groups = "drop")

lonely <- psu_per_strata %>% filter(n_PSU == 1)
cat(sprintf("  Strata dengan 1 PSU (lonely): %d strata\n", nrow(lonely)))
if (nrow(lonely) > 0) {
  cat("  Daftar lonely PSU strata:\n")
  print(lonely)
}
cat(sprintf("  → Ditangani dengan: options(survey.lonely.psu = 'adjust')\n\n"))

  Strata dengan 1 PSU (lonely): 0 strata
  → Ditangani dengan: options(survey.lonely.psu = 'adjust')

Objek Desain svydesign

sakernas_coba <- svydesign(
  id      = ~PSU + SSU,
  strata  = ~STRATA,
  data    = data_coba,
  weights = ~SUPASfinal_weightR
)

data_coba$PSU <- as.factor(data_coba$PSU)
data_coba$SSU <- as.factor(data_coba$SSU)

cat("  Ringkasan desain sampling:\n")
print(summary(sakernas_coba))

  Ringkasan desain sampling:
Stratified 2 - level Cluster Sampling design (with replacement)
With (448, 12129) clusters.
svydesign(id = ~PSU + SSU, strata = ~STRATA, data = data_coba, 
    weights = ~SUPASfinal_weightR)
Probabilities:
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0003987 0.0010163 0.0014347 0.0038716 0.0021645 0.2500000 
First-level Stratum Sizes: 
           36 01 36 02
obs         8158  3971
design.PSU   304   144
actual.PSU   304   144
Data variables:
 [1] "tahun"              "urutan"             "SUPASfinal_weightR"
 [4] "KODE_PROV"          "NAMA_PROV"          "KODE_KAB2"         
 [7] "NAMA_KAB"           "KODE_KEC2"          "NAMA_KEC"          
[10] "KLASIFIKASI2"       "no_urutNKS"         "NURT2"             
[13] "B1_R11"             "B4_K8"              "B5_R5A1"           
[16] "B5_R6"              "B5_R15A"            "B5_R15B"           
[19] "B5_R20A"            "STATUS"             "pengangguran"      
[22] "angkatankerja"      "STAT"               "PSU"               
[25] "SSU"                "STRATA"            

Konversi Satuan dan Rename Kolom

cat("================================================================\n")
cat("  POST-PROCESSING HASIL\n")
cat("================================================================\n\n")

hasilku$`pengangguran/angkatankerja`     <- round(
  hasilku$`pengangguran/angkatankerja` * 100, 4)

hasilku$`se.pengangguran/angkatankerja`  <- round(
  hasilku$`se.pengangguran/angkatankerja` * 100, 4)

hasilku$`var.pengangguran/angkatankerja` <- round(
  hasilku$`se.pengangguran/angkatankerja` *
  hasilku$`se.pengangguran/angkatankerja`, 4)

hasilku$ci_l <- round(hasilku$ci_l * 100, 4)
hasilku$ci_u <- round(hasilku$ci_u * 100, 4)
hasilku$DEff <- round(hasilku$DEff, 4)

hasilku[is.na(hasilku)] <- 0

hasilku <- hasilku %>%
  rename(
    TPT      = `pengangguran/angkatankerja`,
    SE       = `se.pengangguran/angkatankerja`,
    VAR      = `var.pengangguran/angkatankerja`,
    ci_lower = ci_l,
    ci_upper = ci_u
  )

hasilku$RSE <- ifelse(
  hasilku$TPT == 0, 0,
  round((hasilku$SE / hasilku$TPT) * 100, 4)
)
hasilku[is.na(hasilku)] <- 0

cat("  Kolom akhir : ", paste(names(hasilku), collapse = ", "), "\n\n")

================================================================
  POST-PROCESSING HASIL
================================================================

  Kolom akhir :  STAT, TPT, SE, ci_lower, ci_upper, VAR, DEff, RSE 

Ringkasan Statistik Hasil

stat_ringkasan <- data.frame(
  Statistik = c("Min", "Median", "Mean*", "Max"),
  `TPT (%)`  = c(min(hasilku$TPT),
                  median(hasilku$TPT),
                  mean(hasilku$TPT[hasilku$TPT > 0]),
                  max(hasilku$TPT)),
  `SE (%)`   = c(min(hasilku$SE),
                  median(hasilku$SE),
                  mean(hasilku$SE[hasilku$SE > 0]),
                  max(hasilku$SE)),
  `RSE (%)`  = c(min(hasilku$RSE),
                  median(hasilku$RSE),
                  mean(hasilku$RSE[hasilku$RSE > 0]),
                  max(hasilku$RSE)),
  DEff       = c(min(hasilku$DEff),
                  median(hasilku$DEff),
                  mean(hasilku$DEff[hasilku$DEff > 0]),
                  max(hasilku$DEff)),
  check.names = FALSE
)

knitr::kable(
  stat_ringkasan,
  caption = "Statistik Ringkasan Hasil Estimasi TPT per Kecamatan  (* Mean dihitung untuk nilai > 0)",
  digits  = 4,
  align   = c("l","r","r","r","r")
)

n_tidak_presisi <- sum(hasilku$RSE > 25 & hasilku$RSE != 0)
cat(sprintf("\nKecamatan RSE > 25%% (tidak presisi) = %d dari %d kecamatan\n",
            n_tidak_presisi, nrow(hasilku)))

Kecamatan RSE > 25% (tidak presisi) = 74 dari 144 kecamatan

knitr::kable(
  head(hasilku, 10),
  caption = "10 Baris Pertama Hasil Estimasi TPT",
  digits  = 4
)

Export ke Excel

cat("\n================================================================\n")
cat("  EXPORT HASIL\n")
cat("================================================================\n")

write_xlsx(hasilku, "hasilku.xlsx", col_names = TRUE)

cat("  ✓ File 'hasilku.xlsx' berhasil disimpan!\n")
cat(sprintf("  Total kecamatan            : %d\n",      nrow(hasilku)))
cat(sprintf("  Rata-rata TPT (TPT > 0)    : %.4f%%\n",  mean(hasilku$TPT[hasilku$TPT > 0])))
cat(sprintf("  Rata-rata RSE (RSE > 0)    : %.4f%%\n",  mean(hasilku$RSE[hasilku$RSE > 0])))
cat(sprintf("  Kecamatan tidak presisi    : %d\n",      n_tidak_presisi))

================================================================
  EXPORT HASIL
================================================================
  ✓ File 'hasilku.xlsx' berhasil disimpan!
  Total kecamatan            : 144
  Rata-rata TPT (TPT > 0)    : 10.3264%
  Rata-rata RSE (RSE > 0)    : 40.2107%
  Kecamatan tidak presisi    : 74

V1 — Tidak Valid (SE = Inf, Lonely PSU)

# ================================================================
# [V1] TIDAK VALID — SE = Inf
# ================================================================
cat(LINE)
cat("  [V1] TIDAK VALID — SE = Inf (Lonely PSU dalam Strata)\n")
cat(LINE)

v1 <- hasilku %>% filter(STATUS_VALIDITAS == "TIDAK VALID")
cat(sprintf("  Jumlah kecamatan: %d\n\n", nrow(v1)))

if (nrow(v1) > 0) {
  cat(sprintf("  %-14s  %-9s  %-9s  %-9s  %-6s  %-8s  %-12s\n",
              "STAT","TPT","SE","DEff","n_PSU","n_AK","Kategori_V1"))
  cat(sprintf("  %s\n", strrep("-", 75)))
  for (i in 1:min(nrow(v1), 15)) {
    kat <- ifelse(v1$n_PSU[i] == 1,
                  "1 PSU murni",
                  paste0(v1$n_PSU[i]," PSU-lonely strata"))
    cat(sprintf("  %-14s  %-9.4f  %-9s  %-9s  %-6d  %-8d  %-12s\n",
      v1$STAT[i], v1$TPT[i],
      ifelse(is.infinite(v1$SE[i]),   "Inf", sprintf("%.4f", v1$SE[i])),
      ifelse(is.infinite(v1$DEff[i]), "Inf", sprintf("%.4f", v1$DEff[i])),
      v1$n_PSU[i], v1$n_sampel_AK[i], kat))
  }
  if (nrow(v1) > 15)
    cat(sprintf("  ... dan %d kecamatan lainnya\n", nrow(v1) - 15))

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  ANALISIS DATA:\n\n")

  v1_1psu    <- v1 %>% filter(n_PSU == 1)
  v1_multpsu <- v1 %>% filter(n_PSU > 1)

  cat("  Breakdown penyebab SE = Inf:\n")
  cat(sprintf("    a) n_PSU = 1 (benar-benar hanya 1 blok sensus)  : %d kecamatan\n",
              nrow(v1_1psu)))
  cat(sprintf("    b) n_PSU > 1 tapi lonely dalam strata           : %d kecamatan\n",
              nrow(v1_multpsu)))
  cat("\n  Penjelasan kasus (b):\n")
  cat("    Kecamatan memiliki >1 PSU, NAMUN dalam strata-nya\n")
  cat("    (kombinasi KODE_PROV × KLASIFIKASI) hanya ada 1 PSU\n")
  cat("    yang ter-sampel → formula Taylor tidak bisa hitung\n")
  cat("    varians antar PSU dalam strata → SE = Inf.\n\n")

  if (nrow(v1_multpsu) > 0) {
    cat("  Detail kecamatan multi-PSU tapi SE=Inf:\n")
    cat(sprintf("  %-14s  %-6s  %-8s  %-9s\n","STAT","n_PSU","n_AK","TPT"))
    cat(sprintf("  %s\n", strrep("-", 42)))
    for (i in 1:min(nrow(v1_multpsu), 10)) {
      cat(sprintf("  %-14s  %-6d  %-8d  %-9.4f\n",
        v1_multpsu$STAT[i], v1_multpsu$n_PSU[i],
        v1_multpsu$n_sampel_AK[i], v1_multpsu$TPT[i]))
    }
  }

  tpt_v1 <- v1$TPT[is.finite(v1$TPT)]
  cat(sprintf("\n  Distribusi TPT kecamatan V1 (n=%d):\n", nrow(v1)))
  cat(sprintf("    Min = %.4f%%  |  Median = %.4f%%  |  Mean = %.4f%%  |  Max = %.4f%%\n",
    min(tpt_v1), median(tpt_v1), mean(tpt_v1), max(tpt_v1)))

  v1_risiko <- v1 %>% filter(TPT > 20) %>% arrange(desc(TPT))
  if (nrow(v1_risiko) > 0) {
    cat(sprintf("\n  ⚠ Kecamatan V1 dengan TPT > 20%% (risiko under-report tinggi):\n"))
    cat(sprintf("  %-14s  %-9s  %-6s  %-8s\n","STAT","TPT","n_PSU","n_AK"))
    cat(sprintf("  %s\n", strrep("-", 42)))
    for (i in 1:nrow(v1_risiko)) {
      cat(sprintf("  %-14s  %-9.4f  %-6d  %-8d\n",
        v1_risiko$STAT[i], v1_risiko$TPT[i],
        v1_risiko$n_PSU[i], v1_risiko$n_sampel_AK[i]))
    }
    cat("  → Kecamatan ini memiliki estimasi TPT tinggi namun\n")
    cat("    tidak dapat disajikan karena masalah desain sampling.\n")
    cat("    Potensi UNDER-COVER jika dikecualikan dari publikasi.\n")
  }

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  KESIMPULAN [V1]:\n")
  cat(sprintf("  ✗ %d kecamatan (%.1f%%) estimasi TIDAK VALID.\n",
              nrow(v1), nrow(v1)/nrow(hasilku)*100))
  cat("  ✗ Lonely PSU dalam strata → SE = Inf → tidak bisa\n")
  cat("    disajikan tanpa perbaikan desain.\n")
  cat("  → REKOMENDASI: Gunakan SAE atau collapsed strata.\n\n")
}

====================================================================
  [V1] TIDAK VALID — SE = Inf (Lonely PSU dalam Strata)
====================================================================
  Jumlah kecamatan: 0

V2 — Perlu Cek (TPT = 0 & SE = 0)

# ================================================================
# [V2] PERLU CEK — TPT = 0 & SE = 0
# ================================================================
cat(LINE)
cat("  [V2] PERLU CEK — TPT = 0 & SE = 0\n")
cat(LINE)

v2 <- hasilku %>% filter(STATUS_VALIDITAS == "PERLU CEK")
cat(sprintf("  Jumlah kecamatan: %d\n\n", nrow(v2)))

if (nrow(v2) > 0) {
  cat(sprintf("  %-14s  %-9s  %-9s  %-10s  %-10s  %-8s  %-12s\n",
              "STAT","TPT","SE","ci_lower","ci_upper","n_AK","n_Nganggur"))
  cat(sprintf("  %s\n", strrep("-", 76)))
  for (i in 1:min(nrow(v2), 15)) {
    cat(sprintf("  %-14s  %-9.4f  %-9.4f  %-10.4f  %-10.4f  %-8d  %-12d\n",
      v2$STAT[i], v2$TPT[i], v2$SE[i],
      v2$ci_lower[i], v2$ci_upper[i],
      v2$n_sampel_AK[i], v2$n_pengangguran[i]))
  }

  v2_no_ak    <- v2 %>% filter(n_sampel_AK    == 0)
  v2_no_ngang <- v2 %>% filter(n_sampel_AK    >  0 & n_pengangguran == 0)

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  ANALISIS DATA:\n\n")
  cat("  Breakdown penyebab TPT=0 & SE=0:\n")
  cat(sprintf("    a) n_AK = 0 (tidak ada AK di sampel)            : %d kecamatan\n",
              nrow(v2_no_ak)))
  cat(sprintf("    b) n_AK > 0, n_Pengangguran = 0                 : %d kecamatan\n",
              nrow(v2_no_ngang)))

  if (nrow(v2_no_ngang) > 0) {
    cat(sprintf("\n  Detail kecamatan dengan n_AK>0 & nol pengangguran:\n"))
    cat(sprintf("  %-14s  %-8s  %-14s  %-10s\n",
                "STAT","n_AK","n_Pengangguran","n_PSU"))
    cat(sprintf("  %s\n", strrep("-", 50)))
    for (i in 1:nrow(v2_no_ngang)) {
      cat(sprintf("  %-14s  %-8d  %-14d  %-10d\n",
        v2_no_ngang$STAT[i], v2_no_ngang$n_sampel_AK[i],
        v2_no_ngang$n_pengangguran[i], v2_no_ngang$n_PSU[i]))
    }
    cat(sprintf("\n  Ukuran sampel AK: Min=%d | Median=%.0f | Max=%d\n",
      min(v2_no_ngang$n_sampel_AK), median(v2_no_ngang$n_sampel_AK),
      max(v2_no_ngang$n_sampel_AK)))
    v2_besar <- v2_no_ngang %>% filter(n_sampel_AK >= 30)
    cat(sprintf("  Kec. dengan n_AK ≥ 30 (nol pengangguran lebih meyakinkan): %d\n",
                nrow(v2_besar)))
  }

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  KESIMPULAN [V2]:\n")
  if (nrow(v2_no_ak) > 0) {
    cat(sprintf("  ✗ %d kecamatan tidak ada sampel AK sama sekali\n", nrow(v2_no_ak)))
    cat("    → estimasi 0 TIDAK VALID, bukan gambaran populasi nyata.\n")
  }
  if (nrow(v2_no_ngang) > 0) {
    cat(sprintf("  ✓ %d kecamatan punya sampel AK (n_AK > 0) tapi\n", nrow(v2_no_ngang)))
    cat("    nol pengangguran teridentifikasi di sampel.\n")
    cat("  ✓ Ini BERPOTENSI VALID: mungkin memang tidak ada pengangguran.\n")
    n_kuat  <- nrow(v2_no_ngang %>% filter(n_sampel_AK >= 30))
    n_lemah <- nrow(v2_no_ngang %>% filter(n_sampel_AK < 15))
    cat(sprintf("    - %d kec. n_AK ≥ 30 → estimasi nol LEBIH DIPERCAYA\n", n_kuat))
    cat(sprintf("    - %d kec. n_AK < 15 → estimasi nol perlu konfirmasi\n", n_lemah))
  }
  cat("  → REKOMENDASI:\n")
  cat("    • Kasus (a) n_AK=0 → masuk kandidat SAE\n")
  cat("    • Kasus (b) n_AK besar → bisa disajikan dengan catatan\n")
  cat("      'tidak terdeteksi pengangguran dalam sampel'\n\n")
}

====================================================================
  [V2] PERLU CEK — TPT = 0 & SE = 0
====================================================================
  Jumlah kecamatan: 8

  STAT            TPT        SE         ci_lower    ci_upper    n_AK      n_Nganggur  
  ----------------------------------------------------------------------------
  36 01 020       0.0000     0.0000     0.0000      0.0000      14        0           
  36 01 050       0.0000     0.0000     0.0000      0.0000      15        0           
  36 02 031       0.0000     0.0000     0.0000      0.0000      14        0           
  36 02 100       0.0000     0.0000     0.0000      0.0000      34        0           
  36 02 121       0.0000     0.0000     0.0000      0.0000      19        0           
  36 03 130       0.0000     0.0000     0.0000      0.0000      23        0           
  36 04 130       0.0000     0.0000     0.0000      0.0000      49        0           
  36 71 050       0.0000     0.0000     0.0000      0.0000      34        0           

  --------------------------------------------------------------------
  ANALISIS DATA:

  Breakdown penyebab TPT=0 & SE=0:
    a) n_AK = 0 (tidak ada AK di sampel)            : 0 kecamatan
    b) n_AK > 0, n_Pengangguran = 0                 : 8 kecamatan

  Detail kecamatan dengan n_AK>0 & nol pengangguran:
  STAT            n_AK      n_Pengangguran  n_PSU     
  --------------------------------------------------
  36 01 020       14        0               1         
  36 01 050       15        0               1         
  36 02 031       14        0               1         
  36 02 100       34        0               2         
  36 02 121       19        0               1         
  36 03 130       23        0               1         
  36 04 130       49        0               3         
  36 71 050       34        0               2         

  Ukuran sampel AK: Min=14 | Median=21 | Max=49
  Kec. dengan n_AK ≥ 30 (nol pengangguran lebih meyakinkan): 3

  --------------------------------------------------------------------
  KESIMPULAN [V2]:
  ✓ 8 kecamatan punya sampel AK (n_AK > 0) tapi
    nol pengangguran teridentifikasi di sampel.
  ✓ Ini BERPOTENSI VALID: mungkin memang tidak ada pengangguran.
    - 3 kec. n_AK ≥ 30 → estimasi nol LEBIH DIPERCAYA
    - 2 kec. n_AK < 15 → estimasi nol perlu konfirmasi
  → REKOMENDASI:
    • Kasus (a) n_AK=0 → masuk kandidat SAE
    • Kasus (b) n_AK besar → bisa disajikan dengan catatan
      'tidak terdeteksi pengangguran dalam sampel'

V3 — Mencurigakan (TPT > 0 & SE = 0)

# ================================================================
# [V3] MENCURIGAKAN — TPT > 0 & SE = 0
# ================================================================
cat(LINE)
cat("  [V3] MENCURIGAKAN — TPT > 0 & SE = 0\n")
cat(LINE)

v3 <- hasilku %>% filter(STATUS_VALIDITAS == "MENCURIGAKAN")
cat(sprintf("  Jumlah kecamatan: %d\n\n", nrow(v3)))

if (nrow(v3) > 0) {
  cat(sprintf("  %-14s  %-9s  %-9s  %-9s  %-8s  %-12s  %-12s  %-15s\n",
              "STAT","TPT","SE","RSE","n_AK","n_Nganggur",
              "n_PSU_ada_ngang","max_ngang_1PSU"))
  cat(sprintf("  %s\n", strrep("-", 88)))
  for (i in 1:nrow(v3)) {
    cat(sprintf("  %-14s  %-9.4f  %-9.4f  %-9.4f  %-8d  %-12d  %-12d  %-15d\n",
      v3$STAT[i], v3$TPT[i], v3$SE[i], v3$RSE[i],
      v3$n_sampel_AK[i], v3$n_pengangguran[i],
      v3$n_PSU_dgn_ngang[i], v3$max_ngang_1psu[i]))
  }

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  ANALISIS DATA:\n\n")
  cat("  Mengapa SE = 0 padahal TPT > 0 dan n_pengangguran > 1?\n\n")
  cat("  Penjelasan teknis:\n")
  cat("  SE Taylor dihitung dari variasi ANTAR PSU dalam strata.\n")
  cat("  SE = 0 terjadi bila semua pengangguran terkonsentrasi\n")
  cat("  di PSU yang SAMA → tidak ada variasi antar PSU → SE = 0.\n\n")

  cat("  Diagnosis per kecamatan:\n")
  cat(sprintf("  %s\n", strrep("-", 68)))
  for (i in 1:nrow(v3)) {
    cat(sprintf("  Kecamatan: %s\n", v3$STAT[i]))
    cat(sprintf("    TPT = %.4f%% | n_AK = %d | n_Pengangguran = %d\n",
      v3$TPT[i], v3$n_sampel_AK[i], v3$n_pengangguran[i]))
    cat(sprintf("    PSU yang memuat pengangguran = %d dari %d PSU total\n",
      v3$n_PSU_dgn_ngang[i], v3$n_PSU[i]))
    cat(sprintf("    Maks pengangguran dalam 1 PSU = %d orang\n",
      v3$max_ngang_1psu[i]))
    if (v3$n_PSU_dgn_ngang[i] == 1 && v3$n_PSU[i] > 1) {
      cat(sprintf("    ⚠ DIAGNOSIS: Semua %d pengangguran ada di 1 PSU saja\n",
                  v3$n_pengangguran[i]))
      cat("       → PSU lain semua nol → tidak ada variasi → SE = 0\n")
    } else if (v3$n_PSU_dgn_ngang[i] == 1 && v3$n_PSU[i] == 1) {
      cat("    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung\n")
    } else {
      cat("    ⚠ DIAGNOSIS: Variasi antar PSU sangat kecil → SE mendekati 0\n")
    }
    cat(sprintf("  %s\n", strrep("-", 68)))
  }

  cat("  KESIMPULAN [V3]:\n")
  cat(sprintf("  ✗ %d kecamatan punya TPT > 0 tapi SE = 0.\n", nrow(v3)))
  cat("  ✗ Semua pengangguran terkonsentrasi di 1 PSU → Taylor SE = 0.\n")
  cat(sprintf("  ✗ TPT berkisar %.4f%% – %.4f%%: tampak besar tapi\n",
              min(v3$TPT), max(v3$TPT)))
  cat("    presisinya sama sekali tidak terukur.\n")
  cat("  → REKOMENDASI: Masuk kandidat SAE. TIDAK BOLEH disajikan.\n\n")
}

====================================================================
  [V3] MENCURIGAKAN — TPT > 0 & SE = 0
====================================================================
  Jumlah kecamatan: 30

  STAT            TPT        SE         RSE        n_AK      n_Nganggur    n_PSU_ada_ngang  max_ngang_1PSU 
  ----------------------------------------------------------------------------------------
  36 01 010       6.5541     0.0000     0.0000     17        1             1             1              
  36 01 031       20.6821    0.0000     0.0000     17        3             1             3              
  36 01 061       20.7317    0.0000     0.0000     15        3             1             3              
  36 01 070       8.7767     0.0000     0.0000     24        2             1             2              
  36 01 072       7.1978     0.0000     0.0000     22        1             1             1              
  36 01 090       4.6243     0.0000     0.0000     20        1             1             1              
  36 01 110       8.1115     0.0000     0.0000     18        1             1             1              
  36 01 111       15.1212    0.0000     0.0000     13        2             1             2              
  36 01 131       6.0776     0.0000     0.0000     32        2             1             2              
  36 01 150       22.9045    0.0000     0.0000     17        4             1             4              
  36 01 161       23.8040    0.0000     0.0000     19        4             1             4              
  36 01 170       3.7254     0.0000     0.0000     22        1             1             1              
  36 01 190       3.5866     0.0000     0.0000     20        1             1             1              
  36 01 192       8.1342     0.0000     0.0000     15        1             1             1              
  36 02 050       7.3188     0.0000     0.0000     19        1             1             1              
  36 02 051       6.6883     0.0000     0.0000     21        1             1             1              
  36 02 080       15.6825    0.0000     0.0000     27        4             1             4              
  36 02 110       4.1319     0.0000     0.0000     23        1             1             1              
  36 02 190       13.8823    0.0000     0.0000     17        2             1             2              
  36 03 060       4.7842     0.0000     0.0000     19        1             1             1              
  36 03 081       6.9307     0.0000     0.0000     15        1             1             1              
  36 03 141       45.0711    0.0000     0.0000     10        4             1             4              
  36 03 151       20.3172    0.0000     0.0000     23        5             1             5              
  36 03 160       20.7888    0.0000     0.0000     19        4             1             4              
  36 03 161       19.7013    0.0000     0.0000     16        3             1             3              
  36 04 010       11.6571    0.0000     0.0000     13        1             1             1              
  36 04 110       10.6004    0.0000     0.0000     19        2             1             2              
  36 04 180       30.4661    0.0000     0.0000     14        4             1             4              
  36 04 251       5.0839     0.0000     0.0000     17        1             1             1              
  36 04 271       71.3973    0.0000     0.0000     20        14            1             14             

  --------------------------------------------------------------------
  ANALISIS DATA:

  Mengapa SE = 0 padahal TPT > 0 dan n_pengangguran > 1?

  Penjelasan teknis:
  SE Taylor dihitung dari variasi ANTAR PSU dalam strata.
  SE = 0 terjadi bila semua pengangguran terkonsentrasi
  di PSU yang SAMA → tidak ada variasi antar PSU → SE = 0.

  Diagnosis per kecamatan:
  --------------------------------------------------------------------
  Kecamatan: 36 01 010
    TPT = 6.5541% | n_AK = 17 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 031
    TPT = 20.6821% | n_AK = 17 | n_Pengangguran = 3
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 3 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 061
    TPT = 20.7317% | n_AK = 15 | n_Pengangguran = 3
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 3 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 070
    TPT = 8.7767% | n_AK = 24 | n_Pengangguran = 2
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 2 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 072
    TPT = 7.1978% | n_AK = 22 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 090
    TPT = 4.6243% | n_AK = 20 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 110
    TPT = 8.1115% | n_AK = 18 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 111
    TPT = 15.1212% | n_AK = 13 | n_Pengangguran = 2
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 2 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 131
    TPT = 6.0776% | n_AK = 32 | n_Pengangguran = 2
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 2 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 150
    TPT = 22.9045% | n_AK = 17 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 161
    TPT = 23.8040% | n_AK = 19 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 170
    TPT = 3.7254% | n_AK = 22 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 190
    TPT = 3.5866% | n_AK = 20 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 01 192
    TPT = 8.1342% | n_AK = 15 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 02 050
    TPT = 7.3188% | n_AK = 19 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 02 051
    TPT = 6.6883% | n_AK = 21 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 02 080
    TPT = 15.6825% | n_AK = 27 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 02 110
    TPT = 4.1319% | n_AK = 23 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 02 190
    TPT = 13.8823% | n_AK = 17 | n_Pengangguran = 2
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 2 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 060
    TPT = 4.7842% | n_AK = 19 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 081
    TPT = 6.9307% | n_AK = 15 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 141
    TPT = 45.0711% | n_AK = 10 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 151
    TPT = 20.3172% | n_AK = 23 | n_Pengangguran = 5
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 5 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 160
    TPT = 20.7888% | n_AK = 19 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 03 161
    TPT = 19.7013% | n_AK = 16 | n_Pengangguran = 3
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 3 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 04 010
    TPT = 11.6571% | n_AK = 13 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 04 110
    TPT = 10.6004% | n_AK = 19 | n_Pengangguran = 2
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 2 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 04 180
    TPT = 30.4661% | n_AK = 14 | n_Pengangguran = 4
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 4 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 04 251
    TPT = 5.0839% | n_AK = 17 | n_Pengangguran = 1
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 1 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  Kecamatan: 36 04 271
    TPT = 71.3973% | n_AK = 20 | n_Pengangguran = 14
    PSU yang memuat pengangguran = 1 dari 1 PSU total
    Maks pengangguran dalam 1 PSU = 14 orang
    ⚠ DIAGNOSIS: Hanya ada 1 PSU total → SE tidak bisa dihitung
  --------------------------------------------------------------------
  KESIMPULAN [V3]:
  ✗ 30 kecamatan punya TPT > 0 tapi SE = 0.
  ✗ Semua pengangguran terkonsentrasi di 1 PSU → Taylor SE = 0.
  ✗ TPT berkisar 3.5866% – 71.3973%: tampak besar tapi
    presisinya sama sekali tidak terukur.
  → REKOMENDASI: Masuk kandidat SAE. TIDAK BOLEH disajikan.

V4 — Perlu Perhatian (CI Lower < 0)

# ================================================================
# [V4] PERLU PERHATIAN — CI Lower < 0
# ================================================================
cat(LINE)
cat("  [V4] PERLU PERHATIAN — CI Lower < 0\n")
cat(LINE)

v4 <- hasilku %>% filter(STATUS_VALIDITAS == "PERLU PERHATIAN")
cat(sprintf("  Jumlah kecamatan: %d\n", nrow(v4)))
if (nrow(v4) == 0) {
  cat("  → Tidak ada kecamatan dalam kategori ini.\n")
  cat("  ✓ Semua CI lower dari estimasi yang lolos V1-V3 bernilai >= 0.\n\n")
} else {
  cat("\n")
}

====================================================================
  [V4] PERLU PERHATIAN — CI Lower < 0
====================================================================
  Jumlah kecamatan: 0
  → Tidak ada kecamatan dalam kategori ini.
  ✓ Semua CI lower dari estimasi yang lolos V1-V3 bernilai >= 0.

V5 — Tidak Presisi (RSE > 25%)

# ================================================================
# [V5] TIDAK PRESISI — RSE > 25%
# ================================================================
cat(LINE)
cat("  [V5] TIDAK PRESISI — RSE > 25%\n")
cat(LINE)

v5 <- hasilku %>% filter(STATUS_VALIDITAS == "TIDAK PRESISI")
cat(sprintf("  Jumlah kecamatan: %d\n\n", nrow(v5)))

if (nrow(v5) > 0) {
  v5_sorted <- v5 %>% arrange(desc(RSE))

  cat(sprintf("  %-14s  %-9s  %-9s  %-9s  %-9s  %-8s\n",
              "STAT","TPT","SE","RSE","DEff","n_AK"))
  cat(sprintf("  %s\n", strrep("-", 62)))
  for (i in 1:min(nrow(v5_sorted), 15)) {
    cat(sprintf("  %-14s  %-9.4f  %-9.4f  %-9.4f  %-9.4f  %-8d\n",
      v5_sorted$STAT[i], v5_sorted$TPT[i], v5_sorted$SE[i],
      v5_sorted$RSE[i], v5_sorted$DEff[i], v5_sorted$n_sampel_AK[i]))
  }
  if (nrow(v5) > 15)
    cat(sprintf("  ... dan %d kecamatan lainnya\n", nrow(v5) - 15))

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  ANALISIS DATA:\n\n")

  z1 <- sum(v5$RSE > 25 & v5$RSE <= 35)
  z2 <- sum(v5$RSE > 35 & v5$RSE <= 50)
  z3 <- sum(v5$RSE > 50)
  cat("  Distribusi zona RSE:\n")
  cat(sprintf("    25%% < RSE ≤ 35%%  : %2d kecamatan — tidak terlalu buruk\n",z1))
  cat(sprintf("    35%% < RSE ≤ 50%%  : %2d kecamatan — buruk\n",              z2))
  cat(sprintf("    RSE > 50%%         : %2d kecamatan — sangat buruk\n\n",      z3))

  grp1 <- v5$RSE[v5$n_sampel_AK < 10]
  grp2 <- v5$RSE[v5$n_sampel_AK >= 10 & v5$n_sampel_AK <= 30]
  grp3 <- v5$RSE[v5$n_sampel_AK > 30]
  cat("  Hubungan ukuran sampel AK vs RSE:\n")
  cat(sprintf("    n_AK < 10   (n=%2d) → median RSE = %s\n",
    length(grp1), ifelse(length(grp1)==0, "N/A (tidak ada kec.)",
                         sprintf("%.2f%%", median(grp1)))))
  cat(sprintf("    n_AK 10–30  (n=%2d) → median RSE = %s\n",
    length(grp2), ifelse(length(grp2)==0, "N/A (tidak ada kec.)",
                         sprintf("%.2f%%", median(grp2)))))
  cat(sprintf("    n_AK > 30   (n=%2d) → median RSE = %s\n",
    length(grp3), ifelse(length(grp3)==0, "N/A (tidak ada kec.)",
                         sprintf("%.2f%%", median(grp3)))))
  cat("  → Semakin kecil n_AK, semakin besar RSE (konfirmasi teori).\n\n")

  cat("  Hubungan nilai TPT vs RSE:\n")
  t1 <- v5$RSE[v5$TPT < 5]
  t2 <- v5$RSE[v5$TPT >= 5 & v5$TPT < 10]
  t3 <- v5$RSE[v5$TPT >= 10]
  cat(sprintf("    TPT < 5%%   (n=%2d) → median RSE = %s\n",
    length(t1), ifelse(length(t1)==0,"N/A",sprintf("%.2f%%",median(t1)))))
  cat(sprintf("    TPT 5–10%%  (n=%2d) → median RSE = %s\n",
    length(t2), ifelse(length(t2)==0,"N/A",sprintf("%.2f%%",median(t2)))))
  cat(sprintf("    TPT ≥ 10%%  (n=%2d) → median RSE = %s\n",
    length(t3), ifelse(length(t3)==0,"N/A",sprintf("%.2f%%",median(t3)))))
  cat("  → Kecamatan dengan TPT rendah cenderung RSE lebih tinggi\n")
  cat("    karena pembilang (SE) relatif besar terhadap TPT kecil.\n\n")

  v5_deff <- v5 %>% filter(DEff > 2)
  cat(sprintf("  Kecamatan DEff > 2 (efisiensi desain rendah): %d\n",
              nrow(v5_deff)))
  if (nrow(v5_deff) > 0) {
    for (i in 1:nrow(v5_deff)) {
      cat(sprintf("    %-14s → DEff=%.4f | RSE=%.2f%% | n_AK=%d\n",
        v5_deff$STAT[i], v5_deff$DEff[i],
        v5_deff$RSE[i], v5_deff$n_sampel_AK[i]))
    }
  }

  cat(sprintf("\n  %s\n", strrep("-", 68)))
  cat("  KESIMPULAN [V5]:\n")
  cat(sprintf("  ✗ %d kecamatan (%.1f%%) RSE > 25%% → tidak presisi (standar BPS).\n",
              nrow(v5), nrow(v5)/nrow(hasilku)*100))
  cat(sprintf("  ✗ %d kec. RSE > 50%% (sangat buruk) → kandidat kuat SAE.\n", z3))
  cat(sprintf("  ✗ %d kec. RSE 25–35%% → bisa disajikan dengan catatan.\n", z1))
  cat("  → REKOMENDASI:\n")
  cat("    • RSE 25–35%% → beri tanda '*' di tabel publikasi\n")
  cat("    • RSE 35–50%% → beri tanda '**', pertimbangkan SAE\n")
  cat("    • RSE > 50%%   → wajib SAE, jangan disajikan langsung\n\n")
}