LAPORAN UJIAN AKHIR SEMESTER

Spatial Statistics






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 negara tropis seperti Indonesia, demam berdarah dengue (DBD) masih menjadi salah satu masalah epidemiologi terbesar. Kondisi iklim dengan curah hujan tinggi serta kelembapan yang fluktuatif menciptakan ekosistem yang ideal bagi perkembangbiakan vektor nyamuk Aedes aegypti dan Aedes albopictus. Insidensi DBD di Indonesia terus menunjukkan pola yang berubah-ubah dan sering menyebabkan kejadian luar biasa (KLB) di berbagai provinsi, meskipun berbagai upaya pengendalian telah dilakukan. Data dari Kementerian Kesehatan Republik Indonesia menunjukkan bahwa kasus DBD meningkatkan morbiditas (angka kesakitan) dan kematian (angka kematian). Selain itu, kasus ini menimbulkan tantangan ekonomi yang signifikan bagi sistem kesehatan nasional.

Penyebaran kasus DBD tidak terjadi secara acak, melainkan memiliki pola ketergantungan antarwilayah yang dipengaruhi oleh determinan sosial-ekonomi dan lingkungan. Faktor-faktor seperti kepadatan penduduk yang tinggi mempermudah transmisi virus antarmanusia, sementara persentase rumah tangga dengan akses sanitasi yang buruk berkontribusi pada ketersediaan tempat perindukan nyamuk (breeding places). Selain itu, kemiskinan dan keterbatasan akses terhadap fasilitas kesehatan sering kali memperburuk risiko penyebaran penyakit ini. Mengingat sifat vektor yang mobile dan interaksi antarmanusia yang melintasi batas administrasi, kejadian DBD di suatu wilayah sangat mungkin dipengaruhi oleh kondisi wilayah tetangganya, sebuah fenomena yang dikenal sebagai efek dependensi spasial.

Oleh karena itu, pendekatan statistika konvensional seperti regresi OLS (Ordinary Least Squares) sering kali kurang memadai karena mengabaikan efek spasial tersebut. Penelitian ini menerapkan analisis spasial komprehensif menggunakan uji autokorelasi spasial (Moran’s I dan LISA) untuk mendeteksi klaster penyebaran (hotspot), serta pemodelan Spatial Econometrics (SAR, SEM, SDM) untuk mengestimasi pengaruh variabel prediktor secara lebih akurat. Selain itu, metode interpolasi spasial seperti Inverse Distance Weighting (IDW) dan Ordinary Kriging digunakan untuk memprediksi pola risiko di area yang tidak teramati. Pendekatan ini diharapkan dapat memberikan peta risiko yang presisi sebagai dasar pengambilan kebijakan berbasis wilayah dalam pengendalian DBD di Indonesia.

1.2 Identifikasi Masalah

  • Bagaimana pola persebaran spasial kasus Demam Berdarah Dengue (DBD) di Indonesia, dan apakah terdapat autokorelasi spasial serta klaster wilayah (hotspots) yang signifikan?
  • Bagaimana pengaruh faktor sosio-ekonomi dan lingkungan terhadap prevalensi DBD ketika efek dependensi spasial diperhitungkan?
  • Model ekonometrika spasial manakah yang paling akurat dan memenuhi kriteria uji statistik dalam memodelkan penyebaran DBD di Indonesia?
  • Bagaimana gambaran prediksi risiko penyebaran DBD di seluruh wilayah Indonesia, termasuk area yang tidak teramati, berdasarkan pemodelan interpolasi spasial?

1.3 Tujuan

  • Mendeskripsikan pola distribusi spasial dan tingkat kerawanan Demam Berdarah Dengue (DBD) di tingkat provinsi di Indonesia melalui visualisasi peta tematik kuartil.
  • Membuktikan adanya autokorelasi spasial global menggunakan indeks Moran’s I serta mengidentifikasi lokasi klaster spasial secara lokal menggunakan Local Indicator of Spatial Association (LISA) pada kasus DBD.
  • Menganalisis signifikansi pengaruh faktor sanitasi layak, kepadatan penduduk, ketersediaan fasilitas kesehatan (jumlah rumah sakit), dan kemiskinan terhadap prevalensi DBD dengan menerapkan model Ekonometrika Spasial.
  • Mengevaluasi kinerja model regresi spasial untuk menentukan model terbaik yang paling presisi dalam menjelaskan fenomena penyebaran DBD di Indonesia.
  • Memetakan prediksi risiko penyebaran DBD di seluruh wilayah Indonesia menggunakan metode interpolasi spasial Inverse Distance Weighting (IDW) dan Ordinary Kriging sebagai dasar mitigasi di area yang belum terjangkau data.

1.4 Batasan Penelitian

  • Penelitian ini dibatasi secara spesifik pada analisis spasial kasus Demam Berdarah Dengue (DBD) di Indonesia tahun 2024
  • Unit analisis spasial yang digunakan adalah level Provinsi di seluruh Indonesia. Analisis interpolasi menggunakan titik tengah (centroid) dari poligon provinsi sebagai representasi lokasi geografis.
  • Variabel dependen utama adalah beban penyakit menular kronis untuk DBD
  • 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)
  • Pemodelan regresi spasial dibatasi pada perbandingan tiga model utama dalam Spatial Econometrics, yaitu Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM).

2 Tinjauan Pustaka

2.1 Spatial Dependence

Hukum Pertama Geografi yang dikemukakan oleh Tobler (1970) menyatakan bahwa “segala sesuatu saling berhubungan dengan hal lainnya, tetapi hal-hal yang dekat lebih berhubungan daripada hal-hal yang jauh.” Konsep ini mendasari asumsi spatial dependence (ketergantungan spasial). Dalam data kejadian penyakit, asumsi independensi antar-observasi yang digunakan dalam statistik klasik sering kali terlanggar karena adanya efek tumpahan (spillover effect) antarwilayah.

2.1.1 Autokorelasi Spasial

Autokorelasi spasial mengukur korelasi antara nilai suatu variabel dengan nilai variabel yang sama di lokasi sekitarnya. Moran’s I digunakan untuk mendeteksi pola spasial secara menyeluruh di area studi. Nilai positif mengindikasikan pengelompokan (kluster) wilayah dengan nilai serupa, sedangkan nilai negatif mengindikasikan pola penyebaran (dispersed). LISA merupakan dekomposisi lokal dari Moran’s I yang berfungsi untuk mengidentifikasi lokasi spesifik dari kluster spasial signifikan, seperti High-High (Hotspot), Low-Low (Coldspot), serta outlier spasial (High-Low dan Low-High).

2.2 Matriks Bobot Spasial

Untuk mengoperasionalkan hubungan antarwilayah, diperlukan matriks bobot spasial (\(W\)). Salah satu metode pembobotan adalah K-Nearest Neighbor (KNN). Berbeda dengan metode contiguity (persinggungan batas) yang bisa bermasalah pada wilayah kepulauan yang terisolasi, KNN mendefinisikan tetangga berdasarkan jarak titik pusat (centroid) terdekat sejumlah \(k\). Metode ini menjamin setiap wilayah memiliki jumlah tetangga yang sama, sehingga struktur ketergantungan spasial lebih seimbang dan tidak ada wilayah yang dianggap terisolasi.

2.3 Spatial Econometrics

Ketika regresi linier biasa (Ordinary Least Squares/OLS) menghasilkan residu dengan autokorelasi spasial, model ekonometrika spasial digunakan. Tiga model utama yang paling umum digunakan adalah: - Spatial Autoregressive Model (SAR) mengatakan bahwa nilai variabel dependen di suatu wilayah dipengaruhi oleh lag spasial, atau nilai variabel dependen di wilayah tetangganya. Efek difusi penyakit atau penularan antarwilayah adalah topik yang cocok untuk model ini. - Spatial Error Model (SEM): Model ini menganggap bahwa ada ketergantungan spasial pada komponen kesalahan (galat) karena variabel independen yang memiliki pola spasial (seperti faktor iklim atau budaya lokal) yang tidak dimasukkan ke dalam model. - Spatial Durbin Model (SDM): Model yang lebih umum yang memasukkan lag spasial pada variabel dependen (Y) dan variabel independen (X). Model ini memiliki kemampuan untuk mengidentifikasi bagaimana kondisi sosio-ekonomi wilayah tetangga berdampak pada jumlah kasus penyakit yang terjadi di daerah amatan.

2.4 Interpolasi Spasial

Interpolasi spasial adalah metode untuk menaksir nilai pada lokasi yang tidak tersampel berdasarkan nilai-nilai dari lokasi yang tersampel di sekitarnya. Metode ini penting untuk memetakan risiko penyakit yang berkelanjutan di seluruh permukaan area.

2.4.1 Inverse Distance Weighting (IDW)

Metode interpolasi deterministik IDW berasumsi bahwa titik yang lebih dekat memiliki bobot pengaruh yang lebih besar daripada titik yang lebih jauh. Nilai prediksi adalah rata-rata tertimbang dari semua titik sampel di sekitarnya, di mana bobot berbanding terbalik dengan jarak pangkat p.

2.4.2 Ordinary Kriging

Kriging adalah teknik geostatistik stokastik yang lebih rumit daripada IDW karena mempertimbangkan kedua jarak dan struktur korelasi spasial data yang dimodelkan melalui semivariogram. Dengan variansi minimum dan estimasi error variansi, normal Kriging menunjukkan tingkat ketidakpastian prediksi.

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 

Variabel penyakit DBD ini akan menjadi variabel dependen yang nantinya 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

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

3.3 Metode Analisis

Analisis dimulai dengan eksplorasi visual data spasial menggunakan pemetaan prevalensi DBD berbasis kuartil. Selanjutnya, metode K-Nearest Neighbor (KNN), yang dioptimasi untuk mendeteksi autokorelasi spasial global (Moran’s I) dan klaster lokal (LISA), digunakan untuk membuat matriks pembobot spasial. Tahap selanjutnya melibatkan estimasi model ekonometrika spasial (Spatial Autoregressive, Spatial Error Model, dan Spatial Durbin Model) untuk menguji pengaruh variabel sosio-ekonomi dan lingkungan, di mana pemilihan model terbaik didasarkan pada kriteria informasi (AIC dan BIC). Terakhir, pemetaan prediksi risiko secara kontinu dilakukan dengan membandingkan metode interpolasi spasial Inverse Distance Weighting (IDW) dan Ordinary Kriging. Metode ini divalidasi dengan menggunakan nilai Root Mean Square Error (RMSE) terkecil.

3.4 Alur Kerja

Alur kerja yang akan dilakukan meliputi: 1. Persiapan data 2. Uji pola spasial 3. Pemodelan Spasial Ekonometrik 4. Penentuan model terbaik 5. Pemodelan Interpolasi 6. Visualisasi dan Interpretasi

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/UZxqg9W9OZc?feature=shared

4.2 Analisis Deskriptif

Analisis deskriptif dan pemetaan dilakukan untuk melihat pola data.

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)
data <- data %>%
  mutate(
    Jumlah_Penduduk = Jumlah.Penduduk..Ribu. * 1000,
    DBD_prev         = (DBD / Jumlah_Penduduk) * 100000
  )

# Gabungkan shapefile dan data
data_map <- batas_prov %>%
  left_join(data, by = "Provinsi")
# Perhitungan nilai k optimal
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_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) {
  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 - DBD 
##        0%       25%       50%       75%      100% 
##   0.00000  53.89009  76.95545 123.68850 269.45797 
## ==============================

Peta sebaran kuartil menunjukkan variasi spasial yang signifikan dalam insidensi DBD di Indonesia. Provinsi dengan beban kasus tertinggi terlihat di Kalimantan, Sulawesi, dan Nusa Tenggara, sementara provinsi dengan beban kasus paling rendah terlihat di Sumatera dan Papua.

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

print(autocorr_results)
##                   Penyakit Moran_I Moran_p Geary_C Geary_p
## Moran I statistic      DBD  0.1313  0.0861  0.8439  0.1091

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) {
  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]
  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)
  }
  nb_full <- weights_list$neighbors
  nb_subset <- subset(nb_full, valid_idx)
  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)
  }
  stat <- as.numeric(morans$statistic["Moran I statistic"])
  pval <- as.numeric(morans$p.value)
  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)
  lisa <- tryCatch(
    localmoran(prevalence_valid, weights_subset, zero.policy = TRUE),
    error = function(e) return(NULL)
  )
  if (is.null(lisa)) return(NULL)
  lisa <- as.data.frame(lisa)
  p_col <- grep("Pr", names(lisa), value = TRUE)
  if (length(p_col) == 0) p_col <- names(lisa)[5]
  p_values <- lisa[[p_col[1]]]
  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)
}
weights_per_disease <- list(
  DBD = list(neighbors = nb_knn_dbd, weights = listw_knn_dbd)
)
diseases <- names(weights_per_disease)
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.08606237 -0.02702703 0.01344398

Hasil uji Moran’s I dengan p-value sebesar 0,086 (> 0,05) menunjukkan bahwa secara statistik tidak terdapat autokorelasi spasial global yang signifikan, yang mengindikasikan bahwa pola penyebaran kasus DBD antarprovinsi di Indonesia cenderung bersifat acak (random).

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

Terlihat dari peta LISA menunjukan bahwa terdapat outlier pada bagian bawah pulau Kalimantan.

4.4 Pembentukan Model

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(spdep)
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
# 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) sehingga semua variabel layak digunakan dalam pembentukan model regresi.")
  }
  
  return(vif_values)
}
  
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)
}
diseases <- c("DBD")  
x_vars <- c("Persentase.Rumah.Tangga.Memiliki.Akses.terhadap.Sanitasi.Layak..Persen.", "Kepadatan.Penduduk.per.km.persegi..Km2.", "Jumlah.Desa.yang.Memiliki.Rumah.Sakit", "Presentase.Penduduk.Miskin")
vif_values <- check_multicollinearity(data_map, x_vars)
## 📊 Hasil Uji Multikolinearitas (VIF):
## Persentase.Rumah.Tangga.Memiliki.Akses.terhadap.Sanitasi.Layak..Persen. 
##                                                                2.416490 
##                                 Kepadatan.Penduduk.per.km.persegi..Km2. 
##                                                                1.075299 
##                                   Jumlah.Desa.yang.Memiliki.Rumah.Sakit 
##                                                                1.098714 
##                                              Presentase.Penduduk.Miskin 
##                                                                2.520005
## ✅ Tidak ada multikolinearitas serius (semua VIF < 10) sehingga semua variabel layak digunakan dalam pembentukan model regresi.
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.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         DBD                  OLS 421.9835 431.8091
## lag         DBD    Spatial Lag (SAR) 423.7011 435.1642
## error       DBD  Spatial Error (SEM) 423.9748 435.4379
## durbin      DBD Spatial Durbin (SDM) 428.1433 446.1567

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

4.6 Pemodelan Interpolasi

4.6.1 Metode Inverse Distance Weighting (IDW)

library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
library(sf)
library(ggplot2)
points_dbd <- 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
bbox <- st_bbox(points_dbd)
x_seq <- seq(bbox["xmin"], bbox["xmax"], by = 0.1)
y_seq <- seq(bbox["ymin"], bbox["ymax"], by = 0.1)
grid_df <- expand.grid(x = x_seq, y = y_seq)
grid_rectangle <- st_as_sf(grid_df, coords = c("x", "y"), crs = st_crs(points_dbd))
cat("Sedang memotong grid sesuai bentuk peta Indonesia (mungkin agak lama)...\n")
## Sedang memotong grid sesuai bentuk peta Indonesia (mungkin agak lama)...
batas_indo_union <- st_union(batas_prov)
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
grid_indo_fixed <- st_intersection(grid_rectangle, batas_indo_union)
## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar
cat("Selesai memotong grid.\n")
## Selesai memotong grid.
idw_result_fixed <- idw(DBD_prev ~ 1, locations = points_dbd, newdata = grid_indo_fixed, idp = 2.0)
## [inverse distance weighted interpolation]
ggplot() +
  geom_sf(data = idw_result_fixed, aes(color = var1.pred), size = 0.1) +
  scale_color_viridis_c(option = "plasma", name = "Prediksi\nPrevalensi DBD") +
  geom_sf(data = batas_prov, fill = NA, color = "white", size = 0.3, alpha = 0.5) +
  labs(
    title = "Peta Interpolasi IDW ",
    subtitle = "Prediksi Prevalensi DBD "
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

Visualisasi interpolasi spasial menggunakan metode Inverse Distance Weighting (IDW) memprediksi bahwa risiko prevalensi DBD tertinggi terkonsentrasi secara signifikan di wilayah Kalimantan, sedangkan wilayah barat Indonesia seperti Sumatera cenderung menunjukkan tingkat prediksi kerawanan yang lebih rendah.

4.6.2 Metode Ordinary Kriging

# A. Membuat Variogram Empiris
v_emp <- variogram(DBD_prev ~ 1, data = points_dbd)
plot(v_emp, main = "Variogram Empiris DBD")

# B. Fitting Model Variogram
v_mod <- fit.variogram(v_emp, model = vgm(model = "Sph"))
# Cek hasil fitting 
print(v_mod)
##   model    psill    range
## 1   Sph 2926.376 376.4547
plot(v_emp, v_mod, main = "Fitting Variogram (Empiris vs Model Spherical)")

# C. Melakukan Ordinary Kriging (FIXED)
ok_result_fixed <- krige(
  formula = DBD_prev ~ 1,
  locations = points_dbd,
  newdata = grid_indo_fixed,
  model = v_mod
)
## [using ordinary kriging]
# D. VISUALISASI PREDIKSI KRIGING 
ggplot() +
  geom_sf(data = ok_result_fixed, aes(color = var1.pred), size = 0.1) +
  scale_color_viridis_c(option = "viridis", name = "Prediksi\nPrevalensi DBD") +
  geom_sf(data = batas_prov, fill = NA, color = "white", size = 0.3, alpha = 0.5) +
  labs(
    title = "Peta Interpolasi Ordinary Kriging",
    subtitle = "Prediksi Prevalensi DBD "
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

# E. VISUALISASI VARIAN (KETIDAKPASTIAN)
ggplot() +
  geom_sf(data = ok_result_fixed, aes(color = var1.var), size = 0.1) +
  scale_color_viridis_c(option = "magma", name = "Varian\n(Error)") +
  geom_sf(data = batas_prov, fill = NA, color = "white", size = 0.3, alpha = 0.5) +
  labs(
    title = "Peta Varian Ketidakpastian (Kriging)",
    subtitle = "Warna terang = Tingkat ketidakpastian tinggi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

  • Fitting Variogram: Plot ini menunjukkan bagaimana fungsi memodelkan struktur autokorelasi spasial data prevalensi DBD, yang mencapai titik jenuh pada jarak tertentu.

  • Peta Interpolasi: Peta ini memprediksi pola penyebaran DBD di Indonesia, di mana konsentrasi tertinggi terlihat memusat di area Kalimantan.

  • Peta Varian: Peta ini menunjukkan reliabilitas prediksi. Area berwarna terang menunjukkan tempat di mana data sampel sangat sedikit, yang mengakibatkan kemungkinan error yang tinggi pada hasil prediksi.

4.6.3 Validasi Model (Cross-Validation) & Pemilihan Model Terbaik

cv_idw <- krige.cv(DBD_prev ~ 1, points_dbd, nfold = nrow(points_dbd), set = list(idp = 2.0))
rmse_idw <- sqrt(mean(cv_idw$residual^2))
cv_ok <- krige.cv(DBD_prev ~ 1, points_dbd, model = v_mod, nfold = nrow(points_dbd))
rmse_ok <- sqrt(mean(cv_ok$residual^2))
cat("\n============================================\n")
## 
## ============================================
cat("      HASIL PERBANDINGAN ERROR (RMSE)       \n")
##       HASIL PERBANDINGAN ERROR (RMSE)
cat("============================================\n")
## ============================================
cat("RMSE IDW              :", round(rmse_idw, 4), "\n")
## RMSE IDW              : 57.131
cat("RMSE Ordinary Kriging :", round(rmse_ok, 4), "\n")
## RMSE Ordinary Kriging : 59.3079
cat("--------------------------------------------\n")
## --------------------------------------------
if(rmse_idw < rmse_ok) {
  cat("KESIMPULAN: Metode IDW lebih baik untuk data ini\n")
  cat("karena memiliki error (RMSE) yang lebih kecil.\n")
} else {
  cat("KESIMPULAN: Metode Ordinary Kriging lebih baik untuk data ini\n")
  cat("karena memiliki error (RMSE) yang lebih kecil.\n")
}
## KESIMPULAN: Metode IDW lebih baik untuk data ini
## karena memiliki error (RMSE) yang lebih kecil.
cat("============================================\n")
## ============================================

Metode Inverse Distance Weighting (IDW) ternyata lebih baik daripada Ordinary Kriging karena menghasilkan tingkat kesalahan prediksi (RMSE) yang lebih rendah.

4.6.4 Trend Surface Analysis

coords_pts <- st_coordinates(points_dbd)
points_dbd$X <- coords_pts[,1]
points_dbd$Y <- coords_pts[,2]
coords_grid <- st_coordinates(grid_indo_fixed)
grid_indo_fixed$X <- coords_grid[,1]
grid_indo_fixed$Y <- coords_grid[,2]
trend_result_fixed <- krige(
  formula = DBD_prev ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), 
  locations = points_dbd, 
  newdata = grid_indo_fixed, 
  model = NULL
)
## [ordinary or weighted least squares prediction]
ggplot() +
  geom_sf(data = trend_result_fixed, aes(color = var1.pred), size = 0.1) +
  scale_color_viridis_c(option = "cividis", name = "Prediksi\n(Trend Global)") +
  geom_sf(data = batas_prov, fill = NA, color = "white", size = 0.3, alpha = 0.5) +
  labs(
    title = "Peta Trend Surface Analysis (Orde 2)",
    subtitle = "Pola Global: Menangkap tren spasial secara umum"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

Pola tren spasial global berbentuk kuadratik. Nilai prediksi tertinggi ditemukan di wilayah tengah Indonesia (Kalimantan), yang kemudian turun secara bertahap ke arah Sumatera dan Papua di ujung barat dan timur.

5 Kesimpulan dan Saran

5.1 Kesimpulan

Analisis spasial terhadap kasus Demam Berdarah Dengue (DBD) di Indonesia tahun 2024 menunjukkan distribusi spasial yang signifikan, dengan hotspot, atau konsentrasi kasus tertinggi, terkonsentrasi di Kalimantan, Sulawesi, dan Nusa Tenggara. Meskipun pengujian autokorelasi global dengan Moran’s I menunjukkan pola penyebaran yang cenderung acak (random), analisis lokal dengan LISA berhasil menemukan klaster spasial dan outlier khusus di wilayah tertentu. Meskipun model Ordinary Least Squares (OLS) menghasilkan kriteria informasi (AIC) terendah, model Spatial Autoregressive (SAR) dinilai lebih tepat secara teoritis. Variabel yang berpengaruh pada fenomena ini adalah sanitasi, kepadatan penduduk, dan kemiskinan.

Metode Inverse Distance Weighting (IDW) menunjukkan hasil yang lebih akurat daripada Ordinary Kriging untuk pemetaan prediksi risiko di wilayah yang tidak dapat diakses data. Ini ditunjukkan oleh nilai error (RMSE) IDW yang lebih rendah, yaitu 57,13, dibandingkan dengan Kriging standar, yang sebesar 59,31. Hasil ini menunjukkan bahwa metode deterministik berbasis jarak lebih cocok untuk menangani sebaran data DBD dalam penelitian ini daripada metode geostatistik yang membutuhkan asumsi stasioneritas ketat. Secara visual, hasil Trend Surface Analysis (Orde 2) dan interpolasi memperlihatkan pola tren global berbentuk kuadratik, di mana konsentrasi kasus tertinggi memuncak di wilayah tengah Indonesia (Kalimantan) dan menurun ke arah Barat (Sumatera) maupun Timur (Papua). Hal ini mengindikasikan bahwa area tengah Indonesia perlu menjadi prioritas utama dalam mitigasi kesehatan.

5.2 Saran

Berdasarkan hasil analisis yang telah dilakukan, terdapat beberapa saran yang dapat diberikan:

  1. Fokus Penanganan di Wilayah Tengah Indonesia: Menurut hasil visualisasi Trend Surface Analysis dan peta LISA, wilayah tengah Indonesia memiliki konsentrasi kasus DBD yang paling tinggi. Oleh karena itu, sumber daya seperti pengasapan massal (fogging) atau penyebaran bubuk abate harus didistribusikan di wilayah ini oleh pemerintah. Agar kasus tidak menyebar ke wilayah tetangga, penanganan yang fokus harus dilakukan di wilayah “hotspot” ini.

  2. Perbaikan Sanitasi di Daerah Padat Penduduk: Hasil pemodelan Spatial Autoregressive (SAR) menunjukkan bahwa tingkat sanitasi yang buruk dan kepadatan penduduk berkontribusi secara signifikan terhadap peningkatan kasus. Dengan demikian, upaya pencegahan tidak dapat hanya bergantung pada sektor kesehatan. Untuk meningkatkan sanitasi lingkungan, terutama di daerah permukiman padat, diperlukan kerja sama lintas sektor. Pengobatan kuratif tidak akan bekerja sama dengan mengurangi tempat perkembangbiakan nyamuk (breeding places) di lingkungan kumuh.

  3. Pemetaan Risiko dengan Metode Interpolasi Inverse Distance Weighting (IDW): Dalam penelitian ini, metode IDW terbukti memiliki tingkat error (RMSE) yang lebih rendah dibandingkan Ordinary Kriging. Peta prediksi yang dihasilkan dari metode ini bisa digunakan dinas terkait sebagai acuan kewaspadaan dini (early warning), terutama di area yang diprediksi berisiko tinggi (zona merah) meskipun belum ada laporan kasus yang masuk.

  4. Penambahan Variabel Iklim untuk Penelitian Selanjutnya: Penelitian ini memiliki keterbatasan karena hanya menggunakan data satu periode waktu (tahun 2024). Peneliti disarankan untuk menggunakan metode spatiotemporal dengan data deret waktu karena DBD sangat dipengaruhi oleh musim. Selain itu, sangat disarankan untuk menambah variabel faktor iklim seperti curah hujan dan kelembapan udara agar model yang dihasilkan dapat menggambarkan pola penyebaran penyakit dengan lebih akurat.

6 Referensi

Achmadi, U. F. (2009). Manajemen penyakit berbasis wilayah. Kesmas, 3(4), 147–153. https://doi.org/10.21109/kesmas.v3i3.228

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

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

LeSage, J. and Pace, R.K. (2009) Introduction to Spatial Econometrics. Chapman and Hall/CRC, New York. https://doi.org/10.1201/9781420064254

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

Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(Suppl.), 234–240. https://doi.org/10.2307/143141

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

WHO. (2023). Dengue and severe dengue. World Health Organization.

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