ANALISIS EPIDEMIOLOGI KASUS PNEUMONIA PADA BALITA DI INDONESIA TAHUN 2024
Gusti Ayu Anisa Eldina
2025-12-22
- 1 Pendahuluan
- 2 Tinjauan Pustaka
- 3 Metodologi
- 4 Analisis dan Hasil
- 5 Pembahasan Hasil
Analisis
- 5.1 Gambaran Umum Pneumonia Balita Nasional
- 5.2 Distribusi Insidensi Pneumonia Balita Antarprovinsi
- 5.3 Pola Case Fatality Rate (CFR) Pneumonia Balita
- 5.4 Angka Kematian Pneumonia Balita
- 5.5 Analisis Spasial Global: Moran’s I
- 5.6 Analisis Spasial Lokal (LISA)
- 5.7 Standardized Morbidity Ratio (SMR) dan Empirical Bayes Smoothing
- 5.8 Pemodelan Regresi Poisson dan Negative Binomial
- 5.9 Implikasi Epidemiologis dan Kebijakan
- 5.10 Keterbatasan Analisis
- 6 Kesimpulan dan Saran
- 7 Daftar Pustaka
- 8 Lampiran
Program Studi
S-1 Statistika
NPM
140610230042
Mata Kuliah
Epidemiologi
Dosen Pembimbing Mata Kuliah
I Gede Nyoman Mindra Jaya, P.hd
Abstrak
Pneumonia merupakan salah satu penyebab utama morbiditas dan mortalitas pada balita di Indonesia, dengan distribusi kejadian yang bervariasi antarwilayah. Pemahaman mengenai pola epidemiologis dan spasial pneumonia balita menjadi penting untuk mendukung perencanaan intervensi kesehatan berbasis wilayah. Penelitian ini bertujuan untuk menganalisis distribusi, risiko, dan pola spasial pneumonia balita di Indonesia tahun 2024 serta mengevaluasi keterkaitannya dengan indikator pelayanan kesehatan pada tingkat provinsi. Penelitian ini menggunakan desain studi ekologi dengan unit analisis provinsi. Analisis dilakukan secara deskriptif untuk menggambarkan insidensi, Case Fatality Rate (CFR), dan mortalitas pneumonia balita. Analisis spasial meliputi pemetaan tematik, perhitungan Standardized Morbidity Ratio (SMR) dan Empirical Bayes (EB) smoothing, uji autokorelasi spasial global Moran’s I, serta analisis Local Indicators of Spatial Association (LISA). Selain itu, dilakukan pemodelan regresi Poisson dan Negative Binomial untuk mengevaluasi hubungan antara indikator pelayanan kesehatan dengan kejadian pneumonia balita.
Hasil penelitian menunjukkan adanya variasi insidensi dan risiko pneumonia balita yang cukup besar antarprovinsi. Beberapa provinsi memiliki risiko kejadian yang lebih tinggi dibandingkan rata-rata nasional. Analisis spasial global tidak menunjukkan autokorelasi spasial yang signifikan, namun analisis lokal mengidentifikasi kluster terbatas pada wilayah tertentu. Pemodelan regresi menunjukkan bahwa hubungan antara indikator pelayanan kesehatan dan kejadian pneumonia balita bersifat kompleks dan dipengaruhi oleh heterogenitas data.
Kesimpulannya, pendekatan epidemiologi dan spasial memberikan gambaran komprehensif mengenai pola pneumonia balita di Indonesia dan dapat mendukung pengambilan keputusan kesehatan masyarakat berbasis wilayah.
1 Pendahuluan
1.1 Latar Belakang
Pneumonia merupakan infeksi paru-paru yang menyebabkan akumulasi cairan di alveoli sehingga mengganggu pertukaran udara dan sistem pernapasan. Penyakit ini menjadi salah satu penyebab utama morbiditas dan mortalitas pada anak di bawah usia lima tahun di seluruh dunia, meskipun banyak kasusnya dapat dicegah dan diobati dengan intervensi kesehatan yang tepat (World Health Organization, 2023; UNICEF, 2025). Secara global, pneumonia masih membunuh lebih banyak anak kecil daripada penyakit infeksi lainnya, dengan sekitar lebih dari 700.000 kematian anak balita setiap tahun atau rata-rata satu anak meninggal setiap 43 detik akibat pneumonia (UNICEF, 2025; Kementerian Kesehatan Republik Indonesia, 2024).
Di Indonesia, beban penyakit pneumonia juga tergolong tinggi dan menjadi masalah kesehatan masyarakat yang signifikan, terutama pada kelompok balita yang memiliki sistem kekebalan tubuh yang belum matang. Pemerintah Indonesia mengidentifikasi pneumonia sebagai salah satu penyebab utama kematian balita dan menempatkannya sebagai prioritas intervensi kesehatan publik. Selain memengaruhi angka kematian, pneumonia juga memberikan beban ekonomi yang substansial. Data BPJS Kesehatan mencatat bahwa pneumonia menempati peringkat tertinggi dalam biaya pengobatan penyakit pernapasan pada tahun 2023 (Kementerian Kesehatan Republik Indonesia, 2024).
Selain dampak langsung pada kesehatan, pneumonia memiliki kaitan erat dengan determinan sosial dan lingkungan seperti status gizi, polusi udara, dan akses terhadap layanan kesehatan dasar. Faktor-faktor ini sering kali memperburuk risiko infeksi serta keparahan penyakit, terutama di wilayah dengan infrastruktur dan akses kesehatan yang kurang memadai. Dengan karakteristik geografi, ekonomi, dan pelayanan kesehatan yang bervariasi di tingkat provinsi, analisis yang mempertimbangkan dimensi spasial menjadi sangat relevan untuk memahami distribusi penyakit secara kontekstual dan merancang intervensi yang tepat sasaran.
Relevansi analisis spasial dalam konteks Indonesia diperkuat oleh komitmen terhadap pencapaian target Sustainable Development Goals (SDGs), khususnya SDG 3.2, yang menargetkan pengurangan angka kematian anak balita secara substansial hingga mencapai tingkat terendah yang memungkinkan (United Nations, 2015; World Health Organization, 2023). Analisis epidemiologi berbasis spasial dapat mengidentifikasi wilayah-wilayah dengan risiko tinggi, mengungkap kluster penyakit, serta membantu pembuat kebijakan dalam mengalokasikan sumber daya kesehatan secara efisien untuk mengurangi kesenjangan regional dalam beban pneumonia balita.
1.2 Rumusan Masalah
Keterbatasan informasi mengenai pola spasial pneumonia balita menyulitkan upaya identifikasi wilayah dengan beban penyakit tinggi (high-risk areas) yang memerlukan prioritas intervensi. Tanpa analisis epidemiologi berbasis wilayah, kebijakan kesehatan berisiko bersifat umum dan kurang tepat sasaran, sehingga tidak optimal dalam mendukung pencapaian target penurunan kematian balita sesuai dengan Sustainable Development Goals (SDG 3.2).
Berdasarkan kondisi tersebut, rumusan masalah dalam penelitian ini difokuskan pada pertanyaan utama berikut :
Bagaimana distribusi dan besaran kejadian pneumonia balita antar provinsi di Indonesia pada tahun 2024?
Apakah terdapat pola spasial atau kluster pneumonia balita di Indonesia pada tingkat provinsi
Bagaimana hubungan antara indikator pelayanan kesehatan (cakupan antibiotik dan pengobatan sesuai standar) dengan kejadian dan keparahan pneumonia balita pada tingkat provinsi?
1.3 Tujuan Penelitian
Penelitian ini bertujuan untuk menganalisis distribusi dan pola pneumonia balita di Indonesia tahun 2024 secara epidemiologis dan spasial untuk mendukung pengambilan keputusan kesehatan masyarakat berbasis wilayah. Secara lebih lengkap tujuan penelitian diuraikan menjadi :
Menggambarkan distribusi pneumonia balita antar provinsi di Indonesia secara kuantitatif.
Mengidentifikasi variasi risiko serta kemungkinan kluster pneumonia balita antar provinsi menggunakan pendekatan analisis spasial.
Mengevaluasi keterkaitan indikator pelayanan kesehatan, khususnya pemberian antibiotik dan pengobatan sesuai standar, dengan beban kasus dan keparahan pneumonia balita pada tingkat provinsi.
Menghasilkan informasi berbasis data yang dapat mendukung perencanaan dan pengambilan keputusan kesehatan oleh masyarakat dan pemerintah.
2 Tinjauan Pustaka
2.1 Konsep Epidemiologi Pneumonia Balita
2.1.1 Agent–Host–Environment
Pneumonia pada balita merupakan penyakit infeksi saluran pernapasan bawah yang bersifat multifaktorial dan dipengaruhi oleh interaksi antara agen penyebab (agent), faktor individu (host), dan lingkungan (environment). Kerangka agent–host–environment merupakan konsep dasar epidemiologi yang digunakan untuk memahami proses terjadinya penyakit serta mengidentifikasi titik intervensi yang efektif dalam pencegahan dan pengendalian penyakit.
Agent (Agen Penyebab)
Agen penyebab pneumonia balita umumnya berupa mikroorganisme infeksius, baik virus maupun bakteri. Secara global dan nasional, agen yang paling sering terlibat antara lain Streptococcus pneumoniae, Haemophilus influenzae, Respiratory Syncytial Virus (RSV), serta virus influenza. Agen-agen ini dapat menyebabkan peradangan jaringan paru sehingga mengganggu pertukaran oksigen dan menimbulkan gejala klinis pneumonia (World Health Organization \[WHO\], 2023).
Di Indonesia, pneumonia balita sebagian besar ditegakkan berdasarkan diagnosis klinis di fasilitas pelayanan kesehatan primer, sehingga spektrum agen yang terlibat dapat bervariasi antar wilayah. Keterlambatan diagnosis atau tatalaksana yang tidak sesuai standar dapat meningkatkan risiko pneumonia berat dan kematian, terutama pada balita dengan kondisi rentan (Kementerian Kesehatan Republik Indonesia, 2024).
Host (Faktor Individu Balita)
Faktor host memainkan peran penting dalam menentukan kerentanan terhadap pneumonia serta keparahan penyakit. Balita, khususnya bayi usia kurang dari satu tahun, memiliki sistem imun yang belum berkembang sempurna sehingga lebih rentan terhadap infeksi saluran pernapasan. Selain usia, status gizi yang buruk, berat badan lahir rendah, riwayat prematuritas, serta komorbiditas tertentu dapat meningkatkan risiko terjadinya pneumonia dan komplikasinya (UNICEF, 2025).
Di Indonesia, faktor host juga berkaitan dengan praktik pengasuhan, seperti keterlambatan mencari pertolongan medis dan rendahnya pemanfaatan layanan kesehatan. Perbedaan karakteristik host ini berkontribusi terhadap variasi kejadian dan keparahan pneumonia balita antar wilayah, sebagaimana tercermin dalam perbedaan angka kasus dan case fatality rate antar provinsi.
Environment (Faktor Lingkungan dan Sistem Pelayanan)
Lingkungan fisik dan sosial berperan besar dalam penularan dan keparahan pneumonia balita. Faktor lingkungan yang relevan meliputi kepadatan hunian, ventilasi rumah yang buruk, paparan asap rokok dan polusi udara, serta kondisi sosial ekonomi yang rendah. Lingkungan dengan kualitas udara yang buruk dapat meningkatkan kerentanan saluran pernapasan balita terhadap infeksi (WHO, 2023).
Selain lingkungan fisik, sistem pelayanan kesehatan juga merupakan bagian dari komponen lingkungan dalam konteks kesehatan masyarakat. Akses terhadap fasilitas kesehatan, ketersediaan tenaga medis, serta kualitas tatalaksana klinis—termasuk pemberian antibiotik dan pengobatan sesuai standar—sangat memengaruhi luaran pneumonia balita. Di Indonesia, variasi akses dan kualitas layanan kesehatan antar provinsi dapat menyebabkan perbedaan dalam deteksi dini, penanganan kasus, dan risiko kematian akibat pneumonia balita (Kementerian Kesehatan Republik Indonesia, 2024).
Sintesis Kerangka Agent–Host–Environment
Interaksi antara agent, host, dan environment membentuk dasar terjadinya pneumonia balita. Agen infeksius yang beredar di masyarakat dapat menyebabkan penyakit terutama pada balita dengan faktor host yang rentan, dan risiko tersebut diperberat oleh kondisi lingkungan yang tidak mendukung kesehatan serta keterbatasan sistem pelayanan. Oleh karena itu, pendekatan epidemiologi yang mempertimbangkan ketiga komponen ini sangat penting untuk memahami distribusi pneumonia balita di Indonesia dan merancang intervensi yang komprehensif serta berbasis wilayah.
2.1.2 Ukuran Epidemiologi
Ukuran epidemiologi digunakan untuk menggambarkan besarnya masalah kesehatan, membandingkan risiko antar kelompok atau wilayah, serta menilai dampak faktor risiko terhadap kejadian penyakit. Dalam konteks pneumonia balita di Indonesia, ukuran epidemiologi yang relevan meliputi ukuran frekuensi, ukuran asosiasi, dan ukuran dampak yang disesuaikan dengan karakteristik data agregat wilayah.
Ukuran Frekuensi
Ukuran frekuensi menggambarkan besarnya kejadian pneumonia
balita dalam suatu populasi dan periode tertentu. Pada penelitian ini,
data yang tersedia berupa jumlah kasus pneumonia balita selama satu
tahun (2024) dan jumlah balita per provinsi, sehingga ukuran frekuensi
yang paling tepat adalah insidensi kumulatif.
Insidensi kumulatif (Cumulative Incidence)
Insidensi kumulatif didefinisikan sebagai proporsi individu yang
mengalami penyakit dalam suatu periode waktu tertentu. Dalam studi ini,
insidensi kumulatif pneumonia balita dihitung sebagai jumlah kasus
pneumonia balita dibagi dengan jumlah balita pada masing-masing provinsi
dan dinyatakan per 1.000 balita. Ukuran ini menggambarkan risiko
rata-rata seorang balita untuk mengalami pneumonia selama tahun
2024.
Selain insidensi, ukuran frekuensi lain yang relevan adalah Case Fatality Rate (CFR). CFR menunjukkan proporsi kematian di antara kasus pneumonia yang terlaporkan dan digunakan sebagai indikator keparahan penyakit serta kualitas tatalaksana klinis. CFR dihitung sebagai jumlah kematian akibat pneumonia dibagi dengan jumlah kasus pneumonia balita dalam periode yang sama dan dinyatakan dalam persentase (Gordis, 2014).
Ukuran seperti prevalensi tidak dapat dihitung secara valid dalam penelitian ini karena data yang tersedia tidak membedakan antara kasus baru dan kasus yang sedang sakit pada suatu titik waktu tertentu. Oleh karena itu, analisis difokuskan pada insidensi dan CFR yang sesuai dengan struktur data tahunan.
Ukuran Asosiasi
Ukuran asosiasi digunakan untuk membandingkan risiko atau laju
kejadian penyakit antara kelompok yang berbeda. Dalam studi epidemiologi
berbasis wilayah dengan data agregat, ukuran asosiasi yang paling sesuai
adalah Risk Ratio (RR) atau Incidence Rate
Ratio (IRR) yang diperoleh dari model regresi Poisson atau
Negative Binomial.
Risk Ratio (RR / IRR)
Risk ratio menggambarkan perbandingan risiko atau laju kejadian
pneumonia balita antara dua kelompok wilayah, misalnya provinsi dengan
cakupan pengobatan sesuai standar yang tinggi dibandingkan dengan
provinsi dengan cakupan yang lebih rendah. Nilai RR lebih besar dari
satu menunjukkan risiko yang lebih tinggi pada kelompok terpapar,
sedangkan nilai kurang dari satu menunjukkan risiko yang lebih rendah
(Rothman, Greenland, & Lash, 2008).
Dalam penelitian ini, RR diestimasi menggunakan model regresi Poisson dengan offset jumlah balita untuk memperhitungkan perbedaan ukuran populasi antar provinsi. Pendekatan ini memungkinkan perbandingan laju kasus pneumonia secara adil antar wilayah.
Penggunaan odds ratio (OR) tidak direkomendasikan dalam studi ini karena OR umumnya digunakan pada desain case-control atau data individu, dan dapat menimbulkan interpretasi yang menyesatkan bila diterapkan pada data agregat wilayah.
Ukuran Dampak (Attributable Measures)
Selain ukuran asosiasi, epidemiologi juga menggunakan ukuran
dampak untuk mengestimasi kontribusi faktor risiko terhadap kejadian
penyakit.
Attributable Risk (AR)
Attributable risk menunjukkan selisih risiko penyakit antara kelompok
terpapar dan tidak terpapar. Ukuran ini menggambarkan besarnya risiko
pneumonia balita yang dapat dikaitkan dengan suatu faktor pada kelompok
terpapar.
Attributable Fraction (AF)
Attributable fraction menyatakan proporsi kejadian penyakit pada
kelompok terpapar yang dapat dicegah apabila paparan tersebut
dihilangkan, dengan asumsi hubungan kausal. Dalam konteks analisis
berbasis RR, AF dapat dihitung dari nilai RR yang diperoleh dari model
regresi.
Pada tingkat populasi, Population Attributable Fraction (PAF) dapat digunakan untuk memperkirakan proporsi kasus pneumonia balita di seluruh provinsi yang berpotensi dicegah jika faktor risiko tertentu dapat dikendalikan. Namun, interpretasi ukuran dampak dalam studi ekologi harus dilakukan secara hati-hati karena keterbatasan inferensi kausal pada data agregat.
Relevansi Ukuran Epidemiologi dalam Analisis Wilayah
Penggunaan ukuran frekuensi, asosiasi, dan dampak secara
bersamaan memberikan gambaran komprehensif mengenai pneumonia balita di
Indonesia. Insidensi dan CFR menggambarkan besarnya masalah dan tingkat
keparahan penyakit, sementara RR dan ukuran dampak membantu memahami
perbedaan risiko antar wilayah serta potensi perbaikan melalui
intervensi layanan kesehatan. Pendekatan ini mendukung analisis
epidemiologi berbasis wilayah yang sejalan dengan tujuan pengendalian
pneumonia balita dan pencapaian target SDG 3.2.
2.1.3 Desain Studi Epidemiologi
Desain studi epidemiologi merupakan kerangka metodologis yang digunakan untuk mengkaji hubungan antara paparan dan kejadian penyakit. Pemilihan desain studi harus disesuaikan dengan tujuan penelitian, karakteristik data, serta sumber daya yang tersedia. Dalam analisis pneumonia balita di Indonesia tahun 2024, desain yang relevan adalah studi potong lintang (cross-sectional) dengan pendekatan studi ekologi.
Cross-sectional (Studi Potong Lintang)
Cross-sectional merupakan desain observasional yang mengukur paparan dan outcome pada waktu yang sama dalam suatu populasi. Desain ini banyak digunakan untuk menggambarkan distribusi penyakit dan mengidentifikasi hubungan awal antara faktor risiko dan outcome (Gordis, 2014). Pada konteks pneumonia balita, studi potong lintang memungkinkan peneliti untuk menilai besarnya kejadian pneumonia dan indikator pelayanan kesehatan dalam satu periode waktu tertentu, yaitu tahun 2024.
Kelebihan Cross-sectional :
Efisien dari segi waktu dan biaya karena pengumpulan data dilakukan pada satu periode.
Cocok untuk analisis deskriptif dan eksploratif, termasuk perhitungan insidensi kumulatif dan CFR.
Dapat digunakan untuk perbandingan antar wilayah dalam skala nasional.
Keterbatasan studi potong lintang:
Tidak dapat memastikan urutan sebab-akibat antara paparan dan penyakit (temporal ambiguity).
Rentan terhadap reverse causality, di mana outcome dapat memengaruhi paparan.
Kurang sesuai untuk mengkaji penyakit dengan durasi sangat singkat atau fluktuasi musiman tanpa data tambahan.
Studi Ekologi (Ecological Study)
Studi ekologi merupakan desain epidemiologi yang menggunakan
data agregat kelompok atau wilayah sebagai unit analisis, bukan
individu. Dalam penelitian ini, unit analisis adalah provinsi di
Indonesia, dengan variabel yang direkap pada tingkat provinsi. Studi
ekologi sangat relevan untuk analisis berbasis wilayah dan pemetaan
penyakit, terutama ketika data individu tidak tersedia (Rothman et al.,
2008).
Kelebihan studi ekologi :
Memungkinkan analisis distribusi penyakit pada skala geografis yang luas.
Mendukung identifikasi variasi risiko dan kluster penyakit antar wilayah.
Sangat berguna untuk perencanaan kebijakan dan alokasi sumber daya kesehatan masyarakat.
Memungkinkan integrasi analisis spasial dan pemodelan penyakit berbasis wilayah.
Keterbatasan studi ekologi :
Rentan terhadap ecological fallacy, yaitu kesalahan menarik kesimpulan pada tingkat individu berdasarkan data kelompok.
Variabel perancu pada tingkat individu tidak dapat dikendalikan secara langsung.
Kualitas hasil sangat bergantung pada kelengkapan dan konsistensi data agregat.
Relevansi Desain Studi dengan Analisis Pneumonia Balita
Kombinasi desain potong lintang dan pendekatan ekologi sesuai
untuk tujuan penelitian ini, yaitu mendeskripsikan distribusi pneumonia
balita, mengidentifikasi pola spasial, serta mengeksplorasi hubungan
indikator pelayanan kesehatan dengan beban penyakit pada tingkat
provinsi. Meskipun desain ini memiliki keterbatasan dalam inferensi
kausal, hasilnya tetap bernilai penting sebagai dasar perencanaan dan
evaluasi program kesehatan masyarakat, khususnya dalam konteks
pencapaian target penurunan kematian balita sesuai SDG 3.2.
2.2 Kerangka Konseptual
Kerangka konseptual dalam penelitian ini menggambarkan bahwa kejadian pneumonia pada balita merupakan hasil interaksi antara faktor agen, host, lingkungan, serta sistem pelayanan kesehatan. Agen penyebab berupa virus dan bakteri pernapasan dapat menginfeksi balita, terutama ketika didukung oleh kondisi lingkungan yang tidak sehat, seperti kepadatan hunian tinggi, ventilasi rumah yang buruk, paparan asap rokok atau polusi, serta kondisi sosial ekonomi dan sanitasi yang rendah. Di sisi lain, faktor host seperti usia balita yang masih sangat muda, imunitas yang belum optimal, status gizi buruk, serta adanya penyakit penyerta meningkatkan kerentanan balita untuk mengalami pneumonia setelah terpapar agen infeksi.
Setelah pneumonia terjadi, luaran penyakit sangat dipengaruhi oleh sistem pelayanan kesehatan. Keterbatasan akses layanan dan keterlambatan deteksi dini dapat meningkatkan risiko pneumonia berkembang menjadi pneumonia berat dan berujung pada kematian. Sebaliknya, intervensi pelayanan kesehatan yang baik, seperti pemberian antibiotik yang tepat dan pengobatan sesuai standar, berperan penting dalam menurunkan keparahan penyakit dan risiko kematian (CFR). Dengan demikian, kerangka ini menunjukkan bahwa upaya pengendalian pneumonia balita tidak hanya berfokus pada pengobatan, tetapi juga memerlukan perbaikan kondisi lingkungan, penguatan faktor host, serta peningkatan akses dan kualitas layanan kesehatan, terutama pada wilayah dengan beban penyakit yang tinggi.
3 Metodologi
3.1 Sumber Data dan Pemilihan Kasus
Penelitian ini menggunakan data sekunder yang bersumber dari Profil Kesehatan Indonesia Tahun 2024 yang diterbitkan oleh Kementerian Kesehatan Republik Indonesia. Secara khusus, data yang digunakan berasal dari Lampiran 58: Penemuan Kasus Pneumonia Balita Menurut Kelompok Umur, Jenis Kelamin, dan Provinsi Tahun 2024, yang memuat rekapitulasi jumlah kasus pneumonia balita pada tingkat provinsi di seluruh Indonesia. Dataset ini mencakup informasi jumlah balita, perkiraan kasus pneumonia, jumlah kasus pneumonia dan pneumonia berat menurut kelompok umur dan jenis kelamin, jumlah kematian akibat pneumonia, serta indikator pelayanan kesehatan seperti jumlah kasus yang mendapatkan antibiotik dan persentase pengobatan sesuai standar Profil Kesehatan Indonesia 2024.
Pemilihan data pneumonia balita didasarkan pada pertimbangan bahwa pneumonia merupakan salah satu penyebab utama morbiditas dan mortalitas pada anak di bawah usia lima tahun di Indonesia serta menjadi indikator penting dalam pemantauan pencapaian target Sustainable Development Goals (SDG) 3.2, yaitu penurunan angka kematian balita. Data dari Profil Kesehatan Indonesia dipilih karena bersifat resmi, nasional, terstandar, dan dikompilasi secara rutin melalui sistem pelaporan kesehatan, sehingga memungkinkan perbandingan antar provinsi secara konsisten.
Kasus yang dianalisis dalam penelitian ini adalah seluruh kasus pneumonia balita (usia <5 tahun) yang tercatat pada tahun 2024, tanpa pembatasan tambahan berdasarkan fasilitas kesehatan atau karakteristik individu. Unit analisis yang digunakan adalah provinsi, sehingga penelitian ini menggunakan pendekatan studi ekologi dengan data agregat wilayah. Pendekatan ini dipilih karena sesuai dengan tujuan penelitian untuk menggambarkan distribusi pneumonia balita secara spasial, mengidentifikasi variasi risiko antar provinsi, serta mendukung analisis pemetaan penyakit dan pemodelan epidemiologi berbasis wilayah.
3.2 Variabel Penelitian
3.2.1 Variabel Outcome
Variabel outcome dalam penelitian ini merepresentasikan besaran kejadian dan keparahan pneumonia balita pada tingkat provinsi. Variabel-variabel tersebut digunakan untuk menggambarkan beban penyakit serta sebagai keluaran utama dalam analisis epidemiologi dan pemodelan.
1. Jumlah kasus pneumonia balita
(case_tot)
Jumlah kasus pneumonia balita merupakan total kasus pneumonia yang
tercatat pada balita usia <5 tahun (laki-laki dan perempuan) di
setiap provinsi selama tahun 2024. Variabel ini menjadi outcome
utama dalam penelitian dan digunakan sebagai dasar perhitungan
ukuran frekuensi (insidensi kumulatif) serta sebagai variabel respon
dalam model regresi Poisson atau Negative Binomial. Karena berupa data
hitungan (count data), variabel ini dianalisis dengan
mempertimbangkan jumlah balita sebagai offset.
2. Jumlah kematian akibat pneumonia balita
(death_tot)
Jumlah kematian akibat pneumonia balita adalah total kematian balita
yang disebabkan oleh pneumonia di setiap provinsi selama tahun 2024.
Variabel ini digunakan untuk menilai dampak fatal pneumonia balita serta
sebagai komponen dalam perhitungan case fatality rate. Selain
itu, variabel ini juga digunakan untuk menghitung mortality
rate spesifik pneumonia balita pada tingkat provinsi.
3. Case Fatality Rate (CFR)
Case Fatality Rate (CFR) merupakan proporsi kematian di antara
seluruh kasus pneumonia balita yang terlaporkan dalam periode tertentu.
CFR dihitung sebagai jumlah kematian akibat pneumonia balita dibagi
dengan jumlah kasus pneumonia balita, kemudian dinyatakan dalam
persentase. CFR digunakan sebagai indikator keparahan penyakit dan
kualitas tatalaksana klinis, karena nilai CFR yang lebih tinggi dapat
mencerminkan keterlambatan diagnosis, keterbatasan akses layanan, atau
penanganan yang kurang optimal.
4. Ukuran turunan outcome
Selain variabel outcome utama, penelitian ini juga menggunakan beberapa
ukuran turunan untuk memperkaya analisis, antara lain :
Insidensi kumulatif pneumonia balita, yaitu jumlah kasus pneumonia balita per 1.000 balita di setiap provinsi.
Mortality rate pneumonia balita, yaitu jumlah kematian akibat pneumonia balita per 100.000 balita.
3.2.2 Variabel Paparan dan Pendukung
Variabel paparan dan pendukung dalam penelitian ini digunakan untuk mengeksplorasi keterkaitan antara indikator pelayanan kesehatan serta karakteristik populasi dengan kejadian dan keparahan pneumonia balita pada tingkat provinsi. Variabel-variabel ini berperan sebagai explanatory variables dalam analisis epidemiologi dan pemodelan.
1. Persentase pengobatan sesuai standar
(tx_std_pct)
Variabel tx_std_pct merepresentasikan persentase kasus
pneumonia balita yang mendapatkan pengobatan sesuai dengan standar
pelayanan kesehatan yang ditetapkan. Variabel ini digunakan sebagai
indikator kualitas tatalaksana klinis pada tingkat provinsi. Cakupan
pengobatan sesuai standar yang lebih tinggi diharapkan berkaitan dengan
penurunan keparahan penyakit dan risiko kematian akibat pneumonia
balita. Dalam analisis, variabel ini digunakan sebagai variabel paparan
utama untuk mengevaluasi hubungan kualitas pelayanan kesehatan dengan
luaran penyakit.
2. Jumlah kasus pneumonia yang mendapatkan antibiotik
(abx_n)
Variabel abx_n menunjukkan jumlah absolut kasus pneumonia
balita yang menerima terapi antibiotik di setiap provinsi selama tahun
2024. Variabel ini mencerminkan cakupan intervensi pengobatan terhadap
kasus pneumonia yang terlaporkan. Namun, karena nilai abx_n
secara langsung dipengaruhi oleh jumlah kasus pneumonia, interpretasinya
dilakukan secara hati-hati, terutama dalam analisis asosiasi.
3. Cakupan pemberian antibiotik
(abx_cov)
Untuk meningkatkan interpretabilitas analisis, variabel
abx_cov digunakan sebagai ukuran turunan yang dihitung
sebagai persentase kasus pneumonia balita yang mendapatkan antibiotik
terhadap total kasus pneumonia balita. Variabel ini memberikan gambaran
proporsional mengenai cakupan terapi antibiotik dan digunakan sebagai
variabel paparan alternatif dalam pemodelan, sehingga mengurangi potensi
bias akibat hubungan mekanis dengan jumlah kasus.
4. Jumlah balita (balita_n)
Variabel balita_n merepresentasikan jumlah balita usia
<5 tahun di setiap provinsi dan digunakan sebagai
denominator dalam perhitungan ukuran frekuensi, seperti
insidensi kumulatif dan mortality rate. Dalam pemodelan regresi
Poisson atau Negative Binomial, variabel ini digunakan sebagai
offset (dalam bentuk logaritmik) untuk memperhitungkan
perbedaan ukuran populasi antar provinsi.
5. Variabel pendukung lainnya
Selain variabel utama di atas, penelitian ini juga mempertimbangkan
variabel pendukung berupa ukuran turunan seperti insidensi kumulatif dan
CFR untuk keperluan deskriptif, visualisasi, serta interpretasi
epidemiologis. Variabel-variabel pendukung ini tidak digunakan sebagai
paparan langsung, tetapi berperan dalam memberikan konteks terhadap
beban penyakit dan variasi antar wilayah.
3.2.3 Populasi dan Unit Analisis
Populasi target dalam penelitian ini adalah seluruh balita berusia kurang dari lima tahun (<5 tahun) yang tinggal di Indonesia pada tahun 2024. Populasi ini dipilih karena balita merupakan kelompok usia yang paling rentan terhadap pneumonia akibat sistem kekebalan tubuh yang belum berkembang secara optimal serta tingginya paparan faktor risiko lingkungan dan sosial.
Populasi terobservasi dalam penelitian ini adalah seluruh balita yang tercatat dalam sistem pelaporan kesehatan dan direkap dalam Profil Kesehatan Indonesia Tahun 2024. Data yang digunakan mencerminkan hasil kompilasi laporan rutin pelayanan kesehatan yang dikumpulkan oleh Kementerian Kesehatan Republik Indonesia pada tingkat nasional.
Unit analisis yang digunakan adalah provinsi di Indonesia. Setiap provinsi diperlakukan sebagai satu unit observasi dengan variabel-variabel yang direpresentasikan dalam bentuk data agregat wilayah, seperti jumlah balita, jumlah kasus pneumonia balita, jumlah kematian akibat pneumonia, serta indikator pelayanan kesehatan. Pemilihan unit analisis provinsi memungkinkan analisis distribusi penyakit secara spasial, perbandingan risiko antar wilayah, serta penerapan pemodelan epidemiologi dalam konteks pemetaan penyakit.
Penggunaan populasi dan unit analisis ini sesuai dengan pendekatan studi ekologi cross-sectional, di mana interpretasi hasil difokuskan pada perbedaan antar wilayah, bukan pada tingkat individu. Oleh karena itu, temuan penelitian ini dimaksudkan untuk mendukung perencanaan dan evaluasi kebijakan kesehatan masyarakat berbasis wilayah, dengan tetap mempertimbangkan keterbatasan inferensi kausal pada tingkat individu.
3.3 Metode Analisis
Analisis data dalam penelitian ini dilakukan secara bertahap untuk menggambarkan besarnya beban pneumonia balita, mengidentifikasi pola spasial antar provinsi, serta mengeksplorasi hubungan antara indikator pelayanan kesehatan dan kejadian pneumonia balita. Seluruh analisis dilakukan menggunakan software R.
3.3.1 Analisis Deskriptif
Analisis deskriptif digunakan untuk menggambarkan karakteristik data pada tingkat provinsi, meliputi jumlah balita, jumlah kasus pneumonia balita, jumlah kematian akibat pneumonia, serta indikator pelayanan kesehatan. Hasil analisis deskriptif disajikan dalam bentuk tabel ringkasan, peta tematik, dan grafik perbandingan antar provinsi.
3.3.2 Ukuran Epidemiologi
a. Insidensi Kumulatif
Insidensi kumulatif pneumonia balita dihitung untuk menggambarkan risiko terjadinya pneumonia selama tahun 2024 pada setiap provinsi. Rumus insidensi kumulatif adalah:
\[ CI_i = \frac{C_i}{N_i} \]
dengan :
- \(C_i\) = jumlah kasus pneumonia
balita di provinsi ke-\(i\)
- \(N_i\) = jumlah balita di provinsi
ke-\(i\)
Insidensi kumulatif kemudian dinyatakan per 1.000 balita sebagai :
\[ CI_{1000,i} = \frac{C_i}{N_i} \times 1000 \]
b. Case Fatality Rate (CFR)
Case Fatality Rate (CFR) digunakan untuk menggambarkan tingkat keparahan pneumonia balita dan dihitung sebagai :
\[ CFR_i = \frac{D_i}{C_i} \times 100 \]
dengan :
- \(D_i\) = jumlah kematian akibat
pneumonia balita di provinsi ke-\(i\)
c. Mortality Rate Spesifik Pneumonia Balita
Angka kematian spesifik pneumonia balita dihitung sebagai :
\[ MR_{100000,i} = \frac{D_i}{N_i} \times 100000 \]
Ukuran ini menggambarkan risiko kematian akibat pneumonia balita pada populasi balita di setiap provinsi.
3.3.3 Analisis Spasial
Analisis spasial dilakukan untuk mengidentifikasi pola geografis dan kemungkinan kluster pneumonia balita antar provinsi.
a. Global Moran’s I
Global Moran’s I digunakan untuk menilai apakah distribusi insidensi pneumonia balita menunjukkan pola acak atau mengelompok secara spasial. Rumus Moran’s I adalah :
\[ I = \frac{n}{\sum_{i}\sum_{j} w_{ij}} \cdot \frac{\sum_{i}\sum_{j} w_{ij}(x_i - \bar{x})(x_j - \bar{x})} {\sum_{i}(x_i - \bar{x})^2} \]
dengan :
- \(x_i\) = nilai insidensi pneumonia
balita di wilayah \(i\)
- \(\bar{x}\) = rata-rata
insidensi
- \(w_{ij}\) = bobot spasial antara
wilayah \(i\) dan \(j\)
- \(n\) = jumlah wilayah
b. Local Indicators of Spatial Association (LISA)
Local Indicators of Spatial Association (LISA) digunakan untuk mengidentifikasi kluster lokal pneumonia balita, seperti hotspot dan coldspot. Nilai Local Moran’s I dihitung untuk setiap wilayah guna mengetahui kontribusi lokal terhadap autokorelasi spasial global.
3.3.4 Pemodelan Epidemiologi
a. Standardized Morbidity Ratio (SMR)
Dalam konteks pemetaan penyakit, risiko relatif antar provinsi dihitung menggunakan Standardized Morbidity Ratio (SMR) :
\[ SMR_i = \frac{O_i}{E_i} \]
dengan :
- \(O_i\) = jumlah kasus pneumonia
balita yang terobservasi
- \(E_i\) = jumlah kasus pneumonia
balita yang diharapkan
Jumlah kasus yang diharapkan dihitung sebagai :
\[ E_i = N_i \times r \]
di mana \(r\) merupakan laju insidensi nasional pneumonia balita.
b. Regresi Poisson
Untuk mengevaluasi hubungan antara indikator pelayanan kesehatan dan jumlah kasus pneumonia balita, digunakan model regresi Poisson dengan offset jumlah balita :
\[ Y_i \sim Poisson(\mu_i) \]
\[ \log(\mu_i) = \log(N_i) + \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} \]
dengan :
- \(Y_i\) = jumlah kasus pneumonia
balita
- \(N_i\) = jumlah balita
(offset)
- \(X_{1i}\) = persentase pengobatan
sesuai standar
- \(X_{2i}\) = cakupan pemberian
antibiotik
Hasil model dilaporkan dalam bentuk Incidence Rate Ratio (IRR) :
\[ IRR = e^{\beta} \]
c. Regresi Negative Binomial
Apabila ditemukan overdispersi, yaitu varians data lebih besar dibandingkan nilai rataan, maka model ditingkatkan menjadi regresi Negative Binomial untuk memperoleh estimasi parameter yang lebih stabil dan reliabel.
d. Empirical Bayes (EB) Smoothing
Untuk menstabilkan estimasi risiko pada provinsi dengan populasi balita kecil, dilakukan pendekatan Empirical Bayes smoothing terhadap nilai SMR. Pendekatan ini menghasilkan estimasi risiko relatif tersmoothing yang lebih andal untuk keperluan pemetaan penyakit.
3.3.5 Evaluasi Model
Evaluasi kinerja model dilakukan melalui beberapa pendekatan, meliputi: 1. Pemeriksaan overdispersi menggunakan rasio deviance terhadap derajat bebas. 2. Perbandingan model menggunakan Akaike Information Criterion (AIC). 3. Analisis residual dan uji Moran’s I terhadap residual untuk menilai adanya struktur spasial yang belum tertangkap oleh model.
3.4 Alur Kerja Analisis
Diagram alur kerja tersebut menunjukkan tahapan analisis pneumonia balita yang dilakukan secara sistematis dari awal hingga akhir penelitian. Proses dimulai dengan penggunaan data resmi dari Profil Kesehatan Indonesia Tahun 2024, kemudian data tersebut disusun menjadi dataset tingkat provinsi. Setelah itu dilakukan pra-pemrosesan untuk memastikan kualitas data, seperti pengecekan nilai yang hilang, kesesuaian tipe data, dan keseragaman penamaan provinsi, sehingga data siap digunakan untuk analisis lebih lanjut.
Tahap berikutnya adalah perhitungan ukuran epidemiologi untuk menggambarkan besarnya masalah pneumonia balita, seperti insidensi, tingkat kematian kasus (CFR), dan angka kematian spesifik. Data kemudian dipetakan dengan menggabungkan informasi kesehatan dan peta batas provinsi untuk menampilkan sebaran penyakit secara visual. Analisis spasial digunakan untuk melihat apakah terdapat pola pengelompokan kasus antar wilayah, sementara pemodelan dilakukan untuk memperkirakan risiko relatif dan hubungan antara faktor layanan kesehatan dengan jumlah kasus. Seluruh hasil analisis selanjutnya dievaluasi dan diinterpretasikan secara epidemiologis, lalu disajikan dalam bentuk laporan RMarkdown/RPubs yang mudah dibaca dan digunakan sebagai dasar pengambilan keputusan di bidang kesehatan.
4 Analisis dan Hasil
4.1 Deskripsi Kasus dan Proses Penyakit
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.5.1## Warning: package 'ggplot2' was built under R version 4.5.1## Warning: package 'tibble' was built under R version 4.5.1## Warning: package 'tidyr' was built under R version 4.5.1## Warning: package 'readr' was built under R version 4.5.1## Warning: package 'purrr' was built under R version 4.5.1## Warning: package 'dplyr' was built under R version 4.5.1## Warning: package 'stringr' was built under R version 4.5.1## Warning: package 'forcats' was built under R version 4.5.1## Warning: package 'lubridate' was built under R version 4.5.1## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(scales)## Warning: package 'scales' was built under R version 4.5.1##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(knitr)## Warning: package 'knitr' was built under R version 4.5.1library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.5.1##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowspne <- readr::read_csv(
"C:/Users/user/OneDrive/Documents/Epidem_DashboardUAS/Epidem_DataPneumonia.csv"
)## Rows: 38 Columns: 25
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): prov
## dbl (24): balita_n, pneu_est, pneu_u1_m, pneu_u1_f, pneu_1to4_m, pneu_1to4_f...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Pastikan tipe numerik benar
num_cols <- pne %>% dplyr::select(-prov) %>% names()
pne <- pne %>%
mutate(across(all_of(num_cols), ~as.numeric(.x)))
# Hitung ukuran turunan utama
pne <- pne %>%
mutate(
# Insidensi
ci = if_else(balita_n > 0, case_tot / balita_n, NA_real_),
ci_1000 = ci * 1000,
# Case Fatality Rate
cfr_calc = if_else(case_tot > 0, (death_tot / case_tot) * 100, NA_real_),
# Mortality rate per 100.000 balita
mort_100k = if_else(balita_n > 0, (death_tot / balita_n) * 100000, NA_real_),
# Cakupan antibiotik
abx_cov = if_else(case_tot > 0, (abx_n / case_tot) * 100, NA_real_)
)
head(pne)## # A tibble: 6 × 30
## prov balita_n pneu_est pneu_u1_m pneu_u1_f pneu_1to4_m pneu_1to4_f sev_u1_m
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aceh 548866 24479 403 286 1254 1000 75
## 2 Sumate… 1560648 46663 863 636 1824 1605 455
## 3 Sumate… 565241 22101 1151 842 3281 2871 68
## 4 Riau 644531 17209 319 198 654 557 32
## 5 Jambi 368234 11599 306 221 573 459 14
## 6 Sumate… 923818 33350 845 826 2587 2281 118
## # ℹ 22 more variables: sev_u1_f <dbl>, sev_1to4_m <dbl>, sev_1to4_f <dbl>,
## # case_m <dbl>, case_f <dbl>, case_tot <dbl>, case_pct <dbl>,
## # death_u1_m <dbl>, death_u1_f <dbl>, death_1to4_m <dbl>, death_1to4_f <dbl>,
## # death_m <dbl>, death_f <dbl>, death_tot <dbl>, cfr <dbl>, abx_n <dbl>,
## # tx_std_pct <dbl>, ci <dbl>, ci_1000 <dbl>, cfr_calc <dbl>, mort_100k <dbl>,
## # abx_cov <dbl>4.2 Perhitungan Ukuran Epidemiologi
4.2.1 Ukuran Frekuensi
Ringkasan Data Nasional
nasional <- pne %>%
summarise(
balita_nasional = sum(balita_n, na.rm = TRUE),
kasus_nasional = sum(case_tot, na.rm = TRUE),
kematian_nasional = sum(death_tot, na.rm = TRUE),
CI_1000_nasional = (kasus_nasional / balita_nasional) * 1000,
CFR_nasional = (kematian_nasional / kasus_nasional) * 100,
Mort_100k_nasional = (kematian_nasional / balita_nasional) * 100000
)
nasional %>%
mutate(across(where(is.numeric), ~round(.x, 2))) %>%
kable(caption = "Ringkasan Nasional Pneumonia Balita (Indonesia, 2024)") %>%
kable_styling(full_width = FALSE)| balita_nasional | kasus_nasional | kematian_nasional | CI_1000_nasional | CFR_nasional | Mort_100k_nasional |
|---|---|---|---|---|---|
| 26071952 | 530641 | 635 | 20.35 | 0.12 | 2.44 |
Ukuran Frekuensi per Provinsi
freq_prov <- pne %>%
transmute(
prov,
balita_n,
case_tot,
death_tot,
ci_1000 = round(ci_1000, 2),
cfr_pct = round(cfr_calc, 2),
mort_100k = round(mort_100k, 2),
tx_std_pct = round(tx_std_pct, 2),
abx_cov = round(abx_cov, 2)
) %>%
arrange(desc(ci_1000))
freq_prov %>%
kable(caption = "Ukuran Frekuensi Pneumonia Balita per Provinsi (2024)") %>%
kable_styling(full_width = FALSE) %>%
scroll_box(height = "400px")| prov | balita_n | case_tot | death_tot | ci_1000 | cfr_pct | mort_100k | tx_std_pct | abx_cov |
|---|---|---|---|---|---|---|---|---|
| DI Yogyakarta | 198343 | 12780 | 11 | 64.43 | 0.09 | 5.55 | 0.86 | 85.99 |
| DKI Jakarta | 891491 | 41032 | 0 | 46.03 | 0.00 | 0.00 | 1.00 | 100.00 |
| Kalimantan Selatan | 422507 | 13952 | 0 | 33.02 | 0.00 | 0.00 | 0.98 | 97.89 |
| Papua Selatan | 58181 | 1873 | 23 | 32.19 | 1.23 | 39.53 | 0.73 | 73.14 |
| Nusa Tenggara Barat | 564602 | 17883 | 9 | 31.67 | 0.05 | 1.59 | 0.99 | 99.22 |
| Banten | 1168992 | 35808 | 4 | 30.63 | 0.01 | 0.34 | 0.95 | 95.06 |
| Sulawesi Tengah | 294550 | 8840 | 5 | 30.01 | 0.06 | 1.70 | 0.92 | 91.89 |
| Jawa Timur | 3537130 | 99628 | 10 | 28.17 | 0.01 | 0.28 | 0.98 | 98.03 |
| Jawa Barat | 4960090 | 133926 | 91 | 27.00 | 0.07 | 1.83 | 0.98 | 97.57 |
| Gorontalo | 118164 | 2843 | 1 | 24.06 | 0.04 | 0.85 | 0.96 | 96.17 |
| Papua Barat | 50702 | 1203 | 7 | 23.73 | 0.58 | 13.81 | 0.96 | 95.84 |
| Jawa Tengah | 3547812 | 82862 | 275 | 23.36 | 0.33 | 7.75 | 0.97 | 97.49 |
| Kalimantan Utara | 73996 | 1524 | 3 | 20.60 | 0.20 | 4.05 | 0.99 | 98.75 |
| Bali | 323353 | 6055 | 11 | 18.73 | 0.18 | 3.40 | 0.98 | 97.79 |
| Kalimantan Timur | 365934 | 5799 | 1 | 15.85 | 0.02 | 0.27 | 0.99 | 99.40 |
| Kep. Bangka Belitung | 145333 | 2240 | 3 | 15.41 | 0.13 | 2.06 | 0.99 | 98.88 |
| Sumatera Barat | 565241 | 8576 | 1 | 15.17 | 0.01 | 0.18 | 0.96 | 95.64 |
| Papua | 92253 | 1153 | 8 | 12.50 | 0.69 | 8.67 | 0.92 | 91.67 |
| Sulawesi Selatan | 908767 | 9052 | 17 | 9.96 | 0.19 | 1.87 | 0.96 | 95.98 |
| Sulawesi Barat | 142964 | 1381 | 22 | 9.66 | 1.59 | 15.39 | 0.87 | 86.75 |
| Nusa Tenggara Timur | 540594 | 4966 | 2 | 9.19 | 0.04 | 0.37 | 1.00 | 100.00 |
| Lampung | 918823 | 7537 | 0 | 8.20 | 0.00 | 0.00 | 0.99 | 99.08 |
| Sumatera Selatan | 923818 | 6896 | 0 | 7.46 | 0.00 | 0.00 | 0.90 | 90.41 |
| Maluku Utara | 124666 | 870 | 2 | 6.98 | 0.23 | 1.60 | 0.96 | 96.21 |
| Papua Barat Daya | 56757 | 386 | 14 | 6.80 | 3.63 | 24.67 | 1.00 | 100.00 |
| Maluku | 184774 | 1121 | 11 | 6.07 | 0.98 | 5.95 | 0.98 | 98.04 |
| Kepulauan Riau | 207944 | 1255 | 3 | 6.04 | 0.24 | 1.44 | 0.97 | 97.37 |
| Sulawesi Tenggara | 273700 | 1647 | 3 | 6.02 | 0.18 | 1.10 | 1.00 | 99.57 |
| Aceh | 548866 | 3296 | 2 | 6.01 | 0.06 | 0.36 | 0.90 | 90.17 |
| Kalimantan Tengah | 275864 | 1480 | 6 | 5.36 | 0.41 | 2.17 | 1.00 | 100.00 |
| Jambi | 368234 | 1620 | 0 | 4.40 | 0.00 | 0.00 | 0.98 | 97.90 |
| Kalimantan Barat | 544043 | 2270 | 25 | 4.17 | 1.10 | 4.60 | 0.97 | 96.92 |
| Sumatera Utara | 1560648 | 6328 | 58 | 4.05 | 0.92 | 3.72 | 0.99 | 99.27 |
| Riau | 644531 | 1853 | 0 | 2.87 | 0.00 | 0.00 | 0.94 | 93.63 |
| Bengkulu | 207056 | 463 | 0 | 2.24 | 0.00 | 0.00 | 1.00 | 99.78 |
| Sulawesi Utara | 261229 | 243 | 7 | 0.93 | 2.88 | 2.68 | 0.99 | 98.77 |
| Papua Tengah | 0 | 0 | 0 | NA | NA | NA | 0.00 | NA |
| Papua Pegunungan | 0 | 0 | 0 | NA | NA | NA | 0.00 | NA |
Top 10 dan Bottom 10 Insidensi (per 1000)
top10_ci <- freq_prov %>% slice(1:10)
bottom10_ci <- freq_prov %>% arrange(ci_1000) %>% slice(1:10)
top10_ci %>%
kable(caption = "Top 10 Provinsi Insidensi Tertinggi (CI per 1.000)") %>%
kable_styling(full_width = FALSE)| prov | balita_n | case_tot | death_tot | ci_1000 | cfr_pct | mort_100k | tx_std_pct | abx_cov |
|---|---|---|---|---|---|---|---|---|
| DI Yogyakarta | 198343 | 12780 | 11 | 64.43 | 0.09 | 5.55 | 0.86 | 85.99 |
| DKI Jakarta | 891491 | 41032 | 0 | 46.03 | 0.00 | 0.00 | 1.00 | 100.00 |
| Kalimantan Selatan | 422507 | 13952 | 0 | 33.02 | 0.00 | 0.00 | 0.98 | 97.89 |
| Papua Selatan | 58181 | 1873 | 23 | 32.19 | 1.23 | 39.53 | 0.73 | 73.14 |
| Nusa Tenggara Barat | 564602 | 17883 | 9 | 31.67 | 0.05 | 1.59 | 0.99 | 99.22 |
| Banten | 1168992 | 35808 | 4 | 30.63 | 0.01 | 0.34 | 0.95 | 95.06 |
| Sulawesi Tengah | 294550 | 8840 | 5 | 30.01 | 0.06 | 1.70 | 0.92 | 91.89 |
| Jawa Timur | 3537130 | 99628 | 10 | 28.17 | 0.01 | 0.28 | 0.98 | 98.03 |
| Jawa Barat | 4960090 | 133926 | 91 | 27.00 | 0.07 | 1.83 | 0.98 | 97.57 |
| Gorontalo | 118164 | 2843 | 1 | 24.06 | 0.04 | 0.85 | 0.96 | 96.17 |
bottom10_ci %>%
kable(caption = "Bottom 10 Provinsi Insidensi Terendah (CI per 1.000)") %>%
kable_styling(full_width = FALSE)| prov | balita_n | case_tot | death_tot | ci_1000 | cfr_pct | mort_100k | tx_std_pct | abx_cov |
|---|---|---|---|---|---|---|---|---|
| Sulawesi Utara | 261229 | 243 | 7 | 0.93 | 2.88 | 2.68 | 0.99 | 98.77 |
| Bengkulu | 207056 | 463 | 0 | 2.24 | 0.00 | 0.00 | 1.00 | 99.78 |
| Riau | 644531 | 1853 | 0 | 2.87 | 0.00 | 0.00 | 0.94 | 93.63 |
| Sumatera Utara | 1560648 | 6328 | 58 | 4.05 | 0.92 | 3.72 | 0.99 | 99.27 |
| Kalimantan Barat | 544043 | 2270 | 25 | 4.17 | 1.10 | 4.60 | 0.97 | 96.92 |
| Jambi | 368234 | 1620 | 0 | 4.40 | 0.00 | 0.00 | 0.98 | 97.90 |
| Kalimantan Tengah | 275864 | 1480 | 6 | 5.36 | 0.41 | 2.17 | 1.00 | 100.00 |
| Aceh | 548866 | 3296 | 2 | 6.01 | 0.06 | 0.36 | 0.90 | 90.17 |
| Sulawesi Tenggara | 273700 | 1647 | 3 | 6.02 | 0.18 | 1.10 | 1.00 | 99.57 |
| Kepulauan Riau | 207944 | 1255 | 3 | 6.04 | 0.24 | 1.44 | 0.97 | 97.37 |
Top 10 CFR Tertinggi
pne %>%
transmute(prov, case_tot, death_tot, cfr_pct = round(cfr_calc, 2)) %>%
arrange(desc(cfr_pct)) %>%
slice(1:10) %>%
kable(caption = "Top 10 Provinsi CFR Tertinggi (2024)") %>%
kable_styling(full_width = FALSE)| prov | case_tot | death_tot | cfr_pct |
|---|---|---|---|
| Papua Barat Daya | 386 | 14 | 3.63 |
| Sulawesi Utara | 243 | 7 | 2.88 |
| Sulawesi Barat | 1381 | 22 | 1.59 |
| Papua Selatan | 1873 | 23 | 1.23 |
| Kalimantan Barat | 2270 | 25 | 1.10 |
| Maluku | 1121 | 11 | 0.98 |
| Sumatera Utara | 6328 | 58 | 0.92 |
| Papua | 1153 | 8 | 0.69 |
| Papua Barat | 1203 | 7 | 0.58 |
| Kalimantan Tengah | 1480 | 6 | 0.41 |
Grafik Distribusi CI dan CFR
library(ggplot2)
ggplot(freq_prov, aes(x = reorder(prov, ci_1000), y = ci_1000, fill = ci_1000)) +
geom_col(show.legend = FALSE) +
coord_flip() +
scale_fill_viridis_c(option = "plasma") +
labs(
title = "Insidensi Pneumonia Balita per Provinsi (2024)",
subtitle = "Jumlah kasus per 1.000 balita",
x = NULL,
y = "Kasus per 1.000 balita"
) +
theme_minimal(base_size = 11) +
theme(
axis.text.y = element_text(size = 5), # Label provinsi
axis.text.x = element_text(size = 9),
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
panel.grid.major.y = element_blank()
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_col()`).
library(ggplot2)
p_cfr <- pne %>%
transmute(prov, cfr_pct = cfr_calc) %>%
filter(is.finite(cfr_pct))
ggplot(p_cfr, aes(x = reorder(prov, cfr_pct), y = cfr_pct, fill = cfr_pct)) +
geom_col(show.legend = FALSE) +
coord_flip() +
scale_fill_viridis_c(option = "inferno") +
labs(
title = "Case Fatality Rate (CFR) Pneumonia Balita per Provinsi (2024)",
subtitle = "Persentase kematian dari seluruh kasus pneumonia balita",
x = NULL,
y = "CFR (%)"
) +
theme_minimal(base_size = 11) +
theme(
axis.text.y = element_text(size = 5), # Label provinsi
axis.text.x = element_text(size = 9),
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
panel.grid.major.y = element_blank()
)
4.2.2 Ukuran Asosiasi
library(tidyverse)
library(broom)## Warning: package 'broom' was built under R version 4.5.1library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select# Variabel turunan untuk asosiasi:
pne <- pne %>%
mutate(
abx_cov = if_else(case_tot > 0, abx_n / case_tot * 100, NA_real_),
log_balita = log(balita_n)
)
# Buat paparan kategorik agar RR mudah diinterpretasi
med_tx <- median(pne$tx_std_pct, na.rm = TRUE)
pne <- pne %>%
mutate(
tx_high = if_else(tx_std_pct > med_tx, 1, 0) # 1=di atas median
)
head(pne)## # A tibble: 6 × 32
## prov balita_n pneu_est pneu_u1_m pneu_u1_f pneu_1to4_m pneu_1to4_f sev_u1_m
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aceh 548866 24479 403 286 1254 1000 75
## 2 Sumate… 1560648 46663 863 636 1824 1605 455
## 3 Sumate… 565241 22101 1151 842 3281 2871 68
## 4 Riau 644531 17209 319 198 654 557 32
## 5 Jambi 368234 11599 306 221 573 459 14
## 6 Sumate… 923818 33350 845 826 2587 2281 118
## # ℹ 24 more variables: sev_u1_f <dbl>, sev_1to4_m <dbl>, sev_1to4_f <dbl>,
## # case_m <dbl>, case_f <dbl>, case_tot <dbl>, case_pct <dbl>,
## # death_u1_m <dbl>, death_u1_f <dbl>, death_1to4_m <dbl>, death_1to4_f <dbl>,
## # death_m <dbl>, death_f <dbl>, death_tot <dbl>, cfr <dbl>, abx_n <dbl>,
## # tx_std_pct <dbl>, ci <dbl>, ci_1000 <dbl>, cfr_calc <dbl>, mort_100k <dbl>,
## # abx_cov <dbl>, log_balita <dbl>, tx_high <dbl>RR/IRR dari Regresi Poisson (Offset Populasi)
m_pois <- glm(
case_tot ~ tx_high + abx_cov,
offset = log_balita,
family = poisson,
data = pne
)
irr_pois <- tidy(m_pois, exponentiate = TRUE, conf.int = TRUE) %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value)
irr_pois## # A tibble: 3 × 5
## term estimate conf.low conf.high p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.105 0.0940 0.118 0
## 2 tx_high 1.30 1.29 1.31 0
## 3 abx_cov 0.982 0.980 0.983 5.92e-207Overdispersi
disp <- sum(residuals(m_pois, type="pearson")^2) / df.residual(m_pois)
disp## [1] 5162.44Jika disp > 1.5
m_nb <- glm.nb(
case_tot ~ tx_high + abx_cov + offset(log_balita),
data = pne
)
irr_nb <- tidy(m_nb, exponentiate = TRUE, conf.int = TRUE) %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value)
irr_nb## # A tibble: 3 × 5
## term estimate conf.low conf.high p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.76 0.00732 1799. 0.849
## 2 tx_high 1.21 0.610 2.34 0.580
## 3 abx_cov 0.951 0.883 1.01 0.120Attributable Fraction (AF) dan Population Attributable Fraction (PAF)
# Pilih RR dari model utama (default: Poisson)
RR <- exp(coef(m_pois)["tx_high"])
Pe <- mean(pne$tx_high, na.rm = TRUE)
AF_e <- (RR - 1) / RR
PAF <- (Pe * (RR - 1)) / (Pe * (RR - 1) + 1)
tibble(
RR_tx_high = RR,
Pe = Pe,
AF_e = AF_e,
PAF = PAF
)## # A tibble: 1 × 4
## RR_tx_high Pe AF_e PAF
## <dbl> <dbl> <dbl> <dbl>
## 1 1.30 0.5 0.228 0.129Alternatif : RR untuk perubahan 10% pada tx_std_pct
(kontinu)
pne <- pne %>%
mutate(tx_std_10 = tx_std_pct / 10)
m_pois2 <- glm(
case_tot ~ tx_std_10 + abx_cov,
offset = log_balita,
family = poisson,
data = pne
)
tidy(m_pois2, exponentiate = TRUE, conf.int = TRUE) %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value)## # A tibble: 3 × 5
## term estimate conf.low conf.high p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.32e- 3 8.46e- 3 1.03e- 2 0
## 2 tx_std_10 1.83e+101 3.08e+97 1.09e+105 0
## 3 abx_cov 7.98e- 1 7.91e- 1 8.05e- 1 04.3 Visualisasi dan Interpretasi
library(sf)## Warning: package 'sf' was built under R version 4.5.1## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(dplyr)
library(leaflet)## Warning: package 'leaflet' was built under R version 4.5.1library(spdep)## Warning: package 'spdep' was built under R version 4.5.1## Loading required package: spData## Warning: package 'spData' was built under R version 4.5.1## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`# Path shapefile
shp_path <- "C:/Users/user/OneDrive/Documents/Epidem_DashboardUAS/Batas_Provinsi_2024.shp"
# Baca shapefile menjadi objek sf
prov_sf <- st_read(shp_path, quiet = TRUE)library(dplyr)
library(sf)
names(prov_sf)## [1] "KDPPUM" "WADMPR" "METADATA" "UPDATED" "geometry"# 1) Buat key di shapefile (pakai WADMPR)
prov_sf <- prov_sf %>%
mutate(
prov_key = WADMPR,
prov_key = case_when(
prov_key == "Daerah Istimewa Yogyakarta" ~ "DI Yogyakarta",
prov_key %in% c("Bangka Belitung", "Kepulauan Bangka Belitung") ~ "Kep. Bangka Belitung",
prov_key %in% c("DKI Jakarta", "Daerah Khusus Ibukota Jakarta") ~ "DKI Jakarta",
TRUE ~ prov_key
)
)
# 2) Buat key di data (pne)
pne <- pne %>%
mutate(
prov_key = case_when(
# Bangka Belitung
prov %in% c("Kep. Bangka Belitung", "Bangka Belitung", "Kepulauan Bangka Belitung")
~ "Kep. Bangka Belitung",
# Jakarta
prov %in% c("DKI Jakarta", "Jakarta Raya", "Daerah Khusus Ibukota Jakarta")
~ "DKI Jakarta",
# Yogyakarta
prov %in% c("DI Yogyakarta", "Yogyakarta", "D.I. Yogyakarta", "Daerah Istimewa Yogyakarta")
~ "DI Yogyakarta",
# Papua pemekaran (biarkan apa adanya)
prov %in% c("Papua Tengah", "Papua Pegunungan", "Papua Selatan", "Papua Barat Daya")
~ prov,
# default
TRUE ~ prov
)
)
# 3) Join ulang
map_df <- prov_sf %>%
left_join(pne, by = "prov_key")
# 4) Cek provinsi yang masih NA setelah join
map_df %>%
st_drop_geometry() %>%
filter(is.na(case_tot) | is.na(balita_n)) %>%
dplyr::select(prov_key, case_tot, balita_n)## [1] prov_key case_tot balita_n
## <0 rows> (or 0-length row.names)anti_join(
map_df %>% st_drop_geometry() %>% distinct(prov_key),
pne %>% distinct(prov_key),
by = "prov_key"
)## [1] prov_key
## <0 rows> (or 0-length row.names)anti_join(
pne %>% distinct(prov_key),
map_df %>% st_drop_geometry() %>% distinct(prov_key),
by = "prov_key"
)## # A tibble: 0 × 1
## # ℹ 1 variable: prov_key <chr>map_df %>%
st_drop_geometry() %>%
summarise(
n_total = n(),
n_na_case = sum(is.na(case_tot)),
n_na_balita = sum(is.na(balita_n))
)## n_total n_na_case n_na_balita
## 1 38 0 04.3.1 Peta Distribusi Penyakit
Peta Insidensi Pneumonia Balita (per 1000)
Case Fatality Rate (CFR)
# Pastikan CFR numerik
map_df <- map_df %>%
mutate(
cfr_calc = as.numeric(cfr_calc)
)# Paket minimal biar ringan
library(sf)
library(dplyr)
library(ggplot2)
library(scales)
library(viridisLite)## Warning: package 'viridisLite' was built under R version 4.5.1# Pastikan CRS (biar stabil)
# (Kalau sudah benar, ini aman dan tidak mengubah geometri)
map_df <- st_make_valid(map_df)
# Pastikan CFR numerik & finite
map_df <- map_df %>%
mutate(cfr_calc = as.numeric(cfr_calc)) %>%
mutate(cfr_calc = if_else(is.finite(cfr_calc), cfr_calc, NA_real_))
# --- Klasifikasi CFR (quantile) biar choroplethnya rapi ---
# fallback kalau ada banyak nilai sama / data kecil
breaks_q <- tryCatch(
quantile(map_df$cfr_calc, probs = seq(0, 1, 0.2), na.rm = TRUE),
error = function(e) NULL
)
# Kalau quantile gagal / terlalu banyak duplikat, pakai pretty breaks
if (is.null(breaks_q) || length(unique(breaks_q)) < 3) {
breaks_q <- pretty(map_df$cfr_calc, n = 5)
}
breaks_q <- unique(as.numeric(breaks_q))
breaks_q <- breaks_q[is.finite(breaks_q)]
map_df <- map_df %>%
mutate(
cfr_cat = cut(
cfr_calc,
breaks = breaks_q,
include.lowest = TRUE,
dig.lab = 5
)
)
# --- Split untuk inset (Papua) supaya peta tidak "ketarik" jauh ---
papua_keys <- c(
"Papua", "Papua Barat",
"Papua Tengah", "Papua Pegunungan", "Papua Selatan", "Papua Barat Daya"
)
map_main <- map_df %>% filter(!prov_key %in% papua_keys)
map_papua <- map_df %>% filter(prov_key %in% papua_keys)ggplot(map_df) +
geom_sf(aes(fill = cfr_cat), color = "white", linewidth = 0.15) +
scale_fill_viridis_d(
option = "magma",
direction = -1,
na.value = "grey90",
name = "CFR (%)\n(kuantil)"
) +
labs(
title = "Peta Case Fatality Rate (CFR) Pneumonia Balita per Provinsi — Indonesia, 2024",
subtitle = "Choropleth kategori kuantil; warna lebih gelap = CFR lebih tinggi",
caption = "Sumber: Data Pneumonia Balita 2024 (olah penulis)."
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
panel.grid = element_blank()
)
Standardized Morbidity Ratio (SMR)
rate_nat <- sum(pne$case_tot, na.rm = TRUE) / sum(pne$balita_n, na.rm = TRUE)
map_df <- map_df %>%
mutate(
expected = as.numeric(balita_n) * rate_nat,
smr = if_else(expected > 0, as.numeric(case_tot) / expected, NA_real_)
)Peta Empirical Bayes (EB) Smoothed Relative Risk
library(sf)
library(spdep)
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched off# geometri super-ringan hanya untuk membentuk nb
map_nb <- map_df
st_geometry(map_nb) <- st_simplify(
st_make_valid(st_geometry(map_df)),
dTolerance = 0.05, # naikkan = makin cepat
preserveTopology = TRUE
)## Warning in st_simplify.sfc(st_make_valid(st_geometry(map_df)), dTolerance =
## 0.05, : st_simplify does not correctly simplify longitude/latitude data,
## dTolerance needs to be in decimal degreessystem.time(nb <- poly2nb(map_nb, queen = TRUE))## although coordinates are longitude/latitude, st_intersects assumes that they
## are planar## Warning in poly2nb(map_nb, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(map_nb, queen = TRUE): neighbour object has 13 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.## user system elapsed
## 1.48 0.13 1.65system.time(lw <- nb2listw(nb, style = "W", zero.policy = TRUE))## user system elapsed
## 0 0 0sf::sf_use_s2(TRUE)## Spherical geometry (s2) switched onlength(nb) # 38## [1] 38length(lw$neighbours) # 38## [1] 38summary(card(nb)) # jumlah tetangga tiap provinsi## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 2.000 1.737 3.000 4.000library(dplyr)
# Pastikan expected sudah ada
# rate_nat <- sum(pne$case_tot, na.rm = TRUE) / sum(pne$balita_n, na.rm = TRUE)
# map_df <- map_df %>% mutate(expected = as.numeric(balita_n) * rate_nat)
# 1) Siapkan O dan E (buang E=0)
eb_dat <- map_df %>%
mutate(
O = as.numeric(case_tot),
E = as.numeric(expected)
) %>%
mutate(
smr = if_else(E > 0, O / E, NA_real_)
)
ok <- is.finite(eb_dat$smr) & eb_dat$E > 0
O <- eb_dat$O[ok]
E <- eb_dat$E[ok]
smr <- eb_dat$smr[ok]
# 2) Rata-rata nasional (mu = total O / total E)
mu <- sum(O) / sum(E)
# 3) Estimasi heterogenitas tau^2 (method of moments)
# Var(SMR_i) ≈ mu/E_i + tau^2
tau2 <- var(smr) - mean(mu / E)
tau2 <- max(tau2, 1e-8) # jaga agar >0
# 4) Parameter prior Gamma(alpha, beta) dengan mean = mu dan var = tau2
alpha <- (mu^2) / tau2
beta <- mu / tau2
# 5) Posterior mean EB untuk RR
rr_eb_global <- (O + alpha) / (E + beta)
# 6) Simpan balik ke map_df (provinsi E=0 -> NA)
map_df$rr_eb <- NA_real_
map_df$rr_eb[which(ok)] <- rr_eb_global
# Pastikan rr_eb sudah ada dari cell permodelan (Global EB)
# (opsional) cek cepat
class(map_df$rr_eb)## [1] "numeric"summary(map_df$rr_eb)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.04608 0.29664 0.55198 0.81754 1.21813 3.16469 2# Domain bersih untuk palet
domain_eb <- map_df$rr_eb
domain_eb <- domain_eb[is.finite(domain_eb)]
length(domain_eb)## [1] 36summary(domain_eb)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04608 0.29664 0.55198 0.81754 1.21813 3.16469library(sf)
library(dplyr)
library(ggplot2)
library(scales)
library(viridisLite)
# Pastikan valid & rr_eb numeric
map_df <- map_df %>%
mutate(rr_eb = as.numeric(rr_eb)) %>%
mutate(rr_eb = if_else(is.finite(rr_eb), rr_eb, NA_real_))
map_df <- st_make_valid(map_df)
# Kategori kuantil untuk RR EB (lebih stabil & enak dibaca)
breaks_eb <- tryCatch(
quantile(map_df$rr_eb, probs = seq(0, 1, 0.2), na.rm = TRUE),
error = function(e) NULL
)
if (is.null(breaks_eb) || length(unique(breaks_eb)) < 3) {
breaks_eb <- pretty(map_df$rr_eb, n = 5)
}
breaks_eb <- unique(as.numeric(breaks_eb))
breaks_eb <- breaks_eb[is.finite(breaks_eb)]
map_df <- map_df %>%
mutate(
rr_eb_cat = cut(
rr_eb,
breaks = breaks_eb,
include.lowest = TRUE,
dig.lab = 5
)
)
# Inset Papua biar peta nggak ketarik panjang
papua_keys <- c(
"Papua", "Papua Barat",
"Papua Tengah", "Papua Pegunungan", "Papua Selatan", "Papua Barat Daya"
)
map_main <- map_df %>% filter(!prov_key %in% papua_keys)
map_papua <- map_df %>% filter(prov_key %in% papua_keys)
# Untuk judul/caption: ringkas range RR EB (opsional)
rr_rng <- range(map_df$rr_eb, na.rm = TRUE)
rr_rng## [1] 0.0460832 3.1646914# Pakai palet yang "diverging vibe" tapi tetap ringan (viridis option 'turbo' / 'rocket' tidak ada di viridisLite base)
# Jadi kita gunakan viridis_d option "C" (default viridis) + arah warna, dan beri garis referensi di caption.
p_main <- ggplot() +
geom_sf(
data = map_main,
aes(fill = rr_eb_cat),
color = "white",
linewidth = 0.15
) +
scale_fill_viridis_d(
option = "viridis",
direction = -1,
na.value = "grey90",
name = "RR EB\n(kuantil)"
) +
labs(
title = "Peta Empirical Bayes (EB) Smoothed Relative Risk Pneumonia Balita — Indonesia, 2024",
subtitle = "RR EB dibanding rata-rata nasional (RR=1 ≈ setara nasional); kategori berdasarkan kuantil",
caption = "Sumber: Data Pneumonia Balita 2024 (olah penulis). EB smoothing: prior global (Gamma) dari O dan E."
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
legend.position = "right",
panel.grid = element_blank()
)
# Inset Papua (opsional, ringan)
if (nrow(map_papua) > 0) {
p_inset <- ggplot() +
geom_sf(
data = map_papua,
aes(fill = rr_eb_cat),
color = "white",
linewidth = 0.15
) +
scale_fill_viridis_d(
option = "viridis",
direction = -1,
na.value = "grey90"
) +
theme_void() +
theme(legend.position = "none")
library(grid)
print(p_main)
vp <- viewport(
x = 0.86, y = 0.18, width = 0.28, height = 0.28,
just = c("center", "center")
)
print(p_inset, vp = vp)
} else {
print(p_main)
}
4.3.2 Pola Spasial dan Cluster Penyakit
Global Moran’s I (Autokorelasi Spasial)
Pada analisis ini, variabel yang diuji adalah insidensi kumulatif per
1.000 balita (ci_1000).
Hipotesis uji Moran’s I :
- H0 : tidak ada autokorelasi spasial (pola acak)
- H1 : ada autokorelasi spasial (pola mengelompok/tersebar)
library(spdep)
# Pastikan tidak ada NA di indikator yang diuji
x <- map_df$ci_1000
ok <- is.finite(x)
# Tetangga & bobot
# nb <- poly2nb(map_df)
# lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
moran_global <- moran.test(x[ok], nb2listw(subset(nb, ok), style="W", zero.policy=TRUE), zero.policy = TRUE)## Warning in subset.nb(nb, ok): subsetting caused increase in subgraph countmoran_global##
## Moran I test under randomisation
##
## data: x[ok]
## weights: nb2listw(subset(nb, ok), style = "W", zero.policy = TRUE)
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 1.0761, p-value = 0.1409
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1601943 -0.0400000 0.0346122LISA (Local Moran’s I) dan Klasifikasi Hotspot/Coldspot
library(dplyr)
library(sf)
library(spdep)
# 1) pastikan numeric
map_df <- map_df %>%
mutate(ci_1000 = as.numeric(ci_1000))
# 2) buat flag valid (provinsi yang bisa dihitung LISA)
map_df <- map_df %>%
mutate(ci_valid = is.finite(ci_1000))
# 3) z-score hanya dari yang valid
mu <- mean(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
sdv <- sd(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
map_df <- map_df %>%
mutate(
ci_z = (ci_1000 - mu) / sdv,
# agar panjang tetap sama, isi 0 untuk yang invalid (mean z-score)
ci_z = if_else(ci_valid, ci_z, 0)
)
# 4) spatial lag (sekarang aman karena tidak ada NA/Inf)
map_df$lag_ci_z <- lag.listw(lw, map_df$ci_z, zero.policy = TRUE)
# 5) Local Moran (jalan karena input finite semua)
lisa <- localmoran(map_df$ci_z, lw, zero.policy = TRUE)
colnames(lisa)## [1] "Ii" "E.Ii" "Var.Ii" "Z.Ii"
## [5] "Pr(z != E(Ii))"pcol <- grep("^Pr\\(", colnames(lisa), value = TRUE)[1]
map_df$lisa_I <- lisa[, "Ii"]
map_df$lisa_p <- lisa[, pcol]
# 6) klasifikasi cluster (hanya untuk yang valid, sisanya "Tidak dihitung")
alpha <- 0.05
map_df <- map_df %>%
mutate(
cluster = case_when(
!ci_valid ~ "Tidak Dihitung (CI NA/Inf)",
lisa_p < alpha & ci_z > 0 & lag_ci_z > 0 ~ "Hotspot (High-High)",
lisa_p < alpha & ci_z < 0 & lag_ci_z < 0 ~ "Coldspot (Low-Low)",
lisa_p < alpha & ci_z > 0 & lag_ci_z < 0 ~ "Outlier (High-Low)",
lisa_p < alpha & ci_z < 0 & lag_ci_z > 0 ~ "Outlier (Low-High)",
TRUE ~ "Tidak signifikan"
)
)
# 7) ringkasan
map_df %>%
st_drop_geometry() %>%
count(cluster, sort = TRUE)## cluster n
## 1 Tidak signifikan 35
## 2 Tidak Dihitung (CI NA/Inf) 2
## 3 Hotspot (High-High) 1Peta Cluster LISA
library(dplyr)
library(sf)
library(spdep)
library(leaflet)
# =========================================================
# Hitung LISA
# =========================================================
map_df <- map_df %>%
mutate(ci_1000 = as.numeric(ci_1000),
ci_valid = is.finite(ci_1000))
mu <- mean(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
sdv <- sd(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
map_df <- map_df %>%
mutate(
ci_z = (ci_1000 - mu) / sdv,
ci_z = if_else(ci_valid, ci_z, 0)
)
map_df$lag_ci_z <- lag.listw(lw, map_df$ci_z, zero.policy = TRUE)
lisa <- localmoran(map_df$ci_z, lw, zero.policy = TRUE)
pcol <- grep("^Pr\\(", colnames(lisa), value = TRUE)[1]
map_df$lisa_p <- lisa[, pcol]
alpha <- 0.05
map_df <- map_df %>%
mutate(
cluster = case_when(
!ci_valid ~ "Tidak dihitung",
lisa_p < alpha & ci_z > 0 & lag_ci_z > 0 ~ "Hotspot (HH)",
lisa_p < alpha & ci_z < 0 & lag_ci_z < 0 ~ "Coldspot (LL)",
lisa_p < alpha & ci_z > 0 & lag_ci_z < 0 ~ "Outlier (HL)",
lisa_p < alpha & ci_z < 0 & lag_ci_z > 0 ~ "Outlier (LH)",
TRUE ~ "Tidak signifikan"
)
)
# =========================================================
# Objek view
# =========================================================
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched offmap_view <- map_df
st_geometry(map_view) <- st_simplify(
st_make_valid(st_geometry(map_df)),
dTolerance = 0.08, # lebih ringan (kalau masih berat: 0.1)
preserveTopology = TRUE
)## Warning in st_simplify.sfc(st_make_valid(st_geometry(map_df)), dTolerance =
## 0.08, : st_simplify does not correctly simplify longitude/latitude data,
## dTolerance needs to be in decimal degreessf::sf_use_s2(TRUE)## Spherical geometry (s2) switched onlibrary(sf)
library(dplyr)
library(spdep)
library(ggplot2)
# =========================================================
# HITUNG LISA (kalau kolom 'cluster' belum ada / mau rerun)
# =========================================================
map_df <- map_df %>%
mutate(
ci_1000 = as.numeric(ci_1000),
ci_valid = is.finite(ci_1000)
)
mu <- mean(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
sdv <- sd(map_df$ci_1000[map_df$ci_valid], na.rm = TRUE)
map_df <- map_df %>%
mutate(
ci_z = (ci_1000 - mu) / sdv,
ci_z = if_else(ci_valid, ci_z, 0) # agar localmoran tidak error
)
# Spatial lag z-score
map_df$lag_ci_z <- lag.listw(lw, map_df$ci_z, zero.policy = TRUE)
# Local Moran
lisa <- localmoran(map_df$ci_z, lw, zero.policy = TRUE)
pcol <- grep("^Pr\\(", colnames(lisa), value = TRUE)[1]
map_df$lisa_p <- lisa[, pcol]
alpha <- 0.05
map_df <- map_df %>%
mutate(
cluster = case_when(
!ci_valid ~ "Tidak dihitung",
lisa_p < alpha & ci_z > 0 & lag_ci_z > 0 ~ "Hotspot (HH)",
lisa_p < alpha & ci_z < 0 & lag_ci_z < 0 ~ "Coldspot (LL)",
lisa_p < alpha & ci_z > 0 & lag_ci_z < 0 ~ "Outlier (HL)",
lisa_p < alpha & ci_z < 0 & lag_ci_z > 0 ~ "Outlier (LH)",
TRUE ~ "Tidak signifikan"
),
cluster = factor(cluster, levels = c(
"Hotspot (HH)",
"Coldspot (LL)",
"Outlier (HL)",
"Outlier (LH)",
"Tidak signifikan",
"Tidak dihitung"
))
)
# Ringkasan cepat (opsional, buat narasi hasil)
tab_lisa <- map_df %>% st_drop_geometry() %>% count(cluster, sort = TRUE)
tab_lisa## cluster n
## 1 Tidak signifikan 35
## 2 Tidak dihitung 2
## 3 Hotspot (HH) 1# =========================================================
# GEOMETRI RINGAN UNTUK PLOT (biar RPubs aman)
# =========================================================
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched offmap_view <- map_df
st_geometry(map_view) <- st_simplify(
st_make_valid(st_geometry(map_df)),
dTolerance = 0.08, # kalau masih berat: 0.10
preserveTopology = TRUE
)## Warning in st_simplify.sfc(st_make_valid(st_geometry(map_df)), dTolerance =
## 0.08, : st_simplify does not correctly simplify longitude/latitude data,
## dTolerance needs to be in decimal degreessf::sf_use_s2(TRUE)## Spherical geometry (s2) switched on# =========================================================
# PALET WARNA (manual biar cluster gampang dibaca)
# =========================================================
pal_cluster <- c(
"Hotspot (HH)" = "#D73027", # merah
"Coldspot (LL)" = "#4575B4", # biru
"Outlier (HL)" = "#FC8D59", # oranye
"Outlier (LH)" = "#91BFDB", # biru muda
"Tidak signifikan" = "grey85",
"Tidak dihitung" = "grey95"
)
# =========================================================
# PLOT
# =========================================================
ggplot(map_view) +
geom_sf(aes(fill = cluster), color = "white", linewidth = 0.15) +
scale_fill_manual(values = pal_cluster, drop = FALSE, name = "Kelas LISA") +
labs(
title = "Peta LISA (Local Moran’s I) — Klaster Insidensi Pneumonia Balita (CI per 1.000), Indonesia 2024",
subtitle = paste0(
"Alpha = ", alpha,
" | Hotspot (HH) = tinggi dikelilingi tinggi; Coldspot (LL) = rendah dikelilingi rendah"
),
caption = "Metode: Local Moran’s I pada z-score CI per provinsi dengan bobot spasial (queen contiguity)."
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
panel.grid = element_blank(),
legend.position = "right"
)
4.3.3 Perbandingan Wilayah
Top–Bottom Provinsi Berdasarkan Insidensi (CI per 1.000)
library(dplyr)
library(knitr)
library(kableExtra)
tb_ci <- pne %>%
transmute(
prov,
ci_1000 = round(ci_1000, 2),
case_tot,
balita_n
) %>%
arrange(desc(ci_1000))
top_ci <- tb_ci %>% slice(1:10)
bot_ci <- tb_ci %>% arrange(ci_1000) %>% slice(1:10)
top_ci %>%
kable(caption = "Top 10 Provinsi Insidensi Pneumonia Balita Tertinggi (per 1.000)") %>%
kable_styling(full_width = FALSE)| prov | ci_1000 | case_tot | balita_n |
|---|---|---|---|
| DI Yogyakarta | 64.43 | 12780 | 198343 |
| DKI Jakarta | 46.03 | 41032 | 891491 |
| Kalimantan Selatan | 33.02 | 13952 | 422507 |
| Papua Selatan | 32.19 | 1873 | 58181 |
| Nusa Tenggara Barat | 31.67 | 17883 | 564602 |
| Banten | 30.63 | 35808 | 1168992 |
| Sulawesi Tengah | 30.01 | 8840 | 294550 |
| Jawa Timur | 28.17 | 99628 | 3537130 |
| Jawa Barat | 27.00 | 133926 | 4960090 |
| Gorontalo | 24.06 | 2843 | 118164 |
bot_ci %>%
kable(caption = "Bottom 10 Provinsi Insidensi Pneumonia Balita Terendah (per 1.000)") %>%
kable_styling(full_width = FALSE)| prov | ci_1000 | case_tot | balita_n |
|---|---|---|---|
| Sulawesi Utara | 0.93 | 243 | 261229 |
| Bengkulu | 2.24 | 463 | 207056 |
| Riau | 2.87 | 1853 | 644531 |
| Sumatera Utara | 4.05 | 6328 | 1560648 |
| Kalimantan Barat | 4.17 | 2270 | 544043 |
| Jambi | 4.40 | 1620 | 368234 |
| Kalimantan Tengah | 5.36 | 1480 | 275864 |
| Aceh | 6.01 | 3296 | 548866 |
| Sulawesi Tenggara | 6.02 | 1647 | 273700 |
| Kepulauan Riau | 6.04 | 1255 | 207944 |
Top–Bottom Provinsi Berdasarkan Case Fatality Rate (CFR)
tb_cfr <- pne %>%
transmute(
prov,
cfr_pct = round(cfr_calc, 2),
case_tot,
death_tot
) %>%
filter(is.finite(cfr_pct)) %>%
arrange(desc(cfr_pct))
top_cfr <- tb_cfr %>% slice(1:10)
bot_cfr <- tb_cfr %>% arrange(cfr_pct) %>% slice(1:10)
top_cfr %>%
kable(caption = "Top 10 Provinsi CFR Pneumonia Balita Tertinggi (%)") %>%
kable_styling(full_width = FALSE)| prov | cfr_pct | case_tot | death_tot |
|---|---|---|---|
| Papua Barat Daya | 3.63 | 386 | 14 |
| Sulawesi Utara | 2.88 | 243 | 7 |
| Sulawesi Barat | 1.59 | 1381 | 22 |
| Papua Selatan | 1.23 | 1873 | 23 |
| Kalimantan Barat | 1.10 | 2270 | 25 |
| Maluku | 0.98 | 1121 | 11 |
| Sumatera Utara | 0.92 | 6328 | 58 |
| Papua | 0.69 | 1153 | 8 |
| Papua Barat | 0.58 | 1203 | 7 |
| Kalimantan Tengah | 0.41 | 1480 | 6 |
bot_cfr %>%
kable(caption = "Bottom 10 Provinsi CFR Pneumonia Balita Terendah (%)") %>%
kable_styling(full_width = FALSE)| prov | cfr_pct | case_tot | death_tot |
|---|---|---|---|
| Riau | 0.00 | 1853 | 0 |
| Jambi | 0.00 | 1620 | 0 |
| Sumatera Selatan | 0.00 | 6896 | 0 |
| Bengkulu | 0.00 | 463 | 0 |
| Lampung | 0.00 | 7537 | 0 |
| DKI Jakarta | 0.00 | 41032 | 0 |
| Kalimantan Selatan | 0.00 | 13952 | 0 |
| Sumatera Barat | 0.01 | 8576 | 1 |
| Jawa Timur | 0.01 | 99628 | 10 |
| Banten | 0.01 | 35808 | 4 |
Perbandingan Risiko Relatif (SMR / EB-smoothed RR)
# Pilih indikator risiko relatif yang tersedia
rr_var <- if ("rr_eb" %in% names(map_df)) "rr_eb" else "smr"
tb_rr <- map_df %>%
st_drop_geometry() %>%
transmute(
prov,
rr = round(.data[[rr_var]], 2)
) %>%
filter(is.finite(rr)) %>%
arrange(desc(rr))
top_rr <- tb_rr %>% slice(1:10)
bot_rr <- tb_rr %>% arrange(rr) %>% slice(1:10)
top_rr %>%
kable(caption = "Top 10 Provinsi Risiko Relatif Tertinggi") %>%
kable_styling(full_width = FALSE)| prov | rr |
|---|---|
| DI Yogyakarta | 3.16 |
| DKI Jakarta | 2.26 |
| Kalimantan Selatan | 1.62 |
| Papua Selatan | 1.58 |
| Nusa Tenggara Barat | 1.56 |
| Banten | 1.50 |
| Sulawesi Tengah | 1.47 |
| Jawa Timur | 1.38 |
| Jawa Barat | 1.33 |
| Gorontalo | 1.18 |
bot_rr %>%
kable(caption = "Bottom 10 Provinsi Risiko Relatif Terendah") %>%
kable_styling(full_width = FALSE)| prov | rr |
|---|---|
| Sulawesi Utara | 0.05 |
| Bengkulu | 0.11 |
| Riau | 0.14 |
| Sumatera Utara | 0.20 |
| Kalimantan Barat | 0.21 |
| Jambi | 0.22 |
| Kalimantan Tengah | 0.26 |
| Aceh | 0.30 |
| Kepulauan Riau | 0.30 |
| Sulawesi Tenggara | 0.30 |
Grafik Distribusi Insidensi & CFR
library(ggplot2)
ggplot(top_ci, aes(x = reorder(prov, ci_1000), y = ci_1000)) +
geom_col() +
coord_flip() +
labs(
title = "Top 10 Provinsi Insidensi Pneumonia Balita Tertinggi",
x = NULL, y = "Kasus per 1.000 balita"
) +
theme_minimal()
ggplot(top_cfr, aes(x = reorder(prov, cfr_pct), y = cfr_pct)) +
geom_col() +
coord_flip() +
labs(
title = "Top 10 Provinsi CFR Pneumonia Balita Tertinggi",
x = NULL, y = "CFR (%)"
) +
theme_minimal()
4.4 Pemodelan Epidemiologi
library(tidyverse)
library(MASS)
library(broom)
library(spdep)
library(sf)
# Baca data
pne <- read_csv("Epidem_DataPneumonia.csv")## Rows: 38 Columns: 25
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): prov
## dbl (24): balita_n, pneu_est, pneu_u1_m, pneu_u1_f, pneu_1to4_m, pneu_1to4_f...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Pastikan numerik
num_cols <- pne %>% dplyr::select(-prov) %>% names()
pne <- pne %>% mutate(across(all_of(num_cols), as.numeric))
# Variabel turunan
pne <- pne %>%
mutate(
abx_cov = if_else(case_tot > 0, abx_n / case_tot * 100, NA_real_),
log_balita = log(balita_n)
)4.4.1 Pemilihan Model dan Asumsi
Model Standardized Morbidity Ratio (SMR)
# Laju insidensi nasional
r_nasional <- sum(pne$case_tot, na.rm = TRUE) /
sum(pne$balita_n, na.rm = TRUE)
# Hitung expected cases & SMR
pne <- pne %>%
mutate(
expected = balita_n * r_nasional,
smr = case_tot / expected
)
# Ringkasan SMR
pne %>%
summarise(
smr_min = min(smr, na.rm = TRUE),
smr_mean = mean(smr, na.rm = TRUE),
smr_max = max(smr, na.rm = TRUE)
)## # A tibble: 1 × 3
## smr_min smr_mean smr_max
## <dbl> <dbl> <dbl>
## 1 0.0457 0.817 3.17Model Regresi Poisson (Offset Populasi)
m_pois <- glm(
case_tot ~ tx_std_pct + abx_cov,
offset = log_balita,
family = poisson,
data = pne
)
summary(m_pois)##
## Call:
## glm(formula = case_tot ~ tx_std_pct + abx_cov, family = poisson,
## data = pne, offset = log_balita)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.675513 0.049376 -94.69 <2e-16 ***
## tx_std_pct 23.316495 0.443381 52.59 <2e-16 ***
## abx_cov -0.225133 0.004478 -50.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 198606 on 35 degrees of freedom
## Residual deviance: 195518 on 33 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 195891
##
## Number of Fisher Scoring iterations: 5IRR (Incidence Rate Ratio)
irr_pois <- tidy(m_pois, exponentiate = TRUE, conf.int = TRUE)
irr_pois %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value)## # A tibble: 3 × 5
## term estimate conf.low conf.high p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.32e- 3 8.46e-3 1.03e- 2 0
## 2 tx_std_pct 1.34e+10 5.61e+9 3.19e+10 0
## 3 abx_cov 7.98e- 1 7.91e-1 8.05e- 1 0Cek Overdispersi → Model Regresi Negative Binomial
dispersion <- sum(residuals(m_pois, type = "pearson")^2) /
df.residual(m_pois)
dispersion## [1] 5329.5# Jika `dispersion > 1.5`, lanjutkan dengan Negative Binomial:
m_nb <- glm.nb(
case_tot ~ tx_std_pct + abx_cov + offset(log_balita),
data = pne
)
summary(m_nb)##
## Call:
## glm.nb(formula = case_tot ~ tx_std_pct + abx_cov + offset(log_balita),
## data = pne, init.theta = 1.582465775, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4420 2.4190 -0.183 0.855
## tx_std_pct -17.7201 55.6937 -0.318 0.750
## abx_cov 0.1389 0.5588 0.248 0.804
##
## (Dispersion parameter for Negative Binomial(1.5825) family taken to be 1)
##
## Null deviance: 42.454 on 35 degrees of freedom
## Residual deviance: 39.653 on 33 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 701.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.582
## Std. Err.: 0.341
##
## 2 x log-likelihood: -693.195IRR Negative Binomial
irr_nb <- tidy(m_nb, exponentiate = TRUE, conf.int = TRUE)
irr_nb %>%
dplyr::select(term, estimate, conf.low, conf.high, p.value)## # A tibble: 3 × 5
## term estimate conf.low conf.high p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.643 8.57e- 3 1.36e 2 0.855
## 2 tx_std_pct 0.0000000202 8.66e-60 3.44e42 0.750
## 3 abx_cov 1.15 3.61e- 1 3.76e 0 0.804# Bandingkan Model (AIC)
AIC(m_pois, m_nb)## df AIC
## m_pois 3 195890.682
## m_nb 4 701.195Model Empirical Bayes (EB) Smoothing untuk SMR
library(dplyr)
# Pastikan expected sudah ada
# rate_nat <- sum(pne$case_tot, na.rm = TRUE) / sum(pne$balita_n, na.rm = TRUE)
# map_df <- map_df %>% mutate(expected = as.numeric(balita_n) * rate_nat)
# 1) Siapkan O dan E (buang E=0)
eb_dat <- map_df %>%
mutate(
O = as.numeric(case_tot),
E = as.numeric(expected)
) %>%
mutate(
smr = if_else(E > 0, O / E, NA_real_)
)
ok <- is.finite(eb_dat$smr) & eb_dat$E > 0
O <- eb_dat$O[ok]
E <- eb_dat$E[ok]
smr <- eb_dat$smr[ok]
# 2) Rata-rata nasional (mu = total O / total E)
mu <- sum(O) / sum(E)
# 3) Estimasi heterogenitas tau^2 (method of moments)
# Var(SMR_i) ≈ mu/E_i + tau^2
tau2 <- var(smr) - mean(mu / E)
tau2 <- max(tau2, 1e-8) # jaga agar >0
# 4) Parameter prior Gamma(alpha, beta) dengan mean = mu dan var = tau2
alpha <- (mu^2) / tau2
beta <- mu / tau2
# 5) Posterior mean EB untuk RR
rr_eb_global <- (O + alpha) / (E + beta)
# 6) Simpan balik ke map_df (provinsi E=0 -> NA)
map_df$rr_eb <- NA_real_
map_df$rr_eb[which(ok)] <- rr_eb_global
summary(map_df$rr_eb)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.04608 0.29664 0.55198 0.81754 1.21813 3.16469 2Output
map_df %>%
st_drop_geometry() %>%
transmute(
prov,
case_tot,
expected,
smr = round(smr, 2),
rr_eb = round(rr_eb, 2)
) %>%
arrange(desc(rr_eb)) %>%
head(10)## prov case_tot expected smr rr_eb
## 1 DI Yogyakarta 12780 4036.864 3.17 3.16
## 2 DKI Jakarta 41032 18144.467 2.26 2.26
## 3 Kalimantan Selatan 13952 8599.262 1.62 1.62
## 4 Papua Selatan 1873 1184.155 1.58 1.58
## 5 Nusa Tenggara Barat 17883 11491.313 1.56 1.56
## 6 Banten 35808 23792.430 1.51 1.50
## 7 Sulawesi Tengah 8840 5994.960 1.47 1.47
## 8 Jawa Timur 99628 71991.012 1.38 1.38
## 9 Jawa Barat 133926 100952.438 1.33 1.33
## 10 Gorontalo 2843 2404.985 1.18 1.184.4.2 Hasil Model
Hasil Regresi Poisson
library(knitr)
library(kableExtra)
irr_pois %>%
mutate(
estimate = round(estimate, 3),
conf.low = round(conf.low, 3),
conf.high = round(conf.high, 3),
p.value = signif(p.value, 3)
) %>%
rename(
Variabel = term,
IRR = estimate,
`95% CI (Lower)` = conf.low,
`95% CI (Upper)` = conf.high,
`p-value` = p.value
) %>%
kable(
caption = "Hasil Regresi Poisson: Incidence Rate Ratio (IRR)"
) %>%
kable_styling(full_width = FALSE)| Variabel | IRR | std.error | statistic | p-value | 95% CI (Lower) | 95% CI (Upper) |
|---|---|---|---|---|---|---|
| (Intercept) | 9.000000e-03 | 0.0493763 | -94.69139 | 0 | 8.00000e-03 | 1.000000e-02 |
| tx_std_pct | 1.337288e+10 | 0.4433813 | 52.58790 | 0 | 5.60837e+09 | 3.188972e+10 |
| abx_cov | 7.980000e-01 | 0.0044781 | -50.27461 | 0 | 7.91000e-01 | 8.050000e-01 |
Evaluasi Overdispersi
dispersion## [1] 5329.5Hasil Regresi Negative Binomial
irr_nb %>%
mutate(
estimate = round(estimate, 3),
conf.low = round(conf.low, 3),
conf.high = round(conf.high, 3),
p.value = signif(p.value, 3)
) %>%
rename(
Variabel = term,
IRR = estimate,
`95% CI (Lower)` = conf.low,
`95% CI (Upper)` = conf.high,
`p-value` = p.value
) %>%
kable(
caption = "Hasil Regresi Negative Binomial: Incidence Rate Ratio (IRR)"
) %>%
kable_styling(full_width = FALSE)| Variabel | IRR | std.error | statistic | p-value | 95% CI (Lower) | 95% CI (Upper) |
|---|---|---|---|---|---|---|
| (Intercept) | 0.643 | 2.4189718 | -0.1827095 | 0.855 | 0.009 | 1.355820e+02 |
| tx_std_pct | 0.000 | 55.6937369 | -0.3181696 | 0.750 | 0.000 | 3.437266e+42 |
| abx_cov | 1.149 | 0.5587873 | 0.2484925 | 0.804 | 0.361 | 3.760000e+00 |
Perbandingan Model
aic_tbl <- tibble(
Model = c("Poisson", "Negative Binomial"),
AIC = c(AIC(m_pois), AIC(m_nb))
)
aic_tbl %>%
kable(caption = "Perbandingan AIC Model Poisson dan Negative Binomial") %>%
kable_styling(full_width = FALSE)| Model | AIC |
|---|---|
| Poisson | 195890.682 |
| Negative Binomial | 701.195 |
Hasil Disease Mapping : SMR dan EB-Smoothed RR
map_df %>%
st_drop_geometry() %>%
transmute(
Provinsi = prov,
Kasus_Observed = case_tot,
Kasus_Expected = round(expected, 1),
SMR = round(smr, 2),
EB_RR = round(rr_eb, 2)
) %>%
arrange(desc(EB_RR)) %>%
slice(1:10) %>%
kable(
caption = "Top 10 Provinsi Berdasarkan EB-Smoothed Relative Risk"
) %>%
kable_styling(full_width = FALSE)| Provinsi | Kasus_Observed | Kasus_Expected | SMR | EB_RR |
|---|---|---|---|---|
| DI Yogyakarta | 12780 | 4036.9 | 3.17 | 3.16 |
| DKI Jakarta | 41032 | 18144.5 | 2.26 | 2.26 |
| Kalimantan Selatan | 13952 | 8599.3 | 1.62 | 1.62 |
| Papua Selatan | 1873 | 1184.2 | 1.58 | 1.58 |
| Nusa Tenggara Barat | 17883 | 11491.3 | 1.56 | 1.56 |
| Banten | 35808 | 23792.4 | 1.51 | 1.50 |
| Sulawesi Tengah | 8840 | 5995.0 | 1.47 | 1.47 |
| Jawa Timur | 99628 | 71991.0 | 1.38 | 1.38 |
| Jawa Barat | 133926 | 100952.4 | 1.33 | 1.33 |
| Gorontalo | 2843 | 2405.0 | 1.18 | 1.18 |
5 Pembahasan Hasil Analisis
Berikut merupakan hasil analisis epidemiologi pneumonia balita di Indonesia tahun 2024 berdasarkan pendekatan deskriptif, spasial, dan pemodelan statistik, serta pembahasan implikasi epidemiologis dan kebijakan dari temuan tersebut.
5.1 Gambaran Umum Pneumonia Balita Nasional
Berdasarkan hasil analisis, jumlah balita di Indonesia tahun 2024 tercatat sebanyak 26.071.952 anak, dengan total 530.641 kasus pneumonia balita dan 635 kematian. Nilai cumulative incidence (CI) nasional sebesar 20,35 per 1.000 balita, dengan Case Fatality Rate (CFR) nasional 0,12% dan angka kematian sebesar 2,44 per 100.000 balita.
Secara umum, beban pneumonia balita masih menjadi masalah kesehatan masyarakat yang signifikan, terutama mengingat tingginya jumlah kasus absolut dan adanya variasi besar antarprovinsi baik dari sisi insidensi maupun fatalitas.
5.2 Distribusi Insidensi Pneumonia Balita Antarprovinsi
Analisis insidensi menunjukkan variasi spasial yang mencolok antarprovinsi. Sepuluh provinsi dengan insidensi tertinggi (per 1.000 balita) adalah :
- DI Yogyakarta (64,43)
- DKI Jakarta (46,03)
- Kalimantan Selatan (33,02)
- Papua Selatan (32,19)
- Nusa Tenggara Barat (31,67)
- Banten (30,63)
- Sulawesi Tengah (30,01)
- Jawa Timur (28,17)
- Jawa Barat (27,00)
- Gorontalo (24,06)
Sebaliknya, provinsi dengan insidensi terendah antara lain Sulawesi Utara (0,93 per 1.000), Bengkulu, dan Riau. Perbedaan ini mencerminkan heterogenitas risiko pneumonia balita yang kemungkinan dipengaruhi oleh faktor kepadatan penduduk, akses layanan kesehatan, cakupan pengobatan standar, serta kualitas sistem surveilans.
5.3 Pola Case Fatality Rate (CFR) Pneumonia Balita
Distribusi CFR pneumonia balita menunjukkan bahwa provinsi dengan insidensi rendah tidak selalu memiliki CFR rendah. Provinsi dengan CFR tertinggi antara lain :
- Papua Barat Daya (3,63%)
- Sulawesi Utara (2,88%)
- Sulawesi Barat (1,59%)
- Papua Selatan (1,23%)
- Kalimantan Barat (1,10%)
Tingginya CFR di wilayah-wilayah tersebut mengindikasikan kemungkinan keterlambatan diagnosis, keterbatasan fasilitas rujukan, atau akses layanan kesehatan yang belum optimal, khususnya di wilayah timur Indonesia.
5.4 Angka Kematian Pneumonia Balita
Angka kematian per 100.000 balita tertinggi ditemukan di Papua Selatan (39,53) dan Papua Barat Daya (24,67). Meskipun jumlah kasus absolut di wilayah ini relatif kecil, proporsi kematian yang tinggi menunjukkan kerentanan sistem kesehatan terhadap kasus pneumonia berat pada balita.
Sebaliknya, beberapa provinsi dengan insidensi tinggi seperti DKI Jakarta dan Kalimantan Selatan menunjukkan angka kematian mendekati nol, yang mengindikasikan sistem deteksi dan tata laksana kasus yang lebih baik.
5.5 Analisis Spasial Global: Moran’s I
Uji Global Moran’s I terhadap insidensi pneumonia balita menghasilkan nilai Moran’s I sebesar 0,160 dengan p-value 0,141, yang menunjukkan tidak terdapat autokorelasi spasial global yang signifikan.
Hasil ini mengindikasikan bahwa secara nasional, pola insidensi pneumonia balita belum membentuk klaster spasial yang kuat, melainkan lebih bersifat lokal dan dipengaruhi oleh karakteristik spesifik wilayah.
5.6 Analisis Spasial Lokal (LISA)
Analisis Local Indicator of Spatial Association (LISA) mengidentifikasi satu provinsi yang tergolong hotspot (High–High), yaitu Jawa Tengah, yang memiliki insidensi tinggi dan dikelilingi oleh provinsi dengan insidensi tinggi pula.
Sebagian besar provinsi (35 dari 38 provinsi) berada dalam kategori tidak signifikan, sedangkan dua provinsi (Papua Tengah dan Papua Pegunungan) tidak dapat dianalisis karena tidak memiliki kasus dan populasi berisiko.
Temuan ini menunjukkan bahwa klaster pneumonia balita bersifat sangat lokal dan tidak merata di seluruh wilayah Indonesia.
5.7 Standardized Morbidity Ratio (SMR) dan Empirical Bayes Smoothing
Perhitungan Standardized Morbidity Ratio (SMR) menunjukkan nilai SMR berkisar antara 0,05 hingga 3,17, dengan rata-rata nasional sebesar 0,82. Nilai tertinggi ditemukan di DI Yogyakarta (SMR = 3,17) dan DKI Jakarta (SMR = 2,26).
Penerapan Empirical Bayes (EB) smoothing menghasilkan nilai EB-smoothed Relative Risk yang hampir identik dengan SMR pada provinsi berpopulasi besar, namun lebih stabil pada provinsi dengan populasi kecil. Sepuluh provinsi dengan EB-RR tertinggi adalah :
- DI Yogyakarta (3,16)
- DKI Jakarta (2,26)
- Kalimantan Selatan (1,62)
- Papua Selatan (1,58)
- Nusa Tenggara Barat (1,56)
- Banten (1,50)
- Sulawesi Tengah (1,47)
- Jawa Timur (1,38)
- Jawa Barat (1,33)
- Gorontalo (1,18)
Hasil ini menegaskan bahwa EB smoothing efektif untuk mengurangi fluktuasi risiko akibat ukuran populasi kecil tanpa menghilangkan pola risiko utama.
5.8 Pemodelan Regresi Poisson dan Negative Binomial
Model regresi Poisson menunjukkan bahwa variabel cakupan pengobatan standar (tx_std_pct) berhubungan sangat kuat dengan peningkatan angka kejadian pneumonia balita (IRR sangat besar), sementara cakupan antibiotik (abx_cov) berhubungan negatif dengan kejadian.
Namun, nilai dispersion parameter yang sangat tinggi (>5000) menunjukkan adanya overdispersi berat, sehingga model Poisson tidak sesuai.
Model Negative Binomial menunjukkan penurunan drastis nilai AIC (701,2 vs 195.890,7 pada Poisson), mengindikasikan kesesuaian model yang jauh lebih baik. Meski demikian, pada model ini tidak ditemukan hubungan yang signifikan secara statistik antara variabel independen dan kejadian pneumonia balita, yang menunjukkan adanya faktor lain di luar variabel yang dianalisis.
5.9 Implikasi Epidemiologis dan Kebijakan
Hasil analisis menunjukkan bahwa beban pneumonia balita di Indonesia bersifat heterogen secara spasial. Wilayah dengan insidensi tinggi tidak selalu memiliki fatalitas tinggi, dan sebaliknya. Hal ini menegaskan pentingnya pendekatan kebijakan yang spesifik wilayah, bukan kebijakan nasional yang seragam.
Wilayah dengan EB-RR dan CFR tinggi, khususnya di Papua dan Sulawesi, memerlukan prioritas dalam penguatan sistem rujukan, peningkatan kualitas layanan pneumonia berat, serta peningkatan kapasitas tenaga kesehatan.
Sementara itu, wilayah dengan insidensi tinggi namun fatalitas rendah menunjukkan keberhasilan tata laksana klinis yang dapat dijadikan pembelajaran bagi wilayah lain.
5.10 Keterbatasan Analisis
Beberapa keterbatasan dalam analisis ini antara lain penggunaan data agregat provinsi yang berpotensi menyembunyikan variasi intra-provinsi, serta keterbatasan variabel kovariat yang tersedia. Selain itu, wilayah dengan nol populasi berisiko tidak dapat dianalisis secara spasial.
6 Kesimpulan dan Saran
6.1 Kesimpulan
Berdasarkan hasil analisis deskriptif, spasial, dan pemodelan statistik, dapat disimpulkan beberapa temuan utama sebagai berikut :
6.2 6.1 Kesimpulan
Penelitian ini bertujuan untuk menganalisis distribusi dan pola pneumonia balita di Indonesia tahun 2024 secara epidemiologis dan spasial guna mendukung pengambilan keputusan kesehatan masyarakat berbasis wilayah. Berdasarkan hasil analisis yang telah dilakukan, diperoleh beberapa kesimpulan utama sebagai jawaban dari tujuan penelitian yang sudah dibuat, yaitu :
Distribusi pneumonia balita antar provinsi di Indonesia menunjukkan variasi yang sangat lebar.
Secara kuantitatif, insidensi pneumonia balita berbeda signifikan antarprovinsi, dengan nilai tertinggi tercatat di DI Yogyakarta, DKI Jakarta, Kalimantan Selatan, dan Papua Selatan, sementara beberapa provinsi di kawasan timur dan Sumatra menunjukkan insidensi yang relatif lebih rendah. Variasi ini mengindikasikan adanya perbedaan karakteristik populasi, sistem pelaporan, serta kapasitas pelayanan kesehatan antarwilayah.Risiko pneumonia balita antarprovinsi tidak merata dan bersifat heterogen secara spasial.
Analisis risiko menggunakan Standardized Morbidity Ratio (SMR) dan Empirical Bayes (EB) smoothing menunjukkan bahwa sejumlah provinsi memiliki risiko kejadian pneumonia balita yang secara konsisten lebih tinggi dibandingkan rata-rata nasional. Metode EB smoothing memberikan estimasi risiko yang lebih stabil, terutama pada provinsi dengan populasi balita kecil, sehingga memperkuat interpretasi perbedaan risiko antarwilayah.Analisis spasial global tidak menunjukkan autokorelasi spasial yang signifikan, namun analisis lokal mengungkap adanya kluster terbatas.
Uji Moran’s I global menunjukkan tidak adanya pola spasial yang kuat secara nasional, namun analisis Local Indicators of Spatial Association (LISA) mengidentifikasi kluster hotspot terbatas, khususnya di Provinsi Jawa Tengah. Temuan ini menunjukkan bahwa pola pneumonia balita di Indonesia cenderung bersifat lokal dan kontekstual, sehingga intervensi perlu disesuaikan dengan karakteristik wilayah tertentu.Keterkaitan antara indikator pelayanan kesehatan dengan beban pneumonia balita tidak bersifat linier dan sederhana.
Hasil pemodelan regresi menunjukkan bahwa meskipun cakupan pengobatan sesuai standar dan pemberian antibiotik berkaitan dengan variasi kejadian pneumonia balita pada analisis tertentu, hubungan tersebut menjadi tidak signifikan setelah mempertimbangkan overdispersi dan heterogenitas data melalui model Negative Binomial. Hal ini mengindikasikan bahwa beban pneumonia balita dipengaruhi oleh faktor multidimensional di luar indikator pelayanan kesehatan yang dianalisis.Informasi epidemiologis dan spasial yang dihasilkan dapat menjadi dasar penting bagi pengambilan keputusan kesehatan berbasis wilayah.
Pemetaan distribusi, estimasi risiko berbasis EB smoothing, serta identifikasi kluster lokal memberikan gambaran yang lebih komprehensif mengenai beban pneumonia balita di Indonesia. Informasi ini dapat dimanfaatkan oleh pemerintah dan pemangku kepentingan untuk menyusun prioritas intervensi, penguatan sistem pelayanan kesehatan, serta perencanaan program pencegahan yang lebih tepat sasaran.
6.3 Saran
Saran Kebijakan
Penguatan intervensi berbasis wilayah perlu diprioritaskan, khususnya di provinsi dengan insidensi dan EB-smoothed Relative Risk tinggi seperti DI Yogyakarta, DKI Jakarta, Kalimantan Selatan, Papua Selatan, dan Nusa Tenggara Barat.
Wilayah dengan CFR dan angka kematian tinggi, terutama di kawasan Papua dan Sulawesi, memerlukan peningkatan kapasitas layanan kesehatan rujukan, deteksi dini pneumonia berat, serta ketersediaan tenaga kesehatan dan fasilitas penunjang.
Peningkatan kualitas sistem surveilans pneumonia balita diperlukan untuk memastikan konsistensi pelaporan antarprovinsi, terutama di wilayah dengan angka nol kasus atau populasi berisiko nol.
Integrasi pendekatan analisis spasial dan risiko berbasis data dalam perencanaan program kesehatan balita di tingkat provinsi dan nasional dapat membantu penentuan prioritas intervensi secara lebih tepat sasaran.
Saran Penelitian Lanjutan
Penelitian selanjutnya disarankan menggunakan data tingkat kabupaten/kota atau individu untuk menangkap variasi risiko intra-provinsi yang tidak dapat terdeteksi pada analisis agregat provinsi.
Penambahan variabel lain seperti status gizi balita, kepadatan hunian, kualitas udara, tingkat kemiskinan, dan akses sanitasi diharapkan dapat memberikan pemahaman yang lebih komprehensif mengenai determinan pneumonia balita.
Pengembangan model spasial lanjutan, seperti spatial regression atau Bayesian hierarchical models, dapat digunakan untuk mengakomodasi ketergantungan spasial dan heterogenitas wilayah secara lebih optimal.
Studi longitudinal atau analisis tren multi-tahun diperlukan untuk menilai dinamika pneumonia balita dari waktu ke waktu dan mengevaluasi dampak kebijakan kesehatan secara berkelanjutan.
7 Daftar Pustaka
Kementerian Kesehatan Republik Indonesia. (2024). Pneumonia terus
ancam anak-anak (profil kesehatan Indonesia 2024).
United Nations. (2015). Transforming our world: The 2030 Agenda for
Sustainable Development.
UNICEF. (2025). Pneumonia in children statistics.
World Health Organization. (2023). Pneumonia in children.
Kementerian Kesehatan Republik Indonesia. (2024). Profil Kesehatan
Indonesia 2024. Jakarta: Kemenkes RI.
UNICEF. (2025). Pneumonia in children. https://www.unicef.org
World Health Organization. (2023). Pneumonia. https://www.who.int
Gordis, L. (2014). Epidemiology (5th ed.). Philadelphia, PA:
Elsevier Saunders.
Rothman, K. J., Greenland, S., & Lash, T. L. (2008). Modern
epidemiology (3rd ed.). Philadelphia, PA: Lippincott Williams &
Wilkins.
Anselin, L. (1995). Local indicators of spatial association—LISA.
Geographical Analysis, 27(2), 93–115.
Lawson, A. B. (2018). Bayesian disease mapping: Hierarchical
modeling in spatial epidemiology (2nd ed.). Boca Raton, FL: CRC
Press.
Waller, L. A., & Gotway, C. A. (2004). Applied spatial
statistics for public health data. Hoboken, NJ: Wiley.
8 Lampiran
8.1 Lampiran Video Pembelajaran
https://youtube.com/@gustiayuanisaeldina?si=BVvmlCzTvAdLgnv4
8.2 Lampiran Link Dashboard
https://nisadashboard.shinyapps.io/NisaEpidemiologi_DashboardUAS/