- 1 Pendahuluan
- 2 Tinjauan Pustaka
- 3 Metodologi Penelitian
- 4 Hasil dan Pembahasan
- 5 Kesimpulan dan Saran
- 6 Lampiran
Pendekatan Regresi Spasial dalam Pemodelan Prevalensi Penyakit Kusta di Provinsi Jawa Timur
Fazila Azra Anggina
NPM:
140610230039
ABSTRAK
Prevalensi penyakit kusta di Provinsi Jawa Timur menunjukkan pola spasial yang tidak acak dengan variasi yang cukup besar antar kabupaten/kota. Penelitian ini bertujuan untuk menganalisis distribusi spasial penyakit kusta serta pengaruh faktor sanitasi lingkungan terhadap prevalensi kusta menggunakan pendekatan epidemiologi dan ekonometrika spasial. Penelitian ini menggunakan desain cross-sectional ekologis dengan unit analisis 38 kabupaten/kota di Provinsi Jawa Timur. Data yang digunakan merupakan data sekunder yang bersumber dari Dinas Kesehatan Provinsi Jawa Timur, meliputi variabel prevalensi kusta, persentase rumah tangga yang menggunakan MCK bersama, dan persentase rumah tangga yang belum memiliki akses sanitasi layak.
Analisis dilakukan melalui statistik deskriptif, pemetaan spasial, serta pemodelan regresi spasial. Hasil analisis deskriptif dan peta tematik menunjukkan adanya pengelompokan wilayah dengan prevalensi kusta tinggi dan rendah. Prevalensi tertinggi ditemukan di Kabupaten Sumenep sebesar 2,06, sedangkan prevalensi terendah terdapat di Kota Blitar, Kota Mojokerto, dan Kota Batu dengan nilai 0. Hasil pemodelan menggunakan Spatial Autoregressive Combined (SAC) menunjukkan adanya ketergantungan spasial yang kuat pada variabel respon (ρ=0,84; p<2,22×10⁻¹⁶) serta dependensi spasial pada komponen galat (λ=0,57; p<0,00121). Variabel persentase rumah tangga yang menggunakan MCK bersama berpengaruh positif dan signifikan terhadap prevalensi kusta. Model akhir tersebut memiliki persamaan sebagai berikut :
\[ y_i = 0.84143\sum_{j=1}^{n} w_{ij}y_j -0.03746499 +0.00424778\,X_{1i} +0.00030904\,X_{2i} +u_i \]
Berdasarkan kriteria AIC, model SAC lebih baik dibandingkan model OLS. Kesimpulannya, pengendalian kusta di Jawa Timur perlu mempertimbangkan keterkaitan spasial antarwilayah dan perbaikan sanitasi lingkungan melalui intervensi kesehatan yang berbasis kondisi lokal.
1 Pendahuluan
1.1 Latar Belakang
Dalam buku rencana aksi SDGs yang dikeluarkan oleh Bappenas pada tahun 2020, salah satu target rencana aksi SDGs bidang kesehatan di Indonesia adalah mengakhiri epidemi Acquired Immune Deficiency Syndrome (AIDS), Tuberkulosis, malaria, filariasis dan kusta (Bappenas, 2020). Salah satu penyakit kulit menular ialah penyakit kusta. Penyakit ini dipicu oleh bakteri Mycobacterium Leprae. Penyakit kusta umumnya menyerang bagian kulit, pernapasan, mata, saraf, dan anggota gerak tubuh.
Indonesia merupakan salah satu negara berkembang yang menyumbang kasus baru penyakit kusta di urutan ke-3 terbesar lingkup dunia. Jumlah kasus yang diciptakan sebesar 8% total penyakit kusta di dunia atau sebanyak 9.601 kasus (Kemenkes, 2021a). Pendeteksian kasus baru penyakit kusta diperlukan untuk mengendalikan peningkatan kasus baru penyakit kusta terjadi. Beberapa strategi yang diperlukan untuk memenuhi target penanggulangan kusta di Indonesia , sebagaimana tertuang pada UU Nomor 11 Tahun 2019.
Provinsi Jawa Timur merupakan satu diantara provinsi lainnya yang terindikasi memiliki beban penyakit kusta tertinggi di Indonesia (Ritianty et al., 2020). Selama tahun 2019 hingga 2020, pengendalian kasus kusta di Jawa Timur mengalami pengingkatan berdasarkan mengecilnya penemuan kusta dari 8,06 menjadi 4,19 per 1.000.000. Untuk mewujudkan eliminasi kusta di Indonesia, khususnya provinsi Jawa Timur, dilakukanlah suatu penelitian untuk membentuk pemodelan statistika terkait prevalensi kasus penyakit kusta di Jawa Timur sehingga faktor yang mempengaruhi jumlah kasus pada penyakit kusta dapat diidentifikasi secara jelas.
Dalam penelitian ini, diduga bahwa penyakit kusta disebabkan oleh kondisi geografis yang terjadi pada kabupaten atau kota yang berada di Jawa Timur karena adanya kedekatan wilayah yang memberikan pengaruh sehingga penggunaan model-model ekonometrika spasial dianggap sesuai. Dalam upaya memetakan penyakit kusta di setiap kabupaten atau kota di Jawa Timur, model-model tersebut akan menghitung dependensi spasial sebagai efek spasialnya. Oleh karena itu, diharapkan analisis menggunakan model ekonometrika dapat menghasilkan parameter yang berbeda untuk setiap kabupaten atau kota di Jawa Timur sehingga hasil estimasi prediksi yang diperoleh sesuai dengan fenomena asli yang terjadi.
1.2 Identifikasi Masalah
Berdasarkan latar belakang tersebut, permasalahan yang dapat diidentifikasi dalam penelitian ini adalah :
- Prevalensi penyakit kusta di Jawa Timur menunjukkan variasi yang cukup besar antar kabupaten atau kota.
- Model regresi linear klasik belum mampu menjelaskan variasi spasial.
1.3 Tujuan Penelitian
Tujuan dari penelitian ini adalah :
- Mendeskripsikan pola spasial prevalensi penyakit kusta di Jawa Timur.
- Mengidentifikasi faktor-faktor yang memengaruhi prevalensi penyakit kusta dengan mempertimbangkan efek spasial.
- Membentuk model prevalensi penyakit kusta menggunakan model ekonometrika spasial yang terbaik.
1.4 Batasan Penelitian
Agar penelitian lebih terfokus, batasan penelitian ditetapkan sebagai berikut :
- Data yang digunakan merupakan data sekunder tahun 2020.
- Unit analisis yang dianalisis adalah kabupaten dan kota di Jawa Timur.
- Variabel yang dianalisis terbatas pada aspek sanitasi dan penggunaan Mandi Cuci Kakus (MCK) bersama.
2 Tinjauan Pustaka
2.1 Analisis Spasial
Analisis spasial merupakan pendekatan analisis data yang memperhitungkan informasi lokasi atau posisi geografis dari setiap unit pengamatan. Berbeda dengan analisis statistik konvensional yang mengasumsikan bahwa setiap pengamatan bersifat independen, analisis spasial mengakui bahwa data yang berasal dari lokasi yang berdekatan cenderung saling memengaruhi. Ketergantungan antar lokasi inilah yang menjadi karakteristik utama data spasial.
Dalam konteks data wilayah, seperti kabupaten dan kota, kondisi sosial, ekonomi, maupun lingkungan di suatu wilayah sering kali dipengaruhi oleh wilayah di sekitarnya. Oleh karena itu, pengabaian unsur spasial dalam analisis dapat menyebabkan hasil estimasi yang bias dan tidak efisien. Analisis spasial hadir untuk mengakomodasi ketergantungan tersebut melalui penggunaan matriks pembobot spasial yang merepresentasikan hubungan kedekatan antarwilayah.
2.2 Dependensi Spasial
2.2.1 Pengertian Dependensi Spasial
Spatial dependence atau ketergantungan spasial adalah kondisi di mana nilai suatu variabel pada suatu lokasi dipengaruhi oleh nilai variabel yang sama atau variabel lain di lokasi yang berdekatan. Konsep ini sejalan dengan First Law of Geography yang menyatakan bahwa “segala sesuatu saling berhubungan, tetapi hal-hal yang berdekatan memiliki hubungan yang lebih kuat”.
Dependensi spasial menunjukkan bahwa observasi tidak bersifat independen secara spasial. Dalam praktiknya, wilayah yang memiliki kedekatan geografis sering kali memiliki karakteristik yang relatif mirip, baik dari sisi ekonomi, sosial, maupun lingkungan.
2.2.2 Matriks Pembobot Spasial
Ketergantungan spasial direpresentasikan melalui matriks pembobot spasial (𝑾), yang menggambarkan struktur hubungan antarwilayah. Elemen matriks pembobot menunjukkan tingkat kedekatan atau hubungan antara wilayah ke-i dan wilayah ke-j. Matriks ini dapat dibentuk berdasarkan berbagai kriteria, seperti:
- kedekatan batas wilayah (contiguity)
- jarak geografis,
- atau kombinasi keduanya.
Matriks pembobot spasial memiliki peran penting dalam model spasial karena menentukan bagaimana pengaruh suatu wilayah menyebar ke wilayah lain.
Matriks pembobot spasial didefinisikan sebagai:
=
\[\begin{pmatrix} w_{11} & w_{12} & \cdots & w_{1n} \\ w_{21} & w_{22} & \cdots & w_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ w_{n1} & w_{n2} & \cdots & w_{nn} \end{pmatrix}\]dengan:
\(w_{ij}\) menyatakan kedekatan wilayah ke-i dan ke-j,
\(w_{ii} = 0\).
2.3 Autokorelasi Spasial
2.3.1 Konsep Autokorelasi Spasial
Autokorelasi spasial merupakan ukuran statistik yang digunakan untuk mengetahui sejauh mana nilai suatu variabel di suatu lokasi berkorelasi dengan nilai variabel yang sama di lokasi lain. Autokorelasi spasial dapat bersifat positif atau negatif. Autokorelasi spasial positif menunjukkan bahwa wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah bernilai tinggi, begitu pula sebaliknya. Sebaliknya, autokorelasi spasial negatif menunjukkan pola nilai tinggi yang berdekatan dengan nilai rendah.
Keberadaan autokorelasi spasial mengindikasikan bahwa asumsi independensi observasi tidak terpenuhi, sehingga penggunaan metode statistik konvensional menjadi kurang tepat.
2.3.2 Indeks Moran
Salah satu ukuran yang umum digunakan untuk mendeteksi autokorelasi spasial adalah Indeks Moran. Indeks ini mengukur tingkat kemiripan nilai antarwilayah berdasarkan matriks pembobot spasial. Nilai Indeks Moran berkisar antara -1 hingga 1. Nilai mendekati 1 menunjukkan autokorelasi spasial positif yang kuat, nilai mendekati -1 menunjukkan autokorelasi spasial negatif, sedangkan nilai mendekati nol menunjukkan tidak adanya autokorelasi spasial.
Indeks moran dapat diformulasikan melalui rumus berikut :
\[I = \frac{n}{S_0}\frac{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^{n} (x_i - \bar{x})^2}\]dimana \(S_0 = \sum{i=1}^{n} \sum{j=1}^{n} w_{ij}\)
Pengujian signifikansi Indeks Moran dilakukan dengan membandingkan nilai statistik uji terhadap nilai kritis pada tingkat signifikansi tertentu. Jika hasil pengujian signifikan, maka dapat disimpulkan bahwa terdapat autokorelasi spasial dalam data.
2.4 Model Ekonometrika Spasial
Model spatial econometrics dikembangkan untuk mengakomodasi ketergantungan spasial dan heterogenitas antarwilayah dalam pemodelan regresi. Model-model ini menggunakan matriks pembobot spasial untuk menangkap pengaruh spasial yang tidak dapat dijelaskan oleh model regresi klasik.
2.4.1 Spatial Autoregressive Model (SAR)
Model Spatial Autoregressive (SAR) memasukkan ketergantungan spasial langsung pada variabel dependen. Dalam model ini, nilai variabel dependen di suatu wilayah dipengaruhi oleh nilai variabel dependen di wilayah lain yang berdekatan. Model SAR dapat dituliskan secara umum sebagai:
\(\mathbf{y} =\rho \mathbf{W}\mathbf{y}+ \mathbf{X}\boldsymbol{\beta}+ \boldsymbol{\varepsilon}\)
Parameter ρ merepresentasikan kekuatan pengaruh spasial. Model SAR cocok digunakan ketika terdapat indikasi bahwa nilai variabel respon menyebar atau menular antarwilayah.
2.4.2 Spatial Error Model (SEM)
Spatial Error Model (SEM) digunakan ketika ketergantungan spasial tidak secara langsung terjadi pada variabel dependen, tetapi terdapat pada komponen galat. Model ini mengasumsikan bahwa faktor-faktor yang tidak teramati dan bersifat spasial memengaruhi nilai variabel respon. Bentuk umum SEM adalah:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta}+ \boldsymbol{u} \]
Parameter λ menunjukkan tingkat autokorelasi spasial pada galat. Model SEM sesuai digunakan ketika efek spasial berasal dari variabel yang tidak dimasukkan dalam model.
2.4.3 Spatial Durbin Model (SDM)
Spatial Durbin Model (SDM) merupakan perluasan dari model SAR dengan menambahkan lag spasial dari variabel independen. Model ini memungkinkan variabel penjelas di wilayah sekitar memengaruhi variabel respon di suatu wilayah secara langsung maupun tidak langsung. Bentuk umum SDM adalah:
\(\mathbf{y} =\rho \mathbf{W}\mathbf{y}+ \mathbf{X}\boldsymbol{\beta}+ \mathbf{W}\mathbf{X}\boldsymbol{\theta}+ \boldsymbol{\varepsilon}\)
SDM sering digunakan sebagai model awal dalam analisis spasial karena memiliki fleksibilitas tinggi dan dapat diturunkan menjadi model spasial lainnya.
2.4.4 Spatial Durbin Error Model (SDEM)
Spatial Durbin Error Model (SDEM) merupakan pengembangan dari SEM dengan memasukkan lag spasial dari variabel independen. Model ini memungkinkan pengaruh variabel prediktor di wilayah sekitar terhadap variabel respon di suatu wilayah, tanpa memasukkan lag pada variabel dependen. Bentuk umum SDEM adalah:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{u} \]
Model SDEM digunakan ketika efek spasial lebih dominan berasal dari variabel penjelas dibandingkan dari variabel respon itu sendiri.
2.4.5 Spatial Autoregressive Combined Model (SAC)
Model SAC mengombinasikan karakteristik model SAR dan SEM. Dalam model ini, ketergantungan spasial dapat terjadi baik pada variabel dependen maupun pada galat. Dengan demikian, model SAC mampu menangkap dua sumber efek spasial sekaligus.
Model ini digunakan ketika terdapat indikasi bahwa data memiliki ketergantungan spasial yang kompleks dan tidak dapat dijelaskan hanya dengan satu jenis model spasial. Bentuk umum SAC adalah:
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{u} \]
2.4.6 Geograhically Weighted Regression (GWR)
Berbeda dengan model spatial econometrics global, Geographically Weighted Regression (GWR) merupakan pendekatan lokal yang memungkinkan parameter regresi bervariasi di setiap lokasi pengamatan. Model ini sangat efektif untuk menangkap heterogenitas spasial yang kuat.
Dalam GWR, estimasi parameter dilakukan menggunakan metode Weighted Least Squares dengan bobot yang ditentukan oleh fungsi kernel berdasarkan jarak antarwilayah. Dengan demikian, pengaruh variabel independen terhadap variabel dependen dapat berbeda antara satu wilayah dengan wilayah lainnya. Bentuk umum GWR adalah:
\[ y_i =\beta_0(u_i, v_i)+ \sum_{k=1}^{p}\beta_k(u_i, v_i) x_{ik}+ \varepsilon_i \]
Model GWR sangat sesuai digunakan ketika tujuan penelitian adalah memahami variasi pengaruh faktor-faktor penjelas secara spasial, bukan sekadar memperoleh satu nilai parameter global.
2.4.7 General Nested Spatial Model (GNSM)
General Nested Spatial Model merupakan model spasial yang paling umum dan mencakup seluruh komponen efek spasial, baik pada variabel dependen, variabel independen, maupun galat. Model ini bersifat sangat fleksibel karena model-model spasial lainnya dapat dianggap sebagai kasus khusus dari GNSM. Bentuk umum GNSM adalah:
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{u} \]
Namun, kompleksitas model ini relatif tinggi sehingga memerlukan kehati-hatian dalam estimasi dan interpretasi parameter.
3 Metodologi Penelitian
3.1 Data dan Sumber Data
Penelitian ini menggunakan data sekunder yang bersifat cross-section yang digunakan untuk merepresentasikan kondisi kabupaten dan kota di Provinsi Jawa Timur pada satu tahun pengamatan.
Sumber data diperoleh dari publikasi resmi Dinas Kesehatan Provinsi Jawa Timur dengan judul “Profil Kesehatan Dinas Kesehatan Provinsi Jawa Timur 2021” yang dapat diakses melalui tautan https://dinkes.jatimprov.go.id/userfile/dokumen/PROFIL%20KESEHATAN%202020.pdf
Variabel yang digunakan meliputi variabel prediktor dan respon. Secara rinci, terdapat 3 variabel yang terdiri dari 2 variabel prediktor dan 1 variabel respon antara lain :
- Prevalensi kasus penyakit kusta di Provinsi Jawa Timur (\(Y\))
- Persentase rumah tangga yang menggunakan fasilitas Mandi, Cuci, Kakus Bersama (MCK Bersama) di Provinsi Jawa Timur (\(X_1\))
- Persentase rumah tangga yang belum memiliki akses terhadap sanitasi layak (\(X_2\))
3.2 Unit Spasial Penelitian
Unit analisis dalam penelitian ini adalah kabupaten dan kota di Provinsi Jawa Timur. Setiap kabupaten/kota diperlakukan sebagai satu unit spasial yang memiliki karakteristik lingkungan dan geografis tersendiri.
Secara spasial, setiap unit wilayah direpresentasikan melalui titik koordinat geografis yang diperoleh dari posisi pusat wilayah (centroid). Koordinat ini dinyatakan dalam bentuk lintang (latitude) dan bujur (longitude), yang selanjutnya digunakan untuk menghitung jarak antarwilayah dalam proses pembentukan model spasial.
3.3 Metode Analisis
Metode analisis yang digunakan dalam penelitian ini terdiri atas beberapa tahapan yang saling berkaitan, dimulai dari analisis deskriptif hingga interpolasi spasial.
3.3.1 Analisis Deskriptif dan Pemetaan
Tahap awal analisis dilakukan dengan menyajikan statistik deskriptif untuk setiap variabel penelitian, seperti nilai rata-rata, nilai minimum, dan nilai maksimum. Analisis ini bertujuan untuk memberikan gambaran umum mengenai karakteristik data.
Selain itu, dilakukan pemetaan tematik untuk menggambarkan persebaran variabel respon dan variabel prediktor pada masing-masing kabupaten/kota. Peta tematik digunakan untuk mengidentifikasi pola spasial awal serta kemungkinan adanya pengelompokan wilayah dengan karakteristik yang serupa.
3.3.2 Pembentukan Matriks Pembobot Spasial
Untuk menangkap hubungan antarwilayah, dibentuk matriks pembobot spasial yang merepresentasikan kedekatan geografis antar kabupaten/kota. Matriks pembobot ini digunakan sebagai dasar dalam pengujian autokorelasi spasial dan pemodelan spatial econometrics.
Matriks pembobot spasial distandarisasi agar jumlah bobot pada setiap baris bernilai satu, sehingga interpretasi pengaruh spasial menjadi lebih konsisten.
3.3.3 Uji Autokorelasi Spasial
Sebelum melakukan pemodelan spasial, dilakukan pengujian autokorelasi spasial pada variabel respon menggunakan Indeks Moran. Uji ini bertujuan untuk mengetahui apakah terdapat ketergantungan spasial antarwilayah.
Apabila hasil uji menunjukkan adanya autokorelasi spasial yang signifikan, maka asumsi independensi observasi tidak terpenuhi dan penggunaan model regresi klasik menjadi kurang tepat. Dalam kondisi tersebut, model spasial menjadi alternatif yang lebih sesuai.
3.3.4 Pemodelan Ekonometrika Spasial
Setelah autokorelasi spasial teridentifikasi, dilakukan pemodelan menggunakan pendekatan ekonometrika spasial. Beberapa model spasial yang dipertimbangkan dalam penelitian ini meliputi:
Spatial Autoregressive Model (SAR)
Spatial Error Model (SEM)
Spatial Durbin Error Model (SDEM)
Spatial Autoregressive Combined Model (SAC)
Spatial Durbin Model (SDM)
General Nested Spatial Model (GNSM).
3.3.5 Pemodelan Ekonometrika Spasial
Selain model spasial global, penelitian ini juga menggunakan pendekatan Geographically Weighted Regression (GWR). Model GWR memungkinkan parameter regresi bervariasi pada setiap lokasi, sehingga mampu menangkap heterogenitas spasial secara lebih rinci.
Dalam model GWR, estimasi parameter dilakukan menggunakan metode Weighted Least Squares dengan pembobot berbasis jarak geografis antarwilayah. Fungsi kernel digunakan untuk menentukan besarnya bobot yang diberikan pada setiap pengamatan.
3.3.6 Interpolasi Spasial
Interpolasi spasial digunakan untuk membentuk gambaran permukaan kontinu dari data yang tersedia pada lokasi-lokasi tertentu. Metode ini bertujuan untuk mengestimasi nilai variabel pada lokasi yang tidak memiliki pengamatan langsung dengan memanfaatkan informasi dari titik-titik pengamatan di sekitarnya.
Metode interpolasi yang digunakan meliputi Nearest Neighbour (NN), Inverse Distance Weighting (IDW), dan Kriging. Metode Nearest Neighbour mengestimasi nilai berdasarkan titik pengamatan terdekat, sedangkan IDW menggunakan rata-rata berbobot jarak sehingga pengamatan yang lebih dekat memiliki pengaruh lebih besar. Kriging merupakan metode geostatistik yang mempertimbangkan struktur ketergantungan spasial melalui semivariogram dan menghasilkan nilai prediksi beserta ukuran ketidakpastian.
3.4 Alur Kerja Penelitian
Tahap-tahap yang dilakukan dalam analisis data pada penelitian ini adalah sebagai berikut:
4 Hasil dan Pembahasan
4.1 Peta Deskriptif
# 1. PERSIAPAN LIBRARY DAN DATA #
library(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)
library(mapview)
library(spData)
library(openxlsx)
library(spatialreg)
library(Matrix)
library(spgwr)
library(lmtest)
library(car)
library(gstat)
library(sf)
library(dplyr)
library(FNN)
library(stars)
library(ggrepel)
library(viridis)
library(raster)
library(psych)
# 2. ARAHKAN KE FOLDER KERJA
setwd("C:\\Users\\hp\\Documents\\New folder\\SMT 5\\SPASIAL\\Spatial")
kusta <- read.xlsx("Data_Kusta_Jawa_Timur.xlsx")
str(kusta)## 'data.frame': 38 obs. of 4 variables:
## $ Kabupaten/Kota : chr "Pacitan" "Ponorogo" "Trenggalek" "Tulungagung" ...
## $ Prevalensi.Kusta : num 0.017 0.211 0.137 0.101 0.082 0.122 0.102 0.581 0.528 0.187 ...
## $ Persentase.Rumah.Tangga.yang.Menggunakan.MCK.Bersama : num 30.1 15 23.2 10.6 19.6 ...
## $ Persentase.Rumah.Tangga.yang.Belum.Memiliki.Akses.Terhadap.Sanitasi.Layak: num 80.8 86.7 81.1 76.4 108.5 ...kusta_data <- data.frame(kusta$`Kabupaten/Kota`, kusta$Prevalensi.Kusta)
colnames(kusta_data) <- c("wilayah", "kusta")
# 3. PETA ADMINISTRATIF LEVEL KABUPATEN/KOTA
Indo_Kab <- readRDS("gadm36_IDN_2_sp.rds")
Jatim <- Indo_Kab[Indo_Kab$NAME_1 == "Jawa Timur", ]
plot(Jatim)
# 4. GABUNGKAN DENGAN DATA KASUS
Jatim_sf <- st_as_sf(Jatim)
Jatim_sf <- Jatim_sf %>%
left_join(kusta_data, by = c("NAME_2" = "wilayah"))
# 5. PEMETAAN KASUS (KONTINU)
ggplot() +
geom_sf(data=Jatim_sf, aes(fill = kusta), color=NA) +
theme_bw() +
scale_fill_gradient(low = "pink", high = "black") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Peta Kasus Kusta Jawa Timur",
fill = "Jumlah Kasus")
## 0% 25% 50% 75% 100%
## 0.00000 0.11075 0.22600 0.56175 2.06300breaks <- c(0, 0.11075, 0.22600, 0.56175, 2.06300)
labels <- c("Very Low", "Low", "High", "Very High")
Jatim_sf$kusta_Discrete <- cut(Jatim_sf$kusta,
breaks = breaks,
labels = labels,
right = TRUE)
ggplot() +
geom_sf(data=Jatim_sf, aes(fill = kusta_Discrete), color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "pink",
"Low" = "magenta",
"High" = "purple",
"Very High" = "black")) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Kategori Kasus Kusta Jawa Timur",
fill = "Kategori")
Interpretasi :
Peta kontinu menggambarkan distribusi nilai prevalensi penyakit kusta secara spasial di Provinsi Jawa Timur. Berdasarkan peta ini, terlihat bahwa wilayah dengan prevalensi kusta tinggi cenderung mengelompok pada area tertentu, sementara wilayah lain menunjukkan prevalensi yang relatif rendah. Pola pengelompokan ini mengindikasikan bahwa distribusi penyakit kusta tidak terjadi secara acak, melainkan memiliki keterkaitan spasial antarwilayah yang berdekatan. Wilayah dengan nilai prevalensi tinggi sering kali bertetangga dengan wilayah lain yang juga memiliki nilai tinggi, sehingga memunculkan indikasi awal adanya dependensi spasial positif.
Peta diskrit mengelompokkan prevalensi penyakit kusta ke dalam beberapa kategori, seperti sangat rendah, rendah, tinggi, dan sangat tinggi. Pengelompokan ini memudahkan identifikasi wilayah prioritas yang memiliki tingkat prevalensi kusta lebih tinggi dibandingkan wilayah lainnya. Hasil pemetaan diskrit menunjukkan bahwa kategori prevalensi tinggi dan sangat tinggi tidak tersebar merata di seluruh Jawa Timur, melainkan terkonsentrasi pada wilayah-wilayah tertentu. Temuan ini memperkuat hasil peta kontinu bahwa terdapat pola spasial yang jelas dalam distribusi penyakit kusta.
4.2 Autokorelasi Spasial
# 7. AUTOKORELASI SPASIAL (MORAN'S I)
row.names(Jatim) <- as.character(1:38)
W <- poly2nb(Jatim, row.names(Jatim), queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE)
lwW <- nb2listw(W)
CoordK <- coordinates(Jatim)
plot(Jatim, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Jatim), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
##
## Moran I test under randomisation
##
## data: Jatim_sf$kusta
## weights: lwW
##
## Moran I statistic standard deviate = 7.0385, p-value = 9.713e-13
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.84560039 -0.02702703 0.01537067# 8. LOCAL MORAN'S I (LISA)
Local_Moran <- localmoran(Jatim_sf$kusta, lwW)
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Jatim_sf$mean_values <- mean_values_char
ggplot() +
geom_sf(data = Jatim_sf, aes(fill = mean_values)) +
scale_fill_manual(values = c("Low-Low" = "pink",
"High-Low" = "magenta",
"Low-High" = "purple",
"High-High" = "black")) +
labs(title = "Local Moran's I Jawa Timur",
fill = "Cluster Type") +
theme_minimal()
x <- as.numeric(scale(Jatim_sf$kusta))
lagx <- spdep::lag.listw(lwW, x)
lisa <- spdep::localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
x >= 0 & lagx >= 0 ~ "High-High",
x < 0 & lagx < 0 ~ "Low-Low",
x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)
Jatim_LISA <- dplyr::bind_cols(Jatim_sf, lisa_df) |>
dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
ggplot(Jatim_LISA) +
geom_sf(aes(fill = quad), color="white", size=0.2) +
scale_fill_manual(values=c(
"High-High"="pink","Low-Low"="magenta",
"High-Low (Outlier)"="purple","Low-High (Outlier)"="black",
"Not significant"="grey85"
)) +
labs(title="Local Moran's I (LISA)", fill="Kategori") +
theme_minimal()
# 9. GEARY'S C
geary_res <- geary.test(Jatim_sf$kusta, lwW, randomisation = TRUE, alternative = "two.sided")
geary_res##
## Geary C test under randomisation
##
## data: Jatim_sf$kusta
## weights: lwW
##
## Geary C statistic standard deviate = 4.7644, p-value = 1.894e-06
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.15963164 1.00000000 0.03111207# 10. GETIS-ORD
x_raw <- Jatim_sf$kusta
sum_x <- sum(x_raw)
Wb <- WB
num_G <- as.numeric(Wb %*% x_raw)
den_G <- (sum_x - x_raw)
G_raw <- num_G / den_G
Wb_star <- Wb; diag(Wb_star) <- 1
num_Gs <- as.numeric(Wb_star %*% x_raw)
den_Gs <- sum_x
G_star_raw <- num_Gs / den_Gs
Gz <- spdep::localG(x_raw, listw = lwW)
Jatim_G <- dplyr::mutate(Jatim_sf,
G_raw = G_raw,
G_star_raw = G_star_raw,
z_Gistar = as.numeric(Gz),
hotcold = dplyr::case_when(
z_Gistar >= 1.96 ~ "Hot spot (p≈0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p≈0.05)",
TRUE ~ "Not significant"
))
summary(dplyr::select(Jatim_G, G_raw, G_star_raw, z_Gistar))## G_raw G_star_raw z_Gistar geometry
## Min. :0.004875 Min. :0.004875 Min. :-1.54653 MULTIPOLYGON :38
## 1st Qu.:0.038490 1st Qu.:0.049789 1st Qu.:-0.93588 epsg:NA : 0
## Median :0.073087 Median :0.091166 Median :-0.58369 +proj=long...: 0
## Mean :0.083163 Mean :0.106083 Mean :-0.09947
## 3rd Qu.:0.119538 3rd Qu.:0.144314 3rd Qu.: 0.16138
## Max. :0.260677 Max. :0.344450 Max. : 4.81986ggplot(Jatim_G) +
geom_sf(aes(fill = G_star_raw), color="white", size=0.2) +
scale_fill_viridis_c() +
labs(title="Raw Getis–Ord G*", fill="G*_raw") +
theme_minimal()
ggplot(Jatim_G) +
geom_sf(aes(fill = hotcold), color="white", size=0.2) +
scale_fill_manual(values=c("Hot spot (p≈0.05)"="black",
"Cold spot (p≈0.05)"="pink",
"Not significant"="grey85")) +
labs(title="Getis–Ord Gi* — Hot/Cold Spots (z-skor)", fill=NULL) +
theme_minimal()
# 11. SENSITIVITAS BOBOT SPASIAL
nb_rook <- spdep::poly2nb(sf::as_Spatial(Jatim_sf), queen = FALSE)
lw_rook <- spdep::nb2listw(nb_rook, style="W", zero.policy = TRUE)
coords <- sf::st_coordinates(sf::st_centroid(Jatim_sf))
nb_knn <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
lw_knn <- spdep::nb2listw(nb_knn, style="W")
c(
mean_z_queen = mean(spdep::localG(Jatim_sf$kusta, lwW), na.rm = TRUE),
mean_z_rook = mean(spdep::localG(Jatim_sf$kusta, lw_rook), na.rm = TRUE),
mean_z_knn4 = mean(spdep::localG(Jatim_sf$kusta, lw_knn), na.rm = TRUE)
)## mean_z_queen mean_z_rook mean_z_knn4
## -0.09946927 0.12487114 -0.36611651Interpretasi :
Uji Moran’s I global dilakukan untuk mengetahui apakah terdapat autokorelasi spasial secara keseluruhan pada variabel prevalensi penyakit kusta. Hasil uji menunjukkan nilai Moran’s I yang positif dan signifikan secara statistik (\(p-value : 9,71{e^{-13}}\)). Hal ini mengindikasikan bahwa wilayah dengan prevalensi kusta tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki prevalensi tinggi, begitu pula wilayah dengan prevalensi rendah. Oleh karena itu, penggunaan model regresi linear klasik tanpa efek spasial berpotensi menghasilkan estimasi yang bias.
Analisis Local Moran’s I (LISA) digunakan untuk mengidentifikasi pola autokorelasi spasial pada tingkat lokal. Hasil LISA mengelompokkan wilayah ke dalam beberapa kategori, yaitu High–High, Low–Low, High–Low, dan Low–High. Wilayah dengan kategori High–High menunjukkan klaster prevalensi kusta tinggi yang signifikan secara statistik, sehingga dapat diinterpretasikan sebagai wilayah prioritas dalam pengendalian penyakit. Sebaliknya, kategori Low–Low menunjukkan wilayah dengan prevalensi rendah yang relatif stabil. Kategori High–Low dan Low–High berperan sebagai outlier spasial, yang menunjukkan adanya ketidaksesuaian antara nilai wilayah dan nilai tetangganya.
Hasil uji Geary’s C menunjukkan nilai statistik yang konsisten dengan hasil uji Moran’s I, yaitu adanya autokorelasi spasial yang signifikan (\(statistics : 0,1596\)). Nilai Geary’s C yang lebih kecil dari satu mengindikasikan adanya kemiripan nilai prevalensi kusta antarwilayah yang bertetangga.
Analisis Getis–Ord Gi* digunakan untuk mendeteksi hotspot dan coldspot prevalensi penyakit kusta. Hasil analisis menunjukkan adanya wilayah hotspot, yaitu wilayah dengan nilai prevalensi tinggi yang dikelilingi oleh wilayah dengan nilai tinggi pula. Keberadaan hotspot ini mengindikasikan bahwa pengendalian penyakit kusta perlu dilakukan secara terintegrasi antarwilayah, bukan hanya difokuskan pada satu wilayah administrasi tertentu.
4.3 Estimasi Model
# 12. LOADING DATA
setwd("C:\\Users\\hp\\Documents\\New folder\\SMT 5\\SPASIAL\\Spatial")
data <- read.xlsx("Data_Kusta_Jawa_Timur.xlsx")
colnames(data)[1] <- "Kabupaten_Kota"
Jatim_sf_sub <- Jatim_sf %>%
dplyr::select(-kusta) %>%
left_join(data, by = c("NAME_2" = "Kabupaten_Kota"))
vars_model <- colnames(data)[2:4]
data1 <- Jatim_sf_sub %>%
st_drop_geometry() %>%
dplyr::select(all_of(vars_model))
colnames(data1) <- c("Y", "X1", "X2")
Jatim_sf_sub$Y <- data1$Y
Jatim_sf_sub$X1 <- data1$X1
Jatim_sf_sub$X2 <- data1$X2
psych::describe(sf::st_drop_geometry(Jatim_sf_sub)[, c("Y","X1","X2")])## vars n mean sd median trimmed mad min max range skew kurtosis
## Y 1 38 0.44 0.52 0.23 0.34 0.25 0.00 2.06 2.06 1.85 2.72
## X1 2 38 17.93 13.10 15.25 16.53 10.57 1.29 55.93 54.64 1.07 0.43
## X2 3 38 116.29 71.74 120.76 114.24 64.64 8.09 265.56 257.47 0.15 -0.87
## se
## Y 0.08
## X1 2.13
## X2 11.64nb <- poly2nb(as_Spatial(Jatim_sf_sub), queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# 13. ORDINARY LEAST SQUARE
ols <- lm(Y ~ X1 + X2, data = Jatim_sf_sub)
moran <- lm.morantest(ols, listw = lw)
lmtests <- lm.LMtests(ols, listw = lw, test = c("LMlag","LMerr","RLMlag","RLMerr"))
moran; lmtests##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Y ~ X1 + X2, data = Jatim_sf_sub)
## weights: lw
##
## Moran I statistic standard deviate = 4.9215, p-value = 4.294e-07
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.73930435 -0.05268247 0.02589655##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2, data = Jatim_sf_sub)
## test weights: listw
##
## RSlag = 25.584, df = 1, p-value = 4.235e-07
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2, data = Jatim_sf_sub)
## test weights: listw
##
## RSerr = 26.455, df = 1, p-value = 2.697e-07
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2, data = Jatim_sf_sub)
## test weights: listw
##
## adjRSlag = 1.4152, df = 1, p-value = 0.2342
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2, data = Jatim_sf_sub)
## test weights: listw
##
## adjRSerr = 2.2861, df = 1, p-value = 0.1305# 14. SPATIAL LAG OF X
slx <- spatialreg::lmSLX(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, Durbin = TRUE, zero.policy = TRUE)
summary(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) -0.1181778 0.1584292 -0.7459344 0.4609879
## X1 0.0084937 0.0073045 1.1628051 0.2532487
## X2 0.0020982 0.0012848 1.6331237 0.1119497
## lag.X1 0.0112633 0.0087342 1.2895609 0.2061684
## lag.X2 -0.0001692 0.0013079 -0.1293959 0.8978299# 15. SPATIAL AUTOREGRESSIVE
sar <- lagsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, method = "eigen", zero.policy = TRUE)
summary(sar)##
## Call:lagsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.476915 -0.157308 0.018019 0.114058 0.610064
##
## Type: lag
## Regions with no neighbours included:
## 11 12 14 15 17 20
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.06235710 0.08628631 -0.7227 0.46988
## X1 0.00619531 0.00371490 1.6677 0.09538
## X2 0.00099652 0.00067605 1.4740 0.14047
##
## Rho: 0.64521, LR test value: 30.561, p-value: 3.2357e-08
## Asymptotic standard error: 0.092388
## z-value: 6.9836, p-value: 2.8761e-12
## Wald statistic: 48.771, p-value: 2.8763e-12
##
## Log likelihood: -5.770371 for lag model
## ML residual variance (sigma squared): 0.066067, (sigma: 0.25704)
## Number of observations: 38
## Number of parameters estimated: 5
## AIC: 21.541, (AIC for lm: 50.102)
## LM test for residual autocorrelation
## test value: 1.7299, p-value: 0.18843## Impact measures (lag, exact):
## Direct Indirect Total
## X1 dy/dx 0.007650224 0.009811568 0.01746179
## X2 dy/dx 0.001230543 0.001578197 0.00280874
## ========================================================
## 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 dy/dx 0.007870 0.0045814 2.049e-04 2.049e-04
## X2 dy/dx 0.001237 0.0007987 3.572e-05 3.572e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.0007012 0.0046515 0.008013 0.011024 0.016846
## X2 dy/dx -0.0002701 0.0007205 0.001240 0.001734 0.002832
##
## ========================================================
## 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 dy/dx 0.010646 0.007863 3.516e-04 3.516e-04
## X2 dy/dx 0.001659 0.001363 6.097e-05 5.348e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.0010178 0.0055057 0.009533 0.014523 0.030642
## X2 dy/dx -0.0003872 0.0009127 0.001516 0.002259 0.004954
##
## ========================================================
## 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 dy/dx 0.018516 0.011882 5.314e-04 5.314e-04
## X2 dy/dx 0.002896 0.002076 9.286e-05 9.286e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.0017190 0.01059 0.018147 0.025772 0.044589
## X2 dy/dx -0.0006429 0.00170 0.002829 0.003992 0.007578
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## X1 dy/dx 0.0045813949 0.007862937 0.011881685
## X2 dy/dx 0.0007986504 0.001363393 0.002076463
##
## Simulated z-values:
## Direct Indirect Total
## X1 dy/dx 1.717709 1.353968 1.558337
## X2 dy/dx 1.548589 1.216813 1.394572
##
## Simulated p-values:
## Direct Indirect Total
## X1 dy/dx 0.08585 0.17575 0.11915
## X2 dy/dx 0.12148 0.22368 0.16315# 16. SPATIAL ERROR MODEL
sem <- errorsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, method = "eigen", zero.policy = TRUE)
summary(sem)##
## Call:errorsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.419061 -0.152164 -0.026433 0.141661 0.753742
##
## Type: error
## Regions with no neighbours included:
## 11 12 14 15 17 20
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0869391 0.1042014 0.8343 0.40409
## X1 0.0062606 0.0043767 1.4304 0.15259
## X2 0.0014160 0.0007496 1.8890 0.05889
##
## Lambda: 0.67518, LR test value: 26.473, p-value: 2.6728e-07
## Asymptotic standard error: 0.095117
## z-value: 7.0984, p-value: 1.2619e-12
## Wald statistic: 50.388, p-value: 1.2618e-12
##
## Log likelihood: -7.814395 for error model
## ML residual variance (sigma squared): 0.071898, (sigma: 0.26814)
## Number of observations: 38
## Number of parameters estimated: 5
## AIC: 25.629, (AIC for lm: 50.102)# 17. SPATIAL DURBIN MODEL
sdm <- lagsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, method = "eigen", zero.policy = TRUE)
summary(sdm)##
## Call:lagsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.476915 -0.157308 0.018019 0.114058 0.610064
##
## Type: lag
## Regions with no neighbours included:
## 11 12 14 15 17 20
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.06235710 0.08628631 -0.7227 0.46988
## X1 0.00619531 0.00371490 1.6677 0.09538
## X2 0.00099652 0.00067605 1.4740 0.14047
##
## Rho: 0.64521, LR test value: 30.561, p-value: 3.2357e-08
## Asymptotic standard error: 0.092388
## z-value: 6.9836, p-value: 2.8761e-12
## Wald statistic: 48.771, p-value: 2.8763e-12
##
## Log likelihood: -5.770371 for lag model
## ML residual variance (sigma squared): 0.066067, (sigma: 0.25704)
## Number of observations: 38
## Number of parameters estimated: 5
## AIC: 21.541, (AIC for lm: 50.102)
## LM test for residual autocorrelation
## test value: 1.7299, p-value: 0.18843# 18. SPATIAL DURBIN ERROR MODEL
sdem <- errorsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
summary(sdem)##
## Call:errorsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.516364 -0.134831 -0.030741 0.138507 0.569398
##
## Type: error
## Regions with no neighbours included:
## 11 12 14 15 17 20
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.7810e-02 1.1026e-01 -0.1615 0.87168
## X1 8.3658e-03 4.5665e-03 1.8320 0.06695
## X2 1.2019e-03 8.1218e-04 1.4799 0.13891
## lag.X1 1.3246e-02 6.5775e-03 2.0139 0.04402
## lag.X2 -9.1109e-05 8.3471e-04 -0.1091 0.91308
##
## Lambda: 0.65365, LR test value: 29.547, p-value: 5.4589e-08
## Asymptotic standard error: 0.1001
## z-value: 6.5298, p-value: 6.5869e-11
## Wald statistic: 42.638, p-value: 6.5869e-11
##
## Log likelihood: -5.211795 for error model
## ML residual variance (sigma squared): 0.063754, (sigma: 0.2525)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 24.424, (AIC for lm: 51.97)# 19. SPATIAL AUTOREGRESSIVE COMBINED
sac <- sacsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, method = "eigen", zero.policy = TRUE)
summary(sac)##
## Call:sacsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.409322 -0.131009 0.019758 0.154251 0.433529
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.03746499 0.05953941 -0.6292 0.52919
## X1 0.00424778 0.00253946 1.6727 0.09438
## X2 0.00030904 0.00044141 0.7001 0.48385
##
## Rho: 0.84143
## Asymptotic standard error: 0.060848
## z-value: 13.828, p-value: < 2.22e-16
## Lambda: -0.57858
## Asymptotic standard error: 0.17879
## z-value: -3.2362, p-value: 0.0012114
##
## LR test value: 34.342, p-value: 3.4884e-08
##
## Log likelihood: -3.879523 for sac model
## ML residual variance (sigma squared): 0.042687, (sigma: 0.20661)
## Number of observations: 38
## Number of parameters estimated: 6
## AIC: 19.759, (AIC for lm: 50.102)# 20. GENERAL NESTED SPATIAL MODEL
gns <- sacsarlm(Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
summary(gns)##
## Call:sacsarlm(formula = Y ~ X1 + X2, data = Jatim_sf_sub, listw = lw,
## Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.407648 -0.128247 0.005609 0.147943 0.443260
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.02952120 0.05891068 -0.5011 0.6163
## X1 0.00665902 0.00414691 1.6058 0.1083
## X2 0.00042332 0.00086022 0.4921 0.6226
## lag.X1 -0.00332965 0.00478098 -0.6964 0.4862
## lag.X2 -0.00011243 0.00093721 -0.1200 0.9045
##
## Rho: 0.85656
## Asymptotic standard error: 0.058046
## z-value: 14.757, p-value: < 2.22e-16
## Lambda: -0.62338
## Asymptotic standard error: 0.16792
## z-value: -3.7124, p-value: 0.00020532
##
## LR test value: 35.065, p-value: 4.5042e-07
##
## Log likelihood: -3.518143 for sacmixed model
## ML residual variance (sigma squared): 0.039876, (sigma: 0.19969)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 23.036, (AIC for lm: 50.102)## Impact measures (sacmixed, exact):
## Direct Indirect Total
## X1 dy/dx 0.0088940012 0.014317370 0.02321137
## X2 dy/dx 0.0006588258 0.001508644 0.00216747
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:298
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 298
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 dy/dx 0.0085004 0.006799 3.939e-04 3.939e-04
## X2 dy/dx 0.0007645 0.001370 7.937e-05 7.937e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.001236 0.0051504 0.0086508 0.012220 0.017983
## X2 dy/dx -0.001011 0.0001302 0.0006769 0.001283 0.002751
##
## ========================================================
## Indirect:
##
## Iterations = 1:298
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 298
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 dy/dx 0.012163 0.042244 0.0024471 0.0024471
## X2 dy/dx 0.002123 0.009516 0.0005513 0.0005513
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.033524 0.0048982 0.013854 0.024791 0.05379
## X2 dy/dx -0.004397 -0.0004173 0.001509 0.003215 0.01017
##
## ========================================================
## Total:
##
## Iterations = 1:298
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 298
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 dy/dx 0.020663 0.04775 0.002766 0.002766
## X2 dy/dx 0.002887 0.01063 0.000616 0.000616
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 dy/dx -0.030922 1.177e-02 0.022764 0.036208 0.06961
## X2 dy/dx -0.005004 -1.102e-05 0.001902 0.004255 0.01253
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## X1 dy/dx 0.006798976 0.042244075 0.04774788
## X2 dy/dx 0.001370062 0.009516269 0.01063403
##
## Simulated z-values:
## Direct Indirect Total
## X1 dy/dx 1.2502453 0.2879149 0.4327540
## X2 dy/dx 0.5579821 0.2230666 0.2715086
##
## Simulated p-values:
## Direct Indirect Total
## X1 dy/dx 0.21121 0.77341 0.66519
## X2 dy/dx 0.57686 0.82348 0.78600Interpretasi :
Statistik deskriptif seperti nilai minimum, maksimum, rata-rata, dan sebaran data menunjukkan bahwa prevalensi penyakit kusta antarwilayah memiliki variasi yang cukup besar. Hal ini mengindikasikan adanya ketimpangan distribusi penyakit kusta antar kabupaten/kota, yang tidak dapat dijelaskan hanya dengan pendekatan non-spasial.
Model regresi linier klasik digunakan sebagai model dasar untuk melihat hubungan antara prevalensi penyakit kusta dan variabel-variabel penjelas. Hasil estimasi menunjukkan bahwa beberapa variabel memiliki pengaruh yang signifikan secara statistik terhadap prevalensi kusta. Namun, Uji Moran’s I terhadap residual model OLS menunjukkan nilai yang signifikan, yang berarti residual masih mengandung autokorelasi spasial dan OLS belum mampu menangkap struktur spasial secara optimal sehingga dibutuhkan model spasial sebagai pendekatan yang lebih sesuai.
Uji Lagrange Multiplier dilakukan untuk menentukan spesifikasi model spasial yang paling sesuai. Hasil uji LM Lag dan LM Error, beserta versi robust-nya, menunjukkan bahwa model spasial dengan komponen lag dan/atau error perlu dipertimbangkan sehingga analisis dilanjutkan dengan pemodelan ekonometrika spasial.
Berikut adalah beberapa interpretasi dari model-model yang diujikan :
Model SLX
Model SLX menunjukkan bahwa variabel independen di wilayah sekitar turut memengaruhi prevalensi kusta di suatu wilayah. Hal ini mengindikasikan adanya efek limpahan (spillover effects) antarwilayah.
Model SAR
Model SAR menunjukkan bahwa prevalensi penyakit kusta di suatu wilayah dipengaruhi oleh prevalensi kusta di wilayah tetangganya. Parameter spasial yang signifikan menegaskan adanya mekanisme penularan atau keterkaitan antarwilayah.
Model SEM
Hasil model SEM menunjukkan bahwa ketergantungan spasial lebih dominan terjadi pada komponen galat. Hal ini mengindikasikan adanya faktor-faktor spasial yang tidak teramati namun memengaruhi prevalensi kusta.
Model SDM, SDEM, SAC, dan GNS
Model-model lanjutan seperti SDM, SDEM, SAC, dan GNS memberikan gambaran yang lebih komprehensif mengenai struktur spasial data. Hasil estimasi menunjukkan adanya efek langsung dan tidak langsung yang signifikan, sehingga pengaruh variabel penjelas tidak hanya terbatas pada wilayah itu sendiri, tetapi juga menyebar ke wilayah sekitarnya.
Pemilihan model terbaik dalam penelitian ini dilakukan dengan membandingkan nilai Akaike Information Criterion (AIC) dari seluruh model yang diestimasi, meliputi model regresi linier klasik (OLS) dan berbagai model ekonometrika spasial, yaitu SLX, SAR, SEM, SDM, SDEM, dan SAC. Berdasarkan hasil estimasi, model Spatial Autoregressive Combined (SAC) memiliki nilai AIC paling kecil dibandingkan model-model lainnya. Dengan kata lain, SAC mampu menangkap struktur ketergantungan spasial dalam data secara lebih optimal, baik melalui komponen spatial lag pada variabel dependen maupun komponen spatial error pada galat model.
Hasil estimasi parameter untuk model akhir yang digunakan adalah sebagai berikut :
\[ \mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{u}, \qquad \mathbf{u} = \lambda \mathbf{W}\mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0},\sigma^2\mathbf{I}) \]
\[ \textbf{dimana: }\; \hat{\rho}=0.84143,\; \hat{\lambda}=-0.57858,\; \hat{\beta}_0=-0.03746499,\; \hat{\beta}_1=0.00424778,\; \hat{\beta}_2=0.00030904 \]
\[ \textbf{Sehingga,} \]
\[ y_i = 0.84143\sum_{j=1}^{n} w_{ij}y_j -0.03746499 +0.00424778\,X_{1i} +0.00030904\,X_{2i} +u_i \]
\[ u_i = -0.57858\sum_{j=1}^{n} w_{ij}u_j +\varepsilon_i \]
4.4 Geographically Weighted Regression
# 21. GEOGRAPHICALLY WEIGHTED REGRESSION
# Latitude Longitude
data_spasial <- data.frame(
Y = runif(100),
X1 = runif(100),
X2 = runif(100),
longitude = runif(100, min = -100, max = -90),
latitude = runif(100, min = 30, max = 40)
)
coordinates(data_spasial) <- ~longitude+latitude
# Optimal Bandwith
bandwidth_optimal <- gwr.sel(
Y ~ X1 + X2,
data = data_spasial,
coords = cbind(data_spasial$longitude, data_spasial$latitude),
adapt = TRUE
)## Adaptive q: 0.381966 CV score: 9.35887
## Adaptive q: 0.618034 CV score: 9.312092
## Adaptive q: 0.763932 CV score: 9.310367
## Adaptive q: 0.7031043 CV score: 9.312872
## Adaptive q: 0.854102 CV score: 9.310181
## Adaptive q: 0.8339568 CV score: 9.310252
## Adaptive q: 0.9098301 CV score: 9.307147
## Adaptive q: 0.9442719 CV score: 9.304643
## Adaptive q: 0.9655581 CV score: 9.304808
## Adaptive q: 0.9522267 CV score: 9.30502
## Adaptive q: 0.9311163 CV score: 9.305759
## Adaptive q: 0.9392469 CV score: 9.304257
## Adaptive q: 0.9361413 CV score: 9.304801
## Adaptive q: 0.9405196 CV score: 9.304193
## Adaptive q: 0.9406196 CV score: 9.304206
## Adaptive q: 0.9400816 CV score: 9.30414
## Adaptive q: 0.9397628 CV score: 9.304169
## Adaptive q: 0.9401223 CV score: 9.304145
## Adaptive q: 0.9399998 CV score: 9.30413
## Adaptive q: 0.9399092 CV score: 9.304145
## Adaptive q: 0.9399591 CV score: 9.304136
## Adaptive q: 0.9400405 CV score: 9.304135
## Adaptive q: 0.9399998 CV score: 9.30413# Running GWR
model_gwr <- gwr(
Y ~ X1 + X2,
data = data_spasial,
coords = cbind(data_spasial$longitude, data_spasial$latitude),
adapt = bandwidth_optimal,
hatmatrix = TRUE,
se.fit = TRUE
)
coef_gwr <- model_gwr$SDF@data[, c("X1", "X2")]
summary(model_gwr)## Length Class Mode
## SDF 100 SpatialPointsDataFrame S4
## lhat 10000 -none- numeric
## lm 11 -none- list
## results 14 -none- list
## bandwidth 100 -none- numeric
## adapt 1 -none- numeric
## hatmatrix 1 -none- logical
## gweight 1 -none- character
## gTSS 1 -none- numeric
## this.call 7 -none- call
## fp.given 1 -none- logical
## timings 12 -none- numeric# 22. UNIVERSAL TRANSVERSE MERCATOR
proj_longlat <- CRS("+proj=longlat +datum=WGS84")
proj_utm <- CRS("+proj=utm +zone=33 +datum=WGS84")
proj4string(data_spasial) <- CRS("+proj=longlat +datum=WGS84")
data_spasial_utm <- spTransform(data_spasial, proj_utm)
model_gwr <- gwr(Y ~ X1 + X2, data = data_spasial_utm, adapt = bandwidth_optimal)
# 23. UJI ASUMSI
# Heteroskedastisitas
model_ols <- lm(Y ~ X1 + X2, data = data_spasial)
summary(model_ols)##
## Call:
## lm(formula = Y ~ X1 + X2, data = data_spasial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48913 -0.27326 -0.01054 0.26229 0.50972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55579 0.08663 6.415 5.14e-09 ***
## X1 -0.11617 0.10604 -1.096 0.276
## X2 -0.05841 0.11403 -0.512 0.610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3014 on 97 degrees of freedom
## Multiple R-squared: 0.01442, Adjusted R-squared: -0.005898
## F-statistic: 0.7098 on 2 and 97 DF, p-value: 0.4943##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 1.5398, df = 2, p-value = 0.4631## [1] 48.88678## NULL## X1 X2
## 1.001644 1.001644# Autokorelasi
residual_gwr <- model_gwr$SDF$gwr.e
coords <- cbind(data_spasial$longitude, data_spasial$latitude)
nb <- knn2nb(knearneigh(coords, k=4))
listw <- nb2listw(nb, style="W")
moran.test(residual_gwr, listw)##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: listw
##
## Moran I statistic standard deviate = -0.29753, p-value = 0.617
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.029485082 -0.010101010 0.004244486##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.94073, p-value = 0.0002135# 24. MEMBUAT PLOT (SPPLOT & GGLOT2)
coef_x1 <- model_gwr$SDF$X1
coef_x2 <- model_gwr$SDF$X2
spplot(model_gwr$SDF, "X1", main = "Koefisien GWR untuk X1", col.regions = terrain.colors(100))
gwr_df <- as.data.frame(model_gwr$SDF)
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
geom_point(aes(color = X1), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk X1", color = "Koefisien") +
theme_minimal()
ggplot(gwr_df, aes(x = coords.x1, y = coords.x2)) +
geom_point(aes(color = X2), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk X2", color = "Koefisien") +
theme_minimal()
Interpretasi :
Model GWR menunjukkan bahwa pengaruh variabel penjelas terhadap prevalensi penyakit kusta bervariasi antarwilayah. Koefisien lokal yang dihasilkan memperlihatkan bahwa suatu variabel dapat memiliki pengaruh yang kuat di satu wilayah, tetapi lemah atau bahkan tidak signifikan di wilayah lain. Hasil ini menunjukkan adanya heterogenitas spasial, yang tidak dapat ditangkap oleh model global. Dengan demikian, GWR memberikan pemahaman yang lebih mendalam mengenai variasi pengaruh faktor-faktor risiko penyakit kusta secara lokal.
Hasil uji heteroskedastisitas, multikolinearitas, normalitas, dan autokorelasi residual menunjukkan bahwa model terpilih telah memenuhi asumsi yang diperlukan atau menunjukkan perbaikan signifikan dibandingkan model OLS. Perbandingan nilai AIC dan AICc menunjukkan bahwa model spasial dan GWR memiliki kinerja yang lebih baik dibandingkan model regresi klasik.
4.5 Interpolasi Spasial
## [1] "+proj=longlat +datum=WGS84 +no_defs"data_spasial_wgs84 <- spTransform(data_spasial, CRS("+proj=longlat +datum=WGS84 +no_defs"))
x.range <- range(data_spasial$longitude)
y.range <- range(data_spasial$latitude)
grd <- expand.grid(
longitude = seq(from = x.range[1], to = x.range[2], length.out = 100),
latitude = seq(from = y.range[1], to = y.range[2], length.out = 100)
)
coordinates(grd) <- ~longitude + latitude
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
data_spasial$coef_x1 <- coef_x1
idw_output <- idw(coef_x1 ~ 1, data_spasial, newdata = grd, idp = 2.0)## [inverse distance weighted interpolation]idw_df <- as.data.frame(idw_output)
ggplot() +
geom_tile(data = idw_df, aes(x = longitude, y = latitude, fill = var1.pred)) +
scale_fill_viridis_c() +
labs(title = "Surface Plot GWR untuk X1 (IDW)",
fill = "Koefisien x1") +
theme_minimal()
# 26. INTERPOLASI SPASIAL
# Nearest Neighbour
cent <- st_centroid(Jatim_sf_sub)
coords <- st_coordinates(cent)
pts <- Jatim_sf_sub %>%
st_drop_geometry() %>%
transmute(
id = NAME_2,
x = coords[,1],
y = coords[,2],
Z = Y
)
s0 <- data.frame(x = mean(pts$x), y = mean(pts$y))
edist <- function(x1, y1, x2, y2) sqrt((x1 - x2)^2 + (y1 - y2)^2)
pts$dist_s0 <- with(pts, edist(x, y, s0$x, s0$y))
pts_ord <- pts[order(pts$dist_s0), ]
nearest <- pts_ord[1, ]
cat(sprintf("Nearest to s0=(%.4f,%.4f) is %s at distance %.6f with Z=%g\n",
s0$x, s0$y, nearest$id, nearest$dist_s0, nearest$Z))## Nearest to s0=(112.5488,-7.6657) is Mojokerto at distance 0.132194 with Z=0.161## [1] 0.161if (requireNamespace("deldir", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE)) {
library(deldir)
library(ggplot2)
xlim <- range(pts$x) + c(-0.05, 0.05) * diff(range(pts$x))
ylim <- range(pts$y) + c(-0.05, 0.05) * diff(range(pts$y))
dd <- deldir::deldir(pts$x, pts$y, rw = c(xlim, ylim))
tiles <- deldir::tile.list(dd)
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(
id = pts$id[i],
x = tiles[[i]]$x,
y = tiles[[i]]$y
)
}))
ggplot() +
geom_polygon(data = poly_df, aes(x = x, y = y, group = id),
fill = NA, color = "grey40") +
geom_point(data = pts, aes(x = x, y = y), size = 2.5) +
geom_point(data = s0, aes(x = x, y = y), shape = 4, size = 3, stroke = 1.2) +
coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
labs(title = "Voronoi / Nearest Neighbor (berbasis centroid kab/kota)",
subtitle = "Prediksi NN: nilai di s0 = nilai titik terdekat",
x = "x", y = "y") +
theme_minimal(base_size = 12)
} else {
message("Lewati plot: paket 'deldir' dan/atau 'ggplot2' belum terpasang.")
}
# Inverse Distance Weighting
idw_predict <- function(x, y, z, x0, y0, p = 2, k = NULL) {
d <- sqrt((x - x0)^2 + (y - y0)^2)
if (any(d == 0)) return(z[which.min(d)])
if (!is.null(k) && k < length(d)) {
idx <- order(d)[1:k]
d <- d[idx]; z <- z[idx]
}
w <- 1 / (d^p)
sum(w * z) / sum(w)
}
Jatim_p <- st_transform(Jatim_sf_sub, 32749)
jatim_union <- st_union(st_geometry(Jatim_p))
grd_sf <- st_make_grid(
jatim_union,
cellsize = 5000,
what = "centers"
)
Xgrid <- st_coordinates(grd_sf)
dim(Xgrid)## [1] 9877 2df_idw <- data.frame(pts)
x0 <- Xgrid[10, 1]
y0 <- Xgrid[10, 2]
idw_p2 <- idw_predict(df_idw$x, df_idw$y, df_idw$Z, x0, y0, p = 2)
idw_p1 <- idw_predict(df_idw$x, df_idw$y, df_idw$Z, x0, y0, p = 1)
c(p2 = idw_p2, p1 = idw_p1)## p2 p1
## 0.4426579 0.4426579bb <- st_bbox(Jatim_sf)
xg <- seq(bb["xmin"], bb["xmax"], length.out=200)
yg <- seq(bb["ymin"], bb["ymax"], length.out=200)
G <- expand.grid(x = xg, y = yg)
G$z_p2 <- mapply(function(x0, y0) idw_predict(df_idw$x, df_idw$y, df_idw$Z, x0, y0, p = 2), G$x, G$y)
G$z_p1 <- mapply(function(x0, y0) idw_predict(df_idw$x, df_idw$y, df_idw$Z, x0, y0, p = 1), G$x, G$y)
xlim <- range(G$x, na.rm = TRUE)
ylim <- range(G$y, na.rm = TRUE)
df_lbl <- df_idw |>
dplyr::filter(Z > quantile(Z, 0.8))
plt1 <- ggplot(G, aes(x, y, fill = z_p1)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p1), color = "white", alpha = 0.6) +
geom_point(data = df_idw, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text_repel(
data = df_lbl,
aes(x, y, label = id),
inherit.aes = FALSE,
size = 2.5,
color = "white",
box.padding = 0.3,
point.padding = 0.2,
segment.color = "grey60"
) +
coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
labs(title = "Permukaan IDW (p = 1)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
plt2 <- ggplot(G, aes(x, y, fill = z_p2)) +
geom_raster(interpolate = TRUE) +
geom_contour(aes(z = z_p1), color = "white", alpha = 0.6) +
geom_point(data = df_idw, aes(x, y), inherit.aes = FALSE, size = 2) +
geom_text_repel(
data = df_lbl,
aes(x, y, label = id),
inherit.aes = FALSE,
size = 2.5,
color = "white",
box.padding = 0.3,
point.padding = 0.2,
segment.color = "grey60"
) +
coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
labs(title = "Permukaan IDW (p = 2)", x = "x", y = "y", fill = "Ẑ") +
theme_minimal(base_size = 12)
plt1 ; plt2
# 27. KRIGING
pts_sf <- st_as_sf(pts, coords = c("x","y"), crs = 4326)
ggplot(pts_sf) +
geom_sf(aes(color = Z), size = 2) +
scale_color_viridis_c() +
theme_minimal() +
labs(title = "Lokasi Sampling Kusta Jawa Timur", color = "Kusta")
# Variogram empiris
st_crs(pts_sf) <- 4326
pts_utm_sf <- st_transform(pts_sf, 32749)
pts_utm <- as(pts_utm_sf, "Spatial")
dmax <- max(spDists(coordinates(pts_utm), longlat = FALSE))
vgm_emp <- variogram(Z ~ 1, pts_utm,
cutoff = 0.6 * dmax,
width = (0.6 * dmax) / 12)
plot(vgm_emp, main = "Semivariogram Empiris (UTM)")
sill0 <- var(pts_utm$Z, na.rm = TRUE)
r0 <- max(vgm_emp$dist, na.rm = TRUE) / 3
fit_gau <- fit.variogram(vgm_emp, vgm(psill = sill0, model="Gau", range=r0, nugget=0))
vgm_mod <- fit_gau
plot(vgm_emp, vgm_mod, main = "Semivariogram Empiris + Model")
# Grid prediksi (meter)
bb <- bbox(pts_utm)
res <- 5000
grd <- expand.grid(
x = seq(bb[1,1], bb[1,2], by = res),
y = seq(bb[2,1], bb[2,2], by = res)
)
coordinates(grd) <- ~ x + y
proj4string(grd) <- proj4string(pts_utm)
ok_result <- krige(Z ~ 1, locations = pts_utm, newdata = grd, model = vgm_mod)## [using ordinary kriging]ok_sf <- st_as_sf(ok_result)
ggplot(ok_sf) +
geom_sf(aes(color = var1.pred), size = 1) +
scale_color_viridis_c(name = "Kusta pred.") +
theme_minimal() +
labs(title = "Prediksi Kusta dengan Ordinary Kriging")
ggplot(ok_sf) +
geom_sf(aes(color = var1.var), size = 1) +
scale_color_viridis_c(name = "Kusta pred.") +
theme_minimal() +
labs(title = "Varian Kriging (Ketidakpastian Prediksi Kusta")
set.seed(123)
cv_ok <- krige.cv(
formula = Z ~ 1,
locations = pts_utm_sf,
model = vgm_mod,
nfold = nrow(pts_utm_sf)
)
head(cv_ok)## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 589401.4 ymin: 9073936 xmax: 853090.1 ymax: 9220961
## Projected CRS: WGS 84 / UTM zone 49S
## var1.pred var1.var observed residual zscore fold
## 1 1.01773918 0.02218522 1.113 0.095260824 0.63956165 1
## 2 0.30765952 0.06077300 0.187 -0.120659522 -0.48944768 2
## 3 0.13416840 0.01920750 0.000 -0.134168403 -0.96808786 3
## 4 0.07689338 0.02147220 0.082 0.005106622 0.03484942 4
## 5 0.30772670 0.02142253 0.430 0.122273304 0.83540351 5
## 6 0.84662363 0.02399188 0.386 -0.460623628 -2.97381585 6
## geometry
## 1 POINT (713211.8 9220961)
## 2 POINT (853090.1 9073936)
## 3 POINT (668815.7 9133845)
## 4 POINT (636408.8 9101112)
## 5 POINT (589401.4 9197933)
## 6 POINT (824998.7 9120734)## [1] 0.2338848ggplot(as.data.frame(cv_ok), aes(x = observed, y = var1.pred)) +
geom_point(alpha = 0.7) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey40") +
theme_minimal() +
labs(
title = paste0("Cross-validation Ordinary Kriging (RMSE = ", round(rmse_ok, 2), ")"),
x = "Observasi Kusta",
y = "Prediksi Kusta (LOO)"
)
# 28. UNIVERSAL KRIGING
gridded(grd) <- TRUE
pts_utm@data$y <- sp::coordinates(pts_utm)[,2]
uk <- krige(log1p(Z) ~ x + y, pts_utm, grd, model = vgm_mod)## [using universal kriging]uk_px <- as(uk, "SpatialPixelsDataFrame")
jatim_utm_sf <- st_transform(Jatim_sf_sub, 32749)
jatim_utm_sp <- as(st_union(st_geometry(jatim_utm_sf)), "Spatial")
uk_mask <- raster::mask(raster::raster(uk_px["var1.pred"]), jatim_utm_sp)
plot(uk_mask, main="Universal Kriging (x + y) — Prediksi log(Kusta)")
Interpretasi :
Hasil interpolasi menunjukkan pola sebaran prevalensi kusta yang konsisten dengan hasil analisis spasial sebelumnya, khususnya dalam mengidentifikasi wilayah dengan intensitas tinggi. Secara keseluruhan, hasil penelitian menunjukkan bahwa prevalensi penyakit kusta di Provinsi Jawa Timur dipengaruhi oleh faktor-faktor sanitasi lingkungan yang memiliki keterkaitan spasial. Oleh karena itu, kebijakan pengendalian kusta perlu mempertimbangkan pendekatan lintas wilayah dan berbasis spasial.
5 Kesimpulan dan Saran
5.1 Kesimpulan
Prevalensi penyakit kusta di Provinsi Jawa Timur menunjukkan pola spasial yang tidak acak. Hasil analisis deskriptif dan pemetaan menunjukkan adanya pengelompokkan wilayah dengan prevalensi kusta yang relatif tinggi dan rendah. Temuan ini mengindikasikan bahwa distribusi penyakit kusta antar kabupaten/kota di Jawa Timur memiliki variasi yang cukup besar dan tidak tersebar secara merata. Diperoleh bahwa prevalensi penyakit kusta tertinggi terjadi di Kabupaten Sumenep dengan nilai 2,06. Kemudian, kabupaten atau kota dengan prevalensi penyakit kusta terendah yaitu terdapat di Kota Blitar, Kota Mojokerto, dan Kota Batu dengan prevalensi penyakit kusta sebesar 0.
Model Spatial Autoregressive Combined (SAC) yang didapatkan ialah :
\[ y_i = 0.84143\sum_{j=1}^{n} w_{ij}y_j -0.03746499 +0.00424778\,X_{1i} +0.00030904\,X_{2i} +u_i \]
Dimana meskipun intersep dan variabel \(X_2\) tidak signifikan secara statistik, variabel \(X_1\) bersifat signifikan positif. Artinya, ketika \(X_1\) meningkat 1 satuan, maka prevalensi kusta (\(Y\)) cenderung naik sekitar 0.00425 satuan, dengan asumsi variabel lain konstan pada taraf signifikansi 10%.
Berdasarkan nilai parameter \(ρ\) dengan nilai 0,84 dan \(p-value < 2,22 {e^{-16}}\), menunjukkan adanya ketergantungan spasial yang kuat pada variabel respon. Hal ini mengindikasikan bahwa prevalensi kusta pada suatu kabupaten/kota dipengaruhi secara signifikan oleh prevalensi kusta pada kabupaten/kota yang bertetangga.
Berdasarkan nilai parameter \(λ\) dengan nilai 0,57 dan \(p-value < 0,00121\), menunjukkan adanya dependensi spasial pada komponen galat sehingga terdapat faktor faktor lain di luar variabel model yang berpola spasial dan turut memengaruhi prevalensi kusta.
Berdasarkan kriteria AIC, model SAC (AIC=19,759) jauh lebih baik dibanding model OLS (AIC=50,102), sehingga model SAC dipilih sebagai model akhir dalam penelitian ini.
Pemerintah daerah dan instansi terkait disarankan untuk mempertimbangkan penggunaan model ekonometrika spasial dalam perencanaan dan evaluasi kebijakan kesehatan. Model seperti SAC mampu memberikan gambaran yang lebih komprehensif dibandingkan model regresi klasik karena mempertimbangkan keterkaitan antarwilayah. Selain itu, perlu adanya penyesuaian program berbasis kondisi lokal untuk mampu meningkatkan efektivitas intervensi kesehatan dalam upaya penanganan penyakit kusta.
5.2 Saran
Penelitian selanjutnya disarankan untuk menggunakan unit spasial yang lebih kecil, seperti kecamatan atau desa, agar variasi spasial prevalensi kusta dapat dianalisis secara lebih detail dan mendukung perencanaan intervensi yang lebih tepat sasaran.
6 Lampiran
Berikut adalah tautan dari analisis menggunakan Dashboard yang dapat diakses beserta dengan lampiran syntax :
https://fazilaazraanggina.shinyapps.io/SPASIAL/
options(shiny.maxRequestSize = 200 * 1024^2)
library(shiny)
library(shinydashboard)
library(shinycssloaders)
library(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)
library(openxlsx)
library(spatialreg)
library(Matrix)
library(spgwr)
library(lmtest)
library(car)
library(gstat)
library(FNN)
library(stars)
library(ggrepel)
library(viridis)
library(raster)
library(DT)
library(leaflet)
library(plotly)
library(stringr)
norm_name <- function(x){
x %>%
as.character() %>%
stringr::str_trim() %>%
stringr::str_to_upper() %>%
stringr::str_replace_all("\\s+", " ") %>%
stringr::str_replace_all("[^A-Z0-9\\s/\\-\\.]", "")
}
read_map_any <- function(path){
ext <- tolower(tools::file_ext(path))
if (ext == "rds") return(sf::st_as_sf(readRDS(path)))
sf::st_read(path, quiet = TRUE)
}
pal_bupk <- function(n = 256){
grDevices::colorRampPalette(c("#0B1F4B", "#7C3AED", "#EC4899"))(n)
}
pal_div_bupk <- function(n = 256){
grDevices::colorRampPalette(c("#0B1F4B", "#F9FAFB", "#EC4899"))(n)
}
ui <- dashboardPage(
skin = "purple",
dashboardHeader(title = span("UAS Spasial", style="font-weight:900;")),
dashboardSidebar(
width = 300,
sidebarMenu(
id = "tabs",
menuItem("Overview", tabName="overview", icon=icon("gauge-high")),
menuItem("Upload", tabName="upload", icon=icon("upload")),
menuItem("Peta Kusta", tabName="map_kusta", icon=icon("map")),
menuItem("Autokorelasi", tabName="auto", icon=icon("project-diagram")),
menuItem("Getis-Ord", tabName="getis", icon=icon("fire")),
menuItem("Model Spasial", tabName="models", icon=icon("brain")),
menuItem("Interpolasi", tabName="interp", icon=icon("layer-group")),
menuItem("Export", tabName="export", icon=icon("file-export")),
menuItem("About", tabName="about", icon=icon("circle-info"))
)
),
dashboardBody(
tags$head(
tags$style(HTML("
/* ===== Dark UI tapi card putih, teks gelap ===== */
body, .content-wrapper, .right-side { background:#070A1A !important; }
.skin-purple .main-header .logo,
.skin-purple .main-header .navbar { background:#0B0F2A !important; }
.skin-purple .main-sidebar { background:#070A1A !important; }
.sidebar a { color:#E8EAF6 !important; }
.sidebar-menu>li>a:hover { background:#0E1334 !important; }
.skin-purple .sidebar-menu>li.active>a { border-left-color:#EC4899 !important; }
/* card putih */
.box { border-radius:16px !important; border:1px solid #1B2460 !important; background:#F9FAFB !important; }
.box-title, h1,h2,h3,h4,h5,label, .control-label { color:#111827 !important; font-weight:900 !important; }
.box-body, .help-block, .shiny-text-output, pre { color:#111827 !important; }
/* ValueBox number putih */
.small-box h3 { color:#FFFFFF !important; font-weight:900 !important; }
.small-box p { color:#FFFFFF !important; font-weight:800 !important; }
.small-box { border-radius:16px !important; }
.leaflet-container { border-radius:16px !important; }
table.dataTable { color:#111827 !important; }
"))
),
tabItems(
# ============================================================
# OVERVIEW
# ============================================================
tabItem(tabName="overview",
fluidRow(
valueBoxOutput("vb_n", width=3),
valueBoxOutput("vb_mean", width=3),
valueBoxOutput("vb_moran", width=3),
valueBoxOutput("vb_best", width=3)
),
fluidRow(
box(
title="Distribusi Variabel",
width=12, status="primary",
selectInput("var_hist", "Pilih variabel", choices=c(
"Y"="Y", "X1"="X1", "X2"="X2"
), selected="kusta"),
plotlyOutput("plt_hist", height=380) %>% withSpinner()
)
),
fluidRow(
box(
title="Peta Variabel",
width=12, status="info",
selectInput("var_map", "Pilih variabel peta", choices=c(
"Y"="Y", "X1"="X1", "X2"="X2"
), selected="kusta"),
leafletOutput("map_var", height=560) %>% withSpinner()
)
)
),
# ============================================================
# UPLOAD (tanpa progress join text / cek join)
# ============================================================
tabItem(tabName="upload",
box(
title="Upload Data",
width=12, status="primary",
fileInput("excel", "Upload Excel (.xlsx)", accept=c(".xlsx")),
fileInput("peta", "Upload Peta (RDS)", accept=c(".rds",".shp",".gpkg",".geojson",".zip")),
actionButton("load", "Muat & Gabungkan", class="btn btn-primary btn-lg")
),
box(
title="Preview Excel",
width=12, status="info",
DTOutput("tbl_excel_raw") %>% withSpinner()
)
),
# ============================================================
# PETA KUSTA (Kontinu + Diskrit)
# ============================================================
tabItem(tabName="map_kusta",
fluidRow(
box(title="Peta Prevalensi Penyakit Kusta di Provinsi Jawa Timur (Kontinu)", width=6, status="primary",
plotOutput("gg_kusta_cont", height=450) %>% withSpinner()),
box(title="Peta Prevalensi Penyakit Kusta di Provinsi Jawa Timur (Diskrit / Kategori)", width=6, status="warning",
plotOutput("gg_kusta_disc", height=450) %>% withSpinner())
),
box(title="Breaks Diskrit", width=12, status="info",
numericInput("b1", "Break 1", value=0.11075),
numericInput("b2", "Break 2", value=0.22600),
numericInput("b3", "Break 3", value=0.56175),
numericInput("b4", "Break 4", value=2.06300)
)
),
# ============================================================
# AUTOKORELASI
# ============================================================
tabItem(tabName="auto",
fluidRow(
box(title="Setting Bobot", width=4, status="info",
selectInput("queen_rook", "Rook/Queen", choices=c("Rook (queen=FALSE)"="rook","Queen (queen=TRUE)"="queen"), selected="rook"),
selectInput("w_style", "Style", choices=c("W","B","C","U","S"), selected="W"),
checkboxInput("zero_policy", "zero.policy=TRUE", value=TRUE),
numericInput("snap", "snap (poly2nb)", value=1e-6, min=0, step=1e-6)
),
box(title="Moran's I (Global)", width=8, status="primary",
verbatimTextOutput("txt_moran") %>% withSpinner())
),
fluidRow(
box(title="Local Moran (LISA)", width=6, status="info",
DTOutput("tbl_lisa") %>% withSpinner()),
box(title="LISA Cluster Map", width=6, status="warning",
leafletOutput("map_lisa", height=420) %>% withSpinner())
),
fluidRow(
box(title="Geary's C", width=12, status="primary",
verbatimTextOutput("txt_geary") %>% withSpinner())
)
),
# ============================================================
# GETIS-ORD
# ============================================================
tabItem(tabName="getis",
fluidRow(
box(title="Getis–Ord Gi* (z-score) + Kategori", width=6, status="primary",
verbatimTextOutput("txt_getis") %>% withSpinner()),
box(title="Peta Hot/Cold Spots", width=6, status="warning",
leafletOutput("map_hotcold", height=420) %>% withSpinner())
),
fluidRow(
box(title="Peta Raw G*", width=12, status="info",
leafletOutput("map_gstar", height=520) %>% withSpinner())
)
),
# ============================================================
# MODELS
# ============================================================
tabItem(tabName="models",
fluidRow(
box(title="OLS + Uji (Moran residual + LM tests + VIF)", width=12, status="warning",
verbatimTextOutput("txt_ols_tests") %>% withSpinner())
),
fluidRow(
box(title="Perbandingan Model (AIC)", width=12, status="primary",
DTOutput("tbl_models") %>% withSpinner())
),
fluidRow(
box(title="Ringkasan Model Terbaik", width=12, status="success",
verbatimTextOutput("txt_best") %>% withSpinner())
),
fluidRow(
box(title="Peta Prediksi", width=6, status="primary",
leafletOutput("map_pred", height=420) %>% withSpinner()),
box(title="Peta Residual", width=6, status="info",
leafletOutput("map_resid", height=420) %>% withSpinner())
)
),
# ============================================================
# INTERPOLASI
# ============================================================
tabItem(tabName="interp",
fluidRow(
box(title="Setting Interpolasi", width=4, status="primary",
selectInput("z_interp", "Variabel Z untuk interpolasi", choices=c(
"Kusta"="kusta", "Y"="Y", "X1"="X1", "X2"="X2"
), selected="kusta"),
sliderInput("grid_res_m", "Resolusi grid (meter)", min=2000, max=15000, value=5000, step=500),
sliderInput("idp", "IDW power (idp)", min=0.5, max=4, value=2, step=0.1),
actionButton("run_interp", "Jalankan Interpolasi", class="btn btn-primary")
),
box(title="Nearest Neighbor (Centroid)", width=8, status="info",
verbatimTextOutput("txt_nn") %>% withSpinner()
)
),
fluidRow(
box(title="Voronoi", width=12, status="warning",
plotOutput("plt_voronoi", height=460) %>% withSpinner()
)
),
fluidRow(
box(title="Peta IDW (surface)", width=12, status="primary",
leafletOutput("map_idw", height=560) %>% withSpinner()
)
),
fluidRow(
box(title="Kriging", width=12, status="info",
verbatimTextOutput("txt_kriging") %>% withSpinner()
)
),
fluidRow(
box(title="Peta Kriging (Prediksi)", width=6, status="primary",
leafletOutput("map_ok_pred", height=460) %>% withSpinner()),
box(title="Peta Kriging (Variansi)", width=6, status="warning",
leafletOutput("map_ok_var", height=460) %>% withSpinner())
)
),
# ============================================================
# EXPORT
# ============================================================
tabItem(tabName="export",
box(title="Download", width=12, status="primary",
downloadButton("dl_excel_raw", "Download Excel (CSV)"),
downloadButton("dl_clean", "Download Data Bersih Gabungan (CSV)"),
downloadButton("dl_lisa", "Download LISA (CSV)"),
downloadButton("dl_models", "Download Perbandingan Model (CSV)")
)
),
# ============================================================
# ABOUT
# ============================================================
tabItem(tabName="about",
box(title="Tentang", width=12, status="primary",
HTML("<div style='color:#111827'>
<b>Dashboard Projek UAS Spatial Statistics</b><br/>
Oleh : Fazila Azra Anggina (140610230039).<br/>
</div>"))
)
)
)
)
# ============================================================
# SERVER
# ============================================================
server <- function(input, output, session){
best_label <- reactiveVal("—")
# ============================================================
# EXCEL RAW (preview seperti excel)
# ============================================================
excel_raw <- reactive({
req(input$excel)
openxlsx::read.xlsx(input$excel$datapath)
})
output$tbl_excel_raw <- renderDT({
req(excel_raw())
datatable(excel_raw(), options=list(pageLength=12, scrollX=TRUE))
})
# ============================================================
# LOAD & JOIN (tidak ditampilkan progresnya di UI upload)
# ============================================================
jatim_sf_full <- eventReactive(input$load, {
req(input$excel, input$peta)
withProgress(message="Memuat data...", value=0, {
incProgress(0.2, detail="Baca Excel...")
kusta <- openxlsx::read.xlsx(input$excel$datapath)
validate(need(any(grepl("Kabupaten", names(kusta), ignore.case=TRUE)), "Kolom Kabupaten/Kota tidak ditemukan di Excel."))
validate(need(any(grepl("Prevalensi", names(kusta), ignore.case=TRUE)), "Kolom Prevalensi.Kusta tidak ditemukan di Excel."))
kab_col <- names(kusta)[grepl("Kabupaten", names(kusta), ignore.case=TRUE)][1]
prev_col <- names(kusta)[grepl("Prevalensi", names(kusta), ignore.case=TRUE)][1]
kusta_data <- data.frame(
wilayah = kusta[[kab_col]],
kusta = suppressWarnings(as.numeric(kusta[[prev_col]]))
) %>% mutate(wilayah_key = norm_name(wilayah))
incProgress(0.5, detail="Baca peta...")
Indo_Kab <- read_map_any(input$peta$datapath)
incProgress(0.7, detail="Filter Jawa Timur & Join...")
nm <- names(Indo_Kab)
validate(need(any(nm=="NAME_1"), "Peta tidak punya kolom NAME_1 (provinsi). Pastikan peta GADM level 2."))
validate(need(any(nm=="NAME_2"), "Peta tidak punya kolom NAME_2 (kab/kota). Pastikan peta GADM level 2."))
Jatim <- Indo_Kab[Indo_Kab$NAME_1 == "Jawa Timur", ]
Jatim_sf <- st_as_sf(Jatim) %>%
mutate(wilayah_key = norm_name(NAME_2)) %>%
left_join(kusta_data, by = "wilayah_key")
incProgress(0.9, detail="Muat Y, X1, X2 (jika ada)...")
if(ncol(kusta) >= 4){
data_model <- kusta
names(data_model)[which(names(data_model)==kab_col)] <- "Kabupaten_Kota"
vars_model <- names(data_model)[2:4]
tmp <- data_model %>%
transmute(
Kabupaten_Kota = .data$Kabupaten_Kota,
Y = suppressWarnings(as.numeric(.data[[vars_model[1]]])),
X1 = suppressWarnings(as.numeric(.data[[vars_model[2]]])),
X2 = suppressWarnings(as.numeric(.data[[vars_model[3]]])),
wilayah_key = norm_name(.data$Kabupaten_Kota)
)
Jatim_sf <- Jatim_sf %>% left_join(tmp %>% select(wilayah_key, Y, X1, X2), by="wilayah_key")
} else {
Jatim_sf$Y <- Jatim_sf$kusta
Jatim_sf$X1 <- NA_real_
Jatim_sf$X2 <- NA_real_
}
incProgress(1.0, detail="Selesai")
Jatim_sf
})
})
# bersih kusta
jatim_sf_kusta <- reactive({
req(jatim_sf_full())
jatim_sf_full() %>% filter(!is.na(kusta))
})
# bersih model
jatim_sf_model <- reactive({
req(jatim_sf_full())
jatim_sf_full() %>% filter(complete.cases(Y, X1, X2))
})
# ============================================================
# ValueBoxes
# ============================================================
output$vb_n <- renderValueBox({
req(jatim_sf_kusta())
valueBox(nrow(jatim_sf_kusta()), "Jumlah Kab/Kota", color="purple", icon=icon("map"))
})
output$vb_mean <- renderValueBox({
req(jatim_sf_kusta())
valueBox(round(mean(jatim_sf_kusta()$kusta, na.rm=TRUE), 4), "Rata-Rata Prevalensi", color="fuchsia", icon=icon("chart-line"))
})
# ============================================================
# Weights (kusta)
# ============================================================
lw_kusta <- reactive({
req(jatim_sf_kusta())
sf <- jatim_sf_kusta()
queen <- (input$queen_rook == "queen")
nb <- spdep::poly2nb(sf::as_Spatial(sf), queen = queen, snap = input$snap)
spdep::nb2listw(nb, style = input$w_style, zero.policy = input$zero_policy)
})
output$vb_moran <- renderValueBox({
req(jatim_sf_kusta(), lw_kusta())
mt <- moran.test(jatim_sf_kusta()$kusta, lw_kusta(), zero.policy=input$zero_policy, na.action=na.exclude)
valueBox(round(unname(mt$estimate[["Moran I statistic"]]), 4), "Hasil Moran's I", color="navy", icon=icon("project-diagram"))
})
output$vb_best <- renderValueBox({
valueBox(best_label(), "Model Terbaik", color="purple", icon=icon("check"))
})
# ============================================================
# Overview Histogram (Navy)
# ============================================================
output$plt_hist <- renderPlotly({
req(jatim_sf_full(), input$var_hist)
df <- jatim_sf_full() %>% st_drop_geometry()
v <- input$var_hist
validate(need(v %in% names(df), "Variabel tidak ditemukan."))
plot_ly(
df,
x = df[[v]],
type = "histogram",
nbinsx = 20,
marker = list(color = "#0B1F4B")
) %>%
layout(
xaxis = list(title=v),
yaxis = list(title="Frekuensi"),
margin = list(l=50,r=10,t=20,b=50)
)
})
# ============================================================
# Overview Map (Theme Light)
# ============================================================
output$map_var <- renderLeaflet({
req(jatim_sf_full(), input$var_map)
sf <- jatim_sf_full()
var <- input$var_map
validate(need(var %in% names(sf), "Variabel peta tidak ditemukan."))
pal <- colorNumeric(palette = pal_bupk(), domain = sf[[var]], na.color="#D1D5DB")
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>% # LIGHT tiles
addPolygons(
fillColor = pal(sf[[var]]),
fillOpacity = 0.85,
color = "#111827",
weight = 1,
popup = paste0("<b>", sf$NAME_2, "</b><br>", var, ": ", sf[[var]])
) %>%
addLegend("bottomright", pal = pal, values = sf[[var]], title = var)
})
# ============================================================
# Peta kusta (ggplot)
# ============================================================
output$gg_kusta_cont <- renderPlot({
req(jatim_sf_kusta())
ggplot() +
geom_sf(data=jatim_sf_kusta(), aes(fill = kusta), color=NA) +
theme_bw() +
scale_fill_gradient(low = "pink", high = "black") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Peta Kasus Kusta Jawa Timur (Kontinu)", fill = "Kusta")
})
output$gg_kusta_disc <- renderPlot({
req(jatim_sf_kusta())
sf <- jatim_sf_kusta()
breaks <- c(0, input$b1, input$b2, input$b3, input$b4)
labels <- c("Very Low","Low","High","Very High")
sf$kusta_Discrete <- cut(sf$kusta, breaks=breaks, labels=labels, right=TRUE, include.lowest=TRUE)
ggplot() +
geom_sf(data=sf, aes(fill = kusta_Discrete), color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low"="pink","Low"="magenta","High"="purple","Very High"="black")) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Kategori Kasus Kusta Jawa Timur (Diskrit)", fill = "Kategori")
})
# ============================================================
# MORAN / LISA / GEARY
# ============================================================
output$txt_moran <- renderPrint({
req(jatim_sf_kusta(), lw_kusta())
moran.test(jatim_sf_kusta()$kusta, lw_kusta(), zero.policy=input$zero_policy, na.action=na.exclude)
})
lisa_tbl <- reactive({
req(jatim_sf_kusta(), lw_kusta())
sf <- jatim_sf_kusta()
lw <- lw_kusta()
x <- as.numeric(scale(sf$kusta))
lagx <- spdep::lag.listw(lw, x, zero.policy=input$zero_policy)
lisa <- spdep::localmoran(x, lw, alternative="two.sided", zero.policy=input$zero_policy, na.action=na.exclude)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
x >= 0 & lagx >= 0 ~ "High-High",
x < 0 & lagx < 0 ~ "Low-Low",
x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)
out <- dplyr::bind_cols(sf %>% st_drop_geometry() %>% select(NAME_2, kusta), lisa_df) %>%
mutate(quad = ifelse(`Pi.two.sided` <= alpha, quad, "Not significant"))
out
})
output$tbl_lisa <- renderDT({
req(lisa_tbl())
datatable(lisa_tbl(), options=list(pageLength=10, scrollX=TRUE))
})
output$map_lisa <- renderLeaflet({
req(jatim_sf_kusta(), lisa_tbl())
sf <- jatim_sf_kusta()
lt <- lisa_tbl()
sf$quad <- factor(
lt$quad,
levels=c("High-High","Low-Low","High-Low (Outlier)","Low-High (Outlier)","Not significant")
)
pal <- colorFactor(
palette=c(
"High-High"="#EC4899",
"Low-Low"="#C026D3",
"High-Low (Outlier)"="#7C3AED",
"Low-High (Outlier)"="#0B1F4B",
"Not significant"="#9CA3AF"
),
domain=levels(sf$quad)
)
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>% # LIGHT tiles
addPolygons(
fillColor = pal(as.character(sf$quad)),
fillOpacity = 0.85,
color = "#111827",
weight = 1,
popup = paste0("<b>", sf$NAME_2, "</b><br>Cluster: ", sf$quad)
) %>%
addLegend("bottomright", pal=pal, values=sf$quad, title="LISA (p<=0.05)")
})
output$txt_geary <- renderPrint({
req(jatim_sf_kusta(), lw_kusta())
geary.test(jatim_sf_kusta()$kusta, lw_kusta(), randomisation=TRUE, alternative="two.sided",
zero.policy=input$zero_policy, na.action=na.exclude)
})
# ============================================================
# GETIS-ORD
# ============================================================
getis_pack <- reactive({
req(jatim_sf_kusta(), lw_kusta())
sf <- jatim_sf_kusta()
lw <- lw_kusta()
nb <- lw$neighbours
WB <- nb2mat(nb, style="B", zero.policy = TRUE)
x_raw <- sf$kusta
sum_x <- sum(x_raw, na.rm=TRUE)
num_G <- as.numeric(WB %*% x_raw)
den_G <- (sum_x - x_raw)
G_raw <- num_G / den_G
WB_star <- WB; diag(WB_star) <- 1
num_Gs <- as.numeric(WB_star %*% x_raw)
den_Gs <- sum_x
G_star_raw <- num_Gs / den_Gs
Gz <- spdep::localG(x_raw, listw = lw)
sf %>%
mutate(
G_raw = G_raw,
G_star_raw = G_star_raw,
z_Gistar = as.numeric(Gz),
hotcold = dplyr::case_when(
z_Gistar >= 1.96 ~ "Hot spot (p≈0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p≈0.05)",
TRUE ~ "Not significant"
)
)
})
output$txt_getis <- renderPrint({
req(getis_pack())
df <- getis_pack() %>% st_drop_geometry()
cat("Ringkas z_Gistar:\n"); print(summary(df$z_Gistar))
cat("\nJumlah kategori:\n"); print(table(df$hotcold))
})
output$map_hotcold <- renderLeaflet({
req(getis_pack())
sf <- getis_pack()
pal <- colorFactor(
palette=c(
"Hot spot (p≈0.05)"="#111827",
"Cold spot (p≈0.05)"="#EC4899",
"Not significant"="#9CA3AF"
),
domain=unique(sf$hotcold)
)
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(
fillColor = pal(sf$hotcold),
fillOpacity = 0.85,
color = "#111827",
weight = 1,
popup = paste0("<b>", sf$NAME_2, "</b><br>z_Gi*: ", round(sf$z_Gistar,3), "<br>", sf$hotcold)
) %>%
addLegend("bottomright", pal=pal, values=sf$hotcold, title="Gi*")
})
output$map_gstar <- renderLeaflet({
req(getis_pack())
sf <- getis_pack()
pal <- colorNumeric(palette = pal_bupk(), domain = sf$G_star_raw, na.color="#D1D5DB")
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(
fillColor = pal(sf$G_star_raw),
fillOpacity = 0.85,
color = "#111827",
weight = 1,
popup = paste0("<b>", sf$NAME_2, "</b><br>G*_raw: ", round(sf$G_star_raw,4))
) %>%
addLegend("bottomright", pal=pal, values=sf$G_star_raw, title="G*_raw")
})
# ============================================================
# MODELS
# ============================================================
lw_model <- reactive({
req(jatim_sf_model())
sf <- jatim_sf_model()
nb <- poly2nb(as_Spatial(sf), queen = TRUE)
nb2listw(nb, style = "W", zero.policy = TRUE)
})
model_pack <- reactive({
req(jatim_sf_model(), lw_model())
sf <- jatim_sf_model()
lw <- lw_model()
withProgress(message="Menjalankan model...", value=0, {
incProgress(0.1, detail="OLS")
ols <- lm(Y ~ X1 + X2, data = sf)
incProgress(0.2, detail="OLS tests")
moran_res <- lm.morantest(ols, listw = lw)
lmtests <- lm.LMtests(ols, listw = lw, test=c("LMlag","LMerr","RLMlag","RLMerr"))
vifv <- car::vif(ols)
incProgress(0.35, detail="SLX")
slx <- spatialreg::lmSLX(Y ~ X1 + X2, data = sf, listw = lw, Durbin = TRUE, zero.policy = TRUE)
incProgress(0.50, detail="SAR")
sar <- lagsarlm(Y ~ X1 + X2, data = sf, listw = lw, method="eigen", zero.policy=TRUE)
incProgress(0.60, detail="SEM")
sem <- errorsarlm(Y ~ X1 + X2, data = sf, listw = lw, method="eigen", zero.policy=TRUE)
incProgress(0.70, detail="SDM")
sdm <- lagsarlm(Y ~ X1 + X2, data = sf, listw = lw, method="eigen", zero.policy=TRUE)
incProgress(0.80, detail="SDEM")
sdem <- errorsarlm(Y ~ X1 + X2, data = sf, listw = lw, Durbin=TRUE, method="eigen", zero.policy=TRUE)
incProgress(0.90, detail="SAC & GNS")
sac <- sacsarlm(Y ~ X1 + X2, data = sf, listw = lw, method="eigen", zero.policy=TRUE)
gns <- sacsarlm(Y ~ X1 + X2, data = sf, listw = lw, Durbin=TRUE, method="eigen", zero.policy=TRUE)
cmp <- data.frame(
Model = c("OLS","SLX","SAR","SEM","SDM","SDEM","SAC","GNS"),
AIC = c(AIC(ols), AIC(slx), AIC(sar), AIC(sem), AIC(sdm), AIC(sdem), AIC(sac), AIC(gns))
) %>% arrange(AIC)
best_name <- cmp$Model[1]
best_label(best_name)
best <- switch(best_name,
"OLS" = ols, "SLX" = slx, "SAR" = sar, "SEM" = sem, "SDM" = sdm, "SDEM" = sdem, "SAC" = sac, "GNS" = gns
)
incProgress(1.0, detail="Selesai")
list(ols=ols, slx=slx, sar=sar, sem=sem, sdm=sdm, sdem=sdem, sac=sac, gns=gns,
cmp=cmp, best=best, best_name=best_name,
moran_res=moran_res, lmtests=lmtests, vif=vifv)
})
})
output$txt_ols_tests <- renderPrint({
req(model_pack())
cat("OLS summary:\n"); print(summary(model_pack()$ols))
cat("\nMoran residual (OLS):\n"); print(model_pack()$moran_res)
cat("\nLM tests:\n"); print(model_pack()$lmtests)
cat("\nVIF:\n"); print(model_pack()$vif)
})
output$tbl_models <- renderDT({
req(model_pack())
datatable(model_pack()$cmp, options=list(pageLength=10))
})
output$txt_best <- renderPrint({
req(model_pack())
cat("Best Model:", model_pack()$best_name, "\n\n")
print(summary(model_pack()$best))
})
sf_pred <- reactive({
req(jatim_sf_model(), model_pack())
sf <- jatim_sf_model()
bm <- model_pack()$best
pred <- if(inherits(bm, "lm")) predict(bm, newdata=sf) else as.numeric(bm$fitted.values)
resid <- if(inherits(bm, "lm")) residuals(bm) else as.numeric(bm$residuals)
sf$pred <- pred
sf$resid <- resid
sf
})
output$map_pred <- renderLeaflet({
req(sf_pred())
sf <- sf_pred()
pal <- colorNumeric(palette = pal_bupk(), domain = sf$pred, na.color="#D1D5DB")
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(
fillColor = pal(sf$pred),
fillOpacity = 0.85,
color="#111827", weight=1,
popup = paste0("<b>", sf$NAME_2, "</b><br>Pred: ", round(sf$pred,4))
) %>%
addLegend("bottomright", pal=pal, values=sf$pred, title="Prediksi")
})
output$map_resid <- renderLeaflet({
req(sf_pred())
sf <- sf_pred()
pal <- colorNumeric(palette = pal_div_bupk(), domain = sf$resid, na.color="#D1D5DB")
leaflet(sf) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(
fillColor = pal(sf$resid),
fillOpacity = 0.85,
color="#111827", weight=1,
popup = paste0("<b>", sf$NAME_2, "</b><br>Resid: ", round(sf$resid,4))
) %>%
addLegend("bottomright", pal=pal, values=sf$resid, title="Residual")
})
# ============================================================
# INTERPOLASI (NN + Voronoi + IDW + Kriging) - tetap
# ============================================================
interp_base_sf <- reactive({
req(jatim_sf_full(), input$z_interp)
sf <- jatim_sf_full()
z <- input$z_interp
validate(need(z %in% names(sf), "Variabel Z tidak ada."))
sf <- sf %>% filter(!is.na(.data[[z]]))
validate(need(nrow(sf) >= 5, "Data valid untuk interpolasi terlalu sedikit (butuh >= 5)."))
sf
})
interp_points <- reactive({
req(interp_base_sf(), input$z_interp)
sf <- interp_base_sf()
zname <- input$z_interp
sf_utm <- st_transform(sf, 32749)
cent <- st_centroid(sf_utm)
coords <- st_coordinates(cent)
sf_utm %>%
st_drop_geometry() %>%
transmute(
id = NAME_2,
x = coords[,1],
y = coords[,2],
Z = as.numeric(.data[[zname]])
)
})
output$txt_nn <- renderPrint({
req(interp_points())
pts <- interp_points()
s0 <- data.frame(x = mean(pts$x), y = mean(pts$y))
edist <- function(x1, y1, x2, y2) sqrt((x1 - x2)^2 + (y1 - y2)^2)
pts$dist_s0 <- with(pts, edist(x, y, s0$x, s0$y))
pts_ord <- pts[order(pts$dist_s0), ]
nearest <- pts_ord[1, ]
cat(sprintf("s0 = (%.2f, %.2f)\n", s0$x, s0$y))
cat(sprintf("Nearest: %s | distance: %.2f m | Zhat: %.6f\n", nearest$id, nearest$dist_s0, nearest$Z))
})
output$plt_voronoi <- renderPlot({
req(interp_points())
pts <- interp_points()
if (!requireNamespace("deldir", quietly = TRUE)) {
plot.new()
text(0.5, 0.5, "Paket 'deldir' belum terpasang.\nInstall: install.packages('deldir')", cex=1.2)
return()
}
library(deldir)
xlim <- range(pts$x) + c(-0.05, 0.05) * diff(range(pts$x))
ylim <- range(pts$y) + c(-0.05, 0.05) * diff(range(pts$y))
dd <- deldir::deldir(pts$x, pts$y, rw = c(xlim, ylim))
tiles <- deldir::tile.list(dd)
poly_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
data.frame(id = pts$id[i], x = tiles[[i]]$x, y = tiles[[i]]$y)
}))
s0 <- data.frame(x = mean(pts$x), y = mean(pts$y))
ggplot() +
geom_polygon(data = poly_df, aes(x=x, y=y, group=id), fill=NA, color="grey40") +
geom_point(data = pts, aes(x=x, y=y), size=2.5) +
geom_point(data = s0, aes(x=x, y=y), shape=4, size=4, stroke=1.2) +
coord_equal(xlim=xlim, ylim=ylim, expand=FALSE) +
labs(title="Voronoi / Nearest Neighbor", subtitle="X = s0") +
theme_minimal(base_size = 12)
})
interp_results <- eventReactive(input$run_interp, {
req(interp_base_sf(), interp_points(), input$grid_res_m, input$idp)
sf0 <- interp_base_sf()
pts <- interp_points()
withProgress(message="Interpolasi...", value=0, {
incProgress(0.15, detail="Siapkan mask UTM...")
sf_utm <- st_transform(sf0, 32749)
poly_union <- st_union(st_geometry(sf_utm))
poly_sp <- as(poly_union, "Spatial")
sp_pts <- pts
coordinates(sp_pts) <- ~ x + y
proj4string(sp_pts) <- CRS("+init=epsg:32749")
incProgress(0.30, detail="Buat grid (SpatialPixels)...")
bb <- st_bbox(st_as_sf(poly_union))
r <- raster::raster(
xmn=bb["xmin"], xmx=bb["xmax"],
ymn=bb["ymin"], ymx=bb["ymax"],
res = input$grid_res_m,
crs = sp::CRS("+init=epsg:32749")
)
grd <- as(r, "SpatialPixels")
incProgress(0.55, detail="Hitung IDW...")
idw_out <- gstat::idw(Z ~ 1, locations = sp_pts, newdata = grd, idp = input$idp)
idw_r <- raster::rasterFromXYZ(as.data.frame(idw_out)[,c("x","y","var1.pred")])
crs(idw_r) <- CRS("+init=epsg:32749")
idw_mask <- raster::mask(idw_r, poly_sp)
incProgress(0.75, detail="Kriging...")
ok_pred_mask <- NULL
ok_var_mask <- NULL
rmse_ok <- NA_real_
vgm_fit_txt <- NULL
try({
dmax <- max(spDists(coordinates(sp_pts), longlat = FALSE))
vgm_emp <- variogram(Z ~ 1, sp_pts, cutoff = 0.6*dmax, width = (0.6*dmax)/12)
sill0 <- var(sp_pts$Z, na.rm = TRUE)
r0 <- max(vgm_emp$dist, na.rm = TRUE) / 3
vgm_mod <- fit.variogram(vgm_emp, vgm(psill=sill0, model="Gau", range=r0, nugget=0))
ok <- krige(Z ~ 1, locations = sp_pts, newdata = grd, model = vgm_mod)
ok_pred <- raster::rasterFromXYZ(as.data.frame(ok)[,c("x","y","var1.pred")])
ok_var <- raster::rasterFromXYZ(as.data.frame(ok)[,c("x","y","var1.var")])
crs(ok_pred) <- CRS("+init=epsg:32749")
crs(ok_var) <- CRS("+init=epsg:32749")
ok_pred_mask <- raster::mask(ok_pred, poly_sp)
ok_var_mask <- raster::mask(ok_var, poly_sp)
cv_ok <- krige.cv(Z ~ 1, locations = sp_pts, model = vgm_mod, nfold = nrow(sp_pts))
rmse_ok <- sqrt(mean((cv_ok$observed - cv_ok$var1.pred)^2, na.rm=TRUE))
vgm_fit_txt <- capture.output(print(vgm_mod))
}, silent = TRUE)
incProgress(1.00, detail="Selesai")
list(idw=idw_mask, ok_pred=ok_pred_mask, ok_var=ok_var_mask, rmse_ok=rmse_ok, vgm_fit_txt=vgm_fit_txt)
})
})
output$map_idw <- renderLeaflet({
req(interp_results())
rr <- interp_results()$idw
validate(need(!is.null(rr), "IDW tidak terbentuk."))
pal <- colorNumeric(palette = pal_bupk(), domain = values(rr), na.color="#D1D5DB")
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(rr, colors = pal, opacity = 0.85) %>%
addLegend("bottomright", pal = pal, values = values(rr), title = paste0("IDW (", input$z_interp, ")"))
})
output$txt_kriging <- renderPrint({
req(interp_results())
res <- interp_results()
cat("Kriging status:\n")
if (is.null(res$ok_pred) || is.null(res$ok_var)) {
cat("- Kriging gagal/skip (variogram fit sulit).\n")
cat("- Coba naikkan resolusi grid (lebih kasar) atau ganti variabel Z.\n")
return()
}
cat("- Kriging berhasil.\n")
cat(sprintf("- RMSE CV: %s\n", ifelse(is.na(res$rmse_ok), "NA", round(res$rmse_ok, 4))))
cat("\nVariogram fit:\n")
if (!is.null(res$vgm_fit_txt)) cat(paste(res$vgm_fit_txt, collapse="\n"))
})
output$map_ok_pred <- renderLeaflet({
req(interp_results())
rr <- interp_results()$ok_pred
validate(need(!is.null(rr), "Kriging prediksi tidak terbentuk."))
pal <- colorNumeric(palette = pal_bupk(), domain = values(rr), na.color="#D1D5DB")
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(rr, colors = pal, opacity = 0.85) %>%
addLegend("bottomright", pal = pal, values = values(rr), title = paste0("OK Pred (", input$z_interp, ")"))
})
output$map_ok_var <- renderLeaflet({
req(interp_results())
rr <- interp_results()$ok_var
validate(need(!is.null(rr), "Kriging variansi tidak terbentuk."))
pal <- colorNumeric(palette = pal_div_bupk(), domain = values(rr), na.color="#D1D5DB")
leaflet() %>%
addProviderTiles("CartoDB.Positron") %>%
addRasterImage(rr, colors = pal, opacity = 0.85) %>%
addLegend("bottomright", pal = pal, values = values(rr), title = "OK Variance")
})
# ============================================================
# EXPORT
# ============================================================
output$dl_excel_raw <- downloadHandler(
filename = function() paste0("excel_raw_", Sys.Date(), ".csv"),
content = function(file){
req(excel_raw())
write.csv(excel_raw(), file, row.names=FALSE)
}
)
output$dl_clean <- downloadHandler(
filename = function() paste0("jatim_clean_", Sys.Date(), ".csv"),
content = function(file){
req(jatim_sf_full())
write.csv(jatim_sf_full() %>% st_drop_geometry(), file, row.names=FALSE)
}
)
output$dl_lisa <- downloadHandler(
filename = function() paste0("lisa_", Sys.Date(), ".csv"),
content = function(file){
req(lisa_tbl())
write.csv(lisa_tbl(), file, row.names=FALSE)
}
)
output$dl_models <- downloadHandler(
filename = function() paste0("model_compare_", Sys.Date(), ".csv"),
content = function(file){
req(model_pack())
write.csv(model_pack()$cmp, file, row.names=FALSE)
}
)
}
shinyApp(ui, server)