Survival Analysis dengan R

Survival Analysis dengan R

Halo teman-teman, kembali lagi bersama jokoding.com, sebuah situs belajar dan berbagi analisis data menggunakan aplikasi olah data. Kita kemarin telah membahas mengenai beberapa teori mengenai machine learning dalam data science, seperti supervised learning, unsupervised learning, serta reinforcement learning. Dalam kesempatan kali ini, kita coba kembali meneruskan pembahasan mengenai pemodelan termasuk bagaimana mengaplikasikannya dengan R.

Baik, kali ini kita akan membahas sebuah alat analisis yang sangat populer di bidang kesehatan teman-teman, yaitu survival analysis. Mungkin sebagian dari kita pernah mendengar kata survive atau yang artinya bertahan ini. Jadi, di survival analysis ini kita akan mengenal bagaimana cara menganalisis waktu atau lama hingga sebuah kejadian benar-benar terjadi. Dalam survival analysis, kita akan mengenal 2 outcome utama, yaitu waktu (T), serta event (e). Event ini merupakan variabel biner 0 dan 1 dengan nilai 1 menunjukkan terjadinya sebuah kejadian, sedangkan 0 tidak terjadinya sebuah kejadian.

Dalam survival analysis ini pula, kita juga akan mengenal istilah tersensor (cencored), yaitu ketika nilai dari event (e) bernilai 0. Namun yang perlu dicatat bahwa tersensor dalam hal ini tidak menunjukkan besaran peluang dari e itu sendiri.

Dalam survival analysis, fungsi yang digunakan di dalam pemodelan dikenal dengan fungsi survival dengan notasi sebagai berikut:

S(t) = P(T > t)

yang dapat diterjemahkan sebagai peluang kejadian (hidup) saat melampaui waktu sebesar t

Selain itu, dalam survival analysis kita juga akan mengenal fungsi Hazard yang dinotasikan sebagai berikut:

h(t) = lim dt --> 0 {P(t <= T <= (t + dt))/dt}

kalau kita amati, fungsi Hazard di atas merupakan fungsi kelajuan sesaat yang diterjemahkan sebagai laju penurunan peluang sebuah kejadian sesaat setelah waktu T. Dalam dunia kesehatan atau medis, kita dapat menerjemahkan seberapa besar peluang seseorang meninggal sesaat setelah ia hidup selama T satuan (hari, bulan, tahun). Kelajuan berdasarkan fungsi Hazard di atas menunjukkan sebuah grafik atau kurva yang menurun. Hal ini sesuai dengan situasi dalam bidang kesehatan, semakin tua seseorang, maka kualitas fisiknya juga mengalami penurunan hingga ia telah dinyatakan meninggal.

Di dalam survival analysis sendiri, setidaknya kita dapat berkenalan dengan 3 jenis model, yaitu Kaplan-Meier, model Eksponensial (exponential model), serta Cox Proportional Hazard (Coxph). Sekilas cara membedakan ketiga jenis model tersebut bisa kita amati dari grafiknya. Kaplan-Meier biasanya digambarkan dengan bentuk grafik fungsi tangga (floor function), sementara model eksponensial digambarkan dengan bentuk grafik kontinu yang mulus (smooth), sedangkan model Cox Proportional Hazard merupakan model dengan mengombinasikan model Kaplan-Meier dan Eksponensial.

Sebelum praktikum, agar pemahaman kita tidak tipis-tipis, ada baiknya kita melakukan sebuah simulasi manual terlebih dahulu. Misalkan kita meneliti ketahanan seseorang terkena penyakit kanker (cancer) sampai ia meninggal (e = 1). Dalam penelitian kita misalkan melakukan wawancara dengan jumlah sampel sebanyak 9 orang dan kita amati selama 5 tahun. Dari hasil wawancara tersebut, kemudian kita buat sebuah oretan sederhana sebagai berikut:

Hasil pengamatan 9 responden

Responden 1 ditemukan meninggal akibat penyakit kanker pada tahun ke-3 dari acuan awal penelitian. Kemudian responden 2 terpilih sampel pada tahun ke-2, namun pada tahun ke-3 ia dinyatakan meninggal akibat kanker. Responden ke-3, dia tidak tahu apakah dirinya mengidap kanker atau tidak, dia bergabung dalam penelitian pada tahun ke-2 dan ternyata ia meninggal pada tahun ke-4. Sedangkan responden ke-4, ia bergabung pada tahun ke-2 tetapi pada tahun ke-4 ia tidak lagi mau diwawancarai sehingga dianggap tersensor (e = 0). Responden ke-5, ia sejak awal penelitian tidak tahu pasti apakah ia mengidap kanker atau tidak, dan pada tahun ke-2 dia ternyata non-respon karena pindah tempat tinggal atau tinggal di luar area penelitian. Dan seterusnya sampai responden ke-9, ia baru bergabung pada tahun ke-4 dan ternyata dinyatakan meninggal akibat kanker di akhir masa penelitian (e = 1). Dalam hal inilah terlihat bahwa situasi tersensor tidak ada hubungan dengan peluang seseorang meninggal akibat penyakit kanker itu sendiri, karena bisa jadi ia tersensor karena pindah tempat, tidak mau diwawancarai, atau hal lain sehingga responden tidak tercatat lagi dalam selang periode penelitian.

Dari ke-9 data penelitian ini, selanjutnya kita coba tabulasikan menjadi life table survival sebagai berikut:

Membangun life table survival

Setidaknya terdapat dua penekanan dalam pembentukan life table survival. Pertama, life table hanya mengakomodir waktu-waktu dengan event (e) yang tidak tersensor secara berurutan. Bila terdapat waktu yang berurutan tersensor (e = 0), maka dapat diabaikan, misal waktu T-1 e = 0, T-2 e = 0, T-3 e = 1, maka langsung saja buat life table dari T-3. Kedua, bila kondisinya T-1 e = 1, T-2 e = 0, T-3 e = 1, maka life table dimulai dari T-1 tetapi T-2 kendati tersensor, ia tetap diperhitungkan karena berkaitan dengan jumlah responden berisiko dalam life table survival.

Setelah life table survival terbentuk, berikutnya adalah kita akan membuat grafik Kaplan-Meiernya. Dalam visualisasi grafis Kaplan-Meier, titik mulainya dari peluang survival 100 persen atau 1. Kemudian menurun seperti fungsi tangga sampai batas akhir penelitian. Adapun grafik Kaplan-Meier kasus di atas dapat kita visualisasikan berikut:

Grafik Kaplan-Meier

Itu tadi sekilas mengenai bagaimana kita melakukan survival analysis paling sederhana dengan menggunakan penghitungan manual. Kita akan coba melanjutkan dengan penerapan survival analysis dengan menggunakan R. Untuk praktik survival analysis kali ini, kita akan mengangkat topik mengenai ketahanan seseorang terhadap penyakit kanker darah atau yang biasa kita sebut sebagai leukemia. Leukemia ini kalau ingat ingat pengertian secara medis merupakan kondisi seseorang yang mengalami gangguan akibat produksi sel darah putihnya terlalu berlebih dan biasa terjadi dalam tulang sumsum. Selain terjadi pada orang dewasa, penyakit ini juga kerapkali ditemukan pada anak-anak. Saking berbahayanya, orang yang mengidap penyakit ini memiliki umur yang tidak panjang.

Dalam praktik kali ini, kita akan mencoba memodelkan survival dalam bentuk sederhana mirip dengan pendahuluan kita di atas tadi dan model survival dengan menggunakan 1 kovariat. Untuk itu, teman-teman perlu mengunduh datanya pada tautan berikut. Setelah datanya diunduh, pemodelan survival analysis dapat menggunakan beberapa code berikut:

Code:

#Instal dan aktivasi package
install.packages("survival")
library(survival)

#Import Data Leukimia
library(readxl)
leukimia <- read_excel("E:/R/Survival Analysis/leukimia.xlsx")

#Mengattach data
attach(leukimia)

#melihat Variabel dalam Data
names(leukimia)

Hasil:

# A tibble: 42 x 5
   Pasien waktu sensor   TRT   WBC
    <dbl> <dbl> <dbl> <dbl> <dbl>
1      1     6      1     1     1
2      2     6      1     1     2
3      3     6      1     1     2
4      4     7      1     1     2
5      5    10      1     1     2
6      6    13      1     1     2
7      7    16      1     1     2
8      8    22      1     1     1
9      9    23      1     1     2
10     10     6      0     1     2
# ... with 32 more rows

[1] "Pasien" "waktu" "sensor" "TRT" "WBC"

Code:

#angka 1 itu kalau tidak ada covariat atau variabel X
km.model1 <- survfit(Surv(waktu, sensor) ~ 1, data = leukimia, type = "kaplan-meier")

#Melihat ringkasan model
summary(km.model1)

Hasil:

Call: survfit(formula = Surv(waktu, sensor) ~ 1, data = leukimia, type = "kaplan-meier")

time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1     42       2    0.952 0.0329       0.8901        1.000
    2     40       2    0.905 0.0453       0.8202        0.998
    3     38       1    0.881 0.0500       0.7883        0.985
    4     37       2    0.833 0.0575       0.7279        0.954
    5     35       2    0.786 0.0633       0.6709        0.920
    6     33       3    0.714 0.0697       0.5899        0.865
    7     29       1    0.690 0.0715       0.5628        0.845
    8     28       4    0.591 0.0764       0.4588        0.762
   10     23       1    0.565 0.0773       0.4325        0.739
   11     21       2    0.512 0.0788       0.3783        0.692
   12     18       2    0.455 0.0796       0.3227        0.641
   13     16       1    0.426 0.0795       0.2958        0.615
   15     15       1    0.398 0.0791       0.2694        0.588
   16     14       1    0.369 0.0784       0.2437        0.560
   17     13       1    0.341 0.0774       0.2186        0.532
   22      9       2    0.265 0.0765       0.1507        0.467
   23      7       2    0.189 0.0710       0.0909        0.395

Code:

#Plot model tanpa covariat
plot(km.model1, conf.int = F, xlab = "waktu (t)", ylab = "%Hidup = S(t)", las = 1,
col = "blue", mark.time = T, main = "Survival Analysis Sederhana")
abline(h = 0.5, col = "red")

Hasil:

Survival Analysis tanpa kovariat

Code:

#Data leukimia dengan 1 covariat TRT
#Covariat TRT (1, 0) yang menggambarkan ada 2 grup, yaitu grup 1 dan 0
#sensor 0 berarti tidak ada kejadian kematian selama penelitian dilakukan (sensoring)
km.model2 <- survfit(Surv(waktu, sensor) ~ TRT, type = "kaplan-meier")

#Ringkasan model dengan 1 covariat TRT
summary(km.model2)

Hasil:

Call: survfit(formula = Surv(waktu, sensor) ~ TRT, type = "kaplan-meier")

                TRT=0
time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1     21       2   0.9048 0.0641      0.78754        1.000
    2     19       2   0.8095 0.0857      0.65785        0.996
    3     17       1   0.7619 0.0929      0.59988        0.968
    4     16       2   0.6667 0.1029      0.49268        0.902
    5     14       2   0.5714 0.1080      0.39455        0.828
    8     12       4   0.3810 0.1060      0.22085        0.657
   11      8       2   0.2857 0.0986      0.14529        0.562
   12      6       2   0.1905 0.0857      0.07887        0.460
   15      4       1   0.1429 0.0764      0.05011        0.407
   17      3       1   0.0952 0.0641      0.02549        0.356
   22      2       1   0.0476 0.0465      0.00703        0.322
   23      1       1   0.0000     NaN           NA           NA

                TRT=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
    6     21       3    0.857 0.0764        0.720        1.000
    7     17       1    0.807 0.0869        0.653        0.996
   10     15       1    0.753 0.0963        0.586        0.968
   13     12       1    0.690 0.1068        0.510        0.935
   16     11       1    0.627 0.1141        0.439        0.896
   22      7       1    0.538 0.1282        0.337        0.858
   23      6       1    0.448 0.1346        0.249        0.807

Code:

#Uji signifikansi dengan perlakuan dan tanpa perlakuan
#H0 tidak terdapat berbedaan fungsi survival antara dengan perlakuan dan tanpa perlakuan
#Jika alpha < 0,05 tolak H0
survdiff(Surv(waktu, sensor) ~ TRT)

Hasil:

Call:
survdiff(formula = Surv(waktu, sensor) ~ TRT)

       N Observed Expected (O-E)^2/E (O-E)^2/V
TRT=0 21       21     10.7      9.77      16.8
TRT=1 21        9     19.3      5.46      16.8

Chisq= 16.8 on 1 degrees of freedom, p= 4e-05

Terlihat bahwa dengan perlakuan menunjukkan efek positif signifikan karena mampu menambah atau memperpanjang waktu bertahan hidup seorang pengidap leukemia

Code:

#plot model dengan covariat
plot(km.model2, conf.int = F, xlab = "waktu (t)", ylab = "%Hidup = S(t)", las = 1,
col = c("red","blue"), mark.time = T, main = "Survival Analysis 1 Kovariat")
abline(h = 0.5, col = "green")
text(0,0.55, c("Median"))
legend(29,1, legend = c("tanpa perlakuan", "dengan perlakuan"),
lwd = 1, lty = 2, col = c("red", "blue"), cex = 0.6)

Hasil: