Fungsi Daya Tahan
Laily Nissa Atul Mualifah
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.
- Jika dalam pengamatan tidak ada data tersensor, penduga empiris dari fungsi daya tahan dapat dihitung dari proporsi individu yang masih bertahan pada waktu \(t\), yaitu:
\(\hat S(t)=\frac{individu\;yang\;bertahan\;di\;waktu\;t}{total\;individu}\)
- Jika dalam pengamatan ada data tersensor, penduga empiris dari fungsi daya tahan dapat dihitung menggunakan pendekatan nonparametrik (Kaplan eier Life Table). Untuk Kaplan Meier fungsi daya tahan dihitung dengan formula sebagai berikut:
\(\hat S(t)=\Pi(1-\frac{d_i}{n_i})\)
\(d_i\): banyaknya kejadian pada waktu t
\(n_i\): banyaknya individu pada waktu t
Ilustrasi (Kaplan Meier)
Pada ilustrasi ini akan kita gunakan data percobaan acak terapi untuk menghentikan kebiasaan merokok. Data memuat peubah sebagai berikut:
id: id responden
ttr: waktu (hari) sampai individu kembali merokok
relapse: indikator kembali merokok (0: Tidak/tersensor, 1:Ya/kejadian)
grp: jenis terapi yang diberikan (combination atau patchOnly)
age: usia (tahun) saat percobaan
gender: jenis kelamin (Female atau Male)
race: ras (black, hispanic, white atau other)
employment: pekerjaan (ft: full-time, pt: part-time, other)
yearsSmoking: lama merokok (tahun)
levelSmoking: kategori perokok (heavy atau light)
ageGroup2: kelompok umur dengan kategori 21-49 atau 50+
ageGroup4: kelompok umur dengan kategori 21-34, 35-49, 50-64 atau 65+
priorAttempts: banyaknya percobaan yang telah dilakukan untuk berhenti merokok
longestNoSmoke: jangka waktu terlama (hari) yang pernah dilakukan tanpa merokok
library(asaur)
library(dplyr)
library(survival)
library(survminer)
library(ggsurvfit)data("pharmacoSmoking")
str(pharmacoSmoking)## 'data.frame': 125 obs. of 14 variables:
## $ id : int 21 113 39 80 87 29 16 35 54 70 ...
## $ ttr : int 182 14 5 16 0 182 14 77 2 0 ...
## $ relapse : int 0 1 1 1 1 0 1 1 1 1 ...
## $ grp : Factor w/ 2 levels "combination",..: 2 2 1 1 1 1 2 2 2 2 ...
## $ age : int 36 41 25 54 45 43 66 78 40 38 ...
## $ gender : Factor w/ 2 levels "Female","Male": 2 2 1 2 2 2 2 1 1 2 ...
## $ race : Factor w/ 4 levels "black","hispanic",..: 4 4 4 4 4 2 1 1 1 1 ...
## $ employment : Factor w/ 3 levels "ft","other","pt": 1 2 2 1 2 1 3 2 1 1 ...
## $ yearsSmoking : int 26 27 12 39 30 30 54 56 25 23 ...
## $ levelSmoking : Factor w/ 2 levels "heavy","light": 1 1 1 1 1 1 1 2 1 2 ...
## $ ageGroup2 : Factor w/ 2 levels "21-49","50+": 1 1 1 2 1 1 2 2 1 1 ...
## $ ageGroup4 : Factor w/ 4 levels "21-34","35-49",..: 2 2 1 3 2 2 4 4 2 2 ...
## $ priorAttempts : int 0 3 3 0 0 2 0 10 4 10 ...
## $ longestNoSmoke: int 0 90 21 0 0 1825 0 15 7 90 ...Pertama kita akan melihat perbedaan peluang daya tahan perokok berdasarkan jenis terapi yang diberikan.
surv_obj=Surv(pharmacoSmoking$ttr,pharmacoSmoking$relapse)
fit=survfit(surv_obj~grp, data=pharmacoSmoking)
summary(fit, times=head(fit$time,10))## Call: survfit(formula = surv_obj ~ grp, data = pharmacoSmoking)
##
## grp=combination
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 61 4 0.934 0.0317 0.874 0.999
## 2 57 3 0.885 0.0408 0.809 0.969
## 4 54 1 0.869 0.0432 0.788 0.958
## 5 53 2 0.836 0.0474 0.748 0.934
## 8 51 2 0.803 0.0509 0.709 0.909
## 10 49 1 0.787 0.0524 0.691 0.897
## 12 48 1 0.770 0.0538 0.672 0.884
## 14 47 1 0.754 0.0551 0.653 0.870
## 15 46 2 0.721 0.0574 0.617 0.843
## 16 44 1 0.705 0.0584 0.599 0.829
##
## grp=patchOnly
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 64 8 0.875 0.0413 0.798 0.960
## 2 51 8 0.750 0.0541 0.651 0.864
## 4 47 3 0.703 0.0571 0.600 0.824
## 5 45 0 0.703 0.0571 0.600 0.824
## 8 43 3 0.656 0.0594 0.550 0.784
## 10 42 0 0.656 0.0594 0.550 0.784
## 12 42 1 0.641 0.0600 0.533 0.770
## 14 41 6 0.547 0.0622 0.438 0.684
## 15 35 2 0.516 0.0625 0.407 0.654
## 16 33 0 0.516 0.0625 0.407 0.654Berdasarkan hasil di atas, perokok yang memperoleh terapi combination memiliki daya tahan yang lebih baik dibanding yang memperoleh terapi patchOnly, hal ini dapat dilihat dari nilai peluang bertahan (survival) dari kedua kelompok. Misalnya untuk waktu 14 hari, dari kelompok terapi combination peluang bertahan untuk tidak merokok setelah 14 hari adalah 75.4% sementara dari kelompok terapi patchOnly peluang bertahan untuk tidak merokok setelah 14 hari 54.7%. Sehingga dapat dikatakan bahwa terapi combination lebih efektif untuk digunakan.
Plot Kaplan Meier untuk hasil di atas adalah sebagai berikut.
ggsurvplot(fit,
conf.int = T,
risk.table = T,
risk.table.col='strata',
surv.median.line = 'hv',
ggtheme = theme_bw())
Garis hitam pada plot merupakan nilai median waktu bertahan untuk kedua jenis terapi, yaitu 65 hari untuk terapi combination dan 23 hari untuk terapi patchOnly.
fit## Call: survfit(formula = surv_obj ~ grp, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## grp=combination 61 37 65 50 NA
## grp=patchOnly 64 52 23 14 56Kita akan memperbaiki tampilan plot di atas dengan menambahkan info banyaknya pengamatan tersensor.
ggsurvplot(fit,
conf.int = T,
conf.int.style='step',
break.time.by=25,
xlab='Waktu (hari)',
xlim=c(0,190),
fontsize=3,
risk.table = 'abs_pct',
risk.table.y.text.col=T,
risk.table.y.text=F,
ncensor.plot=T,
surv.median.line = 'hv',
ggtheme = theme_light())
Dari tabel banyaknya pengamatan tersensor tampak bahwa untuk kedua jenis terapi tidak ada pengamatan tersensor yang waktu bertahannya kurang dari masa studi (182 hari), sehingga dapat kita duga bahwa pengamatan-pengamatan yang tersensor adalah tersensor kanan (sampai akhir masa studi tidak merokok kembali).
Selanjutnya kita lihat perbedaan peluang daya tahan perokok berdasarkan semua peubah kategorik yang ada dalam data, yaitu gender, race, employment, levelsmoking, dan agegroup4
grup=c('gender','race','employment','levelSmoking','ageGroup4')
formula=sapply(grup,
function(x) as.formula(paste('Surv(ttr,relapse)~',x)))
models=lapply(formula, function(x) survfit(x, data=pharmacoSmoking))
models## $gender
## Call: survfit(formula = x, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## gender=Female 81 61 49 21 75
## gender=Male 44 28 40 14 NA
##
## $race
## Call: survfit(formula = x, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## race=black 38 29 22.5 12 77
## race=hispanic 8 5 94.5 2 NA
## race=other 2 1 80.0 80 NA
## race=white 77 54 56.0 25 140
##
## $employment
## Call: survfit(formula = x, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## employment=ft 72 49 56 30 140
## employment=other 39 28 25 12 110
## employment=pt 14 12 18 14 NA
##
## $levelSmoking
## Call: survfit(formula = x, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## levelSmoking=heavy 89 62 49.0 15 80
## levelSmoking=light 36 27 50.5 28 140
##
## $ageGroup4
## Call: survfit(formula = x, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## ageGroup4=21-34 16 14 53 5 140
## ageGroup4=35-49 50 42 18 12 63
## ageGroup4=50-64 48 25 122 49 NA
## ageGroup4=65+ 11 8 45 20 NAsplot=lapply(models, function(x) {ggsurvfit(x)+add_confidence_interval(type='lines')+
add_risktable(risktable_stats ='n.risk',
risktable_group='risktable_stats',
size=3)+
add_quantile()})
splot## $gender
##
## $race
##
## $employment
##
## $levelSmoking
##
## $ageGroup4
Dari hasil di atas, untuk peubah jenis kelamin, perempuan (female) memiliki median waktu bertahan 49 hari, sementara laki-laki (male) 40 hari. Jika kita perhatikan plot peluang daya tahan untuk jenis kelamin di atas, peluang daya tahan perempuan sebelum 60 hari relatif sama dengan laki-laki, akan tetapi setelah 60 hari peluang daya tahan perempuan mengecil sedangkan peluang daya tahan laki-laki cukup konstan. Peluang daya tahan perempuan dan laki-laki setelah 140 hari berturut-turut adalah 27.2% dan 38.6%.
Sementara untuk jenis pekerjaan (employment), tampak bahwa pekerja penuh waktu (fulltime) memiliki peluang daya tahan yang lebih besar dibanding paruh waktu (parttime) atau lainnya (other), sedangkan pekerja paruh waktu memiliki peluang daya tahan paling kecil,dengan median waktu daya tahan untuk fulltime, parttime dan other berturut-turut 56 hari, 18 hari, dan 25 hari.
Jika kita lihat berdasarkan kategori perokok, antara perokok berat dan perokok ringan, peluang daya tahan kedua tidak berbeda secara signifikan.
Interpretasi untuk peubah kategorik yang lain silakan anda interpretasikan.
Selanjutnya kita akan mencoba melihat peluang daya tahan jika kita gunakan beberapa peubah kategorik, yang dalam hal ini kita gunakan peubah jenis kelamin dan pekerjaan.
fit2<-survfit(surv_obj~grp+employment, data=pharmacoSmoking)
fit2## Call: survfit(formula = surv_obj ~ grp + employment, data = pharmacoSmoking)
##
## n events median 0.95LCL 0.95UCL
## grp=combination, employment=ft 34 21 64.0 50 NA
## grp=combination, employment=other 20 11 125.0 20 NA
## grp=combination, employment=pt 7 5 63.0 8 NA
## grp=patchOnly, employment=ft 38 28 42.5 14 140
## grp=patchOnly, employment=other 19 17 8.0 3 77
## grp=patchOnly, employment=pt 7 7 15.0 1 NAggsurv=ggsurvplot(fit2,
conf.int = T,
conf.int.style='step',
surv.median.line = 'hv',
ggtheme = theme_bw())
ggsurv$plot+
theme_bw()+
theme(legend.position = 'right')+
facet_grid(rows =vars(grp))
Untuk jenis terapi combination, peluang daya tahan pekerja fulltime dan other relatif sama, tetapi pekerja other memiliki median waktu daya tahan yang lebih besar yaitu 125 hari, sedangkan median waktu bertahan antara pekerja fulltime atau parttime hampir sama yaitu 64 dan 63 hari. Sementara untuk jenis terapi patchOnly, pekerja fulltime memiliki peluang daya tahan yang paling besar dibanding parttime dan other dengan median waktu bertahan berturut-turut 43 hari, 25 hari, dan 8 hari.
Untuk model yang lebih kompleks lagi kita gunakan 3 peubah yaitu jenis kelamin, pekerjaan, dan jenis terapi
fit3<-survfit(surv_obj~grp+employment+gender, data=pharmacoSmoking)
fit3## Call: survfit(formula = surv_obj ~ grp + employment + gender, data = pharmacoSmoking)
##
## n events median 0.95LCL
## grp=combination, employment=ft , gender=Female 19 12 60.0 42
## grp=combination, employment=ft , gender=Male 15 9 140.0 16
## grp=combination, employment=other, gender=Female 14 7 110.0 20
## grp=combination, employment=other, gender=Male 6 4 85.0 0
## grp=combination, employment=pt , gender=Female 6 4 69.0 14
## grp=combination, employment=pt , gender=Male 1 1 8.0 NA
## grp=patchOnly, employment=ft , gender=Female 26 23 34.0 14
## grp=patchOnly, employment=ft , gender=Male 12 5 NA 6
## grp=patchOnly, employment=other, gender=Female 10 9 7.5 2
## grp=patchOnly, employment=other, gender=Male 9 8 14.0 2
## grp=patchOnly, employment=pt , gender=Female 6 6 18.0 1
## grp=patchOnly, employment=pt , gender=Male 1 1 14.0 NA
## 0.95UCL
## grp=combination, employment=ft , gender=Female NA
## grp=combination, employment=ft , gender=Male NA
## grp=combination, employment=other, gender=Female NA
## grp=combination, employment=other, gender=Male NA
## grp=combination, employment=pt , gender=Female NA
## grp=combination, employment=pt , gender=Male NA
## grp=patchOnly, employment=ft , gender=Female 84
## grp=patchOnly, employment=ft , gender=Male NA
## grp=patchOnly, employment=other, gender=Female NA
## grp=patchOnly, employment=other, gender=Male NA
## grp=patchOnly, employment=pt , gender=Female NA
## grp=patchOnly, employment=pt , gender=Male NAggsurv=ggsurvplot(fit3,
conf.int = T,
conf.int.style='step',
ggtheme = theme_bw())
ggsurv$plot+
theme_bw()+
theme(legend.position = 'bottom',
legend.text = element_text(size=4),
legend.title = element_text(size=4))+
facet_grid(cols =vars(employment),
rows=vars(gender))
Silakan anda interpretasikan dan coba untuk peubah yang lain.