library(survival)
library(survminer)
leuk <- read.csv('dados.csv', sep = ";")
leuk$group <- factor(leuk$group, levels = c(1, 0))
head(leuk)
##   id time status sex logWBC group
## 1  1   35      0   1   1.45     0
## 2  2   34      0   1   1.47     0
## 3  3   32      0   1   2.20     0
## 4  4   32      0   1   2.53     0
## 5  5   25      0   1   1.78     0
## 6  6   23      1   1   2.57     0

Model Eksponensial

Default software R adalah model AFT, artinya output di dibawah merupakan koefisien regresi dari model AFT.

summary(model_exp)

model_exp <- survreg(Surv(time, status) ~ group, dist = "exponential", data = leuk)
summary(model_exp)
## 
## Call:
## survreg(formula = Surv(time, status) ~ group, data = leuk, dist = "exponential")
##             Value Std. Error    z       p
## (Intercept) 2.159      0.218 9.90 < 2e-16
## group0      1.527      0.398 3.83 0.00013
## 
## Scale fixed at 1 
## 
## Exponential distribution
## Loglik(model)= -108.5   Loglik(intercept only)= -116.8
##  Chisq= 16.49 on 1 degrees of freedom, p= 4.9e-05 
## Number of Newton-Raphson Iterations: 4 
## n= 42

scale fixed at 1 adalah nilai p untuk weibull untuk membentuk distribusi eksponensial

Karena model diatas merupakan AFT maka koefisiennya merupakan \(\hat{\alpha}_0\) dan \(\hat{\alpha}_1\).

Dari model diatas kita juga dapat menghitung estimasi faktor acceleratednya dengan formula berikut :

\[ \hat{\gamma} = e^{\hat{\alpha}_1} \]

Accelerated Factor

alpha_0 <- unname(model_exp$coefficients[1])
alpha_1 <- unname(model_exp$coefficients[2])

gamma_1 <- exp(alpha_1)
gamma_1
## [1] 4.602564
lower_gamma <- exp(alpha_1 - 1.96 * 0.398)
upper_gamma <- exp(alpha_1 + 1.96 * 0.398)

c(lower_gamma, upper_gamma)
## [1]  2.109674 10.041169
df <- data.frame(q = c(0.25, 0.5, 0.75))
df$t_group0 <- -log(df$q) * exp(alpha_0)
df$t_group1 <- -log(df$q) * exp(alpha_0 + alpha_1)
df$gamma <- df$t_group1 / df$t_group0

df
##      q  t_group0 t_group1    gamma
## 1 0.25 12.014551 55.29774 4.602564
## 2 0.50  6.007276 27.64887 4.602564
## 3 0.75  2.493245 11.47532 4.602564

Nilai q berapapun, menghasilkan nilai faktor accelerated yang sama yaitu 4.602. Artinya variabel TRT memenuhi asumsi AFT

Treatment efektif menunda waktu kekambuhan leukimia dengan membentangkan waktu survive menjadi 4.62 kali lebih panjang dari mereka yang tidak mendapatkan treatment.

Hazard Ratio

(hr1 <- exp(-alpha_1)) #hr1 sama dengan hr
## [1] 0.2172702
hr <- 1 / gamma_1
hr
## [1] 0.2172702
lower_hr <- exp(-alpha_1 - 1.96 * 0.398)
upper_hr <- exp(-alpha_1 + 1.96 * 0.398)

c(lower_hr, upper_hr)
## [1] 0.0995900 0.4740068

interval HR yang signifikan di bawah 1, dapat disimpulkan bahwa hazard pada placebo group (0) lebih rendah dibanding hazard pada treatment group (1), artinya treatment yang dilakukan signifikan menurunkan resiko kambuh pasien kanker leukimia.

Resiko seseorang yang group 0 (treatment) untuk mengalami kambuh 0.22 kali lebih rendah dibandingkan seseorang yang tidak diberi treatment (placebo) (referensi)

HR = Hazard Placebo / Hazard Treatment (referensi)

Model Cox PH

model_cox <- coxph(Surv(time, status) ~ group, data = leuk)
summary(model_cox)
## Call:
## coxph(formula = Surv(time, status) ~ group, data = leuk)
## 
##   n= 42, number of events= 30 
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)    
## group0 -1.5721    0.2076   0.4124 -3.812 0.000138 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## group0    0.2076      4.817   0.09251    0.4659
## 
## Concordance= 0.69  (se = 0.041 )
## Likelihood ratio test= 16.35  on 1 df,   p=5e-05
## Wald test            = 14.53  on 1 df,   p=1e-04
## Score (logrank) test = 17.25  on 1 df,   p=3e-05

Model Weibull

ketika model weibull memiliki p = 1 maka akan menjadi model eksponensial

model_wei <- survreg(Surv(time, status) ~ group, dist = "weibull", data = leuk)
summary(model_wei)
## 
## Call:
## survreg(formula = Surv(time, status) ~ group, data = leuk, dist = "weibull")
##              Value Std. Error     z       p
## (Intercept)  2.248      0.166 13.55 < 2e-16
## group0       1.267      0.311  4.08 4.5e-05
## Log(scale)  -0.312      0.147 -2.12   0.034
## 
## Scale= 0.732 
## 
## Weibull distribution
## Loglik(model)= -106.6   Loglik(intercept only)= -116.4
##  Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06 
## Number of Newton-Raphson Iterations: 5 
## n= 42

Accelerated Factor

(alpha_01 <- unname(model_wei$coefficients[1]))
## [1] 2.248352
(alpha_11 <- unname(model_wei$coefficients[2]))
## [1] 1.267335
(gamma_11 <- exp(alpha_11))
## [1] 3.551374
lower_gamma1 <- exp(alpha_11 - 1.96 * 0.311)
upper_gamma1 <- exp(alpha_11 + 1.96 * 0.311)

c(lower_gamma1, upper_gamma1)
## [1] 1.930491 6.533185

Hazard Ratio

p1 <- model_wei$scale
hr1 <- exp(-alpha_11 * p1)
hr1
## [1] 0.3953692
lower_hr1 <- exp(-alpha_11*p1 - 1.96 * 0.311)
upper_hr1 <- exp(-alpha_11*p1 + 1.96 * 0.311)

c(lower_hr1, upper_hr1)
## [1] 0.2149187 0.7273298
ggsurvplot(
  survfit(Surv(time, status) ~ group, data = leuk),
  data = leuk,
  fun = 'cloglog',
  linetype = 1,
) %++%
  scale_x_continuous(
    'log(t)',
    trans = scales::log_trans(),
    labels = function(x) log(x)
  ) %++%
  geom_point(aes(col = group)) 
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

Likelihood Ratio Test

H0 : Model exponensial lebih baik

H1 : Model weibull lebih baik

lmtest::lrtest(model_exp, model_wei)
## Likelihood ratio test
## 
## Model 1: Surv(time, status) ~ group
## Model 2: Surv(time, status) ~ group
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1   2 -108.52                       
## 2   3 -106.58  1 3.8891     0.0486 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keputusan Tolak Ho yang berarti bahwa model weibull lebih baik dibandingkan model exponensial untuk data tersebut.

AIC(model_wei)
## [1] 219.159
AIC(model_exp)
## [1] 221.0481