Email: elvrianae@gmail.com
RPubs: https://rpubs.com/elvriana_e/

1 Apa itu Analisis Survival?

Analisis survival atau analisis kelangsungan hidup sesuai dengan kumpulan dari pendekatan statistik yang digunakan untuk menyelidiki waktu yang dibutuhkan agar sebuah peristiwa menarik terjadi. Itu disebut juga data Waktu ke peristiwa (time-to-event) yang terdiri dari waktu mulai dan waktu berakhir yang berbeda.

Analisis kelangsungan hidup digunakan dalam berbagai bidang seperti:

Studi kanker, untuk analisis waktu kelangsungan hidup pasien,
Sosiologi, untuk “analisis peristiwa-sejarah”,
Teknik, untuk “analisis kegagalan-waktu”.

Dalam studi kanker, pertanyaan penelitian yang khas seperti:

Apa dampak karakteristik klinis tertentu terhadap kelangsungan hidup pasien. Berapa probabilitas seseorang bertahan hidup selama 3 tahun? Adakah perbedaan dalam kelangsungan hidup antara kelompok pasien?

1.1 Contoh dari Bidang Lain

Data waktu ke peristiwa termasuk umum di berbagai bidang, tetapi tidak terbatas untuk:

Waktu dari infeksi HIV sampai berkembang menjadi AIDS
Waktu serangan jantung
Waktu mulai penyalahgunaan zat
Waktu inisiasi aktivitas seksual
Waktu kerusakan mesin

1.2 Nama Lain untuk Analisis Survival

Karena analisis survival umum di banyak bidang lain, itu juga menggunakan nama-nama lain, seperti:

Analisis reliabilitas
Analisis durasi
Analisis riwayat peristiwa
Analisis waktu ke peristiwa
Model Kelangsungan Hidup (Survival Model)

1.3 Tujuan Analisis Survival

Tujuan dari bab ini adalah untuk mendeskripsikan konsep dasar dari analisis kelangsungan hidup. Dalam studi kanker, sebagian besar analisis kelangsungan hidup menggunakan metode-metode seperti berikut:

Plot Kaplan-Meier untuk memvisualisasikan kurva kelangsungan hidup
Tes log-rank untuk membandingkan kurva kelangsungan hidup dari dua atau lebih kelompok
Regresi risiko proporsional Cox untuk menggambarkan pengaruh variabel terhadap kelangsungan hidup. Model Cox dibahas dalam bab berikutnya: Model risiko proporsional Cox.

Di sini, kita akan mulai dengan menjelaskan konsep penting dari analisis survival (kelangsungan hidup), termasuk:

Bagaimana caranya untuk menghasilkan dan menafsirkan kurva kelangsungan hidup, dan
Bagaimana caranya untuk mengukur dan menguji perbedaan kelangsungan hidup antara dua atau lebih kelompok pasien.

Kemudian, kita akan melanjutkan dengan mendeskripsikan analisis multivariat menggunakan model risiko proporsional Cox.

2 Komponen dari Survival

Untuk subjek \(i\) :

Waktu peristiwa \(T_i\)
Menyensor waktu \(C_i\)
Indikator peristiwa \(\delta_i\)
- 1 jika peristiwa diamati (yaitu \(T_i \leq C_i\))
- 0 jika disensor (yaitu \(T_i \geq C_i\))
Waktu yang diamati \(Y_i = \min(T_i, C_i)\)

Waktu yang diamati dan indikator peristiwa disediakan dalam data lung

Waktu : Waktu bertahan hidup dalam beberapa hari (\(Y_i\))
Status : Status penyensoran, 1 = disensor, 2 = mati (\(\delta_i\))

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

Waktu : Waktu bertahan hidup dalam beberapa hari
Status : menyensor status 1 = disensor, 2 = mati
Jenis Kelamin : Pria = 1, Wanita = 2

##   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 Penyensoran?

Referensinya adalah: 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 Penyensoran

Sebuah subjek dapat disensor karena:

Tidak ditindak lanjut
Penarikan dari penelitian
Tidak ada peristiwa hingga akhir masa penelitian

Secara khusus ini adalah contoh sensor kanan.

Sensor kiri dan sensor interval juga mungkin terjadi, dan ada juga metode untuk menganalisis jenis data ini, tetapi akan dibatasi pada sensor kanan.

3.2 Data Survival yang Disensor

Dalam contoh ini, bagaimana kita menghitung proporsi yang bebas peristiwa pada 10 tahun?

Subjek 6 dan 7 bebas acara pada 10 tahun. Subjek 2, 9, dan 10 memiliki peristiwa tersebut sebelum 10 tahun. Subjek 1, 3, 4, 5, dan 8 disensor sebelum 10 tahun, jadi kita tidak tahu apakah mereka memiliki peristiwa tersebut atau tidak dalam 10 tahun - bagaimana kita memasukkan subjek ini ke dalam perkiraan kita?

3.3 Distribusi Tindak Lanjut

Subjek yang disensor masih memberikan informasi sehingga harus dimasukkan secara tepat dalam analisis
Distribusi waktu tindak lanjut tidak tepat, dan mungkin berbeda antara pasien yang disensor dan mereka yang mengalami peristiwa
Waktu tindak lanjut selalu positif

4 Berurusan dengan Tanggal di R

Data sering datang dengan tanggal mulai dan akhir daripada waktu kelangsungan hidup yang telah dihitung sebelumnya. Langkah pertama adalah memastikan ini diformat sebagai tanggal (as dates) di R.

Mari kita buat contoh dataset kecil dengan variabel sx_date untuk tanggal operasi dan last_fup_date untuk tanggal tindak lanjut terakhir.

library(tibble)
date_ex <- tibble(
  sx_date = c("2007-06-22", "2004-02-13", "2010-10-27"),
  last_fup_date = c("2017-04-15", "2018-07-04", "2016-10-31")
    )

date_ex

## # 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 melihat ini adalah kedua variabel karakter, yang akan sering terjadi, tetapi kita membutuhkannya untuk diformat sebagai tanggal. Cara lain untuk memformat tanggal menggunakan R:

library(dplyr)
date_ex %>% 
  mutate(
    sx_date = as.Date(sx_date, format = "%Y-%m-%d"), 
    last_fup_date = as.Date(last_fup_date, format = "%Y-%m-%d") 
    )

## # 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 di R, formatnya harus menyertakan pemisah serta simbol. misalnya jika tanggal Anda dalam format m/d/Y maka Anda memerlukan format = "%m/%d/%Y"
Lihat daftar lengkap simbol format tanggal di Link Ini

Kita juga dapat menggunakan paket lubridate untuk memformat tanggal. Dalam kasus ini, kita menggunakan fungsi ymd

library(lubridate)
date_ex %>% 
  mutate(
    sx_date = ymd(sx_date), 
    last_fup_date = ymd(last_fup_date)
    )

## # 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, pemisahnya tidak perlu ditentukan
Halaman bantuan untuk ?dmy akan menampilkan semua opsi format.

4.1 Menghitung Waktu Survival

Sekarang setelah tanggal diformat, kita perlu menghitung perbedaan antara waktu mulai dan berakhir dalam beberapa unit, biasanya bulan atau tahun.

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

date_ex %>% 
  mutate(
    os_yrs = 
      as.numeric(
        difftime(last_fup_date, 
                 sx_date, 
                 units = "days")) / 365.25
    )

## # 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

Dengan menggunakan paket lubridate, operator %--% menetapkan interval waktu, yang kemudian diubah menjadi jumlah detik yang telah berlalu menggunakan fungsi as.duration dan akhirnya diubah menjadi tahun dengan membaginya dengan dyears(1), yang memberikan jumlah detik dalam tahun.

date_ex %>% 
  mutate(
    os_yrs = 
      as.duration(sx_date %--% last_fup_date) / dyears(1)
    )

## # 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 memuat paket lubridate menggunakan panggilan ke perpustakaan agar dapat mengakses operator khusus (mirip dengan situasi dengan pipes)

4.2 Indikator Peristiwa

Untuk komponen data survival, disini saya menyebutkan indikator peristiwa:

Indikator peristiwa \(\delta_i\):

1 jika peristiwa diamati (yaitu \(T_i \leq C_i\))
0 jika disensor (yaitu \(T_i > C_i\))

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

Dalam data lung, kita memiliki:

status: status penyensoran 1 = disensor, 2 = mati

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 \leq t)\): fungsi distribusi kumulatif

Secara teori, fungsi bertahan hidup itu lancar, dalam praktiknya kita mengamati peristiwa pada skala waktu yang berbeda.

4.4 Probabilitas Survival

Probabilitas Survival pada waktu tertentu, \(S(t)\), adalah probabilitas bersyarat untuk bertahan hidup setelah waktu itu, mengingat bahwa seseorang telah bertahan sebelum waktu itu.
Dapat diperkirakan sebagai jumlah pasien yang hidup tanpa tindak lanjut pada saat itu, dibagi dengan jumlah pasien yang hidup sebelum waktu itu
Estimasi Kaplan-Meier dari probabilitas survival adalah produk dari probabilitas bersyarat hingga saat itu
Pada waktu 0, probabilitas survival adalah 1, yaitu \(S(t_0) = 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 paket 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 disensor. Mari kita lihat 10 observasi pertama:

library(survival)
Surv(lung$time, lung$status)[1:10]

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

4.6 Memperkirakan Kurva Survival

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

f1 <- survfit(Surv(time, status) ~ 1, data = lung)
names(f1)

##  [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 ini yang akan digunakan untuk membuat kurva survival antara lain:

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

4.7 Plot Survival Kaplan-Meier

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

plot(survfit(Surv(time, status) ~ 1, data = lung), 
     xlab = "Days", 
     ylab = "Overall survival probability")

Plot default di R menunjukkan fungsi langkah (garis penuh) dengan interval kepercayaan terkait (garis putus-putus)
Garis horizontal mewakili durasi survival untuk interval tersebut
Interval diakhiri oleh suatu peristiwa
Ketinggian garis vertikal menunjukkan perubahan probabilitas kumulatif
Pengamatan yang disensor, ditunjukkan dengan tanda centang, mengurangi kelangsungan hidup 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 paket survminer dibangun di atas ggplot2, dan dapat digunakan untuk membuat plot Kaplan-Meier. Lihat lembar cheatsheet pada Link Ini untuk paket survminer.

library(survminer)
ggsurvplot(
    fit = survfit(Surv(time, status) ~ 1, data = lung), 
    xlab = "Days", 
    ylab = "Overall survival probability")

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

4.8 Memperkirakan 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 hidup hingga \(1\) tahun, gunakan summary dengan argumen times (Perhatikan variabel times di data lung sebenarnya dalam hari, jadi kita perlu menggunakan times = 365,25)

summary(survfit(Surv(time, status) ~ 1, data = lung), 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 menemukan 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 \(x\)-tahun Probabilitas Survival

Probabilitas kelangsungan hidup \(1\) tahun adalah titik pada sumbu \(y\) yang sesuai dengan \(1\) tahun pada sumbu \(x\) untuk kurva kelangsungan hidup.

Apa yang terjadi jika Anda menggunakan perkiraan “naive”?

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

Anda mendapatkan perkiraan yang salah tentang kemungkinan bertahan hidup \(1\) tahun ketika Anda mengabaikan fakta bahwa 42 pasien disensor sebelum \(1\) tahun.
- Ingat perkiraan yang benar dari kemungkinan 1 tahun untuk bertahan hidup adalah 41%.

5 Periksa Mengabaikan Sensor

Bayangkan dua studi, masing-masing dengan 228 subjek. Ada 165 kematian di setiap penelitian. Tidak ada sensor di satu (garis oranye), 63 pasien disensor di satu (garis biru)
Mengabaikan penyensoran mengarah pada perkiraan yang terlalu tinggi dari probabilitas kelangsungan hidup secara keseluruhan, karena subjek yang disensor hanya menyumbangkan informasi untuk sebagian dari waktu tindak lanjut, dan kemudian keluar dari set risiko, sehingga menurunkan probabilitas kumulatif untuk bertahan hidup

6 Memperkirakan Survival Median

Kuantitas lain yang sering menarik dalam analisis kelangsungan hidup 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 kita

survfit(Surv(time, status) ~ 1, data = lung)

## 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 bertahan hidup adalah 310 hari Batas bawah dan atas dari interval kepercayaan 95% juga ditampilkan.

6.1 Waktu dan Probabilitas

Kelangsungan hidup median adalah waktu yang sesuai dengan probabilitas kelangsungan hidup 0,5:

6.2 Dampak dari Mengabaikan Sensor

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

lung %>% 
  filter(status == 2) %>% 
  summarize(median_surv = median(time))

##   median_surv
## 1         226

Anda mendapatkan perkiraan waktu kelangsungan hidup rata-rata yang salah yaitu 226 hari ketika Anda mengabaikan fakta bahwa pasien yang disensor juga berkontribusi pada waktu tindak lanjut.
Ingat perkiraan yang benar dari waktu bertahan hidup rata-rata adalah 310 hari.
Mengabaikan penyensoran menciptakan kurva kelangsungan hidup yang diturunkan secara artifisial karena waktu tindak lanjut yang berkontribusi pada pasien yang disensor tidak termasuk (garis ungu)
Kurva kelangsungan hidup sebenarnya untuk data lung ditampilkan dengan warna biru untuk perbandingan

7 Membandingkan Waktu Survival

Kita dapat melakukan uji signifikansi antarkelompok menggunakan tes log-rank
Tes log-rank sama-sama membobotkan pengamatan selama seluruh waktu tindak lanjut dan merupakan cara paling umum untuk membandingkan waktu bertahan hidup antar kelompok
Ada versi yang lebih membebani tindak lanjut awal atau akhir yang bisa lebih sesuai tergantung pada pertanyaan penelitian (lihat ?Survdiff untuk opsi tes yang berbeda)

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

survdiff(Surv(time, status) ~ sex, 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 cara melakukannya

sd <- survdiff(Surv(time, status) ~ sex, data = lung)
1 - pchisq(sd$chisq, length(sd$n) - 1)

## [1] 0.001311165

Atau ada fungsi sdp di paket ezfun, yang dapat Anda instal menggunakan devtools::install_github(“zabore/ezfun”). Ini mengembalikan nilai p yang sudah diformat

ezfun::sdp(sd)

## [1] 0.001

9 Model Regresi Cox

Kita mungkin ingin mengukur ukuran efek untuk satu variabel, atau menyertakan lebih dari satu variabel ke dalam model regresi untuk memperhitungkan efek 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)\) :risiko, atau kecepatan sesaat di mana peristiwa terjadi
\(h_0(t)\):risiko dasar yang mendasari

Beberapa asumsi utama model:

sensor non-informatif
risiko proporsional

Catatan: model regresi parametrik untuk hasil kelangsungan hidup juga tersedia, tetapi model tersebut tidak akan dibahas dalam pelatihan 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.

coxph(Surv(time, status) ~ sex, data = lung)

## 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 keluaran yang rapi menggunakan fungsi tidy dari paket broom:

library(kableExtra)
broom::tidy(
  coxph(Surv(time, status) ~ sex, data = lung), 
  exp = TRUE
  ) %>% 
  kable()

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

Atau gunakan tbl_regression dari paket gtsummary

coxph(Surv(time, status) ~ sex, data = lung) %>% 
  gtsummary::tbl_regression(exp = TRUE)

Characteristic	HR¹	95% CI¹	p-value
sex	0.59	0.42, 0.82	0.001
¹ HR = Hazard Ratio, CI = Confidence Interval

11 Rasio Risiko (Hazard Ratio)

Kuantitas yang menarik dari model regresi Cox adalah rasio risiko atau Hazard Ratio (HR). HR mewakili rasio bahaya 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 ditafsirkan seperti itu.
Jika Anda memiliki parameter regresi \(\beta\) (dari estimasi kolom di coxph kita) maka \(HR = \exp(\beta)\).
\(HR <1\) menunjukkan pengurangan bahaya kematian sedangkan \(HR> 1\) menunjukkan peningkatan bahaya kematian.
Jadi HR = 0,59 kita menyiratkan bahwa sekitar \(0,6\) kali lebih banyak wanita yang meninggal daripada pria, pada waktu tertentu.

Pengenalan Analisis Survival

Konsep Dasar dan Analisis Awal

Elvriana Elvani

February 16, 2021