Ringkasan dan Cara Menjalankan

Dokumen ini menerapkan model disease mapping Bayesian (Poisson log-linear) dengan efek spasial BYM2 menggunakan INLA pada data kasus DBD 27 kabupaten/kota di Jawa Barat. Alur dokumen mengikuti Bagian 7 modul:

Bagian modul	Isi pada dokumen ini
7.1 Model dasar disease mapping	Notasi \(y_i\), \(E_i\), \(\lambda_i\) (RR), SMR
7.2 Model spasial (IID, ICAR, BYM, BYM2)	Matriks ketetanggaan, matriks presisi ICAR, Moran’s I
7.3 Konsep dasar INLA	Latent Gaussian Model dan ringkasan pendekatan INLA
7.4 Prior dan PC priors	PC prior untuk \(\sigma\) dan \(\phi\)
7.5 Implementasi R-INLA	Model Poisson, IID, BYM, dan BYM2; posterior hiperparameter
7.6 Contoh dengan peta area nyata	Peta posterior RR, ketidakpastian, dan exceedance probability
7.7 Perbandingan model	IID vs BYM2 (DIC, WAIC, LCPO, perbandingan RR)

Syarat menjalankan

Data kasus DBD sudah diinput langsung di dalam dokumen (chunk data), sehingga file Data_Spasial.csv tidak diperlukan.
Koneksi internet diperlukan untuk (a) mengunduh peta batas kabupaten/kota dari BIG dan (b) memasang paket yang belum ada. Peta yang sudah terunduh disimpan sebagai jabar_2026.geojson dan dipakai ulang pada knit berikutnya. Jika unduhan gagal, letakkan sendiri file GeoJSON batas kabupaten/kota Jawa Barat dengan nama tersebut di folder yang sama.
Klik Knit (HTML). Seluruh keluaran tabel (CSV) tersimpan di folder Output_Bayesian_DBD_Jabar.

Persiapan Paket dan Pengaturan Awal

Paket CRAN yang belum terpasang akan dipasang otomatis. Paket INLA tidak ada di CRAN sehingga dipasang dari repositori khusus INLA (sesuai catatan 7.5.1 pada modul).

# Pastikan repositori CRAN terdefinisi saat knit non-interaktif
if (identical(unname(getOption("repos")["CRAN"]), "@CRAN@") ||
    is.null(getOption("repos")["CRAN"]) || is.na(getOption("repos")["CRAN"])) {
  options(repos = c(CRAN = "https://cloud.r-project.org"))
}

pkg_cran <- c("sf", "dplyr", "readr", "spdep", "ggplot2", "knitr", "tidyr")
pkg_hilang <- pkg_cran[!vapply(pkg_cran, requireNamespace, logical(1), quietly = TRUE)]
if (length(pkg_hilang) > 0) install.packages(pkg_hilang)

if (!requireNamespace("INLA", quietly = TRUE)) {
  install.packages(
    "INLA",
    repos = c(getOption("repos"),
              INLA = "https://inla.r-inla-download.org/R/stable"),
    dep = TRUE
  )
}

library(sf)
library(dplyr)
library(readr)
library(spdep)
library(ggplot2)
library(tidyr)
library(INLA)

set.seed(123)

output_dir <- "Output_Bayesian_DBD_Jabar"
if (!dir.exists(output_dir)) dir.create(output_dir)

7.1 Model Dasar Disease Mapping

Untuk setiap wilayah \(i = 1, \dots, n\) (di sini \(n = 27\)):

\(y_i\): jumlah kasus DBD yang diamati,
\(E_i\): expected count, yaitu jumlah kasus yang diharapkan jika risiko wilayah \(i\) sama dengan risiko rata-rata Jawa Barat,
\(\lambda_i\): risiko relatif (relative risk, RR) wilayah \(i\).

Model Poisson log-linear:

\[ y_i \mid \lambda_i \sim \text{Poisson}(E_i \lambda_i), \qquad \log(\lambda_i) = \eta_i = \alpha + u_i + v_i . \]

Karena data tidak memuat kovariat, bagian \(\mathbf{x}_i^\top \boldsymbol\beta\) pada modul tidak dipakai (model hanya memuat intercept \(\alpha\) dan efek acak area). Pada BYM2, \(u_i + v_i\) diganti oleh satu efek total \(b_i\) (lihat 7.2.4).

Data dan Perhitungan Expected Count

File CSV berisi proyeksi penduduk 2025 (ribu jiwa), angka kesakitan DBD per 100.000 penduduk, dan jumlah kasus DBD (ribu kasus). Nilai \(E_i\) dihitung dengan standardisasi internal:

\[ E_i = n_i \times \frac{\sum_i y_i}{\sum_i n_i}, \qquad \text{SMR}_i = \frac{y_i}{E_i}. \]

# Data diinput manual (sumber: Data_Spasial.csv), sehingga tidak perlu file eksternal.
# Kolom: wilayah, proyeksi penduduk 2025 (ribu jiwa), angka kesakitan DBD per 100.000
# penduduk, dan jumlah kasus DBD (ribu kasus).
data_raw <- data.frame(
  wilayah = c(
  "Bogor", "Sukabumi", "Cianjur", "Bandung", "Garut", "Tasikmalaya", "Ciamis",
  "Kuningan", "Cirebon", "Majalengka", "Sumedang", "Indramayu", "Subang",
  "Purwakarta", "Karawang", "Bekasi", "Bandung Barat", "Pangandaran", "Kota Bogor",
  "Kota Sukabumi", "Kota Bandung", "Kota Cirebon", "Kota Bekasi", "Kota Depok",
  "Kota Cimahi", "Kota Tasikmalaya", "Kota Banjar"
  ),
  pop_ribu = c(
  5721.62, 2852.11, 2610.32, 3783.21, 2748.7, 1933.81, 1266.33, 1225.49, 2413.81,
  1363.76, 1195.53, 1932.52, 1675.52, 1062.97, 2581.25, 3303.02, 1907.82, 436.49,
  1083.78, 370.68, 2548.78, 347.54, 2648.27, 2167.96, 606.38, 759.37, 211.96
  ),
  ir = c(
  39, 10, 48, 49, 90, 44, 58, 105, 48, 86, 121, 43, 75, 44, 294, 84, 97, 58, 72,
  274, 142, 374, 153, 142, 63, 92, 46
  ),
  kasus_ribu = c(
  2.2314318, 0.285211, 1.2529536, 1.8537729, 2.47383, 0.8508764, 0.7344714,
  1.2867645, 1.1586288, 1.1728336, 1.4465913, 0.8309836, 1.25664, 0.4677068,
  7.588875, 2.7745368, 1.8505854, 0.2531642, 0.7803216, 1.0156632, 3.6192676,
  1.2997996, 4.0518531, 3.0785032, 0.3820194, 0.6986204, 0.0975016
  ),
  stringsAsFactors = FALSE
)

stopifnot(ncol(data_raw) >= 4)
data_raw <- as.data.frame(data_raw[, 1:4])
names(data_raw) <- c("wilayah", "pop_ribu", "ir", "kasus_ribu")

# Fungsi penyeragaman nama wilayah agar penggabungan dengan peta tidak gagal
norm_nama <- function(x) {
  x <- toupper(trimws(as.character(x)))
  x <- gsub("^KABUPATEN\\s+|^KAB\\.?\\s+", "", x)
  x <- gsub("[^A-Z ]", " ", x)
  gsub("\\s+", " ", trimws(x))
}

dat <- data_raw %>%
  mutate(
    wilayah    = trimws(wilayah),
    pop_ribu   = as.numeric(pop_ribu),
    ir         = as.numeric(ir),
    kasus_ribu = as.numeric(kasus_ribu),
    id         = row_number(),
    kunci      = norm_nama(wilayah),
    pop        = pop_ribu * 1000,
    y          = round(kasus_ribu * 1000),
    E          = pop * (sum(y) / sum(pop)),
    SMR        = y / E
  )

stopifnot(nrow(dat) == 27, !any(is.na(dat$y)), !any(is.na(dat$pop)))

dat %>%
  select(id, wilayah, pop, ir, y, E, SMR) %>%
  mutate(across(c(E, SMR), ~ round(.x, 3))) %>%
  knitr::kable(caption = "Data DBD per kabupaten/kota (y = jumlah kasus, E = expected count)")

Data DBD per kabupaten/kota (y = jumlah kasus, E = expected count)
id	wilayah	pop	ir	y	E	SMR
1	Bogor	5721620	39	2231	5049.576	0.442
2	Sukabumi	2852110	10	285	2517.110	0.113
3	Cianjur	2610320	48	1253	2303.720	0.544
4	Bandung	3783210	49	1854	3338.845	0.555
5	Garut	2748700	90	2474	2425.846	1.020
6	Tasikmalaya	1933810	44	851	1706.670	0.499
7	Ciamis	1266330	58	734	1117.591	0.657
8	Kuningan	1225490	105	1287	1081.548	1.190
9	Cirebon	2413810	48	1159	2130.291	0.544
10	Majalengka	1363760	86	1173	1203.577	0.975
11	Sumedang	1195530	121	1447	1055.107	1.371
12	Indramayu	1932520	43	831	1705.532	0.487
13	Subang	1675520	75	1257	1478.718	0.850
14	Purwakarta	1062970	44	468	938.117	0.499
15	Karawang	2581250	294	7589	2278.064	3.331
16	Bekasi	3303020	84	2775	2915.057	0.952
17	Bandung Barat	1907820	97	1851	1683.733	1.099
18	Pangandaran	436490	58	253	385.221	0.657
19	Kota Bogor	1083780	72	780	956.482	0.815
20	Kota Sukabumi	370680	274	1016	327.141	3.106
21	Kota Bandung	2548780	142	3619	2249.408	1.609
22	Kota Cirebon	347540	374	1300	306.719	4.238
23	Kota Bekasi	2648270	153	4052	2337.212	1.734
24	Kota Depok	2167960	142	3079	1913.318	1.609
25	Kota Cimahi	606380	63	382	535.156	0.714
26	Kota Tasikmalaya	759370	92	699	670.177	1.043
27	Kota Banjar	211960	46	98	187.064	0.524

Peta Batas Wilayah (BIG)

url_big <- paste0(
  "https://geoservices.big.go.id/rbi/rest/services/",
  "BATASWILAYAH/BATAS_KABKOTA_AR/MapServer/0/query?",
  "where=WADMPR%3D%27Jawa%20Barat%27&outFields=*&maxAllowableOffset=0.001&f=geojson"
)
file_geo <- "jabar_2026.geojson"

if (!file.exists(file_geo)) {
  options(timeout = 300)
  ok <- tryCatch({
    download.file(url_big, destfile = file_geo, mode = "wb", quiet = TRUE)
    TRUE
  }, error = function(e) {
    message("Unduhan peta gagal: ", conditionMessage(e))
    FALSE
  }, warning = function(w) {
    message("Peringatan saat mengunduh peta: ", conditionMessage(w))
    FALSE
  })
  if (!ok && file.exists(file_geo)) file.remove(file_geo)
}

if (!file.exists(file_geo)) {
  stop("Peta tidak tersedia. Unduh GeoJSON batas kabupaten/kota Jawa Barat (BIG) dan simpan sebagai '",
       file_geo, "' di folder: ", getwd())
}

jabar_shp <- tryCatch(
  st_read(file_geo, quiet = TRUE),
  error = function(e) {
    file.remove(file_geo)
    stop("File GeoJSON tidak valid (sudah dihapus). Jalankan ulang agar diunduh kembali. Detail: ",
         conditionMessage(e))
  }
)
if (is.na(st_crs(jabar_shp))) st_crs(jabar_shp) <- 4326

# Cari kolom nama kabupaten/kota yang paling cocok dengan data (biasanya WADMKK)
kolom_chr <- names(jabar_shp)[vapply(st_drop_geometry(jabar_shp), is.character, logical(1))]
kolom_chr <- c(intersect("WADMKK", kolom_chr), setdiff(kolom_chr, "WADMKK"))
skor <- vapply(kolom_chr, function(k) sum(unique(norm_nama(jabar_shp[[k]])) %in% dat$kunci), numeric(1))
kolom_nama <- names(skor)[which.max(skor)]

if (max(skor) < nrow(dat)) {
  tdk_cocok <- setdiff(dat$kunci, norm_nama(jabar_shp[[kolom_nama]]))
  stop("Nama wilayah pada peta tidak cocok dengan data untuk: ", paste(tdk_cocok, collapse = ", "),
       ". Kolom nama peta yang dipakai: ", kolom_nama)
}
cat("Kolom nama pada peta yang dipakai:", kolom_nama, "\n")

## Kolom nama pada peta yang dipakai: WADMKK

sf_use_s2(FALSE)

jabar_shp$kunci <- norm_nama(jabar_shp[[kolom_nama]])
jabar_shp <- jabar_shp %>% filter(kunci %in% dat$kunci) %>% select(kunci)

# Jika satu wilayah terdiri atas beberapa fitur, gabungkan
if (anyDuplicated(jabar_shp$kunci) > 0) {
  jabar_shp <- jabar_shp %>% group_by(kunci) %>% summarise(.groups = "drop")
}

jabar_utm <- jabar_shp %>%
  inner_join(dat, by = "kunci") %>%
  arrange(id) %>%
  mutate(idx_area = row_number()) %>%
  st_transform(32748) %>%
  st_make_valid() %>%
  st_buffer(0)

stopifnot(nrow(jabar_utm) == 27)

# Versi lon/lat untuk pemetaan
jabar <- st_transform(jabar_utm, 4326)

ggplot(jabar) +
  geom_sf(fill = "grey90", color = "black", linewidth = 0.3) +
  labs(title = "Batas kabupaten/kota Jawa Barat (BIG)", x = "Longitude", y = "Latitude") +
  theme_bw()

7.2 Struktur Spasial: Ketetanggaan, ICAR, dan Eksplorasi

Matriks Ketetanggaan (Queen Contiguity)

Matriks \(W = (w_{ij})\) dengan \(w_{ij} = 1\) bila wilayah \(i\) dan \(j\) berbatasan. Karena peta dari BIG disederhanakan (maxAllowableOffset), batas bersama dua wilayah kadang tidak persis berimpit. Toleransi (snap) karena itu dinaikkan bertahap sampai seluruh wilayah terhubung menjadi satu komponen.

snap_kandidat <- c(1, 50, 100, 250, 500, 1000)
for (s in snap_kandidat) {
  nb <- poly2nb(jabar_utm, queen = TRUE, snap = s)
  ncomp <- n.comp.nb(nb)$nc
  if (ncomp == 1) break
}
cat("Toleransi snap yang dipakai (meter):", s, "\n")

## Toleransi snap yang dipakai (meter): 1

cat("Jumlah komponen graf:", ncomp, "\n")

## Jumlah komponen graf: 1

if (ncomp > 1) {
  warning("Graf ketetanggaan belum terhubung penuh. BYM2 tetap dapat dijalankan, tetapi periksa peta.")
}

W <- nb2mat(nb, style = "B", zero.policy = TRUE)
stopifnot(isTRUE(all.equal(W, t(W), check.attributes = FALSE)))

# Visualisasi jaringan ketetanggaan
pusat <- suppressWarnings(st_coordinates(st_centroid(st_geometry(jabar_utm))))
plot(st_geometry(jabar_utm), border = "grey50", main = "Jaringan ketetanggaan (queen)")
plot(nb, pusat, add = TRUE, col = "red", pch = 19, cex = 0.6)

# File graf untuk INLA (modul 7.5.3, ekstensi .graph)
file_adj <- file.path(output_dir, "Adjacency_Graph.graph")
nb2INLA(file_adj, nb)
g <- inla.read.graph(file_adj)
cat("Jumlah node pada graf INLA:", g$n, "\n")

## Jumlah node pada graf INLA: 27

stopifnot(g$n == nrow(jabar_utm))

Matriks Presisi ICAR

Prior ICAR (7.2.2) memakai matriks presisi \(Q = D - W\) dengan \(D = \text{diag}(m_i)\), \(m_i\) = jumlah tetangga wilayah \(i\). Matriks ini singular, sehingga diperlukan constraint \(\sum_i u_i = 0\) (pada INLA ditangani lewat constr = TRUE).

D <- diag(rowSums(W))
Q <- D - W
cat("Dimensi Q:", nrow(Q), "x", ncol(Q), "\n")

## Dimensi Q: 27 x 27

cat("Rank Q:", qr(Q)$rank, "(jika graf terhubung penuh, rank = n - 1 =", nrow(Q) - 1, ")\n")

## Rank Q: 26 (jika graf terhubung penuh, rank = n - 1 = 26 )

cat("Jumlah tetangga: min =", min(rowSums(W)), ", median =", median(rowSums(W)),
    ", maks =", max(rowSums(W)), "\n")

## Jumlah tetangga: min = 1 , median = 4 , maks = 8

Eksplorasi: Autokorelasi Spasial dan Overdispersi

lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
moran_smr     <- moran.mc(jabar_utm$SMR, listw = lw, nsim = 999, zero.policy = TRUE)
moran_log_smr <- moran.mc(log(jabar_utm$SMR), listw = lw, nsim = 999, zero.policy = TRUE)

model_glm    <- glm(y ~ 1 + offset(log(E)), family = poisson(), data = st_drop_geometry(jabar_utm))
pearson_disp <- sum(residuals(model_glm, type = "pearson")^2) / df.residual(model_glm)

print(moran_smr)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  jabar_utm$SMR 
## weights: lw  
## number of simulations + 1: 1000 
## 
## statistic = -0.17646, observed rank = 105, p-value = 0.895
## alternative hypothesis: greater

print(moran_log_smr)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  log(jabar_utm$SMR) 
## weights: lw  
## number of simulations + 1: 1000 
## 
## statistic = -0.26742, observed rank = 46, p-value = 0.954
## alternative hypothesis: greater

cat("Pearson dispersion model Poisson tanpa efek acak:", round(pearson_disp, 3), "\n")

## Pearson dispersion model Poisson tanpa efek acak: 1024.579

Pearson dispersion jauh di atas 1 menunjukkan overdispersi, sehingga model Poisson sederhana tidak memadai dan efek acak area (\(v_i\) atau \(b_i\)) diperlukan. Moran’s I pada SMR bernilai -0.176 (p-value Monte Carlo = 0.895) dan pada log-SMR bernilai -0.267 (p-value = 0.954).

Peta SMR

label_df <- jabar %>%
  st_point_on_surface() %>%
  st_coordinates() %>%
  as.data.frame() %>%
  setNames(c("X", "Y")) %>%
  mutate(id = jabar$id)

tema_peta <- theme_bw(base_family = "serif") +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(face = "bold", size = 12),
    legend.key = element_rect(color = "black", linewidth = 0.4)
  )

ggplot(jabar) +
  geom_sf(aes(fill = SMR), color = "black", linewidth = 0.35) +
  geom_text(data = label_df, aes(x = X, y = Y, label = id), size = 2.8, fontface = "bold") +
  scale_fill_gradient2(low = "#2ca25f", mid = "#ffea00", high = "#b30000",
                       midpoint = 1, name = "SMR") +
  labs(title = "Peta Standardized Morbidity Ratio (SMR) DBD Jawa Barat",
       x = "Longitude", y = "Latitude") +
  tema_peta

7.3 Konsep Dasar INLA (Ringkasan)

Model di atas adalah Latent Gaussian Model (LGM):

Likelihood: \(y_i \mid \eta_i \sim \text{Poisson}(E_i e^{\eta_i})\).
Latent field Gaussian: \(\mathbf{x} \mid \boldsymbol\theta \sim \text{N}(\mathbf{0}, Q(\boldsymbol\theta)^{-1})\), dengan \(\mathbf{x}\) memuat intercept dan efek acak spasial.
Hiperparameter: \(\boldsymbol\theta = (\tau, \phi)\) dengan prior \(\pi(\boldsymbol\theta)\).

INLA mengaproksimasi posterior hiperparameter \(\tilde\pi(\boldsymbol\theta \mid \mathbf{y})\) dengan Laplace approximation, mengintegrasikannya secara numerik pada ruang berdimensi rendah, dan menghitung posterior marginal \(\pi(x_k \mid \mathbf{y})\) secara langsung tanpa sampling MCMC.

7.4 Prior: PC Priors untuk BYM2

Model BYM2 (7.2.4) menulis efek total sebagai

\[ b_i = \frac{1}{\sqrt{\tau}}\left(\sqrt{1-\phi}\,\theta_i + \sqrt{\phi}\, s_i\right), \]

dengan \(\theta_i \sim \text{N}(0,1)\) (IID terstandar), \(s_i\) ICAR terstandar (varians marginal 1), \(\tau\) presisi global, dan \(\phi \in [0,1]\) proporsi variasi spasial terstruktur. Prior yang dipakai:

PC prior untuk \(\sigma = 1/\sqrt{\tau}\): \(\Pr(\sigma > U_\sigma) = \alpha_\sigma\), di sini \(\Pr(\sigma > 1) = 0{,}01\).
PC prior untuk \(\phi\): base model \(\phi = 0\) (tanpa struktur spasial), dengan \(\Pr(\phi < 0{,}5) = 0{,}5\).

U_SIGMA <- 1;   A_SIGMA <- 0.01   # Pr(sigma > 1) = 0.01
U_PHI   <- 0.5; A_PHI   <- 0.5    # Pr(phi < 0.5) = 0.5

7.5 Implementasi R-INLA

Sesuai contoh pada modul, ukuran eksposur dimasukkan lewat argumen E = E, dan indeks area idx_area harus berurutan sama dengan urutan node pada graf.

dat_inla <- as.data.frame(st_drop_geometry(jabar_utm))
dat_inla$idx_area <- seq_len(nrow(dat_inla))
stopifnot(max(dat_inla$idx_area) == g$n)

ctrl_compute <- list(dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE)

Model 1: Poisson tanpa efek acak (pembanding)

fit_pois <- inla(
  y ~ 1,
  family = "poisson", data = dat_inla, E = dat_inla$E,
  control.predictor = list(compute = TRUE),
  control.compute = ctrl_compute
)
summary(fit_pois)

## Time used:
##     Pre = 0.545, Running = 0.491, Post = 0.103, Total = 1.14 
## Fixed effects:
##             mean    sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept)    0 0.005     -0.009        0      0.009    0   0
## 
## Deviance Information Criterion (DIC) ...............: -19029.59
## Deviance Information Criterion (DIC, saturated) ....: -19154.63
## Effective number of parameters .....................: -20490.73
## 
## Watanabe-Akaike information criterion (WAIC) ...: 160490.84
## Effective number of parameters .................: 79481.77
## 
## Marginal log-Likelihood:  -10980.37 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Model 2: Poisson-IID (efek non-spasial, 7.2.1)

fit_iid <- inla(
  y ~ 1 + f(idx_area, model = "iid",
            hyper = list(prec = list(prior = "pc.prec", param = c(U_SIGMA, A_SIGMA)))),
  family = "poisson", data = dat_inla, E = dat_inla$E,
  control.predictor = list(compute = TRUE),
  control.compute = ctrl_compute
)
summary(fit_iid)

## Time used:
##     Pre = 0.267, Running = 0.989, Post = 0.0709, Total = 1.33 
## Fixed effects:
##              mean   sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) -0.12 0.14     -0.397    -0.12      0.157 -0.12   0
## 
## Random effects:
##   Name     Model
##     idx_area IID model
## 
## Model hyperparameters:
##                        mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for idx_area 1.99 0.516       1.12     1.95       3.13 1.85
## 
## Deviance Information Criterion (DIC) ...............: 293.59
## Deviance Information Criterion (DIC, saturated) ....: 53.92
## Effective number of parameters .....................: 26.92
## 
## Watanabe-Akaike information criterion (WAIC) ...: 285.45
## Effective number of parameters .................: 13.54
## 
## Marginal log-Likelihood:  -224.03 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Model 3: Poisson-BYM klasik (7.2.3 dan 7.5.4)

fit_bym <- tryCatch(
  inla(
    y ~ 1 + f(idx_area, model = "bym", graph = g, scale.model = TRUE,
              hyper = list(
                prec.unstruct = list(prior = "pc.prec", param = c(U_SIGMA, A_SIGMA)),
                prec.spatial  = list(prior = "pc.prec", param = c(U_SIGMA, A_SIGMA))
              )),
    family = "poisson", data = dat_inla, E = dat_inla$E,
    control.predictor = list(compute = TRUE),
    control.compute = ctrl_compute
  ),
  error = function(e) {
    message("Model BYM klasik dilewati: ", conditionMessage(e))
    NULL
  }
)
if (!is.null(fit_bym)) summary(fit_bym)

## Time used:
##     Pre = 0.299, Running = 0.854, Post = 0.0922, Total = 1.25 
## Fixed effects:
##              mean   sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) -0.12 0.13      -0.38    -0.12       0.14 -0.12   0
## 
## Random effects:
##   Name     Model
##     idx_area BYM model
## 
## Model hyperparameters:
##                                               mean       sd 0.025quant 0.5quant
## Precision for idx_area (iid component)        2.12 6.31e-01       1.10     2.04
## Precision for idx_area (spatial component) 1129.75 1.32e+04       2.42    66.85
##                                            0.975quant mode
## Precision for idx_area (iid component)           3.56 1.90
## Precision for idx_area (spatial component)    6868.37 4.36
## 
## Deviance Information Criterion (DIC) ...............: 293.60
## Deviance Information Criterion (DIC, saturated) ....: 53.93
## Effective number of parameters .....................: 26.93
## 
## Watanabe-Akaike information criterion (WAIC) ...: 285.45
## Effective number of parameters .................: 13.54
## 
## Marginal log-Likelihood:  -204.25 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Model 4: Poisson-BYM2 (model utama, 7.2.4 dan 7.5.5)

formula_bym2 <- y ~ 1 +
  f(idx_area, model = "bym2", graph = g, scale.model = TRUE, constr = TRUE,
    hyper = list(
      prec = list(prior = "pc.prec", param = c(U_SIGMA, A_SIGMA)),
      phi  = list(prior = "pc",      param = c(U_PHI, A_PHI))
    ))

fit_bym2 <- inla(
  formula_bym2,
  family = "poisson", data = dat_inla, E = dat_inla$E,
  control.predictor = list(compute = TRUE),
  control.compute = c(ctrl_compute, list(return.marginals.predictor = TRUE))
)
summary(fit_bym2)

## Time used:
##     Pre = 16.7, Running = 0.995, Post = 0.196, Total = 17.9 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) -0.12 0.125      -0.37    -0.12      0.129 -0.12   0
## 
## Random effects:
##   Name     Model
##     idx_area BYM2 model
## 
## Model hyperparameters:
##                         mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for idx_area 1.983 0.522      1.132    1.924       3.17 1.818
## Phi for idx_area       0.202 0.178      0.011    0.146       0.67 0.029
## 
## Deviance Information Criterion (DIC) ...............: 293.60
## Deviance Information Criterion (DIC, saturated) ....: 53.94
## Effective number of parameters .....................: 26.93
## 
## Watanabe-Akaike information criterion (WAIC) ...: 285.46
## Effective number of parameters .................: 13.54
## 
## Marginal log-Likelihood:  -205.43 
## CPO, PIT is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Posterior Hiperparameter: \(\sigma\) dan \(\phi\)

Presisi global \(\tau\) ditransformasi menjadi simpangan baku total \(\sigma = 1/\sqrt{\tau}\). Parameter \(\phi\) menunjukkan proporsi variasi spasial terstruktur (jawaban untuk pertanyaan kontribusi spatial structured effect): \(\phi \approx 0\) berarti variasi dominan non-spasial, \(\phi \approx 1\) berarti hampir seluruhnya spasial.

marg_sigma <- inla.tmarginal(function(x) 1 / sqrt(x),
                             fit_bym2$marginals.hyperpar[["Precision for idx_area"]])
marg_phi   <- fit_bym2$marginals.hyperpar[["Phi for idx_area"]]

ringkas_marg <- function(m) {
  c(Mean     = inla.emarginal(identity, m),
    Median   = inla.qmarginal(0.5, m),
    Lower_95 = inla.qmarginal(0.025, m),
    Upper_95 = inla.qmarginal(0.975, m))
}

hyper_summary <- data.frame(
  Parameter = c("Sigma (Total SD)", "Phi (Proporsi Spasial)"),
  rbind(ringkas_marg(marg_sigma), ringkas_marg(marg_phi)),
  row.names = NULL
)
hyper_summary[, -1] <- round(hyper_summary[, -1], 4)
knitr::kable(hyper_summary, caption = "Ringkasan posterior hiperparameter BYM2")

Ringkasan posterior hiperparameter BYM2
Parameter	Mean	Median	Lower_95	Upper_95
Sigma (Total SD)	0.7284	0.7204	0.5624	0.9373
Phi (Proporsi Spasial)	0.2025	0.1470	0.0114	0.6675

write_csv(hyper_summary, file.path(output_dir, "01b_Hyperparameter_BYM2.csv"))

df_hyper <- bind_rows(
  data.frame(Parameter = "Sigma (Total SD)",         x = marg_sigma[, 1], dens = marg_sigma[, 2]),
  data.frame(Parameter = "Phi (Proporsi Spasial)",   x = marg_phi[, 1],   dens = marg_phi[, 2])
)

ggplot(df_hyper, aes(x, dens)) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~ Parameter, scales = "free") +
  labs(x = NULL, y = "Posterior density",
       title = "Posterior hiperparameter BYM2") +
  theme_bw()

Rata-rata posterior \(\phi\) sebesar 0.202 (interval kredibel 95%: 0.011 sampai 0.667). Interval yang lebar menandakan data (27 area) hanya memberi informasi terbatas tentang \(\phi\).

7.7 Perbandingan Model

Kriteria: DIC, WAIC, dan LCPO (\(-\sum_i \log \text{CPO}_i\)); nilai yang lebih kecil lebih baik. Nilai MargLogLik adalah log marginal likelihood dari INLA (semakin besar semakin baik).

get_metrics <- function(m, nama) {
  cpo <- m$cpo$cpo
  data.frame(
    Model       = nama,
    DIC         = round(m$dic$dic, 2),
    WAIC        = round(m$waic$waic, 2),
    LCPO        = round(-sum(log(cpo), na.rm = TRUE), 2),
    CPO_gagal   = sum(m$cpo$failure > 0, na.rm = TRUE),
    MargLogLik  = round(m$mlik[1], 2)
  )
}

daftar_metrik <- list(
  get_metrics(fit_pois, "Poisson"),
  get_metrics(fit_iid,  "Poisson-IID")
)
if (!is.null(fit_bym)) daftar_metrik[[length(daftar_metrik) + 1]] <- get_metrics(fit_bym, "Poisson-BYM")
daftar_metrik[[length(daftar_metrik) + 1]] <- get_metrics(fit_bym2, "Poisson-BYM2")

model_comp <- bind_rows(daftar_metrik) %>% arrange(WAIC)
knitr::kable(model_comp, caption = "Perbandingan model (diurutkan berdasarkan WAIC)")

Perbandingan model (diurutkan berdasarkan WAIC)
Model	DIC	WAIC	LCPO	MargLogLik
Poisson-IID	293.59	285.45	174.39	-223.65
Poisson-BYM	293.59	285.45	174.28	-204.57
Poisson-BYM2	293.60	285.46	174.31	-205.78
Poisson	-19029.59	160490.84	6131.92	-10980.37

write_csv(model_comp, file.path(output_dir, "01_Model_Comparison.csv"))

Model terbaik berdasarkan WAIC: Poisson-IID.

Analisis Sensitivitas Prior BYM2

Pengaruh pilihan prior diuji dengan memvariasikan \(U_\sigma \in \{0{,}5;\,1;\,2\}\) dan \(\alpha_\phi \in \{1/3;\,1/2;\,2/3\}\) pada prior PC.

prior_grid <- expand.grid(U_SIGMA_s = c(0.5, 1, 2), A_PHI_s = c(1/3, 1/2, 2/3))

sens_table <- bind_rows(lapply(seq_len(nrow(prior_grid)), function(i) {
  u_s <- prior_grid$U_SIGMA_s[i]
  a_p <- prior_grid$A_PHI_s[i]
  f_s <- y ~ 1 + f(idx_area, model = "bym2", graph = g, scale.model = TRUE, constr = TRUE,
                   hyper = list(
                     prec = list(prior = "pc.prec", param = c(u_s, A_SIGMA)),
                     phi  = list(prior = "pc",      param = c(U_PHI, a_p))
                   ))
  m_s <- inla(f_s, family = "poisson", data = dat_inla, E = dat_inla$E,
              control.compute = list(dic = TRUE, waic = TRUE))
  data.frame(
    U_SIGMA  = u_s,
    A_PHI    = round(a_p, 3),
    DIC      = round(m_s$dic$dic, 2),
    WAIC     = round(m_s$waic$waic, 2),
    Phi_mean = round(inla.emarginal(identity, m_s$marginals.hyperpar[["Phi for idx_area"]]), 3)
  )
}))

knitr::kable(sens_table, caption = "Sensitivitas DIC, WAIC, dan rata-rata posterior phi terhadap prior")

Sensitivitas DIC, WAIC, dan rata-rata posterior phi terhadap prior
U_SIGMA	A_PHI	DIC	WAIC	Phi_mean
0.5	0.333	293.61	285.50	0.206
1.0	0.333	293.60	285.46	0.226
2.0	0.333	293.60	285.45	0.239
0.5	0.500	293.60	285.49	0.187
1.0	0.500	293.60	285.46	0.202
2.0	0.500	293.60	285.44	0.213
0.5	0.667	293.60	285.49	0.165
1.0	0.667	293.60	285.46	0.176
2.0	0.667	293.60	285.44	0.183

write_csv(sens_table, file.path(output_dir, "02_Sensitivity_Analysis.csv"))

Rentang WAIC antar skenario prior adalah 0.06 dan rentang rata-rata posterior \(\phi\) adalah 0.165 sampai 0.239.

7.6 Posterior Relative Risk dan Klasifikasi Risiko

Karena eksposur dimasukkan lewat E, prediktor linear dari INLA adalah \(\eta_i = \log \lambda_i\) sehingga \(\text{RR}_i = \exp(\eta_i)\). Selain rata-rata dan interval kredibel 95%, dihitung exceedance probability \(P(\text{RR}_i > 1)\) dan \(P(\text{RR}_i > 1{,}2)\).

n_area <- nrow(dat_inla)
marg_eta <- fit_bym2$marginals.linear.predictor
stopifnot(length(marg_eta) >= n_area)

post_rr <- bind_rows(lapply(seq_len(n_area), function(i) {
  m <- marg_eta[[i]]
  data.frame(
    id          = dat_inla$id[i],
    RR_mean     = inla.emarginal(exp, m),
    RR_lower    = exp(inla.qmarginal(0.025, m)),
    RR_upper    = exp(inla.qmarginal(0.975, m)),
    P_exceed    = 1 - inla.pmarginal(0, m),         # P(RR > 1)  <=> P(eta > 0)
    P_exceed_12 = 1 - inla.pmarginal(log(1.2), m)   # P(RR > 1.2)
  )
})) %>%
  mutate(RR_lebar = RR_upper - RR_lower)

map_result <- jabar %>%
  left_join(post_rr, by = "id") %>%
  mutate(
    kategori = case_when(
      P_exceed >= 0.95 ~ "Risiko Tinggi (P >= 0.95)",
      P_exceed >= 0.80 ~ "Cenderung Tinggi",
      P_exceed <= 0.05 ~ "Risiko Rendah (P <= 0.05)",
      P_exceed <= 0.20 ~ "Cenderung Rendah",
      TRUE             ~ "Tidak Berbeda Nyata"
    ),
    kategori = factor(kategori, levels = c(
      "Risiko Tinggi (P >= 0.95)", "Cenderung Tinggi", "Tidak Berbeda Nyata",
      "Cenderung Rendah", "Risiko Rendah (P <= 0.05)"
    ))
  )

hasil_tabel <- map_result %>%
  st_drop_geometry() %>%
  select(id, wilayah, y, E, SMR, RR_mean, RR_lower, RR_upper, P_exceed, P_exceed_12, kategori) %>%
  arrange(desc(RR_mean)) %>%
  mutate(across(c(E, SMR, RR_mean, RR_lower, RR_upper, P_exceed, P_exceed_12), ~ round(.x, 3)))

knitr::kable(hasil_tabel, caption = "Posterior relative risk BYM2 per kabupaten/kota (urut RR tertinggi)")

Posterior relative risk BYM2 per kabupaten/kota (urut RR tertinggi)
id	wilayah	y	E	SMR	RR_mean	RR_lower	RR_upper	P_exceed	P_exceed_12	kategori
22	Kota Cirebon	1300	306.719	4.238	4.230	4.005	4.465	1.000	1.000	Risiko Tinggi (P >= 0.95)
15	Karawang	7589	2278.064	3.331	3.330	3.256	3.406	1.000	1.000	Risiko Tinggi (P >= 0.95)
20	Kota Sukabumi	1016	327.141	3.106	3.098	2.912	3.293	1.000	1.000	Risiko Tinggi (P >= 0.95)
23	Kota Bekasi	4052	2337.212	1.734	1.733	1.680	1.787	1.000	1.000	Risiko Tinggi (P >= 0.95)
24	Kota Depok	3079	1913.318	1.609	1.609	1.553	1.666	1.000	1.000	Risiko Tinggi (P >= 0.95)
21	Kota Bandung	3619	2249.408	1.609	1.608	1.556	1.661	1.000	1.000	Risiko Tinggi (P >= 0.95)
11	Sumedang	1447	1055.107	1.371	1.370	1.301	1.442	1.000	1.000	Risiko Tinggi (P >= 0.95)
8	Kuningan	1287	1081.548	1.190	1.189	1.126	1.255	1.000	0.368	Risiko Tinggi (P >= 0.95)
17	Bandung Barat	1851	1683.733	1.099	1.099	1.050	1.150	1.000	0.000	Risiko Tinggi (P >= 0.95)
26	Kota Tasikmalaya	699	670.177	1.043	1.042	0.967	1.121	0.859	0.000	Cenderung Tinggi
5	Garut	2474	2425.846	1.020	1.020	0.980	1.060	0.831	0.000	Cenderung Tinggi
10	Majalengka	1173	1203.577	0.975	0.974	0.920	1.031	0.182	0.000	Cenderung Rendah
16	Bekasi	2775	2915.057	0.952	0.952	0.917	0.988	0.005	0.000	Risiko Rendah (P <= 0.05)
13	Subang	1257	1478.718	0.850	0.850	0.804	0.898	0.000	0.000	Risiko Rendah (P <= 0.05)
19	Kota Bogor	780	956.482	0.815	0.816	0.760	0.874	0.000	0.000	Risiko Rendah (P <= 0.05)
25	Kota Cimahi	382	535.156	0.714	0.715	0.646	0.789	0.000	0.000	Risiko Rendah (P <= 0.05)
18	Pangandaran	253	385.221	0.657	0.658	0.581	0.742	0.000	0.000	Risiko Rendah (P <= 0.05)
7	Ciamis	734	1117.591	0.657	0.657	0.611	0.706	0.000	0.000	Risiko Rendah (P <= 0.05)
4	Bandung	1854	3338.845	0.555	0.556	0.531	0.581	0.000	0.000	Risiko Rendah (P <= 0.05)
9	Cirebon	1159	2130.291	0.544	0.545	0.514	0.577	0.000	0.000	Risiko Rendah (P <= 0.05)
3	Cianjur	1253	2303.720	0.544	0.544	0.515	0.575	0.000	0.000	Risiko Rendah (P <= 0.05)
27	Kota Banjar	98	187.064	0.524	0.528	0.432	0.638	0.000	0.000	Risiko Rendah (P <= 0.05)
14	Purwakarta	468	938.117	0.499	0.500	0.457	0.547	0.000	0.000	Risiko Rendah (P <= 0.05)
6	Tasikmalaya	851	1706.670	0.499	0.499	0.467	0.534	0.000	0.000	Risiko Rendah (P <= 0.05)
12	Indramayu	831	1705.532	0.487	0.488	0.456	0.522	0.000	0.000	Risiko Rendah (P <= 0.05)
1	Bogor	2231	5049.576	0.442	0.442	0.424	0.461	0.000	0.000	Risiko Rendah (P <= 0.05)
2	Sukabumi	285	2517.110	0.113	0.115	0.102	0.129	0.000	0.000	Risiko Rendah (P <= 0.05)

write_csv(hasil_tabel, file.path(output_dir, "03_Posterior_Relative_Risk.csv"))

as.data.frame(table(Kategori = map_result$kategori)) %>%
  rename(Jumlah_Wilayah = Freq) %>%
  knitr::kable(caption = "Jumlah wilayah per kategori risiko")

Jumlah wilayah per kategori risiko
Kategori	Jumlah_Wilayah
Risiko Tinggi (P >= 0.95)	9
Cenderung Tinggi	2
Tidak Berbeda Nyata	0
Cenderung Rendah	1
Risiko Rendah (P <= 0.05)	15

Wilayah dengan risiko tinggi yang bertahan setelah smoothing (P(RR > 1) \(\ge\) 0,95): Kota Cirebon, Karawang, Kota Sukabumi, Kota Bekasi, Kota Depok, Kota Bandung, Sumedang, Kuningan, Bandung Barat.

Peta Posterior RR, Ketidakpastian, dan Exceedance Probability

ggplot(map_result) +
  geom_sf(aes(fill = RR_mean), color = "black", linewidth = 0.35) +
  geom_text(data = label_df, aes(x = X, y = Y, label = id), size = 2.8, fontface = "bold") +
  scale_fill_gradient2(low = "#2ca25f", mid = "#ffea00", high = "#b30000",
                       midpoint = 1, name = "Posterior RR") +
  labs(title = "Peta Posterior Mean Relative Risk (BYM2) DBD Jawa Barat",
       x = "Longitude", y = "Latitude") +
  tema_peta

ggplot(map_result) +
  geom_sf(aes(fill = RR_lebar), color = "black", linewidth = 0.35) +
  geom_text(data = label_df, aes(x = X, y = Y, label = id), size = 2.8, fontface = "bold") +
  scale_fill_gradient(low = "#ffffcc", high = "#800026", name = "Lebar CrI 95%\n(RR upper - lower)") +
  labs(title = "Peta Ketidakpastian Posterior RR (Lebar Interval Kredibel 95%)",
       x = "Longitude", y = "Latitude") +
  tema_peta

ggplot(map_result) +
  geom_sf(aes(fill = P_exceed), color = "black", linewidth = 0.35) +
  geom_text(data = label_df, aes(x = X, y = Y, label = id), size = 2.8, fontface = "bold") +
  scale_fill_gradient2(low = "#31a354", mid = "#ffea00", high = "#a50f15",
                       midpoint = 0.5, limits = c(0, 1), name = "P(RR > 1)") +
  labs(title = "Peta Exceedance Probability P(RR > 1) DBD Jawa Barat",
       x = "Longitude", y = "Latitude") +
  tema_peta

ggplot(map_result) +
  geom_sf(aes(fill = kategori), color = "black", linewidth = 0.35) +
  geom_text(data = label_df, aes(x = X, y = Y, label = id), size = 2.8, fontface = "bold") +
  scale_fill_manual(
    values = c("Risiko Tinggi (P >= 0.95)" = "#a50f15",
               "Cenderung Tinggi"          = "#fb6a4a",
               "Tidak Berbeda Nyata"       = "#f0f0f0",
               "Cenderung Rendah"          = "#a1d99b",
               "Risiko Rendah (P <= 0.05)" = "#238b45"),
    name = "Kategori risiko", drop = FALSE
  ) +
  labs(title = "Klasifikasi Risiko DBD Berdasarkan Exceedance Probability",
       x = "Longitude", y = "Latitude") +
  tema_peta

7.7 (Lanjutan) Perbandingan Risiko: IID-only vs BYM2

Seperti pada modul, RR dari model IID dibandingkan dengan RR dari BYM2. Titik di sekitar garis diagonal berarti kedua model memberi estimasi serupa; selisih menunjukkan efek smoothing spasial tambahan pada BYM2.

rr_banding <- data.frame(
  id      = dat_inla$id,
  wilayah = dat_inla$wilayah,
  SMR     = dat_inla$SMR,
  RR_iid  = exp(fit_iid$summary.linear.predictor$mean[seq_len(n_area)]),
  RR_bym2 = exp(fit_bym2$summary.linear.predictor$mean[seq_len(n_area)])
)

ggplot(rr_banding, aes(x = RR_iid, y = RR_bym2)) +
  geom_point(alpha = 0.7) +
  geom_text(aes(label = id), size = 3, vjust = -0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey30") +
  labs(x = "RR (IID-only)", y = "RR (BYM2)",
       title = "Perbandingan Relative Risk: IID-only vs BYM2") +
  theme_minimal()

rr_banding %>%
  mutate(across(c(SMR, RR_iid, RR_bym2), ~ round(.x, 3)),
         Selisih_BYM2_IID = round(RR_bym2 - RR_iid, 3)) %>%
  arrange(desc(abs(Selisih_BYM2_IID))) %>%
  knitr::kable(caption = "SMR, RR IID, dan RR BYM2 (diurutkan berdasarkan selisih terbesar)")

SMR, RR IID, dan RR BYM2 (diurutkan berdasarkan selisih terbesar)
id	wilayah	SMR	RR_iid	RR_bym2	Selisih_BYM2_IID
22	Kota Cirebon	4.238	4.226	4.229	0.003
27	Kota Banjar	0.524	0.527	0.525	-0.002
5	Garut	1.020	1.020	1.019	-0.001
18	Pangandaran	0.657	0.657	0.656	-0.001
20	Kota Sukabumi	3.106	3.096	3.097	0.001
1	Bogor	0.442	0.442	0.442	0.000
2	Sukabumi	0.113	0.115	0.115	0.000
3	Cianjur	0.544	0.544	0.544	0.000
4	Bandung	0.555	0.555	0.555	0.000
6	Tasikmalaya	0.499	0.499	0.499	0.000
7	Ciamis	0.657	0.657	0.657	0.000
8	Kuningan	1.190	1.189	1.189	0.000
9	Cirebon	0.544	0.544	0.544	0.000
10	Majalengka	0.975	0.974	0.974	0.000
11	Sumedang	1.371	1.370	1.370	0.000
12	Indramayu	0.487	0.488	0.488	0.000
13	Subang	0.850	0.850	0.850	0.000
14	Purwakarta	0.499	0.500	0.500	0.000
15	Karawang	3.331	3.330	3.330	0.000
16	Bekasi	0.952	0.952	0.952	0.000
17	Bandung Barat	1.099	1.099	1.099	0.000
19	Kota Bogor	0.815	0.815	0.815	0.000
21	Kota Bandung	1.609	1.608	1.608	0.000
23	Kota Bekasi	1.734	1.733	1.733	0.000
24	Kota Depok	1.609	1.608	1.608	0.000
25	Kota Cimahi	0.714	0.714	0.714	0.000
26	Kota Tasikmalaya	1.043	1.042	1.042	0.000

Kesimpulan

Overdispersi dan autokorelasi. Pearson dispersion model Poisson sederhana sebesar 1024.58 menunjukkan overdispersi yang kuat, sehingga efek acak area diperlukan. Moran’s I pada log-SMR bernilai -0.267 (p = 0.954).
Pemilihan model. Berdasarkan WAIC, model terbaik adalah Poisson-IID (WAIC = 285.45, DIC = 293.59).
Kontribusi efek spasial. Rata-rata posterior \(\phi\) sebesar 0.202 (CrI 95%: 0.011 sampai 0.667), dan simpangan baku total \(\sigma\) rata-rata 0.728.
Wilayah berisiko tinggi setelah smoothing. Wilayah dengan P(RR > 1) \(\ge\) 0,95 tercantum pada tabel dan peta klasifikasi di atas. Peta ketidakpastian membantu menilai wilayah mana yang estimasinya kurang pasti.
Sensitivitas prior. Hasil pada analisis sensitivitas menunjukkan seberapa stabil kesimpulan terhadap pilihan prior PC.

Catatan: Model Leroux (7.2.5) hanya dibahas konseptual pada modul dan tidak diestimasi di sini. Data tidak memuat kovariat, sehingga suku \(\mathbf{x}_i^\top\boldsymbol\beta\) tidak dimasukkan. Jika kovariat tersedia (misalnya curah hujan atau kepadatan penduduk), cukup tambahkan ke formula seperti contoh Cases ~ 1 + Rainfall + f(...) pada modul 7.6.

Informasi Sesi

sessionInfo()

## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Bangkok
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] INLA_26.08.07 Matrix_1.7-5  tidyr_1.3.2   ggplot2_4.0.3 spdep_1.4-2  
## [6] spData_2.3.5  readr_2.2.0   dplyr_1.2.1   sf_1.1-3     
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.62          bslib_0.12.0       lattice_0.22-9    
##  [5] tzdb_0.5.0         vctrs_0.7.3        tools_4.6.1        generics_0.1.4    
##  [9] parallel_4.6.1     tibble_3.3.1       proxy_0.4-29       pkgconfig_2.0.3   
## [13] KernSmooth_2.23-26 RColorBrewer_1.1-3 S7_0.2.2           lifecycle_1.0.5   
## [17] compiler_4.6.1     farver_2.1.2       deldir_2.0-4       fmesher_0.8.0     
## [21] textshaping_1.0.5  MatrixModels_0.5-4 htmltools_0.5.9    class_7.3-23      
## [25] sass_0.4.10        yaml_2.3.12        crayon_1.5.3       pillar_1.11.1     
## [29] jquerylib_0.1.4    classInt_0.4-11    cachem_1.1.0       wk_0.9.5          
## [33] boot_1.3-32        Deriv_4.3.6        tidyselect_1.2.1   digest_0.6.39     
## [37] purrr_1.2.2        labeling_0.4.3     splines_4.6.1      fastmap_1.2.0     
## [41] grid_4.6.1         cli_3.6.6          magrittr_2.0.5     e1071_1.7-17      
## [45] withr_3.0.3        scales_1.4.0       sp_2.2-4           bit64_4.8.6       
## [49] rmarkdown_2.32     bit_4.6.0          igraph_2.3.4       otel_0.2.0        
## [53] ragg_1.5.2         hms_1.1.4          evaluate_1.0.5     knitr_1.52        
## [57] s2_1.1.13          rlang_1.3.0        Rcpp_1.1.2         glue_1.8.1        
## [61] DBI_1.3.0          rstudioapi_0.19.0  vroom_1.7.1        jsonlite_2.0.0    
## [65] R6_2.6.1           systemfonts_1.3.2  units_1.0-1

Bayesian Spatial Disease Mapping DBD Jawa Barat dengan BYM2 dan INLA

Proyek Spatial Statistics — mengacu pada Bagian 7: Analisis Spasial dengan Pendekatan Bayesian - INLA

13 November 2026