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.
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:
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:
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:
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.
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).
Autokorelasi spasial mengukur kemiripan nilai antar unit spasial berdasarkan bobot \(\mathbf{W}\).
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).
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).
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).
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).
Matriks bobot spasial \(\mathbf{W}\) merepresentasikan struktur ketetanggaan antarwilayah. Pada data poligon kabupaten/kota, bobot sering dibangun dari kontiguitas:
Bobot kemudian dinormalisasi (mis. row-standardized, style = “W”) agar jumlah bobot per baris menjadi 1, sehingga interpretasi spatial lag lebih stabil.
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}\).
\[ \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).
\[ \mathbf{y} = \mathbf{X}\beta + \mathbf{WX}\theta + \varepsilon. \]
Model SLX menangkap spillover kovariat dari wilayah tetangga melalui \(\mathbf{WX}\) (LeSage & Pace, 2009).
\[ \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).
\[ \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).
\[ \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).
\[ \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).
\[ \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).
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).
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).
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).
Penelitian ini menggunakan dua sumber data utama:
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.
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:
Tahap pengolahan data (pra-analisis) meliputi:
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.
Penelitian dilakukan melalui alur kerja berikut:
Persiapan data
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:
Autokorelasi spasial global (menggunakan Y_rate)
Menyusun bobot spasial berbasis kontiguitas:
Menguji autokorelasi global:
Membuat Moran scatterplot untuk interpretasi arah hubungan (cluster vs dispersed).
Autokorelasi spasial lokal (LISA)
Menghitung Local Moran’s I (LISA) menggunakan Y_rate terstandarisasi.
Mengklasifikasikan pola wilayah menjadi:
Memetakan hasil LISA dalam peta interaktif.
Hotspot / Coldspot (Getis–Ord Gi*)
Sensitivitas bobot spasial
Membandingkan hasil Moran’s I untuk tiga definisi bobot:
Menilai apakah kesimpulan autokorelasi stabil terhadap perubahan bobot.
Pemodelan spatial econometrics (Y_rate sebagai respon)
Menjalankan OLS sebagai baseline, termasuk:
Mengestimasi model spasial berikut:
(Opsional) Menghitung dampak (impacts) untuk model yang relevan (SAR/SDM/SAC/GNS).
Pemilihan model terbaik
Membandingkan model spasial menggunakan:
Menentukan model terbaik (utama berdasarkan AIC terendah) dan menguji autokorelasi residual model terbaik.
Pemetaan hasil model terbaik
Membuat peta interaktif:
Menyertakan popup ringkas (wilayah, Y_obs, Y_fitted, residual).
Interpolasi spasial (tambahan)
Mengubah poligon menjadi titik representatif (centroid aman).
Membuat grid prediksi di Jawa Barat dan melakukan masking (hanya titik dalam boundary).
Melakukan:
Validasi silang Leave-One-Out (LOO) dan membandingkan RMSE IDW vs Kriging.
Menyajikan peta prediksi dan peta varian (ketidakpastian) untuk Kriging.
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
## 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.
## 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.
## 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.
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.
##
## 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.
##
## 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.
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.
##
## 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.
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.
##
## 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.
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...
## 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.
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.
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.
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]
##
## 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).
## 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.
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.
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.
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.
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.
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:
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.
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.
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.
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.
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.
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.
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.
Beberapa saran untuk pengembangan analisis dan peningkatan kualitas hasil penelitian adalah:
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.
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.
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.
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.
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).
Dashboard dapat diakses melalui: 🔗 Klik untuk membuka dashboard