Email: theodoraputrina@gmail.com
RPubs: https://rpubs.com/Theodora

1 Apa itu Analisis Survival?

1.1 Contoh Dari Bidang Lain

Data waktu ke acara umum yang termasuk pada banyak bidang, tetapi tidak terbatas pada:

Waktu dari Infeksi HIV hingga perkembangan AIDS
Waktu hingga serangan jantung
Waktu hingga penyalahgunaan zat
Waktu untuk memulai aktivitas seksual
Waktu hingga kerusakan mesin

1.2 Sebutan lain Analisis Survival

Karena Analisis Survival umum pada banyak bidang lainnya, maka ada juga beberapa sebutan lainnya

Analisis Keandalan
Analisis Durasi
Analisis Riwayat Peristiwa
Waktu untuk Analisis Acara
Model SUrvival

1.3 Tujuan Model Bertahan Hidup

tujuan dari bab ini adalah untuk mendeskripsikan konsep dasar dari analisis survival. dalam studi kanker, sebagian besar analisis survival menggunakan metode berikut:

Plot Kaplan - Meler untuk memvisualisasikan kurva survival
Tes log-rank untuk membandingkan kurva survival dari dua atau lebih kelompok
regresi bahay proporsional COx untuk menggambarkan pengaruh variabel pada survival. Model Cox dibahas dalam bab berikutnya/: Model bahay proprorsional COx.

Disini, kita akan mulai dengan menjelaskan konsep esensial analisis kelangsungan hidup, termasuk:

Bagaimana membuat dan menafsirkan kurva survival, dan
Bagaiman mengukur dan menguji pebedaan kelangsungan hidup antara dua atau lebih kelompok pasien.

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

2 Komponen Survival

Untuk subjek \(i\)

Waktu acara \(T_i\)
Menyensor Waktu \(C_i\)
Indikator acara \(δ_i\) * 1 jika peristiwa diamati (i.e \(T_i≤C_i\)) * 0 jika disensor (i.e \(T_i>C_i\))
Waktu obesitas \(Y_i=min(T_i,C_i)\)

Waktu yang diamati dan indikator kejadian disediakan dalam data lung

waktu: waktu survival pada hari \((Y_i)\)
Status: menyesnor status 1=censored, 2=dead \((δ_i)\)

dataset lung tersedia dari survival paket pada R. Data berisis 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

Waktu: waktu bertahan hidup dalam hari
status: menyensor status 1=censored, 2=dead
jenis kelamin: Male=1, Female=2

3 Apa itu Menyensor?

Refernsinya 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 Tipe-tipe Menyensor

Suatu subjek dapat disensor karena:

Mangkir
Penarikan dari studi
Tidak ada acara hingga akhir periode studi tetap

Secara khusu ini adalah contoh rights sensor.

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

3.2 Data Kelangsungan Hidup yang Disensor

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

Subjek 6 dan 7 bebas acara pada 10 tahun. subjek 2, 9, dan 10 memiliki kejadian tersebut sebelum 10 tahun. Subjek 1, 3, 4, 5, dan 8 *disensor sebelum 10 tahun, jadi kita tidak tahu apakah mereka memiliki acara tersebut atau tidak hingga 10 tahun - bagaiman kita bisa 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 miring, dan mungkin berbeda anatara pasien yang disensor dan mereka yang memiliki kejadian waktu
Tindak lanjut positif

#Berurusan dengan tanggal di R{.tabset .tabset-fade .tabset-pills}

Data sering kali datang dengan tanggal awal dan akhir daripada waktu kelangsungan hidup yang telah dihitung sebelumnya. Langkah pertaa adalah memastikan ini diformat sebagi tanggal di R.

Mari kita buat contoh kecil dataset dengan variabel sx_date untuk tanggal opersi dan last_fup_date untuk tanggal tindak lanjur 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 di 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 basis R format harus meneyrtakan pemisah serta simbol. mislanya jika tanggal anda dalam format m/d/Y maka anda akan membutuhkan format="%m/%d/%Y"
Lihat daftar lengkap simbol format tanggal di https://www.statmethods.net/input/dates.html

Kita juga bisa menggunakan lubdridate untuk memformat tanggal. Dalam kasus ini, gunakan 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, epmisah tidka perlu ditentukan
Halaman bantuan untuk ?Dmy akan menampilkan semua ospi format

3.4 Menghitung waktu bertahan hidup

Sekarang setelah tanggal diformat, kita perlu menghitung perbedaan antara waktu mulai dan berkahir dalam beberpa unit, biasanya 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 ubahke tahun dengan membaginya dengan 365,25, jumlah rata-rat 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 menajdi jumlah detik yang telah berlalu menggunakan as.duration dan akhirnya diubah menjadi tahun dengan membagi 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 lubridate menggunakan panggilan ke library agar dapat mengakses operator khusus (mirip dengan situasi dengan pipa)

3.5 Indikator Kejadian

Untuk komponen data survival saya sebutkan indikator kejadian:

Indikator kejadian \(δ_i\): * 1 jika kejadian diamati (yaitu \(T_i≤C_i\)) * 0 jika disnsor (yaitu \(T_i>C_i\))

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

Pada data lung, kita tau:

status: menyensor status 1=censored, 2=dead

3.6 Fungsi bertahan hidup

Probabilitas bahwa subjek akan bertahan melebihi waktu tertentu

\[S(t) = Pr(T>t) = 1 - F(t)\] \(S(t)\): fungsi bertahan hidup \(F(t)=Pr(T≤t)\): distribusi kumulatif

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

3.7 Probabilitas kelangsungan hidup

Probabilitas kelangsungan hidup pada waktu tertentu, \(S(t)\) adalah probabilitas bersyarat untuk berthana 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 prbabilita skelangsungan hidup adalah produk dari probabilitas bersyarat hingga saat itu
Pada waktu 0, probabilitas kelangsungan hidup adalah 1, ayitu \(S(t_0)\)=1

3.8 Membuat objek bertahan hidup

Metode Kaplan-Meier adalah cara paling umum untuk memperkirakan waktu kelangsungan hidup dan probabilitas. Ini adalah pendekatan non-parametik yang menghasilkan fungsi langkah, dimana ada penuruan 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

3.9 Memperkirakan kurva kelangsungan hidup

Fungsi survfit menciptakan kurva kelangsungan hidup berdasarkan rumus. Mari buat kurva kelangsungan hidup keseluruhan untuk seluruh kelompok, tetapkan ke objek f1, dan lihat name 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 meliputi:

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

3.10 Kaplan-Meier plot bertahan hidup

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 langkah (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 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 catatan 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

3.11 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

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.

3.12 \(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 “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%.

4 Periksa yang menghindari penyensoran

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

5 Memperkirkaan median kelangsungan hidup

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

Kita dapat memperoleh ini langsung dari objek survfit kami

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.

5.1 Waktu dan Probabilitas

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

5.2 Dampak dari Mengabaikan Penyensran

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 paru-paru ditampilkan dengan warna biru untuk perbandingan

6 Membandingkan model kelangsungan hidup

Kita dapat melakukan uji signifikansi antarkelompok menggunakan uji 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, 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

7 Mengekstraksi Infromasi

Sebenarnya agak merepotkan untuk mengekstrak nilai-p dari hasil survdiff. Berikut adalah baris kode untuk 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 diformat

ezfun::sdp(sd)

## [1] 0.001

8 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)\): hazard, atau kecepatan sesaat di mana peristiwa terjadi \(h_0(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

9 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

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

10 Rasio Hazard

Kuantitas yang menarik dari model regresi Cox adalah rasio hazard (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 \(β\) (dari estimasi kolom di coxph kita) maka HR = \(exp (β)\).
HR <1 menunjukkan pengurangan bahaya kematian sedangkan HR> 1 menunjukkan peningkatan bahaya kematian.
Jadi HR kita = 0,59 menyiratkan bahwa sekitar 0,6 kali lebih banyak wanita yang meninggal daripada pria, pada waktu tertentu.

Survival Model

Terjemahan

Theodora Putrina Gea

February 14, 2021