Software 10.5.1 (A Weibull model in both languages). survreg fits AFT models; its coefficients are on the log time scale. The Weibull fit for lung, its time ratios and scale, and the hazard ratios from Theorem 10.2.1 beside those of coxph:

> library(survival); options(show.signif.stars = FALSE)
> lg <- na.omit(lung[, c("time", "status", "age", "sex", "ph.ecog")])
> w <- survreg(Surv(time, status) ~ age + sex + ph.ecog, data = lg,
+              dist = "weibull")
> round(cbind(coef = coef(w)[-1], time_ratio = exp(coef(w)[-1])), 4)
           coef time_ratio
age     -0.0075     0.9926
sex      0.4011     1.4935
ph.ecog -0.3396     0.7120
> round(w$scale, 4)
[1] 0.7311
> hr_from_aft <- exp(-coef(w)[-1] / w$scale)     # HR = exp(-b / sigma)
> cx <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lg)
> round(rbind(weibull_HR = hr_from_aft, cox_HR = exp(coef(cx))), 4)
              age    sex ph.ecog
weibull_HR 1.0103 0.5778  1.5913
cox_HR     1.0111 0.5754  1.5900

Reading the output. Women’s times are \(1.49\) times men’s; the scale \(0.731\) is the reciprocal of the Weibull shape, \(k = 1.37\), an increasing hazard. Converted, the hazard ratios agree with Cox’s to two decimals.

Python: lifelines.WeibullAFTFitter; its rho_ parameter is the shape \(k\).

Software 10.5.2 (Comparing AFT families). The veteran model fitted with three distributions and compared by AIC; then the time ratios of the log-logistic with \(95\%\) intervals, and the ratio for ten Karnofsky points:

> vet <- veteran; vet$trt <- factor(vet$trt, labels = c("standard", "test"))
> f <- Surv(time, status) ~ trt + celltype + karno
> aic <- sapply(c("weibull", "lognormal", "loglogistic"),
+               function(d) AIC(survreg(f, data = vet, dist = d)))
> round(aic, 2)
    weibull   lognormal loglogistic
    1445.71     1445.81     1439.10
> ll <- survreg(f, data = vet, dist = "loglogistic")
> round(ll$scale, 4)
[1] 0.5811
> round(exp(cbind(TR = coef(ll), confint(ll)))[-1, ], 4)   # time ratios
                      TR  2.5 % 97.5 %
trttest           0.9479 0.6705 1.3401
celltypesmallcell 0.4966 0.3070 0.8033
celltypeadeno     0.4605 0.2737 0.7748
celltypelarge     0.9671 0.5765 1.6224
karno             1.0366 1.0276 1.0456
> round(exp(10 * coef(ll)["karno"]), 4)          # per ten Karnofsky points
 karno
1.4325

Reading the output. The log-logistic has the lowest AIC, by more than \(6\). Treatment has no clear effect (the interval for its time ratio contains \(1\)); small-cell and adenocarcinoma tumours halve survival times compared with squamous; ten Karnofsky points lengthen them by \(43\%\).

Python: LogLogisticAFTFitter and LogNormalAFTFitter in lifelines, with an AIC_ attribute.

Software 10.5.3 (Predictions and Cox–Snell residuals). Predicted medians at three Karnofsky scores; then the Cox–Snell residuals and the cumulative hazard of their Kaplan–Meier curve at \(r = 0.5\), \(1\) and \(2\):

> nd <- data.frame(trt = "standard", celltype = "squamous",
+                  karno = c(40, 60, 80))
> round(predict(ll, newdata = nd, type = "quantile", p = 0.5), 1)
    1     2     3
 52.3 107.2 220.0
> z <- (log(vet$time) - predict(ll, type = "lp")) / ll$scale
> r <- log(1 + exp(z))                           # Cox-Snell: H(t_i | x_i)
> km <- survfit(Surv(r, vet$status) ~ 1)         # should be exponential(1)
> H <- -log(summary(km, times = c(0.5, 1, 2))$surv)
> round(rbind(r = c(0.5, 1, 2), H_hat = H), 3)
       [,1]  [,2]  [,3]
r     0.500 1.000 2.000
H_hat 0.466 1.066 1.925

Reading the output. The medians are \(52\), \(107\) and \(220\) days, each \(1.43^2 = 2.05\) times the one before. The residuals’ cumulative hazard is close to \(r\) at each point, so the model is adequate.

Python: predict_median on a fitted lifelines AFT model.