Github: https://github.com/sofia3484

1 Apa itu Analisis Survival?

Analisis survival berkoresponden dengan seperangkat pendekatan statistik yang digunakan untuk menyelidiki suatu waktu yang dibutuhkan untuk suatu peristiwa yang menarik terjadi. Ini disebut juga sebagai waktu ke peristiwa (time-to-event) yang terdiri dari waktu mulai dan waktu berakhir yang pastinya berbeda.

Analisis Survaival digunakan dalam berbagai bidang seperti:

Studi penyakit kanker untuk menganalisis waktu kelangsungan hidup atau waktu survival pasien;
Studi sosiologi untukk “analisis peristiwa sejarah”;
Di bidang teknik untuk mempelajari “analisis waktu kegagalan” atau " failure-time-analysis";

Pada studi kanker, terdapat beberapa pertanyaan yang khas seperti:

Apa dampak karakteristik klinis tertentu terhadap kelangsungan hidup atau survival hidup pasien. Berapa probabilitas seseorang bertahan hidup dalam tiga tahun? Adakah perbedaan dalam kelangsungan hidup atau survival antar kelompok pasien?

1.1 Contoh dari Bidang Lain

Data waktu ke peristiwa (time-to-event) sangatlah umum di banyak bidang termasuk, tetapi tidak terbatas pada

Waktu untuk infeksi HIV sampai berkembang menjadi AIDS;
Waktu untuk serangan jantung;
Waktu untuk mulainya penyalahgunaan zat;
Waktu untuk inisiasi aktivitas seksual;
Waktu untuk mesin rusak

1.2 Alias untuk Analisis survival

Karena analiss survival atau analisis kelangsungan hidup sangatlah umum di banyak bidang, analisis ini juga dinamakan sebagai:

Analisis Reabilitas;
Analisis Durasi;
analisis peristiwa lampau (Event History Analysis);
Analisis Waktu ke Peristiwa (Time-to-Event Analysis);
Model Bertahan Hidup

1.3 1.3 Tujuan Analisis Kelangsungan Hidup atau 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 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, kami akan mulai dengan menjelaskan konsep penting dari analisis kelangsungan hidup, termasuk:

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

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

2 Komponen dari Kelangsungan Hidup

untuk subjek $i$:

1. Waktu Perisiwa $T_i$
1. Waktu Sensoring $C_i$
1. Indikator Peristiwa: $δ_i$:
- 1 jika peristiwa diamati (yakni $T_i$ ≤ $C_i$)
- 0 jika tersensor (yakni T_i$ > $C_i$)
1. Waktu Observasi $Y_i = min(T_i, C_i$)

Waktu yang diamati dan indikator kejadian disediakan dalam data lung

time: Waktu bertahan dalam hari ($Y_i$)
status: status sensor 1 = tersensor, 2 = meninggal ($δ_i$)

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

time: Waktu bertahan dalam hari
status: menyensor status 1 = tersensor, 2 = meninggal
sex: 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 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

Sebuah subjek dapat tersensor karena:

Tidak lanjut
Penarikan dari studi
Tidak ada peristiwa hingga akhir masa studi tetap

Secara khusus ini adalah contoh sensor kanan.

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

3.2 Sensoring pada Data Survival

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

Subjek 6 dan 7 adalah bebas peristiwa pada 10 tahun. Subjek 2, 9, dan 10 mengalami kejadian tersebut sebelum 10 tahun. Subjek 1, 3, 4, 5, dan 8 tersensor sebelum tahun ke-10, jadi kami tidak tahu apakah mereka mengalami peristiwa tersebut atau tidak hingga tahun ke-10 - bagaimana kami memasukkan subjek ini ke dalam perkiraan kami?

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 kejadiannya
Waktu tindak lanjut (follow-up) selalu positif

4 Berurusan dengan Tanggal di R

Data sering kali dimulai dan diakhiri dengan tanggal daripada menggunakan waktu kelangsungan hidup yang telah dihitung sebelumnya. Langkah pertama adalah memastikan ini diformat sebagai tanggal di R.

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

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

Kami melihat ini adalah kedua variabel karakter, yang akan sering terjadi, tetapi kami membutuhkannya untuk diformat sebagai tanggal. Cara lain untuk memformat tanggal menggunakan R:

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 basis format R 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 https://www.statmethods.net/input/dates.html

Kami juga dapat menggunakan paket lubridate untuk memformat tanggal. Dalam kasus ini, gunakan fungsi ymd

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

4.1 Mengkalkulasikan Survival Times atau Waktu Kelangsungan Hidup

Sekarang setelah tanggal diformat, kita perlu menghitung perbedaan antara waktu mulai dan akhir dalam beberapa unit, biasanya dalam bulan atau tahun. Dalam basis R, gunakan difftime untuk menghitung jumlah hari antara dua tanggal kita dan mengubahnya menjadi nilai numerik menggunakan 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 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 menjalankan paket lubridate menggunakan panggilan ke library agar dapat mengakses operator khusus (mirip dengan situasi dengan pipes)

4.2 Indikator Peristiwa

Untuk komponen data survival yangsaya sebutkan, indikator peristiwanya:

Indikator Peristiwa: $δ_i$:

* 1 jika peristiwa diamati (yakni $T_i$ \ ≤ \ $C_i$)

* 0 jika tersensor (yakni T_i$ \ > \ $C_i$)

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

Dalam data lung, kami memiliki:

status: menyensor status 1 = disensor, 2 = mati

4.3 Fungsi Kelangsungan Hidup atau Fungsi Survival

Probabilitas bahwa subjek akan bertahan melebihi waktu tertentu

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

$S (t)$: fungsi kelangsungan hidup $F (t) = Pr (T≤t)$: fungsi distribusi kumulatif

Secara teori, fungsi bertahan hidup itu mulus; dalam praktiknya kami mengamati peristiwa pada skala waktu yang berbeda.

4.4 Probabilitas Kelangsungan Hidup

Probabilitas kelangsungan hidup pada waktu tertentu, S (t), adalah probabilitas bersyarat untuk bertahan hidup setelah waktu tersebut, mengingat bahwa seseorang telah bertahan sebelum waktu itu.
Dapat diperkirakan sebagai jumlah pasien yang hidup tanpa mangkir pada saat itu, dibagi dengan jumlah pasien yang hidup sebelum waktu itu
Estimasi Kaplan-Meier dari probabilitas kelangsungan hidup adalah produk dari probabilitas bersyarat hingga saat itu
Pada waktu 0, probabilitas kelangsungan hidup 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 paket survival adalah 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:

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

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

4.6 Memperkirakan Kurva Bertahan Hidup

Fungsi survfit adalah menciptakan kurva kelangsungan hidup berdasarkan rumus. Mari buat kurva kelangsungan hidup keseluruhan untuk seluruh kelompok, tetapkan ke objek f1, dan lihat names dari objek-objek tersebut:

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 meliputi:

waktu, yang berisi titik awal dan akhir setiap interval waktu
surv, yang berisi probabilitas kelangsungan hidup yang berkorespondensi dengan setiap waktu

4.7 Plot Kelangsungan Hidup Kaplan-Meier

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

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

Plot default di basis R menunjukkan fungsi step (garis penuh) dengan interval kepercayaan terkait (garis putus-putus)
Garis horizontal mewakili durasi kelangsungan hidup untuk interval tersebut
Interval diakhiri oleh suatu peristiwa
Ketinggian garis vertikal menunjukkan perubahan probabilitas kumulatif
Pengamatan sensor, 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 contekan untuk paket 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, dalam contoh ini garis tersebut agak dikaburkan , dan dapat ditekan menggunakan opsi censor = FALSE

4.8 Memperkirakan Kelangsungan Hidup 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 time 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

Kami 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 Kemungkinan Bertahan Hidup

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 “naif”?

121 dari 228 pasien meninggal dalam $1$ tahun jadi:

\[(1− \frac{121}{228})×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 Cek 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 Median Survival

Kuantitas lain yang sering menarik dalam analisis kelangsungan hidup adalah rata-rata waktu bertahan hidup, yang kami ukur menggunakan median. Waktu bertahan tidak diharapkan terdistribusi secara normal sehingga rata-ratanya bukan ringkasan yang tepat.

Kami dapat memperoleh ini langsung dari objek survfit

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

Kami melihat median waktu bertahan hidup adalah 310 hari. Batas bawah dan atas dari interval kepercayaan 95% juga ditampilkan.

6.1 Waktu dan Kemungkinan

Median Kelangsungan hidup adalah waktu yang sesuai dengan probabilitas kelangsungan hidup $0.5$:

knitr::include_graphics("C:/Users/rony/Downloads/SM8.png")

6.2 Dampak Mengabaikan Sensor

Apa yang terjadi jika anda menggunakan perkiraan “naif”? 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 Bertahan Hidup

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

Kami mendapatkan nilai p log-rank menggunakan fungsi survdiff. Misalnya, kami 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 merepotkan untuk mengekstrak nilai-p dari hasil survdiff. Berikut di bawah in adalah

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 diinstal menggunakan devtools::install_github("zabore/ezfun"). Ini mengembalikan format nilai p

ezfun::sdp(sd)

## [1] 0.001

9 Model Regresi Cox

Kami 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(β_1X_{i1}+⋯+β_pX_{ip})\]

$h(t)$: bahaya (hazard), atau kecepatan sesaat di mana peristiwa terjadi $h0(t)$: bahaya dasar yang mendasari

Beberapa asumsi utama model:

sensor non-informatif
bahaya 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 lebih 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 gtsummarry

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 Bahaya

Kuantitas yang menarik dari model regresi Cox adalah rasio bahaya (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 dalam peristiwa tersebut. Ini bukan risiko, meskipun secara umum ditafsirkan seperti itu.
- Jika anda memiliki parameter regresi $β$ (dari kolom estimate di coxph kami) maka $HR = exp (β)$.
- HR < 1 menunjukkan pengurangan bahaya kematian sedangkan HR > 1 menunjukkan peningkatan bahaya kematian.
- Jadi HR = 0,59 menyiratkan bahwa sekitar 0,6 kali lebih banyak wanita yang meninggal daripada pria, pada waktu tertentu.

Pengenalan Analis Survival

Konsep Dasar dan Analisis Pertama

Sofia Nicklaus S.

February 07, 2021