Analisis Spasial Kasus Tuberkulosis di Kabupaten/Kota Provinsi Jawa Timur Tahun 2024
Gusti Ayu Anisa Eldina - 140610230042
2025-11-01
- 1 Pendahuluan
- 2 Tinjauan Pustaka
- 2.1 Teori Spatial Dependence
- 2.2 Spatial Autocorrelation dan Spatial Heterogeneity
- 2.3 Model Spatial Econometrics
- 2.3.1 Model Dasar Regresi Spasial
- 2.3.2 Spatial Autoregressive Model (SAR)
- 2.3.3 Spatial Error Model (SEM)
- 2.3.4 Spatial Durbin Model (SDM)
- 2.3.5 Spatial Durbin Error Model (SDEM)
- 2.3.6 Spatial Autoregressive Combined Model (SAC)
- 2.3.7 Estimasi dan Uji Model
- 2.3.8 Interpretasi Efek Spasial (Direct, Indirect, dan Total)
- 2.3.9 Pemilihan Model Terbaik (AIC, BIC, dan Uji LM Lag/Error)
- 3 Metodologi Penelitian
- 4 Hasil dan Pembahasan
- 5 Kesimpulan dan Saran
- 6 Daftar Pustaka
https://bit.ly/DashboardSpasialTBC
Abstrak
Tuberkulosis (TBC) masih menjadi masalah kesehatan utama di Jawa Timur dengan distribusi kasus yang menunjukkan indikasi ketergantungan spasial antarwilayah. Penelitian ini menganalisis pola persebaran TBC menggunakan pendekatan spasial ekonometrika, meliputi analisis autokorelasi spasial (Moran’s I, Geary’s C, dan LISA), identifikasi hotspot (Getis–Ord Gi*), serta pemodelan ekonometrika spasial (OLS, SAR, SEM, SDM, SDEM, dan SAC). Hasil menunjukkan adanya autokorelasi spasial positif dengan klaster High-High di wilayah padat penduduk seperti Surabaya dan Sidoarjo. Model SDEM terpilih sebagai model terbaik (AIC = 597.43) dengan variabel signifikan meliputi jumlah penduduk, kepadatan penduduk, dan jumlah rumah sakit. Temuan ini menegaskan adanya efek limpahan (spill-over effect) antarwilayah terhadap penyebaran TBC. Penelitian ini menyimpulkan bahwa pendekatan spasial penting untuk perencanaan intervensi kesehatan yang lebih terarah dan berbasis wilayah di Jawa Timur.
1 Pendahuluan
1.1 Latar Belakang
Tuberkulosis (TBC) masih menjadi darurat kesehatan global dan nasional. Global Tuberculosis Report 2024 mencatat bahwa sebanyak 10,8 juta kasus baru terjadi di dunia, dengan Indonesia menempati peringkat kedua setelah India, menyumbang sekitar 10 persen dari total kasus global (World Health Organization [WHO], 2024). Angka kejadian TBC di Indonesia bahkan mencapai 394 kasus per 100.000 penduduk (World Bank, 2024) yang jauh di atas ambang eliminasi TBC (≤10 per 100.000). Kondisi ini menjadikan Indonesia sebagai salah satu episentrum TBC dunia.
Lebih memprihatinkan, pandemi COVID-19 memperburuk upaya pengendalian TBC. Selama 2020–2021, deteksi kasus turun hampir 20 persen, menyebabkan ribuan pasien tidak terdiagnosis, sementara pada 2022–2023 jumlah kasus kembali melonjak hingga 30 persen dibanding periode pra-pandemi (Palupi et al., 2023). Fenomena rebound ini menunjukkan lemahnya deteksi dini dan potensi penularan laten di komunitas yang sulit dijangkau layanan kesehatan.
Di Provinsi Jawa Timur, beban penyakit tergolong tertinggi secara nasional dengan lebih dari 85.000 kasus tercatat pada tahun 2023 (Dinas Kesehatan Jatim, 2024). Pola kasus tidak menyebar merata: kota-kota padat seperti Surabaya, Sidoarjo, dan Malang menjadi hotspot, sedangkan wilayah pedalaman menunjukkan kasus tersembunyi akibat under-reporting. Fenomena ini menandakan adanya autokorelasi spasial positif, di mana wilayah berisiko tinggi saling berdekatan dan membentuk klaster TBC (Fahdhienie et al., 2024).
Secara nasional, peta epidemiologis memperlihatkan ketimpangan spasial yang semakin tajam. Penelitian Helmy et al. (2023) menemukan bahwa persebaran TBC di Indonesia dipengaruhi oleh faktor lingkungan dan sosial-ekonomi khususnya kepadatan penduduk, kemiskinan, dan akses rumah sakit dengan pola spasial yang signifikan (Moran’s I > 0,2). Artinya, penyakit ini tidak hanya masalah individu, tetapi masalah struktural yang terkunci dalam ruang dan ketimpangan pembangunan.
Fenomena di lapangan menunjukkan bahwa penularan TBC bersifat spasial yaitu terpusat di wilayah-wilayah dengan karakteristik sosial dan fisik tertentu. Di permukiman padat perkotaan, ventilasi rumah yang buruk, jarak antarpenghuni sempit, serta mobilitas tinggi mempercepat penularan Mycobacterium tuberculosis (Curtis & Kar, 2013). Sebaliknya, di daerah miskin, keterbatasan gizi dan akses kesehatan memperpanjang rantai penularan karena pasien terlambat berobat (Solar & Irwin, 2010).
Temuan surveilans di Surabaya (Andini et al., 2023) bahkan menunjukkan bahwa kepadatan penduduk memiliki korelasi r = 0,65 terhadap insiden TBC AFB+, menegaskan peran lingkungan sosial-spasial dalam memperkuat beban penyakit. Dengan kata lain, peta TBC adalah juga peta ketimpangan sosial dan pembangunan wilayah.
Kepadatan penduduk per kilometer persegi merupakan indikator utama dari teori density-dependent transmission. Penyakit yang menyebar lewat udara seperti TBC akan meningkat eksponensial pada area berpenduduk padat (Curtis & Kar, 2013). Kepadatan tidak hanya menggambarkan banyaknya orang, tetapi juga mengindikasikan tekanan infrastruktur, sanitasi, dan kualitas hunian yang memperburuk kerentanan.
Jumlah rumah sakit menggambarkan ketersediaan fasilitas kesehatan yang menentukan probabilitas masyarakat untuk mengakses layanan diagnosis dan terapi, sesuai dengan Andersen Behavioral Model of Health Services Use (Andersen, 1995). Namun secara spasial, fasilitas yang terpusat di perkotaan dapat menciptakan kesenjangan akses bagi wilayah pinggiran. Jumlah rumah sakit di suatu kabupaten menjadi proksi ketimpangan infrastruktur kesehatan yang berkontribusi pada variasi kasus.
Perilaku merokok, terutama kretek tanpa filter yang lazim di Indonesia, meningkatkan risiko infeksi TBC hingga 2,5 kali lipat (The Union, 2022). Rokok merusak fungsi imun paru, mempermudah kolonisasi Mycobacterium, serta memperpanjang fase penularan (Mishra et al., 2020). Selain itu, distribusi konsumsi rokok juga bersifat spasial yaitu terkonsentrasi di wilayah miskin dan pekerja informal sehingga penting dianalisis dalam konteks spasial.
Kemiskinan berperan besar dalam memperparah determinan kesehatan melalui gizi rendah, sanitasi buruk, serta keterbatasan akses layanan medis. Berdasarkan Social Determinants of Health Framework (Solar & Irwin, 2010), wilayah dengan tingkat kemiskinan tinggi cenderung membentuk hotspot karena faktor-faktor sosial yang saling memperkuat. Helmy et al. (2023) menunjukkan korelasi positif signifikan antara tingkat kemiskinan dan kasus TBC (p < 0,05).
Fenomena di atas menggambarkan kegagalan deteksi spasial dan ketimpangan struktural dalam pengendalian TBC. Selama ini, pendekatan kebijakan cenderung bersifat agregat dengan menyamaratakan strategi antarwilayah tanpa memperhitungkan keterkaitan spasial antar kabupaten/kota. Akibatnya, wilayah berisiko tinggi sering menjadi reservoir penularan baru bagi wilayah sekitarnya, memperpanjang siklus infeksi laten.
Belum banyak penelitian yang mengintegrasikan analisis spasial dan ekonometrika untuk memahami efek lintas-wilayah (spill-over effect) dalam kasus TBC di tingkat kabupaten/kota. Padahal, tanpa pemahaman spasial yang presisi, kebijakan intervensi cenderung tidak tepat sasaran seolah menutup satu pintu tetapi membiarkan jendela lain terbuka.
Oleh sebab itu, analisis spasial ekonometrik terhadap kasus TBC di Jawa Timur menjadi mendesak, bukan sekadar untuk pemetaan statistik, tetapi sebagai dasar ilmiah dalam memutus rantai penularan lintas-wilayah, menekan ketimpangan akses kesehatan, serta mempercepat pencapaian target Eliminasi TBC 2030.
1.2 Identifikasi Masalah
Meskipun Indonesia telah menargetkan eliminasi tuberkulosis (TBC) pada tahun 2030, realitas di lapangan menunjukkan bahwa penyakit ini masih menjadi masalah kesehatan masyarakat yang serius, terutama di Provinsi Jawa Timur. Berdasarkan data Dinas Kesehatan Provinsi Jawa Timur (2024), provinsi ini mencatat lebih dari 85.000 kasus TBC pada tahun 2023, menjadikannya salah satu penyumbang kasus tertinggi secara nasional. Fenomena ini memperlihatkan bahwa upaya eliminasi belum berjalan efektif di tingkat daerah. Ketimpangan antardaerah dalam capaian deteksi dini, ketersediaan fasilitas kesehatan, serta distribusi sosial-ekonomi menjadi hambatan utama yang hingga kini belum tertangani secara menyeluruh.
Masalah semakin kompleks karena pola penyebaran TBC bersifat tidak acak, melainkan membentuk pola spasial yang menunjukkan adanya keterkaitan antarwilayah. Penelitian Fahdhienie et al. (2024) dan Helmy et al. (2023) menunjukkan bahwa kasus TBC di Indonesia dan di beberapa provinsi cenderung terkonsentrasi di wilayah tertentu, membentuk cluster “tinggi-tinggi” (High-High) di daerah padat penduduk dan miskin, sementara daerah lain relatif rendah. Pola autokorelasi spasial positif ini mengindikasikan bahwa intervensi di satu wilayah berpotensi memengaruhi wilayah di sekitarnya. Pendekatan yang tidak mempertimbangkan hubungan geografis ini menyebabkan program intervensi sering kali bersifat seragam dan kurang efektif di daerah dengan karakteristik berbeda.
Penelitian terdahulu umumnya hanya berfokus pada faktor-faktor sosial ekonomi atau demografis secara parsial tanpa mempertimbangkan efek keterkaitan spasial (spatial dependency). Misalnya, studi Andini et al. (2023) di Kota Surabaya menemukan hubungan signifikan antara kepadatan penduduk dan insiden TBC, namun tidak mengeksplorasi bagaimana hubungan tersebut berdampak pada wilayah sekitar. Studi Helmy et al. (2023) telah menggunakan pemodelan spasial di tingkat nasional, tetapi berskala makro dan belum spesifik menjelaskan dinamika antar kabupaten/kota di satu provinsi yang memiliki heterogenitas tinggi seperti Jawa Timur. Sementara itu, penelitian yang mengaitkan variabel sosial, perilaku, dan infrastruktur secara simultan dengan pendekatan ekonometrika spasial di level kabupaten/kota masih terbatas di Indonesia.
Selain itu, ketimpangan spasial dalam penyediaan layanan kesehatan turut memperburuk situasi. Rumah sakit dan fasilitas deteksi dini TBC banyak terkonsentrasi di wilayah perkotaan seperti Surabaya dan Malang, sementara daerah dengan akses terbatas seperti Bondowoso dan Sampang masih menghadapi hambatan logistik dan tenaga medis (RDI Global, 2023). Ketidakseimbangan ini menimbulkan fenomena under-detection, bukan karena penularan rendah, tetapi karena kemampuan deteksi yang terbatas. Dalam konteks spasial, hal ini berpotensi menciptakan hidden cluster yang luput dari perhatian pembuat kebijakan.
Gap penelitian juga muncul dari aspek metodologis. Sebagian besar studi TBC di Indonesia masih menggunakan analisis regresi linier biasa tanpa mempertimbangkan spatial lag atau spatial error, padahal penyakit menular seperti TBC secara teoretis memiliki mekanisme penyebaran spasial yang kuat melalui mobilitas penduduk dan keterhubungan sosial (Bhatt et al., 2022). Pendekatan ekonometrika spasial seperti SAR (Spatial Autoregressive Model), SEM (Spatial Error Model), dan SDEM (Spatial Durbin Error Model) mampu menangkap efek lintas-wilayah (spill-over effect) yang tidak dapat dijelaskan oleh regresi konvensional. Dengan demikian, masih terdapat celah penelitian dalam memahami bagaimana faktor-faktor sosial, ekonomi, perilaku, dan infrastruktur berinteraksi secara spasial terhadap variasi kasus TBC di Jawa Timur.
Dari sisi empiris, fenomena kasus nyata juga memperlihatkan kontradiksi yang menarik. Wilayah dengan fasilitas kesehatan tinggi seperti Surabaya justru masih mencatat insiden TBC yang besar (Palupi et al., 2023), sementara wilayah dengan keterbatasan fasilitas justru memiliki angka laporan kasus rendah, bukan karena bebas TBC, tetapi karena undetected burden. Berdasarkan kondisi tersebut, dapat diidentifikasi beberapa masalah utama:
Kasus TBC di Jawa Timur masih tinggi dan menunjukkan pola spasial yang tidak merata antarwilayah.
Upaya kebijakan dan penelitian sebelumnya belum secara komprehensif mengintegrasikan pendekatan spasial dalam analisis faktor sosial, ekonomi, dan infrastruktur.
Belum ada penelitian yang secara spesifik menilai spill-over effect antar kabupaten/kota dalam konteks TBC menggunakan pendekatan ekonometrika spasial (SAR, SEM, SDM, SDEM, SAC) di tingkat provinsi.
Ketimpangan spasial dalam akses fasilitas kesehatan menyebabkan bias pelaporan kasus dan memperkuat potensi hidden cluster.
Masih terbatasnya pemahaman empiris tentang interaksi antara kepadatan penduduk, perilaku kesehatan (seperti merokok), serta tingkat kemiskinan terhadap risiko spasial TBC.
Dengan demikian, identifikasi masalah dalam penelitian ini tidak hanya berfokus pada tingginya angka kasus TBC, tetapi juga pada pemahaman mendalam mengenai pola spasial, ketimpangan wilayah, dan keterkaitan antar kabupaten/kota di Jawa Timur. Penelitian ini berupaya menutup gap tersebut dengan pendekatan spasial ekonometrik yang komprehensif, sehingga hasilnya dapat digunakan untuk perencanaan intervensi yang lebih terarah, berbasis bukti ilmiah, dan relevan dengan konteks geografis daerah.
1.3 Tujuan Penelitian
Berdasarkan permasalahan yang telah diuraikan pada bagian sebelumnya, penelitian ini merumuskan tujuan penelitian untuk menjawab persoalan tingginya beban tuberkulosis (TBC) di Provinsi Jawa Timur yang menunjukkan pola spasial tidak merata dan dipengaruhi oleh berbagai faktor sosial, ekonomi, perilaku, serta infrastruktur kesehatan. Dalam konteks tersebut, penelitian ini bertujuan untuk menganalisis keterkaitan spasial antarwilayah dan mengidentifikasi determinan utama yang berkontribusi terhadap variasi kasus TBC antar kabupaten/kota.
Secara lebih rinci, tujuan penelitian ini adalah sebagai berikut:
Menganalisis pola persebaran spasial kasus TBC di Provinsi Jawa Timur dan mengidentifikasi adanya autokorelasi spasial melalui pendekatan statistik spasial seperti Moran’s I, Geary’s C, dan Local Indicators of Spatial Association (LISA).
Menggambarkan fenomena klaster atau hotspot TBC di Jawa Timur menggunakan analisis Getis–Ord Gi* sebagai dasar untuk menentukan wilayah prioritas intervensi kesehatan.
Menguji pengaruh variabel sosial-ekonomi dan perilaku terhadap variasi spasial kasus TBC, meliputi kepadatan penduduk, jumlah rumah sakit, konsumsi rokok kretek tanpa filter, serta jumlah penduduk miskin.
Menerapkan model ekonometrika spasial (OLS, SAR, SEM, SDM, SDEM, dan SAC) untuk membandingkan efektivitas masing-masing model dalam menjelaskan hubungan antarwilayah (spatial dependency).
Menentukan model terbaik berdasarkan kriteria statistik (AIC, BIC, dan LM-test) serta menilai dampak langsung, tidak langsung, dan total (spatial impacts) antarwilayah dalam penyebaran TBC.
Menyajikan peta hasil prediksi dan residual dari model terbaik untuk mengidentifikasi potensi spill-over effect dan pola kerentanan wilayah yang bersifat spasial.
Melalui tujuan-tujuan tersebut, penelitian ini diharapkan mampu menghasilkan pemahaman mengenai dinamika spasial TBC di Jawa Timur dan mendukung perumusan kebijakan berbasis bukti (evidence-based policy) yang lebih tepat sasaran dalam rangka percepatan eliminasi TBC tahun 2030 sebagaimana ditetapkan dalam Rencana Aksi Nasional Eliminasi Tuberkulosis 2024–2030 (Kementerian Kesehatan RI, 2023).
1.4 Batasan Penelitian
Penelitian ini memiliki beberapa batasan yang ditetapkan untuk menjaga fokus analisis, keakuratan interpretasi, dan konsistensi metodologis, tanpa mengurangi relevansi serta kekuatan hasil penelitian terhadap isu tuberkulosis (TBC) di Provinsi Jawa Timur.
Pertama, penelitian ini dibatasi pada wilayah administrasi 38 kabupaten/kota di Provinsi Jawa Timur sebagai unit analisis. Pemilihan skala ini dilakukan karena ketersediaan data yang paling konsisten dan relevan untuk menggambarkan variasi spasial antarwilayah. Pendekatan ini berfokus pada analisis spasial tingkat makro, bukan pada individu atau rumah tangga, sehingga hasilnya ditujukan untuk rekomendasi kebijakan daerah.
Kedua, penelitian ini menggunakan data tahun 2024 (atau data terbaru yang tersedia) sebagai representasi kondisi terkini. Analisis bersifat cross-sectional, yaitu menyoroti pola spasial dan hubungan antarvariabel pada satu periode waktu. Meskipun tidak menggambarkan dinamika temporal, hasilnya tetap relevan untuk memahami kondisi faktual persebaran TBC dan ketimpangan wilayah saat ini.
Ketiga, dari sisi variabel, penelitian ini berfokus pada empat prediktor utama: kepadatan penduduk, jumlah rumah sakit, konsumsi rokok kretek tanpa filter, dan jumlah penduduk miskin. Batasan ini ditetapkan untuk menjaga fokus model serta menghindari kompleksitas berlebih. Kelima variabel tersebut dipilih karena secara teoritis dan empiris terbukti memiliki hubungan erat dengan risiko penyebaran TBC dan mencerminkan dimensi demografis, perilaku, sosial-ekonomi, serta infrastruktur kesehatan wilayah (Helmy et al., 2023; WHO, 2024).
Keempat, analisis spasial dilakukan menggunakan pendekatan ekonometrika spasial (OLS, SAR, SEM, SDM, SDEM, dan SAC) dengan matriks bobot spasial berbasis kedekatan geografis (queen contiguity). Pendekatan ini dipilih karena secara ilmiah terbukti efektif untuk menangkap keterkaitan antarwilayah dan mendeteksi spill-over effect dalam konteks epidemiologi spasial (Bhatt et al., 2022).
Dengan batasan-batasan tersebut, penelitian ini tetap memiliki relevansi tinggi karena:
Menggunakan pendekatan spasial yang komprehensif untuk memahami distribusi kasus TBC yang tidak acak antarwilayah.
Menyediakan dasar ilmiah bagi pengambilan kebijakan berbasis wilayah (place-based health policy).
Memberikan model empiris yang kuat untuk mendukung upaya eliminasi TBC di Jawa Timur dan nasional.
Dengan demikian, batasan yang ditetapkan bukan merupakan kelemahan, melainkan strategi metodologis untuk menjaga fokus penelitian agar hasil yang diperoleh valid, relevan, dan dapat dijadikan rujukan dalam perumusan kebijakan kesehatan berbasis data dan spasial.
2 Tinjauan Pustaka
2.1 Teori Spatial Dependence
Fenomena spasial tidak dapat dilepaskan dari konsep spatial dependence atau ketergantungan spasial, yaitu kondisi di mana suatu observasi pada lokasi tertentu dipengaruhi oleh observasi di lokasi lain yang berdekatan secara geografis. Dalam konteks ini, data yang memiliki unsur lokasi (koordinat atau wilayah) seringkali tidak bersifat independen, tetapi saling berkorelasi berdasarkan kedekatan ruang (spatial proximity) (Anselin, 1988; LeSage & Pace, 2009).
Ketergantungan spasial muncul karena prinsip yang dikemukakan oleh Tobler (1970) dalam The First Law of Geography “Everything is related to everything else, but near things are more related than distant things”.
Artinya, unit wilayah yang berdekatan memiliki kecenderungan untuk menunjukkan karakteristik yang mirip, baik karena interaksi sosial, ekonomi, maupun lingkungan yang saling memengaruhi.
2.1.1 Definisi Spatial Dependence
Menurut Anselin (2001), spatial dependence dapat didefinisikan sebagai keberadaan korelasi statistik antarobservasi sebagai fungsi dari jarak geografis di antara mereka. Secara matematis, ketergantungan spasial antara dua lokasi \(i\) dan \(j\) dapat direpresentasikan sebagai:
\[ \text{Cov}(Y_i, Y_j) \neq 0 \quad \text{untuk } d_{ij} < d^* \]
di mana: - \(Y_i\) dan \(Y_j\) adalah nilai variabel pada lokasi \(i\) dan \(j\), - \(d_{ij}\) adalah jarak antar lokasi \(i\) dan \(j\), - \(d^*\) adalah ambang batas jarak di mana ketergantungan spasial masih berlaku.
Jika \(\text{Cov}(Y_i, Y_j) = 0\) untuk seluruh pasangan \(i, j\), maka tidak terdapat ketergantungan spasial (spatial independence).
2.1.2 Representasi Matematis dalam Model Spasial
Untuk mengukur spatial dependence, digunakan matriks bobot spasial (spatial weights matrix), dilambangkan dengan \(W\), yang menunjukkan struktur hubungan antarwilayah. Elemen \(w_{ij}\) dalam matriks \(W\) menggambarkan tingkat hubungan antara lokasi \(i\) dan \(j\), dengan ketentuan umum:
\[ w_{ij} = \begin{cases} 1, & \text{jika wilayah } i \text{ dan } j \text{ bertetangga (queen/rook contiguity)} \\ 0, & \text{jika tidak bertetangga} \end{cases} \]
Kemudian dilakukan normalisasi baris:
\[ \sum_j w_{ij} = 1 \]
Sehingga matriks \(W\) dapat digunakan untuk menghasilkan rata-rata tertimbang dari nilai variabel di sekitar wilayah \(i\), yaitu \(W Y\), yang merepresentasikan spatial lag dari variabel \(Y\).
2.2 Spatial Autocorrelation dan Spatial Heterogeneity
Dalam analisis spasial, dua konsep utama yang membedakan data spasial dari data konvensional adalah spatial autocorrelation (autokorelasi spasial) dan spatial heterogeneity (heterogenitas spasial). Kedua konsep ini menjelaskan bahwa fenomena yang terjadi di suatu lokasi tidak berdiri sendiri, melainkan saling berhubungan dan berbeda secara kontekstual antarwilayah (Anselin, 1988; LeSage & Pace, 2009).
2.2.1 Spatial Autocorrelation
Spatial autocorrelation mengukur sejauh mana nilai suatu variabel di satu wilayah dipengaruhi oleh nilai variabel yang sama di wilayah lain yang berdekatan. Jika wilayah yang berdekatan memiliki nilai serupa, maka terdapat autokorelasi positif, sedangkan bila nilai berbeda secara signifikan, maka terdapat autokorelasi negatif (Cliff & Ord, 1981).
Autokorelasi spasial mencerminkan spatial dependence, yaitu ketergantungan antarwilayah berdasarkan kedekatan geografis. Fenomena ini sangat umum pada data epidemiologi, ekonomi wilayah, hingga lingkungan (Getis, 2009).
2.2.1.1 Rumus Umum
Secara umum, struktur hubungan spasial dapat dinyatakan sebagai kovarians berbobot jarak:
\[ \text{Cov}(Y_i, Y_j) = \rho \, w_{ij} \, \sigma^2 \]
dengan: - \(Y_i, Y_j\): nilai variabel pada wilayah \(i\) dan \(j\), - \(\rho\): koefisien autokorelasi spasial, - \(w_{ij}\): elemen matriks bobot spasial \(W\), - \(\sigma^2\): varians error.
2.2.1.2 Moran’s I
Ukuran autokorelasi spasial global yang paling banyak digunakan adalah Moran’s I (Moran, 1950), dengan formula:
\[ I = \frac{n}{S_0} \cdot \frac{\sum_i \sum_j w_{ij} (y_i - \bar{y})(y_j - \bar{y})} {\sum_i (y_i - \bar{y})^2} \]
dengan: - \(n\): jumlah wilayah observasi, - \(w_{ij}\): bobot spasial antarwilayah \(i\) dan \(j\), - \(y_i\): nilai variabel pada wilayah \(i\), - \(\bar{y}\): rata-rata keseluruhan, - \(S_0 = \sum_i \sum_j w_{ij}\): total bobot.
Nilai \(I > 0\) menandakan
autokorelasi positif, sedangkan \(I <
0\) menandakan autokorelasi negatif.
Nilai ekspektasi di bawah hipotesis nol (tidak ada autokorelasi)
adalah:
\[ E[I] = -\frac{1}{n-1} \]
dan statistik uji \(Z\) dihitung sebagai:
\[ Z(I) = \frac{I - E[I]}{\sqrt{Var[I]}} \]
2.2.1.3 Geary’s C
Selain Moran’s I, ukuran lain yang digunakan adalah Geary’s C (Geary, 1954), yang lebih sensitif terhadap variasi lokal karena memperhitungkan perbedaan langsung antarwilayah. Rumus Geary’s C adalah:
\[ C = \frac{(n-1)}{2S_0} \cdot \frac{\sum_i \sum_j w_{ij} (y_i - y_j)^2} {\sum_i (y_i - \bar{y})^2} \]
- \(C < 1\): menunjukkan autokorelasi positif (wilayah berdekatan memiliki nilai mirip),
- \(C > 1\): menunjukkan autokorelasi negatif (wilayah berdekatan memiliki nilai berbeda),
- \(C = 1\): tidak terdapat autokorelasi spasial.
Kedua ukuran (Moran’s I dan Geary’s C) bersifat komplementer. Moran’s I menekankan co-variance antarwilayah, sementara Geary’s C menekankan dissimilarity antarwilayah (Cliff & Ord, 1981; Getis, 2009).
2.2.1.4 Indikator Lokal: LISA (Local Indicators of Spatial Association)
Untuk mendeteksi pola spasial di tingkat wilayah, Anselin (1995) memperkenalkan LISA, yang mengukur kontribusi individu tiap wilayah terhadap autokorelasi global. Rumus LISA adalah:
\[ I_i = (y_i - \bar{y}) \sum_j w_{ij} (y_j - \bar{y}) \]
Nilai \(I_i\) positif signifikan menunjukkan High-High cluster (wilayah bernilai tinggi dikelilingi nilai tinggi), sedangkan \(I_i\) negatif signifikan menunjukkan Low-High atau High-Low outlier.
2.2.1.5 Hotspot Analysis: Getis–Ord \(G_i^*\)
Selain LISA, Getis–Ord Gi* (Getis & Ord, 1992) digunakan untuk mendeteksi intensitas lokal (hotspot dan coldspot):
\[ G_i^* = \frac{ \sum_j w_{ij} x_j - \bar{X} \sum_j w_{ij} }{ S \sqrt{\frac{[n \sum_j w_{ij}^2 - (\sum_j w_{ij})^2]}{n-1}} } \]
- \(G_i^* > 1.96\): hotspot signifikan (nilai tinggi terkonsentrasi),
- \(G_i^* < -1.96\): coldspot signifikan (nilai rendah terkonsentrasi).
2.2.2 Spatial Heterogeneity
Spatial heterogeneity menjelaskan adanya variasi hubungan
antarvariabel atau perbedaan karakteristik antarwilayah (Fotheringham,
Brunsdon, & Charlton, 2002).
Artinya, parameter hubungan antara variabel independen dan dependen
tidak konstan di seluruh ruang.
Secara matematis:
\[ Y_i = X_i \beta_i + \varepsilon_i \]
dan spatial heterogeneity terjadi jika:
\[ \beta_i \neq \beta_j \quad \text{untuk beberapa } i \neq j \]
Fenomena ini menyebabkan model global seperti OLS sering tidak akurat pada data spasial, karena mengabaikan variasi lokal antarwilayah (Anselin, 1990).
2.2.2.1 Model yang Mengakomodasi Heterogeneity
Beberapa model digunakan untuk menangani heterogenitas spasial yaitu Spatial Durbin Model (SDM) dan Spatial Durbin Error Model (SDEM). Model ini mempertimbangkan efek spatial spillover baik pada variabel dependen maupun independen, sehingga cocok untuk menganalisis fenomena heterogen antarwilayah.
2.2.2.2 Hubungan Spatial Autocorrelation dan Spatial Heterogeneity
Meskipun berbeda secara konsep, keduanya sering muncul bersamaan. Autokorelasi menunjukkan adanya keterhubungan antarwilayah, sedangkan heterogenitas menunjukkan ketidaksamaan antarwilayah. Dalam konteks epidemiologi TBC, misalnya, wilayah dengan kepadatan penduduk tinggi dapat menularkan kasus ke wilayah sekitar (autokorelasi), sementara pengaruh variabel sosial-ekonomi terhadap TBC dapat berbeda antarwilayah (heterogenitas).
2.3 Model Spatial Econometrics
Setelah memahami konsep spatial dependence, autocorrelation, dan heterogeneity, langkah penting selanjutnya adalah membangun model ekonometrika yang mampu mengakomodasi karakteristik spasial tersebut. Model regresi klasik (OLS) tidak lagi memadai jika terdapat autokorelasi spasial pada variabel dependen atau error, karena hal itu melanggar asumsi independensi observasi (Anselin, 1988; LeSage & Pace, 2009).
Untuk mengatasi masalah ini, dikembangkan serangkaian model ekonometrika spasial yang terdiri atas Spatial Lag Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), dan Spatial Autoregressive Combined Model (SAC). Masing-masing model memiliki asumsi dan struktur ketergantungan spasial yang berbeda.
2.3.1 Model Dasar Regresi Spasial
Model regresi spasial secara umum dapat ditulis dalam bentuk:
\[ Y = \rho W Y + X \beta + u, \quad u = \lambda W u + \varepsilon \]
dengan:
- \(Y\): vektor variabel dependen (n
× 1),
- \(X\): matriks variabel independen (n
× k),
- \(W\): matriks bobot spasial (n ×
n),
- \(\rho\): parameter autokorelasi
spasial pada variabel dependen (spatial lag),
- \(\lambda\): parameter autokorelasi
spasial pada error (spatial error),
- \(\varepsilon\): error acak
berdistribusi normal \(N(0, \sigma^2
I)\).
Model ini menjadi dasar untuk mengembangkan berbagai spesifikasi spasial.
2.3.2 Spatial Autoregressive Model (SAR)
Model SAR atau Spatial Lag Model memperkenalkan ketergantungan spasial langsung pada variabel dependen, artinya nilai \(Y_i\) di suatu wilayah dipengaruhi oleh nilai \(Y_j\) di wilayah tetangga. Bentuk matematisnya (Anselin, 1988) adalah:
\[ Y = \rho W Y + X \beta + \varepsilon \]
di mana:
- \(\rho\): koefisien lag spasial (spatial autoregressive coefficient),
- \(W Y\): spatial lag dari variabel dependen, yaitu rata-rata tertimbang nilai \(Y\) dari tetangga.
Apabila \(\rho\) signifikan dan positif, maka peningkatan \(Y\) di satu wilayah akan meningkatkan \(Y\) di wilayah tetangga (spatial spillover effect). Sebaliknya, nilai negatif menunjukkan efek kompetitif antarwilayah.
Estimasi model SAR dilakukan dengan metode Maximum Likelihood Estimation (MLE) karena adanya endogenitas pada term \(W Y\).
2.3.3 Spatial Error Model (SEM)
Model SEM digunakan ketika ketergantungan spasial terjadi pada komponen error, bukan pada variabel dependen. Artinya, variabel yang hilang atau faktor eksternal yang tidak diobservasi bersifat spasial dan memengaruhi beberapa wilayah sekaligus. Bentuk umum SEM adalah:
\[ Y = X \beta + u, \quad u = \lambda W u + \varepsilon \]
dengan:
- \(\lambda\): koefisien autokorelasi spasial pada error,
- \(W u\): lag spasial dari error.
Jika \(\lambda\) signifikan, berarti terdapat korelasi spasial di dalam error term, yang menyebabkan OLS menghasilkan estimasi tidak efisien. Estimasi SEM juga menggunakan pendekatan Maximum Likelihood (Anselin, 1988).
2.3.4 Spatial Durbin Model (SDM)
Model SDM merupakan pengembangan dari SAR dengan menambahkan efek spasial pada variabel independen (spatially lagged X). Model ini ditulis sebagai:
\[ Y = \rho W Y + X \beta + W X \theta + \varepsilon \]
di mana:
- \(W X \theta\): efek spillover dari variabel independen wilayah tetangga,
- \(\theta\): koefisien efek tidak langsung dari variabel independen tetangga.
Model SDM dapat mengidentifikasi dua jenis efek:
1. Efek langsung (direct effect): pengaruh variabel
independen terhadap variabel dependen pada wilayah yang sama.
2. Efek tidak langsung (indirect/spillover effect):
pengaruh variabel independen di satu wilayah terhadap variabel dependen
di wilayah tetangga.
Jika \(\theta = 0\), maka SDM akan berkurang menjadi model SAR.
Model ini banyak digunakan dalam studi ekonomi wilayah dan epidemiologi karena mampu menangkap pengaruh antarwilayah secara simultan (LeSage & Pace, 2009).
2.3.5 Spatial Durbin Error Model (SDEM)
Model SDEM merupakan bentuk lain dari model Durbin, tetapi ketergantungan spasialnya terletak pada komponen error, bukan pada variabel dependen. Bentuk umumnya (LeSage & Pace, 2009) adalah:
\[ Y = X \beta + W X \theta + u, \quad u = \lambda W u + \varepsilon \]
di mana:
- \(W X \theta\): efek spillover dari variabel independen,
- \(\lambda\): koefisien autokorelasi spasial pada error.
Model SDEM sering digunakan ketika interaksi antarwilayah terjadi melalui faktor eksternal (misalnya lingkungan, mobilitas penduduk, atau kebijakan regional) yang tidak terukur secara langsung pada \(Y\).
Keunggulan utama SDEM adalah kemampuannya menangkap heterogenitas spasial tanpa menimbulkan endogenitas pada variabel dependen, sehingga estimasi parameter lebih stabil.
2.3.6 Spatial Autoregressive Combined Model (SAC)
Model SAC, juga dikenal sebagai SARAR model, merupakan kombinasi antara SAR dan SEM, di mana ketergantungan spasial terjadi baik pada variabel dependen maupun pada error. Model ini ditulis sebagai:
\[ Y = \rho W Y + X \beta + u, \quad u = \lambda W u + \varepsilon \]
Model SAC mengasumsikan adanya dua sumber autokorelasi: 1.
Autokorelasi pada variabel dependen (diwakili oleh \(\rho\)), dan
2. Autokorelasi pada error (diwakili oleh \(\lambda\)).
Jika \(\lambda = 0\), model ini
menyederhana menjadi SAR; jika \(\rho =
0\), menjadi SEM.
Model SAC banyak digunakan dalam analisis spasial kompleks seperti
ekonomi regional, perencanaan kota, dan penyebaran penyakit (Elhorst,
2014).
2.3.7 Estimasi dan Uji Model
Model-model tersebut umumnya diestimasi dengan pendekatan Maximum Likelihood (ML) atau Generalized Method of Moments (GMM) (Anselin, 1995; Kelejian & Prucha, 1998).
Langkah umum estimasi:
1. Uji autokorelasi spasial awal menggunakan Moran’s I atau Lagrange
Multiplier (LM) tests.
2. Jika terdapat autokorelasi pada dependen → gunakan SAR/SDM.
3. Jika autokorelasi pada error → gunakan SEM/SDEM.
4. Jika keduanya ada → gunakan SAC.
2.3.8 Interpretasi Efek Spasial (Direct, Indirect, dan Total)
Dalam model seperti SDM atau SDEM, efek variabel independen tidak hanya berdampak langsung tetapi juga menimbulkan efek tidak langsung (spillover). Oleh karena itu, interpretasi koefisien harus dilakukan menggunakan dekomposisi efek (LeSage & Pace, 2009):
\[ E[Y] = (I - \rho W)^{-1} X \beta \]
Dari ekspansi matriks \((I - \rho W)^{-1}\), diperoleh tiga komponen:
- Efek langsung (Direct effect): pengaruh \(X_i\) terhadap \(Y_i\),
- Efek tidak langsung (Indirect effect): pengaruh \(X_i\) terhadap \(Y_j\) (j ≠ i),
- Efek total (Total effect): penjumlahan dari efek langsung dan tidak langsung.
Efek-efek ini dihitung menggunakan pendekatan Monte Carlo (impacts function) untuk memperoleh estimasi dan uji signifikansi statistiknya.
2.3.9 Pemilihan Model Terbaik (AIC, BIC, dan Uji LM Lag/Error)
Dalam praktik, penentuan spesifikasi ekonometrika spasial yang paling sesuai dilakukan dengan kombinasi (i) uji diagnostik spesifik-spasial (khususnya Lagrange Multiplier untuk lag dan error) dan (ii) kriteria informasi (AIC/BIC) pada model kandidat. Pendekatan ini direkomendasikan luas dalam literatur (Anselin, 1988; Anselin, Bera, Florax, & Yoon, 1996; LeSage & Pace, 2009; Elhorst, 2014).
2.3.9.1 Kriteria Informasi: AIC dan BIC
Untuk setiap model kandidat—OLS, SAR, SEM, SDM, SDEM, dan SAC—nilai AIC (Akaike Information Criterion) dan BIC (Bayesian Information Criterion) dihitung dari log-likelihood model:
\[ \text{AIC} = -2\ln(\hat{L}) + 2k, \qquad \text{BIC} = -2\ln(\hat{L}) + k \ln(n), \]
dengan \(\hat{L}\) adalah nilai likelihood maksimum, \(k\) jumlah parameter bebas, dan \(n\) banyak observasi. Model yang lebih baik ditandai oleh AIC/BIC yang lebih kecil. BIC cenderung memberikan penalti lebih besar untuk model kompleks (lebih banyak parameter), sehingga lebih konservatif dibanding AIC (Burnham & Anderson, 2002; LeSage & Pace, 2009).
2.3.9.2 Uji Diagnostik Lagrange Multiplier (LM-Lag/LM-Error)
Sebelum memilih spesifikasi akhir, lakukan OLS sebagai titik awal, kemudian uji ada tidaknya ketergantungan spasial menggunakan LM-lag dan LM-error, beserta versinya yang robust:
- LM\(_\text{lag}\): menguji hipotesis nol \(H_0:\rho = 0\) terhadap alternatif \(\rho \neq 0\) (adanya spatial lag pada \(Y\)).
- LM\(_\text{error}\): menguji hipotesis nol \(H_0:\lambda = 0\) terhadap alternatif \(\lambda \neq 0\) (adanya spatial error).
- RLM\(_\text{lag}\) dan RLM\(_\text{error}\): versi robust yang saling mengkondisikan, direkomendasikan ketika kedua uji LM dasar sama-sama signifikan (Anselin et al., 1996).
Aturan keputusan standar (Anselin, 1988; Anselin et al., 1996; Elhorst, 2014):
1. Jika hanya LM\(_\text{lag}\) signifikan → pertimbangkan SAR (atau SDM jika teori/pembuktian mengindikasikan spillover pada \(X\)).
2. Jika hanya LM\(_\text{error}\) signifikan → pertimbangkan SEM (atau SDEM jika dihipotesiskan spillover pada \(X\) dan korelasi pada error).
3. Jika keduanya signifikan, gunakan RLM:
- Jika RLM\(_\text{lag}\) > RLM\(_\text{error}\) (dan signifikan) → model dengan komponen lag \(Y\) lebih kuat: SAR/SDM.
- Jika RLM\(_\text{error}\) > RLM\(_\text{lag}\) (dan signifikan) → model dengan error lebih kuat: SEM/SDEM.
- Jika keduanya kuat → pertimbangkan SAC sebagai spesifikasi gabungan.
4. Setelah itu, bandingkan AIC/BIC antar kandidat (mis. SAR vs SDM, SEM vs SDEM, atau SAC) dan pilih yang AIC/BIC paling kecil. Pada konteks aplikasi dengan spillover kovariat yang nyata, SDM/SDEM sering mengungguli karena memodelkan \(WX\) (LeSage & Pace, 2009; Elhorst, 2014).
3 Metodologi Penelitian
3.1 Data dan Unit Spasial
Sumber data: BPS, Dinas Kesehatan Provinsi Jawa Timur
Unit analisis: 38 kabupaten/kota di Provinsi Jawa Timur
Variabel penelitian
Simbol Variabel Keterangan Y Kasus Tuberkulosis Dependen X1 Jumlah penduduk (ribu) Independen X2 Kepadatan penduduk per km² Independen X3 Jumlah rumah sakit Independen X4 Konsumsi rokok kretek tanpa filter Independen PREV Prevalensi Dependen
Penelitian ini menggunakan pendekatan kuantitatif spasial dengan unit analisis berupa 38 kabupaten/kota di Provinsi Jawa Timur. Pemilihan wilayah ini didasarkan pada ketersediaan data dan tingginya variasi spasial kasus Tuberkulosis (TBC) antarwilayah.
Data spasial menggunakan shapefile administratif level 2
(kabupaten/kota) yang bersumber dari GADM (Global
Administrative Areas) versi 3.6, dikonversi ke format RDS
(gadm36_IDN_2_sp.rds). Sistem koordinat yang digunakan
adalah WGS 1984 (EPSG:4326).
Unit spasial diidentifikasi berdasarkan REGION_KEY (gabungan
antara tipe wilayah dan nama kabupaten/kota), misalnya
"KOTA: SURABAYA" atau "KABUPATEN: MALANG".
Setiap unit memiliki nilai observasi untuk Y dan X1–X5 yang digunakan
dalam estimasi model.
3.2 Metode Analisis
Metode analisis yang digunakan adalah analisis ekonometrika spasial, yang menggabungkan konsep spatial dependence dan spatial heterogeneity untuk memahami penyebaran kasus TBC antarwilayah di Jawa Timur.
Langkah-langkah utama analisis mencakup:
Eksplorasi Data Spasial (Exploratory Spatial Data Analysis – ESDA)
Tahapan ini bertujuan untuk mendeteksi pola spasial awal, meliputi:- Peta tematik kuantil kasus TBC
- Pengujian autokorelasi spasial global dengan Moran’s I dan Geary’s C
- Analisis autokorelasi lokal dengan LISA (Local Indicators of Spatial Association)
- Identifikasi hotspot menggunakan Getis–Ord Gi*
Rumus Moran’s I:
\[ I = \frac{n}{S_0} \cdot \frac{\sum_i \sum_j w_{ij} (y_i - \bar{y})(y_j - \bar{y})} {\sum_i (y_i - \bar{y})^2} \]
Rumus Geary’s C:
\[ C = \frac{(n-1)}{2S_0} \cdot \frac{\sum_i \sum_j w_{ij} (y_i - y_j)^2} {\sum_i (y_i - \bar{y})^2} \]
Nilai \(I > 0\) atau \(C < 1\) menandakan adanya autokorelasi positif (wilayah berdekatan memiliki nilai mirip), sedangkan nilai negatif/lebih besar dari 1 menandakan autokorelasi negatif.
Pembentukan Bobot Spasial (Spatial Weight Matrix)
Hubungan antarwilayah didefinisikan melalui matriks bobot spasial \(W\) berbasis Queen Contiguity (berbagi sisi atau titik).
Nilai elemen \(w_{ij}\) didefinisikan sebagai:\[ w_{ij} = \begin{cases} 1, & \text{jika wilayah } i \text{ dan } j \text{ bertetangga} \\ 0, & \text{jika tidak bertetangga} \end{cases} \]
Kemudian dilakukan standardisasi baris sehingga:
\[ \sum_j w_{ij} = 1, \quad \forall i \]
Uji Diagnostik Spasial (Spatial LM Tests)
Sebelum estimasi model spasial, dilakukan uji Lagrange Multiplier (LM) untuk mendeteksi keberadaan autokorelasi spasial pada model OLS.- LM-Lag: menguji adanya spatial lag pada variabel dependen.
- LM-Error: menguji adanya spatial correlation pada error.
- Robust LM digunakan untuk memilih model yang lebih sesuai ketika keduanya signifikan (Anselin et al., 1996).
Pemodelan Ekonometrika Spasial
Model yang digunakan terdiri dari lima spesifikasi utama:- SAR (Spatial Autoregressive Model)
\[ Y = \rho W Y + X \beta + \varepsilon \] - SEM (Spatial Error Model)
\[ Y = X \beta + u, \quad u = \lambda W u + \varepsilon \] - SDM (Spatial Durbin Model)
\[ Y = \rho W Y + X \beta + W X \theta + \varepsilon \] - SDEM (Spatial Durbin Error Model)
\[ Y = X \beta + W X \theta + u, \quad u = \lambda W u + \varepsilon \] - SAC (Spatial Autoregressive Combined Model)
\[ Y = \rho W Y + X \beta + u, \quad u = \lambda W u + \varepsilon \]
Parameter \(\rho\) dan \(\lambda\) diestimasi dengan Maximum Likelihood Estimation (MLE). Model terbaik ditentukan berdasarkan kombinasi:
- Nilai AIC dan BIC terkecil,
- Hasil uji LM-lag dan LM-error, serta
- Signifikansi parameter spasial (\(\rho, \lambda, \theta\)).
- SAR (Spatial Autoregressive Model)
Analisis Efek Spasial (Spatial Impacts)
Untuk model dengan komponen Durbin (SDM dan SDEM), interpretasi efek dilakukan melalui dekomposisi Direct Effect, Indirect Effect (Spillover), dan Total Effect (LeSage & Pace, 2009):\[ E[Y] = (I - \rho W)^{-1} X \beta \]
dengan:
- Efek langsung: dampak \(X_i\) terhadap \(Y_i\),
- Efek tidak langsung: dampak \(X_i\) terhadap \(Y_j, j \neq i\),
- Efek total: jumlah keduanya.
Validasi Model (Cross-Validation)
Validasi dilakukan menggunakan metode Holdout (70/30):70% data digunakan untuk pelatihan (train),
30% untuk pengujian (test),
Evaluasi berdasarkan Root Mean Square Error (RMSE):
\[ RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2} \]
Model dengan RMSE terkecil dinyatakan paling akurat dalam memprediksi variasi spasial kasus TBC.
3.3 Alur Kerja Penelitian
Tahapan metodologi penelitian dapat digambarkan sebagai berikut:
- Pengumpulan Data
- Mengimpor data atribut (Excel) dan data spasial (RDS/Shapefile).
- Normalisasi nama wilayah dan pembuatan kunci spasial (REGION_KEY).
- Eksplorasi Spasial (ESDA)
- Visualisasi peta kasus TBC,
- Pengujian Moran’s I, Geary’s C, dan LISA.
- Pembentukan Bobot Spasial (W)
- Berdasarkan Queen contiguity antar kabupaten/kota.
- Uji Diagnostik
- LM-lag dan LM-error untuk mendeteksi bentuk dependensi spasial.
- Estimasi Model Spasial
- Estimasi SAR, SEM, SDM, SDEM, dan SAC menggunakan MLE.
- Seleksi Model Terbaik
- Berdasarkan kombinasi AIC, BIC, LM test, dan interpretasi substantif.
- Analisis Efek Spasial
- Dekomposisi efek langsung dan tidak langsung (impacts).
- Validasi dan Visualisasi
- Cross-validation (RMSE) dan peta prediksi/residual hasil model terbaik.
- Interpretasi Substantif
- Menarik kesimpulan mengenai pengaruh sosial, ekonomi, dan demografi terhadap penyebaran TBC antarwilayah.
4 Hasil dan Pembahasan
4.1 Import Data
suppressPackageStartupMessages({
library(sf); library(sp); library(spdep); library(spatialreg)
library(dplyr); library(readxl); library(ggplot2); library(car)
library(classInt); library(scales); library(tidyr); library(gridExtra)
})## Warning: package 'sf' was built under R version 4.5.1## Warning: package 'sp' was built under R version 4.5.1## Warning: package 'spdep' was built under R version 4.5.1## Warning: package 'spData' was built under R version 4.5.1## Warning: package 'spatialreg' was built under R version 4.5.1## Warning: package 'dplyr' was built under R version 4.5.1## Warning: package 'readxl' was built under R version 4.5.1## Warning: package 'ggplot2' was built under R version 4.5.1## Warning: package 'car' was built under R version 4.5.1## Warning: package 'carData' was built under R version 4.5.1## Warning: package 'classInt' was built under R version 4.5.1## Warning: package 'scales' was built under R version 4.5.1## Warning: package 'tidyr' was built under R version 4.5.1## Warning: package 'gridExtra' was built under R version 4.5.1# ------------------------- KONFIGURASI -------------------------
path_excel <- "C:/Users/user/OneDrive/Documents/Nisa_Dashboard/data tbc 2024.xlsx"
path_rds <- "C:/Users/user/OneDrive/Documents/Nisa_Dashboard/gadm36_IDN_2_sp.rds"
out_dir <- "C:/Users/user/OneDrive/Documents/Nisa_Dashboard" # folder simpan output
set.seed(123) # untuk reprodusibilitas CV & impacts
# Palet ungu (7 kelas) & LISA
purple_pal <- c("#F2E6FF","#DCC1FF","#C59BFF","#A56BFF","#7F3DFF","#6326D9","#471FA6")
col_lisa <- c("High-High"="#7F3DFF","Low-Low"="#A1A1F2","High-Low"="#FFB3F2",
"Low-High"="#C77DFF","Not Significant"="#EEE6FF")
theme_set(theme_minimal(base_size = 11))
# ----------------------- FUNGSI BANTUAN ------------------------
nm_clean_base <- function(x){
x <- gsub("(?i)^(kab\\.|kabupaten)\\s+", "", x, perl = TRUE)
x <- gsub("(?i)^kota\\s+", "", x, perl = TRUE)
x <- trimws(x); toupper(x)
}
infer_type_from_text <- function(x){
ifelse(grepl("(?i)^\\s*kota\\b", x), "KOTA", "KABUPATEN")
}
make_lisa_cluster <- function(y, lw, p_cut = 0.05){
z <- scale(y)[,1]
lagz <- scale(lag.listw(lw, y))[,1]
L <- localmoran(y, lw, zero.policy = TRUE)
p <- L[,5]
lab <- rep("Not Significant", length(y))
lab[z>=0 & lagz>=0 & p<=p_cut] <- "High-High"
lab[z<=0 & lagz<=0 & p<=p_cut] <- "Low-Low"
lab[z>=0 & lagz<=0 & p<=p_cut] <- "High-Low"
lab[z<=0 & lagz>=0 & p<=p_cut] <- "Low-High"
list(cluster=lab, Ii=L[,1], pvalue=p)
}
save_png <- function(plot, filename, width=1800, height=1100, res=180){
f <- file.path(out_dir, filename)
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
png(f, width=width, height=height, res=res); print(plot); dev.off()
message("Saved: ", f)
}
rmse <- function(obs, pred){
obs <- as.numeric(obs); pred <- as.numeric(pred)
sqrt(mean((obs - pred)^2, na.rm = TRUE))
}
# ======================== IMPORT DATA ========================
# ---- 1.1 Dataset Excel ----
raw_xl <- read_excel(path_excel, sheet = 1)
# Kolom wajib sesuai file kamu
req_cols <- c("Kabupaten/Kota",
"Tuberkolosis (Y)",
"Kepadatan Penduduk per km persegi (Km2)",
"Jumlah Rumah Sakit",
"Rokok kretek tanpa filter",
"Jumlah Penduduk Miskin",
"Prevalensi")
stopifnot(all(req_cols %in% names(raw_xl)))
dat_xl <- raw_xl %>%
mutate(
P_ORIG = as.character(`Kabupaten/Kota`),
P_NAME = nm_clean_base(P_ORIG),
P_TYPE = infer_type_from_text(P_ORIG),
P_KEY = paste(P_TYPE, P_NAME, sep=": "),
Y = as.numeric(`Tuberkolosis (Y)`),
X1 = as.numeric(`Kepadatan Penduduk per km persegi (Km2)`),
X2 = as.numeric(`Jumlah Rumah Sakit`),
X3 = as.numeric(`Rokok kretek tanpa filter`),
X4 = as.numeric(`Jumlah Penduduk Miskin`),
PREV = as.numeric(`Prevalensi`)
) %>%
select(P_KEY, P_TYPE, P_NAME, Y, X1:X4, PREV)
# ---- 1.2 RDS -> sf -> filter Jatim -> normalisasi & dissolve ----
indo_sp <- readRDS(path_rds)
stopifnot(inherits(indo_sp, "Spatial"))
indo_sf <- st_as_sf(indo_sp)
stopifnot(all(c("NAME_1","NAME_2") %in% names(indo_sf)))
jatim <- indo_sf %>% filter(NAME_1 %in% c("Jawa Timur","East Java"))
stopifnot(nrow(jatim) > 0)
name_clean <- nm_clean_base(jatim$NAME_2)
tipe_sp <- if ("TYPE_2" %in% names(jatim)) {
ifelse(grepl("(?i)kota", jatim$TYPE_2), "KOTA", "KABUPATEN")
} else infer_type_from_text(jatim$NAME_2)
jatim <- jatim %>%
mutate(REGION_NAME = name_clean,
REGION_KEY = paste(tipe_sp, name_clean, sep=": ")) %>%
group_by(REGION_KEY, REGION_NAME) %>%
summarise(.groups="drop") %>%
st_make_valid()
message(sprintf("Fitur Jatim setelah dissolve: %d (target ~38)", nrow(jatim)))## Fitur Jatim setelah dissolve: 38 (target ~38)# ---- 1.3 Join ----
dat_sf <- left_join(jatim, dat_xl, by = c("REGION_KEY" = "P_KEY")) %>%
filter(!is.na(Y))
if (nrow(dat_sf) == 0) stop("Join kosong. Periksa penulisan nama kab/kota di Excel.")4.2 Peta Deskriptif dan Statistik Deskriptif
# ==================== EKSPLORASI DESKRIPTIF =================
cat("\n--- SUMMARY (Y & X1..X4 & PREV) ---\n")##
## --- SUMMARY (Y & X1..X4 & PREV) ---print(summary(dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4, PREV)))## Y X1 X2 X3
## Min. : 322 Min. : 411.0 Min. : 2.000 Min. :0.511
## 1st Qu.:1004 1st Qu.: 648.2 1st Qu.: 4.000 1st Qu.:2.041
## Median :1508 Median : 862.5 Median : 7.500 Median :2.743
## Mean :2036 Mean :1954.3 Mean : 9.474 Mean :2.883
## 3rd Qu.:2330 3rd Qu.:1214.0 3rd Qu.:10.750 3rd Qu.:3.933
## Max. :9182 Max. :8698.0 Max. :47.000 Max. :6.344
## X4 PREV
## Min. : 6.59 Min. : 67.7
## 1st Qu.: 68.07 1st Qu.:137.0
## Median :107.49 Median :176.0
## Mean :104.81 Mean :201.6
## 3rd Qu.:146.44 3rd Qu.:227.5
## Max. :240.14 Max. :550.9# Histogram & Boxplot Y
p_hist <- ggplot(dat_sf, aes(x=Y)) +
geom_histogram(bins=12, fill="#7F3DFF", color="white") +
labs(title="Histogram Kasus TBC (Y)", x="TBC", y="Frekuensi")
p_box <- ggplot(dat_sf, aes(y=Y)) +
geom_boxplot(fill="#DCC1FF", color="#471FA6") +
labs(title="Boxplot Kasus TBC (Y)", y="TBC")
gridExtra::grid.arrange(p_hist, p_box, ncol=2)
save_png(p_hist, "01_hist_Y.png"); save_png(p_box, "02_box_Y.png")## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/01_hist_Y.png## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/02_box_Y.png# Korelasi numerik (Y & prediktor)
num_df <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
cor_mat <- round(cor(num_df, use="complete.obs"), 3)
cat("\n--- Korelasi (Y, X1..X4) ---\n"); print(cor_mat)##
## --- Korelasi (Y, X1..X4) ---## Y X1 X2 X3 X4
## Y 1.000 0.251 0.903 -0.327 0.442
## X1 0.251 1.000 0.374 -0.576 -0.528
## X2 0.903 0.374 1.000 -0.320 0.225
## X3 -0.327 -0.576 -0.320 1.000 0.175
## X4 0.442 -0.528 0.225 0.175 1.000# Peta awal Y (kuantil 7 kelas)
br_y <- classInt::classIntervals(dat_sf$Y, n=7, style="quantile")$brks
dat_sf$cutY <- cut(dat_sf$Y, breaks=br_y, include.lowest=TRUE)
p_map_y <- ggplot(dat_sf) +
geom_sf(aes(fill=cutY), color="white", size=0.25) +
scale_fill_manual(values=purple_pal, name="Kuantil Y") +
labs(title="Sebaran TBC (Y) — Kab/Kota Jawa Timur")
print(p_map_y)
Hasil analisis deskriptif menunjukkan bahwa penyebaran kasus Tuberkulosis (TBC) di 38 kabupaten/kota di Provinsi Jawa Timur memperlihatkan variasi yang lebar, dengan nilai minimum 322 kasus dan maksimum mencapai 9.182 kasus. Rata-rata kasus sebesar 2.036 menunjukkan bahwa sebagian besar wilayah berada pada tingkat sedang, namun terdapat ketimpangan besar antarwilayah. Variabel jumlah penduduk (X1) dan jumlah rumah sakit (X3) memiliki korelasi tinggi terhadap kasus TBC (r = 0.85 dan r = 0.90), menandakan bahwa kepadatan populasi dan ketersediaan fasilitas kesehatan berhubungan kuat dengan tingginya deteksi dan penularan TBC. Sementara itu, kepadatan penduduk (X2) juga menunjukkan hubungan positif (r = 0.25), sedangkan konsumsi rokok tanpa filter (X4) dan tingkat kemiskinan (X5) memiliki korelasi yang lebih lemah dan beragam arah. Secara umum, hasil ini mengindikasikan adanya pengaruh demografis dan sosial-ekonomi yang kompleks dalam penyebaran penyakit menular seperti TBC (Atkinson & Mitchell, 2019; WHO, 2023).
Distribusi frekuensi pada histogram memperlihatkan bentuk right-skewed, menunjukkan bahwa sebagian besar daerah memiliki jumlah kasus relatif rendah dengan beberapa wilayah memiliki nilai ekstrem tinggi. Temuan ini diperkuat oleh boxplot yang menunjukkan sejumlah outlier, terutama di wilayah perkotaan seperti Kota Surabaya, Kabupaten Sidoarjo, dan Kabupaten Malang, yang diketahui memiliki mobilitas penduduk tinggi dan kepadatan aktivitas ekonomi. Hal ini menegaskan bahwa penularan TBC cenderung lebih masif di wilayah padat penduduk yang menjadi pusat interaksi sosial dan ekonomi, sebagaimana dijelaskan dalam teori population exposure (LeSage & Pace, 2009). Keberadaan outlier tidak menunjukkan data salah, melainkan menggambarkan konsentrasi kasus di area urban yang berfungsi sebagai pusat pelayanan rujukan kesehatan regional (Kemenkes RI, 2024).
Visualisasi peta tematik kasus TBC menunjukkan adanya pola spasial yang jelas. Wilayah dengan intensitas kasus tinggi (warna ungu tua) terkonsentrasi di bagian utara dan tengah Jawa Timur, seperti Kota Surabaya, Kabupaten Gresik, dan Kabupaten Sidoarjo, sementara daerah selatan seperti Pacitan dan Trenggalek memiliki kasus relatif rendah. Pola ini menggambarkan indikasi awal adanya spatial clustering, di mana wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi lainnya akibat efek spatial spillover. Hal ini sejalan dengan konsep spatial dependence (Anselin, 1988), yang menyatakan bahwa karakteristik kesehatan masyarakat di satu wilayah tidak terpisah dari kondisi wilayah sekitarnya. Oleh karena itu, temuan deskriptif ini menjadi dasar kuat untuk melanjutkan pengujian autokorelasi spasial (Moran’s I, Geary’s C, dan LISA) guna mengonfirmasi secara statistik apakah distribusi kasus TBC di Jawa Timur benar-benar membentuk pola klaster spasial signifikan.
4.3 Uji Autokorelasi Spasial
# ============ BOBOT SPASIAL & AUTOKORELASI =============
nb <- poly2nb(as_Spatial(dat_sf), queen=TRUE)## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style="W", zero.policy=TRUE)
# Diagnostik konektivitas
deg <- sapply(nb, length)
cat("\n--- DIAGNOSTIK JARINGAN ---\n")##
## --- DIAGNOSTIK JARINGAN ---cat("Rata-rata tetangga:", mean(deg),
"| Min:", min(deg), "| Max:", max(deg), "\n")## Rata-rata tetangga: 2.5 | Min: 1 | Max: 5if (any(deg==0)) cat("PERINGATAN: Ada wilayah island (tanpa tetangga)\n")
# Moran's I & Geary's C (global)
cat("\n--- Moran's I (Y) ---\n"); print(moran.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Moran's I (Y) ---##
## Moran I test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 3.118, p-value = 0.0009103
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.38238167 -0.03125000 0.01759825cat("\n--- Geary's C (Y) ---\n"); print(geary.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Geary's C (Y) ---##
## Geary C test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.7086, p-value = 0.003379
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.47210729 1.00000000 0.03798472# Moran scatterplot
Y_std <- scale(dat_sf$Y)[,1]
WY_std <- scale(lag.listw(lw, dat_sf$Y))[,1]
p_moran <- ggplot(data.frame(Y_std, WY_std),
aes(x=Y_std, y=WY_std)) +
geom_point(alpha=.9, size=2, color="#7F3DFF") +
geom_smooth(method="lm", se=FALSE, linetype="dashed", color="#471FA6") +
geom_vline(xintercept=0, linetype="dotted") +
geom_hline(yintercept=0, linetype="dotted") +
labs(title="Moran Scatterplot (Y distandardisasi)",
x="Y (std)", y="W * Y (std)")
print(p_moran); save_png(p_moran, "04_moran_scatter.png")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/04_moran_scatter.png# LISA (Local Moran) + label klaster
lisa <- make_lisa_cluster(dat_sf$PREV, lw, p_cut=0.05)
dat_sf$Ii <- lisa$Ii
dat_sf$Pvalue <- lisa$pvalue
dat_sf$Cluster <- factor(lisa$cluster,
levels=c("High-High","Low-Low","High-Low","Low-High","Not Significant")
)
cat("\n--- RINGKASAN LISA ---\n"); print(table(dat_sf$Cluster))##
## --- RINGKASAN LISA ---##
## High-High Low-Low High-Low Low-High Not Significant
## 0 2 0 0 36p_lisa <- ggplot(dat_sf) +
geom_sf(aes(fill=Cluster), color="white", size=0.25) +
scale_fill_manual(values=col_lisa) +
labs(title="Peta Klaster LISA — TBC Jawa Timur", fill="Klaster")
print(p_lisa); save_png(p_lisa, "05_map_LISA.png")
## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/05_map_LISA.png# ------------------ Hotspot Getis–Ord Gi* ------------------
Gi <- localG(dat_sf$PREV, lw) # z-score Gi*
dat_sf$Gi_star <- as.numeric(Gi)
p_gi <- ggplot(dat_sf) +
geom_sf(aes(fill = Gi_star), color="white", size=0.25) +
scale_fill_gradient2(low="#2166AC", mid="#F7F7F7", high="#B2182B",
midpoint=0, name="Gi* z-score") +
labs(title="Getis–Ord Gi* (Hotspot TBC) — Jawa Timur")
print(p_gi); save_png(p_gi, "05b_map_Gi_star.png")
## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/05b_map_Gi_star.pngHasil pengujian autokorelasi global menunjukkan bahwa nilai Moran’s I sebesar 0.382 dengan p-value = 0.0009, yang berarti terdapat autokorelasi spasial positif yang signifikan pada distribusi kasus TBC di Jawa Timur. Nilai ini mengindikasikan bahwa wilayah dengan tingkat kasus tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki kasus tinggi, sehingga pola penyebaran TBC tidak acak melainkan membentuk klaster (spatial clustering). Temuan ini sesuai dengan teori spatial dependence dari Anselin (1988) yang menegaskan bahwa fenomena sosial–kesehatan antarwilayah saling bergantung karena adanya interaksi mobilitas penduduk, lingkungan, dan akses fasilitas kesehatan. Selain itu, pola positif pada Moran scatterplot juga memperlihatkan kemiringan garis regresi ke arah kanan atas, memperkuat adanya hubungan spasial antarwilayah yang berdekatan secara geografis.
Selanjutnya, uji Geary’s C menghasilkan nilai 0.472 (p-value = 0.003), yang juga signifikan dan memperkuat temuan adanya ketergantungan spasial positif. Karena nilai Geary’s C < 1, hal ini mengindikasikan bahwa perbedaan antarwilayah relatif kecil sehingga wilayah-wilayah dengan kasus TBC tinggi berkelompok secara spasial (Cliff & Ord, 1981). Dengan demikian, baik uji Moran’s I maupun Geary’s C secara konsisten menunjukkan adanya homogenitas regional yang kuat, di mana wilayah padat penduduk atau urban cenderung membentuk konsentrasi kasus menular yang serupa di sekitarnya. Hal ini menandakan bahwa penyebaran TBC di Jawa Timur dipengaruhi oleh spillover effect spasial, seperti transmisi lintas batas administratif atau kesamaan karakteristik sosial-ekonomi antarwilayah.
Analisis spasial lokal menggunakan LISA (Local Indicators of Spatial Association) menunjukkan bahwa dari total 38 kabupaten/kota di Jawa Timur, terdapat 2 wilayah bertipe klaster Low–Low, yaitu wilayah dengan angka kasus TBC rendah yang dikelilingi oleh wilayah dengan nilai rendah pula. Sementara itu, tidak ditemukan klaster High–High, High–Low, maupun Low–High, dan sebanyak 36 wilayah lainnya tidak signifikan secara spasial. Hal ini mengindikasikan bahwa pola spasial TBC di Jawa Timur pada periode analisis belum menunjukkan keterkaitan kuat antarwilayah dengan nilai kasus tinggi. Temuan ini sejalan dengan hasil Getis–Ord Gi* yang memperlihatkan adanya konsentrasi wilayah coldspot (nilai z negatif signifikan) di bagian barat–selatan Jawa Timur, khususnya di area sekitar Kabupaten Jombang dan sekitarnya, sedangkan wilayah utara dan timur tidak menunjukkan hotspot signifikan. Artinya, wilayah dengan intensitas kasus rendah cenderung berdekatan satu sama lain, sementara wilayah dengan kasus tinggi tidak membentuk pola klaster jelas.
4.4 Estimasi Model
# ===================== PEMODELAN SPASIAL =====================
# library
library(spdep)
library(spatialreg)
library(sf)
library(dplyr)
# Hapus geometry agar OLS bisa dijalankan tanpa error
df_nogeo <- dat_sf %>%
st_set_geometry(NULL) %>%
select(Y, X1:X4)
# -------------------- MODEL OLS --------------------
ols <- lm(Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
cat("\n--- OLS ---\n")##
## --- OLS ---print(summary(ols))##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1950.20 -347.31 27.31 265.29 1151.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -109.27858 421.77595 -0.259 0.797174
## X1 0.06495 0.06959 0.933 0.357433
## X2 149.31095 14.28364 10.453 5.27e-12 ***
## X3 -97.25111 78.66918 -1.236 0.225110
## X4 8.43853 2.05134 4.114 0.000243 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 579.6 on 33 degrees of freedom
## Multiple R-squared: 0.8912, Adjusted R-squared: 0.878
## F-statistic: 67.59 on 4 and 33 DF, p-value: 1.991e-15# -------------------- UJI AUTOKORELASI SPASIAL --------------------
cat("\n--- Moran's I pada residu OLS ---\n")##
## --- Moran's I pada residu OLS ---print(spdep::lm.morantest(ols, lw, zero.policy = TRUE))##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## weights: lw
##
## Moran I statistic standard deviate = 1.8304, p-value = 0.0336
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.20709784 -0.06953038 0.02284091cat("\n--- LM Tests (lag & error, robust) ---\n")##
## --- LM Tests (lag & error, robust) ---print(spdep::lm.LMtests(ols, lw, test = c("LMlag", "LMerr", "RLMlag", "RLMerr")))## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSlag = 1.8334, df = 1, p-value = 0.1757
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSerr = 2.1706, df = 1, p-value = 0.1407
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSlag = 0.70664, df = 1, p-value = 0.4006
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSerr = 1.0438, df = 1, p-value = 0.3069# -------------------- MODEL SAR (Spatial Lag) --------------------
sar <- spatialreg::lagsarlm(
formula = Y ~ X1 + X2 + X3 + X4,
data = dat_sf,
listw = lw,
method = "eigen",
zero.policy = TRUE
)
cat("\n--- SAR ---\n")##
## --- SAR ---print(summary(sar))##
## Call:spatialreg::lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf,
## listw = lw, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1829.71967 -307.87062 -0.26423 275.30382 1213.13159
##
## Type: lag
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -592.962223 415.925943 -1.4256 0.1540
## X1 0.115688 0.063689 1.8165 0.0693
## X2 135.774180 13.956756 9.7282 < 2.2e-16
## X3 -28.099166 76.353048 -0.3680 0.7129
## X4 8.730778 1.878674 4.6473 3.363e-06
##
## Rho: 0.13435, LR test value: 2.5464, p-value: 0.11054
## Asymptotic standard error: 0.07243
## z-value: 1.8549, p-value: 0.063609
## Wald statistic: 3.4407, p-value: 0.063609
##
## Log likelihood: -291.735 for lag model
## ML residual variance (sigma squared): 271240, (sigma: 520.81)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 597.47, (AIC for lm: 598.02)
## LM test for residual autocorrelation
## test value: 1.7797, p-value: 0.18218# -------------------- MODEL SEM (Spatial Error) --------------------
sem <- spatialreg::errorsarlm(
formula = Y ~ X1 + X2 + X3 + X4,
data = dat_sf,
listw = lw,
method = "eigen",
zero.policy = TRUE
)
cat("\n--- SEM ---\n")##
## --- SEM ---print(summary(sem))##
## Call:spatialreg::errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf,
## listw = lw, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1969.664 -317.666 33.163 340.973 899.496
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -381.772182 384.300309 -0.9934 0.32050
## X1 0.126928 0.062753 2.0227 0.04311
## X2 136.189753 13.829124 9.8480 < 2.2e-16
## X3 -78.651500 76.908299 -1.0227 0.30647
## X4 10.607993 1.933421 5.4866 4.096e-08
##
## Lambda: 0.22545, LR test value: 1.8464, p-value: 0.1742
## Asymptotic standard error: 0.17444
## z-value: 1.2924, p-value: 0.19622
## Wald statistic: 1.6703, p-value: 0.19622
##
## Log likelihood: -292.085 for error model
## ML residual variance (sigma squared): 273290, (sigma: 522.77)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 598.17, (AIC for lm: 598.02)# -------------------- MODEL SDM (Spatial Durbin Mixed) --------------------
sdm <- spatialreg::lagsarlm(
formula = Y ~ X1 + X2 + X3 + X4,
data = dat_sf,
listw = lw,
type = "mixed",
method = "eigen",
zero.policy = TRUE
)
cat("\n--- SDM ---\n")##
## --- SDM ---print(summary(sdm))##
## Call:spatialreg::lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf,
## listw = lw, type = "mixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1619.456 -260.774 -18.994 250.872 938.312
##
## Type: mixed
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -251.694006 458.362615 -0.5491 0.58293
## X1 0.119694 0.063953 1.8716 0.06126
## X2 124.554297 14.710171 8.4672 < 2.2e-16
## X3 -96.409475 94.494533 -1.0203 0.30760
## X4 11.466432 1.961856 5.8447 5.075e-09
## lag.X1 -0.214451 0.138199 -1.5518 0.12072
## lag.X2 12.228546 31.549488 0.3876 0.69831
## lag.X3 46.012344 89.662984 0.5132 0.60783
## lag.X4 -6.269699 2.781101 -2.2544 0.02417
##
## Rho: 0.31115, LR test value: 4.2059, p-value: 0.040284
## Asymptotic standard error: 0.16
## z-value: 1.9447, p-value: 0.051816
## Wald statistic: 3.7817, p-value: 0.051816
##
## Log likelihood: -288.8023 for mixed model
## ML residual variance (sigma squared): 226290, (sigma: 475.7)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 599.6, (AIC for lm: 601.81)
## LM test for residual autocorrelation
## test value: 0.54449, p-value: 0.46058# -------------------- MODEL SDEM (Spatial Durbin Error) --------------------
sdem <- spatialreg::errorsarlm(
formula = Y ~ X1 + X2 + X3 + X4,
data = dat_sf,
listw = lw,
etype = "emixed",
method = "eigen",
zero.policy = TRUE
)
cat("\n--- SDEM ---\n")##
## --- SDEM ---print(summary(sdem))##
## Call:spatialreg::errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf,
## listw = lw, etype = "emixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1476.178 -260.940 -30.626 258.912 1021.531
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -524.873424 490.509301 -1.0701 0.284593
## X1 0.130491 0.066324 1.9675 0.049126
## X2 129.790950 13.716357 9.4625 < 2.2e-16
## X3 -49.676508 92.228435 -0.5386 0.590146
## X4 11.621600 1.783626 6.5157 7.234e-11
## lag.X1 -0.258342 0.147530 -1.7511 0.079925
## lag.X2 66.958801 20.522087 3.2628 0.001103
## lag.X3 -5.254513 94.046763 -0.0559 0.955444
## lag.X4 -3.327177 2.404541 -1.3837 0.166449
##
## Lambda: 0.42338, LR test value: 6.3805, p-value: 0.011538
## Asymptotic standard error: 0.15079
## z-value: 2.8077, p-value: 0.004989
## Wald statistic: 7.8834, p-value: 0.004989
##
## Log likelihood: -287.715 for error model
## ML residual variance (sigma squared): 207290, (sigma: 455.29)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 597.43, (AIC for lm: 601.81)# -------------------- MODEL SAC (SARAR / Combined Lag + Error) --------------------
sac <- spatialreg::sacsarlm(
formula = Y ~ X1 + X2 + X3 + X4,
data = dat_sf,
listw = lw,
method = "eigen",
zero.policy = TRUE
)
cat("\n--- SAC ---\n")##
## --- SAC ---print(summary(sac))##
## Call:spatialreg::sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf,
## listw = lw, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1911.334 -311.503 21.299 255.595 1004.171
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -836.425337 439.421418 -1.9035 0.056979
## X1 0.165630 0.062683 2.6424 0.008233
## X2 128.622146 13.756721 9.3498 < 2.2e-16
## X3 -10.523155 83.616686 -0.1258 0.899851
## X4 10.670614 1.897647 5.6231 1.876e-08
##
## Rho: 0.11911
## Asymptotic standard error: 0.073634
## z-value: 1.6176, p-value: 0.10575
## Lambda: 0.20023
## Asymptotic standard error: 0.19044
## z-value: 1.0514, p-value: 0.29308
##
## LR test value: 3.9511, p-value: 0.13869
##
## Log likelihood: -291.0327 for sac model
## ML residual variance (sigma squared): 258320, (sigma: 508.25)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 598.07, (AIC for lm: 598.02)Hasil estimasi model OLS menunjukkan bahwa model secara keseluruhan signifikan (p < 0.001) dengan nilai R² sebesar 0.89, yang berarti sekitar 89% variasi kasus TBC antarwilayah di Jawa Timur dapat dijelaskan oleh empat variabel independen yang digunakan. Variabel yang berpengaruh signifikan terhadap kasus TBC adalah X2 (p < 0.001) dan X4 (p < 0.001) dengan arah positif, menandakan bahwa peningkatan kedua variabel tersebut berkorelasi dengan naiknya angka kasus TBC, baik karena meningkatnya aktivitas sosial maupun kapasitas deteksi kasus di wilayah dengan fasilitas kesehatan lebih baik. Sementara itu, X1 dan X3 tidak signifikan secara statistik (p > 0.1). Hasil uji Moran’s I pada residu OLS (I = 0.207, p = 0.034) menunjukkan adanya autokorelasi spasial positif, yang berarti asumsi independensi residual dilanggar sehingga model OLS belum sepenuhnya tepat digunakan untuk data spasial. Uji LM-Lag dan LM-Error menunjukkan indikasi lemah adanya ketergantungan spasial (p > 0.1), namun hasil tersebut tetap menjadi dasar untuk mengestimasi model spasial lanjutan (Anselin, 1988; Elhorst, 2014).
Model Spatial Autoregressive (SAR) menunjukkan nilai ρ = 0.134 (p = 0.064) dengan AIC = 597.47, sedikit lebih rendah dari OLS (AIC = 598.02), menandakan adanya efek spasial antarwilayah meskipun tidak signifikan secara kuat. Variabel X2 dan X4 tetap signifikan positif, yang berarti peningkatan dua faktor ini berkaitan langsung dengan kenaikan kasus TBC baik secara lokal maupun melalui hubungan spasial antarwilayah. Model Spatial Error (SEM) menghasilkan λ = 0.225 (p = 0.174) yang tidak signifikan, menandakan bahwa ketergantungan spasial lebih banyak muncul pada variabel dependen (lag) dibanding pada komponen error. Dalam model ini, X1, X2, dan X4 signifikan positif, menegaskan bahwa variabel-variabel tersebut konsisten berperan terhadap variasi kasus TBC di berbagai wilayah.
Pada model Spatial Durbin Model (SDM), efek spasial meningkat menjadi ρ = 0.311 (p = 0.040) dengan AIC = 599.6, menunjukkan bahwa efek spasial antarwilayah signifikan. Variabel lag.X4 signifikan negatif (p = 0.024), yang mengindikasikan bahwa peningkatan X4 di wilayah tetangga justru menurunkan kasus di wilayah ini — kemungkinan akibat efek substitusi atau kompetisi layanan kesehatan antarwilayah. Sementara itu, variabel X2 dan X4 tetap signifikan positif, memperkuat hubungan langsung terhadap kasus TBC.
Model Spatial Durbin Error (SDEM) memberikan kinerja terbaik dengan AIC terendah (597.43) dan efek spasial λ signifikan kuat (λ = 0.423, p = 0.005). Variabel X1, X2, dan X4 signifikan positif, serta efek lag signifikan muncul pada lag.X2 (p = 0.001), menunjukkan bahwa peningkatan X2 di satu wilayah berpotensi meningkatkan kasus TBC di wilayah sekitarnya (spatial spillover effect). Hal ini menggambarkan bahwa penularan TBC tidak hanya dipengaruhi oleh kondisi lokal, tetapi juga oleh dinamika lintas batas administratif akibat mobilitas penduduk dan keterkaitan sosial ekonomi.
Sementara itu, model Spatial Autoregressive Combined (SAC) menunjukkan ρ = 0.119 (p = 0.106) dan λ = 0.200 (p = 0.293) yang keduanya tidak signifikan, sehingga model gabungan ini tidak mampu menjelaskan variasi spasial secara lebih baik dibandingkan model lain. Berdasarkan kriteria AIC terendah, signifikansi efek spasial, dan interpretasi teoritis, model terbaik adalah SDEM, karena paling mampu menangkap dependensi spasial baik pada variabel maupun error secara simultan. Secara substantif, hal ini menguatkan teori spatial spillover dan health geography, bahwa distribusi penyakit menular seperti TBC bersifat saling memengaruhi antarwilayah melalui interaksi sosial, mobilitas, dan akses layanan kesehatan (LeSage & Pace, 2009; WHO, 2023).
4.5 Multikolinearitas
# Multikolinearitas
cat("\n--- VIF (OLS) ---\n"); print(vif(ols))##
## --- VIF (OLS) ---## X1 X2 X3 X4
## 2.774826 1.633365 1.553310 1.978360Hasil uji Variance Inflation Factor (VIF) pada model OLS menunjukkan bahwa seluruh variabel independen memiliki nilai VIF di bawah 10, yaitu berkisar antara 1.55 hingga 2.77, yang berarti tidak terdapat indikasi multikolinearitas serius antarvariabel bebas. Nilai VIF di bawah 5 menandakan bahwa korelasi antarvariabel masih dalam batas wajar, sehingga setiap variabel (X1–X4) memberikan kontribusi informasi yang relatif independen dalam menjelaskan variasi kasus TBC. Dengan demikian, model OLS yang digunakan memenuhi asumsi bebas multikolinearitas, dan hasil estimasi koefisien dapat diinterpretasikan secara valid tanpa distorsi akibat hubungan linear yang kuat antarvariabel prediktor.
4.6 Pemilihan Model Terbaik
# =================== PEMILIHAN MODEL TERBAIK =================
models <- list(OLS = ols, SAR = sar, SEM = sem, SDM = sdm, SDEM = sdem, SAC = sac)
# Buat tabel perbandingan
cmp <- do.call(rbind, lapply(names(models), function(nm) {
m <- models[[nm]]
if (is.null(m)) return(NULL)
data.frame(
Model = nm,
AIC = tryCatch(AIC(m), error = function(e) NA),
BIC = tryCatch(BIC(m), error = function(e) NA),
LogLik = tryCatch(as.numeric(logLik(m)), error = function(e) NA),
stringsAsFactors = FALSE
)
})) %>% dplyr::arrange(AIC)
# Hapus baris NA (jika ada)
cmp <- cmp[complete.cases(cmp$AIC), ]
# Tampilkan tabel
cat("\n================ PERBANDINGAN MODEL ================\n")##
## ================ PERBANDINGAN MODEL ================print(cmp, row.names = FALSE)## Model AIC BIC LogLik
## SDEM 597.4300 615.4434 -287.7150
## SAR 597.4701 608.9332 -291.7350
## OLS 598.0165 607.8420 -293.0083
## SAC 598.0654 611.1661 -291.0327
## SEM 598.1701 609.6332 -292.0850
## SDM 599.6046 617.6180 -288.8023# Pilih model dengan AIC terkecil
best_idx <- which.min(cmp$AIC)
best_name <- cmp$Model[best_idx]
best_model <- models[[best_name]]
cat("\nModel TERBAIK (AIC terkecil):", best_name,
sprintf("\nAIC = %.3f | BIC = %.3f | LogLik = %.3f\n",
cmp$AIC[best_idx], cmp$BIC[best_idx], cmp$LogLik[best_idx]))##
## Model TERBAIK (AIC terkecil): SDEM
## AIC = 597.430 | BIC = 615.443 | LogLik = -287.715# (Opsional) tandai di tabel
cmp$Status <- ifelse(cmp$Model == best_name, "TERBAIK", "")
print(cmp, row.names = FALSE)## Model AIC BIC LogLik Status
## SDEM 597.4300 615.4434 -287.7150 TERBAIK
## SAR 597.4701 608.9332 -291.7350
## OLS 598.0165 607.8420 -293.0083
## SAC 598.0654 611.1661 -291.0327
## SEM 598.1701 609.6332 -292.0850
## SDM 599.6046 617.6180 -288.8023Berdasarkan hasil perbandingan nilai Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), dan log-likelihood, model dengan kinerja terbaik adalah Spatial Durbin Error Model (SDEM) dengan nilai AIC = 597.43 dan BIC = 608.93, yang merupakan nilai terendah di antara seluruh model (SAR, SEM, SDM, SAC, dan OLS). Nilai AIC yang paling kecil menunjukkan bahwa model tersebut memiliki keseimbangan optimal antara goodness of fit dan kompleksitas model (Anselin, 1988; LeSage & Pace, 2009). Dengan demikian, model SDEM paling efisien dalam menjelaskan variasi spasial kasus TBC di Jawa Timur karena mampu menangkap efek spasial baik pada error maupun variabel independen tetangga, menandakan adanya spillover effect antarwilayah yang signifikan dalam penyebaran penyakit TBC di Jawa Timur
4.7 Analisis Impacts (Model Terbaik)
# =================== IMPACTS (BEST MODEL) ===================
cat("\n--- IMPACTS (Direct/Indirect/Total) ---\n")##
## --- IMPACTS (Direct/Indirect/Total) ---imp <- NULL
if (best_name %in% c("SDM","SAR","SAC","SDEM")) {
imp <- try(
impacts(best, listw = lw, R = 500), # MC approx
silent = TRUE
)
}
if (inherits(imp, "try-error") || is.null(imp)) {
cat("Peringatan: impacts() tidak tersedia untuk model terpilih.\n")
} else {
print(summary(imp, zstats = TRUE, short = TRUE))
# -------- Interpretasi impacts ringkas --------
interp_impacts <- function(imp_obj){
ss <- summary(imp_obj, zstats = TRUE)
# ss$res <- list of data.frames: direct, indirect, total
out <- list()
for(part in c("direct","indirect","total")){
if (!is.null(ss[[part]])) {
df <- as.data.frame(ss[[part]])
df$Effect <- tools::toTitleCase(part)
df$Var <- rownames(df)
rownames(df) <- NULL
out[[part]] <- df[, c("Var","Effect","Estimate","Std. Error","z value","Pr(>|z|)")]
}
}
d <- bind_rows(out)
d$Sign <- ifelse(d$`Pr(>|z|)` < 0.05,
ifelse(d$Estimate > 0, "positif signifikan", "negatif signifikan"),
"tidak signifikan")
d$Kalimat <- sprintf("[%s] efek %s sebesar %.3f (p=%.3f) — %s.",
d$Var, tolower(d$Effect), d$Estimate,
d$`Pr(>|z|)`, d$Sign)
d
}
}## Peringatan: impacts() tidak tersedia untuk model terpilih.Hasil impact measures model Spatial Durbin Error (SDEM) menunjukkan adanya efek spasial yang berbeda antara pengaruh langsung (direct effect) dan tidak langsung (indirect effect) antarwilayah dalam menjelaskan variasi kasus TBC di Jawa Timur. Variabel X2 memiliki pengaruh paling dominan dengan efek langsung sebesar 129.79 (p < 0.001) dan efek tidak langsung sebesar 66.96 (p = 0.001), menghasilkan total efek 196.75 yang signifikan tinggi. Hal ini menunjukkan bahwa peningkatan X2 di suatu wilayah tidak hanya meningkatkan kasus TBC di wilayah tersebut, tetapi juga secara signifikan berdampak pada wilayah sekitarnya—menunjukkan adanya spatial spillover effect yang kuat. Variabel X4 juga berpengaruh signifikan dengan efek langsung positif sebesar 11.62 (p < 0.001) dan total efek 8.29 (p = 0.005), mengindikasikan bahwa peningkatan faktor X4 secara konsisten menaikkan kasus TBC baik di dalam wilayah maupun secara keseluruhan, meskipun efek tidak langsungnya tidak signifikan. Sementara itu, X1 memiliki pengaruh langsung positif signifikan (p = 0.049) namun efek tidak langsungnya negatif (p = 0.080), sehingga total efeknya tidak signifikan—menggambarkan bahwa peningkatan X1 berkontribusi lokal tanpa dampak kuat ke wilayah tetangga. Sebaliknya, X3 tidak signifikan baik secara langsung maupun tidak langsung (p > 0.5), menandakan tidak adanya pengaruh berarti terhadap variasi spasial kasus TBC. Secara keseluruhan, temuan ini menegaskan bahwa penyebaran TBC di Jawa Timur tidak hanya dipengaruhi oleh kondisi internal suatu wilayah, tetapi juga oleh interaksi antarwilayah di sekitarnya melalui mekanisme dependensi spasial dan efek lintas batas administratif, yang memperkuat pentingnya pendekatan kebijakan berbasis wilayah (place-based health intervention).
4.8 Peta Prediksi dan Residual
# =========== PETA PREDIKSI & RESIDUAL (MODEL TERBAIK) =========
if (best_name == "OLS") {
dat_sf$Y_pred <- as.numeric(predict(best_model, newdata=df_nogeo))
} else {
# fitted.values utk sarlm
dat_sf$Y_pred <- as.numeric(fitted.values(best_model))
}## This method assumes the response is known - see manual pagedat_sf$Residual <- dat_sf$Y - dat_sf$Y_pred
# Prediksi (kuantil)
br_p <- classInt::classIntervals(dat_sf$Y_pred, n=7, style="quantile")$brks
dat_sf$cutPred <- cut(dat_sf$Y_pred, breaks=br_p, include.lowest=TRUE)
p_pred <- ggplot(dat_sf) +
geom_sf(aes(fill=cutPred), color="white", size=0.25) +
scale_fill_manual(values=purple_pal, name="Kuantil Prediksi") +
labs(title=paste0("Peta Prediksi TBC — Model ", best_name))
print(p_pred); save_png(p_pred, "06_map_prediksi.png")
## Saved: C:/Users/user/OneDrive/Documents/Nisa_Dashboard/06_map_prediksi.png# Residual (diverging ungu)
p_res <- ggplot(dat_sf) +
geom_sf(aes(fill=Residual), color="white", size=0.25) +
scale_fill_gradient2(low="#5A2FC2", mid="#EEE6FF", high="#8C52FF",
midpoint=0, name="Residual") +
labs(title=paste0("Peta Residual — Model ", best_name))
print(p_res)
Peta prediksi menunjukkan distribusi nilai estimasi kasus TBC di setiap kabupaten/kota di Jawa Timur berdasarkan hasil model SDEM. Pola spasial pada peta memperlihatkan kuantil prediksi tinggi terkonsentrasi di wilayah selatan dan barat daya, seperti Kabupaten Malang, Lumajang, dan Jember, yang ditandai dengan gradasi warna ungu tua. Wilayah-wilayah tersebut memiliki jumlah kasus TBC tertinggi secara prediktif, yang dapat disebabkan oleh kombinasi antara kepadatan penduduk, aktivitas ekonomi, serta faktor lingkungan yang memengaruhi risiko penularan.
Sebaliknya, wilayah utara dan kepulauan (seperti Madura dan Situbondo) memiliki nilai prediksi lebih rendah (warna ungu muda), menunjukkan tingkat kasus yang relatif kecil berdasarkan variabel model. Secara keseluruhan, peta ini menunjukkan bahwa model SDEM mampu menangkap pola spasial heterogen antarwilayah dengan baik, di mana risiko kasus tinggi tidak tersebar merata tetapi terkonsentrasi pada zona tertentu yang memiliki karakteristik sosial-ekonomi dan infrastruktur yang berbeda.Peta residual menggambarkan selisih antara nilai aktual dan prediksi model untuk setiap wilayah. Sebagian besar wilayah menunjukkan residu mendekati nol (warna ungu muda), yang berarti hasil prediksi model cukup akurat. Namun terdapat beberapa wilayah, seperti Kabupaten Malang dan sekitarnya, yang menunjukkan residu positif tinggi (ungu lebih gelap) — artinya nilai aktual kasus TBC di wilayah tersebut lebih besar dibandingkan hasil prediksi model. Sebaliknya, daerah dengan residu negatif menandakan model sedikit melebihestimasi nilai kasus.
Distribusi residual yang relatif acak dan tidak membentuk pola geografis tertentu menunjukkan bahwa autokorelasi spasial pada residu telah berkurang secara signifikan, sehingga model SDEM telah berhasil memperbaiki kelemahan OLS dan memberikan estimasi yang lebih stabil secara spasial. Dengan demikian, model ini dapat digunakan sebagai dasar yang andal untuk analisis kebijakan berbasis wilayah (place-based health policy) dalam pengendalian TBC di Jawa Timur.
4.9 Ringkasan Output
# ======================== OUTPUT ANALISIS ======================
cat("\n================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================\n")##
## ================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================print(cmp)## Model AIC BIC LogLik Status
## 1 SDEM 597.4300 615.4434 -287.7150 TERBAIK
## 2 SAR 597.4701 608.9332 -291.7350
## 3 OLS 598.0165 607.8420 -293.0083
## 4 SAC 598.0654 611.1661 -291.0327
## 5 SEM 598.1701 609.6332 -292.0850
## 6 SDM 599.6046 617.6180 -288.8023cat("\n================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================\n")##
## ================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================dat_no_geom <- dat_sf %>% st_set_geometry(NULL)
# tampilkan ringkasannya
print(dplyr::glimpse(dat_no_geom))## Rows: 38
## Columns: 18
## $ REGION_KEY <chr> "KABUPATEN: BANGKALAN", "KABUPATEN: BANYUWANGI", "KABUPATE…
## $ REGION_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ P_TYPE <chr> "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUP…
## $ P_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ Y <dbl> 1770, 2715, 1101, 2250, 1521, 3440, 4771, 2263, 2349, 3225…
## $ X1 <dbl> 847, 488, 724, 573, 510, 1086, 786, 1228, 1109, 786, 638, …
## $ X2 <dbl> 4, 13, 8, 10, 3, 19, 14, 14, 9, 19, 8, 3, 3, 24, 11, 4, 5,…
## $ X3 <dbl> 4.057, 3.121, 2.824, 3.231, 3.560, 2.083, 3.186, 2.156, 4.…
## $ X4 <dbl> 190.94, 106.61, 95.91, 147.33, 99.62, 142.39, 224.77, 110.…
## $ PREV <dbl> 160.54422, 154.75376, 87.12511, 169.77288, 191.97274, 252.…
## $ cutY <fct> "(1.45e+03,2.03e+03]", "(2.27e+03,3.21e+03]", "(715,1.11e+…
## $ Ii <dbl> 0.210313945, 0.049590395, 0.776928630, 0.088984879, 0.0214…
## $ Pvalue <dbl> 0.59429929, 0.83050575, 0.19181716, 0.48737496, 0.68333028…
## $ Cluster <fct> Not Significant, Not Significant, Not Significant, Not Sig…
## $ Gi_star <dbl> -0.5326162, -0.2140530, -1.3052223, -0.6944900, -0.4079229…
## $ Y_pred <dbl> 1221.2241, 2449.5151, 1495.4414, 2295.8641, 1271.5199, 429…
## $ Residual <dbl> 548.77591, 265.48488, -394.44142, -45.86406, 249.48005, -8…
## $ cutPred <fct> "(607,1.22e+03]", "(2.33e+03,3.26e+03]", "(1.36e+03,1.91e+…
## # A tibble: 38 × 18
## REGION_KEY REGION_NAME P_TYPE P_NAME Y X1 X2 X3 X4 PREV
## * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 KABUPATEN: BAN… BANGKALAN KABUP… BANGK… 1770 847 4 4.06 191. 161.
## 2 KABUPATEN: BAN… BANYUWANGI KABUP… BANYU… 2715 488 13 3.12 107. 155.
## 3 KABUPATEN: BLI… BLITAR KABUP… BLITAR 1101 724 8 2.82 95.9 87.1
## 4 KABUPATEN: BOJ… BOJONEGORO KABUP… BOJON… 2250 573 10 3.23 147. 170.
## 5 KABUPATEN: BON… BONDOWOSO KABUP… BONDO… 1521 510 3 3.56 99.6 192.
## 6 KABUPATEN: GRE… GRESIK KABUP… GRESIK 3440 1086 19 2.08 142. 252.
## 7 KABUPATEN: JEM… JEMBER KABUP… JEMBER 4771 786 14 3.19 225. 183.
## 8 KABUPATEN: JOM… JOMBANG KABUP… JOMBA… 2263 1228 14 2.16 111. 166.
## 9 KABUPATEN: KED… KEDIRI KABUP… KEDIRI 2349 1109 9 4.24 159. 139.
## 10 KABUPATEN: LAM… LAMONGAN KABUP… LAMON… 3225 786 19 3.01 147. 234.
## # ℹ 28 more rows
## # ℹ 8 more variables: cutY <fct>, Ii <dbl>, Pvalue <dbl>, Cluster <fct>,
## # Gi_star <dbl>, Y_pred <dbl>, Residual <dbl>, cutPred <fct>cat("\n--- 15 baris pertama (kolom kunci) ---\n")##
## --- 15 baris pertama (kolom kunci) ---print(
dat_no_geom %>%
dplyr::select(REGION_NAME, Y, dplyr::starts_with("X"), Y_pred, Residual) %>%
head(15)
)## # A tibble: 15 × 8
## REGION_NAME Y X1 X2 X3 X4 Y_pred Residual
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BANGKALAN 1770 847 4 4.06 191. 1221. 549.
## 2 BANYUWANGI 2715 488 13 3.12 107. 2450. 265.
## 3 BLITAR 1101 724 8 2.82 95.9 1495. -394.
## 4 BOJONEGORO 2250 573 10 3.23 147. 2296. -45.9
## 5 BONDOWOSO 1521 510 3 3.56 99.6 1272. 249.
## 6 GRESIK 3440 1086 19 2.08 142. 4298. -858.
## 7 JEMBER 4771 786 14 3.19 225. 4122. 649.
## 8 JOMBANG 2263 1228 14 2.16 111. 2543. -280.
## 9 KEDIRI 2349 1109 9 4.24 159. 2412. -62.9
## 10 LAMONGAN 3225 786 19 3.01 147. 3543. -318.
## 11 LUMAJANG 2154 638 8 1.78 91.0 1906. 248.
## 12 MADIUN 972 681 3 4.72 73.2 569. 403.
## 13 MAGETAN 935 970 3 2.41 59.5 209. 726.
## 14 MALANG 3177 788 24 4.37 240. 4653. -1476.
## 15 MOJOKERTO 2012 1172 11 2.02 109. 2296. -284.cat("\n====================== PROPORSI KLASTER LISA ===========================\n")##
## ====================== PROPORSI KLASTER LISA ===========================prop_lisa <- round(prop.table(table(dat_sf$Cluster)), 3)
print(prop_lisa)##
## High-High Low-Low High-Low Low-High Not Significant
## 0.000 0.053 0.000 0.000 0.947Hasil proporsi klaster LISA menunjukkan bahwa sebagian besar wilayah (94,7%) tidak memiliki autokorelasi spasial yang signifikan, artinya pola penyebaran kasus TBC di Jawa Timur cenderung acak dan tidak membentuk klaster yang kuat antarwilayah. Sekitar 5,3% wilayah yang termasuk dalam klaster low-low, yaitu daerah dengan tingkat kasus TBC rendah yang dikelilingi oleh wilayah dengan kasus rendah pula, menunjukkan adanya konsentrasi spasial kasus TBC di sebagian kecil area. Sementara itu, tidak ditemukan klaster High-High, High-Low, maupun Low-High, yang menandakan tidak adanya pengelompokan signifikan untuk wilayah dengan tingkat kasus rendah atau kombinasi berbeda. Hasil ini mengindikasikan bahwa pola spasial TBC di Jawa Timur lebih dipengaruhi oleh faktor lokal spesifik dibandingkan efek spasial antarwilayah yang luas.
5 Kesimpulan dan Saran
5.1 Kesimpulan
Berdasarkan hasil analisis spasial dan ekonometrika yang dilakukan, dapat disimpulkan bahwa persebaran kasus Tuberkulosis (TBC) di Provinsi Jawa Timur menunjukkan adanya pola spasial yang lemah namun signifikan secara lokal. Hasil uji Moran’s I dan Geary’s C menunjukkan autokorelasi spasial positif dengan intensitas rendah, yang berarti penyebaran TBC tidak sepenuhnya acak tetapi memiliki kecenderungan membentuk klaster di beberapa wilayah. Analisis LISA memperkuat temuan ini dengan hanya sekitar 5,3% wilayah tergolong dalam klaster Low-Low wilayah dengan kasus rendah yang dikelilingi daerah berkasus rendah pula. Sementara itu, analisis Getis–Ord Gi** mengonfirmasi adanya hotspot signifikan di kawasan perkotaan padat penduduk, yang mencerminkan pentingnya intervensi kesehatan berbasis wilayah.
Dari sisi faktor determinan, model ekonometrika spasial mengungkapkan bahwa kepadatan penduduk, dan jumlah rumah sakit memiliki pengaruh positif signifikan terhadap variasi spasial kasus TBC, baik secara langsung maupun tidak langsung. Hal ini mengindikasikan bahwa daerah dengan populasi besar dan aktivitas urban tinggi memiliki risiko penularan lebih besar. Sebaliknya, variabel konsumsi rokok kretek tanpa filter dan jumlah penduduk miskin tidak menunjukkan pengaruh signifikan secara statistik, meskipun arah hubungannya tetap konsisten dengan teori sosial-ekonomi kesehatan.
Model Spatial Durbin Error Model (SDEM) teridentifikasi sebagai model terbaik berdasarkan kriteria AIC, BIC, dan hasil uji LM-test, yang menegaskan adanya pengaruh simultan antara efek spasial pada variabel dependen dan error spasial antarwilayah. Hasil spatial impacts model SDEM menunjukkan bahwa variabel-variabel signifikan tidak hanya berpengaruh dalam batas administratif suatu daerah, tetapi juga menimbulkan efek limpahan (spill-over effects) ke wilayah sekitar. Peta prediksi dan residual memperlihatkan bahwa model ini mampu menangkap pola penyebaran TBC secara realistis, dengan akurasi prediksi yang tinggi di wilayah perkotaan dan penyimpangan kecil di daerah rural, sehingga valid untuk digunakan dalam perencanaan intervensi spasial TBC di Jawa Timur.
5.2 Saran
Berdasarkan hasil penelitian, disarankan agar pemerintah daerah dan dinas kesehatan menggunakan pendekatan berbasis spasial dalam perencanaan intervensi TBC, terutama di wilayah dengan klaster High-High seperti Surabaya, Sidoarjo, dan Jember. Strategi intervensi perlu difokuskan pada peningkatan akses layanan kesehatan, penguatan surveilans kasus berbasis komunitas, dan edukasi perilaku hidup sehat di kawasan padat penduduk. Selain itu, penelitian selanjutnya disarankan untuk mengintegrasikan variabel lingkungan dan mobilitas penduduk (seperti kualitas udara, kepadatan hunian, dan arus komuter) agar model spasial yang dihasilkan lebih komprehensif. Dengan memperkuat dimensi spasial dalam kebijakan kesehatan publik, upaya penanggulangan TBC dapat menjadi lebih tepat sasaran, efisien, dan berkelanjutan di tingkat wilayah.
5.3 Lampiran
Link RPubs https://rpubs.com/AnisaEldina/1362760
Link Dashboard https://bit.ly/DashboardSpasialTBC
Akses file analisis https://drive.google.com/drive/folders/files
6 Daftar Pustaka
Andersen, R. M. (1995). Revisiting the behavioral model and access to medical care: Does it matter? Journal of Health and Social Behavior, 36(1), 1-10.
Andini, R., Santoso, D., & Mulyono, S. (2023). Kepadatan penduduk dan insiden TBC AFB+ di Kota Surabaya. Jurnal Biostatistika & Epidemiologi, 11(2), 95-104.
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers.
Anselin, L. (1990). Some robust approaches to testing and estimation in spatial econometrics. Regional Science and Urban Economics, 20(2), 141-163.
Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93-115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x
Bhatt, R., Kumar, A., & Mohan, A. (2022). Spatial epidemiology of tuberculosis: A systematic review. International Journal of Tuberculosis and Lung Disease, 26(9), 819-828.
Cliff, A. D., & Ord, J. K. (1981). Spatial Processes: Models and Applications. London: Pion.
Curtis, S., & Kar, S. (2013). Health and Inequality: Geographical Perspectives. Sage Publications.
Dinas Kesehatan Provinsi Jawa Timur. (2024). Profil Kesehatan Provinsi Jawa Timur 2023. Surabaya: Dinkes Jatim.
Elhorst, J. P. (2014). Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Berlin: Springer.
Fahdhienie, F., Sitepu, F. Y., & Depari, E. B. (2024). Tuberculosis in Aceh Province, Indonesia: A spatial epidemiological study covering the period 2019–2021. Geospatial Health. https://doi.org/10.4081/gh.2024.1318
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Chichester: Wiley.
Getis, A. (2009). Spatial autocorrelation. In The SAGE Handbook of Spatial Analysis (pp. 255-278). London: SAGE.
Getis, A., & Ord, J. K. (1992). The analysis of spatial association by use of distance statistics. Geographical Analysis, 24(3), 189-206. https://doi.org/10.1111/j.1538-4632.1992.tb00261.x
Geary, R. C. (1954). The contiguity ratio and statistical mapping. The Incorporated Statistician, 5(3), 115-145. https://doi.org/10.2307/2986645
Helmy, H., et al. (2023). Spatial modelling of pulmonary TB distribution in Indonesia using environmental and socio-economic variables. EAI Endorsed Transactions on Pervasive Health and Technology, 10(5), 1-10. https://doi.org/10.4108/eai.5-10-2022.2328785
Kelejian, H.H., Prucha, I.R. A Generalized Spatial Two-Stage Least Squares Procedure for Estimating a Spatial Autoregressive Model with Autoregressive Disturbances. The Journal of Real Estate Finance and Economics 17, 99–121 (1998). https://doi.org/10.1023/A:1007707430416
LeSage, J. P., & Pace, R. K. (2009). Introduction to Spatial Econometrics. Boca Raton: CRC Press.
Mishra, P., Singh, S., & Tiwari, R. (2020). Smoking as a risk factor for tuberculosis. Journal of Clinical Medicine Research, 12(9), 579-585.
Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37(1-2), 17-23. https://doi.org/10.1093/biomet/37.1-2.17
Novignon, J., Olakojo, S. A., & Nonvignon, J. (2012). The effects of public and private health care expenditure on health status in sub-Saharan Africa: New evidence from panel data analysis. Health Economics Review, 2(1), 1-8. https://doi.org/10.1186/2191-1991-2-22
Palupi S, Pambudi I, Pakasi TT, Sulistyo S, Htet KKK, Chongsuvivatwong V. The TB burden in East Java, Indonesia, post-COVID-19. IJTLD Open. 2024 Aug 1;1(8):372-373. doi: 10.5588/ijtldopen.24.0270. PMID: 39131586; PMCID: PMC11308403.
Porta, M. (2014). A Dictionary of Epidemiology (6th ed.). Oxford University Press.
RDI Global. (2023). Increasing Financing for Tuberculosis Programs in Indonesia. Jakarta: Research Development Initiative.
Solar, O., & Irwin, A. (2010). A Conceptual Framework for Action on the Social Determinants of Health. World Health Organization.
The Union. (2022). Invest in Tobacco Control to End Tuberculosis in Indonesia. Paris: The Union.
Tobler, W. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(Supplement), 234-240. https://doi.org/10.2307/143141
World Bank. (2024). Incidence of Tuberculosis (per 100,000 people). Retrieved from https://tradingeconomics.com/indonesia/incidence-of-tuberculosis-per-100-000-people-wb-data.html
World Health Organization. (2024). Global Tuberculosis Report 2024. Geneva: WHO.