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}\)

Seperti dalam sebaran Weibull, dalam sebaran log-logistik \(\lambda\) diparameterisasi sebagai \(\lambda^{1/p}=exp(\alpha_0+\sum_{i=1}^p \alpha_iX_i)\) (untuk model AFT). Sementara shape parameter menentukan bentuk dari fungsi hazard, jika \(p\le 1\) maka fungsi hazard menurun seiring meningkatnya waktu, dan jika \(p>1\) maka fungsi hazard meningkat sampai suatu titik waktu tertentu kemudian menurun seiring meningkatnya waktu (kondisi ini disebut unimodal).

Sebaran log-logistik memenuhi asumsi model AFT dan model PO (Proportional Odd). Model PO daya tahan merupakan model dengan rasio odd dari fungsi daya tahan diasumsikan konstan (jika dalam model PH rasio dari fungsi hazard yang bernilai konstan). Asumsi log-logistik dapat diperiksa melalui plot antara \(log(t)\) dan \(log(\frac{\hat{S}(t)}{1-\hat{S}(t)})\). Jika plot berupa garis lurus (linier) maka asumsi log-logistik terpenuhi. Selain itu, jika plot antar kelompok paralel maka asumsi PO terpenuhi. Hubungan nilai koefisien dari parameterisasi model PO dan AFT adalah sebagai berikut: \(\beta_j=-\alpha_j p\).

Ilustrasi

Pada ilustrasi ini digunakan data uji coba pasien kanker payudara tahun 1984-1989 yang dilakukan oleh German Breast Cancer Study Group (GBSG).

library(survival)
library(dplyr)
library(bshazard)
library(SurvRegCensCov)
data("gbsg")
str(gbsg)
## 'data.frame':    686 obs. of  11 variables:
##  $ pid    : int  132 1575 1140 769 130 1642 475 973 569 1180 ...
##  $ age    : int  49 55 56 45 65 48 48 37 67 45 ...
##  $ meno   : int  0 1 1 0 1 0 0 0 1 0 ...
##  $ size   : int  18 20 40 25 30 52 21 20 20 30 ...
##  $ grade  : int  2 3 3 3 2 2 3 2 2 2 ...
##  $ nodes  : int  2 16 3 1 5 11 8 9 1 1 ...
##  $ pgr    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ er     : int  0 0 0 4 36 0 0 0 0 0 ...
##  $ hormon : int  0 0 0 0 1 0 0 1 1 0 ...
##  $ rfstime: int  1838 403 1603 177 1855 842 293 42 564 1093 ...
##  $ status : int  0 1 0 0 0 1 1 0 1 1 ...

Yang menjadi interest kita adalah membandingkan pasien yang memperoleh terapi hormonal dan tidak memperoleh terapi hormonal. Pertama kita lihat plot fungsi hazard (untuk melihat perbedaan pola fungsi hazard antara distribusi Eksponensial, Weibull dan Log-Logistik).

hazard=bshazard(Surv(rfstime,status)~as.factor(hormon), data=gbsg, verbose = F)
## NOTE: entry.status has been set to 0 for all.
plot(hazard,
     conf.int = F, overall = F,
     lty=1:2)
legend(2000,0.0002,
       c('No','Yes')[1:2],
       lty=1:2)

Plot fungsi hazard memiliki pola yang fluktiatif (tidak konstan atau monoton naik), sehingga digunakan distribusi Log-logistik untuk data tersebut.

Selanjutnya diperiksa plot untuk asumsi log-logistik dan PO.

plot(survfit(Surv(rfstime,status)~as.factor(hormon), data=gbsg), 
     col=c("blue", "red"), lty=1:2,
     fun= function(s) log(s/(1-s)),log = 'x',
     xlab='log(t)', ylab='log(S(t)/(1-S(t)))')
legend(10,2, c('No','Yes')[1:2], lty=1:2, col=c("blue", "red"))

Berdasarkan plot di atas tampak pola plot yang cenderung paralel dan linier, sehingg asumsi log-logistik dan PO terpenuhi, yang berarti bahwa asumsi AFT juga terpenuhi.

Selanjutnya pemodelan regresi loglogistik menggunakan fungsi survreg dengan argument dist='loglogistic'.

fit=survreg(Surv(rfstime,status)~as.factor(hormon), data=gbsg, dist='loglogistic')
summary(fit)
## 
## Call:
## survreg(formula = Surv(rfstime, status) ~ as.factor(hormon), 
##     data = gbsg, dist = "loglogistic")
##                      Value Std. Error      z      p
## (Intercept)         7.2842     0.0624 116.79 <2e-16
## as.factor(hormon)1  0.3291     0.1035   3.18 0.0015
## Log(scale)         -0.4400     0.0487  -9.04 <2e-16
## 
## Scale= 0.644 
## 
## Log logistic distribution
## Loglik(model)= -2622.8   Loglik(intercept only)= -2627.9
##  Chisq= 10.22 on 1 degrees of freedom, p= 0.0014 
## Number of Newton-Raphson Iterations: 4 
## n= 686

Ingat kembali bahwa hasil diatas merupakan hasil model AFT. Berdasarkan hasil di atas diperoleh nilai faktor akselerasi adalah 1.38967 dan shape parameter yaitu 1.552745. Sementara nilai koefisien dari parameterisasi model PO dapat dihitung berdasarkan hubungan koefisien parameterisasi dengan model AFT \(\beta_j=-\alpha_j p\), yaitu \(\beta_1=-0.51096\), sehingga diperoleh odd rasio kegagalan sebesar 0.59992. Nilai rasio odd kegagalan model PO dan faktor akselesasi model AFT dapat diakses menggunakan fungsi ConvertWeibull.

ConvertWeibull(fit, conf.level = 0.95)
## $vars
##                         Estimate           SE
## lambda              1.224278e-05 6.503753e-06
## gamma               1.552745e+00 7.561913e-02
## as.factor(hormon)1 -5.109608e-01 1.613529e-01
## 
## $HR
##                           HR       LB        UB
## as.factor(hormon)1 0.5999189 0.437269 0.8230693
## 
## $ETR
##                         ETR       LB       UB
## as.factor(hormon)1 1.389674 1.134532 1.702195

Berdasarkan hasil di atas nilai gamma adalah nilai shape parameter. Pada komponen HR mmerupakan nilai rasio odd kegagalan, sementara komponen ETR merupakan nilai faktor akselerasi, dengan LB dan UB masing-masing merupakan batas bawah dan batas atas. Nilai rasio odd kegagalan<1 berarti bahwa odd kegagalan pasien yang menerima terapi hormon lebih kecil dibanding pasien yang tidak menerima terapi hormon (0.59992 kali ), atau odd kegagalan pasien yang tidak menerima terapi hormon (1.66689 kali) lebih besar dibanding pasien yang menerima terapi hormon. Sementara nilai faktor akselerasi>1 berarti bahwa daya tahan pasien yang menerima terapi hormon lebih besar (1.38967 kali) daripada pasien yang tidak menerima terapi hormon.

Selanjutnya dilakukan pemodelan dengan beberapa peubah, dengan terlebih dahulu dilakukan seleksi peubah berdasarkan uji Wald.

library(rms)
gbsg2=gbsg%>%
  mutate(meno=as.factor(meno),
         hormon=as.factor(hormon))
all=psm(Surv(rfstime,status)~age+meno+size+grade+nodes+pgr+er+hormon, 
        data=gbsg2, dist='loglogistic')
anova(all)
##                 Wald Statistics          Response: Surv(rfstime, status) 
## 
##  Factor     Chi-Square d.f. P     
##  age          3.47     1    0.0624
##  meno         4.08     1    0.0435
##  size         4.33     1    0.0374
##  grade        7.89     1    0.0050
##  nodes       37.20     1    <.0001
##  pgr         16.98     1    <.0001
##  er           0.03     1    0.8616
##  hormon      10.92     1    0.0010
##  TOTAL      111.27     8    <.0001

Semua peubah signifikan kecuali age dan er.

fit2=survreg(Surv(rfstime,status)~meno+size+grade+nodes+pgr+hormon, 
             data=gbsg2, dist='loglogistic')
ConvertWeibull(fit2, conf.level = 0.95)
## $vars
##              Estimate           SE
## lambda   7.040853e-07 5.201077e-07
## gamma    1.746885e+00 8.467109e-02
## meno1    1.549095e-01 1.634834e-01
## size     1.177461e-02 5.672894e-03
## grade    4.131035e-01 1.439063e-01
## nodes    8.786546e-02 1.494529e-02
## pgr     -2.878002e-03 6.729780e-04
## hormon1 -5.754230e-01 1.705786e-01
## 
## $HR
##                HR        LB        UB
## meno1   1.1675522 0.8474596 1.6085466
## size    1.0118442 1.0006562 1.0231573
## grade   1.5115015 1.1400274 2.0040192
## nodes   1.0918412 1.0603226 1.1242967
## pgr     0.9971261 0.9958118 0.9984422
## hormon1 0.5624669 0.4026245 0.7857669
## 
## $ETR
##               ETR        LB        UB
## meno1   0.9151406 0.7617244 1.0994559
## size    0.9932823 0.9869950 0.9996097
## grade   0.7894016 0.6720606 0.9272302
## nodes   0.9509457 0.9353529 0.9667983
## pgr     1.0016489 1.0008954 1.0024029
## hormon1 1.3901331 1.1488596 1.6820768

Nilai rasio odd kegagalan antara pasien postmenopausal (meno1) dan premenopausal (meno0) sebesar 1.16755 berarti bahwa odd kegagalan pasien postmenopausal lebih besar 1.16755 kali dibanding pasien premenopausal dengan peubah lain dianggap konstan. Hal ini sejalan dengan nilai rasio daya tahan (faktor akselerasi) nya, yaitu 0.91514 yang berarti bahwa daya tahan pasien postmenopausal lebih kecil 0.91514 kali dibanding pasien premenopausal. Sementara untuk peubah size, nilai faktor akselerasi sebesar 0.99328 berarti bahwa setiap perubahan ukuran tumor 1 mm maka daya tahan pasien akan berubah 0.99328. Dan untuk peubah nodes, nilai rasio odd kegagalan 1.09184 berarti bahwa odd kegagalan pasien berubah sebesar 1.09184 untuk setiap bertambahnya jumlah positive lymph nodes satu satuan.