Pemodelan Hazard Proporsional dengan Perkalian Gamma Frailty Menggunakan Pendekatan Bayesian 1 1,2,3
Ismi Try Amalia Jaya, 2Armin Lawi, dan 3Sri Astuti Thamrin
Jurusan Matematika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Universitas Hasanuddin, Makassar, Indonesia
Abstrak Analisis survival merupakan analisis yang mempelajari waktu kelangsungan hidup suatu individu terhadap suatu kejadian. Salah satu model yang sering digunakan dalam analisis survival adalah model frailty. Model frailty merupakan perluasan dari model Cox proportional hazard. Penelitian ini bertujuan untuk mendapatkan model hazard proporsional dengan perkalian gamma frailty melalui penaksiran parameter dengan pendekatan Bayesian. Untuk mendapatkan nilai estimasi parameter digunakan algoritma Gibbs Sampling. Data yang digunakan berasal dari penelitian Mc. Gilchrist dan Aisbett tentang kekambuhan penyakit ginjal. Hasil yang diperoleh menunjukkan bahwa pasien wanita memiliki peluang yang lebih kecil terkena infeksi ginjal dibandingkan pasien pria. Kata kunci: Analisis survival, Bayesian, Distribusi Weibull, Gamma frailty, Gibbs sampling, Hazard proporsional.
1. Pendahuluan Analisis yang digunakan untuk menganalisis data waktu hidup disebut analisis data kelangsungan hidup (survival). Seringkali dalam analisis data survival, waktu kelangsungan hidup dalam kelompok yang sama saling berkorelasi karena kovariat tidak teramati. Oleh karena itu, dikembangkan sebuah pemodelan data yang merupakan salah satu cara agar kovariat dapat dimasukkan dalam model yaitu dengan mempertimbangkan kovariat tersebut sebagai kelemahan (frailty), (Clayton(1978), Oakes(1982), serta Clayton dan Cuzick(1985)). Model hazard dasar yang diusulkan adalah model hazard dasar berdistribusi Weibull karena distribusi Weibull dapat menyajikan keakuratan kegagalan dengan sampel yang sangat kecil (Hossain dan Zimmer, 2003). Distribusi gamma pada frailty dibutuhkan agar parameter dari model dapat diidentifikasi (Ibrahim et al., 2001). Oleh karena itu, dalam penulisan ini akan dikaji model hazard dasar yang diberikan frailty pada data survival multivariat. Model waktu kelangsungan hidup (survival times) yang akan digunakan untuk fungsi hazard dasar berdistribusi Weibull dengan perkalian gamma frailty.
2. Tinjauan Pustaka 2.1 Analisis Kelangsungan Hidup (Survival Analysis) Analisis survival berkaitan erat dengan waktu kelangsungan hidup (survival time). Waktu kelangsungan hidup merupakan jangka waktu dari awal pengamatan sampai terjadinya suatu peristiwa. Peristiwa itu dapat berupa kegagalan, kematian, respon, timbulnya gejala dan lain-lain (Lee, 1992). Akan tetapi, sering dijumpai suatu individu tidak mengalami kegagalan sampai batas waktu penelitian. Hal ini mengakibatkan ketidaklengkapan data kegagalan (failure time) yang sering disebut sensoring (censoring). Menurut Cox dan Oakes (1984), terdapat tiga hal yang harus diperhatikan dalam menentukan waktu survival secara tepat: a. Waktu awal yang tidak ambigu yang berarti tidak ada dua pengertian atau lebih. b. Definisi terjadinya kegagalan secara keseluruhan harus jelas. c. Skala waktu sebagai satuan pengukuran harus jelas. Pada analisis kelangsungan hidup, ada 2 hal yang mendasar yaitu fungsi survival, S(t) dan fungsi hazard, h(t). Fungsi survival merupakan dasar dari analisis ini, karena meliputi probabilitas kelangsungan hidup dari waktu yang berbeda-beda yang memberikan informasi penting tentang data survival. Berbeda dengan fungsi survival yang fokus pada tidak terjadinya peristiwa, fungsi hazard fokus pada terjadinya peristiwa. Oleh karena itu fungsi hazard dapat dipandang sebagai pemberi informasi yang berlawanan dengan fungsi survival. 2.2 Model Gamma Frailty Model frailty adalah perluasan dari model Cox hazard proporsional yang dikenal sebagai Model Cox (Cox, 1972). Model Cox proportional hazard diberikan sebagai berikut : ( )
( )
(
)
(2.1)
dimana ( ) merupakan fungsi hazard untuk setiap objek dengan variabel penduga Z , sehingga ( ) dikatakan fungsi hazard dasar. (Collet, 1994), adalah parameter regresi, dan adalah kovariat. Distribusi Gamma merupakan distribusi yang digunakan pada pemodelan fungsi hazard yang mengandung frailty karena distribusi Gamma dianggap sebagai distribusi yang paling umum digunakan (Sahu et al., 1997). Dengan mengasumsikan bahwa waktu survival dari subjek ke-i (i=1,…,n), dalam kelompok ke-j (j=1,…,m), dinotasikan dengan tij dan diberikan wi (untuk kelompok ke-i) sebagai parameter frailty yang tidak teramati, maka model Gamma frailty diberikan sebagai berikut :
(
|
)
(
)
(
)
(2.2)
dimana:
( )
= resiko gagal pada waktu tertentu dan merupakan vektor kovariat p x 1untuk subyek ke- dalam kelompok ke- , dan waktunya terikat. = vektor p x 1 dari koefisien regresi yang tidak diketahui = fungsi hazard dasar (baseline hazard function) ( ), yang berdistribusi Weibull diberikan sebagai
Fungsi hazard dasar, berikut: (
)
Dimana adalah parameter bentuk, waktu kelangsungan hidup.
(2.3) adalah parameter skala, dan
adalah
2.3 Pendekatan Bayesian Analisis Bayesian memiliki keuntungan dalam melakukan pemodelan statistik dan analisis data. Dalam pendekatan klasik, parameter merupakan suatu nilai yang konstan, dimana parameter adalah sebuah peluang sampel dari populasi (Raudenbush et al., 2002). Pada estimasi parameter dengan menggunakan metode Bayesian, terdapat tiga komponen, yaitu distribusi prior, fungsi likelihood, dan distribusi posterior (Nurmalitasari, 2011). Menurut Soejoeti dan Soebanar (1988), distribusi posterior dapat ditulis dalam bentuk distribusi prior dan fungsi likelihood sebagaimana diberikan sebagai berikut: ( ) ( ) ( ) (2.4) ( ) ( ) dimana merupakan fungsi likelihood dan merupakan distribusi prior. 2.3.1
Gibbs Sampling Gibbs Sampling diperkenalkan oleh Geman dan Geman (1984). Gibbs Sampling adalahsuatu teknik simulasi untuk membangkitkan variabel acak dari suatu distribusi tertentu secara langsung, tanpa perlu menghitung densitasnya (Casella, et al., 1992). 2.3.2
WinBugs WinBUGS adalah sebuah perangkat lunak untuk melakukan analisis Bayesian menggunakan metode Markov Chain Monte Carlo(MCMC). WinBUGS merupakan sebuah perangkat lunak berbasis BUGS (Bayesian Using Gibbs Sampling). 3. Metodologi Penelitian Data yang digunakan untuk penelitian ini adalah data waktu kekambuhan infeksi pada titik penyisipan kateter untuk pasien ginjal menggunakan peralatan dialisi portable yang bersumber dari penelitian McGilchrist dan Aisbett (1991). Data dapat diunggah pada situs http://lib.stat.cmu.edu/datasets/kidney
Adapun langkah-langkah yang akan dilakukan dalam penelitian ini adalah : 1. Menentukan fungsi hazard dasar menggunakan distribusi Weibull dengan ( ) menggunakan rumus ( ) , dimana f(t) adalah fungsi kepadatan ( ) peluang dari distribusi Weibull dan S(t) adalah fungsi survival dari distribusi Weibull. 2. Memodelkan fungsi hazard yang diberikan frailty dengan menambahkan parameter frailty, , kedalam model Cox proporsional hazard. 3. Menentukan distribusi prior menggunakan pemilihan prior noninformatif. 4. Menentukan likelihood dari fungsi hazard yang dibangun berdasarkan model proporsional hazard yang diberi frailty dengan distribusi Weibull. 5. Menentukan distribusi posterior dengan perkalian proporsional distribusi prior dengan fungsi likelihood 6. Mengestimasi parameter dari posterior yang diperoleh menggunakan algoritma Gibbs Sampling. 7. Menerapkan model proporsional hazard berdistribusi Weibull dengan frailty pada data waktu kekambuhan infeksi pada titik penyisipan kateter untuk pasien ginjal. 4. Hasil dan Pembahasan Berdasarkan data pada penelitian McGilchrist dan Aisbett (1991), diketahui jumlah pasien yang diobservasi berjumlah 38 orang. adalah waktu terjadinya failure event dalam skala hari. Status pasien pada akhir penelitian bernilai 1 jika mengalami kejadian (kekambuhan) dan bernilai 0 jika tidak mengalami kejadian.Untuk jenis kelamin diberi nilai 1 untuk laki-laki dan nilai 2 untuk perempuan. Kejadian atau penyakit yang diderita pasien diberi nilai 0 untuk jenis penyakit Glomerulo Nephritis, bernilai 1 untuk jenis penyakit Acute Nephritis, bernilai 2 untuk jenis penyakit Polycystic Kidney Disease dan bernilai 3 untuk yang lainnya. Pengolahan data pada penelitian ini menggunakan bantuan software WinBUGS 1.4.3. Iterasi pertama dilakukan sebanyak 1000 iterasi, lalu dilanjutkan hingga mencapai konvergen yaitu sebanyak 100.000 iterasi. Model dikatakan konvergen pada iterasi ke-100.000 dengan melihat pada Gambar 4.1.
(a) (b) Gambar 4.1. (a) Trace Plot βage, (b) Trace Plot βsex, (c) Trace Plot σ2
(c) Gambar 4.1. (a) Trace Plot βage, (b) Trace Plot βsex, (c) Trace Plot σ2 Pada Gambar 4.1 menunjukkan sudah tidak adanya lagi pola yang terbentuk , baik itu trend naik maupun trendturun, pada grafik trace plot. Selain melihat pada trace plot, konvergensi model juga dapat ditunjukkan melalui gambar history plot. Jika hasil history plot tampak rapat dan mampu merespon semua parameter yang diestimasi, maka model telah konvergen pada nilai mean. History plot ditunjukkan pada Gambar 4.2.
(a)
(b)
(c) Gambar 4.2. (a) History Plot βage,(b)History Plot βsex, (c) History Plot σ2
Berdasarkan Gambar 4.2 terlihat bahwa history plot tiap parameter telah rapat dan dapat merespon semua parameter. Selanjutnya untuk lebih memastikan bahwa model telah konvergen, dilihat pula grafik density tiap parameter. Gambar density tiap parameter ditunjukkan pada Gambar 4.3.
(a)
(b)
(c) Gambar 4.3 (a) Grafik Density βage, (b) Grafik Density βsex, (c) Grafik Density σ2 Gambar 4.3 menunjukkan grafik density tiap parameter telah halus sehingga model dikatakan konvergen. Setelah kondisi konvergen telah terpenuhi, maka langkah selanjutnya adalah mencari nilai estimasi tiap parameter. Hasil estimasi tiap parameter ditunjukkan pada Tabel 4.1. Tabel 4.1. Hasil estimasi parameter
,
,dan
untuk data waktu
kekambuhan infeksi pada titik penyisipan kateter untuk pasien ginjal. Parameter
Rata-Rata Posterior
Standar Deviasi
(95% Credible Interval)
-1.924
0.5875
(-2.929, -0.9938)
0.008024
0.01326
(-0.01365, 0.02984)
0.6003
0.317
(0.1589, 1.184)
Dari Tabel 4.1 diperoleh hasil taksiran parameter adalah -1.924. Berdasarkan rata-rata nilai Gibbs sampler yang diperoleh setelah 100.000 iterasi, dihasilkan fungsi densitas pada Gambar 4.3(a) sehingga interval kredibel 95% untuk taksiran adalah (-2.929, -0.9938). Sedangkan hasil taksiran parameter
adalah 0.008024. Berdasarkan rata-rata nilai Gibbs sampler yang diperoleh setelah 100.000 iterasi, dihasilkan fungsi densitas pada Gambar 4.3(b) sehingga interval kredibel 95% untuk taksiran parameter adalah (-0.01365, 0.02984). Hasil taksiran yang diperoleh untuk parameter adalah 0.6003. Berdasarkan ratarata nilai Gibbs sampler yang diperoleh setelah 100.000 iterasi pada Tabel 4.1, dihasilkan fungsi densitas pada Gambar 4.3(c) sehingga interval kredibel 95% untuk taksiran parameter adalah (0.1589, 1.184). Berdasarkan Tabel 4.1, terlihat bahwa variabel yang mempengaruhi waktu terjadinya infeksi adalah . Hal ini dapat diamati dari nilai interval kepercayaan untuk masing-masing parameter. Parameter yang interval kepercayaannya mengandung nilai nol tidak mempengaruhi waktu terjadinya infeksi. Sebaliknya, jika nilai interval kepercayaannya tidak mengandung nilai nol maka mempengaruhi waktu terjadinya infeksi seorang pasien. Nilai negatif pada interval kepercayaan yang ditunjukkan pula oleh grafik density pada Gambar 4.3, cenderung ke sebelah kiri 0, artinya pasien perempuan memiliki resiko yang lebih kecil terkena infeksi dibandingkan dengan pasien laki-laki. Nilai yang lebih dari 0 yang berada pada interval (0.1589, 1.184) menujukkan bahwa adanya heterogenitas pada data sehingga mempengaruhi waktu kelangsungan hidup pasien. Nilai estimasi parameter , dan ditunjukkan pada Tabel 4.2. Tabel 4.2. Hasil estimasi parameter dan untuk data waktu kekambuhan infeksi pada titik penyisipan kateter untuk pasien ginjal Rata-Rata Parameter Standar Deviasi (95% CI) Posterior 1.235
0.1724
(0.9598, 1.526)
0.01651
0.01587
(0.00278, 0.04533)
Pada Tabel 4.2, diperoleh hasil taksiran parameter adalah 1.235 dengan interval kredibel untuk taksiran parameter adalah (0.9598, 1.526), sedangkan hasil taksiran parameter adalah 0.01651 dengan interval kredibel parameter adalah (0.00278, 0.04533). Diketahui parameter merupakan parameter bentuk dari distribusi Weibull, sehingga berdasarkan Tabel 4.2 tersebut, nilai menujukkan bahwa laju kegagalan meningkat seiring berjalannya waktu. Dari Tabel 4.1 diketahui nilai taksiran untuk parameter βsex dan βage, sehingga jika disubtitusi kedalam Persamaan (4.1), maka akan diperoleh model hazard proporsional dengan perkalian gamma frailty berikut : (
|
)
(
)
(
)
)
5. Kesimpulan Berdasarkan hasil penulisan ini, maka dapat disimpulkan bahwa: 1. Nilai estimasi untuk paramaeter βsex adalah -1.924 dan βage adalah 0.008024, sehingga diperoleh model hazard proporsional dengan perkalian Gamma frailty berikut : ( | ) ( ) ( ) ) 2. Estimasi parameter pada model menggunakan pendekatan Bayesian dilakukan melalui tiga tahap. Tahap pertama ialah menentukan distribusi awal sebelum melakukan analisis data atau yang disebut distribusi prior. Prior yang digunakan merupakan prior non-informatif. Tahap selanjutnya ialah mengkontruksi fungsi likelihood yang dibentuk melalui model hazard proporsional dengan perkalian gamma frailty. Tahap akhir ialah menentukan distribusi posterior dengan cara proporsional fungsi likelihood dikalikan dengan distribusi prior. 3. Dari hasil pengolahan data diperoleh bahwa hanya jenis kelamin dan frailty yang mempengaruhi waktu kelangsungan hidup pasien. Untuk jenis kelamin, diketahui bahwa pasien wanita memiliki resiko yang lebih rendah terkena infeksi dibanding pasien pria. Sedangkan untuk variansi frailty yang bernilai lebih dari nol mengindikasikan bahwa terdapat heterogenitas pada data yang dapat mempengaruhi waktu kelangsungan hidup pasien. 6. Daftar Pustaka Anifa, Moch. Abdul Mukid, Agus Rusgiyono. 2012. Simulasi Stokastik Menggunakan Algoritma Gibbs Sampling. Jurnal Gaussian.1 : 21-30. Bain, L. J. dan Engelhardt, M. 1992. Introduction to Probability and Mathematical Statistics Second Edition. Duxbury Press. California. Box, G.E.P dan Tiao, G.C. 1973.Bayesian Inference In Statistical Analysis. Addision-Wesley Publishing Company, Inc; Philippines Casella, G. dan George, E.I. 1992.Explaining Gibbs Sampler.Journal of The American Statistical Association, vol. 46.167 – 174. Clayton, D. G. 1978. AModel for Association in Bivariate Life Tables and its Application in Epidemiological Studies of Familial Tendency in Chronic Disease Incidence. Biometrika,65: 141-15. Clayton, D. G. dan J. Cuzick. 1985. Multivariate Generalizations of the Proportional Hazards Model (with discussion). Journal of Royal Statistical Society Ser. A,148:82-117. Cox, D. R. 1972. Regression Models and Life Tables (with discussion), Journal of The Royal Statistical Society B,34:187-220. Cox dan Oakes, D. 1984. Analysis of Survival Data. Chapman and Hall. London. Ebeling, Charles. 1997. An Introduction to reliability and maintainability Engineering. Singspore : The MC. Graw Hill Companier Inc, New York. Evans, M., Hastings, N. dan Peacock, B. 2000. Statistical Distribution. John Wiley and Sons, Inc. United State of America.
Geman, S. dan Geman, D. 1984. Stochastic Relaxation, Gibbs Distribution, and the Bayesian Restoration of Images. IEE Transaction on Pattern Analysis and Machine Intelligence. Gilks, W.R., Richardson, S. dan Spiegelhalter, D. J. 1996. Markov Chain Monte Carlo in Practice.Chapman and Hall. New York. Gill, Jeff.2002. Bayesian Methods: A Social dan Behavioral Sciences Approach. A CRC Press Company, Washington D.C. Hossain,A.M., Zimmer, W.J., 2003. Comparison of estimation methods for Weibull parameters: complete and censored samples. J. Statist. Comput. Simulation, 73 (2), 145–153. Ibrahim, J,.Chen, M dan Sinha, D. 2001. Bayesian Survival Analysis. Springer Verlag, New York. Klein, J.P. dan Moeschberger, M.L. 2003. Survival Anlysis : Techniques for Censored dan Truncated Data. Springer, New York. Kleinbaum, D. G. dan Klein, M. 2005. Survival Analysis. Springer Science Business Media, Inc. New York. Lee, E.T. 1992. Statistical Methods for Survival Data Analysis. John Willey and Sons, Inc. New York. Le, C. T. 2003. Introductory Biostatistics. John Wiley and Sons, Inc. Lunn, D. et al. 2003. WinBugs User Manual, Version 1.4. Department of Epidemiologi & Public Health, London : http://www.mrc-bsu.cam.ac.uk/bugs Nurmalitasari. 2011. Estimasi Parameter Model Inar(1) Menggunakan Metode Bayes. Skripsi. Surakarta: Fakultas Matematika dan Ilmu Pengetahuan Alam, Universitas Sebelas Maret. Ntzoufras, I. 2009. Bayesian Modeling Using WinBUGS. Greece. John Wiley. Oakes, D. 1989. Bivariate Survival Models Induced by Frailties. Journal of the American Statistical Assosiation, 84, 487493. Raudenbush, S.W. dan Bryk, A.S. 2002.Hierarchical Linear Models, Applications and Data Analysis Methods, Second Edition. Sage Publication. Sahu, S.K., Dey, D.K. dan Aslanidou, H., Sinha, D. 1997. A Weibull Regression Model Model with Gamma Frailty for Multivariate Survival Data. Life Time Data Analysis, vol. 3. 123-137. Soejoeti, Z dan Soebanar. 1988. Inferensi Bayesian. Karunia Universitas Terbuka. Jakarta. Vaupel, J.W. Manton, K.G. Stallard, E. 1979. The Impact of Heterogeneity in Individual Frailty on the Dynamics of Mortality. Demography 16, 439 – 454.
