Survival Analysis with R

April 25, 2017
By

(This article was first published on R Views, and kindly contributed to R-bloggers)

With roots dating back to at least 1662 when John Graunt, a London merchant, published an extensive set of inferences based on mortality records, Survival Analysis is one of the oldest subfields of Statistics [1]. Basic life-table methods, including techniques for dealing with censored data, were known before 1700 [2]. In the early eighteenth century, the old masters, de Moivre working on annuities and Daniel Bernoulli studying competing risks for his work on smallpox inoculation, developed the foundations of time-to-event modeling [2]. Today, survival analysis models are important in Engineering, Insurance, Marketing and Medicine and many more application areas. So, it is not surprising that the R Task View on Survival Analysis, a curated, organized and annotated list of relevant R packages and functions, is formidable.

Looking at the Task View on a small screen is a bit like standing too close to a brick wall – left-right, up-down, bricks all around. It is a fantastic edifice that gives some idea of the significant contributions R developers have made both to the theory and practice of Survival Analysis. As well-organized as it is, however, I imagine that even survival analysis experts need some time to find their way around this task view. (I would be remiss not to mention that we all owe a great deal of gratitude to Arthur Allignol and Aurielien Latouche, the task view maintainers.) Newcomers, people either new to R or new to survival analysis or both, must find it overwhelming. So, it is with newcomers in mind that I offer the following slim trajectory through the task view that relies on just a few packages: survival, KMsurv, Oisurv and ranger

The survival package, which began life as an S package in the late ’90s, is the cornerstone of the entire R Survival Analysis edifice. Not only is the package itself rich in features, but the object created by the Surv() function, which contains failure time and censoring information, is the basic survival analysis data structure in R.

KMsurv contains some interesting data sets from John Klein and Melvin Moeschberger’s classic text, Survival Analysis Techniques for Censored and Truncated Data.

My main reason for selecting the OIsurv package was to draw attention the very helpful guide to Survival Analysis in R, produced by the folks at OpenIntro.

ranger might be the surprise in my very short list of survival packages. The ranger() function is well-known for being a fast implementation of the Random Forests algorithm for building ensembles of classification and regression trees. But it was news to me that ranger() also builds survival models. Benchmarks indicate that ranger() is suitable for building time-to-event models with the large, high dimensional data sets important to internet marketing applications. Since ranger() uses standard Surv() survival objects, it’s an ideal tool for getting acquainted with survival analysis in the in this machine learning age.

*** Load and transform the data This first block of code loads the required packages along with the bone marrow transplant data frame from the KMsurv package. I chose this because has number of covariates and no missing values. The only data preparation is to make the appropriate variables into factors.

library(survival)
library(dplyr)
library(OIsurv) # Aumatically loads KMsurv
library(ranger)
library(ggplot2)

#------------
data(bmt)
# sapply(bmt,class)

# Prepare new data frame for modeling
bmt2 <- select(bmt,-c(t2,d2,d3))
bmt2  <- mutate(bmt2,
                group = as.factor(group),
                d1 = as.factor(d1),
                da = as.factor(da),
                dc = as.factor(dc),
                dp = as.factor(dp),
                z3 = as.factor(z3),
                z4 = as.factor(z4),
                z5 = as.factor(z5),
                z6 = as.factor(z6),
                z8 = as.factor(z8),
                z9 = as.factor(z9),
                z10 = as.factor(z10)
)
head(bmt2)
##   group   t1 d1   ta da  tc dc tp dp z1 z2 z3 z4 z5 z6   z7 z8 z9 z10
## 1     1 2081  0   67  1 121  1 13  1 26 33  1  0  1  1   98  0  1   0
## 2     1 1602  0 1602  0 139  1 18  1 21 37  1  1  0  0 1720  0  1   0
## 3     1 1496  0 1496  0 307  1 12  1 26 35  1  1  1  0  127  0  1   0
## 4     1 1462  0   70  1  95  1 13  1 17 21  0  1  0  0  168  0  1   0
## 5     1 1433  0 1433  0 236  1 12  1 32 36  1  1  1  1   93  0  1   0
## 6     1 1377  0 1377  0 123  1 12  1 22 31  1  1  1  1 2187  0  1   0

Kaplan Meier Analysis

The first thing to do is to use Surv() to build the standard survival object. The variable t1 records the time to death or the censored time; d1 indicates that the patient died (d1 = 1) or that the patient survived until the end of the study (d1 = 0). Note that a “+” after the time in the print out of y_bmt indicates censoring. The formula y_bmt ~ 1 instructs the survfit() function to fit a model with intercept only, and produces the Kaplan-Meier estimate.

The confBands() function from the OIsurv package estimates Hall-Wellner confidence bands. These are large-sample, simultaneous estimates and are plotted in red. They contrast with the point-wise confidence bands rendered as black dashed lines.

# Kaplan Meier Survival Curve
y_bmt <- Surv(bmt$t1, bmt$d1)
y_bmt
##   [1] 2081+ 1602+ 1496+ 1462+ 1433+ 1377+ 1330+  996+  226+ 1199+ 1111+
##  [12]  530+ 1182+ 1167+  418   417   276   156   781   172   487   716 
##  [23]  194   371   526   122  1279   110   243    86   466   262   162 
##  [34]  262     1   107   269   350  2569+ 2506+ 2409+ 2218+ 1857+ 1829+
##  [45] 1562+ 1470+ 1363+ 1030+  860+ 1258+ 2246+ 1870+ 1799+ 1709+ 1674+
##  [56] 1568+ 1527+ 1324+  957+  932+  847+  848+ 1850+ 1843+ 1535+ 1447+
##  [67] 1384+  414  2204  1063   481   105   641   390   288   522    79 
##  [78] 1156   583    48   431  1074   393    10    53    80    35  1499+
##  [89]  704   653   222  1356+ 2640+ 2430+ 2252+ 2140+ 2133+ 1238+ 1631+
## [100] 2024+ 1345+ 1136+  845+  491   162  1298   121     2    62   265 
## [111]  547   341   318   195   469    93   515   183   105   128   164 
## [122]  129   122    80   677    73   168    74    16   248   732   105 
## [133]  392    63    97   153   363
fit1_bmt <- survfit(y_bmt ~ 1)
summary(fit1_bmt)
## Call: survfit(formula = y_bmt ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    137       1    0.993 0.00727        0.979        1.000
##     2    136       1    0.985 0.01025        0.966        1.000
##    10    135       1    0.978 0.01250        0.954        1.000
##    16    134       1    0.971 0.01438        0.943        0.999
##    35    133       1    0.964 0.01602        0.933        0.995
##    48    132       1    0.956 0.01748        0.923        0.991
##    53    131       1    0.949 0.01881        0.913        0.987
##    62    130       1    0.942 0.02003        0.903        0.982
##    63    129       1    0.934 0.02117        0.894        0.977
##    73    128       1    0.927 0.02222        0.884        0.972
##    74    127       1    0.920 0.02322        0.875        0.966
##    79    126       1    0.912 0.02415        0.866        0.961
##    80    125       2    0.898 0.02588        0.848        0.950
##    86    123       1    0.891 0.02668        0.840        0.944
##    93    122       1    0.883 0.02744        0.831        0.939
##    97    121       1    0.876 0.02817        0.822        0.933
##   105    120       3    0.854 0.03017        0.797        0.915
##   107    117       1    0.847 0.03078        0.788        0.909
##   110    116       1    0.839 0.03137        0.780        0.903
##   121    115       1    0.832 0.03193        0.772        0.897
##   122    114       2    0.818 0.03300        0.755        0.885
##   128    112       1    0.810 0.03350        0.747        0.879
##   129    111       1    0.803 0.03399        0.739        0.872
##   153    110       1    0.796 0.03445        0.731        0.866
##   156    109       1    0.788 0.03490        0.723        0.860
##   162    108       2    0.774 0.03575        0.707        0.847
##   164    106       1    0.766 0.03615        0.699        0.841
##   168    105       1    0.759 0.03653        0.691        0.834
##   172    104       1    0.752 0.03690        0.683        0.828
##   183    103       1    0.745 0.03726        0.675        0.821
##   194    102       1    0.737 0.03760        0.667        0.815
##   195    101       1    0.730 0.03793        0.659        0.808
##   222    100       1    0.723 0.03825        0.651        0.802
##   243     98       1    0.715 0.03856        0.644        0.795
##   248     97       1    0.708 0.03886        0.636        0.788
##   262     96       2    0.693 0.03943        0.620        0.775
##   265     94       1    0.686 0.03969        0.612        0.768
##   269     93       1    0.678 0.03995        0.604        0.761
##   276     92       1    0.671 0.04019        0.597        0.755
##   288     91       1    0.664 0.04042        0.589        0.748
##   318     90       1    0.656 0.04063        0.581        0.741
##   341     89       1    0.649 0.04084        0.574        0.734
##   350     88       1    0.642 0.04104        0.566        0.727
##   363     87       1    0.634 0.04122        0.558        0.720
##   371     86       1    0.627 0.04140        0.551        0.713
##   390     85       1    0.619 0.04156        0.543        0.706
##   392     84       1    0.612 0.04172        0.535        0.699
##   393     83       1    0.605 0.04186        0.528        0.693
##   414     82       1    0.597 0.04199        0.520        0.686
##   417     81       1    0.590 0.04212        0.513        0.679
##   418     80       1    0.583 0.04223        0.505        0.671
##   431     79       1    0.575 0.04234        0.498        0.664
##   466     78       1    0.568 0.04243        0.490        0.657
##   469     77       1    0.560 0.04252        0.483        0.650
##   481     76       1    0.553 0.04259        0.476        0.643
##   487     75       1    0.546 0.04266        0.468        0.636
##   491     74       1    0.538 0.04271        0.461        0.629
##   515     73       1    0.531 0.04276        0.453        0.622
##   522     72       1    0.524 0.04280        0.446        0.615
##   526     71       1    0.516 0.04282        0.439        0.607
##   547     69       1    0.509 0.04285        0.431        0.600
##   583     68       1    0.501 0.04287        0.424        0.593
##   641     67       1    0.494 0.04288        0.416        0.585
##   653     66       1    0.486 0.04288        0.409        0.578
##   677     65       1    0.479 0.04286        0.402        0.571
##   704     64       1    0.471 0.04284        0.394        0.563
##   716     63       1    0.464 0.04281        0.387        0.556
##   732     62       1    0.456 0.04277        0.380        0.548
##   781     61       1    0.449 0.04272        0.372        0.541
##  1063     52       1    0.440 0.04276        0.364        0.533
##  1074     51       1    0.432 0.04278        0.355        0.524
##  1156     48       1    0.423 0.04282        0.346        0.515
##  1279     42       1    0.413 0.04297        0.336        0.506
##  1298     41       1    0.402 0.04308        0.326        0.496
##  2204      9       1    0.358 0.05696        0.262        0.489
cb <- confBands(y_bmt, type = "hall")
plot(fit1_bmt,
    main = 'Kaplan Meyer Plot with confidence bands')
lines(cb, col = "red",lty = 3)
legend(1000, 0.99, legend = c('K-M survival estimate',
'pointwise intervals', 'Hall-Werner conf bands'), lty = 1:3)

Cox Proportional Hazards Model

Next, I fit a Cox Proportional Hazards model, which makes use of several of the variables contained in the data set. I don’t pretend to have the best model here, or even a very good one. My intention is just to show the survival package’s coxph() function, and show how the covariates enter the formula. Note, however, that this model does achieve an R2 value of 0.8, and that three variables – ta, tc, and dc1 – show up as significant.

# Fit Cox Model
form <- formula(y_bmt ~ group + ta + tc + tp + dc + dp + 
                   z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z10)

cox_bmt <- coxph(form,data = bmt2)
summary(cox_bmt)
## Call:
## coxph(formula = form, data = bmt2)
## 
##   n= 137, number of events= 81 
## 
##              coef  exp(coef)   se(coef)      z Pr(>|z|)    
## group2  0.2540084  1.2891826  0.4342118  0.585   0.5586    
## group3  0.5280287  1.6955865  0.4197516  1.258   0.2084    
## ta     -0.0017962  0.9982054  0.0003736 -4.808 1.52e-06 ***
## tc     -0.0060005  0.9940175  0.0010512 -5.708 1.14e-08 ***
## tp     -0.0005709  0.9994293  0.0012622 -0.452   0.6510    
## dc1    -3.5832062  0.0277865  0.4938320 -7.256 3.99e-13 ***
## dp1    -1.2266468  0.2932744  0.4832251 -2.538   0.0111 *  
## z1      0.0109077  1.0109674  0.0230293  0.474   0.6358    
## z2     -0.0155304  0.9845895  0.0200202 -0.776   0.4379    
## z31     0.0162146  1.0163468  0.2480277  0.065   0.9479    
## z41     0.0334753  1.0340419  0.2722755  0.123   0.9021    
## z51     0.0792485  1.0824733  0.2856339  0.277   0.7814    
## z61    -0.0702048  0.9322029  0.2557463 -0.275   0.7837    
## z7      0.0000440  1.0000440  0.0004235  0.104   0.9173    
## z81     0.4618598  1.5870228  0.3386689  1.364   0.1726    
## z101    0.5145619  1.6729055  0.3330479  1.545   0.1223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## group2   1.28918     0.7757   0.55044   3.01937
## group3   1.69559     0.5898   0.74478   3.86023
## ta       0.99821     1.0018   0.99747   0.99894
## tc       0.99402     1.0060   0.99197   0.99607
## tp       0.99943     1.0006   0.99696   1.00190
## dc1      0.02779    35.9887   0.01056   0.07314
## dp1      0.29327     3.4098   0.11375   0.75613
## z1       1.01097     0.9892   0.96635   1.05764
## z2       0.98459     1.0157   0.94670   1.02399
## z31      1.01635     0.9839   0.62506   1.65258
## z41      1.03404     0.9671   0.60642   1.76319
## z51      1.08247     0.9238   0.61842   1.89474
## z61      0.93220     1.0727   0.56470   1.53887
## z7       1.00004     1.0000   0.99921   1.00087
## z81      1.58702     0.6301   0.81716   3.08218
## z101     1.67291     0.5978   0.87092   3.21338
## 
## Concordance= 0.926  (se = 0.034 )
## Rsquare= 0.806   (max possible= 0.995 )
## Likelihood ratio test= 224.6  on 16 df,   p=0
## Wald test            = 107.5  on 16 df,   p=1.332e-15
## Score (logrank) test = 211.7  on 16 df,   p=0
cox_fit_bmt <- survfit(cox_bmt)
# plot(cox_fit_bmt)

Random Forests Model

To give some idea of the scope of R’s capabilities to work with time to event data, I use the ranger() function to fit a Random Forests Ensemble model to the data. Note that ranger() builds a model for each observation in the data set. The next block of code builds the model using the same variables used the the Cox model above, and each of the 137 survival curves computed for the bmt data set, along with a curve of average values.

# ranger model

# ranger model
r_fit_bmt <- ranger(form,
                    data = bmt2,
                    importance = "permutation",
                    seed = 1234)

# Average the survival models
death_times <- r_fit_bmt$unique.death.times
surv_prob <- data.frame(r_fit_bmt$survival)
avg_prob <- sapply(surv_prob,mean)

# Plot the survival models for each patient
plot(r_fit_bmt$unique.death.times,r_fit_bmt$survival[1,], type = "l", 
     ylim = c(0,1),
     col = "red",
     xlab = "death times",
     ylab = "survival",
     main = "Patient Survival Curves")

for(n in c(2:137)){
  lines(r_fit_bmt$unique.death.times, r_fit_bmt$survival[n,], type = "l", col = "red")
}
lines(death_times, avg_prob, lwd = 2)
legend(100, 0.2, legend = c('Averages - black'))

Here, we show the ranking of variable importance computed by the permutation method, which is ranger()’s default for survival data. Note that ta, tc, and dc are the same top three variables flagged in the Cox model.

Also listed is a measure of prediction error calculated from Harrell’s c-index. This index is defined as “… the proportion of all usable patient pairs in which the predictions and outcomes are concordant” (cf. [8] p 370), where predictions for pairs are concordant if predicted survival times are larger for patients who lived longer. Note that Harrell’s c-index may be thought of as a generalization of finding the are under an ROC curve. (For binary outcomes Harrell’s c-index reduces to the Wilcoxon-Mann-Whitney statistic which, in turn, is equivalent to computing the area under the ROC curve.)

vi <- data.frame(sort(round(r_fit_bmt$variable.importance, 4), decreasing = TRUE))
names(vi) <- "importance"
head(vi)
##    importance
## ta     0.1259
## tc     0.0688
## dc     0.0190
## tp     0.0117
## dp     0.0092
## z2     0.0046
cat("Prediction Error = 1 - Harrell's c-index = ", r_fit_bmt$prediction.error)
## Prediction Error = 1 - Harrell's c-index =  0.09304771

Finally, we plot the survival curves computed for all three models on the same graph. Note that the “ad hoc” curve of average survival curves computed by the ranger model tracks the Kaplan-Meier curve fairly well.

# Set up for ggplot
km <- rep("KM", length(fit1_bmt$time))
km_df <- data.frame(fit1_bmt$time,fit1_bmt$surv,km)
names(km_df) <- c("Time","Surv","Model")

cox <- rep("Cox",length(cox_fit_bmt$time))
cox_df <- data.frame(cox_fit_bmt$time,cox_fit_bmt$surv,cox)
names(cox_df) <- c("Time","Surv","Model")

rf <- rep("RF",length(r_fit_bmt$unique.death.times))
rf_df <- data.frame(r_fit_bmt$unique.death.times,avg_prob,rf)
names(rf_df) <- c("Time","Surv","Model")
plot_df <- rbind(km_df,cox_df,rf_df)

p <- ggplot(plot_df, aes(x = Time, y = Surv, color = Model))
p + geom_line() + ggtitle("Comparison of Survival Curves") 

For a very nice exposition of the sort of predictive survival analysis modeling that can be done with ranger, be sure to have a look at Manuel Amunategui’s post and video.

This four-package excursion only hints at the Survival Analysis tools that are available in R, but it does illustrate some of the richness of the R platform which has been under continuous development and improvement for nearly twenty years. The use of the Surv() function shows how open source code allows generations of developers to build on the work of their predecessors. The ranger packages provides a practical example of how R can incorporate fast C++ code and adapt to the world of machine learning applications, and the incidental use of options such as Hall-Wellner Confidence bands and Harrell’s c-index gives some idea of the statistical depth that underlies almost everything R.

References

For convenience, I have collected the references used throughout the post here.

[1] Hacking, Ian. (2006) The Emergence of Probability: A Philosophical Study of Early Ideas about Probability Induction and Statistical Inference. Cambridge University Press, 2nd ed. p11
[2] Andersen, P.K., Keiding, N. (1998) Survival analysis (overview) Encyclopedia of Biostatistics 6. Wiley, p 4452-4461
[3] Kaplan, E.L. & Meier, P. (1958). Non-parametric estimation from incomplete observations, J American Stats Assn. 53, 457–481, 562–563.
[4] Cox, D.R. (1972). Regression models and life-tables (with discussion), Journal of the Royal Statistical Society (B) 34, 187–220.
[5] Diez, David. Survival Analysis in R. OpenIntro [6] Klein, John P and Moeschberger, Melvin L. (1997) Survival Analysis Techniques for Censored and Truncated Data, Springer.
[7] Wright, Marvin & Ziegler, Andreas. (2017) * [ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R] (https://www.jstatsoft.org/article/view/v077i01)*, JSS Vol 77, Issue 1.
[8] Harrell, Frank, Lee, Kerry & Mark, Daniel. (1996) Multivariable Prognostic Models: Issues in Developing Models, Evaluating Assumptions and Adequacy, and Measuring and Reducing Errors. Statistics in Medicine, Vol 15, 361-387
[9] Amunategui, Manuel. Survival Ensembles: Survival Plus Classification for Improved Time-Based Predictions in R

To leave a comment for the author, please follow the link and comment on their blog: R Views.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)