ABSTRAK

Kemiskinan adalah cerminan dari kecilnya pendapatan seseorang sehingga standar hidup yang layak tidak dapat diakses olehnya. Kemisikinan telah menjadi tantangan global sejak bertahun-tahun lamanya, termasuk bagi pemerintah di Indonesia. Persentase kemiskinan di Jawa Barat tahun 2024 sendiri menduduki urutan ketiga tertinggi di Indonesia. Setelah dilakukan pemetaan persentase kemiskinan, terlihat adanya indikasi dependensi spasial antarwilayah. Tujuan penelitian ini adalah memperoleh model terbaik yang dapat menjelaskan hubungan antarwilayah serta faktor sosial-ekonomi utama yang berpengaruh terhadap kemiskinan di Jawa Barat. Data yang digunakan pada penelitian ini bersumber dari Badan Pusat Statistik dimana Persentase Kemiskinan menjadi variabel respons, lalu Tingkat Pengangguran Terbuka (TPT), Produk Domestik Regional Bruto (PDRB) per kapita, Indeks Pembangunan Manusia (IPM), dan Gini Ratio menjadi variabel prediktornya. Berdasarkan hasil analisis, didapatkan model terbaik yaitu SAR. Variabel yang signfikan pada taraf signifikansi 5% adalah IPM, dimana setiap kenaikan satu unit pada IPM akan menurunkan persentase penduduk miskin sebesar 0.4132%. Sedangkan pada taraf signfikansi 1%, setiap kenaikan satu unit pada Gini Ratio akan meningkatkan persentase penduduk miskin sebesar 12.705%. Selain itu, variabel lain belum menunjukkan signfikansi yang berarti pada penelitian ini.

BAB 1 PENDAHULUAN

1.1 Latar Belakang Masalah

Kemiskinan menunjukkan kecilnya pendapatan sehingga tidak memungkinkan bagi seseorang untuk mendapatkan standar hidup yang layak (Corral et al., 2020). Kemiskinan yang berkepanjangan akan berdampak pada berbagai aspek, diantaranya melebarnya ketimpangan sosial, menghambat pembangunan, atau bahkan memicu konflik horizontal dan meningkatkan tingkat kriminalitas suatu negara (UN, 2023). Apabila tidak ditangani dengan tepat, kemiskinan akan membentuk siklus antargenerasi yang sulit diselesaikan (Rivera & Currais (1999); Benhabib & Spiegel (1994; Lucas (1990)).

Meskipun tingkat kemiskinan di Indonesia, khususnya Jawa Barat, terus menurun hingga 7.46% pada tahun 2024, persentase masyarakat miskin di Jawa Barat sejak 2020-2024 cenderung stagnan.

Hal ini menunjukkan adanya urgensi untuk memahami faktor-faktor pengaruh kemiskinan dan pola spasialnya demi mendukung perencanaan pembangunan yang menyeluruh dan pengentasan kemiskinan di Jawa Barat.

1.2 Identifikasi Masalah

Pemetaan persentase kemiskinan di kabupaten/kota Provinsi Jawa Barat menunjukkan adanya indikasi tidak independen secara spasial, dimana nilai kemiskinan di satu kab/kota memiliki keterkaitan dengan nilai observasi di kab/kota lain yang berdekatan. Kondisi ini menunjukkan adanya spatial dependence.

Adanya dependensi spasial ini mengindikasikan bahwa asumsi independensi antar observasi yang digunakan pada model regresi klasik (OLS) tidak terpenuhi. Jika kondisi tersebut diabaikan, maka hasil estimasi dapat menjadi bias dan tidak efisien. Oleh karena itu, perlu digunakan pendekatan ekonometrika spasial untuk menangkap pengaruh spasial secara lebih akurat.

1.3 Maksud dan Tujuan Penelitian

Maksud dari penelitian ini adalah untuk menerapkan pendekatan spasial untuk menganalisis faktor-faktor yang mempengaruhi kemiskinan yang terjadi di Jawa Barat pada tahun 2024 secara signifikan. Melalui analisis spasial, penelitian ini bertujuan untuk memperoleh model terbaik yang dapat menjelaskan hubungan antarwilayah serta faktor-faktor utama yang berpengaruh terhadap kemiskinan di Jawa Barat.

1.4 Batasan Penelitian

Penelitian ini dibatasi pada ruang lingkup wilayah Provinsi Jawa Barat dengan unit analisis kabupaten/kota. Analisis difokuskan pada pola spasial kemiskinan tanpa membahas aspek penyebab non-spasial secara mendalam. Selain itu, penelitian ini hanya menggunakan data tahun 2024 sehingga hasilnya menggambarkan kondisi pada periode tertentu dan tidak mencakup dinamika perubahan antarwaktu.

BAB 2 TINJAUAN PUSTAKA

2.1 Spatial Autocorrelation

Dalam analisis spasial, keterkaitan antar wilayah secara spasial diukur menggunakan Moran’s I. Nilai I > 0 menunjukkan adanya autokorelasi spasial positif, yaitu wilayah dengan nilai mirip (high-high atau low-low) cenderung berdekatan. Sebaliknya, I < 0 menandakan pola checkerboard (high-low), sedangkan I \(\approx\) 0 menunjukkan pola yang acak. Dalam perhitungannya, Moran’s I memanfaatkan matriks bobot spasial (W) yang menentukan seberapa kuat pengaruh satu wilayah terhadap wilayah lain berdasarkan kriteria kedekatan tertentu, yakni Rook, Queen, ataupun Jarak.

2.2 Ordinary Least Square (OLS)

\[Y = X\beta + \epsilon\]

OLS adalah metode estimasi parameter untuk regresi klasik. Metode ini membutuhkan pemenuhan pengujian asumsi klasik, diantaranya pengujian multikolinearitas, residual berdistribusi normal, dan homoskedastisitas. Metode ini mengasumsikan bahwa error antarwilayah bersifat independen dan homogen, sehingga tidak dapat menangkap spillover atau pengaruh kemiskinan dari satu daerah ke daerah lain. Hal ini menyebabkan hasil estimasi menjadi bias dan kurang akurat. Untuk mengatasi keterbatasan tersebut, digunakan model ekonometrika spasial yang mampu memperhitungkan spatial dependence antarwilayah.

2.3 Uji Lagrange Multiplier

Uji ini digunakan untuk mengidentifikasi ada tidaknya dependensi spasial pada model. Jika hasil uji menunjukkan adanya dependensi lag, maka digunakan model Spatial Autoregressive (SAR). Jika yang muncul adalah dependensi galat, maka digunakan model Spatial Error Model (SEM). Uji ini membantu menentukan model spasial yang paling sesuai sebelum analisis lanjut.

2.4 Spatial Autoregressive (SAR)

Menurut Anselin (2003), Spatial Autoregressive Model (SAR) merupakan pengembangan dari model OLS dengan menambahkan pengaruh lag pada variabel dependen. Model ini digunakan ketika \(\lambda = 0\), sehingga hanya mempertimbangkan proses autoregresif pada variabel respon. Dengan kata lain, SAR menunjukkan bahwa nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah-wilayah tetangganya, sesuai dengan hubungan spasial yang terbentuk dalam model.

\[Y = \rho Wy + X\beta + \epsilon\]

Asumsi penting dalam model SAR adalah bahwa residual tidak mengalami autokorelasi spasial, sehingga pengaruh spasial telah sepenuhnya ditangkap oleh komponen lag variabel dependen. Jika asumsi ini terpenuhi, maka model SAR dianggap tepat untuk menggambarkan struktur keterkaitan spasial dalam data.

2.5 AIC

Menurut Suryaningrat, S. H. (2021) Akaike Information Criterion (AIC) merupakan metode analisis yang digunakan untuk memperoleh model yang terbaik dengan menggunakan estimasi maximum likelihood sebagai perhitungan yang sesuai.. AIC mempertimbangkan keseimbangan antara kompleksitas model (jumlah parameter) dan kecocokan model terhadap data (likelihood). AIC dihitung dengan rumus: \[AIC=2k-2ln(\hat{L})\] di mana k adalah jumlah parameter dalam model dan \(\hat{L}\) adalah nilai maksimum dari fungsi likelihood. Model dengan nilai AIC yang lebih rendah dianggap lebih baik karena menunjukkan bahwa model tersebut memiliki kompleksitas yang wajar tetapi tetap mampu menjelaskan data dengan baik.

BAB 3 METODOLOGI PENELITIAN

3.1 Data

Penelitian ini menggunakan data sekunder yang diperoleh dari Badan Pusat Statistik (BPS) tahun 2024 dengan unit analisis sebanyak 27 kabupaten/kota di Provinsi Jawa Barat. Data yang digunakan bersifat cross section dan mencerminkan kondisi sosial ekonomi daerah yang diduga berpengaruh terhadap tingkat kemiskinan.

Variabel yang digunakan terdiri atas satu variabel dependen yaitu kemiskinan, yang diukur melalui persentase penduduk miskin menurut kabupaten/kota, serta empat variabel independen yaitu Tingkat Pengangguran Terbuka (TPT), Produk Domestik Regional Bruto (PDRB) per kapita, Indeks Pembangunan Manusia (IPM), dan Gini Ratio. Keempat variabel tersebut dipilih karena secara teori memiliki hubungan erat dengan tingkat kesejahteraan dan ketimpangan sosial ekonomi masyarakat.

Selain data atribut, penelitian ini juga menggunakan data spasial berupa peta batas administrasi kabupaten/kota di Jawa Barat yang diperoleh dari geoBoundaries Indonesia Administrative Level 2 (ADM2) melalui situs www.geoboundaries.org dalam format shapefile (.shp). Data ini digunakan untuk membentuk hubungan antarwilayah dalam analisis spasial.

Tabel 1. Definisi Operasional Variabel
Variabel Simbol Definisi Operasional Satuan
Kemiskinan Y Persentase penduduk miskin yang berada di bawah garis kemiskinan. %
TPT X1 Persentase jumlah pengangguran terhadap angkatan kerja. %
PDRB X2 Nilai total produksi barang dan jasa wilayah per jumlah penduduk berdasarkan harga berlaku. Juta Rupiah
IPM X3 Ukuran yang menunjukkan tingkat pencapaian pembangunan manusia di suatu wilayah. Indeks
Gini Ratio X4 Ukuran kesenjangan distribusi pendapatan antar kelompok penduduk. Indeks

3.2 Metode Analisis

Analisis dilakukan menggunakan pendekatan statistik spasial dengan bantuan perangkat lunak R. Pendekatan ini bertujuan untuk mengidentifikasi pengaruh faktor sosial ekonomi terhadap tingkat kemiskinan serta mendeteksi adanya ketergantungan spasial antarwilayah di Jawa Barat. Analisis diawali dengan uji autokorelasi spasial global (Moran’s I, Geary’s C) dan lokal (LISA serta Getis-Ord Gi*) untuk melihat pola keterkaitan antarwilayah.

Selanjutnya digunakan model regresi OLS sebagai dasar analisis hubungan antarvariabel, diikuti uji Lagrange Multiplier (LM) untuk menentukan model spasial yang tepat. Beberapa model yang dipertimbangkan mencakup SAR, SEM, SDM, SDEM, SAC/SARAR, dan GNS. Pemilihan model terbaik dilakukan berdasarkan nilai AIC, signifikansi parameter spasial (ρ dan λ), serta kesesuaian teori yang mendukung hubungan spasial antarwilayah.

3.3 Alur Kerja Penelitian

Penelitian ini dilakukan dengan langkah-langkah yang ditampilkan dalam diagram alir berikut.

BAB 4 HASIL DAN PEMBAHASAN

4.1 Analisis Deskriptif

Analisis deskriptif dilakukan untuk memberikan gambaran umum mengenai karakteristik data penelitian. Hasil analisis deskriptif data dapat dilihat pada tabel berikut:

Tabel 2. Hasil Analisis Deskriptif
Variabel Mean Minimum Maximum Simpangan_Baku
Kemiskinan 8.01 2.34 11.93 2.65
TPT 6.52 1.58 8.97 1.71
PDRB 53.11 24.91 147.08 33.26
IPM 74.68 68.89 83.75 4.33
Gini Ratio 0.37 0.29 0.48 0.05

Selanjutnya, sebagai gambaran awal dilakukan pula pemetaan untuk penyebaraan masing-masing variabel penelitian sebagai berikut.

Penyebaran variabel pada peta menunjukkan adanya pola spasial yang saling berkaitan antarwilayah di Jawa Barat. Daerah dengan tingkat kemiskinan tinggi cenderung berkelompok di bagian timur laut, yang juga memiliki tingkat pengangguran (TPT) tinggi dan PDRB per kapita rendah. Wilayah selatan dengan IPM lebih rendah menunjukkan kemiskinan yang cukup tinggi. Pola serupa juga terlihat pada Rasio Gini, dimana ketimpangan relatif tinggi muncul di wilayah yang sama dengan PDRB tinggi. Kecenderungan pengelompokan nilai-nilai serupa ini mengindikasikan adanya dependensi spasial, sehingga observasi antarwilayah tidak bersifat independen.

4.2 Ordinary Least Square (OLS)

Hasil analisis dengan R menghasilkan persamaan model regresi linear berikut.

\[ Y = 44.770 - 0.0879X_1 - 0.0046X_2 - 0.5521X_3^{***} + 14.2366X_4 \] Model yang terbentuk memiliki nilai R Squared sebesar 57.82% dengan p-value uji signifikansi simultan sebesar < 0.001 yang artinya variabel TPT, PDRB, Gini Ratio dan IPM secara bersama-sama berpengaruh terhadap Kemiskinan. Adapun pengaruh parsial dari IPM signifikan pada taraf signifikansi 5%, sedangkan pengaruh dari variabel lainnya tidak signifikan. Untuk menghasilkan model yang valid dan reliabel, model OLS perlu diuji asumsinya yang hasilnya disajikan pada Tabel 3. berikut.

Tabel 3. Hasil Uji Asumsi Klasik
Jenis Uji Statistik Statistik Hitung p-value Kesimpulan
Linearitas (Ramsey RESET) RESET 0.47600 0.6281 variabel independen terhadap variabel dependen membentuk model linier
Normalitas (Shapiro Wilk) W 0.95524 0.2865 residual berdistribusi normal
Homoskedastisitas (Breusch–Pagan) BP 2.56000 0.6339 tidak terdapat gejala heteroskedastisitas

Lalu untuk pengujian asumsi multikolinearitas, dilakukan dengan melihat nilai VIF. Didapatkan nilai VIF setiap variabel x < 10, yang artinya tidak terjadi multikolinearitas pada variabel independen. Diperoleh kesimpulan bahwa seluruh asumsi klasik pada model OLS terpenuhi.

Selanjutnya dilakukan Uji Moran’s I residual untuk melihat apakah residual pad model regresi OLS bersifat acak atau tidak. Hasil pengujian ditampilkan pada tabel berikut:

Tabel 4. Hasil Uji Moran’s I
I p (Monte Carlo) Keputusan
2.3277 0.015 Tolak H0

4.3 Autokorelasi Spasial

Uji autokorelasi spasial dilakukan pada setiap variabel untuk mengetahui ada tidaknya ketergantungan spasial, baik pada variabel dependen maupun variabel independen.

4.3.1 Moran’s I

Tabel 5. Hasil Uji Moran’s I per Variabel
Variabel Moran’s I Z p (Monte Carlo) Keputusan
Kemiskinan 0.283 3.65 0.001 Signifikan
TPT 0.511 4.07 0.001 Signifikan
PDRB 0.211 1.86 0.102 Tidak Signifikan
IPM 0.275 2.24 0.032 Signifikan
Gini Ratio 0.299 2.40 0.036 Signifikan

Berdasarkan hasil uji Moran’s I pada Tabel 5., variabel Kemiskinan, TPT, IPM, dan Gini Ratiomenunjukkan adanya autokorelasi spasial yang signifikan. Artinya, nilai variabel-variabel tersebut di suatu wilayah cenderung dipengaruhi oleh wilayah sekitarnya. Sebaliknya, variabel PDRB tidak menunjukkan pola spasial yang signifikan.

4.3.2 Geary’s C

Tabel 6. Hasil Uji Geary’s C per Variabel
Variabel Geary’s C Z p (Monte Carlo) Keputusan
Kemiskinan 0.410 3.76 0.001 Signifikan
TPT 0.465 2.92 0.001 Signifikan
PDRB 0.660 1.82 0.031 Signifikan
IPM 0.607 2.42 0.005 Signifikan
GINI 0.504 3.07 0.004 Signifikan

Hasil menunjukkan bahwa dengan taraf signifikansi 5%, seluruh variabel dependen maupun independen memiliki autokorelasi spasial positif yang signifikan.

4.3.3 Local Moran’s I (LISA)

Hasil analisis Local Moran’s I (LISA) pada Gambar 4. menunjukkan bahwa variabel kemiskinanmembentuk klaster High-High di wilayah Cirebon dan Majalengka serta klaster Low-Low di wilayah Bogor dan Kota Bekasi. Variabel TPT, IPM, dan Gini Ratio juga menunjukkan pola pengelompokan lokal, sedangkan PDRB tidak membentuk klaster signifikan.

4.3.4 Getis-Ord \(G_{i}^{*}\)

Hasil analisis Getis–Ord G*i pada Gambar 5. menunjukkan pola hotspot dan coldspot yang konsisten dengan hasil LISA. Wilayah Cirebon dan Majalengka teridentifikasi sebagai hotspot, yang berarti tingkat kemiskinannya tinggi dan dikelilingi oleh wilayah dengan tingkat kemiskinan tinggi. Sebaliknya, Bogor dan Kota Bekasi teridentifikasi sebagai coldspot.

4.4 Uji Lagrange Multiplier

Tabel 7. Hasil Uji Lagrange Multiplier
Statistik Uji Value p-value Keputusan
LM Error (RSer) 3.4944 0.0616 Tidak Signifikan
LM Lag (RSlag) 8.8489 0.0029 Signifikan
Robust LM Error (AdjRSer) 1.0849 0.2976 Tidak Signifikan
Robust LM Lag (AdjRSlag) 6.4394 0.0112 Signifikan
LM SARMA (Gabungan) 9.9338 0.0070 Signifikan

Berdasarkan hasil uji Lagrange Multiplier (LM), diketahui bahwa nilai LM Lag (p = 0.0029) dan Adjusted LM Lag (p = 0.0112) signifikan pada α = 0.05, sedangkan LM Error dan Adjusted LM Error tidak signifikan. Hal ini menunjukkan adanya pengaruh spasial pada variabel dependen (lag spasial), bukan pada error-nya. Dengan demikian, salah satu model spasial yang paling tepat digunakan adalah model Spatial Autoregressive (SAR). Nilai LM SARMA (p = 0.007) juga signifikan, menunjukkan adanya dependensi spasial secara keseluruhan dalam model.

4.5 Perbandingan Model

Hasil uji Moran’s I pada residual model OLS menunjukkan adanya autokorelasi spasial pada error, sementara hasil uji LM-Lag (Adjusted) menunjukkan signifikansi, sedangkan uji LM-Error (Adjusted) tidak signifikan. Hal ini mengindikasikan adanya endogenous spatial dependence pada variabel dependen. Berdasarkan temuan ini, dapat disimpulkan bahwa model spasial yang tidak mengakomodasi adanya endogenous spatial dependence (seperti SLX, SEM, dan SDEM) tidak sesuai untuk digunakan. Oleh karena itu, empat model kandidat yang lebih cocok untuk analisis ini adalah SDM, GNS, SAC/SARAR, dan SAR. Adapun hasil pemodelan disajikan dalam Tabel 8. berikut ini.

Tabel 8. Perbandingan Model Spasial
Model AIC BIC p uji ρ p uji λ p PDRB p GINI p TPT p IPM
SAR 101.953 111 4e-04
0.294 0.097 0.51 2.67e-07
SDM 102.327 117 0.005
0.694* 0.473* 0.141* 0.776*
GNS 103.903 119 0.59 0.517 0.904* 0.585* 0.116* 0.567*
SAC/SARAR 102.898 113 2.9e-07 0.361 0.388 0.101 0.596 1.33e-07

4.6 Estimasi Model SAR

Berdasarkan hasil perbandingan model, model Spatial Autoregressive (SAR) dipilih sebagai model terbaik. Dengan menggunakan metode Maximum Likelihood Estimation (MLE), model SAR untuk kemiskinan di Jawa Barat dapat dirumuskan sebagai:

\[ y_i = \rho \sum_{j=1}^{n} w_{ij} y_j + 28.449^{***} + 0.1132 (TPT)_i - 0.0087 (PDRB)_i - 0.4132 (IPM)_i^{***} + 12.705 (GINI)_i^{*} \]

Adapun hasil uji Likelihood Ratio disajikan pada Tabel 9. berikut.

Tabel 9. Hasil Uji LR (Likelihood Ratio)
ρ LR test value p-value Keputusan
0.65269 12.428 0.00042 Signifikan

Berdasarkan hasil uji Likelihood Ratio, nilai ρ sebesar 0.65269 menunjukkan adanya hubungan positif antara kemiskinan di suatu wilayah dengan kemiskinan di wilayah sekitarnya. Nilai p-value yang lebih kecil dari α = 0.05, mengindikasikan bahwa pengaruh spasial dalam model ini signifikan. Sehingga dapat disimpulkan bahwa jika tingkat kemiskinan di suatu kabupaten/kota meningkat, maka kemiskinan di kabupaten/kota tetangga cenderung ikut meningkat. Hasil perhitungan impact pada model SAR disajikan pada Tabel 12. berikut.

Tabel 10. Hasil Estimasi Direct dan Indirect Impact
Variabel Jenis Impact Impact Measures p-value Keputusan
IPM Direct -0.4810 0.018 Signifikan
IPM Indirect -0.7080 0.800 Tidak Signifikan
TPT Direct 0.1320 0.553 Tidak Signifikan
TPT Indirect 0.1940 0.874 Tidak Signifikan
PDRB Direct -0.0101 0.351 Tidak Signifikan
PDRB Indirect -0.0150 0.822 Tidak Signifikan
Gini Ratio Direct 14.8060 0.232 Tidak Signifikan
Gini Ratio Indirect 21.7740 0.855 Tidak Signifikan

Berdasarkan hasil perhitungan impact, variabel IPM menunjukkan pengaruh langsung yang signifikan terhadap tingkat kemiskinan. Hal ini mengindikasikan bahwa peningkatan IPM di suatu wilayah akan menurunkan tingkat kemiskinan di wilayah tersebut. Namun, pengaruh tidak langsung dari IPM terhadap kemiskinan tidak signifikan. Hal ini mungkin terjadi karena spillover dari wilayah tetangga tidak cukup kuat atau tidak terdeteksi dengan baik dalam model ini. Selain itu, pengaruh IPM terhadap kemiskinan bisa lebih terasa langsung daripada melalui jalur spasial atau mediasi, yang menjelaskan mengapa efek tidak langsungnya tidak signifikan.

Variabel Gini Ratio, PDRB, dan TPT tidak signifikan baik secara langsung maupun tidak langsung karena kemungkinan pengaruhnya terhadap kemiskinan tidak bersifat linier atau terhambat oleh faktor lain yang tidak tercakup dalam model ini. Misalnya, ketimpangan pendapatan (Gini) dan tingkat pengangguran (TPT) mungkin dipengaruhi oleh faktor-faktor sosial atau kebijakan yang lebih kompleks yang tidak termodel dengan baik, sementara pengaruh PDRB per kapita mungkin lebih terasa dalam jangka panjang, bukan secara langsung di tahun 2024.

4.6 Uji Asumsi Model SAR

Untuk memastikan apakah model yang dibentuk dapat diandalkan dilakukan uji asumsi yang hasilnya sebagai berikut.

Tabel 11. Hasil Uji Diagnostik Model
Jenis Uji Statistik Statistik Hitung p-value Keputusan
Residual Autokorelasi (Moran’s I) I -0.11100 0.7077 tidak terdapat autokorelasi spasial pada residual
Linearitas (Ramsey RESET) RESET 0.02100 0.9790 terdapat pola linear antara variabel dependen dan independen dalam model
Normalitas (Shapiro Wilk) W 0.95524 0.2865 residual berdistribusi normal
Homoskedastisitas (Breusch-Pagan) BP 3.01690 0.5550 tidak terdapat heteroskedastisitas

Berdasarkan hasil uji asumsi, diperoleh kesimpulan bahwa seluruh asumsi model terpenuhi sehingga model SAR yang telah diestimasi dapat dianggap valid dan dapat digunakan untuk menganalisis pengaruh faktor ekonomi dan sosial terhadap tingkat kemiskinan di Jawa Barat.

4.7 Uji Sensivitas Bobot

Tabel 12. Hasil Perbandingan Berdasarkan Bobot Spasial
Bobot ρ AIC LogLik
Queen 0.6527 101.9525 -43.97625
Rook 0.6527 101.9525 -43.97625
KNN-4 0.4871 105.5565 -45.77824

Uji sensitivitas bobot spasial dilakukan untuk melihat stabilitas dari model SAR sebelumnya. Berdasarkan tabel bobot di atas, diketahui bahwa Queen dan Rook memberikan nila , AIC, dan hasil ρ log-likelihood yang sama, meskipun terjadi perbedaan pada bobot KNN-4. Hal ini menunjukkan bahwa spatial dependence di Jawa Barat lebih dipengaruhi oleh kedekatan geografis administratif dan bukan jarak nyata.

4.8 Interpretasi Hasil

Berdasarkan peta prediksi, terlihat bahwa Kota Bandung, Kota Bekasi, dan Depok diprediksi memiliki tingkat kemiskinan yang rendah, sementara Kabupaten Bekasi diprediksi dengan nilai menengah-rendah. Sebagian besar wilayah lainnya menunjukkan prediksi dengan tingkat kemiskinan menengah-atas. Hasil ini menggambarkan pola kemiskinan yang lebih merata di sebagian besar wilayah, dengan konsentrasi kemiskinan yang relatif rendah di daerah perkotaan besar seperti Bandung, Bekasi, dan Depok.

Namun, pada peta residual, beberapa wilayah menunjukkan adanya ketidaksesuaian antara nilai prediksi dan nilai aktual. Kota Tasikmalaya, misalnya, menunjukkan perbedaan yang cukup besar antara prediksi (7.782) dan nilai aktual (11.100), yang menandakan bahwa model ini tidak sepenuhnya akurat di wilayah tersebut. Kabupaten Sukabumi juga menunjukkan deviasi yang signifikan dengan nilai z-score sebesar -2,391, menandakan underprediction yang cukup besar. Di sisi lain, Kabupaten Indramayu menunjukkan residual yang lebih moderat dengan nilai z-score 1,65, menunjukkan bahwa prediksinya cukup mendekati nilai aktual. Sementara itu, untuk sebagian besar wilayah lainnya, distribusi residual tampak lebih merata dan akurat, mengindikasikan bahwa model SAR memberikan prediksi yang lebih stabil dan sesuai untuk sebagian besar kabupaten/kota di Jawa Barat.

Model SAR yang telah terbentuk untuk kemiskinan di Jawa Barat 2024 secara matematis dapat dituliskan sebagai:

\[ y_i = \rho \sum_{j=1}^{n} w_{ij} y_j + 28.449^{***} + 0.1132 (TPT)_i - 0.0087 (PDRB)_i - 0.4132 (IPM)_i^{***} + 12.705 (GINI)_i^{*} \]

Hasil menunjukkan bahwa IPM berpengaruh signifikan pada taraf signifikansi 5%, Gini Ratio signifikan pada taraf signifikansi 1%, sedangkan TPT, PDRB tidak berpengaruh signifikan terhadap Persentase Penduduk Miskin. Berdasarkan nilai koefisien hasil estimasi, setiap kenaikan satu unit pada IPM akan menurunkan persentase penduduk miskin sebesar 0.4132%, yang menunjukkan bahwa perbaikan dalam kualitas hidup dan pembangunan manusia dapat mengurangi tingkat kemiskinan. Di sisi lain, setiap kenaikan satu unit pada Gini Ratio akan meningkatkan persentase penduduk miskin sebesar 12.705%, yang mengindikasikan bahwa semakin tinggi ketimpangan pendapatan, semakin besar pula tingkat kemiskinan. Hasil ini menunjukkan bahwa IPM berkontribusi untuk mengurangi kemiskinan, sementara GINI Ratio justru berkontribusi pada peningkatan kemiskinan.

BAB 5 KESIMPULAN

5.1 Kesimpulan

Berdasarkan hasil analisis yang telah dilakukan untuk melihat faktor ekonomi-spasial untuk menganalisis persentase kemiskinan di Jawa Barat pada Tahun 2024, didapatkan model terbaik yaitu SAR. Berdasarkan model terbaik tersebut, disarankan kepada pemerintah ataupun stake holder terkait untuk fokus pada peningkatan Indeks Pembangunan Manusia (IPM), karena IPM sendiri sudah meliputi ketiga dimensi kehidupan dasar yaitu Pendidikan, Kesehatan, dan Standar Hidup Layak. Langkah perbaikan dengan fokus kepada IPM dinilai sebagai langkah yang bijak karena dapat memberikan dampak yang positif terhadap penurunan persentase kemisikanan yang terjadi di Jawa Barat.

Selain itu, langkah untuk mengurangi atau memperpendek ketimpangan ekonomi antar masyarakat yang kaya dan miskin juga dinilai sebagai langkah yang cukup menjanjikan, karena semakin sempit jarak ketimpangan, peningkatan masyarakat yang tergolong miskin akan semakin rendah pula. Oleh karena itu, berdasarkan penelitian ini dapat disimpulkan bahwa yang dapat mempengaruhi dan yang seharusnya menjadi fokus pemerintah dalam mengentaskan kemiskinan di Jawa Barat adalah aspek Pendidikan, Kesehatan, Ekonomi (dengan menjadikan adanya standar hidup layak).

5.2 Saran

Untuk penelitian selanjutnya, disarankan untuk melakukan eksplorasi variabel lebih banyak untuk melihat faktor-faktor yang mempengaruhi persentase kemiskinan di Jawa Barat. Hal ini dapat dilakukan dengan melihat faktor sosial, ekonomi, budaya, kesehatan, pemerintah setempat, yang mungkin saja memiliki pengaruh yang belum pernah diungkap sebelumnya. Pendekatan spasial dapat terus dikembangkan untuk penelitian berikutnya seperti diterapkan pada daerah lain yang memiliki pola atau fenomena spasial yang lebih merinci.

DAFTAR PUSTAKA

Corral. P. et al.. 2020. Fragility and Conflict: On the Front Lines of the Fight against Poverty. Washington DC: World Bank Publications

United Nations Development Programme. 2022. Sustainable Development Goals. [Online] Available at: https://www.undp.org/sustain able-development-goals [Accessed 13 OKT 2024].

LeSage, J. P., & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press. DOI: 10.1201/9781420064244

Anselin, L. (2001). Spatial Econometrics. In B. H. Baltagi (Ed.), A Companion to Theoretical Econometrics (pp. 310–330). Blackwell Publishing. DOI: 10.1002/9780470996249.ch15

Ayudia, N., & Putri, D. (2024). Pemodelan Spatial Autoregressive Model (SAR) pada Kasus Kemiskinan di Jawa Timur. Jurnal Gaussian, 13(2), 308–318. DOI: 10.14710/j.gauss.13.2.308-318

Noviyanti, N. K., Setiawan, M. Y., & Wijayanti, S. N. (2023). Spatial Regression and Spatial Autocorrelation Analysis of the Determinants of Poverty in Indonesia in 2022. Jurnal Ilmiah Aplikatif: Statistika, 7(2), 273–288. DOI: 10.31092/jia.v7i2.2318

Arsyad, A. M., & Riani, P. (2020). Analisis Faktor-faktor yang Mempengaruhi Kemiskinan di Jawa Barat Menggunakan Regresi Data Panel Spasial. Jurnal Statistik: Jurnal Ilmu Statistik dan Aplikasi, 7(1), 25–36. DOI: 10.29244/jst.v7i1.621

Curtis, K. J., Dorius, S. F., & Haghish, E. F. (2018). Spatial Regression Analysis of Poverty in R. Journal of Geographical Systems, 20(4), 383–404. DOI: 10.1007/s41322-018-0081-3

Widiyanto, P. S., Rahmawati, Y., & Rahayu, S. M. (2021). Spatial Econometric Model for Mapping Poverty Area in West Sulawesi. Journal of Physics: Conference Series, 1752(1), 012048. DOI: 10.1088/1742-6596/1752/1/012048

Purwono, R., & Prawoto, N. (2019). Analisis Spasial Faktor-Faktor yang Mempengaruhi Kemiskinan di Indonesia. Jurnal Ekonomi dan Bisnis, 22(2), 195–211. DOI: 10.24914/jeb.v22i2.2743

PSKK UGM. (2018, 8 Juni). Upaya penanggulangan kemiskinan dari masa ke masa. Diakses dari https://cpps.ugm.ac.id/upaya-penanggulangan-kemiskinan-dari-masa-ke-masa/

Haryanto, T. (2020). Faktor-Faktor yang Mempengaruhi Tingkat Kemiskinan di Jawa Tengah: Pendekatan Regresi Spasial. Jurnal Ilmiah Ekonomi, 13(2), 193–204. DOI: 10.22219/jiek.v13i2.14488

Yusuf, A. A., & Rabbani, A. F. (2018). The Impact of Social and Economic Factors on Poverty in Indonesia: A Spatial Econometrics Approach. Economics and Finance in Indonesia, 64(2), 133–152. DOI: 10.47291/efi.v64i2.607

Fitri, F., & Sitorus, S. R. (2019). Analisis Faktor Penentu Kemiskinan di Pulau Sumatera dengan Pendekatan Ekonometrika Spasial. Jurnal Pengelolaan Sumberdaya Alam dan Lingkungan (JPSL), 9(1), 58–66. DOI: 10.29244/jpsl.9.1.58-66

Dini, M. S., Rosadi, D., & Supartono. (2020). Spatial Durbin Model untuk Pemodelan Kemiskinan di Jawa Tengah. Jurnal Matematika, Statistika & Komputasi, 17(2), 185–196. DOI: 10.20956/jmsk.v17i2.10901

Siregar, R., & Saragih, H. M. (2021). Analisis Faktor-Faktor yang Mempengaruhi Kemiskinan di Indonesia dengan Pendekatan Spatial Error Model (SEM). Jurnal Pendidikan Tambusai, 5(2), 4867–4877. DOI: 10.31004/jptam.v5i2.1428

Suryaningrat, S. H. (2021). Aplikasi Akaike Information Criterion (AIC) pada Perhitungan Efisiensi Teknis Perikanan Pukat Cincin di Indonesia. Jurnal Penyuluhan Perikanan dan Kelautan, 15(2). Diakses dari: https://doi.org/10.29244/jmf.v11i2.38550.

Tursinawati, Y., Prawoto, N., & Sunarti, S. (2020). Model Spasial Kemiskinan di Provinsi Jawa Tengah: Pengaruh IPM dan Belanja Pemerintah. Jurnal Ekonomi Pembangunan, 21(2), 101–114. DOI: 10.14710/jep.21.2.101-114

LAMPIRAN

Laporan dalam Rpubs dapat diakses melalui: RPUBS_WEST_JAVA_POVERTY*

Dashboard dapat diakses melalui: DASHBOARD_WEST_JAVA_POVERTY*

Syntax R yang digunakan pada penelitian ini adalah sebagai berikut:

library(sf); library(dplyr); library(stringr)
library(spdep); library(spatialreg)
library(tmap); library(broom); library(ggplot2)
library(psych); library(readxl); library(viridis)   
library(scales); library(tibble); library(purrr)
library(forcats); library(htmltools)

# Import shapefile
setwd("C:\\Users\\ACER\\Documents\\SEMESTER 5\\SPATIAL\\UTS\\geoBoundaries-IDN-ADM2-all")
indo_adm2 <- st_read("geoBoundaries-IDN-ADM2.shp")
kab_jabar <- c("Bandung","Bandung Barat","Bekasi","Bogor","Ciamis","Cianjur","Cirebon",
               "Garut","Indramayu","Karawang","Kuningan","Majalengka","Pangandaran",
               "Purwakarta","Subang","Sukabumi","Sumedang","Tasikmalaya",
               "Kota Bandung","Kota Bekasi","Kota Bogor","Kota Cimahi",
               "Kota Cirebon","Kota Depok","Kota Sukabumi","Kota Tasikmalaya",
               "Kota Banjar")
kab_jabar
# normalisasi kolom shapeName
indo_adm2 <- indo_adm2 %>% mutate(shapeName_clean = str_to_lower(str_squish(shapeName)))
# normalisasi daftar kabupaten agar sama format
kab_jabar_clean <- str_to_lower(kab_jabar)
# subset Jawa Barat
jabar <- indo_adm2 %>% filter(shapeName_clean %in% kab_jabar_clean)
# jika masih kosong, coba pendekatan fuzzy match (string distance)
if(nrow(jabar) == 0){
  message("Tidak ada match exact; mencoba fuzzy join...")
  # minta user menyesuaikan nama; contoh pendekatan:
  # tampilkan beberapa nama dari shapefile untuk manual mapping
  print(unique(indo_adm2$shapeName)[1:80])
  stop("Periksa & adaptasi nama kab/kota lalu jalankan ulang.")
}

# Import atribut 
attr <- read_xlsx("C:/Users/ACER/Documents/SEMESTER 5/SPATIAL/UTS/DATA SPASIAL.xlsx")
# nama kolom kab/kota di CSV (sesuaikan jika berbeda)
# misal kolom bernama "KAB.KOTA" atau "KAB/KOTA" atau "KAB_KOTA"
names(attr)
# contoh: ubah nama menjadi 'kabupaten' dan normalisasi
attr <- attr %>%
  rename(kabupaten = starts_with("KAB")) %>%  # adapt jika perlu
  mutate(kab_clean = str_to_lower(str_squish(kabupaten)))
# cek beberapa baris
head(attr)

# Join atribut ke shapefile 
# join berdasarkan nama kabupaten (clean)
jabar <- jabar %>% left_join(attr, by = c("shapeName_clean" = "kab_clean"))
jabar
attach(jabar)

# Statistika Deskriptif ----
data_spasial <- jabar %>%
  sf::st_drop_geometry() %>%                         
  dplyr::select(KEMISKINAN, TPT, PDRB, IPM, GINI)
psych::describe(data_spasial)

# Histogram Setiap Variabel ----
# Histogram KEMISKINAN
ggplot(data_spasial, aes(x = KEMISKINAN)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  labs(title = "Distribusi Tingkat Kemiskinan", x = "Kemiskinan (%)", y = "Frekuensi")

# Histogram TPT
ggplot(data_spasial, aes(x = TPT)) +
  geom_histogram(bins = 10, fill = "lightgreen", color = "black") +
  labs(title = "Distribusi Tingkat Pengangguran Terbuka", x = "TPT (%)", y = "Frekuensi")

# Histogram PDRB
ggplot(data_spasial, aes(x = PDRB)) +
  geom_histogram(bins = 10, fill = "lightcoral", color = "black") +
  labs(title = "Distribusi PDRB per Kapita", x = "PDRB (juta rupiah)", y = "Frekuensi")

# Histogram IPM
ggplot(data_spasial, aes(x = IPM)) +
  geom_histogram(bins = 10, fill = "gold", color = "black") +
  labs(title = "Distribusi Indeks Pembangunan Manusia", x = "IPM", y = "Frekuensi")

# Histogram GINI
ggplot(data_spasial, aes(x = GINI)) +
  geom_histogram(bins = 10, fill = "plum", color = "black") +
  labs(title = "Distribusi Rasio Gini", x = "GINI", y = "Frekuensi")

# Boxplot Setiap Variabel ----
# Boxplot KEMISKINAN
ggplot(data_spasial, aes(y = KEMISKINAN)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  labs(title = "Boxplot Tingkat Kemiskinan", y = "Kemiskinan (%)", x = "") +
  theme_minimal()

# Boxplot TPT
ggplot(data_spasial, aes(y = TPT)) +
  geom_boxplot(fill = "lightgreen", color = "black") +
  labs(title = "Boxplot Tingkat Pengangguran Terbuka", y = "TPT (%)", x = "") +
  theme_minimal()

# Boxplot PDRB
ggplot(data_spasial, aes(y = PDRB)) +
  geom_boxplot(fill = "lightcoral", color = "black") +
  labs(title = "Boxplot PDRB per Kapita", y = "PDRB (juta rupiah)", x = "") +
  theme_minimal()

# Boxplot IPM
ggplot(data_spasial, aes(y = IPM)) +
  geom_boxplot(fill = "gold", color = "black") +
  labs(title = "Boxplot Indeks Pembangunan Manusia", y = "IPM", x = "") +
  theme_minimal()

# Boxplot GINI
ggplot(data_spasial, aes(y = GINI)) +
  geom_boxplot(fill = "plum", color = "black") +
  labs(title = "Boxplot Rasio Gini", y = "GINI", x = "") +
  theme_minimal()

# Plot Antar variabel ----
# Scatter Plot KEMISKINAN vs TPT
ggplot(data_spasial, aes(x = TPT, y = KEMISKINAN)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Kemiskinan dan TPT",
       x = "Tingkat Pengangguran Terbuka (%)",
       y = "Kemiskinan (%)") +
  theme_minimal()

# Scatter Plot KEMISKINAN vs PDRB
ggplot(data_spasial, aes(x = PDRB, y = KEMISKINAN)) +
  geom_point(color = "seagreen", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Kemiskinan dan PDRB per Kapita",
       x = "PDRB per Kapita (juta rupiah)",
       y = "Kemiskinan (%)") +
  theme_minimal()

# Scatter Plot KEMISKINAN vs IPM
ggplot(data_spasial, aes(x = IPM, y = KEMISKINAN)) +
  geom_point(color = "goldenrod", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Kemiskinan dan IPM",
       x = "Indeks Pembangunan Manusia",
       y = "Kemiskinan (%)") +
  theme_minimal()

# Scatter Plot KEMISKINAN vs GINI
ggplot(data_spasial, aes(x = GINI, y = KEMISKINAN)) +
  geom_point(color = "purple", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Kemiskinan dan Rasio Gini",
       x = "Rasio Gini",
       y = "Kemiskinan (%)") +
  theme_minimal()

# Choropleth ----
# Tema peta yang konsisten
theme_map <- function() {
  theme_minimal(base_size = 11) +
    theme(
      axis.title = element_blank(),
      axis.text  = element_blank(),
      panel.grid = element_blank(),
      legend.position = "right",
      plot.title = element_text(face = "bold")
    )
}

# 1) Kemiskinan
p_kemiskinan <- ggplot(jabar) +
  geom_sf(aes(fill = KEMISKINAN), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "Kemiskinan (%)", option = "C", direction = -1) +
  labs(title = "Penyebaran Kemiskinan (%) – Jawa Barat") +
  coord_sf() + theme_map()

# 2) TPT
p_tpt <- ggplot(jabar) +
  geom_sf(aes(fill = TPT), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "TPT (%)", option = "C", direction = -1) +
  labs(title = "Penyebaran TPT (%) – Jawa Barat") +
  coord_sf() + theme_map()

# 3) PDRB (biasanya miring kanan → gunakan skala log agar kontrasnya lebih informatif)
p_pdrb <- ggplot(jabar) +
  geom_sf(aes(fill = PDRB), color = "white", size = 0.2) +
  scale_fill_viridis_c(
    name = "PDRB per Kapita\n(juta rupiah, log10)",
    option = "C", direction = -1,
    trans = "log10",
    labels = label_number(accuracy = 0.1)
  ) +
  labs(title = "Penyebaran PDRB per Kapita – Jawa Barat") +
  coord_sf() + theme_map()

# 4) IPM
p_ipm <- ggplot(jabar) +
  geom_sf(aes(fill = IPM), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "IPM", option = "C", direction = -1) +
  labs(title = "Penyebaran IPM – Jawa Barat") +
  coord_sf() + theme_map()

# 5) GINI
p_gini <- ggplot(jabar) +
  geom_sf(aes(fill = GINI), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "Rasio Gini", option = "C", direction = -1) +
  labs(title = "Penyebaran Rasio Gini – Jawa Barat") +
  coord_sf() + theme_map()

# Tampilkan satu per satu (setiap baris akan memunculkan satu plot)
p_kemiskinan
p_tpt
p_pdrb
p_ipm
p_gini


# Autokorelasi Spatial
# Moran’s I ----
# 1) Weights Queen
nb_queen  <- poly2nb(jabar, queen = TRUE)
lw_queen  <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

# 2) Variabel yang diuji
vars <- c("KEMISKINAN", "TPT", "PDRB", "IPM", "Gini Ratio")

# 3) Helper: Moran's I (Asimtotik + Permutasi) aman dari NA
safe_moran_row_both <- function(y, varname, nb_full, nsim = 999, seed = 77) {
  # Tangani NA: subset data & struktur tetangga
  ok <- is.finite(y)
  y_use <- y[ok]
  nb_sub <- if (all(ok)) nb_full else spdep::subset.nb(nb_full, ok)
  lw_sub <- nb2listw(nb_sub, style = "W", zero.policy = TRUE)
  
  ## --- Asimtotik ---
  mt <- moran.test(y_use, lw_sub, zero.policy = TRUE)
  est <- mt$estimate
  
  I_asym  <- suppressWarnings(as.numeric(est[["Moran I statistic"]]))
  EI      <- suppressWarnings(as.numeric(est[["Expectation"]]))
  VarI    <- est[["Variance"]]
  if (is.null(VarI)) VarI <- est[["variance of Moran's I"]]
  VarI    <- suppressWarnings(as.numeric(VarI))
  Z_asym  <- suppressWarnings(as.numeric(mt$statistic))
  SD_asym <- if (!is.null(VarI) && is.finite(VarI)) sqrt(VarI) else NA_real_
  p_asym  <- mt$p.value
  
  ## --- Permutasi (Monte Carlo) ---
  set.seed(seed)
  mmc <- moran.mc(y_use, lw_sub, nsim = nsim, zero.policy = TRUE, alternative = "two.sided")
  I_perm <- suppressWarnings(as.numeric(mmc$statistic))
  p_perm <- mmc$p.value
  
  tibble(
    Variabel    = varname,
    Moran_I     = round(I_asym, 6),
    Expected_I  = round(EI, 6),
    SD          = round(SD_asym, 6),
    Z_value     = round(Z_asym, 6),
    P_asym      = signif(p_asym, 6),
    I_perm      = round(I_perm, 6),
    P_perm      = signif(p_perm, 6),
    nsim        = nsim
  )
}

# 4) Jalankan untuk semua variabel dan gabungkan ke satu tabel
vars <- names(jabar)[sapply(jabar, is.numeric)]

hasil_moran_both <- map_dfr(vars, ~ safe_moran_row_both(jabar[[.x]], .x, nb_queen))

# 5) Lihat tabel hasil
hasil_moran_both


# Geary C ----
safe_geary_row <- function(y, varname, nb_full) {
  # Tangani NA: subset data & struktur tetangga
  ok <- is.finite(y)
  if (!all(ok)) {
    # subset nb (fungsi subset.nb ada di spdep)
    nb_sub <- spdep::subset.nb(nb_full, ok)
  } else {
    nb_sub <- nb_full
  }
  lw_sub <- nb2listw(nb_sub, style = "W", zero.policy = TRUE)
  
  # Asimtotik 
  gt <- geary.test(y[ok], lw_sub, randomisation = TRUE, zero.policy = TRUE)
  est <- gt$estimate
  
  # Ambil nilai C & ekspektasi dengan fallback nama
  C_asym <- suppressWarnings(
    as.numeric(est[["Geary C statistic"]])
  )
  if (is.na(C_asym)) {
    C_asym <- suppressWarnings(as.numeric(est[["Geary's C"]]))
  }
  E_C <- suppressWarnings(as.numeric(est[["Expectation"]]))
  if (is.na(E_C)) E_C <- 1  # default teoritis Geary's C
  
  Z_asym <- suppressWarnings(as.numeric(gt$statistic))
  p_asym <- gt$p.value
  
  # Permutasi (Monte Carlo) 
  set.seed(77)
  gmc <- geary.mc(y[ok], lw_sub, nsim = 999, zero.policy = TRUE)
  C_perm <- suppressWarnings(as.numeric(gmc$statistic))
  p_perm <- gmc$p.value
  
  tibble(
    Variabel   = varname,
    Geary_C    = round(C_asym, 6),
    Expected_C = round(E_C, 6),
    Z_value    = round(Z_asym, 6),
    P_asym     = signif(p_asym, 6),
    C_perm     = round(C_perm, 6),
    P_perm     = signif(p_perm, 6)
  )
}

# 3) Hitung untuk semua variabel dan gabungkan ke tabel
hasil_geary <- map_dfr(vars, ~ safe_geary_row(jabar[[.x]], .x, nb_queen))

# 4) Lihat tabel hasil
hasil_geary

# Interpetasi
# Geary’s C < 1 → kecenderungan clustering (mirip Moran’s I positif).
# Geary’s C ≈ 1 → acak (tidak ada autokorelasi spasial).
# Geary’s C > 1 → dispersion (tetangga cenderung berbeda/kontras).
# Gunakan P_asym (uji asimtotik) dan P_perm (Monte Carlo) untuk signifikansi. Biasanya saya laporkan keduanya; jika hasilnya konsisten, kesimpulan lebih kuat.

# Local Moran’s I (LISA) -----
# 3️⃣ Fungsi untuk menghitung LISA & klasifikasi kuadran
safe_lisa <- function(y, varname, sf_obj, lw, sig_thr = 0.05) {
  ok <- is.finite(y)
  z <- as.numeric(scale(y[ok]))
  lag_z <- lag.listw(lw, z, zero.policy = TRUE)
  
  lisa_df <- as.data.frame(localmoran(z, lw, zero.policy = TRUE))
  
  # Hitung p-value dua sisi (karena versi spdep bisa beda nama)
  if ("Pr(z > 0)" %in% names(lisa_df)) {
    p_right <- lisa_df[,"Pr(z > 0)"]
    p_left  <- 1 - p_right
    p_two   <- 2*pmin(p_right, p_left)
  } else if ("Pr(z <= 0)" %in% names(lisa_df)) {
    p_left  <- lisa_df[,"Pr(z <= 0)"]
    p_right <- 1 - p_left
    p_two   <- 2*pmin(p_right, p_left)
  } else if ("Pr(z != E(Ii))" %in% names(lisa_df)) {
    p_two <- lisa_df[,"Pr(z != E(Ii))"]
  } else {
    zval <- lisa_df[,"Z.Ii"]
    p_two <- 2*pnorm(-abs(zval))
  }
  
  # Klasifikasi kuadran
  quad_raw <- ifelse(z >= 0 & lag_z >= 0, "High-High",
                     ifelse(z < 0 & lag_z < 0, "Low-Low",
                            ifelse(z >= 0 & lag_z < 0, "High-Low", "Low-High")))
  cluster <- ifelse(p_two <= sig_thr, quad_raw, "Not Significant")
  
  # Gabungkan ke sf
  sf_sub <- sf_obj[ok, ] %>%
    mutate(
      z_score  = z,
      lag_z    = lag_z,
      Ii       = lisa_df$Ii,
      Z_Ii     = lisa_df$Z.Ii,
      p_lisa   = p_two,
      LISA_cat = factor(cluster,
                        levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")),
      Variabel = varname
    )
  
  # Tabel ringkasan kategori
  tab <- sf_sub %>%
    st_drop_geometry() %>%
    group_by(LISA_cat) %>%
    summarise(
      Kabupaten = paste(kabupaten, collapse = ", "),
      Jumlah = n(),
      .groups = "drop"
    )
  
  list(sf = sf_sub, tabel = tab)
}

# 4️⃣ Jalankan untuk semua variabel
lisa_results <- map(vars, ~ safe_lisa(jabar[[.x]], .x, jabar, lw_queen))

# 5️⃣ Gabungkan tabel hasil
lisa_tables <- map2_df(lisa_results, vars, ~ mutate(.x$tabel, Variabel = .y))
lisa_tables

# 6️⃣ Plot peta cluster untuk tiap variabel
plots <- map2(lisa_results, vars, ~ {
  ggplot(.x$sf) +
    geom_sf(aes(fill = LISA_cat), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "High-High" = "#d7191c",
      "Low-Low"   = "#2c7bb6",
      "High-Low"  = "#fdae61",
      "Low-High"  = "#abd9e9",
      "Not Significant" = "grey80"
    ), drop = FALSE, name = "Cluster") +
    labs(title = paste0("Peta Cluster LISA - ", .y),
         subtitle = "Metode Local Moran’s I (Queen Contiguity)") +
    theme_minimal() +
    theme(axis.text = element_blank(),
          axis.title = element_blank(),
          panel.grid = element_blank())
})

# 7️⃣ Tampilkan tabel & peta contoh
print(lisa_tables)
plots[[1]]  # tampilkan peta KEMISKINAN
plots[[2]]  # tampilkan peta TPT
plots[[3]]  # tampilkan peta PDRB
plots[[4]]  # tampilkan peta IPM
plots[[5]]  # tampilkan peta GINI
# Jika ingin semua peta ditampilkan sekaligus:
# library(patchwork)
# wrap_plots(plots, ncol = 2)

# Getis–Ord 𝐺∗𝑖----
# === Palet kategori (sesuai contohmu) ===
cols_gstar <- c("Hot spot" = "#D7191C", "Cold spot" = "#1B9E77", "Not significant" = "grey99")

# === Helper: hitung Gi* satu variabel, hasilkan sf + tabel ringkasan ===
compute_gistar <- function(sf_obj, varname, lw) {
  x_raw <- sf_obj[[varname]]
  ok    <- is.finite(x_raw)
  
  # Subset bobot & data bila ada NA
  nb_sub <- if (all(ok)) lw$neighbours else spdep::subset.nb(lw$neighbours, ok)
  lw_sub <- nb2listw(nb_sub, style = "W", zero.policy = TRUE)
  
  # Matriks bobot (untuk G, G*)
  Wb      <- spdep::listw2mat(lw_sub)
  sum_x   <- sum(x_raw[ok])
  
  # --- G_i (tanpa i) ---
  num_G   <- as.numeric(Wb %*% x_raw[ok])
  den_G   <- (sum_x - x_raw[ok])
  G_raw   <- num_G / den_G
  
  # --- G_i* (dengan i): set w_ii = 1 ---
  Wb_star <- Wb; diag(Wb_star) <- 1
  num_Gs  <- as.numeric(Wb_star %*% x_raw[ok])
  den_Gs  <- sum_x
  G_star_raw <- num_Gs / den_Gs
  
  # --- Z-score Gi* (uji lokal) ---
  Gz <- spdep::localG(x_raw[ok], listw = lw_sub, zero.policy = TRUE)
  
  # Susun vektor full-length (NA di posisi non-ok)
  z_Gistar_full   <- rep(NA_real_, nrow(sf_obj)); z_Gistar_full[ok]   <- as.numeric(Gz)
  G_raw_full      <- rep(NA_real_, nrow(sf_obj)); G_raw_full[ok]      <- G_raw
  G_star_raw_full <- rep(NA_real_, nrow(sf_obj)); G_star_raw_full[ok] <- G_star_raw
  
  # Kategori
  GSTAR_cat <- case_when(
    z_Gistar_full >=  1.96 ~ "Hot spot",
    z_Gistar_full <= -1.96 ~ "Cold spot",
    TRUE                   ~ "Not significant"
  )
  
  # Gabungkan ke sf
  sf_out <- sf_obj |>
    mutate(
      !!paste0("G_raw_",   varname) := G_raw_full,
      !!paste0("Gstar_",   varname) := G_star_raw_full,
      !!paste0("z_Gi_",    varname) := z_Gistar_full,
      !!paste0("GSTAR_",   varname) := factor(GSTAR_cat,
                                              levels = c("Hot spot","Cold spot","Not significant")),
      Variabel = varname
    )
  
  # Tabel ringkasan hotspot/coldspot
  tab <- sf_out |>
    st_drop_geometry() |>
    transmute(
      Variabel = varname,
      Kabupaten = kabupaten,
      Kategori  = .data[[paste0("GSTAR_", varname)]],
      Gi_star   = round(.data[[paste0("z_Gi_", varname)]], 3)
    ) |>
    filter(Kategori %in% c("Hot spot","Cold spot")) |>
    arrange(Kategori, dplyr::desc(if_else(Kategori == "Hot spot", Gi_star, -Gi_star)))
  
  list(sf = sf_out, table = tab)
}

# === Jalankan untuk semua variabel ===
gistar_results <- map(vars, ~ compute_gistar(jabar, .x, lw_queen))

# === Tabel gabungan hotspot/coldspot semua variabel ===
gistar_tables <- map_dfr(gistar_results, "table")
gistar_tables
# -> berisi kolom: Variabel, Kabupaten, Kategori ("Hot spot"/"Cold spot"), Gi_star (z-score)

# === Plot peta cluster per variabel
# --- palet lebih kontras ---
cols_gstar <- c(
  "Hot spot"        = "#D7191C",
  "Cold spot"       = "#1B9E77",
  "Not significant" = "#EAEAEA"  # bukan putih; abu muda biar batas terlihat
)

# --- buat layer batas (garis) yang tajam ---
sf_plot <- gistar_results[[which(vars == "KEMISKINAN")]]$sf
batas_prov <- sf::st_as_sf(sf::st_cast(sf::st_union(sf_plot), "MULTILINESTRING"))

fill_col <- paste0("GSTAR_", "KEMISKINAN")
sf_plot <- sf_plot |>
  mutate(
    !!fill_col := as.character(.data[[fill_col]]),
    !!fill_col := ifelse(is.na(.data[[fill_col]]), "Not significant", .data[[fill_col]]),
    !!fill_col := forcats::fct(.data[[fill_col]],
                               levels = c("Hot spot","Cold spot","Not significant"))
  )

gistar_plots <- purrr::map2(gistar_results, vars, ~{
  sf_plot <- .x$sf
  fill_col <- paste0("GSTAR_", .y)
  
  batas_prov <- sf::st_as_sf(sf::st_cast(sf::st_union(sf_plot), "MULTILINESTRING"))
  
  sf_plot <- sf_plot |>
    mutate(
      !!fill_col := as.character(.data[[fill_col]]),
      !!fill_col := ifelse(is.na(.data[[fill_col]]), "Not significant", .data[[fill_col]]),
      !!fill_col := forcats::fct(.data[[fill_col]],
                                 levels = c("Hot spot","Cold spot","Not significant"))
    )
  
  ggplot(sf_plot) +
    geom_sf(aes(fill = .data[[fill_col]]), color = "#FFFFFF", size = 0.3) +
    geom_sf(data = sf::st_boundary(sf_plot), color = "#6B7280", size = 0.35, fill = NA) +
    geom_sf(data = batas_prov, color = "#374151", size = 0.6, fill = NA) +
    scale_fill_manual(values = cols_gstar,
                      limits = c("Hot spot","Cold spot","Not significant"),
                      drop = FALSE, name = "Kategori") +
    labs(title = paste0("Getis–Ord G* – ", .y)) +
    coord_sf(expand = FALSE) +
    theme_minimal(base_size = 12) +
    theme(axis.title = element_blank(),
          axis.text  = element_blank(),
          panel.grid = element_blank())
})
gistar_plots[[1]]  # KEMISKINAN
gistar_plots[[2]]  # TPT
gistar_plots[[3]]  # PDRB
gistar_plots[[4]]  # IPM
gistar_plots[[5]]  # GINI

# UJI OLS ----
ols_jabar <- lm(KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
summary(ols_jabar)

# Uji Moran’s I residual OLS ----
moran.test(residuals(ols_jabar), lw_queen)
moran.mc(residuals(ols_jabar), lw_queen, nsim = 999)

# Uji Lagrange Multiplier (LM) ----
lm.LMtests(ols_jabar, lw_queen, test = "all")

# 1. SDM ----
# Estimasi SDM (Durbin = TRUE mengaktifkan WX)
sdm_jabar <- lagsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  Durbin = TRUE,             # aktifkan Spatial Lag of X
  method = "eigen",          # metode umum untuk estimasi SDM
  zero.policy = TRUE
)
# Tampilkan ringkasan hasil
summary(sdm_jabar)
# Efek langsung, tidak langsung, dan total
impacts_sdm <- impacts(sdm_jabar, listw = lw_queen, R = 999)
summary(impacts_sdm, zstats = TRUE, short = FALSE)

# 2. GNS -----
# Estimasi model GNS
gns_jabar <- sacsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  Durbin = TRUE,          # aktifkan Spatial Lag of X
  method = "eigen",       # metode numerik untuk dekomposisi spasial
  zero.policy = TRUE
)

# Ringkasan hasil estimasi
summary(gns_jabar)
impacts_gns <- impacts(gns_jabar, listw = lw_queen, R = 999)
summary(impacts_gns, zstats = TRUE, short = TRUE)

# 3. SAC/SARAR ----
# Estimasi model SAC / SARAR
sac_jabar <- sacsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  method = "eigen",       # metode numerik dekomposisi eigen
  zero.policy = TRUE
)

# Ringkasan hasil model
summary(sac_jabar)

# 4. SAR ----
# Estimasi model SAR
sar_jabar <- lagsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  method = "eigen",
  zero.policy = TRUE
)

# Ringkasan hasil estimasi
summary(sar_jabar)
# Hitung impacts (dengan simulasi Monte Carlo R=500)
imp_sar <- impacts(
  sar_jabar,
  listw = lw_queen,
  R = 500,             # jumlah replikasi simulasi (biasanya 500–1000)
  zero.policy = TRUE
)

# Tampilkan hasil ringkas
summary(imp_sar, zstats = TRUE, short = TRUE)

# SLX ----

# Estimasi model SLX
slx_jabar <- lmSLX(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  Durbin = TRUE,        # aktifkan lag pada variabel X (WX)
  zero.policy = TRUE
)

# Ringkasan hasil model
summary(slx_jabar)
# Efek langsung, tidak langsung, dan total
impacts_slx <- impacts(slx_jabar, listw = lw_queen, R = 999)
summary(impacts_slx, zstats = TRUE, short = FALSE)

# SEM ----
# Estimasi model SEM
sem_jabar <- errorsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  method = "eigen",
  zero.policy = TRUE
)

# Ringkasan hasil model
summary(sem_jabar)

# SDEM ----
# Estimasi model SDEM
sdem_jabar <- errorsarlm(
  KEMISKINAN ~ TPT + PDRB + IPM + GINI,
  data = jabar,
  listw = lw_queen,
  Durbin = TRUE,         # aktifkan lag pada variabel X
  method = "eigen",
  zero.policy = TRUE
)

# Ringkasan hasil model
summary(sdem_jabar)

# BANDINGIN MODEL ----
compare_stats <- data.frame(
  Model = c("OLS","SLX","SAR","SDM","GNS","SAC/SARAR","SEM","SDEM"),
  LogLik = c(logLik(ols_jabar), logLik(slx_jabar), logLik(sar_jabar),
             logLik(sdm_jabar), logLik(gns_jabar), logLik(sac_jabar),
             logLik(sem_jabar), logLik(sdem_jabar)),
  AIC = c(AIC(ols_jabar), AIC(slx_jabar), AIC(sar_jabar),
          AIC(sdm_jabar), AIC(gns_jabar), AIC(sac_jabar),
          AIC(sem_jabar), AIC(sdem_jabar))
)
compare_stats

# VISUALISASI MODEL ----
# Prediksi in-sample & residual
jabar_pred <- jabar %>%
  mutate(
    y_hat   = as.numeric(fitted(sar_jabar)),       # prediksi SAR
    resid   = as.numeric(residuals(sar_jabar)),    # residual
    z_resid = as.numeric(scale(resid))             # residual distandarkan (opsional)
  )
# Peta Prediksi (nilai taksiran kemiskinan)
p_pred <- ggplot(jabar_pred) +
  geom_sf(aes(fill = y_hat), color = "white", size = 0.2) +
  scale_fill_viridis_c(name = "Prediksi\nKemiskinan (%)",
                       option = "C", direction = -1) +
  labs(title = "Peta Prediksi (SAR) – Kemiskinan Jawa Barat") +
  theme_minimal(base_size = 12) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank())
p_pred
# Peta Kemiskinan Aktual
p_kemiskinan

# Peta Residual (standar, pusat di 0, divergen)
# Pilih limit simetris agar mudah baca deviasi +/- (opsional)
lim_resid <- max(abs(jabar_pred$z_resid), na.rm = TRUE)

p_resid <- ggplot(jabar_pred) +
  geom_sf(aes(fill = z_resid), color = "white", size = 0.2) +
  scale_fill_gradient2(
    name = "Residual (z)",
    low = "#2c7bb6", mid = "white", high = "#d7191c",
    midpoint = 0, limits = c(-lim_resid, lim_resid), oob = squish
  ) +
  labs(title = "Peta Residual Standar (SAR) – Kemiskinan Jawa Barat",
       subtitle = "Merah: terprediksi terlalu rendah (underfit) • Biru: terprediksi terlalu tinggi (overfit)") +
  theme_minimal(base_size = 12) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank())
p_resid

# Uji autokorelasi residual + peta klasifikasi residual
# Uji Moran pada residual SAR (harusnya TIDAK signifikan untuk model yang baik)
moran_res_sar <- moran.test(jabar_pred$resid, lw_queen, zero.policy = TRUE)
moran_res_sar

# Klasifikasi sederhana residual z (untuk legenda diskrit)
jabar_pred <- jabar_pred %>%
  mutate(res_cat = cut(z_resid,
                       breaks = c(-Inf, -1.96, -1, 1, 1.96, Inf),
                       labels = c("<= -1.96", "(-1.96,-1]", "(-1,1]", "(1,1.96]", ">= 1.96"),
                       include.lowest = TRUE))

p_resid_disc <- ggplot(jabar_pred) +
  geom_sf(aes(fill = res_cat), color = "white", size = 0.2) +
  scale_fill_brewer(name = "Residual (z)", palette = "RdBu", direction = -1, drop = FALSE) +
  labs(title = "Peta Residual Standar (SAR) – Klasifikasi Diskrit") +
  theme_minimal(base_size = 12) +
  theme(axis.title = element_blank(),
        axis.text  = element_blank(),
        panel.grid = element_blank())
p_resid_disc

# Peta interaktif (leaflet)
# jabar_pred: sf dengan kolom kabupaten, y_hat, KEMISKINAN, z_resid

# 1) siapkan VEKTOR label (satu elemen per baris)
labels <- sprintf(
  "<b>%s</b><br/>Prediksi: %.2f<br/>Aktual: %.2f<br/>Residual (z): %.2f",
  jabar_pred$kabupaten,
  jabar_pred$y_hat,
  jabar_pred$KEMISKINAN,
  jabar_pred$z_resid
) |> lapply(HTML)  # <- list of HTML snippets, panjang = nrow(jabar_pred)

# 2) palet
pal_pred <- colorNumeric("viridis", domain = jabar_pred$y_hat)

# 3) plot leaflet (label per-fitur)
leaflet(jabar_pred) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor   = ~pal_pred(y_hat),
    color       = "white", weight = 1, fillOpacity = 0.85,
    label       = labels,                         # <- pakai vektor 'labels', bukan ~HTML(...)
    labelOptions = labelOptions(direction = "auto"),
    highlightOptions = highlightOptions(weight = 2, bringToFront = TRUE)
  ) %>%
  addLegend("bottomright", pal = pal_pred, values = ~y_hat, title = "Prediksi (%)")

# Uji Bobot W ----
# --- 1. Definisikan beberapa matriks bobot ---
nb_queen <- poly2nb(jabar, queen = TRUE)
lw_queen <- nb2listw(nb_queen, style = "W")

nb_rook  <- poly2nb(jabar, queen = FALSE)
lw_rook  <- nb2listw(nb_rook, style = "W")

coords <- sf::st_coordinates(sf::st_centroid(jabar))
nb_knn4 <- knn2nb(knearneigh(coords, k = 4))
lw_knn4 <- nb2listw(nb_knn4, style = "W")

# --- 2. Jalankan model SAR untuk setiap W ---
sar_queen <- lagsarlm(KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, listw = lw_queen)
sar_rook  <- lagsarlm(KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, listw = lw_rook)
sar_knn4  <- lagsarlm(KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, listw = lw_knn4)

# --- 3. Bandingkan hasil utama ---
data.frame(
  Bobot = c("Queen", "Rook", "KNN-4"),
  Rho   = c(sar_queen$rho, sar_rook$rho, sar_knn4$rho),
  AIC   = c(AIC(sar_queen), AIC(sar_rook), AIC(sar_knn4)),
  LogLik = c(logLik(sar_queen), logLik(sar_rook), logLik(sar_knn4))
)


# Uji Asumsi Klasik
model_klasik <- lagsarlm(KEMISKINAN ~ TPT + PDRB + IPM + GINI)
summary(model_klasik)

library(car)
uji_multikol <- vif(model_klasik)
uji_multikol

# UJI NORMALITAS RESIDUAL
shapiro.test(residuals(model_klasik))

# Plot visual (opsional tapi disarankan)
par(mfrow=c(1,2))
hist(residuals(model_klasik), main="Histogram Residual", xlab="Residuals")
qqnorm(residuals(model_klasik)); qqline(residuals(model_klasik))

# UJI HETEROSKEDASTISITAS
library(lmtest)
bptest(model_klasik)