Survival analysis with R

The R framework has (strictly speaking) no abilities to perform survival analysis. The survival package fill in this gap. You can also perform “modern” graphics using the survminer package. To be able to use them, you must load it each time your launch a new session by invoking

install.packages(c("survival", "survminer"))##need to be done once

library(survival)
library(survminer)

Loading required package: ggplot2
Loading required package: ggpubr

Attaching package: ‘survminer’

The following object is masked from ‘package:survival’:

    myeloma

One main characteristics of survival data is censoring. How to deal with such data using the survival package?

The main function to handle such data is the Surv function which apparently is doing nothing except that now data belongs to the class Surv. It is thus mandatory to use it extensively in your analysis. A typical use is

duration <- 1:5##this are your observed duration eventually right censored
delta <- c(0, 1, 1, 1, 0)##this the binary vector specifying if an observation is right censored or not (0: censored, 1: event)
Surv(duration, delta)

[1] 1+ 2  3  4  5+

1. Generator fans

We are going to work with the genfan data included in the survival package.

• Load the data genfan and have a look at its documentation.
• Read the documentation of the survfit function to be sure how to use it properly.
• Fit a survival curve using the Kaplan–Meier estimator.
• What are the outputs of the summary and print functions when applied to your just fitted Kaplan–Meier estimates?
• What about the ggsurvplot function?

2. Motor insulation

We are now going to work with the imotor data included in the survival package.

• Load the data imotor and have a look at its documentation.
• For each temperature, fit a survival function using the Kaplan–Meier estimator.
• If, in the previous question, you make several call to the survfit function, try to make use of an R formula to do it in a single call.
• Analyze the outputs of the summary and print functions when applied to your just fitted Kaplan–Meier estimates?
• Plot the estimated survival curves (you may want to pass conf.int = 0.95 to get confidence intervals)

3. Confidence intervals

Recall that a widely use confidence interval for some parameter \(\theta\) is a symmetric one which is typically of the form \[ \left[\hat{\theta} - z_{1 - \alpha / 2} \sqrt{\mbox{Var}(\hat{\theta})}, \hat{\theta} - z_{1 - \alpha / 2} \sqrt{\mbox{Var}(\hat{\theta})} \right], \]
where \(\hat{\theta}\) is an estimator of \(\theta\) which is supposed here to be at least asymptotically normal and whose standard error is \(\sqrt{\mbox{Var}(\hat{\theta})}\).

• Have a look at the lecture note about confidence intervals for survival curves and make sure you understand why we introduce the *log-log survivor ship.
• Have a look at the argument of the survfit function and learn how to use it.
• On the imotor data set, compare the different type of confidence intervals.

5. How dumb see survivals? (Optional)

In this exercise, we are going to see what happens if we were very dumb and completely ignore censoring. To this aim we will work on simulated data to “know the truth” (and thus see how dumb we are!). Recall that right censoring corresponds to the situation where we observe realizations from the random variable \[T = \min(C, T_*),\] where \(C\) is a random variable related to censoring and \(T_*\) is the time to failure (and only the latter is of interest). Throughout this exercise, we will assume that both \(T_*\) and \(C\) have an exponential distribution with parameter \(\lambda_*\) and \(\lambda_c\) respectively and \(C\) and \(T_*\) are independent.

Without loss of generality, in our simulation study, we will assume that \(\lambda_c = 1\) and will only vary \(\lambda_*\).

• Generate 100 independent realizations of \(T\) with a \(\lambda_*\) value of your choice and make sure to keep track of the indicator variable indicating censoring.
• Estimate the survival curve on this data set and compare it to the true one.
• Estimate the survival curve ignoring censoring, i.e., we have no censored observation, and compare it to the theoretical surrvival curve.
• Play a bit with the value of \(\lambda_*\) to see how it impacts results and number of censored observations. Try to make sense of it.
• We now turn to a theoretical derivation of what we experienced.
  • What is the distribution of \(T\)? *Hint: You may want to compute its survivor function
  • What is the probability of being censored?
  • Are those results in agreement with your conclusions?
• We now switch back to simulation and focus on median survival time.
  • Learn how to retrieve the median survival time estimate using the quantile function. *{Hint: you may want to have a look at the documentation of quantile.survfit
  • Generate a bunch of samples of size 100, say \(N\), and estimate the mean square error for the median survival time, i.e., \[\widehat{MSE}(\lambda_*) = \frac{1}{N} \sum_{i=1}^N \left(\hat t_{0.5,i} - \frac{\log 2}{\lambda_*} \right)^2, \] where \(\hat t_{0.5,i}\) is the estimated median survival time for the \(i\)–th sample and \(\log 2 / \lambda_*\) is the theoretical one (check it!).
  • Plot the evolution of these mean square error as \(\lambda_*\) changes for both the smart way, i.e., taking into account censoring, and the dumb way, i.e., ignoring censoring.

6. Tongue cancer

A study was conducted on the effects of ploidy on the prognosis of patients with cancers of the mouth. Patients were selected who had a paraffin-embedded sample of the cancerous tissue taken at the time of surgery. Follow-up survival data was obtained on each patient. The tissue samples were examined using a flow cytometer to determine if the tumor had an aneuploid (abnormal) or diploid (normal) DNA profile using a technique discussed in Sickle–Santanello et al. (1988). Times are in weeks.

• The dataset is available from the KMsurv package under the name *tongue. (You may want to install it and load the library)
• Is there an impact of the tumour DNA profile in survival?

7. Cox’s regression: the basics

In this exercise we will learn how to make use of the *Cox’s proportional hazards model.

• Have a look at the documentation of the coxph function and make sure you understand its basic use.
• Does it make sense to use Cox’s model on the genfan data set?
• Run the following piece of code and try to understand what’s going on:

fit <- coxph(Surv(hours, status) ~ 1, data = genfan)
fit

Call:  coxph(formula = Surv(hours, status) ~ 1, data = genfan)

Null model
  log likelihood= -45.45218 
  n= 70 

plot(survfit(fit))

• We now focus on the imotor dataset. Using Cox’s model, how would you check if the temperature has an impact on survival? Is this in agreement with the log-rank test?
• Read the documentation of the residuals.coxph function and make sure to make the connection with the lecture notes.
• Perform a residual analysis to check if everything is OK with your fitted model. Hint: the function lowess might be useful to plot a local polynomial curve to your scatterplots.
• Read the documentation of the function cox.zph and use it to assess whether the proportional hazards assumption is sensible on the imotor dataset.
• Predict the survival for temperature of \(250, 280\) and \(310\). *Hint: we already saw how to predict survival curve from a fitted Cox’s model in this exercise

8. Start working on your project (if grading is about project not exam)

This exercise is a actually not an exercise but explain the expectation for grading this lecture. You will have to conduct a whole survival analysis on a dataset chosen from the ones listed below. You should write a technical report (and include your R script as a separate file) and upload it using the dedicated place on Hippocampus.

• melanoma, available from the boot package, consists of measurements made on patients with malignant melanoma;
• bfeed, available from the KMsurv package, consists of the information from 927 first–born children to mothers who chose to breast feed their children;
• burn, available from the KMsurv package, consists of measurements made on burn patients and the time until staphylococcus infection was recorded (in addition to other features);
• a dataset of your choice which is ``more related to your future job’’ ;-)