Overplotting Non-parametric vs Parametric survival models

References

• How to plot the survival curve generated by survreg (package survival of R)? on StackOverflow
• http://stackoverflow.com/questions/9151591/how-to-plot-the-survival-curve-generated-by-survreg-package-survival-of-r

Prepare data

## Load necessary packages
library(survival)
library(rms)

## Load data
data(lung, package = "survival")

## Create Surv object
lung$SurvObj <- with(lung, Surv(time, status))

## Check first 6 rows
head(lung)

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss SurvObj
1    3  306      2  74   1       1       90       100     1175      NA    306 
2    3  455      2  68   1       0       90        90     1225      15    455 
3    3 1010      1  56   1       0       90        90       NA      15   1010+
4    5  210      2  57   1       1       90        60     1150      11    210 
5    1  883      2  60   1       0      100        90       NA       0    883 
6   12 1022      1  74   1       1       50        80      513       0   1022+

Survival curves without covariates

## Kaplan-Meier estimator without grouping
km.null <- survfit(data = lung, SurvObj ~ 1)
survplot(km.null, conf = "none")

## Overplot estimation from Cox regression by Efron method
cox.null <- coxph(data = lung, SurvObj ~ 1)
lines(survfit(cox.null), col = "green", mark.time = FALSE)

## Parametric estimation with Weibull distribution
weibull.null <- survreg(data = lung, SurvObj ~ 1, dist = "weibull")
lines(x = predict(weibull.null, type = "quantile", p = seq(0.01, 0.99, by=.01))[1,],
      y = rev(seq(0.01, 0.99, by = 0.01)),
      col = "red")

## Parametric estimation with log-logistic distribution
loglogistic.null <- survreg(data = lung, SurvObj ~ 1, dist = "loglogistic")
lines(x = predict(loglogistic.null, type = "quantile", p = seq(0.01, 0.99, by=.01))[1,],
      y = rev(seq(0.01, 0.99, by = 0.01)),
      col = "blue")

## Add legends
legend(x = "topright",
       legend = c("Kaplan-Meier", "Cox (Efron)", "Weibull", "Log-logistic"),
       lwd = 2, bty = "n",
       col = c("black", "green", "red", "blue"))

Survival curves with covariates

Sex is added as a predictor.

## Load plyr for List to (none) PLY (l_ply) for looping
library(plyr)

## Kaplan-Meier estimator
km.sex <- survfit(data = lung, SurvObj ~ sex)

survplot(km.sex, conf = "none")


## Cox regression by Efron method
cox.sex <- coxph(data = lung, SurvObj ~ sex)

l_ply(1:2, function(X) {
    lines(survfit(cox.sex, newdata = data.frame(sex = X)), col = "green", lty = X, mark.time = FALSE)
})


## Parametric with Weibull distribution
weibull.sex <- survreg(data = lung, SurvObj ~ sex, dist = "weibull")

l_ply(1:2, function(X) {
    lines(x = predict(weibull.sex, newdata = data.frame(sex = X), type = "quantile", p = seq(0.01, 0.99, by=.01)),
          y = rev(seq(0.01, 0.99, by = 0.01)),
          col = "red", lty = X)
})


## Parametric with log-logistic distribution
loglogistic.sex <- survreg(data = lung, SurvObj ~ sex, dist = "loglogistic")

l_ply(1:2, function(X) {
    lines(x = predict(loglogistic.sex, newdata = data.frame(sex = X), type = "quantile", p = seq(0.01, 0.99, by=.01)),
          y = rev(seq(0.01, 0.99, by = 0.01)),
          col = "blue", lty = X)
})

## Add legends
legend(x = "topright",
       legend = c("Kaplan-Meier", "Cox (Efron)", "Weibull", "Log-logistic"),
       lwd = 2, bty = "n",
       col = c("black", "green", "red", "blue"))