Analisis Epidemiologi pada Penyakit DBD

LAPORAN UJIAN AKHIR SEMESTER

Epidemiologi

Disusun oleh:

Andini Dwi Kurnia Putri (140610230082)

Dosen Pengampu :

Dr. I Gede Nyoman Mindra Jaya, M.Si

PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

1 Pendahuluan

1.1 Latar Belakang

Di Indonesia, Demam Berdarah Dengue (DBD) masih menjadi masalah besar bagi kesehatan masyarakat, dengan insidensi yang meningkat dan menyebar di seluruh provinsi. Siklus hidup nyamuk Aedes aegypti sangat didukung oleh iklim tropis Indonesia, tetapi ada banyak faktor lain yang meningkatkan risiko penularan penyakit ini. Jumlah kasus yang berbeda di berbagai wilayah menunjukkan bahwa DBD bukan hanya masalah biologis; itu adalah hasil dari interaksi kompleks antara agen penyakit, inang, dan lingkungan. Untuk menentukan apakah suatu wilayah merupakan pusat penyebaran penyakit yang memerlukan intervensi prioritas, sangat penting untuk memahami pola distribusi penyakit secara geografis.

Menurut epidemiologi lingkungan, faktor sosial dan infrastruktur dasar sangat memengaruhi penyebaran penyakit menular. Faktor demografi seperti kepadatan penduduk mempercepat penyebaran virus antarmanusia, dan faktor lingkungan seperti akses sanitasi yang buruk membantu menyediakan tempat perindukan vektor. Sebaliknya, kapasitas suatu wilayah untuk mencegah dan menanggulangi wabah dipengaruhi oleh kerentanan sosial, yang ditunjukkan oleh persentase penduduk miskin dan ketersediaan layanan kesehatan, seperti jumlah rumah sakit. Faktor risiko lingkungan dan sosial ini seringkali bersifat lokal dan spesifik, sehingga kesimpulan yang mengabaikan aspek kewilayahan dapat menjadi kurang akurat.

Karena efek autokorelasi spasial, di mana kejadian penyakit di suatu lokasi cenderung dipengaruhi oleh kondisi di sekitarnya (hukum Tobler), pendekatan statistik konvensional sering kali tidak cukup untuk memodelkan penyebaran penyakit menular. Tujuan dari penelitian ini adalah untuk mengeksplorasi faktor epidemiologi DBD di Indonesia dengan menggunakan metode spasial. Diharapkan dapat dibuat estimasi yang lebih akurat tentang hubungan antara sanitasi, kepadatan penduduk, kemiskinan, dan fasilitas kesehatan terhadap insidensi DBD melalui pemetaan distribusi penyakit, identifikasi klaster wilayah (LISA), dan pemodelan regresi spasial yang diharapkan akan membantu menciptakan kebijakan kesehatan berbasis wilayah yang lebih efektif.

1.2 Identifikasi Masalah

1. Bagaimana distribusi geografis dan insidensi Demam Berdarah Dengue (DBD) di berbagai provinsi di Indonesia, serta daerah mana saja yang diidentifikasi sebagai hotspot atau coldspot?
2. Apakah ada autokorelasi spasial yang signifikan yang menunjukkan ketergantungan antarwilayah dalam penyebaran kasus DBD?
3. Bagaimana besar risiko (asosiasi) antara kondisi sanitasi lingkungan yang buruk terhadap peningkatan kejadian DBD berdasarkan perhitungan ukuran epidemiologi?
4. Model manakah yang paling cocok untuk menjelaskan bagaimana faktor-faktor seperti sanitasi, kepadatan penduduk, kemiskinan, dan fasilitas kesehatan memengaruhi insidensi DBD?

1.3 Tujuan Penelitian

1. Mendeskripsikan pola epidemiologi dan distribusi spasial penyakit DBD di tingkat provinsi pada tahun 2024
2. Menghitung ukuran-ukuran epidemiologi
3. Mengidentifikasi autokorelasi spasial global dan lokal di seluruh provinsi Indonesia
4. Menguji signifikansi autokorelasi spasial global untuk membuktikan apakah kejadian DBD di suatu provinsi dipengaruhi oleh kondisi di provinsi tetangganya
5. Membangun dan membandingkan kinerja model guna menentukan model terbaik yang mampu menjelaskan pengaruh faktor lingkungan dan sosial terhadap DBD.

1.4 Batasan Penelitian

1. Penelitian ini dikhususkan pada analisis spasial penyakit Demam Berdarah Dengue (DBD).
2. Unit analisis adalah provinsi di Indonesia tahun 2024
3. Variabel dependen yang dimodelkan adalah Angka Insidensi (Incidence Rate), yaitu jumlah kasus DBD per 100.000 penduduk, bukan sekadar jumlah kasus absolut, untuk menghindari bias akibat perbedaan jumlah populasi antarprovinsi.
4. Variabel independen adalah Persentase akses sanitasi layak, kepadatan penduduk, ketersediaan fasilitas kesehatan, dan persentase penduduk miskin. 5.Matriks pembobotan spasial menggunakan metode k-Nearest Neighbors (KNN) dengan penentuan nilai k optimal berdasarkan nilai Moran’s I tertinggi. Analisis asosiasi risiko (Odds Ratio) menggunakan kategorisasi data berdasarkan nilai median (cut-off median).

2 Tinjauan Pustaka

2.1 Demam Berdarah Dengue (DBD)

Demam Berdarah Dengue (DBD) adalah penyakit infeksi akut yang disebabkan oleh virus Dengue dan ditularkan terutama melalui gigitan nyamuk Aedes aegypti dan Aedes albopictus. DBD sangat umum di Indonesia karena iklimnya yang tropis. Indikator utama yang digunakan untuk mengukur beban penyakit ini di suatu wilayah adalah Incidence Rate (Angka Insidensi), yang menggambarkan jumlah kasus baru per populasi berisiko (biasanya per 100.000 penduduk) dalam periode waktu tertentu. Karena penyebaran DBD seringkali tidak acak, melainkan membentuk pola spasial yang dipengaruhi oleh kondisi lokal, pemahaman tentang distribusi insidensi sangat penting.

2.2 Determinan Agent-Host-Environment pada DBD

Agent adalah mikroorganisme atau pathogen seperti virus, bakteri, parasit, dan sebagainya yang dapat menginfeksi manusia atau makhluk hidup. Host merupakan manusia atau makhluk hidup yang dapat diinfeksi oleh agent. Semua faktor eksternal yang mempengaruhi kemungkinan terjadinya kontak antara host dan agen disebut environment.

Pada penelitian ini, agennya adalah virus Flavivirus yang dibawa oleh vektor Aedes Aegypti dan beberapa spesies nyamuk lainnya dalam genus Aedes. Environment pada penelitian ini adalah persentase rumah tangga dengan akses sanitasi layak, kepadatan penduduk per km2, jumlah desa dengan fasilitas rumah sakit, dan persentase penduduk miskin. Risiko paparan terhadap agen penyakit berbasis lingkungan meningkat karena sanitasi yang tidak memadai; kepadatan penduduk yang tinggi mempercepat penularan antar orang; kekurangan fasilitas kesehatan memperlambat deteksi dan penanganan kasus; dan kemiskinan meningkatkan kerentanan masyarakat terhadap infeksi.

2.3 Desain Studi

Metode utama untuk menguji hipotesis tentang hubungan antara paparan dan penyakit adalah pengertian dari desain studi dalam konteks epidemiologi. Studi biasanya dikategorikan menjadi dua kategori: observasional dan eksperimental. Studi observasional analitik termasuk Studi Cross-Sectional (mengukur paparan dan penyakit secara bersamaan), Studi Kasus-Kontrol (membandingkan riwayat paparan pada kelompok sakit dan tidak sakit), dan Studi Kohort (melacak perkembangan penyakit di kelompok terpapar dan tidak terpapar). Setiap desain memiliki kelebihan dan kelemahan, serta bias, seperti bias mengingat atau kehilangan mengingat, yang dapat memengaruhi validitas kesimpulan kausalitas. Akibatnya, memilih desain yang tepat sangat penting.

2.4 Ukuran Frekuensi Penyakit

Dalam epidemiologi deskriptif, pengukuran frekuensi penyakit merupakan langkah fundamental untuk mengkarakterisasi pola kejadian penyakit di suatu populasi.

2.4.1 Incidence Rate (Angka Insidensi)

Laju kejadian kasus baru penyakit dalam suatu populasi berisiko disebut insidensi. Penggunaan angka mutlak, atau jumlah kasus saja, akan bias karena variasi populasi di seluruh Indonesia. Oleh karena itu, untuk memungkinkan perbandingan yang adil antarwilayah, digunakan standarisasi per konstanta penduduk (k), biasanya per 100.000 orang. \[ \text{Insidensi} = \frac{\text{Jumlah Kasus Baru}}{\text{Populasi}} \times k \]

2.4.2 Prevalensi

Prosentase individu dalam populasi yang terkena penyakit pada rentang waktu tertentu disebut prevalensi. Diperlukan jumlah kasus yang dilakukan selama periode penelitian untuk nilai prevalensi. \[ \text{Prevalensi} = \frac{\text{Jumlah Total Kasus}}{\text{Jumlah Populasi}} \]

2.4.3 Case Fatality Rate (CFR)

CFR adalah persentase kasus yang meninggal dari total kasus. \[ CFR = \frac{\text{Jumlah kematian akibat penyakit}}{\text{Jumlah kasus penyakit}} \times 100\% \]

2.4.4 Attack Rate

Jumlah orang yang terkena penyakit selama wabah atau epidemi dikenal sebagai tingkat infeksi. \[ \text{Attack Rate} = \frac{\text{Jumlah kasus selama wabah}}{\text{Populasi}} \times 100\% \]

2.4.5 Mortality Rate

Mortality rate adalah jumlah kematian per populasi.

2.5 Ukuran Asosiasi

Tabel kontingensi 2x2 digunakan untuk mengukur hubungan deterministik antara faktor paparan (eksposur) dan kasus penyakit.

2.5.1 Odds Ratio (OR)

Kemungkinan terjadinya penyakit pada kelompok terpapar dibandingkan dengan kemungkinan terjadinya penyakit pada kelompok tidak terpapar dikenal sebagai odds ratio. OR digunakan untuk menentukan seberapa besar kemungkinan suatu wilayah dengan sanitasi buruk (terpapar) memiliki insidensi DBD yang tinggi dibandingkan dengan wilayah dengan sanitasi baik. \[ OR = \frac{a \times d}{b \times c} \] Dimana jika nilai \(OR > 1\), maka faktor paparan dianggap sebagai faktor risiko yang meningkatkan kejadian penyakit.

2.5.2 Risk Ratio (RR)

Risiko rasio, juga dikenal sebagai risiko relatif, menilai kemungkinan (risiko) munculnya penyakit pada kelompok terpapar dibandingkan dengan kelompok tidak terpapar. \[ RR = \frac{Insidensi_{terpapar}}{Insidensi_{tidakterpapar}} \]

2.5.3 Attributable Risk

AR menghitung proporsi kasus yang dapat dihindari dalam kelompok terpapar jika paparan dicegah.

2.6 Eksplorasi Data Spasial (ESDA)

2.6.1 Indeks Moran’s I (Global)

Global Moran’s I adalah statistik uji korelasi yang mengukur autokorelasi spasial umum di wilayah studi. Nilai indeks antara -1 dan +1. Nilai yang positif menunjukkan pengelompokan (clustering) area dengan nilai sejenis yang tinggi atau rendah, sedangkan nilai yang negatif menunjukkan pola penyebaran.

2.6.2 Local Indicators of Spatial Association (LISA)

Karena Global Moran’s I hanya memberikan satu nilai rangkuman untuk seluruh wilayah, digunakan statistik lokal (Local Moran’s I) untuk mengidentifikasi lokasi spesifik dari klaster spasial. Peta LISA membagi wilayah menjadi empat kuadran: 1. High-High (Hotspot): Wilayah insidensi tinggi dikelilingi wilayah insidensi tinggi. 2. Low-Low (Coldspot): Wilayah insidensi rendah dikelilingi wilayah insidensi rendah. 3. High-Low & Low-High (Spatial Outliers): Wilayah yang nilainya berbeda signifikan dengan tetangganya.

2.7 Matriks Pembobotan Spasial (Spatial Weights Matrix)

2.7.1 k-Nearest Neighbors (KNN)

Studi ini menggunakan metode k-Nearest Neighbors (KNN). KNN mendefinisikan tetangga berdasarkan jarak titik pusat (centroid) terdekat sejumlah k. Metode ini menjamin bahwa setiap wilayah memiliki jumlah tetangga yang sama, sehingga struktur konektivitas spasial lebih seimbang. Ini berbeda dengan pembobotan berbasis persinggungan batas (contiguity) yang dapat bias pada wilayah kepulauan yang terpisah laut (seperti provinsi di Indonesia).

2.8 Pemodelan Ekonometrika Spasial

2.8.1 Spatial Autoregressive Model (SAR)

Pada Model SAR, yang juga dikenal sebagai Model Lag Spatial, variabel dependen di suatu wilayah dipengaruhi secara langsung oleh variabel dependen di wilayah sekitarnya. Variabel dependen ini dikenal sebagai efek difusi atau penularan. \[ y = \rho W y + X \beta + \epsilon \] Dimana \(\rho\) adalah koefisien lag spasial.

2.8.2 Spatial Error Model (SEM)

Model SEM mengasumsikan bahwa ketergantungan spasial tidak terjadi pada variabel dependen, melainkan pada komponen error. Hal ini biasanya disebabkan oleh adanya variabel pengganggu yang memiliki pola spasial namun tidak dimasukkan ke dalam model (omitted spatial variables). \[ y = X \beta + u, \quad u = \lambda W u + \epsilon \] Dimana \(\lambda\) adalah koefisien error spasial.

3 Metodologi

3.1 Sumber Data

Penelitian ini menggunakan data dari Profil Kesehatan Indonesia tahun 2024 dari Kementerian Kesehatan Republik Indonesia. Semua 38 provinsi Indonesia adalah unit observasi dalam penelitian ini. Berikut ini adalah variabel yang digunakan:

data <-read.csv("D:/SEMESTER 5/Spasial/Data Spasial Epidem Andin.csv")
library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data <- data %>%
  rename(
    Akses_Sanitasi = Persentase.Rumah.Tangga.Memiliki.Akses.terhadap.Sanitasi.Layak..Persen.,
    Kepadatan_Penduduk = Kepadatan.Penduduk.per.km.persegi..Km2.,
    Jumlah_Rumah_Sakit = Jumlah.Desa.yang.Memiliki.Rumah.Sakit,
    Penduduk_Miskin = Presentase.Penduduk.Miskin,
    Jumlah_Penduduk = Jumlah.Penduduk..Ribu.
  )
data

Model dan mapping akan dibuat untuk variabel penyakit DBD ini. Selanjutnya, Presentase Rumah Tangga yang Memiliki Akses Sanitasi Layak, Kepadatan Penduduk, Jumlah Desa yang Memiliki Rumah Sakit, dan Presentase Penduduk Miskin akan menjadi variabel independen yang sekiranya akan dianalisis apakah mempengaruhi dan ada aspek autokorelasi spasial.

4 Hasil dan Pembahasan

4.1 Dashboard dan Video Paparan Project

4.1.1 Dashboard

Berikut adalah link untuk dashboard yang telah dibuat : https://dashboard-dbd-indonesia-hvdrv2f5jonwr2wayd4jmy.streamlit.app/ Terdapat beberapa keterbatasan, salah satunya peta yang dapat tersajikan masih yang 34 provinsi Indonesia.

4.1.2 Video Paparan Project

Berikut adalah link untuk video paparan project https://youtu.be/5y5KSttquGI?feature=shared

4.2 Ukuran Epidemiologi

Karena keterbatasan akses data, ukuran epidemiologi hanyalah insidensi. Data ini hanya mencakup total jumlah kasus pada tahun 2024 dan jumlah penduduk per provinsi di Indonesia.

4.2.1 Insidensi

options(scipen=10)
data <- data %>%
  mutate(
    Jumlah_Penduduk = Jumlah_Penduduk * 1000,
    DBD_insid        = (DBD / Jumlah_Penduduk) * 100000
    )
insid <- data %>%
  select(Provinsi, DBD_insid)
insid

Nilai Insidensi ini merepresentasikan jumlah kasus baru DBD yang ditemukan per 100.000 penduduk di setiap provinsi. Semakin tinggi angka insidensi, semakin tinggi pula laju penularan dan risiko kesehatan di wilayah tersebut. Penggunaan insidensi (bukan jumlah kasus mutlak) bertujuan untuk menstandarisasi data, sehingga perbandingan antarprovinsi menjadi tidak bias(menghindari asumsi bahwa provinsi padat pasti memiliki risiko penyakit yang lebih tinggi)

library(sf)

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(dplyr)
library(spdep)

## Warning: package 'spdep' was built under R version 4.3.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.3.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(ggplot2)
library(viridis)

## Warning: package 'viridis' was built under R version 4.3.3

## Loading required package: viridisLite

library(readr)

## Warning: package 'readr' was built under R version 4.3.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.3

## Warning: package 'tidyr' was built under R version 4.3.3

## Warning: package 'forcats' was built under R version 4.3.3

## Warning: package 'lubridate' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

batas_prov <- readRDS("D:/SEMESTER 5/Spasial/Batas_Provinsi_Indonesia.rds")
batas_prov <- batas_prov %>%
  mutate(Provinsi = as.character(WADMPR)) %>%
  st_make_valid() %>%
  st_zm(drop = TRUE, what = "ZM") %>%
  st_buffer(0)
batas_prov <- batas_prov %>%
  mutate(Provinsi = as.character(WADMPR)) %>%
  st_make_valid() %>%
  st_zm(drop = TRUE, what = "ZM") %>%
  st_buffer(0)
# Gabungkan shapefile dan data
data_map <- batas_prov %>%
  left_join(data, by = "Provinsi")
# Menentukan K untuk Bobot Spasial KNN
sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

4.2.2 Ukuran Asosiasi

library(epitools)
# 1. Kategorisasi Data (Cut-off menggunakan Median)
median_sanitasi <- median(data_map$Akses_Sanitasi, na.rm = TRUE)
median_dbd      <- median(data_map$DBD_insid, na.rm = TRUE) 
data_epi <- data_map %>%
  st_drop_geometry() %>%
  mutate(
    Exposure_Sanitasi = ifelse(Akses_Sanitasi < median_sanitasi, "Buruk (Exposed)", "Baik (Unexposed)"),
    Outcome_DBD       = ifelse(DBD_insid > median_dbd, "Kasus Tinggi", "Kasus Rendah")
  )
# 2. Membuat Tabel 2x2
tabel_2x2 <- table(data_epi$Exposure_Sanitasi, data_epi$Outcome_DBD)
tabel_2x2 <- tabel_2x2[c("Buruk (Exposed)", "Baik (Unexposed)"), c("Kasus Tinggi", "Kasus Rendah")]
print(tabel_2x2)

##                   
##                    Kasus Tinggi Kasus Rendah
##   Buruk (Exposed)             9           10
##   Baik (Unexposed)           10            9

# 3. Odds Ratio (OR)
hasil_or <- oddsratio(tabel_2x2, method = "wald")
cat("\n--- Odds Ratio (OR) ---\n")

## 
## --- Odds Ratio (OR) ---

print(hasil_or$measure)

##                   odds ratio with 95% C.I.
##                    estimate     lower   upper
##   Buruk (Exposed)      1.00        NA      NA
##   Baik (Unexposed)     0.81 0.2266658 2.89457

# 4. Risk Ratio (RR)
hasil_rr <- riskratio(tabel_2x2, method = "wald")
cat("\n--- Risk Ratio (RR) ---\n")

## 
## --- Risk Ratio (RR) ---

print(hasil_rr$measure)

##                   risk ratio with 95% C.I.
##                    estimate     lower    upper
##   Buruk (Exposed)       1.0        NA       NA
##   Baik (Unexposed)      0.9 0.4756749 1.702844

or_val <- hasil_or$measure[2, 1]
cat("\nInterpretasi: Wilayah dengan sanitasi buruk memiliki peluang kejadian DBD tinggi sebesar", 
    round(or_val, 2), "kali lipat dibandingkan wilayah dengan sanitasi baik.\n")

## 
## Interpretasi: Wilayah dengan sanitasi buruk memiliki peluang kejadian DBD tinggi sebesar 0.81 kali lipat dibandingkan wilayah dengan sanitasi baik.

4.3 Analisis Deskriptif Spasial

centroids <- st_centroid(data_map)

## Warning: st_centroid assumes attributes are constant over geometries

## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data

coords <- st_coordinates(centroids)
# Fungsi untuk hitung Moran's I berdasarkan k dan variabel target
calculate_moran <- function(k, var) {
  nb <- knn2nb(knearneigh(coords, k = k))
  lw <- nb2listw(nb, style = "W")
  moran_res <- moran.test(var, lw, zero.policy = TRUE)
  return(moran_res$estimate["Moran I statistic"])
}
# Fungsi untuk menentukan k optimal dengan fallback ke 3 jika hasil = 1
find_optimal_k <- function(var, disease_name) {
  k_range <- 1:15
  moran_values <- sapply(k_range, calculate_moran, var = var)
  k_optimal <- k_range[which.max(moran_values)]
  if (k_optimal == 1) {
    message("⚠️  k optimal untuk ", disease_name, " = 1, diganti menjadi 3 (default).")
    k_optimal <- 3
  } else {
    message("✅ k optimal untuk ", disease_name, " =", k_optimal)
  }
  # Plot nilai Moran’s I terhadap k
  plot(k_range, moran_values, type = "b", pch = 19,
       xlab = "Nilai k (jumlah tetangga)",
       ylab = "Nilai Moran's I",
       main = paste("Pemilihan Nilai k Optimal -", disease_name))
  abline(v = k_optimal, col = "red", lty = 2)
  text(k_optimal, moran_values[which.max(moran_values)],
       labels = paste("k =", k_optimal), pos = 4)
  return(k_optimal)
}
k_dbd        <- find_optimal_k(data_map$DBD, "DBD")

## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 12
## sub-graphs

## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points
## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points
## Warning in knearneigh(coords, k = k): k greater than one-third of the number of
## data points

## ✅ k optimal untuk DBD =3

nb_knn_dbd        <- knn2nb(knearneigh(coords, k = k_dbd))
listw_knn_dbd     <- nb2listw(nb_knn_dbd, style = "W", zero.policy = TRUE)
listw <- nb2listw(nb_knn_dbd, style = "W", zero.policy = TRUE)
weights <- listw_knn_dbd
cat("   - DBD: k =", k_dbd, "\n")

##    - DBD: k = 3

diseases <- c("DBD")
for (d in diseases) {
  insid_col <- paste0(d, "_insid")
  p <- ggplot(data_map) +
    geom_sf(aes(fill = .data[[insid_col]]), color = "white", size = 0.3) +
    scale_fill_viridis_c(
      option = "inferno",
      na.value = "grey90",
      direction = -1,
      name = paste("Insidensi", d, "(per 100.000)")
    ) +
    labs(
      title = paste("Peta Insidensi", d, "di Indonesia"),
      subtitle = "Gradasi warna menunjukkan intensitas kasus per 100.000 penduduk",
      caption = "Sumber data: Laporan Kesehatan Tahun 2024 Kementerian Kesehatan Republik Indonesia "
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 11, hjust = 0.5),
      legend.position = "right"
    )
  print(p)
  ggsave(
    filename = paste0("Peta_Gradasi_", d, ".png"),
    plot = p, width = 8, height = 5, dpi = 300
  )
}

Secara visual, peta ini menunjukkan pola penyebaran DBD, di mana Kalimantan terlihat dengan insidensi paling ekstrem dibanding wilayah timur Indonesia yang cenderung lebih rendah.

diseases <- c("DBD")
for (d in diseases) {
  insid_col <- paste0(d, "_insid")
  # Hitung nilai kuantil (Q0, Q1, Q2, Q3, Q4)
  quantile_values <- quantile(
    data_map[[insid_col]],
    probs = seq(0, 1, 0.25),
    na.rm = TRUE
  )
  # Cetak ke konsol
  cat("\n==============================\n")
  cat("📊 Nilai Kuantil Insidensi -", d, "\n")
  print(quantile_values)
  cat("==============================\n\n")
  # Membuat kategori kuartil (Q1 - Q4)
  data_map <- data_map %>%
    mutate(
      quartile = cut(
        .data[[insid_col]],
        breaks = quantile_values,
        include.lowest = TRUE,
        labels = c("Q1 (terendah)", "Q2", "Q3", "Q4 (tertinggi)")
      )
    )
  # Plot dengan warna diskrit berdasarkan kuartil
  p <- ggplot(data_map) +
    geom_sf(aes(fill = quartile), color = "white", size = 0.3) +
    scale_fill_viridis_d(
      option = "inferno",
      direction = -1,
      name = paste("Kuartil Insidensi", d)
    ) +
    labs(
      title = paste("Peta Kuartil Insidensi", d, "di Indonesia"),
      subtitle = "Warna menunjukkan kelompok kuartil jumlah kasus per 100.000 penduduk",
      caption = "Sumber data: Laporan Kesehatan Tahun 2024 Kementerian Kesehatan Republik Indonesia"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 11, hjust = 0.5),
      legend.position = "right"
    )
  print(p)
  ggsave(
    filename = paste0("Peta_Kuartil_", d, ".png"),
    plot = p, width = 8, height = 5, dpi = 300
  )
}

## 
## ==============================
## 📊 Nilai Kuantil Insidensi - DBD 
##        0%       25%       50%       75%      100% 
##   0.00000  53.89009  76.95545 123.68850 269.45797 
## ==============================

Peta kuartil memperlihatkan wilayah Kalimantan teridentifikasi sebagai zona kerentanan tinggi (Q4) dengan insidensi melampaui 123 kasus per 100.000 penduduk, berbeda dengan wilayah timur yang mayoritas berada di tingkat terendah.

4.4 Uji Autokorelasi Spasial

diseases <- c("DBD")
# Data frame hasil
autocorr_results <- data.frame(
  Penyakit = character(),
  Moran_I = numeric(),
  Moran_p = numeric(),
  Geary_C = numeric(),
  Geary_p = numeric(),
  stringsAsFactors = FALSE
)
# Loop tiap penyakit
for (d in diseases) {
  insid_col <- paste0(d, "_insid")
  values <- data_map[[insid_col]]
  # Hilangkan NA
  valid_idx <- !is.na(values)
  values <- values[valid_idx]
  if (length(values) > 3) {
    # Moran's I
    moran_res <- moran.test(values, listw, zero.policy = TRUE)
    # Geary's C
    geary_res <- geary.test(values, listw, zero.policy = TRUE)
    autocorr_results <- rbind(
      autocorr_results,
      data.frame(
        Penyakit = d,
        Moran_I = round(moran_res$estimate[1], 4),
        Moran_p = round(moran_res$p.value, 4),
        Geary_C = round(geary_res$estimate[1], 4),
        Geary_p = round(geary_res$p.value, 4)
      )
    )
  }
}
print(autocorr_results)

##                   Penyakit Moran_I Moran_p Geary_C Geary_p
## Moran I statistic      DBD  0.1271  0.0918  0.8534  0.1255

Dari hasil pengujian autokorelasi spasial didapatkan bahwa DBD memiliki autokorelasi spasial yang lemah. Hal ini menunjukkan bahwa DBD dapat menyebar melalui wilayah.

calculate_morans_i <- function(data, disease_var, weights_list) {
  insid_var <- paste0(disease_var, "_insid")
  insid_values <- st_drop_geometry(data)[[insid_var]]
  valid_idx <- !is.na(insid_values)
  insid_values <- insid_values[valid_idx]
  # Pastikan selalu mengembalikan kolom yang sama
  out_template <- data.frame(
    Disease = disease_var,
    Statistic = NA_real_,
    P_value = NA_real_,
    Expected = NA_real_,
    Variance = NA_real_,
    stringsAsFactors = FALSE
  )
  if (length(insid_values) < 3) {
    return(out_template)
  }
  # Ambil nb dari objek weights_list (diasumsikan struktur: list(neighbors = nb_obj, weights = listw_obj))
  nb_full <- weights_list$neighbors
  # subset nb sesuai indeks valid
  nb_subset <- subset(nb_full, valid_idx)
  # buat listw dari nb_subset
  weights_subset <- nb2listw(nb_subset, style = "W", zero.policy = TRUE)
  morans <- tryCatch(
    moran.test(insid_values, weights_subset, zero.policy = TRUE),
    error = function(e) return(NULL)
  )
  if (is.null(morans)) {
    return(out_template)
  }
  # Ekstraksi hasil secara aman
  stat <- as.numeric(morans$statistic["Moran I statistic"])
  pval <- as.numeric(morans$p.value)
  # estimate bisa berupa vektor bernama; ambil elemen sesuai nama jika ada
  expected <- if (!is.null(morans$estimate["Expectation"])) {
    as.numeric(morans$estimate["Expectation"])
  } else if (length(morans$estimate) >= 2) {
    as.numeric(morans$estimate[2])
  } else NA_real_
  variance <- if (!is.null(morans$estimate["Variance"])) {
    as.numeric(morans$estimate["Variance"])
  } else if (length(morans$estimate) >= 3) {
    as.numeric(morans$estimate[3])
  } else NA_real_
  out_template$Statistic <- stat
  out_template$P_value <- pval
  out_template$Expected <- expected
  out_template$Variance <- variance
  return(out_template)
}
calculate_lisa <- function(data, disease_var, weights_list) {
  data_df <- st_drop_geometry(data)
  insid_var <- paste0(disease_var, "_insid")
  insid <- data_df[[insid_var]]
  valid_idx <- !is.na(insid)
  insid_valid <- insid[valid_idx]
  if (length(insid_valid) < 3) return(NULL)
  nb_full <- weights_list$neighbors
  nb_subset <- subset(nb_full, valid_idx)
  weights_subset <- nb2listw(nb_subset, style = "W", zero.policy = TRUE)
  # Jalankan Local Moran’s I dengan proteksi error
  lisa <- tryCatch(
    localmoran(insid_valid, weights_subset, zero.policy = TRUE),
    error = function(e) return(NULL)
  )
  if (is.null(lisa)) return(NULL)
  # Konversi ke data.frame agar mudah diakses
  lisa <- as.data.frame(lisa)
  # Coba deteksi kolom p-value secara otomatis
  p_col <- grep("Pr", names(lisa), value = TRUE)
  if (length(p_col) == 0) p_col <- names(lisa)[5]  # fallback ke kolom ke-5
  p_values <- lisa[[p_col[1]]]
  # Buat output sf
  data_subset <- data[valid_idx, ]
  data_subset$insid <- insid_valid
  data_subset$Ii <- lisa[, "Ii"]
  data_subset$p_value <- p_values
  data_subset$insid_std <- as.numeric(scale(insid_valid)[, 1])
  lag_vals <- lag.listw(weights_subset, insid_valid, zero.policy = TRUE)
  data_subset$Lag_std <- as.numeric(scale(lag_vals)[, 1])
  # Klasifikasi cluster
  data_subset$Cluster <- "Not Significant"
  sig_mask <- data_subset$p_value < 0.05
  data_subset$Cluster[sig_mask & data_subset$insid_std > 0 & data_subset$Lag_std > 0] <- "HH (Hotspot)"
  data_subset$Cluster[sig_mask & data_subset$insid_std < 0 & data_subset$Lag_std < 0] <- "LL (Coldspot)"
  data_subset$Cluster[sig_mask & data_subset$insid_std > 0 & data_subset$Lag_std < 0] <- "HL (Outlier)"
  data_subset$Cluster[sig_mask & data_subset$insid_std < 0 & data_subset$Lag_std > 0] <- "LH (Outlier)"
  data_subset$Disease <- disease_var
  return(data_subset)
}
# Contoh: jika Anda punya hanya tiga disease KNN
weights_per_disease <- list(
  DBD = list(neighbors = nb_knn_dbd, weights = listw_knn_dbd))
diseases <- names(weights_per_disease)  # gunakan nama di weights_per_disease
# --- Global Moran’s I (gabungkan hasil menjadi dataframe)
morans_results <- lapply(diseases, function(d) {
  calculate_morans_i(data_map, d, weights_per_disease[[d]])
})
morans_df <- do.call(rbind, morans_results)
print(morans_df)

##   Disease Statistic    P_value    Expected   Variance
## 1     DBD        NA 0.09181663 -0.02702703 0.01343638

lisa_results <- lapply(diseases, function(d) {
  calculate_lisa(data_map, d, weights_per_disease[[d]])
})
names(lisa_results) <- diseases
# === Visualisasi Peta LISA ===
for (d in diseases) {
  lisa_data <- lisa_results[[d]]
  if (is.null(lisa_data)) {
    message("Tidak cukup data LISA untuk ", d)
    next
  }
  p <- ggplot(lisa_data) +
    geom_sf(aes(fill = Cluster), color = "white", size = 0.3) +
    scale_fill_manual(
      values = c(
        "HH (Hotspot)" = "#E31A1C",
        "LL (Coldspot)" = "#1F78B4",
        "HL (Outlier)" = "#FEB24C",
        "LH (Outlier)" = "#A1D99B",
        "Not Significant" = "grey90"
      )
    ) +
    labs(
      title = paste("Peta LISA (Local Moran’s I):", d),
      subtitle = "Hotspot dan Coldspot berdasarkan insidensi penyakit",
      fill = "Cluster"
    ) +
    theme_minimal()
  print(p)
}

Hasil analisis LISA mendeteksi spatial outlier (Low-High) di sebagian Kalimantan dan Sulawesi, yang artinya wilayah-wilayah tersebut memiliki insidensi rendah padahal dikelilingi oleh tetangga dengan insidensi tinggi.

4.5 Modeling

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.3.3

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)

## Warning: package 'car' was built under R version 4.3.3

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

equation <- DBD_insid ~ Akses_Sanitasi + Kepadatan_Penduduk + Jumlah_Rumah_Sakit + Penduduk_Miskin
# 2. Model OLS 
model_ols <- lm(equation, data = data_map)
summary(model_ols)

## 
## Call:
## lm(formula = equation, data = data_map)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -85.08 -32.59 -15.92  11.80 161.66 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)
## (Intercept)         59.189998 100.186456   0.591    0.559
## Akses_Sanitasi       0.668889   0.969524   0.690    0.495
## Kepadatan_Penduduk   0.001280   0.003735   0.343    0.734
## Jumlah_Rumah_Sakit  -0.103004   0.114367  -0.901    0.374
## Penduduk_Miskin     -1.451954   2.343634  -0.620    0.540
## 
## Residual standard error: 57.18 on 33 degrees of freedom
## Multiple R-squared:  0.1101, Adjusted R-squared:  0.002215 
## F-statistic: 1.021 on 4 and 33 DF,  p-value: 0.4111

# Cek Multikolinearitas
cat("\n--- Cek VIF ---\n")

## 
## --- Cek VIF ---

print(vif(model_ols))

##     Akses_Sanitasi Kepadatan_Penduduk Jumlah_Rumah_Sakit    Penduduk_Miskin 
##           2.416490           1.075299           1.098714           2.520005

# 3. Uji Lagrange Multiplier 
lm_test <- lm.LMtests(model_ols, listw_knn_dbd, test = c("LMerr","LMlag"))

## Please update scripts to use lm.RStests in place of lm.LMtests

print(lm_test)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = equation, data = data_map)
## test weights: listw
## 
## RSerr = 0.0012473, df = 1, p-value = 0.9718
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = equation, data = data_map)
## test weights: listw
## 
## RSlag = 0.19895, df = 1, p-value = 0.6556

# 4. Model SAR (Spatial Lag)
model_sar <- lagsarlm(equation, data = data_map, listw = listw_knn_dbd)
summary(model_sar)

## 
## Call:lagsarlm(formula = equation, data = data_map, listw = listw_knn_dbd)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.517 -32.126 -18.783  12.878 162.596 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)        42.8886456 97.3272562  0.4407   0.6595
## Akses_Sanitasi      0.7057031  0.9036962  0.7809   0.4349
## Kepadatan_Penduduk  0.0012745  0.0034699  0.3673   0.7134
## Jumlah_Rumah_Sakit -0.0941956  0.1061519 -0.8874   0.3749
## Penduduk_Miskin    -1.1451480  2.2065940 -0.5190   0.6038
## 
## Rho: 0.097165, LR test value: 0.21216, p-value: 0.64508
## Asymptotic standard error: 0.1978
##     z-value: 0.49124, p-value: 0.62326
## Wald statistic: 0.24131, p-value: 0.62326
## 
## Log likelihood: -204.8857 for lag model
## ML residual variance (sigma squared): 2816.2, (sigma: 53.068)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 423.77, (AIC for lm: 421.98)
## LM test for residual autocorrelation
## test value: 3.1955, p-value: 0.073843

# 5. Model SEM (Spatial Error)
model_sem <- errorsarlm(equation, data = data_map, listw = listw_knn_dbd)
summary(model_sem)

## 
## Call:errorsarlm(formula = equation, data = data_map, listw = listw_knn_dbd)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -85.370 -32.498 -16.397  12.006 161.872 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)        58.9930891 93.4924840  0.6310   0.5280
## Akses_Sanitasi      0.6671698  0.9037136  0.7383   0.4604
## Kepadatan_Penduduk  0.0012603  0.0034800  0.3622   0.7172
## Jumlah_Rumah_Sakit -0.1017279  0.1069442 -0.9512   0.3415
## Penduduk_Miskin    -1.4325322  2.1926588 -0.6533   0.5135
## 
## Lambda: 0.0090664, LR test value: 0.0015099, p-value: 0.969
## Asymptotic standard error: 0.21134
##     z-value: 0.042899, p-value: 0.96578
## Wald statistic: 0.0018403, p-value: 0.96578
## 
## Log likelihood: -204.991 for error model
## ML residual variance (sigma squared): 2838.7, (sigma: 53.279)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 423.98, (AIC for lm: 421.98)

# 6. Pemilihan Model Terbaik
aic_vals <- AIC(model_ols, model_sar, model_sem)
print(aic_vals)

##           df      AIC
## model_ols  6 421.9835
## model_sar  7 423.7714
## model_sem  7 423.9820

best_mod <- rownames(aic_vals)[which.min(aic_vals$AIC)]
cat("\nModel terbaik berdasarkan AIC adalah:", best_mod, "\n")

## 
## Model terbaik berdasarkan AIC adalah: model_ols

Hasil pemodelan menunjukkan bahwa model terbaik (model dengan AIC terkecil) adalah model OLS, namun penyakit DBD ini memiliki korelasi spasial, sehingga model yang digunakan adalah SAR. Hal ini menunjukkan bahwa penyakit DBD erat kaitannya dengan faktor-faktor resiko yang ada.

4.5.1 Pemilihan Model Terbaik

Model Ordinary Least Squares (OLS) menghasilkan nilai terkecil yang menunjukkan kecocokan statistik yang baik, menurut hasil perbandingan nilai Akaike Information Criterion (AIC). Namun, keputusan untuk menggunakan model ini tidak hanya didasarkan pada metrik statistik,namun ada juga pertimbangan dari penyebaran Demam Berdarah Dengue (DBD) memiliki hubungan yang kuat antara wilayah. Karena itu, model Spatial Autoregressive (SAR) adalah model yang paling cocok untuk penelitian ini.

5 Kesimpulan dan Saran

5.1 Kesimpulan

Secara deskriptif, penyebaran kasus DBD di Indonesia terlihat tidak merata. Kawasan timur Indonesia cenderung memiliki insidensi yang lebih rendah, tetapi wilayah Kalimantan menonjol karena memiliki tingkat insidensi yang sangat tinggi, yang ditunjukkan pada peta kuartil. Analisis LISA juga menemukan adanya spatial outlier di sekitar Kalimantan dan Sulawesi, yang menandakan adanya provinsi dengan kasus rendah padahal dikelilingi oleh tetangga yang kasusnya tinggi. Kemudian hasil uji statistik menggunakan Indeks Moran’s I mengonfirmasi bahwa terdapat autokorelasi positif pada penyebaran DBD, meskipun nilainya tergolong lemah. Hal ini membuktikan bahwa tingginya kasus DBD di suatu provinsi tidak berdiri sendiri, melainkan ada pengaruh dari kondisi provinsi di sekitarnya.

Di bagian faktor resiko, ukuran asosiasi (Odds Ratio), sanitasi lingkungan terbukti memiliki hubungan dengan kejadian penyakit. Wilayah dengan akses sanitasi yang kurang layak memiliki risiko lebih besar untuk mengalami lonjakan kasus DBD dibandingkan wilayah dengan sanitasi yang baik.

Pada penentuan model terbaik, dipilih model Spatial Autoregressive (SAR). Meskipun nilai AIC model OLS sedikit lebih kecil, uji diagnostik menunjukkan bahwa OLS melanggar asumsi kebebasan sisaan. Akibatnya, model SAR dianggap lebih masuk akal secara statistik dan teoritis untuk menjelaskan penyebaran penyakit menular seperti DBD, yang dapat menyebar antarwilayah.

5.2 Saran

1. Bagi Pemerintah: Hasil pemetaan menunjukkan bahwa Kalimantan adalah daerah yang sangat rentan. Oleh karena itu, pemerintah harus memprioritaskan anggaran kesehatan dan program pencegahan, seperti fogging atau vaksinasi massal, serta meningkatkan koordinasi penanganan penyakit antarprovinsi.
2. Bagi Masyarakat: Hasil analisis menunjukkan bahwa tingkat sanitasi yang buruk berkorelasi positif dengan risiko terkena DBD. Oleh karena itu, kesadaran akan kebersihan lingkungan dan pemberantasan sarang nyamuk harus ditingkatkan (3M Plus).
3. Bagi Penelitian Selanjutnya: Agar pola penyebaran lebih jelas dan bias dikurangi, saya menyarankan agar penelitian berikutnya menggunakan unit analisis yang lebih kecil, seperti tingkat kabupaten atau kota. Sangat disarankan untuk menambah faktor iklim seperti suhu udara dan curah hujan karena faktor-faktor ini sangat memengaruhi perkembangbiakan nyamuk.

6 Referensi

Achmadi, U. F. (2009). Manajemen penyakit berbasis wilayah. Kesmas, 3(4), 147–153.

Ariyadi, T., Parasitologi, L., Keperawatan, I., Kesehatan, D., Mulia Ningsih, A., Santosa, B., Kesehatan, A., Universitas, K., Semarang, M., & Studi, P. (2020). UJI SKRINING FILARIASIS DI DESA JATIBARANG LOR KECAMATAN JATIBARANG KABUPATEN BREBES. Jurnal Labora Medika, 4, 20–24.

Beier, J. C. (1998). MALARIA PARASITE DEVELOPMENT IN MOSQUITOES. Annual Review of Entomology, 43(1), 519–543. https://doi.org/10.1146/annurev.ento.43.1.519

Er, A. C., Rosli, M. H., Asmahani, A., Mohamad Naim, M. R., & Harsuzilawati, M. (2010). Spatial mapping of dengue incidence: A case study in Hulu Langat District, Selangor, Malaysia. International Journal of Earth, Energy and Environmental Sciences, 4(7).

Kementerian Kesehatan RI. (2024). Profil Kesehatan Indonesia Tahun 2023. Jakarta: Kemenkes RI.

Loka Laboratorium Kesehatan Masyarakat Pangandaran. (2022, Juni 12). Segitiga penyakit tular vektor. https://labkesmaspangandaran.id/segitiga-penyakit-tular-vektor/

Septian Maksum, T., Amalan Tomia, Mk., & Ayu Rofia Nurfadillah, Ms. (n.d.). ENTOMOLOGI & PENGENDALIAN VEKTOR PENYAKIT.

Sobari, M., Jaya, I. G. N. M., & Ruchjana, B. N. (2023). Spatial Analysis of Dengue Disease in Jakarta Province. CAUCHY: Jurnal Matematika Murni Dan Aplikasi, 7(4), 535–547. https://doi.org/10.18860/ca.v7i4.17423

Tansil, M. G., Rampengan, N. H., & Wilar, R. (2021). Faktor Risiko Terjadinya Kejadian Demam Berdarah Dengue Pada Anak. Jurnal Biomedik:JBM, 13(1), 90. https://doi.org/10.35790/jbm.13.1.2021.31760