LAPORAN UJIAN TENGAH SEMESTER

Epidemiologi

Disusun oleh:

Fara Dhiyaa Ramadhani (140610230066)

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

Penyakit menular melalui hewan seperti Demam Berdarah Dengue (DBD), malaria, dan filariasis menjadi salah satu masalah kesehatan masyarakat di Indonesia sebagai negara tropis. Ketiga penyakit ini disebabka oleh vektor nyamuk menimbulkan morbiditas dan mortalitas yang tinggi sehingga penting untuk dilakukan berbagai upaya pengendalian.

Penularan ketiga penyakit ini sangat dipengaruhi oleh faktor lingkungan dan faktor sosial-ekonomi. Menjadi sangat penting bagaimana masyarakat mendapat akses terhadap sanitasi layak. Lingkungan dengan sanitasi buruk, dengan genangan air, menjadi tempat berkembang biak vektor yang menjadi transmitor Aedes aegypti, Culex sp., dan Anopheles sp..

Selain itu, dengan kepadatan penduduk yang tinggi dapat meningkatkan kontak antara manusia dan vektor, serta memperluaas wilayah penyebaran kasus. Jumlah fasilitas kesehatan seperti rumah sakit menjadi salah satu hal yang dibutuhkan untuk melakukan deteksi dini, pengobatan, dan pengendalian penyakit-penyakit tersebut. Daerah dengan keterbatasan rumah sakit seringkali menghadapi keterlambatan diagnosis dan penangan kasus.

Faktor kemiskinan juga memiliki peran besar, masyarakat dengan tingkat kesejahteraan rendah umumnya memiliki akses terbatas terhadap air bersih, perumahan layak, dan layanan kesehatan, sehingga lebih rentan terhadap infeksi penyakit.

1.2 Identifikasi Masalah

Bagaimana pola epidemiologi dan distribusi spasial penyakit tular vektor (DBD, filariasis, dan malaria) di Indonesia pada tingkat provinisi tahun 2024
Apakah terdapat autokorelasi spasial antar-provinisi dalam penyebaran penyakit DBD, Filariasis, dan Malaria di Indonesia
Sejauh mana faktor akses sanitasi layak, kepadatan penduduk, jumlah rumah sakit, dan tingkat kemiskinan berpengaruh terhadap kejadian penyakit.

1.3 Tujuan

Mendeskripsikan pola epidemiologi dan distribusi spasial penyakit DBD, malaria, dan filariasis di tingkat provinsi pada tahun 2024
Menghitung ukuran-ukuran epidemiologi untuk masing-masing penyakit (ukuran frekuensi dan asosiasi)
Mengidentifikasi adanya autokorelasi spasial global dan lokal di seluruh provinsi Indonesia
Menganalisis pengaruh faktor-faktor terhadap kejadian penyakit DBD, malaria, dan filariasis
Memberikan interpretasi epidemiologis terhadap hasil visualisasi spasial dan ukuran asosiasi, serta menyusun rekomendasi kebijakan untuk pengendalian dan pencegahan penyakit menular tersebut.

1.4 Batasan Penelitian

Penelitian hanya terfokus pada tiga penyakit menular : DBD, Malaria, Filariasis.
Unit analisis adalah provinsi di Indonesia tahun 2024
Variabel dependen utama adalah beban penyakit menular kronis untuk malaria, DBD, dan filariasis.
Variabel independen (penjelas) dibatasi pada empat faktor utama: Persentase Rumah Tangga yang Memiliki Akses terhadap Sanitasi Layak, Kepadatan Penduduk per-km, Jumlah Desa yang Memiliki Rumah Sakit, dan Persentase Penduduk Miskin.
Model spasial menggunakan matriks bobot spasial KNN (K- Nearest Neighbor)

2 Tinjauan Pustaka

2.1 Penyakit Tular Vektor (Vector-borne diseases)

Penyakit Tular Vektor adalah penyakit yang ditularkan melalui vektor dan hewan pembawa penyakit. Salah satu vektor yang menjadi perantara biologis kepada manusia adalah nyamuk. Penyakit yang menjdai fokus pada penelitian ini adalah deman berdarah dengue (DBD), malaria, dan filariasis. Ketiga penyakit ini menular melalui vektor nyamuk dan menurut Kemenkes Labkesmas Pangandaran penyakit tular vektor ini masih menjadi permasalahan kesehatan di Indonesia.

2.2 Faktor Agent-Host-Environment

Agent merupakan micro-organisme atau pathogen seperti virus, bakteri, parasit, dsb yang terinfeksi dan menyebabkan penyakit muncul. Host merupakan manusia atau makhluk hidup yang dapat diinfeksi oleh agen. Environment adalah semua faktor eksternal yang mempengaruhi peluang terjadinya kontak antara host dan agen.

Pada penelitian ini, masing-masing penyakit memiliki agen nya tersendiri. Pada penyakit demam berdarah dengue (DBD) agennya adalah virus Flavivirus yang terbawa oleh vektor nyamuk Aedes Aegypti dan beberapa nyamuk lainnya dari genus Aedes. Penyakit Malaria menular melalui parasit Plasmodium yang terbawa oleh nyamuk Anopheles betina. Filariasis limfatik atau kaki gajah menular melalui infeksi cacing filaria yang dibawa oleh nyamuk Culex, Anopheles, Mansonia, Aedes.

Lingkungan pada penelitian ini menjadi faktor risiko dari penyakit. Faktor risikonya pada penelitian ini persentase rumah tangga dengan akses terhadap sanitasi layak, kepadatan penduduk per km\(^2\) , jumlah desa yang memiliki fasilitas rumah sakit, serta persentase penduduk miskin, berperan penting dalam menentukan tingkat risiko dan penyebaran penyakit menular. Sanitasi yang tidak memadai meningkatkan risiko paparan terhadap agen penyakit berbasis lingkungan; kepadatan penduduk yang tinggi mempercepat penularan antarindividu; keterbatasan fasilitas kesehatan memperlambat deteksi dan penanganan kasus; sementara kemiskinan meningkatkan kerentanan masyarakat terhadap infeksi.

**Gambar 1.** Grafik Hubungan Agent-Host-Environment
*Sumber:* Microbe Notes

2.3 Desain Studi

Desain studi dalam konteks epidemiologi adalah metodologi utama untuk menguji hipotesis tentang hubungan antara paparan dan penyakit. Secara umum, studi dibagi 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 tertentu, serta kemungkinan bias (seperti bias mengingat atau kehilangan mengingat) yang dapat memengaruhi validitas kesimpulan kausalitas. Oleh karena itu, pemilihan desain yang tepat sangat penting.

2.4 Ukuran Frekuensi Epidemiologi

2.4.1 Incidence Rate (Insidensi)

Insidensi merupakan jumlah kasus yang muncul selama periode waktu tertentu dibagai populasi.

\[ \text{Insidensi} = \frac{\text{Jumlah Kasus Baru}}{\text{Populasi}} \times k \]

k adalah perbandingan jumlah orang yang diinginkan.

2.4.2 Prevalensi

Prevalensi merupakan proporsi individu dalam populasi yang terkena penyakit pada rentang waktu tertentu

\[ \text{Prevalensi}= \frac{\text{Jumlah Total Kasus}}{\text{Jumlah Populasi}} \]

Untuk nilai prevalensi dibutuhkan total kasus pada periode penelitian dilakukan (tidak hanya kasus baru).

2.4.3 Attack Rate

Attack rate merupakan proporsi populasi yang terkena penyakit selama adanya kejadian wabah/epidemik singkat.

\[ \text{Attack Rate} = \frac{\text{Jumlah kasus selama wabah}}{\text{Populasi}}\times 100\% \]

2.4.4 Case Fatality Rate (CFR)

CFR adalah proporsi kasus yang meninggal dari keseluruhan jumlah kasus

\[ CFR = \frac{\text{Jumlah kematian akibat penyakit}}{\text{Jumlah kasus penyakit}} \times 100\% \]

2.4.5 Mortality Rate

Mortality Rate berbeda dikit dengan CFR, kalau CFR merupakan proporsi tetapi mortality rate adalah jumlah kematian per populasi

2.5 Ukuran Asosiasi Epidemiologi

Ukuran asosiasi adlah alat untuk mengukur kekuatan hubungan antara faktor paparan dan hasil penyakit. Ukuran asosiasi membanidngkan kelompok terpapar dan kelompok tidak terpapar.

2.5.1 Risk Ratio

Risk ratio merupakan perbandingan insidensi antara kelompok terpapar dengan tidak terapar.

\[ RR = \frac{Insidensi_{terpapar}}{Insidensi_{tidak terpapar}} \]

2.5.2 Odds Ratio

Odds ratio merupakan perbandingan odds paparan antara kasus dan kontrol. Odds ratio dapat dihitung ketika data dapat disajikan dalam tabel frekuensi dan sering umumnya digunakan di studi case-control.

2.5.3 Attributable Risk

AR mengukur seberapa besar proporsi kasus pada kelompok terpapar yang dapat dihindarkan jika paparan dicegah.

2.6 Spatial Dependence

Dalam statistik klasik, asumsi independensi observasi tidak berlaku, namun, dependensi spasial terjadi ketika nilai variabel di suatu provinsi dipengaruhi oleh variabel yang sama di provinsi tetangganya. Dalam kasus penyakit menular kronis seperti malaria, diare, dan filariasis, ketergantungan ini dapat disebabkan oleh penyebaran penyakit antar-wilayah (spatial lag) atau kesamaan faktor sosio-ekonomi dan lingkungan yang tidak teramati di daerah yang berdekatan. Pemodelan ini penting karena wilayah dengan beban penyakit tinggi cenderung membentuk kluster geografis, yang memerlukan pengujian autokorelasi spasial untuk mengetahui kekuatan dan karakteristik penyebarannya antar-provinsi.

2.6.1 Autokorelasi Spasial

Salah satu ukuran statistik yang dikenal sebagai autokorelasi spasial menunjukkan seberapa sistematis nilai suatu fenomena di suatu lokasi bergantung pada nilai fenomena yang sama di lokasi-lokasi sekitarnya. Kondisi ini dapat bersifat positif, yaitu daerah di sekitarnya cenderung memiliki nilai variabel yang sebanding (kluster nilai tinggi atau rendah), atau negatif, yaitu daerah di sekitarnya cenderung memiliki nilai variabel yang negatif. Untuk mempelajari pola ini, digunakan metrik global seperti Indeks Moran’s I untuk menguji keberadaan autokorelasi spasial umum di seluruh wilayah, dan metrik lokal seperti Local Indicators of Spatial Association (LISA), yang termasuk Local Moran’s I untuk menentukan lokasi kluster (misalnya, tinggi atau rendah) dan outlier spasial.

2.7 Matriks Bobot Spasial

Matriks bobot spasial merupakan komponen penting dalam merepresentasikan hubungan ketergantungan antarwilayah atau antarunit observasi. Matriks bobot spasisal berfungsi untuk menunjukkan sejauh mana suatu wilayah dipengaruhi oleh wilayah lain berdasarkan kedekatan geografis, jarak, atau hal lainnya. Bobot ini dapat ditentukan berdasarkan kedekatan fisik, jarak geografis, atau tetangga terdekat, hal ini dapat dipilih tergantung penelitian dan karakteristik data.

2.7.1 K-Nearest Neighbor (KNN)

Metode K-Nearest Neighbor (KNN) merupakan pendekatan pembobotan spasial berdasarkan jarak Euclidean sejumlah tetangga terdekat (k). Dengan metode ini, setiap wilayah akan memiliki jumlah hubungan spasial yang sama sehingga menghasilkan matriks bobot yang lebih seimbang.

3 Metodologi

3.1 Sumber Data

Data yang digunakan sebagai populasi dalam penelitian ini berasal dari laporan resmi Kementerian Kesehatan Republik Indonesia, terutama Profil Kesehatan Indonesia selama periode tahun 2024. Data yang digunakan sebagai unit observasi dalam penelitian mencakup semua 38 provinsi di Indonesia. Variabel-variabel yang digunakan adalah sebagai berikut:

data <-read.csv("C:/KAMPUS CHECK!/SPASIAL/DATA SPASIAL EPIDEM - VARIABEL DEPENDEN.csv")
library(dplyr)

## 
## 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

Variabel Malaria, DBD, dan Filariasis merupakan jumlah kasus yang terjadi di masing-masing provinsi. Variabel penyakit ini akan menjadi variabel dependen yang masing-masing penyakit ini akan dibuatkan model dan dilakukan mapping. Kemudian 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.

3.2 Unit Spasial

Unit spasial (spatial unit) yang digunakan dalam penelitian ini adalah 38 Provinsi di Indonesia. Setiap provinsi dilihat sebagai satuan geografis yang mewakili satu observasi dalam analisis spasial. Data kasus yang digunakan berasal dari Kementerian Kesehatan dan Badan Pusat Statistik (BPS) yang tersedia dan disagregasi pada tingkat administrasi.

4 Hasil dan Pembahasan

4.1 Dashboard

Berikut adalah link untuk dashboard yang telah dibuat : https://spatial-analysis-for-zoonosis-diseases-in-indonesia.streamlit.app/

Terdapat beberapa keterbatasan, salah satunya peta yang dapat tersajikan masih yang 34 provinsi Indonesia.

4.2 Ukuran Epidemiologi

Pada ukuran epidemiologi ini, karena ada keterbatasan pada akses data sehingga ukuran epidemiologi yang dapat dibentuk hanyalah prevalensi. Data hanya terbatas pada total jumlah kasus pada tahun 2024 dan jumlah penduduk pada masing-masing provinsi.

4.2.1 Prevalensi

options(scipen=10)
data <- data %>%
  mutate(
    Jumlah_Penduduk = Jumlah_Penduduk * 1000,
    Malaria_prev     = (Malaria / Jumlah_Penduduk) * 100000,
    DBD_prev        = (DBD / Jumlah_Penduduk) * 100000,
    Filariasis_prev  = (Filariasis / Jumlah_Penduduk) * 100000
  )
prev <- data %>%
  select(Provinsi, DBD_prev, Malaria_prev, Filariasis_prev)
prev

Nilai prevalensi ini menunjukkan berapa banyak kasus yang sekiranya akan ada pada 100000 orang di masing-masing provinsi. Semakin tinggi prevalensi pada sebuah wilayah menunjukkan bahwa penyakit tersebut semakin banyak dan semakin berpotensi untuk menular pada wilayah tersebut. Selain itu setiap perhitung selanjutnya akan digunakan nilai prevalensi karena jika menggunakan jumlah kasus ada kemungkinan beberapa wilayah yang memiliki jumlah populasi yang tinggi akan tinggi pula, padahal sebenarnya insidensinya rendah.

4.3 Analisis Deskriptif Spasial

library(sf)

## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE

library(dplyr)
library(spdep)

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

## Loading required package: spData

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

## 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.5.1

## Loading required package: viridisLite

library(readr)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.3.0
## ✔ purrr     1.0.4     ✔ 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

# Membaca data spasial
batas_prov <- readRDS("C:/KAMPUS CHECK!/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)

# 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

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)
}

# Hitung untuk setiap penyakit
k_malaria    <- find_optimal_k(data_map$Malaria, "Malaria")

## 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 Malaria = 1, diganti menjadi 3 (default).

k_dbd        <- find_optimal_k(data_map$DBD, "DBD")

## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 12
## sub-graphs
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points

## ✅ k optimal untuk DBD =3

k_filariasis <- find_optimal_k(data_map$Filariasis, "Filariasis")

## Warning in knn2nb(knearneigh(coords, k = k)): neighbour object has 12
## sub-graphs
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points
## Warning in knn2nb(knearneigh(coords, k = k)): k greater than one-third of the
## number of data points

## ⚠️  k optimal untuk Filariasis = 1, diganti menjadi 3 (default).

nb_knn_malaria    <- knn2nb(knearneigh(coords, k = k_malaria))
listw_knn_malaria <- nb2listw(nb_knn_malaria, style = "W", zero.policy = TRUE)

nb_knn_dbd        <- knn2nb(knearneigh(coords, k = k_dbd))
listw_knn_dbd     <- nb2listw(nb_knn_dbd, style = "W", zero.policy = TRUE)

nb_knn_filariasis <- knn2nb(knearneigh(coords, k = k_filariasis))
listw_knn_filariasis <- nb2listw(nb_knn_filariasis, style = "W", zero.policy = TRUE)
listw <- nb2listw(nb_knn_dbd, style = "W", zero.policy = TRUE)
weights <- listw_knn_dbd


cat("✅ Bobot spasial KNN telah dibuat untuk ketiga penyakit.\n")

## ✅ Bobot spasial KNN telah dibuat untuk ketiga penyakit.

cat("   - Malaria: k =", k_malaria, "\n")

##    - Malaria: k = 3

cat("   - DBD: k =", k_dbd, "\n")

##    - DBD: k = 3

cat("   - Filariasis: k =", k_filariasis, "\n")

##    - Filariasis: k = 3

# ==============================================================
# 6. PEMETAAN GRADASI PENYAKIT
# ==============================================================
diseases <- c("Malaria", "DBD", "Filariasis")

for (d in diseases) {
  prev_col <- paste0(d, "_prev")
  
  p <- ggplot(data_map) +
    geom_sf(aes(fill = .data[[prev_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
  )
}

# ==============================================================
# 6. PEMETAAN PENYAKIT BERDASARKAN KUARTIL + CETAK NILAI KUANTIL
# ==============================================================

diseases <- c("Malaria", "DBD", "Filariasis")

for (d in diseases) {
  prev_col <- paste0(d, "_prev")
  
  # Hitung nilai kuantil (Q0, Q1, Q2, Q3, Q4)
  quantile_values <- quantile(
    data_map[[prev_col]],
    probs = seq(0, 1, 0.25),
    na.rm = TRUE
  )
  
  # Cetak ke konsol
  cat("\n==============================\n")
  cat("📊 Nilai Kuantil Prevalensi -", d, "\n")
  print(quantile_values)
  cat("==============================\n\n")
  
  # Membuat kategori kuartil (Q1 - Q4)
  data_map <- data_map %>%
    mutate(
      quartile = cut(
        .data[[prev_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 Prevalensi - Malaria 
##           0%          25%          50%          75%         100% 
##     0.142032     1.761651    16.820487    45.138604 39707.274927 
## ==============================

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

## 
## ==============================
## 📊 Nilai Kuantil Prevalensi - Filariasis 
##          0%         25%         50%         75%        100% 
##   0.0000000   0.3847621   1.5811700   4.9331521 176.7755314 
## ==============================

4.4 Uji Autokorelasi Spasial

# Daftar penyakit yang diuji
diseases <- c("Malaria", "DBD", "Filariasis")

# 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) {
  prev_col <- paste0(d, "_prev")
  values <- data_map[[prev_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)
      )
    )
  }
}

# ===============================
# 5. TAMPILKAN HASIL
# ===============================
print(autocorr_results)

##                      Penyakit Moran_I Moran_p Geary_C Geary_p
## Moran I statistic     Malaria  0.4098  0.0000  0.6305  0.0120
## Moran I statistic1        DBD  0.1313  0.0861  0.8439  0.1091
## Moran I statistic2 Filariasis  0.2797  0.0000  0.7173  0.0613

Dari hasil pengujian autokorelasi spasial didapatkan bahwa yang memang memiliki autokorelasi spasial adalah penyakit Malaria dan Filariasis sedangkan penyakit DBD tidak secara signifikan memiliki autokorelasi spasial atau autokorelasinya sangat lemah.

calculate_morans_i <- function(data, disease_var, weights_list) {
  prevalence_var <- paste0(disease_var, "_prev")
  prevalence_values <- st_drop_geometry(data)[[prevalence_var]]
  valid_idx <- !is.na(prevalence_values)
  prevalence_values <- prevalence_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(prevalence_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(prevalence_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)
  prevalence_var <- paste0(disease_var, "_prev")
  prevalence <- data_df[[prevalence_var]]

  valid_idx <- !is.na(prevalence)
  prevalence_valid <- prevalence[valid_idx]

  if (length(prevalence_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(prevalence_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$Prevalence <- prevalence_valid
  data_subset$Ii <- lisa[, "Ii"]
  data_subset$p_value <- p_values
  data_subset$Prevalence_std <- as.numeric(scale(prevalence_valid)[, 1])
  lag_vals <- lag.listw(weights_subset, prevalence_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$Prevalence_std > 0 & data_subset$Lag_std > 0] <- "HH (Hotspot)"
  data_subset$Cluster[sig_mask & data_subset$Prevalence_std < 0 & data_subset$Lag_std < 0] <- "LL (Coldspot)"
  data_subset$Cluster[sig_mask & data_subset$Prevalence_std > 0 & data_subset$Lag_std < 0] <- "HL (Outlier)"
  data_subset$Cluster[sig_mask & data_subset$Prevalence_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(
  Malaria = list(neighbors = nb_knn_malaria, weights = listw_knn_malaria),
  DBD = list(neighbors = nb_knn_dbd, weights = listw_knn_dbd),
  Filariasis = list(neighbors = nb_knn_filariasis, weights = listw_knn_filariasis)
)

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    Malaria        NA 0.0000003285410 -0.02702703 0.007713138
## 2        DBD        NA 0.0860623657305 -0.02702703 0.013443977
## 3 Filariasis        NA 0.0000008182231 -0.02702703 0.004095232

# --- Local Moran’s I (LISA) per penyakit (kembali sebagai list of sf / NULL)
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 prevalensi penyakit",
      fill = "Cluster"
    ) +
    theme_minimal()

  print(p)
}

Dari LISA Maps menunjukkan pada penyakit Malaria terdapat hotspot pada Pulau Papua, untuk penyakit DBD terdapat outlier pada sekitar pulau Kalimantan, dan untuk penyakit Filariasis menunjukkan adanya hotspot juga di Pulau Papua.

5 Kesimpulan

Berdasarkan hasil analisis epidemiologi dan spasial, didapatkan bahwa distribusi dari penyakit DBD cukup merata pada provinsi-provinsi di Indonesia tetapi yang paling banyak adalah di Kalimantan, untuk penyakit Malaria penyebaran penyakit cenderung banyak di Pulau Papua dan begitupun penyakit Filariasis. Hasil analisis autokorelasi spasial global menunjukkan bahwa penyakit Malaria dan FIlariasis memiliki hubungan antar wilayahnya sehingga ketika satu provinsi memiliki tingkat penyakit yang tinggi maka daerah sekitarnya juga tinggi. Namun hal ini tidak dengan penyakit DBD.

Hasil penelitian ini menegaskan bahwa peningkatan akses sanitasi, pengurangan kemiskinan, pemerataan akses fasilitas sehatan menjadi kunci untuk menurunkan risiko penyakit DBD, malaria, dan filariasis. Provinisi-provinisi yang menjadi hotspot harus diperhatian agar penyakit menular ini tidak semakin menyebar bahkan ke wilayah lain.

6 Referensi

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

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.

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

Analisis Epidemiologi pada Penyakit DBD, Malaria, dan Filariasis