Email:
RPubs: hhttps://rpubs.com/imelda123

1 Apa itu Analisis Survival?

Analisis Survival adalah suatu studi tentang metode untuk menganalisis data, dimana variabel responnya adalah waktu sampai terjadinya suatu peristiwa tertentu. Dapat disebut juga sebagai data waktu ke peristiwa yang terdiri dari waktu mulai dan waktu akhir yang berbeda.

Analisis Survival secara umum banyak digunakan di bidang kesehatan, namun sekarang metode ini sudah banyak diaplikasikan di beberapa bidang, contohnya:

  • Bidang Studi Kanker untuk menghitung taraf kehidupan pasien,
  • Bidang Sosiologi, untuk menganalisis peristiwa sejarah yang terjadi,
  • Bidang Teknik, untuk menganalisis kegagalan yang mungkin dapat terjadi.

Dalam Studi Kanker, pertanyaan yang biasanya di tanyakan adalah sebagai berikut:

  1. Apa dampak yang paling berpengaruh terhadap jangka waktu hidup pasien?
  2. Berapa peluang pasien bertahan hidup selama 3 tahun?
  3. Adakah perbedaan jangka waktu hidup suatu pasien dengan pasien lainnya?

1.1 Contoh dalam Bidang Lainnya

Metode analisis data waktu ke peristiwa tertentu, banyak dimanfaatkan di berbagai bidang, tidak terkecuali dalam bidang:

  • Waktu seseorang terinfeksi HIV, sampai berkembang menjadi AIDS,
  • Waktu saat terjadinya serangan jantung,
  • Waktu saat seseorang menyalahgunakan suatu zat,
  • Waktu saat seseorang inisiasi aktivitas seksual,
  • Waktu untuk melihat waktu kerusakan mesin terjadi.

1.2 Nama lain dari Analisis Survival

Karena Analisis Survival banyak digunakan diberbagai bidang, maka namanya pun berbeda ditiap bidangnya, seperti:

  • Analisis Reliabilitas
  • Analisis Durasi
  • Analisis Riwayat Peristiwa
  • Analisis Waktu-ke-peristiwa
  • Model Survival (Model Bertahap Hidup)

1.3 Tujuan Analisis Survival

Tujuan dari pembelajaran ini adalah untuk mendeskripsikan konsep dasar dari analisis survival. Dalam studi kanker, kebanyakan metode analisis survival yang digunakan adalah:

  1. Plot Kaplan-Meier, untuk memvisualisasikan kurva survival (kurva kelangsungan hidup).
  2. Tes Log-Rank, untuk membandingkan kurva survival dari dua atau lebih kelompok.
  3. Regresi Proporsional hazards Cox, untuk menggambarkan pengaruh variabel terhadap survival seseorang. Model ini dibahas dalam bab berikutnya: Model proporsional hazard Cox.

Saya akan menjelaskan konsep penting dari analisis survival, termasuk:

  • Bagaimana untuk menghasilkan atau menginterpretasikan kurva survival,
  • Bagaimana untuk mengukur dan menguji perbedaan survival antara dua pasien atau lebih kelompok pasien.

Kemudian, saya akan melanjutkan dengan mendeskripsikan analisis multivariat menggunakan model proporsional hazard Cox.

2 Komponen dari Survival

Dengan \(i\) sebagai subjek:

  1. Waktu Peristiwa Terjadi \(T_i\)

  2. Waktu Sensoring \(C_i\)

  3. Indikator Peristiwa \(\delta_i\):

    • Bernilai 1 jika peristiwa berhasil diamati (observasi) \(T_i \leq C_i\)
    • Bernilai 0 jika peristiwa tidak selesai diamati (sensoring) \(T_i > C_i\)

  4. Waktu Observasi \(Y_i = \min(T_i, C_i)\)

Contohnya:

Dalam data lung, sepanjang waktu yang diamati dan indikator peristiwa disediakan.

  • Waktu : waktu bertahan hidup (waktu survival) dalam hitungan hari (\(Y_i\))
  • Status : 1 = Sensor, 2 = Meninggal (\(\delta_i\))

Dataset lung tersedia dalam package survival dalam R. Data ini berisi subjek dengan kondisi stadium lanjut dari North Central Cancer Treatment Group. Beberapa variabel yang akan saya gunakan dalam mendemostrasikan metode ini adalah:

  • Waktu : waktu survive dinilai dalam harian
  • Status : 1 = Sensor, 2 = Meninggal
  • Jenis Kelamin : 1 = Pria, 2 = Wanita
##   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       0

3 Apa itu Sensoring?

Referensi : RICH JT, NEELY JG, PANIELLO RC, VOELKER CCJ, NUSSENBAUM B, WANG EW. A PRACTICAL GUIDE TO UNDERSTANDING KAPLAN-MEIER CURVES. Otolaryngology head and neck surgery: official journal of American Academy of Otolaryngology Head and Neck Surgery. 2010;143(3):331-336. doi:10.1016/j.otohns.2010.05.007.

3.1 Tipe-tipe Sensoring

Suatu subjek dapat disensor dikarenakan:

  1. Sampel hilang (seorang pasien yang melarikan diri dari rumah sakit, ketika rumah sakit tersebut sedang meneliti pasien tersebut dan pasien lainnya)
  2. Penarikan dari studi (seorang pasien tidak lagi menjadi bahan penelitian)
  3. Tidak ada peristiwa yang terjadi hingga akhir masa penelitian (penelitian yang dilakukan telah melewati batas waktu yang ditentukan, sehingga penelitian terhadap seorang pasien bisa di hentikan)

Secara khusus, hal diatas ini adalah contoh dari sensoring di bagian kanan.

Sensoring bagian kiri dan sensoring interval juga dapat terjadi. Jenis data yang memiliki sensoring bagian kiri dan sensoring interval juga ada metode untuk menganalisisnya. Namun, analisis nya akan dibatasi pada sensoring kanan.

3.2 Data Survival Sensoring

Dalam contoh diatas, bagaimana kita akan menghitung proporsi bebas dari peristiwa dalam 10 tahun?

Subjek 6 dan 7 bebas dari peristiwa dalam 10 tahun, yang artinya subjek 6 dan 7 tidak meninggal ataupun tidak sembuh dalam penelitian yang berjangka 10 tahun.

Subjek 2,9, dan 10 memiliki peristiwa dalam 10 tahun, yang artinya subjek 2,9 dan 10 meninggal/sembuh dalam penelitian yang berjangka 10 tahun.

Subjek 1,3,4,5, dan 8 sensoring sebelum 10 tahun, artinya adalah kita tidak tahu apakah mereka meninggal atau tetap hidup dalam penelitian berjangka 10 tahun itu. Lalu yang menjadi pertanyaan adalah bagaimana kita bisa memasukkan subjek ini ke dalam perkiraan kita?

3.3 Distribusi dari Tindak Lanjut

  • Subjek yang disensor masih memberikan informasi sehingga harus dimasukkan secara tepat dalam analisis.
  • Distribusi dari tindak lanjut waktu tidak tepat, dan mungkin akan berbeda antara pasien yang disensor dengan pasien yang mengalami peristiwa (meninggal/sembuh)
  • Tindak lanjut waktu selalu positif

4 Pembahasan mengenai “waktu” dalam R

Data yang kita miliki sering sekali dalam kondisi tanggal mulai - tanggal akhir, daripada perhitungan waktu survive nya. Langkah pertama untuk kita adalah memastikan formatnya berbentuk format tanggal dalam R.

Berikut ini adalah contoh dari dataset kecil dengan variabel sx_date sebagai tanggal operasi dan last_fup_date sebagai tanggal tidak lanjut.

## # A tibble: 3 x 2
##   sx_date    last_fup_date
##   <chr>      <chr>        
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

Kita bisa melihat bahwa kedua variabel ini jenisnya adalah karakter, kejadian ini sering kali terjadi. Kita butuh untuk membuat formatnya menjadi format tanggal. Cara untuk mengubahnya adalah sebagai berikut:

## # A tibble: 3 x 2
##   sx_date    last_fup_date
##   <date>     <date>       
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

  • Catat! Dalam R format harus menyertakan pemisah serta simbol. Contohnya, jika tanggal anda dalam format m/d/Y maka anda memerlukan format = "% m /% d /% Y"

  • Klik link ini untuk melihat daftar lengkap format simbol tanggal https://www.statmethods.net/input/dates.html

Kita juga bisa menggunakan package lubridate untuk format tanggal. Kasus dibawah ini, menggunakan fungsi ymd

## # A tibble: 3 x 2
##   sx_date    last_fup_date
##   <date>     <date>       
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

  • Catat! Tidak seperti opsi dasar R, pemisah tidak perlu ditentukan

  • Halaman bantuan untuk ?Dmy akan menampilkan semua opsi format

4.1 Menghitung Waktu Survival

Setelah membuat format tanggal, kita perlu untuk menghitung perbedaan antara waktu mulai dan waktu berakhir dalam beberapa unit, biasanya dalam bulan atau tahun.

Dalam R, gunakan difftime untuk menghitung jumlah hari antara dua tanggal, dan mengubahnya menjadi nilai numerik menggunakan as.numeric. Kemudian ubah ke tahun dengan membaginya dengan 365,25 (jumlah rata-rata hari dalam setahun).

## # A tibble: 3 x 3
##   sx_date    last_fup_date os_yrs
##   <chr>      <chr>          <dbl>
## 1 2007-06-22 2017-04-15      9.82
## 2 2004-02-13 2018-07-04     14.4 
## 3 2010-10-27 2016-10-31      6.01

Menggunakan package lubridate, operator %--% akan menetapkan interval waktu, yang mana akan di ubah menjadi jumlah detik yang telah berlalu menggunakan as.duration, dan akhirnya diubah lagi menjadi tahun, dengan membaginya menggunakan dyears(1), yang memberikan jumlah detik dalam tahun.

## # A tibble: 3 x 3
##   sx_date    last_fup_date os_yrs
##   <chr>      <chr>          <dbl>
## 1 2007-06-22 2017-04-15      9.82
## 2 2004-02-13 2018-07-04     14.4 
## 3 2010-10-27 2016-10-31      6.01
  • Catat!: kita perlu package lubridate lalu menjalankannya dengan library(lubridate) untuk dapat mengakses operator khusus (lakukan hal yang sama untuk setiap kasus ingin menggunakan package tertentu)

4.2 Indikator Peristiwa

Untuk komponen dari survival data, saya menjelaskan dua indikator peristiwa:

Indikator peristiwa \(\delta_i\):

  • 1, Jika peristiwa berhasil di obervasi (i.e. \(T_i \leq C_i\))
  • 0, Jika terjadi sensoring (i.e. \(T_i > C_i\))

Namun, dalam R fungsi dari Surv dapat ditulisakan sebagai TRUE/FALSE (TRUE = peristiwa) atau bisa juga di tuliskan 1/2 (2 = peristiwa).

Dalam data lung, kita memiliki:

  • Status : 1 = Sensoring, 2 = Meninggal

4.3 Fungsi Survival

Dalam fungsi matematika peluang bahwa subjek akan survive (bertahan) melebihi waktu yang ditentukan, dirumuskan sebagai berikut:

\[S(t) = Pr(T>t) = 1 - F(t)\]

dimana:

  • \(S(t)\) : Fungsi Survival
  • \(F(t) = Pr(T \leq t)\): Fungsi Distribusi Kumulatif (CDF)

Secara teori, fungsi survive (bertahan hidup) itu mulus; dalam praktiknya saya mengamati peristiwa pada skala waktu yang berbeda.

4.4 Peluang Survival

  • Peluang Survival pada waktu tertentu (\(S(t)\)), adalah peluang bersyarat untuk survive (bertahan hidup) setelah waktu tersebut. Ingat selalu bahwa seseorang telah bertahan sebelum waktu itu.
  • Dapat di perkirakan sebagai jumlah pasien (subjek yang diteliti) tanpa ada yang kabur/ melarikan diri pada saat penelitian, lalu dibagi dengan jumlah pasien yang hidup sebelum waktu itu.
  • Estimasi Kaplan-Meier dari peluang survive (bertahan hidup) adalah produk dari peluang bersyarat hingga saat waktu itu.
  • Pada waktu 0, peluang survive adalah 1, yaitu \(S(t_0)=1\)

4.5 Membuat Object Survival

Metode Kaplan-Meier adalah cara yang paling banyak digunakan untuk memperkirakan waktu survive seseorang dan juga peluangnya. Metode ini merupakan pendekatan non-parametrik yang menghasilkan fungsi langkah, dimana ada penurunan dalam setiap kali suatu peristiwa terjadi.

  • Fungsi Surv dari package survival membuat objek survival yang digunakan sebagai respons dalam perumusan modelnya. Akan ada satu data untuk setiap subjek yaitu survival time (waktu bertahan hidup), yang di ikuti dengan + jika subjek mengalami sensoring. Mari lihat 10 observasi pertama:
##  [1]  306   455  1010+  210   883  1022+  310   361   218   166

4.6 Estimasi Kurva Survival

  • Fungsi survfit digunakan untuk membuat kurva survival berdasarkan dari formula yang ada. Mari buat kurva survival (kelangsungan hidup) secara keseluruhan untuk semua kelompok, tetapkan ke objek f1, lalu lihat names dari object tersebut:
##  [1] "n"         "time"      "n.risk"    "n.event"   "n.censor"  "surv"     
##  [7] "std.err"   "cumhaz"    "std.chaz"  "type"      "logse"     "conf.int" 
## [13] "conf.type" "lower"     "upper"     "call"

Beberapa komponen kunci dari objek survfit, akan digunakan untuk membuat kurva survival antara lain:

  • time, yang berisi titik awal dan akhir dari setiap interval waktu
  • surv, yang berisi peluang survive (bertahan hidup) yang sesuai untuk setiap time

4.7 Plot Survival Kaplan-Meier

Sekarang kita plot objek survfit dalam R untuk mendapatkan plot Kaplan-Meier.

  • Plot yang di berikan oleh R menunjukkan fungsi langkah (garis penuh) dengan interval kepercayaan (garis putus-putus).
  • Garis horizontal mewakili durasi kelangsungan hidup untuk interval tersebut, dalam plot tersebut durasi nya dalam hari.
  • Interval akan diakhiri oleh suatu peristiwa.
  • Ketinggian dari garis vertikal menunjukkan perubahan nilai kumulatif peluangnya.
  • Obervasi yang disensor, ditunjukkan dengan tanda centang, mengurangi survive kumulatif di antara interval. (Perhatikan tanda centang untuk pasien yang disensor tidak ditampilkan secara default, tetapi bisa ditambahkan menggunakan opsi mark.time = TRUE)

Sebagai alternatif, fungsi ggsurvplot dari package survminer yang dibangun diatas ggplot2 dan dapat digunakan untuk membuat plot Kaplan-Meier. Lihat cheatsheet untuk package survminer.

  • Plot default yang muncul menggunakan ggsurvplot menunjukkan fungsi langkah (garis solid) dengan pita kepercayaan yang terkait (area yang berbayang).
  • Tanda centang untuk pasien yang disensor ditampilkan secara default, dalam contoh ini agak mengaburkan garis itu sendiri, dan hal tersebut dapat ditekan menggunakan opsi censor=FALSE

4.8 Estimasi Survival x-year

Satu kuantitas yang sering menjadi perhatian dalam analisis survival adalah peluang untuk survive melebihi jumlah (\(x\)) tahun tertentu.

Sebagai contoh, untuk mengestimasi peluang survive sampai 1 tahun, gunakan summary dengan argumen times (Perhatikan variabel time di data lung sebenarnya dalam hari, jadi kita perlu menggunakan times = 365,25)

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   365     65     121    0.409  0.0358        0.345        0.486

Kita bisa melihat bahwa peluang survival (kelangsungan hidup) selama 1 tahun dalam penelitian ini adalah sebesar 41%

Batas bawah dan atas yang terkait dari interval kepercayaan 95% juga ditampilkan.

4.9 Peluang Survival x-year

Peluang survival (kelangsungan hidup) 1-tahun adalah titik pada sumbu y yang sesuai dengan 1 tahun pada sumbu x untuk kurva survival (kelangsungan hidup)

Apa yang terjadi jika kita menggunakan perkiraan “naif”?

121 dari 228 pasien meninggal dalam 1 tahun, jadi: \[\Big(1 - \frac{121}{228}\Big) \times 100 = 47\%\]

  • Kamu akan mendapatkan perkiraan yang salah tentang kemungkinan survive selama 1 tahun, ketika kamu mengabaikan fakta bahwa 42 pasien sensoring (hilang/tidak diketahui kondisinya) sebelum 1 tahun.

    • Ingat perkiraan yang benar dari kemungkinan 1 tahun untuk survive adalah 41%

5 Periksa Sensoring yang diabaikan

  • Misalkan ada 2 studi, masing-masing memiliki 228 subjek. Ada 165 kematian di setiap penelitian. Tidak ada sensor di studi 1 (garis pink), 63 pasien di sensor di studi 2 (garis biru).
  • Mengabaikan penyensoran mengarah pada perkiraan yang terlalu tinggi dari peluang survive secara keseluruhan, karena subjek yang disensor hanya menyumbangkan informasi untuk sebagian dari waktu tindak lanjut (penanganan), dan kemudian keluar dari set risiko, sehingga menurunkan peluang kumulatif untuk bertahan hidup.

6 Memperkirakan Nilai Tengah (Rata-rata) Survival

Hal penting lain yang menarik dalam analisis survival adalah waktu survive rata-rata, yang kita ukur menggunakan median. Waktu survive tidak diharapkan terdistribusi secara normal sehingga nilai meannya bukan ringkasan yang tepat.

Kita dapat memperoleh ini langsung dari objek survfit kita dengan memakai data lung:

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##       n  events  median 0.95LCL 0.95UCL 
##     228     165     310     285     363

Kita bisa melihat median waktu survive adalah 310 hari batas bawah dan atas dari interval kepercayaan 95%.

6.1 Waktu dan Peluang

Median survival adalah waktu yang sesuai dengan peluang survive 0,5:

Dalam plot diatas, kita bisa melihat bahwa rata-rata survival dari subjek yang diteliti adalah dibawah 12 bulan. Dengan kata lain, peluang seorang untuk bertahan hidup (survive) di atas 12 bulan bernilai kurang dari sama dengan 0,5 (\(\leq 0.5\))

6.2 Dampak dari mengabaikan Sensoring

Apa yang terjadi jika kita menggunakan perkiraan “naif”? Rangkuman waktu median survival di antara 165 pasien yang meninggal

##   median_surv
## 1         226
  • Kita mendapatkan perkiraan waktu median survival yang salah yaitu 226 hari, ketika kita mengabaikan fakta bahwa pasien yang disensor juga berkontribusi pada waktu tindak lanjut (penanganan).
  • Ingat bahwa perkiraan yang benar dari waktu median survival (waktu rata-rata bertahan hidup) adalah 310 hari.
  • Mengabaikan sensoring akan menciptakan kurva survival yang diturunkan secara artifisial karena waktu tindak lanjut (penanganan) yang berkontribusi pada pasien yang disensor tidak termasuk (garis pink).
  • Kurva survival yang benar untuk data lung di presentasikan oleh warna biru sebagai perbandingan terhadap kurva yang mengabaikan sensoring.

7 Membandingkan Waktu Survival

  • Kita dapat melakukan uji signifikansi antarkelompok menggunakan uji log-rank.
  • Tes log-rank memiliki bobot pengamatan yang sama selama seluruh waktu tindak lanjut (penanganan) dan merupakan cara paling umum untuk membandingkan waktu survive antar kelompok.
  • Ada versi yang lebih lengkap mengenai waktu tindak lanjut awal atau akhir yang bisa lebih sesuai tergantung pada pertanyaan penelitian (cek ?Survdiff untuk melihat opsi tes lainnya yang berbeda)

Kita akan mendapatkan nilai p-value dari tes log-rank menggunakan fungsi survdiff. Misalnya, kita dapat menguji apakah ada perbedaan waktu survive (bertahan hidup) menurut jenis kelamin di data lung

## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138      112     91.6      4.55      10.3
## sex=2  90       53     73.4      5.68      10.3
## 
##  Chisq= 10.3  on 1 degrees of freedom, p= 0.001

8 Mengekstraksi Informasi

Untuk mengekstraksi nilai p-value dari hasil survdiff sebenarnya agak merepotkan. Berikut adalah cara untuk melakukannya:

## [1] 0.001311165

Atau juga kita dapat menggunakan fungsi sdp dalam package ezfun, yang mana bisa kita install menggunakan devtools::install_github("zabore/ezfun"). Ini berguna untuk mengembalikan nilai p-value yang di format

## [1] 0.001

9 Model Regresi Cox

Kita mungkin ingin mengukur ukuran efek untuk satu variabel, atau memasukkan lebih dari satu variabel ke dalam model regresi untuk memperhitungkan efek dari beberapa variabel.

Model regresi Cox merupakan model semi parametrik yang dapat digunakan untuk menyesuaikan model regresi univariabel dan multivariabel yang memiliki hasil survival.

\[h(t|X_i) = h_0(t) \exp(\beta_1 X_{i1} + \cdots + \beta_p X_{ip})\]

  • \(h(t)\) : Hazard (kecepatan sesaat dimana peristiwa terjadi)
  • \(h_0(t)\) : Sumber awal bahaya

Beberapa poin asumsi dari model tersebut:

  • Sensoring tidak informatif
  • Bahaya yang proporsional

Catat! : model regresi parametrik untuk hasil survival juga tersedia, tetapi model tersebut tidak akan dibahas dalam pembahasan ini.

Kita dapat menyesuaikan model regresi untuk data survival menggunakan fungsi coxph, yang mengambil objek Surv di sisi kiri dan memiliki sintaks standar untuk rumus regresi dalam R di sisi kanan.

## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
## 
##        coef exp(coef) se(coef)      z       p
## sex -0.5310    0.5880   0.1672 -3.176 0.00149
## 
## Likelihood ratio test=10.63  on 1 df, p=0.001111
## n= 228, number of events= 165

10 Format Regresi Cox

Kita dapat melihat output versi tidy menggunakan fungsi tidy dari package broom:

term estimate std.error statistic p.value
sex 0.5880028 0.1671786 -3.176385 0.0014912

Atau kita bisa menggunakan tbl_regression dari package gtsummary:

Characteristic HR1 95% CI1 p-value
sex 0.59 0.42, 0.82 0.001

1 HR = Hazard Ratio, CI = Confidence Interval

11 Rasio Hazard (Rasio Bahaya)

  • Kuantitas yang menarik dari model regresi Cox adalah rasio bahaya/ hazard ratio(HR). HR mewakili rasio bahaya antara dua kelompok pada titik waktu tertentu.

    • HR diartikan sebagai tingkat peristiwa sesaat dari peristiwa yang menarik bagi mereka yang masih berisiko untuk peristiwa tersebut. Ini bukan risiko, meskipun secara umum ditafsirkan seperti itu

    • Jika kita memiliki parameter regresi \(\beta\) (dari kolom estimate dalam coxph kami) maka HR = exp (\(\beta\)).

    • HR bernilai < 1 menunjukkan pengurangan bahaya kematian sedangkan HR bernilai > 1 menunjukkan peningkatan bahaya kematian.

    • Jadi HR kita adalah = 0,59 menyiratkan bahwa sekitar 0,6 kali lebih banyak wanita yang meninggal daripada pria, pada waktu tertentu.