Female Breadwinners
Data Preprocessing & Analisis Regresi Logistik Biner | Firth Logistic Regression
Tentang Dokumen Ini
Dokumen ini mendokumentasikan seluruh alur kerja analisis data Female Breadwinners menggunakan data Susenas 2025. Proses mencakup tiga tahap utama: (1) Data Preprocessing — mulai dari filter, join, hingga rekonstruksi variabel; (2) Regresi Logistik Biner sebagai model baseline; dan (3) Firth Logistic Regression sebagai model utama yang lebih robust terhadap kondisi data imbalance dan separation yang umum ditemui pada data kemiskinan.
📦 Load Library
Sebelum memulai analisis, seluruh paket yang dibutuhkan dimuat terlebih dahulu. Paket
dplyrmenjadi tulang punggung manipulasi data berbasis tidy syntax,foreigndigunakan untuk membaca file berformat.dbfdari hasil ekspor data survei, danreadxldigunakan untuk mengimpor data Garis Kemiskinan yang tersimpan dalam format Excel.
🔷 Bagian 1 — Female Breadwinner + KRT + Cerai Hidup/Mati
Tujuan Bagian Ini: Mengidentifikasi dan menyiapkan data individu perempuan yang secara spesifik merupakan Kepala Rumah Tangga (KRT) dengan status perkawinan cerai hidup (kode 3) atau cerai mati (kode 4). Kelompok ini disebut sebagai Female Breadwinner tipe IND1, karena definisinya lebih ketat — hanya KRT, bukan sekadar anggota rumah tangga.
1. Filter Female Breadwinner (IND1)
Prosedur: Data individu utama (
Data_ami) difilter menggunakan tiga kondisi sekaligus: jenis kelamin perempuan (r405 == 2), status sebagai KRT (r403 == 1), dan status perkawinan cerai (r404 %in% c(3,4)).Teknis: Fungsi
filter()daridplyrmengevaluasi ketiga kondisi secara bersamaan dengan operatorANDimplisit, sehingga hanya baris yang memenuhi seluruh syarat yang dipertahankan.
# Female + KRT + Cerai Hidup/Mati
library(haven)
Data_ami <- read_sav("C:/Users/User/Downloads/Susenas Maret 2025 - Permintaan Data STIS.sav")
data1 <- Data_ami %>%
filter(
r405 == 2, # Perempuan
r403 == 1, # KRT
r404 %in% c(3,4) # Cerai hidup/mati
)2. Import Data Individu (IND1)
Prosedur: File
.dbfyang merupakan data individu dari basis data Susenas diimpor menggunakan fungsiread.dbf()dari paketforeign.Teknis: Argumen
as.is = TRUEmencegah R secara otomatis mengubah kolom bertipe karakter menjadifactor, sehingga tipe data asli tetap dipertahankan dan prosesmutatelebih mudah dikontrol pada tahap berikutnya.
data_ind1 <- read.dbf(
"C:/Users/User/Downloads/bv2_2026_01_30_14_04_11_ssn202503_kor_ind1.dbf",
as.is = TRUE
)3. Filter Individu Sesuai Data IND1
Prosedur: Data individu (
data_ind1) difilter agar hanya memuat individu yang beririsan dengandata1— yakni perempuan KRT cerai. Proses ini menggunakan teknik semi-join.Teknis:
semi_join()mempertahankan baris dari tabel kiri (data_ind1) yang memiliki pasangan di tabel kanan (data1), berdasarkan kunci gabunganURUT(ID rumah tangga) danR403/r403(status KRT). Konversias.numeric()dilakukan terlebih dahulu untuk memastikan tipe data dari kedua sisi join kompatibel — penting karena file.dbfsering membaca kolom numerik sebagai karakter.
data_ind1_filtered <- data_ind1 %>%
mutate(
URUT = as.numeric(URUT),
R403 = as.numeric(R403)
) %>%
semi_join(
data1 %>%
mutate(
id_ruta = as.numeric(id_ruta),
r403 = as.numeric(r403)
),
by = c("URUT" = "id_ruta",
"R403" = "r403")
)4. Import Data Rumah Tangga
Prosedur: Data level rumah tangga diimpor dari file
.dbfterpisah. Data ini diperlukan untuk mengambil informasi jumlah anggota rumah tangga (R301) yang tidak tersedia di file data individu.Teknis: Struktur file data rumah tangga (KOR-RT) berisi satu baris per rumah tangga (
URUT), sehingga bisa dijadikan referensileft_joinke data individu menggunakan kunciURUT.
data_rt <- read.dbf(
"C:/Users/User/Downloads/bv2_2026_01_30_16_10_01_ssn202503_kor_rt.dbf",
as.is = TRUE
)5. Tambah Jumlah ART ke Level Individu
Prosedur: Informasi jumlah Anggota Rumah Tangga (
R301) dari data level rumah tangga ditambahkan ke data individu menggunakanleft_join. Kemudian dibuat variabel binerkategori_artuntuk menandai rumah tangga dengan lebih dari 4 ART.Teknis:
left_join()menyatukan dua tabel berdasarkan kunciURUT. Karena setiap individu berada dalam satu rumah tangga, join ini bersifat many-to-one (banyak individu ke satu rumah tangga). Variabelkategori_artdibuat menggunakanas.integer(R301 > 4), yang menghasilkan nilai1jika kondisi TRUE dan0jika FALSE — output langsung bertipe integer.
data_ind1_filtered <- data_ind1_filtered %>%
left_join(
data_rt %>%
mutate(URUT = as.numeric(URUT)) %>%
select(URUT, R301),
by = "URUT"
) %>%
mutate(
kategori_art = as.integer(R301 > 4)
)6. Import Blok 43 (Data Pengeluaran)
Prosedur: Data pengeluaran rumah tangga berasal dari tiga file
.dbfterpisah (Blok 43), yang mencakup berbagai kelompok komoditas (kode 11–97). Ketiga file digabung menjadi satu data frame secara vertikal.Teknis:
lapply()mengiterasi setiap path file dan menerapkanread.dbf(), menghasilkan sebuah list of data frames.bind_rows()kemudian menggabungkan semuanya secara vertikal (row-binding). Hasilnya adalah satu data frame terpadu yang merepresentasikan seluruh pengeluaran rumah tangga per unit (URUT).
files_43 <- c(
"C:/Users/User/Downloads/bv2_2026_01_30_17_45_51_ssn202503_kp_blok43_51_97.dbf",
"C:/Users/User/Downloads/bv2_2026_01_30_17_42_55_ssn202503_kp_blok43_32_36.dbf",
"C:/Users/User/Downloads/bv2_2026_01_30_17_39_45_ssn202503_kp_blok43_11_31.dbf"
)
data_kp_43 <- bind_rows(
lapply(files_43, read.dbf, as.is = TRUE)
) %>%
mutate(
R101 = as.numeric(R101),
R102 = as.numeric(R102),
URUT = as.numeric(URUT)
)7. Hitung Pengeluaran Per Kapita
Prosedur: Kolom jumlah ART (
R301) dari data rumah tangga digabungkan ke data pengeluaran, kemudian dihitung pengeluaran per kapita dengan membagi total pengeluaran rumah tangga (EXPEND) dengan jumlah ART (R301). Kode kabupaten juga dikonstruksi untuk keperluan join garis kemiskinan.Teknis: Penggabungan dilakukan dengan tiga kunci (
R101,R102,URUT) untuk menghindari ambiguitas lintas provinsi dan kabupaten. Kode kabupaten (kodekab) dikonstruksi dengan menggabungkan kode provinsi (R101) dan kode kabupaten dua digit (sprintf("%02d", R102)) menggunakanpaste0(), kemudian dikonversi ke numerik agar kompatibel dengan data Garis Kemiskinan.
data_kp_43 <- data_kp_43 %>%
left_join(
data_rt %>%
mutate(
R101 = as.numeric(R101),
R102 = as.numeric(R102),
URUT = as.numeric(URUT)
) %>%
select(R101, R102, URUT, R301),
by = c("R101", "R102", "URUT")
) %>%
mutate(
EXPEND_KAPITA = EXPEND / R301,
kodekab = as.numeric(
paste0(R101, sprintf("%02d", R102))
)
)8. Join Garis Kemiskinan & Tentukan Status Miskin
Prosedur: Data Garis Kemiskinan per kabupaten (
GK_Kab) diimpor dari Excel, kemudian digabungkan ke data pengeluaran. Status kemiskinan ditentukan dengan membandingkan pengeluaran per kapita terhadap Garis Kemiskinan kabupaten setempat.Teknis:
if_else()digunakan (bukanifelse()) karena lebih ketat dalam penanganan tipe data — memastikan output selalu bertipeinteger(1L/0L). Argumenmissing = 0Lmenangani baris di mana nilaiGKadalahNA(rumah tangga tanpa pasangan di data GK), sehingga tidak adaNAyang masuk ke variabel dependen model.
GK_Kab <- read_excel("C:/Users/User/Downloads/GK Kab.xlsx") %>%
mutate(Kode = as.numeric(Kode))
data_kp_43 <- data_kp_43 %>%
left_join(GK_Kab, by = c("kodekab" = "Kode")) %>%
mutate(
status_miskin = if_else(EXPEND_KAPITA < GK, 1L, 0L, missing = 0L)
)9. Join Status Kemiskinan ke Level Individu
Prosedur: Variabel hasil komputasi di level rumah tangga — yakni
status_miskin, Garis Kemiskinan (GK), total pengeluaran (EXPEND), dan pengeluaran per kapita (EXPEND_KAPITA) — dipindahkan kembali ke level individu melalui join.Teknis: Join tiga kunci (
R101,R102,URUT) digunakan untuk memastikan setiap individu mendapatkan status kemiskinan rumah tangganya masing-masing secara presisi. Ini adalah teknik disaggregasi data rumah tangga ke level individu yang umum digunakan dalam analisis mikro berbasis survei.
data_ind1_filtered <- data_ind1_filtered %>%
mutate(
R101 = as.numeric(R101),
R102 = as.numeric(R102),
URUT = as.numeric(URUT)
) %>%
left_join(
data_kp_43 %>%
select(R101, R102, URUT,
status_miskin, GK,
EXPEND, EXPEND_KAPITA),
by = c("R101", "R102", "URUT")
)🔷 Bagian 2 — Female Breadwinner + Cerai Hidup/Mati (IND2)
Tujuan Bagian Ini: Membuat dataset alternatif (IND2) dengan definisi Female Breadwinner yang lebih longgar: perempuan yang bercerai (hidup atau mati), tanpa syarat harus menjadi KRT. Pendekatan ini menangkap lebih banyak perempuan yang secara de facto menanggung keluarga, meskipun secara administratif mungkin tidak tercatat sebagai KRT. Bagian ini juga menangani isu duplikat yang muncul akibat definisi yang lebih luas ini.
1. Filter Female Breadwinner Cerai (IND2)
Prosedur: Dari data
Data_ami, dipilih perempuan dengan status cerai hidup atau cerai mati, tanpa pembatasan status KRT.Teknis: Filter yang lebih longgar ini menghasilkan potensi duplikat pada kolom
URUTkarena dalam satu rumah tangga bisa terdapat lebih dari satu perempuan cerai. Hal ini akan dideteksi dan ditangani pada sub-bagian berikutnya.
2. Filter Data Individu (IND2)
Prosedur: Data individu (
data_ind1) difilter menggunakansemi_joinberdasarkan empat kunci sekaligus:URUT,R403(posisi dalam rumah tangga),R404(status perkawinan), danR405(jenis kelamin). Penambahan kunciR404danR405dibanding IND1 bertujuan mengidentifikasi individu spesifik, bukan sekadar rumah tangganya.Teknis: Multi-key join meminimalkan false positive, yakni individu dari rumah tangga yang sama namun bukan perempuan cerai yang dimaksud. Semua variabel kunci dikonversi ke numerik sebelum join untuk menghindari mismatch tipe data antara file
.dbfdan tabel referensi.
data_ind2_filtered <- data_ind1 %>%
mutate(
URUT = as.numeric(URUT),
R403 = as.numeric(R403),
R404 = as.numeric(R404),
R405 = as.numeric(R405)
) %>%
semi_join(
data2 %>%
mutate(
id_ruta = as.numeric(id_ruta),
r403 = as.numeric(r403),
r404 = as.numeric(r404),
r405 = as.numeric(r405)
),
by = c("URUT" = "id_ruta",
"R403" = "r403",
"R404" = "r404",
"R405" = "r405")
)Cek Duplikat Awal
Prosedur: Setelah filter IND2 yang lebih longgar, perlu dipastikan apakah terdapat nilai
URUTyang muncul lebih dari sekali — yang mengindikasikan satu rumah tangga memiliki lebih dari satu perempuan cerai.Teknis:
count(URUT)menghitung frekuensi setiap nilaiURUT, kemudianfilter(n > 1)mengekstrak hanya yang duplikat. Output berupa tabel dan jumlah baris duplikat.
#> [1] 32# Lihat detail baris duplikat untuk inspeksi manual
data_ind2_filtered %>%
filter(URUT %in% duplikat_urut$URUT) %>%
arrange(URUT)3. Tambah Jumlah ART ke IND2
Prosedur dan Teknis: Sama seperti IND1, kolom jumlah ART (
R301) dari data rumah tangga digabungkan ke data IND2 menggunakan tiga kunci (R101,R102,URUT). Variabelkategori_artdibuat sebagai penanda rumah tangga besar (> 4 ART).
data_ind2_filtered <- data_ind2_filtered %>%
mutate(
R101 = as.numeric(R101),
R102 = as.numeric(R102),
URUT = as.numeric(URUT)
) %>%
left_join(
data_rt %>%
select(R101, R102, URUT, R301),
by = c("R101", "R102", "URUT")
) %>%
mutate(
kategori_art = as.integer(R301 > 4)
)4. Join Status Kemiskinan ke IND2
Prosedur dan Teknis: Status kemiskinan dan variabel ekonomi rumah tangga (pengeluaran, GK, bobot) digabungkan ke data IND2 menggunakan join tiga kunci. Kolom
WEIND(bobot individu) juga disertakan untuk keperluan analisis berbasis survei yang memperhitungkan desain sampling.
data_ind2_filtered <- data_ind2_filtered %>%
left_join(
data_kp_43 %>%
select(R101, R102, URUT,
status_miskin, GK,
EXPEND, EXPEND_KAPITA,
WEIND),
by = c("R101", "R102", "URUT")
)5. Penanganan Duplikat — 2 Pendekatan
Prosedur: Karena duplikat pada IND2 dapat menyebabkan pengamatan yang sama rumah tangganya terhitung lebih dari sekali dalam model regresi, perlu dipilih satu representasi per rumah tangga. Dua strategi digunakan:
Pendekatan 1 — Hapus Semua Duplikat: Membuang seluruh baris dari rumah tangga yang memiliki lebih dari satu perempuan cerai. Pendekatan ini lebih konservatif tetapi kehilangan observasi.
Pendekatan 2 — Pertahankan Umur Tertinggi: Dari setiap kelompok duplikat, dipilih satu individu dengan umur tertinggi (variabel
R407), dengan asumsi individu tertua lebih merepresentasikan peran kepala de facto rumah tangga.Teknis (Pendekatan 2):
group_by(URUT)mengelompokkan baris per rumah tangga, kemudianslice_max(order_by = R407, n = 1, with_ties = FALSE)mengambil tepat satu baris dengan nilaiR407tertinggi. Argumenwith_ties = FALSEmemastikan tidak ada seri — jika ada dua individu dengan umur sama, hanya satu yang dipertahankan (baris pertama).
# Identifikasi ulang URUT yang duplikat
duplikat_urut <- data_ind2_filtered %>%
count(URUT) %>%
filter(n > 1)
# ── PENDEKATAN 1: HAPUS SEMUA DUPLIKAT ─────────────────────────────
data_ind2_nodup <- data_ind2_filtered %>%
filter(!URUT %in% duplikat_urut$URUT)
cat("Pendekatan 1 - Hapus semua duplikat:\n")#> Pendekatan 1 - Hapus semua duplikat:#> Sebelum : 28892 baris#> Sesudah : 28828 baris#> Dihapus : 64 baris# ── PENDEKATAN 2: PERTAHANKAN UMUR TERTINGGI (R407) ────────────────
data_ind2_maxage <- data_ind2_filtered %>%
mutate(R407 = as.numeric(R407)) %>%
group_by(URUT) %>%
slice_max(order_by = R407, n = 1, with_ties = FALSE) %>%
ungroup()
cat("Pendekatan 2 - Pertahankan umur tertinggi (R407):\n")#> Pendekatan 2 - Pertahankan umur tertinggi (R407):#> Sebelum : 28892 baris#> Sesudah : 28860 baris#> Dihapus : 32 baris# ── VERIFIKASI: pastikan tidak ada duplikat tersisa ─────────────────
cat("Cek duplikat Pendekatan 1:", sum(duplicated(data_ind2_nodup$URUT)), "duplikat tersisa\n")#> Cek duplikat Pendekatan 1: 0 duplikat tersisa#> Cek duplikat Pendekatan 2: 0 duplikat tersisa🔷 Bagian 3 — Ekspansi Dataset
Tujuan Bagian Ini: Dari 3 dataset awal, dibuat total 12 dataset melalui dua lapisan filter tambahan: filter jumlah ART dan filter provinsi (Sumatera Barat). Ekspansi ini bertujuan menguji apakah temuan model bersifat konsisten di berbagai sub-populasi dan definisi operasional Female Breadwinner.
1. Skema Ekspansi 12 Dataset
Struktur ekspansi dilakukan secara bertingkat (layered) seperti pada tabel berikut:
| Layer | Keterangan |
|---|---|
| Layer 1 | 3 dataset awal (IND1 Filtered, IND2 No-Duplikat, IND2 Max-Age) |
| Layer 2 | Layer 1 + Filter R301 > 1 (minimal 2
ART) → 6 dataset |
| Layer 3 | Layer 1 & 2 + Filter R101 == 13
(Sumatera Barat) → 12 dataset |
2. Buat 12 Dataset
Prosedur: Setiap layer dataset dibuat dengan menambahkan filter menggunakan
%>% filter(). Filter ART (R301 > 1) bertujuan mengecualikan rumah tangga yang hanya terdiri dari satu anggota — yang secara konsep tidak relevan dengan peran breadwinner. Filter provinsi (R101 == 13) membatasi analisis pada Sumatera Barat sebagai studi kasus regional.Teknis: Karena semua 12 dataset merupakan subset dari dataset yang sudah ada, operasi filter cukup dijalankan satu kali. Tidak diperlukan impor data ulang. Sebuah tabel
ringkasandibuat menggunakantibble()untuk mendokumentasikan jumlah observasi (N) setiap dataset sebelum analisis dimulai.
# ── LAYER 2: + FILTER R301 > 1 (minimal 2 ART) → menjadi 6 ─────────
data_ind1_filtered_art2 <- data_ind1_filtered %>% filter(R301 > 1)
data_ind2_nodup_art2 <- data_ind2_nodup %>% filter(R301 > 1)
data_ind2_maxage_art2 <- data_ind2_maxage %>% filter(R301 > 1)
# ── LAYER 3: + FILTER R101 == 13 (Sumatera Barat) → menjadi 12 ──────
# Dari dataset awal (tanpa filter ART)
data_ind1_filtered_smb <- data_ind1_filtered %>% filter(R101 == 13)
data_ind2_nodup_smb <- data_ind2_nodup %>% filter(R101 == 13)
data_ind2_maxage_smb <- data_ind2_maxage %>% filter(R101 == 13)
# Dari dataset + filter ART
data_ind1_filtered_art2_smb <- data_ind1_filtered_art2 %>% filter(R101 == 13)
data_ind2_nodup_art2_smb <- data_ind2_nodup_art2 %>% filter(R101 == 13)
data_ind2_maxage_art2_smb <- data_ind2_maxage_art2 %>% filter(R101 == 13)
# ── RINGKASAN 12 DATASET ─────────────────────────────────────────────
ringkasan <- tibble(
Dataset = c(
"data_ind1_filtered",
"data_ind2_nodup",
"data_ind2_maxage",
"data_ind1_filtered_art2",
"data_ind2_nodup_art2",
"data_ind2_maxage_art2",
"data_ind1_filtered_smb",
"data_ind2_nodup_smb",
"data_ind2_maxage_smb",
"data_ind1_filtered_art2_smb",
"data_ind2_nodup_art2_smb",
"data_ind2_maxage_art2_smb"
),
Filter_Cerai = c("Cerai + KRT", "Cerai", "Cerai",
"Cerai + KRT", "Cerai", "Cerai",
"Cerai + KRT", "Cerai", "Cerai",
"Cerai + KRT", "Cerai", "Cerai"),
Duplikat = c("Awal", "Dihapus", "Max Umur",
"Awal", "Dihapus", "Max Umur",
"Awal", "Dihapus", "Max Umur",
"Awal", "Dihapus", "Max Umur"),
Filter_ART = c("—", "—", "—",
">1", ">1", ">1",
"—", "—", "—",
">1", ">1", ">1"),
Filter_Provinsi = c("—", "—", "—",
"—", "—", "—",
"Sumbar", "Sumbar", "Sumbar",
"Sumbar", "Sumbar", "Sumbar"),
N = c(
nrow(data_ind1_filtered),
nrow(data_ind2_nodup),
nrow(data_ind2_maxage),
nrow(data_ind1_filtered_art2),
nrow(data_ind2_nodup_art2),
nrow(data_ind2_maxage_art2),
nrow(data_ind1_filtered_smb),
nrow(data_ind2_nodup_smb),
nrow(data_ind2_maxage_smb),
nrow(data_ind1_filtered_art2_smb),
nrow(data_ind2_nodup_art2_smb),
nrow(data_ind2_maxage_art2_smb)
)
)
print(ringkasan)#> # A tibble: 12 × 6
#> Dataset Filter_Cerai Duplikat Filter_ART Filter_Provinsi N
#> <chr> <chr> <chr> <chr> <chr> <int>
#> 1 data_ind1_filtered Cerai + KRT Awal — — 27537
#> 2 data_ind2_nodup Cerai Dihapus — — 28828
#> 3 data_ind2_maxage Cerai Max Umur — — 28860
#> 4 data_ind1_filtered_ar… Cerai + KRT Awal >1 — 16035
#> 5 data_ind2_nodup_art2 Cerai Dihapus >1 — 17326
#> 6 data_ind2_maxage_art2 Cerai Max Umur >1 — 17358
#> 7 data_ind1_filtered_smb Cerai + KRT Awal — Sumbar 1135
#> 8 data_ind2_nodup_smb Cerai Dihapus — Sumbar 1195
#> 9 data_ind2_maxage_smb Cerai Max Umur — Sumbar 1196
#> 10 data_ind1_filtered_ar… Cerai + KRT Awal >1 Sumbar 677
#> 11 data_ind2_nodup_art2_… Cerai Dihapus >1 Sumbar 737
#> 12 data_ind2_maxage_art2… Cerai Max Umur >1 Sumbar 738🔷 Bagian 4 — Pembuatan Variabel
Tujuan Bagian Ini: Mengkonstruksi variabel prediktor (
X1–X11) dari kolom-kolom mentah data survei menjadi variabel kategorik biner atau ordinal yang siap dimasukkan ke dalam model regresi. Semua variabel dikonstruksi secara konsisten menggunakan fungsi terpusatbuat_variabel()sehingga dapat diterapkan ke seluruh 12 dataset dengan satu perintah.
1. Definisi Variabel Prediktor
Berikut adalah pemetaan variabel prediktor beserta sumber kolom asli dan kategorisasinya. Tanda * menunjukkan kategori referensi dalam model regresi logistik.
| Variabel | Nama | Sumber | Kategorisasi |
|---|---|---|---|
| X1 | Pendidikan | R615 | 2 = Dasar (≤SMP), 1 = Menengah (SMA/K), 0 = Tinggi* (>SMA) |
| X2 | Umur | R407 | 1 = ≥45 tahun, 0 = <45 tahun* |
| X3 | Sektor Pekerjaan | R706 | 1 = Informal, 0 = Formal* |
| X4 | Lapangan Usaha | R705 | 2 = Tersier, 1 = Sekunder, 0 = Primer* |
| X5 | Jam Kerja | R708 | 1 = <35 jam (tidak penuh), 0 = ≥35 jam (penuh)* |
| X6 | Wilayah | R105 | 1 = Perdesaan, 0 = Perkotaan* |
| X7 | Jumlah ART | R301 | 1 = >4 orang, 0 = ≤4 orang* |
| X8 | Jaminan Kesehatan | R1101_X | 1 = Tidak memiliki, 0 = Memiliki* |
| X9 | Keluhan Kesehatan | R1102 | 1 = Ada keluhan, 0 = Tidak ada keluhan* |
| X10 | Kepemilikan Rumah | R801 | 1 = Bukan milik sendiri, 0 = Milik sendiri* |
| X11 | Penggunaan Internet | R808 | 1 = Tidak menggunakan, 0 = Menggunakan* |
2. Fungsi Pembuatan Variabel
Prosedur: Fungsi
buat_variabel()menerima sebuah data frame dan mengembalikannya dengan tambahan kolomX1hinggaX11. Semua transformasi dilakukan dalam satu blokmutate()untuk efisiensi.Teknis:
case_when()digunakan untuk kategorisasi multi-level (X1, X4) yang tidak dapat diekspresikan hanya dengan satu kondisi biner.if_else()(bukanifelse()) digunakan untuk kondisi biner karena lebih strict dalam tipe data — menghindari promosi tipe yang tidak diinginkan.- Semua output dikodekan sebagai
integer(1L/0L) menggunakan suffixL, bukandouble, karena nilai kategorikal tidak membutuhkan presisi desimal dan lebih ringan secara memori.as.numeric()selalu diterapkan sebelum perbandingan untuk menghindari mismatch tipe karakter dari file.dbf.
# ============================================================
# FUNGSI PEMBUATAN VARIABEL BARU
# ============================================================
buat_variabel <- function(df) {
df %>%
mutate(
# ── 1. PENDIDIKAN (R615) ───────────────────────────────
# 2 = Dasar ≤SMP (Kode 0–9, 25)
# 1 = Menengah SMA/SMK (Kode 10–16)
# 0 = Tinggi >SMA/SMK* (Kode 17–24)
X1 = case_when(
as.numeric(R615) %in% c(0:9, 25) ~ 2L,
as.numeric(R615) %in% 10:16 ~ 1L,
as.numeric(R615) %in% 17:24 ~ 0L
),
# ── 2. UMUR (R407) ────────────────────────────────────
# 1 = ≥45 tahun
# 0 = <45 tahun*
X2 = if_else(as.numeric(R407) >= 45, 1L, 0L),
# ── 3. SEKTOR PEKERJAAN (R706) ────────────────────────
# 1 = Informal (Kode 1, 2, 5, 6)
# 0 = Formal* (Kode 3, 4)
X3 = case_when(
as.numeric(R706) %in% c(1, 2, 5, 6) ~ 1L,
as.numeric(R706) %in% c(3, 4) ~ 0L
),
# ── 4. LAPANGAN USAHA (R705) ──────────────────────────
# 2 = Tersier (Kode 1–7)
# 1 = Sekunder (Kode 8–11)
# 0 = Primer* (Kode 12–26)
X4 = case_when(
as.numeric(R705) %in% 1:7 ~ 2L,
as.numeric(R705) %in% 8:11 ~ 1L,
as.numeric(R705) %in% 12:26 ~ 0L
),
# ── 5. JAM KERJA (R708) ───────────────────────────────
# 1 = Tidak penuh waktu (<35 jam)
# 0 = Penuh waktu* (≥35 jam)
X5 = if_else(as.numeric(R708) < 35, 1L, 0L),
# ── 6. KLASIFIKASI WILAYAH (R105) ─────────────────────
# 1 = Perdesaan (Kode 2)
# 0 = Perkotaan* (Kode 1)
X6 = if_else(as.numeric(R105) == 2, 1L, 0L),
# ── 7. JUMLAH ART (R301) ──────────────────────────────
# 1 = >4 orang
# 0 = ≤4 orang*
X7 = if_else(as.numeric(R301) > 4, 1L, 0L),
# ── 8. JAMINAN KESEHATAN (R1101_X) ────────────────────
# 1 = Tidak memiliki (semua R1101 bernilai "X")
# 0 = Memiliki* (ada ≥1 R1101 bukan "X")
X8 = if_else(R1101_X == "X", 1L, 0L, missing = 0L),
# ── 9. KELUHAN KESEHATAN (R1102) ──────────────────────
# 1 = Ada keluhan (Kode 1)
# 0 = Tidak keluhan* (Kode 5)
X9 = case_when(
as.numeric(R1102) == 1 ~ 1L,
as.numeric(R1102) == 5 ~ 0L
),
# ── 10. KEPEMILIKAN RUMAH (R801) ──────────────────────
# 1 = Bukan milik sendiri (Kode 1)
# 0 = Milik sendiri* (Selain 1)
X10 = if_else(as.numeric(R801) == 1, 1L, 0L),
# ── 11. PENGGUNAAN INTERNET (R808) ────────────────────
# 1 = Tidak menggunakan internet (Kode 5)
# 0 = Menggunakan internet* (Kode 1)
X11 = case_when(
as.numeric(R808) == 5 ~ 1L,
as.numeric(R808) == 1 ~ 0L
)
)
}3. Terapkan ke Semua 12 Dataset
Prosedur: Fungsi
buat_variabel()diterapkan ke seluruh 12 dataset menggunakan pola yang konsisten. Setiap dataset dikembalikan dengan kolomX1–X11yang baru.Teknis: Setelah penerapan, dilakukan verifikasi cepat menggunakan
summary()pada kolomXuntuk memastikan distribusi nilai berada dalam rentang yang diharapkan (tidak ada nilai di luar 0, 1, atau 2) dan mengecek apakah terdapatNAberlebih.
# ── Layer 1 — Nasional tanpa filter tambahan ──────────────────────
data_ind1_filtered <- buat_variabel(data_ind1_filtered)
data_ind2_nodup <- buat_variabel(data_ind2_nodup)
data_ind2_maxage <- buat_variabel(data_ind2_maxage)
# ── Layer 2 — Nasional + R301 > 1 ────────────────────────────────
data_ind1_filtered_art2 <- buat_variabel(data_ind1_filtered_art2)
data_ind2_nodup_art2 <- buat_variabel(data_ind2_nodup_art2)
data_ind2_maxage_art2 <- buat_variabel(data_ind2_maxage_art2)
# ── Layer 3 — Sumatera Barat tanpa filter ART ─────────────────────
data_ind1_filtered_smb <- buat_variabel(data_ind1_filtered_smb)
data_ind2_nodup_smb <- buat_variabel(data_ind2_nodup_smb)
data_ind2_maxage_smb <- buat_variabel(data_ind2_maxage_smb)
# ── Layer 3 — Sumatera Barat + R301 > 1 ──────────────────────────
data_ind1_filtered_art2_smb <- buat_variabel(data_ind1_filtered_art2_smb)
data_ind2_nodup_art2_smb <- buat_variabel(data_ind2_nodup_art2_smb)
data_ind2_maxage_art2_smb <- buat_variabel(data_ind2_maxage_art2_smb)
# ── Verifikasi cepat (cek 1 dataset sebagai sampel) ───────────────
data_ind1_filtered %>%
select(starts_with("X")) %>%
summary()#> X1 X2 X3 X4
#> Min. :0.000 Min. :0.0000 Min. :0.0000 Min. :0.0000
#> 1st Qu.:2.000 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000
#> Median :2.000 Median :1.0000 Median :1.0000 Median :1.0000
#> Mean :1.672 Mean :0.8057 Mean :0.7052 Mean :0.9359
#> 3rd Qu.:2.000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:2.0000
#> Max. :2.000 Max. :1.0000 Max. :1.0000 Max. :2.0000
#> NA's :1145 NA's :1145
#> X5 X6 X7 X8
#> Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.0000
#> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000
#> Median :0.0000 Median :1.0000 Median :0.00000 Median :0.0000
#> Mean :0.4189 Mean :0.5464 Mean :0.04761 Mean :0.1787
#> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:0.0000
#> Max. :1.0000 Max. :1.0000 Max. :1.00000 Max. :1.0000
#>
#> X9 X10 X11
#> Min. :0.0000 Min. :0.0000 Min. :0.0000
#> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
#> Median :0.0000 Median :1.0000 Median :0.0000
#> Mean :0.4296 Mean :0.7082 Mean :0.4872
#> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
#> Max. :1.0000 Max. :1.0000 Max. :1.0000
#> 🔷 Bagian 5 — Model Baseline: Regresi Logistik Biner
Tujuan Bagian Ini: Membangun model baseline menggunakan Regresi Logistik Biner standar dengan estimasi Maximum Likelihood (MLE). Model ini dijalankan pada seluruh 12 dataset sebagai pembanding terhadap model Firth yang lebih robust.
Catatan Metodologis: Regresi logistik biner klasik (via
glm()) dapat menghasilkan estimasi yang biased ketika terdapat kondisi data imbalance (proporsi kelas minoritas sangat kecil) atau complete/quasi-complete separation (variabel prediktor dapat memisahkan kelas respons secara sempurna). Kedua kondisi ini umum ditemui pada data kemiskinan dengan prevalensi rendah. Model Firth (Bagian 7) hadir untuk mengatasi limitasi ini.
1. Fungsi Utama Regresi Logistik Biner
Prosedur: Fungsi
analisis_reglog()menjalankan pipeline analisis lengkap dalam satu pemanggilan: (1) Validasi data, (2) EDA distribusi Y dan X, (3) Pemodelanglm(), (4) Uji simultan G², (5) Uji parsial Wald, (6) Goodness of Fit Hosmer-Lemeshow, dan (7) Interpretasi Odds Ratio.Teknis — Uji Simultan (G²): Statistik G² dihitung sebagai selisih deviance model null dan model penuh:
G² = -2LL(null) - (-2LL(full)). Distribusi sampling G² mengikuti Chi-Kuadrat dengan derajat bebas = jumlah prediktor.Teknis — Uji Parsial (Wald): Statistik Wald dihitung sebagai
W = (β̂/SE)², yang mengikuti distribusi Chi-Kuadrat dengan df = 1 di bawah H₀: βⱼ = 0.Teknis — Hosmer-Lemeshow: Uji kesesuaian model dengan data dilakukan dengan membagi observasi ke dalam 10 desil probabilitas prediksi, kemudian membandingkan frekuensi observasi dan ekspektasi model menggunakan statistik Chi-Kuadrat. H₀: model fit, sehingga kita berharap p-value > 0.05.
Teknis — Odds Ratio:
exp(β̂)menginterpretasikan efek setiap prediktor sebagai odds ratio relatif terhadap kategori referensi. CI 95% dihitung menggunakan metode Wald:exp(β̂ ± 1.96 × SE).
# ============================================================
# FUNGSI UTAMA ANALISIS REGRESI LOGISTIK BINER
# ============================================================
library(dplyr)
library(ResourceSelection) # Hosmer-Lemeshow
library(lmtest) # lrtest (Uji Simultan)
analisis_reglog <- function(data, nama_dataset, drop_var_na = FALSE) {
cat("\n")
cat("╔══════════════════════════════════════════════════════╗\n")
cat(" DATASET:", nama_dataset, "\n")
cat("╚══════════════════════════════════════════════════════╝\n\n")
xvars <- c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10","X11")
semua_vars <- c("status_miskin", xvars)
missing_vars <- semua_vars[!semua_vars %in% names(data)]
if (length(missing_vars) > 0) {
cat("⚠️ Variabel tidak ditemukan:", paste(missing_vars, collapse = ", "), "\n")
return(invisible(NULL))
}
df <- data %>%
select(all_of(semua_vars)) %>%
mutate(across(everything(), as.numeric))
# ── Drop variabel mengandung NA (opsional) ──────────────────────
if (drop_var_na) {
var_na <- xvars[sapply(xvars, function(v) any(is.na(df[[v]])))]
if (length(var_na) > 0)
cat("⚠️ Variabel mengandung NA (dikeluarkan):",
paste(var_na, collapse = ", "), "\n\n")
xvars_model <- setdiff(xvars, var_na)
} else {
xvars_model <- xvars
}
df <- df %>%
select(all_of(c("status_miskin", xvars_model))) %>%
na.omit()
cat("── Observasi valid:", nrow(df), "baris\n\n")
# ── Cek variabel konstan ────────────────────────────────────────
konstan <- xvars_model[sapply(xvars_model, function(v) length(unique(df[[v]])) < 2)]
if (length(konstan) > 0)
cat("⚠️ Variabel KONSTAN (dikeluarkan):",
paste(konstan, collapse = ", "), "\n\n")
xvars_model <- setdiff(xvars_model, konstan)
# ── 1. EDA ─────────────────────────────────────────────────────
cat("▶ 1. EDA\n")
cat("────────────────────────────────────────────────────────\n")
tbl_y <- table(df$status_miskin)
prop_y <- prop.table(tbl_y) * 100
cat(" [Y] Distribusi status_miskin:\n")
cat(" Tidak Miskin (0):", tbl_y["0"], "obs (", round(prop_y["0"], 1), "%)\n")
cat(" Miskin (1):", tbl_y["1"], "obs (", round(prop_y["1"], 1), "%)\n\n")
cat(" [X] Distribusi & Proporsi Variabel:\n")
for (v in xvars_model) {
cat("\n ──", v, "──\n")
tab <- table(df[[v]])
prop <- prop.table(tab) * 100
for (lvl in names(tab))
cat(" ", lvl, ":", tab[lvl], "obs (", round(prop[lvl], 2), "%)\n")
}
cat("\n")
# ── 2. PEMODELAN ────────────────────────────────────────────────
cat("▶ 2. PEMODELAN REGRESI LOGISTIK BINER\n")
cat("────────────────────────────────────────────────────────\n")
formula_model <- as.formula(
paste("status_miskin ~", paste(xvars_model, collapse = " + "))
)
model <- glm(formula_model, data = df, family = binomial(link = "logit"))
model_null <- glm(status_miskin ~ 1, data = df, family = binomial(link = "logit"))
cat(" Model: logit[P(Y=1)] = β0 +",
paste(paste0("β·", xvars_model), collapse = " + "), "\n\n")
# ── 3. UJI SIMULTAN (G²) ────────────────────────────────────────
cat("▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)\n")
cat("────────────────────────────────────────────────────────\n")
G2 <- model_null$deviance - model$deviance
df_G2 <- length(xvars_model)
chi_kritis <- qchisq(0.95, df = df_G2)
pval_G2 <- pchisq(G2, df = df_G2, lower.tail = FALSE)
cat(" H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0\n")
cat(" -2 Log Likelihood (Model Null) :", round(model_null$deviance, 3), "\n")
cat(" -2 Log Likelihood (Model Penuh) :", round(model$deviance, 3), "\n")
cat(" G² hitung :", round(G2, 3), "\n")
cat(" df :", df_G2, "\n")
cat(" χ² tabel (α=5%, df=", df_G2, ") :", round(chi_kritis, 3), "\n")
cat(" p-value :", round(pval_G2, 4), "\n")
cat(ifelse(G2 > chi_kritis,
" ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan\n\n",
" ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan\n\n"))
# ── 4. UJI PARSIAL (Wald) ───────────────────────────────────────
cat("▶ 4. UJI PARSIAL (Wald Test)\n")
cat("────────────────────────────────────────────────────────\n")
s <- summary(model)
koef <- s$coefficients
chi_kritis1 <- qchisq(0.95, df = 1)
cat(" H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel =",
round(chi_kritis1, 3), ")\n\n")
cat(sprintf(" %-20s %8s %8s %10s %8s %10s\n",
"Variabel", "B", "SE", "Wald", "p-value", "Kesimpulan"))
cat(" ", paste(rep("-", 70), collapse = ""), "\n")
for (i in 1:nrow(koef)) {
wald_i <- (koef[i, 1] / koef[i, 2])^2
pval_i <- koef[i, 4]
ket <- ifelse(pval_i < 0.05, "Signifikan ✅", "Tdk Signifikan ❌")
cat(sprintf(" %-20s %8.3f %8.3f %10.3f %8.4f %s\n",
rownames(koef)[i], koef[i, 1], koef[i, 2], wald_i, pval_i, ket))
}
cat("\n")
# ── 5. GOODNESS OF FIT (Hosmer-Lemeshow) ────────────────────────
cat("▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)\n")
cat("────────────────────────────────────────────────────────\n")
hl_test <- tryCatch(
hoslem.test(df$status_miskin, fitted(model), g = 10),
error = function(e) NULL
)
if (!is.null(hl_test)) {
chi_kritis_hl <- qchisq(0.95, df = 8)
cat(" H0: Model fit dengan data\n")
cat(" χ² Hosmer-Lemeshow :", round(hl_test$statistic, 3), "\n")
cat(" df :", hl_test$parameter, "\n")
cat(" χ² tabel (df=8) :", round(chi_kritis_hl, 3), "\n")
cat(" p-value :", round(hl_test$p.value, 4), "\n")
cat(ifelse(hl_test$p.value > 0.05,
" ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data\n\n",
" ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data\n\n"))
} else {
cat(" ⚠️ Hosmer-Lemeshow tidak dapat dihitung\n\n")
}
# ── 6. INTERPRETASI PARAMETER (Odds Ratio) ──────────────────────
cat("▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))\n")
cat("────────────────────────────────────────────────────────\n")
or <- exp(coef(model))
ci <- exp(confint.default(model))
cat(sprintf(" %-20s %8s %10s %s\n",
"Variabel", "exp(B)", "95% CI", "Interpretasi Singkat"))
cat(" ", paste(rep("-", 75), collapse = ""), "\n")
for (i in seq_along(or)) {
var_nm <- names(or)[i]
if (var_nm == "(Intercept)") {
cat(sprintf(" %-20s %8.3f [%5.3f, %5.3f] Odds baseline\n",
var_nm, or[i], ci[i, 1], ci[i, 2]))
} else {
arah <- ifelse(or[i] > 1, "↑ lebih berisiko miskin", "↓ lebih kecil risiko miskin")
cat(sprintf(" %-20s %8.3f [%5.3f, %5.3f] %s\n",
var_nm, or[i], ci[i, 1], ci[i, 2], arah))
}
}
cat("\n")
cat("══════════════════════════════════════════════════════\n\n")
return(invisible(model))
}🔷 Bagian 6 — Hasil: Regresi Logistik Biner (Baseline)
Berikut ringkasan 12 dataset yang akan dianalisis, diikuti hasil model regresi logistik biner untuk setiap dataset. Perhatikan distribusi
status_miskindi setiap dataset — kondisi imbalance yang teridentifikasi di sini menjadi dasar argumentasi penggunaan model Firth pada Bagian 8.
#> # A tibble: 12 × 6
#> Dataset Filter_Cerai Duplikat Filter_ART Filter_Provinsi N
#> <chr> <chr> <chr> <chr> <chr> <int>
#> 1 data_ind1_filtered Cerai + KRT Awal — — 27537
#> 2 data_ind2_nodup Cerai Dihapus — — 28828
#> 3 data_ind2_maxage Cerai Max Umur — — 28860
#> 4 data_ind1_filtered_ar… Cerai + KRT Awal >1 — 16035
#> 5 data_ind2_nodup_art2 Cerai Dihapus >1 — 17326
#> 6 data_ind2_maxage_art2 Cerai Max Umur >1 — 17358
#> 7 data_ind1_filtered_smb Cerai + KRT Awal — Sumbar 1135
#> 8 data_ind2_nodup_smb Cerai Dihapus — Sumbar 1195
#> 9 data_ind2_maxage_smb Cerai Max Umur — Sumbar 1196
#> 10 data_ind1_filtered_ar… Cerai + KRT Awal >1 Sumbar 677
#> 11 data_ind2_nodup_art2_… Cerai Dihapus >1 Sumbar 737
#> 12 data_ind2_maxage_art2… Cerai Max Umur >1 Sumbar 738Dataset 1 — Nasional | IND1 Filtered
Dataset ini mencakup seluruh perempuan KRT cerai di tingkat nasional. Variabel yang seluruh nilainya
NA(seperti X3 pada beberapa dataset) akan otomatis dikeluarkan karena argumendrop_var_na = TRUE.
model_ind1_filtered <- analisis_reglog(data_ind1_filtered,
"Nasional | IND1 Filtered", drop_var_na = TRUE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND1 Filtered
#> ╚══════════════════════════════════════════════════════╝
#>
#> ⚠️ Variabel mengandung NA (dikeluarkan): X3, X4
#>
#> ── Observasi valid: 27537 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 26599 obs ( 96.6 %)
#> Miskin (1): 938 obs ( 3.4 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2174 obs ( 7.89 %)
#> 1 : 4692 obs ( 17.04 %)
#> 2 : 20671 obs ( 75.07 %)
#>
#> ── X2 ──
#> 0 : 5351 obs ( 19.43 %)
#> 1 : 22186 obs ( 80.57 %)
#>
#> ── X5 ──
#> 0 : 16003 obs ( 58.11 %)
#> 1 : 11534 obs ( 41.89 %)
#>
#> ── X6 ──
#> 0 : 12491 obs ( 45.36 %)
#> 1 : 15046 obs ( 54.64 %)
#>
#> ── X7 ──
#> 0 : 26226 obs ( 95.24 %)
#> 1 : 1311 obs ( 4.76 %)
#>
#> ── X8 ──
#> 0 : 22615 obs ( 82.13 %)
#> 1 : 4922 obs ( 17.87 %)
#>
#> ── X9 ──
#> 0 : 15708 obs ( 57.04 %)
#> 1 : 11829 obs ( 42.96 %)
#>
#> ── X10 ──
#> 0 : 8034 obs ( 29.18 %)
#> 1 : 19503 obs ( 70.82 %)
#>
#> ── X11 ──
#> 0 : 14122 obs ( 51.28 %)
#> 1 : 13415 obs ( 48.72 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 8183.689
#> -2 Log Likelihood (Model Penuh) : 7062.246
#> G² hitung : 1121.442
#> df : 9
#> χ² tabel (α=5%, df= 9 ) : 16.919
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.885 0.182 456.384 0.0000 Signifikan ✅
#> X1 0.496 0.086 33.515 0.0000 Signifikan ✅
#> X2 -1.256 0.083 228.375 0.0000 Signifikan ✅
#> X5 0.200 0.072 7.664 0.0056 Signifikan ✅
#> X6 0.342 0.077 19.624 0.0000 Signifikan ✅
#> X7 2.288 0.084 742.425 0.0000 Signifikan ✅
#> X8 0.083 0.087 0.916 0.3386 Tdk Signifikan ❌
#> X9 -0.418 0.076 30.493 0.0000 Signifikan ✅
#> X10 -0.401 0.093 18.727 0.0000 Signifikan ✅
#> X11 0.663 0.106 39.415 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 14.855
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.062
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.021 [0.014, 0.029] Odds baseline
#> X1 1.643 [1.389, 1.944] ↑ lebih berisiko miskin
#> X2 0.285 [0.242, 0.335] ↓ lebih kecil risiko miskin
#> X5 1.221 [1.060, 1.406] ↑ lebih berisiko miskin
#> X6 1.408 [1.210, 1.638] ↑ lebih berisiko miskin
#> X7 9.858 [8.362, 11.622] ↑ lebih berisiko miskin
#> X8 1.087 [0.916, 1.289] ↑ lebih berisiko miskin
#> X9 0.659 [0.568, 0.764] ↓ lebih kecil risiko miskin
#> X10 0.670 [0.558, 0.803] ↓ lebih kecil risiko miskin
#> X11 1.941 [1.578, 2.388] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 2 — Nasional | IND2 No-Duplikat
Dataset IND2 dengan penanganan duplikat menggunakan pendekatan penghapusan total. Cakupan observasi lebih kecil namun bersih dari potensi pseudo-replication.
model_ind2_nodup <- analisis_reglog(data_ind2_nodup,
"Nasional | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND2 No-Duplikat
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 27617 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 26613 obs ( 96.4 %)
#> Miskin (1): 1004 obs ( 3.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2358 obs ( 8.54 %)
#> 1 : 4933 obs ( 17.86 %)
#> 2 : 20326 obs ( 73.6 %)
#>
#> ── X2 ──
#> 0 : 5971 obs ( 21.62 %)
#> 1 : 21646 obs ( 78.38 %)
#>
#> ── X3 ──
#> 0 : 8437 obs ( 30.55 %)
#> 1 : 19180 obs ( 69.45 %)
#>
#> ── X4 ──
#> 0 : 13547 obs ( 49.05 %)
#> 1 : 2817 obs ( 10.2 %)
#> 2 : 11253 obs ( 40.75 %)
#>
#> ── X5 ──
#> 0 : 16911 obs ( 61.23 %)
#> 1 : 10706 obs ( 38.77 %)
#>
#> ── X6 ──
#> 0 : 12724 obs ( 46.07 %)
#> 1 : 14893 obs ( 53.93 %)
#>
#> ── X7 ──
#> 0 : 25988 obs ( 94.1 %)
#> 1 : 1629 obs ( 5.9 %)
#>
#> ── X8 ──
#> 0 : 22671 obs ( 82.09 %)
#> 1 : 4946 obs ( 17.91 %)
#>
#> ── X9 ──
#> 0 : 16056 obs ( 58.14 %)
#> 1 : 11561 obs ( 41.86 %)
#>
#> ── X10 ──
#> 0 : 7686 obs ( 27.83 %)
#> 1 : 19931 obs ( 72.17 %)
#>
#> ── X11 ──
#> 0 : 14711 obs ( 53.27 %)
#> 1 : 12906 obs ( 46.73 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 8626.444
#> -2 Log Likelihood (Model Penuh) : 7508.501
#> G² hitung : 1117.943
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.928 0.170 534.864 0.0000 Signifikan ✅
#> X1 0.489 0.080 37.821 0.0000 Signifikan ✅
#> X2 -1.234 0.080 236.275 0.0000 Signifikan ✅
#> X3 0.043 0.083 0.272 0.6022 Tdk Signifikan ❌
#> X4 0.068 0.044 2.432 0.1189 Tdk Signifikan ❌
#> X5 0.204 0.071 8.221 0.0041 Signifikan ✅
#> X6 0.261 0.080 10.532 0.0012 Signifikan ✅
#> X7 2.108 0.080 694.184 0.0000 Signifikan ✅
#> X8 0.060 0.084 0.518 0.4718 Tdk Signifikan ❌
#> X9 -0.370 0.073 25.510 0.0000 Signifikan ✅
#> X10 -0.344 0.091 14.327 0.0002 Signifikan ✅
#> X11 0.674 0.103 42.511 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 19.773
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.0112
#> ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.020 [0.014, 0.027] Odds baseline
#> X1 1.631 [1.395, 1.906] ↑ lebih berisiko miskin
#> X2 0.291 [0.249, 0.341] ↓ lebih kecil risiko miskin
#> X3 1.044 [0.888, 1.228] ↑ lebih berisiko miskin
#> X4 1.070 [0.983, 1.166] ↑ lebih berisiko miskin
#> X5 1.226 [1.067, 1.409] ↑ lebih berisiko miskin
#> X6 1.298 [1.109, 1.520] ↑ lebih berisiko miskin
#> X7 8.228 [7.034, 9.624] ↑ lebih berisiko miskin
#> X8 1.062 [0.901, 1.252] ↑ lebih berisiko miskin
#> X9 0.691 [0.599, 0.797] ↓ lebih kecil risiko miskin
#> X10 0.709 [0.593, 0.847] ↓ lebih kecil risiko miskin
#> X11 1.963 [1.603, 2.404] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 3 — Nasional | IND2 Max-Age
Dataset IND2 dengan penanganan duplikat menggunakan seleksi individu tertua per rumah tangga. Mempertahankan lebih banyak observasi dibanding pendekatan penghapusan.
model_ind2_maxage <- analisis_reglog(data_ind2_maxage,
"Nasional | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND2 Max-Age
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 27642 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 26633 obs ( 96.3 %)
#> Miskin (1): 1009 obs ( 3.7 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2359 obs ( 8.53 %)
#> 1 : 4937 obs ( 17.86 %)
#> 2 : 20346 obs ( 73.61 %)
#>
#> ── X2 ──
#> 0 : 5980 obs ( 21.63 %)
#> 1 : 21662 obs ( 78.37 %)
#>
#> ── X3 ──
#> 0 : 8450 obs ( 30.57 %)
#> 1 : 19192 obs ( 69.43 %)
#>
#> ── X4 ──
#> 0 : 13559 obs ( 49.05 %)
#> 1 : 2824 obs ( 10.22 %)
#> 2 : 11259 obs ( 40.73 %)
#>
#> ── X5 ──
#> 0 : 16926 obs ( 61.23 %)
#> 1 : 10716 obs ( 38.77 %)
#>
#> ── X6 ──
#> 0 : 12736 obs ( 46.07 %)
#> 1 : 14906 obs ( 53.93 %)
#>
#> ── X7 ──
#> 0 : 26002 obs ( 94.07 %)
#> 1 : 1640 obs ( 5.93 %)
#>
#> ── X8 ──
#> 0 : 22691 obs ( 82.09 %)
#> 1 : 4951 obs ( 17.91 %)
#>
#> ── X9 ──
#> 0 : 16073 obs ( 58.15 %)
#> 1 : 11569 obs ( 41.85 %)
#>
#> ── X10 ──
#> 0 : 7694 obs ( 27.83 %)
#> 1 : 19948 obs ( 72.17 %)
#>
#> ── X11 ──
#> 0 : 14728 obs ( 53.28 %)
#> 1 : 12914 obs ( 46.72 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 8661.053
#> -2 Log Likelihood (Model Penuh) : 7535.555
#> G² hitung : 1125.497
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.924 0.170 533.704 0.0000 Signifikan ✅
#> X1 0.495 0.080 38.723 0.0000 Signifikan ✅
#> X2 -1.233 0.080 236.547 0.0000 Signifikan ✅
#> X3 0.038 0.082 0.213 0.6442 Tdk Signifikan ❌
#> X4 0.071 0.043 2.673 0.1020 Tdk Signifikan ❌
#> X5 0.197 0.071 7.710 0.0055 Signifikan ✅
#> X6 0.246 0.080 9.443 0.0021 Signifikan ✅
#> X7 2.110 0.080 700.525 0.0000 Signifikan ✅
#> X8 0.066 0.084 0.616 0.4324 Tdk Signifikan ❌
#> X9 -0.367 0.073 25.232 0.0000 Signifikan ✅
#> X10 -0.350 0.091 14.864 0.0001 Signifikan ✅
#> X11 0.676 0.103 42.786 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 19.91
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.0107
#> ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.020 [0.014, 0.028] Odds baseline
#> X1 1.640 [1.403, 1.917] ↑ lebih berisiko miskin
#> X2 0.291 [0.249, 0.341] ↓ lebih kecil risiko miskin
#> X3 1.039 [0.884, 1.221] ↑ lebih berisiko miskin
#> X4 1.074 [0.986, 1.169] ↑ lebih berisiko miskin
#> X5 1.217 [1.060, 1.398] ↑ lebih berisiko miskin
#> X6 1.279 [1.093, 1.497] ↑ lebih berisiko miskin
#> X7 8.247 [7.054, 9.641] ↑ lebih berisiko miskin
#> X8 1.068 [0.907, 1.258] ↑ lebih berisiko miskin
#> X9 0.693 [0.601, 0.800] ↓ lebih kecil risiko miskin
#> X10 0.705 [0.590, 0.842] ↓ lebih kecil risiko miskin
#> X11 1.966 [1.606, 2.408] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 4 — Nasional + ART>1 | IND1 Filtered
Subset dari Dataset 1 yang hanya menyertakan rumah tangga dengan minimal 2 ART, untuk mengecualikan rumah tangga “tunggal” yang memiliki karakteristik berbeda.
model_ind1_filtered_art2 <- analisis_reglog(data_ind1_filtered_art2,
"Nasional + ART>1 | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND1 Filtered
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 15421 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 14558 obs ( 94.4 %)
#> Miskin (1): 863 obs ( 5.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1595 obs ( 10.34 %)
#> 1 : 3457 obs ( 22.42 %)
#> 2 : 10369 obs ( 67.24 %)
#>
#> ── X2 ──
#> 0 : 4649 obs ( 30.15 %)
#> 1 : 10772 obs ( 69.85 %)
#>
#> ── X3 ──
#> 0 : 5194 obs ( 33.68 %)
#> 1 : 10227 obs ( 66.32 %)
#>
#> ── X4 ──
#> 0 : 8185 obs ( 53.08 %)
#> 1 : 1597 obs ( 10.36 %)
#> 2 : 5639 obs ( 36.57 %)
#>
#> ── X5 ──
#> 0 : 10418 obs ( 67.56 %)
#> 1 : 5003 obs ( 32.44 %)
#>
#> ── X6 ──
#> 0 : 7434 obs ( 48.21 %)
#> 1 : 7987 obs ( 51.79 %)
#>
#> ── X7 ──
#> 0 : 14160 obs ( 91.82 %)
#> 1 : 1261 obs ( 8.18 %)
#>
#> ── X8 ──
#> 0 : 12770 obs ( 82.81 %)
#> 1 : 2651 obs ( 17.19 %)
#>
#> ── X9 ──
#> 0 : 9716 obs ( 63 %)
#> 1 : 5705 obs ( 37 %)
#>
#> ── X10 ──
#> 0 : 3107 obs ( 20.15 %)
#> 1 : 12314 obs ( 79.85 %)
#>
#> ── X11 ──
#> 0 : 9806 obs ( 63.59 %)
#> 1 : 5615 obs ( 36.41 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 6652.957
#> -2 Log Likelihood (Model Penuh) : 5741.038
#> G² hitung : 911.919
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.699 0.188 386.537 0.0000 Signifikan ✅
#> X1 0.496 0.088 31.728 0.0000 Signifikan ✅
#> X2 -0.895 0.083 115.568 0.0000 Signifikan ✅
#> X3 -0.007 0.092 0.006 0.9385 Tdk Signifikan ❌
#> X4 0.099 0.048 4.211 0.0402 Signifikan ✅
#> X5 0.338 0.077 19.317 0.0000 Signifikan ✅
#> X6 0.247 0.091 7.370 0.0066 Signifikan ✅
#> X7 1.807 0.086 445.193 0.0000 Signifikan ✅
#> X8 0.099 0.093 1.139 0.2859 Tdk Signifikan ❌
#> X9 -0.264 0.081 10.731 0.0011 Signifikan ✅
#> X10 -0.501 0.100 25.270 0.0000 Signifikan ✅
#> X11 0.757 0.110 47.781 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 9.555
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.2976
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [0.017, 0.036] Odds baseline
#> X1 1.643 [1.382, 1.953] ↑ lebih berisiko miskin
#> X2 0.409 [0.347, 0.481] ↓ lebih kecil risiko miskin
#> X3 0.993 [0.829, 1.190] ↓ lebih kecil risiko miskin
#> X4 1.104 [1.004, 1.214] ↑ lebih berisiko miskin
#> X5 1.402 [1.206, 1.630] ↑ lebih berisiko miskin
#> X6 1.280 [1.071, 1.530] ↑ lebih berisiko miskin
#> X7 6.092 [5.151, 7.206] ↑ lebih berisiko miskin
#> X8 1.105 [0.920, 1.326] ↑ lebih berisiko miskin
#> X9 0.768 [0.656, 0.899] ↓ lebih kecil risiko miskin
#> X10 0.606 [0.498, 0.737] ↓ lebih kecil risiko miskin
#> X11 2.133 [1.720, 2.643] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 5 — Nasional + ART>1 | IND2 No-Duplikat
model_ind2_nodup_art2 <- analisis_reglog(data_ind2_nodup_art2,
"Nasional + ART>1 | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND2 No-Duplikat
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 16646 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 15676 obs ( 94.2 %)
#> Miskin (1): 970 obs ( 5.8 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1840 obs ( 11.05 %)
#> 1 : 3816 obs ( 22.92 %)
#> 2 : 10990 obs ( 66.02 %)
#>
#> ── X2 ──
#> 0 : 5336 obs ( 32.06 %)
#> 1 : 11310 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 5850 obs ( 35.14 %)
#> 1 : 10796 obs ( 64.86 %)
#>
#> ── X4 ──
#> 0 : 9011 obs ( 54.13 %)
#> 1 : 1772 obs ( 10.65 %)
#> 2 : 5863 obs ( 35.22 %)
#>
#> ── X5 ──
#> 0 : 11326 obs ( 68.04 %)
#> 1 : 5320 obs ( 31.96 %)
#>
#> ── X6 ──
#> 0 : 8162 obs ( 49.03 %)
#> 1 : 8484 obs ( 50.97 %)
#>
#> ── X7 ──
#> 0 : 15017 obs ( 90.21 %)
#> 1 : 1629 obs ( 9.79 %)
#>
#> ── X8 ──
#> 0 : 13773 obs ( 82.74 %)
#> 1 : 2873 obs ( 17.26 %)
#>
#> ── X9 ──
#> 0 : 10616 obs ( 63.78 %)
#> 1 : 6030 obs ( 36.22 %)
#>
#> ── X10 ──
#> 0 : 3282 obs ( 19.72 %)
#> 1 : 13364 obs ( 80.28 %)
#>
#> ── X11 ──
#> 0 : 10753 obs ( 64.6 %)
#> 1 : 5893 obs ( 35.4 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 7397.045
#> -2 Log Likelihood (Model Penuh) : 6461.919
#> G² hitung : 935.126
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.686 0.171 466.758 0.0000 Signifikan ✅
#> X1 0.500 0.079 40.184 0.0000 Signifikan ✅
#> X2 -0.901 0.079 130.132 0.0000 Signifikan ✅
#> X3 0.002 0.086 0.001 0.9773 Tdk Signifikan ❌
#> X4 0.077 0.045 2.916 0.0877 Tdk Signifikan ❌
#> X5 0.354 0.073 23.650 0.0000 Signifikan ✅
#> X6 0.270 0.084 10.349 0.0013 Signifikan ✅
#> X7 1.643 0.079 429.606 0.0000 Signifikan ✅
#> X8 0.059 0.088 0.454 0.5007 Tdk Signifikan ❌
#> X9 -0.247 0.076 10.574 0.0011 Signifikan ✅
#> X10 -0.458 0.095 23.253 0.0000 Signifikan ✅
#> X11 0.781 0.104 56.720 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 7.728
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.4605
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [0.018, 0.035] Odds baseline
#> X1 1.649 [1.413, 1.925] ↑ lebih berisiko miskin
#> X2 0.406 [0.348, 0.474] ↓ lebih kecil risiko miskin
#> X3 1.002 [0.848, 1.185] ↑ lebih berisiko miskin
#> X4 1.080 [0.989, 1.180] ↑ lebih berisiko miskin
#> X5 1.424 [1.235, 1.642] ↑ lebih berisiko miskin
#> X6 1.310 [1.111, 1.544] ↑ lebih berisiko miskin
#> X7 5.171 [4.427, 6.041] ↑ lebih berisiko miskin
#> X8 1.061 [0.893, 1.260] ↑ lebih berisiko miskin
#> X9 0.781 [0.673, 0.906] ↓ lebih kecil risiko miskin
#> X10 0.633 [0.525, 0.762] ↓ lebih kecil risiko miskin
#> X11 2.183 [1.782, 2.675] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 6 — Nasional + ART>1 | IND2 Max-Age
model_ind2_maxage_art2 <- analisis_reglog(data_ind2_maxage_art2,
"Nasional + ART>1 | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND2 Max-Age
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 16671 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 15696 obs ( 94.2 %)
#> Miskin (1): 975 obs ( 5.8 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1841 obs ( 11.04 %)
#> 1 : 3820 obs ( 22.91 %)
#> 2 : 11010 obs ( 66.04 %)
#>
#> ── X2 ──
#> 0 : 5345 obs ( 32.06 %)
#> 1 : 11326 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 5863 obs ( 35.17 %)
#> 1 : 10808 obs ( 64.83 %)
#>
#> ── X4 ──
#> 0 : 9023 obs ( 54.12 %)
#> 1 : 1779 obs ( 10.67 %)
#> 2 : 5869 obs ( 35.2 %)
#>
#> ── X5 ──
#> 0 : 11341 obs ( 68.03 %)
#> 1 : 5330 obs ( 31.97 %)
#>
#> ── X6 ──
#> 0 : 8174 obs ( 49.03 %)
#> 1 : 8497 obs ( 50.97 %)
#>
#> ── X7 ──
#> 0 : 15031 obs ( 90.16 %)
#> 1 : 1640 obs ( 9.84 %)
#>
#> ── X8 ──
#> 0 : 13793 obs ( 82.74 %)
#> 1 : 2878 obs ( 17.26 %)
#>
#> ── X9 ──
#> 0 : 10633 obs ( 63.78 %)
#> 1 : 6038 obs ( 36.22 %)
#>
#> ── X10 ──
#> 0 : 3290 obs ( 19.73 %)
#> 1 : 13381 obs ( 80.27 %)
#>
#> ── X11 ──
#> 0 : 10770 obs ( 64.6 %)
#> 1 : 5901 obs ( 35.4 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 7427.86
#> -2 Log Likelihood (Model Penuh) : 6486.52
#> G² hitung : 941.34
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -3.681 0.171 465.742 0.0000 Signifikan ✅
#> X1 0.506 0.079 41.125 0.0000 Signifikan ✅
#> X2 -0.900 0.079 130.470 0.0000 Signifikan ✅
#> X3 -0.003 0.085 0.001 0.9763 Tdk Signifikan ❌
#> X4 0.081 0.045 3.208 0.0733 Tdk Signifikan ❌
#> X5 0.347 0.073 22.901 0.0000 Signifikan ✅
#> X6 0.254 0.084 9.253 0.0024 Signifikan ✅
#> X7 1.646 0.079 434.039 0.0000 Signifikan ✅
#> X8 0.065 0.087 0.557 0.4556 Tdk Signifikan ❌
#> X9 -0.243 0.076 10.307 0.0013 Signifikan ✅
#> X10 -0.464 0.095 23.927 0.0000 Signifikan ✅
#> X11 0.782 0.104 57.013 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 6.627
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.5774
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [0.018, 0.035] Odds baseline
#> X1 1.659 [1.421, 1.936] ↑ lebih berisiko miskin
#> X2 0.406 [0.348, 0.474] ↓ lebih kecil risiko miskin
#> X3 0.997 [0.844, 1.179] ↓ lebih kecil risiko miskin
#> X4 1.084 [0.992, 1.184] ↑ lebih berisiko miskin
#> X5 1.415 [1.227, 1.631] ↑ lebih berisiko miskin
#> X6 1.290 [1.095, 1.519] ↑ lebih berisiko miskin
#> X7 5.184 [4.441, 6.053] ↑ lebih berisiko miskin
#> X8 1.067 [0.900, 1.266] ↑ lebih berisiko miskin
#> X9 0.784 [0.676, 0.910] ↓ lebih kecil risiko miskin
#> X10 0.629 [0.522, 0.757] ↓ lebih kecil risiko miskin
#> X11 2.186 [1.785, 2.679] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 7 — Sumatera Barat | IND1 Filtered
Analisis dimulai pada sub-populasi Sumatera Barat. Ukuran sampel yang lebih kecil membuat risiko separation meningkat, sehingga perbandingan dengan model Firth menjadi semakin relevan.
model_ind1_filtered_smb <- analisis_reglog(data_ind1_filtered_smb,
"Sumatera Barat | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND1 Filtered
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1069 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1057 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 101 obs ( 9.45 %)
#> 1 : 266 obs ( 24.88 %)
#> 2 : 702 obs ( 65.67 %)
#>
#> ── X2 ──
#> 0 : 210 obs ( 19.64 %)
#> 1 : 859 obs ( 80.36 %)
#>
#> ── X3 ──
#> 0 : 287 obs ( 26.85 %)
#> 1 : 782 obs ( 73.15 %)
#>
#> ── X4 ──
#> 0 : 487 obs ( 45.56 %)
#> 1 : 93 obs ( 8.7 %)
#> 2 : 489 obs ( 45.74 %)
#>
#> ── X5 ──
#> 0 : 662 obs ( 61.93 %)
#> 1 : 407 obs ( 38.07 %)
#>
#> ── X6 ──
#> 0 : 450 obs ( 42.1 %)
#> 1 : 619 obs ( 57.9 %)
#>
#> ── X7 ──
#> 0 : 979 obs ( 91.58 %)
#> 1 : 90 obs ( 8.42 %)
#>
#> ── X8 ──
#> 0 : 864 obs ( 80.82 %)
#> 1 : 205 obs ( 19.18 %)
#>
#> ── X9 ──
#> 0 : 613 obs ( 57.34 %)
#> 1 : 456 obs ( 42.66 %)
#>
#> ── X10 ──
#> 0 : 176 obs ( 16.46 %)
#> 1 : 893 obs ( 83.54 %)
#>
#> ── X11 ──
#> 0 : 608 obs ( 56.88 %)
#> 1 : 461 obs ( 43.12 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 131.615
#> -2 Log Likelihood (Model Penuh) : 103.216
#> G² hitung : 28.399
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.0028
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -5.268 1.404 14.087 0.0002 Signifikan ✅
#> X1 0.098 0.631 0.024 0.8768 Tdk Signifikan ❌
#> X2 -2.648 0.798 10.996 0.0009 Signifikan ✅
#> X3 0.070 0.757 0.008 0.9268 Tdk Signifikan ❌
#> X4 0.219 0.412 0.282 0.5956 Tdk Signifikan ❌
#> X5 0.848 0.651 1.697 0.1926 Tdk Signifikan ❌
#> X6 0.037 0.712 0.003 0.9589 Tdk Signifikan ❌
#> X7 2.208 0.681 10.498 0.0012 Signifikan ✅
#> X8 0.518 0.670 0.598 0.4392 Tdk Signifikan ❌
#> X9 -1.103 0.817 1.820 0.1773 Tdk Signifikan ❌
#> X10 0.406 0.918 0.195 0.6584 Tdk Signifikan ❌
#> X11 1.866 0.955 3.815 0.0508 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 10.918
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.2064
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.005 [0.000, 0.081] Odds baseline
#> X1 1.103 [0.320, 3.797] ↑ lebih berisiko miskin
#> X2 0.071 [0.015, 0.339] ↓ lebih kecil risiko miskin
#> X3 1.072 [0.243, 4.731] ↑ lebih berisiko miskin
#> X4 1.244 [0.555, 2.790] ↑ lebih berisiko miskin
#> X5 2.335 [0.652, 8.362] ↑ lebih berisiko miskin
#> X6 1.037 [0.257, 4.191] ↑ lebih berisiko miskin
#> X7 9.097 [2.392, 34.589] ↑ lebih berisiko miskin
#> X8 1.679 [0.452, 6.238] ↑ lebih berisiko miskin
#> X9 0.332 [0.067, 1.647] ↓ lebih kecil risiko miskin
#> X10 1.501 [0.248, 9.068] ↑ lebih berisiko miskin
#> X11 6.460 [0.994, 41.992] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 8 — Sumatera Barat | IND2 No-Duplikat
model_ind2_nodup_smb <- analisis_reglog(data_ind2_nodup_smb,
"Sumatera Barat | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND2 No-Duplikat
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1124 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1112 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 116 obs ( 10.32 %)
#> 1 : 283 obs ( 25.18 %)
#> 2 : 725 obs ( 64.5 %)
#>
#> ── X2 ──
#> 0 : 235 obs ( 20.91 %)
#> 1 : 889 obs ( 79.09 %)
#>
#> ── X3 ──
#> 0 : 310 obs ( 27.58 %)
#> 1 : 814 obs ( 72.42 %)
#>
#> ── X4 ──
#> 0 : 525 obs ( 46.71 %)
#> 1 : 98 obs ( 8.72 %)
#> 2 : 501 obs ( 44.57 %)
#>
#> ── X5 ──
#> 0 : 699 obs ( 62.19 %)
#> 1 : 425 obs ( 37.81 %)
#>
#> ── X6 ──
#> 0 : 483 obs ( 42.97 %)
#> 1 : 641 obs ( 57.03 %)
#>
#> ── X7 ──
#> 0 : 1017 obs ( 90.48 %)
#> 1 : 107 obs ( 9.52 %)
#>
#> ── X8 ──
#> 0 : 906 obs ( 80.6 %)
#> 1 : 218 obs ( 19.4 %)
#>
#> ── X9 ──
#> 0 : 651 obs ( 57.92 %)
#> 1 : 473 obs ( 42.08 %)
#>
#> ── X10 ──
#> 0 : 180 obs ( 16.01 %)
#> 1 : 944 obs ( 83.99 %)
#>
#> ── X11 ──
#> 0 : 653 obs ( 58.1 %)
#> 1 : 471 obs ( 41.9 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 132.825
#> -2 Log Likelihood (Model Penuh) : 106.157
#> G² hitung : 26.668
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.0052
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -5.509 1.408 15.307 0.0001 Signifikan ✅
#> X1 0.297 0.628 0.223 0.6364 Tdk Signifikan ❌
#> X2 -2.618 0.805 10.591 0.0011 Signifikan ✅
#> X3 0.121 0.744 0.026 0.8709 Tdk Signifikan ❌
#> X4 0.185 0.397 0.217 0.6415 Tdk Signifikan ❌
#> X5 0.775 0.647 1.434 0.2312 Tdk Signifikan ❌
#> X6 -0.005 0.673 0.000 0.9942 Tdk Signifikan ❌
#> X7 1.984 0.676 8.618 0.0033 Signifikan ✅
#> X8 0.455 0.660 0.475 0.4907 Tdk Signifikan ❌
#> X9 -1.065 0.815 1.709 0.1911 Tdk Signifikan ❌
#> X10 0.352 0.913 0.149 0.6997 Tdk Signifikan ❌
#> X11 1.833 0.957 3.670 0.0554 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 4.779
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.7809
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.004 [0.000, 0.064] Odds baseline
#> X1 1.346 [0.393, 4.612] ↑ lebih berisiko miskin
#> X2 0.073 [0.015, 0.353] ↓ lebih kecil risiko miskin
#> X3 1.129 [0.262, 4.854] ↑ lebih berisiko miskin
#> X4 1.203 [0.552, 2.621] ↑ lebih berisiko miskin
#> X5 2.171 [0.610, 7.724] ↑ lebih berisiko miskin
#> X6 0.995 [0.266, 3.722] ↓ lebih kecil risiko miskin
#> X7 7.275 [1.934, 27.369] ↑ lebih berisiko miskin
#> X8 1.576 [0.432, 5.745] ↑ lebih berisiko miskin
#> X9 0.345 [0.070, 1.702] ↓ lebih kecil risiko miskin
#> X10 1.422 [0.237, 8.522] ↑ lebih berisiko miskin
#> X11 6.252 [0.959, 40.786] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 9 — Sumatera Barat | IND2 Max-Age
model_ind2_maxage_smb <- analisis_reglog(data_ind2_maxage_smb,
"Sumatera Barat | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND2 Max-Age
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1125 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1113 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 116 obs ( 10.31 %)
#> 1 : 283 obs ( 25.16 %)
#> 2 : 726 obs ( 64.53 %)
#>
#> ── X2 ──
#> 0 : 235 obs ( 20.89 %)
#> 1 : 890 obs ( 79.11 %)
#>
#> ── X3 ──
#> 0 : 311 obs ( 27.64 %)
#> 1 : 814 obs ( 72.36 %)
#>
#> ── X4 ──
#> 0 : 525 obs ( 46.67 %)
#> 1 : 99 obs ( 8.8 %)
#> 2 : 501 obs ( 44.53 %)
#>
#> ── X5 ──
#> 0 : 700 obs ( 62.22 %)
#> 1 : 425 obs ( 37.78 %)
#>
#> ── X6 ──
#> 0 : 483 obs ( 42.93 %)
#> 1 : 642 obs ( 57.07 %)
#>
#> ── X7 ──
#> 0 : 1018 obs ( 90.49 %)
#> 1 : 107 obs ( 9.51 %)
#>
#> ── X8 ──
#> 0 : 907 obs ( 80.62 %)
#> 1 : 218 obs ( 19.38 %)
#>
#> ── X9 ──
#> 0 : 651 obs ( 57.87 %)
#> 1 : 474 obs ( 42.13 %)
#>
#> ── X10 ──
#> 0 : 180 obs ( 16 %)
#> 1 : 945 obs ( 84 %)
#>
#> ── X11 ──
#> 0 : 654 obs ( 58.13 %)
#> 1 : 471 obs ( 41.87 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 132.847
#> -2 Log Likelihood (Model Penuh) : 106.158
#> G² hitung : 26.689
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.0051
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -5.509 1.408 15.307 0.0001 Signifikan ✅
#> X1 0.297 0.628 0.223 0.6365 Tdk Signifikan ❌
#> X2 -2.619 0.805 10.594 0.0011 Signifikan ✅
#> X3 0.121 0.744 0.026 0.8707 Tdk Signifikan ❌
#> X4 0.185 0.397 0.217 0.6415 Tdk Signifikan ❌
#> X5 0.775 0.647 1.434 0.2311 Tdk Signifikan ❌
#> X6 -0.005 0.673 0.000 0.9942 Tdk Signifikan ❌
#> X7 1.985 0.676 8.619 0.0033 Signifikan ✅
#> X8 0.455 0.660 0.475 0.4907 Tdk Signifikan ❌
#> X9 -1.066 0.815 1.710 0.1910 Tdk Signifikan ❌
#> X10 0.352 0.913 0.149 0.6997 Tdk Signifikan ❌
#> X11 1.833 0.957 3.671 0.0554 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 4.78
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.7808
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.004 [0.000, 0.064] Odds baseline
#> X1 1.346 [0.393, 4.612] ↑ lebih berisiko miskin
#> X2 0.073 [0.015, 0.353] ↓ lebih kecil risiko miskin
#> X3 1.129 [0.262, 4.855] ↑ lebih berisiko miskin
#> X4 1.203 [0.552, 2.621] ↑ lebih berisiko miskin
#> X5 2.171 [0.610, 7.725] ↑ lebih berisiko miskin
#> X6 0.995 [0.266, 3.722] ↓ lebih kecil risiko miskin
#> X7 7.276 [1.934, 27.371] ↑ lebih berisiko miskin
#> X8 1.576 [0.432, 5.745] ↑ lebih berisiko miskin
#> X9 0.345 [0.070, 1.702] ↓ lebih kecil risiko miskin
#> X10 1.422 [0.237, 8.522] ↑ lebih berisiko miskin
#> X11 6.254 [0.959, 40.790] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 10 — Sumbar + ART>1 | IND1 Filtered
model_ind1_filtered_art2_smb <- analisis_reglog(data_ind1_filtered_art2_smb,
"Sumbar + ART>1 | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND1 Filtered
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 650 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 639 obs ( 98.3 %)
#> Miskin (1): 11 obs ( 1.7 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 80 obs ( 12.31 %)
#> 1 : 209 obs ( 32.15 %)
#> 2 : 361 obs ( 55.54 %)
#>
#> ── X2 ──
#> 0 : 201 obs ( 30.92 %)
#> 1 : 449 obs ( 69.08 %)
#>
#> ── X3 ──
#> 0 : 212 obs ( 32.62 %)
#> 1 : 438 obs ( 67.38 %)
#>
#> ── X4 ──
#> 0 : 349 obs ( 53.69 %)
#> 1 : 65 obs ( 10 %)
#> 2 : 236 obs ( 36.31 %)
#>
#> ── X5 ──
#> 0 : 469 obs ( 72.15 %)
#> 1 : 181 obs ( 27.85 %)
#>
#> ── X6 ──
#> 0 : 305 obs ( 46.92 %)
#> 1 : 345 obs ( 53.08 %)
#>
#> ── X7 ──
#> 0 : 560 obs ( 86.15 %)
#> 1 : 90 obs ( 13.85 %)
#>
#> ── X8 ──
#> 0 : 543 obs ( 83.54 %)
#> 1 : 107 obs ( 16.46 %)
#>
#> ── X9 ──
#> 0 : 400 obs ( 61.54 %)
#> 1 : 250 obs ( 38.46 %)
#>
#> ── X10 ──
#> 0 : 65 obs ( 10 %)
#> 1 : 585 obs ( 90 %)
#>
#> ── X11 ──
#> 0 : 458 obs ( 70.46 %)
#> 1 : 192 obs ( 29.54 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 111.552
#> -2 Log Likelihood (Model Penuh) : 89.867
#> G² hitung : 21.685
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.0269
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -4.774 1.450 10.839 0.0010 Signifikan ✅
#> X1 0.145 0.626 0.054 0.8166 Tdk Signifikan ❌
#> X2 -2.264 0.800 8.010 0.0047 Signifikan ✅
#> X3 0.203 0.770 0.070 0.7916 Tdk Signifikan ❌
#> X4 0.181 0.425 0.182 0.6694 Tdk Signifikan ❌
#> X5 0.832 0.671 1.536 0.2152 Tdk Signifikan ❌
#> X6 -0.065 0.727 0.008 0.9290 Tdk Signifikan ❌
#> X7 1.877 0.672 7.810 0.0052 Signifikan ✅
#> X8 0.315 0.752 0.176 0.6753 Tdk Signifikan ❌
#> X9 -0.806 0.828 0.947 0.3304 Tdk Signifikan ❌
#> X10 0.049 1.003 0.002 0.9614 Tdk Signifikan ❌
#> X11 1.723 0.948 3.305 0.0691 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 6.028
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.6441
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.008 [0.000, 0.145] Odds baseline
#> X1 1.156 [0.339, 3.944] ↑ lebih berisiko miskin
#> X2 0.104 [0.022, 0.498] ↓ lebih kecil risiko miskin
#> X3 1.226 [0.271, 5.544] ↑ lebih berisiko miskin
#> X4 1.199 [0.522, 2.755] ↑ lebih berisiko miskin
#> X5 2.297 [0.617, 8.558] ↑ lebih berisiko miskin
#> X6 0.937 [0.225, 3.897] ↓ lebih kecil risiko miskin
#> X7 6.537 [1.752, 24.390] ↑ lebih berisiko miskin
#> X8 1.371 [0.314, 5.988] ↑ lebih berisiko miskin
#> X9 0.447 [0.088, 2.263] ↓ lebih kecil risiko miskin
#> X10 1.050 [0.147, 7.497] ↑ lebih berisiko miskin
#> X11 5.599 [0.874, 35.863] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 11 — Sumbar + ART>1 | IND2 No-Duplikat
model_ind2_nodup_art2_smb <- analisis_reglog(data_ind2_nodup_art2_smb,
"Sumbar + ART>1 | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND2 No-Duplikat
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 705 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 694 obs ( 98.4 %)
#> Miskin (1): 11 obs ( 1.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 95 obs ( 13.48 %)
#> 1 : 226 obs ( 32.06 %)
#> 2 : 384 obs ( 54.47 %)
#>
#> ── X2 ──
#> 0 : 226 obs ( 32.06 %)
#> 1 : 479 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 235 obs ( 33.33 %)
#> 1 : 470 obs ( 66.67 %)
#>
#> ── X4 ──
#> 0 : 387 obs ( 54.89 %)
#> 1 : 70 obs ( 9.93 %)
#> 2 : 248 obs ( 35.18 %)
#>
#> ── X5 ──
#> 0 : 506 obs ( 71.77 %)
#> 1 : 199 obs ( 28.23 %)
#>
#> ── X6 ──
#> 0 : 338 obs ( 47.94 %)
#> 1 : 367 obs ( 52.06 %)
#>
#> ── X7 ──
#> 0 : 598 obs ( 84.82 %)
#> 1 : 107 obs ( 15.18 %)
#>
#> ── X8 ──
#> 0 : 585 obs ( 82.98 %)
#> 1 : 120 obs ( 17.02 %)
#>
#> ── X9 ──
#> 0 : 438 obs ( 62.13 %)
#> 1 : 267 obs ( 37.87 %)
#>
#> ── X10 ──
#> 0 : 69 obs ( 9.79 %)
#> 1 : 636 obs ( 90.21 %)
#>
#> ── X11 ──
#> 0 : 503 obs ( 71.35 %)
#> 1 : 202 obs ( 28.65 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 113.354
#> -2 Log Likelihood (Model Penuh) : 92.69
#> G² hitung : 20.664
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.037
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -4.971 1.442 11.886 0.0006 Signifikan ✅
#> X1 0.324 0.623 0.270 0.6036 Tdk Signifikan ❌
#> X2 -2.259 0.804 7.895 0.0050 Signifikan ✅
#> X3 0.238 0.755 0.100 0.7523 Tdk Signifikan ❌
#> X4 0.153 0.411 0.138 0.7108 Tdk Signifikan ❌
#> X5 0.762 0.670 1.296 0.2550 Tdk Signifikan ❌
#> X6 -0.092 0.689 0.018 0.8934 Tdk Signifikan ❌
#> X7 1.700 0.665 6.538 0.0106 Signifikan ✅
#> X8 0.236 0.740 0.102 0.7498 Tdk Signifikan ❌
#> X9 -0.766 0.824 0.863 0.3530 Tdk Signifikan ❌
#> X10 -0.029 0.993 0.001 0.9766 Tdk Signifikan ❌
#> X11 1.684 0.952 3.127 0.0770 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 10.845
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.2106
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.007 [0.000, 0.117] Odds baseline
#> X1 1.382 [0.407, 4.689] ↑ lebih berisiko miskin
#> X2 0.104 [0.022, 0.505] ↓ lebih kecil risiko miskin
#> X3 1.269 [0.289, 5.573] ↑ lebih berisiko miskin
#> X4 1.165 [0.520, 2.609] ↑ lebih berisiko miskin
#> X5 2.143 [0.577, 7.962] ↑ lebih berisiko miskin
#> X6 0.912 [0.236, 3.520] ↓ lebih kecil risiko miskin
#> X7 5.473 [1.487, 20.141] ↑ lebih berisiko miskin
#> X8 1.266 [0.297, 5.397] ↑ lebih berisiko miskin
#> X9 0.465 [0.092, 2.339] ↓ lebih kecil risiko miskin
#> X10 0.971 [0.139, 6.806] ↓ lebih kecil risiko miskin
#> X11 5.388 [0.833, 34.842] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 12 — Sumbar + ART>1 | IND2 Max-Age
model_ind2_maxage_art2_smb <- analisis_reglog(data_ind2_maxage_art2_smb,
"Sumbar + ART>1 | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND2 Max-Age
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 706 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 695 obs ( 98.4 %)
#> Miskin (1): 11 obs ( 1.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 95 obs ( 13.46 %)
#> 1 : 226 obs ( 32.01 %)
#> 2 : 385 obs ( 54.53 %)
#>
#> ── X2 ──
#> 0 : 226 obs ( 32.01 %)
#> 1 : 480 obs ( 67.99 %)
#>
#> ── X3 ──
#> 0 : 236 obs ( 33.43 %)
#> 1 : 470 obs ( 66.57 %)
#>
#> ── X4 ──
#> 0 : 387 obs ( 54.82 %)
#> 1 : 71 obs ( 10.06 %)
#> 2 : 248 obs ( 35.13 %)
#>
#> ── X5 ──
#> 0 : 507 obs ( 71.81 %)
#> 1 : 199 obs ( 28.19 %)
#>
#> ── X6 ──
#> 0 : 338 obs ( 47.88 %)
#> 1 : 368 obs ( 52.12 %)
#>
#> ── X7 ──
#> 0 : 599 obs ( 84.84 %)
#> 1 : 107 obs ( 15.16 %)
#>
#> ── X8 ──
#> 0 : 586 obs ( 83 %)
#> 1 : 120 obs ( 17 %)
#>
#> ── X9 ──
#> 0 : 438 obs ( 62.04 %)
#> 1 : 268 obs ( 37.96 %)
#>
#> ── X10 ──
#> 0 : 69 obs ( 9.77 %)
#> 1 : 637 obs ( 90.23 %)
#>
#> ── X11 ──
#> 0 : 504 obs ( 71.39 %)
#> 1 : 202 obs ( 28.61 %)
#>
#> ▶ 2. PEMODELAN REGRESI LOGISTIK BINER
#> ────────────────────────────────────────────────────────
#> Model: logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#>
#> ▶ 3. UJI SIMULTAN (Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 113.386
#> -2 Log Likelihood (Model Penuh) : 92.691
#> G² hitung : 20.694
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.0367
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Wald Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#>
#> Variabel B SE Wald p-value Kesimpulan
#> ----------------------------------------------------------------------
#> (Intercept) -4.971 1.442 11.886 0.0006 Signifikan ✅
#> X1 0.323 0.623 0.269 0.6038 Tdk Signifikan ❌
#> X2 -2.260 0.804 7.899 0.0049 Signifikan ✅
#> X3 0.239 0.755 0.100 0.7519 Tdk Signifikan ❌
#> X4 0.153 0.411 0.137 0.7108 Tdk Signifikan ❌
#> X5 0.762 0.670 1.296 0.2550 Tdk Signifikan ❌
#> X6 -0.092 0.689 0.018 0.8933 Tdk Signifikan ❌
#> X7 1.700 0.665 6.540 0.0105 Signifikan ✅
#> X8 0.236 0.740 0.102 0.7498 Tdk Signifikan ❌
#> X9 -0.766 0.824 0.864 0.3528 Tdk Signifikan ❌
#> X10 -0.029 0.993 0.001 0.9766 Tdk Signifikan ❌
#> X11 1.685 0.952 3.129 0.0769 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 10.849
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.2104
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B))
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.007 [0.000, 0.117] Odds baseline
#> X1 1.382 [0.407, 4.688] ↑ lebih berisiko miskin
#> X2 0.104 [0.022, 0.505] ↓ lebih kecil risiko miskin
#> X3 1.269 [0.289, 5.574] ↑ lebih berisiko miskin
#> X4 1.165 [0.520, 2.609] ↑ lebih berisiko miskin
#> X5 2.143 [0.577, 7.963] ↑ lebih berisiko miskin
#> X6 0.912 [0.236, 3.519] ↓ lebih kecil risiko miskin
#> X7 5.474 [1.487, 20.145] ↑ lebih berisiko miskin
#> X8 1.266 [0.297, 5.398] ↑ lebih berisiko miskin
#> X9 0.465 [0.092, 2.339] ↓ lebih kecil risiko miskin
#> X10 0.971 [0.139, 6.806] ↓ lebih kecil risiko miskin
#> X11 5.390 [0.834, 34.851] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════🔷 Bagian 7 — Model Utama: Firth Logistic Regression
Tujuan & Latar Belakang: Firth Logistic Regression (Firth, 1993) merupakan metode estimasi berbasis Penalized Maximum Likelihood (PML) yang menambahkan penalti Jeffreys sebagai koreksi bias ke dalam fungsi log-likelihood standar:
\[\ell^*(\beta) = \ell(\beta) + \frac{1}{2} \log |I(\beta)|\]
di mana \(I(\beta)\) adalah matriks informasi Fisher. Penalti ini menghasilkan estimator yang bias-reduced — secara asimptotik lebih dekat ke nilai populasi dibanding MLE biasa, khususnya ketika:
- Data imbalance: proporsi kelas minoritas (miskin) sangat kecil
- Complete/quasi-complete separation: prediktor dapat memisahkan Y=0 dan Y=1 secara sempurna
- Sampel kecil: terutama pada dataset level provinsi (Sumatera Barat)
Perbedaan Teknis Utama dibanding Model Baseline:
| Aspek | Logistik Biasa (glm) |
Firth (logistf) |
|---|---|---|
| Estimasi | MLE biasa | Penalized MLE |
| Uji Parsial | Wald Test (W²) | Profile Penalized LRT |
| Confidence Interval | Wald CI | Profile Penalized CI |
| Fitted values | fitted(model) |
model$predict |
| Log-likelihood | -model$deviance/2 |
model$loglik[2] |
| Cocok untuk | Data seimbang, n besar | Data imbalance, separation ✅ |
1. Fungsi Utama Firth Logistic Regression
Prosedur: Fungsi
analisis_firth()mengikuti struktur yang sama persis dengananalisis_reglog(), namun menggantiglm()denganlogistf(). Seluruh komponen uji juga disesuaikan dengan output objeklogistf.Teknis —
logistf(firth = TRUE, pl = TRUE): Argumenfirth = TRUEmengaktifkan penalti Jeffreys. Argumenpl = TRUE(profile likelihood) mengaktifkan penghitungan CI berbasis profile likelihood, yang lebih akurat dibanding CI Wald terutama ketika distribusi koefisien tidak simetris — kondisi yang umum pada data dengan separation.Teknis — Uji Simultan Penalized: G² dihitung dari log-likelihood model penuh dan null:
G² = 2 × (loglik_full − loglik_null), menggunakanmodel$loglik[2]yang menyimpan nilai log-likelihood setelah estimasi.Teknis — Uji Parsial Profile Penalized LRT: p-value parsial diambil langsung dari
model$prob, yang merupakan hasil profile penalized likelihood ratio test untuk setiap koefisien — bukan Wald, sehingga lebih reliable untuk kondisi separation dan imbalance.Teknis — CI Profile Penalized:
model$ci.lowerdanmodel$ci.uppermenyimpan batas CI skala log-odds yang dihitung via profile likelihood. Hasilexp()menghasilkan CI Odds Ratio yang lebih akurat.
# ============================================================
# FUNGSI UTAMA ANALISIS FIRTH LOGISTIC REGRESSION
# (Bias-Reduced Logistic Regression — Menangani Data Imbalance)
# ============================================================
library(dplyr)
library(logistf) # Firth Logistic Regression
library(ResourceSelection) # Hosmer-Lemeshow
analisis_firth <- function(data, nama_dataset, drop_var_na = FALSE) {
cat("\n")
cat("╔══════════════════════════════════════════════════════╗\n")
cat(" DATASET:", nama_dataset, "\n")
cat(" MODEL : Firth Logistic Regression (Bias Reduction)\n")
cat("╚══════════════════════════════════════════════════════╝\n\n")
xvars <- c("X1","X2","X3","X4","X5","X6","X7","X8","X9","X10","X11")
semua_vars <- c("status_miskin", xvars)
missing_vars <- semua_vars[!semua_vars %in% names(data)]
if (length(missing_vars) > 0) {
cat("⚠️ Variabel tidak ditemukan:", paste(missing_vars, collapse = ", "), "\n")
return(invisible(NULL))
}
df <- data %>%
select(all_of(semua_vars)) %>%
mutate(across(everything(), as.numeric))
# ── Drop variabel mengandung NA (opsional) ──────────────────────
if (drop_var_na) {
var_na <- xvars[sapply(xvars, function(v) any(is.na(df[[v]])))]
if (length(var_na) > 0)
cat("⚠️ Variabel mengandung NA (dikeluarkan):",
paste(var_na, collapse = ", "), "\n\n")
xvars_model <- setdiff(xvars, var_na)
} else {
xvars_model <- xvars
}
df <- df %>%
select(all_of(c("status_miskin", xvars_model))) %>%
na.omit()
cat("── Observasi valid:", nrow(df), "baris\n\n")
# ── Cek variabel konstan ────────────────────────────────────────
konstan <- xvars_model[sapply(xvars_model, function(v) length(unique(df[[v]])) < 2)]
if (length(konstan) > 0)
cat("⚠️ Variabel KONSTAN (dikeluarkan):",
paste(konstan, collapse = ", "), "\n\n")
xvars_model <- setdiff(xvars_model, konstan)
# ── 1. EDA ─────────────────────────────────────────────────────
cat("▶ 1. EDA\n")
cat("────────────────────────────────────────────────────────\n")
tbl_y <- table(df$status_miskin)
prop_y <- prop.table(tbl_y) * 100
cat(" [Y] Distribusi status_miskin:\n")
cat(" Tidak Miskin (0):", tbl_y["0"], "obs (", round(prop_y["0"], 1), "%)\n")
cat(" Miskin (1):", tbl_y["1"], "obs (", round(prop_y["1"], 1), "%)\n\n")
cat(" [X] Distribusi & Proporsi Variabel:\n")
for (v in xvars_model) {
cat("\n ──", v, "──\n")
tab <- table(df[[v]])
prop <- prop.table(tab) * 100
for (lvl in names(tab))
cat(" ", lvl, ":", tab[lvl], "obs (", round(prop[lvl], 2), "%)\n")
}
cat("\n")
# ── 2. PEMODELAN FIRTH ──────────────────────────────────────────
cat("▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION\n")
cat("────────────────────────────────────────────────────────\n")
formula_model <- as.formula(
paste("status_miskin ~", paste(xvars_model, collapse = " + "))
)
model <- logistf(formula_model, data = df, firth = TRUE, pl = TRUE)
model_null <- logistf(status_miskin ~ 1, data = df, firth = TRUE, pl = TRUE)
cat(" Model : logit[P(Y=1)] = β0 +",
paste(paste0("β·", xvars_model), collapse = " + "), "\n")
cat(" Metode: Penalized Maximum Likelihood (Firth, 1993)\n\n")
# ── 3. UJI SIMULTAN (Penalized LRT / G²) ───────────────────────
cat("▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)\n")
cat("────────────────────────────────────────────────────────\n")
loglik_null <- model_null$loglik[2]
loglik_full <- model$loglik[2]
G2 <- 2 * (loglik_full - loglik_null)
df_G2 <- length(xvars_model)
chi_kritis <- qchisq(0.95, df = df_G2)
pval_G2 <- pchisq(G2, df = df_G2, lower.tail = FALSE)
cat(" H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0\n")
cat(" -2 Log Likelihood (Model Null) :", round(-2 * loglik_null, 3), "\n")
cat(" -2 Log Likelihood (Model Penuh) :", round(-2 * loglik_full, 3), "\n")
cat(" G² hitung :", round(G2, 3), "\n")
cat(" df :", df_G2, "\n")
cat(" χ² tabel (α=5%, df=", df_G2, ") :", round(chi_kritis, 3), "\n")
cat(" p-value :", round(pval_G2, 4), "\n")
cat(ifelse(G2 > chi_kritis,
" ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan\n\n",
" ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan\n\n"))
# ── 4. UJI PARSIAL (Profile Penalized LRT) ──────────────────────
cat("▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)\n")
cat("────────────────────────────────────────────────────────\n")
koef <- model$coefficients
se_koef <- sqrt(diag(vcov(model)))
pval_koef <- model$prob
chi_kritis1 <- qchisq(0.95, df = 1)
cat(" H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel =",
round(chi_kritis1, 3), ")\n")
cat(" Statistik uji: Profile Penalized LRT (bukan Wald)\n\n")
cat(sprintf(" %-20s %8s %8s %10s %8s\n",
"Variabel", "B", "SE", "p-value", "Kesimpulan"))
cat(" ", paste(rep("-", 65), collapse = ""), "\n")
for (nm in names(koef)) {
ket <- ifelse(pval_koef[nm] < 0.05, "Signifikan ✅", "Tdk Signifikan ❌")
cat(sprintf(" %-20s %8.3f %8.3f %10.4f %s\n",
nm, koef[nm], se_koef[nm], pval_koef[nm], ket))
}
cat("\n")
# ── 5. GOODNESS OF FIT (Hosmer-Lemeshow) ────────────────────────
cat("▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)\n")
cat("────────────────────────────────────────────────────────\n")
fitted_prob <- model$predict
hl_test <- tryCatch(
hoslem.test(df$status_miskin, fitted_prob, g = 10),
error = function(e) NULL
)
if (!is.null(hl_test)) {
chi_kritis_hl <- qchisq(0.95, df = 8)
cat(" H0: Model fit dengan data\n")
cat(" χ² Hosmer-Lemeshow :", round(hl_test$statistic, 3), "\n")
cat(" df :", hl_test$parameter, "\n")
cat(" χ² tabel (df=8) :", round(chi_kritis_hl, 3), "\n")
cat(" p-value :", round(hl_test$p.value, 4), "\n")
cat(ifelse(hl_test$p.value > 0.05,
" ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data\n\n",
" ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data\n\n"))
} else {
cat(" ⚠️ Hosmer-Lemeshow tidak dapat dihitung\n\n")
}
# ── 6. INTERPRETASI PARAMETER (Odds Ratio — Profile CI) ─────────
cat("▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)\n")
cat("────────────────────────────────────────────────────────\n")
or <- exp(koef)
ci_low <- exp(model$ci.lower)
ci_high <- exp(model$ci.upper)
cat(sprintf(" %-20s %8s %10s %s\n",
"Variabel", "exp(B)", "95% CI", "Interpretasi Singkat"))
cat(" ", paste(rep("-", 75), collapse = ""), "\n")
for (nm in names(or)) {
if (nm == "(Intercept)") {
cat(sprintf(" %-20s %8.3f [%6.3f, %6.3f] Odds baseline\n",
nm, or[nm], ci_low[nm], ci_high[nm]))
} else {
arah <- ifelse(or[nm] > 1,
"↑ lebih berisiko miskin",
"↓ lebih kecil risiko miskin")
cat(sprintf(" %-20s %8.3f [%6.3f, %6.3f] %s\n",
nm, or[nm], ci_low[nm], ci_high[nm], arah))
}
}
cat("\n")
cat("══════════════════════════════════════════════════════\n\n")
return(invisible(model))
}🔷 Bagian 8 — Hasil: Firth Logistic Regression
Hasil analisis Firth Logistic Regression untuk seluruh 12 dataset disajikan di bawah. Bandingkan hasil uji parsial dan nilai Odds Ratio di sini dengan hasil Bagian 6 untuk mengidentifikasi apakah koreksi bias mengubah signifikansi atau arah efek variabel prediktor.
#> # A tibble: 12 × 6
#> Dataset Filter_Cerai Duplikat Filter_ART Filter_Provinsi N
#> <chr> <chr> <chr> <chr> <chr> <int>
#> 1 data_ind1_filtered Cerai + KRT Awal — — 27537
#> 2 data_ind2_nodup Cerai Dihapus — — 28828
#> 3 data_ind2_maxage Cerai Max Umur — — 28860
#> 4 data_ind1_filtered_ar… Cerai + KRT Awal >1 — 16035
#> 5 data_ind2_nodup_art2 Cerai Dihapus >1 — 17326
#> 6 data_ind2_maxage_art2 Cerai Max Umur >1 — 17358
#> 7 data_ind1_filtered_smb Cerai + KRT Awal — Sumbar 1135
#> 8 data_ind2_nodup_smb Cerai Dihapus — Sumbar 1195
#> 9 data_ind2_maxage_smb Cerai Max Umur — Sumbar 1196
#> 10 data_ind1_filtered_ar… Cerai + KRT Awal >1 Sumbar 677
#> 11 data_ind2_nodup_art2_… Cerai Dihapus >1 Sumbar 737
#> 12 data_ind2_maxage_art2… Cerai Max Umur >1 Sumbar 738Dataset 1 — Nasional | IND1 Filtered
Dataset terbesar dengan cakupan nasional dan definisi IND1 (KRT). Merupakan benchmark utama analisis.
drop_var_na = FALSEdigunakan di sini — variabel dengan NA tidak otomatis dibuang, melainkan baris yang mengandung NA pada variabel yang dipakai yang dihapus.
model_ind1_filtered <- analisis_firth(data_ind1_filtered,
"Nasional | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND1 Filtered
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 26392 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 25495 obs ( 96.6 %)
#> Miskin (1): 897 obs ( 3.4 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2113 obs ( 8.01 %)
#> 1 : 4574 obs ( 17.33 %)
#> 2 : 19705 obs ( 74.66 %)
#>
#> ── X2 ──
#> 0 : 5284 obs ( 20.02 %)
#> 1 : 21108 obs ( 79.98 %)
#>
#> ── X3 ──
#> 0 : 7781 obs ( 29.48 %)
#> 1 : 18611 obs ( 70.52 %)
#>
#> ── X4 ──
#> 0 : 12721 obs ( 48.2 %)
#> 1 : 2642 obs ( 10.01 %)
#> 2 : 11029 obs ( 41.79 %)
#>
#> ── X5 ──
#> 0 : 16003 obs ( 60.64 %)
#> 1 : 10389 obs ( 39.36 %)
#>
#> ── X6 ──
#> 0 : 11996 obs ( 45.45 %)
#> 1 : 14396 obs ( 54.55 %)
#>
#> ── X7 ──
#> 0 : 25131 obs ( 95.22 %)
#> 1 : 1261 obs ( 4.78 %)
#>
#> ── X8 ──
#> 0 : 21668 obs ( 82.1 %)
#> 1 : 4724 obs ( 17.9 %)
#>
#> ── X9 ──
#> 0 : 15156 obs ( 57.43 %)
#> 1 : 11236 obs ( 42.57 %)
#>
#> ── X10 ──
#> 0 : 7511 obs ( 28.46 %)
#> 1 : 18881 obs ( 71.54 %)
#>
#> ── X11 ──
#> 0 : 13764 obs ( 52.15 %)
#> 1 : 12628 obs ( 47.85 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 7823.275
#> -2 Log Likelihood (Model Penuh) : 7767.17
#> G² hitung : 56.105
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.917 0.186 0.0000 Signifikan ✅
#> X1 0.473 0.088 0.0000 Signifikan ✅
#> X2 -1.240 0.084 0.0000 Signifikan ✅
#> X3 0.021 0.088 0.8100 Tdk Signifikan ❌
#> X4 0.088 0.046 0.0559 Tdk Signifikan ❌
#> X5 0.178 0.075 0.0177 Signifikan ✅
#> X6 0.252 0.086 0.0034 Signifikan ✅
#> X7 2.314 0.086 0.0000 Signifikan ✅
#> X8 0.104 0.088 0.2449 Tdk Signifikan ❌
#> X9 -0.383 0.077 0.0000 Signifikan ✅
#> X10 -0.381 0.095 0.0000 Signifikan ✅
#> X11 0.648 0.109 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 11.653
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.1674
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.020 [ 0.014, 0.029] Odds baseline
#> X1 1.605 [ 1.353, 1.916] ↑ lebih berisiko miskin
#> X2 0.289 [ 0.245, 0.341] ↓ lebih kecil risiko miskin
#> X3 1.021 [ 0.860, 1.217] ↑ lebih berisiko miskin
#> X4 1.092 [ 0.998, 1.197] ↑ lebih berisiko miskin
#> X5 1.194 [ 1.031, 1.383] ↑ lebih berisiko miskin
#> X6 1.286 [ 1.086, 1.526] ↑ lebih berisiko miskin
#> X7 10.112 [ 8.541, 11.951] ↑ lebih berisiko miskin
#> X8 1.110 [ 0.930, 1.317] ↑ lebih berisiko miskin
#> X9 0.682 [ 0.586, 0.793] ↓ lebih kecil risiko miskin
#> X10 0.683 [ 0.566, 0.822] ↓ lebih kecil risiko miskin
#> X11 1.912 [ 1.544, 2.367] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 2 — Nasional | IND2 No-Duplikat
model_ind2_nodup <- analisis_firth(data_ind2_nodup,
"Nasional | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND2 No-Duplikat
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 27617 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 26613 obs ( 96.4 %)
#> Miskin (1): 1004 obs ( 3.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2358 obs ( 8.54 %)
#> 1 : 4933 obs ( 17.86 %)
#> 2 : 20326 obs ( 73.6 %)
#>
#> ── X2 ──
#> 0 : 5971 obs ( 21.62 %)
#> 1 : 21646 obs ( 78.38 %)
#>
#> ── X3 ──
#> 0 : 8437 obs ( 30.55 %)
#> 1 : 19180 obs ( 69.45 %)
#>
#> ── X4 ──
#> 0 : 13547 obs ( 49.05 %)
#> 1 : 2817 obs ( 10.2 %)
#> 2 : 11253 obs ( 40.75 %)
#>
#> ── X5 ──
#> 0 : 16911 obs ( 61.23 %)
#> 1 : 10706 obs ( 38.77 %)
#>
#> ── X6 ──
#> 0 : 12724 obs ( 46.07 %)
#> 1 : 14893 obs ( 53.93 %)
#>
#> ── X7 ──
#> 0 : 25988 obs ( 94.1 %)
#> 1 : 1629 obs ( 5.9 %)
#>
#> ── X8 ──
#> 0 : 22671 obs ( 82.09 %)
#> 1 : 4946 obs ( 17.91 %)
#>
#> ── X9 ──
#> 0 : 16056 obs ( 58.14 %)
#> 1 : 11561 obs ( 41.86 %)
#>
#> ── X10 ──
#> 0 : 7686 obs ( 27.83 %)
#> 1 : 19931 obs ( 72.17 %)
#>
#> ── X11 ──
#> 0 : 14711 obs ( 53.27 %)
#> 1 : 12906 obs ( 46.73 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 8619.569
#> -2 Log Likelihood (Model Penuh) : 8562.028
#> G² hitung : 57.541
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.919 0.169 0.0000 Signifikan ✅
#> X1 0.487 0.079 0.0000 Signifikan ✅
#> X2 -1.233 0.080 0.0000 Signifikan ✅
#> X3 0.042 0.082 0.6102 Tdk Signifikan ❌
#> X4 0.068 0.043 0.1191 Tdk Signifikan ❌
#> X5 0.204 0.071 0.0042 Signifikan ✅
#> X6 0.261 0.080 0.0011 Signifikan ✅
#> X7 2.106 0.080 0.0000 Signifikan ✅
#> X8 0.062 0.084 0.4612 Tdk Signifikan ❌
#> X9 -0.369 0.073 0.0000 Signifikan ✅
#> X10 -0.343 0.091 0.0001 Signifikan ✅
#> X11 0.674 0.103 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 19.971
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.0104
#> ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.020 [ 0.014, 0.028] Odds baseline
#> X1 1.627 [ 1.396, 1.906] ↑ lebih berisiko miskin
#> X2 0.291 [ 0.249, 0.341] ↓ lebih kecil risiko miskin
#> X3 1.043 [ 0.888, 1.228] ↑ lebih berisiko miskin
#> X4 1.070 [ 0.983, 1.166] ↑ lebih berisiko miskin
#> X5 1.226 [ 1.066, 1.408] ↑ lebih berisiko miskin
#> X6 1.298 [ 1.109, 1.520] ↑ lebih berisiko miskin
#> X7 8.212 [ 7.016, 9.597] ↑ lebih berisiko miskin
#> X8 1.064 [ 0.901, 1.251] ↑ lebih berisiko miskin
#> X9 0.692 [ 0.599, 0.797] ↓ lebih kecil risiko miskin
#> X10 0.710 [ 0.593, 0.847] ↓ lebih kecil risiko miskin
#> X11 1.963 [ 1.602, 2.403] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 3 — Nasional | IND2 Max-Age
model_ind2_maxage <- analisis_firth(data_ind2_maxage,
"Nasional | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional | IND2 Max-Age
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 27642 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 26633 obs ( 96.3 %)
#> Miskin (1): 1009 obs ( 3.7 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 2359 obs ( 8.53 %)
#> 1 : 4937 obs ( 17.86 %)
#> 2 : 20346 obs ( 73.61 %)
#>
#> ── X2 ──
#> 0 : 5980 obs ( 21.63 %)
#> 1 : 21662 obs ( 78.37 %)
#>
#> ── X3 ──
#> 0 : 8450 obs ( 30.57 %)
#> 1 : 19192 obs ( 69.43 %)
#>
#> ── X4 ──
#> 0 : 13559 obs ( 49.05 %)
#> 1 : 2824 obs ( 10.22 %)
#> 2 : 11259 obs ( 40.73 %)
#>
#> ── X5 ──
#> 0 : 16926 obs ( 61.23 %)
#> 1 : 10716 obs ( 38.77 %)
#>
#> ── X6 ──
#> 0 : 12736 obs ( 46.07 %)
#> 1 : 14906 obs ( 53.93 %)
#>
#> ── X7 ──
#> 0 : 26002 obs ( 94.07 %)
#> 1 : 1640 obs ( 5.93 %)
#>
#> ── X8 ──
#> 0 : 22691 obs ( 82.09 %)
#> 1 : 4951 obs ( 17.91 %)
#>
#> ── X9 ──
#> 0 : 16073 obs ( 58.15 %)
#> 1 : 11569 obs ( 41.85 %)
#>
#> ── X10 ──
#> 0 : 7694 obs ( 27.83 %)
#> 1 : 19948 obs ( 72.17 %)
#>
#> ── X11 ──
#> 0 : 14728 obs ( 53.28 %)
#> 1 : 12914 obs ( 46.72 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 8654.173
#> -2 Log Likelihood (Model Penuh) : 8596.573
#> G² hitung : 57.6
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.915 0.169 0.0000 Signifikan ✅
#> X1 0.493 0.079 0.0000 Signifikan ✅
#> X2 -1.232 0.080 0.0000 Signifikan ✅
#> X3 0.037 0.082 0.6524 Tdk Signifikan ❌
#> X4 0.071 0.043 0.1022 Tdk Signifikan ❌
#> X5 0.197 0.071 0.0056 Signifikan ✅
#> X6 0.246 0.080 0.0020 Signifikan ✅
#> X7 2.108 0.079 0.0000 Signifikan ✅
#> X8 0.067 0.083 0.4227 Tdk Signifikan ❌
#> X9 -0.366 0.073 0.0000 Signifikan ✅
#> X10 -0.349 0.091 0.0001 Signifikan ✅
#> X11 0.676 0.103 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 19.888
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.0108
#> ❌ KEPUTUSAN: Tolak H0 → Model TIDAK fit dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.020 [ 0.014, 0.028] Odds baseline
#> X1 1.637 [ 1.404, 1.917] ↑ lebih berisiko miskin
#> X2 0.292 [ 0.249, 0.341] ↓ lebih kecil risiko miskin
#> X3 1.038 [ 0.884, 1.221] ↑ lebih berisiko miskin
#> X4 1.073 [ 0.986, 1.169] ↑ lebih berisiko miskin
#> X5 1.217 [ 1.059, 1.398] ↑ lebih berisiko miskin
#> X6 1.279 [ 1.093, 1.497] ↑ lebih berisiko miskin
#> X7 8.231 [ 7.036, 9.614] ↑ lebih berisiko miskin
#> X8 1.070 [ 0.906, 1.257] ↑ lebih berisiko miskin
#> X9 0.694 [ 0.601, 0.800] ↓ lebih kecil risiko miskin
#> X10 0.706 [ 0.590, 0.842] ↓ lebih kecil risiko miskin
#> X11 1.966 [ 1.605, 2.407] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 4 — Nasional + ART>1 | IND1 Filtered
model_ind1_filtered_art2 <- analisis_firth(data_ind1_filtered_art2,
"Nasional + ART>1 | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND1 Filtered
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 15421 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 14558 obs ( 94.4 %)
#> Miskin (1): 863 obs ( 5.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1595 obs ( 10.34 %)
#> 1 : 3457 obs ( 22.42 %)
#> 2 : 10369 obs ( 67.24 %)
#>
#> ── X2 ──
#> 0 : 4649 obs ( 30.15 %)
#> 1 : 10772 obs ( 69.85 %)
#>
#> ── X3 ──
#> 0 : 5194 obs ( 33.68 %)
#> 1 : 10227 obs ( 66.32 %)
#>
#> ── X4 ──
#> 0 : 8185 obs ( 53.08 %)
#> 1 : 1597 obs ( 10.36 %)
#> 2 : 5639 obs ( 36.57 %)
#>
#> ── X5 ──
#> 0 : 10418 obs ( 67.56 %)
#> 1 : 5003 obs ( 32.44 %)
#>
#> ── X6 ──
#> 0 : 7434 obs ( 48.21 %)
#> 1 : 7987 obs ( 51.79 %)
#>
#> ── X7 ──
#> 0 : 14160 obs ( 91.82 %)
#> 1 : 1261 obs ( 8.18 %)
#>
#> ── X8 ──
#> 0 : 12770 obs ( 82.81 %)
#> 1 : 2651 obs ( 17.19 %)
#>
#> ── X9 ──
#> 0 : 9716 obs ( 63 %)
#> 1 : 5705 obs ( 37 %)
#>
#> ── X10 ──
#> 0 : 3107 obs ( 20.15 %)
#> 1 : 12314 obs ( 79.85 %)
#>
#> ── X11 ──
#> 0 : 9806 obs ( 63.59 %)
#> 1 : 5615 obs ( 36.41 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 6646.254
#> -2 Log Likelihood (Model Penuh) : 6590.296
#> G² hitung : 55.958
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.688 0.187 0.0000 Signifikan ✅
#> X1 0.494 0.088 0.0000 Signifikan ✅
#> X2 -0.894 0.083 0.0000 Signifikan ✅
#> X3 -0.008 0.092 0.9271 Tdk Signifikan ❌
#> X4 0.099 0.048 0.0401 Signifikan ✅
#> X5 0.338 0.077 0.0000 Signifikan ✅
#> X6 0.246 0.091 0.0065 Signifikan ✅
#> X7 1.805 0.085 0.0000 Signifikan ✅
#> X8 0.101 0.093 0.2799 Tdk Signifikan ❌
#> X9 -0.263 0.080 0.0010 Signifikan ✅
#> X10 -0.499 0.099 0.0000 Signifikan ✅
#> X11 0.757 0.109 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 9.658
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.2898
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [ 0.017, 0.036] Odds baseline
#> X1 1.638 [ 1.383, 1.953] ↑ lebih berisiko miskin
#> X2 0.409 [ 0.348, 0.482] ↓ lebih kecil risiko miskin
#> X3 0.992 [ 0.829, 1.190] ↓ lebih kecil risiko miskin
#> X4 1.104 [ 1.004, 1.214] ↑ lebih berisiko miskin
#> X5 1.402 [ 1.205, 1.629] ↑ lebih berisiko miskin
#> X6 1.279 [ 1.071, 1.529] ↑ lebih berisiko miskin
#> X7 6.079 [ 5.136, 7.182] ↑ lebih berisiko miskin
#> X8 1.107 [ 0.920, 1.325] ↑ lebih berisiko miskin
#> X9 0.769 [ 0.656, 0.899] ↓ lebih kecil risiko miskin
#> X10 0.607 [ 0.499, 0.737] ↓ lebih kecil risiko miskin
#> X11 2.132 [ 1.719, 2.641] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 5 — Nasional + ART>1 | IND2 No-Duplikat
model_ind2_nodup_art2 <- analisis_firth(data_ind2_nodup_art2,
"Nasional + ART>1 | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND2 No-Duplikat
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 16646 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 15676 obs ( 94.2 %)
#> Miskin (1): 970 obs ( 5.8 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1840 obs ( 11.05 %)
#> 1 : 3816 obs ( 22.92 %)
#> 2 : 10990 obs ( 66.02 %)
#>
#> ── X2 ──
#> 0 : 5336 obs ( 32.06 %)
#> 1 : 11310 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 5850 obs ( 35.14 %)
#> 1 : 10796 obs ( 64.86 %)
#>
#> ── X4 ──
#> 0 : 9011 obs ( 54.13 %)
#> 1 : 1772 obs ( 10.65 %)
#> 2 : 5863 obs ( 35.22 %)
#>
#> ── X5 ──
#> 0 : 11326 obs ( 68.04 %)
#> 1 : 5320 obs ( 31.96 %)
#>
#> ── X6 ──
#> 0 : 8162 obs ( 49.03 %)
#> 1 : 8484 obs ( 50.97 %)
#>
#> ── X7 ──
#> 0 : 15017 obs ( 90.21 %)
#> 1 : 1629 obs ( 9.79 %)
#>
#> ── X8 ──
#> 0 : 13773 obs ( 82.74 %)
#> 1 : 2873 obs ( 17.26 %)
#>
#> ── X9 ──
#> 0 : 10616 obs ( 63.78 %)
#> 1 : 6030 obs ( 36.22 %)
#>
#> ── X10 ──
#> 0 : 3282 obs ( 19.72 %)
#> 1 : 13364 obs ( 80.28 %)
#>
#> ── X11 ──
#> 0 : 10753 obs ( 64.6 %)
#> 1 : 5893 obs ( 35.4 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 7390.228
#> -2 Log Likelihood (Model Penuh) : 7332.86
#> G² hitung : 57.368
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.677 0.170 0.0000 Signifikan ✅
#> X1 0.498 0.079 0.0000 Signifikan ✅
#> X2 -0.899 0.079 0.0000 Signifikan ✅
#> X3 0.001 0.085 0.9863 Tdk Signifikan ❌
#> X4 0.077 0.045 0.0880 Tdk Signifikan ❌
#> X5 0.353 0.073 0.0000 Signifikan ✅
#> X6 0.269 0.084 0.0013 Signifikan ✅
#> X7 1.641 0.079 0.0000 Signifikan ✅
#> X8 0.061 0.087 0.4896 Tdk Signifikan ❌
#> X9 -0.246 0.076 0.0010 Signifikan ✅
#> X10 -0.456 0.095 0.0000 Signifikan ✅
#> X11 0.781 0.103 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 7.737
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.4596
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [ 0.018, 0.035] Odds baseline
#> X1 1.646 [ 1.414, 1.926] ↑ lebih berisiko miskin
#> X2 0.407 [ 0.349, 0.475] ↓ lebih kecil risiko miskin
#> X3 1.001 [ 0.848, 1.186] ↑ lebih berisiko miskin
#> X4 1.080 [ 0.989, 1.180] ↑ lebih berisiko miskin
#> X5 1.424 [ 1.235, 1.641] ↑ lebih berisiko miskin
#> X6 1.309 [ 1.111, 1.544] ↑ lebih berisiko miskin
#> X7 5.162 [ 4.417, 6.025] ↑ lebih berisiko miskin
#> X8 1.063 [ 0.893, 1.259] ↑ lebih berisiko miskin
#> X9 0.782 [ 0.673, 0.906] ↓ lebih kecil risiko miskin
#> X10 0.634 [ 0.525, 0.762] ↓ lebih kecil risiko miskin
#> X11 2.183 [ 1.781, 2.673] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 6 — Nasional + ART>1 | IND2 Max-Age
model_ind2_maxage_art2 <- analisis_firth(data_ind2_maxage_art2,
"Nasional + ART>1 | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Nasional + ART>1 | IND2 Max-Age
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 16671 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 15696 obs ( 94.2 %)
#> Miskin (1): 975 obs ( 5.8 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 1841 obs ( 11.04 %)
#> 1 : 3820 obs ( 22.91 %)
#> 2 : 11010 obs ( 66.04 %)
#>
#> ── X2 ──
#> 0 : 5345 obs ( 32.06 %)
#> 1 : 11326 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 5863 obs ( 35.17 %)
#> 1 : 10808 obs ( 64.83 %)
#>
#> ── X4 ──
#> 0 : 9023 obs ( 54.12 %)
#> 1 : 1779 obs ( 10.67 %)
#> 2 : 5869 obs ( 35.2 %)
#>
#> ── X5 ──
#> 0 : 11341 obs ( 68.03 %)
#> 1 : 5330 obs ( 31.97 %)
#>
#> ── X6 ──
#> 0 : 8174 obs ( 49.03 %)
#> 1 : 8497 obs ( 50.97 %)
#>
#> ── X7 ──
#> 0 : 15031 obs ( 90.16 %)
#> 1 : 1640 obs ( 9.84 %)
#>
#> ── X8 ──
#> 0 : 13793 obs ( 82.74 %)
#> 1 : 2878 obs ( 17.26 %)
#>
#> ── X9 ──
#> 0 : 10633 obs ( 63.78 %)
#> 1 : 6038 obs ( 36.22 %)
#>
#> ── X10 ──
#> 0 : 3290 obs ( 19.73 %)
#> 1 : 13381 obs ( 80.27 %)
#>
#> ── X11 ──
#> 0 : 10770 obs ( 64.6 %)
#> 1 : 5901 obs ( 35.4 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 7421.037
#> -2 Log Likelihood (Model Penuh) : 7363.61
#> G² hitung : 57.428
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0
#> ✅ KEPUTUSAN: Tolak H0 → Model signifikan secara simultan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -3.672 0.170 0.0000 Signifikan ✅
#> X1 0.504 0.079 0.0000 Signifikan ✅
#> X2 -0.899 0.079 0.0000 Signifikan ✅
#> X3 -0.003 0.085 0.9674 Tdk Signifikan ❌
#> X4 0.080 0.045 0.0734 Tdk Signifikan ❌
#> X5 0.347 0.072 0.0000 Signifikan ✅
#> X6 0.254 0.083 0.0023 Signifikan ✅
#> X7 1.644 0.079 0.0000 Signifikan ✅
#> X8 0.067 0.087 0.4455 Tdk Signifikan ❌
#> X9 -0.242 0.076 0.0012 Signifikan ✅
#> X10 -0.462 0.095 0.0000 Signifikan ✅
#> X11 0.782 0.103 0.0000 Signifikan ✅
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 6.649
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.575
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.025 [ 0.018, 0.035] Odds baseline
#> X1 1.655 [ 1.422, 1.937] ↑ lebih berisiko miskin
#> X2 0.407 [ 0.349, 0.475] ↓ lebih kecil risiko miskin
#> X3 0.997 [ 0.844, 1.179] ↓ lebih kecil risiko miskin
#> X4 1.084 [ 0.992, 1.184] ↑ lebih berisiko miskin
#> X5 1.415 [ 1.227, 1.630] ↑ lebih berisiko miskin
#> X6 1.289 [ 1.095, 1.519] ↑ lebih berisiko miskin
#> X7 5.175 [ 4.431, 6.036] ↑ lebih berisiko miskin
#> X8 1.069 [ 0.899, 1.265] ↑ lebih berisiko miskin
#> X9 0.785 [ 0.676, 0.909] ↓ lebih kecil risiko miskin
#> X10 0.630 [ 0.522, 0.757] ↓ lebih kecil risiko miskin
#> X11 2.186 [ 1.784, 2.676] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 7 — Sumatera Barat | IND1 Filtered
Pada dataset level provinsi, ukuran sampel jauh lebih kecil sehingga risiko quasi-complete separation meningkat. Keunggulan Firth dalam kondisi ini paling terasa.
model_ind1_filtered_smb <- analisis_firth(data_ind1_filtered_smb,
"Sumatera Barat | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND1 Filtered
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1069 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1057 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 101 obs ( 9.45 %)
#> 1 : 266 obs ( 24.88 %)
#> 2 : 702 obs ( 65.67 %)
#>
#> ── X2 ──
#> 0 : 210 obs ( 19.64 %)
#> 1 : 859 obs ( 80.36 %)
#>
#> ── X3 ──
#> 0 : 287 obs ( 26.85 %)
#> 1 : 782 obs ( 73.15 %)
#>
#> ── X4 ──
#> 0 : 487 obs ( 45.56 %)
#> 1 : 93 obs ( 8.7 %)
#> 2 : 489 obs ( 45.74 %)
#>
#> ── X5 ──
#> 0 : 662 obs ( 61.93 %)
#> 1 : 407 obs ( 38.07 %)
#>
#> ── X6 ──
#> 0 : 450 obs ( 42.1 %)
#> 1 : 619 obs ( 57.9 %)
#>
#> ── X7 ──
#> 0 : 979 obs ( 91.58 %)
#> 1 : 90 obs ( 8.42 %)
#>
#> ── X8 ──
#> 0 : 864 obs ( 80.82 %)
#> 1 : 205 obs ( 19.18 %)
#>
#> ── X9 ──
#> 0 : 613 obs ( 57.34 %)
#> 1 : 456 obs ( 42.66 %)
#>
#> ── X10 ──
#> 0 : 176 obs ( 16.46 %)
#> 1 : 893 obs ( 83.54 %)
#>
#> ── X11 ──
#> 0 : 608 obs ( 56.88 %)
#> 1 : 461 obs ( 43.12 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 129.121
#> -2 Log Likelihood (Model Penuh) : 117.344
#> G² hitung : 11.777
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.3807
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.566 1.118 0.0001 Signifikan ✅
#> X1 0.059 0.509 0.9194 Tdk Signifikan ❌
#> X2 -2.455 0.646 0.0010 Signifikan ✅
#> X3 -0.034 0.606 0.9609 Tdk Signifikan ❌
#> X4 0.182 0.335 0.6393 Tdk Signifikan ❌
#> X5 0.784 0.535 0.1999 Tdk Signifikan ❌
#> X6 0.039 0.579 0.9534 Tdk Signifikan ❌
#> X7 2.076 0.562 0.0023 Signifikan ✅
#> X8 0.522 0.550 0.4142 Tdk Signifikan ❌
#> X9 -0.898 0.631 0.1881 Tdk Signifikan ❌
#> X10 0.275 0.738 0.7345 Tdk Signifikan ❌
#> X11 1.724 0.776 0.0530 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 5.82
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.6674
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.010 [ 0.001, 0.101] Odds baseline
#> X1 1.061 [ 0.344, 3.651] ↑ lebih berisiko miskin
#> X2 0.086 [ 0.019, 0.373] ↓ lebih kecil risiko miskin
#> X3 0.966 [ 0.262, 4.367] ↓ lebih kecil risiko miskin
#> X4 1.199 [ 0.572, 2.705] ↑ lebih berisiko miskin
#> X5 2.190 [ 0.654, 7.448] ↑ lebih berisiko miskin
#> X6 1.040 [ 0.283, 4.127] ↑ lebih berisiko miskin
#> X7 7.972 [ 2.209, 27.461] ↑ lebih berisiko miskin
#> X8 1.685 [ 0.453, 5.528] ↑ lebih berisiko miskin
#> X9 0.407 [ 0.076, 1.505] ↓ lebih kecil risiko miskin
#> X10 1.317 [ 0.276, 7.974] ↑ lebih berisiko miskin
#> X11 5.605 [ 0.978, 35.042] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 8 — Sumatera Barat | IND2 No-Duplikat
model_ind2_nodup_smb <- analisis_firth(data_ind2_nodup_smb,
"Sumatera Barat | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND2 No-Duplikat
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1124 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1112 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 116 obs ( 10.32 %)
#> 1 : 283 obs ( 25.18 %)
#> 2 : 725 obs ( 64.5 %)
#>
#> ── X2 ──
#> 0 : 235 obs ( 20.91 %)
#> 1 : 889 obs ( 79.09 %)
#>
#> ── X3 ──
#> 0 : 310 obs ( 27.58 %)
#> 1 : 814 obs ( 72.42 %)
#>
#> ── X4 ──
#> 0 : 525 obs ( 46.71 %)
#> 1 : 98 obs ( 8.72 %)
#> 2 : 501 obs ( 44.57 %)
#>
#> ── X5 ──
#> 0 : 699 obs ( 62.19 %)
#> 1 : 425 obs ( 37.81 %)
#>
#> ── X6 ──
#> 0 : 483 obs ( 42.97 %)
#> 1 : 641 obs ( 57.03 %)
#>
#> ── X7 ──
#> 0 : 1017 obs ( 90.48 %)
#> 1 : 107 obs ( 9.52 %)
#>
#> ── X8 ──
#> 0 : 906 obs ( 80.6 %)
#> 1 : 218 obs ( 19.4 %)
#>
#> ── X9 ──
#> 0 : 651 obs ( 57.92 %)
#> 1 : 473 obs ( 42.08 %)
#>
#> ── X10 ──
#> 0 : 180 obs ( 16.01 %)
#> 1 : 944 obs ( 83.99 %)
#>
#> ── X11 ──
#> 0 : 653 obs ( 58.1 %)
#> 1 : 471 obs ( 41.9 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 130.331
#> -2 Log Likelihood (Model Penuh) : 118.415
#> G² hitung : 11.917
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.3699
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.817 1.122 0.0000 Signifikan ✅
#> X1 0.252 0.508 0.6617 Tdk Signifikan ❌
#> X2 -2.444 0.653 0.0013 Signifikan ✅
#> X3 0.024 0.598 0.9723 Tdk Signifikan ❌
#> X4 0.150 0.325 0.6880 Tdk Signifikan ❌
#> X5 0.718 0.533 0.2397 Tdk Signifikan ❌
#> X6 -0.012 0.550 0.9851 Tdk Signifikan ❌
#> X7 1.881 0.558 0.0050 Signifikan ✅
#> X8 0.462 0.543 0.4622 Tdk Signifikan ❌
#> X9 -0.871 0.632 0.2028 Tdk Signifikan ❌
#> X10 0.230 0.735 0.7760 Tdk Signifikan ❌
#> X11 1.709 0.779 0.0568 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 5.999
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.6473
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.008 [ 0.001, 0.079] Odds baseline
#> X1 1.287 [ 0.422, 4.424] ↑ lebih berisiko miskin
#> X2 0.087 [ 0.019, 0.387] ↓ lebih kecil risiko miskin
#> X3 1.024 [ 0.283, 4.517] ↑ lebih berisiko miskin
#> X4 1.162 [ 0.568, 2.542] ↑ lebih berisiko miskin
#> X5 2.051 [ 0.613, 6.958] ↑ lebih berisiko miskin
#> X6 0.988 [ 0.288, 3.592] ↓ lebih kecil risiko miskin
#> X7 6.562 [ 1.828, 22.620] ↑ lebih berisiko miskin
#> X8 1.588 [ 0.433, 5.105] ↑ lebih berisiko miskin
#> X9 0.418 [ 0.078, 1.548] ↓ lebih kecil risiko miskin
#> X10 1.259 [ 0.266, 7.593] ↑ lebih berisiko miskin
#> X11 5.523 [ 0.952, 34.660] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 9 — Sumatera Barat | IND2 Max-Age
model_ind2_maxage_smb <- analisis_firth(data_ind2_maxage_smb,
"Sumatera Barat | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumatera Barat | IND2 Max-Age
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 1125 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 1113 obs ( 98.9 %)
#> Miskin (1): 12 obs ( 1.1 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 116 obs ( 10.31 %)
#> 1 : 283 obs ( 25.16 %)
#> 2 : 726 obs ( 64.53 %)
#>
#> ── X2 ──
#> 0 : 235 obs ( 20.89 %)
#> 1 : 890 obs ( 79.11 %)
#>
#> ── X3 ──
#> 0 : 311 obs ( 27.64 %)
#> 1 : 814 obs ( 72.36 %)
#>
#> ── X4 ──
#> 0 : 525 obs ( 46.67 %)
#> 1 : 99 obs ( 8.8 %)
#> 2 : 501 obs ( 44.53 %)
#>
#> ── X5 ──
#> 0 : 700 obs ( 62.22 %)
#> 1 : 425 obs ( 37.78 %)
#>
#> ── X6 ──
#> 0 : 483 obs ( 42.93 %)
#> 1 : 642 obs ( 57.07 %)
#>
#> ── X7 ──
#> 0 : 1018 obs ( 90.49 %)
#> 1 : 107 obs ( 9.51 %)
#>
#> ── X8 ──
#> 0 : 907 obs ( 80.62 %)
#> 1 : 218 obs ( 19.38 %)
#>
#> ── X9 ──
#> 0 : 651 obs ( 57.87 %)
#> 1 : 474 obs ( 42.13 %)
#>
#> ── X10 ──
#> 0 : 180 obs ( 16 %)
#> 1 : 945 obs ( 84 %)
#>
#> ── X11 ──
#> 0 : 654 obs ( 58.13 %)
#> 1 : 471 obs ( 41.87 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 130.353
#> -2 Log Likelihood (Model Penuh) : 118.437
#> G² hitung : 11.916
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.37
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.817 1.122 0.0000 Signifikan ✅
#> X1 0.252 0.508 0.6618 Tdk Signifikan ❌
#> X2 -2.444 0.652 0.0013 Signifikan ✅
#> X3 0.024 0.598 0.9721 Tdk Signifikan ❌
#> X4 0.150 0.325 0.6880 Tdk Signifikan ❌
#> X5 0.718 0.533 0.2397 Tdk Signifikan ❌
#> X6 -0.012 0.550 0.9851 Tdk Signifikan ❌
#> X7 1.881 0.558 0.0050 Signifikan ✅
#> X8 0.462 0.543 0.4622 Tdk Signifikan ❌
#> X9 -0.872 0.632 0.2027 Tdk Signifikan ❌
#> X10 0.230 0.735 0.7760 Tdk Signifikan ❌
#> X11 1.709 0.779 0.0567 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 5.982
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.6492
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.008 [ 0.001, 0.079] Odds baseline
#> X1 1.287 [ 0.422, 4.424] ↑ lebih berisiko miskin
#> X2 0.087 [ 0.019, 0.387] ↓ lebih kecil risiko miskin
#> X3 1.024 [ 0.283, 4.518] ↑ lebih berisiko miskin
#> X4 1.162 [ 0.568, 2.542] ↑ lebih berisiko miskin
#> X5 2.051 [ 0.613, 6.958] ↑ lebih berisiko miskin
#> X6 0.988 [ 0.288, 3.592] ↓ lebih kecil risiko miskin
#> X7 6.563 [ 1.828, 22.622] ↑ lebih berisiko miskin
#> X8 1.588 [ 0.433, 5.105] ↑ lebih berisiko miskin
#> X9 0.418 [ 0.078, 1.548] ↓ lebih kecil risiko miskin
#> X10 1.259 [ 0.266, 7.593] ↑ lebih berisiko miskin
#> X11 5.524 [ 0.952, 34.664] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 10 — Sumbar + ART>1 | IND1 Filtered
model_ind1_filtered_art2_smb <- analisis_firth(data_ind1_filtered_art2_smb,
"Sumbar + ART>1 | IND1 Filtered", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND1 Filtered
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 650 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 639 obs ( 98.3 %)
#> Miskin (1): 11 obs ( 1.7 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 80 obs ( 12.31 %)
#> 1 : 209 obs ( 32.15 %)
#> 2 : 361 obs ( 55.54 %)
#>
#> ── X2 ──
#> 0 : 201 obs ( 30.92 %)
#> 1 : 449 obs ( 69.08 %)
#>
#> ── X3 ──
#> 0 : 212 obs ( 32.62 %)
#> 1 : 438 obs ( 67.38 %)
#>
#> ── X4 ──
#> 0 : 349 obs ( 53.69 %)
#> 1 : 65 obs ( 10 %)
#> 2 : 236 obs ( 36.31 %)
#>
#> ── X5 ──
#> 0 : 469 obs ( 72.15 %)
#> 1 : 181 obs ( 27.85 %)
#>
#> ── X6 ──
#> 0 : 305 obs ( 46.92 %)
#> 1 : 345 obs ( 53.08 %)
#>
#> ── X7 ──
#> 0 : 560 obs ( 86.15 %)
#> 1 : 90 obs ( 13.85 %)
#>
#> ── X8 ──
#> 0 : 543 obs ( 83.54 %)
#> 1 : 107 obs ( 16.46 %)
#>
#> ── X9 ──
#> 0 : 400 obs ( 61.54 %)
#> 1 : 250 obs ( 38.46 %)
#>
#> ── X10 ──
#> 0 : 65 obs ( 10 %)
#> 1 : 585 obs ( 90 %)
#>
#> ── X11 ──
#> 0 : 458 obs ( 70.46 %)
#> 1 : 192 obs ( 29.54 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 109.15
#> -2 Log Likelihood (Model Penuh) : 98.038
#> G² hitung : 11.112
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.4339
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.150 1.150 0.0006 Signifikan ✅
#> X1 0.098 0.499 0.8649 Tdk Signifikan ❌
#> X2 -2.035 0.635 0.0042 Signifikan ✅
#> X3 0.100 0.611 0.8868 Tdk Signifikan ❌
#> X4 0.156 0.341 0.6949 Tdk Signifikan ❌
#> X5 0.780 0.544 0.2157 Tdk Signifikan ❌
#> X6 -0.049 0.583 0.9429 Tdk Signifikan ❌
#> X7 1.736 0.545 0.0075 Signifikan ✅
#> X8 0.369 0.605 0.5994 Tdk Signifikan ❌
#> X9 -0.614 0.634 0.3816 Tdk Signifikan ❌
#> X10 -0.015 0.805 0.9863 Tdk Signifikan ❌
#> X11 1.601 0.760 0.0713 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 4.677
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.7915
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.016 [ 0.001, 0.170] Odds baseline
#> X1 1.103 [ 0.364, 3.758] ↑ lebih berisiko miskin
#> X2 0.131 [ 0.028, 0.531] ↓ lebih kecil risiko miskin
#> X3 1.105 [ 0.293, 5.009] ↑ lebih berisiko miskin
#> X4 1.168 [ 0.540, 2.645] ↑ lebih berisiko miskin
#> X5 2.182 [ 0.621, 7.459] ↑ lebih berisiko miskin
#> X6 0.952 [ 0.251, 3.749] ↓ lebih kecil risiko miskin
#> X7 5.677 [ 1.640, 19.172] ↑ lebih berisiko miskin
#> X8 1.447 [ 0.326, 5.219] ↑ lebih berisiko miskin
#> X9 0.541 [ 0.101, 2.032] ↓ lebih kecil risiko miskin
#> X10 0.985 [ 0.173, 6.605] ↓ lebih kecil risiko miskin
#> X11 4.956 [ 0.865, 28.555] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 11 — Sumbar + ART>1 | IND2 No-Duplikat
model_ind2_nodup_art2_smb <- analisis_firth(data_ind2_nodup_art2_smb,
"Sumbar + ART>1 | IND2 No-Duplikat", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND2 No-Duplikat
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 705 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 694 obs ( 98.4 %)
#> Miskin (1): 11 obs ( 1.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 95 obs ( 13.48 %)
#> 1 : 226 obs ( 32.06 %)
#> 2 : 384 obs ( 54.47 %)
#>
#> ── X2 ──
#> 0 : 226 obs ( 32.06 %)
#> 1 : 479 obs ( 67.94 %)
#>
#> ── X3 ──
#> 0 : 235 obs ( 33.33 %)
#> 1 : 470 obs ( 66.67 %)
#>
#> ── X4 ──
#> 0 : 387 obs ( 54.89 %)
#> 1 : 70 obs ( 9.93 %)
#> 2 : 248 obs ( 35.18 %)
#>
#> ── X5 ──
#> 0 : 506 obs ( 71.77 %)
#> 1 : 199 obs ( 28.23 %)
#>
#> ── X6 ──
#> 0 : 338 obs ( 47.94 %)
#> 1 : 367 obs ( 52.06 %)
#>
#> ── X7 ──
#> 0 : 598 obs ( 84.82 %)
#> 1 : 107 obs ( 15.18 %)
#>
#> ── X8 ──
#> 0 : 585 obs ( 82.98 %)
#> 1 : 120 obs ( 17.02 %)
#>
#> ── X9 ──
#> 0 : 438 obs ( 62.13 %)
#> 1 : 267 obs ( 37.87 %)
#>
#> ── X10 ──
#> 0 : 69 obs ( 9.79 %)
#> 1 : 636 obs ( 90.21 %)
#>
#> ── X11 ──
#> 0 : 503 obs ( 71.35 %)
#> 1 : 202 obs ( 28.65 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 110.951
#> -2 Log Likelihood (Model Penuh) : 99.732
#> G² hitung : 11.218
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.4252
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.354 1.144 0.0003 Signifikan ✅
#> X1 0.270 0.498 0.6356 Tdk Signifikan ❌
#> X2 -2.044 0.640 0.0046 Signifikan ✅
#> X3 0.140 0.601 0.8388 Tdk Signifikan ❌
#> X4 0.127 0.332 0.7418 Tdk Signifikan ❌
#> X5 0.717 0.544 0.2554 Tdk Signifikan ❌
#> X6 -0.082 0.556 0.8983 Tdk Signifikan ❌
#> X7 1.589 0.539 0.0133 Signifikan ✅
#> X8 0.294 0.597 0.6701 Tdk Signifikan ❌
#> X9 -0.587 0.633 0.4040 Tdk Signifikan ❌
#> X10 -0.088 0.797 0.9211 Tdk Signifikan ❌
#> X11 1.578 0.766 0.0780 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 6.48
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.5936
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.013 [ 0.001, 0.136] Odds baseline
#> X1 1.310 [ 0.435, 4.451] ↑ lebih berisiko miskin
#> X2 0.130 [ 0.028, 0.537] ↓ lebih kecil risiko miskin
#> X3 1.151 [ 0.312, 5.095] ↑ lebih berisiko miskin
#> X4 1.135 [ 0.537, 2.510] ↑ lebih berisiko miskin
#> X5 2.049 [ 0.582, 7.022] ↑ lebih berisiko miskin
#> X6 0.921 [ 0.259, 3.360] ↓ lebih kecil risiko miskin
#> X7 4.899 [ 1.422, 16.475] ↑ lebih berisiko miskin
#> X8 1.341 [ 0.308, 4.726] ↑ lebih berisiko miskin
#> X9 0.556 [ 0.104, 2.090] ↓ lebih kecil risiko miskin
#> X10 0.916 [ 0.163, 6.071] ↓ lebih kecil risiko miskin
#> X11 4.845 [ 0.831, 28.262] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Dataset 12 — Sumbar + ART>1 | IND2 Max-Age
model_ind2_maxage_art2_smb <- analisis_firth(data_ind2_maxage_art2_smb,
"Sumbar + ART>1 | IND2 Max-Age", drop_var_na = FALSE)#>
#> ╔══════════════════════════════════════════════════════╗
#> DATASET: Sumbar + ART>1 | IND2 Max-Age
#> MODEL : Firth Logistic Regression (Bias Reduction)
#> ╚══════════════════════════════════════════════════════╝
#>
#> ── Observasi valid: 706 baris
#>
#> ▶ 1. EDA
#> ────────────────────────────────────────────────────────
#> [Y] Distribusi status_miskin:
#> Tidak Miskin (0): 695 obs ( 98.4 %)
#> Miskin (1): 11 obs ( 1.6 %)
#>
#> [X] Distribusi & Proporsi Variabel:
#>
#> ── X1 ──
#> 0 : 95 obs ( 13.46 %)
#> 1 : 226 obs ( 32.01 %)
#> 2 : 385 obs ( 54.53 %)
#>
#> ── X2 ──
#> 0 : 226 obs ( 32.01 %)
#> 1 : 480 obs ( 67.99 %)
#>
#> ── X3 ──
#> 0 : 236 obs ( 33.43 %)
#> 1 : 470 obs ( 66.57 %)
#>
#> ── X4 ──
#> 0 : 387 obs ( 54.82 %)
#> 1 : 71 obs ( 10.06 %)
#> 2 : 248 obs ( 35.13 %)
#>
#> ── X5 ──
#> 0 : 507 obs ( 71.81 %)
#> 1 : 199 obs ( 28.19 %)
#>
#> ── X6 ──
#> 0 : 338 obs ( 47.88 %)
#> 1 : 368 obs ( 52.12 %)
#>
#> ── X7 ──
#> 0 : 599 obs ( 84.84 %)
#> 1 : 107 obs ( 15.16 %)
#>
#> ── X8 ──
#> 0 : 586 obs ( 83 %)
#> 1 : 120 obs ( 17 %)
#>
#> ── X9 ──
#> 0 : 438 obs ( 62.04 %)
#> 1 : 268 obs ( 37.96 %)
#>
#> ── X10 ──
#> 0 : 69 obs ( 9.77 %)
#> 1 : 637 obs ( 90.23 %)
#>
#> ── X11 ──
#> 0 : 504 obs ( 71.39 %)
#> 1 : 202 obs ( 28.61 %)
#>
#> ▶ 2. PEMODELAN FIRTH LOGISTIC REGRESSION
#> ────────────────────────────────────────────────────────
#> Model : logit[P(Y=1)] = β0 + β·X1 + β·X2 + β·X3 + β·X4 + β·X5 + β·X6 + β·X7 + β·X8 + β·X9 + β·X10 + β·X11
#> Metode: Penalized Maximum Likelihood (Firth, 1993)
#>
#> ▶ 3. UJI SIMULTAN (Penalized Likelihood Ratio Test / Uji G²)
#> ────────────────────────────────────────────────────────
#> H0: semua β = 0 vs H1: minimal ada 1 β ≠ 0
#> -2 Log Likelihood (Model Null) : 110.982
#> -2 Log Likelihood (Model Penuh) : 99.766
#> G² hitung : 11.216
#> df : 11
#> χ² tabel (α=5%, df= 11 ) : 19.675
#> p-value : 0.4253
#> ❌ KEPUTUSAN: Gagal Tolak H0 → Model tidak signifikan
#>
#> ▶ 4. UJI PARSIAL (Profile Penalized Likelihood Ratio Test)
#> ────────────────────────────────────────────────────────
#> H0: βj = 0 vs H1: βj ≠ 0 (α = 5%, χ² tabel = 3.841 )
#> Statistik uji: Profile Penalized LRT (bukan Wald)
#>
#> Variabel B SE p-value Kesimpulan
#> -----------------------------------------------------------------
#> (Intercept) -4.354 1.144 0.0003 Signifikan ✅
#> X1 0.270 0.498 0.6358 Tdk Signifikan ❌
#> X2 -2.044 0.640 0.0046 Signifikan ✅
#> X3 0.141 0.601 0.8384 Tdk Signifikan ❌
#> X4 0.127 0.333 0.7418 Tdk Signifikan ❌
#> X5 0.717 0.544 0.2553 Tdk Signifikan ❌
#> X6 -0.082 0.556 0.8981 Tdk Signifikan ❌
#> X7 1.589 0.539 0.0133 Signifikan ✅
#> X8 0.294 0.597 0.6701 Tdk Signifikan ❌
#> X9 -0.587 0.633 0.4037 Tdk Signifikan ❌
#> X10 -0.088 0.797 0.9211 Tdk Signifikan ❌
#> X11 1.579 0.765 0.0778 Tdk Signifikan ❌
#>
#> ▶ 5. GOODNESS OF FIT (Hosmer-Lemeshow Test)
#> ────────────────────────────────────────────────────────
#> H0: Model fit dengan data
#> χ² Hosmer-Lemeshow : 6.482
#> df : 8
#> χ² tabel (df=8) : 15.507
#> p-value : 0.5934
#> ✅ KEPUTUSAN: Gagal Tolak H0 → Model FIT dengan data
#>
#> ▶ 6. INTERPRETASI PARAMETER (Odds Ratio / exp(B) — Profile Penalized CI)
#> ────────────────────────────────────────────────────────
#> Variabel exp(B) 95% CI Interpretasi Singkat
#> ---------------------------------------------------------------------------
#> (Intercept) 0.013 [ 0.001, 0.136] Odds baseline
#> X1 1.310 [ 0.435, 4.450] ↑ lebih berisiko miskin
#> X2 0.129 [ 0.028, 0.537] ↓ lebih kecil risiko miskin
#> X3 1.151 [ 0.312, 5.096] ↑ lebih berisiko miskin
#> X4 1.135 [ 0.537, 2.510] ↑ lebih berisiko miskin
#> X5 2.049 [ 0.582, 7.023] ↑ lebih berisiko miskin
#> X6 0.921 [ 0.259, 3.359] ↓ lebih kecil risiko miskin
#> X7 4.900 [ 1.422, 16.478] ↑ lebih berisiko miskin
#> X8 1.341 [ 0.308, 4.727] ↑ lebih berisiko miskin
#> X9 0.556 [ 0.104, 2.089] ↓ lebih kecil risiko miskin
#> X10 0.916 [ 0.163, 6.071] ↓ lebih kecil risiko miskin
#> X11 4.848 [ 0.832, 28.270] ↑ lebih berisiko miskin
#>
#> ══════════════════════════════════════════════════════Catatan Akhir: Seluruh analisis dalam dokumen ini menggunakan data Survei Sosial Ekonomi Nasional (Susenas) Maret 2025. Garis Kemiskinan yang digunakan merupakan GK per kabupaten/kota dari BPS. Model Firth direkomendasikan sebagai model utama mengingat karakteristik data yang menunjukkan ketidakseimbangan kelas pada variabel dependen
status_miskin. Packagelogistf(Heinze & Schemper, 2002) digunakan untuk implementasi di R.