ANALISIS SPASIAL DAN DETERMINAN EPIDEMIOLOGI
PENYAKIT PNEUMONIA DI PROVINSI DKI JAKARTA
TAHUN 2017–2024

Ditujukan untuk Memenuhi Nilai Mata Kuliah Epidemiologi

Dosen Pengampu:

Dr. I Gede Nyoman Mindra Jaya, S.Si., M.Si.

Disusun Oleh:

Nabiel Alfallah Herdiana (140610230018)

PROGRAM STUDI STATISTIKA
DEPARTEMEN STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR, KABUPATEN SUMEDANG
TAHUN 2025

ABSTRAK

Pneumonia merupakan infeksi saluran pernapasan akut yang menjadi salah satu penyebab utama morbiditas dan mortalitas di Indonesia, khususnya di wilayah urban padat penduduk seperti DKI Jakarta. Penelitian ini bertujuan untuk menganalisis karakteristik epidemiologi Pneumonia secara komprehensif, memetakan distribusi spasial penyakit, serta mengestimasi pengaruh faktor risiko lingkungan (kepadatan penduduk) terhadap kejadian penyakit. Data yang digunakan adalah data sekunder periode 2017–2024 yang mencakup seluruh wilayah administrasi DKI Jakarta. Metode analisis meliputi perhitungan indikator epidemiologi (Incidence, Mortality, CFR, Recovery), pemetaan spasial (Disease Mapping), dan inferensi statistik menggunakan Regresi Data Panel. Hasil penelitian menunjukkan bahwa Jakarta Pusat memiliki beban penyakit tertinggi. Secara temporal, terjadi tren peningkatan kasus yang signifikan di seluruh wilayah, dengan lonjakan drastis yang anomali di Kepulauan Seribu pada tahun 2023–2024. Berdasarkan pemodelan Random Effect, ditemukan bahwa kepadatan penduduk tidak berpengaruh signifikan (p > 0.05) terhadap variasi Incidence Rate, melainkan faktor tren waktu yang menjadi determinan dominan. Temuan ini mengindikasikan bahwa pengendalian penyakit tidak cukup hanya dengan intervensi fisik lingkungan, tetapi memerlukan penguatan sistem surveilans menghadapi tren kenaikan alami.

Kata Kunci: Pneumonia, Epidemiologi Spasial, Regresi Data Panel, Random Effect, DKI Jakarta.

BAB 1. PENDAHULUAN

1.1 Latar Belakang Masalah

Pneumonia masih menjadi penyebab utama kematian balita dan lansia di Indonesia. DKI Jakarta, dengan kepadatan penduduk ekstrem (>16.000 jiwa/km²) dan mobilitas tinggi, menjadi lingkungan berisiko tinggi bagi transmisi penyakit airborne. Kompleksitas pengendalian penyakit ini semakin diperberat oleh polusi udara persisten serta disparitas lingkungan yang mencolok antar-wilayah administrasi (seperti Jakarta Pusat vs Kepulauan Seribu).

Guna mendukung kebijakan kesehatan berbasis bukti (evidence-based policy), analisis data agregat semata tidaklah cukup. Penelitian ini menerapkan pendekatan Epidemiologi Spasial untuk memetakan hotspot risiko secara visual, serta Statistika Inferensia (Regresi Data Panel) untuk mengukur secara presisi pengaruh faktor risiko kepadatan penduduk terhadap laju insidensi penyakit dalam dimensi ruang dan waktu.

1.2 Rumusan Masalah

Berdasarkan latar belakang tersebut, permasalahan utama yang dikaji meliputi: (1) Gambaran epidemiologi (Incidence, Mortality, CFR, Recovery) periode 2017–2024; (2) Pola distribusi spasial dan tren temporal; serta (3) Hubungan asosiasi antara kepadatan penduduk dengan Incidence Rate menggunakan Regresi Data Panel.

1.3 Tujuan Penelitian

Penelitian ini bertujuan mendeskripsikan tren penyakit secara komprehensif, memvisualisasikan wilayah berisiko tinggi (hotspot), dan mengestimasi besaran pengaruh kepadatan penduduk terhadap kejadian Pneumonia menggunakan model statistik yang valid dan robust.

BAB 2. TINJAUAN PUSTAKA

2.1 Konsep Dasar Epidemiologi

Epidemiologi adalah ilmu yang mempelajari distribusi dan determinan keadaan atau kejadian yang berkaitan dengan kesehatan pada populasi tertentu. Kejadian penyakit Pneumonia dapat dijelaskan melalui interaksi Segitiga Epidemiologi (Epidemiologic Triangle), yang meliputi:

  • Agent (Agen): Penyebab utama infeksi, yang paling umum adalah bakteri Streptococcus pneumoniae atau virus Influenza/RSV.
  • Host (Pejamu): Manusia yang rentan. Faktor risiko internal meliputi usia (balita dan lansia lebih rentan), status gizi, dan status imunisasi.
  • Environment (Lingkungan): Faktor eksternal yang mempengaruhi kontak antara agen dan pejamu. Dalam konteks perkotaan, kepadatan penduduk merupakan variabel lingkungan fisik yang krusial karena meningkatkan frekuensi kontak antar individu (contact rate) dan mempermudah transmisi penyakit.

2.2 Ukuran-Ukuran Epidemiologi

Penelitian ini menggunakan empat indikator kunci untuk memotret beban penyakit secara utuh, yaitu:

  1. Incidence Rate (IR): Risiko kasus baru (\(\frac{\text{Kasus Baru}}{\text{Populasi}} \times 100.000\)).
  2. Mortality Rate (MR): Beban kematian (\(\frac{\text{Kematian}}{\text{Populasi}} \times 100.000\)).
  3. Case Fatality Rate (CFR): Tingkat fatalitas (\(\frac{\text{Kematian}}{\text{Kasus Sakit}} \times 100\%\)).
  4. Recovery Rate (RR): Tingkat kesembuhan (\(\frac{\text{Sembuh}}{\text{Kasus Sakit}} \times 100\%\)).

2.3 Analisis Regresi Data Panel

Regresi data panel menggabungkan data cross-section dan time-series, memiliki keunggulan utama dalam mengontrol heterogenitas (karakteristik unik) antar-wilayah yang tidak teramati. Terdapat tiga pendekatan estimasi:

  • Common Effect Model (CEM): Mengabaikan dimensi waktu/ruang dan menggunakan metode OLS biasa.
  • Fixed Effect Model (FEM): Mengasumsikan adanya perbedaan intercept yang bersifat tetap (fixed) untuk menangkap karakteristik spesifik tiap wilayah.
  • Random Effect Model (REM): Mengasumsikan perbedaan antar-wilayah bersifat acak dan masuk ke dalam komponen galat (error term).

Penentuan Model Terbaik (Uji Diagnostik):

Nama Uji Perbandingan Model Keputusan (jika p < 0.05)
Uji Chow CEM vs FEM Pilih FEM
Uji Hausman FEM vs REM Pilih FEM
Uji Lagrange Multiplier CEM vs REM Pilih REM

BAB 3. METODOLOGI

3.1 Sumber Data

Studi ini menggunakan data sekunder agregat periode 2017–2024 yang diperoleh dari sumber resmi berikut:

  1. Seksi Surveilans Dinkes DKI Jakarta: Data epidemiologi utama meliputi jumlah kasus pneumonia dirawat, meninggal (death), dan sembuh (discharge).
  2. BPS Provinsi DKI Jakarta: Data demografi (populasi pertengahan tahun) dan luas wilayah sebagai basis perhitungan Incidence Rate dan kepadatan penduduk.
  3. Data Spasial (Shapefile): Peta digital batas wilayah administrasi untuk keperluan visualisasi pemetaan (choropleth map).

3.2 Desain Penelitian

Studi ini menerapkan Desain Ekologi Longitudinal (Data Panel) dengan unit analisis wilayah (Kabupaten/Kota), bukan individu. Pendekatan ini dipilih untuk menganalisis korelasi faktor lingkungan dan penyakit pada tingkat populasi menggunakan data sekunder agregat.

3.2.1 Variabel dan Sampel

  1. Variabel: Terdiri dari Incidence Rate (Dependen), Kepadatan Penduduk (Independen), dan Tahun (Kontrol).
  2. Populasi & Sampel: Populasi mencakup seluruh penduduk DKI Jakarta (2017–2024). Teknik Total Sampling digunakan terhadap 6 wilayah administrasi (5 Kota, 1 Kabupaten) sehingga data representatif secara provinsi.

3.2.2 Potensi Bias

Mengingat desain ekologi, keterbatasan yang perlu diperhatikan meliputi: Ecological Fallacy: Simpulan level wilayah tidak serta-merta berlaku pada level individu. Under-reporting Bias: Fluktuasi kasus sangat bergantung pada sensitivitas sistem surveilans (artefak data). Omitted Variable Bias: Model belum mengontrol variabel cuaca dan polusi udara spesifik karena keterbatasan data.

3.3 Tahapan Analisis Data

Analisis dilakukan sistematis menggunakan RStudio melalui tahapan:

  1. Pra-pemrosesan: Data cleaning, penggabungan data spasial-demografi-kasus, dan agregasi tahunan.
  2. Epidemiologi Deskriptif: Perhitungan indikator (IR, MR, CFR, RR) serta visualisasi tren temporal dan peta sebaran (disease mapping).
  3. Inferensia (Data Panel): Estimasi model (CEM, FEM, REM), pemilihan model terbaik (Uji Chow, Hausman, LM), dan interpretasi signifikansi statistik.

BAB 4. HASIL DAN PEMBAHASAN

4.1 Gambaran Umum Epidemiologi

Bagian ini menyajikan hasil perhitungan statistik deskriptif untuk empat indikator epidemiologi utama di 6 wilayah administrasi DKI Jakarta selama periode pengamatan (2017–2024).

Tabel 4.1 Profil Epidemiologi Pneumonia DKI Jakarta (Metode Agregat 2017-2024)
Wilayah Administrasi Incidence Rate (Mean) Mortality Rate (Mean) CFR (%) Recovery Rate (%)
JAKARTA UTARA 10.37 0.03 0.33 16.55
JAKARTA PUSAT 9.38 0.05 0.57 17.45
JAKARTA BARAT 9.12 0.04 0.41 7.27
JAKARTA TIMUR 6.16 0.03 0.51 13.73
JAKARTA SELATAN 5.24 0.04 0.72 27.05
KEPULAUAN SERIBU 4.60 0.12 2.54 5.08

Interpretasi Tabel 4.1: Berdasarkan tabel di atas, terlihat adanya disparitas epidemiologi yang signifikan antar wilayah administrasi di DKI Jakarta:

Berdasarkan tabel di atas, terlihat adanya disparitas epidemiologi yang signifikan antar wilayah administrasi di DKI Jakarta:

  1. Episentrum Jakarta Utara: Mencatat Incidence Rate tertinggi (10.40 per 100k). Risiko tinggi ini berkorelasi dengan densitas pemukiman padat pesisir dan aktivitas industri yang mempermudah transmisi airborne.
  2. Anomali Kepulauan Seribu: Menunjukkan pola Low Incidence, High Fatality. Meski kasus terendah, wilayah ini memiliki CFR Tertinggi (2.54%). Hal ini mengindikasikan fenomena gunung es akibat hambatan akses rujukan medis (referral barriers) dan keterlambatan deteksi dini.
  3. Kinerja Kesembuhan: Jakarta Selatan memiliki Recovery Rate terbaik (27.0%). Rendahnya angka kesembuhan agregat DKI (<30%) diduga kuat akibat under-reporting pada pasien rawat jalan yang sembuh mandiri.

4.2 Analisis Tren Temporal dan Spasial

Analisis ini bertujuan untuk melihat dinamika penyakit dari waktu ke waktu dan sebarannya secara geografis.

4.2.1 Tren Waktu (Time Series)

Grafik berikut menampilkan pergerakan Incidence Rate dari tahun 2017 hingga 2024.

Gambar 4.1 Tren Incidence Rate Pneumonia (2017-2024)

Interpretasi Tren: Secara umum, terdapat pola fluktuatif dengan kecenderungan tren kenaikan (uptrend) yang jelas pada tahun-tahun terakhir, khususnya pasca-pandemi COVID-19. Fenomena ini bisa jadi disebabkan oleh “immunity debt” (penurunan kekebalan populasi akibat isolasi panjang) atau perbaikan sistem pelaporan surveilans yang semakin sensitif mendeteksi kasus.

4.2.2 Analisis Spasial (Disease Mapping)

Peta di bawah ini memvisualisasikan konsentrasi kasus di wilayah daratan DKI Jakarta.

Gambar 4.2 Peta Sebaran Incidence Rate (Daratan Jakarta)

Gambar 4.2 Peta Sebaran Incidence Rate (Daratan Jakarta)

Berdasarkan peta kloroplet di atas, zona dengan warna merah pekat (insidensi sangat tinggi) terlihat jelas terkonsentrasi di Jakarta Utara. Hal ini membentuk klaster risiko yang konsisten di wilayah pesisir. Tingginya kasus di Jakarta Utara sangat mungkin berkaitan dengan faktor lingkungan spesifik seperti kepadatan pemukiman kumuh, polusi udara dari aktivitas industri/pelabuhan, serta sanitasi yang kurang memadai di beberapa titik. Wilayah ini menjadi prioritas utama intervensi kesehatan.

4.2.3 Anomali Lonjakan Kasus di Kepulauan Seribu

Mengingat keterbatasan visualisasi peta untuk wilayah kepulauan, berikut disajikan analisis khusus untuk Kepulauan Seribu. Meskipun rata-rata insidensinya rendah (Tabel 4.1), data tahunan menunjukkan fenomena yang mengkhawatirkan.

Tabel 4.2 Tren Incidence Rate di Kepulauan Seribu (per 100k)
Indikator 2017 2018 2019 2020 2021 2022 2023 2024
Incidence Rate 2.09 1.38 4.46 4.5 0.89 2.06 8.47 11.87

Analisis Lonjakan Kasus: Tabel 4.2 memperlihatkan anomali yang signifikan. Setelah bertahun-tahun stabil di angka rendah (0-4 kasus per 100k), terjadi lonjakan ekstrem pada tahun 2023 (8.47) dan 2024 (11.87). Peningkatan insidensi hampir 6 kali lipat dibanding tahun 2022 ini mengindikasikan adanya perubahan pola epidemiologi di wilayah kepulauan. Hal ini bisa disebabkan oleh peningkatan pariwisata pasca-pandemi yang membawa agen penyakit dari daratan, atau adanya outbreak lokal yang belum tertangani maksimal.

4.3 Analisis Asosiasi dan Pemodelan Statistik

Bagian ini menganalisis pengaruh faktor risiko Kepadatan Penduduk terhadap Incidence Rate Pneumonia dengan mengontrol faktor tren waktu, menggunakan pendekatan Regresi Data Panel.

4.3.1 Seleksi Model Terbaik

Penentuan model estimasi dilakukan melalui tiga tahapan uji diagnostik. Berikut ringkasan hasil pengujian:

Tabel 4.3 Ringkasan Hasil Uji Diagnostik Pemilihan Model
Jenis Uji Perbandingan Model Hasil (P-Value) Keputusan
Uji Chow CEM vs FEM 0.0002 (< 0.05) Pilih Fixed Effect
Uji Hausman FEM vs REM 0.9931 (> 0.05) Pilih Random Effect
Uji LM CEM vs REM < 0.05 (Signifikan) Pilih Random Effect

Kesimpulan: Berdasarkan konsistensi hasil uji Hausman dan Lagrange Multiplier, Random Effect Model (REM) ditetapkan sebagai model estimasi terbaik.

4.3.2 Evaluasi Asumsi Klasik

Diagnostik model Random Effect menunjukkan hasil sebagai berikut:

  1. Multikolinearitas: Aman. Nilai VIF (1.001) jauh di bawah 10, mengindikasikan tidak ada korelasi antar prediktor.
  2. Heteroskedastisitas: Terdeteksi (\(p_{BP} = 0.0027\)). Varian residual tidak konstan akibat karakteristik wilayah yang beragam (heterogen).
  3. Autokorelasi: Terindikasi pada taraf 10% (\(p = 0.0666\)), wajar pada data runtun waktu.
  4. Normalitas: Residual tidak berdistribusi normal (\(p < 0.05\)), namun dapat diabaikan berdasar Central Limit Theorem pada jumlah observasi yang memadai.

Tindakan Perbaikan (Robust Correction): Mengingat adanya pelanggaran asumsi heteroskedastisitas dan autokorelasi, estimasi standar error biasa menjadi bias. Oleh karena itu, model akhir dikoreksi menggunakan Robust Standard Errors (HAC Estimator - Metode Arellano) untuk menjamin validitas uji signifikansi (uji-t).

4.3.3 Hasil Estimasi Akhir (Robust Correction)

Berikut adalah hasil estimasi model Random Effect setelah dikoreksi dengan Robust Standard Errors, yang merinci pengaruh Kepadatan Penduduk dan Tren Waktu (Tahun) terhadap Incidence Rate:

Tabel 4.4 Hasil Estimasi Model Random Effect dengan Robust Standard Errors
Var. Prediktor Est. Koefisien Std. Error P-Value Tingkat Signifikansi
(Intercept) 3.3836 1.4446 0.0243 Signifikan
K. Penduduk 0.000174 0.000071 0.0191 Signifikan
Tahun 2018 -0.1229 0.5731 0.8313 Tidak Signifikan
Tahun 2019 0.6679 1.0372 0.5233 Tidak Signifikan
Tahun 2020 0.5485 0.7570 0.4729 Tidak Signifikan
Tahun 2021 -2.5304 0.8514 0.0050 Signifikan
Tahun 2022 -0.7823 0.8850 0.3821 Tidak Signifikan
Tahun 2023 5.4307 2.2018 0.0181 Signifikan
Tahun 2024 8.9353 2.2349 0.0002 Signifikan

Keterangan: Tahun 2017 menjadi kategori referensi (basis). Sumber: Output Analisis R (coeftest with vcovHC).

Interpretasi Hasil Estimasi:

  1. Kepadatan Penduduk (\(p = 0.0191\)): Terbukti berpengaruh positif dan signifikan. Setiap peningkatan densitas populasi meningkatkan risiko transmisi pneumonia. Temuan ini mengonfirmasi bahwa interaksi fisik yang rapat mempermudah penyebaran agen infeksius airborne/droplet.

  2. Dinamika Tren Waktu: Variabel dummy waktu mengungkap pola fluktuasi ekstrem dibanding tahun dasar (2017):

    • Fase Stabil (2018-2020, 2022): Tidak ada perubahan signifikan pada pola insidensi (\(p > 0.05\)).
    • Fase Penurunan (2021): Terjadi penurunan signifikan (Koef: -2.53, \(p < 0.01\)) akibat restriksi mobilitas (PPKM) saat gelombang Delta.
    • Fase Rebound (2023-2024): Terjadi lonjakan tajam pasca-pandemi , memuncak di 2024 (Koef: 8.93, \(p < 0.001\)). Ini mengindikasikan fenomena immunity debt dan dampak normalisasi aktivitas (polusi) yang memperburuk kesehatan paru.

BAB 5. KESIMPULAN DAN SARAN

5.1 Kesimpulan

Berdasarkan analisis spasial, temporal, dan pemodelan statistik Random Effect dengan koreksi Robust terhadap kasus Pneumonia di DKI Jakarta (2017–2024), diperoleh kesimpulan utama sebagai berikut:

  1. Pola Spasial (Hotspot Wilayah): Beban penyakit tertinggi (Incidence Rate) secara rata-rata terkonsentrasi di Jakarta Utara (10.40 per 100k). Hal ini mengindikasikan bahwa wilayah pesisir dengan karakteristik pemukiman padat dan aktivitas pelabuhan/industri memiliki kerentanan lingkungan yang lebih tinggi terhadap penyebaran pneumonia dibandingkan wilayah administratif lainnya.

  2. Disparitas Layanan Kesehatan di Kepulauan Seribu: Ditemukan ketimpangan indikator klinis yang ekstrem di Kepulauan Seribu. Meskipun laju insidensinya relatif rendah dibandingkan daratan, wilayah ini mencatatkan Case Fatality Rate (CFR) Tertinggi (2.54%) dan angka kesembuhan (Recovery Rate) yang tercatat sangat rendah. Tingginya fatalitas ini menjadi indikator kuat adanya hambatan akses rujukan (referral barriers) dan keterlambatan penanganan medis di wilayah kepulauan.

  3. Determinan Penyakit (Hasil Uji Statistik Robust): Berbeda dengan dugaan awal, hasil uji statistik membuktikan bahwa:

    • Kepadatan Penduduk Berpengaruh Signifikan Positif (\(p = 0.0191\)): Wilayah dengan densitas penduduk yang tinggi terbukti secara statistik memiliki laju kasus pneumonia yang lebih tinggi. Hal ini mengonfirmasi mekanisme transmisi droplet yang semakin mudah terjadi pada jarak interaksi fisik yang rapat.
    • Dinamika Tren Waktu (\(p < 0.05\)): Risiko pneumonia tidak statis. Terjadi penurunan signifikan pada tahun 2021 (efek pembatasan sosial pandemi), namun diikuti oleh lonjakan ekstrem (Rebound) pada tahun 2023 dan puncaknya di 2024.

5.2 Saran

Berdasarkan temuan di atas, direkomendasikan langkah-langkah strategis sebagai berikut:

5.2.1 Rekomendasi Kebijakan Kesehatan

  1. Intervensi Berbasis Wilayah:
    • Jakarta Utara & Wilayah Padat: Fokus pada program penyehatan lingkungan pemukiman (ventilasi rumah sehat) dan pengurangan polusi indoor di kecamatan dengan kepadatan tinggi, mengingat variabel kepadatan terbukti signifikan memicu kasus.
    • Kepulauan Seribu: Dinas Kesehatan DKI Jakarta perlu memprioritaskan audit sistem rujukan gawat darurat (evakuasi medis) dan peningkatan fasilitas perawatan intensif di RSUD setempat untuk menekan angka kematian (CFR).
  2. Kewaspadaan Dini Pasca-Pandemi: Mengingat adanya lonjakan kasus yang signifikan pada tahun 2024 (rebound phenomenon), sistem surveilans perlu meningkatkan sensitivitas deteksi dini terhadap klaster pneumonia baru, terutama pada populasi rentan yang mungkin mengalami penurunan imunitas pasca-pandemi (immunity debt).

5.2.2 Rekomendasi Studi Lanjutan

Penelitian ini memiliki keterbatasan dalam ketersediaan data lingkungan spesifik per kecamatan. Peneliti selanjutnya disarankan untuk: 1. Mengintegrasikan data Kualitas Udara (PM 2.5) level mikro dan data Iklim (Kelembaban/Curah Hujan) ke dalam model, karena faktor ini diduga kuat menjadi variabel perancu yang memperbesar efek tren waktu. 2. Melakukan studi kualitatif atau survei primer terkait perilaku pencarian pengobatan (health seeking behavior) masyarakat di Kepulauan Seribu untuk menjelaskan tingginya fatalitas.

DAFTAR PUSTAKA

  1. Badan Pusat Statistik (BPS) Provinsi DKI Jakarta. (2024). Provinsi DKI Jakarta Dalam Angka 2024. Jakarta: BPS Provinsi DKI Jakarta.
  2. Baltagi, B. H. (2013). Econometric Analysis of Panel Data (5th ed.). Chichester: John Wiley & Sons.
  3. Bonita, R., Beaglehole, R., & Kjellstrom, T. (2006). Basic Epidemiology (2nd ed.). Geneva: World Health Organization.
  4. Dinas Kesehatan Provinsi DKI Jakarta. (2023). Profil Kesehatan Provinsi DKI Jakarta Tahun 2022. Jakarta: Dinas Kesehatan Provinsi DKI Jakarta.
  5. Gordis, L. (2014). Epidemiology (5th ed.). Philadelphia: Elsevier Saunders.
  6. Greene, W. H. (2012). Econometric Analysis (7th ed.). Upper Saddle River: Pearson Education.
  7. Kementerian Kesehatan Republik Indonesia. (2022). Pedoman Pengendalian Penyakit Infeksi Saluran Pernapasan Akut. Jakarta: Kemenkes RI.
  8. Lawson, A. B. (2018). Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology (3rd ed.). Boca Raton: CRC Press.
  9. Pemerintah Provinsi DKI Jakarta. (2024). Portal Data Terpadu Pemprov DKI Jakarta (Open Data Jakarta). Diakses dari https://data.jakarta.go.id/.
  10. Riskesdas. (2018). Laporan Nasional Riset Kesehatan Dasar 2018. Jakarta: Badan Penelitian dan Pengembangan Kesehatan.
  11. Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data (2nd ed.). Cambridge: MIT Press.
  12. World Health Organization. (2021). Pneumonia: Key Facts. Diakses dari https://www.who.int/news-room/fact-sheets/detail/pneumonia.

ACKNOWLEDGEMENT

Sebagai bentuk integritas akademik, saya ingin menyatakan secara terbuka mengenai penggunaan Artificial Intelligence (ChatGPT dan Google Gemini) dalam pengerjaan tugas UAS ini.

Dalam proses penyusunan proyek ini, AI digunakan sebagai alat bantu (assistive tool) dengan batasan peran sebagai berikut:

  1. Technical Support: Berperan sebagai debugging partner saat saya menemui kendala teknis (error) pada penulisan kode R, terutama untuk fitur interaktif di R Shiny dan visualisasi peta Leaflet.
  2. Writing Assistant: Membantu memperbaiki struktur kalimat dan tata bahasa agar narasi laporan ini lebih runut dan enak dibaca (human-readable), tanpa mengubah substansi ide asli.
  3. Concept Validator: Berfungsi sebagai mitra diskusi untuk memvalidasi pemahaman saya terkait interpretasi output statistik (khususnya pada model Random Effect dan Robust Standard Errors).

Meskipun menggunakan bantuan alat tersebut, saya menjamin bahwa seluruh verifikasi data, logika analisis, dan pengambilan kesimpulan akhir telah saya lakukan secara mandiri dan sepenuhnya menjadi tanggung jawab saya sebagai penulis.

LAMPIRAN

Lampiran A: Dashboard Analisis Mandiri (R Shiny)

Dashboard interaktif untuk memvisualisasikan data dan peta secara real-time dapat diakses melalui: KLIK DISINI: Link Dashboard Pneumonia DKI Jakarta

Lampiran B: Video Paparan Project

Video presentasi yang menjelaskan metodologi, temuan, dan demo aplikasi dapat dilihat pada tautan berikut: KLIK DISINI: Video Presentasi YouTube

Lampiran C: Source Code Analisis Utama (R)

Berikut adalah cuplikan kode R yang digunakan untuk analisis statistik dan pemodelan:

# ==============================================================================
# SCRIPT ANALISIS UAS EPIDEMIOLOGI
# ==============================================================================

# 1. SETUP
rm(list = ls()) 
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, sf, tmap, plm, zoo, scales, viridis, grid, units)

# --- SESUAIKAN FOLDER ---
setwd("C:/Users/Nabie/OneDrive - Universitas Padjadjaran/Documents/Statistika/Semester 5/Epidem/uas/opsi 4")

# 2. LOAD DATA
tryCatch({
  df_pop_raw       <- read_csv("penduduk_dki_2017_2025.csv", show_col_types = FALSE)
  df_dirawat_raw   <- read_csv("pneumonia_dki_2017_2025_dirawat.csv", show_col_types = FALSE)
  df_meninggal_raw <- read_csv("pneumonia_dki_2017_2025_meninggal.csv", show_col_types = FALSE)
  df_sembuh_raw    <- read_csv("pneumonia_dki_2017_2025_sembuh.csv", show_col_types = FALSE)
  shp_jakarta      <- st_read("Kabupaten-Kota (Provinsi Jakarta).shp", quiet = TRUE) 
}, error = function(e) {
  stop("ERROR: File CSV/SHP tidak ditemukan di folder kerja!")
})

# 3. DATA CLEANING & PREP
# A. Peta & Luas
shp_clean <- shp_jakarta %>%
  mutate(Luas_m2 = st_area(.), Luas_km2 = as.numeric(set_units(Luas_m2, km^2)),
         KUNCI_RAHASIA = case_when(
           str_detect(NAME_2, regex("Seribu", ignore_case = T)) ~ "KEPULAUAN SERIBU",
           str_detect(NAME_2, regex("Pusat", ignore_case = T)) ~ "JAKARTA PUSAT",
           str_detect(NAME_2, regex("Utara", ignore_case = T)) ~ "JAKARTA UTARA",
           str_detect(NAME_2, regex("Barat", ignore_case = T)) ~ "JAKARTA BARAT",
           str_detect(NAME_2, regex("Selatan", ignore_case = T)) ~ "JAKARTA SELATAN",
           str_detect(NAME_2, regex("Timur", ignore_case = T)) ~ "JAKARTA TIMUR",
           TRUE ~ NA_character_
         )) %>% filter(!is.na(KUNCI_RAHASIA)) %>% select(KUNCI_RAHASIA, Luas_km2) %>% st_drop_geometry()

# B. Populasi
pop_clean <- df_pop_raw %>% filter(year <= 2023, area != "DKI Jakarta") %>%
  rename(Wilayah = area, Tahun = year, Populasi = population)
pop_growth <- pop_clean %>% filter(Tahun %in% c(2022, 2023)) %>%
  pivot_wider(names_from = Tahun, values_from = Populasi, names_prefix = "Pop_") %>%
  mutate(Growth_Rate = (Pop_2023 - Pop_2022) / Pop_2022)
pop_2024_proj <- pop_growth %>% mutate(Tahun = 2024, Populasi = round(Pop_2023 * (1 + Growth_Rate), 0)) %>%
  select(Wilayah, Tahun, Populasi)
pop_final <- bind_rows(pop_clean, pop_2024_proj)

# C. Penyakit
clean_disease_data <- function(df) {
  df %>% filter(year < 2025) %>% mutate(across(`1`:`12`, as.numeric)) %>%
    pivot_longer(cols = `1`:`12`, names_to = "Bulan", values_to = "Jumlah") %>%
    mutate(Bulan = as.numeric(Bulan),
           KUNCI_RAHASIA = case_when(
             str_detect(Nama, regex("Pusat", ignore_case=T)) ~ "JAKARTA PUSAT",
             str_detect(Nama, regex("Utara", ignore_case=T)) ~ "JAKARTA UTARA",
             str_detect(Nama, regex("Barat", ignore_case=T)) ~ "JAKARTA BARAT",
             str_detect(Nama, regex("Selatan", ignore_case=T)) ~ "JAKARTA SELATAN",
             str_detect(Nama, regex("Timur", ignore_case=T)) ~ "JAKARTA TIMUR",
             str_detect(Nama, regex("Seribu", ignore_case=T)) ~ "KEPULAUAN SERIBU",
             TRUE ~ "LAINNYA")) %>%
    select(KUNCI_RAHASIA, Wilayah_Asli = Nama, Tahun = year, Bulan, Jumlah)
}
data_dirawat <- clean_disease_data(df_dirawat_raw)
data_meninggal <- clean_disease_data(df_meninggal_raw)
data_sembuh <- clean_disease_data(df_sembuh_raw)

# 4. MERGE & HITUNG RUMUS
data_epi <- data_dirawat %>%
  left_join(data_meninggal, by = c("KUNCI_RAHASIA", "Tahun", "Bulan")) %>%
  left_join(data_sembuh, by = c("KUNCI_RAHASIA", "Tahun", "Bulan")) %>%
  rename(Wilayah = Wilayah_Asli.x, Kasus_Dirawat = Jumlah.x, Kasus_Meninggal = Jumlah.y, Kasus_Sembuh = Jumlah) %>%
  select(Wilayah, KUNCI_RAHASIA, Tahun, Bulan, Kasus_Dirawat, Kasus_Meninggal, Kasus_Sembuh)

data_final <- data_epi %>%
  left_join(pop_final %>% mutate(KUNCI_RAHASIA = case_when(
    str_detect(Wilayah, regex("Seribu", ignore_case=T)) ~ "KEPULAUAN SERIBU",
    str_detect(Wilayah, regex("Pusat", ignore_case=T)) ~ "JAKARTA PUSAT",
    str_detect(Wilayah, regex("Utara", ignore_case=T)) ~ "JAKARTA UTARA",
    str_detect(Wilayah, regex("Barat", ignore_case=T)) ~ "JAKARTA BARAT",
    str_detect(Wilayah, regex("Selatan", ignore_case=T)) ~ "JAKARTA SELATAN",
    str_detect(Wilayah, regex("Timur", ignore_case=T)) ~ "JAKARTA TIMUR",
    TRUE ~ NA_character_)) %>% select(-Wilayah), by = c("KUNCI_RAHASIA", "Tahun")) %>%
  left_join(shp_clean, by = "KUNCI_RAHASIA") %>%
  mutate(
    Tanggal = as.yearmon(paste(Tahun, Bulan, sep = "-")),
    Incidence_Rate = (Kasus_Dirawat / Populasi) * 100000,   
    Mortality_Rate = (Kasus_Meninggal / Populasi) * 100000, 
    CFR = (Kasus_Meninggal / Kasus_Dirawat) * 100,          
    Recovery_Rate = (Kasus_Sembuh / Kasus_Dirawat) * 100,   
    Kepadatan_Penduduk = Populasi / Luas_km2 
  ) %>% filter(!is.na(Populasi)) 

write_csv(data_final, "Hasil_Analisis_Epidemiologi_Pneumonia_Lengkap.csv")

mean(data_final$Recovery_Rate, na.rm=T)

# 5. VISUALISASI GRAFIK TREN
# A. Incidence Rate
plot_ir <- ggplot(data_final, aes(x = as.Date(Tanggal), y = Incidence_Rate, color = Wilayah)) +
  geom_line(linewidth = 1) + geom_point(size = 1.5, alpha = 0.5) +
  labs(title = "Tren Incidence Rate (Risiko Penularan)", x = "Tahun", y = "IR (per 100k)") +
  theme_minimal() + scale_x_date(date_breaks = "1 year", date_labels = "%Y") + theme(legend.position = "bottom")
ggsave("Grafik_Tren_IR.png", plot_ir, width = 10, height = 6)

# B. CFR (Severity)
plot_cfr <- ggplot(data_final, aes(x = as.Date(Tanggal), y = CFR, color = Wilayah)) +
  geom_smooth(se = FALSE, span = 0.3) + geom_point(alpha = 0.3) +
  labs(title = "Tren Case Fatality Rate (Tingkat Fatalitas)", x = "Tahun", y = "CFR (%)") +
  theme_minimal() + theme(legend.position = "bottom")
ggsave("Grafik_Tren_CFR.png", plot_cfr, width = 10, height = 6)

# C. Mortality Rate (Beban Kematian)
plot_mort <- ggplot(data_final, aes(x = as.Date(Tanggal), y = Mortality_Rate, color = Wilayah)) +
  geom_line(linewidth = 1) + geom_point(alpha = 0.5) +
  labs(title = "Tren Mortality Rate (Angka Kematian Penduduk)", x = "Tahun", y = "Mortality Rate (per 100k)") +
  theme_minimal() + theme(legend.position = "bottom")
ggsave("Grafik_Tren_Mortality.png", plot_mort, width = 10, height = 6)

# D. Recovery Rate (Tingkat Kesembuhan)
plot_rec <- ggplot(data_final, aes(x = as.Date(Tanggal), y = Recovery_Rate, color = Wilayah)) +
  geom_smooth(se = FALSE, span = 0.3) + geom_point(alpha = 0.3) +
  labs(title = "Tren Recovery Rate (Tingkat Kesembuhan)", x = "Tahun", y = "Recovery Rate (%)") +
  theme_minimal() + theme(legend.position = "bottom")
ggsave("Grafik_Tren_Recovery.png", plot_rec, width = 10, height = 6)

# 6. PETA SPASIAL (TAHUNAN / FACET MAP)
graphics.off() 
library(sf)
library(tidyverse)
library(tmap)

tmap_mode("plot")

# --- A. PERSIAPAN DATA ---
# 1. Siapkan Shapefile Mentah & Tambah KUNCI_RAHASIA
shp_clean_map <- shp_jakarta %>%
  mutate(KUNCI_RAHASIA = case_when(
    str_detect(NAME_2, regex("Seribu", ignore_case=T)) ~ "KEPULAUAN SERIBU",
    str_detect(NAME_2, regex("Pusat", ignore_case=T)) ~ "JAKARTA PUSAT",
    str_detect(NAME_2, regex("Utara", ignore_case=T)) ~ "JAKARTA UTARA",
    str_detect(NAME_2, regex("Barat", ignore_case=T)) ~ "JAKARTA BARAT",
    str_detect(NAME_2, regex("Selatan", ignore_case=T)) ~ "JAKARTA SELATAN",
    str_detect(NAME_2, regex("Timur", ignore_case=T)) ~ "JAKARTA TIMUR",
    TRUE ~ NA_character_
  )) %>%
  filter(!is.na(KUNCI_RAHASIA))

# 2. Gabungkan dengan Data Angka
data_angka <- data_final %>% 
  st_drop_geometry() %>% 
  select(KUNCI_RAHASIA, Tahun, Incidence_Rate)

shp_lengkap <- shp_clean_map %>%
  left_join(data_angka, by = "KUNCI_RAHASIA") %>%
  filter(!is.na(Tahun)) 

# --- B. FILTER: HANYA DARATAN ---
shp_darat <- shp_lengkap %>% filter(KUNCI_RAHASIA != "KEPULAUAN SERIBU")

# --- C. LABELING ---
shp_darat <- shp_darat %>%
  mutate(
    Label_Text = paste0(KUNCI_RAHASIA, "\n", sprintf("%.1f", Incidence_Rate))
  )

# --- D. VISUALISASI PETA (DARATAN SAJA) ---
min_val <- floor(min(data_final$Incidence_Rate, na.rm=T))
max_val <- ceiling(max(data_final$Incidence_Rate, na.rm=T))
my_breaks <- pretty(c(min_val, max_val), n = 5)

peta_final <- tm_shape(shp_darat) +
  tm_polygons("Incidence_Rate", 
              breaks = my_breaks,
              palette = "Reds", 
              title = "Incidence Rate (per 100k)") +
  tm_text("Label_Text", size = 0.5, 
          col = "black", fontface = "bold", 
          bg.color = "white", bg.alpha = 0.5) +
  tm_facets(by = "Tahun", ncol = 4, free.coords = FALSE) +
  tm_layout(
    main.title = "Evolusi Sebaran Pneumonia DKI Jakarta (2017-2024)",
    main.title.size = 1.2,
    panel.labels = as.character(unique(shp_lengkap$Tahun)),
    panel.label.bg.color = "white",
    frame = FALSE,
    legend.outside = TRUE,
    legend.outside.position = "bottom",
    legend.text.size = 1.1,
    legend.title.size = 1.3
  )

tmap_save(peta_final, "Peta_Sebaran_Tahunan_Lengkap.png", width = 3200, height = 1600, dpi = 200)

# --- E. OUTPUT DATA KEPULAUAN SERIBU (RATA-RATA PER TAHUN) ---

cat("\n=================================================\n")
cat(" RATA-RATA INCIDENCE RATE KEP. SERIBU (PER TAHUN) \n")
cat(" (Silakan salin angka ini ke naskah laporan)\n")
cat("=================================================\n")

# 1. Hitung Rata-rata dulu (Summarise)
stats_kepser_agg <- data_final %>%
  filter(KUNCI_RAHASIA == "KEPULAUAN SERIBU") %>%
  group_by(Tahun) %>%  # KUNCI: Kelompokkan per tahun
  summarise(Mean_IR = mean(Incidence_Rate, na.rm = TRUE)) %>% # Hitung Rata2
  ungroup()

# 2. Print Loop Hasil Agregat
for(i in 1:nrow(stats_kepser_agg)) {
  cat(sprintf("Tahun %d : %.2f\n", stats_kepser_agg$Tahun[i], stats_kepser_agg$Mean_IR[i]))
}

# 3. Hitung Rata-rata Total Seluruh Periode
total_avg <- mean(stats_kepser_agg$Mean_IR, na.rm=T)

cat("-------------------------------------------------\n")
cat("Rata-rata Total (2017-2024): ", sprintf("%.2f", total_avg), "\n")
cat("=================================================\n")

# ==============================================================================
# 7. MODELING & EVALUASI (MODEL SELECTION + ASSUMPTION CHECKS)
# ==============================================================================

if (!require("lmtest")) install.packages("lmtest")
if (!require("car")) install.packages("car")
library(lmtest)
library(car)

# --- A. PERSIAPAN DATA PANEL ---
data_model <- data_final %>%
  group_by(Wilayah, Tahun) %>%
  summarise(
    Incidence_Rate = mean(Incidence_Rate, na.rm = TRUE),
    Kepadatan_Penduduk = mean(Kepadatan_Penduduk, na.rm = TRUE),
    .groups = "drop"
  )

# Definisikan Struktur Panel
pdata <- pdata.frame(
  data_model,
  index = c("Wilayah", "Tahun")
)

# --- B. ESTIMASI 3 JENIS MODEL (UNTUK PERBANDINGAN) ---
# 1. Common Effect (Pooled OLS) - Asumsi: Antar wilayah dianggap sama aja
model_ols <- plm(Incidence_Rate ~ Kepadatan_Penduduk + Tahun, 
                 data = pdata, model = "pooling")

# 2. Fixed Effect (FE) - Asumsi: Tiap wilayah punya karakteristik unik yg tetap
model_fe <- plm(Incidence_Rate ~ Kepadatan_Penduduk + Tahun, 
                data = pdata, model = "within")

# 3. Random Effect (RE) - Asumsi: Perbedaan antar wilayah bersifat acak
model_re <- plm(Incidence_Rate ~ Kepadatan_Penduduk + Tahun, 
                data = pdata, model = "random")

# --- C. UJI PEMILIHAN MODEL (ALASAN LOGIS MEMILIH MODEL) ---
sink("Hasil_Uji_Pemilihan_Model_dan_Asumsi.txt")

cat("=== 1. UJI CHOW (Common Effect vs Fixed Effect) ===\n")
# H0: Common Effect lebih baik (Wilayah dianggap sama)
# H1: Fixed Effect lebih baik (Ada perbedaan karakteristik antar wilayah)
chow_test <- pooltest(model_ols, model_fe)
print(chow_test)
cat("Interp: Jika p-value < 0.05, maka pilih Fixed Effect.\n\n")

cat("=== 2. UJI HAUSMAN (Fixed Effect vs Random Effect) ===\n")
# H0: Random Effect lebih baik (Konsisten & Efisien)
# H1: Fixed Effect lebih baik (Konsisten)
hausman_test <- phtest(model_fe, model_re)
print(hausman_test)
cat("Interp: Jika p-value > 0.05, pilih Random Effect (Aman pakai RE).\n")
cat("        Jika p-value < 0.05, pilih Fixed Effect.\n\n")

cat("=== 3. UJI LM BREUSCH-PAGAN (Common Effect vs Random Effect) ===\n")
# H0: Tidak ada efek panel (Common Effect cukup)
# H1: Ada efek panel (Perlu Random Effect)
lm_test <- plmtest(model_ols, type="bp")
print(lm_test)
cat("Interp: Jika p-value < 0.05, maka Random Effect lebih baik dari OLS biasa.\n\n")

# --- D. UJI ASUMSI KLASIK (DIAGNOSTIK MODEL TERPILIH) ---
model_terpilih <- model_re 

cat("=== 4. UJI NORMALITAS RESIDUAL ===\n")
# Menggunakan Shapiro-Wilk (Jika data < 5000 baris)
shapiro_res <- shapiro.test(residuals(model_terpilih))
print(shapiro_res)
cat("Interp: H0 = Residual Normal. Jika P > 0.05 (Normal). Jika P < 0.05 (Tidak Normal).\n")
cat("Note: Pada data panel besar, asumsi normalitas seringkali dilanggar tapi dimaklumi (CLT).\n\n")

cat("=== 5. UJI MULTIKOLINEARITAS (VIF) ===\n")
# VIF dihitung menggunakan model pooling (OLS) karena struktur formulanya sama
vif_val <- vif(lm(Incidence_Rate ~ Kepadatan_Penduduk + Tahun, data = data_model))
print(vif_val)
cat("Interp: Jika nilai VIF < 10, maka TIDAK terjadi multikolinearitas (Aman).\n\n")

cat("=== 6. UJI HETEROSKEDASTISITAS (Breusch-Pagan) ===\n")
# H0: Varian error konstan (Homoskedastisitas/Bagus)
# H1: Varian error berubah-ubah (Heteroskedastisitas)
bptest_res <- bptest(Incidence_Rate ~ Kepadatan_Penduduk + Tahun, data = data_model, studentize=F)
print(bptest_res)
cat("Interp: Jika P > 0.05 maka Aman (Homoskedastis). Jika P < 0.05 maka Hetero.\n")
cat("Solusi jika Hetero: Gunakan Robust Standard Errors di interpretasi akhir.\n\n")

cat("=== 7. UJI AUTOKORELASI (Serial Correlation) ===\n")
# Wooldridge Test untuk autokorelasi pada data panel
pbg_test <- pbgtest(model_terpilih)
print(pbg_test)
cat("Interp: H0 = Tidak ada autokorelasi. Jika P < 0.05 ada autokorelasi.\n\n")

cat("=== 8. HASIL AKHIR MODEL TERPILIH (SUMMARY) ===\n")
print(summary(model_terpilih))

sink() # Tutup simpan file

# --- E. VISUALISASI EVALUASI (Plot Residual) ---
data_model$Fitted <- as.numeric(fitted(model_terpilih))
data_model$Residuals <- as.numeric(residuals(model_terpilih))

# 1. Histogram Residual (Cek Normalitas Visual)
plot_hist <- ggplot(data_model, aes(x = Residuals)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "skyblue", color = "black") +
  geom_density(color = "red", size = 1) +
  labs(title = "Uji Normalitas: Histogram Residual", x = "Residuals", y = "Density") +
  theme_minimal()
ggsave("Evaluasi_Histogram_Residual.png", plot_hist, width = 8, height = 6)

# 2. Residual vs Fitted (Cek Heteroskedastisitas Visual)
plot_resfit <- ggplot(data_model, aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.5, color = "darkblue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Uji Heteroskedastisitas: Residual vs Fitted",
       subtitle = "Jika menyebar acak tanpa pola = Homoskedastis (Bagus)",
       x = "Fitted Values", y = "Residuals") +
  theme_minimal()
ggsave("Evaluasi_Residual_vs_Fitted.png", plot_resfit, width = 8, height = 6)

print("SUKSES! Semua uji asumsi & pemilihan model selesai.")
print("Buka file 'Hasil_Uji_Pemilihan_Model_dan_Asumsi.txt' untuk melihat angka p-value dan interpretasinya.")

# ==============================================================================
# 8. PERBAIKAN MODEL (ROBUST STANDARD ERRORS)
# ==============================================================================
# Karena uji asumsi menunjukkan adanya Heteroskedastisitas & Autokorelasi,
# maka digunakan Robust Standard Errors (HAC Estimator)
# agar kesimpulan p-value valid.

if (!require("sandwich")) install.packages("sandwich")
library(sandwich)

# A. Tentukan Model Terpilih (Misal: Random Effect dari langkah sebelumnya)
# Pastikan 'model_terpilih' sudah didefinisikan dari langkah 7 (misal: model_re)

# B. Hitung Koefisien dengan Robust Standard Errors (Arellano/White Method)
# vcovHC adalah fungsi untuk menghitung Heteroskedasticity-Consistent Covariance Matrix
robust_model <- coeftest(model_terpilih, vcov = vcovHC(model_terpilih, method = "arellano", type = "HC1"))

# C. Simpan Hasil Robust ke File Teks
sink("Hasil_Akhir_Robust_Correction.txt")

cat("==============================================================\n")
cat(" HASIL REGRESI DENGAN KOREKSI ROBUST (HAC Standard Errors) \n")
cat("==============================================================\n")
cat("Digunakan karena terdeteksi Heteroskedastisitas & Autokorelasi.\n")
cat("Metode Koreksi: Arellano (Panel Data Robust Covariance)\n\n")

print(robust_model)

cat("\n--------------------------------------------------------------\n")
cat("INTERPRETASI:\n")
cat("Lihat kolom 'Pr(>|z|)'. Jika nilai ini < 0.05, maka variabel signifikan berpengaruh\n")
cat("walaupun terdapat masalah asumsi klasik.\n")
cat("==============================================================\n")

sink()

print("SUKSES! Koreksi Robust SE selesai. Cek file 'Hasil_Akhir_Robust_Correction.txt'.")

print("SUKSES! Semua Grafik & Peta Tahunan sudah jadi. Cek folder kerjamu.")

Lampiran D: Source Code Dashboard (Rshiny)

Berikut adalah cuplikan kode R yang digunakan untuk pembuatan dashboard:

# ==============================================================================
# DASHBOARD EPIDEMIOLOGI PNEUMONIA DKI JAKARTA (FINAL FIX COLUMN NAME)
# ==============================================================================

# 1. LOAD LIBRARIES
library(shiny)
library(shinydashboard)
library(tidyverse)
library(leaflet)
library(plotly)
library(DT)
library(sf)
library(scales)
library(shinycssloaders)
library(DiagrammeR)
library(shinyWidgets)
library(shinyjs)

# --- 2. LOAD DATA & HANDLING ---
data_epi <- tryCatch({
  read_csv("Hasil_Analisis_Epidemiologi_Pneumonia_Lengkap.csv", show_col_types = FALSE)
}, error = function(e) { NULL })

shp_jakarta <- tryCatch({
  st_read("Kabupaten-Kota (Provinsi Jakarta).shp", quiet = TRUE)
}, error = function(e) { NULL })

map_ready <- FALSE
if (!is.null(data_epi) && !is.null(shp_jakarta)) {
  shp_clean <- shp_jakarta %>%
    mutate(KUNCI_RAHASIA = case_when(
      str_detect(NAME_2, regex("Seribu", ignore_case = T)) ~ "KEPULAUAN SERIBU",
      str_detect(NAME_2, regex("Pusat", ignore_case = T)) ~ "JAKARTA PUSAT",
      str_detect(NAME_2, regex("Utara", ignore_case = T)) ~ "JAKARTA UTARA",
      str_detect(NAME_2, regex("Barat", ignore_case = T)) ~ "JAKARTA BARAT",
      str_detect(NAME_2, regex("Selatan", ignore_case = T)) ~ "JAKARTA SELATAN",
      str_detect(NAME_2, regex("Timur", ignore_case = T)) ~ "JAKARTA TIMUR",
      TRUE ~ NA_character_
    )) %>%
    filter(!is.na(KUNCI_RAHASIA))
  map_ready <- TRUE
}

# ==============================================================================
# UI (USER INTERFACE)
# ==============================================================================

ui <- dashboardPage(
  title = "Pneumonia DKI Jakarta",
  skin = "black",
  
  # HEADER
  dashboardHeader(
    title = span(tagList(icon("lungs"), "Pneumonia DKI"), 
                 style = "font-weight: 700; font-family: 'Poppins', sans-serif; letter-spacing: 1px;"),
    titleWidth = 260
  ),
  
  # SIDEBAR
  dashboardSidebar(
    width = 260,
    sidebarMenu(
      menuItem("Dashboard Eksekutif", tabName = "dashboard", icon = icon("chart-pie")),
      menuItem("Laporan Epidemiologi", tabName = "epidem_report", icon = icon("file-medical-alt")),
      menuItem("Database", tabName = "data", icon = icon("database")),
      hr(),
      
      # TOMBOL RPUBS (SIDEBAR)
      div(style = "text-align: center; margin-bottom: 15px;",
          actionButton("btn_rpubs_sidebar", 
                       label = span(icon("external-link-alt"), " Buka Laporan (RPubs)", style="color: white; font-weight: bold;"),
                       style = "background-color: #e74c3c; border-color: #c0392b; width: 85%; border-radius: 20px; box-shadow: 0 4px 6px rgba(0,0,0,0.3);",
                       onclick ="window.open('https://rpubs.com/nabielherdiana/pneumonia-dki-jakarta', '_blank')") 
      ),
      hr(),
      
      div(style = "padding: 0 15px 15px 15px;",
          h5("FILTER GLOBAL", style = "color:#7f8c8d; font-weight:700; font-size:11px; margin-bottom:10px;"),
          pickerInput("filter_wilayah", "Wilayah:", choices = c("SEMUA WILAYAH", unique(data_epi$KUNCI_RAHASIA)),
                      selected = "SEMUA WILAYAH", options = list(`style` = "btn-default")),
          br(),
          sliderTextInput("year_range", "Rentang Waktu:", choices = 2017:2024, selected = c(2017, 2024), grid = TRUE)
      )
    )
  ),
  
  # BODY
  dashboardBody(
    tags$head(
      tags$link(rel = "stylesheet", href = "https://fonts.googleapis.com/css2?family=Poppins:wght@300;400;500;600;700&display=swap"),
      tags$style(HTML('
        body, h1, h2, h3, h4, h5, h6, .main-header .logo { font-family: "Poppins", sans-serif !important; }
        .content-wrapper { background-color: #f4f6f9; transition: background-color 0.3s; }
        .box { border-radius: 12px; border: none; box-shadow: 0 3px 12px rgba(0,0,0,0.05); margin-bottom: 20px; transition: background-color 0.3s; }
        .box-header { border-bottom: 1px solid #f4f4f4; padding: 15px; }
        .small-box { border-radius: 12px; transition: transform 0.3s; box-shadow: 0 4px 10px rgba(0,0,0,0.1); }
        .small-box:hover { transform: translateY(-5px); }
        .table-responsive { overflow-x: auto; }
        table { width: 100% !important; font-size: 13px; }
        th { white-space: nowrap; }
        .skin-black .main-header .navbar { background: linear-gradient(to right, #2c3e50, #3498db); }
        .skin-black .main-header .logo { background-color: #2c3e50; color: white; }
        .skin-black .main-sidebar { background-color: #2c3e50; }
      '))
    ),
    
    tabItems(
      
      # --- TAB 1: DASHBOARD ---
      tabItem(tabName = "dashboard",
              fluidRow(column(12, h3("Monitoring Kesehatan Masyarakat: Pneumonia DKI Jakarta", style="font-weight:700; margin-top:0; color:#2c3e50;"))),
              fluidRow(
                valueBoxOutput("vbox_kasus", width = 3),
                valueBoxOutput("vbox_ir", width = 3),
                valueBoxOutput("vbox_mortality", width = 3),
                valueBoxOutput("vbox_cfr", width = 3)
              ),
              fluidRow(
                box(width = 8, title = "Peta Risiko Spasial (Time-Lapse)", status = "danger", solidHeader = TRUE,
                    div(style = "padding: 10px 20px; background-color: #fff3cd; border-radius: 8px; margin-bottom: 15px; color:#856404;",
                        sliderTextInput("map_year", "Geser untuk animasi tahun:", choices = 2017:2024, selected = 2024,
                                        animate = animationOptions(interval = 1500), grid = TRUE, width = "100%")
                    ),
                    withSpinner(leafletOutput("epid_map", height = "450px"), type = 6, color = "#e74c3c")
                ),
                column(width = 4,
                       box(width = 12, title = "Top 3 Wilayah (High Risk)", status = "warning", solidHeader = TRUE,
                           div(style = "overflow-x: auto;", tableOutput("top_rank_table"))
                       ),
                       box(width = 12, title = "Proporsi Kasus", status = "primary", solidHeader = TRUE,
                           plotlyOutput("donut_chart", height = "220px")
                       )
                )
              ),
              fluidRow(
                box(width = 12, title = "Tren Utama: Incidence Rate", status = "primary", solidHeader = TRUE,
                    withSpinner(plotlyOutput("plot_trend_ir_main", height = "350px"), type = 4, color = "#3498db")
                )
              )
      ),
      
      # --- TAB 2: LAPORAN EPIDEMIOLOGI ---
      tabItem(tabName = "epidem_report",
              fluidRow(
                column(12, 
                       div(style="display: flex; justify-content: space-between; align-items: center; border-bottom:2px solid #ddd; padding-bottom:10px; margin-bottom:20px;",
                           h3("Laporan Analisis & Metodologi", style="margin:0; font-weight:700;"),
                           
                           # TOMBOL RPUBS (MAIN)
                           actionButton("btn_rpubs_main", 
                                        label = span(icon("book-open"), " BACA FULL REPORT (RPUBS)", style="color: white; font-weight: bold;"),
                                        class = "btn-danger", 
                                        onclick ="window.open('https://rpubs.com/nabielherdiana/pneumonia-dki-jakarta', '_blank')") 
                       )
                )
              ),
              fluidRow(
                column(6,
                       box(width = 12, title = "1. Segitiga Epidemiologi", status = "info", solidHeader = TRUE,
                           grVizOutput("triad_diagram", height = "400px"), 
                           p("Diagram interaksi dinamis antara Agen, Pejamu, dan Lingkungan.", style="text-align:center; font-size:12px; color:#777; margin-top:10px;")
                       ),
                       box(width = 12, title = "2. Ukuran Epidemiologi", status = "success", solidHeader = TRUE,
                           tags$table(class = "table table-striped",
                                      tags$thead(tags$tr(tags$th("Indikator"), tags$th("Rumus"))),
                                      tags$tbody(
                                        tags$tr(tags$td("Incidence Rate"), tags$td("(Kasus Baru / Populasi) x 100k")),
                                        tags$tr(tags$td("Mortality Rate"), tags$td("(Kematian / Populasi) x 100k")),
                                        tags$tr(tags$td("Case Fatality Rate"), tags$td("(Meninggal / Sakit) x 100%")),
                                        tags$tr(tags$td("Recovery Rate"), tags$td("(Sembuh / Sakit) x 100%"))
                                      )
                           )
                       )
                ),
                column(6,
                       box(width = 12, title = "3. Desain Studi & Metodologi", status = "warning", solidHeader = TRUE,
                           p(strong("Desain: Studi Ekologi (Ecological Time-Series)")),
                           p("Studi ini menganalisis korelasi populasi menggunakan data agregat wilayah.", style="font-size:13px;"),
                           br(),
                           strong("Rincian Variabel:"),
                           tags$ul(style="font-size:13px; padding-left:20px;",
                                   tags$li(strong("Unit Analisis:"), " 6 Wilayah Kota/Kab Administrasi."),
                                   tags$li(strong("Variabel Y:"), " Incidence Rate Pneumonia."),
                                   tags$li(strong("Variabel X:"), " Kepadatan Penduduk (jiwa/km2).")
                           ),
                           br(),
                           strong("Potensi Bias (Keterbatasan):"),
                           tags$ul(style="font-size:13px; padding-left:20px; color:#c0392b;",
                                   tags$li(strong("Ecological Fallacy:"), " Bias kesimpulan individu dari data agregat."),
                                   tags$li(strong("Confounding:"), " Faktor tak teramati (iklim/polusi) tidak masuk model.")
                           )
                       ),
                       box(width = 12, title = "4. Hasil Modeling & Evaluasi", status = "danger", solidHeader = TRUE,
                           p("Model Terpilih: ", strong("Random Effect Model"), " (Hausman p > 0.05) dengan Robust SE."),
                           tableOutput("robust_summary"),
                           hr(),
                           div(style="background-color: #fceceb; padding: 15px; border-radius: 5px; border-left: 5px solid #e74c3c;",
                               h5(icon("lightbulb"), strong("Kesimpulan Analisis (Update):"), style="color:#c0392b; margin-top:0;"),
                               tags$ul(style="padding-left: 20px; color:#333; line-height:1.6; font-size:13px;",
                                       tags$li(strong("Kepadatan Penduduk SIGNIFIKAN (+):"), 
                                               " Terbukti secara statistik (p < 0.05) bahwa wilayah padat memiliki risiko pneumonia lebih tinggi."),
                                       tags$li(strong("Tren Waktu (Rebound):"), 
                                               " Kasus turun signifikan saat pandemi (2021), namun terjadi lonjakan drastis pada 2023-2024 (Rebound Phenomenon).")
                               )
                           )
                       )
                )
              )
      ),
      
      # --- TAB 3: DATABASE ---
      tabItem(tabName = "data",
              fluidRow(
                box(width = 12, title = "Database Lengkap", status = "primary", solidHeader = TRUE,
                    p("Gunakan filter di setiap kolom untuk pencarian spesifik."),
                    dataTableOutput("raw_table")
                )
              )
      )
    )
  )
)

# ==============================================================================
# SERVER (LOGIC)
# ==============================================================================

server <- function(input, output) {
  
  filtered_data_range <- reactive({
    req(data_epi)
    df <- data_epi %>% filter(Tahun >= input$year_range[1] & Tahun <= input$year_range[2])
    if (input$filter_wilayah != "SEMUA WILAYAH") df <- df %>% filter(KUNCI_RAHASIA == input$filter_wilayah)
    df
  })
  
  output$vbox_kasus <- renderValueBox({
    valueBox(comma(sum(filtered_data_range()$Kasus_Dirawat, na.rm=T)), "Total Kasus", icon = icon("hospital-user"), color = "aqua")
  })
  output$vbox_ir <- renderValueBox({
    valueBox(sprintf("%.2f", mean(filtered_data_range()$Incidence_Rate, na.rm=T)), "Avg Incidence Rate", icon = icon("chart-line"), color = "blue")
  })
  output$vbox_mortality <- renderValueBox({
    valueBox(sprintf("%.3f", mean(filtered_data_range()$Mortality_Rate, na.rm=T)), "Avg Mortality Rate", icon = icon("skull"), color = "red")
  })
  output$vbox_cfr <- renderValueBox({
    num <- sum(filtered_data_range()$Kasus_Meninggal, na.rm=T)
    den <- sum(filtered_data_range()$Kasus_Dirawat, na.rm=T)
    val <- ifelse(den > 0, (num/den)*100, 0)
    valueBox(paste0(sprintf("%.2f", val), "%"), "Case Fatality Rate", icon = icon("procedures"), color = "orange")
  })
  
  output$epid_map <- renderLeaflet({
    if (!map_ready) return(leaflet() %>% addTiles())
    data_year <- data_epi %>% filter(Tahun == input$map_year)
    shp_joined <- shp_clean %>% left_join(data_year, by = "KUNCI_RAHASIA")
    pal <- colorNumeric(palette = "Reds", domain = data_epi$Incidence_Rate)
    
    leaflet(shp_joined) %>%
      addProviderTiles(providers$CartoDB.Positron) %>%
      addPolygons(fillColor = ~pal(Incidence_Rate), weight = 2, opacity = 1, color = "white", fillOpacity = 0.8,
                  highlight = highlightOptions(weight = 4, color = "#e74c3c", bringToFront = TRUE),
                  label = ~KUNCI_RAHASIA,
                  popup = ~paste0("<b>", KUNCI_RAHASIA, "</b><br>IR: ", sprintf("%.2f", Incidence_Rate))) %>%
      addLegend("bottomright", pal = pal, values = data_epi$Incidence_Rate, title = "IR")
  })
  
  output$top_rank_table <- renderTable({
    filtered_data_range() %>% 
      group_by(KUNCI_RAHASIA) %>% summarise(IR = mean(Incidence_Rate, na.rm=T)) %>% 
      arrange(desc(IR)) %>% head(3) %>% 
      select(Wilayah = KUNCI_RAHASIA, `Avg IR` = IR)
  }, width = "100%", striped=TRUE)
  
  output$donut_chart <- renderPlotly({
    df <- filtered_data_range() %>% group_by(KUNCI_RAHASIA) %>% summarise(Total = sum(Kasus_Dirawat, na.rm=T))
    plot_ly(df, labels = ~KUNCI_RAHASIA, values = ~Total, type = 'pie', hole = 0.5, textinfo = 'label+percent', showlegend=FALSE) %>% 
      layout(margin = list(t=0,b=0,l=0,r=0), paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)')
  })
  
  output$plot_trend_ir_main <- renderPlotly({
    df_trend <- filtered_data_range() %>% group_by(Tahun, KUNCI_RAHASIA) %>% summarise(Val = mean(Incidence_Rate, na.rm=T), .groups="drop")
    p <- ggplot(df_trend, aes(x = Tahun, y = Val, color = KUNCI_RAHASIA)) +
      geom_line(linewidth = 1) + geom_point(size = 2) +
      theme_minimal() + scale_color_brewer(palette = "Set1") +
      labs(y = "Incidence Rate", x = "Tahun", color = "Wilayah") +
      scale_x_continuous(breaks = seq(min(df_trend$Tahun), max(df_trend$Tahun), 1))
    
    ggplotly(p) %>% config(displayModeBar = FALSE) %>% layout(paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)')
  })
  
  output$triad_diagram <- renderGrViz({
    grViz("
      digraph epidemiological_triad {
        graph [layout = neato, overlap = false, bgcolor = 'transparent']
        node [fontname = 'Poppins', style=filled, penwidth=0]
        
        # AGENT (Merah, Teks Putih)
        AGENT [pos='0,2.2!', shape=circle, width=1.8, fillcolor='#e74c3c', fontcolor='white',
               label='AGENT\n\n- Streptococcus\n- Influenza Virus\n- RSV']
               
        # HOST (Biru, Teks Putih)
        HOST [pos='-2.2,-1.5!', shape=circle, width=1.8, fillcolor='#3498db', fontcolor='white',
              label='HOST\n\n- Balita & Lansia\n- Imunitas Rendah\n- Status Gizi']
              
        # ENVIRONMENT (Hijau, Teks Putih)
        ENV [pos='2.2,-1.5!', shape=circle, width=1.8, fillcolor='#2ecc71', fontcolor='white',
             label='ENVIRONMENT\n\n- Kepadatan Hunian\n- Polusi Udara\n- Ventilasi Buruk']
             
        # PENYAKIT (Kuning, Teks HITAM)
        node [shape=doublecircle, fillcolor='#f1c40f', fontcolor='#2c3e50', width=1.5]
        DISEASE [pos='0,0.2!', label='PNEUMONIA\n(Kejadian Sakit)']
        
        edge [color='#7f8c8d', penwidth=2, arrowhead=normal]
        AGENT -> DISEASE; HOST -> DISEASE; ENV -> DISEASE
        edge [style=dashed, color='#bdc3c7', penwidth=1, arrowhead=none]
        AGENT -> HOST; HOST -> ENV; ENV -> AGENT
      }
    ")
  })
  
  output$robust_summary <- renderTable({
    data.frame(
      Variabel = c("(Intercept)", "Kepadatan Penduduk", 
                   "Tahun 2018", "Tahun 2019", "Tahun 2020", 
                   "Tahun 2021", "Tahun 2022", "Tahun 2023", "Tahun 2024"),
      
      Koefisien = c("3.3836", "0.000174", 
                    "-0.1229", "0.6679", "0.5485", 
                    "-2.5304", "-0.7823", "5.4307", "8.9353"),
      
      `Std. Error` = c("1.4446", "0.000071", 
                       "0.5731", "1.0372", "0.7570", 
                       "0.8514", "0.8850", "2.2018", "2.2349"),
      
      `P-Value` = c("0.024", "0.019", 
                    "0.831", "0.523", "0.473", 
                    "0.005", "0.382", "0.018", "0.000"),
      
      Signifikansi = c("Signifikan", "Signifikan", 
                       "Tidak Signifikan", "Tidak Signifikan", "Tidak Signifikan", 
                       "Signifikan", "Tidak Signifikan", "Signifikan", "Signifikan")
    )
  }, striped = TRUE, width = "100%", align = "l")
  
  # --- DATABASE FIX ---
  output$raw_table <- renderDataTable({
    df <- filtered_data_range()
    
    cols_to_show <- c("KUNCI_RAHASIA", "Tahun", "Bulan", "Incidence_Rate", "Mortality_Rate", "Kasus_Dirawat")
    available_cols <- cols_to_show[cols_to_show %in% names(df)]
    
    bulan_indo <- c("Januari", "Februari", "Maret", "April", "Mei", "Juni", 
                    "Juli", "Agustus", "September", "Oktober", "November", "Desember")
    
    # 1. Pilih kolom dulu
    df_temp <- df %>% 
      select(all_of(available_cols)) %>%
      mutate(across(any_of("Bulan"), ~ {
        if(is.numeric(.) || all(str_detect(as.character(.), "^[0-9]+$"))) {
          idx <- as.numeric(.)
          ifelse(idx >= 1 & idx <= 12, bulan_indo[idx], .)
        } else { . }
      }))
    
    # 2. RENAME KOLOM SECARA PAKSA (Base R - Anti Gagal)
    # Jika kolom KUNCI_RAHASIA ada, ubah namanya jadi Kabupaten/Kota
    if("KUNCI_RAHASIA" %in% names(df_temp)) {
      colnames(df_temp)[colnames(df_temp) == "KUNCI_RAHASIA"] <- "Kabupaten/Kota"
    }
    
    # 3. Rename sisanya (underscore jadi spasi)
    df_final <- df_temp %>%
      rename_with(~ str_replace_all(., "_", " "))
    
    datatable(df_final, 
              filter = 'top', 
              extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                pageLength = 10,
                scrollX = TRUE,
                autoWidth = TRUE
              ),
              rownames = FALSE
    ) %>% 
      formatRound(columns = c("Incidence Rate", "Mortality Rate"), digits = 2)
  })
}

shinyApp(ui, server)