Analisis Survival sesuai dengan serangkaian pendekatan statistik yang digunakan untuk menyelidiki waktu yang diperlukan untuk hal yang terjadi. Juga disebut Data Peristiwa yang terdiri dari waktu mulai sampai dengan waktu akhir yang berbeda.
Analisis kelangsungan hidup digunakan dalam berbagai bidang seperti:
Dalam Studi kanker, tipe pertanyaan penelitian seperti :
Apa Karakteristik dampak tertentu pada kelangsungan hidup pasien yang kemungkinan bertahan selama 3 tahun ? Apakah ada perbedaan dalam kelangsungan hidup antara kelompok pasien?
Waktu untuk data peristiwa pada umumnya terjadi di banyak bidang, tetapi tidak terbatas pada
Karena analisis kelangsungan hidup adalah umum digunakan di banyak bidang lain dengan nama lain :
Tujuan dari bab ini adalah untuk menggambarkan konsep dasar analisis kelangsungan hidup. Dalam studi kanker, sebagian besar analisis kelangsungan hidup menggunakan metode berikut:
Di sini, kita akan mulai dengan menjelaskan konsep penting analisis kelangsungan hidup, termasuk:
Kemudian, kita akan melanjutkan dengan menggambarkan analisis multivariat menggunakan model bahaya proporsional Cox.
untuk subjek \(i\) :
Waktu Observasi dan indikator peristiwa diberikan oleh data lung
Dataset lung tersedia dari package survival dalam R. Data tersebut mengandung subjek dengan kanker paru-paru stadium lanjut dari North Central Cancer Treatment Group. Beberapa variabel yang akan kami gunakan untuk menunjukkan metode ini termasuk :
## 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
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.
Subjek dapat disensor karena:
Secara khusus ini adalah contoh pensesoran kanan.
Penyensoran kiri dan sensor interval juga mungkin terjadi, dan metode yang ada untuk menganalisis jenis data ini, tetapi pelatihan ini akan terbatas pada sensor kanan.
Dalam contoh ini, bagaimana kita menghitung proporsi siapa yang bebas acara pada 10 tahun?
Subyek 6 dan 7 bebas acara pada usia 10 tahun. Subyek 2, 9, dan 10 memiliki acara sebelum 10 tahun. Subjek 1, 3, 4, 5, dan 8 disensor sebelum 10 tahun, jadi kita tidak tahu apakah mereka memiliki acara atau tidak oleh 10 tahun - bagaimana kita memasukkan subjek ini ke dalam perkiraan kita?
Data akan sering datang dengan tanggal mulai dan berakhir daripada waktu bertahan hidup yang telah dihitung sebelumnya. Langkah pertama adalah memastikan ini diformat sebagai tanggal dalam R.
Mari membuat contoh kecil dari dataset dengan variabel sx_date untuk tanggal operasi dan last_fup_date untuk tanggal tindak lanjut terakhir.
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.4 v stringr 1.4.0
## v tidyr 1.1.2 v forcats 0.5.0
## v readr 1.4.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
##
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice
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
Kami melihat ini adalah kedua karakter variabel yang 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
Kita juga dapat menggunakan paket lubridate untuk memformat tanggal. Dalam hal ini, gunakan fungsi ?dmy untuk opsi lebih lanjut.
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
## # 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
R, pemisah tidak perlu ditentukanSetelah tanggal diformat, kita perlu menghitung perbedaan antara waktu mulai dan berakhir di beberapa unit, biasanya bulan atau tahun. Di basis R, gunakan difftime untuk menghitung jumlah hari di antara dua tanggal dan mengonversinya menjadi nilai numerik menggunakan as.numerik. Kemudian konversi ke tahun dengan membagi sebesar 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 %--% menunjuk interval waktu, yang kemudian dikonversi ke jumlah detik yang berlalu menggunakan as.duration dan akhirnya dikonversi ke tahun dengan membagi dengan dyears(1), yang memberikan jumlah detik 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
Catatan: kita perlu memuat paket lubridate menggunakan panggilan ke perpustakaan agar dapat mengakses operator khusus (mirip dengan situasi dengan pipa)
Untuk komponen data kelangsungan hidup saya menyebutkan indikator peristiwa: Indikator Peristiwa \(\delta_i\):
Bagaimana pun, dalam R fungsi Surv akan menerima TRUE/FALSE (TRUE = peristiwa) atau 1/2 (2 = peristiwa).
dalam data lung, kami memiliki:
Peluang bahwa subjek akan bertahan di luar waktu yang ditentukan
\[S(t)=Pr(T>t)=1−F(t)\]
\(S(t)\) : fungsi kelangsungan hidup \(F(t)=Pr(T≤t)\): fungsi kumulatif distribusi
Secara teori fungsi kelangsungan hidup lancar; dalam praktiknya kami mengamati peristiwa dalam skala waktu diskrit.
Metode Kaplan-Meier adalah cara paling umum untuk memperkirakan waktu dan probabilitas kelangsungan hidup. Ini adalah pendekatan non-parametrik yang menghasilkan fungsi langkah, di mana ada langkah mundur setiap kali peristiwa terjadi.
Surv dari paket survival menciptakan objek bertahan hidup untuk digunakan sebagai respons dalam rumus model. Akan ada satu entri untuk setiap subjek yang merupakan waktu bertahan hidup, yang diikuti oleh + jika subjek disensor. Mari kita lihat 10 pengamatan pertama:## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166
survfit membuat kurva survival pada rumus, Mari kita hasilkan kurva kelangsungan hidup keseluruhan untuk seluruh kohort, tetapkan ke objek f1, dan lihat names objek itu:## [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 utama objek survfit ini yang akan digunakan untuk membuat kurva bertahan hidup meliputi:
time, yang berisi titik awal dan titik akhir dari setiap interval waktusurv, yang berisi probabilitas kelangsungan hidup yang sesuai dengan masing-masing timeSekarang kita plot objek survfit di R dasar untuk mendapatkan plot Kaplan-Meier.
plot(survfit(Surv(time, status) ~ 1, data = lung),
xlab = "Days",
ylab = "Overall survival probability")mark.time = TRUE)Atau, fungsi ggsurvplot dari paket survminer dibangun di atas ggplot2 , dan dapat digunakan untuk membuat plot Kaplan-Meier. Lihat lembar cheet untuk paket survminer.
## Loading required package: ggpubr
ggsurvplot(
fit = survfit(Surv(time, status) ~ 1, data = lung),
xlab = "Days",
ylab = "Overall survival probability")ggsurvplot menunjukkan fungsi langkah (garis solid) dengan pita kepercayaan terkait (area berbayang).censor= FALSESatu kuantitas yang menarik dalam analisis kelangsungan hidup adalah kemungkinan bertahan hidup di luar angka tertentu \((x)\) tahun.
sebagai contoh untuk mengestimasi kelangsungan hidup pada tahun ke 1 menggunakan summary dengan argument times (Perhatikan variabel waktu dalam data lung sebenarnya dalam beberapa 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
Kami menemukan bahwa probabilitas 1 tahun untuk bertahan hidup dalam penelitian ini adalah 41%.
Batas bawah dan atas terkait dari interval kepercayaan 95% juga ditampilkan.
Probabilitas kelangsungan hidup 1 tahun adalah titik pada sumbu y yang sesuai dengan 1 tahun pada sumbu x untuk kurva bertahan hidup.
Apa yang terjadi jika kamu menggunakan estimasi “Naive”?
121 dari 228 pasien meninggal 1 tahun jadi: \((1-\frac{121}{228})x100=47%\)
Kuantitas lain yang sering menarik dalam analisis kelangsungan hidup adalah waktu bertahan hidup rata-rata, yang kami kuantifikasi menggunakan median. Waktu bertahan hidup tidak diharapkan untuk didistribusikan secara normal sehingga rata-rata bukan ringkasan yang tepat.
Kita bisa mendapatkan 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
median Kelangsungan hidup adalah waktu yang sesuai dengan probabilitas kelangsungan hidup 0,5:
Apa yang terjadi jika Anda menggunakan perkiraan “naif”? Meringkas waktu kelangsungan hidup median di antara 165 pasien yang meninggal
## median_surv
## 1 226
lung ditunjukkan dalam warna biru untuk perbandingansurvdiff untuk opsi pengujian yang berbeda)Kami mendapatkan nilai p peringkat log menggunakan fungsi survdiff. Misalnya, kita dapat menguji apakah ada perbedaan waktu kelangsungan 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
Ini sebenarnya agak rumit untuk mengekstrak nilai p dari hasil survdiff. Berikut adalah baris kode untuk melakukannya
## [1] 0.001311165
terdapat fungsi sdp dalam paket ezfun , yang mana dapat diinstall menggunakan devtools::install_github("zabore/ezfun"). Mengembalikan pengformatan p-value
## [1] 0.001
Kami mungkin ingin mengukur ukuran efek untuk satu variabel, atau menyertakan lebih dari satu variabel ke dalam model regresi untuk memperhitungkan efek dari beberapa variabel.
Model regresi Cox adalah model semi-parametrik yang dapat digunakan untuk menyesuaikan model regresi yang tidak dapat divariasi dan multivariasi yang memiliki hasil bertahan hidup. \[h(t|X_i)=h_0(t)exp(β_1X_{i1}+⋯+β_{p}X_{ip})\]
\(h(t)\): bahaya, atau laju seketika di mana peristiwa terjadi
\(h_0(t)\): bahaya dasar yang mendasari
Beberapa asumsi utama model:
Catatan: model regresi parametrik untuk hasil bertahan hidup juga tersedia, tetapi 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.
## 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
Kita dapat melihat versi rapi dari output menggunakan fungsi tidy dari paket broom
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| sex | 0.5880028 | 0.1671786 | -3.176385 | 0.0014912 |
estimate di coxph) maka \(HR = exp(\beta)\).