Ujian Akhir Semester - Epidemiologi
Analisis Epidemiologi Kasus Kronis Filariasis di Jawa Barat Tahun 2024
Ditujukan Untuk Memenuhi Tugas UAS Mata Kuliah Epidemiologi
Dosen Pengampu:
Dr. I Gede Nyoman Mindra Jaya, M.Si
Disusun oleh:
Nadwah Khairunnisa (140610230034)
PROGRAM STUDI S1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

Abstrak

Filariasis merupakan penyakit menular yang dapat menyebabkan kecacatan kronis dan berdampak pada kualitas hidup serta produktivitas masyarakat. Studi ini melakukan mini project analisis epidemiologi filariasis kronis di Provinsi Jawa Barat pada unit analisis kabupaten/kota (\(n = 27\)). Analisis mencakup: (i) ukuran frekuensi kejadian dan keparahan (count, rate per 10.000 penduduk, mortalitas, CFR), (ii) ukuran asosiasi sederhana berbasis tabel \(2\times2\) (RR, OR, AR/RD, AF) untuk faktor sosial-lingkungan terpilih, (iii) Standardized Incidence Ratio (SIR) untuk membandingkan observed vs expected, (iv) pemetaan tematik dan zona risiko terpadu, (v) pemodelan count (Poisson vs Negative Binomial) serta (vi) uji dan pemodelan spasial (Moran, Geary, LISA, Getis–Ord \(Gi^*\), OLS dan model spasial—SLX/SAR/SEM/SDM/SDEM/SAC/GNS).

Hasil ringkas tingkat provinsi menunjukkan total populasi 50.345.190, total kasus kronis 262, dan total kematian 65 dengan rate 0,052 per 10.000 serta CFR 24,8%. Pada analisis spasial global untuk rate filariasis, Moran’s \(I\) menunjukkan kecenderungan pengelompokan (\(p \approx 0{,}052\)) dan Geary’s \(C\) signifikan (\(p \approx 0{,}038\)). Perbandingan model spasial berbasis AIC mengindikasikan model terbaik adalah SEM. Temuan ini menegaskan pentingnya pendekatan epidemiologi-spasial untuk memahami heterogenitas risiko antarwilayah dan mendukung penentuan prioritas intervensi.


Bab 1. Pendahuluan

1.1 Latar Belakang

Filariasis limfatik (lymphatic filariasis), sering dikenal sebagai elephantiasis, adalah penyakit parasit yang ditularkan melalui gigitan nyamuk yang membawa larva cacing filaria. Infeksi umumnya terjadi sejak masa kanak-kanak dan dapat menimbulkan kerusakan sistem limfatik yang bersifat “tersembunyi” dalam jangka panjang; manifestasi yang terlihat seperti limfedema/elefantiasis atau hidrokel pada laki-laki sering muncul kemudian pada usia dewasa. Kondisi-kondisi tersebut berpotensi menyebabkan pembengkakan menetap, keterbatasan aktivitas, nyeri berulang, infeksi sekunder, dan disabilitas kronis. Dampak lanjutannya tidak hanya medis, tetapi juga sosial dan ekonomi: stigma, penurunan partisipasi pendidikan/kerja, hingga beban biaya rumah tangga.

Dari perspektif kesehatan masyarakat, beban filariasis kronis memiliki dua karakter penting. Pertama, filariasis merupakan NTD yang kuat keterkaitannya dengan kemiskinan, lingkungan permukiman, dan ketersediaan infrastruktur sanitasi. Wilayah dengan kualitas perumahan rendah, kepadatan tertentu, drainase buruk, serta akses air bersih dan sanitasi yang tidak memadai cenderung menyediakan kondisi yang lebih mendukung habitat vektor serta meningkatkan peluang paparan gigitan nyamuk. Kedua, kasus kronis adalah indikator “beban morbiditas” yang persisten: kasus kronis dapat tetap ada meskipun transmisi menurun, sehingga membutuhkan strategi pelayanan kesehatan yang berkelanjutan melalui morbidity management and disability prevention (MMDP) untuk mengurangi disabilitas, mencegah komplikasi, dan meningkatkan kualitas hidup.

Di Jawa Barat, heterogenitas antar kabupaten/kota berpotensi menciptakan perbedaan risiko yang nyata. Perbedaan tersebut dapat berkaitan dengan indikator determinan kesehatan seperti:

  • Sanitasi layak (X1) yang mencerminkan kualitas lingkungan dan potensi pengelolaan tempat perindukan vektor secara tidak langsung,
  • Rumah layak huni (X2) yang mencerminkan kualitas hunian dan perlindungan terhadap paparan vektor (misalnya ventilasi, kerapatan bangunan, kondisi dinding/atap),
  • Rasio puskesmas terhadap penduduk (X3) sebagai proksi akses layanan kesehatan yang berpengaruh pada deteksi, pencatatan, tata laksana, dan pelaporan kasus kronis,
  • Persentase penduduk miskin (X4) sebagai determinan sosial yang berkaitan dengan kerentanan, keterbatasan pencegahan, dan keterbatasan akses layanan.

Walaupun program eliminasi filariasis menekankan penghentian transmisi, dalam praktik kebijakan kesehatan daerah masih dibutuhkan informasi yang lebih operasional: wilayah mana yang relatif lebih berisiko dan determinannya apa. Analisis epidemiologi berbasis wilayah (kab/kota) memberikan manfaat karena:

  1. Memetakan beban (berapa banyak kasus dan laju per penduduk),
  2. Mengidentifikasi determinan wilayah yang mungkin dapat menjadi prioritas intervensi lintas sektor (sanitasi, perumahan, kemiskinan),
  3. Mengarahkan layanan kesehatan untuk penemuan kasus, perawatan, dan manajemen disabilitas pada area yang lebih membutuhkan.

Dalam kerangka pembangunan berkelanjutan, filariasis termasuk NTD yang secara eksplisit terkait dengan SDG 3 (Good Health and Well-being), khususnya target 3.3 yang menekankan pengendalian penyakit menular termasuk NTDs. (UN SDG 3.3) Selain itu, determinan yang diuji (sanitasi dan perumahan) memiliki irisan kebijakan dengan SDG 6 (air bersih dan sanitasi) dan SDG 1 (pengentasan kemiskinan) sebagai akar struktural kerentanan. Dari sudut pandang pembangunan, NTDs sering disebut sebagai “penyakit kemiskinan” karena berhubungan dengan akses layanan, lingkungan, dan kondisi sosial-ekonomi yang tertinggal. (World Bank Blogs)

Berdasarkan kebutuhan praktis dan perencanaan program, terdapat setidaknya tiga celah yang umum muncul pada analisis penyakit menular berbasis wilayah:

  1. Dominasi pelaporan deskriptif: sering berhenti pada “jumlah kasus per wilayah” tanpa mengaitkan determinan sosial–lingkungan secara sistematis.
  2. Kurangnya standardisasi risiko: perbandingan antarwilayah sering bias oleh ukuran populasi; SIR membantu membandingkan kejadian teramati terhadap ekspektasi berdasarkan tingkat referensi.
  3. Keterbatasan pemeriksaan pola spasial: wilayah bertetangga berpotensi berbagi karakter lingkungan dan mobilitas; menguji apakah laju kasus membentuk klaster dapat memperkuat prioritisasi intervensi berbasis area.

Dengan demikian, penelitian ini berupaya menyusun laporan epidemiologi yang lebih komprehensif (teori → metodologi → hasil pada Bab 4), dengan tetap menjaga interpretasi sesuai desain area-level (ekologis) dan menghindari klaim kausal individual.

1.2 Identifikasi Masalah

  1. Apakah Rate Kasus Kronis Filariasis menunjukkan autokorelasi spasial pada tingkat kabupaten/kota di Jawa Barat?
  2. Wilayah mana yang membentuk cluster (High–High atau Low–Low) serta outlier (High–Low atau Low–High) berdasarkan analisis lokal?
  3. Wilayah mana yang teridentifikasi sebagai hotspot dan coldspot berdasarkan Getis–Ord Gi*?
  4. Model spatial econometrics mana yang paling sesuai untuk memodelkan hubungan Rate Kasus Kronis Filariasis dengan Persentase Sanitasi Layak, Persentase Rumah Layak Huni, Rasio Puskesmas terhadap Penduduk, dan Persentase Penduduk Miskin, setelah mempertimbangkan struktur spasial?
  5. Seberapa baik pendekatan interpolasi (IDW vs Ordinary Kriging) dalam membentuk permukaan prediksi Rate Kasus Kronis Filariasis?

1.3 Tujuan Penelitian

  1. Memetakan distribusi Jumlah dan Rate Kasus Kronis Filariasis, dan variabel penjelas yaitu Persentase Sanitasi Layak, Persentase Rumah Layak Huni, Rasio Puskesmas terhadap Penduduk, dan Persentase Penduduk Miskinpada kabupaten/kota di Jawa Barat.
  2. Menguji autokorelasi spasial global menggunakan Moran’s I dan Geary’s C.
  3. Mengidentifikasi cluster dan outlier menggunakan LISA (Local Moran).
  4. Mengidentifikasi hotspot/coldspot menggunakan Getis–Ord Gi*.
  5. Mengestimasi dan membandingkan beberapa model spatial econometrics (SLX, SAR, SEM, SDM, SDEM, SAC, GNS) serta menentukan model terbaik berdasarkan AIC/BIC/LogLik.
  6. Memetakan nilai prediksi (fitted) dan residual model terbaik.
  7. Melakukan interpolasi spasial (IDW dan Ordinary Kriging) serta mengevaluasi akurasi melalui LOO-CV (RMSE).

1.4 Batasan Penelitian

  • Unit analisis dibatasi pada 27 kabupaten/kota di Provinsi Jawa Barat, sesuai cakupan shapefile RBI 50K 2023.
  • Variabel respon utama untuk analisis spasial dan pemodelan adalah Rate Kasus Kronis Filariasis.
  • Variabel Jumlah Kasus Kronis Filariasis digunakan untuk distribusi peta spasial dan deskriptif awal.
  • Bobot spasial utama menggunakan Queen contiguity dengan analisis sensitivitas pada Rook dan kNN-4.
  • Interpolasi dilakukan pada titik representatif wilayah (centroid di dalam poligon) dan dibatasi pada batas administrasi gabungan Jawa Barat.

Bab 2. TINJAUAN PUSTAKA

2.1 Teori Dasar Spatial Dependence

Spatial dependence menyatakan bahwa observasi antar lokasi tidak independen; nilai variabel pada lokasi \(i\) dapat berkaitan dengan nilai pada lokasi tetangga \(j\) melalui struktur kedekatan yang dirangkum oleh matriks bobot spasial \(\mathbf{W}\). Ketergantungan ini dapat muncul karena proses difusi, kemiripan lingkungan, atau keterhubungan sosial-ekonomi antarwilayah. Dalam analisis penyakit menular, fenomena seperti mobilitas penduduk, kesamaan kondisi sanitasi, dan kesenjangan akses layanan kesehatan dapat menimbulkan pola pengelompokan spasial (Cliff & Ord, 1981; Anselin, 1988).

2.1 Autokorelasi Spasial (Global & Lokal)

Autokorelasi spasial mengukur kemiripan nilai antar unit spasial berdasarkan bobot \(\mathbf{W}\).

2.1.1 Moran’s I (Global)

Salah satu bentuk umum Moran’s I global: \[ I = \frac{n}{S_0}\,\frac{\sum_i\sum_j w_{ij}(y_i-\bar y)(y_j-\bar y)}{\sum_i (y_i-\bar y)^2}, \quad S_0 = \sum_i\sum_j w_{ij}, \] dengan \(n\) jumlah wilayah dan \(w_{ij}\) elemen matriks bobot spasial. Nilai \(I>0\) mengindikasikan cluster (nilai mirip berdekatan), sedangkan \(I<0\) mengindikasikan pola menyebar (Anselin, 1988).

2.1.2 Geary’s C (Global)

Bentuk umum Geary’s C: \[ C = \frac{(n-1)}{2S_0}\,\frac{\sum_i\sum_j w_{ij}(y_i-y_j)^2}{\sum_i (y_i-\bar y)^2}. \]

Nilai \(C<1\) menunjukkan autokorelasi positif, \(C>1\) autokorelasi negatif. Geary’s C relatif lebih sensitif pada variasi lokal dibanding Moran’s I (Cliff & Ord, 1981).

2.1.3 LISA (Local Moran)

Local Moran untuk wilayah \(i\): \[ I_i = z_i \sum_j w_{ij} z_j, \quad z_i = \frac{y_i-\bar y}{s_y}. \]

Nilai \(I_i\) dan signifikansinya digunakan untuk mengklasifikasikan pola lokal: High–High, Low–Low, High–Low, Low–High (Anselin, 1995).

2.1.4 Getis–Ord Gi*

Getis–Ord \(G_i^*\) (versi z-score) digunakan untuk mendeteksi hotspot/coldspot berdasarkan konsentrasi nilai tinggi/rendah di sekitar \(i\) (Getis & Ord, 1992; Ord & Getis, 1995).

2.2 Matriks Bobot Spasial (\(\mathbf{W}\))

Matriks bobot spasial \(\mathbf{W}\) merepresentasikan struktur ketetanggaan antarwilayah. Pada data poligon kabupaten/kota, bobot sering dibangun dari kontiguitas:

  • Queen contiguity: bertetangga jika berbagi sisi atau titik sudut.
  • Rook contiguity: bertetangga jika berbagi sisi.
  • k-nearest neighbors (kNN): setiap wilayah dihubungkan dengan \(k\) tetangga terdekat berdasarkan jarak (Bivand, Pebesma, & Gómez-Rubio, 2013).

Bobot kemudian dinormalisasi (mis. row-standardized, style = “W”) agar jumlah bobot per baris menjadi 1, sehingga interpretasi spatial lag lebih stabil.

2.3 Model Spatial Econometrics (Definisi + Model Matematis)

Misalkan \(\mathbf{y}\) adalah vektor respon \((n\times 1)\), \(\mathbf{X}\) matriks kovariat \((n\times k)\), \(\beta\) koefisien, dan \(\varepsilon\) error i.i.d. dengan varian \(\sigma^2\). Matriks bobot spasial adalah \(\mathbf{W}\).

Notasi ringkas: spatial lag pada variabel \(\mathbf{y}\) ditulis \(\mathbf{Wy}\), dan spatial lag pada kovariat ditulis \(\mathbf{WX}\).

2.3.1 Ordinary Least Square (OLS)

\[ \mathbf{y} = \mathbf{X}\beta + \varepsilon. \]

OLS sebagai baseline mengasumsikan error independen. Jika terdapat autokorelasi spasial pada \(\varepsilon\) atau struktur dependensi pada \(\mathbf{y}\), OLS menjadi tidak efisien dan inferensi dapat bias (Anselin, 1988; LeSage & Pace, 2009).

2.3.2 SLX (Spatial Lag of X)

\[ \mathbf{y} = \mathbf{X}\beta + \mathbf{WX}\theta + \varepsilon. \]

Model SLX menangkap spillover kovariat dari wilayah tetangga melalui \(\mathbf{WX}\) (LeSage & Pace, 2009).

2.3.3 SAR / SLM (Spatial Autoregressive / Spatial Lag Model)

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\beta + \varepsilon, \quad |\rho|<1. \]

Parameter \(\rho\) mengukur pengaruh rata-rata nilai respon wilayah tetangga terhadap wilayah \(i\) (Anselin, 1988; LeSage & Pace, 2009).

2.3.4 SEM (Spatial Error Model)

\[ \mathbf{y} = \mathbf{X}\beta + \mathbf{u},\quad \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \varepsilon, \quad |\lambda|<1. \]

Dependensi spasial masuk melalui error \(\mathbf{u}\), cocok bila ada faktor tak-teramati yang berpola spasial (Anselin, 1988; Elhorst, 2014).

2.3.5 SDM (Spatial Durbin Model)

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\beta + \mathbf{WX}\theta + \varepsilon. \]

SDM menggabungkan spatial lag respon dan spatial lag kovariat sehingga fleksibel dalam menangkap direct dan indirect effects (LeSage & Pace, 2009; Elhorst, 2014).

2.3.6 SDEM (Spatial Durbin Error Model)

\[ \mathbf{y} = \mathbf{X}\beta + \mathbf{WX}\theta + \mathbf{u},\quad \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \varepsilon. \]

SDEM adalah kombinasi SEM + lag kovariat (Elhorst, 2014).

2.3.7 SAC / SARAR (Spatial Autoregressive Combined)

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\beta + \mathbf{u},\quad \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \varepsilon. \]

SAC mengizinkan dependensi pada \(\mathbf{y}\) sekaligus dependensi pada error (Anselin, 1988; LeSage & Pace, 2009).

2.3.8 GNS (General Nesting Spatial Model)

Salah satu bentuk umum GNS: \[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\beta + \mathbf{WX}\theta + \mathbf{u},\quad \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \varepsilon. \]

GNS mencakup struktur SDM + error spasial. Secara praktis, identifikasi model dapat menantang dan sering memerlukan kehati-hatian dalam interpretasi (LeSage & Pace, 2009; Elhorst, 2014).

2.4 Interpolasi Spasial (IDW & Ordinary Kriging)

2.4.1 IDW

Prediksi IDW pada lokasi \(s_0\): \[ \hat z(s_0) = \frac{\sum_{i=1}^{n} w_i(s_0)\,z(s_i)}{\sum_{i=1}^{n} w_i(s_0)},\quad w_i(s_0)=\frac{1}{d(s_0,s_i)^p}, \] dengan \(p\) parameter power dan \(d(\cdot)\) jarak (Wackernagel, 2003).

2.4.2 Ordinary Kriging

Ordinary Kriging memodelkan korelasi spasial melalui variogram: \[ \gamma(h)=\frac{1}{2}\mathrm{Var}\big(Z(s)-Z(s+h)\big), \] dan prediksi dibentuk sebagai kombinasi linier tak-bias: \[ \hat z(s_0)=\sum_{i=1}^{n}\lambda_i z(s_i),\quad \sum_i\lambda_i=1, \] dengan \(\lambda_i\) ditentukan dari sistem Kriging berdasarkan model variogram (Cressie, 1993).

Bab. 3 METODOLOGI PENELITIAN

3.1 Data

Penelitian ini menggunakan dua sumber data utama:

  1. Data spasial (shapefile) Shapefile batas administrasi kabupaten/kota Provinsi Jawa Barat dari RBI 50K 2023. Data ini digunakan sebagai dasar pembuatan peta tematik dan pembentukan struktur ketetanggaan (spatial weights) untuk analisis autokorelasi dan pemodelan spasial.

  2. Dataset berisi indikator filariasis dan variabel penjelas yang dianalisis yang didapat dari sumber resmi Badan Pusat Statistik Provinsi Jawa Barat dalam publikasi berjudul Profil Kesehatan Jawa Barat 2025, dengan keterangan variabel sebagai berikut:

    • Y (Jumlah): Jumlah Kasus Kronis Filariasis
    • Y (Rate): Rate Kasus Kronis Filariasis
    • X1: Persentase Sanitasi Layak
    • X2: Persentase Rumah Layak Huni
    • X3: Rasio Puskesmas terhadap Penduduk
    • X4: Persentase Penduduk Miskin

Tahap pengolahan data (pra-analisis) meliputi:

  • Pembersihan format numerik (mengatasi koma desimal menjadi titik).
  • Normalisasi nama kabupaten/kota agar konsisten antara shapefile dan Excel (mis. penyeragaman “KAB.”/“KABUPATEN”/“KOTA” dan spasi).
  • Join data spasial dan tabular menggunakan join_key hasil normalisasi.
  • Validasi hasil join (memastikan tidak ada nilai Y_rate yang hilang akibat mismatch nama wilayah).

3.2 Unit Spasial

Unit spasial penelitian adalah 27 kabupaten/kota di Provinsi Jawa Barat. Setiap unit spasial direpresentasikan sebagai poligon pada shapefile RBI 50K 2023.

Dalam analisis interpolasi, setiap poligon direpresentasikan sebagai satu titik observasi menggunakan titik representatif di dalam poligon (st_point_on_surface) agar titik selalu berada di dalam wilayah administratif.

3.3 Metode Analisis dan Alur Kerja Penelitian

Penelitian dilakukan melalui alur kerja berikut:

  1. Persiapan data

    • Memuat shapefile RBI 50K 2023 Jawa Barat.
    • Memuat data Excel indikator filariasis dan variabel penjelas.
    • Cleaning numerik dan normalisasi nama wilayah.
    • Join shapefile dan data Excel, lalu validasi kecocokan.
  2. Eksplorasi deskriptif dan pemetaan

    • Statistik ringkas (min, max, mean, kuartil) untuk Y_jumlah dan X1–X4.

    • Visualisasi histogram dan boxplot Y_jumlah.

    • Peta tematik interaktif (tmap mode view) untuk:

      • Y_jumlah, Y_rate
      • X1, X2, X3, X4
  3. Autokorelasi spasial global (menggunakan Y_rate)

    • Menyusun bobot spasial berbasis kontiguitas:

      • Utama: Queen contiguity
    • Menguji autokorelasi global:

      • Moran’s I
      • Geary’s C
    • Membuat Moran scatterplot untuk interpretasi arah hubungan (cluster vs dispersed).

  4. Autokorelasi spasial lokal (LISA)

    • Menghitung Local Moran’s I (LISA) menggunakan Y_rate terstandarisasi.

    • Mengklasifikasikan pola wilayah menjadi:

      • High–High, Low–Low (cluster)
      • High–Low, Low–High (outlier)
      • Not significant
    • Memetakan hasil LISA dalam peta interaktif.

  5. Hotspot / Coldspot (Getis–Ord Gi*)

    • Menghitung z-score Gi* dan p-value (dua sisi).
    • Mengklasifikasikan wilayah menjadi kategori hotspot/coldspot berdasarkan ambang signifikansi (p < 0.10, p < 0.05, p < 0.01).
    • Memetakan hasil hotspot/coldspot dalam peta interaktif.
  6. Sensitivitas bobot spasial

    • Membandingkan hasil Moran’s I untuk tiga definisi bobot:

      • Queen contiguity
      • Rook contiguity
      • kNN-4 (nearest neighbors)
    • Menilai apakah kesimpulan autokorelasi stabil terhadap perubahan bobot.

  7. Pemodelan spatial econometrics (Y_rate sebagai respon)

    • Menjalankan OLS sebagai baseline, termasuk:

      • Uji multikolinearitas (VIF)
      • Moran’s I residual OLS
      • LM tests (Rao Score) untuk indikasi kebutuhan model spasial
    • Mengestimasi model spasial berikut:

      • SLX, SAR, SEM, SDM, SDEM, SAC, GNS
    • (Opsional) Menghitung dampak (impacts) untuk model yang relevan (SAR/SDM/SAC/GNS).

  8. Pemilihan model terbaik

    • Membandingkan model spasial menggunakan:

      • AIC
      • BIC
      • LogLik
    • Menentukan model terbaik (utama berdasarkan AIC terendah) dan menguji autokorelasi residual model terbaik.

  9. Pemetaan hasil model terbaik

    • Membuat peta interaktif:

      • Prediksi (fitted) model terbaik
      • Residual model terbaik
    • Menyertakan popup ringkas (wilayah, Y_obs, Y_fitted, residual).

  10. Interpolasi spasial (tambahan)

  • Mengubah poligon menjadi titik representatif (centroid aman).

  • Membuat grid prediksi di Jawa Barat dan melakukan masking (hanya titik dalam boundary).

  • Melakukan:

    • IDW
    • Ordinary Kriging dengan tuning variogram (banding spherical vs exponential)
  • Validasi silang Leave-One-Out (LOO) dan membandingkan RMSE IDW vs Kriging.

  • Menyajikan peta prediksi dan peta varian (ketidakpastian) untuk Kriging.


Bab 4. HASIL DAN PEMBAHASAN

Catatan variabel model
- Y = Rate Kasus Kronis Filariasis
- X1 = Persentase Sanitasi Layak
- X2 = Persentase Rumah Layak Huni
- X3 = Rasio Puskesmas terhadap Penduduk
- X4 = Persentase Penduduk Miskin

4.1 Persiapan Paket, Data, dan Join Spasial

## Jumlah entitas spasial (SHP): 27
## Jumlah baris Excel: 27
##   total missing_Yrate
## 1    27             0

Berdasarkan hasil pemeriksaan, jumlah entitas spasial pada shapefile sama dengan jumlah observasi dataset (27 wilayah), serta tidak terdapat missing join (missing_Yrate = 0). Hal ini menunjukkan bahwa integrasi data spasial dan tabular telah berhasil sehingga analisis dapat dilanjutkan.

4.2 Statistik Deskriptif (EDA)

##     Y_jumlah          Y_rate              X1               X2       
##  Min.   : 0.000   Min.   :0.00000   Min.   : 34.56   Min.   :32.83  
##  1st Qu.: 4.000   1st Qu.:0.03064   1st Qu.: 83.67   1st Qu.:43.49  
##  Median : 6.000   Median :0.05104   Median : 95.78   Median :61.91  
##  Mean   : 9.704   Mean   :0.05794   Mean   : 89.45   Mean   :58.72  
##  3rd Qu.:15.500   3rd Qu.:0.08107   3rd Qu.: 98.55   3rd Qu.:72.64  
##  Max.   :29.000   Max.   :0.16405   Max.   :100.00   Max.   :85.78  
##        X3              X4        
##  Min.   :16.21   Min.   : 2.340  
##  1st Qu.:35.04   1st Qu.: 6.360  
##  Median :42.80   Median : 8.410  
##  Mean   :43.54   Mean   : 8.014  
##  3rd Qu.:53.06   3rd Qu.:10.185  
##  Max.   :66.42   Max.   :11.930

Berdasarkan hasil statistik deskriptif, jumlah kasus kronis filariasis (Y_jumlah) menunjukkan variasi antarwilayah, dengan beberapa kabupaten/kota memiliki kasus 0 dan wilayah lain mencapai hingga 29 kasus. Variasi juga terlihat pada Y_rate (min 0 hingga max 0.16405), sehingga patut diuji apakah terdapat pola spasial.

4.3 Pembentukan Bobot Spasial (Queen Contiguity)

## Rata-rata tetangga: 3.925926  | Min: 1  | Max: 8

Bobot spasial Queen contiguity digunakan karena mempertimbangkan perbatasan maupun titik sudut yang bersinggungan. Bobot ini relevan untuk data kabupaten/kota yang memiliki variasi bentuk wilayah.

4.4 Peta Distribusi Spasial

Peta distribusi spasial menunjukkan bahwa jumlah kasus kronis filariasis tidak tersebar merata di Provinsi Jawa Barat. Beberapa wilayah tampak memiliki konsentrasi kasus yang lebih tinggi dibandingkan wilayah lainnya, sehingga mengindikasikan potensi adanya pola spasial.

4.5 OLS Baseline dan Diagnostik Spasial

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.044760 -0.026955 -0.007218  0.022057  0.081844 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.0509175  0.0601620   0.846    0.406
## X1           0.0004647  0.0005466   0.850    0.404
## X2          -0.0002185  0.0005028  -0.435    0.668
## X3          -0.0009928  0.0006043  -1.643    0.115
## X4           0.0026838  0.0028642   0.937    0.359
## 
## Residual standard error: 0.03671 on 22 degrees of freedom
## Multiple R-squared:  0.1703, Adjusted R-squared:  0.01943 
## F-statistic: 1.129 on 4 and 22 DF,  p-value: 0.3686
##       X1       X2       X3       X4 
## 1.278364 1.425374 1.132313 1.109951
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## weights: lw_q
## 
## Moran I statistic standard deviate = 2.0712, p-value = 0.01917
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.17507357      -0.10020524       0.01766398
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## test weights: lw_q
## 
## RSlag = 1.1618, df = 1, p-value = 0.2811
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## test weights: lw_q
## 
## RSerr = 1.3351, df = 1, p-value = 0.2479
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## test weights: lw_q
## 
## adjRSlag = 0.063562, df = 1, p-value = 0.801
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_model)
## test weights: lw_q
## 
## adjRSerr = 0.23693, df = 1, p-value = 0.6264

Hasil estimasi model Ordinary Least Squares (OLS) menunjukkan bahwa seluruh variabel penjelas (X1–X4) tidak berpengaruh signifikan terhadap Y_rate kasus filariasis pada taraf signifikansi 5%. Selain itu, nilai Adjusted R-squared yang rendah mengindikasikan bahwa model OLS memiliki kemampuan penjelasan yang terbatas terhadap variasi data antarwilayah.

Meskipun demikian, uji Moran’s I terhadap residual OLS menunjukkan adanya autokorelasi spasial yang signifikan (p = 0,01917). Hal ini menandakan bahwa residual model OLS masih mengandung pola spasial, sehingga asumsi independensi residual pada model OLS tidak terpenuhi.

Oleh karena itu, meskipun model OLS tidak signifikan, keberadaan autokorelasi spasial pada residual menunjukkan bahwa pendekatan non-spasial belum memadai. Untuk menangkap ketergantungan spasial antarwilayah yang tersisa, analisis kemudian dilanjutkan menggunakan model ekonometrika spasial.

4.6 Uji Autokorelasi Spasial Global (Moran’s I & Geary’s C)

## 
##  Moran I test under randomisation
## 
## data:  Yv  
## weights: lw_q    
## 
## Moran I statistic standard deviate = 1.6282, p-value = 0.05174
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.18381007       -0.03846154        0.01863615
## 
##  Geary C test under randomisation
## 
## data:  Yv 
## weights: lw_q   
## 
## Geary C statistic standard deviate = 1.7717, p-value = 0.03823
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.68535304        1.00000000        0.03154205

Hasil uji Moran’s I menunjukkan nilai positif (I ≈ 0,184) dengan nilai p = 0,05174 yang mendekati taraf signifikansi 5%, sedangkan uji Geary’s C menghasilkan hasil signifikan pada α = 5% (p = 0,03823). Hal ini menunjukkan adanya autokorelasi spasial positif lemah, yaitu wilayah dengan rate kasus yang serupa cenderung berdekatan secara geografis.

4.7 Moran Scatterplot (Global Pattern)

Moran scatterplot menunjukkan kecenderungan garis regresi yang berkemiringan positif, konsisten dengan Moran’s I global yang bernilai positif. Hal ini mendukung adanya kecenderungan pengelompokan nilai Y_rate yang mirip, meskipun kekuatannya lemah.

4.8 Autokorelasi Spasial Lokal (LISA)

## 
## Not significant         Low-Low 
##              25               2

Hasil analisis LISA menunjukkan bahwa mayoritas wilayah tidak signifikan secara lokal, dengan ringkasan Not significant = 25 dan hanya Low–Low = 2. Artinya, terdapat klaster lokal terbatas berupa wilayah dengan rate filariasis rendah yang dikelilingi wilayah bertetangga dengan nilai rendah pula.

Peta LISA

Berdasarkan peta LISA dan ringkasan output, klaster signifikan bertipe Low–Low teridentifikasi pada Kabupaten Bandung dan Kota Bandung. Hal ini mengindikasikan adanya pengelompokan wilayah dengan rate kasus rendah pada area tersebut.

4.9 Getis–Ord Gi* (Hotspot / Coldspot)

## 
##      Not significant Cold Spot (p < 0.05) Cold Spot (p < 0.10) 
##                   24                    2                    1

Analisis Getis–Ord Gi* menunjukkan bahwa sebagian besar wilayah tidak signifikan (24 wilayah). Ditemukan Cold Spot (p < 0.05) sebanyak 2 wilayah dan Cold Spot (p < 0.10) sebanyak 1 wilayah, sementara tidak ditemukan hotspot signifikan. Ini menunjukkan konsentrasi spasial nilai rendah pada area tertentu, sedangkan wilayah bernilai tinggi tidak membentuk klaster kuat.

Peta Getis–Ord Gi*

Sesuai dengan hasil LISA, peta Gi* menunjukkan keberadaan cold spot yang terkonsentrasi di area Bandung. Hal ini memperkuat indikasi bahwa wilayah dengan rate rendah cenderung berdekatan pada area tersebut.

## Simple feature collection with 3 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 107.251 ymin: -7.309771 xmax: 107.9384 ymax: -6.812845
## Geodetic CRS:  WGS 84
##       join_key              hotspot      Gi_z       p_gi
## 1 KOTA BANDUNG Cold Spot (p < 0.05) -2.216871 0.02663191
## 2      BANDUNG Cold Spot (p < 0.05) -2.174471 0.02966976
## 3  KOTA CIMAHI Cold Spot (p < 0.10) -1.883786 0.05959393
##                         geometry
## 1 MULTIPOLYGON (((107.6238 -6...
## 2 MULTIPOLYGON (((107.9147 -6...
## 3 MULTIPOLYGON (((107.5708 -6...

4.11 Model Spasial Terbaik (SEM), Diagnostik Residual, dan Interpolasi (IDW vs Ordinary Kriging)

4.11.1 Estimasi Model Spasial dan Perbandingan (AIC/BIC/LogLik)

##   Model       AIC       BIC   LogLik
## 1   SEM -95.31806 -86.24720 54.65903
## 2   SAR -94.97028 -85.89942 54.48514
## 3   SAC -93.38232 -83.01562 54.69116
## 4   SLX -90.23484 -77.27647 55.11742
## 5   SDM -89.55009 -75.29589 55.77505
## 6  SDEM -89.24950 -74.99529 55.62475
## 7   GNS -87.56272 -72.01267 55.78136
## 
## Model terbaik (AIC terendah): SEM

Berdasarkan tabel perbandingan model, model terbaik dengan AIC terendah adalah SEM (Spatial Error Model). Pemilihan SEM menunjukkan bahwa dependensi spasial pada data lebih dominan muncul pada komponen error (gangguan), bukan pada lag variabel dependen secara langsung.

4.11.2 Diagnostik Residual Model Terpilih

best_res <- tryCatch(residuals(best_model), error = function(e) NULL) moran_best_res <- if (!is.null(best_res)) moran.test(best_res, lw_q, zero.policy = TRUE) else NULL moran_best_res

Hasil uji Moran’s I untuk residual model spasial terbaik menunjukkan p-value = 0.3014 (tidak signifikan), yang berarti bahwa SEM berhasil menghilangkan autokorelasi spasial yang sebelumnya masih tersisa pada residual OLS. Dengan demikian, model SEM dapat dianggap lebih sesuai dibandingkan OLS.

4.11.3 Peta Residual (Model Terpilih: SEM)

Peta residual membantu melihat wilayah mana yang masih memiliki deviasi relatif besar antara nilai observasi dan prediksi model. Karena residual model SEM tidak lagi berautokorelasi secara spasial, pola residual yang muncul tidak menunjukkan klaster spasial yang konsisten.

4.12 Interpolasi Spasial

Interpolasi spasial dilakukan untuk menggambarkan pola kontinu sebaran rate kasus kronis filariasis (Y_rate) di Provinsi Jawa Barat berdasarkan data titik kabupaten/kota. Metode yang digunakan meliputi Inverse Distance Weighting (IDW) dan Ordinary Kriging (OK), di mana Ordinary Kriging juga menghasilkan informasi ketidakpastian prediksi (variansi kriging). Ketiga peta interpolasi tersebut memberikan gambaran yang saling melengkapi mengenai tingkat risiko dan keandalan estimasi spasial.

## [inverse distance weighted interpolation]
## [using ordinary kriging]

4.12.1 Variogram

## 
## RMSE variogram Sph: 0.03923
## RMSE variogram Exp: 0.03909
## Variogram terpilih : Exp

Grafik variogram empiris menggambarkan hubungan antara jarak (sumbu X) dan semivarians (sumbu Y). Semivarians yang meningkat seiring jarak menunjukkan bahwa lokasi yang semakin jauh cenderung semakin tidak mirip nilainya. Titik-titik pada variogram empiris adalah ringkasan variasi data pada beberapa “kelas jarak” (bin) yang ditentukan oleh cutoff dan width (tuned). Tuning ini penting karena menentukan seberapa jauh jarak maksimum yang dipakai dan bagaimana jarak dikelompokkan, sehingga pola semivarians bisa terbaca lebih stabil.

Pada grafik fitting variogram — model terpilih, garis kurva menunjukkan model teoretis (misalnya Spherical/Exponential) yang dipasang untuk mengikuti titik-titik empiris. Secara konsep, ada tiga komponen kunci: 1. Nugget: nilai semivarians saat jarak mendekati 0 (atau intersep awal). Nugget merepresentasikan variasi mikro-skala dan/atau error pengukuran—artinya ada variasi yang tidak bisa dijelaskan hanya oleh kedekatan spasial. 2. Sill: nilai plateau/ambang atas yang didekati kurva ketika jarak makin besar. Sill merepresentasikan total variasi yang “stabil” ketika lokasi sudah tidak lagi berkorelasi kuat. 3. Range: jarak ketika kurva mendekati sill (atau ketika korelasi spasial praktis hilang). Di bawah range, lokasi relatif masih berkorelasi; di atas range, tambahan jarak tidak banyak mengubah kemiripan (nilai dianggap hampir independen).

Dari bentuk kurva yang naik cepat di jarak dekat lalu melandai, interpretasinya adalah: ketergantungan spasial paling kuat terjadi pada jarak pendek, kemudian melemah sampai akhirnya cenderung stabil. Model variogram terpilih inilah yang menjadi “mesin” OK untuk menentukan bobot dan sekaligus menghasilkan prediksi serta variansi prediksi. Jika fitting cukup baik (kurva masuk akal mengikuti pola titik), maka prediksi kriging biasanya lebih reliabel; jika titik-titik empiris menyebar jauh dari kurva, interpretasi dan hasil kriging perlu lebih hati-hati (misalnya mempertimbangkan model variogram lain, transformasi data, atau penambahan titik).

4.12.2 Cross Validation

##       ME      MAE     RMSE     MSSE 
## 0.004366 0.029412 0.041272      NaN
##       ME      MAE     RMSE     MSSE 
## 0.002542 0.029866 0.039090 1.009781

Hasil cross-validation menunjukkan bahwa Ordinary Kriging (OK) memiliki kinerja lebih baik dibanding IDW karena RMSE lebih kecil (OK RMSE = 0.03909 vs IDW RMSE = 0.04127). Selain itu, MSSE OK ≈ 1.0098 mendekati 1, yang mengindikasikan bahwa estimasi ketidakpastian kriging terkalibrasi dengan baik. Nilai MSSE pada IDW tidak tersedia (NaN) karena IDW tidak menyediakan varian prediksi.

4.12.3 Interpolasi Ordinary Kriging — Prediksi Rate

Peta prediksi Ordinary Kriging menunjukkan bahwa nilai Y_rate filariasis di Provinsi Jawa Barat memiliki pola spasial yang tidak merata dan membentuk beberapa zona konsentrasi. Area dengan nilai prediksi lebih tinggi (ditunjukkan oleh warna hangat seperti kuning tua hingga merah) tampak terkonsentrasi pada beberapa wilayah tertentu, yang mengindikasikan potensi wilayah berisiko tinggi terhadap kasus filariasis kronis. Sebaliknya, wilayah dengan nilai prediksi rendah (warna hijau hingga biru) mencerminkan area dengan tingkat kejadian filariasis yang relatif lebih rendah.

Metode Ordinary Kriging menghasilkan permukaan prediksi yang relatif halus karena mempertimbangkan struktur korelasi spasial antarwilayah melalui variogram. Hal ini menyebabkan transisi nilai antarwilayah berlangsung lebih gradual dan mencerminkan asumsi adanya ketergantungan spasial pada data Y_rate. Dengan demikian, peta ini memberikan estimasi spasial yang secara teoritis lebih konsisten dibandingkan metode berbasis jarak semata.

4.12.4 Interpolasi Ordinary Kriging — Variansi Prediksi (Ketidakpastian)

Peta variansi prediksi Ordinary Kriging menggambarkan tingkat ketidakpastian estimasi pada setiap lokasi interpolasi. Nilai variansi yang lebih tinggi (warna merah tua) umumnya muncul pada wilayah yang lebih jauh dari titik observasi atau berada di bagian tepi wilayah studi. Hal ini menunjukkan bahwa prediksi di area tersebut memiliki tingkat kepercayaan yang lebih rendah akibat keterbatasan informasi data di sekitarnya.

Sebaliknya, wilayah dengan variansi prediksi yang lebih rendah (warna kuning terang) menunjukkan area yang relatif dekat dengan lokasi pengamatan, sehingga estimasi nilai Y_rate lebih andal. Informasi variansi ini penting dalam analisis spasial karena memungkinkan peneliti tidak hanya menilai besarnya prediksi, tetapi juga memahami tingkat keandalan hasil interpolasi, terutama dalam konteks perencanaan kebijakan kesehatan berbasis wilayah.

4.12.5 Interpolasi IDW — Prediksi Rate

Peta interpolasi menggunakan metode Inverse Distance Weighting (IDW) menunjukkan pola sebaran Y_rate yang cenderung lebih lokal dan tajam dibandingkan Ordinary Kriging. Nilai prediksi yang tinggi terlihat sangat dipengaruhi oleh keberadaan titik observasi dengan nilai Y_rate besar, sehingga terbentuk zona-zona konsentrasi yang lebih kontras di sekitar titik tersebut.

Karakteristik IDW yang hanya mempertimbangkan jarak tanpa memperhitungkan struktur korelasi spasial menyebabkan hasil interpolasi menjadi lebih sensitif terhadap nilai ekstrem lokal. Akibatnya, variasi spasial yang dihasilkan oleh IDW tampak lebih kasar dan kurang halus dibandingkan Ordinary Kriging. Meskipun demikian, IDW tetap berguna sebagai pendekatan eksploratif awal karena sifatnya yang sederhana dan intuitif.

4.12.6 Perbandingan dan Implikasi Analisis

Secara keseluruhan, Ordinary Kriging menghasilkan permukaan prediksi yang lebih stabil dan disertai informasi ketidakpastian, sehingga lebih sesuai digunakan sebagai dasar analisis lanjutan dan pengambilan keputusan. Sementara itu, IDW memberikan gambaran yang lebih menekankan pengaruh lokal dari titik pengamatan, namun kurang mempertimbangkan struktur spasial global.

Keberadaan peta variansi pada Ordinary Kriging memberikan nilai tambah yang penting, karena memungkinkan identifikasi wilayah yang memerlukan kehati-hatian lebih tinggi dalam interpretasi hasil. Dengan demikian, kombinasi hasil interpolasi ini memberikan pemahaman yang komprehensif mengenai pola spasial filariasis di Jawa Barat, baik dari sisi estimasi nilai maupun tingkat kepercayaannya.

Bab 5. KESIMPULAN DAN SARAN

5.1 Kesimpulan

Berdasarkan seluruh hasil analisis spasial kabupaten/kota di Provinsi Jawa Barat (n = 27) dengan variabel respon Rate Kasus Kronis Filariasis (Y_rate) serta variabel penjelas X1–X4, diperoleh kesimpulan berikut:

  1. Model OLS sebagai baseline belum mampu menjelaskan variasi Y_rate dengan baik. Nilai Adjusted R-squared sebesar 0,01943 (R² = 0,1703) dan uji F tidak signifikan (p-value = 0,3686). Secara parsial, tidak ada variabel yang signifikan pada taraf 5% pada model OLS.

  2. Tidak terdapat indikasi multikolinearitas yang bermasalah pada variabel penjelas. Seluruh nilai VIF berada pada rentang rendah (X1 = 1,278; X2 = 1,425; X3 = 1,132; X4 = 1,110), sehingga hubungan antar variabel independen tidak mengganggu estimasi secara serius.

  3. Residual OLS masih mengandung autokorelasi spasial yang signifikan. Uji Moran’s I untuk residual OLS menunjukkan z = 2,0712 dengan p-value = 0,01917, sehingga terdapat ketergantungan spasial yang belum ditangkap oleh model OLS. Temuan ini menguatkan alasan penggunaan model spasial.

  4. Uji LM (Rao’s score) tidak memberikan sinyal kuat untuk memilih SAR atau SEM secara langsung, karena baik uji lag maupun error tidak signifikan (RSlag p = 0,2811; RSerr p = 0,2479), termasuk versi adjusted. Namun, mengingat Moran residual OLS signifikan, pemodelan spasial tetap relevan untuk dicoba dan dibandingkan.

  5. Model spasial terbaik berdasarkan AIC adalah SEM (Spatial Error Model). Dari perbandingan AIC, model SEM memperoleh nilai AIC terendah (AIC = -95,318) dibanding model lain (SAR, SAC, SLX, SDM, SDEM, GNS). Ini menunjukkan bahwa ketergantungan spasial pada kasus ini lebih efektif ditangkap melalui komponen error spasial.

  6. Setelah menggunakan model terbaik (SEM), autokorelasi spasial residual menjadi tidak signifikan. Moran’s I residual SEM menghasilkan p-value = 0,3014, yang mengindikasikan bahwa model SEM sudah cukup mengakomodasi struktur spasial yang sebelumnya muncul pada residual OLS.

  7. Interpolasi spasial menunjukkan Ordinary Kriging lebih akurat dibanding IDW untuk membentuk permukaan prediksi Y_rate. Hasil validasi silang LOO menunjukkan RMSE IDW = 0,041272 sedangkan RMSE Ordinary Kriging (model spherical) = 0,039230, sehingga Kriging memberikan galat prediksi yang lebih kecil.

5.2 Saran

Beberapa saran untuk pengembangan analisis dan peningkatan kualitas hasil penelitian adalah:

  1. Perkuat kualitas dan konsistensi definisi Y_rate. Pastikan satuan rate (misalnya per 1.000/10.000 penduduk) konsisten antar wilayah dan tahun pengamatan jelas, karena interpretasi spasial sangat bergantung pada validitas definisi indikator.

  2. Tambahkan variabel penjelas yang lebih dekat dengan mekanisme kejadian filariasis. Misalnya kepadatan penduduk, kualitas drainase/lingkungan, curah hujan, elevasi, proporsi wilayah rawan genangan, kepadatan vektor, atau indikator akses layanan kesehatan yang lebih rinci.

  3. Lakukan uji sensitivitas tambahan pada spesifikasi model. Misalnya mencoba transformasi Y_rate (log/Box-Cox bila sesuai), menguji outlier, serta uji heteroskedastisitas/robust standard errors agar interpretasi koefisien lebih stabil.

  4. Uji sensitivitas resolusi grid interpolasi dan parameter model. Untuk interpolasi, resolusi grid (misalnya 2 km, 5 km, 10 km) dan parameter IDW (idp, nmax) maupun variogram pada Kriging sebaiknya diuji agar hasil peta prediksi tidak bergantung pada satu pilihan parameter saja.

  5. Laporkan peta fitted dan residual sebagai bahan evaluasi model. Peta fitted membantu melihat wilayah dengan prediksi tinggi/rendah, sedangkan peta residual membantu mengidentifikasi wilayah yang masih “menyimpang” dari model (under/over-prediction) untuk evaluasi lanjutan (misalnya variabel yang belum dimasukkan).

Daftar Pustaka

  1. Anselin, L. (1988). Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers.
  2. Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
  3. Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). New York: Springer.
  4. Cliff, A. D., & Ord, J. K. (1981). Spatial Processes: Models & Applications. London: Pion.
  5. Cressie, N. (1993). Statistics for Spatial Data (Revised ed.). New York: Wiley.
  6. Elhorst, J. P. (2014). Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Heidelberg: Springer.
  7. Getis, A., & Ord, J. K. (1992). The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis, 24(3), 189–206.
  8. LeSage, J., & Pace, R. K. (2009). Introduction to Spatial Econometrics. Boca Raton: CRC Press.
  9. Ord, J. K., & Getis, A. (1995). Local Spatial Autocorrelation Statistics: Distributional Issues and an Application. Geographical Analysis, 27(4), 286–306.
  10. Wackernagel, H. (2003). Multivariate Geostatistics: An Introduction with Applications (3rd ed.). Berlin: Springer.

Lampiran

Dashboard Analisis Spasial

Dashboard dapat diakses melalui: 🔗 Klik untuk membuka dashboard