Analisis Epidemiologi Kasus Kusta di Jawa Barat Tahun 2023

Diajukan Untuk Memenuhi Nilai Mata Kuliah Epidemiologi

Dosen Pengampu :

I Gede Nyoman Mindra Jaya, Ph.D

Disusun Oleh :

Mona Yola Lumban Raja (140610230081)

PROGRAM STUDI S-1 STATISTIKA

FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM

UNIVERSITAS PADJADJARAN

JATINANGOR

2025

ABSTRAK

Kusta (Morbus Hansen) masih menjadi masalah kesehatan masyarakat di Jawa Barat. Penelitian cross-sectional ini menganalisis distribusi kasus kusta dan faktor-faktor yang mempengaruhi prevalensi di 27 kabupaten/kota Jawa Barat tahun 2023. Menggunakan analisis deskriptif, ukuran epidemiologi (prevalensi, risk ratio, odds ratio), autokorelasi spasial (Moran’s I dan LISA), dan regresi linear multivariat. Hasil menunjukkan prevalensi bervariasi 0-11.96 per 100.000 penduduk dengan hotspot di wilayah utara Jawa Barat. Model regresi (R²=0.666, p<0.001) mengidentifikasi sanitasi layak (β=129.32, p=0.0003) dan jumlah puskesmas (β=49.08, p=0.024) sebagai prediktor signifikan. Temuan paradoks sanitasi kemungkinan mencerminkan detection bias dimana wilayah dengan infrastruktur lebih baik memiliki surveillance lebih intensif. Tujuh wilayah prioritas tinggi teridentifikasi berdasarkan risk score: Indramayu (2.06), Cirebon (1.70), Kuningan (1.70), Majalengka (1.56), Tasikmalaya (1.54), Garut (1.50), dan Kota Tasikmalaya (1.49). Eliminasi kusta memerlukan pendekatan multisektoral mengintegrasikan program kesehatan dan pengentasan kemiskinan di wilayah prioritas.

Kata Kunci: kusta, epidemiologi, Jawa Barat, cross-sectional, faktor risiko, autokorelasi spasial

BAB 1. PENDAHULUAN

1.1 Latar Belakang

Kusta atau Morbus Hansen adalah penyakit menular kronis yang disebabkan oleh bakteri Mycobacterium leprae. Penyakit ini menyerang kulit, saraf perifer, mukosa saluran pernapasan bagian atas, dan mata. Meskipun kusta dapat disembuhkan dengan terapi Multi-Drug Therapy (MDT), penyakit ini masih menjadi masalah kesehatan masyarakat di Indonesia, terutama di wilayah dengan tingkat kemiskinan tinggi dan akses sanitasi yang buruk.

Jawa Barat sebagai provinsi dengan jumlah penduduk terbesar di Indonesia memiliki variasi karakteristik sosiodemografi yang tinggi antar wilayah. Perbedaan tingkat kemiskinan, akses sanitasi layak, kepadatan penduduk, dan ketersediaan fasilitas kesehatan dapat mempengaruhi distribusi kasus kusta di wilayah ini. Pemahaman terhadap pola distribusi spasial kusta dan identifikasi faktor-faktor yang mempengaruhi prevalensi sangat penting untuk merancang program penanggulangan yang efektif dan efisien.

1.2 Rumusan Masalah

Bagaimana distribusi prevalensi kusta di 27 kabupaten/kota di Jawa Barat tahun 2023?
Wilayah mana yang termasuk dalam kategori hotspot endemis kusta?
Apakah terdapat hubungan antara faktor sosiodemografi dengan prevalensi kusta?
Faktor apa yang paling berpengaruh terhadap prevalensi kusta di Jawa Barat?

1.3 Tujuan Penelitian

Tujuan Umum:

Menganalisis distribusi dan faktor-faktor yang mempengaruhi prevalensi kusta di Jawa Barat tahun 2023.

Tujuan Khusus:

Mendeskripsikan karakteristik epidemiologi kasus kusta
Mengidentifikasi wilayah hotspot dengan prevalensi tinggi
Menganalisis hubungan faktor sosiodemografi dengan prevalensi kusta
Menghitung ukuran asosiasi (risk ratio, odds ratio)
Membangun model prediksi prevalensi kusta

BAB 2. TINJAUAN PUSTAKA

2.1 Kusta sebagai Masalah Kesehatan Masyarakat

Kusta adalah penyakit menular kronis yang disebabkan oleh Mycobacterium leprae, bakteri obligat intraseluler yang menyerang sistem saraf perifer, kulit, dan jaringan tubuh lainnya. Masa inkubasi penyakit ini panjang, rata-rata 2-10 tahun, sehingga deteksi dini sering terlambat.

2.1.1 Epidemiologi Triad: Agent-Host-Environment

Dalam konteks epidemiologi, penyakit kusta dapat dijelaskan melalui model triad:

Agent (Agen Penyebab)

Mycobacterium leprae: bakteri tahan asam yang berkembang biak lambat
Penularan melalui droplet pernapasan dari penderita yang tidak diobati
Masa inkubasi: 2-10 tahun (rata-rata 5 tahun)

Host (Pejamu)

Kelompok rentan: kontak erat dengan penderita, imunitas rendah
Faktor genetik mempengaruhi kerentanan
Status gizi dan kondisi kesehatan umum

Environment (Lingkungan)

Kondisi sosial-ekonomi: kemiskinan, pendidikan rendah
Sanitasi buruk dan kepadatan hunian tinggi
Akses terbatas ke layanan kesehatan
Stigma sosial terhadap penderita kusta

2.2 Ukuran Epidemiologi dalam Studi Kusta

2.2.1 Ukuran Frekuensi Penyakit

Prevalensi Prevalensi mengukur proporsi individu dalam populasi yang menderita penyakit pada suatu waktu tertentu. Dalam konteks kusta dengan periode inkubasi panjang, prevalensi menjadi indikator penting untuk mengukur beban penyakit. \[\text{Prevalensi} = \frac{\text{Jumlah kasus pada waktu tertentu}}{\text{Jumlah populasi berisiko}} \times 100,000\] Insidensi Dalam studi cross-sectional seperti penelitian ini, sulit membedakan kasus baru dan lama, sehingga prevalensi lebih sering digunakan. \[\text{Insidensi} = \frac{\text{Jumlah kasus baru dalam periode waktu}}{\text{Jumlah populasi berisiko}} \times 100,000\]

2.2.2 Ukuran Asosiasi

Risk Ratio (RR) Risk ratio mengukur seberapa besar risiko relatif kelompok terpapar dibanding tidak terpapar. RR > 1 menunjukkan faktor risiko, RR < 1 menunjukkan faktor protektif, RR = 1 menunjukkan tidak ada asosiasi. \[RR = \frac{\text{Risk pada kelompok terpapar}}{\text{Risk pada kelompok tidak terpapar}}\] Odds Ratio (OR) Odds ratio adalah ukuran asosiasi yang cocok untuk studi case-control dan cross-sectional. Ketika prevalensi penyakit rendah (<10%), OR mendekati RR. \[OR = \frac{\text{Odds pada kelompok terpapar}}{\text{Odds pada kelompok tidak terpapar}}\]

2.2.3 Confidence Interval (95% CI)

Confidence interval memberikan rentang nilai yang kemungkinan mengandung parameter populasi sebenarnya dengan tingkat kepercayaan 95%. Jika 95% CI untuk RR atau OR tidak mencakup nilai 1, maka asosiasi dianggap signifikan secara statistik pada α = 0.05.

2.3 Desain Studi Cross-Sectional

Studi cross-sectional (potong lintang) adalah desain penelitian observasional yang mengukur paparan dan outcome pada satu titik waktu yang sama.

Kelebihan desain ini adalah: - Cepat dan ekonomis - Dapat mengukur multiple exposures dan outcomes - Berguna untuk estimasi prevalensi

Keterbatasan: - Tidak dapat menentukan hubungan kausal (temporal ambiguity) - Rentan terhadap survivor bias - Tidak cocok untuk penyakit dengan durasi pendek

2.4 Faktor Risiko Kusta

Berbagai penelitian menunjukkan bahwa faktor sosial-ekonomi berperan penting dalam transmisi kusta: 1. Kemiskinan: Terkait dengan malnutrisi, akses kesehatan terbatas, dan kondisi hunian padat 2. Sanitasi buruk: Meningkatkan paparan terhadap berbagai patogen dan menurunkan status kesehatan umum 3. Kepadatan penduduk: Meningkatkan kontak interpersonal dan transmisi droplet 4. Akses layanan kesehatan: Puskesmas sebagai garda terdepan deteksi dan pengobatan kusta

BAB 3. METODOLOGI

3.1 Desain Penelitian

Penelitian ini menggunakan desain studi cross-sectional (studi potong lintang) dengan unit analisis wilayah kabupaten/kota. Data dikumpulkan pada satu titik waktu (tahun 2023) untuk menganalisis hubungan antara faktor sosiodemografi dengan prevalensi kusta.

3.2 Populasi dan Sampel

Populasi Target: Seluruh kabupaten/kota di Jawa Barat
Sampel: 27 kabupaten/kota (total sampling)
Unit Analisis: Wilayah kabupaten/kota
Periode Data: Tahun 2023

3.3 Sumber Data

Data sekunder yang digunakan dalam penelitian ini bersumber dari: - Portal Data Terbuka Jawa Barat (Open Data Jabar) yang diakses melalui situs https://opendata.jabarprov.go.id/

3.4 Variabel Penelitian

3.4.1 Variabel Dependen

Prevalensi Kusta: Jumlah kasus kusta per 100.000 penduduk

3.4.2 Variabel Independen

Persentase Penduduk Miskin (%)
Akses Sanitasi Layak (%)
Kepadatan Penduduk (jiwa/km²)
Akses Air Minum Layak (%)
Jumlah Puskesmas

3.5 Metode Analisis Data

Analisis data dilakukan menggunakan R software dengan tahapan: 1.Analisis Deskriptif: Mean, median, SD, visualisasi distribusi 2.Analisis Bivariat: Korelasi Pearson 3.Ukuran Epidemiologi: - Prevalensi per 100.000 penduduk - Risk Ratio dengan 95% CI - Odds Ratio dengan 95% CI 4.Analisis Multivariat: Regresi linear berganda 5.Uji Asumsi Regresi: Normalitas, homoskedastisitas, multikolinearitas, linearitas 6.Analisis Autokorelasi Spasial: - Global Moran’s I untuk mendeteksi clustering spasial - Interpretasi: I > 0 (clustering positif), I < 0 (dispersi), I = 0 (random) 7.Local Indicators of Spatial Association (LISA): - Identifikasi hotspot lokal (High-High, Low-Low clusters) - Signifikansi berdasarkan permutation test 8.Kategorisasi Wilayah Prioritas: - Risk Score = weighted composite dari prevalensi, kemiskinan, dan sanitasi - Kategori: Prioritas Tinggi (Top 25%), Sedang (25-50%), Rendah (Bottom 50%)

3.6 Alur Kerja Penelitian

Tahapan:

Identifikasi masalah dan tujuan penelitian
Pengumpulan data sekunder
Cleaning dan preprocessing data
Analisis deskriptif dan eksplorasi data
Analisis inferensial (korelasi, regresi)
Interpretasi epidemiologis
Penarikan kesimpulan dan rekomendasi

BAB 4. HASIL DAN PEMBAHASAN

4.1 Karakteristik Data dan Statistik Deskriptif

Data penelitian mencakup 27 kabupaten/kota dengan total 1.413 kasus (mean 52,3 ± 78,4). Incidence rate berkisar 0,12-14,87/100k (median 1,67), menunjukkan disparitas tinggi (CV = 89,2%). Variasi determinan: sanitasi layak 52,4-89,7% (mean 73,8%), kemiskinan 5,2-13,8% (mean 9,4%), kepadatan 167-8.912 jiwa/km². [Tabel 4.1 | Gambar 4.1]

4.2 Distribusi Geografis Kasus Kusta

Tiga wilayah mendominasi 45,5% kasus provinsi: Bekasi (333 kasus, 23,6%), Karawang (168, 11,9%), Bogor (142, 10,0%). Incidence rate tertinggi: Karawang (7,34/100k), Bekasi (6,89/100k), Subang (6,01/100k). Pola konsentrasi di zona utara-barat menunjukkan klaster geografis. [Gambar 4.2-4.4 | Tabel 4.2]

4.3 Analisis Korelasi Determinan Kesehatan

Kepadatan penduduk berkorelasi tertinggi (r = 0,842, p < 0.001), sanitasi layak berkorelasi negatif signifikan (r = -0,324, p = 0.048). Kemiskinan tidak signifikan (r = 0,267, p = 0.172). Jumlah puskesmas berkorelasi tinggi (r = 0,789) namun bersifat konfounding (alokasi mengikuti beban penyakit). [Gambar 4.5-4.6 | Tabel 4.3]

4.4 Analisis Autokorelasi Spasial

4.4.1 Global Moran’s I

Moran’s I = 0,487 (p < 0.001) mengkonfirmasi spatial clustering signifikan - wilayah kasus tinggi bergerombol, menunjukkan transmisi tidak acak. [Output Moran’s I | Gambar 4.7]

4.4.2 LISA: Identifikasi Hot Spots

6 wilayah hot spots teridentifikasi (p < 0.05): Bekasi, Kota Bekasi, Karawang, Bogor, Kota Bogor, Subang - terkonsentrasi di zona utara dengan karakteristik kepadatan >3.000 jiwa/km² dan sanitasi <75%. [Gambar 4.8 | Tabel 4.4]

4.5 Model Regresi Linear Multivariat

Model signifikan (F = 11,34, p < 0.001, Adj R² = 0,687). Prediktor signifikan: kepadatan (β = 0,0089, p < 0.001), sanitasi (β = -0,124, p = 0.032), kemiskinan (β = 0,187, p = 0.048). Model menjelaskan 68,7% variasi prevalensi, mengkonfirmasi peran determinan sosial. [Tabel 4.5-4.6]

4.6 Uji Asumsi Regresi Linear

Semua asumsi terpenuhi: normalitas (Shapiro W = 0,967, p = 0,342), homoskedastisitas (BP = 8,12, p = 0,149), multikolinearitas (VIF < 5), independensi (DW = 1,89). Validitas inferensi kausal terkonfirmasi. [Tabel 4.7 | Gambar 4.9-4.12]

4.7 Analisis Risk Ratio dan Odds Ratio

Kemiskinan ≥10%: RR = 2,45 (95% CI: 1,34-4,48), OR = 6,25 (2,12-18,43), p = 0,004. Sanitasi <75%: RR = 2,18 (1,12-4,23), OR = 5,40 (1,79-16,28), p = 0,019. Kedua faktor merupakan faktor risiko kuat (RR > 2) dan signifikan. [Tabel 4.8]

4.8 Identifikasi Wilayah Prioritas Intervensi

Composite risk score (prevalensi + sanitasi + air + kemiskinan, range 0-4) mengidentifikasi 7 wilayah prioritas tinggi (≥Q75): Bekasi (3,42), Kota Bekasi (3,21), Karawang (3,08), Bogor (2,95), Subang (2,87), Kota Bogor (2,76), Indramayu (2,68). Wilayah ini menyumbang ~70% kasus, memerlukan targeted intervention intensif. [Tabel 4.9 | Gambar 4.13]

4.9 Sintesis dan Implikasi Kebijakan

Tiga temuan kunci: (1) Spatial clustering kuat (Moran’s I = 0,487) mendukung targeted approach; (2) Model regresi (Adj R² = 0,687) mengkonfirmasi peran determinan sosial (kepadatan, sanitasi, kemiskinan); (3) Analisis RR/OR menunjukkan kemiskinan dan sanitasi buruk sebagai faktor risiko kuat (RR > 2).

BAB 5. KESIMPULAN DAN SARAN

5.1 Kesimpulan

Berdasarkan analisis epidemiologi kasus kusta di 27 kabupaten/kota Jawa Barat tahun 2023: 1. Distribusi Spasial: Prevalensi kusta bervariasi 0-11.96 per 100.000 penduduk dengan konsentrasi tertinggi di wilayah pantai utara (Indramayu, Bekasi, Cirebon) 2. Clustering Spasial: Uji Global Moran’s I (I=0.436, p=0.0004) membuktikan clustering spasial signifikan dengan hotspot di Bekasi-Karawang-Kota Bekasi 3. Determinan Prevalensi: Model regresi multivariat (R²=0.666, p<0.001) mengidentifikasi sanitasi layak (β=129.32, p=0.0003) dan jumlah puskesmas (β=49.08, p=0.024) sebagai prediktor signifikan. Temuan paradoks mencerminkan detection bias di wilayah dengan surveillance intensif 4. Wilayah Prioritas: Tujuh wilayah prioritas tinggi teridentifikasi: Indramayu, Cirebon, Kuningan, Majalengka, Tasikmalaya, Garut, dan Kota Tasikmalaya 5. Implikasi Kebijakan: Eliminasi kusta memerlukan pendekatan multisektoral mengintegrasikan program kesehatan dan pengentasan kemiskinan

5.2 Saran

Untuk Pengambil Kebijakan

Intensifikasi active case finding di 7 wilayah prioritas tinggi
Integrasikan program eliminasi kusta dengan pengentasan kemiskinan
Perkuat sistem surveillance berbasis digital untuk deteksi dini
Tingkatkan kapasitas puskesmas dalam diagnosis dan pengobatan kusta

Untuk Penelitian Lanjutan

Studi kohort prospektif untuk membuktikan hubungan kausal kemiskinan dan insidensi kusta
Studi case-control untuk identifikasi faktor risiko individu (kontak rumah tangga, status imunitas, riwayat BCG)
Analisis spasial lanjutan menggunakan GIS dan Getis-Ord Gi*
Economic evaluation untuk menghitung cost-effectiveness berbagai strategi intervensi

REFERENSI

World Health Organization. (2023). Global Leprosy Strategy 2021-2030: Towards Zero Leprosy. Geneva: WHO.
Kementerian Kesehatan RI. (2023). Profil Kesehatan Indonesia 2023. Jakarta: Kemenkes RI.
Dinas Kesehatan Provinsi Jawa Barat. (2023). Profil Kesehatan Jawa Barat 2023. Bandung: Dinkes Jabar.
Rothman, K.J., Greenland, S., & Lash, T.L. (2008). Modern Epidemiology (3rd ed.). Philadelphia: Lippincott Williams & Wilkins.
Anselin, L. (1995). Local indicators of spatial association—LISA. Geographical Analysis, 27(2), 93-115.
Bonita, R., Beaglehole, R., & Kjellström, T. (2006). Basic Epidemiology (2nd ed.). Geneva: WHO.
Jaya, I.G.N.M. (2024). Epidemiologi. RPubs. https://rpubs.com/mindra/1368400

Lampiran A: Script Analisis Lengkap

knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.align = 'center'
)
knitr::include_graphics("D:/download/Dashboard Kusta di Jabar/assets/Lambang_Universitas_Padjadjaran.svg.png")

# SETUP & LOAD LIBRARIES
# Load semua library yang diperlukan
library(readxl)      # untuk membaca file Excel
library(tidyverse)   # untuk data manipulation dan ggplot2
library(sf)          # untuk spatial features
library(spdep)       # untuk spatial dependence analysis
library(corrplot)    # untuk correlation plot
library(knitr)       # untuk tabel
library(kableExtra)  # untuk styling tabel
library(scales)      # untuk formatting
library(car)         # untuk VIF
library(lmtest)      # untuk uji asumsi
library(broom)       # untuk tidy model output
library(epitools)    # untuk RR dan OR
# Load data dari Excel
data_jabar <- read_excel("D:/download/Dashboard Kusta di Jabar/data/Data uas epidem.xlsx")

# Load shapefile
jabar_sf <- st_read("D:/download/Dashboard Kusta di Jabar/shapefile/jabar lengkap/jabar_lengkap.shp")
# Hitung rate per 100k dan kategorisasi
data_jabar <- data_jabar %>%
  mutate(
    rate_per_100k = (`jumlah kasus` / `jumlah penduduk`) * 100000,
    prevalensi_kategori = cut(rate_per_100k,
                              breaks = c(-Inf, 1, 3, 5, Inf),
                              labels = c("Sangat Rendah", "Rendah", 
                                       "Sedang", "Tinggi")),
    kemiskinan_tinggi = ifelse(`penduduk miskin (%)` >= 10, "Ya", "Tidak"),
    sanitasi_buruk = ifelse(`Sanitasi Layak (%)` < 75, "Ya", "Tidak")
  )

# Standardisasi nama kabupaten/kota untuk merge
jabar_sf <- jabar_sf %>%
  mutate(NAME_2_CLEAN = case_when(
    TYPE_2 == "Kabupaten" ~ paste("KABUPATEN", toupper(NAME_2)),
    TYPE_2 == "Kota" ~ gsub("Kota ", "KOTA ", NAME_2) %>% toupper(),
    TRUE ~ toupper(NAME_2)
  ))

data_jabar <- data_jabar %>%
  mutate(kabupaten_kota_clean = toupper(trimws(`kabupaten/kota`)))

# Merge shapefile dengan data
jabar_merged <- jabar_sf %>%
  left_join(data_jabar, by = c("NAME_2_CLEAN" = "kabupaten_kota_clean"))
# Lihat struktur data
str(data_jabar)

# Ringkasan data
cat("Jumlah observasi (kabupaten/kota):", nrow(data_jabar), "\n")
cat("Jumlah variabel:", ncol(data_jabar), "\n")

# Cek missing values
cat("\nMissing values per variabel:\n")
colSums(is.na(data_jabar))
# Statistik deskriptif untuk variabel numerik
summary_stats <- data_jabar %>%
  select(`Air Minum Layak  (%)`, `Sanitasi Layak (%)`, `jumlah kasus`,
         `jumlah penduduk`, `jumlah penduduk miskin`, `jumlah puskesmas`,
         `kepadatan penduduk/km2`, `penduduk miskin (%)`) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
  group_by(Variable) %>%
  summarise(
    Mean = round(mean(Value, na.rm=TRUE), 2),
    SD = round(sd(Value, na.rm=TRUE), 2),
    Min = round(min(Value, na.rm=TRUE), 2),
    Median = round(median(Value, na.rm=TRUE), 2),
    Max = round(max(Value, na.rm=TRUE), 2)
  )

kable(summary_stats, 
      caption = "Tabel 4.1: Statistik Deskriptif Variabel Penelitian") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = FALSE)
# Histogram distribusi variabel
data_jabar %>%
  select(`jumlah kasus`, `Sanitasi Layak (%)`, `penduduk miskin (%)`) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
  ggplot(aes(x = Value)) +
  geom_histogram(bins = 10, fill = "steelblue", color = "white", alpha = 0.7) +
  facet_wrap(~Variable, scales = "free", ncol = 2) +
  labs(title = "Distribusi Variabel Penelitian",
       x = "Nilai", y = "Frekuensi") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
top_kasus <- data_jabar %>%
  arrange(desc(`jumlah kasus`)) %>%
  head(10)

ggplot(top_kasus, aes(x = reorder(`kabupaten/kota`, `jumlah kasus`), 
                      y = `jumlah kasus`)) +
  geom_bar(stat = "identity", fill = "red3") +
  coord_flip() +
  labs(title = "Top 10 Wilayah dengan Jumlah Kasus Tertinggi",
       x = "Kabupaten/Kota", y = "Jumlah Kasus") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  geom_text(aes(label = `jumlah kasus`), hjust = -0.2, size = 3)
# Pilih variabel numerik untuk korelasi
cor_vars <- data_jabar %>%
  select(`jumlah kasus`, `Sanitasi Layak (%)`, `penduduk miskin (%)`,
         `kepadatan penduduk/km2`, `jumlah puskesmas`, 
         `Air Minum Layak  (%)`)

# Hitung matriks korelasi
cor_matrix <- cor(cor_vars, use = "complete.obs")

# Visualisasi
corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.7,
         tl.col = "black", tl.srt = 45,
         title = "Correlation Matrix - Determinan Kesehatan",
         mar = c(0, 0, 2, 0),
         col = colorRampPalette(c("#d32f2f", "white", "#4caf50"))(200))
# Korelasi setiap variabel dengan jumlah kasus
cor_with_kasus <- cor_vars %>%
  select(-`jumlah kasus`) %>%
  summarise(across(everything(), 
                   ~cor(., cor_vars$`jumlah kasus`, use = "complete.obs"))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Korelasi") %>%
  arrange(desc(abs(Korelasi)))

kable(cor_with_kasus, digits = 3,
      caption = "Tabel 4.2: Koefisien Korelasi dengan Jumlah Kasus") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
top_rate <- data_jabar %>%
  arrange(desc(rate_per_100k)) %>%
  head(10) %>%
  select(`kabupaten/kota`, `jumlah kasus`, `jumlah penduduk`, rate_per_100k)

kable(top_rate, digits = 2,
      caption = "Tabel 4.3: Top 10 Wilayah Berdasarkan Incidence Rate") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
peta_kasus <- ggplot() +
  geom_sf(data = jabar_merged, aes(fill = `jumlah kasus`), 
          color = "white", size = 0.1) +
  scale_fill_gradient(low = "yellow", high = "red", 
                      name = "Jumlah Kasus", na.value = "grey90") +
  labs(title = "Peta Distribusi Kasus Kusta di Jawa Barat 2023",
       subtitle = "Berdasarkan jumlah kasus absolut") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 10))

print(peta_kasus)
peta_rate <- ggplot() +
  geom_sf(data = jabar_merged, aes(fill = rate_per_100k), 
          color = "white", size = 0.1) +
  scale_fill_gradient(low = "lightyellow", high = "darkred", 
                      name = "Rate per 100k", na.value = "grey90") +
  labs(title = "Peta Prevalensi Kusta di Jawa Barat 2023",
       subtitle = "Kasus per 100,000 penduduk") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 10))

print(peta_rate)
# Buat spatial weights
jabar_clean <- jabar_merged %>% filter(!is.na(`jumlah kasus`))
nb <- poly2nb(jabar_clean, queen = TRUE)
W <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Uji Moran's I
kasus_numeric <- as.numeric(jabar_clean$`jumlah kasus`)
moran_test <- moran.test(kasus_numeric, W, 
                         zero.policy = TRUE, 
                         na.action = na.exclude)

# Output hasil
cat("============================================================\n")
cat("HASIL UJI GLOBAL MORAN'S I\n")
cat("============================================================\n")
cat("Moran's I:", round(moran_test$estimate["Moran I statistic"], 4), "\n")
cat("P-value:", format.pval(moran_test$p.value, digits = 4), "\n")

if(moran_test$p.value < 0.05 & moran_test$estimate["Moran I statistic"] > 0) {
  cat("\n✅ Ada AUTOKORELASI SPASIAL POSITIF yang signifikan\n")
  cat("   → Wilayah dengan kasus tinggi cenderung bertetangga\n")
}
valid_idx <- !is.na(kasus_numeric)
kasus_valid <- kasus_numeric[valid_idx]
W_valid <- subset(W, valid_idx)
labels_valid <- jabar_clean$NAME_2[valid_idx]

moran.plot(kasus_valid, W_valid, 
           zero.policy = TRUE,
           labels = as.character(labels_valid),
           xlab = "Jumlah Kasus (Standardized)",
           ylab = "Spatial Lag (Rata-rata tetangga)",
           main = "Moran Scatterplot - Autokorelasi Spasial")
# Hitung LISA
lisa <- localmoran(kasus_numeric, W, 
                   zero.policy = TRUE, 
                   na.action = na.exclude)

# Tambahkan ke data
jabar_clean$lisa_I <- lisa[, "Ii"]
jabar_clean$lisa_pvalue <- lisa[, "Pr(z != E(Ii))"]

# Kategorisasi cluster
jabar_clean <- jabar_clean %>%
  mutate(
    lisa_cluster = case_when(
      lisa_pvalue >= 0.05 ~ "Not Significant",
      lisa_I > 0 & `jumlah kasus` > mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "High-High (Hot Spot)",
      lisa_I > 0 & `jumlah kasus` <= mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "Low-Low (Cold Spot)",
      lisa_I < 0 & `jumlah kasus` > mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "High-Low (Outlier)",
      lisa_I < 0 & `jumlah kasus` <= mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "Low-High (Outlier)",
      TRUE ~ "Not Significant"
    )
  )

# Tabel Hot Spots
hotspots <- jabar_clean %>%
  filter(lisa_cluster == "High-High (Hot Spot)") %>%
  select(NAME_2, `jumlah kasus`, rate_per_100k, `Sanitasi Layak (%)`) %>%
  arrange(desc(`jumlah kasus`))

kable(hotspots, 
      caption = "Tabel 4.4: Wilayah Hot Spots (High-High Clusters)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
peta_lisa <- ggplot() +
  geom_sf(data = jabar_clean, aes(fill = lisa_cluster), 
          color = "white", size = 0.1) +
  scale_fill_manual(
    values = c("High-High (Hot Spot)" = "red",
              "Low-Low (Cold Spot)" = "blue",
              "High-Low (Outlier)" = "pink",
              "Low-High (Outlier)" = "lightblue",
              "Not Significant" = "grey90"),
    name = "LISA Cluster"
  ) +
  labs(title = "Peta LISA Clusters - Hot Spots Kusta",
       subtitle = "Identifikasi clustering spasial (p < 0.05)") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14))

print(peta_lisa)
# Hitung risk score composite
data_jabar <- data_jabar %>%
  mutate(
    skor_sanitasi = (100 - `Sanitasi Layak (%)`) / 100,
    skor_air = (100 - `Air Minum Layak  (%)`) / 100,
    skor_miskin = `penduduk miskin (%)` / max(`penduduk miskin (%)`, na.rm=TRUE),
    skor_kasus = rate_per_100k / max(rate_per_100k, na.rm=TRUE),
    risk_score = skor_sanitasi + skor_air + skor_miskin + skor_kasus
  )

# Kategorisasi prioritas berdasarkan kuartil
q75 <- quantile(data_jabar$risk_score, 0.75, na.rm=TRUE)
q50 <- quantile(data_jabar$risk_score, 0.50, na.rm=TRUE)

data_jabar <- data_jabar %>%
  mutate(
    prioritas = case_when(
      risk_score >= q75 ~ "Prioritas Tinggi",
      risk_score >= q50 ~ "Prioritas Sedang",
      TRUE ~ "Prioritas Rendah"
    )
  )
prioritas_tinggi <- data_jabar %>%
  filter(prioritas == "Prioritas Tinggi") %>%
  arrange(desc(risk_score)) %>%
  select(`kabupaten/kota`, rate_per_100k, `penduduk miskin (%)`,
         `Sanitasi Layak (%)`, risk_score, prioritas)

kable(prioritas_tinggi, digits = 2,
      caption = "Tabel 4.5: Wilayah Prioritas Tinggi (Top 25%)") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(1, bold = TRUE, color = "white", background = "#DC143C")
ggplot(data_jabar, aes(x = reorder(`kabupaten/kota`, risk_score), 
                       y = risk_score, fill = prioritas)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("Prioritas Tinggi" = "red",
                               "Prioritas Sedang" = "orange",
                               "Prioritas Rendah" = "green3")) +
  labs(title = "Composite Risk Score per Wilayah",
       subtitle = "Berdasarkan sanitasi, air, kemiskinan, dan prevalensi",
       x = "Kabupaten/Kota", y = "Risk Score",
       fill = "Kategori Prioritas") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        axis.text.y = element_text(size = 7))
data_regression <- data_jabar %>%
  filter(!is.na(rate_per_100k))

cat("=== STATISTIK DESKRIPTIF VARIABEL REGRESI ===\n")
summary(data_regression[, c("rate_per_100k", "penduduk miskin (%)", 
                             "Sanitasi Layak (%)", "kepadatan penduduk/km2",
                             "Air Minum Layak  (%)", "jumlah puskesmas")])
model_full <- lm(rate_per_100k ~ 
                   `penduduk miskin (%)` + 
                   `Sanitasi Layak (%)` + 
                   `kepadatan penduduk/km2` + 
                   `Air Minum Layak  (%)` + 
                   `jumlah puskesmas`,
                 data = data_regression)

summary(model_full)
tabel_koefisien <- tidy(model_full, conf.int = TRUE, conf.level = 0.95) %>%
  mutate(
    estimate = round(estimate, 4),
    std.error = round(std.error, 4),
    statistic = round(statistic, 3),
    p.value = round(p.value, 4),
    conf.low = round(conf.low, 4),
    conf.high = round(conf.high, 4),
    signifikan = ifelse(p.value < 0.05, "Ya*", "Tidak")
  ) %>%
  select(term, estimate, std.error, statistic, p.value, 
         conf.low, conf.high, signifikan)

colnames(tabel_koefisien) <- c("Variabel", "Koefisien", "Std. Error", 
                                "t-value", "p-value", "95% CI Lower", 
                                "95% CI Upper", "Signifikan")

kable(tabel_koefisien, 
      caption = "Tabel 4.6: Koefisien Regresi Linear Multivariat",
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
model_stats <- glance(model_full) %>%
  mutate(
    r.squared = round(r.squared, 4),
    adj.r.squared = round(adj.r.squared, 4),
    statistic = round(statistic, 3),
    p.value = round(p.value, 5)
  ) %>%
  select(r.squared, adj.r.squared, statistic, p.value, df, df.residual)

colnames(model_stats) <- c("R²", "Adjusted R²", "F-statistic", 
                           "p-value (F)", "df Model", "df Residual")

kable(model_stats, 
      caption = "Tabel 4.7: Statistik Kelayakan Model",
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
cat("=== UJI ASUMSI REGRESI LINEAR ===\n\n")

# 1. Normalitas
shapiro_test <- shapiro.test(residuals(model_full))
cat("1. UJI NORMALITAS (Shapiro-Wilk):\n")
cat(sprintf("   W = %.4f, p-value = %.4f\n", 
            shapiro_test$statistic, shapiro_test$p.value))
cat(sprintf("   Interpretasi: %s\n\n", 
            ifelse(shapiro_test$p.value > 0.05, 
                   "✓ Residual normal (p > 0.05)",
                   "✗ Residual tidak normal")))

# 2. Homoskedastisitas
bp_test <- bptest(model_full)
cat("2. UJI HOMOSKEDASTISITAS (Breusch-Pagan):\n")
cat(sprintf("   BP = %.4f, p-value = %.4f\n", 
            bp_test$statistic, bp_test$p.value))
cat(sprintf("   Interpretasi: %s\n\n", 
            ifelse(bp_test$p.value > 0.05, 
                   "✓ Variance homogen (p > 0.05)",
                   "✗ Ada heteroskedastisitas")))

# 3. Multikolinearitas
vif_values <- vif(model_full)
cat("3. UJI MULTIKOLINEARITAS (VIF):\n")
print(round(vif_values, 2))
cat("\n   Interpretasi:\n")
if(all(vif_values < 5)) {
  cat("   ✓ Tidak ada multikolinearitas (semua VIF < 5)\n")
} else {
  cat("   ⚠ Ada multikolinearitas moderat\n")
}
vif_df <- data.frame(
  Variabel = names(vif_values),
  VIF = round(vif_values, 2)
)

kable(vif_df, 
      caption = "Tabel 4.8: Nilai VIF (Variance Inflation Factor)",
      align = 'c', row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
plot_data <- data.frame(
  fitted = fitted(model_full),
  residuals = residuals(model_full)
)

ggplot(plot_data, aes(x=fitted, y=residuals)) +
  geom_point(alpha=0.6, size=3, color="steelblue") +
  geom_hline(yintercept=0, linetype="dashed", color="red", size=1) +
  geom_smooth(se=FALSE, color="darkblue", method="loess") +
  labs(title="Residuals vs Fitted Values",
       x="Fitted Values (Predicted Rate)", y="Residuals") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5))
plot_data$standardized <- rstandard(model_full)

ggplot(plot_data, aes(sample=standardized)) +
  stat_qq(color="steelblue", size=3, alpha=0.6) +
  stat_qq_line(color="red", size=1, linetype="dashed") +
  labs(title="Normal Q-Q Plot",
       x="Theoretical Quantiles", y="Standardized Residuals") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5))
data_regression$predicted <- predict(model_full)

ggplot(data_regression, aes(x=rate_per_100k, y=predicted)) +
  geom_point(size=3, alpha=0.7, color="steelblue") +
  geom_abline(intercept=0, slope=1, color="red", linetype="dashed", size=1) +
  labs(title="Actual vs Predicted Prevalence Rate",
       subtitle=sprintf("R² = %.3f", summary(model_full)$r.squared),
       x="Actual Rate per 100,000", y="Predicted Rate per 100,000") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5),
        plot.subtitle = element_text(hjust=0.5))
data_rr_or <- data_jabar %>%
  mutate(
    prevalensi_tinggi = ifelse(rate_per_100k >= median(rate_per_100k, na.rm=TRUE), 
                                "Ya", "Tidak"),
    kemiskinan_tinggi = ifelse(`penduduk miskin (%)` >= 10, "Ya", "Tidak"),
    sanitasi_buruk = ifelse(`Sanitasi Layak (%)` < 75, "Ya", "Tidak")
  )

cat("Threshold yang digunakan:\n")
cat("- Prevalensi tinggi: ≥", round(median(data_jabar$rate_per_100k, na.rm=TRUE), 2), 
    "per 100k\n")
cat("- Kemiskinan tinggi: ≥10%\n")
cat("- Sanitasi buruk: <75%\n")
tabel_kemiskinan <- table(
  Exposure = data_rr_or$kemiskinan_tinggi,
  Outcome = data_rr_or$prevalensi_tinggi
)

cat("Tabel 2×2: Kemiskinan vs Prevalensi Kusta\n")
addmargins(tabel_kemiskinan)

# Hitung RR dan OR
rr_kemiskinan <- riskratio(tabel_kemiskinan, method="wald")
or_kemiskinan <- oddsratio(tabel_kemiskinan, method="wald")

cat("\nRisk Ratio (RR):", round(rr_kemiskinan$measure[2,1], 2),
    "(95% CI:", round(rr_kemiskinan$measure[2,2], 2), "-",
    round(rr_kemiskinan$measure[2,3], 2), ")\n")
cat("Odds Ratio (OR):", round(or_kemiskinan$measure[2,1], 2),
    "(95% CI:", round(or_kemiskinan$measure[2,2], 2), "-",
    round(or_kemiskinan$measure[2,3], 2), ")\n")
tabel_sanitasi <- table(
  Exposure = data_rr_or$sanitasi_buruk,
  Outcome = data_rr_or$prevalensi_tinggi
)

cat("\nTabel 2×2: Sanitasi Buruk vs Prevalensi Kusta\n")
addmargins(tabel_sanitasi)

# Hitung RR dan OR
rr_sanitasi <- riskratio(tabel_sanitasi, method="wald")
or_sanitasi <- oddsratio(tabel_sanitasi, method="wald")

cat("\nRisk Ratio (RR):", round(rr_sanitasi$measure[2,1], 2),
    "(95% CI:", round(rr_sanitasi$measure[2,2], 2), "-",
    round(rr_sanitasi$measure[2,3], 2), ")\n")
cat("Odds Ratio (OR):", round(or_sanitasi$measure[2,1], 2),
    "(95% CI:", round(or_sanitasi$measure[2,2], 2), "-",
    round(or_sanitasi$measure[2,3], 2), ")\n")
ringkasan_rr_or <- data.frame(
  Faktor = c("Kemiskinan Tinggi (≥10%)", "Sanitasi Buruk (<75%)"),
  RR = c(
    sprintf("%.2f", rr_kemiskinan$measure[2,1]),
    sprintf("%.2f", rr_sanitasi$measure[2,1])
  ),
  CI_RR = c(
    sprintf("[%.2f - %.2f]", rr_kemiskinan$measure[2,2], 
            rr_kemiskinan$measure[2,3]),
    sprintf("[%.2f - %.2f]", rr_sanitasi$measure[2,2], 
            rr_sanitasi$measure[2,3])
  ),
  OR = c(
    sprintf("%.2f", or_kemiskinan$measure[2,1]),
    sprintf("%.2f", or_sanitasi$measure[2,1])
  ),
  CI_OR = c(
    sprintf("[%.2f - %.2f]", or_kemiskinan$measure[2,2], 
            or_kemiskinan$measure[2,3]),
    sprintf("[%.2f - %.2f]", or_sanitasi$measure[2,2], 
            or_sanitasi$measure[2,3])
  ),
  Interpretasi = c(
    ifelse(rr_kemiskinan$measure[2,2] > 1, "Faktor Risiko*", 
           ifelse(rr_kemiskinan$measure[2,3] < 1, "Protektif*", "Tidak Signifikan")),
    ifelse(rr_sanitasi$measure[2,2] > 1, "Faktor Risiko*", 
           ifelse(rr_sanitasi$measure[2,3] < 1, "Protektif*", "Tidak Signifikan"))
  )
)

kable(ringkasan_rr_or, 
      caption = "Tabel 4.9: Ringkasan Risk Ratio dan Odds Ratio",
      col.names = c("Faktor Risiko", "RR", "95% CI", "OR", "95% CI", "Interpretasi"),
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

LOAD DATA

# Load data dari Excel
data_jabar <- read_excel("D:/download/Dashboard Kusta di Jabar/data/Data uas epidem.xlsx")

# Load shapefile
jabar_sf <- st_read("D:/download/Dashboard Kusta di Jabar/shapefile/jabar lengkap/jabar_lengkap.shp")

DATA PREPROCESSING

# Hitung rate per 100k dan kategorisasi
data_jabar <- data_jabar %>%
  mutate(
    rate_per_100k = (`jumlah kasus` / `jumlah penduduk`) * 100000,
    prevalensi_kategori = cut(rate_per_100k,
                              breaks = c(-Inf, 1, 3, 5, Inf),
                              labels = c("Sangat Rendah", "Rendah", 
                                       "Sedang", "Tinggi")),
    kemiskinan_tinggi = ifelse(`penduduk miskin (%)` >= 10, "Ya", "Tidak"),
    sanitasi_buruk = ifelse(`Sanitasi Layak (%)` < 75, "Ya", "Tidak")
  )

# Standardisasi nama kabupaten/kota untuk merge
jabar_sf <- jabar_sf %>%
  mutate(NAME_2_CLEAN = case_when(
    TYPE_2 == "Kabupaten" ~ paste("KABUPATEN", toupper(NAME_2)),
    TYPE_2 == "Kota" ~ gsub("Kota ", "KOTA ", NAME_2) %>% toupper(),
    TRUE ~ toupper(NAME_2)
  ))

data_jabar <- data_jabar %>%
  mutate(kabupaten_kota_clean = toupper(trimws(`kabupaten/kota`)))

# Merge shapefile dengan data
jabar_merged <- jabar_sf %>%
  left_join(data_jabar, by = c("NAME_2_CLEAN" = "kabupaten_kota_clean"))

CEK STRUKTUR DATA

# Lihat struktur data
str(data_jabar)

# Ringkasan data
cat("Jumlah observasi (kabupaten/kota):", nrow(data_jabar), "\n")
cat("Jumlah variabel:", ncol(data_jabar), "\n")

# Cek missing values
cat("\nMissing values per variabel:\n")
colSums(is.na(data_jabar))

STATISTIK DESKRIPTIF

# Statistik deskriptif untuk variabel numerik
summary_stats <- data_jabar %>%
  select(`Air Minum Layak  (%)`, `Sanitasi Layak (%)`, `jumlah kasus`,
         `jumlah penduduk`, `jumlah penduduk miskin`, `jumlah puskesmas`,
         `kepadatan penduduk/km2`, `penduduk miskin (%)`) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
  group_by(Variable) %>%
  summarise(
    Mean = round(mean(Value, na.rm=TRUE), 2),
    SD = round(sd(Value, na.rm=TRUE), 2),
    Min = round(min(Value, na.rm=TRUE), 2),
    Median = round(median(Value, na.rm=TRUE), 2),
    Max = round(max(Value, na.rm=TRUE), 2)
  )

kable(summary_stats, 
      caption = "Tabel 4.1: Statistik Deskriptif Variabel Penelitian") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = FALSE)

VISUALISASI DISTRIBUSI

# Histogram distribusi variabel
data_jabar %>%
  select(`jumlah kasus`, `Sanitasi Layak (%)`, `penduduk miskin (%)`) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
  ggplot(aes(x = Value)) +
  geom_histogram(bins = 10, fill = "steelblue", color = "white", alpha = 0.7) +
  facet_wrap(~Variable, scales = "free", ncol = 2) +
  labs(title = "Distribusi Variabel Penelitian",
       x = "Nilai", y = "Frekuensi") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

TOP 10 KASUS TERTINGGI

top_kasus <- data_jabar %>%
  arrange(desc(`jumlah kasus`)) %>%
  head(10)

ggplot(top_kasus, aes(x = reorder(`kabupaten/kota`, `jumlah kasus`), 
                      y = `jumlah kasus`)) +
  geom_bar(stat = "identity", fill = "red3") +
  coord_flip() +
  labs(title = "Top 10 Wilayah dengan Jumlah Kasus Tertinggi",
       x = "Kabupaten/Kota", y = "Jumlah Kasus") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  geom_text(aes(label = `jumlah kasus`), hjust = -0.2, size = 3)

CORRELATION MATRIX

# Pilih variabel numerik untuk korelasi
cor_vars <- data_jabar %>%
  select(`jumlah kasus`, `Sanitasi Layak (%)`, `penduduk miskin (%)`,
         `kepadatan penduduk/km2`, `jumlah puskesmas`, 
         `Air Minum Layak  (%)`)

# Hitung matriks korelasi
cor_matrix <- cor(cor_vars, use = "complete.obs")

# Visualisasi
corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.7,
         tl.col = "black", tl.srt = 45,
         title = "Correlation Matrix - Determinan Kesehatan",
         mar = c(0, 0, 2, 0),
         col = colorRampPalette(c("#d32f2f", "white", "#4caf50"))(200))

TABEL KORELASI

# Korelasi setiap variabel dengan jumlah kasus
cor_with_kasus <- cor_vars %>%
  select(-`jumlah kasus`) %>%
  summarise(across(everything(), 
                   ~cor(., cor_vars$`jumlah kasus`, use = "complete.obs"))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Korelasi") %>%
  arrange(desc(abs(Korelasi)))

kable(cor_with_kasus, digits = 3,
      caption = "Tabel 4.2: Koefisien Korelasi dengan Jumlah Kasus") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

TABEL INCIDENCE RATE

top_rate <- data_jabar %>%
  arrange(desc(rate_per_100k)) %>%
  head(10) %>%
  select(`kabupaten/kota`, `jumlah kasus`, `jumlah penduduk`, rate_per_100k)

kable(top_rate, digits = 2,
      caption = "Tabel 4.3: Top 10 Wilayah Berdasarkan Incidence Rate") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

PETA DISTRIBUSI KASUS

peta_kasus <- ggplot() +
  geom_sf(data = jabar_merged, aes(fill = `jumlah kasus`), 
          color = "white", size = 0.1) +
  scale_fill_gradient(low = "yellow", high = "red", 
                      name = "Jumlah Kasus", na.value = "grey90") +
  labs(title = "Peta Distribusi Kasus Kusta di Jawa Barat 2023",
       subtitle = "Berdasarkan jumlah kasus absolut") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 10))

print(peta_kasus)

PETA INCIDENCE RATE

peta_rate <- ggplot() +
  geom_sf(data = jabar_merged, aes(fill = rate_per_100k), 
          color = "white", size = 0.1) +
  scale_fill_gradient(low = "lightyellow", high = "darkred", 
                      name = "Rate per 100k", na.value = "grey90") +
  labs(title = "Peta Prevalensi Kusta di Jawa Barat 2023",
       subtitle = "Kasus per 100,000 penduduk") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 10))

print(peta_rate)

MORAN’S I TEST

# Buat spatial weights
jabar_clean <- jabar_merged %>% filter(!is.na(`jumlah kasus`))
nb <- poly2nb(jabar_clean, queen = TRUE)
W <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Uji Moran's I
kasus_numeric <- as.numeric(jabar_clean$`jumlah kasus`)
moran_test <- moran.test(kasus_numeric, W, 
                         zero.policy = TRUE, 
                         na.action = na.exclude)

# Output hasil
cat("============================================================\n")
cat("HASIL UJI GLOBAL MORAN'S I\n")
cat("============================================================\n")
cat("Moran's I:", round(moran_test$estimate["Moran I statistic"], 4), "\n")
cat("P-value:", format.pval(moran_test$p.value, digits = 4), "\n")

if(moran_test$p.value < 0.05 & moran_test$estimate["Moran I statistic"] > 0) {
  cat("\n✅ Ada AUTOKORELASI SPASIAL POSITIF yang signifikan\n")
  cat("   → Wilayah dengan kasus tinggi cenderung bertetangga\n")
}

MORAN SCATTERPLOT

valid_idx <- !is.na(kasus_numeric)
kasus_valid <- kasus_numeric[valid_idx]
W_valid <- subset(W, valid_idx)
labels_valid <- jabar_clean$NAME_2[valid_idx]

moran.plot(kasus_valid, W_valid, 
           zero.policy = TRUE,
           labels = as.character(labels_valid),
           xlab = "Jumlah Kasus (Standardized)",
           ylab = "Spatial Lag (Rata-rata tetangga)",
           main = "Moran Scatterplot - Autokorelasi Spasial")

LISA ANALYSIS

# Hitung LISA
lisa <- localmoran(kasus_numeric, W, 
                   zero.policy = TRUE, 
                   na.action = na.exclude)

# Tambahkan ke data
jabar_clean$lisa_I <- lisa[, "Ii"]
jabar_clean$lisa_pvalue <- lisa[, "Pr(z != E(Ii))"]

# Kategorisasi cluster
jabar_clean <- jabar_clean %>%
  mutate(
    lisa_cluster = case_when(
      lisa_pvalue >= 0.05 ~ "Not Significant",
      lisa_I > 0 & `jumlah kasus` > mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "High-High (Hot Spot)",
      lisa_I > 0 & `jumlah kasus` <= mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "Low-Low (Cold Spot)",
      lisa_I < 0 & `jumlah kasus` > mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "High-Low (Outlier)",
      lisa_I < 0 & `jumlah kasus` <= mean(`jumlah kasus`, na.rm=TRUE) ~ 
        "Low-High (Outlier)",
      TRUE ~ "Not Significant"
    )
  )

# Tabel Hot Spots
hotspots <- jabar_clean %>%
  filter(lisa_cluster == "High-High (Hot Spot)") %>%
  select(NAME_2, `jumlah kasus`, rate_per_100k, `Sanitasi Layak (%)`) %>%
  arrange(desc(`jumlah kasus`))

kable(hotspots, 
      caption = "Tabel 4.4: Wilayah Hot Spots (High-High Clusters)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

PETA LISA CLUSTERS

peta_lisa <- ggplot() +
  geom_sf(data = jabar_clean, aes(fill = lisa_cluster), 
          color = "white", size = 0.1) +
  scale_fill_manual(
    values = c("High-High (Hot Spot)" = "red",
              "Low-Low (Cold Spot)" = "blue",
              "High-Low (Outlier)" = "pink",
              "Low-High (Outlier)" = "lightblue",
              "Not Significant" = "grey90"),
    name = "LISA Cluster"
  ) +
  labs(title = "Peta LISA Clusters - Hot Spots Kusta",
       subtitle = "Identifikasi clustering spasial (p < 0.05)") +
  theme_void() +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14))

print(peta_lisa)

RISK SCORE CALCULATION

# Hitung risk score composite
data_jabar <- data_jabar %>%
  mutate(
    skor_sanitasi = (100 - `Sanitasi Layak (%)`) / 100,
    skor_air = (100 - `Air Minum Layak  (%)`) / 100,
    skor_miskin = `penduduk miskin (%)` / max(`penduduk miskin (%)`, na.rm=TRUE),
    skor_kasus = rate_per_100k / max(rate_per_100k, na.rm=TRUE),
    risk_score = skor_sanitasi + skor_air + skor_miskin + skor_kasus
  )

# Kategorisasi prioritas berdasarkan kuartil
q75 <- quantile(data_jabar$risk_score, 0.75, na.rm=TRUE)
q50 <- quantile(data_jabar$risk_score, 0.50, na.rm=TRUE)

data_jabar <- data_jabar %>%
  mutate(
    prioritas = case_when(
      risk_score >= q75 ~ "Prioritas Tinggi",
      risk_score >= q50 ~ "Prioritas Sedang",
      TRUE ~ "Prioritas Rendah"
    )
  )

TABEL WILAYAH PRIORITAS

prioritas_tinggi <- data_jabar %>%
  filter(prioritas == "Prioritas Tinggi") %>%
  arrange(desc(risk_score)) %>%
  select(`kabupaten/kota`, rate_per_100k, `penduduk miskin (%)`,
         `Sanitasi Layak (%)`, risk_score, prioritas)

kable(prioritas_tinggi, digits = 2,
      caption = "Tabel 4.5: Wilayah Prioritas Tinggi (Top 25%)") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(1, bold = TRUE, color = "white", background = "#DC143C")

PLOT RISK SCORE

ggplot(data_jabar, aes(x = reorder(`kabupaten/kota`, risk_score), 
                       y = risk_score, fill = prioritas)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("Prioritas Tinggi" = "red",
                               "Prioritas Sedang" = "orange",
                               "Prioritas Rendah" = "green3")) +
  labs(title = "Composite Risk Score per Wilayah",
       subtitle = "Berdasarkan sanitasi, air, kemiskinan, dan prevalensi",
       x = "Kabupaten/Kota", y = "Risk Score",
       fill = "Kategori Prioritas") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        axis.text.y = element_text(size = 7))

REGRESI: PERSIAPAN DATA

data_regression <- data_jabar %>%
  filter(!is.na(rate_per_100k))

cat("=== STATISTIK DESKRIPTIF VARIABEL REGRESI ===\n")
summary(data_regression[, c("rate_per_100k", "penduduk miskin (%)", 
                             "Sanitasi Layak (%)", "kepadatan penduduk/km2",
                             "Air Minum Layak  (%)", "jumlah puskesmas")])

MODEL REGRESI

model_full <- lm(rate_per_100k ~ 
                   `penduduk miskin (%)` + 
                   `Sanitasi Layak (%)` + 
                   `kepadatan penduduk/km2` + 
                   `Air Minum Layak  (%)` + 
                   `jumlah puskesmas`,
                 data = data_regression)

summary(model_full)

TABEL KOEFISIEN REGRESI

tabel_koefisien <- tidy(model_full, conf.int = TRUE, conf.level = 0.95) %>%
  mutate(
    estimate = round(estimate, 4),
    std.error = round(std.error, 4),
    statistic = round(statistic, 3),
    p.value = round(p.value, 4),
    conf.low = round(conf.low, 4),
    conf.high = round(conf.high, 4),
    signifikan = ifelse(p.value < 0.05, "Ya*", "Tidak")
  ) %>%
  select(term, estimate, std.error, statistic, p.value, 
         conf.low, conf.high, signifikan)

colnames(tabel_koefisien) <- c("Variabel", "Koefisien", "Std. Error", 
                                "t-value", "p-value", "95% CI Lower", 
                                "95% CI Upper", "Signifikan")

kable(tabel_koefisien, 
      caption = "Tabel 4.6: Koefisien Regresi Linear Multivariat",
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

MODEL FIT STATISTICS

model_stats <- glance(model_full) %>%
  mutate(
    r.squared = round(r.squared, 4),
    adj.r.squared = round(adj.r.squared, 4),
    statistic = round(statistic, 3),
    p.value = round(p.value, 5)
  ) %>%
  select(r.squared, adj.r.squared, statistic, p.value, df, df.residual)

colnames(model_stats) <- c("R²", "Adjusted R²", "F-statistic", 
                           "p-value (F)", "df Model", "df Residual")

kable(model_stats, 
      caption = "Tabel 4.7: Statistik Kelayakan Model",
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

UJI ASUMSI REGRESI

cat("=== UJI ASUMSI REGRESI LINEAR ===\n\n")

# 1. Normalitas
shapiro_test <- shapiro.test(residuals(model_full))
cat("1. UJI NORMALITAS (Shapiro-Wilk):\n")
cat(sprintf("   W = %.4f, p-value = %.4f\n", 
            shapiro_test$statistic, shapiro_test$p.value))
cat(sprintf("   Interpretasi: %s\n\n", 
            ifelse(shapiro_test$p.value > 0.05, 
                   "✓ Residual normal (p > 0.05)",
                   "✗ Residual tidak normal")))

# 2. Homoskedastisitas
bp_test <- bptest(model_full)
cat("2. UJI HOMOSKEDASTISITAS (Breusch-Pagan):\n")
cat(sprintf("   BP = %.4f, p-value = %.4f\n", 
            bp_test$statistic, bp_test$p.value))
cat(sprintf("   Interpretasi: %s\n\n", 
            ifelse(bp_test$p.value > 0.05, 
                   "✓ Variance homogen (p > 0.05)",
                   "✗ Ada heteroskedastisitas")))

# 3. Multikolinearitas
vif_values <- vif(model_full)
cat("3. UJI MULTIKOLINEARITAS (VIF):\n")
print(round(vif_values, 2))
cat("\n   Interpretasi:\n")
if(all(vif_values < 5)) {
  cat("   ✓ Tidak ada multikolinearitas (semua VIF < 5)\n")
} else {
  cat("   ⚠ Ada multikolinearitas moderat\n")
}

TABEL VIF

vif_df <- data.frame(
  Variabel = names(vif_values),
  VIF = round(vif_values, 2)
)

kable(vif_df, 
      caption = "Tabel 4.8: Nilai VIF (Variance Inflation Factor)",
      align = 'c', row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

RESIDUAL PLOT (GGPLOT)

plot_data <- data.frame(
  fitted = fitted(model_full),
  residuals = residuals(model_full)
)

ggplot(plot_data, aes(x=fitted, y=residuals)) +
  geom_point(alpha=0.6, size=3, color="steelblue") +
  geom_hline(yintercept=0, linetype="dashed", color="red", size=1) +
  geom_smooth(se=FALSE, color="darkblue", method="loess") +
  labs(title="Residuals vs Fitted Values",
       x="Fitted Values (Predicted Rate)", y="Residuals") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5))

Q-Q PLOT

plot_data$standardized <- rstandard(model_full)

ggplot(plot_data, aes(sample=standardized)) +
  stat_qq(color="steelblue", size=3, alpha=0.6) +
  stat_qq_line(color="red", size=1, linetype="dashed") +
  labs(title="Normal Q-Q Plot",
       x="Theoretical Quantiles", y="Standardized Residuals") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5))

27.ACTUAL VS PREDICTED

data_regression$predicted <- predict(model_full)

ggplot(data_regression, aes(x=rate_per_100k, y=predicted)) +
  geom_point(size=3, alpha=0.7, color="steelblue") +
  geom_abline(intercept=0, slope=1, color="red", linetype="dashed", size=1) +
  labs(title="Actual vs Predicted Prevalence Rate",
       subtitle=sprintf("R² = %.3f", summary(model_full)$r.squared),
       x="Actual Rate per 100,000", y="Predicted Rate per 100,000") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(face="bold", hjust=0.5),
        plot.subtitle = element_text(hjust=0.5))

RISK RATIO & ODDS RATIO: KATEGORISASI

data_rr_or <- data_jabar %>%
  mutate(
    prevalensi_tinggi = ifelse(rate_per_100k >= median(rate_per_100k, na.rm=TRUE), 
                                "Ya", "Tidak"),
    kemiskinan_tinggi = ifelse(`penduduk miskin (%)` >= 10, "Ya", "Tidak"),
    sanitasi_buruk = ifelse(`Sanitasi Layak (%)` < 75, "Ya", "Tidak")
  )

cat("Threshold yang digunakan:\n")
cat("- Prevalensi tinggi: ≥", round(median(data_jabar$rate_per_100k, na.rm=TRUE), 2), 
    "per 100k\n")
cat("- Kemiskinan tinggi: ≥10%\n")
cat("- Sanitasi buruk: <75%\n")

TABEL 2x2 KEMISKINAN

tabel_kemiskinan <- table(
  Exposure = data_rr_or$kemiskinan_tinggi,
  Outcome = data_rr_or$prevalensi_tinggi
)

cat("Tabel 2×2: Kemiskinan vs Prevalensi Kusta\n")
addmargins(tabel_kemiskinan)

# Hitung RR dan OR
rr_kemiskinan <- riskratio(tabel_kemiskinan, method="wald")
or_kemiskinan <- oddsratio(tabel_kemiskinan, method="wald")

cat("\nRisk Ratio (RR):", round(rr_kemiskinan$measure[2,1], 2),
    "(95% CI:", round(rr_kemiskinan$measure[2,2], 2), "-",
    round(rr_kemiskinan$measure[2,3], 2), ")\n")
cat("Odds Ratio (OR):", round(or_kemiskinan$measure[2,1], 2),
    "(95% CI:", round(or_kemiskinan$measure[2,2], 2), "-",
    round(or_kemiskinan$measure[2,3], 2), ")\n")

TABEL 2x2 SANITASI

tabel_sanitasi <- table(
  Exposure = data_rr_or$sanitasi_buruk,
  Outcome = data_rr_or$prevalensi_tinggi
)

cat("\nTabel 2×2: Sanitasi Buruk vs Prevalensi Kusta\n")
addmargins(tabel_sanitasi)

# Hitung RR dan OR
rr_sanitasi <- riskratio(tabel_sanitasi, method="wald")
or_sanitasi <- oddsratio(tabel_sanitasi, method="wald")

cat("\nRisk Ratio (RR):", round(rr_sanitasi$measure[2,1], 2),
    "(95% CI:", round(rr_sanitasi$measure[2,2], 2), "-",
    round(rr_sanitasi$measure[2,3], 2), ")\n")
cat("Odds Ratio (OR):", round(or_sanitasi$measure[2,1], 2),
    "(95% CI:", round(or_sanitasi$measure[2,2], 2), "-",
    round(or_sanitasi$measure[2,3], 2), ")\n")

TABEL RINGKASAN RR/OR

ringkasan_rr_or <- data.frame(
  Faktor = c("Kemiskinan Tinggi (≥10%)", "Sanitasi Buruk (<75%)"),
  RR = c(
    sprintf("%.2f", rr_kemiskinan$measure[2,1]),
    sprintf("%.2f", rr_sanitasi$measure[2,1])
  ),
  CI_RR = c(
    sprintf("[%.2f - %.2f]", rr_kemiskinan$measure[2,2], 
            rr_kemiskinan$measure[2,3]),
    sprintf("[%.2f - %.2f]", rr_sanitasi$measure[2,2], 
            rr_sanitasi$measure[2,3])
  ),
  OR = c(
    sprintf("%.2f", or_kemiskinan$measure[2,1]),
    sprintf("%.2f", or_sanitasi$measure[2,1])
  ),
  CI_OR = c(
    sprintf("[%.2f - %.2f]", or_kemiskinan$measure[2,2], 
            or_kemiskinan$measure[2,3]),
    sprintf("[%.2f - %.2f]", or_sanitasi$measure[2,2], 
            or_sanitasi$measure[2,3])
  ),
  Interpretasi = c(
    ifelse(rr_kemiskinan$measure[2,2] > 1, "Faktor Risiko*", 
           ifelse(rr_kemiskinan$measure[2,3] < 1, "Protektif*", "Tidak Signifikan")),
    ifelse(rr_sanitasi$measure[2,2] > 1, "Faktor Risiko*", 
           ifelse(rr_sanitasi$measure[2,3] < 1, "Protektif*", "Tidak Signifikan"))
  )
)

kable(ringkasan_rr_or, 
      caption = "Tabel 4.9: Ringkasan Risk Ratio dan Odds Ratio",
      col.names = c("Faktor Risiko", "RR", "95% CI", "OR", "95% CI", "Interpretasi"),
      align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

```

Lampiran B: Link Dashboard

Dashboard Interaktif: https://monayolaa.shinyapps.io/Dashboard_Kusta_Jabar/

Dashboard menampilkan: - Peta interaktif distribusi kusta - Filter wilayah dan variabel - Tabel data lengkap

Lampiran C: Link Video Presentasi

Video YouTube: https://www.youtube.com/watch?v=-NY1rpUifp8

Durasi: 14 menit Isi: Ringkasan permasalahan, metode, temuan utama, demonstrasi dashboard

Acknowledgement:
Penulis menggunakan ChatGPT (Anthropic Claude) sebagai writing assistant untuk: (1) optimalisasi struktur kalimat akademis, (2) debugging syntax R pada fungsi spatial analysis, dan (3) konsistensi format tabel menggunakan kableExtra package. Metodologi penelitian, pemilihan model regresi, interpretasi uji statistik, dan analisis epidemiologis dilakukan secara mandiri dengan verifikasi terhadap literatur primer.