LAPORAN UJIAN TENGAH SEMESTER

Spatial Statistics

Disusun oleh:

Fara Dhiyaa Ramadhani (140610230066)
Andini Dwi Kurnia Putri (140610230082)
Jesslyn Alvina Gunawan (140610230085)

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

World Health Organization (WHO) menjelaskan bahwa Indonesia, sebagai negara kepulauan dengan kompleksitas geografis dan keragaman sosio-ekonomi yang tinggi, sedang menghadapi tantangan kesehatan masyarakat yang memerlukan pendekatan analitis yang mendalam. Di berbagai provinsi, beban penyakit menular kronis seperti malaria, demam berdarah dengue (DBD), dan filariasis masih menjadi masalah utama kesehatan masyarakat pada tahun 2024. Penyakit-penyakit tersebut merupakan penyakit yang ditularkan melalui vektor, dan penyebarannya tidak hanya dipengaruhi oleh faktor biologis, tetapi juga oleh faktor sosial, ekonomi, dan lingkungan yang saling berinteraksi.

Faktor-faktor seperti 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. Oleh karena itu, analisis hubungan spasial antara faktor-faktor tersebut dengan beban penyakit menjadi penting untuk memahami dinamika penularan di tingkat wilayah.

Dalam konteks tersebut, pendekatan One Health memiliki peran yang sangat penting. One Health menekankan adanya keterkaitan erat antara kesehatan manusia, hewan, dan lingkungan, dengan tujuan akhir untuk meningkatkan kesejahteraan manusia secara berkelanjutan melalui upaya pengendalian penyakit yang ditularkan oleh hewan maupun vektor. Pendekatan ini menuntut adanya kerja sama lintas sektor, melibatkan bidang kesehatan, pertanian, lingkungan, dan sosial, agar penanganan penyakit dapat dilakukan secara menyeluruh dan terpadu. Di Indonesia, penerapan prinsip One Health menjadi semakin relevan karena penyakit seperti malaria, DBD, dan filariasis sangat dipengaruhi oleh kondisi lingkungan, perilaku manusia, serta dinamika ekosistem vektor yang saling berinteraksi secara kompleks.

1.2 Identifikasi masalah

Apakah faktor sosio-ekonomi dan lingkungan/infrastruktur kesehatan (Sanitasi Layak, Kepadatan Penduduk, Jumlah Rumah Sakit, dan Kemiskinan) secara signifikan mempengaruhi beban penyakit menular kronis (malaria, DBD, dan filariasis) di tingkat provinsi?
Apakah terdapat ketergantungan spasial (antar-provinsi) dalam penyebaran penyakit menular kronis, dan seberapa kuat pengaruhnya?
Model statistik manakah yang paling akurat dan kuat dalam memodelkan hubungan ini untuk rekomendasi kebijakan?

1.3 Tujuan

Mendeskripsikan pola dan distribusi spasial kasus penyakit menular kronis (malaria, DBD, dan filariasis) di seluruh provinsi di Indonesia pada tahun 2024 melalui visualisasi peta tematik.
Menguji adanya autokorelasi spasial global dan lokal pada kasus penyakit menular kronis (malaria, DBD, dan filariasis).
Menganalisis dampak faktor lingkungan, sosio-ekonomi, dan infrastruktur kesehatan (Sanitasi Layak, Kepadatan Penduduk, Jumlah Rumah Sakit, dan Kemiskinan) terhadap kasus penyakit dengan mempertimbangkan efek spasial.
Mengidentifikasi model regresi spasial yang paling cocok untuk menjelaskan beban penyakit menular kronis di Indonesia melalui perbandingan kekuatan model.

1.4 Batasan Penelitian

Penelitian ini difokuskan pada kasus penyakit menular kronis (malaria, DBD, dan filariasis) yang menjadi beban utama masyarakat di tingkat provinsi di seluruh Indonesia pada tahun 2024.
Analisis spasial dilakukan pada tingkat Provinsi secara keseluruhan.
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 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.1.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.2 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.2.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.

2.3 Spatial Econometrics

Model Ekonometrik Spasial digunakan untuk menganalisis data antar-wilayah (seperti provinsi) di mana nilai suatu provinsi dipengaruhi oleh provinsi tetangganya (disebut dependensi spasial). Pendekatan ini digunakan ketika model regresi konvensional tidak dapat (Ordinary Least Squares) tidak mampu menangkap keterkaitan spasial yang terjadi. Ada tigamodel utama yaitu model SAR untuk melihat apakah jumlah kasus penyakit di suatu provinsi dipengaruhi langsung oleh jumlah kasus penyakit yang sama di provinsi tetangganya (efek penyebaran). Sementara itu, model SEM untuk mencari apakah ada faktor tersembunyi (seperti kesamaan iklim atau budaya) yang belum kita masukkan dalam model, namun memiliki pola spasial dan memengaruhi provinsi tetangga secara bersamaan. SDM memperluas SAR dengan memasukkan pengaruh variabel independen dari wilayah tetangga. Ketiga model ini biasanya diestimasi dengan metode Maximum Likelihood (MLE).

3 Metodologi

3.1 Sumber Data

Data yang digunakan 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.

3.3 Metode Analisis

Metode analisis yang dilakukan pada penelitian ini meliputi eksplorasi data secara spasial untuk melihat penyebaran data. Kemudian akan dilakukan perhitungan autokorelasi spasial menggunakan indeks seperti Morans I dan LISA. Selanjutnya, model-model spasial ekonometrik akan dibandingkan dan akan dipilih model terbaik paling cocok dan sesuai dengan data untuk memprediksi data.

3.4 Alur Kerja

Alur kerja yang akan dilakukan adalah :

Persiapan data
Uji pola spasial
Pemodelan Spatial Ekonometrik
Penentuan model terbaik
Interpretasi

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 Analisis Deskriptif

Untuk melihat pola dari data yang ada dilakukan analisis deskriptif dan pemetaan untuk masing-masing penyakit.

# ==============================================================
# 1. LOAD LIBRARIES
# ==============================================================
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

# ==============================================================
# 2. BACA DAN PERBAIKI 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)

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
  )

# Gabungkan shapefile dan data
data_map <- batas_prov %>%
  left_join(data, by = "Provinsi")

# ==============================================================
# 4. PERHITUNGAN NILAI k OPTIMAL UNTUK KNN (PER PENYAKIT)
# ==============================================================
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).

# ==============================================================
# 5. BUAT BOBOT SPASIAL KNN UNTUK SETIAP PENYAKIT
# ==============================================================
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 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.3 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, DBD, dan Filariasis. Hal ini menunjukkan bahwa penyakit DBD, Malaria, dan Filariasis dapat menyebar melalui wilayah.

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 pemanggilan untuk semua penyakit ===
# Pastikan 'weights_per_disease' berisi list per penyakit, misal:
# 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),

# 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 3.285410e-07 -0.02702703 0.007713138
## 2        DBD        NA 8.606237e-02 -0.02702703 0.013443977
## 3 Filariasis        NA 8.182231e-07 -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.

4.4 Pembentukan Model

library(spatialreg)

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

## Loading required package: Matrix

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

## 
## 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(spdep)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

library(car)

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

# Uji VIF antar semua X (pakai model dummy)
check_multicollinearity <- function(data, x_vars) {
  data_clean <- na.omit(st_drop_geometry(data)[, x_vars])
  
  # Buat model dengan satu variabel dummy
  formula_obj <- as.formula(paste("dummy ~", paste(x_vars, collapse = " + ")))
  data_clean$dummy <- rnorm(nrow(data_clean))  # variabel acak
  
  lm_model <- lm(formula_obj, data = data_clean)
  vif_values <- vif(lm_model)
  
  cat("📊 Hasil Uji Multikolinearitas (VIF):\n")
  print(vif_values)
  
  if (any(vif_values > 10)) {
    message("⚠️  Terdapat indikasi multikolinearitas (VIF > 10). Pertimbangkan menghapus variabel tersebut.")
  } else {
    message("✅ Tidak ada multikolinearitas serius (semua VIF < 10).")
  }
  
  return(vif_values)
}
  

# ===============================================================
# 🧮 Fungsi untuk membandingkan model spasial
# ===============================================================
fit_spatial_models <- function(data, disease_var, x_vars, weights) {
  data_clean <- st_drop_geometry(data)
  prevalence_var <- paste0(disease_var, "_prev")
  valid_idx <- !is.na(data_clean[[prevalence_var]]) & rowSums(is.na(data_clean[, x_vars])) == 0
  data_valid <- data_clean[valid_idx, ]
  
  if (nrow(data_valid) < 5) return(NULL)
  
  formula_str <- paste(prevalence_var, "~", paste(x_vars, collapse = " + "))
  formula_obj <- as.formula(formula_str)
  
  models <- list()
  
  # Model dasar OLS
  models$ols <- lm(formula_obj, data = data_valid)
  
  # Cek apakah bobot spasial valid
  if (!inherits(weights, "listw")) {
    stop("❌ Object 'weights' harus berupa objek listw hasil nb2listw().")
  }
  
  # Coba model spasial dengan perlindungan terhadap singularitas
  models$lag <- tryCatch(
    lagsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, tol.solve = 1e-10),
    error = function(e) NULL
  )
  models$error <- tryCatch(
    errorsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, tol.solve = 1e-10),
    error = function(e) NULL
  )
  models$durbin <- tryCatch(
    lagsarlm(formula_obj, data = data_valid, listw = weights, type = "mixed", zero.policy = TRUE, tol.solve = 1e-10),
    error = function(e) NULL
  )
  
  # Bandingkan berdasarkan AIC & BIC
  aic_values <- sapply(models, function(m) if (!is.null(m)) AIC(m) else NA)
  bic_values <- sapply(models, function(m) if (!is.null(m)) BIC(m) else NA)
  
  df <- data.frame(
    Penyakit = disease_var,
    Model = c("OLS", "Spatial Lag (SAR)", "Spatial Error (SEM)", "Spatial Durbin (SDM)"),
    AIC = aic_values,
    BIC = bic_values
  )
  
  return(df)
}

# ===============================================================
# ⚙️ Jalankan analisis
# ===============================================================
diseases <- c("Malaria", "DBD", "Filariasis")  
x_vars <- c("Akses_Sanitasi", "Kepadatan_Penduduk", "Jumlah_Rumah_Sakit", "Penduduk_Miskin")

# 1️⃣ Uji multikolinearitas
vif_values <- check_multicollinearity(data_map, x_vars)

## 📊 Hasil Uji Multikolinearitas (VIF):
##     Akses_Sanitasi Kepadatan_Penduduk Jumlah_Rumah_Sakit    Penduduk_Miskin 
##           2.416490           1.075299           1.098714           2.520005

## ✅ Tidak ada multikolinearitas serius (semua VIF < 10).

# 2️⃣ Jalankan model spasial
results <- lapply(diseases, function(d) fit_spatial_models(data_map, d, x_vars, weights))

## Warning in lagsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.19194e-16 - using numerical Hessian.

## Warning in errorsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.7155e-16 - using numerical Hessian.

## Warning in lagsarlm(formula_obj, data = data_valid, listw = weights, type = "mixed", : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.51383e-16 - using numerical Hessian.

## Warning in lagsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.42383e-11 - using numerical Hessian.

## Warning in errorsarlm(formula_obj, data = data_valid, listw = weights, zero.policy = TRUE, : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.44832e-11 - using numerical Hessian.

## Warning in lagsarlm(formula_obj, data = data_valid, listw = weights, type = "mixed", : inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 
##   reciprocal condition number = 2.28711e-11 - using numerical Hessian.

model_compare <- do.call(rbind, results)

print(model_compare)

##           Penyakit                Model      AIC      BIC
## ols        Malaria                  OLS 793.2026 803.0281
## lag        Malaria    Spatial Lag (SAR) 792.4623 803.9254
## error      Malaria  Spatial Error (SEM) 793.9999 805.4630
## durbin     Malaria Spatial Durbin (SDM) 795.2483 813.2618
## ols1           DBD                  OLS 421.9835 431.8091
## lag1           DBD    Spatial Lag (SAR) 423.7011 435.1642
## error1         DBD  Spatial Error (SEM) 423.9748 435.4379
## durbin1        DBD Spatial Durbin (SDM) 428.1433 446.1567
## ols2    Filariasis                  OLS 366.3758 376.2013
## lag2    Filariasis    Spatial Lag (SAR) 368.0396 379.5027
## error2  Filariasis  Spatial Error (SEM) 368.3621 379.8252
## durbin2 Filariasis Spatial Durbin (SDM) 371.8888 389.9023

Dari hasil pemodelan, model dengan AIC terkecil adalah model OLS tetapi dengan ketiga penyakit ini memiliki korelasi spasial sehingga model yang digunakan adalah model spasial ekonometrik. Sehingga untuk penyakit DBD, Malaria, dan Filariasis model spasial ekonometrik terbaik adalah model SAR (Spasial Autoregressive). Hal ini menunjukkan bahwa ketiga penyakit ini erat kaitannya dengan faktor-faktor risiko yang ada.

5 Kesimpulan dan Saran

5.1 Kesimpulan

Berdasarkan hasil uji autokorelasi spasial menggunakan Moran’s I dan Geary C didapatkan bahwa tiga penyakit (Malaria, DBD, dan Filariasis) memiliki autokorelasi spasial yang memiliki pola persebaran yang mengelompok antarwilayah. Temuan ini mengindikasikan bahwa kasus penyakit Malaria, DBD, dan Filariasis pada suatu wilayah dipengaruhi oleh kondisi wilayah-wilayah sekitarnya.

Pembentukan model ekonometrik spasial dilakukan untuk menjelaskan variasi prevalensi penyakit. Untuk penyakit DBD, Malaria, dan Filariasis model terbaiknya adalah SAR dan penyakit-penyakit ini dapat dijelaskan oleh lingkungan dan faktor-faktor risiko.

5.2 Saran

Implementasi Kebijakan Kesehatan Berbasis Klaster Spasial
Pemerintah dan dinas kesehatan perlu melakukan intervensi lintas-wilayah terutama pada provinsi yang teridentifikasi sebagai hotspot. Penanganan penyakit yang memiliki efek spasial langsung, seperti Malaria dan Filariasis, sebaiknya dilakukan dengan intervensi serentak antarwilayah yang saling berdekatan.

Fokus pada Peningkatan Sanitasi dan Pengentasan Kemiskinan
Karena sanitasi dan kemiskinan terbukti sebagai determinan paling signifikan, program peningkatan sanitasi lingkungan, pemenuhan air bersih, serta program pemberdayaan ekonomi perlu diprioritaskan untuk menurunkan risiko penyakit menular kronis secara berkelanjutan.

Optimalisasi Program Pencegahan DBD Berbasis Perubahan Perilaku
Mengingat faktor spasial pada DBD lebih banyak dipengaruhi oleh unsur tidak teramati, seperti perilaku masyarakat dan kebersihan lingkungan, maka edukasi dan pemberdayaan masyarakat melalui kampanye PSN, 3M Plus, dan gerakan kesehatan berbasis komunitas perlu ditingkatkan.

Penguatan Sistem Informasi dan Monitoring Spasial Penyakit
Pemerintah disarankan mengembangkan sistem pemantauan (surveillance) berbasis GIS dan dashboard spasial yang terintegrasi untuk memantau pergerakan penyakit secara real time dan memetakan wilayah prioritas intervensi secara cepat dan akurat.

Penelitian Lanjutan dengan Pendekatan Spatio-Temporal
Untuk mendapatkan pemahaman lebih mendalam, penelitian selanjutnya dapat menggunakan model Panel Data Spasial agar dapat menangkap dinamika spasial dan temporal secara bersamaan, serta menambahkan variabel lain seperti iklim, mobilitas penduduk, dan kualitas layanan kesehatan. :::::::::::

6 Referensi

Aurora, W. I. D., Maria, I., Kusdiyah, E., Darmawan, A., Nuriyah, & Mulyadi, D. (2023). Spatial Analysis of Filariasis and Malaria Determinants as Neglected Tropical Diseases in Indonesia (pp. 1010–1022). https://doi.org/10.2991/978-2-38476-110-4_97

Naim, S., Hastono, S. P., Rahayu, S., & Wangi, M. P. (2022). A Spatial Analysis of Dengue Hemorrhagic Fever (DHF), Hygiene, and Latrines in Depok City in 2020. JURNAL KESEHATAN LINGKUNGAN, 14(2), 122–129. https://doi.org/10.20473/jkl.v14i2.2022.122-129

Rahmawati, Y., Jamil, I. R., & Hidayah, I. (2025). A Spatial Analysis on Heterogenous Determinant of Dengue Fever Cases in Indonesia. Journal of Geovisualization and Spatial Analysis, 9(1), 10. https://doi.org/10.1007/s41651-024-00212-1

Schmidt, W.-P., Suzuki, M., Dinh Thiem, V., White, R. G., Tsuzuki, A., Yoshida, L.-M., Yanai, H., Haque, U., Huu Tho, L., Anh, D. D., & Ariyoshi, K. (2011). Population Density, Water Supply, and the Risk of Dengue Fever in Vietnam: Cohort Study and Spatial Analysis. PLoS Medicine, 8(8), e1001082. https://doi.org/10.1371/journal.pmed.1001082

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

Wang, W., Xiao, X., Qian, J., Chen, S., Liao, F., Yin, F., Zhang, T., Li, X., & Ma, Y. (2022). Reclaiming independence in spatial‐clustering datasets: A series of data‐driven spatial weights matrices. Statistics in Medicine, 41(15), 2939–2956. https://doi.org/10.1002/sim.9395

Widawati, M., Ipa, M., Astuti, E. P., Wahono, T., & Yuliasih, Y. (2022). The Activities on Prevention of Malaria and Filariasis Vector Bites among Indonesian Society: A Nationwide Disease Prevention Survey. Indonesian Journal of Tropical and Infectious Disease, 10(2), 104–112. https://doi.org/10.20473/ijtid.v10i2.36053

Zeanova, H., Taniwan, P., Aulia Putri, R., & Yusti Faidah, D. (2024). ANALISIS FAKTOR PENYEBAB PENYAKIT TBC DI JAWA BARAT MENGGUNAKAN REGRESI BINOMIAL NEGATIF. Jurnal Lebesgue : Jurnal Ilmiah Pendidikan Matematika, Matematika Dan Statistika, 5(3), 2284–2302. https://doi.org/10.46306/lb.v5i3.859

Zebua, H. I., & Jaya, I. G. N. M. (2022). Spatial Autoregressive Model of Tuberculosis Cases in Central Java Province 2019. CAUCHY: Jurnal Matematika Murni Dan Aplikasi, 7(2), 240–248. https://doi.org/10.18860/ca.v7i2.13451

Analisis Spasial pada Penyakit DBD, Malaria, dan Filariasis