# ================================================================
# PEMODELAN BAYESIAN SPASIO-TEMPORAL RISIKO KUSTA
# KABUPATEN/KOTA PROVINSI JAWA BARAT TAHUN 2019-2024
# ================================================================

# ------------------------------------------------
# 1. FOLDER KERJA DAN PACKAGE
# ------------------------------------------------

setwd("D:/Semester 2/Epidemiologi")

library(readxl)
library(dplyr)
library(stringr)
library(ggplot2)
library(sf)
library(spdep)
library(tmap)
library(INLA)
library(openxlsx)
library(tibble)

# ------------------------------------------------
# 2. MEMBACA DATA
# ------------------------------------------------

kusta <- read_excel("Penyakit Kusta Jabar.xlsx")

peta_jabar <- st_read(
  "geojson_jawa_barat_adm2.geojson",
  quiet = TRUE
)

# ------------------------------------------------
# 3. FUNGSI MEMBERSIHKAN NAMA WILAYAH
# ------------------------------------------------

bersihkan_nama <- function(x) {
  x <- as.character(x)
  x <- str_to_upper(x)
  x <- str_trim(x)
  x <- str_replace_all(x, "\\s+", " ")
  x
}

# ------------------------------------------------
# 4. MENYIAPKAN DATA ANALISIS
# ------------------------------------------------

data_model <- kusta %>%
  mutate(
    nama_kabupaten_kota = bersihkan_nama(nama_kabupaten_kota),
    tahun = as.integer(tahun),
    jumlah_kasus = as.numeric(jumlah_kasus),
    
    # Jumlah penduduk pada data dalam satuan ribu jiwa
    penduduk = as.numeric(jumlah_penduduk) * 1000,
    
    jumlah_faskes = as.numeric(jumlah_faskes),
    jumlah_nakes = as.numeric(jumlah_nakes),
    sanitasi_layak = as.numeric(sanitasi_layak),
    air_minum_layak = as.numeric(air_minum_layak)
  ) %>%
  filter(tahun >= 2019, tahun <= 2024) %>%
  select(
    nama_kabupaten_kota,
    tahun,
    jumlah_kasus,
    penduduk,
    jumlah_faskes,
    jumlah_nakes,
    sanitasi_layak,
    air_minum_layak
  )

# Pastikan data lengkap: 27 wilayah x 6 tahun = 162 baris
stopifnot(nrow(data_model) == 162)

Abstrak. Kusta masih menjadi salah satu masalah kesehatan masyarakat karena dapat menyebabkan kecacatan apabila tidak ditemukan dan ditangani secara dini. Distribusi kasus kusta antarwilayah tidak selalu merata serta dapat berubah dari waktu ke waktu, sehingga diperlukan analisis yang mempertimbangkan aspek spasial dan temporal. Penelitian ini bertujuan untuk memetakan dan memodelkan risiko relatif kusta di kabupaten/kota Provinsi Jawa Barat tahun 2019–2024 menggunakan model Bayesian spasio-temporal dengan pendekatan Poisson BYM2. Data yang digunakan meliputi jumlah kasus kusta, jumlah penduduk, persentase sanitasi layak, persentase akses air minum layak, jumlah fasilitas kesehatan, dan jumlah tenaga kesehatan pada 27 kabupaten/kota di Jawa Barat. Jumlah penduduk digunakan untuk menghitung insidensi kusta dan membentuk expected cases sebagai offset dalam model. Model memasukkan efek spasial terstruktur dan tidak terstruktur melalui BYM2, efek temporal menggunakan random walk orde satu (RW1), serta interaksi spasio-temporal. Hasil analisis menunjukkan bahwa insidensi kusta di Jawa Barat selama tahun 2019–2024 berfluktuasi, dengan insidensi tertinggi pada tahun 2019 sebesar 4,26 per 100.000 penduduk dan menurun menjadi 2,39 per 100.000 penduduk pada tahun 2024. Pemetaan insidensi, Standardized Incidence Ratio (SIR), dan relative risk menunjukkan adanya ketidakmerataan risiko kusta antarwilayah. Model menghasilkan estimasi risiko relatif yang lebih stabil karena telah mempertimbangkan perbedaan populasi, keterkaitan antarwilayah, dan perubahan antarwaktu. Hasil penelitian diharapkan dapat menjadi dasar dalam menentukan wilayah prioritas untuk penguatan surveilans, penemuan kasus aktif, pelacakan kontak, serta intervensi pengendalian kusta di Jawa Barat.

Kata kunci: kusta, risiko relatif, Bayesian spasio-temporal, BYM2, Poisson, Jawa Barat.

1 PENDAHULUAN

1.1 Latar Belakang

Kusta (leprosy atau penyakit Hansen) merupakan penyakit infeksi kronis yang disebabkan oleh Mycobacterium leprae. Penyakit ini terutama menyerang kulit dan saraf perifer; apabila terlambat ditemukan serta diobati, kusta dapat menimbulkan gangguan saraf, kecacatan menetap, dan stigma yang memengaruhi perilaku pencarian pengobatan maupun kehidupan sosial penderita. Walaupun pengobatan multidrug therapy (MDT) tersedia dan efektif, deteksi dini serta pengobatan segera tetap diperlukan untuk mencegah kecacatan dan memutus rantai penularan (World Health Organization [WHO], 2026).

Urgensi pengendalian kusta masih tinggi pada tingkat global maupun nasional. Pada tahun 2024, Brazil, India, dan Indonesia masing-masing masih melaporkan lebih dari 10.000 kasus baru kusta. Indonesia juga masih melaporkan sekitar 14.000–15.000 kasus baru setiap tahun, yang menunjukkan bahwa penularan belum sepenuhnya terkendali dan pengendalian kusta memerlukan strategi yang lebih terarah menuju target zero leprosy tahun 2030 (WHO, 2026; Sebong et al., 2025).

Pada tingkat daerah, besarnya masalah kusta tidak cukup dilihat melalui jumlah kasus atau angka insidensi tingkat provinsi saja. Risiko kusta dapat berbeda antar kabupaten/kota karena perbedaan ukuran populasi, kondisi lingkungan, akses air minum dan sanitasi, serta ketersediaan fasilitas dan tenaga kesehatan. Selain itu, wilayah yang berdekatan secara geografis berpotensi memiliki kondisi sosial, lingkungan, mobilitas penduduk, dan akses layanan yang serupa. Konsekuensinya, risiko kusta pada suatu kabupaten/kota dapat berkaitan dengan risiko di wilayah tetangganya.

Data penelitian pada 27 kabupaten/kota di Jawa Barat selama 2019–2024 menunjukkan bahwa insidensi kusta berfluktuasi. Insidensi provinsi tercatat sebesar 4,26 kasus per 100.000 penduduk pada tahun 2019, menurun pada 2020–2021, meningkat kembali pada 2022–2023, dan menjadi 2,39 kasus per 100.000 penduduk pada tahun 2024. Pola tersebut menunjukkan perubahan beban kusta antarwaktu, tetapi angka agregat provinsi belum dapat menunjukkan kabupaten/kota yang mempunyai risiko lebih tinggi setelah perbedaan jumlah penduduk diperhitungkan. Oleh karena itu, diperlukan pemetaan risiko pada tingkat kabupaten/kota untuk mengidentifikasi variasi risiko yang tersembunyi di balik angka rata-rata provinsi.

Insidensi kasar berguna untuk menggambarkan kejadian penyakit, tetapi nilainya dapat tidak stabil pada wilayah dengan populasi kecil atau jumlah kasus sedikit. Dalam pemetaan penyakit, jumlah kasus lazim dimodelkan menggunakan distribusi Poisson dengan jumlah kasus harapan sebagai offset, sehingga perbedaan ukuran populasi antarwilayah dapat diperhitungkan dan hasilnya dinyatakan sebagai relative risk. Pendekatan ini membantu membedakan wilayah yang memiliki banyak kasus karena populasinya besar dari wilayah yang benar-benar memiliki risiko penyakit lebih tinggi.

Model Bayesian spasio-temporal dengan pendekatan Poisson BYM2 digunakan untuk mengakomodasi kebutuhan tersebut. Model BYM2 menggabungkan komponen spasial terstruktur—yang merepresentasikan kemiripan risiko antarwilayah bertetangga—dan komponen spasial tidak terstruktur untuk menangkap variasi spesifik wilayah. Parameterisasi BYM2 juga dirancang agar komponen spasial dapat diskalakan dan diinterpretasikan secara lebih konsisten dalam disease mapping (Riebler et al., 2016). Efek temporal dan interaksi spasio-temporal ditambahkan untuk menangkap perubahan risiko dari tahun ke tahun serta kemungkinan perubahan risiko yang berbeda pada setiap kabupaten/kota.

Berdasarkan uraian tersebut, penelitian ini bertujuan untuk memetakan persebaran insidensi kusta dan memodelkan risiko relatif kusta pada 27 kabupaten/kota di Jawa Barat selama tahun 2019–2024 menggunakan model Bayesian spasio-temporal Poisson BYM2. Model memasukkan jumlah kasus kusta sebagai respons, jumlah penduduk untuk membentuk kasus harapan sebagai offset, serta kovariat sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, dan rasio tenaga kesehatan. Hasil penelitian diharapkan menghasilkan informasi wilayah dengan risiko relatif tinggi sebagai dasar penguatan surveilans, penemuan kasus aktif, pemeriksaan kontak, dan penetapan prioritas pengendalian kusta di Jawa Barat.

1.2 Identifikasi Masalah

Berdasarkan latar belakang penelitian, identifikasi masalah dalam penelitian ini adalah sebagai berikut.

Kusta masih menjadi masalah kesehatan masyarakat karena keterlambatan penemuan dan pengobatan kasus dapat menyebabkan kecacatan permanen serta stigma pada penderita.
Insidensi kusta di Jawa Barat selama tahun 2019–2024 menunjukkan pola fluktuatif, sehingga perlu dianalisis perubahan risiko kusta dari waktu ke waktu.
Nilai insidensi tingkat provinsi belum mampu menggambarkan perbedaan risiko kusta pada setiap kabupaten/kota di Jawa Barat.
Jumlah kasus kusta antarwilayah tidak dapat dibandingkan secara langsung karena jumlah penduduk pada setiap kabupaten/kota berbeda.
Risiko kusta diduga memiliki keterkaitan spasial, yaitu kabupaten/kota yang bertetangga berpotensi mempunyai tingkat risiko yang serupa.
Kondisi sanitasi layak, akses air minum layak, fasilitas kesehatan, dan tenaga kesehatan diduga berkaitan dengan variasi risiko kusta antarwilayah.
Diperlukan metode pemodelan yang mampu memperhitungkan jumlah kasus, perbedaan populasi, ketergantungan spasial, perubahan waktu, serta faktor-faktor pendukung secara simultan.

1.3 Rumusan Masalah

Bagaimana pola insidensi kusta di kabupaten/kota Jawa Barat selama tahun 2019–2024?
Bagaimana persebaran spasial insidensi dan Standardized Incidence Ratio (SIR) kusta di Jawa Barat?
Bagaimana estimasi risiko relatif kusta pada setiap kabupaten/kota di Jawa Barat tahun 2019–2024 menggunakan model Bayesian spasio-temporal Poisson BYM2?
Bagaimana pengaruh sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, dan rasio tenaga kesehatan terhadap risiko relatif kusta di Jawa Barat?
Kabupaten/kota mana yang memiliki risiko relatif kusta tertinggi dan dapat menjadi wilayah prioritas pengendalian kusta?

1.4 Tujuan Penelitian

Mendeskripsikan tren insidensi kusta di Jawa Barat selama tahun 2019–2024.
Memetakan persebaran insidensi kusta dan SIR pada kabupaten/kota di Jawa Barat.
Mengestimasi risiko relatif kusta pada setiap kabupaten/kota di Jawa Barat menggunakan model Bayesian spasio-temporal Poisson BYM2.
Menganalisis hubungan sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, dan rasio tenaga kesehatan dengan risiko relatif kusta.
Mengidentifikasi kabupaten/kota dengan risiko relatif kusta tinggi sebagai wilayah prioritas pengendalian kusta.

2 Data dan Metode Penelitian

2.1 Data dan Variabel Penelitian

Penelitian ini menggunakan data sekunder berbentuk data panel pada 27 kabupaten/kota di Provinsi Jawa Barat selama periode 2019–2024. Unit analisis penelitian adalah kabupaten/kota per tahun, sehingga jumlah observasi yang digunakan adalah:

\[ 27 \text{ kabupaten/kota} \times 6 \text{ tahun} = 162 \text{ observasi} \]

Data penelitian meliputi jumlah kasus kusta, jumlah penduduk, sanitasi layak, akses air minum layak, jumlah fasilitas kesehatan, jumlah tenaga kesehatan, serta peta batas administrasi kabupaten/kota di Jawa Barat. Jumlah kasus kusta digunakan sebagai variabel respons pada model Poisson BYM2. Jumlah penduduk tidak digunakan sebagai kovariat, tetapi digunakan untuk menghitung insidensi, rasio fasilitas kesehatan, rasio tenaga kesehatan, dan expected cases sebagai offset dalam model.

Data dan Variabel Penelitian
No	Data/Variabel	Simbol	Sumber Data	Definisi Operasional	Satuan/Skala	Peran dalam Analisis
1	Jumlah kasus kusta	\(Y_{it}\)	Dinas Kesehatan Provinsi Jawa Barat / Profil Kesehatan Jawa Barat	Jumlah kasus kusta yang tercatat pada kabupaten/kota ke-i pada tahun ke-t.	Kasus; rasio	Variabel respons pada model Poisson BYM2
2	Jumlah penduduk	\(N_{it}\)	Badan Pusat Statistik Provinsi Jawa Barat	Jumlah penduduk pada kabupaten/kota ke-i pada tahun ke-t.	Jiwa; rasio	Menghitung insidensi, rasio layanan kesehatan, dan expected cases; bukan kovariat
3	Insidensi kusta	\(IR_{it}\)	Hasil perhitungan dari jumlah kasus dan jumlah penduduk	Jumlah kasus kusta per 100.000 penduduk.	Kasus per 100.000 penduduk; rasio	Analisis deskriptif dan peta insidensi
4	Expected cases	\(E_{it}\)	Hasil perhitungan dari jumlah kasus dan jumlah penduduk	Jumlah kasus kusta yang diharapkan berdasarkan jumlah penduduk wilayah dan laju kasus Jawa Barat pada tahun yang sama.	Kasus; rasio	Offset pada model Poisson BYM2
5	Risiko relatif	\(RR_{it}\)	Hasil estimasi model Bayesian spasio-temporal BYM2	Perbandingan risiko kusta suatu kabupaten/kota terhadap risiko rata-rata Jawa Barat setelah penyesuaian model.	Rasio	Keluaran utama model
6	Sanitasi layak	\(X_{1it}\)	Open Data Jawa Barat / Dinas terkait	Persentase penduduk atau rumah tangga yang memiliki akses sanitasi layak.	Persen; rasio	Kovariat lingkungan
7	Akses air minum layak	\(X_{2it}\)	Open Data Jawa Barat / Dinas terkait	Persentase penduduk atau rumah tangga yang memiliki akses air minum layak.	Persen; rasio	Kovariat lingkungan
8	Jumlah fasilitas kesehatan	\(X_{3it}\)	Dinas Kesehatan Provinsi Jawa Barat / Profil Kesehatan Jawa Barat	Jumlah fasilitas kesehatan pada setiap kabupaten/kota.	Fasilitas; rasio	Dibentuk menjadi rasio faskes per 100.000 penduduk sebagai kovariat
9	Jumlah tenaga kesehatan	\(X_{4it}\)	Dinas Kesehatan Provinsi Jawa Barat / Profil Kesehatan Jawa Barat	Jumlah tenaga kesehatan pada setiap kabupaten/kota.	Tenaga; rasio	Dibentuk menjadi rasio nakes per 100.000 penduduk sebagai kovariat
10	Wilayah	\(i\)	GeoJSON batas administrasi Jawa Barat	Identitas kabupaten/kota di Provinsi Jawa Barat.	Nominal	Indeks spasial pada model
11	Tahun	\(t\)	Data penelitian	Tahun pengamatan penelitian, yaitu 2019–2024.	Interval	Indeks temporal pada model

2.2 Metode Analisis Data

Analisis data dilakukan secara bertahap, meliputi persiapan data, analisis deskriptif, perhitungan ukuran epidemiologi, analisis spasial eksploratif, pembentukan matriks ketetanggaan, pemodelan Bayesian spasio-temporal Poisson BYM2, evaluasi model, dan pemetaan risiko relatif. Tahapan analisis disajikan pada Tabel berikut:

Tahapan Metode Analisis Data
No	Metode Analisis	Tujuan	Output
1	Persiapan dan pembersihan data	Memastikan data lengkap, konsisten, dan siap dianalisis.	Data panel 27 kabupaten/kota selama 2019–2024.
2	Analisis deskriptif	Menggambarkan karakteristik dan perkembangan data kusta serta kovariat penelitian.	Tabel ringkasan statistik dan grafik tren.
3	Perhitungan insidensi kusta	Menggambarkan tingkat kejadian kusta pada setiap kabupaten/kota.	Insidensi kusta per 100.000 penduduk.
4	Perhitungan expected cases dan SIR	Membandingkan jumlah kasus aktual dengan jumlah kasus yang diharapkan berdasarkan populasi.	Expected cases dan Standardized Incidence Ratio (SIR).
5	Analisis spasial eksploratif	Menggambarkan pola persebaran geografis insidensi dan SIR kusta.	Peta insidensi dan peta SIR.
6	Pembentukan matriks ketetanggaan	Menentukan hubungan spasial antarwilayah kabupaten/kota.	Matriks ketetanggaan queen contiguity.
7	Pemodelan Bayesian spasio-temporal Poisson BYM2	Mengestimasi risiko relatif kusta dengan mempertimbangkan populasi, spasial, temporal, dan kovariat.	Estimasi relative risk, efek kovariat, efek spasial, dan efek temporal.
8	Evaluasi model	Menilai kecocokan model dan ketidakpastian estimasi parameter.	DIC, WAIC, dan credible interval 95%.
9	Pemetaan risiko relatif dan wilayah prioritas	Mengidentifikasi kabupaten/kota dengan risiko relatif kusta tinggi.	Peta relative risk dan tabel wilayah prioritas.

2.2.1 Persiapan dan Pembersihan Data

Tahap persiapan data dilakukan dengan memeriksa kelengkapan data, konsistensi nama kabupaten/kota, dan nilai hilang pada setiap variabel. Nama kabupaten/kota pada data tabular disamakan dengan nama wilayah pada peta administrasi Jawa Barat. Data kemudian disusun dalam bentuk data panel dengan unit analisis kabupaten/kota per tahun.

Penelitian menggunakan 27 kabupaten/kota selama enam tahun, yaitu tahun 2019–2024, sehingga jumlah observasi adalah:

\[ 27 \text{ kabupaten/kota} \times 6 \text{ tahun} = 162 \text{ observasi} \]

2.2.2 Analisis Deskriptif

Analisis deskriptif dilakukan untuk menggambarkan jumlah kasus kusta, jumlah penduduk, sanitasi layak, akses air minum layak, fasilitas kesehatan, dan tenaga kesehatan. Analisis juga dilakukan untuk melihat tren insidensi kusta Jawa Barat selama tahun 2019–2024.

2.2.3 Perhitungan Insidensi Kusta

Insidensi kusta dihitung untuk menggambarkan tingkat kejadian kusta pada setiap kabupaten/kota. Rumus insidensi kusta adalah:

\[ IR_{it} = \frac{Y_{it}}{N_{it}} \times 100.000 \]

dengan:

\(IR_{it}\) adalah insidensi kusta pada wilayah ke-\(i\) tahun ke-\(t\);
\(Y_{it}\) adalah jumlah kasus kusta pada wilayah ke-\(i\) tahun ke-\(t\);
\(N_{it}\) adalah jumlah penduduk pada wilayah ke-\(i\) tahun ke-\(t\).

2.2.4 Perhitungan Expected Cases

Expected cases atau kasus harapan adalah jumlah kasus yang diperkirakan terjadi pada suatu wilayah apabila wilayah tersebut memiliki tingkat risiko yang sama dengan rata-rata Jawa Barat pada tahun yang sama.

\[ E_{it} = N_{it} \times \left( \frac{\sum_i Y_{it}}{\sum_i N_{it}} \right) \]

dengan \(E_{it}\) adalah jumlah kasus kusta yang diharapkan pada wilayah ke-\(i\) tahun ke-\(t\).

2.2.5 Standardized Incidence Ratio (SIR)

SIR membandingkan jumlah kasus aktual dengan jumlah kasus harapan, nilai Standardized Incidence Ratio (SIR) dihitung menggunakan rumus:

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

Interpretasi nilai SIR adalah sebagai berikut:

\(SIR_{it} > 1\) menunjukkan jumlah kasus aktual lebih tinggi daripada jumlah kasus harapan;
\(SIR_{it} = 1\) menunjukkan jumlah kasus aktual sama dengan jumlah kasus harapan;
\(SIR_{it} < 1\) menunjukkan jumlah kasus aktual lebih rendah daripada jumlah kasus harapan.

2.2.6 Analisis Spasial Eksploratif

Analisis spasial eksploratif dilakukan melalui peta tematik insidensi kusta dan SIR. Peta insidensi digunakan untuk menggambarkan tingkat kejadian kusta per 100.000 penduduk, sedangkan peta SIR digunakan untuk menunjukkan wilayah dengan jumlah kasus aktual yang lebih tinggi atau lebih rendah dibandingkan kasus harapan.

2.2.7 Pembentukan Matriks Ketetanggaan Spasial

Hubungan spasial antar kabupaten/kota di Jawa Barat dibentuk menggunakan metode queen contiguity. Dalam metode ini, dua wilayah dikategorikan sebagai bertetangga apabila memiliki batas sisi atau titik sudut yang bersinggungan.

Matriks ketetanggaan spasial dinyatakan sebagai:

\[ w_{ij} = \begin{cases} 1, & \text{jika wilayah } i \text{ dan wilayah } j \text{ bertetangga}, \\ 0, & \text{jika wilayah } i \text{ dan wilayah } j \text{ tidak bertetangga}. \end{cases} \]

dengan:

\(w_{ij}\) adalah nilai hubungan ketetanggaan antara wilayah ke-\(i\) dan wilayah ke-\(j\);
\(w_{ij}=1\) menunjukkan bahwa kedua wilayah bertetangga;
\(w_{ij}=0\) menunjukkan bahwa kedua wilayah tidak bertetangga;
\(i,j = 1,2,\ldots,27\) menunjukkan indeks kabupaten/kota di Jawa Barat.

Matriks ketetanggaan ini digunakan untuk membentuk komponen spasial terstruktur dalam model BYM2. Komponen tersebut memungkinkan estimasi risiko kusta pada suatu kabupaten/kota dipengaruhi oleh pola risiko pada kabupaten/kota yang berbatasan langsung dengannya. Penggunaan struktur ketetanggaan dalam pemetaan penyakit membantu menghasilkan estimasi risiko yang lebih stabil dibandingkan penggunaan rasio kasus mentah, terutama pada wilayah dengan jumlah penduduk atau jumlah kasus yang kecil.

2.2.8 7. Pemodelan Bayesian Spasio-Temporal Poisson BYM2

Jumlah kasus kusta dimodelkan menggunakan distribusi Poisson:

\[ Y_{it} \sim \text{Poisson}(\mu_{it}) \]

Nilai rata-rata jumlah kasus dinyatakan sebagai:

\[ \mu_{it} = E_{it} \times RR_{it} \]

Model risiko relatif kusta dirumuskan sebagai berikut:

\[ \log(RR_{it}) = \alpha + \beta_1 \text{Sanitasi}_{it} + \beta_2 \text{AirMinum}_{it} + \beta_3 \text{Faskes}_{it} + \beta_4 \text{Nakes}_{it} + u_i + v_i + \gamma_t + \delta_{it} \]

dengan:

\(\alpha\) adalah intercept;
\(\beta_1\), \(\beta_2\), \(\beta_3\), dan \(\beta_4\) adalah koefisien kovariat;
\(u_i\) adalah efek spasial terstruktur;
\(v_i\) adalah efek spasial tidak terstruktur;
\(\gamma_t\) adalah efek temporal menggunakan Random Walk orde satu (RW1);
\(\delta_{it}\) adalah efek interaksi spasio-temporal;
\(RR_{it}\) adalah risiko relatif kusta pada wilayah ke-\(i\) tahun ke-\(t\).

Komponen \(u_i\) dan \(v_i\) dimodelkan menggunakan pendekatan BYM2. Jumlah penduduk tidak digunakan sebagai kovariat, tetapi digunakan untuk membentuk expected cases sebagai offset dalam model Poisson.

2.2.9 8. Evaluasi Model

Evaluasi model dilakukan menggunakan nilai Deviance Information Criterion (DIC), Watanabe-Akaike Information Criterion (WAIC), dan credible interval 95%.

Nilai DIC dan WAIC digunakan untuk menilai kecocokan model. Nilai yang lebih kecil menunjukkan kecocokan model yang lebih baik apabila dilakukan perbandingan dengan model lain. Credible interval 95% digunakan untuk menggambarkan ketidakpastian estimasi parameter model.

2.2.10 9. Pemetaan Risiko Relatif dan Identifikasi Wilayah Prioritas

Hasil utama model adalah estimasi relative risk kusta pada setiap kabupaten/kota dan tahun pengamatan. Risiko relatif dihitung sebagai:

\[ RR_{it} = \frac{\mu_{it}}{E_{it}} \]

dengan \(\mu_{it}\) merupakan estimasi jumlah kasus dari model dan \(E_{it}\) merupakan jumlah kasus harapan.

Kabupaten/kota dengan nilai:

\[ RR_{it} > 1 \]

diinterpretasikan memiliki risiko kusta lebih tinggi dibandingkan risiko rata-rata Jawa Barat. Wilayah dengan nilai relative risk tertinggi pada tahun 2024 diidentifikasi sebagai wilayah prioritas untuk penguatan surveilans, penemuan kasus aktif, pemeriksaan kontak, dan pengendalian kusta.

3 Hasil dan Pembahasan

3.1 1. Persiapan Data

Data penelitian terdiri atas 27 kabupaten/kota di Jawa Barat selama tahun 2019–2024. Dengan enam tahun pengamatan, jumlah observasi yang dianalisis adalah 162 observasi. Tahap persiapan data dilakukan melalui pemeriksaan struktur data, kelengkapan data, duplikasi, dan konsistensi data kabupaten/kota-tahun.

Pemeriksaan Struktur Data Penelitian
Keterangan	Nilai
Jumlah observasi	162
Jumlah kabupaten/kota	27
Jumlah tahun	6
Rentang tahun	2019–2024
Jumlah data duplikat	0

3.2 Statistik Deskriptif

Statistik Deskriptif Variabel Penelitian di Jawa Barat Tahun 2019–2024
Variabel	Jumlah Data	Rata_rata	Standar Deviasi	Minimum	Kuartil 1	Median	Kuartil 3	Maksimum
Akses Air Minum Layak (%)	162	93.75	6.02	76.18	91.17	95.86	98.18	100
Jumlah Fasilitas Kesehatan	162	40.69	21.06	10.00	22.00	38.00	50.00	101
Jumlah Kasus Kusta	162	58.74	82.50	0.00	5.00	18.00	89.25	382
Jumlah Penduduk	162	1825432.84	1198583.83	183100.00	1037945.00	1768145.00	2521475.00	5965400
Jumlah Tenaga Kesehatan	162	98.70	59.05	1.00	58.50	88.00	120.00	310
Sanitasi Layak (%)	162	84.89	18.79	0.00	76.62	90.06	96.33	173

Berdasarkan Tabel Statistik Deskriptif, penelitian ini menggunakan 162 observasi yang berasal dari 27 kabupaten/kota di Jawa Barat selama periode 2019–2024. Akses air minum layak memiliki rata-rata sebesar 93,75% dengan standar deviasi 6,02%, nilai minimum 76,18%, dan nilai maksimum 100%. Hasil tersebut menunjukkan bahwa akses air minum layak di Jawa Barat secara umum sudah cukup tinggi dan relatif merata, meskipun masih terdapat perbedaan antarwilayah dan antarwaktu.

Jumlah fasilitas kesehatan memiliki rata-rata 40,69 unit dengan standar deviasi 21,06 unit. Nilai minimum sebesar 10 unit dan maksimum 101 unit menunjukkan adanya variasi ketersediaan fasilitas kesehatan yang cukup besar antar kabupaten/kota. Variasi tersebut dapat berkaitan dengan perbedaan jumlah penduduk, luas wilayah, tingkat urbanisasi, dan kapasitas pelayanan kesehatan di masing-masing wilayah.

Jumlah kasus kusta memiliki rata-rata 58,74 kasus dengan standar deviasi 82,50 kasus. Nilai standar deviasi yang lebih besar daripada rata-rata menunjukkan bahwa jumlah kasus kusta memiliki keragaman yang tinggi antarwilayah dan antarwaktu. Jumlah kasus minimum sebesar 0 kasus dan maksimum sebesar 382 kasus memperlihatkan bahwa terdapat wilayah dengan tidak adanya kasus tercatat, sedangkan wilayah lain memiliki jumlah kasus yang jauh lebih tinggi. Oleh karena itu, jumlah kasus mentah tidak dapat digunakan secara langsung untuk membandingkan risiko antarwilayah, sehingga diperlukan perhitungan insidensi, SIR, dan risiko relatif hasil model.

Jumlah penduduk memiliki rata-rata sekitar 1.825.432 jiwa dengan standar deviasi sekitar 1.198.584 jiwa. Nilai minimum sebesar 183.100 jiwa dan maksimum sebesar 5.965.400 jiwa menunjukkan adanya perbedaan ukuran populasi yang sangat besar antar kabupaten/kota. Oleh sebab itu, jumlah penduduk digunakan untuk menghitung insidensi dan expected cases sebagai penyesuaian dalam model Poisson, bukan dimasukkan sebagai kovariat penelitian.

Jumlah tenaga kesehatan memiliki rata-rata 98,70 orang dengan standar deviasi 59,05 orang. Nilai minimum sebesar 1 orang dan maksimum sebesar 310 orang menunjukkan adanya ketimpangan ketersediaan tenaga kesehatan antarwilayah. Agar perbandingan antarwilayah menjadi lebih proporsional, jumlah fasilitas kesehatan dan tenaga kesehatan selanjutnya dikonversi menjadi rasio per 100.000 penduduk sebelum digunakan sebagai kovariat dalam model.

Sanitasi layak memiliki rata-rata sebesar 84,89% dengan standar deviasi 18,79%. Namun, nilai maksimum sebesar 173% menunjukkan adanya nilai yang perlu diperiksa kembali karena secara konsep persentase seharusnya berada pada rentang 0% hingga 100%. Nilai tersebut dapat berasal dari perbedaan definisi indikator, kesalahan pencatatan, atau ketidaksesuaian satuan data. Oleh karena itu, data sanitasi layak perlu divalidasi atau dibersihkan terlebih dahulu sebelum digunakan dalam pemodelan Bayesian spasio-temporal.

3.3 Perhitungan Ukuran Epidemiologi

##  [1] "nama_kabupaten_kota"  "tahun"                "jumlah_kasus"        
##  [4] "penduduk"             "jumlah_faskes"        "jumlah_nakes"        
##  [7] "sanitasi_layak"       "air_minum_layak"      "total_kasus_jabar"   
## [10] "total_penduduk_jabar" "incidence_rate"       "rate_jabar"          
## [13] "expected_cases"       "SIR"                  "rasio_faskes_100rb"  
## [16] "rasio_nakes_100rb"

3.3.1 Ukuran kejadian kusta

Ukuran Kejadian Kusta di Jawa Barat Tahun 2019–2024
Indikator	Nilai
Total Kasus Kusta	9516.00
Total Penduduk	295720120.00
Rata-rata Insidensi per 100.000 Penduduk	2.86
Standar Deviasi Insidensi	3.37
Minimum Insidensi	0.00
Kuartil 1 Insidensi	0.51
Median Insidensi	1.41
Kuartil 3 Insidensi	4.10
Maksimum Insidensi	20.12

Berdasarkan Tabel Ukuran Kejadian Kusta di Jawa Barat Tahun 2019–2024, tercatat total 9.516 kasus kusta selama periode penelitian. Rata-rata insidensi kusta sebesar 2,86 kasus per 100.000 penduduk, dengan median sebesar 1,41 kasus per 100.000 penduduk. Perbedaan antara rata-rata dan median menunjukkan bahwa sebagian besar kabupaten/kota memiliki insidensi yang relatif rendah, tetapi terdapat beberapa wilayah dengan insidensi yang lebih tinggi. Nilai minimum insidensi sebesar 0,00 menunjukkan adanya kabupaten/kota-tahun yang tidak mencatat kasus kusta, sedangkan nilai maksimum sebesar 20,12 kasus per 100.000 penduduk menunjukkan adanya wilayah dengan tingkat kejadian yang jauh lebih tinggi dibandingkan rata-rata.

Standar deviasi insidensi sebesar 3,37 menunjukkan adanya variasi insidensi yang cukup besar antar kabupaten/kota dan antar tahun. Kuartil pertama sebesar 0,51, median sebesar 1,41, dan kuartil ketiga sebesar 4,10 memperlihatkan bahwa 75% observasi memiliki insidensi tidak lebih dari 4,10 kasus per 100.000 penduduk, sementara sebagian kecil observasi memiliki nilai yang jauh lebih tinggi. Pola distribusi tersebut menunjukkan adanya heterogenitas risiko kusta di Jawa Barat, sehingga analisis menggunakan model Bayesian spasio-temporal BYM2 diperlukan untuk mengidentifikasi wilayah dengan risiko relatif tinggi secara lebih stabil setelah memperhitungkan perbedaan jumlah penduduk, hubungan spasial antarwilayah, dan perubahan risiko antarwaktu.

3.3.2 2. Kasus harapan dan SIR

Kasus Harapan dan Standardized Incidence Ratio Kusta Tahun 2019–2024
Indikator	Nilai
Rata-rata Expected Cases	58.74
Minimum Expected Cases	5.01
Maksimum Expected Cases	254.14
Rata-rata SIR	0.89
Standar Deviasi SIR	1.03
Minimum SIR	0.00
Kuartil 1 SIR	0.17
Median SIR	0.43
Kuartil 3 SIR	1.34
Maksimum SIR	5.76

Berdasarkan tabel kasus harapan dan Standardized Incidence Ratio (SIR), rata-rata jumlah kasus harapan kusta pada setiap kabupaten/kota-tahun adalah 58,74 kasus, dengan nilai minimum 5,01 kasus dan maksimum 254,14 kasus. Perbedaan tersebut menunjukkan bahwa jumlah kasus yang diharapkan pada setiap wilayah tidak sama karena dipengaruhi oleh perbedaan jumlah penduduk dan laju kasus kusta Jawa Barat pada tahun yang bersangkutan. Wilayah dengan jumlah penduduk lebih besar cenderung memiliki jumlah kasus harapan yang lebih tinggi.

Rata-rata SIR sebesar 0,89 menunjukkan bahwa secara keseluruhan jumlah kasus aktual cenderung sedikit lebih rendah dibandingkan jumlah kasus yang diharapkan berdasarkan rata-rata Jawa Barat. Namun, standar deviasi SIR sebesar 1,03 serta rentang nilai dari 0,00 hingga 5,76 menunjukkan adanya variasi yang besar antarwilayah dan antarwaktu. Median SIR sebesar 0,43 berarti lebih dari separuh observasi memiliki jumlah kasus aktual di bawah kasus harapan, sedangkan kuartil ketiga sebesar 1,34 menunjukkan bahwa sekitar 25% observasi memiliki SIR di atas 1,34. Nilai maksimum SIR sebesar 5,76 mengindikasikan adanya wilayah tertentu dengan jumlah kasus aktual hingga sekitar 5,76 kali lebih tinggi daripada jumlah kasus harapan. Variasi SIR yang cukup besar ini mendukung penggunaan model Bayesian spasio-temporal untuk memperoleh estimasi risiko relatif yang lebih stabil.

3.3.3 3. Ketersediaan layanan kesehatan

Rasio Ketersediaan Layanan Kesehatan di Jawa Barat Tahun 2019–2024
Indikator	Nilai
Rata-rata Rasio Faskes per 100.000 Penduduk	2.68
Standar Deviasi Rasio Faskes	1.11
Minimum Rasio Faskes	1.20
Median Rasio Faskes	2.41
Maksimum Rasio Faskes	6.89
Rata-rata Rasio Nakes per 100.000 Penduduk	7.41
Standar Deviasi Rasio Nakes	5.95
Minimum Rasio Nakes	0.05
Median Rasio Nakes	6.11
Maksimum Rasio Nakes	38.13

Berdasarkan tabel rasio ketersediaan layanan kesehatan, rata-rata rasio fasilitas kesehatan di Jawa Barat selama tahun 2019–2024 adalah 2,68 fasilitas kesehatan per 100.000 penduduk. Median sebesar 2,41 menunjukkan bahwa sebagian besar kabupaten/kota memiliki rasio faskes yang berada di sekitar nilai tersebut. Namun, nilai minimum 1,20 dan maksimum 6,89, dengan standar deviasi 1,11, menunjukkan adanya perbedaan ketersediaan fasilitas kesehatan antarwilayah. Wilayah dengan rasio faskes rendah dapat memiliki keterbatasan akses layanan dibandingkan wilayah dengan rasio yang lebih tinggi.

Rata-rata rasio tenaga kesehatan adalah 7,41 tenaga kesehatan per 100.000 penduduk, dengan median 6,11. Variasi rasio tenaga kesehatan lebih besar dibandingkan rasio fasilitas kesehatan, yang terlihat dari standar deviasi sebesar 5,95, nilai minimum 0,05, dan maksimum 38,13. Hal ini menunjukkan adanya ketimpangan yang cukup besar dalam distribusi tenaga kesehatan antar kabupaten/kota dan antar tahun. Perbedaan ketersediaan fasilitas dan tenaga kesehatan tersebut relevan dimasukkan sebagai kovariat dalam model Bayesian spasio-temporal, karena dapat berkaitan dengan kemampuan deteksi, pelaporan, diagnosis, dan penanganan kasus kusta di setiap wilayah.

3.4 Tren Insidensi Kusta di Jawa Barat

tren_insidensi_hasil <- data_model %>%
  group_by(tahun) %>%
  summarise(
    total_kasus = sum(jumlah_kasus, na.rm = TRUE),
    total_penduduk = sum(penduduk, na.rm = TRUE),
    insidensi_jabar = ifelse(
      total_penduduk > 0,
      (total_kasus / total_penduduk) * 100000,
      NA_real_
    ),
    .groups = "drop"
  )

ggplot(tren_insidensi_hasil, aes(x = tahun, y = insidensi_jabar)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_text(
    aes(label = round(insidensi_jabar, 2)),
    vjust = -0.8,
    size = 3.5
  ) +
  scale_x_continuous(breaks = 2019:2024) +
  labs(
    x = "Tahun",
    y = "Insidensi per 100.000 Penduduk"
  ) +
  theme_minimal()

Grafik menunjukkan bahwa insidensi kusta di Jawa Barat selama tahun 2019–2024 mengalami pola yang berfluktuasi. Insidensi tertinggi terjadi pada tahun 2019, yaitu sebesar 4,26 kasus per 100.000 penduduk. Setelah itu, insidensi menurun cukup tajam menjadi 2,92 pada tahun 2020 dan kembali menurun menjadi 2,70 pada tahun 2021. Penurunan pada periode 2020–2021 dapat mencerminkan berkurangnya kasus yang tercatat, tetapi juga perlu ditafsirkan secara hati-hati karena dapat dipengaruhi oleh perubahan akses pelayanan, penemuan kasus, atau pelaporan selama periode tersebut.

Pada tahun 2022, insidensi kembali meningkat menjadi 3,49 kasus per 100.000 penduduk dan relatif stabil pada tahun 2023 sebesar 3,55 kasus per 100.000 penduduk. Selanjutnya, insidensi menurun menjadi 2,39 kasus per 100.000 penduduk pada tahun 2024, yang merupakan nilai terendah selama periode pengamatan. Secara keseluruhan, grafik menunjukkan kecenderungan penurunan insidensi dibandingkan tahun awal penelitian, tetapi perubahan antar tahun yang tidak konsisten menunjukkan bahwa risiko kusta masih bervariasi secara temporal. Oleh karena itu, analisis spasio-temporal diperlukan untuk mengidentifikasi apakah perubahan tersebut terjadi merata di seluruh kabupaten/kota atau hanya terkonsentrasi pada wilayah tertentu.

3.5 Analisis Spasial Eksploratif

3.5.1 Peta Insidensi Kusta Tahun 2019–2024

Peta Insidensi Kusta di Jawa Barat Tahun 2019–2024

Peta insidensi menunjukkan bahwa persebaran kejadian kusta di Jawa Barat selama tahun 2019–2024 tidak merata antar kabupaten/kota. Wilayah dengan kategori insidensi tinggi hingga sangat tinggi, yang ditunjukkan oleh warna biru tua, cenderung terkonsentrasi di bagian utara dan timur Jawa Barat. Pola tersebut tampak berulang pada hampir seluruh tahun pengamatan, terutama pada tahun 2019, 2022, 2023, dan 2024. Sementara itu, sebagian wilayah di bagian tengah dan selatan Jawa Barat cenderung berada pada kategori insidensi rendah, yang ditunjukkan oleh warna lebih terang.

Secara temporal, peta memperlihatkan adanya perubahan intensitas insidensi dari tahun ke tahun, tetapi pola konsentrasi wilayah dengan insidensi tinggi relatif tetap pada beberapa wilayah tertentu. Tahun 2020 dan 2021 terlihat memiliki intensitas insidensi yang cenderung lebih rendah dibandingkan tahun 2019, 2022, dan 2023. Pada tahun 2024, beberapa wilayah masih menunjukkan kategori insidensi tinggi meskipun rata-rata insidensi Jawa Barat secara keseluruhan menurun. Temuan ini menunjukkan adanya heterogenitas spasial dan temporal pada kejadian kusta, sehingga penggunaan model Bayesian spasio-temporal BYM2 relevan untuk mengestimasi risiko relatif yang lebih stabil dengan mempertimbangkan keterkaitan antarwilayah dan perubahan risiko antarwaktu.

3.6 Peta SIR

Peta Standardized Incidence Ratio Kusta di Jawa Barat Tahun 2019–2024

Peta Standardized Incidence Ratio (SIR) menunjukkan perbandingan antara jumlah kasus kusta aktual dan jumlah kasus yang diharapkan pada setiap kabupaten/kota setelah mempertimbangkan perbedaan jumlah penduduk. Wilayah dengan warna biru tua memiliki nilai SIR lebih tinggi, yang berarti jumlah kasus aktual di wilayah tersebut lebih besar dibandingkan jumlah kasus yang diharapkan berdasarkan rata-rata Jawa Barat. Sebaliknya, wilayah dengan warna lebih terang memiliki nilai SIR lebih rendah, yang menunjukkan jumlah kasus aktual lebih rendah daripada kasus harapan. Selama tahun 2019–2024, wilayah dengan SIR tinggi tampak cenderung terkonsentrasi pada bagian utara hingga timur Jawa Barat, sedangkan beberapa wilayah di bagian tengah dan selatan relatif berada pada kategori SIR rendah.

Secara temporal, pola SIR terlihat relatif konsisten dari tahun ke tahun, meskipun terdapat perubahan tingkat intensitas pada beberapa kabupaten/kota. Pada tahun 2019, 2022, 2023, dan 2024, beberapa wilayah di bagian utara dan timur secara berulang berada pada kategori SIR tinggi. Hal ini menunjukkan bahwa wilayah tersebut memiliki jumlah kasus kusta yang secara konsisten lebih tinggi dibandingkan yang diharapkan. Namun, SIR masih merupakan ukuran kasar yang dapat berfluktuasi, terutama pada wilayah dengan jumlah penduduk atau kasus kecil. Oleh karena itu, peta SIR digunakan sebagai eksplorasi awal, sedangkan penentuan risiko relatif dan wilayah prioritas utama dilakukan berdasarkan hasil model Bayesian spasio-temporal BYM2.

3.7 Matrix Bobot Spasial

peta_sp_hasil <- as(peta_jabar, "Spatial")

nb_queen_hasil <- poly2nb(peta_sp_hasil, queen = TRUE)

ringkasan_tetangga_hasil <- data.frame(
  Keterangan = c(
    "Jumlah wilayah",
    "Jumlah hubungan ketetanggaan",
    "Rata-rata jumlah tetangga per wilayah"
  ),
  Nilai = c(
    length(nb_queen_hasil),
    sum(card(nb_queen_hasil)),
    round(mean(card(nb_queen_hasil)), 2)
  )
)

kable(
  ringkasan_tetangga_hasil,
  caption = "Ringkasan Matriks Ketetanggaan Queen Contiguity",
  booktabs = TRUE,
  align = c("l", "r")
) %>%
  kable_styling(
    full_width = FALSE,
    latex_options = "hold_position"
  )

Ringkasan Matriks Ketetanggaan Queen Contiguity
Keterangan	Nilai
Jumlah wilayah	27.00
Jumlah hubungan ketetanggaan	106.00
Rata-rata jumlah tetangga per wilayah	3.93

Berdasarkan ringkasan matriks ketetanggaan queen contiguity, analisis spasial melibatkan 27 kabupaten/kota di Jawa Barat. Metode queen contiguity menetapkan dua wilayah sebagai bertetangga apabila keduanya berbagi batas sisi maupun titik sudut. Dari pembentukan matriks tersebut diperoleh 106 hubungan ketetanggaan, yang menunjukkan bahwa antarwilayah di Jawa Barat memiliki keterhubungan spasial yang cukup luas.

Rata-rata jumlah tetangga per wilayah adalah 3,93 wilayah. Artinya, setiap kabupaten/kota rata-rata memiliki sekitar tiga hingga empat wilayah yang berbatasan langsung. Struktur ketetanggaan ini digunakan dalam model BYM2 untuk membentuk efek spasial terstruktur, sehingga estimasi risiko kusta pada suatu wilayah dapat memanfaatkan informasi dari wilayah-wilayah di sekitarnya. Dengan demikian, model dapat menangkap kemungkinan adanya kemiripan risiko kusta pada kabupaten/kota yang saling berdekatan.

3.8 Hasil Pemodelan Bayessian Spatio-Temporal Poisson (BYM2)

Hasil Estimasi Parameter Model Bayesian Spasio-Temporal Poisson BYM2 dengan Kovariat
Variabel	Estimasi	Batas Bawah 95%	Batas Atas 95%	Relative Risk	RR Bawah 95%	RR Atas 95%
Intercept	-0.8436	-1.0660	-0.6197	0.430	0.344	0.538
Sanitasi Layak	-0.0393	-0.1235	0.0435	0.961	0.884	1.044
Akses Air Minum Layak	0.0579	-0.0927	0.2080	1.060	0.911	1.231
Rasio Fasilitas Kesehatan	-0.0930	-0.5467	0.3515	0.911	0.579	1.421
Rasio Tenaga Kesehatan	-0.0869	-0.2248	0.0529	0.917	0.799	1.054

3.8.1 Model Estimasi Bayesian Spasio-Temporal Poisson BYM2

Berdasarkan hasil estimasi parameter model Bayesian spasio-temporal Poisson BYM2 dengan kovariat, model estimasi risiko kusta dapat dituliskan sebagai berikut:

\[ \begin{aligned} \log(\hat{\mu}_{it}) =\;& \log(E_{it}) - 0{,}8436 \\ &- 0{,}0393(\text{Sanitasi}_{it}) \\ &+ 0{,}0579(\text{AirMinum}_{it}) \\ &- 0{,}0930(\text{RasioFaskes}_{it}) \\ &- 0{,}0869(\text{RasioNakes}_{it}) \\ &+ b_i + \gamma_t + \delta_{it}. \end{aligned} \]

Model tersebut dapat juga dinyatakan dalam bentuk risiko relatif sebagai berikut:

\[ \begin{aligned} RR_{it} = \exp \Big[ &-0{,}8436 -0{,}0393(\text{Sanitasi}_{it}) \\ &+0{,}0579(\text{AirMinum}_{it}) \\ &-0{,}0930(\text{RasioFaskes}_{it}) \\ &-0{,}0869(\text{RasioNakes}_{it}) \\ &+b_i+\gamma_t+\delta_{it} \Big]. \end{aligned} \]

dengan:

\(\hat{\mu}_{it}\) adalah estimasi jumlah kasus kusta pada kabupaten/kota ke-\(i\) pada tahun ke-\(t\);
\(E_{it}\) adalah jumlah kasus harapan (expected cases) pada kabupaten/kota ke-\(i\) pada tahun ke-\(t\);
\(\log(E_{it})\) merupakan offset untuk menyesuaikan perbedaan jumlah penduduk antarwilayah;
\(\text{Sanitasi}_{it}\) adalah nilai sanitasi layak yang telah distandardisasi;
\(\text{AirMinum}_{it}\) adalah nilai akses air minum layak yang telah distandardisasi;
\(\text{RasioFaskes}_{it}\) adalah rasio fasilitas kesehatan per 100.000 penduduk yang telah distandardisasi;
\(\text{RasioNakes}_{it}\) adalah rasio tenaga kesehatan per 100.000 penduduk yang telah distandardisasi;
\(b_i\) adalah efek spasial wilayah menggunakan model BYM2;
\(\gamma_t\) adalah efek temporal pada tahun ke-\(t\);
\(\delta_{it}\) adalah efek interaksi spasio-temporal;
\(RR_{it}\) adalah risiko relatif kusta pada wilayah ke-\(i\) dan tahun ke-\(t\).

Karena seluruh kovariat telah distandardisasi menggunakan z-score, maka setiap koefisien menggambarkan perubahan risiko relatif kusta akibat peningkatan satu standar deviasi pada masing-masing kovariat. Misalnya, koefisien sanitasi layak sebesar \(-0{,}0393\) menunjukkan bahwa peningkatan satu standar deviasi sanitasi layak cenderung menurunkan log risiko relatif kusta, sedangkan koefisien akses air minum layak sebesar \(0{,}0579\) menunjukkan kecenderungan peningkatan log risiko relatif kusta. Namun, interpretasi pengaruh kovariat tetap perlu mempertimbangkan interval kredibel 95% dari masing-masing parameter.

Hasil estimasi menunjukkan bahwa sanitasi layak memiliki relative risk sebesar 0,961. Artinya, peningkatan satu standar deviasi pada persentase sanitasi layak cenderung berkaitan dengan penurunan risiko relatif kusta sekitar 3,9% (1−0,961)×100%. Namun, interval kredibel 95% untuk relative risk berada pada rentang 0,884–1,044 dan masih memuat nilai 1. Dengan demikian, belum terdapat bukti posterior yang cukup kuat bahwa sanitasi layak berhubungan dengan risiko relatif kusta pada data penelitian ini.

Akses air minum layak memiliki relative risk sebesar 1,060, yang menunjukkan kecenderungan peningkatan risiko relatif kusta sekitar 6,0% untuk setiap peningkatan satu standar deviasi akses air minum layak. Akan tetapi, interval kredibel 95% sebesar 0,911–1,231 juga memuat nilai 1. Oleh karena itu, hubungan akses air minum layak dengan risiko relatif kusta belum dapat dinyatakan meyakinkan. Arah hubungan yang positif ini tidak langsung berarti bahwa air minum layak meningkatkan risiko kusta, melainkan dapat mencerminkan perbedaan pelaporan, penemuan kasus, karakteristik wilayah, atau faktor lain yang belum tercakup dalam model.

Rasio fasilitas kesehatan memiliki relative risk sebesar 0,911, yang menunjukkan kecenderungan penurunan risiko relatif sekitar 8,9% untuk setiap peningkatan satu standar deviasi rasio fasilitas kesehatan. Namun, interval kredibel 95% sebesar 0,579–1,421 masih memuat nilai 1. Rasio tenaga kesehatan juga memiliki relative risk sebesar 0,917, yang mengarah pada kecenderungan penurunan risiko relatif sekitar 8,3%, tetapi interval kredibel 95% sebesar 0,799–1,054 masih memuat nilai 1. Dengan demikian, kedua variabel layanan kesehatan tersebut belum menunjukkan hubungan yang meyakinkan secara posterior dengan risiko relatif kusta.

Secara keseluruhan, seluruh kovariat memiliki interval kredibel 95% relative risk yang memuat nilai 1. Artinya, berdasarkan data tahun 2019–2024 dan struktur model yang digunakan, belum terdapat bukti kuat bahwa sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, maupun rasio tenaga kesehatan berhubungan secara langsung dengan risiko relatif kusta. Meskipun demikian, kovariat tetap relevan dipertahankan dalam model sebagai faktor kontekstual, sedangkan variasi risiko kusta antarwilayah dan antarwaktu lebih lanjut ditangkap oleh komponen spasial, temporal, dan interaksi spasio-temporal pada model BYM2.

3.9 Evaluasi Kecocokan Model

evaluasi_model_hasil <- data.frame(
  Indikator = c("DIC", "WAIC"),
  Nilai = c(
    model_bym2_st$dic$dic,
    model_bym2_st$waic$waic
  )
) %>%
  mutate(
    Nilai = round(Nilai, 3)
  )

kable(
  evaluasi_model_hasil,
  caption = "Evaluasi Kecocokan Model Bayesian Spasio-Temporal Poisson BYM2",
  booktabs = TRUE,
  align = c("l", "r")
) %>%
  kable_styling(
    full_width = FALSE,
    latex_options = "hold_position"
  )

Evaluasi Kecocokan Model Bayesian Spasio-Temporal Poisson BYM2
Indikator	Nilai
DIC	1072.650
WAIC	1073.839

Nilai DIC sebesar 1072,651 dan WAIC sebesar 1073,840 merupakan ukuran evaluasi kecocokan model Bayesian Spasio-Temporal Poisson BYM2. Kedua indikator tersebut digunakan untuk menilai keseimbangan antara kemampuan model dalam menjelaskan data kasus kusta dan kompleksitas model yang dibentuk. Nilai DIC dan WAIC tidak memiliki batas absolut untuk menyatakan bahwa suatu model sudah baik atau buruk. Oleh karena itu, nilai tersebut tidak diinterpretasikan secara terpisah, tetapi umumnya digunakan untuk membandingkan beberapa model dengan data yang sama.

Evaluasi Kecocokan Model Bayesian Spasio-Temporal Poisson BYM2
Indikator	Nilai
Mean Absolute Error (MAE)	1.783
Root Mean Squared Error (RMSE)	2.492

Nilai MAE sebesar 1,783 menunjukkan bahwa rata-rata selisih absolut antara jumlah kasus kusta aktual dan jumlah kasus hasil prediksi model adalah sekitar 1,78 kasus pada setiap observasi kabupaten/kota-tahun. Nilai RMSE sebesar 2,492 menunjukkan bahwa kesalahan prediksi model, dengan penekanan lebih besar pada kesalahan yang besar, berada pada kisaran 2,49 kasus. Nilai RMSE yang sedikit lebih besar daripada MAE menunjukkan adanya beberapa observasi dengan kesalahan prediksi yang lebih tinggi, tetapi perbedaannya tidak besar sehingga kesalahan prediksi ekstrem tidak dominan.

Grafik perbandingan kasus aktual dan prediksi menunjukkan bahwa hampir seluruh titik berada sangat dekat dengan garis diagonal. Hal ini berarti nilai prediksi model BYM2 secara umum sangat mendekati jumlah kasus aktual. Dengan demikian, berdasarkan MAE, RMSE, dan pola grafik aktual–prediksi, model Bayesian Spasio-Temporal Poisson BYM2 memiliki kecocokan prediksi yang baik terhadap data kasus kusta tahun 2019–2024.

3.10 Peta Risiko Relatif

data_peta_rr_semua_tahun <- peta_jabar %>%
  left_join(
    data_model %>%
      select(nama_kabupaten_kota, tahun, relative_risk_model),
    by = "nama_kabupaten_kota"
  )

tm_shape(data_peta_rr_semua_tahun) +
  tm_polygons(
    col = "relative_risk_model",
    style = "quantile",
    n = 5,
    title = "Relative Risk"
  ) +
  tm_borders() +
  tm_facets(
    by = "tahun",
    ncol = 3,
    free.scales = FALSE
  ) +
  tm_layout(
    frame = FALSE,
    legend.outside = TRUE
  )

Peta risiko relatif hasil model Bayesian Spasio-Temporal Poisson BYM2 menunjukkan bahwa risiko kusta di Jawa Barat tidak tersebar merata antar kabupaten/kota maupun antar tahun. Wilayah dengan warna biru tua memiliki risiko relatif lebih tinggi, yaitu jumlah kasus kusta yang diperkirakan model lebih besar dibandingkan jumlah kasus harapan. Sebaliknya, wilayah dengan warna lebih terang memiliki risiko relatif lebih rendah. Selama periode 2019–2024, wilayah dengan risiko relatif tinggi cenderung terkonsentrasi di bagian utara hingga timur Jawa Barat, sedangkan sebagian wilayah di bagian tengah dan selatan cenderung memiliki risiko relatif rendah.

Secara temporal, pola wilayah berisiko tinggi terlihat cukup konsisten dari tahun ke tahun, terutama pada beberapa kabupaten/kota di bagian utara dan timur Jawa Barat yang berulang kali berada pada kategori risiko relatif tinggi. Meskipun demikian, terdapat perubahan intensitas risiko pada beberapa wilayah di setiap tahun. Peta risiko relatif ini lebih stabil dibandingkan peta insidensi dan SIR karena estimasinya telah mempertimbangkan perbedaan jumlah penduduk melalui expected cases, hubungan ketetanggaan antarwilayah melalui komponen BYM2, perubahan antar tahun, serta pengaruh kovariat. Oleh karena itu, wilayah dengan risiko relatif tinggi dapat dijadikan prioritas dalam penguatan penemuan kasus, pemeriksaan kontak, edukasi, dan pelayanan pengendalian kusta.

3.11 Wilayah Prioritas

Peringkat Risiko Relatif Kusta Seluruh Kabupaten/Kota di Jawa Barat Tahun 2024
Peringkat	nama_kabupaten_kota	jumlah_kasus	penduduk	incidence_rate	SIR	relative_risk_model	rr_lower	rr_upper
1	KABUPATENUPATEN INDRAMAYU	169	1914040	8.829	3.698	3.697	3.184	4.268
2	KABUPATENUPATEN BEKASI	198	3273870	6.048	2.533	2.515	2.190	2.875
3	KABUPATENUPATEN SUBANG	87	1663160	5.231	2.191	2.172	1.769	2.638
4	KABUPATENUPATEN CIREBON	120	2387960	5.025	2.105	2.109	1.769	2.494
5	KABUPATENUPATEN KARAWANG	111	2554380	4.345	1.820	1.827	1.523	2.173
6	KABUPATENUPATEN KUNINGAN	44	1213930	3.625	1.518	1.537	1.168	1.986
7	KOTA BEKASI	93	2644060	3.517	1.473	1.450	1.189	1.752
8	KOTA CIREBON	10	344850	2.900	1.215	1.259	0.758	1.970
9	KABUPATENUPATEN BOGOR	158	5682300	2.781	1.165	1.173	1.006	1.360
10	KABUPATENUPATEN PURWAKARTA	27	1050340	2.571	1.077	1.078	0.768	1.473
11	KOTA DEPOK	46	2163640	2.126	0.890	0.824	0.622	1.073
12	KABUPATENUPATEN MAJALENGKA	15	1352540	1.109	0.465	0.568	0.383	0.809
13	KABUPATENUPATEN SUKABUPATENUMI	41	2828020	1.450	0.607	0.538	0.397	0.715
14	KOTA BOGOR	7	1078350	0.649	0.272	0.334	0.201	0.522
15	KABUPATENUPATEN SUMEDANG	8	1187130	0.674	0.282	0.331	0.203	0.507
16	KABUPATENUPATEN PANGANDARAN	3	434100	0.691	0.289	0.329	0.168	0.583
17	KABUPATENUPATEN TASIKMALAYA	12	1920920	0.625	0.262	0.283	0.182	0.419
18	KABUPATENUPATEN CIANJUR	18	2584990	0.696	0.292	0.267	0.177	0.388
19	KOTA TASIKMALAYA	5	750730	0.666	0.279	0.242	0.130	0.413
20	KOTA SUKABUPATENUMI	0	365740	0.000	0.000	0.181	0.080	0.351
21	KABUPATENUPATEN CIAMIS	3	1259230	0.238	0.100	0.168	0.093	0.279
22	KOTA BANJAR	1	209790	0.477	0.200	0.150	0.050	0.351
23	KABUPATENUPATEN GARUT	12	2716950	0.442	0.185	0.149	0.092	0.231
24	KABUPATENUPATEN BANDUNG BARAT	5	1884190	0.265	0.111	0.098	0.054	0.165
25	KABUPATENUPATEN BANDUNG	5	3753120	0.133	0.056	0.051	0.028	0.086
26	KOTA BANDUNG	4	2528160	0.158	0.066	0.050	0.025	0.088
27	KOTA CIMAHI	0	598700	0.000	0.000	0.035	0.011	0.085

Pada tahun 2024, Kabupaten Indramayu merupakan wilayah dengan risiko relatif kusta tertinggi, yaitu 3,697 kali dibandingkan rata-rata Jawa Barat, diikuti Kabupaten Bekasi, Subang, Cirebon, dan Karawang. Wilayah tersebut memiliki interval kredibel 95% yang seluruhnya berada di atas 1, sehingga dapat ditetapkan sebagai wilayah prioritas utama pengendalian kusta.

Sebaliknya, Kabupaten Bandung, Kota Bandung, dan Kota Cimahi memiliki risiko relatif paling rendah. Secara umum, risiko kusta tinggi cenderung terkonsentrasi di wilayah utara dan timur Jawa Barat, sehingga perlu diprioritaskan untuk penemuan kasus aktif, pemeriksaan kontak, serta penguatan layanan diagnosis dan pengobatan.

4 Kesimpulan

Berdasarkan analisis yang telah dilakukan dapat disimpulkan bahwa: 1. Pola insidensi kusta di kabupaten/kota Jawa Barat selama tahun 2019–2024 mengalami fluktuasi. Insidensi tertinggi terjadi pada tahun 2019 sebesar 4,26 kasus per 100.000 penduduk, kemudian menurun pada tahun 2020 dan 2021, meningkat kembali pada tahun 2022–2023, serta menurun menjadi 2,39 kasus per 100.000 penduduk pada tahun 2024. Secara umum, terdapat kecenderungan penurunan insidensi dibandingkan tahun awal penelitian, meskipun perubahan antar tahun tidak berlangsung secara konsisten.

Persebaran spasial insidensi dan Standardized Incidence Ratio (SIR) kusta menunjukkan bahwa kejadian kusta tidak tersebar merata di Jawa Barat. Wilayah dengan insidensi dan SIR tinggi cenderung terkonsentrasi pada bagian utara dan timur Jawa Barat, sedangkan sebagian wilayah tengah dan selatan cenderung memiliki nilai yang lebih rendah. Pola tersebut terlihat berulang pada beberapa tahun pengamatan, sehingga menunjukkan adanya heterogenitas spasial dalam kejadian kusta.
Model Bayesian Spasio-Temporal Poisson BYM2 berhasil mengestimasi risiko relatif kusta pada setiap kabupaten/kota dan tahun pengamatan dengan mempertimbangkan perbedaan jumlah penduduk, hubungan ketetanggaan antarwilayah, perubahan risiko antar tahun, serta interaksi spasio-temporal. Hasil model menunjukkan bahwa risiko relatif kusta memiliki pola spasial dan temporal yang tidak merata, dengan beberapa wilayah secara konsisten memiliki risiko lebih tinggi dibandingkan rata-rata Jawa Barat.
Sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, dan rasio tenaga kesehatan belum menunjukkan pengaruh yang meyakinkan terhadap risiko relatif kusta. Hal tersebut ditunjukkan oleh seluruh interval kredibel 95% relative risk yang masih mencakup nilai 1. Dengan demikian, variasi risiko kusta pada penelitian ini lebih banyak dijelaskan oleh komponen spasial, temporal, dan interaksi spasio-temporal dibandingkan oleh kovariat yang digunakan.
Pada tahun 2024, Kabupaten Indramayu merupakan wilayah dengan risiko relatif kusta tertinggi, yaitu 3,697 kali dibandingkan rata-rata Jawa Barat. Wilayah prioritas berikutnya adalah Kabupaten Bekasi, Kabupaten Subang, Kabupaten Cirebon, Kabupaten Karawang, Kabupaten Kuningan, Kota Bekasi, dan Kabupaten Bogor. Wilayah-wilayah tersebut perlu diprioritaskan dalam penguatan surveilans, penemuan kasus aktif, pemeriksaan kontak serumah, edukasi masyarakat, serta peningkatan akses diagnosis dan pengobatan kusta.

5 Keterbatasan Penelitian

Berikut merupakan keterbatasan penelitian dan rekomendasi yang dapat diberikan berdasarkan hasil analisis risiko relatif kusta di Jawa Barat tahun 2019–2024. Keterbatasan perlu diperhatikan dalam menafsirkan hasil penelitian, sedangkan rekomendasi diharapkan dapat menjadi masukan bagi penelitian selanjutnya maupun program pengendalian kusta.

Data penelitian menggunakan data agregat tingkat kabupaten/kota, sehingga hasil tidak dapat menggambarkan risiko kusta pada tingkat individu.
Periode data hanya mencakup tahun 2019–2024, sehingga pola risiko kusta dalam jangka waktu yang lebih panjang belum dapat dianalisis.
Kovariat yang digunakan terbatas pada sanitasi layak, akses air minum layak, rasio fasilitas kesehatan, dan rasio tenaga kesehatan.
Faktor lain yang berpotensi memengaruhi risiko kusta belum dimasukkan, seperti kepadatan penduduk, kemiskinan, kondisi perumahan, mobilitas penduduk, cakupan pemeriksaan kontak, penemuan kasus, dan kepatuhan pengobatan.
Jumlah kasus tercatat dapat dipengaruhi oleh kualitas surveilans, kemampuan deteksi kasus, akses layanan kesehatan, serta pencatatan dan pelaporan pada masing-masing wilayah.
Penelitian hanya menggunakan satu spesifikasi model BYM2, sehingga belum dilakukan perbandingan kecocokan dengan model spasial atau temporal alternatif.

6 Rekomendasi

Penelitian selanjutnya dapat menggunakan periode data yang lebih panjang agar perubahan risiko kusta antarwaktu dapat diamati dengan lebih baik.
Analisis dapat dilakukan pada tingkat wilayah yang lebih rinci, seperti kecamatan atau desa, agar identifikasi wilayah berisiko tinggi lebih spesifik.
Penelitian berikutnya disarankan menambahkan kovariat seperti kepadatan penduduk, tingkat kemiskinan, kondisi rumah, cakupan pemeriksaan kontak, angka penemuan kasus, proporsi kasus multibasiler, dan aksesibilitas layanan kesehatan.
Perlu dilakukan peningkatan kualitas surveilans, penemuan kasus aktif, serta pencatatan dan pelaporan kasus kusta di seluruh kabupaten/kota.
Kabupaten Indramayu, Kabupaten Bekasi, Kabupaten Subang, Kabupaten Cirebon, dan Kabupaten Karawang perlu diprioritaskan dalam penemuan kasus aktif, pemeriksaan kontak serumah, edukasi masyarakat, serta penguatan layanan diagnosis dan pengobatan.
Peta risiko relatif hasil model BYM2 dapat digunakan sebagai dasar awal dalam penentuan prioritas program, alokasi sumber daya, dan pemantauan pengendalian kusta di Jawa Barat.
Penelitian selanjutnya dapat membandingkan model BYM2 dengan model lain, seperti model Poisson tanpa efek spasial, model temporal, atau model spasial lainnya, untuk memperoleh model dengan kecocokan terbaik.

7 DAFTAR PUSTAKA

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

Badan Pusat Statistik Provinsi Jawa Barat. (2025). Kasus penyakit menurut kabupaten/kota dan jenis penyakit di Provinsi Jawa Barat, 2024. Badan Pusat Statistik Provinsi Jawa Barat.

Dinas Kesehatan Provinsi Jawa Barat. (2025). Jumlah kasus penderita kusta berdasarkan kabupaten/kota di Jawa Barat. Open Data Jabar.

Pemerintah Provinsi Jawa Barat. (2025). Open Data Jabar: Portal data terbuka Provinsi Jawa Barat. Open Data Jabar.

Martino, S., & Riebler, A. (2019). Integrated nested Laplace approximations (INLA). arXiv. https://arxiv.org/abs/1907.01248

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.

Sebong, P. H., et al. (2025). Participatory development of Indonesia’s national action plan for leprosy toward zero leprosy 2030. Frontiers in Public Health.

World Health Organization. (2026, January 23). Leprosy (Hansen disease). WHO leprosy fact

Pemetaan dan Pemodelan Risiko Relatif Kusta Menggunakan Model Bayesian Spasio-Temporal BYM2 di Jawa Barat Tahun 2019–2024

Linda Apriliana

2026-06-13

1 PENDAHULUAN

1.1 Latar Belakang

1.2 Identifikasi Masalah

1.3 Rumusan Masalah

1.4 Tujuan Penelitian

2 Data dan Metode Penelitian

2.1 Data dan Variabel Penelitian

2.2 Metode Analisis Data

2.2.1 Persiapan dan Pembersihan Data

2.2.2 Analisis Deskriptif

2.2.3 Perhitungan Insidensi Kusta

2.2.4 Perhitungan Expected Cases

2.2.5 Standardized Incidence Ratio (SIR)

2.2.6 Analisis Spasial Eksploratif

2.2.7 Pembentukan Matriks Ketetanggaan Spasial

2.2.8 7. Pemodelan Bayesian Spasio-Temporal Poisson BYM2

2.2.9 8. Evaluasi Model

2.2.10 9. Pemetaan Risiko Relatif dan Identifikasi Wilayah Prioritas

3 Hasil dan Pembahasan

3.1 1. Persiapan Data

3.2 Statistik Deskriptif

3.3 Perhitungan Ukuran Epidemiologi

3.3.1 Ukuran kejadian kusta

3.3.2 2. Kasus harapan dan SIR

3.3.3 3. Ketersediaan layanan kesehatan

3.4 Tren Insidensi Kusta di Jawa Barat

3.5 Analisis Spasial Eksploratif

3.5.1 Peta Insidensi Kusta Tahun 2019–2024

3.6 Peta SIR

3.7 Matrix Bobot Spasial

3.8 Hasil Pemodelan Bayessian Spatio-Temporal Poisson (BYM2)

3.8.1 Model Estimasi Bayesian Spasio-Temporal Poisson BYM2

3.9 Evaluasi Kecocokan Model

3.10 Peta Risiko Relatif

3.11 Wilayah Prioritas

4 Kesimpulan

5 Keterbatasan Penelitian

6 Rekomendasi

7 DAFTAR PUSTAKA