1 Pendahuluan

1.1 Latar Belakang

Mengapa DBD di Jawa Barat perlu analisis spasial?
Selama 2018–2024, Jawa Barat mencatat 203.058 kasus DBD dengan distribusi yang sangat tidak merata secara geografis — Kota Sukabumi (IR rata-rata 211,3/100.000) vs Kabupaten Tasikmalaya (13,8/100.000): perbedaan 15× lipat.

Demam Berdarah Dengue (DBD) yang ditularkan melalui vektor nyamuk Aedes aegypti dan Aedes albopictus tetap menjadi salah satu beban kesehatan masyarakat paling kritis di Indonesia (Kemenkes RI, 2024). Sebagai provinsi berpenduduk terbanyak (>50 juta jiwa), Jawa Barat secara konsisten mencatat beban kasus tahunan yang masif.

Ketimpangan geografis tersebut mengonfirmasi adanya ketergantungan spasial (spatial dependence): estimasi risiko penularan suatu daerah bertaut erat dengan karakteristik lingkungan dan pola mobilitas wilayah tetangganya (Anselin, 1995). Pendekatan regresi Poisson standar mengasumsikan setiap unit wilayah bersifat independen, menghasilkan standar eror yang terlalu rendah dan parameter yang bias (Lawson, 2018).

1.2 Alasan Pendekatan Spasial BYM2

Penelitian ini menerapkan model BYM2 (Besag-York-Mollié 2; Riebler et al., 2016) karena:

  1. Menangani autokorelasi spasial — model secara eksplisit memodelkan struktur tetangga via ICAR
  2. Estimasi stabil untuk small area — meminjam kekuatan statistik dari wilayah bertetangga
  3. Parameter φ (phi) yang interpretatif — memisahkan klaster geografis dari heterogenitas acak dalam satu parameter tunggal
  4. Bayesian via INLA — distribusi posterior penuh, termasuk Exceedance Probability P(RR>1), dengan komputasi jauh lebih cepat dibanding MCMC (Rue et al., 2009)
  5. Integrasi temporal (RW1) — mendeteksi pergeseran pola risiko dari tahun ke tahun

1.3 Tujuan Analisis

Penelitian ini bertujuan untuk:

  1. Mengidentifikasi dan memetakan pola autokorelasi spasial distribusi kasus DBD di Jawa Barat 2018–2024
  2. Mengestimasi relative risk (RR) DBD per kabupaten/kota menggunakan model BYM2 Bayesian Spatiotemporal
  3. Mengidentifikasi wilayah hotspot (RR > 1) dan coldspot (RR < 1)
  4. Menganalisis pengaruh faktor sosial, ekonomi, dan kesehatan terhadap risiko DBD
  5. Menganalisis dinamika temporal risiko DBD selama 2018–2024

2 Data dan Metode

2.1 Sumber Data & Unit Observasi

Unit observasi: 27 kabupaten/kota Jawa Barat × 7 tahun (2018–2024) = 189 observasi panel

Sumber Data Link
Open Data Jabar Kasus DBD, Faskes, Sanitasi, Air Bersih opendata.jabarprov.go.id
BPS Jabar Kepadatan, IPM, Garis Kemiskinan, Penduduk opendata.jabarprov.go.id
GADM v4.1 Shapefile Indonesia Level 2 gadm.org

Catatan: Kabupaten Pangandaran (kode 3218) tidak tersedia dalam GADM 4.1. Geometri Kabupaten Ciamis digunakan sebagai proksi,ini merupakan keterbatasan analisis.

2.2 Variabel Penelitian

vars <- data.frame(
  No  = 1:8,
  Variabel = c("Kasus DBD (Y_it)", "IPM (%)", "Air Bersih (%)",
               "Sanitasi (%)", "Kepadatan Penduduk",
               "Garis Kemiskinan", "Faskes Rate", "Jumlah Penduduk"),
  Satuan = c("Orang/tahun", "Persen", "Persen", "Persen",
             "Jiwa/km² (log)", "Rp/kapita/bulan (log)",
             "Unit/10.000 pddk", "Jiwa"),
  Keterangan = c(
    "Variabel outcome; kasus terkonfirmasi per tahun per wilayah",
    "Indeks Pembangunan Manusia; prediktor kualitas hidup",
    "Proporsi RT dengan akses air bersih",
    "Proporsi RT dengan sanitasi layak; di-cap 100%",
    "Log-transformasi kepadatan penduduk per luas wilayah",
    "Log-transformasi garis kemiskinan; proksi daya beli",
    "Ketersediaan fasilitas kesehatan per kapita",
    "Dari satuan ribu orang ×1000; digunakan sebagai offset Poisson"
  )
)
kable(vars, caption = "Tabel 1. Variabel Penelitian", align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  column_spec(1, bold = TRUE, width = "3em") %>%
  column_spec(4, width = "35em") %>%
  footnote(general = "Semua prediktor distandarisasi (z-score) sebelum masuk model.",
           general_title = "Catatan:")
Tabel 1. Variabel Penelitian
No Variabel Satuan Keterangan
1 Kasus DBD (Y_it) Orang/tahun Variabel outcome; kasus terkonfirmasi per tahun per wilayah
2 IPM (%) Persen Indeks Pembangunan Manusia; prediktor kualitas hidup
3 Air Bersih (%) Persen Proporsi RT dengan akses air bersih
4 Sanitasi (%) Persen Proporsi RT dengan sanitasi layak; di-cap 100%
5 Kepadatan Penduduk Jiwa/km² (log) Log-transformasi kepadatan penduduk per luas wilayah
6 Garis Kemiskinan Rp/kapita/bulan (log) Log-transformasi garis kemiskinan; proksi daya beli
7 Faskes Rate Unit/10.000 pddk Ketersediaan fasilitas kesehatan per kapita
8 Jumlah Penduduk Jiwa Dari satuan ribu orang ×1000; digunakan sebagai offset Poisson
Catatan:
Semua prediktor distandarisasi (z-score) sebelum masuk model.

2.3 Spesifikasi Model BYM2

2.3.1 Likelihood Model

\[Y_{it} \sim \text{Poisson}(E_{it} \times RR_{it})\] \[\log(RR_{it}) = \alpha + \sum_{k=1}^{p}\beta_k X_{kit} + \underbrace{u_i + v_i}_{\text{BYM2 spasial}} + \underbrace{\gamma_t + \delta_t}_{\text{temporal}} + \phi_{it}\]

Di mana: - \(E_{it} = N_{it} \times \rho_{global}\) = expected count (offset Poisson)
- \(\rho_{global}\) = rate global Jawa Barat = 58,9593 per 100.000
- \(RR_{it} > 1\) → wilayah i tahun t berisiko di atas rata-rata Jawa Barat

2.3.2 Komponen Model

comp <- data.frame(
  Komponen   = c("u_i", "v_i", "γ_t", "δ_t", "φ_it", "β'X_it"),
  Nama       = c("Efek spasial terstruktur (ICAR)",
                 "Efek spasial tidak terstruktur (IID)",
                 "Efek temporal terstruktur (RW1)",
                 "Efek temporal tidak terstruktur (IID)",
                 "Interaksi ruang-waktu (Tipe I IID)",
                 "Fixed effects prediktor"),
  Interpretasi = c(
    "Klaster geografis: wilayah bertetangga cenderung berisiko mirip",
    "Heterogenitas wilayah yang tidak terkait lokasi",
    "Tren risiko yang berubah secara gradual antar tahun",
    "Fluktuasi acak per tahun; menangkap outbreak atau anomali",
    "Variasi residual yang tidak tertangkap komponen lain",
    "Pengaruh linear 6 prediktor (z-score)"
  )
)
kable(comp, caption = "Tabel 2. Komponen Model BYM2 Spatiotemporal") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = TRUE) %>%
  column_spec(1, bold = TRUE, color = "#1F3864", width = "5em")
Tabel 2. Komponen Model BYM2 Spatiotemporal
Komponen Nama Interpretasi
u_i Efek spasial terstruktur (ICAR) Klaster geografis: wilayah bertetangga cenderung berisiko mirip
v_i Efek spasial tidak terstruktur (IID) Heterogenitas wilayah yang tidak terkait lokasi
γ_t Efek temporal terstruktur (RW1) Tren risiko yang berubah secara gradual antar tahun
δ_t Efek temporal tidak terstruktur (IID) Fluktuasi acak per tahun; menangkap outbreak atau anomali
φ_it Interaksi ruang-waktu (Tipe I IID) Variasi residual yang tidak tertangkap komponen lain
β’X_it Fixed effects prediktor Pengaruh linear 6 prediktor (z-score)

2.3.3 Parameter BYM2: Phi (φ)

\[\theta_i = \frac{1}{\sqrt{\tau}} \times \left(\sqrt{\varphi}\, u_i^* + \sqrt{1-\varphi}\, v_i\right)\]

φ ∈ [0,1]: proporsi variabilitas yang disebabkan klaster spasial (ICAR).
- φ → 1: dominasi klaster geografis
- φ → 0: heterogenitas acak dominan

2.3.4 Prior PC (Penalized Complexity)

# Prior PC mengikuti Simpson et al. (2017)
# Mendorong parsimoni: penalti pada kompleksitas di luar null model

prior_bym2 <- list(
  prec = list(prior = "pc.prec", param = c(1.0, 0.01)),  # P(σ > 1.0) = 0.01
  phi  = list(prior = "pc",      param = c(0.5, 2/3))    # P(φ > 0.5) = 2/3
)
prior_rw1  <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))
prior_iid  <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))

3 Preprocessing Data

# ── 1. Load Data ────────────────────────────────────────────────────────────
data_raw <- read_excel("data_dbd_jabar__2018-2024.xlsx")

data_raw <- data_raw %>%
  rename(
    kode_kab     = kode_kabupaten_kota,
    nama_kab     = nama_kabupaten_kota,
    dbd          = `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(%)`,
    penduduk     = `jumlah penduduk (ribu orang)`
  )

# ── 2. Konversi & Cleaning ──────────────────────────────────────────────────
data_raw <- data_raw %>%
  mutate(
    air_bersih = as.numeric(air_bersih),
    sanitasi   = as.numeric(sanitasi),
    penduduk   = penduduk * 1000,          # KRITIS: ribu → jiwa
    sanitasi   = ifelse(!is.na(sanitasi) & sanitasi > 100, 100, sanitasi)  # Winsorization
  ) %>%
  group_by(kode_kab) %>%
  mutate(
    air_bersih = ifelse(is.na(air_bersih), mean(air_bersih, na.rm=TRUE), air_bersih),
    sanitasi   = ifelse(is.na(sanitasi),   mean(sanitasi,   na.rm=TRUE), sanitasi)
  ) %>%
  ungroup()

# ── 3. Feature Engineering ──────────────────────────────────────────────────
data_clean <- data_raw %>%
  mutate(
    ir_dbd           = (dbd / penduduk) * 100000,
    log_kepadatan    = log(kepadatan),
    log_garis_miskin = log(garis_miskin),
    faskes_rate      = (faskes / penduduk) * 10000,
    rate_global      = sum(dbd) / sum(penduduk),
    E_i              = penduduk * rate_global,
    SIR              = dbd / E_i,
    t_idx            = tahun - min(tahun) + 1
  )

# ── 4. Standardisasi (z-score) ──────────────────────────────────────────────
data_clean <- data_clean %>%
  mutate(
    z_ipm    = as.vector(scale(ipm)),
    z_air    = as.vector(scale(air_bersih)),
    z_san    = as.vector(scale(sanitasi)),
    z_kep    = as.vector(scale(log_kepadatan)),
    z_mis    = as.vector(scale(log_garis_miskin)),
    z_fas    = as.vector(scale(faskes_rate))
  )

cat("Rate global:", round(unique(data_clean$rate_global) * 1e5, 4), "per 100.000\n")
cat("Total obs:", nrow(data_clean), "(", n_distinct(data_clean$kode_kab), "× 7 tahun)\n")
## Verifikasi data:
##   • Sanitasi > 100%: Kota Sukabumi 2019 (152%) dan Kota Depok 2019 (173%) → di-cap ke 100%
##   • 1 NA air_bersih → imputasi mean wilayah
##   • 1 NA sanitasi → imputasi mean wilayah
##   • Rate global Jawa Barat: 58.9593 per 100.000 penduduk per tahun

4 Statistik Deskriptif

4.1 Ringkasan Per Tahun

tbl <- data.frame(
  Tahun = 2018:2024,
  `Total Kasus`    = c("12.492","25.282","24.471","23.454","36.608","19.328","61.423"),
  `Pddk Jabar`     = c("48.683.700","49.316.700","48.152.270","48.738.830",
                        "49.306.800","49.860.330","50.345.190"),
  `IR per 100rb`   = c("25,66","51,26","50,82","48,12","74,25","38,76","122,00"),
  `Wilayah Tertinggi` = c("Kota Bandung (112,9)","Kota Sukabumi (239,1)",
                           "Kota Tasikmalaya (197,4)","Kota Bandung (152,1)",
                           "Kota Sukabumi (289,2)","Kota Bogor (137,7)",
                           "Kota Sukabumi (445,9)"),
  `Morans I (p-perm)` = c("-0,049 (0,686) ns","-0,002 (0,482) ns",
                            "-0,020 (0,576) ns","-0,028 (0,605) ns",
                            "-0,143 (0,929) ns","-0,258 (0,995) ns",
                            "-0,228 (0,989) ns"),
  check.names = FALSE
)

kable(tbl,
      caption = "Tabel 3.1. Rekap Kasus DBD Jawa Barat 2018–2024",
      align   = c("c","r","r","r","l","c")) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed","bordered"),
                full_width = TRUE, font_size = 13) %>%
  row_spec(c(1,3,5,7), background = "#EBF3FB") %>%
  row_spec(7, bold = TRUE, background = "#FFE0E0") %>%
  footnote(
    general = "Penduduk dalam jiwa (satuan asli ribu × 1000). Moran's I: permutation test 9.999 replikasi, seed=2024.",
    general_title = "Catatan:"
  )
Tabel 3.1. Rekap Kasus DBD Jawa Barat 2018–2024
Tahun Total Kasus Pddk Jabar IR per 100rb Wilayah Tertinggi Morans I (p-perm)
2018 12.492 48.683.700 25,66 Kota Bandung (112,9) -0,049 (0,686) ns
2019 25.282 49.316.700 51,26 Kota Sukabumi (239,1) -0,002 (0,482) ns
2020 24.471 48.152.270 50,82 Kota Tasikmalaya (197,4) -0,020 (0,576) ns
2021 23.454 48.738.830 48,12 Kota Bandung (152,1) -0,028 (0,605) ns
2022 36.608 49.306.800 74,25 Kota Sukabumi (289,2) -0,143 (0,929) ns
2023 19.328 49.860.330 38,76 Kota Bogor (137,7) -0,258 (0,995) ns
2024 61.423 50.345.190 122,00 Kota Sukabumi (445,9) -0,228 (0,989) ns
Catatan:
Penduduk dalam jiwa (satuan asli ribu × 1000). Moran’s I: permutation test 9.999 replikasi, seed=2024.

Interpretasi: IR pooled global 58,96/100.000. Lonjakan terbesar 2024 (IR 122,00, meningkat 215% dari 2023). Semua nilai Moran’s I tidak signifikan — pola dispersi pada IR kasar; struktur spasial terdeteksi pada level residual model (φ = 0,644).

4.2 Plot Eksplorasi Deskriptif

knitr::include_graphics(.find_fig("fig1_deskriptif.png"))
Gambar 1. Eksplorasi deskriptif DBD Jawa Barat 2018–2024: (A) distribusi IR seluruh obs, (B) tren total kasus, (C) scatter IPM vs IR, (D) 10 wilayah IR tertinggi rata-rata.

Gambar 1. Eksplorasi deskriptif DBD Jawa Barat 2018–2024: (A) distribusi IR seluruh obs, (B) tren total kasus, (C) scatter IPM vs IR, (D) 10 wilayah IR tertinggi rata-rata.

Temuan utama dari plot:

  • Panel A (Distribusi IR): Right-skewed, mean ~70/100K; sebagian besar obs berIR rendah tetapi ada outlier ekstrem (Kota Sukabumi 2024 IR ~446/100K)
  • Panel B (Tren): Pola U-shape dengan nadir 2023, puncak 2024; penurunan 2020–2021 kemungkinan terkait pandemi COVID-19
  • Panel C (IPM vs IR): Korelasi positif — wilayah urban (IPM tinggi) cenderung berisiko lebih tinggi
  • Panel D (Top 10): Kota Sukabumi dominan (IR rata-rata 211,3/100K), diikuti Kota Bandung (172,9) dan Kota Tasikmalaya (141,9)

4.3 Matriks Korelasi

knitr::include_graphics(.find_fig("fig2_korelasi.png"))
Gambar 2. Matriks korelasi variabel penelitian. Korelasi tinggi antara IPM–Log Kepadatan (r=0,88) dan IPM–Log Garis Miskin (r=0,81) perlu diperhatikan.

Gambar 2. Matriks korelasi variabel penelitian. Korelasi tinggi antara IPM–Log Kepadatan (r=0,88) dan IPM–Log Garis Miskin (r=0,81) perlu diperhatikan.

catatan Multikolinearitas: IPM berkorelasi tinggi dengan Log Kepadatan (r=0,88) dan Log Garis Miskin (r=0,81). Hal ini menjelaskan mengapa Sanitasi tidak signifikan dalam model — kemungkinan variance inflation akibat kolinearitas dengan Air Bersih (r≈0,58). Model BYM2 dengan prior PC sebagian memitigasi masalah ini melalui regularisasi Bayesian.

5 Analisis Autokorelasi Spasial

5.1 Struktur Tetangga (Queen Contiguity)

knitr::include_graphics(.find_fig("fig_adjacency_structure.png"))
Gambar 3. Struktur tetangga queen contiguity — matriks bobot W untuk 27 kab/kota Jawa Barat. Total 238 pasangan tetangga (rata-rata 8,8 tetangga per wilayah).

Gambar 3. Struktur tetangga queen contiguity — matriks bobot W untuk 27 kab/kota Jawa Barat. Total 238 pasangan tetangga (rata-rata 8,8 tetangga per wilayah).

Matriks bobot W dibangun menggunakan kriteria queen contiguity (berbagi ≥1 titik batas), kemudian distandardisasi per baris (row-standardized) sehingga \(\sum_j w_{ij} = 1\).

# Membangun adjacency matrix dari shapefile
shp <- st_read("jabar.geojson", quiet = TRUE)
nb  <- poly2nb(shp, queen = TRUE)   # queen contiguity

cat("Jumlah wilayah:", length(nb), "\n")
cat("Rata-rata tetangga:", round(mean(card(nb)), 2), "\n")

# Simpan graph untuk INLA
tmp <- tempfile(fileext = ".graph")
nb2INLA(tmp, nb)

# Matriks bobot row-standardized untuk Moran's I
W_listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

5.2 Uji Global Moran’s I

5.2.1 Formula Moran’s I

\[I = \frac{n}{S_0} \cdot \frac{\sum_i \sum_j w_{ij} z_i z_j}{\sum_i z_i^2}\]

di mana \(z_i\) = deviasi IR wilayah \(i\) dari rata-rata global, \(S_0\) = jumlah semua bobot

# Uji Moran's I per tahun (permutation test 9.999 replikasi)
set.seed(2024)
moran_results <- lapply(2018:2024, function(yr) {
  dy <- filter(data_clean, tahun == yr) %>% arrange(kode_kab)
  mi <- moran.mc(dy$ir_dbd, W_listw, nsim = 9999,
                 alternative = "greater", zero.policy = TRUE)
  data.frame(tahun = yr, I = round(mi$statistic, 4),
             p_value = round(mi$p.value, 4),
             sig = ifelse(mi$p.value < 0.05, "***", "ns"))
})
moran_df <- do.call(rbind, moran_results)
moran_df <- data.frame(
  Tahun      = 2018:2024,
  `Moran I`  = c(-0.0490,-0.0016,-0.0204,-0.0276,-0.1427,-0.2580,-0.2280),
  `p-value`  = c(0.6857,0.4816,0.5762,0.6054,0.9287,0.9951,0.9889),
  Signifikan = rep("ns (tidak signifikan)", 7),
  check.names = FALSE
)
kable(moran_df, caption = "Tabel 3.2. Hasil Uji Global Moran's I (permutation test, 9.999 replikasi)", align="c") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  column_spec(4, color = "#6c757d", italic = TRUE)
Tabel 3.2. Hasil Uji Global Moran’s I (permutation test, 9.999 replikasi)
Tahun Moran I p-value Signifikan
2018 -0.0490 0.6857 ns (tidak signifikan)
2019 -0.0016 0.4816 ns (tidak signifikan)
2020 -0.0204 0.5762 ns (tidak signifikan)
2021 -0.0276 0.6054 ns (tidak signifikan)
2022 -0.1427 0.9287 ns (tidak signifikan)
2023 -0.2580 0.9951 ns (tidak signifikan)
2024 -0.2280 0.9889 ns (tidak signifikan)

Interpretasi: Moran’s I yang tidak signifikan menunjukkan tidak adanya autokorelasi spasial yang kuat pada IR DBD mentah. Sebaliknya, nilai φ = 0,644 menunjukkan bahwa setelah pengaruh kovariat dikendalikan, masih terdapat ketergantungan spasial yang cukup kuat pada risiko DBD.

5.3 Moran Scatterplot & LISA

knitr::include_graphics(.find_fig("fig4_moran_scatterplot.png"))
Gambar 4. Moran Scatterplot IR DBD 2022. Kemiringan regresi negatif (I = −0,143) mengkonfirmasi dispersi. Kota Sukabumi sebagai outlier ekstrem (HL quadrant).

Gambar 4. Moran Scatterplot IR DBD 2022. Kemiringan regresi negatif (I = −0,143) mengkonfirmasi dispersi. Kota Sukabumi sebagai outlier ekstrem (HL quadrant). 

knitr::include_graphics(.find_fig("fig5_peta_lisa.png"))
Gambar 5. Peta LISA DBD Jawa Barat 2022. Seluruh wilayah 'Tidak Signifikan' konsisten dengan Moran's I global yang tidak signifikan.

Gambar 5. Peta LISA DBD Jawa Barat 2022. Seluruh wilayah ‘Tidak Signifikan’ konsisten dengan Moran’s I global yang tidak signifikan.

5.4 Peta IR & SIR Observasi

knitr::include_graphics(.find_fig("fig3_peta_ir_2022.png"))
Gambar 6. Incidence Rate DBD Jawa Barat 2022. Klaster urban di kawasan Bandung Raya dan beberapa kota (Sukabumi, Kuningan) menonjol.

Gambar 6. Incidence Rate DBD Jawa Barat 2022. Klaster urban di kawasan Bandung Raya dan beberapa kota (Sukabumi, Kuningan) menonjol. 

knitr::include_graphics(.find_fig("fig3b_peta_sir_2022.png"))
Gambar 7. Standardized Incidence Ratio (SIR) DBD Jawa Barat 2022. Merah = kasus melebihi ekspektasi (SIR ≥ 1,25); biru = kasus di bawah ekspektasi berdasarkan rate global.

Gambar 7. Standardized Incidence Ratio (SIR) DBD Jawa Barat 2022. Merah = kasus melebihi ekspektasi (SIR ≥ 1,25); biru = kasus di bawah ekspektasi berdasarkan rate global.

6 Pemodelan BYM2 Spatiotemporal

6.1 Script Estimasi R-INLA

# ── Area index untuk INLA ───────────────────────────────────────────────────
kode_unik <- sort(unique(data_clean$kode_kab))
idx_map   <- data.frame(kode_kab = kode_unik, area_idx = seq_along(kode_unik))
data_clean <- data_clean %>% left_join(idx_map, by = "kode_kab")

data_clean <- data_clean %>%
  mutate(
    as = area_idx,          # spatial ICAR index
    ts = t_idx,             # temporal RW1 index
    ti = t_idx,             # temporal IID index
    si = area_idx           # spatial IID index (space-time interaction)
  )

# ── Prior PC ────────────────────────────────────────────────────────────────
prior_bym2 <- list(
  prec = list(prior = "pc.prec", param = c(1.0, 0.01)),
  phi  = list(prior = "pc",      param = c(0.5, 2/3))
)
prior_rw1 <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))
prior_iid <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))

# ── Formula BYM2 Spatiotemporal ──────────────────────────────────────────────
formula_bym2 <- outcome ~
  z_ipm + z_air + z_san + z_kep + z_mis + z_fas +
  f(as, model = "bym2",  graph = graph_file,
    scale.model = TRUE, constr = TRUE, hyper = prior_bym2) +
  f(ts, model = "rw1",   scale.model = TRUE,
    constr = TRUE, hyper = prior_rw1) +
  f(ti, model = "iid",   hyper = prior_iid) +
  f(si, model = "iid",   hyper = prior_iid)

# ── Estimasi INLA ────────────────────────────────────────────────────────────
mod <- inla(
  formula_bym2,
  family             = "poisson",
  data               = data_clean,
  E                  = E_i,
  control.compute    = list(dic = TRUE, waic = TRUE, cpo = TRUE),
  control.predictor  = list(compute = TRUE, link = 1),
  control.inla       = list(strategy = "adaptive"),
  verbose            = FALSE
)

# ── Ekstrak posterior ────────────────────────────────────────────────────────
data_clean$RR    <- mod$summary.fitted.values$mean
data_clean$RR_lo <- mod$summary.fitted.values$`0.025quant`
data_clean$RR_hi <- mod$summary.fitted.values$`0.975quant`

# Exceedance probability P(RR > 1)
data_clean$exc <- suppressWarnings(sapply(seq_len(nrow(data_clean)), function(i)
  tryCatch(1 - inla.pmarginal(1, mod$marginals.fitted.values[[i]]),
           error = function(e) 0.5)))

# Goodness-of-fit
cat("DIC  :", round(mod$dic$dic, 1), "\n")
cat("WAIC :", round(mod$waic$waic, 1), "\n")
cat("CPO valids:", sum(!is.na(mod$cpo$cpo)), "dari", nrow(data_clean), "\n")

6.2 Parameter φ (Phi BYM2)

6.2.1 Hasil Estimasi φ = 0,644 (95% CI: 0,521–0,762)

Interpretasi: 64,4% variabilitas risiko DBD bersifat spasial terstruktur — artinya, wilayah yang berdekatan secara geografis cenderung memiliki pola risiko yang mirip, bahkan setelah semua prediktor dikendalikan.

Faktor-faktor yang mendorong klaster geografis ini meliputi: - Kepadatan vektor Aedes aegypti yang bergantung pada kondisi lingkungan lokal - Pola mobilitas penduduk lintas wilayah (commuting) - Karakteristik infrastruktur perkotaan (drainase, permukiman padat)

Mengapa φ tinggi tapi Moran’s I tidak signifikan?

Ini bukan kontradiksi. Moran’s I mengukur autokorelasi IR mentah. BYM2 mengukur autokorelasi pada residual (setelah IPM, sanitasi, kepadatan, dan variabel lain dikontrol). Efek confounding spasial yang tersembunyi di balik prediktor inilah yang ditangkap φ = 0,644.

6.3 Kriteria Fit Model

fit_df <- data.frame(
  Kriteria = c("DIC", "WAIC", "CPO Valid"),
  Nilai    = c("-13.993,3", "56.768,8", "189 dari 189 (100%)"),
  Keterangan = c(
    "Nilai DIC yang lebih kecil menunjukkan model memiliki kecocokan yang lebih baik terhadap data",
    "Nilai WAIC menunjukkan performa model yang baik dengan tetap mempertimbangkan kompleksitas model",
    "Seluruh observasi dapat dijelaskan dengan baik oleh model dan tidak terdapat data yang terlalu menyimpang"
  )
)
kable(fit_df, caption = "Tabel 3.3. Kriteria Goodness-of-Fit Model BYM2") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = TRUE)
Tabel 3.3. Kriteria Goodness-of-Fit Model BYM2
Kriteria Nilai Keterangan
DIC -13.993,3 Nilai DIC yang lebih kecil menunjukkan model memiliki kecocokan yang lebih baik terhadap data
WAIC 56.768,8 Nilai WAIC menunjukkan performa model yang baik dengan tetap mempertimbangkan kompleksitas model
CPO Valid 189 dari 189 (100%) Seluruh observasi dapat dijelaskan dengan baik oleh model dan tidak terdapat data yang terlalu menyimpang

7 Hasil Model BYM2

7.1 Peta Risiko Posterior

knitr::include_graphics(.find_fig("fig6_peta_risiko_bym2.png"))
Gambar 8. Hasil model BYM2 Spatiotemporal: (A) IR observasi rata-rata, (B) Posterior Mean RR, (C) Exceedance Probability P(RR>1), (D) Efek spasial terstruktur ICAR (u_i). Panel (B) dan (C) adalah output utama model Bayesian.

Gambar 8. Hasil model BYM2 Spatiotemporal: (A) IR observasi rata-rata, (B) Posterior Mean RR, (C) Exceedance Probability P(RR>1), (D) Efek spasial terstruktur ICAR (u_i). Panel (B) dan (C) adalah output utama model Bayesian.

Interpretasi panel:

Panel Isi Interpretasi
A — IR Observasi Rata-rata IR 2018–2024 Data kasar; terpengaruh populasi kecil (Kota Sukabumi)
B — Posterior Mean RR Estimasi risiko setelah smoothing Bayesian Lebih stabil; “meminjam kekuatan” dari tetangga
C — P(RR>1) Exceedance probability Bukti Bayesian kekuatan hotspot
D — Efek Spasial (u_i) Komponen ICAR murni Mengisolasi efek “klaster geografis” dari efek lain

7.2 Fixed Effects (Rate Ratio)

fe <- data.frame(
  Prediktor  = c("IPM","Air Bersih","Sanitasi",
                 "Log Kepadatan","Log Garis Miskin","Faskes Rate"),
  RR         = c(1.203, 0.915, 0.959, 0.820, 1.119, 0.856),
  CI_lo      = c(1.132, 0.866, 0.895, 0.753, 1.038, 0.790),
  CI_hi      = c(1.279, 0.966, 1.027, 0.894, 1.205, 0.928),
  Signifikan = c("✓ Ya","✓ Ya","✗ Tidak","✓ Ya","✓ Ya","✓ Ya"),
  Arah       = c("Risiko (+)","Protektif (−)","Protektif (−)",
                 "Protektif (−)","Risiko (+)","Protektif (−)")
)

kable(fe, digits = 3,
      col.names = c("Prediktor","RR","CI 2,5%","CI 97,5%","Signifikan","Arah"),
      caption   = "Tabel 3.4. Estimasi Fixed Effects Model BYM2 (Rate Ratio = exp(β))") %>%
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = TRUE) %>%
  column_spec(2, bold = TRUE,
              color = ifelse(fe$RR > 1, "#CC0000", "#006600")) %>%
  column_spec(5, bold = TRUE,
              color = ifelse(fe$Signifikan == "✓ Ya", "#006600", "#6c757d")) %>%
  footnote(general = "(*) Signifikan jika 95% CI tidak melewati RR = 1. Semua prediktor dalam z-score.",
           general_title = "Catatan:")
Tabel 3.4. Estimasi Fixed Effects Model BYM2 (Rate Ratio = exp(β))
Prediktor RR CI 2,5% CI 97,5% Signifikan Arah
IPM 1.203 1.132 1.279 ✓ Ya Risiko (+)
Air Bersih 0.915 0.866 0.966 ✓ Ya Protektif (−)
Sanitasi 0.959 0.895 1.027 ✗ Tidak Protektif (−)
Log Kepadatan 0.820 0.753 0.894 ✓ Ya Protektif (−)
Log Garis Miskin 1.119 1.038 1.205 ✓ Ya Risiko (+)
Faskes Rate 0.856 0.790 0.928 ✓ Ya Protektif (−)
Catatan:
(*) Signifikan jika 95% CI tidak melewati RR = 1. Semua prediktor dalam z-score. 
knitr::include_graphics(.find_fig("fig8_forest_plot.png"))
Gambar 9. Forest plot Rate Ratio fixed effects BYM2. Titik = posterior mean; garis = 95% Credible Interval. Merah = faktor risiko (RR>1); biru = protektif (RR<1). Ukuran titik ∝ presisi estimasi.

Gambar 9. Forest plot Rate Ratio fixed effects BYM2. Titik = posterior mean; garis = 95% Credible Interval. Merah = faktor risiko (RR>1); biru = protektif (RR<1). Ukuran titik ∝ presisi estimasi.

Interpretasi kunci (per 1 SD kenaikan prediktor):

  • IPM (RR=1,203): Urbanisasi meningkatkan risiko 20,3% — wilayah dengan tingkat pembangunan yang lebih tinggi cenderung memiliki risiko DBD yang lebih besar, yang dapat dikaitkan dengan karakteristik wilayah perkotaan seperti kepadatan aktivitas dan mobilitas penduduk yang tinggi.
  • Air Bersih (RR=0,915): Akses air bersih menurunkan risiko DBD sebesar 8,5%. Hal ini menunjukkan bahwa akses air bersih yang lebih baik berperan sebagai faktor protektif terhadap kejadian DBD.
  • Log Kepadatan (RR=0,820):Kepadatan penduduk menurunkan risiko DBD sebesar 18,0%. Hasil ini mengindikasikan bahwa setelah pengaruh variabel lain dikendalikan, wilayah yang lebih padat tidak selalu memiliki risiko DBD yang lebih tinggi.
  • Log Garis Kemiskinan (RR=1,119): Log garis kemiskinan meningkatkan risiko DBD sebesar 11,9%. Hasil ini menunjukkan adanya hubungan positif antara tingkat kemiskinan dan risiko DBD.
  • Faskes Rate (RR=0,856): Jumlah fasilitas kesehatan menurunkan risiko DBD sebesar 14,4%. Hal ini menunjukkan bahwa ketersediaan fasilitas kesehatan yang lebih baik berkontribusi dalam menurunkan risiko DBD.
  • Sanitasi: Tidak signifikan (p>0,05) — kemungkinan collinear dengan Air Bersih

7.3 Profil Risiko Wilayah

knitr::include_graphics(.find_fig("fig9_caterpillar.png"))
Gambar 10. Caterpillar plot posterior RR per wilayah (rata-rata 2018–2024). Titik merah: RR>1 (di atas rata-rata Jawa Barat). Ukuran titik ∝ P(RR>1). Kota Bandung dan Kota Sukabumi: RR tertinggi dengan bukti Bayesian kuat.

Gambar 10. Caterpillar plot posterior RR per wilayah (rata-rata 2018–2024). Titik merah: RR>1 (di atas rata-rata Jawa Barat). Ukuran titik ∝ P(RR>1). Kota Bandung dan Kota Sukabumi: RR tertinggi dengan bukti Bayesian kuat. 

profil <- data.frame(
  Kategori   = c("Sangat Tinggi (RR > 1,2)", "Tinggi (1,0–1,2)", "Sedang (0,85–1,0)", "Rendah (< 0,85)"),
  Wilayah    = c("Kota Sukabumi, Kota Bandung, Kota Tasikmalaya",
                 "Kota Cimahi, Kab. Sumedang, Kota Depok, Kota Bogor",
                 "Kab. Bandung, Kab. Kuningan, Kota Bekasi",
                 "Kab. Tasikmalaya, Kab. Sukabumi, Kab. Indramayu, Kab. Bekasi"),
  `P(RR>1) > 0.8` = c("Ya (≥0,887)","Sebagian","Tidak","Tidak"),
  `Prioritas`     = c("🔴 Sangat Tinggi","🟠 Tinggi","🟡 Sedang","🟢 Aman"),
  check.names = FALSE
)
kable(profil, caption = "Tabel 3.5. Klasifikasi Wilayah Berdasarkan Posterior RR") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = TRUE) %>%
  row_spec(1, background = "#FFE0E0", bold = TRUE) %>%
  row_spec(4, background = "#D4EDDA")
Tabel 3.5. Klasifikasi Wilayah Berdasarkan Posterior RR
Kategori Wilayah P(RR>1) > 0.8 Prioritas
Sangat Tinggi (RR > 1,2) Kota Sukabumi, Kota Bandung, Kota Tasikmalaya Ya (≥0,887) 🔴 Sangat Tinggi |
Tinggi (1,0–1,2) Kota Cimahi, Kab. Sumedang, Kota Depok, Kota Bogor Sebagian 🟠 Tinggi |
Sedang (0,85–1,0) Kab. Bandung, Kab. Kuningan, Kota Bekasi Tidak 🟡 Sedang |
Rendah (< 0,85) Kab. Tasikmalaya, Kab. Sukabumi, Kab. Indramayu, Kab. Bekasi Tidak 🟢 Aman |

7.4 Efek Temporal RW1

knitr::include_graphics(.find_fig("fig7_efek_temporal.png"))
Gambar 11. Efek temporal terstruktur RW1 (γ_t). Garis biru = posterior mean; area = 95% Credible Interval. γ_t > 0 = tahun lebih berisiko dari rata-rata setelah kontrol semua faktor.

Gambar 11. Efek temporal terstruktur RW1 (γ_t). Garis biru = posterior mean; area = 95% Credible Interval. γ_t > 0 = tahun lebih berisiko dari rata-rata setelah kontrol semua faktor. 

temp_df <- data.frame(
  Tahun  = 2018:2024,
  `γ_t (posterior mean)` = c(0.42, 0.51, 0.27, 0.05, -0.14, -0.68, -0.46),
  Interpretasi = c(
    "Risiko intrinsik di atas rata-rata",
    "Puncak residual risiko (tertinggi)",
    "Di atas rata-rata (menurun)",
    "Mendekati rata-rata",
    "Di bawah rata-rata",
    "Terendah (konsisten pandemi)",
    "Mulai naik kembali"
  ),
  check.names = FALSE
)
kable(temp_df, caption = "Tabel 3.6. Efek Temporal Terstruktur RW1 (γ_t)") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(2, bold = TRUE, background = "#FFE0E0") %>%
  row_spec(6, background = "#D4EDDA")
Tabel 3.6. Efek Temporal Terstruktur RW1 (γ_t)
Tahun γ_t (posterior mean) Interpretasi
2018 0.42 Risiko intrinsik di atas rata-rata
2019 0.51 Puncak residual risiko (tertinggi)
2020 0.27 Di atas rata-rata (menurun)
2021 0.05 Mendekati rata-rata
2022 -0.14 Di bawah rata-rata
2023 -0.68 Terendah (konsisten pandemi)
2024 -0.46 Mulai naik kembali

Catatan: 2018–2019 berisiko tinggi secara intrinsik. 2020–2021: turun signifikan (COVID-19 restriction). 2022–2024: pemulihan dan eskalasi risiko — konsisten dengan siklus dengue multi-tahunan dan mobilitas pasca-pandemi.

7.5 Exceedance Probability

exc_df <- data.frame(
  Wilayah = c("Kota Sukabumi","Kota Bandung","Kota Tasikmalaya",
              "Kota Depok","Kab. Sumedang","Kota Cimahi"),
  `RR rata-rata` = c(1.55, 1.45, 1.22, 1.07, 1.10, 1.11),
  `P(RR>1)`      = c(0.991, 0.981, 0.887, 0.673, 0.727, 0.739),
  `IR rata-rata` = c("211,3","172,9","141,9","104,3","107,9","101,4"),
  `Prioritas`    = c("🔴 P1","🔴 P1","🔴 P1","🟠 P2","🟠 P2","🟠 P2"),
  check.names = FALSE
)
kable(exc_df, digits = 3,
      caption = "Tabel 3.7. Exceedance Probability Hotspot Prioritas") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(1:3, background = "#FFE0E0", bold = TRUE)
Tabel 3.7. Exceedance Probability Hotspot Prioritas
Wilayah RR rata-rata P(RR>1) IR rata-rata Prioritas
Kota Sukabumi 1.55 0.991 211,3 🔴 P1 |
Kota Bandung 1.45 0.981 172,9 🔴 P1 |
Kota Tasikmalaya 1.22 0.887 141,9 🔴 P1 |
Kota Depok 1.07 0.673 104,3 🟠 P2 |
Kab. Sumedang 1.10 0.727 107,9 🟠 P2 |
Kota Cimahi 1.11 0.739 101,4 🟠 P2 |

3 Hotspot Prioritas Tinggi (P(RR>1) > 0,88):

  1. Kota Sukabumi — P(RR>1) = 0,991 → bukti Bayesian sangat kuat
  2. Kota Bandung — P(RR>1) = 0,981 → populasi terbesar, beban absolut tertinggi
  3. Kota Tasikmalaya — P(RR>1) = 0,887 → klaster selatan yang perlu perhatian

8 Diagnostik Model

8.1 PIT Histogram

knitr::include_graphics(.find_fig("fig10_pit_histogram.png"))
Gambar 12. PIT (Probability Integral Transform) Histogram. Pola U-shape: spike di 0 dan 1 mengindikasikan model sedikit conservatif untuk nilai-nilai ekstrem. Secara keseluruhan, fit dapat diterima untuk data count epidemiologi.

Gambar 12. PIT (Probability Integral Transform) Histogram. Pola U-shape: spike di 0 dan 1 mengindikasikan model sedikit conservatif untuk nilai-nilai ekstrem. Secara keseluruhan, fit dapat diterima untuk data count epidemiologi.

Interpretasi PIT: Histogram PIT yang seragam (uniform) menunjukkan fit sempurna. Pola U-shape yang terlihat adalah umum untuk model penyakit dengan over-dispersion atau observasi ekstrem — tidak mengindikasikan kegagalan model secara substansial.

8.2 Observed vs Fitted

knitr::include_graphics(.find_fig("fig11_obs_vs_fitted.png"))
Gambar 13. Observed vs Fitted (skala log-log). Garis merah = garis sempurna (obs=fitted). Kepadatan titik di sekitar garis diagonal mengkonfirmasi model menangkap pola kasus dengan baik untuk 189 observasi.

Gambar 13. Observed vs Fitted (skala log-log). Garis merah = garis sempurna (obs=fitted). Kepadatan titik di sekitar garis diagonal mengkonfirmasi model menangkap pola kasus dengan baik untuk 189 observasi.

Kualitas fit: Titik-titik tersebar secara simetris di kedua sisi garis merah tanpa pola sistematik, menunjukkan tidak ada systematic bias. Beberapa titik yang menyimpang di pojok kiri bawah adalah observasi dengan kasus sangat kecil (rare events) — wajar untuk model Poisson.

8.3 Ringkasan Diagnostik

diag <- data.frame(
  Metrik    = c("DIC","WAIC","CPO valid","PIT shape","Obs vs Fitted"),
  Nilai     = c("-13.993,3","56.768,8","189/189 (100%)","U-shape ringan","R² ≈ 0,94 (log-log)"),
  Interpretasi = c(
    "Menunjukkan model memiliki kecocokan yang baik terhadap data",
    "Menunjukkan performa model yang baik dengan mempertimbangkan kompleksitas model",
    "Seluruh observasi dapat dijelaskan dengan baik oleh model",
    "Model sedikit kurang akurat pada nilai ekstrem, tetapi masih dapat diterima",
    "Model mampu menjelaskan sekitar 94% variasi kasus DBD"
  )
)
kable(diag, caption = "Tabel 3.8. Ringkasan Diagnostik Model BYM2") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = TRUE) %>%
  column_spec(3, color = "#006600")
Tabel 3.8. Ringkasan Diagnostik Model BYM2
Metrik Nilai Interpretasi
DIC -13.993,3 Menunjukkan model memiliki kecocokan yang baik terhadap data
WAIC 56.768,8 Menunjukkan performa model yang baik dengan mempertimbangkan kompleksitas model
CPO valid 189/189 (100%) Seluruh observasi dapat dijelaskan dengan baik oleh model
PIT shape U-shape ringan Model sedikit kurang akurat pada nilai ekstrem, tetapi masih dapat diterima
Obs vs Fitted R² ≈ 0,94 (log-log) Model mampu menjelaskan sekitar 94% variasi kasus DBD

9 Kesimpulan

5 Kesimpulan Utama Penelitian:

  1. Beban DBD masif dan tidak merata — 203.058 kasus 2018–2024 dengan disparitas 15× antara wilayah tertinggi dan terendah. Tren memuncak pada 2024 (IR 122,00/100.000; perlu verifikasi kelengkapan data).

  2. Struktur spasial terdeteksi pada level model — Moran’s I global tidak signifikan pada IR kasar, namun φ = 0,644 dari BYM2 mengkonfirmasi 64,4% variabilitas risiko residual bersifat spasial terstruktur. Ini memvalidasi keunggulan model spasial atas regresi Poisson biasa.

  3. Hotspot prioritas dengan bukti Bayesian kuat — Kota Sukabumi [P(RR>1)=0,991], Kota Bandung [0,981], Kota Tasikmalaya [0,887] harus mendapat prioritas tertinggi dalam program PSN, fogging, dan surveilans entomologi.

  4. Prediktor signifikan: efek urbanisasi dominan — IPM dan garis kemiskinan meningkatkan risiko (urbanisasi); air bersih, kepadatan, dan faskes bersifat protektif. Sanitasi tidak signifikan (kemungkinan kolinearitas dengan air bersih).

  5. Dinamika temporal mengkonfirmasi dampak pandemi — Penurunan γ_t pada 2020–2021 (COVID-19), pemulihan dan eskalasi 2022–2024. Pola ini konsisten dengan siklus dengue multi-tahunan dan peningkatan mobilitas pasca-pandemi.

9.1 Rekomendasi Kebijakan

Prioritas Wilayah Intervensi
🔴 Sangat Tinggi Kota Sukabumi, Kota Bandung, Kota Tasikmalaya Intensifikasi PSN, fogging terjadwal, surveilans entomologi
🟠 Tinggi Kota Cimahi, Kab. Sumedang, Kota Depok Pemantauan rutin, edukasi masyarakat
🟢 Fokus preventif Kab. Tasikmalaya, Kab. Sukabumi Pertahankan program air bersih & faskes

9.2 Keterbatasan

  1. Geometri Kab. Pangandaran menggunakan proksi Kab. Ciamis (GADM 4.1 belum update)
  2. Data 2024 kemungkinan belum tahun penuh (lonjakan 215% perlu verifikasi)
  3. Variabel iklim (curah hujan, suhu) tidak dimasukkan
  4. Mobilitas penduduk lintas wilayah tidak dimodelkan secara eksplisit

10 Daftar Pustaka

Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93–115.

Badan Pusat Statistik Provinsi Jawa Barat. (n.d.). Garis kemiskinan berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. 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. 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. 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. 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. https://opendata.jabarprov.go.id/id/dataset/persentase-rumah-tangga-yang-memiliki-akses-terhadap-sumber-air-minum-layak-berdasarkan-kabupatenkota-di-jawa-barat

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.

Bivand, R., Pebesma, E., & Gomez-Rubio, V. (2013). Applied spatial data analysis with R (2nd ed.). Springer.

Blangiardo, M., & Cameletti, M. (2015). Spatial and spatio-temporal Bayesian models with R-INLA. John Wiley & Sons.

Dinas Kesehatan Provinsi Jawa Barat. (n.d.). Jumlah kasus penyakit demam berdarah dengue (DBD) berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/jumlah-kasus-penyakit-demam-berdarah-dengue-dbd-berdasarkan-kabupatenkota-di-jawa-barat

Dinas Kesehatan Provinsi Jawa Barat. (n.d.). Jumlah fasilitas kesehatan berdasarkan jenis di Jawa Barat. Open Data Jabar. https://opendata.jabarprov.go.id/en/dataset/jumlah-fasilitas-kesehatan-berdasarkan-jenis-di-jawa-barat

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

Riebler, A., Sørbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145–1165. https://doi.org/10.1177/0962280216660421

Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2), 319–392. https://doi.org/10.1111/j.1467-9868.2008.00700.x

Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sørbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science, 32(1), 1–28. https://doi.org/10.1214/16-STS576

World Health Organization. (2023). Dengue and severe dengue. https://www.who.int/news-room/fact-sheets/detail/dengue-and-severe-dengue

11 Lampiran: Script R Lengkap

11.1 Formula BYM2

# Spesifikasi lengkap formula model BYM2 Spatiotemporal

formula_bym2 <- outcome ~
  # Fixed effects (prediktor z-score)
  z_ipm + z_air + z_san + z_kep + z_mis + z_fas +
  
  # BYM2: efek spasial (terstruktur ICAR + tak terstruktur IID)
  f(as, model     = "bym2",
    graph         = graph_file,
    scale.model   = TRUE,
    constr        = TRUE,
    hyper         = prior_bym2) +
  
  # RW1: efek temporal terstruktur (Random Walk orde 1)
  f(ts, model     = "rw1",
    scale.model   = TRUE,
    constr        = TRUE,
    hyper         = prior_rw1) +
  
  # IID temporal tak terstruktur
  f(ti, model     = "iid", hyper = prior_iid) +
  
  # Interaksi ruang-waktu Tipe I
  f(si, model     = "iid", hyper = prior_iid)

11.2 Prior PC

# Prior PC (Penalized Complexity) — Simpson et al. (2017)
# Filosofi: penalti kompleksitas → mendorong parsimoni → data lebih dominan

# BYM2: precision τ & mixing φ
prior_bym2 <- list(
  prec = list(prior = "pc.prec", param = c(1.0, 0.01)),  # P(σ > 1.0) = 0.01
  phi  = list(prior = "pc",      param = c(0.5, 2/3))    # P(φ > 0.5) = 2/3
)

# RW1 temporal terstruktur
prior_rw1 <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))

# IID temporal & interaksi ruang-waktu
prior_iid <- list(prec = list(prior = "pc.prec", param = c(1.0, 0.01)))

11.3 Ekstraksi Posterior

# Ekstrak komponen posterior model
n  <- nrow(data_clean)
na <- length(unique(data_clean$kode_kab))

# Fitted values (posterior RR)
data_clean$RR    <- mod$summary.fitted.values$mean
data_clean$RR_lo <- mod$summary.fitted.values$`0.025quant`
data_clean$RR_hi <- mod$summary.fitted.values$`0.975quant`

# Exceedance probability P(RR > 1)
data_clean$exc <- suppressWarnings(sapply(seq_len(n), function(i)
  tryCatch(1 - inla.pmarginal(1, mod$marginals.fitted.values[[i]]),
           error = function(e) 0.5)))

# Efek spasial terstruktur u_i (BYM2 structured component)
rf_spatial  <- mod$summary.random$as
spat_effect <- data.frame(
  kode_kab = sort(unique(data_clean$kode_kab)),
  spat_u   = rf_spatial$mean[(na+1):(2*na)]  # Indeks na+1:2na = komponen ICAR
)

# Efek temporal γ_t (RW1)
nt <- length(unique(data_clean$tahun))
rf_temporal <- mod$summary.random$ts
temp_effect <- data.frame(
  tahun   = sort(unique(data_clean$tahun)),
  gamma_t = rf_temporal$mean[1:nt],
  gam_lo  = rf_temporal$`0.025quant`[1:nt],
  gam_hi  = rf_temporal$`0.975quant`[1:nt]
)

# Phi (mixing parameter BYM2)
phi_val <- mod$summary.hyperpar["Phi for as", "mean"]
phi_ci  <- mod$summary.hyperpar["Phi for as", c("0.025quant","0.975quant")]
cat("phi =", round(phi_val, 3),
    "95% CI: [", round(phi_ci[1], 3), ",", round(phi_ci[2], 3), "]\n")

*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.