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 menganalisis variasi spasial dan temporal secara simultan, digunakan pendekatan spatio-temporal melalui Generalized Linear Mixed Model (GLMM). Pada penelitian ini, kabupaten/kota dimodelkan sebagai random effect untuk menangkap heterogenitas antarwilayah, sedangkan tahun dimasukkan sebagai fixed effect untuk mengendalikan variasi antarperiode pengamatan. Pendekatan ini sesuai untuk data panel epidemiologi yang memiliki pengamatan berulang pada setiap wilayah (Zuur et al., 2009).
Namun, data kasus TB sebagai data hitungan (count data) 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 equidispersion (kesamaan antara nilai rata-rata dan varians). Oleh karena itu, digunakan pendekatan Negative Binomial Generalized Linear Mixed Model (GLMM-NB) sebagai model yang lebih fleksibel dalam menangani overdispersion. Pada penelitian ini, model memasukkan random intercept pada tingkat kabupaten/kota untuk menangkap heterogenitas antarwilayah, sedangkan tahun dimasukkan sebagai fixed effect untuk mengendalikan perubahan jumlah kasus pada setiap periode pengamatan (Hilbe, 2011).
Dengan demikian, pendekatan Negative Binomial Generalized Linear Mixed Model (GLMM-NB) dalam penelitian ini diharapkan mampu memberikan hasil pemodelan yang lebih akurat dalam menjelaskan faktor-faktor yang memengaruhi jumlah kasus Tuberkulosis, menjelaskan variasi kasus antar kabupaten/kota melalui random effect, serta mengendalikan perubahan jumlah kasus pada setiap periode pengamatan melalui fixed effect tahun di Provinsi Jawa Barat selama periode 2021-2025.
Berdasarkan latar belakang yang telah diuraikan, maka rumusan masalah dalam penelitian ini adalah sebagai berikut:
Bagaimana gambaran distribusi kasus Tuberkulosis (TB) di Provinsi Jawa Barat selama periode 2021-2025 berdasarkan ukuran epidemiologi (Incidence Rate dan Standardized Incidence Ratio)?
Bagaimana pemetaan risiko Tuberkulosis di Provinsi Jawa Barat menggunakan pendekatan disease mapping berbasis Empirical Bayes?
Apakah terdapat pola autokorelasi spasial (cluster) pada risiko Tuberkulosis di Provinsi Jawa Barat berdasarkan analisis Global Moran’s I dan Local Indicators of Spatial Association (LISA)?
Bagaimana pengaruh faktor sosial, ekonomi, dan lingkungan terhadap jumlah kasus Tuberkulosis di Provinsi Jawa Barat menggunakan pendekatan Generalized Linear Mixed Model (GLMM) Negative Binomial?
Berdasarkan rumusan masalah yang telah disusun, maka tujuan penelitian ini adalah sebagai berikut:
Menganalisis distribusi kasus Tuberkulosis (TB) di Provinsi Jawa Barat selama periode 2021-2025 berdasarkan ukuran epidemiologi (Incidence Rate dan Standardized Incidence Ratio).
Melakukan pemetaan risiko Tuberkulosis di Provinsi Jawa Barat menggunakan pendekatan disease mapping berbasis Empirical Bayes.
Mengidentifikasi pola autokorelasi spasial (cluster) pada risiko Tuberkulosis di Provinsi Jawa Barat menggunakan analisis Global Moran’s I dan Local Indicators of Spatial Association (LISA).
Menganalisis pengaruh faktor sosial, ekonomi, dan lingkungan terhadap jumlah kasus Tuberkulosis di Provinsi Jawa Barat menggunakan pendekatan Generalized Linear Mixed Model (GLMM) Negative Binomial.
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\) merupakan 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.
Generalized Linear Mixed Model (GLMM) merupakan pengembangan dari Generalized Linear Model (GLM) yang memungkinkan dimasukkannya efek acak (random effect) untuk mengakomodasi korelasi antar pengamatan yang berasal dari kelompok yang sama. Pada penelitian ini, GLMM digunakan untuk memodelkan jumlah kasus Tuberkulosis yang bersifat data cacah (count data) dan mengalami overdispersi. Oleh karena itu, digunakan distribusi Negative Binomial tipe nbinom2 melalui paket glmmTMB.
Model GLMM dalam penelitian ini memasukkan kabupaten/kota sebagai random intercept untuk mengakomodasi heterogenitas antarwilayah yang tidak dapat dijelaskan oleh variabel prediktor, sedangkan variabel tahun dimasukkan sebagai fixed effect untuk mengendalikan perubahan jumlah kasus Tuberkulosis pada setiap periode pengamatan. Selain itu, digunakan offset berupa logaritma jumlah penduduk sehingga model mengestimasi laju kejadian (incidence rate) setelah memperhitungkan ukuran populasi (Hardin & Hilbe, 2007).
\[ Var(TB_{it}\mid u_i) = \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:
\[ \begin{aligned} \log(\mu_{it}) =\;& \beta_0 +\beta_1X_{1it} +\beta_2X_{2it} +\beta_3X_{3it} +\beta_4X_{4it} +\beta_5X_{5it} \\ & +\beta_6Tahun_{it} +\log(Penduduk_{it}) +u_i \end{aligned} \]
Dengan:
\(TB_{it}\) : jumlah kasus Tuberkulosis pada kabupaten/kota ke-\(i\) dan tahun ke-\(t\)
\(\mu_{it}\) : nilai harapan jumlah kasus Tuberkulosis
\(\beta_0\) : intersep model
\(\beta_1,\ldots,\beta_6\) : koefisien regresi (fixed effect)
\(X_{1it}\) : persentase penduduk miskin
\(X_{2it}\) : kepadatan penduduk
\(X_{3it}\) : persentase rumah tangga dengan akses air minum layak
\(X_{4it}\) : persentase rumah tangga dengan akses sanitasi layak
\(X_{5it}\) : jumlah kasus HIV
\(Tahun_{it}\) : variabel indikator tahun pengamatan (fixed effect)
\(\log(Penduduk_{it})\) : offset berupa logaritma jumlah penduduk
\(u_i\) : random intercept kabupaten/kota dengan \[ u_i\sim N(0,\sigma_u^2) \] Pada penelitian ini, kabupaten/kota diperlakukan sebagai random effect (random intercept) untuk menangkap heterogenitas antarwilayah yang tidak sepenuhnya dapat dijelaskan oleh variabel prediktor dalam model. Sementara itu, tahun dimasukkan sebagai fixed effect karena periode pengamatan (2021-2025) merupakan rentang waktu yang telah ditetapkan dalam penelitian, sehingga pengaruh setiap tahun diestimasi secara langsung melalui parameter tetap. Pendekatan ini memungkinkan model menangkap variasi antarwilayah sekaligus mengendalikan perubahan jumlah kasus Tuberkulosis pada setiap periode pengamatan.
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
)
| 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
)
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 penggunaan random intercept pada tingkat kabupaten/kota untuk menangkap heterogenitas antarwilayah, sedangkan variasi antarperiode dikendalikan melalui variabel tahun sebagai fixed effect.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 Kabupaten/Kota. Variansi yang teramati pada efek acak ini menunjukkan bahwa terdapat keragaman risiko dasar ( baseline risk) yang menunjukkan adanya keragaman risiko dasar antar kabupaten/kota yang belum dapat dijelaskan oleh variabel prediktor tetap 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. Estimasi Simpangan Baku Random Effect Wilayah
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 estimasi simpangan baku random intercept pada tingkat kabupaten/kota sebesar 0,3059. Nilai tersebut menunjukkan adanya variasi karakteristik dasar antar kabupaten/kota yang belum dapat dijelaskan oleh variabel tetap (fixed effects) dalam model. Hal ini menunjukkan bahwa setiap kabupaten/kota memiliki perbedaan karakteristik yang memengaruhi jumlah kasus TB, sehingga penggunaan random effect pada tingkat wilayah diperlukan untuk mengakomodasi heterogenitas tersebut.
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:
Interpretasi Nilai AIC dan BIC, Nilai AIC dan BIC berfungsi sebagai indikator goodness-of-fit yang menghukum model atas kompleksitas yang berlebihan (parsimony). Tabel tersebut menunjukkan bahwa model Negative Binomial GLMM menghasilkan nilai AIC dan BIC yang paling rendah dibandingkan dengan model OLS maupun model Poisson. Perlu ditekankan bahwa perbandingan yang paling relevan secara metodologis adalah antara model Poisson dan Negative Binomial GLMM, mengingat kedua model tersebut memiliki struktur yang identik (fixed effects yang sama, offset populasi, serta random intercept wilayah) dan hanya berbeda pada asumsi distribusi terhadap varians. Penurunan AIC dan BIC yang nyata pada model NB GLMM dibandingkan Poisson GLMM dengan demikian secara langsung mengonfirmasi bahwa penambahan parameter dispersi mampu memperbaiki kecocokan model secara substansial, bukan sekadar artefak dari penambahan kompleksitas. Sementara itu, nilai AIC/BIC model OLS yang jauh lebih tinggi terutama berfungsi sebagai baseline pembanding yang menegaskan bahwa pendekatan linear sederhana sejak awal tidak sesuai untuk memodelkan data cacah (count data) seperti kasus TB, sehingga perbandingannya dengan dua model lain bersifat ilustratif belaka, bukan perbandingan yang setara (apple-to-apple).
Analisis Log-Likelihood, Nilai Log-Likelihood yang lebih tinggi pada model Negative Binomial GLMM dibandingkan model Poisson GLMM memperkuat temuan di atas, karena kenaikan likelihood ini terjadi pada struktur model yang sebanding, sehingga dapat diatribusikan secara spesifik pada perbaikan asumsi distribusi (Negative Binomial) dalam menangani overdispersi, bukan pada perbedaan spesifikasi model yang mendasar.
Justifikasi Pemilihan Model, Ketidaksesuaian model OLS terhadap data cacah serta kegagalan model Poisson dalam memenuhi asumsi equidispersion - sebagaimana telah dibuktikan secara kuantitatif melalui uji dispersi Pearson pada Tabel 17 - menjadi dasar utama mengapa kedua model tersebut menghasilkan nilai AIC/BIC yang lebih tinggi. Bukti uji overdispersi dan bukti penurunan AIC/BIC pada model NB GLMM ini bersifat saling melengkapi: uji dispersi menjelaskan mengapa asumsi Poisson dilanggar, sementara kriteria informasi menunjukkan secara empiris seberapa besar perbaikan yang diperoleh setelah beralih ke distribusi Negative Binomial. Dengan memasukkan parameter dispersi tambahan, random intercept pada tingkat kabupaten/kota, serta tahun sebagai fixed effect, model Negative Binomial GLMM memberikan kecocokan model yang lebih baik dalam menjelaskan variasi jumlah kasus Tuberkulosis antar kabupaten/kota sekaligus mengendalikan perubahan jumlah kasus pada setiap periode pengamatan di Jawa Barat.
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.
Meskipun penelitian ini telah berhasil menjawab tujuan yang dirumuskan, terdapat beberapa keterbatasan yaitu sebagai berikut:
Keterbatasan Unit Analisis (Ecological Fallacy), Data yang digunakan bersifat agregat tingkat kabupaten/kota, bukan data individu, sehingga berisiko menimbulkan ecological fallacy-hubungan yang teramati pada tingkat wilayah belum tentu berlaku pada tingkat individu. Misalnya, hubungan positif kepadatan penduduk dengan kasus TB tidak berarti setiap individu di wilayah padat otomatis berisiko lebih tinggi, karena variasi risiko intra-wilayah tidak tertangkap data agregat.
Tidak Tersedianya Data Karakteristik Individu, Variabel level individu seperti usia, status gizi, riwayat kontak, dan riwayat pengobatan tidak tersedia, sehingga model belum dapat mengidentifikasi kelompok berisiko tinggi (high-risk subgroups) secara spesifik untuk intervensi yang lebih tepat sasaran.
Variabel Perilaku dan Kualitas Pelayanan Kesehatan Belum Terakomodasi, Faktor perilaku kesehatan (health-seeking behavior, kepatuhan pengobatan) serta indikator kualitas layanan seperti case detection rate dan treatment success rate belum dimasukkan ke dalam model akibat keterbatasan data sekunder, padahal keduanya berpotensi memengaruhi variasi kasus antarwilayah.
Potensi Bias Pelaporan (Reporting Bias), Jumlah kasus tercatat sangat bergantung pada kapasitas surveilans masing-masing wilayah. Wilayah dengan akses diagnostik lebih baik (misalnya TCM) berpotensi melaporkan kasus lebih tinggi bukan karena insidensi riil yang lebih besar, melainkan kemampuan deteksi yang lebih unggul-sehingga estimasi risiko spasial dapat bias.
Periode Waktu dan Cakupan Spasial Terbatas, Penelitian hanya mencakup lima tahun (2021-2025) dan 27 kabupaten/kota di Jawa Barat, sehingga generalisasi ke wilayah lain atau proyeksi jangka panjang perlu kehati-hatian. Periode ini juga beririsan dengan masa pemulihan pasca-pandemi COVID-19 yang dapat memengaruhi pola pelaporan secara tidak proporsional.
Keterbatasan Pendekatan Spasial yang Digunakan, Pendekatan Empirical Bayes Global memperlakukan setiap wilayah secara independen tanpa mempertimbangkan struktur ketetanggaan spasial, berbeda dengan model CAR atau BYM yang mengakomodasi dependensi spasial secara langsung. Estimasi risiko EB pada penelitian ini berpotensi belum optimal menangkap pengaruh spasial lokal.
Berdasarkan hasil analisis epidemiologi, disease mapping, dan pemodelan spasio-temporal Tuberkulosis di Provinsi Jawa Barat periode 2021-2025, dapat ditarik kesimpulan sebagai berikut:
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 peningkatan jumlah kasus yang dapat dipengaruhi oleh intensifikasi penemuan kasus maupun peningkatan penularan di masyarakat. Nilai Koefisien Variasi (CV) yang berkisar 72-80% menunjukkan disparitas beban TB yang sangat tinggi antarwilayah.
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.
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.
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 memasukkan random intercept pada tingkat kabupaten/kota serta tahun sebagai fixed effect mampu menjelaskan variasi jumlah kasus TB antarwilayah sekaligus mengendalikan perubahan jumlah kasus pada setiap periode pengamatan. Hasil estimasi Incidence Rate Ratio (IRR) menunjukkan bahwa variabel HIV, kepadatan penduduk, dan persentase penduduk miskin berpengaruh positif terhadap laju kejadian TB, sedangkan akses air minum layak dan sanitasi layak berperan sebagai faktor protektif yang menurunkan laju kejadian 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 dapat dipertimbangkan pada penelitian selanjutnya untuk memodelkan ketergantungan spasial secara lebih eksplisit serta menghasilkan estimasi risiko yang lebih stabil. Selain itu, pengembangan model yang mengintegrasikan komponen spasio-temporal juga dapat dilakukan apabila tersedia data dengan periode pengamatan yang lebih panjang.
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.