epidemmmmmmmmmm.knit

Studi Epidemiologi: Analisis Kasus Campak di Kabupaten/Kota
Provinsi Jawa Timur Tahun 2022

Ditujukan Guna Memenuhi UAS Mata Kuliah Epidemiologi

Disusun oleh:
Tiffanny Chantika Silalahi (140610230027)

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

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

ABSTRAK

Campak tetap menjadi tantangan kesehatan publik yang serius di Indonesia, khususnya di Jawa Timur. Studi ini mengevaluasi karakteristik epidemiologi kasus campak di Jawa Timur selama tahun 2022 dengan fokus pada distribusi kasus, indikator epidemiologi, dan pola penyebaran geografis. Analisis ini menggunakan data agregat tingkat kabupaten/kota yang mencakup jumlah kasus campak dan data populasi. Pendekatan epidemiologi meliputi pemetaan distribusi kasus, penghitungan incidence rate per 100.000 populasi dan prevalensi, pengujian klasterisasi spasial menggunakan indeks Moran’s I, serta identifikasi hotspot lokal melalui Local Indicator of Spatial Autocorrelation (LISA). Kami juga menghitung Standardized Morbidity Ratio (SMR) dengan confidence interval 95% dan menerapkan metode Empirical Bayes untuk memperbaiki estimasi risiko di wilayah berpopulasi kecil. Kejadian campak di Jawa Timur menunjukkan distribusi yang heterogen antarwilayah. Indeks Moran’s I yang signifikan mengonfirmasi adanya klasterisasi spasial positif, yang diperkuat oleh identifikasi zona risiko tinggi melalui analisis LISA. SMR mengungkap variasi risiko geografis yang substansial, sementara koreksi Empirical Bayes berhasil meredam volatilitas estimasi risiko. Penyebaran campak di Jawa Timur menunjukkan pola epidemiologi yang terstruktur secara spasial, bukan distribusi acak. Temuan ini memberikan landasan empiris untuk merancang intervensi kesehatan masyarakat yang berbasis geografis dan lebih efektif.

Kata kunci: campak; epidemiologi spasial; incidence rate; prevalensi; autokorelasi spasial; Moran’s I; LISA; Standardized Morbidity Ratio; Empirical Bayes

BAB I

PENDAHULUAN

Latar Belakang

Epidemiologi merupakan disiplin ilmu kesehatan masyarakat yang mempelajari distribusi dan determinan penyakit dalam populasi untuk mendukung pencegahan serta pengendalian penyakit secara efektif[1]. Campak masih menjadi masalah kesehatan global dengan angka kejadian dan kematian yang meningkat secara signifikan pada periode 2021–2022 [2], menunjukkan tantangan terkini dalam pengendalian penyakit ini di banyak negara berkembang.

Campak dikenal memiliki tingkat penularan yang sangat tinggi—dengan potensi satu kasus dapat menulari banyak individu lain, dan penularan ini sangat dipengaruhi oleh tingkat kekebalan kelompok (herd immunity)[3]. Cakupan imunisasi yang turun di berbagai wilayah telah dikaitkan dengan munculnya kembali kejadian luar biasa campak, sehingga pemahaman epidemiologis terhadap distribusi kasus sangat penting untuk intervensi kesehatan masyarakat[3].

Dalam kajian epidemiologi, desain studi potong-lintang (cross-sectional) digunakan untuk menggambarkan sebaran penyakit dan mengukur ukuran frekuensi seperti insidensi dan prevalensi pada satu titik waktu tertentu [4], sehingga memungkinkan perbandingan risiko antar wilayah tanpa menilai hubungan kausal temporal. Desain ini sangat sesuai untuk menilai variasi epidemiologi penyakit seperti campak dalam satu tahun pengamatan.

Rumusan Masalah

Bagaimana distribusi angka insidensi campak di kabupaten/kota Provinsi Jawa Timur tahun 2022?
Apakah terdapat autokorelasi spasial kasus campak antar wilayah?
Wilayah mana yang berpotensi menjadi hotspot atau coldspot campak?
Bagaimana estimasi risiko relatif campak setelah dilakukan stabilisasi menggunakan Empirical Bayes?

Tujuan

Penelitian ini bertujuan untuk menganalisis karakteristik epidemiologi penyakit campak di Provinsi Jawa Timur tahun 2022 melalui perhitungan ukuran epidemiologi dasar dan estimasi risiko wilayah.

BAB II

TINJAUAN PUSTAKA

Campak dalam Perspektif Epidemiologi

Campak merupakan penyakit infeksi virus yang disebabkan oleh Measles virus dari genus Morbillivirus [5]. Penyakit ini memiliki masa inkubasi sekitar 10–14 hari dan ditandai oleh gejala demam, batuk, pilek, konjungtivitis, serta ruam makulopapular [6]. Campak dapat menimbulkan komplikasi serius seperti pneumonia dan ensefalitis, terutama pada anak-anak dengan status gizi buruk atau imunitas rendah [7].

Dalam epidemiologi, kejadian penyakit campak dipahami sebagai hasil interaksi antara agent, host, dan environment. Kerangka ini dikenal sebagai segitiga epidemiologi dan digunakan secara luas untuk menjelaskan dinamika penyakit menular.

Agent-Host-Environment

Agent

Agent penyebab campak adalah virus campak yang memiliki tingkat infektivitas sangat tinggi. Virus ini mudah ditularkan melalui udara dan droplet, sehingga kontak singkat dengan individu terinfeksi dapat menyebabkan penularan [8]. Tidak adanya reservoir non-manusia menjadikan manusia sebagai satu-satunya host alami, sehingga eliminasi penyakit sangat bergantung pada tingkat kekebalan populasi melalui imunisasi[9].

Host

Host dalam epidemiologi campak adalah manusia yang rentan terhadap infeksi. Faktor host yang berperan meliputi status imunisasi, usia, status gizi, dan kondisi imunitas. Individu yang belum mendapatkan imunisasi campak atau memiliki kekebalan yang tidak adekuat memiliki risiko lebih tinggi untuk terinfeksi [10]. Dalam penelitian ini, faktor host direpresentasikan secara agregat melalui cakupan imunisasi campak berdasarkan jenis kelamin di tingkat kabupaten/kota.

Environment

Lingkungan berperan dalam memfasilitasi atau menghambat transmisi campak. Kepadatan penduduk, kondisi sanitasi, tingkat kemiskinan, serta akses terhadap pelayanan kesehatan merupakan faktor lingkungan yang dapat memengaruhi peluang penularan [11]. Lingkungan dengan kondisi sosial ekonomi yang kurang baik cenderung memiliki risiko penularan yang lebih tinggi [12]. Dalam studi ini, faktor lingkungan direpresentasikan oleh persentase rumah tangga dengan akses sanitasi layak dan persentase penduduk miskin.

Ukuran Epidemiologi

Ukuran epidemiologi digunakan untuk menggambarkan frekuensi dan risiko penyakit dalam populasi. Angka insidensi menggambarkan jumlah kasus baru dalam suatu periode tertentu [13], sedangkan Standardized Morbidity Ratio digunakan untuk membandingkan kejadian penyakit antar wilayah dengan memperhitungkan perbedaan ukuran populasi [14]. Pendekatan Empirical Bayes digunakan untuk menstabilkan estimasi risiko pada wilayah dengan jumlah penduduk kecil.

Desain Studi

Penelitian ini menggunakan desain studi potong lintang (cross-sectional) berbasis wilayah. Pada desain ini, pengukuran paparan dan kejadian penyakit dilakukan pada satu periode waktu yang sama, yaitu tahun 2022. Desain cross-sectional sesuai untuk menggambarkan distribusi penyakit dan membandingkan risiko antar wilayah, namun tidak dapat digunakan untuk menilai hubungan sebab-akibat secara temporal seperti pada studi kohort [16].

BAB III

METODOLOGI PENELITIAN

Sumber dan Jenis Data

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

Variabel Penelitian

Variabel utama dalam penelitian ini adalah jumlah kasus campak, yang merepresentasikan kejadian penyakit. Jumlah penduduk digunakan sebagai penyebut dalam perhitungan ukuran frekuensi penyakit. Variabel host meliputi persentase penduduk yang mendapatkan imunisasi campak pada laki-laki dan perempuan. Variabel lingkungan meliputi persentase penduduk miskin dan persentase rumah tangga dengan akses sanitasi layak.

Ukuran Epidemiologi dan Analisis

Incidense Rate (IR)

Angka insidensi digunakan untuk menggambarkan jumlah kasus baru campak dalam populasi berisiko pada periode tertentu.

Angka insidensi didefinisikan sebagai jumlah kasus baru campak dalam populasi berisiko pada periode tertentu:

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

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

Prevalensi

Prevalensi digunakan untuk menggambarkan proporsi penduduk yang mengalami campak pada tahun pengamatan.

Prevalensi menunjukkan proporsi penduduk yang mengalami campak pada tahun pengamatan:

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

Standardized Morbidity Ratio (SMR)

SMR digunakan untuk membandingkan kejadian campak di suatu wilayah terhadap kejadian yang diharapkan berdasarkan rata-rata populasi.

Standardized Morbidity Ratio (SMR) digunakan untuk membandingkan jumlah kasus teramati dengan jumlah kasus yang diharapkan:

\[ \text{SMR}_i = \frac{O_i}{E_i} \]

dengan: - $O_i$ : jumlah kasus campak teramati di wilayah ke-$i$ - $E_i$ : jumlah kasus campak yang diharapkan di wilayah ke-$i$

Empirical Bayes Estimation

Empirical Bayes Estimation digunakan untuk menstabilkan nilai risiko wilayah dengan jumlah penduduk kecil dengan mengombinasikan SMR lokal dan rata-rata global.

Empirical Bayes Estimation digunakan untuk menstabilkan estimasi risiko wilayah:

\[ \theta_i^{EB} = w_i \cdot \theta_i + (1 - w_i) \cdot \bar{\theta} \]

dengan: - $\theta_i$ : SMR wilayah ke-$i$ - $\bar{\theta}$ : rata-rata SMR seluruh wilayah - $w_i$ : bobot keandalan wilayah ke-$i$, $0 \leq w_i \leq 1$

BAB IV

HASIL DAN PEMBAHASAN

4.1 Peta Prevalensi Kasus Campak di Jawa Timur 2022

options(stringsAsFactors = FALSE)

library(sf)
library(dplyr)
library(readxl)
library(ggplot2)
library(spdep)
library(spatialreg)
library(gstat)
library(sp)
library(spgwr)
library(car)
library(lmtest)
library(viridis)
library(DT)


## PEMANGGILAN DATA 
# === GANTI PATH SESUAI LOKASI FILE ===
shp_path   <- "C:/Users/ACER/Downloads/Spasial_Epidem_UAS/GADM_JATIM"
excel_path <- "C:/Users/ACER/Downloads/Spasial_Epidem_UAS/epidem_spasial_campak_2022.xlsx"

shp <- st_read(shp_path, quiet = TRUE) %>%
  filter(tolower(NAME_1) == "jawa timur") %>%
  st_make_valid()

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

df <- read_excel(excel_path)
names(df) <- trimws(gsub("\\s+", "_", names(df)))


## 🔗 JOIN DATA & PEMBENTUKAN VARIABEL
clean_name <- function(x) tolower(gsub("\\s+|kabupaten|kota|,|\\.", "", x))

shp <- shp %>% mutate(key = clean_name(NAME_2))
df  <- df  %>% mutate(key = clean_name(Kabupaten_Kota)) %>%
  filter(Tahun == 2022)

dat <- left_join(shp, df, by = "key") %>%
  mutate(
    IR_100k = Kasus_Campak / Penduduk * 100000
  )

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

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

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

Peta persebaran kasus menunjukkan bahwa kejadian campak di Provinsi Jawa Timur tahun 2022 tidak tersebar merata antar kabupaten/kota. Beberapa wilayah terlihat memiliki jumlah kasus yang jauh lebih tinggi dibandingkan wilayah lain. Pola ini mengindikasikan bahwa risiko campak sangat dipengaruhi oleh kondisi lokal, seperti kepadatan penduduk, cakupan imunisasi, dan faktor lingkungan setempat. Wilayah dengan warna lebih gelap pada peta menandakan prevalensi kasus yang lebih tinggi, sehingga dapat dianggap sebagai daerah prioritas dalam pengendalian penyakit campak.

4.2 Ukuran Epidemiologi (Prevalensi)

library(knitr)

## Warning: package 'knitr' was built under R version 4.3.3

library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.3.3

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

##  Incidence rate (100.000)
dat %>%
  st_drop_geometry() %>%
  select(NAME_2, Kasus_Campak, Penduduk, IR_100k) %>%
  datatable()

##  Prevalensi per Kabupaten/Kota

dat %>%
  st_drop_geometry() %>%
  mutate(Prevalensi = Kasus_Campak / Penduduk) %>%
  select(NAME_2, Prevalensi) %>%
  datatable()

Berdasarkan tabel ukuran epidemiologi, terlihat bahwa kabupaten/kota dengan jumlah kasus tinggi tidak selalu memiliki risiko tertinggi. Incidence Rate (IR per 100.000 penduduk) dan prevalensi memberikan gambaran yang lebih jelas karena mempertimbangkan jumlah penduduk.

Kabupaten Pamekasan memiliki IR tertinggi yaitu sebesar 3.59. Hal ini menunjukkan bahwa risiko campak di Pamekasan relatif lebih besar dibandingkan wilayah lain, bukan hanya karena jumlah penduduknya. Prevalensi per kabupaten/kota juga menunjukkan Kabupaten Pamekasa memiliki nilai tertinggi, yakni sebesar 0.000036. hal ini menunjukkan proporsi penduduk yang terdampak campak dalam satu periode waktu tertentu.

4.3 Auto Korelasi Spasial

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

Hasil uji Moran’s I menunjukkan nilai statistik sebesar 0,694 dengan p-value < 0,001, yang berarti terdapat autokorelasi spasial positif yang sangat signifikan. Artinya, wilayah dengan IR campak tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki IR tinggi, dan sebaliknya untuk wilayah dengan IR rendah.

Peta LISA menunjukkan adanya klaster High–High, terutama di wilayah utara Jawa Timur. Klaster ini berarti bahwa kabupaten/kota tersebut memiliki nilai IR tinggi dan dikelilingi oleh wilayah lain yang juga memiliki IR tinggi. Kondisi ini mengindikasikan adanya hotspot epidemiologis yang berpotensi menjadi pusat penularan campak. Wilayah lain sebagian besar berada pada kategori tidak signifikan, yang menunjukkan bahwa tidak semua daerah dengan IR rendah atau tinggi membentuk klaster yang bermakna secara statistik.

4.4 Standardized Morbidity Ratio (SMR)

rbar <- sum(dat$Kasus_Campak)/sum(dat$Penduduk)

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

datatable(dat_smr %>% select(NAME_2, O, E, SMR, LCL, UCL))

Hasil SMR menunjukkan adanya perbedaan besar antara jumlah kasus yang teramati (O) dan yang diharapkan (E) di beberapa wilayah. Kabupaten Pamekasan, memiliki nilai SMR lebih dari 3 dengan batas bawah dan atas interval kepercayaan 95% yang berada di atas 1. Hal ini menandakan bahwa kejadian campak di wilayah tersebut secara signifikan lebih tinggi dari yang diharapkan berdasarkan rata-rata provinsi.

Sebaliknya, beberapa wilayah memiliki SMR jauh di bawah 1, yang menunjukkan kejadian campak lebih rendah dari ekspektasi. SMR membantu mengidentifikasi daerah dengan risiko relatif tinggi atau rendah, sehingga sangat berguna dalam pemantauan dan evaluasi program pengendalian penyakit.

4.5 Empirical Bayes (EB)

dat_smr$EB <- (1 + dat_smr$O)/(1 + dat_smr$E)

ggplot(dat_smr, aes(SMR, EB)) +
  geom_point() +
  geom_abline(linetype=2) +
  theme_minimal()

Hasil scatter plot antara SMR dan estimasi Empirical Bayes menunjukkan bahwa metode EB menarik nilai SMR ekstrem mendekati rata-rata. Wilayah dengan populasi kecil dan SMR sangat tinggi atau sangat rendah mengalami penyesuaian yang lebih besar dibandingkan wilayah berpopulasi besar.

Pendekatan ini penting karena SMR mentah pada wilayah dengan jumlah penduduk kecil cenderung tidak stabil dan mudah dipengaruhi oleh fluktuasi acak. Dengan EB smoothing, estimasi risiko menjadi lebih stabil dan lebih dapat diandalkan untuk pengambilan keputusan kebijakan kesehatan.

BAB V

KESIMPULAN DAN SARAN

Kesimpulan

Berdasarkan hasil analisis epidemiologi campak di Provinsi Jawa Timur tahun 2022, dapat disimpulkan bahwa kejadian campak tidak tersebar secara merata antar kabupaten/kota, melainkan membentuk pola spasial yang jelas. Beberapa wilayah menunjukkan konsentrasi kasus dan tingkat risiko yang lebih tinggi dibandingkan wilayah lainnya, sehingga kejadian campak tidak dapat dianggap terjadi secara acak.

Penggunaan ukuran epidemiologi berbasis penduduk, seperti incidence rate dan prevalensi, memberikan gambaran risiko yang lebih akurat dibandingkan hanya melihat jumlah kasus absolut. Wilayah dengan penduduk lebih kecil dapat memiliki risiko yang lebih tinggi, meskipun jumlah kasusnya tidak paling besar. Hal ini menegaskan pentingnya penggunaan indikator epidemiologi yang mempertimbangkan ukuran populasi dalam penilaian risiko penyakit.

Hasil uji Moran’s I menunjukkan adanya autokorelasi spasial positif yang signifikan, yang berarti wilayah dengan tingkat kejadian campak tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki tingkat kejadian tinggi. Temuan ini diperkuat oleh analisis LISA, yang mengidentifikasi klaster wilayah berisiko tinggi (High–High), terutama di wilayah Pulau Madura. Klaster ini menunjukkan adanya daerah yang berpotensi menjadi pusat penularan campak.

Analisis Standardized Morbidity Ratio (SMR) mengungkapkan bahwa beberapa kabupaten/kota memiliki kejadian campak yang secara signifikan lebih tinggi dari yang diharapkan berdasarkan rata-rata provinsi. Selain itu, penerapan Empirical Bayes smoothing berhasil menstabilkan estimasi risiko, terutama pada wilayah dengan populasi kecil, sehingga memberikan gambaran risiko yang lebih andal untuk keperluan perencanaan kesehatan masyarakat.

Secara keseluruhan, hasil penelitian ini menunjukkan bahwa pendekatan epidemiologi yang dikombinasikan dengan analisis spasial sangat efektif dalam mengidentifikasi wilayah berisiko tinggi dan memahami pola penyebaran campak di tingkat kabupaten/kota.

Saran

Berdasarkan hasil dan kesimpulan yang diperoleh, disarankan agar intervensi pengendalian campak difokuskan pada wilayah yang teridentifikasi sebagai klaster risiko tinggi, terutama daerah dengan nilai incidence rate, prevalensi, dan SMR yang tinggi. Pendekatan berbasis wilayah ini diharapkan lebih efektif dibandingkan intervensi yang dilakukan secara merata di seluruh provinsi.

Dinas kesehatan perlu meningkatkan cakupan dan pemerataan imunisasi campak, khususnya di wilayah yang menunjukkan konsentrasi kasus tinggi dan risiko relatif yang signifikan. Upaya ini dapat dikombinasikan dengan edukasi masyarakat, peningkatan surveilans penyakit, serta perbaikan akses layanan kesehatan di daerah tersebut.

Penggunaan metode Empirical Bayes disarankan untuk digunakan secara rutin dalam analisis penyakit berbasis wilayah, terutama ketika terdapat perbedaan besar jumlah penduduk antar kabupaten/kota. Metode ini membantu mengurangi bias akibat fluktuasi acak dan menghasilkan estimasi risiko yang lebih stabil.

Untuk penelitian selanjutnya, disarankan agar analisis epidemiologi ini dikembangkan dengan memasukkan variabel penjelas tambahan, seperti cakupan imunisasi, kepadatan penduduk, dan kondisi sosial ekonomi, serta dikombinasikan dengan model regresi spasial dan GWR. Pendekatan tersebut akan memberikan pemahaman yang lebih mendalam mengenai faktor-faktor yang memengaruhi penyebaran campak dan mendukung perumusan kebijakan kesehatan berbasis bukti yang lebih kuat.

Daftar Pustaka

[1] Penulis, A. A. (Tahun). Judul buku: Subjudul buku (Edisi jika ada). Springer Publishing. https://connect.springerpub.com/content/book/978-0-8261-6694-4

[2] Minta, A. A., Ferrari, M., Antoni, S., et al. (2023). Progress toward measles elimination—Worldwide, 2000–2022. MMWR Morbidity and Mortality Weekly Report, 72(46), 1262–1268. https://doi.org/10.15585/mmwr.mm7246a3

[3] World Health Organization. Measles [Internet]. Geneva: WHO; 2023 [cited 2024 May 22]. Available from: https://www.who.int/news-room/fact-sheets/detail/measles

[4] Zuleika, P., & Legiran. (2022). Cross-sectional study as research design in medicine. Archives of The Medicine and Case Reports, 3(2), 256–259. https://doi.org/10.37275/amcr.v3i2.193

[5] Nwalozie, R., Uzoechi, M., Esiere, R. K., & Nnokam, B. A. (2023). Biology of measles virus: Epidemiology and clinical manifestations. International Journal of Pathogen Research, 12(4), 1–10. https://doi.org/10.9734/ijpr/2023/v12i4231

[6] Misin, A., Antonello, R. M., Di Bella, S., Campisciano, G., Zanotta, N., Giacobbe, D. R., Comar, M., & Luzzati, R. (2020). Measles: An overview of a re-emerging disease in children and immunocompromised patients. Microorganisms, 8(2), 276. https://doi.org/10.3390/microorganisms8020276

[7] Perry, R. T., & Halsey, N. A. (2004). The clinical significance of measles: A review. Journal of Infectious Diseases, 189(Supplement_1), S4–S16. https://doi.org/10.1086/377712

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

[9] Centers for Disease Control and Prevention. (2024). Measles (rubeola): Epidemiology and prevention of vaccine-preventable diseases (The Pink Book, Chapter 13). U.S. Department of Health and Human Services. https://www.cdc.gov/pinkbook/hcp/table-of-contents/chapter-13-measles.html

[10] World Health Organization. (2023). Measles fact sheet. World Health Organization. https://www.who.int/news-room/fact-sheets/detail/measles

[11] World Health Organization. (2020). Framework for strengthening the immunization and surveillance system for measles and rubella elimination. WHO. https://apps.who.int/iris/handle/10665/332240

[12] Cutts, F. T., Ferrari, M. J., Krause, L. K., Tatem, A. J., & Mosser, J. F. (2021). Vaccination strategies for measles control and elimination: Time to strengthen local initiatives. BMC Medicine, 19(1), 2–9. https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-020-01836-6

[13] Centers for Disease Control and Prevention. (2022). Principles of epidemiology in public health practice (Lesson 3). CDC. https://www.cdc.gov/csels/dsepd/ss1978/lesson3/section2.html

[14] Lawson, A. B. (2021). Bayesian disease mapping: Hierarchical modeling in spatial epidemiology (2nd ed.). CRC Press. https://www.taylorfrancis.com/books/oa-mono/10.1201/9780429056069

[15] Martuzzi, M., & Elliott, P. (1996). Empirical Bayes estimation of small area prevalence of non-rare conditions. Statistics in Medicine, 15(17–18), 1867–1873. https://onlinelibrary.wiley.com/doi/10.1002/(SICI)1097-0258(19960930)15:17/18 <1867::AID-SIM362>3.0.CO;2-A

[16] Levin, K. A. (2006). Study design III: Cross-sectional studies. Evidence-Based Dentistry, 7(1), 24–25. https://www.nature.com/articles/6400375

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

Lampiran

Lampiran Script Analisis R

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

📊 Dashboard Interaktif Analisis Spasial Campak (Jawa Timur) 2022

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

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

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

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

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

🧭 THEME

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

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

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

🔧 Helper

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

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

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

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

🧱 UI

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

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

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

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

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

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

⚙️ SERVER

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

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

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

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

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

shp

})

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

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

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

})

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

req(shp, df)

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

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

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

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

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

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

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

joined

})

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

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

}

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

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

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

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

})

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

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

})

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

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

})

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

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

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

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

})

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

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

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

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

})

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

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

})

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

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

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

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

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

})

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

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

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

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

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

grid_sp <- as(grid_sf, "Spatial")

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

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

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

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

bb <- st_bbox(sf_use)

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

})

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

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

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

v <- variogram(yvar ~ 1, pts_sp)

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

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

})

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

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

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

as(grid_sf, "Spatial")

}

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

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

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

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

bb <- st_bbox(obj$sf_use)

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

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

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

})

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

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

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

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

bb <- st_bbox(obj$sf_use)

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

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

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

})

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

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

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

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

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

bb <- st_bbox(obj$sf_use)

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

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

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

})

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

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

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

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

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

bb <- st_bbox(obj$sf_use)

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

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

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

})

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

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

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

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

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

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

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

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

})

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

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

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

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

})

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

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

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

})

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

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

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

})

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

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

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

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

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

})

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

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

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

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

})

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

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

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

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

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

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

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

})

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

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

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

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

})

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

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

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

})

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

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

})

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

vif_val <- car::vif(model)

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

})

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

bp <- bptest(model)

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

})

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

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

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

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

lw <- make_listw(sf_use, input$wtype)

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

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

})

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

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

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

})

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

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

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

}

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

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

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

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

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

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

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

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

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

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

})

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

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

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

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

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

res

})

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

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

test <- shapiro.test(res)

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

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

})

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

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

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

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

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

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

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

})

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

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

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

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

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

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

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

})

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

sf_use$pred <- gwr_res$SDF$pred

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

}) }

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

🚀 Run App

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

shinyApp(ui, server)

Lampiran Link Dashboard

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

Ditujukan Guna Memenuhi UAS Mata Kuliah Epidemiologi

Latar Belakang

Rumusan Masalah

Tujuan

Campak dalam Perspektif Epidemiologi

Agent-Host-Environment

Agent

Host

Environment

Ukuran Epidemiologi

Desain Studi

Sumber dan Jenis Data

Variabel Penelitian

Ukuran Epidemiologi dan Analisis

Incidense Rate (IR)

Prevalensi

Standardized Morbidity Ratio (SMR)

Empirical Bayes Estimation

Kesimpulan

Saran

Lampiran Script Analisis R

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

📊 Dashboard Interaktif Analisis Spasial Campak (Jawa Timur) 2022

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

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

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

🧭 THEME

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

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

🔧 Helper

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

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

🧱 UI

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

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

⚙️ SERVER

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

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

🚀 Run App

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

Lampiran Link Dashboard

Lampiran Link Video Social Media