Abstrak

Bayi Berat Lahir Rendah (BBLR) merupakan salah satu indikator penting kesehatan ibu dan anak yang masih menjadi masalah kesehatan masyarakat di Indonesia. Penelitian ini bertujuan menganalisis pola persebaran kejadian BBLR antar kabupaten/kota di Provinsi Jawa Barat menggunakan pendekatan disease mapping dan analisis spasio-temporal. Data yang digunakan merupakan data sekunder 27 kabupaten/kota di Jawa Barat periode 2019–2023 yang diperoleh dari Open Data Jabar, BPS Provinsi Jawa Barat, dan data batas administrasi wilayah. Analisis dilakukan melalui perhitungan Incidence Rate (IR), Standardized Incidence Ratio (SIR), Empirical Bayes smoothing, uji autokorelasi spasial Moran’s I, serta pemodelan spasio-temporal Bayesian menggunakan metode Integrated Nested Laplace Approximation (INLA) dengan komponen spasial BYM2 dan temporal Random Walk Order 1 (RW1).

Hasil analisis menunjukkan bahwa kejadian BBLR memiliki pola persebaran yang tidak merata antarwilayah dan terdapat autokorelasi spasial positif pada sebagian besar tahun pengamatan. Peta risiko menunjukkan beberapa wilayah secara konsisten memiliki risiko lebih tinggi dibandingkan wilayah lainnya, terutama Kabupaten Tasikmalaya dan Kabupaten Kuningan. Analisis spasio-temporal juga menunjukkan adanya variasi kejadian BBLR menurut dimensi ruang dan waktu. Temuan ini menunjukkan bahwa pendekatan spasial dan spasio-temporal dapat mendukung identifikasi wilayah prioritas sehingga perencanaan program kesehatan ibu dan anak dapat dilakukan secara lebih tepat sasaran.

Kata kunci: BBLR, disease mapping, spasio-temporal, INLA, Jawa Barat.

1. Pendahuluan

1.1. Latar Belakang

Berat Badan Lahir Rendah (BBLR) merupakan kondisi ketika bayi lahir dengan berat badan kurang dari 2.500 gram tanpa memandang usia kehamilan (World Health Organization, 2023). BBLR merupakan salah satu indikator penting kesehatan ibu dan anak karena berkaitan dengan peningkatan risiko morbiditas dan mortalitas neonatal (World Health Organization, 2023). Bayi yang lahir dengan BBLR juga memiliki risiko lebih tinggi mengalami gangguan pertumbuhan dan perkembangan pada masa anak-anak hingga dewasa (UNICEF & WHO, 2023).

Secara global, sekitar 15–20% bayi lahir dengan kondisi BBLR atau setara dengan lebih dari 20 juta kelahiran setiap tahunnya (UNICEF & WHO, 2023). BBLR masih menjadi salah satu penyebab utama kematian neonatal di berbagai negara (World Health Organization, 2023). Selain itu, BBLR berhubungan dengan meningkatnya risiko stunting, gangguan kognitif, dan penyakit kronis pada masa mendatang (Blencowe dkk., 2019).

Di Indonesia, prevalensi BBLR masih berada pada kisaran 6–7% dari total kelahiran hidup (Kementerian Kesehatan RI, 2023). BBLR juga menjadi salah satu penyumbang utama kematian neonatal di Indonesia (Kementerian Kesehatan RI, 2023). Faktor-faktor yang memengaruhi kejadian BBLR di Indonesia tidak hanya berasal dari kondisi biologis ibu, tetapi juga dipengaruhi oleh faktor sosial ekonomi, pendidikan, dan akses pelayanan kesehatan (Aini & Kurniawan, 2023).

Provinsi Jawa Barat sebagai provinsi dengan jumlah penduduk terbesar di Indonesia masih menghadapi permasalahan BBLR yang cukup tinggi (BPS Provinsi Jawa Barat, 2023). Jumlah kasus BBLR di Jawa Barat pada tahun 2022 tercatat sekitar 19.971 kasus (BPS Provinsi Jawa Barat, 2023). Tingginya jumlah kasus tersebut menunjukkan bahwa upaya pencegahan dan pengendalian BBLR masih perlu menjadi perhatian dalam pembangunan kesehatan daerah (Dinas Kesehatan Provinsi Jawa Barat, 2023).

Kejadian BBLR dipengaruhi oleh berbagai faktor yang bersifat multidimensional (Tshotetsi dkk., 2019). Status gizi ibu selama kehamilan berpengaruh terhadap risiko terjadinya BBLR (Tshotetsi et al., 2019). Cakupan pemberian Tablet Tambah Darah (TTD) berperan dalam pencegahan anemia pada ibu hamil yang merupakan salah satu faktor risiko BBLR (Kementerian Kesehatan RI, 2023). Ketersediaan tenaga kesehatan, kondisi sanitasi, akses air minum layak, tingkat pendidikan, dan kondisi sosial ekonomi juga dilaporkan berkaitan dengan kejadian BBLR (Aini & Kurniawan, 2023). Selain dipengaruhi oleh berbagai faktor risiko, kejadian BBLR juga berpotensi menunjukkan pola spasial antarwilayah. Wilayah yang berdekatan cenderung memiliki karakteristik sosial, ekonomi, dan lingkungan yang serupa sehingga risiko kesehatan dapat membentuk pola pengelompokan geografis (Waller & Gotway, 2004). Oleh karena itu, pendekatan spasial diperlukan untuk mengidentifikasi wilayah dengan risiko BBLR yang relatif tinggi maupun rendah.

Data jumlah kasus BBLR merupakan data cacah (count data) yang dapat dimodelkan menggunakan regresi Poisson atau Negative Binomial (Cameron & Trivedi, 2013). Pada data kesehatan masyarakat sering ditemukan overdispersi sehingga model Negative Binomial menjadi alternatif yang lebih sesuai dibandingkan model Poisson standar (Hilbe, 2011). Selain itu, keberadaan variasi spasial dan temporal memerlukan pendekatan pemodelan yang mampu mengakomodasi ketergantungan antarwilayah dan antarwaktu secara simultan (Knorr-Held, 2000).

Berdasarkan karakteristik tersebut, penelitian ini menggunakan pendekatan Bayesian melalui metode Integrated Nested Laplace Approximation (INLA) dengan komponen spasial Besag–York–Mollié 2 (BYM2) dan efek temporal Random Walk Order 1 (RW1) (Rue dkk., 2009; Riebler dkk., 2016). Pendekatan ini banyak digunakan dalam analisis disease mapping karena mampu menghasilkan estimasi risiko yang lebih stabil dan akurat pada data penyakit yang memiliki variasi spasial dan temporal (Lawson, 2018).

1.2. Rumusan Masalah

Bagaimana pola persebaran kejadian BBLR antar kabupaten/kota di Provinsi Jawa Barat berdasarkan pendekatan disease mapping?
Bagaimana pola kejadian BBLR berdasarkan pendekatan analisis spasio-temporal di Provinsi Jawa Barat?

1.3. Tujuan Penelitian

Menganalisis pola persebaran kejadian Bayi Berat Lahir Rendah (BBLR) antar kabupaten/kota di Provinsi Jawa Barat menggunakan pendekatan disease mapping.
Menganalisis pola kejadian Bayi Berat Lahir Rendah (BBLR) di Provinsi Jawa Barat berdasarkan dimensi ruang dan waktu menggunakan pendekatan analisis spasio-temporal.

2. Data dan Metode

2.1. Sumber Data dan Variabel Penelitian

Penelitian ini menggunakan data sekunder yang berasal dari Open Data Jabar, Badan Pusat Statistik (BPS) Provinsi Jawa Barat, dan data spasial batas administrasi kabupaten/kota Provinsi Jawa Barat dalam bentuk shapefile. Unit analisis penelitian adalah 27 kabupaten/kota di Provinsi Jawa Barat selama periode tahun 2019–2023 sehingga membentuk data panel spasial-temporal.

Data kesehatan yang digunakan meliputi jumlah kasus Bayi Berat Lahir Rendah (BBLR), jumlah bayi lahir, persentase penerima Tablet Tambah Darah (TTD), jumlah bidan puskesmas, dan jumlah ibu hamil yang memperoleh tablet zat besi. Selain itu digunakan pula variabel lingkungan dan sosial ekonomi berupa persentase sanitasi layak, akses air minum layak, angka melek aksara, dan tingkat pengangguran terbuka. Data spasial digunakan untuk membentuk hubungan ketetanggaan antarwilayah dalam analisis spasial dan spasio-temporal.

Tabel 1. Deskripsi Variabel

Variabel	Deskripsi Variabel	Satuan	Sumber Data
Jumlah BBLR	Jumlah bayi yang lahir dengan berat badan kurang dari 2.500 gram pada setiap kabupaten/kota di Provinsi Jawa Barat. Variabel ini digunakan sebagai variabel respon dalam penelitian.	Orang	Open Data Jabar
Persentase Penerima TTD	Persentase ibu hamil yang menerima Tablet Tambah Darah (TTD) sesuai standar pelayanan kesehatan. Variabel ini menggambarkan cakupan intervensi pencegahan anemia pada ibu hamil.	Persen (%)	Open Data Jabar
Persentase Keluarga dengan Akses Sanitasi Layak	Persentase keluarga yang memiliki akses terhadap fasilitas sanitasi layak (jamban sehat) sesuai standar kesehatan lingkungan. Variabel ini menggambarkan kondisi sanitasi masyarakat.	Persen (%)	Open Data Jabar
Jumlah Bidan Puskesmas	Jumlah tenaga kebidanan yang bertugas pada fasilitas puskesmas di suatu wilayah. Variabel ini menggambarkan ketersediaan tenaga kesehatan maternal dan neonatal.	Orang	Open Data Jabar
Persentase Rumah Tangga yang Memiliki Akses terhadap Sumber Air Minum Layak	Persentase rumah tangga yang memiliki akses terhadap sumber air minum layak dan aman untuk dikonsumsi. Variabel ini menggambarkan kondisi kesehatan lingkungan masyarakat.	Persen (%)	Open Data Jabar
Angka Melek Aksara	Persentase penduduk usia 15 tahun ke atas yang mampu membaca dan menulis. Variabel ini menggambarkan tingkat pendidikan masyarakat.	Persen (%)	Open Data Jabar dan BPS Jawa Barat
Tingkat Pengangguran Terbuka	Persentase angkatan kerja yang tidak bekerja dan sedang mencari pekerjaan terhadap total angkatan kerja. Variabel ini menggambarkan kondisi ketenagakerjaan suatu wilayah.	Persen (%)	Open Data Jabar
Jumlah Ibu Hamil Mendapatkan Tablet Zat Besi (FE3-90 Tablet)	Jumlah ibu hamil yang memperoleh minimal 90 tablet zat besi selama masa kehamilan sesuai standar pelayanan kesehatan. Variabel ini menggambarkan cakupan suplementasi zat besi pada ibu hamil.	Orang	Open Data Jabar

Data spasial yang digunakan berupa batas administrasi kabupaten/kota Provinsi Jawa Barat dalam format shapefile yang diperoleh dari Pusat Pemetaan Batas Wilayah (PPBW) Badan Informasi Geospasial (BIG) tahun 2024. Data spasial tersebut digunakan untuk membentuk struktur ketetanggaan antarwilayah, menyusun matriks bobot spasial, serta melakukan visualisasi peta risiko pada analisis disease mapping dan pemodelan spasio-temporal.

2.2. Ukuran Epidemiologi

2.2.1. Incidence Rate

Incidence Rate (IR) untuk menggambarkan tingkat kejadian kasus pada setiap wilayah dan waktu pengamatan. IR dihitung sebagai:

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

dengan:

\(IR_{it}\) = Incidence Rate pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(O_{it}\) = jumlah kasus pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(n_{it}\) = populasi berisiko pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(k\) = konstanta pengali (misalnya 1.000 atau 100.000)

Nilai IR menunjukkan banyaknya kasus yang terjadi per sejumlah populasi berisiko. Semakin besar nilai IR, semakin tinggi tingkat kejadian kasus pada wilayah dan waktu tersebut.

2.2.2. Standardized Incidence Ratio

Standardized Incidence Ratio (SIR) untuk membandingkan jumlah kasus yang diamati dengan jumlah kasus yang diharapkan berdasarkan tingkat kejadian rata-rata. Jumlah kasus yang diharapkan dihitung sebagai:

\[ E_{it} = n_{it} \left( \frac{\sum_{i=1}^{n}O_{it}} {\sum_{i=1}^{n}n_{it}} \right) \]

dengan:

\(E_{it}\) = jumlah kasus yang diharapkan pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(\sum_{i=1}^{n}O_{it}\) = total kasus pada waktu ke-\(t\)
\(\sum_{i=1}^{n}n_{it}\) = total populasi berisiko pada waktu ke-\(t\)

Selanjutnya, nilai SIR dihitung menggunakan rumus:

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

dengan:

\(SIR_{it}\) = Standardized Incidence Ratio pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(O_{it}\) = jumlah kasus yang diamati pada wilayah ke-\(i\) dan waktu ke-\(t\)
\(E_{it}\) = jumlah kasus yang diharapkan pada wilayah ke-\(i\) dan waktu ke-\(t\)

Interpretasi nilai SIR:

\(SIR_{it} > 1\) menunjukkan jumlah kasus yang diamati lebih tinggi daripada jumlah kasus yang diharapkan, sehingga wilayah tersebut memiliki risiko relatif lebih tinggi dibandingkan rata-rata.
\(SIR_{it} = 1\) menunjukkan jumlah kasus yang diamati sama dengan jumlah kasus yang diharapkan.
\(SIR_{it} < 1\) menunjukkan jumlah kasus yang diamati lebih rendah daripada jumlah kasus yang diharapkan, sehingga wilayah tersebut memiliki risiko relatif lebih rendah dibandingkan rata-rata.

2.3. Empirical Bayes Smoothing

Peta risiko berdasarkan Incidence Rate (IR) kasar (crude rate) sering kali menunjukkan variasi yang tinggi, terutama pada wilayah dengan jumlah populasi berisiko yang kecil. Pada wilayah tersebut, perubahan jumlah kasus yang relatif sedikit dapat menghasilkan nilai risiko yang sangat tinggi atau sangat rendah akibat variasi acak (random fluctuation). Kondisi ini dapat menyebabkan interpretasi risiko antarwilayah menjadi kurang stabil (Clayton & Kaldor, 1987).

Untuk mengurangi pengaruh variasi acak tersebut digunakan metode Empirical Bayes (EB) smoothing. Metode ini menghasilkan estimasi risiko yang lebih stabil dengan mengombinasikan risiko yang diamati pada suatu wilayah dengan rata-rata risiko seluruh wilayah melalui proses shrinkage. Wilayah dengan populasi kecil akan mengalami penyesuaian yang lebih besar dibandingkan wilayah dengan populasi besar sehingga estimasi risiko menjadi lebih stabil (Marshall, 1991).

Pada penelitian ini, smoothing dilakukan terhadap Incidence Rate BBLR menggunakan pendekatan Empirical Bayes untuk data Poisson. Misalkan jumlah kasus pada wilayah ke-\(i\) mengikuti distribusi:

\[ O_i \sim Poisson(E_i\theta_i) \]

dengan:

\(O_i\) = jumlah kasus BBLR yang diamati pada wilayah ke-\(i\),
\(E_i\) = jumlah bayi lahir pada wilayah ke-\(i\),
\(\theta_i\) = risiko relatif wilayah ke-\(i\).

Risiko kasar dihitung sebagai:

\[ r_i= \frac{O_i}{E_i} \]

Estimasi Empirical Bayes diperoleh dengan mengombinasikan risiko kasar dan rata-rata risiko seluruh wilayah:

\[ r_i^{EB}= w_i r_i + (1-w_i)\bar r \]

dengan:

\(r_i^{EB}\) = risiko hasil Empirical Bayes,
\(r_i\) = risiko kasar wilayah ke-\(i\),
\(\bar r\) = rata-rata risiko seluruh wilayah,
\(w_i\) = bobot yang bergantung pada ukuran populasi wilayah.

Wilayah dengan jumlah bayi lahir yang besar akan memiliki nilai \(w_i\) yang mendekati satu sehingga estimasi risiko hasil smoothing relatif dekat dengan nilai risiko yang diamati. Sebaliknya, wilayah dengan jumlah bayi lahir yang kecil akan memiliki nilai \(w_i\) yang lebih kecil sehingga estimasi risiko akan lebih mendekati rata-rata keseluruhan wilayah. Dengan demikian, metode Empirical Bayes mampu menghasilkan peta risiko yang lebih stabil dan mengurangi pengaruh fluktuasi acak pada wilayah dengan populasi kecil.

2.4. Autokorelasi Spasial

Hubungan spasial antarwilayah dibentuk menggunakan matriks bobot spasial berbasis queen contiguity. Pada pendekatan ini dua wilayah dianggap bertetangga apabila memiliki sisi maupun titik sudut yang saling bersinggungan. Autokorelasi spasial global diuji menggunakan statistik Moran’s I dengan rumus:

\[ I= \frac{n} {\sum_i\sum_jw_{ij}} \frac{ \sum_i\sum_j w_{ij} (x_i-\bar{x}) (x_j-\bar{x}) } { \sum_i (x_i-\bar{x})^2 } \]

dengan:

\(n\) = jumlah wilayah,
\(w_{ij}\) = bobot spasial antara wilayah ke-\(i\) dan ke-\(j\),
\(x_i\) = nilai variabel pada wilayah ke-\(i\),
\(\bar{x}\) = rata-rata seluruh wilayah.

Interpretasi Moran’s I adalah:

\(I>0\) menunjukkan autokorelasi spasial positif.
\(I<0\) menunjukkan autokorelasi spasial negatif.
\(I \approx 0\) menunjukkan pola acak.

2.5. Regresi Poisson

Regresi Poisson merupakan salah satu metode regresi yang digunakan untuk memodelkan variabel respon berbentuk data cacah (count data) yang menyatakan jumlah kejadian suatu peristiwa dalam periode atau wilayah tertentu. Model ini banyak digunakan dalam bidang epidemiologi untuk menganalisis jumlah kasus penyakit, kematian, maupun kejadian kesehatan lainnya (Cameron & Trivedi, 2013). Misalkan variabel respon \(Y_i\) menyatakan jumlah kejadian pada unit pengamatan ke-\(i\), maka model Poisson mengasumsikan bahwa:

\[ Y_i \sim Poisson(\mu_i) \]

dengan fungsi peluang:

\[ P(Y_i=y_i)= \frac{e^{-\mu_i}\mu_i^{y_i}} {y_i!}, \qquad y_i=0,1,2,\ldots \]

di mana:

\(Y_i\) = jumlah kejadian pada unit ke-\(i\),
\(\mu_i\) = nilai harapan (mean) jumlah kejadian,
\(y_i\) = nilai observasi.

Hubungan antara rata-rata kejadian dan variabel prediktor dimodelkan menggunakan fungsi hubungan logaritmik (log link function):

\[ \log(\mu_i) = \beta_0 + \beta_1X_{i1} + \beta_2X_{i2} + \cdots + \beta_pX_{ip} \]

atau dapat dituliskan sebagai:

\[ \mu_i = \exp \left(\beta_0 + \sum_{k=1}^{p} \beta_kX_{ik} \right) \]

dengan:

\(\beta_0\) = intersep,
\(\beta_k\) = koefisien regresi variabel ke-\(k\),
\(X_{ik}\) = nilai variabel prediktor ke-\(k\) pada unit ke-\(i\).

Salah satu karakteristik utama distribusi Poisson adalah asumsi equidispersion, yaitu nilai rata-rata dan varians memiliki nilai yang sama:

\[ E(Y_i)=Var(Y_i)=\mu_i \]

Apabila varians data lebih besar daripada rata-ratanya, maka terjadi kondisi overdispersi yang menyebabkan model Poisson menjadi kurang sesuai digunakan (Hilbe, 2011).

2.6. Pemeriksaan Overdispersi

Sebelum membangun model spasio-temporal dilakukan pemeriksaan overdispersi untuk mengevaluasi kesesuaian asumsi regresi Poisson. Overdispersi diperiksa menggunakan rasio devians terhadap derajat bebas:

\[ \phi= \frac{\text{Deviance}} {\text{df}} \]

dengan:

\(\phi\) = parameter dispersi,
df = derajat bebas residual.

Interpretasi nilai dispersi adalah:

\(\phi \approx 1\) menunjukkan asumsi equidispersion terpenuhi.
\(\phi > 1\) menunjukkan adanya overdispersi.
\(\phi < 1\) menunjukkan adanya underdispersi.

Menurut Hilbe (2011), overdispersi menyebabkan standar error pada model Poisson menjadi terlalu kecil sehingga hasil pengujian parameter dapat menjadi kurang akurat. Oleh karena itu, apabila ditemukan overdispersi maka digunakan model Negative Binomial yang memiliki parameter dispersi tambahan. Pada model Negative Binomial, varians dinyatakan sebagai:

\[ Var(Y_i) = \mu_i + \frac{\mu_i^2}{\theta} \]

dengan:

\(\mu_i\) = rata-rata,
\(\theta\) = parameter dispersi.

Model ini memungkinkan varians lebih besar daripada rata-rata sehingga lebih sesuai untuk data cacah yang mengalami overdispersi.

2.7. Model Spasio-Temporal Bayesian INLA

Analisis spasio-temporal dilakukan menggunakan pendekatan Bayesian melalui metode Integrated Nested Laplace Approximation (INLA) yang diperkenalkan oleh Rue, Martino, dan Chopin (2009). Metode INLA digunakan untuk memperoleh estimasi posterior secara efisien pada model hierarkis Bayesian tanpa memerlukan simulasi Markov Chain Monte Carlo (MCMC). Pada penelitian ini digunakan model Negative Binomial dengan komponen efek spasial dan temporal. Secara umum model dapat dituliskan sebagai:

\[ Y_{it} \sim NegBin(\mu_{it},\theta) \]

dengan:

\[ \log(\mu_{it}) = \log(E_{it}) + \beta_0 + \sum_{k=1}^{p} \beta_kX_{kit} + u_i + v_i + \gamma_t \]

dengan:

\(Y_{it}\) = jumlah kasus BBLR pada wilayah ke-\(i\) dan tahun ke-\(t\),
\(E_{it}\) = jumlah kasus harapan (expected cases),
\(X_{kit}\) = variabel prediktor ke-\(k\),
\(\beta_k\) = koefisien regresi,
\(u_i\) = efek spasial terstruktur,
\(v_i\) = efek spasial tidak terstruktur,
\(\gamma_t\) = efek temporal.

Komponen \(\log(E_{it})\) digunakan sebagai offset sehingga model yang dihasilkan mengukur risiko relatif terhadap jumlah kasus yang diharapkan.

2.7.1. Efek Spasial BYM2

Efek spasial dimodelkan menggunakan pendekatan Besag–York–Mollié 2 (BYM2) yang dikembangkan oleh Riebler dkk. (2016). Model BYM2 menggabungkan dua komponen efek acak, yaitu efek spasial terstruktur yang menggambarkan kemiripan risiko antarwilayah bertetangga dan efek spasial tidak terstruktur yang menangkap variasi acak yang tidak dijelaskan oleh lokasi geografis. Hubungan ketetanggaan antarwilayah dibentuk menggunakan matriks ketetanggaan berbasis queen contiguity, yaitu dua wilayah dianggap bertetangga apabila memiliki sisi atau titik sudut yang saling bersinggungan.

2.7.2. Efek Temporal Random Walk Order 1 (RW1)

Efek temporal dimodelkan menggunakan Random Walk Order 1 (RW1) yang mengasumsikan bahwa kondisi pada suatu tahun dipengaruhi oleh kondisi pada tahun sebelumnya.

Secara matematis:

\[ \gamma_t = \gamma_{t-1} + \varepsilon_t \]

dengan:

\[ \varepsilon_t \sim N(0,\sigma^2) \]

Pendekatan ini menghasilkan pola perubahan risiko yang lebih halus dari waktu ke waktu dan banyak digunakan dalam pemodelan spasio-temporal penyakit (Knorr-Held, 2000).

2.8. Pemilihan Model

Pemilihan model dilakukan menggunakan nilai Deviance Information Criterion (DIC) dan Widely Applicable Information Criterion (WAIC).

2.8.1. Deviance Information Criterion (DIC)

DIC merupakan ukuran pemilihan model Bayesian yang mempertimbangkan kecocokan model dan kompleksitas parameter (Spiegelhalter dkk., 2002).

\[ DIC= \bar D + p_D \]

dengan:

\(\bar D\) = rata-rata deviance posterior,
\(p_D\) = jumlah parameter efektif.

Model dengan nilai DIC yang lebih kecil dianggap memiliki kecocokan yang lebih baik.

2.8.2. Widely Applicable Information Criterion (WAIC)

WAIC merupakan ukuran pemilihan model Bayesian yang digunakan untuk mengevaluasi kemampuan prediksi model terhadap data baru (Watanabe, 2010).

\[ WAIC = -2 \left(lppd - p_{WAIC}\right) \]

dengan:

\(lppd\) = log pointwise predictive density,
\(p_{WAIC}\) = jumlah parameter efektif.

Model dengan nilai WAIC yang lebih kecil dianggap memiliki kemampuan prediksi yang lebih baik. Pada penelitian ini, model terbaik dipilih berdasarkan nilai DIC dan WAIC terkecil.

3. Hasil dan Pembahasan

3.1. Cek Kualitas Data

Persiapan dan Cleaning Data

setwd("G:/My Drive/Semester 2/Epidemiologi/Tugas Individu")

# =========================================================
# IMPORT DATA
# =========================================================

data_awal = import("Data Lengkap.xlsx")

jabar_shp <- st_read("Batas_Kabupaten_Jawa_Barat.shp")

# =========================================================
# PENYESUAIAN NAMA WILAYAH
# =========================================================

data_awal <- data_awal %>%
  mutate(
    WADMKK = nama_kabupaten_kota %>%
      str_to_title() %>%
      str_replace("^Kabupaten ", "") %>%
      str_replace("^Kota ", "Kota ")
  )

Cek Missing Value dan Nilai Nol

# =========================================================
# CEK MISSING VALUE
# =========================================================

tabel_missing <- data.frame(
  Variabel = names(data_awal),
  Missing = colSums(is.na(data_awal))
)

# =========================================================
# CEK NILAI 0
# =========================================================

tabel_nol <- data_awal %>%
  
  filter(jumlah_bblr == 0) %>%
  
  select(
    WADMKK,
    tahun,
    jumlah_bblr
  )

Imputasi Data

# =========================================================
# IMPUTASI DENGAN RATA-RATA HISTORIS
# =========================================================

data_imputasi <- data_awal %>%
  
  arrange(WADMKK, tahun) %>%
  
  group_by(WADMKK) %>%
  
  mutate(
    
    # ubah 0 jadi NA
    jumlah_bblr = ifelse(
      jumlah_bblr == 0,
      NA,
      jumlah_bblr
    ),
    
    # imputasi rata-rata historis wilayah
    jumlah_bblr_imputasi = ifelse(
      
      is.na(jumlah_bblr),
      
      round(
        mean(
          jumlah_bblr,
          na.rm = TRUE
        )
      ),
      
      jumlah_bblr
    )
  ) %>%
  
  ungroup()

# =========================================================
# TABEL HASIL IMPUTASI
# =========================================================

tabel_imputasi <- data_imputasi %>%

  filter(is.na(jumlah_bblr)) %>%

  select(
    WADMKK,
    tahun,
    jumlah_bblr,
    jumlah_bblr_imputasi
  ) %>%

  arrange(WADMKK, tahun)

kable(
  tabel_imputasi,
  caption = "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 2. Data Asli dan Hasil Imputasi
  </div>",
  escape = FALSE
) %>%

  kable_styling(
    full_width = FALSE,
    position = "center"
  )

Tabel 2. Data Asli dan Hasil Imputasi
WADMKK	tahun	jumlah_bblr	jumlah_bblr_imputasi
Ciamis	2023	NA	821
Kota Bandung	2023	NA	778
Kota Cirebon	2023	NA	237
Subang	2023	NA	549

Pada proses pemeriksaan data tidak ditemukan missing value namun ditemukan nilai nol pada variabel jumlah BBLR tahun 2023 di Kabupaten Ciamis, Kabupaten Subang, Kota Bandung, dan Kota Cirebon. Nilai nol tersebut diperlakukan sebagai data hilang (missing value) atau kemungkinan kesalahan pencatatan karena tidak konsisten dengan pola historis data. Sehingga dilakukan imputasi menggunakan metode rata-rata historis wilayah, yaitu mengganti nilai yang hilang dengan rata-rata jumlah kasus BBLR tahun 2019–2022 pada wilayah yang sama. Hasil imputasi menghasilkan estimasi jumlah kasus BBLR sebesar 821 kasus di Kabupaten Ciamis, 549 kasus di Kabupaten Subang, 778 kasus di Kota Bandung, dan 237 kasus di Kota Cirebon.

3.2. Statistik Desktriptif dan Exploratory Data Analysis

3.2.1. Statistik Deskriptif

tabel_ringkas <- data_imputasi %>%

group_by(tahun) %>%

summarise(

Total = sum(
  jumlah_bblr_imputasi,
  na.rm = TRUE
),

`Rata-rata` = round(
  mean(
    jumlah_bblr_imputasi,
    na.rm = TRUE
  ),
  2
),

Median = round(
  median(
    jumlah_bblr_imputasi,
    na.rm = TRUE
  ),
  2
),

Minimum = min(
  jumlah_bblr_imputasi,
  na.rm = TRUE
),

Maksimum = max(
  jumlah_bblr_imputasi,
  na.rm = TRUE
),

`Standar Deviasi` = round(
  sd(
    jumlah_bblr_imputasi,
    na.rm = TRUE
  ),
  2
),

`Standard Error` = round(
  sd(
    jumlah_bblr_imputasi,
    na.rm = TRUE
  ) /
  sqrt(
    sum(
      !is.na(jumlah_bblr_imputasi)
    )
  ),
  2
),

)

kable(
tabel_ringkas,

caption = "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 3. Ringkasan Statistik Kasus BBLR per Tahun
  </div>",

align = "c"
) %>%

kable_styling(
full_width = FALSE,
position = "center"
)

Tabel 3. Ringkasan Statistik Kasus BBLR per Tahun
tahun	Total	Rata-rata	Median	Minimum	Maksimum	Standar Deviasi	Standard Error
2019	21744	805.33	723	129	2033	549.59	105.77
2020	21568	798.81	723	122	2033	538.98	103.73
2021	20861	772.63	548	136	2015	540.75	104.07
2022	21130	782.59	637	118	1965	548.43	105.55
2023	24330	901.11	778	132	2672	680.84	131.03

Berdasarkan tabel ringkasan statistik setelah imputasi, jumlah kasus BBLR di Jawa Barat selama periode 2019–2023 menunjukkan pola yang berfluktuasi. Total kasus tertinggi terjadi pada tahun 2023 yaitu sebanyak 24.330 kasus dengan rata-rata 901,11 kasus per kabupaten/kota, sedangkan total terendah terjadi pada tahun 2021 yaitu sebanyak 20.861 kasus. Nilai maksimum kasus tertinggi juga terjadi pada tahun 2023 yaitu sebesar 2.672 kasus.

Nilai standar deviasi yang cukup besar pada setiap tahun menunjukkan adanya variasi jumlah kasus BBLR antar kabupaten/kota di Jawa Barat. Standar deviasi tertinggi terdapat pada tahun 2023 yaitu sebesar 680,84 yang mengindikasikan perbedaan jumlah kasus antarwilayah relatif tinggi.

3.2.2. Grafik Tren Kasus BBLR

# ====================================
# GRAFIK TREN KASUS BBLR
# ====================================

tabel_tren <- data_imputasi %>%
  
  group_by(tahun) %>%
  
  summarise(
    Total = sum(
      jumlah_bblr_imputasi,
      na.rm = TRUE
    )
  )

ggplot(
  tabel_tren,
  aes(x = tahun, y = Total)
) +
  
  geom_area(
    fill = "#F4A6A6",
    alpha = 0.35
  ) +
  
  geom_line(
    color = "#C1121F",
    linewidth = 1.8
  ) +
  
  geom_point(
    color = "#C1121F",
    size = 4
  ) +
  
  geom_text(
    aes(label = comma(Total)),
    vjust = -1,
    fontface = "bold",
    size = 4
  ) +
  
  coord_cartesian(
    ylim = c(20000, 25000)
  ) +
  
  scale_y_continuous(
    labels = comma
  ) +
  
  labs(
    title = "Tren Kasus BBLR di Jawa Barat",
    subtitle = "Periode Tahun 2019–2023",
    x = "Tahun",
    y = "Jumlah Kasus BBLR"
  ) +
  
  theme_minimal(base_size = 14) +
  
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 20
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40",
      size = 12
    ),
    
    axis.title = element_text(
      face = "bold",
      size = 14
    ),
    
    axis.text = element_text(
      size = 12
    ),
    
    panel.grid.minor = element_blank(),
    
    panel.grid.major.x = element_blank()
  )

Gambar 1. Tren Kasus BBLR Jawa Barat

Berdasarkan grafik tren kasus BBLR di Jawa Barat tahun 2019–2023, jumlah kasus BBLR cenderung mengalami penurunan pada periode 2019–2021. Penurunan ini diduga berkaitan dengan kondisi pandemi COVID-19 yang terjadi sejak tahun 2020. Pada masa pandemi, terjadi pembatasan aktivitas masyarakat dan perubahan sistem pelayanan kesehatan sehingga pelaporan maupun akses pemeriksaan ke fasilitas kesehatan dapat mengalami gangguan. Kondisi tersebut berpotensi memengaruhi pencatatan kasus BBLR di beberapa wilayah.

# ====================================
# GRAFIK TREN BBLR PER DAERAH
# ====================================

# Hitung total kasus per daerah
top5_daerah <- data_imputasi %>%
  
  group_by(WADMKK) %>%
  
  summarise(
    Total_BBLR = sum(
      jumlah_bblr_imputasi,
      na.rm = TRUE
    ),
    .groups = "drop"
  ) %>%
  
  arrange(desc(Total_BBLR)) %>%
  
  slice(1:5) %>%
  
  pull(WADMKK)

# Data tren
tabel_tren_daerah <- data_imputasi %>%
  
  group_by(WADMKK, tahun) %>%
  
  summarise(
    Total = sum(
      jumlah_bblr_imputasi,
      na.rm = TRUE
    ),
    .groups = "drop"
  ) %>%
  
  mutate(
    kategori = ifelse(
      WADMKK %in% top5_daerah,
      WADMKK,
      "Lainnya"
    )
  )

# Grafik
ggplot(
  tabel_tren_daerah,
  aes(
    x = tahun,
    y = Total,
    group = WADMKK,
    color = kategori
  )
) +
  
  geom_line(
    aes(
      linewidth = kategori != "Lainnya"
    ),
    alpha = 0.9
  ) +
  
  geom_point(
    aes(
      size = kategori != "Lainnya"
    )
  ) +
  
  scale_x_continuous(
    breaks = 2019:2023
  ) +
  
  scale_linewidth_manual(
    values = c(
      "TRUE" = 1.4,
      "FALSE" = 0.5
    ),
    guide = "none"
  ) +
  
scale_size_manual(
  values = c(
    "TRUE" = 2.5,
    "FALSE" = 1
  ),
  guide = "none"
) +

scale_color_manual(
  values = c(
    "Lainnya" = "grey70",
    setNames(
      scales::hue_pal()(length(top5_daerah)),
      top5_daerah
    )
  )
) +
  
  labs(
    title = "Tren Kasus BBLR per Kabupaten/Kota di Jawa Barat",
    subtitle = "Sorotan pada 5 daerah dengan total kasus tertinggi (2019–2023)",
    x = "Tahun",
    y = "Jumlah Kasus BBLR",
    color = "Kabupaten/Kota"
  ) +
  
  theme_minimal(base_size = 13) +
  
  theme(
  plot.title = element_text(
    hjust = 0.5,
    face = "bold",
    size = 18
  ),
  
  plot.subtitle = element_text(
    hjust = 0.5,
    color = "grey40"
  ),
  
  axis.title = element_text(
    face = "bold"
  ),
  
  legend.position = "right",
  
  legend.title = element_text(
    face = "bold"
  ),
  
  legend.text = element_text(
    size = 8
  ),
  
  panel.grid.minor = element_blank(),
  
  plot.margin = margin(
    t = 20,
    r = 20,
    b = 10,
    l = 20
  )
)

Gambar 2. Tren Kasus BBLR Jawa Barat per Kabupaten/Kota

Grafik menunjukkan bahwa tren kasus BBLR di kabupaten/kota Jawa Barat selama tahun 2019–2023 cenderung berfluktuasi antarwilayah. Kabupaten Bandung mengalami peningkatan paling tajam pada tahun 2023 dan menjadi daerah dengan jumlah kasus tertinggi. Kabupaten Bogor dan Sukabumi juga menunjukkan jumlah kasus yang relatif tinggi dengan kecenderungan meningkat pada akhir periode pengamatan, sedangkan Kabupaten Garut mengalami kenaikan sejak tahun 2021 dan kemudian relatif stabil. Daerah lainnya memiliki jumlah kasus yang lebih rendah dan cenderung stabil dibandingkan lima daerah dengan kasus tertinggi.

3.3. Ukuran Epidemiologi

3.3.1. Perhitungan IR dan SIR

# =========================================================
# JOIN DATA SPASIAL
# =========================================================

# gunakan shapefile sebagai base object
jabar_join <- jabar_shp %>%
  
  left_join(
    data_imputasi,
    by = "WADMKK"
  )

# pastikan object tetap sf
jabar_join <- st_as_sf(jabar_join)

# perbaiki geometry jika ada yang invalid
jabar_join <- st_make_valid(jabar_join)

# =========================================================
# PERHITUNGAN IR DAN SIR
# =========================================================

jabar_join <- jabar_join %>%
  
  group_by(tahun) %>%
  
  mutate(
    
    # =====================================
    # INCIDENCE RATE (per 1.000 kelahiran)
    # =====================================
    
    IR = (
      jumlah_bblr_imputasi /
        jumlah_bayi_lahir
    ) * 1000,
    
    # =====================================
    # EXPECTED CASES
    # berdasarkan rata-rata provinsi
    # pada tahun yang sama
    # =====================================
    
    expected_cases =
      jumlah_bayi_lahir *
      
      (
        sum(
          jumlah_bblr_imputasi,
          na.rm = TRUE
        ) /
        
        sum(
          jumlah_bayi_lahir,
          na.rm = TRUE
        )
      ),
    
    # =====================================
    # STANDARDIZED INCIDENCE RATIO
    # =====================================
    
    SIR =
      jumlah_bblr_imputasi /
      expected_cases
  ) %>%
  
  ungroup()

# =========================================================
# TABEL IR DAN SIR
# =========================================================

tabel_sir <- jabar_join %>%
  
  st_drop_geometry() %>%
  
  select(
    WADMKK,
    tahun,
    jumlah_bblr_imputasi,
    jumlah_bayi_lahir,
    IR,
    expected_cases,
    SIR
  ) %>%
  
  mutate(
    IR = round(IR, 2),
    expected_cases = round(expected_cases, 2),
    SIR = round(SIR, 3)
  )

kable(
  head(tabel_sir, 15),
  caption = "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 4. Contoh Hasil Perhitungan IR dan SIR Kasus BBLR
  </div>"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )

Tabel 4. Contoh Hasil Perhitungan IR dan SIR Kasus BBLR
WADMKK	tahun	jumlah_bblr_imputasi	jumlah_bayi_lahir	IR	expected_cases	SIR
Bogor	2019	1633	116630	14.00	2764.12	0.591
Bogor	2020	1790	117040	15.29	2870.69	0.624
Bogor	2021	1383	139565	9.91	3299.41	0.419
Bogor	2022	1927	108215	17.81	2862.28	0.673
Bogor	2023	2200	107699	20.43	3214.45	0.684
Sukabumi	2019	2033	41442	49.06	982.17	2.070
Sukabumi	2020	2033	47162	43.11	1156.76	1.757
Sukabumi	2021	2015	45196	44.58	1068.46	1.886
Sukabumi	2022	1965	42772	45.94	1131.32	1.737
Sukabumi	2023	2240	42995	52.10	1283.25	1.746
Cianjur	2019	1143	42702	26.77	1012.03	1.129
Cianjur	2020	1143	40408	28.29	991.10	1.153
Cianjur	2021	1085	42438	25.57	1003.26	1.081
Cianjur	2022	1220	41700	29.26	1102.96	1.106
Cianjur	2023	1125	39495	28.48	1178.79	0.954

Pada penelitian ini digunakan dua ukuran epidemiologi, yaitu Incidence Rate (IR) dan Standardized Incidence Ratio (SIR). Incidence Rate digunakan untuk menggambarkan besarnya kejadian kasus BBLR di antara seluruh kelahiran hidup pada setiap kabupaten/kota dan periode waktu tertentu. Secara epidemiologis, denominator yang seharusnya digunakan dalam perhitungan BBLR adalah jumlah bayi lahir hidup (live births). Namun, karena keterbatasan ketersediaan data, penelitian ini menggunakan variabel jumlah bayi baru lahir sebagai pendekatan populasi berisiko. Rumus Incidence Rate adalah:

\[ IR = \frac{\text{Jumlah Kasus BBLR}}{\text{Jumlah Bayi Baru Lahir}} \times 1000 \]

Sebagai contoh, Kabupaten Sukabumi tahun 2019 memiliki nilai IR sebesar 49,06 per 1.000 kelahiran, yang menunjukkan bahwa dari setiap 1.000 bayi lahir terdapat sekitar 49 bayi dengan kondisi BBLR.

Selain itu, digunakan Standardized Incidence Ratio (SIR) untuk membandingkan jumlah kasus BBLR yang diamati pada suatu wilayah dengan jumlah kasus yang diharapkan berdasarkan rata-rata kejadian di seluruh Jawa Barat. Perhitungan SIR diawali dengan menghitung expected cases menggunakan rumus:

\[ E_i = n_i \times \frac{\sum O_i}{\sum n_i} \]

dengan:

\(E_i\) = jumlah kasus BBLR yang diharapkan pada wilayah ke-\(i\),
\(n_i\) = jumlah bayi baru lahir pada wilayah ke-\(i\),
\(\sum O_i\) = total kasus BBLR seluruh kabupaten/kota di Jawa Barat,
\(\sum n_i\) = total jumlah bayi baru lahir seluruh kabupaten/kota di Jawa Barat.

Selanjutnya, nilai SIR dihitung menggunakan rumus:

\[ SIR_i = \frac{O_i}{E_i} \]

dengan:

\(O_i\) = jumlah kasus BBLR yang diamati pada wilayah ke-\(i\),
\(E_i\) = jumlah kasus BBLR yang diharapkan pada wilayah ke-\(i\).

Nilai SIR lebih dari 1 menunjukkan bahwa risiko BBLR pada suatu wilayah lebih tinggi dari kasus yang diharapkan, sedangkan nilai SIR kurang dari 1 menunjukkan risiko yang lebih rendah dibandingkan kasus yang diharapkan.

3.3.2. Visualisasi Incidence Rate per Tahun

# =========================================================
# BOXPLOT IR PER TAHUN
# =========================================================

ggplot(
  jabar_join,
  aes(
    x = factor(tahun),
    y = IR
  )
) +
  
  geom_boxplot(
    fill = "#A8DADC",
    color = "#1D3557",
    alpha = 0.8
  ) +
  
  geom_jitter(
    width = 0.15,
    alpha = 0.5,
    color = "#1D3557"
  ) +
  
  labs(
    title = "Distribusi Incidence Rate (IR) BBLR per Tahun",
    subtitle = "Kabupaten/Kota di Jawa Barat Tahun 2019–2023",
    x = "Tahun",
    y = "IR per 1.000 Bayi Baru Lahir"
  ) +
  
  theme_minimal(base_size = 14) +
  
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 20
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40"
    ),
    
    axis.title = element_text(
      face = "bold"
    ),
    
    panel.grid.minor = element_blank()
  )

Gambar 3. Boxplot IR per Tahun

Boxplot menunjukkan distribusi nilai Incidence Rate (IR) kasus BBLR antar kabupaten/kota pada setiap tahun pengamatan. Semakin tinggi median boxplot menunjukkan semakin tinggi risiko BBLR secara umum pada tahun tersebut, sedangkan titik ekstrem menunjukkan adanya wilayah dengan risiko jauh lebih tinggi dibanding wilayah lainnya. Berdasarkan boxplot distribusi Incidence Rate (IR) kasus BBLR tahun 2019–2023, terlihat bahwa nilai IR antar kabupaten/kota di Jawa Barat cenderung berfluktuasi setiap tahun. Median IR mengalami peningkatan pada tahun 2022–2023 yang menunjukkan adanya peningkatan risiko BBLR secara umum pada sebagian besar wilayah. Selain itu, lebar boxplot yang relatif besar mengindikasikan adanya variasi tingkat kejadian BBLR antarwilayah yang cukup tinggi. Grafik juga menunjukkan adanya beberapa titik ekstrem (outlier) pada setiap tahun pengamatan, terutama pada tahun 2021 dan 2022. Hal tersebut menandakan bahwa terdapat beberapa kabupaten/kota dengan nilai IR jauh lebih tinggi dibandingkan wilayah lainnya.

3.4. Disease Mapping

Disease mapping dilakukan untuk menggambarkan distribusi risiko kasus BBLR antar kabupaten/kota di Jawa Barat menggunakan peta crude risk berdasarkan Incidence Rate (IR) dan metode Empirical Bayes smoothing.

3.4.1. Analisis Autokorelasi Spasial Global Moran’s I Berdasarkan IR

Analisis Moran’s I dilakukan untuk mengetahui pola autokorelasi spasial berdasarkan nilai Incidence Rate (IR) kasus BBLR antar kabupaten/kota di Jawa Barat periode 2019–2023.

# =========================================================
# SPATIAL WEIGHT MATRIX
# =========================================================

# ambil satu tahun
# karena batas wilayah tetap
data_sp <- jabar_join %>%
  filter(tahun == 2019)

# ubah ke spatial object
data_sp <- as(
  data_sp,
  "Spatial"
)

# queen contiguity
nb <- poly2nb(
  data_sp,
  queen = TRUE
)

# spatial weights
lw <- nb2listw(
  nb,
  style = "W",
  zero.policy = TRUE
)

# =========================================================
# MORAN'S I BERDASARKAN IR
# =========================================================

hasil_moran_ir <- data.frame()

for(th in sort(unique(jabar_join$tahun))){
  
  # filter data per tahun
  data_th <- jabar_join %>%
    filter(tahun == th)
  
  # hitung Moran's I
  moran_ir <- moran.test(
    data_th$IR,
    lw,
    zero.policy = TRUE
  )
  
  # simpan hasil
  hasil_moran_ir <- rbind(
    hasil_moran_ir,
    
    data.frame(
      Tahun = th,
      
      Moran_I = as.numeric(
        moran_ir$estimate[1]
      ),
      
      Expected_I = as.numeric(
        moran_ir$estimate[2]
      ),
      
      Variance = as.numeric(
        moran_ir$estimate[3]
      ),
      
      P_Value = moran_ir$p.value
    )
  )
}

# =========================================================
# TABEL HASIL MORAN'S I IR
# =========================================================

hasil_moran_ir_tabel <- hasil_moran_ir %>%
  
  mutate(
    Moran_I = round(Moran_I, 4),
    Expected_I = round(Expected_I, 4),
    Variance = round(Variance, 4),
    P_Value = round(P_Value, 4),
    
    Interpretasi = ifelse(
      P_Value < 0.05 & Moran_I > 0,
      "Autokorelasi Positif",
      
      ifelse(
        P_Value < 0.05 & Moran_I < 0,
        "Autokorelasi Negatif",
        "Tidak Signifikan"
      )
    )
  )

kable(
  hasil_moran_ir_tabel,
  caption = "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 5. Hasil Global Moran’s I Berdasarkan IR Kasus BBLR Tahun 2019–2023
  </div>",
  align = "c"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )

Tabel 5. Hasil Global Moran’s I Berdasarkan IR Kasus BBLR Tahun 2019–2023
Tahun	Moran_I	Expected_I	Variance	P_Value	Interpretasi
2019	0.6581	-0.0385	0.0195	0.0000	Autokorelasi Positif
2020	0.6723	-0.0385	0.0199	0.0000	Autokorelasi Positif
2021	0.1008	-0.0385	0.0110	0.0917	Tidak Signifikan
2022	0.1951	-0.0385	0.0081	0.0047	Autokorelasi Positif
2023	0.5808	-0.0385	0.0198	0.0000	Autokorelasi Positif

Berdasarkan hasil Global Moran’s I, nilai Incidence Rate (IR) kasus BBLR di Jawa Barat menunjukkan adanya autokorelasi spasial positif yang signifikan pada tahun 2019, 2020, 2022, dan 2023. Hal ini ditunjukkan oleh nilai Moran’s I yang positif dengan p-value < 0,05, yang mengindikasikan bahwa kabupaten/kota dengan nilai IR tinggi cenderung berdekatan dengan wilayah yang juga memiliki nilai IR tinggi, begitu pula wilayah dengan nilai IR rendah cenderung berdekatan dengan wilayah bernilai rendah. Pada tahun 2021, tidak menunjukkan autokorelasi spasial yang signifikan karena memiliki p-value > 0,05 sehingga pola persebaran IR kasus BBLR cenderung lebih acak antarwilayah.

# =========================================================
# MORAN SCATTERPLOT IR
# =========================================================

par(mfrow = c(2,3))

for(th in sort(unique(jabar_join$tahun))){
  
  data_th <- jabar_join %>%
    filter(tahun == th)
  
  moran.plot(
    data_th$IR,
    lw,
    zero.policy = TRUE,
    labels = data_th$WADMKK,
    pch = 19,
    col = "darkblue",
    main = paste(
      "Moran Scatterplot IR",
      th
    )
  )
}

Gambar 4. Moran Scatter Plot

Berdasarkan Moran Scatterplot tahun 2019–2023, seluruh tahun menunjukkan kecenderungan autokorelasi spasial positif pada incidence rate (IR) BBLR. Hal ini menandakan bahwa wilayah dengan IR tinggi cenderung berdekatan dengan wilayah dengan IR tinggi, wilayah dengan IR rendah berdekatan dengan wilayah IR rendah.

3.4.2. Peta Crude Risk Berdasarkan IR

# =========================================================
# PETA CRUDE RISK 
# =========================================================


breaks_ir <- quantile(
  jabar_join$IR,
  
  probs = seq(0, 1, 0.2),
  
  na.rm = TRUE
)

# hindari duplikat break
breaks_ir <- unique(breaks_ir)

# =========================================================
# BUAT LABEL RENTANG ANGKA
# =========================================================

label_ir <- paste0(
  
  round(
    head(breaks_ir, -1),
    1
  ),
  
  " - ",
  
  round(
    tail(breaks_ir, -1),
    1
  )
)

# =========================================================
# KATEGORI IR
# =========================================================

jabar_join <- jabar_join %>%
  
  mutate(
    
    kategori_IR = cut(
      IR,
      
      breaks = breaks_ir,
      
      include.lowest = TRUE,
      
      labels = label_ir
    )
  )

# =========================================================
# PLOT PETA
# =========================================================

ggplot(jabar_join) +
  
  geom_sf(
    aes(fill = kategori_IR),
    
    color = "grey95",
    
    linewidth = 0.4
  ) +
  
  facet_wrap(
    ~tahun,
    ncol = 3
  ) +
  
  scale_fill_manual(
    
    values = c(
      "#FFF5EB",
      "#FDBE85",
      "#FD8D3C",
      "#E6550D",
      "#A63603"
    ),
    
    name = "IR per 1.000\nkelahiran"
  ) +
  
  labs(
    title = "Peta Crude Risk BBLR Berdasarkan Incidence Rate (IR)",
    
    subtitle = "Kabupaten/Kota di Jawa Barat Tahun 2019–2023",
    
  ) +
  
  theme_minimal(base_size = 14) +
  
  theme(
    
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 20
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40",
      size = 12
    ),
    
    strip.text = element_text(
      face = "bold",
      size = 12
    ),
    
    legend.title = element_text(
      face = "bold"
    ),
    
    legend.text = element_text(
      size = 10
    ),
    
    axis.text = element_blank(),
    
    axis.title = element_blank(),
    
    panel.grid = element_blank(),
    
    plot.caption = element_text(
      hjust = 0.5,
      color = "grey40"
    )
  )

Gambar 5. Peta Crude Risk

Berdasarkan peta crude risk Incidence Rate (IR) kasus BBLR tahun 2023, terlihat bahwa distribusi risiko BBLR antar kabupaten/kota di Jawa Barat tidak merata. Wilayah dengan warna lebih gelap menunjukkan nilai IR yang lebih tinggi, yang berarti memiliki kejadian kasus BBLR lebih besar dibandingkan wilayah lainnya. Beberapa wilayah di bagian timur dan barat Jawa Barat terlihat memiliki risiko BBLR relatif tinggi, sedangkan wilayah dengan warna lebih terang menunjukkan nilai IR yang lebih rendah. Pola ini mengindikasikan adanya variasi risiko BBLR antarwilayah yang kemungkinan dipengaruhi oleh faktor sosial, ekonomi, maupun akses pelayanan kesehatan.

3.4.3. Empirical Bayes Smoothing

# =========================================================
# EMPIRICAL BAYES SMOOTHING BERDASARKAN IR
# =========================================================

library(spdep)
library(dplyr)
library(sf)

# simpan hasil semua tahun
hasil_eb <- list()

# =========================================================
# LOOP SETIAP TAHUN
# =========================================================

for(th in sort(unique(jabar_join$tahun))){
  
  # data per tahun
  data_th <- jabar_join %>%
    
    filter(tahun == th)
  
  # =====================================================
  # EMPIRICAL BAYES SMOOTHING
  # =====================================================
  
  eb <- EBest(
    
    n = data_th$jumlah_bblr_imputasi,
    
    x = data_th$jumlah_bayi_lahir
  )
  
  # crude IR
  data_th$IR_crude <- (
    
    data_th$jumlah_bblr_imputasi /
      data_th$jumlah_bayi_lahir
    
  ) * 1000
  
  # empirical bayes smoothed IR
  data_th$EB <- eb$est * 1000
  
  # shrinkage
  data_th$shrinkage <- abs(
    
    data_th$IR_crude -
      data_th$EB
  )
  
  hasil_eb[[as.character(th)]] <- data_th
}

# =========================================================
# GABUNGKAN SEMUA TAHUN
# =========================================================

map_eb_all <- do.call(
  rbind,
  hasil_eb
)

# pastikan object sf
map_eb_all <- st_as_sf(map_eb_all)

# =========================================================
# RINGKASAN EMPIRICAL BAYES PER TAHUN
# =========================================================

ringkasan_eb <- map_eb_all %>%
  
  st_drop_geometry() %>%
  
  group_by(tahun) %>%
  
  summarise(
    
    `Range IR Kasar` = paste0(
      round(min(IR_crude, na.rm = TRUE), 2),
      " - ",
      round(max(IR_crude, na.rm = TRUE), 2)
    ),
    
    `Range EB IR` = paste0(
      round(min(EB, na.rm = TRUE), 2),
      " - ",
      round(max(EB, na.rm = TRUE), 2)
    ),
    
    `Shrinkage Rata-rata` = round(
      mean(shrinkage, na.rm = TRUE),
      4
    ),
    
    `Shrinkage Maksimum` = round(
      max(shrinkage, na.rm = TRUE),
      4
    )
  )

#ringkasan_eb

kable(
  ringkasan_eb,
  
  caption =
    "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 6. Ringkasan Empirical Bayes Smoothing
  </div>",
  
  align = "c",
  
  digits = 3
) %>%
  
  kable_styling(
    bootstrap_options = c(
      "striped",
      "hover",
      "condensed"
    ),
    
    full_width = FALSE,
    
    position = "center"
  )

Tabel 6. Ringkasan Empirical Bayes Smoothing
tahun	Range IR Kasar	Range EB IR	Shrinkage Rata-rata	Shrinkage Maksimum
2019	4.64 - 59.51	4.66 - 59.27	0.120	0.787
2020	4.85 - 58.18	4.88 - 57.95	0.145	0.864
2021	4.49 - 140.05	4.51 - 138.57	0.140	1.475
2022	3.3 - 186.92	3.34 - 185.2	0.137	1.722
2023	5.58 - 68.94	5.65 - 68.65	0.137	0.757

Berdasarkan ringkasan hasil Empirical Bayes smoothing, rentang nilai Empirical Bayes IR pada setiap tahun sedikit lebih sempit dibandingkan IR kasar. Nilai shrinkage rata-rata berada pada kisaran 0,120–0,145, sedangkan shrinkage maksimum tertinggi terjadi pada tahun 2022. Hasil tersebut menunjukkan adanya perubahan nilai risiko setelah dilakukan proses smoothing pada masing-masing wilayah.

# =========================================================
# PETA HASIL EMPIRICAL BAYES
# =========================================================

breaks_eb <- quantile(
  map_eb_all$IR_crude,
  
  probs = seq(0, 1, 0.2),
  
  na.rm = TRUE
)

# hindari duplikat
breaks_eb <- unique(breaks_eb)

# =========================================================
# LABEL RENTANG ANGKA
# =========================================================

label_eb <- paste0(
  
  round(
    head(breaks_eb, -1),
    1
  ),
  
  " - ",
  
  round(
    tail(breaks_eb, -1),
    1
  )
)

# =========================================================
# KATEGORI EMPIRICAL BAYES
# =========================================================

map_eb_all <- map_eb_all %>%
  
  mutate(
    
    kategori_EB = cut(
      EB,
      
      breaks = breaks_eb,
      
      include.lowest = TRUE,
      
      labels = label_eb
    )
  )

# =========================================================
# PLOT PETA EMPIRICAL BAYES
# =========================================================

ggplot(map_eb_all) +
  
  geom_sf(
    aes(fill = kategori_EB),
    
    color = "grey95",
    
    linewidth = 0.4
  ) +
  
  facet_wrap(
    ~tahun,
    ncol = 3
  ) +
  
  scale_fill_manual(
    
    values = c(
      "#FFF5EB",
      "#FDBE85",
      "#FD8D3C",
      "#E6550D",
      "#A63603"
    ),
    
    name = "EB Risk\nper 1.000"
  ) +
  
  labs(
    title = "Peta Empirical Bayes Smoothed Risk BBLR",
    
    subtitle = "Kabupaten/Kota di Jawa Barat Tahun 2019–2023",
    
  ) +
  
  theme_minimal(base_size = 14) +
  
  theme(
    
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 20
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40",
      size = 12
    ),
    
    strip.text = element_text(
      face = "bold",
      size = 12
    ),
    
    legend.title = element_text(
      face = "bold"
    ),
    
    legend.text = element_text(
      size = 10
    ),
    
    axis.text = element_blank(),
    
    axis.title = element_blank(),
    
    panel.grid = element_blank(),
    
    plot.caption = element_text(
      hjust = 0.5,
      color = "grey40"
    )
  )

Gambar 6. Peta Empirical Bayes

Berdasarkan peta Empirical Bayes smoothed risk kasus BBLR tahun 2019–2023, terlihat adanya perbedaan tingkat risiko antar kabupaten/kota di Jawa Barat. Wilayah dengan warna lebih gelap menunjukkan nilai Empirical Bayes risk yang lebih tinggi, sedangkan wilayah dengan warna lebih terang menunjukkan nilai risiko yang lebih rendah. Secara teoritis, metode Empirical Bayes akan memberikan koreksi yang lebih besar pada wilayah dengan populasi berisiko kecil karena estimasi risikonya cenderung lebih dipengaruhi oleh fluktuasi acak, sedangkan pada wilayah dengan populasi besar estimasi risiko cenderung stabil sehingga nilai hasil smoothing tetap mendekati nilai risiko kasar (Wakefield, 2020). Oleh karena itu, kemiripan antara peta crude risk dan peta Empirical Bayes risk.

3.5. Analisis Spasio-temporal

3.5.1. Regresi Poisson

# ====================================
# BUAT ID WILAYAH TETAP
# ====================================

mapping_wilayah <- data.frame(
  WADMKK = data_sp$WADMKK,
  id_wilayah = seq_len(length(data_sp$WADMKK))
)

# ====================================
# DATA MODEL
# ====================================

data_model <- jabar_join %>%
  
  st_drop_geometry() %>%
  
  left_join(
    mapping_wilayah,
    by = "WADMKK"
  ) %>%
  
  mutate(
    id_tahun = tahun - min(tahun) + 1
  )
# =========================================================
# POISSON
# =========================================================

model_pois <- glm(

  jumlah_bblr_imputasi ~

    persentase_penerima_ttd +
    persentase_sanitasi_layak +
    scale(jumlah_bidan_puskesmas) +
    air_minum_layak +
    angka_melek_aksara +
    tingkat_pengangguran_terbuka +
    scale(jumlah_ibu_hamil_zat_besi),

  family = poisson(),

  offset = log(expected_cases),

  data = data_model
)

3.5.2. Overdispersi

# =========================================================
# OVERDISPERSION
# =========================================================

dispersion <-

  model_pois$deviance /

  model_pois$df.residual

dispersion

## [1] 155.9748

Nilai overdispersi pada model Poisson sebesar 155.9748 menunjukkan bahwa varians data jauh lebih besar dibandingkan nilai rata-ratanya sehingga terjadi overdispersi pada data. Menurut Rahayu (2021), nilai statistik deviance atau Pearson Chi-Square yang dibagi dengan derajat bebas lebih besar dari 1 mengindikasikan adanya overdispersi pada data. Negative Binomial merupakan salah satu metode untuk mengatasi overdispersi pada regresi Poisson, dimana regresi dengan distribusi Negative Binomial tidak mengharuskan adanya kondisi equidispersion seperti pada regresi Poisson (Suryadi, 2023).

# =========================================================
# ADJACENCY MATRIX
# =========================================================

jabar_map <- jabar_shp

nb <- poly2nb(
  jabar_map,
  queen = TRUE
)

nb2INLA(
  "jabar.adj",
  nb
)

3.5.3. Analisis Spatio Temporal dengan INLA

Code Model Poisson INLA

# =========================================================
# MODEL POISSON INLA
# =========================================================

model_inla_pois <- inla(

  jumlah_bblr_imputasi ~

    persentase_penerima_ttd +
    persentase_sanitasi_layak +
    scale(jumlah_bidan_puskesmas) +
    air_minum_layak +
    angka_melek_aksara +
    tingkat_pengangguran_terbuka +
    scale(jumlah_ibu_hamil_zat_besi) +

    f(
  id_wilayah,
  model = "bym2",
  graph = "jabar.adj",
  scale.model = TRUE,

  hyper = list(

    prec = list(
      prior = "pc.prec",
      param = c(0.5, 0.01)
    ),

    phi = list(
      prior = "pc",
      param = c(0.5, 0.5)
    )
  )
) +

    f(
      id_tahun,
      model = "rw1"
    ),

  family = "poisson",

  offset = log(expected_cases),

  data = data_model,

  control.predictor =
    list(
      compute = TRUE
    ),

  control.compute =
    list(
      dic = TRUE,
      waic = TRUE
    )
)

Code Model Negative Binomial INLA

# =========================================================
# MODEL NEGATIVE BINOMIAL INLA
# =========================================================

model_inla_nb <- inla(

  jumlah_bblr_imputasi ~

    persentase_penerima_ttd +
    persentase_sanitasi_layak +
    scale(jumlah_bidan_puskesmas) +
    air_minum_layak +
    angka_melek_aksara +
    tingkat_pengangguran_terbuka +
    scale(jumlah_ibu_hamil_zat_besi) +

    f(
  id_wilayah,
  model = "bym2",
  graph = "jabar.adj",
  scale.model = TRUE,

  hyper = list(

    prec = list(
      prior = "pc.prec",
      param = c(0.5, 0.01)
    ),

    phi = list(
      prior = "pc",
      param = c(0.5, 0.5)
    )
  )
) +

    f(
      id_tahun,
      model = "rw1"
    ),

  family = "nbinomial",

  offset = log(expected_cases),

  data = data_model,

  control.predictor =
    list(
      compute = TRUE
    ),

  control.compute =
    list(
      dic = TRUE,
      waic = TRUE,
      cpo = TRUE
    )
)

Analisis spasial-temporal pada penelitian ini menggunakan pendekatan Bayesian melalui metode Integrated Nested Laplace Approximation (INLA). Model yang digunakan terdiri atas model Poisson dan Negative Binomial dengan komponen efek spasial dan temporal. Variabel respon berupa jumlah kasus BBLR hasil imputasi, sedangkan variabel prediktor meliputi persentase penerima tablet tambah darah, persentase sanitasi layak, jumlah bidan puskesmas, akses air minum layak, angka melek aksara, tingkat pengangguran terbuka, dan jumlah ibu hamil yang memperoleh zat besi.

Hubungan spasial antarwilayah dibentuk berdasarkan matriks ketetanggaan menggunakan pendekatan queen contiguity, dimana dua wilayah dianggap bertetangga apabila memiliki sisi maupun titik sudut yang saling bersinggungan. Struktur spasial kemudian dimasukkan ke dalam model menggunakan pendekatan Besag-York-Mollié 2 (BYM2), yaitu model yang menggabungkan efek spasial terstruktur dan tidak terstruktur. Efek spasial terstruktur digunakan untuk menangkap kemiripan risiko antarwilayah yang berdekatan, sedangkan efek spasial tidak terstruktur digunakan untuk menangkap variasi acak yang tidak dipengaruhi oleh lokasi geografis.

Pada model BYM2 digunakan hyperparameter berupa precision dan phi dengan pendekatan Penalized Complexity Prior (PC Prior). Hyperparameter precision menggunakan prior PC Prior dengan parameter (0,5; 0,01) untuk mengontrol tingkat keragaman efek acak spasial, dimana nilai precision yang lebih besar menunjukkan variasi spasial yang lebih kecil. Sementara itu, hyperparameter phi menggunakan prior PC Prior dengan parameter (0,5; 0,5) untuk mengatur proporsi kontribusi antara efek spasial terstruktur dan tidak terstruktur dalam model. Penggunaan PC Prior bertujuan untuk menghindari model yang terlalu kompleks dan menghasilkan estimasi yang lebih stabil (Simpson dkk., 2017).

Model juga memasukkan efek temporal menggunakan random walk orde satu. Pendekatan ini mengasumsikan bahwa kondisi pada suatu tahun dipengaruhi oleh kondisi pada tahun sebelumnya sehingga mampu menggambarkan perubahan risiko BBLR dari waktu ke waktu secara lebih halus. Analisis juga menggunakan offset berdasarkan jumlah kasus harapan.

# =========================================================
# DIC & WAIC
# =========================================================

perbandingan_model <- data.frame(

  Model = c(
    "Poisson",
    "Negative Binomial"
  ),

  DIC = c(
    model_inla_pois$dic$dic,
    model_inla_nb$dic$dic
  ),

  WAIC = c(
    model_inla_pois$waic$waic,
    model_inla_nb$waic$waic
  )
)

kable(
  perbandingan_model,
  caption = 
    "<div style='text-align:center; color:black; font-weight:bold;'>
Tabel 7. Perbandingan DIC dan WAIC  
  </div>"
)

Table:

Tabel 7. Perbandingan DIC dan WAIC

Model	DIC	WAIC
Poisson	-926.0363	6073.631
Negative Binomial	1816.1443	1831.573

Pemilihan model terbaik dilakukan menggunakan nilai Deviance Information Criterion (DIC) dan Widely Applicable Information Criterion (WAIC), dimana model dengan nilai yang lebih kecil dianggap memiliki kecocokan yang lebih baik (Spiegelhalter dkk., 2002). Hasil perbandingan menunjukkan bahwa model Poisson memiliki nilai DIC yang lebih rendah dibandingkan model Negative Binomial, sedangkan model Negative Binomial memiliki nilai WAIC yang jauh lebih rendah dibandingkan model Poisson. Karena data BBLR juga mengalami overdispersi sehingga model Negative Binomial lebih sesuai digunakan karena mampu mengakomodasi varians yang lebih besar daripada mean (Hilbe, 2011).

# =========================================================
# PARAMETER MODEL
# =========================================================

hasil_fixed <-
model_inla_nb$summary.fixed %>%

select(
mean,
sd,
`0.025quant`,
`0.5quant`,
`0.975quant`,
mode
)

kable(
round(
hasil_fixed,
4
),
caption = 
  "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 8. Parameter Fixed Effect
  </div>"
)

Table:

Tabel 8. Parameter Fixed Effect

	mean	sd	0.025quant	0.5quant	0.975quant	mode
(Intercept)	-5.4978	4.7154	-14.8823	-5.4527	3.6346	-5.4515
persentase_penerima_ttd	0.0027	0.0034	-0.0039	0.0027	0.0093	0.0027
persentase_sanitasi_layak	0.0006	0.0017	-0.0028	0.0006	0.0039	0.0006
scale(jumlah_bidan_puskesmas)	0.1436	0.0958	-0.0458	0.1440	0.3306	0.1440
air_minum_layak	0.0032	0.0100	-0.0162	0.0032	0.0228	0.0032
angka_melek_aksara	0.0516	0.0437	-0.0330	0.0511	0.1387	0.0511
tingkat_pengangguran_terbuka	-0.0100	0.0191	-0.0476	-0.0100	0.0275	-0.0100
scale(jumlah_ibu_hamil_zat_besi)	-0.1602	0.1020	-0.3630	-0.1591	0.0377	-0.1588

Berdasarkan hasil estimasi parameter fixed effect pada model spasial-temporal Negative Binomial INLA, tidak terdapat variabel yang signifikan secara statistik karena seluruh interval kredibel 95% masih mencakup nilai nol. Nilai koefisien positif menunjukkan kecenderungan peningkatan risiko BBLR, sedangkan koefisien negatif menunjukkan kecenderungan penurunan risiko BBLR. Variabel persentase penerima tablet tambah darah memiliki koefisien sebesar 0,0027, persentase sanitasi layak sebesar 0,0006, jumlah bidan puskesmas sebesar 0,1436, akses air minum layak sebesar 0,0032, dan angka melek aksara sebesar 0,0516 yang menunjukkan kecenderungan peningkatan risiko BBLR. Sementara itu, tingkat pengangguran terbuka memiliki koefisien sebesar -0,0100 dan jumlah ibu hamil yang memperoleh zat besi sebesar -0,1603 yang menunjukkan kecenderungan penurunan risiko BBLR.

# ====================================
# RELATIVE RISK PER TAHUN
# ====================================

# fitted value dari model
fitted_mean <- 
  model_inla_nb$
  summary.fitted.values$mean

# hitung RR
data_rr <- data_model %>%
  
  mutate(
    
    RR =
      fitted_mean /
      expected_cases
  )

# ambil variabel yang diperlukan
rr_map <- data_rr %>%
  
  select(
    WADMKK,
    tahun,
    RR
  )

# ====================================
# GABUNGKAN KE SHAPEFILE
# ====================================

map_final <- jabar_join %>%
  
  left_join(
    rr_map,
    by = c(
      "WADMKK",
      "tahun"
    )
  )

# ====================================
# PETA RR
# ====================================

ggplot(map_final) +

  geom_sf(
    aes(fill = RR),
    color = "white"
  ) +

  facet_wrap(
    ~tahun,
    ncol = 3
  ) +

  scale_fill_viridis_c() +

  labs(
    title = "Peta Relative Risk BBLR",
    subtitle = "Kabupaten/Kota di Jawa Barat Tahun 2019–2023",
    fill = "RR"
  ) +

  theme_minimal(base_size = 14) +

  theme(
    
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 18
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40"
    ),
    
    strip.text = element_text(
      face = "bold"
    ),
    
    axis.text = element_blank(),
    
    axis.ticks = element_blank(),
    
    axis.title = element_blank(),
    
    panel.grid = element_blank()
  )

Gambar 7. Peta RR Model Spasio-Temporal

Berdasarkan peta relative risk (RR) BBLR hasil model spasial-temporal INLA Negative Binomial tahun 2019–2023, terlihat adanya variasi risiko kejadian BBLR antar kabupaten/kota di Jawa Barat. Wilayah dengan RR lebih dari 1 menunjukkan jumlah kasus BBLR yang lebih tinggi dibandingkan kasus yang diharapkan, sedangkan RR kurang dari 1 menunjukkan risiko yang lebih rendah. Secara umum, pola risiko relatif BBLR cenderung konsisten dari tahun ke tahun, dimana beberapa wilayah di bagian timur dan tenggara Jawa Barat memiliki risiko relatif lebih tinggi dibandingkan wilayah lainnya. Dengan demikian, wilayah dengan nilai RR yang tinggi dapat menjadi prioritas dalam upaya pencegahan dan penanganan kasus BBLR. Tingginya risiko BBLR berkaitan dengan faktor sosio demografi (Aini dan Kurniawan, 2023)

# =========================================================
# HEATMAP RELATIVE RISK (RR)
# =========================================================

library(ggplot2)
library(dplyr)
library(forcats)

# ====================================
# DATA RR DARI MODEL INLA
# ====================================

fitted_mean <-
  model_inla_nb$
  summary.fitted.values$mean

heatmap_rr <- data_model %>%

  mutate(
    RR = fitted_mean / expected_cases
  ) %>%

  select(
    WADMKK,
    tahun,
    RR
  )

# ====================================
# URUTKAN DAERAH BERDASARKAN
# RR RATA-RATA
# ====================================

urutan_daerah <- heatmap_rr %>%

  group_by(WADMKK) %>%

  summarise(
    RR_mean = mean(RR)
  ) %>%

  arrange(desc(RR_mean)) %>%

  pull(WADMKK)

heatmap_rr$WADMKK <-
  factor(
    heatmap_rr$WADMKK,
    levels = urutan_daerah
  )

# ====================================
# HEATMAP
# ====================================

ggplot(
  heatmap_rr,
  aes(
    x = factor(tahun),
    y = WADMKK,
    fill = RR
  )
) +

  geom_tile(
    color = "white",
    linewidth = 0.3
  ) +

  scale_fill_viridis_c(
    option = "magma",
    name = "Relative\nRisk (RR)"
  ) +

  labs(
    title = "Heatmap Relative Risk (RR) BBLR",
    subtitle = "Hasil Model Spasio-Temporal Negative Binomial INLA",
    x = "Tahun",
    y = "Kabupaten/Kota"
  ) +

  theme_minimal(base_size = 13) +

  theme(

    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 18
    ),

    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40"
    ),

    axis.title = element_text(
      face = "bold"
    ),

    axis.text.y = element_text(
      size = 8
    ),

    panel.grid = element_blank()
  )

Gambar 8. Heatmap RR BBLR

Berdasarkan heatmap Relative Risk (RR) BBLR tahun 2019–2023, terlihat bahwa risiko BBLR berbeda antar kabupaten/kota di Jawa Barat dan cenderung stabil dari tahun ke tahun. Kabupaten Tasikmalaya dan Kabupaten Kuningan secara konsisten memiliki RR tertinggi, diikuti oleh Indramayu, Pangandaran, Majalengka, Ciamis, dan Kota Cirebon. Sebaliknya, Kabupaten Bekasi, Kota Bekasi, Kota Depok, dan Bogor menunjukkan RR yang relatif lebih rendah. Temuan ini menunjukkan adanya wilayah prioritas yang memerlukan perhatian lebih dalam upaya pencegahan dan penanggulangan BBLR.

# =========================================================
# TOP 5 WILAYAH DENGAN RR TERTINGGI
# SETIAP TAHUN
# =========================================================


# fitted mean dari model INLA
fitted_mean <- 
  model_inla_nb$
  summary.fitted.values$mean

# data RR
data_rr <- data_model %>%
  
  mutate(
    
    RR =
      fitted_mean /
      expected_cases
  )

# =========================================================
# AMBIL TOP 5 RR SETIAP TAHUN
# =========================================================

top5_rr <- data_rr %>%
  
  group_by(tahun) %>%
  
  slice_max(
    order_by = RR,
    n = 5
  ) %>%
  
  ungroup() %>%
  
  arrange(
    tahun,
    desc(RR)
  ) %>%
  
  mutate(
    
    RR = round(RR, 3)
  ) %>%
  
  select(
    tahun,
    WADMKK,
    RR
  )

# =========================================================
# VISUALISASI TOP 5 RR
# =========================================================

ggplot(
  top5_rr,
  
  aes(
    x = fct_reorder(
      WADMKK,
      RR
    ),
    
    y = RR,
    
    fill = factor(tahun)
  )
) +
  
  geom_col(
    alpha = 0.9
  ) +
  
  coord_flip() +
  
  facet_wrap(
    ~tahun,
    scales = "free_y"
  ) +
  
  geom_hline(
    yintercept = 1,
    linetype = "dashed",
    linewidth = 1
  ) +
  
  geom_text(
    aes(
      label = round(RR, 2)
    ),
    
    hjust = -0.1,
    
    size = 3.5,
    
    fontface = "bold"
  ) +
  
  labs(
    
    title =
      "Top 5 Wilayah dengan Relative Risk (RR) Tertinggi",
    
    subtitle =
      "Kasus BBLR Kabupaten/Kota di Jawa Barat Tahun 2019–2023",
    
    x =
      "Kabupaten/Kota",
    
    y =
      "Relative Risk (RR)",
    
    fill =
      "Tahun"
  ) +
  
  theme_minimal(
    base_size = 13
  ) +
  
  theme(
    
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      size = 18
    ),
    
    plot.subtitle = element_text(
      hjust = 0.5,
      color = "grey40"
    ),
    
    axis.title = element_text(
      face = "bold"
    ),
    
    legend.position = "right",
    
    legend.title = element_text(
      face = "bold"
    ),
    
    panel.grid.minor =
      element_blank()
  )

Gambar 9. Top 5 Wilayah dengan RR RRTertinggi

Berdasarkan gambar Top 5 Wilayah dengan Relative Risk (RR) Tertinggi tahun 2019–2023, Kabupaten Tasikmalaya dan Kabupaten Kuningan secara konsisten memiliki nilai RR tertinggi dengan nilai lebih dari 2 pada hampir seluruh periode pengamatan. Hal ini menunjukkan bahwa risiko kejadian BBLR di kedua wilayah tersebut lebih dari dua kali lipat dibandingkan risiko rata-rata Provinsi Jawa Barat. Selain itu, Kabupaten Indramayu, Pangandaran, Majalengka, Kota Cirebon, dan Ciamis juga beberapa kali muncul sebagai wilayah dengan RR tinggi. Temuan ini mengindikasikan adanya wilayah-wilayah yang secara konsisten memiliki risiko BBLR lebih tinggi sehingga perlu menjadi prioritas dalam upaya pencegahan dan penanganan BBLR di Jawa Barat.

Hasil penelitian menunjukkan bahwa kejadian BBLR di Jawa Barat masih memiliki variasi antarwilayah, sehingga diperlukan intervensi yang lebih terfokus pada daerah dengan risiko tinggi. Penguatan pelayanan kesehatan ibu dan anak, pemantauan kehamilan berisiko, serta optimalisasi pemberian tablet tambah darah dan suplementasi zat besi perlu terus ditingkatkan untuk menurunkan kejadian BBLR (World Health Organization, 2023). Selain itu, pemanfaatan analisis spasial dan spasio-temporal dapat mendukung penentuan wilayah prioritas sehingga program kesehatan dan alokasi sumber daya dapat dilakukan secara lebih tepat sasaran (MacNab, 2022).

3.6. Uji Autokorelasi Spasial Residual Model

# =========================================================
# UJI MORAN'S I PADA RESIDUAL MODEL NEGATIVE BINOMIAL INLA
# =========================================================

# residual model
data_model$residual <-

  data_model$jumlah_bblr_imputasi -

  model_inla_nb$summary.fitted.values$mean

# =========================================================
# MORAN'S I RESIDUAL PER TAHUN
# =========================================================

hasil_moran_res <- data.frame()

for(th in sort(unique(data_model$tahun))){

  data_th <- data_model %>%
    
    filter(
      tahun == th
    )

  moran_res <- moran.test(
    
    data_th$residual,
    
    lw,
    
    zero.policy = TRUE
  )

  hasil_moran_res <- rbind(
    
    hasil_moran_res,
    
    data.frame(
      
      Tahun = th,
      
      Moran_I = as.numeric(
        moran_res$estimate[1]
      ),
      
      Expected_I = as.numeric(
        moran_res$estimate[2]
      ),
      
      Variance = as.numeric(
        moran_res$estimate[3]
      ),
      
      P_Value = moran_res$p.value
    )
  )
}

# =========================================================
# TABEL HASIL MORAN'S I RESIDUAL
# =========================================================

hasil_moran_res_tabel <- hasil_moran_res %>%
  
  mutate(
    
    Moran_I = round(
      Moran_I,
      4
    ),
    
    Expected_I = round(
      Expected_I,
      4
    ),
    
    Variance = round(
      Variance,
      4
    ),
    
    P_Value = round(
      P_Value,
      4
    ),
    
    Interpretasi =
      
      ifelse(
        
        P_Value < 0.05 &
          Moran_I > 0,
        
        "Autokorelasi Positif",
        
        ifelse(
          
          P_Value < 0.05 &
            Moran_I < 0,
          
          "Autokorelasi Negatif",
          
          "Tidak Signifikan"
        )
      )
  )

kable(
  
  hasil_moran_res_tabel,
  
  caption =
    "<div style='text-align:center; color:black; font-weight:bold;'>
  Tabel 9. Hasil Uji Moran's I Residual Model Negative Binomial INLA
  </div>",
  
  align = "c"
  
) %>%
  
  kable_styling(
    
    full_width = FALSE,
    
    position = "center"
  )

Tabel 9. Hasil Uji Moran’s I Residual Model Negative Binomial INLA
Tahun	Moran_I	Expected_I	Variance	P_Value	Interpretasi
2019	0.1707	-0.0385	0.0165	0.0518	Tidak Signifikan
2020	-0.0235	-0.0385	0.0135	0.4488	Tidak Signifikan
2021	-0.2685	-0.0385	0.0173	0.9598	Tidak Signifikan
2022	0.0468	-0.0385	0.0151	0.2438	Tidak Signifikan
2023	-0.1098	-0.0385	0.0153	0.7181	Tidak Signifikan

Berdasarkan hasil uji Moran’s I terhadap residual model Negative Binomial INLA, seluruh nilai p-value lebih besar dari 0,05 sehingga tidak ditemukan autokorelasi spasial yang signifikan pada residual untuk seluruh tahun pengamatan (2019–2023). Hasil ini menunjukkan bahwa komponen spasial BYM2 yang dimasukkan ke dalam model telah mampu mengakomodasi ketergantungan spasial antar kabupaten/kota di Jawa Barat. Dengan demikian, pola spasial yang sebelumnya terdapat pada data BBLR telah berhasil dijelaskan oleh model, sehingga residual yang tersisa cenderung tersebar secara acak dan tidak menunjukkan pengelompokan spasial yang signifikan.

4. Implikasi

4.1. Temuan

Kejadian BBLR di Jawa Barat menunjukkan pola persebaran yang tidak merata antarwilayah dan adanya pengelompokan spasial pada beberapa tahun pengamatan. Hal ini menunjukkan bahwa risiko BBLR dipengaruhi oleh karakteristik wilayah sehingga memerlukan pendekatan penanganan yang mempertimbangkan aspek geografis.

Berdasarkan hasil pemetaan risiko relatif, Kabupaten Tasikmalaya dan Kabupaten Kuningan termasuk wilayah yang secara konsisten memiliki risiko BBLR lebih tinggi dibandingkan wilayah lainnya. Oleh karena itu, kedua wilayah tersebut dapat menjadi prioritas dalam upaya pencegahan dan pengendalian BBLR.

Penelitian ini memiliki beberapa keterbatasan. Pertama, analisis menggunakan data agregat tingkat kabupaten/kota sehingga belum dapat menggambarkan variasi risiko pada tingkat wilayah yang lebih kecil. Kedua, terdapat kemungkinan underreporting atau ketidaklengkapan pelaporan kasus BBLR pada beberapa wilayah, yang dapat memengaruhi estimasi risiko yang dihasilkan. Ketiga, variabel yang digunakan masih terbatas pada data sekunder yang tersedia sehingga belum mencakup seluruh faktor risiko individu yang berkaitan dengan kejadian BBLR.

Pemanfaatan disease mapping dan analisis spasio-temporal dapat mendukung pemantauan BBLR secara berkala untuk mengidentifikasi wilayah prioritas. Selain itu, peningkatan kualitas pencatatan dan pelaporan data kesehatan diperlukan agar perencanaan program kesehatan dapat dilakukan secara lebih tepat sasaran.

4.2. Rekomendasi Kebijakan

Penguatan pelayanan antenatal terpadu, terutama melalui deteksi dini anemia, kekurangan energi kronis (KEK), dan komplikasi kehamilan yang merupakan faktor risiko utama terjadinya BBLR (World Health Organization, 2023).
Peningkatan cakupan dan kepatuhan konsumsi tablet tambah darah serta suplementasi zat besi pada ibu hamil untuk mencegah anemia dan menurunkan risiko kelahiran BBLR.
Peningkatan edukasi gizi selama kehamilan serta pemantauan status gizi ibu hamil secara berkala guna mendukung pertumbuhan janin yang optimal dan mengurangi risiko BBLR (UNICEF & WHO, 2023).
Optimalisasi kunjungan antenatal sesuai standar agar kehamilan berisiko dapat diidentifikasi dan ditangani lebih dini melalui intervensi yang tepat.
Penguatan sistem surveilans kesehatan berbasis spasial melalui pemanfaatan disease mapping dan analisis spasio-temporal untuk memantau wilayah berisiko tinggi dan mendukung penentuan prioritas program kesehatan ibu dan anak.

Kesimpulan

Berdasarkan hasil analisis disease mapping, kejadian Bayi Berat Lahir Rendah (BBLR) di Provinsi Jawa Barat menunjukkan pola persebaran yang tidak merata antar kabupaten/kota. Hasil Moran’s I menunjukkan adanya autokorelasi spasial positif pada sebagian besar tahun pengamatan, yang mengindikasikan bahwa wilayah dengan kejadian BBLR tinggi cenderung berdekatan dengan wilayah yang juga memiliki kejadian tinggi. Peta risiko dan hasil Empirical Bayes smoothing juga menunjukkan adanya variasi risiko BBLR antarwilayah, dengan beberapa kabupaten/kota secara konsisten berada pada kelompok risiko yang relatif lebih tinggi.

Berdasarkan analisis spasio-temporal menggunakan model Bayesian INLA dengan komponen BYM2 dan Random Walk Order 1 (RW1), diperoleh bahwa kejadian BBLR di Jawa Barat memiliki variasi menurut dimensi ruang dan waktu. Hasil pemetaan relative risk menunjukkan bahwa Kabupaten Tasikmalaya dan Kabupaten Kuningan termasuk wilayah yang secara konsisten memiliki risiko lebih tinggi dibandingkan wilayah lainnya selama periode 2019–2023. Selain itu, hasil uji Moran’s I residual menunjukkan tidak adanya autokorelasi spasial yang signifikan pada residual model, sehingga model yang digunakan telah mampu mengakomodasi pola spasial yang terdapat pada data.

Hasil penelitian ini menunjukkan bahwa upaya penanggulangan BBLR perlu dilakukan secara lebih terarah pada wilayah dengan risiko tinggi melalui penguatan pelayanan kesehatan ibu dan anak, peningkatan pemantauan kehamilan berisiko, serta optimalisasi program pencegahan BBLR. Selain itu, pendekatan disease mapping dan analisis spasio-temporal dapat digunakan sebagai dasar dalam menentukan wilayah prioritas sehingga perencanaan program kesehatan dan alokasi sumber daya dapat dilakukan secara lebih efektif dan tepat sasaran.

Daftar Pustaka

Aini, Y., & Kurniawan, F. (2023). The maternal sociodemographic determinants of low birth weight in Indonesia. Jurnal Kesehatan Masyarakat, 18(4), 536–545. https://doi.org/10.15294/kemas.v18i4.42006
Blencowe, H., Krasevec, J., de Onis, M., Black, R. E., An, X., Stevens, G. A., Borghi, E., Hayashi, C., Estevez, D., Cegolon, L., Shiekh, S., Ponce Hardy, V., Lawn, J. E., & Cousens, S. (2019). National, regional, and worldwide estimates of low birthweight in 2015, with trends from 2000: A systematic analysis. The Lancet Global Health, 7(7), e849–e860. https://doi.org/10.1016/S2214-109X(18)30565-5
Badan Pusat Statistik Provinsi Jawa Barat. (2023). Jawa Barat dalam angka 2023. Badan Pusat Statistik Provinsi Jawa Barat.
Cameron, A. C., & Trivedi, P. K. (2013). Regression analysis of count data (2nd ed.). Cambridge University Press.
Dinas Kesehatan Provinsi Jawa Barat. (2023). Profil kesehatan Provinsi Jawa Barat tahun 2022. Dinas Kesehatan Provinsi Jawa Barat.
Hilbe, J. M. (2011). Negative binomial regression (2nd ed.). Cambridge University Press.
Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19(17–18), 2555–2567. https://doi.org/10.1002/1097-0258(20000915/30)19:17/18<2555::AID-SIM587>3.0.CO;2-#
Kementerian Kesehatan Republik Indonesia. (2023). Profil kesehatan Indonesia tahun 2022. Kementerian Kesehatan Republik Indonesia.
Lawson, A. B. (2018). Bayesian disease mapping: Hierarchical modeling in spatial epidemiology (3rd ed.). CRC Press.
Marshall, R. J. (1991). Mapping disease and mortality rates using empirical Bayes estimators. Journal of the Royal Statistical Society: Series C (Applied Statistics), 40(2), 283–294.
Rahayu, A. (2020). Model-model regresi untuk mengatasi masalah overdispersi pada regresi Poisson. Peqguruang: Conference Series, 2(1), 1–6.
Riebler, A., Sørbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145–1165. https://doi.org/10.1177/0962280216660421
Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2), 319–392. https://doi.org/10.1111/j.1467-9868.2008.00700.x
Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sørbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science, 32(1), 1–28. https://doi.org/10.1214/16-STS576
Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583–639. https://doi.org/10.1111/1467-9868.00353
Suryadi, F., Jonathan, S., Jonatan, K., & Ohyver, M. (2023). Handling overdispersion in Poisson regression using negative binomial regression for poverty case in West Java. Procedia Computer Science, 216, 517–523. https://doi.org/10.1016/j.procs.2023.01.319
Tshotetsi, L., Dzikiti, L., Hajison, P., & Feresu, S. (2019). Maternal factors contributing to low birth weight deliveries in Tshwane District, South Africa. PLoS ONE, 14(3), e0213058. https://doi.org/10.1371/journal.pone.0213058
UNICEF, & World Health Organization. (2023). Low birth weight estimates: Levels and trends 2000–2020. UNICEF & WHO.
Waller, L. A., & Gotway, C. A. (2004). Applied spatial statistics for public health data. John Wiley & Sons.
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571–3594.
World Health Organization. (2023). Low birth weight estimates: Levels and trends 2000–2020. World Health Organization.
Wakefield, J. (2020). Small Area Estimation for Disease Prevalence Mapping. International Statistical Review, 88(S1), S40–S66.
MacNab, Y. C. (2022). Bayesian disease mapping: Past, present, and future. Spatial Statistics, 50, 100593.
Badan Informasi Geospasial. (2024). Batas administrasi kabupaten/kota Provinsi Jawa Barat (shapefile). Pusat Pemetaan Batas Wilayah (PPBW), Badan Informasi Geospasial.
Open Data Jabar. (2024). Jumlah bayi berat badan lahir rendah (BBLR) berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/jumlah-bayi-berat-badan-lahir-rendah-bblr-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Persentase ibu hamil penerima tablet tambah darah (TTD) berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/persentase-ibu-hamil-penerima-tablet-tambah-darah-ttd-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Persentase keluarga dengan akses sanitasi layak (jamban sehat) berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/persentase-keluarga-dengan-akses-sanitasi-layak-jamban-sehat-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Jumlah tenaga kebidanan berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/jumlah-tenaga-kebidanan-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Persentase rumah tangga yang memiliki akses terhadap sumber air minum layak berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/persentase-rumah-tangga-yang-memiliki-akses-terhadap-sumber-air-minum-layak-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Tingkat pengangguran terbuka berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/en/dataset/tingkat-pengangguran-terbuka-berdasarkan-kabupatenkota-di-jawa-barat
Open Data Jabar. (2024). Jumlah ibu hamil mendapatkan tablet zat besi (FE3-90 tablet) berdasarkan kabupaten/kota di Jawa Barat. Pemerintah Provinsi Jawa Barat. https://opendata.jabarprov.go.id/id/dataset/jumlah-ibu-hamil-mendapatkan-tablet-zat-besi-fe3-90-tablet-berdasarkan-kabupatenkota-di-jawa-barat
Badan Pusat Statistik Provinsi Jawa Barat. (2024). Angka Melek Aksara Penduduk Umur 15 Tahun ke Atas Menurut Kabupaten/Kota di Provinsi Jawa Barat (Percent). Badan Pusat Statistik Provinsi Jawa Barat. https://jabar.bps.go.id/en/statistics-table/2/OTU2IzI=/angka-melek-aksara-penduduk-umur-15-tahun-ke-atas-menurut-kabupaten-kota-di-provinsi-jawa-barat.html

Pemetaan Risiko dan Analisis Spasio-Temporal Bayi Berat Lahir Rendah (BBLR) di Provinsi Jawa Barat Tahun 2019–2023

Nama : Rana Aribah NPM : 140720260012 Dosen Pengampu : Prof. Dr. I Gede Nyoman Mindra Jaya, M.Si.

2026-06-16

Abstrak

1. Pendahuluan

1.1. Latar Belakang

1.2. Rumusan Masalah

1.3. Tujuan Penelitian

2. Data dan Metode

2.1. Sumber Data dan Variabel Penelitian

2.2. Ukuran Epidemiologi

2.2.1. Incidence Rate

2.2.2. Standardized Incidence Ratio

2.3. Empirical Bayes Smoothing

2.4. Autokorelasi Spasial

2.5. Regresi Poisson

2.6. Pemeriksaan Overdispersi

2.7. Model Spasio-Temporal Bayesian INLA

2.7.1. Efek Spasial BYM2

2.7.2. Efek Temporal Random Walk Order 1 (RW1)

2.8. Pemilihan Model

2.8.1. Deviance Information Criterion (DIC)

2.8.2. Widely Applicable Information Criterion (WAIC)

3. Hasil dan Pembahasan

3.1. Cek Kualitas Data

3.2. Statistik Desktriptif dan Exploratory Data Analysis

3.2.1. Statistik Deskriptif

3.2.2. Grafik Tren Kasus BBLR

3.3. Ukuran Epidemiologi

3.3.1. Perhitungan IR dan SIR

3.3.2. Visualisasi Incidence Rate per Tahun

3.4. Disease Mapping

3.4.1. Analisis Autokorelasi Spasial Global Moran’s I Berdasarkan IR

3.4.2. Peta Crude Risk Berdasarkan IR

3.4.3. Empirical Bayes Smoothing

3.5. Analisis Spasio-temporal

3.5.1. Regresi Poisson

3.5.2. Overdispersi

3.5.3. Analisis Spatio Temporal dengan INLA

3.6. Uji Autokorelasi Spasial Residual Model

4. Implikasi

4.1. Temuan

4.2. Rekomendasi Kebijakan

Kesimpulan

Daftar Pustaka

Nama : Rana Aribah
NPM : 140720260012
Dosen Pengampu : Prof. Dr. I Gede Nyoman Mindra Jaya, M.Si.