knitr::include_graphics("C:/Users/Admin/Downloads/KULIAH/Semester 5/Epidem/UAS & Project (Rshiny)/Cover Rpubs UAS Epidem.png")

Abstrak

Lonjakan kasus campak global sebesar 20% menjadi ancaman serius bagi pencapaian target Sustainable Development Goals (SDGs) poin 3.2 dan 3.b, menempatkan Jawa Barat sebagai wilayah dengan kerentanan tinggi akibat adanya celah imunitas pasca-pandemi. Penelitian ini bertujuan menganalisis pola sebaran spasial dan determinan risiko kejadian campak di Jawa Barat tahun 2024 dengan menggunakan variabel imunisasi, sanitasi layak, dan kepadatan penduduk melalui komparasi metode Ordinary Least Squares (OLS) dan Ekonometrika Spasial. Berdasarkan hasil analisis, ditemukan bahwa meskipun secara visual kasus terkonsentrasi di wilayah urban utara, uji diagnostik spasial (Moran’s I) membuktikan pola penyebaran bersifat acak (random) tanpa adanya dependensi antar-wilayah yang signifikan, sehingga model OLS terbukti lebih efisien (AIC terendah) dibandingkan model spasial. Secara statistik, variabel kepadatan penduduk dan sanitasi terkonfirmasi sebagai faktor pendorong utama peningkatan kasus, mengindikasikan bahwa strategi pengendalian wabah tidak cukup hanya mengandalkan target vaksinasi, tetapi harus disertai intervensi perbaikan kualitas lingkungan di permukiman padat penduduk.

Bab 1. Pendahuluan

1.1 Latar Belakang

Penyakit campak (Measles) kembali menjadi perhatian serius bagi kesehatan global. Laporan terbaru dari organisasi kesehatan dunia (World Health Organization) dan CDC pada November 2024 menunjukkan data yang cukup mengkhawatirkan: kasus campak di seluruh dunia naik sekitar 20% pada tahun 2023, dengan perkiraan total mencapai 10,3 juta kasus [1]. Walaupun angka kematian sedikit turun, penyakit ini masih menyebabkan sekitar 107.500 orang meninggal, dan sebagian besarnya adalah anak-anak balita [1]. Kondisi ini tentu menghambat pencapaian target Sustainable Development Goals (SDGs) poin 3.2 tentang mengakhiri kematian balita yang bisa dicegah, serta poin 3.b tentang akses vaksin untuk semua [2].

Indonesia pun termasuk negara yang rentan. Data WHO memasukkan Indonesia ke dalam 10 besar negara dengan kasus campak terbanyak, sejajar dengan negara seperti Nigeria dan Pakistan [3]. Cakupan imunisasi dasar yang belum merata setelah pandemi membuat adanya celah kekebalan (immunity gap) yang berisiko. Di tingkat daerah, Jawa Barat sebagai provinsi dengan penduduk terbanyak punya risiko penularan yang lebih tinggi. Masalah penularan campak di sini tidak cuma soal virusnya, tapi juga dipengaruhi lingkungan dan kondisi sosial. Karena itu, penelitian ini melihat tiga faktor utama: (1) Persentase Imunisasi, sebagai pertahanan utama herd immunity untuk memutus rantai penularan; (2) Sanitasi Layak, yang bisa menggambarkan kondisi permukiman padat atau kumuh di mana penyakit menular lebih mudah menyebar [4]; dan (3) Kepadatan Penduduk, yang secara teori membuat interaksi antar-orang makin sering sehingga penyakit yang menular lewat udara seperti campak lebih mudah berpindah [5].

Akan tetapi, analisis biasa menggunakan regresi linier standar (OLS) sering kali kurang tepat kalau dipakai untuk data wilayah. Masalahnya, data antar-wilayah sering kali saling berhubungan dan tidak berdiri sendiri. Hal ini sangat terasa di Jawa Barat, terutama di daerah penyangga ibu kota (Bodebek). Data BPS tahun 2024 mencatat ada sekitar 7,59 juta pekerja yang setiap hari pergi-pulang antar kota/kabupaten [6]. Tingginya pergerakan orang ini membuat jumlah penduduk di suatu wilayah bisa berubah drastis pada waktu yang berbeda, sehingga virus sangat mungkin terbawa lintas wilayah.

Karena data ini punya keterkaitan antar-lokasi (spatial dependence), metode Ekonometrika Spasial sangat perlu digunakan. Uji Moran’s I dan Local Indicator of Spatial Association (LISA) dibutuhkan untuk melihat apakah penyebaran penyakit ini mengelompok di daerah tertentu (cluster) atau menyebar secara acak. Setelah itu, model Spatial Autoregressive Model (SAR) atau Spatial Error Model (SEM) akan digunakan untuk menghitung risiko dengan lebih akurat dibanding metode biasa, karena metode ini memperhitungkan pengaruh wilayah tetangga atau hubungan antar-wilayah lainnya [7]. Harapannya, pendekatan ini bisa menghasilkan peta risiko yang lebih jelas sebagai dasar penanganan kesehatan masyarakat yang lebih tepat sasaran.

1.2 Rumusan Masalah

  1. Bagaimana pola distribusi spasial dan tingkat kerawanan kasus penyakit campak di Kabupaten/Kota Provinsi Jawa Barat tahun 2024?
  2. Apakah terdapat autokorelasi spasial global (Moran’s I) dan lokal (LISA) yang signifikan pada insidensi campak, yang mengindikasikan adanya pengelompokan (clustering) wilayah wabah?
  3. Bagaimana pengaruh variabel prediktor (persentase imunisasi, akses sanitasi layak, dan kepadatan penduduk) terhadap peningkatan kasus campak?
  4. Model manakah di antara Ordinary Least Square (OLS), Spatial Autoregressive Model (SAR), atau Spatial Error Model (SEM) yang memiliki performa terbaik (berdasarkan nilai AIC dan uji diagnostik) dalam menjelaskan fenomena penyebaran campak di Jawa Barat?

1.3 Tujuan Penelitian

  1. Mendeskripsikan karakteristik epidemiologi dan peta sebaran kasus campak di Jawa Barat secara visual dan statistik.
  2. Mengidentifikasi keberadaan ketergantungan spasial dan lokasi hotspot (wilayah High-High) penyebaran campak menggunakan uji statistik spasial.
  3. Menganalisis signifikansi dan besaran pengaruh faktor risiko imunisasi, sanitasi, dan kepadatan penduduk terhadap insidensi penyakit.
  4. Mengevaluasi dan memilih model regresi spasial terbaik yang mampu menangkap efek interaksi antar-wilayah, guna memberikan rekomendasi kebijakan pengendalian penyakit yang berbasis bukti (evidence-based policy).

Bab 2. Tinjauan Pustaka

2.1 Agent-Host-Environment

Dalam memahami pola penyebaran dan penyebab penyakit, epidemiologi menggunakan konsep dasar yang dikenal sebagai Segitiga Epidemiologi (Epidemiologic Triad). Model ini menggambarkan bahwa penyakit muncul akibat adanya interaksi antara tiga komponen utama, yaitu Agent (penyebab penyakit), Host (inang/manusia), dan Environment (lingkungan) [8], [9].

Kondisi kesehatan sangat bergantung pada keseimbangan ketiga faktor tersebut. Jika keseimbangan ini terganggu misalnya agen penyakit menjadi lebih ganas, daya tahan tubuh menurun, atau kondisi lingkungan mendukung penularan maka penyakit akan terjadi [8].

2.1.1 Agen (Agent)

Agen merupakan penyebab utama yang keberadaannya mutlak diperlukan agar suatu penyakit bisa muncul. Secara umum, agen penyebab penyakit ini bisa dibagi menjadi beberapa kelompok [9], [10]:

  1. Biologis: Meliputi mikroorganisme seperti virus, bakteri, jamur, parasit, dan protozoa.
  2. Kimia: Berupa zat-zat seperti limbah beracun, obat-obatan, atau bahan kimia berbahaya lainnya.
  3. Fisik: Termasuk faktor fisik seperti radiasi, benturan atau trauma, suhu yang terlalu ekstrem, serta kebisingan.
  4. Nutrisi: Berkaitan dengan masalah gizi, baik itu kekurangan maupun kelebihan (contohnya defisiensi vitamin atau obesitas).

Selain jenisnya, karakteristik agen juga sangat menentukan, seperti infektivitas (kemampuannya untuk masuk dan berkembang biak dalam tubuh), patogenisitas (kemampuannya menyebabkan sakit), serta virulensi (seberapa parah penyakit yang ditimbulkannya) [9].

2.1.2 Penjamu (Host)

Penjamu atau host adalah sebutan untuk manusia atau organisme yang rentan terkena penyakit. Karakteristik yang dimiliki oleh penjamu ini sangat menentukan seberapa besar risiko mereka tertular serta tingkat keparahan penyakit yang akan dialami. Faktor-faktor tersebut meliputi [8], [10]:

  • Demografis: Mencakup aspek-aspek dasar kependudukan seperti usia, jenis kelamin, ras, dan etnis.
  • Biologis dan Genetik: Berkaitan dengan sistem kekebalan tubuh, kondisi fisik, serta faktor keturunan. Isu interaksi antara gen dan lingkungan (gene-environment interaction) juga menjadi bahasan penting dalam aspek ini [8].
  • Perilaku: Berhubungan dengan kebiasaan hidup sehari-hari, misalnya pola makan, aktivitas fisik, serta kebersihan diri.

2.1.3 Lingkungan (Environment)

Lingkungan merupakan faktor eksternal yang mempengaruhi kondisi agen maupun penjamu. Secara umum, lingkungan terbagi menjadi beberapa aspek [9], [10]:

  • Fisik: Meliputi iklim, letak geografis, kualitas udara, serta kondisi perumahan.
  • Biologis: Mencakup keberadaan flora, fauna, dan vektor pembawa penyakit seperti nyamuk.
  • Sosial-Ekonomi: Termasuk faktor pekerjaan, pendidikan, kepadatan penduduk, dan akses ke fasilitas kesehatan.

2.2 Ukuran Frekuensi

Dalam epidemiologi, langkah paling dasar adalah mengukur kejadian penyakit serta hubungannya dengan paparan tertentu. Bagian ini akan membahas tiga ukuran utama: Prevalensi sebagai ukuran frekuensi, serta Prevalence Ratio (PR) dan Prevalence Odds Ratio (POR) sebagai ukuran asosiasi, terutama dalam studi cross-sectional [8], [9].

Untuk mempermudah perhitungan rumus selanjutnya, data diasumsikan tersusun dalam tabel kontingensi \(2 \times 2\) sebagai berikut:

Sakit (Kasus) Tidak Sakit Total
Terpapar \(a\) \(b\) \(a+b\)
Tidak Terpapar \(c\) \(d\) \(c+d\)

2.2.1 Prevalensi (Prevalence)

Prevalensi menggambarkan proporsi individu dalam populasi yang sedang mengalami penyakit atau kondisi tertentu, baik pada satu titik waktu (point prevalence) maupun dalam rentang waktu tertentu (period prevalence). Berbeda dengan insidensi yang hanya mencatat kasus baru, prevalensi mencakup gabungan antara kasus baru dan kasus lama yang masih ada [8], [10].

Secara matematis, prevalensi (\(P\)) dapat dihitung menggunakan rumus:

\[ P = \frac{\text{Jumlah kasus yang ada (lama + baru)}}{\text{Total populasi berisiko}} \]

Atau berdasarkan notasi tabel di atas, prevalensi penyakit dalam total populasi studi adalah:

\[ P = \frac{a + c}{a + b + c + d} \]

Prevalensi dipengaruhi oleh dua faktor utama: laju insidensi (munculnya kasus baru) dan durasi penyakit. Jika penyakit bersifat kronis dan tidak mematikan dengan cepat, prevalensi akan cenderung tinggi [10].

2.2.2 Prevalensi Ratio (Prevalence Ratio)

Prevalence Ratio (PR) merupakan ukuran asosiasi yang membandingkan prevalensi penyakit pada kelompok terpapar dengan kelompok tidak terpapar. Ukuran ini lebih sering digunakan dalam studi cross-sectional ketika penyakit yang diteliti cukup umum (prevalensinya tinggi), karena memberikan gambaran risiko yang lebih mudah dipahami dibandingkan Odds Ratio [11].

Adapun rumus untuk menghitung Prevalence Ratio adalah:

\[ PR = \frac{\text{Prevalensi pada Kelompok Terpapar}}{\text{Prevalensi pada Kelompok Tidak Terpapar}} \]

Dalam notasi aljabar:

\[ PR = \frac{a / (a+b)}{c / (c+d)} \]

Interpretasi nilai PR: * \(PR = 1\): Tidak ada hubungan antara paparan dan penyakit. * \(PR > 1\): Paparan merupakan faktor risiko (meningkatkan prevalensi). * \(PR < 1\): Paparan bersifat protektif (menurunkan prevalensi).

2.2.3 Prevalensi Odds Ratio (Prevalence Odds Ratio)

Prevalence Odds Ratio (POR) adalah perbandingan odds penyakit antara kelompok yang terpapar dengan kelompok yang tidak terpapar. Walaupun perhitungan rumusnya sama persis dengan Odds Ratio (OR) dalam studi kasus-kontrol, POR secara khusus dihitung menggunakan data prevalensi yang diperoleh dari studi cross-sectional [8], [11].

Rumus untuk menghitung Prevalence Odds Ratio adalah: \[ POR = \frac{\text{Odds Sakit pada Terpapar}}{\text{Odds Sakit pada Tidak Terpapar}} = \frac{a/b}{c/d} \]

Disederhanakan menjadi rumus silang:

\[ POR = \frac{ad}{bc} \]

Penting untuk dicatat bahwa jika penyakit yang diteliti jarang terjadi (prevalensi < 10%), nilai POR akan mendekati nilai PR. Namun, jika prevalensi penyakit tinggi, POR cenderung memberikan angka yang lebih ekstrem (lebih jauh dari angka 1) dibandingkan PR, sehingga dapat menaksir risiko secara berlebihan (overestimation) [2], [4].

2.3 Desain Studi

Studi cross-sectional atau potong lintang merupakan desain penelitian observasional di mana pengukuran paparan (faktor risiko) dan penyakit (outcome) dilakukan secara bersamaan pada satu waktu. Karena karakteristiknya ini, studi tersebut sering diibaratkan sebagai “potret” (snapshot) kondisi populasi. Unit analisisnya umumnya adalah individu, dengan tujuan utama mengestimasi prevalensi penyakit serta menjajaki hubungan awal antar variabel [8], [9]. Karena tidak melibatkan pemantauan berkelanjutan (follow-up), ukuran frekuensi yang dihasilkan adalah prevalensi, sedangkan ukuran asosiasinya menggunakan Prevalence Ratio (PR) atau Prevalence Odds Ratio (POR) [11].

Keunggulan utama dari desain ini adalah efisiensi dari segi waktu dan biaya, serta sangat berguna untuk memetakan beban penyakit sebagai dasar perencanaan layanan kesehatan, misalnya pada survei nasional [8]. Namun, kelemahan signifikannya adalah ketidakmampuan memastikan urutan waktu (temporal relationship) antara sebab dan akibat. Selain itu, desain ini rentan terhadap bias prevalensi karena cenderung hanya menjaring kasus yang bertahan hidup lama (survivor bias) dan melewatkan kasus yang cepat sembuh atau fatal [9], [11].

Bab 3. Metodologi Penelitian

Penelitian ini menggunakan pendekatan kuantitatif dengan metode deskriptif dan analitik untuk mengkaji distribusi spasial dan faktor risiko kasus campak di Provinsi Jawa Barat pada tahun 2024. Analisis dilakukan terhadap data sekunder dari 27 kabupaten/kota.

3.1 Sumber Data

Penelitian ini menggunakan jenis data sekunder yang diperoleh dari portal resmi Open Data Jabar (Pemerintah Provinsi Jawa Barat). Data yang dikumpulkan mencakup data penyakit (kasus campak), data demografi (kependudukan), dan data faktor risiko terkait di level Kabupaten/Kota [12]-[16].

3.2 Variabel Penelitian

Variabel dalam penelitian ini diklasifikasikan berdasarkan peranannya dalam analisis asosiasi dan perhitungan frekuensi. Variabel Jumlah Penduduk ditempatkan sebagai Variabel Z (Variabel Tambahan/Denominator) yang berfungsi untuk menstandarisasi jumlah kasus menjadi ukuran frekuensi (Prevalensi/Insidensi) agar dapat dibandingkan antarwilayah.

Berikut adalah tabel operasional variabel:

Kategori Variabel Nama Variabel Definisi Operasional Peran dalam Analisis Skala Data
Variabel Dependen (Y) Jumlah Kasus Campak Jumlah kejadian kasus campak yang tercatat di fasilitas kesehatan pada periode waktu tertentu. Outcome (Akibat) Rasio (Diskrit)
Variabel Independen (X1) Cakupan Imunisasi Persentase target populasi yang menerima vaksinasi campak di suatu wilayah. Exposure (Paparan/Faktor Risiko) Rasio (Kontinu)
Variabel Independen (X2) Persentase Akses Sanitasi Layak Persentase rumah tangga yang memiliki akses terhadap sanitasi layak Exposure (Faktor Lingkungan) Rasio (Kontinu)
Variabel Independen (X3) Kepadatan Penduduk Jumlah penduduk dibagi dengan luas wilayah (\(jiwa/km^2\)). Exposure (Faktor Lingkungan) Rasio (Kontinu)
Variabel Tambahan (Z) Jumlah Penduduk Total populasi yang berisiko (population at risk) di wilayah tersebut. Denominator (Pembagi) untuk menghitung Rate/Proporsi Rasio (Diskrit)

3.3 Metode Analisis

Analisis data dilakukan secara bertahap menggunakan perangkat lunak statistik (R/RStudio). Tahapan analisis meliputi analisis deskriptif, perhitungan ukuran epidemiologis, hingga pemodelan ekonometrika spasial.

3.3.1 Statistika Deskriptif dan Visualisasi

Langkah pertama adalah melakukan eksplorasi data untuk memahami karakteristik awalnya. Tahapan ini mencakup tiga bagian utama:

  1. Statistika Deskriptif: Menghitung nilai rata-rata (mean), standar deviasi, serta nilai minimum dan maksimum. Tujuannya adalah untuk memberikan gambaran umum mengenai distribusi data pada variabel dependen maupun independen.
  2. Peta Sebaran (Choropleth Map): Membuat visualisasi peta untuk melihat sebaran prevalensi campak dan faktor risikonya di tingkat kabupaten/kota, sehingga pola spasialnya dapat terlihat dengan jelas.
  3. Visualisasi Tren: Menggunakan grafik garis untuk mengamati pola pergerakan jumlah kasus atau tingkat prevalensi antar wilayah.

3.3.2 Ukuran Epidemiologis

Pada tahap ini, dilakukan perhitungan statistik menggunakan data agregat yang mencakup:

  1. Ukuran Frekuensi: Menghitung nilai Prevalensi, yaitu proporsi jumlah kasus campak terhadap total populasi yang berisiko.
  2. Ukuran Asosiasi: Mengestimasi hubungan antar variabel menggunakan Prevalence Ratio (PR) dan Prevalence Odds Ratio (POR). Kedua ukuran ini digunakan untuk melihat seberapa besar peluang risiko paparan terhadap kejadian penyakit pada data cross-sectional.

3.3.3 Pemodelan Spasial

Analisis ini dilakukan untuk menangkap efek ketergantungan antarwilayah atau spatial dependence, hal yang sering kali luput dalam analisis regresi biasa.

  1. Uji Autokorelasi Spasial Global (Global Moran’s I) Langkah ini dipakai untuk melihat apakah ada autokorelasi spasial pada variabel terikat di seluruh wilayah Jawa Barat secara umum. Tujuannya untuk menguji hipotesis apakah penyebaran kasusnya bersifat acak (random), mengelompok (clustered), atau justru menyebar (dispersed).

  2. Uji Autokorelasi Spasial Lokal (LISA) Local Indicators of Spatial Association atau LISA digunakan untuk melacak lokasi-lokasi spesifik yang punya hubungan spasial kuat. Hasilnya nanti dibagi menjadi empat kategori:

    • High-High (Hotspots): Wilayah dengan kasus tinggi dikelilingi wilayah yang juga tinggi.
    • Low-Low (Coldspots): Wilayah dengan kasus rendah dikelilingi wilayah yang rendah.
    • High-Low dan Low-High: Wilayah yang menjadi pencilan spasial (spatial outliers).

  3. Estimasi Model Regresi Ada tiga model yang akan dibangun dan dibandingkan dalam penelitian ini:

    • Ordinary Least Squares (OLS): Ini adalah regresi linier biasa yang belum memperhitungkan efek spasial, dipakai sebagai pembanding (baseline).
    • Spatial Autoregressive Model (SAR): Model ini berasumsi kalau ketergantungan spasial ada di variabel dependennya (spatial lag). Simpelnya, kejadian penyakit di satu daerah dipengaruhi oleh kejadian di daerah tetangganya.
    • Spatial Error Model (SEM): Kalau model ini, asumsinya ketergantungan spasial ada di error atau sisaannya. Biasanya terjadi karena ada faktor lain yang tidak kita teliti tapi polanya mirip antarwilayah.

  4. Pemilihan Model Terbaik Untuk menentukan mana model yang paling pas, dilakukan beberapa tes diagnostik:

    • Uji Lagrange Multiplier (LM): Menggunakan uji LM-lag dan LM-error (termasuk versi Robust-nya) untuk melihat apakah model spasial memang lebih bagus dari OLS, serta memberi petunjuk apakah lebih cocok pakai SAR atau SEM.
    • Akaike Information Criterion (AIC): Kita akan memilih model berdasarkan nilai AIC. Model yang dianggap terbaik (best fit) adalah yang memiliki nilai AIC terendah dengan tetap mempertimbangkan nilai \(R^2\) yang baik.

3.4 Alur Penelitian

Secara garis besar, alur pengerjaan penelitian ini dibagi menjadi beberapa tahapan:

  1. Identifikasi Masalah: Menentukan fokus analisis dan variabel yang akan digunakan terkait kasus campak.
  2. Pengumpulan Data: Mengambil data sekunder (data kasus, populasi, dan variabel prediktor) dari portal Open Data Jabar.
  3. Pengolahan Data (Data Preprocessing): Menyiapkan data agar siap diolah, meliputi pembersihan data (cleaning), sinkronisasi kode wilayah, dan penggabungan dataset (merging).
  4. Analisis Deskriptif: Melakukan eksplorasi data awal dan visualisasi untuk melihat gambaran umum data.
  5. Perhitungan Ukuran Asosiasi: Menghitung nilai-nilai dasar seperti Prevalensi, Prevalence Ratio (PR), dan Prevalence Odds Ratio (POR).
  6. Uji Asumsi Spasial: Mengecek apakah ada ketergantungan antar-lokasi menggunakan uji Moran’s I (Global) dan melihat pola penyebarannya dengan LISA (Lokal).
  7. Pemodelan Statistik: Membangun model regresi mulai dari metode standar (OLS) hingga model spasial (SAR dan SEM).
  8. Evaluasi dan Pemilihan Model: Menentukan model terbaik (best fit) menggunakan uji Lagrange Multiplier (LM) dan kriteria AIC, dilanjutkan dengan interpretasi hasil.

Bab 4. Hasil dan Pembahasan

4.1 Statistika Deskriptif dan Visualisasi

4.1.1 Tabel Hasil Data

library(sf)
library(dplyr)
library(stringr)
library(tmap)
library(spdep)
library(sp)
library(stars)
library(gstat)
library(ggplot2)

# 1. Menyiapkan Data Penduduk (Denominator)
data <- read.csv("C:/Users/Admin/Downloads/KULIAH/Semester 5/Epidem/UAS & Project (Rshiny)/File Bahan Dashboard/(CSV) Data UAS & Project Epidem 2025.csv")
jabar_shp <- st_read("C:/Users/Admin/Downloads/KULIAH/Semester 5/Epidem/UAS & Project (Rshiny)/File Bahan Dashboard/RBI_50K_2023_Jawa Barat.shp")
## Reading layer `RBI_50K_2023_Jawa Barat' from data source 
##   `C:\Users\Admin\Downloads\KULIAH\Semester 5\Epidem\UAS & Project (Rshiny)\File Bahan Dashboard\RBI_50K_2023_Jawa Barat.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3703 ymin: -7.82099 xmax: 108.8468 ymax: -5.806538
## Geodetic CRS:  WGS 84
penduduk_df <- data.frame(
  nama_kabupaten_kota = c(
    "KABUPATEN BOGOR","KABUPATEN SUKABUMI","KABUPATEN CIANJUR",
    "KABUPATEN BANDUNG","KABUPATEN GARUT","KABUPATEN TASIKMALAYA",
    "KABUPATEN CIAMIS","KABUPATEN KUNINGAN","KABUPATEN CIREBON",
    "KABUPATEN MAJALENGKA","KABUPATEN SUMEDANG","KABUPATEN INDRAMAYU",
    "KABUPATEN SUBANG","KABUPATEN PURWAKARTA","KABUPATEN KARAWANG",
    "KABUPATEN BEKASI","KABUPATEN BANDUNG BARAT","KABUPATEN PANGANDARAN",
    "KOTA BOGOR","KOTA SUKABUMI","KOTA BANDUNG",
    "KOTA CIREBON","KOTA BEKASI","KOTA DEPOK",
    "KOTA CIMAHI","KOTA TASIKMALAYA","KOTA BANJAR"
  ),
  jumlah_penduduk = c(
    56823,28282,25849,37532,27169,19202,12593,12133,23876,
    13524,11873,19144,16634,10504,25548,32739,18841,
    4341,10785,3657,25286,3447,26445,21637,5982,7501,2093
  )
)

# 2. Cleaning Shapefile Jawa Barat
jabar_shp <- jabar_shp %>%
  mutate(
    KabKota = case_when(
      grepl("^Kota ", WADMKK) ~ gsub("^Kota ", "", WADMKK),
      grepl("^Kabupaten ", WADMKK) ~ gsub("^Kabupaten ", "", WADMKK),
      TRUE ~ WADMKK
    ),
    KabKota = trimws(KabKota)
  )

# 3. Filter Data Tahun 2024 dan Cleaning Nama Wilayah
data_2024 <- data %>%
  filter(tahun == 2024) %>%
  mutate(
    WADMKK = case_when(
      grepl("^KOTA", nama_kabupaten_kota) ~ 
        str_to_title(sub("^KOTA ", "Kota ", nama_kabupaten_kota)),
      grepl("^KABUPATEN", nama_kabupaten_kota) ~ 
        str_to_title(sub("^KABUPATEN ", "", nama_kabupaten_kota)),
      TRUE ~ str_to_title(nama_kabupaten_kota)
    ),
    WADMKK = str_trim(WADMKK)
  )

# 4. Join Data Spasial dan Tabular
jabar_2024_sf <- jabar_shp %>%
  left_join(data_2024, by = "WADMKK") %>%
  rename(
    jumlah_kasus = jumlah_kasus..Y.,
    imunisasi    = jumlah_bayi_imunisasi..X1.,
    sanitasi     = persentase_sanitasi_layak..X2.,
    kepadatan    = kepadatan_penduduk..X3.
  ) %>%
  left_join(penduduk_df, by = "nama_kabupaten_kota")

# 5. Menghitung Prevalensi per 100.000 Penduduk
jabar_2024_sf <- jabar_2024_sf %>%
  mutate(
    prevalensi = jumlah_kasus / jumlah_penduduk,
    prevalensi_100k = prevalensi * 100000
  )

# Menampilkan Tabel 5 Wilayah dengan Prevalensi Tertinggi
prev_table <- jabar_2024_sf %>%
  st_drop_geometry() %>%
  select(nama_kabupaten_kota, jumlah_kasus, jumlah_penduduk, prevalensi_100k) %>%
  arrange(desc(prevalensi_100k)) %>%
  head(5)

print(prev_table)
##   nama_kabupaten_kota jumlah_kasus jumlah_penduduk prevalensi_100k
## 1        KOTA CIREBON          107            3447        3104.149
## 2    KOTA TASIKMALAYA          152            7501        2026.396
## 3       KOTA SUKABUMI           66            3657        1804.758
## 4          KOTA BOGOR          194           10785        1798.795
## 5          KOTA DEPOK          388           21637        1793.225

Berdasarkan tabel, terlihat jelas bahwa jumlah kasus absolut yang besar belum tentu menunjukkan tingkat risiko yang paling parah. Data menunjukkan Kota Cirebon justru memiliki rasio kasus per penduduk (prevalensi) tertinggi meskipun jumlah pasiennya sedikit, sangat kontras dengan Kota Depok yang memiliki jumlah pasien terbanyak namun rasio risikonya paling rendah di kelompok lima besar ini.

4.1.2 Statistika Deskriptif Kasus Campak

library(dplyr)
library(sf) # Diperlukan karena data Anda formatnya spasial (sf)

# Membuat tabel ringkasan
stat_deskriptif <- jabar_2024_sf %>%
  st_drop_geometry() %>% # Hilangkan geometri peta agar proses hitung fokus ke angka
  summarise(
    Total_Kasus = sum(jumlah_kasus),
    Rata_rata = mean(jumlah_kasus),
    Median = median(jumlah_kasus),
    Std_Dev = sd(jumlah_kasus),
    Min = min(jumlah_kasus),
    Max = max(jumlah_kasus),
    Rentang_Data = max(jumlah_kasus) - min(jumlah_kasus)
  )

# Menampilkan hasil
print(stat_deskriptif)
##   Total_Kasus Rata_rata Median  Std_Dev Min Max Rentang_Data
## 1        3682  136.3704    107 99.94081   0 413          413

Berdasarkan tabel deskriptif, terlihat adanya ketimpangan yang cukup mencolok pada pola penyebaran data. Nilai rata-ratanya (136) tercatat lebih tinggi dibandingkan nilai tengah (107). Ini menunjukkan kalau distribusinya tidak normal, kemungkinan besar karena adanya wilayah dengan kasus yang sangat ekstrem hingga 413 kasus (outlier). Variasi datanya juga sangat luas, terlihat dari standar deviasi yang besar (hampir 100), yang artinya perbedaan jumlah kasus antar satu wilayah dengan wilayah lainnya sangat tajam.

4.1.3 Visualisasi Peta Sebaran (Choropleth Map)

# Ubah mode ke "plot" agar aman untuk PDF/Word
tmap_mode("plot") 
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(jabar_2024_sf) +
  tm_fill(
    col = "jumlah_kasus",          # Variabel yang akan dipetakan
    style = "quantile",            # Metode klasifikasi (quantile)
    n = 5,                         # Jumlah kelas
    palette = "Reds",              # Skala warna merah
    title = "Jumlah Kasus (2024)", # Judul Legenda
    popup.vars = c(                # Info yang muncul saat diklik (hanya di mode view)
      "Kab/Kota"   = "WADMKK",
      "Kasus"      = "jumlah_kasus",
      "Sanitasi %" = "sanitasi",
      "Kepadatan"  = "kepadatan",
      "Prev per 100k" = "prevalensi_100k"
    )
  ) +
  tm_borders(lwd = 0.5) +          # Garis batas wilayah
  tm_layout(
    title = "Peta Kasus Penyakit Jawa Barat Tahun 2024",
    title.position = c("center", "top"),
    legend.position = c("right", "bottom"),
    frame = TRUE
  ) +
  tm_compass(position = c("left", "top")) +
  tm_scale_bar(position = c("left", "bottom"))
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

Berdasarkan peta, distribusi data terlihat sangat tidak merata (heterogen). Terdapat penumpukan nilai ekstrem di wilayah utara (merah pekat) dengan kasus di atas 193, yang sangat kontras dibandingkan wilayah lainnya. Ketimpangan ini menunjukkan bahwa datanya memiliki variasi spasial yang tinggi dan cenderung mengelompok di area tertentu saja.

4.1.4 Tren Kejadian Penyakit

# Persiapan Data Tren
data_tren <- data %>%
  rename(
    jumlah_kasus = jumlah_kasus..Y.,
    imunisasi    = jumlah_bayi_imunisasi..X1.,
    sanitasi     = persentase_sanitasi_layak..X2.,
    kepadatan    = kepadatan_penduduk..X3.
  ) %>%
  mutate(
    tahun = as.integer(tahun),
    jumlah_kasus = as.numeric(jumlah_kasus),
    imunisasi = as.numeric(imunisasi),
    sanitasi = as.numeric(sanitasi),
    kepadatan = as.numeric(kepadatan)
  )

# Agregasi Total Kasus per Tahun
tren_global <- data_tren %>%
  group_by(tahun) %>%
  summarise(total_kasus = sum(jumlah_kasus, na.rm = TRUE))

# Plotting Grafik Garis
ggplot(tren_global, aes(x = tahun, y = total_kasus)) +
  geom_line(linewidth = 1.2, color = "red") +
  geom_point(size = 2) +
  labs(
    title = "Tren Total Kasus Penyakit Jawa Barat",
    x = "Tahun",
    y = "Total Kasus"
  ) +
  theme_minimal()

Berdasarkan grafik, data menunjukkan pola fluktuasi yang ekstrem, di mana terjadi lonjakan nilai yang sangat drastis pada tahun 2023 (mencapai puncak) setelah sebelumnya cenderung stabil dan rendah. Meskipun tren di tahun 2024 sudah mengalami penurunan tajam, posisi datanya masih berada di atas rata-rata periode awal (2018-2021), yang menandakan bahwa kondisi data saat ini belum sepenuhnya kembali stabil ke level normal.

4.2 Perhitungan Ukuran Asosiasi dan Frekuensi

4.2.1 Ukuran Frekuensi (Prevalensi)

# Menggunakan data 'jabar_2024_sf' yang sudah dihitung pada subbab 4.1

# Menampilkan ringkasan statistik Prevalensi
summary_prev <- summary(jabar_2024_sf$prevalensi_100k)

# Menampilkan wilayah dengan prevalensi terendah dan tertinggi
min_prev <- jabar_2024_sf %>% 
  arrange(prevalensi_100k) %>% 
  slice(1) %>% 
  select(WADMKK, prevalensi_100k) %>% 
  st_drop_geometry()

max_prev <- jabar_2024_sf %>% 
  arrange(desc(prevalensi_100k)) %>% 
  slice(1) %>% 
  select(WADMKK, prevalensi_100k) %>% 
  st_drop_geometry()

print("Ringkasan Statistik Prevalensi per 100.000 Penduduk:")
## [1] "Ringkasan Statistik Prevalensi per 100.000 Penduduk:"
print(summary_prev)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   534.1   678.1   954.6  1311.7  3104.1
print(paste("Prevalensi Terendah:", min_prev$WADMKK, "(", round(min_prev$prevalensi_100k, 2), ")"))
## [1] "Prevalensi Terendah: Sukabumi ( 0 )"
print(paste("Prevalensi Tertinggi:", max_prev$WADMKK, "(", round(max_prev$prevalensi_100k, 2), ")"))
## [1] "Prevalensi Tertinggi: Kota Cirebon ( 3104.15 )"

Berdasarkan ringkasan statistik, data menunjukkan ketimpangan yang cukup ekstrem karena nilai rata-rata (954,6) posisinya jauh lebih tinggi dibandingkan nilai tengah (678,1). Kondisi ini menandakan bahwa sebaran data “tertarik” ke atas oleh nilai ekstrem di Kota Cirebon (3.104) yang angkanya melompat sangat jauh meninggalkan mayoritas wilayah lain, sehingga menciptakan rentang data yang sangat lebar dari 0 hingga 3.000-an.

4.2.2 Ukuran Asosiasi (Faktor Cakupan Imunisasi)

# 1. Menentukan Cut-off (Rata-rata)
mean_kasus <- mean(jabar_2024_sf$jumlah_kasus, na.rm = TRUE)
mean_imun <- mean(jabar_2024_sf$imunisasi, na.rm = TRUE)

# 2. Membentuk Tabel Kontingensi 2x2
# a: Exposed (Imun Rendah) & Sakit (Kasus Tinggi)
a <- nrow(filter(jabar_2024_sf, imunisasi < mean_imun & jumlah_kasus >= mean_kasus))

# b: Exposed (Imun Rendah) & Tidak Sakit (Kasus Rendah)
b <- nrow(filter(jabar_2024_sf, imunisasi < mean_imun & jumlah_kasus < mean_kasus))

# c: Unexposed (Imun Tinggi) & Sakit (Kasus Tinggi)
c <- nrow(filter(jabar_2024_sf, imunisasi >= mean_imun & jumlah_kasus >= mean_kasus))

# d: Unexposed (Imun Tinggi) & Tidak Sakit (Kasus Rendah)
d <- nrow(filter(jabar_2024_sf, imunisasi >= mean_imun & jumlah_kasus < mean_kasus))

# Membuat Matrix Tampilan
tabel_2x2_imun <- matrix(c(a, c, b, d), nrow = 2, byrow = TRUE)
colnames(tabel_2x2_imun) <- c("Kasus Tinggi", "Kasus Rendah")
rownames(tabel_2x2_imun) <- c("Imunisasi Rendah (Risk)", "Imunisasi Tinggi (Ref)")

print("Tabel Kontingensi 2x2 (Imunisasi vs Kasus):")
## [1] "Tabel Kontingensi 2x2 (Imunisasi vs Kasus):"
print(tabel_2x2_imun)
##                         Kasus Tinggi Kasus Rendah
## Imunisasi Rendah (Risk)            4            9
## Imunisasi Tinggi (Ref)            12            2
# 3. Menghitung PR dan POR
# Prevalence Ratio
PR_imun <- (a / (a + b)) / (c / (c + d))

# Prevalence Odds Ratio
POR_imun <- (a * d) / (b * c)

print(paste("Prevalence Ratio (PR) Imunisasi:", round(PR_imun, 3)))
## [1] "Prevalence Ratio (PR) Imunisasi: 0.306"
print(paste("Prevalence Odds Ratio (POR) Imunisasi:", round(POR_imun, 3)))
## [1] "Prevalence Odds Ratio (POR) Imunisasi: 0.074"

Berdasarkan hasil analisis statistik, data memperlihatkan hubungan yang terbalik (anomali) dari dugaan umum. Wilayah dengan status Imunisasi Rendah justru tercatat memiliki risiko yang jauh lebih kecil (nilai rasio 0,3) untuk mengalami kejadian “Kasus Tinggi” dibandingkan wilayah dengan Imunisasi Tinggi. Secara statistik, data ini menunjukkan bahwa kejadian luar biasa (kasus tinggi) justru menumpuk di kelompok wilayah yang cakupan imunisasinya baik. Kemungkinan besar karena wilayah tersebut adalah kota padat penduduk yang memudahkan penularan sementara wilayah dengan imunisasi rendah justru didominasi oleh kasus yang sedikit.

4.2.3 Ukuran Asosiasi (Faktor Sanitasi Layak)

# 1. Menentukan Cut-off
mean_sanitasi <- mean(jabar_2024_sf$sanitasi, na.rm = TRUE)

# 2. Membentuk Tabel Kontingensi 2x2
# a: Exposed (Sanitasi Buruk) & Sakit (Kasus Tinggi)
a_san <- nrow(filter(jabar_2024_sf, sanitasi < mean_sanitasi & jumlah_kasus >= mean_kasus))

# b: Exposed (Sanitasi Buruk) & Tidak Sakit (Kasus Rendah)
b_san <- nrow(filter(jabar_2024_sf, sanitasi < mean_sanitasi & jumlah_kasus < mean_kasus))

# c: Unexposed (Sanitasi Baik) & Sakit (Kasus Tinggi)
c_san <- nrow(filter(jabar_2024_sf, sanitasi >= mean_sanitasi & jumlah_kasus >= mean_kasus))

# d: Unexposed (Sanitasi Baik) & Tidak Sakit (Kasus Rendah)
d_san <- nrow(filter(jabar_2024_sf, sanitasi >= mean_sanitasi & jumlah_kasus < mean_kasus))

# Membuat Matrix Tampilan
tabel_2x2_san <- matrix(c(a_san, c_san, b_san, d_san), nrow = 2, byrow = TRUE)
colnames(tabel_2x2_san) <- c("Kasus Tinggi", "Kasus Rendah")
rownames(tabel_2x2_san) <- c("Sanitasi Buruk (Risk)", "Sanitasi Baik (Ref)")

print("Tabel Kontingensi 2x2 (Sanitasi vs Kasus):")
## [1] "Tabel Kontingensi 2x2 (Sanitasi vs Kasus):"
print(tabel_2x2_san)
##                       Kasus Tinggi Kasus Rendah
## Sanitasi Buruk (Risk)            8            5
## Sanitasi Baik (Ref)              7            7
# 3. Menghitung PR dan POR
# Prevalence Ratio
PR_san <- (a_san / (a_san + b_san)) / (c_san / (c_san + d_san))

# Prevalence Odds Ratio
POR_san <- (a_san * d_san) / (b_san * c_san)

print(paste("Prevalence Ratio (PR) Sanitasi:", round(PR_san, 3)))
## [1] "Prevalence Ratio (PR) Sanitasi: 1.28"
print(paste("Prevalence Odds Ratio (POR) Sanitasi:", round(POR_san, 3)))
## [1] "Prevalence Odds Ratio (POR) Sanitasi: 1.6"

Berdasarkan hasil analisis statistik, data menunjukkan pola hubungan yang positif (sejalan) di mana faktor lingkungan berperan nyata. Wilayah dengan status Sanitasi Buruk tercatat memiliki risiko 1,28 kali lebih besar untuk mengalami kejadian “Kasus Tinggi” dibandingkan wilayah dengan sanitasi baik. Angka ini mengonfirmasi secara data bahwa buruknya kondisi kebersihan berbanding lurus dengan peningkatan peluang terjadinya lonjakan kasus di suatu wilayah.

4.2.4 Ukuran Asosiasi (Faktor Kepadatan Penduduk)

# 1. Menentukan Cut-off
mean_kepadatan <- mean(jabar_2024_sf$kepadatan, na.rm = TRUE)

# 2. Membentuk Tabel Kontingensi 2x2
# a: Exposed (Padat) & Sakit (Kasus Tinggi)
a_kep <- nrow(filter(jabar_2024_sf, kepadatan >= mean_kepadatan & jumlah_kasus >= mean_kasus))

# b: Exposed (Padat) & Tidak Sakit (Kasus Rendah)
b_kep <- nrow(filter(jabar_2024_sf, kepadatan >= mean_kepadatan & jumlah_kasus < mean_kasus))

# c: Unexposed (Jarang) & Sakit (Kasus Tinggi)
c_kep <- nrow(filter(jabar_2024_sf, kepadatan < mean_kepadatan & jumlah_kasus >= mean_kasus))

# d: Unexposed (Jarang) & Tidak Sakit (Kasus Rendah)
d_kep <- nrow(filter(jabar_2024_sf, kepadatan < mean_kepadatan & jumlah_kasus < mean_kasus))

# Membuat Matrix Tampilan
tabel_2x2_kep <- matrix(c(a_kep, c_kep, b_kep, d_kep), nrow = 2, byrow = TRUE)
colnames(tabel_2x2_kep) <- c("Kasus Tinggi", "Kasus Rendah")
rownames(tabel_2x2_kep) <- c("Kepadatan Tinggi (Risk)", "Kepadatan Rendah (Ref)")

print("Tabel Kontingensi 2x2 (Kepadatan vs Kasus):")
## [1] "Tabel Kontingensi 2x2 (Kepadatan vs Kasus):"
print(tabel_2x2_kep)
##                         Kasus Tinggi Kasus Rendah
## Kepadatan Tinggi (Risk)            5            8
## Kepadatan Rendah (Ref)             3           11
# 3. Menghitung PR dan POR
# Prevalence Ratio
PR_kep <- (a_kep / (a_kep + b_kep)) / (c_kep / (c_kep + d_kep))

# Prevalence Odds Ratio
POR_kep <- (a_kep * d_kep) / (b_kep * c_kep)

print(paste("Prevalence Ratio (PR) Kepadatan:", round(PR_kep, 3)))
## [1] "Prevalence Ratio (PR) Kepadatan: 1.484"
print(paste("Prevalence Odds Ratio (POR) Kepadatan:", round(POR_kep, 3)))
## [1] "Prevalence Odds Ratio (POR) Kepadatan: 2.292"

Berdasarkan hasil analisis risiko, data menunjukkan pola hubungan yang sejalan di mana kepadatan penduduk menjadi faktor pendorong kenaikan kasus. Wilayah dengan status kepadatan tinggi tercatat memiliki risiko hampir 1,5 kali lebih besar untuk mengalami kejadian kasus tinggi dibandingkan wilayah yang lebih lengang. Secara statistik, angka ini mengonfirmasi bahwa semakin padat populasi di suatu area, semakin tinggi pula probabilitas data kasusnya untuk melonjak.

4.3 Pemodelan Spasial

4.3.1 Matriks Pembobot dan Uji Autokorelasi Global (Moran’s I)

library(spdep)
library(spatialreg) # Diperlukan untuk model SAR/SEM

# 1. Membuat Matriks Pembobot Spasial (Queen Contiguity)
# Menggunakan data 'jabar_2024_sf' dari subbab sebelumnya
nb <- poly2nb(jabar_2024_sf, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# 2. Uji Global Moran's I pada Variabel Kasus
moran_result <- moran.test(jabar_2024_sf$jumlah_kasus, lw, randomisation = FALSE)

print("Hasil Uji Global Moran's I (Variabel Kasus):")
## [1] "Hasil Uji Global Moran's I (Variabel Kasus):"
print(moran_result)
## 
##  Moran I test under normality
## 
## data:  jabar_2024_sf$jumlah_kasus  
## weights: lw    
## 
## Moran I statistic standard deviate = 1.0192, p-value = 0.1541
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.10366301       -0.03846154        0.01944673

Berdasarkan hasil uji statistik, data menghasilkan angka indeks yang kecil yaitu 0,10 dengan nilai peluang (p-value) sebesar 0,15. Karena nilai peluang ini melebihi batas standar 0,05, maka pola pengelompokan wilayah yang terlihat sebelumnya dianggap tidak terbukti secara nyata. Kesimpulannya, tingginya kasus di suatu tempat secara statistik terjadi secara acak saja dan bukan disebabkan oleh pengaruh penularan yang kuat dari wilayah tetangganya.

4.3.2 Uji Autokorelasi Lokal (LISA)

# 1. Menghitung Local Moran's I
local_moran <- localmoran(jabar_2024_sf$jumlah_kasus, lw)
jabar_2024_sf$local_moran_pval <- local_moran[, 5] # Ambil p-value

# 2. Standarisasi Data (Z-score) untuk menentukan kuadran
# Centering variable (x - mean)
x <- scale(jabar_2024_sf$jumlah_kasus) %>% as.vector()
# Centering lag spatial (lag - mean)
y <- scale(lag.listw(lw, jabar_2024_sf$jumlah_kasus)) %>% as.vector()

# 3. Klasifikasi Kuadran LISA
# Buat kolom baru, default diisi NA
jabar_2024_sf$quadrant <- NA

# Tentukan Significance Level
sig_level <- 0.05

# Logika Pengisian (Hanya yang signifikan yang diberi label)
# High-High (Hotspot)
jabar_2024_sf$quadrant[x > 0 & y > 0 & jabar_2024_sf$local_moran_pval <= sig_level] <- "High-High"
# Low-Low (Coldspot)
jabar_2024_sf$quadrant[x < 0 & y < 0 & jabar_2024_sf$local_moran_pval <= sig_level] <- "Low-Low"
# Low-High (Outlier)
jabar_2024_sf$quadrant[x < 0 & y > 0 & jabar_2024_sf$local_moran_pval <= sig_level] <- "Low-High"
# High-Low (Outlier)
jabar_2024_sf$quadrant[x > 0 & y < 0 & jabar_2024_sf$local_moran_pval <= sig_level] <- "High-Low"

# Ubah menjadi Factor agar urutan legenda rapi
jabar_2024_sf$quadrant <- factor(jabar_2024_sf$quadrant, 
                                 levels = c("High-High", "High-Low", "Low-High", "Low-Low"))

# 4. Visualisasi Peta LISA
tmap_mode("plot") 
## ℹ tmap modes "plot" - "view"
tm_shape(jabar_2024_sf) +
  tm_fill(
    col = "quadrant", 
    # Warna sesuai urutan levels factor di atas: HH(Merah), HL(Pink), LH(Biru Muda), LL(Biru)
    palette = c("red", "lightpink", "lightblue", "blue"), 
    style = "cat",
    title = "Kategori LISA",
    # Bagian ini yang mengatasi double label:
    colorNA = "white",           # Warna untuk yang tidak signifikan (NA)
    textNA = "Tidak Signifikan"  # Label untuk yang tidak signifikan (NA)
  ) +
  tm_borders(lwd = 0.5) +
  tm_layout(
    title = "Peta Cluster LISA (Kasus Campak)",
    title.position = c("center", "top"),
    legend.position = c("left", "bottom"),
    frame = TRUE
  ) +
  tm_compass(position = c("right", "top")) +
  tm_scale_bar(position = c("right", "bottom"))
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "cat"`, use fill.scale =
## `tm_scale_categorical()`.
## ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'colorNA' (rename to
##   'value.na'), 'textNA' (rename to 'label.na') to
##   'tm_scale_categorical(<HERE>)'
## For small multiples, specify a 'tm_scale_' for each multiple, and put them in a
## list: 'fill'.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.

Berdasarkan peta, hasil analisis lokal ini mempertegas bahwa penyebaran kasus di Jawa Barat mayoritas bersifat acak karena hampir seluruh wilayah berwarna putih atau tidak memiliki hubungan spasial yang signifikan. Hanya terdapat satu area berwarna merah (High-High) di wilayah utara yang terdeteksi sebagai pusat pengelompokan (hotspot) di mana wilayah dengan kasus tinggi saling berdekatan, serta satu titik biru muda (Low-High) sebagai anomali data, sehingga secara statistik pengaruh antar-tetangga ini sifatnya sangat lokal dan tidak terjadi secara merata di seluruh provinsi.

4.3.3 Estimasi Model OLS dan Seleksi Model (Uji LM)

# 1. Model OLS (Baseline)
# Model: Kasus ~ Imunisasi + Sanitasi + Kepadatan
model_ols <- lm(jumlah_kasus ~ imunisasi + sanitasi + kepadatan, data = jabar_2024_sf)

print("Ringkasan Model OLS:")
## [1] "Ringkasan Model OLS:"
summary(model_ols)
## 
## Call:
## lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan, 
##     data = jabar_2024_sf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -137.03  -48.76  -11.84   41.15  212.82 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 99.2113253 44.9953728   2.205  0.03773 * 
## imunisasi    0.0027138  0.0008117   3.343  0.00282 **
## sanitasi    -0.9983793  0.5485434  -1.820  0.08179 . 
## kepadatan    0.0065807  0.0034855   1.888  0.07170 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 81.65 on 23 degrees of freedom
## Multiple R-squared:  0.4095, Adjusted R-squared:  0.3325 
## F-statistic: 5.318 on 3 and 23 DF,  p-value: 0.006232
# 2. Uji Lagrange Multiplier (LM Test) untuk Pemilihan Model
# Menguji apakah error OLS memiliki ketergantungan spasial
lm_test <- lm.LMtests(model_ols, lw, test = "all")
## Please update scripts to use lm.RStests in place of lm.LMtests
print("Hasil Uji Lagrange Multiplier (Diagnostik Spasial):")
## [1] "Hasil Uji Lagrange Multiplier (Diagnostik Spasial):"
print(lm_test)
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan,
## data = jabar_2024_sf)
## test weights: listw
## 
## RSerr = 0.69395, df = 1, p-value = 0.4048
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan,
## data = jabar_2024_sf)
## test weights: listw
## 
## RSlag = 0.00018578, df = 1, p-value = 0.9891
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan,
## data = jabar_2024_sf)
## test weights: listw
## 
## adjRSerr = 2.5706, df = 1, p-value = 0.1089
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan,
## data = jabar_2024_sf)
## test weights: listw
## 
## adjRSlag = 1.8769, df = 1, p-value = 0.1707
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan,
## data = jabar_2024_sf)
## test weights: listw
## 
## SARMA = 2.5708, df = 2, p-value = 0.2765

Berdasarkan hasil pemodelan OLS, variabel-variabel yang digunakan secara bersama-sama mampu menjelaskan pola naik-turunnya kasus campak sebesar 41% (R-squared 0,409), yang berarti model ini cukup layak digunakan. Secara statistik, variabel Imunisasi justru menunjukkan anomali karena berhubungan positif (kasus tetap tinggi meski imunisasi tinggi), sedangkan variabel Sanitasi dan Kepadatan Penduduk memiliki arah hubungan yang sesuai logika, di mana perbaikan sanitasi menurunkan kasus dan kepadatan penduduk menaikkan risiko kasus.

Lalu berdasarkan hasil uji diagnostik spasial, seluruh metode pengujian (Lagrange Multiplier) menunjukkan nilai peluang (p-value) yang jauh di atas 0,05. Hasil ini menyimpulkan bahwa secara statistik tidak ditemukan bukti adanya pengaruh hubungan antar-wilayah atau efek spasial dalam data ini. Dengan demikian, model regresi standar (OLS) yang sudah dibuat sebelumnya dinyatakan sudah valid dan cukup untuk menjelaskan data tanpa perlu menggunakan model spasial yang lebih rumit.

4.3.4 Estimasi Model Spasial (SAR & SEM) dan Perbandingan AIC

# 1. Estimasi Spatial Autoregressive Model (SAR/Lag Model)
model_sar <- lagsarlm(jumlah_kasus ~ imunisasi + sanitasi + kepadatan, 
                      data = jabar_2024_sf, listw = lw)

# 2. Estimasi Spatial Error Model (SEM/Error Model)
model_sem <- errorsarlm(jumlah_kasus ~ imunisasi + sanitasi + kepadatan, 
                        data = jabar_2024_sf, listw = lw)

# Menampilkan Ringkasan Singkat
print("--- Ringkasan Model SAR ---")
## [1] "--- Ringkasan Model SAR ---"
summary(model_sar)
## 
## Call:lagsarlm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan, 
##     data = jabar_2024_sf, listw = lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -136.984  -48.944  -11.850   41.165  212.749 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) 99.49870434 45.54872872  2.1844 0.0289295
## imunisasi    0.00271426  0.00075354  3.6020 0.0003158
## sanitasi    -0.99867453  0.51301463 -1.9467 0.0515733
## kepadatan    0.00659381  0.00349362  1.8874 0.0591082
## 
## Rho: -0.0023258, LR test value: 0.000151, p-value: 0.9902
## Asymptotic standard error: 0.21
##     z-value: -0.011075, p-value: 0.99116
## Wald statistic: 0.00012266, p-value: 0.99116
## 
## Log likelihood: -155.0128 for lag model
## ML residual variance (sigma squared): 5679.1, (sigma: 75.36)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 322.03, (AIC for lm: 320.03)
## LM test for residual autocorrelation
## test value: 2.5556, p-value: 0.1099
print("--- Ringkasan Model SEM ---")
## [1] "--- Ringkasan Model SEM ---"
summary(model_sem)
## 
## Call:errorsarlm(formula = jumlah_kasus ~ imunisasi + sanitasi + kepadatan, 
##     data = jabar_2024_sf, listw = lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -134.903  -53.427   10.014   40.287  207.582 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  1.0447e+02  3.6867e+01  2.8336 0.0046024
## imunisasi    2.7647e-03  7.2705e-04  3.8026 0.0001432
## sanitasi    -1.1167e+00  4.8245e-01 -2.3146 0.0206328
## kepadatan    6.5425e-03  3.0221e-03  2.1649 0.0303991
## 
## Lambda: -0.15622, LR test value: 0.53138, p-value: 0.46603
## Asymptotic standard error: 0.24853
##     z-value: -0.62858, p-value: 0.52962
## Wald statistic: 0.39511, p-value: 0.52962
## 
## Log likelihood: -154.7472 for error model
## ML residual variance (sigma squared): 5535.9, (sigma: 74.404)
## Number of observations: 27 
## Number of parameters estimated: 6 
## AIC: 321.49, (AIC for lm: 320.03)
# 3. Perbandingan AIC
aic_comparison <- data.frame(
  Model = c("OLS", "SAR", "SEM"),
  AIC = c(AIC(model_ols), AIC(model_sar), AIC(model_sem))
)

print("Perbandingan Nilai AIC:")
## [1] "Perbandingan Nilai AIC:"
print(aic_comparison)
##   Model      AIC
## 1   OLS 320.0258
## 2   SAR 322.0256
## 3   SEM 321.4944
# Menentukan Model Terbaik
best_model <- aic_comparison[which.min(aic_comparison$AIC), ]
print(paste("Model Terbaik berdasarkan AIC terendah adalah:", best_model$Model))
## [1] "Model Terbaik berdasarkan AIC terendah adalah: OLS"

Berdasarkan tabel perbandingan model, nilai efisiensi terbaik (AIC terendah) justru ditemukan pada model OLS standar (320,02), yang angkanya lebih kecil dibandingkan kedua model spasial SAR maupun SEM. Secara statistik, hasil ini menegaskan bahwa penambahan unsur spasial yang rumit tidak memberikan perbaikan apa pun pada kualitas prediksi, sehingga model regresi biasa (OLS) adalah pilihan terbaik dan paling efisien untuk data ini karena mampu bekerja maksimal tanpa perlu modifikasi tambahan.

4.4 Contoh Penerapan Desain Studi: Cross-Sectional

Penelitian ini menerapkan desain studi Cross-Sectional dengan pendekatan ekologi (Ecological Study), di mana pengukuran variabel paparan dan outcome dilakukan secara serentak pada unit analisis wilayah (Kabupaten/Kota) di Provinsi Jawa Barat pada tahun 2024.

4.4.1 Variabel Penelitian

Berdasarkan data yang digunakan, variabel diklasifikasikan sebagai berikut:

  • Variabel Dependen (Outcome): Kejadian penyakit Campak, yang diukur menggunakan indikator Prevalensi Campak (jumlah kasus per total populasi berisiko di setiap wilayah).
  • Variabel Independen (Exposure): Faktor-faktor risiko lingkungan dan pelayanan kesehatan yang diukur pada waktu yang sama, meliputi:
    1. Cakupan Imunisasi: Persentase bayi yang mendapatkan imunisasi dasar/campak.
    2. Kepadatan Penduduk: Jumlah penduduk per kilometer persegi sebagai indikator kerentanan transmisi.
    3. Sanitasi Layak: Persentase akses sanitasi yang memenuhi syarat kesehatan.

4.4.2 Populasi dan Sampel

  • Populasi Target: Seluruh penduduk di Provinsi Jawa Barat.
  • Unit Analisis: Karena data berbentuk agregat, unit analisis adalah 27 Kabupaten dan Kota di Jawa Barat.
  • Teknik Sampling: Menggunakan metode Total Sampling (Sensus Wilayah), di mana seluruh Kabupaten/Kota yang terdaftar di database Open Data Jabar dilibatkan dalam analisis tanpa terkecuali.

4.4.3 Potensi Bias

Mengingat desain studi ini bersifat observasional potong lintang pada level agregat, terdapat beberapa potensi bias yang perlu diperhatikan dalam interpretasi hasil:

  1. Ecological Fallacy (Kekeliruan Ekologis): Bias ini terjadi ketika kesimpulan mengenai individu ditarik langsung dari data agregat wilayah. Contohnya, wilayah dengan cakupan imunisasi tinggi mungkin memiliki kasus campak rendah, namun tidak menjamin bahwa individu yang sakit campak adalah individu yang tidak diimunisasi.
  2. Reporting Bias (Bias Pelaporan): Mengingat data bersumber dari laporan surveilans fasilitas kesehatan, terdapat risiko underreporting (kasus tidak terlaporkan) terutama di wilayah dengan akses kesehatan sulit atau kesadaran masyarakat yang rendah.
  3. Temporal Ambiguity: Pada desain cross-sectional, sulit untuk memastikan apakah faktor risiko (misalnya kepadatan penduduk atau kondisi sanitasi buruk) mendahului terjadinya lonjakan kasus campak, atau apakah kondisi tersebut kebetulan terjadi bersamaan (co-occurrence).

Bab 5. Kesimpulan dan Saran

5.1 Kesimpulan

Berdasarkan hasil analisis secara keseluruhan, penelitian ini menyimpulkan bahwa sebaran penyakit campak di Jawa Barat tahun 2024 menunjukkan pola ketimpangan data yang unik namun tidak memiliki ketergantungan antar-wilayah yang kuat. Meskipun secara visual terlihat adanya penumpukan kasus di kawasan urban utara, pengujian statistik membuktikan bahwa pola tersebut bersifat acak (random), sehingga model regresi standar (OLS) justru terbukti lebih valid dan efisien untuk menjelaskan fenomena ini dibandingkan model spasial yang rumit. Variabel kepadatan penduduk dan sanitasi terkonfirmasi sebagai faktor penentu utama yang sejalan dengan data kasus, sedangkan anomali pada data imunisasi mengisyaratkan bahwa risiko penularan di kota padat sangat tinggi terlepas dari status cakupan vaksinnya.

5.2 Saran

Berdasarkan temuan bahwa model spasial ternyata tidak terbukti signifikan dan model regresi standar hanya mampu menjelaskan sekitar 41% pola data, peneliti selanjutnya disarankan untuk tidak lagi memaksakan analisis berbasis wilayah makro, melainkan beralih mengeksplorasi data yang lebih detail di tingkat individu atau perilaku masyarakat. Rendahnya akurasi model saat ini mengindikasikan adanya variabel tersembunyi (lurking variable) di luar faktor lingkungan fisik seperti pola interaksi sosial atau mobilitas harian yang belum tertangkap, sehingga penambahan variabel baru yang lebih spesifik sangat diperlukan untuk mendapatkan gambaran statistik yang lebih utuh.

Berdasarkan data yang mengonfirmasi bahwa buruknya sanitasi dan tingginya kepadatan penduduk berbanding lurus dengan lonjakan kasus, pemerintah disarankan untuk tidak hanya terpaku pada target angka vaksinasi semata, tetapi mulai memprioritaskan intervensi perbaikan kualitas lingkungan di wilayah padat penduduk (urban) seperti Cirebon dan Bodebek. Data membuktikan bahwa wilayah dengan cakupan imunisasi tinggi pun bisa tetap mengalami ledakan kasus (“kebobolan”) jika kondisi lingkungannya mendukung penularan, sehingga masyarakat juga perlu sadar bahwa menjaga kebersihan lingkungan hunian adalah benteng pertahanan yang sama pentingnya dengan imunisasi untuk memutus rantai penularan yang sifatnya acak ini.

Bab 6. Daftar Pustaka

[1] World Health Organization, “Measles cases surge worldwide, infecting 10.3 million people in 2023,” WHO News Release, Nov. 14, 2024. [Online]. Available: https://www.who.int/news/item/14-11-2024-measles-cases-surge-worldwide--infecting-10.3-million-people-in-2023.

[2] United Nations, “Goal 3: Ensure healthy lives and promote well-being for all at all ages,” Sustainable Development Goals, 2023. [Online]. Available: https://sdgs.un.org/goals/goal3.

[3] Centers for Disease Control and Prevention (CDC), “Global Measles Outbreaks,” CDC Global Health, Dec. 10, 2025. [Online]. Available: https://www.cdc.gov/global-measles-vaccination/data-research/global-measles-outbreaks/index.html.

[4] N. M. R. Riastini dan I. M. Sutarga, “Hubungan Sanitasi Lingkungan dan Kejadian Penyakit Menular: Studi Literatur,” Jurnal Kesehatan Lingkungan, vol. 12, no. 1, pp. 20-25, 2020.

[5] L. Andriani, “Hubungan Karakteristik Balita dan Kepadatan Hunian dengan Kejadian Campak Klinis,” Jurnal Berkala Epidemiologi, vol. 5, no. 3, pp. 315-326, 2017.

[6] Badan Pusat Statistik, “Statistik Komuter Jabodetabek: Hasil Survei Komuter Jabodetabek 2024,” Jakarta: BPS RI, Nov. 2024.

[7] J. LeSage and R. K. Pace, Introduction to Spatial Econometrics. Boca Raton: CRC Press, 2009.

[8] Wassertheil-Smoller, S., & Smoller, J. (2015). Biostatistics and Epidemiology: A Primer for Health and Biomedical Professionals (4th ed.). New York, NY: Springer.

[9] Gordis, L. (2014). Epidemiology (5th ed.). Philadelphia, PA: Elsevier Saunders.

[10] Centers for Disease Control and Prevention. (2012). Principles of Epidemiology in Public Health Practice: An Introduction to Applied Epidemiology and Biostatistics (3rd ed.). Atlanta, GA: U.S. Department of Health and Human Services.

[11] Rothman, K. J., Greenland, S., & Lash, T. L. (2008). Modern Epidemiology (3rd ed.). Philadelphia, PA: Lippincott Williams & Wilkins.

[12] Pemerintah Provinsi Jawa Barat, “Jumlah Kasus Penyakit Campak Berdasarkan Kabupaten/Kota di Jawa Barat,” *Open Data Jabar* , 2024. [Online]. Available: https://opendata.jabarprov.go.id/id/dataset/jumlah-kasus-penyakit-campak-berdasarkan-kabupatenkota-di-jawa-barat.

[13] Pemerintah Provinsi Jawa Barat, “Jumlah Bayi yang Diimunisasi Campak Berdasarkan Kabupaten/Kota di Jawa Barat,” *Open Data Jabar*, 2024. [Online]. Available: https://opendata.jabarprov.go.id/id/dataset/jumlah-bayi-yang-diimunisasi-campak-berdasarkan-kabupatenkota-di-jawa-barat.

[14] Pemerintah Provinsi Jawa Barat, “Persentase Keluarga dengan Akses Sanitasi Layak (Jamban Sehat) Berdasarkan Kabupaten/Kota di Jawa Barat,” *Open Data Jabar*, 2024. [Online]. Available: https://opendata.jabarprov.go.id/id/dataset/persentase-keluarga-dengan-akses-sanitasi-layak-jamban-sehat-berdasarkan-kabupatenkota-di-jawa-barat.

[15] Pemerintah Provinsi Jawa Barat, “Kepadatan Penduduk Berdasarkan Kabupaten/Kota di Jawa Barat,” *Open Data Jabar*, 2024. [Online]. Available: https://opendata.jabarprov.go.id/id/dataset/kepadatan-penduduk-berdasarkan-kabupatenkota-di-jawa-barat.

[16] Pemerintah Provinsi Jawa Barat, “Proyeksi Penduduk Berdasarkan Kabupaten/Kota di Jawa Barat,” *Open Data Jabar*, 2024. [Online]. Available: https://opendata.jabarprov.go.id/id/dataset/proyeksi-penduduk-berdasarkan-kabupatenkota-di-jawa-barat.