########### Accelerated Failure Time Models ##########

library(survival)
library(MASS)

data(VA)
attach(VA)
# VA data 
VA <- data.frame(stime,status,treat=factor(treat),age,Karn,diag.time,
                 cell=factor(cell), prior=factor(prior))
# Let's fit the VA data using different models and then compare results 
# Fit a Cox proportional hazards model
va.cox <- coxph(Surv(stime, status) ~ Karn + cell, data=VA)
# Fit a Weibull model, which is an accelerated failure time model 
va.wei <- survreg(Surv(stime, status) ~ Karn + cell, data=VA, dist='weibull')
# Fit a lognormal model, which is also an accelerated failure time model 
va.lnorm <- survreg(Surv(stime, status) ~ Karn + cell, data=VA, dist='lognormal') 
# Now, let's compare the results 
> summary(va.cox)
Call:
coxph(formula = Surv(stime, status) ~ Karn + cell, data = VA)

  n= 137 
         coef exp(coef) se(coef)     z       p
Karn  -0.0311      0.97  0.00518 -6.00 2.0e-09
cell2  0.7153      2.04  0.25269  2.83 4.6e-03
cell3  1.1577      3.18  0.29294  3.95 7.7e-05
cell4  0.3256      1.38  0.27668  1.18 2.4e-01

      exp(coef) exp(-coef) lower .95 upper .95
Karn       0.97      1.032     0.960      0.98
cell2      2.04      0.489     1.246      3.36
cell3      3.18      0.314     1.792      5.65
cell4      1.38      0.722     0.805      2.38

Rsquare= 0.352   (max possible= 0.999 )
Likelihood ratio test= 59.4  on 4 df,   p=3.93e-12
Wald test            = 61.3  on 4 df,   p=1.58e-12
Score (logrank) test = 64  on 4 df,   p=4.31e-13

> summary(va.wei)

Call:
survreg(formula = Surv(stime, status) ~ Karn + cell, data = VA, 
    dist = "weibull")
              Value Std. Error      z        p
(Intercept)  3.4806    0.34048 10.223 1.57e-24
Karn         0.0292    0.00462  6.310 2.78e-10
cell2       -0.7082    0.22609 -3.132 1.74e-03
cell3       -1.1085    0.25269 -4.387 1.15e-05
cell4       -0.3220    0.25002 -1.288 1.98e-01
Log(scale)  -0.0642    0.06596 -0.973 3.30e-01

Scale= 0.938 

Weibull distribution
Loglik(model)= -716.5   Loglik(intercept only)= -748.1
        Chisq= 63.15 on 4 degrees of freedom, p= 6.3e-13 
Number of Newton-Raphson Iterations: 5 
n= 137 
# Notice that Cox model does not have an intercept, while 
# the AFT models have an intercept term. The estimates from 
# the Cox model are different from that from the AFT model 
# because their interpretations are different. In the Weibull 
# model, the estimated scale is near 1, so an exponential 
# model may fit the data.  

> summary(va.lnorm)

Call:
survreg(formula = Surv(stime, status) ~ Karn + cell, data = VA, 
    dist = "lognormal")
              Value Std. Error      z        p
(Intercept)  2.2613     0.3430  6.592 4.33e-11
Karn         0.0374     0.0048  7.801 6.16e-15
cell2       -0.5334     0.2455 -2.173 2.98e-02
cell3       -0.6543     0.2805 -2.332 1.97e-02
cell4        0.1126     0.2808  0.401 6.88e-01
Log(scale)   0.0707     0.0626  1.130 2.59e-01

Scale= 1.07 

Log Normal distribution
Loglik(model)= -716.2   Loglik(intercept only)= -749.5
        Chisq= 66.62 on 4 degrees of freedom, p= 1.2e-13 
Number of Newton-Raphson Iterations: 4 
n= 137 
# The Weibull model and the lognormal model seem to produce 
# somewhat different estimates, but the overall conclusions 
# do not differ much. Model diagnostics are needed to 
# check which model fits better. 

# If the accelerated-life model holds, T*exp(-beta*x) has the same distribution 
# for all subjects. 
ntimes <- VA$stime * exp(-va.wei$linear.predictors)
plot(survfit(Surv(ntimes,VA$status)~1), log=T)
# The figure is plausibly linear, confirming the suitability of a Weibull 
# model (exponential model) 

# Motors data 
data(motors)
attach(motors)
motor.cox <- coxph(Surv(time, cens) ~ temp, motors)
motor.wei <- survreg(Surv(time,cens) ~ temp, motors)
> summary(motor.wei)

Call:
survreg(formula = Surv(time, cens) ~ temp, data = motors)
              Value Std. Error     z         p
(Intercept) 16.3185    0.62296  26.2 3.03e-151
temp        -0.0453    0.00319 -14.2  6.74e-46
Log(scale)  -1.0956    0.21480  -5.1  3.38e-07

Scale= 0.334 

Weibull distribution
Loglik(model)= -147.4   Loglik(intercept only)= -169.5
        Chisq= 44.32 on 1 degrees of freedom, p= 2.8e-11 
Number of Newton-Raphson Iterations: 7 
n= 40 

# prediction for the centre and quantiles 
predict(motor.wei, data.frame(temp=130), type="quantile", p=c(0.5,0.1))
[1] 29913.58 15934.59
>