Analisis Spasial Kasus Tuberkulosis di Kabupaten/Kota Provinsi Jawa
Timur Tahun 2024
Analisis Spasial Kasus Tuberkulosis di Kabupaten/Kota Provinsi Jawa Timur Tahun 2024
Fazila Azra Anggina / 140610230039
Abstrak Tuberkulosis (TB) merupakan salah satu penyakit menular yang masih menjadi masalah kesehatan global, terutama di Indonesia. Jawa Timur menempati peringkat kedua dengan jumlah kasus TB tertinggi di Indonesia pada tahun 2024. Penelitian ini bertujuan untuk menganalisis faktor-faktor yang berhubungan dengan penyebaran kasus TB di Provinsi Jawa Timur menggunakan pendekatan ekonometrika spasial. Data yang digunakan merupakan data sekunder tahun 2024 yang mencakup 38 kabupaten/kota, dengan variabel dependen jumlah kasus TB (Y) dan lima variabel independen yaitu kepadatan penduduk (X1), jumlah rumah sakit (X2), konsumsi rokok kretek tanpa filter (X3), dan jumlah penduduk miskin (X4). Hasil uji Moran’s I sebesar 0,37 dan Geary’s C sebesar 0,53 menunjukkan adanya autokorelasi spasial positif antarwilayah. Berdasarkan nilai Akaike Information Criterion, model terbaik yang diperoleh adalah SDEM sebesar 564,469. Model ini menunjukkan bahwa variabel jumlah rumah sakit dan jumlah penduduk miskin berpengaruh signifikan terhadap penyebaran kasus TB, baik secara langsung maupun melalui pengaruh tetangga (spatial lag). Selain itu, koefisien error spasial (\(\lambda = 0,356\)) signifikan, mengindikasikan adanya dependensi spasial pada residual. Hasil uji Moran dan Geary terhadap residual menunjukkan tidak adanya autokorelasi spasial tersisa (p-value > 0,05), sehingga model SDEM dinilai efektif dalam menjelaskan variasi spasial penyebaran TB di Jawa Timur.
1 Pendahuluan
1.1 Latar Belakang
Tuberkulosis (TB) merupakan penyakit infeksi pada saluran pernapasan yang disebabkan oleh bakteri Mycobacterium tuberculosis dan menular dari seorang individu kepada individu lainnya melalui droplet udara yang terhirup sehingga terperangkap dalam saluran napas bagian atas (Aulia et al., 2024). Menurut laporan dari World Health Organization (WHO), TB merupakan salah satu dari 10 penyebab utama kematian akibat penyakit infeksi. Jumlah kematian yang disebabkannya pada tahun 2017 mencapai angka 1,3 juta jiwa dan bahkan tercatat muncul 6,4 juta kasus baru di seluruh dunia (WHO, 2017). Kasus dari penyakit TB pun kian meningkat setiap harinya, terbukti dengan tertularnya 1 orang dan meninggalnya 8-12 orang setiap 30 detiknya (Dwi et al., 2021).
Di Indonesia, penyakit TB menjadi penyakit yang harus ditangani. Hasil riset Kementerian Kesehatan RI 2018 menyebutkan beberapa provinsi dengan angka prevalensi TB yang berada di atas angka nasional di antaranya adalah Provinsi Aceh, DKI Jakarta, Daerah Istimewa Yogyakarta, Sumatera Barat, Kepulauan Riau, Nusa Tenggara Barat, Nusa Tenggara, Sulawesi Selatan, Sulawesi Tengah, dan Daerah Timur Indonesia (Riskesdas, 2018). Namun, Rahmawati et al. (2024) menambahkan bahwa salah satu provinsi lainnya yang juga menyumbang angka TB terbanyak di Indonesia adalah Jawa Timur, bahkan menempati peringkat kedua di Indonesia pada tahun 2024. Proporsi kasus TB pada laki-laki lebih tinggi dibandingkan perempuan, yakni sebesar 51.108 kasus (56,2%) untuk laki-laki dan 39.822 (43,8%) untuk perempuan. Tingginya mobilitas dan faktor risiko seperti kebiasaan merokok serta konsumsi alkohol turut berpengaruh terhadap persentase jumlah tersebut (Profil Kesehatan Jawa Timur, 2024).
Pertumbuhan penduduk yang semakin pesat merupakan salah satu tantangan di Indonesia. Bertambahkan jumlah penduduk tidak menjamin bahwa hasilnya akan secara otomatis meningkatkan kualitas hidup. Hal ini karena pertambahan penduduk yang cukup tinggi akan menimbulkan masalah baru apabila tidak diikuti dengan daya dukung dan tampung lingkungan yang ideal (Sari et al., 2023), contohnya kepadatan penduduk. Kepadatan penduduk sendiri adalah perbandingan antara antara jumlah penduduk dengan luas wilayah (Rao & Jhonson, 2021). Semakin banyak jumlah penduduknya, maka semakin tinggi pula kepadatan penduduk di wilayah tersebut. Salah satu dampaknya dari kondisi pemukiman yang memiliki hunian padat dapat menghalangi pencahayaan untuk masuk ke dalam rumah. Padahal, faktor pencahayaan merupakan faktor risiko yang cukup signifikan karena cahaya matahari merupakan salah satu faktor yang mampu membunuh bakteri Mycobacterium tuberculosis (Asrijun, 2018). Demikian, kepadatan penduduk yang tinggi memungkinkan terjadinya kasus penularan TB yang semakin besar.
Berdasarkan laporan Kementerian Kesehatan RI 2022 mengenai program penanggulangan TB, penemuan kasus TB dapat dilakukan secara aktif masif di masyarakat dan pasif intensif di fasilitas pelayanan kesehatan (fasyankes), dimana ini merupakan upaya menemukan terduga TB melalui skrining. Seluruh fasyankes yang memberi pelayanan kesehatan dan berpotensi menemukan terduga TB perlu terlibat dalam program penanggulangan. Upaya intensifikasi pelibatan fasyankes dalam program TB terbukti menghasilkan peningkatan penemuan kasus secara keseluruhan.
Selain itu, kebiasaan buruk seperti merokok pun turut berpotensi dalam menaikkan risiko terjangkit TB sebagaimana dijelaskan oleh dr. Agi Hidjri Tarigan, M.Ked (Paru), Sp. P. Kandungan zat berbahaya dalam rokok dapat melemahkan sistem kekebalan tubuh dan merusak silia di saluran pernapasan yang berfungsi mengeluarkan kuman, bakteri, dan virus. Dengan kata lain, rokok bisa menjadi pintu masuk bagi bakteri yang menyebabkan penyakit TB (PPTI, 2024). Tidak hanya kebiasaan hidup, keadaan sosial ekonomi yang rendah memiliki pengaruh yang positif signifikan terhadap kasus Tuberkulosis (Sihaloho et al., 2019). Hal ini dapat disebabkan oleh ketidakmampuan mengakses fasilitas kesehatan (Giyarsih, 2014). Karena keterbatasan dalam memenuhi syarat-syarat kesehatan dan lingkungan yang baik, maka keadaan sosial bisa menyebabkan terjadinya kasus TB.
Dari uraian latar belakang di atas, penelitian ini dimaksudkan untuk menganalisis faktor-faktor yang berhubungan dengan penyebaran kasus TB di Indonesia, khususnya dengan pendekatan spasial, karena fenomena ini tidak hanya dipengaruhi oleh faktor individu, melainkan juga oleh kondisi spasial seperti kepadatan penduduk, perilaku masyarakat, serta ketersediaan fasilitas kesehatan di sekitar wilayah tempat tinggal.
1.2 Identifikasi Masalah
Berdasarkan latar belakang yang telah dijelaskan sebelumnya, dapat diidentifikasi beberapa permasalahan yang menjadi dasar penelitian ini, antara lain sebagai berikut :
- Apakah terdapat autokorelasi spasial pada kasus TB di Provinsi Jawa Timur tahun 2024?
- Faktor-faktor apa saja yang secara signifikan memengaruhi tingkat penyebaran kasus TB di Provinsi Jawa Timur tahun 2024?
- Model spasial apa yang paling sesuai untuk menjelaskan hubungan antara faktor-faktor penjelas terhadap kasus TB di Provinsi Jawa Timur tahun 2024?
1.3 Tujuan Penelitian
Tujuan dari penelitian ini dapat dirumuskan ke dalam beberapa poin sebagai berikut :
- Menganalisis pola persebaran spasial kasus TB di Provinsi Jawa Timur tahun 2024 untuk mengetahui apakah terdapat dependensi spasial antar wilayah.
- Mengidentifikasi faktor-faktor yang berpengaruh signifikan terhadap tingkat penyebaran kasus TB di Provinsi Jawa Timur tahun 2024.
- Menentukan model spasial yang paling sesuai dalam menjelaskan hubungan antara faktor-faktor penjelas dengan penyebaran kasus TB di Provinsi Jawa Timur tahun 2024.
1.4 Batasan Penelitian
Penelitian ini hanya mencakup 38 kabupaten/kota di Provinsi Jawa Timur menggunakan data sekunder tahun 2024, sehingga hasil analisis yang diperoleh tidak dapat digeneralisasikan untuk provinsi lain di Indonesia. Penelitian ini juga bersifat cross-sectional dan dilakukan tanpa mempertimbangkan dinamika waktu kasus TB dari waktu ke waktu.
2 Tinjauan Pustaka
2.1 Dependensi Spasial
Dependensi spasial menyatakan bahwa nilai suatu variabel pada satu wilayah tidak independen dari nilai variabel yang sama pada wilayah tetangga. Hal ini mencerminkan bahwa “barang-apa yang berada di dekat lebih berkaitan dibanding yang jauh” sesuai dengan Tobler’s First Law of Geography. Pada penelitian ini, dependensi spasial dapat muncul karena jumlah penduduk, kondisi lingkungan, akses layanan kesehatan antar daerah, dan kebiasaan hidup.
Beberapa implikasi dependensi spasial diantara lain :
- Variabel dalam satu wilayah dapat memengaruhi variabel di wilayah tetangga (spillover).
- Observasi pada satu unit wilayah tidak independen, maka terjadi pelanggaran asumsi klasik regresi.
- Memerlukan pendekatan analisis spasial untuk menangkap efek ini.
2.2 Autokorelasi Spasial
Autokorelasi spasial mengukur sejauh mana suatu wilayah memiliki kesamaan nilai dengan wilayah di sekitarnya. Autokorelasi positif menunjukkan bahwa wilayah dengan nilai tinggi (atau rendah) cenderung berdekatan dengan wilayah lain yang juga memiliki nilai serupa (cluster), sedangkan autokorelasi negatif menunjukkan pola penyebaran yang saling berbeda (dispersion).
Ukuran utama yang digunakan untuk mengukur autokorelasi spasial adalah Moran’s I dan Local Indicators of Spatial Association (LISA) dengan formulasi sebagai berikut.
2.2.1 Rumus Moran’s I (Global Spatial Autocorrelation)
\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij}(y_i - \bar{y})(y_j - \bar{y})} {\sum_i (y_i - \bar{y})^2} \]
Keterangan :
\(n\) = jumlah observasi
\(y_i, y_j\) = nilai variabel pada lokasi \(i\) dan \(j\)
\(\bar{y}\) = rata-rata keseluruhan
\(w_{ij}\) = bobot spasial antara wilayah \(i\) dan \(j\)
Nilai I > 0 menunjukkan autokorelasi spasial positif, I < 0 menunjukkan negatif, dan I = 0 berarti tidak ada autokorelasi (Anselin, 1995).
2.2.2 Rumus Geary’s C
\[ C = \frac{(n-1)\sum_i\sum_j w_{ij}(y_i - y_j)^2} {2\,\sum_i\sum_j w_{ij}\,\sum_i(y_i - \bar y)^2} \]
Nilai \(C\) < 1 menunjukkan bahwa terdapat autokorelasi positif, \(C\) > 1 menunjukkan bahwa terdapat autokorelasi negatif.
2.2.3 Rumus LISA (Local Indicators of Spatial Association)
\[ I_i = (y_i - \bar{y}) \sum_j w_{ij}(y_j - \bar{y}) \]
LISA digunakan untuk mengidentifikasi klaster lokal, seperti daerah High-High (tinggi dikelilingi tinggi), Low-Low, High-Low, dan Low-High.
2.2.4 Scatterplot Moran & Hotspot
\[ G_i^* = \frac{\sum_j w_{ij} x_j - \bar X \sum_j w_{ij}} {S \sqrt{[ n\sum_j w_{ij}^2 - (\sum_j w_{ij})^2 ] / (n-1) }} \]
dimana \(x_{j}\) adalah nilai variabel tetangga dan \(S\) standar deviasi. Hotspot menandakan klaster tinggi yang signifikan secara spasial.
2.3 Model Ekonometrika Spasial
Model ekonometrik spasial digunakan untuk mengatasi pelanggaran asumsi independensi residual dalam model regresi konvensional. Menurut Anselin (1988), model-model ini memperhitungkan adanya efek spasial yang mungkin berasal dari pengaruh tetangga (spatial lag) atau dari error yang berkorelasi secara spasial (spatial error).
Beberapa model utama yang digunakan meliputi:
2.3.1 Spatial Autoregressive Model (SAR)
\[ y = \rho W y + X\beta + \varepsilon \]
Keterangan :
\(\rho\) = koefisien autokorelasi spasial pada variabel dependen
\(Wy\) = lag spasial dari \(y\)
\(X\beta\) = variabel independen dan koefisiennya
\(\varepsilon\) = error
Model ini menunjukkan bahwa nilai \(y\) pada suatu lokasi dipengaruhi oleh nilai \(y\) pada lokasi di sekitarnya.
2.3.2 Spatial Error Model (SEM)
\[ y = X\beta + u, \quad u = \lambda W u + \varepsilon \]
Keterangan :
\(y\) = vektor variabel dependen
\(X\) = matriks variabel independen
\(\beta\) = vektor koefisien regresi
\(u\) = komponen error yang mengandung dependensi spasial
\(\lambda\) = koefisien spatial error dependence
\(W\) = matriks bobot spasial
\(\varepsilon\) = error
Model ini digunakan ketika efek spasial muncul pada residual, bukan pada variabel dependen. Nilai \(\lambda\) menunjukkan tingkat ketergantungan error spasial.
2.3.3 Spatial Durbin Model (SDM)
\[ y = \rho W y + X\beta + W X\theta + \varepsilon \]
Keterangan :
\(\rho\) = koefisien autokorelasi spasial pada variabel dependen
\(Wy\) = lag spasial dari \(y\)
\(\theta\) = vektor koefisien untuk lag spasial dari variabel independen
\(WX\) = lag spasial dari variabel independen
\(\varepsilon\) = error
Model ini memperhitungkan efek lag pada variabel dependen maupun independen, sehingga lebih fleksibel dalam menangkap pengaruh spasial.
2.3.4 Spatial Durbin Error Model (SDEM)
\[ y = X\beta + W X\theta + u, \quad u = \lambda W u + \varepsilon \]
Keterangan :
\(X\beta\) = pengaruh langsung dari variabel independen
\(WX\theta\) = pengaruh tidak langsung dari variabel independen tetangga
\(u\) = komponen error yang juga memiliki autokorelasi spasial
\(\lambda\) = koefisien dependensi spasial pada error
SDEM menggabungkan lag spasial dari variabel independen serta error yang bersifat spasial, menjadikannya cocok untuk fenomena lingkungan dan sosial seperti penyebaran penyakit.
3 Metodologi Penelitian
3.1 Data dan Unit Spasial
Penelitian ini menggunakan data sekunder yang bersumber dari Badan Pusat Statistik dan Profil Kesehatan Jawa Timur tahun 2024 yang dapat diakses melalui https://dinkes.jatimprov.go.id/userfile/dokumen/PROFIL%20KESEHATAN%20PROVINSI%20JAWA%20TIMUR%20TAHUN%202024.pdf. Unit analisis yang dihimpun mencakup 38 kabupaten/kota di Provinsi Jawa Timur.
Dari penelitian ini pula diperoleh 1 variabel dependen dan 5 variabel independen yang dijelaskan dalam tabel keterangan variabel di bawah ini :
| Simbol | Variabel | Keterangan |
|---|---|---|
| Y | Kasus Tuberkulosis | Dependen |
| X1 | Kepadatan penduduk per km² | Independen |
| X2 | Jumlah rumah sakit | Independen |
| X3 | Konsumsi rokok kretek tanpa filter | Independen |
| X4 | Jumlah penduduk miskin | Independen |
3.2 Metode Analisis
Analisis dilakukan dengan pendekatan ekonometrika spasial yang digunakan untuk mengetahui hubungan antara variabel dependen (jumlah kasus TB) dengan berbagai faktor sosial dan ekonomi yang memengaruhinya pada tingkat kabupaten/kota di Provinsi Jawa Timur.
Tahapan analisis dimulai dengan menyajikan statistik deskriptif dan peta awal untuk memahami pola spasial dari data. Selanjutnya, akan dilakukan pengujian autokorelasi spasial global menggunakan Global Moran’s I untuk melihat apakah terdapat keterkaitan spasial secara keseluruhan. Apabila hasil uji menunjukkan adanya autokorelasi spasial yang signifikan, maka dilakukan analisis spasial lokal menggunakan Local Indicators of Spatial Association untuk mengidentifikasi klaster wilayah dengan kasus tinggi atau rendah.
Selanjutnya, dilakukan pemodelan menggunakan beberapa pendekatan, yaitu :
Model Spatial Autoregressive (SAR) untuk menguji pengaruh spasial pada variabel dependen.
Model Spatial Error (SEM) untuk menguji pengaruh spasial pada error.
Model Spatial Durbin Error (SDEM) untuk menguji pengaruh spasial, baik pada variabel independen maupun error.
Pemilihan model terbaik dilakukan berdasarkan nilai Akaike Information Criterion (AIC) dan uji signifikansi koefisien spasial.
3.3 Alur Kerja Penelitian
Alur kerja penelitian ini meliputi beberapa tahapan yang saling berurutan sebagaimana berikut :
Pengumpulan Data
Data yang digunakan meliputi data kasus TB, Jumlah Penduduk, Kepadan Penduduk per \(Km^2\) , Jumlah Rumah Sakit, Perokok Tanpa Filter, Jumlah Penduduk Miskin, dan data spasial wilayah administrasi Provinsi Jawa Timur.
Eksplorasi dan Pemetaan Data
Melakukan statistik deskriptif dan membuat peta sebaran kasus TBC untuk mengidentifikasi pola awal distribusi spasial.
Uji Autokorelasi Spasial (Global Moran’s I)
Menguji apakah terdapat hubungan spasial antar wilayah berdasarkan distribusi kasus TB.
Uji Spasial Lokal (LISA)
Mengidentifikasi wilayah-wilayah yang membentuk klaster High-High, Low-Low, atau Low-High dan High-Low.
Uji Getis–Ord Gi
Mengidentifikasi wilayah yang termasuk hotspot dan coldspot.
Pemodelan Spasial (SAR, SEM, SDEM)
Mengidentifikasi model spasial dan menentukan model terbaik berdasarkan nilai AIC.
Interpretasi
Memahami dan menjelaskan hubungan antar variabel, baik dalam wilayah itu sendiri maupun pengaruhnya terhadap wilayah tetangga.
Cross-Validation (Holdout 70/30) menggunakan RMSE per model
Mengevaluasi performa model spasial dengan menguji seberapa baik model mampu memprediksi data yang belum pernah dilihat sebelumnya.
Peta Prediksi dan Residual
Menunjukkan hasil estimasi kasus TB berdasarkan model spasial di setiap kabupaten/kota.
4 Hasil dan Pembahasan
4.1 Import Data
suppressPackageStartupMessages({
library(sf); library(sp); library(spdep); library(spatialreg)
library(dplyr); library(readxl); library(ggplot2); library(car)
library(classInt); library(scales); library(tidyr); library(gridExtra)
})## Warning: package 'sf' was built under R version 4.4.3## Warning: package 'sp' was built under R version 4.4.3## Warning: package 'spdep' was built under R version 4.4.3## Warning: package 'spData' was built under R version 4.4.3## Warning: package 'spatialreg' was built under R version 4.4.3## Warning: package 'readxl' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'car' was built under R version 4.4.3## Warning: package 'carData' was built under R version 4.4.3## Warning: package 'classInt' was built under R version 4.4.3## Warning: package 'gridExtra' was built under R version 4.4.3# ------------------------- KONFIGURASI -------------------------
path_excel <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/data tbc 2024.xlsx"
path_rds <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/gadm36_IDN_2_sp.rds"
out_dir <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL" # folder simpan output
set.seed(123) # untuk reprodusibilitas CV & impacts
# Palet biru (7 kelas) & LISA
blue_pal <- c("#E6F0FF", "#B3D1FF", "#80B3FF", "#4D94FF", "#1A75FF", "#0052CC", "#003380")
col_lisa <- c(
"High-High" = "#084594", # biru tua (indikasi hotspot kuat)
"Low-Low" = "#9ECAE1", # biru muda
"High-Low" = "#6BAED6", # biru menengah
"Low-High" = "#C6DBEF", # biru sangat muda
"Not Significant" = "#F7FBFF" # hampir putih
)
theme_set(theme_minimal(base_size = 11))
# ----------------------- FUNGSI BANTUAN ------------------------
nm_clean_base <- function(x){
x <- gsub("(?i)^(kab\\.|kabupaten)\\s+", "", x, perl = TRUE)
x <- gsub("(?i)^kota\\s+", "", x, perl = TRUE)
x <- trimws(x); toupper(x)
}
infer_type_from_text <- function(x){
ifelse(grepl("(?i)^\\s*kota\\b", x), "KOTA", "KABUPATEN")
}
make_lisa_cluster <- function(y, lw, p_cut = 0.05){
z <- scale(y)[,1]
lagz <- scale(lag.listw(lw, y))[,1]
L <- localmoran(y, lw, zero.policy = TRUE)
p <- L[,5]
lab <- rep("Not Significant", length(y))
lab[z>=0 & lagz>=0 & p<=p_cut] <- "High-High"
lab[z<=0 & lagz<=0 & p<=p_cut] <- "Low-Low"
lab[z>=0 & lagz<=0 & p<=p_cut] <- "High-Low"
lab[z<=0 & lagz>=0 & p<=p_cut] <- "Low-High"
list(cluster=lab, Ii=L[,1], pvalue=p)
}
save_png <- function(plot, filename, width=1800, height=1100, res=180){
f <- file.path(out_dir, filename)
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
png(f, width=width, height=height, res=res); print(plot); dev.off()
message("Saved: ", f)
}
rmse <- function(obs, pred){
obs <- as.numeric(obs); pred <- as.numeric(pred)
sqrt(mean((obs - pred)^2, na.rm = TRUE))
}
# ======================== IMPORT DATA ========================
# ---- 1.1 Dataset Excel ----
raw_xl <- read_excel(path_excel, sheet = 1)
# Kolom wajib sesuai file kamu
req_cols <- c("Kabupaten/Kota",
"Tuberkolosis (Y)",
"Kepadatan Penduduk per km persegi (Km2)",
"Jumlah Rumah Sakit",
"Rokok kretek tanpa filter",
"Jumlah Penduduk Miskin",
"Prevalensi")
stopifnot(all(req_cols %in% names(raw_xl)))
dat_xl <- raw_xl %>%
mutate(
P_ORIG = as.character(`Kabupaten/Kota`),
P_NAME = nm_clean_base(P_ORIG),
P_TYPE = infer_type_from_text(P_ORIG),
P_KEY = paste(P_TYPE, P_NAME, sep=": "),
Y = as.numeric(`Tuberkolosis (Y)`),
X1 = as.numeric(`Kepadatan Penduduk per km persegi (Km2)`),
X2 = as.numeric(`Jumlah Rumah Sakit`),
X3 = as.numeric(`Rokok kretek tanpa filter`),
X4 = as.numeric(`Jumlah Penduduk Miskin`),
prev = as.numeric(`Prevalensi`)
) %>%
select(P_KEY, P_TYPE, P_NAME, Y, X1:prev)
# ---- 1.2 RDS -> sf -> filter Jatim -> normalisasi & dissolve ----
indo_sp <- readRDS(path_rds)
stopifnot(inherits(indo_sp, "Spatial"))
indo_sf <- st_as_sf(indo_sp)
stopifnot(all(c("NAME_1","NAME_2") %in% names(indo_sf)))
jatim <- indo_sf %>% filter(NAME_1 %in% c("Jawa Timur","East Java"))
stopifnot(nrow(jatim) > 0)
name_clean <- nm_clean_base(jatim$NAME_2)
tipe_sp <- if ("TYPE_2" %in% names(jatim)) {
ifelse(grepl("(?i)kota", jatim$TYPE_2), "KOTA", "KABUPATEN")
} else infer_type_from_text(jatim$NAME_2)
jatim <- jatim %>%
mutate(REGION_NAME = name_clean,
REGION_KEY = paste(tipe_sp, name_clean, sep=": ")) %>%
group_by(REGION_KEY, REGION_NAME) %>%
summarise(.groups="drop") %>%
st_make_valid()
message(sprintf("Fitur Jatim setelah dissolve: %d (target ~38)", nrow(jatim)))## Fitur Jatim setelah dissolve: 38 (target ~38)# ---- 1.3 Join ----
dat_sf <- left_join(jatim, dat_xl, by = c("REGION_KEY" = "P_KEY")) %>%
filter(!is.na(Y))
if (nrow(dat_sf) == 0) stop("Join kosong. Periksa penulisan nama kab/kota di Excel.")4.2 Peta Deskriptif dan Statistik Deskriptif
# ==================== EKSPLORASI DESKRIPTIF =================
cat("\n--- SUMMARY (Y & X1..X4) ---\n")##
## --- SUMMARY (Y & X1..X4) ---print(summary(dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)))## Y X1 X2 X3
## Min. : 322 Min. : 411.0 Min. : 2.000 Min. :0.511
## 1st Qu.:1004 1st Qu.: 648.2 1st Qu.: 4.000 1st Qu.:2.041
## Median :1508 Median : 862.5 Median : 7.500 Median :2.743
## Mean :2036 Mean :1954.3 Mean : 9.474 Mean :2.883
## 3rd Qu.:2330 3rd Qu.:1214.0 3rd Qu.:10.750 3rd Qu.:3.933
## Max. :9182 Max. :8698.0 Max. :47.000 Max. :6.344
## X4
## Min. : 6.59
## 1st Qu.: 68.07
## Median :107.49
## Mean :104.81
## 3rd Qu.:146.44
## Max. :240.14# Histogram & Boxplot Y
p_hist <- ggplot(dat_sf, aes(x=Y)) +
geom_histogram(aes(y=..count..), bins=12, fill="#1A75FF", color="white") +
geom_text(
stat="bin",
bins=12,
aes(label=..count.., y=..count..),
vjust=-0.3,
color="black",
size=3
) +
labs(title="Histogram Kasus TBC (Y)", x="TBC", y="Frekuensi") +
theme_minimal()
median_y <- median(dat_sf$Y, na.rm = TRUE)
p_box <- ggplot(dat_sf, aes(y = Y)) +
geom_boxplot(fill = "#B3D1FF", color = "#003380") +
geom_text(
aes(y = median_y, label = paste("Median =", round(median_y, 1))),
x = 1.05,
hjust = 0,
size = 3.5,
color = "black"
) +
labs(title = "Boxplot Kasus TBC (Y)", y = "TBC") +
theme_minimal() +
coord_cartesian(clip = "off")
gridExtra::grid.arrange(p_hist, p_box, ncol=2)## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in geom_text(aes(y = median_y, label = paste("Median =", round(median_y, : All aesthetics have length 1, but the data has 38 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
save_png(p_hist, "01_hist_Y.png"); save_png(p_box, "02_box_Y.png")## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/01_hist_Y.png## Warning in geom_text(aes(y = median_y, label = paste("Median =", round(median_y, : All aesthetics have length 1, but the data has 38 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/02_box_Y.png# Korelasi numerik (Y & prediktor)
num_df <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
cor_mat <- round(cor(num_df, use="complete.obs"), 3)
cat("\n--- Korelasi (Y, X1..X4) ---\n"); print(cor_mat)##
## --- Korelasi (Y, X1..X4) ---## Y X1 X2 X3 X4
## Y 1.000 0.251 0.903 -0.327 0.442
## X1 0.251 1.000 0.374 -0.576 -0.528
## X2 0.903 0.374 1.000 -0.320 0.225
## X3 -0.327 -0.576 -0.320 1.000 0.175
## X4 0.442 -0.528 0.225 0.175 1.000# Peta awal Y (kuantil 7 kelas)
br_y <- classInt::classIntervals(dat_sf$Y, n=7, style="quantile")$brks
dat_sf$cutY <- cut(dat_sf$Y, breaks=br_y, include.lowest=TRUE)
p_map_y <- ggplot(dat_sf) +
geom_sf(aes(fill=cutY), color="white", size=0.25) +
scale_fill_manual(values=blue_pal, name="Kuantil Y") +
labs(title="Sebaran TBC (Y) — Kab/Kota Jawa Timur")
print(p_map_y)
library(dplyr)
top5_tb <- dat_sf %>%
st_set_geometry(NULL) %>%
select(REGION_KEY, Y) %>%
arrange(desc(Y)) %>%
head(5)
cat("\n 5 Kabupaten/Kota dengan Kasus TBC Tertinggi : \n")##
## 5 Kabupaten/Kota dengan Kasus TBC Tertinggi :print(top5_tb)## # A tibble: 5 × 2
## REGION_KEY Y
## <chr> <dbl>
## 1 KOTA: SURABAYA 9182
## 2 KABUPATEN: SIDOARJO 5409
## 3 KABUPATEN: JEMBER 4771
## 4 KABUPATEN: GRESIK 3440
## 5 KABUPATEN: PASURUAN 32964.3 Multikolinearitas
# Multikolinearitas
df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
ols <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
cat("\n--- VIF (OLS) ---\n"); print(vif(ols))##
## --- VIF (OLS) ---## X1 X2 X3 X4
## 2.774826 1.633365 1.553310 1.9783604.4 Uji Autokorelasi Spasial
# ============ BOBOT SPASIAL & AUTOKORELASI =============
nb <- poly2nb(as_Spatial(dat_sf), queen=TRUE)## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style="W", zero.policy=TRUE)
# Diagnostik konektivitas
deg <- sapply(nb, length)
cat("\n--- DIAGNOSTIK JARINGAN ---\n")##
## --- DIAGNOSTIK JARINGAN ---cat("Rata-rata tetangga:", mean(deg),
"| Min:", min(deg), "| Max:", max(deg), "\n")## Rata-rata tetangga: 2.5 | Min: 1 | Max: 5if (any(deg==0)) cat("PERINGATAN: Ada kabupaten/kota (tanpa tetangga)\n")
# Moran's I & Geary's C (global)
cat("\n--- Moran's I (Y) ---\n"); print(moran.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Moran's I (Y) ---##
## Moran I test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 3.118, p-value = 0.0009103
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.38238167 -0.03125000 0.01759825cat("\n--- Geary's C (Y) ---\n"); print(geary.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Geary's C (Y) ---##
## Geary C test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.7086, p-value = 0.003379
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.47210729 1.00000000 0.03798472# Moran scatterplot
Y_std <- scale(dat_sf$Y)[,1]
WY_std <- scale(lag.listw(lw, dat_sf$Y))[,1]
p_moran <- ggplot(data.frame(Y_std, WY_std),
aes(x=Y_std, y=WY_std)) +
geom_point(alpha=.9, size=2, color="#1A75FF") +
geom_smooth(method="lm", se=FALSE, linetype="dashed", color="#003380") +
geom_vline(xintercept=0, linetype="dotted") +
geom_hline(yintercept=0, linetype="dotted") +
labs(title="Moran Scatterplot (Y distandardisasi)",
x="Y (std)", y="W * Y (std)")
print(p_moran); save_png(p_moran, "04_moran_scatter.png")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/04_moran_scatter.png# LISA (Local Moran) + label klaster
lisa <- make_lisa_cluster(dat_sf$prev, lw, p_cut=0.05)
dat_sf$Ii <- lisa$Ii
dat_sf$Pvalue <- lisa$pvalue
dat_sf$Cluster <- factor(lisa$cluster,
levels=c("High-High","Low-Low","High-Low","Low-High","Not Significant")
)
cat("\n--- RINGKASAN LISA ---\n"); print(table(dat_sf$Cluster))##
## --- RINGKASAN LISA ---##
## High-High Low-Low High-Low Low-High Not Significant
## 0 2 0 0 36p_lisa <- ggplot(dat_sf) +
geom_sf(aes(fill=Cluster), color="white", size=0.25) +
scale_fill_manual(values=col_lisa) +
labs(title="Peta Klaster LISA — TBC Jawa Timur", fill="Klaster")
print(p_lisa); save_png(p_lisa, "05_map_LISA.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/05_map_LISA.png# ------------------ Hotspot Getis–Ord Gi* ------------------
Gi <- localG(dat_sf$prev, lw) # z-score Gi*
dat_sf$Gi_star <- as.numeric(Gi)
p_gi <- ggplot(dat_sf) +
geom_sf(aes(fill = Gi_star), color="white", size=0.25) +
scale_fill_gradient2(low="#2166AC", mid="#F7F7F7", high="#B2182B",
midpoint=0, name="Gi* z-score") +
labs(title="Getis–Ord Gi* (Hotspot TBC) — Jawa Timur")
print(p_gi); save_png(p_gi, "05b_map_Gi_star.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/05b_map_Gi_star.png4.5 Estimasi Model
# ===================== PEMODELAN SPASIAL ====================
df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
# MODEL OLS
ols <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
cat("\n--- OLS ---\n"); print(summary(ols))##
## --- OLS ---##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1950.20 -347.31 27.31 265.29 1151.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -109.27858 421.77595 -0.259 0.797174
## X1 0.06495 0.06959 0.933 0.357433
## X2 149.31095 14.28364 10.453 5.27e-12 ***
## X3 -97.25111 78.66918 -1.236 0.225110
## X4 8.43853 2.05134 4.114 0.000243 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 579.6 on 33 degrees of freedom
## Multiple R-squared: 0.8912, Adjusted R-squared: 0.878
## F-statistic: 67.59 on 4 and 33 DF, p-value: 1.991e-15# Autokorelasi residu OLS + LM-tests
cat("\n--- Moran's I pada residu OLS ---\n")##
## --- Moran's I pada residu OLS ---print(lm.morantest(ols, lw, zero.policy=TRUE))##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## weights: lw
##
## Moran I statistic standard deviate = 1.8304, p-value = 0.0336
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.20709784 -0.06953038 0.02284091cat("\n--- LM Tests (lag & error, robust) ---\n")##
## --- LM Tests (lag & error, robust) ---print(lm.LMtests(ols, lw, test=c("LMlag","LMerr","RLMlag","RLMerr")))## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSlag = 1.8334, df = 1, p-value = 0.1757
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSerr = 2.1706, df = 1, p-value = 0.1407
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSlag = 0.70664, df = 1, p-value = 0.4006
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSerr = 1.0438, df = 1, p-value = 0.3069# MODEL SAR (Spatial Lag)
sar <- lagsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SAR ---\n"); print(summary(sar))##
## --- SAR ---##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1829.71967 -307.87062 -0.26423 275.30382 1213.13159
##
## Type: lag
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -592.962223 415.925943 -1.4256 0.1540
## X1 0.115688 0.063689 1.8165 0.0693
## X2 135.774180 13.956756 9.7282 < 2.2e-16
## X3 -28.099166 76.353048 -0.3680 0.7129
## X4 8.730778 1.878674 4.6473 3.363e-06
##
## Rho: 0.13435, LR test value: 2.5464, p-value: 0.11054
## Asymptotic standard error: 0.07243
## z-value: 1.8549, p-value: 0.063609
## Wald statistic: 3.4407, p-value: 0.063609
##
## Log likelihood: -291.735 for lag model
## ML residual variance (sigma squared): 271240, (sigma: 520.81)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 597.47, (AIC for lm: 598.02)
## LM test for residual autocorrelation
## test value: 1.7797, p-value: 0.18218# MODEL SEM (Spatial Error)
sem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SEM ---\n"); print(summary(sem))##
## --- SEM ---##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1969.664 -317.666 33.163 340.973 899.496
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -381.772182 384.300309 -0.9934 0.32050
## X1 0.126928 0.062753 2.0227 0.04311
## X2 136.189753 13.829124 9.8480 < 2.2e-16
## X3 -78.651500 76.908299 -1.0227 0.30647
## X4 10.607993 1.933421 5.4866 4.096e-08
##
## Lambda: 0.22545, LR test value: 1.8464, p-value: 0.1742
## Asymptotic standard error: 0.17444
## z-value: 1.2924, p-value: 0.19622
## Wald statistic: 1.6703, p-value: 0.19622
##
## Log likelihood: -292.085 for error model
## ML residual variance (sigma squared): 273290, (sigma: 522.77)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 598.17, (AIC for lm: 598.02)# MODEL SDM (Spatial Durbin Mixed: lag Y + lag X)
sdm <- lagsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, type="mixed",
method="eigen", zero.policy=TRUE)
cat("\n--- SDM ---\n"); print(summary(sdm))##
## --- SDM ---##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## type = "mixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1619.456 -260.774 -18.994 250.872 938.312
##
## Type: mixed
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -251.694006 458.362615 -0.5491 0.58293
## X1 0.119694 0.063953 1.8716 0.06126
## X2 124.554297 14.710171 8.4672 < 2.2e-16
## X3 -96.409475 94.494533 -1.0203 0.30760
## X4 11.466432 1.961856 5.8447 5.075e-09
## lag.X1 -0.214451 0.138199 -1.5518 0.12072
## lag.X2 12.228546 31.549488 0.3876 0.69831
## lag.X3 46.012344 89.662984 0.5132 0.60783
## lag.X4 -6.269699 2.781101 -2.2544 0.02417
##
## Rho: 0.31115, LR test value: 4.2059, p-value: 0.040284
## Asymptotic standard error: 0.16
## z-value: 1.9447, p-value: 0.051816
## Wald statistic: 3.7817, p-value: 0.051816
##
## Log likelihood: -288.8023 for mixed model
## ML residual variance (sigma squared): 226290, (sigma: 475.7)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 599.6, (AIC for lm: 601.81)
## LM test for residual autocorrelation
## test value: 0.54449, p-value: 0.46058# MODEL SDEM (Spatial Durbin Error: error spatial + lag X)
sdem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, etype="emixed",
method="eigen", zero.policy=TRUE)
cat("\n--- SDEM ---\n"); print(summary(sdem))##
## --- SDEM ---##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## etype = "emixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1476.178 -260.940 -30.626 258.912 1021.531
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -524.873424 490.509301 -1.0701 0.284593
## X1 0.130491 0.066324 1.9675 0.049126
## X2 129.790950 13.716357 9.4625 < 2.2e-16
## X3 -49.676508 92.228435 -0.5386 0.590146
## X4 11.621600 1.783626 6.5157 7.234e-11
## lag.X1 -0.258342 0.147530 -1.7511 0.079925
## lag.X2 66.958801 20.522087 3.2628 0.001103
## lag.X3 -5.254513 94.046763 -0.0559 0.955444
## lag.X4 -3.327177 2.404541 -1.3837 0.166449
##
## Lambda: 0.42338, LR test value: 6.3805, p-value: 0.011538
## Asymptotic standard error: 0.15079
## z-value: 2.8077, p-value: 0.004989
## Wald statistic: 7.8834, p-value: 0.004989
##
## Log likelihood: -287.715 for error model
## ML residual variance (sigma squared): 207290, (sigma: 455.29)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 597.43, (AIC for lm: 601.81)# MODEL SAC (SARAR: lag + error)
sac <- sacsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SAC ---\n"); print(summary(sac))##
## --- SAC ---##
## Call:sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1911.334 -311.503 21.299 255.595 1004.171
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -836.425337 439.421418 -1.9035 0.056979
## X1 0.165630 0.062683 2.6424 0.008233
## X2 128.622146 13.756721 9.3498 < 2.2e-16
## X3 -10.523155 83.616686 -0.1258 0.899851
## X4 10.670614 1.897647 5.6231 1.876e-08
##
## Rho: 0.11911
## Asymptotic standard error: 0.073634
## z-value: 1.6176, p-value: 0.10575
## Lambda: 0.20023
## Asymptotic standard error: 0.19044
## z-value: 1.0514, p-value: 0.29308
##
## LR test value: 3.9511, p-value: 0.13869
##
## Log likelihood: -291.0327 for sac model
## ML residual variance (sigma squared): 258320, (sigma: 508.25)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 598.07, (AIC for lm: 598.02)4.6 Pemilihan Model Terbaik
# =================== PEMILIHAN MODEL TERBAIK =================
models <- list(OLS=ols, SAR=sar, SEM=sem, SDM=sdm, SDEM=sdem, SAC=sac)
get_wald_pvalue <- function(model) {
summ <- summary(model)
pval <- tryCatch({
if (!is.null(summ$Wald$p.value)) {
summ$Wald$p.value
} else if (!is.null(summ$Wald_p.value)) {
summ$Wald_p.value
} else {
txt <- capture.output(summary(model))
line <- grep("Wald statistic", txt, value = TRUE)
if (length(line) > 0) {
as.numeric(sub(".*p-value:\\s*([0-9.]+).*", "\\1", line))
} else {
NA
}
}
}, error = function(e) NA)
if (length(pval) == 0) pval <- NA
return(pval)
}
cmp <- do.call(rbind, lapply(names(models), function(nm) {
m <- models[[nm]]
data.frame(
Model = nm,
AIC = AIC(m),
BIC = BIC(m),
LogLik = as.numeric(logLik(m)),
Wald_p = get_wald_pvalue(m),
stringsAsFactors = FALSE
)
}))
cmp <- cmp %>%
filter(!is.na(Wald_p) & Wald_p < 0.05) %>%
arrange(AIC, BIC)
# ===================== CETAK HASIL =====================
cat("\n================ PERBANDINGAN MODEL (p < 0.05) ================\n")##
## ================ PERBANDINGAN MODEL (p < 0.05) ================print(cmp)## Model AIC BIC LogLik Wald_p
## Wald statistic3 SDEM 597.43 615.4434 -287.715 0.004989006if (nrow(cmp) > 0) {
best_name <- cmp$Model[1]
best <- models[[best_name]]
cat("\nModel TERBAIK (signifikan & AIC terendah):", best_name, "\n")
cat("AIC =", round(cmp$AIC[1], 3), "| p-value model =", round(cmp$Wald_p[1], 5), "\n")
} else {
cat("\nTidak ada model signifikan (p < 0.05)\n")
}##
## Model TERBAIK (signifikan & AIC terendah): SDEM
## AIC = 597.43 | p-value model = 0.004994.7 Uji Autokorelasi Spasial Model Terbaik
sdem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, etype="emixed",
method="eigen", zero.policy=TRUE)
res_sdem <- residuals (sdem)
# Moran's I & Geary's C (global)
cat("\n--- Moran's I (Y) ---\n"); print(moran.test(res_sdem, lw, zero.policy=TRUE))##
## --- Moran's I (Y) ---##
## Moran I test under randomisation
##
## data: res_sdem
## weights: lw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.11792, p-value = 0.5469
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04891648 -0.03125000 0.02244706cat("\n--- Geary's C (Y) ---\n"); print(geary.test(res_sdem, lw, zero.policy=TRUE))##
## --- Geary's C (Y) ---##
## Geary C test under randomisation
##
## data: res_sdem
## weights: lw
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = -2.0711, p-value = 0.9808
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 1.33773251 1.00000000 0.026592794.8 Peta Prediksi dan Residual
# =========== PETA PREDIKSI & RESIDUAL (MODEL TERBAIK) =========
if (best_name == "OLS") {
dat_sf$Y_pred <- as.numeric(predict(best, newdata=df_nogeo))
} else {
# fitted.values utk sarlm
dat_sf$Y_pred <- as.numeric(fitted.values(best))
}## This method assumes the response is known - see manual pagedat_sf$Residual <- dat_sf$Y - dat_sf$Y_pred
# Prediksi (kuantil)
br_p <- classInt::classIntervals(dat_sf$Y_pred, n=7, style="quantile")$brks
dat_sf$cutPred <- cut(dat_sf$Y_pred, breaks=br_p, include.lowest=TRUE)
p_pred <- ggplot(dat_sf) +
geom_sf(aes(fill=cutPred), color="white", size=0.25) +
scale_fill_manual(values=blue_pal, name="Kuantil Prediksi") +
labs(title=paste0("Peta Prediksi TBC — Model ", best_name))
print(p_pred); save_png(p_pred, "06_map_prediksi.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/06_map_prediksi.png# Residual
p_res <- ggplot(dat_sf) +
geom_sf(aes(fill=Residual), color="white", size=0.25) +
scale_fill_gradient2(low="#084594", mid="#F7FBFF", high="#2171B5",
midpoint=0, name="Residual") +
labs(title=paste0("Peta Residual — Model ", best_name))
print(p_res)
# Plot perbandingan
p_lin <- ggplot(dat_sf, aes(x = Y_pred, y = Y)) +
geom_point(color = "#1A75FF", size = 2, alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, color = "red", linewidth = 0.8) +
geom_abline(slope = 1, intercept = 0, color = "black", linetype = "dashed") +
labs(
title = paste0("Plot Prediksi vs Aktual ", best_name),
x = "Prediksi Kasus TBC",
y = "Aktual Kasus TBC"
) +
theme_minimal()
print(p_lin)## `geom_smooth()` using formula = 'y ~ x'
4.9 Ringkasan Output
# ======================== OUTPUT ANALISIS ======================
cat("\n================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================\n")##
## ================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================print(cmp)## Model AIC BIC LogLik Wald_p
## Wald statistic3 SDEM 597.43 615.4434 -287.715 0.004989006cat("\n================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================\n")##
## ================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================dat_no_geom <- dat_sf %>% st_set_geometry(NULL)
print(dplyr::glimpse(dat_no_geom))## Rows: 38
## Columns: 18
## $ REGION_KEY <chr> "KABUPATEN: BANGKALAN", "KABUPATEN: BANYUWANGI", "KABUPATE…
## $ REGION_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ P_TYPE <chr> "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUP…
## $ P_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ Y <dbl> 1770, 2715, 1101, 2250, 1521, 3440, 4771, 2263, 2349, 3225…
## $ X1 <dbl> 847, 488, 724, 573, 510, 1086, 786, 1228, 1109, 786, 638, …
## $ X2 <dbl> 4, 13, 8, 10, 3, 19, 14, 14, 9, 19, 8, 3, 3, 24, 11, 4, 5,…
## $ X3 <dbl> 4.057, 3.121, 2.824, 3.231, 3.560, 2.083, 3.186, 2.156, 4.…
## $ X4 <dbl> 190.94, 106.61, 95.91, 147.33, 99.62, 142.39, 224.77, 110.…
## $ prev <dbl> 160.54422, 154.75376, 87.12511, 169.77288, 191.97274, 252.…
## $ cutY <fct> "(1.45e+03,2.03e+03]", "(2.27e+03,3.21e+03]", "(715,1.11e+…
## $ Ii <dbl> 0.210313945, 0.049590395, 0.776928630, 0.088984879, 0.0214…
## $ Pvalue <dbl> 0.59429929, 0.83050575, 0.19181716, 0.48737496, 0.68333028…
## $ Cluster <fct> Not Significant, Not Significant, Not Significant, Not Sig…
## $ Gi_star <dbl> -0.5326162, -0.2140530, -1.3052223, -0.6944900, -0.4079229…
## $ Y_pred <dbl> 1221.2241, 2449.5151, 1495.4414, 2295.8641, 1271.5199, 429…
## $ Residual <dbl> 548.77591, 265.48488, -394.44142, -45.86406, 249.48005, -8…
## $ cutPred <fct> "(607,1.22e+03]", "(2.33e+03,3.26e+03]", "(1.36e+03,1.91e+…
## # A tibble: 38 × 18
## REGION_KEY REGION_NAME P_TYPE P_NAME Y X1 X2 X3 X4 prev
## * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 KABUPATEN: BAN… BANGKALAN KABUP… BANGK… 1770 847 4 4.06 191. 161.
## 2 KABUPATEN: BAN… BANYUWANGI KABUP… BANYU… 2715 488 13 3.12 107. 155.
## 3 KABUPATEN: BLI… BLITAR KABUP… BLITAR 1101 724 8 2.82 95.9 87.1
## 4 KABUPATEN: BOJ… BOJONEGORO KABUP… BOJON… 2250 573 10 3.23 147. 170.
## 5 KABUPATEN: BON… BONDOWOSO KABUP… BONDO… 1521 510 3 3.56 99.6 192.
## 6 KABUPATEN: GRE… GRESIK KABUP… GRESIK 3440 1086 19 2.08 142. 252.
## 7 KABUPATEN: JEM… JEMBER KABUP… JEMBER 4771 786 14 3.19 225. 183.
## 8 KABUPATEN: JOM… JOMBANG KABUP… JOMBA… 2263 1228 14 2.16 111. 166.
## 9 KABUPATEN: KED… KEDIRI KABUP… KEDIRI 2349 1109 9 4.24 159. 139.
## 10 KABUPATEN: LAM… LAMONGAN KABUP… LAMON… 3225 786 19 3.01 147. 234.
## # ℹ 28 more rows
## # ℹ 8 more variables: cutY <fct>, Ii <dbl>, Pvalue <dbl>, Cluster <fct>,
## # Gi_star <dbl>, Y_pred <dbl>, Residual <dbl>, cutPred <fct>cat("\n--- 15 baris pertama (kolom kunci) ---\n")##
## --- 15 baris pertama (kolom kunci) ---print(
dat_no_geom %>%
dplyr::select(REGION_NAME, Y, dplyr::starts_with("X"), Y_pred, Residual) %>%
head(15)
)## # A tibble: 15 × 8
## REGION_NAME Y X1 X2 X3 X4 Y_pred Residual
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BANGKALAN 1770 847 4 4.06 191. 1221. 549.
## 2 BANYUWANGI 2715 488 13 3.12 107. 2450. 265.
## 3 BLITAR 1101 724 8 2.82 95.9 1495. -394.
## 4 BOJONEGORO 2250 573 10 3.23 147. 2296. -45.9
## 5 BONDOWOSO 1521 510 3 3.56 99.6 1272. 249.
## 6 GRESIK 3440 1086 19 2.08 142. 4298. -858.
## 7 JEMBER 4771 786 14 3.19 225. 4122. 649.
## 8 JOMBANG 2263 1228 14 2.16 111. 2543. -280.
## 9 KEDIRI 2349 1109 9 4.24 159. 2412. -62.9
## 10 LAMONGAN 3225 786 19 3.01 147. 3543. -318.
## 11 LUMAJANG 2154 638 8 1.78 91.0 1906. 248.
## 12 MADIUN 972 681 3 4.72 73.2 569. 403.
## 13 MAGETAN 935 970 3 2.41 59.5 209. 726.
## 14 MALANG 3177 788 24 4.37 240. 4653. -1476.
## 15 MOJOKERTO 2012 1172 11 2.02 109. 2296. -284.cat("\n====================== PROPORSI KLASTER LISA ===========================\n")##
## ====================== PROPORSI KLASTER LISA ===========================prop_lisa <- round(prop.table(table(dat_sf$Cluster)), 3)
print(prop_lisa)##
## High-High Low-Low High-Low Low-High Not Significant
## 0.000 0.053 0.000 0.000 0.9475 Kesimpulan dan Saran
Diketahui persebaran kasus TB di Jawa Timur, dimana kasus prevalensi TB tertinggi di Provinsi Jawa Timur tahun 2024 terjadi di Kabupaten Sidoarjo, Kabupaten Jember, Kabupaten Gresik, Kabupaten Pasuruan, dan Kabupaten Lamongan.
Perhitungan Indeks Moran dan Geary dilakukan pada banyaknya kasus TB yang terjadi sebagai variabel dependen untuk mengetahui apakah terdapat indikasi dependensi spasial atau tidak. Nilai yang diperoleh untuk keduanya berada pada rentang 0 < I < 1 (0,37 untuk Moran dan 0,53 untuk Geary), yang menunjukkan penyebaran kasus TB di Jawa Timur memiliki autokorelasi spasial yang positif. Dengan kata lain, antar kabupaten/kota yang jaraknya berdekatan memiliki nilai persebaran yang hampir sama dan cenderung berkelompok secara spasial. Hal ini turut diperkuat oleh perolehan p-value yang lebih kecil daripada \(\alpha\) (0,05) untuk kedua metode, menandakan bahwa dengan menggunakan tingkat kepercayaan 95%, terdapat dependensi spasial pada kasus TB di Provinsi Jawa Timur tahun 2024.
Berdasarkan uji LM, diperoleh bahwa tidak terdapat dependensi spasial error dan spasial lag. Akan tetapi, karena terindikasi pola autokorelasi pada residual model OLS, maka pendekatan secara spasial akan tetap dilakukan. Setelah mencobakan beberapa model, diperoleh bahwa model terbaik berdasarkan AIC terkecil adalah model SDEM yaitu sebesar 564,469.
Setelah melihat uji simultan model, diperoleh nilai p-value sebesar 0,0309. Jika dibandingkan dengan nilai \(\alpha\) (0,05), maka diperoleh kesimpulan bahwa data yang ada mendukung untuk menolak \(H_0\) yang berarti variabel independen berpengaruh secara simultan terhadap variabel dependen. Sementara untuk uji parsial, diperoleh bahwa koefisien \(\beta_2\), \(\beta_4\), dan \(\theta_2\) secara signifikan berpengaruh terhadap variabel dependen.
Sementara itu, uji asumsi residual untuk model SDEM sendiri mengindikasikan bahwa sudah tidak ada lagi dependensi spasial yang terdeteksi. Dapat terlihat setelah dilakukan uji Moran dan Geary, diperoleh p-value > \(\alpha\) (0,05) atau dengan kata lain, model SDEM sudah berhasil menangkap seluruh pola spasial yang relevan di dalam data.
Maka, persamaan untuk model spasial ekonometrika dalam penelitian ini adalah sebagai berikut : \[ \begin{equation} Y_i = 54.773 + 111.28X_{2i} + 9.893X_{4i} + 53.079(WX_{2i}) + u_i, \end{equation} \] \[ \begin{equation} u_i = 0.356 \sum_{j} w_{ij}u_j + \varepsilon_i \end{equation} \]
6 Daftar Pustaka
Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Springer.
Asrijun. (2018). Pengaruh Kondisi Lingkungan Terhadap Penularan Tuberkulosis Paru. Jurnal Kesehatan Lingkungan.
Aulia, K. . T., Yuniati, S. K., & Seto, A. A. B. (2024). The Effect of Population Density on The Risk of Tuberculosis in Densely Populated Environments.Journal of Diverse Medical Research: Medicosphere,1(2), 7–12.
Dinas Kesehatan Provinsi Jawa Timur. (2025). Profil Kesehatan Provinsi Jawa Timur 2024. Dinas Kesehatan Provinsi Jawa Timur.
Dwi, I. K., Tintin, S., & Makhfudli. (2021). Gambaran Perilaku Pengawas Minum Obat (PMO) Terhadap Sikap, Kepatuhan Minum Obat Dan Kualitas Hidup Pasien TB Paru. Jurnal Penelitian Kesehatan Suara Forikes, 12(1), 39–42.
Giyarsih, Sri Rum. 2014. “Pengentasan Kemiskinan Yang Komprehensif Di Bagian Wilayah Terluar Di Indonesia - Kasus Kabupaten Nunukan Provinsi Kalimantan Utara.” Jurnal Manusia Dan Lingkungan 21 (2): 239–46.
Kementerian Kesehatan Republik Indonesia. (2022). Program penanggulangan tuberkulosis. Kementerian Kesehatan RI.
PPTI. (2024). Rokok dan TBC: Dua ancaman serius bagi Gen-Z. Diakses dari https://ppti.id/rokok-dan-tbc-dua-ancaman-serius-bagi-gen-z/
Rahmawati, N., Karno, F., & Hermanto, E. M. P. (2024). Analisis Penyakit Tuberkulosis (TBC) pada Provinsi Jawa Timur Tahun 2021 Menggunakan Geographically Weighted Regression (GWR). Indonesian Journal of Applied Statistics, 6(2), 116. https://doi.org/10.13057/ijas.v6i2.78593
Rao, M., & Johnson, A. (2021). Impact of Population Density and Elevation on Tuberculosis Spread and Transmission in Maharashtra, India. Journal of Emerging Inversigators. https://doi.org/10.59720/21056
Riskesdas. (2018). Hasil Utama Riset Kesehatan Dasar. Kementerian Kesehatan Republik Indonesia. 1-100.
Sari, A. P., Rahmadini, G., Carlina, H., Ramadan, M. I., & Pradani, Z. E. (2023). Analisis masalah kependudukan di Indonesia. Journal of Economic Education (JEecE), 2(1), 29–37. e-ISSN: 2964-559X.
Sihaloho, E.D., Alfarizy, I.L., & Sagala, E.B. (2019). “Indikator Ekonomi Dan Angka Tuberkulosis Di Kabupaten Kota Jawa Barat.” Jurnal Ilmu Ekonomi Dan Pembangunan 19 (2): 136-46. https://doi.org/10.20961/jiep.v19i2.33698.
WHO. (2017). Global Tuberculosis Report. Geneva : World Health Organization
7 Lampiran
Laporan dalam Rpubs ini dapat diakses melalui link berikut ini :
https://rpubs.com/fazilazranggina/UTSProjekStatistikaSpasial
Laporan dashboard dapat diakses melalui link berikut ini :
https://fazilaazraanggina.shinyapps.io/DASHBOARD/
Lampiran syntax dapat dilihat melalui uraian di bawah ini :
suppressPackageStartupMessages({
library(sf); library(sp); library(spdep); library(spatialreg)
library(dplyr); library(readxl); library(ggplot2); library(car)
library(classInt); library(scales); library(tidyr); library(gridExtra)
})
# ------------------------- KONFIGURASI -------------------------
path_excel <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/data tbc 2024.xlsx"
path_rds <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/gadm36_IDN_2_sp.rds"
out_dir <- "C:/Users/hp/Documents/New folder/SMT 5/SPASIAL" # folder simpan output
set.seed(123) # untuk reprodusibilitas CV & impacts
# Palet biru (7 kelas) & LISA
blue_pal <- c("#E6F0FF", "#B3D1FF", "#80B3FF", "#4D94FF", "#1A75FF", "#0052CC", "#003380")
col_lisa <- c(
"High-High" = "#084594", # biru tua (indikasi hotspot kuat)
"Low-Low" = "#9ECAE1", # biru muda
"High-Low" = "#6BAED6", # biru menengah
"Low-High" = "#C6DBEF", # biru sangat muda
"Not Significant" = "#F7FBFF" # hampir putih
)
theme_set(theme_minimal(base_size = 11))
# ----------------------- FUNGSI BANTUAN ------------------------
nm_clean_base <- function(x){
x <- gsub("(?i)^(kab\\.|kabupaten)\\s+", "", x, perl = TRUE)
x <- gsub("(?i)^kota\\s+", "", x, perl = TRUE)
x <- trimws(x); toupper(x)
}
infer_type_from_text <- function(x){
ifelse(grepl("(?i)^\\s*kota\\b", x), "KOTA", "KABUPATEN")
}
make_lisa_cluster <- function(y, lw, p_cut = 0.05){
z <- scale(y)[,1]
lagz <- scale(lag.listw(lw, y))[,1]
L <- localmoran(y, lw, zero.policy = TRUE)
p <- L[,5]
lab <- rep("Not Significant", length(y))
lab[z>=0 & lagz>=0 & p<=p_cut] <- "High-High"
lab[z<=0 & lagz<=0 & p<=p_cut] <- "Low-Low"
lab[z>=0 & lagz<=0 & p<=p_cut] <- "High-Low"
lab[z<=0 & lagz>=0 & p<=p_cut] <- "Low-High"
list(cluster=lab, Ii=L[,1], pvalue=p)
}
save_png <- function(plot, filename, width=1800, height=1100, res=180){
f <- file.path(out_dir, filename)
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
png(f, width=width, height=height, res=res); print(plot); dev.off()
message("Saved: ", f)
}
rmse <- function(obs, pred){
obs <- as.numeric(obs); pred <- as.numeric(pred)
sqrt(mean((obs - pred)^2, na.rm = TRUE))
}
# ======================== IMPORT DATA ========================
# ---- 1.1 Dataset Excel ----
raw_xl <- read_excel(path_excel, sheet = 1)
# Kolom wajib sesuai file kamu
req_cols <- c("Kabupaten/Kota",
"Tuberkolosis (Y)",
"Kepadatan Penduduk per km persegi (Km2)",
"Jumlah Rumah Sakit",
"Rokok kretek tanpa filter",
"Jumlah Penduduk Miskin",
"Prevalensi")
stopifnot(all(req_cols %in% names(raw_xl)))
dat_xl <- raw_xl %>%
mutate(
P_ORIG = as.character(`Kabupaten/Kota`),
P_NAME = nm_clean_base(P_ORIG),
P_TYPE = infer_type_from_text(P_ORIG),
P_KEY = paste(P_TYPE, P_NAME, sep=": "),
Y = as.numeric(`Tuberkolosis (Y)`),
X1 = as.numeric(`Kepadatan Penduduk per km persegi (Km2)`),
X2 = as.numeric(`Jumlah Rumah Sakit`),
X3 = as.numeric(`Rokok kretek tanpa filter`),
X4 = as.numeric(`Jumlah Penduduk Miskin`),
prev = as.numeric(`Prevalensi`)
) %>%
select(P_KEY, P_TYPE, P_NAME, Y, X1:prev)
# ---- 1.2 RDS -> sf -> filter Jatim -> normalisasi & dissolve ----
indo_sp <- readRDS(path_rds)
stopifnot(inherits(indo_sp, "Spatial"))
indo_sf <- st_as_sf(indo_sp)
stopifnot(all(c("NAME_1","NAME_2") %in% names(indo_sf)))
jatim <- indo_sf %>% filter(NAME_1 %in% c("Jawa Timur","East Java"))
stopifnot(nrow(jatim) > 0)
name_clean <- nm_clean_base(jatim$NAME_2)
tipe_sp <- if ("TYPE_2" %in% names(jatim)) {
ifelse(grepl("(?i)kota", jatim$TYPE_2), "KOTA", "KABUPATEN")
} else infer_type_from_text(jatim$NAME_2)
jatim <- jatim %>%
mutate(REGION_NAME = name_clean,
REGION_KEY = paste(tipe_sp, name_clean, sep=": ")) %>%
group_by(REGION_KEY, REGION_NAME) %>%
summarise(.groups="drop") %>%
st_make_valid()
message(sprintf("Fitur Jatim setelah dissolve: %d (target ~38)", nrow(jatim)))## Fitur Jatim setelah dissolve: 38 (target ~38)# ---- 1.3 Join ----
dat_sf <- left_join(jatim, dat_xl, by = c("REGION_KEY" = "P_KEY")) %>%
filter(!is.na(Y))
if (nrow(dat_sf) == 0) stop("Join kosong. Periksa penulisan nama kab/kota di Excel.")
# ==================== EKSPLORASI DESKRIPTIF =================
cat("\n--- SUMMARY (Y & X1..X4) ---\n")##
## --- SUMMARY (Y & X1..X4) ---print(summary(dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)))## Y X1 X2 X3
## Min. : 322 Min. : 411.0 Min. : 2.000 Min. :0.511
## 1st Qu.:1004 1st Qu.: 648.2 1st Qu.: 4.000 1st Qu.:2.041
## Median :1508 Median : 862.5 Median : 7.500 Median :2.743
## Mean :2036 Mean :1954.3 Mean : 9.474 Mean :2.883
## 3rd Qu.:2330 3rd Qu.:1214.0 3rd Qu.:10.750 3rd Qu.:3.933
## Max. :9182 Max. :8698.0 Max. :47.000 Max. :6.344
## X4
## Min. : 6.59
## 1st Qu.: 68.07
## Median :107.49
## Mean :104.81
## 3rd Qu.:146.44
## Max. :240.14# Histogram & Boxplot Y
p_hist <- ggplot(dat_sf, aes(x=Y)) +
geom_histogram(aes(y=..count..), bins=12, fill="#1A75FF", color="white") +
geom_text(
stat="bin",
bins=12,
aes(label=..count.., y=..count..),
vjust=-0.3,
color="black",
size=3
) +
labs(title="Histogram Kasus TBC (Y)", x="TBC", y="Frekuensi") +
theme_minimal()
median_y <- median(dat_sf$Y, na.rm = TRUE)
p_box <- ggplot(dat_sf, aes(y = Y)) +
geom_boxplot(fill = "#B3D1FF", color = "#003380") +
geom_text(
aes(y = median_y, label = paste("Median =", round(median_y, 1))),
x = 1.05,
hjust = 0,
size = 3.5,
color = "black"
) +
labs(title = "Boxplot Kasus TBC (Y)", y = "TBC") +
theme_minimal() +
coord_cartesian(clip = "off")
gridExtra::grid.arrange(p_hist, p_box, ncol=2)## Warning in geom_text(aes(y = median_y, label = paste("Median =", round(median_y, : All aesthetics have length 1, but the data has 38 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
save_png(p_hist, "01_hist_Y.png"); save_png(p_box, "02_box_Y.png")## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/01_hist_Y.png## Warning in geom_text(aes(y = median_y, label = paste("Median =", round(median_y, : All aesthetics have length 1, but the data has 38 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/02_box_Y.png# Korelasi numerik (Y & prediktor)
num_df <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
cor_mat <- round(cor(num_df, use="complete.obs"), 3)
cat("\n--- Korelasi (Y, X1..X4) ---\n"); print(cor_mat)##
## --- Korelasi (Y, X1..X4) ---## Y X1 X2 X3 X4
## Y 1.000 0.251 0.903 -0.327 0.442
## X1 0.251 1.000 0.374 -0.576 -0.528
## X2 0.903 0.374 1.000 -0.320 0.225
## X3 -0.327 -0.576 -0.320 1.000 0.175
## X4 0.442 -0.528 0.225 0.175 1.000# Peta awal Y (kuantil 7 kelas)
br_y <- classInt::classIntervals(dat_sf$Y, n=7, style="quantile")$brks
dat_sf$cutY <- cut(dat_sf$Y, breaks=br_y, include.lowest=TRUE)
p_map_y <- ggplot(dat_sf) +
geom_sf(aes(fill=cutY), color="white", size=0.25) +
scale_fill_manual(values=blue_pal, name="Kuantil Y") +
labs(title="Sebaran TBC (Y) — Kab/Kota Jawa Timur")
print(p_map_y)
library(dplyr)
top5_tb <- dat_sf %>%
st_set_geometry(NULL) %>%
select(REGION_KEY, Y) %>%
arrange(desc(Y)) %>%
head(5)
cat("\n 5 Kabupaten/Kota dengan Kasus TBC Tertinggi : \n")##
## 5 Kabupaten/Kota dengan Kasus TBC Tertinggi :print(top5_tb)## # A tibble: 5 × 2
## REGION_KEY Y
## <chr> <dbl>
## 1 KOTA: SURABAYA 9182
## 2 KABUPATEN: SIDOARJO 5409
## 3 KABUPATEN: JEMBER 4771
## 4 KABUPATEN: GRESIK 3440
## 5 KABUPATEN: PASURUAN 3296# Multikolinearitas
df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
ols <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
cat("\n--- VIF (OLS) ---\n"); print(vif(ols))##
## --- VIF (OLS) ---## X1 X2 X3 X4
## 2.774826 1.633365 1.553310 1.978360# ============ BOBOT SPASIAL & AUTOKORELASI =============
nb <- poly2nb(as_Spatial(dat_sf), queen=TRUE)## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(as_Spatial(dat_sf), queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style="W", zero.policy=TRUE)
# Diagnostik konektivitas
deg <- sapply(nb, length)
cat("\n--- DIAGNOSTIK JARINGAN ---\n")##
## --- DIAGNOSTIK JARINGAN ---cat("Rata-rata tetangga:", mean(deg),
"| Min:", min(deg), "| Max:", max(deg), "\n")## Rata-rata tetangga: 2.5 | Min: 1 | Max: 5if (any(deg==0)) cat("PERINGATAN: Ada kabupaten/kota (tanpa tetangga)\n")
# Moran's I & Geary's C (global)
cat("\n--- Moran's I (Y) ---\n"); print(moran.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Moran's I (Y) ---##
## Moran I test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 3.118, p-value = 0.0009103
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.38238167 -0.03125000 0.01759825cat("\n--- Geary's C (Y) ---\n"); print(geary.test(dat_sf$Y, lw, zero.policy=TRUE))##
## --- Geary's C (Y) ---##
## Geary C test under randomisation
##
## data: dat_sf$Y
## weights: lw
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = 2.7086, p-value = 0.003379
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.47210729 1.00000000 0.03798472# Moran scatterplot
Y_std <- scale(dat_sf$Y)[,1]
WY_std <- scale(lag.listw(lw, dat_sf$Y))[,1]
p_moran <- ggplot(data.frame(Y_std, WY_std),
aes(x=Y_std, y=WY_std)) +
geom_point(alpha=.9, size=2, color="#1A75FF") +
geom_smooth(method="lm", se=FALSE, linetype="dashed", color="#003380") +
geom_vline(xintercept=0, linetype="dotted") +
geom_hline(yintercept=0, linetype="dotted") +
labs(title="Moran Scatterplot (Y distandardisasi)",
x="Y (std)", y="W * Y (std)")
print(p_moran); save_png(p_moran, "04_moran_scatter.png")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/04_moran_scatter.png# LISA (Local Moran) + label klaster
lisa <- make_lisa_cluster(dat_sf$prev, lw, p_cut=0.05)
dat_sf$Ii <- lisa$Ii
dat_sf$Pvalue <- lisa$pvalue
dat_sf$Cluster <- factor(lisa$cluster,
levels=c("High-High","Low-Low","High-Low","Low-High","Not Significant")
)
cat("\n--- RINGKASAN LISA ---\n"); print(table(dat_sf$Cluster))##
## --- RINGKASAN LISA ---##
## High-High Low-Low High-Low Low-High Not Significant
## 0 2 0 0 36p_lisa <- ggplot(dat_sf) +
geom_sf(aes(fill=Cluster), color="white", size=0.25) +
scale_fill_manual(values=col_lisa) +
labs(title="Peta Klaster LISA — TBC Jawa Timur", fill="Klaster")
print(p_lisa); save_png(p_lisa, "05_map_LISA.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/05_map_LISA.png# ------------------ Hotspot Getis–Ord Gi* ------------------
Gi <- localG(dat_sf$prev, lw) # z-score Gi*
dat_sf$Gi_star <- as.numeric(Gi)
p_gi <- ggplot(dat_sf) +
geom_sf(aes(fill = Gi_star), color="white", size=0.25) +
scale_fill_gradient2(low="#2166AC", mid="#F7F7F7", high="#B2182B",
midpoint=0, name="Gi* z-score") +
labs(title="Getis–Ord Gi* (Hotspot TBC) — Jawa Timur")
print(p_gi); save_png(p_gi, "05b_map_Gi_star.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/05b_map_Gi_star.png# ===================== PEMODELAN SPASIAL ====================
df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% select(Y, X1:X4)
# MODEL OLS
ols <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
cat("\n--- OLS ---\n"); print(summary(ols))##
## --- OLS ---##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1950.20 -347.31 27.31 265.29 1151.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -109.27858 421.77595 -0.259 0.797174
## X1 0.06495 0.06959 0.933 0.357433
## X2 149.31095 14.28364 10.453 5.27e-12 ***
## X3 -97.25111 78.66918 -1.236 0.225110
## X4 8.43853 2.05134 4.114 0.000243 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 579.6 on 33 degrees of freedom
## Multiple R-squared: 0.8912, Adjusted R-squared: 0.878
## F-statistic: 67.59 on 4 and 33 DF, p-value: 1.991e-15# Autokorelasi residu OLS + LM-tests
cat("\n--- Moran's I pada residu OLS ---\n")##
## --- Moran's I pada residu OLS ---print(lm.morantest(ols, lw, zero.policy=TRUE))##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## weights: lw
##
## Moran I statistic standard deviate = 1.8304, p-value = 0.0336
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.20709784 -0.06953038 0.02284091cat("\n--- LM Tests (lag & error, robust) ---\n")##
## --- LM Tests (lag & error, robust) ---print(lm.LMtests(ols, lw, test=c("LMlag","LMerr","RLMlag","RLMerr")))## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSlag = 1.8334, df = 1, p-value = 0.1757
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## RSerr = 2.1706, df = 1, p-value = 0.1407
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSlag = 0.70664, df = 1, p-value = 0.4006
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
##
## adjRSerr = 1.0438, df = 1, p-value = 0.3069# MODEL SAR (Spatial Lag)
sar <- lagsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SAR ---\n"); print(summary(sar))##
## --- SAR ---##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1829.71967 -307.87062 -0.26423 275.30382 1213.13159
##
## Type: lag
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -592.962223 415.925943 -1.4256 0.1540
## X1 0.115688 0.063689 1.8165 0.0693
## X2 135.774180 13.956756 9.7282 < 2.2e-16
## X3 -28.099166 76.353048 -0.3680 0.7129
## X4 8.730778 1.878674 4.6473 3.363e-06
##
## Rho: 0.13435, LR test value: 2.5464, p-value: 0.11054
## Asymptotic standard error: 0.07243
## z-value: 1.8549, p-value: 0.063609
## Wald statistic: 3.4407, p-value: 0.063609
##
## Log likelihood: -291.735 for lag model
## ML residual variance (sigma squared): 271240, (sigma: 520.81)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 597.47, (AIC for lm: 598.02)
## LM test for residual autocorrelation
## test value: 1.7797, p-value: 0.18218# MODEL SEM (Spatial Error)
sem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SEM ---\n"); print(summary(sem))##
## --- SEM ---##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1969.664 -317.666 33.163 340.973 899.496
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -381.772182 384.300309 -0.9934 0.32050
## X1 0.126928 0.062753 2.0227 0.04311
## X2 136.189753 13.829124 9.8480 < 2.2e-16
## X3 -78.651500 76.908299 -1.0227 0.30647
## X4 10.607993 1.933421 5.4866 4.096e-08
##
## Lambda: 0.22545, LR test value: 1.8464, p-value: 0.1742
## Asymptotic standard error: 0.17444
## z-value: 1.2924, p-value: 0.19622
## Wald statistic: 1.6703, p-value: 0.19622
##
## Log likelihood: -292.085 for error model
## ML residual variance (sigma squared): 273290, (sigma: 522.77)
## Number of observations: 38
## Number of parameters estimated: 7
## AIC: 598.17, (AIC for lm: 598.02)# MODEL SDM (Spatial Durbin Mixed: lag Y + lag X)
sdm <- lagsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, type="mixed",
method="eigen", zero.policy=TRUE)
cat("\n--- SDM ---\n"); print(summary(sdm))##
## --- SDM ---##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## type = "mixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1619.456 -260.774 -18.994 250.872 938.312
##
## Type: mixed
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -251.694006 458.362615 -0.5491 0.58293
## X1 0.119694 0.063953 1.8716 0.06126
## X2 124.554297 14.710171 8.4672 < 2.2e-16
## X3 -96.409475 94.494533 -1.0203 0.30760
## X4 11.466432 1.961856 5.8447 5.075e-09
## lag.X1 -0.214451 0.138199 -1.5518 0.12072
## lag.X2 12.228546 31.549488 0.3876 0.69831
## lag.X3 46.012344 89.662984 0.5132 0.60783
## lag.X4 -6.269699 2.781101 -2.2544 0.02417
##
## Rho: 0.31115, LR test value: 4.2059, p-value: 0.040284
## Asymptotic standard error: 0.16
## z-value: 1.9447, p-value: 0.051816
## Wald statistic: 3.7817, p-value: 0.051816
##
## Log likelihood: -288.8023 for mixed model
## ML residual variance (sigma squared): 226290, (sigma: 475.7)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 599.6, (AIC for lm: 601.81)
## LM test for residual autocorrelation
## test value: 0.54449, p-value: 0.46058# MODEL SDEM (Spatial Durbin Error: error spatial + lag X)
sdem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, etype="emixed",
method="eigen", zero.policy=TRUE)
cat("\n--- SDEM ---\n"); print(summary(sdem))##
## --- SDEM ---##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## etype = "emixed", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1476.178 -260.940 -30.626 258.912 1021.531
##
## Type: error
## Regions with no neighbours included:
## 12 31 32 35 37
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -524.873424 490.509301 -1.0701 0.284593
## X1 0.130491 0.066324 1.9675 0.049126
## X2 129.790950 13.716357 9.4625 < 2.2e-16
## X3 -49.676508 92.228435 -0.5386 0.590146
## X4 11.621600 1.783626 6.5157 7.234e-11
## lag.X1 -0.258342 0.147530 -1.7511 0.079925
## lag.X2 66.958801 20.522087 3.2628 0.001103
## lag.X3 -5.254513 94.046763 -0.0559 0.955444
## lag.X4 -3.327177 2.404541 -1.3837 0.166449
##
## Lambda: 0.42338, LR test value: 6.3805, p-value: 0.011538
## Asymptotic standard error: 0.15079
## z-value: 2.8077, p-value: 0.004989
## Wald statistic: 7.8834, p-value: 0.004989
##
## Log likelihood: -287.715 for error model
## ML residual variance (sigma squared): 207290, (sigma: 455.29)
## Number of observations: 38
## Number of parameters estimated: 11
## AIC: 597.43, (AIC for lm: 601.81)# MODEL SAC (SARAR: lag + error)
sac <- sacsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
cat("\n--- SAC ---\n"); print(summary(sac))##
## --- SAC ---##
## Call:sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, listw = lw,
## method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1911.334 -311.503 21.299 255.595 1004.171
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -836.425337 439.421418 -1.9035 0.056979
## X1 0.165630 0.062683 2.6424 0.008233
## X2 128.622146 13.756721 9.3498 < 2.2e-16
## X3 -10.523155 83.616686 -0.1258 0.899851
## X4 10.670614 1.897647 5.6231 1.876e-08
##
## Rho: 0.11911
## Asymptotic standard error: 0.073634
## z-value: 1.6176, p-value: 0.10575
## Lambda: 0.20023
## Asymptotic standard error: 0.19044
## z-value: 1.0514, p-value: 0.29308
##
## LR test value: 3.9511, p-value: 0.13869
##
## Log likelihood: -291.0327 for sac model
## ML residual variance (sigma squared): 258320, (sigma: 508.25)
## Number of observations: 38
## Number of parameters estimated: 8
## AIC: 598.07, (AIC for lm: 598.02)## Pemilihan Model Terbaik
# =================== PEMILIHAN MODEL TERBAIK =================
models <- list(OLS=ols, SAR=sar, SEM=sem, SDM=sdm, SDEM=sdem, SAC=sac)
get_wald_pvalue <- function(model) {
summ <- summary(model)
pval <- tryCatch({
if (!is.null(summ$Wald$p.value)) {
summ$Wald$p.value
} else if (!is.null(summ$Wald_p.value)) {
summ$Wald_p.value
} else {
txt <- capture.output(summary(model))
line <- grep("Wald statistic", txt, value = TRUE)
if (length(line) > 0) {
as.numeric(sub(".*p-value:\\s*([0-9.]+).*", "\\1", line))
} else {
NA
}
}
}, error = function(e) NA)
if (length(pval) == 0) pval <- NA
return(pval)
}
cmp <- do.call(rbind, lapply(names(models), function(nm) {
m <- models[[nm]]
data.frame(
Model = nm,
AIC = AIC(m),
BIC = BIC(m),
LogLik = as.numeric(logLik(m)),
Wald_p = get_wald_pvalue(m),
stringsAsFactors = FALSE
)
}))
cmp <- cmp %>%
filter(!is.na(Wald_p) & Wald_p < 0.05) %>%
arrange(AIC, BIC)
# ===================== CETAK HASIL =====================
cat("\n================ PERBANDINGAN MODEL (p < 0.05) ================\n")##
## ================ PERBANDINGAN MODEL (p < 0.05) ================print(cmp)## Model AIC BIC LogLik Wald_p
## Wald statistic3 SDEM 597.43 615.4434 -287.715 0.004989006if (nrow(cmp) > 0) {
best_name <- cmp$Model[1]
best <- models[[best_name]]
cat("\nModel TERBAIK (signifikan & AIC terendah):", best_name, "\n")
cat("AIC =", round(cmp$AIC[1], 3), "| p-value model =", round(cmp$Wald_p[1], 5), "\n")
} else {
cat("\nTidak ada model signifikan (p < 0.05)\n")
}##
## Model TERBAIK (signifikan & AIC terendah): SDEM
## AIC = 597.43 | p-value model = 0.00499## Uji Autokorelasi Spasial Model Terbaik
sdem <- errorsarlm(Y ~ X1 + X2 + X3 + X4,
data=dat_sf, listw=lw, etype="emixed",
method="eigen", zero.policy=TRUE)
res_sdem <- residuals (sdem)
# Moran's I & Geary's C (global)
cat("\n--- Moran's I (Y) ---\n"); print(moran.test(res_sdem, lw, zero.policy=TRUE))##
## --- Moran's I (Y) ---##
## Moran I test under randomisation
##
## data: res_sdem
## weights: lw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.11792, p-value = 0.5469
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04891648 -0.03125000 0.02244706cat("\n--- Geary's C (Y) ---\n"); print(geary.test(res_sdem, lw, zero.policy=TRUE))##
## --- Geary's C (Y) ---##
## Geary C test under randomisation
##
## data: res_sdem
## weights: lw
## n reduced by no-neighbour observations
##
## Geary C statistic standard deviate = -2.0711, p-value = 0.9808
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 1.33773251 1.00000000 0.02659279## Peta Prediksi dan Residual
# =========== PETA PREDIKSI & RESIDUAL (MODEL TERBAIK) =========
if (best_name == "OLS") {
dat_sf$Y_pred <- as.numeric(predict(best, newdata=df_nogeo))
} else {
# fitted.values utk sarlm
dat_sf$Y_pred <- as.numeric(fitted.values(best))
}## This method assumes the response is known - see manual pagedat_sf$Residual <- dat_sf$Y - dat_sf$Y_pred
# Prediksi (kuantil)
br_p <- classInt::classIntervals(dat_sf$Y_pred, n=7, style="quantile")$brks
dat_sf$cutPred <- cut(dat_sf$Y_pred, breaks=br_p, include.lowest=TRUE)
p_pred <- ggplot(dat_sf) +
geom_sf(aes(fill=cutPred), color="white", size=0.25) +
scale_fill_manual(values=blue_pal, name="Kuantil Prediksi") +
labs(title=paste0("Peta Prediksi TBC — Model ", best_name))
print(p_pred); save_png(p_pred, "06_map_prediksi.png")
## Saved: C:/Users/hp/Documents/New folder/SMT 5/SPASIAL/06_map_prediksi.png# Residual
p_res <- ggplot(dat_sf) +
geom_sf(aes(fill=Residual), color="white", size=0.25) +
scale_fill_gradient2(low="#084594", mid="#F7FBFF", high="#2171B5",
midpoint=0, name="Residual") +
labs(title=paste0("Peta Residual — Model ", best_name))
print(p_res)
# Plot perbandingan
p_lin <- ggplot(dat_sf, aes(x = Y_pred, y = Y)) +
geom_point(color = "#1A75FF", size = 2, alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, color = "red", linewidth = 0.8) +
geom_abline(slope = 1, intercept = 0, color = "black", linetype = "dashed") +
labs(
title = paste0("Plot Prediksi vs Aktual ", best_name),
x = "Prediksi Kasus TBC",
y = "Aktual Kasus TBC"
) +
theme_minimal()
print(p_lin)## `geom_smooth()` using formula = 'y ~ x'
## Ringkasan Output
# ======================== OUTPUT ANALISIS ======================
cat("\n================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================\n")##
## ================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================print(cmp)## Model AIC BIC LogLik Wald_p
## Wald statistic3 SDEM 597.43 615.4434 -287.715 0.004989006cat("\n================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================\n")##
## ================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================dat_no_geom <- dat_sf %>% st_set_geometry(NULL)
print(dplyr::glimpse(dat_no_geom))## Rows: 38
## Columns: 18
## $ REGION_KEY <chr> "KABUPATEN: BANGKALAN", "KABUPATEN: BANYUWANGI", "KABUPATE…
## $ REGION_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ P_TYPE <chr> "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUP…
## $ P_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ Y <dbl> 1770, 2715, 1101, 2250, 1521, 3440, 4771, 2263, 2349, 3225…
## $ X1 <dbl> 847, 488, 724, 573, 510, 1086, 786, 1228, 1109, 786, 638, …
## $ X2 <dbl> 4, 13, 8, 10, 3, 19, 14, 14, 9, 19, 8, 3, 3, 24, 11, 4, 5,…
## $ X3 <dbl> 4.057, 3.121, 2.824, 3.231, 3.560, 2.083, 3.186, 2.156, 4.…
## $ X4 <dbl> 190.94, 106.61, 95.91, 147.33, 99.62, 142.39, 224.77, 110.…
## $ prev <dbl> 160.54422, 154.75376, 87.12511, 169.77288, 191.97274, 252.…
## $ cutY <fct> "(1.45e+03,2.03e+03]", "(2.27e+03,3.21e+03]", "(715,1.11e+…
## $ Ii <dbl> 0.210313945, 0.049590395, 0.776928630, 0.088984879, 0.0214…
## $ Pvalue <dbl> 0.59429929, 0.83050575, 0.19181716, 0.48737496, 0.68333028…
## $ Cluster <fct> Not Significant, Not Significant, Not Significant, Not Sig…
## $ Gi_star <dbl> -0.5326162, -0.2140530, -1.3052223, -0.6944900, -0.4079229…
## $ Y_pred <dbl> 1221.2241, 2449.5151, 1495.4414, 2295.8641, 1271.5199, 429…
## $ Residual <dbl> 548.77591, 265.48488, -394.44142, -45.86406, 249.48005, -8…
## $ cutPred <fct> "(607,1.22e+03]", "(2.33e+03,3.26e+03]", "(1.36e+03,1.91e+…
## # A tibble: 38 × 18
## REGION_KEY REGION_NAME P_TYPE P_NAME Y X1 X2 X3 X4 prev
## * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 KABUPATEN: BAN… BANGKALAN KABUP… BANGK… 1770 847 4 4.06 191. 161.
## 2 KABUPATEN: BAN… BANYUWANGI KABUP… BANYU… 2715 488 13 3.12 107. 155.
## 3 KABUPATEN: BLI… BLITAR KABUP… BLITAR 1101 724 8 2.82 95.9 87.1
## 4 KABUPATEN: BOJ… BOJONEGORO KABUP… BOJON… 2250 573 10 3.23 147. 170.
## 5 KABUPATEN: BON… BONDOWOSO KABUP… BONDO… 1521 510 3 3.56 99.6 192.
## 6 KABUPATEN: GRE… GRESIK KABUP… GRESIK 3440 1086 19 2.08 142. 252.
## 7 KABUPATEN: JEM… JEMBER KABUP… JEMBER 4771 786 14 3.19 225. 183.
## 8 KABUPATEN: JOM… JOMBANG KABUP… JOMBA… 2263 1228 14 2.16 111. 166.
## 9 KABUPATEN: KED… KEDIRI KABUP… KEDIRI 2349 1109 9 4.24 159. 139.
## 10 KABUPATEN: LAM… LAMONGAN KABUP… LAMON… 3225 786 19 3.01 147. 234.
## # ℹ 28 more rows
## # ℹ 8 more variables: cutY <fct>, Ii <dbl>, Pvalue <dbl>, Cluster <fct>,
## # Gi_star <dbl>, Y_pred <dbl>, Residual <dbl>, cutPred <fct>cat("\n--- 15 baris pertama (kolom kunci) ---\n")##
## --- 15 baris pertama (kolom kunci) ---print(
dat_no_geom %>%
dplyr::select(REGION_NAME, Y, dplyr::starts_with("X"), Y_pred, Residual) %>%
head(15)
)## # A tibble: 15 × 8
## REGION_NAME Y X1 X2 X3 X4 Y_pred Residual
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BANGKALAN 1770 847 4 4.06 191. 1221. 549.
## 2 BANYUWANGI 2715 488 13 3.12 107. 2450. 265.
## 3 BLITAR 1101 724 8 2.82 95.9 1495. -394.
## 4 BOJONEGORO 2250 573 10 3.23 147. 2296. -45.9
## 5 BONDOWOSO 1521 510 3 3.56 99.6 1272. 249.
## 6 GRESIK 3440 1086 19 2.08 142. 4298. -858.
## 7 JEMBER 4771 786 14 3.19 225. 4122. 649.
## 8 JOMBANG 2263 1228 14 2.16 111. 2543. -280.
## 9 KEDIRI 2349 1109 9 4.24 159. 2412. -62.9
## 10 LAMONGAN 3225 786 19 3.01 147. 3543. -318.
## 11 LUMAJANG 2154 638 8 1.78 91.0 1906. 248.
## 12 MADIUN 972 681 3 4.72 73.2 569. 403.
## 13 MAGETAN 935 970 3 2.41 59.5 209. 726.
## 14 MALANG 3177 788 24 4.37 240. 4653. -1476.
## 15 MOJOKERTO 2012 1172 11 2.02 109. 2296. -284.cat("\n====================== PROPORSI KLASTER LISA ===========================\n")##
## ====================== PROPORSI KLASTER LISA ===========================prop_lisa <- round(prop.table(table(dat_sf$Cluster)), 3)
print(prop_lisa)##
## High-High Low-Low High-Low Low-High Not Significant
## 0.000 0.053 0.000 0.000 0.947# ======================== OUTPUT ANALISIS ======================
cat("\n================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================\n")##
## ================ PERBANDINGAN MODEL (AIC, BIC, LogLik) ================print(cmp)## Model AIC BIC LogLik Wald_p
## Wald statistic3 SDEM 597.43 615.4434 -287.715 0.004989006cat("\n================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================\n")##
## ================= DATA GABUNGAN + PREDIKSI/RESIDUAL ====================dat_no_geom <- dat_sf %>% st_set_geometry(NULL)
print(dplyr::glimpse(dat_no_geom))## Rows: 38
## Columns: 18
## $ REGION_KEY <chr> "KABUPATEN: BANGKALAN", "KABUPATEN: BANYUWANGI", "KABUPATE…
## $ REGION_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ P_TYPE <chr> "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUPATEN", "KABUP…
## $ P_NAME <chr> "BANGKALAN", "BANYUWANGI", "BLITAR", "BOJONEGORO", "BONDOW…
## $ Y <dbl> 1770, 2715, 1101, 2250, 1521, 3440, 4771, 2263, 2349, 3225…
## $ X1 <dbl> 847, 488, 724, 573, 510, 1086, 786, 1228, 1109, 786, 638, …
## $ X2 <dbl> 4, 13, 8, 10, 3, 19, 14, 14, 9, 19, 8, 3, 3, 24, 11, 4, 5,…
## $ X3 <dbl> 4.057, 3.121, 2.824, 3.231, 3.560, 2.083, 3.186, 2.156, 4.…
## $ X4 <dbl> 190.94, 106.61, 95.91, 147.33, 99.62, 142.39, 224.77, 110.…
## $ prev <dbl> 160.54422, 154.75376, 87.12511, 169.77288, 191.97274, 252.…
## $ cutY <fct> "(1.45e+03,2.03e+03]", "(2.27e+03,3.21e+03]", "(715,1.11e+…
## $ Ii <dbl> 0.210313945, 0.049590395, 0.776928630, 0.088984879, 0.0214…
## $ Pvalue <dbl> 0.59429929, 0.83050575, 0.19181716, 0.48737496, 0.68333028…
## $ Cluster <fct> Not Significant, Not Significant, Not Significant, Not Sig…
## $ Gi_star <dbl> -0.5326162, -0.2140530, -1.3052223, -0.6944900, -0.4079229…
## $ Y_pred <dbl> 1221.2241, 2449.5151, 1495.4414, 2295.8641, 1271.5199, 429…
## $ Residual <dbl> 548.77591, 265.48488, -394.44142, -45.86406, 249.48005, -8…
## $ cutPred <fct> "(607,1.22e+03]", "(2.33e+03,3.26e+03]", "(1.36e+03,1.91e+…
## # A tibble: 38 × 18
## REGION_KEY REGION_NAME P_TYPE P_NAME Y X1 X2 X3 X4 prev
## * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 KABUPATEN: BAN… BANGKALAN KABUP… BANGK… 1770 847 4 4.06 191. 161.
## 2 KABUPATEN: BAN… BANYUWANGI KABUP… BANYU… 2715 488 13 3.12 107. 155.
## 3 KABUPATEN: BLI… BLITAR KABUP… BLITAR 1101 724 8 2.82 95.9 87.1
## 4 KABUPATEN: BOJ… BOJONEGORO KABUP… BOJON… 2250 573 10 3.23 147. 170.
## 5 KABUPATEN: BON… BONDOWOSO KABUP… BONDO… 1521 510 3 3.56 99.6 192.
## 6 KABUPATEN: GRE… GRESIK KABUP… GRESIK 3440 1086 19 2.08 142. 252.
## 7 KABUPATEN: JEM… JEMBER KABUP… JEMBER 4771 786 14 3.19 225. 183.
## 8 KABUPATEN: JOM… JOMBANG KABUP… JOMBA… 2263 1228 14 2.16 111. 166.
## 9 KABUPATEN: KED… KEDIRI KABUP… KEDIRI 2349 1109 9 4.24 159. 139.
## 10 KABUPATEN: LAM… LAMONGAN KABUP… LAMON… 3225 786 19 3.01 147. 234.
## # ℹ 28 more rows
## # ℹ 8 more variables: cutY <fct>, Ii <dbl>, Pvalue <dbl>, Cluster <fct>,
## # Gi_star <dbl>, Y_pred <dbl>, Residual <dbl>, cutPred <fct>cat("\n--- 15 baris pertama (kolom kunci) ---\n")##
## --- 15 baris pertama (kolom kunci) ---print(
dat_no_geom %>%
dplyr::select(REGION_NAME, Y, dplyr::starts_with("X"), Y_pred, Residual) %>%
head(15)
)## # A tibble: 15 × 8
## REGION_NAME Y X1 X2 X3 X4 Y_pred Residual
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BANGKALAN 1770 847 4 4.06 191. 1221. 549.
## 2 BANYUWANGI 2715 488 13 3.12 107. 2450. 265.
## 3 BLITAR 1101 724 8 2.82 95.9 1495. -394.
## 4 BOJONEGORO 2250 573 10 3.23 147. 2296. -45.9
## 5 BONDOWOSO 1521 510 3 3.56 99.6 1272. 249.
## 6 GRESIK 3440 1086 19 2.08 142. 4298. -858.
## 7 JEMBER 4771 786 14 3.19 225. 4122. 649.
## 8 JOMBANG 2263 1228 14 2.16 111. 2543. -280.
## 9 KEDIRI 2349 1109 9 4.24 159. 2412. -62.9
## 10 LAMONGAN 3225 786 19 3.01 147. 3543. -318.
## 11 LUMAJANG 2154 638 8 1.78 91.0 1906. 248.
## 12 MADIUN 972 681 3 4.72 73.2 569. 403.
## 13 MAGETAN 935 970 3 2.41 59.5 209. 726.
## 14 MALANG 3177 788 24 4.37 240. 4653. -1476.
## 15 MOJOKERTO 2012 1172 11 2.02 109. 2296. -284.cat("\n====================== PROPORSI KLASTER LISA ===========================\n")##
## ====================== PROPORSI KLASTER LISA ===========================prop_lisa <- round(prop.table(table(dat_sf$Cluster)), 3)
print(prop_lisa)##
## High-High Low-Low High-Low Low-High Not Significant
## 0.000 0.053 0.000 0.000 0.947