Analisis Spasio-Temporal Risiko Diare di Provinsi Jawa Barat Tahun 2020–2024
paket <- c(
"readxl", "dplyr", "janitor", "ggplot2", "writexl", "scales",
"stringr", "sf", "spdep", "broom", "viridis", "tidyr", "forcats",
"knitr", "MASS"
)
paket_belum_ada <- paket[!paket %in% installed.packages()[, "Package"]]
if(length(paket_belum_ada) > 0){
install.packages(paket_belum_ada)
}
library(readxl)
library(dplyr)
library(janitor)
library(ggplot2)
library(writexl)
library(scales)
library(stringr)
library(sf)
library(spdep)
library(broom)
library(viridis)
library(tidyr)
library(forcats)
library(knitr)
# MASS tidak dipanggil dengan library(MASS) agar tidak menimpa dplyr::select().
# Fungsi Negative Binomial dipanggil sebagai MASS::glm.nb().standar_wilayah <- function(x){
x %>%
as.character() %>%
stringr::str_to_upper() %>%
stringr::str_replace_all("\\.", "") %>%
stringr::str_replace_all("KAB\\s+", "KABUPATEN ") %>%
stringr::str_replace_all("KABUPATEN\\s+KABUPATEN", "KABUPATEN") %>%
stringr::str_replace_all("KOTA\\s+KOTA", "KOTA") %>%
stringr::str_squish()
}
fmt_angka <- function(x, digits = 0){
formatC(x, format = "f", digits = digits, big.mark = ".", decimal.mark = ",")
}
tema_publikasi <- function(){
ggplot2::theme_minimal(base_size = 12) +
ggplot2::theme(
plot.title = ggplot2::element_text(face = "bold", color = "#1e40af", size = 15),
plot.subtitle = ggplot2::element_text(color = "#475569"),
axis.title = ggplot2::element_text(color = "#1f2937"),
panel.grid.minor = ggplot2::element_blank(),
legend.title = ggplot2::element_text(face = "bold"),
legend.position = "right"
)
}Nama: Mahdayani Putri Yunizar
NPM: 140720250004
Mata Kuliah: Epidemiologi
Program Studi Magister Statistika Terapan
Fakultas Matematika dan Ilmu Pengetahuan Alam
Universitas Padjadjaran
1 Abstrak
Diare merupakan masalah kesehatan masyarakat yang memiliki variasi risiko antarwilayah dan antarwaktu. Penelitian ini bertujuan menganalisis risiko diare di Provinsi Jawa Barat tahun 2020–2024 melalui pendekatan ukuran epidemiologi, disease mapping, dan model spasio-temporal. Data yang digunakan merupakan data agregat 27 kabupaten/kota selama lima tahun, meliputi jumlah kasus diare, populasi berisiko, serta kovariat wilayah seperti kepadatan penduduk, kemiskinan, curah hujan, kelembapan, suhu, sanitasi, fasilitas kesehatan, dan tenaga kesehatan. Metode analisis meliputi perhitungan incidence rate, expected cases, standardized incidence ratio (SIR), pemetaan crude incidence rate dan crude SIR, smoothing Empirical Bayes, serta model spasio-temporal Negative Binomial dengan efek wilayah dan waktu. Hasil menunjukkan bahwa kasus diare relatif tinggi pada tahun 2020–2023 dan menurun pada tahun 2024. Wilayah prioritas tahun 2024 adalah Kota Cirebon, Kota Bogor, Kota Sukabumi, Kota Tasikmalaya, dan Kabupaten Majalengka. Kesimpulannya, risiko diare di Jawa Barat tidak tersebar merata, sehingga monitoring dan intervensi kesehatan masyarakat perlu difokuskan pada wilayah prioritas.
2 Pendahuluan
Diare merupakan salah satu masalah kesehatan masyarakat yang masih perlu mendapatkan perhatian karena dapat menimbulkan beban kesakitan yang besar dan distribusinya dapat berbeda antarwilayah serta antarwaktu. Studi berbasis Global Burden of Disease 2021 menunjukkan bahwa beban diare di Asia mengalami perubahan dalam jangka panjang, tetapi penyakit ini tetap menjadi masalah penting karena masih terdapat variasi beban menurut negara, kelompok umur, dan tingkat sosial ekonomi. Bahkan, proyeksi untuk beberapa negara Asia, termasuk Indonesia, menunjukkan kecenderungan peningkatan angka insidensi terstandar pada periode mendatang (Liu et al., 2025). Hal ini menunjukkan bahwa diare masih relevan untuk dianalisis, terutama dengan pendekatan yang mampu menggambarkan perbedaan risiko menurut wilayah dan waktu.
Risiko diare tidak hanya dipengaruhi oleh faktor biologis, tetapi juga berkaitan dengan kondisi lingkungan dan perilaku kesehatan. Adria et al. (2023) menunjukkan bahwa cakupan air minum layak dan cakupan jamban sehat berhubungan signifikan dengan kejadian diare. Penelitian tersebut juga menunjukkan bahwa pemetaan spasial dapat membantu menggambarkan wilayah dengan kejadian diare yang lebih tinggi, terutama pada daerah dengan cakupan air minum layak yang lebih rendah. Temuan ini menunjukkan bahwa diare perlu dianalisis sebagai masalah kesehatan yang berkaitan erat dengan kondisi lingkungan wilayah.
Selain sanitasi dan akses air minum, faktor iklim juga berperan dalam perubahan risiko diare. Geremew et al. (2024) melalui systematic review dan meta-analysis menemukan bahwa peningkatan suhu dan curah hujan berhubungan dengan peningkatan kejadian diare pada anak. Penelitian tersebut juga menegaskan bahwa kondisi water, sanitation, and hygiene (WASH) dapat memodifikasi hubungan antara faktor iklim dan kejadian diare. Dengan demikian, variabel lingkungan seperti curah hujan, kelembapan, dan suhu penting untuk dipertimbangkan dalam analisis diare, terutama pada wilayah tropis seperti Indonesia.
Dalam konteks Provinsi Jawa Barat, pendekatan spasial menjadi penting karena setiap kabupaten/kota memiliki karakteristik kepadatan penduduk, sanitasi, akses fasilitas kesehatan, dan kondisi lingkungan yang berbeda. Wahyudi et al. (2024) melakukan analisis autokorelasi spasial kasus diare pada 27 kabupaten/kota di Jawa Barat tahun 2023 dan menunjukkan bahwa analisis spasial dapat membantu mengidentifikasi pola sebaran serta wilayah berisiko tinggi. Penelitian tersebut memperkuat pentingnya pemetaan risiko diare di Jawa Barat agar wilayah prioritas dapat ditentukan secara lebih objektif.
Analisis diare tidak cukup hanya menggunakan jumlah kasus absolut karena wilayah dengan jumlah penduduk besar cenderung memiliki kasus yang lebih banyak. Oleh karena itu, diperlukan ukuran epidemiologi yang mempertimbangkan populasi berisiko, seperti incidence rate, expected cases, dan standardized incidence ratio (SIR). Incidence rate digunakan untuk menggambarkan angka kejadian diare per 100.000 penduduk, sedangkan SIR digunakan untuk membandingkan jumlah kasus teramati dengan jumlah kasus yang diharapkan berdasarkan rata-rata provinsi. Nilai SIR lebih dari 1 menunjukkan bahwa suatu wilayah memiliki risiko lebih tinggi dibandingkan rata-rata provinsi pada tahun yang sama.
Disease mapping digunakan untuk memvisualisasikan variasi risiko penyakit antarwilayah secara lebih komunikatif. Coly et al. (2021) menjelaskan bahwa disease mapping bertujuan memperkirakan risiko penyakit yang mendasari data epidemiologi dan menyajikannya dalam bentuk peta risiko yang diperhalus. Penelitian tersebut juga menjelaskan bahwa model spasio-temporal dapat digunakan untuk mengakomodasi dimensi waktu, kovariat, dan overdispersion. Dalam konteks data penyakit menular, Coly et al. (2021) menunjukkan bahwa model Negative Binomial dapat menjadi alternatif yang lebih sesuai dibandingkan Poisson ketika data mengalami overdispersion.
Berdasarkan latar belakang tersebut, analisis ini dilakukan untuk mengkaji risiko diare di Provinsi Jawa Barat tahun 2020–2024 dengan pendekatan epidemiologi spasial dan spasio-temporal. Tujuan analisis ini adalah: (1) mendeskripsikan kasus diare dan populasi berisiko menurut tahun; (2) menghitung incidence rate, expected cases, dan SIR pada setiap kabupaten/kota; (3) memetakan variasi risiko diare menggunakan peta crude incidence rate, crude SIR, dan smoothing Empirical Bayes; (4) menganalisis perubahan risiko diare menurut wilayah dan waktu menggunakan model Negative Binomial; serta (5) mengidentifikasi wilayah prioritas sebagai dasar rekomendasi monitoring dan intervensi kesehatan masyarakat.
3 Data dan Metode
3.1 Sumber Data
link_data_geojson <- "https://drive.google.com/drive/folders/1n-LfSxZWyIeYisQTiex7atgh7qUCtFuV?usp=sharing"
kandidat_excel <- c(
"diare.xlsx",
"diare(1).xlsx",
"diare(2).xlsx"
)
kandidat_geojson <- c(
"geojson_jawa_barat_adm2.geojson",
"jawa_barat_adm2.geojson",
"jabar_adm2.geojson"
)
file_excel <- kandidat_excel[file.exists(kandidat_excel)][1]
file_geojson <- kandidat_geojson[file.exists(kandidat_geojson)][1]
if(is.na(file_excel)){
stop(
"File Excel tidak ditemukan. Simpan salah satu file berikut di folder yang sama dengan Rmd: ",
paste(kandidat_excel, collapse = ", ")
)
}
if(is.na(file_geojson)){
stop(
"File GeoJSON tidak ditemukan. Simpan geojson_jawa_barat_adm2.geojson di folder yang sama dengan Rmd."
)
}
data_raw <- readxl::read_excel(file_excel)
peta_raw <- sf::st_read(file_geojson, quiet = TRUE)Penelitian ini menggunakan data agregat kasus diare tingkat kabupaten/kota di Provinsi Jawa Barat tahun 2020–2024. Unit analisis dalam penelitian ini adalah kabupaten/kota per tahun. Dengan demikian, data berbentuk panel wilayah-waktu yang terdiri atas 27 kabupaten/kota selama 5 tahun pengamatan.
Data utama mencakup jumlah kasus diare, jumlah penduduk, serta beberapa kovariat wilayah yang diduga berkaitan dengan risiko diare. Selain data tabular, penelitian ini juga menggunakan data spasial berupa file GeoJSON batas administrasi kabupaten/kota di Provinsi Jawa Barat. Data dan GeoJSON yang digunakan dalam analisis ini tersedia pada tautan berikut: https://drive.google.com/drive/folders/1n-LfSxZWyIeYisQTiex7atgh7qUCtFuV?usp=sharing.
Data yang digunakan bersifat agregat wilayah dan tidak memuat identitas individu. Oleh karena itu, hasil analisis digunakan untuk menggambarkan variasi risiko diare antarwilayah dan antarwaktu, bukan untuk menyimpulkan risiko pada tingkat individu.
data <- data_raw %>%
janitor::clean_names() %>%
dplyr::rename(
wilayah = wilayah,
tahun = tahun,
kasus = diare_kasus,
kepadatan = kepadatan_jiwa_km2,
kemiskinan = kemiskinan_persen,
curah_hujan = curah_hujan_mm,
kelembapan = kelembapan_persen,
suhu = suhu_c,
sanitasi = sanitasi,
fasilitas_kesehatan = jumlah_fasilitas_kesehatan_unit,
tenaga_kesehatan = jumlah_tenaga_kesehatan_orang,
penduduk_ribu = jumlah_penduduk_ribu,
total_kasus_jabar_file = total_kasus_jabar_per_tahun,
total_penduduk_jabar_ribu_file = total_penduduk_jabar_ribu
) %>%
dplyr::mutate(
wilayah = standar_wilayah(wilayah),
tahun = as.integer(tahun),
kasus = as.numeric(kasus),
kepadatan = as.numeric(kepadatan),
kemiskinan = as.numeric(kemiskinan),
curah_hujan = as.numeric(curah_hujan),
kelembapan = as.numeric(kelembapan),
suhu = as.numeric(suhu),
sanitasi = as.numeric(sanitasi),
fasilitas_kesehatan = as.numeric(fasilitas_kesehatan),
tenaga_kesehatan = as.numeric(tenaga_kesehatan),
penduduk_ribu = as.numeric(penduduk_ribu),
populasi_berisiko = penduduk_ribu * 1000
)
kolom_wilayah_peta <- if("wilayah" %in% names(peta_raw)) "wilayah" else "shapeName"
peta <- peta_raw %>%
sf::st_make_valid() %>%
dplyr::mutate(
wilayah = standar_wilayah(.data[[kolom_wilayah_peta]])
)
jumlah_observasi <- nrow(data)
jumlah_wilayah <- dplyr::n_distinct(data$wilayah)
jumlah_tahun <- dplyr::n_distinct(data$tahun)
tahun_awal <- min(data$tahun, na.rm = TRUE)
tahun_akhir <- max(data$tahun, na.rm = TRUE)Struktur data yang dianalisis:
Data terdiri atas 135 observasi wilayah-tahun, mencakup 27
kabupaten/kota di Provinsi Jawa Barat selama 5 tahun pengamatan, yaitu
2020–2024.
3.2 Variabel Penelitian
Variabel respons dalam penelitian ini adalah jumlah kasus diare pada setiap kabupaten/kota dan tahun. Populasi berisiko dihitung dari jumlah penduduk pada masing-masing wilayah. Karena kolom jumlah penduduk dalam data menggunakan satuan ribu, nilai penduduk dikalikan 1.000 untuk memperoleh jumlah populasi berisiko dalam satuan orang.
tabel_variabel <- tibble::tribble(
~`Kelompok Variabel`, ~`Nama Variabel`, ~`Keterangan`,
"Variabel respons", "Kasus diare", "Jumlah kasus diare pada setiap kabupaten/kota dan tahun",
"Populasi berisiko", "Jumlah penduduk", "Jumlah penduduk wilayah, digunakan sebagai denominator",
"Waktu", "Tahun", "Tahun pengamatan, yaitu 2020--2024",
"Wilayah", "Kabupaten/kota", "Unit spasial analisis di Provinsi Jawa Barat",
"Kovariat sosial", "Kepadatan penduduk, kemiskinan", "Menggambarkan karakteristik sosial-demografis wilayah",
"Kovariat lingkungan", "Curah hujan, kelembapan, suhu", "Menggambarkan kondisi iklim dan lingkungan",
"Kovariat sanitasi", "Sanitasi", "Menggambarkan kondisi sanitasi wilayah",
"Kovariat pelayanan kesehatan", "Fasilitas kesehatan, tenaga kesehatan", "Menggambarkan ketersediaan layanan kesehatan"
)
knitr::kable(tabel_variabel, caption = "Variabel yang digunakan dalam analisis.")| Kelompok Variabel | Nama Variabel | Keterangan |
|---|---|---|
| Variabel respons | Kasus diare | Jumlah kasus diare pada setiap kabupaten/kota dan tahun |
| Populasi berisiko | Jumlah penduduk | Jumlah penduduk wilayah, digunakan sebagai denominator |
| Waktu | Tahun | Tahun pengamatan, yaitu 2020–2024 |
| Wilayah | Kabupaten/kota | Unit spasial analisis di Provinsi Jawa Barat |
| Kovariat sosial | Kepadatan penduduk, kemiskinan | Menggambarkan karakteristik sosial-demografis wilayah |
| Kovariat lingkungan | Curah hujan, kelembapan, suhu | Menggambarkan kondisi iklim dan lingkungan |
| Kovariat sanitasi | Sanitasi | Menggambarkan kondisi sanitasi wilayah |
| Kovariat pelayanan kesehatan | Fasilitas kesehatan, tenaga kesehatan | Menggambarkan ketersediaan layanan kesehatan |
3.3 Ukuran Epidemiologi
Ukuran epidemiologi yang digunakan meliputi incidence rate, expected cases, dan standardized incidence ratio (SIR). Ukuran ini digunakan karena jumlah kasus absolut tidak cukup untuk membandingkan risiko antarwilayah.
\[ Incidence\ Rate = \frac{Jumlah\ Kasus}{Populasi\ Berisiko} \times 100.000 \]
\[ Expected\ Cases = Populasi\ Wilayah \times \frac{Total\ Kasus\ Provinsi}{Total\ Populasi\ Provinsi} \]
\[ SIR = \frac{Kasus\ Teramati}{Expected\ Cases} \]
Interpretasi SIR dilakukan sebagai berikut: SIR > 1 menunjukkan risiko lebih tinggi daripada rata-rata provinsi pada tahun yang sama; SIR = 1 menunjukkan risiko setara dengan rata-rata provinsi; dan SIR < 1 menunjukkan risiko lebih rendah daripada rata-rata provinsi.
3.4 Model Spasial dan Spasio-Temporal
Analisis spasial dilakukan melalui disease mapping dengan menggabungkan data epidemiologi dan GeoJSON kabupaten/kota. Peta yang dibuat terdiri atas peta crude incidence rate, peta crude SIR, dan peta hasil smoothing Empirical Bayes. Estimasi risiko kasar dapat menjadi tidak stabil, terutama pada wilayah dengan populasi kecil atau jumlah kasus yang fluktuatif. Oleh karena itu, metode Empirical Bayes digunakan untuk menstabilkan estimasi risiko.
Model spasio-temporal digunakan untuk menganalisis perubahan risiko diare menurut wilayah dan waktu. Model Poisson digunakan sebagai pembanding awal untuk memeriksa overdispersion. Jika variasi data lebih besar daripada asumsi Poisson, model Negative Binomial lebih sesuai karena memiliki parameter tambahan untuk menangani overdispersion. Model utama yang digunakan dalam penelitian ini adalah Negative Binomial dengan efek wilayah, efek tahun, kovariat wilayah, serta offset log populasi berisiko.
\[ Y_{it} \sim Negative\ Binomial(\mu_{it}) \]
\[ \log(\mu_{it}) = \log(Populasi_{it}) + \alpha + Wilayah_i + Tahun_t + \beta X_{it} \]
Keterangan: \(Y_{it}\) adalah jumlah kasus diare pada wilayah ke-\(i\) dan tahun ke-\(t\), \(\mu_{it}\) adalah nilai harapan kasus, \(Populasi_{it}\) adalah populasi berisiko, \(Wilayah_i\) adalah efek wilayah, \(Tahun_t\) adalah efek waktu, dan \(X_{it}\) adalah kovariat wilayah.
4 Hasil dan Pembahasan
Bagian ini menyajikan tabel, grafik, peta, hasil model, interpretasi, dan implikasi epidemiologis. Setiap tabel dan gambar utama dilengkapi interpretasi agar hasil analisis lebih mudah dipahami.
# A. Ringkasan deskriptif
ringkasan_tahun <- data %>%
dplyr::group_by(tahun) %>%
dplyr::summarise(
total_kasus = sum(kasus, na.rm = TRUE),
total_populasi = sum(populasi_berisiko, na.rm = TRUE),
incidence_rate_provinsi = total_kasus / total_populasi * 100000,
rata_rata_kasus_wilayah = mean(kasus, na.rm = TRUE),
median_kasus_wilayah = median(kasus, na.rm = TRUE),
minimum_kasus_wilayah = min(kasus, na.rm = TRUE),
maksimum_kasus_wilayah = max(kasus, na.rm = TRUE),
.groups = "drop"
)
ringkasan_wilayah <- data %>%
dplyr::group_by(wilayah) %>%
dplyr::summarise(
total_kasus_semua_tahun = sum(kasus, na.rm = TRUE),
rata_rata_kasus_tahunan = mean(kasus, na.rm = TRUE),
total_populasi_semua_tahun = sum(populasi_berisiko, na.rm = TRUE),
incidence_rate_gabungan = total_kasus_semua_tahun / total_populasi_semua_tahun * 100000,
.groups = "drop"
) %>%
dplyr::arrange(dplyr::desc(incidence_rate_gabungan))
# B. Ukuran epidemiologi
total_tahun <- data %>%
dplyr::group_by(tahun) %>%
dplyr::summarise(
total_kasus_tahun = sum(kasus, na.rm = TRUE),
total_populasi_tahun = sum(populasi_berisiko, na.rm = TRUE),
rate_tahun = total_kasus_tahun / total_populasi_tahun,
.groups = "drop"
)
hasil_ukuran <- data %>%
dplyr::left_join(total_tahun, by = "tahun") %>%
dplyr::mutate(
incidence_rate = kasus / populasi_berisiko * 100000,
expected_cases = populasi_berisiko * rate_tahun,
sir = kasus / expected_cases,
interpretasi_sir = dplyr::case_when(
sir > 1 ~ "Risiko lebih tinggi daripada rata-rata provinsi",
sir == 1 ~ "Risiko sama dengan rata-rata provinsi",
sir < 1 ~ "Risiko lebih rendah daripada rata-rata provinsi",
TRUE ~ NA_character_
),
kategori_sir = dplyr::case_when(
sir >= 2 ~ "Sangat tinggi",
sir >= 1.5 ~ "Tinggi",
sir > 1 ~ "Di atas rata-rata",
sir == 1 ~ "Sama dengan rata-rata",
sir < 1 ~ "Di bawah rata-rata",
TRUE ~ NA_character_
)
)
tahun_peta <- max(hasil_ukuran$tahun, na.rm = TRUE)
prioritas_sir_tahun_terakhir <- hasil_ukuran %>%
dplyr::filter(tahun == tahun_peta) %>%
dplyr::arrange(dplyr::desc(sir)) %>%
dplyr::mutate(peringkat_sir = dplyr::row_number()) %>%
dplyr::select(
peringkat_sir, wilayah, tahun, kasus, populasi_berisiko,
incidence_rate, expected_cases, sir, kategori_sir, interpretasi_sir
)
# C. Peta dan Empirical Bayes
data_peta_tahun <- hasil_ukuran %>%
dplyr::filter(tahun == tahun_peta)
peta_tahun <- peta %>%
dplyr::left_join(data_peta_tahun, by = "wilayah")
cek_join_gagal <- peta_tahun %>%
sf::st_drop_geometry() %>%
dplyr::filter(is.na(kasus)) %>%
dplyr::select(wilayah)
peta_eb <- peta_tahun %>%
dplyr::filter(!is.na(kasus), !is.na(populasi_berisiko), populasi_berisiko > 0)
nb <- spdep::poly2nb(peta_eb, queen = TRUE)
hasil_eb <- spdep::EBlocal(
peta_eb$kasus,
peta_eb$populasi_berisiko,
nb,
zero.policy = TRUE
)
rate_provinsi_tahun_peta <- sum(peta_eb$kasus, na.rm = TRUE) /
sum(peta_eb$populasi_berisiko, na.rm = TRUE)
peta_eb$eb_rate <- hasil_eb$est
peta_eb$eb_incidence_rate <- peta_eb$eb_rate * 100000
peta_eb$eb_sir <- peta_eb$eb_rate / rate_provinsi_tahun_peta
tabel_eb <- peta_eb %>%
sf::st_drop_geometry() %>%
dplyr::select(
wilayah, tahun, kasus, populasi_berisiko,
incidence_rate, expected_cases, sir,
eb_incidence_rate, eb_sir,
kategori_sir, interpretasi_sir
) %>%
dplyr::arrange(dplyr::desc(eb_sir))
# D. Model Negative Binomial
data_model <- hasil_ukuran %>%
dplyr::mutate(
wilayah_f = factor(wilayah),
tahun_f = factor(tahun),
log_populasi = log(populasi_berisiko),
kepadatan_z = as.numeric(scale(kepadatan)),
kemiskinan_z = as.numeric(scale(kemiskinan)),
curah_hujan_z = as.numeric(scale(curah_hujan)),
kelembapan_z = as.numeric(scale(kelembapan)),
suhu_z = as.numeric(scale(suhu)),
sanitasi_z = as.numeric(scale(sanitasi)),
fasilitas_kesehatan_z = as.numeric(scale(fasilitas_kesehatan)),
tenaga_kesehatan_z = as.numeric(scale(tenaga_kesehatan))
)
formula_model <- kasus ~ wilayah_f + tahun_f +
kepadatan_z + kemiskinan_z + curah_hujan_z + kelembapan_z + suhu_z +
sanitasi_z + fasilitas_kesehatan_z + tenaga_kesehatan_z +
offset(log_populasi)
model_pois_pembanding <- stats::glm(
formula_model,
family = stats::poisson(link = "log"),
data = data_model
)
nilai_overdispersion_poisson <- sum(stats::residuals(model_pois_pembanding, type = "pearson")^2) /
stats::df.residual(model_pois_pembanding)
model_nb <- MASS::glm.nb(
formula_model,
data = data_model,
control = stats::glm.control(maxit = 100)
)
ringkasan_model_nb <- broom::tidy(model_nb) %>%
dplyr::mutate(
rate_ratio = exp(estimate),
ci_lower = exp(estimate - 1.96 * std.error),
ci_upper = exp(estimate + 1.96 * std.error),
interpretasi = dplyr::case_when(
rate_ratio > 1 ~ "Rate ratio > 1",
rate_ratio < 1 ~ "Rate ratio < 1",
TRUE ~ "Rate ratio = 1"
)
)
rate_referensi_semua_tahun <- sum(data_model$kasus, na.rm = TRUE) /
sum(data_model$populasi_berisiko, na.rm = TRUE)
hasil_model_nb_spasiotemporal <- data_model %>%
dplyr::mutate(
prediksi_kasus_nb = as.numeric(stats::predict(model_nb, type = "response")),
prediksi_rate_nb = prediksi_kasus_nb / populasi_berisiko,
prediksi_incidence_rate_nb = prediksi_rate_nb * 100000,
rr_model_nb = prediksi_rate_nb / rate_referensi_semua_tahun
) %>%
dplyr::select(
wilayah, tahun, kasus, populasi_berisiko,
incidence_rate, expected_cases, sir,
prediksi_kasus_nb, prediksi_incidence_rate_nb, rr_model_nb,
kepadatan, kemiskinan, curah_hujan, kelembapan, suhu,
sanitasi, fasilitas_kesehatan, tenaga_kesehatan
)
prioritas_wilayah_nb <- hasil_model_nb_spasiotemporal %>%
dplyr::filter(tahun == tahun_peta) %>%
dplyr::left_join(
tabel_eb %>% dplyr::select(wilayah, eb_incidence_rate, eb_sir),
by = "wilayah"
) %>%
dplyr::mutate(
skor_prioritas = dplyr::coalesce(eb_sir, sir) + rr_model_nb,
kriteria_prioritas = dplyr::case_when(
sir > 1 & eb_sir > 1 & rr_model_nb > 1 ~ "Prioritas tinggi",
sir > 1 & eb_sir > 1 ~ "Prioritas",
sir > 1 | eb_sir > 1 | rr_model_nb > 1 ~ "Perlu monitoring",
TRUE ~ "Relatif lebih rendah"
)
) %>%
dplyr::arrange(dplyr::desc(skor_prioritas)) %>%
dplyr::mutate(peringkat_prioritas = dplyr::row_number()) %>%
dplyr::select(
peringkat_prioritas, wilayah, tahun, kasus, populasi_berisiko,
incidence_rate, sir, eb_incidence_rate, eb_sir,
prediksi_incidence_rate_nb, rr_model_nb,
skor_prioritas, kriteria_prioritas
)
top5_wilayah <- prioritas_wilayah_nb %>%
dplyr::slice_head(n = 5) %>%
dplyr::pull(wilayah) %>%
stringr::str_to_title() %>%
paste(collapse = ", ")4.1 Analisis Deskriptif
Analisis deskriptif dilakukan untuk menggambarkan perkembangan kasus diare dan populasi berisiko di Provinsi Jawa Barat selama tahun 2020–2024. Ringkasan kasus dan populasi menurut tahun disajikan pada Tabel 1.
tabel_ringkasan_tahun <- ringkasan_tahun %>%
dplyr::mutate(
`Total Kasus` = fmt_angka(total_kasus, 0),
`Total Populasi` = fmt_angka(total_populasi, 0),
`IR Provinsi per 100.000` = fmt_angka(incidence_rate_provinsi, 2),
`Rata-rata Kasus per Wilayah` = fmt_angka(rata_rata_kasus_wilayah, 0)
) %>%
dplyr::select(
Tahun = tahun,
`Total Kasus`,
`Total Populasi`,
`IR Provinsi per 100.000`,
`Rata-rata Kasus per Wilayah`
)
knitr::kable(tabel_ringkasan_tahun, caption = "Ringkasan kasus diare dan populasi berisiko menurut tahun.")| Tahun | Total Kasus | Total Populasi | IR Provinsi per 100.000 | Rata-rata Kasus per Wilayah |
|---|---|---|---|---|
| 2020 | 1.305.102 | 48.274.400 | 2.703,51 | 48.337 |
| 2021 | 1.353.480 | 48.782.500 | 2.774,52 | 50.129 |
| 2022 | 1.367.256 | 49.404.000 | 2.767,50 | 50.639 |
| 2023 | 1.346.230 | 49.860.200 | 2.700,01 | 49.860 |
| 2024 | 672.287 | 50.345.200 | 1.335,35 | 24.900 |
Interpretasi Tabel 1. Tabel ini menunjukkan bahwa total kasus diare di Jawa Barat relatif tinggi pada periode 2020–2023, kemudian menurun tajam pada tahun 2024. Nilai incidence rate provinsi juga menurun pada tahun 2024, sehingga penurunan kasus tidak hanya terjadi secara absolut, tetapi juga setelah mempertimbangkan jumlah populasi berisiko.
ggplot2::ggplot(ringkasan_tahun, ggplot2::aes(x = tahun, y = total_kasus)) +
ggplot2::geom_line(linewidth = 1.1, color = "#2563eb") +
ggplot2::geom_point(size = 3, color = "#1d4ed8") +
ggplot2::scale_y_continuous(labels = scales::comma) +
ggplot2::scale_x_continuous(breaks = sort(unique(ringkasan_tahun$tahun))) +
ggplot2::labs(
title = "Tren Total Kasus Diare Menurut Tahun",
subtitle = "Provinsi Jawa Barat berdasarkan penjumlahan seluruh kabupaten/kota",
x = "Tahun",
y = "Total Kasus Diare"
) +
tema_publikasi()
Gambar 1. Tren total kasus diare di Jawa Barat tahun 2020–2024.
Interpretasi Gambar 1. Grafik menunjukkan bahwa total kasus diare relatif tinggi pada tahun 2020–2023, kemudian menurun tajam pada tahun 2024. Penurunan ini dapat menunjukkan perbaikan situasi, tetapi tetap perlu ditafsirkan hati-hati karena dapat dipengaruhi oleh perubahan pencatatan, pelaporan, atau akses layanan kesehatan.
4.2 Ukuran Epidemiologi
Ukuran epidemiologi yang dihitung meliputi incidence rate, expected cases, dan SIR. Tabel 2 menampilkan 10 wilayah dengan SIR tertinggi pada tahun 2024.
tabel_sir_10 <- prioritas_sir_tahun_terakhir %>%
dplyr::slice_head(n = 10) %>%
dplyr::mutate(
Wilayah = stringr::str_to_title(wilayah),
Kasus = fmt_angka(kasus, 0),
`Populasi Berisiko` = fmt_angka(populasi_berisiko, 0),
`Incidence Rate` = fmt_angka(incidence_rate, 2),
`Expected Cases` = fmt_angka(expected_cases, 2),
SIR = fmt_angka(sir, 2)
) %>%
dplyr::select(
Peringkat = peringkat_sir,
Wilayah,
Kasus,
`Populasi Berisiko`,
`Incidence Rate`,
`Expected Cases`,
SIR,
Kategori = kategori_sir
)
knitr::kable(tabel_sir_10, caption = paste("Sepuluh wilayah dengan SIR diare tertinggi tahun", tahun_peta, "."))| Peringkat | Wilayah | Kasus | Populasi Berisiko | Incidence Rate | Expected Cases | SIR | Kategori |
|---|---|---|---|---|---|---|---|
| 1 | Kota Cirebon | 15.271 | 344.900 | 4.427,66 | 4.605,64 | 3,32 | Sangat tinggi |
| 2 | Kota Bogor | 34.656 | 1.078.400 | 3.213,65 | 14.400,47 | 2,41 | Sangat tinggi |
| 3 | Kota Sukabumi | 11.385 | 365.700 | 3.113,21 | 4.883,39 | 2,33 | Sangat tinggi |
| 4 | Kabupaten Majalengka | 32.382 | 1.352.500 | 2.394,23 | 18.060,67 | 1,79 | Tinggi |
| 5 | Kota Tasikmalaya | 17.845 | 750.700 | 2.377,11 | 10.024,51 | 1,78 | Tinggi |
| 6 | Kota Banjar | 4.096 | 209.800 | 1.952,34 | 2.801,57 | 1,46 | Di atas rata-rata |
| 7 | Kabupaten Cirebon | 46.223 | 2.388.000 | 1.935,64 | 31.888,27 | 1,45 | Di atas rata-rata |
| 8 | Kabupaten Purwakarta | 20.193 | 1.050.300 | 1.922,59 | 14.025,23 | 1,44 | Di atas rata-rata |
| 9 | Kabupaten Pangandaran | 8.003 | 434.100 | 1.843,58 | 5.796,77 | 1,38 | Di atas rata-rata |
| 10 | Kabupaten Subang | 28.398 | 1.663.200 | 1.707,43 | 22.209,62 | 1,28 | Di atas rata-rata |
Interpretasi Tabel 2. Tabel ini memperlihatkan 10 wilayah dengan risiko relatif diare tertinggi berdasarkan SIR pada tahun 2024. Kota Cirebon, Kota Bogor, dan Kota Sukabumi menjadi tiga wilayah teratas. Seluruh wilayah dalam tabel memiliki SIR lebih dari 1, yang berarti jumlah kasus teramati lebih tinggi daripada jumlah kasus yang diharapkan berdasarkan rata-rata provinsi.
prioritas_sir_tahun_terakhir %>%
dplyr::slice_head(n = 10) %>%
dplyr::mutate(wilayah = forcats::fct_reorder(stringr::str_to_title(wilayah), sir)) %>%
ggplot2::ggplot(ggplot2::aes(x = wilayah, y = sir)) +
ggplot2::geom_col(fill = "#93c5fd", color = "#1d4ed8", linewidth = 0.25) +
ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "#334155") +
ggplot2::coord_flip() +
ggplot2::labs(
title = paste("10 Wilayah dengan SIR Diare Tertinggi Tahun", tahun_peta),
subtitle = "Garis putus-putus menunjukkan SIR = 1 atau rata-rata provinsi",
x = "Wilayah",
y = "SIR"
) +
tema_publikasi()
Gambar 2. Sepuluh wilayah dengan SIR diare tertinggi tahun 2024.
Interpretasi Gambar 2. Semua wilayah dalam 10 besar memiliki SIR di atas 1, sehingga berisiko lebih tinggi dibandingkan rata-rata provinsi. Kota Cirebon merupakan wilayah paling menonjol karena nilai SIR-nya paling tinggi, diikuti Kota Bogor dan Kota Sukabumi.
4.3 Disease Mapping
Hasil pencocokan data Excel dan GeoJSON menunjukkan 0 wilayah yang gagal dicocokkan. Dengan demikian, seluruh kabupaten/kota pada data dapat dipetakan.
ggplot2::ggplot(peta_tahun) +
ggplot2::geom_sf(ggplot2::aes(fill = incidence_rate), color = "white", linewidth = 0.25) +
ggplot2::scale_fill_viridis_c(option = "viridis", na.value = "grey90", labels = scales::comma) +
ggplot2::labs(
title = paste("Peta Crude Incidence Rate Diare Tahun", tahun_peta),
subtitle = "Incidence rate per 100.000 penduduk",
fill = "IR per 100.000"
) +
tema_publikasi()
Gambar 3. Peta crude incidence rate diare tahun 2024.
Interpretasi Gambar 3. Peta crude incidence rate menunjukkan bahwa angka kejadian diare per 100.000 penduduk tidak tersebar merata antarwilayah. Wilayah dengan warna lebih pekat memiliki angka kejadian lebih tinggi dan memerlukan perhatian dalam monitoring kesehatan masyarakat.
ggplot2::ggplot(peta_tahun) +
ggplot2::geom_sf(ggplot2::aes(fill = sir), color = "white", linewidth = 0.25) +
ggplot2::scale_fill_viridis_c(option = "magma", na.value = "grey90") +
ggplot2::labs(
title = paste("Peta Crude SIR Diare Tahun", tahun_peta),
subtitle = "SIR > 1 menunjukkan risiko lebih tinggi daripada rata-rata provinsi",
fill = "SIR"
) +
tema_publikasi()
Gambar 4. Peta crude SIR diare tahun 2024.
Interpretasi Gambar 4. Peta crude SIR menunjukkan wilayah yang memiliki risiko relatif lebih tinggi dibandingkan rata-rata Provinsi Jawa Barat. Wilayah yang menonjol pada peta ini konsisten dengan tabel SIR tertinggi, terutama Kota Cirebon, Kota Bogor, Kota Sukabumi, Kota Tasikmalaya, dan Kabupaten Majalengka.
ggplot2::ggplot(peta_eb) +
ggplot2::geom_sf(ggplot2::aes(fill = eb_sir), color = "white", linewidth = 0.25) +
ggplot2::scale_fill_viridis_c(option = "plasma", na.value = "grey90") +
ggplot2::labs(
title = paste("Peta SIR Diare Hasil Empirical Bayes Tahun", tahun_peta),
subtitle = "Estimasi risiko distabilkan menggunakan smoothing Empirical Bayes",
fill = "EB SIR"
) +
tema_publikasi()
Gambar 5. Peta SIR diare hasil Empirical Bayes tahun 2024.
Interpretasi Gambar 5. Peta EB SIR memperhalus estimasi risiko sehingga nilai yang ekstrem akibat fluktuasi populasi atau kasus menjadi lebih stabil. Wilayah yang tetap tinggi setelah smoothing dapat dianggap memiliki indikasi risiko yang lebih kuat.
4.4 Model Spasio-Temporal Negative Binomial
Model Poisson digunakan sebagai pembanding untuk memeriksa overdispersion. Nilai overdispersion Poisson yang diperoleh adalah 1349.19, jauh lebih besar dari 1,5. Oleh karena itu, model utama yang digunakan adalah Negative Binomial.
tabel_model <- tibble::tribble(
~Komponen, ~Hasil,
"Model pembanding", "Poisson",
"Nilai overdispersion Poisson", fmt_angka(nilai_overdispersion_poisson, 2),
"Model utama", "Negative Binomial",
"Efek wilayah", "Ada",
"Efek waktu/tahun", "Ada",
"Offset", "Log populasi berisiko",
"Kovariat", "Kepadatan, kemiskinan, curah hujan, kelembapan, suhu, sanitasi, fasilitas kesehatan, tenaga kesehatan"
)
knitr::kable(tabel_model, caption = "Ringkasan model spasio-temporal.")| Komponen | Hasil |
|---|---|
| Model pembanding | Poisson |
| Nilai overdispersion Poisson | 1.349,19 |
| Model utama | Negative Binomial |
| Efek wilayah | Ada |
| Efek waktu/tahun | Ada |
| Offset | Log populasi berisiko |
| Kovariat | Kepadatan, kemiskinan, curah hujan, kelembapan, suhu, sanitasi, fasilitas kesehatan, tenaga kesehatan |
Interpretasi Tabel 3. Tabel ini menunjukkan bahwa model Poisson hanya digunakan sebagai pembanding untuk memeriksa overdispersion. Nilai overdispersion yang tinggi menunjukkan bahwa variasi kasus diare lebih besar daripada asumsi Poisson, sehingga model Negative Binomial dipilih sebagai model utama spasio-temporal.
model_terms_penting <- ringkasan_model_nb %>%
dplyr::filter(
stringr::str_detect(term, "^tahun_f") |
term %in% c(
"kepadatan_z", "kemiskinan_z", "curah_hujan_z", "kelembapan_z",
"suhu_z", "sanitasi_z", "fasilitas_kesehatan_z", "tenaga_kesehatan_z"
)
) %>%
dplyr::mutate(
Variabel = dplyr::case_when(
stringr::str_detect(term, "^tahun_f") ~ stringr::str_replace(term, "tahun_f", "Tahun "),
term == "kepadatan_z" ~ "Kepadatan penduduk terstandar",
term == "kemiskinan_z" ~ "Kemiskinan terstandar",
term == "curah_hujan_z" ~ "Curah hujan terstandar",
term == "kelembapan_z" ~ "Kelembapan terstandar",
term == "suhu_z" ~ "Suhu terstandar",
term == "sanitasi_z" ~ "Sanitasi terstandar",
term == "fasilitas_kesehatan_z" ~ "Fasilitas kesehatan terstandar",
term == "tenaga_kesehatan_z" ~ "Tenaga kesehatan terstandar",
TRUE ~ term
),
`Rate Ratio` = fmt_angka(rate_ratio, 2),
`CI 95%` = paste0(fmt_angka(ci_lower, 2), " -- ", fmt_angka(ci_upper, 2)),
`p-value` = ifelse(p.value < 0.001, "<0,001", fmt_angka(p.value, 4)),
Interpretasi = dplyr::case_when(
rate_ratio > 1 ~ "Berasosiasi dengan peningkatan rate",
rate_ratio < 1 ~ "Berasosiasi dengan penurunan rate",
TRUE ~ "Tidak berubah"
)
) %>%
dplyr::select(Variabel, `Rate Ratio`, `CI 95%`, `p-value`, Interpretasi)
knitr::kable(model_terms_penting, caption = "Hasil utama model Negative Binomial.")| Variabel | Rate Ratio | CI 95% | p-value | Interpretasi |
|---|---|---|---|---|
| Tahun 2021 | 1,22 | 1,04 – 1,43 | 0,0121 | Berasosiasi dengan peningkatan rate |
| Tahun 2022 | 1,14 | 0,99 – 1,31 | 0,0631 | Berasosiasi dengan peningkatan rate |
| Tahun 2023 | 0,97 | 0,77 – 1,23 | 0,8167 | Berasosiasi dengan penurunan rate |
| Tahun 2024 | 0,52 | 0,43 – 0,64 | <0,001 | Berasosiasi dengan penurunan rate |
| Kepadatan penduduk terstandar | 1,94 | 0,61 – 6,13 | 0,2611 | Berasosiasi dengan peningkatan rate |
| Kemiskinan terstandar | 0,44 | 0,27 – 0,72 | 0,0010 | Berasosiasi dengan penurunan rate |
| Curah hujan terstandar | 0,95 | 0,86 – 1,04 | 0,2634 | Berasosiasi dengan penurunan rate |
| Kelembapan terstandar | 1,12 | 0,88 – 1,42 | 0,3652 | Berasosiasi dengan peningkatan rate |
| Suhu terstandar | 1,44 | 0,83 – 2,48 | 0,1932 | Berasosiasi dengan peningkatan rate |
| Sanitasi terstandar | 0,99 | 0,94 – 1,05 | 0,8166 | Berasosiasi dengan penurunan rate |
| Fasilitas kesehatan terstandar | 0,35 | 0,15 – 0,78 | 0,0105 | Berasosiasi dengan penurunan rate |
| Tenaga kesehatan terstandar | 0,91 | 0,84 – 0,99 | 0,0281 | Berasosiasi dengan penurunan rate |
Interpretasi Tabel 4. Tabel ini menyajikan rate ratio dari model Negative Binomial. Nilai rate ratio lebih dari 1 menunjukkan peningkatan rate kasus, sedangkan nilai kurang dari 1 menunjukkan penurunan rate kasus setelah mempertimbangkan efek wilayah, efek tahun, kovariat, dan populasi berisiko. Interpretasi model perlu dilakukan hati-hati karena hasilnya merupakan hubungan simultan dalam model, bukan hubungan sebab-akibat langsung.
hasil_ukuran %>%
dplyr::mutate(
tahun = factor(tahun),
wilayah = forcats::fct_reorder(stringr::str_to_title(wilayah), sir, .fun = mean, .desc = FALSE)
) %>%
ggplot2::ggplot(ggplot2::aes(x = tahun, y = wilayah, fill = sir)) +
ggplot2::geom_tile(color = "white") +
ggplot2::scale_fill_viridis_c(option = "magma") +
ggplot2::labs(
title = "Perubahan Risiko Diare Menurut Wilayah dan Tahun",
subtitle = "Berdasarkan crude SIR",
x = "Tahun",
y = "Wilayah",
fill = "SIR"
) +
tema_publikasi() +
ggplot2::theme(axis.text.y = ggplot2::element_text(size = 8))
Gambar 6. Heatmap perubahan risiko diare menurut wilayah dan tahun berdasarkan crude SIR.
Interpretasi Gambar 6. Heatmap crude SIR memperlihatkan perubahan risiko relatif setiap kabupaten/kota dari tahun 2020 sampai 2024. Warna yang lebih kuat menunjukkan nilai SIR yang lebih tinggi, sehingga memudahkan identifikasi wilayah dan tahun dengan risiko relatif menonjol.
hasil_model_nb_spasiotemporal %>%
dplyr::mutate(
tahun = factor(tahun),
wilayah = forcats::fct_reorder(stringr::str_to_title(wilayah), rr_model_nb, .fun = mean, .desc = FALSE)
) %>%
ggplot2::ggplot(ggplot2::aes(x = tahun, y = wilayah, fill = rr_model_nb)) +
ggplot2::geom_tile(color = "white") +
ggplot2::scale_fill_viridis_c(option = "plasma") +
ggplot2::labs(
title = "Perubahan Risiko Diare Menurut Wilayah dan Tahun",
subtitle = "Berdasarkan fitted relative risk model Negative Binomial",
x = "Tahun",
y = "Wilayah",
fill = "RR NB"
) +
tema_publikasi() +
ggplot2::theme(axis.text.y = ggplot2::element_text(size = 8))
Gambar 7. Heatmap perubahan risiko berdasarkan model Negative Binomial.
Interpretasi Gambar 7. Heatmap RR Negative Binomial menunjukkan pola risiko hasil model setelah mempertimbangkan efek wilayah, efek tahun, kovariat, dan populasi berisiko. Visualisasi ini melengkapi heatmap SIR karena risiko tidak hanya dilihat dari kasus teramati, tetapi juga dari prediksi model.
peta_model_nb_tahun <- peta %>%
dplyr::left_join(
hasil_model_nb_spasiotemporal %>% dplyr::filter(tahun == tahun_peta),
by = "wilayah"
)
ggplot2::ggplot(peta_model_nb_tahun) +
ggplot2::geom_sf(ggplot2::aes(fill = rr_model_nb), color = "white", linewidth = 0.25) +
ggplot2::scale_fill_viridis_c(option = "plasma", na.value = "grey90") +
ggplot2::labs(
title = paste("Peta Relative Risk Model Negative Binomial Tahun", tahun_peta),
subtitle = "Model dengan efek wilayah, efek tahun, kovariat, dan offset populasi",
fill = "RR NB"
) +
tema_publikasi()
Gambar 8. Peta relative risk model Negative Binomial diare tahun 2024.
Interpretasi Gambar 8. Peta RR Negative Binomial menggambarkan risiko relatif berdasarkan hasil model pada tahun 2024. Peta ini berguna untuk melihat wilayah yang tetap memiliki risiko relatif lebih tinggi setelah perbedaan populasi dan kovariat diperhitungkan.
4.5 Wilayah Prioritas dan Implikasi Epidemiologis
Wilayah prioritas ditentukan berdasarkan kombinasi crude SIR, EB SIR, dan RR model Negative Binomial. Pendekatan gabungan ini digunakan agar penentuan prioritas tidak hanya bergantung pada risiko kasar.
tabel_prioritas <- prioritas_wilayah_nb %>%
dplyr::slice_head(n = 10) %>%
dplyr::mutate(
Wilayah = stringr::str_to_title(wilayah),
Kasus = fmt_angka(kasus, 0),
`Incidence Rate` = fmt_angka(incidence_rate, 2),
SIR = fmt_angka(sir, 2),
`EB SIR` = fmt_angka(eb_sir, 2),
`RR NB` = fmt_angka(rr_model_nb, 2)
) %>%
dplyr::select(
Peringkat = peringkat_prioritas,
Wilayah,
Kasus,
`Incidence Rate`,
SIR,
`EB SIR`,
`RR NB`,
Kriteria = kriteria_prioritas
)
knitr::kable(tabel_prioritas, caption = paste("Sepuluh wilayah prioritas tahun", tahun_peta, "berdasarkan SIR, EB SIR, dan model Negative Binomial."))| Peringkat | Wilayah | Kasus | Incidence Rate | SIR | EB SIR | RR NB | Kriteria |
|---|---|---|---|---|---|---|---|
| 1 | Kota Cirebon | 15.271 | 4.427,66 | 3,32 | 3,31 | 0,90 | Prioritas |
| 2 | Kota Bogor | 34.656 | 3.213,65 | 2,41 | 2,41 | 0,82 | Prioritas |
| 3 | Kota Sukabumi | 11.385 | 3.113,21 | 2,33 | 2,33 | 0,80 | Prioritas |
| 4 | Kota Tasikmalaya | 17.845 | 2.377,11 | 1,78 | 1,78 | 0,88 | Prioritas |
| 5 | Kabupaten Majalengka | 32.382 | 2.394,23 | 1,79 | 1,79 | 0,76 | Prioritas |
| 6 | Kota Banjar | 4.096 | 1.952,34 | 1,46 | 1,46 | 0,71 | Prioritas |
| 7 | Kabupaten Cirebon | 46.223 | 1.935,64 | 1,45 | 1,45 | 0,70 | Prioritas |
| 8 | Kabupaten Pangandaran | 8.003 | 1.843,58 | 1,38 | 1,38 | 0,67 | Prioritas |
| 9 | Kabupaten Purwakarta | 20.193 | 1.922,59 | 1,44 | 1,44 | 0,59 | Prioritas |
| 10 | Kabupaten Sumedang | 18.361 | 1.546,71 | 1,16 | 1,16 | 0,75 | Prioritas |
Interpretasi Tabel 5. Tabel ini menunjukkan wilayah prioritas berdasarkan kombinasi crude SIR, EB SIR, dan RR model Negative Binomial. Lima wilayah prioritas utama pada tahun 2024 adalah Kota Cirebon, Kota Bogor, Kota Sukabumi, Kota Tasikmalaya, Kabupaten Majalengka. Wilayah tersebut perlu menjadi fokus monitoring karena menunjukkan risiko relatif yang tinggi, terutama berdasarkan SIR dan EB SIR.
Implikasi epidemiologis dari hasil ini adalah perlunya intervensi yang lebih terarah secara spasial. Wilayah prioritas dapat menjadi sasaran penguatan surveilans diare, peningkatan akses sanitasi, edukasi perilaku hidup bersih dan sehat, serta pemantauan faktor lingkungan seperti curah hujan, kelembapan, dan suhu.
5 Kesimpulan
Berdasarkan hasil analisis, kasus diare di Provinsi Jawa Barat tahun 2020–2024 menunjukkan variasi menurut wilayah dan waktu. Secara deskriptif, total kasus diare relatif tinggi pada periode 2020–2023, kemudian mengalami penurunan pada tahun 2024. Penurunan kasus pada tahun 2024 juga terlihat dari nilai incidence rate provinsi yang lebih rendah dibandingkan tahun-tahun sebelumnya.
Hasil perhitungan ukuran epidemiologi menunjukkan bahwa beberapa kabupaten/kota memiliki nilai SIR lebih dari 1. Pada tahun 2024, wilayah dengan risiko tertinggi berdasarkan hasil analisis adalah Kota Cirebon, Kota Bogor, Kota Sukabumi, Kota Tasikmalaya, Kabupaten Majalengka. Wilayah-wilayah tersebut menjadi prioritas karena menunjukkan risiko diare yang relatif lebih tinggi berdasarkan incidence rate, SIR, smoothing Empirical Bayes, dan hasil model spasio-temporal Negative Binomial.
Hasil disease mapping menunjukkan bahwa pemetaan crude incidence rate, crude SIR, dan EB SIR dapat membantu menggambarkan pola spasial risiko diare secara lebih jelas. Model spasio-temporal Negative Binomial digunakan karena data kasus diare menunjukkan overdispersion, yaitu variasi kasus yang lebih besar daripada asumsi model Poisson.
Keterbatasan analisis ini adalah data yang digunakan bersifat agregat wilayah, sehingga hasilnya tidak dapat digunakan untuk menyimpulkan risiko pada tingkat individu. Selain itu, hasil analisis bergantung pada kualitas pencatatan kasus, ketepatan data populasi, dan kelengkapan kovariat. Beberapa faktor penting seperti akses air minum aman, perilaku higiene, kejadian banjir, kualitas sanitasi rumah tangga, dan kualitas pelaporan fasilitas kesehatan belum seluruhnya tercakup dalam model.
Rekomendasi utama adalah memperkuat surveilans diare pada wilayah prioritas, meningkatkan akses sanitasi, memperkuat edukasi perilaku hidup bersih dan sehat, memantau faktor lingkungan, serta meningkatkan kualitas deteksi dan pelaporan kasus. Analisis serupa sebaiknya dilakukan secara berkala dengan data tahunan atau bulanan agar perubahan risiko dapat dipantau lebih dini.
6 Daftar Pustaka
Adria, B. P., Wulandari, R. A., & Kurniasari, F. (2023). Spatial analysis of diarrhea incidences, environmental influences, and behavioral factors in an ecological study. Media Kesehatan Masyarakat Indonesia, 19(4), 185–192. https://scholarhub.unhas.ac.id/mkmi/vol19/iss4/2/
Badan Informasi Geospasial. (2024). Area batas wilayah administrasi kabupaten/kota. Geoservices BIG. https://geoservices.big.go.id/rbi/rest/services/BATASWILAYAH/BATAS_WILAYAH/MapServer/layers
Badan Pusat Statistik. (2024). Jumlah penduduk berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/proyeksi-penduduk-berdasarkan-kabupatenkota-di-jawa-barat
Coly, S., Garrido, M., Abrial, D., & Yao, A. F. (2021). Bayesian hierarchical models for disease mapping applied to contagious pathologies. PLOS ONE, 16(1), e0222898. https://doi.org/10.1371/journal.pone.0222898
Dinas Kesehatan Provinsi Jawa Barat. (2024). Jumlah kasus penyakit diare berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/en/dataset/jumlah-kasus-penyakit-diare-berdasarkan-kabupatenkota-di-jawa-barat
Dinas Kesehatan Provinsi Jawa Barat. (2024). Jumlah fasilitas kesehatan berdasarkan jenis di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/jumlah-fasilitas-kesehatan-berdasarkan-jenis-di-jawa-barat
Dinas Kesehatan Provinsi Jawa Barat. (2024). Jumlah tenaga kesehatan masyarakat berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/en/dataset/jumlah-tenaga-kesehatan-masyarakat-berdasarkan-kabupatenkota-di-jawa-barat
Dinas Kesehatan Provinsi Jawa Barat. (2025). Persentase keluarga dengan akses sanitasi layak (jamban sehat) berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/persentase-keluarga-dengan-akses-sanitasi-layak-jamban-sehat-berdasarkan-kabupatenkota-di-jawa-barat
Geremew, G., et al. (2024). Rainfall and temperature influences on childhood diarrhea and the effect modification role of water and sanitation conditions: A systematic review and meta-analysis. International Journal of Environmental Research and Public Health, 21(7), 823. https://doi.org/10.3390/ijerph21070823
Liu, S., Wu, Z., An, Q., Wu, J., Bai, J., Sun, W., Guo, L., & Gong, L. (2025). The burden and etiologies of diarrhea in Asia and its countries from 1990 to 2021 and the forecast to 2040: Analyses informed by the Global Burden of Disease Study 2021. Frontiers in Public Health. https://www.frontiersin.org/journals/public-health/articles/10.3389/fpubh.2025.1651315/full
Pemerintah Provinsi Jawa Barat. (2024). Open Data Jabar: Portal data terbuka Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/
Wahyudi, T., Yuniar, P., & Rahmaniati, M. (2024). Spatial autocorrelation of diarrhea cases in West Java Province in 2023. Media Publikasi Promosi Kesehatan Indonesia, 7(9), 2368–2376. https://doi.org/10.56338/mppki.v7i9.5973
7 Tautan Data dan GeoJSON
Data tabular dan GeoJSON yang digunakan untuk kebutuhan
reproduksibilitas analisis tersedia pada folder Google Drive
berikut:
https://drive.google.com/drive/folders/1n-LfSxZWyIeYisQTiex7atgh7qUCtFuV?usp=sharing