Time-to-event data that consists of a distinct start time and end time. These data are common in the medical field.

Examples from cancer:

• Time from surgery to death
• Time from start of treatment to progression
• Time from response to recurrence

Time-to-event data is common in many fields including, but not limited to:

• Time from HIV infection to development of AIDS
• Time to heart attack
• Time to onset of substance abuse
• Time to initiation of sexual activity
• Time to machine malfunction
• Time to employee resignation/termination

Because survival analysis is common in many other fields, it also goes by other names:

• Reliability analysis
• Duration analysis
• Event history analysis
• Time-to-event analysis

The two most common questions I encounter related to survival analysis are:

1. What is the probability of survival to a certain point in time?
2. What is the average survival time?

In survival analysis it is common for the exact event time to be unknown, or unobserved, which is called censoring. A subject may be censored due to:

• Loss to follow-up
• Withdrawal from study
• No event by end of fixed study period

Specifically these are examples of right censoring. Other common types of censoring include:

• Left
• Interval

When the exact event time is unknown then some patients are censored, and survival analysis methods are required.

In this example, how would we compute the proportion who are event-free at 10 years?

• Subjects 2, 3, 5, 6, 8, 9, and 10 were all event-free at 10 years.
• Subjects 4 and 7 had the event before 10 years.
• Subject 1 was censored before 10 years, so we don't know whether they had the event or not by 10 years. How do we incorporate this subject into our estimate?

Toy example of a Kaplan-Meier curve for this simple data (details to follow):

• Horizontal lines represent survival duration for the interval
• An interval is terminated by an event
• The height of vertical lines show the change in cumulative probability
• Censored observations, indicated by tick marks, reduces the cumulative survival between intervals

For subject \(i\):

1. Event time \(T_i\)
2. Censoring time \(C_i\)

3. Event indicator \(\delta_i\):

   • 1 if event observed (i.e. \(T_i \leq C_i\))
   • 0 if censored (i.e. \(T_i > C_i\))

4. Observed time \(Y_i = \min(T_i, C_i)\)

Organ Transplantation is limited by the supply of cadaveric and living donors.

Unfortunately, many people die while waiting on the list.

What types of characteristics lead to an increase or decrease likelihood of dying on the waiting list?

data("transplant")

• 815 of people on the waiting list
• Outcome: Death
• Predictors: age, sex, Blood type, year added to list

transplant %>% 
    head

##   age sex abo year futime    event
## 1  47   m   B 1994   1197    death
## 2  55   m   A 1991     28      ltx
## 3  52   m   B 1996     85      ltx
## 4  40   f   O 1995    231      ltx
## 5  70   m   O 1996   1271 censored
## 6  66   f   O 1996     58      ltx

• "futime": the observed time \(Y_i = min(T_i, C_i)\)
• "event": the event factor variable can take on censored, death, liver transpant, or withdraw
• "abo": blood type A, B, AB, or O
• sex: Male or Female
• age: age at addition to the waiting list

Most functions used in survival analysis will also require a binary indicator of event that is:

• 0 for no event
• 1 for event

Currently our data example contains a factor variable indicating whehter the patient has died or if they have received a liver transpant, been withdrawn, or are censored.

transplant$delta<-0
transplant$delta[transplant$event=="death"]<-1

Recall the questions of interest:

1. What is the probability of surviving to a certain point in time?
2. What is the average survival time?

The Kaplan-Meier method is the most common way to estimate survival times and probabilities. It is a non-parametric approach that results in a step function, where there is a step down each time an event occurs.

• The Surv function from the survival package creates a survival object for use as the response in a model formula. There will be one entry for each subject that is the survival time, which is followed by a + if the subject was censored. Let's look at the first 10 observations:

survival::Surv(transplant$futime, transplant$delta)[1:10]

##  [1] 1197    28+   85+  231+ 1271+   58+  392+   30+   12+  139+

• The survival::survfit function creates survival curves based on a formula. Let's generate the overall survival curve for the entire cohort, assign it to object f1, and look at the names of that object:

f1 <- survfit(Surv(transplant$futime, transplant$delta) ~ 1, data = transplant)
names(f1)

##  [1] "n"         "time"      "n.risk"    "n.event"   "n.censor" 
##  [6] "surv"      "type"      "std.err"   "lower"     "upper"    
## [11] "conf.type" "conf.int"  "call"

Some key components of this survfit object that will be used to create survival curves include:

• time, which contains the start and endpoints of each time interval
• surv, which contains the survival probability corresponding to each time

Now we plot the survfit object in base R to get the Kaplan-Meier plot:

plot(f1, 
     xlab = "Days", 
     ylab = "Overall survival probability")

• The default plot in base R shows the step function (solid line) with associated confidence intervals (dotted lines). Note that the tick marks for censored patients are not shown by default, but could be added using mark.time = TRUE

One quantity often of interest in a survival analysis is the probability of surviving a certain number (\(x\)) of years.

For example, to estimate the probability of survivng to 5 years, use summary with the times argument:

summary(f1, times = 5*365)

## Call: survfit(formula = Surv(transplant$futime, transplant$delta) ~ 
##     1, data = transplant)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##  1825      2      66     0.49   0.206        0.215            1

We find that the 5-year probability of survival in this study is 49%. The associated lower and upper bounds of the 95% confidence interval are also displayed.

What happens if you use a "naive" estimate?

66 of the 815 patients died by 5 years so:

\[\Big(1 - \frac{66}{815}\Big) \times 100 = 92\%\]

• You get an incorrect estimate of the 5-year probability of survival when you ignore the fact that 747 patients were censored before 5 years.

• Recall the correct estimate of the 5-year probability of survival was 49%.

Another quantity often of interest in a survival analysis is the average survival time, which we quantify using the median (survival times are not expected to be normally distributed so the mean is not an appropriate summary).

We can obtain this directly from our survfit object:

survival::survfit(survival::Surv(futime, delta) ~ 1, data = transplant)

## Call: survfit(formula = survival::Surv(futime, delta) ~ 1, data = transplant)
## 
##       n  events  median 0.95LCL 0.95UCL 
##     815      66    1422    1422      NA

We see the median survival time is 3.9 years. The lower and upper bounds of the 95% confidence interval are also displayed.

What happens if you use a "naive" estimate?

Summarize the median survival time among the 66 patients who died:

transplant$futime[transplant$delta == 1] %>% 
    median

## [1] 65.5

• You get an incorrect estimate of median survival time of 0.2 years when you ignore the fact that censored patients also contribute follow-up time.

• Recall the correct estimate of median survival time is 3.9 years.

Is there a difference in survival probability between groups?

From our example: does the probability of survival differ according to gender among liver transplant patients?

We can add a covariate to the right-hand side of the survival::survfit object to obtain a stratified Kaplan-Meier plot.

Let's also look at some other customization we can do with survminer::ggsurvplot.

survminer::ggsurvplot(
    fit = survival::survfit(survival::Surv(futime, delta) ~ abo, data = transplant), 
    xlab = "Days",
    ylab = "Overall survival probability",
    legend.title = "Blood Type",
    legend.labs = c("A", "B", "AB","O"),
    break.x.by = 100, 
    censor = FALSE,
    risk.table = TRUE,
    risk.table.y.text = FALSE)

• By looking at the curves, we can see that Blood type B patients have the lowest overall survival probability.
• The risk table below the plot shows us the number of patients at risk at certain time points, which can give an idea of how much information is being used to calculate the estimates at each time

• By looking at the curves, we can see that Male patients have the lowest overall survival probability.

We may want to quantify an effect size for a single variable, or include more than one variable into a regression model to account for the effects of multiple variables.

The Cox regression model is a semi-parametric model that can be used to fit univariable and multivariable regression models that have survival outcomes.

Some key assumptions of the model:

• non-informative censoring
• proportional hazards

We can fit regression models for survival data using the survival::coxph function, which takes a survival::Surv object on the left hand side and has standard syntax for regression formulas in R on the right hand side.

survival::coxph(survival::Surv(futime, delta) ~ factor(abo), data = transplant)

We can see a tidy version of the output using the tidy function from the broom package:

broom::tidy(survival::coxph(survival::Surv(futime, delta) ~ factor(abo), 
                            data = transplant))

## # A tibble: 3 x 7
##   term          estimate std.error statistic p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 factor(abo)B    0.183      0.386    0.474    0.636   -0.573     0.938
## 2 factor(abo)AB   0.195      0.618    0.316    0.752   -1.02      1.41 
## 3 factor(abo)O   -0.0110     0.284   -0.0388   0.969   -0.568     0.546

The quantity of interest from a Cox regression model is a hazard ratio (HR).

If you have a regression parameter \(\beta\) (from column estimate in our survival::coxph) then HR = \(\exp(\beta)\).

For example, from our example we obtain the regression parameter \(\beta_1=0.195\) for AB vs B blood type, so we have HR = \(\exp(\beta_1)=1.21\).

A HR < 1 indicates reduced hazard of death whereas a HR > 1 indicates an increased hazard of death.

So we would say that a person with AB blood type has 1.21 times increased hazard of death as compared to a person with A blood type.

• Time-to-event data is common
• Survival analysis techniques are required to account for censored data
• The survival package provides tools for survival analysis, including the Surv and survfit functions
• The survminer package allows for customization of Kaplan-Meier plots based on ggplot2
• Between-group comparisons can be made with the log-rank test using survival::survdiff
• Multiavariable Cox regression analysis can be accomplished using survival::coxph

# First let's import the data 
q = read.csv('https://3jd8gl2iires146kaw2hgqy9-wpengine.netdna-ssl.com/wp-content/uploads/2018/01/turnover.csv', header = TRUE, sep = ",", na.strings = c("",NA))

Our case uses data of 1785 employees.

$ exp - length of employment in the company

$ event - event (1 - terminated, 0 - currently employed)

$ branch - branch

$ pipeline - source of recruitment

Please note that the data is already prepared for survival analysis. Moreover, length of employment is counted in months up to two decimal places, according to the following formula: (date fire - date hire) / (365.25 / 12).

An accelerated failure-time (AFT) model is a parametric model with covariates and failure times following the survival function of the form \(S(x|Z) = S_0 (x exp[\theta Z])\), where \(S_0\) is a function for the baseline survival rate. The term \(exp[\theta Z]\) is called the acceleration factor.

The AFT model uses covariates to place individuals on different time scales { note the scaling by the covariates in \(S(t|Z)\) via \(exp[\theta Z]\). The AFT model can be rewritten in a log-linear form, where the log of failure time (call this logX) is linearly related to the mean \(\mu\), the acceleration factor, and an error term \(\sigma W\): \[log X = \theta'Z+\sigma W\]

There are serveral AFT models depending on the assumption about the distribution of the error \(W\).

The Weibull distribution has a survival function equal to \[S(t)=exp(-(\lambda t)^\rho)\]

and a hazard function equal to \[\Lambda(t)=\rho \lambda(\lambda t)^{\rho-1}\]

where \(\lambda>0\) and \(\rho>0\). A special case of the Weibul functions is the exponential distribution when \(\rho=1\).

If \(\rho > 1\), then the risk increases over time.

If \(\rho < 1\), then the risk decreases over time.

The exponential model assumes a constant effect of time.

We can perform a likelihood ratio test between the exponential and Weibul model to see if AFT is neccessary.

library(lmtest)
lrtest(w2,w3)

## Likelihood ratio test
## 
## Model 1: Surv(exp, event) ~ branch + pipeline
## Model 2: Surv(exp, event) ~ branch + pipeline
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1  10 -2712.4                        
## 2   9 -2721.2 -1 17.63  2.683e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Although there is little difference in the parameter estimates, the likelihood ratio test suggest we should use the Weibul model over the exponential.

• Machine learning methods are improving on survival models
• We have only discussed X's that are constant over time, but what about time-varying covariates (or worse yet time-dependent covariates)
• There are also issues of competing risks