Analisis Epidemiologi Penyakit Kusta dan Faktor Sanitasi di Provinsi Jawa Timur

Fazila Azra Anggina
NPM: 140610230039

ABSTRAK

Penyakit kusta masih menjadi permasalahan kesehatan masyarakat di Provinsi Jawa Timur dengan distribusi kasus yang bervariasi antar kabupaten/kota. Penelitian ini bertujuan untuk menggambarkan kondisi epidemiologi penyakit kusta serta menganalisis faktor sanitasi yang berhubungan dengan kejadian kusta di Provinsi Jawa Timur. Penelitian menggunakan desain cross-sectional ekologis dengan unit analisis 38 kabupaten/kota. Data yang digunakan merupakan data sekunder yang mencakup prevalensi kusta, jumlah kasus baru, jumlah penduduk, persentase rumah tangga yang menggunakan MCK bersama, dan persentase rumah tangga tanpa akses sanitasi layak.

Hasil analisis deskriptif menunjukkan bahwa prevalensi kusta memiliki nilai rata-rata sebesar 0,44 per 1.000 penduduk, dengan nilai tertinggi di Kabupaten Sumenep (2,06) dan beberapa wilayah dengan prevalensi 0, seperti Kota Blitar, Kota Mojokerto, dan Kota Batu. Incidence rate kusta bervariasi antar wilayah, yang mengindikasikan masih terjadinya penularan aktif di beberapa kabupaten/kota. Nilai release from treatment (RFT) umumnya tinggi, dengan sebagian besar kabupaten/kota mencapai nilai mendekati 100%, meskipun variasi lebih besar ditemukan pada kusta tipe multibacillary.

Hasil analisis inferensial menunjukkan bahwa persentase rumah tangga yang menggunakan MCK bersama berhubungan positif dengan jumlah kasus baru kusta. Analisis spasial menggunakan Moran’s I menunjukkan adanya autokorelasi spasial yang signifikan, yang menandakan pengelompokan wilayah dengan prevalensi kusta tinggi. Standardisasi berdasarkan jenis kelamin menunjukkan perbedaan yang kecil antara crude rate dan standardized rate, sehingga struktur jenis kelamin bukan merupakan faktor utama dalam variasi kejadian kusta. Oleh karena itu, dapat disimpulkan bahwa distribusi penyakit kusta di Jawa Timur dipengaruhi oleh faktor lingkungan dan keterkaitan spasial antar wilayah.

---

1 Pendahuluan

1.1 Latar Belakang

Dalam buku rencana aksi SDGs yang dikeluarkan oleh Bappenas pada tahun 2020, salah satu target rencana aksi SDGs bidang kesehatan di Indonesia adalah mengakhiri epidemi Acquired Immune Deficiency Syndrome (AIDS), Tuberkulosis, malaria, filariasis dan kusta (Bappenas, 2020). Salah satu penyakit kulit menular ialah penyakit kusta. Penyakit ini dipicu oleh bakteri Mycobacterium Leprae. Penyakit kusta umumnya menyerang bagian kulit, pernapasan, mata, saraf, dan anggota gerak tubuh.

Indonesia merupakan salah satu negara berkembang yang menyumbang kasus baru penyakit kusta di urutan ke-3 terbesar lingkup dunia. Jumlah kasus yang diciptakan sebesar 8% total penyakit kusta di dunia atau sebanyak 9.601 kasus (Kemenkes, 2021a). Pendeteksian kasus baru penyakit kusta diperlukan untuk mengendalikan peningkatan kasus baru penyakit kusta terjadi. Beberapa strategi yang diperlukan untuk memenuhi target penanggulangan kusta di Indonesia , sebagaimana tertuang pada UU Nomor 11 Tahun 2019.

Provinsi Jawa Timur merupakan satu diantara provinsi lainnya yang terindikasi memiliki beban penyakit kusta tertinggi di Indonesia (Ritianty et al., 2020). Selama tahun 2019 hingga 2020, pengendalian kasus kusta di Jawa Timur mengalami pengingkatan berdasarkan mengecilnya penemuan kusta dari 8,06 menjadi 4,19 per 1.000.000. Untuk mewujudkan eliminasi kusta di Indonesia, khususnya provinsi Jawa Timur, dilakukanlah suatu penelitian untuk membentuk pemodelan statistika terkait prevalensi kasus penyakit kusta di Jawa Timur sehingga faktor yang mempengaruhi jumlah kasus pada penyakit kusta dapat diidentifikasi secara jelas.

1.2 Identifikasi Masalah

Berdasarkan latar belakang tersebut, permasalahan yang dapat diidentifikasi dalam penelitian ini adalah :

1. Distribusi prevalensi penyakit kusta di Jawa Timur cenderung terkonsentrasi pada wilayah tertentu.
2. Faktor lingkungan dan sanitasi diduga berperan dalam tingginya prevalensi penyakit kusta di beberapa daerah.

1.3 Tujuan Penelitian

Tujuan dari penelitian ini adalah :

1. Mendeskripsikan gambaran epidemiologi prevalensi penyakit kusta di Jawa Timur.
2. Mengidentifikasi faktor-faktor yang memengaruhi prevalensi penyakit kusta dengan mempertimbangkan efek spasial.

2 Tinjauan Pustaka

2.1 Agent-Host-Environment

library(DiagrammeR)

grViz("
digraph kusta {
  graph [rankdir = LR]
  A [label='Lingkungan: Sanitasi buruk']
  B [label='Perilaku/Higiene menurun\\n& kepadatan/kontak meningkat']
  C [label='Pajanan M. leprae']
  D [label='Infeksi kusta']
  E [label='Kasus baru meningkat']
  F [label='Prevalensi meningkat']

  A -> B -> C -> D -> E -> F
  M [label='MCK bersama', shape=box]
  M -> B
}
")

Berdasarkan konsep Agent–Host–Environment (AHE), penyakit kusta disebabkan oleh agen biologis berupa Mycobacterium leprae. Faktor host yang berperan meliputi kondisi higiene yang kurang baik serta tingginya intensitas kontak antar individu. Sementara itu, faktor lingkungan yang berkontribusi terhadap penularan kusta antara lain kondisi sanitasi yang buruk dan penggunaan fasilitas Mandi, Cuci, Kakus (MCK) bersama, yang dapat meningkatkan peluang pajanan bakteri penyebab kusta. Hubungan antar komponen tersebut digambarkan melalui diagram kausalitas yang menunjukkan alur dari kondisi lingkungan hingga peningkatan prevalensi penyakit kusta.

2.2 Ukuran Asosiasi

2.2.1 Incidence Rate (IR)

Incidence rate merupakan ukuran frekuensi kejadian penyakit yang menggambarkan kecepatan munculnya kasus baru dalam suatu populasi berisiko selama periode waktu tertentu. Ukuran ini dihitung sebagai perbandingan antara jumlah kasus baru dengan total person-time yang diobservasi.

\[ IR_i = \frac{K_i}{PT_i} \]

2.2.2 Release From Treatment (RFT)

Release From Treatment (RFT) merupakan indikator program pengendalian penyakit yang menunjukkan persentase pasien yang telah menyelesaikan pengobatan sesuai standar nasional maupun WHO (World Health Organization).

\[ RFT_{MB,i} = \frac{SB_{MB,i}}{JK_{MB,i}} \times 100\% \]

2.3 Desain Studi

Penelitian ini menggunakan desain cross-sectional, dengan unit analisis berupa seluruh kabupaten/kota di Provinsi Jawa Timur yang berjumlah 38 wilayah. Pada desain ini, pengukuran variabel paparan dan outcome dilakukan pada waktu yang sama menggunakan data agregat wilayah. Outcome utama dalam penelitian ini adalah prevalensi penyakit kusta, yang pada beberapa analisis juga dinyatakan dalam bentuk insidensi per 1.000 penduduk.

Variabel paparan yang diteliti meliputi persentase rumah tangga yang menggunakan MCK bersama serta persentase rumah tangga yang belum memiliki akses terhadap sanitasi layak, yang merepresentasikan kondisi lingkungan dan sanitasi masyarakat. Populasi target dalam penelitian ini adalah seluruh penduduk yang berada di kabupaten/kota Provinsi Jawa Timur pada tahun 2020. Metode pengambilan sampel yang digunakan adalah total sampling, yaitu seluruh kabupaten/kota yang tersedia dalam dataset dimasukkan sebagai unit observasi, sehingga hasil analisis diharapkan dapat memberikan gambaran menyeluruh mengenai kondisi epidemiologi penyakit kusta di Provinsi Jawa Timur.

3 Metodologi Penelitian

3.1 Sumber Data

Penelitian ini menggunakan data sekunder yang bersifat cross-section yang digunakan untuk merepresentasikan kondisi kabupaten dan kota di Provinsi Jawa Timur pada satu tahun pengamatan.

Sumber data diperoleh dari publikasi resmi Dinas Kesehatan Provinsi Jawa Timur dengan judul “Profil Kesehatan Dinas Kesehatan Provinsi Jawa Timur 2021” yang dapat diakses melalui tautan https://dinkes.jatimprov.go.id/userfile/dokumen/PROFIL%20KESEHATAN%202020.pdf

3.2 Variabel Penelitian

Variabel yang digunakan meliputi variabel prediktor dan respon. Secara rinci, terdapat 3 variabel yang terdiri dari 2 variabel prediktor dan 1 variabel respon antara lain :

1. Prevalensi kasus penyakit kusta di Provinsi Jawa Timur (\(Y\))
2. Persentase rumah tangga yang menggunakan fasilitas Mandi, Cuci, Kakus Bersama (MCK Bersama) di Provinsi Jawa Timur (\(X_1\))
3. Persentase rumah tangga yang belum memiliki akses terhadap sanitasi layak (\(X_2\))

3.3 Metode Analisis

Metode analisis yang digunakan dalam penelitian ini terdiri atas beberapa tahapan yang saling berkaitan, dimulai dari analisis deskriptif hingga standardisasi angka penyakit kusta.

3.3.1 Analisis Deskriptif

Analisis deskriptif dilakukan untuk menggambarkan karakteristik data prevalensi penyakit kusta dan faktor lingkungan pada tingkat kabupaten/kota di Provinsi Jawa Timur. Statistik deskriptif yang dihitung meliputi nilai rata-rata, simpangan baku, nilai minimum, dan nilai maksimum untuk variabel prevalensi kusta, persentase rumah tangga yang menggunakan MCK bersama, serta persentase rumah tangga yang belum memiliki akses sanitasi layak. Selain itu, dilakukan pengurutan wilayah berdasarkan prevalensi kusta tertinggi dan terendah untuk mengidentifikasi daerah dengan beban penyakit yang paling menonjol. Distribusi prevalensi kusta juga divisualisasikan menggunakan histogram.

3.3.2 Perhitungan Ukuran Epidemiologi

Ukuran epidemiologi yang dihitung meliputi Incidence Rate (IR) dan Release From Treatment (RFT). Incidence rate dihitung sebagai perbandingan antara jumlah kasus baru dengan jumlah penduduk (dalam ribu jiwa) pada masing-masing kabupaten/kota, sebagai pendekatan person-time. Nilai IR digunakan untuk menggambarkan laju kejadian kasus baru penyakit kusta.

Selain itu, dihitung indikator Release From Treatment (RFT) untuk kusta tipe paucibacillary (PB) dan multibacillary (MB), yang dinyatakan sebagai persentase pasien yang menyelesaikan pengobatan dibandingkan dengan jumlah total kasus pada masing-masing tipe. Distribusi IR serta RFT divisualisasikan menggunakan histogram untuk melihat sebaran nilai antar wilayah.

3.3.3 Analisis Inferensial

Analisis inferensial dilakukan untuk mengkaji hubungan antara faktor lingkungan dan jumlah kasus baru kusta. Model regresi Poisson digunakan sebagai model awal karena variabel respon berupa data cacah (jumlah kasus baru). Untuk mengevaluasi adanya overdispersion, dihitung rasio dispersi berdasarkan residual Pearson.

Apabila terdeteksi overdispersion, analisis dilanjutkan menggunakan regresi binomial negatif. Hasil model disajikan dalam bentuk Risk Ratio (RR) beserta interval kepercayaan, yang diperoleh melalui eksponensiasi koefisien regresi. Variabel paparan yang dianalisis meliputi persentase rumah tangga yang menggunakan MCK bersama dan persentase rumah tangga tanpa akses sanitasi layak, dengan jumlah penduduk digunakan sebagai offset dalam model.

3.3.4 Pemetaan Penyakit

Pemetaan penyakit dilakukan untuk menggambarkan distribusi spasial prevalensi kusta di Provinsi Jawa Timur. Data prevalensi kusta digabungkan dengan peta batas administrasi kabupaten/kota. Peta tematik disusun dalam dua bentuk, yaitu peta kontinu yang menampilkan gradasi nilai prevalensi dan peta diskrit yang mengelompokkan wilayah ke dalam kategori tingkat prevalensi berdasarkan kuartil. Pemetaan ini bertujuan untuk mengidentifikasi pola sebaran dan wilayah dengan beban penyakit yang relatif tinggi.

3.3.5 Analisis Autokorelasi Spasial

Untuk menilai adanya ketergantungan spasial antar wilayah, dilakukan analisis Global Moran’s I menggunakan matriks pembobot spasial berbasis contiguity. Analisis ini digunakan untuk menguji apakah distribusi prevalensi kusta bersifat acak atau menunjukkan pengelompokan spasial. Selanjutnya, dilakukan analisis Local Moran’s I (LISA) untuk mengidentifikasi klaster lokal, yang diklasifikasikan ke dalam kelompok High-High, Low-Low, High-Low, dan Low-High. Hasil analisis LISA divisualisasikan dalam bentuk peta klaster untuk menunjukkan pola spasial lokal yang signifikan.

3.3.6 Standardisasi Angka Kejadian Penyakit

Untuk membandingkan tingkat kejadian kusta antar wilayah dengan struktur penduduk yang berbeda, dilakukan perhitungan Directly Standardized Rate (DSR) berdasarkan jenis kelamin. DSR dihitung dengan menggunakan populasi standar yang diperoleh dari total populasi penelitian. Selain itu, dihitung pula crude rate sebagai pembanding.

3.4 Alur Kerja Penelitian

Tahap-tahap yang dilakukan dalam analisis data pada penelitian ini adalah sebagai berikut:

4 Hasil dan Pembahasan

4.1 Analisis Deskriptif

# 1. PERSIAPAN LIBRARY DAN DATA #
options(tibble.width = Inf)
library(readxl)
library(dplyr)
library(janitor)
library(ggplot2)
library(sp)
library(sf)
library(MASS)
library(spdep)
library(mapview)
library(spData)
library(openxlsx)
library(spatialreg)
library(Matrix)
library(spgwr)
library(lmtest)
library(car)
library(gstat)
library(dplyr)
library(FNN)
library(stars)
library(ggrepel)
library(viridis)
library(raster)
library(tidyr)

# 2. ARAHKAN KE FOLDER KERJA
setwd("C:\\Users\\hp\\Documents\\New folder\\SMT 5\\EPIDEMIOLOGI") 
df <- read_excel("Data_Epidemiologi_Kusta.xlsx") %>%
  clean_names()

# 2. STATISTIK DESKRIPTIF #
df %>%
  summarise(
    n = n(),
    prev_mean = mean(prevalensi_kusta, na.rm = TRUE),
    prev_sd   = sd(prevalensi_kusta, na.rm = TRUE),
    prev_min  = min(prevalensi_kusta, na.rm = TRUE),
    prev_max  = max(prevalensi_kusta, na.rm = TRUE),
    
    mck_mean  = mean(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
    mck_sd    = sd(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
    mck_min   = min(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
    mck_max   = max(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
    
    san_mean  = mean(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
    san_sd    = sd(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
    san_min   = min(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
    san_max   = max(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE)
  )

## # A tibble: 1 × 13
##       n prev_mean prev_sd prev_min prev_max mck_mean mck_sd mck_min mck_max
##   <int>     <dbl>   <dbl>    <dbl>    <dbl>    <dbl>  <dbl>   <dbl>   <dbl>
## 1    38     0.443   0.519        0     2.06     17.9   13.1    1.29    55.9
##   san_mean san_sd san_min san_max
##      <dbl>  <dbl>   <dbl>   <dbl>
## 1     116.   71.7    8.09    266.

# Urutan prevalensi tertinggi dan terendah
top10 <- df %>%
  dplyr::arrange(desc(prevalensi_kusta)) %>%
  dplyr::select(kabupaten_kota, prevalensi_kusta) %>%
  slice(1:10) %>%
  rename(
    kabupaten_top = kabupaten_kota,
    prevalensi_top = prevalensi_kusta
  )
bottom10 <- df %>%
  dplyr::arrange(desc(prevalensi_kusta)) %>%
  dplyr::select(kabupaten_kota, prevalensi_kusta) %>%
  slice(29:38) %>%
  rename(
    kabupaten_bottom = kabupaten_kota,
    prevalensi_bottom = prevalensi_kusta
  )
tabel_banding <- bind_cols(top10, bottom10) ; tabel_banding

## # A tibble: 10 × 4
##    kabupaten_top prevalensi_top kabupaten_bottom prevalensi_bottom
##    <chr>                  <dbl> <chr>                        <dbl>
##  1 Sumenep                2.06  Kota Malang                  0.107
##  2 Pamekasan              1.91  Kota Kediri                  0.105
##  3 Sampang                1.82  Malang                       0.102
##  4 Bangkalan              1.11  Tulungagung                  0.101
##  5 Situbondo              1.09  Madiun                       0.094
##  6 Probolinggo            0.711 Blitar                       0.082
##  7 Tuban                  0.676 Pacitan                      0.017
##  8 Pasuruan               0.585 Kota Blitar                  0    
##  9 Lumajang               0.581 Kota Mojokerto               0    
## 10 Lamongan               0.573 Batu                         0

while (!is.null(dev.list())) dev.off()
graphics.off()

# Visualisasi Histogram
ggplot(df, aes(x = prevalensi_kusta)) +
  geom_histogram(
    bins = 10,
    fill = "pink",
    color = "white",
    alpha = 0.85
  ) +
  labs(
    title = "Distribusi Prevalensi Kusta",
    subtitle = "Per 1.000 penduduk, Provinsi Jawa Timur",
    x = "Prevalensi Kusta",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(size = 11),
    axis.title = element_text(face = "bold")
  )

Interpretasi :

Hasil statistik deskriptif menunjukkan bahwa prevalensi kusta antar kabupaten/kota memiliki variasi yang cukup besar, dengan nilai minimum sebesar 0 dan maksimum sebesar 2,06 per 1.000 penduduk. Nilai rata-rata prevalensi yang relatif kecil namun dengan simpangan baku yang cukup besar mengindikasikan adanya ketimpangan distribusi penyakit kusta antar wilayah.

Pengurutan wilayah berdasarkan prevalensi menunjukkan bahwa Kabupaten Sumenep, Pamekasan, dan Sampang merupakan wilayah dengan prevalensi tertinggi, sedangkan beberapa kota seperti Kota Blitar, Kota Mojokerto, dan Kota Batu memiliki prevalensi nol. Histogram prevalensi memperlihatkan distribusi yang tidak simetris dan cenderung right-skewed, yang menunjukkan bahwa sebagian besar kabupaten/kota memiliki prevalensi rendah, sementara hanya sedikit wilayah yang memiliki prevalensi sangat tinggi. Pola ini mengindikasikan bahwa beban penyakit kusta terkonsentrasi pada wilayah tertentu.

4.2 Ukuran Epidemiologi

# 3. UKURAN STATISTIK #
data <- df %>%
  dplyr::select(
    kabupaten_kota,
    kasus_baru,
    populasi_ribu,
    jumlah_kasus_pb,
    selesai_berobat_pb_8,
    jumlah_kasus_mb,
    selesai_berobat_pb_10
  )
names(data) <- c("kabupaten_kota", "kasus_baru", "populasi_ribu", "jumlah_kasus_pb", "selesai_berobat_pb", "jumlah_kasus_mb", "selesai_berobat_mb")

# Incidence Rate
ir <- data %>%
  mutate(
    person_time = as.numeric(populasi_ribu) * 1,
    incidence_rate = kasus_baru / person_time
    ) %>%
  dplyr::select(kabupaten_kota, kasus_baru, populasi_ribu, incidence_rate) ; ir

## # A tibble: 38 × 4
##    kabupaten_kota kasus_baru populasi_ribu incidence_rate
##    <chr>               <dbl>         <dbl>          <dbl>
##  1 Pacitan                 1          586.        0.00171
##  2 Ponorogo               20          949.        0.0211 
##  3 Trenggalek             10          731.        0.0137 
##  4 Tulungagung            11         1090.        0.0101 
##  5 Blitar                 10         1224.        0.00817
##  6 Kediri                 20         1635.        0.0122 
##  7 Malang                 27         2654.        0.0102 
##  8 Lumajang               65         1119.        0.0581 
##  9 Jember                134         2537.        0.0528 
## 10 Banyuwangi             32         1708.        0.0187 
## # ℹ 28 more rows

ggplot(ir, aes(x = incidence_rate)) +
  geom_histogram(
    bins = 10,
    fill = "pink",
    color = "white",
    alpha = 0.85
  ) +
  labs(
    title = "Incidence Rate Penyakit Kusta",
    subtitle = "di Provinsi Jawa Timur",
    x = "Incidence Rate Kusta",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(size = 11),
    axis.title = element_text(face = "bold")
  )

# Release From Treatment
rft <- data %>%
  mutate(
    rft_pb = case_when(
      jumlah_kasus_pb == 0 & selesai_berobat_pb == 0 ~ 100,
      jumlah_kasus_pb > 0 ~ (selesai_berobat_pb / jumlah_kasus_pb) * 100,
      TRUE ~ 0
    ),
    rft_mb = case_when(
      jumlah_kasus_mb == 0 & selesai_berobat_mb == 0 ~ 100,
      jumlah_kasus_mb > 0 ~ (selesai_berobat_mb / jumlah_kasus_mb) * 100,
      TRUE ~ 0
    )
  ) %>%
  dplyr::select(kabupaten_kota, rft_pb, rft_mb) ; rft

## # A tibble: 38 × 3
##    kabupaten_kota rft_pb rft_mb
##    <chr>           <dbl>  <dbl>
##  1 Pacitan         100    100  
##  2 Ponorogo        100    100  
##  3 Trenggalek       80     85.7
##  4 Tulungagung     100    100  
##  5 Blitar          100     85  
##  6 Kediri          100    100  
##  7 Malang          100     78.3
##  8 Lumajang        100     89.6
##  9 Jember           96.2   93.0
## 10 Banyuwangi      100     97.2
## # ℹ 28 more rows

ggplot(rft, aes(x = rft_pb)) +
  geom_histogram(
    bins = 10,
    fill = "pink",
    color = "white",
    alpha = 0.85
  ) +
  labs(
    title = " Selesai Berobat Penyakit Kusta",
    subtitle = "Tipe Paucibacillary di Provinsi Jawa Timur",
    x = "RFT Kusta PB",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(size = 11),
    axis.title = element_text(face = "bold")
  )

ggplot(rft, aes(x = rft_mb)) +
  geom_histogram(
    bins = 10,
    fill = "pink",
    color = "white",
    alpha = 0.85
  ) +
  labs(
    title = "Selesai Berobat Penyakit Kusta",
    subtitle = "Tipe Multibacillary di Provinsi Jawa Timur",
    x = "RFT Kusta MB",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(size = 11),
    axis.title = element_text(face = "bold")
  )

Interpretasi :

1. Incidence Rate (IR)

   Hasil perhitungan incidence rate menunjukkan adanya perbedaan laju kejadian kasus baru kusta antar kabupaten/kota. Histogram incidence rate memperlihatkan bahwa sebagian besar wilayah memiliki nilai IR yang rendah, namun terdapat beberapa kabupaten dengan nilai IR yang jauh lebih tinggi dibandingkan wilayah lain.

   Nilai incidence rate yang tinggi mencerminkan masih berlangsungnya transmisi penyakit di masyarakat, sehingga wilayah tersebut memerlukan perhatian khusus dalam hal deteksi dini, penguatan surveilans, dan intervensi kesehatan masyarakat.

2. Release From Treatment (TFT)

   Hasil analisis RFT menunjukkan bahwa sebagian besar kabupaten/kota memiliki nilai RFT yang tinggi, baik untuk tipe paucibacillary (PB) maupun multibacillary (MB). Histogram RFT PB menunjukkan dominasi nilai mendekati 100%, yang mengindikasikan keberhasilan pengobatan pada sebagian besar wilayah. Namun demikian, pada RFT MB terlihat variasi yang lebih besar, dengan beberapa wilayah memiliki nilai RFT yang lebih rendah.

   Perbedaan ini menunjukkan bahwa pengobatan kusta tipe MB yang memerlukan durasi lebih lama dan kepatuhan pengobatan yang lebih tinggi masih menghadapi tantangan di beberapa daerah.

4.3 Analisis Inferensial

# 4. ANALISIS INFERENSIAL
# Regresi Poisson
m_pois <- glm(
  df$kasus_baru ~ df$persentase_rumah_tangga_yang_menggunakan_mck_bersama + df$persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak + offset(log(df$populasi_ribu)),
  family = poisson(),
  data = df
)
disp <- sum(residuals(m_pois, type="pearson")^2) / df.residual(m_pois) ; disp

## [1] 38.09261

# Negatif Binomial
m_nb <- glm.nb(
  df$kasus_baru ~ df$persentase_rumah_tangga_yang_menggunakan_mck_bersama + df$persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak + offset(log(df$populasi_ribu)),
  data = df
)
exp(cbind(RR = coef(m_nb), confint(m_nb)))

##                                                                                       RR
## (Intercept)                                                                  0.009707724
## df$persentase_rumah_tangga_yang_menggunakan_mck_bersama                      1.030329630
## df$persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak 1.006393070
##                                                                                    2.5 %
## (Intercept)                                                                  0.005329685
## df$persentase_rumah_tangga_yang_menggunakan_mck_bersama                      1.008582915
## df$persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak 1.002198653
##                                                                                  97.5 %
## (Intercept)                                                                  0.01798434
## df$persentase_rumah_tangga_yang_menggunakan_mck_bersama                      1.05448130
## df$persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak 1.01070093

Interpretasi :

Analisis inferensial dilakukan menggunakan regresi Poisson untuk memodelkan hubungan antara faktor sanitasi dan jumlah kasus baru kusta. Hasil uji dispersi menunjukkan adanya overdispersion, sehingga model regresi binomial negatif digunakan sebagai model alternatif yang lebih sesuai.

Hasil regresi binomial negatif menunjukkan bahwa persentase rumah tangga yang menggunakan MCK bersama berhubungan positif dengan jumlah kasus baru kusta. Hal ini mengindikasikan bahwa peningkatan penggunaan MCK bersama berpotensi meningkatkan risiko penularan kusta, yang sejalan dengan teori epidemiologi mengenai peran sanitasi dan higiene dalam penularan penyakit menular. Sementara itu, variabel persentase rumah tangga tanpa akses sanitasi layak menunjukkan hubungan yang lebih lemah, yang mengindikasikan kemungkinan adanya faktor lain yang turut berperan.

4.4 Disease Mapping

# 5. DISEASE MAPPING #
setwd("C:\\Users\\hp\\Documents\\New folder\\SMT 5\\EPIDEMIOLOGI") 
kusta_data <- data.frame(df$kabupaten_kota, df$prevalensi_kusta)
colnames(kusta_data) <- c("wilayah", "kusta") 

# Peta administratif
Indo_Kab <- readRDS("gadm36_IDN_2_sp.rds")
Jatim <- Indo_Kab[Indo_Kab$NAME_1 == "Jawa Timur", ]
plot(Jatim)

Jatim_sf <- st_as_sf(Jatim)
Jatim_sf <- Jatim_sf %>%
  left_join(kusta_data, by = c("NAME_2" = "wilayah"))

cat("\nHasil gabung (Kolom Wilayah dan Kusta):\n")

## 
## Hasil gabung (Kolom Wilayah dan Kusta):

Jatim_sf %>%
  sf::st_drop_geometry() %>%
  dplyr::select(NAME_2, kusta) %>%
  print()

##              NAME_2 kusta
## 1         Bangkalan 1.113
## 2        Banyuwangi 0.187
## 3              Batu 0.000
## 4            Blitar 0.082
## 5        Bojonegoro 0.430
## 6         Bondowoso 0.386
## 7            Gresik 0.458
## 8            Jember 0.528
## 9           Jombang 0.402
## 10           Kediri 0.122
## 11      Kota Blitar 0.000
## 12      Kota Kediri 0.105
## 13      Kota Madiun 0.205
## 14      Kota Malang 0.107
## 15   Kota Mojokerto 0.000
## 16    Kota Pasuruan 0.481
## 17 Kota Probolinggo 0.292
## 18         Lamongan 0.573
## 19         Lumajang 0.581
## 20           Madiun 0.094
## 21          Magetan 0.343
## 22           Malang 0.102
## 23        Mojokerto 0.161
## 24          Nganjuk 0.181
## 25            Ngawi 0.241
## 26          Pacitan 0.017
## 27        Pamekasan 1.906
## 28         Pasuruan 0.585
## 29         Ponorogo 0.211
## 30      Probolinggo 0.711
## 31          Sampang 1.825
## 32         Sidoarjo 0.158
## 33        Situbondo 1.093
## 34          Sumenep 2.063
## 35         Surabaya 0.164
## 36       Trenggalek 0.137
## 37            Tuban 0.676
## 38      Tulungagung 0.101

cat("\nJumlah wilayah tanpa data Kusta: ", sum(is.na(Jatim_sf$kusta)), "\n")

## 
## Jumlah wilayah tanpa data Kusta:  0

# Pemetaan kasus kontinu
ggplot() +
  geom_sf(data=Jatim_sf, aes(fill = kusta), color=NA) +
  theme_bw() +
  scale_fill_gradient(low = "pink", high = "black") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.position = "right") +
  labs(title = "Peta Kasus Kusta Jawa Timur",
       fill = "Jumlah Kasus")

# Pemetaan kasus diskrit
quantile(Jatim_sf$kusta, probs = c(0, 0.25, 0.5, 0.75, 1))

##      0%     25%     50%     75%    100% 
## 0.00000 0.11075 0.22600 0.56175 2.06300

breaks <- c(0, 0.11075, 0.22600, 0.56175, 2.06300)
labels <- c("Very Low", "Low", "High", "Very High")

Jatim_sf$kusta_Discrete <- cut(Jatim_sf$kusta, 
                               breaks = breaks, 
                               labels = labels, 
                               right = TRUE)

ggplot() +
  geom_sf(data=Jatim_sf, aes(fill = kusta_Discrete), color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "pink", 
                               "Low" = "magenta",
                               "High" = "purple",
                               "Very High" = "black")) +
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.position = "right") +
  labs(title = "Kategori Kasus Kusta Jawa Timur",
       fill = "Kategori")

# Autokorelasi spasial (Moran's I)
row.names(Jatim) <- as.character(1:38)
W <- poly2nb(Jatim, row.names(Jatim), queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE)  
lwW <- nb2listw(W)

CoordK <- coordinates(Jatim)
plot(Jatim, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Jatim), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)

Global_Moran <- moran.test(Jatim_sf$kusta, lwW)
Global_Moran

## 
##  Moran I test under randomisation
## 
## data:  Jatim_sf$kusta  
## weights: lwW    
## 
## Moran I statistic standard deviate = 7.0385, p-value = 9.713e-13
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.84560039       -0.02702703        0.01537067

# Local Moran's I (LISA)
Local_Moran <- localmoran(Jatim_sf$kusta, lwW)

mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Jatim_sf$mean_values <- mean_values_char

ggplot() +
  geom_sf(data = Jatim_sf, aes(fill = mean_values)) +
  scale_fill_manual(values = c("Low-Low" = "pink", 
                               "High-Low" = "magenta", 
                               "Low-High" = "purple", 
                               "High-High" = "black")) +
  labs(title = "Local Moran's I Jawa Timur",
       fill = "Cluster Type") +
  theme_minimal()

x <- as.numeric(scale(Jatim_sf$kusta))
lagx <- spdep::lag.listw(lwW, x)

lisa <- spdep::localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
  x >= 0 & lagx >= 0 ~ "High-High",
  x <  0 & lagx <  0 ~ "Low-Low",
  x >= 0 & lagx <  0 ~ "High-Low (Outlier)",
  x <  0 & lagx >= 0 ~ "Low-High (Outlier)"
)
Jatim_LISA <- dplyr::bind_cols(Jatim_sf, lisa_df) |>
  dplyr::mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))

ggplot(Jatim_LISA) +
  geom_sf(aes(fill = quad), color="white", size=0.2) +
  scale_fill_manual(values=c(
    "High-High"="pink","Low-Low"="magenta",
    "High-Low (Outlier)"="purple","Low-High (Outlier)"="black",
    "Not significant"="grey85"
  )) +
  labs(title="Local Moran's I (LISA)", fill="Kategori") +
  theme_minimal()

Interpretasi :

Pemetaan penyakit menunjukkan bahwa prevalensi kusta di Jawa Timur tidak terdistribusi secara merata. Peta kontinu memperlihatkan gradasi warna yang lebih gelap pada wilayah dengan prevalensi tinggi, khususnya di wilayah Madura dan sebagian pesisir utara Jawa Timur. Sementara itu, peta diskrit menunjukkan adanya pengelompokan wilayah dengan kategori prevalensi tinggi dan sangat tinggi yang berdekatan secara geografis.

Hasil uji Global Moran’s I menunjukkan nilai yang signifikan secara statistik, yang menandakan adanya autokorelasi spasial positif pada prevalensi kusta. Artinya, kabupaten/kota dengan prevalensi kusta tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki prevalensi tinggi, demikian pula sebaliknya. Analisis Local Moran’s I (LISA) mengidentifikasi klaster High-High dan Low-Low yang signifikan. Klaster High-High menunjukkan wilayah prioritas yang memerlukan intervensi khusus karena tingginya prevalensi kusta yang diperkuat oleh pengaruh wilayah sekitarnya.

4.5 Standardized Morbidity/Mortality Ratio (SMR)

## Direct Standardization
dsr_direct_sex <- function(data,
                           area_col = "kabupaten_kota",
                           pop_m_col = "pop_laki_laki",
                           pop_f_col = "pop_perempuan",
                           case_m_col = "kasus_baru_laki_laki",
                           case_f_col = "kasus_baru_perempuan",
                           N_std = NULL,
                           scale = 1000) {
  
  area <- data[[area_col]]
  popM <- as.numeric(data[[pop_m_col]])
  popF <- as.numeric(data[[pop_f_col]])
  caseM <- as.numeric(data[[case_m_col]])
  caseF <- as.numeric(data[[case_f_col]])
  
  if (is.null(N_std)) {
    N_std <- c(sum(popM, na.rm = TRUE), sum(popF, na.rm = TRUE))
  }
  w <- N_std / sum(N_std)
  
  rM <- ifelse(popM > 0, caseM / popM, NA_real_)
  rF <- ifelse(popF > 0, caseF / popF, NA_real_)
  
  DSR <- (w[1] * rM + w[2] * rF)
  DSR_scaled <- DSR * scale
  
  var_DSR <- (w[1]^2 * rM / popM) + (w[2]^2 * rF / popF)
  se <- sqrt(var_DSR)
  
  lcl <- (DSR - 1.96 * se) * scale
  ucl <- (DSR + 1.96 * se) * scale
  
  out <- tibble::tibble(
    !!area_col := area,
    pop_m = popM, pop_f = popF,
    cases_m = caseM, cases_f = caseF,
    r_m = rM * scale, r_f = rF * scale,
    DSR = DSR_scaled,
    DSR_LCL = lcl,
    DSR_UCL = ucl
  )
  out
}

dsr_jatim <- dsr_direct_sex(
  data = df,
  area_col = "kabupaten_kota",
  pop_m_col = "populasi_laki_laki",
  pop_f_col = "populasi_perempuan",
  case_m_col = "kasus_baru_laki_laki",
  case_f_col = "kasus_baru_perempuan",
  N_std = NULL,
  scale = 1000
)

print(n=38, dsr_jatim %>% arrange(desc(DSR)) %>% dplyr::select(kabupaten_kota, DSR, DSR_LCL, DSR_UCL))

## # A tibble: 38 × 4
##    kabupaten_kota       DSR   DSR_LCL DSR_UCL
##    <chr>              <dbl>     <dbl>   <dbl>
##  1 Sumenep          0.207    0.180    0.234  
##  2 Pamekasan        0.190    0.161    0.219  
##  3 Sampang          0.183    0.156    0.210  
##  4 Bangkalan        0.111    0.0913   0.132  
##  5 Situbondo        0.110    0.0851   0.135  
##  6 Probolinggo      0.0715   0.0560   0.0870 
##  7 Tuban            0.0676   0.0529   0.0823 
##  8 Pasuruan         0.0585   0.0467   0.0703 
##  9 Lumajang         0.0582   0.0440   0.0723 
## 10 Lamongan         0.0573   0.0445   0.0701 
## 11 Jember           0.0528   0.0439   0.0618 
## 12 Kota Pasuruan    0.0481   0.0183   0.0778 
## 13 Gresik           0.0457   0.0341   0.0573 
## 14 Bojonegoro       0.0429   0.0317   0.0542 
## 15 Jombang          0.0401   0.0293   0.0509 
## 16 Bondowoso        0.0388   0.0249   0.0527 
## 17 Magetan          0.0345   0.0204   0.0486 
## 18 Kota Probolinggo 0.0292   0.00756  0.0508 
## 19 Ngawi            0.0242   0.0139   0.0346 
## 20 Ponorogo         0.0211   0.0118   0.0303 
## 21 Kota Madiun      0.0207   0.000411 0.0410 
## 22 Banyuwangi       0.0187   0.0122   0.0252 
## 23 Nganjuk          0.0180   0.0101   0.0259 
## 24 Surabaya         0.0164   0.0117   0.0211 
## 25 Mojokerto        0.0161   0.00864  0.0235 
## 26 Sidoarjo         0.0158   0.0104   0.0211 
## 27 Trenggalek       0.0136   0.00517  0.0220 
## 28 Kediri           0.0121   0.00682  0.0175 
## 29 Kota Malang      0.0107   0.00371  0.0177 
## 30 Kota Kediri      0.0104  -0.00137  0.0222 
## 31 Malang           0.0101   0.00632  0.0140 
## 32 Tulungagung      0.0101   0.00413  0.0161 
## 33 Madiun           0.00941  0.00244  0.0164 
## 34 Blitar           0.00811  0.00308  0.0131 
## 35 Pacitan          0.00170 -0.00163  0.00503
## 36 Kota Blitar      0        0        0      
## 37 Kota Mojokerto   0        0        0      
## 38 Batu             0        0        0

## Crude rate
df_crude <- df %>%
  mutate(
    pop_total = populasi_laki_laki + populasi_perempuan,
    kasus_total = kasus_baru_laki_laki + kasus_baru_perempuan,
    crude_ir = (kasus_total / pop_total) * 1000
  )

df_crude %>%
  dplyr::select(kabupaten_kota, crude_ir) %>%
  dplyr::arrange(desc(crude_ir))

## # A tibble: 38 × 2
##    kabupaten_kota crude_ir
##    <chr>             <dbl>
##  1 Sumenep          0.206 
##  2 Pamekasan        0.189 
##  3 Sampang          0.183 
##  4 Bangkalan        0.111 
##  5 Situbondo        0.109 
##  6 Probolinggo      0.0711
##  7 Tuban            0.0676
##  8 Pasuruan         0.0585
##  9 Lumajang         0.0581
## 10 Lamongan         0.0573
## # ℹ 28 more rows

## Tabel perbandingan crude rate vs DSR
compare_df <- df_crude %>%
  left_join(
    dsr_jatim %>% dplyr::select(kabupaten_kota, DSR),
    by = "kabupaten_kota"
  )

compare_df %>%
  dplyr::select(kabupaten_kota, crude_ir, DSR) %>%
  dplyr::arrange(desc(crude_ir))

## # A tibble: 38 × 3
##    kabupaten_kota crude_ir    DSR
##    <chr>             <dbl>  <dbl>
##  1 Sumenep          0.206  0.207 
##  2 Pamekasan        0.189  0.190 
##  3 Sampang          0.183  0.183 
##  4 Bangkalan        0.111  0.111 
##  5 Situbondo        0.109  0.110 
##  6 Probolinggo      0.0711 0.0715
##  7 Tuban            0.0676 0.0676
##  8 Pasuruan         0.0585 0.0585
##  9 Lumajang         0.0581 0.0582
## 10 Lamongan         0.0573 0.0573
## # ℹ 28 more rows

hist_df <- compare_df %>%
  dplyr::select(kabupaten_kota, crude_ir, DSR) %>%
  tidyr::pivot_longer(
    cols = c(crude_ir, DSR),
    names_to = "jenis",
    values_to = "rate"
  ) %>%
  mutate(
    jenis = dplyr::recode(jenis,
                   crude_ir = "Crude Rate",
                   DSR = "DSR (Standardized)")
  )

ggplot(hist_df, aes(x = rate, fill = jenis)) +
  geom_histogram(
    bins = 10,
    alpha = 0.7,
    color = "white"
  ) +
  facet_wrap(~ jenis, scales = "free_x") +
  labs(
    title = "Perbandingan Distribusi Crude Rate \ndan DSR Kusta",
    subtitle = "Kabupaten/Kota di Jawa Timur",
    x = "Rate per 1.000 penduduk",
    y = "Jumlah Kabupaten/Kota"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold")
  )

Interpretasi :

Hasil perhitungan Directly Standardized Rate (DSR) berdasarkan jenis kelamin menunjukkan bahwa nilai DSR relatif tidak jauh berbeda dibandingkan dengan crude rate pada sebagian besar kabupaten/kota di Provinsi Jawa Timur. Perbedaan antara kedua ukuran tersebut cenderung kecil dan tidak mengubah pola distribusi maupun peringkat wilayah dengan kejadian kusta tinggi.

Temuan ini mengindikasikan bahwa struktur jenis kelamin bukan merupakan faktor utama yang memengaruhi variasi kejadian penyakit kusta antar wilayah. Dengan demikian, perbedaan beban penyakit kusta di Provinsi Jawa Timur lebih kemungkinan dipengaruhi oleh faktor lain, seperti kondisi lingkungan, sanitasi, kepadatan penduduk, serta keterkaitan spasial antar wilayah, dibandingkan oleh perbedaan komposisi jenis kelamin penduduk.

5 Kesimpulan dan Saran

5.1 Kesimpulan

Berdasarkan hasil analisis epidemiologi dan spasial penyakit kusta di Provinsi Jawa Timur, dapat disimpulkan beberapa hal sebagai berikut :

1. Distribusi prevalensi penyakit kusta antar kabupaten/kota menunjukkan variasi yang cukup besar dan tidak tersebar secara merata. Wilayah dengan prevalensi kusta tinggi cenderung terkonsentrasi pada daerah tertentu, sementara beberapa kabupaten/kota lainnya memiliki prevalensi yang sangat rendah bahkan nol. Hal ini menunjukkan bahwa beban penyakit kusta di Jawa Timur bersifat heterogen antar wilayah.
2. Hasil perhitungan ukuran epidemiologi menunjukkan bahwa incidence rate kusta masih bervariasi antar kabupaten/kota, mengindikasikan bahwa penularan penyakit kusta masih aktif terjadi di beberapa wilayah. Sementara itu, indikator Release From Treatment (RFT) menunjukkan bahwa sebagian besar wilayah telah mencapai tingkat keberhasilan pengobatan yang tinggi, meskipun pada tipe multibacillary masih terdapat variasi antar daerah.
3. Analisis inferensial menunjukkan bahwa faktor lingkungan, khususnya persentase rumah tangga yang menggunakan MCK bersama, berhubungan dengan peningkatan jumlah kasus baru kusta. Hasil ini memperkuat peran sanitasi dan higiene sebagai determinan penting dalam penularan penyakit kusta. Selain itu, hasil analisis autokorelasi spasial menunjukkan adanya ketergantungan spasial yang signifikan, di mana prevalensi kusta pada suatu wilayah dipengaruhi oleh wilayah-wilayah yang bertetangga.
4. Hasil standardisasi angka kejadian penyakit menunjukkan bahwa perbedaan antara crude rate dan Directly Standardized Rate (DSR) berdasarkan jenis kelamin relatif kecil. Temuan ini mengindikasikan bahwa struktur jenis kelamin bukan merupakan faktor utama yang menjelaskan variasi kejadian penyakit kusta antar wilayah. Dengan demikian, variasi spasial penyakit kusta di Jawa Timur lebih dipengaruhi oleh faktor lingkungan dan keterkaitan spasial antar wilayah dibandingkan faktor demografis jenis kelamin.
5. Bagi pemerintah daerah dan dinas kesehatan, hasil penelitian ini dapat digunakan sebagai dasar dalam perencanaan program pengendalian kusta yang lebih terarah dan berbasis wilayah, khususnya pada daerah-daerah dengan prevalensi dan incidence rate yang tinggi. Upaya peningkatan sanitasi lingkungan, terutama pengurangan penggunaan MCK bersama dan peningkatan akses terhadap sanitasi layak, perlu menjadi prioritas dalam strategi pencegahan dan pengendalian penyakit kusta.

5.2 Saran

Berdasarkan kesimpulan penelitian, beberapa saran yang dapat diberikan adalah sebagai berikut :

1. Penelitian selanjutnya disarankan untuk memasukkan variabel lain yang belum tercakup dalam penelitian ini, seperti kepadatan penduduk, tingkat kemiskinan, akses pelayanan kesehatan, dan faktor sosial ekonomi, serta menggunakan data longitudinal agar dinamika penularan kusta dapat dianalisis secara lebih mendalam.
2. Pendekatan analisis spasial disarankan untuk terus digunakan dalam kajian epidemiologi penyakit menular, karena mampu memberikan gambaran yang lebih komprehensif dibandingkan analisis non-spasial dalam mengidentifikasi wilayah prioritas intervensi.

6 Lampiran

Berikut adalah tautan dari analisis menggunakan Dashboard yang dapat diakses beserta dengan lampiran syntax :

https://fazilaazraanggina.shinyapps.io/EPIDEMIOLOGI/

---

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

# ============================================================
# LIBRARY
# ============================================================
library(shiny)
library(shinydashboard)
library(shinycssloaders)

library(readxl)
library(janitor)

library(sf)
library(sp)
library(spdep)
library(ggplot2)
library(MASS)

library(DT)
library(leaflet)
library(plotly)
library(stringr)
library(openxlsx)

# IMPORTANT: dplyr & tidyr taruh bawah biar ga ketiban
library(dplyr)
library(tidyr)

# ============================================================
# HELPER
# ============================================================
norm_name <- function(x) {
  x %>%
    as.character() %>%
    stringr::str_trim() %>%
    stringr::str_to_upper() %>%
    stringr::str_replace_all("\\s+", " ") %>%
    stringr::str_replace_all("[^A-Z0-9\\s/\\-\\.]", "")
}

read_map_any <- function(path) {
  ext <- tolower(tools::file_ext(path))
  if (ext == "rds") return(sf::st_as_sf(readRDS(path)))
  sf::st_read(path, quiet = TRUE)
}

pal_bupk <- function(n = 256) {
  grDevices::colorRampPalette(c("#0B1F4B", "#7C3AED", "#EC4899"))(n)
}

pal_div_bupk <- function(n = 256) {
  grDevices::colorRampPalette(c("#0B1F4B", "#F9FAFB", "#EC4899"))(n)
}

# DSR function (sesuai syntax kamu)
dsr_direct_sex <- function(
    data,
    area_col = "kabupaten_kota",
    pop_m_col = "populasi_laki_laki",
    pop_f_col = "populasi_perempuan",
    case_m_col = "kasus_baru_laki_laki",
    case_f_col = "kasus_baru_perempuan",
    N_std = NULL,
    scale = 1000
) {
  area  <- data[[area_col]]
  popM  <- as.numeric(data[[pop_m_col]])
  popF  <- as.numeric(data[[pop_f_col]])
  caseM <- as.numeric(data[[case_m_col]])
  caseF <- as.numeric(data[[case_f_col]])
  
  if (is.null(N_std)) {
    N_std <- c(sum(popM, na.rm = TRUE), sum(popF, na.rm = TRUE))
  }
  w <- N_std / sum(N_std)
  
  rM <- ifelse(popM > 0, caseM / popM, NA_real_)
  rF <- ifelse(popF > 0, caseF / popF, NA_real_)
  
  DSR        <- (w[1] * rM + w[2] * rF)
  DSR_scaled <- DSR * scale
  
  var_DSR <- (w[1]^2 * rM / popM) + (w[2]^2 * rF / popF)
  se      <- sqrt(var_DSR)
  
  lcl <- (DSR - 1.96 * se) * scale
  ucl <- (DSR + 1.96 * se) * scale
  
  tibble::tibble(
    !!area_col := area,
    DSR = DSR_scaled,
    DSR_LCL = lcl,
    DSR_UCL = ucl
  )
}

# ============================================================
# UI (tetap mirip template kamu)
# ============================================================
ui <- dashboardPage(
  skin = "purple",
  dashboardHeader(title = span("UAS Epidemiologi", style = "font-weight:900;")),
  dashboardSidebar(
    width = 300,
    sidebarMenu(
      id = "tabs",
      menuItem("Overview", tabName = "overview", icon = icon("gauge-high")),
      menuItem("Upload", tabName = "upload", icon = icon("upload")),
      menuItem("Deskriptif", tabName = "desc", icon = icon("table")),
      menuItem("Ukuran Statistik", tabName = "measures", icon = icon("calculator")),
      menuItem("Regresi", tabName = "reg", icon = icon("chart-line")),
      menuItem("Peta Kusta", tabName = "map_kusta", icon = icon("map")),
      menuItem("Autokorelasi", tabName = "auto", icon = icon("project-diagram")),
      menuItem("Standardisasi", tabName = "dsr", icon = icon("scale-balanced")),
      menuItem("Export", tabName = "export", icon = icon("file-export")),
      menuItem("About", tabName = "about", icon = icon("circle-info"))
    )
  ),
  dashboardBody(
    tags$head(
      tags$style(
        HTML("
        body, .content-wrapper, .right-side { background:#070A1A !important; }
        .skin-purple .main-header .logo,
        .skin-purple .main-header .navbar { background:#0B0F2A !important; }
        .skin-purple .main-sidebar { background:#070A1A !important; }

        .sidebar a { color:#E8EAF6 !important; }
        .sidebar-menu>li>a:hover { background:#0E1334 !important; }
        .skin-purple .sidebar-menu>li.active>a { border-left-color:#EC4899 !important; }

        .box { border-radius:16px !important; border:1px solid #1B2460 !important; background:#F9FAFB !important; }
        .box-title, h1,h2,h3,h4,h5,label, .control-label { color:#111827 !important; font-weight:900 !important; }
        .box-body, .help-block, .shiny-text-output, pre { color:#111827 !important; }

        .small-box h3 { color:#FFFFFF !important; font-weight:900 !important; }
        .small-box p  { color:#FFFFFF !important; font-weight:800 !important; }
        .small-box { border-radius:16px !important; }

        .leaflet-container { border-radius:16px !important; }
        table.dataTable { color:#111827 !important; }
      ")
      )
    ),
    tabItems(
      # OVERVIEW
      tabItem(
        tabName = "overview",
        fluidRow(
          valueBoxOutput("vb_n", width = 3),
          valueBoxOutput("vb_prev_mean", width = 3),
          valueBoxOutput("vb_disp", width = 3),
          valueBoxOutput("vb_moran", width = 3)
        ),
        fluidRow(
          box(
            title = "Histogram Prevalensi Kusta",
            width = 12,
            status = "primary",
            sliderInput("bins_prev", "Bins", min = 5, max = 40, value = 10),
            plotlyOutput("plt_prev_hist", height = 380) %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Peta Prevalensi Kusta",
            width = 12,
            status = "info",
            radioButtons(
              "map_mode",
              "Mode peta:",
              choices = c("Kontinu" = "cont", "Diskrit" = "disc"),
              inline = TRUE
            ),
            leafletOutput("map_prev", height = 560) %>% withSpinner()
          )
        )
      ),
      
      # UPLOAD
      tabItem(
        tabName = "upload",
        fluidRow(
          box(
            title = "Upload Data",
            width = 12,
            status = "primary",
            fileInput("excel", "Upload Data (.xlsx)", accept = c(".xlsx")),
            fileInput("peta", "Upload Peta (RDS)", accept = c(".rds", ".shp", ".gpkg", ".geojson", ".zip")),
            actionButton("load", "Muat & Gabungkan", class = "btn btn-primary btn-lg"),
            tags$hr()
          )
        ),
        fluidRow(
          box(
            title = "Preview Excel",
            width = 12,
            status = "info",
            DTOutput("tbl_excel_raw") %>% withSpinner()
          )
        )
      ),
      
      # DESKRIPTIF
      tabItem(
        tabName = "desc",
        fluidRow(
          box(
            title = "Statistik Deskriptif",
            width = 12,
            status = "primary",
            DTOutput("tbl_desc") %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "10 Kabupaten/Kota dengan Prevalensi Tertinggi",
            width = 6,
            status = "warning",
            DTOutput("tbl_top10") %>% withSpinner()
          ),
          box(
            title = "10 Kabupaten/Kota dengan Prevalensi Terendah",
            width = 6,
            status = "success",
            DTOutput("tbl_bottom10") %>% withSpinner()
          )
        )
      ),
      
      # UKURAN STATISTIK
      tabItem(
        tabName = "measures",
        fluidRow(
          box(
            title = "Incidence Rate",
            width = 12,
            status = "primary",
            sliderInput("bins_ir", "Bins IR", min = 5, max = 40, value = 10),
            plotlyOutput("plt_ir_hist", height = 320) %>% withSpinner(),
            DTOutput("tbl_ir") %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Release From Treatment Tipe Paucibacillary",
            width = 6,
            status = "warning",
            sliderInput("bins_rft_pb", "Bins RFT PB", min = 5, max = 40, value = 10),
            plotlyOutput("plt_rft_pb", height = 300) %>% withSpinner()
          ),
          box(
            title = "Release From Treatment Tipe Multibacillary",
            width = 6,
            status = "warning",
            sliderInput("bins_rft_mb", "Bins RFT MB", min = 5, max = 40, value = 10),
            plotlyOutput("plt_rft_mb", height = 300) %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Tabel Perbandingan RFT (Paucibacillary & Multibacillary)",
            width = 12,
            status = "info",
            DTOutput("tbl_rft") %>% withSpinner()
          )
        )
      ),
      
      # REGRESI
      tabItem(
        tabName = "reg",
        fluidRow(
          box(
            title = "Analisis Regresi Poisson",
            width = 6,
            status = "primary",
            verbatimTextOutput("txt_pois") %>% withSpinner()
          ),
          box(
            title = "Analisis Negatif Binomial",
            width = 6,
            status = "primary",
            DTOutput("tbl_nb_rr") %>% withSpinner()
          )
        )
      ),
      
      # PETA
      tabItem(
        tabName = "map_kusta",
        fluidRow(
          box(
            title = "Peta Kasus Kusta Jawa Timur (Kontinu)",
            width = 6,
            status = "primary",
            plotOutput("gg_kusta_cont", height = 450) %>% withSpinner()
          ),
          box(
            title = "Peta Kasus Kusta Jawa Timur (Diskrit)",
            width = 6,
            status = "warning",
            plotOutput("gg_kusta_disc", height = 450) %>% withSpinner()
          )
        ),
        box(
          title = "Breaks Diskrit",
          width = 12,
          status = "info",
          numericInput("b0", "Break 0", value = 0),
          numericInput("b1", "Break 1", value = 0.11075),
          numericInput("b2", "Break 2", value = 0.22600),
          numericInput("b3", "Break 3", value = 0.56175),
          numericInput("b4", "Break 4", value = 2.06300)
        )
      ),
      
      # AUTOKORELASI
      tabItem(
        tabName = "auto",
        fluidRow(
          box(
            title = "Setting Bobot",
            width = 4,
            status = "info",
            selectInput(
              "queen_rook",
              "Rook/Queen",
              choices = c("Rook (queen=FALSE)" = "rook", "Queen (queen=TRUE)" = "queen"),
              selected = "rook"
            ),
            selectInput("w_style", "Style", choices = c("W", "B", "C", "U", "S"), selected = "W"),
            checkboxInput("zero_policy", "zero.policy=TRUE", value = TRUE),
            numericInput("snap", "snap (poly2nb)", value = 1e-6, min = 0, step = 1e-6)
          ),
          box(
            title = "Global Moran's I",
            width = 8,
            status = "primary",
            verbatimTextOutput("txt_moran") %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Local Moran (LISA)",
            width = 6,
            status = "info",
            DTOutput("tbl_lisa") %>% withSpinner()
          ),
          box(
            title = "LISA Cluster Map",
            width = 6,
            status = "warning",
            leafletOutput("map_lisa", height = 420) %>% withSpinner()
          )
        )
      ),
      
      # DSR
      tabItem(
        tabName = "dsr",
        fluidRow(
          box(
            title = "Direct Standardization",
            width = 12,
            status = "primary",
            DTOutput("tbl_dsr") %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Tabel Perbandingan (Crude Rate & DSR",
            width = 12,
            status = "info",
            DTOutput("tbl_compare") %>% withSpinner()
          )
        ),
        fluidRow(
          box(
            title = "Histogram Perbandingan (Crude Rate & DSR)",
            width = 12,
            status = "warning",
            sliderInput("bins_dsr", "Bins", min = 5, max = 40, value = 10),
            plotlyOutput("plt_hist_dsr", height = 360) %>% withSpinner()
          )
        )
      ),
      
      # EXPORT
      tabItem(
        tabName = "export",
        box(
          title = "Download",
          width = 12,
          status = "primary",
          downloadButton("dl_excel_raw", "Download Excel Raw (CSV)"),
          downloadButton("dl_clean", "Download Data Gabungan (CSV)"),
          downloadButton("dl_lisa", "Download LISA (CSV)"),
          downloadButton("dl_dsr", "Download DSR (CSV)")
        )
      ),
      
      # ABOUT
      tabItem(
        tabName = "about",
        box(
          title = "Tentang",
          width = 12,
          status = "primary",
          HTML("<div style='color:#111827'>
              <b>Dashboard Projek UAS Epidemiologi</b><br/>
              Oleh : Fazila Azra Anggina (140610230039).<br/>
            </div>")
        )
      )
    )
  )
)

# ============================================================
# SERVER
# ============================================================
server <- function(input, output, session) {
  
  # Excel RAW preview
  excel_raw <- reactive({
    req(input$excel)
    openxlsx::read.xlsx(input$excel$datapath)
  })
  
  output$tbl_excel_raw <- renderDT({
    req(excel_raw())
    DT::datatable(excel_raw(), options = list(pageLength = 12, scrollX = TRUE))
  })
  
  # Load & Join
  dat_pack <- eventReactive(input$load, {
    req(input$excel, input$peta)
    
    withProgress(message = "Memuat & memproses data...", value = 0, {
      
      incProgress(0.2, detail = "Baca Excel + clean_names()...")
      df <- readxl::read_excel(input$excel$datapath) %>% janitor::clean_names()
      
      need_cols <- c(
        "kabupaten_kota",
        "prevalensi_kusta",
        "kasus_baru",
        "populasi_ribu",
        "jumlah_kasus_pb",
        "selesai_berobat_pb_8",
        "jumlah_kasus_mb",
        "selesai_berobat_pb_10",
        "persentase_rumah_tangga_yang_menggunakan_mck_bersama",
        "persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak",
        "populasi_laki_laki",
        "populasi_perempuan",
        "kasus_baru_laki_laki",
        "kasus_baru_perempuan"
      )
      miss <- setdiff(need_cols, names(df))
      validate(need(length(miss) == 0, paste0("Kolom wajib hilang: ", paste(miss, collapse = ", "))))
      
      df <- df %>%
        dplyr::mutate(wilayah_key = norm_name(kabupaten_kota))
      
      incProgress(0.55, detail = "Baca peta...")
      Indo_Kab <- read_map_any(input$peta$datapath)
      nm <- names(Indo_Kab)
      validate(need(any(nm == "NAME_1"), "Peta tidak punya kolom NAME_1."))
      validate(need(any(nm == "NAME_2"), "Peta tidak punya kolom NAME_2."))
      
      incProgress(0.75, detail = "Filter Jawa Timur & join...")
      Jatim <- Indo_Kab[Indo_Kab$NAME_1 == "Jawa Timur", ]
      Jatim_sf <- sf::st_as_sf(Jatim) %>%
        dplyr::mutate(wilayah_key = norm_name(NAME_2)) %>%
        dplyr::left_join(
          df %>%
            dplyr::transmute(
              wilayah_key,
              kabupaten_kota,
              prevalensi_kusta = as.numeric(prevalensi_kusta),
              kasus_baru = as.numeric(kasus_baru),
              populasi_ribu = as.numeric(populasi_ribu),
              jumlah_kasus_pb = as.numeric(jumlah_kasus_pb),
              selesai_berobat_pb = as.numeric(selesai_berobat_pb_8),
              jumlah_kasus_mb = as.numeric(jumlah_kasus_mb),
              selesai_berobat_mb = as.numeric(selesai_berobat_pb_10),
              mck = as.numeric(persentase_rumah_tangga_yang_menggunakan_mck_bersama),
              san = as.numeric(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak),
              pop_m = as.numeric(populasi_laki_laki),
              pop_f = as.numeric(populasi_perempuan),
              case_m = as.numeric(kasus_baru_laki_laki),
              case_f = as.numeric(kasus_baru_perempuan)
            ),
          by = "wilayah_key"
        ) %>%
        dplyr::mutate(
          incidence_rate = kasus_baru / populasi_ribu,
          rft_pb = dplyr::case_when(
            jumlah_kasus_pb == 0 & selesai_berobat_pb == 0 ~ 100,
            jumlah_kasus_pb > 0 ~ (selesai_berobat_pb / jumlah_kasus_pb) * 100,
            TRUE ~ 0
          ),
          rft_mb = dplyr::case_when(
            jumlah_kasus_mb == 0 & selesai_berobat_mb == 0 ~ 100,
            jumlah_kasus_mb > 0 ~ (selesai_berobat_mb / jumlah_kasus_mb) * 100,
            TRUE ~ 0
          )
        )
      
      incProgress(1.0, detail = "Selesai")
      list(df = df, jatim_sf = Jatim_sf)
    })
  })
  
  df <- reactive({ req(dat_pack()); dat_pack()$df })
  jatim_sf_full <- reactive({ req(dat_pack()); dat_pack()$jatim_sf })
  
  jatim_sf_kusta <- reactive({
    req(jatim_sf_full())
    jatim_sf_full() %>% dplyr::filter(!is.na(prevalensi_kusta))
  })
  
  # =========================
  # DESKRIPTIF
  # =========================
  desc_tbl <- reactive({
    req(df())
    df() %>%
      dplyr::summarise(
        n = dplyr::n(),
        prev_mean = mean(prevalensi_kusta, na.rm = TRUE),
        prev_sd   = sd(prevalensi_kusta, na.rm = TRUE),
        prev_min  = min(prevalensi_kusta, na.rm = TRUE),
        prev_max  = max(prevalensi_kusta, na.rm = TRUE),
        
        mck_mean  = mean(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
        mck_sd    = sd(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
        mck_min   = min(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
        mck_max   = max(persentase_rumah_tangga_yang_menggunakan_mck_bersama, na.rm = TRUE),
        
        san_mean  = mean(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
        san_sd    = sd(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
        san_min   = min(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE),
        san_max   = max(persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak, na.rm = TRUE)
      )
  })
  
  top10_tbl <- reactive({
    req(df())
    df() %>%
      dplyr::arrange(dplyr::desc(prevalensi_kusta)) %>%
      dplyr::select(kabupaten_kota, prevalensi_kusta) %>%
      dplyr::slice(1:10)
  })
  
  bottom10_tbl <- reactive({
    req(df())
    df() %>%
      dplyr::arrange(prevalensi_kusta) %>%
      dplyr::select(kabupaten_kota, prevalensi_kusta) %>%
      dplyr::slice(1:10)
  })
  
  # =========================
  # IR (Incidence Rate)
  # =========================
  ir_tbl <- reactive({
    req(df())
    df() %>%
      dplyr::select(kabupaten_kota, kasus_baru, populasi_ribu) %>%
      dplyr::mutate(
        person_time = as.numeric(populasi_ribu) * 1,
        incidence_rate = kasus_baru / person_time
      ) %>%
      dplyr::select(kabupaten_kota, kasus_baru, populasi_ribu, incidence_rate)
  })
  
  # =========================
  # RFT
  # =========================
  rft_tbl <- reactive({
    req(df())
    tmp <- df() %>%
      dplyr::select(
        kabupaten_kota,
        jumlah_kasus_pb,
        selesai_berobat_pb_8,
        jumlah_kasus_mb,
        selesai_berobat_pb_10
      )
    names(tmp) <- c("kabupaten_kota", "jumlah_kasus_pb", "selesai_berobat_pb", "jumlah_kasus_mb", "selesai_berobat_mb")
    
    tmp %>%
      dplyr::mutate(
        rft_pb = dplyr::case_when(
          jumlah_kasus_pb == 0 & selesai_berobat_pb == 0 ~ 100,
          jumlah_kasus_pb > 0 ~ (selesai_berobat_pb / jumlah_kasus_pb) * 100,
          TRUE ~ 0
        ),
        rft_mb = dplyr::case_when(
          jumlah_kasus_mb == 0 & selesai_berobat_mb == 0 ~ 100,
          jumlah_kasus_mb > 0 ~ (selesai_berobat_mb / jumlah_kasus_mb) * 100,
          TRUE ~ 0
        )
      ) %>%
      dplyr::select(kabupaten_kota, rft_pb, rft_mb)
  })
  
  # =========================
  # REGRESI Poisson + Dispersion
  # =========================
  pois_pack <- reactive({
    req(df())
    d <- df()
    m_pois <- glm(
      kasus_baru ~ mck + san + offset(log(populasi_ribu)),
      family = poisson(),
      data = dplyr::mutate(
        d,
        mck = persentase_rumah_tangga_yang_menggunakan_mck_bersama,
        san = persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak
      )
    )
    disp <- sum(residuals(m_pois, type = "pearson")^2) / df.residual(m_pois)
    list(m = m_pois, disp = disp)
  })
  
  # NegBin RR + CI
  nb_rr_tbl <- reactive({
    req(df())
    d <- df() %>%
      dplyr::mutate(
        mck = persentase_rumah_tangga_yang_menggunakan_mck_bersama,
        san = persentase_rumah_tangga_yang_belum_memiliki_akses_terhadap_sanitasi_layak
      )
    
    m_nb <- MASS::glm.nb(
      kasus_baru ~ mck + san + offset(log(populasi_ribu)),
      data = d
    )
    
    rr <- exp(cbind(RR = coef(m_nb), confint(m_nb)))
    rr_df <- as.data.frame(rr)
    rr_df$term <- rownames(rr_df)
    
    rr_df %>% dplyr::select(term, RR, `2.5 %`, `97.5 %`)
  })
  
  # =========================
  # WEIGHTS + MORAN + LISA
  # =========================
  lw_kusta <- reactive({
    req(jatim_sf_kusta())
    sf0 <- jatim_sf_kusta()
    queen <- (input$queen_rook == "queen")
    nb <- spdep::poly2nb(sf::as_Spatial(sf0), queen = queen, snap = input$snap)
    spdep::nb2listw(nb, style = input$w_style, zero.policy = input$zero_policy)
  })
  
  moran_global <- reactive({
    req(jatim_sf_kusta(), lw_kusta())
    spdep::moran.test(
      jatim_sf_kusta()$prevalensi_kusta,
      lw_kusta(),
      zero.policy = input$zero_policy,
      na.action = na.exclude
    )
  })
  
  lisa_tbl <- reactive({
    req(jatim_sf_kusta(), lw_kusta())
    sf0 <- jatim_sf_kusta()
    lw <- lw_kusta()
    
    x <- as.numeric(scale(sf0$prevalensi_kusta))
    lagx <- spdep::lag.listw(lw, x, zero.policy = input$zero_policy)
    lisa <- spdep::localmoran(
      x,
      lw,
      alternative = "two.sided",
      zero.policy = input$zero_policy,
      na.action = na.exclude
    )
    
    lisa_df <- as.data.frame(lisa)
    names(lisa_df) <- c("Ii", "Ei", "Vi", "Zi", "Pi.two.sided")
    
    alpha <- 0.05
    quad <- dplyr::case_when(
      x >= 0 & lagx >= 0 ~ "High-High",
      x <  0 & lagx <  0 ~ "Low-Low",
      x >= 0 & lagx <  0 ~ "High-Low (Outlier)",
      x <  0 & lagx >= 0 ~ "Low-High (Outlier)"
    )
    
    out <- dplyr::bind_cols(
      sf0 %>% sf::st_drop_geometry() %>% dplyr::select(NAME_2, prevalensi_kusta),
      lisa_df
    ) %>%
      dplyr::mutate(quad = ifelse(`Pi.two.sided` <= alpha, quad, "Not significant"))
    
    out
  })
  
  # =========================
  # DSR + crude
  # =========================
  dsr_tbl <- reactive({
    req(df())
    dsr_direct_sex(
      data = df(),
      area_col = "kabupaten_kota",
      pop_m_col = "populasi_laki_laki",
      pop_f_col = "populasi_perempuan",
      case_m_col = "kasus_baru_laki_laki",
      case_f_col = "kasus_baru_perempuan",
      N_std = NULL,
      scale = 1000
    ) %>% dplyr::arrange(dplyr::desc(DSR))
  })
  
  compare_tbl <- reactive({
    req(df(), dsr_tbl())
    df_crude <- df() %>%
      dplyr::mutate(
        pop_total = populasi_laki_laki + populasi_perempuan,
        kasus_total = kasus_baru_laki_laki + kasus_baru_perempuan,
        crude_ir = (kasus_total / pop_total) * 1000
      ) %>%
      dplyr::select(kabupaten_kota, crude_ir)
    
    df_crude %>%
      dplyr::left_join(dsr_tbl() %>% dplyr::select(kabupaten_kota, DSR), by = "kabupaten_kota") %>%
      dplyr::arrange(dplyr::desc(crude_ir))
  })
  
  # ============================================================
  # OUTPUTS
  # ============================================================
  output$vb_n <- renderValueBox({
    req(df())
    valueBox(nrow(df()), "Jumlah Kab/Kota", color = "purple", icon = icon("list"))
  })
  
  output$vb_prev_mean <- renderValueBox({
    req(desc_tbl())
    valueBox(round(desc_tbl()$prev_mean, 4), "Rata-Rata Prevalensi", color = "fuchsia", icon = icon("virus"))
  })
  
  output$vb_disp <- renderValueBox({
    req(pois_pack())
    valueBox(round(pois_pack()$disp, 4), "Dispersion (Poisson)", color = "navy", icon = icon("wave-square"))
  })
  
  output$vb_moran <- renderValueBox({
    req(moran_global())
    val <- unname(moran_global()$estimate[["Moran I statistic"]])
    valueBox(round(val, 4), "Hasil Moran's I", color = "purple", icon = icon("project-diagram"))
  })
  
  output$tbl_desc <- renderDT({
    req(desc_tbl())
    DT::datatable(desc_tbl(), options = list(dom = "t", scrollX = TRUE))
  })
  
  output$tbl_top10 <- renderDT({
    req(top10_tbl())
    DT::datatable(top10_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$tbl_bottom10 <- renderDT({
    req(bottom10_tbl())
    DT::datatable(bottom10_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$plt_prev_hist <- renderPlotly({
    req(df(), input$bins_prev)
    p <- ggplot(df(), aes(x = prevalensi_kusta)) +
      geom_histogram(bins = input$bins_prev, fill = "pink", color = "white", alpha = 0.85) +
      labs(title = "Distribusi Prevalensi Kusta", x = "Prevalensi", y = "Jumlah Kab/Kota") +
      theme_minimal(base_size = 13)
    plotly::ggplotly(p, tooltip = c("x", "count"))
  })
  
  output$tbl_ir <- renderDT({
    req(ir_tbl())
    DT::datatable(ir_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$plt_ir_hist <- renderPlotly({
    req(ir_tbl(), input$bins_ir)
    p <- ggplot(ir_tbl(), aes(x = incidence_rate)) +
      geom_histogram(bins = input$bins_ir, fill = "pink", color = "white", alpha = 0.85) +
      labs(title = "Incidence Rate Kusta", x = "Incidence Rate", y = "Jumlah Kab/Kota") +
      theme_minimal(base_size = 12)
    plotly::ggplotly(p, tooltip = c("x", "count"))
  })
  
  output$tbl_rft <- renderDT({
    req(rft_tbl())
    DT::datatable(rft_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$plt_rft_pb <- renderPlotly({
    req(rft_tbl(), input$bins_rft_pb)
    p <- ggplot(rft_tbl(), aes(x = rft_pb)) +
      geom_histogram(bins = input$bins_rft_pb, fill = "pink", color = "white", alpha = 0.85) +
      labs(title = "RFT PB", x = "RFT PB", y = "Jumlah Kab/Kota") +
      theme_minimal(base_size = 12)
    plotly::ggplotly(p, tooltip = c("x", "count"))
  })
  
  output$plt_rft_mb <- renderPlotly({
    req(rft_tbl(), input$bins_rft_mb)
    p <- ggplot(rft_tbl(), aes(x = rft_mb)) +
      geom_histogram(bins = input$bins_rft_mb, fill = "pink", color = "white", alpha = 0.85) +
      labs(title = "RFT MB", x = "RFT MB", y = "Jumlah Kab/Kota") +
      theme_minimal(base_size = 12)
    plotly::ggplotly(p, tooltip = c("x", "count"))
  })
  
  output$txt_pois <- renderPrint({
    req(pois_pack())
    cat("Model Poisson:\n")
    print(summary(pois_pack()$m))
    cat("\nDispersion (Pearson Chi^2 / df): ", round(pois_pack()$disp, 4), "\n")
  })
  
  output$tbl_nb_rr <- renderDT({
    req(nb_rr_tbl())
    DT::datatable(nb_rr_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$gg_kusta_cont <- renderPlot({
    req(jatim_sf_full())
    ggplot() +
      geom_sf(data = jatim_sf_full(), aes(fill = prevalensi_kusta), color = NA) +
      theme_bw() +
      scale_fill_gradient(low = "pink", high = "black") +
      theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
      labs(title = "Peta Kusta Jawa Timur (Kontinu)", fill = "Prevalensi")
  })
  
  output$gg_kusta_disc <- renderPlot({
    req(jatim_sf_full())
    sf0 <- jatim_sf_full()
    breaks <- c(input$b0, input$b1, input$b2, input$b3, input$b4)
    labels <- c("Very Low", "Low", "High", "Very High")
    sf0$kusta_disc <- cut(sf0$prevalensi_kusta, breaks = breaks, labels = labels, right = TRUE, include.lowest = TRUE)
    
    ggplot() +
      geom_sf(data = sf0, aes(fill = kusta_disc), color = NA) +
      theme_bw() +
      scale_fill_manual(values = c("Very Low" = "pink", "Low" = "magenta", "High" = "purple", "Very High" = "black")) +
      theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
      labs(title = "Kategori Kusta Jawa Timur (Diskrit)", fill = "Kategori")
  })
  
  output$map_prev <- renderLeaflet({
    req(jatim_sf_full(), input$map_mode)
    sf0 <- jatim_sf_full()
    
    breaks <- c(0, 0.11075, 0.22600, 0.56175, 2.06300)
    labels <- c("Very Low", "Low", "High", "Very High")
    sf0$kusta_disc <- cut(sf0$prevalensi_kusta, breaks = breaks, labels = labels, right = TRUE, include.lowest = TRUE)
    
    if (input$map_mode == "cont") {
      pal <- leaflet::colorNumeric(palette = pal_bupk(), domain = sf0$prevalensi_kusta, na.color = "#D1D5DB")
      leaflet(sf0) %>%
        addProviderTiles("CartoDB.Positron") %>%
        addPolygons(
          fillColor = ~pal(prevalensi_kusta),
          fillOpacity = 0.85,
          color = "#111827",
          weight = 1,
          popup = ~paste0("<b>", NAME_2, "</b><br>Prevalensi: ", prevalensi_kusta)
        ) %>%
        addLegend("bottomright", pal = pal, values = ~prevalensi_kusta, title = "Prevalensi")
    } else {
      pal <- leaflet::colorFactor(
        palette = c("Very Low" = "#EC4899", "Low" = "#C026D3", "High" = "#7C3AED", "Very High" = "#0B1F4B"),
        domain = sf0$kusta_disc,
        na.color = "#D1D5DB"
      )
      leaflet(sf0) %>%
        addProviderTiles("CartoDB.Positron") %>%
        addPolygons(
          fillColor = ~pal(kusta_disc),
          fillOpacity = 0.85,
          color = "#111827",
          weight = 1,
          popup = ~paste0(
            "<b>", NAME_2, "</b><br>Kategori: ", kusta_disc,
            "<br>Prevalensi: ", prevalensi_kusta
          )
        ) %>%
        addLegend("bottomright", pal = pal, values = ~kusta_disc, title = "Kategori")
    }
  })
  
  output$txt_moran <- renderPrint({
    req(moran_global())
    moran_global()
  })
  
  output$tbl_lisa <- renderDT({
    req(lisa_tbl())
    DT::datatable(lisa_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$map_lisa <- renderLeaflet({
    req(jatim_sf_kusta(), lisa_tbl())
    sf0 <- jatim_sf_kusta()
    lt <- lisa_tbl()
    
    sf0$quad <- factor(
      lt$quad,
      levels = c("High-High", "Low-Low", "High-Low (Outlier)", "Low-High (Outlier)", "Not significant")
    )
    
    pal <- leaflet::colorFactor(
      palette = c(
        "High-High" = "#EC4899",
        "Low-Low" = "#C026D3",
        "High-Low (Outlier)" = "#7C3AED",
        "Low-High (Outlier)" = "#0B1F4B",
        "Not significant" = "#9CA3AF"
      ),
      domain = levels(sf0$quad)
    )
    
    leaflet(sf0) %>%
      addProviderTiles("CartoDB.Positron") %>%
      addPolygons(
        fillColor = ~pal(as.character(quad)),
        fillOpacity = 0.85,
        color = "#111827",
        weight = 1,
        popup = ~paste0(
          "<b>", NAME_2, "</b><br>Cluster: ", quad,
          "<br>Prevalensi: ", prevalensi_kusta
        )
      ) %>%
      addLegend("bottomright", pal = pal, values = ~quad, title = "LISA (p<=0.05)")
  })
  
  output$tbl_dsr <- renderDT({
    req(dsr_tbl())
    DT::datatable(dsr_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$tbl_compare <- renderDT({
    req(compare_tbl())
    DT::datatable(compare_tbl(), options = list(pageLength = 10, scrollX = TRUE))
  })
  
  output$plt_hist_dsr <- renderPlotly({
    req(compare_tbl(), input$bins_dsr)
    hist_df <- compare_tbl() %>%
      dplyr::select(kabupaten_kota, crude_ir, DSR) %>%
      tidyr::pivot_longer(cols = c(crude_ir, DSR), names_to = "jenis", values_to = "rate") %>%
      dplyr::mutate(jenis = dplyr::recode(jenis, crude_ir = "Crude Rate", DSR = "DSR (Standardized)"))
    
    p <- ggplot(hist_df, aes(x = rate, fill = jenis)) +
      geom_histogram(bins = input$bins_dsr, alpha = 0.7, color = "white") +
      facet_wrap(~jenis, scales = "free_x") +
      labs(title = "Crude Rate vs DSR", x = "Rate per 1.000", y = "Jumlah Kab/Kota") +
      theme_minimal(base_size = 12) +
      theme(legend.position = "none")
    
    plotly::ggplotly(p, tooltip = c("x", "count"))
  })
  
  # EXPORT
  output$dl_excel_raw <- downloadHandler(
    filename = function() paste0("excel_raw_", Sys.Date(), ".csv"),
    content = function(file) {
      req(excel_raw())
      write.csv(excel_raw(), file, row.names = FALSE)
    }
  )
  
  output$dl_clean <- downloadHandler(
    filename = function() paste0("jatim_join_", Sys.Date(), ".csv"),
    content = function(file) {
      req(jatim_sf_full())
      write.csv(jatim_sf_full() %>% sf::st_drop_geometry(), file, row.names = FALSE)
    }
  )
  
  output$dl_lisa <- downloadHandler(
    filename = function() paste0("lisa_", Sys.Date(), ".csv"),
    content = function(file) {
      req(lisa_tbl())
      write.csv(lisa_tbl(), file, row.names = FALSE)
    }
  )
  
  output$dl_dsr <- downloadHandler(
    filename = function() paste0("dsr_", Sys.Date(), ".csv"),
    content = function(file) {
      req(dsr_tbl())
      write.csv(dsr_tbl(), file, row.names = FALSE)
    }
  )
}

shinyApp(ui, server)