Laporan Ujian Akhir Semester

ANALISIS SPASIAL PENEMUAN KASUS BARU KUSTA DI JAWA BARAT 2024: INTEGRASI INTERPOLASI INVERSE DISTANCE WEIGHTING, SPATIAL AUTOREGRESSIVE MODEL, DAN GEOGRAPHICALLY WEIGHTED REGRESSION

Ditujukan Untuk Memenuhi Nilai Mata Kuliah Spatial Statistics TA. 2025/2026

Dosen Pengampu :

Dr. I Gede Nyoman Mindra Jaya, S. Si., M.Si.

Disusun Oleh :

Humaira Zeanova           140610230030

PROGRAM STUDI S1 STATISTIKA

FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM

UNIVERSITAS PADJADJARAN

JATINANGOR

2025

DASHBOARD

Dashboard analisis dapat diakses melalui link berikut: https://spatial-leprosy.streamlit.app/

dengan script yang dapat dilihat melalui: https://github.com/zeanovaa/spatial-leprosy

ataupun pada lampiran laporan ini.

VIDEO PEMBAHASAN

Video pembahasan dapat diakses melalui link berikut: https://youtu.be/DC4-1tEpjas?si=C1LxaxUmbR7F2FLB

Link alternatif: https://youtu.be/DC4-1tEpjas

Judul Video: Analisis Spasial Penemuan Kasus Baru Kusta Di Jawa Barat 2024 - Humaira Zeanova

Nama Channel: Humaira Zeanova

ABSTRAK

Penelitian ini bertujuan menganalisis pola distribusi spasial dan faktor penentu penemuan kasus baru kusta di Provinsi Jawa Barat tahun 2024 dengan mengintegrasikan metode interpolasi spasial dan regresi spasial. Data kabupaten/kota dianalisis menggunakan Inverse Distance Weighting (IDW) untuk menangani nilai hilang, regresi Ordinary Least Squares (OLS) sebagai model awal, Spatial Autoregressive Model (SAR) untuk menangkap ketergantungan spasial global, serta Geographically Weighted Regression (GWR) untuk mengidentifikasi heterogenitas spasial. Hasil analisis menunjukkan adanya autokorelasi spasial positif yang signifikan pada kasus baru kusta (Moran’s I = 0,402; p-value < 0,01). Model OLS memiliki performa rendah dengan R² sebesar 10,47% dan AIC 151,02. Model SAR memberikan peningkatan signifikan dengan parameter spasial ρ = 0,643 (p-value = 0,0038), R² sebesar 44,44%, dan AIC 144,64, namun seluruh variabel penjelas tidak signifikan secara global. Model GWR menunjukkan kinerja terbaik dengan nilai AIC terendah sebesar 122,80 dan quasi-global R² sebesar 69,42%, serta MAE 1,318 dan RMSE 1,926. Hasil GWR mengungkap heterogenitas spasial yang kuat, di mana pengaruh kemiskinan signifikan secara lokal terutama di wilayah timur Jawa Barat. Temuan ini menegaskan pentingnya pendekatan spasial berbasis wilayah dalam perumusan kebijakan pengendalian kusta.

BAB I PENDAHULUAN

1.1 Latar Belakang Masalah

Kusta masih menjadi masalah kesehatan masyarakat yang serius di Indonesia. Bahkan, Indonesia menempati peringkat tiga dunia dalam jumlah penemuan kasus baru kusta, yaitu sebesar 12798 kasus baru, pada tahun 2023. Sebagai penyakit infeksi kronis, kusta berisiko menimbulkan kecacatan permanen serta dampak sosial dan ekonomi yang berkepanjangan apabila tidak ditangani secara dini (Kementerian Kesehatan, 2025). Kondisi ini menunjukkan bahwa upaya eliminasi kusta memerlukan pendekatan yang lebih terarah dan berbasis wilayah, bukan sekadar analisis agregat.

Provinsi Jawa Barat merupakan salah satu provinsi dengan kontribusi kasus kusta yang relatif tinggi secara nasional. Kondisi ini tidak terlepas dari karakteristik wilayah Jawa Barat yang heterogen, baik dari sisi kepadatan penduduk, kondisi sosial ekonomi, maupun perilaku hidup bersih dan sehat. Data penemuan kasus baru kusta tahun 2024 menunjukkan adanya variasi yang cukup besar antar kabupaten/kota di Jawa Barat, yang mengindikasikan bahwa risiko penularan dan penemuan kasus tidak tersebar secara merata.

Gambar 1. Peta Penyebaran Penemuan Kasus Baru Kusta Jawa Barat 2024

Hal tersebut terlihat pada Gambar 1, yang menyajikan peta distribusi spasial penemuan kasus baru kusta di Jawa Barat tahun 2024. Pada peta tersebut tampak bahwa beberapa kabupaten/kota membentuk konsentrasi wilayah dengan nilai penemuan kasus yang relatif tinggi, sementara wilayah lainnya menunjukkan nilai yang jauh lebih rendah. Pola ini mengindikasikan bahwa distribusi penemuan kasus baru kusta bersifat tidak acak dan berpotensi dipengaruhi oleh interaksi antarwilayah yang berdekatan, sesuai dengan prinsip kedekatan geografis dalam analisis spasial (Anselin, 1988).

Keberadaan pola spasial tersebut menunjukkan bahwa analisis non-spasial yang mengasumsikan independensi antarwilayah berpotensi menghasilkan kesimpulan yang kurang tepat. Ketergantungan spasial antar kabupaten/kota dapat menyebabkan bias estimasi parameter apabila tidak dimodelkan secara eksplisit. Oleh karena itu, pendekatan ekonometrika spasial menjadi relevan untuk menangkap efek limpahan (spatial spillover) dalam penemuan kasus baru kusta (LeSage & Pace, 2009).

Selain ketergantungan spasial, hubungan antara penemuan kasus baru kusta dan faktor-faktor penentunya diduga tidak bersifat homogen di seluruh wilayah. Faktor kepadatan penduduk, perilaku hidup bersih dan sehat (PHBS), serta tingkat kemiskinan dapat memiliki pengaruh yang berbeda antar kabupaten/kota, tergantung pada kondisi lokal masing-masing wilayah. Untuk menangkap variasi lokal tersebut, diperlukan pendekatan Geographically Weighted Regression (GWR) yang memungkinkan koefisien regresi bervariasi secara spasial (Fotheringham et al., 2002).

Di sisi lain, data penemuan kasus baru kusta pada tingkat kabupaten/kota tahun 2024 masih menghadapi permasalahan nilai hilang (missing value) pada beberapa wilayah. Keterbatasan ini berpotensi menghambat proses pemetaan risiko dan analisis spasial lanjutan. Oleh karena itu, metode interpolasi spasial Inverse Distance Weighting (IDW) digunakan untuk mengestimasi nilai yang hilang berdasarkan kedekatan geografis antarwilayah. Metode ini dipilih karena bersifat sederhana, deterministik, dan sesuai untuk data agregat wilayah dengan jumlah titik pengamatan yang terbatas (Burrough & McDonnell, 1998).

Berdasarkan uraian tersebut, penelitian ini mengintegrasikan interpolasi IDW untuk melengkapi data, model ekonometrika spasial untuk menganalisis pengaruh global faktor penentu, serta GWR untuk mengkaji heterogenitas spasial pengaruh faktor-faktor tersebut terhadap penemuan kasus baru kusta di Jawa Barat tahun 2024.

1.2 Identifikasi Masalah

Berdasarkan latar belakang yang telah diuraikan, meskipun berbagai upaya pengendalian kusta telah dilakukan, penemuan kasus baru kusta di Jawa Barat tahun 2024 masih menunjukkan distribusi yang tidak merata secara geografis. Kondisi ini mengindikasikan adanya kesenjangan antara pendekatan analisis konvensional yang bersifat agregat dengan realitas spasial penularan kusta yang dipengaruhi oleh interaksi antarwilayah dan karakteristik lokal. Oleh karena itu, permasalahan penelitian ini dapat diidentifikasi sebagai berikut:

1. Distribusi penemuan kasus baru kusta di Jawa Barat tahun 2024 menunjukkan variasi antar kabupaten/kota yang cukup besar, sehingga pola risiko spasial belum dapat dijelaskan secara memadai hanya dengan statistik deskriptif atau analisis non-spasial.

2. Terdapat indikasi ketergantungan spasial antarwilayah dalam penemuan kasus baru kusta, yang berpotensi menyebabkan estimasi parameter menjadi bias apabila dianalisis menggunakan model regresi global tanpa mempertimbangkan efek spasial.

3. Hubungan antara penemuan kasus baru kusta dan faktor-faktor penentunya diduga tidak bersifat homogen di seluruh wilayah Jawa Barat, sehingga diperlukan pendekatan yang mampu menangkap heterogenitas spasial pengaruh variabel antar lokasi.

4. Informasi mengenai wilayah dengan tingkat risiko relatif tinggi masih belum terpetakan secara kontinu dalam bentuk permukaan spasial, sehingga penentuan prioritas intervensi berbasis wilayah belum optimal.

1.3. Rumusan Masalah

Berdasarkan teori autokorelasi spasial dan temuan empiris sebelumnya yang menunjukkan bahwa kasus kusta memiliki pola distribusi tidak acak serta dipengaruhi oleh ketergantungan dan heterogenitas spasial, maka rumusan masalah dalam penelitian ini adalah sebagai berikut:

1. Bagaimana pola distribusi spasial penemuan kasus baru kusta di Jawa Barat tahun 2024 berdasarkan eksplorasi peta dan interpolasi spasial

2. Apakah terdapat autokorelasi spasial dalam penemuan kasus baru kusta antar kabupaten/kota di Jawa Barat tahun 2024

3. Apakah model regresi spasial mampu menjelaskan penemuan kasus baru kusta dengan lebih baik dibandingkan model non-spasial

4. Apakah pengaruh faktor-faktor penentu penemuan kasus baru kusta menunjukkan variasi spasial antar wilayah di Jawa Barat

1.4 Tujuan Penelitian

Penelitian ini bertujuan untuk :

1. Menggambarkan pola distribusi spasial penemuan kasus baru kusta di Jawa Barat tahun 2024 melalui peta tematik dan interpolasi spasial

2. Mengidentifikasi keberadaan autokorelasi spasial penemuan kasus baru kusta menggunakan indeks autokorelasi spasial

3. Memodelkan penemuan kasus baru kusta menggunakan model regresi spasial dan membandingkannya dengan model non-spasial

4. Menganalisis variasi lokal pengaruh faktor-faktor penentu penemuan kasus baru kusta melalui pendekatan Geographically Weighted Regression

1.5 Batasan Penelitian

Penelitian ini memiliki beberapa batasan yang perlu diperhatikan dalam interpretasi hasil, yaitu sebagai berikut :

1. Analisis dilakukan menggunakan data cross section tahun 2024 sehingga hasil penelitian hanya menggambarkan kondisi spasial pada periode tersebut. Implikasinya, kesimpulan yang diperoleh tidak dapat digunakan untuk menilai tren penularan kusta dari waktu ke waktu.
2. Unit analisis regresi spasial dibatasi pada tingkat kabupaten/kota di Jawa Barat karena keterbatasan ketersediaan data pada tingkat kecamatan. Implikasinya, hasil estimasi regresi merefleksikan hubungan antarvariabel pada skala makro wilayah dan belum tentu sepenuhnya menggambarkan variasi kondisi pada tingkat lokal yang lebih kecil.
3. Interpolasi risiko kusta dilakukan menggunakan metode Inverse Distance Weighting (IDW), yang bersifat deterministik dan mengasumsikan pengaruh jarak terhadap nilai prediksi. Oleh karena itu, hasil interpolasi merupakan estimasi spasial dan tidak dapat disamakan dengan data observasi nyata di setiap lokasi.
4. Model regresi spasial yang digunakan dibatasi pada Spatial Econometrics Model dan Geographically Weighted Regression berdasarkan hasil diagnosa autokorelasi spasial sehingga masih ada kemungkinan spesifikasi model lain yang relevan belum dieksplorasi dalam penelitian ini.
5. Variabel penjelas yang digunakan hanya mencakup kepadatan penduduk, PHBS, dan tingkat kemiskinan, sehingga faktor lain seperti akses pelayanan kesehatan, mobilitas penduduk, dan karakteristik lingkungan belum tercakup, dan dapat memengaruhi distribusi penemuan kasus kusta secara tidak langsung.

Batasan-batasan tersebut menegaskan bahwa hasil penelitian ini perlu ditafsirkan secara kontekstual dan digunakan sebagai dasar perencanaan intervensi berbasis wilayah, bukan sebagai penilaian absolut terhadap risiko kusta di setiap kabupaten/kota.

BAB II TINJAUAN PUSTAKA

2.1 Spatial Dependence

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

2.2 Autokorelasi Spasial

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

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

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

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

2.3 Ordinary Least Squares

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

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

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

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

2.4 Uji Lagrange Multiplier

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

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

2.5 Spatial Autoregressive Model

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

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

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

2.6 Akaike Information Criterion

Menurut Suryaningrat, S. H. (2021) Akaike Information Criterion (AIC) merupakan metode analisis yang digunakan untuk memperoleh model yang terbaik dengan menggunakan estimasi maximum likelihood sebagai perhitungan yang sesuai. AIC mempertimbangkan keseimbangan antara kompleksitas model (jumlah parameter) dan kecocokan model terhadap data (likelihood). AIC dihitung dengan rumus:

\[ AIC = 2k - 2 \ln(\hat{L}) \]

di mana k adalah jumlah parameter dalam model dan \(\hat{L}\) adalah nilai maksimum dari fungsi likelihood. Model dengan nilai AIC yang lebih rendah dianggap lebih baik karena menunjukkan bahwa model tersebut memiliki kompleksitas yang wajar tetapi tetap mampu menjelaskan data dengan baik.

2.7 Inverse Distance Weighting (IDW)

Inverse Distance Weighting (IDW) merupakan metode interpolasi deterministik yang menentukan nilai pada lokasi yang tidak teramati berdasarkan rata-rata tertimbang dari titik-titik sampel di sekitarnya. Penggunaan IDW dalam analisis data spasial didasarkan pada asumsi bahwa nilai pada lokasi yang berdekatan memiliki kemiripan yang lebih besar dibandingkan lokasi yang jauh. Bobot (\(w\)) dalam IDW dihitung berdasarkan fungsi jarak (\(d\)) dengan parameter kuasa (power) \(p\) yang dioptimalkan untuk meminimalkan tingkat kesalahan estimasi (Shepard, 1968). Formulasi dasarnya adalah:

\[ \hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i Z(s_i)}{\sum_{i=1}^{n} w_i} ; w_i = \frac{1}{d_i^p} \]

Penerapan IDW mensyaratkan diagnosis terhadap karakteristik data untuk memastikan akurasi interpolasi. Pertama, data diasumsikan memiliki stasioneritas, di mana varians data cenderung konstan di berbagai wilayah, yang dalam syntax dideteksi melalui uji Levene dan variogram cloud. Kedua, IDW bersifat isotropik, yang berarti pengaruh jarak dianggap sama ke segala arah. Oleh karena itu, diagnosis anisotropi melalui directional variogram penting dilakukan untuk memastikan tidak ada arah tertentu yang memiliki korelasi lebih kuat. Ketiga, keberadaan tren spasial (baik linear maupun kuadratik) harus diidentifikasi karena IDW paling efektif digunakan pada data yang tidak memiliki tren global yang ekstrem (Lloyd, 2010).

2.8 Geographically Weighted Regression (GWR)

Geographically Weighted Regression (GWR) adalah pengembangan dari regresi linear yang digunakan untuk memodelkan hubungan antarvariabel yang bersifat lokal. Berbeda dengan OLS yang bersifat global, GWR mengakomodasi adanya heterogenitas spasial atau non-stasioneritas, di mana hubungan antara variabel independen terhadap variabel dependen berubah-ubah antar lokasi (Fotheringham et al., 2002). Model GWR dituliskan sebagai berikut:

\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) x_{ik} + \varepsilon_i \]

Dalam model ini, setiap lokasi pengamatan \((u_i, v_i)\) memiliki koefisien regresi sendiri. Penentuan bobot pada GWR dilakukan melalui fungsi kernel dengan parameter bandwidth. Berdasarkan syntax analisis, penggunaan adaptive bandwidth lebih disarankan jika distribusi titik observasi (pusat massa kabupaten/kota) tidak merata, sehingga radius pencarian tetangga akan menyesuaikan dengan kerapatan data di wilayah tersebut (Brunsdon et al., 1996). Validasi model GWR dilakukan dengan membandingkan nilai AIC-nya dengan model OLS serta melakukan uji signifikansi koefisien lokal menggunakan statistik \(t\) untuk memetakan wilayah mana saja yang memiliki pengaruh variabel signifikan secara spasial.

BAB III METODOLOGI PENELITIAN

3.1 Data

Data penelitian diperoleh dari Badan Pusat Statistik (BPS) terkait penemuan kasus baru kusta dan faktor-faktor determinan di Jawa Barat tahun 2024. Variabel yang digunakan dijelaskan pada Tabel 1. berikut :

Tabel 1. Definisi Operasional Variabel

Variabel	Satuan	Definisi Operasional
KUSTA	Jumlah kasus per 100.000 penduduk	Jumlah penemuan kasus baru kusta di masing-masing kabupaten/kota
KEPADATAN	Jiwa/km²	Jumlah penduduk per km² di tiap kabupaten/kota
PHBS	Persentase (%)	Persentase Rumah Tangga yang menerapkan Perilaku Hidup Bersih dan Sehat
MISKIN	Persentase (%)	Persentase penduduk miskin di tiap kabupaten/kota

3.2 Unit Spasial

Unit analisis dalam penelitian ini adalah kabupaten/kota di Provinsi Jawa Barat, sebanyak 27 wilayah administratif. Analisis dilakukan pada tingkat kabupaten/kota karena ketersediaan data variabel penjelas hanya pada skala ini.

3.3 Metode Analisis dan Alur Kerja Penelitian

Penelitian ini menggunakan kombinasi metode Inverse Distance Weighting (IDW), Spatial Autoregressive Model, dan Geographically Weighted Regression (GWR) untuk menganalisis distribusi dan faktor determinan penemuan kasus baru kusta. Langkah-langkah analisis dijelaskan sebagai berikut :

1. Analisis Deskriptif dan Eksplorasi Spasial

   Analisis diawali dengan statistik deskriptif dan pemetaan tematik untuk menggambarkan pola sebaran penemuan kasus baru kusta dan variabel penjelasnya. Tahap ini digunakan untuk memperoleh gambaran awal distribusi spasial, ketidakseragaman antarwilayah, serta indikasi awal adanya ketergantungan spasial.

2. Diagnosis Data Spasial

   Diagnosis data spasial dilakukan untuk mengevaluasi karakteristik statistik dan spasial data sebagai dasar pemilihan metode interpolasi. Tahapan ini meliputi uji normalitas, uji stasioneritas rata-rata dan varians, pengujian autokorelasi spasial, pemeriksaan anisotropi, serta evaluasi kerapatan dan sebaran data. Berdasarkan hasil diagnosis, metode Inverse Distance Weighting (IDW) dinilai paling sesuai dan dipilih sebagai metode interpolasi.

3. Interpolasi Spasial dengan Inverse Distance Weighting (IDW)

   Interpolasi spasial dilakukan menggunakan metode IDW untuk mengimputasi nilai penemuan kasus baru kusta yang hilang. Tahapan meliputi optimasi parameter power dan jumlah tetangga menggunakan pendekatan Leave-One-Out Cross Validation (LOOCV), dilanjutkan dengan proses prediksi, imputasi nilai hilang, dan visualisasi hasil interpolasi.

4. Ordinary Least Square

   Model regresi Ordinary Least Squares (OLS) dibangun sebagai model awal untuk menganalisis hubungan antara penemuan kasus baru kusta dan variabel penjelas. Selanjutnya dilakukan pemeriksaan asumsi klasik dan uji autokorelasi spasial pada residual model untuk mengidentifikasi keterbatasan model global nonspasial.

5. Analisis Autokorelasi Spasial

   Analisis autokorelasi spasial dilakukan secara global menggunakan indeks Moran’s I dan Geary’s C, serta secara lokal menggunakan Local Moran’s I (LISA) dan statistik Getis–Ord Gi*. Tahap ini bertujuan untuk mengidentifikasi pola klaster spasial dan memperkuat indikasi perlunya pendekatan pemodelan spasial.

6. Pemodelan Ekonometrika Spasial

   Pemodelan ekonometrika spasial diawali dengan uji Lagrange Multiplier untuk menentukan kandidat model yang sesuai. Berdasarkan hasil pengujian dan perbandingan nilai Akaike Information Criterion (AIC), Spatial Lag Model (SAR) dipilih sebagai model ekonometrika spasial utama. Selanjutnya dilakukan estimasi parameter, analisis direct dan indirect effects, pemeriksaan asumsi model, uji sensitivitas terhadap matriks bobot spasial, serta visualisasi hasil prediksi model.

7. Geographically Weighted Regression (GWR)

   Untuk mengidentifikasi heterogenitas spasial hubungan antarvariabel, dilakukan analisis Geographically Weighted Regression (GWR). Tahap ini diawali dengan perbandingan nilai AIC antara model OLS dan GWR, di mana nilai AIC yang lebih kecil pada GWR menunjukkan adanya variasi spasial lokal. Selanjutnya dibangun model GWR, dilakukan visualisasi koefisien lokal, serta pemeriksaan asumsi model.

Seluruh tahapan analisis tersebut disusun secara sistematis dan saling berurutan, mulai dari eksplorasi spasial hingga estimasi dan prediksi spasial.

BAB IV HASIL DAN PEMBAHASAN

4.1 Analisis Deskriptif

Analisis deskriptif dilakukan untuk memberikan gambaran awal mengenai karakteristik penemuan kasus baru kusta serta variabel-variabel penjelas yang digunakan dalam penelitian ini, yaitu kepadatan penduduk, persentase rumah tangga ber-PHBS, dan persentase penduduk miskin pada tingkat kabupaten/kota di Provinsi Jawa Barat tahun 2024. Ringkasan statistika deskriptif masing-masing variabel disajikan pada Tabel 2.

Tabel 2. Statistika Deskriptif

Variabel	n	Mean	SD	Varians	Median	Minimum	Maksimum	Skewness	Kurtosis
KUSTA	25	3.29	3.61	13.03110	1.22	0.16	14.98	1.43	1.86
KEPADATAN	27	3910.93	4668.12	21791310	1468.00	385.00	15176.00	1.22	-0.07
PHBS	27	67.87	10.03	100.66610	68.00	45.48	84.67	-0.51	-0.52
MISKIN	27	8.01	2.65	7.014903	8.41	2.34	11.93	-0.36	-0.97

Berdasarkan Tabel 2, penemuan kasus baru kusta per 100.000 penduduk memiliki nilai rata-rata sebesar 3,29 dengan simpangan baku 3,61. Nilai median yang relatif lebih rendah, yaitu 1,22, serta nilai skewness positif sebesar 1,43 menunjukkan bahwa distribusi kasus baru kusta cenderung right skewed. Hal ini mengindikasikan bahwa sebagian besar kabupaten/kota memiliki tingkat penemuan kasus yang rendah, sementara hanya beberapa wilayah yang memiliki nilai relatif tinggi. Rentang nilai yang cukup lebar, dari 0,16 hingga 14,98, semakin menegaskan adanya variasi yang cukup besar antarwilayah.

Variabel kepadatan penduduk menunjukkan rata-rata sebesar 3.910,93 jiwa per km² dengan simpangan baku yang sangat besar, yaitu 4.668,12. Nilai maksimum mencapai 15.176 jiwa per km², sedangkan nilai minimum hanya 385 jiwa per km². Kondisi ini mencerminkan heterogenitas karakteristik wilayah Jawa Barat, mulai dari daerah perkotaan dengan kepadatan sangat tinggi hingga wilayah perdesaan dengan kepadatan rendah. Skewness positif sebesar 1,22 menunjukkan bahwa sebagian besar wilayah memiliki kepadatan penduduk di bawah rata-rata, dengan beberapa wilayah ekstrem berpenduduk sangat padat.

Persentase rumah tangga ber-PHBS memiliki rata-rata sebesar 67,87% dengan simpangan baku 10,03 dan median 68,00. Nilai skewness yang sedikit negatif (-0,51) menunjukkan bahwa distribusi PHBS relatif simetris dengan kecenderungan sedikit lebih banyak wilayah yang memiliki nilai di atas rata-rata. Rentang nilai PHBS yang cukup lebar, yaitu antara 45,48% hingga 84,67%, menunjukkan adanya perbedaan tingkat perilaku hidup bersih dan sehat antar kabupaten/kota.

Sementara itu, persentase penduduk miskin memiliki nilai rata-rata sebesar 8,01% dengan simpangan baku 2,65. Nilai minimum tercatat sebesar 2,34% dan maksimum 11,93%. Skewness negatif sebesar -0,36 menunjukkan bahwa sebagian besar wilayah memiliki tingkat kemiskinan yang relatif berada di sekitar atau sedikit di atas rata-rata, dengan beberapa wilayah yang memiliki tingkat kemiskinan rendah.

Adapun penyebaran spasial dari masing-masing variabel dapat dilihat pada Gambar 2.

Gambar 2. Peta Penyebaran Variabel

Berdasarkan peta penemuan kasus baru kusta, terlihat bahwa nilai kasus tidak tersebar secara merata di seluruh wilayah Jawa Barat. Beberapa kabupaten/kota membentuk klaster dengan tingkat penemuan kasus relatif tinggi, sementara wilayah lain menunjukkan nilai yang jauh lebih rendah. Pola ini mengindikasikan adanya konsentrasi spasial kasus kusta dan menguatkan dugaan bahwa distribusi kasus bersifat tidak acak secara geografis. Selain itu, terlihat beberapa wilayah yang ditampilkan dengan warna putih, menunjukkan adanya missing value pada data penemuan kasus baru kusta. Keberadaan nilai hilang ini menandakan keterbatasan ketersediaan data pada wilayah tertentu dan berpotensi menghambat analisis spasial lanjutan apabila tidak ditangani secara khusus.

Peta kepadatan penduduk menunjukkan kontras yang jelas antara wilayah perkotaan dan perdesaan. Kabupaten/kota dengan karakteristik perkotaan tampak memiliki kepadatan penduduk yang jauh lebih tinggi dibandingkan wilayah lainnya. Perbedaan spasial ini berpotensi memengaruhi dinamika penularan penyakit menular, termasuk kusta, melalui intensitas interaksi penduduk.

Pada peta persentase rumah tangga ber-PHBS, terlihat variasi tingkat PHBS antarwilayah dengan pola yang tidak seragam. Beberapa wilayah menunjukkan capaian PHBS yang relatif tinggi, sementara wilayah lainnya masih berada pada tingkat menengah hingga rendah. Variasi ini mengindikasikan perbedaan kondisi perilaku kesehatan masyarakat yang berpotensi berpengaruh terhadap risiko penularan dan penemuan kasus kusta.

Sementara itu, peta persentase penduduk miskin menunjukkan adanya konsentrasi wilayah dengan tingkat kemiskinan relatif tinggi, khususnya pada beberapa kabupaten tertentu. Pola ini mengindikasikan bahwa faktor sosial ekonomi memiliki dimensi spasial yang kuat dan berpotensi berinteraksi dengan faktor lingkungan dan kesehatan dalam memengaruhi penemuan kasus kusta.

Secara keseluruhan, hasil analisis deskriptif dan peta sebaran menunjukkan adanya variasi dan pola spasial yang jelas pada variabel-variabel penelitian. Temuan ini menjadi dasar awal yang kuat untuk melakukan analisis autokorelasi spasial dan pemodelan regresi spasial pada tahap selanjutnya, guna menangkap ketergantungan antarwilayah serta heterogenitas pengaruh faktor-faktor penentu penemuan kasus baru kusta di Jawa Barat.

4.2 Diagnosis Data Interpolasi

Analisis statistika deskriptif menunjukkan adanya missing value pada data penemuan kasus baru kusta, sehingga diperlukan penanganan melalui proses interpolasi spasial. Sebelum metode interpolasi ditentukan, dilakukan serangkaian diagnosis data untuk mengidentifikasi karakteristik statistik dan spasial data penemuan kasus baru kusta. Hasil diagnosis tersebut digunakan sebagai dasar pemilihan metode interpolasi yang paling sesuai. Diagnosis ini mencakup uji normalitas, identifikasi tren spasial, autokorelasi spasial, uji stasioneritas dalam rata-rata maupun varians, memeriksa anisotropi, serta kerapatan data observasi.

4.2.1 Uji Normalitas

Uji normalitas dilakukan untuk mengevaluasi apakah data penemuan kasus baru kusta berdistribusi normal, yang penting dalam analisis variogram dan pemodelan spasial.

Hipotesis
H₀ : Data penemuan kasus baru kusta berdistribusi normal
H₁ : Data penemuan kasus baru kusta tidak berdistribusi normal

Kriteria Uji
Tolak H₀ jika p-value < α (0,05)

Hasil Uji

Diperoleh p-value = 0,0002

Berdasarkan hasil uji Shapiro–Wilk terhadap data asli penemuan kasus baru kusta, diperoleh nilai p-value kurang dari taraf signifikansi 0,05 sehingga H₀ ditolak. Dengan demikian, data asli penemuan kasus baru kusta tidak berdistribusi normal. Dengan kondisi tersebut, metode interpolasi yang dipertimbangkan perlu tidak bergantung pada asumsi kenormalan data. Metode deterministik seperti Inverse Distance Weighting (IDW) dan Spline dapat digunakan karena lebih menekankan pada kedekatan spasial antar titik pengamatan. Sebaliknya, metode geostatistik seperti Kriging memerlukan asumsi distribusi residual yang mendekati normal, sehingga kurang sesuai tanpa transformasi data. Oleh karena itu, metode interpolasi deterministik dinilai lebih relevan untuk data kasus kusta ini.

4.2.2 Uji Stasioneritas Rata-Rata (Tren Spasial)

Uji stasioneritas rata-rata bertujuan untuk mengidentifikasi apakah nilai harapan penemuan kasus baru kusta bersifat konstan di seluruh wilayah atau menunjukkan kecenderungan perubahan sistematis terhadap arah geografis tertentu. Identifikasi tren spasial dilakukan menggunakan model regresi polinomial dengan koordinat spasial (X dan Y) sebagai variabel penjelas, yaitu model tren orde pertama (linear) dan orde kedua (kuadratik). Hasil estimasi parameter untuk masing-masing model disajikan pada Tabel 3 dan Tabel 4.

Tabel 3. Hasil Estimasi Model Tren Spasial Orde Pertama

Variabel	Koefisien	Std. Error	p-value
Intersep	-287,87	82,97	0,002
X	3,20	0,81	< 0,001
Y	7,84	1,30	< 0,001

Model tren orde pertama menunjukkan bahwa koordinat X (arah Barat–Timur) dan Y (arah Utara–Selatan) berpengaruh signifikan terhadap penemuan kasus baru kusta. Nilai koefisien yang positif mengindikasikan adanya kecenderungan peningkatan nilai kusta seiring perubahan posisi geografis tertentu.

Tabel 4. Hasil Estimasi Model Tren Spasial Orde Kedua

Variabel	Koefisien	Std. Error	p-value
Intersep	26.392,59	12.283,66	0,045
X	-556,42	232,71	0,027
Y	-1.000,52	234,55	< 0,001
X²	2,96	1,11	0,015
Y²	15,48	3,32	< 0,001
X·Y	11,34	2,50	< 0,001

Pada model tren orde kedua, seluruh komponen polinomial terbukti signifikan secara statistik. Hal ini menunjukkan bahwa pola perubahan penemuan kasus kusta tidak hanya bersifat linier, tetapi juga dipengaruhi oleh interaksi dan efek nonlinier antar arah geografis.

Untuk menentukan model tren yang paling sesuai, dilakukan perbandingan kriteria informasi Akaike (AIC) sebagaimana disajikan pada Tabel 5.

Tabel 5. Perbandingan Kriteria AIC Model Tren Spasial

Model	AIC
Tren Orde Pertama	116,94
Tren Orde Kedua	102,29

Berdasarkan Tabel 5, model tren orde kedua memiliki nilai AIC yang lebih kecil dibandingkan model orde pertama. Temuan ini menunjukkan bahwa pola tren spasial bersifat nonlinier dan tidak dapat dijelaskan secara memadai oleh model linier sederhana. Dengan demikian, secara rata-rata, data penemuan kasus baru kusta tidak memenuhi asumsi stasioneritas mean.

Visualisasi tren spasial pada arah Barat–Timur dan Utara–Selatan ditampilkan pada Gambar 3. Pada kedua arah tersebut terlihat adanya kecenderungan peningkatan nilai kasus, terutama pada arah Utara–Selatan, yang semakin menguatkan indikasi adanya tren spasial.

Gambar 3. Tren Spasial Penemuan Kasus Baru Kusta Arah Barat–Timur dan Utara–Selatan

Selain itu, pola semivariogram eksperimental yang disajikan pada Gambar 4. menunjukkan bahwa nilai semivarian meningkat seiring bertambahnya jarak dan tidak mencapai titik sill atau kondisi mendatar. Pola semivariogram yang cenderung meningkat secara eksponensial ini merupakan indikasi kuat bahwa data mengandung tren dan tidak stasioner terhadap rata-rata.

Gambar 4. Semivariogram Eksperimental Penemuan Kasus Baru Kusta

Berdasarkan hasil estimasi model tren, perbandingan nilai AIC, visualisasi tren arah geografis, serta pola semivariogram, dapat disimpulkan bahwa data penemuan kasus baru kusta tidak stasioner secara rata-rata, sehingga asumsi stasioneritas mean tidak terpenuhi.

Metode interpolasi yang mengasumsikan mean konstan, seperti Ordinary Kriging, menjadi kurang sesuai tanpa proses detrending terlebih dahulu. Pada data yang tidak stasioner, metode interpolasi deterministik seperti Inverse Distance Weighting (IDW) lebih tepat digunakan karena tidak bergantung pada asumsi stasioneritas, atau alternatif lain adalah pendekatan kriging dengan pemodelan tren eksplisit. Namun, tanpa pemodelan tren yang kuat, metode deterministik dinilai lebih robust untuk kondisi data ini.

4.2.3 Uji Stasioneritas Varians

Uji stasioneritas varians dilakukan untuk mengevaluasi apakah varians penemuan kasus baru kusta bersifat homogen di seluruh wilayah pengamatan. Pengujian dilakukan menggunakan uji Levene dengan pembagian wilayah berdasarkan arah Barat–Timur dan Utara–Selatan.

Hipotesis

H₀ : Varians penemuan kasus baru kusta pada wilayah arah Barat–Timur (atau Utara–Selatan) adalah sama

H₁ : Varians penemuan kasus baru kusta pada wilayah arah Barat–Timur (atau Utara–Selatan) tidak sama

Kriteria Uji
Tolak H₀ jika p-value < α (0,05)

Tabel 6. Hasil Uji Levene Stasioneritas Varians

Arah Pembagian	Statistik Levene	F-value	p-value
Barat–Timur	L₁	0,9207	0,3473
Utara–Selatan	L₂	5,169	0,03265

Pada pembagian arah Barat–Timur, nilai p-value > 0,05 sehingga hipotesis nol gagal ditolak dan varians dapat dianggap homogen. Sebaliknya, pada pembagian arah Utara–Selatan, nilai p-value lebih kecil dari 0,05 sehingga hipotesis nol ditolak.

Hasil ini menunjukkan bahwa varians penemuan kasus baru kusta bersifat heterogen pada arah Utara–Selatan, sehingga data tidak sepenuhnya memenuhi asumsi stasioneritas varians dan mengindikasikan adanya perbedaan tingkat variasi kasus antar wilayah dengan karakteristik geografis yang berbeda. Ketidakstasioneran varians ini berpotensi memengaruhi kinerja metode interpolasi geostatistik yang mengasumsikan varians konstan, seperti Ordinary Kriging, sehingga metode tersebut kurang sesuai tanpa penyesuaian atau transformasi lebih lanjut. Adapun metode interpolasi deterministik seperti Inverse Distance Weighting (IDW) lebih robust karena tidak bergantung pada asumsi homogenitas varians.

4.2.4 Autokorelasi Spasial

Diagnosis autokorelasi spasial dilakukan untuk mengetahui apakah penemuan kasus baru kusta di suatu kabupaten/kota berkorelasi dengan wilayah lain yang berdekatan secara geografis. Keberadaan autokorelasi spasial menunjukkan bahwa informasi spasial pada kasus baru kusta mengandung pola yang sistematis, sehingga nilai pada suatu lokasi dapat dimanfaatkan untuk memperkirakan nilai pada lokasi lain di sekitarnya.

Pengujian autokorelasi spasial dilakukan menggunakan indeks Moran’s I dengan matriks pembobot spasial queen contiguity.

Hipotesis
H₀ : Tidak terdapat autokorelasi spasial pada penemuan kasus baru kusta
H₁ : Terdapat autokorelasi spasial pada penemuan kasus baru kusta

Kriteria Uji
Tolak H₀ jika p-value < α (0,05)

Ringkasan hasil uji Moran’s I disajikan pada Tabel 6.

Tabel 7. Hasil Uji Autokorelasi Spasial Moran’s I

Statistik	Nilai
Moran’s I	0,3605
Expected I	-0,0435
Variance	0,0209
Z-score	2,7942
p-value	0,0026

Berdasarkan Tabel 7., hasil pengujian autokorelasi spasial menunjukkan adanya autokorelasi spasial positif yang signifikan, yang mengindikasikan bahwa wilayah dengan tingkat penemuan kasus kusta yang tinggi cenderung berdekatan dengan wilayah lain dengan tingkat kasus yang juga tinggi, demikian pula sebaliknya. Temuan ini menegaskan bahwa data penemuan kasus baru kusta memiliki struktur spasial yang jelas dan layak untuk dilakukan interpolasi spasial.

4.2.5 Anisotropi

Diagnosis anisotropi dilakukan untuk mengetahui apakah struktur variabilitas spasial berbeda menurut arah tertentu. Analisis ini dilakukan menggunakan variogram terarah pada sudut 0°, 45°, 90°, dan 135°.

Gambar 5. Variogram Terarah

Hasil variogram terarah pada Gambar 5. menunjukkan perbedaan pola semivarian antar arah, baik dari sisi tingkat kenaikan semivarian maupun jarak jangkauan (range). Perbedaan ini mengindikasikan bahwa struktur spasial penemuan kasus baru kusta bersifat anisotropik, sehingga pengaruh jarak antarwilayah tidak sama ke segala arah.

4.2.6 Kerapatan Data

Kerapatan dan sebaran titik observasi merupakan aspek penting dalam menentukan keandalan interpolasi spasial. Diagnosis dilakukan melalui variogram eksperimental, peta densitas titik observasi, serta evaluasi jarak rata-rata antar titik terdekat.

Gambar 6. Variogram Eksperimental

Berdasarkan variogram eksperimental yang ditampilkan pada Gambar 6, terlihat bahwa nilai semivarians tidak meningkat secara stabil seiring bertambahnya jarak, melainkan berfluktuasi naik dan turun pada jarak tertentu. Pola ini mengindikasikan adanya hole effect, yaitu kondisi ketika kemiripan spasial melemah pada jarak menengah namun kembali muncul pada jarak yang lebih jauh. Hole effect mencerminkan struktur spasial yang bersifat periodik atau berulang, sehingga hubungan spasial tidak mengikuti pola sederhana. Kondisi ini menyulitkan pembentukan variogram teoretis yang stabil dan menunjukkan bahwa interpolasi geostatistik standar perlu diterapkan dengan kehati-hatian.

Gambar 7. Peta Densitas Titik Observasi

Sementara itu, peta sebaran dan densitas titik observasi pada Gambar 7. memperlihatkan bahwa distribusi titik kabupaten/kota tidak sepenuhnya merata di seluruh wilayah Jawa Barat. Beberapa wilayah menunjukkan konsentrasi titik observasi yang lebih tinggi, sedangkan wilayah lainnya relatif jarang memiliki titik pengamatan. Ketimpangan sebaran ini berpotensi memengaruhi tingkat ketelitian interpolasi, khususnya pada wilayah yang relatif jauh dari titik observasi.

Untuk memberikan ringkasan kuantitatif mengenai kerapatan titik observasi, dihitung rata-rata jarak tetangga terdekat (mean nearest neighbor distance), yaitu jarak dari setiap titik observasi kabupaten/kota ke titik observasi terdekat lainnya, kemudian dirata-ratakan secara keseluruhan. Ringkasan hasil analisis tersebut disajikan pada Tabel 8.

Tabel 8. Ringkasan Kerapatan Titik Observasi Kabupaten/Kota

Indikator Kerapatan	Nilai
Jumlah titik observasi	25 kabupaten/kota
Rata-rata jarak tetangga terdekat	0,208
Interpretasi umum	Titik observasi relatif rapat

Nilai rata-rata jarak tetangga terdekat sebesar 0,208 menunjukkan bahwa secara umum jarak antar titik observasi kabupaten/kota di Jawa Barat relatif berdekatan. Kondisi ini mengindikasikan bahwa informasi spasial dari wilayah sekitar masih cukup representatif untuk digunakan dalam proses interpolasi berbasis jarak.

4.2.7 Identifikasi Metode

Berdasarkan seluruh rangkaian diagnosis data interpolasi, dapat disimpulkan bahwa penemuan kasus baru kusta memiliki karakteristik spasial yang kompleks. Data tidak berdistribusi normal, analisis tren spasial menunjukkan adanya ketidakstasioneran rata-rata dengan pola nonlinier yang signifikan, didukung oleh perbandingan nilai AIC serta pola semivariogram yang tidak mencapai sill. Selain itu, uji stasioneritas varians mengindikasikan adanya heterogenitas varians pada arah Utara–Selatan, sementara analisis anisotropi menunjukkan bahwa struktur spasial berbeda menurut arah tertentu. Di sisi lain, hasil uji Moran’s I membuktikan adanya autokorelasi spasial positif yang signifikan, serta analisis kerapatan titik observasi menunjukkan bahwa distribusi dan jarak antar titik relatif memadai untuk interpolasi berbasis jarak. Kombinasi kondisi data yang tidak sepenuhnya memenuhi asumsi stasioneritas global, keberadaan anisotropi, serta keterbatasan jumlah titik observasi menjadikan metode Inverse Distance Weighting (IDW) sebagai pendekatan yang paling sesuai, karena IDW tidak mensyaratkan stasioneritas ketat, tidak memerlukan pemodelan variogram teoritis, dan mampu memanfaatkan kedekatan spasial antarwilayah secara langsung dalam mengestimasi nilai pada lokasi yang memiliki data hilang.

4.3 Inverse Distance Weighting (IDW)

Metode Inverse Distance Weighting (IDW) digunakan untuk melakukan interpolasi spasial dan imputasi missing value pada data penemuan kasus baru kusta di Jawa Barat, khususnya pada Kabupaten Ciamis dan Kota Cimahi. Sebelum dilakukan interpolasi akhir, dilakukan optimasi parameter power (p) dan jumlah tetangga (k) menggunakan pendekatan Leave-One-Out Cross Validation (LOOCV) guna memperoleh kombinasi parameter yang menghasilkan galat prediksi minimum.

4.3.1 Optimasi Parameter Power (p)

Optimasi parameter power (p) dilakukan untuk menentukan tingkat pengaruh jarak terhadap bobot interpolasi. Nilai p yang semakin besar akan meningkatkan dominasi titik-titik yang lebih dekat terhadap lokasi prediksi. Ringkasan hasil optimasi parameter p disajikan pada Tabel 9.

Tabel 9. Hasil Optimisasi Parameter Power (p)

Power (p)	RMSE	MAE
1	3,251	2,463
2	2,990	2,093
3	2,851	1,794
4	2,805	1,668
5	2,801	1,612
6	2,816	1,582

Berdasarkan Tabel 9., nilai RMSE mengalami penurunan seiring peningkatan nilai p hingga mencapai minimum pada p = 5. Meskipun nilai MAE sedikit lebih kecil pada p = 6, perbedaan yang diperoleh relatif kecil dan diiringi dengan peningkatan kembali nilai RMSE. Oleh karena itu, parameter p = 5 dipilih sebagai nilai optimal karena memberikan keseimbangan terbaik antara akurasi dan stabilitas prediksi.

4.3.2 Optimasi Parameter Tetangga (k)

Selain parameter p, jumlah tetangga terdekat (k) yang digunakan dalam proses interpolasi juga memengaruhi hasil estimasi. Nilai k yang terlalu kecil dapat menyebabkan hasil prediksi menjadi tidak stabil, sedangkan nilai k yang terlalu besar dapat mengaburkan pengaruh lokal. Optimasi dilakukan dengan menguji beberapa nilai k pada parameter p = 5. Hasil evaluasi disajikan pada Tabel 10.

Tabel 10. Hasil Optimasi Parameter Tetangga (k)

Jumlah Tetangga (k)	RMSE	MAE
3	2,865	1,576
5	2,782	1,604
8	2,796	1,618
10	2,792	1,612
15	2,800	1,613
20	2,801	1,612

Hasil pada Tabel 10. menunjukkan bahwa nilai RMSE terendah diperoleh pada k = 5. Meskipun nilai MAE terendah diperoleh pada k = 3, penggunaan jumlah tetangga yang terlalu sedikit berpotensi meningkatkan ketidakstabilan prediksi. Dengan mempertimbangkan akurasi dan representativitas spasial, nilai k = 5 dipilih sebagai parameter tetangga optimal.

4.3.3 Prediksi dan Imputasi Nilai Hilang

Berdasarkan hasil optimasi, interpolasi IDW dilakukan menggunakan parameter optimal p = 5 dan k = 5. Parameter ini kemudian digunakan untuk mengestimasi nilai penemuan kasus baru kusta pada wilayah yang memiliki missing value. Hasil imputasi disajikan pada Tabel 11.

Tabel 11. Nilai Penemuan Kasus Baru Kusta Hasil Imputasi IDW

Wilayah	Estimasi Kusta
Kabupaten Ciamis	0,733
Kota Cimahi	0,219

Hasil ini menunjukkan bahwa metode IDW mampu memberikan estimasi numerik pada wilayah yang sebelumnya tidak memiliki data, sehingga dataset menjadi lengkap dan dapat digunakan untuk analisis spasial lanjutan.

4.3.4 Peta Residual Interpolasi IDW

Untuk mengevaluasi kinerja spasial hasil interpolasi, dilakukan analisis residual menggunakan nilai hasil prediksi LOOCV. Residual dihitung sebagai selisih antara nilai aktual dan nilai hasil prediksi IDW. Peta sebaran residual ditampilkan pada Gambar 8.

Gambar 8. Peta Residual Interpolasi IDW Penemuan Kasus Baru Kusta

Berdasarkan Gambar 8, residual dengan nilai positif (ditampilkan dengan warna merah) menunjukkan kondisi underestimate, sedangkan residual bernilai negatif (warna ungu) menunjukkan overestimate. Secara umum, sebagian besar wilayah memiliki residual yang relatif kecil, namun terdapat satu titik dengan residual positif paling besar yang berada di Kabupaten Indramayu. Hal ini menunjukkan bahwa hasil interpolasi IDW di wilayah tersebut secara nyata lebih rendah dibandingkan nilai observasi aktual, sehingga terjadi kesalahan prediksi yang cukup besar. Kondisi ini diduga dipengaruhi oleh variasi lokal yang kuat atau keterbatasan titik pengamatan di sekitar Indramayu yang tidak sepenuhnya tertangkap oleh pembobotan jarak pada metode IDW, meskipun secara keseluruhan pola residual tidak membentuk klaster spasial yang kuat.

4.3.5 Visualisasi Permukaan Interpolasi

Hasil akhir interpolasi IDW divisualisasikan dalam bentuk permukaan kontinu menggunakan raster dan garis kontur untuk menggambarkan pola spasial estimasi penemuan kasus baru kusta di Jawa Barat. Visualisasi permukaan interpolasi ditampilkan pada Gambar 9.

Gambar 9. Permukaan Interpolasi Penemuan Kasus Baru Kusta di Jawa Barat Menggunakan Metode IDW

Gambar 9. menunjukkan adanya variasi spasial yang jelas pada estimasi penemuan kasus baru kusta. Area dengan intensitas warna paling terang berada di Kabupaten Indramayu, yang mengindikasikan estimasi penemuan kasus baru kusta tertinggi dibandingkan wilayah sekitarnya. Pola permukaan interpolasi yang halus mencerminkan karakteristik metode IDW yang menekankan pengaruh jarak dan kedekatan spasial antarwilayah, sehingga nilai pada suatu lokasi sangat dipengaruhi oleh titik pengamatan di sekitarnya. Pola ini menunjukkan bahwa konsentrasi kasus yang tinggi bersifat lokal dan menurun secara bertahap seiring bertambahnya jarak dari pusat wilayah dengan estimasi tinggi.

4.2 Ordinary Least Square untuk Model Klasik

Sebagai pendekatan awal dalam menganalisis hubungan antara penemuan kasus baru kusta dan faktor-faktor penjelasnya, digunakan model regresi linier dengan metode Ordinary Least Squares (OLS). Model ini berfungsi sebagai model pembanding sebelum diterapkan pendekatan yang mempertimbangkan ketergantungan spasial antarwilayah.

Model regresi linier yang digunakan dapat dituliskan sebagai berikut :

\[ \begin{equation} KUSTA_i = -4.912 + 4.06 \times 10^{-5} \, KEPADATAN_i + 0.05774 \, PHBS_i + 0.4883 \, MISKIN_i + \varepsilon_i, \end{equation} \]

Selanjutnya, dilakukan pengujian asumsi-asumsi klasik OLS secara keseluruhan untuk mengevaluasi kelayakan model. Ringkasan hasil uji asumsi disajikan pada Tabel 12.

Tabel 12. Hasil Uji Asumsi Model OLS

Asumsi	Uji yang Digunakan	Hasil	Kesimpulan
Multikolinearitas	Variance Inflation Factor (VIF)	VIF semua X < 10	Terpenuhi
Homoskedastisitas	Breusch–Pagan Test	p-value = 0,4861	Terpenuhi
Normalitas residual	Kolmogorov–Smirnov Test	p-value = 0,0986	Terpenuhi
Linearitas	Ramsey RESET Test	p-value = 0,2980	Terpenuhi

Berdasarkan Tabel 12, seluruh asumsi klasik OLS secara statistik terpenuhi. Tidak ditemukan indikasi multikolinearitas serius antarvariabel penjelas, varians residual bersifat homogen, residual berdistribusi normal, serta bentuk fungsional model tidak menunjukkan kesalahan spesifikasi yang signifikan.

Meskipun demikian, pemenuhan asumsi klasik OLS belum cukup untuk menjamin kelayakan model dalam konteks data spasial. Diperoleh juga nilai Multiple R-squared sebesar 0,1047 atau 10,47% dan AIC sebesar 151,02. Hal ini menunjukkan kemampuan model OLS dalam menjelaskan variasi penemuan kasus baru kusta masih relatif rendah. Oleh karena itu, dilakukan pengujian autokorelasi spasial terhadap residual model OLS menggunakan indeks Moran’s I dengan matriks pembobot spasial queen contiguity.

4.3 Uji Moran’s I Residual OLS

Hasil uji Moran’s I terhadap residual model OLS disajikan pada tabel berikut.

Tabel 13. Hasil Uji Moran’s I

Indikator	Nilai
Moran’s I	0,3724
Expected Moran’s I	−0,0385
Varians	0,0180
Z-statistic	3,0603
p-value	0,0011

Diperoleh nilai statistik Moran’s I yang positif dan signifikan secara statistik (p-value < 0,05). Hal ini mengindikasikan adanya autokorelasi spasial positif pada residual, yang berarti bahwa kesalahan prediksi pada suatu kabupaten/kota berkorelasi dengan kesalahan pada wilayah lain yang berdekatan secara geografis.

Keberadaan autokorelasi spasial pada residual melanggar asumsi independensi error dalam model OLS dan menunjukkan bahwa model linier global ini belum mampu menangkap struktur spasial yang melekat pada data penemuan kasus baru kusta. Dengan demikian, meskipun asumsi klasik OLS terpenuhi secara statistik, model OLS dinilai tidak memadai untuk digunakan sebagai model utama. Oleh karena itu, perlu dilakukan pemodelan menggunakan odel yang dapat mengakomodasi spasial pada data.

4.4 Autokorelasi Spasial

Analisis autokorelasi spasial dilakukan untuk mengidentifikasi pola keterkaitan antarwilayah yang bertetangga secara geografis. Pengujian menggunakan matriks bobot spasial queen contiguity dan mencakup autokorelasi global (Moran’s I dan Geary’s C) serta autokorelasi lokal (Local Moran’s I/LISA dan Getis–Ord).

4.4.1 Moran’s I

Hasil uji Moran’s I global disajikan pada Tabel 14. berikut ini.

Tabel 14. Hasil Uji Moran’s I Setiap Variabel

Variabel	Moran’s I	p-value (Randomisasi)	p-value (Monte Carlo)	Kesimpulan
Kasus Baru Kusta	0,402	0,000	0,003	Autokorelasi spasial positif signifikan
Kepadatan Penduduk	0,258	0,016	0,025	Autokorelasi spasial positif signifikan
PHBS	−0,095	0,657	0,625	Tidak terdapat autokorelasi spasial
Persentase Penduduk Miskin	0,478	0,000	0,001	Autokorelasi spasial positif signifikan

Berdasarkan tabel tersebut, Kasus Baru Kusta, Kepadatan Penduduk, dan Persentase Penduduk Miskin memiliki nilai Moran’s I positif dan signifikan secara statistik (p-value < 0,05), baik pada uji randomisasi maupun uji Monte Carlo. Hal ini menunjukkan adanya autokorelasi spasial positif, di mana wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi, demikian pula untuk wilayah bernilai rendah. Sebaliknya, variabel PHBS tidak menunjukkan autokorelasi spasial yang signifikan, sehingga distribusinya cenderung acak secara geografis.

4.4.2 Geary’s C

Sebagai pelengkap Moran’s I, uji Geary’s C digunakan untuk mendeteksi autokorelasi spasial dengan sensitivitas yang lebih tinggi terhadap variasi lokal. Hasil uji Geary’s C yang disajikan pada Tabel 15. berikut ini.

Tabel 15. Hasil Uji Geary’s C

Variabel	Geary’s C	p-value	Kesimpulan
Kasus Baru Kusta	0,521	0,008	Autokorelasi spasial positif signifikan
Kepadatan Penduduk	0,594	0,008	Autokorelasi spasial positif signifikan
PHBS	1,033	0,579	Tidak terdapat autokorelasi spasial
Persentase Penduduk Miskin	0,410	< 0,001	Autokorelasi spasial positif signifikan

Hasil menunjukkan pola yang konsisten dengan hasil uji Moran’s I. Nilai Geary’s C yang secara signifikan lebih kecil dari satu pada variabel Kasus Baru Kusta, Kepadatan Penduduk, dan Persentase Penduduk Miskin mengindikasikan adanya kemiripan nilai antarwilayah yang bertetangga, yang mencerminkan autokorelasi spasial positif. Sebaliknya, variabel PHBS memiliki nilai Geary’s C yang tidak signifikan dan mendekati satu, sejalan dengan hasil Moran’s I, sehingga menunjukkan bahwa distribusi PHBS cenderung acak secara geografis dan tidak memiliki ketergantungan spasial yang berarti.

4.4.3 Local Moran’s I (LISA)

Analisis Local Moran’s I (LISA) digunakan untuk mengidentifikasi pola klaster spasial lokal serta spatial outlier pada masing-masing wilayah. Kategori klaster High–High menunjukkan wilayah dengan nilai variabel tinggi yang dikelilingi oleh wilayah bertetangga dengan nilai tinggi pula, sedangkan Low–Low menunjukkan wilayah bernilai rendah yang dikelilingi oleh wilayah bernilai rendah. Sementara itu, High–Low dan Low–High merepresentasikan spatial outlier, yaitu wilayah yang nilainya kontras dengan wilayah sekitarnya.

Gambar 10. Hasil Klaster Local Moran’s I (LISA)

Berdasarkan peta klaster LISA pada Gambar 10., untuk variabel Kasus Baru Kusta, teridentifikasi klaster High–High di Kabupaten Cirebon. Hal ini menunjukkan bahwa Kabupaten Cirebon memiliki jumlah kasus baru kusta yang relatif tinggi dan dikelilingi oleh kabupaten/kota dengan tingkat kasus yang juga relatif tinggi. Kondisi ini mengindikasikan adanya konsentrasi spasial penemuan kasus kusta dan potensi pengaruh faktor regional, seperti mobilitas penduduk atau kesamaan kondisi sosial-lingkungan antarwilayah bertetangga.

Pada variabel Persentase Penduduk Miskin, klaster High–High ditemukan di Kabupaten Cirebon dan Kabupaten Majalengka, yang menunjukkan bahwa wilayah-wilayah tersebut memiliki tingkat kemiskinan tinggi dan berada dalam lingkungan spasial dengan karakteristik serupa. Sebaliknya, klaster Low–Low teridentifikasi di Kabupaten Bogor dan Kota Bekasi, yang merepresentasikan wilayah dengan tingkat kemiskinan relatif rendah dan dikelilingi oleh wilayah dengan kondisi sosial ekonomi yang juga lebih baik.

Sementara itu, variabel Kepadatan Penduduk dan PHBS tidak menunjukkan klaster LISA yang signifikan, yang mengindikasikan bahwa meskipun terdapat variasi antarwilayah, tidak ditemukan pola klaster lokal yang kuat secara statistik.

4.4.4 Getis–Ord 𝐺∗𝑖

Analisis Getis–Ord G*i​ bertujuan untuk mengidentifikasi wilayah yang berperan sebagai hot spot (konsentrasi nilai tinggi) dan cold spot (konsentrasi nilai rendah). Berbeda dengan LISA yang membedakan tipe klaster dan outlier, metode ini lebih menekankan pada intensitas pengelompokan nilai ekstrem secara lokal.

Gambar 11. Hasil Klaster Getis–Ord 𝐺∗𝑖

Berdasarkan Gambar 11, untuk variabel Kasus Baru Kusta, Kabupaten Cirebon teridentifikasi sebagai hot spot yang signifikan. Hal ini mengindikasikan bahwa wilayah tersebut berada dalam area dengan konsentrasi nilai kasus kusta yang tinggi dibandingkan wilayah sekitarnya, sehingga berperan sebagai pusat konsentrasi spasial penemuan kasus.

Pada variabel Persentase Penduduk Miskin, Kabupaten Cirebon dan Kabupaten Majalengka teridentifikasi sebagai hot spot, sedangkan Kabupaten Bogor dan Kota Bekasi sebagai cold spot. Temuan ini konsisten dengan hasil LISA dan menguatkan adanya ketimpangan spasial sosial ekonomi antarwilayah di Provinsi Jawa Barat. Untuk variabel Kepadatan Penduduk dan PHBS, tidak ditemukan hot spot maupun cold spot yang signifikan, yang menunjukkan tidak adanya konsentrasi spasial nilai ekstrem.

4.5 Uji Lagrange Multiplier

Setelah model Ordinary Least Squares (OLS) diestimasi dan dilakukan pengujian asumsi klasik, tahap selanjutnya adalah menguji apakah terdapat ketergantungan spasial yang tidak tertangkap oleh model OLS. Untuk tujuan tersebut digunakan Uji Lagrange Multiplier (LM), yang bertujuan mendeteksi keberadaan dependensi spasial baik pada komponen error maupun pada variabel dependen. Ringkasan hasil uji disajikan pada Tabel 16.

Tabel 16. Hasil Uji Lagrange Multiplier

Jenis Uji	Statistik Uji	Derajat Bebas	p-value	Kesimpulan
LM Error (RSerr)	6,0408	1	0,01398	Signifikan
LM Lag (RSlag)	6,4478	1	0,01111	Signifikan
Robust LM Error (adjRSerr)	0,0135	1	0,9076	Tidak signifikan
Robust LM Lag (adjRSlag)	0,4205	1	0,5167	Tidak signifikan
SARMA	6,4613	2	0,03953	Signifikan

Hasil uji LM menunjukkan bahwa LM Error dan LM Lag sama-sama signifikan pada taraf signifikansi 5%. Temuan ini mengindikasikan bahwa model OLS mengandung ketergantungan spasial, baik yang bersumber dari komponen error maupun dari variabel dependen. Dengan kata lain, model linier klasik belum mampu menangkap sepenuhnya struktur spasial yang terdapat dalam data kasus baru kusta.

Namun, ketika dilakukan pengujian versi robust, baik Robust LM Error maupun Robust LM Lag tidak menunjukkan hasil yang signifikan. Kondisi ini mengindikasikan bahwa sulit untuk menentukan secara tegas apakah bentuk ketergantungan spasial yang dominan berasal dari mekanisme lag atau error, karena keduanya saling memengaruhi. Uji SARMA yang signifikan memperkuat temuan tersebut, yang menunjukkan adanya kemungkinan kombinasi ketergantungan spasial pada variabel dependen dan error secara simultan.

Dengan mempertimbangkan hasil Uji Lagrange Multiplier dan uji SARMA, maka model kandidat yang layak dipertimbangkan dalam analisis lanjutan meliputi Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Autoregressive Combined Model (SAC).

4.6 Perbandingan Model

4.6.1 Perbandingan Model Berdasarkan AIC

Perbandingan awal antar model dilakukan menggunakan Akaike Information Criterion (AIC) dan (Bayesian Information Criterion) BIC untuk menentukan model dengan keseimbangan terbaik antara kecocokan dan kompleksitas.

Tabel 17. Hasil Pemodelan Model Spasial

Model	AIC	BIC
SAR	144.64	152.42
SEM	144.69	152.46
SAC	145.86	154.93
SDM	148.77	160.44

Berdasarkan Tabel 17, model SAR memiliki nilai AIC dan BIC terendah sehingga model tersebut dipilih sebagai model utama.

4.6.2 Uji Likelihood Ratio Model Terpilih

Setelah model terbaik ditentukan berdasarkan kriteria AIC, dilakukan uji Likelihood Ratio (LR) untuk menguji apakah model SAR memberikan peningkatan kecocokan yang signifikan dibandingkan model OLS.

Hipotesis​​

H0 : Model OLS cukup, tidak diperlukan efek spasial

H1 : Model SAR lebih baik dibandingkan model OLS​

Kriteria Uji

Tolak H0 jika p-value ≤ α

Tabel 18. Hasil Uji Likelihood Ratio

Statistik	Nilai
Rho (ρ)	0.6431
Likelihood Ratio (LR)	8.3787
p-value	0.0038

Hasil uji Likelihood Ratio menunjukkan nilai p-value lebih kecil dari taraf signifikansi 5%, sehingga hipotesis nol ditolak. Hal ini mengindikasikan bahwa model SAR memberikan peningkatan kecocokan yang signifikan dibandingkan model OLS. Dengan demikian, keberadaan ketergantungan spasial pada variabel Kasus Baru Kusta perlu dimodelkan secara eksplisit menggunakan pendekatan regresi spasial.

4.7 Estimasi Model SAR

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

\[ y_i = 0.6431^{***} \sum_{j=1}^{n} w_{ij} y_j - 2.9603 - 3.2007 \times 10^{-6} (KEPADATAN)_i + 0.0223 (PHBS)_i + 0.3210 (MISKIN)_i \]

Model ini memasukkan komponen lag spasial dari variabel dependen yang direpresentasikan oleh parameter ρ. Hasil estimasi menunjukkan bahwa parameter autokorelasi spasial ρ bernilai positif dan signifikan secara statistik, yang mengindikasikan adanya pengaruh spasial yang kuat. Artinya, jumlah kasus baru kusta di suatu wilayah dipengaruhi oleh jumlah kasus di wilayah-wilayah yang bertetangga. Sementara itu, koefisien dari seluruh variabel penjelas tidak signifikan secara statistik, yang menunjukkan bahwa variasi kasus baru kusta lebih dominan dijelaskan oleh mekanisme ketergantungan spasial antarwilayah dibandingkan oleh pengaruh langsung variabel-variabel tersebut.

Adapun hasil perhitungan impact pada model Spatial Autoregressive (SAR) disajikan pada Tabel 19. berikut.

Tabel 19. Hasil Perhitungan Impact Model SAR

Variabel	Jenis Impact	Nilai Impact	p-value	Keputusan
Kepadatan Penduduk	Direct	-3.71×10⁻⁶	0.940	Tidak Signifikan
	Indirect	-5.26×10⁻⁶	0.957	Tidak Signifikan
	Total	-8.97×10⁻⁶	0.956	Tidak Signifikan
PHBS	Direct	0.0258	0.924	Tidak Signifikan
	Indirect	0.0367	0.961	Tidak Signifikan
	Total	0.0625	0.960	Tidak Signifikan
Penduduk Miskin	Direct	0.3718	0.671	Tidak Signifikan
	Indirect	0.5275	0.944	Tidak Signifikan
	Total	0.8993	0.932	Tidak Signifikan

Berdasarkan hasil perhitungan impact pada model Spatial Autoregressive (SAR), tidak terdapat variabel independen yang menunjukkan pengaruh signifikan, baik secara langsung (direct impact) maupun tidak langsung (indirect impact), terhadap jumlah kasus kusta pada tingkat signifikansi 5%.

Variabel kepadatan penduduk (jiwa/km²) memiliki nilai direct impact negatif, yang mengindikasikan bahwa peningkatan kepadatan penduduk cenderung diikuti oleh penurunan jumlah kasus kusta di wilayah yang sama. Namun, pengaruh tersebut tidak signifikan secara statistik. Efek tidak langsung kepadatan penduduk juga tidak signifikan, yang menunjukkan bahwa perubahan kepadatan penduduk di wilayah tetangga tidak memberikan limpahan pengaruh yang berarti terhadap kasus kusta di suatu wilayah.

Variabel PHBS (persentase rumah tangga yang menerapkan Perilaku Hidup Bersih dan Sehat) menunjukkan nilai direct impact dan indirect impact yang positif, meskipun tidak signifikan. Hal ini mengindikasikan bahwa peningkatan PHBS belum dapat menunjukkan pengaruh yang terukur secara statistik terhadap penurunan kasus kusta, baik secara langsung maupun melalui mekanisme spasial antarwilayah.

Sementara itu, variabel persentase penduduk miskin memiliki nilai direct, indirect, dan total impact yang relatif lebih besar dibandingkan variabel lain, namun seluruh pengaruh tersebut juga tidak signifikan secara statistik. Temuan ini mengindikasikan bahwa meskipun kemiskinan berpotensi berkaitan dengan penyebaran kusta, pengaruhnya dalam model ini belum cukup kuat untuk terdeteksi secara signifikan, baik pada tingkat wilayah itu sendiri maupun melalui efek spillover spasial.

Secara keseluruhan, hasil impact analysis menunjukkan bahwa ketergantungan spasial dalam model SAR lebih banyak ditangkap melalui parameter lag spasial (\(\rho\)), dibandingkan melalui efek langsung maupun tidak langsung dari variabel-variabel penjelas. Hal ini menegaskan bahwa dinamika spasial kasus kusta lebih dipengaruhi oleh interaksi antarwilayah daripada oleh variasi karakteristik sosial ekonomi yang dimodelkan secara eksplisit.

4.7.1 Uji Asumsi Model

Untuk memastikan kelayakan model SAR, dilakukan pengujian asumsi terhadap residual model. Ringkasan hasil uji asumsi disajikan pada Tabel 20.

Tabel 20. Hasil Uji Asumsi Model SAR

Asumsi	Uji	p-value	Kesimpulan
Tidak ada autokorelasi spasial residual	Moran’s I	0,190	Terpenuhi
Normalitas residual	Kolmogorov–Smirnov	0,530	Terpenuhi
Homoskedastisitas	Breusch–Pagan	0,383	Terpenuhi
Linearitas	RESET	0,458	Terpenuhi

Seluruh uji asumsi menunjukkan bahwa model SAR memenuhi kriteria diagnostik yang diperlukan dan tidak menunjukkan pelanggaran asumsi yang signifikan. Hal ini berlaku juga untuk independensi residual, dimana sudah tidak terdapat autokorelasi spasial yang signifikan pada residual, sehingga ketergantungan spasial yang sebelumnya terdeteksi telah berhasil diakomodasi oleh komponen lag spasial dalam model SAR.

4.7.2 Uji Sensitivitas Bobot

Uji sensitivitas bobot dilakukan untuk menilai kestabilan hasil model SAR terhadap perubahan definisi matriks bobot spasial. Bobot yang digunakan meliputi Queen, Rook, dan KNN-4.

Tabel 21. Hasil Uji Sensitivitas Bobot

Jenis Bobot	Rho	AIC	Log-Likelihood
Queen	0,643	144,64	-66,32
Rook	0,643	144,64	-66,32
KNN-4	0,589	145,49	-66,74

Hasil uji sensitivitas bobot menunjukkan bahwa nilai parameter ketergantungan spasial (ρ) dan kriteria kecocokan model (AIC dan log-likelihood) relatif konsisten ketika menggunakan matriks bobot Queen dan Rook, yang menandakan bahwa struktur ketetanggaan berbasis batas wilayah tidak secara signifikan memengaruhi hasil estimasi model. Pada bobot KNN-4, meskipun nilai ρ sedikit lebih rendah dan AIC meningkat, perubahannya tidak bersifat drastis. Temuan ini mengindikasikan bahwa estimasi parameter spasial model SAR cukup stabil terhadap variasi definisi bobot spasial, sehingga hasil model dapat dianggap robust dan tidak sensitif terhadap pemilihan matriks bobot tertentu.

4.7.3 Hasil Prediksi Model

Secara keseluruhan, hasil analisis menunjukkan bahwa kasus baru kusta di Jawa Barat memiliki ketergantungan spasial yang kuat, di mana kondisi di suatu kabupaten/kota dipengaruhi oleh kondisi wilayah sekitarnya. Model SAR terbukti lebih mampu menangkap pola tersebut dibandingkan model nonspasial.

Nilai pseudo-R squared sebesar 0,444 menunjukkan bahwa hampir 44% variasi kasus baru kusta dapat dijelaskan oleh model SAR, yang mencerminkan peningkatan R Squared sebesar 33,97% dibandingkan OLS.

Gambar 12. Peta Prediksi dan Residual Model SAR 

Selanjutnya, peta prediksi dan residual pada Gambar 12. memperlihatkan bahwa model SAR mampu merepresentasikan pola spasial kasus baru kusta dengan lebih baik. Peta residual model SAR menunjukkan bahwa residual cenderung lebih besar pada wilayah bagian utara dibandingkan wilayah selatan, dengan Kabupaten Indramayu sebagai lokasi dengan residual positif paling tinggi. Pola ini mengindikasikan bahwa meskipun ketergantungan spasial global telah dimodelkan melalui SAR, masih terdapat variasi sistematis antarwilayah yang belum sepenuhnya terjelaskan. Perbedaan sebaran residual antara wilayah utara dan selatan menunjukkan adanya heterogenitas spasial, di mana hubungan antara variabel dependen dan variabel penjelas tidak bersifat konstan secara geografis. Oleh karena itu, analisis selanjutnya dilakukan menggunakan Geographically Weighted Regression (GWR) untuk menangkap variasi hubungan lokal tersebut.

4.11 Estimasi Model Geographically Weighted Regression (GWR)

Model Spatial Autoregressive (SAR) sebelumnya merupakan model global yang mengasumsikan pengaruh variabel penjelas bersifat seragam di seluruh wilayah. Namun, karena adanya potensi heterogenitas spasial, maka digunakan model Geographically Weighted Regression (GWR).

Estimasi model GWR dilakukan menggunakan adaptive bandwidth 0,1873 yang berarti bahwa pada setiap lokasi, estimasi koefisien GWR dilakukan dengan menggunakan sekitar 18,7% dari seluruh titik pengamatan terdekat. Karena jumlah total unit analisis adalah 27 wilayah, proporsi tersebut setara dengan ±5 tetangga terdekat di sekitar setiap kabupaten/kota. Dengan pendekatan ini, setiap regresi lokal dibangun menggunakan jumlah tetangga yang sama, tetapi dengan radius jarak yang dapat berbeda antarwilayah, menyesuaikan kerapatan spasial data. Adapun ringkasan koefisien lokal yang dihasilkan disajikan pada Tabel 22.

Tabel 22. Ringkasan Estimasi Koefisien Lokal GWR

Variabel	Min	Median	Max	Global
Intercept	-21.435	-7.409	10.397	-4.9118
Kepadatan Penduduk	-0.000481	-0.000044	0.000339	0.0000
PHBS	0.01357	0.07666	0.15215	0.0577
Penduduk Miskin	-0.74677	0.44109	1.5519	0.4883

Perbedaan nilai minimum, median, dan maksimum pada koefisien lokal menunjukkan adanya heterogenitas spasial. Selisih antara koefisien lokal dan koefisien global mengindikasikan bahwa model global kurang mampu menangkap variasi hubungan lokal. Pengaruh setiap variabel berbeda-beda di tiap kabupaten/kota, yang secara visual ditunjukkan pada peta sebaran koefisien pada Gambar 13.

Gambar 13. Peta Sebaran Koefisien Lokal GWR

Peta koefisien GWR menunjukkan bahwa pengaruh Penduduk Miskin, Kepadatan Penduduk, dan PHBS terhadap penemuan kasus baru kusta bervariasi secara spasial. Koefisien cenderung lebih besar dan positif pada wilayah bagian timur dan utara, sementara wilayah barat dan selatan menunjukkan pengaruh yang lebih lemah atau mendekati nol. Pola ini menegaskan bahwa hubungan antarvariabel tidak bersifat global dan seragam, melainkan dipengaruhi oleh konteks lokal masing-masing wilayah, sehingga mendukung keberadaan heterogenitas spasial dalam proses penularan dan penemuan kasus kusta.

Gambar 14. Signifikansi Koefisien Lokal Terhadap Kusta

Berdasarkan hasil uji-t lokal yang disajikan pada Gambar 14., terlihat bahwa variabel Kepadatan Penduduk dan PHBS sebagian besar tidak signifikan di tingkat lokal (Kepadatan hanya signifikan di 1 wilayah yaitu Kota Bekasi, sedangkan PHBS tidak signifikan di wilayah manapun). Sebaliknya, variabel Penduduk Miskin menjadi variabel yang paling dominan dengan pengaruh signifikan secara statistik di 12 kabupaten/kota, diantaranya Subang, Indramayu, Kabupaten Cirebon, Kota Cirebon, Majalengka, Sumedang, Kuningan, Ciamis, Kota Tasikmalaya, Kabupaten Tasikmalaya, Pangandaran, dan Kota Banjar, yang seluruhnya mayoritas berada di bagian timur Jawa Barat. Hal ini menunjukkan bahwa kebijakan penanganan kusta di Jawa Barat perlu memperhatikan karakteristik spesifik ekonomi lokal di setiap wilayah tersebut.

4.12 Uji Asumsi Model GWR

Untuk menjamin reliabilitas estimasi lokal, model Geographically Weighted Regression (GWR) dievaluasi melalui serangkaian uji diagnostik residual.

Tabel 23. Ringkasan Diagnostik Residual Model GWR

Asumsi	Metode Uji	Hasil Statistik	Interpretasi
Normalitas	Kolmogorov–Smirnov	\(p\)-value = 0,4504	Residual berdistribusi normal
Homoskedastisitas	Breusch–Pagan	\(p\)-value = 0,4861	Varians residual bersifat homogen
Independensi Spasial	Moran’s I	\(p\)-value = 0,2506	Tidak ada autokorelasi pada residual
Multikolinearitas	VIF	Maks VIF = 1,99	Tidak ada hubungan antar-variabel independen

Seluruh rangkaian uji diagnostik menunjukkan bahwa model Geographically Weighted Regression (GWR) telah memenuhi semua asumsi statistik yang diperlukan, sehingga parameter lokal yang dihasilkan bersifat valid, tidak bias, dan dapat diandalkan untuk pengambilan keputusan. Keandalan model ini semakin dipertegas dengan terpenuhinya asumsi independensi spasial. Hasil uji Moran’s I menunjukkan bahwa residual bersifat acak secara geografis (p-value = 0,2506), yang berarti tidak ada lagi informasi spasial yang tertinggal atau gagal dijelaskan oleh model. Dengan kata lain, model GWR telah secara sempurna mengakomodasi pola keterkaitan antarwilayah yang sebelumnya menjadi kendala pada model regresi klasik.

4.13 Interpretasi Hasil

Untuk mengevaluasi kinerja model regresi yang digunakan, dilakukan perbandingan ukuran kebaikan model antara OLS, SAR, dan GWR. Hasilnya disajikan pada tabel berikut.

Tabel 24. Ukuran Kebaikan Model OLS, SAR dan GWR

Model	AIC	R-Squared
OLS	151.021	10.47%
SAR	144.642	44.44%
GWR	122.797	69.42%

Berdasarkan Tabel 23., GWR ditetapkan sebagai model akhir terbaik untuk memodelkan penemuan kasus baru kusta di Jawa Barat karena memiliki tingkat akurasi dan efisiensi informasi yang paling unggul dengan nilai AIC terendah dan kemampuan model dalam menjelaskan variasi data tertinggi dengan nilai Quasi-global R-squared mencapai 69,42%. Keunggulan ini juga didukung oleh nilai galat prediksi yang relatif rendah, dengan MAE sebesar 1,318 dan RMSE sebesar 1,926. Analisis terhadap performa lokal model GWR dapat dilihat lebih jelas melalui peta sebaran Local R-Square dan peta distribusi residual berikut ini.

Gambar 15. Peta Local R-Square dan Peta Distribusi Residual

Peta Local R-Square GWR menunjukkan bahwa kemampuan model dalam menjelaskan variasi penemuan kasus baru kusta berbeda antarwilayah, dengan nilai R-square lokal yang umumnya lebih tinggi pada sebagian besar wilayah, terutama di bagian tengah dan timur Jawa Barat. Hal ini mengindikasikan bahwa model GWR mampu menangkap hubungan lokal antarvariabel dengan baik pada wilayah-wilayah tersebut. Sementara itu, peta distribusi residual GWR memperlihatkan bahwa sebagian besar wilayah memiliki residual yang relatif kecil dan tersebar lebih acak dibandingkan model global. Namun demikian, residual dengan intensitas warna paling gelap terlihat di Kabupaten Indramayu, yang menunjukkan masih adanya kesalahan prediksi yang relatif besar di wilayah tersebut. Kondisi ini mengindikasikan bahwa meskipun GWR secara umum meningkatkan performa model, terdapat karakteristik lokal di Indramayu yang belum sepenuhnya tertangkap oleh variabel penjelas yang digunakan.

BAB V KESIMPULAN

5.1 Kesimpulan

Hasil eksplorasi spasial menunjukkan bahwa penemuan kasus baru kusta di Jawa Barat memiliki distribusi yang tidak merata antar kabupaten/kota, dengan terbentuknya konsentrasi kasus pada wilayah tertentu. Pola ini menegaskan bahwa risiko kusta tidak tersebar secara acak secara geografis. Pengujian autokorelasi spasial global menggunakan Moran’s I dan Geary’s C membuktikan adanya autokorelasi spasial positif yang signifikan, yang diperkuat oleh analisis lokal (LISA dan Getis–Ord G*i) melalui identifikasi klaster hot spot, terutama di wilayah Cirebon dan sekitarnya.

Pemodelan menunjukkan bahwa model OLS tidak mampu menangkap struktur spasial data secara memadai, meskipun memenuhi asumsi klasik, karena masih menyisakan autokorelasi spasial pada residual. Model Spatial Autoregressive (SAR) memberikan peningkatan kinerja yang signifikan dan menunjukkan parameter autokorelasi spasial (ρ) yang positif dan signifikan, menegaskan adanya pengaruh antarwilayah. Namun, seluruh variabel penjelas tidak signifikan secara global, mengindikasikan ketidakhomogenan hubungan spasial.

Model Geographically Weighted Regression (GWR) memberikan performa terbaik dengan nilai AIC terendah dan Quasi-global R² sebesar 69,42%, serta mampu mengungkap heterogenitas spasial yang kuat. Faktor kemiskinan muncul sebagai determinan paling dominan secara lokal, terutama di wilayah Timur Jawa Barat. Secara keseluruhan, hasil penelitian menegaskan bahwa penemuan kasus baru kusta dipengaruhi oleh ketergantungan spasial dan variasi lokal, sehingga pendekatan spasial, khususnya GWR, menjadi paling tepat untuk analisis dan perumusan kebijakan berbasis wilayah.

5.2 Saran

5.2.1 Implikasi Kebijakan

Hasil penelitian menunjukkan bahwa pengendalian kusta di Jawa Barat perlu diarahkan secara berbasis wilayah (spatially targeted intervention). Keberadaan klaster spasial dan hot spot kasus kusta mengindikasikan bahwa intervensi yang seragam berpotensi kurang efektif. Oleh karena itu, wilayah dengan konsentrasi kasus tinggi, khususnya di kawasan Timur Jawa Barat, perlu diprioritaskan dalam kegiatan surveilans aktif, skrining dini, dan penguatan layanan kesehatan. Selain itu, temuan bahwa faktor kemiskinan berpengaruh signifikan secara lokal menegaskan pentingnya integrasi kebijakan kesehatan dengan kebijakan sosial ekonomi, seperti pengentasan kemiskinan, peningkatan akses layanan dasar, serta perbaikan kondisi hunian dan lingkungan. Ketergantungan spasial antarwilayah juga menunjukkan perlunya koordinasi lintas kabupaten/kota, sehingga penanganan kusta tidak dilakukan secara terfragmentasi berdasarkan batas administrasi semata.

5.2. Saran Pengembangan Analisis

Penelitian selanjutnya dapat memperkaya analisis dengan menambahkan variabel penjelas yang lebih spesifik, seperti akses layanan kesehatan, mobilitas penduduk, dan faktor lingkungan. Pengembangan analisis spasial-temporal berbasis data panel juga penting untuk mengkaji dinamika kusta dan efektivitas intervensi dari waktu ke waktu. Selain itu, penerapan metode lanjutan seperti Spatial Durbin Model berbasis panel atau Multiscale GWR dapat digunakan untuk menangkap variasi skala pengaruh antar faktor penentu. Pada tahap interpolasi, perbandingan IDW dengan metode geostatistik seperti kriging juga dapat dipertimbangkan apabila asumsi stasioneritas data terpenuhi, guna meningkatkan ketelitian estimasi.

ACKNOWLEDGEMENT

Penulis menyampaikan terima kasih kepada dosen pengampu mata kuliah Spatial Statistics, Dr. I Gede Nyoman Mindra Jaya, M.Si., atas arahan dan materi perkuliahan yang menjadi landasan utama dalam pelaksanaan analisis dan penyusunan laporan ini. Dalam proses pengerjaan tugas akhir semester ini, penulis menggunakan bantuan artificial intelligence (AI) sebagai alat pendukung, khususnya untuk membantu penyusunan sintaks pemrograman analisis spasial, pengembangan dashboard interaktif, serta perapihan dan kejelasan penulisan laporan. Seluruh hasil analisis, interpretasi, dan kesimpulan yang disajikan merupakan hasil pemahaman dan tanggung jawab penuh penulis, serta telah diverifikasi melalui perhitungan dan analisis mandiri sesuai dengan prinsip etika akademik yang berlaku.

DAFTAR PUSTAKA

Anselin, L. (1988). Spatial econometrics: Methods and models. Springer.

Barros, I. D. C. A., Sousa, C. D. C. M., Silva, N. D., & Mascarenhas, M. (2024). Characterization of cases and epidemiological and operational indicators of leprosy: Analysis of time series and spatial distribution, Piauí state, Brazil, 2007–2021. Epidemiologia e Serviços de Saúde, 33, e2023090. https://doi.org/10.1590/s2237-96222024v33e2023090.en

Brunsdon, C., Fotheringham, A. S., & Charlton, M. E. (1996). Geographically weighted regression: A method for exploring spatial nonstationarity. Geographical Analysis, 28(4), 281–298. https://doi.org/10.1111/j.1538-4632.1996.tb00936.x

Burrough, P. A., & McDonnell, R. A. (1998). Principles of geographical information systems. Oxford University Press.

Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically weighted regression: The analysis of spatially varying relationships. Wiley.

Indonesia percepat eliminasi kusta dan filariasis, target bebas NTDs pada 2030. (n.d.). Kementerian Kesehatan Republik Indonesia. Retrieved December 21, 2025, from https://kemkes.go.id/id/indonesia-percepat-eliminasi-kusta-dan-filariasis-target-bebas-ntds-pada-2030

LeSage, J., & Pace, R. K. (2009). Introduction to spatial econometrics (1st ed.). Chapman & Hall/CRC. https://doi.org/10.1201/9781420064254

Lloyd, C. D. (2010). Local models for spatial analysis (2nd ed.). CRC Press.

Lu, B., Charlton, M., Harris, P., & Fotheringham, A. S. (2014). Geographically weighted regression with a non-Euclidean distance metric: A case study using hedonic house price data. International Journal of Geographical Information Science, 28(4), 660–681. https://doi.org/10.1080/13658816.2013.865739

Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 23rd ACM National Conference, 517–524. https://doi.org/10.1145/800186.810616

Sintia, I., Suyitno, S., & Hayati, M. (2022). Geographically weighted Poisson regression model with adaptive bisquare weighting function (Case study: Data on number of leprosy cases in Indonesia 2020). Jurnal Matematika, Statistika dan Komputasi, 19(1). https://doi.org/10.20956/j.v19i1.21879

LAMPIRAN

Script Analisis

Berikut script analisis R yang digunakan dalam analisis ini.

rm(list = ls())
library(sf); library(dplyr); library(stringr); library(spgwr)
library(spdep); library(spatialreg); library(car); library(gstat)
library(tmap); library(broom); library(ggplot2); library(sp)
library(psych); library(readxl); library(viridis); library(spatstat)   
library(scales); library(tibble); library(purrr); library(MASS)
library(forcats); library(htmltools); library(lmtest); library(viridis)
library(gridExtra); library(ggthemes); library(patchwork);library(nortest); library(GWmodel)

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

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

kab_jabar <- c("Bandung","Bandung Barat","Bekasi","Bogor","Ciamis","Cianjur","Cirebon",
               "Garut","Indramayu","Karawang","Kuningan","Majalengka","Pangandaran",
               "Purwakarta","Subang","Sukabumi","Sumedang","Tasikmalaya",
               "Kota Bandung","Kota Bekasi","Kota Bogor","Kota Cimahi",
               "Kota Cirebon","Kota Depok","Kota Sukabumi","Kota Tasikmalaya",
               "Kota Banjar")
kab_jabar

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

indo_adm2 <- indo_adm2 %>% mutate(shapeName_clean = str_to_lower(str_squish(shapeName)))
kab_jabar_clean <- str_to_lower(kab_jabar)
jabar <- indo_adm2 %>% filter(shapeName_clean %in% kab_jabar_clean)
if(nrow(jabar) == 0){
  message("Tidak ada match exact; mencoba fuzzy join...")
  print(unique(indo_adm2$shapeName)[1:80])
  stop("Periksa & adaptasi nama kab/kota lalu jalankan ulang.")
}
attr <- read_xlsx("D:/OneDrive - Universitas Padjadjaran/SMT 5/SPASIAL/PROJECT UAS/DATA SPASIAL.xlsx")
names(attr)

## [1] "KAB/KOTA"  "KUSTA"     "KEPADATAN" "PHBS"      "MISKIN"

attr <- attr %>%
  rename(kabupaten = starts_with("KAB")) %>%
  mutate(kab_clean = str_to_lower(str_squish(kabupaten)))
head(attr)

## # A tibble: 6 × 6
##   kabupaten   KUSTA KEPADATAN  PHBS MISKIN kab_clean  
##   <chr>       <dbl>     <dbl> <dbl>  <dbl> <chr>      
## 1 Bogor        4.07      1899  52.7   7.05 bogor      
## 2 Sukabumi     1.22       679  68     6.87 sukabumi   
## 3 Cianjur      0.73       712  76.3  10.1  cianjur    
## 4 Bandung      0.16      2156  64.2   6.19 bandung    
## 5 Garut        0.61       876  60.2   9.68 garut      
## 6 Tasikmalaya  0.93       710  48.2  10.2  tasikmalaya

jabar <- jabar %>% left_join(attr, by = c("shapeName_clean" = "kab_clean"))
jabar

## Simple feature collection with 27 features and 11 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 KUSTA KEPADATAN  PHBS MISKIN
## 1          bandung       Bandung  0.16      2156 64.24   6.19
## 2    bandung barat Bandung Barat  0.41      1468 64.93  10.49
## 3           bekasi        Bekasi  7.41      2617 80.48   4.80
## 4            bogor         Bogor  4.07      1899 52.66   7.05
## 5           ciamis        Ciamis    NA       789 73.68   7.39
## 6          cianjur       Cianjur  0.73       712 76.27  10.14
## 7          cirebon       Cirebon  7.56      2228 75.11  11.00
## 8            garut         Garut  0.61       876 60.23   9.68
## 9        indramayu     Indramayu 14.98       922 61.63  11.93
## 10        karawang      Karawang  7.03      1335 64.48   7.86
##                          geometry
## 1  MULTIPOLYGON (((107.75 -6.8...
## 2  MULTIPOLYGON (((107.525 -6....
## 3  MULTIPOLYGON (((107.0982 -5...
## 4  MULTIPOLYGON (((106.994 -6....
## 5  MULTIPOLYGON (((108.3857 -7...
## 6  MULTIPOLYGON (((107.3118 -7...
## 7  MULTIPOLYGON (((108.685 -6....
## 8  MULTIPOLYGON (((107.7202 -7...
## 9  MULTIPOLYGON (((108.539 -6....
## 10 MULTIPOLYGON (((107.1717 -6...

# STAT DESK ----
data_spasial <- jabar %>%
  sf::st_drop_geometry() %>%                         
  dplyr::select(KUSTA, KEPADATAN, PHBS, MISKIN)
psych::describe(data_spasial)

##           vars  n    mean      sd  median trimmed     mad    min      max
## KUSTA        1 25    3.29    3.61    1.22    2.82    1.57   0.16    14.98
## KEPADATAN    2 27 3910.93 4668.12 1468.00 3271.52 1120.85 385.00 15176.00
## PHBS         3 27   67.87   10.03   68.00   68.37    9.44  45.48    84.67
## MISKIN       4 27    8.01    2.65    8.41    8.10    2.79   2.34    11.93
##              range  skew kurtosis     se
## KUSTA        14.82  1.43     1.86   0.72
## KEPADATAN 14791.00  1.22    -0.07 898.38
## PHBS         39.19 -0.51    -0.52   1.93
## MISKIN        9.59 -0.36    -0.97   0.51

apply(data_spasial, 2, var, na.rm = TRUE)

##        KUSTA    KEPADATAN         PHBS       MISKIN 
## 1.303110e+01 2.179131e+07 1.006661e+02 7.014903e+00

# Histogram 
ggplot(data_spasial, aes(x = KUSTA)) +
  geom_histogram(bins = 10, fill = "skyblue", color = "black") +
  labs(title = "Distribusi Penemuan Kasus Baru Kusta", x = "Kasus per 100.000 Penduduk", y = "Frekuensi")

ggplot(data_spasial, aes(x = KEPADATAN)) +
  geom_histogram(bins = 10, fill = "lightgreen", color = "black") +
  labs(title = "Distribusi Kepadatan Penduduk", x = "Kepadatan Penduduk (jiwa/km²)", y = "Frekuensi")

ggplot(data_spasial, aes(x = PHBS)) +
  geom_histogram(bins = 10, fill = "lightcoral", color = "black") +
  labs(title = "Distribusi Rumah Tangga Ber-PHBS", x = "Persentase (%)", y = "Frekuensi")

ggplot(data_spasial, aes(x = MISKIN)) +
  geom_histogram(bins = 10, fill = "gold", color = "black") +
  labs(title = "Distribusi Penduduk Miskin", x = "Persentase (%)", y = "Frekuensi")

# Boxplot
ggplot(data_spasial, aes(y = KUSTA)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  labs(title = "Boxplot Kasus Baru Kusta", y = "Kasus per 100.000 Penduduk", x = "") +
  theme_minimal()

ggplot(data_spasial, aes(y = KEPADATAN)) +
  geom_boxplot(fill = "lightgreen", color = "black") +
  labs(title = "Boxplot Kepadatan Penduduk", y = "Kepadatan Penduduk (jiwa/km²)", x = "") +
  theme_minimal()

ggplot(data_spasial, aes(y = PHBS)) +
  geom_boxplot(fill = "lightcoral", color = "black") +
  labs(title = "Boxplot Rumah Tangga Ber-PHBS", y = "Persentase (%)", x = "") +
  theme_minimal()

ggplot(data_spasial, aes(y = MISKIN)) +
  geom_boxplot(fill = "gold", color = "black") +
  labs(title = "Boxplot Penduduk Miskin", y = "Persentase (%)", x = "") +
  theme_minimal()

# Peta Sebaran
plot_kusta <- ggplot(jabar) +
  geom_sf(aes(fill = KUSTA), color = "darkgray", size = 0.2) +
  scale_fill_viridis_c(option = "turbo", na.value = "grey90") +
  labs(title = "Penemuan Kasus Baru Kusta",
       fill = "Kasus per\n100.000 Penduduk") +
  theme_minimal(); plot_kusta

plot_kepadatan <- ggplot(jabar) +
  geom_sf(aes(fill = KEPADATAN), color = "darkgray", size = 0.2) +
  scale_fill_viridis_c(option = "turbo", na.value = "grey90") +
  labs(title = "Kepadatan Penduduk",
       fill = "Jiwa/km²") +
  theme_minimal(); plot_kepadatan

plot_phbs <- ggplot(jabar) +
  geom_sf(aes(fill = PHBS), color = "darkgray", size = 0.2) +
  scale_fill_viridis_c(option = "turbo", na.value = "grey90") +
  labs(title = "Persentase Rumah Tangga Ber-PHBS",
       fill = "%") +
  theme_minimal(); plot_phbs

plot_miskin <- ggplot(jabar) +
  geom_sf(aes(fill = MISKIN), color = "darkgray", size = 0.2) +
  scale_fill_viridis_c(option = "turbo") +
  labs(title = "Penduduk Miskin", fill = "%") +
  theme_minimal(); plot_miskin

grid.arrange(plot_kusta, plot_kepadatan, plot_phbs, plot_miskin, nrow = 2, ncol = 2)

# Plot Antar Variabel
# Scatter Plot KUSTA vs KEPADATAN
ggplot(data_spasial, aes(x = KEPADATAN, y = KUSTA)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Kepadatan dan Kusta",
       x = "Kepadatan Penduduk (Jiwa/km2)",
       y = "Kasus Kusta per 100.000 penduduk") +
  theme_minimal()

# Scatter Plot KUSTA vs PHBS
ggplot(data_spasial, aes(x = PHBS, y = KUSTA)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara PHBS dan Kusta",
       x = "Persentase Rumah Tangga Ber-PHBS",
       y = "Kasus Kusta per 100.000 penduduk") +
  theme_minimal()

# Scatter Plot KUSTA vs Penduduk Miskin
ggplot(data_spasial, aes(x = MISKIN, y = KUSTA)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(title = "Hubungan antara Penduduk Miskin dan Kusta",
       x = "Persentase Penduduk Miskin",
       y = "Kasus Kusta per 100.000 penduduk") +
  theme_minimal()

# DIAGNOSIS ----
# Persiapan Data (Hapus NA)
coords_all <- st_coordinates(st_centroid(st_geometry(jabar)))
jabar$X <- coords_all[, 1]
jabar$Y <- coords_all[, 2]
jabar_ada <- jabar %>% filter(!is.na(KUSTA))
jabar_miss <- jabar %>% filter(is.na(KUSTA))

# ==============================================================================
# 1. DIAGNOSIS NORMALITAS
# ==============================================================================
# Uji Shapiro-Wilk (Data Asli)
shapiro.test(jabar_ada$KUSTA)

## 
##  Shapiro-Wilk normality test
## 
## data:  jabar_ada$KUSTA
## W = 0.79675, p-value = 0.0002008

# ==============================================================================
# 2. DIAGNOSIS TREND SPASIAL
# ==============================================================================
# Model Linear (Orde 1)
model_orde1 <- lm(KUSTA ~ X + Y, data = jabar_ada)
summary(model_orde1)

## 
## Call:
## lm(formula = KUSTA ~ X + Y, data = jabar_ada)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7473 -1.3819  0.0641  1.2216  6.9493 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -287.8703    82.9699  -3.470 0.002177 ** 
## X              3.2032     0.8051   3.979 0.000635 ***
## Y              7.8444     1.3011   6.029 4.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.279 on 22 degrees of freedom
## Multiple R-squared:  0.6345, Adjusted R-squared:  0.6013 
## F-statistic:  19.1 on 2 and 22 DF,  p-value: 1.553e-05

# Model Kuadratik (Orde 2)
model_orde2 <- lm(KUSTA ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), data = jabar_ada)
summary(model_orde2)

## 
## Call:
## lm(formula = KUSTA ~ X + Y + I(X^2) + I(Y^2) + I(X * Y), data = jabar_ada)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6399 -0.4864 -0.0915  0.2624  3.6821 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26392.592  12283.659   2.149 0.044770 *  
## X            -556.417    232.706  -2.391 0.027307 *  
## Y           -1000.517    234.552  -4.266 0.000418 ***
## I(X^2)          2.959      1.110   2.666 0.015275 *  
## I(Y^2)         15.477      3.323   4.657 0.000172 ***
## I(X * Y)       11.338      2.500   4.536 0.000226 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.623 on 19 degrees of freedom
## Multiple R-squared:   0.84,  Adjusted R-squared:  0.7979 
## F-statistic: 19.95 on 5 and 19 DF,  p-value: 5.757e-07

# Perbandingan AIC (Pilih nilai terkecil)
AIC(model_orde1, model_orde2)

##             df      AIC
## model_orde1  4 116.9448
## model_orde2  7 102.2933

# Visualisasi Tren (Metode Proyeksi UTM Zone 48S - Jawa Barat)
jabar_utm <- st_transform(jabar_ada, 32748)
coords_centroid <- st_coordinates(st_centroid(st_geometry(jabar_utm)))
jabar_utm$X_coord <- coords_centroid[,1]
jabar_utm$Y_coord <- coords_centroid[,2]

par(mfrow=c(1,2))
# Tren Barat - Timur
plot(jabar_utm$X_coord, jabar_utm$KUSTA, 
     main="Tren Barat-Timur", 
     xlab="Longitude (X)", ylab="Kusta", 
     pch=19, col="darkblue")
abline(lm(KUSTA ~ X_coord, data=jabar_utm), col="red", lwd=2)

# Tren Utara - Selatan
plot(jabar_utm$Y_coord, jabar_utm$KUSTA, 
     main="Tren Utara-Selatan", 
     xlab="Latitude (Y)", ylab="Kusta", 
     pch=19, col="darkgreen")
abline(lm(KUSTA ~ Y_coord, data=jabar_utm), col="red", lwd=2)

par(mfrow=c(1,1))

# ==============================================================================
# 3. DIAGNOSIS AUTOKORELASI SPASIAL
# ==============================================================================
nb_ada <- poly2nb(jabar_ada, queen = TRUE)
lw_ada <- nb2listw(nb_ada, style = "W", zero.policy = TRUE)
moran.test(jabar_ada$KUSTA, lw_ada, zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  jabar_ada$KUSTA  
## weights: lw_ada  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 2.7942, p-value = 0.002601
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.36051967       -0.04347826        0.02090447

# ==============================================================================
# 4. DIAGNOSIS STASIONERITAS
# ==============================================================================
# Semivariogram Cloud
sp_jabar <- as(jabar_ada, "Spatial")
v_cloud <- variogram(KUSTA ~ 1, sp_jabar, cloud = TRUE)
plot(v_cloud, main = "Semivariogram Cloud", 
     xlab = "Jarak (Distance)", ylab = "Semivariance")

# Uji Levene (Barat vs Timur)
jabar_ada$kelompok <- ifelse(jabar_ada$X > mean(jabar_ada$X), "Timur", "Barat")
leveneTest(KUSTA ~ as.factor(kelompok), data = jabar_ada)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.9207 0.3473
##       23

# Uji Levene (Utara vs Selatan)
jabar_ada$kelompok_lat <- ifelse(jabar_ada$Y > mean(jabar_ada$Y), "Utara", "Selatan")
leveneTest(KUSTA ~ as.factor(kelompok_lat), data = jabar_ada)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1   5.169 0.03265 *
##       23                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ==============================================================================
# 5. DIAGNOSIS ANISOTROPI
# ==============================================================================
v_dir <- variogram(KUSTA ~ 1, sp_jabar, alpha = c(0, 45, 90, 135))
plot(v_dir, main = "Directional Variogram")

# ==============================================================================
# 6. DIAGNOSIS KERAPATAN DATA (HEURISTIK)
# ==============================================================================
# Variogram Eksperimental
v_fit <- variogram(KUSTA ~ 1, sp_jabar)
plot(
  v_fit,
  main = "Variogram Eksperimental",
  xlab = "Jarak (Distance)",
  ylab = "Semivariance",
  pch  = 19,     # titik solid
  cex  = 1     # ukuran titik (makin besar = makin tebal)
)

# Peta Densitas Titik
ggplot(jabar) +
  geom_sf(fill = "white") +
  stat_density_2d(data = jabar_ada, aes(x = X, y = Y, fill = ..level..), 
                  geom = "polygon", alpha = 0.3) +
  geom_point(data = jabar_ada, aes(x = X, y = Y), color = "red", size = 1) +
  scale_fill_viridis_c() +
  labs(title = "Peta Densitas Titik Observasi") +
  theme_minimal()

# Range vs Mean Nearest Neighbor Distance
coords_dist <- cbind(jabar_ada$X, jabar_ada$Y)
dist_vector <- nndist(coords_dist)
mean_nn_dist <- mean(dist_vector)
print(paste("Rata-rata Jarak Antar Kab/Kota:", mean_nn_dist))

## [1] "Rata-rata Jarak Antar Kab/Kota: 0.207987466884834"

# IDW -----
# 1. DEFINISI FUNGSI IDW
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
  d <- sqrt((x - x0)^2 + (y - y0)^2)
  if (any(d == 0)) return(z[which.min(d)])
  
  if (!is.null(k)) {
    k <- min(k, length(d))
    idx <- order(d)[1:k]
    d <- d[idx]; z <- z[idx]
  }
  w <- 1 / (d^p)
  return(sum(w * z, na.rm = TRUE) / sum(w, na.rm = TRUE))
}

# Fungsi Evaluasi (LOOCV) untuk menghitung RMSE & MAE
calc_metrics_idw <- function(p_val, k_val = NULL) {
  n <- nrow(jabar_ada)
  preds <- sapply(1:n, function(i) {
    idw_predict(jabar_ada$X[-i], jabar_ada$Y[-i], jabar_ada$KUSTA[-i], 
                jabar_ada$X[i], jabar_ada$Y[i], p = p_val, k = k_val)
  })
  actual <- jabar_ada$KUSTA
  rmse <- sqrt(mean((actual - preds)^2, na.rm = TRUE))
  mae <- mean(abs(actual - preds), na.rm = TRUE)
  return(c(RMSE = rmse, MAE = mae))
}

# ==============================================================================
# 2. OPTIMASI PARAMETER POWER (P)
# ==============================================================================

uji_p <- data.frame(Power = 1:6)
metrics_p <- t(sapply(uji_p$Power, function(p) calc_metrics_idw(p_val = p)))
uji_p <- cbind(uji_p, metrics_p)
print(uji_p)

##   Power     RMSE      MAE
## 1     1 3.251234 2.463341
## 2     2 2.990373 2.092577
## 3     3 2.851183 1.793995
## 4     4 2.805015 1.667966
## 5     5 2.801238 1.612128
## 6     6 2.816423 1.581781

# ==============================================================================
# 3. OPTIMASI PARAMETER TETANGGA (K)
# ==============================================================================
uji_k <- data.frame(k_val = c(3, 5, 8, 10, 15, 20))
metrics_k <- t(sapply(uji_k$k_val, function(k) calc_metrics_idw(p_val = 5, k_val = k)))
uji_k <- cbind(uji_k, metrics_k)
print(uji_k)

##   k_val     RMSE      MAE
## 1     3 2.865247 1.576177
## 2     5 2.781910 1.604221
## 3     8 2.795545 1.617883
## 4    10 2.791622 1.611995
## 5    15 2.799791 1.612698
## 6    20 2.800517 1.611939

# ==============================================================================
# 4. PREDIKSI FINAL & ANALISIS RESIDUAL
# ==============================================================================
p_opt <- 5
k_opt <- 5

# Imputasi data kosong (NA)
jabar$KUSTA_Final <- jabar$KUSTA

if(nrow(jabar_miss) > 0) {
  pred_miss <- sapply(1:nrow(jabar_miss), function(i) {
    idw_predict(jabar_ada$X, jabar_ada$Y, jabar_ada$KUSTA, 
                jabar_miss$X[i], jabar_miss$Y[i], p = p_opt, k = k_opt)
  })
  jabar$KUSTA_Final[is.na(jabar$KUSTA)] <- pred_miss
  ringkasan <- jabar %>% 
    filter(is.na(KUSTA)) %>% 
    st_drop_geometry() %>%
    dplyr::select(shapeName, KUSTA_Final)
  
  cat("\n--- Daftar Nilai Kusta Hasil Imputasi IDW ---\n")
  print(ringkasan)
}

## 
## --- Daftar Nilai Kusta Hasil Imputasi IDW ---
##     shapeName KUSTA_Final
## 1      Ciamis   0.7332690
## 2 Kota Cimahi   0.2194462

# Diagnostik: Hitung Residual LOOCV pada data yang ada
jabar_ada$pred_loocv <- sapply(1:nrow(jabar_ada), function(i) {
  idw_predict(jabar_ada$X[-i], jabar_ada$Y[-i], jabar_ada$KUSTA[-i], 
              jabar_ada$X[i], jabar_ada$Y[i], p = p_opt, k = k_opt)
})
jabar_ada$residual <- jabar_ada$KUSTA - jabar_ada$pred_loocv

# ==============================================================================
# 5. VISUALISASI PETA RESIDUAL (ERROR SPASIAL)
# ==============================================================================

ggplot() +
  geom_sf(data = jabar, fill = "white", color = "gray80") +
  geom_point(data = jabar_ada, aes(x = X, y = Y, color = residual, size = abs(residual))) +
  scale_color_gradient2(low = "blue", mid = "bisque", high = "red") +
  labs(title = "Peta Residual Interpolasi IDW Kusta",
       color = "Residual", size = "Besar Error") +
  theme_minimal()

# ==============================================================================
# 6. VISUALISASI PERMUKAAN INTERPOLASI (RASTER & CONTOUR)
# ==============================================================================

grid_bounds <- st_bbox(jabar)
xg <- seq(grid_bounds["xmin"] - 0.05, grid_bounds["xmax"] + 0.05, length.out = 100)
yg <- seq(grid_bounds["ymin"] - 0.05, grid_bounds["ymax"] + 0.05, length.out = 100)
G_jabar <- expand.grid(X = xg, Y = yg)
G_jabar$Estimasi <- mapply(function(x0, y0) {
  idw_predict(jabar_ada$X, jabar_ada$Y, jabar_ada$KUSTA, x0, y0, p = p_opt, k = k_opt)
}, G_jabar$X, G_jabar$Y)

ggplot() +
  geom_raster(data = G_jabar, aes(X, Y, fill = Estimasi), interpolate = TRUE) +
  geom_contour(data = G_jabar, aes(X, Y, z = Estimasi), color = "white", alpha = 0.3) +
  geom_sf(data = jabar, fill = NA, color = "white", size = 0.2, alpha = 0.5) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Permukaan Interpolasi IDW Kusta",
       x = "Longitude", y = "Latitude", fill = "Estimasi Kusta") +
  theme_minimal()

# EKMET ----
# OLS
ols_jabar_imp <- lm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar)
summary(ols_jabar_imp)

## 
## Call:
## lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.742 -2.245 -1.500  2.164 10.471 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.912e+00  6.954e+00  -0.706    0.487
## KEPADATAN    4.059e-05  2.045e-04   0.199    0.844
## PHBS         5.774e-02  7.344e-02   0.786    0.440
## MISKIN       4.883e-01  3.738e-01   1.306    0.204
## 
## Residual standard error: 3.571 on 23 degrees of freedom
## Multiple R-squared:  0.1047, Adjusted R-squared:  -0.01211 
## F-statistic: 0.8963 on 3 and 23 DF,  p-value: 0.4581

str(jabar$KUSTA_Final)

##  num [1:27] 0.16 0.41 7.41 4.07 0.733 ...

# Uji Asumsi Model OLS
vif(ols_jabar_imp)

## KEPADATAN      PHBS    MISKIN 
##  1.857542  1.107173  1.998637

bptest(ols_jabar_imp)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_jabar_imp
## BP = 2.4407, df = 3, p-value = 0.4861

resettest(ols_jabar_imp, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  ols_jabar_imp
## RESET = 1.2831, df1 = 2, df2 = 21, p-value = 0.298

res <- residuals(ols_jabar_imp)
ks.test(res, "pnorm", mean(res), sd(res))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res
## D = 0.22952, p-value = 0.09859
## alternative hypothesis: two-sided

# Matriks Bobot Spasial Queen
nb_queen <- poly2nb(jabar, queen = TRUE)
lw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

# Uji Moran's I OLS
moran.test(residuals(ols_jabar_imp), lw_queen)

## 
##  Moran I test under randomisation
## 
## data:  residuals(ols_jabar_imp)  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 3.0603, p-value = 0.001105
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.37239424       -0.03846154        0.01802346

# Autokorelasi Spasial ----
vars_test <- c("KUSTA_Final", "KEPADATAN", "PHBS", "MISKIN")

# Moran's I & Permutasi
for (v in vars_test) {
  cat("\n--- Moran's I & Monte Carlo Test:", v, "---\n")
  print(moran.test(jabar[[v]], lw_queen))
  print(moran.mc(jabar[[v]], lw_queen, nsim = 999))
}

## 
## --- Moran's I & Monte Carlo Test: KUSTA_Final ---
## 
##  Moran I test under randomisation
## 
## data:  jabar[[v]]  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 3.3683, p-value = 0.0003781
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.40185371       -0.03846154        0.01708823 
## 
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar[[v]] 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.40185, observed rank = 997, p-value = 0.003
## alternative hypothesis: greater
## 
## 
## --- Moran's I & Monte Carlo Test: KEPADATAN ---
## 
##  Moran I test under randomisation
## 
## data:  jabar[[v]]  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 2.1464, p-value = 0.01592
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.25844576       -0.03846154        0.01913399 
## 
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar[[v]] 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.25845, observed rank = 972, p-value = 0.028
## alternative hypothesis: greater
## 
## 
## --- Moran's I & Monte Carlo Test: PHBS ---
## 
##  Moran I test under randomisation
## 
## data:  jabar[[v]]  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = -0.40289, p-value = 0.6565
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.09478443       -0.03846154        0.01954318 
## 
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar[[v]] 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = -0.094784, observed rank = 351, p-value = 0.649
## alternative hypothesis: greater
## 
## 
## --- Moran's I & Monte Carlo Test: MISKIN ---
## 
##  Moran I test under randomisation
## 
## data:  jabar[[v]]  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 3.6535, p-value = 0.0001293
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.47755328       -0.03846154        0.01994829 
## 
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar[[v]] 
## weights: lw_queen  
## number of simulations + 1: 1000 
## 
## statistic = 0.47755, observed rank = 997, p-value = 0.003
## alternative hypothesis: greater

# Geary's C
for (v in vars_test) {
  cat("\n--- Geary's C Test:", v, "---\n")
  print(geary.test(jabar[[v]], lw_queen, randomisation = TRUE))
}

## 
## --- Geary's C Test: KUSTA_Final ---
## 
##  Geary C test under randomisation
## 
## data:  jabar[[v]] 
## weights: lw_queen   
## 
## Geary C statistic standard deviate = 2.4056, p-value = 0.008074
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.52095981        1.00000000        0.03965615 
## 
## 
## --- Geary's C Test: KEPADATAN ---
## 
##  Geary C test under randomisation
## 
## data:  jabar[[v]] 
## weights: lw_queen   
## 
## Geary C statistic standard deviate = 2.3889, p-value = 0.008451
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##         0.5936679         1.0000000         0.0289324 
## 
## 
## --- Geary's C Test: PHBS ---
## 
##  Geary C test under randomisation
## 
## data:  jabar[[v]] 
## weights: lw_queen   
## 
## Geary C statistic standard deviate = -0.19959, p-value = 0.5791
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.03266735        1.00000000        0.02678747 
## 
## 
## --- Geary's C Test: MISKIN ---
## 
##  Geary C test under randomisation
## 
## data:  jabar[[v]] 
## weights: lw_queen   
## 
## Geary C statistic standard deviate = 3.7556, p-value = 8.646e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.41019155        1.00000000        0.02466388

# LISA
safe_lisa_perm <- function(y, varname, sf_obj, lw,
                           sig_thr = 0.05, nsim = 999) {
  ok <- is.finite(y)
  z  <- as.numeric(scale(y[ok]))
  lag_z <- lag.listw(lw, z, zero.policy = TRUE)
  lisa <- localmoran_perm(
    z, lw,
    nsim = nsim,
    zero.policy = TRUE
  )
  p_two <- lisa[, "Pr(z != E(Ii))"]
  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")
  sf_sub <- sf_obj[ok, ] %>%
    mutate(
      z_score  = z,
      lag_z    = lag_z,
      Ii       = lisa[, "Ii"],
      p_lisa   = p_two,
      LISA_cat = factor(cluster,
                        levels = c("High-High","Low-Low",
                                   "High-Low","Low-High",
                                   "Not Significant")),
      Variabel = varname)
  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)
}
lisa_results <- map(
  vars_test,
  ~ safe_lisa_perm(jabar[[.x]], .x, jabar, lw_queen))
lisa_tables <- map2_df(
  lisa_results, vars_test,
  ~ mutate(.x$tabel, Variabel = .y))
plots <- map2(lisa_results, vars_test, ~ {
  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(.y)) +
    theme_minimal() +
    theme(axis.text = element_blank(),
          axis.title = element_blank(),
          panel.grid = element_blank())
})
print(lisa_tables)

## # A tibble: 7 × 4
##   LISA_cat        Kabupaten                                      Jumlah Variabel
##   <fct>           <chr>                                           <int> <chr>   
## 1 High-High       Cirebon                                             1 KUSTA_F…
## 2 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciamis…     26 KUSTA_F…
## 3 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciamis…     27 KEPADAT…
## 4 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciamis…     27 PHBS    
## 5 High-High       Cirebon, Majalengka                                 2 MISKIN  
## 6 Low-Low         Kota Bekasi                                         1 MISKIN  
## 7 Not Significant Bandung, Bandung Barat, Bekasi, Bogor, Ciamis…     24 MISKIN

(plots[[1]] | plots[[2]]) / (plots[[3]] | plots[[4]])

# Getis–Ord Gi*
cols_gstar <- c("Hot spot" = "#D7191C", "Cold spot" = "#1B9E77", "Not significant" = "#EAEAEA")
compute_gistar <- function(sf_obj, varname, lw) {
  x <- sf_obj[[varname]]
  Gz <- localG(x, lw, zero.policy = TRUE)
  z_val <- as.numeric(Gz)
  cat_name <- paste0("GSTAR_", varname)
  z_name <- paste0("z_Gi_", varname)
  sf_out <- sf_obj %>%
    mutate(
      !!z_name := z_val,
      !!cat_name := case_when(
        z_val >= 1.96 ~ "Hot spot",
        z_val <= -1.96 ~ "Cold spot",
        TRUE ~ "Not significant"
      ) %>% factor(levels = c("Hot spot", "Cold spot", "Not significant"))
    )
  tab <- sf_out %>% 
    st_drop_geometry() %>%
    transmute(Variabel = varname, Kabupaten = shapeName, 
              Kategori = .data[[cat_name]], 
              Gi_star = round(z_val, 3)) %>%
    filter(Kategori != "Not significant") %>%
    arrange(Kategori, desc(Gi_star))
  list(sf = sf_out, table = tab)
}

# === Eksekusi Analisis ===
gistar_results <- map(vars_test, ~ compute_gistar(jabar, .x, lw_queen))
gistar_tables <- map_dfr(gistar_results, "table")
print(gistar_tables)

##      Variabel   Kabupaten  Kategori Gi_star
## 1 KUSTA_Final     Cirebon  Hot spot   2.653
## 2      MISKIN     Cirebon  Hot spot   2.489
## 3      MISKIN  Majalengka  Hot spot   2.485
## 4      MISKIN       Bogor Cold spot  -2.122
## 5      MISKIN Kota Bekasi Cold spot  -2.459

gistar_plots <- map2(gistar_results, vars_test, ~{
  sf_p <- .x$sf
  fill_col <- paste0("GSTAR_", .y)
  batas_prov <- st_as_sf(st_cast(st_union(sf_p), "MULTILINESTRING"))
  ggplot(sf_p) +
    geom_sf(aes(fill = .data[[fill_col]]), color = "#FFFFFF", size = 0.2) +
    geom_sf(data = st_boundary(sf_p), color = "#6B7280", size = 0.3, fill = NA) +
    geom_sf(data = batas_prov, color = "#374151", size = 0.6, fill = NA) +
    scale_fill_manual(values = cols_gstar, drop = FALSE, name = "Kategori") +
    labs(title = paste0(.y)) +
    theme_minimal() +
    theme(axis.title = element_blank(),
          axis.text = element_blank(),
          panel.grid = element_blank(),
          legend.position = "right")
})

(gistar_plots[[1]] | gistar_plots[[2]]) / (gistar_plots[[3]] | gistar_plots[[4]])

# Uji LM
lm.LMtests(ols_jabar_imp, lw_queen, test = "all")

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data =
## jabar)
## test weights: listw
## 
## RSerr = 6.0408, df = 1, p-value = 0.01398
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data =
## jabar)
## test weights: listw
## 
## RSlag = 6.4478, df = 1, p-value = 0.01111
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data =
## jabar)
## test weights: listw
## 
## adjRSerr = 0.013473, df = 1, p-value = 0.9076
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data =
## jabar)
## test weights: listw
## 
## adjRSlag = 0.42052, df = 1, p-value = 0.5167
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data =
## jabar)
## test weights: listw
## 
## SARMA = 6.4613, df = 2, p-value = 0.03953

# Kesimpulan : SAR, SEM, SDM, SAC

# Modeling ----
# 0. Model OLS 
model_ols <- lm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar)
summary(model_ols)

## 
## Call:
## lm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.742 -2.245 -1.500  2.164 10.471 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.912e+00  6.954e+00  -0.706    0.487
## KEPADATAN    4.059e-05  2.045e-04   0.199    0.844
## PHBS         5.774e-02  7.344e-02   0.786    0.440
## MISKIN       4.883e-01  3.738e-01   1.306    0.204
## 
## Residual standard error: 3.571 on 23 degrees of freedom
## Multiple R-squared:  0.1047, Adjusted R-squared:  -0.01211 
## F-statistic: 0.8963 on 3 and 23 DF,  p-value: 0.4581

AIC(model_ols)

## [1] 151.0209

# 1. Model SLX (X)
model_slx <- lmSLX(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                   data = jabar, 
                   listw = lw_queen,
                   Durbin = TRUE)
summary(model_slx)

## 
## 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)    -1.561e+01   1.182e+01  -1.320e+00   2.016e-01
## KEPADATAN      -6.022e-06   2.185e-04  -2.756e-02   9.783e-01
## PHBS            6.715e-02   8.260e-02   8.129e-01   4.259e-01
## MISKIN          5.487e-01   5.701e-01   9.626e-01   3.473e-01
## lag.KEPADATAN   2.947e-04   5.130e-04   5.744e-01   5.721e-01
## lag.PHBS        1.352e-01   1.430e-01   9.457e-01   3.556e-01
## lag.MISKIN      2.833e-03   7.836e-01   3.616e-03   9.972e-01

# 2. Model SAR (Y)
model_sar <- lagsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                      data = jabar, 
                      listw = lw_queen,
                      method = "eigen")
summary(model_sar)

## 
## Call:lagsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     listw = lw_queen, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.77144 -1.37245 -0.62919  1.11116  9.87840 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -2.9603e+00  5.1643e+00 -0.5732   0.5665
## KEPADATAN   -3.2007e-06  1.5160e-04 -0.0211   0.9832
## PHBS         2.2316e-02  5.4608e-02  0.4087   0.6828
## MISKIN       3.2099e-01  2.8359e-01  1.1319   0.2577
## 
## Rho: 0.64307, LR test value: 8.3787, p-value: 0.0037964
## Asymptotic standard error: 0.14489
##     z-value: 4.4384, p-value: 9.063e-06
## Wald statistic: 19.699, p-value: 9.063e-06
## 
## Log likelihood: -66.32109 for lag model
## ML residual variance (sigma squared): 7.0036, (sigma: 2.6464)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 144.64, (AIC for lm: 151.02)
## LM test for residual autocorrelation
## test value: 1.1886, p-value: 0.27561

# 3. Model SEM (E)
model_sem <- errorsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                        data = jabar, 
                        listw = lw_queen,
                        method = "eigen")
summary(model_sem)

## 
## Call:errorsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
##     data = jabar, listw = lw_queen, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.47100 -1.40098 -0.55661  0.84505  9.93585 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept)  0.28509232  5.98311386  0.0476   0.9620
## KEPADATAN   -0.00004716  0.00015442 -0.3054   0.7601
## PHBS         0.00485236  0.05105844  0.0950   0.9243
## MISKIN       0.36893102  0.34853889  1.0585   0.2898
## 
## Lambda: 0.66112, LR test value: 8.3333, p-value: 0.0038926
## Asymptotic standard error: 0.14328
##     z-value: 4.6141, p-value: 3.948e-06
## Wald statistic: 21.29, p-value: 3.948e-06
## 
## Log likelihood: -66.34381 for error model
## ML residual variance (sigma squared): 6.9516, (sigma: 2.6366)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 144.69, (AIC for lm: 151.02)

# 4. Model SDM (X Y)
model_sdm <- lagsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                      data = jabar, 
                      listw = lw_queen, 
                      Durbin = TRUE, method = "eigen")
summary(model_sdm)

## 
## Call:lagsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     listw = lw_queen, Durbin = TRUE, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.47595 -1.09319 -0.45089  0.47587  9.65647 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   -1.0949e+01  8.2633e+00 -1.3250   0.1852
## KEPADATAN     -2.9806e-05  1.5189e-04 -0.1962   0.8444
## PHBS           3.8449e-02  5.7505e-02  0.6686   0.5037
## MISKIN         4.5654e-01  3.9569e-01  1.1538   0.2486
## lag.KEPADATAN  1.6586e-04  3.5634e-04  0.4655   0.6416
## lag.PHBS       1.0936e-01  9.9350e-02  1.1008   0.2710
## lag.MISKIN    -2.1195e-01  5.4969e-01 -0.3856   0.6998
## 
## Rho: 0.63659, LR test value: 8.2597, p-value: 0.0040535
## Asymptotic standard error: 0.14969
##     z-value: 4.2526, p-value: 2.1127e-05
## Wald statistic: 18.085, p-value: 2.1127e-05
## 
## Log likelihood: -65.38745 for mixed model
## ML residual variance (sigma squared): 6.5563, (sigma: 2.5605)
## Number of observations: 27 
## Number of parameters estimated: 9 
## AIC: 148.77, (AIC for lm: 155.03)
## LM test for residual autocorrelation
## test value: 1.1459, p-value: 0.28441

# 5. Model SDEM (X E)
model_sdem <- errorsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                         data = jabar, 
                         listw = lw_queen, 
                         Durbin = TRUE, method = "eigen")
summary(model_sdem)

## 
## Call:errorsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
##     data = jabar, listw = lw_queen, Durbin = TRUE, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.11214 -1.11924 -0.42298  0.55852  9.44075 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   -1.6469e+01  1.3119e+01 -1.2554   0.2093
## KEPADATAN      2.3356e-05  1.6320e-04  0.1431   0.8862
## PHBS           5.2599e-02  6.3153e-02  0.8329   0.4049
## MISKIN         4.6725e-01  3.9104e-01  1.1949   0.2321
## lag.KEPADATAN  2.2514e-04  3.8711e-04  0.5816   0.5608
## lag.PHBS       1.3192e-01  1.0427e-01  1.2651   0.2058
## lag.MISKIN     3.9227e-01  6.2884e-01  0.6238   0.5328
## 
## Lambda: 0.64918, LR test value: 8.3137, p-value: 0.0039347
## Asymptotic standard error: 0.14659
##     z-value: 4.4286, p-value: 9.4842e-06
## Wald statistic: 19.613, p-value: 9.4842e-06
## 
## Log likelihood: -65.36045 for error model
## ML residual variance (sigma squared): 6.5028, (sigma: 2.5501)
## Number of observations: 27 
## Number of parameters estimated: 9 
## AIC: 148.72, (AIC for lm: 155.03)

# 6. Model SAC / SARAR (Y E)
model_sac <- sacsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                      data = jabar, 
                      listw = lw_queen,
                      method = "eigen")
summary(model_sac)

## 
## Call:sacsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     listw = lw_queen, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.84527 -1.25275 -0.40739  0.92490 10.01172 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -1.4607e+00  5.8474e+00 -0.2498   0.8027
## KEPADATAN   -2.3125e-05  1.5725e-04 -0.1471   0.8831
## PHBS         9.3274e-03  5.4834e-02  0.1701   0.8649
## MISKIN       3.5883e-01  3.9719e-01  0.9034   0.3663
## 
## Rho: 0.40069
## Asymptotic standard error: 1.1492
##     z-value: 0.34868, p-value: 0.72733
## Lambda: 0.40406
## Asymptotic standard error: 1.1633
##     z-value: 0.34735, p-value: 0.72833
## 
## LR test value: 9.1646, p-value: 0.010231
## 
## Log likelihood: -65.92814 for sac model
## ML residual variance (sigma squared): 7.082, (sigma: 2.6612)
## Number of observations: 27 
## Number of parameters estimated: 7 
## AIC: 145.86, (AIC for lm: 151.02)

# 7. Model GNS (X Y E)
model_gns <- sacsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, 
                      data = jabar, 
                      listw = lw_queen, 
                      Durbin = TRUE, method = "eigen")
summary(model_gns)

## 
## Call:sacsarlm(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     listw = lw_queen, Durbin = TRUE, method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.29824 -1.09217 -0.51047  0.42367  9.70860 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value Pr(>|z|)
## (Intercept)   -1.3038e+01  2.3809e+01 -0.5476   0.5840
## KEPADATAN     -4.8836e-06  2.4752e-04 -0.0197   0.9843
## PHBS           4.2669e-02  1.0850e-01  0.3933   0.6941
## MISKIN         4.3629e-01  4.6695e-01  0.9343   0.3501
## lag.KEPADATAN  2.0107e-04  4.2127e-04  0.4773   0.6332
## lag.PHBS       1.1492e-01  1.6185e-01  0.7100   0.4777
## lag.MISKIN     7.4216e-02  1.7472e+00  0.0425   0.9661
## 
## Rho: 0.37824
## Asymptotic standard error: 2.7921
##     z-value: 0.13547, p-value: 0.89224
## Lambda: 0.40565
## Asymptotic standard error: 2.7359
##     z-value: 0.14827, p-value: 0.88213
## 
## LR test value: 10.926, p-value: 0.052862
## 
## Log likelihood: -65.04732 for sacmixed model
## ML residual variance (sigma squared): 6.6661, (sigma: 2.5819)
## Number of observations: 27 
## Number of parameters estimated: 10 
## AIC: 150.09, (AIC for lm: 151.02)

perbandingan_ic <- data.frame(
  Model = c("OLS", "SLX", "SAR", "SEM", "SDM", "SDEM", "SAC", "GNS"),
  AIC = c(AIC(model_ols), AIC(model_slx), AIC(model_sar),
          AIC(model_sem), AIC(model_sdm), AIC(model_sdem),
          AIC(model_sac), AIC(model_gns)),
  BIC = c(BIC(model_ols), BIC(model_slx), BIC(model_sar),
          BIC(model_sem), BIC(model_sdm), BIC(model_sdem),
          BIC(model_sac), BIC(model_gns)))
perbandingan_ic <- perbandingan_ic[order(perbandingan_ic$AIC), ]
print(perbandingan_ic)

##   Model      AIC      BIC
## 3   SAR 144.6422 152.4172
## 4   SEM 144.6876 152.4626
## 7   SAC 145.8563 154.9271
## 6  SDEM 148.7209 160.3834
## 5   SDM 148.7749 160.4374
## 8   GNS 150.0946 163.0530
## 1   OLS 151.0209 157.5001
## 2   SLX 155.0346 165.4013

# Uji Asumsi Model SAR
# Autokorelasi Residual
res_sar <- residuals(model_sar)
moran.test(res_sar, lw_queen, zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  res_sar  
## weights: lw_queen    
## 
## Moran I statistic standard deviate = 0.87629, p-value = 0.1904
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.06977038       -0.03846154        0.01525520

# Normalitas Residual
ks.test(res_sar, "pnorm", mean(res_sar), sd(res_sar))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  res_sar
## D = 0.14986, p-value = 0.5304
## alternative hypothesis: two-sided

# Homoskedastisitas 
x_vars <- model_sar$X[, -1, drop = FALSE]
df_bp <- data.frame(res_sar = res_sar, as.data.frame(x_vars))
bptest(res_sar ~ ., data = df_bp)

## 
##  studentized Breusch-Pagan test
## 
## data:  res_sar ~ .
## BP = 3.0536, df = 3, p-value = 0.3834

plot(fitted(model_sar), residuals(model_sar),
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Uji Visual Homoskedastisitas (SAR)")
abline(h = 0, col = "red")

# Uji RESET
sar_residuals <- residuals(model_sar)
sar_fitted <- fitted(model_sar)
jabar$sar_fitted <- sar_fitted
jabar$sar_fitted_sq <- sar_fitted^2
jabar$sar_fitted_cub <- sar_fitted^3
reset_model <- lm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN + sar_fitted_sq + sar_fitted_cub, data = jabar)
resettest(reset_model, power = 2:3)

## 
##  RESET test
## 
## data:  reset_model
## RESET = 0.81363, df1 = 2, df2 = 19, p-value = 0.4581

# Uji Sensitivitas Bobot Model SAR
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 <- st_coordinates(st_centroid(jabar))
nb_knn4 <- knn2nb(knearneigh(coords, k = 4))
lw_knn4 <- nb2listw(nb_knn4, style = "W")
sar_queen <- lagsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, listw = lw_queen)
sar_rook  <- lagsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, listw = lw_rook)
sar_knn4  <- lagsarlm(KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, listw = lw_knn4)
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.6430721 144.6422 -66.32109
## 2  Rook 0.6430721 144.6422 -66.32109
## 3 KNN-4 0.5890997 145.4894 -66.74471

# Impacts
imp_sar <- impacts(model_sar,listw = lw_queen, R = 500, zero.policy = TRUE)
summary(imp_sar, zstats = TRUE, short = TRUE)

## Impact measures (lag, exact):
##                        Direct      Indirect         Total
## KEPADATAN dy/dx -3.707143e-06 -5.260350e-06 -8.967493e-06
## PHBS dy/dx       2.584696e-02  3.667624e-02  6.252320e-02
## MISKIN dy/dx     3.717744e-01  5.275393e-01  8.993137e-01
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                      Direct     Indirect       Total
## KEPADATAN dy/dx 0.000193101 0.0008351966 0.000972386
## PHBS dy/dx      0.067631662 0.3516171845 0.395504844
## MISKIN dy/dx    0.349815227 2.0311320206 2.231459220
## 
## Simulated z-values:
##                      Direct    Indirect       Total
## KEPADATAN dy/dx 0.008651646 -0.04147131 -0.03390223
## PHBS dy/dx      0.350189440  0.11849065  0.16522488
## MISKIN dy/dx    1.118191221  0.35220252  0.49587738
## 
## Simulated p-values:
##                 Direct  Indirect Total  
## KEPADATAN dy/dx 0.99310 0.96692  0.97296
## PHBS dy/dx      0.72620 0.90568  0.86877
## MISKIN dy/dx    0.26349 0.72469  0.61998

# Pseudo-R²
y <- jabar$KUSTA_Final
R2_sar <- 1 - (model_sar$s2 / var(y)) ; R2_sar

## [1] 0.4440053

# Prediksi (fitted values)
jabar$pred_sar <- as.numeric(fitted(model_sar))

# Residual
jabar$resid_sar <- as.numeric(residuals(model_sar))

p_pred <- ggplot(jabar) +
  geom_sf(aes(fill = pred_sar), color = "white", size = 0.2) +
  scale_fill_viridis_c(
    name = "Prediksi Kusta\n per 100.000 Penduduk",
    option = "plasma") +
  labs(
    title = "Peta Prediksi") +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank())

p_resid <- ggplot(jabar) +
  geom_sf(aes(fill = resid_sar), color = "white", size = 0.2) +
  scale_fill_gradient2(
    name = "Residual",
    low = "#2c7bb6",
    mid = "white",
    high = "#d7191c",
    midpoint = 0) +
  labs(
    title = "Peta Residual") +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank())

(p_pred | p_resid)

# GWR
# Centroid poligon
coords_jabar <- st_coordinates(st_centroid(jabar))

# Bandwidth Optimal (Adaptive)
bandwidth_optimal <- gwr.sel(
  KUSTA_Final ~ KEPADATAN + PHBS + MISKIN,
  data = jabar,
  coords = coords_jabar,
  adapt = TRUE)

## Adaptive q: 0.381966 CV score: 289.4378 
## Adaptive q: 0.618034 CV score: 323.9785 
## Adaptive q: 0.236068 CV score: 254.0271 
## Adaptive q: 0.145898 CV score: 250.8523 
## Adaptive q: 0.170955 CV score: 244.6338 
## Adaptive q: 0.1869381 CV score: 242.8551 
## Adaptive q: 0.1956288 CV score: 243.4253 
## Adaptive q: 0.1867075 CV score: 242.8574 
## Adaptive q: 0.1874272 CV score: 242.8534 
## Adaptive q: 0.1873547 CV score: 242.8534 
## Adaptive q: 0.187314 CV score: 242.8534 
## Adaptive q: 0.1873547 CV score: 242.8534

# Model GWR ----
model_gwr <- gwr(
  KUSTA_Final ~ KEPADATAN + PHBS + MISKIN,
  data = jabar,
  coords = coords_jabar,
  adapt = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit = TRUE)
summary(model_gwr)

##           Length Class                  Mode     
## SDF        27    SpatialPointsDataFrame S4       
## lhat      729    -none-                 numeric  
## lm         11    -none-                 list     
## results    14    -none-                 list     
## bandwidth  27    -none-                 numeric  
## adapt       1    -none-                 numeric  
## hatmatrix   1    -none-                 logical  
## gweight     1    -none-                 character
## gTSS        1    -none-                 numeric  
## this.call   7    -none-                 call     
## fp.given    1    -none-                 logical  
## timings    12    -none-                 numeric

print(model_gwr)

## Call:
## gwr(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     coords = coords_jabar, adapt = bandwidth_optimal, hatmatrix = TRUE, 
##     se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.1873547 (about 5 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -2.1435e+01 -1.5809e+01 -7.4085e+00  3.7827e+00  1.0397e+01
## KEPADATAN    -4.8161e-04 -3.0132e-04 -4.4672e-05  2.3995e-04  3.3970e-04
## PHBS          1.3572e-02  3.0588e-02  7.6660e-02  1.0392e-01  1.5215e-01
## MISKIN       -7.4677e-01 -3.5952e-01  4.4109e-01  1.1268e+00  1.5519e+00
##               Global
## X.Intercept. -4.9118
## KEPADATAN     0.0000
## PHBS          0.0577
## MISKIN        0.4883
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 13.53502 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 13.46498 
## Sigma (residual: 2traceS - traceS'S): 2.727213 
## Effective number of parameters (model: traceS): 10.78312 
## Effective degrees of freedom (model: traceS): 16.21688 
## Sigma (model: traceS): 2.485068 
## Sigma (ML): 1.925928 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 156.7706 
## AIC (GWR p. 96, eq. 4.22): 122.7978 
## Residual sum of squares: 100.1483 
## Quasi-global R2: 0.6942126

# uji t
jabar$t_kepadatan <- model_gwr$SDF$KEPADATAN / model_gwr$SDF$KEPADATAN_se
sum(abs(jabar$t_kepadatan) > 1.96)

## [1] 1

jabar$t_phbs <- model_gwr$SDF$PHBS / model_gwr$SDF$PHBS_se
sum(abs(jabar$t_phbs) > 1.96)

## [1] 0

jabar$t_miskin <- model_gwr$SDF$MISKIN / model_gwr$SDF$MISKIN_se
sum(abs(jabar$t_miskin) > 1.96)

## [1] 12

# Status signifikansi (|t| > 1.96)
jabar$sig_kepadatan <- ifelse(abs(jabar$t_kepadatan) > 1.96,
                              "Signifikan", "Tidak Signifikan")

jabar$sig_phbs <- ifelse(abs(jabar$t_phbs) > 1.96,
                         "Signifikan", "Tidak Signifikan")

jabar$sig_miskin <- ifelse(abs(jabar$t_miskin) > 1.96,
                           "Signifikan", "Tidak Signifikan")

# Plot Miskin yang signifikan
plot_sig <- function(data, var_sig, judul) {
  ggplot(data) +
    geom_sf(aes(fill = .data[[var_sig]]),
            color = "white", size = 0.2) +
    scale_fill_manual(
      values = c("Signifikan" = "#d73027",
                 "Tidak Signifikan" = "#eeeeee"),
      drop = FALSE
    ) +
    labs(title = judul, fill = "Status") +
    theme_minimal() +
    theme(
      axis.text = element_blank(),
      axis.title = element_blank(),
      panel.grid = element_blank()
    )
}

p_kepadatan <- plot_sig(jabar, "sig_kepadatan", "Kepadatan Penduduk")
p_phbs      <- plot_sig(jabar, "sig_phbs","PHBS")
p_miskin    <- plot_sig(jabar, "sig_miskin","Penduduk Miskin")

(p_kepadatan | p_phbs) / (p_miskin   | plot_spacer())

# Visualisasi
df_gwr <- as.data.frame(model_gwr$SDF)
coef_gwr <- df_gwr %>%
  dplyr::select(
    Intercept      = `X.Intercept.`, 
    coef_Kepadatan = KEPADATAN,
    coef_PHBS      = PHBS,
    coef_Miskin    = MISKIN,
    R2_Lokal       = localR2,
    Prediksi_GWR   = pred
  )
jabar_gwr <- cbind(jabar, coef_gwr)
koefmiskin = ggplot(jabar_gwr) +
  geom_sf(aes(fill = coef_Miskin)) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Penduduk Miskin",
       fill = "Koefisien") +
  theme_minimal()
koefpadat = ggplot(jabar_gwr) +
  geom_sf(aes(fill = coef_Kepadatan)) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Kepadatan Penduduk",
       fill = "Koefisien") +
  theme_minimal()
koefphbs = ggplot(jabar_gwr) +
  geom_sf(aes(fill = coef_PHBS)) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "PHBS",
       fill = "Koefisien") +
  theme_minimal()


(koefmiskin | koefpadat) / (koefphbs   | plot_spacer())

# Uji Asumsi GWR
# Heterogeneity
bptest(ols_jabar_imp)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_jabar_imp
## BP = 2.4407, df = 3, p-value = 0.4861

# AIC
AIC(ols_jabar_imp)

## [1] 151.0209

print(model_gwr)

## Call:
## gwr(formula = KUSTA_Final ~ KEPADATAN + PHBS + MISKIN, data = jabar, 
##     coords = coords_jabar, adapt = bandwidth_optimal, hatmatrix = TRUE, 
##     se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.1873547 (about 5 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -2.1435e+01 -1.5809e+01 -7.4085e+00  3.7827e+00  1.0397e+01
## KEPADATAN    -4.8161e-04 -3.0132e-04 -4.4672e-05  2.3995e-04  3.3970e-04
## PHBS          1.3572e-02  3.0588e-02  7.6660e-02  1.0392e-01  1.5215e-01
## MISKIN       -7.4677e-01 -3.5952e-01  4.4109e-01  1.1268e+00  1.5519e+00
##               Global
## X.Intercept. -4.9118
## KEPADATAN     0.0000
## PHBS          0.0577
## MISKIN        0.4883
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 13.53502 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 13.46498 
## Sigma (residual: 2traceS - traceS'S): 2.727213 
## Effective number of parameters (model: traceS): 10.78312 
## Effective degrees of freedom (model: traceS): 16.21688 
## Sigma (model: traceS): 2.485068 
## Sigma (ML): 1.925928 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 156.7706 
## AIC (GWR p. 96, eq. 4.22): 122.7978 
## Residual sum of squares: 100.1483 
## Quasi-global R2: 0.6942126

# Multikolinearity
vif(ols_jabar_imp)

## KEPADATAN      PHBS    MISKIN 
##  1.857542  1.107173  1.998637

# Spatial Autokorelasi
residual_gwr <- model_gwr$SDF$gwr.e
nb <- knn2nb(knearneigh(coords_jabar, k=4))
listw <- nb2listw(nb, style="W")
moran.test(residual_gwr, listw)

## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw    
## 
## Moran I statistic standard deviate = 0.67254, p-value = 0.2506
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.03360239       -0.03846154        0.01148136

# Normality
ks.test(residual_gwr, "pnorm", mean(residual_gwr), sd(residual_gwr))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  residual_gwr
## D = 0.15967, p-value = 0.4504
## alternative hypothesis: two-sided

# Ekstrak koefisien lokal untuk variabel X
coef_x1 <- model_gwr$SDF$KEPADATAN
coef_x2 <- model_gwr$SDF$PHBS
coef_x3 <- model_gwr$SDF$MISKIN

# Plot koefisien regresi lokal untuk variabel X
spplot(model_gwr$SDF, "KEPADATAN", main = "Koefisien GWR untuk x1", col.regions = terrain.colors(100))

spplot(model_gwr$SDF, "PHBS", main = "Koefisien GWR untuk x2", col.regions = terrain.colors(100))

spplot(model_gwr$SDF, "MISKIN", main = "Koefisien GWR untuk x3", col.regions = terrain.colors(100))

# Konversi SpatialPointsDataFrame menjadi data frame biasa
gwr_df <- as.data.frame(model_gwr$SDF)

# Plot koefisien lokal (titik)
ggplot(gwr_df, aes(x = X, y = Y)) +
  geom_point(aes(color = KEPADATAN), size = 2) +
  scale_color_viridis_c() +
  labs(title = "Koefisien GWR untuk x1", color = "Koefisien") +
  theme_minimal()

ggplot(gwr_df, aes(x = X, y = Y)) +
  geom_point(aes(color = PHBS), size = 2) +
  scale_color_viridis_c() +
  labs(title = "Koefisien GWR untuk x2", color = "Koefisien") +
  theme_minimal()

ggplot(gwr_df, aes(x = X, y = Y)) +
  geom_point(aes(color = MISKIN), size = 2) +
  scale_color_viridis_c() +
  labs(title = "Koefisien GWR untuk x3", color = "Koefisien") +
  theme_minimal()

# Plot Peta Local R-Square
jabar$local_R2 <- gwr_df$localR2
localr <- ggplot(jabar) +
  geom_sf(aes(fill = local_R2), color = "white", size = 0.2) +
  scale_fill_viridis_c(option = "mako", name = expression(R^2)) +
  labs(title = "Peta Local R-Square GWR") +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank())


# Peta Residual
jabar$res_gwr <- model_gwr$SDF$gwr.e
residgwr <- ggplot(jabar) +
  geom_sf(aes(fill = res_gwr), color = "white", size = 0.2) +
  scale_fill_gradient2(
    low = "#2166ac",
    mid = "#f7f7f7",
    high = "#b2182b",
    midpoint = 0,
    name = "Residual") +
  labs(title = "Peta Distribusi Residual GWR") +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank())

localr | residgwr

# EVALUASI AKURASI MODEL GWR
actual   <- jabar$KUSTA_Final
predicted <- coef_gwr$Prediksi_GWR
mae_gwr <- mean(abs(actual - predicted))
rmse_gwr <- sqrt(mean((actual - predicted)^2))
mae_gwr

## [1] 1.318003

rmse_gwr

## [1] 1.925928

Script Dashboard

Berikut script dashboard yang digunakan untuk membentuk https://spatial-leprosy.streamlit.app/

import streamlit as st
import pandas as pd
import geopandas as gpd
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio 
from scipy.spatial import distance_matrix
from libpysal.weights import Queen, Rook, lag_spatial
from esda.moran import Moran, Moran_Local
from esda.geary import Geary
from esda.getisord import G_Local
from spreg import OLS, ML_Lag, ML_Error, GM_Lag, GM_Error
from mgwr.gwr import GWR
from mgwr.sel_bw import Sel_BW
from fpdf import FPDF
import os
import joblib
from scipy.stats import norm, chi2, shapiro
from streamlit_option_menu import option_menu
import importlib

# load modules explicitly
spreg_out = importlib.import_module("spreg.output")
ml_error_mod = importlib.import_module("spreg.ml_error")

_old_nonspat_top = spreg_out._nonspat_top

def _patched_nonspat_top(reg, ml=False):
    # reg.sig2 kadang array di Py3.13+numpy2.x; paksa jadi float
    try:
        reg.sig2 = float(np.asarray(reg.sig2).squeeze())
    except Exception:
        pass
    return _old_nonspat_top(reg, ml=ml)

# Patch both: the original module + the copied reference in ml_error module
spreg_out._nonspat_top = _patched_nonspat_top
ml_error_mod._nonspat_top = _patched_nonspat_top

if 'run_ekmet' not in st.session_state:
    st.session_state['run_ekmet'] = False

if 'ekmet_models' not in st.session_state:
    st.session_state['ekmet_models'] = None


# --- KONFIGURASI HALAMAN ---
st.set_page_config(page_title="UAS Spasial - Dashboard Kusta", layout="wide")

# --- FUNGSI UTILITAS ---

def impute_idw(gdf, target_col, p=2):
    """Mengisi nilai NaN menggunakan IDW."""
    if 'x' not in gdf.columns:
        gdf['x'] = gdf.geometry.centroid.x
        gdf['y'] = gdf.geometry.centroid.y
    
    known = gdf[gdf[target_col].notna()].copy()
    unknown = gdf[gdf[target_col].isna()].copy()
    
    if unknown.empty:
        return gdf
    
    dist_mtx = distance_matrix(unknown[['x', 'y']], known[['x', 'y']])
    dist_mtx[dist_mtx == 0] = 1e-10
    
    weights = 1.0 / (dist_mtx ** p)
    weighted_sum = np.sum(weights * known[target_col].values, axis=1)
    sum_weights = np.sum(weights, axis=1)
    imputed_values = weighted_sum / sum_weights
    
    gdf.loc[gdf[target_col].isna(), target_col] = imputed_values
    return gdf

def calculate_idw_loocv(gdf, target_col, p=2, k=5):
    """Menghitung RMSE dan MAE menggunakan LOOCV."""
    coords = np.array(list(zip(gdf.geometry.centroid.x, gdf.geometry.centroid.y)))
    values = gdf[target_col].values
    n = len(gdf)
    
    predictions = []
    
    for i in range(n):
        train_coords = np.delete(coords, i, axis=0)
        train_values = np.delete(values, i)
        test_coord = coords[i].reshape(1, -1)
        
        dists = distance_matrix(test_coord, train_coords)[0]
        
        if k is not None and k < len(dists):
            idx_k = np.argsort(dists)[:k]
            d_k = dists[idx_k]
            z_k = train_values[idx_k]
        else:
            d_k = dists
            z_k = train_values
            
        d_k[d_k == 0] = 1e-10
        weights = 1 / (d_k ** p)
        pred = np.sum(weights * z_k) / np.sum(weights)
        predictions.append(pred)
        
    predictions = np.array(predictions)
    residuals = values - predictions
    mse = np.mean(residuals**2)
    rmse = np.sqrt(mse)
    mae = np.mean(np.abs(residuals))
    
    result_df = gdf.copy()
    result_df['Predicted'] = predictions
    result_df['Residual'] = residuals
    
    return rmse, mae, result_df

# Fungsi PDF Generator
def create_pdf_report(model_name, summary_text):
    pdf = FPDF()
    pdf.add_page()
    pdf.set_font("Arial", 'B', 16)
    pdf.cell(0, 10, txt=f"Laporan Analisis Spasial: {model_name}", ln=1, align='C')
    
    pdf.set_font("Arial", size=10)
    pdf.ln(10)
    pdf.cell(0, 10, txt="Ringkasan Statistik Model:", ln=1, align='L')
    
    # Gunakan font Monospace agar tabel di summary rapi
    pdf.set_font("Courier", size=8) 
    # Encoding handling
    clean_summary = summary_text.encode('latin-1', 'replace').decode('latin-1')
    pdf.multi_cell(0, 5, txt=clean_summary)
    
    return pdf.output(dest='S').encode('latin-1')

@st.cache_data
def load_data():
    """Load dan bersihkan data GeoJSON + Excel."""
    # PENGATURAN GEOJSON (Solusi untuk GitHub)
    file_path = "peta_jabar_kecil.geojson"
    
    if not os.path.exists(file_path):
        st.error(f"⚠️ File '{file_path}' tidak ditemukan! Pastikan Anda sudah menjalankan script penyederhanaan peta dan mengupload file .geojson tersebut ke GitHub.")
        st.stop()

    try:
        gdf = gpd.read_file(file_path)
    except Exception as e:
        st.error(f"Gagal membaca file peta: {e}")
        st.stop()
        
    gdf['shapeName_clean'] = gdf['shapeName'].str.strip().str.lower()
    
    # Filter Jawa Barat (Tetap dilakukan untuk memastikan data bersih)
    kab_jabar = [
        "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"
    ]
    jabar_gdf = gdf[gdf['shapeName_clean'].isin(kab_jabar)].copy()
    
    # Load Excel
    try:
        df_attr = pd.read_excel("DATA SPASIAL.xlsx")
    except Exception as e:
        st.error(f"Gagal membaca Excel: {e}")
        st.stop()
        
    kab_col_candidates = [c for c in df_attr.columns if "KAB" in c.upper()]
    if not kab_col_candidates:
        st.error("Kolom 'KABUPATEN' tidak ditemukan di Excel.")
        st.stop()
    
    kab_col_name = kab_col_candidates[0]
    df_attr['kab_clean'] = df_attr[kab_col_name].astype(str).str.strip().str.lower()
    
    merged = jabar_gdf.merge(df_attr, left_on='shapeName_clean', right_on='kab_clean')
    
    # Koordinat Centroid untuk Analisis
    merged['centroid_x'] = merged.geometry.centroid.x
    merged['centroid_y'] = merged.geometry.centroid.y
    merged['x'] = merged.geometry.centroid.x
    merged['y'] = merged.geometry.centroid.y

    # Imputasi NaN
    merged.loc[merged['shapeName_clean'] == 'kota cimahi', 'KUSTA'] = merged.loc[merged['shapeName_clean'] == 'kota cimahi', 'KUSTA'].fillna(0.219)
    merged.loc[merged['shapeName_clean'] == 'ciamis', 'KUSTA'] = merged.loc[merged['shapeName_clean'] == 'ciamis', 'KUSTA'].fillna(0.733)
    nan_count = merged['KUSTA'].isna().sum()
    
    return merged, nan_count

# --- LOAD DATA ---
data, nan_count_awal = load_data()
vars_list = ['KUSTA', 'KEPADATAN', 'PHBS', 'MISKIN']

# --- SIDEBAR NAVIGASI ---
st.sidebar.title("Analysis Dashboard")
menu = st.sidebar.radio(
    "Pilih halaman:",
    [
        "📊 Data and Methods",
        "🔍 Exploratory Data Analysis",
        "🌍 Spatial Autocorrelation",
        "🗺️ Spatial Interpolation",
        "📐 Spatial Econometrics",
        "🧭 GWR"
    ]
)


st.sidebar.markdown("---")
st.sidebar.subheader("Pengaturan Global")

# 1. Matriks Bobot
weight_type = st.sidebar.selectbox("Matriks Bobot Spasial", ["Queen", "Rook"])

import plotly.io as pio
import streamlit as st

# 2. Theme Selector (Default: Dark)
theme_mode = st.sidebar.selectbox(
    "Theme",
    ["Dark Mode", "Light Mode"],
    index=0  
)

# Theme Logic
if theme_mode == "Dark Mode":
    # Plotly Dark
    pio.templates.default = "plotly_dark"
    st.markdown(
        """
        <style>
        /* Main App */
        .stApp {
            background-color: #0E1117;
            color: #FAFAFA;
        }

        /* Sidebar */
        [data-testid="stSidebar"] {
            background-color: #262730;
        }

        /* Text */
        .stMarkdown, p, h1, h2, h3, h4, h5, h6, label, span {
            color: #FAFAFA !important;
        }

        /* Table / DataFrame */
        div[data-testid="stTable"] {
            color: #FAFAFA;
        }
        </style>
        """,
        unsafe_allow_html=True
    )

else:
    # Plotly Light
    pio.templates.default = "plotly_white"
    st.markdown(
        """
        <style>
        .stApp {
            background-color: #FFFFFF;
            color: #000000;
        }

        [data-testid="stSidebar"] {
            background-color: #F0F2F6;
        }

        .stMarkdown, p, h1, h2, h3, h4, h5, h6, label, span {
            color: #000000 !important;
        }
        </style>
        """,
        unsafe_allow_html=True
    )

# Inisialisasi Matriks Bobot berdasarkan pilihan
if weight_type == "Queen":
    w = Queen.from_dataframe(data)
else:
    w = Rook.from_dataframe(data)
w.transform = 'R'

# --- HEADER UTAMA ---
st.title("🩺 Analisis Spasial Penemuan Kasus Baru Kusta")
st.caption(f"**Oleh: Humaira Zeanova - 140610230030**")

# ==============================================================================
# 0. DATA AND METHODS
# ==============================================================================
if menu == "📊 Data and Methods":
    st.header("📊 Data and Methods")

    col_info, col_map = st.columns([1.5, 1])
    
    with col_info:
        st.subheader("Deskripsi Dashboard")
        st.write("""
        Dashboard ini dikembangkan sebagai alat bantu analisis spasial untuk memetakan dan menganalisis faktor-faktor yang mempengaruhi 
        angka penemuan kasus baru kusta di Provinsi Jawa Barat. Sistem ini mengintegrasikan berbagai teknik geostatistik 
        dan ekonometrika spasial untuk memberikan wawasan yang komprehensif bagi pembuat kebijakan kesehatan.
        """)
        
        st.subheader("Tujuan Analisis")
        st.markdown("""
        1. **Menganalisis Pola Spasial**: Mengidentifikasi apakah kasus kusta tersebar secara acak atau membentuk pola pengelompokan (clustering).
        2. **Estimasi Nilai**: Menggunakan teknik interpolasi untuk memprediksi angka kusta di wilayah yang memiliki keterbatasan data.
        3. **Identifikasi Faktor Risiko**: Memodelkan hubungan antara kepadatan penduduk, PHBS, dan kemiskinan terhadap angka kusta baik secara global maupun lokal.
        """)

    with col_map:
        st.subheader("Peta Administrasi Jawa Barat")
        fig_admin = px.choropleth(
            data, geojson=data.geometry, locations=data.index,
            hover_name="shapeName", color_discrete_sequence=["#e0e0e0"]
        )
        fig_admin.update_geos(fitbounds="locations", visible=False)
        fig_admin.update_layout(margin={"r":0,"t":0,"l":0,"b":0}, height=300)
        st.plotly_chart(fig_admin, use_container_width=True)
        st.caption("Cakupan wilayah penelitian: 27 Kabupaten/Kota di Jawa Barat Tahun 2024.")

    st.markdown("---")
    
    st.subheader("Metode Analisis")
    m_col1, m_col2, m_col3 = st.columns(3)
    with m_col1:
        st.info("**Interpolasi (IDW)**\nMengestimasi nilai di lokasi tanpa sampel berdasarkan rata-rata tertimbang dari titik-titik di sekitarnya.")
    with m_col2:
        st.info("**Autokorelasi Spasial**\nMenggunakan Moran's I, Geary's C, LISA dan Getis Ord* untuk mendeteksi korelasi nilai antar lokasi yang bertetangga.")
    with m_col3:
        st.info("**Ekonometrika Spasial**\nModel SLX, SAR, SEM, SDM, SDEM, SAC, dan GNS yang memperhitungkan efek ketergantungan spasial dalam data.")

    st.markdown("---")
    
    st.subheader("Definisi Variabel & Data")
    
    with st.expander("📌 Lihat Definisi Variabel"):
        st.markdown("""
        | Variabel | Definisi Operasional | Satuan |
        | :--- | :--- | :--- |
        | **KUSTA** | Angka penemuan kasus baru kusta (New Case Detection Rate). | Per 100.000 Penduduk |
        | **KEPADATAN** | Jumlah penduduk per unit luas wilayah. | Jiwa / km² |
        | **PHBS** | Rumah tangga yang berperilaku hidup bersih dan sehat. | Persentase (%) |
        | **MISKIN** | Penduduk dengan rata-rata pengeluaran per bulan di bawah garis kemiskinan. | Persentase (%) |
        """)
    
    # Tabel Data
    display_cols = ['shapeName'] + vars_list
    st.write("**Tabel Data yang Digunakan:**")
    st.dataframe(data[display_cols].rename(columns={'shapeName': 'Kabupaten/Kota'}), use_container_width=True)
    st.warning(f"⚠️ Terdapat 2 data Kusta NaN (Ciamis dan Cimahi) yang diimputasi dengan IDW.")
    st.markdown("**Sumber Data:** *Badan Pusat Statistik (BPS) Provinsi Jawa Barat.*")

# ==============================================================================
# 1. EXPLORATORY DATA ANALYSIS (EDA)
# ==============================================================================
elif menu == "🔍 Exploratory Data Analysis":
    st.header("🔍 Exploratory Data Analysis")
    
    # Pemilihan variabel (Multiselect)
    selected_vars = st.multiselect("Pilih Variabel Eksplorasi:", vars_list, default=[vars_list[0]])
    
    if not selected_vars:
        st.error("⚠️ Pilih minimal satu variabel untuk memulai eksplorasi.")
    else:
        # 1. TABEL STATISTIKA DESKRIPTIF (WIDE / MEMANJANG)
        st.subheader("Statistika Deskriptif")
        desc_df = data[selected_vars].describe().T
        desc_df['variance'] = data[selected_vars].var()
        desc_df = desc_df[['mean', 'std', 'variance', 'min', '25%', '50%', '75%', 'max']]
        desc_df.columns = ['Mean', 'Std', 'Var', 'Min', 'Q1', 'Median', 'Q3', 'Max']
        st.dataframe(desc_df, use_container_width=True)

        st.markdown("---")
        
        # 2. VISUALISASI DISTRIBUSI (TABS)
        tab_hist, tab_box, tab_map = st.tabs(["Histogram", "Boxplot", "Peta Sebaran Data Asli"])
        
        with tab_hist:
            st.subheader("Histogram")
            rows_hist = (len(selected_vars) + 1) // 2
            for i in range(rows_hist):
                cols = st.columns(2)
                for j in range(2):
                    idx = i * 2 + j
                    if idx < len(selected_vars):
                        var = selected_vars[idx]
                        fig_hist = px.histogram(data, x=var, nbins=15, title=f"Histogram {var}", 
                                                color_discrete_sequence=['skyblue'], template=pio.templates.default)
                        cols[j].plotly_chart(fig_hist, use_container_width=True)

        with tab_box:
            st.subheader("Boxplot")
            rows_box = (len(selected_vars) + 1) // 2
            for i in range(rows_box):
                cols = st.columns(2)
                for j in range(2):
                    idx = i * 2 + j
                    if idx < len(selected_vars):
                        var = selected_vars[idx]
                        fig_box = px.box(data, y=var, title=f"Boxplot {var}", points="all", 
                                         color_discrete_sequence=['lightcoral'], template=pio.templates.default)
                        cols[j].plotly_chart(fig_box, use_container_width=True)
        
        with tab_map:
            st.subheader("Peta Sebaran Data")
            map_var_eda = st.radio("Pilih Variabel untuk Peta:", selected_vars, horizontal=True, key="map_asli_radio")
            
            # Gunakan variabel standar (yang sudah diimputasi jika ada NA)
            target_col_eda = map_var_eda
            
            # Tentukan label legend berdasarkan variabel yang dipilih
            legend_label = "Nilai"
            if map_var_eda == "KUSTA":
                legend_label = "Kasus Per 100.000 Penduduk"
            elif map_var_eda in ["PHBS", "MISKIN"]:
                legend_label = "Persentase (%)"
            elif map_var_eda == "KEPADATAN":
                legend_label = "Jiwa/km²"

            # Verifikasi kolom ada di data sebelum plotting
            if target_col_eda in data.columns:
                fig_map_asli = px.choropleth(
                    data, geojson=data.geometry, locations=data.index,
                    color=target_col_eda, hover_name="shapeName",
                    color_continuous_scale="Viridis", 
                    labels={target_col_eda: legend_label},
                    title=f"Peta Sebaran {map_var_eda}"
                )
                fig_map_asli.update_geos(fitbounds="locations", visible=False)
                st.plotly_chart(fig_map_asli, use_container_width=True)
            else:
                st.error(f"Kolom '{target_col_eda}' tidak ditemukan. Harap bersihkan cache Streamlit (Menu > Clear Cache) dan muat ulang halaman.")

        st.markdown("---")
        
        # 3. ANALISIS KORELASI (SCATTER PLOT)
        st.subheader("Analisis Korelasi")
        col_x, col_y = st.columns(2)
        with col_x:
            var_x = st.selectbox("Sumbu X:", vars_list, index=0, key="scatter_x")
        with col_y:
            var_y = st.selectbox("Sumbu Y:", vars_list, index=1, key="scatter_y")
            
        if var_x == var_y:
            st.warning("⚠️ Pilih variabel Sumbu X dan Sumbu Y yang berbeda.")
        else:
            fig_scatter = px.scatter(
                data, x=var_x, y=var_y, trendline="ols",
                hover_name="shapeName", title=f"Korelasi: {var_x} vs {var_y}"
            )
            st.plotly_chart(fig_scatter, use_container_width=True)

# ==============================================================================
# 2. SPATIAL AUTOCORRELATION
# ==============================================================================
elif menu == "🌍 Spatial Autocorrelation":
    st.header(f"🌍 Spatial Autocorrelation - {weight_type}")
    var_auto = st.selectbox("Pilih Variabel:", vars_list)
    y = data[var_auto].values
    np.random.seed(42)
    
    st.subheader(f"1. Global Moran's I")
    mi = Moran(y, w)
    col1, col2, col3 = st.columns(3)
    col1.metric("Moran's I", f"{mi.I:.4f}")
    col2.metric("Z-Score", f"{mi.z_norm:.4f}")
    col3.metric("P-Value", f"{mi.p_sim:.4f}")
    
    # Interpretasi Moran's I
    if mi.p_sim < 0.05:
        if mi.I > 0:
            st.success(f"**Interpretasi:** Terdapat autokorelasi spasial positif (Clustered/Berkelompok). Nilai indeks {mi.I:.4f} > 0 dan p-value < 0.05.")
        else:
            st.success(f"**Interpretasi:** Terdapat autokorelasi spasial negatif (Dispersed/Menyebar). Nilai indeks {mi.I:.4f} < 0 dan p-value < 0.05.")
    else:
        st.info(f"**Interpretasi:** Tidak terdapat autokorelasi spasial yang signifikan (Pola Acak/Random). P-value {mi.p_sim:.4f} >= 0.05.")
    
    st.markdown("---")
    st.subheader(f"2. Geary's C")
    gc = Geary(y, w)
    col1, col2, col3 = st.columns(3)
    col1.metric("Geary's C", f"{gc.C:.4f}")
    col2.metric("Z-Score", f"{gc.z_norm:.4f}")
    col3.metric("P-Value", f"{gc.p_sim:.4f}")

    # Interpretasi Geary's C
    if gc.p_sim < 0.05:
        if gc.C < 1:
            st.success(f"**Interpretasi:** Terdapat autokorelasi spasial positif (Clustered/Berkelompok). Nilai indeks {gc.C:.4f} < 1 dan p-value < 0.05.")
        else:
            st.success(f"**Interpretasi:** Terdapat autokorelasi spasial negatif (Dispersed/Menyebar). Nilai indeks {gc.C:.4f} > 1 dan p-value < 0.05.")
    else:
        st.info(f"**Interpretasi:** Tidak terdapat autokorelasi spasial yang signifikan (Pola Acak/Random). P-value {gc.p_sim:.4f} >= 0.05.")

    st.markdown("---")
    st.subheader("3. LISA")
    with joblib.parallel_backend('threading'):
        lisa = Moran_Local(y, w)
        p_two = lisa.p_sim
        
    sig = p_two < 0.05
    hotspot = sig & (lisa.q==1)
    coldspot = sig & (lisa.q==2)
    doughnut = sig & (lisa.q==3)
    diamond = sig & (lisa.q==4)
    
    labels = pd.Series(["Not Significant"] * len(data), index=data.index)
    labels[hotspot] = "High-High"
    labels[coldspot] = "Low-Low"
    labels[doughnut] = "Low-High"
    labels[diamond] = "High-Low"
    data['lisa_label'] = labels
    
    fig_lisa = px.choropleth(
        data, geojson=data.geometry, locations=data.index,
        color='lisa_label',
        color_discrete_map={
            "High-High": "#d7191c", "Low-Low": "#2c7bb6",
            "High-Low": "#fdae61", "Low-High": "#abd9e9",
            "Not Significant": "#eeeeee"
        },
        hover_name="shapeName", title=f"LISA Map: {var_auto}"
    )
    fig_lisa.update_geos(fitbounds="locations", visible=False)
    st.plotly_chart(fig_lisa, use_container_width=True)

    st.markdown("---")
    st.subheader("4. Getis-Ord Gi* (Hotspot Analysis)")
    
    with joblib.parallel_backend('threading'):
        g_local = G_Local(y, w, transform='R', star=True)
        
    gi_sig = g_local.p_sim < 0.05
    gi_labels = pd.Series(["Not Significant"] * len(data), index=data.index)
    gi_labels[(g_local.Zs > 1.96) & gi_sig] = "Hotspot"
    gi_labels[(g_local.Zs < -1.96) & gi_sig] = "Coldspot"
    
    data['gi_label'] = gi_labels
    data['gi_zscore'] = g_local.Zs
    data['gi_pvalue'] = g_local.p_sim

    only_sig = st.checkbox("Hanya tampilkan wilayah signifikan (p < 0.05)", value=False)
    plot_data = data.copy()
    if only_sig:
        plot_data.loc[plot_data['gi_label'] == "Not Significant", 'gi_label'] = None

    fig_gi = px.choropleth(
        plot_data, geojson=plot_data.geometry, locations=plot_data.index,
        color='gi_label',
        color_discrete_map={"Hotspot": "#d73027", "Coldspot": "#4575b4", "Not Significant": "#f0f0f0"},
        hover_name="shapeName",
        hover_data={'gi_zscore':':.2f', 'gi_pvalue':':.4f'},
        title=f"Getis-Ord Gi* Map: {var_auto}"
    )
    fig_gi.update_geos(fitbounds="locations", visible=False)
    st.plotly_chart(fig_gi, use_container_width=True)
    st.info("ℹ️ **Hotspot**: Wilayah nilai tinggi dikelilingi nilai tinggi.")

# ==============================================================================
# 3. SPATIAL INTERPOLATION
# ==============================================================================
elif menu == "🗺️ Spatial Interpolation":
    st.header("🗺️ Spatial Interpolation")
    tab_diag, tab_idw = st.tabs(["Diagnosis Spasial & Tren", "Inverse Distance Weighting (IDW)"])
    
    with tab_diag:
        st.subheader("Visualisasi Tren & Semivariogram Cloud")
        diag_var = st.selectbox("Pilih Variabel untuk Diagnosis:", vars_list, key="diag_var")
        col_trend1, col_trend2 = st.columns(2)
        with col_trend1:
            fig_bx = px.scatter(data, x="centroid_x", y=diag_var, trendline="ols", title=f"Tren Barat - Timur", labels={"centroid_x": "Longitude"})
            st.plotly_chart(fig_bx, use_container_width=True)
        with col_trend2:
            fig_uy = px.scatter(data, x="centroid_y", y=diag_var, trendline="ols", title=f"Tren Utara - Selatan", labels={"centroid_y": "Latitude"})
            st.plotly_chart(fig_uy, use_container_width=True)
        
        st.markdown("##### Semivariogram Cloud")
        coords = list(zip(data.centroid_x, data.centroid_y))
        vals = data[diag_var].values
        d_ij = distance_matrix(coords, coords)
        val_matrix = np.abs(np.subtract.outer(vals, vals))
        gamma_ij = 0.5 * (val_matrix ** 2)
        iu = np.triu_indices(len(vals), k=1)
        distances = d_ij[iu]
        semivars = gamma_ij[iu]
        df_semi = pd.DataFrame({"Jarak": distances, "Semivariance": semivars})
        fig_semi = px.scatter(df_semi, x="Jarak", y="Semivariance", opacity=0.6, title=f"Semivariogram Cloud: {diag_var}")
        st.plotly_chart(fig_semi, use_container_width=True)

    with tab_idw:
        st.subheader("Konfigurasi & Prediksi IDW")
        col_conf1, col_conf2 = st.columns([1, 2])
        with col_conf1:
            target_idw = st.selectbox("Variabel Target:", vars_list, key="target_idw")
            p_val = st.slider("Power (p):", 1.0, 10.0, 2.0, 0.5)
            k_val = st.slider("Jumlah Tetangga (k):", 3, len(data)-1, 5)
            rmse_val, mae_val, res_df = calculate_idw_loocv(data, target_idw, p=p_val, k=k_val)
            st.metric("RMSE", f"{rmse_val:.4f}")
            st.metric("MAE", f"{mae_val:.4f}")
        with col_conf2:
            viz_type = st.radio("Pilih Tampilan:", ["Peta Residual", "Peta Permukaan"], horizontal=True)
            if viz_type == "Peta Residual":
                fig_res = px.choropleth(res_df, geojson=res_df.geometry, locations=res_df.index, color='Residual', color_continuous_scale="RdBu", title="Peta Residual IDW")
                fig_res.update_geos(fitbounds="locations", visible=False)
                st.plotly_chart(fig_res, use_container_width=True)
            else:
                x_min, x_max = data.centroid_x.min(), data.centroid_x.max()
                y_min, y_max = data.centroid_y.min(), data.centroid_y.max()
                grid_x = np.linspace(x_min-0.1, x_max+0.1, 50)
                grid_y = np.linspace(y_min-0.1, y_max+0.1, 50)
                XX, YY = np.meshgrid(grid_x, grid_y)
                known_coords = np.column_stack((data.centroid_x, data.centroid_y))
                known_z = data[target_idw].values
                grid_flat = np.column_stack((XX.ravel(), YY.ravel()))
                dists_g = distance_matrix(grid_flat, known_coords)
                grid_z = []
                for i in range(len(grid_flat)):
                    idx_s = np.argsort(dists_g[i])[:k_val]
                    d_s = dists_g[i][idx_s]
                    d_s[d_s == 0] = 1e-10
                    w_s = 1 / (d_s ** p_val)
                    grid_z.append(np.sum(w_s * known_z[idx_s]) / np.sum(w_s))
                ZZ = np.array(grid_z).reshape(XX.shape)
                fig_surf = go.Figure(data=go.Contour(z=ZZ, x=grid_x, y=grid_y, colorscale='Viridis'))
                fig_surf.update_layout(title=f"Permukaan Interpolasi {target_idw}")
                st.plotly_chart(fig_surf, use_container_width=True)
        
        st.subheader("Kalkulator Prediksi Titik Baru")
        min_x, max_x = 105.5, 108.8
        min_y, max_y = -8.0, -6.0
        col_calc1, col_calc2 = st.columns(2)
        with col_calc1:
            input_lon = st.number_input("Longitude (BT):", value=float(data.centroid_x.mean()), format="%.6f", min_value=105.0, max_value=109.0)
            st.caption(f"Wilayah Jabar: {min_x}° s/d {max_x}°")
            if not (min_x <= input_lon <= max_x):
                st.warning("⚠️ Koordinat diluar zona analisis. Estimasi mungkin berupa ekstrapolasi.")
                        
        with col_calc2:
            input_lat = st.number_input("Latitude (LS):", value=float(data.centroid_y.mean()), format="%.6f", min_value=-9.0, max_value=-5.0)
            st.caption(f"Wilayah Jabar: {max_y}° s/d {min_y}°") # min_y lebih kecil (lebih selatan) dari max_y
            if not (min_y <= input_lat <= max_y):
                st.warning("⚠️ Koordinat diluar zona analisis. Estimasi mungkin berupa ekstrapolasi.")
                
        if st.button("Hitung Prediksi"):
            known_coords = np.column_stack((data.centroid_x, data.centroid_y))
            known_z = data[target_idw].values
            target_point = np.array([[input_lon, input_lat]])
                    
            dists_t = distance_matrix(target_point, known_coords)[0]
                    
            idx_k = np.argsort(dists_t)[:k_val]
            d_k = dists_t[idx_k]
            z_k = known_z[idx_k]
            d_k[d_k == 0] = 1e-10
                    
            weights_k = 1 / (d_k ** p_val)
            prediction = np.sum(weights_k * z_k) / np.sum(weights_k)
                    
            st.success(f"Prediksi Nilai {target_idw} di ({input_lon}, {input_lat}): **{prediction:.4f}**")

# ==============================================================================
# 4. SPATIAL ECONOMETRICS (REVISI ALUR & 7 MODEL)
# ==============================================================================
elif menu == "📐 Spatial Econometrics":
    st.header(f"📐 Spatial Econometrics - {weight_type}")
    col1, col2 = st.columns(2)
    with col1:
        y_var = st.selectbox("Variabel Dependen (Y):", vars_list, index=0)
    with col2:
        x_options = [v for v in vars_list if v != y_var]
        x_vars = st.multiselect("Variabel Independen (X):", x_options, default=x_options)

    if st.button("Run Model"):
        if not x_vars:
            st.error("Pilih minimal satu variabel X.")
        else:
            # Simpan state bahwa analisis sudah dijalankan
            st.session_state['run_ekmet'] = True
            st.session_state['y_var'] = y_var
            st.session_state['x_vars'] = x_vars
            
            with st.spinner("Mengestimasi berbagai model spasial..."):
                Y = data[y_var].astype("float64").to_numpy().reshape(-1, 1)
                X = data[x_vars].astype("float64").to_numpy()
                
                # --- A. FIT MODEL ---
                models = {}
                
                # 1. OLS (Model Klasik)
                ols = OLS(Y, X, w=w, spat_diag=True, moran=True, name_y=y_var, name_x=x_vars, name_ds="Jabar")
                models['OLS'] = ols
                
                # 2. SAR (Spatial Autoregressive) - ML
                sar = ML_Lag(Y, X, w=w, name_y=y_var, name_x=x_vars)
                models['SAR'] = sar

                # 3. SEM (Spatial Error Model) - ML
                sem = ML_Error(Y, X, w=w, name_y=y_var, name_x=x_vars)
                models['SEM'] = sem
                
                # --- PREPARE LAGGED X (WX) FOR DURBIN & SLX ---
                # Lag Spatial untuk setiap kolom X
                WX = np.hstack([lag_spatial(w, X[:, i].reshape(-1, 1)) for i in range(X.shape[1])])
                X_durbin = np.hstack((X, WX))
                name_wx = [f"W_{col}" for col in x_vars]
                name_x_durbin = x_vars + name_wx
                
                # 4. SLX (Spatial Lag of X) - Dijalankan via OLS dengan X + WX
                slx = OLS(Y, X_durbin, w=w, name_y=y_var, name_x=name_x_durbin, name_ds="Jabar (SLX)")
                models['SLX'] = slx
                
                # 5. SDM (Spatial Durbin Model) - SAR dengan X + WX
                sdm = ML_Lag(Y, X_durbin, w=w, name_y=y_var, name_x=name_x_durbin)
                models['SDM'] = sdm
                
                # 6. SDEM (Spatial Durbin Error Model) - SEM dengan X + WX
                sdem = ML_Error(Y, X_durbin, w=w, name_y=y_var, name_x=name_x_durbin)
                models['SDEM'] = sdem
                
                # Simpan models ke session state
                st.session_state['ekmet_models'] = models

    # --- TAMPILKAN HASIL JIKA SUDAH DI-RUN ---
    if st.session_state.get('run_ekmet'):
        models = st.session_state['ekmet_models']
        ols = models['OLS']
        
        # 1. HASIL LM TEST
        st.subheader("1. Uji Lagrange Multiplier (LM Test)")
        try:
            lm_df = pd.DataFrame({
                "Uji Diagnostik": ["LM Error (SEM)", "LM Lag (SAR)", "Robust LM Error", "Robust LM Lag"],
                "Statistik": [ols.lm_error[0], ols.lm_lag[0], ols.rlm_error[0], ols.rlm_lag[0]],
                "P-Value": [ols.lm_error[1], ols.lm_lag[1], ols.rlm_error[1], ols.rlm_lag[1]]
            })
            st.table(lm_df.style.format({"Statistik": "{:.4f}", "P-Value": "{:.4f}"}))
            
            # Rekomendasi Sederhana
            sar_sig = ols.lm_lag[1] < 0.05
            sem_sig = ols.lm_error[1] < 0.05
            if sar_sig and sem_sig:
                if ols.rlm_lag[0] > ols.rlm_error[0]: rek = "SAR / SDM"
                else: rek = "SEM / SDEM"
                info = "Keduanya signifikan, cek Robust LM terbesar."
            elif sar_sig: 
                rek = "SAR"
                info = "Hanya Lag yang signifikan."
            elif sem_sig: 
                rek = "SEM"
                info = "Hanya Error yang signifikan."
            else: 
                rek = "OLS / SLX"
                info = "Tidak ada autokorelasi spasial signifikan."
                
            st.info(f"💡 Rekomendasi Awal: **{rek}** ({info})")
        except:
            st.warning("LM Test tidak tersedia.")
            
        # 2. PERBANDINGAN MODEL
        st.markdown("---")
        st.subheader("2. Perbandingan Performa Model")
        
        comp_list = []
        for name, mod in models.items():
            comp_list.append({
                "Model": name,
                "AIC": mod.aic,
                "Log-Likelihood": mod.logll,
                "R-Squared": mod.r2 if hasattr(mod, 'r2') else (mod.pr2 if hasattr(mod, 'pr2') else np.nan)
            })
            
        df_comp = pd.DataFrame(comp_list).sort_values("AIC")
        st.dataframe(df_comp.style.format({"AIC": "{:.2f}", "Log-Likelihood": "{:.2f}", "R-Squared": "{:.4f}"}))
        
        # Logika Pencarian Nilai Ekstrem
        min_aic_row = df_comp.loc[df_comp['AIC'].idxmin()]
        max_ll_row = df_comp.loc[df_comp['Log-Likelihood'].idxmax()]
        max_r2_row = df_comp.loc[df_comp['R-Squared'].idxmax()]
        
        # Menampilkan 3 Keterangan
        st.success(f"📉 **AIC Terendah:** {min_aic_row['Model']} ({min_aic_row['AIC']:.2f})")
        st.info(f"📈 **Log-Likelihood Tertinggi:** {max_ll_row['Model']} ({max_ll_row['Log-Likelihood']:.2f})")
        st.success(f"📈 **R-Squared Tertinggi:** {max_r2_row['Model']} ({max_r2_row['R-Squared']:.4f})")
        
        # 3. MODEL SUMMARY INTERAKTIF
        st.markdown("---")
        st.subheader("3. Detail Ringkasan Model")
        
        selected_model_name = st.selectbox("Pilih Model untuk melihat detail:", df_comp['Model'].tolist())
        
        if selected_model_name:
            st.markdown(f"**Analisis Detail untuk Model: {selected_model_name}**")
            model_selected = models[selected_model_name]
            
            # --- TABS ORGANISASI BARU ---
            tab_summary, tab_assumptions, tab_vis = st.tabs([
                "📊 Ringkasan Estimasi", 
                "✅ Uji Asumsi Klasik", 
                "🗺️ Visualisasi & Prediksi"
            ])
            
            # ========================
            # TAB 1: RINGKASAN ESTIMASI
            # ========================
            with tab_summary:
                # 1. TABEL EFEK LANGSUNG
                st.markdown("##### Tabel Direct Effect (Koefisien Regresi)")
                
                betas = model_selected.betas
                std_err = model_selected.std_err
                params_list = []
                
                def get_interp(p): return "✅ Signifikan" if p < 0.05 else "❌ Tidak Sig."
                
                # Logic Model Extraction (Direct Effects)
                coeffs = betas.flatten()
                ses = np.array(std_err).flatten()
                var_names = model_selected.name_x 
                
                # Determine P-values based on model type
                if hasattr(model_selected, 't_stat'): # OLS
                    pvals = np.array(model_selected.t_stat)[:, 1]
                elif hasattr(model_selected, 'z_stat'): # ML Models
                    pvals = np.array(model_selected.z_stat)[:, 1]
                else: # Fallback
                    t_stats = coeffs / (ses + 1e-10)
                    pvals = 2 * (1 - norm.cdf(np.abs(t_stats)))

                for i, name in enumerate(var_names):
                    params_list.append({
                        "Variabel": name,
                        "Koefisien": coeffs[i],
                        "Std. Error": ses[i],
                        "P-Value": pvals[i],
                        "Interpretasi": get_interp(pvals[i])
                    })

                # Display Direct Effects Table
                df_direct = pd.DataFrame(params_list)
                st.dataframe(df_direct.style.format({"Koefisien": "{:.4f}", "Std. Error": "{:.4f}", "P-Value": "{:.4f}"}))
                
                # 2. TABEL EFEK TIDAK LANGSUNG (Spatial Params)
                st.markdown("##### Tabel Parameter Spasial (Rho / Lambda)")
                spatial_params = []
                
                # Check for Lambda (SEM/SDEM)
                if hasattr(model_selected, 'lam'):
                    lam = model_selected.lam
                    if hasattr(model_selected, 'lam_se'):
                         lam_se = model_selected.lam_se
                         lam_z = lam / lam_se
                         lam_pval = 2 * (1 - norm.cdf(np.abs(lam_z)))
                         spatial_params.append({
                            "Parameter": "Lambda (λ) - Spatial Error",
                            "Nilai": lam,
                            "Std. Error": lam_se,
                            "P-Value": lam_pval,
                            "Interpretasi": get_interp(lam_pval)
                        })
                    else:
                        spatial_params.append({
                            "Parameter": "Lambda (λ) - Spatial Error",
                            "Nilai": lam,
                            "Std. Error": "-",
                            "P-Value": "-",
                            "Interpretasi": "Lihat LM Test"
                        })

                # Check for Rho (SAR/SDM)
                if hasattr(model_selected, 'rho'):
                    rho = model_selected.rho
                    if hasattr(model_selected, 'rho_se'):
                         rho_se = model_selected.rho_se
                         rho_z = rho / rho_se
                         rho_pval = 2 * (1 - norm.cdf(np.abs(rho_z)))
                         spatial_params.append({
                            "Parameter": "Rho (ρ) - Spatial Lag",
                            "Nilai": rho,
                            "Std. Error": rho_se,
                            "P-Value": rho_pval,
                            "Interpretasi": get_interp(rho_pval)
                        })
                    else:
                        spatial_params.append({
                            "Parameter": "Rho (ρ) - Spatial Lag",
                            "Nilai": rho,
                            "Std. Error": "-",
                            "P-Value": "-",
                            "Interpretasi": "-"
                        })
                
                if spatial_params:
                    df_spatial = pd.DataFrame(spatial_params)
                    st.dataframe(df_spatial.style.format({"Nilai": "{:.4f}"}))
                else:
                    st.info("Tidak ada parameter spasial (Rho/Lambda) untuk model ini (Model OLS/SLX).")

            # ========================
            # TAB 2: UJI ASUMSI KLASIK
            # ========================
            with tab_assumptions:
                st.markdown("##### 🔍 Diagnostik & Uji Asumsi Model")
                
                # Get Residuals
                if hasattr(model_selected, 'u'):
                    resid = model_selected.u.flatten()
                    
                    # 1. Multikolinearitas (VIF)
                    st.markdown("**1. Multikolinearitas (VIF)**")
                    st.caption("Memeriksa korelasi antar variabel independen. VIF < 10 aman.")
                    
                    # Prepare Data for VIF (Calculate from correlations)
                    x_vars = st.session_state['x_vars']
                    if len(x_vars) > 1:
                        df_x = data[x_vars]
                        corr_matrix = df_x.corr()
                        try:
                            inv_corr = np.linalg.inv(corr_matrix.values)
                            vif_vals = np.diag(inv_corr)
                            df_vif = pd.DataFrame({"Variabel": x_vars, "VIF": vif_vals})
                            st.dataframe(df_vif.style.format({"VIF": "{:.2f}"}))
                            if any(vif_vals > 10):
                                st.warning("⚠️ Terdapat Variabel dengan VIF > 10. Indikasi Multikolinearitas Tinggi.")
                            else:
                                st.success("✅ Tidak terdeteksi Multikolinearitas serius (Semua VIF < 10).")
                        except:
                            st.error("Gagal menghitung VIF (Singular Matrix).")
                    else:
                        st.info("Hanya 1 variabel independen, VIF tidak diperlukan.")
                        
                    st.markdown("---")
                    
                    # 2. Normalitas Residual (Shapiro-Wilk)
                    st.markdown("**2. Normalitas Residual**")
                    col_norm1, col_norm2 = st.columns([1, 2])
                    
                    with col_norm1:
                        stat_sw, p_sw = shapiro(resid)
                        st.metric("Shapiro-Wilk Stat", f"{stat_sw:.4f}")
                        st.metric("P-Value", f"{p_sw:.4f}")
                        
                        if p_sw > 0.05:
                            st.success("✅ Residual berdistribusi Normal (P > 0.05).")
                        else:
                            st.warning("⚠️ Residual TIDAK berdistribusi Normal (P < 0.05).")
                    
                    with col_norm2:
                        fig_hist_res = px.histogram(resid, title="Histogram Residual", template="plotly_white")
                        st.plotly_chart(fig_hist_res, use_container_width=True)
                        
                    st.markdown("---")
                    
                    # 3. Homoskedastisitas (Breusch-Pagan)
                    st.markdown("**3. Homoskedastisitas (Uji Breusch-Pagan)**")
                    
                    # Manual BP Test calculation for consistency across ML/OLS models
                    # LM = n * R2 from regression of u^2 on X
                    try:
                        u2 = resid ** 2
                        X_mat = data[x_vars].values
                        # Add constant for auxiliary regression
                        X_aux = np.column_stack([np.ones(len(u2)), X_mat])
                        
                        # OLS u2 on X
                        coeffs_bp = np.linalg.inv(X_aux.T @ X_aux) @ X_aux.T @ u2
                        u2_pred = X_aux @ coeffs_bp
                        
                        # R2 Calculation
                        ss_tot = np.sum((u2 - np.mean(u2))**2)
                        ss_res = np.sum((u2 - u2_pred)**2)
                        r2_bp = 1 - (ss_res / ss_tot)
                        
                        lm_bp = len(u2) * r2_bp
                        df_bp = len(x_vars)
                        p_bp = 1 - chi2.cdf(lm_bp, df_bp)
                        
                        col_bp1, col_bp2 = st.columns(2)
                        col_bp1.metric("BP Statistic", f"{lm_bp:.4f}")
                        col_bp2.metric("P-Value", f"{p_bp:.4f}")
                        
                        if p_bp > 0.05:
                            st.success("✅ Homoskedastisitas Terpenuhi (Varian residual konstan).")
                        else:
                            st.warning("⚠️ Terdeteksi Heteroskedastisitas (P < 0.05).")
                            
                    except Exception as e:
                        st.error(f"Gagal melakukan uji BP: {e}")
                        
                    st.markdown("---")
                    
                    # 4. Autokorelasi Spasial pada Residual (Moran's I)
                    st.markdown("**4. Autokorelasi Spasial Residual (Global Moran's I)**")
                    mi_res = Moran(resid, w)
                    
                    col_mi1, col_mi2, col_mi3 = st.columns(3)
                    col_mi1.metric("Moran's I Residual", f"{mi_res.I:.4f}")
                    col_mi2.metric("P-Value", f"{mi_res.p_sim:.4f}")
                    col_mi3.metric("Z-Score", f"{mi_res.z_norm:.4f}")
                    
                    if mi_res.p_sim > 0.05:
                        st.success("✅ Tidak ada autokorelasi spasial tersisa pada residual (Model sudah baik menangkap efek spasial).")
                    else:
                        st.warning("⚠️ Masih terdapat autokorelasi spasial pada residual. Pertimbangkan model spasial lain (misal: SDM/SDEM).")

                else:
                    st.warning("Residual tidak tersedia untuk model ini, uji asumsi tidak dapat dilakukan.")

            # ========================
            # TAB 3: VISUALISASI & PREDIKSI
            # ========================
            with tab_vis:
                # --- PLOT RESIDUAL ---
                st.markdown("##### Visualisasi Residual Model")
                
                if hasattr(model_selected, 'u'):
                    data['model_resid'] = model_selected.u
                    
                    col_res1, col_res2 = st.columns([1.5, 1])
                    
                    with col_res1:
                        fig_res_map = px.choropleth(
                            data, geojson=data.geometry, locations=data.index,
                            color='model_resid',
                            color_continuous_scale="RdBu",
                            title=f"Peta Residual: {selected_model_name}",
                            hover_name="shapeName"
                        )
                        fig_res_map.update_geos(fitbounds="locations", visible=False)
                        st.plotly_chart(fig_res_map, use_container_width=True)
                    
                    with col_res2:
                        # Scatter Fitted vs Residual (jika ada predicted value / predy)
                        if hasattr(model_selected, 'predy'):
                            fitted_vals = model_selected.predy.flatten()
                            resid_vals = model_selected.u.flatten()
                            df_scatter_res = pd.DataFrame({'Fitted': fitted_vals, 'Residuals': resid_vals})
                            
                            fig_scatter_res = px.scatter(
                                df_scatter_res, x='Fitted', y='Residuals',
                                title=f"Fitted vs Residuals",
                                trendline="ols"
                            )
                            fig_scatter_res.add_hline(y=0, line_dash="dash", line_color="red")
                            st.plotly_chart(fig_scatter_res, use_container_width=True)
                
                # --- SIMULASI PREDIKSI ---
                st.markdown("---")
                st.markdown("##### Simulasi Prediksi Y (Hypothetical)")
                st.info("Prediksi Y mempertimbangkan **Nilai X**, **Lokasi Spasial**, dan **Efek Tetangga** (Nearest Neighbors = 5).")

                # A. Input Variabel X
                st.markdown("**1. Input Nilai Variabel Independen (X)**")
                x_vars = st.session_state['x_vars']
                pred_cols = st.columns(len(x_vars))
                input_vals = {}
                for i, col in enumerate(pred_cols):
                    with col:
                        default_val = float(data[x_vars[i]].mean())
                        input_vals[x_vars[i]] = st.number_input(f"{x_vars[i]}", value=default_val, format="%.2f")

                # B. Input Lokasi (Kontekstualisasi Spasial)
                st.markdown("**2. Input Lokasi Target (Longitude & Latitude)**")
                col_loc1, col_loc2 = st.columns(2)
                with col_loc1:
                    input_sim_lon = st.number_input("Longitude Target:", value=107.6, format="%.4f")
                with col_loc2:
                    input_sim_lat = st.number_input("Latitude Target:", value=-6.9, format="%.4f")

                if st.button("Hitung Prediksi Spasial"):
                    # 1. Cari 5 Tetangga Terdekat dari data asli
                    known_coords = np.column_stack((data.centroid_x, data.centroid_y))
                    d_sim = np.sqrt((known_coords[:,0] - input_sim_lon)**2 + (known_coords[:,1] - input_sim_lat)**2)
                    
                    k_neighbors = 5
                    idx_sim = np.argsort(d_sim)[:k_neighbors] # Indeks 5 tetangga terdekat
                    d_k = d_sim[idx_sim]
                    d_k[d_k == 0] = 1e-10
                    
                    # Hitung Bobot Tetangga (Inverse Distance)
                    w_sim = (1 / d_k**2)
                    w_sim = w_sim / np.sum(w_sim) # Normalize weights to sum to 1
                    
                    # 2. Hitung Nilai Lag Tetangga (Spatially weighted values of neighbors)
                    # Ambil data tetangga
                    neighbors_df = data.iloc[idx_sim]
                    
                    # Lag Y (Spatially lagged dependent variable) -> Untuk model SAR, SDM
                    y_col = st.session_state['y_var']
                    neighbors_y = neighbors_df[y_col].values
                    lag_y_sim = np.sum(w_sim * neighbors_y)
                    
                    # Lag X (Spatially lagged independent variables) -> Untuk model SLX, SDM, SDEM
                    lag_x_sim = {}
                    for var in x_vars:
                        neighbors_x = neighbors_df[var].values
                        lag_x_sim[var] = np.sum(w_sim * neighbors_x)

                    # 3. Hitung Prediksi Berdasarkan Rumus Model
                    # Y_pred = Intercept + X*Beta + WX*Theta + Rho*Wy
                    
                    betas = model_selected.betas.flatten()
                    var_names = model_selected.name_x
                    
                    # Inisialisasi Prediksi
                    y_pred = 0
                    components_text = []
                    
                    for idx, var_name in enumerate(var_names):
                        coeff = betas[idx]
                        
                        if var_name == "CONSTANT" or var_name == "intercept":
                            y_pred += coeff
                            
                        elif var_name in input_vals:
                            val = input_vals[var_name]
                            term = coeff * val
                            y_pred += term
                            
                        elif var_name.startswith("W_"):
                            original_var = var_name[2:] 
                            if original_var in lag_x_sim:
                                val = lag_x_sim[original_var]
                                term = coeff * val
                                y_pred += term
                                
                    # Tambahkan Efek Lag Y (Rho) jika ada
                    if hasattr(model_selected, 'rho'):
                        rho = model_selected.rho
                        term = rho * lag_y_sim
                        y_pred += term
                        
                    st.success(f"Prediksi Nilai **{y_col}** di lokasi tersebut: **{y_pred:.4f}**")
                    st.info(f"Model menggunakan nilai {y_col} rata-rata tertimbang dari 5 kabupaten tetangga terdekat (Lag Y = {lag_y_sim:.2f}) sebagai referensi efek spasial.")

            # --- 5. DOWNLOAD REPORT ---
            st.markdown("---")
            st.subheader("Ekspor Laporan")
            
            if st.button("Siapkan File PDF"):
                summary_txt = str(model_selected.summary)
                pdf_bytes = create_pdf_report(selected_model_name, summary_txt)
                
                st.download_button(
                    label="⬇️ Unduh Laporan PDF",
                    data=pdf_bytes,
                    file_name=f"Laporan_Spasial_{selected_model_name}.pdf",
                    mime="application/pdf"
                )

# ==============================================================================
# 5. GWR
# ==============================================================================
elif menu == "🧭 GWR":
    st.header("🧭 Geographically Weighted Regression")
    
    col_gwr1, col_gwr2 = st.columns(2)
    with col_gwr1:
        # Pilihan Y Bebas
        y_var_gwr = st.selectbox("Variabel Dependen (Y):", vars_list, index=0)
    with col_gwr2:
        # Pilihan X Bebas (termasuk KUSTA jika Y bukan KUSTA)
        x_options_gwr = [v for v in vars_list if v != y_var_gwr]
        x_vars_gwr = st.multiselect("Variabel Independen (X):", x_options_gwr, default=[x_options_gwr[0]])
    if st.button("Run Model"):
        if not x_vars_gwr:
             st.error("Pilih minimal satu variabel X.")
        else:
            with st.spinner("Mencari bandwidth optimal dan menghitung GWR..."):
                coords = list(zip(data.centroid_x, data.centroid_y))
                y = data[y_var_gwr].values.reshape((-1, 1))
                X = data[x_vars_gwr].values
                
                # 1. Bandwidth & GWR
                sel = Sel_BW(coords, y, X)
                bw = sel.search(bw_min=2)
                model = GWR(coords, y, X, bw)
                res = model.fit()
                
                # 2. OLS Global untuk perbandingan
                ols_global = OLS(y, X, name_y=y_var_gwr, name_x=x_vars_gwr)

                # Simpan ke Session State
                st.session_state['gwr_run'] = True
                st.session_state['gwr_res'] = res
                st.session_state['gwr_bw'] = bw
                st.session_state['gwr_ols'] = ols_global
                st.session_state['gwr_x_vars'] = x_vars_gwr

    # LOGIKA TAMPILAN
    if st.session_state.get('gwr_run'):
        res = st.session_state['gwr_res']
        ols = st.session_state['gwr_ols']
        x_vars = st.session_state['gwr_x_vars']
        bw = st.session_state['gwr_bw']
        
        # 1. METRIK PERBANDINGAN
        st.markdown("### 1. Perbandingan Performa Model (GWR vs OLS)")
        col_m1, col_m2, col_m3 = st.columns(3)
        col_m1.metric("Optimal Bandwidth", int(bw))
        delta_aicc = res.aicc - ols.aic 
        col_m2.metric("AICc (Lower is better)", f"{res.aicc:.2f}", delta=f"{delta_aicc:.2f}", delta_color="inverse")
        delta_r2 = res.R2 - ols.r2
        col_m3.metric("R-Square (Higher is better)", f"{res.R2:.4f}", delta=f"{delta_r2:.4f}", delta_color="normal")
        
        # 2. VISUALISASI
        st.markdown("### 2. Visualisasi Lokal")
        
        viz_type = st.radio("Pilih Jenis Peta:", ["Koefisien Lokal", "Signifikansi Lokal (t-test)", "Local R-Square"], horizontal=True)
        
        if viz_type == "Local R-Square":
            data['local_R2'] = res.localR2
            fig = px.choropleth(
                data, geojson=data.geometry, locations=data.index,
                color='local_R2', color_continuous_scale="Viridis",
                title="Sebaran Local R-Square (Akurasi Lokal)",
                hover_name="shapeName"
            )
            fig.update_geos(fitbounds="locations", visible=False)
            st.plotly_chart(fig, use_container_width=True)
            
        else:
            param_options = ["Intercept"] + x_vars
            selected_param = st.selectbox(f"Pilih Variabel untuk {viz_type}:", param_options)
            idx = 0 if selected_param == "Intercept" else x_vars.index(selected_param) + 1
            
            if viz_type == "Koefisien Lokal":
                data['gwr_coef'] = res.params[:, idx]
                fig = px.choropleth(
                    data, geojson=data.geometry, locations=data.index,
                    color='gwr_coef', color_continuous_scale="RdBu_r",
                    title=f"Koefisien Lokal: {selected_param}",
                    hover_name="shapeName"
                )
                fig.update_geos(fitbounds="locations", visible=False)
                st.plotly_chart(fig, use_container_width=True)
                
            elif viz_type == "Signifikansi Lokal (t-test)":
                t_vals = res.params[:, idx] / res.bse[:, idx]
                sig_labels = []
                for t in t_vals:
                    if t > 1.96: sig_labels.append("Signifikan Positif")
                    elif t < -1.96: sig_labels.append("Signifikan Negatif")
                    else: sig_labels.append("Tidak Signifikan")
                
                data['gwr_sig'] = sig_labels
                color_map = {"Signifikan Positif": "red", "Signifikan Negatif": "blue", "Tidak Signifikan": "#eeeeee"}
                
                fig = px.choropleth(
                    data, geojson=data.geometry, locations=data.index,
                    color='gwr_sig', color_discrete_map=color_map,
                    title=f"Signifikansi Lokal (t-stat > 1.96): {selected_param}",
                    hover_name="shapeName"
                )
                fig.update_geos(fitbounds="locations", visible=False)
                st.plotly_chart(fig, use_container_width=True)