Survival Model
UTS
| Kontak | : \(\downarrow\) |
| yosia.yosia@student.matanauniversity.ac.id | |
| yyosia | |
| RPubs | https://rpubs.com/yosia/ |
Analisis Survival 1: pengantar singkat ke dalam Kurva Kaplan-Meier
Analisis waktu bertahan hidup diperlukan dalam penelitian apa pun yang menyelidiki waktu untuk hasil tertentu yang menarik. Studi kanker dalam bidang kedokteran dan kegagalan pertama mobil dalam bidang teknik (analisis waktu kegagalan) adalah contoh yang baik. Hasil yang diinginkan dapat berupa kematian, remisi untuk kambuh, perkembangan, atau kegagalan. Titik waktu untuk mencapai hasil tersebut umumnya disebut peristiwa. Syukurlah, tidak semua “event” berakibat fatal 😃, tetapi kadang-kadang bahkan bisa menjadi hasil yang menguntungkan seperti keluar dari rumah sakit. Dan dengan demikian, analisis ketahanan hidup juga merupakan istilah umum, karena ini bukan hanya tentang kelangsungan hidup.
Mungkinkah ada sesuatu yang lebih mengerikan, selain kecelakaan kapal Titanic dengan banyak korban jiwa? Sayangnya, ya! Yakni, bukan kematian yang langsung. Bayangkan orang-orang yang selamat di tengah lautan gelap yang dingin, sebagian menggunakan jaket pelampung, sebagian lagi tidak, berpegangan pada sisa-sisa kapal, dengan kepanikan yang luar biasa dan sedikit harapan untuk diselamatkan. Tidak akan ada helikopter yang bisa menyelamatkan mereka (ini tahun 1912). Mungkin tidak akan ada satu kapal pun di perairan itu selama berbulan-bulan. Berapa lama mereka akan bertahan? Apakah waktu hanya akan memperpanjang penderitaan mereka atau meningkatkan kemungkinan bertahan hidup, karena pada akhirnya mereka akan ditemukan oleh beberapa kapal secara acak? Apakah wanita memiliki peluang lebih tinggi untuk bertahan hidup dibandingkan pria? Bagaimana dengan penumpang kaya vs penumpang biasa? Nah, dalam artikel ini kita akan menemukan jawaban untuk semua pertanyaan ini dan sambil belajar bagaimana analisis kelangsungan hidup bekerja. Secara khusus, kita akan:
- belajar tentang konsep-konsep terpenting dalam analisis survival: kurva survival, penyensoran, dan uji log-rank,
- menghitung secara manual, kemudian mengkomputasi dan menginterpretasikan kurva survival, dan
- menguji perbedaan kelangsungan hidup antara dua atau lebih kelompok, misalnya perempuan vs laki-laki
Mengapa kita membutuhkan survival analysis?
library(tidyverse)
d <- tibble(
time = c(1,1,1,2,2,3,4,5,7,10),
status = c(1,1,1,1,0,1,0,0,1,0)
)
d %>% kable()| time | status |
|---|---|
| 1 | 1 |
| 1 | 1 |
| 1 | 1 |
| 2 | 1 |
| 2 | 0 |
| 3 | 1 |
| 4 | 0 |
| 5 | 0 |
| 7 | 1 |
| 10 | 0 |
Untuk menjawab pertanyaan ini, mari kita mulai dengan contoh kecil dari korban Titanic yang selamat: hanya 10 orang dan hanya 10 hari di lautan setelah kecelakaan dan menganalisisnya dengan metode statistik klasik, regresi linier dan logistik. Kami tertarik pada berapa banyak orang dan berapa lama selamat dan berapa banyak yang meninggal pada hari tertentu. Dengan demikian, dua parameter utama yang menarik bagi kami adalah (1) waktu (dalam hari) dan (2) “status hidup” orang, di mana “status = 1” berarti kematian dan “status = 0” berarti tidak meninggal. Dan dengan memiliki status orang, kita bisa saja menghitung orang yang selamat dan tidak selamat pada setiap hari tertentu untuk memahami apa yang sedang terjadi, bukan?
ggplot(d, aes(time, fill = factor(status)))+
geom_bar(position = position_dodge())+
theme_bw()
tidak sepenuhnya. Selain itu, fakta bahwa plot dengan hitungan kami tidak terlalu mengungkapkan, kami memiliki dua masalah dengan itu:
Yang pertama adalah bahwa “status = 0” tidak selalu berarti selamat. Pertama, orang tersebut mungkin terbawa ombak dan menghilang. Kita kehilangan orang tersebut dan kita tidak tahu statusnya. Mungkin sudah meninggal, mungkin sudah diselamatkan oleh kapal nelayan atau masih hidup di suatu tempat di lautan. Kami tidak yakin. Satu-satunya hal yang kami yakini adalah bahwa orang ini tidak dapat dihitung sebagai orang yang meninggal atau selamat. Kedua, jika ada orang yang selamat selama 10 hari, mereka mungkin akan meninggal pada hari ke-11. Jadi, waktu memainkan peran di sini dan membawa kita ke masalah kedua.
Kita tidak memiliki data untuk setiap 10 hari tersebut. Hanya untuk 7 hari. Selain itu, jika saya tidak makan apa pun selama 5 hari, probabilitas kelangsungan hidup saya pada hari ke-5, di mana saya pasti selamat, tidak 100%. Itu lebih rendah, karena rasa lapar terakumulasi dari waktu ke waktu. Demikian pula, penurunan probabilitas kelangsungan hidup saya juga terakumulasi dari waktu ke waktu, yang tidak dapat ditentukan oleh penghitungan sederhana dari yang selamat vs yang tidak selamat pada hari tertentu.
Meskipun ada masalah ini, kita masih dapat menganalisis data ini. Kita hanya perlu menemukan metode yang tepat untuk itu.
- Melihat variabel numerik “waktu”, kita mungkin tergoda untuk menggunakan regresi linier untuk memodelkan waktu bertahan hidup untuk dua kelompok status yang berbeda, 0 & 1:
m <- lm(time ~ status, d %>% mutate(status = factor(status)))
library(effects)
plot(allEffects(m))
Namun, ketika kita memplot hasil modelnya, kita menyadari bahwa ada sesuatu yang tidak beres di sana. Semua informasi yang kita miliki untuk memodelkan waktu bertahan hidup adalah angka nol dan satu. Dan waktu rata-rata kematian (status = 1) sebesar 2,5 hari melewatkan banyak informasi. Sebagai contoh, setengah dari orang-orang (3 orang) meninggal pada hari pertama! Dan semakin jauh kita melangkah ke depan, semakin sedikit orang yang meninggal, yang dapat kita lihat pada plot hitungan pertama di atas. Jadi, waktu bertahan hidup tidak benar-benar linier dan tentu saja tidak dapat digambarkan hanya dengan angka nol dan satu! Jadi, model linier tampaknya tidak berguna di sini. Lalu apa yang bisa digunakan?
- Karena kita memiliki angka nol dan satu pada kolom status, kita dapat menggunakan regresi logistik, bukan?! Ya! Terutama karena kurva probabilitasnya tidak linier dan akan menangkap tren kelangsungan hidup dari waktu ke waktu. Keren! Mari kita mulai:
m <- glm(status ~ time, d, family = binomial)
library(sjPlot)
plot_model(m, type = "pred", ci.lvl = NA)## $time
Sekarang kita memiliki probabilitas kelangsungan hidup non-linear :), yang merupakan peningkatan besar dibandingkan dengan regresi linier. Namun, pengamatan kelima memiliki status 0 pada hari ke-2. Ini berarti bahwa orang tersebut pasti bertahan hidup selama 2 hari, tetapi hilang setelahnya. “Lalu kenapa?”, - Anda mungkin bertanya. Probabilitas dalam regresi logistik akan menghitung orang ini sebagai 100% selamat setelah hari kedua hanya karena statusnya tidak “=1”, yang salah (atau bias), karena kita tidak tahu apakah orang ini selamat. Dengan demikian, regresi logistik melebih-lebihkan probabilitas kelangsungan hidup dan oleh karena itu juga merupakan alat yang tidak tepat untuk menganalisis data kelangsungan hidup.
Jadi, sepertinya kita tidak dapat menganalisis data ketahanan hidup dengan metode klasik. Dan itulah mengapa kita membutuhkan metode analisis survival.
Sementara model regresi linier menggambarkan waktu, tetapi melewatkan probabilitas kelangsungan hidup non-linier, regresi logistik menangkap tren non-linier, tetapi melebih-lebihkan probabilitas kelangsungan hidup. Sebaliknya, analisis kelangsungan hidup memecahkan kedua masalah tersebut, karena model ini memodelkan probabilitas kelangsungan hidup non-linear dari waktu ke waktu sambil memperhitungkan subjek penelitian yang hilang, baik yang meninggal maupun yang selamat.
Kematian bukanlah satu-satunya pilihan! Atau: “Apakah yang dimaksud dengan peristiwa?”
Analisis waktu bertahan hidup diperlukan dalam penelitian apa pun yang menyelidiki waktu untuk hasil tertentu yang menarik. Penelitian kanker dalam bidang kedokteran dan kegagalan pertama mobil dalam bidang teknik (analisis waktu kegagalan) adalah contoh yang baik. Hasil yang diinginkan dapat berupa kematian, remisi untuk kambuh, perkembangan, atau kegagalan. Titik waktu untuk mencapai hasil tersebut umumnya disebut event**. Syukurlah, tidak semua “kejadian” berakibat fatal 😃, tetapi kadang-kadang bahkan bisa menjadi hasil yang menguntungkan seperti keluar dari rumah sakit. Dan dengan demikian, analisis ketahanan hidup juga merupakan istilah umum, karena ini bukan hanya tentang kelangsungan hidup.
Kejadian, sebagai titik waktu yang pasti diperlukan karena kita tidak dapat mengamati atau bereksperimen selamanya. Jika kita mempelajari kanker paru-paru, kita tidak bisa menunggu sampai pasien meninggal karena penyebab lain, misalnya karena sudah tua. Hanya ketika beberapa pasien selamat dari kanker, atau beberapa mobil tidak rusak pada suatu waktu (kejadian) tertentu, kita bisa mendapatkan wawasan yang berharga. Misalnya, berapa lama waktu dari awal pengobatan hingga perkembangannya, atau berapa probabilitas untuk bertahan hidup dari kanker paru-paru setelah tepat 1 tahun? Selain itu, hal ini membuka kemungkinan untuk membandingkan waktu bertahan hidup atau kegagalan di antara kelompok yang berbeda, misalnya kanker paru-paru antara perokok dan non-perokok, atau waktu kerusakan antara mobil Jerman dan Korea.
Jadi, waktu bertahan hidup memiliki titik awal dan titik akhir - peristiwa. Namun, bagaimana jika seorang pasien mengundurkan diri dari penelitian karena alasan pribadi, atau mobilnya dicuri 1 hari sebelum acara? Apakah mereka akan “bertahan hidup” pada saat acara berlangsung? Kita tidak tahu! Tapi yang pasti mereka bertahan dari awal penelitian hingga saat kami kehilangan mereka. Dan ini adalah informasi berharga yang tentunya ingin kami sertakan dalam analisis kami! Tapi bagaimana kita membedakan penyintas yang “hilang” dari penyintas yang “nyata”? Yah, kita sebut saja mereka dengan nama baru - disensor. Kata yang keras untuk diterapkan pada seseorang, bukan? Namun, ketika saya berpikir tentang informasi yang tidak pantas dalam sebuah buku atau sumpah serapah di TV, yang sering kali disensor, konsep menyensor pasien menjadi lebih mudah dicerna oleh saya. Lihat, karena kita tidak dapat mengatakan bahwa orang itu telah meninggal atau masih hidup, hanya karena kita TIDAK MEMILIKI INFORMASI YANG CUKUP, kedua kesimpulan itu akan menjadi tidak pantas, dan dengan demikian, kita menyensor orang ini. Konsep penyensoran ini sangat penting, sehingga perlu satu bab tersendiri.
Censoring
Kebutuhan untuk menyensor muncul dari fakta bahwa waktu bertahan hidup yang sebenarnya tidak akan diketahui oleh beberapa individu. Ada banyak cara untuk menghilangkan subjek penelitian:
- kita bisa menghilangkan individu tanpa alasan, ketika individu itu tidak muncul lagi
- individu tersebut mungkin mengalami peristiwa lain, misalnya kematian atau kecelakaan
- bahkan pasien yang bertahan hidup sampai akhir penelitian dapat diperlakukan sebagai disensor karena waktu bertahan hidup yang tidak diketahui (saya pribadi lebih suka melihat pasien-pasien ini sebagai survived)
Semua contoh di atas dianggap sebagai sensor kanan karena arahnya dari kiri ke kanan pada sumbu waktu. Sebagian besar data ketahanan hidup disensor ke kanan. Ada dua jenis penyensoran lainnya. **Penyensoran kiri muncul jika kita tidak tahu di mana penyakit dimulai, sedangkan penyensoran interval terjadi ketika waktu pasti hilangnya pasien tidak diketahui, tetapi hanya jendela waktu. Baik penyensoran kiri maupun penyensoran interval jarang terjadi, sulit untuk dianalisis, dan sering kali merupakan hasil dari desain penelitian yang buruk, dan dengan demikian, tidak akan dibahas di sini.
Yang penting untuk diingat adalah bahwa pengamatan yang disensor masih memberikan informasi yang berguna! Itulah mengapa mereka perlu dimasukkan ke dalam analisis.
Cara menghitung kurva survival Kaplan-Meier “secara manual” selangkah demi selangkah
Bayangkan hari yang cerah dan damai tepat sebelum kecelakaan Titanic terjadi. Tidak ada yang meninggal, tetapi semua 10 orang berisiko tinggi (bahaya) meninggal, hanya saja mereka belum mengetahuinya. Coba tebak, berapa banyak orang yang mungkin akan meninggal sehari sebelum kecelakaan? Jawabannya adalah - mungkin nol. Secara harfiah, kemungkinan kematian adalah nol, karena kita tahu bahwa kecelakaan tidak akan terjadi sehari sebelum tabrakan yang sebenarnya. Untuk menjawab pertanyaan ini dengan lebih tepat, kita perlu mengingat definisi probabilitas: Probabilitas adalah rasio dari sesuatu yang terjadi, terhadap segala sesuatu yang dapat terjadi. Jadi, jika tidak ada seorang pun dari 10 orang yang meninggal sehari sebelum kecelakaan, maka probabilitas kematiannya adalah 0%, sedangkan probabilitas selamat adalah 100%:
\[ probability \ of \ dying = \frac{0}{10} = 0\]
\[ probability \ of \ surviving = \frac{10}{10} = 1 = 100\%\]
Untuk representasi yang lebih baik, mari kita letakkan semua angka dalam satu tabel:
tribble(
~time_in_days, ~N_at_risk, ~N_died, ~P_dying, ~`P_surviving`,
0, 10, 0, "0/10 = 0", "10/10 = 1"
) %>% kable()| time_in_days | N_at_risk | N_died | P_dying | P_surviving |
|---|---|---|---|---|
| 0 | 10 | 0 | 0/10 = 0 | 10/10 = 1 |
Kemudian, pada hari kecelakaan, di mana ke-10 orang tersebut sekarang mengetahui bahwa mereka berisiko tinggi untuk meninggal, 3 di antaranya benar-benar meninggal. Probabilitas kematian pada hari pertama adalah \(3/10 = 0.3 \approx 30\%\) dan probabilitas bertahan hidup adalah \(7/10 = 0.7 \approx 70\%\). Seperti yang Anda lihat, kelangsungan hidup juga dapat dihitung dari probabilitas kematian \(1-0.3 = 0.7\), yang terkadang sangat berguna.
Survival probability
Satu momen penting di sini adalah bahwa kemungkinan untuk bertahan hidup bersifat kumulatif, yang berarti terakumulasi dari hari ke hari. Akumulasi ini terjadi melalui mengalikan probabilitas baru untuk bertahan hidup pada hari ke-1 \((1 - \frac{d_i}{n_i}) = (1 - \frac{3}{10}) = 70\%\) dengan probabilitas lama untuk bertahan hidup sepanjang waktu sebelumnya, yang dalam kasus kami adalah sehari sebelum kecelakaan, atau - hari ke-nol \(S(t_{i-1}) = 100\%\):
\[ S(t_i) = S(t_{i-1})*(1 - \frac{d_i}{n_i}) = 1 * (1 - \frac{3}{10}) = 0.7 \]
Dimana,
- \(S(t_{i−1})\) = probabilitas untuk hidup di \(t_{i−1}\)
- \(n_i\) = jumlah pasien yang masih hidup sebelum \(t_i\)
- \(d_i\) = Jumlah kejadian saat \(t_i\)
- \(t_0\) = 0, \(S(0)\) = 1
Probabilitas bertahan hidup pada waktu tertentu, \(S(t)\), adalah syarat karena orang tersebut harus bertahan hidup setelah waktu tertentu, misalnya hari ke-nol (itulah syaratnya) untuk tetap berada dalam eksperimen untuk hari pertama.
tribble(
~time_in_days, ~N_at_risk, ~N_died, ~P_dying, ~`P_surviving`,
0, 10, 0, "0/10 = 0", "10/10 = 1",
1, 10, 3, "3/10 = 0.3", "1 * (7/10) = 0.7",
) %>% kable()| time_in_days | N_at_risk | N_died | P_dying | P_surviving |
|---|---|---|---|---|
| 0 | 10 | 0 | 0/10 = 0 | 10/10 = 1 |
| 1 | 10 | 3 | 3/10 = 0.3 | 1 * (7/10) = 0.7 |
Hari kedua sedikit lebih menarik, karena satu orang meninggal dan satu orang hilang begitu saja (status = 0). Orang yang hilang ini kemungkinan besar terbawa ombak dan mudah-mudahan bisa diselamatkan. Harapannya besar, karena semakin jauh orang tersebar di lautan, semakin besar kemungkinan salah satu dari mereka akan ditemukan dalam keadaan hidup. Dan hal ini akan memberi tahu dunia bahwa ada orang lain yang masih berada di luar sana. Jadi, meskipun orang hilang terlihat buruk, hal ini bisa jadi merupakan hal yang baik. Jadi, kita tentu tidak bisa menganggap orang yang hilang sebagai orang yang sudah mati. Selain itu, seperti kata pepatah Jerman: “Harapan adalah yang terakhir mati” 😉.
Tetapi karena kita juga tidak yakin apakah orang ini masih hidup, kita tidak bisa mengatakan - orang tersebut selamat. Hal ini menimbulkan dilema: meskipun tidak meninggal, orang tersebut bukan bagian dari eksperimen kami lagi dan kami harus menghapusnya dari jumlah orang yang berisiko. Jadi, sementara hari kedua akan memiliki 7 orang yang berisiko yang tersisa (karena 3 orang telah meninggal di hari pertama), hari ketiga akan tersisa hanya 5 orang, karena 1 orang meninggal dan satu orang menghilang (disensor) dibandingkan dengan hari kedua. Sama halnya dengan hari pertama, probabilitas kelangsungan hidup bersifat kumulatif dan selalu menyertakan probabilitas hari sebelumnya.
tribble(
~time_in_days, ~N_at_risk, ~N_died, ~P_dying, ~`P_surviving`,
0, 10, 0, "0/10 = 0", "10/10 = 1",
1, 10, 3, "3/10 = 0.3", "1 * (7/10) = 0.7",
2, 7, 1, "1/7 = 0.14", "0.7 * (6/7) = 0.6",
3, 5, 1, "1/5 = 0.20", "0.6 * (4/5) = 0.48"
) %>% kable()| time_in_days | N_at_risk | N_died | P_dying | P_surviving |
|---|---|---|---|---|
| 0 | 10 | 0 | 0/10 = 0 | 10/10 = 1 |
| 1 | 10 | 3 | 3/10 = 0.3 | 1 * (7/10) = 0.7 |
| 2 | 7 | 1 | 1/7 = 0.14 | 0.7 * (6/7) = 0.6 |
| 3 | 5 | 1 | 1/5 = 0.20 | 0.6 * (4/5) = 0.48 |
Dua hari berikutnya, dua orang menghilang, jadi kita menyensor mereka. Tidak perlu menghitung sesuatu untuk data yang disensor dalam analisis kelangsungan hidup, karena kita tidak tahu apakah mereka selamat atau tidak. Dan karena kita tertarik pada kelangsungan hidup atau kematian, tabel akhir di bawah ini berisi kasus yang hanya mati. Pada hari ke-6 tidak ada yang mati atau menghilang, oleh karena itu kami juga tidak menghitung apa pun untuk hari itu. Orang terakhir dalam eksperimen kami meninggal pada hari ke-7. Jadi, “N_meninggal” akan menjadi 1, dan “N_berisiko” akan menjadi 2, karena 1 dari 5 orang yang tersisa telah meninggal dan 2 orang disensor. Dengan demikian, tabel yang dihitung secara “manual” akan terlihat seperti ini:
tribble(
~time_in_days, ~N_at_risk, ~N_died, ~P_dying, ~`P_surviving`,
0, 10, 0, "0/10 = 0", "10/10 = 1",
1, 10, 3, "3/10 = 0.3", "1 * (7/10) = 0.7",
2, 7, 1, "1/7 = 0.14", "0.7 * (6/7) = 0.6",
3, 5, 1, "1/5 = 0.20", "0.6 * (4/5) = 0.48",
7, 2, 1, "1/2 = 0.50", "0.48 * (1/2) = 0.24"
) %>% kable()| time_in_days | N_at_risk | N_died | P_dying | P_surviving |
|---|---|---|---|---|
| 0 | 10 | 0 | 0/10 = 0 | 10/10 = 1 |
| 1 | 10 | 3 | 3/10 = 0.3 | 1 * (7/10) = 0.7 |
| 2 | 7 | 1 | 1/7 = 0.14 | 0.7 * (6/7) = 0.6 |
| 3 | 5 | 1 | 1/5 = 0.20 | 0.6 * (4/5) = 0.48 |
| 7 | 2 | 1 | 1/2 = 0.50 | 0.48 * (1/2) = 0.24 |
Namun, apa yang terjadi jika kita mengabaikan penyensoran dan hanya menghitung probabilitas untuk tidak mati? Nah, karena 4 orang tidak mati pada hari ke-7, probabilitas untuk bertahan hidup akan menjadi \(\frac{4}{10} = 40\%\) yang hampir dua kali lipat lebih tinggi, dan dengan demikian kita akan melebih-lebihkan probabilitas untuk bertahan hidup, dibandingkan dengan probabilitas sebenarnya yaitu 0,24 yang memperhitungkan penyensoran!
Story time:
Hanya memperhatikan orang yang selamat bahkan memiliki nama - bias penyintas, dan ada cerita kecil di baliknya. Selama Perang Dunia Kedua, beberapa pesawat kembali dari medan perang dengan banyak kerusakan akibat peluru. Mereka hampir tidak bisa terbang, tapi mereka tetap kembali. Jadi, pihak militer memutuskan untuk melindungi pesawat dengan lebih banyak pelindung di bagian yang paling banyak lubang peluru, seperti sayap, dan mengurangi pelindung di bagian yang tidak berlubang peluru. Anehnya, persentase pesawat yang kembali tidak meningkat. Para insinyur bingung! Hingga seorang ahli matematika, Abraham Walt, yang diundang untuk memecahkan masalah ini berkata: “Pasang lebih banyak lapis baja di tempat-tempat yang tidak memiliki lubang peluru, karena jika tempat ini ditembak, pesawat tidak akan kembali.” Dan ketika yang lain memikirkannya, mereka menyadari bahwa semua pesawat yang kembali tidak memiliki lubang peluru di kokpit atau mesin. Lubang peluru menunjukkan semua tempat di mana pesawat bisa ditembak tapi masih bisa kembali, atau selamat! Itulah yang disebut dengan bias bertahan hidup. Jadi, menurut saya, regresi logistik memiliki survivorship bias dibandingkan dengan analisis ketahanan hidup jika kita ingin menganalisis data ketahanan hidup.
Cara menghitung kurva survival Kaplan-Meier
Muat semua paket yang dibutuhkan sekaligus untuk menghindari interupsi.
library(tidyverse) # data wrangling and visualization
library(knitr) # beautifying tables
library(car) # for checking assumptions, e.g. vif etc.
library(broom) # for tidy model output
library(sjPlot) # for plotting results of log.regr.
library(sjmisc) # for plotting results of log.regr.
library(effects) # for probability output and plotsPertama-tama, pastikan Anda mengetahui dengan pasti apa arti angka 0
dan 1 dalam data Anda! Karena dalam regresi logistik 1 adalah
kelangsungan hidup, sedangkan dalam analisis ketahanan hidup 1
adalah kematian! Kedua, kita harus membedakan kasus tersensor dengan
cara tertentu, misalnya kita dapat menandainya dengan tanda
plus pada data atau pada plot. Data di bawah ini menunjukkan
bahwa orang-orang disensor pada hari ke-2+, 4+, 5+, dan 10+. Instal dan
muat package survival untuk dapat menjalankan kode di bawah
ini:
# install.packages("survival")
library(survival)
Surv(time = d$time, event = d$status)## [1] 1 1 1 2 2+ 3 4+ 5+ 7 10+Fungsi Surv menyatukan data “waktu” dan “status” ke
dalam satu objek “survival” tunggal, yang memungkinkan untuk
“memperhitungkan pengamatan yang disensor”. Objek ini kemudian dapat
digunakan untuk memodelkan probabilitas kelangsungan hidup dengan fungsi
survfit. Kita memodelkan kelangsungan hidup dengan
menambahkan “~1” pada objek tersebut, di mana 1 berarti tidak ada
variabel yang dapat mempengaruhi kelangsungan hidup. Objek
kelangsungan hidup kita di sisi kiri tilde adalah variabel
respons dan di sisi kanan ~ (tilde) adalah prediktor, atau
dalam kasus “1” - tidak ada. Mari kita buat model survival pertama kita
dan lihatlah keluaran modelnya:
survival_model <- survfit(Surv(time, status) ~ 1, data = d)
# small summary
survival_model## Call: survfit(formula = Surv(time, status) ~ 1, data = d)
##
## n events median 0.95LCL 0.95UCL
## [1,] 10 6 3 1 NAwhere we have:
- n - jumlah peserta,
- event - jumlah kematian, yang secara umum disebut event karena Anda dapat mempelajari jenis-jenis “deadlines” lainnya, seperti kambuhnya penyakit, atau kerusakan pertama pada mobil,
- perkiraan waktu bertahan hidup median, bukan rata-rata karena sifat non-parametrik (dijelaskan nanti) dari analisis ketahanan hidup dan
- Interval kepercayaan 95%** untuk waktu ketahanan hidup median
Lebih banyak informasi dapat dinilai dengan fungsi
ringkasan, yang secara mengejutkan 😉 memberikan tabel yang
sama dengan yang baru saja kita hitung secara manual di atas:
# big summary
summary(survival_model)## Call: survfit(formula = Surv(time, status) ~ 1, data = d)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 10 3 0.70 0.145 0.4665 1.000
## 2 7 1 0.60 0.155 0.3617 0.995
## 3 5 1 0.48 0.164 0.2458 0.938
## 7 2 1 0.24 0.188 0.0515 1.000Sekarang Anda bisa yakin, perhitungan kami benar dan Anda telah
berhasil mempelajari cara melakukan analisis survival. Selamat! Sehingga
mulai sekarang Anda dapat mulai menggunakan perangkat lunak. Perangkat
lunak ini memiliki banyak keuntungan! Sebagai contoh, fungsi
summary, yang hanya menampilkan hasil dari kasus yang mati,
dapat diisi dengan argumen tambahan censored = TRUE, yang
menampilkan semua kasus, baik yang mati maupun yang disensor:
summary(survival_model, censored = T)## Call: survfit(formula = Surv(time, status) ~ 1, data = d)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 10 3 0.70 0.145 0.4665 1.000
## 2 7 1 0.60 0.155 0.3617 0.995
## 3 5 1 0.48 0.164 0.2458 0.938
## 4 4 0 0.48 0.164 0.2458 0.938
## 5 3 0 0.48 0.164 0.2458 0.938
## 7 2 1 0.24 0.188 0.0515 1.000
## 10 1 0 0.24 0.188 0.0515 1.000Bagus! Benar kan? Tetapi bagaimana jika Anda memiliki ribuan hari
(tabel Anda akan sangat besar!), tetapi Anda hanya tertarik pada
beberapa hari saja? Nah, Anda dapat menentukan hasil dari hari mana yang
ingin Anda lihat dengan menambahkan argumen times = ... ke
fungsi summary. Hal yang keren 😎 tentangnya adalah Anda
bahkan dapat meminta hari-hari yang tidak tersedia datanya, misalnya
hari ke-8 atau ke-9 dalam contoh sederhana kita. Sebagai contoh,
probabilitas kelangsungan hidup setelah hari ke-8 adalah 24%:
summary(survival_model, times = c(8, 9))## Call: survfit(formula = Surv(time, status) ~ 1, data = d)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 8 1 6 0.24 0.188 0.0515 1
## 9 1 0 0.24 0.188 0.0515 1Sekarang, karena kita mengetahui probabilitas yang tepat untuk
bertahan hidup pada setiap hari, bahkan pada hari-hari
dimana kita tidak memiliki datanya, kita dapat memvisualisasikan hasil
model. Instal dan muat paket survminer untuk dapat
menjalankan kode di bawah ini:
# install.packages("survminer")
library(survminer)
ggsurvplot(
survival_model,
conf.int = FALSE,
surv.median.line = "hv",
xlab = "Days",
ylab = "Survival probability",
break.time.by = 1,
risk.table = T
)
Terlepas dari kenyataan bahwa kita dapat dengan mudah menghitung
probabilitas untuk setiap hari, Anda tidak akan pernah bisa melakukannya
dengan tangan. Bisa jadi ribuan hari. Oleh karena itu, Anda akan
membiarkan perangkat lunak melakukan pekerjaannya dan Anda akan
mendapatkan lebih banyak hasil daripada yang dapat dan ingin Anda hitung
sendiri, misalnya interval kepercayaan (CI) untuk kelangsungan hidup
sehari-hari atau plot kelangsungan hidup. CI divisualisasikan secara
default, tetapi kita dapat menghapusnya jika diperlukan melalui perintah
conf.int = FALSE. Baris terakhir pada kode di atas
menampilkan tabel “Number at risk”, yang juga telah kita hitung secara
manual di atas. Jadi, mengapa kita menghitung sesuatu secara manual jika
kita bisa mendapatkan semuanya dan lebih banyak lagi dari perangkat
lunak? Nah, penting untuk melalui proses perhitungan langkah demi
langkah untuk meningkatkan intuisi Anda tentang kurva analisis
survival!
Menariknya, kurva ini digambarkan secara independen oleh dua ilmuwan
yang berbeda pada waktu yang sama. Edward Kaplan dan Paul Meier kemudian
mempublikasikan temuan mereka bersama-sama sebagai penaksir
Kaplan-Meier (KM) pada tahun 1958.1 Dan fungsi
survfit cocok dengan kurva `type = “kaplan-meier” secara
default.
Sebagian besar data survival menunjukkan banyak kejadian di awal dan jumlah kejadian yang lebih rendah sepanjang waktu. Hal ini membuat kurva survival menjadi tidak linier dan sebagian besar data survival terdistribusi secara miring atau tidak normal (tidak berbentuk lonceng). Ini adalah alasan mengapa metode KM mengestimasi rata-rata waktu bertahan hidup dan bukannya mean sebagai ukuran kecenderungan sentral. Ketidaklinieran dan “stepiness” dari kurva KM membuatnya tidak mungkin untuk meringkas data ke dalam satu parameter tunggal, misalnya slope, yang membuat metode KM menjadi non-parametrik.
Interpretasi Kurva Kaplan-Meier
Sumbu x mewakili waktu dalam hari, dan sumbu
y menunjukkan probabilitas bertahan hidup atau proporsi orang
yang bertahan hidup. Jadi, kurva itu sendiri
menunjukkan kemungkinan bertahan hidup yang tepat dari waktu ke
waktu. Sebuah turunan vertikal pada kurva
menunjukkan setidaknya satu peristiwa. Ketinggian
turunan vertikal menunjukkan perubahan probabilitas
kelangsungan hidup kumulatif. Bagian horizontal dari
kurva menunjukkan durasi ketahanan hidup untuk interval waktu
tertentu, yang diakhiri oleh kejadian berikutnya (dan penurunan
kurva). Kurva KM terlihat seperti tangga yang aneh
dengan anak tangga yang tidak rata, di mana
probabilitas kelangsungan hidup konstan di antara berbagai kejadian, dan
oleh karena itu merupakan fungsi langkah yang hanya
berubah nilainya pada saat setiap kejadian. Dengan cara ini, setiap
pasien menyumbangkan informasi yang berharga untuk perhitungan selama
pasien tersebut masih hidup. Orang yang disensor ditampilkan
persis seperti pada objek bertahan hidup, dengan plus,
seperti yang Anda lihat pada hari ke-2. Namun, sebagian besar nilai
tambah terlihat seperti tanda centang vertikal, karena
mereka berada di bagian horizontal kurva, yaitu hari ke-4 dan ke-5.
Tanda centang ditampilkan secara default, namun dapat ditekan dengan
menggunakan argumen censor = FALSE. Tabel risiko di bawah
plot menunjukkan jumlah orang yang berisiko, yang
sebenarnya semua orang yang “benar-benar hidup” dalam
percobaan** yang tidak mengalami kejadian atau penyensoran pada titik
waktu tertentu. Garis putus-putus menunjukkan rata-rata
waktu bertahan hidup** yang sesuai dengan kemungkinan bertahan
hidup 50%. Dan jika kita mengabaikan penyensoran dan hanya
memperkirakan median waktu (hanya untuk orang yang meninggal), kita akan
mendapatkan 1,5, bukan 3, yang akan meningkatkan provabilitas
kelangsungan hidup dari 48% menjadi 70%. Sekali lagi, mengabaikan
penyensoran akan menghasilkan estimasi yang terlalu tinggi dari
probabilitas kelangsungan hidup karena adanya bias kelangsungan
hidup.
d %>%
filter(status == 1) %>%
summarize(median_survival_whithout_censoring = median(time))Membandingkan survival of groups
2 groups
Sekarang, mari kita ikuti 100 orang setelah kecelakaan Titanic, bukan
10 orang, dan lihat kelangsungan hidup mereka menggunakan kurva
Kaplan-Meier di bawah ini. Angka-angka dalam tabel risiko semakin besar
dan kita sebaiknya menampilkan persentase risiko dalam tanda kurung di
dekat nilai absolut. Argumen risk.table = "abs_pct" akan
membantu dalam hal ini. Contoh “bulat” dari 100 orang memberikan kita
persentase yang sama persis dengan angka absolut, namun dengan jumlah
yang kurang “bulat” atau beberapa kelompok, persentase akan menjadi
sangat berguna. Rata-rata waktu bertahan hidup manusia di lautan dingin
adalah sekitar 70 hari dan interval kepercayaan tidak terlalu lebar,
sehingga kita bisa cukup percaya diri dengan angka yang kita
dapatkan:
set.seed(999) # for reproducible example
d <- ggstatsplot::Titanic_full %>%
mutate(survived = ifelse(Survived == "No", 1, 0),
time = runif(n=2201, min=1, max=100)) %>%
sample_n(100)
m <- survfit(Surv(time, survived) ~ 1, data = d)
ggsurvplot(m,
conf.int = TRUE,
risk.table = "abs_pct",
surv.median.line = "hv")
Tapi bisakah kita memiliki dua kurva kelangsungan hidup pada plot
yang sama dan membandingkannya? Kasar! Sebenarnya, itulah saat di mana
kesenangan dimulai. Sebagai contoh, kita bisa membandingkan probabilitas
kelangsungan hidup jantan vs betina. Untuk itu, kita hanya perlu
mengganti ~1 pada rumus model dengan nama variabel
kategorikal yang diinginkan, misalnya sex. Plot survival
kemudian menampilkan kurva Kaplan-Meier untuk setiap kategori variabel
anda:
m <- survfit(Surv(time, survived) ~ Sex, data = d)
ggsurvplot(m,
pval = TRUE,
conf.int = TRUE,
risk.table = "abs_pct",
surv.median.line = "hv") 
Interpretasi perbandingan kelompok dengan menggunakan 4 tolok ukur
perbandingan visual kurva: hanya dengan melihat dua kelompok, kita dapat mengatakan apakah ada perbedaan. Sebagai contoh, plot kami menunjukkan bahwa betina memiliki probabilitas bertahan hidup yang lebih tinggi hampir di seluruh periode waktu (sumbu x). Namun, apakah perbedaan ini signifikan secara statistik memerlukan uji statistik formal. Kita dapat melangkah lebih jauh dan melihat angka-angka pada kurva, misalnya: pada hari yang cerah dan indah sebelum tabrakan, probabilitas kelangsungan hidup kedua kelompok adalah 1,0 (atau 100% penumpang masih hidup). Pada hari ke-50, probabilitas kelangsungan hidup betina adalah sekitar 0,75 (atau 75%) dan hanya sekitar 0,65 (atau 65%) untuk jantan. Pada hari ke-75, tingkat kelangsungan hidup adalah sekitar 60 dan sekitar 30%.
Perbandingan confidence interval (CI) : CI yang tumpang tindih menunjukkan bahwa kelangsungan hidup jantan dan betina tidak terlalu berbeda dan dapat disebabkan oleh peluang hingga sekitar 80 hari. Setelah 80 hari, CI tidak lagi tumpang tindih yang menunjukkan perbedaan signifikan dalam kelangsungan hidup pada titik waktu tertentu, misalnya 90 hari.
Perkiraan waktu kelangsungan hidup rata-rata menunjukkan lebih dari sekadar interval kepercayaan. Waktu kelangsungan hidup rata-rata untuk setiap kelompok mewakili waktu di mana probabilitas kelangsungan hidup adalah 50%. Sebagai contoh, median kelangsungan hidup betina jauh lebih tinggi (87 hari) daripada jantan (65 hari). Perbedaan yang sangat besar sebesar 23 hari terdengar signifikan bagi saya, karena betina memiliki waktu 23 hari lebih lama untuk ditemukan oleh kapal penangkap ikan. Di sini, sekali lagi, hanya tes yang dapat memberi tahu dan itulah mengapa tolok ukur terakhir untuk membandingkan dua kelompok adalah p-value yang diestimasi oleh Tes Log-Rank (Peto dkk, 1977), yang sangat penting, sehingga perlu mendapat satu bab tambahan.
m## Call: survfit(formula = Surv(time, survived) ~ Sex, data = d)
##
## n events median 0.95LCL 0.95UCL
## Sex=Female 31 10 87.0 60.9 NA
## Sex=Male 69 54 64.5 53.6 72.9Log-Rank test
- Uji statistik non-parametrik Log-Rank (kadang-kadang disebut Mantel-Haenszel) membandingkan median waktu kelangsungan hidup kelompok. Uji Log-Rank mirip dengan uji Wilcoxon-Rank non parametrik (1), yang juga membandingkan median dengan menggunakan peringkat dan juga mirip dengan uji Chi-square (2), di mana kita membandingkan jumlah kejadian yang diamati dengan yang diharapkan dengan menghitung statistik Chi-square:
\[ X^2 = \sum_{i = 1}^{g} \frac{(O_i - E_i)^2}{E_i}\]
Jumlah kejadian yang diharapkan dihitung untuk setiap titik waktu dan setiap kelompok dibandingkan dengan titik waktu sebelumnya. Nilai-nilai ini kemudian dijumlahkan pada semua titik waktu untuk memberikan jumlah total kejadian yang diharapkan pada setiap kelompok.
Sifat non-parametrik dari tes ini membuat tidak ada asumsi tentang distribusi kelangsungan hidup. Jadi, apakah waktu kelangsungan hidup normal (“lonceng melengkung”)? Bisa saja, tetapi tidak harus! Faktanya, data kelangsungan hidup sangat jarang berdistribusi normal, tetapi sering kali miring karena biasanya banyak kejadian awal dan relatif sedikit kejadian akhir.
Hipotesis nol dari uji Log-Rank adalah tidak ada perbedaan dalam ketahanan hidup antara kedua kelompok. Nilai p-value sebesar 0,011 memungkinkan kita untuk menolak hipotesis nol dan mengindikasikan adanya perbedaan rata-rata yang signifikan dalam waktu kelangsungan hidup antara perempuan dan laki-laki.
Uji Log-Rank sangat banyak digunakan untuk membandingkan dua atau lebih kurva kelangsungan hidup, sehingga Anda harus benar-benar membenarkan penggunaan uji lainnya.
Fungsi survdiff() menghitung uji Log-Rank dan
mengembalikan komponen-komponen berikut:
- jumlah subjek dalam setiap kelompok.
- jumlah kejadian teramati tertimbang **dalam setiap kelompok.
- jumlah kejadian tertimbang ** yang diharapkan** dalam setiap kelompok.
- Statistik chi-square untuk uji kesetaraan
- dan p-value untuk perbedaan kelangsungan hidup di antara kelompok
survdiff(Surv(time, survived) ~ Sex, data = d)## Call:
## survdiff(formula = Surv(time, survived) ~ Sex, data = d)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Sex=Female 31 10 19.3 4.45 6.46
## Sex=Male 69 54 44.7 1.92 6.46
##
## Chisq= 6.5 on 1 degrees of freedom, p= 0.01Apakah tes Log-Rank sempurna? Sayangnya, tidak. Salah satu masalahnya adalah: uji ini tidak memberikan ukuran efek yang membuat kita hanya memiliki p-value sebagai ukuran perbedaan. Masalah lain dengan uji Log-Rank (dan metode Kaplan-Meier secara umum) adalah bahwa uji ini tidak mengizinkan adanya perancu, uji ini mengabaikan faktor/variabel lain. Untungnya, kita masih bisa membandingkan lebih dari dua kelompok dengan satu variabel 😉.
> 2 kelompok dan beberapa uji Log-Rank berpasangan (post-hoc)
Uji Log-Rank dapat membandingkan lebih dari dua kelompok dan menyatakan apakah ada perbedaan yang signifikan (p-value < 0,05) di antara kelompok-kelompok tersebut. Hasil uji yang ditampilkan pada tabel dan kurva KM di bawah ini menunjukkan bahwa terdapat perbedaan yang signifikan dalam kelangsungan hidup (p-value = 0,024) di antara kelompok orang di kelas tiket yang berbeda. Namun, seperti kebanyakan tes lainnya (misalnya ANOVA), tes ini tidak menjelaskan secara pasti kelompok mana yang mana. Itulah mengapa kita membutuhkan analisis tambahan yang membandingkan setiap kelompok dengan kelompok lainnya secara berpasangan. Analisis seperti ini sering disebut post-hoc.
survdiff(Surv(time, survived) ~ Class, data = d)## Call:
## survdiff(formula = Surv(time, survived) ~ Class, data = d)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Class=1st 14 4 11.03 4.477 5.532
## Class=2nd 13 5 8.49 1.432 1.666
## Class=3rd 37 30 21.81 3.074 4.760
## Class=Crew 36 25 22.68 0.238 0.383
##
## Chisq= 9.5 on 3 degrees of freedom, p= 0.02m <- survfit(Surv(time, survived) ~ Class, data = d)
ggsurvplot(m,
pval = TRUE,
conf.int = FALSE,
risk.table = "abs_pct",
surv.median.line = "hv",
ncensor.plot = TRUE,
break.time.by = 20,
risk.table.y.text.col = TRUE, risk.table.y.text = FALSE) 
Tiga argumen plot yang berguna lebih lanjut yang dapat membantu memvisualisasikan data ketahanan hidup dengan lebih baik adalah:
break.time.by = 20memecah sumbu x menjadi interval waktu yang diinginkan.risk.table.y.text.col = TRUEdanrisk.table.y.text = FALSEmemplot batang, bukan nama dalam anotasi teks legenda tabel risikoncensor.plot = TRUEmenampilkan jumlah subjek yang disensor, yang membantu untuk memahami apa yang menjadi penyebab angka risiko menjadi lebih kecil: peristiwa atau penyensoran.
Interpretasi
Rata-rata kelangsungan hidup** adalah sekitar 90 hari untuk penumpang
di kelas 1 dan 2, 58 hari untuk kelas 3 dan 65 hari untuk kru Titanic,
menunjukkan kelangsungan hidup yang baik bagi orang-orang kaya,
dibandingkan dengan yang lainnya. Nilai p-value yang rendah (p = 0,024)
menunjukkan bahwa ada perbedaan yang signifikan dalam kelangsungan hidup
di antara kelompok-kelompok tersebut. Di antaranya, hanya post-hoc yang
dapat menunjukkannya. Hasil dari analisis post-hoc
Log-Rank dapat dilakukan dengan fungsi
pairwise_survdiff dari paket survminer dan
ditampilkan di bawah ini. Seperti biasa, jika kita memiliki
banyak perbandingan, kita menghadapi risiko membuat
penemuan yang salah atau kehilangan penemuan penting. Oleh karena itu,
kita harus menyesuaikan nilai-p untuk mengurangi
kemungkinan membuat kesalahan. Saya pribadi lebih memilih metode
penyesuaian Benjamini & Hochberg (1995) daripada metode
Bonferroni** yang terkenal namun konservatif:
pairwise_survdiff(
formula = Surv(time, survived) ~ Class, data = d, p.adjust.method = "fdr"
)##
## Pairwise comparisons using Log-Rank test
##
## data: d and Class
##
## 1st 2nd 3rd
## 2nd 0.414 - -
## 3rd 0.028 0.139 -
## Crew 0.095 0.317 0.414
##
## P value adjustment method: fdrHasil analisis post-hoc mengungkapkan bahwa kelangsungan hidup kelas 1 secara signifikan lebih tinggi dibandingkan dengan kelas 3 dan Awak. Dan terlepas dari kenyataan bahwa rata-rata kelangsungan hidup penumpang kelas 2 sangat mirip dengan kelas 1, interval kepercayaan kelangsungan hidup ini sangat lebar dan banyak tumpang tindih dengan kelompok lain (tidak ditunjukkan untuk menghindari kekacauan). Itulah mengapa penumpang kelas 2 tidak umumnya memiliki waktu bertahan hidup yang lebih tinggi secara signifikan meskipun memiliki waktu bertahan hidup yang jauh lebih tinggi.
Kurva survival berganda
Seperti yang telah disebutkan di atas, uji Log-Rank tidak dapat diterapkan pada beberapa variabel. Namun, kita masih bisa memplot kurva survival dari beberapa variabel untuk mendapatkan intuisi untuk data kita. Metode lebih lanjut, misalnya model Cox, dapat mengestimasi perbedaan di antara beberapa variabel, namun metode ini lebih kompleks dibandingkan metode KM dan akan dibahas di tulisan selanjutnya.
m2 <- survfit( Surv(time, survived) ~ Sex + Class, data = d )
ggsurv <- ggsurvplot(m2, conf.int = TRUE)
ggsurv$plot +theme_bw() +
theme (legend.position = "right")+
facet_grid(Sex ~ .)
Intuisi yang dapat kita peroleh dari plot di atas adalah bahwa wanita kaya (kelas tiket 1 dan 2) dan (hanya 2) wanita dari kru kapal memiliki probabilitas 100% untuk bertahan hidup, sementara wanita murni di kelas 3 akan mati dengan kepastian yang sama dengan pria. CI yang sangat tumpang tindih dari kelangsungan hidup pria menunjukkan bahwa semua pria pada akhirnya akan mati setelah kecelakaan Titanic.
Kesimpulan
Analisis ketahanan hidup menyelidiki waktu yang dibutuhkan untuk suatu peristiwa yang menarik untuk terjadi. Sebagian besar analisis ketahanan hidup univariat (variabel tunggal) menggunakan plot Kaplan-Meier untuk memvisualisasikan kurva ketahanan hidup dan Uji Log-Rank untuk membandingkan kurva ketahanan hidup dua atau lebih kelompok.
- Keuntungan**:
Keuntungan penting dari kurva probabilitas kelangsungan hidup vs kurva regresi logistik adalah memperhitungkan data tersensor, yang tidak mati atau hidup. Regresi logistik memperlakukan semua orang yang tidak meninggal sebagai orang yang selamat, dan ini salah, karena kita tidak tahu status kelangsungan hidup orang yang tidak meninggal. Pasien meninggalkan penelitian karena dua alasan utama, mereka merasa sangat tidak nyaman, sehingga mereka tidak peduli lagi dengan eksperimen Anda, atau mereka merasa jauh lebih baik dan melupakan penelitian Anda. Beberapa pasien mungkin hanya pindah ke kota lain tanpa mengatakan apa-apa. Menghitung mereka semua sebagai pasien yang selamat, seperti yang dilakukan oleh regresi logistik, akan melebih-lebihkan probabilitas kelangsungan hidup (bias kelangsungan hidup) dan meremehkan bahaya kematian. Oleh karena itu, analisis ketahanan hidup lebih tepat dibandingkan dengan regresi logistik, meskipun masih menangkap tren non-linear dalam probabilitas.
Keuntungan lainnya adalah sifat non-parametrik dari metode Kaplan-Meier yang tidak memiliki terlalu banyak asumsi. Bahkan satu-satunya asumsi yang penting adalah bahwa penyensoran harus tidak informatif. Mengapa? Lebih banyak informasi selalu lebih baik, bukan? Ya! Tetapi jika kita tahu mengapa orang meninggalkan penelitian, kita dapat menggunakan informasi ini sebagai variabel baru dan mempelajari pengaruhnya terhadap kelangsungan hidup. Metode KM menjadi tidak tepat, karena akan kehilangan informasi ini. Namun, seringkali kita tidak tahu (no-info) mengapa orang keluar (disensor). Dan dalam hal ini metode KM dapat menarik kesimpulan yang paling banyak dari data yang tidak informatif tersebut.
- Kekurangan**:
“Kurva” Kaplan-Meier tidak benar-benar terlihat seperti kurva. Berlawanan dengan regresi logistik, metode Kaplan-Meier tidak dapat digambarkan sebagai fungsi (kurva) yang mulus oleh beberapa parameter, misalnya kemiringan atau rasio peluang. Itulah mengapa kita membutuhkan tabel hasil atau grafik.
Metode Kaplan-Meier tidak dapat memodelkan variabel numerik, tetapi hanya kategorik.
Metode Kaplan-Meier tidak dapat memasukkan banyak variabel penjelas. Hal ini buruk, karena membandingkan kelompok dalam hal ketahanan hidup dapat melewatkan efek dari faktor lain, yang dikenal sebagai kovariat atau perancu, yang berpotensi memengaruhi waktu ketahanan hidup kelompok tertentu.
- Rekomendasi**:
Selalu menampilkan ketidakpastian statistik dengan menyertakan 95% CI atau/dan nilai p dari uji Log-Rank. Menampilkan CI pada beberapa titik waktu penting pada plot untuk setiap kelompok perlakuan terkadang lebih jelas daripada menampilkannya untuk semua titik waktu. CI yang tidak tumpang tindih menunjukkan perbedaan yang signifikan antar kelompok.
Pertimbangkan untuk memotong sumbu x. Mengapa? Nah, mata secara alami tertarik ke bagian kanan plot, di mana waktu berakhir. Namun, bagian akhir plot mengandung informasi yang paling sedikit dan ketidakpastian yang paling besar karena jumlah peserta yang tersisa sedikit. Seberapa jauh waktu yang dibutuhkan untuk memperpanjang plot? Terserah Anda.
Selalu tampilkan tabel risiko yang menunjukkan jumlah pasien yang bebas dari kejadian dan masih dalam masa tindak lanjut di setiap kelompok pengobatan pada titik waktu yang relevan.
Analisis Survival 2: Survival model parametrik
Metode non-parametrik Kaplan-Meier (KM) tidak dapat menggambarkan probabilitas kelangsungan hidup dengan fungsi yang mulus, yang berarti tidak dapat memprediksi apa pun. Model parametrik (misalnya Eksponensial, Weibull, dll.) dapat melakukannya! Selain itu, dalam kasus di mana model parametrik sesuai, model parametrik lebih tepat, lebih efektif dan lebih informatif daripada KM atau Cox. Namun, sayangnya, langkah ini sering ditinggalkan karena penggunaan model parametrik di belakang. Dalam tulisan ini kami akan mencoba menutup celah ini.
Mengapa kita memerlukan parametrik survival models
Kerugian utama dari metode Kaplan-Meier (KM) non parametrik yang ditunjukkan pada gambar di atas adalah bahwa metode ini tidak dapat menggambarkan probabilitas kelangsungan hidup dengan fungsi yang mulus, yang berarti tidak dapat memprediksi apa pun. Model parametrik** (misalnya Eksponensial, Weibull, dll.) dapat melakukannya! Selain itu, model parametrik adalah langkah logis dalam perjalanan dari KM ke model semi-parametrik Cox, karena model ini dengan indah menghubungkan titik-titik antara model KM dan Cox dan dengan demikian sangat meningkatkan pemahaman tentang analisis survival. Selain itu, jika model parametrik sesuai, model ini lebih tepat, lebih efektif dan lebih informatif daripada KM atau Cox. Namun, sayangnya, langkah ini sering ditinggalkan karena penggunaan model parametrik di belakang. Dalam tulisan ini kami akan mencoba menutup celah ini.
Apakah Waktu adalah Variabel atau Konstanta?
Jadi, bagaimana kita dapat menggambarkan kelangsungan hidup dengan fungsi mulus? Untuk menjawab pertanyaan ini, pertama-tama mari kita gambarkan fungsi NOT-survival, yang merupakan Bahaya mati.
Saya menggambarkan kematian dengan dua hal: peristiwa kematian itu sendiri dan titik waktu tertentu saat kematian terjadi. Dua hal ini selalu menggambarkan satu peristiwa, karena seseorang hanya mati 💀sekali. Namun, jika beberapa orang mati 💀💀💀💀, mereka tidak akan mati pada saat yang sama, bukan? Tidak. Dengan demikian, waktu kematian akan berubah dan jumlah kematian akan bertambah seiring berjalannya waktu, yang akan membuat waktu itu sendiri menjadi variabel. Begitulah beberapa peristiwa kematian pada titik waktu yang berbeda memungkinkan kita untuk mengekspresikan kematian dalam dua cara:
Melalui jumlah kejadian yang berbeda per unit waktu yang tetap, yang sering disebut sebagai Risiko (Bahaya) untuk mati, atau secara sederhana Bahaya (\(\lambda\) - lambda). Hal ini membuat jumlah kejadian menjadi variabel, dan Hazard menjadi tingkat kematian, misalnya, per hari. Atau,
melalui rentang waktu yang berbeda per jumlah kejadian yang tetap. Interval waktu biasanya diukur sampai kejadian berikutnya terjadi. Hal ini membuat waktu menjadi variabel.
Namun, apa pun yang berubah, jumlah kejadian per satuan waktu, atau waktu per satuan kejadian, perubahan itu sendiri merupakan kunci di sini, dan ada berbagai jenis perubahan.
Perubahan yang stabil dalam bahaya dan kelangsungan hidup
Bayangkan bahwa setiap hari tepat 3 dari 10 orang meninggal di air laut yang dingin setelah kecelakaan Titanic. Ini adalah tingkat kematian yang cukup stabil, atau Bahaya, sebesar 30%. Jika bahaya untuk mati akan terus bertambah dengan tingkat yang sama, maka probabilitas untuk bertahan hidup akan terus berkurang dengan tingkat yang sama. Dengan demikian, Bahaya dan Kelangsungan Hidup dapat dinyatakan dalam istilah satu sama lain. Secara khusus, Bahaya kematian selama t waktu dapat dilihat sebagai F kelangsungan hidup yang gagal (\(F(t)\) dalam rumus kiri di bawah). Atau, S bertahan hidup selama t waktu (\(S(t)\) pada rumus kanan di bawah ini) dapat dilihat sebagai Bahaya tidak mati, atau hanya Bahaya negatif, yang secara matematis dapat diekspresikan sebagai tanda minus “-”, atau “1 -”, di depan Bahaya. Kedua fungsi tersebut menghasilkan garis lurus, di mana Bahaya terus meningkat dan kelangsungan hidup terus menurun (plot di bawah):
\[ F(t) = Hazard * t \]
\[ S(t) = 1 - Hazard * t \]
Peningkatan atau penurunan yang stabil dalam bahaya atau kelangsungan hidup seperti itu tidak alami, tidak realistis. Agak sulit untuk merencanakan dan memantau kematian 😉, kecuali jika Anda seorang pembunuh berantai. Bahaya biasanya meningkat secara eksponensial, misalnya dalam kasus kelaparan, atau menurun secara eksponensial, misalnya dalam kasus pandemi setelah vaksin ditemukan. Oleh karena itu, mari kita lihat kedua perubahan eksponensial tersebut.
Perubahan eksponensial positif dalam bahaya dan kelangsungan hidup
Pikirkan sejenak tentang rasa lapar. Rasa lapar itu terakumulasi, bukan? Dan semakin lama kita menahan lapar, semakin tinggi probabilitas (risiko) kita mati. Plot kiri di bawah ini menunjukkan perkembangan seperti itu, di mana Bahaya kematian \(F(t)\) yang dinyatakan dalam probabilitas kecil pada awalnya (sedikit kematian), tetapi tumbuh secara eksponensial seiring waktu (semakin banyak kematian). Saya suka menyebut tren seperti ini sebagai perubahan eksponensial yang positif, atau semakin cepat. “exp” pada rumus sebelah kiri di bawah ini adalah yang perlu kita tambahkan untuk memplot tren seperti itu.
Sekali lagi, karena kelangsungan hidup dapat dilihat sebagai bahaya negatif, kita dapat mengekspresikan kelangsungan hidup \(S(t)\) hanya dengan menggunakan tanda minus (atau “1-”) di depan Bahaya. Plot di sebelah kanan menampilkan hasil dari fungsi kelangsungan hidup tersebut. Cukup jelas, bahwa jika bahaya kematian rendah di awal “waktu”, maka probabilitas kelangsungan hidupnya tinggi. Pada akhir “waktu”, kelangsungan hidup turun secara eksponensial karena peningkatan eksponensial dari Hazard.
\[ F(t)= exp^{Hazard * t} \]
\[ S(t) = 1 - exp^{Hazard * t} \]
Perubahan eksponensial negatif dalam bahaya dan kelangsungan hidup
Namun, dalam beberapa kasus, lebih banyak orang meninggal di awal “waktu”, di mana setelah itu tingkat kematian menurun seiring berjalannya waktu. Pikirkan tentang pandemi, atau kecelakaan besar. Perubahannya masih bersifat eksponensial, dan menuju ke arah yang sama: naik, dalam kasus bahaya, dan turun, dalam kasus kelangsungan hidup. Namun, perubahan eksponensial tersebut berbelok (berbelok) ke dalam. Saya suka menyebutnya sebagai perubahan eksponensial yang negatif, atau melambat.
“Pembengkokan ke dalam” seperti itu dapat dicapai dengan menyanyikan minus “-” di depan Hazard dan fungsi eksponensial itu sendiri (persamaan kiri dan gambar di bawah). Minus ini tidak menggambarkan bahaya negatif, seperti pada contoh di atas. Ini hanya mengubah kurva dari percepatan menjadi perlambatan. Untuk mendapatkan kelangsungan hidup untuk fungsi ini, kita juga, seperti pada contoh di atas, harus menggunakan tanda “-” (atau “1-”) di depan seluruh ekspresi bahaya eksponensial. Menariknya, dua tanda minus di depan “exp” pada rumus Survival akan saling menetralkan dan menjadi nilai tambah, sehingga menghasilkan \(exp^{-Hazard * t}\) (persamaan kanan dan gambar di bawah).
\[ F(t)= - exp^{-Hazard * t} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \]
\[ S(t) = 1 - ( - exp^{-Hazard * t}) = exp^{-Hazard * t} \]
B(u)ilding Exponential model … finally 🥳
Sekarang kita akhirnya memiliki fungsi kelangsungan hidup mulus yang ingin kita gambarkan di awal. Mari kita lihat sekali lagi dan memplotnya di atas fungsi langkah Kaplan-Meier yang tidak mulus.
Karena fungsi survival \(S(t)\) menunjukkan probabilitas kelangsungan hidup melewati titik waktu tertentu**,
\[ S(t) = P(T > t) \]
dan fungsi eksponensial \(F(t)\) menunjukkan probabilitas gagal bertahan melewati titik waktu tertentu,
\[ F(t) = P(T \leq t) = 1 - exp^{-Hazard * t} \]
fungsi survival kemudian dapat dinyatakan dalam bentuk fungsi bukan survival (eksponensial):
\[ S(t) = P(T > t) = 1 - P(T \leq t) = 1 - F(t) = 1 - [1 - exp^{-Hazard * t}] = exp^{-Hazard * t} \]
Sehingga fungsi survival, yang menunjukkan tingkat penurunan, adalah fungsi hazard yang dibalik, yang menunjukkan tingkat kenaikan. Dengan demikian, kita dapat menulis ulang \(S(t)\) sebagai:
\[S(t) = exp^{- Hazard * t} = exp^{- \lambda* t} \]
di mana Bahaya dan waktu adalah dua parameter yang menggambarkan perubahan eksponensial dalam probabilitas kelangsungan hidup kita, itulah sebabnya model seperti itu dinamakan model eksponensial parametrik. Dan karena Hazard adalah negatif, dan fungsi eksponensial tidak linier (alias. curvy) - model ini menghasilkan kurva eksponensial negatif (lihat di bawah). Laju penurunan yang halus seperti itu menggambarkan probabilitas kelangsungan hidup jauh lebih baik daripada metode Kaplan-Meier, yang secara tiba-tiba (secara bertahap) menurunkan probabilitas hanya setelah suatu kejadian, sambil menjaga probabilitas tetap konstan di antara kejadian-kejadian tersebut.
Bagaimana menghitung model parametrik
library(flexsurv) # for Parametric Survival Modelling
ex <- flexsurvreg(Surv(time, status) ~ 1, data = d2, dist="exponential")
we <- flexsurvreg(Surv(time, status) ~ 1, data = d2, dist="weibull")
ggsurvplot(
ex,
conf.int = FALSE,
surv.median.line = "hv",
xlab = "Days",
ylab = "Survival probability",
break.time.by = 1,
risk.table = F
)
Parameter kurva (\(-\lambda * t\)) memungkinkan kita untuk memodelkan dan memprediksi kelangsungan hidup dan bahaya dari waktu ke waktu, yang merupakan keuntungan utama dari model eksponensial dibandingkan dengan metode Kaplan-Meier, yang tidak bersifat parametrik dan oleh karena itu tidak dapat dimodelkan. Namun, non-linearitas sering kali merepotkan, dan kami lebih suka menggunakan konsep regresi linear, yang meringkas (meregresikan) banyak angka menjadi beberapa angka, seperti intersep (\(\beta_0\)) dan kemiringan (\(\beta_1\)). Untungnya, kurva non-linier dapat dengan mudah “dilinierkan” melalui logaritma natural. Untuk melakukan ini, kita bahkan tidak perlu memahami cara kerja logaritma atau fungsi eksponensial, kita hanya perlu mengetahui bahwa keduanya saling menetralkan satu sama lain. Selain itu, menggunakan “log” (melogaritmakan kedua sisi persamaan di bawah ini) menghasilkan tiga efek samping yang positif:
- Pertama, di sisi kanan persamaan, ini akan memindahkan kurva kita (\(-\lambda * t\)) menjadi sebuah garis (\(b_0 + b_1x_1 + ... + b_kx_k\)), di mana kita akan dapat memiliki intersep dan koefisien \(\beta\) seperti pada regresi logistik linier biasa \[Bahaya = exp^{b_0 + b_1x_1 + ... + b_kx_k}\]
atau
\[ log(Hazard) = b_0 + b_1x_1 + ... + b_kx_k \]
Model survival seperti itu sebenarnya adalah model Poisson. Oleh karena itu, akan sangat membantu jika Anda sudah mengetahui distribusi Poisson. Jika belum, tidak apa-apa, Anda tidak perlu memahami Poisson sebelumnya (kutipan dari video Prof. Marin, lihat referensi).
kedua, ini akan membantu kita untuk terhubung ke model-model selanjutnya, seperti model Weibull dan Cox, karena perbedaan di antara keduanya terutama terletak pada intersep “\(b_0\)”;
dan terakhir, ini akan sangat meningkatkan kemampuan interpretasi, karena Hazard-Ratios (HRs) (yang dihasilkan oleh model eksponensial) dapat diinterpretasikan persis seperti Odds-Ratios (ORs) dari regresi logistik (seperti yang dijelaskan dalam posting saya tentang regresi logistik). Sama halnya dengan Odds dalam regresi logistik, Hazard itu sendiri, yang merupakan probabilitas meninggal SEKARANG, kurang berguna dibandingkan dengan Hazard-Ratio. Hazard-Ratio adalah rasio bahaya seseorang yang terpapar (sakit) terhadap seseorang yang tidak terpapar (sehat). Misalnya jika HR = 2, risiko kematian seseorang yang terpapar adalah dua kali lipat dibandingkan dengan seseorang yang tidak terpapar.
Pikiran akhir
Apakah model parametrik berguna? Tentu saja!
Tiga kurva di atas yang menggambarkan berbagai jenis perubahan dalam kelangsungan hidup dan bahaya dari waktu ke waktu dimaksudkan untuk mengatakan bahwa distribusi bisa sangat berbeda. Jika distribusi yang sesuai dapat ditemukan, model parametrik lebih informatif daripada model KM atau Cox
model ini dapat memprediksi probabilitas ketahanan hidup pada titik waktu tertentu, hazard kejadian, rata-rata dan median waktu ketahanan hidup sudah tersedia
model ini juga sedikit lebih efisien dan menghasilkan estimasi yang lebih tepat karena kecocokan yang lebih baik (lihat gambar di atas)
Apakah model parametrik itu sempurna? Tentu saja tidak!
model parametrik perlu menentukan distribusi, yang mungkin sulit untuk diidentifikasi
model ini juga secara matematis lebih kompleks, misalnya KM, dan oleh karena itu lebih jarang digunakan, sehingga mengurangi perbandingan hasil antar studi. Karena popularitasnya yang lebih rendah, saya tidak akan membahas lebih dalam tentang model parametrik tertentu: Model Weibull, Gompertz, model Waktu Kegagalan yang Dipercepat, dll.
Jadi, terlepas dari kenyataan bahwa model parametrik adalah alternatif yang baik untuk model regresi KM dan Cox (yang tidak perlu menentukan distribusi apa pun), KM dan Cox tetap menjadi metode yang paling populer untuk menganalisis data survival. Dan itulah mengapa langkah logis berikutnya dalam perjalanan statistik Anda adalah mempelajari Model Cox Proportional Hazard (sedang berlangsung).
Reference
Kaplan, E.L. and Meier, P. (1958) Nonparametric Estimation from Incomplete Observations. Journal of the American Statistical Association, 53, 457-481. http://dx.doi.org/10.1080/01621459.1958.10501452↩︎