PEMILIHAN MODEL TERBAIK DALAM MENGANALISIS FAKTOR EKONOMI DAN SOSIAL TERHADAP TINGKAT KEMISKINAN DI JAWA BARAT TAHUN 2024
Diajukan untuk Memenuhi Nilai Mata Kuliah Spatial Statistics TA. 2025/2026
ARA FIGURA BILBINA_ 140610230004_STATISTIKA UNIVERSITAS PADJADJARAN
2025-10-31
BAB I. PENDAHULUAN
1.1 Latar Belakang Masalah
Kemiskinan mencerminkan rendahnya tingkat pendapatan seseorang
sehingga tidak mampu memenuhi standar hidup yang layak (Corral et al.,
2020). Kondisi ini, jika berlangsung dalam jangka panjang, dapat
menimbulkan berbagai dampak sosial dan ekonomi seperti meningkatnya
ketimpangan sosial, terhambatnya pembangunan, munculnya konflik
horizontal, serta naiknya tingkat kriminalitas (UN, 2023). Apabila tidak
diatasi secara tepat, kemiskinan berpotensi menciptakan siklus
antargenerasi yang sulit diputus (Rivera & Currais, 1999; Benhabib
& Spiegel, 1994; Lucas, 1990).
Meskipun tingkat kemiskinan di Indonesia, khususnya di Provinsi Jawa Barat, terus mengalami penurunan hingga mencapai 7,46% pada tahun 2024, persentase penduduk miskin dalam periode 2020–2024 cenderung stagnan. Hal ini menunjukkan perlunya analisis lebih mendalam terhadap faktor-faktor yang memengaruhi kemiskinan agar kebijakan penanggulangannya di Jawa Barat dapat dirancang secara lebih efektif dan berkelanjutan.
1.2 Identifikasi Masalah
Pemetaan persentase kemiskinan di kabupaten/kota Provinsi Jawa Barat menunjukkan adanya pola spasial yang jelas, di mana daerah dengan tingkat kemiskinan tinggi cenderung berdekatan dengan daerah lain yang juga memiliki tingkat kemiskinan tinggi, begitu pula sebaliknya. Misalnya, beberapa wilayah di bagian selatan Jawa Barat seperti Kabupaten Garut, Tasikmalaya, dan Cianjur memperlihatkan tingkat kemiskinan yang relatif lebih tinggi dibandingkan wilayah utara seperti Kota Bandung atau Kota Bekasi. Pola pengelompokan ini mencerminkan adanya pola spasial, yaitu kecenderungan nilai suatu wilayah dipengaruhi oleh nilai wilayah tetangganya.
Fenomena tersebut dapat terjadi karena adanya kesamaan karakteristik sosial-ekonomi, seperti keterbatasan akses terhadap infrastruktur, lapangan kerja, dan kualitas pendidikan di wilayah yang berdekatan. Selain itu, interaksi antarwilayah seperti arus tenaga kerja, perdagangan lokal, maupun keterkaitan geografis juga dapat memperkuat keterhubungan antar tingkat kemiskinan. Dengan demikian, kondisi ini menunjukkan bahwa data kemiskinan di Jawa Barat tidak bersifat acak melainkan saling bergantung secara spasial (spatial dependence).
Ketika dependensi spasial ini diabaikan dan model regresi klasik (OLS) tetap digunakan, hasil estimasi yang diperoleh berpotensi bias dan tidak efisien karena pelanggaran asumsi independensi antar observasi. Oleh sebab itu, penggunaan pendekatan ekonometrika spasial menjadi penting untuk ditinjau sehingga akan menggambarkan hubungan spasial secara lebih akurat serta memahami bagaimana kemiskinan di suatu daerah dipengaruhi oleh kondisi di daerah sekitarnya.
1.3 Maksud dan Tujuan Penelitian
Maksud dari penelitian ini adalah untuk 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 kabupaten/kota. Analisis difokuskan pada pola spasial kemiskinan tanpa membahas aspek penyebab non-spasial secara mendalam. Selain itu, penelitian ini hanya menggunakan data tahun 2024 sehingga hasilnya menggambarkan kondisi pada periode tertentu dan tidak mencakup dinamika perubahan antarwaktu.
BAB II. TINJAUAN PUSTAKA
2.1 Spatial Autocorrelation
Dalam analisis spasial, keterkaitan antar wilayah secara spasial diukur menggunakan Moran’s I. Nilai I > 0 menunjukkan adanya autokorelasi spasial positif, yaitu wilayah dengan nilai mirip (high-high atau low-low) cenderung berdekatan. Sebaliknya, I < 0 menandakan pola checkerboard (high-low), sedangkan I ≈ 0 menunjukkan pola yang acak. Dalam perhitungannya, Moran’s I memanfaatkanmatriks bobot spasial (W) yang menentukan seberapa kuat pengaruh satu wilayah terhadap wilayah lain berdasarkan kriteria kedekatan tertentu, yakni Rook, Queen, ataupun Jarak. Selain Moran’s I, terdapat metode pengukuran lain seperti Greary’s C untuk ukuran Global, serta Local Moran’s I dan Getis-Ord G* untuk ukuran lokal.
2.2 OLS (Ordinary Least Square)
OLS adalah metode estimasi parameter untuk regresi klasik. Metode ini membutuhkan pemenuhan pengujian asumsi klasik, diantaranya pengujian multikolinearitas, residual berdistribusi normal, dan homoskedastisitas. Metode ini mengasumsikan bahwa error antarwilayah bersifat independen dan homogen, sehingga tidak dapat menangkap spillover atau pengaruh kemiskinan dari satu daerah ke daerah lain. Hal ini menyebabkan hasil estimasi menjadi bias dan kurang akurat. Untuk mengatasi keterbatasan tersebut, digunakan model ekonometrika spasial yang mampu memperhitungkan spatial dependence antarwilayah.
Adapun model matematis OLS adalah sebagai berikut: \[ Y = X\beta + \varepsilon \]
2.3 Uji Lagrange Multiplier
Uji ini digunakan untuk mengidentifikasi ada tidaknya dependensi spasial pada model. Jika hasil uji menunjukkan adanya dependensi lag, maka digunakan model Spatial Autoregressive (SAR). Jika yang muncul adalah dependensi galat, maka digunakan model Spatial Error Model (SEM). Uji ini membantu menentukan model spasial yang paling sesuai sebelum analisis lanjut.
Uji Pengganda Lagrange pada Lag
\[ H_0 : \rho = 0, \quad \text{(tidak terdapat dependensi spasial pada lag)} \] \[ H_1 : \rho \ne 0, \quad \text{(terdapat dependensi spasial pada lag)} \]
Uji Pengganda Lagrange pada Galat
\[ H_0 : \lambda = 0, \quad \text{(tidak terdapat dependensi spasial pada galat)} \] \[ H_1 : \lambda \ne 0, \quad \text{(terdapat dependensi spasial pada galat)} \]
2.4 Spatial Autoregressive (SAR)
Menurut Anselin (2003), Spatial Autoregressive Model (SAR) merupakan pengembangan dari model OLS dengan menambahkan pengaruh lag pada variabel dependen. Model ini digunakan ketika λ = 0, sehingga hanya mempertimbangkan proses autoregresif pada variabel respon. Dengan kata lain, SAR menunjukkan bahwa nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah-wilayah tetangganya, sesuai dengan hubungan spasial yang terbentuk dalam model.
\[ Y = \rho W Y + X\beta + \varepsilon \]
Asumsi penting dalam model SAR adalah bahwa residual tidak mengalami autokorelasi spasial, sehingga pengaruh spasial telah sepenuhnya ditangkap oleh komponen lag variabel dependen. Jika asumsi ini terpenuhi, maka model SAR dianggap tepat untuk menggambarkan struktur keterkaitan spasial dalam data.
2.5 Akaike Information Criterion (AIC)
Menurut Suryaningrat, S. H. (2021) Akaike Information Criterion (AIC) merupakan metode analisis yang digunakan untuk memperoleh model yang terbaik dengan menggunakan estimasi maximum likelihood sebagai perhitungan yang sesuai. AIC mempertimbangkan keseimbangan antara kompleksitas model (jumlah parameter) dan kecocokan model terhadap data (likelihood).
AIC dihitung dengan rumus:
\[ AIC = 2k - 2\ln(\hat{L}) \]
di mana k adalah jumlah parameter dalam model dan L adalah nilai maksimum dari fungsi likelihood. Model dengan nilai AIC yang lebih rendah dianggap lebih baik karena menunjukkan kompleksitas yang wajar namun 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 mencerminkan kondisi sosial ekonomi daerah yang diduga berpengaruh terhadap tingkat kemiskinan.
Variabel yang digunakan terdiri atas satu variabel dependen yaitu kemiskinan, yang diukur melalui persentase penduduk miskin menurut kabupaten/kota, serta empat variabel independen yaitu Tingkat Pengangguran Terbuka (TPT), Produk Domestik Regional Bruto (PDRB) per kapita, Indeks Pembangunan Manusia (IPM), dan Gini Ratio. Keempat variabel tersebut dipilih karena secara teori memiliki hubungan erat dengan tingkat kesejahteraan dan ketimpangan sosial ekonomi masyarakat.
Selain data atribut, penelitian ini juga menggunakan data spasial berupa peta batas administrasi kabupaten/kota di Jawa Barat yang diperoleh dari geoBoundaries Indonesia Administrative Level 2 (ADM2) melalui situs www.geoboundaries.org dalam format shapefile (.shp). Data ini digunakan untuk membentuk hubungan antarwilayah dalam analisis spasial.
3.2 Metode Analisis
Analisis dilakukan menggunakan pendekatan statistik spasial dengan bantuan perangkat lunak R. Pendekatan ini bertujuan untuk mengidentifikasi pengaruh faktor sosial ekonomi terhadap tingkat kemiskinan serta mendeteksi adanya ketergantungan spasial antarwilayah di Jawa Barat. Analisis diawali dengan uji autokorelasi spasial global (Moran’s I, Geary’s C) dan lokal (LISA serta Getis-Ord Gi*) untuk melihat pola keterkaitan antarwilayah.
Selanjutnya digunakan model regresi OLS sebagai dasar analisis hubungan antarvariabel, diikuti uji Lagrange Multiplier (LM) untuk menentukan model spasial yang tepat. Beberapa model yang dipertimbangkan mencakup SAR, SDM, SAC/SARAR, dan GNS. Pemilihan model terbaik dilakukan berdasarkan nilai AIC, signifikansi parameter spasial (ρ dan λ), serta kesesuaian teori yang mendukung hubungan spasial antarwilayah.
3.3 Alur Kerja Penelitian
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 berikut:
| Variabel | Rata_rata | Minimum | Maksimum | 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 |
Berdasarkan hasil analisis deskriptif, rata-rata tingkat kemiskinan di Jawa Barat sebesar 8,01%, dengan nilai terendah 2,34% dan tertinggi 11,93%, serta simpangan baku sebesar 2,65, yang menunjukkan adanya variasi cukup besar antarwilayah. Variasi ini tampak jelas pada peta persebaran, di mana wilayah selatan seperti Kabupaten Garut, Tasikmalaya, dan Cianjur memiliki tingkat kemiskinan yang lebih tinggi dibandingkan wilayah perkotaan seperti Kota Bandung, Kota Bekasi, dan Kota Depok.
Tingkat Pengangguran Terbuka (TPT) memiliki rata-rata
6,52%, dengan sebaran yang juga menunjukkan ketimpangan
antarwilayah, di mana daerah perkotaan cenderung memiliki TPT lebih
tinggi akibat tingginya kompetisi kerja di sektor formal.
PDRB per kapita rata-rata mencapai 53,11 juta rupiah,
dengan rentang nilai yang lebar (24,91–107,08 juta) dan simpangan baku
tinggi (33,26), mengindikasikan adanya kesenjangan
ekonomi antar daerah. Wilayah seperti Kota Bekasi dan Kota Bandung
memiliki nilai PDRB yang jauh lebih tinggi dibandingkan daerah pedesaan
di selatan.
Sementara itu, IPM rata-rata sebesar 74,68,
menandakan bahwa sebagian besar kabupaten/kota di Jawa Barat telah
mencapai tingkat pembangunan manusia yang cukup baik. Namun demikian,
masih terdapat ketimpangan antarwilayah, dengan nilai IPM lebih rendah
di bagian selatan provinsi.
Rasio Gini yang rata-rata sebesar 0,37 menunjukkan
tingkat ketimpangan pendapatan yang tergolong sedang, dengan variasi
kecil antarwilayah (simpangan baku 0,05).
Secara keseluruhan, hasil deskriptif dan peta tematik menunjukkan adanya perbedaan yang nyata antarwilayah dalam aspek sosial dan ekonomi. Pola spasial yang muncul, seperti pengelompokan wilayah dengan tingkat kemiskinan tinggi di selatan dan tingkat kemiskinan rendah di utara, mengindikasikan adanya spatial dependence atau keterkaitan spasial antarwilayah di Jawa Barat. Hal ini menjadi dasar penting untuk melakukan analisis lanjutan menggunakan model ekonometrika spasial.
4.2 Ordinary Least Square untuk Model Klasik
\[ Y = 44.770 - 0.0879(\text{TPT}) - 0.0046(\text{PDRB}) - 0.5521(\text{IPM})^{***} + 14.2366(\text{GINI}) \]
4.2.1 Uji Linearitas
| Statistik_Uji | F | p | Ket |
|---|---|---|---|
| Ramsey RESET | 0.476 | 0.6281 | Hubungan linier |
4.2.2 Uji Normalitas
| Statistik_Uji | W | p | Ket |
|---|---|---|---|
| Shapiro–Wilk | 0.9552 | 0.2865 | Data berdistribusi normal |
4.2.3 Uji Homoskedastisitas
| Statistik_Uji | BP | p | Ket |
|---|---|---|---|
| Breusch-Pagan Test | 2.56 | 0.6339 | Tidak terdapat heteroskedastisitas pada data. |
4.2.4 Uji Moran’s I Residual
| Statistik_Uji | I | p_Monte_Carlo | Ket |
|---|---|---|---|
| Moran’s I | 2.3277 | 0.015 | Terdapat autokorelasi spasial pada error |
Berdasarkan hasil estimasi model Ordinary Least Square (OLS), seluruh asumsi klasik pada model telah terpenuhi, meliputi uji linearitas, normalitas, dan homoskedastisitas. Hasil uji Ramsey RESET menunjukkan nilai p-value sebesar 0,6281, yang mengindikasikan bahwa hubungan antara variabel dependen dan independen bersifat linier. Uji Shapiro–Wilk menghasilkan p-value 0,2865, sehingga residual berdistribusi normal. Selain itu, uji Breusch–Pagan dengan p-value 0,6339 menunjukkan tidak adanya masalah heteroskedastisitas pada model.
Namun, hasil uji Moran’s I pada residual memberikan nilai I = 2,3277 dengan p-value 0,015, yang berarti terdapat autokorelasi spasial pada residual. Dengan demikian, meskipun model OLS memenuhi asumsi klasik, adanya indikasi dependensi spasial menandakan bahwa model ini belum sepenuhnya sesuai digunakan untuk data yang memiliki keterkaitan antarwilayah. Oleh karena itu, perlu dilakukan analisis lanjutan menggunakan model ekonometrika spasial seperti SAR, SAC/SARAR, GNS atau SDM untuk menangkap efek spasial secara lebih akurat.
4.3 Autokorelasi Spasial
4.3.1 Moran’s I
| Variabel | Morans_I | Z | p_permutasi | Ket |
|---|---|---|---|---|
| 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 | 0.299 | 2.40 | 0.036 | Signifikan |
Berdasarkan hasil uji Moran’s I untuk masing-masing variabel, dapat diketahui bahwa sebagian besar variabel menunjukkan adanya autokorelasi spasial yang signifikan, yang berarti nilai suatu wilayah cenderung dipengaruhi oleh nilai di wilayah sekitarnya.
Nilai Moran’s I untuk variabel kemiskinan sebesar 0,283 (p = 0,001) menunjukkan adanya keterkaitan spasial positif yang signifikan. Artinya, kabupaten/kota dengan tingkat kemiskinan tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki tingkat kemiskinan tinggi, dan sebaliknya. Pola ini menegaskan adanya pengelompokan (spatial clustering) kemiskinan di Jawa Barat.
Variabel Tingkat Pengangguran Terbuka (TPT) memiliki nilai Moran’s I paling tinggi, yaitu 0,511 (p = 0,001), yang mengindikasikan keterkaitan spasial yang sangat kuat. Hal ini menunjukkan bahwa wilayah dengan tingkat pengangguran tinggi cenderung berdekatan dengan wilayah lain yang memiliki karakteristik serupa, kemungkinan karena kesamaan struktur ekonomi dan kesempatan kerja antarwilayah berdekatan.
Variabel PDRB per kapita memiliki nilai Moran’s I 0,211 dengan p-value 0,102, yang berarti tidak signifikan secara spasial. Artinya, distribusi PDRB per kapita antarwilayah di Jawa Barat relatif acak dan tidak menunjukkan pola pengelompokan spasial yang kuat — kemungkinan karena perbedaan sektor unggulan antarwilayah (misalnya industri di utara dan pertanian di selatan).
Sementara itu, variabel IPM (0,275; p = 0,032) dan Rasio Gini (0,299; p = 0,036) sama-sama menunjukkan autokorelasi spasial positif yang signifikan. Ini menandakan bahwa tingkat pembangunan manusia dan ketimpangan pendapatan di suatu kabupaten/kota cenderung menyerupai wilayah yang berdekatan, yang bisa disebabkan oleh kesamaan dalam akses pendidikan, kesehatan, serta distribusi sumber daya ekonomi antarwilayah.
Secara keseluruhan, hasil ini memperkuat bukti adanya dependensi spasial antarwilayah di Jawa Barat, terutama pada variabel kemiskinan, TPT, IPM, dan Gini Ratio, sehingga analisis spasial lebih lanjut sangat diperlukan untuk memodelkan hubungan antarvariabel secara lebih akurat.
4.3.2 Geary’s C
| Variabel | Gearys_C |
|---|---|
| Kemiskinan | 0.410 |
| TPT | 0.465 |
| PDRB | 0.660 |
| IPM | 0.607 |
| GINI | 0.504 |
Dari hasil tersebut, variabel kemiskinan (0,410) dan TPT (0,465) memiliki nilai Geary’s C yang cukup jauh di bawah 1, menunjukkan autokorelasi spasial positif yang kuat. Artinya, wilayah dengan tingkat kemiskinan atau pengangguran tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki nilai tinggi.
Variabel IPM (0,607) dan Gini Ratio (0,504) juga menunjukkan autokorelasi spasial positif, meskipun dengan tingkat yang lebih moderat dibanding dua variabel sebelumnya. Hal ini menandakan bahwa pembangunan manusia dan ketimpangan pendapatan memiliki pola pengelompokan antarwilayah yang masih cukup nyata.
Sementara itu, nilai Geary’s C PDRB (0,660) mendekati 1, mengindikasikan autokorelasi spasial lemah atau hampir tidak ada keterkaitan spasial pada distribusi PDRB antar kabupaten/kota di Jawa Barat.
Secara keseluruhan, hasil Geary’s C ini konsisten dengan uji Moran’s I sebelumnya, di mana sebagian besar variabel (kecuali PDRB) menunjukkan adanya pola spasial positif, sehingga memperkuat indikasi adanya pengaruh spasial yang perlu diperhatikan dalam pemodelan ekonometrika spasial.
4.3.3 Local Moran’s I (LISA)
Hasil analisis Local Moran’s I (LISA) menunjukkan adanya pola pengelompokan spasial pada variabel sosial-ekonomi di Jawa Barat.
Untuk kemiskinan, terdapat klaster High-High di wilayah Garut dan Tasikmalaya, yang berarti daerah dengan kemiskinan tinggi berdekatan dengan daerah berkemiskinan tinggi lainnya. Sementara itu, klaster Low-Low muncul di wilayah utara, seperti Kota Bandung dan sekitarnya, yang menunjukkan daerah dengan kemiskinan rendah cenderung berdekatan dengan wilayah serupa.
TPT, IPM, dan Rasio Gini menunjukkan pola Low-Low, menandakan adanya pengelompokan wilayah dengan pengangguran, pembangunan manusia, dan ketimpangan pendapatan yang relatif rendah. Sedangkan PDRB per kapita tidak menunjukkan klaster signifikan, sehingga distribusinya bersifat acak.
Secara keseluruhan, hasil LISA mengonfirmasi adanya pola spasial lokal pada sebagian besar variabel, terutama kemiskinan, yang memperkuat adanya keterkaitan antarwilayah di Jawa Barat.
4.3.4 Getis–Ord 𝐺∗𝑖
Hasil analisis Getis-Ord G* menunjukkan bahwa kemiskinan membentuk klaster hot spot di wilayah Bogor dan sekitarnya (high-high) serta cold spot di bagian barat Jawa Barat (low-low). Variabel TPT, IPM, dan GINI menunjukkan pola cold spot di wilayah selatan yang menandakan pengangguran, pembangunan manusia, dan ketimpangan yang relatif rendah secara berkelompok. Sementara itu, PDRB tidak menunjukkan pola spasial yang signifikan, menandakan distribusinya relatif acak di seluruh wilayah Jawa Barat.
4.4 Uji Lagrange Multiplier
| Statistik_Uji | Value | p | Ket |
|---|---|---|---|
| LM Error (RSerr) | 3.4944 | 0.0616 | not sig. |
| LM Lag (RSlag) | 8.8489 | 0.0029 | sig. |
| Robust LM Error (AdjRSerr) | 1.0849 | 0.2976 | not sig. |
| Robust LM Lag (AdjRSlag) | 6.4394 | 0.0112 | sig. |
| LM SARMA (Gabungan) | 9.9338 | 0.0070 | sig. |
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, model spasial yang paling tepat digunakan adalah model (SAR). Hasil ini juga diperkuat oleh nilai LM SARMA (p = 0.007) yang signifikan, menegaskan adanya dependensi spasial secara keseluruhan dalam model.
4.5 Perbandingan Model
| Model | LogLik | AIC | Ket |
|---|---|---|---|
| SAR (Spatial Autoregressive) | -43.976 | 101.953 | Memasukkan lag pada variabel dependen; AIC terendah. Menunjukkan peningkatan kecocokan yang signifikan. |
| SDM (Spatial Durbin Model) | -40.163 | 102.327 | Memasukkan lag pada dependen dan independen; AIC rendah dan residual tidak berautokorelasi. |
| GNS (General Nesting Spatial) | -39.951 | 103.903 | Model paling kompleks; meski log-likelihood tertinggi, AIC lebih besar karena banyak parameter. |
| SAC/SARAR | -43.449 | 102.898 | Menggabungkan lag dan error spasial; AIC sedikit lebih tinggi dari SDM. |
Mempertimbangkan hasil uji LM, nilai Log-Likelihood, AIC, dan dependensi spasial pada setiap variabel, akan dibandingkan model Spatial Durbin dan model Spatial Autoregressive (SAR).
4.6 Model SDM
Berdasarkan hasil pengujian autokorelasi spasial, diketahui bahwa terdapat dependensi spasial, baik pada variabel X maupun variabel Y, maka model SDM cocok digunakan karena dapat mengakomodasi hal tersebut. Dengan estimasi parameter model SDM menggunakan MLE, kemiskinan di Jawa Barat dapat dimodelkan sebagai:
\[ \text{Kemiskinan}_i = 0.55 \sum_{j} w_{ij} y_j + 39.14^{***} + 0.26(\text{TPT})_i - 0.0033(\text{PDRB})_i - 0.42(\text{IPM})_i^{***} + 18.98(\text{GINI})_i^{***} \\ \quad + 0.46 \sum_{j} w_{ij} (\text{TPT})_j + 0.01 \sum_{j} w_{ij} (\text{PDRB})_j - 0.06 \sum_{j} w_{ij} (\text{IPM})_j - 13.86 \sum_{j} w_{ij} (\text{GINI})_j \]
| Statistik_Uji | rho | LR | p_value | Ket |
|---|---|---|---|---|
| Likelihood Ratio | 0.555 | 8.029 | 0.0046 | sig. |
Signifikansi nilai menunjukkan bahwa kemiskinan di suatu wilayah dipengaruhi oleh kemiskinan di wilayah sekitarnya. Jika kemiskinan di suatu kabupaten meningkat, maka kemiskinan di kabupaten tetangga cenderung ikut meningkat. Sama halnya dengan nilai variabel lain (TPT, PDRB, IPM, dan Gini Ratio).
| Variabel | Jenis_Impact | Z | P_value | Ket |
|---|---|---|---|---|
| IPM | Direct | -4.626 | 0.0000 | sig |
| Gini Ratio | Direct | 1.819 | 0.0689 | sig (=0.01) |
IPM memiliki dampak langsung negatif dan signifikan terhadap kemiskinan, yang berarti peningkatan pembangunan manusia berperan penting dalam menekan kemiskinan secara lokal, namun tidak menimbulkan efek spillover ke daerah sekitar. Sedangkan Gini ratio signifikan (dengan α = 0.01) mempengaruhi Kemiskinan secara langsung.
4.7 Uji Asumsi Model SDM
Berdasarkan uji asumsi, diketahui bahwa model SDM memenuhi semua uji asumsi. Hasil uji dapat dilihat pada tabel berikut:
| Uji | Statistik | P_value | Ket |
|---|---|---|---|
| Residual Autokorelasi | 0.5686 | 0.2848 | tidak terdapat autokorelasi spasial |
| Normalitas | 0.9845 | 0.9466 | berdistribusi normal |
| Homoskedastisitas | 9.3304 | 0.3152 | tidak terdapat heteroskedastisitas |
4.8 Model SAR
Berdasarkan hasil perbandingan model, model Spatial Autoregressive (SAR) dipilih menjadi model terbaik yang dapat menjelaskan adanya spatial dependence antar kabupaten/kota dan model ini juga mampu menjelaskan faktor-faktor yang mempengaruhi persentase kemiskinan di Jawa Barat. Dengan estimasi model SAR menggunakan MLE, kemiskinan di Jawa Barat dapat dimodelkan sebagai:
\[ y_i = \sum_{j=1}^{n} w_{ij} y_j + 28.449^{***} + 0.1132(\text{TPT})_i - 0.0087(\text{PDRB})_i - 0.4132(\text{IPM})_i^{***} + 12.705(\text{GINI})_i^{*} \]
| Statistik_Uji | rho | LR | p_value | Ket |
|---|---|---|---|---|
| Likelihood Ratio | 0.6527 | 12.428 | 4e-04 | sig. |
Signifikansi nilai menunjukkan bahwa kemiskinan di suatu wilayah dipengaruhi oleh kemiskinan di wilayah sekitarnya. Jika kemiskinan di suatu kabupaten meningkat, maka kemiskinan di kabupaten tetangga cenderung ikut meningkat.
| Variabel | Jenis_Impact | Koefisien | P_value | Ket |
|---|---|---|---|---|
| IPM | Direct | -0.4816 | 0.000 | sig |
| Indirect | -0.7083 | 0.504 | not.sig | |
| Total | -1.1899 | 0.339 | not.sig |
Keterangan:
Direct effect = pengaruh langsung variabel terhadap kemiskinan di wilayahnya sendiri.
Indirect effect = pengaruh tidak langsung (spillover) terhadap wilayah tetangga.
Total effect = gabungan dari efek langsung dan tidak langsung.
IPM memiliki dampak langsung negatif dan signifikan terhadap kemiskinan, yang berarti peningkatan pembangunan manusia berperan penting dalam menekan kemiskinan secara lokal, namun tidak menimbulkan efek ke daerah sekitar.
4.9 Uji Asumsi Model SAR
| Uji | Statistik | P_value | Ket |
|---|---|---|---|
| Residual Autokorelasi | -0.1110 | 0.7077 | not sig. |
| Normalitas | 0.9552 | 0.2865 | sig. |
| Homoskedastisitas | 3.0169 | 0.5550 | not sig. |
Asumsi linearitas pada model SAR dilakukan secara visual melalui pola residual. Berdasarkan plot residual terhadap nilai prediksi yang acak di sekitar nol, sehingga hubungan antara variabel dalam model SAR dapat dikatakan linier.
4.10 Uji Sensitivitas Bobot
| Bobot | rho | AIC | LogLik |
|---|---|---|---|
| Queen | 0.6527 | 101.9525 | -43.9762 |
| Rook | 0.6527 | 101.9525 | -43.9762 |
| KNN-4 | 0.4871 | 105.5565 | -45.7782 |
Uji sensivitas bobot spasial dilakukan untuk melihat stabilitas dari model SAR sebelumnya. Berdasarkan tabel bobot di atas, diketahui bahwa Queen dan Rook memberikan nilai , 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.11 Interpretasi Hasil
Berdasarkan peta prediksi kemiskinan (SDM dan SAR), terlihat bahwa kedua model sama-sama menunjukkan konsentrasi kemiskinan yang relatif tinggi di wilayah Bogor, Sukabumi, dan Cianjur bagian utara, serta Subang dan Indramayu. Sementara itu, wilayah dengan tingkat kemiskinan yang lebih rendah tampak di Bandung Raya (Kota Bandung, Kota Cimahi, Kabupaten Bandung, dan Kabupaten Bandung Barat) serta wilayah timur seperti Sumedang, Majalengka, dan Kuningan. Jika dibandingkan, model SAR (Spatial Autoregressive) memberikan hasil prediksi yang lebih halus dan konsisten di seluruh kabupaten/kota, sedangkan model SDM (Spatial Durbin Model) masih menunjukkan variasi ekstrem antarwilayah, misalnya nilai prediksi yang sangat tinggi di sebagian kecil Kabupaten Bogor dan terlalu rendah di wilayah selatan. Hal ini menunjukkan bahwa model SAR lebih mampu menangkap keterkaitan spasial antarwilayah secara proporsional tanpa menghasilkan prediksi yang berlebihan. | Hasil ini semakin diperkuat oleh peta residual standar. Pada model SDM, masih terlihat adanya pola spasial pada residual, terutama di wilayah Bekasi dan Bogor bagian utara yang cenderung underfit (berwarna merah), serta Tasikmalaya dan Ciamis yang menunjukkan overfit (berwarna biru). Ini menandakan bahwa model SDM belum sepenuhnya mengakomodasi pengaruh spasial antarwilayah tersebut. Sebaliknya, pada model SAR, distribusi residual tampak lebih acak dan intensitas warna lebih lembut di hampir seluruh kabupaten/kota Jawa Barat, termasuk di Cirebon, Garut, dan Sukabumi. Kondisi ini menandakan bahwa model SAR lebih baik dalam menjelaskan variasi spasial kemiskinan dan menghasilkan prediksi yang lebih stabil dan akurat. Dengan demikian, berdasarkan hasil visualisasi residual dan peta prediksi, model SAR dapat dinyatakan sebagai model terbaik untuk menggambarkan pola kemiskinan di Provinsi Jawa Barat.
BAB V. KESIMPULAN DAN SARAN
Berdasarkan hasil perbandingan model, model Spatial Autoregressive (SAR) terpilih sebagai model terbaik dengan nilai AIC terendah (101.95) dan hasil prediksi yang baik. Model SAR untuk kemiskinan secara matematis dapat dituliskan sebagai:
\[ y_i = \sum_{j=1}^{n} w_{ij} y_j + 28.449^{***} + 0.1132(\text{TPT})_i - 0.0087(\text{PDRB})_i - 0.4132(\text{IPM})_i^{***} + 12.705(\text{GINI})_i^{*} \]
Model ini paling sesuai untuk menganalisis faktor ekonomi dan sosial terhadap kemiskinan di Jawa Barat tahun 2024. Hasil menunjukkan bahwa IPM berpengaruh negatif dan signifikan terhadap kemiskinan, Gini Ratio signifikan pada taraf signifikansi 0.01, sedangkan TPT, PDRB tidak signifikan. Pola spasial juga terlihat jelas, di mana wilayah utara dan barat laut memiliki tingkat kemiskinan lebih tinggi dibanding wilayah tengah dan selatan.
Berdasarkan hasil analisis yang sudah dilakukan, dapat disimpulkan bahwa peningkatan IPM perlu menjadi fokus utama kebijakan pengentasan kemiskinan, dimana IPM ini mencakup tiga dimensi utama yaitu Pendidikan, Kesehatan, dan Standar Hidup yang Layak. Selain itu, fokus berdasarkan variabel Gini, pemerintah diharapkan juga bisa mengentaskan ketimpangan ekonomi (sehingga jarak antara yang kaya dan miskin tidak terlalu jauh).
Untuk penelitian berikutnya, disarankan untuk menambahkan variabel lain pada penelitian selanjutnya agar hasil analisis lebih komprehensif, dalam hal ini misalnya variabel akses terhadap air bersih, akses listrik, perumahan, dan lainnya.
DAFTAR PUSTAKA
PSKK UGM. (2018, 8 Juni). Upaya penanggulangan kemiskinan dari masa ke masa. Diakses darihttps://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.
Prasetya, I. P. G. I. B., Baharuddin, B., & Wibawa, G. N. A. (2024). Pemodelan Regresi Spasial untuk Menentukan Faktor-Faktor yang Berpengaruh terhadap Tingkat Kriminalitas di Provinsi Bali dan Jawa Timur. Jurnal Syntax Admiration, 5(6), 2033-2046.
Astutik, D. (2020). Analisis Faktor yang Mempengaruhi Kemiskinan di Jawa Timur (Pendekatan Spasial). Jurnal Ilmiah Mahasiswa FEB, 9(1).
Putri, A. A., & Sanusi, W. (2018). Model Regresi Spasial dan Aplikasinya pada Kasus Tingkat Kemiskinan Kabupaten Soppeng. Indonesian Journal of Fundamental Sciences Vol, 4(2).
LAMPIRAN ANALISIS
library(sf)## 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 TRUElibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(stringr)## Warning: package 'stringr' was built under R version 4.4.3library(spdep)## 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')`library(spatialreg)## 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.ZeroPolicyOptionlibrary(tmap) ## Warning: package 'tmap' was built under R version 4.4.3library(broom)## Warning: package 'broom' was built under R version 4.4.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3library(psych) ## Warning: package 'psych' was built under R version 4.4.3##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(viridis) ## Warning: package 'viridis' was built under R version 4.4.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.4.3library(scales)## 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, rescalelibrary(tibble)
library(purrr)## Warning: package 'purrr' was built under R version 4.4.3##
## Attaching package: 'purrr'## The following object is masked from 'package:scales':
##
## discardlibrary(forcats) ## Warning: package 'forcats' was built under R version 4.4.3library(htmltools)## Warning: package 'htmltools' was built under R version 4.4.3library(leaflet)## Warning: package 'leaflet' was built under R version 4.4.3# Import shapefile
setwd("C:/Users/Owner/Downloads/geoBoundaries-IDN-ADM2-all")
indo_adm2 <- st_read("geoBoundaries-IDN-ADM2.shp")## Reading layer `geoBoundaries-IDN-ADM2' from data source
## `C:\Users\Owner\Downloads\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 84kab_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("C:/Users/Owner/Downloads/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)
# Histogram Setiap Variabel ----
# Histogram KEMISKINAN
ggplot(data_spasial, aes(x = KEMISKINAN)) +
geom_histogram(
binwidth = 1, # atur lebar tiap batang agar batas jelas
fill = "skyblue",
color = "black",
boundary = 0, # memastikan bar mulai dari angka bulat
closed = "left" # agar sisi kiri dan kanan bar tertutup rapat
) +
scale_x_continuous(
name = "Kemiskinan (%)",
breaks = seq(floor(min(data_spasial$KEMISKINAN)), ceiling(max(data_spasial$KEMISKINAN)), by = 1)
) +
scale_y_continuous(name = "Frekuensi") +
labs(title = "Distribusi Tingkat Kemiskinan") +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.minor = element_blank()
)
# 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 GINIplots[[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# 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.01909999moran.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 = FALSE)## 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):
## Direct:
##
## Iterations = 1:996
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 996
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## TPT dy/dx 0.197685 0.238979 0.0075723 0.0075723
## PDRB dy/dx -0.001126 0.009938 0.0003149 0.0003149
## IPM dy/dx -0.494379 0.099348 0.0031480 0.0031480
## GINI dy/dx 18.715194 11.851586 0.3755319 0.3755319
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## TPT dy/dx -0.24152 0.064503 0.187069 0.333676 0.6056
## PDRB dy/dx -0.02078 -0.007163 -0.001676 0.004569 0.0187
## IPM dy/dx -0.69410 -0.548654 -0.489306 -0.433094 -0.3240
## GINI dy/dx 0.43026 12.454428 18.413817 25.042559 39.3885
##
## ========================================================
## Indirect:
##
## Iterations = 1:996
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 996
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## TPT dy/dx -0.6179 2.34520 0.074310 0.074310
## PDRB dy/dx 0.0249 0.08373 0.002653 0.002774
## IPM dy/dx -0.7752 0.90498 0.028675 0.025409
## GINI dy/dx -8.1287 140.80942 4.461717 4.461717
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## TPT dy/dx -2.6061 -1.1240 -0.62917 -0.21004 0.9927
## PDRB dy/dx -0.1016 -0.0167 0.01897 0.06029 0.1690
## IPM dy/dx -2.5393 -1.0601 -0.63136 -0.31984 0.2710
## GINI dy/dx -110.1849 -35.6953 -6.97074 22.39558 116.7446
##
## ========================================================
## Total:
##
## Iterations = 1:996
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 996
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## TPT dy/dx -0.42019 2.50394 0.079340 0.079340
## PDRB dy/dx 0.02378 0.09007 0.002854 0.002964
## IPM dy/dx -1.26953 0.95517 0.030266 0.027210
## GINI dy/dx 10.58653 149.66139 4.742203 4.742203
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## TPT dy/dx -2.7095 -0.99241 -0.46956 0.07018 1.4540
## PDRB dy/dx -0.1139 -0.02062 0.01918 0.06090 0.1863
## IPM dy/dx -3.0767 -1.56041 -1.10231 -0.79602 -0.1988
## GINI dy/dx -108.2936 -19.85880 11.57316 43.07807 148.5133
##
## ========================================================
## 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# 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)# BANDINGKAN 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# 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
# 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 geometriesnb_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.77824LAMPIRAN DASHBOARD
\[ # =========================================================# app.R — West Java Poverty Dashboard 2024 (VERSI 4.1 - FIX)# =========================================================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# -------------------------# Poppins font (UI & plot)try({ font_add_google("Poppins", "Poppins", regular.wt = 400, bold.wt = 700) showtext_auto(enable = TRUE)}, silent = TRUE)theme_set(theme_minimal(base_size = 12, base_family = "Poppins"))# -------------------------# 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")stats_limits <- function(x) { x <- x[is.finite(x)] mu <- mean(x) sig <- sd(x) lo <- mu - 3*sig hi <- mu + 3*sig lo <- max(min(x), lo) hi <- min(max(x), hi) c(lo, hi)}# Weightsnb_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; 2*pmin(p_right, p_left) } else if ("Pr(z <= 0)" %in% names(lisa)) { p_left <- lisa[,"Pr(z <= 0)"]; p_right <- 1 - p_left; 2*pmin(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 = "#b05d39", txt_color = "#3e2723") { 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:Poppins,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 = "logo_unpad.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("About", tabName = "about", icon = icon("info-circle")) ) ), dashboardBody( tags$head( tags$style(HTML(" /* Import Poppins Font */ @import url('https://fonts.googleapis.com/css2?family=Poppins:wght@400;600;700&display=swap'); :root{ --bg: #fdfaf6; /* Light Cream */ --header: #5d4037; /* Dark Warm Brown */ --sidebar: #5d4037; /* COKLAT TERANG (YANG DIGANTI) */ --sidebar-hover: #5d4037; /* Coklat Lebih Terang - Hover (YANG DIGANTI) */ --box-primary: #b05d39; /* Earthy Rust/Terracotta */ --text: #3e2723; /* Darkest Brown */ --text-inv: #fdfaf6; /* Light Cream (Teks sidebar tetap terang) */ --ink: #3e2723; /* Darkest Brown */ } /* Background dengan highlight lembut yang bergerak */ body, .content-wrapper, .right-side { background-color: var(--bg); font-family:'Poppins', 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(93,64,55,0.06), transparent 60%), radial-gradient(700px 450px at 80% 20%, rgba(176,93,57,0.05), transparent 60%), radial-gradient(800px 500px at 30% 85%, rgba(78,52,46,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(62,39,35,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); } /* PENYESUAIAN CSS SIDEBAR (SESUAI WARNA BARU) */ .skin-blue .main-sidebar { background-color: var(--sidebar); } .skin-blue .sidebar a { color: var(--text-inv); transition: color .2s ease; } /* Teks pakai variabel */ .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; /* Teks tetap putih/terang saat hover */ 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; /* Indikator putih cocok di coklat terang */ 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: none; /* Hapus semua border bawaan */ 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); } /* Hanya tambahkan border-top ke box non-solid */ .box:not(.box-solid) { border-top: 3px solid var(--box-primary); } .box:hover{ transform: translateY(-2px); box-shadow: 0 14px 30px rgba(0,0,0,0.12); } /* Pastikan header solid box masih memiliki border radius yang pas */ .box.box-solid.box-primary>.box-header { color:#fff; background: var(--box-primary); position: relative; overflow: hidden; /* Pastikan radius atas sesuai dengan box */ border-top-left-radius: 10px; border-top-right-radius: 10px; } /* 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(176,93,57,0.08); } /* Leaflet: kontrol & tooltip dengan efek halus */ .leaflet-container { font-family: 'Poppins', 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(62,39,35,0.15); color: var(--ink); transition: opacity .2s ease; } /* Scrollbar halus */ ::-webkit-scrollbar{ width: 10px; height:10px; } ::-webkit-scrollbar-thumb{ background: rgba(62,39,35,0.2); border-radius: 8px; transition: background .2s ease; } ::-webkit-scrollbar-thumb:hover{ background: rgba(62,39,35,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 ini TIDAK solid, jadi akan punya border-top box(width = 3, title = "Download Data", status = "primary", solidHeader = FALSE, 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 ini solid, jadi TIDAK akan punya border-top 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 = 10, value = 7, step = 1), 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( "<b>Catatan:</b> R² pada OLS adalah koefisien determinasi konvensional dari regresi linear. Untuk model spasial, nilai R² yang ditampilkan merupakan <i>Pseudo-R²</i> 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 = "About", status = "primary", solidHeader = TRUE, HTML(" <h4><b>Judul Analisis</b></h4> <p><b>PEMILIHAN MODEL TERBAIK DALAM MENGANALISIS FAKTOR EKONOMI DAN SOSIAL TERHADAP TINGKAT KEMISKINAN DI JAWA BARAT TAHUN 2024</b></p> <h4><b>Dibuat pada</b></h4> <p>17 Oktober 2025</p> <h4><b>Dibuat oleh</b></h4> <p>Ara Figura Bilbina — 140610230004</p> <p><b>Mata Kuliah:</b> Spatial Statistics</p> <p><b>Dosen Pengampu:</b> Dr. I Gede Nyoman Mindra Jaya, S. Si., M.Si.</p> <p><b>Instansi:</b> Program Studi Statistika, FMIPA, Universitas Padjadjaran.</p> <p>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.</p> ") ) ) ) ) ))# =========================================================# 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 (Earth Tone Cerah - Opsi 1: Terracotta & Sage) col_hist <- "#D9804C" # Terracotta cerah col_box <- "#9DB589" # Sage green col_outline<- "#3e2723" # Cokelat Tua (kontras) col_point <- "#9DB589" # Sage green col_smooth <- "#f2a900" # Amber output$hist_plot <- renderPlot({ v <- req(input$var_desc) ggplot(jabar, aes(x = .data[[v]])) + geom_histogram(bins = input$bins, color = col_outline, fill = col_hist, closed = "right") + scale_x_continuous(expand = expansion(mult = 0.02)) + labs(title = paste("Histogram", var_labels[[v]]), x = paste0(var_labels[[v]]), y = "Frekuensi") + theme_minimal(base_size = 12, base_family = "Poppins") }) output$box_plot <- renderPlot({ v <- req(input$var_desc) ggplot(jabar, aes(y = .data[[v]])) + geom_boxplot(fill = col_box, color = col_outline) + labs(title = paste("Boxplot", var_labels[[v]]), y = var_labels[[v]], x = "") + theme_minimal(base_size = 12, base_family = "Poppins") }) 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 = "Poppins") }) 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]] lims <- stats_limits(xv) brks <- pretty(lims, n = isolate(input$bins)) p1 <- ggplot(jabar, aes(x = .data[[v]])) + geom_histogram(breaks = brks, color = col_outline, fill = col_hist, closed = "right") + scale_x_continuous(limits = lims, breaks = brks, expand = expansion(mult = 0.02)) + labs(title = paste("Histogram", var_labels[[v]]), x = paste0(var_labels[[v]]), y = "Frekuensi") + theme_minimal(base_size = 12, base_family = "Poppins") p2 <- ggplot(jabar, aes(y = .data[[v]])) + geom_boxplot(fill = col_box, color = col_outline) + labs(title = paste("Boxplot", var_labels[[v]]), y = var_labels[[v]], x = "") + theme_minimal(base_size = 12, base_family = "Poppins") 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 = "Poppins") 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( "<b>%s</b><br/>%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 = "Poppins") + 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("<b>%s</b><br/>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("<b>%s</b><br/>Kategori: %s<br/>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("<b>%s</b><br/>%s: %s<br/>Moran's I: %s<br/>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 # ==================================================================== # === PERBAIKAN BUG ADA DI SINI === # ==================================================================== res <- safe_geary_both_raw(jabar[[v]], nb_queen) # MENGGUNAKAN NAMA YANG BENAR: "Geary C statistic" (tanpa 'R' dan tanpa apostrof) Geary_C <- round(unname(res$gt$estimate[["Geary C statistic"]]), 3) # Baris 'if (is.na...)' tidak diperlukan lagi karena nama di atas sudah benar # 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("<b>%s</b><br/>%s: %s<br/>Geary's C: %s<br/>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 = "Poppins") + 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 = "Poppins") + 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 = "Poppins") + 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 = "Poppins") + 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 = "#fdfaf6"), `Log-Likelihood_num` = `Log-Likelihood`, `AIC (bar)` = vapply(AIC, bar_cell, character(1), minv = rng_aic[1], maxv = rng_aic[2], color = "#b05d39", txt_color = "#fdfaf6"), AIC_num = AIC, `R² (bar)` = vapply(`R²`, bar_cell, character(1), minv = rng_r2[1], maxv = rng_r2[2], color = "#f2a900", txt_color = "#3e2723"), 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( "<b>%s</b><br/>Prediksi: %.3f<br/>Aktual: %.3f<br/>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 = "Poppins") + 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 = "#b05d39") + labs(title = "AIC per Bobot W") + theme_minimal(base_size = 12, base_family = "Poppins") 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 = "Poppins") 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 = "#b05d39") + labs(title = "AIC per Bobot W") + theme_minimal(base_size = 12, base_family = "Poppins") 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 = "Poppins") g <- gridExtra::arrangeGrob(pAIC, pRho, ncol = 2) ggsave(file, g, width = 10, height = 5.5, dpi = 300) } )}# =========================================================# RUN# =========================================================shinyApp(ui, server) \]