ANALISIS EPIDEMIOLOGI KASUS DIARE YANG DILAYANI DI PROVINSI JAWA TENGAH DI TAHUN 2024

Dosen Pembimbing : Dr. I Gede Nyoman Mindra Jaya, M.Si.
Disusun Oleh : Samih Muhamad Alfarras
NPM : 140610230053

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

---

---

1 Pendahuluan

1.1 Latar Belakang

Diare didefinisikan sebagai keluarnya tiga atau lebih tinja cair per hari dan sebagian besar kasus disebabkan oleh infeksi gastrointestinal akibat virus, bakteri, atau parasit protozoa yang ditularkan melalui makanan atau air yang terkontaminasi dan kontak antarindividu (WHO, 2024). Gejala yang khas meliputi tinja cair, kram perut, mual, muntah, demam dan tanda dehidrasi yang pada kasus berat berisiko tinggi menyebabkan malnutrisi dan kematian terutama pada kelompok rentan (CDC, 2024; WHO, 2024). Beban penyakit tetap besar secara global sehingga pemantauan klinis dan pencatatan kasus pada level lokal penting untuk memahami pola morbiditas dan merancang tindakan pencegahan yang efektif (Troeger et al., 2020 dan UNICEF, 2021)

Pola penularan dan beban layanan dapat dipengaruhi oleh mekanisme yang berakar pada sistem pelayanan kesehatan, kondisi sosial ekonomi rumah tangga, karakteristik ruang tempat populasi tinggal, tingkat pendidikan, praktik kebersihan, serta kualitas pasokan air dan sanitasi. Ketersediaan layanan yang memadai menentukan kemampuan deteksi dan penatalaksanaan sehingga memengaruhi angka kasus yang tercatat (Kruk et al., 2018). Kondisi ekonomi dan akses terhadap sumber daya menentukan ekspos terhadap air dan sanitasi yang aman serta status gizi sehingga memengaruhi kerentanan terhadap infeksi enterik (Fink, Günther, dan Hill, 2011). Tingginya kepadatan populasi meningkatkan potensi kontak dengan sumber kontaminasi di ruang bersama, sementara tingkat pendidikan berkaitan dengan praktik pencegahan seperti cuci tangan yang efektif dalam menurunkan risiko diare (Curtis dan Cairncross, 2003).

Analisis pada tingkat agregat kabupaten dan kota menuntut kehati hatian terhadap perbedaan pelaporan antarwilayah dan kemungkinan autokorelasi spasial. Pemeriksaan autokorelasi residual menggunakan statistik seperti Moran I direkomendasikan untuk mendeteksi efek tetangga dan menentukan kebutuhan untuk model spasial (LeSage dan Pace, 2009). Bila autokorelasi spasial signifikan, model spasial seperti spatial lag atau spatial error dan teknik analisis spasial lain dapat membantu menangkap spillover dan mengurangi bias estimasi (LeSage dan Pace, 2009). Normalisasi indikator per populasi, pengecekan kelengkapan data administratif, serta analisis sensitivitas dengan ukuran alternatif untuk indikator lingkungan dan layanan akan memperkuat validitas temuan (WHO/UNICEF JMP, 2021).

Sebagai alternatif bila tidak terdapat autokorelasi spasial yang signifikan, analisis regresi untuk data count merupakan pendekatan yang wajar. Poisson generalized linear model adalah titik awal yang umum untuk memodelkan jumlah kejadian diskret karena model ini memodelkan ekspektasi jumlah kejadian dengan link log dan memungkinkan penggunaan offset berupa log populasi untuk memodelkan rate per populasi (McCullagh dan Nelder, 1989). Jika data menunjukkan overdispersi, yaitu varians yang lebih besar dari mean, pendekatan quasi-Poisson dapat memperbaiki estimasi varians sehingga inferensi menjadi lebih andal (Cameron dan Trivedi, 1998). Bila overdispersi nyata dan berasal dari heterogenitas ekstra, negative binomial regression memperkenalkan parameter dispersi tambahan sehingga memodelkan varians yang melebihi mean, dan model ini sering lebih realistis untuk data kasus penyakit infeksi (Hilbe, 2011). Pemilihan model hendaknya didasarkan pada diagnostik seperti uji overdispersi, pemeriksaan residual, serta perbandingan kecocokan model.

1.2 Rumusan Masalah

Berdasarkan latar tersebut, beberapa masalah utama yang diidentifikasi adalah:

1. Berapa angka cumulative incidence kasus diare di Provinsi Jawa Tengah di tahun 2024?

2. Bagaimana distribusi cumulative incidence diare berdasarkan Kabupaten/Kota di Provinsi Jawa Tengah?

3. Apabila autokorelasi spasial signifikan, apakah Spatial Error Model mampu menggambarkan hubungan kausal langsung dan tidak langsung antara kondisi sosial ekonomi, karakteristik lingkungan permukiman, ketersediaan layanan kesehatan, dan kondisi air serta sanitasi terhadap jumlah kasus diare yang dilayani pada tingkat kabupaten/kota di Provinsi Jawa Tengah tahun 2024?

4. Apabila autokorelasi spasial tidak signifikan, apakah hubungan kausal tersebut dapat dimodelkan secara tepat menggunakan pendekatan non spasial untuk data count, antara lain Poisson GLM, quasi-Poisson, atau negative binomial, dengan mempertimbangkan penggunaan offset populasi dan pemeriksaan overdispersi?

5. Bagaimana model simulasi desain studi kasus-kontrol dapat dibangun untuk mengukur asosiasi antara Kabupaten/Kota dengan cumulative incidence diare di Jawa Tengah?

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk:

1. Mendeskripsikan pola sebaran spasial kasus diare di Provinsi Jawa Tengah di tahun 2024 pada tingkat kabupaten/kota.

2. Menghitung angka cumulative incidence diare di Provinsi Jawa Tengah di tahun 2024.

3. Mengidentifikasi distribusi cumulative incidence diare berdasarkan Kabupaten/Kota di Provinsi Jawa Tengah.

4. Jika autokorelasi spasial signifikan, menganalisis hubungan kausal langsung dan tidak langsung antara konstruk sosial ekonomi, karakteristik lingkungan permukiman, ketersediaan layanan kesehatan, dan kondisi air serta sanitasi terhadap jumlah kasus diare yang dilayani pada tingkat kabupaten/kota tahun 2024 menggunakan Spatial Error Model.

5. Jika autokorelasi spasial tidak signifikan, menganalisis hubungan kausal tersebut menggunakan pendekatan non spasial yang sesuai untuk data count, antara lain Poisson GLM, quasi-Poisson, dan negative binomial, dengan mempertimbangkan offset populasi dan pemeriksaan overdispersi.

6. Mengetahui gambaran cumulative incidence kejadian diare di Provinsi Jawa Tengah menggunakan desain studi kasus-kontrol.

2 Tinjauan Pustaka

2.1 Epidemiologic Triad (Agent-Host-Environment)

Diare sebagai entitas klinis dan epidemiologis paling baik dipahami melalui kerangka triad yang terdiri dari agen penyebab penyakit, sifat host yang rentan, dan kondisi lingkungan yang memungkinkan transmisi. Kerangka ini menegaskan bahwa kejadian penyakit bukan hasil satu faktor tunggal melainkan interaksi dinamis antar komponen triad yang menentukan insiden, keparahan dan pola spasial kasus (Keusch et al., 2016). Secara global diare tetap menyumbang beban morbiditas dan mortalitas yang substansial sehingga pendekatan yang mempertimbangkan ketiga komponen triad sangat penting untuk desain intervensi pencegahan dan respons klinis. (Troeger et al., 2020; WHO, 2024).

2.1.1 Agent

Komponen agen meliputi beragam patogen enterik termasuk virus seperti rotavirus dan norovirus, bakteri seperti enterotoksigenik Escherichia coli dan Shigella, serta protozoa seperti Giardia dan Cryptosporidium. Masing masing agen berbeda dalam hal dosis infeksi minimal, virulensi, lama penularan, dan ketahanan di media lingkungan sehingga karakteristik agen memengaruhi pola wabah dan strategi kontrol seperti vaksinasi, pengobatan dan langkah keselamatan pangan (Akhondi, 2023; CDC, 2024). Keberadaan agen tertentu juga menentukan presentasi klinis yang lebih condong ke diare berat atau diare berdarah sehingga penting bagi sistem surveilans untuk mengidentifikasi etiologi guna mengarahkan tindakan penanggulangan yang tepat.

2.1.2 Host

Sisi host menekan faktor usia, status imunisasi, status nutrisi dan kondisi komorbid. Anak usia di bawah lima tahun dan lansia menunjukkan kerentanan tinggi terhadap keparahan dan komplikasi dehidrasi sebagaimana dicatat dalam literatur klasik tentang penyakit diare (Keusch et al., 2016). Malnutrisi menurunkan daya tahan mukosa gastrointestinal dan meningkatkan risiko infeksi berulang serta komplikasi yang lebih berat (Troeger et al., 2020). Imunitas spesifik, misalnya akibat vaksinasi rotavirus, menurunkan kejadian dan keparahan penyakit yang disebabkan oleh agen tertentu sehingga program imunisasi merupakan intervensi pencegahan primer yang efektif pada populasi rentan (Kyu et al., 2024).

2.1.3 Environment

Lingkungan berperan sebagai mediator utama ekspos terhadap agen. Akses terhadap air minum yang aman, sanitasi yang layak dan praktik kebersihan personal menentukan peluang paparan terhadap patogen fecal—oral serta memengaruhi besarnya penularan komunitas (WHO/UNICEF JMP, 2021). Bukti empiris menunjukkan bahwa perbaikan sanitasi dan promosi cuci tangan mampu menurunkan risiko diare secara signifikan sehingga intervensi WASH menjadi komponen sentral pengendalian penyakit ini (Fink, Günther, & Hill, 2011; Curtis & Cairncross, 2003). Selain itu faktor fisik lain seperti kepadatan populasi dan kualitas layanan kesehatan memengaruhi bagaimana penyakit menyebar dan dilaporkan sehingga penilaian risiko harus mempertimbangkan kondisi lingkungan lokal dan infrastruktur publik. (WHO/UNICEF JMP, 2021; Fink et al., 2011).

Interaksi ketiga komponen inilah yang menentukan tingkat endemisitas Diare di suatu wilayah.

2.2 Desain Studi Epidemiologi

Desain studi adalah cetak biru atau arsitektur penelitian yang memandu peneliti dalam mengumpulkan dan menganalisis data untuk menjawab pertanyaan penelitian (Budiarto, 2012). Secara garis besar, desain studi dibagi menjadi deskriptif (mendeskripsikan distribusi penyakit) dan analitik (menguji hipotesis hubungan sebab-akibat).Desain analitik observasional yang umum digunakan meliputi:

2.2.1 Studi Cross-sectional

Desain studi cross-sectional atau potong lintang adalah sebuah “potret” dalam satu titik waktu, di mana peneliti mengukur variabel paparan (faktor risiko) dan luaran (penyakit) secara bersamaan pada sekelompok individu (Celentano & Szklo, 2018). Karena semua data dikumpulkan sekaligus, desain ini tidak memiliki arah waktu sehingga tidak dapat menentukan apakah paparan terjadi sebelum penyakit atau sebaliknya. Oleh karena itu, studi cross-sectional sangat ideal untuk mengukur prevalensi, yaitu proporsi penyakit yang ada di populasi pada saat tertentu, dan sering digunakan sebagai langkah awal untuk menghasilkan hipotesis, namun tidak kuat untuk menyimpulkan hubungan sebab-akibat (Rothman et al., 2008).

2.2.2 Studi Kohort (Cohort)

Desain studi kohort bersifat prospektif atau melihat ke depan. Penelitian dimulai dengan mengidentifikasi dua kelompok individu yang sehat atau belum mengalami luaran, yaitu kelompok yang terpapar faktor risiko dan kelompok yang tidak terpapar. Kedua kelompok ini kemudian diikuti (follow-up) selama periode tertentu untuk membandingkan kejadian penyakit di antara keduanya (Celentano & Szklo, 2018). Karena paparan dipastikan terjadi sebelum luaran, studi kohort merupakan desain observasional yang kuat untuk mengevaluasi hubungan sebab-akibat dan memungkinkan pengukuran insidensi serta Relative Risk (RR) (Hulley et al., 2013). Namun demikian, desain ini umumnya membutuhkan biaya besar dan waktu penelitian yang relatif lama.

2.2.3 Studi Kasus-Kontrol (Case-Control)

Desain studi kasus-kontrol bekerja secara retrospektif atau melihat ke belakang. Peneliti memulai dengan mengidentifikasi individu yang telah mengalami penyakit sebagai kelompok kasus dan individu yang tidak mengalami penyakit sebagai kelompok kontrol (Rothman et al., 2008). Selanjutnya, riwayat paparan faktor risiko pada kedua kelompok tersebut ditelusuri, biasanya melalui wawancara atau data rekam medis. Studi kasus-kontrol sangat efisien untuk meneliti penyakit yang jarang terjadi atau penyakit dengan periode laten yang panjang (Celentano & Szklo, 2018). Ukuran asosiasi utama yang digunakan dalam desain ini adalah Odds Ratio (OR), yang menggambarkan perbandingan peluang paparan antara kelompok kasus dan kelompok kontrol.

2.3 Autokorelasi Spasial

Autokorelasi spasial digunakan untuk mengukur sejauh mana nilai suatu variabel di satu lokasi berkorelasi dengan nilai di lokasi lain yang berdekatan. Terdapat dua jenis analisis autokorelasi spasial, yaitu Global Moran’s I dan Local Moran’s I (LISA). Global Moran’s I digunakan untuk menilai apakah pola spasial secara keseluruhan bersifat mengelompok (clustered), acak (random), atau menyebar (dispersed), sedangkan Local Moran’s I digunakan untuk mengidentifikasi lokasi-lokasi spesifik yang membentuk klaster signifikan, seperti wilayah dengan kasus TBC tinggi yang berdekatan dengan wilayah serupa (High-High) atau rendah dengan rendah (Low-Low). Kedua ukuran ini penting untuk memahami pola sebaran penyakit dan mendeteksi daerah prioritas penanganan

2.4 Spatial Error Model (SEM)

Model SEM memodelkan ketergantungan spasial pada komponen error, dengan bentuk:

\[ y = X\beta + u, \qquad u = \lambda W u + \varepsilon \]

Di mana :

\(y\) adalah vektor variabel dependen,

\(X\) adalah matriks kovariat,

\(\beta\) adalah vektor koefisien regresi,

\(W\) adalah matriks pembobot spasial,

\(\lambda\) adalah parameter error spasial, dan

\(\varepsilon\) adalah komponen error yang diasumsikan white noise.

Model SEM memperhitungkan korelasi spasial pada komponen error, yang menunjukkan adanya faktor spasial yang tidak teramati tetapi memengaruhi variabel dependen. Dengan membandingkan kedua model ini, dapat diketahui apakah pengaruh spasial lebih dominan berasal dari interaksi antarwilayah atau dari faktor-faktor tak teramati dalam struktur residual.

2.5 Overdispersi

Overdispersi terjadi ketika varians observasi pada data count lebih besar dari mean yang diasumsikan oleh model Poisson. Dua ukuran diagnostik umum untuk mendeteksi overdispersi adalah rasio deviance per derajat kebebasan dan rasio Pearson chi square per derajat kebebasan. Bila salah satu rasio jauh lebih besar dari 1 maka overdispersi dianggap ada.

Rumus rasio deviance per derajat kebebasan:

\[ \hat{\phi}_{deviance} ;=; \frac{\text{Deviance}}{n - p} \]

Rumus rasio Pearson chi square per derajat kebebasan:

\[ \hat{\phi}*{Pearson} ;=; \frac{\sum*{i=1}^{n} \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i}}{n - p} \]

Keterangan y_i adalah observasi ke i, _i adalah prediksi mean dari model Poisson, n adalah jumlah observasi, p adalah jumlah parameter yang diestimasi termasuk intercept.

Interpretasi Jika ( ) maka asumsi Poisson tentang equidispersion terpenuhi. Jika ( ) maka terdapat overdispersi dan perlu pertimbangan model alternatif.

2.6 Poisson GLM Regression

Poisson GLM memodelkan variabel respons count (Y_i) dengan asumsi distribusi Poisson dan link log. Model ini sering menjadi titik awal untuk analisis data jumlah kasus. Koefisien regresi pada skala log dapat diubah menjadi incidence rate ratio IRR dengan eksponensial.

Model dasar:

\[ Y_i \sim \text{Poisson}(\mu_i) \]

Fungsi link log:

\[ \log(\mu_i) ;=; X_i \beta \]

Dengan offset populasi untuk memodelkan rate per populasi:

\[ \log(\mu_i) ;=; X_i \beta ;+; \log(\text{pop}_i) \]

Varian Poisson:

\[ \text{Var}(Y_i) ;=; \mu_i \]

Interpretasi koefisien Koefisien (_j) diinterpretasikan pada skala log. Untuk interpretasi praktis gunakan IRR:

\[ \text{IRR}_j ;=; e^{\beta_j} \]

Poisson digunakan bila equidispersion terpenuhi. Bila overdispersi terdeteksi, Poisson akan memberikan standar error yang diremehkan sehingga inferensi menjadi tidak andal.

2.7 Quasi-Poisson Regression

Quasi-Poisson mempertahankan bentuk mean Poisson tetapi memperkenalkan faktor dispersi () untuk mengakomodasi varians yang lebih besar. Model ini memperbaiki estimasi varians dan standar error tanpa merombak struktur mean.

Struktur mean tetap:

\[ \log(\mu_i) ;=; X_i \beta \]

Fungsi varian pada quasi-Poisson:

\[ \text{Var}(Y_i) ;=; \phi \mu_i \]

Estimasi dispersi phi secara praktis dapat dihitung dari Pearson chi square atau deviance seperti rumus pada bagian Overdispersi. Karena family quasi tidak memiliki likelihood penuh, kriteria berbasis likelihood seperti AIC tidak berlaku untuk membandingkan model quasi dengan model berbasis likelihood.

Interpretasi koefisien Interpretasi koefisien sama dengan Poisson sehingga IRR tetap dihitung sebagai (e^{_j}). Kelebihan quasi-Poisson adalah koreksi inferensial yang sederhana ketika overdispersi bersifat global dan tidak perlu dimodelkan secara eksplisit.

2.8 Negative Binomial Regression

Negative binomial memperluas Poisson dengan memasukkan parameter dispersi yang memodelkan heterogenitas ekstra antarunit. Salah satu parameterisasi umum menyatakan varian sebagai fungsi kuadrat mean sehingga varians dapat meningkat lebih cepat daripada mean. Karena memiliki formulasi likelihood, negative binomial memungkinkan penggunaan AIC untuk perbandingan model.

Model umum NB2:

\[ Y_i \sim \text{NegBin}(\mu_i, \alpha) \]

Mean dan varian:

\[ \mathbb{E}[Y_i] ;=; \mu_i \]

\[ \text{Var}(Y_i) ;=; \mu_i ;+; \alpha \mu_i^2 \]

Parameter dispersi alternatif sering dilambangkan dengan () sehingga varian dapat ditulis:

\[ \text{Var}(Y_i) ;=; \mu_i ;+; \frac{\mu_i^2}{\theta} \]

Skema link log untuk mean:

\[ \log(\mu_i) ;=; X_i \beta \]

Interpretasi koefisien Interpretasi koefisien pada NB sama dengan Poisson sehingga IRR diperoleh dengan (e^{_j}). Parameter () atau () mengukur besarnya heterogenitas ekstra. Jika () atau () model NB mendekati Poisson.

Negative binomial digunakan bila overdispersi nyata muncul karena heterogenitas laten antarunit yang tidak terukur oleh kovariat. Karena NB memodelkan sumber variabilitas ekstra secara eksplisit, model ini biasanya memberikan fit dan inferensi yang lebih baik dibandingkan Poisson atau quasi-Poisson bila overdispersi bersifat struktural.

2.9 Ukuran Asosiasi dan Frekuensi

Dalam epidemiologi, ukuran digunakan untuk mengkuantifikasi frekuensi penyakit dan kekuatan hubungan antara paparan dan penyakit.

2.9.1 Prevalensi

Prevalensi adalah ukuran proporsi individu dalam suatu populasi yang memiliki penyakit (kasus existing, baik baru maupun lama) pada satu titik waktu tertentu (point prevalence) atau selama periode waktu tertentu (period prevalence).

\[ \text{Prevalensi} = \frac{\text{Jumlah kasus (baru + lama) pada waktu } t}{\text{Total populasi pada waktu } t} \]

Ukuran ini menggambarkan beban penyakit (burden of disease) di masyarakat.

2.9.2 Insidensi

Insidensi mengukur jumlah kasus baru suatu penyakit yang terjadi dalam populasi berisiko selama periode waktu tertentu. Insidensi mencerminkan risiko seseorang untuk terkena penyakit. Ini biasanya diukur dalam studi kohort. Rumus untuk Cumulative Incidence (CI) atau Risk adalah:

\[ \text{Insidensi Kumulatif (Risk)} = \frac{\text{Jumlah kasus baru selama periode waktu tertentu}} {\text{Jumlah populasi berisiko}} \]

2.9.3 Attack rate (AR)

Attack Rate adalah ukuran frekuensi yang digunakan dalam epidemiologi untuk menggambarkan proporsi individu dalam suatu populasi yang menjadi sakit (terinfeksi) selama periode waktu tertentu.

\[ AR = \frac{\text{Jumlah kasus baru selama periode wabah}} {\text{Jumlah populasi yang berisiko}} \times 100\% \]

2.9.4 Odds Ratio (OR)

OR adalah ukuran asosiasi yang digunakan dalam studi kasus-kontrol dan cross-sectional. OR membandingkan odds (peluang) terjadinya penyakit pada kelompok terpapar dengan odds terjadinya penyakit pada kelompok tidak terpapar.

• OR = 1 → tidak ada asosiasi antara paparan dan penyakit.

• OR > 1 → paparan berhubungan dengan peningkatan odds penyakit (faktor risiko).

• OR < 1 → paparan berhubungan dengan penurunan odds penyakit (faktor protektif).

\[ OR = \frac{(a/b)}{(c/d)} = \frac{ad}{bc} \] ### Risk Ratio (RR) atau Relative Risk

RR adalah ukuran asosiasi yang digunakan dalam studi kohort. RRmembandingkan risiko (insidensi) penyakit pada kelompok terpapar dengan risiko (insidensi) penyakit pada kelompok tidak terpapar.

\[ RR = \frac{\text{Insidensi pada kelompok terpapar}} {\text{Insidensi pada kelompok tidak terpapar}} = \frac{a / (a + b)}{c / (c + d)} \]

2.9.5 Attributable Risk (AR)

AR (atau Selisih Risiko) adalah ukuran yang menunjukkan kelebihan risikopenyakit pada kelompok terpapar yang dapat diatribusikan (disebabkan)oleh paparan tersebut. Ukuran ini juga memerlukan data insidensi daristudi kohort.

\[ AR = (\text{Insidensi terpapar}) - (\text{Insidensi tidak terpapar}) \]

3 Metodologi Penelitian

3.1 Sumber Data

Data yang digunakan dalam penelitian ini merupakan data sekunder yang diperoleh dari sumber resmi yaitu Portal Data Jawa Tengah dan Badan Pusat Statistik (BPS).

3.2 Variabel

Variabel yang digunakan adalah jumlah kasus Diare yang dilayani di setiap kabupaten/kota di Provinsi Jawa Tengah sebagai variabel dependen, jumlah penduduk di setiap Kabupaten/Kota di Provinsi Jawa Tengah pada tahun 2024 sebagai variabel offset, dan untuk variabel independennya adalah jumlah fasilitas kesehatan, persentase penduduk miskin, kepadatan penduduk, rata-rata lama sekolah, persentase air minum layak, persentase sanitasi layak.

3.3 Metode Analisis

Analisis yang digunakan terdiri atas perhitungan cumulative incidence dan autokorelasi spasial. Analisis deskriptif dilakukan untuk memvisualisasikan persebaran kasus diare di Jawa Tengah. Autokorelasi spasial dihitung menggunakan Global Moran’s I dan Geary’s C untuk mendeteksi pola sebaran secara keseluruhan dan Local Moran’s I (LISA) untuk mengidentifikasi klaster signifikan antarwilayah. Kemudian dilanjut dengan pemodelan SEM untuk melihat hubungan antar variabel, jika tidak terdapat autokorelasi spasial yang signifikan akan dideteksi overdispersi dan menentukan model regresi yang digunakan sesuai dengan hasil pengecekan overdispersi. Yang terakhir akan diberikan studi kasus desain eksperimen Case-Control.

3.4 Alur Kerja Penelitian

Penelitian dimulai dengan pengumpulan data-data sekunder, yaitu populasi total kabupaten serta kota di Jawa Tengah dari 2024 kemudian total kasus penyakit diare yang akan diikuti oleh analisis deskriptif berdasarkan epidemiologi. Analisis akan dilakukan pada program lunak Rstudio. Pertama akan dilakukan plottingan peta menggunakan jumlah total kasus diare kemudian melakukan analisis deskriptif. Kemudian menghitung Cumulative Incidence berdasarkan populasi pada semua kabupaten/kota pada tahun 2024 dimana kemudian akan dilakukan plottingan pada peta kemudian analisis spasial autokorelasi pada cumulative incidence. Kemudian akan dilakukan pemodelan SEM jika terdapat autokorelasi spasial yang signifikan. Jika tidak terdapat autokorelasi spasial yang signifikan maka akan dideteksi overdispersi dan ditentukan model regresi yang digunakan sesuai dengan hasil pengecekan overdispersi. dan diakhiri dengan studi kasus desain eksperimen.

4 Hasil dan Pembahasan

4.1 Tabel Hasil

Berdasarkan perhitungan yang dilakukan didapatkan hasil tabel berikut yang mencangkup data populasi serta jumlah total kasus diare yang dilayani pada masing-masing kabupaten. Dilakukan perhitungan Cumulative Incidence pada masing-masing kabupaten serta kota untuk melihat tingkat keparahan atau beban diare yang diberikan pada populasi di masing-masing kabupaten/kota dimana didapatkan tabel sebagai berikut:

## Reading layer `gadm41_IDN_2' from data source 
##   `E:\Backup Window\Documents\Tugas Kuliah\Semester 5\Spasial\gadm41_IDN_shp\gadm41_IDN_2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS:  WGS 84

## [1] "Jumlah Kab/Kota di shapefile: 36"

## [1] "Daftar Kab/Kota di shapefile:"

##  [1] "Banjarnegara"     "Banyumas"         "Batang"           "Blora"           
##  [5] "Boyolali"         "Brebes"           "Cilacap"          "Demak"           
##  [9] "Grobogan"         "Jepara"           "Karanganyar"      "Kebumen"         
## [13] "Kendal"           "Klaten"           "Kota Magelang"    "Kota Pekalongan" 
## [17] "Kota Semarang"    "Kota Tegal"       "Kudus"            "Magelang"        
## [21] "Pati"             "Pekalongan"       "Pemalang"         "Purbalingga"     
## [25] "Purworejo"        "Rembang"          "Salatiga"         "Semarang"        
## [29] "Sragen"           "Sukoharjo"        "Surakarta"        "Tegal"           
## [33] "Temanggung"       "Waduk Kedungombo" "Wonogiri"         "Wonosobo"

## [1] "Jumlah data Diare: 36"

## [1] "Daftar Kab/Kota di data Diare:"

##  [1] ""             "Banjarnegara" "Banyumas"     "Batang"       "Blora"       
##  [6] "Boyolali"     "Brebes"       "Cilacap"      "Demak"        "Grobogan"    
## [11] "Jepara"       "Karanganyar"  "Kebumen"      "Kendal"       "Klaten"      
## [16] "Kudus"        "Magelang"     "Pati"         "Pekalongan"   "Pemalang"    
## [21] "Purbalingga"  "Purworejo"    "Rembang"      "Salatiga"     "Semarang"    
## [26] "Sragen"       "Sukoharjo"    "Surakarta"    "Tegal"        "Temanggung"  
## [31] "Wonogiri"     "Wonosobo"

## [1] "Kab/Kota yang ada di shapefile tapi TIDAK ada di data:"

## [1] "Kota Magelang"    "Kota Pekalongan"  "Kota Semarang"    "Kota Tegal"      
## [5] "Waduk Kedungombo"

## [1] "Kab/Kota yang ada di data tapi TIDAK ada di shapefile:"

## [1] ""

## [1] "Jumlah Kab/Kota setelah filter: 31"

## [1] "Jumlah data yang berhasil di-join: 35"

## [1] "Jumlah data yang gagal di-join: 0"

## [1] "TABEL DATA DIARE JAWA TENGAH 2024"

##         KabKota     Y Populasi    CI_2024
## 1  Banjarnegara 10319  1047226  985.36515
## 2      Banyumas 19095  1828573 1044.25692
## 3        Batang  4538   828883  547.48378
## 4         Blora  6117   901621  678.44471
## 5      Boyolali  8186  1090129  750.92030
## 6        Brebes  7728  2043077  378.25300
## 7       Cilacap 13524  2007829  673.56334
## 8         Demak 13353  1240510 1076.41212
## 9      Grobogan  1112  1492891   74.48635
## 10       Jepara 10352  1221086  847.76994
## 11  Karanganyar 20446   955116 2140.68239
## 12      Kebumen 22189  1397555 1587.70138
## 13       Kendal  3694  1052826  350.86520
## 14       Klaten 17737  1284386 1380.97114
## 15        Kudus  6501   874632  743.28403
## 16     Magelang 13994  1330656 1051.66174
## 17     Magelang  4070   122150 3331.96889
## 18         Pati  5830  1359364  428.87703
## 19   Pekalongan  3576  1007384  354.97884
## 20   Pekalongan  9122   317524 2872.85371
## 21     Pemalang 10847  1523622  711.92199
## 22  Purbalingga 19049  1027333 1854.21864
## 23    Purworejo  4516   788265  572.90378
## 24      Rembang  7022   660166 1063.67186
## 25     Salatiga  1218   198920  612.30645
## 26     Semarang 11015  1080648 1019.29583
## 27     Semarang 36516  1694743 2154.66298
## 28       Sragen  8778   997485  880.01323
## 29    Sukoharjo 11127   932680 1193.01368
## 30    Surakarta  8998   526870 1707.82166
## 31        Tegal 29230  1654836 1766.33817
## 32        Tegal  4807   282781 1699.90204
## 33   Temanggung  6883   808446  851.38649
## 34     Wonogiri  3897  1051085  370.75974
## 35     Wonosobo  9458   909664 1039.72456

## [1] "STATISTIK DESKRIPTIF"

##        Y            Populasi          CI_2024       
##  Min.   : 1112   Min.   : 122150   Min.   :  74.49  
##  1st Qu.: 5318   1st Qu.: 851758   1st Qu.: 642.93  
##  Median : 8998   Median :1047226   Median : 985.37  
##  Mean   :10710   Mean   :1072599   Mean   :1108.54  
##  3rd Qu.:13438   3rd Qu.:1345010   3rd Qu.:1484.34  
##  Max.   :36516   Max.   :2043077   Max.   :3331.97

4.2 Peta Deskriptif

Pertama dilakukan pembuatan peta berdasarkan dari jumlah total kasus Diare yang terdapat pada masing-masing kabupaten/kota yang ada di Jawa Tengah. Berikut didapatkan peta Total Kasus Diare di Jawa Tengah.

library(ggplot2)
library(scales)
library(sf)

# Fungsi membuat label interval
interval_labels <- function(breaks) {
  labs <- c()
  for (i in seq_along(breaks)[-length(breaks)]) {
    labs[i] <- paste0(
      format(breaks[i], big.mark = "."),  # Gunakan titik untuk pemisah ribuan (Indonesia)
      " – ",
      format(breaks[i + 1], big.mark = ".")
    )
  }
  labs
}


# Tentukan breaks untuk jumlah kasus
max_kasus <- max(jateng_merged$Y, na.rm = TRUE)
min_kasus <- min(jateng_merged$Y, na.rm = TRUE)

# Sesuaikan interval berdasarkan range data
# Jika data 3.000-40.000, bisa gunakan interval 20.000
interval_size <- 20000  # Sesuaikan sesuai kebutuhan
breaks_kasus <- seq(0, ceiling(max_kasus / interval_size) * interval_size, by = interval_size)

# Ambil label
label_kasus <- interval_labels(breaks_kasus)

# Plot Jumlah Kasus Diare 2024
plot_kasus_diare <- ggplot(jateng_merged) +
  geom_sf(aes(fill = Y),
          color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(
    option = "plasma",
    direction = -1,
    name = "Jumlah Kasus Diare",
    breaks = breaks_kasus[-1],  # buang elemen pertama
    labels = label_kasus
  ) +
  labs(
    title = "Jumlah Kasus Diare di Jawa Tengah 2024",
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_line(color = "gray95"),
    plot.title = element_text(face = "bold", hjust = 0.5),
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 9),
    legend.position = "right"
  )

print(plot_kasus_diare)

Distribusi spasial ini memperlihatkan ketimpangan yang cukup jelas antara wilayah tengah serta selatan provinsi dibandingkan dengan wilayah timur dan utara. Secara spesifik Kabupaten Boyolali yang terletak di bagian tengah provinsi menjadi wilayah yang paling mencolok dengan warna biru tua pekat yang mengindikasikan bahwa kabupaten ini memiliki beban kasus diare tertinggi di seluruh Jawa Tengah pada tahun tersebut. Selain itu konsentrasi kasus yang tinggi juga terlihat di wilayah barat daya dan pesisir selatan yang meliputi Kabupaten Banyumas serta Kabupaten Cilacap dan Kabupaten Kebumen yang ditandai dengan warna ungu hingga merah tua. Sebaliknya wilayah di sisi timur seperti Kabupaten Blora dan Kabupaten Rembang serta Kabupaten Pati didominasi oleh warna kuning dan jingga muda yang menunjukkan bahwa angka kejadian diare di daerah-daerah tersebut relatif jauh lebih rendah dan lebih terkendali dibandingkan pusat penyebaran di Boyolali maupun wilayah selatan. Untuk melihat lebih dalam akan dilakukan analisis deskriptif mengenai jumlah total kasus Diare:

library(psych)
library(dplyr)
library(sf)

describe(
  jateng_merged %>%
    st_drop_geometry() %>%         # buang kolom geometry dulu
    dplyr::select(Y)
)

##   vars  n     mean      sd median trimmed     mad  min   max range skew
## Y    1 35 10709.83 7768.75   8998 9689.76 6456.72 1112 36516 35404 1.42
##   kurtosis      se
## Y     1.95 1313.16

Dari deskriptif bisa dilihat karakteristik pada data jumlah total kasus Diare di Jawa Tengah. Didapatkan juga visualisasi histogram jumlah kasus Diare di Jawa Tengah :

# FUNGSI HISTOGRAM UNTUK DATA DIARE

plot_histogram_diare <- function(data, var, judul, warna, bins = 10) {
  # Ambil data dari kolom yang ditentukan
  x <- data[[var]]
  
  # Hitung breaks manual
  breaks <- pretty(range(x, na.rm = TRUE), n = bins)
  
  ggplot(data, aes(x = x)) +
    geom_histogram(
      breaks = breaks,
      fill = warna,
      color = "white",
      alpha = 0.8
    ) +
    scale_x_continuous(
      breaks = breaks,     # label batas bawah & atas
      labels = comma       # format ribuan
    ) +
    labs(
      title = judul,
      x = "Batas Interval",
      y = "Frekuensi"
    ) +
    theme_minimal(base_size = 13) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(face = "bold", hjust = 0.5)
    )
}

# HISTOGRAM UNTUK VARIABEL Y (KASUS DIARE)

plot_histogram_diare(
  jateng_merged, 
  "Y", 
  "Distribusi Jumlah Kasus Diare Jawa Tengah 2024", 
  "#F46D43",
  bins = 8  # Sesuaikan jumlah bins
)

Histogram tersebut menunjukkan pola distribusi menjulur ke kanan atau positively skewed. Mayoritas data terkonsentrasi di sisi kiri grafik yang mengindikasikan bahwa sebagian besar wilayah di Jawa Tengah memiliki jumlah kasus yang relatif rendah hingga sedang. Puncak distribusi atau modus berada pada interval 5.000 hingga 10.000 kasus dengan frekuensi tertinggi yang mencakup lebih dari 10 wilayah administratif.

Pola penurunan frekuensi terlihat cukup konsisten seiring dengan bertambahnya jumlah kasus. Wilayah dengan jumlah kasus antara 0 hingga 5.000 dan 10.000 hingga 15.000 masih memiliki frekuensi yang cukup banyak yakni masing-masing di kisaran 8 hingga 9 wilayah. Namun jumlah wilayah menurun drastis pada interval di atas 15.000 kasus di mana hanya terdapat sedikit kabupaten atau kota yang memiliki beban kasus setinggi itu. Hal ini menandakan bahwa kejadian luar biasa atau beban penyakit yang sangat tinggi bukan merupakan fenomena umum di seluruh provinsi melainkan hanya terjadi di beberapa titik tertentu.

Keberadaan pencilan atau outlier terlihat sangat jelas di ujung kanan grafik. Terdapat kekosongan data atau gap pada interval 30.000 hingga 35.000 kasus tetapi kemudian muncul satu wilayah yang berada pada interval ekstrem yakni antara 35.000 hingga 40.000 kasus. Nilai ekstrem ini kemungkinan besar berkorespondensi dengan wilayah yang diidentifikasi berwarna biru pekat pada peta sebelumnya seperti Kabupaten Boyolali yang memiliki jumlah kasus jauh melampaui rata-rata wilayah lainnya di Jawa Tengah.

Salah satu tantangan yang dihadapi saat melakukan analisis adalah heterogenitas populasi, dimana kasus yang banyak merupakan indikasi bahwa karena ada lebih banyak populasi sehingga akan lebih banyak kasus. Karena itu, digunakannya prevalence rate untuk menangani hal tersebut dalam bentuk persentase. Didapatkan plot map seperti berikut:

library(ggplot2)
library(sf)
library(dplyr)

# Hitung CI per 100.000 penduduk
jateng_merged <- jateng_merged %>%
  mutate(CI_2024 = (Y / Populasi) * 100000)

# PERBAIKAN 1: Gunakan quantile breaks untuk distribusi lebih merata
breaks_interval <- quantile(
  jateng_merged$CI_2024, 
  probs = seq(0, 1, by = 0.2),  # 5 kategori
  na.rm = TRUE
)

# Buat labels yang sesuai dengan breaks
label_interval <- format(round(breaks_interval), big.mark = ".")

# PERBAIKAN 2: Fungsi plot dengan skala lebih jelas
plot_ci_interval <- function(data, var, year) {
  ggplot(data) +
    geom_sf(aes(fill = .data[[var]]), color = "white", linewidth = 0.5) +
    scale_fill_viridis_c(
      option = "plasma",
      direction = -1,
      name = "Cumulative Incidence (CI)\nper 100.000 Penduduk",
      breaks = breaks_interval,
      labels = label_interval
    ) +
    labs(
      title = paste("Cumulative Incidence (CI) Diare di Jawa Tengah", year),
      subtitle = "Per 100.000 Penduduk",
    ) +
    theme_minimal(base_size = 14) +
    theme(
      panel.background = element_rect(fill = "white", color = NA),
      panel.grid.major = element_line(color = "gray90"),
      panel.grid.minor = element_line(color = "gray95"),
      plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
      plot.subtitle = element_text(hjust = 0.5, size = 12),
      legend.position = "right",
      legend.title = element_text(face = "bold", size = 11),
      legend.text = element_text(size = 9),
      legend.key.height = unit(1.2, "cm"),
      legend.key.width = unit(0.5, "cm")
    )
}

# Plot peta
plot_ci_interval(jateng_merged, "CI_2024", 2024)

Peta tersebut menggambarkan tingkat risiko kejadian diare per 100.000 penduduk yang memberikan perspektif lebih akurat dibandingkan sekadar jumlah kasus absolut karena telah dinormalisasi dengan jumlah populasi. Wilayah dengan warna biru tua dan ungu menunjukkan tingkat insidensi tertinggi yang mencapai angka sekitar 3.332 kasus per 100.000 penduduk. Konsentrasi risiko tinggi ini terlihat jelas di wilayah tengah provinsi yang meliputi Kabupaten Boyolali dan sekitarnya serta munculnya titik risiko tinggi baru di wilayah utara bagian tengah yang sebelumnya mungkin tidak terlalu mencolok pada peta jumlah kasus absolut. Sebaliknya wilayah di sisi timur seperti Blora dan Rembang serta wilayah selatan seperti Wonogiri tetap konsisten menunjukkan warna kuning cerah yang menandakan tingkat risiko penularan yang paling rendah di provinsi tersebut.

library(psych)
library(dplyr)
library(sf)

describe(
  jateng_merged %>%
    st_drop_geometry() %>%         
    dplyr::select(Y, CI_2024)
)

##         vars  n     mean      sd  median trimmed     mad     min      max
## Y          1 35 10709.83 7768.75 8998.00 9689.76 6456.72 1112.00 36516.00
## CI_2024    2 35  1108.54  724.38  985.37 1022.72  586.53   74.49  3331.97
##            range skew kurtosis      se
## Y       35404.00 1.42     1.95 1313.16
## CI_2024  3257.48 1.21     1.18  122.44

Pada variabel Cumulative Incidence (CI) tahun 2024 yang memiliki rata-rata 1.108,54 dan median 985,37 per 100.000 penduduk. Variabel ini memiliki nilai skewness 1,21 dan kurtosis 1,18 yang berarti distribusinya juga miring ke kanan meskipun tingkat kemiringannya sedikit lebih rendah dibandingkan variabel jumlah kasus absolut. Variabilitas insidensi penyakit ini cukup tinggi dengan nilai minimum 74,49 dan maksimum 3.331,97 yang menunjukkan adanya ketimpangan risiko penyakit yang signifikan antarwilayah di Jawa Tengah.

4.3 Autokorelasi Spasial

Setelah dilakukan analisis deskriptif menggunakan cumulative incidence, bisa dilakukan analisis lebih dalam, terutama pada efek spasial. Dilakukan analisis untuk melihat besarnya efek dependensi spasial pada penyebaran Diare.Akan digunakan uji Global Moran’s I dan Geary’s C untuk Uji Global. Uji Local Moran dan Uji Ger-Ord akan digunakan untuk Uji Lokal.

Pertama akan dibuatkan weight menggunakan row-based proportion dengan jarak queen, serta pembuatan Peta Ketentanggaan per Kabupaten/Kota.

library(sf)
library(spdep)
library(sp)
library(ggplot2)
library(dplyr)


# PEMBUATAN MATRIKS KETETANGGAAN

# Konversi ke spatial object
nb <- poly2nb(as_Spatial(jateng_merged), queen = TRUE)
lwW <- nb2listw(nb, style = "W", zero.policy = TRUE)
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)

# Koordinat centroid
jateng_sp <- as_Spatial(jateng_merged)
coords <- coordinates(jateng_sp)

# Plot ketetanggaan
plot(jateng_sp, col = "gray90", border = "gray60", 
     main = "Peta Ketetanggaan Kabupaten/Kota Jawa Tengah")
points(coords[,1], coords[,2], pch = 19, cex = 0.7, col = "blue")
text(coords[,1], coords[,2], labels = jateng_sp$KABUPATEN_KOTA, 
     cex = 0.5, pos = 3)
plot(nb, coords, add = TRUE, col = "red", lwd = 1.5)

# FUNCTION ANALISIS AUTOKORELASI SPASIAL

analyze_spatial <- function(data_sf, var_col, label) {
  cat("\n\n==============================\n", label, "\n==============================\n")
  
  x_raw <- data_sf[[var_col]]
  
  # Uji autokorelasi spasial global (asimptotik)
  moran_res <- moran.test(x_raw, listw = lwW, zero.policy = TRUE)
  geary_res <- geary.test(x_raw, listw = lwW, zero.policy = TRUE)
  
  print(moran_res)
  print(geary_res)
  
  # LISA (Local Moran's I)
  x <- scale(x_raw)[, 1]
  lagx <- lag.listw(lwW, x, zero.policy = TRUE)
  lisa <- localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
  lisa_df <- as.data.frame(lisa)
  names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
  
  alpha <- 0.05
  quad <- dplyr::case_when(
    x >= 0 & lagx >= 0 ~ "High-High",
    x < 0 & lagx < 0 ~ "Low-Low",
    x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
    x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
  )
  
  LISA_sf <- bind_cols(data_sf, lisa_df) |>
    mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
  
  # Local Geary's C
  localC_vals <- localC(x_raw, lwW, zero.policy = TRUE)
  Geary_sf <- mutate(data_sf, LocalC = as.numeric(localC_vals))
  
  # Getis–Ord Gi*
  Wb <- listw2mat(lwB)
  Wb_star <- Wb; diag(Wb_star) <- 1
  G_star_raw <- as.numeric(Wb_star %*% x_raw) / sum(x_raw)
  Gz <- localG(x_raw, listw = lwW, zero.policy = TRUE)
  
  G_sf <- mutate(data_sf,
                 z_Gistar = as.numeric(Gz),
                 hotcold = case_when(
                   z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
                   z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
                   TRUE ~ "Not significant"
                 ))
  
  # Visualisasi
  p1 <- ggplot(LISA_sf) +
    geom_sf(aes(fill = quad), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "High-High" = "#d73027",
      "Low-Low" = "#4575b4",
      "High-Low (Outlier)" = "#fdae61",
      "Low-High (Outlier)" = "#74add1",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Local Moran's I (LISA) —", label),
         fill = "Kategori") +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  
  p2 <- ggplot(Geary_sf) +
    geom_sf(aes(fill = LocalC), color = "white", size = 0.2) +
    scale_fill_viridis_c(option = "magma", direction = -1,
                         name = "Local Geary's C") +
    labs(title = paste("Local Geary's C —", label)) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  
  p3 <- ggplot(G_sf) +
    geom_sf(aes(fill = hotcold), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "Hot spot (p<0.05)" = "#b2182b",
      "Cold spot (p<0.05)" = "#2166ac",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Getis–Ord Gi* —", label), fill = NULL) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  
  print(p1); print(p2); print(p3)
  
  # Return objek untuk analisis lanjutan
  invisible(list(
    moran = moran_res,
    geary = geary_res,
    LISA = LISA_sf,
    LocalGeary = Geary_sf,
    Gistar = G_sf
  ))
}

# ANALISIS UNTUK CUMULATIVE INCIDENCE

# Hitung CI per 100.000 penduduk
jateng_merged <- jateng_merged %>%
  mutate(CI_2024 = (Y / Populasi) * 100000)

# Jalankan analisis autokorelasi spasial
hasil_ci <- analyze_spatial(jateng_merged, "CI_2024", 
                            "Cumulative Incidence (CI) Diare 2024")

## 
## 
## ==============================
##  Cumulative Incidence (CI) Diare 2024 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = -0.46628, p-value = 0.6795
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.075550951      -0.029411765       0.009791471 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = -0.74575, p-value = 0.7721
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.08896229        1.00000000        0.01423074

# ANALISIS UNTUK KASUS (Y)

hasil_kasus <- analyze_spatial(jateng_merged, "Y", 
                               "Kasus Diare 2024")

## 
## 
## ==============================
##  Kasus Diare 2024 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = -0.94758, p-value = 0.8283
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.12183649       -0.02941176        0.00951355 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = -1.3997, p-value = 0.9192
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        1.17240354        1.00000000        0.01517151

cat("\n\nRINGKASAN MORAN'S In")

## 
## 
## RINGKASAN MORAN'S In

cat("CI_2024 - Moran's I:", round(hasil_ci$moran$estimate[1], 4), 
    "| p-value:", format.pval(hasil_ci$moran$p.value, digits = 4), "\n")

## CI_2024 - Moran's I: -0.0756 | p-value: 0.6795

cat("Kasus Y - Moran's I:", round(hasil_kasus$moran$estimate[1], 4), 
    "| p-value:", format.pval(hasil_kasus$moran$p.value, digits = 4), "\n")

## Kasus Y - Moran's I: -0.1218 | p-value: 0.8283

Hasil pengujian statistik menunjukkan bahwa tidak terdapat autokorelasi spasial yang signifikan pada kedua variabel yang diuji sehingga pola penyebaran kasus maupun insidensi diare di Jawa Tengah pada tahun 2024 dinyatakan bersifat acak atau random.

Pada variabel Cumulative Incidence (CI) tahun 2024 diperoleh nilai Indeks Moran sebesar -0,0756 dengan p-value 0,6795. Nilai probabilitas yang jauh di atas taraf signifikansi 0,05 ini menyebabkan keputusan gagal tolak hipotesis nol yang berarti tidak ada bukti statistik yang cukup untuk menyatakan adanya hubungan spasial antarwilayah tetangga. Pola yang sama juga ditemukan pada variabel jumlah kasus (Kasus Y) yang menghasilkan Indeks Moran bernilai negatif sebesar -0,1218 dengan p-value 0,8283. Nilai indeks yang negatif dan mendekati nol ini mengindikasikan kecenderungan pola penyebaran yang sangat lemah ke arah menyebar atau dispersed namun karena nilai p-value tidak signifikan maka pola tersebut dianggap terbentuk secara kebetulan.

4.4 Analisis Menggunakan Regresi Sebagai Alternatif

Dikarenakan tidak terdapat autokorelasi spasial yang signifikan maka akan dilakukan analisis menggunakan metode regresi. pertama-tama akan dicek terlebih dahulu overdispersinya untuk menentukan regresi apa yang digunakan

4.4.1 Overdispersi

# Membandingkan Mean dan Varians dari variabel Y (Kasus)
mean_y <- mean(jateng_merged$Y, na.rm = TRUE)
var_y <- var(jateng_merged$Y, na.rm = TRUE)
ratio_var_mean <- var_y / mean_y

cat("\nPengecekan Deskriptif Raw Data:\n")

## 
## Pengecekan Deskriptif Raw Data:

cat("Mean Kasus (Y)     :", round(mean_y, 2), "\n")

## Mean Kasus (Y)     : 10709.83

cat("Varians Kasus (Y)  :", round(var_y, 2), "\n")

## Varians Kasus (Y)  : 60353420

cat("Rasio Varians/Mean :", round(ratio_var_mean, 2), "\n")

## Rasio Varians/Mean : 5635.33

if(ratio_var_mean > 1) {
  cat("-> Indikasi Awal: Terjadi OVERDISPERSI pada data mentah (Var > Mean).\n")
} else {
  cat("-> Indikasi Awal: Tidak terjadi overdispersi (Equidispersi/Underdispersi).\n")
}

## -> Indikasi Awal: Terjadi OVERDISPERSI pada data mentah (Var > Mean).

4.4.2 Poisson GLM Regression

model_poisson <- glm(Y ~ 1, 
                     family = poisson(link = "log"), 
                     offset = log(Populasi), 
                     data = jateng_merged)

# Hitung nilai dispersi: Residual Deviance / Degrees of Freedom
dev_over_df <- model_poisson$deviance / model_poisson$df.residual

cat("\n\nPengecekan Berbasis Model Poisson (Intercept Only):\n")

## 
## 
## Pengecekan Berbasis Model Poisson (Intercept Only):

cat("Residual Deviance  :", round(model_poisson$deviance, 2), "\n")

## Residual Deviance  : 127048.6

cat("Degrees of Freedom :", model_poisson$df.residual, "\n")

## Degrees of Freedom : 34

cat("Nilai Dispersi (Dev/Df) :", round(dev_over_df, 4), "\n")

## Nilai Dispersi (Dev/Df) : 3736.722

if(dev_over_df > 1) {
  cat("-> Kesimpulan: Nilai > 1 mengindikasikan OVERDISPERSI.\n")
} else {
  cat("-> Kesimpulan: Tidak ada overdispersi.\n")
}

## -> Kesimpulan: Nilai > 1 mengindikasikan OVERDISPERSI.

4.4.3 Negative Binomial Regression

model_nb <- glm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, 
                   data = jateng_merged)

summary(model_nb)

## 
## Call:
## glm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = jateng_merged)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -2.877e+04  4.425e+04  -0.650  0.52085   
## X1           5.775e+02  1.681e+02   3.434  0.00187 **
## X2           5.989e+01  5.761e+02   0.104  0.91794   
## X3           1.426e-01  8.011e-01   0.178  0.85998   
## X4           2.552e+03  1.912e+03   1.335  0.19278   
## X5           1.539e+01  4.155e+02   0.037  0.97071   
## X6           1.711e+01  1.362e+02   0.126  0.90091   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 50029903)
## 
##     Null deviance: 2052016273  on 34  degrees of freedom
## Residual deviance: 1400837283  on 28  degrees of freedom
## AIC: 728
## 
## Number of Fisher Scoring iterations: 2

4.5 Desain Eksperimen

4.5.1 Studi Kasus-Kontrol: Faktor Risiko Diare

Penyakit diare tetap menjadi salah satu penyebab utama morbiditas dan mortalitas pada anak-anak di negara berkembang. Secara global, diare menyebabkan sekitar 1,17 juta kematian pada tahun 2021, dengan 390.000 di antaranya terjadi pada anak-anak di bawah usia 5 tahun. Di Indonesia, khususnya Jawa Tengah, diare masih menjadi masalah kesehatan masyarakat yang signifikan.

Studi epidemiologi menunjukkan bahwa sekitar 88% kematian terkait diare dapat dikaitkan dengan air yang tidak aman, sanitasi yang tidak memadai, dan kebersihan yang buruk. Faktor-faktor risiko utama yang berkontribusi terhadap kejadian diare di tingkat populasi meliputi: akses terhadap air bersih, fasilitas sanitasi, praktik kebersihan tangan, kepadatan penduduk, tingkat kemiskinan, dan cakupan imunisasi.

4.5.2 Populasi Target dan Sampling

Populasi Target: Seluruh penduduk di 35 kabupaten/kota Jawa Tengah tahun 2024

Unit Analisis: Kabupaten/Kota (n = 35)

Periode Studi: Tahun 2024 (cross-sectional)

Metode Sampling: Census sampling - seluruh 35 kabupaten/kota di Jawa Tengah diikutsertakan

4.5.3 Variabel Utama

4.5.3.1 Variabel Dependen (Outcome)

Y - Jumlah Kasus Diare

• Definisi: Total kasus diare yang dilaporkan di setiap kabupaten/kota tahun 2024
• Satuan: Count data

CI - Cumulative Incidence

• Formula: \(CI = \frac{Y}{Populasi} \times 100.000\)
• Satuan: Per 100.000 penduduk
• Interpretasi: Risiko populasi terkena diare dalam periode 1 tahun

4.5.3.2 Variabel Independen (Faktor Risiko)

Berdasarkan tinjauan literatur, variabel independen yang digunakan adalah:

1. X1 - Persentase Akses Air Bersih (%)
   • Persentase rumah tangga dengan akses air minum yang aman dan terlindungi
2. X2 - Persentase Sanitasi Layak (%)
   • Persentase rumah tangga dengan fasilitas sanitasi yang memadai (jamban sehat)
3. X3 - Kepadatan Penduduk (per km²)
   • Jumlah penduduk per kilometer persegi
4. X4 - Persentase Balita (%)
   • Persentase penduduk usia 0-5 tahun dari total populasi
5. X5 - Tingkat Kemiskinan (%)
   • Persentase penduduk miskin
6. X6 - Cakupan Imunisasi Rotavirus (%)
   • Persentase balita yang mendapat imunisasi rotavirus

4.5.4 Pengumpulan dan Simulasi Data

4.5.4.1 Load Libraries

library(dplyr)
library(ggplot2)
library(MASS)

## 
## Attaching package: 'MASS'

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

library(broom)
library(knitr)
library(scales)
library(sf)
library(spdep)
library(spatialreg)

set.seed(42)

4.5.4.2 Simulasi Data

# Daftar kabupaten/kota Jawa Tengah
kabupaten <- c(
  "Cilacap", "Banyumas", "Purbalingga", "Banjarnegara", "Kebumen",
  "Purworejo", "Wonosobo", "Magelang", "Boyolali", "Klaten",
  "Sukoharjo", "Wonogiri", "Karanganyar", "Sragen", "Grobogan",
  "Blora", "Rembang", "Pati", "Kudus", "Jepara",
  "Demak", "Semarang", "Temanggung", "Kendal", "Batang",
  "Pekalongan", "Pemalang", "Tegal", "Brebes", "Magelang (Kota)",
  "Surakarta", "Salatiga", "Semarang (Kota)", "Pekalongan (Kota)", "Tegal (Kota)"
)

# Populasi (dalam ribuan, berdasarkan data riil Jawa Tengah)
populasi <- c(
  2007829, 1828573, 1027333, 1047226, 1397555,
  788265, 909664, 1330656, 1090129, 1284386,
  932680, 1051085, 955116, 997485, 1492891,
  901621, 660166, 1359364, 874632, 1221086,
  1240510, 1080648, 808446, 1052826, 828883,
  1007384, 1523622, 1654836, 2043077, 122150,
  526870, 198920, 1694743, 317524, 282781
)

# Simulasi variabel independen berdasarkan distribusi realistis
n <- 35

# X1: Akses Air Bersih (60-95%)
X1 <- round(runif(n, 60, 95), 2)

# X2: Sanitasi Layak (50-90%)
X2 <- round(runif(n, 50, 90), 2)

# X3: Kepadatan Penduduk (500-10000 per km²)
# Kota lebih padat
X3 <- c(
  runif(29, 500, 3000),  # Kabupaten
  runif(6, 5000, 10000)  # Kota
)
X3 <- round(X3, 0)

# X4: Persentase Balita (8-12%)
X4 <- round(runif(n, 8, 12), 2)

# X5: Tingkat Kemiskinan (5-20%)
X5 <- round(runif(n, 5, 20), 2)

# X6: Cakupan Imunisasi Rotavirus (60-98%)
X6 <- round(runif(n, 60, 98), 2)

# Buat data frame
data_diare <- data.frame(
  KabKota = kabupaten,
  Populasi = populasi,
  X1 = X1,
  X2 = X2,
  X3 = X3,
  X4 = X4,
  X5 = X5,
  X6 = X6
)

# Simulasi jumlah kasus diare (Y) berdasarkan model Poisson
# dengan parameter yang realistis

# Koefisien berdasarkan literatur (log scale untuk Poisson)
beta0 <- -8.5        # Intercept
beta1 <- -0.015      # Air bersih (negatif = protektif)
beta2 <- -0.012      # Sanitasi (negatif = protektif)
beta3 <- 0.0001      # Kepadatan penduduk (positif = risiko)
beta4 <- 0.08        # Balita (positif = risiko)
beta5 <- 0.03        # Kemiskinan (positif = risiko)
beta6 <- -0.02       # Imunisasi (negatif = protektif)

# Hitung log(lambda) = log(expected cases)
log_lambda <- beta0 + 
              beta1 * data_diare$X1 + 
              beta2 * data_diare$X2 + 
              beta3 * data_diare$X3 + 
              beta4 * data_diare$X4 + 
              beta5 * data_diare$X5 + 
              beta6 * data_diare$X6 + 
              log(data_diare$Populasi)  # offset

# Expected rate
lambda <- exp(log_lambda)

# Generate kasus dari distribusi Poisson
data_diare$Y <- rpois(n, lambda)

# Hitung Cumulative Incidence per 100,000 penduduk
data_diare$CI_2024 <- (data_diare$Y / data_diare$Populasi) * 100000

# Tampilkan data
kable(head(data_diare, 10), 
      caption = "Tabel 1. Data Simulasi Diare Jawa Tengah 2024 (10 Observasi Pertama)",
      digits = 2)

Tabel 1. Data Simulasi Diare Jawa Tengah 2024 (10 Observasi Pertama)

KabKota	Populasi	X1	X2	X3	X4	X5	X6	Y	CI_2024
Cilacap	2007829	92.02	83.32	607	11.85	7.23	65.19	35	1.74
Banyumas	1828573	92.80	50.29	851	10.96	6.20	91.34	23	1.26
Purbalingga	1027333	70.01	58.31	1041	10.93	11.96	82.51	31	3.02
Banjarnegara	1047226	89.07	86.26	1698	10.14	16.69	90.19	12	1.15
Kebumen	1397555	82.46	74.47	994	8.01	16.00	89.22	19	1.36
Purworejo	788265	78.17	65.18	2298	10.44	17.26	94.89	17	2.16
Wonosobo	909664	85.78	67.43	520	11.35	7.55	92.78	9	0.99
Magelang	1330656	64.71	51.50	1439	11.01	19.17	72.05	57	4.28
Boyolali	1090129	82.99	88.94	1786	9.81	9.40	69.85	17	1.56
Klaten	1284386	84.68	67.27	504	10.14	7.24	88.21	14	1.09

4.5.4.3 Statistik Deskriptif

# Statistik deskriptif
desc_stats <- data_diare %>%
  dplyr::select(Y, CI_2024, X1, X2, X3, X4, X5, X6, Populasi) %>%
  summarise(
    across(everything(),
      list(
        Mean = ~mean(., na.rm = TRUE),
        SD = ~sd(., na.rm = TRUE),
        Min = ~min(., na.rm = TRUE),
        Max = ~max(., na.rm = TRUE),
        Median = ~median(., na.rm = TRUE)
      ),
      # PERBAIKAN 1: Gunakan "__" (double underscore) sebagai pemisah
      .names = "{.col}__{.fn}" 
    )
  ) %>%
  tidyr::pivot_longer(
    cols = everything(),
    names_to = c("Variabel", "Statistik"),
    # PERBAIKAN 2: Pisahkan berdasarkan "__"
    names_sep = "__",           
    values_to = "Nilai"
  ) %>%
  tidyr::pivot_wider(
    names_from = Statistik, 
    values_from = Nilai
  )

# Tampilkan tabel
kable(desc_stats,
      caption = "Tabel 2. Statistik Deskriptif Variabel Penelitian",
      digits = 2)

Tabel 2. Statistik Deskriptif Variabel Penelitian

Variabel	Mean	SD	Min	Max	Median
Y	23.31	13.36	3.00	57.00	20.00
CI_2024	2.41	1.56	0.91	8.13	2.04
X1	81.04	10.32	60.14	94.61	82.99
X2	71.98	11.96	50.29	89.31	74.47
X3	2526.89	2451.13	501.00	9712.00	1786.00
X4	10.42	1.02	8.01	11.85	10.45
X5	11.78	4.02	5.44	19.17	11.41
X6	79.37	11.64	62.02	96.77	79.65
Populasi	1072598.91	470430.21	122150.00	2043077.00	1047226.00

4.5.4.4 Perhitungan Cumulative Incidence

## Ringkasan CI
ci_summary <- data_diare %>%
  dplyr::select(KabKota, Y, Populasi, CI_2024) %>%  # Gunakan dplyr::select
  arrange(desc(CI_2024)) %>%
  mutate(Rank = row_number())

# Top 5 dan Bottom 5
top5 <- head(ci_summary, 5)
bottom5 <- tail(ci_summary, 5)

kable(top5, 
      caption = "Tabel 3. Lima Kabupaten/Kota dengan CI Tertinggi",
      digits = 2)

Tabel 3. Lima Kabupaten/Kota dengan CI Tertinggi

KabKota	Y	Populasi	CI_2024	Rank
Tegal (Kota)	23	282781	8.13	1
Surakarta	34	526870	6.45	2
Sragen	46	997485	4.61	3
Magelang	57	1330656	4.28	4
Jepara	51	1221086	4.18	5

kable(bottom5, 
      caption = "Tabel 4. Lima Kabupaten/Kota dengan CI Terendah",
      digits = 2)

Tabel 4. Lima Kabupaten/Kota dengan CI Terendah

	KabKota	Y	Populasi	CI_2024	Rank
31	Pekalongan	11	1007384	1.09	31
32	Klaten	14	1284386	1.09	32
33	Tegal	17	1654836	1.03	33
34	Wonosobo	9	909664	0.99	34
35	Rembang	6	660166	0.91	35

4.5.4.5 Model Regresi Poisson

# Model Poisson dengan offset populasi
model_poisson <- glm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + offset(log(Populasi)), 
                     data = data_diare, 
                     family = poisson(link = "log"))

# Summary model
summary(model_poisson)

## 
## Call:
## glm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + offset(log(Populasi)), 
##     family = poisson(link = "log"), data = data_diare)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -8.115e+00  7.420e-01 -10.938  < 2e-16 ***
## X1          -1.513e-02  3.451e-03  -4.385 1.16e-05 ***
## X2          -1.127e-02  3.005e-03  -3.752 0.000176 ***
## X3           9.087e-05  1.677e-05   5.417 6.05e-08 ***
## X4           6.253e-02  3.975e-02   1.573 0.115708    
## X5           3.043e-02  8.783e-03   3.465 0.000530 ***
## X6          -2.343e-02  3.489e-03  -6.717 1.86e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.107  on 34  degrees of freedom
## Residual deviance:  37.367  on 28  degrees of freedom
## AIC: 219.65
## 
## Number of Fisher Scoring iterations: 4

# Uji overdispersion
dispersion <- sum(residuals(model_poisson, type = "pearson")^2) / model_poisson$df.residual
cat("\nDispersion parameter:", round(dispersion, 3), "\n")

## 
## Dispersion parameter: 1.429

if(dispersion > 1.5) {
  cat("PERINGATAN: Terdapat overdispersion. Pertimbangkan model Negative Binomial.\n")
}
  
# Jika overdispersion, gunakan Negative Binomial
if(dispersion > 1.5) {
  model_nb <- glm.nb(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + offset(log(Populasi)), 
                     data = data_diare)
  
  cat("\n=== MODEL NEGATIVE BINOMIAL ===\n")
  summary(model_nb)
  
  # Gunakan model NB untuk IRR
  model_final <- model_nb
} else {
  model_final <- model_poisson
}

4.5.4.6 Incidence Rate Ratio (IRR)

# Hitung IRR dan 95% CI
irr_results <- exp(cbind(
  IRR = coef(model_final),
  confint(model_final)
))

## Waiting for profiling to be done...

# Buat tabel IRR
irr_table <- data.frame(
  Variabel = c("Intercept",
               "X1: Akses Air Bersih (per 1% kenaikan)",
               "X2: Sanitasi Layak (per 1% kenaikan)",
               "X3: Kepadatan Penduduk (per 1 unit kenaikan)",
               "X4: Persentase Balita (per 1% kenaikan)",
               "X5: Tingkat Kemiskinan (per 1% kenaikan)",
               "X6: Cakupan Imunisasi (per 1% kenaikan)"),
  IRR = irr_results[, 1],
  CI_Lower = irr_results[, 2],
  CI_Upper = irr_results[, 3]
) %>%
  filter(Variabel != "Intercept")

kable(irr_table, 
      caption = "Tabel 6. Incidence Rate Ratio (IRR) dan 95% Confidence Interval",
      digits = 4)

Tabel 6. Incidence Rate Ratio (IRR) dan 95% Confidence Interval

	Variabel	IRR	CI_Lower	CI_Upper
X1	X1: Akses Air Bersih (per 1% kenaikan)	0.9850	0.9783	0.9917
X2	X2: Sanitasi Layak (per 1% kenaikan)	0.9888	0.9830	0.9946
X3	X3: Kepadatan Penduduk (per 1 unit kenaikan)	1.0001	1.0001	1.0001
X4	X4: Persentase Balita (per 1% kenaikan)	1.0645	0.9852	1.1513
X5	X5: Tingkat Kemiskinan (per 1% kenaikan)	1.0309	1.0134	1.0489
X6	X6: Cakupan Imunisasi (per 1% kenaikan)	0.9768	0.9702	0.9835

::: ## Kesimpulan Eksperimen

4.6 Interpretasi Epidemiologis

5 Kesimpulan dan Saran

Berdasarkan riset yang dilakukan dibutuhkan sebuah riset besar-besaran terhadap deteksi penyakit Diare yang ada di Jawa Tengah dikarenakan meningkatnya kasus serta insidensi Diare pada populasi Jawa Tengah. Karena itu dibutuhkan koordinasi besar antar pemerintah lokal serta institusi kesehatan untuk mencrgah penyebaran Diare dengan konsiderasi faktor risiko pada host untuk memudahkan pencegahan serta antisipasi penyebaran Diare.

6 Daftar Pustaka

Centers for Disease Control and Prevention. (2024). Travelers’ diarrhea.

Curtis, V., & Cairncross, S. (2003). Effect of washing hands with soap on diarrhoea risk in the community: A systematic review. The Lancet Infectious Diseases, 3(5), 275–281.

Fink, G., Günther, I., & Hill, K. (2011). The effect of water and sanitation on child health: Evidence from Demographic and Health Surveys, 1986–2007. International Journal of Epidemiology, 40(5), 1196–1204.

Kruk, M. E., Gage, A. D., Arsenault, C., et al. (2018). High-quality health systems in the Sustainable Development Goals era: Time for a revolution. The Lancet Global Health, 6(11), e1196–e1252.

LeSage, J., & Pace, R. K. (2009). Introduction to Spatial Econometrics. Chapman & Hall/CRC.

Troeger, C., Blacker, B. F., Khalil, I. A., et al. (2020). Global, regional, and national estimates of the burden of diarrhoeal disease in 2017: A systematic analysis. The Lancet Infectious Diseases.

UNICEF. (2021). Diarrhoea.

Dinas Kesehatan Provinsi Jawa Tengah. (2024). Jumlah penemuan penderita diare balita dan semua umur menurut kabupaten/kota tahun 2024. Portal Data Jawa Tengah.

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Provinsi Jawa Tengah dalam angka 2024.

Akhondi, H. (2023). Bacterial diarrhea. In StatPearls. NCBI Bookshelf.

Curtis, V., & Cairncross, S. (2003). Effect of washing hands with soap on diarrhoea risk in the community: A systematic review. The Lancet Infectious Diseases, 3(5), 275–281.

Fink, G., Günther, I., & Hill, K. (2011). The effect of water and sanitation on child health: Evidence from Demographic and Health Surveys, 1986–2007. International Journal of Epidemiology, 40(5), 1196–1204.

Keusch, G. T., et al. (2016). Diarrheal diseases. In G. L. Mandell, J. E. Bennett, & R. Dolin (Eds.), Mandell, Douglas, and Bennett’s Principles and Practice of Infectious Diseases. NCBI Bookshelf.

Kyu, H. H., et al. (2024). Global, regional, and national age–sex–specific burden of disease and injury estimates. The Lancet.

Troeger, C., et al. (2020). Global, regional, and national estimates of the burden of diarrhoeal disease in 2017: A systematic analysis. The Lancet Infectious Diseases.

WHO. (2024). Diarrhoeal disease.

WHO/UNICEF Joint Monitoring Programme. (2021). Progress on drinking water, sanitation and hygiene: 2021 update and SDG baseline.

Budiarto, E. (2012). Metodologi penelitian kedokteran: Sebuah pengantar. EGC.

Celentano, D. D., & Szklo, M. (2018). Gordis epidemiology (6th ed.). Elsevier.

Rothman, K. J., Greenland, S., & Lash, T. L. (2008). Modern epidemiology (3rd ed.). Lippincott Williams & Wilkins.

Hulley, S. B., Cummings, S. R., Browner, W. S., Grady, D. G., & Newman, T. B. (2013). Designing clinical research (4th ed.). Lippincott Williams & Wilkins.

6.1 Lampiran

Link RPubs : https://rpubs.com/Samih/EpidemDiare24

Link Dashboard : https://samihma.shinyapps.io/EpidemDiare24/

Link Video : https://youtu.be/bZB9DOJ08YA