Spasiaalll.knit

Analisis Spasial Kasus Campak di Kabupaten/Kota
Provinsi Jawa Timur Tahun 2022

Ditujukan Guna Memenuhi UAS Mata Kuliah Spatial Statistics

Disusun oleh:
Tiffanny Chantika Silalahi (140610230027)

Dosen Pengampu:
I Gede Nyoman Mindra Jaya, Ph.D

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

PENDAHULUAN

Analisis spasial menawarkan pendekatan yang kuat untuk memahami distribusi geografis penyakit menular dan determinan yang mendasarinya. Penelitian ini mengeksplorasi pola spasial kejadian campak di Jawa Timur tahun 2022 dan mengidentifikasi model regresi spasial yang paling sesuai untuk menjelaskan dependensi geografis antarwilayah. Data agregat kabupaten/kota—meliputi incidence rate campak, persentase kemiskinan, dan akses sanitasi layak—dianalisis menggunakan beberapa teknik spasial. Kami menguji autokorelasi spasial global dengan Moran’s I dan Geary’s C, serta melakukan analisis hotspot lokal menggunakan LISA dan Getis–Ord Gi*. Pemodelan dimulai dengan regresi linier klasik (OLS), dilanjutkan dengan enam model regresi spasial: Spatial Lag Model (SAR), Spatial Error Model (SEM), Spatial Durbin Model (SDM), Spatial Durbin Error Model (SDEM), Spatial Lag of X Model (SLX), dan Spatial Autoregressive Combined Model (SAC). Seleksi model optimal didasarkan pada Akaike Information Criterion (AIC) dan diagnostik residual. Kami juga menerapkan Geographically Weighted Regression (GWR) untuk mendeteksi variasi hubungan lokal, serta tiga metode interpolasi spasial—Inverse Distance Weighting (IDW), Ordinary Kriging, dan Universal Kriging—untuk eksplorasi pola risiko. Analisis mengonfirmasi klasterisasi spasial yang signifikan pada kasus campak. Model SAC terpilih sebagai model terbaik berdasarkan AIC terendah dan residual yang bebas autokorelasi spasial. GWR mengungkap heterogenitas spasial dalam efek faktor sosial-lingkungan terhadap kejadian campak. Peta interpolasi menunjukkan pola risiko yang konsisten dengan hasil klasterisasi dan pemodelan regresi. Studi ini menggarisbawahi pentingnya integrasi perspektif spasial dalam analisis penyakit menular untuk mendukung kebijakan kesehatan berbasis bukti dan spesifik lokasi. Kata kunci: analisis spasial; campak; regresi spasial; Geographically Weighted Regression; autokorelasi spasial; Moran’s I; LISA; kriging; Akaike Information Criterion

BAB I

PENDAHULUAN

Latar Belakang

Penyakit campak adalah penyakit infeksi akut yang disebabkan oleh virus campak (Measles virus) dan ditularkan melalui droplet udara [1]. Meskipun program imunisasi telah lama diterapkan, kasus campak masih ditemukan di berbagai wilayah Indonesia dengan distribusi yang tidak merata secara geografis [2]. Ketimpangan ini mengindikasikan adanya faktor spasial yang memengaruhi penyebaran penyakit.

Provinsi Jawa Timur merupakan salah satu provinsi dengan jumlah penduduk besar dan heterogenitas sosial ekonomi yang tinggi [3], sehingga berpotensi menunjukkan variasi spasial kasus campak antar kabupaten/kota. Analisis statistik konvensional sering kali mengabaikan ketergantungan spasial antar wilayah, yang dapat menyebabkan estimasi parameter yang bias dan interpretasi yang kurang tepat [4].

Oleh karena itu, penelitian ini menggunakan pendekatan Spatial Statistics untuk mengeksplorasi pola spasial prevalensi campak, menguji adanya autokorelasi spasial, serta membandingkan berbagai model spatial econometrics dan GWR guna memperoleh model yang paling sesuai secara statistik dan substantif.

Identifikasi Masalah

Berdasarkan latar belakang tersebut, tersusun beberapa pertanyaan berikut:

Apakah prevalensi campak di Provinsi Jawa Timur menunjukkan pola spasial tertentu?
Apakah terdapat ketergantungan spasial antar kabupaten/kota dalam prevalensi campak?
Model regresi spasial apa yang paling sesuai untuk menjelaskan variasi prevalensi campak antar wilayah?
Bagaimana interpretasi hasil model spasial dalam konteks kesehatan masyarakat?

Tujuan Penelitian

Tujuan penelitian ini adalah:

Memetakan prevalensi campak secara spasial pada tingkat kabupaten/kota.
Menguji keberadaan autokorelasi spasial global dan lokal pada prevalensi campak.
Mengestimasi dan membandingkan beberapa model regresi spasial.
Menginterpretasikan hasil pemodelan spasial untuk mendukung pengambilan kebijakan kesehatan masyarakat.

Batasan Penelitian

Batasan dalam penelitian ini meliputi:

Data yang digunakan bersifat agregat pada tingkat kabupaten/kota.
Analisis dilakukan untuk satu tahun pengamatan.
Metode interpolasi spasial bersifat eksploratif dan tidak digunakan untuk inferensi kausal.
Penelitian tidak membahas aspek klinis pada tingkat individu.

BAB II

TINJAUAN PUSTAKA

**Ketergantungan Spasial (Spatial Dependence)**

Ketergantungan spasial merujuk pada kondisi di mana nilai suatu variabel pada suatu lokasi dipengaruhi oleh nilai variabel yang sama di lokasi sekitarnya. Konsep ini sejalan dengan Tobler’s First Law of Geography, yang menyatakan bahwa “segala sesuatu saling berhubungan, tetapi hal yang berdekatan memiliki hubungan yang lebih kuat” [5]. Dalam konteks epidemiologi, ketergantungan spasial muncul akibat mobilitas penduduk, kesamaan lingkungan, dan akses fasilitas kesehatan yang berdekatan.

Autokorelasi Spasial

Autokorelasi spasial menggambarkan hubungan antara nilai variabel pada suatu lokasi dengan nilai di lokasi lain berdasarkan kedekatan geografis [6]. Autokorelasi spasial dapat bersifat:

Positif, jika wilayah dengan nilai tinggi (atau rendah) cenderung berdekatan [6].
Negatif, jika wilayah dengan nilai tinggi berdekatan dengan wilayah bernilai rendah [6].

Pengukuran autokorelasi spasial global umumnya dilakukan menggunakan Moran’s I dan Geary’s C, sedangkan pola lokal diidentifikasi menggunakan Local Indicators of Spatial Association (LISA) dan Getis–Ord Gi*.

Model Ekonometrika Spasial

Model ekonometrika spasial digunakan untuk menangkap pengaruh ketergantungan spasial secara global, dengan asumsi bahwa hubungan antara variabel bersifat stasioner di seluruh wilayah penelitian [7]. Model-model yang digunakan dalam penelitian ini meliputi:

Ordinary Least Squares (OLS) sebagai model dasar tanpa mempertimbangkan efek spasial.
Spatial Lag of X (SLX) untuk menangkap pengaruh variabel independen dari wilayah sekitar.
Spatial Autoregressive Model (SAR) yang memasukkan pengaruh variabel dependen dari wilayah tetangga.
Spatial Error Model (SEM) yang mengakomodasi ketergantungan spasial pada komponen galat.
Spatial Durbin Model (SDM) dan Spatial Durbin Error Model (SDEM) sebagai model gabungan antara lag dependen dan lag independen.
Spatial Autocorrelation Model (SAC) sebagai bentuk paling umum dari model regresi spasial global.

Model-model tersebut digunakan untuk menjelaskan variasi prevalensi campak dengan mempertimbangkan ketergantungan spasial antar kabupaten/kota secara global.

Geographically Weighted Regression (GWR)

Geographically Weighted Regression (GWR) merupakan metode regresi spasial lokal yang digunakan untuk menangkap ketidakstasioneran spasial, yaitu kondisi ketika hubungan antara variabel dependen dan variabel independen berbeda antar lokasi [8]. Berbeda dengan regresi global yang mengasumsikan koefisien konstan, GWR memungkinkan setiap koefisien bervariasi sesuai lokasi geografis [9].

Pendekatan GWR menggunakan pembobotan spasial berbasis jarak, di mana observasi yang lebih dekat dengan lokasi analisis memiliki pengaruh yang lebih besar dibandingkan observasi yang lebih jauh [10]. Metode ini bersifat eksploratif dan digunakan untuk mengidentifikasi variasi lokal hubungan antar variabel.

Interpolasi Spasial

Interpolasi spasial merupakan metode statistik yang digunakan untuk memperkirakan nilai suatu variabel pada lokasi yang tidak teramati berdasarkan nilai pada lokasi yang teramati dengan mempertimbangkan kedekatan spasial. Dalam penelitian epidemiologi spasial, interpolasi digunakan untuk menggambarkan permukaan risiko penyakit secara kontinu sehingga pola sebaran spasial dapat dianalisis secara lebih komprehensif dibandingkan data diskret wilayah administratif. Pada penelitian ini, interpolasi spasial diterapkan untuk memetakan distribusi kejadian penyakit campak di Provinsi Jawa Timur.

Inverse Distance Weighting (IDW)

Metode Inverse Distance Weighting (IDW) merupakan teknik interpolasi deterministik yang mengasumsikan bahwa pengaruh suatu titik pengamatan terhadap lokasi prediksi berbanding terbalik dengan jaraknya. Semakin dekat suatu titik pengamatan dengan lokasi yang diprediksi, semakin besar bobot kontribusinya terhadap nilai estimasi. IDW tidak mempertimbangkan struktur variabilitas spasial secara eksplisit, sehingga cocok digunakan sebagai pendekatan awal untuk eksplorasi pola spasial penyakit.

Kriging

Kriging merupakan metode interpolasi geostatistik yang bersifat stokastik dan mempertimbangkan struktur ketergantungan spasial melalui fungsi semivariogram. Berbeda dengan IDW, Kriging menghasilkan estimasi yang bersifat Best Linear Unbiased Estimator (BLUE) dan mampu memberikan ukuran ketidakpastian prediksi. Dalam epidemiologi spasial, Kriging banyak digunakan untuk memodelkan distribusi risiko penyakit yang dipengaruhi oleh korelasi spasial antar lokasi. Terdapat 2 jenis kriging, sebagai berikut:

Ordinary Kriging

Ordinary Kriging mengasumsikan bahwa rata-rata proses spasial bersifat konstan tetapi tidak diketahui. Estimasi nilai pada lokasi yang tidak teramati diperoleh sebagai kombinasi linear dari nilai pengamatan di sekitarnya dengan bobot yang ditentukan oleh semivariogram.

Universal Kriging

Universal Kriging memperluas Ordinary Kriging dengan mengakomodasi adanya tren spasial yang bersifat sistematis. Metode ini mengasumsikan bahwa nilai variabel terdiri dari komponen tren deterministik dan residual yang bersifat stokastik. Universal Kriging sesuai digunakan ketika terdapat variasi spasial yang dipengaruhi oleh faktor geografis atau lingkungan tertentu.

BAB III

METODOLOGI PENELITIAN

Sumber dan Jenis Data

Data yang digunakan dalam penelitian ini bersumber dari Badan Pusat Statistik (BPS) Provinsi Jawa Timur melalui publikasi “Jumlah Jenis Penyakit Tetanus, Campak, Diare, DBD, IMS Menurut Kabupaten/Kota di Provinsi Jawa Timur, 2022” [11]. Data berbentuk tabel statistik yang memuat jumlah kasus penyakit menular pada tingkat kabupaten/kota untuk tahun 2022, dengan variabel utama berupa jumlah kasus campak. Selain itu, data juga mencakup jumlah kasus tetanus, diare, demam berdarah dengue (DBD), dan infeksi menular seksual (IMS). Seluruh data bersifat cross-sectional dan dibedakan berdasarkan wilayah administratif kabupaten/kota, sehingga sesuai untuk analisis spasial guna mengidentifikasi pola distribusi dan ketergantungan spasial penyakit campak di Provinsi Jawa Timur.

Unit Spasial

Unit analisis dalam penelitian ini adalah seluruh kabupaten/kota di Provinsi Jawa Timur. Pemilihan unit ini didasarkan pada ketersediaan data dan relevansinya dalam perumusan kebijakan kesehatan tingkat regional.

Metode Analisis

Prevalensi

Perhitungan prevalensi dilakukan untuk menstandarkan jumlah kasus campak terhadap jumlah penduduk di setiap kabupaten/kota. Langkah ini bertujuan agar perbandingan antar wilayah menjadi lebih adil dan tidak dipengaruhi oleh perbedaan ukuran populasi.

Prevalensi campak dihitung sebagai:

\[ \text{Prevalensi}_i = \frac{C_i}{P_i} \times 100.000 \]

dengan: - $C_i$ : jumlah kasus campak di wilayah ke-$i$ - $P_i$ : jumlah penduduk di wilayah ke-$i$

Autokorelasi Spasial Global (Moran’s I)

Uji autokorelasi spasial global dilakukan menggunakan statistik Moran’s I untuk mengetahui apakah prevalensi campak menunjukkan pola spasial secara keseluruhan. Uji ini menentukan apakah distribusi penyakit bersifat acak, mengelompok, atau menyebar secara spasial.

Statistik Moran’s I dirumuskan sebagai:

\[ I = \frac{n}{\sum_{i}\sum_{j} w_{ij}} \frac{\sum_{i}\sum_{j} w_{ij}(y_i - \bar{y})(y_j - \bar{y})} {\sum_{i}(y_i - \bar{y})^2} \]

dengan: - $y_i$ : nilai prevalensi di wilayah ke-$i$ - $w_{ij}$ : bobot spasial antara wilayah $i$ dan $j$ - $n$ : jumlah wilayah

Autokorelasi Spasial Lokal (LISA)

Uji autokorelasi spasial lokal dilakukan menggunakan Local Indicators of Spatial Association (LISA) untuk mengidentifikasi pola spasial pada tingkat wilayah individual. Metode ini memungkinkan pendeteksian klaster lokal seperti wilayah dengan prevalensi tinggi yang dikelilingi wilayah serupa.

Indikator LISA untuk wilayah ke-$i$ dinyatakan sebagai:

\[ I_i = (y_i - \bar{y}) \sum_{j} w_{ij}(y_j - \bar{y}) \]

Model Regresi Linier OLS

Regresi Ordinary Least Squares (OLS) digunakan sebagai model dasar untuk menjelaskan hubungan antara prevalensi campak dan variabel penjelas tanpa mempertimbangkan efek spasial. Model ini berfungsi sebagai pembanding awal sebelum menerapkan model spasial.

Model regresi linear dasar (OLS) dirumuskan sebagai:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

dengan: - $\mathbf{y}$ : vektor variabel dependen - $\mathbf{X}$ : matriks variabel independen - $\boldsymbol{\beta}$ : vektor parameter - $\boldsymbol{\varepsilon}$ : vektor galat

Model Regresi Spasial Global

Model regresi spasial global seperti SAR, SEM, SLX, SDM, dan SAC digunakan untuk menangkap ketergantungan spasial antar wilayah yang tidak dapat dijelaskan oleh OLS. Model-model ini memungkinkan pengaruh wilayah sekitar dimasukkan secara eksplisit ke dalam pemodelan.

SLX

Model SLX memasukkan lag spasial dari variabel independen:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \]

dengan: - $\mathbf{W}$ : matriks bobot spasial - $\boldsymbol{\theta}$ : parameter lag spasial variabel independen

SAR

Model SAR dirumuskan sebagai:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

dengan: - $\rho$ : koefisien autoregresif spasial

SEM

Model SEM dinyatakan sebagai:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{u} \]

\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]

dengan: - $\lambda$ : koefisien ketergantungan galat spasial

SDM

Model SDM merupakan perluasan SAR dengan lag variabel independen:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon} \]

SDEM

Model SDEM merupakan perluasan SEM dengan lag variabel independen:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \mathbf{u} \]

\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]

SAC

Model SAC menggabungkan ketergantungan pada variabel dependen dan galat spasial:

\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{u} \]

\[ \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon} \]

Pemilihan Model terbaik berdasarkan AIC

Akaike Information Criterion (AIC) digunakan untuk membandingkan kinerja beberapa model regresi. Model dengan nilai AIC terkecil dipilih sebagai model terbaik karena memberikan keseimbangan antara kecocokan model dan kompleksitas parameter.

Geographically Weighted Regression (GWR)

GWR digunakan untuk menganalisis hubungan lokal antara prevalensi campak dan variabel penjelas dengan memperbolehkan koefisien regresi bervariasi antar wilayah. Metode ini bertujuan untuk mengidentifikasi adanya heterogenitas spasial yang tidak dapat ditangkap oleh model global.

Model Geographically Weighted Regression (GWR):

\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i)x_{ik} + \varepsilon_i \]

Interpolasi Spasial

Inverse Distance Weighting (IDW)

Rumus matematis IDW sebagai berikut:

\[ \hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i(s_0) \, Z(s_i)}{\sum_{i=1}^{n} w_i(s_0)} \]

\[ w_i(s_0) = \frac{1}{d(s_0, s_i)^p} \]

Kriging

Ordinary Kriging

Rumus matematis Ordinary Kriging :

\[ \hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i \, Z(s_i) \]

\[ \sum_{i=1}^{n} \lambda_i = 1 \]

Universal Kriging

Rumus matematis Universal Kriging :

\[ Z(s) = \mu(s) + \varepsilon(s) \]

\[ \mu(s) = \sum_{k=0}^{p} \beta_k f_k(s) \]

\[ \hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i) \]

BAB IV

HASIL DAN PEMBAHASAN

4.1 Peta Persebaran Kasus Campak di Jawa Timur 2022

# Pastikan kolom IR_100k tersedia
stopifnot("IR_100k" %in% names(dat))

ggplot(dat) +
  geom_sf(
    aes(fill = IR_100k),
    color = "grey60",
    linewidth = 0.2
  ) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.3,
    check_overlap = TRUE
  ) +
  scale_fill_viridis_c(
    option = "magma",
    na.value = "grey90",
    labels = scales::number_format(accuracy = 0.01)
  ) +
  theme_minimal() +
  labs(
    title = "Peta Distribusi IR Campak per 100.000 Penduduk",
    subtitle = "Provinsi Jawa Timur, Tahun 2022",
    fill = "IR per 100.000"
  ) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(size = 11),
    legend.position = "right"
  )

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Analisis spasial terhadap kejadian campak di Provinsi Jawa Timur tahun 2022 menunjukkan bahwa distribusi kasus tidak tersebar secara merata antar kabupaten/kota. Peta persebaran incidence rate campak memperlihatkan adanya konsentrasi kejadian yang relatif lebih tinggi pada wilayah tertentu, khususnya di Pulau Madura dan beberapa wilayah pesisir. Pola ini mengindikasikan bahwa faktor lokasi dan kedekatan geografis berperan dalam dinamika penyebaran penyakit campak.

4.2 Autokorelasi Spasial

## 
##  Moran I test under randomisation
## 
## data:  dat$IR_100k  
## weights: lw    
## 
## Moran I statistic standard deviate = 6.0635, p-value = 6.66e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.69383782       -0.02702703        0.01413388

Hasil uji autokorelasi spasial global menggunakan Moran’s I menunjukkan nilai Moran’s I sebesar 0,694 dengan p-value yang sangat kecil (p < 0,001). Nilai ini mengindikasikan adanya autokorelasi spasial positif yang kuat, yang berarti wilayah dengan tingkat kejadian campak tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki tingkat kejadian tinggi, demikian pula sebaliknya. Dengan demikian, hipotesis bahwa distribusi kejadian campak bersifat acak secara spasial dapat ditolak.

Analisis autokorelasi lokal melalui Local Indicators of Spatial Association (LISA) memberikan gambaran yang lebih rinci mengenai lokasi klaster spasial. Peta LISA menunjukkan adanya klaster High–High yang signifikan secara statistik di wilayah utara Jawa Timur. Klaster ini merepresentasikan kabupaten/kota dengan incidence rate campak tinggi yang dikelilingi oleh wilayah dengan nilai serupa. Sementara itu, sebagian besar wilayah Jawa Timur lainnya tergolong ke dalam kategori tidak signifikan, yang menunjukkan tidak adanya pola klaster lokal yang kuat di wilayah tersebut. Keberadaan klaster High–High ini mengindikasikan bahwa wilayah Utara merupakan area prioritas dalam pengendalian dan pencegahan campak.

gi <- localG(dat$IR_100k, lw, zero.policy = TRUE)
dat$Gi <- as.numeric(gi)

ggplot(dat) +
  geom_sf(aes(fill = Gi)) +
  scale_fill_viridis_c() +
  theme_minimal()

Temuan tersebut diperkuat oleh hasil Getis–Ord Gi* yang mengidentifikasi hotspot kejadian campak pada wilayah yang sama. Nilai statistik Gi* yang tinggi dan signifikan menunjukkan bahwa wilayah tersebut tidak hanya memiliki nilai kejadian tinggi, tetapi juga dikelilingi oleh wilayah dengan nilai kejadian tinggi secara konsisten. Konsistensi hasil antara Moran’s I, LISA, dan Getis–Ord Gi* memperkuat bukti adanya konsentrasi spasial kejadian campak di wilayah tertentu dan menegaskan bahwa pola yang diamati bukan merupakan hasil variasi acak semata.

4.3 Regresi Linier OLS

## 
## Call:
## lm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, data = df_mod)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9907 -0.9345 -0.5086  0.4958  7.2294 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     -0.390118   3.275428  -0.119   0.9059  
## Penduduk_Miskin  0.178558   0.091393   1.954   0.0588 .
## Sanitasi_Layak  -0.005312   0.031595  -0.168   0.8674  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.874 on 35 degrees of freedom
## Multiple R-squared:  0.163,  Adjusted R-squared:  0.1152 
## F-statistic: 3.409 on 2 and 35 DF,  p-value: 0.04441

dat$pred_ols <- as.numeric(predict(ols))

ggplot(dat) +
  geom_sf(aes(fill = pred_ols), color = "white", linewidth = 0.2) +
  geom_sf_text(aes(label = NAME_2),
               color = "white", size = 2.2, check_overlap = TRUE) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(
    title = "Peta Prediksi OLS",
    fill  = "Prediksi IR"
  )

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Analisis regresi linier global (OLS) dilakukan untuk mengevaluasi hubungan antara faktor sosial-lingkungan dengan kejadian campak. Hasil OLS menunjukkan bahwa variabel penduduk miskin memiliki koefisien positif terhadap incidence rate campak, meskipun signifikansinya berada pada tingkat marginal. Hal ini mengindikasikan bahwa wilayah dengan proporsi penduduk miskin yang lebih tinggi cenderung memiliki risiko kejadian campak yang lebih besar. Sebaliknya, variabel sanitasi layak tidak menunjukkan hubungan yang signifikan secara statistik. Nilai koefisien determinasi yang relatif rendah menunjukkan bahwa model OLS belum sepenuhnya mampu menjelaskan variasi kejadian campak, serta belum memperhitungkan ketergantungan spasial antar wilayah.

4.4 Regresi Spasial Global

plot_spatial_pred <- function(model, sfdata, title_txt) {
  sfdata$pred <- as.numeric(fitted(model))
  
  ggplot(sfdata) +
    geom_sf(aes(fill = pred), color = "white", linewidth = 0.2) +
    geom_sf_text(
      aes(label = NAME_2),
      color = "white",
      size = 2.2,
      check_overlap = TRUE
    ) +
    scale_fill_viridis_c(option = "magma") +
    theme_minimal() +
    labs(
      title = title_txt,
      fill  = "Prediksi IR"
    )
}

SLX

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                      Estimate   Std. Error  t value    Pr(>|t|) 
## (Intercept)          -5.673396   3.838957   -1.477848   0.148930
## Penduduk_Miskin       0.074047   0.090862    0.814942   0.420948
## Sanitasi_Layak        0.015202   0.034806    0.436778   0.665117
## lag.Penduduk_Miskin   0.359524   0.120419    2.985611   0.005299
## lag.Sanitasi_Layak    0.008407   0.038072    0.220819   0.826594

SAR

## 
## Call:lagsarlm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, 
##     data = df_mod, listw = lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.43790 -0.60734 -0.41821 -0.16403  5.64768 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.8705100  2.5051506  0.3475   0.7282
## Penduduk_Miskin  0.0297172  0.0695489  0.4273   0.6692
## Sanitasi_Layak  -0.0092636  0.0241090 -0.3842   0.7008
## 
## Rho: 0.54087, LR test value: 14.04, p-value: 0.00017895
## Asymptotic standard error: 0.13676
##     z-value: 3.955, p-value: 7.6548e-05
## Wald statistic: 15.642, p-value: 7.6548e-05
## 
## Log likelihood: -69.20517 for lag model
## ML residual variance (sigma squared): 2.0307, (sigma: 1.425)
## Number of observations: 38 
## Number of parameters estimated: 5 
## AIC: 148.41, (AIC for lm: 160.45)
## LM test for residual autocorrelation
## test value: 12.333, p-value: 0.00044512

SEM

## 
## Call:errorsarlm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, 
##     data = df_mod, listw = lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.49053 -0.56054 -0.36528 -0.15272  5.27696 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      3.449620   2.767269  1.2466   0.2126
## Penduduk_Miskin -0.055963   0.074876 -0.7474   0.4548
## Sanitasi_Layak  -0.024645   0.027427 -0.8986   0.3689
## 
## Lambda: 0.60327, LR test value: 14.21, p-value: 0.00016348
## Asymptotic standard error: 0.1258
##     z-value: 4.7956, p-value: 1.622e-06
## Wald statistic: 22.998, p-value: 1.622e-06
## 
## Log likelihood: -69.12015 for error model
## ML residual variance (sigma squared): 1.9656, (sigma: 1.402)
## Number of observations: 38 
## Number of parameters estimated: 5 
## AIC: 148.24, (AIC for lm: 160.45)

SDM

## 
## Call:lagsarlm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, 
##     data = df_mod, listw = lw, Durbin = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.795640 -0.637211 -0.301283  0.064816  4.591499 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -3.3406354  3.0857062 -1.0826  0.27898
## Penduduk_Miskin     -0.0191967  0.0753940 -0.2546  0.79902
## Sanitasi_Layak      -0.0027388  0.0278362 -0.0984  0.92162
## lag.Penduduk_Miskin  0.2406886  0.0977530  2.4622  0.01381
## lag.Sanitasi_Layak   0.0197396  0.0302737  0.6520  0.51438
## 
## Rho: 0.45491, LR test value: 9.6862, p-value: 0.0018566
## Asymptotic standard error: 0.14557
##     z-value: 3.125, p-value: 0.0017779
## Wald statistic: 9.7658, p-value: 0.0017779
## 
## Log likelihood: -66.48914 for mixed model
## ML residual variance (sigma squared): 1.8157, (sigma: 1.3475)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 146.98, (AIC for lm: 154.66)
## LM test for residual autocorrelation
## test value: 5.1484, p-value: 0.023268

SDEM

## 
## Call:errorsarlm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, 
##     data = df_mod, listw = lw, Durbin = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.716564 -0.658722 -0.336928 -0.091758  4.642811 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                        Estimate  Std. Error z value Pr(>|z|)
## (Intercept)         -2.92502648  5.06344848 -0.5777  0.56348
## Penduduk_Miskin      0.02940490  0.07831376  0.3755  0.70731
## Sanitasi_Layak       0.00032807  0.02931902  0.0112  0.99107
## lag.Penduduk_Miskin  0.22256542  0.10947135  2.0331  0.04204
## lag.Sanitasi_Layak   0.01237187  0.03438305  0.3598  0.71898
## 
## Lambda: 0.49246, LR test value: 7.5067, p-value: 0.006147
## Asymptotic standard error: 0.14621
##     z-value: 3.3682, p-value: 0.00075651
## Wald statistic: 11.345, p-value: 0.00075651
## 
## Log likelihood: -67.57887 for error model
## ML residual variance (sigma squared): 1.8989, (sigma: 1.378)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 149.16, (AIC for lm: 154.66)

SAC

## 
## Call:sacsarlm(formula = IR_100k ~ Penduduk_Miskin + Sanitasi_Layak, 
##     data = df_mod, listw = lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.338872 -0.388472 -0.290519 -0.099105  3.749014 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)      1.0924780  2.3157454  0.4718   0.6371
## Penduduk_Miskin -0.0147878  0.0505720 -0.2924   0.7700
## Sanitasi_Layak   0.0066207  0.0185526  0.3569   0.7212
## 
## Rho: -0.77824
## Asymptotic standard error: 0.13313
##     z-value: -5.8457, p-value: 5.0431e-09
## Lambda: 0.89189
## Asymptotic standard error: 0.05214
##     z-value: 17.106, p-value: < 2.22e-16
## 
## LR test value: 24.31, p-value: 5.2633e-06
## 
## Log likelihood: -64.0705 for sac model
## ML residual variance (sigma squared): 0.97098, (sigma: 0.98539)
## Number of observations: 38 
## Number of parameters estimated: 6 
## AIC: 140.14, (AIC for lm: 160.45)

Untuk mengakomodasi ketergantungan spasial tersebut, dilakukan estimasi menggunakan berbagai model regresi spasial global, yaitu SLX, SAR, SEM, SDM, SDEM, dan SAC. Hasil model SLX menunjukkan bahwa efek spasial dari variabel penduduk miskin (lag variabel) signifikan, yang mengindikasikan bahwa kondisi sosial di wilayah sekitar turut memengaruhi kejadian campak di suatu wilayah. Model SAR dan SEM menunjukkan adanya parameter spasial (rho dan lambda) yang signifikan, menandakan bahwa baik efek ketergantungan spasial langsung maupun korelasi error spasial berperan dalam model.

4.5 Geographically Weighted Regression (GWR)

dat_utm <- st_transform(dat, 32749)
pts <- as(st_centroid(dat_utm), "Spatial")

## Warning: st_centroid assumes attributes are constant over geometries

bw <- gwr.sel(IR_100k ~ Penduduk_Miskin + Sanitasi_Layak,
              data = df_mod,
              coords = coordinates(pts),
              adapt = TRUE)

## Adaptive q: 0.381966 CV score: 127.6561 
## Adaptive q: 0.618034 CV score: 133.643 
## Adaptive q: 0.236068 CV score: 122.4254 
## Adaptive q: 0.145898 CV score: 121.2911 
## Adaptive q: 0.1271733 CV score: 123.273 
## Adaptive q: 0.1851989 CV score: 121.4115 
## Adaptive q: 0.1609096 CV score: 120.7172 
## Adaptive q: 0.1646473 CV score: 120.7879 
## Adaptive q: 0.1596755 CV score: 120.6985 
## Adaptive q: 0.154413 CV score: 120.8321 
## Adaptive q: 0.1576654 CV score: 120.686 
## Adaptive q: 0.1575233 CV score: 120.6921 
## Adaptive q: 0.1585334 CV score: 120.6836 
## Adaptive q: 0.1589696 CV score: 120.689 
## Adaptive q: 0.1582173 CV score: 120.6799 
## Adaptive q: 0.1581515 CV score: 120.6791 
## Adaptive q: 0.1579658 CV score: 120.677 
## Adaptive q: 0.1578511 CV score: 120.6781 
## Adaptive q: 0.1580065 CV score: 120.6775 
## Adaptive q: 0.1579251 CV score: 120.6766 
## Adaptive q: 0.1579251 CV score: 120.6766

gwr_mod <- gwr(IR_100k ~ Penduduk_Miskin + Sanitasi_Layak,
               data = df_mod,
               coords = coordinates(pts),
               adapt = bw)

summary(gwr_mod)

##           Length Class                  Mode     
## SDF       38     SpatialPointsDataFrame S4       
## lhat       1     -none-                 logical  
## lm        11     -none-                 list     
## results    0     -none-                 NULL     
## bandwidth 38     -none-                 numeric  
## adapt      1     -none-                 numeric  
## hatmatrix  1     -none-                 logical  
## gweight    1     -none-                 character
## gTSS       1     -none-                 numeric  
## this.call  5     -none-                 call     
## fp.given   1     -none-                 logical  
## timings    8     -none-                 numeric

# Proyeksi ke CRS metrik
dat_utm <- st_transform(dat, 32749)

# Titik centroid
pts_sf <- st_centroid(dat_utm)
pts_sp <- as(pts_sf, "Spatial")

# Data model
df_gwr <- dat_utm %>%
  st_drop_geometry() %>%
  select(IR_100k, Penduduk_Miskin, Sanitasi_Layak)

# Bandwidth
bw <- gwr.sel(
  IR_100k ~ Penduduk_Miskin + Sanitasi_Layak,
  data   = df_gwr,
  coords = coordinates(pts_sp),
  adapt  = TRUE
)

## Adaptive q: 0.381966 CV score: 127.6561 
## Adaptive q: 0.618034 CV score: 133.643 
## Adaptive q: 0.236068 CV score: 122.4254 
## Adaptive q: 0.145898 CV score: 121.2911 
## Adaptive q: 0.1271733 CV score: 123.273 
## Adaptive q: 0.1851989 CV score: 121.4115 
## Adaptive q: 0.1609096 CV score: 120.7172 
## Adaptive q: 0.1646473 CV score: 120.7879 
## Adaptive q: 0.1596755 CV score: 120.6985 
## Adaptive q: 0.154413 CV score: 120.8321 
## Adaptive q: 0.1576654 CV score: 120.686 
## Adaptive q: 0.1575233 CV score: 120.6921 
## Adaptive q: 0.1585334 CV score: 120.6836 
## Adaptive q: 0.1589696 CV score: 120.689 
## Adaptive q: 0.1582173 CV score: 120.6799 
## Adaptive q: 0.1581515 CV score: 120.6791 
## Adaptive q: 0.1579658 CV score: 120.677 
## Adaptive q: 0.1578511 CV score: 120.6781 
## Adaptive q: 0.1580065 CV score: 120.6775 
## Adaptive q: 0.1579251 CV score: 120.6766 
## Adaptive q: 0.1579251 CV score: 120.6766

# GWR
gwr_mod <- gwr(
  IR_100k ~ Penduduk_Miskin + Sanitasi_Layak,
  data   = df_gwr,
  coords = coordinates(pts_sp),
  adapt  = bw
)

dat_utm$gwr_pred <- gwr_mod$SDF$pred

ggplot(dat_utm) +
  geom_sf(aes(fill = gwr_pred), color = "white", linewidth = 0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.2,
    check_overlap = TRUE
  ) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(
    title = "Peta Prediksi Geographically Weighted Regression (GWR)",
    fill = "Prediksi IR"
  )

dat_utm$coef_miskin <- gwr_mod$SDF$Penduduk_Miskin

ggplot(dat_utm) +
  geom_sf(aes(fill = coef_miskin), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(
    title = "Koefisien Lokal GWR — Penduduk Miskin",
    fill = "Koefisien"
  )

dat_utm$coef_sanit <- gwr_mod$SDF$Sanitasi_Layak

ggplot(dat_utm) +
  geom_sf(aes(fill = coef_sanit), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(
    title = "Koefisien Lokal GWR — Sanitasi Layak",
    fill = "Koefisien"
  )

Analisis Geographically Weighted Regression (GWR) selanjutnya dilakukan untuk mengeksplorasi variasi lokal hubungan antara variabel independen dan kejadian campak. Hasil GWR menunjukkan bahwa pengaruh penduduk miskin dan sanitasi layak tidak bersifat homogen di seluruh wilayah. Peta koefisien lokal menunjukkan bahwa pengaruh penduduk miskin terhadap kejadian campak cenderung lebih kuat di wilayah timur dan selatan Jawa Timur, sedangkan di beberapa wilayah barat pengaruhnya relatif lebih lemah. Demikian pula, pengaruh sanitasi layak menunjukkan variasi arah dan besaran koefisien antar wilayah. Temuan ini mengindikasikan adanya heterogenitas spasial, sehingga pendekatan lokal seperti GWR memberikan informasi tambahan yang tidak dapat diperoleh dari model global.

4.6 Pemilihan Model Terbaik

aic_tbl <- data.frame(
  Model = c("OLS","SLX","SAR","SEM","SDM","SDEM","SAC"),
  AIC   = c(
    AIC(ols),
    AIC(SLX),
    AIC(SAR),
    AIC(SEM),
    AIC(SDM),
    AIC(SDEM),
    AIC(SAC)
  )
)

aic_tbl <- aic_tbl[order(aic_tbl$AIC), ]
aic_tbl

best_model_name <- aic_tbl$Model[1]

best_model <- switch(
  best_model_name,
  "OLS"  = ols,
  "SLX"  = SLX,
  "SAR"  = SAR,
  "SEM"  = SEM,
  "SDM"  = SDM,
  "SDEM" = SDEM,
  "SAC"  = SAC
)

best_model_name

## [1] "SAC"

# Residual
residuals_best <- if (best_model_name == "OLS") {
  residuals(best_model)
} else {
  residuals(best_model)
}

# Moran's I residual
moran_res <- moran.test(
  residuals_best,
  lw,
  zero.policy = TRUE
)

moran_res

## 
##  Moran I test under randomisation
## 
## data:  residuals_best  
## weights: lw    
## 
## Moran I statistic standard deviate = -0.094816, p-value = 0.5378
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.03830106       -0.02702703        0.01413831

Di antara seluruh model yang diuji, model Spatial Autoregressive Combined (SAC) memiliki nilai AIC terendah, sehingga dipilih sebagai model terbaik. Model ini mampu menangkap ketergantungan spasial baik pada komponen lag maupun error secara simultan. Hasil uji Moran’s I terhadap residual model SAC menunjukkan tidak adanya autokorelasi spasial yang signifikan, yang menandakan bahwa struktur spasial dalam data telah berhasil dijelaskan oleh model. Dengan demikian, model SAC memberikan estimasi yang paling sesuai untuk menggambarkan hubungan antara faktor sosial-lingkungan dan kejadian campak di Jawa Timur.

4.7 Interpolasi Spasial

pts$y <- dat$IR_100k
vg <- variogram(y ~ 1, pts)
plot(vg)

# Proyeksi metrik
dat_utm <- st_transform(dat, 32749)

# Titik centroid
pts_sf <- st_point_on_surface(dat_utm)

## Warning: st_point_on_surface assumes attributes are constant over geometries

pts_sp <- as(pts_sf, "Spatial")
pts_sp$yvar <- dat_utm$IR_100k

# Grid
grid_sf <- st_make_grid(dat_utm, n = c(120,120), what = "centers") |>
  st_as_sf() |>
  st_intersection(dat_utm)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

grid_sp <- as(grid_sf, "Spatial")

# IDW
idw_res <- idw(
  yvar ~ 1,
  locations = pts_sp,
  newdata   = grid_sp,
  idp       = 2
)

## [inverse distance weighted interpolation]

# Data frame + koordinat (ANTI ERROR)
idw_df <- as.data.frame(idw_res)
xy     <- coordinates(idw_res)
idw_df$x <- xy[,1]
idw_df$y <- xy[,2]

ggplot(idw_df) +
  geom_raster(aes(x = x, y = y, fill = var1.pred), interpolate = TRUE) +
  geom_contour(aes(x = x, y = y, z = var1.pred),
               color = "white", bins = 10, linewidth = 0.3) +
  scale_fill_viridis_c(option = "magma") +
  coord_equal() +
  theme_minimal() +
  labs(
    title = "Interpolasi Spasial IDW (IR Campak)",
    fill = "Prediksi"
  )

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

vgm_fit <- fit.variogram(vg, vgm("Sph"))

## Warning in fit.variogram(vg, vgm("Sph")): No convergence after 200 iterations:
## try different initial values?

## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial
## values?

ok <- krige(
  yvar ~ 1,
  locations = pts_sp,
  newdata   = grid_sp,
  model     = vgm_fit
)

## [using ordinary kriging]

ok_df <- as.data.frame(ok)
xy_ok <- coordinates(ok)
ok_df$x <- xy_ok[,1]
ok_df$y <- xy_ok[,2]

ggplot(ok_df) +
  geom_raster(aes(x = x, y = y, fill = var1.pred), interpolate = TRUE) +
  scale_fill_viridis_c(option = "magma") +
  coord_equal() +
  theme_minimal() +
  labs(
    title = "Ordinary Kriging (OK)",
    fill = "Prediksi"
  )

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

coords_pts <- coordinates(pts_sp)
pts_sp$x <- coords_pts[,1]
pts_sp$y <- coords_pts[,2]

coords_grd <- coordinates(grid_sp)
grid_sp$x <- coords_grd[,1]
grid_sp$y <- coords_grd[,2]

uk <- krige(
  yvar ~ x + y,
  locations = pts_sp,
  newdata   = grid_sp,
  model     = vgm_fit
)

## [using universal kriging]

uk_df <- as.data.frame(uk)
xy_uk <- coordinates(uk)
uk_df$x <- xy_uk[,1]
uk_df$y <- xy_uk[,2]

ggplot(uk_df) +
  geom_raster(aes(x = x, y = y, fill = var1.pred), interpolate = TRUE) +
  scale_fill_viridis_c(option = "plasma") +
  coord_equal() +
  theme_minimal() +
  labs(
    title = "Universal Kriging (UK)",
    fill = "Prediksi"
  )

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

Selain analisis regresi, dilakukan pula pendekatan interpolasi spasial untuk memvisualisasikan pola risiko kejadian campak secara kontinu. Variogram empiris menunjukkan adanya struktur spasial yang jelas, yang menjadi dasar dalam penerapan metode kriging. Hasil Interpolasi Inverse Distance Weighting (IDW) menggambarkan area dengan risiko tinggi yang terkonsentrasi di sekitar wilayah dengan nilai kejadian tinggi, namun dengan pola yang relatif lebih kasar. Sementara itu, hasil Ordinary Kriging dan Universal Kriging menghasilkan permukaan prediksi yang lebih halus dan stabil, karena mempertimbangkan struktur variogram dan, pada Universal Kriging, adanya tren spasial. Meskipun metode interpolasi ini bersifat eksploratif karena menggunakan data agregat wilayah, hasilnya tetap berguna untuk mengidentifikasi area berisiko tinggi secara visual di luar batas administratif. Secara keseluruhan, hasil analisis spasial menunjukkan bahwa kejadian campak di Jawa Timur memiliki pola ketergantungan spasial yang kuat, dipengaruhi oleh faktor sosial-lingkungan, serta menunjukkan variasi pengaruh antar wilayah. Oleh karena itu, pendekatan spasial dan spasial-ekonometrik memberikan pemahaman yang lebih komprehensif dibandingkan analisis non-spasial semata dalam mengkaji distribusi dan determinan kejadian campak.

BAB V

KESIMPULAN DAN SARAN

Kesimpulan

Analisis spasial kejadian campak di Provinsi Jawa Timur tahun 2022 menunjukkan bahwa distribusi penyakit tidak bersifat acak, melainkan membentuk pola spasial yang jelas. Hasil uji autokorelasi global (Moran’s I) dan analisis lokal (LISA dan Getis–Ord Gi*) secara konsisten mengidentifikasi adanya klaster kejadian campak, khususnya di wilayah Madura. Temuan ini menegaskan bahwa faktor kedekatan geografis berperan penting dalam dinamika penyebaran campak antar kabupaten/kota.

Hasil pemodelan regresi menunjukkan bahwa model linier global (OLS) belum mampu menangkap ketergantungan spasial dalam data secara optimal. Penerapan model regresi spasial global memperlihatkan peningkatan kinerja model, dengan model Spatial Autoregressive Combined (SAC) sebagai model terbaik berdasarkan nilai AIC terendah. Model SAC mampu menjelaskan ketergantungan spasial baik pada komponen lag maupun error, serta menghasilkan residual yang tidak lagi menunjukkan autokorelasi spasial, sehingga memberikan estimasi yang lebih reliabel.

Analisis Geographically Weighted Regression (GWR) mengungkap adanya heterogenitas spasial dalam pengaruh faktor sosial-lingkungan terhadap kejadian campak. Variabel penduduk miskin dan sanitasi layak memiliki pengaruh yang berbeda antar wilayah, menunjukkan bahwa determinan campak bersifat kontekstual dan tidak seragam di seluruh Jawa Timur. Selain itu, pendekatan interpolasi spasial (IDW dan Kriging) memperkuat temuan adanya konsentrasi risiko di wilayah tertentu dan memberikan gambaran permukaan risiko yang bersifat kontinu secara geografis.

Saran

Berdasarkan temuan tersebut, disarankan agar upaya pengendalian dan pencegahan campak di Jawa Timur dilakukan secara berbasis wilayah, dengan prioritas pada daerah yang teridentifikasi sebagai klaster risiko tinggi, khususnya wilayah Madura. Pendekatan ini lebih efektif dibandingkan strategi yang bersifat seragam, karena mempertimbangkan karakteristik spasial dan keterkaitan antar wilayah.

Pemerintah daerah dan dinas kesehatan disarankan untuk mengintegrasikan analisis spasial dalam sistem surveilans penyakit menular. Informasi mengenai klaster, hotspot, dan variasi lokal determinan penyakit dapat digunakan sebagai dasar perencanaan intervensi, seperti peningkatan cakupan imunisasi, perbaikan sanitasi, dan penguatan layanan kesehatan di wilayah dengan risiko lebih tinggi.

Untuk penelitian selanjutnya, disarankan penggunaan data dengan resolusi yang lebih rinci, baik secara spasial maupun temporal, serta penambahan variabel lain yang relevan seperti kepadatan penduduk, mobilitas penduduk, dan cakupan imunisasi. Selain itu, pengembangan analisis spasial dinamis dan pendekatan Bayesian dapat memberikan pemahaman yang lebih mendalam mengenai perubahan pola risiko campak dari waktu ke waktu.

Daftar Pustaka

[1] Reynard, O., et al. (2022). Nebulized fusion inhibitory peptide protects cynomolgus macaques from measles virus infection. Nature Communications, 13, Article 33832. https://doi.org/10.1038/s41467-022-33832-6

[2] Wijaya, S. (2018). Pengaruh cakupan imunisasi campak terhadap incidence rate penyakit campak di Indonesia tahun 2016. Journal of Health Sciences, 11(2), 159–166. https://doi.org/10.33086/jhs.v11i2.108

[3] Aryanto, R. P., Nilogiri, A., & Wardoyo, A. E. (2024). Klasterisasi jumlah penduduk Provinsi Jawa Timur tahun 2021–2023 menggunakan algoritma K-Means. JISKA (Jurnal Informatika Sunan Kalijaga), 9(2), 134–146. https://doi.org/10.14421/jiska.2024.9.2.134-146

[4] Anselin, L. (1988). Spatial econometrics: Methods and models. Springer. https://link.springer.com/book/10.1007/978-94-015-7799-1

[5] Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46(2), 234–240. https://www.tandfonline.com/doi/10.2307/143141

[6] Cliff, A. D., & Ord, J. K. (1981). Spatial processes: Models & applications. Pion. https://onlinelibrary.wiley.com/doi/book/10.1002/978047052400

[7] LeSage, J. P., & Pace, R. K. (2009). Introduction to spatial econometrics. CRC Press. https://www.taylorfrancis.com/books/oa-mono/10.1201/9781420064254

[8] Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically weighted regression: The analysis of spatially varying relationships. Wiley. https://onlinelibrary.wiley.com/doi/book/10.1002/0470017331

[9] Brunsdon, C., Fotheringham, A. S., & Charlton, M. (1996). Geographically weighted regression: A method for exploring spatial nonstationarity. Geographical Analysis, 28(4), 281–298. https://onlinelibrary.wiley.com/doi/10.1111/j.1538-4632.1996.tb00936.x

[10] Charlton, M., Fotheringham, A. S., & Brunsdon, C. (2009). Geographically weighted regression. International Encyclopedia of Human Geography, 1–7. https://doi.org/10.1016/B978-008044910-4.00487-6

[11] Badan Pusat Statistik Provinsi Jawa Timur. (2022). Jumlah jenis penyakit tetanus, campak, diare, DBD, dan IMS menurut kabupaten/kota di Provinsi Jawa Timur. https://jatim.bps.go.id/id/statistics-table/1/Mjk3NSMx/jumlah-jenis-penyakit-tetanus--campak--diare--dbd--ims-menurut-kabupaten-kota-di-provinsi-jawa-timur--2022.html

Lampiran

Lampiran Script Analisis R

📊 Dashboard Interaktif Analisis Spasial Campak (Jawa Timur) 2022

Materi: Peta, Moran, Geary, LISA, Gi*, IDW, Variogram, OK, UK, Cross-Variogram (LMC)

options(shiny.maxRequestSize = 500*1024^2)

library(shiny) library(bslib) library(sf) library(dplyr) library(readxl) library(ggplot2) library(DT) library(spdep) library(gstat) library(sp) library(viridis) library(spatialreg) # SAR, SEM, SDM, SDEM, SAC library(spgwr) # GWR library(broom) library(car) library(lmtest) library(gridExtra)

============================================================

🧭 THEME

============================================================

my_theme <- bs_theme( version = 5, bootswatch = “flatly”, primary = “#1565C0” )

============================================================

🔧 Helper

============================================================

clean_name <- function(x) { tolower(gsub(“\s+|kabupaten|kota|,|\.”, ““, x)) }

safe_numeric <- function(x) suppressWarnings(as.numeric(x)) fmt3 <- function(x) round(x, 3)

============================================================

🧱 UI

============================================================

ui <- fluidPage( theme = my_theme, titlePanel( h2(“📊 Dashboard Analisis Spasial Campak — Jawa Timur (2022)”, style=“font-weight:700; color:#1565C0;”) ),

sidebarLayout( sidebarPanel( style=“background:#f8f9fa; border-radius:12px; padding:16px;”, h4(“⚙️ Input”, style=“font-weight:700; color:#1565C0;”), hr(),

  fileInput("shp_zip", "Upload shapefile GADM level-2 (.zip)",
            accept = c(".zip")),
  helpText("Zip harus berisi: .shp .shx .dbf .prj (dan file terkait)."),
  
  fileInput("excel", "Upload data Campak 2022 (.xlsx)",
            accept = c(".xlsx")),
  
  hr(),
  h4("🧪 Pengaturan Analisis", style="font-weight:700; color:#1565C0;"),
  
  selectInput(
    "outcome",
    "Variabel yang dianalisis (untuk Moran/LISA/Gi*/Interpolasi/Kriging):",
    choices = c("Kasus_Campak", "IR_100k"),
    selected = "IR_100k"
  ),
  
  sliderInput("alpha", "Alpha (signifikansi LISA):",
              min = 0.001, max = 0.10, value = 0.05, step = 0.001),
  
  selectInput(
    "wtype",
    "Bobot spasial (neighbors):",
    choices = c("Queen (default)" = "queen",
                "Rook" = "rook",
                "kNN (k=4)" = "knn4"),
    selected = "queen"
  ),
  
  hr(),
  h4("🌐 Interpolasi & Kriging", style="font-weight:700; color:#1565C0;"),
  sliderInput("grid_n", "Resolusi grid (n x n) (lebih besar = lebih lambat):",
              min = 60, max = 200, value = 120, step = 10),
  
  sliderInput("idw_p", "IDW power (p):",
              min = 1, max = 4, value = 2, step = 0.5),
  
  hr(),
  actionButton("run", "Jalankan / Refresh Analisis", class="btn btn-primary"),
  br(), br(),
  
  tags$small(
    "Catatan: kriging & interpolasi bersifat eksploratif karena data agregat wilayah (centroid kab/kota)."
  )
),

mainPanel(
  tabsetPanel(
    tabPanel("🧾 Data",
             h4("Data Gabungan (Shapefile + Excel)"),
             DTOutput("tbl_join"),
             br(),
             verbatimTextOutput("notes_data")
    ),
    
    tabPanel("🗺️ Peta Choropleth",
             fluidRow(
               column(12, plotOutput("map_choro", height = 650))
             )
    ),
    
    tabPanel("📌 Autokorelasi Global",
             h4("Moran’s I & Geary’s C"),
             DTOutput("tbl_global"),
             br(),
             plotOutput("moran_scatter", height = 450)
    ),
    
    tabPanel("🧩 LISA",
             h4("Local Moran’s I (LISA) Cluster Map"),
             plotOutput("map_lisa", height = 650),
             br(),
             DTOutput("tbl_lisa")
    ),
    
    tabPanel("🔥 Hotspot (Gi*)",
             h4("Getis–Ord Gi* (Hot/Cold Spots)"),
             plotOutput("map_gi", height = 650),
             br(),
             DTOutput("tbl_gi")
    ),
    
    tabPanel("🌈 IDW (Cepat)",
             h4("Interpolasi IDW via gstat (lebih cepat dari mapply)"),
             plotOutput("map_idw", height = 650),
             br(),
             DTOutput("tbl_idw_summary")
    ),
    
    tabPanel("📈 Variogram",
             h4("Variogram Empiris & Fit Model"),
             plotOutput("plot_variogram", height = 520),
             br(),
             DTOutput("tbl_vgmfit")
    ),
    
    tabPanel("📍 Kriging (OK/UK)",
             h4("Ordinary Kriging (OK) & Universal Kriging (UK)"),
             fluidRow(
               column(6, plotOutput("map_ok", height = 520)),
               column(6, plotOutput("map_ok_var", height = 520))
             ),
             hr(),
             fluidRow(
               column(6, plotOutput("map_uk", height = 520)),
               column(6, plotOutput("map_uk_var", height = 520))
             )
    ),
    
    tabPanel("🔁 Cross-Variogram (LMC)",
             h4("Variogram & Cross-Variogram (Campak vs Variabel Sekunder)"),
             selectInput(
               "secvar",
               "Pilih variabel sekunder (untuk cross-variogram):",
               choices = c("Penduduk_Miskin", "Sanitasi_Layak",
                           "Imunisasi_Laki", "Imunisasi_Perempuan"),
               selected = "Penduduk_Miskin"
             ),
             plotOutput("plot_crossvgm", height = 560),
             br(),
             verbatimTextOutput("notes_lmc")
    ),
    
    tabPanel("🧬 Epidemiologi",
             
             h4("Agent – Host – Environment"),
             verbatimTextOutput("notes_ahe"),
             br(),
             
             h4("Ukuran Epidemiologi Dasar"),
             DTOutput("tbl_epi_basic"),
             br(),
             
             h4("Prevalensi Campak per Kabupaten/Kota"),
             DTOutput("tbl_prevalensi_kab"),
             br(),
             
             h4("Standardized Morbidity Ratio (SMR)"),
             DTOutput("tbl_smr"),
             br(),
             
             h4("Empirical Bayes (EB) Smoothing"),
             plotOutput("plot_eb", height = 420),
             br(),
    ),
    
    tabPanel("📊 OLS",
             
             h4("Ordinary Least Squares (OLS)"),
             
             plotOutput("map_ols", height = 520),
             hr(),
             
             h4("Ringkasan Model"),
             verbatimTextOutput("summary_ols"),
             hr(),
             
             h4("Kriteria Informasi"),
             DTOutput("tbl_ols_aic"),
             hr(),
             
             h4("Uji Multikolinearitas (VIF)"),
             DTOutput("tbl_vif"),
             hr(),
             
             h4("Uji Heteroskedastisitas (Breusch–Pagan)"),
             DTOutput("tbl_bp")
    ),
    
    tabPanel("📐 Spatial Regression",
             h4("Model Econometric Spasial"),
             DTOutput("tbl_aic"),
             hr(),
             h5("SLX"),
             plotOutput("map_slx", height = 450),
             verbatimTextOutput("summary_slx"),
             hr(),
             h5("SAR"),
             plotOutput("map_sar", height = 450),
             verbatimTextOutput("summary_sar"),
             hr(),
             h5("SEM"),
             plotOutput("map_sem", height = 450),
             verbatimTextOutput("summary_sem"),
             hr(),
             h5("SDM"),
             plotOutput("map_sdm", height = 450),
             verbatimTextOutput("summary_sdm"),
             hr(),
             h5("SDEM"),
             plotOutput("map_sdem", height = 450),
             verbatimTextOutput("summary_sdem"),
             hr(),
             h5("SAC"),
             plotOutput("map_sac", height = 450),
             verbatimTextOutput("summary_sac")
    ),
    
    tabPanel("🌍 GWR",
             h4("Geographically Weighted Regression"),
             plotOutput("map_gwr_pred", height = 520),
             hr(),
             plotOutput("map_gwr_coef", height = 520),
             hr(),
             verbatimTextOutput("summary_gwr"),
             hr(),
             h4("Uji Normalitas Residual GWR"),
             DTOutput("tbl_gwr_normality"),
             br(),
             h4("Moran’s I Residual GWR"),
             DTOutput("tbl_gwr_moran"),
    )
  )
)

============================================================

⚙️ SERVER

============================================================

server <- function(input, output, session) {

# ————————— # 1) Read shapefile zip # ————————— shp_sf <- eventReactive(input$run, { req(input$shp_zip) tmpdir <- tempdir() unzip(input$shp_zip$datapath, exdir = tmpdir) shp_file <- list.files(tmpdir, “\.shp$”, full.names = TRUE)[1] shp <- st_read(shp_file, quiet = TRUE)

# filter Jawa Timur (umumnya NAME_1)
if ("NAME_1" %in% names(shp)) {
  shp <- shp %>% filter(tolower(NAME_1) == "jawa timur")
}

shp <- st_make_valid(shp)
if (is.na(st_crs(shp))) st_crs(shp) <- 4326

shp

})

# ————————— # 2) Read excel # ————————— df_xlsx <- eventReactive(input$run, { req(input$excel) df <- read_excel(input$excel$datapath)

# normalize names
names(df) <- trimws(gsub("\\s+", "_", names(df)))

# Ensure required cols exist
# Expected: Kabupaten_Kota, Tahun, Kasus_Campak, Penduduk, Sanitasi_Layak, Penduduk_Miskin, Imunisasi_Laki, Imunisasi_Perempuan
df

})

# ————————— # 3) Join shapefile + excel (only year 2022) # ————————— jatim_join <- reactive({ shp <- shp_sf() df <- df_xlsx()

req(shp, df)

# Filter year 2022 if column exists
if ("Tahun" %in% names(df)) {
  df <- df %>% filter(Tahun == 2022)
}

# Key from shapefile
if (!("NAME_2" %in% names(shp))) {
  stop("Shapefile tidak punya kolom NAME_2 (nama kab/kota).")
}

shp <- shp %>% mutate(.key = clean_name(NAME_2))

# Key from excel
if (!("Kabupaten_Kota" %in% names(df))) {
  stop("Excel harus punya kolom 'Kabupaten_Kota'.")
}

df <- df %>%
  mutate(.key = clean_name(Kabupaten_Kota)) %>%
  mutate(
    Kasus_Campak = safe_numeric(Kasus_Campak),
    Penduduk = safe_numeric(Penduduk),
    Sanitasi_Layak = safe_numeric(Sanitasi_Layak),
    Penduduk_Miskin = safe_numeric(Penduduk_Miskin),
    Imunisasi_Laki = safe_numeric(Imunisasi_Laki),
    Imunisasi_Perempuan = safe_numeric(Imunisasi_Perempuan)
  )

joined <- left_join(shp, df, by = ".key")

# compute IR per 100k
joined <- joined %>%
  mutate(
    IR_100k = ifelse(!is.na(Penduduk) & Penduduk > 0,
                     Kasus_Campak / Penduduk * 100000,
                     NA_real_)
  )

joined

})

# ————————— # 4) Neighbors / weights # ————————— make_listw <- function(sfobj, wtype = “queen”) { sfobj <- sfobj %>% filter(!st_is_empty(geometry))

if (wtype == "queen") {
  nb <- poly2nb(sfobj, queen = TRUE)
  nb2listw(nb, style = "W", zero.policy = TRUE)
} else if (wtype == "rook") {
  nb <- poly2nb(sfobj, queen = FALSE)
  nb2listw(nb, style = "W", zero.policy = TRUE)
} else {
  # kNN (k=4) from centroids
  ctr <- st_coordinates(st_centroid(sfobj))
  nb <- knn2nb(knearneigh(ctr, k = 4))
  nb2listw(nb, style = "W", zero.policy = TRUE)
}

}

# ————————— # 5) Data table # ————————— output$tbl_join <- renderDT({ req(input$run) dat <- jatim_join() %>% st_drop_geometry() datatable(dat, options = list(pageLength = 15, scrollX = TRUE)) })

output$notes_data <- renderPrint({ req(input$run) cat(“Catatan data:\n”) cat(“- Data digabung berdasarkan nama kab/kota (NAME_2 di shapefile vs Kabupaten_Kota di Excel).”) cat(“- IR_100k dihitung: Kasus_Campak / Penduduk * 100000.”) cat(“- Analisis spasial berbasis kab/kota (area), interpolasi & kriging menggunakan centroid (eksploratif).”) })

# ————————— # 6) Choropleth map (Cases / IR) # ————————— output$map_choro <- renderPlot({ req(input$run) sfdat <- jatim_join()

var <- input$outcome
if (!(var %in% names(sfdat))) {
  ggplot() + theme_void() + ggtitle(paste("Kolom tidak ditemukan:", var))
} else {
  ggplot(sfdat) +
    geom_sf(aes(fill = .data[[var]]), color = "grey60", linewidth = 0.2) +
    geom_sf_text(
      aes(label = NAME_2),
      color = "white",
      size = 2.2,
      check_overlap = TRUE
    ) +
    scale_fill_viridis_c(option = "magma", na.value = "grey95") +
    theme_minimal() +
    labs(
      title = "Peta Choropleth Campak — Jawa Timur (2022)",
      subtitle = paste("Variabel:", var),
      fill = var
    )
}

})

# ————————— # 7) Global Moran & Geary + Moran scatter # ————————— output$tbl_global <- renderDT({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(datatable(data.frame(Message = "Data valid terlalu sedikit untuk uji autokorelasi."), options = list(dom='t')))
}

lw <- make_listw(sf_use, input$wtype)

mor <- moran.test(sf_use[[var]], lw, zero.policy = TRUE)
gea <- geary.test(sf_use[[var]], lw, zero.policy = TRUE)

out <- data.frame(
  Statistik = c("Moran's I", "Moran p-value", "Geary's C", "Geary p-value"),
  Nilai = c(as.numeric(mor$estimate[1]),
            as.numeric(mor$p.value),
            as.numeric(gea$estimate[1]),
            as.numeric(gea$p.value))
)

datatable(out, options = list(dom='t'))

})

output$moran_scatter <- renderPlot({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) return(NULL)

lw <- make_listw(sf_use, input$wtype)

x <- sf_use[[var]]
zx <- as.numeric(scale(x))
lagx <- lag.listw(lw, zx, zero.policy = TRUE)

dfp <- data.frame(z = zx, lagz = lagx)

ggplot(dfp, aes(z, lagz)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(
    title = "Moran Scatter Plot",
    subtitle = paste("Variabel:", var),
    x = "Z-score",
    y = "Spatial lag (W*Z)"
  )

})

# ————————— # 8) LISA map + table # ————————— output$map_lisa <- renderPlot({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(ggplot() + theme_void() + ggtitle("Data valid terlalu sedikit untuk LISA."))
}

lw <- make_listw(sf_use, input$wtype)

x <- sf_use[[var]]
zx <- as.numeric(scale(x))
lagz <- lag.listw(lw, zx, zero.policy = TRUE)

lisa <- localmoran(zx, lw, zero.policy = TRUE)
pval <- lisa[,5]

alpha <- input$alpha
quad <- case_when(
  zx >= 0 & lagz >= 0 ~ "High-High",
  zx <  0 & lagz <  0 ~ "Low-Low",
  zx >= 0 & lagz <  0 ~ "High-Low (Outlier)",
  zx <  0 & lagz >= 0 ~ "Low-High (Outlier)"
)

sf_use$LISA <- ifelse(pval <= alpha, quad, "Not significant")
sf_use$LISA <- factor(sf_use$LISA,
                      levels = c("High-High","Low-Low","High-Low (Outlier)","Low-High (Outlier)","Not significant"))

ggplot(sf_use) +
  geom_sf(aes(fill = LISA), color="white", linewidth=0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.2,
    check_overlap = TRUE
  ) +
  scale_fill_manual(values = c(
    "High-High"="#d73027",
    "Low-Low"="#4575b4",
    "High-Low (Outlier)"="#fdae61",
    "Low-High (Outlier)"="#74add1",
    "Not significant"="grey85"
  )) +
  theme_minimal() +
  labs(title = "LISA Cluster Map (Local Moran’s I)",
       subtitle = paste("Alpha =", alpha, "| Variabel:", var),
       fill = NULL)

})

output$tbl_lisa <- renderDT({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(datatable(data.frame(Message="Data valid terlalu sedikit."), options = list(dom='t')))
}

lw <- make_listw(sf_use, input$wtype)

x <- sf_use[[var]]
zx <- as.numeric(scale(x))
lagz <- lag.listw(lw, zx, zero.policy = TRUE)

lisa <- localmoran(zx, lw, zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","P")

alpha <- input$alpha
quad <- case_when(
  zx >= 0 & lagz >= 0 ~ "High-High",
  zx <  0 & lagz <  0 ~ "Low-Low",
  zx >= 0 & lagz <  0 ~ "High-Low (Outlier)",
  zx <  0 & lagz >= 0 ~ "Low-High (Outlier)"
)

out <- sf_use %>%
  st_drop_geometry() %>%
  mutate(z = zx, lagz = lagz) %>%
  bind_cols(lisa_df) %>%
  mutate(cluster = ifelse(P <= alpha, quad, "Not significant")) %>%
  select(NAME_2, all_of(var), Ii, Zi, P, cluster)

datatable(out, options = list(pageLength=10, scrollX=TRUE))

})

# ————————— # 9) Getis-Ord Gi* hotspot # ————————— output$map_gi <- renderPlot({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(ggplot() + theme_void() + ggtitle("Data valid terlalu sedikit untuk Gi*."))
}

lw <- make_listw(sf_use, input$wtype)

gi <- localG(sf_use[[var]], lw, zero.policy = TRUE)
zgi <- as.numeric(gi)

sf_use$hotcold <- case_when(
  zgi >= 1.96 ~ "Hot spot (p≈0.05)",
  zgi <= -1.96 ~ "Cold spot (p≈0.05)",
  TRUE ~ "Not significant"
)

ggplot(sf_use) +
  geom_sf(aes(fill = hotcold), color="white", linewidth=0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.2,
    check_overlap = TRUE
  ) + 
  scale_fill_manual(values = c(
    "Hot spot (p≈0.05)"="#b2182b",
    "Cold spot (p≈0.05)"="#2166ac",
    "Not significant"="grey85"
  )) +
  theme_minimal() +
  labs(title = "Getis–Ord Gi* Hot/Cold Spots",
       subtitle = paste("Variabel:", var),
       fill = NULL)

})

output$tbl_gi <- renderDT({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(datatable(data.frame(Message="Data valid terlalu sedikit."), options = list(dom='t')))
}

lw <- make_listw(sf_use, input$wtype)
gi <- localG(sf_use[[var]], lw, zero.policy = TRUE)
zgi <- as.numeric(gi)

out <- sf_use %>%
  st_drop_geometry() %>%
  mutate(z_gi = zgi) %>%
  mutate(kategori = case_when(
    z_gi >= 1.96 ~ "Hot spot",
    z_gi <= -1.96 ~ "Cold spot",
    TRUE ~ "Not significant"
  )) %>%
  select(NAME_2, all_of(var), z_gi, kategori)

datatable(out, options = list(pageLength=10, scrollX=TRUE))

})

# ————————— # 10) IDW via gstat (FINAL FIX) # ————————— output$map_idw <- renderPlot({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% dplyr::filter(!is.na(.data[[var]]))
if (nrow(sf_use) < 5) {
  return(ggplot() + theme_void() +
           ggtitle("Data valid terlalu sedikit untuk IDW."))
}

# =====================================================
# IMPORTANT: Project ke CRS METRIK (WAJIB untuk IDW)
# =====================================================
sf_use <- st_transform(sf_use, 32749)  # UTM 49S (Jawa Timur)

# centroid points
pts_sf <- st_point_on_surface(sf_use)
pts_sp <- as(pts_sf, "Spatial")
pts_sp$yvar <- sf_use[[var]]

# grid (centers)
n <- input$grid_n
grid_sf <- st_make_grid(sf_use, n = c(n, n), what = "centers") |>
  st_as_sf() |>
  st_intersection(sf_use)

grid_sp <- as(grid_sf, "Spatial")

# IDW
idw_res <- gstat::idw(
  yvar ~ 1,
  locations = pts_sp,
  newdata   = grid_sp,
  idp       = input$idw_p
)

# =====================================================
# AMBIL KOORDINAT SECARA EKSPLISIT (ANTI x1/x2 ERROR)
# =====================================================
idw_df <- as.data.frame(idw_res)
xy_idw <- sp::coordinates(idw_res)

idw_df$x <- xy_idw[, 1]
idw_df$y <- xy_idw[, 2]

df_obs <- data.frame(
  x = sp::coordinates(pts_sp)[, 1],
  y = sp::coordinates(pts_sp)[, 2]
)

bb <- st_bbox(sf_use)

ggplot() +
  geom_raster(
    data = idw_df,
    aes(x = x, y = y, fill = var1.pred),
    interpolate = TRUE
  ) +
  geom_contour(
    data = idw_df,
    aes(x = x, y = y, z = var1.pred),
    color = "white",
    linewidth = 0.35,
    bins = 10
  ) +
  geom_point(
    data = df_obs,
    aes(x = x, y = y),
    shape = 21,
    size = 2.5,
    stroke = 0.8,
    fill = "yellow",
    color = "black"
  ) +
  coord_equal(
    xlim = c(bb["xmin"], bb["xmax"]),
    ylim = c(bb["ymin"], bb["ymax"])
  ) + 
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(
    title = "Interpolasi IDW (gstat)",
    subtitle = paste("Variabel:", var, "| power (p) =", input$idw_p),
    fill = "Prediksi"
  )

})

# ————————— # 11) Variogram & fit (OK base) # ————————— vgm_fit_obj <- reactive({ req(input$run) sfdat <- jatim_join() var <- input$outcome

sf_use <- sfdat %>% filter(!is.na(.data[[var]]))
validate(need(nrow(sf_use) >= 8, "Data valid terlalu sedikit untuk variogram/kriging."))

pts_sf <- st_point_on_surface(sf_use)
pts_sp <- as(pts_sf, "Spatial")
pts_sp$yvar <- sf_use[[var]]

v <- variogram(yvar ~ 1, pts_sp)

# initial guess (robust-ish)
vfit <- fit.variogram(
  v,
  model = vgm(
    psill = var(pts_sp$yvar, na.rm = TRUE),
    model = "Sph",
    range = max(v$dist, na.rm = TRUE) / 2,
    nugget = 0
  )
)

list(v = v, vfit = vfit, pts_sp = pts_sp, sf_use = sf_use)

})

output$plot_variogram <- renderPlot({ req(input$run) obj <- vgm_fit_obj() plot(obj$v, obj$vfit, main = “Variogram Empiris & Model Fit (Spherical)”) })

output$tbl_vgmfit <- renderDT({ req(input$run) obj <- vgm_fit_obj() datatable(as.data.frame(obj$vfit), options = list(dom = ‘t’)) })

# ————————— # 12) OK & UK Kriging maps – FIX: UTM + aman dari x1/x2 # ————————— make_kriging_grid <- function(sf_use, n) { # sf_use sudah UTM grid_sf <- st_make_grid(sf_use, n = c(n, n), what = “centers”) |> st_as_sf() |> st_intersection(sf_use)

as(grid_sf, "Spatial")

}

output$map_ok <- renderPlot({ req(input$run) obj <- vgm_fit_obj()

n <- input$grid_n
grid_sp <- make_kriging_grid(obj$sf_use, n)

ok <- gstat::krige(yvar ~ 1, locations = obj$pts_sp, newdata = grid_sp, model = obj$vfit)

ok_df <- as.data.frame(ok)
xy_ok <- sp::coordinates(ok)
ok_df$x <- xy_ok[, 1]
ok_df$y <- xy_ok[, 2]

bb <- st_bbox(obj$sf_use)

df_obs <- data.frame(
  x = sp::coordinates(obj$pts_sp)[, 1],
  y = sp::coordinates(obj$pts_sp)[, 2]
)

contour_levels <- pretty(range(ok_df$var1.pred, na.rm = TRUE), n = 10)

ggplot() +
  geom_raster(data = ok_df, aes(x = x, y = y, fill = var1.pred), interpolate = TRUE) +
  geom_contour(data = ok_df, aes(x = x, y = y, z = var1.pred),
               breaks = contour_levels, color = "white", alpha = 0.9, linewidth = 0.35) +
  geom_point(data = df_obs, aes(x = x, y = y),
             shape = 21, size = 2.8, stroke = 0.9, color = "black", fill = "yellow") +
  coord_equal(xlim = c(bb["xmin"], bb["xmax"]), ylim = c(bb["ymin"], bb["ymax"])) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(title = "Ordinary Kriging (OK) — Prediksi",
       subtitle = paste("Variabel:", input$outcome),
       fill = "Pred")

})

output$map_ok_var <- renderPlot({ req(input$run) obj <- vgm_fit_obj()

n <- input$grid_n
grid_sp <- make_kriging_grid(obj$sf_use, n)

ok <- gstat::krige(yvar ~ 1, locations = obj$pts_sp, newdata = grid_sp, model = obj$vfit)

ok_df <- as.data.frame(ok)
xy_ok <- sp::coordinates(ok)
ok_df$x <- xy_ok[, 1]
ok_df$y <- xy_ok[, 2]

bb <- st_bbox(obj$sf_use)

df_obs <- data.frame(
  x = sp::coordinates(obj$pts_sp)[, 1],
  y = sp::coordinates(obj$pts_sp)[, 2]
)

contour_levels <- pretty(range(ok_df$var1.var, na.rm = TRUE), n = 10)

ggplot() +
  geom_raster(data = ok_df, aes(x = x, y = y, fill = var1.var), interpolate = TRUE) +
  geom_contour(data = ok_df, aes(x = x, y = y, z = var1.var),
               breaks = contour_levels, color = "white", alpha = 0.9, linewidth = 0.35) +
  geom_point(data = df_obs, aes(x = x, y = y),
             shape = 21, size = 2.8, stroke = 0.9, color = "black", fill = "yellow") +
  coord_equal(xlim = c(bb["xmin"], bb["xmax"]), ylim = c(bb["ymin"], bb["ymax"])) +
  scale_fill_viridis_c(option = "viridis") +
  theme_minimal() +
  labs(title = "Ordinary Kriging (OK) — Varians Prediksi",
       subtitle = "Semakin besar = semakin tidak pasti",
       fill = "Var")

})

output$map_uk <- renderPlot({ req(input$run) obj <- vgm_fit_obj()

# drift x+y untuk UK
coords <- sp::coordinates(obj$pts_sp)
obj$pts_sp$x <- coords[, 1]
obj$pts_sp$y <- coords[, 2]

n <- input$grid_n
grid_sp <- make_kriging_grid(obj$sf_use, n)
coords_g <- sp::coordinates(grid_sp)
grid_sp$x <- coords_g[, 1]
grid_sp$y <- coords_g[, 2]

uk <- gstat::krige(yvar ~ x + y, locations = obj$pts_sp, newdata = grid_sp, model = obj$vfit)

uk_df <- as.data.frame(uk)
xy_uk <- sp::coordinates(uk)
uk_df$x <- xy_uk[, 1]
uk_df$y <- xy_uk[, 2]

bb <- st_bbox(obj$sf_use)

df_obs <- data.frame(
  x = sp::coordinates(obj$pts_sp)[, 1],
  y = sp::coordinates(obj$pts_sp)[, 2]
)

contour_levels <- pretty(range(uk_df$var1.pred, na.rm = TRUE), n = 10)

ggplot() +
  geom_raster(data = uk_df, aes(x = x, y = y, fill = var1.pred), interpolate = TRUE) +
  geom_contour(data = uk_df, aes(x = x, y = y, z = var1.pred),
               breaks = contour_levels, color = "white", alpha = 0.9, linewidth = 0.35) +
  geom_point(data = df_obs, aes(x = x, y = y),
             shape = 21, size = 2.8, stroke = 0.9, color = "black", fill = "yellow") +
  coord_equal(xlim = c(bb["xmin"], bb["xmax"]), ylim = c(bb["ymin"], bb["ymax"])) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(title = "Universal Kriging (UK) — Prediksi",
       subtitle = "Drift: x + y",
       fill = "Pred")

})

output$map_uk_var <- renderPlot({ req(input$run) obj <- vgm_fit_obj()

coords <- sp::coordinates(obj$pts_sp)
obj$pts_sp$x <- coords[, 1]
obj$pts_sp$y <- coords[, 2]

n <- input$grid_n
grid_sp <- make_kriging_grid(obj$sf_use, n)
coords_g <- sp::coordinates(grid_sp)
grid_sp$x <- coords_g[, 1]
grid_sp$y <- coords_g[, 2]

uk <- gstat::krige(yvar ~ x + y, locations = obj$pts_sp, newdata = grid_sp, model = obj$vfit)

uk_df <- as.data.frame(uk)
xy_uk <- sp::coordinates(uk)
uk_df$x <- xy_uk[, 1]
uk_df$y <- xy_uk[, 2]

bb <- st_bbox(obj$sf_use)

df_obs <- data.frame(
  x = sp::coordinates(obj$pts_sp)[, 1],
  y = sp::coordinates(obj$pts_sp)[, 2]
)

contour_levels <- pretty(range(uk_df$var1.var, na.rm = TRUE), n = 10)

ggplot() +
  geom_raster(data = uk_df, aes(x = x, y = y, fill = var1.var), interpolate = TRUE) +
  geom_contour(data = uk_df, aes(x = x, y = y, z = var1.var),
               breaks = contour_levels, color = "white", alpha = 0.9, linewidth = 0.35) +
  geom_point(data = df_obs, aes(x = x, y = y),
             shape = 21, size = 2.8, stroke = 0.9, color = "black", fill = "yellow") +
  coord_equal(xlim = c(bb["xmin"], bb["xmax"]), ylim = c(bb["ymin"], bb["ymax"])) +
  scale_fill_viridis_c(option = "viridis") +
  theme_minimal() +
  labs(title = "Universal Kriging (UK) — Varians Prediksi",
       subtitle = "Ketidakpastian model UK",
       fill = "Var")

})

# ————————— # 13) Cross-Variogram / LMC (Exploratory) # ————————— output$plot_crossvgm <- renderPlot({ req(input$run) sfdat <- jatim_join()

# use only non-NA
sec <- input$secvar
mainv <- input$outcome

sf_use <- sfdat %>%
  filter(!is.na(.data[[mainv]]), !is.na(.data[[sec]]))

validate(need(nrow(sf_use) >= 8, "Data valid terlalu sedikit untuk cross-variogram."))

pts_sf <- st_point_on_surface(sf_use)
pts_sp <- as(pts_sf, "Spatial")

pts_sp$campak <- sf_use[[mainv]]
pts_sp$secvar <- sf_use[[sec]]

gs <- gstat(NULL, id = "campak", formula = campak ~ 1, data = pts_sp)
gs <- gstat(gs,   id = "secvar", formula = secvar ~ 1, data = pts_sp)

v.cross <- variogram(gs)
plot(v.cross, main = paste("Variogram & Cross-Variogram:", mainv, "vs", sec))

})

output$notes_lmc <- renderPrint({ cat(“Catatan LMC/Cokriging (eksploratif):”) cat(“- Cross-variogram mengevaluasi keterkaitan spasial antara campak dan variabel sekunder.”) cat(“- Untuk UAS, cukup tampilkan plot variogram & cross-variogram sebagai bukti konsep.”) cat(“- Prediksi cokriging/LMC penuh dapat ditambahkan bila dibutuhkan, namun tidak wajib.”) })

# ============================================================ # 14) EPIDEMIOLOGI (NON-SPASIAL) # ============================================================ # ————————— # Agent – Host – Environment # ————————— output$notes_ahe <- renderPrint({ cat(“Agent:”) cat(“- Virus campak (Measles virus, genus Morbillivirus)”) cat(“- Sangat menular melalui droplet udara”)

cat("Host:\n")
cat("- Anak-anak tidak imun\n")
cat("- Cakupan imunisasi rendah (Imunisasi_Laki / Imunisasi_Perempuan)\n")
cat("- Status gizi dan kepadatan penduduk\n\n")

cat("Environment:\n")
cat("- Sanitasi lingkungan (Sanitasi_Layak)\n")
cat("- Kemiskinan (Penduduk_Miskin)\n")
cat("- Akses layanan kesehatan\n")

})

# tabel utama output$tbl_epi_basic <- renderDT({ req(input$run) dat <- jatim_join() %>% st_drop_geometry()

out <- dat %>%
  summarise(
    Total_Kasus = sum(Kasus_Campak, na.rm = TRUE),
    Total_Penduduk = sum(Penduduk, na.rm = TRUE),
    IR_per_100k = (Total_Kasus / Total_Penduduk) * 100000
  )

datatable(
  fmt3(out),
  options = list(dom = "t")
)

})

# ————————— # Prevalensi per Kabupaten/Kota # ————————— output$tbl_prevalensi_kab <- renderDT({ req(input$run)

dat <- jatim_join() %>% 
  st_drop_geometry() %>%
  filter(
    !is.na(Kasus_Campak),
    !is.na(Penduduk),
    Penduduk > 0
  ) %>%
  mutate(
    Prevalensi = Kasus_Campak / Penduduk
  ) %>%
  select(
    Kabupaten_Kota = NAME_2,
    Kasus_Campak,
    Penduduk,
    Prevalensi
  ) %>%
  mutate(
    across(where(is.numeric), ~ round(.x, 8))
  )

datatable(
  dat,
  options = list(
    pageLength = 10,
    scrollX = TRUE
  )
)

})

# ————————— # SMR + CI 95% (Poisson) # ————————— output$tbl_smr <- renderDT({ req(input$run) dat <- jatim_join() %>% st_drop_geometry()

validate(need(all(c("Kasus_Campak","Penduduk") %in% names(dat)),
              "Variabel kasus & penduduk diperlukan"))

rbar <- sum(dat$Kasus_Campak, na.rm = TRUE) /
  sum(dat$Penduduk, na.rm = TRUE)

dat <- dat %>%
  mutate(
    E = rbar * Penduduk,
    O = Kasus_Campak,
    SMR = O / E,
    LCL = ifelse(O > 0, qchisq(0.025, 2*O)/(2*E), 0),
    UCL = qchisq(0.975, 2*(O+1))/(2*E)
  )

datatable(
  dat %>% select(NAME_2, O, E, SMR, LCL, UCL),
  options = list(pageLength = 10, scrollX = TRUE)
)

})

# ————————— # Empirical Bayes (EB) # ————————— output$plot_eb <- renderPlot({ req(input$run) dat <- jatim_join() %>% st_drop_geometry()

rbar <- sum(dat$Kasus_Campak, na.rm = TRUE) /
  sum(dat$Penduduk, na.rm = TRUE)

dat <- dat %>%
  mutate(
    E = rbar * Penduduk,
    O = Kasus_Campak,
    SMR = O / E,
    EB = (1 + O) / (1 + E)
  )

ggplot(dat, aes(SMR, EB)) +
  geom_point(color = "#1565C0", size = 2) +
  geom_abline(linetype = 2, color = "grey40") +
  theme_minimal() +
  labs(
    title = "Empirical Bayes Smoothing",
    subtitle = "EB mengecilkan ekstrem SMR pada area berpopulasi kecil",
    x = "SMR (Raw)",
    y = "EB Estimate"
  )

})

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

ols_model <- reactive({ req(input$run)

sfdat <- jatim_join()
var   <- input$outcome

sf_use <- sfdat %>%
  filter(
    !is.na(.data[[var]]),
    !is.na(Penduduk_Miskin),
    !is.na(Sanitasi_Layak)
  )

validate(need(nrow(sf_use) >= 10, "Data terlalu sedikit untuk OLS"))

df <- sf_use %>%
  st_drop_geometry() %>%
  select(all_of(var), Penduduk_Miskin, Sanitasi_Layak)

lm(
  as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
  data = df
)

})

output$map_ols <- renderPlot({ model <- ols_model() sf_use <- jatim_join()

sf_use <- sf_use %>%
  filter(
    !is.na(.data[[input$outcome]]),
    !is.na(Penduduk_Miskin),
    !is.na(Sanitasi_Layak)
  )

sf_use$pred <- round(as.numeric(predict(model)), 3)

ggplot(sf_use) +
  geom_sf(aes(fill = pred), color = "white", linewidth = 0.2) +
  geom_sf_text(aes(label = NAME_2), color = "white", size = 2) +
  scale_fill_viridis_c(option = "magma", labels = scales::number_format(accuracy = 0.001)) +
  theme_minimal() +
  labs(title = "Prediksi OLS", fill = "Prevalensi")

})

output$summary_ols <- renderPrint({ summary(ols_model()) })

output$tbl_ols_aic <- renderDT({ model <- ols_model()

datatable(
  data.frame(
    Model = "OLS",
    AIC = round(AIC(model), 3)
  ),
  options = list(dom = "t")
)

})

output$tbl_ols_aic <- renderDT({ model <- ols_model()

datatable(
  data.frame(
    Model = "OLS",
    AIC = round(AIC(model), 3)
  ),
  options = list(dom = "t")
)

})

output$tbl_vif <- renderDT({ model <- ols_model()

vif_val <- car::vif(model)

datatable(
  data.frame(
    Variabel = names(vif_val),
    VIF = round(as.numeric(vif_val), 3)
  ),
  options = list(dom = "t")
)

})

output$tbl_bp <- renderDT({ model <- ols_model()

bp <- bptest(model)

datatable(
  data.frame(
    Statistik = c("BP statistic", "df", "P-value"),
    Nilai = round(c(bp$statistic, bp$parameter, bp$p.value), 3)
  ),
  options = list(dom = "t")
)

})

# ============================================================ # 14) Spatial Econometric Models # ============================================================ spatial_models <- reactive({ req(input$run)

sfdat <- jatim_join()
var   <- input$outcome

sf_use <- sfdat %>%
  filter(
    !is.na(.data[[var]]),
    !is.na(Penduduk_Miskin),
    !is.na(Sanitasi_Layak)
  )

validate(need(nrow(sf_use) >= 10, "Data terlalu sedikit untuk model spasial"))

lw <- make_listw(sf_use, input$wtype)

df <- sf_use %>%
  st_drop_geometry() %>%
  select(all_of(var), Penduduk_Miskin, Sanitasi_Layak)

list(
  sf = sf_use,
  
  SLX  = lmSLX(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw
  ),
  
  SAR  = lagsarlm(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw
  ),
  
  SEM  = errorsarlm(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw
  ),
  
  SDM  = lagsarlm(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw,
    Durbin = TRUE
  ),
  
  SDEM = errorsarlm(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw,
    Durbin = TRUE
  ),
  
  SAC  = sacsarlm(
    as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
    data = df,
    listw = lw
  )
)

})

output$tbl_aic <- renderDT({ obj <- spatial_models()

aic_df <- data.frame(
  Model = c("SLX","SAR","SEM","SDM","SDEM","SAC"),
  AIC = fmt3(c(
    AIC(obj$SLX),
    AIC(obj$SAR),
    AIC(obj$SEM),
    AIC(obj$SDM),
    AIC(obj$SDEM),
    AIC(obj$SAC)
  ))
)

datatable(aic_df %>% arrange(AIC), options = list(dom = "t"))

})

# OUTPUT SEMUA MODEL output$summary_slx <- renderPrint({ summary(spatial_models()$SLX) }) output$summary_sar <- renderPrint({ summary(spatial_models()$SAR) }) output$summary_sem <- renderPrint({ summary(spatial_models()$SEM) }) output$summary_sdm <- renderPrint({ summary(spatial_models()$SDM) }) output$summary_sdem <- renderPrint({ summary(spatial_models()$SDEM) }) output$summary_sac <- renderPrint({ summary(spatial_models()$SAC) })

plot_spatial_pred <- function(model_obj, title_txt) { sf_use <- spatial_models()$sf sf_use$pred <- as.numeric(fitted(model_obj))

ggplot(sf_use) +
  geom_sf(aes(fill = pred), color = "white", linewidth = 0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.2,
    check_overlap = TRUE
  ) + 
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(title = title_txt, fill = "Prediksi")

}

output$map_slx <- renderPlot({ plot_spatial_pred(spatial_models()$SLX, “Prediksi SLX”) }) output$map_sar <- renderPlot({ plot_spatial_pred(spatial_models()$SAR, “Prediksi SAR”) }) output$map_sem <- renderPlot({ plot_spatial_pred(spatial_models()$SEM, “Prediksi SEM”) }) output$map_sdm <- renderPlot({ plot_spatial_pred(spatial_models()$SDM, “Prediksi SDM”) }) output$map_sdem <- renderPlot({ plot_spatial_pred(spatial_models()$SDEM, “Prediksi SDEM”) }) output$map_sac <- renderPlot({ plot_spatial_pred(spatial_models()$SAC, “Prediksi SAC”) })

# ============================================================ # 15) Geographically Weighted Regression (GWR) # ============================================================

gwr_model <- reactive({ req(input$run)

sfdat <- jatim_join()
var   <- input$outcome

sf_use <- sfdat %>%
  filter(
    !is.na(.data[[var]]),
    !is.na(Penduduk_Miskin),
    !is.na(Sanitasi_Layak)
  )

validate(need(nrow(sf_use) >= 15, "Data terlalu sedikit untuk GWR"))

# PROYEKSI KE METRIK (WAJIB)
sf_use <- st_transform(sf_use, 32749)
pts_sp <- as(st_centroid(sf_use), "Spatial")

df <- sf_use %>%
  st_drop_geometry() %>%
  select(all_of(var), Penduduk_Miskin, Sanitasi_Layak)

bw <- gwr.sel(
  as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
  data = df,
  coords = coordinates(pts_sp),
  adapt = TRUE
)

gwr(
  as.formula(paste(var, "~ Penduduk_Miskin + Sanitasi_Layak")),
  data = df,
  coords = coordinates(pts_sp),
  adapt = bw
)

})

# ============================================================ # Residual GWR + Bobot Spasial (kNN) # ============================================================

gwr_residuals <- reactive({ gwr_res <- gwr_model()

# ambil residual GWR
res <- as.numeric(gwr_res$SDF$gwr.e)

# buang NA
res <- res[!is.na(res)]

validate(
  need(length(res) >= 3, "Residual GWR terlalu sedikit untuk uji asumsi")
)

res

})

# ============================================================ # Uji Normalitas Residual GWR (AMAN) # ============================================================

output$tbl_gwr_normality <- renderDT({ res <- gwr_residuals()

test <- shapiro.test(res)

out <- data.frame(
  Test = "Shapiro-Wilk",
  Statistic = round(test$statistic, 3),
  P_value = round(test$p.value, 3),
  Keputusan = ifelse(test$p.value > 0.05,
                     "Normal",
                     "Tidak normal")
)

datatable(out, options = list(dom = "t"))

})

# ============================================================ # Moran’s I Residual GWR # ============================================================

output$tbl_gwr_moran <- renderDT({ res <- gwr_residuals()

sf_use <- st_transform(jatim_join(), 32749)
coords <- st_coordinates(st_centroid(sf_use))

nb <- knn2nb(knearneigh(coords, k = 4))
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

mor <- moran.test(res, listw, zero.policy = TRUE)

out <- data.frame(
  Statistik = c("Moran's I", "p-value"),
  Nilai = round(c(mor$estimate[1], mor$p.value), 3)
)

datatable(out, options = list(dom = "t"))

})

output$summary_gwr <- renderPrint({ summary(gwr_model()) })

output$map_gwr_coef <- renderPlot({ gwr_res <- gwr_model() sf_use <- st_transform(jatim_join(), 32749)

gwr_df <- as.data.frame(gwr_res$SDF)

# pastikan urutan cocok
sf_use$coef_miskin <- gwr_df$Penduduk_Miskin
sf_use$coef_sanit  <- gwr_df$Sanitasi_Layak

p1 <- ggplot(sf_use) +
  geom_sf(aes(fill = coef_miskin), color = "white", linewidth = 0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.1,
    check_overlap = TRUE
  ) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(
    title = "GWR — Koefisien Lokal Penduduk Miskin",
    fill = "Koefisien"
  )

p2 <- ggplot(sf_use) +
  geom_sf(aes(fill = coef_sanit), color = "white", linewidth = 0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.1,
    check_overlap = TRUE
  ) +
  scale_fill_viridis_c(option = "viridis") +
  theme_minimal() +
  labs(
    title = "GWR — Koefisien Lokal Sanitasi Layak",
    fill = "Koefisien"
  )

gridExtra::grid.arrange(p1, p2, ncol = 2)

})

output$map_gwr_pred <- renderPlot({ gwr_res <- gwr_model() sf_use <- st_transform(jatim_join(), 32749)

sf_use$pred <- gwr_res$SDF$pred

ggplot(sf_use) +
  geom_sf(aes(fill = pred), color = "white", linewidth = 0.2) +
  geom_sf_text(
    aes(label = NAME_2),
    color = "white",
    size = 2.2,
    check_overlap = TRUE
  ) +
  scale_fill_viridis_c(option = "magma") +
  theme_minimal() +
  labs(
    title = "Prediksi GWR",
    fill = "Prediksi"
  )

}) }

============================================================

🚀 Run App

============================================================

shinyApp(ui, server)

Lampiran Link Dashboard

https://tiffannyspatialstatistics.shinyapps.io/spasial_epidem_uas/

Ditujukan Guna Memenuhi UAS Mata Kuliah Spatial Statistics

Latar Belakang

Identifikasi Masalah

Tujuan Penelitian

Batasan Penelitian

Ketergantungan Spasial (Spatial Dependence)

Autokorelasi Spasial

Model Ekonometrika Spasial

Geographically Weighted Regression (GWR)

Interpolasi Spasial

Inverse Distance Weighting (IDW)

Kriging

Ordinary Kriging

Universal Kriging

Sumber dan Jenis Data

Unit Spasial

Metode Analisis

Prevalensi

Autokorelasi Spasial Global (Moran’s I)

Autokorelasi Spasial Lokal (LISA)

Model Regresi Linier OLS

Model Regresi Spasial Global

SLX

SAR

SEM

SDM

SDEM

SAC

Pemilihan Model terbaik berdasarkan AIC

Geographically Weighted Regression (GWR)

Interpolasi Spasial

Inverse Distance Weighting (IDW)

Kriging

Ordinary Kriging

Universal Kriging

SLX

SAR

SEM

SDM

SDEM

SAC

Kesimpulan

Saran

Lampiran Script Analisis R

📊 Dashboard Interaktif Analisis Spasial Campak (Jawa Timur) 2022

Materi: Peta, Moran, Geary, LISA, Gi*, IDW, Variogram, OK, UK, Cross-Variogram (LMC)

============================================================

🧭 THEME

============================================================

============================================================

🔧 Helper

============================================================

============================================================

🧱 UI

============================================================

============================================================

⚙️ SERVER

============================================================

============================================================

🚀 Run App

============================================================

Lampiran Link Dashboard

Lampiran Link Video Social Media

**Ketergantungan Spasial (Spatial Dependence)**