Sensoring terjadi ketika kita mempunyai informasi mengenai waktu survival individu, namun kita tidak mengetahui dengan pasti waktu survivalnya (kejadian yang diamati tidak terjadi). Sensoring dapat dianggap sebagai data hilang (missing data) dalam analisis daya tahan. Beberapa penyebab data tersensor adalah:
Subyek tidak mengalami kejadian sampai akhir studi
Subyek tidak dapat di follow-up (pindah; tidak hadir saat pengumpulan data; tidak dapat dikontak)
Subyek mengundurkan diri dari studi
Subyek meninggal (disebabkan oleh hal lain yang bukan merupakan interest dalam kajian)
Terdapat beberapa jenis sensoring, yaitu: sensoring kanan (tetap dan acak), sensoring kiri, dan sensoring interval. Kita akan fokus pada jenis sensoring kanan (tetap ataupun acak).
Fungsi daya tahan menunjukkan peluang subjek bertahan setelah waktu tertentu. Secara teoritis, fungsi daya tahan didefinisikan sebagai berikut:
\(S(t)=P(T>t)=1-P(T<=t)=1-F(t)\)
dengan \(F(t)\) adalah peluang kumulatif.
Untuk melakukan analisis daya tahan, kita harus memiliki data terkait waktu bertahan subjek (\(Y_i\)) dan indikator kejadian (apakah subjek tersensor atau mengalami kejadian).
Pada ilustrasi ini akan kita gunakan data pasien kanker paru-paru yang berasal dari North Central Cancer Treatment Group. Data memuat peubah sebagai berikut:
nst: Institution code
time: Survival time in days
status: censoring status 1=censored, 2=dead
age: Age in years
sex: Male=1 Female=2
ph.ecog: ECOG performance score as rated by the physician. (0=asymptomatic, 1= symptomatic but completely ambulatory, 2= in bed <50% of the day, 3= in bed > 50% of the day but not bedbound, 4 = bedbound)
ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
pat.karno: Karnofsky performance score as rated by patient
meal.cal: Calories consumed at meals
wt.loss: Weight loss in last six months (pounds)
library(dplyr)
library(survival)
library(ggsurvfit)
data("cancer")
str(cancer)
## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 2 2 1 2 2 1 2 2 2 2 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...
Struktur data menunjukkan bahwa peubah status memiliki
nilai 1 dan 2, yang mana nilai tersebut tidak dalam bentuk baku. Fungsi
Surv() secara default menerima nilai dari peubah indikator
berupa nilai True/False atau 1/0, dengan 1
berarti kejadian dan 0 berarti tersensor. Sehingga kita perlu mengubah
nilai dari peubah status ke dalam bentuk yang dapat
diterima oleh fungsi Surv().
cancer_df=cancer%>%
mutate(status=recode(status,'1'=0, '2'=1))
head(cancer_df)
Peubah yang akan digunakan adalah time dan
status.
cancer_df2=cancer_df%>%
select(time,status)
head(cancer_df2)
Untuk melakukan survival analisis, langkah awal kita membuat objek survival.
surv_obj=Surv(cancer_df2$time,cancer_df2$status)
head(surv_obj)
## [1] 306 455 1010+ 210 883 1022+
Perhatikan hasil dari objek survival di atas. Nilai data ada yang hanya berupa angka, dan ada yang berupa angka diikuti tanda plus (+). Nilai yang diikuti tanda \(+\) merupakan subjek yang tersensor (0), sementara nilai yang tidak diikuti tanda \(+\) adalah responden yang mengalami kejadian (1).
Selanjutnya kita lakukan analisis daya tahan menggunakan metode Kaplan Meier.
fit=survfit(surv_obj~1)
summary(fit, times=head(fit$time))
## Call: survfit(formula = surv_obj ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 228 1 0.996 0.00438 0.987 1.000
## 11 227 3 0.982 0.00869 0.966 1.000
## 12 224 1 0.978 0.00970 0.959 0.997
## 13 223 2 0.969 0.01142 0.947 0.992
## 15 221 1 0.965 0.01219 0.941 0.989
## 26 220 1 0.961 0.01290 0.936 0.986
Berdasarkan hasil di atas, interpretasi untuk baris pertama adalah untuk waktu 5 hari, dari 228 orang ada 1 orang yang mengalami kejadian, sehingga peluang seseorang yang terkena kanker paru-paru dan bertahan setelah 5 hari adalah 0.996 atau 99.6%. Nilai peluang ini dapat kita kitung sebagai berikut:
\(S(t)=1-\frac{1}{228}=0.996\)
Sementara untuk baris kedua, yaitu untuk waktu 11 hari, dari 227 orang ada 3 orang yang mengalami kejadian, sehingga peluang seseorang yang terkena kanker paru-paru dan bertahan setelah 11 hari adalah 0.982 atau 98.2%. Nilai peluang ini dapat kita kitung sebagai berikut:
\(S(t)=(1-\frac{1}{228})(1-\frac{3}{227})=0.9956\times0.9868=0.982\)
Hal ini analog untuk baris keenam, yaitu untuk waktu 26 hari, dari 220 orang ada 1 orang yang mengalami kejadian, sehingga peluang seseorang yang terkena kanker paru-paru dan bertahan setelah 26 hari adalah 0.961 atau 96.1%. Nilai peluang ini dapat kita kitung sebagai berikut:
\(S(t)=(1-\frac{1}{228})(1-\frac{3}{227})...(1-\frac{1}{220})=0.9956\times0.9868\times...\times0.9955=0.961\)
Untuk mempermudah melihat peluang bertahan hidup secara keseluruhan, kita buat plot untuk hasil analisis Kaplan Meier.
plot(fit, xlab='Hari', ylab='Peluang Survive')
Selain menggunakan fungsi plot, plot dari analisis Kaplan Meier dapat
kita buat menggunakan fungsi survfit.
survfit2(surv_obj~1)%>%
ggsurvfit()+
labs(x='Hari',
y='Peluang Survive')+
add_confidence_interval()+
add_risktable()
Grafik tersebut menunjukkan bahwa peluang seseorang untuk bertahan semakin menurun seiring bertambahnya waktu (\(t\) membesar). Di bawah plot terdapat tabel risiko yang berisi angka banyaknya individu yang berisiko (angka di atas) dan banyaknya kejadian kumulatif (angka di bawah).
Selanjutnya kita akan menduga peluang seseorang akan bertahan selama \(x\) tahun.
# menduga seseorang bertahan selama 1 tahun
summary(fit, times=365)
## Call: survfit(formula = surv_obj ~ 1)
##
## 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
Berdasarkan hasil kita peroleh bahwa peluang seseorang bertahan selama 1 tahun (365 hari) adalah 0.41 atau 41%.
## menduga seseorang bertahan selama 2 tahun
summary(fit, times=730)
## Call: survfit(formula = surv_obj ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 730 13 159 0.116 0.0283 0.0716 0.187
Sedangkan peluang seseorang bertahan selama 2 tahun (730 hari) adalah 0.116 atau 11.6%.
Selanjutnya untuk melihat rata-rata waktu bertahan pasien kanker paru-paru, kita gunakan nilai mediannya.
fit
## Call: survfit(formula = surv_obj ~ 1)
##
## n events median 0.95LCL 0.95UCL
## [1,] 228 165 310 285 363
Rata-rata waktu bertahan pasien kanker paru-paru adalah selama 310 hari.