BAB I PENDAHULUAN
1.1 Latar Belakang Masalah
1.2 Identifikasi Masalah
1.3. Maksud dan Tujuan Penelitian
1.4 Batasan Penelitian
BAB II TINJAUAN PUSTAKA
2.1 Spatial Dependence
2.2 Autokorelasi Spasial
2.3 Ordinary Least Squares
2.4 Uji Lagrange Multiplier
2.5 Spatial Autoregressive Model
2.6 Akaike Information Criterion
BAB III METODOLOGI PENELITIAN
3.1 Data
3.2 Metode Analisis
BAB IV HASIL DAN PEMBAHASAN
4.1 Analisis Deskriptif
4.2 Ordinary Least Square untuk Model Klasik
4.3 Uji Moran’s I Residual OLS
4.4 Autokorelasi Spasial
- 4.4.1 Moran’s I
- 4.4.2 Geary’s C
- 4.4.3 Local Moran’s I (LISA)
- 4.4.4 Getis–Ord 𝐺∗𝑖
4.5 Uji Lagrange Multiplier
4.6 Perbandingan Model
4.7 Estimasi Model SAR
4.8 Uji Asumsi Model SAR
4.9 Uji Sensitivitas Bobot
4.10 Interpretasi Hasil
BAB V KESIMPULAN
5.1 Kesimpulan
5.2 Saran
DAFTAR PUSTAKA
LAMPIRAN
Script Analisis
Script Dashboard

Link Dashboard: https://zeanova.shinyapps.io/WestJavaPovertyDashboard/

BAB I PENDAHULUAN

1.1 Latar Belakang Masalah

Kemiskinan mencerminkan rendahnya kemampuan individu dalam memenuhi standar hidup yang layak akibat keterbatasan pendapatan (Corral et al., 2020). Kondisi ini, apabila berlangsung secara berkelanjutan, dapat menimbulkan berbagai dampak sosial dan ekonomi seperti meningkatnya ketimpangan, terhambatnya pembangunan, munculnya konflik horizontal, hingga naiknya tingkat kriminalitas (UN, 2022). Jika tidak diatasi secara tepat, kemiskinan juga berpotensi menciptakan siklus antargenerasi yang sulit diputus.

Secara nasional, tingkat kemiskinan di Jawa Barat, menunjukkan tren penurunan hingga mencapai 7,27% pada tahun 2024. Namun, dalam kurun 2020–2024, penurunannya relatif lambat dan cenderung stagnan sebagaimana yang dapat dilihat pada Gambar 1. berikut.

Gambar 1. Persentase Penduduk Miskin di Jawa Barat 2002-2024

Kondisi kemiskinan ini diduga dipengaruhi oleh dinamika ekonomi dan sosial, seperti tingkat pengangguran terbuka (TPT), pertumbuhan ekonomi yang direpresentasikan melalui Produk Domestik Regional Bruto (PDRB), kualitas pembangunan manusia yang tercermin dalam Indeks Pembangunan Manusia (IPM), serta ketimpangan pendapatan yang diukur melalui Gini Ratio. Oleh karena itu, diperlukan analisis komprehensif untuk mengidentifikasi model terbaik yang mampu merepresentasikan hubungan antara variabel ekonomi dan sosial tersebut terhadap tingkat kemiskinan, sehingga kebijakan pengentasan kemiskinan di Jawa Barat dapat disusun secara lebih efektif, tepat sasaran, dan berkelanjutan.

1.2 Identifikasi Masalah

Pemetaan persentase kemiskinan di kabupaten/kota Provinsi Jawa Barat menunjukkan adanya indikasi ketidakterpisahan antar wilayah secara spasial, di mana tingkat kemiskinan pada suatu wilayah cenderung dipengaruhi oleh kondisi wilayah di sekitarnya. Indikasi tersebut dapat dilihat melalui Gambar 2. berikut ini.

Gambar 2. Peta Persentase Kemiskinan di Jawa Barat 2024

Hal ini mengindikasikan adanya spatial dependence atau ketergantungan spasial antar unit observasi. Dalam kondisi seperti ini, asumsi independensi antar observasi yang digunakan dalam model regresi klasik (OLS) tidak terpenuhi. Jika dependensi spasial tersebut diabaikan, maka hasil estimasi dapat menjadi bias dan tidak efisien.

Selain itu, tingkat kemiskinan diduga dipengaruhi oleh faktor ekonomi dan sosial, seperti Tingkat Pengangguran Terbuka (TPT), Produk Domestik Regional Bruto (PDRB), Indeks Pembangunan Manusia (IPM), dan Gini Ratio. Oleh karena itu, diperlukan pemodelan spasial yang mampu menangkap keterkaitan antar wilayah sekaligus mengukur pengaruh variabel-variabel tersebut secara akurat. Mengingat terdapat beberapa alternatif model spasial seperti SAR, SEM, SDM, SDEM, SLX, SAC, dan GNS, maka perlu dilakukan identifikasi model spasial terbaik yang paling sesuai untuk menggambarkan faktor ekonomi dan sosial yang memengaruhi tingkat kemiskinan di Jawa Barat.

1.3. Maksud dan Tujuan Penelitian

Maksud dari penelitian ini adalah menerapkan pendekatan spasial untuk menganalisis faktor yang mempengaruhi kemiskinan di Jawa Barat pada tahun 2024. 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 wilayah Provinsi Jawa Barat dengan unit analisis berupa 27 kabupaten/kota. Analisis difokuskan pada keterkaitan spasial tingkat kemiskinan antarwilayah sehingga aspek penyebab non-spasial tidak dibahas secara mendalam di luar variabel ekonomi dan sosial yang digunakan dalam model (TPT, PDRB, IPM, dan Gini Ratio). Selain itu, penelitian ini menggunakan data cross section tahun 2024 sehingga hasil yang diperoleh hanya mencerminkan kondisi pada periode tersebut dan tidak memperhatikan dinamika perubahan kemiskinan antarwaktu.

BAB II TINJAUAN PUSTAKA

2.1 Spatial Dependence

Ketergantungan spasial (spatial dependence) merujuk pada kondisi ketika suatu fenomena di suatu wilayah dipengaruhi oleh fenomena yang terjadi di wilayah lain yang berdekatan secara geografis. Menurut Anselin (1988), ketergantungan spasial muncul karena adanya the First Law of Geography, yang dikemukakan oleh Tobler (1970), yaitu: “Everything is related to everything else, but near things are more related than distant things.” Artinya, wilayah yang berdekatan cenderung memiliki karakteristik yang mirip karena adanya interaksi sosial-ekonomi, keterhubungan fisik, atau kedekatan budaya. Dalam konteks kemiskinan, misalnya, meningkatnya tingkat kemiskinan pada suatu wilayah dapat berkorelasi dengan wilayah sekitar karena adanya mobilitas tenaga kerja, transmisi ekonomi, dan kesamaan kondisi sosial.

2.2 Autokorelasi Spasial

Autokorelasi spasial adalah bentuk khusus dari ketergantungan spasial yang menunjukkan sejauh mana nilai suatu variabel pada satu wilayah berkorelasi dengan nilai variabel yang sama pada wilayah lain yang berdekatan. Dalam analisis spasial, autokorelasi spasial global dapat diukur menggunakan statistik Moran’s I yang dihitung dengan rumus sebagai berikut.

\[ I = \frac{n}{S_0} \cdot \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij} (x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^{n} (x_i - \bar{x})^2} \]

Nilai Moran’s I yang lebih besar dari nol (I > 0) mengindikasikan adanya autokorelasi spasial positif, yaitu wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi lainnya (high-high), demikian pula untuk nilai rendah (low-low). Sebaliknya, nilai Moran’s I yang kurang dari nol (I < 0) menunjukkan autokorelasi negatif atau pola checkerboard, di mana wilayah dengan nilai tinggi dikelilingi wilayah bernilai rendah atau sebaliknya. Jika Moran’s I mendekati nol (I ≈ 0), maka distribusi spasial diasumsikan acak dan tidak terdapat ketergantungan spasial yang signifikan.

Perhitungan Moran’s I membutuhkan matriks bobot spasial (W) yang merepresentasikan kedekatan antarwilayah. Matriks ini dapat disusun berdasarkan kontiguitas wilayah seperti rook contiguity (berbagi sisi), queen contiguity (berbagi sisi atau sudut), atau berdasarkan jarak geografis (distance-based weights). Selain Moran’s I, ukuran lain yang dapat digunakan untuk mengukur autokorelasi spasial meliputi Geary’s C untuk tingkat global, serta Local Moran’s I (LISA) dan Getis-Ord G* untuk mendeteksi pola klaster spasial pada tingkat lokal.

2.3 Ordinary Least Squares

Ordinary Least Square (OLS) adalah metode estimasi parameter untuk regresi klasik yang bertujuan untuk meminimalkan jumlah kuadrat galat. Model regresi linier klasik dapat dituliskan dalam bentuk umum sebagai:

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

di mana Y merupakan variabel dependen, X adalah matriks variabel independen, β adalah vektor parameter, dan ε menunjukkan komponen galat (error).

Metode ini membutuhkan pemenuhan pengujian asumsi klasik, diantaranya pengujian multikolinearitas, normalitas residual, dan homoskedastisitas. Metode ini mengasumsikan bahwa error antarwilayah bersifat independen dan homogen, sehingga tidak dapat menangkap spillover atau pengaruh suatu variabel dari satu daerah ke daerah lain. Hal ini menyebabkan hasil estimasi menjadi bias dan kurang akurat apabila diterapkan pada data dengan spatial dependence. Untuk mengatasi keterbatasan tersebut, digunakan model ekonometrika spasial yang mampu memperhitungkan spatial dependence antarwilayah.

2.4 Uji Lagrange Multiplier

Uji Lagrange Multiplier (LM) digunakan untuk mendeteksi keberadaan dependensi spasial dalam model regresi linier klasik (OLS). Uji ini membantu menentukan jenis keterkaitan spasial yang dominan, apakah terjadi pada bentuk spatial lag (ketergantungan pada variabel dependen wilayah tetangga) atau spatial error (ketergantungan pada komponen galat antarwilayah) (Anselin, 1988). Jika LM-Lag signifikan, maka terdapat indikasi bahwa model Spatial Autoregressive (SAR) lebih sesuai. Sebaliknya, apabila LM-Error signifikan, maka model Spatial Error Model (SEM) lebih direkomendasikan.

Uji LM terbagi menjadi dua jenis, yaitu LM biasa (LM-lag dan LM-error) dan Robust LM (RLM-lag dan RLM-error). Uji LM biasa hanya menguji masing-masing bentuk dependensi secara terpisah, sehingga dapat memberikan indikasi yang bias apabila terdapat dua jenis ketergantungan spasial secara simultan. Untuk mengatasi hal tersebut, digunakan Robust LM yang telah mengoreksi pengaruh jenis ketergantungan spasial lainnya, sehingga memberikan hasil yang lebih akurat dalam menentukan model spasial yang lebih dominan.

2.5 Spatial Autoregressive Model

Menurut Anselin (2003), Spatial Autoregressive Model (SAR) merupakan model regresi spasial yang mengakomodasi adanya pengaruh spasial pada variabel dependen. Model ini merupakan pengembangan dari regresi OLS dengan menambahkan komponen lag spasial dari variabel dependen sehingga nilai suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah-wilayah tetangganya. Bentuk umum model SAR adalah:

\[ Y = \rho W y + X \beta + \varepsilon \]

di mana Y adalah variabel dependen, X adalah matriks variabel independen, W merupakan matriks bobot spasial, ρ adalah koefisien autoregresif spasial yang mengukur kekuatan pengaruh wilayah tetangga, dan ε adalah error yang diasumsikan bersifat identik, independen, dan berdistribusi normal dengan varian konstan. Dalam model ini diasumsikan bahwa error tidak mengalami autokorelasi spasial, sehingga seluruh ketergantungan antarwilayah sudah tercermin melalui komponen lag spasial pada variabel dependen. Jika asumsi tersebut terpenuhi dan parameter ρ signifikan, maka model SAR dianggap tepat untuk digunakan.

2.6 Akaike Information Criterion

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 - 2 \ln(\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 III 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 merepresentasikan kondisi sosial ekonomi wilayah dalam satu periode waktu tertentu.

Variabel yang digunakan dalam penelitian ini terdiri atas satu variabel dependen, yaitu tingkat kemiskinan yang diukur melalui persentase penduduk miskin di setiap kabupaten/kota. Variabel independen yang digunakan meliputi: (1) Tingkat Pengangguran Terbuka (TPT), yang mencerminkan ketidaktersediaan lapangan kerja yang dapat menurunkan pendapatan rumah tangga; (2) Produk Domestik Regional Bruto (PDRB) per kapita sebagai indikator pertumbuhan dan kapasitas ekonomi wilayah; (3) Indeks Pembangunan Manusia (IPM) yang merefleksikan kualitas hidup melalui aspek pendidikan, kesehatan, dan ekonomi; serta (4) Gini Ratio yang mengukur tingkat ketimpangan pendapatan dan potensi kesenjangan sosial ekonomi. Pemilihan keempat variabel tersebut didasarkan pada teori bahwa kemiskinan dipengaruhi oleh kondisi pasar tenaga kerja, pertumbuhan ekonomi, kualitas sumber daya manusia, dan distribusi pendapatan.

Selain data atribut, penelitian ini juga memanfaatkan data spasial berupa peta batas administrasi wilayah kabupaten/kota di Jawa Barat yang diperoleh dari geoBoundaries Indonesia Administrative Level 2 (ADM2) dalam format shapefile (.shp) melalui situs www.geoboundaries.org. Data spasial ini digunakan untuk membentuk hubungan kedekatan antarwilayah yang direpresentasikan dalam matriks bobot spasial pada analisis ekonometrika spasial. Definisi operasional masing-masing variabel disajikan pada Tabel 1. berikut untuk mempermudah pemahaman dan sebagai acuan dalam pengolahan data.

Tabel 1. Definisi Operasional Variabel

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

3.2 Metode Analisis

Analisis dilakukan menggunakan pendekatan ekonometrika spasial dengan bantuan perangkat lunak R. Proses dimulai dari input data atribut dan data spasial yang kemudian dianalisis secara deskriptif untuk menggambarkan karakteristik masing-masing variabel. Selanjutnya dilakukan estimasi model OLS beserta uji asumsi klasik sebagai model dasar. Untuk mendeteksi keterkaitan spasial antarwilayah, dilakukan uji autokorelasi spasial global dan lokal, diikuti uji Lagrange Multiplier (LM) sebagai dasar penentuan model spasial yang sesuai.

Berdasarkan hasil uji LM, dilakukan estimasi beberapa model kandidat, seperti SAR, SDM, SAC, dan GNS. Pemilihan model terbaik dilakukan berdasarkan nilai AIC, signifikansi parameter spasial, serta tidak adanya autokorelasi residual. Model terpilih kemudian diuji kembali asumsi spasialnya dan dianalisis sensitivitas bobot spasialnya. Tahap akhir berupa interpretasi hasil yang meliputi direct effect, indirect effect (spillover), dan total effect. Alur selengkapnya disajikan pada Gambar 3.

Gambar 3. Alur Kerja Penelitian

BAB IV 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 2.

Tabel 2. Statistika Deskriptif

Variabel	Rata-rata	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

Adapun penyebaran dari setiap variabel dapat dilihat melalui Gambar 3. berikut.

Gambar 3. Peta Penyebaran Variabel

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 untuk Model Klasik

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 Model OLS

Jenis Uji	Statistik	Statistik Hitung	p-value	Kesimpulan
Linearitas (Ramsey RESET)	RESET	0.4760	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.56	0.6339	tidak terdapat gejala heteroskedastisitas

Adapun pengujian asumsi multikolinearitas dilakukan dengan melihat nilai VIF. Diketahui bahwa nilai VIF (X1) =1.67, VIF (X2) = 1.49, VIF (X3) = 2.17, VIF (X4) = 2.43. Semua nilai VIF < 5 yang berarti bahwa tidak terjadi multikolinearitas pada variabel independen. Diperoleh kesimpulan bahwa seluruh asumsi klasik pada model OLS terpenuhi.

4.3 Uji Moran’s I Residual OLS

Tabel 4. Hasil Uji Moran’s I

I	p (Monte Carlo)	Keputusan
2.3277	0.015	Tolak $H_{0}$

Berdasarkan hasil uji Moran’s I residual diperoleh kesimpulan bahwa residual pada model regresi OLS tidak bersifat acak, melainkan memiliki pola keterkaitan spasial yang signifikan. Dengan demikian, terdapat autokorelasi spasial pada error (residual), sehingga asumsi independensi residual dalam regresi klasik (OLS) tidak terpenuhi. Kondisi ini mengindikasikan bahwa model OLS tidak mampu menangkap pengaruh keterkaitan antarwilayah secara memadai, sehingga diperlukan pendekatan model ekonometrika spasial untuk menggambarkan hubungan spasial secara lebih akurat.

4.4 Autokorelasi Spasial

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

4.4.1 Moran’s I

Tabel 5. Hasil Uji Moran’s I Setiap 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 Ratio menunjukkan 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.4.2 Geary’s C

Tabel 6. Hasil Uji Geary’s C

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.4.3 Local Moran’s I (LISA)

Gambar 4. Hasil Klaster LISA

Hasil analisis Local Moran’s I (LISA) pada Gambar 4. menunjukkan bahwa variabel kemiskinan membentuk 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.4.4 Getis–Ord 𝐺∗𝑖

Gambar 5. Hasil Klaster Getis–Ord 𝐺∗𝑖

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.5 Uji Lagrange Multiplier

Tabel 7. Hasil Uji Lagrange Multiplier

Statistik Uji	Value	p-value	Keputusan
LM Error (RSer r)	3.4944	0.0616	Tidak Signifikan
LM Lag (RSlag)	8.8489	0.0029	Signifikan
Robust LM Error (AdjRSer r)	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.6 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. Hasil Pemodelan Model Spasial

Model	AIC	BIC	p uji ρ	p uji λ	p PDRB	p GINI	p TPT	p IPM
SAR	101.953	111	0.0004	-	0.294	0.097	0.510	2.67e-07
SDM	102.327	117	0.005	-	0.694*	0.473*	0.141*	0.776*
GNS	103.903	119	0.590	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

Catatan: * = p-value dari lag variabel

Berdasarkan hasil analisis pada Tabel 8., model terbaik yang dipilih adalah model SAR, yang ditunjukkan oleh nilai AIC dan BIC terendah, serta p-value uji Rho yang signifikan (0.0004), yang mengindikasikan adanya ketergantungan spasial yang kuat dalam data. Model SAR juga lebih efektif dalam menjelaskan pengaruh spasial dari Gini Ratio dan IPM terhadap tingkat kemiskinan di Jawa Barat, di mana kedua variabel tersebut memiliki p-value yang signifikan.

Model SDM dan GNS tidak dipilih karena semua variabel lag pada model ini tidak signifikan, yang menunjukkan bahwa efek spasial pada variabel independen (X) tidak memberikan kontribusi yang berarti dalam memprediksi kemiskinan di Jawa Barat. Selain itu, nilai Rho dan Lambda pada model GNS juga tidak signifikan. Ketidaksignifikanan kedua nilai ini menunjukkan bahwa model GNS tidak cukup efektif dalam menangkap ketergantungan spasial pada variabel dependen atau error, sehingga model ini kurang relevan untuk menggambarkan hubungan spasial dalam data kemiskinan di Jawa Barat. Model SAC/SARAR juga menunjukkan nilai Lambda yang tidak signifikan, berarti bahwa error dalam model tidak menunjukkan pola spasial yang signifikan, yang berarti tidak ada ketergantungan spasial pada error yang perlu diperhitungkan.

Oleh karena itu, model SAR menjadi pilihan terbaik karena mampu menjelaskan ketergantungan spasial yang ada dengan baik, sekaligus memenuhi prinsip kesederhanaan (parsimonious) dalam model, dengan mempertimbangkan nilai AIC dan BIC yang lebih rendah serta signifikansi yang kuat pada variabel dependen.

4.7 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.113(X_{1})_{i} - 0.0087(X_{2})_{i} - 0.413(X_{3})_{i}^{***} + 12.705(X_{4})_{i}^{*} \]

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

Tabel 9. Hasil Uji 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 Perhitungan Impact Model SAR

Variabel	Jenis Impact	Impact Measures	p-value	Keputusan
IPM	Direct	-0.481	0.018	Signifikan
	Indirect	-0.708	0.800	Tidak Signifikan
TPT	Direct	0.132	0.553	Tidak Signifikan
	Indirect	0.194	0.874	Tidak Signifikan
PDRB	Direct	-0.0101	0.351	Tidak Signifikan
	Indirect	-0.015	0.822	Tidak Signifikan
Gini Ratio	Direct	14.806	0.232	Tidak Signifikan
	Indirect	21.774	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.8 Uji Asumsi Model SAR

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

Tabel 11. Hasil Uji Asumsi Model SAR

Jenis Uji	Statistik	Statistik Hitung	p-value	Keputusan
Residual Autokorelasi (Moran’s I)	I	-0.111	0.7077	tidak terdapat autokorelasi spasial pada residual
Linearitas (Ramsey RESET)	RESET	0.021	0.979	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.0169	0.555	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.9 Uji Sensitivitas Bobot

Tabel 12. Hasil Uji Sensitivitas Bobot

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.10 Interpretasi Hasil

Hasil prediksi model SAR disajikan dalam gambar berikut.

Gambar 6. Peta Prediksi dan Residual Model SAR

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 V KESIMPULAN

5.1 Kesimpulan

Berdasarkan hasil analisis pemilihan model terbaik untuk kemiskinan di Jawa Barat 2024, disarankan agar pemerintah fokus pada peningkatan Indeks Pembangunan Manusia (IPM) sebagai langkah strategis untuk mengurangi tingkat kemiskinan. Peningkatan kualitas hidup melalui akses yang lebih baik terhadap pendidikan, kesehatan, dan infrastruktur dapat memberikan dampak positif dalam mengurangi kemiskinan. Selain itu, upaya untuk mengurangi ketimpangan pendapatan yang tercermin dari GINI Ratio juga sangat penting, karena semakin tinggi ketimpangan, semakin besar dampaknya terhadap peningkatan kemiskinan. Oleh karena itu, kebijakan yang mendorong pemerataan ekonomi, peningkatan kesejahteraan sosial, dan pengurangan ketimpangan pendapatan antar wilayah harus menjadi prioritas dalam perencanaan pembangunan. Pemerintah perlu menerapkan kebijakan yang lebih inklusif, memperkuat jaringan pengaman sosial, serta mendorong distribusi ekonomi yang lebih merata untuk menciptakan kemiskinan yang lebih rendah dan berkelanjutan di Jawa Barat.

5.2 Saran

Untuk penelitian selanjutnya, disarankan agar mempertimbangkan variabel-variabel lain yang dapat mempengaruhi tingkat kemiskinan di Jawa Barat, seperti faktor-faktor sosial, budaya, kesehatan, atau kebijakan pemerintah yang lebih spesifik. Penelitian juga dapat memperluas cakupan waktu dan wilayah dengan data panel untuk menganalisis dinamika jangka panjang dan variasi antar daerah. Akhirnya, penelitian lebih lanjut dapat mempertimbangkan pendekatan kualitatif untuk memahami faktor-faktor yang tidak terungkap oleh model kuantitatif, serta untuk mengeksplorasi kebijakan-kebijakan lokal yang mungkin mempengaruhi kemiskinan secara lebih spesifik.

DAFTAR PUSTAKA

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/

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 16 OKT 2024].

World Population Review. 2022. Poorest Countries in the World 2022. [Online] Available at: https:// worldpopulationreview.com/country-rankings /poorest-countries-in-the-world [Accessed 2 November 2022].

Pohan, H. L. M., & Vitale, J. D. (2016). Overcoming the poverty trap through education: An intergenerational study on Indonesia. Journal of Indonesian Applied Economics, 6(1), 1–21.

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.

LAMPIRAN

Script Analisis

Berikut script analisis R yang digunakan dalam analisis ini.

library(sf); library(dplyr); library(stringr)

## Warning: package 'sf' was built under R version 4.4.3

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spdep); library(spatialreg)

## Warning: package 'spdep' was built under R version 4.4.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.4.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Warning: package 'spatialreg' was built under R version 4.4.3

## Loading required package: Matrix

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(tmap); library(broom); library(ggplot2)

## Warning: package 'tmap' was built under R version 4.4.3

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

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

library(psych); library(readxl); library(viridis)

## Warning: package 'psych' was built under R version 4.4.3

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

## Warning: package 'readxl' was built under R version 4.4.3

## Warning: package 'viridis' was built under R version 4.4.3

## Loading required package: viridisLite

library(scales); library(tibble); library(purrr)

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

## 
## Attaching package: 'scales'

## The following object is masked from 'package:viridis':
## 
##     viridis_pal

## The following objects are masked from 'package:psych':
## 
##     alpha, rescale

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

## 
## Attaching package: 'purrr'

## The following object is masked from 'package:scales':
## 
##     discard

library(forcats); library(htmltools); library(lmtest); library(leaflet)

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

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.4.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Warning: package 'leaflet' was built under R version 4.4.3

# Import shapefile
setwd("D:/OneDrive - Universitas Padjadjaran/SMT 5/SPASIAL/PROJECT UTS/geoBoundaries-IDN-ADM2-all")
indo_adm2 <- st_read("geoBoundaries-IDN-ADM2.shp")

## Reading layer `geoBoundaries-IDN-ADM2' from data source 
##   `D:\OneDrive - Universitas Padjadjaran\SMT 5\SPASIAL\PROJECT UTS\geoBoundaries-IDN-ADM2-all\geoBoundaries-IDN-ADM2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 519 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84

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

##  [1] "Bandung"          "Bandung Barat"    "Bekasi"           "Bogor"           
##  [5] "Ciamis"           "Cianjur"          "Cirebon"          "Garut"           
##  [9] "Indramayu"        "Karawang"         "Kuningan"         "Majalengka"      
## [13] "Pangandaran"      "Purwakarta"       "Subang"           "Sukabumi"        
## [17] "Sumedang"         "Tasikmalaya"      "Kota Bandung"     "Kota Bekasi"     
## [21] "Kota Bogor"       "Kota Cimahi"      "Kota Cirebon"     "Kota Depok"      
## [25] "Kota Sukabumi"    "Kota Tasikmalaya" "Kota Banjar"

# 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("D:/OneDrive - Universitas Padjadjaran/SMT 5/SPASIAL/PROJECT 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)

## [1] "KAB/KOTA"   "KEMISKINAN" "IPM"        "GINI"       "RLS"       
## [6] "TPT"        "PDRB"       "PPP"        "AHH"

# 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)

## # A tibble: 6 × 10
##   kabupaten   KEMISKINAN   IPM  GINI   RLS   TPT  PDRB   PPP   AHH kab_clean  
##   <chr>            <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      
## 1 Bogor             7.05  73.6 0.398  8.39  7.34  54.9 11.6   72.2 bogor      
## 2 Sukabumi          6.87  70.2 0.342  7.34  7.11  31.2  9.82  72.1 sukabumi   
## 3 Cianjur          10.1   68.9 0.338  7.33  5.99  24.9  9.03  71.0 cianjur    
## 4 Bandung           6.19  74.6 0.364  9.15  6.36  44.1 11.4   74.6 bandung    
## 5 Garut             9.68  69.9 0.324  7.85  6.96  29.0  9.17  72.3 garut      
## 6 Tasikmalaya      10.2   70.0 0.363  7.97  3.74  25.9  8.96  70.4 tasikmalaya

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

## Simple feature collection with 27 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3703 ymin: -7.820979 xmax: 108.8469 ymax: -5.918148
## Geodetic CRS:  WGS 84
## First 10 features:
##        shapeName shapeISO                 shapeID shapeGroup shapeType
## 1        Bandung     <NA> 22746128B73297609500416        IDN      ADM2
## 2  Bandung Barat     <NA> 22746128B79756617117888        IDN      ADM2
## 3         Bekasi     <NA> 22746128B74302729610221        IDN      ADM2
## 4          Bogor     <NA> 22746128B99423643145371        IDN      ADM2
## 5         Ciamis     <NA> 22746128B17121530580660        IDN      ADM2
## 6        Cianjur     <NA> 22746128B34542084533901        IDN      ADM2
## 7        Cirebon     <NA> 22746128B77544549346240        IDN      ADM2
## 8          Garut     <NA> 22746128B66421464971893        IDN      ADM2
## 9      Indramayu     <NA> 22746128B59462811669242        IDN      ADM2
## 10      Karawang     <NA> 22746128B53840275251202        IDN      ADM2
##    shapeName_clean     kabupaten KEMISKINAN   IPM  GINI  RLS  TPT      PDRB
## 1          bandung       Bandung       6.19 74.59 0.364 9.15 6.36  44.12583
## 2    bandung barat Bandung Barat      10.49 70.77 0.400 8.24 6.70  32.50007
## 3           bekasi        Bekasi       4.80 76.80 0.367 9.76 8.82 128.69290
## 4            bogor         Bogor       7.05 73.63 0.398 8.39 7.34  54.85684
## 5           ciamis        Ciamis       7.39 73.64 0.326 8.10 3.37  34.33605
## 6          cianjur       Cianjur      10.14 68.89 0.338 7.33 5.99  24.91437
## 7          cirebon       Cirebon      11.00 72.30 0.374 7.65 6.74  28.13950
## 8            garut         Garut       9.68 69.91 0.324 7.85 6.96  29.01236
## 9        indramayu     Indramayu      11.93 70.72 0.306 6.95 6.25  55.16555
## 10        karawang      Karawang       7.86 73.82 0.376 8.05 8.04 121.29850
##       PPP      AHH                       geometry
## 1  11.392 74.55279 MULTIPOLYGON (((107.75 -6.8...
## 2   9.583 73.35373 MULTIPOLYGON (((107.525 -6....
## 3  12.500 74.62647 MULTIPOLYGON (((107.0982 -5...
## 4  11.563 72.15454 MULTIPOLYGON (((106.994 -6....
## 5  10.085 72.92110 MULTIPOLYGON (((108.3857 -7...
## 6   9.026 70.95969 MULTIPOLYGON (((107.3118 -7...
## 7  11.529 73.06450 MULTIPOLYGON (((108.685 -6....
## 8   9.168 72.28017 MULTIPOLYGON (((107.7202 -7...
## 9  11.010 72.80831 MULTIPOLYGON (((108.539 -6....
## 10 12.942 73.17286 MULTIPOLYGON (((107.1717 -6...

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

##            vars  n  mean    sd median trimmed   mad   min    max  range  skew
## KEMISKINAN    1 27  8.01  2.65   8.41    8.10  2.79  2.34  11.93   9.59 -0.36
## TPT           2 27  6.52  1.71   6.73    6.67  0.99  1.58   8.97   7.39 -1.08
## PDRB          3 27 53.11 33.26  39.70   48.14 15.84 24.91 147.08 122.17  1.51
## IPM           4 27 74.68  4.33  73.82   74.36  4.42 68.89  83.75  14.86  0.73
## GINI          5 27  0.37  0.05   0.37    0.37  0.05  0.29   0.48   0.19  0.28
##            kurtosis   se
## KEMISKINAN    -0.97 0.51
## TPT            0.89 0.33
## PDRB           1.19 6.40
## IPM           -0.59 0.83
## GINI          -0.68 0.01

# 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()

## `geom_smooth()` using formula = 'y ~ x'

# 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()

## `geom_smooth()` using formula = 'y ~ x'

# 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()

## `geom_smooth()` using formula = 'y ~ x'

# 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()

## `geom_smooth()` using formula = 'y ~ x'

# 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")

# 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
hasil_moran_both <- map_dfr(vars, ~ safe_moran_row_both(jabar[[.x]], .x, nb_queen))

# 5) Lihat tabel hasil
hasil_moran_both

## # A tibble: 5 × 9
##   Variabel   Moran_I Expected_I    SD Z_value    P_asym I_perm P_perm  nsim
##   <chr>        <dbl>      <dbl> <dbl>   <dbl>     <dbl>  <dbl>  <dbl> <dbl>
## 1 KEMISKINAN   0.478    -0.0385 0.141    3.65 0.000129   0.478  0       999
## 2 TPT          0.511    -0.0385 0.135    4.07 0.0000240  0.511  0       999
## 3 PDRB         0.211    -0.0385 0.134    1.86 0.0314     0.211  0.102   999
## 4 IPM          0.275    -0.0385 0.140    2.24 0.0127     0.275  0.032   999
## 5 GINI         0.299    -0.0385 0.140    2.40 0.00811    0.299  0.036   999

# 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

## # A tibble: 5 × 7
##   Variabel   Geary_C Expected_C Z_value    P_asym C_perm P_perm
##   <chr>        <dbl>      <dbl>   <dbl>     <dbl>  <dbl>  <dbl>
## 1 KEMISKINAN   0.410          1    3.76 0.0000865  0.410  0.001
## 2 TPT          0.465          1    2.92 0.00175    0.465  0.001
## 3 PDRB         0.660          1    1.82 0.0344     0.660  0.031
## 4 IPM          0.607          1    2.42 0.00787    0.607  0.005
## 5 GINI         0.504          1    3.07 0.00105    0.504  0.004

# 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

## # A tibble: 10 × 4
##    LISA_cat        Kabupaten                                     Jumlah Variabel
##    <fct>           <chr>                                          <int> <chr>   
##  1 High-High       Cirebon, Majalengka                                2 KEMISKI…
##  2 Low-Low         Bogor, Kota Bekasi                                 2 KEMISKI…
##  3 Not Significant Bandung, Bandung Barat, Bekasi, Ciamis, Cian…     23 KEMISKI…
##  4 Low-Low         Ciamis, Kota Tasikmalaya, Pangandaran, Tasik…      4 TPT     
##  5 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Cianj…     23 TPT     
##  6 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     27 PDRB    
##  7 Low-Low         Sumedang                                           1 IPM     
##  8 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     26 IPM     
##  9 Low-Low         Sumedang, Tasikmalaya                              2 GINI    
## 10 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     25 GINI

# 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)

## # A tibble: 10 × 4
##    LISA_cat        Kabupaten                                     Jumlah Variabel
##    <fct>           <chr>                                          <int> <chr>   
##  1 High-High       Cirebon, Majalengka                                2 KEMISKI…
##  2 Low-Low         Bogor, Kota Bekasi                                 2 KEMISKI…
##  3 Not Significant Bandung, Bandung Barat, Bekasi, Ciamis, Cian…     23 KEMISKI…
##  4 Low-Low         Ciamis, Kota Tasikmalaya, Pangandaran, Tasik…      4 TPT     
##  5 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Cianj…     23 TPT     
##  6 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     27 PDRB    
##  7 Low-Low         Sumedang                                           1 IPM     
##  8 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     26 IPM     
##  9 Low-Low         Sumedang, Tasikmalaya                              2 GINI    
## 10 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciami…     25 GINI

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

##      Variabel        Kabupaten  Kategori Gi_star
## 1  KEMISKINAN          Cirebon  Hot spot   2.489
## 2  KEMISKINAN       Majalengka  Hot spot   2.485
## 3  KEMISKINAN      Kota Bekasi Cold spot  -2.459
## 4  KEMISKINAN            Bogor Cold spot  -2.122
## 5         TPT      Pangandaran Cold spot  -3.278
## 6         TPT      Tasikmalaya Cold spot  -3.175
## 7         TPT           Ciamis Cold spot  -3.119
## 8         TPT Kota Tasikmalaya Cold spot  -2.513
## 9         IPM         Sumedang Cold spot  -2.051
## 10       GINI      Tasikmalaya Cold spot  -2.414
## 11       GINI         Sumedang Cold spot  -2.052

# -> 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)

## 
## Call:
## lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2564 -0.8408 -0.2461  0.8028  3.9446 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.770143   6.655302   6.727 9.23e-07 ***
## TPT         -0.087919   0.255841  -0.344    0.734    
## PDRB        -0.004582   0.012384  -0.370    0.715    
## IPM         -0.552077   0.114684  -4.814 8.28e-05 ***
## GINI        14.236635  11.360662   1.253    0.223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.72 on 22 degrees of freedom
## Multiple R-squared:  0.6431, Adjusted R-squared:  0.5782 
## F-statistic: 9.912 on 4 and 22 DF,  p-value: 9.655e-05

# Multikolinearitas → VIF
car::vif(ols_jabar) # Terpenuhi < 5

##      TPT     PDRB      IPM     GINI 
## 1.675131 1.491281 2.170410 2.438083

# Normalitas residual
shapiro.test(residuals(ols_jabar)) # Terpenuhi terima H0

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(ols_jabar)
## W = 0.96628, p-value = 0.5071

# Homoskedastisitas (Breusch–Pagan)
lmtest::bptest(ols_jabar) # Terpenuhi terima H0

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_jabar
## BP = 2.56, df = 4, p-value = 0.6339

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

## 
##  Moran I test under randomisation
## 
## data:  residuals(ols_jabar)  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 2.3277, p-value = 0.009964
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.28323118       -0.03846154        0.01909999

# Uji Permutasi
moran.mc(residuals(ols_jabar), lw_queen, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  residuals(ols_jabar) 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.28323, observed rank = 985, p-value = 0.015
## alternative hypothesis: greater

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

## Please update scripts to use lm.RStests in place of lm.LMtests

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## test weights: listw
## 
## RSerr = 3.4944, df = 1, p-value = 0.06158
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## test weights: listw
## 
## RSlag = 8.8489, df = 1, p-value = 0.002933
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## test weights: listw
## 
## adjRSerr = 1.0849, df = 1, p-value = 0.2976
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## test weights: listw
## 
## adjRSlag = 6.4394, df = 1, p-value = 0.01116
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar)
## test weights: listw
## 
## SARMA = 9.9338, df = 2, p-value = 0.006965

# 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)

## 
## Call:lagsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.089226 -0.655366  0.098144  0.608557  2.156242 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  39.1412149  13.4707343  2.9056  0.003665
## TPT           0.2570906   0.1778978  1.4452  0.148413
## PDRB         -0.0033225   0.0080008 -0.4153  0.677948
## IPM          -0.4249462   0.0895674 -4.7444 2.091e-06
## GINI         18.9755585   7.6787113  2.4712  0.013466
## lag.TPT      -0.4591547   0.3122474 -1.4705  0.141431
## lag.PDRB      0.0100997   0.0256549  0.3937  0.693819
## lag.IPM      -0.0686563   0.2408683 -0.2850  0.775616
## lag.GINI    -13.8560728  19.2998455 -0.7179  0.472796
## 
## Rho: 0.55505, LR test value: 8.029, p-value: 0.0046033
## Asymptotic standard error: 0.15968
##     z-value: 3.476, p-value: 0.00050893
## Wald statistic: 12.083, p-value: 0.00050893
## 
## Log likelihood: -40.16333 for mixed model
## ML residual variance (sigma squared): 1.048, (sigma: 1.0237)
## Number of observations: 27 
## Number of parameters estimated: 11 
## AIC: 102.33, (AIC for lm: 108.36)
## LM test for residual autocorrelation
## test value: 0.35755, p-value: 0.54987

# Efek langsung, tidak langsung, dan total
impacts_sdm <- impacts(sdm_jabar, listw = lw_queen, R = 999)
summary(impacts_sdm, zstats = TRUE, short = TRUE)

## Impact measures (mixed, exact):
##                  Direct    Indirect       Total
## TPT dy/dx   0.197308186 -0.65143378 -0.45412560
## PDRB dy/dx -0.001762873  0.01699438  0.01523151
## IPM dy/dx  -0.482473829 -0.62686447 -1.10933830
## GINI dy/dx 18.347667907 -6.84196900 11.50569891
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                  Direct     Indirect        Total
## TPT dy/dx   0.238978566   2.34519929   2.50393527
## PDRB dy/dx  0.009938453   0.08372882   0.09007449
## IPM dy/dx   0.099347883   0.90498164   0.95517030
## GINI dy/dx 11.851585543 140.80941724 149.66139389
## 
## Simulated z-values:
##                Direct    Indirect       Total
## TPT dy/dx   0.8272080 -0.26346185 -0.16781008
## PDRB dy/dx -0.1132518  0.29741656  0.26396809
## IPM dy/dx  -4.9762426 -0.85654032 -1.32911632
## GINI dy/dx  1.5791299 -0.05772811  0.07073656
## 
## Simulated p-values:
##            Direct    Indirect Total  
## TPT dy/dx  0.40812   0.79219  0.86673
## PDRB dy/dx 0.90983   0.76615  0.79180
## IPM dy/dx  6.483e-07 0.39170  0.18381
## GINI dy/dx 0.11431   0.95397  0.94361

# 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)

## 
## Call:sacsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01845 -0.73957  0.16932  0.66633  2.08399 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  50.8859986  29.1817765  1.7438   0.08120
## TPT           0.2518314   0.1839399  1.3691   0.17097
## PDRB         -0.0026465   0.0081610 -0.3243   0.74572
## IPM          -0.4399850   0.0956475 -4.6001 4.224e-06
## GINI         19.2987749   7.8683875  2.4527   0.01418
## lag.TPT      -0.5364431   0.3414370 -1.5711   0.11615
## lag.PDRB      0.0031261   0.0260318  0.1201   0.90441
## lag.IPM      -0.1942414   0.3390701 -0.5729   0.56674
## lag.GINI    -10.2259445  18.7534928 -0.5453   0.58556
## 
## Rho: 0.31514
## Asymptotic standard error: 0.58522
##     z-value: 0.5385, p-value: 0.59023
## Lambda: 0.38332
## Asymptotic standard error: 0.59239
##     z-value: 0.64707, p-value: 0.51758
## 
## LR test value: 20.478, p-value: 0.0022761
## 
## Log likelihood: -39.9514 for sacmixed model
## ML residual variance (sigma squared): 1.0572, (sigma: 1.0282)
## Number of observations: 27 
## Number of parameters estimated: 12 
## AIC: 103.9, (AIC for lm: 112.38)

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

## Impact measures (sacmixed, exact):
##                 Direct     Indirect         Total
## TPT dy/dx   0.21185842 -0.627435256 -0.4155768354
## PDRB dy/dx -0.00244606  0.003146329  0.0007002691
## IPM dy/dx  -0.46909794 -0.456970178 -0.9260681209
## GINI dy/dx 18.93636004 -5.688633168 13.2477268723
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                 Direct   Indirect       Total
## TPT dy/dx   0.34865335   6.640903   6.9247158
## PDRB dy/dx  0.02644265   0.579019   0.6036337
## IPM dy/dx   0.25034767   5.554890   5.7850484
## GINI dy/dx 15.14916279 281.032374 293.3327719
## 
## Simulated z-values:
##                 Direct    Indirect        Total
## TPT dy/dx   0.65242103 -0.15511006 -0.115903972
## PDRB dy/dx -0.04622883  0.01188410  0.009374409
## IPM dy/dx  -1.93974093 -0.13540182 -0.213957039
## GINI dy/dx  1.28430861 -0.02003272  0.047135396
## 
## Simulated p-values:
##            Direct   Indirect Total  
## TPT dy/dx  0.514130 0.87673  0.90773
## PDRB dy/dx 0.963128 0.99052  0.99252
## IPM dy/dx  0.052411 0.89229  0.83058
## GINI dy/dx 0.199034 0.98402  0.96241

# 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)

## 
## Call:sacsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.618260 -0.433617 -0.067514  0.402969  3.049638 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 27.3556558  5.3130193  5.1488 2.622e-07
## TPT          0.0828665  0.1562761  0.5303    0.5959
## PDRB        -0.0066943  0.0077519 -0.8636    0.3878
## IPM         -0.4018809  0.0762069 -5.2735 1.338e-07
## GINI        11.7918439  7.2013261  1.6375    0.1015
## 
## Rho: 0.73341
## Asymptotic standard error: 0.14298
##     z-value: 5.1295, p-value: 2.9057e-07
## Lambda: -0.31297
## Asymptotic standard error: 0.34317
##     z-value: -0.91202, p-value: 0.36176
## 
## LR test value: 13.483, p-value: 0.0011812
## 
## Log likelihood: -43.44892 for sac model
## ML residual variance (sigma squared): 1.194, (sigma: 1.0927)
## Number of observations: 27 
## Number of parameters estimated: 8 
## AIC: 102.9, (AIC for lm: 112.38)

# 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)

## 
## Call:lagsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.811768 -0.603120  0.058127  0.410091  3.318332 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 28.4991886  5.1628398  5.5201 3.389e-08
## TPT          0.1132664  0.1721913  0.6578   0.51067
## PDRB        -0.0087101  0.0083088 -1.0483   0.29450
## IPM         -0.4132629  0.0803241 -5.1449 2.676e-07
## GINI        12.7051666  7.6638468  1.6578   0.09736
## 
## Rho: 0.65269, LR test value: 12.428, p-value: 0.00042297
## Asymptotic standard error: 0.12453
##     z-value: 5.2411, p-value: 1.5964e-07
## Wald statistic: 27.469, p-value: 1.5964e-07
## 
## Log likelihood: -43.97625 for lag model
## ML residual variance (sigma squared): 1.3317, (sigma: 1.154)
## Number of observations: 27 
## Number of parameters estimated: 7 
## AIC: 101.95, (AIC for lm: 112.38)
## LM test for residual autocorrelation
## test value: 1.4493, p-value: 0.22865

# 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)

## Impact measures (lag, exact):
##                 Direct    Indirect       Total
## TPT dy/dx   0.13200312  0.19411954  0.32612265
## PDRB dy/dx -0.01015097 -0.01492769 -0.02507866
## IPM dy/dx  -0.48162546 -0.70826290 -1.18988836
## GINI dy/dx 14.80687435 21.77451312 36.58138746
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                 Direct    Indirect       Total
## TPT dy/dx  0.205683618  0.50647847  0.68261387
## PDRB dy/dx 0.009929494  0.03226359  0.03991287
## IPM dy/dx  0.103413999  0.83325294  0.89829356
## GINI dy/dx 9.366448647 38.55651486 45.32990789
## 
## Simulated z-values:
##               Direct   Indirect      Total
## TPT dy/dx   0.660727  0.4926766  0.5646395
## PDRB dy/dx -1.058324 -0.6203655 -0.7647617
## IPM dy/dx  -4.779253 -1.0765981 -1.5488481
## GINI dy/dx  1.630268  0.7295042  0.9573582
## 
## Simulated p-values:
##            Direct     Indirect Total  
## TPT dy/dx  0.50879    0.62224  0.57232
## PDRB dy/dx 0.28991    0.53502  0.44441
## IPM dy/dx  1.7595e-06 0.28166  0.12142
## GINI dy/dx 0.10304    0.46569  0.33839

# Koefisien Determinasi dari SAR
# Ambil variabel dependen (Y)
y <- jabar$KEMISKINAN
# Hitung pseudo-R² (versi LeSage & Pace)
R2_sar <- 1 - (sar_jabar$s2 / var(y))
R2_sar

## [1] 0.8101626

# 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)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##              Estimate   Std. Error  t value    Pr(>|t|) 
## (Intercept)  62.758787  16.696300    3.758844   0.001437
## TPT           0.188540   0.264058    0.714009   0.484378
## PDRB          0.001629   0.011827    0.137714   0.891995
## IPM          -0.494633   0.131822   -3.752286   0.001458
## GINI         19.022804  11.413348    1.666715   0.112876
## lag.TPT      -0.517804   0.459869   -1.125981   0.274963
## lag.PDRB     -0.002175   0.037836   -0.057478   0.954798
## lag.IPM      -0.264498   0.342100   -0.773158   0.449464
## lag.GINI     -9.208540  28.632934   -0.321607   0.751455

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

## Impact measures (SlX, glht, n-k):
##                  Direct     Indirect         Total
## TPT dy/dx   0.188539631 -0.517803935 -0.3292643039
## PDRB dy/dx  0.001628748 -0.002174735 -0.0005459867
## IPM dy/dx  -0.494633172 -0.264497733 -0.7591309044
## GINI dy/dx 19.022803508 -9.208539596  9.8142639114
## ========================================================
## Standard errors:
##                 Direct    Indirect       Total
## TPT dy/dx   0.26405783  0.45986930  0.50045586
## PDRB dy/dx  0.01182703  0.03783591  0.03690502
## IPM dy/dx   0.13182183  0.34210045  0.29496886
## GINI dy/dx 11.41334794 28.63293396 29.69946736
## ========================================================
## Z-values:
##                Direct    Indirect       Total
## TPT dy/dx   0.7140089 -1.12598064 -0.65792876
## PDRB dy/dx  0.1377140 -0.05747805 -0.01479437
## IPM dy/dx  -3.7522856 -0.77315810 -2.57359678
## GINI dy/dx  1.6667155 -0.32160657  0.33045252
## 
## p-values:
##            Direct     Indirect Total   
## TPT dy/dx  0.47522169 0.26017  0.510584
## PDRB dy/dx 0.89046643 0.95416  0.988196
## IPM dy/dx  0.00017523 0.43943  0.010065
## GINI dy/dx 0.09557100 0.74775  0.741058

# 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)

## 
## Call:errorsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.16279 -0.86976 -0.11366  0.80556  3.20803 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 33.649799   5.679634  5.9246 3.130e-09
## TPT          0.303626   0.192064  1.5809   0.11391
## PDRB        -0.005517   0.008960 -0.6157   0.53807
## IPM         -0.456813   0.092540 -4.9364 7.959e-07
## GINI        19.147958   8.015995  2.3887   0.01691
## 
## Lambda: 0.76303, LR test value: 8.3202, p-value: 0.0039206
## Asymptotic standard error: 0.1119
##     z-value: 6.8192, p-value: 9.1565e-12
## Wald statistic: 46.501, p-value: 9.1563e-12
## 
## Log likelihood: -46.03007 for error model
## ML residual variance (sigma squared): 1.4495, (sigma: 1.204)
## Number of observations: 27 
## Number of parameters estimated: 7 
## AIC: 106.06, (AIC for lm: 112.38)

# 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)

## 
## Call:errorsarlm(formula = KEMISKINAN ~ TPT + PDRB + IPM + GINI, data = jabar, 
##     listw = lw_queen, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01471 -0.68286  0.15237  0.55950  2.17018 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 64.6681636 12.4764414  5.1832 2.181e-07
## TPT          0.2309217  0.1778903  1.2981   0.19425
## PDRB        -0.0028092  0.0080739 -0.3479   0.72789
## IPM         -0.4642902  0.0819782 -5.6636 1.482e-08
## GINI        19.0523236  7.9231049  2.4047   0.01619
## lag.TPT     -0.5258451  0.3602129 -1.4598   0.14434
## lag.PDRB    -0.0021198  0.0238506 -0.0889   0.92918
## lag.IPM     -0.3421141  0.2055820 -1.6641   0.09609
## lag.GINI    -4.9824738 17.1255537 -0.2909   0.77110
## 
## Lambda: 0.6422, LR test value: 8.0242, p-value: 0.0046156
## Asymptotic standard error: 0.14849
##     z-value: 4.3249, p-value: 1.526e-05
## Wald statistic: 18.705, p-value: 1.526e-05
## 
## Log likelihood: -40.16573 for error model
## ML residual variance (sigma squared): 1.0095, (sigma: 1.0047)
## Number of observations: 27 
## Number of parameters estimated: 11 
## AIC: 102.33, (AIC for lm: 108.36)

# 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

##       Model    LogLik      AIC
## 1       OLS -50.19019 112.3804
## 2       SLX -44.17784 108.3557
## 3       SAR -43.97625 101.9525
## 4       SDM -40.16333 102.3267
## 5       GNS -39.95140 103.9028
## 6 SAC/SARAR -43.44892 102.8978
## 7       SEM -46.03007 106.0601
## 8      SDEM -40.16573 102.3315

# Menghitung BIC untuk setiap model spasial

# SAR Model
bic_sar <- BIC(sar_jabar)
# SDM Model
bic_sdm <- BIC(sdm_jabar)
# GNS Model
bic_gns <- BIC(gns_jabar)
# SAC Model
bic_sac <- BIC(sac_jabar)

# Menampilkan hasil BIC dari keempat model
bic_results <- tibble::tibble(
  Model = c("SAR", "SDM", "GNS", "SAC"),
  BIC = c(bic_sar, bic_sdm, bic_gns, bic_sac)
)

# Menampilkan hasil BIC
print(bic_results)

## # A tibble: 4 × 2
##   Model   BIC
##   <chr> <dbl>
## 1 SAR    111.
## 2 SDM    117.
## 3 GNS    119.
## 4 SAC    113.

# Prediksi nilai fitted
pred_sar <- sar_jabar$fitted.values
pred_sdm <- sdm_jabar$fitted.values

# Hitung MAPE masing-masing model
mape_sar <- mean(abs((jabar$KEMISKINAN - pred_sar) / jabar$KEMISKINAN)) * 100
mape_sdm <- mean(abs((jabar$KEMISKINAN - pred_sdm) / jabar$KEMISKINAN)) * 100

# Bandingkan hasil
cat("MAPE SAR:", mape_sar, "\n")

## MAPE SAR: 12.06973

cat("MAPE SDM:", mape_sdm, "\n")

## MAPE SDM: 11.40382

# 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)
  )

## This method assumes the response is known - see manual page

# 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

## 
##  Moran I test under randomisation
## 
## data:  jabar_pred$resid  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = -0.54681, p-value = 0.7077
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.11174814       -0.03846154        0.01796299

# 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))

## Warning: st_centroid assumes attributes are constant over geometries

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))
)

##   Bobot       Rho      AIC    LogLik
## 1 Queen 0.6526877 101.9525 -43.97625
## 2  Rook 0.6526877 101.9525 -43.97625
## 3 KNN-4 0.4871373 105.5565 -45.77824

# UJI ASUMSI SAR -----
# Residual SAR tidak ber-autokorelasi spasial
res_sar <- residuals(sar_jabar)
moran.test(res_sar, lw_queen, zero.policy = TRUE)  # terpenuhi p > 0.05

## 
##  Moran I test under randomisation
## 
## data:  res_sar  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = -0.54681, p-value = 0.7077
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.11174814       -0.03846154        0.01796299

## Normalitas residual SAR
shapiro.test(res_sar)  # p > 0.05 ~ normal

## 
##  Shapiro-Wilk normality test
## 
## data:  res_sar
## W = 0.95524, p-value = 0.2865

## Homoskedastisitas pada SAR
library(lmtest)
# residual SAR
res_sar <- residuals(sar_jabar)
# matriks X dari objek sarlm (kolom pertama = intercept)
x_vars <- sar_jabar$X[, -1, drop = FALSE]   # buang intercept
# siapkan data.frame untuk BP test
df_bp <- data.frame(res_sar = res_sar, as.data.frame(x_vars))
# Breusch–Pagan (uji homoskedastisitas)
bptest(res_sar ~ ., data = df_bp)

## 
##  studentized Breusch-Pagan test
## 
## data:  res_sar ~ .
## BP = 3.0169, df = 4, p-value = 0.555

# Secara Visual
plot(fitted(sar_jabar), residuals(sar_jabar),
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Uji Visual Homoskedastisitas (SAR)")

## This method assumes the response is known - see manual page

abline(h = 0, col = "red")

# Ekstrak residuals dan nilai yang diprediksi dari model SAR
sar_residuals <- residuals(sar_jabar)
sar_fitted <- fitted(sar_jabar)

## This method assumes the response is known - see manual page

# Tambahkan variabel kuadrat dan kubik dari prediksi ke dalam data
jabar$sar_fitted <- sar_fitted
jabar$sar_fitted_sq <- sar_fitted^2
jabar$sar_fitted_cub <- sar_fitted^3
# Melakukan regresi OLS untuk uji RESET
reset_model <- lm(KEMISKINAN ~ TPT + PDRB + IPM + GINI + sar_fitted_sq + sar_fitted_cub, data = jabar)

# Uji RESET
reset_test_sar <- resettest(reset_model, power = 2:3)

# Menampilkan hasil uji RESET
print(reset_test_sar)

## 
##  RESET test
## 
## data:  reset_model
## RESET = 0.021031, df1 = 2, df2 = 18, p-value = 0.9792

Script Dashboard

Berikut script dashboard yang digunakan dalam project ini:

atau untuk tampilan yang lebih baik dapat mengakses script pada tautan berikut: https://bit.ly/ScriptWestJavaPovertyDashboard

=========================================================

app.R — West Java Poverty Dashboard 2024 (VERSI 2)

(Navy theme, Inter font, Unpad logo, proper sorting with in-cell bars)

=========================================================

suppressPackageStartupMessages({ library(shiny) library(shinydashboard) library(DT) library(sf) library(dplyr) library(stringr) library(spdep) library(spatialreg) library(ggplot2) library(viridis) library(scales) library(psych) library(forcats) library(purrr) library(lmtest) library(car) library(readxl) library(leaflet) library(htmltools) library(gridExtra) library(openxlsx) library(showtext) library(sysfonts) })

————————-

GLOBAL: Fonts & Assets

————————-

if (!dir.exists(“www”)) dir.create(“www”) src_logo <- “D:/OneDrive - Universitas Padjadjaran/SMT 5/SPASIAL/PROJECT UTS/dashboard/www/logo-unpad.png” dst_logo <- “www/unpad_logo.png” if (file.exists(src_logo)) file.copy(src_logo, dst_logo, overwrite = TRUE)

Inter font (UI & plot)

try({ font_add_google(“Inter”, “Inter”, regular.wt = 400, bold.wt = 700) showtext_auto(enable = TRUE) }, silent = TRUE) theme_set(theme_minimal(base_size = 12, base_family = “Inter”))

————————-

GLOBAL: Data & helpers

————————-

shp_path <- “data/geoBoundaries-IDN-ADM2.shp” excel_path <- “data/DATA SPASIAL.xlsx”

indo_adm2 <- st_read(shp_path, quiet = TRUE) 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” ) indo_adm2 <- indo_adm2 %>% mutate(shapeName_clean = str_to_lower(str_squish(shapeName))) jabar <- indo_adm2 %>% filter(shapeName_clean %in% str_to_lower(kab_jabar))

attr <- read_xlsx(excel_path) %>% rename(kabupaten = dplyr::matches(“^KAB|KOTA”, ignore.case = TRUE)) %>% mutate(kabupaten = as.character(kabupaten), shapeName_clean = str_to_lower(str_squish(kabupaten))) jabar <- jabar %>% left_join(attr, by = “shapeName_clean”)

vars <- c(“KEMISKINAN”,“TPT”,“PDRB”,“IPM”,“GINI”) var_labels <- c( KEMISKINAN = “KEMISKINAN”, TPT = “TPT”, PDRB = “PDRB”, IPM = “IPM”, GINI = “GINI RATIO” )

Weights

nb_queen <- poly2nb(jabar, queen = TRUE) lw_queen <- nb2listw(nb_queen, style = “W”, zero.policy = TRUE) nb_rook <- poly2nb(jabar, queen = FALSE) lw_rook <- nb2listw(nb_rook, style = “W”, zero.policy = TRUE) coords <- sf::st_coordinates(sf::st_centroid(jabar)) nb_knn4 <- knn2nb(knearneigh(coords, k = 4)) lw_knn4 <- nb2listw(nb_knn4, style = “W”, zero.policy = TRUE)

round3_df <- function(df) df %>% mutate(across(where(is.numeric), ~ round(., 3)))

desc_table <- function(df_num) { d <- psych::describe(df_num) %>% as.data.frame() %>% tibble::rownames_to_column(“Variabel”) %>% select(Variabel, mean, sd, median, min, max, skew, kurtosis) d$var <- d$sd^2 d %>% select(Variabel, mean, sd, var, median, min, max, skew, kurtosis) %>% round3_df() }

safe_moran_both_raw <- function(y, nb, nsim = 999, seed = 77) { ok <- is.finite(y); y_use <- y[ok] nb_sub <- if (all(ok)) nb else spdep::subset.nb(nb, ok) lw_sub <- nb2listw(nb_sub, style = “W”, zero.policy = TRUE) mt <- moran.test(y_use, lw_sub, zero.policy = TRUE) set.seed(seed); mmc <- moran.mc(y_use, lw_sub, nsim = nsim, zero.policy = TRUE) list(mt = mt, mmc = mmc) } safe_geary_both_raw <- function(y, nb, nsim = 999, seed = 77) { ok <- is.finite(y); y_use <- y[ok] nb_sub <- if (all(ok)) nb else spdep::subset.nb(nb, ok) lw_sub <- nb2listw(nb_sub, style = “W”, zero.policy = TRUE) gt <- geary.test(y_use, lw_sub, randomisation = TRUE, zero.policy = TRUE) set.seed(seed); gmc <- geary.mc(y_use, lw_sub, nsim = nsim, zero.policy = TRUE) list(gt = gt, gmc = gmc) }

compute_lisa <- function(y, nb, sig_thr = 0.05) { ok <- is.finite(y); y_use <- y[ok] lw_sub <- nb2listw(if (all(ok)) nb else spdep::subset.nb(nb, ok), style = “W”, zero.policy = TRUE) z <- as.numeric(scale(y_use)) lz <- lag.listw(lw_sub, z, zero.policy = TRUE) lisa <- as.data.frame(localmoran(z, lw_sub, zero.policy = TRUE)) p_two <- if (“Pr(z > 0)” %in% names(lisa)) { p_right <- lisa[,“Pr(z > 0)”]; p_left <- 1 - p_right; 2pmin(p_right, p_left) } else if (“Pr(z <= 0)” %in% names(lisa)) { p_left <- lisa[,“Pr(z <= 0)”]; p_right <- 1 - p_left; 2pmin(p_right, p_left) } else if (“Pr(z != E(Ii))” %in% names(lisa)) { lisa[,“Pr(z != E(Ii))”] } else 2*pnorm(-abs(lisa[,“Z.Ii”])) quad <- ifelse(z >= 0 & lz >= 0, “High-High”, ifelse(z < 0 & lz < 0, “Low-Low”, ifelse(z >= 0 & lz < 0, “High-Low”, “Low-High”))) cl <- ifelse(p_two <= sig_thr, quad, “Not Significant”) list(z = z, lagz = lz, Ii = lisa$Ii, Zi = lisa$Z.Ii, p = p_two, cluster = cl, ok = ok) }

compute_gistar <- function(x, nb) { ok <- is.finite(x); x_use <- x[ok] nb_sub <- if (all(ok)) nb else spdep::subset.nb(nb, ok) lw_sub <- nb2listw(nb_sub, style = “W”, zero.policy = TRUE) Gz <- spdep::localG(x_use, listw = lw_sub, zero.policy = TRUE) z_full <- rep(NA_real_, length(x)); z_full[ok] <- as.numeric(Gz) cat <- case_when( z_full >= 1.96 ~ “Hot spot”, z_full <= -1.96 ~ “Cold spot”, TRUE ~ “Not significant” ) list(z = z_full, cat = factor(cat, levels = c(“Hot spot”,“Cold spot”,“Not significant”))) }

fit_models <- function(sfobj, yname, xnames, lw) { fml <- as.formula(paste(yname, “~”, paste(xnames, collapse = ” + “))) ols <- lm(fml, data = sfobj) slx <- lmSLX(fml, data = sfobj, listw = lw, Durbin = TRUE, zero.policy = TRUE) sar <- lagsarlm(fml, data = sfobj, listw = lw, method =”eigen”, zero.policy = TRUE) sdm <- lagsarlm(fml, data = sfobj, listw = lw, Durbin = TRUE, method = “eigen”, zero.policy = TRUE) gns <- sacsarlm(fml, data = sfobj, listw = lw, Durbin = TRUE, method = “eigen”, zero.policy = TRUE) sac <- sacsarlm(fml, data = sfobj, listw = lw, method = “eigen”, zero.policy = TRUE) sem <- errorsarlm(fml, data = sfobj, listw = lw, method = “eigen”, zero.policy = TRUE) sdem <- errorsarlm(fml, data = sfobj, listw = lw, Durbin = TRUE, method = “eigen”, zero.policy = TRUE) list(ols = ols, slx = slx, sar = sar, sdm = sdm, gns = gns, sac = sac, sem = sem, sdem = sdem) }

model_long_names <- c( “OLS” = “OLS (Ordinary Least Square)”, “SLX” = “SLX (Spatial Lag of X Model)”, “SAR” = “SAR (Spatial Autoregressive)”, “SDM” = “SDM (Spatial Durbin Model)”, “GNS” = “GNS (General Nesting Spatial)”, “SAC/SARAR” = “SAC/SARAR (Spatial Autoregressive Combined)”, “SEM” = “SEM (Spatial Error Model)”, “SDEM” = “SDEM (Spatial Durbin Error Model)” )

model_stats <- function(models, y) { get_loglik <- function(m) tryCatch(as.numeric(logLik(m)), error = function(e) NA_real_) get_aic <- function(m) tryCatch(AIC(m), error = function(e) NA_real_) get_r2 <- function(m, y) { if (inherits(m, “lm”)) return(summary(m)$r.squared) if (!is.null(m$s2)) return(1 - (m$s2 / var(y))) NA_real_ } tibble::tibble( ModelKey = c("OLS","SLX","SAR","SDM","GNS","SAC/SARAR","SEM","SDEM"), `Log-Likelihood` = c(get_loglik(models$ols), get_loglik(models$slx), get_loglik(models$sar), get_loglik(models$sdm), get_loglik(models$gns), get_loglik(models$sac), get_loglik(models$sem), get_loglik(models$sdem)), AIC = c(get_aic(models$ols), get_aic(models$slx), get_aic(models$sar), get_aic(models$sdm), get_aic(models$gns), get_aic(models$sac), get_aic(models$sem), get_aic(models$sdem)), `R²` = c(get_r2(models$ols, y), get_r2(models$slx, y), get_r2(models$sar, y), get_r2(models$sdm, y), get_r2(models$gns, y), get_r2(models$sac, y), get_r2(models$sem, y), get_r2(models$sdem, y)) ) %>% mutate(Model = unname(model_long_names[ModelKey])) %>% select(ModelKey, Model, Log-Likelihood, AIC, R²) %>% round3_df() }

bar_cell <- function(value, minv, maxv, color = “#1f2d5a”, txt_color = “#0d1b2a”) { if (is.na(value)) return(““) pct <- if (maxv == minv) 100 else scales::rescale(value, to = c(8, 100), from = c(minv, maxv)) lbl <- format(round(value, 3), nsmall = 3) htmltools::tags$div(style = sprintf(”background: linear-gradient(to right, %s %0.2f%%, rgba(0,0,0,0) %0.2f%%); padding:3px 8px; border-radius:6px; color:%s; font-weight:600; font-family:Inter,system-ui,sans-serif;“, color, pct, pct, txt_color ), lbl) |> as.character() }

=========================================================

UI

=========================================================

ui <- dashboardPage( dashboardHeader(title = “West Java Poverty Dashboard 2024”, tags$li(class = "dropdown", tags$a(href = “https://unpad.ac.id”, target = “_blank”, tags$img(src = "unpad_logo.png", height = "30px", style = "margin-top:10px;") ) ) ), dashboardSidebar(collapsed = FALSE, sidebarMenu(id = "tabs", menuItem("Data", tabName = "data", icon = icon("table")), menuItem("Statistika Deskriptif", tabName = "des", icon = icon("chart-bar")), menuItem("Peta", tabName = "peta", icon = icon("map")), menuItem("Autokorelasi Spasial", tabName = "auto", icon = icon("project-diagram")), menuItem("Pemodelan", tabName = "model", icon = icon("cogs")), menuItem("Tentang Kami", tabName = "about", icon = icon("info-circle")) ) ), dashboardBody( tags$head( tags$style(HTML(” @import url(‘https://fonts.googleapis.com/css2?family=Inter:wght@400;600;700&display=swap’);
/* Warna asli dipertahankan */
:root{
 --bg:#f5f7fb;
 --header:#0d1b2a;
 --sidebar:#14213d;
 --sidebar-hover:#1b2a49;
 --box-primary:#1f2d5a;
 --text:#0d1b2a;
 --text-inv:#ffffff;
 --ink:#2a2a2a;
}

/* Background dengan highlight lembut yang bergerak */
body, .content-wrapper, .right-side {
 background-color: var(--bg);
 font-family:'Inter', system-ui, -apple-system, Segoe UI, Roboto, Arial, sans-serif;
 position: relative;
 overflow-x: hidden;
}
body::before{
 content:'';
 position: fixed; inset:-20% -20% -20% -20%;
 background:
 radial-gradient(600px 400px at 15% 25%, rgba(13,27,42,0.06), transparent 60%),
 radial-gradient(700px 450px at 80% 20%, rgba(27,42,73,0.05), transparent 60%),
 radial-gradient(800px 500px at 30% 85%, rgba(31,45,90,0.05), transparent 60%);
 animation: floatBg 18s ease-in-out infinite alternate;
 pointer-events:none;
 z-index:-1;
}
@keyframes floatBg{ 
 0%{ transform: translate3d(0,0,0) scale(1); filter: blur(0px); }
 100%{ transform: translate3d(-2%, -1%, 0) scale(1.02); filter: blur(0.2px); }
}

/* Header & navbar dengan animasi halus */
.skin-blue .main-header .logo,
.skin-blue .main-header .navbar {
 background-color: var(--header);
 color: var(--text-inv);
 box-shadow: 0 6px 24px rgba(13,27,42,0.18);
 transition: box-shadow .25s ease, transform .2s ease, filter .25s ease;
}
.skin-blue .main-header .logo:hover{ transform: translateY(-1px); filter: brightness(1.04); }

/* Sidebar dengan transisi */
.skin-blue .main-sidebar {
 background-color: var(--sidebar);
}
.skin-blue .sidebar a { color: #dce6f2; transition: color .2s ease; }
.skin-blue .sidebar-menu>li>a{
 position: relative;
 transition: background-color .2s ease, color .2s ease, transform .15s ease;
}
.skin-blue .sidebar-menu>li.active>a, 
.skin-blue .sidebar-menu>li:hover>a {
 background: var(--sidebar-hover);
 color: #fff;
 transform: translateX(2px);
}
/* Indikator aktif animatif di sisi kiri item */
.skin-blue .sidebar-menu>li.active>a::before,
.skin-blue .sidebar-menu>li:hover>a::before{
 content:''; position:absolute; left:0; top:8px; bottom:8px; width:3px;
 background:#ffffff; border-radius:2px;
 animation: glow 2.4s ease-in-out infinite;
}
@keyframes glow{
 0%,100%{ box-shadow: 0 0 0 rgba(255,255,255,0.0); }
 50%{ box-shadow: 0 0 10px rgba(255,255,255,0.35); }
}

/* Kartu/box dengan efek hover */
.box {
 border-top: 3px solid var(--box-primary);
 border-radius: 10px;
 background: #fff;
 transition: transform .18s ease, box-shadow .25s ease, border-color .25s ease;
 box-shadow: 0 6px 20px rgba(0,0,0,0.06);
}
.box:hover{ transform: translateY(-2px); box-shadow: 0 14px 30px rgba(0,0,0,0.12); }
.box.box-solid.box-primary>.box-header { 
 color:#fff; background: var(--box-primary);
 position: relative; overflow: hidden;
}
/* sheen animasi lembut pada header box */
.box.box-solid.box-primary>.box-header::after{
 content:''; position:absolute; top:0; bottom:0; left:-40%;
 width:40%; background: linear-gradient(90deg, rgba(255,255,255,0.2), transparent);
 transform: skewX(-20deg);
 animation: sheen 4s linear infinite;
}
@keyframes sheen{ 
 0%{ left:-40%; } 
 100%{ left:120%; } 
}

/* Tombol dengan ripple & hover */
.btn, .btn-default, .btn-primary {
 border-radius: 10px;
 transition: transform .08s ease, box-shadow .2s ease, filter .2s ease;
 box-shadow: 0 6px 18px rgba(0,0,0,0.08);
}
.btn:hover{ filter: brightness(1.04); box-shadow: 0 10px 22px rgba(0,0,0,0.12); }
.btn:active{ transform: scale(0.98); }

/* Konten */
.content { padding-top: 10px; }

/* DataTable: highlight baris, transisi */
table.dataTable { color: var(--ink); }
.dataTables_wrapper .dataTables_paginate .paginate_button { 
 border-radius:6px;
 transition: background-color .2s ease, transform .1s ease;
}
.dataTables_wrapper .dataTables_paginate .paginate_button:hover{
 background: rgba(0,0,0,0.06); transform: translateY(-1px);
}
table.dataTable tbody tr{ transition: background-color .15s ease; }
table.dataTable tbody tr:hover{ background-color: rgba(31,45,90,0.06); }

/* Leaflet: kontrol & tooltip dengan efek halus */
.leaflet-container { font-family: 'Inter', sans-serif; }
.leaflet-bar a{
 transition: transform .12s ease, box-shadow .2s ease, filter .2s ease;
 box-shadow: 0 4px 12px rgba(0,0,0,0.08);
}
.leaflet-bar a:hover{ transform: translateY(-1px); filter: brightness(1.04); }
.leaflet-tooltip{
 background: rgba(255,255,255,0.95);
 border: 1px solid rgba(13,27,42,0.15);
 color: var(--ink);
 transition: opacity .2s ease;
}

/* Scrollbar halus */
::-webkit-scrollbar{ width: 10px; height:10px; }
::-webkit-scrollbar-thumb{
 background: rgba(13,27,42,0.2); border-radius: 8px;
 transition: background .2s ease;
}
::-webkit-scrollbar-thumb:hover{ background: rgba(13,27,42,0.32); }
“)), # Ripple kecil untuk tombol (tanpa ubah warna) tags$script(HTML(" (function(){ function ripple(e){ var btn = e.currentTarget; var circle = document.createElement('span'); circle.className = 'ripple-span'; var d = Math.max(btn.clientWidth, btn.clientHeight); circle.style.width = circle.style.height = d+'px'; var rect = btn.getBoundingClientRect(); circle.style.left = (e.clientX - rect.left - d/2)+'px'; circle.style.top = (e.clientY - rect.top - d/2)+'px'; btn.appendChild(circle); setTimeout(function(){ circle.remove(); }, 450); } document.addEventListener('DOMContentLoaded', function(){ document.querySelectorAll('.btn, .btn-default, .btn-primary').forEach(function(b){ b.style.position='relative'; b.style.overflow='hidden'; b.addEventListener('click', ripple); }); }); })(); ")), tags$style(HTML(” .ripple-span{ position:absolute; border-radius:50%; transform: scale(0); background: rgba(0,0,0,0.12); animation: ripple .45s ease-out; pointer-events:none; opacity:1; } @keyframes ripple{ to{ transform: scale(2.6); opacity:0; } } “)) ), tabItems( # ———— DATA ———— tabItem(tabName =”data”, fluidRow( box(width = 3, title = “Download Data”, status = “primary”, solidHeader = TRUE, tags$p(tags$em(“Sumber data: Badan Pusat Statistik Provinsi Jawa Barat 2024”)), downloadButton(“download_data_csv”, “Download CSV”), br(), br(), downloadButton(“download_data_xlsx”, “Download XLSX”) ), box(width = 9, title = “Data Overview”, status = “primary”, solidHeader = TRUE, DTOutput(“data_table”) ) ) ),
 # ------------ DESKRIPTIF ------------
 tabItem(tabName = "des",
 fluidRow(
 box(width = 3, title = "Pengaturan", status = "primary", solidHeader = TRUE,
 selectInput("var_desc", "Pilih variabel:", choices = setNames(vars, var_labels), selected = "KEMISKINAN"),
 sliderInput("bins", "Jumlah bins histogram:", min = 5, max = 30, value = 10),
 sliderInput("pointSize", "Ukuran titik scatter:", min = 1, max = 6, value = 3),
 selectInput("x_for_scatter", "Pilih variabel X untuk scatter plot:",
 choices = setNames(vars[vars != "KEMISKINAN"], var_labels[vars != "KEMISKINAN"]),
 selected = "IPM"),
 br(),
 downloadButton("download_plots", "Download Plot (PNG)")
 ),
 box(width = 9, title = "Tabel Statistika Deskriptif", status = "primary", solidHeader = TRUE,
 DTOutput("desc_table"),
 hr(),
 fluidRow(
 column(6, plotOutput("hist_plot", height = 340)),
 column(6, plotOutput("box_plot", height = 340))
 ),
 fluidRow(
 column(12, plotOutput("scatter_plot", height = 520))
 )
 )
 )
 ),
 
 # ------------ PETA (LEAFLET + download PNG statis) ------------
 tabItem(tabName = "peta",
 fluidRow(
 box(width = 3, title = "Pengaturan", status = "primary", solidHeader = TRUE,
 selectInput("var_map_left", "Peta Kiri: variabel", choices = setNames(vars, var_labels), selected = "KEMISKINAN"),
 selectInput("var_map_right", "Peta Kanan: variabel", choices = setNames(vars, var_labels), selected = "IPM")
 ),
 box(width = 9, title = "Peta Penyebaran Data", status = "primary", solidHeader = TRUE,
 fluidRow(
 column(6,
 leafletOutput("map_left", height = 520),
 br(), downloadButton("dl_map_left_png", "Download PNG (Peta Kiri)")
 ),
 column(6,
 leafletOutput("map_right", height = 520),
 br(), downloadButton("dl_map_right_png", "Download PNG (Peta Kanan)")
 )
 )
 )
 )
 ),
 
 # ------------ AUTOKORELASI SPASIAL ------------
 tabItem(tabName = "auto",
 fluidRow(
 box(width = 3, title = "Pengaturan", status = "primary", solidHeader = TRUE,
 selectInput("var_auto", "Pilih variabel:", choices = setNames(vars, var_labels), selected = "KEMISKINAN"),
 selectInput("auto_method", "Metode:",
 choices = c("Moran's I","Geary's C","LISA","Getis-Ord G*"),
 selected = "Moran's I"),
 conditionalPanel(
 "input.auto_method == 'LISA'",
 checkboxGroupInput("lisa_filter", "Filter kategori:",
 choices = c("High-High","Low-Low","High-Low","Low-High","Not Significant"),
 selected = c("High-High","Low-Low","High-Low","Low-High","Not Significant"))
 ),
 conditionalPanel(
 "input.auto_method == 'Getis-Ord G*'",
 checkboxGroupInput("gstar_filter", "Filter kategori:",
 choices = c("Hot spot","Cold spot","Not significant"),
 selected = c("Hot spot","Cold spot","Not significant"))
 ),
 br(),
 downloadButton("dl_auto_map_png", "Download PNG (Peta)")
 ),
 box(width = 9, title = "Hasil Analisis", status = "primary", solidHeader = TRUE,
 DTOutput("auto_table"),
 br(),
 leafletOutput("auto_map", height = 520)
 )
 )
 ),
 
 # ------------ PEMODELAN ------------
 tabItem(tabName = "model",
 fluidRow(
 box(width = 3, title = "Pengaturan", status = "primary", solidHeader = TRUE,
 selectInput("y_model", "Pilih Variabel Y:", choices = setNames(vars, var_labels), selected = "KEMISKINAN"),
 checkboxGroupInput("x_model", "Pilih Variabel X:",
 choices = setNames(vars, var_labels),
 selected = c("TPT","PDRB","IPM","GINI")),
 selectInput("w_choice", "Pilih jenis matriks bobot spasial:",
 choices = c("Queen","Rook","KNN-4"), selected = "Queen"),
 selectInput("model_pred", "Pilih model untuk peta residual:",
 choices = c("SAR","SDM","SAC/SARAR","GNS","SLX","SEM","SDEM"),
 selected = "SAR"),
 actionButton("fit_btn", "Fit / Refit Model")
 ),
 box(width = 9, title = "Ukuran Kebaikan Model", status = "primary", solidHeader = TRUE,
 DTOutput("model_table"),
 uiOutput("best_model_cards"),
 tags$small(HTML(
 "Catatan: R² pada OLS adalah koefisien determinasi konvensional dari regresi linear. 
 Untuk model spasial, nilai R² yang ditampilkan merupakan Pseudo-R² 
 yang dihitung sebagai 1 − (σ²_model / Var(y))"
 )),
 hr(),
 h4("Peta Residual"),
 leafletOutput("resid_leaf", height = 520),
 br(), downloadButton("dl_resid_map_png", "Download PNG (Residual)"),
 hr(),
 h4("Uji Sensitivitas Matriks Pembobot Spasial (berdasarkan pilihan Y dan X)"),
 DTOutput("sens_table"),
 plotOutput("sens_bar", height = 360),
 br(), downloadButton("dl_sens_bar", "Download PNG (Bar Chart Sensitivitas)")
 )
 )
 ),
 
 # ------------ TENTANG KAMI ------------
 tabItem(tabName = "about",
 fluidRow(
 box(width = 12, title = "Tentang Kami", status = "primary", solidHeader = TRUE,
 HTML("
 <h4>Judul Analisis</h4>
 PEMILIHAN MODEL TERBAIK DALAM MENGANALISIS FAKTOR EKONOMI DAN SOSIAL 
 TERHADAP TINGKAT KEMISKINAN DI JAWA BARAT TAHUN 2024
 <h4>Dibuat pada</h4>
 17 Oktober 2025
 <h4>Dibuat oleh</h4>
 Humaira Zeanova — 140610230030
 Dosen Pengampu: Dr. I Gede Nyoman Mindra Jaya, S. Si., M.Si.
 Instansi: Program Studi Statistika, FMIPA, Universitas Padjadjaran.
 Dashboard ini dirancang untuk membantu eksplorasi data kemiskinan tingkat kabupaten/kota di Jawa Barat,
 analisis autokorelasi spasial, pemodelan spasial lanjutan, serta evaluasi sensitivitas matriks bobot. 
 Pengguna dapat menyesuaikan variabel, model, serta parameter visualisasi secara interaktif.
 ")
 )
 )
 )
)
) )

=========================================================

SERVER

=========================================================

server <- function(input, output, session){

observeEvent(input$y_model, { if (input$y_model %in% input$x_model) { updateCheckboxGroupInput(session, "x_model", selected = setdiff(input$x_model, input$y_model) ) } }) observeEvent(input$x_model, { if (input$y_model %in% input$x_model) { updateCheckboxGroupInput(session, “x_model”, selected = setdiff(input$x_model, input$y_model) ) } })

# ———- DATA data_df <- reactive({ jabar |> st_drop_geometry() |> transmute( Kabupaten = kabupaten, KEMISKINAN (%) = KEMISKINAN, TPT (%) = TPT, PDRB (juta) = PDRB, IPM = IPM, GINI RATIO = GINI ) |> round3_df() }) output$data_table <- renderDT({ DT::datatable(data_df(), options = list(pageLength = 15), rownames = FALSE) }) output$download_data_csv <- downloadHandler( filename = function() “data_spasial_jabar.csv”, content = function(file) write.csv(data_df(), file, row.names = FALSE, na = ““) ) output$download_data_xlsx <- downloadHandler( filename = function()”data_spasial_jabar.xlsx”, content = function(file) openxlsx::write.xlsx(data_df(), file, overwrite = TRUE) )

# ———- DESKRIPTIF df_num <- reactive({ jabar |> sf::st_drop_geometry() |> dplyr::select(all_of(vars)) }) output$desc_table <- renderDT({ DT::datatable(desc_table(df_num()), rownames = FALSE, options = list(pageLength = 8)) })

# Palet plot tema col_hist <- “#1f2d5a” # navy col_box <- “#3a6ea5” # steel blue col_outline<- “#0d1b2a” # whisker/outline col_point <- “#3a6ea5” # scatter + boxplot jitter col_smooth <- “#f2a900” # amber

output$hist_plot <- renderPlot({ v <- req(input$var_desc) x <- jabar[[v]] brks <- pretty(range(x, na.rm = TRUE), n = input$bins) ggplot(jabar, aes(x = .data[[v]])) + geom_histogram(breaks = brks, color = "black", fill = col_hist, closed = "right") + scale_x_continuous(breaks = brks) + labs(title = paste("Histogram", var_labels[[v]]), x = paste0(var_labels[[v]]), y = "Frekuensi") + theme_minimal(base_size = 12, base_family = "Inter") }) output$box_plot <- renderPlot({ v <- req(input$var_desc) ggplot(jabar, aes(y = .data[[v]])) + geom_boxplot(fill = "#8B3A3A", color = "black") + labs(title = paste("Boxplot", var_labels[[v]]), y = var_labels[[v]], x = "") + theme_minimal(base_size = 12) }) output$scatter_plot <- renderPlot({ x <- req(input$x_for_scatter) ggplot(jabar |> st_drop_geometry(), aes(x = .data[[x]], y = KEMISKINAN)) + geom_point(size = input$pointSize, color = col_point) + geom_smooth(method = “lm”, se = FALSE, linetype = “dashed”, color = col_smooth) + labs(title = paste(“KEMISKINAN vs”, var_labels[[x]]), x = var_labels[[x]], y = “KEMISKINAN”) + theme_minimal(base_size = 12, base_family = “Inter”) }) output$download_plots <- downloadHandler( filename = function() "deskriptif_plots.png", content = function(file) { v <- isolate(input$var_desc); x <- isolate(input$x_for_scatter) xv <- jabar[[v]] brks <- pretty(range(xv, na.rm = TRUE), n = isolate(input$bins)) p1 <- ggplot(jabar, aes(x = .data[[v]])) + geom_histogram(breaks = brks, color = “black”, fill = col_hist, closed = “right”) + scale_x_continuous(breaks = brks) + labs(title = paste(“Histogram”, var_labels[[v]]), x = paste0(var_labels[[v]]), y = “Frekuensi”) + theme_minimal(base_size = 12, base_family = “Inter”) p2 <- ggplot(jabar, aes(y = .data[[v]])) + geom_boxplot(fill = “#8B3A3A”, color = “black”) + labs(title = paste(“Boxplot”, var_labels[[v]]), y = var_labels[[v]], x = ““) + theme_minimal(base_size = 12) p3 <- ggplot(jabar |> st_drop_geometry(), aes(x = .data[[x]], y = KEMISKINAN)) + geom_point(size = isolate(input$pointSize), color = col_point) + geom_smooth(method =”lm”, se = FALSE, linetype = “dashed”, color = col_smooth) + labs(title = paste(“KEMISKINAN vs”, var_labels[[x]]), x = var_labels[[x]], y = “KEMISKINAN”) + theme_minimal(base_size = 12, base_family = “Inter”) g <- gridExtra::arrangeGrob(p1, p2, p3, ncol = 2) ggsave(file, g, width = 12, height = 8, dpi = 300) } )

# ———- PETA pal_fun <- function(x) colorNumeric(“viridis”, domain = x) render_leaflet_map <- function(varname) { x <- jabar[[varname]] pal <- pal_fun(x); cols <- pal(x) labels <- sprintf( “%s
%s: %s”, jabar$kabupaten, var_labels[[varname]], format(round(x, 3), nsmall = 3) ) |> lapply(HTML) leaflet(jabar) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons( fillColor = cols, color = “white”, weight = 1, fillOpacity = 0.85, label = labels, labelOptions = labelOptions(direction = “auto”), highlightOptions = highlightOptions(weight = 2, bringToFront = TRUE) ) %>% addLegend(“bottomright”, pal = pal, values = x, title = var_labels[[varname]]) } output$map_left <- renderLeaflet({ render_leaflet_map(req(input$var_map_left)) }) output$map_right <- renderLeaflet({ render_leaflet_map(req(input$var_map_right)) })

make_static_choropleth <- function(varname) { ggplot(jabar) + geom_sf(aes(fill = .data[[varname]]), color = “white”, size = 0.25) + scale_fill_viridis_c(name = var_labels[[varname]], option = “C”, direction = -1) + labs(title = paste(“Peta”, var_labels[[varname]])) + theme_minimal(base_size = 12, base_family = “Inter”) + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } output$dl_map_left_png <- downloadHandler( filename = function() sprintf("peta_kiri_%s.png", input$var_map_left), content = function(file) { p <- make_static_choropleth(isolate(input$var_map_left)) ggsave(file, p, width = 7, height = 6, dpi = 300) } ) output$dl_map_right_png <- downloadHandler( filename = function() sprintf(“peta_kanan_%s.png”, input$var_map_right), content = function(file) { p <- make_static_choropleth(isolate(input$var_map_right)) ggsave(file, p, width = 7, height = 6, dpi = 300) } )

# ———- AUTOKORELASI output$auto_table <- renderDT({ v <- req(input$var_auto) if (input$auto_method == "Moran's I") { res <- safe_moran_both_raw(jabar[[v]], nb_queen) Moran_I <- unname(res$mt$estimate[["Moran I statistic"]]) Expected_I <- unname(res$mt$estimate[["Expectation"]]) Z_value <- unname(res$mt$statistic) P_perm <- res$mmc$p.value tab <- tibble::tibble( `Moran's I` = Moran_I, `Expected I` = Expected_I, `Z-value` = Z_value, `p-value Uji Permutasi` = P_perm, `Keputusan` = ifelse(P_perm < 0.05, "Signifikan", "Tidak signifikan") ) %>% round3_df() DT::datatable(tab, options = list(dom = 't')) } else if (input$auto_method == “Geary’s C”) { res <- safe_geary_both_raw(jabar[[v]], nb_queen) estC <- res$gt$estimate[[“Geary C statistic”]]; if (is.null(estC)) estC <- res$gt$estimate[[“Geary’s C”]] Expected_C <- unname(res$gt$estimate[[“Expectation”]]) Z_value <- unname(res$gt$statistic) P_perm <- res$gmc$p.value tab <- tibble::tibble( Geary's C = unname(estC), Expected C = Expected_C, Z-value = Z_value, p-value Uji Permutasi = P_perm, Keputusan = ifelse(P_perm < 0.05, “Signifikan”, “Tidak signifikan”) ) %>% round3_df() DT::datatable(tab, options = list(dom = ‘t’)) } else if (input$auto_method == "LISA") { lis <- compute_lisa(jabar[[v]], nb_queen) sf_sub <- jabar sf_sub$LISA_cat <- factor(“Not Significant”, levels = c(“High-High”,“Low-Low”,“High-Low”,“Low-High”,“Not Significant”)) sf_sub$LISA_cat[lis$ok] <- factor(lis$cluster, levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")) tab <- sf_sub |> st_drop_geometry() |> transmute( Kabupaten = kabupaten, Kategori = as.character(LISA_cat) ) keep <- input$lisa_filter if (!is.null(keep)) tab <- tab %>% filter(Kategori %in% keep) DT::datatable(tab, options = list(pageLength = 15)) } else { # Getis-Ord G* gs <- compute_gistar(jabar[[v]], nb_queen) tab <- jabar |> st_drop_geometry() |> transmute( Kabupaten = kabupaten, Z(G*) = round(gs$z, 3), Kategori = as.character(gs$cat) ) keep <- input$gstar_filter if (!is.null(keep)) tab <- tab %>% filter(Kategori %in% keep) DT::datatable(tab, options = list(pageLength = 15)) } })

build_auto_leaflet <- function(method, v) { if (method == “LISA”) { lis <- compute_lisa(jabar[[v]], nb_queen) sf_sub <- jabar sf_sub$LISA_cat <- factor("Not Significant", levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")) sf_sub$LISA_cat[lis$ok] <- factor(lis$cluster, levels = c(“High-High”,“Low-Low”,“High-Low”,“Low-High”,“Not Significant”)) # Warna LISA tegas (non-pastel) cols_lisa <- c(“High-High”=“#E53935”,“Low-Low”=“#1E88E5”,“High-Low”=“#FB8C00”,“Low-High”=“#43A047”,“Not Significant”=“#B0BEC5”) pal <- colorFactor(cols_lisa, levels = names(cols_lisa)) labels <- sprintf(“%s
Kategori: %s”, sf_sub$kabupaten, as.character(sf_sub$LISA_cat)) |> lapply(HTML) leaflet(sf_sub) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(fillColor = pal(sf_sub$LISA_cat), color=“white”, weight=1, fillOpacity=0.85, label=labels, labelOptions = labelOptions(direction=“auto”), highlightOptions = highlightOptions(weight = 2, bringToFront=TRUE)) %>% addLegend(“bottomright”, pal = pal, values = sf_sub$LISA_cat, title = "LISA") } else if (method == "Getis-Ord G*") { gs <- compute_gistar(jabar[[v]], nb_queen) sf_g <- jabar %>% mutate(GSTAR = gs$cat, z_g = gs$z) cols_g <- c("Hot spot"="#E53935","Cold spot"="#1E88E5","Not significant"="#CFD8DC") pal <- colorFactor(cols_g, levels = names(cols_g)) labels <- sprintf("%s Kategori: %s Z(G*): %s", sf_g$kabupaten, as.character(sf_g$GSTAR), format(round(sf_g$z_g,3), nsmall=3)) |> lapply(HTML) leaflet(sf_g) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(fillColor = pal(sf_g$GSTAR), color=“white”, weight=1, fillOpacity=0.85, label=labels, labelOptions = labelOptions(direction=“auto”), highlightOptions = highlightOptions(weight=2, bringToFront=TRUE)) %>% addLegend(“bottomright”, pal = pal, values = sf_g$GSTAR, title = "Getis–Ord G*") } else if (method == "Moran's I") { res <- safe_moran_both_raw(jabar[[v]], nb_queen) Moran_I <- round(unname(res$mt$estimate[["Moran I statistic"]]), 3) P_perm <- round(res$mmc$p.value, 3) x <- jabar[[v]] pal <- colorNumeric("viridis", domain = x) labels <- sprintf("%s %s: %s Moran's I: %s p (perm): %s", jabar$kabupaten, var_labels[[v]], format(round(x,3), nsmall=3), Moran_I, P_perm) |> lapply(HTML) leaflet(jabar) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(fillColor = pal(x), color="white", weight=1, fillOpacity=0.85, label=labels, labelOptions = labelOptions(direction="auto"), highlightOptions = highlightOptions(weight=2, bringToFront=TRUE)) %>% addLegend("bottomright", pal = pal, values = x, title = var_labels[[v]]) } else { # Geary res <- safe_geary_both_raw(jabar[[v]], nb_queen) Geary_C <- round(unname(res$gt$estimate[["Geary C statistic"]]), 3) if (is.na(Geary_C)) Geary_C <- round(unname(res$gt$estimate[["Geary's C"]]), 3) P_perm <- round(res$gmc$p.value, 3) x <- jabar[[v]] pal <- colorNumeric("viridis", domain = x) labels <- sprintf("%s %s: %s Geary's C: %s p (perm): %s", jabar$kabupaten, var_labels[[v]], format(round(x,3), nsmall=3), Geary_C, P_perm) |> lapply(HTML) leaflet(jabar) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(fillColor = pal(x), color="white", weight=1, fillOpacity=0.85, label=labels, labelOptions = labelOptions(direction="auto"), highlightOptions = highlightOptions(weight=2, bringToFront=TRUE)) %>% addLegend("bottomright", pal = pal, values = x, title = var_labels[[v]]) } } output$auto_map <- renderLeaflet({ build_auto_leaflet(input$auto_method, req(input$var_auto)) })

make_static_auto_plot <- function(method, v) { if (method == “LISA”) { lis <- compute_lisa(jabar[[v]], nb_queen) sf_sub <- jabar sf_sub$LISA_cat <- factor("Not Significant", levels = c("High-High","Low-Low","High-Low","Low-High","Not Significant")) sf_sub$LISA_cat[lis$ok] <- factor(lis$cluster, levels = c(“High-High”,“Low-Low”,“High-Low”,“Low-High”,“Not Significant”)) cols_lisa <- c(“High-High”=“#E53935”,“Low-Low”=“#1E88E5”,“High-Low”=“#FB8C00”,“Low-High”=“#43A047”,“Not Significant”=“#B0BEC5”) ggplot(sf_sub) + geom_sf(aes(fill = LISA_cat), color = “white”, size = 0.25) + scale_fill_manual(values = cols_lisa, drop = FALSE, name = “Kategori”) + labs(title = paste(“LISA —”, var_labels[[v]])) + theme_minimal(base_size = 12, base_family = “Inter”) + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } else if (method == “Getis-Ord G*”) { gs <- compute_gistar(jabar[[v]], nb_queen) sf_g <- jabar %>% mutate(GSTAR = gs\(cat) cols_g <- c("Hot spot"="#E53935","Cold spot"="#1E88E5","Not significant"="#CFD8DC") ggplot(sf_g) + geom_sf(aes(fill = GSTAR), color = "white", size = 0.25) + scale_fill_manual(values = cols_g, drop = FALSE, name = "Kategori") + labs(title = paste("Getis–Ord G* —", var_labels[[v]])) + theme_minimal(base_size = 12, base_family = "Inter") + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } else if (method == "Moran's I") { ggplot(jabar) + geom_sf(aes(fill = .data[[v]]), color = "white", size = 0.25) + scale_fill_viridis_c(name = var_labels[[v]], option = "C", direction = -1) + labs(title = paste("Moran's I —", var_labels[[v]])) + theme_minimal(base_size = 12, base_family = "Inter") + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } else { ggplot(jabar) + geom_sf(aes(fill = .data[[v]]), color = "white", size = 0.25) + scale_fill_viridis_c(name = var_labels[[v]], option = "C", direction = -1) + labs(title = paste("Geary's C —", var_labels[[v]])) + theme_minimal(base_size = 12, base_family = "Inter") + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } } output\)dl_auto_map_png <- downloadHandler( filename = function() sprintf(“autokorelasi_%s_%s.png”, input$auto_method, input$var_auto), content = function(file) { p <- make_static_auto_plot(isolate(input$auto_method), isolate(input$var_auto)) ggsave(file, p, width = 7.2, height = 6.2, dpi = 300) } )

# ———- PEMODELAN current_w <- reactive({ switch(input$w_choice, “Queen” = lw_queen, “Rook” = lw_rook, “KNN-4” = lw_knn4 ) })

models_reactive <- eventReactive(input$fit_btn, { y <- req(input$y_model); x <- req(input$x_model) validate(need(length(x) >= 1, “Pilih minimal satu X”)) fit_models(jabar, y, x, current_w()) }, ignoreInit = TRUE)

output$model_table <- renderDT({ req(models_reactive()) yv <- jabar[[input$y_model]] df <- model_stats(models_reactive(), yv)
rng_aic <- range(df$AIC, na.rm = TRUE)
rng_ll <- range(df$`Log-Likelihood`, na.rm = TRUE)
rng_r2 <- range(df$`R²`, na.rm = TRUE)

df_render <- df %>%
 transmute(
 Model,
 `Log-Likelihood (bar)` = vapply(`Log-Likelihood`, bar_cell, character(1),
 minv = rng_ll[1], maxv = rng_ll[2],
 color = "#2e7d32", txt_color = "#0d1b2a"),
 `Log-Likelihood_num` = `Log-Likelihood`,
 `AIC (bar)` = vapply(AIC, bar_cell, character(1),
 minv = rng_aic[1], maxv = rng_aic[2],
 color = "#1f2d5a", txt_color = "#e3eaf2"),
 AIC_num = AIC,
 `R² (bar)` = vapply(`R²`, bar_cell, character(1),
 minv = rng_r2[1], maxv = rng_r2[2],
 color = "#f2a900", txt_color = "#0d1b2a"),
 R2_num = `R²`
 )

DT::datatable(
 df_render,
 escape = FALSE,
 rownames = FALSE,
 options = list(
 pageLength = 8,
 dom = 'tip',
 columnDefs = list(
 list(targets = 1, orderData = 2), # sort LL by hidden numeric
 list(visible = FALSE, targets = 2),
 list(targets = 3, orderData = 4), # sort AIC by hidden numeric
 list(visible = FALSE, targets = 4),
 list(targets = 5, orderData = 6), # sort R2 by hidden numeric
 list(visible = FALSE, targets = 6)
 )
 )
)
})

output$best_model_cards <- renderUI({ req(models_reactive()) yv <- jabar[[input$y_model]] df <- model_stats(models_reactive(), yv) best_aic <- df$Model[which.min(df$AIC)] best_r2 <- df$Model[which.max(df$R²)] best_ll <- df$Model[which.max(df$Log-Likelihood)] tagList( tags$b(sprintf("Model dengan AIC terendah: %s", best_aic)), tags$br(), tags$b(sprintf("Model dengan R² tertinggi: %s", best_r2)), tags$br(), tags$b(sprintf(“Model dengan Log-Likelihood tertinggi: %s”, best_ll)) ) })

pred_resid_from_model <- reactive({ req(models_reactive()) mods <- models_reactive() mdl <- switch(input$model_pred, "SAR" = mods$sar, “SDM” = mods$sdm, "SAC/SARAR" = mods$sac, “GNS” = mods$gns, "SLX" = mods$slx, “SEM” = mods$sem, "SDEM" = mods$sdem ) data.frame( yhat = as.numeric(fitted(mdl)), resid = as.numeric(residuals(mdl)) ) })

build_resid_leaflet <- reactive({ df <- req(pred_resid_from_model()) jabar_res <- jabar %>% mutate( y_hat = round(df$yhat, 3), z_resid = round(as.numeric(scale(df$resid)), 3) ) rng <- range(jabar_res$z_resid, na.rm = TRUE) lim <- max(abs(rng)) pal_res <- colorNumeric(c("#1E88E5","white","#E53935"), domain = c(-lim, lim)) labels <- sprintf( "%s Prediksi: %.3f Aktual: %.3f Residual (z): %.3f", jabar_res$kabupaten, jabar_res$y_hat, jabar_res$KEMISKINAN, jabar_res$z_resid ) |> lapply(HTML) leaflet(jabar_res) %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons( fillColor = pal_res(jabar_res$z_resid), color = "white", weight = 1, fillOpacity = 0.85, label = labels, labelOptions = labelOptions(direction="auto"), highlightOptions = highlightOptions(weight = 2, bringToFront = TRUE) ) %>% addLegend("bottomright", colors = c("#1E88E5","white","#E53935"), labels = c("Negatif besar", "0", "Positif besar"), title = "Residual (z)") }) output$resid_leaf <- renderLeaflet({ build_resid_leaflet() })

make_static_resid_plot <- function() { df <- req(pred_resid_from_model()) sfp <- jabar %>% mutate(z_resid = round(as.numeric(scale(df$resid)), 3)) lim_res <- max(abs(sfp$z_resid), na.rm = TRUE) ggplot(sfp) + geom_sf(aes(fill = z_resid), color = “white”, size = 0.25) + scale_fill_gradient2(low = “#1E88E5”, mid = “white”, high = “#E53935”, midpoint = 0, limits = c(-lim_res, lim_res), oob = scales::squish, name = “Residual (z)”) + labs(title = paste(“Peta Residual —”, input$model_pred)) + theme_minimal(base_size = 12, base_family = "Inter") + theme(axis.title = element_blank(), axis.text = element_blank(), panel.grid = element_blank()) } output$dl_resid_map_png <- downloadHandler( filename = function() sprintf(“peta_residual_%s.png”, input$model_pred), content = function(file) { p <- make_static_resid_plot() ggsave(file, p, width = 7.2, height = 6.2, dpi = 300) } )

# ———- SENSITIVITAS W output$sens_table <- renderDT({ y <- req(input$y_model); x <- req(input$x_model) validate(need(length(x) >= 1, "Pilih minimal satu X lalu klik Fit.")) fml <- as.formula(paste(y, "~", paste(x, collapse = " + "))) mQ <- lagsarlm(fml, data = jabar, listw = lw_queen, method = "eigen", zero.policy = TRUE) mR <- lagsarlm(fml, data = jabar, listw = lw_rook, method = "eigen", zero.policy = TRUE) mK <- lagsarlm(fml, data = jabar, listw = lw_knn4, method = "eigen", zero.policy = TRUE) tab <- data.frame( Bobot = c("Queen","Rook","KNN-4"), Rho = round(c(mQ$rho, mR$rho, mK$rho), 3), AIC = round(c(AIC(mQ), AIC(mR), AIC(mK)), 3), Log-Likelihood = round(c(as.numeric(logLik(mQ)), as.numeric(logLik(mR)), as.numeric(logLik(mK))), 3) ) DT::datatable(tab, rownames = FALSE, options = list(pageLength = 5)) })

output$sens_bar <- renderPlot({ y <- req(input$y_model); x <- req(input$x_model) validate(need(length(x) >= 1, "Pilih minimal satu X lalu klik Fit.")) fml <- as.formula(paste(y, "~", paste(x, collapse = " + "))) mQ <- lagsarlm(fml, data = jabar, listw = lw_queen, method = "eigen", zero.policy = TRUE) mR <- lagsarlm(fml, data = jabar, listw = lw_rook, method = "eigen", zero.policy = TRUE) mK <- lagsarlm(fml, data = jabar, listw = lw_knn4, method = "eigen", zero.policy = TRUE) df <- data.frame( Bobot = factor(c("Queen","Rook","KNN-4"), levels = c("Queen","Rook","KNN-4")), Rho = c(mQ$rho, mR$rho, mK$rho), AIC = c(AIC(mQ), AIC(mR), AIC(mK)) ) pAIC <- ggplot(df, aes(x = Bobot, y = AIC)) + geom_col(width = 0.6, fill = “#1f2d5a”) + labs(title = “AIC per Bobot W”) + theme_minimal(base_size = 12, base_family = “Inter”) pRho <- ggplot(df, aes(x = Bobot, y = Rho)) + geom_col(width = 0.6, fill = “#43A047”) + labs(title = “Rho (ρ) per Bobot W”) + theme_minimal(base_size = 12, base_family = “Inter”) gridExtra::grid.arrange(pAIC, pRho, ncol = 2) }) output$dl_sens_bar <- downloadHandler( filename = function() "sensitivitas_W.png", content = function(file) { y <- isolate(input$y_model); x <- isolate(input$x_model) fml <- as.formula(paste(y, "~", paste(x, collapse = " + "))) mQ <- lagsarlm(fml, data = jabar, listw = lw_queen, method = "eigen", zero.policy = TRUE) mR <- lagsarlm(fml, data = jabar, listw = lw_rook, method = "eigen", zero.policy = TRUE) mK <- lagsarlm(fml, data = jabar, listw = lw_knn4, method = "eigen", zero.policy = TRUE) df <- data.frame( Bobot = factor(c("Queen","Rook","KNN-4"), levels = c("Queen","Rook","KNN-4")), Rho = c(mQ$rho, mR$rho, mK$rho), AIC = c(AIC(mQ), AIC(mR), AIC(mK)) ) pAIC <- ggplot(df, aes(x = Bobot, y = AIC)) + geom_col(width = 0.6, fill = “#1f2d5a”) + labs(title = “AIC per Bobot W”) + theme_minimal(base_size = 12, base_family = “Inter”) pRho <- ggplot(df, aes(x = Bobot, y = Rho)) + geom_col(width = 0.6, fill = “#43A047”) + labs(title = “Rho (ρ) per Bobot W”) + theme_minimal(base_size = 12, base_family = “Inter”) g <- gridExtra::arrangeGrob(pAIC, pRho, ncol = 2) ggsave(file, g, width = 10, height = 5.5, dpi = 300) } ) }

=========================================================

RUN

=========================================================

shinyApp(ui, server)

PEMODELAN EKONOMETRIKA SPASIAL TERHADAP FAKTOR SOSIAL-EKONOMI PENENTU TINGKAT KEMISKINAN DI JAWA BARAT

Humaira Zeanova

BAB I PENDAHULUAN

1.1 Latar Belakang Masalah

1.2 Identifikasi Masalah

1.3. Maksud dan Tujuan Penelitian

1.4 Batasan Penelitian

BAB II TINJAUAN PUSTAKA

2.1 Spatial Dependence

2.2 Autokorelasi Spasial

2.3 Ordinary Least Squares

2.4 Uji Lagrange Multiplier

2.5 Spatial Autoregressive Model

2.6 Akaike Information Criterion

BAB III METODOLOGI PENELITIAN

3.1 Data

3.2 Metode Analisis

BAB IV HASIL DAN PEMBAHASAN

4.1 Analisis Deskriptif

4.2 Ordinary Least Square untuk Model Klasik

4.3 Uji Moran’s I Residual OLS

4.4 Autokorelasi Spasial

4.4.1 Moran’s I

4.4.2 Geary’s C

4.4.3 Local Moran’s I (LISA)

4.4.4 Getis–Ord 𝐺∗𝑖

4.5 Uji Lagrange Multiplier

4.6 Perbandingan Model

4.7 Estimasi Model SAR

4.8 Uji Asumsi Model SAR

4.9 Uji Sensitivitas Bobot

4.10 Interpretasi Hasil

BAB V KESIMPULAN

5.1 Kesimpulan

5.2 Saran

DAFTAR PUSTAKA

LAMPIRAN

Script Analisis

Script Dashboard

=========================================================

app.R — West Java Poverty Dashboard 2024 (VERSI 2)

(Navy theme, Inter font, Unpad logo, proper sorting with in-cell bars)

=========================================================

————————-

GLOBAL: Fonts & Assets

————————-

Inter font (UI & plot)

————————-

GLOBAL: Data & helpers

————————-

Weights

=========================================================

UI

=========================================================

=========================================================

SERVER

=========================================================

=========================================================

RUN

=========================================================