Dalam fungsi daya tahan (survivor function) kita dapat membandingkan daya tahan antar kelompok Sebelumnya (dalam Kaplan Meier), kita dapat membandingkan peluang daya tahan antar kelompok yang berbeda berdasarkan kurva daya tahannya. Selanjutnya kita perlu melakukan pengujian untuk mengetahui apakah daya tahan antar kelompok tersebut berbeda secara signifikan.

Uji Log-Rank merupakan uji statistik untuk membandingkan kurva daya tahan dari dua kelompok atau lebih. Hipotesis yang digunakan dalam Uji Log-Rank adalah:

\(H_0:tidak\;ada\;perbedaan\;daya\;tahan\;antar\;kelompok\)

\(H_1:terdapat\;perbedaan\;daya\;tahan\;antar\;kelompok\)

Statistik uji Log-Rank mengikuti sebaran Khi-Kuadrat dengan derajat bebas \(k-1\), dimana \(k\) adalah banyaknya kelompok yang dibandingkan.

\(\chi^2=\sum_{i=1}^m\frac{(\sum O_{it}-\sum E_{it})^2}{\sum E_{it}}\)

dengan,

\(E_{it}=\frac{n_{it}\times O_t}{n_t}\)

dimana,

\(O_{it}\): banyaknya kejadian pada kelompok i waktu t

\(E_{it}\): banyaknya kejadian yang diharapkan pada kelompok i waktu t

\(n_{it}\): banyaknya individu berisiko pada kelompok i

\(O_t\): total banyaknya kejadian waktu t

\(n_t\): total banyaknya individu berisiko waktu t

Ilustrasi

Pada ilustrasi ini akan kita gunakan data percobaan acak terapi untuk menghentikan kebiasaan merokok (data yang sama dengan ilustrasi Kaplan Meier). Kejadian yang diamati adalah perokok kembali merokok.

library(asaur)
library(dplyr)
library(survival)
library(survminer)
library(ggsurvfit)
data("pharmacoSmoking")
glimpse(pharmacoSmoking)
## Rows: 125
## Columns: 14
## $ id             <int> 21, 113, 39, 80, 87, 29, 16, 35, 54, 70, 84, 85, 25, 47…
## $ ttr            <int> 182, 14, 5, 16, 0, 182, 14, 77, 2, 0, 12, 182, 21, 3, 1…
## $ relapse        <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1…
## $ grp            <fct> patchOnly, patchOnly, combination, combination, combina…
## $ age            <int> 36, 41, 25, 54, 45, 43, 66, 78, 40, 38, 64, 51, 37, 65,…
## $ gender         <fct> Male, Male, Female, Male, Male, Male, Male, Female, Fem…
## $ race           <fct> white, white, white, white, white, hispanic, black, bla…
## $ employment     <fct> ft, other, other, ft, other, ft, pt, other, ft, ft, oth…
## $ yearsSmoking   <int> 26, 27, 12, 39, 30, 30, 54, 56, 25, 23, 30, 35, 23, 50,…
## $ levelSmoking   <fct> heavy, heavy, heavy, heavy, heavy, heavy, heavy, light,…
## $ ageGroup2      <fct> 21-49, 21-49, 21-49, 50+, 21-49, 21-49, 50+, 50+, 21-49…
## $ ageGroup4      <fct> 35-49, 35-49, 21-34, 50-64, 35-49, 35-49, 65+, 65+, 35-…
## $ priorAttempts  <int> 0, 3, 3, 0, 0, 2, 0, 10, 4, 10, 12, 1, 5, 6, 5, 2, 1, 1…
## $ longestNoSmoke <int> 0, 90, 21, 0, 0, 1825, 0, 15, 7, 90, 365, 7, 1095, 180,…

Pertama kita akan membandingkan peluang daya tahan perokok berdasarkan jenis terapi yang diberikan (combination dan patchOnly). Hipotesis yang kita gunakan adalah:

\(H_0:\) daya tahan perokok yang mendapat terapi combination dan patchOnly sama

\(H_1:\) daya tahan perokok yang mendapat terapi combination dan patchOnly tidak sama

uji_logrank=survdiff(Surv(ttr,relapse)~grp, data=pharmacoSmoking)
uji_logrank
## Call:
## survdiff(formula = Surv(ttr, relapse) ~ grp, data = pharmacoSmoking)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## grp=combination 61       37     49.9      3.36      8.03
## grp=patchOnly   64       52     39.1      4.29      8.03
## 
##  Chisq= 8  on 1 degrees of freedom, p= 0.005

Berdasarkan hasil di atas diperoleh statistik uji \(\chi^2=8>\chi^2_{tabel}=3.84\) dengan \(p-value=0.005\), maka kita tolak \(H_0\). Jadi pada taraf nyata 5% kita memiliki cukup bukti untuk mengatakan bahwa daya tahan perokok yang mendapat terapi combination dan patchOnly tidak sama.

Berikut plot kurva daya tahan untuk kedua jenis terapi beserta p-value uji Log-Rank.

fit=survfit(Surv(ttr,relapse)~grp, data=pharmacoSmoking)
ggsurvplot(fit,
           pval = T,
           conf.int = T,
           conf.int.style='step',
           break.time.by=25,
           xlab='Waktu (hari)',
           xlim=c(0,190),
           surv.median.line = 'hv',
           ggtheme = theme_bw())

Berdasarkan plot di atas perokok yang memperoleh terapi combination memiliki daya tahan yang lebih tinggi dengan median 65 hari dibanding perokok dengan terapi patchOnly dengan median 23 hari.

Kedua kita akan membandingkan peluang daya tahan perokok berdasarkan kelompok umur (21-34, 35-49, 50-64 dan 65+).

uji_logrank2=survdiff(Surv(ttr,relapse)~ageGroup4, data=pharmacoSmoking)
uji_logrank2
## Call:
## survdiff(formula = Surv(ttr, relapse) ~ ageGroup4, data = pharmacoSmoking)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## ageGroup4=21-34 16       14    10.43   1.22342  1.45e+00
## ageGroup4=35-49 50       42    30.46   4.37639  6.95e+00
## ageGroup4=50-64 48       25    40.07   5.66881  1.09e+01
## ageGroup4=65+   11        8     8.04   0.00025  2.86e-04
## 
##  Chisq= 11.9  on 3 degrees of freedom, p= 0.008

Berdasarkan hasil di atas diperoleh statistik uji \(\chi^2=11.9>\chi^2_{tabel}=7.82\) dengan \(p-value=0.005\), maka kita tolak \(H_0\). Jadi pada taraf nyata 5% kita memiliki cukup bukti untuk mengatakan bahwa terdapat perbedaan daya tahan perokok antar kelompok umur.

Selanjutnya kita lakukan uji lanjut/ perbandingan berganda untuk mengetahui kelompok umur mana saja yang memiliki daya tahan berbeda.

res_lanjut=pairwise_survdiff(Surv(ttr,relapse)~ageGroup4, 
                  data=pharmacoSmoking, p.adjust.method='hommel')
res_lanjut
## 
##  Pairwise comparisons using Log-Rank test 
## 
## data:  pharmacoSmoking and ageGroup4 
## 
##       21-34  35-49  50-64 
## 35-49 0.9967 -      -     
## 50-64 0.1018 0.0078 -     
## 65+   0.9967 0.7987 0.7987
## 
## P value adjustment method: hommel

Nilai-nilai yang diperoleh merupakan p-value untuk perbandingan setiap 2 kelompok berbeda. Misalnya, baris pertama perbandingan antara kelompok umur 21-34 dan 35-49 diperoleh p-value 0.9967, sehingga dapat disimpulkan bahwa daya tahan perokok antara kelompok umur 21-34 dan 35-49 tidak berbeda secara signifikan. Satu-satunya perbandingan yang signifikan adalah perbandingan antara kelompok umur 35-49 dan 50-64 yang mana diperoleh p-value 0.0078.

Untuk mempermudah dalam melihat perbandingan mana saja yang signifikan, kita ubah p-value yang diperoleh dengan simbol sesuai yang kita inginkan.

symnum(res_lanjut$p.value, cutpoints = c(0, 0.0001, 0.001, 
                                         0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")
##       21-34 35-49 50-64
## 35-49                  
## 50-64       **         
## 65+                    
## attr(,"legend")
## [1] 0 '****' 1e-04 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1 \t    ## NA: ''

Terlihat bahwa perbandingan antara kelompok umur 35-49 dan 50-64 signifikan pada taraf nyata 5% (0.05).

Berikut plot kurva daya tahan berdasarkan kelompok umur beserta p-value uji Log-Rank.

fit2=survfit(Surv(ttr,relapse)~ageGroup4, data=pharmacoSmoking)
ggsurvplot(fit2,
           pval = T,
           conf.int = T,
           conf.int.style='step',
           break.time.by=25,
           xlab='Waktu (hari)',
           xlim=c(0,190),
           surv.median.line = 'hv',
           ggtheme = theme_bw())

Berdasarkan plot di atas perokok pada kelompok usia 50-64 memiliki daya tahan yang lebih tinggi dengan median 122 hari dibanding perokok pada usia 35-49 dengan median 18 hari.

Terakhir kita bandingkan peluang daya tahan perokok berdasarkan jenis terapi dan jenis pekerjaan (full-time, part-time dan lainnya).

uji_logrank3=survdiff(Surv(ttr,relapse)~grp+employment, 
                      data=pharmacoSmoking)
uji_logrank3
## Call:
## survdiff(formula = Surv(ttr, relapse) ~ grp + employment, data = pharmacoSmoking)
## 
##                                    N Observed Expected (O-E)^2/E (O-E)^2/V
## grp=combination, employment=ft    34       21    28.07    1.7804   2.71076
## grp=combination, employment=other 20       11    16.77    1.9880   2.55117
## grp=combination, employment=pt     7        5     5.10    0.0021   0.00231
## grp=patchOnly, employment=ft      38       28    26.62    0.0714   0.10596
## grp=patchOnly, employment=other   19       17     8.93    7.2827   8.53696
## grp=patchOnly, employment=pt       7        7     3.50    3.5067   3.83448
## 
##  Chisq= 15.6  on 5 degrees of freedom, p= 0.008

Berdasarkan hasil di atas diperoleh statistik uji \(\chi^2=15.6>\chi^2_{tabel}=11.07\) dengan \(p-value=0.005\), maka kita tolak \(H_0\). Jadi pada taraf nyata 5% kita memiliki cukup bukti untuk mengatakan bahwa terdapat perbedaan daya tahan perokok berdasarkan jenis terapi dan jenis pekerjaan.

res_lanjut2=pairwise_survdiff(Surv(ttr,relapse)~grp+employment, 
                  data=pharmacoSmoking, p.adjust.method='hommel')
symnum(res_lanjut2$p.value, cutpoints = c(0, 0.0001, 0.001, 
                                         0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")
##                                   grp=combination, employment=ft   
## grp=combination, employment=other                                  
## grp=combination, employment=pt                                     
## grp=patchOnly, employment=ft                                       
## grp=patchOnly, employment=other   *                                
## grp=patchOnly, employment=pt                                       
##                                   grp=combination, employment=other
## grp=combination, employment=other                                  
## grp=combination, employment=pt                                     
## grp=patchOnly, employment=ft                                       
## grp=patchOnly, employment=other   +                                
## grp=patchOnly, employment=pt                                       
##                                   grp=combination, employment=pt   
## grp=combination, employment=other                                  
## grp=combination, employment=pt                                     
## grp=patchOnly, employment=ft                                       
## grp=patchOnly, employment=other                                    
## grp=patchOnly, employment=pt                                       
##                                   grp=patchOnly, employment=ft   
## grp=combination, employment=other                                
## grp=combination, employment=pt                                   
## grp=patchOnly, employment=ft                                     
## grp=patchOnly, employment=other                                  
## grp=patchOnly, employment=pt                                     
##                                   grp=patchOnly, employment=other
## grp=combination, employment=other                                
## grp=combination, employment=pt                                   
## grp=patchOnly, employment=ft                                     
## grp=patchOnly, employment=other                                  
## grp=patchOnly, employment=pt                                     
## attr(,"legend")
## [1] 0 '****' 1e-04 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1 \t    ## NA: ''

Pada taraf nyata 1% (0.01) daya tahan perokok yang memperoleh terapi combination dan bekerja full-time berbeda dengan perokok yang memperoleh terapi patchOnly dan bekerja lainnya. Sementara pada taraf nyata 10% (0.1) daya tahan perokok yang memperoleh terapi combination dan bekerja lainnya berbeda dengan perokok yang memperoleh terapi patchOnly dan bekerja lainnya. Sedangkan untuk kombinasi jenis terapi dan pekerjaan yang lain, daya tahannya tidak berbeda secara signifikan.

fit3=survfit(Surv(ttr,relapse)~grp+employment, data=pharmacoSmoking)
ggplot=ggsurvplot(fit3,
                  pval = T,
                  conf.int = T,
                  conf.int.style='step',
                  break.time.by=25,
                  xlab='Waktu (hari)',
                  xlim=c(0,190),
                  surv.median.line = 'hv',
                  ggtheme = theme_bw())
ggplot$plot+
  theme(legend.position = 'bottom',
        legend.text = element_text(size=7),
        legend.title = element_text(size=7))