Berikut adalah fungsi daya tahan dan fungsi hazard untuk beberapa sebaran:

Distribusi \(S(t)\) \(h(t)\)
Eksponensial \(exp(-\lambda t)\) \(\lambda\)
Weibull \(exp(-\lambda t^p)\) \(\lambda p t^{p-1}\)
Log-logistik \(\frac{1}{1+\lambda t^p}\) \(\frac{\lambda p t^{p-1}}{1+\lambda t^p}\)

Dalam sebaran Weibull parameter \(\lambda\) merupakan reparameterisasi dari kovariat dan \(p\) adalah shape parameter yang menentukan bentuk dari fungsi hazard. Jika \(p>1\) maka fungsi hazard meningkat seiring meningkatnya waktu, jika \(p=1\) maka fungsi hazard konstan sehingga menjadi sebaran eksponensial, jika \(p<1\) maka fungsi hazard menurun seiring meningkatnya waktu. Jika asumsi AFT dalam sebaran Weibull terpenuhi maka asumsi PH juga terpenuhi, dan sebaliknya (berlaku jika nilai \(p\) sama untuk semua kategori). Sifat unik lain dari sebaran Weibull adalah log-log plot, yaitu plot antara \(log(t)\) dan \(-log(-log(S(t)))\), linier.

Ilustrasi

Pada ilustrasi ini akan kita gunakan data kelangsungan hidup pada pasien dengan kanker paru stadium lanjut dari North Central Cancer Treatment Group.

#install.packages('SurvRegCensCov')
library(SurvRegCensCov)
library(survival)
library(dplyr)
library(bshazard)
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 ...

Pertama akan diperiksa terlebih dahulu log-log plot.

plot(survfit(Surv(time,status)~as.factor(sex), data=cancer), 
     col=c("black", "red"), fun="cloglog",
     xlab='log(t)', ylab='-log(-log(S(t)))')

Berdasarkan plot log-log tampak pola plot yang paralel dan linier, sehingg asumsi PH terpenuhi.

Selanjutnya dilakukan pemodelan regresi Weibull.

cancer=cancer%>%
  mutate(sex=as.factor(sex))%>%
  mutate(status=recode(status,'1'=0, '2'=1))
fit=survreg(Surv(time,status)~sex, data=cancer, dist='weibull')
summary(fit)
## 
## Call:
## survreg(formula = Surv(time, status) ~ sex, data = cancer, dist = "weibull")
##               Value Std. Error     z       p
## (Intercept)  5.8842     0.0720 81.76 < 2e-16
## sex2         0.3956     0.1276  3.10  0.0019
## Log(scale)  -0.2809     0.0619 -4.54 5.7e-06
## 
## Scale= 0.755 
## 
## Weibull distribution
## Loglik(model)= -1148.7   Loglik(intercept only)= -1153.9
##  Chisq= 10.4 on 1 degrees of freedom, p= 0.0013 
## Number of Newton-Raphson Iterations: 5 
## n= 228

Ingat kembali bahwa hasil diatas merupakan hasil model AFT. Berdasarkan hasil di atas diperoleh nilai faktor akselerasi adalah 1.485 dan hazard ratio 0.592, sementara shape parameter yaitu 1.324. Nilai shape parameter lebih dari 1 yang dapat diartikan bahwa fungsi hazard meningkat seiring dengan meningkatnya waktu. Hal ini bisa dikonfirmasi dari plot fungsi hazard.

hazard=bshazard(Surv(time,status)~sex, data=cancer, verbose = F)
## NOTE: entry.status has been set to 0 for all.
plot(hazard,
     conf.int = F, overall = F,
     lty = 1:2, ylim = c(0,0.008))

Untuk mempermudah melihat nilai parameter asli dapat digunakan fungsi ConvertWeibull

ConvertWeibull(fit, conf.level = 0.95)
## $vars
##             Estimate           SE
## lambda  0.0004127439 0.0002079792
## gamma   1.3243489259 0.0820015808
## sex2   -0.5238833093 0.1667827835
## 
## $HR
##             HR        LB        UB
## sex2 0.5922163 0.4270853 0.8211947
## 
## $ETR
##           ETR       LB       UB
## sex2 1.485242 1.156491 1.907447

Berdasarkan hasil di atas nilai gamma adalah nilai shape parameter. HR merupakan nilai Hazard Ratio sementara ETR merupakan nilai perubahan fungsi daya tahan (faktor akselerasi), dengan LB dan UB masing-masing merupakan batas bawah dan batas atas. Nilai hazard ratio<1 berarti bahwa laju bahaya pasien perempuan (sex2) lebih kecil dibanding pasien laki-laki (0.592 kali ), atau laju bahaya pasien laki-laki (1.689 kali) lebih besar dibanding pasien perempuan. Sementara nilai faktor akselerasi>1 berarti bahwa daya tahan pasien perempuan lebih besar (1.485 kali) daripada pasien laki-laki.

Selanjutnya dilakukan pemodelan dengan beberapa kovariat. Untuk menentukan kovariat mana saja yang signifikan dilakukan uji Wald.

#install.packages('rms')
library(rms)
cancer=cancer%>%
  mutate(ph.ecog=as.factor(ph.ecog))
all=psm(Surv(time,status)~sex+age+ph.ecog+pat.karno+ph.karno+meal.cal+wt.loss, 
        data=cancer, dist='weibull')
anova(all)
##                 Wald Statistics          Response: Surv(time, status) 
## 
##  Factor     Chi-Square d.f. P     
##  sex         7.38      1    0.0066
##  age         0.51      1    0.4764
##  ph.ecog    12.96      3    0.0047
##  pat.karno   1.87      1    0.1712
##  ph.karno    4.48      1    0.0344
##  meal.cal    0.01      1    0.9373
##  wt.loss     3.18      1    0.0748
##  TOTAL      30.23      9    0.0004

Berdasarkan Uji Wald di atas kovariat yang signifikan adalah `sex`,`ph.ecog`, dan `ph.karno`. Ketiga kovariat tersebut selanjutnya digunakan untuk pemodelan regresi Weibull.


```r
fit2=survreg(Surv(time,status)~sex+ph.ecog+ph.karno, data=cancer, dist='weibull')
ConvertWeibull(fit2, conf.level = 0.95)
## $vars
##               Estimate           SE
## lambda    5.736305e-05 6.328531e-05
## gamma     1.385144e+00 8.612418e-02
## sex2     -5.609872e-01 1.689725e-01
## ph.ecog1  5.770130e-01 2.346729e-01
## ph.ecog2  1.299742e+00 3.551540e-01
## ph.ecog3  2.398262e+00 1.076346e+00
## ph.karno  1.266276e-02 9.515096e-03
## 
## $HR
##                  HR        LB         UB
## sex2      0.5706454 0.4097667  0.7946868
## ph.ecog1  1.7807115 1.1241897  2.8206391
## ph.ecog2  3.6683504 1.8287872  7.3583163
## ph.ecog3 11.0040327 1.3346422 90.7274864
## ph.karno  1.0127433 0.9940314  1.0318074
## 
## $ETR
##                ETR         LB        UB
## sex2     1.4993069 1.17768905 1.9087562
## ph.ecog1 0.6593026 0.47404783 0.9169537
## ph.ecog2 0.3912750 0.23752408 0.6445499
## ph.ecog3 0.1770333 0.03921629 0.7991776
## ph.karno 0.9908998 0.97771748 1.0042599

Berdasarkan tabel HR, Nilai Hazard Ratio antara pasien laki-laki dan perempuan adalah 0.571, sementara nilai 1.781, 3.668, dan 11.004 berturut-turut dapat diartikan sebagai nilai hazard ratio antara pasien yang memiliki skor ECOG 1 (symptomatic but completely ambulatory), pasien yang memiliki skor ECOG 2 (in bed <50% of the day), pasien yang memiliki skor ECOG 3 (in bed > 50% of the day but not bedbound) dengan pasien yang memiliki skor ECOG 0 (symptomatic), dan nilai 1.0127 merupakan perubahan laju hazard setiap kenaikan satu satuan ph.karno.