Email :
RPubs : https://rpubs.com/Nicholas321
instagram : https://www.instagram.com/nicholasandrian/
linkedin : https://www.linkedin.com/in/nicholas-andrian/

1 Apa itu Analisis Survival?

Analisis Survival merupakan sebuah pendekatan statistik yang digunakan untuk menyelidiki berapa lama waktu yang dibutuhkan untuk suatu peristiwa akan benar-benar terjadi. Atau disebut juga sebagai Waktu Menuju Peristiwa (Time-to-event) yang terdiri dari waktu awal analisis dan waktu akhir analisis.

Analisis Survival dapat digunakan di berbagai bidang, seperti:

  • Menganalisis waktu bertahan hidup dari pasien kanker,
  • Menganalisis peristiwa sejarah dari bidang Sosiologi,
  • Menganalisis waktu terjadinya sebuah kerusakan dari bidang teknik.

Beberapa pertanyaan di dalam suatu penelitian mengenai kanker adalah sebagai berikut:

  • Apa dampak dari perbedaan karakteristik klinis tertentu terhadap kelangsungan hidup pasien?
  • Berapa kemungkinan seorang pasien akan bertahan hidup selama 3 tahun?
  • Adakah perbedaan waktu kelangsungan hidup antara kelompok pasien?

1.1 Beberapa Contoh Bidang Lainnya

Banyak bidang umum yang memiliki data Waktu Menuju Peristiwa, contoh lainnya sepertinya:

  • Waktu dari terinfeksi HIV sampai berkembang menjadi AIDS,
  • Waktu akan terjadinya serangan jantung,
  • Waktu timbulnya penyelahgunaan zat
  • Waktu untuk memulai aktivitas seksual,
  • Waktu akan terjadinya kerusakan mesin.

1.2 Nama Lain untuk Analisis Survival

Karena analisis survival ini umum digunakan di berbagai bidang lain, maka analisis survival ini juga memiliki berbagai nama:

  • Analisis Reliabilitas
  • Analisis Durasi
  • Analisis Riwayat Kejadian
  • Analisis Waktu Menuju Peristiwa (Time-to-event)
  • Model Survival

1.3 Tujuan Analisis Survival

Tujuan dari bab ini adalah untuk menggambarkan konsep dasar dari analisis survival. Dalam studi mengenai kanker, sebagian besar analisis survival menggunakan metode berikut:

  • Plot Kaplan-Meier untuk memvisualisasikan kurva survival
  • Uji log-rank untuk membandingkan kurva survival dari dua kelompok atau lebih
  • Regresi proporsional hazard Cox untuk menggambarkan pengaruh variabel terhadap survival. Model Cox akan dibahas dalam bab berikutnya: Model proporsional hazard Cox.

Di sini, kita akan memulainya dengan menjelaskan konsep penting dari analisis survival, seperti:

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

Kemudian, kita akan melanjutkannya dengan menggambarkan analisis multivariat menggunakan model proporsional hazard Cox.

2 Komponen Survival

Untuk subjek \(i\):

  1. Waktu kejadian \(T_i\)
  2. Waktu censoring \(C_i\)
  3. Indikator kejadian \(δ_i\):
  • 1 jika terjadi suatu kejadian (yaitu \(T_i\)≤\(C_i\))
  • 0 jika sensor (yaitu \(T_i\)>\(C_i\))
  1. Waktu yang diamati \(Y_i=min(T_i,C_i)\)

Waktu yang diamati dan indikator kejadian tersedia pada data lung.

  • time: Waktu bertahan hidup dalam hari \((Y_i)\)
  • status: status censoring, 1 = sensor, 2 = meninggal \((δ_i)\)

Data set lung tersedia di R dalam package survival. Data ini berisi subjek dengan kanker paru-paru stadium lanjut dari North Central Cancer Treatment Group. Beberapa variabel yang akan kita gunakan untuk mendemonstrasikan metode ini adalah:

  • time: Waktu bertahan hidup dalam hari
  • status: status censoring, 1 = sensor, 2 = meninggal
  • sex: 1 = Laki-laki, 2 = Perempuan

3 Apa itu Censoring?

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 Jenis-jenis Censoring

Sebuah subjek kemungkinan terjadi sensor karena:

  • Tidak di follow-up
  • Ditarik dari penelitian
  • Tidak terjadi apa-apa hingga akhir masa penelitian yang ditentukan

Contoh di atas merupakan contoh sensor sebelah kanan.

Sensor sebelah kiri dan sensor interval juga kemungkinan dapat terjadi, dan ada metode untuk menganalisis jenis data ini, tetapi pada latihan kali ini hanya dibatasi pada sensor sebelah kanan.

Dalam contoh ini, bagaimana cara kita menghitung proporsi yang tidak mengalami kejadian setelah 10 tahun waktu penelitian?

Subjek 6 dan 7 tidak mengalami kejadian setelah 10 tahun. Subjek 2, 9, dan 10 mengalami kejadian sebelum 10 tahun. Subjek 1, 3, 4, 5, dan 8 terjadi sensor sebelum 10 tahun, jadi kita tidak tahu apakah mereka mengalami kejadian atau tidak selama 10 tahun - bagaimana cara untuk memasukkan subjek ini ke dalam perkiraan kita?

3.2 Distribusi dari Follow-up

  • Subjek yang mengalami sensor masih dapat memberikan informasi sehingga harus dimasukkan secara tepat dalam analisis
  • Distribusi waktu follow-up mengalami kecondongan, dan mungkin berbeda antara pasien yang mengalami sensor dengan pasien yang mengalami kejadian (meninggal)
  • Waktu follow-up selalu positif

4 Membuat Tanggal di R

Biasanya data tidak akan menampilkan secara langsung waktu survival yang telah dihitung sebelumnya, melainkan data akan sering muncul dengan tanggal mulai dan tanggal akhir dari survival. Langkah pertamanya adalah memastikan tanggal ini dirofmat sebagai dates di R.

Mari kita buat contoh data set kecil dengan variabel sx_date untuk tanggal operasi dan last_fup_date untuk tanggal terakhir follow-up.

## # 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 perhatikan keduanya adalah variabel karakter, tetapi kita membutuhkannya dalam bentuk format tanggal. Cara lainnya untuk memformat tanggal menggunakan R:

## # 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
  • Perhatikan bahwa format tanggal di R harus menyertakan pemisah dalam bentuk simbol, misal. jika tanggal Anda dalam format month/day/year, maka Anda harus menggunakan "format = %m/%d/%Y`
  • Lihat daftarlengkap simbol untuk format tanggal di https://www.statmethods.net/input/dates.html

Kita juga dapat menggunakan package lubridate untuk memformat tanggal. Dalam kasus ini, gunakan 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
  • Perhatikan bahwa tidak seperti opsi dasar R, pemisah tidak perlu disertakan.
  • Halaman bantuan untuk ?dmy akan menampilkan semua opsi format.

4.1 Menghitung Waktu Survival

Setelah memformat tangga, sekarang kita perlu menghitung perbedaan antara waktu mulai dan berakhir dalam beberapa unit, biasanya dalam bulan atau tahun. Dalam basis R, gunakan difftime untuk menghitung jumlah hari antara dua tanggal mulai dan akhir serta 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 kemudian diubah menjadi jumlah detik yang telah berlalu menggunakan as.duration dan yang terakhir diubah menjadi tahun dengan membaginya dengan dyears(1) , yang memberikan jumlah detik dalam satu 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
  • Catatan: kita perlu mengaktifkan package lubridate dengan menggunakan library agar dapat mengakses operator khusus (mirip dengan penggunaan pipa)

4.2 Indikator Kejadian

Untuk komponen data survival, indikator kejadian \(δ_i\):

  • 1 jika terjadi suatu kejadian (yaitu \(T_i≤C_i\))
  • 0 jika sensor (yaitu \(T_i>C_i\))

Namun, di R fungsi Surv juga akan menerima TRUE/FALSE (TRUE = event) atau 1/2 (2 = event).

Dalam data lung, kita memiliki:

  • status: status censoring, 1 = sensor, 2 = meninggal

4.3 Fungsi Survival

Probabilitas bahwa subjek akan bertahan melebihi waktu tertentu

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

\(S(t)\): fungsi survival

\(F(t)=Pr(T≤t)\): fungsi distribusi kumulatif

Secara teori, fungsi survival itu mulus; dalam praktiknya kita akan mengamati peristiwa pada skala waktu yang berbeda.

4.4 Probabilitas Survival

  • Probabilitas survival pada waktu tertentu, \(S(t)\) , adalah probabilitas bersyarat untuk survival setelah waktu tersebut, mengingat bahwa seseorang telah bertahan sebelum waktu itu.

  • Dapat diperkirakan sebagai jumlah pasien yang hidup tanpa kehilangan follow-up pada saat itu, dibagi dengan jumlah pasien yang hidup sebelum waktu tersebut.

  • Estimasi Kaplan-Meier dari probabilitas survival adalah produk dari probabilitas bersyarat hingga saat itu.

  • Pada waktu 0, probabilitas survival adalah 1, yaitu \(S(t0)=1\).

4.5 Membuat Objek Survival

Metode Kaplan-Meier adalah cara paling umum untuk memperkirakan waktu kelangsungan hidup dan probabilitas. Ini adalah pendekatan non-parametrik yang menghasilkan fungsi langkah, di mana ada penurunan setiap kali suatu peristiwa terjadi.

Fungsi Surv dari package survival membuat objek survival untuk digunakan sebagai respons dalam rumus model. Akan ada satu entri untuk setiap subjek yaitu survival time, yang diikuti dengan + jika subjek mengalami sensor. Mari kita lihat 10 observasi pertama:

##  [1]  306   455  1010+  210   883  1022+  310   361   218   166

4.6 Perkiraan Kurva Survival

Fungsi survfit membuat kurva survival berdasarkan rumus. Mari kita buat kurva survival keseluruhan untuk seluruh kelompok, tetapkan ke objek f1, dan lihat names objek tersebut:

##  [1] "n"          "time"       "n.risk"     "n.event"    "n.censor"  
##  [6] "surv"       "std.err"    "cumhaz"     "std.chaz"   "start.time"
## [11] "type"       "logse"      "conf.int"   "conf.type"  "lower"     
## [16] "upper"      "call"

Beberapa komponen kunci dari objek survfit ini yang akan digunakan untuk membuat kurva survival, yaitu:

  • time, yang berisi titik awal dan akhir setiap interval waktu
  • surv, yang berisi probabilitas survival yang sesuai untuk setiap time.

4.7 Plot Survival Kaplan-Meier

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

  • Plot default di basis R menunjukkan fungsi langkah (garis utuh) dengan interval kepercayaan terkait (garis putus-putus)
  • Garis horizontal mewakili durasi survival untuk interval tersebut
  • Interval diakhiri oleh suatu kejadian
  • Ketinggian garis vertikal menunjukkan perubahan probabilitas kumulatif
  • Pengamatan yang mengalami sensor, ditunjukkan dengan tanda centang, mengurangi survival 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 berasal dari ggplot2, dan dapat digunakan untuk membuat plot Kaplan-Meier. Lihatlah cheatsheet ini untuk package survminer.

Plot default menggunakan ggsurvplot menunjukkan fungsi langkah (garis utuh) dengan pita kepercayaan terkait (area berbayang).

Tanda centang untuk pasien yang mengalami sensor ditampilkan secara default, agak mengaburkan garis itu sendiri dalam contoh ini, dan dapat ditekan menggunakan opsi censor = FALSE.

4.8 Perkiraan Survival x-tahun

Satu kuantitas yang sering menjadi perhatian dalam analisis kelangsungan hidup adalah probabilitas untuk bertahan hidup melebihi jumlah (x) tahun tertentu.

Misalnya, untuk mengestimasi probabilitas bertahan hingga 1 tahun, gunakan summary dengan argumen times (Perhatikan variabel times 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 temukan bahwa probabilitas kelangsungan hidup 1 tahun dalam penelitian ini adalah 41%.

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

4.9 Kemungkinan Survival x-tahun

Probabilitas survival 1 tahun adalah titik pada sumbu y yang sesuai dengan 1 tahun pada sumbu x untuk kurva survival.

Apa yang terjadi jika Anda menggunakan perkiraan “naive”?

121 dari 228 pasien meninggal dalam 1 tahun, jadi:

\(\left(1-\frac{121}{228} \right)\times 100 = 47\%\)

Anda akan mendapatkan perkiraan yang salah mengenai kemungkinan bertahan hidup 1 tahun ketika Anda mengabaikan fakta bahwa 42 pasien mengalami sensor sebelum 1 tahun.

Ingatlah perkiraan yang benar dari kemungkinan 1 tahun untuk survival adalah 41%.

5 Memeriksa Mengabaikan Censoring

Bayangkan dua penelitian, masing-masing dengan 228 subjek. Ada 165 kematian di setiap penelitian. Tidak ada sensor di suatu garis (garis oranye), 63 pasien sensor di suatu garis lainnya (garis biru).

Mengabaikan censoring mengarah pada perkiraan yang terlalu tinggi dari probabilitas survival secara keseluruhan, karena subjek yang disensor hanya menyumbangkan informasi untuk sebagian dari waktu follow-up, dan kemudian keluar dari set risiko, sehingga menurunkan probabilitas kumulatif untuk bertahan hidup.

6 Memperkirakan Median Survival

Kuantitas lain yang sering digunakan dalam analisis survival adalah waktu bertahan hidup rata-rata, yang kita ukur menggunakan median. Waktu bertahan tidak diharapkan terdistribusi secara normal sehingga meannya bukan ringkasan yang tepat.

Kita dapat memperoleh ini langsung dari objek survfit

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

Kita melihat median waktu survival adalah 310 hari. Batas bawah dan atas dari interval kepercayaan 95% juga ditampilkan.

6.1 Waktu dan Kemungkinan

Median Survival adalah waktu yang sesuai dengan probabilitas survival 0,5:

6.2 Dampak Mengabaikan Sensor

Apa yang terjadi jika Anda menggunakan perkiraan “naive”? Rangkum waktu bertahan hidup rata-rata di antara 165 pasien yang meninggal.

##   median_surv
## 1         226

Anda mendapatkan perkiraan waktu survival rata-rata yang salah yaitu 226 hari ketika Anda mengabaikan fakta bahwa pasien yang disensor juga berkontribusi pada waktu follow-up.

Ingat perkiraan yang benar dari waktu survival rata-rata adalah 310 hari.

Mengabaikan penyensoran membuat kurva survival yang diturunkan secara artifisial karena waktu follow-up yang berkontribusi pada pasien yang disensor tidak termasuk (garis ungu).

Kurva survival untuk data lung ditampilkan dengan warna biru untuk perbandingan.

7 Membandingkan Waktu Survival

Kita dapat melakukan tes signifikansi antar-grup menggunakan uji log-rank.

Uji log-rank sama-sama membobotkan pengamatan selama seluruh waktu follow-up dan merupakan cara paling umum untuk membandingkan waktu bertahan hidup antar kelompok.

Ada versi yang lebih menekankan pada follow-up awal atau akhir yang mungkin lebih sesuai tergantung pada pertanyaan penelitian (lihat ?Survdiff untuk opsi pengujian yang berbeda).

Kita mendapatkan nilai p log-rank menggunakan fungsi survdiff. Misalnya, kita dapat menguji apakah ada perbedaan waktu bertahan hidup menurut jenis kelamin dalam 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 Mengekstrak informasi

Sebenarnya agak rumit untuk mengekstrak nilai-p dari hasil survdiff. Berikut adalah baris kode untuk melakukannya

## [1] 0.001311165

Atau ada fungsi sdp dalam package ezfun, yang dapat Anda instal menggunakan devtools::install_github("zabore/ezfun"). Ini menampilkan nilai p yang diformat

## [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 pengaruh 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)\text{exp}(\beta_1X_{i1}+⋯+\beta_pX{ip})\]

\(h(t)\) : hazard, atau kurs sesaat di mana peristiwa terjadi \(h_0(t)\) : bahaya dasar yang mendasarinya

Beberapa asumsi utama model:

  • sensor non-informatif
  • hazard proporsional

Catatan: model regresi parametrik untuk hasil survival juga tersedia, tetapi model tersebut tidak akan dibahas pada latihan ini.

Kita dapat menyesuaikan model regresi untuk data kelangsungan hidup 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 Memformat Regresi Cox

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

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

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

Atau 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

Jumlah minat dari model regresi Cox adalah rasio hazard (HR). HR mewakili rasio hazard antara dua kelompok pada titik waktu tertentu.

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

Jika Anda memiliki parameter regresi β (dari kolom perkiraan di coxph) maka HR = \(exp(β)\).

HR < 1 menunjukkan pengurangan hazard kematian sedangkan HR > 1 menunjukkan peningkatan hazard kematian.

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