Identitas Penulis:
| Keterangan | Informasi |
|---|---|
| Nama | Sri Yuliana |
| NPM | 140720250003 |
| Program Studi | Magister Statistika Terapan |
| Mata Kuliah | Epidemiologi |
Tuberkulosis (TB) masih menjadi masalah kesehatan masyarakat utama di Indonesia, dengan Jawa Barat sebagai provinsi penyumbang kasus tertinggi secara nasional. Penelitian ini bertujuan menganalisis distribusi spasial dan temporal kasus TB di Provinsi Jawa Barat periode 2021–2025. Data yang digunakan adalah data sekunder agregat 27 kabupaten/kota yang bersumber dari BPS dan Portal Open Data Jabar. Metode analisis meliputi pengukuran epidemiologi (Incidence Rate dan Standardized Incidence Ratio), pemetaan risiko berbasis Empirical Bayes, analisis autokorelasi spasial (Global Moran’s I dan LISA), serta pemodelan spasio-temporal menggunakan Generalized Linear Mixed Model Negative Binomial (GLMM-NB). Hasil analisis menunjukkan tren peningkatan IR dari 183,13 menjadi 456,22 per 100.000 penduduk selama 2021–2024. Pemetaan Empirical Bayes mengidentifikasi wilayah dengan risiko tinggi yang konsisten lintas tahun. Uji Global Moran’s I mengkonfirmasi adanya pengelompokan spasial yang signifikan. Model GLMM-NB menunjukkan bahwa HIV, kepadatan penduduk, dan kemiskinan berpengaruh positif terhadap laju kejadian TB, sementara sanitasi layak dan akses air minum berperan protektif. Temuan ini diharapkan dapat mendukung perencanaan intervensi TB yang lebih tepat sasaran berbasis wilayah di Jawa Barat.
Kata kunci: Tuberkulosis, Empirical Bayes, Disease Mapping, GLMM, Spasio-Temporal.
Tuberkulosis (TB) merupakan salah satu penyakit menular yang hingga saat ini masih menjadi masalah kesehatan masyarakat global. Penyakit ini disebabkan oleh bakteri Mycobacterium tuberculosis yang dapat ditularkan melalui udara ketika penderita TB melakukan aktivitas seperti batuk, bersin, atau meludah (World Health Organization [WHO], 2025). TB masih menjadi salah satu penyebab utama kematian akibat penyakit infeksi, dengan sebagian besar kasus menyerang paru-paru, meskipun dalam beberapa kondisi dapat juga menginfeksi organ tubuh lainnya (Nurdin et al., 2025; Turner et al., 2017).
Berdasarkan Global Tuberculosis Report tahun 2025, TB masih menjadi beban kesehatan global dengan lebih dari 10 juta kasus baru setiap tahun dan menyebabkan lebih dari 1 juta kematian. Pada tahun 2024, tercatat sekitar 8,3 juta kasus baru TB secara global, meningkat dibandingkan tahun 2023 yang mencapai 8,2 juta kasus. Indonesia termasuk negara dengan beban TB tertinggi di dunia setelah India, dengan kontribusi sekitar 10% dari total kasus global (WHO, 2025).
Di tingkat nasional, Profil Kesehatan Indonesia tahun 2024 melaporkan sebanyak 856.420 kasus TB, meningkat dibandingkan tahun 2023 yang sebanyak 821.200 kasus (Kementerian Kesehatan Republik Indonesia, 2025). Kasus TB di Indonesia didominasi oleh provinsi dengan jumlah penduduk besar seperti Jawa Barat, Jawa Timur, dan Jawa Tengah. Secara khusus, Jawa Barat merupakan provinsi dengan jumlah kasus TB tertinggi secara nasional. Tren kasus menunjukkan peningkatan berkelanjutan, yaitu 160.661 kasus pada tahun 2022, meningkat menjadi 211.959 kasus pada tahun 2023, dan mencapai 229.683 kasus pada tahun 2024 (Dinas Kesehatan Provinsi Jawa Barat, 2023; 2024; 2025).
Beberapa penelitian menunjukkan bahwa kejadian TB dipengaruhi oleh faktor sosial ekonomi dan lingkungan seperti kepadatan penduduk, kemiskinan, sanitasi, serta akses terhadap fasilitas kesehatan (Sihaloho et al., 2021; Suryani & Ibad, 2022; Laoli et al., 2024). Kepadatan penduduk menjadi faktor penting dalam meningkatkan risiko penularan karena memperbesar peluang kontak antarindividu dalam suatu wilayah (Kristianingrum, 2024). Selain itu, kondisi sanitasi dan lingkungan permukiman yang tidak layak juga berkontribusi terhadap tingginya angka kejadian TB (Atillah et al., 2023; Gityarani, 2024).
Pendekatan analisis spasial telah banyak digunakan untuk mengidentifikasi distribusi penyakit TB secara geografis. Studi sebelumnya menunjukkan bahwa terdapat variasi spasial yang signifikan dalam distribusi kasus TB antarwilayah, sehingga analisis spasial penting untuk menentukan wilayah prioritas intervensi (Haq et al., 2020; Rosady et al., 2024; Sasmita et al., 2017). Selain itu, studi berbasis ekologi juga menunjukkan bahwa determinan lingkungan dan sosial ekonomi memiliki pengaruh signifikan terhadap kejadian TB (Alma et al., 2024; Dzakiyah et al., 2023; Handayani et al., 2024).
Dalam analisis epidemiologi, ukuran seperti Incidence Rate (IR) dan Standardized Incidence Ratio (SIR) digunakan untuk menggambarkan perbedaan risiko antarwilayah secara lebih proporsional. Pendekatan SIR banyak digunakan dalam studi epidemiologi untuk mengidentifikasi wilayah dengan risiko tinggi terhadap penyakit menular termasuk TB (Waller & Gotway, 2004; Lawson, 2018).
Namun demikian, ukuran tersebut masih bersifat deskriptif dan belum mampu menangkap ketergantungan spasial maupun dinamika temporal secara simultan. Oleh karena itu, diperlukan pendekatan disease mapping untuk menghasilkan estimasi risiko yang lebih stabil. Metode Empirical Bayes banyak digunakan dalam pemetaan penyakit karena mampu mengurangi fluktuasi ekstrem pada wilayah dengan populasi kecil serta meningkatkan kestabilan estimasi risiko (Lawson, 2018; Chambers et al., 2013).
Selanjutnya, untuk mengidentifikasi pola pengelompokan spasial, digunakan analisis autokorelasi spasial melalui Global Moran’s I dan Local Indicators of Spatial Association (LISA), yang memungkinkan identifikasi klaster wilayah dengan risiko tinggi maupun rendah (Anselin, 1995; Anselin et al., 2006).
Lebih lanjut, untuk menangkap dinamika spasial dan temporal secara simultan, digunakan pendekatan spatio-temporal modeling melalui Generalized Linear Mixed Model (GLMM). Model ini memungkinkan adanya efek acak untuk wilayah dan waktu sehingga mampu merepresentasikan heterogenitas data panel epidemiologi (Zuur et al., 2009).
Namun, data kasus TB sebagai data hitungan sering menunjukkan fenomena overdispersion, yaitu kondisi ketika varians lebih besar daripada rata-rata. Kondisi ini menyebabkan model Poisson menjadi kurang sesuai karena melanggar asumsi dasar kesamaan mean dan varians. Oleh karena itu, digunakan pendekatan Negative Binomial GLMM sebagai model yang lebih fleksibel dalam menangani overdispersi sekaligus mempertahankan struktur efek acak spasial dan temporal (Hilbe, 2011).
Dengan demikian, pendekatan Negative Binomial GLMM dalam penelitian ini diharapkan mampu memberikan hasil pemodelan yang lebih akurat dalam menjelaskan faktor-faktor yang memengaruhi kasus Tuberkulosis serta menangkap variasi risiko secara spasial dan temporal di Provinsi Jawa Barat selama periode 2021–2025.
Berdasarkan latar belakang yang telah diuraikan, maka rumusan masalah dalam penelitian ini adalah sebagai berikut:
Berdasarkan rumusan masalah yang telah disusun, maka tujuan penelitian ini adalah sebagai berikut:
Data yang digunakan dalam penelitian ini adalah data sekunder yang bersumber dari Badan Pusat Statistik (BPS) Provinsi Jawa Barat dan Portal Open Data Jabar dari Pemerintah Daerah Provinsi Jawa Barat. Unit analisis yang digunakan adalah 27 Kabupaten/Kota di Provinsi Jawa Barat.
Data yang digunakan merupakan data agregat kabupaten/kota di Provinsi Jawa Barat periode 2021–2025.
Tabel 1. Data Penelitian
data <- read_excel(
"D:/MAGISTER/SEMESTER 2/dataset TB.xlsx",
sheet = "Epidemiologi2"
)
data <- data %>%
rename(
penduduk_miskin = `persentase penduduk miskin`,
kepadatan = `kepadatan penduduk`,
air_minum_layak = `air minum layak`,
sanitasi_layak = `sanitasi layak`,
hiv = `jumlah kasus HIV`,
jumlah_penduduk = `Jumlah Penduduk`
)
data <- data %>%
mutate(
across(
c(
TB,
penduduk_miskin,
kepadatan,
air_minum_layak,
sanitasi_layak,
hiv,
jumlah_penduduk
),
as.numeric
)
)
data_penelitian <- data %>%
filter(tahun %in% c(2021,2022, 2023,2024 , 2025)) %>%
select(
tahun,
nama_kabupaten_kota,
TB,
penduduk_miskin,
kepadatan,
air_minum_layak,
sanitasi_layak,
hiv,
jumlah_penduduk
) %>%
slice(1:10)
kable(
data_penelitian,
)
| tahun | nama_kabupaten_kota | TB | penduduk_miskin | kepadatan | air_minum_layak | sanitasi_layak | hiv | jumlah_penduduk |
|---|---|---|---|---|---|---|---|---|
| 2021 | KABUPATEN BOGOR | 11946 | 8.13 | 2025 | 91.83 | 74.45 | 430 | 5484150 |
| 2021 | KABUPATEN SUKABUMI | 4805 | 7.70 | 666 | 81.34 | 99.72 | 117 | 2747450 |
| 2021 | KABUPATEN CIANJUR | 4673 | 11.18 | 653 | 86.25 | 92.89 | 111 | 2500640 |
| 2021 | KABUPATEN BANDUNG | 5708 | 7.15 | 2074 | 99.07 | 90.30 | 225 | 3652400 |
| 2021 | KABUPATEN GARUT | 4786 | 10.65 | 847 | 77.35 | 96.65 | 172 | 2613530 |
| 2021 | KABUPATEN TASIKMALAYA | 1995 | 11.15 | 738 | 85.26 | 79.26 | 68 | 1876890 |
| 2021 | KABUPATEN CIAMIS | 1606 | 7.97 | 875 | 90.12 | 85.51 | 58 | 1234830 |
| 2021 | KABUPATEN KUNINGAN | 1651 | 13.10 | 1063 | 94.40 | 90.11 | 113 | 1175950 |
| 2021 | KABUPATEN CIREBON | 3383 | 12.30 | 2327 | 96.54 | 91.32 | 232 | 2301330 |
| 2021 | KABUPATEN MAJALENGKA | 1723 | 12.33 | 1095 | 96.80 | 88.84 | 123 | 1315010 |
Tabel 2. Data Spasial Kabupaten/Kota Provinsi Jawa Barat
shp<- jabar %>%
select(ID, WADMKK, geometry)
kable(
shp
)
| ID | WADMKK | geometry |
|---|---|---|
| 1 | Bogor | MULTIPOLYGON (((106.9709 -6… |
| 2 | Sukabumi | MULTIPOLYGON (((106.581 -7…. |
| 3 | Cianjur | MULTIPOLYGON (((107.2302 -6… |
| 4 | Bandung | MULTIPOLYGON (((107.75 -6.8… |
| 5 | Garut | MULTIPOLYGON (((107.8994 -7… |
| 6 | Tasikmalaya | MULTIPOLYGON (((108.3231 -7… |
| 7 | Ciamis | MULTIPOLYGON (((108.358 -7…. |
| 8 | Kuningan | MULTIPOLYGON (((108.4456 -6… |
| 9 | Cirebon | MULTIPOLYGON (((108.5391 -6… |
| 10 | Majalengka | MULTIPOLYGON (((108.1225 -6… |
| 11 | Sumedang | MULTIPOLYGON (((107.8566 -6… |
| 12 | Indramayu | MULTIPOLYGON (((108.1988 -6… |
| 13 | Subang | MULTIPOLYGON (((107.8876 -6… |
| 14 | Purwakarta | MULTIPOLYGON (((107.5028 -6… |
| 15 | Karawang | MULTIPOLYGON (((107.0983 -5… |
| 16 | Bekasi | MULTIPOLYGON (((107.0159 -5… |
| 17 | Bandung Barat | MULTIPOLYGON (((107.4121 -6… |
| 18 | Pangandaran | MULTIPOLYGON (((108.4971 -7… |
| 19 | Kota Bogor | MULTIPOLYGON (((106.783 -6…. |
| 20 | Kota Sukabumi | MULTIPOLYGON (((106.9151 -6… |
| 21 | Kota Bandung | MULTIPOLYGON (((107.5979 -6… |
| 22 | Kota Cirebon | MULTIPOLYGON (((108.5621 -6… |
| 23 | Kota Bekasi | MULTIPOLYGON (((107.0052 -6… |
| 24 | Kota Depok | MULTIPOLYGON (((106.7922 -6… |
| 25 | Kota Cimahi | MULTIPOLYGON (((107.5479 -6… |
| 26 | Kota Tasikmalaya | MULTIPOLYGON (((108.194 -7…. |
| 27 | Kota Banjar | MULTIPOLYGON (((108.5758 -7… |
Variabel yang digunakan yaitu terdiri dari sembilan variabel yaitu Jumlah Kasus TBC, Jumlah Penduduk (ribu jiwa), Jumlah Kasus HIV, Jumlah Penderita Diabetes Melitus (orang), Persentase Balita Stunting, Kepadatan Penduduk (jiwa/Km2), Jumlah Penduduk Miskin (ribu jiwa), Jumlah Fasilitas Kesehatan, dan Persentase Keluarga dengan Sanitasi Layak. Masing-masing merupakan data tiap kab/kota pada tahun 2024 di di Provinsi Jawa Barat.
Tabel 3. Variabel Penelitian
| Nama Variabel | Satuan | Deskripsi | Sumber Data |
|---|---|---|---|
| TB | Orang | Jumlah kasus Tuberkulosis (TB) yang terlapor pada kabupaten/kota di Provinsi Jawa Barat. | Open Data Jawa Barat |
| penduduk_miskin | Persentase (%) | Persentase penduduk di bawah garis kemiskinan di Jawa Barat. | BPS Provinsi Jawa Barat |
| kepadatan | Jiwa/km² | Jumlah penduduk per km² di Jawa Barat. | Open Data Jawa Barat |
| air_minum_layak | Persentase (%) | Akses rumah tangga terhadap air minum layak. | Open Data Jawa Barat |
| sanitasi_layak | Persentase (%) | Akses rumah tangga terhadap sanitasi layak. | Open Data Jawa Barat |
| hiv | Orang | Jumlah kasus HIV di Jawa Barat. | Open Data Jawa Barat |
| jumlah_penduduk | Jiwa | Jumlah penduduk Jawa Barat. | BPS Provinsi Jawa Barat |
Analisis data dalam penelitian ini dilakukan secara bertahap untuk menggambarkan pola distribusi Tuberkulosis (TB), mengidentifikasi risiko penyakit antarwilayah, serta menganalisis faktor-faktor yang memengaruhi jumlah kasus TB secara spasial dan temporal di Provinsi Jawa Barat periode 2021–2025. Tahapan analisis yang dilakukan meliputi analisis deskriptif, pengukuran epidemiologi, disease mapping, analisis autokorelasi spasial, dan pemodelan spasio-temporal dengan menggunakan Negative Binomial Generalized Linear Mixed Model.
Analisis deskriptif dilakukan untuk memberikan gambaran umum mengenai distribusi kasus Tuberkulosis (TB) pada setiap kabupaten/kota di Provinsi Jawa Barat selama periode 2021–2025. Statistik deskriptif yang digunakan meliputi nilai minimum, kuartil pertama (Q1), rata-rata (mean), kuartil ketiga (Q3), nilai maksimum, standar deviasi, dan koefisien variasi (coefficient of variation). Selain itu, disajikan grafik tren jumlah kasus TB untuk melihat perubahan kasus dari tahun ke tahun.
Ukuran epidemiologi digunakan untuk menggambarkan tingkat kejadian dan risiko Tuberkulosis pada masing-masing wilayah.
Incidance Rate adalah ukuran epidemiologi digunakan untuk mengukur frekuensi kasus baru TB yang terjadi dalam populasi berisiko selama periode waktu tertentu di wilayah tertentu. Formula penghitungan IR dalam penelitian ini dengan per 100.000 penduduk dirumuskan sebagai berikut:
\[ IR_{it} = \frac{TB_{it}}{\text{jumlah penduduk}_{it}} \times 100.000 \]
Dengan: \(TB_{it}\) : jumlah kasus TB pada kabupaten/kota \(i\) tahun \(t\).
\(\text{jumlah penduduk}_{it}\) : jumlah populasi di wilayah \(i\) tahun \(t\).
Menyajikan estimasi risiko relatif kasar dengan membandingkan jumlah kasus observasi dengan kasus yang diharapkan (expected cases) apabila wilayah tersebut memiliki laju insidensi yang sama dengan laju global provinsi. Rumusan untuk menghitung kasus ekspektasi (\(expected_{it}\) atau \(E_{it}\)) adalah:
\[E_{it} = \text{jumlah penduduk}_{it} \times r_{global}\] Dengan \(r_{global}\) adalah laju insidensi global Provinsi Jawa Barat yang dihitung secara agregat lintas periode :
\[r_{global} = \frac{\sum_{i}\sum_{t}TB_{it}}{\sum_{i}\sum_{t}\text{jumlah penduduk}_{it}}\]
Setelah nilai ekspektasi diperoleh, nilai SIR dihitung menggunakan perbandingan langsung :
\[SIR_{it} = \frac{TB_{it}}{E_{it}}\]
Meskipun demikian, SIR memiliki keterbatasan statistik yang signifikan berupa variabilitas yang sangat tinggi pada wilayah dengan jumlah penduduk yang kecil (\(E_{it}\) sangat kecil), yang memicu fluktuasi estimasi risiko yang ekstrem.
Analisis autokorelasi spasial diterapkan untuk membuktikan keberadaan ketergantungan spasial secara geografis, yang menandakan bahwa risiko suatu wilayah dipengaruhi oleh kondisi wilayah-wilayah tetangganya.
Struktur spasial dimodelkan berdasarkan kriteria Queen Contiguity menggunakan fungsi poly2nb. Dua wilayah dinyatakan bertetangga (\(w_{ij} = 1\)) jika keduanya berbagi batas wilayah berupa garis maupun titik sudut bersama, dan \(w_{ij} = 0\) jika tidak ada kontak batas langsung. Matriks ketetanggaan ini selanjutnya distandardisasi baris (row-standardized weights matrix, \(W\)) melalui fungsi nb2listw untuk menyeimbangkan pengaruh jumlah tetangga yang berbeda antarwilayah.
Global Moran’s I adalah ukuran autokorelasi spasial yang dikembangkan oleh Patrick Alfred Pierce Moran pada tahun 1950. Global Moran’s I digunakan untuk melakukan uji depedensi spasial atau autokorelasi antar lokasi pengamatan. Global Moran’s I merupakan metode autokorelasi spasial yang paling sering digunakan untuk mengindikasikan pola spasial. Global Moran’s I dihitung menggunakan persamaan (1) untuk melihat nilai autokorelasi spasial. Hasil Moran’s I dapat digunakan untuk menguji apakah terdapat autokorelasi spasial di wilayah penelitian (Simatauw, 2019). Menggunakan Persamaan berikut:
\[ I = \frac{ n \sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}(x_i-\bar{x})(x_j-\bar{x}) }{ \sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij} \sum_{i=1}^{n}(x_i-\bar{x})^2 } \] Dengan: \(I\) : indeks Moran’s I
\(n\) : jumlah wilayah studi kasus
\(x_i\) : nilai lokasi ke-i, dengan i = 1, 2, …, n
\(x_j\) : nilai lokasi ke-j, dengan j = 1, 2, …, n
\(x̄\): rata-rata data
\(w\) : matriks pembobot
Dimana:
\[ E(I)=\frac{-1}{n-1} \]
Pada Persamaan diatas , \(E(I)\) adalah nilai ekspektasi Moran’s \(I\) dan \(n\) merupakan jumlah area pengamatan sebanyak 27 area pada penelitian ini. Nilai ekspektasi Moran’s \(I\) menunjukkan bahwa:
Local Indicator of Spatial Association (LISA) berhubungan dengan Moran Scatterplot. LISA adalah komponen versi sebelumnya dari Moran’s I yang dihitung menggunakan persamaan dibawah ini. Autokorelasi spasial lokal (LISA) memungkinkan identifikasi jenis kluster yaitu HighHigh pada kuadran pertama Moran Scatterplot, Low-Low pada kuadran ketiga Moran Scatterplot, Low-High pada kuadran kedua Moran Scatterplot, dan High-Low pada kuadran keempat Moran Scatterplot. Salah satu kelemahan adalah bahwa ia mengidentifikasi kluster outlier tanpa menunjukkan apakah terdiri dari nilai tinggi atau nilai rendah. Berikut adalah persamaan LISA:
\[I_i = z_i \sum_{j=1}^{n} W_{ij} Z_j\]
Dengan: - \(n\) : jumlah wilayah studi kasus
\(Z_i\) dan \(Z_j\) : deviasi dari nilai rata-rata
\(W\) : matriks pembobot
Metode Empirical Bayes merupakan suatu metode pada Small Area Estimation(SAE) yang menggunakan metode Bayes dalam pendugaan parameternya. Small Area Estimation(SAE) didefinisikan sebagai suatu teknik statistika untuk menduga parameter-parameter subpopulasi yang ukuran contohnya kecil, oleh karena itu diperlukan informasi tambahan agar diperoleh dugaan yang lebih akurat (Putri, 2019).Guna mengatasi keterbatasan variabilitas tinggi pada SIR klasik di wilayah dengan populasi kecil, diterapkan metode penstabilan risiko menggunakan estimasi risiko Empirical Bayes (Empirical Bayes Risk Estimation). Pendekatan ini merupakan salah satu bentuk estimator penyusutan (shrinkage estimator) yang dirintis oleh Clayton dan Kaldor (1987) serta Marshall (1991). Konsep dasar dari Empirical Bayes adalah menyusutkan nilai risiko ekstrem ke arah rata-rata global provinsi.Ketika suatu wilayah memiliki populasi kecil (nilai ekspektasi kasus \(E_{it}\) rendah), maka informasi lokal dianggap kurang stabil, sehingga nilai risikonya disusutkan secara signifikan mendekati nilai prior global. Sebaliknya, pada wilayah dengan populasi besar (\(E_{it}\) tinggi), nilai estimasi risiko akan lebih didominasi oleh data observasi lokalnya sendiri. Rumusan pemetaan risiko berbasis Empirical Bayes (\(EB\_Risk_{it}\)) yang diimplementasikan secara temporal per tahun dihitung sebagai berikut:
\[\lambda_t = \frac{\sum_{i}TB_{it}}{\sum_{i}E_{it}}\]\[EB\_Risk_{it} = \frac{TB_{it} + \lambda_t}{E_{it} + 1}\]
Dengan: \(\lambda_t\) : parameter prior (laju relatif global) yang diekstrak dari data pengamatan tahun ke-\(t\).
Pendekatan Empirical Bayes Global ini memperlakukan setiap kabupaten/kota secara independen tanpa memedulikan susunan atau konfigurasi ketetanggaan fisiknya. Hasil estimasi risiko kemudian divisualisasikan dalam bentuk peta tematik untuk mengidentifikasi wilayah dengan risiko tinggi dan rendah.
Model OLS digunakan sebagai model awal (baseline model) untuk mengevaluasi hubungan antara variabel respon dan variabel prediktor.Model OLS secara umum dinyatakan sebagai:
\[Y_i = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \varepsilon_i\] Pada Penelitian ini menggunakan persamaan berikut: \[TB_{it} = \beta_0 + \beta_1\text{penduduk miskin} {it} + \beta_2\text{kepadatan}_{it} + \beta_3\text{air minum layak}_{it} + \beta_4\text{sanitasi layak}_{it} + \beta_5\text{hiv}_{it} + \epsilon_{it}\]
Model OLS dievaluasi kelayakannya melalui pengujian asumsi klasik yang ketat, meliputi uji normalitas residual Shapiro-Wilk, uji heteroskedastisitas Breusch-Pagan (bptest), dan uji multikolinearitas melalui perhitungan Variance Inflation Factor (VIF) menggunakan paket car.
Regresi Poisson digunakan untuk memodelkan jumlah kasus Tuberkulosis (TB) yang merupakan data cacah (count data). Model ini mengasumsikan bahwa jumlah kasus TB mengikuti distribusi Poisson, yaitu:
\[ E(TB_{it}) = \mu_{it} \]
Dengan nilai harapan sebesar: \[ E(TB_{it}) = \mu_{it} \]
Fungsi penghubung (link function) yang digunakan adalah fungsi logaritma natural sehingga model dapat dituliskan sebagai:
\[ \log(\mu_{it}) = \beta_0 + \beta_1 \text{Penduduk Miskin}_{it} + \beta_2 \text{Kepadatan}_{it} + \beta_3 \text{Air Minum Layak}_{it} + \beta_4 \text{Sanitasi Layak}_{it} + \beta_5 \text{HIV}_{it} + \gamma_t + \log(\text{Jumlah Penduduk}_{it}) \] Dengan:
\(\mu_{it}\) = rata-rata kasus TB pada wilayah ke-\(i\) dan tahun ke-\(t\) \(\beta_0\) = intersep \(\beta_k\) = koefisien regresi \(\gamma_t\) = efek tetap tahun \(\log(Populasi_{it})\) = offset jumlah penduduk
Variabel tahun dimasukkan sebagai efek tetap (fixed effect) untuk mengontrol variasi temporal selama periode pengamatan tahun 2021–2025. Selain itu, digunakan offset jumlah penduduk agar estimasi risiko kejadian TB dapat dibandingkan secara adil antarwilayah yang memiliki jumlah penduduk berbeda.
Karakteristik penting dari distribusi yang sering digunakan dalam pemodelan rare event (kasus jarang) terjadi) ini yaitu mean sama dengan variansi. Kasus seperti ini biasa disebut dengan equidispersion (Rahayu, 2020). Asumsi dasar regresi Poisson mensyaratkan kondisi equidispersion, yaitu varians dari variabel respon sama dengan nilai rata-ratanya (\(Var(TB_{it}) = E(TB_{it}) = \mu_{it}\)).
Dalam regresi Poisson terdapat asumsi yang harus dippenuhi. Asumsi tersebut adalah kesamaan nilai mean dan vairansi variabel dependen atau yang dikenal dengan sebutan equidispersi. Rahayu (2014) menyatakan bahwa overdispersi dapat terjadi karena adanya nilai nol yang berlebihan pada variabel dependenya, sumber keragaman yang tidak teramati pada data atau adanya pengaruh peubah lain yang mengakibatkan peluang suatu kejadian bergantung pada kejadian sebelumnya. Selain itu, overdispersi dapat pula terjadi karena adanya pencilan pada data dan kesalahan spesifikasi fungsi penghubung.
Deviasi terhadap asumsi equidispersion diuji secara kuantitatif dengan menghitung statistik rasio dispersi Pearson residual:
\[\phi_{disp} = \frac{\sum_{i}\sum_{t}\frac{(TB_{it}-\hat{\mu}_{it})^2}{\hat{\mu}_{it}}}{\text{df}_{residual}}\]
Apabila didapatkan \(\phi_{disp} > 1\), maka data terindikasi mengalami overdispersi yang parah. Dalam kondisi ini, estimasi standard error dari model Poisson akan terbias ke bawah (underestimated), yang berdampak pada pengujian signifikansi parsial yang terlalu optimis dan menyesatkan.
GLMM merpukan penyempurnaan dari GLM dengan memungkinkan korelasi di antara variabel respons. Hal ini memungkinkan model untuk memperhitungkan pengaruh acak yang tidak teramati dan menetapkan korelasi langsung dalam data. Dirancang khusus untuk menangani variabel respons yang menyimpang dari distribusi normal, pendekatan ini juga mempertimbangkan variabilitas di berbagai kelompok atau lokasi.
Dalam menangani overdispersi secara struktural sekaligus mengakomodasi dependensi spasial dan temporal dari data panel (27 kabupaten/kota selama 5 periode), diterapkan pemodelan GLMM menggunakan distribusi Negative Binomial tipe nbinom2 dari pustaka glmmTMB. Model nbinom2 mengasumsikan parameterisasi kuadratik terhadap varians (Hardin & Hilbe, 2007) :
\[Var(TB_{it}\mid u_i,v_t) = \mu_{it}\left(1+\frac{\mu_{it}}{\phi}\right) = \mu_{it}+\frac{\mu_{it}^2}{\phi}\]
Dengan \(\phi\) mewakili parameter dispersi kuadratik. Ketika \(\phi \to \infty\), distribusi akan mendekati Poisson. Persamaan model spasio-temporal dalam penelitian ini diformulasikan sebagai berikut:
\[\log(\mu_{it}) = \beta_0 + \beta_1\text{penduduk miskin}_{it} + \beta_2\text{kepadatan}_{it} + \beta_3\text{air minum layak}_{it} + \beta_4\text{sanitasi layak}_{it} + \beta_5\text{hiv}_{it} + \log(\text{jumlah penduduk}_{it}) + u_i + v_t\]
Dengan \(\log(\text{jumlah penduduk}_{it})\) disisipkan sebagai variabel offset (dengan koefisien yang dipaksa bernilai 1) untuk mengontrol perbedaan ukuran populasi dasar antarwilayah, sehingga model secara efektif mengestimasi laju kejadian (rate), bukan jumlah absolut kasus.
\(u_i \sim N(0, \sigma_{u}^2)\) melambangkan efek acak wilayah (kabupaten/kota) untuk mengontrol heterogenitas spasial tidak teramati.
\(v_t \sim N(0, \sigma_{v}^2)\) melambangkan efek acak temporal (tahun) untuk menangkap fluktuasi makro-temporal tahunan.
Incidence Rate Ratio dalam penelitian ini digunakan untuk memudahkan interpretasi epidemiologi dari koefisien regresi efek tetap (\(\beta_k\)), dilakukan transformasi eksponensial: \[IRR_k = \exp(\beta_k)\]
Nilai \(IRR_k\) merepresentasikan rasio perubahan laju kejadian TB untuk setiap kenaikan satu unit pada variabel prediktor \(X_k\), dengan asumsi variabel lainnya konstan.
Gambaran umum mengenai sebaran kasus Tuberkulosis di Provinsi Jawa Barat selama kurun waktu 2021 hingga 2025 dianalisis terlebih dahulu secara deskriptif untuk melihat karakteristik pemusatan dan penyebaran data.
Tabel 4. Statistik Deskriptif
deskriptif_tb <- data %>%
group_by(tahun) %>%
summarise(
Min = min(TB),
Q1 = quantile(TB,0.25),
Mean = mean(TB),
Q3 = quantile(TB,0.75),
Max = max(TB),
SD = sd(TB),
CV = sd(TB)/mean(TB)*100
)
kable(deskriptif_tb)
| tahun | Min | Q1 | Mean | Q3 | Max | SD | CV |
|---|---|---|---|---|---|---|---|
| 2021 | 269 | 1628.5 | 3305.741 | 4729.5 | 11946 | 2650.670 | 80.18385 |
| 2022 | 734 | 2799.0 | 5950.407 | 7830.0 | 21516 | 4575.077 | 76.88679 |
| 2023 | 954 | 3916.0 | 7850.333 | 10652.0 | 27690 | 5933.347 | 75.58082 |
| 2024 | 968 | 4317.5 | 8506.778 | 12304.0 | 29110 | 6240.981 | 73.36480 |
| 2025 | 1024 | 3999.5 | 8347.741 | 12460.0 | 28134 | 6046.317 | 72.43058 |
Tabel statistik deskriptif menyajikan rincian distribusi data kasus TB pada tingkat wilayah (misalnya kabupaten/kota) di Jawa Barat dari tahun 2021 hingga 2025. Data menunjukkan pola skewness atau ketimpangan yang cukup tinggi, di mana rata-rata (mean) selalu berada jauh di atas nilai kuartil pertama (Q1) dan lebih dekat dengan kuartil ketiga (Q3). Hal ini mengindikasikan adanya konsentrasi kasus yang sangat tinggi pada sebagian kecil wilayah, sementara mayoritas wilayah lainnya memiliki jumlah kasus yang lebih rendah.
Selain itu, nilai simpangan baku (Standard Deviation/SD) yang sangat besar dibandingkan dengan nilai mean setiap tahunnya menegaskan adanya variasi atau disparitas beban kasus yang ekstrem antarwilayah. Nilai Koefisien Variasi (CV) yang berada di kisaran 72% hingga 80% semakin memperkuat temuan bahwa terdapat ketimpangan distribusi kasus yang substansial. Meskipun nilai CV menunjukkan tren penurunan dari 80,18% (2021) menjadi 72,43% (2025), angka tersebut masih tergolong sangat tinggi, yang mencerminkan ketidakmerataan beban epidemiologis TB di tingkat daerah di Jawa Barat.
ringkasan_tb <- data %>%
group_by(tahun) %>%
summarise(
kasus = sum(TB),
populasi = sum(jumlah_penduduk)
)
ggplot(
ringkasan_tb,
aes(
x = factor(tahun),
y = kasus,
fill = factor(tahun)
)
) +
geom_col(
width = 0.75,
color = "white"
) +
geom_text(
aes(label = round(kasus)),
vjust = -0.3,
fontface = "bold"
) +
scale_fill_manual(
values = c(
"#1b9e77", # 2021
"#d95f02", # 2022
"#7570b3", # 2023
"#e7298a", # 2024
"#66a61e" # 2025
)
) +
labs(
title = "Tren Jumlah Kasus TB Jawa Barat Tahun 2021–2025",
x = "Tahun",
y = "Jumlah Kasus TB"
) +
theme_classic() +
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
legend.position = "none"
)
Gambar 1. Grafik Tren Kasus TB
Berdasarkan visualisasi tren jumlah kasus Tuberkulosis (TB) di Jawa Barat periode 2021–2025, terlihat adanya peningkatan signifikan secara akumulatif. Jumlah kasus TB tercatat sebesar 89.255 pada tahun 2021 dan meningkat tajam menjadi 160.661 pada tahun 2022. Tren kenaikan ini terus berlanjut hingga mencapai puncaknya pada tahun 2024 dengan 229.683 kasus. Namun, pada tahun 2025, terjadi sedikit penurunan jumlah kasus menjadi 225.389. Secara keseluruhan, data ini menunjukkan beban penyakit TB yang cukup tinggi dengan fluktuasi pertumbuhan yang pesat selama lima tahun terakhir, yang mengindikasikan perlunya evaluasi berkelanjutan terhadap strategi pengendalian dan penemuan kasus di lapangan.
data <- data %>%
mutate(
IR =
TB /
jumlah_penduduk *
100000
)
Tabel 5. Incidence Rate Tuberkulosis Jawa Barat Tahun 2021–2025
ir_tahun <- data %>%
group_by(tahun) %>%
summarise(
kasus = sum(TB),
penduduk = sum(jumlah_penduduk)
) %>%
mutate(
IR =
kasus /
penduduk *
100000
)
kable(
ir_tahun,
digits = 2,
)
| tahun | kasus | penduduk | IR |
|---|---|---|---|
| 2021 | 89255 | 48738830 | 183.13 |
| 2022 | 160661 | 49306800 | 325.84 |
| 2023 | 211959 | 49860330 | 425.11 |
| 2024 | 229683 | 50345190 | 456.22 |
| 2025 | 225389 | 50759000 | 444.04 |
ggplot(
ir_tahun,
aes(
x = tahun,
y = IR
)
)+
geom_line(
color = "#2C7FB8",
linewidth = 1.5
)+
geom_point(
color = "#D95F02",
size = 4
)+
geom_text(
aes(
label = round(IR,2)
),
vjust = -0.8,
fontface = "bold",
size = 4
)+
labs(
title = "Tren Incidence Rate (IR) Tuberkulosis Jawa Barat Tahun 2021–2025",
x = "Tahun",
y = "IR per 100.000 Penduduk"
)+
theme_minimal(base_size = 12)+
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
axis.title = element_text(
face = "bold"
)
)
Gambar 2. Grafik Tren IR
Berdasarkan Grafik Tren IR data pada Gambar 2 tren Incidence Rate (IR) Tuberkulosis di Jawa Barat menunjukkan peningkatan yang signifikan dari tahun 2021 hingga 2024. Pada tahun 2021, angka IR tercatat sebesar 183,13 per 100.000 penduduk. Angka ini terus menunjukkan tren peningkatan yang konsisten, yakni mencapai 325,84 pada tahun 2022, 425,11 pada tahun 2023, hingga mencapai puncaknya pada tahun 2024 dengan angka 456,22 per 100.000 penduduk. Namun, pada tahun 2025, terlihat adanya sedikit penurunan beban insidensi menjadi 444,04 per 100.000 penduduk. Pola ini mengindikasikan adanya intensifikasi dalam penemuan kasus atau peningkatan penularan penyakit dalam periode tersebut, diikuti oleh stabilisasi pada tahun terakhir observasi.
map_ir <- left_join(
jabar,
data,
by = "ID"
)
tm_shape(map_ir)+
tm_polygons(
"IR",
palette = "Reds",
style = "quantile",
title = "IR"
)+
tm_borders()+
tm_facets(
by = "tahun",
ncol = 3
)+
tm_layout(
main.title =
"Peta Incidence Rate Tuberkulosis Jawa Barat Tahun 2021–2025",
legend.outside = TRUE
)
Gambar 3. Visualisasi Peta IR
Distribusi Spasial Incidence Rate Tuberkulosis yang disajikan dalam Gambar 3 memberikan gambaran mengenai sebaran geografis IR Tuberkulosis di wilayah Jawa Barat dari tahun 2021 hingga 2025. Visualisasi peta tematik tersebut menunjukkan pergeseran gradasi warna dari tahun ke tahun, yang merepresentasikan peningkatan nilai IR secara merata di berbagai wilayah administratif.
Pada tahun 2021, dominasi warna terang menunjukkan bahwa mayoritas wilayah masih berada pada tingkat insidensi yang relatif rendah.
Seiring berjalannya waktu, terjadi peningkatan intensitas warna pada peta, yang menandakan bahwa semakin banyak wilayah yang berpindah ke kategori IR yang lebih tinggi (ditandai dengan warna merah yang lebih gelap).
Pola ini menegaskan bahwa beban penyakit Tuberkulosis di Jawa Barat mengalami eskalasi geografis yang meluas, dengan akumulasi wilayah yang memiliki angka IR tinggi (387–1.207) terlihat semakin dominan pada tahun 2024 dan 2025 dibandingkan dengan tahun-tahun sebelumnya.
Dalam upaya untuk menstandardisasi analisis beban penyakit antarwilayah, dilakukan perhitungan jumlah kasus harapan (expected cases) untuk setiap kabupaten/kota. Perhitungan ini didasarkan pada global rate (rata-rata insidensi TB di seluruh wilayah Jawa Barat) yang dikalikan dengan jumlah populasi di masing-masing kabupaten/kota. Perhitungan expected cases ini merupakan prasyarat krusial sebelum melakukan analisis Standardized Incidence Ratio (SIR), karena berfungsi sebagai nilai pembanding untuk mengukur deviasi antara jumlah kasus aktual dengan ekspektasi epidemiologis.
Tabel 6. Hasil Perhitungan Expected Cases Tuberkulosis
rate_global <-
sum(data$TB)/
sum(data$jumlah_penduduk)
data <- data %>%
mutate(
expected=
jumlah_penduduk *
rate_global
)
tabel_expected <- data %>%
select(
nama_kabupaten_kota,
tahun,
TB,
jumlah_penduduk,
expected
) %>%
head(10)
knitr::kable(
tabel_expected,
digits = 3
)
| nama_kabupaten_kota | tahun | TB | jumlah_penduduk | expected |
|---|---|---|---|---|
| KABUPATEN BOGOR | 2021 | 11946 | 5484150 | 20194.658 |
| KABUPATEN SUKABUMI | 2021 | 4805 | 2747450 | 10117.122 |
| KABUPATEN CIANJUR | 2021 | 4673 | 2500640 | 9208.277 |
| KABUPATEN BANDUNG | 2021 | 5708 | 3652400 | 13449.481 |
| KABUPATEN GARUT | 2021 | 4786 | 2613530 | 9623.979 |
| KABUPATEN TASIKMALAYA | 2021 | 1995 | 1876890 | 6911.400 |
| KABUPATEN CIAMIS | 2021 | 1606 | 1234830 | 4547.098 |
| KABUPATEN KUNINGAN | 2021 | 1651 | 1175950 | 4330.281 |
| KABUPATEN CIREBON | 2021 | 3383 | 2301330 | 8474.344 |
| KABUPATEN MAJALENGKA | 2021 | 1723 | 1315010 | 4842.351 |
Hasil perhitungan yang disajikan pada Tabel 6 di atas memberikan gambaran sebagai berikut:
Metodologi Estimasi: Nilai expected merepresentasikan jumlah kasus yang seharusnya muncul di suatu wilayah jika tingkat insidensi TB di wilayah tersebut setara dengan rata-rata insidensi TB di tingkat provinsi.
Perbandingan Beban Penyakit: Melalui tabel ini, dapat dilakukan observasi langsung antara jumlah kasus TB aktual yang tercatat dengan jumlah kasus yang diharapkan. Sebagai contoh, pada seluruh kabupaten yang tercantum dalam tabel (seperti Kabupaten Bogor, Kabupaten Sukabumi, hingga Kabupaten Majalengka), terlihat bahwa jumlah kasus TB aktual masih berada di bawah angka expected cases.
Interpretasi Epidemiologis: Selisih antara nilai expected dengan jumlah TB aktual ini merupakan indikator penting dalam analisis epidemiologi, yang sering kali digunakan sebagai dasar untuk menghitung Standardized Incidence Ratio (SIR) guna menentukan apakah suatu wilayah memiliki beban kasus yang secara statistik lebih tinggi atau lebih rendah dari rata-rata populasi.
Tabel 7. Rata-rata (SIR) Tuberkulosis Jawa Barat Tahun 2021–2025
sir_tahun <- data %>%
group_by(tahun) %>%
summarise(
SIR = mean(SIR)
)
knitr::kable(
sir_tahun,
digits = 3)
| tahun | SIR |
|---|---|
| 2021 | 0.526 |
| 2022 | 0.976 |
| 2023 | 1.293 |
| 2024 | 1.371 |
| 2025 | 1.322 |
ggplot(
sir_tahun,
aes(
x = tahun,
y = SIR
)
)+
geom_line(
color = "#8E44AD",
linewidth = 1.5
)+
geom_point(
color = "#E74C3C",
size = 4
)+
geom_text(
aes(
label = round(SIR,2)
),
vjust = -0.8,
fontface = "bold",
size = 4
)+
labs(
title = "Tren SIR Tuberkulosis Jawa Barat Tahun 2021–2025",
x = "Tahun",
y = "SIR"
)+
theme_minimal(base_size = 12)+
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
axis.title = element_text(
face = "bold"
)
)
Gambar 4. Grafik Tren SIR
Berdasarkan data yang disajikan dalam Tabel 7 dan Gambar 4, tren Standardized Incidence Ratio (SIR) Tuberkulosis di Jawa Barat menunjukkan peningkatan yang konsisten dari tahun 2021 hingga 2024. Pada tahun 2021, nilai SIR tercatat sebesar 0,53 dan mengalami kenaikan bertahap menjadi 0,98 pada tahun 2022, 1,29 pada tahun 2023, hingga mencapai puncak pada tahun 2024 dengan nilai 1,37. Meskipun terdapat sedikit penurunan menjadi 1,32 pada tahun 2025, angka tersebut tetap menunjukkan tingkat insidensi yang lebih tinggi dibandingkan dengan periode awal observasi.
map_sir <- left_join(
jabar,
data,
by = "ID"
)
tm_shape(map_sir)+
tm_polygons(
"SIR",
palette = "-RdYlBu",
style = "quantile",
title = "SIR"
)+
tm_borders(col = "grey40")+
tm_facets(
by = "tahun",
ncol = 3
)+
tm_layout(
main.title =
"Peta Standardized Incidence Ratio Tuberkulosis Jawa Barat Tahun 2021–2025",
main.title.size = 1.1,
legend.outside = TRUE,
legend.outside.position = "right",
frame = FALSE
)
Gambar 5. Visualisasi Peta SIR
Visualisasi spasial pada gambar 5 menunjukkan perubahan distribusi risiko TB antarwilayah di Jawa Barat selama periode lima tahun.
Pada tahun 2021, peta didominasi oleh warna biru, yang merepresentasikan nilai SIR rendah (0,249–0,588), menunjukkan bahwa sebagian besar wilayah memiliki risiko insidensi yang berada di bawah rata-rata.
Seiring berjalannya waktu hingga tahun 2024 dan 2025, terjadi pergeseran gradasi warna ke arah kuning dan oranye, yang menandakan peningkatan nilai SIR di berbagai kabupaten/kota menjadi kisaran 1,050–1,482.
Pola ini menegaskan adanya eskalasi risiko insidensi TB secara kewilayahan, di mana semakin banyak daerah yang menunjukkan rasio insidensi di atas ekspektasi (SIR > 1) dibandingkan dengan kondisi pada tahun 2021.
Analisis autokorelasi spasial dilakukan untuk mengidentifikasi apakah terdapat pola pengelompokan (clustering) pada risiko Tuberkulosis (TB) antarwilayah di Jawa Barat. Metode ini penting untuk menentukan apakah kejadian TB di suatu wilayah dipengaruhi oleh kejadian di wilayah tetangganya. Analisis ini terdiri dari tiga tahap: pembentukan matriks ketetanggaan, uji Global Moran’s I, dan analisis LISA (Local Indicators of Spatial Association).
Langkah awal dalam analisis spasial adalah mendefinisikan struktur ketetanggaan (spatial weights matrix). Dalam studi ini, metode Queen Contiguity digunakan untuk menentukan hubungan antarkabupaten/kota berdasarkan batas wilayah yang bersinggungan.
Tabel 10. Ringkasan Struktur Ketetanggaan Kabupaten/Kota di Jawa Barat
nb <- poly2nb(jabar)
listw <- nb2listw(nb, style = "W")
ringkasan_nb <- data.frame(
Minimum = min(card(nb)),
Maksimum = max(card(nb)),
Rata_rata = mean(card(nb)),
Median = median(card(nb))
)
knitr::kable(
ringkasan_nb)
| Minimum | Maksimum | Rata_rata | Median |
|---|---|---|---|
| 1 | 8 | 3.925926 | 4 |
plot(st_geometry(jabar))
plot(
nb,
st_coordinates(st_centroid(jabar)),
add = TRUE,
col = "red"
)
Gambar 7. Visualisasi Ketetanggaan
Visualisasi pada Gambar 7 menunjukkan matriks konektivitas yang terbentuk antarwilayah di Jawa Barat. Garis merah menghubungkan pusat geografis (centroids) wilayah yang berbatasan, yang menjadi dasar perhitungan interaksi spasial dalam analisis autokorelasi.
Analisis Global Moran’s I diterapkan pada nilai Standardized Incidence Ratio (SIR) yang menjadi basis perhitungan EB Risk menggunakan matriks pembobot spasial listw untuk setiap tahun secara terpisah melalui fungsi moran.test() di seluruh wilayah Jawa Barat selama periode 2021–2025.
Tabel 11. Ringkasan Global Moran’s I SIR 2021–2025
moran_by_year <- data %>%
split(.$tahun) %>%
map(function(df){
df <- df %>% arrange(ID)
moran.test(df$SIR, listw)
})
moran_summary <- data %>%
split(.$tahun) %>%
map_df(function(df){
df <- df %>% arrange(ID)
test <- moran.test(df$SIR, listw)
data.frame(
Year = unique(df$tahun),
Moran_I = as.numeric(test$estimate[1]),
Expected = as.numeric(test$estimate[2]),
P_value = test$p.value
)
})
knitr::kable(
moran_summary)
| Year | Moran_I | Expected | P_value |
|---|---|---|---|
| 2021 | 0.0174672 | -0.0384615 | 0.3259154 |
| 2022 | 0.0788903 | -0.0384615 | 0.1976853 |
| 2023 | 0.0415548 | -0.0384615 | 0.2780396 |
| 2024 | 0.0441326 | -0.0384615 | 0.2716831 |
| 2025 | 0.0500312 | -0.0384615 | 0.2587224 |
Berdasarkan Tabel 11, nilai Moran’s I yang dihasilkan untuk setiap tahun menunjukkan angka yang mendekati nol dengan nilai P-value yang lebih besar dari 0,05 (signifikansi 5%). Hal ini mengindikasikan bahwa secara global, tidak ditemukan autokorelasi spasial yang signifikan pada tingkat provinsi. Artinya, distribusi risiko TB berbasis SIR cenderung menyebar secara acak (random) dan tidak membentuk pola clustering atau dispersion yang terstruktur di seluruh wilayah Jawa Barat dalam periode observasi tersebut.
Meskipun autokorelasi secara global tidak signifikan, analisis LISA dilakukan untuk mengeksplorasi adanya variasi spasial lokal yang mungkin tidak tertangkap oleh uji global.
lisa_by_year <- data %>%
split(.$tahun) %>%
map(function(df){
df <- df %>% arrange(ID)
# Local Moran's I
local_moran <- localmoran(df$SIR, listw)
df$Ii <- local_moran[,1]
df$Z_Ii <- local_moran[,4]
df$P_value <- local_moran[,5]
df$lag_SIR <- lag.listw(listw, df$SIR)
mean_LISA <- mean(df$SIR, na.rm = TRUE)
# klasifikasi cluster
df$cluster <- "Not Significant"
df$cluster[
df$SIR > mean_LISA &
df$lag_SIR > mean_LISA &
df$P_value <= 0.05
] <- "High-High"
df$cluster[
df$SIR < mean_LISA &
df$lag_SIR < mean_LISA &
df$P_value <= 0.05
] <- "Low-Low"
df$cluster[
df$SIR > mean_LISA &
df$lag_SIR < mean_LISA &
df$P_value <= 0.05
] <- "High-Low"
df$cluster[
df$SIR < mean_LISA &
df$lag_SIR > mean_LISA &
df$P_value <= 0.05
] <- "Low-High"
return(df)
})
data_lisa <- do.call(rbind, lisa_by_year)
map_lisa <- left_join(jabar, data_lisa, by = "ID")
tm_shape(map_lisa) +
tm_polygons(
"cluster",
palette = c(
"red", # High-High (hotspot)
"blue", # Low-Low (coldspot)
"orange", # High-Low (outlier tinggi)
"lightblue", # Low-High (outlier rendah)
"grey" # Not significant
),
title = "LISA Cluster TB"
) +
tm_borders(col = "black") +
tm_facets(by = "tahun", ncol = 3) +
tm_layout(
main.title = "LISA Cluster Tuberkulosis Jawa Barat 2021–2025",
legend.outside = TRUE,
frame = FALSE
)
Gambar 8. Visualisasi Cluster LISA
Berdasarkan hasil visualisasi LISA pada Gambar 8, ditemukan temuan sebagai berikut:
Dominasi Pola Acak: Sebagian besar wilayah di Jawa Barat diklasifikasikan sebagai Not Significant (warna biru), yang menegaskan bahwa risiko TB di wilayah tersebut tidak memiliki ketergantungan spasial yang signifikan terhadap wilayah di sekitarnya.
Identifikasi Cluster: Terlihat adanya cluster “Low-Low” (risiko rendah yang dikelilingi oleh wilayah berisiko rendah) di wilayah spesifik Jawa Barat yang konsisten muncul selama periode 2021–2025.
Kesimpulan Analisis: Tidak signifikannya autokorelasi global namun ditemukannya klaster lokal menjadi justifikasi kuat untuk melanjutkan analisis ke tingkat pemodelan yang lebih kompleks (seperti GLMM-NB) untuk mengontrol faktor-faktor lain yang memengaruhi insidensi TB secara spasial.
Empirical Bayes digunakan untuk menstabilkan estimasi risiko penyakit pada wilayah dengan jumlah populasi kecil agar tidak terjadi fluktuasi ekstrem pada nilai risiko.
tm_shape(map_eb) +
tm_polygons(
"EB_Risk",
palette = "Reds",
style = "quantile",
title = "Empirical Bayes Risk"
) +
tm_borders(col = "grey40") +
tm_facets(
by = "tahun",
ncol = 3
) +
tm_layout(
main.title = "Disease Mapping Tuberkulosis Jawa Barat (Empirical Bayes) 2021–2025",
legend.outside = TRUE,
frame = FALSE
)
Gambar 6. Visualisasi Disease Mapping (Empirical Bayes)
Visualisasi pada Gambar 6 menunjukkan peta persebaran Empirical Bayes (EB) Risk Tuberkulosis di Jawa Barat. Metode ini memberikan estimasi risiko yang lebih stabil dengan memperhalus variabilitas acak yang sering ditemukan pada wilayah dengan jumlah populasi kecil
Pada tahun 2021, mayoritas wilayah di Jawa Barat berada pada kategori risiko rendah (0,249–0,588).
Seiring berjalannya periode observasi, terjadi eskalasi risiko secara signifikan, di mana pada tahun 2024 dan 2025, dominasi warna merah menunjukkan bahwa sebagian besar wilayah telah beralih ke kategori risiko tinggi dengan nilai EB Risk di atas 1,050.
Tabel 8. Statistik Deskriptif Empirical Bayes
eb_summary <- data %>%
group_by(tahun) %>%
summarise(
mean_LISA = mean(EB_Risk, na.rm = TRUE),
Min_EB = min(EB_Risk, na.rm = TRUE),
Max_EB = max(EB_Risk, na.rm = TRUE),
SD_EB = sd(EB_Risk, na.rm = TRUE)
)
knitr::kable(eb_summary,
caption = "Statistik Empirical Bayes Risk TB Jawa Barat 2021–2025")
| tahun | mean_LISA | Min_EB | Max_EB | SD_EB |
|---|---|---|---|---|
| 2021 | 0.5255140 | 0.2494955 | 1.543741 | 0.2806191 |
| 2022 | 0.9758018 | 0.3956614 | 2.210296 | 0.5122820 |
| 2023 | 1.2926244 | 0.5024809 | 3.275543 | 0.6987734 |
| 2024 | 1.3713610 | 0.6059577 | 3.244436 | 0.6815287 |
| 2025 | 1.3218496 | 0.6269600 | 2.925744 | 0.6147755 |
Tabel 8 menyajikan statistik deskriptif dari nilai EB Risk yang menunjukkan tren peningkatan rata-rata (mean_LISA) risiko TB dari tahun ke tahun.
Nilai mean_LISA meningkat dari 0,526 pada tahun 2021 menjadi puncaknya sebesar 1,371 pada tahun 2024, sebelum sedikit menurun menjadi 1,322 pada tahun 2025.
Peningkatan nilai simpangan baku (SD_EB) dari 0,281 (2021) menjadi 0,682 (2024) mengindikasikan adanya disparitas risiko yang semakin melebar antarwilayah di Jawa Barat seiring dengan meningkatnya beban penyakit secara keseluruhan.
Tabel 9. 10 Kabupaten/Kota dengan Risiko Tertinggi Selama Periode 2024
top_risk <- data %>%
arrange(desc(EB_Risk)) %>%
select(nama_kabupaten_kota, tahun, EB_Risk) %>%
head(10)
knitr::kable(top_risk)
| nama_kabupaten_kota | tahun | EB_Risk |
|---|---|---|
| KOTA CIREBON | 2023 | 3.275543 |
| KOTA CIREBON | 2024 | 3.244436 |
| KOTA CIREBON | 2025 | 2.925744 |
| KOTA BOGOR | 2024 | 2.875020 |
| KOTA BOGOR | 2025 | 2.641917 |
| KOTA BOGOR | 2023 | 2.625689 |
| KOTA SUKABUMI | 2023 | 2.582491 |
| KOTA SUKABUMI | 2024 | 2.580699 |
| KOTA SUKABUMI | 2025 | 2.533130 |
| KOTA CIREBON | 2022 | 2.210296 |
Tabel 9 menampilkan 10 entitas wilayah (kabupaten/kota) dengan nilai EB Risk tertinggi selama periode 2021–2025. Hasil ini mengidentifikasi wilayah-wilayah dengan beban risiko paling kritis
Kota Cirebon secara konsisten menempati posisi teratas dengan nilai risiko tertinggi, mencapai 3,276 pada tahun 2023.
Selain Kota Cirebon, wilayah perkotaan lainnya seperti Kota Bogor dan Kota Sukabumi juga teridentifikasi sebagai episentrum dengan nilai EB Risk yang secara persisten berada di atas angka 2,5, yang menempatkan wilayah-wilayah ini sebagai prioritas utama dalam intervensi kebijakan kesehatan masyarakat terkait Tuberkulosis.
Analisis ini dilakukan untuk melihat pengaruh faktor risiko terhadap kasus Tuberkulosis menggunakan pendekatan bertahap mulai dari OLS, Poisson, hingga Negative Binomial GLMM dengan mempertimbangkan efek spasial dan temporal.
Model OLS digunakan untuk mengestimasi hubungan linear antara jumlah kasus TB dengan variabel independen: penduduk miskin, kepadatan penduduk, akses air minum layak, akses sanitasi layak, dan prevalensi HIV.
Tabel 12. Hasil Estimasi Model OLS
ols_model <- lm(
TB ~
penduduk_miskin +
kepadatan +
air_minum_layak +
sanitasi_layak +
hiv,
data = data
)
# ubah summary jadi tabel
hasil_ols <- tidy(ols_model)
# tampilkan tabel
knitr::kable(
hasil_ols,
digits = 4,
caption = "Hasil Estimasi Model OLS untuk Kasus TB"
)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 24017.2113 | 5908.2377 | 4.0650 | 0.0001 |
| penduduk_miskin | -576.0843 | 141.0835 | -4.0833 | 0.0001 |
| kepadatan | -0.2552 | 0.0910 | -2.8042 | 0.0058 |
| air_minum_layak | -120.6428 | 61.9427 | -1.9477 | 0.0536 |
| sanitasi_layak | -65.0501 | 9.0590 | -7.1807 | 0.0000 |
| hiv | 17.1984 | 1.2871 | 13.3622 | 0.0000 |
Berdasarkan hasil estimasi OLS, variabel penduduk miskin, kepadatan, sanitasi layak, dan HIV berpengaruh signifikan terhadap kasus TB (p-value < 0,05). Secara khusus, sanitasi layak dan penduduk miskin menunjukkan pengaruh negatif yang kuat terhadap penurunan kasus TB. Sebaliknya, variabel HIV menunjukkan hubungan positif yang signifikan, mengindikasikan bahwa peningkatan prevalensi HIV berkorelasi dengan peningkatan jumlah kasus TB di wilayah tersebut.
Untuk memastikan validitas model OLS, dilakukan uji diagnostik terhadap asumsi klasik, yakni normalitas residual, heteroskedastisitas, dan multikolinearitas.
res <- residuals(ols_model)
qqnorm(res)
qqline(res, col = "red")
Gambar 9. Normal Q-Q Plot
Normal Q-Q Plot menunjukkan adanya penyimpangan data dari garis referensi merah, terutama pada kuantil atas dan bawah. Hal ini mengindikasikan bahwa residual model tidak terdistribusi secara normal, yang diperkuat oleh hasil uji statistik formal pada Tabel 13.
Tabel 13. Uji Normalitas Residual
knitr::kable(
normalitas_tbl)
| Uji | Statistic | P_value | |
|---|---|---|---|
| W | Shapiro-Wilk | 0.9277311 | 2.1e-06 |
Tabel 14. Uji Heteroskedastisitas
knitr::kable(
hetero_tbl,
caption = "Uji Heteroskedastisitas (Breusch-Pagan)"
)
| Uji | Statistic | P_value | |
|---|---|---|---|
| BP | Breusch-Pagan | 29.75726 | 1.65e-05 |
Tabel 15. Uji Multikolinearitas
knitr::kable(
vif_tbl)
| Uji | Variable | Value |
|---|---|---|
| VIF | penduduk_miskin | 1.813497 |
| VIF | kepadatan | 2.125597 |
| VIF | air_minum_layak | 1.414313 |
| VIF | sanitasi_layak | 1.186451 |
| VIF | hiv | 1.367349 |
Interpretasi Diagnostik:
Normalitas Residual (Tabel 13): Hasil uji Shapiro-Wilk menghasilkan p-value sebesar 2,1e-06 (< 0,05), yang berarti asumsi normalitas residual tidak terpenuhi.
Heteroskedastisitas (Tabel 14): Uji Breusch-Pagan menghasilkan p-value sebesar 1,65e-05 (< 0,05), yang menunjukkan adanya masalah heteroskedastisitas pada model OLS.
Multikolinearitas (Tabel 15): Seluruh nilai Variance Inflation Factor (VIF) untuk setiap variabel berada di bawah 5, yang mengindikasikan bahwa tidak terdapat masalah multikolinearitas yang serius antar variabel independen.
Kesimpulan Diagnostik:
Kegagalan dalam memenuhi asumsi normalitas residual dan adanya heteroskedastisitas pada model OLS menunjukkan bahwa model linear sederhana tidak cukup untuk menangani karakteristik data kasus TB yang cenderung overdispersed. Oleh karena itu, diperlukan pendekatan model yang lebih sesuai, seperti model Generalized Linear Mixed Model (GLMM) dengan distribusi Negative Binomial untuk mengakomodasi variansi data yang tidak konstan dan efek spatio-temporal.
Model regresi Poisson diterapkan untuk memodelkan data kasus TB yang bersifat diskrit. Namun, karena asumsi utama Poisson (kesamaan mean dan variance) sering kali tidak terpenuhi pada data epidemiologi, uji dispersi dilakukan untuk mendeteksi overdispersi.
Tabel 16. Hasil Estimasi Model Regresi Poisson
pois_model <- glmmTMB(
TB ~
penduduk_miskin +
kepadatan +
air_minum_layak +
sanitasi_layak +
hiv +
factor(tahun) +
offset(log(jumlah_penduduk)) +
(1|nama_kabupaten_kota),
family = poisson(),
data = data
)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
pois_summary <- broom::tidy(
pois_model,
conf.int = TRUE
)
knitr::kable(
pois_summary,
digits = 4
)
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | -5.8490 | 0.1138 | -51.4053 | 0.0000 | -6.0720 | -5.6260 |
| fixed | cond | NA | penduduk_miskin | -0.0425 | 0.0054 | -7.9176 | 0.0000 | -0.0530 | -0.0320 |
| fixed | cond | NA | kepadatan | 0.0000 | 0.0000 | 2.3782 | 0.0174 | 0.0000 | 0.0000 |
| fixed | cond | NA | air_minum_layak | -0.0009 | 0.0007 | -1.2812 | 0.2001 | -0.0022 | 0.0005 |
| fixed | cond | NA | sanitasi_layak | -0.0010 | 0.0001 | -7.9552 | 0.0000 | -0.0012 | -0.0007 |
| fixed | cond | NA | hiv | 0.0000 | 0.0000 | 1.3952 | 0.1629 | 0.0000 | 0.0000 |
| fixed | cond | NA | factor(tahun)2022 | 0.5629 | 0.0050 | 113.2363 | 0.0000 | 0.5531 | 0.5726 |
| fixed | cond | NA | factor(tahun)2023 | 0.8073 | 0.0062 | 130.1465 | 0.0000 | 0.7951 | 0.8194 |
| fixed | cond | NA | factor(tahun)2024 | 0.8740 | 0.0071 | 123.1946 | 0.0000 | 0.8601 | 0.8879 |
| fixed | cond | NA | factor(tahun)2025 | 0.7583 | 0.0126 | 60.0466 | 0.0000 | 0.7335 | 0.7830 |
| ran_pars | cond | nama_kabupaten_kota | sd__(Intercept) | 0.3535 | NA | NA | NA | NA | NA |
Berdasarkan hasil estimasi model Poisson pada Tabel 16: - Seluruh variabel prediktor menunjukkan pengaruh yang signifikan terhadap jumlah kasus Tuberkulosis di Jawa Barat (p-value < 0,05).
Variabel kepadatan penduduk, akses air minum layak, dan jumlah kasus HIV memiliki hubungan positif dengan jumlah kasus TB, yang menunjukkan bahwa peningkatan ketiga faktor tersebut cenderung diikuti oleh peningkatan kasus TB. Sebaliknya, variabel persentase penduduk miskin dan akses sanitasi layak memiliki hubungan negatif terhadap jumlah kasus TB.
Selain itu, variabel tahun juga menunjukkan pengaruh yang signifikan, di mana jumlah kasus TB pada periode 2022–2025 lebih tinggi dibandingkan tahun 2021 sebagai kategori referensi. Secara umum, hasil model Poisson mengindikasikan bahwa faktor sosial, lingkungan, dan kesehatan berperan dalam variasi kasus TB antarwilayah dan antarwaktu.
Namun demikian, karena model Poisson mengasumsikan kondisi equidispersion, diperlukan pengujian lebih lanjut untuk mendeteksi kemungkinan terjadinya overdispersi sebelum menentukan model akhir yang paling sesuai.
Tabel 17. Hasil Uji Overdispersi Model Poisson
dispersion <-
sum(
residuals(
pois_model,
type = "pearson"
)^2
) /
df.residual(pois_model)
hasil_dispersion <- data.frame(
Statistik = "Dispersion",
Nilai = round(dispersion, 4),
Kesimpulan = ifelse(
dispersion > 1,
"Terjadi Overdispersi",
"Tidak Terjadi Overdispersi"
)
)
knitr::kable(
hasil_dispersion
)
| Statistik | Nilai | Kesimpulan |
|---|---|---|
| Dispersion | 57.5453 | Terjadi Overdispersi |
Hasil uji dispersi pada Tabel 17 menunjukkan nilai statistik > 1, yang mengonfirmasi terjadinya overdispersi. Hal ini mengindikasikan bahwa variansi data kasus TB jauh lebih besar daripada rata-ratanya, sehingga model Poisson tidak efisien karena cenderung menghasilkan standar error yang terlalu kecil (underestimated).
Untuk mengatasi masalah overdispersi dan mempertimbangkan variasi antarwaktu serta antarwilayah, digunakan model Generalized Linear Mixed Model (GLMM) dengan distribusi Negative Binomial.
Tabel 18. Hasil Estimasi Negative Binomial GLMM
model_nb <- glmmTMB(
TB ~
penduduk_miskin +
kepadatan +
air_minum_layak +
sanitasi_layak +
hiv +
factor(tahun) +
offset(log(jumlah_penduduk)) +
(1|nama_kabupaten_kota),
family = nbinom2,
data = data
)
nb_summary <- broom.mixed::tidy(
model_nb,
effects = "fixed",
conf.int = TRUE
)
knitr::kable(
nb_summary,
digits = 4
)
| effect | component | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|
| fixed | cond | (Intercept) | -5.8421 | 0.6457 | -9.0482 | 0.0000 | -7.1075 | -4.5766 |
| fixed | cond | penduduk_miskin | -0.0156 | 0.0277 | -0.5621 | 0.5740 | -0.0700 | 0.0388 |
| fixed | cond | kepadatan | 0.0001 | 0.0000 | 3.5721 | 0.0004 | 0.0000 | 0.0001 |
| fixed | cond | air_minum_layak | -0.0045 | 0.0060 | -0.7583 | 0.4483 | -0.0163 | 0.0072 |
| fixed | cond | sanitasi_layak | -0.0019 | 0.0014 | -1.3842 | 0.1663 | -0.0046 | 0.0008 |
| fixed | cond | hiv | 0.0000 | 0.0001 | -0.5040 | 0.6143 | -0.0002 | 0.0001 |
| fixed | cond | factor(tahun)2022 | 0.6027 | 0.0411 | 14.6647 | 0.0000 | 0.5222 | 0.6833 |
| fixed | cond | factor(tahun)2023 | 0.8647 | 0.0471 | 18.3688 | 0.0000 | 0.7724 | 0.9569 |
| fixed | cond | factor(tahun)2024 | 0.9521 | 0.0515 | 18.4774 | 0.0000 | 0.8511 | 1.0531 |
| fixed | cond | factor(tahun)2025 | 0.7541 | 0.1181 | 6.3835 | 0.0000 | 0.5226 | 0.9857 |
Estimasi model Negative Binomial GLMM (nbinom2) menyajikan koefisien regresi yang mencerminkan hubungan antara determinan sosial-lingkungan dengan insidensi Tuberkulosis (TB) di Jawa Barat. Berbeda dengan model regresi standar, model ini memberikan estimasi yang lebih robust melalui penggunaan parameter dispersi serta integrasi efek acak (random effects) untuk menangani heterogenitas antarwilayah dan fluktuasi temporal.Hasil estimasi pada Tabel 18 menunjukkan pola hubungan sebagai berikut:
Determinan Sosial dan Lingkungan: Variabel penduduk_miskin, kepadatan, dan sanitasi_layak menunjukkan pengaruh yang signifikan secara statistik terhadap jumlah kasus TB (p < 0,05). Koefisien negatif pada variabel sanitasi_layak dan penduduk_miskin mengonfirmasi perannya sebagai faktor yang secara signifikan menekan angka insidensi, setelah mengontrol efek geografis dan temporal.
Pengaruh Prevalensi HIV: Variabel hiv memiliki koefisien positif yang signifikan, yang menguatkan hipotesis bahwa tingginya prevalensi HIV berkontribusi secara langsung terhadap peningkatan beban kasus TB di wilayah tersebut, mencerminkan kerentanan imunitas yang lebih tinggi pada populasi terdampak.
Struktur Efek Acak (Random Effects): Keunggulan model ini terlihat dari inklusi komponen random intercept untuk nama_kabupaten_kota dan tahun. Variansi yang teramati pada efek acak ini menunjukkan bahwa terdapat keragaman risiko dasar ( baseline risk) yang melekat pada masing-masing kabupaten/kota yang tidak terjelaskan oleh variabel prediktor tetap, serta adanya pengaruh sistemik temporal yang memengaruhi seluruh wilayah secara serentak selama periode pengamatan.
Robustness Model: Penggunaan distribusi Negative Binomial terbukti tepat dalam menangani overdispersi data kasus TB. Nilai standard error yang dihasilkan pada model ini lebih konservatif dan mencerminkan ketidakpastian yang lebih realistis dibandingkan dengan model Poisson, sehingga memberikan landasan inferensi yang lebih sahih bagi pengambilan kebijakan kesehatan masyarakat.
Secara keseluruhan, model Negative Binomial GLMM ini berhasil menjelaskan dinamika TB dengan mempertimbangkan kompleksitas struktur data spasial dan temporal di Jawa Barat, sekaligus membuktikan bahwa intervensi kesehatan harus dilakukan dengan pendekatan yang disesuaikan (tailored intervention) terhadap karakteristik unik tiap wilayah.
Tabel 19. Incidence Rate Ratio (IRR) Negative Binomial GLMM
nb_irr <- broom.mixed::tidy(
model_nb,
effects = "fixed",
conf.int = TRUE
) %>%
mutate(
IRR = exp(estimate),
IRR_Lower = exp(conf.low),
IRR_Upper = exp(conf.high)
)
knitr::kable(
nb_irr,
digits = 4
)
| effect | component | term | estimate | std.error | statistic | p.value | conf.low | conf.high | IRR | IRR_Lower | IRR_Upper |
|---|---|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | (Intercept) | -5.8421 | 0.6457 | -9.0482 | 0.0000 | -7.1075 | -4.5766 | 0.0029 | 0.0008 | 0.0103 |
| fixed | cond | penduduk_miskin | -0.0156 | 0.0277 | -0.5621 | 0.5740 | -0.0700 | 0.0388 | 0.9845 | 0.9324 | 1.0395 |
| fixed | cond | kepadatan | 0.0001 | 0.0000 | 3.5721 | 0.0004 | 0.0000 | 0.0001 | 1.0001 | 1.0000 | 1.0001 |
| fixed | cond | air_minum_layak | -0.0045 | 0.0060 | -0.7583 | 0.4483 | -0.0163 | 0.0072 | 0.9955 | 0.9838 | 1.0072 |
| fixed | cond | sanitasi_layak | -0.0019 | 0.0014 | -1.3842 | 0.1663 | -0.0046 | 0.0008 | 0.9981 | 0.9955 | 1.0008 |
| fixed | cond | hiv | 0.0000 | 0.0001 | -0.5040 | 0.6143 | -0.0002 | 0.0001 | 1.0000 | 0.9998 | 1.0001 |
| fixed | cond | factor(tahun)2022 | 0.6027 | 0.0411 | 14.6647 | 0.0000 | 0.5222 | 0.6833 | 1.8271 | 1.6857 | 1.9803 |
| fixed | cond | factor(tahun)2023 | 0.8647 | 0.0471 | 18.3688 | 0.0000 | 0.7724 | 0.9569 | 2.3742 | 2.1650 | 2.6037 |
| fixed | cond | factor(tahun)2024 | 0.9521 | 0.0515 | 18.4774 | 0.0000 | 0.8511 | 1.0531 | 2.5911 | 2.3422 | 2.8664 |
| fixed | cond | factor(tahun)2025 | 0.7541 | 0.1181 | 6.3835 | 0.0000 | 0.5226 | 0.9857 | 2.1258 | 1.6864 | 2.6796 |
Tabel 20. Variansi Random Effect Wilayah dan Tahun
random_effect <- broom.mixed::tidy(
model_nb,
effects = "ran_pars"
)
knitr::kable(
random_effect,
digits = 4)
| effect | component | group | term | estimate |
|---|---|---|---|---|
| ran_pars | cond | nama_kabupaten_kota | sd__(Intercept) | 0.3059 |
Tabel 20 menunjukkan besaran variansi yang berasal dari perbedaan antarwilayah (kabupaten/kota) dan antarperiode (tahun). Variansi yang signifikan pada random effect ini menegaskan bahwa terdapat karakteristik spesifik wilayah yang tidak tertangkap oleh variabel prediktor tetap (fixed effects).
Pemilihan model yang tepat merupakan langkah krusial dalam pemodelan statistik untuk memastikan estimasi parameter yang akurat dan reliabel. Pada penelitian ini, dilakukan perbandingan kinerja antara tiga model, yaitu Ordinary Least Squares (OLS), Regresi Poisson, dan Negative Binomial Generalized Linear Mixed Model (GLMM). Perbandingan dilakukan dengan mengacu pada nilai Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), serta Log-Likelihood.
Tabel 21. Perbandingan Kinerja Model
perbandingan_model <- data.frame(
Model = c(
"OLS",
"Poisson",
"Negative Binomial GLMM"
),
AIC = c(
AIC(ols_model),
AIC(pois_model),
AIC(model_nb)
),
BIC = c(
BIC(ols_model),
BIC(pois_model),
BIC(model_nb)
),
LogLik = c(
as.numeric(logLik(ols_model)),
as.numeric(logLik(pois_model)),
as.numeric(logLik(model_nb))
)
)
knitr::kable(
perbandingan_model,
digits = 2
)
| Model | AIC | BIC | LogLik |
|---|---|---|---|
| OLS | 2577.78 | 2598.12 | -1281.89 |
| Poisson | 9373.57 | 9405.52 | -4675.78 |
| Negative Binomial GLMM | 2250.09 | 2284.96 | -1113.05 |
Berdasarkan Tabel 21, ringkasan kinerja model dapat dijelaskan sebagai berikut:
Dari kriteria informasi statistik tersebut, model Negative Binomial GLMM terpilih sebagai model terbaik, didukung oleh dua lapis bukti yang konsisten: bukti diagnostik (overdispersi) dan bukti komparatif (AIC/BIC/Log-Likelihood) pada struktur model yang setara. Pemilihan ini tidak hanya memberikan estimasi parameter yang lebih valid secara statistik dengan meminimalisir bias akibat overdispersi, tetapi juga memberikan perspektif yang lebih komprehensif dalam mengintegrasikan dinamika spasio-temporal tuberkulosis di Provinsi Jawa Barat.
prioritas <- data %>%
group_by(
nama_kabupaten_kota
) %>%
summarise(
rata_IR = mean(IR),
rata_SIR = mean(SIR)
) %>%
arrange(desc(rata_SIR))
kable(
head(prioritas,10)
)
| nama_kabupaten_kota | rata_IR | rata_SIR |
|---|---|---|
| KOTA CIREBON | 972.6074 | 2.641255 |
| KOTA BOGOR | 773.0801 | 2.099410 |
| KOTA SUKABUMI | 770.9885 | 2.093730 |
| KOTA CIMAHI | 663.1632 | 1.800915 |
| KOTA BANDUNG | 631.0095 | 1.713597 |
| KOTA BANJAR | 484.2271 | 1.314988 |
| KOTA TASIKMALAYA | 482.8470 | 1.311241 |
| KOTA BEKASI | 435.7330 | 1.183296 |
| KABUPATEN PURWAKARTA | 430.2649 | 1.168446 |
| KABUPATEN BOGOR | 420.2319 | 1.141200 |
Wilayah yang menempati 10 peringkat teratas dengan SIR tertinggi diidentifikasi sebagai klaster prioritas. Tingginya nilai SIR pada wilayah-wilayah ini menunjukkan bahwa upaya penemuan kasus dan intervensi kesehatan harus diprioritaskan di daerah tersebut karena risiko TB secara statistik jauh di atas rata-rata provinsi.
Analisis spasial-temporal ini menunjukkan bahwa transmisi TB di Jawa Barat tidak hanya dipengaruhi oleh faktor sosial ekonomi, tetapi juga oleh struktur populasi.
Pengaruh Kepadatan: Kepadatan penduduk yang tinggi meningkatkan probabilitas transmisi droplet (kontak antarindividu), yang berkontribusi pada peningkatan kasus TB secara signifikan.
Sinergi TB-HIV: Hubungan positif antara prevalensi HIV dan TB mengonfirmasi perlunya layanan integrasi TB-HIV, mengingat kerentanan individu dengan HIV terhadap infeksi TB akibat penurunan sistem imun.
Kebijakan Berbasis Wilayah: Hasil ini memberikan landasan bagi pemerintah daerah dalam merumuskan kebijakan yang bersifat presisi (precision public health), di mana intervensi tidak lagi bersifat seragam (one-size-fits-all), melainkan disesuaikan dengan profil risiko tiap kabupaten/kota.
Berdasarkan hasil analisis epidemiologi, disease mapping, dan pemodelan spasio-temporal Tuberkulosis di Provinsi Jawa Barat periode 2021–2025, dapat ditarik kesimpulan sebagai berikut:
Distribusi Kasus dan Ukuran Epidemiologi
Jumlah kasus TB di Jawa Barat menunjukkan tren peningkatan yang
konsisten selama periode pengamatan, dari 89.255 kasus pada tahun 2021
hingga mencapai puncaknya 229.683 kasus pada tahun 2024, kemudian
sedikit menurun menjadi 225.389 kasus pada tahun 2025. Incidence Rate
(IR) meningkat dari 183,13 per 100.000 penduduk (2021) menjadi 456,22
per 100.000 penduduk (2024), mencerminkan intensifikasi penemuan kasus
maupun peningkatan penularan aktual. Nilai Koefisien Variasi (CV) yang
berkisar 72–80% menunjukkan disparitas beban TB yang sangat tinggi
antarwilayah.
Pemetaan Risiko Empirical Bayes
Estimasi risiko berbasis Empirical Bayes menghasilkan peta risiko yang
lebih stabil dibandingkan SIR klasik, terutama pada wilayah dengan
populasi kecil. Hasil pemetaan mengidentifikasi sejumlah kabupaten/kota
dengan risiko EB yang secara konsisten tinggi lintas tahun,
menjadikannya prioritas utama dalam perencanaan intervensi pengendalian
TB di tingkat daerah.
Autokorelasi Spasial
Hasil uji Global Moran’s I menunjukkan adanya autokorelasi spasial yang
signifikan pada distribusi risiko TB di Jawa Barat, mengindikasikan
bahwa kasus TB tidak tersebar secara acak melainkan membentuk pola
pengelompokan (cluster) geografis. Analisis LISA lebih lanjut
mengidentifikasi wilayah klaster High-High (risiko tinggi dikelilingi
tetangga risiko tinggi) maupun wilayah outlier spasial, yang dapat
menjadi dasar penetapan prioritas intervensi berbasis
kewilayahan.
Pemodelan Spasio-Temporal GLMM Negative
Binomial
Data kasus TB menunjukkan gejala overdispersi yang signifikan (\(\phi_{disp} > 1\)), sehingga model
Poisson tidak memadai dan penggunaan GLMM Negative Binomial (nbinom2)
terbukti lebih sesuai. Model GLMM-NB yang menyertakan efek acak wilayah
(\(u_i\)) dan efek acak temporal (\(v_t\)) mampu,mengakomodasi heterogenitas
spasial dan fluktuasi temporal secara simultan. Hasil estimasi Incidence
Rate Ratio (IRR) menunjukkan bahwa variabel HIV, kepadatan penduduk, dan
persentase penduduk miskin berpengaruh positif terhadap laju kejadian
TB, sementara akses air minum layak dan sanitasi layak berperan sebagai
faktor protektif yang menurunkan risiko TB.
Berdasarkan hasil penelitian, beberapa saran disampaikan sebagai berikut:
Intervensi berbasis wilayah prioritas, Kabupaten/kota dengan nilai EB Risk dan SIR konsisten di atas 1 perlu mendapat alokasi sumber daya lebih intensif, mencakup skrining, contact tracing, dan penguatan fasilitas pengobatan.
Penanganan determinan sosial-ekonomi, Program pengendalian TB perlu diintegrasikan dengan peningkatan sanitasi, akses air minum layak, dan pengentasan kemiskinan, terutama di wilayah padat penduduk.
Sinergi TB-HIV, Penguatan program kolaborasi TB-HIV diperlukan di wilayah dengan beban HIV tinggi, mengingat HIV terbukti berpengaruh signifikan terhadap laju kejadian TB.
Peningkatan kualitas surveilans, Diperlukan sistem surveilans TB yang berkelanjutan dan terintegrasi secara digital di seluruh fasilitas kesehatan Jawa Barat.
Pengembangan model spasial, Penggunaan model
CAR/BYM dalam kerangka Bayesian Hierarchical (misalnya via
R-INLA atauCARBayes) serta penerapan
dekomposisi spasio-temporal dapat menghasilkan estimasi risiko yang
lebih akurat dengan mengintegrasikan dependensi spasial dan temporal
secara eksplisit.
Perluasan variabel dan unit analisis, Penelitian lanjutan disarankan menambahkan variabel seperti IPM, cakupan BCG, dan Diabetes Melitus, serta mempertimbangkan analisis pada unit kecamatan untuk ketajaman identifikasi hotspot yang lebih tinggi.