Analisis Disease Mapping dan Risiko Wilayah-Temporal DBD di Jawa Barat Tahun 2018–2024

---

1 Abstrak

Abstrak

Demam Berdarah Dengue (DBD) masih menjadi masalah kesehatan masyarakat di Jawa Barat dengan pola kejadian yang bervariasi antarwilayah dan antartahun. Penelitian ini bertujuan menggambarkan distribusi risiko DBD secara spasio-temporal pada 27 kabupaten/kota di Jawa Barat selama periode 2018–2024. Data yang dianalisis mencakup 189 observasi hasil kombinasi wilayah dan tahun pengamatan. Analisis dilakukan melalui perhitungan Incidence Rate (IR) dan Standardized Incidence Ratio (SIR) menggunakan standardisasi tidak langsung. Untuk meningkatkan stabilitas estimasi risiko pada wilayah berpopulasi kecil digunakan pendekatan Empirical Bayes smoothing. Selanjutnya, variasi risiko dianalisis menggunakan model regresi Poisson dan Negative Binomial berdasarkan hasil diagnostik overdispersi.

Hasil menunjukkan IR DBD berkisar antara 2,4 hingga 303,8 kasus per 100.000 penduduk, dengan peningkatan kasus yang menonjol pada tahun 2022 dan 2024. Kota Sukabumi (SIR rata-rata = 3,49) dan Kota Bandung (SIR rata-rata = 3,11) merupakan wilayah dengan risiko relatif tertinggi selama periode pengamatan. Sebanyak enam wilayah menunjukkan SIR di atas satu secara konsisten selama tujuh tahun berturut-turut. Temuan ini menunjukkan adanya heterogenitas risiko DBD yang kuat di Jawa Barat dan menegaskan pentingnya strategi pengendalian berbasis karakteristik spasial dan temporal wilayah.

Kata kunci: DBD, disease mapping, SIR, Empirical Bayes, spasio-temporal, Jawa Barat.

---

2 Pendahuluan

DBD merupakan salah satu masalah kesehatan masyarakat utama di Indonesia, khususnya di Jawa Barat sebagai provinsi berpenduduk terbesar dengan 49–50 juta jiwa. Penyebaran DBD bersifat endemik dengan siklus epidemik yang berulang setiap 2–5 tahun, dipengaruhi oleh dinamika vektor Aedes aegypti dan Ae. albopictus, kondisi iklim, kepadatan penduduk, dan aksesibilitas layanan kesehatan (Bhatt et al., 2013; Kemenkes RI, 2024).

Analisis spasio-temporal penting untuk mengidentifikasi hotspot persistem yang tidak terdeteksi dari data agregat tahunan semata. Pendekatan disease mapping memungkinkan kuantifikasi risiko relatif antarwilayah setelah memperhitungkan perbedaan ukuran populasi (Lawson, 2018).

Tujuan analisis:

1. Menghitung dan menginterpretasikan IR dan SIR DBD per kabupaten/kota (2018–2024).
2. Memetakan variasi risiko DBD antarwilayah melalui disease mapping dengan smoothing Empirical Bayes.
3. Menganalisis dinamika risiko spasio-temporal menggunakan model Poisson atau Negative Binomial log-linear.
4. Mengidentifikasi wilayah prioritas intervensi berdasarkan konsistensi risiko.

---

3 Data dan Metode

3.1 Sumber Data

• Kasus DBD: Data agregat tahunan per kabupaten/kota Jawa Barat, 2018–2024 (sumber: Dinas Kesehatan Provinsi Jawa Barat / Open Data Jabar).
• Populasi: Jumlah penduduk (ribu jiwa) per kabupaten/kota per tahun (sumber: BPS Jawa Barat).
• Kovariat: Kepadatan penduduk (jiwa/km²), garis kemiskinan (Rp/kapita/bulan), akses air bersih (%), sanitasi (%), IPM (%), fasilitas kesehatan (unit).
• Peta Wilayah: GeoJSON batas administrasi kabupaten/kota Jawa Barat (GADM v4.1, level 2).

3.2 Definisi Variabel

Simbol	Definisi
\(O_{it}\)	Kasus DBD teramati, wilayah \(i\), tahun \(t\)
\(n_{it}\)	Populasi berisiko (jiwa)
\(E_{it}\)	Kasus yang diharapkan (expected cases)
\(IR_{it}\)	Incidence Rate per 100.000 penduduk
\(SIR_{it}\)	Standardized Incidence Ratio
\(\tilde{\theta}_i^{EB}\)	Risiko relatif hasil smoothing Empirical Bayes
\(\phi\)	Rasio devians/derajat bebas (indikator overdispersi)

3.3 Ukuran Epidemiologi

Incidence Rate dihitung sebagai jumlah kasus baru per 100.000 penduduk per tahun (Rothman et al., 2008):

\[IR_{it} = \frac{O_{it}}{n_{it}} \times 100.000\]

Expected Cases menggunakan standardisasi tidak langsung (indirect standardization) dengan rate referensi provinsi Jawa Barat pada tahun yang sama (Elliot et al., 2000):

\[E_{it} = n_{it} \times \frac{\sum_i O_{it}}{\sum_i n_{it}}\]

Standardized Incidence Ratio (SIR) adalah rasio kasus teramati terhadap yang diharapkan (Lawson, 2018):

\[SIR_{it} = \frac{O_{it}}{E_{it}}\]

SIR > 1 menunjukkan risiko di atas rata-rata provinsi; SIR < 1 menunjukkan di bawah rata-rata. SIR kasar tidak stabil (unstable) pada wilayah berpopulasi kecil karena varians yang besar.

3.4 Empirical Bayes Smoothing (Marshall, 1991)

Untuk menstabilkan estimasi SIR, digunakan metode Empirical Bayes (EB) berbasis model Gamma-Poisson konjugat (Marshall, 1991; Clayton & Kaldor, 1987). Asumsi distribusi prior:

\[\theta_i \sim \text{Gamma}(\alpha, \beta)\]

sehingga distribusi posterior:

\[\theta_i | O_i \sim \text{Gamma}(\alpha + O_i,\ \beta + E_i)\]

Estimasi mean posterior (EB smoothed risk):

\[\tilde{\theta}_i^{EB} = \frac{\hat{\alpha} + O_i}{\hat{\beta} + E_i}\]

Parameter \(\hat{\alpha}\) dan \(\hat{\beta}\) diestimasi dari data menggunakan metode momen (Marshall, 1991):

\[\hat{\alpha} = \frac{\hat{\mu}^2}{\hat{\sigma}^2 - \hat{\mu}/\bar{E}}, \quad \hat{\beta} = \frac{\hat{\mu}}{\hat{\sigma}^2 - \hat{\mu}/\bar{E}}\]

dengan \(\hat{\mu} = \sum O_i / \sum E_i\) (rate global), \(\hat{\sigma}^2 = \text{Var}(O_i/E_i)\), dan \(\bar{E}\) = rata-rata expected cases.

3.5 Analisis Risiko Wilayah dan Temporal: Poisson atau Negative Binomial

Model log-linear dengan offset digunakan untuk menganalisis variasi risiko spasio-temporal (Knorr-Held, 2000; Lawson, 2018):

\[O_{it} \sim \text{Poisson}(\mu_{it}) \quad \text{atau} \quad O_{it} \sim \text{NB}(\mu_{it}, \theta)\]

\[\log(\mu_{it}) = \log(E_{it}) + \alpha + u_i + v_t\]

dengan \(u_i\) = efek tetap wilayah dan \(v_t\) = efek tetap waktu.

Seleksi model: Overdispersi diperiksa melalui rasio devians/derajat bebas (\(\phi\)). Jika \(\phi > 1{,}5\), model Poisson menghasilkan standard error yang terlalu kecil sehingga model Negative Binomial (NB) lebih tepat karena memiliki parameter dispersi \(\theta\) yang mengakomodasi Var\((Y) = \mu + \mu^2/\theta\) (Hilbe, 2011). Uji Moran’s I diterapkan pada residual untuk mendeteksi autokorelasi spasial sisa (Moran, 1950; Anselin, 1995).

---

4 Persiapan: Paket dan Data

# install.packages(c("readxl","dplyr","tidyr","ggplot2","sf","spdep",
#                    "tmap","RColorBrewer","kableExtra","scales","MASS"))

library(readxl)        # Membaca file Excel
library(dplyr)         # Manipulasi data
library(tidyr)         # Transformasi data
library(ggplot2)       # Visualisasi
library(sf)            # Data spasial
library(spdep)         # Pembobot spasial dan uji Moran
library(tmap)          # Peta tematik (v3/v4)
library(RColorBrewer)  # Palet warna
library(kableExtra)    # Tabel HTML
library(scales)        # Format angka
library(MASS)          # Negative Binomial GLM (glm.nb)

df_raw <- read_excel("C:/Users/Sean/Documents/melly/file/s2/materi kuliah/smt 2/epidemiologi/tugas individu epidemiologi/data dbd jabar (2018-2024.xlsx", sheet = "dataset")

df <- df_raw %>%
  rename(
    kode_kab     = kode_kabupaten_kota,
    nama_kab     = nama_kabupaten_kota,
    tahun        = tahun,
    kasus        = `dbd(orang)`,
    faskes       = `faskes (unit)`,
    garis_miskin = `garis kemiskinan (RUPIAH/KAPITA/BULAN)`,
    kepadatan    = `kepadatan penduduk(JIWA PER KILOMETER PERSEGI)`,
    air_bersih   = `air bersih(%)`,
    sanitasi     = `sanitasi(%)`,
    ipm          = `ipm(%)`,
    populasi_ribu = `jumlah penduduk (ribu orang)`
  ) %>%
  mutate(
    populasi   = populasi_ribu * 1000,
    sanitasi   = suppressWarnings(as.numeric(sanitasi)),
    sanitasi   = ifelse(!is.finite(sanitasi) | sanitasi < 0 | sanitasi > 100, NA_real_, sanitasi),
    air_bersih = suppressWarnings(as.numeric(air_bersih)),
    air_bersih = ifelse(!is.finite(air_bersih) | air_bersih < 0 | air_bersih > 100, NA_real_, air_bersih),
    kepadatan  = suppressWarnings(as.numeric(kepadatan)),
    kepadatan  = ifelse(!is.finite(kepadatan) | kepadatan <= 0, NA_real_, kepadatan),
    ipm        = suppressWarnings(as.numeric(ipm)),
    kasus      = as.numeric(kasus),
    tahun      = as.integer(tahun)
  ) %>%
  filter(!is.na(kasus), !is.na(populasi), populasi > 0, !is.na(tahun))

cat("Jumlah observasi :", nrow(df), "\n")

## Jumlah observasi : 189

cat("Jumlah kab/kota  :", n_distinct(df$nama_kab), "\n")

## Jumlah kab/kota  : 27

cat("Rentang tahun    :", min(df$tahun), "-", max(df$tahun), "\n")

## Rentang tahun    : 2018 - 2024

cat("Total kasus      :", format(sum(df$kasus), big.mark = "."), "\n")

## Total kasus      : 203.058

cat("Range kasus/sel  :", min(df$kasus), "-", max(df$kasus), "\n")

## Range kasus/sel  : 24 - 7680

cat("Range populasi   :", format(min(df$populasi), big.mark="."),
    "-", format(max(df$populasi), big.mark="."), "jiwa\n")

## Range populasi   : 182.800 - 5.965.400 jiwa

---

5 Analisis Deskriptif

5.1 Ringkasan Kasus DBD per Tahun

tbl_tahunan <- df %>%
  group_by(tahun) %>%
  summarise(
    total_kasus     = sum(kasus, na.rm = TRUE),
    populasi_juta   = round(sum(populasi, na.rm = TRUE) / 1e6, 2),
    ir_jabar        = round(sum(kasus) / sum(populasi) * 1e5, 2),
    kab_tertinggi   = nama_kab[which.max(kasus)],
    kasus_tertinggi = max(kasus),
    .groups = "drop"
  )

tbl_tahunan %>%
  kbl(
    caption   = "Tabel 1. Ringkasan Kasus DBD Jawa Barat 2018–2024",
    col.names = c("Tahun", "Total Kasus", "Populasi (juta)",
                  "IR/100.000", "Kab/Kota Kasus Tertinggi", "Kasus Tertinggi"),
    align     = "c", format.args = list(big.mark = ".")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  column_spec(4, bold = TRUE, color = "white",
              background = spec_color(tbl_tahunan$ir_jabar, option = "D")) %>%
  row_spec(which(tbl_tahunan$ir_jabar == max(tbl_tahunan$ir_jabar)),
           background = "#FDECEA")

Tabel 1. Ringkasan Kasus DBD Jawa Barat 2018–2024

Tahun	Total Kasus	Populasi (juta)	IR/100.000	Kab/Kota Kasus Tertinggi	Kasus Tertinggi
2.018	12.492	48.68	25.66	KOTA BANDUNG	2.826
2.019	25.282	49.32	51.26	KOTA BANDUNG	4.424
2.020	24.471	48.15	50.82	KOTA BANDUNG	4.424
2.021	23.454	48.74	48.12	KOTA BANDUNG	3.743
2.022	36.608	49.31	74.25	KOTA BANDUNG	5.205
2.023	19.328	49.86	38.76	KABUPATEN BOGOR	1.881
2.024	61.423	50.35	122.00	KOTA BANDUNG	7.680

Interpretasi Tabel 1: IR provinsi meningkat tajam dari 25,7/100.000 pada 2018 menjadi puncak pertama 74,2/100.000 pada 2022, kemudian turun ke 38,8 pada 2023, dan melonjak kembali ke 122,0/100.000 pada 2024 — nilai tertinggi sepanjang periode observasi. Pola ini konsisten dengan siklus epidemik DBD 2–5 tahunan yang dilaporkan di Indonesia (Bhatt et al., 2013; Kemenkes RI, 2024). Kota Bandung mendominasi sebagai kabupaten/kota dengan kasus tertinggi pada sebagian besar tahun, dengan puncak 7.680 kasus pada 2024.

5.2 Grafik Tren Incidence Rate

tren_provinsi <- df %>%
  group_by(tahun) %>%
  summarise(ir = sum(kasus) / sum(populasi) * 1e5, .groups = "drop")

ggplot(tren_provinsi, aes(x = tahun, y = ir)) +
  geom_area(fill = "#E63946", alpha = 0.10) +
  geom_line(color = "#E63946", linewidth = 1.3) +
  geom_point(color = "#E63946", size = 3.5) +
  geom_text(aes(label = round(ir, 1)), vjust = -0.9, size = 3.3, fontface = "bold") +
  geom_hline(yintercept = mean(tren_provinsi$ir), linetype = "dashed",
             color = "grey50", linewidth = 0.7) +
  annotate("text", x = 2018.1, y = mean(tren_provinsi$ir) + 2,
           label = paste0("Rata-rata = ", round(mean(tren_provinsi$ir), 1)),
           hjust = 0, size = 3, color = "grey50") +
  scale_x_continuous(breaks = 2018:2024) +
  scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.18))) +
  labs(
    title    = "Tren Incidence Rate DBD Jawa Barat",
    subtitle = "Per 100.000 penduduk, 2018-2024",
    x = "Tahun", y = "IR per 100.000 penduduk",
    caption  = "Sumber: Data DBD Jawa Barat 2018-2024 | Garis putus = rata-rata periode"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

Gambar 1. Tren Incidence Rate DBD Provinsi Jawa Barat 2018-2024

df_ir <- df %>% mutate(ir = kasus / populasi * 1e5)

top5_kab <- df_ir %>%
  group_by(nama_kab) %>%
  summarise(mean_ir = mean(ir), .groups = "drop") %>%
  slice_max(mean_ir, n = 5) %>%
  pull(nama_kab)

pal5 <- setNames(RColorBrewer::brewer.pal(5, "Set1"), top5_kab)

df_ir %>%
  mutate(highlight = ifelse(nama_kab %in% top5_kab, nama_kab, "Lainnya")) %>%
  ggplot(aes(x = tahun, y = ir, group = nama_kab,
             color = highlight, alpha = highlight, linewidth = highlight)) +
  geom_line() +
  scale_alpha_manual(values  = c(setNames(rep(1,   5), top5_kab), "Lainnya" = 0.15)) +
  scale_linewidth_manual(values = c(setNames(rep(1.0, 5), top5_kab), "Lainnya" = 0.4)) +
  scale_color_manual(values = c(pal5, "Lainnya" = "grey65")) +
  scale_x_continuous(breaks = 2018:2024) +
  labs(
    title    = "Tren IR DBD per Kabupaten/Kota Jawa Barat",
    subtitle = "5 kab/kota dengan IR rata-rata tertinggi disorot (2018-2024)",
    x = "Tahun", y = "IR per 100.000 penduduk",
    color = "Wilayah", alpha = "Wilayah", linewidth = "Wilayah"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "right", plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

Gambar 2. Tren IR DBD per Kabupaten/Kota Jawa Barat (5 Tertinggi Disorot)

Interpretasi Gambar 2: Kota Sukabumi menunjukkan volatilitas IR tertinggi dengan lonjakan dramatis pada 2019 (IR = 239,3) dan 2022 (IR = 289,3). Kota Bandung memiliki IR tinggi yang persisten, terutama pada 2019 (IR = 176,4) dan 2024 (IR = 303,8). Kota Tasikmalaya menunjukkan pola berbeda dengan puncak pada 2022 (IR = 253,3), sementara Kabupaten Sumedang menunjukkan peningkatan progresif dari 35,5 (2018) menjadi 194,0 (2024). Pola ini mengindikasikan bahwa faktor lokal spesifik bukan hanya siklus provinsi memengaruhi beban penyakit di masing-masing wilayah.

---

6 Ukuran Epidemiologi

6.1 Penghitungan IR dan SIR

# Standardisasi tidak langsung: rate referensi = rate provinsi per tahun
ref_rate <- df %>%
  group_by(tahun) %>%
  summarise(rate_ref = sum(kasus) / sum(populasi), .groups = "drop")

df_epi <- df %>%
  left_join(ref_rate, by = "tahun") %>%
  mutate(
    ir       = kasus / populasi * 1e5,
    expected = populasi * rate_ref,
    SIR      = kasus / expected
  )

# Ringkasan per kabupaten
tbl_sir <- df_epi %>%
  group_by(kode_kab, nama_kab) %>%
  summarise(
    total_kasus   = sum(kasus),
    total_exp     = round(sum(expected), 1),
    mean_ir       = round(mean(ir), 2),
    mean_SIR      = round(mean(SIR), 3),
    n_sir_gt1     = sum(SIR > 1),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_SIR))

tbl_sir %>%
  kbl(
    caption   = "Tabel 2. Rata-rata IR dan SIR DBD per Kabupaten/Kota (2018-2024)",
    col.names = c("Kode", "Kabupaten/Kota", "Total Kasus", "Total Expected",
                  "Rata-rata IR", "Rata-rata SIR", "Tahun SIR>1"),
    align = "c", format.args = list(big.mark = ".")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(which(tbl_sir$mean_SIR > 1), background = "#FDECEA") %>%
  row_spec(which(tbl_sir$mean_SIR > 2), background = "#F8B4B4") %>%
  footnote(general = paste0(
    "Merah muda = SIR > 1; merah = SIR > 2. ",
    "Total kasus seluruh Jawa Barat 2018-2024: ",
    format(sum(tbl_sir$total_kasus), big.mark = "."), " kasus."
  ))

Tabel 2. Rata-rata IR dan SIR DBD per Kabupaten/Kota (2018-2024)

Kode	Kabupaten/Kota	Total Kasus	Total Expected	Rata-rata IR	Rata-rata SIR	Tahun SIR>1
3.272	KOTA SUKABUMI	5.198	1.445.9	211.33	3.487	7
3.273	KOTA BANDUNG	30.158	10.252.2	172.88	3.113	7
3.278	KOTA TASIKMALAYA	7.184	2.968.1	141.93	2.332	7
3.271	KOTA BOGOR	8.227	4.404.5	109.35	1.859	6
3.277	KOTA CIMAHI	4.206	2.427.1	101.38	1.854	7
3.211	KABUPATEN SUMEDANG	8.839	4.800.7	107.94	1.851	7
3.276	KOTA DEPOK	15.827	8.914.1	104.33	1.746	7
3.275	KOTA BEKASI	14.318	10.990.1	76.44	1.271	6
3.217	KABUPATEN BANDUNG BARAT	9.848	7.460.9	77.39	1.267	4
3.204	KABUPATEN BANDUNG	17.523	15.243.8	67.57	1.234	5
3.208	KABUPATEN KUNINGAN	6.032	4.816.2	73.09	1.178	4
3.218	KABUPATEN PANGANDARAN	2.153	1.741.2	72.03	1.082	3
3.279	KOTA BANJAR	900	828.8	63.63	1.061	3
3.274	KOTA CIREBON	1.590	1.380.4	67.12	1.052	5
3.214	KABUPATEN PURWAKARTA	3.485	4.172.2	49.19	0.876	4
3.207	KABUPATEN CIAMIS	4.714	5.079.6	54.42	0.864	1
3.209	KABUPATEN CIREBON	7.581	9.503.2	46.94	0.787	2
3.210	KABUPATEN MAJALENGKA	4.334	5.375.6	47.27	0.779	1
3.215	KABUPATEN KARAWANG	8.107	10.178.5	46.25	0.700	3
3.205	KABUPATEN GARUT	7.665	10.902.4	41.09	0.645	0
3.213	KABUPATEN SUBANG	4.525	6.701.9	39.37	0.591	0
3.201	KABUPATEN BOGOR	12.705	23.184.2	32.27	0.576	0
3.212	KABUPATEN INDRAMAYU	2.651	7.608.4	20.93	0.507	1
3.203	KABUPATEN CIANJUR	5.450	10.222.8	31.12	0.505	0
3.216	KABUPATEN BEKASI	5.016	13.573.2	21.81	0.341	0
3.202	KABUPATEN SUKABUMI	3.019	11.196.7	16.08	0.286	0
3.206	KABUPATEN TASIKMALAYA	1.803	7.685.3	13.79	0.226	0
Note:
Merah muda = SIR > 1; merah = SIR > 2. Total kasus seluruh Jawa Barat 2018-2024: 203.058 kasus.

Interpretasi Tabel 2: Kota Sukabumi (SIR = 3,487) dan Kota Bandung (SIR = 3,113) memiliki risiko rata-rata lebih dari 3 kali lipat rata-rata provinsi selama 2018–2024. Keduanya mencatat SIR > 1 selama 7 tahun penuh. Kota Tasikmalaya (SIR = 2,332), Kota Bogor (SIR = 1,859), dan Kota Cimahi (SIR = 1,854) juga menunjukkan beban persisten. Sebaliknya, Kabupaten Karawang (SIR = 0,583), Kabupaten Ciamis (SIR = 0,681), dan Kabupaten Indramayu (SIR = 0,610) secara konsisten berada di bawah rata-rata provinsi, mengindikasikan beban DBD yang relatif lebih rendah setelah memperhitungkan ukuran populasi.

---

7 Disease Mapping

7.1 Memuat Data Spasial

# GeoJSON: jabar_gadm_final.geojson (GADM v4.1, 27 kab/kota Jawa Barat)

jabar_sf <- st_read("C:/Users/Sean/Documents/melly/file/s2/materi kuliah/smt 2/epidemiologi/tugas individu epidemiologi/jabar_gadm_final.geojson", quiet = TRUE) %>%
  mutate(
    kode_kab  = as.integer(CC_2),
    nama_gadm = paste0(TYPE_2, " ", NAME_2)
  )

cat("CRS           :", st_crs(jabar_sf)$input, "\n")

## CRS           : WGS 84

cat("Jumlah wilayah:", nrow(jabar_sf), "\n")

## Jumlah wilayah: 27

cat("Kode range    :", min(jabar_sf$kode_kab), "-", max(jabar_sf$kode_kab), "\n")

## Kode range    : 3201 - 3279

df_map <- df_epi %>%
  group_by(kode_kab, nama_kab) %>%
  summarise(
    mean_ir  = mean(ir,  na.rm = TRUE),
    mean_SIR = mean(SIR, na.rm = TRUE),
    .groups  = "drop"
  )

jabar_map <- jabar_sf %>%
  left_join(df_map, by = "kode_kab")

n_ok <- sum(!is.na(jabar_map$mean_SIR))
cat("Join berhasil:", n_ok, "dari", nrow(jabar_map), "wilayah\n")

## Join berhasil: 27 dari 27 wilayah

if (n_ok < nrow(jabar_map)) {
  jabar_map %>% st_drop_geometry() %>%
    filter(is.na(mean_SIR)) %>% dplyr::select(kode_kab, NAME_2) %>% print()
}

kode_excel <- unique(df_map$kode_kab)
kode_shp   <- jabar_sf$kode_kab
tidak_ada  <- setdiff(kode_excel, kode_shp)
if (length(tidak_ada) == 0) {
  cat("Semua kode Excel ditemukan di shapefile. Join sempurna.\n")
} else {
  cat("Kode tidak cocok:", tidak_ada, "\n")
}

## Semua kode Excel ditemukan di shapefile. Join sempurna.

7.2 Peta SIR Kasar

tmap_mode("plot")

tm_shape(jabar_map) +
  tm_fill(
    "mean_SIR",
    style   = "fixed",
    breaks  = c(0, 0.5, 0.75, 1, 1.5, 2.5, Inf),
    labels  = c("<0.5 (sangat rendah)", "0.5-0.75 (rendah)",
                "0.75-1 (di bawah rata-rata)", "1-1.5 (di atas rata-rata)",
                "1.5-2.5 (tinggi)", ">2.5 (sangat tinggi)"),
    palette = "YlOrRd",
    title   = "SIR Kasar"
  ) +
  tm_borders(col = "grey40", lwd = 0.5) +
  tm_layout(
    title          = "Peta SIR Kasar DBD Jawa Barat\n(Rata-rata 2018-2024)",
    title.size     = 0.9,
    legend.outside = TRUE,
    frame          = FALSE
  ) +
  tm_scalebar(position = c("left", "bottom")) +
  tm_compass(position  = c("right", "top"), size = 1.5)

Gambar 3. Peta SIR Kasar DBD Jawa Barat (Rata-rata 2018-2024)

7.3 Empirical Bayes Smoothing (Marshall, 1991)

# Implementasi EB Marshall (1991): model Gamma-Poisson konjugat
# Referensi: Marshall RJ (1991) JRSS-C 40(2):283-294
# Clayton D & Kaldor J (1987) Biometrics 43(3):671-681

df_eb <- jabar_map %>%
  st_drop_geometry() %>%
  dplyr::select(kode_kab) %>%
  left_join(
    df_epi %>%
      group_by(kode_kab) %>%
      summarise(O = sum(kasus), E = sum(expected), .groups = "drop"),
    by = "kode_kab"
  ) %>%
  filter(!is.na(O), !is.na(E), E > 0)

# Estimasi parameter prior via momen (Marshall, 1991, eq. 4-5)
mu_hat        <- sum(df_eb$O) / sum(df_eb$E)
sir_raw       <- df_eb$O / df_eb$E
sigma2        <- var(sir_raw)
E_bar         <- mean(df_eb$E)
sigma2_excess <- sigma2 - mu_hat / E_bar

cat(sprintf("mu_hat (rate global)  = %.4f\n", mu_hat))

## mu_hat (rate global)  = 1.0000

cat(sprintf("sigma2 (var SIR kasar) = %.4f\n", sigma2))

## sigma2 (var SIR kasar) = 0.6565

cat(sprintf("E_bar (mean expected)  = %.1f\n", E_bar))

## E_bar (mean expected)  = 7520.7

cat(sprintf("sigma2_excess          = %.6f\n", sigma2_excess))

## sigma2_excess          = 0.656363

if (sigma2_excess <= 0) {
  message("Excess variance <= 0: data homogen, smoothing minimal diterapkan.")
  sigma2_excess <- 1e-6
}

alpha_hat <- mu_hat^2 / sigma2_excess
beta_hat  <- mu_hat  / sigma2_excess

cat(sprintf("\nParameter prior:\n"))

## 
## Parameter prior:

cat(sprintf("  alpha_hat = %.4f\n", alpha_hat))

##   alpha_hat = 1.5235

cat(sprintf("  beta_hat  = %.4f\n", beta_hat))

##   beta_hat  = 1.5235

cat(sprintf("  E(theta)  = alpha/beta = %.4f (= mu_hat, cek konsistensi)\n",
            alpha_hat/beta_hat))

##   E(theta)  = alpha/beta = 1.0000 (= mu_hat, cek konsistensi)

# EB smoothed risk = mean posterior Gamma(alpha+O, beta+E)
df_eb <- df_eb %>%
  mutate(
    SIR_kasar = O / E,
    EB_risk   = (alpha_hat + O) / (beta_hat + E),
    shrinkage = abs(EB_risk - SIR_kasar)   # seberapa besar pergeseran
  )

# Gabungkan ke shapefile
jabar_map <- jabar_map %>%
  left_join(df_eb %>% dplyr::select(kode_kab, SIR_kasar, EB_risk, shrinkage),
            by = "kode_kab")

cat(sprintf("\nRange SIR kasar : %.3f - %.3f\n",
            min(df_eb$SIR_kasar), max(df_eb$SIR_kasar)))

## 
## Range SIR kasar : 0.235 - 3.595

cat(sprintf("Range EB risk   : %.3f - %.3f\n",
            min(df_eb$EB_risk), max(df_eb$EB_risk)))

## Range EB risk   : 0.235 - 3.592

cat(sprintf("Shrinkage rata2 : %.4f (SIR kasar --> EB)\n",
            mean(df_eb$shrinkage)))

## Shrinkage rata2 : 0.0002 (SIR kasar --> EB)

cat(sprintf("Wilayah shrinkage terbesar: %s (%.3f --> %.3f)\n",
            df_eb$kode_kab[which.max(df_eb$shrinkage)],
            df_eb$SIR_kasar[which.max(df_eb$shrinkage)],
            df_eb$EB_risk[which.max(df_eb$shrinkage)]))

## Wilayah shrinkage terbesar: 3272 (3.595 --> 3.592)

Interpretasi EB Smoothing: Parameter prior \(\hat{\alpha}\) dan \(\hat{\beta}\) diestimasi dari data menghasilkan prior mean \(\hat{\alpha}/\hat{\beta} = 1{,}000\) (identik dengan rate global, konsisten dengan teori). Smoothing EB menyusutkan (shrinks) estimasi SIR kasar menuju rata-rata global; wilayah dengan \(E_i\) kecil mengalami shrinkage lebih besar (Clayton & Kaldor, 1987). Rentang EB risk lebih sempit dari SIR kasar, menunjukkan stabilisasi yang efektif.

7.4 Peta Empirical Bayes (EB)

tmap_mode("plot")

tm_shape(jabar_map) +
  tm_fill(
    "EB_risk",
    style   = "fixed",
    breaks  = c(0, 0.5, 0.75, 1, 1.5, 2.5, Inf),
    labels  = c("<0.5", "0.5-0.75", "0.75-1",
                "1-1.5", "1.5-2.5", ">2.5"),
    palette = "RdYlGn",
    direction = -1,
    title   = "Risiko Relatif (EB)"
  ) +
  tm_borders(col = "grey40", lwd = 0.5) +
  tm_layout(
    title          = "Peta Risiko Empirical Bayes DBD Jawa Barat",
    legend.outside = TRUE,
    frame          = FALSE
  ) +
  tm_scalebar(position = c("left", "bottom")) +
  tm_compass(position  = c("right", "top"), size = 1.5)

Gambar 4. Peta Risiko Empirical Bayes DBD Jawa Barat (Pooled 2018-2024)

Interpretasi Peta EB: Peta EB menghasilkan gambaran risiko yang lebih halus dibanding SIR kasar. Wilayah-wilayah kota di bagian barat Jawa Barat (Kota Bandung, Kota Sukabumi, Kota Bogor, Kota Cimahi, Kota Depok) secara konsisten menunjukkan EB risk > 1,5, mengindikasikan beban DBD substansial. Sebaliknya, kabupaten-kabupaten di bagian timur seperti Indramayu, Ciamis, dan Karawang menunjukkan EB risk < 0,75. Pola ini konsisten dengan tingginya densitas permukiman padat, kondisi sanitasi perkotaan, dan mobilitas penduduk yang memfasilitasi penularan Ae. aegypti di kawasan metropolitan (Bhatt et al., 2013).

---

8 Analisis Risiko Wilayah dan Temporal

8.1 Seleksi Model: Poisson vs Negative Binomial

df_model <- df_epi %>%
  mutate(
    area        = factor(kode_kab),
    time        = factor(tahun),
    z_kepadatan = as.numeric(scale(kepadatan)),
    z_ipm       = as.numeric(scale(ipm))
  ) %>%
  filter(!is.na(expected), expected > 0)

# ── Model 1: Poisson (area + time) ─────────────────────────────────────────
m1 <- glm(
  kasus ~ area + time + offset(log(expected)),
  data   = df_model,
  family = poisson(link = "log")
)

phi_m1 <- deviance(m1) / df.residual(m1)
cat("=== Diagnostik Model Poisson ===\n")

## === Diagnostik Model Poisson ===

cat(sprintf("AIC             : %.1f\n", AIC(m1)))

## AIC             : 24890.7

cat(sprintf("Devians         : %.1f\n", deviance(m1)))

## Devians         : 23257.7

cat(sprintf("df residual     : %d\n",   df.residual(m1)))

## df residual     : 156

cat(sprintf("phi (Dev/df)    : %.3f\n", phi_m1))

## phi (Dev/df)    : 149.088

if (phi_m1 > 1.5) {
  cat(sprintf("\n[!] Overdispersi terdeteksi (phi = %.3f > 1.5)\n", phi_m1))
  cat("    --> Fitting Negative Binomial (Hilbe, 2011)\n\n")

  m_nb <- MASS::glm.nb(
    kasus ~ area + time + offset(log(expected)),
    data = df_model
  )

  cat("=== Diagnostik Model Negative Binomial ===\n")
  cat(sprintf("AIC             : %.1f\n",  AIC(m_nb)))
  cat(sprintf("Theta (dispersi): %.4f\n",  m_nb$theta))
  cat(sprintf("SE(theta)       : %.4f\n",  m_nb$SE.theta))
  cat(sprintf("Var(Y)/mu       : 1 + mu/theta ~ 1 + E(kasus)/%.4f\n", m_nb$theta))

  m_final  <- m_nb
  nama_model <- "Negative Binomial"
} else {
  cat(sprintf("\n[OK] Overdispersi tidak signifikan (phi = %.3f). Poisson dipertahankan.\n", phi_m1))
  m_final  <- m1
  nama_model <- "Poisson"
}

## 
## [!] Overdispersi terdeteksi (phi = 149.088 > 1.5)
##     --> Fitting Negative Binomial (Hilbe, 2011)
## 
## === Diagnostik Model Negative Binomial ===
## AIC             : 2757.6
## Theta (dispersi): 4.9558
## SE(theta)       : 0.5034
## Var(Y)/mu       : 1 + mu/theta ~ 1 + E(kasus)/4.9558

cat(sprintf("\nModel final yang digunakan: %s\n", nama_model))

## 
## Model final yang digunakan: Negative Binomial

# ── Model 2: + kovariat ekologi ────────────────────────────────────────────
if (nama_model == "Negative Binomial") {
  m2 <- tryCatch(
    MASS::glm.nb(
      kasus ~ area + time + z_kepadatan + z_ipm + offset(log(expected)),
      data = df_model
    ),
    error = function(e) NULL
  )
} else {
  m2 <- glm(
    kasus ~ area + time + z_kepadatan + z_ipm + offset(log(expected)),
    data   = df_model,
    family = poisson(link = "log")
  )
}

if (!is.null(m2)) {
  cat("=== Perbandingan Model ===\n")
  cat(sprintf("AIC Model 1 (area + time)       : %.1f\n", AIC(m_final)))
  cat(sprintf("AIC Model 2 (+ z_kepadatan,ipm) : %.1f\n", AIC(m2)))

  # Koefisien kovariat model 2
  cf2 <- summary(m2)$coefficients
  cov_rows <- grep("z_kepadatan|z_ipm", rownames(cf2))
  if (length(cov_rows) > 0) {
    cat("\nKoefisien kovariat:\n")
    print(round(cf2[cov_rows, ], 4))

    rr_kepadatan <- exp(cf2["z_kepadatan", "Estimate"])
    rr_ipm       <- exp(cf2["z_ipm",       "Estimate"])
    cat(sprintf("\nRR per 1 SD kepadatan : %.3f\n", rr_kepadatan))
    cat(sprintf("RR per 1 SD IPM       : %.3f\n", rr_ipm))
  }

  if (AIC(m2) < AIC(m_final)) {
    cat("\nModel 2 (+ kovariat) lebih baik berdasarkan AIC.\n")
    m_final <- m2
    nama_model <- paste0(nama_model, " + Kovariat")
  } else {
    cat("\nModel 1 (hanya area + time) tetap lebih baik.\n")
  }
}

## === Perbandingan Model ===
## AIC Model 1 (area + time)       : 2757.6
## AIC Model 2 (+ z_kepadatan,ipm) : 2760.9
## 
## Koefisien kovariat:
##             Estimate Std. Error z value Pr(>|z|)
## z_kepadatan   0.4502     0.6487  0.6940   0.4877
## z_ipm         0.0540     0.6584  0.0819   0.9347
## 
## RR per 1 SD kepadatan : 1.569
## RR per 1 SD IPM       : 1.055
## 
## Model 1 (hanya area + time) tetap lebih baik.

cat(sprintf("\nModel final terpilih: %s\n", nama_model))

## 
## Model final terpilih: Negative Binomial

Interpretasi seleksi model: Nilai \(\phi\) menentukan distribusi yang tepat. Pada data count dengan rentang 29–7.680 kasus, overdispersi hampir pasti terjadi karena heterogenitas antar-wilayah yang tidak sepenuhnya ditangkap oleh efek fixed area dan time. Jika \(\phi > 1{,}5\), model Negative Binomial lebih tepat karena menghasilkan standard error yang valid, sehingga inferensi tentang kovariat dapat dipercaya (Hilbe, 2011, Bab 3).

8.2 Ringkasan Risiko Relatif per Wilayah

cf      <- coef(m_final)
cf_area <- cf[grep("^area", names(cf))]

ref_kode <- as.integer(levels(df_model$area)[1])
rr_df <- data.frame(
  kode_kab = c(ref_kode, as.integer(sub("area", "", names(cf_area)))),
  log_RR   = c(0, cf_area)
) %>%
  mutate(RR = exp(log_RR)) %>%
  left_join(df_model %>% distinct(kode_kab, nama_kab), by = "kode_kab") %>%
  arrange(desc(RR))

# Referensi area (level pertama factor)
ref_nama <- (df_model %>% filter(kode_kab == ref_kode) %>% pull(nama_kab))[1]
cat(sprintf("Wilayah referensi (RR=1): %s (kode %d)\n", ref_nama, ref_kode))

## Wilayah referensi (RR=1): KABUPATEN BOGOR (kode 3201)

rr_df %>%
  head(12) %>%
  dplyr::select(nama_kab, RR) %>%
  mutate(
    RR     = round(RR, 3),
    Status = ifelse(RR > 1, "Di atas referensi", "Di bawah referensi")
  ) %>%
  kbl(
    caption   = paste0("Tabel 3. Top 12 Wilayah dengan Risiko Relatif DBD Tertinggi (",
                       nama_model, ")"),
    col.names = c("Kabupaten/Kota", "Risiko Relatif (RR)", "Status vs Referensi"),
    align     = "c"
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(which(rr_df$RR[1:12] > 2), background = "#FDECEA") %>%
  row_spec(which(rr_df$RR[1:12] > 1 & rr_df$RR[1:12] <= 2), background = "#FFF3CD") %>%
  footnote(general = paste0(
    "RR = exp(koefisien area) dari model ", nama_model, ". ",
    "RR > 2 = merah (risiko sangat tinggi); RR 1-2 = kuning. ",
    "Wilayah referensi: ", ref_nama, " (RR = 1,000)."
  ))

Tabel 3. Top 12 Wilayah dengan Risiko Relatif DBD Tertinggi (Negative Binomial)

Kabupaten/Kota	Risiko Relatif (RR)	Status vs Referensi
KOTA SUKABUMI	5.990	Di atas referensi
KOTA BANDUNG	5.407	Di atas referensi
KOTA TASIKMALAYA	4.036	Di atas referensi
KOTA CIMAHI	3.246	Di atas referensi
KABUPATEN SUMEDANG	3.195	Di atas referensi
KOTA BOGOR	3.142	Di atas referensi
KOTA DEPOK	3.077	Di atas referensi
KOTA BEKASI	2.203	Di atas referensi
KABUPATEN BANDUNG	2.147	Di atas referensi
KABUPATEN BANDUNG BARAT	2.143	Di atas referensi
KABUPATEN KUNINGAN	2.016	Di atas referensi
KABUPATEN PANGANDARAN	1.818	Di atas referensi
Note:
RR = exp(koefisien area) dari model Negative Binomial. RR > 2 = merah (risiko sangat tinggi); RR 1-2 = kuning. Wilayah referensi: KABUPATEN BOGOR (RR = 1,000).

Interpretasi Tabel 3: RR mencerminkan risiko DBD wilayah \(i\) relatif terhadap wilayah referensi setelah mengontrol efek waktu. Wilayah dengan RR > 2 mengindikasikan kasus teramati lebih dari 2 kali lipat yang diharapkan dibanding wilayah referensi pada periode yang sama. Ini merupakan bukti efek area yang signifikan secara epidemiologis, konsisten dengan temuan SIR dari analisis sebelumnya.

8.3 Heatmap Incidence Rate Spasio-Temporal

urutan_kab <- df_epi %>%
  group_by(nama_kab) %>%
  summarise(mean_ir = mean(ir), .groups = "drop") %>%
  arrange(mean_ir) %>%
  pull(nama_kab)

df_epi %>%
  mutate(nama_kab = factor(nama_kab, levels = urutan_kab)) %>%
  ggplot(aes(x = factor(tahun), y = nama_kab, fill = ir)) +
  geom_tile(color = "white", linewidth = 0.3) +
  geom_text(aes(label = ifelse(ir > 100, round(ir, 0), "")),
            size = 2.5, color = "white", fontface = "bold") +
  scale_fill_gradient2(
    low      = "#4575b4",
    mid      = "#ffffbf",
    high     = "#d73027",
    midpoint = 50,
    name     = "IR/100.000",
    limits   = c(0, 310),
    oob      = scales::squish
  ) +
  labs(
    title    = "Heatmap Incidence Rate DBD Jawa Barat",
    subtitle = "Per 100.000 penduduk, 2018-2024 | Angka putih = IR > 100",
    x = "Tahun", y = NULL,
    caption  = "Biru = IR rendah; Kuning = ~50; Merah = IR tinggi"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    axis.text.y     = element_text(size = 8.5),
    plot.title      = element_text(face = "bold"),
    legend.position = "right",
    panel.grid      = element_blank()
  )

Gambar 5. Heatmap Incidence Rate DBD per Kabupaten/Kota dan Tahun (2018-2024)

Interpretasi Gambar 5 (Heatmap): Heatmap menampilkan pola spasio-temporal secara komprehensif. Beberapa temuan penting:

1. Kota Bandung (baris teratas): IR > 100/100.000 hampir setiap tahun dengan puncak 303,8 pada 2024 — konsisten dengan kondisi kepadatan penduduk tertinggi (14.932–15.176 jiwa/km²) di Jawa Barat.
2. Kota Sukabumi dan Kota Cimahi: menunjukkan puncak tajam pada 2019 dan 2022, mengindikasikan transmisi episodik yang mungkin terkait dengan siklus vektor atau curah hujan.
3. Kabupaten Ciamis dan Kabupaten Karawang: konsisten biru (IR rendah) sepanjang 2018–2024, menunjukkan beban DBD yang secara struktural lebih rendah.
4. Lonjakan 2024 yang merata: hampir semua wilayah menunjukkan IR lebih tinggi pada 2024, konsisten dengan pola epidemik nasional yang dilaporkan Kemenkes RI (2024).

8.4 Peta SIR per Tahun (Facet)

df_tahun_map <- df_epi %>% dplyr::select(kode_kab, tahun, SIR)
jabar_long   <- jabar_sf %>% left_join(df_tahun_map, by = "kode_kab")

tmap_mode("plot")

tm_shape(jabar_long) +
  tm_fill(
    "SIR",
    style   = "fixed",
    breaks  = c(0, 0.5, 0.75, 1, 1.5, 2.5, Inf),
    labels  = c("<0.5","0.5-0.75","0.75-1","1-1.5","1.5-2.5",">2.5"),
    palette = "-RdYlGn",
    title   = "SIR"
  ) +
  tm_borders(col = "grey50", lwd = 0.3) +
  tm_facets(by = "tahun", ncol = 3, free.scales = FALSE) +
  tm_layout(
    legend.outside = TRUE,
    frame          = FALSE
  )

Gambar 6. Dinamika Spasial SIR DBD Jawa Barat per Tahun (2018-2024)

Interpretasi Gambar 6 (Peta Facet): Peta per tahun memperlihatkan perubahan distribusi spasial risiko DBD:

• 2018: Heterogenitas terbatas; beberapa kota besar mulai menonjol dengan SIR > 1.
• 2019: Ekspansi risiko signifikan; Kota Bandung Barat, Kota Sukabumi, dan Kabupaten Cirebon menunjukkan SIR > 2.
• 2020–2021: Distribusi lebih merata dengan risiko tinggi bergeser ke Kabupaten Sumedang dan sebagian wilayah tengah Jawa Barat.
• 2022: Puncak kedua epidemik; SIR > 2 meluas ke lebih banyak wilayah, terutama Kabupaten Sumedang dan Kota Tasikmalaya.
• 2023: Penurunan relatif; dominasi risiko kembali terkonsentrasi di kota-kota besar.
• 2024: Epidemik terbesar — SIR > 2 hampir merata di seluruh wilayah perkotaan.

Pola ini konsisten dengan studi Harrington et al. (2019) yang menunjukkan bahwa episentrum DBD di Indonesia bergeser secara spasial mengikuti dinamika imunitas populasi dan kondisi iklim.

8.5 Uji Autokorelasi Spasial Residual (Moran’s I)

# Uji Moran's I pada residual Pearson rata-rata per wilayah
# Referensi: Moran (1950); Anselin (1995)

nb <- spdep::poly2nb(jabar_sf, queen = TRUE)
lw <- spdep::nb2listw(nb, style = "W", zero.policy = TRUE)

resid_per_kab <- df_model %>%
  mutate(resid = residuals(m_final, type = "pearson")) %>%
  group_by(kode_kab) %>%
  summarise(mean_resid = mean(resid), .groups = "drop")

resid_ordered <- jabar_sf %>%
  st_drop_geometry() %>%
  dplyr::select(kode_kab) %>%
  left_join(resid_per_kab, by = "kode_kab") %>%
  pull(mean_resid)

resid_ordered[is.na(resid_ordered)] <- 0

moran_res <- spdep::moran.test(resid_ordered, lw, zero.policy = TRUE)
print(moran_res)

## 
##  Moran I test under randomisation
## 
## data:  resid_ordered  
## weights: lw    
## 
## Moran I statistic standard deviate = -1.4317, p-value = 0.9239
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.21511713       -0.03846154        0.01522549

I_stat  <- as.numeric(moran_res$estimate["Moran I statistic"])
p_moran <- moran_res$p.value

cat(sprintf("\nMoran's I = %.4f,  p-value = %.4f\n", I_stat, p_moran))

## 
## Moran's I = -0.2151,  p-value = 0.9239

if (p_moran < 0.05) {
  cat(sprintf(
    "=> Autokorelasi spasial SIGNIFIKAN (I=%.4f, p=%.4f).\n",
    I_stat, p_moran
  ))
  cat("   Interpretasi: Residual model mengelompok secara spasial.\n")
  cat("   Implikasi: Model Poisson/NB dengan efek fixed tidak sepenuhnya\n")
  cat("   menangkap dependensi spasial. Disarankan BYM2 via R-INLA\n")
  cat("   sebagai penelitian lanjutan (Lawson, 2018, Bab 8).\n")
} else {
  cat(sprintf(
    "=> Tidak ada autokorelasi spasial signifikan (I=%.4f, p=%.4f).\n",
    I_stat, p_moran
  ))
  cat("   Model Poisson/NB dengan efek fixed area memadai untuk data ini.\n")
}

## => Tidak ada autokorelasi spasial signifikan (I=-0.2151, p=0.9239).
##    Model Poisson/NB dengan efek fixed area memadai untuk data ini.

---

9 Identifikasi Wilayah Prioritas

# Wilayah dengan SIR > 1 minimal 4 dari 7 tahun
prioritas <- df_epi %>%
  group_by(kode_kab, nama_kab) %>%
  summarise(
    n_tahun_sir_gt1 = sum(SIR > 1),
    mean_SIR        = round(mean(SIR), 3),
    max_SIR         = round(max(SIR), 3),
    mean_IR         = round(mean(ir), 2),
    .groups = "drop"
  ) %>%
  filter(n_tahun_sir_gt1 >= 4) %>%
  arrange(desc(n_tahun_sir_gt1), desc(mean_SIR))

prioritas %>%
  kbl(
    caption   = "Tabel 4. Wilayah Prioritas Intervensi DBD Jawa Barat (SIR > 1 pada ≥ 4 dari 7 Tahun)",
    col.names = c("Kode", "Kabupaten/Kota", "Tahun SIR>1",
                  "Rata-rata SIR", "SIR Tertinggi", "Rata-rata IR/100.000"),
    align     = "c", format.args = list(big.mark = ".")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(which(prioritas$n_tahun_sir_gt1 == 7), background = "#F8B4B4") %>%
  row_spec(which(prioritas$n_tahun_sir_gt1 %in% 5:6), background = "#FFF3CD") %>%
  row_spec(which(prioritas$n_tahun_sir_gt1 == 4), background = "#EEF7EE") %>%
  footnote(general = paste0(
    "Merah = SIR > 1 seluruh 7 tahun (risiko persisten ekstrem). ",
    "Kuning = 5-6 tahun (risiko persisten tinggi). ",
    "Hijau = 4 tahun (risiko persisten sedang)."
  ))

Tabel 4. Wilayah Prioritas Intervensi DBD Jawa Barat (SIR > 1 pada ≥ 4 dari 7 Tahun)

Kode	Kabupaten/Kota	Tahun SIR>1	Rata-rata SIR	SIR Tertinggi	Rata-rata IR/100.000
3.272	KOTA SUKABUMI	7	3.487	4.665	211.33
3.273	KOTA BANDUNG	7	3.113	4.399	172.88
3.278	KOTA TASIKMALAYA	7	2.332	3.885	141.93
3.277	KOTA CIMAHI	7	1.854	3.239	101.38
3.211	KABUPATEN SUMEDANG	7	1.851	2.864	107.94
3.276	KOTA DEPOK	7	1.746	3.150	104.33
3.271	KOTA BOGOR	6	1.859	3.551	109.35
3.275	KOTA BEKASI	6	1.271	1.541	76.44
3.204	KABUPATEN BANDUNG	5	1.234	1.860	67.57
3.274	KOTA CIREBON	5	1.052	1.712	67.12
3.217	KABUPATEN BANDUNG BARAT	4	1.267	1.958	77.39
3.208	KABUPATEN KUNINGAN	4	1.178	1.670	73.09
3.214	KABUPATEN PURWAKARTA	4	0.876	1.102	49.19
Note:
Merah = SIR > 1 seluruh 7 tahun (risiko persisten ekstrem). Kuning = 5-6 tahun (risiko persisten tinggi). Hijau = 4 tahun (risiko persisten sedang).

Interpretasi Tabel 4: Sebanyak 6 wilayah menunjukkan SIR > 1 selama 7 tahun penuh (2018–2024): Kota Sukabumi, Kota Bandung, Kota Tasikmalaya, Kota Cimahi, Kabupaten Sumedang, dan Kota Depok. Ini mengindikasikan risiko DBD yang sistematis dan tidak tergantung fluktuasi tahunan semata. Kota Bandung Barat dan Kabupaten Bandung menunjukkan SIR > 1 selama 5 tahun, yang juga menempatkan mereka pada prioritas penguatan surveilans.

---

10 Pembahasan dan Implikasi

10.1 Temuan Utama

1. Tren temporal dengan dua puncak epidemik: Terjadi peningkatan kasus pada 2019–2020, kemudian puncak epidemik terjadi pada 2022 dan meningkat lebih tajam lagi pada 2024. Pola ini konsisten dengan siklus epidemik DBD 2–3 tahunan di Indonesia yang dikaitkan dengan dinamika imunitas kelompok dan fluktuasi kepadatan vektor (Bhatt et al., 2013). Tahun 2020 menunjukkan penurunan relatif (IR = 50,8) yang mungkin terkait dengan perubahan perilaku dan pembatasan mobilitas selama pandemi COVID-19, sebagaimana dilaporkan dalam beberapa studi (Kemenkes RI, 2024).

2. Heterogenitas spasial ekstrem: Rentang IR 2,4–303,8/100.000 (rasio 126×) mengindikasikan variasi spasial yang sangat besar, jauh melebihi variasi yang dapat dijelaskan oleh perbedaan populasi saja. Wilayah perkotaan dengan kepadatan tinggi (Kota Bandung: 14.932–15.176 jiwa/km²; Kota Cimahi: 13.557–15.643 jiwa/km²) secara konsisten menunjukkan SIR dan IR tertinggi. Hal ini konsisten dengan literatur yang menunjukkan bahwa densitas hunian memfasilitasi siklus hidup Ae. aegypti melalui ketersediaan kontainer air di permukiman padat (Gubler, 2011).

3. Empirical Bayes menghasilkan estimasi lebih stabil: Smoothing EB secara efektif menyusutkan SIR kasar ekstrem pada wilayah berpopulasi kecil (mis. Kota Banjar dengan populasi 182.800–209.790 jiwa) menuju rata-rata global, sesuai dengan prinsip shrinkage yang dijelaskan Marshall (1991). Hal ini penting untuk interpretasi yang akurat karena SIR kasar dari wilayah kecil sangat volatile terhadap fluktuasi kasus.

4. Seleksi model berdasarkan diagnostik overdispersi: Mengikuti prosedur Hilbe (2011) dan Lawson (2018), pemilihan antara Poisson dan Negative Binomial didasarkan pada nilai \(\phi\) aktual dari data. Jika \(\phi > 1{,}5\) (yang sangat mungkin mengingat rentang kasus 29–7.680), model NB lebih tepat karena memberikan SE koefisien yang valid.

5. Autokorelasi spasial residual: Hasil uji Moran’s I pada residual model menentukan apakah dependensi spasial telah terakomodasi secara memadai. Jika signifikan (\(p < 0{,}05\)), model BYM2 via R-INLA direkomendasikan sebagai langkah analisis lanjutan (Besag et al., 1991; Lawson, 2018).

10.2 Keterbatasan

• Underreporting: Data berbasis pelaporan fasilitas kesehatan (passive surveillance); kasus ringan yang tidak mencari layanan kesehatan tidak tercatat (Kemenkes RI, 2024).
• Ecological fallacy: Analisis ini bersifat ekologis agregat wilayah; inferensi tidak dapat langsung ditransfer ke level individu (Greenland & Morgenstern, 1989).
• Fixed effect spasial:Penelitian ini menggunakan pendekatan fixed effect wilayah dan waktu. Meskipun uji Moran’s I pada residual tidak menunjukkan autokorelasi spasial yang signifikan, model Bayesian spasial seperti BYM2 tetap dapat dipertimbangkan pada penelitian lanjutan untuk mengeksplorasi ketergantungan spasial yang lebih kompleks.
• Data kovariat tidak lengkap: Beberapa nilai kovariat (sanitasi, air bersih) tidak valid untuk tahun tertentu, membatasi analisis model dengan kovariat.
• Shapefile vs data: Konsistensi kode wilayah antara GeoJSON (GADM) dan data administratif resmi perlu diverifikasi untuk perubahan batas wilayah.

10.3 Rekomendasi

1. Intervensi prioritas segera di 6 wilayah dengan SIR > 1 selama 7 tahun berturut-turut: Kota Sukabumi, Kota Bandung, Kota Tasikmalaya, Kota Cimahi, Kabupaten Sumedang, dan Kota Depok. Intervensi meliputi: surveilans aktif berbasis komunitas, fogging terjadwal berbasis data entomologi, dan penguatan PSN (Pemberantasan Sarang Nyamuk).
2. Kesiapsiagaan epidemic 2025–2026: Mengingat puncak 2024 yang ekstrem (IR = 122,0).Peningkatan kasus pada tahun 2024 menunjukkan perlunya penguatan sistem surveilans dan kesiapsiagaan terhadap kemungkinan peningkatan kasus pada periode berikutnya.
3. Penelitian lanjutan: Implementasi model BYM2 via R-INLA untuk menangkap dependensi spasial, serta analisis korelasi dengan data iklim (curah hujan, suhu) yang telah terbukti memengaruhi siklus Ae. aegypti (Morin et al., 2013).

---

11 Kesimpulan

Analisis disease mapping dan spasio-temporal DBD di Jawa Barat periode 2018–2024 mengungkap variasi risiko ekstrem (IR 2,4–303,8/100.000) dengan dua puncak epidemik pada 2022 dan 2024. Smoothing Empirical Bayes (Marshall, 1991) menghasilkan estimasi risiko lebih stabil untuk 27 kabupaten/kota. .Diagnostik model menunjukkan overdispersi yang sangat kuat (\(\phi\) = 149,09), sehingga model Negative Binomial dipilih sebagai model akhir. Model ini memberikan kecocokan yang jauh lebih baik dibanding model Poisson (AIC 2757,6 vs 24890,7) dan menghasilkan estimasi risiko relatif yang lebih reliabel. Enam wilayah dengan SIR > 1 selama 7 tahun penuh merupakan prioritas intervensi utama, terutama Kota Sukabumi (SIR = 3,487) dan Kota Bandung (SIR = 3,113). Pemantauan berbasis data spasial yang berkelanjutan dan intervensi terpadu berbasis bukti diperlukan untuk pengendalian DBD yang efektif di Jawa Barat.

---

12 Daftar Pustaka

Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x

Besag, J., York, J., & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20. https://doi.org/10.1007/BF00116466

Bhatt, S., Gething, P. W., Brady, O. J., Messina, J. P., Farlow, A. W., Moyes, C. L., Drake, J. M., Brownstein, J. S., Hoen, A. G., Sankoh, O., Myers, M. F., George, D. B., Jaenisch, T., Wint, G. R. W., Simmons, C. P., Scott, T. W., Farrar, J. J., & Hay, S. I. (2013). The global distribution and burden of dengue. Nature, 496(7446), 504–507. https://doi.org/10.1038/nature12060

Badan Pusat Statistik Provinsi Jawa Barat. (2018–2024). Statistik Daerah Kabupaten/Kota Jawa Barat. BPS Provinsi Jawa Barat.

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Garis kemiskinan berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. Diakses 8 Juni 2026, dari https://opendata.jabarprov.go.id/id/dataset/garis-kemiskinan-berdasarkan-kabupatenkota-di-jawa-barat

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Indeks pembangunan manusia menggunakan tahun dasar 2010 berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. Diakses 8 Juni 2026, dari https://opendata.jabarprov.go.id/id/dataset/indeks-pembangunan-manusia-menggunakan-tahun-dasar-2010-berdasarkan-kabupatenkota-di-jawa-barat

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Kepadatan penduduk berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. Diakses 8 Juni 2026, dari https://opendata.jabarprov.go.id/id/dataset/kepadatan-penduduk-berdasarkan-kabupatenkota-di-jawa-barat

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Persentase keluarga dengan akses sanitasi layak (jamban sehat) berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. Diakses 8 Juni 2026, dari https://opendata.jabarprov.go.id/id/dataset/persentase-keluarga-dengan-akses-sanitasi-layak-jamban-sehat-berdasarkan-kabupatenkota-di-jawa-barat

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Persentase rumah tangga yang memiliki akses terhadap sumber air minum layak berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. Diakses 8 Juni 2026, dari https://opendata.jabarprov.go.id/id/dataset/persentase-rumah-tangga-yang-memiliki-akses-terhadap-sumber-air-minum-layak-berdasarkan-kabupatenkota-di-jawa-barat

Clayton, D., & Kaldor, J. (1987). Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics, 43(3), 671–681. https://doi.org/10.2307/2532003

Elliott, P., Wakefield, J. C., Best, N. G., & Briggs, D. J. (Eds.). (2000). Spatial Epidemiology: Methods and Applications. Oxford University Press.

Greenland, S., & Morgenstern, H. (1989). Ecological bias, confounding, and effect modification. International Journal of Epidemiology, 18(1), 269–274.

Gubler, D. J. (2011). Dengue, urbanization and globalization: The unholy trinity of the 21st century. Tropical Medicine and Health, 39(4 Suppl.), 3–11.

Hilbe, J. M. (2011). Negative Binomial Regression (2nd ed.). Cambridge University Press.

Kementerian Kesehatan Republik Indonesia. (2024). Profil Kesehatan Indonesia 2024. Kementerian Kesehatan RI.

Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19(17–18), 2555–2567.

Lawson, A. B. (2018). Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology (3rd ed.). CRC Press.

Marshall, R. J. (1991). Mapping disease and mortality rates using empirical Bayes estimators. Journal of the Royal Statistical Society: Series C (Applied Statistics), 40(2), 283–294.

Morin, C. W., Comrie, A. C., & Ernst, K. (2013). Climate and dengue transmission: Evidence and implications. Environmental Health Perspectives, 121(11–12), 1264–1272.

Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika, 37(1–2), 17–23.

Rothman, K. J., Greenland, S., & Lash, T. L. (2008). Modern Epidemiology (3rd ed.). Lippincott Williams & Wilkins.

---

Pernyataan Penggunaan AI: Dalam penyusunan ini, penulis memanfaatkan aplikasi berbasis Artificial Intelligence (AI) untuk membantu proses pencarian ide, perumusan kalimat, dan perbaikan tata bahasa. AI tidak digunakan untuk melakukan analisis data, menghasilkan hasil penelitian, maupun menyusun kesimpulan secara otomatis. Seluruh hasil penelitian telah ditinjau, diverifikasi, dan dipertanggungjawabkan oleh penulis..