ANALISIS SPASIAL KASUS TUBERKULOSIS DI SUMATERA UTARA PADA 2024
Ditujukan Untuk Memenuhi Nilai Mata Kuliah Spatial Statistics TA. 2025/2026
ANALISIS SPASIAL KASUS TUBERKULOSIS DI SUMATRA UTARA TAHUN 2024
Disusun oleh:
Charles Joshua Nathaniel Waruwu (140610230048)
Dosen Pengampu:
Dr. I Gede Nyoman Mindra Jaya, M.Si
PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025
Pendahuluan
Latar Belakang
Tuberkulosis (TBC) masih menjadi masalah kesehatan masyarakat yang signifikan di berbagai negara, termasuk Indonesia. Selain meningkatnya morbiditas dan mortalitas, penyakit ini menimbulkan tantangan sosial-ekonomi bagi individu dan keluarga yang terdampak. Upaya pengendalian TBC mengarahkan pemahaman yang mendalam mengenai pola sebarannya sehingga intervensi kesehatan dapat diarahkan secara tepat sasaran. Pendekatan spasial memberikan kerangka kerja untuk mengidentifikasi daerah-daerah dengan risiko tinggi dan mendeteksi klaster penyakit yang mungkin luput dari analisis agregat konvensional.
Sebagai salah satu provinsi terbesar di Indonesia, Sumatera Utara memiliki keragaman demografis, infrastruktur kesehatan, serta sanitasi yang berbeda di setiap kabupaten dan kota. Perbedaan tersebut dapat memengaruhi penularan serta keberhasilan program pengobatan TBC. Oleh karena itu, analisis yang memanfaatkan data spasial dan variabel kontekstual seperti angka kesembuhan, akses air minum layak, sanitasi layak, kepadatan penduduk, dan jumlah puskesmas, diperlukan untuk menggambarkan hubungan antar faktor dan distribusi kasus TBC secara geografis pada tahun 2024. Hal ini menjadi dasar bagi pengambilan keputusan program kesehatan di tingkat kabupaten/kota.
Beberapa penelitian menyebutkan bahwa distribusi TBC tidak bersifat acak secara geografis. Penelitian oleh Anselin (1995) memperkenalkan konsep Local Indicators of Spatial Association (LISA) yang banyak digunakan untuk mengidentifikasi klaster penyakit menular, termasuk TBC. Penelitian lanjutan di Indonesia menunjukkan adanya ketimpangan beban TBC antar provinsi maupun antar kabupaten/kota. Hal ini mengindikasikan bahwa analisis statistik non-spasial saja kurang memadai untuk memahami dinamika penularan TBC secara komprehensif.
Analisis spasial menyediakan alat dan metode yang mampu mengidentifikasi konsentrasi penyakit, mendeteksi klaster statistik, serta menguji keterkaitan antara kondisi lingkungan atau layanan kesehatan dengan penyakit TBC. Hasil analisis spasial yang jelas dan terukur memiliki implikasi kebijakan yang konkret. Informasi spasial mendukung evaluasi efektivitas program intervensi dari waktu ke waktu dan membantu arahan upaya pendidikan serta kampanye kesehatan masyarakat pada kelompok dan wilayah yang paling berisiko.
Dengan demikian, penelitian analisis spasial mengenai Tuberkulosis ini diharapkan dapat menghasilkan peta risiko serta bukti empiris yang berguna bagi pengambilan kebijakan kesehatan yang lebih efektif. Penelitian ini juga membuka peluang untuk pengembangan model prediksi dan evaluasi dampak intervensi berbasis lokal, sehingga upaya pengendalian TBC dapat lebih terukur dan berorientasi hasil.
Identifikasi Masalah
Berdasarkan latar belakang tersebut, dapat dilakukan identifikasi masalah dengan rincian sebagai berikut:
TBC masih menjadi masalah kesehatan masyarakat yang signifikan di Indonesia, yang berdampak pada sosial-ekonomi individu dan keluarga.
Distribusi kasus TBC antar wilayah belum sepenuhnya merata, khususnya di Provinsi Sumatera Utara, yang memiliki perbedaan karakteristik demografis, sanitasi, serta infrastruktur kesehatan.
Pendekatan analisis agregat konvensional belum mampu menggambarkan pola sebaran TBC secara rinci.
Belum optimalnya pemanfaatan informasi spasial dalam perencanaan dan evaluasi program pengendalian TBC.
Kebutuhan akan model analisis yang mampu mengintegrasikan data epidemiologu dengan aspek spasial sehingga dinamika penularan TBC dapat dipahami secara mendalam.
Tujuan Penelitian
Tujuan umum dari penelitian ini adalah untuk menganalisis pola sebaran spasial Tuberkulosis (TBC) di Provinsi Sumatera Utara tahun 2024 serta mengidentifikasi wilayah-wilayah yang memiliki risiko tinggi kejadian TBC sebagai dasar perencaan kebijakan yang lebih tepat sasaran. Adapun tujuan khususnya adalah sebagai berikut:
Menggambarkan pola distribusi geografis kasus Tuberkulosis (TBC) di tingkat kabupaten/kota di Provinsi Sumatera Utara tahun 2024.
Mengidentifikasi adanya pengelompokan (klaster) spasial kasus TBC di Provinsi Sumatera Utara menggunakan pendekatan analisis spasial.
Menyusun peta risiko TBC berbasis wilayah sebagai alat bantu visual untuk mendukung pengambilan keputusan dalam program pengendalian TBC.
Memberikan informasi empiris berbasis spasial yang dapat dimanfaatkan oleh pemangku kebijakan dalam penentuan wilayah prioritas intervensi kesehatan.
Batasan Penelitian
Agar penelitian ini terarah dan sesuai dengan tujuan yang ditetapkan, maka batasan penelitian ditentukan sebagai berikut:
Ruang lingkup wilayah dalam penelitian ini dibatasi pada seluruh kabupaten dan kota di Provinsi Sumatera Utara sesuai dengan batas administrasi wilayah yang digunakan dalam data spasial.
Ruang lingkup waktu penelitian dibatasi pada tahun 2024, sehingga analisis yang dilakukan bersifat potret kondisi (cross-sectional) dan tidak membahas perubahan pola sebaran TBC antar tahun.
Objek penelitian dalam studi ini adalah Tuberkulosis (TBC) yang dianalisis berdasarkan data agregat pada tingkat kabupaten/kota, tanpa membedakan jenis TBC (paru atau ekstra paru) maupun karakteristik individu penderita.
Unit analisis spasial dibatasi pada batas administrasi kabupaten/kota, sehingga variasi kejadian TBC pada tingkat yang lebih kecil (misalnya kecamatan atau desa) tidak dibahas dalam penelitian ini.
Model regresi spasial yang digunakan dibatasi pada SAR, SEM, SDM, SAC, serta OLS untuk memudahkan dalam pembelajaran. Masih terdapat beberapa pendekatan model lain yang relevan, namun tidak dilakukan dalam penelitian ini.
Tinjauan Pustaka
Spatial Dependence
Dependensi Spasial merupakan kondisi di mana nilai suatu variabel pada suatu wilayah dipengaruhi oleh nilai variabel yang sama di wilayah sekitaranya, sehingga observasi antarwilayah tidak bersifat independen. Hal ini muncul akibat pendekatan geografis, interaksi sosial, atau kesamaan karakteristik lingkungan yang menyebabkan fenomena di suatu wilayah dapat menyebar ke lokasi lain di sekitarnya. Keberadaan ketergantungan spasial menyebabkan asumsi independensi antarwilayah tidak terpenuhi, sehingga analisis statistik perlu mempertimbangkan aspek lokasi.
Autokorelasi Spasial
Autokorelasi spasial adalah estimasi dari hubungan antara nilai observasi yang berkaitan dengan lokasi geografis pada variabel yang sama. Autokorelasi spasial positif terjadi ketika variabel yang diobservasi pada lokasi yang berdekatan cenderung mirip atau berkorelasi positif (+1). Sementara itu, autokorelasi spasial negatif terjadi ketika variabel pada lokasi yang berdekatan cenderung berbeda atau berkorelasi negatif (-1).
Secara umum, autokorelasi spasial memiliki ukuran global klasik seperti Moran’s I dan Geary’s C, begitu juga analisis lokal menggunakan Local Moran’s I dan Getis-Ord.
1. Moran’s I (Global)
Moran’s I merupakan ukuran global autokorelasi spasial yang memiliki kemiripan dengan koefisien korelasi Pearson, tetapi disesuaikan untuk data spasial. Statistik ini mengukur sejauh mana nilai-nilai yang serupa berkumpul secara spasial. Nilai yang berkisar antara -1 (autokorelasi negatif sempurna) hingga +1 (autokorelasi positif sempurna), dengan nilai mendekati 0 menunjukkan tidak ada autokorelasi.
Secara matematis, dapat dinyatakan sebagai:
\[ I = \frac{N}{S_0} \cdot \frac{ \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x})(x_j - \bar{x}) }{ \sum_{i=1}^N (x_i - \bar{x})^2 } \]
di mana \(N\) adalah jumlah observasi, \(w_{ij}\) adalah bobot spasial antara lokasi \(i\) dan \(j\), \(x_i\) adalah nilai variabel pada lokasi \(i\), dan \(\bar{x}\) adalah rata-rata nilai variabel.
2. Geary’s C (Global)
Geary’s C adalah ukuran global autokorelasi spasial yang berfokus pada perbedaan kuadrat antara nilai-nilai berdekatan. Berbeda dengan Moran’s I, Geary’s C lebih sensitif terhadap perbedaan lokal. Nilai mendekati 0 menunjukkan autokorelasi positif kuat, sekitar 1 menunjukkan tidak ada autokorelasi, dan di atas 1 menunjukkan autokorelasi negatif.
Secara matematis, dapat dinyatakan sebagai:
\[ C = \frac{(N - 1) \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - x_j)^2}{2 S_0 \sum_{i=1}^N (x_i - \bar{x})^2} \]
di mana \(S_0 = \sum_{i=1}^N \sum_{j=1}^N w_{ij}\), dan variabel lain sama seperti sebelumnya.
- Interpretasi: \(C < 1\) menunjukkan autokorelasi positif (tetangga mirip), \(C > 1\) menunjukkan autokorelasi negatid (tetangga berlawanan), serta \(C = 1\) menunjukkan tidak ada autokorelasi (acak).
3. Local Moran’s I_i (LISA)
Local Moran’s I_i adalah versi lokal dari Moran’s I yang digunakan untuk mengidentifikasi cluster atau pencilan spesifik pada tiap lokasi. Hal ini membantu mendeteksi area high-high (cluster nilai tinggi), low-low (cluster nilai rendah), high-low (outlier tinggi di tengah rendah), atau low-high (outlier rendah di tengah tinggi).
Secara matematis, dapat dinyatakan sebagai:
\[ I_i = \frac{x_i - \bar{x}}{m_2} \sum_{j=1}^N w_{ij} (x_j - \bar{x}) \]
Di mana \(m_2 = \frac{1}{N} \sum_{k=1}^N (x_k - \bar{x})^2 , (x_i)\) adalah nilai pada lokasi \(i\), \(\bar{x}\) adalah rata-rata, \(w_{ij}\) adalah bobot spasial, dan \(N\) adalah jumlah observasi.
4. Getis-Ord G dan G*
Getis-Ord Gi dan Gi^* adalah statistik lokal untuk mendeteksi hot spots (area dengan nilai tinggi dikelilingi nilai tinggi) dan cold spots (area dengan nilai rendah dikelilingi nilai rendah). \(G_i\) tidak menyertakan nilai lokasi itu sendiri, sementara \(G_i^*\) menyertakannya untuk pertimbangan lebih inklusif. Untuk suatu band jarak \(d\) (atau bobot ketetaanggaan biner) dengan \(w_{ij}(d) \in \{0, 1\}\):
- \(G_i\) (tanpa nilai sendiri):
\[ G_i(d) = \frac{\sum_{j \neq i} w_{ij}(d) x_j}{\sum_{j \neq i} x_j} \]
Di mana \(w_{ij}(d)\) adalah bobot spasial berdasarkan jarak \(d\), \(x_j\) adalah nilai variabel pada lokasi \(j\), dan penyebut adalah jumlah nilai kecuali lokasi \(i\).
- Contoh: Di data curah hujan Jakarta pada September 2025, \(G_i(d) = 0.7\) di kecamatan Tambora menunjukkan cluster hujan tinggi di sekitar, tapi tidak termasuk data Tambora sendiri.
- \(G_i^*\) (termasuk nilai sendiri):
\[ G_i^*(d) = \frac{\sum_{j=1}^N w_{ij}(d) x_j}{\sum_{j=1}^N x_j}, \quad \text{dengan } w_{ii}(d) = 1 \]
Di mana \(w_{ij}(d)\) adalah bobot spasial, termasuk \(w_{ii}(d) = 1\) untuk lokasi itu sendiri, dan \(x_j\) adalah nilai variabel.
Model Spasial Ekonometrik
Model spasial ekonometrik merupakan pendekatan analisis yang digunakan untuk mengakomodasi adanya ketergantungan spasial antarwilayah dalam suatu model statistik atau ekonometrik. Model ini dikembangkan untuk mengatasi pelanggaran asumsi independensi observasi yang sering terjadi pada data berbasis wilayah, di mana nilai suatu variabel di satu lokasi dapat dipengaruhi oleh lokasi lain di sekitarnya. Terdapat beberapa model spasial yang akan digunakan, seperti:
Ordinary Least Squares (OLS)
Ordinary Least Squares (OLS) merupakan metode regresi linier klasik yang digunakan untuk mengestimasi hubungan antara variabel dependen dan variabel independen dengan meminimalkan jumlah kuadrat galat. Model ini mengasumsikan bahwa setiap observasi bersifat independen antarwilayah, sehingga tidak mempertimbangkan adanya ketergantungan spasial. Oleh karena itu, OLS sering digunakan sebagai model dasar (baseline) sebelum dilakukan pengujian dan penerapan model spasial.
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
Penjelasan Komponen:
\(\mathbf{y}\) : vektor variabel dependen
\(\mathbf{X}\) : matriks variabel independen
\(\boldsymbol{\beta}\) : vektor parameter regresi
\(\boldsymbol{\varepsilon}\) : vektor galat (error)
Spatial Autoregressive (SAR)
Spatial Autoregressive Model (SAR) digunakan ketika nilai variabel dependen pada suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah sekitarnya. Model ini secara eksplisit memasukkan ketergantungan spasial melalui lag spasial dari variabel dependen, sehingga mampu menangkap interaksi antarwilayah.
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
Penjelasan Komponen:
\(\rho\) : koefisien autoregresif spasial
\(\mathbf{W}\) : matriks bobot spasial
\(\mathbf{W}\mathbf{y}\) : lag spasial dari variabel dependen
Komponen lainnya sama dengan model OLS
Spatial Error Model (SEM)
Spatial Error Model (SEM) digunakan ketika ketergantungan spasial muncul pada komponen galat, bukan pada variabel dependen. Model ini sesuai untuk kondisi di mana terdapat faktor-faktor yang tidak teramati namun bersifat spasial sehingga memengaruhi nilai galat antarwilayah.
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{u} \]
\[ \boldsymbol{u} = \lambda \mathbf{W}\boldsymbol{u} + \boldsymbol{\varepsilon} \]
Penjelasan Komponen:
\(\boldsymbol{u}\) : galat yang mengandung ketergantungan spasial
\(\lambda\) : koefisien error spasial
\(\boldsymbol{\varepsilon}\) : galat acak
Spatial Durbin Model (SDM)
Spatial Durbin Model (SDM) merupakan pengembangan dari model SAR dengan menambahkan lag spasial pada variabel independen. Model ini mampu menangkap efek langsung dan tidak langsung (spillover effects) antarwilayah, sehingga memberikan gambaran yang lebih komprehensif mengenai pengaruh spasial.
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \]
Penjelasan Komponen:
\(\mathbf{W}\mathbf{X}\) : lag spasial dari variabel independen
\(\boldsymbol{\theta}\) : parameter efek spillover
Komponen lainnya sama dengan model SAR
Spatial Autoregressive Combined Model (SAC)
Spatial Autoregressive Combined Model (SAC) mengombinasikan karakteristik model SAR dan SEM dengan memasukkan ketergantungan spasial pada variabel dependen dan komponen galat secara simultan. Model ini digunakan ketika proses spasial terjadi melalui interaksi antarwilayah dan pengaruh faktor laten yang bersifat spasial.
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{u} \]
\[ \boldsymbol{u} = \lambda \mathbf{W}\boldsymbol{u} + \boldsymbol{\varepsilon} \]
Penjelasan Komponen:
\(\rho\) : koefisien lag spasial variabel dependen
\(\lambda\) : koefisien error spasial
\(\mathbf{W}\) : matriks bobot spasial
Komponen lainnya sama dengan model SAR dan SEM
Metode Interpolasi - Ordinary Kriging (OK)
Ordinary Kriging (OK) adalah metode interpolasi geostatistik yang menghasilkan taksiran linier tak bias pada lokasi yang belum diamati dengan meminimalkan variansi taksiran. OK memanfaatkan struktur ketergantungan spasial yang diwakili oleh variogram atau fungsi kovarians sehingga bobot observasi disusun berdasarkan hubungan spasial (jarak dan model variogram) antara titik-titik pengamatan dan lokasi target. Keunggulan OK adalah sifatnya sebagai Best Linear Unbiased Estimator (BLUE) — memberikan estimasi dengan variansi minimum di antara semua taksiran linier tak bias — serta kemampuannya menghasilkan peta ketidakpastian (kriging variance) yang berguna untuk interpretasi keandalan taksiran.
Secara praktis, OK mengasumsikan bahwa rata-rata lokal tidak diketahui tetapi konstan dalam neighborhood yang dipertimbangkan (mean lokal konstan namun tidak perlu diketahui nilainya secara global). Karena itu OK menggabungkan constraint tak bias (jumlah bobot = 1) melalui metode Lagrange multiplier saat menyelesaikan sistem persamaan normal yang diperoleh dari peminimalan variansi taksiran. Prosesnya mencakup pemodelan variogram empiris, pemilihan fungsi variogram teoretis (mis. exponential, spherical, gaussian), dan penyelesaian sistem kriging untuk mendapatkan bobot yang kemudian dipakai menghasilkan peta interpolasi dan peta variansnya.
Taksiran Ordinary Kriging pada titik \(s_0\) dinyatakan sebagai kombinasi linier observasi: \[ \hat{Z}(s_0) \;=\; \sum_{i=1}^{n} \lambda_i \, Z(s_i) \]
dengan syarat tak-bias:
\[ \sum_{i=1}^{n} \lambda_i = 1. \]
Bobot \(\lambda_i\) diperoleh dengan menyelesaikan sistem kriging (menggunakan kovarians \(C(\cdot,\cdot)\) atau variogram \(\gamma(\cdot)\)). Dalam bentuk kovarians, sistem linear dengan Lagrange multiplier \(\mu\) adalah:
\[ \begin{cases} \displaystyle \sum_{j=1}^{n} \lambda_j \, C(s_i,s_j) + \mu = C(s_i,s_0), & i=1,\dots,n\\[6pt] \displaystyle \sum_{j=1}^{n} \lambda_j = 1. \end{cases} \]
Atau dalam bentuk variogram (intrinsic hypothesis) menjadi:
\[ \begin{cases} \displaystyle \sum_{j=1}^{n} \lambda_j \, \gamma(s_i,s_j) + \mu = \gamma(s_i,s_0), & i=1,\dots,n\\[6pt] \displaystyle \sum_{j=1}^{n} \lambda_j = 1. \end{cases} \]
Kriging variance (varians taksiran) diberikan oleh:
\[ \sigma^2_K(s_0) \;=\; \operatorname{Var}\big(Z(s_0)-\hat{Z}(s_0)\big) \;=\; \sum_{i=1}^{n} \lambda_i \, C(s_i,s_0) + \mu, \]
atau ekuivalen dinyatakan melalui variogram sesuai bentuk yang dipakai.
Penjelasan Komponen:
\(Z(s_i)\) : nilai observasi pada lokasi \(s_i\).
\(\hat{Z}(s_0)\) : taksiran nilai pada lokasi target \(s_0\).
\(\lambda_i\) : bobot kriging untuk tiap observasi; ditentukan oleh struktur spasial (variogram/kovarians) dan constraint tak-bias.
\(\mathbf{W}\) (implisit) : matriks kovarians/variogram antar titik yang membentuk koefisien sistem linear.
\(\mu\) : Lagrange multiplier yang memastikan kondisi tak-bias (\(\sum \lambda_i = 1\)).
\(C(h)\) atau \(\gamma(h)\) : fungsi kovarians atau variogram (tergantung pendekatan); parameter model meliputi nugget, sill, dan range:
Nugget: komponen variabilitas pada jarak sangat kecil / kesalahan pengukuran.
Sill: nilai variogram pada plateau yang mencerminkan varians total (untuk kovarians, nilai sill berkaitan dengan varians proses).
Range: jarak di mana korelasi spasial menjadi tidak signifikan (batas pengaruh spasial).
Kriging variance \(\sigma^2_K\): ukuran ketidakpastian taksiran pada \(s_0\); berguna untuk menilai keandalan peta interpolasi.
Neighborhood selection: pemilihan titik pengamatan sekitar \(s_0\) (mis. k terdekat atau radius tertentu) yang mempengaruhi bobot dan efisiensi perhitungan.
Metodologi Penelitian
Data Penelitian
Penelitian ini menggunakan data sekunder yang bersumber dari Badan Pusat Statistik (BPS) Sumatra Utara serta instansi resmi terkait. Data yang digunakan bersifat cross-sectional serta menunjukkan kondisi kabupaten/kota di Provinsi Sumatera Utara dengan satu periode pengamatan, yaitu tahun 2024. Variabel-variabel yang digunakan dapat dilihat sebagai berikut:
| No | Nama Variabel | Simbol | Jenis Variabel | Skala Pengukuran | Satuan / Keterangan |
|---|---|---|---|---|---|
| 1 | Tuberkulosis (TBC) | \(Y\) | Dependen | Rasio | Jumlah / angka kejadian TBC per kabupaten/kota |
| 2 | Angka Kesembuhan | \(X_1\) | Independen | Rasio | Persentase (%) |
| 3 | % RT dengan Sumber Air Minum Layak | \(X_2\) | Independen | Rasio | Persentase (%) |
| 4 | % RT dengan Sanitasi Layak | \(X_3\) | Independen | Rasio | Persentase (%) |
| 5 | Kepadatan Penduduk | \(X_4\) | Independen | Rasio | Jiwa per km\(^2\) |
| 6 | Jumlah Puskesmas | \(X_5\) | Independen | Rasio | Unit fasilitas kesehatan |
Unit Spasial
Unit spasial dalam penelitian ini adalah wilayah kabupaten/kota di Provinsi Sumatra Utara. Setiap objek dianggap sebagai satuan unit observasi spasial yang memiliki batas geografis yang jelas sesuai dengan pembagian administrasi resmi. Penggunaan unit spasial ini dipilih karena ketersediaan data epidemiologi TBC dan data pendukung lainnya umumnya disajikan pada tingkat ini, serta reevan dengan pelaksanaan dan evaluasi kebijakan kesehatan daerah.
Metode Analisis
Metode analisis dalam penelitian ini disusun secara bertahap untuk menjawab tujuan penelitian yang telah dirumuskan. Secata umum, analisis mencakup sebagai berikut:
Tahapan pertama yaitu pengumpulan data sekunder. Pada tahap ini dikumpulkan seluruh data sekunder yang diperlukan yaitu peta batas administrasi (shapefile GADM level kabupaten/kota), data kasus Tuberkulosis tahun 2024, dan data pendukung terkait determinan lingkungan/layanan (mis. angka kesembuhan, akses air minum, sanitasi, kepadatan penduduk, jumlah puskesmas). Sumber data dicatat secara rinci (instansi/dataset, tahun, format) dan dilakukan pengecekan awal kelengkapan serta konsistensi (kesesuaian nama wilayah, satuan, dan cakupan spasial) sebelum diproses lebih lanjut.
Tahapan kedua yaitu Eksplorasi data dan analisis awal. Eksplorasi meliputi pemeriksaan kualitas data (missing values, outlier, konsistensi nama kabupaten), transformasi jika perlu (mis. normalisasi atau log-transform untuk variabel skewed), serta deskriptif statistik dan visualisasi awal (histogram, boxplot, peta choropleth). Untuk aspek spasial dilakukan Exploratory Spatial Data Analysis (ESDA) termasuk perhitungan Moran’s I global dan peta LISA untuk mendeteksi autokorelasi spasial dan klaster awal; hasil ESDA memandu keputusan pemilihan model dan perlakuan data selanjutnya.
Tahapan ketiga yaitu Permodelan regresi Non-spasial menggunakan OLS untuk menguji hubungan antara variabel dependen dengan independen. Pada tahap ini dilakukan pengecekan asumsi OLS (linearitas, homoskedastisitas, normalitas residual, multikolinearitas) serta pemeriksaan residual untuk mendeteksi adanya autokorelasi spasial; bila residual OLS menunjukkan pola spasial signifikan, hal tersebut menjadi justifikasi penerapan model spasial pada tahap berikutnya.
Tahapan keempat yaitu Permodelan Regresi Spasial. Langkah ini mengestimasi model-model spasial untuk mengakomodasi ketergantungan ruang yang terdeteksi: SAR ketika efek lag dependen dominan, SEM bila ketergantungan muncul pada error, SDM untuk menangkap spillover variabel independen, dan SAC jika efek lag dependen dan error muncul simultan.
Tahapan kelima yaitu Analisis interpolasi Spasial menggunakan Ordinary Kinging (OK). Hasil OK memberikan gambaran kontinuitas spasial dan ketidakpastian taksiran yang berguna untuk mengidentifikasi hotspot potensial, namun interpretasi kuantitatif harus mempertimbangkan bahwa data penelitian berbentuk agregat kabupaten/kota.
Tahap akhir yaitu menyajikan hasil dalam bentuk peta tematik (choropleth, peta klaster LISA, peta kriging dan peta varians), tabel ringkasan parameter model, serta grafik diagnostik; interpretasi menekankan temuan utama (zona risiko tinggi, efek spillover, signifikansi variabel) dan implikasi kebijakan.
Alur Kerja Penelitian
Penelitian ini diawali dengan pengumpulan data sekunder yang meliputi data kasus Tuberkulosis tahun 2024 dan data spasial batas administrasi kabupaten/kota di Provinsi Sumatera Utara. Selanjutnya dilakukan eksplorasi dan pembersihan data untuk memastikan kualitas serta konsistensi antar sumber. Tahap berikutnya adalah analisis awal melalui statistik deskriptif dan eksplorasi spasial guna mengidentifikasi pola sebaran dan indikasi ketergantungan spasial. Berdasarkan hasil tersebut, dilakukan permodelan regresi non-spasial sebagai model dasar, kemudian dilanjutkan dengan permodelan regresi spasial untuk mengakomodasi efek ketergantungan ruang. Sebagai pelengkap analisis, dilakukan interpolasi spasial untuk memvisualisasikan kontinuitas risiko antarwilayah. Seluruh hasil kemudian divisualisasikan dan diinterpretasikan untuk menarik kesimpulan serta rekomendasi kebijakan pengendalian Tuberkulosis berbasis wilayah.
Hasil dan Pembahasan
Statistika Deskriptif
Terlebih dahulu dilakukan statistika deskriptif untuk memberikan gambaran awal mengenai karakteristik data kasus Tuberkulosis serta variabel penjelas yang digunakan dalam analisis, yaitu Angka Kesembuhan, Sumber Air Minum Layak, Sanitasi Layak, Jumlah Puskesmas, serta Kepadatan Penduduk pada tingkat kabupaten/kota di Provinsi Sumatra Utara.
# Import Data
library(openxlsx)## Warning: package 'openxlsx' was built under R version 4.3.3data = read.xlsx("C:/Users/MSI Cyborg15/OneDrive/Documents/college/Semester 5/Spatial/uasss/data_tbc_sumut.xlsx")
head(data)## Kabupaten.Kota Y X1 X2 X3 X4 X5
## 1 Asahan 1865 1206 93.94 92.13 221 31
## 2 Batu Bara 923 733 97.21 92.54 500 15
## 3 Dairi 1117 996 92.31 94.42 161 18
## 4 Deli Serdang 5972 4640 95.51 96.78 805 37
## 5 Humbang Hasundutan 482 286 94.94 96.99 89 12
## 6 Karo 1222 701 93.13 88.11 196 19# Statistika Deskriptif
summary(data)## Kabupaten.Kota Y X1 X2
## Length:33 Min. : 49 Min. : 106 Min. : 51.33
## Class :character 1st Qu.: 610 1st Qu.: 405 1st Qu.: 75.78
## Mode :character Median : 775 Median : 647 Median : 90.83
## Mean : 1571 Mean :1073 Mean : 85.16
## 3rd Qu.: 1379 3rd Qu.: 993 3rd Qu.: 96.04
## Max. :17161 Max. :8950 Max. :100.00
## X3 X4 X5
## Min. :15.63 Min. : 42 Min. : 5.00
## 1st Qu.:55.37 1st Qu.: 96 1st Qu.:12.00
## Median :88.11 Median : 178 Median :17.00
## Mean :74.53 Mean :1171 Mean :18.91
## 3rd Qu.:92.96 3rd Qu.: 707 3rd Qu.:21.00
## Max. :96.99 Max. :8945 Max. :47.00Berdasarkan ringkasan statistik deskriptif pada 33 kabupaten/kota, terlihat bahwa jumlah kasus TBC (Y) memiliki variasi yang sangat besar, dengan nilai maksimum jauh lebih tinggi dibandingkan median dan rata-ratanya. Hal ini menunjukkan adanya ketimpangan beban TBC antarwilayah, di mana hanya beberapa kabupaten/kota yang memiliki jumlah kasus sangat tinggi. Variabel angka kesembuhan (X1) dan kepadatan penduduk (X4) juga menunjukkan pola serupa, dengan distribusi yang menceng ke kanan akibat keberadaan nilai ekstrem. Sementara itu, akses air minum layak (X2) dan sanitasi layak (X3) pada umumnya berada pada tingkat yang cukup tinggi, meskipun masih terdapat beberapa wilayah dengan capaian yang relatif rendah. Variabel jumlah puskesmas (X5) memiliki sebaran yang lebih moderat dibandingkan variabel lainnya. Secara keseluruhan, data menunjukkan heterogenitas antarwilayah yang cukup kuat, sehingga analisis lanjutan, khususnya dengan pendekatan spasial, menjadi relevan untuk menangkap perbedaan karakteristik antar kabupaten/kota.
# Variabel Numerik
vars = data[, c("Y", "X1", "X2", "X3", "X4", "X5")]
# Layout Panel
par(mfrow = c(2, 3), # 2 baris, 3 kolom
mar = c(4, 4, 2, 1))
# Histogram
for (i in colnames(vars)) {
hist(vars[[i]],
main = paste("Histogram", i),
xlab = i,
col = "lightblue",
border = "white")
}
# Boxplot
for (i in colnames(vars)) {
boxplot(vars[[i]],
main = paste("Histogram", i),
xlab = i,
col = "lightblue",
border = "white")
}
Secara keseluruhan, data menunjukkan distribusi yang tidak merata dan cenderung tidak simetris, dengan sebagian besar nilai terkonsentrasi pada tingkat rendah–menengah dan beberapa nilai ekstrem. Hal ini menandakan adanya perbedaan karakteristik antarwilayah serta potensi outlier, sehingga analisis spasial menjadi penting untuk memahami pola sebarannya.
# Library
library(sp)## Warning: package 'sp' was built under R version 4.3.3library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(spdep)## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tmap)
# Shapefile
Indo = st_read("C:/Users/MSI Cyborg15/OneDrive/Documents/college/Semester 5/Spatial/uasss/gadm41_IDN_2.json")## Reading layer `gadm41_IDN_2' from data source
## `C:\Users\MSI Cyborg15\OneDrive\Documents\college\Semester 5\Spatial\uasss\gadm41_IDN_2.json'
## using driver `GeoJSON'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.0102 ymin: -11.0075 xmax: 141.0194 ymax: 6.0766
## Geodetic CRS: WGS 84# Filter Sumatera Utara
sumut = Indo[Indo$NAME_1 == "SumateraUtara",]
# Hapus Lake Toba (bukan kab/kota)
sumut <- sumut %>%
filter(NAME_2 != "LakeToba")
# Penyesuaian NAma
sumut$NAME_2 <- dplyr::recode(
sumut$NAME_2,
"Asahan" = "Asahan",
"BatuBara" = "Batu Bara",
"Dairi" = "Dairi",
"DeliSerdang" = "Deli Serdang",
"Gunungsitoli" = "Kota Gunungsitoli",
"HumbangHasundutan" = "Humbang Hasundutan",
"Karo" = "Karo",
"KotaBinjai" = "Kota Binjai",
"KotaMedan" = "Kota Medan",
"KotaTanjungbalai" = "Kota Tanjung Balai",
"Labuhanbatu" = "Labuhan Batu",
"LabuhanbatuSelatan" = "Labuhan Batu Selatan",
"LabuhanbatuUtara" = "Labuhan Batu Utara",
"Langkat" = "Langkat",
"MandailingNatal" = "Mandailing Natal",
"Nias" = "Nias",
"NiasBarat" = "Nias Barat",
"NiasSelatan" = "Nias Selatan",
"NiasUtara" = "Nias Utara",
"PadangLawas" = "Padang Lawas",
"PadangLawasUtara" = "Padang Lawas Utara",
"Padangsidimpuan" = "Kota Padangsidimpuan",
"PakpakBarat" = "Pakpak Bharat",
"Pematangsiantar" = "Kota Pematang Siantar",
"Samosir" = "Samosir",
"SerdangBedagai" = "Serdang Bedagai",
"Sibolga" = "Kota Sibolga",
"Simalungun" = "Simalungun",
"TapanuliSelatan" = "Tapanuli Selatan",
"TapanuliTengah" = "Tapanuli Tengah",
"TapanuliUtara" = "Tapanuli Utara",
"Tebingtinggi" = "Kota Tebing Tinggi",
"TobaSamosir" = "Toba Samosir"
)
# Konversi ke sf
sumut_sf <- st_as_sf(sumut)
# Join dengan data TBC
sumut_merged <- sumut_sf %>%
left_join(
data,
by = c("NAME_2" = "Kabupaten.Kota")
)
# Peta Sebaran
vars_num <- sumut_merged %>%
st_drop_geometry() %>%
select(where(is.numeric)) %>%
names()
library(ggplot2)
for (v in vars_num) {
p <- ggplot(sumut_merged) +
geom_sf(aes(fill = .data[[v]]), color = "white", linewidth = 0.2) +
scale_fill_viridis_c(option = "C", na.value = "grey90") +
labs(
title = paste("Sebaran", v, "di Sumatra Utara"),
fill = v
) +
theme_minimal()
print(p)
}
Secara umum, peta sebaran menunjukkan bahwa setiap variabel memiliki pola spasial yang berbeda antarwilayah di Sumatera Utara, dengan adanya konsentrasi nilai tinggi dan rendah yang tidak tersebar secara acak. Beberapa wilayah tampak membentuk klaster nilai serupa, sementara wilayah lain menunjukkan kontras yang cukup jelas dengan daerah sekitarnya. Pola ini mengindikasikan adanya ketimpangan spasial dan memperkuat bahwa karakteristik wilayah berperan penting dalam menjelaskan variasi data, sehingga pendekatan analisis spasial memang relevan untuk digunakan.
Autokorelasi Spasial
# Pastikan CRS planar (bukan longlat)
sumut_merged <- st_transform(sumut_merged, 32748) # UTM zona Sumut
# Neighbour berbasis queen contiguity
nb <- poly2nb(sumut_merged, queen = TRUE)## Warning in poly2nb(sumut_merged, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.# Bobot spasial (row-standardized)
WL <- nb2listw(nb, style = "W", zero.policy = TRUE)
## Moran's I
moran_test <- moran.test(
sumut_merged$Y,
WL,
zero.policy = TRUE
)
moran_test##
## Moran I test under randomisation
##
## data: sumut_merged$Y
## weights: WL
##
## Moran I statistic standard deviate = 5.3952, p-value = 3.422e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.349522477 -0.031250000 0.004980973## Plot
moran.plot(
sumut_merged$Y,
WL,
zero.policy = TRUE,
main = "Moran Scatterplot Kasus TBC\nProvinsi Sumatera Utara, 2024"
)
Hasil uji Moran’s I menunjukkan nilai statistik Moran sebesar 0,35 dengan p-value yang sangat kecil, sehingga dapat disimpulkan bahwa terdapat autokorelasi spasial positif yang signifikan pada data kasus TBC di Provinsi Sumatera Utara. Artinya, wilayah dengan jumlah kasus tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki kasus tinggi, begitu pula wilayah dengan kasus rendah yang saling berdekatan. Pola ini mengindikasikan bahwa sebaran kasus TBC tidak terjadi secara acak, melainkan membentuk pola spasial tertentu.
Hal tersebut diperkuat oleh Moran scatterplot, di mana sebagian besar wilayah berada pada kuadran High–High dan Low–Low. Kondisi ini menunjukkan adanya pengelompokan wilayah dengan karakteristik serupa dalam hal jumlah kasus TBC. Dengan demikian, analisis lanjutan berbasis spasial menjadi penting karena faktor kedekatan wilayah berperan dalam menjelaskan variasi kasus TBC antar daerah.
## Geary's C
geary_test <- geary.test(
sumut_merged$Y,
WL,
zero.policy = TRUE
)
geary_test##
## Geary C test under randomisation
##
## data: sumut_merged$Y
## weights: WL
##
## Geary C statistic standard deviate = 2.2069, p-value = 0.01366
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.33573060 1.00000000 0.09059621Hasil uji Geary’s C menunjukkan nilai statistik sebesar 0,336, yang lebih kecil dari nilai ekspektasinya (1), dengan p-value 0,01366 sehingga signifikan pada taraf 5%. Hal ini mengindikasikan adanya autokorelasi spasial positif pada data kasus TBC, di mana wilayah-wilayah yang berdekatan cenderung memiliki nilai kasus yang mirip. Dengan kata lain, sebaran kasus TBC tidak bersifat acak, melainkan membentuk pola pengelompokan spasial antar wilayah yang saling bertetangga.
## LISA (Local Moran's I)
lisa <- localmoran(
sumut_merged$Y,
WL,
zero.policy = TRUE
)
# Simpan ke data
sumut_merged$Ii <- lisa[, "Ii"]
sumut_merged$P_Ii <- lisa[, "Pr(z != E(Ii))"]
sumut_merged$Z_Ii <- lisa[, "Z.Ii"]
# Standarisasi
y_std <- scale(sumut_merged$Y)
wy_std <- lag.listw(WL, y_std, zero.policy = TRUE)
sumut_merged$LISA_Cluster <- NA
sumut_merged$LISA_Cluster[y_std > 0 & wy_std > 0 & sumut_merged$P_Ii < 0.05] <- "High-High"
sumut_merged$LISA_Cluster[y_std < 0 & wy_std < 0 & sumut_merged$P_Ii < 0.05] <- "Low-Low"
sumut_merged$LISA_Cluster[y_std > 0 & wy_std < 0 & sumut_merged$P_Ii < 0.05] <- "High-Low"
sumut_merged$LISA_Cluster[y_std < 0 & wy_std > 0 & sumut_merged$P_Ii < 0.05] <- "Low-High"
sumut_merged$LISA_Cluster[is.na(sumut_merged$LISA_Cluster)] <- "Not Significant"
tmap_mode("plot")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`# Matikan autoscale info
tmap_options(component.autoscale = FALSE)
tm_shape(sumut_merged) +
tm_polygons(
fill = "LISA_Cluster",
fill.scale = tm_scale_categorical(
values = c(
"High-High" = "#B2182B", # hotspot
"Low-Low" = "#2166AC", # coldspot
"High-Low" = "#EF8A62", # outlier
"Low-High" = "#67A9CF",
"Not Significant" = "grey80"
)
),
col = "white",
lwd = 0.4
) +
tm_layout(
legend.outside = TRUE,
frame = FALSE
)
# Hitung Getis–Ord Gi*
gi_star <- localG(
x = sumut_merged$Y,
listw = WL,
zero.policy = TRUE
)
# Simpan ke data
sumut_merged$GiZ <- as.numeric(gi_star)
sumut_merged <- sumut_merged %>%
mutate(
Gi_Class = case_when(
GiZ >= 2.58 ~ "Hotspot (99%)",
GiZ >= 1.96 ~ "Hotspot (95%)",
GiZ >= 1.65 ~ "Hotspot (90%)",
GiZ <= -2.58 ~ "Coldspot (99%)",
GiZ <= -1.96 ~ "Coldspot (95%)",
GiZ <= -1.65 ~ "Coldspot (90%)",
TRUE ~ "Not Significant"
)
)
library(tmap)
tmap_mode("plot")## ℹ tmap modes "plot" - "view"tmap_options(component.autoscale = FALSE)
tm_shape(sumut_merged) +
tm_polygons(
fill = "Gi_Class",
fill.scale = tm_scale_categorical(
values = c(
"Hotspot (99%)" = "#7F0000",
"Hotspot (95%)" = "#B2182B",
"Hotspot (90%)" = "#EF8A62",
"Coldspot (99%)" = "#08306B",
"Coldspot (95%)" = "#2166AC",
"Coldspot (90%)" = "#67A9CF",
"Not Significant" = "grey85"
)
),
col = "white",
lwd = 0.4
) +
tm_layout(
legend.outside = TRUE,
frame = FALSE
)
Estimasi Model OLS
Estimasi Parameter OLS
ols = lm(
Y ~
Y + X1 + X2 + X3 + X4 + X5,
data = sumut_merged
)## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assignedsummary(ols)##
## Call:
## lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data = sumut_merged)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1549.82 -206.95 54.86 259.86 867.41
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 554.71891 547.54023 1.013 0.3200
## X1 1.79212 0.09987 17.945 <2e-16 ***
## X2 -7.95741 9.84773 -0.808 0.4261
## X3 -1.19752 5.30813 -0.226 0.8232
## X4 0.13435 0.05956 2.256 0.0324 *
## X5 -15.69078 13.81902 -1.135 0.2662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 471 on 27 degrees of freedom
## Multiple R-squared: 0.9791, Adjusted R-squared: 0.9752
## F-statistic: 252.5 on 5 and 27 DF, p-value: < 2.2e-16Hasil estimasi model OLS menunjukkan bahwa secara keseluruhan model memiliki kinerja yang sangat baik. Hal ini ditunjukkan oleh nilai R-squared sebesar 0,979 dan Adjusted R-squared sebesar 0,975, yang berarti sebagian besar variasi pada variabel respon dapat dijelaskan oleh variabel-variabel penjelas dalam model. Selain itu, hasil uji F yang signifikan (p-value < 0,001) mengindikasikan bahwa model secara simultan signifikan dan layak digunakan untuk analisis hubungan antara variabel respon dan variabel independen.
Secara parsial, hanya X1 dan X4 yang berpengaruh signifikan terhadap variabel respon, dengan arah pengaruh positif. Variabel X1 memiliki pengaruh yang sangat kuat dan signifikan, sedangkan X4 berpengaruh signifikan pada taraf 5%. Sementara itu, variabel X2, X3, dan X5 tidak menunjukkan pengaruh yang signifikan secara statistik. Meskipun model OLS menunjukkan hasil yang baik, adanya indikasi autokorelasi spasial pada analisis sebelumnya mengisyaratkan bahwa model spasial dapat memberikan estimasi yang lebih tepat dibandingkan OLS.
Uji Asumsi Klasik
Asumsi Multikolinieritas
library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(ols)## Warning in model.matrix.default(mod, data = structure(list(Y = c(1865, 923, :
## the response appeared on the right-hand side and was dropped## Warning in model.matrix.default(mod, data = structure(list(Y = c(1865, 923, :
## problem with term 1 in model.matrix: no columns are assigned## GVIF Df GVIF^(1/(2*Df))
## Y 21.278460 0 Inf
## X1 3.820684 1 1.954657
## X2 2.746567 1 1.657277
## X3 2.571929 1 1.603724
## X4 2.558103 1 1.599407
## X5 2.944981 1 1.716095Asumsi Normalitas
shapiro.test(residuals(ols))##
## Shapiro-Wilk normality test
##
## data: residuals(ols)
## W = 0.91518, p-value = 0.0135Asumsi Linearitas
library(lmtest)## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericharvtest(ols)## Warning in model.matrix.default(terms(formula), model.frame(formula)): the
## response appeared on the right-hand side and was dropped## Warning in model.matrix.default(terms(formula), model.frame(formula)): problem
## with term 1 in model.matrix: no columns are assigned##
## Harvey-Collier test
##
## data: ols
## HC = 1.178, df = 26, p-value = 0.2495Asumsi Homoskedastisitas
library(lmtest)
bptest(ols)## Warning in model.matrix.default(terms(formula), model.frame(formula)): the
## response appeared on the right-hand side and was dropped## Warning in model.matrix.default(terms(formula), model.frame(formula)): problem
## with term 1 in model.matrix: no columns are assigned##
## studentized Breusch-Pagan test
##
## data: ols
## BP = 12.769, df = 5, p-value = 0.02564Model Ekonometriks Spasial
Uji Lagrange Multiplier (LM)
lm.RStests(ols, listw = WL, zero.policy = TRUE, test = "all")## Warning in model.matrix.default(terms(model), model.frame(model)): the response
## appeared on the right-hand side and was dropped## Warning in model.matrix.default(terms(model), model.frame(model)): problem with
## term 1 in model.matrix: no columns are assigned##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data =
## sumut_merged)
## test weights: WL
##
## RSerr = 3.3478, df = 1, p-value = 0.0673
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data =
## sumut_merged)
## test weights: WL
##
## RSlag = 0.50002, df = 1, p-value = 0.4795
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data =
## sumut_merged)
## test weights: WL
##
## adjRSerr = 2.8483, df = 1, p-value = 0.09147
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data =
## sumut_merged)
## test weights: WL
##
## adjRSlag = 0.00055756, df = 1, p-value = 0.9812
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data =
## sumut_merged)
## test weights: WL
##
## SARMA = 3.3484, df = 2, p-value = 0.1875Hasil uji Lagrange Multiplier (Rao’s score) menunjukkan bahwa indikasi ketergantungan spasial pada model OLS relatif lemah. Uji LM error (RSerr) dan adjusted LM error hanya signifikan pada taraf 10% (p-value ≈ 0,07–0,09), sedangkan LM lag (RSlag) dan adjusted LM lag sama sekali tidak signifikan. Selain itu, uji SARMA juga tidak signifikan. Secara keseluruhan, hasil ini mengindikasikan bahwa tidak terdapat bukti kuat adanya ketergantungan spasial murni dalam bentuk lag maupun error pada model OLS, meskipun terdapat kecenderungan lemah ke arah spatial error, sehingga penggunaan model spasial tetap dapat dipertimbangkan sebagai pembanding atau penguatan analisis.
Pemilihan Model Terbaik
library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.3.3## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption# Spatial Autoregressive Model
sar = lagsarlm(
Y ~
Y + X1 + X2 + X3 + X4 + X5,
data = sumut_merged,
listw = WL,
method = "eigen"
)## Warning in model.matrix.default(mt, mf): the response appeared on the
## right-hand side and was dropped## Warning in model.matrix.default(mt, mf): problem with term 1 in model.matrix:
## no columns are assigned# Spatial Error Model
sem = errorsarlm(
Y ~
Y + X1 + X2 + X3 + X4 + X5,
data = sumut_merged,
listw = WL,
method = "eigen"
)## Warning in model.matrix.default(mt, mf): the response appeared on the
## right-hand side and was dropped
## Warning in model.matrix.default(mt, mf): problem with term 1 in model.matrix:
## no columns are assigned# Spatial Durbin Model
sdm = lagsarlm(
Y ~
Y + X1 + X2 + X3 + X4 + X5,
data = sumut_merged,
listw = WL,
Durbin = TRUE,
method = "eigen"
)## Warning in model.matrix.default(mt, mf): the response appeared on the
## right-hand side and was dropped
## Warning in model.matrix.default(mt, mf): problem with term 1 in model.matrix:
## no columns are assigned# Spatial Autoreggresive Combined Model
sac = sacsarlm(
Y ~
Y + X1 + X2 + X3 + X4 + X5,
data = sumut_merged,
listw = WL,
method = "eigen"
)## Warning in model.matrix.default(mt, mf): the response appeared on the
## right-hand side and was dropped
## Warning in model.matrix.default(mt, mf): problem with term 1 in model.matrix:
## no columns are assigned## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped## Warning in model.matrix.default(mt, mf, contrasts): problem with term 1 in
## model.matrix: no columns are assigned# Perbandingan Model
aic_values = c(
SAR = AIC(sar),
SEM = AIC(sem),
SDM = AIC(sdm),
SAC = AIC(sac)
)
aic_values[order(aic_values)]## SDM SAC SEM SAR
## 487.9999 493.1934 503.2354 508.4059Berdasarkan perbandingan nilai Akaike Information Criterion (AIC), model Spatial Durbin Model (SDM) memiliki nilai AIC paling kecil dibandingkan model lainnya, yaitu sebesar 487,99, diikuti oleh SAC, SEM, dan SAR. Nilai AIC yang lebih kecil menunjukkan bahwa model tersebut memiliki keseimbangan terbaik antara ketepatan model dan kompleksitas parameter, sehingga lebih efisien dalam menjelaskan variasi data.
Dengan demikian, SDM dapat dinyatakan sebagai model permodelan terbaik untuk data yang dianalisis. Hal ini mengindikasikan bahwa selain pengaruh variabel penjelas secara langsung, pengaruh spasial dari wilayah sekitar (spillover effect) juga berperan penting dalam menjelaskan variasi variabel respon. Oleh karena itu, penggunaan SDM lebih tepat dibandingkan model spasial lainnya maupun model OLS dalam menggambarkan pola dan hubungan spasial pada data ini.
Model Terbaik: Spatial Durbin Model (SDM)
summary(sdm)##
## Call:lagsarlm(formula = Y ~ Y + X1 + X2 + X3 + X4 + X5, data = sumut_merged,
## listw = WL, Durbin = TRUE, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -535.614 -140.862 -12.579 133.904 411.961
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1299.672081 523.112447 2.4845 0.01297
## X1 1.871468 0.068809 27.1979 < 2.2e-16
## X2 -9.100091 6.027126 -1.5099 0.13108
## X3 -1.813900 5.021345 -0.3612 0.71792
## X4 0.053820 0.038172 1.4100 0.15855
## X5 -5.072277 7.617558 -0.6659 0.50550
## lag.X1 1.233123 0.255212 4.8318 1.353e-06
## lag.X2 -18.051585 11.005373 -1.6403 0.10095
## lag.X3 8.808842 8.333009 1.0571 0.29046
## lag.X4 -0.490808 0.107132 -4.5814 4.620e-06
## lag.X5 6.504961 11.602570 0.5606 0.57504
##
## Rho: -0.86333, LR test value: 12.596, p-value: 0.00038661
## Asymptotic standard error: 0.13988
## z-value: -6.1721, p-value: 6.7409e-10
## Wald statistic: 38.094, p-value: 6.7409e-10
##
## Log likelihood: -231 for mixed model
## ML residual variance (sigma squared): 57170, (sigma: 239.1)
## Number of observations: 33
## Number of parameters estimated: 13
## AIC: 488, (AIC for lm: 498.6)
## LM test for residual autocorrelation
## test value: 1.6647, p-value: 0.19697Berdasarkan hasil estimasi Spatial Durbin Model (SDM), model ini menunjukkan adanya ketergantungan spasial yang signifikan, yang tercermin dari nilai parameter rho sebesar −0,8633 dengan p-value yang sangat kecil. Hal ini mengindikasikan bahwa nilai variabel respon di suatu wilayah dipengaruhi secara signifikan oleh kondisi wilayah sekitarnya, namun dengan arah hubungan yang berlawanan. Uji Likelihood Ratio dan Wald juga signifikan, sehingga memperkuat bahwa komponen spasial memang penting untuk dimasukkan dalam model. Selain itu, nilai AIC SDM lebih kecil dibandingkan model OLS, yang menandakan bahwa SDM memberikan kinerja pemodelan yang lebih baik.
Dari sisi pengaruh variabel, terdapat pengaruh langsung dan tidak langsung (spillover) yang berbeda-beda. Salah satu variabel menunjukkan pengaruh langsung yang signifikan terhadap variabel respon, serta diperkuat oleh pengaruh tidak langsung dari wilayah tetangga. Sebaliknya, beberapa variabel lain tidak signifikan baik secara langsung maupun spasial, yang menunjukkan bahwa pengaruhnya relatif lemah dalam menjelaskan variasi data. Uji autokorelasi residual yang tidak signifikan menandakan bahwa SDM telah mampu menangkap ketergantungan spasial dengan baik, sehingga residual model tidak lagi mengandung pola spasial yang berarti.
moran.test(residuals(sdm), WL, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: residuals(sdm)
## weights: WL
##
## Moran I statistic standard deviate = -0.60499, p-value = 0.7274
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.11416359 -0.03125000 0.01878268Hasil uji Moran’s I pada residual SDM menunjukkan nilai statistik Moran’s I sebesar −0,114 dengan p-value 0,7274, yang berarti tidak signifikan secara statistik. Hal ini mengindikasikan bahwa tidak terdapat autokorelasi spasial yang tersisa pada residual model, sehingga pola ketergantungan spasial telah berhasil ditangkap oleh Spatial Durbin Model. Dengan demikian, residual bersifat acak secara spasial dan model SDM dapat dinilai telah sesuai serta layak digunakan untuk analisis lebih lanjut.
Ordinary Kringing (OK)
Permodelan Struktur Variogram
Permodelan struktur variogram digunakan untuk menggambarkan dan mengukur ketergantungan spasial suatu data, yaitu sejauh mana nilai pada lokasi yang berdekatan memiliki kemiripan dibandingkan lokasi yang berjauhan. Variogram menunjukkan hubungan antara jarak antar lokasi (lag) dengan tingkat perbedaan nilai (semivarians). Secara umum, semakin dekat jaraknya, semakin kecil perbedaannya, dan semakin jauh jaraknya, perbedaan cenderung membesar hingga mencapai kondisi stabil.
Dalam struktur variogram terdapat tiga komponen utama, yaitu nugget, sill, dan range. Nugget merepresentasikan variasi pada jarak sangat kecil atau kesalahan pengukuran, sill menunjukkan nilai semivarians maksimum ketika pengaruh spasial mulai menghilang, dan range adalah jarak maksimum di mana ketergantungan spasial masih ada. Setelah variogram empiris dibentuk dari data, dilakukan pemilihan model variogram teoretis (seperti spherical, exponential, atau gaussian) yang paling sesuai. Model inilah yang kemudian digunakan dalam proses interpolasi spasial (misalnya kriging) untuk menghasilkan estimasi yang lebih akurat dan konsisten secara spasial.
Uji Anisotropi
# Plot Ketetanggaan
sumut_sp = as_Spatial(sumut_merged)
library(gstat)## Warning: package 'gstat' was built under R version 4.3.3v = variogram(Y ~ 1, data = sumut_sp, alpha = c(0, 45, 90, 135))
plot(v, maun = "Variogram Berdasarkan Arah")
Berdasarkan uji anisotropi melalui variogram arah (0°, 45°, 90°, dan 135°), terlihat bahwa pola kenaikan semivarians terhadap jarak tidak seragam di semua arah. Pada beberapa arah tertentu, semivarians meningkat lebih cepat dan mencapai nilai yang lebih tinggi dibandingkan arah lainnya. Hal ini menunjukkan bahwa tingkat kemiripan nilai antar lokasi berbeda-beda tergantung arah, sehingga terdapat indikasi anisotropi spasial pada data yang dianalisis.
Temuan ini mengimplikasikan bahwa proses spasial yang terjadi tidak bersifat isotropik (seragam ke segala arah), melainkan dipengaruhi oleh arah tertentu, yang bisa berkaitan dengan faktor geografis, lingkungan, atau pola distribusi wilayah. Oleh karena itu, dalam tahap permodelan variogram selanjutnya, penggunaan model variogram anisotropik menjadi lebih tepat dibandingkan model isotropik, agar struktur ketergantungan spasial dapat direpresentasikan secara lebih akurat dan hasil interpolasi menjadi lebih andal.
Fitting Variogram
dist_max = max(spDists(sumut_sp))
v_exp = variogram(Y ~ 1, data = sumut_sp, cutoff = dist_max/2, Width = dist_max/20)## Warning in variogram.default(y, locations, X, trend.beta = beta, grid = grid, :
## the following arguments are ignored: 17538.4085323374# Lag yang tidak stabil
v_exp_cleaned = v_exp[v_exp$np > 10, ]
plot(v_exp_cleaned, main = "Variogram Eksperimental")
Variogram eksperimental tersebut menggambarkan hubungan antara jarak antar lokasi dengan tingkat ketidaksamaan nilai (semivarians). Terlihat bahwa pada jarak yang relatif dekat, nilai semivarians masih rendah, yang menunjukkan bahwa lokasi-lokasi yang berdekatan cenderung memiliki nilai yang mirip. Seiring bertambahnya jarak, semivarians meningkat dan mulai berfluktuasi pada nilai yang lebih tinggi, menandakan bahwa kemiripan antar lokasi semakin berkurang ketika jaraknya makin jauh.
Pola ini menunjukkan adanya ketergantungan spasial hingga jarak tertentu, sebelum akhirnya mencapai kondisi relatif stabil (mendekati sill), meskipun masih tampak variasi karena keterbatasan jumlah data dan sebaran titik. Variogram eksperimental ini menjadi dasar penting untuk proses fitting variogram teoretis (misalnya model spherical, exponential, atau gaussian), agar struktur nugget, sill, dan range dapat ditentukan secara lebih jelas dan selanjutnya digunakan dalam analisis spasial lanjutan seperti interpolasi.
init_nugget <- 10
init_psill <- 15
init_range <- 100
m_sph <- fit.variogram(v_exp_cleaned, vgm(init_psill, "Sph", init_range, init_nugget))## Warning in fit.variogram(v_exp_cleaned, vgm(init_psill, "Sph", init_range, :
## singular model in variogram fitm_exp <- fit.variogram(v_exp_cleaned, vgm(init_psill, "Exp", init_range, init_nugget))## Warning in fit.variogram(v_exp_cleaned, vgm(init_psill, "Exp", init_range, :
## singular model in variogram fitm_gau <- fit.variogram(v_exp_cleaned, vgm(init_psill, "Gau", init_range, init_nugget))## Warning in fit.variogram(v_exp_cleaned, vgm(init_psill, "Gau", init_range, :
## singular model in variogram fitsse = c(
"SSE Spherical" = attr(m_sph, "SSErr"),
"SSE Exponential" = attr(m_exp, "SSErr"),
"SSE Gaussian" = attr(m_gau, "SSErr")
); sse## SSE Spherical SSE Exponential SSE Gaussian
## 2968897 2968897 2968897Hasil perhitungan Sum of Squared Errors (SSE) menunjukkan bahwa nilai SSE untuk model Spherical, Exponential, dan Gaussian adalah sama, yaitu sebesar 2.968.897. Ini berarti ketiga model variogram tersebut memiliki tingkat kecocokan yang setara terhadap variogram eksperimental, sehingga tidak ada satu model pun yang secara statistik lebih unggul dibandingkan yang lain berdasarkan kriteria SSE.
Dengan kondisi ini, pemilihan model variogram teoretis dapat didasarkan pada pertimbangan tambahan, seperti kemudahan interpretasi, karakteristik proses spasial yang dianalisis, atau kestabilan model saat digunakan untuk interpolasi (misalnya kriging). Dalam praktik, model spherical sering dipilih sebagai default karena interpretasinya sederhana dan umum digunakan, namun secara kuantitatif ketiga model pada hasil ini sama-sama layak digunakan.
Interpolasi Ordinary Kriging
Hasil Interpolasi
# Load Peta Kecamatan
indo_kec = readRDS("C:/Users/MSI Cyborg15/OneDrive/Documents/college/Semester 5/Spatial/gadm36_IDN_3_sp.rds")
sumut_kec = indo_kec[indo_kec$NAME_1 == "Sumatera Utara", ]
# Centroid dari Tiap Kabupaten
titik_kab = coordinates(sumut_sp)
sumut_points = SpatialPointsDataFrame(coords = titik_kab, data = sumut_sp@data)
proj4string(sumut_points) = CRS(as.character(NA))
# Centroid dari Tiap Kecamatan
titik_kec = coordinates(sumut_kec)
kec_pred_grid = SpatialPoints(titik_kec)
proj4string(kec_pred_grid) = CRS(as.character(NA))
# Prediksi Kriging ke Kecamatan
hasil_kec = krige(Y ~ 1,
sumut_points,
kec_pred_grid,
model = m_sph)## [using ordinary kriging]# Hasil Interpolasi
hasil_df = as.data.frame(hasil_kec)
hasil = data.frame(
Kabupaten_Kota = sumut_kec$NAME_2,
Kecamatan = sumut_kec$NAME_3,
Longitude = hasil_df$coords.x1,
Latitude = hasil_df$coords.x2,
TBC_Prediksi = hasil_df$var1.pred,
Error_var = hasil_df$var1.var
); head(hasil)## Kabupaten_Kota Kecamatan Longitude Latitude TBC_Prediksi Error_var
## 1 Asahan Aek Ledong 99.64753 2.611061 1570.515 25.75758
## 2 Asahan Aek Kuasan 99.71028 2.667194 1570.515 25.75758
## 3 Asahan Aek Songsongan 99.36953 2.605118 1570.515 25.75758
## 4 Asahan Air Batu 99.62501 2.863546 1570.515 25.75758
## 5 Asahan Air Joman 99.71566 3.012231 1570.515 25.75758
## 6 Asahan Bandar Pasir Mandoge 99.26253 2.735140 1570.515 25.75758library(viridis)## Warning: package 'viridis' was built under R version 4.3.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.3.3kec_sf = st_as_sf(sumut_kec)
st_crs(kec_sf)## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]kec_sf$TBC_Prediksi = hasil$TBC_Prediksi
library(ggplot2)
ggplot(data = kec_sf) +
geom_sf(
aes(fill = TBC_Prediksi),
color = "white",
linewidth = 0.2
) +
scale_fill_viridis_c(option = "C", na.value = "grey90") +
labs(
title = "Sebaran Prediksi TBC Tingkat Kecamatan Sumatera Utara",
fill = "Prediksi TBC"
) +
theme_minimal()
Peta sebaran prediksi TBC tingkat kecamatan di Provinsi Sumatera Utara menunjukkan bahwa nilai prediksi relatif homogen antar wilayah, dengan perbedaan spasial yang tidak terlalu kontras. Hal ini mengindikasikan bahwa model menghasilkan estimasi rata-rata yang cenderung merata, sehingga tidak tampak klaster ekstrem dengan nilai prediksi yang sangat tinggi maupun sangat rendah. Pola ini mencerminkan hasil pemodelan yang telah mengakomodasi ketergantungan spasial sehingga variasi lokal yang tajam menjadi lebih teredam.
Secara substantif, peta ini memberikan gambaran umum tingkat risiko TBC di seluruh kecamatan sebagai baseline prediksi, yang berguna untuk perencanaan kebijakan kesehatan secara makro. Meskipun variasinya tidak mencolok, hasil prediksi tetap dapat dimanfaatkan untuk mengidentifikasi wilayah yang memerlukan perhatian lebih lanjut, terutama jika dikombinasikan dengan informasi tambahan seperti kepadatan penduduk, akses layanan kesehatan, atau faktor lingkungan lainnya.
Kesimpulan dan Saran
Kesimpulan
Berdasarkan seluruh rangkaian analisis, dapat disimpulkan bahwa data menunjukkan adanya ketergantungan spasial yang signifikan, baik secara global maupun lokal, sehingga pendekatan regresi spasial lebih tepat dibandingkan regresi OLS biasa. Hasil pengujian dan perbandingan model menunjukkan bahwa Spatial Durbin Model (SDM) merupakan model terbaik, karena memiliki nilai AIC terendah serta mampu menangkap pengaruh langsung dan pengaruh spasial dari wilayah sekitar. Model ini juga telah memenuhi asumsi yang baik, ditunjukkan oleh tidak adanya autokorelasi spasial pada residual. Peta prediksi TBC memperlihatkan pola sebaran yang relatif merata antar kecamatan di Sumatera Utara, yang mencerminkan hasil pemodelan yang stabil dan telah mengakomodasi efek spasial.
Saran
Untuk penelitian selanjutnya, disarankan agar analisis spasial tuberkulosis diperluas dengan memasukkan variabel determinan tambahan, seperti kondisi sosial ekonomi, kepadatan penduduk, dan akses layanan kesehatan, sehingga pemodelan memberikan gambaran yang lebih utuh mengenai perbedaan kasus antarwilayah di Provinsi Sumatra Utara.
Referensi
Kementerian Kesehatan Republik Indonesia, Profil Kesehatan Indonesia 2024, Jakarta: Kemenkes RI, 2024.
M. Sobari, I. G. N. M. Jaya, and B. N. Ruchjana, “Spatial Analysis of Dengue Disease in Jakarta Province,” CAUCHY – Jurnal Matematika Murni dan Aplikasi, vol.7, no. 4, pp. 535–547, 2023.
L. Anselin, “Spatial Econometrics,” in A Companion to Theoretical Econometrics, B. H. Baltagi, Ed., Oxford: Blackwell Publishing, 2001, pp. 310–330.
H. I. Zebua and I. G. N. M. Jaya, “Spatial Autoregressive Model of Tuberculosis Cases in Central Java Province 2019,” CAUCHY – Jurnal Matematika Murni dan Aplikasi, vol. 7, no. 2, pp. 240–248, 2022.
H. Yasin, A. R. Hakim, and B. Warsito, Regresi Spasial (Aplikasi dengan R), Pekalongan: WADE Group, 2020.
W. Wang et al., “Reclaiming independence in spatial-clustering datasets: A series of data-driven spatial weights matrices,” Statistics in Medicine, vol. 41, no. 15, pp. 2939–2956, 2022.
A. A. Grasa, Econometric Model Selection: A New Approach, Dordrecht: Springer, 2018.
Lampiran
Dashboard: https://wnajochar2005.shinyapps.io/uasss/