Epidemiologi
From Statistical Methods to Public Health Applications
I Gede Nyoman Mindra Jaya, Ph.D
04 December 2025
Kata Pengantar
1 Pendahuluan Epidemiologi
1.1 Sejarah dan Perkembangan Epidemiologi
Epidemiologi berasal dari kata epi (atas), demos (masyarakat), dan logos (ilmu). Secara etimologis, epidemiologi berarti ilmu yang mempelajari penyakit yang menimpa sekelompok masyarakat. Dalam perkembangannya, epidemiologi tidak hanya mempelajari penyakit menular, tetapi juga berbagai kondisi kesehatan lain yang berdampak pada populasi, termasuk penyakit tidak menular, kecelakaan, dan faktor risiko lingkungan maupun perilaku.
1.1.1 Awal Mula (Era Klasik)
- Pada abad ke-5 SM, Hippokrates sering dianggap
sebagai “Bapak Epidemiologi”. Ia berusaha mencari penyebab penyakit
dengan menghubungkan faktor lingkungan (air, udara, tempat tinggal)
dengan kesehatan manusia. Pandangan ini merupakan pergeseran dari
anggapan bahwa penyakit semata-mata disebabkan oleh kutukan atau hal
mistis.
- Pada abad pertengahan, pencatatan kematian mulai dilakukan, salah satunya oleh John Graunt (1620–1674) di Inggris. Ia menganalisis data kelahiran dan kematian di London, menemukan pola mortalitas, dan mengembangkan konsep life table yang menjadi dasar statistik vital.
1.1.2 Era Modern Awal (Abad ke-19)
- Edward Jenner (1749–1823) mengembangkan vaksin
cacar (smallpox) pada akhir abad ke-18, yang menjadi tonggak awal
pencegahan penyakit melalui imunisasi.
- John Snow (1813–1858) dikenal sebagai pelopor
epidemiologi modern. Pada tahun 1854, ia meneliti wabah kolera di London
dan menunjukkan bahwa penyakit tersebut ditularkan melalui air yang
terkontaminasi. Penelitiannya, terutama peta pompa air Broad
Street, dianggap sebagai contoh klasik penerapan metode
epidemiologi.
- William Farr (1807–1883) mengembangkan sistem pencatatan kesehatan yang sistematis di Inggris dan menghubungkan data mortalitas dengan faktor lingkungan serta demografi.
1.1.3 Era Perkembangan Ilmu Kedokteran dan Statistik (Abad ke-20)
- Penemuan mikroorganisme penyebab penyakit oleh
Louis Pasteur dan Robert Koch memperkuat teori germ theory of
disease.
- Perkembangan metode statistik memungkinkan epidemiologi semakin
kuantitatif. Studi kohort dan case-control mulai digunakan untuk menilai
hubungan sebab-akibat.
- Setelah Perang Dunia II, banyak penelitian besar dilakukan, seperti Framingham Heart Study (1948) yang meneliti faktor risiko penyakit jantung koroner, dan studi kohort perokok di Inggris yang menunjukkan hubungan rokok dengan kanker paru.
1.1.4 Era Kontemporer (Abad ke-21)
- Epidemiologi kini tidak hanya terbatas pada penyakit menular, tetapi
juga mencakup penyakit tidak menular (misalnya diabetes, hipertensi,
kanker), kecelakaan, kesehatan kerja, hingga kesehatan mental.
- Perkembangan teknologi informasi dan biostatistik memungkinkan
analisis data skala besar (big data) dan penggunaan metode
komputasi canggih seperti machine learning.
- Epidemiologi modern terkait erat dengan epidemiologi
spasial, epidemiologi molekuler, serta
epidemiologi digital.
- Pandemi global, seperti COVID-19 (2019–sekarang), menegaskan kembali pentingnya epidemiologi dalam surveilans, analisis data, pemodelan prediktif, serta dasar kebijakan kesehatan.
Ringkasnya: sejarah epidemiologi menunjukkan pergeseran paradigma dari pemahaman mistis menuju pendekatan ilmiah berbasis data.
1.2 Definisi dan Ruang Lingkup Epidemiologi
1.2.1 Definisi Epidemiologi
Secara umum, epidemiologi didefinisikan sebagai ilmu yang mempelajari distribusi dan determinan masalah kesehatan dalam populasi, serta penerapan hasil kajiannya untuk mengendalikan masalah kesehatan.
Beberapa definisi penting:
- Hippokrates (400 SM): penyakit tidak muncul secara
acak, tetapi memiliki penyebab yang dapat ditelusuri dari lingkungan,
pola hidup, dan faktor host.
- MacMahon & Pugh (1970): studi distribusi penyakit
dalam populasi manusia serta faktor yang menentukan distribusinya.
- Last (1988): studi distribusi dan determinan
keadaan/kejadian kesehatan dalam populasi tertentu, serta penerapannya
untuk mengendalikan masalah kesehatan.
- CDC (2002): ilmu dasar kesehatan masyarakat yang
menyediakan metode ilmiah untuk menemukan penyebab, pencegahan, serta
dasar kebijakan kesehatan.
Dengan demikian, epidemiologi tidak hanya menghitung angka penyakit, tetapi juga menganalisis siapa, di mana, kapan, dan mengapa penyakit terjadi.
1.2.2 Ruang Lingkup Epidemiologi
- Epidemiologi Deskriptif
- Mempelajari distribusi penyakit menurut orang,
tempat, dan waktu.
- Contoh: distribusi kasus DBD di Kota Bandung berdasarkan kecamatan.
- Mempelajari distribusi penyakit menurut orang,
tempat, dan waktu.
- Epidemiologi Analitik
- Menganalisis faktor risiko dengan membandingkan kelompok terpapar
dan tidak terpapar.
- Metode: studi case-control dan cohort.
- Menganalisis faktor risiko dengan membandingkan kelompok terpapar
dan tidak terpapar.
- Epidemiologi Eksperimental
- Menguji efektivitas vaksin, obat, atau intervensi melalui randomized controlled trial.
- Epidemiologi Terapan
- Digunakan dalam perencanaan, implementasi, dan evaluasi program
kesehatan masyarakat.
- Termasuk surveilans, deteksi dini outbreak, serta kebijakan berbasis bukti.
- Digunakan dalam perencanaan, implementasi, dan evaluasi program
kesehatan masyarakat.
- Bidang Spesifik Epidemiologi Modern
- Penyakit menular: HIV/AIDS, TBC, malaria, COVID-19.
- Penyakit tidak menular: kanker, hipertensi, diabetes.
- Lingkungan & kerja: polusi, kecelakaan kerja.
- Genetik & molekuler: biomarker, interaksi gen-lingkungan.
- Digital & big data: rekam medis elektronik, media sosial, machine learning.
- Penyakit menular: HIV/AIDS, TBC, malaria, COVID-19.
Kesimpulan: epidemiologi adalah fondasi kesehatan masyarakat, dengan ruang lingkup deskriptif, analitik, eksperimental, terapan, hingga bidang modern.
1.3 Peran Epidemiologi dalam Kesehatan Masyarakat
Epidemiologi disebut sebagai ilmu dasar kesehatan masyarakat (basic science of public health). Perannya meliputi:
- Identifikasi Masalah Kesehatan
- Mengukur besarnya masalah kesehatan dengan morbiditas, mortalitas,
prevalensi, dan insidensi.
- Contoh: surveilans kasus COVID-19 untuk memantau tren wilayah.
- Mengukur besarnya masalah kesehatan dengan morbiditas, mortalitas,
prevalensi, dan insidensi.
- Analisis Faktor Risiko dan Penyebab
- Mengidentifikasi determinan kesehatan biologis, lingkungan, sosial,
dan perilaku.
- Contoh: studi kohort tentang rokok sebagai penyebab kanker paru.
- Mengidentifikasi determinan kesehatan biologis, lingkungan, sosial,
dan perilaku.
- Perencanaan dan Evaluasi Program Kesehatan
- Dasar penyusunan kebijakan kesehatan masyarakat.
- Contoh: evaluasi efektivitas vaksin HPV.
- Dasar penyusunan kebijakan kesehatan masyarakat.
- Surveilans dan Deteksi Dini Wabah
- Pengumpulan, analisis, dan interpretasi data kesehatan secara
sistematis.
- Contoh: sistem peringatan dini DBD berbasis iklim.
- Pengumpulan, analisis, dan interpretasi data kesehatan secara
sistematis.
- Dasar Kebijakan Berbasis Bukti
- Mendukung keputusan kesehatan berbasis hasil kajian ilmiah.
- Contoh: pembatasan iklan rokok berdasarkan bukti epidemiologi.
- Mendukung keputusan kesehatan berbasis hasil kajian ilmiah.
- Pengembangan Ilmu Pengetahuan dan Riset
- Menjadi jembatan antar disiplin (kedokteran, statistik, sosiologi,
lingkungan, data science).
- Contoh: pemodelan spasial-temporal dan machine learning untuk prediksi penyakit.
- Menjadi jembatan antar disiplin (kedokteran, statistik, sosiologi,
lingkungan, data science).
Kesimpulan: epidemiologi sangat penting dalam semua aspek kesehatan masyarakat: identifikasi masalah, analisis risiko, perencanaan, surveilans, hingga kebijakan berbasis bukti.
2 Proses Terjadinya Penyakit
2.1 Konsep Agent–Host–Environment
Terjadinya penyakit pada manusia merupakan hasil interaksi antara agent, host, dan environment. Konsep ini dikenal sebagai epidemiologic triad atau segitiga epidemiologi.
- Agent (Agen Penyebab)
Agen adalah faktor yang keberadaannya diperlukan untuk timbulnya penyakit. Agen bisa berupa:- Biologis: bakteri, virus, jamur, parasit.
- Kimia: racun, polutan, asap rokok.
- Fisik: radiasi, trauma, suhu ekstrem.
- Sosial/psikologis: stres, faktor perilaku.
- Biologis: bakteri, virus, jamur, parasit.
- Host (Pejamu)
Host adalah individu yang menjadi tempat agen hidup dan berkembang. Kerentanan host dipengaruhi oleh:- Faktor genetik.
- Umur, jenis kelamin.
- Status gizi, imunitas.
- Perilaku (gaya hidup, pola makan, kebiasaan).
- Faktor genetik.
- Environment (Lingkungan)
Lingkungan mencakup semua kondisi eksternal yang memengaruhi hubungan host dengan agent. Contoh:- Lingkungan fisik: iklim, cuaca, sanitasi.
- Lingkungan biologis: flora-fauna, vektor penyakit.
- Lingkungan sosial: kepadatan penduduk, kondisi ekonomi, budaya, pelayanan kesehatan.
- Lingkungan fisik: iklim, cuaca, sanitasi.
Hubungan ketiganya menentukan apakah penyakit muncul atau tidak. Penyakit timbul ketika ada ketidakseimbangan antara host, agent, dan lingkungan.
2.2 Model Penyebab Penyakit (Multifaktor)
Tidak semua penyakit disebabkan oleh satu agen saja. Banyak penyakit bersifat multifaktorial, artinya terjadi akibat interaksi berbagai faktor risiko.
2.2.1 Model Web of Causation (jaringan sebab-akibat)
- Menjelaskan penyakit sebagai hasil interaksi kompleks berbagai
faktor yang saling terkait.
- Tidak ada satu faktor tunggal yang cukup, melainkan jaringan
penyebab.
- Contoh: penyakit jantung koroner dipengaruhi oleh hipertensi, merokok, obesitas, kolesterol tinggi, stres, dan faktor genetik.
2.2.2 Model Sufficient-Component Cause (Rothman’s Causal Pie Model)
- Penyakit timbul jika semua komponen penyebab
terbentuk seperti “potongan pai” yang lengkap.
- Setiap faktor (misalnya rokok, genetik, polusi) merupakan bagian
dari “pai penyebab”.
- Berbagai kombinasi faktor dapat menghasilkan penyakit yang sama.
2.2.3 Faktor Determinan Kesehatan
WHO membagi determinan kesehatan ke dalam beberapa kelompok:
- Faktor biologis (umur, genetik).
- Faktor perilaku (gaya hidup).
- Faktor sosial-ekonomi (pendidikan, pendapatan, pekerjaan).
- Faktor lingkungan (air bersih, sanitasi, polusi).
- Akses pelayanan kesehatan.
2.3 Proses Transmisi dan Dinamika Penyakit
2.3.2 Dinamika Penyakit dalam Populasi
- Endemi: penyakit yang selalu ada di suatu wilayah
dengan tingkat relatif konstan (contoh: malaria di Papua).
- Epidemi (KLB/outbreak): peningkatan kasus penyakit
di atas normal pada suatu daerah/waktu tertentu.
- Pandemi: epidemi yang meluas ke banyak negara atau
benua (contoh: COVID-19).
- Eliminasi: penurunan kasus penyakit hingga hampir
tidak ada di suatu wilayah.
- Eradikasi: penghilangan total penyakit dari dunia (contoh: eradikasi cacar oleh WHO tahun 1980).
Kesimpulan: Proses terjadinya penyakit merupakan hasil interaksi antara agen, host, dan lingkungan. Penyakit tidak selalu disebabkan oleh satu faktor, tetapi seringkali merupakan kombinasi multifaktor. Dinamika penyakit dalam populasi mencakup proses transmisi, distribusi, hingga upaya eliminasi atau eradikasi.
3 Jenis Kejadian Penyakit dalam Masyarakat
3.1 Epidemi, Endemi, Pandemi, Outbreak
Dalam epidemiologi, ada beberapa istilah penting untuk menggambarkan pola kejadian penyakit di masyarakat:
Kasus Sporadik
Penyakit yang terjadi secara terpisah-pisah, jumlahnya kecil, dan tidak menunjukkan pola tertentu.
Contoh: kasus rabies pada manusia di wilayah tertentu yang muncul sesekali.Endemi
Penyakit yang selalu ada di suatu wilayah dengan frekuensi relatif konstan dalam jangka panjang.
Contoh: malaria di Papua, atau demam berdarah di daerah tropis.Epidemi (Kejadian Luar Biasa / KLB)
Peningkatan jumlah kasus penyakit dalam suatu wilayah dan waktu tertentu, melebihi angka normal atau baseline.
Contoh: wabah campak di suatu kota dengan lonjakan kasus mendadak.Outbreak
Istilah lain untuk epidemi, biasanya digunakan untuk wilayah yang lebih kecil atau terbatas.
Contoh: outbreak keracunan makanan di sekolah atau pesta pernikahan.Pandemi
Epidemi yang meluas hingga ke berbagai negara atau benua, menjangkiti populasi global.
Contoh: pandemi COVID-19 (2019–sekarang), pandemi influenza 1918.
3.2 Kasus, Insidensi, Prevalensi
Selain istilah epidemi, epidemiologi juga menggunakan ukuran dasar untuk menggambarkan kejadian penyakit:
- Kasus
Individu yang teridentifikasi mengalami penyakit atau kondisi kesehatan tertentu.- Kasus baru (new case): pertama kali terdiagnosis pada
periode tertentu.
- Kasus lama (existing case): sudah terdiagnosis sebelumnya dan masih ada pada periode observasi.
- Kasus baru (new case): pertama kali terdiagnosis pada
periode tertentu.
- Insidensi (Incidence)
Jumlah kasus baru dari suatu penyakit yang muncul dalam populasi tertentu selama periode waktu tertentu.- Rumus sederhana:
\[ \text{Insidensi} = \frac{\text{Jumlah kasus baru dalam periode}}{\text{Jumlah penduduk berisiko}} \times 100\% \]
- Contoh: 100 kasus baru DBD dari 10.000 penduduk berisiko dalam 1 tahun → insidensi = 1%.
- Rumus sederhana:
- Prevalensi (Prevalence)
Jumlah seluruh kasus (baru + lama) suatu penyakit dalam populasi tertentu pada periode waktu tertentu.- Rumus sederhana:
\[ \text{Prevalensi} = \frac{\text{Jumlah seluruh kasus (baru + lama)}}{\text{Jumlah penduduk pada periode tertentu}} \times 100\% \]
- Contoh: jika terdapat 500 orang penderita hipertensi dari 10.000 penduduk → prevalensi = 5%.
- Rumus sederhana:
- Perbedaan utama insidensi dan prevalensi:
- Insidensi → fokus pada kasus baru (new
occurrence).
- Prevalensi → mencakup semua kasus, baik baru maupun
lama.
- Prevalensi dipengaruhi oleh insidensi dan lama penyakit (durasi).
- Insidensi → fokus pada kasus baru (new
occurrence).
Kesimpulan:
Jenis kejadian penyakit dalam masyarakat dapat berupa sporadik, endemi, epidemi/outbreak, maupun pandemi. Untuk mengukur frekuensinya, epidemiologi menggunakan ukuran kasus, insidensi, dan prevalensi, yang menjadi dasar analisis lebih lanjut dalam kesehatan masyarakat.
4 Ukuran Statistik dalam Epidemiologi
Dalam epidemiologi, terdapat beberapa ukuran utama yang digunakan untuk menggambarkan dinamika penyakit dan faktor risikonya. Di antaranya adalah incidence rate, cumulative incidence, prevalensi, serta ukuran asosiasi yang melibatkan exposure seperti Relative Risk (RR), Odds Ratio (OR), dan ukuran tambahan seperti Attributable Risk (AR), Population Attributable Risk (PAR), dan Case Fatality Rate (CFR).
4.1 Incidence Rate (Angka Insidensi)
Definisi: Jumlah kasus baru suatu penyakit dalam periode tertentu dibagi dengan total waktu pemajanan (person-time) dari populasi berisiko. Ukuran ini menggambarkan kecepatan munculnya kasus baru.
Rumus:
\[ IR = \frac{\text{Kasus Baru}}{\text{Total Person-Time}} \]
Contoh: Jika terdapat 5 kasus baru TBC dalam 1 tahun pada populasi 1.000 orang, namun total person-time yang terhitung adalah 950 person-years, maka:
\[ IR = \frac{5}{950} = 0.00526 \; \text{per person-year} \]
atau 5,26 per 1.000 person-years.
4.2 Cumulative Incidence (Insidensi Kumulatif)
Definisi: Proporsi individu yang mengalami penyakit baru dalam suatu periode tertentu di antara populasi berisiko pada awal periode. Menggambarkan risiko rata-rata terkena penyakit dalam periode tersebut.
Rumus:
\[ CI = \frac{\text{Kasus Baru dalam Periode}}{\text{Populasi Berisiko Awal}} \]
Contoh: Dari 1.000 orang sehat di awal tahun, terdapat 20 kasus baru diabetes selama 1 tahun. Maka:
\[ CI = \frac{20}{1000} = 0.02 = 2\% \]
Artinya, risiko terkena diabetes dalam setahun adalah 2%.
4.3 Prevalensi
Definisi: Proporsi individu yang memiliki penyakit tertentu (kasus lama + baru) pada suatu titik waktu (point prevalence) atau periode tertentu (period prevalence). Menggambarkan beban penyakit pada populasi.
Rumus:
\[ P = \frac{\text{Kasus (baru + lama)}}{\text{Total Populasi}} \]
Contoh: Pada 1 Januari 2025, di kota berpenduduk 10.000 orang, terdapat 200 penderita hipertensi. Maka:
\[ P = \frac{200}{10000} = 0.02 = 2\% \]
Prevalensi hipertensi = 2%.
4.4 Exposure dan Ukuran Asosiasi
Dalam epidemiologi, kita sering meneliti hubungan antara exposure (paparan faktor risiko) dan outcome (penyakit). Untuk itu digunakan ukuran asosiasi seperti Relative Risk (RR) dan Odds Ratio (OR).
4.4.1 Relative Risk (RR)
Definisi: Perbandingan risiko penyakit antara kelompok yang terpapar dan kelompok yang tidak terpapar. RR hanya dapat dihitung pada studi kohort.
Rumus:
\[ RR = \frac{CI_{exposed}}{CI_{unexposed}} \]
Contoh:
- Dari 100 orang perokok, 20 terkena kanker paru → \(CI_{exposed} = 20/100 = 0.20\).
- Dari 100 orang bukan perokok, 5 terkena kanker paru → \(CI_{unexposed} = 5/100 = 0.05\).
Maka:
\[ RR = \frac{0.20}{0.05} = 4 \]
Artinya, risiko kanker paru pada perokok 4 kali lebih tinggi dibanding bukan perokok.
4.4.2 Odds Ratio (OR)
Definisi: Perbandingan odds (peluang) terjadinya penyakit pada kelompok terpapar dibanding kelompok tidak terpapar. OR biasanya digunakan pada studi kasus-kontrol.
Rumus (dari tabel 2x2):
| Paparan / Penyakit | Penyakit (+) | Tidak Penyakit (-) |
|---|---|---|
| Terpapar (+) | a | b |
| Tidak Terpapar (-) | c | d |
\[ OR = \frac{a/b}{c/d} = \frac{ad}{bc} \]
Contoh:
- Dari 50 pasien kanker paru: 30 perokok (a=30), 20 bukan perokok (c=20).
- Dari 50 kontrol sehat: 10 perokok (b=10), 40 bukan perokok (d=40).
\[ OR = \frac{30 \times 40}{20 \times 10} = \frac{1200}{200} = 6 \]
Artinya, perokok memiliki 6 kali odds lebih tinggi terkena kanker paru dibanding bukan perokok.
4.5 Attributable Risk (AR) / Risk Difference (RD)
Definisi: Selisih risiko penyakit antara kelompok terpapar dan tidak terpapar. Menggambarkan beban absolut risiko yang diatribusikan pada paparan.
Rumus:
\[ AR = CI_{exposed} - CI_{unexposed} \]
Contoh: Jika \(CI_{exposed} = 0.20\) dan \(CI_{unexposed} = 0.05\), maka:
\[ AR = 0.20 - 0.05 = 0.15 \]
Artinya, ada 15 kasus tambahan per 100 orang pada kelompok terpapar akibat paparan.
4.6 Population Attributable Risk (PAR)
Definisi: Proporsi insidensi penyakit dalam populasi yang dapat dicegah bila paparan dihilangkan. Mengukur dampak paparan pada tingkat populasi.
Rumus:
\[ PAR = CI_{population} - CI_{unexposed} \]
atau dalam bentuk proporsi:
\[ PAR\% = \frac{CI_{population} - CI_{unexposed}}{CI_{population}} \times 100\% \]
Contoh: Jika \(CI_{population} = 0.12\) dan \(CI_{unexposed} = 0.05\), maka:
\[ PAR = 0.12 - 0.05 = 0.07 \]
Artinya, 7 kasus per 100 orang dalam populasi dapat dicegah jika paparan dihilangkan.
4.7 Case Fatality Rate (CFR)
Definisi: Proporsi kasus suatu penyakit yang berakhir dengan kematian. Menggambarkan tingkat keparahan penyakit.
Rumus:
\[ CFR = \frac{\text{Jumlah kematian akibat penyakit}}{\text{Jumlah kasus penyakit}} \times 100\% \]
Contoh: Jika dari 100 kasus COVID-19, terdapat 2 kematian, maka:
\[ CFR = \frac{2}{100} \times 100\% = 2\% \]
4.8 Ringkasan
- Incidence Rate = laju kasus baru per person-time.
- Cumulative Incidence = risiko kasus baru dalam periode tertentu.
- Prevalensi = proporsi total kasus (baru + lama).
- Relative Risk = perbandingan risiko antara terpapar vs tidak terpapar.
- Odds Ratio = perbandingan odds penyakit antara terpapar vs tidak terpapar.
- Attributable Risk = selisih risiko antara terpapar dan tidak terpapar.
- Population Attributable Risk = proporsi insidensi populasi yang dapat dicegah bila paparan dihilangkan.
- Case Fatality Rate = proporsi kasus yang berakhir dengan kematian.
Dengan memahami ukuran-ukuran ini, peneliti dan praktisi kesehatan dapat menggambarkan distribusi penyakit, mengukur beban kesehatan, serta menilai dampak faktor risiko baik pada individu maupun populasi.
4.9 Tugas untuk Mahasiswa
4.9.1 Petunjuk
Jawablah soal-soal berikut dengan menunjukkan langkah perhitungan secara jelas. Gunakan data yang diberikan dan tuliskan interpretasi hasilnya.
4.9.2 Soal 1: Incidence Rate
Dalam sebuah studi, ditemukan 8 kasus baru penyakit X dalam 1 tahun pada populasi 1.200 orang. Karena ada migrasi dan kematian, total person-time yang berhasil dicatat adalah 1.050 person-years.
- Hitung Incidence Rate penyakit X.
- Nyatakan hasil dalam bentuk per 1.000 person-years.
4.9.3 Soal 2: Cumulative Incidence
Dalam sebuah penelitian kohort terhadap 500 orang sehat di awal periode, setelah 2 tahun pengamatan ditemukan 25 kasus baru hipertensi.
- Hitung Cumulative Incidence.
- Bagaimana interpretasinya dalam konteks risiko?
4.9.4 Soal 3: Prevalensi
Pada survei kesehatan di kota berpenduduk 20.000 orang, ditemukan 400 orang menderita diabetes pada saat survei dilakukan.
- Hitung prevalensi diabetes di kota tersebut.
- Apa makna angka tersebut dalam konteks kesehatan masyarakat?
4.9.5 Soal 4: Relative Risk dan Attributable Risk
Data kohort:
- Dari 200 perokok, 40 orang menderita penyakit paru kronis.
- Dari 300 bukan perokok, 15 orang menderita penyakit paru kronis.
Tentukan:
- Cumulative Incidence pada kelompok perokok dan bukan perokok.
- Relative Risk (RR).
- Attributable Risk (AR).
- Interpretasikan hasilnya.
4.9.6 Soal 5: Odds Ratio
Penelitian kasus-kontrol memberikan data berikut:
| Paparan / Penyakit | Penyakit (+) | Tidak Penyakit (-) |
|---|---|---|
| Terpapar (+) | 45 | 30 |
| Tidak Terpapar (-) | 20 | 55 |
Hitung Odds Ratio (OR).
Apa makna hasil tersebut terkait hubungan paparan dengan penyakit?
4.9.7 Soal 6: Case Fatality Rate
Pada suatu wabah, ditemukan 250 kasus dengan 10 di antaranya meninggal.
- Hitung CFR.
- Bagaimana interpretasi tingkat keparahan penyakit berdasarkan CFR tersebut?
4.10 Jawaban Tugas
Petunjuk Umum Tuliskan langkah perhitungan dengan jelas dan berikan interpretasi pada setiap jawaban. Hitung pula jika diminta dalam skala per 1.000 atau 100.000 agar mudah diinterpretasikan.
4.10.1 Soal 1: Incidence Rate
Diketahui: 8 kasus baru selama 1 tahun; total person-time = 1.050 person-years (PY).
Rumus umum:
\[ \text{Incidence Rate (IR)} =
\frac{\text{Kasus baru}}{\text{Total person-time}} \; (\text{per PY})
\]
Langkah hitung:
\[ IR = \frac{8}{1050} = 0{.}007619 \;
\text{per PY} \]
Dalam per 1.000 PY:
\[ IR_{(per\,1000)} = 0{.}007619 \times 1000
= 7{.}619 \approx \mathbf{7{.}62\;kasus/1000\;PY} \]
Interpretasi: Terdapat kira-kira 7,62 kasus baru penyakit X per 1.000 person-years.
cases <- 8; py <- 1050
IR <- cases/py
IR_per1000 <- IR*1000
IR; IR_per1000## [1] 0.007619048## [1] 7.6190484.10.2 Soal 2: Cumulative Incidence
Diketahui: Kohort 500 orang sehat; selama 2 tahun muncul 25 kasus baru.
Rumus:
\[ \text{Cumulative Incidence (CI)} =
\frac{\text{Kasus baru pada periode}}{\text{Populasi berisiko di awal}}
\]
Hitung:
\[ CI = \frac{25}{500} = 0{.}05 =
\mathbf{5\%} \;\text{(selama 2 tahun)} \]
Interpretasi: Risiko seseorang di kohort untuk mengalami hipertensi selama 2 tahun adalah 5% (≈ 1 dari 20 orang).
new_cases <- 25; n0 <- 500
CI <- new_cases/n0
CI## [1] 0.054.10.3 Soal 3: Prevalensi
Diketahui: Populasi 20.000; kasus diabetes saat survei = 400.
Rumus:
\[ \text{Prevalensi} = \frac{\text{Kasus
(lama + baru) saat itu}}{\text{Total populasi}} \]
Hitung:
\[ P = \frac{400}{20000} = 0{.}02 =
\mathbf{2\%} \]
Interpretasi: Pada saat survei, 2% penduduk adalah penyandang diabetes (menggambarkan beban penyakit).
prev <- 400/20000
prev## [1] 0.024.10.4 Soal 4: Relative Risk dan Attributable Risk
Data kohort:
- Perokok: 40 sakit dari 200 → \(CI_{exp} =
40/200\)
- Non-perokok: 15 sakit dari 300 → \(CI_{unexp} = 15/300\)
(1) CI Kelompok:
\[ CI_{exp} = \frac{40}{200} = 0{.}20 = 20\%
\]
\[ CI_{unexp} = \frac{15}{300} = 0{.}05 = 5\%
\]
(2) Relative Risk (RR):
\[ RR = \frac{CI_{exp}}{CI_{unexp}} =
\frac{0{.}20}{0{.}05} = \mathbf{4{.}0} \]
(3) Attributable Risk (AR) / Risk Difference:
\[ AR = CI_{exp} - CI_{unexp} = 0{.}20 -
0{.}05 = \mathbf{0{.}15} = \mathbf{15\%\;points} \]
(Opsional) Attributable Fraction among Exposed
(AF_e):
\[ AF_e = \frac{RR-1}{RR} = \frac{3}{4} =
0{.}75 = 75\% \]
Interpretasi: Perokok memiliki risiko 4 kali lebih tinggi mengalami penyakit paru kronis dibanding bukan perokok. Sekitar 15 poin persentase dari risiko pada perokok dapat diatribusikan pada kebiasaan merokok (≈ 75% dari risiko pada perokok terkait paparan).
a <- 40; n_exp <- 200
c_ <- 15; n_unexp <- 300
CI_exp <- a/n_exp
CI_unexp <- c_/n_unexp
RR <- CI_exp/CI_unexp
AR <- CI_exp - CI_unexp
AF_e <- (RR-1)/RR
c(CI_exp=CI_exp, CI_unexp=CI_unexp, RR=RR, AR=AR, AF_e=AF_e)## CI_exp CI_unexp RR AR AF_e
## 0.20 0.05 4.00 0.15 0.754.10.5 Soal 5: Odds Ratio
Tabel kasus-kontrol:
| Paparan / Penyakit | Penyakit (+) | Tidak Penyakit (-) |
|---|---|---|
| Terpapar (+) | 45 | 30 |
| Tidak Terpapar (-) | 20 | 55 |
Notasi: \(a=45,\; b=30,\; c=20,\; d=55\).
Rumus OR:
\[ OR = \frac{a\times d}{b\times c} =
\frac{45\times 55}{30\times 20} = \frac{2475}{600} = \mathbf{4{.}125}
\]
Interpretasi: Peluang (odds) paparan di antara kasus adalah ≈ 4,13 kali peluang paparan di antara kontrol → ada asosiasi positif antara paparan dan penyakit.
a <- 45; b <- 30; c <- 20; d <- 55
OR <- (a*d)/(b*c)
OR## [1] 4.1254.10.6 Soal 6: Case Fatality Rate (CFR)
Diketahui: 250 kasus; 10 meninggal.
Rumus:
\[ CFR =
\frac{\text{Meninggal}}{\text{Kasus}} \times 100\% \]
Hitung:
\[ CFR = \frac{10}{250} \times 100\% = 4\%
\]
Interpretasi: Sekitar 4% dari kasus pada wabah tersebut berujung kematian → menunjukkan tingkat keparahan penyakit ringan–sedang (interpretasi akhir bergantung pada konteks klinis & pembanding historis).
death <- 10; cases <- 250
CFR <- death/cases*100
CFR## [1] 45 Metode Pengumpulan Data Epidemiologi
Pengumpulan data dalam epidemiologi sangat penting untuk memantau kesehatan masyarakat, mendeteksi wabah, serta mengevaluasi program intervensi. Metode yang digunakan dapat berupa surveilans, sumber data primer maupun sekunder, serta teknik lapangan untuk pencatatan dan registrasi kasus.
5.1 Surveilans Epidemiologi
5.1.1 Surveilans Aktif
- Petugas kesehatan secara proaktif mencari informasi kasus di fasilitas pelayanan, masyarakat, atau populasi tertentu.
- Biasanya dilakukan dengan kunjungan rutin,
wawancara, dan pelaporan langsung oleh petugas.
- Kelebihan: Data lebih cepat, akurat, dan
lengkap.
- Kekurangan: Membutuhkan biaya dan tenaga lebih besar.
5.1.2 Surveilans Pasif
- Mengandalkan laporan rutin dari fasilitas pelayanan
kesehatan (puskesmas, rumah sakit, laboratorium).
- Petugas tidak mencari kasus secara langsung, melainkan menunggu
laporan masuk.
- Kelebihan: Murah, praktis, dapat mencakup wilayah
luas.
- Kekurangan: Data cenderung underreporting (kasus tidak semua tercatat).
5.2 Sumber Data Epidemiologi
5.2.1 Data Primer
- Data yang dikumpulkan langsung dari responden atau populasi
sasaran.
- Contoh: wawancara, kuesioner, observasi, pemeriksaan laboratorium,
pemeriksaan fisik.
- Kelebihan: Data sesuai kebutuhan, lebih
spesifik.
- Kekurangan: Biaya, waktu, dan tenaga lebih besar.
5.2.2 Data Sekunder
- Data yang sudah tersedia dan dikumpulkan oleh pihak lain.
- Contoh: catatan medis, registrasi rumah sakit, laporan surveilans,
sensus penduduk, laporan BPS.
- Kelebihan: Cepat diperoleh, murah, cakupan lebih
luas.
- Kekurangan: Kualitas dan kelengkapan tidak selalu sesuai kebutuhan penelitian.
5.3 Teknik Pengumpulan Data Lapangan dan Registrasi Kasus
- Survei lapangan: Pengumpulan data langsung melalui
kunjungan rumah, posyandu, sekolah, atau fasilitas umum.
- Wawancara: Tanya jawab terstruktur atau
semi-terstruktur dengan responden untuk mendapatkan informasi
kesehatan.
- Observasi: Melihat langsung kondisi lingkungan,
perilaku kesehatan, dan gejala penyakit.
- Registrasi kasus: Pencatatan sistematis semua kasus
penyakit yang ditemukan di fasilitas kesehatan atau masyarakat (contoh:
registrasi kanker, registrasi kelahiran & kematian).
- Formulir standar: Penggunaan kuesioner, formulir WHO, atau sistem elektronik (SIRS, e-health) untuk pencatatan yang konsisten.
5.4 Contoh Kasus Surveilans
5.4.1 Surveilans Dengue Disease
- Deskripsi: Dengue adalah penyakit menular yang
ditularkan melalui nyamuk Aedes aegypti.
- Metode Surveilans:
- Pasif: Laporan kasus demam berdarah dari rumah
sakit dan puskesmas yang dikumpulkan ke dinas kesehatan.
- Aktif: Petugas kesehatan melakukan kunjungan ke
rumah pasien dan pemeriksaan jentik nyamuk di lingkungan sekitar.
- Pasif: Laporan kasus demam berdarah dari rumah
sakit dan puskesmas yang dikumpulkan ke dinas kesehatan.
- Manfaat: Data surveilans digunakan untuk mendeteksi kejadian luar biasa (KLB) dan menentukan kapan dilakukan fogging, PSN (Pemberantasan Sarang Nyamuk), serta kampanye kesehatan.
5.4.2 Surveilans HIV
- Deskripsi: HIV adalah penyakit menular seksual yang
memerlukan surveilans jangka panjang untuk memantau penyebaran dan
efektivitas program pencegahan.
- Metode Surveilans:
- Pasif: Pencatatan jumlah pasien HIV/AIDS di rumah
sakit dan klinik VCT (Voluntary Counseling and Testing).
- Aktif: Survei sentinel pada kelompok berisiko
tinggi (misalnya pekerja seks, pengguna narkoba suntik) untuk mendeteksi
tren penularan.
- Pasif: Pencatatan jumlah pasien HIV/AIDS di rumah
sakit dan klinik VCT (Voluntary Counseling and Testing).
- Manfaat: Surveilans HIV digunakan untuk merancang kebijakan pencegahan, alokasi obat ARV, dan program edukasi masyarakat.
Kesimpulan Metode pengumpulan data epidemiologi mencakup surveilans (aktif dan pasif), pemanfaatan sumber data primer dan sekunder, serta penggunaan teknik lapangan dan registrasi kasus. Contoh pada dengue disease dan HIV menunjukkan bagaimana surveilans membantu deteksi dini, pengendalian penyakit, dan perencanaan program kesehatan masyarakat.
5.5 Bagan Alur Surveilans (Aktif vs Pasif)
Di bawah ini disajikan bagan alur (flowchart) untuk
dua penyakit: Dengue dan HIV.
Saya sediakan dua bentuk: ASCII flowchart (pasti
tampil) dan Mermaid/DiagrammeR (opsional untuk HTML
jika paket tersedia).
5.5.1 Dengue — Surveilans Pasif (pelaporan rutin)
Fasilitas Kesehatan (RS/PKM/Lab)
|
v
Formulir Kasus DBD Terisi (standar Kemenkes)
|
v
Sistem Pelaporan Elektronik (SIRS/EL-Surveilans)
|
v
Dinas Kesehatan Kab/Kota --(validasi & agregasi)--> Dinkes Provinsi
| |
v v
Analisis Tren (insidens, CFR, peta) Analisis Provinsi (KLB?)
| |
v v
Umpan Balik ke Fasyankes + Rekomendasi Tindak Lanjut (PSN, fogging terbatas, komunikasi risiko)5.5.2 Dengue — Surveilans Aktif (pencarian kasus & vektor)
Signal Kenaikan Kasus / KLB
|
v
Tim Lapangan Dinkes Turun (PE: Penyelidikan Epidemiologi)
|
v
- Wawancara Kasus Indeks & Kontak
- Verifikasi Diagnosis (Lab RDT/NS1/IgM/IgG)
- Survey Jentik (ABJ, HI, CI, BI)
|
v
Input Temuan ke Sistem + Pemetaan Spasial
|
v
Tindakan Cepat: PSN 3M+, Fogging Fokus, Edukasi Masyarakat
|
v
Monitoring 7–14 hari & Evaluasi Efektivitas5.5.3 HIV — Surveilans Pasif (registrasi program rutin)
Layanan VCT / CST / RS
|
v
Tes & Konseling (Consent) --> Hasil Lab (ELISA/RT-PCR)
|
v
Pencatatan Kasus (register HIV, ART status, CD4/VL)
|
v
Pelaporan Rutin ke Dinkes + Program HIV
|
v
Analisis Tren (insidens, cascade of care, gap ART)
|
v
Perencanaan Logistik ARV & Intervensi Pencegahan (PMTCT, harm reduction)5.5.4 HIV — Surveilans Aktif (survei sentinel/terfokus)
Penentuan Populasi Kunci (WPS, LSL, Waria, Penasun)
|
v
Desain Survei (Sentinel/IBBS) + Etik + Sampling
|
v
Pengumpulan Data (kuesioner perilaku + tes biologis teranonimkan)
|
v
Analisis Tren Prevalensi & Perilaku Berisiko
|
v
Kebijakan: Targeting Intervensi, Edukasi, Distribusi Kondom/Alat Suntik SterilCatatan: Untuk output HTML, Anda bisa mengaktifkan versi diagram dengan Mermaid/DiagrammeR berikut (opsional). Pastikan paket
DiagrammeRterpasang.
# install.packages("DiagrammeR") # bila belum ada
library(DiagrammeR)
mermaid("
flowchart TD
A[Fasyankes (RS/PKM/Lab)] --> B[Formulir Kasus DBD]
B --> C[Sistem Pelaporan Elektronik]
C --> D[Dinkes Kab/Kota]
D --> E{Validasi & Agregasi}
E --> F[Dinkes Provinsi]
D --> G[Analisis Tren]
F --> H[Analisis Provinsi (KLB?)]
G --> I[Umpan Balik & Rekomendasi]
H --> I
")library(DiagrammeR)
mermaid("
flowchart TD
A[Signal Kenaikan Kasus/KLB] --> B[Tim Lapangan (PE)]
B --> C[Wawancara & Verifikasi Lab]
B --> D[Survey Jentik (ABJ/HI/CI/BI)]
C --> E[Input Sistem + Peta]
D --> E
E --> F[Tindakan Cepat (PSN, Fogging Fokus, Edukasi)]
F --> G[Monitoring 7-14 hari & Evaluasi]
")library(DiagrammeR)
mermaid("
flowchart TD
A[Layanan VCT/CST/RS] --> B[Tes & Konseling -> Hasil Lab]
B --> C[Registrasi HIV (ART, CD4, VL)]
C --> D[Pelaporan ke Dinkes/Program HIV]
D --> E[Analisis Tren: Cascade of Care]
E --> F[Perencanaan ARV & Pencegahan]
")library(DiagrammeR)
mermaid("
flowchart TD
A[Populasi Kunci] --> B[Desain Survei + Etik + Sampling]
B --> C[Pengumpulan Data (Perilaku + Tes)]
C --> D[Analisis Tren Prevalensi & Risiko]
D --> E[Kebijakan & Target Intervensi]
")6 Desain studi epidemiologi
Desain studi epidemiologi digunakan untuk menjawab pertanyaan tentang frekuensi, determinans, dan efektivitas intervensi pada masalah kesehatan. Dua kelompok besar desain adalah studi observasional dan studi eksperimental (intervensi).
6.1 Studi Observasional
Pada studi observasional, peneliti tidak memberikan intervensi; hanya mengamati hubungan antara paparan (exposure) dan kejadian penyakit (outcome).
6.1.1 Cross-Sectional Study
Definisi: Mengukur paparan dan outcome pada
satu titik waktu (snapshot).
Tujuan utama: Mengestimasi prevalensi
dan mengeksplorasi asosiasi awal (hipotesis).
Ciri & Analisis:
- Unit analisis: individu/rumah tangga/populasi.
- Ukuran frekuensi: prevalensi (point/period).
- Ukuran asosiasi: prevalence ratio atau prevalence odds ratio.
- Contoh: Survei hipertensi nasional tahun berjalan.
Kelebihan:
- Relatif cepat dan murah.
- Cocok untuk perencanaan layanan (beban penyakit).
Kelemahan:
- Tidak dapat menentukan urutan waktu (sebab–akibat).
- Tidak tepat untuk penyakit yang jarang atau berumur sangat pendek.
6.1.2 Case-Control Study
Definisi: Memilih kasus (sakit) dan
kontrol (tidak sakit), lalu menelusuri
paparan di masa lalu.
Tujuan utama: Menguji hipotesis asosiasi
paparan–penyakit, efisien untuk penyakit langka.
Ciri & Analisis:
- Arah waktu: retrospektif (dari outcome ke paparan).
- Ukuran asosiasi utama: Odds Ratio (OR).
- Teknik: pemilihan kontrol sebanding, matching, dan regresi logistik untuk kontrol confounding.
Kelebihan:
- Efisien untuk penyakit langka atau laten panjang.
- Membutuhkan waktu dan biaya lebih sedikit daripada kohort.
Kelemahan:
- Tidak dapat mengestimasi insidens langsung.
- Rentan bias seleksi dan bias recall.
- Sulit memastikan temporality (paparan benar mendahului penyakit).
6.1.3 Cohort Study
Definisi: Mengikuti kelompok terpapar dan tidak terpapar ke depan waktu untuk melihat kejadian penyakit. Dapat prospektif (ke depan) atau retrospektif (menggunakan data historis).
Ciri & Analisis:
- Mengukur incidence (cumulative incidence, incidence rate).
- Ukuran asosiasi: Relative Risk (RR), Rate Ratio, atau Hazard Ratio (Cox).
- Dapat mempelajari banyak outcome dari satu paparan.
Kelebihan:
- Kuat untuk menilai temporality dan bukti kausal.
- Dapat menghitung risiko, laju, dan RR langsung.
Kelemahan:
- Mahal dan lama (khususnya prospektif).
- Kurang efisien untuk penyakit langka.
- Risiko loss to follow-up dan misclassification paparan/outcome.
6.2 Studi Eksperimental / Intervensi
Pada studi eksperimental, peneliti memberikan intervensi dan mengacak partisipan ke kelompok perlakuan vs kontrol bila memungkinkan.
6.2.1 Randomized Controlled Trial (RCT)
Definisi: Peserta diacak
(randomized) ke kelompok intervensi dan
kelompok kontrol (plasebo/standar).
Tujuan: Menilai efikasi/efektivitas
intervensi (obat, vaksin, program).
Prinsip kunci:
- Randomisasi: menyeimbangkan confounder terukur/tak terukur.
- Blinding (single/double) untuk meminimalkan bias.
- Kontrol: plasebo atau standar terapi.
- Intention-to-treat: analisis berdasarkan alokasi awal.
- Variasi: cluster RCT (acak pada level puskesmas/komunitas), factorial RCT.
Kelebihan:
- Standar emas (gold standard) untuk inferensi kausal.
- Bias minimal bila desain & pelaksanaan baik.
Kelemahan:
- Biaya tinggi, kebutuhan sampel besar, dan isu etika.
- Generalisasi bisa terbatas (setting terkontrol).
6.3 Perbandingan Ringkas Desain
| Aspek | Cross-sectional | Case-control | Cohort | RCT |
|---|---|---|---|---|
| Arah waktu | Sekali ukur | Mundur (retrospektif) | Maju (prospektif) / mundur (retrospektif) | Maju (intervensi) |
| Ukuran utama | Prevalensi | OR | Insidens, RR/HR | Efek kausal (RR, HR) |
| Cocok untuk | Beban penyakit, hipotesis awal | Penyakit langka/latent | Paparan langka | Uji intervensi |
| Kekuatan | Cepat, murah | Efisien untuk penyakit langka | Temporalitas jelas | Bias minimal, kausal kuat |
| Kelemahan | Tidak tahu temporality | Recall/selection bias; tidak ada insidens | Mahal, lama; loss to follow-up | Biaya, etika, generalisasi |
| Bukti kausal | Lemah | Sedang (tergantung bias) | Baik | Sangat kuat |
6.4 Memilih Desain: Panduan Praktis
- Penyakit langka →
case-control.
- Paparan langka (mis. pajanan industri tertentu) →
cohort.
- Uji efektivitas intervensi → RCT/cluster
RCT.
- Estimasi beban penyakit saat ini →
cross-sectional.
- Pertimbangkan etika, biaya/waktu, ketersediaan data, dan feasibility.
6.5 Lampiran: Tabel 2×2 & Rumus OR/RR
Tabel 2×2 (paparan × penyakit):
| Penyakit (+) | Penyakit (−) | Total | |
|---|---|---|---|
| Terpapar (+) | a | b | a+b |
| Tidak terpapar (−) | c | d | c+d |
| Total | a+c | b+d | N |
Odds Ratio (OR): \(OR =
\frac{a\times d}{b\times c}\).
Risk pada terpapar: \(CI_{exp} = a/(a+b)\).
Risk pada tidak terpapar: \(CI_{unexp} = c/(c+d)\).
Relative Risk (RR): \(RR =
CI_{exp}/CI_{unexp}\).
# Contoh kecil menghitung OR dan RR
a <- 45; b <- 30; c <- 20; d <- 55
OR <- (a*d)/(b*c)
CI_exp <- a/(a+b); CI_unexp <- c/(c+d)
RR <- CI_exp/CI_unexp
c(OR=OR, CI_exp=CI_exp, CI_unexp=CI_unexp, RR=RR)## OR CI_exp CI_unexp RR
## 4.1250000 0.6000000 0.2666667 2.25000006.6 Catatan Bias & Confounding (Singkat)
- Bias seleksi: pemilihan partisipan tidak mewakili
populasi target (umum pada case-control).
- Bias informasi: salah klasifikasi paparan/outcome
(recall bias, misclassification).
- Confounding: faktor ketiga yang terkait paparan & outcome → kontrol dengan randomisasi, matching, stratifikasi, regresi multivariat.
6.7 Contoh Case-Control Study
Latar Belakang & Desain Studi Studi case–control tentang hubungan wadah air tidak tertutup di rumah tangga dengan kejadian Demam Berdarah Dengue (DBD) di Kota X (periode Jan–Jun 2025).
- Kasus (case): DBD terkonfirmasi lab
(NS1/IgM/IgG).
- Kontrol: Tetangga tanpa DBD 2 bulan terakhir,
frequency-matched umur & jenis kelamin.
- Paparan: Ada ≥1 wadah air tidak tertutup (validasi observasi lapangan).
Data Ringkas (Tabel 2×2) Berikut data fiktif konsisten dengan contoh:
| Paparan | Kasus (DBD) | Kontrol |
|---|---|---|
| Wadah tak tertutup (Exposed +) | a = 78 | b = 72 |
| Tidak (Exposed −) | c = 42 | d = 168 |
a <- 78; b <- 72; c <- 42; d <- 168
tab <- matrix(c(a,b,c,d), nrow = 2, byrow = TRUE,
dimnames = list(Exposure = c("Exposed +","Exposed -"),
Status = c("Case +","Control -")))
tab## Status
## Exposure Case + Control -
## Exposed + 78 72
## Exposed - 42 168Odds Ratio (OR) & 95% CI
Perhitungan Manual (Wald CI) Rumus OR: \( OR = \frac{a\times
d}{b\times c} \).
SE Wald untuk \(\ln(OR)\): \( SE = \sqrt{1/a + 1/b + 1/c + 1/d}
\).
CI 95% Wald untuk OR: \( \exp\{\ln(OR) \pm 1.96\,SE\} \).
OR <- (a*d)/(b*c)
SE <- sqrt(1/a + 1/b + 1/c + 1/d)
LCI <- exp(log(OR) - 1.96*SE)
UCI <- exp(log(OR) + 1.96*SE)
c(OR = OR, CI95_L = LCI, CI95_U = UCI)## OR CI95_L CI95_U
## 4.333333 2.719828 6.904031Interpretasi singkat: Odds DBD pada rumah tangga dengan wadah air tak tertutup sekitar 4.33 kali dibanding tanpa paparan. CI 95% Wald: 2.72–6.9.
OR & CI dari exact test (Fisher)
fisher.test() memberikan OR dan CI eksak (berbasis
hipergeometrik).
ft <- fisher.test(tab)
c(OR_Fisher = unname(ft$estimate), CI95_L = ft$conf.int[1], CI95_U = ft$conf.int[2], P_value = ft$p.value)## OR_Fisher CI95_L CI95_U P_value
## 4.314131e+00 2.651348e+00 7.101069e+00 2.792061e-10Catatan: Nilai OR eksak bisa sedikit berbeda dari perhitungan Wald, namun interpretasinya sama.
Regresi Logistik
Bentuk Teragregasi (disarankan untuk data 2×2) Kita
memodelkan peluang menjadi kasus sebagai fungsi
paparan menggunakan glm dengan respon
binomial berbobot kasus vs kontrol per
kelompok paparan.
df_agg <- data.frame(
exposed = c(1, 0),
cases = c(a, c),
controls = c(b, d)
)
fit_agg <- glm(cbind(cases, controls) ~ exposed, family = binomial, data = df_agg)
summary(fit_agg)##
## Call:
## glm(formula = cbind(cases, controls) ~ exposed, family = binomial,
## data = df_agg)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3863 0.1725 -8.036 9.30e-16 ***
## exposed 1.4663 0.2376 6.170 6.81e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40.417 on 1 degrees of freedom
## Residual deviance: 0.000 on 0 degrees of freedom
## AIC: 14.821
##
## Number of Fisher Scoring iterations: 3# OR dan CI 95% (Wald) dari koefisien logit
beta <- coef(fit_agg)["exposed"]
se <- sqrt(vcov(fit_agg)["exposed","exposed"])
OR_glm <- exp(beta)
LCI_glm <- exp(beta - 1.96*se)
UCI_glm <- exp(beta + 1.96*se)
c(OR_glm = OR_glm, CI95_L = LCI_glm, CI95_U = UCI_glm)## OR_glm.exposed CI95_L.exposed CI95_U.exposed
## 4.333333 2.719828 6.904031Konsistensi: OR dari regresi logistik (teragregasi) harus sama dengan OR dari perhitungan 2×2 (perbedaan minimal karena pembulatan).
Bentuk Data Individu Kita dapat memperluas data menjadi baris per individu dan menyesuaikan model logistik yang sama.
dat <- rbind(
data.frame(case = 1, exposed = 1, n = a),
data.frame(case = 0, exposed = 1, n = b),
data.frame(case = 1, exposed = 0, n = c),
data.frame(case = 0, exposed = 0, n = d)
)
# ekspansi baris
dat_full <- dat[rep(1:nrow(dat), dat$n), c("case","exposed")]
fit_ind <- glm(case ~ exposed, family = binomial, data = dat_full)
summary(fit_ind)##
## Call:
## glm(formula = case ~ exposed, family = binomial, data = dat_full)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3863 0.1725 -8.036 9.30e-16 ***
## exposed 1.4663 0.2376 6.170 6.81e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 458.29 on 359 degrees of freedom
## Residual deviance: 417.87 on 358 degrees of freedom
## AIC: 421.87
##
## Number of Fisher Scoring iterations: 4exp(cbind(OR = coef(fit_ind), confint.default(fit_ind)))## OR 2.5 % 97.5 %
## (Intercept) 0.250000 0.1782764 0.3505792
## exposed 4.333333 2.7198516 6.9039714Pada data case–control, koefisien
exposeddariglm(logit) mengestimasi log(OR); sehinggaexp(coef)= OR.
Pelaporan
- OR kasar (2×2): 4.33 (CI 95% Wald 2.72–6.90).
- OR (Fisher exact): 4.31 (CI 95% 2.65–7.10, p =
2.79^{-10}).
- Regresi logistik: OR 4.33 (CI 95% 2.72–6.90).
Kesimpulan: Ada asosiasi positif antara wadah air tak tertutup dan DBD pada studi case–control ini; peluang DBD meningkat sekitar 4.33 kali pada rumah tangga dengan paparan.
Catatan Lanjutan
- Untuk matching individual, gunakan
conditional logistic regression
(
survival::clogit).
- Untuk penyesuaian confounder pada data individu, tambahkan
variabel (mis. kepadatan hunian, 3M, repellent, SES) sebagai kovariat di
glm.
- Hindari menafsirkan OR sebagai RR saat outcome tidak jarang; OR mendekati RR hanya jika outcome jarang (\(<10\%\)).
# Analisis Data Epidemiologi
6.8 Pengolahan dan Pembersihan Data Epidemiologi
Pengolahan data epidemiologi merupakan tahap awal yang penting sebelum analisis dilakukan. Beberapa langkah utama meliputi:
- Pemeriksaan kelengkapan data: identifikasi data yang hilang (missing values).
- Pemeriksaan inkonsistensi: periksa adanya duplikasi, kesalahan input, atau outlier yang tidak masuk akal.
- Transformasi data: standardisasi format tanggal, kode wilayah, atau variabel kategorik.
- Pembuatan variabel baru: misalnya kelompok umur, status penyakit, atau paparan faktor risiko.
Tujuan utama pembersihan data adalah memastikan kualitas data agar analisis yang dilakukan valid dan dapat diinterpretasikan dengan benar.
6.9 Analisis Deskriptif dan Inferensial pada Data Epidemiologi
6.9.1 Analisis Deskriptif
Analisis deskriptif bertujuan untuk menggambarkan data epidemiologi secara umum, meliputi:
- Distribusi frekuensi kasus menurut umur, jenis kelamin, lokasi, dan waktu (person, place, time).
- Statistik ringkasan: mean, median, proporsi, prevalensi, dan insidensi.
- Visualisasi data: histogram, bar chart, peta sebaran kasus.
6.9.2 Analisis Inferensial
Analisis inferensial digunakan untuk menguji hubungan atau perbedaan dalam data epidemiologi, meliputi:
- Regresi logistik untuk mengukur odds ratio suatu faktor risiko terhadap penyakit.
- Regresi Poisson atau negatif binomial untuk menghitung rate ratio pada data insidensi.
- Regresi Cox (survival analysis) untuk mengukur hazard ratio.
- Uji asosiasi antara faktor risiko dengan kejadian penyakit menggunakan tabel kontingensi.
6.10 Uji Hipotesis dalam Epidemiologi
Uji hipotesis digunakan untuk menilai apakah perbedaan atau asosiasi yang ditemukan signifikan secara statistik.
- Hipotesis nol (H0): Tidak ada perbedaan/ asosiasi antara faktor risiko dan penyakit.
- Hipotesis alternatif (H1): Ada perbedaan/ asosiasi.
Beberapa uji yang umum digunakan dalam epidemiologi:
- Chi-square test: untuk data kategorik (misalnya hubungan antara paparan dan status penyakit).
- Uji t dan ANOVA: untuk membandingkan rata-rata antar kelompok.
- Uji non-parametrik (Mann-Whitney, Kruskal-Wallis): bila asumsi parametrik tidak terpenuhi.
- Confidence interval: memberikan estimasi rentang nilai parameter populasi.
Interpretasi hasil uji hipotesis harus mempertimbangkan konteks epidemiologi, potensi bias, dan validitas data.
7 Disease Mapping
7.1 Apa itu Disease Mapping
Disease mapping adalah representasi spasial (peta) dari distribusi suatu penyakit, biasanya dalam bentuk peta tematik yang menunjukkan tingkat prevalensi, insidensi, atau risiko penyakit pada wilayah tertentu. Pendekatan ini memadukan epidemiologi dengan geografi dan statistika spasial.
Tujuan Disease Mapping
Mengidentifikasi distribusi geografis penyakit Untuk melihat pola sebaran penyakit di berbagai wilayah, apakah terkonsentrasi (cluster), menyebar merata, atau mengikuti pola tertentu (misalnya di sepanjang sungai, daerah padat penduduk, dll).
Mendeteksi area risiko tinggi (hotspot) Peta dapat membantu menemukan daerah dengan prevalensi tinggi sehingga bisa difokuskan untuk intervensi.
Membandingkan antarwilayah Disease mapping memungkinkan perbandingan prevalensi penyakit antar kabupaten, provinsi, bahkan antarnegara.
Menyediakan dasar perencanaan kebijakan kesehatan Data spasial penyakit sangat penting untuk alokasi sumber daya medis (vaksin, tenaga kesehatan, rumah sakit lapangan).
Mendukung hipotesis epidemiologi Misalnya untuk menduga hubungan penyakit dengan faktor lingkungan (polusi udara, sanitasi, iklim).
Manfaat Disease Mapping dalam Studi Epidemiologi
Visualisasi kompleks menjadi sederhana: Data angka kasus yang besar dapat dipahami lebih cepat jika ditampilkan dalam bentuk peta.
Meningkatkan kewaspadaan: Aparat kesehatan dapat lebih cepat tanggap terhadap wabah bila tahu letak spasial penyebarannya.
Intervensi yang lebih tepat sasaran: Misalnya, program penyemprotan insektisida difokuskan ke daerah dengan peta risiko malaria tertinggi.
Mendukung komunikasi publik: Peta mudah dipahami masyarakat umum, sehingga lebih efektif untuk kampanye kesehatan.
Contoh Disease Mapping
Peta Malaria di Papua Pemetaan insidensi malaria per kabupaten menunjukkan daerah endemis tinggi (Merauke, Mimika) sehingga program distribusi kelambu difokuskan ke sana.
Peta COVID-19 di Indonesia (2020–2022) Pemerintah membuat peta zonasi (merah, oranye, hijau) berdasarkan jumlah kasus aktif. Peta ini digunakan untuk menentukan kebijakan PPKM di tiap daerah.
Peta Stunting di Jawa Barat Data prevalensi stunting divisualisasikan per desa/kecamatan, untuk menentukan desa prioritas dalam program perbaikan gizi anak.
Peta TBC di Kota Besar Cluster pasien TBC di pemukiman padat penduduk sering muncul jelas pada peta, memunculkan hipotesis hubungan dengan kualitas udara dan sanitasi.
Peta Demam Berdarah Dengue (DBD) Kasus DBD sering dipetakan hingga tingkat kelurahan untuk mendukung fogging dan eliminasi sarang nyamuk.
Dalam materi ini kita akan menggunakan data kasus Diare di Kota Bandung untuk contoh pemetaan.
7.2 Persiapan Library dan Data
library(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)
# Arahkan ke folder kerja Anda
setwd("/Users/mindra/Library/CloudStorage/GoogleDrive-mindra@unpad.ac.id/My Drive/Spatial")
# Data kasus
Diare <- read.csv("Diare.csv", sep=";")
# Peta administratif Indonesia level kecamatan (gadm36 level 3)
Indo_Kec <- readRDS('gadm36_IDN_3_sp.rds')
# Ambil hanya Kota Bandung
Bandung <- Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
plot(Bandung)
# Tambahkan id untuk join
Bandung$id <- c(1:30)
# Ubah ke format sf
Bandung_sf <- st_as_sf(Bandung)
# Gabungkan dengan data kasus
Bandung_merged <- Bandung_sf %>%
left_join(Diare, by = "id")7.3 Pemetaan Kasus (Continuous Mapping)
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare), color=NA) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Peta Kasus Diare Kota Bandung",
fill = "Jumlah Kasus")
Interpretasi: Warna merah menandakan kasus lebih tinggi, kuning lebih rendah.
7.4 Pemetaan Diskrit (Discrete Mapping)
# Interval kategori
breaks <- c(-Inf, 550, 800, 1200, Inf)
labels <- c("Very Low", "Low", "High", "Very High")
# Buat variabel diskrit
Bandung_merged$Diare_Discrete <- cut(Bandung_merged$Diare,
breaks = breaks,
labels = labels,
right = TRUE)
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare_Discrete), color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow",
"Low" = "orange",
"High" = "red",
"Very High" = "red3")) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "right") +
labs(title = "Kategori Kasus Diare Kota Bandung",
fill = "Kategori")
Interpretasi: Kategori “High” dan “Very High” dapat mengindikasikan wilayah dengan risiko lebih besar.
7.5 Autokorelasi Spasial (Moran’s I)
Autokorelasi spasial digunakan untuk melihat apakah data memiliki pola keterkaitan antar wilayah tetangga.
# Matriks bobot spasial
row.names(Bandung) <- as.character(1:30)
W <- poly2nb(Bandung, row.names(Bandung), queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE)
WL <- nb2listw(W)
# Plot jaringan tetangga
CoordK <- coordinates(Bandung)
plot(Bandung, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Bandung), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
# Global Moran's I
Global_Moran <- moran.test(Diare$Diare, WL)
Global_Moran##
## Moran I test under randomisation
##
## data: Diare$Diare
## weights: WL
##
## Moran I statistic standard deviate = 0.67226, p-value = 0.2507
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.04045578 -0.03448276 0.01242630Interpretasi: Nilai Moran’s I menunjukkan apakah ada pola clustering (positif) atau sebaran acak.
7.6 Local Moran’s I (LISA)
Local Moran’s I digunakan untuk mengidentifikasi cluster lokal (misalnya high-high, low-low).
Local_Moran <- localmoran(Diare$Diare, WL)
# Ambil informasi quadran
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Bandung_sf$mean_values <- mean_values_char
# Plot hasil LISA
ggplot() +
geom_sf(data = Bandung_sf, aes(fill = mean_values)) +
scale_fill_manual(values = c("Low-Low" = "blue",
"High-Low" = "green",
"Low-High" = "yellow",
"High-High" = "red")) +
labs(title = "Local Moran's I Kota Bandung",
fill = "Cluster Type") +
theme_minimal()
Interpretasi:
High-High (Merah): cluster wilayah dengan kasus tinggi dikelilingi wilayah tinggi.
Low-Low (Biru): cluster wilayah dengan kasus rendah dikelilingi wilayah rendah.
Low-High / High-Low: wilayah dengan nilai berbeda dari tetangganya (potensi outlier spasial).
7.7 Kesimpulan
- Mapping spasial membantu visualisasi distribusi penyakit.
- Continuous mapping menunjukkan gradien kasus.
- Discrete mapping memberi kategori untuk interpretasi lebih jelas.
- Moran’s I (global dan lokal) mendeteksi adanya pola spasial (clustering).
7.8 Standardized Morbidity/Mortality Ratio (SMR)
Materi menjelaskan standardized morbidity/mortality ratio (SMR) serta direct dan indirect standardization dalam epidemiologi. Tujuannya adalah:
- memahami mengapa standardisasi diperlukan saat struktur umur/jenis kelamin berbeda,
- memahami perbedaan metodologis direct vs indirect,
- melakukan perhitungan manual dan replikasi dengan R,
- menyajikan interpretasi dan batasan.
Catatan: SMR pada konteks ini bisa berarti Standardized Morbidity Ratio (morbiditas) atau Standardized Mortality Ratio (mortalitas). Secara matematis identik; yang berbeda hanyalah jenis kejadian (kasus sakit vs kematian).
Mengapa Perlu Standardisasi?
Angka penyakit/kematian mentah sering tidak dapat dibandingkan langsung antar wilayah/kelompok karena perbedaan struktur umur (atau variabel demografis lainnya). Penyakit degeneratif, misalnya, meningkat tajam pada usia lanjut; wilayah yang lebih “tua” akan tampak lebih “sakit” jika kita tidak menyesuaikan perbedaan ini. Standardisasi menghilangkan bias komposisi ini agar perbandingan menjadi adil.
Definisi dan Notasi Dasar
- \(O\): jumlah kasus diamati (observed cases) di populasi yang diteliti.
- \(E\): jumlah kasus yang diharapkan (expected) jika populasi yang diteliti memiliki risiko seperti populasi standar.
- \(N_i\): ukuran populasi yang diteliti pada kelompok umur ke-\(i\).
- \(r_i^{(s)}\): angka kejadian spesifik umur dari populasi standar pada kelompok \(i\) (mis. per 1.000).
- \(c_i\): jumlah kasus pada kelompok \(i\) di populasi yang diteliti.
- \(n_i\): populasi/py pada kelompok \(i\) di populasi yang diteliti.
- \(r_i = c_i/n_i\): rate spesifik umur populasi yang diteliti.
- \(N_i^{(s)}\): populasi standar pada kelompok \(i\).
- \(w_i = N_i^{(s)} / \sum_j N_j^{(s)}\): bobot standar (proporsi standar pada kelompok \(i\)).
7.8.1 Direct vs Indirect Standarization
SMR (Indirect Standardization) didefinisikan sebagai \[ \mathrm{SMR} \;=\; \frac{O}{E}, \qquad E = \sum_i r_i^{(s)} \, N_i. \]
Directly Standardized Rate (DSR) didefinisikan sebagai \[ \mathrm{DSR} \;=\; \sum_i w_i \, r_i \;=\; \frac{\sum_i r_i \, N_i^{(s)}}{\sum_i N_i^{(s)}}. \]
| Aspek | Direct Standardization | Indirect Standardization (SMR) |
|---|---|---|
| Data yang dibutuhkan | Rate spesifik umur (\(r_i\)) stabil di populasi yang diteliti | Rate spesifik umur tersedia di populasi standar (\(r_i^{(s)}\)); di populasi yang diteliti mungkin kecil/tidak stabil |
| Keluaran utama | Angka kejadian terstandar (DSR) yang bisa dibandingkan langsung antar wilayah | Rasio SMR = \(O/E\) yang membandingkan populasi yang diteliti terhadap standar |
| Kegunaan | Membandingkan banyak wilayah secara setara | Menilai apakah satu populasi “lebih/kurang berisiko” daripada standar |
| Stabilitas pada populasi kecil | Bisa tidak stabil (varian besar) | Lebih stabil untuk populasi kecil |
Contoh 1 — Indirect Standardization (SMR)
Misal kita menilai penyakit paru di Daerah X. Kita punya ukuran populasi Daerah X per umur (\(N_i\)) dan rate standar (\(r_i^{(s)}\)) per 1.000 dari populasi nasional. Kita juga punya kasus diamati (\(O_i\)).
indirect_df <- data.frame(
umur = c("20–39", "40–59", "60+"),
N_i = c(4000, 3000, 2000),
r_std_per1000 = c(2, 6, 12),
O_i = c(10, 20, 30)
)
indirect_df$E_i <- indirect_df$r_std_per1000/1000 * indirect_df$N_i
indirect_df## umur N_i r_std_per1000 O_i E_i
## 1 20–39 4000 2 10 8
## 2 40–59 3000 6 20 18
## 3 60+ 2000 12 30 24Perhitungan manual:
O <- sum(indirect_df$O_i)
E <- sum(indirect_df$E_i)
SMR <- O / E
O; E; SMR## [1] 60## [1] 50## [1] 1.2Interpretasi: Jika Daerah X memiliki risiko seperti populasi standar, kita mengharapkan 50 kasus. Nyatanya ada 60 kasus, sehingga SMR = 1.2. Nilai > 1 berarti lebih tinggi dari harapan.
Interval Kepercayaan 95% untuk SMR
Dengan asumsi \(O \sim \mathrm{Poisson}(E \cdot \mathrm{SMR})\), CI eksak (berbasis \(\chi^2\)) untuk SMR adalah: \[ \left[ \frac{\chi^2_{\alpha/2, \, 2O}}{2E}, \; \frac{\chi^2_{1-\alpha/2, \, 2(O+1)}}{2E} \right]. \]
alpha <- 0.05
lower <- if (O > 0) qchisq(alpha/2, 2*O) / (2*E) else 0
upper <- qchisq(1 - alpha/2, 2*(O + 1)) / (2*E)
c(lower=lower, upper=upper)## lower upper
## 0.9157264 1.5446379Pelaporan: SMR = 1.2 (95% CI 0.92–1.54).
Contoh 2 — Direct Standardization (DSR)
Sekarang kita bandingkan Daerah A dan Daerah B. Kita memiliki kasus dan populasi per umur (sehingga dapat menghitung \(r_i\)), serta populasi standar \(N_i^{(s)}\).
direct_df <- data.frame(
umur = c("<40", "40–59", "60+"),
# Populasi & kasus di A
n_A = c(6000, 3000, 1000),
c_A = c(12, 18, 15),
# Populasi & kasus di B
n_B = c(5000, 3500, 1500),
c_B = c(15, 28, 35),
# Populasi standar
N_std = c(5000, 3000, 2000)
)
# Rate spesifik umur (per 1.000 untuk presentasi; gunakan per 1 untuk perhitungan)
direct_df$r_A <- direct_df$c_A / direct_df$n_A
direct_df$r_B <- direct_df$c_B / direct_df$n_B
direct_df$w <- direct_df$N_std / sum(direct_df$N_std)
direct_df$DSR_A <- sum(direct_df$w * direct_df$r_A)
direct_df$DSR_B <- sum(direct_df$w * direct_df$r_B)
# Varians aproksimasi (Greenwood-like): Var(DSR) ≈ sum w_i^2 * r_i / n_i
var_DSR_A <- sum(direct_df$w^2 * direct_df$r_A / direct_df$n_A)
var_DSR_B <- sum(direct_df$w^2 * direct_df$r_B / direct_df$n_B)
se_A <- sqrt(var_DSR_A)
se_B <- sqrt(var_DSR_B)
ci_A <- c(lower = direct_df$DSR_A - 1.96*se_A, upper = direct_df$DSR_A + 1.96*se_A)
ci_B <- c(lower = direct_df$DSR_B - 1.96*se_B, upper = direct_df$DSR_B + 1.96*se_B)
list(
table = direct_df[, c("umur","n_A","c_A","n_B","c_B","N_std","w","r_A","r_B")],
DSR_A = direct_df$DSR_A,
DSR_B = direct_df$DSR_B,
CI_A = ci_A,
CI_B = ci_B
)## $table
## umur n_A c_A n_B c_B N_std w r_A r_B
## 1 <40 6000 12 5000 15 5000 0.5 0.002 0.00300000
## 2 40–59 3000 18 3500 28 3000 0.3 0.006 0.00800000
## 3 60+ 1000 15 1500 35 2000 0.2 0.015 0.02333333
##
## $DSR_A
## [1] 0.0058 0.0058 0.0058
##
## $DSR_B
## [1] 0.008566667 0.008566667 0.008566667
##
## $CI_A
## lower1 lower2 lower3 upper1 upper2 upper3
## 0.003978852 0.003978852 0.003978852 0.007621148 0.007621148 0.007621148
##
## $CI_B
## lower1 lower2 lower3 upper1 upper2 upper3
## 0.006628409 0.006628409 0.006628409 0.010504924 0.010504924 0.010504924Untuk presentasi per 1.000 penduduk:
DSR_A_1000 <- direct_df$DSR_A * 1000
DSR_B_1000 <- direct_df$DSR_B * 1000
CI_A_1000 <- c(lower = (direct_df$DSR_A - 1.96*se_A)*1000,
upper = (direct_df$DSR_A + 1.96*se_A)*1000)
CI_B_1000 <- c(lower = (direct_df$DSR_B - 1.96*se_B)*1000,
upper = (direct_df$DSR_B + 1.96*se_B)*1000)
c(DSR_A_per1000 = DSR_A_1000,
DSR_B_per1000 = DSR_B_1000)## DSR_A_per10001 DSR_A_per10002 DSR_A_per10003 DSR_B_per10001 DSR_B_per10002
## 5.800000 5.800000 5.800000 8.566667 8.566667
## DSR_B_per10003
## 8.566667CI_A_1000; CI_B_1000## lower1 lower2 lower3 upper1 upper2 upper3
## 3.978852 3.978852 3.978852 7.621148 7.621148 7.621148## lower1 lower2 lower3 upper1 upper2 upper3
## 6.628409 6.628409 6.628409 10.504924 10.504924 10.504924Interpretasi: Setelah menyesuaikan struktur umur standar, rate terstandar Daerah B lebih tinggi dari A. Perbedaan ini mencerminkan risiko yang lebih tinggi, bukan sekadar perbedaan struktur umur.
Ringkasan Perbedaan Direct vs Indirect
- Indirect (SMR) membandingkan kasus aktual dengan kasus harapan berdasarkan rate standar: \(\mathrm{SMR} = O/E\).
- Direct (DSR) membandingkan rate antar wilayah setelah menyamakan komposisi umur melalui bobot standar: \(\mathrm{DSR} = \sum w_i r_i\).
- Pilih indirect jika populasi kecil atau
rate umur di populasi yang diteliti tidak stabil/tidak
tersedia.
- Pilih direct jika tersedia data lengkap dan ingin membandingkan beberapa wilayah secara langsung dalam satu skala rate.
Fungsi Bantuan di R
Berikut fungsi generik untuk menghitung SMR (dengan CI eksak) dan DSR (dengan CI normal-approx).
# SMR dengan CI eksak berbasis chi-square
smr_ci <- function(O, E, alpha = 0.05){
lower <- if (O > 0) qchisq(alpha/2, 2*O) / (2*E) else 0
upper <- qchisq(1 - alpha/2, 2*(O + 1)) / (2*E)
list(SMR = O/E, lower = lower, upper = upper)
}
# DSR dengan Var ≈ sum w_i^2 * r_i / n_i
dsr_ci <- function(cases, pop, std_pop, alpha = 0.05){
r <- cases / pop
w <- std_pop / sum(std_pop)
dsr <- sum(w * r)
var <- sum((w^2) * r / pop)
se <- sqrt(var)
z <- qnorm(1 - alpha/2)
list(DSR = dsr, lower = dsr - z*se, upper = dsr + z*se)
}Contoh pakai fungsi:
# Indirect
O <- sum(indirect_df$O_i); E <- sum(indirect_df$E_i)
smr_ci(O, E)## $SMR
## [1] 1.2
##
## $lower
## [1] 0.9157264
##
## $upper
## [1] 1.544638# Direct
dsr_ci(cases = direct_df$c_A, pop = direct_df$n_A, std_pop = direct_df$N_std)## $DSR
## [1] 0.0058
##
## $lower
## [1] 0.003978885
##
## $upper
## [1] 0.007621115dsr_ci(cases = direct_df$c_B, pop = direct_df$n_B, std_pop = direct_df$N_std)## $DSR
## [1] 0.008566667
##
## $lower
## [1] 0.006628445
##
## $upper
## [1] 0.01050489Catatan Praktis dan Batasan
- Konsistensi definisi: pastikan definisi kasus, periode observasi, dan sumber data konsisten antar kelompok.
- Unit rate: jaga konsistensi (per 1, per 1.000, per 100.000) antara \(r_i\), \(E\), dan pelaporan.
- Ketidakstabilan rate pada populasi kecil: pertimbangkan penggabungan umur atau smoothing (mis. Empirical Bayes) bila perlu.
- Pemilihan populasi standar memengaruhi hasil. Laporkan dengan jelas standar yang digunakan (mis. WHO World Standard Population atau populasi nasional).
- Interval kepercayaan DSR di atas memakai pendekatan normal. Untuk kejadian langka, pertimbangkan metode alternatif (mis. gamma-based).
Standardisasi memungkinkan perbandingan yang adil dengan menyesuaikan perbedaan komposisi umur. SMR (indirect) menjawab pertanyaan “apakah suatu populasi lebih/kurang berisiko dibanding standar?”, sedangkan DSR (direct) memberikan rate terstandar yang dapat dibandingkan lintas wilayah.
Rekomendasi pelaporan minimal: sebutkan metode (direct/indirect), standar yang dipakai, tabel input per umur, nilai ringkas (SMR atau DSR) beserta 95% CI, dan interpretasi substantifnya.
7.8.2 Crude Rate vs Standardized Rate
Crude rate (angka kasar) adalah jumlah kejadian dibagi total populasi, tanpa penyesuaian struktur umur/jenis kelamin: \[ \text{Crude Rate} = \frac{\text{Total Kasus}}{\text{Total Populasi}} \times k \] dengan \(k\) lazimnya 1.000 atau 100.000.
Keterbatasan: crude rate bias apabila struktur umur antar wilayah berbeda, terutama untuk penyakit yang sangat dipengaruhi umur. Karena itu, perbandingan antarwilayah dengan crude rate saja bisa menyesatkan.
Mengapa Perlu Standardisasi (Direct & Indirect)
- Direct standardization (DSR): menghasilkan rate terstandar dengan menyamakan bobot umur ke populasi standar. Cocok bila rate per umur pada populasi yang diteliti tersedia dan cukup stabil. Dapat dipakai untuk membandingkan banyak wilayah dalam satu skala rate yang sebanding.
- Indirect standardization (SMR): menghasilkan rasio (O/E) yang membandingkan jumlah kasus aktual terhadap jumlah kasus harapan, ketika kita meminjam rate per umur dari populasi standar. Cocok untuk populasi kecil/angka langka.
Ringkasnya, crude rate menjawab “berapa kejadian apa adanya”, sedangkan standardized rate/ratio menjawab “apakah risikonya berbeda setelah mengendalikan struktur umur”.
Ilustrasi Numerik Kecil
Misalkan ada dua wilayah dengan struktur umur berbeda kuat. Kita bandingkan crude rate vs DSR dan SMR secara ringkas.
crude_demo <- data.frame(
wilayah = c("Kota A", "Kota B"),
kasus_total = c(80, 55),
pop_total = c(10000, 10000)
)
crude_demo$crude_per1000 <- with(crude_demo, kasus_total/pop_total*1000)
crude_demo## wilayah kasus_total pop_total crude_per1000
## 1 Kota A 80 10000 8.0
## 2 Kota B 55 10000 5.5Crude rate memberi kesan Kota A lebih “buruk”. Namun apakah itu murni karena risikonya lebih tinggi, atau karena struktur umur Kota A lebih tua? Jawaban adilnya diperoleh melalui standardisasi (DSR/SMR).
Contoh Perhitungan Terpadu: Crude vs DSR vs SMR
Kita gunakan data pada bagian sebelumnya dan hitung secara berdampingan untuk Daerah A dan Daerah B (DSR) serta Daerah X (SMR).
# Crude untuk A dan B (tanpa umur)
crude_AB <- data.frame(
wilayah = c("A","B"),
kasus_total = c(sum(direct_df$c_A), sum(direct_df$c_B)),
pop_total = c(sum(direct_df$n_A), sum(direct_df$n_B))
)
crude_AB$crude_per1000 <- with(crude_AB, kasus_total/pop_total*1000)
# DSR (sudah dihitung di atas)
DSR_A_1000 <- direct_df$DSR_A * 1000
DSR_B_1000 <- direct_df$DSR_B * 1000
# SMR untuk Daerah X (indirect)
O_X <- sum(indirect_df$O_i)
E_X <- sum(indirect_df$E_i)
SMR_X <- O_X / E_X
list(
crude_AB = crude_AB,
DSR_A_per1000 = DSR_A_1000,
DSR_B_per1000 = DSR_B_1000,
SMR_X = SMR_X
)## $crude_AB
## wilayah kasus_total pop_total crude_per1000
## 1 A 45 10000 4.5
## 2 B 78 10000 7.8
##
## $DSR_A_per1000
## [1] 5.8 5.8 5.8
##
## $DSR_B_per1000
## [1] 8.566667 8.566667 8.566667
##
## $SMR_X
## [1] 1.2Interpretasi umum:
- Crude rate berguna untuk deskripsi cepat, tetapi bisa berat sebelah jika komposisi umur berbeda.
- DSR memungkinkan perbandingan rate yang adil antarwilayah.
- SMR menjawab apakah sebuah populasi memiliki risiko relatif yang lebih tinggi/rendah daripada standar, ketika rate per umur pada populasi kecil tidak stabil.
Visualisasi Sederhana
# install.packages("ggplot2") # jika belum ada
library(ggplot2)
plot_df <- data.frame(
wilayah = factor(rep(c("A","B"), each = 2), levels = c("A","B")),
jenis = factor(rep(c("Crude","DSR"), times = 2), levels = c("Crude","DSR")),
nilai = c(crude_AB$crude_per1000[match("A", crude_AB$wilayah)],
DSR_A_1000,
crude_AB$crude_per1000[match("B", crude_AB$wilayah)],
DSR_B_1000)
)
ggplot(plot_df, aes(x = wilayah, y = nilai, fill = jenis)) +
geom_col(position = position_dodge(width = 0.7), width = 0.6) +
labs(x = NULL, y = "Rate per 1.000", fill = NULL) +
scale_y_continuous(expand = expansion(mult = c(0, .1))) +
theme_minimal(base_size = 12)
Catatan: Pada data sintetis ini, perbedaan antara Crude dan DSR menunjukkan dampak penyesuaian umur. Jika DSR mengecilkan celah antardaerah dibanding crude, berarti sebagian perbedaan tadi hanyalah efek struktur umur.
7.8.3 SMR untuk Pemetaan Penyakit (Disease Mapping) — Model Poisson
Materi ini menyajikan Standardized Morbidity Ratio (SMR) untuk disease mapping dengan asumsi proses Poisson pada area kecil. Fokusnya pada standar internal: \[ E_i = \Big(\frac{\sum_{j=1}^m O_j}{\sum_{j=1}^m N_j}\Big)\,N_i,\qquad SMR_i = \frac{O_i}{E_i},\quad i=1,\dots,m. \]
Diikuti hubungan konsep ini dengan SMR tidak langsung (indirect) pada umumnya, penurunan model probabilistik, cara menghitung CI 95%, serta contoh numerik dan kode R untuk analisis praktis dan flagging hotspot/coldspot.
Notasi
- Memiliki \(m\) lokasi (kecamatan/kabupaten) di suatu region.
- \(O_i\): jumlah kasus teramati di lokasi \(i\) (periode tertentu).
- \(N_i\): populasi (atau person-time) di lokasi \(i\).
- \(\bar r = \frac{\sum_j O_j}{\sum_j N_j}\): crude rate keseluruhan region.
- Expected (standar internal): \(E_i = \bar r\, N_i\).
- SMR (alias SIR) area-\(i\): \(SMR_i = O_i/E_i\).
Intuisi: jika semua lokasi memiliki risiko yang sama dengan rata-rata regional \(\bar r\), maka kasus harapan di lokasi \(i\) adalah \(E_i\). Deviasi \(O_i\) dari \(E_i\) menunjukkan kemungkinan hotspot (\(>1\)) atau coldspot (\(<1\)).
Hubungan dengan SMR Indirect “Umum” SMR tidak langsung
(umum) menggunakan \[
E_i = \sum_a r_a^{(s)}\,N_{i,a},\qquad SMR_i = \frac{O_i}{E_i},
\] dengan \(r_a^{(s)}\) =
rate standar per strata (misalnya kelompok umur) dari
populasi standar (bisa eksternal, seperti nasional/WHO,
atau internal).
Formulasi pada bagian ini adalah kasus khusus: 1) Tidak
ada stratifikasi (satu strata), sehingga \(r^{(s)}\) tunggal.
2) Standar internal: ambil \(r^{(s)} = \bar r = (\sum O)/(\sum
N)\).
Maka \(E_i = r^{(s)}N_i = \bar r N_i\)
dan \(SMR_i = O_i/E_i\) adalah
SMR tak langsung berbasis standar internal per
lokasi.
7.8.4 Model Probabilistik (Poisson) & CI
Model probabilistik lazim untuk small-area disease mapping: \[ O_i \;\sim\; \text{Poisson}(E_i\,\theta_i), \] dengan \(\theta_i\) = relative risk lokasi \(i\). Estimator sederhana: \[ \widehat{\theta}_i = SMR_i = \frac{O_i}{E_i}. \]
95% CI eksak untuk \(\theta_i\) (dan setara untuk \(SMR_i\)) berbasis \(\chi^2\): \[ \left[ \frac{\chi^2_{\alpha/2,\,2O_i}}{2E_i},\; \frac{\chi^2_{1-\alpha/2,\,2(O_i+1)}}{2E_i} \right]. \]
Catatan: pada lokasi dengan \(O_i\) kecil, CI akan lebar; praktik umum adalah melakukan smoothing (Empirical Bayes, BYM/BYM2).
Contoh Numerik Mini (m = 8 lokasi) Kita buat data sintetis agar langkah-langkah transparan.
set.seed(123)
m <- 8
N <- round(runif(m, 8000, 25000)) # populasi area
# Bangun nilai "benar" risiko relatif (beberapa hotspot)
theta_true <- c(1.8, 1.2, 1.0, 0.9, 1.4, 0.8, 1.0, 1.6)
# Tentukan crude rate global (misal 0.6 per 1.000)
rbar <- 0.6/1000
E <- rbar * N
# Simulasi O_i ~ Poisson(E_i * theta_i)
O <- rpois(m, lambda = E * theta_true)
df <- data.frame(
area = paste0("A", 1:m),
N = N, E = E, O = O
)
df$SMR <- with(df, O/E)
df## area N E O SMR
## 1 A1 12889 7.7334 14 1.8103292
## 2 A2 21401 12.8406 22 1.7133156
## 3 A3 14953 8.9718 10 1.1146035
## 4 A4 23011 13.8066 13 0.9415787
## 5 A5 23988 14.3928 25 1.7369796
## 6 A6 8774 5.2644 1 0.1899552
## 7 A7 16978 10.1868 8 0.7853300
## 8 A8 23171 13.9026 24 1.7262958Hitung CI 95% eksak untuk masing-masing area:
smr_ci <- function(O, E, alpha = 0.05){
lower <- if (O > 0) qchisq(alpha/2, 2*O) / (2*E) else 0
upper <- qchisq(1 - alpha/2, 2*(O + 1)) / (2*E)
c(lower = lower, upper = upper)
}
ci_mat <- t(mapply(smr_ci, O = df$O, E = df$E))
df$SMR_L <- ci_mat[, "lower"]
df$SMR_U <- ci_mat[, "upper"]
df## area N E O SMR SMR_L SMR_U
## 1 A1 12889 7.7334 14 1.8103292 0.989723831 3.037425
## 2 A2 21401 12.8406 22 1.7133156 1.073725751 2.593980
## 3 A3 14953 8.9718 10 1.1146035 0.534495719 2.049796
## 4 A4 23011 13.8066 13 0.9415787 0.501350984 1.610128
## 5 A5 23988 14.3928 25 1.7369796 1.124081614 2.564125
## 6 A6 8774 5.2644 1 0.1899552 0.004809249 1.058362
## 7 A7 16978 10.1868 8 0.7853300 0.339049768 1.547413
## 8 A8 23171 13.9026 24 1.7262958 1.106070293 2.568591Visualisasi sederhana SMR + CI per area (tanpa peta geospasial):
op <- par(mar = c(5, 5, 1, 1))
ord <- order(df$SMR, decreasing = TRUE)
y <- df$SMR[ord]
yl <- df$SMR_L[ord]
yu <- df$SMR_U[ord]
labs <- df$area[ord]
plot(seq_along(y), y, ylim = range(c(yl, yu)), xaxt = "n",
xlab = "", ylab = "SMR (O_i / E_i)",
pch = 19)
arrows(x0 = seq_along(y), y0 = yl, x1 = seq_along(y), y1 = yu, angle = 90, code = 3, length = 0.04)
axis(1, at = seq_along(y), labels = labs, las = 2)
abline(h = 1, lty = 2, col = "gray50")
par(op)Flagging hotspot/coldspot (berdasarkan CI vs 1):
df$flag <- ifelse(df$SMR_L > 1, "Hotspot (signif. > 1)",
ifelse(df$SMR_U < 1, "Coldspot (signif. < 1)",
"Tidak signifikan"))
df[, c("area","O","E","SMR","SMR_L","SMR_U","flag")]## area O E SMR SMR_L SMR_U flag
## 1 A1 14 7.7334 1.8103292 0.989723831 3.037425 Tidak signifikan
## 2 A2 22 12.8406 1.7133156 1.073725751 2.593980 Hotspot (signif. > 1)
## 3 A3 10 8.9718 1.1146035 0.534495719 2.049796 Tidak signifikan
## 4 A4 13 13.8066 0.9415787 0.501350984 1.610128 Tidak signifikan
## 5 A5 25 14.3928 1.7369796 1.124081614 2.564125 Hotspot (signif. > 1)
## 6 A6 1 5.2644 0.1899552 0.004809249 1.058362 Tidak signifikan
## 7 A7 8 10.1868 0.7853300 0.339049768 1.547413 Tidak signifikan
## 8 A8 24 13.9026 1.7262958 1.106070293 2.568591 Hotspot (signif. > 1)7.8.5 Empirical Bayes Ringkas (Poisson–Gamma)
Untuk menstabilkan estimasi pada area kecil, pendekatan EB (naive, global) sering dipakai. Intuisinya: shrink \(SMR_i\) ke 1 dengan bobot tergantung “kekuatan data” (besar \(E_i\)). Skema sederhana berikut hanya ilustrasi (bukan turunan rigor mendalam).
# Skema EB global sederhana (ilustrasi):
# Prior Gamma(a, b) untuk theta; posterior ~ Gamma(a + O_i, b + E_i).
# Estimasi Bayes (mean posterior): (a + O_i) / (b + E_i).
# Pilih (a,b) agar mean prior = 1 dan var moderat (misal a=1, b=1) -> prior mean = 1.
a <- 1; b <- 1
df$theta_EB <- (a + df$O) / (b + df$E)
# Bandingkan raw SMR vs EB
head(df[, c("area","SMR","theta_EB")])## area SMR theta_EB
## 1 A1 1.8103292 1.7175441
## 2 A2 1.7133156 1.6617777
## 3 A3 1.1146035 1.1031108
## 4 A4 0.9415787 0.9455243
## 5 A5 1.7369796 1.6891014
## 6 A6 0.1899552 0.3192644Visual perbandingan SMR vs EB:
op <- par(mfrow = c(1,2), mar = c(4,4,2,1))
plot(df$SMR, df$theta_EB, pch = 19, xlab = "SMR (raw)", ylab = "EB (posterior mean)")
abline(0,1,lty=2,col="gray50")
plot(df$E, abs(df$theta_EB - df$SMR), pch = 19, xlab = "E_i", ylab = "|EB - SMR|")
par(op)Praktik lapangan biasanya memakai pendekatan EB lokal atau model Bayesian spasial (BYM/BYM2) yang menggabungkan spatial smoothing melalui conditional autoregressive (CAR) priors—membiarkan area bertetangga saling “berbagi informasi” agar peta risiko lebih halus.
7.8.6 Kaitan dengan Crude Rate & Direct Standardization
- Formulasi \(E_i = (\sum O/\sum
N)N_i\) menggunakan crude rate global \(\bar r\) sebagai standar
internal. Ini setara dengan SMR indirect tanpa
stratifikasi.
- Bila ingin mengontrol struktur umur antarlokasi
(ketika profil umur berbeda tajam), gunakan indirect
“umum” dengan \(E_i = \sum_a
r_a^{(s)} N_{i,a}\) (standar internal/eksternal)
atau pakai direct standardization
(DSR) untuk menghasilkan rate terstandar yang bisa
dibandingkan antarwilayah dalam satu skala.
- Dalam disease mapping, SMR internal praktis sebagai indikator relatif awal; untuk inferensi kuat, pertimbangkan EB/Bayesian.
Ringkasan Praktis
- Hitung \(\bar r=(\sum O)/(\sum
N)\), lalu \(E_i=\bar r
N_i\).
- Hitung \(SMR_i=O_i/E_i\) dan
CI eksak (Poisson–\(\chi^2\)).
- Buat plot SMR+CI; flag area dengan CI di atas/bawah
1.
- Untuk stabilitas pada area kecil: gunakan EB atau
model spasial (BYM/BYM2).
- Laporkan: metode standar (internal/eksternal), periode data, unit rate, CI, dan interpretasi hotspot/coldspot.
Lampiran: Fungsi Bantuan
smr_table <- function(O, N, k = 1000, alpha = 0.05){
stopifnot(length(O) == length(N))
rbar <- sum(O) / sum(N)
E <- rbar * N
SMR <- O / E
lower <- ifelse(O > 0, qchisq(alpha/2, 2*O)/(2*E), 0)
upper <- qchisq(1-alpha/2, 2*(O+1))/(2*E)
data.frame(O = O, N = N, E = E, SMR = SMR, L = lower, U = upper)
}
# Contoh pakai fungsi:
# smr_table(O = df$O, N = df$N)7.9 Poisson-Gamma Model
7.9.1 Konsep Dasar Poisson-Gamma Model
Dalam pemetaan penyakit (disease mapping), kita sering ingin memodelkan jumlah kasus penyakit pada area ke-\(i\), misalnya \(O_i\), dengan mempertimbangkan variasi spasial dan faktor risiko. Model yang paling dasar adalah model Poisson, yaitu:
\[ O_i \mid \lambda_i \sim \text{Poisson}(E_i \lambda_i) \]
di mana:
- \(O_i\): jumlah kasus teramati di area \(i\),
- \(E_i\): jumlah kasus yang diharapkan (expected cases), biasanya berdasarkan populasi atau tingkat dasar (baseline rate),
- \(\lambda_i\): risiko relatif (relative risk) di area ke-\(i\).
Namun, data penyakit sering menunjukkan overdispersion — varians lebih besar dari mean — yang tidak dapat dijelaskan oleh model Poisson murni. Untuk mengakomodasi variasi ini, digunakan model Poisson-Gamma, di mana parameter risiko \(\lambda_i\) dianggap acak dan mengikuti distribusi Gamma:
\[ \lambda_i \sim \text{Gamma}(a, b) \]
Sehingga model hierarkisnya adalah:
\[ O_i \mid \lambda_i \sim \text{Poisson}(E_i \lambda_i), \quad \lambda_i \sim \text{Gamma}(a, b) \]
7.9.2 Distribusi Marginal: Hubungan dengan Negative Binomial
Kita dapat mengintegrasikan \(\lambda_i\) dari model hierarkis di atas untuk mendapatkan distribusi marginal dari \(O_i\).
Distribusi gabungan:
\[ P(O_i = o_i, \lambda_i) = \frac{(E_i \lambda_i)^{o_i} e^{-E_i \lambda_i}}{o_i!} \cdot \frac{b^a}{\Gamma(a)} \lambda_i^{a-1} e^{-b\lambda_i} \]
Marginalisasi terhadap \(\lambda_i\) menghasilkan:
\[ P(O_i = o_i) = \int_0^\infty P(O_i = o_i, \lambda_i) \, d\lambda_i \]
Setelah diselesaikan secara analitik, hasilnya adalah distribusi Negative Binomial (NB):
\[ O_i \sim \text{NegBin}(r = a, p = \frac{b}{b + E_i}) \]
Artinya, model Poisson-Gamma ekuivalen dengan model Negative Binomial dengan:
\[ \mathbb{E}[O_i] = E_i \cdot \frac{a}{b}, \quad \text{Var}(O_i) = E_i \cdot \frac{a}{b} \left(1 + \frac{E_i}{b}\right) \]
Inilah yang membuat model Poisson-Gamma sangat berguna untuk menangani overdispersion: variasi tambahan berasal dari variasi acak pada \(\lambda_i\).
7.9.3 Distribusi Posterior
Kita ingin memperkirakan risiko relatif \(\lambda_i\) setelah mengamati data \(O_i\). Dengan teorema Bayes:
\[ P(\lambda_i \mid O_i) \propto P(O_i \mid \lambda_i) P(\lambda_i) \]
Substitusi:
\[ P(O_i \mid \lambda_i) \propto (E_i \lambda_i)^{O_i} e^{-E_i \lambda_i} \] \[ P(\lambda_i) \propto \lambda_i^{a-1} e^{-b\lambda_i} \]
Maka posteriornya adalah:
\[ \lambda_i \mid O_i \sim \text{Gamma}(a + O_i, b + E_i) \]
Nilai harapan posterior (posterior mean):
\[ \mathbb{E}[\lambda_i \mid O_i] = \frac{a + O_i}{b + E_i} \]
Ini dapat diinterpretasikan sebagai penyusutan Bayes (Bayesian shrinkage): risiko relatif \(\lambda_i\) “ditarik” ke arah rata-rata global \(a/b\), tergantung pada kekuatan informasi data \(O_i\).
7.9.4 Contoh Aplikasi di Disease Mapping
Misalkan kita memiliki data jumlah kasus penyakit di 10 wilayah berikut:
set.seed(123)
area <- 1:10
E <- rep(20, 10) # expected cases
a <- 2; b <- 1 # prior parameters (mean = 2)
# True lambda (unknown risks)
lambda_true <- rgamma(10, a, b)
O <- rpois(10, E * lambda_true)
# Posterior estimation
lambda_post_mean <- (a + O) / (b + E)
data.frame(area, O, E, lambda_true, lambda_post_mean)## area O E lambda_true lambda_post_mean
## 1 1 19 20 0.8920936 1.0000000
## 2 2 67 20 3.3118474 3.2857143
## 3 3 2 20 0.1443750 0.1904762
## 4 4 27 20 1.6625233 1.3809524
## 5 5 91 20 4.3358790 4.4285714
## 6 6 29 20 2.1176157 1.4761905
## 7 7 5 20 0.3507177 0.3333333
## 8 8 2 20 0.1304001 0.1904762
## 9 9 61 20 3.3737820 3.0000000
## 10 10 38 20 1.9730466 1.9047619Interpretasi: nilai lambda_post_mean memberikan estimasi
risiko relatif posterior untuk tiap area. Area dengan sedikit kasus akan
memiliki estimasi yang “menyusut” ke rata-rata global, sementara area
dengan banyak kasus memiliki nilai lebih tinggi.
7.9.5 Hubungan dengan Model Spasial (BYM dan Hierarchical Poisson-Gamma)
Model Poisson-Gamma ini dapat diperluas secara spasial dengan menambahkan struktur dependensi antar area, misalnya dengan prior Conditional Autoregressive (CAR) atau Besag-York-Mollié (BYM). Dalam hal ini, log-risiko dapat dimodelkan sebagai:
\[ \log(\lambda_i) = \alpha + u_i + v_i \]
di mana:
- \(u_i\): efek spasial terstruktur (misal: CAR prior),
- \(v_i\): efek acak tidak terstruktur (iid normal).
Model ini dapat diestimasi dengan INLA (Integrated Nested Laplace Approximation) untuk efisiensi komputasi.
7.9.6 Kesimpulan
Model Poisson-Gamma menyediakan jembatan alami antara Poisson dan Negative Binomial, menjelaskan overdispersion melalui ketidakpastian pada parameter risiko. Posterior dari risiko relatif memiliki bentuk Gamma yang mudah diturunkan secara analitik, menjadikannya model sederhana namun kuat untuk disease mapping.
Untuk kasus spasial yang kompleks, model ini menjadi dasar bagi model hierarkis spasial seperti BYM dan BYM2 yang lebih realistis.
Referensi:
- Clayton, D., & Kaldor, J. (1987). Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics.
- Lawson, A. (2018). Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology. CRC Press.
- Besag, J., York, J., & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics.
8 Bayesian Disease Mapping Model
8.1 Tujuan Pembelajaran
- Memahami konsep Bayesian dalam disease mapping.
- Memahami smoothing (penyugihan/penyeragaman) risiko untuk area kecil.
- Menjelaskan mengapa pendekatan Bayesian diperlukan pada peta penyakit.
- Membedakan Empirical Bayes vs Full Bayes.
- Mengimplementasikan contoh model Poisson–BYM2 dengan R-INLA untuk memetakan relative risk (RR), termasuk peta SMR kasar vs RR tersmoothing, dan peta exceedance probability \(P(\mathrm{RR}>1)\).
Catatan: contoh kode bersifat self-contained (simulasi data + adjacency), sehingga bisa dijalankan tanpa data eksternal.
8.2 Konsep Bayesian dalam Disease Mapping
Pada disease mapping, kita sering memodelkan jumlah kasus \(Y_i\) pada area \(i\) (mis. kecamatan/kabupaten) sebagai \[ Y_i \sim \text{Poisson}(E_i \, \theta_i), \] di mana \(E_i\) adalah ekspektasi kasus jika tidak ada variasi spasial (mis. indirect standardization berdasarkan struktur umur/jenis kelamin) dan \(\theta_i\) adalah relative risk (RR) yang ingin dipetakan.
Pendekatan Bayesian memadukan likelihood dengan prior untuk menghasilkan posterior: \[ p(\boldsymbol{\theta} \mid \mathbf{Y}) \propto p(\mathbf{Y} \mid \boldsymbol{\theta}) \, p(\boldsymbol{\theta}). \] Pada konteks spasial, kita memilih prior yang mendorong kedekatan nilai di area bertetangga (mis. CAR/BYM), sehingga informasi “dipinjam” dari tetangga (borrowing strength) untuk menstabilkan estimasi pada area dengan populasi kecil.
8.3 Apa itu Smoothing? Mengapa Diperlukan?
Smoothing adalah proses menstabilkan estimasi risiko/angka kejadian yang bervariasi liar ketika jumlah populasi atau kasus kecil. Raw rate atau SMR \(\mathrm{SMR}_i = Y_i/E_i\) sering tidak stabil—area dengan \(E_i\) kecil akan memiliki nilai ekstrem karena fluktuasi acak. Bayesian smoothing (atau shrinkage) menarik estimasi menuju rata-rata global maupun rata-rata tetangga (spasial). Hasilnya: - Peta risiko lebih halus, - Pola lebih mudah ditangkap, - Interval kredibel lebih realistis, khususnya untuk area kecil.
8.4 Mengapa Bayesian Penting untuk Disease Mapping?
- Small-area estimation: banyak area berpenduduk kecil → raw rate tidak andal.
- Borrowing strength: informasi dari area tetangga memperbaiki estimasi tiap area.
- Ketidakpastian eksplisit: Bayesian memberi interval kredibel dan kuantitas keputusan seperti probabilitas exceedance \(P(\mathrm{RR}>1)\).
- Model fleksibel: gabungkan structured spatial effect (tetangga) + unstructured heterogeneity (IID), kovariat, efek waktu, dsb.
8.5 Empirical Bayes vs Full Bayes
Empirical Bayes (EB):
- Hyperparameter prior (mis. parameter precision atau parameter Gamma) diestimasi dari data (sering plug-in via MLE atau method of moments), lalu digunakan seolah-olah diketahui.
- Komputasi ringan, cepat, dan sering cukup baik untuk pemetaan ringkas.
- Kekurangan: ketidakpastian hyperparameter tidak sepenuhnya dipropagasi.
Full Bayes:
- Tetapkan prior juga untuk hyperparameter, dan inferensi dilakukan atas seluruh hierarki.
- Ketidakpastian sepenuhnya dipropagasi, sehingga kredibel interval/probabilitas keputusan lebih defensibel.
- Komputasi lebih berat; dengan INLA kita bisa melakukan approximate full Bayes secara efisien untuk model Gaussian tersembunyi (LGM).
8.6 Model Spasial Klasik: Poisson + BYM/BYM2
Struktur umum (untuk area \(i=1,\dots,n\)): \[ \begin{aligned} Y_i &\sim \text{Poisson}(E_i \,\theta_i),\\ \log \theta_i &= \alpha + \beta^\top x_i + u_i + v_i, \end{aligned} \]
dengan:
- \(u_i\): efek spasial terstruktur (mis. prior CAR/Besag: tetangga mirip),
- \(v_i \sim \text{iid Normal}(0, \sigma_v^2)\): heterogenitas tak terstruktur,
- \(\text{BYM}\) menggabungkan \(u_i\) dan \(v_i\); BYM2 adalah reparametrisasi dengan skala yang baik dan parameter mixing \(\phi \in [0,1]\) untuk proporsi varian terstruktur vs tak terstruktur.
- Hyperprior sering memakai PC priors (penalized complexity) agar weakly informative namun stabil.
8.7 Bayesian Disease Mapping with INLA
Pemetaan penyakit (disease mapping) adalah salah satu aplikasi paling populer dari pemodelan spasial dalam epidemiologi. Tujuan utamanya adalah:
- Mengestimasi risiko relatif penyakit pada unit wilayah (misalnya kecamatan, kabupaten, provinsi).
- Menghaluskan (smoothing) angka insiden yang noisy karena ukuran populasi kecil.
- Mengidentifikasi pola klaster spasial yang berpotensi penting untuk kebijakan kesehatan publik.
Dalam pendekatan klasik, disease mapping sering dilakukan dengan membandingkan angka insiden kasar (crude rates) antar wilayah. Namun, angka tersebut bisa sangat tidak stabil jika jumlah penduduk kecil atau kasus jarang. Pendekatan Bayesian mengatasi masalah ini dengan:
- Menggunakan model hierarkis: memodelkan data kasus, risiko yang tidak teramati, dan efek spasial/non-spasial secara eksplisit.
- Menggunakan prior untuk menggabungkan informasi tambahan dan melakukan smoothing.
- Menghasilkan distribusi posterior dari parameter dan risiko relatif, sehingga ketidakpastian dapat dikuantifikasi.
Di antara metode komputasi Bayesian, pendekatan klasik adalah MCMC (Markov Chain Monte Carlo). Namun, untuk model spasial yang kompleks dan dataset besar, MCMC bisa sangat lambat. Di sinilah Integrated Nested Laplace Approximation (INLA) menjadi alternatif yang sangat efisien untuk kelas Latent Gaussian Models (LGM), termasuk sebagian besar model disease mapping.
Dokumen ini membahas:
- Formulasi matematis model disease mapping Bayesian (Poisson log-linear).
- Berbagai spesifikasi efek spasial: IID, ICAR, BYM, BYM2, Leroux.
- Konsep dasar INLA dan bagaimana ia mengaproksimasi posterior.
- Penentuan prior, termasuk Penalized Complexity (PC) priors.
- Implementasi lengkap di R menggunakan paket INLA untuk contoh data simulasi dan contoh realistis.
Tujuan akhirnya adalah menyediakan materi yang cukup lengkap (konseptual, matematis, dan praktis) untuk pembelajaran maupun pengajaran.
8.7.1 Model Dasar Disease Mapping
8.7.2 Notasi dasar
Misalkan kita memiliki \(n\) wilayah (area) yang saling bersebelahan. Untuk setiap wilayah \(i = 1, \dots, n\):
- \(y_i\) : banyaknya kasus penyakit yang diamati.
- \(E_i\) : expected count (jumlah kasus yang diharapkan) jika risiko di area \(i\) sama dengan risiko rata-rata populasi referensi (misalnya berdasarkan struktur umur dan jenis kelamin).
- \(\lambda_i\) : risiko relatif (relative risk) di area \(i\), yaitu apakah area tersebut memiliki risiko lebih tinggi/rendah dari rata-rata.
Model disease mapping standar (Poisson log-linear) adalah:
\[ y_i \mid \lambda_i \sim \text{Poisson}(E_i \lambda_i), \]
dengan
\[ \log(\lambda_i) = \eta_i = \alpha + \mathbf{x}_i^\top \boldsymbol{\beta} + u_i + v_i. \]
Di sini:
- \(\alpha\) : intercept (log-risiko rata-rata).
- \(\mathbf{x}_i\) : vektor kovariat untuk area \(i\) (misalnya indeks kemiskinan, kepadatan penduduk, dll).
- \(\boldsymbol{\beta}\) : koefisien regresi kovariat.
- \(u_i\) : efek spasial terstruktur (area berdekatan cenderung memiliki nilai serupa).
- \(v_i\) : efek tidak terstruktur (IID noise antar area).
Model ini dikenal juga sebagai model BYM (Besag–York–Mollié) ketika \(u_i\) dimodelkan sebagai ICAR (Intrinsic Conditional AutoRegressive) dan \(v_i\) sebagai IID Gaussian.
Dalam bentuk log-mean:
\[ \mu_i = \log(E_i) + \eta_i, \]
sehingga offset \(\log(E_i)\) bisa dimasukkan langsung dalam formulasi model.
8.7.3 Tahap hierarkis
Model Bayesian biasanya ditulis dalam tiga tahap:
Data model (likelihood)
\[ y_i \mid \eta_i \sim \text{Poisson}(\exp(\log(E_i) + \eta_i)). \]Latent field (model untuk \(\eta_i\))
\[ \eta_i = \alpha + \mathbf{x}_i^\top\boldsymbol{\beta} + u_i + v_i. \]Hyperparameters (model untuk parameter varians, precision, dan parameter lain)
Misalnya: \[ \alpha \sim \mathrm{N}(0, \sigma_{\alpha}^2), \quad \boldsymbol{\beta} \sim \mathrm{N}(\mathbf{0}, \Sigma_{\beta}), \] dan prior untuk parameter presisi \(\tau_u, \tau_v\) yang mengontrol variabilitas \(u_i\) dan \(v_i\).
INLA memanfaatkan fakta bahwa model ini adalah Latent Gaussian Model: latent field \((\eta, u, v, \alpha, \beta)\) berdistribusi Gaussian (conditional pada hyperparameters), sementara likelihood dapat berupa salah satu distribusi dari keluarga yang luas (termasuk Poisson).
8.7.4 Model Spasial: IID, ICAR, BYM, BYM2, Leroux
8.7.4.1 Efek IID (non-spasial)
Efek non-struktural (IID) dimodelkan sebagai:
\[ v_i \mid \tau_v \sim \mathrm{N}(0, \tau_v^{-1}), \quad i = 1, \dots, n, \]
dengan:
- \(\tau_v\) : precision (kebalikan variance). Variance = \(1/\tau_v\).
Efek ini tidak mengandung struktur spasial—setiap wilayah bisa memiliki deviasi sendiri tanpa memperhatikan tetangga.
8.7.4.2 Struktur ICAR (Intrinsic CAR)
Model ICAR adalah pondasi untuk banyak model spasial area-level. Asumsinya:
- Kita punya matriks ketetanggaan \(W = (w_{ij})\) dengan \(w_{ij} = 1\) jika area \(i\) dan \(j\) berbatasan, dan 0 jika tidak.
- Derajat area \(i\) : \(m_i = \sum_{j} w_{ij}\).
Model ICAR mendefinisikan distribusi bersama \(u = (u_1, \dots, u_n)\) melalui full conditional distributions:
\[ u_i \mid u_{-i}, \tau_u \sim \mathrm{N}\left( \frac{1}{m_i}\sum_{j \sim i} u_j,\; \frac{1}{\tau_u m_i} \right), \]
di mana \(j \sim i\) berarti \(w_{ij} = 1\) (tetangga).
Dalam bentuk matriks, prior ICAR dapat ditulis sebagai:
\[ \pi(\mathbf{u} \mid \tau_u) \propto \tau_u^{(n-1)/2} \exp\left( -\frac{\tau_u}{2} \mathbf{u}^\top Q \mathbf{u} \right), \]
dengan matriks presisi:
\[ Q_{ij} = \begin{cases} m_i, & \text{jika } i = j, \\ -1, & \text{jika } i \text{ dan } j \text{ bertetangga}, \\ 0, & \text{lainnya}. \end{cases} \]
Matriks ini singular (rank-deficient) karena model hanya mengidentifikasi perbedaan relatif antar area, bukan level absolut. Oleh karena itu prior ICAR disebut intrinsic.
Untuk melakukan inferensi, perlu diberi constraint, misalnya:
\[ \sum_{i=1}^n u_i = 0. \]
INLA menangani hal ini secara otomatis ketika kita menentukan
model = "icar" atau model = "besag"
(tergantung versi paket).
8.7.4.3 Model BYM klasik (Besag–York–Mollié)
Model BYM menggabungkan:
\[ u_i \sim \text{ICAR}(\tau_u), \quad v_i \sim \mathrm{N}(0, \tau_v^{-1}), \]
sehingga:
\[ \eta_i = \alpha + \mathbf{x}_i^\top\boldsymbol{\beta} + u_i + v_i. \]
Model ini melakukan smoothing risiko relatif baik secara spasial maupun non-spasial.
Masalah klasik BYM:
- Interpretabilitas: sulit memisahkan kontribusi spasial dan non-spasial.
- Identifiability: skala \(u_i\) dan \(v_i\) bergantung pada unit dan struktur graf.
Untuk mengatasi hal ini, dikembangkan reparameterization seperti BYM2.
8.7.4.4 Model BYM2
Model BYM2 menulis efek total sebagai kombinasi terstandar:
\[ b_i = \frac{1}{\sqrt{\tau}} \left( \sqrt{1 - \phi}\, \theta_i + \sqrt{\phi}\, s_i \right), \]
dengan:
- \(\theta_i \sim \mathrm{N}(0,1)\) efek IID terstandar.
- \(s_i\) adalah efek spasial ICAR terstandar (varians marjinal 1).
- \(\tau\) adalah precision global (mengontrol total varians).
- \(\phi \in [0,1]\) adalah parameter yang mengontrol proporsi variasi spasial terstruktur.
Interpretasi:
- \(\phi \approx 0\) → mayoritas variasi adalah non-spasial (IID).
- \(\phi \approx 1\) → variasi hampir seluruhnya spasial.
Keuntungan BYM2:
- Skala varians terkendali (tidak bergantung pada struktur graf).
- Memungkinkan penggunaan prior yang bermakna untuk \(\tau\) dan \(\phi\) (misalnya PC priors).
Dalam INLA, model BYM2 diakses dengan
model = "bym2".
8.7.4.5 Model Leroux
Model Leroux adalah model CAR yang memadukan IID dan CAR dalam satu parameter smoothing \(\rho\):
\[ \mathbf{u} \mid \tau, \rho \sim \mathrm{N}\left(\mathbf{0}, \tau^{-1} \left[ (1-\rho)I + \rho R \right]^{-1} \right), \]
di mana \(R\) adalah matriks presisi ICAR (biasanya \(R = D - W\)).
Parameter \(\rho \in [0,1]\):
- \(\rho = 0\): efek IID.
- \(\rho \to 1\): mendekati ICAR.
Model ini memberikan cara lain untuk mengontrol kekuatan smoothing spasial.
8.7.5 Konsep Dasar INLA
8.7.5.1 Latent Gaussian Model (LGM)
INLA dirancang untuk model dengan struktur:
- Likelihood:
\[ y_i \mid \eta_i, \boldsymbol{\theta} \sim \pi(y_i \mid \eta_i, \boldsymbol{\theta}), \quad i = 1,\dots,n. \]
- Linear predictor:
\[ \eta_i = \mu_i + \sum_k z_{ik} x_k, \]
di mana \(x_k\) adalah komponen dari latent field \(\mathbf{x}\) (koefisien regresi, random effects, dsb).
- Latent field Gaussian:
\[ \mathbf{x} \mid \boldsymbol{\theta} \sim \mathrm{N}(\mathbf{0}, Q(\boldsymbol{\theta})^{-1}), \]
dengan \(Q(\boldsymbol{\theta})\) matriks presisi (jarang/sparse).
- Hyperparameters:
\[ \boldsymbol{\theta} \sim \pi(\boldsymbol{\theta}). \]
Tujuan: menghitung posterior marginal:
- \(\pi(x_k \mid \mathbf{y})\) untuk semua komponen \(x_k\).
- \(\pi(\theta_j \mid \mathbf{y})\) untuk semua hyperparameter.
8.7.5.2 Ide utama INLA
INLA menggunakan kombinasi pendekatan:
- Laplace approximation untuk posterior \(\pi(\boldsymbol{\theta} \mid \mathbf{y})\) dan conditional \(\pi(\mathbf{x} \mid \boldsymbol{\theta}, \mathbf{y})\).
- Numerical integration (grid/Adaptive) pada ruang hyperparameter berdimensi rendah (biasanya 1–5 dimensi).
- Memanfaatkan struktur sparse matriks presisi \(Q(\boldsymbol{\theta})\) untuk efisiensi komputasi.
Skema kasar:
- Aproksimasi posterior hyperparameters:
\[ \tilde{\pi}(\boldsymbol{\theta} \mid \mathbf{y}) \propto \frac{ \pi(\mathbf{x}, \boldsymbol{\theta}, \mathbf{y}) \big|_{\mathbf{x} = \mathbf{x}^*(\boldsymbol{\theta})} }{ \tilde{\pi}_G(\mathbf{x} \mid \boldsymbol{\theta}, \mathbf{y}) \big|_{\mathbf{x} = \mathbf{x}^*(\boldsymbol{\theta})} }, \]
di mana \(\mathbf{x}^*(\boldsymbol{\theta})\) adalah modus conditional \(\pi(\mathbf{x} \mid \boldsymbol{\theta}, \mathbf{y})\) dan \(\tilde{\pi}_G\) adalah pendekatan Gaussian.
- Posterior marginal untuk \(x_k\) dihitung sebagai:
\[ \pi(x_k \mid \mathbf{y}) = \int \pi(x_k \mid \boldsymbol{\theta}, \mathbf{y}) \pi(\boldsymbol{\theta} \mid \mathbf{y}) \, d\boldsymbol{\theta} \approx \sum_{m} \tilde{\pi}(x_k \mid \boldsymbol{\theta}_m, \mathbf{y}) \tilde{\pi}(\boldsymbol{\theta}_m \mid \mathbf{y}) \Delta_m. \]
INLA menyediakan tiga tingkat aproksimasi untuk \(\pi(x_k \mid \boldsymbol{\theta}, \mathbf{y})\):
- Gaussian approximation.
- Simplified Laplace.
- Full Laplace.
8.7.5.3 Keunggulan INLA
- Sangat cepat untuk model dengan struktur GMRF (Gaussian Markov Random Field), termasuk model spasial area-level.
- Menghasilkan posterior marginal langsung, bukan sampel yang harus diringkas.
- Mempermudah eksplorasi model (misal penggantian prior, model alternatif).
8.7.6 Prior dan PC Priors dalam Disease Mapping
Pemilihan prior untuk model disease mapping sangat penting. Prior yang kurang tepat dapat menghasilkan smoothing berlebihan atau under-smoothing.
8.7.6.1 Prior standar untuk precision
Pada model klasik, precision \(\tau\) sering diberikan prior Gamma:
\[ \tau \sim \text{Gamma}(a, b). \]
Namun, interpretasi langsung dari parameter \(a, b\) terhadap varians sering tidak intuitif, dan prior ini bisa terlalu informatif pada skala log-varians.
8.7.6.2 Penalized Complexity (PC) priors
PC priors didesain berdasarkan prinsip bahwa:
- Model sederhana (base model) harus lebih disukai secara a priori.
- Deviasi dari base model diberi penalti secara eksponensial.
Untuk varians atau precision, PC prior sering dinyatakan dalam bentuk:
- Prior untuk standard deviation \(\sigma\) dengan tail behavior:
\[ \Pr(\sigma > U) = \alpha, \]
untuk \(U\) dan \(\alpha\) yang ditentukan pengguna. Ini kemudian diubah menjadi prior untuk \(\tau = 1/\sigma^2\).
Contoh: jika kita ingin menyatakan bahwa a priori ada probabilitas 0.01 bahwa standard deviation melebihi 1, maka:
\[ \Pr(\sigma > 1) = 0.01. \]
INLA menyediakan PC priors untuk banyak komponen, termasuk:
- Standard deviation dari random effect (misal
theta1). - Parameter proporsi spasial \(\phi\) pada BYM2.
8.7.6.3 PC prior untuk BYM2
Pada BYM2:
\[ b_i = \frac{1}{\sqrt{\tau}} \left( \sqrt{1 - \phi}\, \theta_i + \sqrt{\phi}\, s_i \right), \]
dengan hyperparameters:
- \(\tau\) (precision global) → PC prior berdasarkan \(\sigma = 1/\sqrt{\tau}\).
- \(\phi\) (proporsi variasi spasial) → PC prior berdasarkan base model \(\phi = 0\) (tidak ada struktur spasial).
Misalnya, kita dapat mengatur:
- \(\Pr(\sigma > 1) = 0.01\).
- \(\Pr(\phi < 0.5) = 0.5\) (prior median berada sekitar 0.5).
Di R-INLA, ini dilakukan melalui argumen hyper di
f(..., model = "bym2", hyper = list(...)).
8.7.7 Implementasi R-INLA untuk Disease Mapping
8.7.7.1 Persiapan R dan paket
Sebelum memulai, pastikan paket INLA sudah terinstal.
Biasanya, INLA dipasang dari repositori khusus.
# install.packages("INLA", repos=c(getOption("repos"),
# INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
library(readxl)
library(INLA)
library(dplyr)
library(ggplot2)
library(sf) # untuk data spasial vector
library(spdep) # untuk adjacency dan nb2INLA
library(MASS)Catatan: Sesuaikan instalasi INLA dengan dokumentasi versi terbaru.
8.7.7.2 Simulasi data disease mapping sederhana
Kita akan mensimulasikan data pada peta grid sederhana (misalnya \(5 \times 5\) = 25 area) dengan struktur ketetanggaan rook (atas–bawah–kiri–kanan).
8.7.7.2.1 Membuat grid dan adjacency
set.seed(123)
# Buat grid 5x5
nrow <- 5
ncol <- 5
n <- nrow * ncol
grid <- expand.grid(
col = 1:ncol,
row = 1:nrow
) %>%
mutate(id = 1:n)
# Buat objek sf titik (hanya untuk contoh)
coords <- cbind(grid$col, grid$row)
grid_sf <- st_as_sf(data.frame(grid), coords = c("col", "row"), crs = 4326)
# Buat neighbor list rook (4-neighbors)
nb <- spdep::cell2nb(nrow, ncol, type = "rook")
W <- spdep::nb2mat(nb, style = "B") # matriks ketetanggaan8.7.7.2.2 Membuat efek spasial terstruktur (ICAR) dan non-spasial
# Fungsi sederhana untuk generate ICAR dari Q
sim_icar <- function(Q, tau = 1) {
Q_eps <- Q + diag(1e-05, nrow(Q)) # tambahkan ridge
Sigma <- solve(Q_eps)
z <- MASS::mvrnorm(1, mu = rep(0, nrow(Q)), Sigma = (1 / tau) * Sigma)
as.numeric(z - mean(z)) # pusatkan ke mean 0
}
# Matriks presisi ICAR Q = D - W
D <- diag(rowSums(W))
Q <- D - W
tau_u <- 3 # precision efek spasial
tau_v <- 10 # precision efek iid
u_true <- sim_icar(Q, tau = tau_u)
v_true <- rnorm(n, 0, sd = 1 / sqrt(tau_v))8.7.7.2.3 Kovariat dan risiko relatif
x <- rnorm(n, mean = 0, sd = 1)
alpha_true <- 1.0
beta_true <- 0.4
eta_true <- alpha_true + beta_true * x + u_true + v_true
# Eksposur (expected counts)
E <- runif(n, min = 50, max = 500)
lambda_true <- exp(eta_true)
mu_true <- E * lambda_true
# Simulasi kasus
y <- rpois(n, lambda = mu_true)
dat <- data.frame(
id = 1:n,
x = x,
E = E,
y = y,
u_true = u_true,
v_true = v_true
)
head(dat)## id x E y u_true v_true
## 1 1 0.25331851 431.35392 369 -0.73764676 -0.53337926
## 2 2 -0.02854676 273.88727 935 -0.02755701 0.26493153
## 3 3 -0.04287046 224.55906 556 -0.13238557 0.04850084
## 4 4 1.36860228 160.90205 468 -0.09587661 -0.35991050
## 5 5 -0.22577099 99.99341 405 0.08289763 0.39649109
## 6 6 1.51647060 225.49750 886 -0.35134275 0.134859838.7.7.3 Menyiapkan struktur INLA: adjacency dan index
# Indeks area untuk INLA
dat$idx_area <- dat$id
# Buat file adjacency untuk INLA
graph_basename <- tempfile(pattern = "grid_adj")
file_adj <- paste0(graph_basename, ".graph")
nb2INLA(file_adj, nb) # pastikan nama file lengkap termasuk .graph
# Pengecekan isi file
first_lines <- readLines(file_adj, n = 5)
print(first_lines)## [1] "25" "1 2 2 6" "2 3 1 3 7" "3 3 2 4 8" "4 3 3 5 9"# Baca graph
g <- inla.read.graph(file_adj)
# Cek hasil
print(g)## $n
## [1] 25
##
## $nnbs
## [1] 2 3 3 3 2 3 4 4 4 3 3 4 4 4 3 3 4 4 4 3 2 3 3 3 2
##
## $nbs
## $nbs[[1]]
## [1] 2 6
##
## $nbs[[2]]
## [1] 1 3 7
##
## $nbs[[3]]
## [1] 2 4 8
##
## $nbs[[4]]
## [1] 3 5 9
##
## $nbs[[5]]
## [1] 4 10
##
## $nbs[[6]]
## [1] 1 7 11
##
## $nbs[[7]]
## [1] 2 6 8 12
##
## $nbs[[8]]
## [1] 3 7 9 13
##
## $nbs[[9]]
## [1] 4 8 10 14
##
## $nbs[[10]]
## [1] 5 9 15
##
## $nbs[[11]]
## [1] 6 12 16
##
## $nbs[[12]]
## [1] 7 11 13 17
##
## $nbs[[13]]
## [1] 8 12 14 18
##
## $nbs[[14]]
## [1] 9 13 15 19
##
## $nbs[[15]]
## [1] 10 14 20
##
## $nbs[[16]]
## [1] 11 17 21
##
## $nbs[[17]]
## [1] 12 16 18 22
##
## $nbs[[18]]
## [1] 13 17 19 23
##
## $nbs[[19]]
## [1] 14 18 20 24
##
## $nbs[[20]]
## [1] 15 19 25
##
## $nbs[[21]]
## [1] 16 22
##
## $nbs[[22]]
## [1] 17 21 23
##
## $nbs[[23]]
## [1] 18 22 24
##
## $nbs[[24]]
## [1] 19 23 25
##
## $nbs[[25]]
## [1] 20 24
##
##
## $cc
## $cc$id
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $cc$n
## [1] 1
##
## $cc$nodes
## $cc$nodes[[1]]
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
##
##
## $cc$mean
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Levels: 1
##
##
## attr(,"class")
## [1] "inla.graph"8.7.7.4 Model BYM klasik dengan INLA
formula_bym <- y ~ 1 + x +
f(idx_area, model = "bym",
graph = g,
scale.model = TRUE,
hyper = list(
prec.unstruct = list(prior = "pc.prec", param = c(0.5, 0.01)),
prec.spatial = list(prior = "pc.prec", param = c(0.5, 0.01))
))
res_bym <- inla(
formula_bym,
family = "poisson",
data = dat,
E = E,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
summary(res_bym)## Time used:
## Pre = 1.2, Running = 2.11, Post = 0.0213, Total = 3.33
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.039 0.037 0.962 1.039 1.116 1.039 388.332
## x 0.455 0.075 0.306 0.455 0.604 0.455 52.568
##
## Random effects:
## Name Model
## idx_area BYM model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for idx_area (iid component) 7864.32 1.46e+05 4.06 152.45
## Precision for idx_area (spatial component) 5.25 2.53e+00 1.58 4.81
## 0.975quant mode
## Precision for idx_area (iid component) 35662.48 6.73
## Precision for idx_area (spatial component) 11.24 3.88
##
## Deviance Information Criterion (DIC) ...............: 258.56
## Deviance Information Criterion (DIC, saturated) ....: 49.74
## Effective number of parameters .....................: 24.70
##
## Watanabe-Akaike information criterion (WAIC) ...: 251.32
## Effective number of parameters .................: 12.58
##
## Marginal log-Likelihood: -168.27
## CPO, PIT is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')8.7.7.5 Model BYM2 (reparameterized)
formula_bym2 <- y ~ 1 + x +
f(idx_area, model = "bym2",
graph = g,
scale.model = TRUE,
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01)),
phi = list(prior = "pc", param = c(0.5, 0.5))
))
res_bym2 <- inla(
formula_bym2,
family = "poisson",
data = dat,
E = E,
control.predictor = list(compute = TRUE),
control.compute = list(
dic = TRUE,
waic = TRUE,
cpo = TRUE,
return.marginals.predictor = TRUE # **ini tambahan penting**
)
)
summary(res_bym2)## Time used:
## Pre = 8.26, Running = 0.227, Post = 0.016, Total = 8.5
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.039 0.054 0.925 1.039 1.152 1.039 177.255
## x 0.453 0.080 0.296 0.454 0.610 0.454 47.076
##
## Random effects:
## Name Model
## idx_area BYM2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for idx_area 5.602 1.55 3.103 5.420 9.144 5.091
## Phi for idx_area 0.617 0.22 0.174 0.641 0.955 0.814
##
## Deviance Information Criterion (DIC) ...............: 258.57
## Deviance Information Criterion (DIC, saturated) ....: 49.75
## Effective number of parameters .....................: 24.72
##
## Watanabe-Akaike information criterion (WAIC) ...: 251.30
## Effective number of parameters .................: 12.57
##
## Marginal log-Likelihood: -170.02
## CPO, PIT is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')8.7.7.5.1 Posterior untuk phi
phi_marg <- res_bym2$marginals.hyperpar[["Phi for idx_area"]]
# Cetak beberapa elemen untuk melihat struktur
str(phi_marg)## num [1:43, 1:2] 0.0191 0.033 0.0612 0.1256 0.1736 ...
## - attr(*, "hyperid")= chr "11002|idx_area"
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "x" "y"head(phi_marg, 3)## x y
## [1,] 0.01910488 0.002186451
## [2,] 0.03301109 0.011683780
## [3,] 0.06120761 0.056509568phi_df <- data.frame(
phi = phi_marg[, 1],
dens = phi_marg[, 2]
)
ggplot(phi_df, aes(phi, dens)) +
geom_line() +
labs(x = expression(phi), y = "Posterior density",
title = "Posterior untuk proporsi variasi spasial (phi)")
8.7.7.6 Menghasilkan peta risiko relatif dari simulasi grid
Pada contoh simulasi grid, kita belum punya polygon area, hanya titik. Secara praktis, untuk analisis nyata, kita akan memakai shapefile polygon. Namun, kita bisa tetap mem-plot nilai RR (Relative Risk) sebagai heatmap titik.
eta_mean <- res_bym2$summary.linear.predictor$mean
lambda_mean <- exp(eta_mean)
dat$lambda_mean <- lambda_mean
# Gabungkan dengan koordinat grid
grid_plot <- cbind(grid, dat["lambda_mean"])
ggplot(grid_plot, aes(x = col, y = row, fill = lambda_mean)) +
geom_tile(color = "grey50") +
scale_fill_viridis_c(option = "C") +
coord_equal() +
labs(fill = "RR (mean)",
title = "Posterior mean Relative Risk (BYM2) pada Grid 5x5") +
theme_minimal()
8.7.7.7 Exceedance Probability
# Setelah fitting res_bym2
# Pastikan marginals tersedia:
# control.compute = list(..., return.marginals.fitted.values = TRUE) was set.
# Atau jika tidak, bisa gunakan marginals.linear.predictor dan transformasi.
# Misalkan threshold
c_threshold <- 1.2
# Ambil marginals untuk fitted values (λ)
marg_list <- res_bym2$marginals.fitted.values
if (is.null(marg_list)) {
# Jika marginals.fitted.values NULL, gunakan marginals.linear.predictor dan transformasi
marg_lp <- res_bym2$marginals.linear.predictor
marg_list <- lapply(marg_lp,
function(m) inla.tmarginal(fun = exp, m = m))
}
# Hitung exceedance probability
exceed_prob <- sapply(marg_list, function(marg) {
1 - inla.pmarginal(q = c_threshold, marginal = marg)
})
# Tambahkan ke dataset
dat$exceed_prob <- exceed_prob
# Gabungkan ke grid untuk plotting
grid_plot2 <- cbind(grid, dat["exceed_prob"])
library(ggplot2)
ggplot(grid_plot2, aes(x = col, y = row, fill = exceed_prob)) +
geom_tile(color = "grey50") +
scale_fill_viridis_c(option = "C", labels = scales::percent_format(accuracy = 1)) +
coord_equal() +
labs(fill = paste0("Prob(λ > ", c_threshold, ")"),
title = paste0("Peta Probabilitas Melebihi RR > ", c_threshold)) +
theme_minimal()
8.7.8 Contoh dengan Peta Area Nyata (Sketsa Konseptual)
Bagian ini sketsa konseptual untuk data nyata; .
library(sf)
library(dplyr)
library(spdep)
library(INLA)
library(ggplot2)
# Jika file ada di folder & nama misalnya "Bandung.shp"
setwd("/Users/mindra/mybook/Epidemiologi/BANDUNG")
Bandung_sf <- st_read("Bandung.shp")## Reading layer `BANDUNG' from data source
## `/Users/mindra/mybook/Epidemiologi/BANDUNG/BANDUNG.shp' using driver `ESRI Shapefile'
## Simple feature collection with 30 features and 12 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 107.5453 ymin: -6.969748 xmax: 107.7395 ymax: -6.840746
## Geodetic CRS: WGS84 Lat/Long's, Degrees, -180 ==> +180Bandung_sf$id<-c(1:30)
Bandung_sf<-st_as_sf(Bandung_sf) # Pastikan kolom geometry ada
df_case<-read_xlsx("DataDBDRainfall.xlsx")
dat2 <- Bandung_sf %>%
st_drop_geometry() %>%
left_join(df_case, by = "id")
# 2. Buat neighbor list & graph untuk INLA
nb2 <- poly2nb(Bandung_sf) # opsi ‘row.names’ agar matching lebih aman
# Cek jumlah area
n_area <- length(nb2)
graph_basename2 <- tempfile(pattern = "kec_adj")
file_graph2 <- paste0(graph_basename2, ".graph")
# Buat graph file secara eksplisit dengan ekstensi
spdep::nb2INLA(file = file_graph2, nb = nb2)
# Cek file
if (!file.exists(file_graph2)) {
stop("File .graph tidak ditemukan: ", file_graph2)
}
# Tampilkan beberapa baris awal dari file, untuk verifikasi
first_lines <- readLines(file_graph2, n = 5)
print(first_lines)## [1] "30" "1 4 3 25 29 30" "2 5 3 21 22 25 30"
## [4] "3 4 1 2 25 30" "4 2 5 6"# Baca graph
g2 <- INLA::inla.read.graph(file_graph2)
# Verifikasi objek graph
print(g2$n) # harus = number of areas## [1] 30str(g2)## List of 4
## $ n : int 30
## $ nnbs: num [1:30] 4 5 4 2 7 4 4 4 3 5 ...
## $ nbs :List of 30
## ..$ : int [1:4] 3 25 29 30
## ..$ : int [1:5] 3 21 22 25 30
## ..$ : int [1:4] 1 2 25 30
## ..$ : int [1:2] 5 6
## ..$ : int [1:7] 4 6 7 11 13 16 18
## ..$ : int [1:4] 4 5 7 23
## ..$ : int [1:4] 5 6 18 23
## ..$ : int [1:4] 9 10 11 28
## ..$ : int [1:3] 8 10 28
## ..$ : int [1:5] 8 9 11 12 24
## ..$ : int [1:5] 5 8 10 12 13
## ..$ : int [1:6] 10 11 13 14 15 24
## ..$ : int [1:7] 5 11 12 14 16 18 27
## ..$ : int [1:7] 12 13 15 17 19 20 27
## ..$ : int [1:5] 12 14 19 21 24
## ..$ : int [1:5] 5 13 17 18 27
## ..$ : int [1:6] 14 16 18 20 23 27
## ..$ : int [1:6] 5 7 13 16 17 23
## ..$ : int [1:5] 14 15 20 21 26
## ..$ : int [1:6] 14 17 19 26 29 30
## ..$ : int [1:5] 2 15 19 26 30
## ..$ : int [1:2] 2 25
## ..$ : int [1:4] 6 7 17 18
## ..$ : int [1:3] 10 12 15
## ..$ : int [1:4] 1 2 3 22
## ..$ : int [1:4] 19 20 21 30
## ..$ : int [1:4] 13 14 16 17
## ..$ : int [1:2] 8 9
## ..$ : int [1:3] 1 20 30
## ..$ : int [1:7] 1 2 3 20 21 26 29
## $ cc :List of 4
## ..$ id : int [1:30] 1 1 1 1 1 1 1 1 1 1 ...
## ..$ n : int 1
## ..$ nodes:List of 1
## .. ..$ : int [1:30] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ mean : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "class")= chr "inla.graph"# 3. Tambahkan index area untuk random effect
dat2 <- dat2 %>%
mutate(
idx_area = as.integer(factor(id)) # 1…n_area
)
stopifnot(max(dat2$idx_area) == n_area)
# 4. Definisikan formula dengan BYM2 + kovariat
formula_kec <- Cases ~ 1 + Rainfall +
f(idx_area, model = "bym2",
graph = g2,
scale.model = TRUE,
constr = TRUE,
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01)),
phi = list(prior = "pc", param = c(0.5, 0.5))
))
# 5. Jalankan INLA
res_kec <- inla(
formula = formula_kec,
family = "poisson",
data = dat2,
E = dat2$E,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
# 6. Ekstrak RR
dat2 <- dat2 %>%
mutate(
eta_mean = res_kec$summary.linear.predictor$mean,
RR_mean = exp(eta_mean)
)
# 7. Gabungkan kembali ke sf untuk plot
library(dplyr)
Bandung_res <- Bandung_sf %>%
left_join(dat2, by = "id")
Bandung_res_2024<-Bandung_res[Bandung_res$Year==2024,]
ggplot(Bandung_res_2024) +
geom_sf(aes(fill = RR_mean)) +
facet_wrap(~Month)+
scale_fill_viridis_c(option = "C") +
labs(
title = "Posterior Mean Relative Risk per Bandungan",
fill = "RR"
) +
theme_minimal()
8.7.9 Perbandingan Model: IID-only vs BYM2
# —————
# Perbandingan Model: IID‐only vs BYM2
# Pastikan data model IID vs BYM2 menggunakan *data yang sama (dat2)*
formula_iid <- Cases ~ 1 +Rainfall +
f(idx_area, model = "iid",
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01))
))
res_iid <- inla(
formula = formula_iid,
family = "poisson",
data = dat2,
E = dat2$E,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
# Ambil DIC/WAIC
model_comp <- tibble(
model = c("IID_only", "BYM2"),
DIC = c(res_iid$dic$dic, res_kec$dic$dic),
WAIC = c(res_iid$waic$waic, res_kec$waic$waic)
)
print(model_comp)## # A tibble: 2 × 3
## model DIC WAIC
## <chr> <dbl> <dbl>
## 1 IID_only 26575. 26267.
## 2 BYM2 26575. 26267.dat2 <- dat2 %>%
mutate(
RR_iid = exp(res_iid$summary.linear.predictor$mean),
RR_bym2 = RR_mean
)
rr_long <- dat2 %>%
dplyr::select(id, RR_iid, RR_bym2) %>%
tidyr::pivot_longer(
cols = starts_with("RR_"),
names_to = "model",
values_to= "RR"
)
library(ggplot2)
ggplot(dat2, aes(x = RR_iid, y = RR_bym2)) +
geom_point(alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey30") +
labs(
x = "RR (IID-only)",
y = "RR (BYM2)",
title = "Perbandingan Relative Risk: IID-only vs BYM2"
) +
theme_minimal()
8.7.10 Garis Besar Ekstensi Spatio-Temporal
# Misal kita punya data panel area-time:
# data_st: id_area, time, y, E, x, ...
data_st <- df_case %>%
mutate(
area_idx = as.integer(factor(id)),
time_idx = as.integer(factor(Month)),
area_time_idx = as.integer(interaction(area_idx, time_idx, drop = TRUE))
)
formula_st <- Cases ~ 1 + Rainfall +
f(area_idx, model = "bym2", graph = g2, scale.model = TRUE, constr = TRUE,
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01)),
phi = list(prior = "pc", param = c(0.5, 0.5))
)) +
f(time_idx, model = "rw1",
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01))
)) +
f(area_time_idx, model = "iid",
hyper = list(
prec = list(prior = "pc.prec", param = c(0.5, 0.01))
))
res_st <- inla(
formula = formula_st,
family = "poisson",
data = data_st,
E = data_st$E,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)Penutup
Materi ini telah mengulas:
- Formulasi dasar disease mapping Bayesian (Poisson log-linear).
- Model spasial area-level: ICAR, BYM, BYM2, dan Leroux.
- Konsep Latent Gaussian Model dan cara kerja INLA (Laplace approximation dan integrasi nested).
- PC priors dan bagaimana menggunakannya dalam konteks disease mapping.
- Implementasi praktis di R-INLA: dari simulasi data, fitting model BYM/BYM2, hingga visualisasi relative risk dan perbandingan model.
- Sketsa singkat perluasan ke model spatio-temporal.
Dengan R Markdown ini, Anda dapat:
- Menggunakannya sebagai modul kuliah atau e-book mini untuk materi Bayesian disease mapping.
- Memodifikasi bagian simulasi menjadi studi kasus real (misalnya data DBD atau TB per kecamatan).
- Menambahkan latihan untuk mahasiswa: misalnya, eksplorasi prior, membandingkan model, dan interpretasi peta risiko.
Selanjutnya, Anda bisa mengintegrasikan materi ini dengan topik lain seperti:
- Model multivariate disease mapping (beberapa penyakit sekaligus).
- Penggunaan SPDE–GMRF untuk spatial smoothing di domain kontinu.
- Integrasi faktor lingkungan (polusi udara, iklim) dalam model risiko penyakit.
9 Analisis Survival
9.1 Konsep Dasar Analisis Survival
Analisis survival (juga disebut time-to-event analysis) mempelajari waktu sampai terjadinya suatu kejadian (event).
Contoh event: - Kematian pasien - Kekambuhan penyakit - Kegagalan mesin - Drop-out dari program
9.1.1 Notasi Dasar
Misalkan: - \(T\) = peubah acak waktu sampai event terjadi, dengan \(T \ge 0\) - \(t\) = realisasi waktu (angka tertentu, misal 12 bulan)
Fungsi-fungsi utama dalam analisis survival:
Fungsi distribusi kumulatif: \[ F(t) = P(T \le t) \] Probabilitas event sudah terjadi paling lambat pada waktu \(t\).
Fungsi survival: \[ S(t) = P(T > t) = 1 - F(t) \] Probabilitas subjek masih bertahan (belum mengalami event) sampai lewat waktu \(t\).
Fungsi densitas peluang (jika kontinu): \[ f(t) = rac{d}{dt}F(t) \]
Fungsi hazard: \[ h(t) = \lim_{\Delta t o 0} rac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} \]
Secara intuitif, \(h(t)\) adalah laju risiko sesaat di waktu \(t\), diberikan bahwa subjek masih hidup/bertahan sampai tepat sebelum \(t\).
Hubungan penting: \[ h(t) = rac{f(t)}{S(t)} \]
9.1.2 Fungsi Hazard Kumulatif
Definisikan hazard kumulatif: \[ H(t) = \int_0^t h(u)\,du \]
Hubungan antara \(S(t)\) dan \(H(t)\): \[ S(t) = \exp\{-H(t)\} \]
Ini sering dipakai dalam teori dan estimasi model survival.
9.1.3 Data Tersensor (Censoring)
Keunikan analisis survival adalah adanya censoring:
- Right censoring: kita tahu subjek belum mengalami
event sampai waktu tertentu, tapi setelah itu tidak tahu.
- Contoh: studi selesai sebelum pasien meninggal.
- Left censoring: event terjadi sebelum waktu observasi awal.
- Interval censoring: hanya tahu bahwa event terjadi di antara dua waktu.
Kasus paling umum dalam praktik: right-censored.
Representasi data: - \(T_i = \min(T_i^st,
C_i)\)
\(T_i^st\) = waktu event
sebenarnya,
\(C_i\) = waktu censoring (misal akhir
studi), - \(\delta_i = 1\) jika event
teramati, \(\delta_i = 0\) jika
tersensor.
9.2 Estimasi Nonparametrik: Kaplan–Meier & Nelson–Aalen
9.2.1 Estimator Kaplan–Meier (Product-Limit)
Kaplan–Meier digunakan untuk mengestimasi fungsi survival \(S(t)\) tanpa asumsi bentuk distribusi.
Misalkan: - Waktu event terurut: \(t_{(1)}
< t_{(2)} < \dots < t_{(k)}\)
(hanya waktu di mana terjadi event, bukan censoring) - \(d_j\) = jumlah event pada waktu \(t_{(j)}\) - \(n_j\) = jumlah individu yang masih
berisiko tepat sebelum \(t_{(j)}\)
Estimator Kaplan–Meier dituliskan sebagai:
\[ \hat{S}(t) = \prod_{t_{(j)} \le t} \left( 1 - \frac{d_j}{n_j} \right) \]
dengan:
- \(d_j\): jumlah event pada waktu \(t_{(j)}\)
- \(n_j\): jumlah individu yang masih berisiko tepat sebelum waktu \(t_{(j)}\) Ciri penting:
- Fungsi survival empiris berbentuk step function (tangga).
- Bisa digambar dengan confidence interval.
9.2.2 Contoh Perhitungan Manual Kaplan–Meier (10 Pasien Kanker)
Contoh ini menggunakan data 10 pasien kanker yang dipantau selama
beberapa bulan.
Status: - 1 = meninggal (event) - 0 =
tersensor (misal hilang follow up / studi berakhir)
9.2.2.1 Data
kaplan_data <- data.frame(
id = 1:10,
waktu = c(3, 4, 5, 6, 6, 10, 15, 15, 16, 28),
status = c(1, 0, 1, 1, 1, 0, 1, 1, 1, 1)
)
kaplan_data## id waktu status
## 1 1 3 1
## 2 2 4 0
## 3 3 5 1
## 4 4 6 1
## 5 5 6 1
## 6 6 10 0
## 7 7 15 1
## 8 8 15 1
## 9 9 16 1
## 10 10 28 1Total ada 10 pasien pada awal studi.
9.2.2.2 Langkah 1: Waktu event dan censoring
- Waktu event (status = 1):
\(3, 5, 6, 6, 15, 15, 16, 28\) - Waktu censoring (status = 0):
\(4, 10\)
Untuk Kaplan–Meier, kita fokus pada waktu event
unik:
\[
3, 5, 6, 15, 16, 28
\]
9.2.2.3 Langkah 2: Hitung \(n_j\) dan \(d_j\)
Definisi: - \(n_j\): banyaknya individu yang masih berisiko tepat sebelum \(t_{(j)}\) - \(d_j\): banyaknya event pada waktu \(t_{(j)}\)
Secara manual, dapat diringkas sebagai berikut:
| Waktu | \(n_j\) | \(d_j\) | Keterangan singkat |
|---|---|---|---|
| 3 | 10 | 1 | 10 berisiko, 1 meninggal |
| 5 | 8 | 1 | 1 dicensor di 4, tersisa 8 berisiko, 1 meninggal di 5 |
| 6 | 7 | 2 | 7 berisiko, 2 meninggal |
| 15 | 4 | 2 | 1 dicensor di 10, 4 berisiko, 2 meninggal |
| 16 | 2 | 1 | 2 berisiko, 1 meninggal |
| 28 | 1 | 1 | 1 berisiko, 1 meninggal |
9.2.2.4 Langkah 3: Hitung faktor \((1 - d_j/n_j)\)
\[ \begin{aligned} t_{(1)} = 3&: && \frac{d_1}{n_1} = \frac{1}{10} = 0.10, && 1 - \frac{d_1}{n_1} = 0.90 \\ t_{(2)} = 5&: && \frac{d_2}{n_2} = \frac{1}{8} = 0.125, && 1 - \frac{d_2}{n_2} = 0.875 \\ t_{(3)} = 6&: && \frac{d_3}{n_3} = \frac{2}{7} \approx 0.28571, && 1 - \frac{d_3}{n_3} \approx 0.71429 \\ t_{(4)} = 15&: && \frac{d_4}{n_4} = \frac{2}{4} = 0.50, && 1 - \frac{d_4}{n_4} = 0.50 \\ t_{(5)} = 16&: && \frac{d_5}{n_5} = \frac{1}{2} = 0.50, && 1 - \frac{d_5}{n_5} = 0.50 \\ t_{(6)} = 28&: && \frac{d_6}{n_6} = \frac{1}{1} = 1.00, && 1 - \frac{d_6}{n_6} = 0 \end{aligned} \]
9.2.2.5 Langkah 4: Hitung \(\hat{S}(t)\) sebagai produk kumulatif
Rumus Kaplan–Meier:
\[ \hat{S}(t) = \prod_{t_{(j)} \le t} \left(1 - \frac{d_j}{n_j}\right) \]
Secara bertahap:
\[ \begin{aligned} \hat{S}(0) &= 1 \\ \hat{S}(3) &= 1 \times 0.90 = 0.90000 \\ \hat{S}(5) &= 0.90000 \times 0.875 = 0.78750 \\ \hat{S}(6) &= 0.78750 \times 0.71429 \approx 0.56250 \\ \hat{S}(15) &= 0.56250 \times 0.50 = 0.28125 \\ \hat{S}(16) &= 0.28125 \times 0.50 = 0.14063 \\ \hat{S}(28) &= 0.14063 \times 0 = 0 \end{aligned} \]
Ringkasan:
| Waktu | \(n_j\) | \(d_j\) | Hazard \(d_j/n_j\) | \(1 - d_j/n_j\) | \(\hat{S}(t)\) |
|---|---|---|---|---|---|
| 0 | 10 | 0 | 0.00000 | 1.00000 | 1.00000 |
| 3 | 10 | 1 | 0.10000 | 0.90000 | 0.90000 |
| 5 | 8 | 1 | 0.12500 | 0.87500 | 0.78750 |
| 6 | 7 | 2 | 0.28571 | 0.71429 | 0.56250 |
| 15 | 4 | 2 | 0.50000 | 0.50000 | 0.28125 |
| 16 | 2 | 1 | 0.50000 | 0.50000 | 0.14063 |
| 28 | 1 | 1 | 1.00000 | 0.00000 | 0.00000 |
Censoring pada waktu 4 dan 10 tidak muncul sebagai baris event, tetapi memengaruhi nilai \(n_j\) berikutnya.
9.2.2.6 Cek dengan Fungsi
survfit di R
km_fit <- survfit(Surv(waktu, status) ~ 1, data = kaplan_data)
summary(km_fit)## Call: survfit(formula = Surv(waktu, status) ~ 1, data = kaplan_data)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3 10 1 0.900 0.0949 0.7320 1.000
## 5 8 1 0.787 0.1340 0.5641 1.000
## 6 7 2 0.562 0.1651 0.3165 1.000
## 15 4 2 0.281 0.1631 0.0903 0.876
## 16 2 1 0.141 0.1286 0.0234 0.844
## 28 1 1 0.000 NaN NA NAPerhatikan kolom time dan surv pada output
summary(km_fit) – nilai surv seharusnya sama
(atau sangat dekat) dengan nilai \(\hat{S}(t)\) hasil perhitungan manual.
9.2.2.7 Plot Kurva Kaplan–Meier
ggsurvplot(
km_fit,
conf.int = TRUE,
risk.table = TRUE,
xlab = "Waktu (bulan)",
ylab = "Probabilitas survival",
surv.median.line = "hv"
)
Kurva survival Kaplan–Meier akan berbentuk tangga, dengan penurunan di setiap waktu event, dan risk table menunjukkan berapa banyak subjek yang masih berisiko pada berbagai waktu.
9.2.3 Estimator Hazard Kumulatif: Nelson–Aalen
Nelson–Aalen mengestimasi fungsi hazard kumulatif \(H(t)\) dengan:
\[ \hat{H}(t) = \sum_{t_{(j)} \le t} \frac{d_j}{n_j} \]
Lalu fungsi survival dapat didekati dengan:
\[ \hat{S}(t) \approx \exp\{-\hat{H}(t)\} \] —
9.3 Uji Log Rank
Dokumen ini menjelaskan perhitungan manual uji
log-rank (perbandingan dua kurva Kaplan–Meier) menggunakan
contoh data ovarian dari paket
survival di R.
Variabel utama:
futime= waktu follow-up (misal dalam hari)fustat= status kejadian1= event terjadi (misal: meninggal)
0= tersensor (masih hidup / lost follow-up)
rx= kelompok perlakuan1= treatment 1
2= treatment 2
Uji log-rank digunakan untuk menguji:
H_0: tidak ada perbedaan fungsi survival antara dua kelompok
H_1: ada perbedaan fungsi survival antara minimal dua kelompok
9.3.1 Uji Log-Rank: Formulasi Matematis Lengkap
Sekarang kita masuk ke rumus log-rank seperti yang muncul di gambar Excel.
9.3.1.1 Notasi 2 Kelompok
Misal ada dua kelompok:
- Kelompok 1: misalnya
rx = 1 - Kelompok 2: misalnya
rx = 2
Untuk setiap waktu event \(t_j\) (j = 1, 2, …, k), didefinisikan:
- \(Y_{1j}\) = jumlah subjek kelompok 1 yang masih at risk tepat sebelum \(t_j\)
- \(Y_{2j}\) = jumlah subjek kelompok 2 yang masih at risk tepat sebelum \(t_j\)
- \(Y_j = Y_{1j} + Y_{2j}\) = total at risk
- \(d_{1j}\) = jumlah event pada kelompok 1 pada waktu \(t_j\)
- \(d_{2j}\) = jumlah event pada kelompok 2 pada waktu \(t_j\)
- \(d_j = d_{1j} + d_{2j}\) = total event pada waktu \(t_j\)
Dalam tabel Excel:
- \(d_{1j}\) muncul sebagai Observed untuk kelompok 1 pada baris \(t_j\).
- \(Y_{1j} - d_{1j}\) adalah Not Observed (masih hidup / belum event pada waktu itu).
- \(Y_{1j}\) adalah M.Sum untuk kelompok 1.
9.3.1.2 Ekspektasi (Expected) Event
Di bawah hipotesis nol (H_0) “tidak ada perbedaan survival antar kelompok”, probabilitas gagal (event) pada \(t_j\) dianggap sama untuk semua subjek di risk set.
Probabilitas gagal untuk satu subjek (tanpa membedakan kelompok) pada \(t_j\) diperkirakan oleh:
\[ \widehat{p}_j = \frac{d_j}{Y_j} \]
Sehingga expected number of failures di kelompok 1 pada waktu \(t_j\) adalah:
\[ E_{1j} = Y_{1j} \widehat{p}_j = \frac{Y_{1j}}{Y_j} d_j \]
Begitu pula untuk kelompok 2:
\[ E_{2j} = Y_{2j} \widehat{p}_j = \frac{Y_{2j}}{Y_j} d_j \]
Jelas bahwa:
\[ E_{1j} + E_{2j} = d_j \]
Total expected event untuk kelompok 1 sepanjang semua waktu event:
\[ E_1 = \sum_{j=1}^k E_{1j} \]
Total observed event di kelompok 1:
\[ O_1 = \sum_{j=1}^k d_{1j} \]
9.3.1.3 Varians Komponen (V_j)
Jika kita menganggap \(d_j\) mengikuti distribusi hiper-geometrik (atau binomial aproksimasi) dengan alokasi ke dua kelompok secara acak menurut ukuran risk set, maka varians dari selisih \(d_{1j} - E_{1j}\) adalah:
\[ V_j = \frac{Y_{1j} Y_{2j} d_j (Y_j - d_j)}{Y_j^2 (Y_j - 1)} \]
Intuisi:
- \(Y_{1j} Y_{2j}\) → “ukuran” kedua kelompok
- \(d_j (Y_j - d_j)\) → variasi karena ada yang gagal dan ada yang tidak
- Penyebut \(Y_j^2 (Y_j - 1)\) → penyesuaian skala dan finite population correction
Total varians untuk perbandingan dua kelompok:
\[ V = \sum_{j=1}^k V_j \]
9.3.1.4 Statistik Uji Log-Rank
Statistik uji log-rank fokus pada selisih kumulatif antara observed dan expected di satu kelompok (misalnya kelompok 1):
\[ Z = \frac{O_1 - E_1}{\sqrt{V}} \]
Di bawah H_0, secara asimtotik:
\[ Z \sim N(0, 1) \]
Sehingga kita dapat juga menggunakan bentuk chi-square:
\[ \chi^2 = Z^2 = \frac{(O_1 - E_1)^2}{V} \]
Dengan derajat bebas 1. p-value dihitung sebagai:
\[ p\text{-value} = P\big(\chi^2\_{(1)} \ge \chi^2\_{\text{hitung}}\big) \]
Dalam R:
p_value <- 1 - pchisq(chisq_manual, df = 1)9.3.1.5 Ringkasan Rumus (versi tabel, 2 kelompok)
Untuk memudahkan pengajaran, berikut ringkasan seluruh rumus dalam bentuk tabel:
- Risk set dan kejadian per waktu event \(t_j\)
\[ \begin{aligned} Y_{1j} &= \text{jumlah at risk di grup 1 sebelum } t_j \\ Y_{2j} &= \text{jumlah at risk di grup 2 sebelum } t_j \\ Y_j &= Y_{1j} + Y_{2j} \\ d_{1j} &= \text{jumlah event di grup 1 pada } t_j \\ d_{2j} &= \text{jumlah event di grup 2 pada } t_j \\ d_j &= d_{1j} + d_{2j} \end{aligned} \]
- Ekspektasi dan varians
\[ E_{1j} = \frac{Y_{1j}}{Y_j} d_j, \qquad V_j = \frac{Y_{1j} Y_{2j} d_j (Y_j - d_j)}{Y_j^2 (Y_j - 1)} \]
- Akumulasi
\[ O_1 = \sum_j d_{1j}, \quad E_1 = \sum_j E_{1j}, \quad V = \sum_j V_j \]
- Statistik uji
\[ \chi^2 = \frac{(O_1 - E_1)^2}{V} \sim \chi^2\_{(1)} \quad \text{di bawah } H_0. \]
9.3.1.6 Data dan Pemilihan Sampel
data(ovarian)
ovarian <- ovarian
# Lihat struktur data
str(ovarian)## 'data.frame': 26 obs. of 6 variables:
## $ futime : num 59 115 156 421 431 448 464 475 477 563 ...
## $ fustat : num 1 1 1 0 1 0 1 1 0 1 ...
## $ age : num 72.3 74.5 66.5 53.4 50.3 ...
## $ resid.ds: num 2 2 2 2 2 1 2 2 2 1 ...
## $ rx : num 1 1 1 2 1 1 2 2 1 2 ...
## $ ecog.ps : num 1 1 2 1 1 2 2 2 1 2 ...# Ringkas per kelompok rx
ovarian %>%
group_by(rx) %>%
summarise(
n = n(),
event = sum(fustat),
censored = n() - sum(fustat)
)## # A tibble: 2 × 4
## rx n event censored
## <dbl> <int> <dbl> <dbl>
## 1 1 13 7 6
## 2 2 13 5 89.3.2 Menyusun Tabel Event Time dan Risk Set
9.3.2.1 Menentukan Waktu Event
Uji log-rank hanya menggunakan waktu ketika event benar-benar terjadi, bukan waktu censoring.
event_times <- ovarian %>%
filter(fustat == 1) %>% # hanya yang mengalami event
arrange(futime) %>%
distinct(futime) %>% # ambil waktu event unik
pull(futime)
event_times## [1] 59 115 156 268 329 353 365 431 464 475 563 6389.3.2.2 Risk Set dan Jumlah Event per Waktu
make_logrank_table <- function(dat) {
times <- dat %>%
filter(fustat == 1) %>%
arrange(futime) %>%
distinct(futime) %>%
pull(futime)
out <- lapply(times, function(tj) {
# siapa saja yang masih at risk tepat sebelum tj?
at_risk <- dat %>%
filter(futime >= tj) # belum event / belum censor sebelum tj
# risk set per grup
n1j <- sum(at_risk$rx == 1)
n2j <- sum(at_risk$rx == 2)
# event pada tj per grup
d1j <- dat %>% filter(futime == tj, fustat == 1, rx == 1) %>% nrow()
d2j <- dat %>% filter(futime == tj, fustat == 1, rx == 2) %>% nrow()
data.frame(
time = tj,
Y1j = n1j,
Y2j = n2j,
d1j = d1j,
d2j = d2j
)
})
out <- bind_rows(out)
out <- out %>%
mutate(
Yj = Y1j + Y2j, # total at risk
dj = d1j + d2j # total event
)
out
}
logrank_tab <- make_logrank_table(ovarian)
kable(logrank_tab, caption = "Tabel Risk Set dan Event per Waktu") %>%
kable_styling(full_width = FALSE)| time | Y1j | Y2j | d1j | d2j | Yj | dj |
|---|---|---|---|---|---|---|
| 59 | 13 | 13 | 1 | 0 | 26 | 1 |
| 115 | 12 | 13 | 1 | 0 | 25 | 1 |
| 156 | 11 | 13 | 1 | 0 | 24 | 1 |
| 268 | 10 | 13 | 1 | 0 | 23 | 1 |
| 329 | 9 | 13 | 1 | 0 | 22 | 1 |
| 353 | 8 | 13 | 0 | 1 | 21 | 1 |
| 365 | 8 | 12 | 0 | 1 | 20 | 1 |
| 431 | 8 | 9 | 1 | 0 | 17 | 1 |
| 464 | 6 | 9 | 0 | 1 | 15 | 1 |
| 475 | 6 | 8 | 0 | 1 | 14 | 1 |
| 563 | 5 | 7 | 0 | 1 | 12 | 1 |
| 638 | 5 | 6 | 1 | 0 | 11 | 1 |
9.3.3 Perhitungan Manual Expected dan Variance
logrank_tab <- logrank_tab %>%
mutate(
# Expected event untuk grup 1 dan 2
E1j = (Y1j / Yj) * dj,
E2j = (Y2j / Yj) * dj,
# Komponen varians untuk grup 1
Vj = (Y1j * Y2j * dj * (Yj - dj)) / (Yj^2 * (Yj - 1))
)
# pastikan menggunakan dplyr::select
logrank_clean <- logrank_tab %>%
dplyr::select(time, Y1j, Y2j, d1j, d2j, Yj, dj, E1j, E2j, Vj)
kable(
logrank_clean,
digits = 4,
caption = "Perhitungan Manual: Observed, Expected, dan Komponen Varians"
) %>%
kable_styling(full_width = FALSE)| time | Y1j | Y2j | d1j | d2j | Yj | dj | E1j | E2j | Vj |
|---|---|---|---|---|---|---|---|---|---|
| 59 | 13 | 13 | 1 | 0 | 26 | 1 | 0.5000 | 0.5000 | 0.2500 |
| 115 | 12 | 13 | 1 | 0 | 25 | 1 | 0.4800 | 0.5200 | 0.2496 |
| 156 | 11 | 13 | 1 | 0 | 24 | 1 | 0.4583 | 0.5417 | 0.2483 |
| 268 | 10 | 13 | 1 | 0 | 23 | 1 | 0.4348 | 0.5652 | 0.2457 |
| 329 | 9 | 13 | 1 | 0 | 22 | 1 | 0.4091 | 0.5909 | 0.2417 |
| 353 | 8 | 13 | 0 | 1 | 21 | 1 | 0.3810 | 0.6190 | 0.2358 |
| 365 | 8 | 12 | 0 | 1 | 20 | 1 | 0.4000 | 0.6000 | 0.2400 |
| 431 | 8 | 9 | 1 | 0 | 17 | 1 | 0.4706 | 0.5294 | 0.2491 |
| 464 | 6 | 9 | 0 | 1 | 15 | 1 | 0.4000 | 0.6000 | 0.2400 |
| 475 | 6 | 8 | 0 | 1 | 14 | 1 | 0.4286 | 0.5714 | 0.2449 |
| 563 | 5 | 7 | 0 | 1 | 12 | 1 | 0.4167 | 0.5833 | 0.2431 |
| 638 | 5 | 6 | 1 | 0 | 11 | 1 | 0.4545 | 0.5455 | 0.2479 |
9.3.3.1 Akumulasi O, E, dan V serta Statistik Uji
O1 <- sum(logrank_tab$d1j) # total observed event di grup 1
E1 <- sum(logrank_tab$E1j) # total expected event di grup 1
V <- sum(logrank_tab$Vj) # total variance
O1; E1; V## [1] 7## [1] 5.233531## [1] 2.936196chisq_manual <- (O1 - E1)^2 / V
chisq_manual## [1] 1.06274p_value <- 1 - pchisq(chisq_manual, df = 1)
p_value## [1] 0.30259119.3.3.2 Membandingkan
dengan survdiff()
fit_logrank <- survdiff(Surv(futime, fustat) ~ rx, data = ovarian)
fit_logrank## Call:
## survdiff(formula = Surv(futime, fustat) ~ rx, data = ovarian)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=1 13 7 5.23 0.596 1.06
## rx=2 13 5 6.77 0.461 1.06
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3fit_logrank$chisq## [1] 1.06274Hasil fit_logrank$chisq harus sama (atau sangat dekat)
dengan chisq_manual yang kita hitung menggunakan rumus:
\[ \chi^2 = \frac{(O_1 - E_1)^2}{V}. \]
9.3.3.3 Tabel Gaya Excel: Observed, Not Observed, dan M.Sum
Untuk menyambungkan ke gambar Excel, kita bentuk tabel berikut.
excel_like <- logrank_tab %>%
transmute(
time,
Observed_rx1 = d1j,
NotObs_rx1 = Y1j - d1j,
MSum_rx1 = Y1j,
Observed_rx2 = d2j,
NotObs_rx2 = Y2j - d2j,
MSum_rx2 = Y2j,
Total_at_risk = Yj,
E_rx1 = E1j,
E_rx2 = E2j,
Vj = Vj,
Obs_minus_E_rx1 = d1j - E1j
)
kable(excel_like, digits = 4,
caption = "Tabel Gaya Excel: Observed, Not Observed, M.Sum, Expected, dan Variance") %>%
kable_styling(full_width = FALSE)| time | Observed_rx1 | NotObs_rx1 | MSum_rx1 | Observed_rx2 | NotObs_rx2 | MSum_rx2 | Total_at_risk | E_rx1 | E_rx2 | Vj | Obs_minus_E_rx1 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 59 | 1 | 12 | 13 | 0 | 13 | 13 | 26 | 0.5000 | 0.5000 | 0.2500 | 0.5000 |
| 115 | 1 | 11 | 12 | 0 | 13 | 13 | 25 | 0.4800 | 0.5200 | 0.2496 | 0.5200 |
| 156 | 1 | 10 | 11 | 0 | 13 | 13 | 24 | 0.4583 | 0.5417 | 0.2483 | 0.5417 |
| 268 | 1 | 9 | 10 | 0 | 13 | 13 | 23 | 0.4348 | 0.5652 | 0.2457 | 0.5652 |
| 329 | 1 | 8 | 9 | 0 | 13 | 13 | 22 | 0.4091 | 0.5909 | 0.2417 | 0.5909 |
| 353 | 0 | 8 | 8 | 1 | 12 | 13 | 21 | 0.3810 | 0.6190 | 0.2358 | -0.3810 |
| 365 | 0 | 8 | 8 | 1 | 11 | 12 | 20 | 0.4000 | 0.6000 | 0.2400 | -0.4000 |
| 431 | 1 | 7 | 8 | 0 | 9 | 9 | 17 | 0.4706 | 0.5294 | 0.2491 | 0.5294 |
| 464 | 0 | 6 | 6 | 1 | 8 | 9 | 15 | 0.4000 | 0.6000 | 0.2400 | -0.4000 |
| 475 | 0 | 6 | 6 | 1 | 7 | 8 | 14 | 0.4286 | 0.5714 | 0.2449 | -0.4286 |
| 563 | 0 | 5 | 5 | 1 | 6 | 7 | 12 | 0.4167 | 0.5833 | 0.2431 | -0.4167 |
| 638 | 1 | 4 | 5 | 0 | 6 | 6 | 11 | 0.4545 | 0.5455 | 0.2479 | 0.5455 |
Interpretasi satu baris pada waktu event \(t_j\):
Observed_rx1= \(d_{1j}\): jumlah event di grup 1.
NotObs_rx1= \(Y_{1j} - d_{1j}\): subjek grup 1 yang masih hidup / belum event pada \(t_j\).
MSum_rx1= \(Y_{1j}\): total at risk di grup 1.
Demikian pula untuk grup 2.
Bagian penting untuk dibawa ke slide/materi kuliah:
Kaplan–Meier:
\[ \widehat{S}(t) = \prod_{t_{(j)} \le t} \left(1 - \frac{d_j}{Y_j}\right) \]Log-rank (2 kelompok)
- Risk set dan event per waktu: \[ Y_{1j}, Y_{2j}, Y_j, d_{1j}, d_{2j}, d_j \]
- Expected: \[ E_{1j} = \frac{Y_{1j}}{Y_j} d_j \]
- Variance: \[ V_j = \frac{Y_{1j} Y_{2j} d_j (Y_j - d_j)}{Y_j^2 (Y_j - 1)} \]
- Akumulasi: \[ O_1 = \sum_j d_{1j}, \quad E_1 = \sum_j E_{1j}, \quad V = \sum_j V_j \]
- Statistik uji: \[ \chi^2 = \frac{(O_1 - E_1)^2}{V} \sim \chi^2_{(1)}. \]
Anda bisa memotong bagian-bagian rumus ini untuk dimasukkan ke PPT, e-book, atau handout Analisis Survival.
9.4 Estimasi Parametrik: Exponential & Weibull
Pendekatan parametrik mengasumsikan bentuk distribusi spesifik untuk \(T\). Keuntungannya: - Model lebih ringkas - Bisa mengekstrapolasi di luar rentang data (dengan hati-hati)
9.4.1 Model Exponential
Asumsi: hazard konstan \[ h(t) = \lambda, \quad \lambda > 0 \]
Maka: \[ S(t) = \exp(-\lambda t) \] \[ f(t) = \lambda \exp(-\lambda t) \]
Likelihood untuk data dengan censoring:
Misalkan: - \(t_i\) = waktu observasi -
\(\delta_i\) = indikator event (1) /
censored (0)
Untuk distribusi eksponensial: \[ L(\lambda) = \prod_{i=1}^n [f(t_i)]^{\delta_i} [S(t_i)]^{1-\delta_i} = \prod_{i=1}^n [\lambda e^{-\lambda t_i}]^{\delta_i} [e^{-\lambda t_i}]^{1-\delta_i} \]
\[ \ln L(\lambda) = \sum_{i=1}^n \delta_i \ln \lambda - \lambda \sum_{i=1}^n t_i \]
Turunan terhadap \(\lambda\) dan set ke nol: \[ \frac{d}{d\lambda} \ln L(\lambda) = \frac{\sum_{i=1}^n \delta_i}{\lambda} - \sum_{i=1}^n t_i = 0 \Rightarrow \hat{\lambda} = \frac{\sum_{i=1}^n \delta_i}{\sum_{i=1}^n t_i} \] Jadi estimasi MLE parameter eksponensial sangat sederhana.
9.4.2 Model Weibull
Model Weibull lebih fleksibel: - Hazard bisa meningkat atau menurun terhadap waktu.
Fungsi hazard: \[ h(t) = \lambda p t^{p-1}, \quad \lambda > 0,\, p > 0 \]
Fungsi survival: \[ S(t) = \exp(-\lambda t^p) \]
Interpretasi: - Jika \(p = 1\): sama dengan eksponensial - Jika \(p > 1\): hazard meningkat seiring waktu - Jika \(p < 1\): hazard menurun seiring waktu
Parameter biasanya diestimasi via maximum likelihood dengan optim numerik (software).
9.5 Model Semiparametrik: Cox Proportional Hazards
Model ini sangat populer karena: - Tidak perlu menspesifikkan bentuk baseline hazard \(h_0(t)\), - Tetap memungkinkan efek kovariat.
9.5.1 Spesifikasi Model Cox
Misalkan \(X_i\) adalah vektor kovariat untuk individu ke-\(i\). Model Cox ditulis sebagai:
\[ h(t \mid X_i) = h_0(t)\,\exp(\eta^\top X_i) \]
dengan:
- \(h_0(t)\): baseline hazard (tidak dispesifikkan bentuknya)
- \(\eta\): parameter regresi
- \(\exp(\eta_j)\): hazard ratio untuk kenaikan satu unit pada kovariat \(X_j\)
9.5.2 Partial Likelihood (Tanpa Bentuk Eksplisit \(h_0\))
Cox mengusulkan partial likelihood sebagai berikut.
Misalkan terdapat event pada waktu \(t_{(1)} < t_{(2)} < \dots < t_{(k)}\), dan individu yang mengalami event ke-\(j\) adalah \(i_j\).
Set risk set \(R_j\) didefinisikan sebagai himpunan individu yang masih berisiko tepat sebelum waktu \(t_{(j)}\).
Partial likelihood ditulis sebagai:
\[ L_p(\eta) = \prod_{j=1}^k \frac{\exp(\eta^\top X_{i_j})}{\sum_{l \in R_j} \exp(\eta^\top X_l)} \]
Log partial likelihood:
\[ \ell_p(\eta) = \sum_{j=1}^k \left[ \eta^\top X_{i_j} - \log \left( \sum_{l \in R_j} \exp(\eta^\top X_l) \right) \right] \]
Estimasi \(\hat{\eta}\) diperoleh dengan memaksimumkan \(\ell_p(\eta)\).
9.5.3 Interpretasi Koefisien Cox
Jika model dituliskan sebagai:
\[ h(t \mid X) = h_0(t)\,\exp(\eta_1 X_1 + \dots + \eta_p X_p) \]
maka hazard ratio untuk kovariat \(X_j\) adalah:
\[ \text{HR}_j = \exp(\eta_j) \]
Interpretasi:
- \(\text{HR}_j > 1\): peningkatan \(X_j\) meningkatkan hazard (risiko event lebih besar).
- \(\text{HR}_j < 1\): peningkatan \(X_j\) menurunkan hazard (efek protektif).
- \(\text{HR}_j = 1.5\): hazard meningkat 50% untuk setiap satu unit kenaikan \(X_j\).
9.6 Contoh Analisis Survival dengan R
Kita gunakan data lung dari paket survival.
Ini adalah data pasien kanker paru.
data("lung", package = "survival")
?lung## Help on topic 'lung' was found in the following packages:
##
## Package Library
## KMsurv /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## survival /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
##
##
## Using the first match ...head(lung)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0summary(lung)## inst time status age
## Min. : 1.00 Min. : 5.0 Min. :1.000 Min. :39.00
## 1st Qu.: 3.00 1st Qu.: 166.8 1st Qu.:1.000 1st Qu.:56.00
## Median :11.00 Median : 255.5 Median :2.000 Median :63.00
## Mean :11.09 Mean : 305.2 Mean :1.724 Mean :62.45
## 3rd Qu.:16.00 3rd Qu.: 396.5 3rd Qu.:2.000 3rd Qu.:69.00
## Max. :33.00 Max. :1022.0 Max. :2.000 Max. :82.00
## NA's :1
## sex ph.ecog ph.karno pat.karno
## Min. :1.000 Min. :0.0000 Min. : 50.00 Min. : 30.00
## 1st Qu.:1.000 1st Qu.:0.0000 1st Qu.: 75.00 1st Qu.: 70.00
## Median :1.000 Median :1.0000 Median : 80.00 Median : 80.00
## Mean :1.395 Mean :0.9515 Mean : 81.94 Mean : 79.96
## 3rd Qu.:2.000 3rd Qu.:1.0000 3rd Qu.: 90.00 3rd Qu.: 90.00
## Max. :2.000 Max. :3.0000 Max. :100.00 Max. :100.00
## NA's :1 NA's :1 NA's :3
## meal.cal wt.loss
## Min. : 96.0 Min. :-24.000
## 1st Qu.: 635.0 1st Qu.: 0.000
## Median : 975.0 Median : 7.000
## Mean : 928.8 Mean : 9.832
## 3rd Qu.:1150.0 3rd Qu.: 15.750
## Max. :2600.0 Max. : 68.000
## NA's :47 NA's :14Penjelasan variabel (singkat): - time: waktu follow-up
(hari) - status: 1 = censored, 2 = death (perlu di-recode)
- age: usia - sex: 1 = male, 2 = female -
ph.ecog: ECOG performance score (semakin tinggi = kondisi
lebih buruk), dsb.
9.6.1 Menyiapkan Objek Surv
Dalam paket survival, kita pakai fungsi
Surv().
lung2 <- lung %>%
mutate(
status = ifelse(status == 2, 1, 0) # 1 = event (death), 0 = censored
)
surv_obj <- with(lung2, Surv(time = time, event = status))
head(surv_obj)## [1] 306 455 1010+ 210 883 1022+Objek Surv menyimpan informasi waktu dan status
censoring.
9.7 Estimasi Kaplan–Meier dan Plot Kurva Survival
9.7.1 Kaplan–Meier Tanpa Kovariat (Overall)
fit_km <- survfit(surv_obj ~ 1, data = lung2)
fit_km## Call: survfit(formula = surv_obj ~ 1, data = lung2)
##
## n events median 0.95LCL 0.95UCL
## [1,] 228 165 310 285 363Menampilkan ringkasan survival:
summary(fit_km)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## 228.00000 228.00000 228.00000 165.00000 376.27475 19.70779 310.00000 285.00000
## 0.95UCL
## 363.00000Plot kurva survival:
ggsurvplot(
fit_km,
data = lung2,
risk.table = TRUE,
conf.int = TRUE,
xlab = "Waktu (hari)",
ylab = "Probabilitas Survival",
title = "Kurva Kaplan–Meier (Overall)",
surv.median.line = "hv"
)
Interpretasi: - Kurva survival menunjukkan probabilitas bertahan hidup seiring waktu. - Garis horizontal/vertikal median menunjukkan waktu ketika survival = 0.5 (median survival). - Risk table menunjukkan berapa banyak pasien yang masih diikuti (at risk) pada tiap titik waktu.
9.7.2 Kaplan–Meier Dibedakan menurut Jenis Kelamin
lung2 <- lung2 %>%
mutate(
sex_factor = factor(sex, labels = c("Laki-laki", "Perempuan"))
)
fit_km_sex <- survfit(surv_obj ~ sex_factor, data = lung2)
ggsurvplot(
fit_km_sex,
data = lung2,
risk.table = TRUE,
conf.int = TRUE,
pval = TRUE,
xlab = "Waktu (hari)",
ylab = "Probabilitas Survival",
title = "Kurva Kaplan–Meier menurut Jenis Kelamin",
legend.title = "Jenis Kelamin",
legend.labs = c("Laki-laki", "Perempuan")
)
Interpretasi: - Kita bandingkan kurva survival
laki-laki vs perempuan. - Nilai pval (log-rank test)
menunjukkan apakah perbedaan survival antara grup signifikan secara
statistik. - Misal: jika p-value < 0.05, ada bukti bahwa survival
laki-laki dan perempuan berbeda.
9.8 Model Cox Proportional Hazards
Mari kita bangun model Cox dengan kovariat: - sex_factor
(laki vs perempuan) - age (usia) - ph.ecog
(kondisi klinis)
9.8.1 Fitting Model Cox
cox_fit <- coxph(
Surv(time, status) ~ sex_factor + age + ph.ecog,
data = lung2
)
summary(cox_fit)## Call:
## coxph(formula = Surv(time, status) ~ sex_factor + age + ph.ecog,
## data = lung2)
##
## n= 227, number of events= 164
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex_factorPerempuan -0.552612 0.575445 0.167739 -3.294 0.000986 ***
## age 0.011067 1.011128 0.009267 1.194 0.232416
## ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex_factorPerempuan 0.5754 1.7378 0.4142 0.7994
## age 1.0111 0.9890 0.9929 1.0297
## ph.ecog 1.5900 0.6289 1.2727 1.9864
##
## Concordance= 0.637 (se = 0.025 )
## Likelihood ratio test= 30.5 on 3 df, p=1e-06
## Wald test = 29.93 on 3 df, p=1e-06
## Score (logrank) test = 30.5 on 3 df, p=1e-06Output utama: - coef = \(\hat{eta}\) - exp(coef) =
hazard ratio - se(coef), z,
Pr(>|z|) = uji signifikansi - Concordance,
Likelihood ratio test, dsb.
9.8.2 Interpretasi Koefisien
Misalkan (contoh pola interpretasi; angka aktual lihat output):
sex_factorPerempuandengan HR < 1:
Berarti perempuan memiliki hazard (risiko kematian) yang lebih rendah dibanding laki-laki.
Misal HR = 0.7 → hazard 30% lebih rendah.agedengan HR sedikit > 1:
Misal HR = 1.02 → setiap kenaikan 1 tahun usia, hazard meningkat 2%.ph.ecogdengan HR > 1, misal HR = 1.5: Setiap peningkatan satu poin ECOG (kondisi lebih buruk), hazard meningkat 50%.
9.8.3 Tabel Hazard Ratio yang Lebih Rapih
# Jika paket broom belum terpasang, jalankan dulu:
# install.packages("broom")
broom::tidy(cox_fit, exponentiate = TRUE, conf.int = TRUE)9.9 Plot Survival Adjusted (Berbasis Model Cox)
Kita bisa menggambar kurva survival yang diprediksi dari model Cox untuk berbagai kombinasi kovariat.
Contoh: prediksi survival untuk laki-laki vs perempuan pada usia dan
ph.ecog tertentu.
# Data baru untuk prediksi
newdata <- data.frame(
sex_factor = factor(c("Laki-laki", "Perempuan"),
levels = levels(lung2$sex_factor)),
age = c(60, 60),
ph.ecog = c(1, 1)
)
cox_fit## Call:
## coxph(formula = Surv(time, status) ~ sex_factor + age + ph.ecog,
## data = lung2)
##
## coef exp(coef) se(coef) z p
## sex_factorPerempuan -0.552612 0.575445 0.167739 -3.294 0.000986
## age 0.011067 1.011128 0.009267 1.194 0.232416
## ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05
##
## Likelihood ratio test=30.5 on 3 df, p=1.083e-06
## n= 227, number of events= 164
## (1 observation deleted due to missingness)fit_surv_cox <- survfit(cox_fit, newdata = newdata)
ggsurvplot(
fit_surv_cox,
data = newdata,
conf.int = TRUE,
legend.labs = c("Laki-laki, 60 th, ECOG 1", "Perempuan, 60 th, ECOG 1"),
legend.title = "Profil Pasien",
xlab = "Waktu (hari)",
ylab = "Probabilitas Survival",
title = "Kurva Survival Diprediksi dari Model Cox"
)
Interpretasi: - Kurva memperlihatkan survival yang diharapkan untuk profil pasien tertentu. - Selisih antara kurva mencerminkan perbedaan hazard ratio antar kovariat.
9.10 Mengecek Asumsi Proportional Hazards
Asumsi utama model Cox: - Hazard ratio konstan terhadap waktu (proportional hazards).
Kita dapat mengecek dengan uji Schoenfeld residuals:
ph_test <- cox.zph(cox_fit)
ph_test## chisq df p
## sex_factor 2.305 1 0.13
## age 0.188 1 0.66
## ph.ecog 2.054 1 0.15
## GLOBAL 4.464 3 0.22plot(ph_test)
- Jika p-value kecil (misal < 0.05), indikasi pelanggaran asumsi PH untuk kovariat tersebut.
- Plot residual terhadap waktu: jika garis mendekati horizontal → asumsi PH cukup baik.
9.11 Model Parametrik
dengan survreg
Sebagai tambahan, kita bisa memodelkan survival secara parametrik
dengan survreg().
Contoh: model Weibull.
weibull_fit <- survreg(
Surv(time, status) ~ sex_factor + age + ph.ecog,
data = lung2,
dist = "weibull"
)
summary(weibull_fit)##
## Call:
## survreg(formula = Surv(time, status) ~ sex_factor + age + ph.ecog,
## data = lung2, dist = "weibull")
## Value Std. Error z p
## (Intercept) 6.67453 0.42740 15.62 < 2e-16
## sex_factorPerempuan 0.40109 0.12373 3.24 0.0012
## age -0.00748 0.00676 -1.11 0.2690
## ph.ecog -0.33964 0.08348 -4.07 4.7e-05
## Log(scale) -0.31319 0.06135 -5.11 3.3e-07
##
## Scale= 0.731
##
## Weibull distribution
## Loglik(model)= -1132.4 Loglik(intercept only)= -1147.4
## Chisq= 29.98 on 3 degrees of freedom, p= 1.4e-06
## Number of Newton-Raphson Iterations: 5
## n=227 (1 observation deleted due to missingness)Catatan: - Parametrisasi survreg menggunakan bentuk
accelerated failure time (AFT). - Interpretasi koefisien: efek
pada log-waktu survival, bukan langsung hazard ratio. - Namun kita bisa
transform untuk mendapatkan parameter hazard/dekadence tertentu.
9.12 Ringkasan & Intuisi Konsep
9.12.1 Intuisi Survival Function
- \(S(t)\) tinggi → sebagian besar subjek masih bertahan pada waktu \(t\).
- Turun tajam di awal → banyak event terjadi pada awal follow-up.
9.12.2 Intuisi Hazard
- \(h(t)\) = “laju risiko sesaat”:
- Jika hazard besar di awal lalu turun: risiko tertinggi terjadi di awal (contoh: operasi besar).
- Jika hazard kecil di awal lalu naik: risiko meningkat seiring waktu (contoh: keausan mesin).
9.12.3 Intuisi Kaplan–Meier
- Kaplan–Meier memperkirakan survival secara empiris, dengan tetap memperhitungkan censoring.
- Setiap kali terjadi event, survival dikalikan dengan faktor:
\[ \left(1 - \frac{d_j}{n_j}\right) \]
9.12.4 Intuisi Cox PH
- Cox memodelkan perbedaan hazard antar individu melalui kovariat, dengan asumsi rasio hazard konstan terhadap waktu.
- Baseline hazard \(h_0(t)\) boleh bentuk apa saja → sangat fleksibel.
- Koefisien \(eta\) memberi hazard ratio yang mudah ditafsirkan.
9.13 Lampiran: Contoh Simulasi Data Survival Sederhana
Sebagai ilustrasi tambahan, berikut contoh simulasi data survival eksponensial dengan satu kovariat biner.
set.seed(123)
n <- 300
x <- rbinom(n, 1, 0.5) # 0 atau 1
lambda0 <- 0.01 # baseline hazard
beta <- log(2) # hazard ratio = 2 jika x = 1
# Hazard individu: lambda_i = lambda0 * exp(beta * x)
lambda_i <- lambda0 * exp(beta * x)
# Generate waktu event dari exponential: T ~ Exp(lambda_i)
# Jika U ~ Uniform(0,1), maka T = -log(U)/lambda_i
u <- runif(n)
T_true <- -log(u) / lambda_i
# Tambahkan censoring acak
C <- rexp(n, rate = 0.005) # censoring times
time_obs <- pmin(T_true, C)
status_obs <- as.integer(T_true <= C)
sim_data <- data.frame(
time = time_obs,
status = status_obs,
x = factor(x, labels = c("Grup 0", "Grup 1"))
)
head(sim_data)## time status x
## 1 24.26128 1 Grup 0
## 2 193.42511 0 Grup 1
## 3 24.96597 1 Grup 0
## 4 15.77729 1 Grup 1
## 5 23.09131 1 Grup 1
## 6 73.20734 1 Grup 0table(sim_data$status)##
## 0 1
## 74 226Fitting Kaplan–Meier dan Cox:
# KM by group
fit_km_sim <- survfit(Surv(time, status) ~ x, data = sim_data)
ggsurvplot(
fit_km_sim,
data = sim_data,
risk.table = TRUE,
pval = TRUE,
title = "Kurva Kaplan–Meier (Data Simulasi)",
xlab = "Waktu",
ylab = "Probabilitas Survival"
)
# Model Cox
cox_sim <- coxph(Surv(time, status) ~ x, data = sim_data)
summary(cox_sim)## Call:
## coxph(formula = Surv(time, status) ~ x, data = sim_data)
##
## n= 300, number of events= 226
##
## coef exp(coef) se(coef) z Pr(>|z|)
## xGrup 1 0.4405 1.5535 0.1341 3.286 0.00102 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## xGrup 1 1.554 0.6437 1.195 2.02
##
## Concordance= 0.565 (se = 0.019 )
## Likelihood ratio test= 10.73 on 1 df, p=0.001
## Wald test = 10.8 on 1 df, p=0.001
## Score (logrank) test = 10.97 on 1 df, p=9e-04Perhatikan: - Hazard ratio hasil model Cox seharusnya mendekati nilai “benar” yaitu 2 (karena kita set \(eta = \log 2\)).
9.14 Penutup
Dalam R Markdown ini kita telah membahas:
- Konsep dasar analisis survival:
- Fungsi survival, hazard, dan hazard kumulatif.
- Data tersensor dan representasinya.
- Metode nonparametrik:
- Estimasi Kaplan–Meier.
- Estimasi hazard kumulatif Nelson–Aalen.
- Metode parametrik:
- Model Eksponensial dan Weibull.
- Estimasi MLE.
- Model semiparametrik Cox PH:
- Partial likelihood.
- Interpretasi hazard ratio.
- Pengecekan asumsi proportional hazards.
- Contoh lengkap dengan data
lungdan simulasi data eksponensial.
Silakan modifikasi, tambahkan penjelasan teoretis yang lebih dalam (misalnya derivasi MLE lengkap, link ke teori likelihood, atau perluasan ke model competing risks) sesuai kebutuhan materi kuliah atau modul pelatihan.
9.15 Appendix Cox Proporsional Hazard
Pada analisis Kaplan–Meier dan uji log-rank, kita hanya
membandingkan kurva survival antar kelompok (misal obat
A vs B) satu per satu.
Jika ingin memasukkan banyak kovariat sekaligus (usia,
jenis kelamin, status merokok, jenis obat, dll.), kita membutuhkan
model regresi untuk waktu survival, yaitu:
Cox Proportional Hazards Regression (Cox PH).
Slide kuliah menggambarkan bahwa:
Kita ingin memodelkan hazard (laju kejadian sesaat) sebagai fungsi dari beberapa faktor risiko, sehingga kita dapat memahami bagaimana kovariat memengaruhi risiko kejadian pada setiap titik waktu.
- Output penting: koefisien regresi dan Hazard Ratio (HR).
- Contoh bentuk model di slide: \[ \lambda(t) = \lambda_0(t)\,\exp(-5\,\text{gender} + 2\,\text{smokingStatus} + 0.01\,\text{Age}) \] di mana \(\lambda_0(t)\) adalah baseline hazard.
Dokumen ini menyusun:
- Definisi fungsi hazard dan survival.
- Bentuk umum model Cox PH.
- Formulasi matematis lengkap: hazard ratio, log-hazard, partial likelihood.
- Contoh perhitungan manual step-by-step dengan data
ovarianuntuk model sederhana. - Perbandingan hasil manual dengan fungsi
coxph()di R.
9.15.1 Dasar Teori Survival dan Hazard
9.15.1.1 Fungsi distribusi, survival, dan hazard
Misalkan \(T\) adalah peubah acak waktu sampai kejadian (event), dengan \(T \ge 0\).
- Fungsi distribusi kumulatif (CDF): \[ F(t) = P(T \le t) \]
- Fungsi survival: \[ S(t) = P(T > t) = 1 - F(t) \]
- Fungsi densitas (jika kontinu): \[ f(t) = \frac{d}{dt} F(t) \]
Fungsi hazard (instantaneous failure rate) didefinisikan sebagai:
\[ h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} \]
Jika \(T\) kontinu, hubungan penting:
\[ h(t) = \frac{f(t)}{S(t)} \]
dan
\[ S(t) = \exp\left(-\int_0^t h(u)\,du\right) \]
9.15.2 Model Cox Proportional Hazards
9.15.2.1 Bentuk umum model
Misalkan kita punya vektor kovariat \(\mathbf{x} = (x_1,\dots,x_p)'\).
Model Cox PH menyatakan:
\[ h(t \mid \mathbf{x}) = h_0(t)\,\exp(\beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p) \]
dengan:
- \(h_0(t)\) = baseline hazard (hazard untuk subjek dengan semua kovariat 0).
- \(\beta_1,\dots,\beta_p\) = koefisien regresi yang mengukur efek kovariat pada hazard.
Model ini semi-parametrik:
- Bagian \(\exp(\beta'x)\) bersifat parametrik.
- Bentuk \(h_0(t)\) tidak perlu dispesifikasikan (nonparametrik).
9.15.2.2 Log-hazard dan interpretasi koefisien
Ambil log di kedua sisi:
\[ \log h(t \mid \mathbf{x}) = \log h_0(t) + \beta_1 x_1 + \cdots + \beta_p x_p \]
Untuk satu kovariat \(x_k\):
- Jika \(x_k\) naik 1 unit (kovariat lain tetap), perubahan log-hazard = \(\beta_k\).
- Hazard ratio (HR) untuk kenaikan 1 unit:
\[ \text{HR}_k = \exp(\beta_k) \]
Jika \(\beta_k > 0\): HR > 1 →
hazard meningkat → survival memburuk.
Jika \(\beta_k < 0\): HR < 1 →
hazard menurun → survival membaik.
Slide memberikan contoh:
\[ \log \lambda(t) = \log \lambda_0(t) -5\,\text{gender} + 2\,\text{smokingStatus} + 0.01\,\text{Age} \]
- Koefisien 2 untuk
smokingStatus→ HR ≈ \(\exp(2) \approx 7.39\)
Artinya perokok memiliki hazard sekitar 7.39 kali hazard non-perokok (ceteris paribus). fileciteturn0file0
9.15.2.3 Prinsip proportional hazards
Untuk dua subjek dengan kovariat \(\mathbf{x}_A\) dan \(\mathbf{x}_B\), hazard ratio:
\[ \frac{h(t \mid \mathbf{x}_A)}{h(t \mid \mathbf{x}_B)} = \frac{h_0(t)\exp(\beta'\mathbf{x}_A)}{h_0(t)\exp(\beta'\mathbf{x}_B)} = \exp\big(\beta'(\mathbf{x}_A - \mathbf{x}_B)\big) \]
Catatan penting:
- Tidak tergantung waktu \(t\) → inilah asumsi proportional hazards.
- HR konstan sepanjang waktu.
9.15.3 Estimasi Parameter: Partial Likelihood
9.15.3.1 Ide dasar
Karena \(h_0(t)\) tidak dispesifikasikan, Cox mengusulkan partial likelihood yang:
- Menggunakan informasi urutan dan posisi waktu event.
- Tidak perlu memodelkan \(h_0(t)\).
Anggap:
- Ada \(n\) subjek.
- Terdapat \(K\) waktu event terurut: \(t_{(1)} < t_{(2)} < \dots < t_{(K)}\).
- Pada setiap waktu event \(t_{(k)}\), tepat satu subjek mengalami event (untuk memudahkan penjelasan; kasus ties bisa ditangani dengan koreksi Efron/Breslow).
Definisikan:
- \(R_k\) = risk set pada waktu \(t_{(k)}\): semua subjek yang masih at risk tepat sebelum \(t_{(k)}\).
- \(i(k)\) = indeks subjek yang mengalami event pada \(t_{(k)}\).
- \(\mathbf{x}_{i(k)}\) = kovariat subjek yang event pada waktu itu.
Probabilitas bahwa subjek \(i(k)\) yang mengalami event di antara semua anggota risk set \(R_k\) (dengan asumsi model Cox benar):
\[ P(\text{subjek } i(k) \text{ yang gagal pada } t_{(k)} \mid \text{seseorang gagal di } R_k) = \frac{h(t_{(k)} \mid \mathbf{x}_{i(k)})}{\sum_{l \in R_k} h(t_{(k)} \mid \mathbf{x}_l)} \]
Karena \(h(t \mid \mathbf{x}) = h_0(t)\exp(\beta'\mathbf{x})\), maka:
\[ \frac{h_0(t_{(k)})\exp(\beta'\mathbf{x}_{i(k)})} {\sum_{l \in R_k} h_0(t_{(k)})\exp(\beta'\mathbf{x}_l)} = \frac{\exp(\beta'\mathbf{x}_{i(k)})} {\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \]
Baseline hazard \(h_0(t)\) hilang dari rasio → inilah trik partial likelihood.
9.15.3.2 Partial likelihood
Dengan mengalikan probabilitas untuk semua waktu event:
\[ L_p(\beta) = \prod_{k=1}^K \frac{\exp(\beta'\mathbf{x}_{i(k)})} {\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \]
Log-partial likelihood:
\[ \ell_p(\beta) = \sum_{k=1}^K \left[ \beta'\mathbf{x}_{i(k)} - \log\left( \sum_{l \in R_k} \exp(\beta'\mathbf{x}_l) \right) \right] \]
Estimasi \(\hat{\beta}\) diperoleh dengan memaksimalkan \(\ell_p(\beta)\).
9.15.3.3 Score function dan matriks informasi
Turunan pertama (score vector):
\[ U(\beta) = \frac{\partial \ell_p(\beta)}{\partial \beta} = \sum_{k=1}^K \left[ \mathbf{x}_{i(k)} - \frac{\sum_{l \in R_k} \mathbf{x}_l \exp(\beta'\mathbf{x}_l)}{\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \right] \]
Istilah kedua adalah rata-rata kovariat berbobot pada risk set.
Turunan kedua (observed information matrix):
\[ I(\beta) = -\frac{\partial^2 \ell_p(\beta)}{\partial \beta \partial \beta'} = \sum_{k=1}^K \left[ \frac{\sum_{l \in R_k} \mathbf{x}_l \mathbf{x}_l' \exp(\beta'\mathbf{x}_l)}{\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} - \left( \frac{\sum_{l \in R_k} \mathbf{x}_l \exp(\beta'\mathbf{x}_l)}{\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \right) \left( \frac{\sum_{l \in R_k} \mathbf{x}_l \exp(\beta'\mathbf{x}_l)}{\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \right)' \right] \]
Estimator Cox PH:
\[ \hat{\beta} = \arg\max_\beta \ell_p(\beta) \]
dan varian asimtotik:
\[ \widehat{\text{Var}}(\hat{\beta}) = I(\hat{\beta})^{-1} \]
9.15.4 Contoh Data:
ovarian
Slide kuliah menggunakan data ovarian untuk ilustrasi
Kaplan–Meier dan log-rank. fileciteturn0file0
Kita juga akan memakainya untuk Cox PH (satu kovariat:
rx).
data(ovarian)
ovarian <- ovarian
head(ovarian)## futime fustat age resid.ds rx ecog.ps
## 1 59 1 72.3315 2 1 1
## 2 115 1 74.4932 2 1 1
## 3 156 1 66.4658 2 1 2
## 4 421 0 53.3644 2 2 1
## 5 431 1 50.3397 2 1 1
## 6 448 0 56.4301 1 1 2summary(ovarian$rx)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 1.0 1.5 1.5 2.0 2.0# Ringkas jumlah event dan censoring per grup
ovarian %>%
group_by(rx) %>%
summarise(
n = n(),
event = sum(fustat),
censored = n() - sum(fustat)
)## # A tibble: 2 × 4
## rx n event censored
## <dbl> <int> <dbl> <dbl>
## 1 1 13 7 6
## 2 2 13 5 8Kita pasang model:
\[ h(t \mid \text{rx}) = h_0(t)\,\exp(\beta\,\text{rx}) \]
dengan:
rx = 1→ kelompok referensi (misalnya obat A)rx = 2→ obat B
9.15.5 Estimasi Cox
dengan coxph()
fit_cox <- coxph(Surv(futime, fustat) ~ rx, data = ovarian)
summary(fit_cox)## Call:
## coxph(formula = Surv(futime, fustat) ~ rx, data = ovarian)
##
## n= 26, number of events= 12
##
## coef exp(coef) se(coef) z Pr(>|z|)
## rx -0.5964 0.5508 0.5870 -1.016 0.31
##
## exp(coef) exp(-coef) lower .95 upper .95
## rx 0.5508 1.816 0.1743 1.74
##
## Concordance= 0.608 (se = 0.07 )
## Likelihood ratio test= 1.05 on 1 df, p=0.3
## Wald test = 1.03 on 1 df, p=0.3
## Score (logrank) test = 1.06 on 1 df, p=0.3Output penting:
- Estimate \(\hat{\beta}\) untuk
rx. exp(coef)= hazard ratio (HR).se(coef)= standar error.zdanPr(>|z|)= uji Wald.- CI 95% untuk HR.
Kita simpan angka-angka tersebut:
beta_hat <- coef(fit_cox)["rx"]
se_beta <- sqrt(vcov(fit_cox)["rx", "rx"])
beta_hat## rx
## -0.59638se_beta## [1] 0.5869895HR <- exp(beta_hat)
CI_HR <- exp(beta_hat + c(-1, 1) * 1.96 * se_beta)
HR## rx
## 0.5508019CI_HR## [1] 0.174317 1.7404089.15.5.1 Perhitungan Manual Partial Likelihood (Model 1 Kovariat)
Di sini kita menunjukkan langkah manual sesuai formulasi partial likelihood.
9.15.5.2 Menyiapkan data terurut
library(survival)
library(dplyr)
data(ovarian)
# Buat data sederhana berisi waktu, status, dan rx
dat <- ovarian %>%
dplyr::select(futime, fustat, rx) %>%
arrange(futime)
head(dat, 10)## futime fustat rx
## 1 59 1 1
## 2 115 1 1
## 3 156 1 1
## 22 268 1 1
## 23 329 1 1
## 24 353 1 2
## 25 365 1 2
## 26 377 0 2
## 4 421 0 2
## 5 431 1 1Kita hanya menggunakan baris dengan event untuk partial likelihood; baris sensored hanya mempengaruhi risk set.
event_idx <- which(dat$fustat == 1)
event_times <- dat$futime[event_idx]
event_times## [1] 59 115 156 268 329 353 365 431 464 475 563 638length(event_times) # jumlah event K## [1] 129.15.5.3 Membentuk risk set dan kontribusi partial likelihood
Untuk setiap event ke-\(k\) pada waktu \(t_k\):
- Tentukan risk set \(R_k = \{ i : T_i \ge t_k \}\).
- Catat kovariat
rxuntuk subjek yang event (numerator). - Hitung jumlah exp(\(\beta \cdot \text{rx}_i\)) di risk set (denominator).
Untuk model 1 kovariat, \(\beta\) skalar.
# Fungsi untuk membangun tabel partial likelihood pada suatu nilai beta
build_partial_table <- function(beta, dat){
# dat harus sudah terurut naik berdasarkan futime
out_list <- list()
ev_rows <- which(dat$fustat == 1)
for(k in seq_along(ev_rows)){
i_event <- ev_rows[k]
t_k <- dat$futime[i_event]
# risk set: semua dengan futime >= t_k
R_k <- which(dat$futime >= t_k)
# kovariat untuk event dan untuk risk set
x_event <- dat$rx[i_event]
x_R <- dat$rx[R_k]
num <- exp(beta * x_event)
den <- sum(exp(beta * x_R))
out_list[[k]] <- data.frame(
k = k,
t_k = t_k,
i_event = i_event,
x_event = x_event,
R_size = length(R_k),
num = num,
den = den,
term = num / den # kontribusi L_p di waktu ke-k
)
}
do.call(rbind, out_list)
}
# Coba dengan beta = 0 (tidak ada efek rx)
partial_beta0 <- build_partial_table(beta = 0, dat = dat)
kable(head(partial_beta0, 10), digits = 4,
caption = "Kontribusi Partial Likelihood untuk \\(\\beta = 0\\) (10 event pertama)") %>%
kable_styling(full_width = FALSE)| k | t_k | i_event | x_event | R_size | num | den | term |
|---|---|---|---|---|---|---|---|
| 1 | 59 | 1 | 1 | 26 | 1 | 26 | 0.0385 |
| 2 | 115 | 2 | 1 | 25 | 1 | 25 | 0.0400 |
| 3 | 156 | 3 | 1 | 24 | 1 | 24 | 0.0417 |
| 4 | 268 | 4 | 1 | 23 | 1 | 23 | 0.0435 |
| 5 | 329 | 5 | 1 | 22 | 1 | 22 | 0.0455 |
| 6 | 353 | 6 | 2 | 21 | 1 | 21 | 0.0476 |
| 7 | 365 | 7 | 2 | 20 | 1 | 20 | 0.0500 |
| 8 | 431 | 10 | 1 | 17 | 1 | 17 | 0.0588 |
| 9 | 464 | 12 | 2 | 15 | 1 | 15 | 0.0667 |
| 10 | 475 | 13 | 2 | 14 | 1 | 14 | 0.0714 |
Untuk \(\beta = 0\):
num = 1, den = ukuran risk set.
Partial likelihood \(L_p(\beta)\) dan log-partial likelihood:
Lp_beta0 <- prod(partial_beta0$term)
logLp_beta0 <- sum(log(partial_beta0$term))
Lp_beta0## [1] 6.400788e-16logLp_beta0## [1] -34.984949.15.5.4 Partial likelihood pada \(\hat{\beta}\)
Sekarang kita gunakan \(\hat{\beta}\) dari coxph()
untuk melihat berapa nilai log-partial likelihood.
partial_betahat <- build_partial_table(beta = beta_hat, dat = dat)
kable(head(partial_betahat, 10), digits = 4,
caption = "Kontribusi Partial Likelihood untuk \\(\\beta = \\hat\\beta\\) (10 event pertama)") %>%
kable_styling(full_width = FALSE)| k | t_k | i_event | x_event | R_size | num | den | term | |
|---|---|---|---|---|---|---|---|---|
| rx | 1 | 59 | 1 | 1 | 26 | 0.5508 | 11.1044 | 0.0496 |
| rx1 | 2 | 115 | 2 | 1 | 25 | 0.5508 | 10.5536 | 0.0522 |
| rx2 | 3 | 156 | 3 | 1 | 24 | 0.5508 | 10.0028 | 0.0551 |
| rx3 | 4 | 268 | 4 | 1 | 23 | 0.5508 | 9.4520 | 0.0583 |
| rx4 | 5 | 329 | 5 | 1 | 22 | 0.5508 | 8.9012 | 0.0619 |
| rx5 | 6 | 353 | 6 | 2 | 21 | 0.3034 | 8.3504 | 0.0363 |
| rx6 | 7 | 365 | 7 | 2 | 20 | 0.3034 | 8.0470 | 0.0377 |
| rx7 | 8 | 431 | 10 | 1 | 17 | 0.5508 | 7.1369 | 0.0772 |
| rx8 | 9 | 464 | 12 | 2 | 15 | 0.3034 | 6.0353 | 0.0503 |
| rx9 | 10 | 475 | 13 | 2 | 14 | 0.3034 | 5.7319 | 0.0529 |
Lp_betahat <- prod(partial_betahat$term)
logLp_betahat <- sum(log(partial_betahat$term))
Lp_betahat## [1] 1.082813e-15logLp_betahat## [1] -34.45921Bandingkan log-partial likelihood:
c(
logL_beta0 = logLp_beta0,
logL_betahat = logLp_betahat,
logL_coxph = fit_cox$loglik[2] # from coxph: loglik at beta_hat
)## logL_beta0 logL_betahat logL_coxph
## -34.98494 -34.45921 -34.45921Perbedaan kecil bisa muncul karena rounding atau penanganan ties (Cox PH default menggunakan metode Efron/Breslow jika ada ties).
9.15.6 Uji Signifikansi dan Interval Kepercayaan
9.15.6.1 Uji Wald
Untuk parameter tunggal \(\beta\):
- H\(_0\): \(\beta = 0\)
- Statistik uji: \[ Z_{\text{Wald}} = \frac{\hat{\beta}}{\text{SE}(\hat{\beta})} \sim N(0, 1) \quad \text{(asimtotik)} \]
p-value:
\[ p = 2 \left[ 1 - \Phi\left( |Z_{\text{Wald}}| \right) \right] \]
Z_wald <- beta_hat / se_beta
p_wald <- 2 * (1 - pnorm(abs(Z_wald)))
Z_wald## rx
## -1.015998p_wald## rx
## 0.309630495% Confidence Interval (CI) untuk \(\beta\):
\[ \hat{\beta} \pm z_{0.975} \,\text{SE}(\hat{\beta}) \]
CI untuk HR:
\[ \text{CI}_{\text{HR}} = \left[ \exp(\hat{\beta} - 1.96\,\text{SE}),\ \exp(\hat{\beta} + 1.96\,\text{SE}) \right] \]
CI_beta <- beta_hat + c(-1, 1) * 1.96 * se_beta
CI_beta## [1] -1.7468795 0.5541194exp(CI_beta)## [1] 0.174317 1.740408Hasilnya harus sama dengan CI yang ditampilkan oleh
summary(fit_cox).
9.15.6.2 Likelihood Ratio Test (LRT)
- Bandingkan:
- Model null: tanpa kovariat (hanya \(h_0(t)\)) → log-likelihood \(\ell_0\).
- Model penuh: dengan kovariat
rx→ log-likelihood \(\ell_1\).
Statistik uji:
\[ \text{LRT} = -2(\ell_0 - \ell_1) \sim \chi^2_{(df)} \]
dengan \(df =\) banyaknya parameter kovariat (di sini 1).
loglik_null <- fit_cox$loglik[1]
loglik_full <- fit_cox$loglik[2]
LRT <- -2 * (loglik_null - loglik_full)
p_LRT <- 1 - pchisq(LRT, df = 1)
c(LRT = LRT, p_LRT = p_LRT)## LRT p_LRT
## 1.0514531 0.30517279.15.6.3 Score test (Log-rank)
Score test pada \(\beta = 0\) untuk
model 1 kovariat setara dengan uji log-rank.
R menghitungnya secara otomatis, tetapi di sini fokus kita adalah:
- Score \(U(0)\) ≈ perbedaan antara observed vs expected event di setiap risk set.
- Varian \(I(0)\) mirip dengan formula varians di uji log-rank.
9.15.7 Interpretasi Hazard Ratio
Dari output:
HR## rx
## 0.5508019CI_HR## [1] 0.174317 1.740408- Jika HR > 1: kelompok
rx = 2memiliki hazard lebih tinggi dibandingrx = 1. - Jika HR < 1: hazard lebih rendah (survival lebih baik) di
rx = 2. - Misal HR = 0.6 → hazard 40% lebih rendah; HR = 1.5 → hazard 50% lebih tinggi.
Slide memberi contoh interpretasi:
“Hazard ratio 0.6 berarti risiko pasien dengan treatment A meninggal 40%
lebih kecil dibanding treatment B.” fileciteturn0file0
9.15.8 Asumsi Model Cox PH
Ringkas dari slide: fileciteturn0file0
- Proportional hazards
Perbandingan hazard (HR) antar kelompok konstan sepanjang waktu: \[ \frac{h(t \mid \mathbf{x}_A)}{h(t \mid \mathbf{x}_B)} = \exp\big(\beta'(\mathbf{x}_A - \mathbf{x}_B)\big) \quad \forall t \] - Independensi survival time antar subjek (tidak saling bergantung).
- Kovariat diukur tanpa error dan hubungannya log-linear terhadap hazard (dalam bentuk \(\beta'x\)).
9.15.9 Ringkasan untuk Bahan Ajar
Poin-poin yang bisa langsung dipindah ke slide/handout:
Model Cox PH: \[ h(t \mid \mathbf{x}) = h_0(t)\,\exp(\beta_1 x_1 + \cdots + \beta_p x_p) \]
Hazard ratio: \[ \text{HR} = \exp(\beta_k) \] untuk kenaikan satu unit kovariat \(x_k\).
Partial likelihood: \[ L_p(\beta) = \prod_{k=1}^K \frac{\exp(\beta'\mathbf{x}_{i(k)})} {\sum_{l \in R_k} \exp(\beta'\mathbf{x}_l)} \]
Log-partial likelihood: \[ \ell_p(\beta) = \sum_{k=1}^K \left[ \beta'\mathbf{x}_{i(k)} - \log\left( \sum_{l \in R_k} \exp(\beta'\mathbf{x}_l) \right) \right] \]
Varian dan uji Wald: \[ \widehat{\text{Var}}(\hat{\beta}) = I(\hat{\beta})^{-1}, \quad Z_{\text{Wald}} = \frac{\hat{\beta}}{\text{SE}(\hat{\beta})} \]
CI 95% untuk HR: \[ \text{CI}_{\text{HR}} = \left[\exp(\hat{\beta} - 1.96\,\text{SE}),\ \exp(\hat{\beta} + 1.96\,\text{SE})\right] \]
Dengan R Markdown ini, Anda bisa:
- Menjelaskan konsep hazard dan Cox PH secara matematis.
- Menunjukkan perhitungan manual partial likelihood untuk satu kovariat.
- Menyandingkan hasil manual dengan
coxph()sehingga mahasiswa melihat jembatan antara rumus dan output software.