AED UAS

# =============================================================================
# EDA LENGKAP: Analisis Salah Sasaran Bansos — Jawa Barat, SUSENAS Maret 2023
# Tujuan analisis:
#   1. Mengidentifikasi inclusion error & exclusion error bansos
#   2. Memotret kelas menengah terjepit (di atas GK, tapi ekonomi tidak aman)
#   3. Membandingkan pola perkotaan vs perdesaan
#   4. Eksplorasi profil demografis KRT terkait salah sasaran
#
# File yang dibutuhkan (letakkan di folder kerja):
#   - 2023 Maret JABAR - SUSENAS KOR Rumah Tangga.csv
#   - 2023 Maret JABAR - SUSENAS KOR INDIVIDU PART1.csv
#   - 2023 Maret JABAR - SUSENAS KP BP 4.3.csv
#
# Catatan data:
#   - KAPITA = pengeluaran per kapita per BULAN (Rp, bukan ribuan)
#   - Garis Kemiskinan Jabar Maret 2023: Rp 479.934/kapita/bulan (BPS)
#   - Kode bansos: 1=aktif, 2=pernah tapi tidak lagi, 5=tidak pernah
#   - Kode aset:   1=punya, 5=tidak punya
#   - R403 == 1   → Kepala Rumah Tangga (KRT)
# =============================================================================

Warna

# =============================================================================
# 0. SETUP & LIBRARY
# =============================================================================

# install.packages(c("tidyverse","janitor","scales","patchwork","ggridges","corrplot","gt"))

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.5.3

## Warning: package 'ggplot2' was built under R version 4.5.2

## Warning: package 'tidyr' was built under R version 4.5.2

## Warning: package 'readr' was built under R version 4.5.2

## Warning: package 'purrr' was built under R version 4.5.2

## Warning: package 'dplyr' was built under R version 4.5.3

## Warning: package 'forcats' was built under R version 4.5.2

## Warning: package 'lubridate' was built under R version 4.5.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.5.1
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(janitor)

## Warning: package 'janitor' was built under R version 4.5.3

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(scales)

## Warning: package 'scales' was built under R version 4.5.3

## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.5.3

library(ggridges)

## Warning: package 'ggridges' was built under R version 4.5.2

library(corrplot)

## Warning: package 'corrplot' was built under R version 4.5.3

## corrplot 0.95 loaded

library(gt)

## Warning: package 'gt' was built under R version 4.5.3

# --- Tema ggplot global -------------------------------------------------------
theme_set(
  theme_minimal(base_size = 12) +
    theme(
      plot.title    = element_text(face = "bold", size = 14),
      plot.subtitle = element_text(color = "grey40", size = 11),
      plot.caption  = element_text(color = "grey55", size = 9),
      legend.position = "bottom",
      strip.text    = element_text(face = "bold")
    )
)

# --- Palet warna --------------------------------------------------------------
warna_wilayah <- c("Perkotaan" = "#2563EB", "Perdesaan" = "#16A34A")
warna_sasaran <- c(
  "Tepat sasaran"                         = "#10ce12",
  "Exclusion error (miskin, tidak dapat)" = "#1324aa",
  "Inclusion error (mampu, tapi dapat)"   = "#ef4c5b",
  "Bukan sasaran"                         = "#94A3B8"
)

# --- Garis kemiskinan BPS Jabar Maret 2023 ------------------------------------
GARIS_KEMISKINAN <- 479934
GK_JT            <- GARIS_KEMISKINAN / 1e6

# =============================================================================
# 1. LOAD DATA
# =============================================================================

kor_rt <- read_csv(
  "C:/Users/hp/Downloads/2023 Maret JABAR - SUSENAS KOR Rumah Tangga.csv",
  show_col_types = FALSE
) |> clean_names()

## New names:
## • `` -> `...1`

kor_ind <- read_csv(
  "C:/Users/hp/Downloads/2023 Maret JABAR - SUSENAS KOR INDIVIDU PART1.csv",
  show_col_types = FALSE,
  col_types = cols(.default = "c")     # baca semua char dulu (mixed types)
) |>
  clean_names() |>
  mutate(across(c(urut, r401, r403, r405, r407, r614, r707, r704), as.numeric))

## New names:
## • `` -> `...1`

kp_b43 <- read_csv(
  "C:/Users/hp/Downloads/2023 Maret JABAR - SUSENAS KP BP 4.3 (1).csv",
  show_col_types = FALSE
) |> clean_names()

## New names:
## • `` -> `...1`

cat("=== Dimensi data ===\n")

## === Dimensi data ===

cat("  KOR RT  :", nrow(kor_rt),  "RT  x", ncol(kor_rt),  "kolom\n")

##   KOR RT  : 25890 RT  x 199 kolom

cat("  KOR IND :", nrow(kor_ind), "ind x", ncol(kor_ind), "kolom\n")

##   KOR IND : 84688 ind x 183 kolom

cat("  KP B43  :", nrow(kp_b43),  "RT  x", ncol(kp_b43),  "kolom\n")

##   KP B43  : 25890 RT  x 20 kolom

# =============================================================================
# 2. SELEKSI VARIABEL, REKODE, & MERGE
# =============================================================================

# --- 2.1 KOR RT ---------------------------------------------------------------
kor_rt_sel <- kor_rt |>
  select(
    urut, r101, r102, r105, r301, fwt,
    # Bansos
    r2202, r2203, r2204a, r2207, r2208a2, r2209a, r2209b,
    # Perumahan
    r1802, r1804, r1806, r1807, r1808, r1809a, r1810a, r1816,
    # Aset
    r2001b, r2001c, r2001f, r2001h, r2001k
  )

# --- 2.2 KP Blok 4.3 — pengeluaran per kapita ---------------------------------
kp_sel <- kp_b43 |>
  select(urut, food, nonfood, expend, kapita, kalori_kap, prote_kap)

# --- 2.3 KOR Individu — ambil KRT (r403 == 1) ---------------------------------
kor_krt <- kor_ind |>
  filter(r403 == 1) |>
  select(urut, r405, r407, r614, r707, r704) |>
  rename(
    jk_krt     = r405,   # jenis kelamin KRT
    umur_krt   = r407,   # umur KRT
    pend_krt   = r614,   # pendidikan tertinggi KRT
    kerja_krt  = r707,   # kedudukan pekerjaan
    sektor_krt = r704    # lapangan usaha
  )

cat("\nJumlah KRT terambil:", nrow(kor_krt), "\n")

## 
## Jumlah KRT terambil: 25890

# --- 2.4 Merge ----------------------------------------------------------------
df_raw <- kor_rt_sel |>
  left_join(kp_sel,  by = "urut") |>
  left_join(kor_krt, by = "urut")

cat("Jumlah RT setelah merge:", nrow(df_raw), "\n")

## Jumlah RT setelah merge: 25890

cat("Missing KAPITA   :", sum(is.na(df_raw$kapita)),   "\n")

## Missing KAPITA   : 0

cat("Missing pend_krt :", sum(is.na(df_raw$pend_krt)), "\n")

## Missing pend_krt : 0

# =============================================================================
# 3. KONSTRUKSI VARIABEL ANALITIK
# =============================================================================

label_kab <- c(
  "1"="Bogor","2"="Sukabumi","3"="Cianjur","4"="Bandung",
  "5"="Garut","6"="Tasikmalaya","7"="Ciamis","8"="Kuningan",
  "9"="Cirebon","10"="Majalengka","11"="Sumedang","12"="Indramayu",
  "13"="Subang","14"="Purwakarta","15"="Karawang","16"="Bekasi",
  "17"="Bandung Barat","18"="Pangandaran",
  "71"="Kota Bogor","72"="Kota Sukabumi","73"="Kota Bandung",
  "74"="Kota Cirebon","75"="Kota Bekasi","76"="Kota Depok",
  "77"="Kota Cimahi","78"="Kota Tasikmalaya","79"="Kota Banjar"
)

df <- df_raw |>
  mutate(

    # ---- WILAYAH ----
    wilayah   = factor(r105, levels = c(1,2), labels = c("Perkotaan","Perdesaan")),
    nama_kab  = label_kab[as.character(r102)],
    jenis_kab = if_else(r102 >= 71, "Kota", "Kabupaten"),

    # ---- BANSOS (kode 1 = aktif) ----
    penerima_pkh      = if_else(r2204a == 1, 1L, 0L, missing = 0L),
    penerima_bpnt     = if_else(r2207  == 1, 1L, 0L, missing = 0L),
    penerima_kks      = if_else(r2202  == 1, 1L, 0L, missing = 0L),
    penerima_blt_bbm  = if_else(r2209a == 1, 1L, 0L, missing = 0L),
    penerima_blt_desa = if_else(r2209b == 1, 1L, 0L, missing = 0L),

    penerima  = if_else(
      penerima_pkh == 1 | penerima_bpnt == 1 | penerima_kks == 1 |
        penerima_blt_bbm == 1 | penerima_blt_desa == 1, 1L, 0L),
    n_program = penerima_pkh + penerima_bpnt + penerima_kks +
      penerima_blt_bbm + penerima_blt_desa,

    # ---- KEMISKINAN & DESIL ----
    layak       = if_else(kapita <= GARIS_KEMISKINAN, 1L, 0L, missing = NA_integer_),
    desil       = ntile(kapita, 10),
    kapita_juta = kapita / 1e6,

    kelompok_ekono = case_when(
      desil %in% 1:2  ~ "Desil 1-2\n(Sangat miskin)",
      desil %in% 3:4  ~ "Desil 3-4\n(Miskin)",
      desil %in% 5:6  ~ "Desil 5-6\n(Rentan)",
      desil %in% 7:8  ~ "Desil 7-8\n(Menengah bawah)",
      desil %in% 9:10 ~ "Desil 9-10\n(Menengah atas)"
    ) |> factor(levels = c(
      "Desil 1-2\n(Sangat miskin)","Desil 3-4\n(Miskin)",
      "Desil 5-6\n(Rentan)","Desil 7-8\n(Menengah bawah)",
      "Desil 9-10\n(Menengah atas)"
    )),

    # ---- JENIS KESALAHAN SASARAN ----
    jenis_sasaran = case_when(
      layak == 1 & penerima == 1 ~ "Tepat sasaran",
      layak == 1 & penerima == 0 ~ "Exclusion error (miskin, tidak dapat)",
      layak == 0 & penerima == 1 ~ "Inclusion error (mampu, tapi dapat)",
      layak == 0 & penerima == 0 ~ "Bukan sasaran",
      TRUE ~ NA_character_
    ) |> factor(levels = names(warna_sasaran)),

    # ---- ASET (kode 1=punya, 5=tidak) ----
    punya_kulkas = if_else(r2001b == 1, 1L, 0L, missing = 0L),
    punya_ac     = if_else(r2001c == 1, 1L, 0L, missing = 0L),
    punya_laptop = if_else(r2001f == 1, 1L, 0L, missing = 0L),
    punya_motor  = if_else(r2001h == 1, 1L, 0L, missing = 0L),
    punya_mobil  = if_else(r2001k == 1, 1L, 0L, missing = 0L),
    skor_aset    = punya_kulkas + punya_ac + punya_laptop + punya_motor + punya_mobil,

    # ---- KONDISI RUMAH ----
    listrik_pln  = if_else(r1816  == 1,         1L, 0L, missing = 0L),
    rumah_milik  = if_else(r1802  == 1,         1L, 0L, missing = 0L),
    lantai_layak = if_else(r1808  != 1,         1L, 0L, missing = 0L),
    punya_bab    = if_else(r1809a == 1,         1L, 0L, missing = 0L),
    air_layak    = if_else(r1810a %in% c(1,2,3),1L, 0L, missing = 0L),

    # ---- DEMOGRAFI KRT ----
    jk_krt_label = factor(jk_krt, levels = c(1,2),
                           labels = c("Laki-laki","Perempuan")),

    kelompok_umur = case_when(
      umur_krt < 30              ~ "< 30 tahun",
      umur_krt < 45              ~ "30-44 tahun",
      umur_krt < 60              ~ "45-59 tahun",
      TRUE                       ~ "≥ 60 tahun"
    ) |> factor(levels = c("< 30 tahun","30-44 tahun","45-59 tahun","≥ 60 tahun")),

    pendidikan_krt = case_when(
      pend_krt %in% 0:2   ~ "Tidak/Belum tamat SD",
      pend_krt %in% 3:7   ~ "SD",
      pend_krt %in% 8:12  ~ "SMP",
      pend_krt %in% 13:17 ~ "SMA/SMK",
      pend_krt %in% 18:25 ~ "Diploma/Sarjana+",
      TRUE ~ NA_character_
    ) |> factor(levels = c(
      "Tidak/Belum tamat SD","SD","SMP","SMA/SMK","Diploma/Sarjana+"
    )),

    status_kerja_krt = case_when(
      kerja_krt == 0          ~ "Tidak bekerja",
      kerja_krt == 1          ~ "Berusaha sendiri",
      kerja_krt %in% c(2,3)  ~ "Berusaha (punya buruh)",
      kerja_krt == 4          ~ "Buruh/Karyawan",
      kerja_krt == 5          ~ "Pekerja bebas",
      kerja_krt == 6          ~ "Pekerja keluarga",
      TRUE ~ NA_character_
    ) |> factor(levels = c(
      "Tidak bekerja","Berusaha sendiri","Berusaha (punya buruh)",
      "Buruh/Karyawan","Pekerja bebas","Pekerja keluarga"
    )),

    sektor_krt_label = case_when(
      sektor_krt == 0 ~ "Tidak bekerja",
      sektor_krt == 1 ~ "Pertanian",
      sektor_krt == 2 ~ "Industri",
      sektor_krt == 3 ~ "Perdagangan/Jasa",
      sektor_krt == 4 ~ "Lainnya",
      TRUE ~ NA_character_
    ) |> factor(levels = c(
      "Tidak bekerja","Pertanian","Industri","Perdagangan/Jasa","Lainnya"
    )),

    # ---- KELAS MENENGAH TERJEPIT ----
    # Definisi: di atas GK, tapi masih desil < 9, dan tidak dapat bansos apapun
    kejepit = if_else(layak == 0 & desil < 9 & penerima == 0, 1L, 0L, missing = 0L)
  )

cat("\n=== Distribusi Wilayah ===\n");        print(table(df$wilayah))

## 
## === Distribusi Wilayah ===

## 
## Perkotaan Perdesaan 
##     16953      8937

cat("\n=== Distribusi Jenis Sasaran ===\n");   print(table(df$jenis_sasaran, useNA="ifany"))

## 
## === Distribusi Jenis Sasaran ===

## 
##                         Tepat sasaran Exclusion error (miskin, tidak dapat) 
##                                   621                                   512 
##   Inclusion error (mampu, tapi dapat)                         Bukan sasaran 
##                                  9048                                 15709

cat("\n=== % Miskin Jabar (GK Rp 479.934) ===\n")

## 
## === % Miskin Jabar (GK Rp 479.934) ===

cat(percent(mean(df$layak, na.rm=TRUE), accuracy=0.01), "\n")

## 4.38%

# =============================================================================
# 4. STATISTIK DESKRIPTIF — TABEL RINGKASAN
# =============================================================================

ringkasan <- df |>
  filter(!is.na(wilayah)) |>
  group_by(wilayah) |>
  summarise(
    n_rt              = n(),
    pct_penerima      = mean(penerima,   na.rm=TRUE),
    pct_layak         = mean(layak,      na.rm=TRUE),
    pct_tepat         = mean(jenis_sasaran == "Tepat sasaran",                         na.rm=TRUE),
    pct_exclusion     = mean(jenis_sasaran == "Exclusion error (miskin, tidak dapat)", na.rm=TRUE),
    pct_inclusion     = mean(jenis_sasaran == "Inclusion error (mampu, tapi dapat)",   na.rm=TRUE),
    pct_kejepit       = mean(kejepit,    na.rm=TRUE),
    median_kapita     = median(kapita,   na.rm=TRUE),
    median_umur_krt   = median(umur_krt, na.rm=TRUE),
    pct_krt_perempuan = mean(jk_krt == 2, na.rm=TRUE),
    .groups = "drop"
  )

ringkasan |>
  gt() |>
  tab_header(
    title    = "Ringkasan EDA Bansos Jawa Barat — SUSENAS Maret 2023",
    subtitle = "Perbandingan Perkotaan vs Perdesaan"
  ) |>
  fmt_number(columns  = c(n_rt, median_kapita, median_umur_krt), decimals=0, use_seps=TRUE) |>
  fmt_percent(columns = starts_with("pct_"), decimals=1) |>
  cols_label(
    wilayah           = "Wilayah",
    n_rt              = "Jumlah RT",
    pct_penerima      = "% Penerima",
    pct_layak         = "% Miskin",
    pct_tepat         = "% Tepat sasaran",
    pct_exclusion     = "% Exclusion error",
    pct_inclusion     = "% Inclusion error",
    pct_kejepit       = "% Kelas menengah terjepit",
    median_kapita     = "Median kapita (Rp/bln)",
    median_umur_krt   = "Median umur KRT (thn)",
    pct_krt_perempuan = "% KRT perempuan"
  )

Ringkasan EDA Bansos Jawa Barat — SUSENAS Maret 2023
Perbandingan Perkotaan vs Perdesaan
Wilayah	Jumlah RT	% Penerima	% Miskin	% Tepat sasaran	% Exclusion error	% Inclusion error	% Kelas menengah terjepit	Median kapita (Rp/bln)	Median umur KRT (thn)	% KRT perempuan
Perkotaan	16,953	31.6%	4.1%	2.1%	2.0%	29.6%	44.5%	1,303,901	49	15.7%
Perdesaan	8,937	48.2%	5.0%	3.0%	2.0%	45.1%	42.2%	1,085,099	50	14.4%

# =============================================================================
# 5. VISUALISASI — PLOT DISTRIBUSI & KOMPOSISI
# =============================================================================

# --------------------------------------------------------------------------
# PLOT 1: Density distribusi pengeluaran per kapita
# --------------------------------------------------------------------------
p1 <- df |>
  filter(!is.na(wilayah), !is.na(penerima)) |>
  mutate(status = if_else(penerima == 1, "Penerima bansos","Bukan penerima")) |>
  ggplot(aes(x=kapita_juta, fill=wilayah, color=wilayah)) +
  geom_density(alpha=0.3, linewidth=0.8) +
  geom_vline(xintercept=GK_JT, linetype="dashed", color="firebrick", linewidth=0.9) +
  annotate("text", x=GK_JT+0.02, y=Inf, vjust=1.5, hjust=0,
           label="GK Jabar\nRp 480 rb", color="firebrick", size=2.8) +
  facet_wrap(~status, ncol=1) +
  scale_x_continuous(limits=c(0,8), labels=label_number(suffix=" jt", accuracy=0.1)) +
  scale_fill_manual(values=warna_wilayah) +
  scale_color_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 1: Distribusi Pengeluaran Per Kapita",
    subtitle = "Penerima vs bukan penerima bansos",
    x="Pengeluaran per kapita per bulan (juta Rp)", y="Kepadatan",
    fill=NULL, color=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p1)

## Warning: Removed 515 rows containing non-finite outside the scale range
## (`stat_density()`).

# --------------------------------------------------------------------------
# PLOT 2: % Penerima bansos per desil
# --------------------------------------------------------------------------
p2 <- df |>
  filter(!is.na(wilayah), !is.na(desil)) |>
  group_by(wilayah, desil) |>
  summarise(pct=mean(penerima, na.rm=TRUE), .groups="drop") |>
  ggplot(aes(x=desil, y=pct, color=wilayah, group=wilayah)) +
  geom_line(linewidth=1.3) + geom_point(size=3.5) +
  geom_vline(xintercept=4.5, linetype="dotted", color="grey50") +
  annotate("text", x=3.5, y=0.97, label="← Miskin",       color="grey40", size=3.2, fontface="italic") +
  annotate("text", x=5.5, y=0.97, label="Tidak miskin →", color="grey40", size=3.2, fontface="italic") +
  scale_x_continuous(breaks=1:10, labels=paste0("D",1:10)) +
  scale_y_continuous(labels=percent_format(accuracy=1), limits=c(0,1)) +
  scale_color_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 2: % Penerima Bansos per Desil Pengeluaran",
    subtitle = "Desil 1 = termiskin | Garis putus = batas miskin",
    x="Desil pengeluaran per kapita", y="Proporsi penerima",
    color=NULL, caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p2)

# --------------------------------------------------------------------------
# PLOT 3: Stacked bar — komposisi jenis sasaran per wilayah
# --------------------------------------------------------------------------
p3 <- df |>
  filter(!is.na(wilayah), !is.na(jenis_sasaran)) |>
  count(wilayah, jenis_sasaran) |>
  group_by(wilayah) |>
  mutate(pct=n/sum(n)) |>
  ggplot(aes(x=wilayah, y=pct, fill=jenis_sasaran)) +
  geom_col(width=0.6, color="white", linewidth=0.4) +
  geom_text(aes(label=if_else(pct>0.03, percent(pct, accuracy=0.1), "")),
            position=position_stack(vjust=0.5),
            size=3.5, color="white", fontface="bold") +
  scale_y_continuous(labels=percent_format()) +
  scale_fill_manual(values=warna_sasaran) +
  labs(
    title    = "Plot 3: Komposisi Ketepatan Sasaran Bansos",
    subtitle = "Perkotaan vs Perdesaan",
    x=NULL, y="Proporsi rumah tangga", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p3)

# --------------------------------------------------------------------------
# PLOT 4: Grouped bar — inclusion vs exclusion error
# --------------------------------------------------------------------------
p4 <- df |>
  filter(!is.na(wilayah), !is.na(layak)) |>
  group_by(wilayah) |>
  summarise(
    `Inclusion error\n(mampu, dapat bansos)` = mean(layak==0 & penerima==1, na.rm=TRUE),
    `Exclusion error\n(miskin, tidak dapat)` = mean(layak==1 & penerima==0, na.rm=TRUE),
    .groups="drop"
  ) |>
  pivot_longer(-wilayah, names_to="jenis_error", values_to="pct") |>
  ggplot(aes(x=jenis_error, y=pct, fill=wilayah)) +
  geom_col(position="dodge", width=0.6) +
  geom_text(aes(label=percent(pct, accuracy=0.1)),
            position=position_dodge(width=0.6),
            vjust=-0.5, size=3.8, fontface="bold", color="grey30") +
  scale_y_continuous(labels=percent_format(), expand=expansion(mult=c(0,0.2))) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 4: Tingkat Inclusion Error vs Exclusion Error",
    subtitle = "Perbandingan perkotaan dan perdesaan",
    x=NULL, y="Proporsi rumah tangga", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p4)

# --------------------------------------------------------------------------
# PLOT 5: Ridge plot — distribusi kapita per jenis sasaran
# --------------------------------------------------------------------------
p5 <- df |>
  filter(!is.na(jenis_sasaran), !is.na(wilayah)) |>
  ggplot(aes(x=kapita_juta, y=jenis_sasaran, fill=wilayah, color=wilayah)) +
  geom_density_ridges(alpha=0.35, scale=0.85, rel_min_height=0.01) +
  geom_vline(xintercept=GK_JT, linetype="dashed", color="firebrick", linewidth=0.8) +
  facet_wrap(~wilayah) +
  scale_x_continuous(limits=c(0,5), labels=label_number(suffix=" jt", accuracy=0.1)) +
  scale_fill_manual(values=warna_wilayah) +
  scale_color_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 5: Distribusi Kapita berdasarkan Status Sasaran",
    subtitle = "Garis merah = garis kemiskinan (Rp 479.934/kapita/bln)",
    x="Pengeluaran per kapita (juta Rp/bln)", y=NULL,
    fill=NULL, color=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  ) + theme(legend.position="none")
print(p5)

## Picking joint bandwidth of 0.0649

## Picking joint bandwidth of 0.0597

## Warning: Removed 1481 rows containing non-finite outside the scale range
## (`stat_density_ridges()`).

# --------------------------------------------------------------------------
# PLOT 6: Cakupan per program bansos per wilayah
# --------------------------------------------------------------------------
p6 <- df |>
  filter(!is.na(wilayah)) |>
  group_by(wilayah) |>
  summarise(
    PKH        = mean(penerima_pkh,      na.rm=TRUE),
    BPNT       = mean(penerima_bpnt,     na.rm=TRUE),
    KKS        = mean(penerima_kks,      na.rm=TRUE),
    `BLT BBM`  = mean(penerima_blt_bbm,  na.rm=TRUE),
    `BLT Desa` = mean(penerima_blt_desa, na.rm=TRUE),
    .groups="drop"
  ) |>
  pivot_longer(-wilayah, names_to="program", values_to="pct") |>
  ggplot(aes(x=reorder(program, pct), y=pct, fill=wilayah)) +
  geom_col(position="dodge", width=0.65) +
  geom_text(aes(label=percent(pct, accuracy=0.1)),
            position=position_dodge(width=0.65),
            hjust=-0.1, size=3.2, color="grey30") +
  coord_flip() +
  scale_y_continuous(labels=percent_format(), limits=c(0,0.40),
                     expand=expansion(mult=c(0,0.05))) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 6: Cakupan per Program Bansos",
    subtitle = "% RT penerima aktif per program",
    x=NULL, y="% penerima", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p6)

# --------------------------------------------------------------------------
# PLOT 7: Cakupan program per desil (kebocoran ke desil atas)
# --------------------------------------------------------------------------
p7 <- df |>
  filter(!is.na(wilayah), !is.na(desil)) |>
  group_by(wilayah, desil) |>
  summarise(
    PKH  = mean(penerima_pkh,  na.rm=TRUE),
    BPNT = mean(penerima_bpnt, na.rm=TRUE),
    BLT  = mean(penerima_blt_bbm + penerima_blt_desa > 0, na.rm=TRUE),
    .groups="drop"
  ) |>
  pivot_longer(c(PKH,BPNT,BLT), names_to="program", values_to="pct") |>
  ggplot(aes(x=desil, y=pct, color=wilayah)) +
  geom_line(linewidth=1) + geom_point(size=2) +
  facet_wrap(~program, ncol=3) +
  scale_x_continuous(breaks=1:10, labels=paste0("D",1:10)) +
  scale_y_continuous(labels=percent_format(accuracy=1)) +
  scale_color_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 7: Cakupan Program per Desil",
    subtitle = "Seberapa banyak desil atas yang masih menerima tiap program?",
    x="Desil pengeluaran", y="% penerima",
    color=NULL, caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p7)

# --------------------------------------------------------------------------
# PLOT 8: Kelas menengah terjepit per desil
# --------------------------------------------------------------------------
p8 <- df |>
  filter(!is.na(wilayah), !is.na(desil)) |>
  group_by(wilayah, desil) |>
  summarise(pct_kejepit=mean(kejepit, na.rm=TRUE), .groups="drop") |>
  ggplot(aes(x=desil, y=pct_kejepit, fill=wilayah)) +
  geom_col(position="dodge", width=0.7) +
  scale_x_continuous(breaks=1:10, labels=paste0("D",1:10)) +
  scale_y_continuous(labels=percent_format(accuracy=1),
                     expand=expansion(mult=c(0,0.1))) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 8: Proporsi Kelas Menengah Terjepit per Desil",
    subtitle = "Di atas GK, desil < 9, tidak menerima bansos apapun",
    x="Desil pengeluaran per kapita", y="% rumah tangga", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p8)

# --------------------------------------------------------------------------
# PLOT 9: Heatmap korelasi
# --------------------------------------------------------------------------
vars_cor <- df |>
  select(kapita, skor_aset, r1804,
         punya_motor, punya_mobil, punya_kulkas, punya_ac,
         listrik_pln, lantai_layak, punya_bab,
         penerima, layak) |>
  mutate(across(everything(), as.numeric)) |>
  drop_na()

mat_cor <- cor(vars_cor, use="complete.obs")
colnames(mat_cor) <- rownames(mat_cor) <- c(
  "Pengeluaran/kapita","Skor aset","Luas lantai",
  "Motor","Mobil","Kulkas","AC",
  "Listrik PLN","Lantai layak","Punya BAB",
  "Penerima bansos","Layak bansos"
)
corrplot(mat_cor, method="color", type="upper",
         addCoef.col="black", number.cex=0.65,
         tl.col="black", tl.cex=0.8,
         title="Plot 9: Korelasi Antar Variabel Utama",
         mar=c(0,0,2,0))

# --------------------------------------------------------------------------
# PLOT 10: Exclusion error per kabupaten/kota
# --------------------------------------------------------------------------
p10 <- df |>
  filter(!is.na(layak), !is.na(nama_kab)) |>
  group_by(nama_kab, wilayah) |>
  summarise(exclusion_err=mean(layak==1 & penerima==0, na.rm=TRUE),
            n=n(), .groups="drop") |>
  filter(n >= 50) |>
  ggplot(aes(x=reorder(nama_kab, exclusion_err), y=exclusion_err, fill=wilayah)) +
  geom_col(width=0.7) +
  geom_text(aes(label=percent(exclusion_err, accuracy=1)),
            hjust=-0.1, size=2.8, color="grey30") +
  coord_flip() +
  scale_y_continuous(labels=percent_format(accuracy=1), limits=c(0,0.55),
                     expand=expansion(mult=c(0,0.08))) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 10: Exclusion Error per Kabupaten/Kota Jawa Barat",
    subtitle = "Proporsi RT miskin yang tidak mendapat bansos apapun",
    x=NULL, y="% exclusion error", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p10)

# --------------------------------------------------------------------------
# PLOT 11: Skor aset per kelompok sasaran & wilayah
# --------------------------------------------------------------------------
p11 <- df |>
  filter(!is.na(jenis_sasaran), !is.na(wilayah)) |>
  group_by(wilayah, jenis_sasaran) |>
  summarise(rata_aset=mean(skor_aset, na.rm=TRUE), .groups="drop") |>
  ggplot(aes(x=jenis_sasaran, y=rata_aset, fill=wilayah)) +
  geom_col(position="dodge", width=0.65) +
  geom_text(aes(label=round(rata_aset,2)),
            position=position_dodge(width=0.65),
            vjust=-0.4, size=3.5, fontface="bold", color="grey30") +
  scale_fill_manual(values=warna_wilayah) +
  scale_x_discrete(labels=function(x) str_wrap(x,15)) +
  scale_y_continuous(expand=expansion(mult=c(0,0.15))) +
  labs(
    title    = "Plot 11: Rata-rata Skor Aset per Kelompok Sasaran",
    subtitle = "Skor aset: jumlah dari 5 aset (motor, mobil, kulkas, AC, laptop)",
    x=NULL, y="Rata-rata skor aset (0–5)", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p11)

# --------------------------------------------------------------------------
# PLOT 12: Boxplot kapita per jenis sasaran & wilayah
# --------------------------------------------------------------------------
p12 <- df |>
  filter(!is.na(jenis_sasaran), !is.na(wilayah),
         jenis_sasaran != "Bukan sasaran", kapita_juta < 8) |>
  ggplot(aes(x=jenis_sasaran, y=kapita_juta, fill=wilayah)) +
  geom_boxplot(alpha=0.7, outlier.size=0.5, outlier.alpha=0.2) +
  geom_hline(yintercept=GK_JT, linetype="dashed", color="firebrick", linewidth=0.8) +
  facet_wrap(~wilayah) +
  scale_x_discrete(labels=function(x) str_wrap(x,15)) +
  scale_y_continuous(limits=c(0,5), labels=label_number(suffix=" jt", accuracy=0.1)) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 12: Boxplot Pengeluaran per Kapita berdasarkan Status Sasaran",
    subtitle = "Garis merah putus = garis kemiskinan Jawa Barat",
    x=NULL, y="Pengeluaran per kapita (juta Rp/bln)", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  ) + theme(legend.position="none", axis.text.x=element_text(size=8))
print(p12)

## Warning: Removed 77 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

# --------------------------------------------------------------------------
# PLOT 13: % penerima berdasarkan pendidikan KRT & wilayah
# --------------------------------------------------------------------------
p13 <- df |>
  filter(!is.na(pendidikan_krt), !is.na(wilayah)) |>
  group_by(wilayah, pendidikan_krt) |>
  summarise(pct_penerima=mean(penerima, na.rm=TRUE), n=n(), .groups="drop") |>
  ggplot(aes(x=pendidikan_krt, y=pct_penerima, fill=wilayah)) +
  geom_col(position="dodge", width=0.7) +
  geom_text(aes(label=percent(pct_penerima, accuracy=0.1)),
            position=position_dodge(width=0.7),
            vjust=-0.4, size=3, color="grey30") +
  scale_fill_manual(values=warna_wilayah) +
  scale_y_continuous(labels=percent_format(), expand=expansion(mult=c(0,0.15))) +
  labs(
    title    = "Plot 13: % Penerima Bansos berdasarkan Pendidikan KRT",
    subtitle = "Apakah pendidikan rendah meningkatkan kemungkinan menerima bansos?",
    x="Pendidikan Kepala Rumah Tangga", y="% penerima bansos",
    fill=NULL, caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p13)

# --------------------------------------------------------------------------
# PLOT 14: Exclusion error berdasarkan pendidikan KRT
# --------------------------------------------------------------------------
p14 <- df |>
  filter(!is.na(pendidikan_krt), !is.na(wilayah), !is.na(layak)) |>
  group_by(wilayah, pendidikan_krt) |>
  summarise(pct_exclusion=mean(layak==1 & penerima==0, na.rm=TRUE),
            n=n(), .groups="drop") |>
  ggplot(aes(x=pendidikan_krt, y=pct_exclusion, color=wilayah, group=wilayah)) +
  geom_line(linewidth=1.2) + geom_point(size=3.5) +
  scale_color_manual(values=warna_wilayah) +
  scale_y_continuous(labels=percent_format(accuracy=1),
                     expand=expansion(mult=c(0.02,0.1))) +
  labs(
    title    = "Plot 14: Exclusion Error berdasarkan Pendidikan KRT",
    subtitle = "RT miskin yang tidak menerima bansos — per tingkat pendidikan KRT",
    x="Pendidikan Kepala Rumah Tangga", y="% exclusion error",
    color=NULL, caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p14)

# --------------------------------------------------------------------------
# PLOT 15: % penerima berdasarkan status pekerjaan KRT
# --------------------------------------------------------------------------
p15 <- df |>
  filter(!is.na(status_kerja_krt), !is.na(wilayah)) |>
  group_by(wilayah, status_kerja_krt) |>
  summarise(pct_penerima=mean(penerima, na.rm=TRUE), n=n(), .groups="drop") |>
  filter(n >= 30) |>
  ggplot(aes(x=reorder(status_kerja_krt, pct_penerima), y=pct_penerima, fill=wilayah)) +
  geom_col(position="dodge", width=0.7) +
  geom_text(aes(label=percent(pct_penerima, accuracy=0.1)),
            position=position_dodge(width=0.7),
            hjust=-0.1, size=3, color="grey30") +
  coord_flip() +
  scale_fill_manual(values=warna_wilayah) +
  scale_y_continuous(labels=percent_format(), limits=c(0,0.65),
                     expand=expansion(mult=c(0,0.05))) +
  labs(
    title    = "Plot 15: % Penerima Bansos berdasarkan Status Pekerjaan KRT",
    subtitle = "Apakah status pekerja informal meningkatkan akses bansos?",
    x=NULL, y="% penerima bansos", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p15)

# --------------------------------------------------------------------------
# PLOT 16: % penerima berdasarkan sektor pekerjaan KRT
# --------------------------------------------------------------------------
p16 <- df |>
  filter(!is.na(sektor_krt_label), !is.na(wilayah)) |>
  group_by(wilayah, sektor_krt_label) |>
  summarise(pct_penerima=mean(penerima, na.rm=TRUE), n=n(), .groups="drop") |>
  filter(n >= 30) |>
  ggplot(aes(x=sektor_krt_label, y=pct_penerima, fill=wilayah)) +
  geom_col(position="dodge", width=0.7) +
  geom_text(aes(label=percent(pct_penerima, accuracy=0.1)),
            position=position_dodge(width=0.7),
            vjust=-0.4, size=3, color="grey30") +
  scale_fill_manual(values=warna_wilayah) +
  scale_y_continuous(labels=percent_format(), expand=expansion(mult=c(0,0.18))) +
  labs(
    title    = "Plot 16: % Penerima Bansos berdasarkan Sektor Pekerjaan KRT",
    subtitle = "Pertanian vs Industri vs Perdagangan/Jasa",
    x=NULL, y="% penerima bansos", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p16)

# --------------------------------------------------------------------------
# PLOT 17: Distribusi umur KRT per jenis sasaran (boxplot)
# --------------------------------------------------------------------------
p17 <- df |>
  filter(!is.na(jenis_sasaran), !is.na(wilayah), !is.na(umur_krt)) |>
  ggplot(aes(x=jenis_sasaran, y=umur_krt, fill=wilayah)) +
  geom_boxplot(alpha=0.7, outlier.size=0.5, outlier.alpha=0.2) +
  facet_wrap(~wilayah) +
  scale_x_discrete(labels=function(x) str_wrap(x,15)) +
  scale_fill_manual(values=warna_wilayah) +
  labs(
    title    = "Plot 17: Distribusi Umur KRT berdasarkan Status Sasaran",
    subtitle = "Apakah KRT yang lebih tua lebih rentan terhadap exclusion error?",
    x=NULL, y="Umur Kepala Rumah Tangga (tahun)", fill=NULL,
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  ) + theme(legend.position="none", axis.text.x=element_text(size=8))
print(p17)

# --------------------------------------------------------------------------
# PLOT 18: Heatmap % penerima per pendidikan KRT x kelompok desil
# --------------------------------------------------------------------------
heatmap_data <- df |>
  filter(!is.na(pendidikan_krt), !is.na(kelompok_ekono), !is.na(wilayah)) |>
  group_by(wilayah, pendidikan_krt, kelompok_ekono) |>
  summarise(pct_penerima=mean(penerima, na.rm=TRUE), n=n(), .groups="drop") |>
  filter(n >= 20)

p18 <- heatmap_data |>
  ggplot(aes(x=pendidikan_krt, y=kelompok_ekono, fill=pct_penerima)) +
  geom_tile(color="white", linewidth=0.5) +
  geom_text(aes(label=percent(pct_penerima, accuracy=1)),
            size=2.8, color="white", fontface="bold") +
  facet_wrap(~wilayah) +
  scale_fill_gradient(low="#FEF3C7", high="#DC2626", labels=percent_format()) +
  scale_x_discrete(labels=function(x) str_wrap(x,10)) +
  labs(
    title    = "Plot 18: Heatmap % Penerima Bansos",
    subtitle = "Pendidikan KRT × Kelompok Desil Pengeluaran",
    x="Pendidikan KRT", y=NULL, fill="% penerima",
    caption="Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  ) + theme(axis.text.x=element_text(size=7.5))
print(p18)

Plot 1

# =============================================================================
# 6. PLOT TEBARAN + MEDIAN LINE (digabung lintas wilayah)
# =============================================================================
# Tujuan: melihat sebaran individu RT, bukan hanya ringkasan agregat.
# Tiga sudut pandang sesuai konteks analisis:
#   Scatter A — kapita × jenis sasaran  (inti analisis salah sasaran)
#   Scatter B — kapita × desil          (kebocoran ke desil atas)
#   Scatter C — kapita × pendidikan KRT (profil demografis)
# =============================================================================

# --- Helper: hitung median + IQR per grup ------------------------------------
median_iqr <- function(data, grup_var, y_var = "kapita_juta") {
  data |>
    group_by(across(all_of(grup_var))) |>
    summarise(
      med = median(.data[[y_var]], na.rm = TRUE),
      q1  = quantile(.data[[y_var]], 0.25, na.rm = TRUE),
      q3  = quantile(.data[[y_var]], 0.75, na.rm = TRUE),
      n   = n(),
      .groups = "drop"
    )
}

# Subsample 40% agar titik tidak terlalu padat (masih ~10.000 titik)
set.seed(42)
df_samp <- df |>
  filter(!is.na(kapita_juta), kapita_juta < 8) |>
  slice_sample(prop = 0.2)

# --------------------------------------------------------------------------
# SCATTER A: Pengeluaran per kapita × Jenis Sasaran Bansos
# Pertanyaan: seberapa bersih pemisahan miskin vs tidak miskin per kategori?
# --------------------------------------------------------------------------
med_A <- median_iqr(
  df |> filter(!is.na(jenis_sasaran), kapita_juta < 8),
  "jenis_sasaran"
)

psA <- df_samp |>
  filter(!is.na(jenis_sasaran)) |>
  ggplot(aes(x = jenis_sasaran, y = kapita_juta, color = jenis_sasaran)) +
  
  # Titik tebaran
   geom_jitter(width = 0.28, alpha = 0.15, size = 0.65, show.legend = FALSE) +

  # Bar IQR (Q1–Q3)
    geom_errorbar(data = med_A,
                aes(x = jenis_sasaran, y = med, ymin = q1, ymax = q3),
                width = 0.25, linewidth = 1.2,
                color = "grey25", inherit.aes = FALSE) +

  # Garis median (crossbar tebal)
   geom_crossbar(data = med_A,
                aes(x = jenis_sasaran, y = med, ymin = med, ymax = med),
                width = 0.55, linewidth = 1.5, fatten = 0,
                color = "black", inherit.aes = FALSE) +

  # Label nilai median
 geom_text(data = med_A,
            aes(x = jenis_sasaran, y = med,
                label = paste0("Md: ", label_number(suffix=" jt", accuracy=0.2)(med))),
            hjust = -0.8, vjust = 0.4, size = 3,
            color = "black", fontface = "bold", inherit.aes = FALSE) +

  # Garis kemiskinan
  geom_hline(yintercept = GK_JT,
             linetype = "dashed", color = "firebrick", linewidth = 0.9) +
  annotate("text", x = 0.55, y = GK_JT + 0.3,
           label = "GK Rp 480 rb", color = "firebrick",
           size = 2.8, hjust = 0, fontface = "italic") +

  scale_color_manual(values = warna_sasaran) +
  scale_x_discrete(labels = function(x) str_wrap(x, 14)) +
  scale_y_continuous(limits = c(0, 6),
                     labels = label_number(suffix = " jt", accuracy = 0.1)) +
  labs(
    title    = "Scatter A: Pengeluaran Per Kapita × Status Sasaran Bansos",
    subtitle = "Titik = RT (subsample 40%) | Garis hitam = median | Bar abu = IQR (Q1–Q3) | Semua wilayah",
    x        = NULL,
    y        = "Pengeluaran per kapita (juta Rp/bulan)",
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(psA)

## Warning: The `fatten` argument of `geom_crossbar()` is deprecated as of ggplot2 4.0.0.
## ℹ Please use the `middle.linewidth` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 84 rows containing missing values or values outside the scale range
## (`geom_point()`).

Plot 2

# --------------------------------------------------------------------------
# SCATTER B: Pengeluaran per kapita × Desil
# Pertanyaan: apakah bansos benar-benar terkonsentrasi di desil bawah?
# --------------------------------------------------------------------------
med_B <- median_iqr(
  df |> filter(!is.na(desil), kapita_juta < 8),
  "desil"
)

df_samp_B <- df_samp |>
  filter(!is.na(desil)) |>
  mutate(status_penerima = factor(penerima, levels=c(1,0),
                                  labels=c("Penerima bansos","Bukan penerima")))

psB <- df_samp_B |>
  ggplot(aes(x = factor(desil), y = kapita_juta)) +

  # Titik tebaran, warna = status penerima
  geom_jitter(aes(color = status_penerima),
              width = 0.3, alpha = 0.7, size = 0.65) +

  # Bar IQR
  geom_errorbar(data = med_B,
                aes(x = factor(desil), y = med, ymin = q1, ymax = q3),
                width = 0.3, linewidth = 1.0,
                color = "grey20", inherit.aes = FALSE) +

  # Garis median sambung antar desil
  geom_line(data = med_B,
            aes(x = factor(desil), y = med, group = 1),
            linewidth = 1.4, color = "black", inherit.aes = FALSE) +
  geom_point(data = med_B,
             aes(x = factor(desil), y = med),
             size = 3.8, color = "black", shape = 18, inherit.aes = FALSE) +

  # Garis kemiskinan
  geom_hline(yintercept = GK_JT,
             linetype = "dashed", color = "firebrick", linewidth = 0.9) +
  annotate("text", x = "9", y = GK_JT + 0.3,
           label = "GK Rp 480 rb", color = "firebrick",
           size = 2.8, hjust = 0.1, fontface = "italic") +

  scale_color_manual(values = c("Penerima bansos"="#1d12e3","Bukan penerima"="#f3ed2e")) +
  scale_x_discrete(labels = paste0("D", 1:10)) +
  scale_y_continuous(limits = c(0, 7),
                     labels = label_number(suffix = " jt", accuracy = 0.1)) +
  guides(color = guide_legend(override.aes = list(size = 4, alpha = 1)))+
  labs(
    title    = "Scatter B: Pengeluaran Per Kapita × Desil",
    subtitle = "Warna = status penerima | ◆ = median desil | Bar abu = IQR | Semua wilayah",
    x        = "Desil pengeluaran per kapita (D1 = termiskin)",
    y        = "Pengeluaran per kapita (juta Rp/bulan)",
    color    = NULL,
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(psB)

## Warning: Removed 38 rows containing missing values or values outside the scale range
## (`geom_point()`).

Plot 3

# --------------------------------------------------------------------------
# SCATTER C: Pengeluaran per kapita × Pendidikan KRT
# Pertanyaan: bagaimana profil ekonomi RT berdasarkan pendidikan KRT,
#             dan apakah yang berpendidikan rendah lebih banyak yang exclusion error?
# --------------------------------------------------------------------------
med_C <- median_iqr(
  df |> filter(!is.na(pendidikan_krt), kapita_juta < 8),
  "pendidikan_krt"
)

psC <- df_samp |>
  filter(!is.na(pendidikan_krt), !is.na(jenis_sasaran)) |>
  ggplot(aes(x = pendidikan_krt, y = kapita_juta, color = jenis_sasaran)) +

  # Titik tebaran, warna = jenis sasaran
  geom_jitter(width = 0.28, alpha = 0.5, size = 0.65) +

  # Bar IQR
  geom_errorbar(data = med_C,
                aes(x = pendidikan_krt, y = med, ymin = q1, ymax = q3),
                width = 0.25, linewidth = 1.1,
                color = "grey25", inherit.aes = FALSE) +

  # Garis median sambung
  geom_line(data = med_C,
            aes(x = pendidikan_krt, y = med, group = 1),
            linewidth = 1.4, color = "black", inherit.aes = FALSE) +
  geom_point(data = med_C,
             aes(x = pendidikan_krt, y = med),
             size = 4, color = "black", shape = 18, inherit.aes = FALSE) +

  # Label nilai median
  geom_text(data = med_C,
            aes(x = pendidikan_krt, y = med+0.12,
                label = label_number(suffix=" jt", accuracy=0.1)(med)),
            nudge_x=0.12, nudge_y = 0.22, size = 3.0, color = "black",
            fontface = "bold", inherit.aes = FALSE) +

  # Garis kemiskinan
  geom_hline(yintercept = GK_JT,
             linetype = "dashed", color = "firebrick", linewidth = 0.9) +
  annotate("text", x = 0.6, y = GK_JT - 0.2,
           label = "GK Rp 480 rb", color = "firebrick",
           size = 2.8, hjust = 0.38, fontface = "italic") +

  scale_color_manual(values = warna_sasaran) +
  scale_x_discrete(labels = function(x) str_wrap(x, 12)) +
  scale_y_continuous(limits = c(0, 6),
                     labels = label_number(suffix = " jt", accuracy = 0.1)) +
  guides(color = guide_legend(override.aes = list(alpha=0.9, size=2.5), nrow=2)) +
  labs(
    subtitle = "Warna = status sasaran bansos | ◆ = median | Bar abu = IQR | Semua wilayah",
    x        = "Pendidikan Kepala Rumah Tangga",
    y        = "Pengeluaran per kapita (juta Rp/bulan)",
    color    = NULL,
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(psC)

## Warning: Removed 84 rows containing missing values or values outside the scale range
## (`geom_point()`).

# =============================================================================
# 7. UJI STATISTIK
# =============================================================================

cat("\n\n========== UJI STATISTIK ==========\n")

## 
## 
## ========== UJI STATISTIK ==========

cat("\n[1] Chi-square: Penerimaan Bansos vs Wilayah\n")

## 
## [1] Chi-square: Penerimaan Bansos vs Wilayah

print(chisq.test(table(df$wilayah, df$penerima)))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df$wilayah, df$penerima)
## X-squared = 681.26, df = 1, p-value < 2.2e-16

cat("\n[2] Chi-square: Jenis Sasaran vs Wilayah\n")

## 
## [2] Chi-square: Jenis Sasaran vs Wilayah

print(chisq.test(table(df$wilayah, df$jenis_sasaran)))

## 
##  Pearson's Chi-squared test
## 
## data:  table(df$wilayah, df$jenis_sasaran)
## X-squared = 689.77, df = 3, p-value < 2.2e-16

cat("\n[3] Chi-square: Jenis Sasaran vs Pendidikan KRT\n")

## 
## [3] Chi-square: Jenis Sasaran vs Pendidikan KRT

print(chisq.test(table(df$pendidikan_krt, df$jenis_sasaran)))

## 
##  Pearson's Chi-squared test
## 
## data:  table(df$pendidikan_krt, df$jenis_sasaran)
## X-squared = 1630.4, df = 12, p-value < 2.2e-16

cat("\n[4] Chi-square: Jenis Sasaran vs Status Pekerjaan KRT\n")

## 
## [4] Chi-square: Jenis Sasaran vs Status Pekerjaan KRT

print(chisq.test(table(df$status_kerja_krt, df$jenis_sasaran)))

## 
##  Pearson's Chi-squared test
## 
## data:  table(df$status_kerja_krt, df$jenis_sasaran)
## X-squared = 733.74, df = 15, p-value < 2.2e-16

cat("\n[5] Mann-Whitney: Kapita ~ Penerima (per Wilayah)\n")

## 
## [5] Mann-Whitney: Kapita ~ Penerima (per Wilayah)

df |>
  filter(!is.na(wilayah)) |>
  group_by(wilayah) |>
  group_split() |>
  walk(function(d) {
    w <- as.character(unique(d$wilayah))
    cat("\nWilayah:", w, "\n")
    mw  <- wilcox.test(kapita ~ penerima, data = d)
    med <- d |> group_by(penerima) |>
      summarise(med = median(kapita, na.rm=TRUE))
    cat("  W =", mw$statistic, "| p-value =", round(mw$p.value, 6), "\n")
    cat("  Median kapita penerima      : Rp", comma(med$med[med$penerima==1]), "\n")
    cat("  Median kapita bukan penerima: Rp", comma(med$med[med$penerima==0]), "\n")
  })

## 
## Wilayah: Perkotaan 
##   W = 42760984 | p-value = 0 
##   Median kapita penerima      : Rp 965,179 
##   Median kapita bukan penerima: Rp 1,552,504 
## 
## Wilayah: Perdesaan 
##   W = 11773980 | p-value = 0 
##   Median kapita penerima      : Rp 988,947 
##   Median kapita bukan penerima: Rp 1,184,944

cat("\n[6] Kruskal-Wallis: Kapita ~ Pendidikan KRT\n")

## 
## [6] Kruskal-Wallis: Kapita ~ Pendidikan KRT

print(kruskal.test(kapita ~ pendidikan_krt, data = df))

## 
##  Kruskal-Wallis rank sum test
## 
## data:  kapita by pendidikan_krt
## Kruskal-Wallis chi-squared = 1903.5, df = 4, p-value < 2.2e-16

cat("\n[7] Tabel Silang Wilayah x Jenis Sasaran\n")

## 
## [7] Tabel Silang Wilayah x Jenis Sasaran

df |>
  filter(!is.na(wilayah), !is.na(jenis_sasaran)) |>
  tabyl(wilayah, jenis_sasaran) |>
  adorn_percentages("row") |>
  adorn_pct_formatting(digits=1) |>
  adorn_ns() |>
  print()

##    wilayah Tepat sasaran Exclusion error (miskin, tidak dapat)
##  Perkotaan    2.1% (352)                            2.0% (337)
##  Perdesaan    3.0% (269)                            2.0% (175)
##  Inclusion error (mampu, tapi dapat)  Bukan sasaran
##                        29.6% (5,013) 66.4% (11,251)
##                        45.1% (4,035) 49.9%  (4,458)

cat("\n[8] Tabel Silang Pendidikan KRT x Jenis Sasaran\n")

## 
## [8] Tabel Silang Pendidikan KRT x Jenis Sasaran

df |>
  filter(!is.na(pendidikan_krt), !is.na(jenis_sasaran)) |>
  tabyl(pendidikan_krt, jenis_sasaran) |>
  adorn_percentages("row") |>
  adorn_pct_formatting(digits=1) |>
  adorn_ns() |>
  print()

##        pendidikan_krt Tepat sasaran Exclusion error (miskin, tidak dapat)
##  Tidak/Belum tamat SD    3.0%  (18)                            3.0%  (18)
##                    SD    3.7% (358)                            2.6% (251)
##                   SMP    2.2%  (98)                            2.2%  (98)
##               SMA/SMK    0.9%  (57)                            1.2%  (80)
##      Diploma/Sarjana+    1.9%  (90)                            1.4%  (65)
##  Inclusion error (mampu, tapi dapat) Bukan sasaran
##                        50.4%   (299) 43.5%   (258)
##                        46.7% (4,514) 47.0% (4,545)
##                        34.0% (1,491) 61.5% (2,700)
##                        21.9% (1,434) 76.0% (4,983)
##                        27.9% (1,310) 68.8% (3,223)

# =============================================================================
# 8. PROFIL KELAS MENENGAH TERJEPIT
# =============================================================================

cat("\n\n========== PROFIL KELAS MENENGAH TERJEPIT ==========\n")

## 
## 
## ========== PROFIL KELAS MENENGAH TERJEPIT ==========

profil_kejepit <- df |>
  filter(kejepit == 1) |>
  group_by(wilayah) |>
  summarise(
    n_rt              = n(),
    median_kapita     = median(kapita,   na.rm=TRUE),
    median_desil      = median(desil,    na.rm=TRUE),
    median_umur_krt   = median(umur_krt, na.rm=TRUE),
    pct_krt_perempuan = mean(jk_krt == 2, na.rm=TRUE),
    pct_motor         = mean(punya_motor,  na.rm=TRUE),
    pct_kulkas        = mean(punya_kulkas, na.rm=TRUE),
    pct_mobil         = mean(punya_mobil,  na.rm=TRUE),
    pct_listrik_pln   = mean(listrik_pln,  na.rm=TRUE),
    .groups = "drop"
  )

profil_kejepit |>
  gt() |>
  tab_header(
    title    = "Profil Kelas Menengah Terjepit — Jawa Barat 2023",
    subtitle = "RT di atas GK, desil < 9, tidak menerima bansos apapun"
  ) |>
  fmt_number(columns  = c(n_rt, median_kapita, median_umur_krt, median_desil),
             decimals=0, use_seps=TRUE) |>
  fmt_percent(columns = starts_with("pct_"), decimals=1)

Profil Kelas Menengah Terjepit — Jawa Barat 2023
RT di atas GK, desil < 9, tidak menerima bansos apapun
wilayah	n_rt	median_kapita	median_desil	median_umur_krt	pct_krt_perempuan	pct_motor	pct_kulkas	pct_mobil	pct_listrik_pln
Perkotaan	7,552	1,155,765	5	47	12.3%	79.5%	76.8%	8.2%	95.8%
Perdesaan	3,768	1,091,487	5	48	12.3%	73.4%	59.0%	5.5%	93.5%

# =============================================================================
# 9. DASHBOARD — GABUNG PLOT UTAMA
# =============================================================================

dashboard_utama <- (p1 | p2) /
  (p3 | p4) /
  (p13 | p15) +
  plot_annotation(
    title    = "EDA Salah Sasaran Bansos — Jawa Barat, SUSENAS Maret 2023",
    subtitle = "25.890 rumah tangga · Perkotaan vs Perdesaan · Termasuk demografi KRT",
    caption  = "Sumber: BPS Jawa Barat | GK Rp 479.934/kapita/bulan",
    theme    = theme(
      plot.title    = element_text(size=16, face="bold"),
      plot.subtitle = element_text(size=12, color="grey40"),
      plot.caption  = element_text(size=9,  color="grey55")
    )
  )

dashboard_scatter <- psA / psB / psC +
  plot_annotation(
    title    = "Plot Tebaran + Median Line — Analisis Bansos Jawa Barat 2023",
    subtitle = "Semua wilayah digabung | ◆ = median | Bar abu = IQR | Garis merah putus = GK",
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat",
    theme    = theme(
      plot.title    = element_text(size=15, face="bold"),
      plot.subtitle = element_text(size=11, color="grey40")
    )
  )

ggsave("dashboard_bansos_jabar2023_utama.png",
       plot=dashboard_utama, width=18, height=22, dpi=150, bg="white")

## Warning: Removed 515 rows containing non-finite outside the scale range
## (`stat_density()`).

ggsave("dashboard_bansos_jabar2023_scatter.png",
       plot=dashboard_scatter, width=14, height=22, dpi=150, bg="white")

## Warning: Removed 84 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 38 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 84 rows containing missing values or values outside the scale range
## (`geom_point()`).

cat("\nDashboard tersimpan:\n")

## 
## Dashboard tersimpan:

cat("  dashboard_bansos_jabar2023_utama.png\n")

##   dashboard_bansos_jabar2023_utama.png

cat("  dashboard_bansos_jabar2023_scatter.png\n")

##   dashboard_bansos_jabar2023_scatter.png

# =============================================================================
# 10. MODEL PREDIKSI EXCLUSION ERROR
# =============================================================================
# Outcome  : RT miskin (layak==1) yang TIDAK menerima bansos → exclusion error
# Metode   : (A) Logistic Regression + (B) Random Forest
# Evaluasi : Confusion Matrix, ROC-AUC, Forest Plot OR, Variable Importance
# =============================================================================

# --- Install & load package tambahan -----------------------------------------
# install.packages(c("broom", "pROC", "ranger", "caret", "gt"))

library(broom)   # tidy() untuk output model → tabel rapi

## Warning: package 'broom' was built under R version 4.5.2

library(pROC)    # ROC curve & AUC

## Warning: package 'pROC' was built under R version 4.5.3

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(ranger)  # Random Forest versi cepat

## Warning: package 'ranger' was built under R version 4.5.3

library(caret)   # confusionMatrix(), createDataPartition()

## Warning: package 'caret' was built under R version 4.5.3

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

set.seed(42)

# =============================================================================
# 10.1 PERSIAPAN DATA MODEL
# =============================================================================
# Filter: HANYA RT miskin (layak == 1)
# Logika: kita ingin tahu FAKTOR APA yang membuat RT miskin tidak dapat bansos
# Jika pakai semua RT, model akan didominasi oleh kapita (bukan itu tujuannya)

df_model <- df |>
  filter(
    layak == 1,                    # hanya RT di bawah garis kemiskinan
    !is.na(pendidikan_krt),
    !is.na(status_kerja_krt),
    !is.na(sektor_krt_label),
    !is.na(umur_krt),
    !is.na(wilayah),
    !is.na(skor_aset)
  ) |>
  mutate(
    # Outcome: 1 = exclusion error (miskin tapi tidak dapat), 0 = tepat sasaran
    exclusion = factor(
      if_else(penerima == 0, "Ya", "Tidak"),
      levels = c("Tidak", "Ya")   # "Tidak" = referensi (tidak exclusion error)
    ),
    # Jumlah ART sebagai prediktor
    jml_art = as.numeric(r301)
  )

# Cek distribusi outcome
cat("=== Distribusi Outcome ===\n")

## === Distribusi Outcome ===

cat("Tepat sasaran (Tidak)  :", sum(df_model$exclusion == "Tidak"),
    sprintf("(%.1f%%)\n", 100*mean(df_model$exclusion == "Tidak")))

## Tepat sasaran (Tidak)  : 621 (54.8%)

cat("Exclusion error (Ya)   :", sum(df_model$exclusion == "Ya"),
    sprintf("(%.1f%%)\n", 100*mean(df_model$exclusion == "Ya")))

## Exclusion error (Ya)   : 512 (45.2%)

# --- Cek class imbalance -----------------------------------------------------
# Jika rasio 1:3 atau lebih parah, pertimbangkan SMOTE / class weight
# Untuk saat ini kita pakai model standar dulu, lihat hasilnya

# --- Formula prediktor -------------------------------------------------------
formula_model <- exclusion ~
  wilayah + pendidikan_krt + umur_krt + jk_krt_label +
  status_kerja_krt + sektor_krt_label + skor_aset +
  listrik_pln + rumah_milik + lantai_layak + punya_bab + air_layak +
  jml_art

# --- Train/Test Split 70:30 --------------------------------------------------
idx_train <- createDataPartition(df_model$exclusion, p = 0.70, list = FALSE)
df_train  <- df_model[ idx_train, ]
df_test   <- df_model[-idx_train, ]

cat("\nTrain:", nrow(df_train), "| Test:", nrow(df_test), "\n")

## 
## Train: 794 | Test: 339

# =============================================================================
# (A) LOGISTIC REGRESSION
# =============================================================================

logit_fit <- glm(
  formula_model,
  data   = df_train,
  family = binomial(link = "logit")   # link = "probit" jika mau probit
)

cat("\n========== RINGKASAN LOGISTIC REGRESSION ==========\n")

## 
## ========== RINGKASAN LOGISTIC REGRESSION ==========

print(summary(logit_fit))

## 
## Call:
## glm(formula = formula_model, family = binomial(link = "logit"), 
##     data = df_train)
## 
## Coefficients:
##                                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                             1.797536   1.183556   1.519  0.12882   
## wilayahPerdesaan                       -0.302210   0.161392  -1.873  0.06113 . 
## pendidikan_krtSD                       -0.220717   0.432585  -0.510  0.60989   
## pendidikan_krtSMP                      -0.037831   0.461869  -0.082  0.93472   
## pendidikan_krtSMA/SMK                   0.125425   0.476752   0.263  0.79249   
## pendidikan_krtDiploma/Sarjana+          0.028227   0.465089   0.061  0.95160   
## umur_krt                               -0.014938   0.006922  -2.158  0.03093 * 
## jk_krt_labelPerempuan                   0.153864   0.265297   0.580  0.56193   
## status_kerja_krtBerusaha sendiri       -0.177374   0.487644  -0.364  0.71605   
## status_kerja_krtBerusaha (punya buruh)  0.235121   0.512008   0.459  0.64608   
## status_kerja_krtBuruh/Karyawan          0.057804   0.501473   0.115  0.90823   
## status_kerja_krtPekerja bebas          -0.044167   0.497584  -0.089  0.92927   
## status_kerja_krtPekerja keluarga       -0.629525   1.273461  -0.494  0.62106   
## sektor_krt_labelPertanian               0.296094   0.580627   0.510  0.61008   
## sektor_krt_labelPerdagangan/Jasa       -0.229111   0.456609  -0.502  0.61583   
## sektor_krt_labelLainnya                -0.179721   0.526034  -0.342  0.73261   
## skor_aset                               0.300024   0.102814   2.918  0.00352 **
## listrik_pln                            -0.145432   0.244712  -0.594  0.55231   
## rumah_milik                            -0.116547   0.203137  -0.574  0.56615   
## lantai_layak                           -0.707498   0.896788  -0.789  0.43016   
## punya_bab                              -0.003708   0.181421  -0.020  0.98369   
## air_layak                               0.247482   0.169484   1.460  0.14423   
## jml_art                                -0.158277   0.058311  -2.714  0.00664 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1093.4  on 793  degrees of freedom
## Residual deviance: 1046.0  on 771  degrees of freedom
## AIC: 1092
## 
## Number of Fisher Scoring iterations: 4

# --- McFadden Pseudo R² (ukuran goodness-of-fit untuk logit) -----------------
logit_null <- glm(exclusion ~ 1, data = df_train, family = binomial)
mcfadden_r2 <- 1 - (logit_fit$deviance / logit_null$deviance)
cat("McFadden Pseudo R²:", round(mcfadden_r2, 4), "\n")

## McFadden Pseudo R²: 0.0433

# Cara lebih simpel:
mcfadden_r2 <- 1 - (logit_fit$deviance / logit_null$deviance)
cat("\nMcFadden Pseudo R²:", round(mcfadden_r2, 4), "\n")

## 
## McFadden Pseudo R²: 0.0433

# Interpretasi: > 0.2 = acceptable, > 0.4 = excellent (McFadden 1977)

# --- Odds Ratio + CI 95% -----------------------------------------------------
tabel_or <- logit_fit |>
  tidy(conf.int = TRUE, exponentiate = TRUE) |>   # exponentiate = TRUE → exp(coef) = OR
  filter(term != "(Intercept)") |>
  mutate(
    signif = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ ".",
      TRUE            ~ ""
    ),
    arah = if_else(estimate > 1, "Meningkatkan risiko", "Menurunkan risiko")
  ) |>
  arrange(desc(estimate))

cat("\n--- Odds Ratio (OR > 1 = meningkatkan risiko exclusion error) ---\n")

## 
## --- Odds Ratio (OR > 1 = meningkatkan risiko exclusion error) ---

print(tabel_or |> select(term, estimate, conf.low, conf.high, p.value, signif))

## # A tibble: 22 × 6
##    term                               estimate conf.low conf.high p.value signif
##    <chr>                                 <dbl>    <dbl>     <dbl>   <dbl> <chr> 
##  1 skor_aset                             1.35     1.10      1.65  0.00352 "**"  
##  2 sektor_krt_labelPertanian             1.34     0.435     4.29  0.610   ""    
##  3 air_layak                             1.28     0.919     1.79  0.144   ""    
##  4 status_kerja_krtBerusaha (punya b…    1.27     0.454     3.42  0.646   ""    
##  5 jk_krt_labelPerempuan                 1.17     0.693     1.97  0.562   ""    
##  6 pendidikan_krtSMA/SMK                 1.13     0.443     2.91  0.792   ""    
##  7 status_kerja_krtBuruh/Karyawan        1.06     0.387     2.80  0.908   ""    
##  8 pendidikan_krtDiploma/Sarjana+        1.03     0.412     2.59  0.952   ""    
##  9 punya_bab                             0.996    0.699     1.42  0.984   ""    
## 10 umur_krt                              0.985    0.972     0.999 0.0309  "*"   
## # ℹ 12 more rows

# --- Forest Plot Odds Ratio --------------------------------------------------
p_or <- tabel_or |>
  ggplot(aes(x = estimate, y = reorder(term, estimate), color = arah)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey60", linewidth = 0.8) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0.3, linewidth = 0.7, color = "grey40"
  ) +
  geom_point(size = 3.5) +
  scale_color_manual(values = c(
    "Meningkatkan risiko" = "#ef4c5b",
    "Menurunkan risiko"   = "#2563EB"
  )) +
  scale_x_log10(
    labels = label_number(accuracy = 0.01)
  ) +
  labs(
    title    = "Forest Plot: Odds Ratio Logistic Regression",
    subtitle = "Prediksi Exclusion Error (RT miskin yang tidak dapat bansos) | CI 95%",
    x        = "Odds Ratio (skala log) | Garis putus = OR 1 (tidak berpengaruh)",
    y        = NULL,
    color    = NULL,
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )

## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(p_or)

## `height` was translated to `width`.

# --- Prediksi & Evaluasi Logit -----------------------------------------------
pred_logit_prob <- predict(logit_fit, newdata = df_test, type = "response")
pred_logit_klas <- factor(
  if_else(pred_logit_prob >= 0.5, "Ya", "Tidak"),
  levels = c("Tidak", "Ya")
)

cat("\n--- Confusion Matrix: Logistic Regression ---\n")

## 
## --- Confusion Matrix: Logistic Regression ---

cm_logit <- confusionMatrix(pred_logit_klas, df_test$exclusion, positive = "Ya")
print(cm_logit)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Tidak  Ya
##      Tidak   133  80
##      Ya       53  73
##                                         
##                Accuracy : 0.6077        
##                  95% CI : (0.5535, 0.66)
##     No Information Rate : 0.5487        
##     P-Value [Acc > NIR] : 0.01632       
##                                         
##                   Kappa : 0.1952        
##                                         
##  Mcnemar's Test P-Value : 0.02417       
##                                         
##             Sensitivity : 0.4771        
##             Specificity : 0.7151        
##          Pos Pred Value : 0.5794        
##          Neg Pred Value : 0.6244        
##              Prevalence : 0.4513        
##          Detection Rate : 0.2153        
##    Detection Prevalence : 0.3717        
##       Balanced Accuracy : 0.5961        
##                                         
##        'Positive' Class : Ya            
## 

# Perhatikan: Sensitivity (recall) penting di sini!
# Sensitivity tinggi = model bagus menangkap exclusion error yang sesungguhnya

roc_logit <- roc(df_test$exclusion, pred_logit_prob,
                 levels = c("Tidak", "Ya"), direction = "<", quiet = TRUE)
cat("\nAUC Logistic Regression:", round(auc(roc_logit), 4), "\n")

## 
## AUC Logistic Regression: 0.6572

# =============================================================================
# (B) RANDOM FOREST
# =============================================================================
# Package: ranger (jauh lebih cepat dari randomForest)
# probability = TRUE agar bisa hitung ROC-AUC

cat("\n========== RANDOM FOREST ==========\n")

## 
## ========== RANDOM FOREST ==========

rf_fit <- ranger(
  formula_model,
  data          = df_train,
  num.trees     = 500,           # 500 tree sudah cukup stabil
  mtry          = 4,             # ≈ sqrt(jumlah prediktor), bisa di-tune
  min.node.size = 10,            # minimum observasi di leaf node
  importance    = "permutation", # untuk variable importance plot
  probability   = TRUE,          # output probabilitas, bukan kelas
  seed          = 42
)

cat("OOB Prediction Error:", round(rf_fit$prediction.error, 4), "\n")

## OOB Prediction Error: 0.2475

# OOB = Out-of-Bag error, estimasi error tanpa test set terpisah

# --- Prediksi & Evaluasi RF --------------------------------------------------
pred_rf_prob <- predict(rf_fit, data = df_test)$predictions[, "Ya"]
pred_rf_klas <- factor(
  if_else(pred_rf_prob >= 0.5, "Ya", "Tidak"),
  levels = c("Tidak", "Ya")
)

cat("\n--- Confusion Matrix: Random Forest ---\n")

## 
## --- Confusion Matrix: Random Forest ---

cm_rf <- confusionMatrix(pred_rf_klas, df_test$exclusion, positive = "Ya")
print(cm_rf)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Tidak  Ya
##      Tidak   132  80
##      Ya       54  73
##                                           
##                Accuracy : 0.6047          
##                  95% CI : (0.5505, 0.6571)
##     No Information Rate : 0.5487          
##     P-Value [Acc > NIR] : 0.02137         
##                                           
##                   Kappa : 0.1897          
##                                           
##  Mcnemar's Test P-Value : 0.03080         
##                                           
##             Sensitivity : 0.4771          
##             Specificity : 0.7097          
##          Pos Pred Value : 0.5748          
##          Neg Pred Value : 0.6226          
##              Prevalence : 0.4513          
##          Detection Rate : 0.2153          
##    Detection Prevalence : 0.3746          
##       Balanced Accuracy : 0.5934          
##                                           
##        'Positive' Class : Ya              
## 

roc_rf <- roc(df_test$exclusion, pred_rf_prob,
              levels = c("Tidak", "Ya"), direction = "<", quiet = TRUE)
cat("\nAUC Random Forest:", round(auc(roc_rf), 4), "\n")

## 
## AUC Random Forest: 0.6212

# --- Variable Importance Plot ------------------------------------------------
vip_df <- tibble(
  variable   = names(rf_fit$variable.importance),
  importance = rf_fit$variable.importance
) |>
  arrange(desc(importance)) |>
  slice_head(n = 15)   # top 15

p_vip <- vip_df |>
  ggplot(aes(x = importance, y = reorder(variable, importance),
             fill = importance)) +
  geom_col(show.legend = FALSE, alpha = 0.85) +
  geom_text(aes(label = round(importance, 4)),
            hjust = -0.1, size = 3, color = "grey30") +
  scale_fill_gradient(low = "#93C5FD", high = "#1D4ED8") +
  scale_x_continuous(expand = expansion(mult = c(0, 0.25))) +
  labs(
    title    = "Variable Importance — Random Forest",
    subtitle = "Permutation importance: seberapa turun akurasi jika variabel di-shuffle",
    x        = "Importance (permutation)",
    y        = NULL,
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p_vip)

# =============================================================================
# (C) PERBANDINGAN MODEL — ROC Curve & Tabel Performa
# =============================================================================

# --- Plot ROC berdampingan ---------------------------------------------------
roc_df <- bind_rows(
  tibble(
    Model = sprintf("Logistic Regression (AUC = %.3f)", auc(roc_logit)),
    FPR   = 1 - roc_logit$specificities,
    TPR   = roc_logit$sensitivities
  ),
  tibble(
    Model = sprintf("Random Forest (AUC = %.3f)", auc(roc_rf)),
    FPR   = 1 - roc_rf$specificities,
    TPR   = roc_rf$sensitivities
  )
)

p_roc <- roc_df |>
  ggplot(aes(x = FPR, y = TPR, color = Model)) +
  geom_line(linewidth = 1.3) +
  geom_abline(intercept = 0, slope = 1,
              linetype = "dashed", color = "grey60", linewidth = 0.7) +
  annotate("text", x = 0.75, y = 0.2,
           label = "Random guess\n(AUC = 0.5)",
           color = "grey50", size = 3, fontface = "italic") +
  scale_color_manual(values = c("#ef4c5b", "#2563EB")) +
  coord_equal() +
  labs(
    title    = "ROC Curve: Logistic Regression vs Random Forest",
    subtitle = "Prediksi Exclusion Error Bansos | Test Set (30%)",
    x        = "False Positive Rate (1 - Specificity)",
    y        = "True Positive Rate (Sensitivity / Recall)",
    color    = NULL,
    caption  = "Sumber: SUSENAS Maret 2023, BPS Jawa Barat"
  )
print(p_roc)

# --- Tabel Ringkasan Performa ------------------------------------------------
tabel_performa <- tibble(
  Model       = c("Logistic Regression", "Random Forest"),
  Accuracy    = c(cm_logit$overall["Accuracy"],    cm_rf$overall["Accuracy"]),
  Sensitivity = c(cm_logit$byClass["Sensitivity"], cm_rf$byClass["Sensitivity"]),
  Specificity = c(cm_logit$byClass["Specificity"], cm_rf$byClass["Specificity"]),
  `F1-Score`  = c(cm_logit$byClass["F1"],          cm_rf$byClass["F1"]),
  AUC         = c(as.numeric(auc(roc_logit)),       as.numeric(auc(roc_rf)))
)

tabel_performa |>
  gt() |>
  tab_header(
    title    = "Perbandingan Performa Model Prediksi Exclusion Error",
    subtitle = "SUSENAS Maret 2023 — Jawa Barat | Test set 30%"
  ) |>
  fmt_percent(columns = -Model, decimals = 1) |>
  tab_style(
    style = list(cell_fill(color = "#DBEAFE"), cell_text(weight = "bold")),
    locations = cells_body(
      rows = AUC == max(AUC)
    )
  ) |>
  cols_label(
    Model       = "Model",
    Accuracy    = "Accuracy",
    Sensitivity = "Sensitivity (Recall)",
    Specificity = "Specificity",
    `F1-Score`  = "F1-Score",
    AUC         = "AUC"
  )

Perbandingan Performa Model Prediksi Exclusion Error
SUSENAS Maret 2023 — Jawa Barat | Test set 30%
Model	Accuracy	Sensitivity (Recall)	Specificity	F1-Score	AUC
Logistic Regression	60.8%	47.7%	71.5%	52.3%	65.7%
Random Forest	60.5%	47.7%	71.0%	52.1%	62.1%

1. Forest Plot Logistic Regression (Image 2) Yang meningkatkan risiko exclusion error (OR > 1, merah):

• skor_aset tinggi → RT miskin yang punya lebih banyak aset justru lebih sering tidak dapat bansos. Kemungkinan karena petugas menganggap mereka “tidak terlihat miskin” meskipun pengeluarannya di bawah GK
• Sektor Pertanian → RT miskin di sektor pertanian lebih berisiko tidak terdaftar, kemungkinan akibat database sasaran yang tidak update dan mobilitas tinggi petani
• KRT Perempuan → kepala RT perempuan lebih rentan exclusion error, konsisten dengan literatur tentang ketimpangan gender dalam akses bansos
• Pendidikan SMA/Diploma+ → RT yang KRT-nya berpendidikan lebih tinggi sering tidak dianggap prioritas meskipun tetap miskin

Yang menurunkan risiko exclusion error (OR < 1, biru):

• Wilayah Perdesaan → RT miskin di desa justru lebih terlindungi dari exclusion error dibanding perkotaan — distribusi bansos di desa lebih merata (ada BLT Desa, dll.)
• Listrik PLN & Rumah Milik → RT yang sudah terdaftar di sistem PLN atau punya rumah lebih mudah diidentifikasi oleh petugas
• Pekerja Keluarga → kelompok ini cenderung lebih “terlihat” dalam jaringan sosial desa sehingga lebih mudah terdata

2. Variable Importance Random Forest (Image 3) RF mengkonfirmasi: umur KRT adalah prediktor terpenting, jauh di atas variabel lain. Artinya:

• KRT yang lebih tua cenderung lebih dikenali sebagai kelompok rentan → lebih kecil kemungkinan exclusion error
• Jumlah anggota RT penting karena RT besar lebih “mencolok” kebutuhannya
• Status kerja & skor aset penting — menunjukkan kondisi ekonomi riil lebih berperan dari variabel perumahan

3. Performa Model (Image 4 & 5) Logistic Regression menang tipis dari Random Forest — ini agak tidak biasa karena biasanya RF lebih unggul. Artinya hubungan antara prediktor dan exclusion error di data ini cenderung linear, sehingga logit sudah cukup menangkapnya.

MetrikMakna di konteks iniAUC 65.7% (Logit) = “Acceptable”, model bisa membedakan exclusion error vs tidak, tapi tidak sempurna

• Sensitivity 47.7% = ⚠️ Model hanya mendeteksi ~48% dari RT miskin yang benar-benar tidak dapat bansos
• Specificity 71.5% = Model lebih baik mengidentifikasi yang tepat sasaran (71%) daripada yang exclusion error
• AUC < 0.7 = keduanyaWajar untuk data survei sosial, ada faktor administratif/politik yang tidak tertangkap oleh variabel SUSENAS

Kesimpulan : Model logistic regression (AUC = 65.7%) sedikit lebih unggul dari random forest (AUC = 62.1%) dalam memprediksi exclusion error bansos di Jawa Barat. Kedua model memiliki akurasi moderat (~60%), yang wajar mengingat mekanisme penetapan penerima bansos juga dipengaruhi faktor administratif yang tidak terobservasi dalam data SUSENAS. Variabel yang paling konsisten memengaruhi risiko exclusion error adalah umur KRT, jumlah anggota RT, status kerja, dan skor aset. RT miskin di perkotaan, yang dikepalai perempuan, atau yang bekerja di sektor pertanian memiliki risiko exclusion error lebih tinggi, temuan yang perlu dipertimbangkan dalam pembaruan Data Terpadu Kesejahteraan Sosial (DTKS).