Analisis Spasial Kasus Leptospirosis di Jawa Tengah

UAS Mata Kuliah Analisis Statistik Spasial

ABSTRAK

Penyakit Leptospirosis merupakan zoonosis yang menjadi permasalahan kesehatan masyarakat persisten di Provinsi Jawa Tengah, terutama dipengaruhi oleh interaksi faktor demografi dan akses fasilitas kesehatan. Penelitian ini secara strategis mendukung pencapaian Sustainable Development Goals (SDGs), khususnya Tujuan 3 mengenai kehidupan sehat dan sejahtera melalui pengendalian penyakit menular, serta Tujuan 6 terkait akses air bersih dan sanitasi layak yang menjadi faktor risiko utama penularan bakteri Leptospira. Penelitian ini bertujuan untuk menganalisis pola persebaran spasial kasus Leptospirosis di 35 kabupaten/kota di Provinsi Jawa Tengah pada periode 2022-2025 serta menentukan model estimasi terbaik dengan mempertimbangkan aspek dependensi spasial. Metodologi yang digunakan mencakup uji autokorelasi spasial melalui statistik Geary’s C dan perbandingan model ekonometrika spasial yang meliputi Ordinary Least Square (OLS), Spatial Autoregressive (SAR), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM).

Hasil analisis menunjukkan bahwa distribusi kasus Leptospirosis cenderung bersifat acak secara lokal dengan nilai Geary’s C sebesar 0,9194 (\(p\text{-value} = 0,2409\)), yang mengindikasikan tidak adanya pola klasterisasi yang signifikan secara statistik pada taraf kepercayaan 95%. Berdasarkan kriteria Akaike Information Criterion (AIC), model OLS terpilih sebagai model paling efisien dengan nilai AIC terendah sebesar 305,9 dibandingkan model spasial lainnya. Uji diagnostik sisaan menunjukkan nilai Moran’s I Residual sebesar 0,0587 (\(p\text{-value} = 0,1925\)), membuktikan bahwa model OLS telah berhasil mengekstrak seluruh informasi tanpa menyisakan autokorelasi spasial pada residual. Temuan ini diharapkan dapat menjadi rujukan dalam perancangan intervensi kesehatan dan perbaikan infrastruktur sanitasi yang lebih tepat sasaran guna mengakselerasi target SDGs di Jawa Tengah.

Kata Kunci: Leptospirosis, Analisis Statistik Spasial, OLS, SDGs, Jawa Tengah

---

1 PENDAHULUAN

1.1 Latar Belakang

Leptospirosis merupakan salah satu penyakit zoonosis yang masih menjadi permasalahan kesehatan masyarakat di Indonesia, khususnya di wilayah dengan kondisi lingkungan yang mendukung perkembangan bakteri Leptospira, seperti daerah dengan curah hujan tinggi, sanitasi yang kurang memadai, serta kepadatan penduduk yang relatif tinggi. Provinsi Jawa Tengah sebagai salah satu wilayah dengan karakteristik geografis dan demografis yang beragam menunjukkan variasi jumlah kasus leptospirosis antar kabupaten/kota, sehingga memerlukan analisis yang mampu menangkap perbedaan tersebut secara komprehensif.

Perbedaan kondisi sosial, lingkungan, dan infrastruktur kesehatan antar wilayah menyebabkan risiko leptospirosis tidak bersifat homogen. Kabupaten/kota yang memiliki masalah sanitasi, pengelolaan limbah, atau akses layanan kesehatan yang terbatas cenderung memiliki kerentanan lebih tinggi terhadap penyebaran penyakit ini. Oleh karena itu, analisis berbasis wilayah menjadi penting agar kebijakan penanggulangan dapat disusun secara lebih tepat sasaran sesuai karakteristik masing-masing daerah.

Selain faktor lokal, aspek keterkaitan antar wilayah juga perlu diperhatikan. Secara geografis, kabupaten/kota di Jawa Tengah saling berdekatan dan memungkinkan adanya pola spasial dalam distribusi kasus leptospirosis. Hal ini mendorong penggunaan pendekatan statistik spasial, tidak hanya untuk mengidentifikasi hubungan antara kasus leptospirosis dengan faktor-faktor penjelas, tetapi juga untuk menguji apakah terdapat pengaruh ketergantungan spasial antar wilayah. Dengan demikian, model regresi spasial digunakan sebagai pengembangan dari regresi linier klasik agar potensi autokorelasi spasial dapat diakomodasi.

Kajian ini menjadi penting dalam konteks pencapaian Sustainable Development Goals (SDGs), khususnya Tujuan 3: Good Health and Well-being, yang menekankan upaya menjamin kehidupan yang sehat dan meningkatkan kesejahteraan bagi seluruh masyarakat. Pengendalian penyakit menular seperti leptospirosis merupakan bagian dari target SDGs untuk menurunkan angka kesakitan dan kematian akibat penyakit. Selain itu, penelitian ini juga berkaitan dengan SDGs Tujuan 6: Clean Water and Sanitation, karena leptospirosis sangat erat hubungannya dengan kualitas sanitasi lingkungan dan akses terhadap air bersih. Melalui pemahaman pola spasial dan faktor-faktor yang memengaruhi kasus leptospirosis, hasil penelitian ini diharapkan dapat menjadi dasar bagi pemerintah daerah dalam merancang intervensi kesehatan dan perbaikan sanitasi yang lebih terarah.

1.2 Identifikasi Masalah

Berdasarkan latar belakang di atas, masalah penelitian dapat dirumuskan sebagai berikut:

1. Bagaimana pola distribusi spasial kasus leptospirosis di Jawa Tengah?
2. Apakah terdapat autokorelasi spasial dalam distribusi kasus leptospirosis?
3. Variabel independen apa saja yang signifikan mempengaruhi kasus leptospirosis dengan mempertimbangkan efek spasial?
4. Model spatial mana yang paling sesuai untuk kasus ini?

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk:

1. Mengidentifikasi pola distribusi spasial [variabel dependen]
2. Mengukur tingkat autokorelasi spasial menggunakan Moran’s I, Geary’s C, dan LISA
3. Menganalisis faktor-faktor yang mempengaruhi kasus leptospirosis menggunakan model spatial econometrics
4. Menentukan model terbaik berdasarkan kriteria statistik dan interpretasi substantif

1.4 Batasan Penelitian

1. Unit analisis: Kabupaten/Kota di Jawa Tengah
2. Periode analisis: Tahun 2022–2025
3. Variabel yang digunakan: Jumlah Penduduk per Kabupaten/Kota, Jumlah Puskesmas, dan Jumlah Kasus Leptospirosis di Jawa Tengah
4. Metode spasial: SAR, SEM, SDM

---

2 TINJAUAN PUSTAKA

2.1 Spatial Dependence dan Hukum Tobler

Hukum pertama geografi Tobler menyatakan bahwa “segala sesuatu berhubungan dengan yang lain, tetapi hal-hal yang dekat lebih berhubungan daripada hal-hal yang jauh” (Tobler, 1970). Prinsip ini menjadi dasar spatial dependence, yaitu nilai variabel di suatu lokasi dipengaruhi oleh nilai variabel di lokasi tetangganya.

2.2 Matriks Pembobot Spasial

Matriks pembobot spasial (W) mendefinisikan struktur ketetanggaan antar wilayah. Penelitian ini menggunakan:

• Queen Contiguity: Dua wilayah dianggap tetangga jika berbagi sisi atau titik sudut
• K-Nearest Neighbors: Setiap wilayah memiliki k tetangga terdekat

2.3 Autokorelasi Spasial

2.3.1 Moran’s I

Moran’s I mengukur korelasi spasial global dengan rumus:

\[I = \frac{n}{\sum_i \sum_j w_{ij}} \frac{\sum_i \sum_j w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_i (x_i - \bar{x})^2}\]

Nilai I > 0 menunjukkan autokorelasi positif, I < 0 autokorelasi negatif, dan I ≈ 0 tidak ada autokorelasi.

2.3.2 Geary’s C

Geary’s C fokus pada perbedaan nilai antar tetangga:

\[C = \frac{(n-1)}{2\sum_i \sum_j w_{ij}} \frac{\sum_i \sum_j w_{ij}(x_i - x_j)^2}{\sum_i (x_i - \bar{x})^2}\]

Nilai C < 1 menunjukkan autokorelasi positif, C > 1 autokorelasi negatif.

2.3.3 LISA (Local Indicators of Spatial Association)

LISA mengidentifikasi hotspot dan coldspot lokal:

\[I_i = \frac{(x_i - \bar{x})}{\sigma^2} \sum_j w_{ij}(x_j - \bar{x})\]

Klasifikasi pola: High-High, Low-Low, High-Low, Low-High.

2.4 Model Spatial Econometrics

2.4.1 Spatial Autoregressive Model (SAR)

Model SAR mengasumsikan spatial dependence pada variabel dependen:

\[y = \rho Wy + X\beta + \epsilon\]

2.4.2 Spatial Error Model (SEM)

Model SEM mengasumsikan spatial dependence pada error term:

\[y = X\beta + u, \quad u = \lambda Wu + \epsilon\]

2.4.3 Spatial Durbin Model (SDM)

Model SDM memasukkan spatial lag pada variabel independen dan dependen:

\[y = \rho Wy + X\beta + WX\theta + \epsilon\]

2.5 Pemilihan Model

Kriteria pemilihan model meliputi:

• AIC (Akaike Information Criterion): Lebih kecil lebih baik
• BIC (Bayesian Information Criterion): Lebih kecil lebih baik
• LM Test: Lagrange Multiplier test untuk SAR dan SEM
• Likelihood Ratio Test: Perbandingan goodness-of-fit

---

3 METODOLOGI PENELITIAN

3.1 Data dan Sumber Data

Data yang digunakan dalam penelitian ini meliputi:

• Unit Analisis: Kabupaten/Kota di Provinsi Jawa Tengah (35 Kab./Kota)
• Tahun: Tahun 2022–2025
• Sumber Data: Portal Data Jawa Tengah dan BPS Jawa Tengah

Variabel Penelitian:

1. Variabel Dependen:
   • Kasus Leptospirosis (Y): Rata-rata jumlah kasus leptospirosis per tahun periode 2022-2025 di setiap kabupaten/kota. Leptospirosis adalah penyakit zoonosis yang disebabkan oleh bakteri Leptospira yang ditularkan melalui air atau tanah yang terkontaminasi urine hewan terinfeksi.
2. Variabel Independen:
   • Jumlah Penduduk (X1): Rata-rata jumlah penduduk dalam ribuan per tahun periode 2022-2025. Populasi yang lebih besar dapat meningkatkan risiko paparan terhadap sumber infeksi leptospirosis.
   • Jumlah Puskesmas (X2): Rata-rata jumlah puskesmas per tahun periode 2022-2025. Ketersediaan fasilitas kesehatan dasar berpengaruh terhadap deteksi, pelaporan, dan penanganan kasus leptospirosis.

3.2 Software dan Packages

# Instalasi packages jika belum terinstal
packages <- c(
  "readxl", "sf", "spdep", "spatialreg", "tmap", "ggplot2",
  "dplyr", "tidyr", "knitr", "kableExtra", "viridis",
  "RColorBrewer", "leaflet", "plotly", "corrplot",
  "scales", "gridExtra", "GWmodel", "stringr", "gstat"
)

# Install jika belum ada
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Load semua packages
invisible(lapply(packages, library, character.only = TRUE))

3.3 Tahapan Analisis

Analisis dilakukan melalui tahapan berikut:

1. Eksplorasi Data Spasial
   • Statistik deskriptif
   • Visualisasi peta distribusi
2. Analisis Autokorelasi Spasial
   • Konstruksi matriks pembobot spasial
   • Perhitungan Moran’s I global
   • Perhitungan Geary’s C
   • Analisis LISA lokal
3. Pemodelan Spatial Econometrics
   • Estimasi model OLS (baseline)
   • Uji spatial dependence (LM test)
   • Estimasi SAR, SEM, SDM
   • Perbandingan model
4. Interpolasi Spasial
   • Kriging atau IDW
   • Validasi cross-validation
5. Interpretasi dan Visualisasi

---

4 HASIL DAN PEMBAHASAN

4.1 Load dan Prepare Data

# Load data Excel
data <- read_excel("D:/Semester 5/Spasial/DATA UAS/Dashboard/Data UAS Statistik Spasial.xlsx")

# Load GeoJSON
indo_sf <- st_read("D:/Semester 5/Spasial/DATA UAS/Dashboard/gadm41_IDN_2.json", 
                   quiet = TRUE)

# Cek struktur data
cat("Struktur Data:\n")

## Struktur Data:

str(data)

## tibble [140 × 5] (S3: tbl_df/tbl/data.frame)
##  $ Kabupaten/Kota            : chr [1:140] "Kab.Cilacap" "Kab.Cilacap" "Kab.Cilacap" "Kab.Cilacap" ...
##  $ Tahun                     : num [1:140] 2022 2023 2024 2025 2022 ...
##  $ Jumlah Penduduk (Ribu)    : num [1:140] 1989 2008 2027 2046 1806 ...
##  $ Jumlah Puskesmas          : num [1:140] 38 38 38 38 40 40 40 40 22 22 ...
##  $ Jumlah Kasus Leptospirosis: num [1:140] 12 37 6 20 34 88 24 108 0 2 ...

cat("\n\nKolom GeoJSON:\n")

## 
## 
## Kolom GeoJSON:

names(indo_sf)

##  [1] "GID_2"     "GID_0"     "COUNTRY"   "GID_1"     "NAME_1"    "NL_NAME_1"
##  [7] "NAME_2"    "VARNAME_2" "NL_NAME_2" "TYPE_2"    "ENGTYPE_2" "CC_2"     
## [13] "HASC_2"    "geometry"

cat("\n\nDimensi Data:\n")

## 
## 
## Dimensi Data:

cat("Data Excel:", nrow(data), "baris,", ncol(data), "kolom\n")

## Data Excel: 140 baris, 5 kolom

cat("Data Spasial:", nrow(indo_sf), "wilayah\n")

## Data Spasial: 502 wilayah

# Preview data
cat("\n\nPreview Data Excel:\n")

## 
## 
## Preview Data Excel:

head(data)

Interpretasi: Data penelitian terdiri dari 140 observasi (35 kabupaten/kota × 4 tahun) dengan 3 variabel. Data spasial berbentuk poligon yang merepresentasikan batas administratif kabupaten/kota di Jawa Tengah.

# Standarisasi nama kolom
data <- data %>%
  rename(
    kabkota = `Kabupaten/Kota`,
    tahun = Tahun,
    jumlah_penduduk = `Jumlah Penduduk (Ribu)`,
    jumlah_puskesmas = `Jumlah Puskesmas`,
    kasus_leptospirosis = `Jumlah Kasus Leptospirosis`
  )

# Agregasi data per kabupaten/kota (rata-rata 2022-2025)
data_avg <- data %>%
  group_by(kabkota) %>%
  summarise(
    y_var = mean(kasus_leptospirosis, na.rm = TRUE),  # Rata-rata kasus
    x1 = mean(jumlah_penduduk, na.rm = TRUE),         # Rata-rata penduduk
    x2 = mean(jumlah_puskesmas, na.rm = TRUE),        # Rata-rata puskesmas
    total_kasus = sum(kasus_leptospirosis, na.rm = TRUE),
    .groups = 'drop'
  )

cat("Wilayah di Data Excel:", nrow(data_avg), "\n")

## Wilayah di Data Excel: 35

cat("Wilayah di GeoJSON (Seluruh Indonesia):", nrow(indo_sf), "\n\n")

## Wilayah di GeoJSON (Seluruh Indonesia): 502

# FILTER HANYA JAWA TENGAH
names(indo_sf)

##  [1] "GID_2"     "GID_0"     "COUNTRY"   "GID_1"     "NAME_1"    "NL_NAME_1"
##  [7] "NAME_2"    "VARNAME_2" "NL_NAME_2" "TYPE_2"    "ENGTYPE_2" "CC_2"     
## [13] "HASC_2"    "geometry"

unique(indo_sf$NAME_1)

##  [1] "Aceh"              "Bali"              "BangkaBelitung"   
##  [4] "Banten"            "Bengkulu"          "Gorontalo"        
##  [7] "JakartaRaya"       "Jambi"             "JawaBarat"        
## [10] "JawaTengah"        "JawaTimur"         "KalimantanBarat"  
## [13] "KalimantanSelatan" "KalimantanTengah"  "KalimantanTimur"  
## [16] "KalimantanUtara"   "KepulauanRiau"     "Lampung"          
## [19] "Maluku"            "MalukuUtara"       "NusaTenggaraBarat"
## [22] "NusaTenggaraTimur" "Papua"             "PapuaBarat"       
## [25] "Riau"              "SulawesiBarat"     "SulawesiSelatan"  
## [28] "SulawesiTengah"    "SulawesiTenggara"  "SulawesiUtara"    
## [31] "SumateraBarat"     "SumateraSelatan"   "SumateraUtara"    
## [34] "Yogyakarta"

indo_sf_jateng <- indo_sf %>%
  filter(NAME_1 == "JawaTengah")

cat("Wilayah di GeoJSON (Jawa Tengah saja):", nrow(indo_sf_jateng), "\n\n")

## Wilayah di GeoJSON (Jawa Tengah saja): 36

# Lihat sample nama di Jawa Tengah
cat("Sample nama wilayah di GeoJSON Jawa Tengah (NAME_2):\n")

## Sample nama wilayah di GeoJSON Jawa Tengah (NAME_2):

print(head(indo_sf_jateng$NAME_2, 10))

##  [1] "Banjarnegara" "Banyumas"     "Batang"       "Blora"        "Boyolali"    
##  [6] "Brebes"       "Cilacap"      "Demak"        "Grobogan"     "Jepara"

cat("\n\nSample nama wilayah di Excel:\n")

## 
## 
## Sample nama wilayah di Excel:

print(head(data_avg$kabkota, 10))

##  [1] "Kab.Banjarnegara" "Kab.Banyumas"     "Kab.Batang"       "Kab.Blora"       
##  [5] "Kab.Boyolali"     "Kab.Brebes"       "Kab.Cilacap"      "Kab.Demak"       
##  [9] "Kab.Grobogan"     "Kab.Jepara"

library(stringr)

data_avg <- data_avg %>%
  mutate(
    nama_clean = str_to_lower(
      str_trim(
        str_replace(
          str_replace(kabkota, "^Kab\\.", ""),
          "^Kota ", ""
        )
      )
    )
  )

cat("\n\nSetelah cleaning:\n")

## 
## 
## Setelah cleaning:

cat("Excel:\n")

## Excel:

print(head(data_avg %>% select(kabkota, nama_clean), 10))

## # A tibble: 10 × 2
##    kabkota          nama_clean  
##    <chr>            <chr>       
##  1 Kab.Banjarnegara banjarnegara
##  2 Kab.Banyumas     banyumas    
##  3 Kab.Batang       batang      
##  4 Kab.Blora        blora       
##  5 Kab.Boyolali     boyolali    
##  6 Kab.Brebes       brebes      
##  7 Kab.Cilacap      cilacap     
##  8 Kab.Demak        demak       
##  9 Kab.Grobogan     grobogan    
## 10 Kab.Jepara       jepara

indo_sf_jateng <- indo_sf_jateng %>%
  mutate(
    nama_clean = str_to_lower(str_trim(NAME_2))
  )

cat("\n\nGeoJSON:\n")

## 
## 
## GeoJSON:

print(head(indo_sf_jateng %>% st_drop_geometry() %>% select(NAME_2, nama_clean), 10))

##          NAME_2   nama_clean
## 1  Banjarnegara banjarnegara
## 2      Banyumas     banyumas
## 3        Batang       batang
## 4         Blora        blora
## 5      Boyolali     boyolali
## 6        Brebes       brebes
## 7       Cilacap      cilacap
## 8         Demak        demak
## 9      Grobogan     grobogan
## 10       Jepara       jepara

head(data_avg$nama_clean, 10)

##  [1] "banjarnegara" "banyumas"     "batang"       "blora"        "boyolali"    
##  [6] "brebes"       "cilacap"      "demak"        "grobogan"     "jepara"

head(indo_sf_jateng$nama_clean, 10)

##  [1] "banjarnegara" "banyumas"     "batang"       "blora"        "boyolali"    
##  [6] "brebes"       "cilacap"      "demak"        "grobogan"     "jepara"

data_sf <- indo_sf_jateng %>%
  left_join(data_avg, by = "nama_clean")

sum(!is.na(data_sf$y_var))

## [1] 35

sum(is.na(data_sf$y_var))

## [1] 5

# Merge berdasarkan nama yang sudah di-clean
data_sf <- indo_sf_jateng %>%
  left_join(data_avg, by = "nama_clean")

# Cek hasil matching
matched <- sum(!is.na(data_sf$y_var))
not_matched <- sum(is.na(data_sf$y_var))

cat("\n\n=== HASIL MATCHING ===\n")

## 
## 
## === HASIL MATCHING ===

cat("✓ Wilayah berhasil match:", matched, "\n")

## ✓ Wilayah berhasil match: 35

cat("✗ Wilayah tidak match   :", not_matched, "\n")

## ✗ Wilayah tidak match   : 5

if(not_matched > 0) {
  cat("\n\nWilayah di GeoJSON yang TIDAK MATCH:\n")
  not_matched_geo <- data_sf %>% 
    filter(is.na(y_var)) %>% 
    st_drop_geometry() %>% 
    select(NAME_2, nama_clean)
  print(not_matched_geo)
  
  cat("\n\nWilayah di Excel yang TIDAK MATCH:\n")
  matched_clean <- data_sf %>% 
    filter(!is.na(y_var)) %>% 
    pull(nama_clean)
  not_matched_excel <- data_avg %>%
    filter(!nama_clean %in% matched_clean) %>%
    select(kabkota, nama_clean)
  print(not_matched_excel)
  
  cat("\n\n💡 Coba cek perbedaan penulisan antara kedua tabel di atas\n")
}

## 
## 
## Wilayah di GeoJSON yang TIDAK MATCH:
##            NAME_2      nama_clean
## 1    KotaMagelang    kotamagelang
## 2  KotaPekalongan  kotapekalongan
## 3    KotaSemarang    kotasemarang
## 4       KotaTegal       kotategal
## 5 WadukKedungombo wadukkedungombo
## 
## 
## Wilayah di Excel yang TIDAK MATCH:
## # A tibble: 0 × 2
## # ℹ 2 variables: kabkota <chr>, nama_clean <chr>
## 
## 
## 💡 Coba cek perbedaan penulisan antara kedua tabel di atas

# Filter hanya yang match
data_sf <- data_sf %>% filter(!is.na(y_var))

cat("\n\n=== DATA FINAL ===\n")

## 
## 
## === DATA FINAL ===

cat("Jumlah wilayah untuk analisis:", nrow(data_sf), "kabupaten/kota\n")

## Jumlah wilayah untuk analisis: 35 kabupaten/kota

cat("Total kasus 2022-2025        :", sum(data_sf$total_kasus), "kasus\n")

## Total kasus 2022-2025        : 2034 kasus

cat("Rata-rata kasus per wilayah  :", round(mean(data_sf$y_var), 2), "kasus/tahun\n")

## Rata-rata kasus per wilayah  : 14.53 kasus/tahun

cat("Rata-rata penduduk           :", round(mean(data_sf$x1), 1), "ribu jiwa\n")

## Rata-rata penduduk           : 1076.4 ribu jiwa

cat("Rata-rata puskesmas          :", round(mean(data_sf$x2), 1), "unit\n")

## Rata-rata puskesmas          : 25.2 unit

# Cek apakah semua 35 wilayah sudah match
if(nrow(data_sf) == 35) {
  cat("\n✅ SEMPURNA! Semua 35 kabupaten/kota Jawa Tengah berhasil di-match!\n")
} else {
  cat("\n⚠️ WARNING: Hanya", nrow(data_sf), "dari 35 wilayah yang match.\n")
  cat("   Mohon periksa wilayah yang tidak match di atas.\n")
}

## 
## ✅ SEMPURNA! Semua 35 kabupaten/kota Jawa Tengah berhasil di-match!

4.2 Statistik Deskriptif

# Variabel yang akan dianalisis
vars_analisis <- c("y_var", "x1", "x2")

# Statistik deskriptif
desc_stats <- data_sf %>%
  st_drop_geometry() %>%
  select(all_of(vars_analisis)) %>%
  summary()

print(desc_stats)

##      y_var              x1               x2       
##  Min.   : 0.000   Min.   : 122.2   Min.   : 5.00  
##  1st Qu.: 0.375   1st Qu.: 854.3   1st Qu.:21.00  
##  Median : 4.000   Median :1052.8   Median :26.00  
##  Mean   :14.529   Mean   :1076.4   Mean   :25.17  
##  3rd Qu.:25.000   3rd Qu.:1348.6   3rd Qu.:30.00  
##  Max.   :63.500   Max.   :2051.6   Max.   :40.00

# Tabel statistik lebih detail
desc_table <- data_sf %>%
  st_drop_geometry() %>%
  select(all_of(vars_analisis)) %>%
  summarise(
    across(everything(), 
           list(
             Mean = ~mean(.x, na.rm = TRUE),
             Median = ~median(.x, na.rm = TRUE),
             SD = ~sd(.x, na.rm = TRUE),
             Min = ~min(.x, na.rm = TRUE),
             Max = ~max(.x, na.rm = TRUE),
             CV = ~(sd(.x, na.rm = TRUE)/mean(.x, na.rm = TRUE))*100
           ))
  ) %>%
  pivot_longer(everything(), 
               names_to = c("Variable", "Statistic"),
               names_sep = "_") %>%
  pivot_wider(names_from = Statistic, values_from = value)

desc_table %>%
  kable(caption = "Tabel 1. Statistik Deskriptif Variabel Penelitian",
        digits = 2,
        format.args = list(big.mark = ",")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Tabel 1. Statistik Deskriptif Variabel Penelitian

Variable	var	Mean	Median	SD	Min	Max	CV
y	14.52857, 4.00000, 18.46596, 0.00000, 63.50000, 127.10103	NULL	NULL	NULL	NULL	NULL	NULL
x1	NULL	1076.432	1052.775	472.4567	122.25	2051.625	43.89099
x2	NULL	25.17143	26	9.040101	5	40	35.91414

Interpretasi:

• Rata-rata kasus leptospirosis di Jawa Tengah adalah 14.53 kasus per kabupaten/kota per tahun dengan standar deviasi 18.47. Coefficient of Variation (CV) sebesar 127.1% menunjukkan variabilitas yang sangat tinggi antar wilayah.

• Rentang kasus dari 0 hingga 63.5 mengindikasikan disparitas spasial yang signifikan dalam distribusi kasus leptospirosis.

• Jumlah penduduk berkisar antara 122 ribu hingga 2052 ribu jiwa, menunjukkan heterogenitas ukuran wilayah.

• Jumlah puskesmas bervariasi dari 5 hingga 40 unit, mengindikasikan perbedaan kapasitas layanan kesehatan dasar antar wilayah.

# Matriks korelasi
cor_matrix <- data_sf %>%
  st_drop_geometry() %>%
  select(all_of(vars_analisis)) %>%
  cor(use = "complete.obs")

# Visualisasi
corrplot(cor_matrix, 
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45,
         diag = FALSE,
         number.cex = 0.8,
         col = colorRampPalette(c("#BB4444", "#FFFFFF", "#4444BB"))(200))

Gambar 1. Matriks Korelasi Antar Variabel

Interpretasi:

Matriks korelasi menunjukkan bahwa variabel X1 memiliki korelasi positif lemah–sedang (r = 0,35) terhadap kasus leptospirosis, sedangkan X2 memiliki korelasi positif lemah (r = 0,18). Korelasi antar variabel independen sebesar 0,42 menunjukkan tidak adanya multikolinearitas yang kuat, sehingga kedua variabel layak digunakan dalam pemodelan.

4.3 Visualisasi Peta Distribusi Spasial

# Peta distribusi kasus leptospirosis
tm_shape(data_sf) +
  tm_polygons("y_var",
              palette = "YlOrRd",
              style = "jenks",
              n = 5,
              title = "Rata-rata Kasus\nLeptospirosis\n(2022-2025)",
              border.col = "white",
              border.alpha = 0.5,
              legend.hist = TRUE) +
  tm_layout(main.title = "Distribusi Spasial Kasus Leptospirosis di Jawa Tengah",
            main.title.size = 1.2,
            legend.outside = TRUE,
            legend.outside.position = "right",
            frame = FALSE) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_compass(position = c("right", "top")) +
  tm_credits("Sumber: Data UAS Statistik Spasial 2025", 
             position = c("left", "bottom"))

Gambar 2. Peta Distribusi Kasus Leptospirosis di Jawa Tengah

Interpretasi:

Peta choropleth rata-rata kasus leptospirosis tahun 2022–2025 menunjukkan adanya variasi spasial yang cukup jelas antar kabupaten/kota di Provinsi Jawa Tengah. Beberapa wilayah tampak memiliki rata-rata kasus yang tinggi hingga sangat tinggi, sementara sebagian besar wilayah lainnya berada pada kategori rendah. Pola ini mengindikasikan bahwa beban leptospirosis cenderung terlokalisasi pada wilayah tertentu, sehingga diperlukan perhatian dan intervensi kesehatan yang lebih terfokus pada daerah-daerah dengan intensitas kasus tinggi.

# Peta untuk semua variabel independen
map1 <- tm_shape(data_sf) +
  tm_polygons("x1",
              palette = "Blues",
              style = "jenks",
              n = 5,
              title = "Penduduk\n(Ribu)",
              border.col = "white",
              border.alpha = 0.3) +
  tm_layout(main.title = "Jumlah Penduduk",
            legend.outside = FALSE,
            frame = FALSE)

map2 <- tm_shape(data_sf) +
  tm_polygons("x2",
              palette = "Greens",
              style = "jenks",
              n = 5,
              title = "Jumlah\nPuskesmas",
              border.col = "white",
              border.alpha = 0.3) +
  tm_layout(main.title = "Jumlah Puskesmas",
            legend.outside = FALSE,
            frame = FALSE)

# Arrange maps
tmap_arrange(map1, map2, ncol = 2)

Gambar 3. Peta Distribusi Variabel Independen

Interpretasi:

Pola Distribusi Jumlah Penduduk (X1)

Peta menunjukkan adanya ketimpangan distribusi penduduk antarwilayah di Jawa Tengah, di mana beberapa kabupaten/kota menjadi pusat konsentrasi populasi, sementara wilayah lainnya berpenduduk relatif lebih jarang. Perbedaan ini berpotensi memengaruhi tingkat risiko dan jumlah absolut kasus leptospirosis, karena wilayah berpenduduk besar cenderung memiliki peluang paparan dan pelaporan kasus yang lebih tinggi.

Pola Distribusi Jumlah Puskesmas (X2):

Peta jumlah puskesmas menunjukkan variasi kapasitas layanan kesehatan primer antarwilayah. Daerah dengan jumlah puskesmas lebih banyak diharapkan memiliki kemampuan deteksi dan penanganan leptospirosis yang lebih baik, sedangkan wilayah dengan fasilitas terbatas berpotensi menghadapi kendala dalam akses layanan kesehatan.

4.4 Konstruksi Matriks Pembobot Spasial

# Cek validitas geometri
data_sf <- st_make_valid(data_sf)

# Konversi ke sp untuk spdep (jika perlu)
# data_sp <- as_Spatial(data_sf)

# Queen contiguity
nb_queen <- poly2nb(data_sf, queen = TRUE)

cat("Struktur Ketetanggaan (Queen Contiguity):\n")

## Struktur Ketetanggaan (Queen Contiguity):

cat("Jumlah wilayah:", length(nb_queen), "\n")

## Jumlah wilayah: 35

cat("Rata-rata tetangga:", mean(card(nb_queen)), "\n")

## Rata-rata tetangga: 5.142857

cat("Min tetangga:", min(card(nb_queen)), "\n")

## Min tetangga: 2

cat("Max tetangga:", max(card(nb_queen)), "\n")

## Max tetangga: 10

# Identifikasi wilayah tanpa tetangga (islands)
islands <- which(card(nb_queen) == 0)
if(length(islands) > 0) {
  cat("\nWilayah tanpa tetangga (islands):", length(islands), "\n")
  cat("Nama wilayah:", data_sf$NAME_2[islands], "\n")
}

# Atasi islands dengan k-nearest neighbors
if(length(islands) > 0) {
  coords <- st_coordinates(st_centroid(data_sf))
  k <- ceiling(mean(card(nb_queen))) + 1
  nb_knn <- knn2nb(knearneigh(coords, k = k))
  
  # Gabungkan queen dan knn untuk islands
  nb_combined <- nb_queen
  nb_combined[islands] <- nb_knn[islands]
  
  cat("\nSetelah penambahan KNN untuk islands:\n")
  cat("Rata-rata tetangga:", mean(card(nb_combined)), "\n")
  
  nb_final <- nb_combined
} else {
  nb_final <- nb_queen
}

# Row-standardized weights
w_list <- nb2listw(nb_final, style = "W", zero.policy = TRUE)

# Visualisasi struktur ketetanggaan
plot(st_geometry(data_sf), border = "lightgray", main = "Struktur Ketetanggaan")
plot(nb_final, coords = st_coordinates(st_centroid(data_sf)), 
     add = TRUE, col = "red", lwd = 0.5)

Interpretasi: Nilai suatu wilayah tidak berdiri sendiri, tetapi berpotensi dipengaruhi oleh wilayah-wilayah yang terhubung dengannya.

4.5 Analisis Autokorelasi Spasial Global

4.5.1 Moran’s I

# Hitung Moran's I
moran_test <- moran.test(data_sf$y_var, w_list, zero.policy = TRUE)

print(moran_test)

## 
##  Moran I test under randomisation
## 
## data:  data_sf$y_var  
## weights: w_list    
## 
## Moran I statistic standard deviate = 0.83258, p-value = 0.2025
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.05473869       -0.02941176        0.01021554

# Ekstrak nilai
moran_i <- moran_test$estimate[1]
expected_i <- moran_test$estimate[2]
variance_i <- moran_test$estimate[3]
p_value <- moran_test$p.value

cat("\n=== RINGKASAN MORAN'S I ===\n")

## 
## === RINGKASAN MORAN'S I ===

cat("Moran's I         :", round(moran_i, 4), "\n")

## Moran's I         : 0.0547

cat("Expected I        :", round(expected_i, 4), "\n")

## Expected I        : -0.0294

cat("Variance          :", round(variance_i, 6), "\n")

## Variance          : 0.010216

cat("Z-score           :", round(moran_test$statistic, 4), "\n")

## Z-score           : 0.8326

cat("P-value           :", format.pval(p_value, digits = 4), "\n")

## P-value           : 0.2025

cat("Kesimpulan        :", 
    ifelse(p_value < 0.05, 
           ifelse(moran_i > 0, "Autokorelasi Positif Signifikan", 
                  "Autokorelasi Negatif Signifikan"),
           "Tidak Ada Autokorelasi Signifikan"), "\n")

## Kesimpulan        : Tidak Ada Autokorelasi Signifikan

Interpretasi:

• [Jika I > 0 dan p < 0.05]: Terdapat autokorelasi spasial positif yang signifikan, artinya wilayah dengan nilai tinggi cenderung bertetangga dengan wilayah nilai tinggi (clustering), dan wilayah nilai rendah bertetangga dengan wilayah nilai rendah.

• [Jika I < 0 dan p < 0.05]: Terdapat autokorelasi spasial negatif yang signifikan, artinya wilayah dengan nilai tinggi cenderung bertetangga dengan wilayah nilai rendah (dispersi).

• [Jika p > 0.05]: Tidak terdapat autokorelasi spasial yang signifikan, distribusi cenderung acak (random).

Berdasarkan perhitungan tersebut, tidak terdapat autokorelasi spasial global yang signifikan pada distribusi kasus leptospirosis antar kabupaten/kota. Pola sebaran kasus tidak berbeda nyata dari pola acak secara global di Jawa Tengah

# Moran scatterplot
moran.plot(data_sf$y_var, w_list, 
           labels = FALSE,
           xlab = "Variabel Y (Standardized)",
           ylab = "Spatial Lag of Y (Standardized)",
           main = "Moran's I Scatterplot",
           zero.policy = TRUE)

Gambar 4. Moran’s I Scatterplot

Interpretasi:

Moran’s I scatterplot menunjukkan kemiringan garis regresi yang kecil dan positif, sejalan dengan nilai Moran’s I sebesar 0,0547. Sebaran titik yang relatif acak di keempat kuadran mengindikasikan bahwa hubungan antara nilai kasus leptospirosis suatu wilayah dengan rata-rata tetangganya bersifat lemah. Hal ini menegaskan bahwa tidak terdapat pola pengelompokan spasial global yang kuat, sesuai dengan hasil uji Moran’s I yang tidak signifikan secara statistik (p-value = 0,2025).

4.5.2 Geary’s C

# Hitung Geary's C
geary_test <- geary.test(data_sf$y_var, w_list, zero.policy = TRUE)

print(geary_test)

## 
##  Geary C test under randomisation
## 
## data:  data_sf$y_var 
## weights: w_list   
## 
## Geary C statistic standard deviate = 0.70337, p-value = 0.2409
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.91942972        1.00000000        0.01312163

geary_c <- geary_test$estimate[1]
geary_p <- geary_test$p.value

cat("\n=== RINGKASAN GEARY'S C ===\n")

## 
## === RINGKASAN GEARY'S C ===

cat("Geary's C         :", round(geary_c, 4), "\n")

## Geary's C         : 0.9194

cat("Expected C        :", 1, "\n")

## Expected C        : 1

cat("Z-score           :", round(geary_test$statistic, 4), "\n")

## Z-score           : 0.7034

cat("P-value           :", format.pval(geary_p, digits = 4), "\n")

## P-value           : 0.2409

cat("Kesimpulan        :", 
    ifelse(geary_p < 0.05,
           ifelse(geary_c < 1, "Autokorelasi Positif Signifikan",
                  "Autokorelasi Negatif Signifikan"),
           "Tidak Ada Autokorelasi Signifikan"), "\n")

## Kesimpulan        : Tidak Ada Autokorelasi Signifikan

Interpretasi:

Tidak terdapat autokorelasi spasial yang signifikan pada variabel kasus leptospirosis berdasarkan uji Geary’s C. Pola nilai antar wilayah cenderung acak, dan kemiripan antar tetangga tidak cukup kuat secara statistik.

4.6 Analisis LISA (Local Indicators of Spatial Association)

# Hitung Local Moran's I
local_moran <- localmoran(data_sf$y_var, w_list, zero.policy = TRUE)

# Tambahkan ke data
data_sf$local_i <- local_moran[, 1]
data_sf$local_p <- local_moran[, 5]

# Standardize variable
data_sf$y_std <- scale(data_sf$y_var)

# Hitung spatial lag
data_sf$lag_y_std <- lag.listw(w_list, data_sf$y_std, zero.policy = TRUE)

# Klasifikasi quadrant
data_sf <- data_sf %>%
  mutate(
    quadrant = case_when(
      y_std > 0 & lag_y_std > 0 ~ "High-High",
      y_std < 0 & lag_y_std < 0 ~ "Low-Low",
      y_std < 0 & lag_y_std > 0 ~ "Low-High",
      y_std > 0 & lag_y_std < 0 ~ "High-Low",
      TRUE ~ "Not Significant"
    ),
    lisa_cluster = ifelse(local_p < 0.05, quadrant, "Not Significant")
  )

# Hitung jumlah per kategori
lisa_summary <- data_sf %>%
  st_drop_geometry() %>%
  count(lisa_cluster) %>%
  mutate(pct = n/sum(n)*100)

cat("=== RINGKASAN LISA CLUSTER ===\n")

## === RINGKASAN LISA CLUSTER ===

print(lisa_summary)

##      lisa_cluster  n       pct
## 1       High-High  2  5.714286
## 2 Not Significant 33 94.285714

# Identifikasi wilayah per cluster
cat("\n\nWilayah High-High (Hotspot):\n")

## 
## 
## Wilayah High-High (Hotspot):

print(data_sf %>% filter(lisa_cluster == "High-High") %>% 
        select(kabkota, y_var, x1, x2) %>% st_drop_geometry())

##       kabkota y_var      x1 x2
## 1 Kab.Cilacap 18.75 2017.55 38
## 2 Kab.Kebumen 47.75 1405.25 35

cat("\n\nWilayah Low-Low (Coldspot):\n")

## 
## 
## Wilayah Low-Low (Coldspot):

print(data_sf %>% filter(lisa_cluster == "Low-Low") %>% 
        select(kabkota, y_var, x1, x2) %>% st_drop_geometry())

## [1] kabkota y_var   x1      x2     
## <0 rows> (or 0-length row.names)

# Peta LISA
tm_shape(data_sf) +
  tm_polygons("lisa_cluster",
              palette = c("High-High" = "#d7191c",
                         "Low-Low" = "#2b83ba",
                         "Low-High" = "#abdda4",
                         "High-Low" = "#fdae61",
                         "Not Significant" = "#f0f0f0"),
              title = "LISA Cluster",
              border.col = "white",
              border.alpha = 0.5) +
  tm_layout(main.title = "Klaster Lokal Kasus Leptospirosis (LISA)",
            main.title.size = 1.2,
            legend.outside = TRUE,
            frame = FALSE) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_compass(position = c("right", "top")) +
  tm_credits("α = 0.05", position = c("right", "bottom"))

Gambar 5. Peta LISA Cluster Kasus Leptospirosis

Interpretasi LISA:

Hasil analisis LISA menunjukkan adanya klaster High–High kasus leptospirosis yang signifikan secara statistik (α = 0,05). Wilayah tersebut memiliki tingkat kasus leptospirosis yang tinggi dan dikelilingi oleh wilayah dengan tingkat kasus yang juga tinggi, sehingga mengindikasikan adanya autokorelasi spasial positif. Sementara itu, sebagian besar wilayah lainnya tidak menunjukkan pola klaster yang signifikan.

# Peta p-value
tm_shape(data_sf) +
  tm_polygons("local_p",
              breaks = c(0, 0.01, 0.05, 0.1, 1),
              palette = "-RdYlGn",
              title = "P-value",
              border.col = "white",
              labels = c("< 0.01", "0.01-0.05", "0.05-0.10", "> 0.10")) +
  tm_layout(main.title = "Signifikansi Statistik Local Moran's I",
            legend.outside = TRUE,
            frame = FALSE)

Gambar 6. Peta Signifikansi Statistik LISA

Interpretasi:

Peta signifikansi Local Moran’s I menunjukkan apakah wilayah tersebut valid secara statistik. Peta tersebut menunjukkan bahwa hanya sebagian wilayah di Jawa Tengah, khususnya di bagian barat–selatan, yang memiliki autokorelasi spasial lokal kasus leptospirosis yang signifikan (p < 0,05). Sebagian besar wilayah lainnya tidak menunjukkan signifikansi statistik, mengindikasikan bahwa pola kejadian leptospirosis di wilayah tersebut bersifat acak.

4.7 Pemodelan Spatial Econometrics

4.7.1 Model OLS (Ordinary Least Squares) - Baseline

# Model OLS tanpa efek spasial
formula_ols <- y_var ~ x1 + x2

model_ols <- lm(formula_ols, data = data_sf)

# Summary model
summary(model_ols)

## 
## Call:
## lm(formula = formula_ols, data = data_sf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.355 -11.768  -3.412   8.109  44.499 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1.441177   9.777788  -0.147   0.8837  
## x1           0.012813   0.007145   1.793   0.0824 .
## x2           0.086507   0.373395   0.232   0.8183  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.85 on 32 degrees of freedom
## Multiple R-squared:  0.121,  Adjusted R-squared:  0.06603 
## F-statistic: 2.202 on 2 and 32 DF,  p-value: 0.1271

# Ekstrak hasil
ols_coef <- coef(model_ols)
ols_r2 <- summary(model_ols)$r.squared
ols_adj_r2 <- summary(model_ols)$adj.r.squared
ols_aic <- AIC(model_ols)

cat("\n=== RINGKASAN MODEL OLS ===\n")

## 
## === RINGKASAN MODEL OLS ===

cat("R-squared        :", round(ols_r2, 4), "\n")

## R-squared        : 0.121

cat("Adjusted R2      :", round(ols_adj_r2, 4), "\n")

## Adjusted R2      : 0.066

cat("AIC              :", round(ols_aic, 2), "\n")

## AIC              : 305.91

cat("Koefisien:\n")

## Koefisien:

print(round(ols_coef, 4))

## (Intercept)          x1          x2 
##     -1.4412      0.0128      0.0865

Interpretasi Model OLS:

Model OLS baseline menghasilkan:

• R² = 0.121: Model menjelaskan 12.1% variasi kasus leptospirosis
• Koefisien x1 (Jumlah Penduduk): Setiap peningkatan 1.000 penduduk berhubungan dengan perubahan 0.0128 kasus leptospirosis, tidak signifikan pada α = 0.05
• Koefisien x2 (Jumlah Puskesmas): Setiap penambahan 1 puskesmas berhubungan dengan perubahan 0.0865 kasus, tidak signifikan pada α = 0.05

Namun, model OLS mengabaikan spatial dependence yang telah terbukti ada (Moran’s I signifikan). Kita perlu menguji apakah spatial lag atau spatial error lebih sesuai.

4.7.2 Uji Spatial Dependence (Lagrange Multiplier Test)

# LM test untuk spatial dependence
lm_tests <- lm.LMtests(model_ols, w_list, test = "all")

print(lm_tests)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = formula_ols, data = data_sf)
## test weights: listw
## 
## RSerr = 0.28134, df = 1, p-value = 0.5958
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = formula_ols, data = data_sf)
## test weights: listw
## 
## RSlag = 0.094877, df = 1, p-value = 0.7581
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = formula_ols, data = data_sf)
## test weights: listw
## 
## adjRSerr = 0.79958, df = 1, p-value = 0.3712
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = formula_ols, data = data_sf)
## test weights: listw
## 
## adjRSlag = 0.61312, df = 1, p-value = 0.4336
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = formula_ols, data = data_sf)
## test weights: listw
## 
## SARMA = 0.89446, df = 2, p-value = 0.6394

# Ekstrak p-values
lm_lag_p <- lm_tests$LMlag$p.value
lm_err_p <- lm_tests$LMerr$p.value
lm_rlag_p <- lm_tests$RLMlag$p.value
lm_rerr_p <- lm_tests$RLMerr$p.value
lm_sarma_p <- lm_tests$SARMA$p.value

cat("\n=== RINGKASAN LM TESTS ===\n")

## 
## === RINGKASAN LM TESTS ===

cat("LM lag (p-value)         :", format.pval(lm_lag_p, digits = 4), 
    ifelse(lm_lag_p < 0.05, " ***", ""), "\n")

## LM lag (p-value)         :

cat("Robust LM lag (p-value)  :", format.pval(lm_rlag_p, digits = 4),
    ifelse(lm_rlag_p < 0.05, " ***", ""), "\n")

## Robust LM lag (p-value)  :

cat("LM error (p-value)       :", format.pval(lm_err_p, digits = 4),
    ifelse(lm_err_p < 0.05, " ***", ""), "\n")

## LM error (p-value)       :

cat("Robust LM error (p-value):", format.pval(lm_rerr_p, digits = 4),
    ifelse(lm_rerr_p < 0.05, " ***", ""), "\n")

## Robust LM error (p-value):

cat("LM SARMA (p-value)       :", format.pval(lm_sarma_p, digits = 4),
    ifelse(lm_sarma_p < 0.05, " ***", ""), "\n")

## LM SARMA (p-value)       : 0.6394

str(lm_tests)

## List of 5
##  $ RSerr   :List of 5
##   ..$ statistic: Named num 0.281
##   .. ..- attr(*, "names")= chr "RSerr"
##   ..$ parameter: Named num 1
##   .. ..- attr(*, "names")= chr "df"
##   ..$ p.value  : Named num 0.596
##   .. ..- attr(*, "names")= chr ""
##   ..$ method   : chr "Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence"
##   ..$ data.name: chr "\nmodel: lm(formula = formula_ols, data = data_sf)\ntest weights: listw\n"
##   ..- attr(*, "class")= chr "htest"
##  $ RSlag   :List of 5
##   ..$ statistic: Named num 0.0949
##   .. ..- attr(*, "names")= chr "RSlag"
##   ..$ parameter: Named num 1
##   .. ..- attr(*, "names")= chr "df"
##   ..$ p.value  : Named num 0.758
##   .. ..- attr(*, "names")= chr ""
##   ..$ method   : chr "Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence"
##   ..$ data.name: chr "\nmodel: lm(formula = formula_ols, data = data_sf)\ntest weights: listw\n"
##   ..- attr(*, "class")= chr "htest"
##  $ adjRSerr:List of 5
##   ..$ statistic: Named num 0.8
##   .. ..- attr(*, "names")= chr "adjRSerr"
##   ..$ parameter: Named num 1
##   .. ..- attr(*, "names")= chr "df"
##   ..$ p.value  : Named num 0.371
##   .. ..- attr(*, "names")= chr ""
##   ..$ method   : chr "Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence"
##   ..$ data.name: chr "\nmodel: lm(formula = formula_ols, data = data_sf)\ntest weights: listw\n"
##   ..- attr(*, "class")= chr "htest"
##  $ adjRSlag:List of 5
##   ..$ statistic: Named num 0.613
##   .. ..- attr(*, "names")= chr "adjRSlag"
##   ..$ parameter: Named num 1
##   .. ..- attr(*, "names")= chr "df"
##   ..$ p.value  : Named num 0.434
##   .. ..- attr(*, "names")= chr ""
##   ..$ method   : chr "Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence"
##   ..$ data.name: chr "\nmodel: lm(formula = formula_ols, data = data_sf)\ntest weights: listw\n"
##   ..- attr(*, "class")= chr "htest"
##  $ SARMA   :List of 5
##   ..$ statistic: Named num 0.894
##   .. ..- attr(*, "names")= chr "SARMA"
##   ..$ parameter: Named num 2
##   .. ..- attr(*, "names")= chr "df"
##   ..$ p.value  : Named num 0.639
##   .. ..- attr(*, "names")= chr ""
##   ..$ method   : chr "Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence"
##   ..$ data.name: chr "\nmodel: lm(formula = formula_ols, data = data_sf)\ntest weights: listw\n"
##   ..- attr(*, "class")= chr "htest"
##  - attr(*, "class")= chr "RStestlist"

names(lm_tests)

## [1] "RSerr"    "RSlag"    "adjRSerr" "adjRSlag" "SARMA"

lm_lag_p   <- lm_tests$RSlag$p.value
lm_err_p   <- lm_tests$RSerr$p.value
lm_rlag_p  <- lm_tests$adjRSlag$p.value
lm_rerr_p  <- lm_tests$adjRSerr$p.value
lm_sarma_p <- lm_tests$SARMA$p.value

cat("\n=== RINGKASAN LM TESTS ===\n")

## 
## === RINGKASAN LM TESTS ===

cat("LM lag (p-value)         :", format.pval(lm_lag_p, digits = 4), "\n")

## LM lag (p-value)         : 0.7581

cat("Robust LM lag (p-value)  :", format.pval(lm_rlag_p, digits = 4), "\n")

## Robust LM lag (p-value)  : 0.4336

cat("LM error (p-value)       :", format.pval(lm_err_p, digits = 4), "\n")

## LM error (p-value)       : 0.5958

cat("Robust LM error (p-value):", format.pval(lm_rerr_p, digits = 4), "\n")

## Robust LM error (p-value): 0.3712

cat("LM SARMA (p-value)       :", format.pval(lm_sarma_p, digits = 4), "\n")

## LM SARMA (p-value)       : 0.6394

cat("\n=== REKOMENDASI MODEL ===\n")

## 
## === REKOMENDASI MODEL ===

if(lm_rlag_p < 0.05 & lm_rerr_p >= 0.05) {
  cat("Rekomendasi: Gunakan SAR (Spatial Lag Model)\n")
} else if(lm_rerr_p < 0.05 & lm_rlag_p >= 0.05) {
  cat("Rekomendasi: Gunakan SEM (Spatial Error Model)\n")
} else if(lm_rlag_p < 0.05 & lm_rerr_p < 0.05) {
  cat("Rekomendasi: Gunakan SDM (Spatial Durbin Model) atau SAC\n")
} else {
  cat("Rekomendasi: Model OLS mungkin sudah cukup\n")
}

## Rekomendasi: Model OLS mungkin sudah cukup

Kesimpulan :

Meskipun hasil Moran’s I, Geary’s C, dan LM tests menunjukkan tidak adanya autokorelasi spasial yang signifikan, penelitian ini tetap mengestimasi model SAR, SEM, dan/atau SDM sebagai analisis pembanding untuk melihat apakah terdapat peningkatan kinerja model dibandingkan OLS

4.7.3 Model SAR (Spatial Autoregressive Model)

# Estimasi SAR
model_sar <- lagsarlm(formula_ols, data = data_sf, listw = w_list, 
                      zero.policy = TRUE)

# Summary
summary(model_sar, Nagelkerke = TRUE)

## 
## Call:lagsarlm(formula = formula_ols, data = data_sf, listw = w_list, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.7142 -11.6589  -3.3297   7.5815  44.4236 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2285806  9.8712126 -0.2258  0.82138
## x1           0.0128092  0.0068202  1.8781  0.06036
## x2           0.0631910  0.3562512  0.1774  0.85921
## 
## Rho: 0.088492, LR test value: 0.11027, p-value: 0.73984
## Asymptotic standard error: 0.24055
##     z-value: 0.36787, p-value: 0.71297
## Wald statistic: 0.13533, p-value: 0.71297
## 
## Log likelihood: -148.9015 for lag model
## ML residual variance (sigma squared): 289.8, (sigma: 17.024)
## Nagelkerke pseudo-R-squared: 0.12374 
## Number of observations: 35 
## Number of parameters estimated: 5 
## AIC: 307.8, (AIC for lm: 305.91)
## LM test for residual autocorrelation
## test value: 1.1023, p-value: 0.29376

# Ekstrak hasil
sar_rho <- model_sar$rho
sar_aic <- AIC(model_sar)
sar_bic <- BIC(model_sar)
sar_loglik <- logLik(model_sar)
sar_pseudo_r2 <- summary(model_sar, Nagelkerke = TRUE)$NK

cat("\n=== RINGKASAN MODEL SAR ===\n")

## 
## === RINGKASAN MODEL SAR ===

cat("Rho (ρ)          :", round(sar_rho, 4), "\n")

## Rho (ρ)          : 0.0885

cat("Pseudo R²        :", round(sar_pseudo_r2, 4), "\n")

## Pseudo R²        : 0.1237

cat("AIC              :", round(sar_aic, 2), "\n")

## AIC              : 307.8

cat("BIC              :", round(sar_bic, 2), "\n")

## BIC              : 315.58

cat("Log-likelihood   :", round(as.numeric(sar_loglik), 2), "\n")

## Log-likelihood   : -148.9

Interpretasi Model SAR:

Model SAR menangkap spatial spillover effect dengan persamaan :

\(y_i = \rho \sum_j w_{ij} y_j + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_i\)

Model Spatial Autoregressive (SAR) digunakan untuk menguji apakah terdapat pengaruh nilai kasus leptospirosis di suatu wilayah terhadap wilayah tetangganya. Hasil model SAR menunjukkan bahwa pengaruh spasial global tidak signifikan. Hal ini menegaskan bahwa kejadian leptospirosis di Jawa Tengah lebih dipengaruhi oleh karakteristik lokal wilayah dibandingkan efek spillover antarwilayah.

4.7.4 Model SEM (Spatial Error Model)

# Estimasi SEM
model_sem <- errorsarlm(formula_ols, data = data_sf, listw = w_list,
                        zero.policy = TRUE)

# Summary
summary(model_sem, Nagelkerke = TRUE)

## 
## Call:errorsarlm(formula = formula_ols, data = data_sf, listw = w_list, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.1271 -11.1673  -2.1158   7.0468  42.6718 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -0.37189480  9.47727598 -0.0392  0.96870
## x1           0.01381132  0.00677398  2.0389  0.04146
## x2          -0.00081442  0.34989298 -0.0023  0.99814
## 
## Lambda: 0.17717, LR test value: 0.37313, p-value: 0.5413
## Asymptotic standard error: 0.23784
##     z-value: 0.74492, p-value: 0.45632
## Wald statistic: 0.55491, p-value: 0.45632
## 
## Log likelihood: -148.7701 for error model
## ML residual variance (sigma squared): 286.22, (sigma: 16.918)
## Nagelkerke pseudo-R-squared: 0.13029 
## Number of observations: 35 
## Number of parameters estimated: 5 
## AIC: 307.54, (AIC for lm: 305.91)

# Ekstrak hasil
sem_lambda <- model_sem$lambda
sem_aic <- AIC(model_sem)
sem_bic <- BIC(model_sem)
sem_loglik <- logLik(model_sem)
sem_pseudo_r2 <- summary(model_sem, Nagelkerke = TRUE)$NK

cat("\n=== RINGKASAN MODEL SEM ===\n")

## 
## === RINGKASAN MODEL SEM ===

cat("Lambda (λ)       :", round(sem_lambda, 4), "\n")

## Lambda (λ)       : 0.1772

cat("Pseudo R²        :", round(sem_pseudo_r2, 4), "\n")

## Pseudo R²        : 0.1303

cat("AIC              :", round(sem_aic, 2), "\n")

## AIC              : 307.54

cat("BIC              :", round(sem_bic, 2), "\n")

## BIC              : 315.32

cat("Log-likelihood   :", round(as.numeric(sem_loglik), 2), "\n")

## Log-likelihood   : -148.77

Interpretasi Model SEM:

Model SEM menangkap spatial dependence pada error term:

\(y_i = \beta_1 x_{1i} + \beta_2 x_{2i} + u_i\) \(u_i = \lambda \sum_j w_{ij} u_j + \epsilon_i\)

Model SEM menunjukkan bahwa tidak terdapat ketergantungan spasial yang signifikan pada komponen error kasus leptospirosis di Jawa Tengah. Meskipun salah satu variabel independen berpengaruh signifikan, model SEM tidak memberikan peningkatan kinerja yang berarti dibandingkan model OLS.

4.7.5 Model SDM (Spatial Durbin Model)

# Estimasi SDM
model_sdm <- lagsarlm(formula_ols, data = data_sf, listw = w_list,
                      type = "mixed", zero.policy = TRUE)

# Summary
summary(model_sdm, Nagelkerke = TRUE)

## 
## Call:lagsarlm(formula = formula_ols, data = data_sf, listw = w_list, 
##     type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.4237 -12.4095  -4.4841  10.0453  41.8234 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -14.6805349  20.0978219 -0.7305  0.46511
## x1            0.0116266   0.0065451  1.7764  0.07567
## x2            0.1662768   0.3409782  0.4876  0.62580
## lag.x1       -0.0289227   0.0139748 -2.0696  0.03849
## lag.x2        1.6621811   0.8390971  1.9809  0.04760
## 
## Rho: 0.1567, LR test value: 0.34751, p-value: 0.55553
## Asymptotic standard error: 0.23494
##     z-value: 0.66698, p-value: 0.50478
## Wald statistic: 0.44486, p-value: 0.50478
## 
## Log likelihood: -146.3513 for mixed model
## ML residual variance (sigma squared): 249.64, (sigma: 15.8)
## Nagelkerke pseudo-R-squared: 0.24257 
## Number of observations: 35 
## Number of parameters estimated: 7 
## AIC: 306.7, (AIC for lm: 305.05)
## LM test for residual autocorrelation
## test value: 1.6173, p-value: 0.20347

# Ekstrak hasil
sdm_aic <- AIC(model_sdm)
sdm_bic <- BIC(model_sdm)
sdm_loglik <- logLik(model_sdm)
sdm_pseudo_r2 <- summary(model_sdm, Nagelkerke = TRUE)$NK

cat("\n=== RINGKASAN MODEL SDM ===\n")

## 
## === RINGKASAN MODEL SDM ===

cat("Pseudo R²        :", round(sdm_pseudo_r2, 4), "\n")

## Pseudo R²        : 0.2426

cat("AIC              :", round(sdm_aic, 2), "\n")

## AIC              : 306.7

cat("BIC              :", round(sdm_bic, 2), "\n")

## BIC              : 317.59

cat("Log-likelihood   :", round(as.numeric(sdm_loglik), 2), "\n")

## Log-likelihood   : -146.35

# Impact estimates (Direct, Indirect, Total)
impacts_sdm <- impacts(model_sdm, listw = w_list, R = 500)
cat("\n=== IMPACT ESTIMATES (SDM) ===\n")

## 
## === IMPACT ESTIMATES (SDM) ===

summary(impacts_sdm, zstats = TRUE)

## Impact measures (mixed, exact):
##        Direct    Indirect       Total
## x1 0.01073168 -0.03124173 -0.02051005
## x2 0.22202595  1.94619221  2.16821816
## ========================================================
## Simulation results ( variance matrix):
## Direct:
## 
## Iterations = 1:500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean       SD  Naive SE Time-series SE
## x1 0.01041 0.006825 0.0003052      0.0002705
## x2 0.22656 0.362600 0.0162160      0.0148897
## 
## 2. Quantiles for each variable:
## 
##         2.5%       25%     50%     75%   97.5%
## x1 -0.003135  0.005753 0.01084 0.01484 0.02353
## x2 -0.467390 -0.017420 0.21435 0.45921 0.97001
## 
## ========================================================
## Indirect:
## 
## Iterations = 1:500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean    SD  Naive SE Time-series SE
## x1 -0.03378 0.022 0.0009838      0.0009838
## x2  2.14726 1.397 0.0624889      0.0574128
## 
## 2. Quantiles for each variable:
## 
##        2.5%      25%      50%      75%    97.5%
## x1 -0.09192 -0.04436 -0.03096 -0.02034 0.003922
## x2 -0.11188  1.37509  2.03176  2.74437 4.880112
## 
## ========================================================
## Total:
## 
## Iterations = 1:500
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean      SD Naive SE Time-series SE
## x1 -0.02337 0.02471 0.001105       0.001028
## x2  2.37381 1.58659 0.070955       0.061748
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%      50%       75%   97.5%
## x1 -0.09540 -0.034 -0.02125 -0.008541 0.01742
## x2 -0.06352  1.447  2.28636  3.087532 5.51719
## 
## ========================================================
## Simulated standard errors
##         Direct   Indirect      Total
## x1 0.006825108 0.02199777 0.02470529
## x2 0.362600108 1.39729482 1.58659409
## 
## Simulated z-values:
##       Direct  Indirect      Total
## x1 1.5250202 -1.535467 -0.9458877
## x2 0.6248096  1.536724  1.4961685
## 
## Simulated p-values:
##    Direct  Indirect Total  
## x1 0.12725 0.12467  0.34421
## x2 0.53210 0.12436  0.13461

Interpretasi Model SDM:

Model SDM adalah model paling general yang mencakup spatial lag pada Y dan X:

\(y_i = \rho \sum_j w_{ij} y_j + \beta_1 x_{1i} + \beta_2 x_{2i} + \theta_1 \sum_j w_{ij} x_{1j} + \theta_2 \sum_j w_{ij} x_{2j} + \epsilon_i\)

Model SDM menunjukkan bahwa meskipun tidak terdapat ketergantungan spasial langsung pada kasus leptospirosis antarwilayah, terdapat efek spillover yang signifikan melalui variabel independen.

4.7.6 Perbandingan Model

# Tabel perbandingan
comparison_table <- data.frame(
  Model = c("OLS", "SAR", "SEM", "SDM"),
  AIC = c(ols_aic, sar_aic, sem_aic, sdm_aic),
  BIC = c(BIC(model_ols), sar_bic, sem_bic, sdm_bic),
  LogLik = c(as.numeric(logLik(model_ols)), 
             as.numeric(sar_loglik),
             as.numeric(sem_loglik),
             as.numeric(sdm_loglik)),
  R2 = c(ols_r2, sar_pseudo_r2, sem_pseudo_r2, sdm_pseudo_r2)
) %>%
  mutate(
    Delta_AIC = AIC - min(AIC),
    Rank_AIC = rank(AIC)
  ) %>%
  arrange(AIC)

comparison_table %>%
  kable(caption = "Tabel 2. Perbandingan Model Spatial Econometrics",
        digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE) %>%
  row_spec(which.min(comparison_table$AIC), bold = TRUE, color = "white", background = "#2b8cbe")

Tabel 2. Perbandingan Model Spatial Econometrics

Model	AIC	BIC	LogLik	R2	Delta_AIC	Rank_AIC
OLS	305.91	312.13	-148.96	0.12	0.00	1
SDM	306.70	317.59	-146.35	0.24	0.79	2
SEM	307.54	315.32	-148.77	0.13	1.63	3
SAR	307.80	315.58	-148.90	0.12	1.89	4

# Visualisasi perbandingan AIC
ggplot(comparison_table, aes(x = reorder(Model, AIC), y = AIC, fill = Model)) +
  geom_col() +
  geom_text(aes(label = round(AIC, 1)), vjust = -0.5) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Perbandingan AIC Antar Model",
       subtitle = "Semakin rendah AIC, semakin baik model",
       x = "Model", y = "AIC") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", size = 14))

Interpretasi Perbandingan Model:

Dibandingkan model OLS, SAR, SEM, dan SDM, model OLS menunjukkan AIC terkecil sehingga kinerja OLS lebih baik dalam menjelaskan variasi kasus leptospirosis.

4.7.7 Visualisasi Residual Model Terpilih

# Pilih model terbaik berdasarkan AIC
best_model_name <- comparison_table$Model[1]

if(best_model_name == "SAR") {
  best_model <- model_sar
  residuals_best <- residuals(model_sar)
} else if(best_model_name == "SEM") {
  best_model <- model_sem
  residuals_best <- residuals(model_sem)
} else if(best_model_name == "SDM") {
  best_model <- model_sdm
  residuals_best <- residuals(model_sdm)
} else {
  best_model <- model_ols
  residuals_best <- residuals(model_ols)
}

# Tambahkan residual ke data
data_sf$residual <- residuals_best

# Plot residual
p1 <- ggplot(data.frame(residual = residuals_best), aes(x = residual)) +
  geom_histogram(bins = 20, fill = "#4292c6", alpha = 0.7, color = "white") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Distribusi Residual", x = "Residual", y = "Frekuensi") +
  theme_minimal()

p2 <- ggplot(data.frame(fitted = fitted(best_model), residual = residuals_best),
             aes(x = fitted, y = residual)) +
  geom_point(alpha = 0.6, color = "#4292c6") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "loess", se = FALSE, color = "blue") +
  labs(title = "Residual vs Fitted", x = "Fitted Values", y = "Residual") +
  theme_minimal()

p3 <- ggplot(data.frame(residual = residuals_best), aes(sample = residual)) +
  stat_qq(color = "#4292c6") +
  stat_qq_line(color = "red", linetype = "dashed") +
  labs(title = "Q-Q Plot", x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()

grid.arrange(p1, p2, p3, ncol = 2, nrow = 2)

Gambar 7. Analisis Residual Model Terpilih

# Test Moran's I pada residual
moran_residual <- moran.test(residuals_best, w_list, zero.policy = TRUE)
cat("\n=== MORAN'S I PADA RESIDUAL ===\n")

## 
## === MORAN'S I PADA RESIDUAL ===

cat("Moran's I :", round(moran_residual$estimate[1], 4), "\n")

## Moran's I : 0.0587

cat("P-value   :", format.pval(moran_residual$p.value, digits = 4), "\n")

## P-value   : 0.1925

cat("Kesimpulan:", ifelse(moran_residual$p.value > 0.05,
                          "Tidak ada autokorelasi spasial pada residual (BAIK)",
                          "Masih ada autokorelasi spasial pada residual"), "\n")

## Kesimpulan: Tidak ada autokorelasi spasial pada residual (BAIK)

Interpretasi :

Hasil ini dikategorikan sebagai BAIK dalam pemodelan spasial. Hal ini membuktikan bahwa model OLS yang Anda gunakan sudah mampu menangkap seluruh pola spasial dalam data Kasus Leptospirosis Jawa Tengah, sehingga sisaan (error) yang dihasilkan bersifat acak (random) dan tidak saling bergantung antar wilayah.

# Peta residual
tm_shape(data_sf) +
  tm_polygons("residual",
              palette = "RdBu",
              style = "jenks",
              midpoint = 0,
              n = 7,
              title = "Residual",
              border.col = "white") +
  tm_layout(main.title = paste("Peta Residual Model", best_model_name),
            legend.outside = TRUE,
            frame = FALSE) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_compass(position = c("right", "top"))

Gambar 8. Peta Residual Model

Interpretasi Peta Residual:

• Residual positif (merah): Wilayah dengan kasus leptospirosis lebih tinggi dari prediksi model (under-prediction)
• Residual negatif (biru): Wilayah dengan kasus lebih rendah dari prediksi (over-prediction)
• Residual mendekati nol (putih): Prediksi model akurat

5 KESIMPULAN DAN SARAN

5.1 Kesimpulan

Berdasarkan rangkaian analisis statistik spasial yang telah dilakukan terhadap data Leptospirosis di Provinsi Jawa Tengah periode 2022-2025, dapat ditarik beberapa kesimpulan utama sebagai berikut:

1. Pola Distribusi Spasial: Persebaran kasus Leptospirosis di Jawa Tengah menunjukkan adanya konsentrasi kasus yang fluktuatif, dengan tahun 2024 menjadi tahun puncak beban penyakit (258 kasus). Identifikasi melalui peta residual dan interpolasi menunjukkan wilayah pesisir utara (seperti Demak) dan wilayah tengah (Banyumas) memiliki risiko relatif lebih tinggi dibandingkan wilayah lainnya.

2. Hasil Uji Autokorelasi: Meskipun terdapat indikasi visual adanya pengelompokan kasus, hasil uji Geary’s C menunjukkan nilai 0,9194 dengan \(p\text{-value}\) 0,2409, yang berarti pola persebaran kasus secara lokal cenderung bersifat acak (random) pada tingkat signifikansi 5%. Hal ini mengindikasikan bahwa keterkaitan antar wilayah bertetangga tidak selalu linear atau seragam di seluruh kabupaten.

3. Pemilihan Model Terbaik: Melalui perbandingan nilai Akaike Information Criterion (AIC), model Ordinary Least Square (OLS) terpilih sebagai model yang paling efisien dengan nilai AIC terendah sebesar 305,9. Penambahan bobot spasial pada model (SAR, SEM, dan SDM) tidak memberikan peningkatan kebaikan suai (goodness of fit) yang signifikan dibandingkan model linear klasik untuk dataset ini.

4. Validitas : Model OLS terbukti valid secara spasial karena uji Moran’s I pada residual menghasilkan nilai 0,0587 dengan \(p\text{-value}\) 0,1925.

5.2 Saran

5.2.1 Saran Akademik dan Penelitian Selanjutnya

• Penambahan Variabel Lingkungan: Penelitian selanjutnya disarankan untuk menyertakan variabel meteorologi dan lingkungan yang lebih spesifik, seperti curah hujan, frekuensi banjir, atau tutupan lahan, guna meningkatkan nilai \(R^2\) dan akurasi model.

• Penggunaan Interpolasi (IDW): Disarankan untuk menerapkan teknik Inverse Distance Weighting (IDW) yang divalidasi dengan Cross Validation guna memetakan gradien risiko secara kontinu yang tidak dibatasi oleh batas administratif kabupaten.

• Analisis Temporal: Perlu dilakukan analisis spasio-temporal (seperti Spatial Panel Data) untuk melihat bagaimana tren perpindahan klaster penyakit berubah dari tahun 2022 hingga 2025 secara simultan.

5.2.2 Saran Kebijakan

1. Pencegahan dan Pengendalian Berbasis Wilayah:
   • Wilayah Hotspot (High-High): Prioritas tertinggi untuk program pencegahan intensif, perbaikan sanitasi lingkungan, pengendalian vektor, dan edukasi masyarakat
   • Wilayah Transisi: Surveilans aktif untuk deteksi dini dan pencegahan perluasan hotspot
2. Penguatan Sistem Surveilans:
   • Peningkatan kapasitas deteksi dan pelaporan kasus di wilayah dengan puskesmas terbatas
   • Implementasi early warning system berbasis spatial analysis
   • Integrasi data lingkungan dan klimatologi untuk prediksi risiko
3. Alokasi Sumber Daya:
   • Alokasi sumber daya kesehatan proporsional dengan tingkat risiko spasial
   • Penambahan puskesmas atau pos kesehatan di wilayah berisiko tinggi dengan akses terbatas
4. Penelitian Lanjutan:
   • Investigasi mendalam pada wilayah outlier (High-Low dan Low-High)
   • Studi kualitatif untuk memahami faktor lokal spesifik
   • Monitoring dan evaluasi dampak intervensi secara spasial

---

DAFTAR PUSTAKA

Anselin, L. (1988). Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers.

Bivand, R.S., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R. 2nd ed. New York: Springer.

Cressie, N. (1993). Statistics for Spatial Data. New York: John Wiley & Sons.

Elhorst, J.P. (2014). Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Berlin: Springer.

Getis, A., & Ord, J.K. (1992). The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis, 24(3), 189-206.

LeSage, J., & Pace, R.K. (2009). Introduction to Spatial Econometrics. Boca Raton: CRC Press.

Tobler, W.R. (1970). A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography, 46(sup1), 234-240.

---

LAMPIRAN

5.3 Lampiran A: Session Info

sessionInfo()

## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gstat_2.1-3         stringr_1.5.1       GWmodel_2.4-1      
##  [4] Rcpp_1.0.14         sp_2.2-0            robustbase_0.99-4-1
##  [7] gridExtra_2.3       scales_1.4.0        corrplot_0.95      
## [10] plotly_4.11.0       leaflet_2.2.3       RColorBrewer_1.1-3 
## [13] viridis_0.6.5       viridisLite_0.4.2   kableExtra_1.4.0   
## [16] knitr_1.49          tidyr_1.3.1         dplyr_1.1.4        
## [19] ggplot2_4.0.0       tmap_4.2            spatialreg_1.3-6   
## [22] Matrix_1.6-1.1      spdep_1.3-10        spData_2.3.4       
## [25] sf_1.0-20           readxl_1.4.5       
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.2.3            deldir_2.0-4         tmaptools_3.3       
##  [4] s2_1.1.7             logger_0.4.1         sandwich_3.1-1      
##  [7] rlang_1.1.4          magrittr_2.0.3       multcomp_1.4-29     
## [10] e1071_1.7-16         compiler_4.3.1       mgcv_1.8-42         
## [13] png_0.1-8            systemfonts_1.2.2    vctrs_0.6.5         
## [16] pkgconfig_2.0.3      wk_0.9.4             fastmap_1.2.0       
## [19] labeling_0.4.3       lwgeom_0.2-14        leafem_0.2.5        
## [22] utf8_1.2.4           rmarkdown_2.29       spacesXYZ_1.6-0     
## [25] purrr_1.0.4          xfun_0.52            cachem_1.1.0        
## [28] jsonlite_1.8.9       terra_1.8-29         parallel_4.3.1      
## [31] LearnBayes_2.15.1    R6_2.6.1             bslib_0.9.0         
## [34] stringi_1.8.4        boot_1.3-28.1        jquerylib_0.1.4     
## [37] cellranger_1.1.0     stars_0.6-8          zoo_1.8-12          
## [40] base64enc_0.1-3      FNN_1.1.4.1          splines_4.3.1       
## [43] tidyselect_1.2.1     rstudioapi_0.17.1    abind_1.4-8         
## [46] yaml_2.3.10          maptiles_0.10.0      codetools_0.2-19    
## [49] intervals_0.15.5     lattice_0.21-8       tibble_3.2.1        
## [52] leafsync_0.1.0       withr_3.0.2          S7_0.2.0            
## [55] coda_0.19-4.1        evaluate_1.0.3       survival_3.5-5      
## [58] units_0.8-7          proxy_0.4-27         xts_0.14.1          
## [61] xml2_1.3.8           pillar_1.10.1        KernSmooth_2.23-21  
## [64] generics_0.1.3       spacetime_1.3-3      class_7.3-23        
## [67] glue_1.8.0           lazyeval_0.2.2       tools_4.3.1         
## [70] leaflegend_1.2.1     data.table_1.16.4    mvtnorm_1.3-3       
## [73] XML_3.99-0.18        grid_4.3.1           crosstalk_1.2.1     
## [76] colorspace_2.1-1     nlme_3.1-162         cols4all_0.9        
## [79] raster_3.6-32        cli_3.6.1            svglite_2.1.3       
## [82] DEoptimR_1.1-4       gtable_0.3.6         sass_0.4.9          
## [85] digest_0.6.33        classInt_0.4-11      TH.data_1.1-4       
## [88] htmlwidgets_1.6.4    farver_2.1.2         htmltools_0.5.8.1   
## [91] lifecycle_1.0.4      httr_1.4.7           microbenchmark_1.5.0
## [94] MASS_7.3-60

5.4 Lampiran B: Link Sumber Daya

• Dashboard Interaktif: Dashboard Leptospirosis Jateng

• Video Presentasi:

• Data Source:

  Variabel Independen : Badan Pusat Statistik (BPS) Provinsi Jawa Tengah

  Variabel Dependen : Portal Data Jawa Tengah

5.5 Lampiran C: Acknowledgement

Analisis ini menggunakan bantuan AI tools (Claude AI by Anthropic) untuk:

- Menyusun struktur kode R dan RMarkdown

- Membantu visualisasi data dan peta

Namun, semua penjelasan, narasi, interpretasi, dan kesimpulan merupakan hasil pemahaman dan pemikiran mandiri penulis.

---

TERIMA KASIH

---

“In God we trust, all others must bring data” - W. Edwards Deming

“Geography matters” - Paul Krugman