Survival Analysis

Survival analysis

So far, have seen:
- response variable counted or measured (regression)
- response variable categorized (logistic regression)
But what if response is time until event (eg. time of survival after surgery)?
Additional complication: event might not have happened at end of study (eg. patient still alive). But knowing that patient has “not died yet” presumably informative. Such data called censored.

… continued

Enter survival analysis, in particular the “Cox proportional hazards model”.
Explanatory variables in this context often called covariates.

Packages

Install package survival if not done. Also use broom and marginaleffects from earlier.

library(tidyverse)
library(survival)
library(broom)
library(marginaleffects)

Content warning

In survival analysis, the “event” of interest is often death.
Many of the studies are of cancers of various types and treatments for them.

Example: cancer of the larynx

90 males were diagnosed with cancer of the larynx during the period 1970-1978 at a Dutch hospital.
Times in time are the intervals (in years) between first treatment and either death or the end of the study (January 1, 1983).
Column status records whether the patient was dead or alive at the end of the study.
Also recorded are the patient’s age at the time of diagnosis and the stage of the patient’s cancer.
The stage is one of A through D, with A being least serious and D being most serious.
Predict time until death from stage and age.

The data

my_url <- "http://datafiles.ritsokiguess.site/larynx.csv"
larynx <- read_csv(my_url)
larynx

Effect of age at diagnosis on survival time?

ggplot(larynx, aes(x = age, y = time, colour = status)) + 
  geom_point()

Comments

Not much effect of age on survival time
- even if you focus only on the patients that died (blue points)
What about effect of stage (categorical)?

ggplot(larynx, aes(x = stage, y = time, fill = status)) + 
  geom_boxplot()

Effect of stage on survival time

Comments

Focus on patients that died (blue boxes)
The later the stage, the less the survival time
Better: analysis that incorporates information from patients that survived as well.

Response variable

to allow for patients that survived to contribute to analysis, make response variable out of both time and status:

with(larynx, Surv(time, status == "dead"))

 [1]  0.6   1.3   2.4   2.5+  3.2   3.2+  3.3   3.3+  3.5 
[10]  3.5   4.0   4.0   4.3   4.5+  4.5+  5.3   5.5+  5.9+
[19]  5.9+  6.0   6.1+  6.2+  6.4   6.5   6.5+  6.7+  7.0+
[28]  7.4   7.4+  8.1+  8.1+  9.6+ 10.7+  0.2   1.8   2.0 
[37]  2.2+  2.6+  3.3+  3.6   3.6+  4.0   4.3+  4.3+  5.0+
[46]  6.2   7.0   7.5+  7.6+  9.3+  0.3   0.3   0.5   0.7 
[55]  0.8   1.0   1.3   1.6   1.8   1.9   1.9   3.2   3.5 
[64]  3.7+  4.5+  4.8+  4.8+  5.0   5.0+  5.1+  6.3   6.4 
[73]  6.5+  7.8   8.0+  9.3+ 10.1+  0.1   0.3   0.4   0.8 
[82]  0.8   1.0   1.5   2.0   2.3   2.9+  3.6   3.8   4.3+

Patients with a + by survival time were still alive at end.

Fit Cox model

larynx.1 <- coxph(Surv(time, status == "dead") ~ age + stage, 
                  data = larynx)

stage is categorical with 4 levels, so to see what we can remove, use drop1:

drop1(larynx.1, test = "Chisq")

Comments

Age is not significant (has nothing to add to prediction over and above stage)
There is definitely a significant effect of stage.
Remove age and re-fit:

larynx.2 <- coxph(Surv(time, status == "dead") ~ stage, 
                  data = larynx)
drop1(larynx.2, test = "Chisq")

Understanding the `stage` effect

summary(larynx.2)

Call:
coxph(formula = Surv(time, status == "dead") ~ stage, data = larynx)

  n= 90, number of events= 50 

          coef exp(coef) se(coef)     z Pr(>|z|)    
stageB 0.06481   1.06696  0.45843 0.141   0.8876    
stageC 0.61481   1.84930  0.35519 1.731   0.0835 .  
stageD 1.73490   5.66838  0.41939 4.137 3.52e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

       exp(coef) exp(-coef) lower .95 upper .95
stageB     1.067     0.9372    0.4344      2.62
stageC     1.849     0.5407    0.9219      3.71
stageD     5.668     0.1764    2.4916     12.90

Concordance= 0.668  (se = 0.037 )
Likelihood ratio test= 16.49  on 3 df,   p=9e-04
Wald test            = 19.24  on 3 df,   p=2e-04
Score (logrank) test = 22.88  on 3 df,   p=4e-05

Comments

Numbers in coef column like regression slopes
For a quantitative variable, a more positive coefficient means that for a larger value, the event is more likely to happen sooner
For a categorical variable, a more positive coefficient means that compared to the baseline, the event is more likely to happen sooner.
Here, the more advanced the stage of the cancer, the more likely the patient is to die sooner: the coefficients become more positive for later stages.

Background

All depends on hazard rate, which is based on probability that event happens in the next short time period, given that event has not happened yet:
\(X\) denotes time to event, \(\delta\) is small time interval:
\(h(t) = P(X \le t + \delta | X \ge t) / \delta\)
if \(h(t)\) large, event likely to happen soon (lifetime short)
if \(h(t)\) small, event unlikely to happen soon (lifetime long).

Modelling lifetime

want to model hazard rate
but hazard rate always positive, so actually model log of hazard rate
modelling how (log-)hazard rate depends on other things eg \(X_1 =\) age, \(X_2 =\) treatment, with the \(\beta\) being regression coefficients:
Cox model \(h(t)=h_0(t)\exp(\beta_0+\beta_1X_1+\beta_2 X_2+\cdots)\), or:
\(\log(h(t)) = \log(h_0(t)) + \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots\)
like a generalized linear model with log link.

Plotting predictions

Uses plot_predictions from marginaleffects.
Recall: in condition, inputs are (in order):
- \(x\)-axis
- colours
- facets
In this kind of analysis, the \(x\)-axis is always time.
- hence put variable of interest second in condition:

plot_predictions(larynx.2, condition = c("time", "stage"),
                 type = "survival", vcov = FALSE)

The plot

Comments

graph shows the probability of surviving (event does not happen) over time
the more time passes, the more likely the event (death) is to happen, regardless of stage of cancer
- so all the trends go downhill
stages A and B are very similar in terms of likely time to death
but stage D has the worst survival (most likely to die soon)
- consistent with coefficients on summary output.

Effect of age should be small

go back to first model (with age in it)
not much difference in survival among ages (older age slightly worse)

Numerical predictions

Graphs above based on predicted survival probability over time for each stage
to get numerical predictions, make a dataframe new with explanatory variable values and times to predict for
let`s do 2, 4, 6 years for stages C and D:

new <- datagrid(model = larynx.2, time = c(2, 4, 6), 
                stage = c("C", "D"))
new

The predictions

cbind(predictions(larynx.2, newdata = new,
                  type = "survival")) %>% 
  select(time, stage, estimate)

The predicted survival probabilities go down over time and are lower for stage D than for stage C.

Example 2: kidney data

People whose kidneys do not work properly may have to undergo dialysis.
Procedure where blood is extracted, filtered, and returned to the body.
Sometimes a catheter is used: a soft tube placed in a large vein.
The catheter can be inserted surgically (with an operation) or percutaneously (through the skin).
The exit site (where the blood comes out) can get infected.
Does the time to infection depend on how the catheter was inserted?

The data

The kidney dataset from the KMsurv package (install first):

library(KMsurv)
data(kidney)
kidney %>% slice_sample(n = 10)

Data processing

Time to infection or to censoring is in months
Numerical codes:
- delta: 1 = infected, 0 not (censored, usually from catheter failure)
- type: 1 = surgical, 2 = percutaneous
Turn numerical codes into something more descriptive.

Data processing cont’d

kidney %>% 
  mutate(infected = recode_values(
    delta,
    0 ~ "no",
    1 ~ "yes"
  )) %>% 
  mutate(type = recode_values(
    type,
    1 ~ "surgical",
    2 ~ "percutaneous"
  )) -> kidney

What we have now (sample of rows)

kidney %>% slice_sample(n = 10)

Fit Cox model

kidney.1 <- coxph(Surv(time, infected == "yes") ~ type, 
                  data = kidney)
drop1(kidney.1, test = "Chisq")

No significant difference between catheter insertion types in terms of time to infection
But… model assumes that proportional hazards is correct, and we need a way to test this.
Function cox.zph. Gives a P-value, and also a graph to help understand why proportional hazards fails (if it does).

Test for proportional hazards

kidney.1.zph <- cox.zph(kidney.1)
kidney.1.zph

       chisq df      p
type    8.62  1 0.0033
GLOBAL  8.62  1 0.0033

Proportional hazards is rejected (P-value 0.0033).

A plot

plot(kidney.1.zph)

Comments

If proportional hazards correct, trace on the plot should go straight across at the coefficient for type in the model, here 0.6:

tidy(kidney.1)

For most times, coefficient is positive (infection more likely to happen sooner with surgical vs percutaneous)
But for small times, coefficient is negative (infection less likely to happen sooner with surgical vs percutaneous)
Under proportional hazards, coefficient cannot depend on time.
Need a “time-dependent covariate” model (beyond our scope).

Test larynx model for proportional hazards

We don’t (yet) know that proportional hazards is good for the larynx cancer data
Explanatory variable stage has 4 levels, so get test and graph for each one vs baseline A:

larynx.2.zph <- cox.zph(larynx.2, terms = FALSE)
larynx.2.zph

         chisq df     p
stageB 0.18227  1 0.669
stageC 3.26317  1 0.071
stageD 0.00393  1 0.950
GLOBAL 3.58948  3 0.309

No problems with proportional hazards here, though stage C is on the edge.

The graphs

par(mfrow = c(2, 2)) # make a 2x2 grid of base R plots
par(mar = c(4, 4, 0.1, 0.1)) # reduce blank areas around them
plot(larynx.2.zph)

Comments

Trends go near enough to straight across, so no time-dependence of hazard.
Stage B (vs A) and D (vs A) definitely have no problems
Stage C (vs A) questionable. This went with P-value 0.071, so not quite bad enough to be significant.

Example 3: lung cancer

survival package has lung: data set measuring survival in patients with advanced lung cancer.
Along with survival time, number of “performance scores” included, measuring how well patients can perform daily activities.
Sometimes high good, but sometimes bad!
Variables below, from the data set help file (?lung).

The variables

Uh oh, missing values

lung %>% select(meal.cal, wt.loss)

A closer look

summary(lung)

      inst            time            status           age       
 Min.   : 1.00   Min.   :   5.0   Min.   :1.000   Min.   :39.00  
 1st Qu.: 3.00   1st Qu.: 166.8   1st Qu.:1.000   1st Qu.:56.00  
 Median :11.00   Median : 255.5   Median :2.000   Median :63.00  
 Mean   :11.09   Mean   : 305.2   Mean   :1.724   Mean   :62.45  
 3rd Qu.:16.00   3rd Qu.: 396.5   3rd Qu.:2.000   3rd Qu.:69.00  
 Max.   :33.00   Max.   :1022.0   Max.   :2.000   Max.   :82.00  
 NAs    :1                                                       
      sex           ph.ecog          ph.karno        pat.karno     
 Min.   :1.000   Min.   :0.0000   Min.   : 50.00   Min.   : 30.00  
 1st Qu.:1.000   1st Qu.:0.0000   1st Qu.: 75.00   1st Qu.: 70.00  
 Median :1.000   Median :1.0000   Median : 80.00   Median : 80.00  
 Mean   :1.395   Mean   :0.9515   Mean   : 81.94   Mean   : 79.96  
 3rd Qu.:2.000   3rd Qu.:1.0000   3rd Qu.: 90.00   3rd Qu.: 90.00  
 Max.   :2.000   Max.   :3.0000   Max.   :100.00   Max.   :100.00  
                 NAs    :1        NAs    :1        NAs    :3       
    meal.cal         wt.loss       
 Min.   :  96.0   Min.   :-24.000  
 1st Qu.: 635.0   1st Qu.:  0.000  
 Median : 975.0   Median :  7.000  
 Mean   : 928.8   Mean   :  9.832  
 3rd Qu.:1150.0   3rd Qu.: 15.750  
 Max.   :2600.0   Max.   : 68.000  
 NAs    :47       NAs    :14

Remove obs with any missing values

lung %>% drop_na() -> lung.complete
lung.complete %>%
  select(meal.cal:wt.loss) %>%
  slice(1:10)

Check!

summary(lung.complete)

      inst            time            status           age       
 Min.   : 1.00   Min.   :   5.0   Min.   :1.000   Min.   :39.00  
 1st Qu.: 3.00   1st Qu.: 174.5   1st Qu.:1.000   1st Qu.:57.00  
 Median :11.00   Median : 268.0   Median :2.000   Median :64.00  
 Mean   :10.71   Mean   : 309.9   Mean   :1.719   Mean   :62.57  
 3rd Qu.:15.00   3rd Qu.: 419.5   3rd Qu.:2.000   3rd Qu.:70.00  
 Max.   :32.00   Max.   :1022.0   Max.   :2.000   Max.   :82.00  
      sex           ph.ecog          ph.karno        pat.karno     
 Min.   :1.000   Min.   :0.0000   Min.   : 50.00   Min.   : 30.00  
 1st Qu.:1.000   1st Qu.:0.0000   1st Qu.: 70.00   1st Qu.: 70.00  
 Median :1.000   Median :1.0000   Median : 80.00   Median : 80.00  
 Mean   :1.383   Mean   :0.9581   Mean   : 82.04   Mean   : 79.58  
 3rd Qu.:2.000   3rd Qu.:1.0000   3rd Qu.: 90.00   3rd Qu.: 90.00  
 Max.   :2.000   Max.   :3.0000   Max.   :100.00   Max.   :100.00  
    meal.cal         wt.loss       
 Min.   :  96.0   Min.   :-24.000  
 1st Qu.: 619.0   1st Qu.:  0.000  
 Median : 975.0   Median :  7.000  
 Mean   : 929.1   Mean   :  9.719  
 3rd Qu.:1162.5   3rd Qu.: 15.000  
 Max.   :2600.0   Max.   : 68.000

No missing values left.

Recoding

All these variables are recorded as numbers, though status, sex are actually categorical:

lung.complete %>% 
  mutate(status = recode_values(status,
                                1 ~ "censored",
                                2 ~ "dead")) %>% 
  mutate(sex = recode_values(sex, 
                             1 ~ "male",
                             2 ~ "female")) -> lung.complete

Check recoding

lung.complete %>% select(-inst) %>% slice_sample(n = 15)

Model 1: use everything except `inst`

names(lung.complete)

 [1] "inst"      "time"      "status"    "age"       "sex"      
 [6] "ph.ecog"   "ph.karno"  "pat.karno" "meal.cal"  "wt.loss"

Event was status being dead:

lung.1 <- coxph(
  Surv(time, status == "dead") ~ . - inst - time - status,
  data = lung.complete
)

“Dot” means “all the other variables”; list ones to explicitly exclude.

Model 1

tidy(lung.1)

Some of these are definitely not significant:

sex and ph.ecog definitely significant here
age, pat.karno and meal.cal definitely not

Using `step`

(Output is long, suppressed)

lung.2 <- step(lung.1, test = "Chisq")

What is left

tidy(lung.2)

Fit model 3

Try taking out those marginally significant ones:

lung.3 <- update(lung.2, . ~ . - ph.karno - pat.karno - wt.loss)

and test whether that was sensible:

anova(lung.3, lung.2)

Only just.

Look at model 3

summary(lung.3)

Call:
coxph(formula = Surv(time, status == "dead") ~ sex + ph.ecog, 
    data = lung.complete)

  n= 167, number of events= 120 

          coef exp(coef) se(coef)     z Pr(>|z|)    
sexmale 0.5101    1.6655   0.1969 2.591 0.009579 ** 
ph.ecog 0.4825    1.6201   0.1323 3.647 0.000266 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

        exp(coef) exp(-coef) lower .95 upper .95
sexmale     1.665     0.6004     1.132      2.45
ph.ecog     1.620     0.6172     1.250      2.10

Concordance= 0.641  (se = 0.031 )
Likelihood ratio test= 19.48  on 2 df,   p=6e-05
Wald test            = 19.35  on 2 df,   p=6e-05
Score (logrank) test = 19.62  on 2 df,   p=5e-05

Commentary

Both remaining variables strongly significant.
Nature of effect on survival time? Consider later.
Picture?

Plotting survival probabilities

Assess (separately) the effect of sex and ph.ecog score using plot_predictions
Don’t forget to add time (here actually called time) to the condition.

Effect of `sex`:

plot_predictions(lung.3, condition = c("time", "sex"), 
                 type = "survival")

Females (sex = 2) have better survival than males.
This graph from a mean ph.ecog score, but the male-female comparison is the same for any score.

Effect of `ph.ecog` score:

plot_predictions(lung.3, condition = c("time", "ph.ecog"), 
                 type = "survival")

Comments

A lower ph.ecog score is better.
For example, a patient with a score of 0 has almost a 50-50 chance of living 500 days, but a patient with a score of 3 has almost no chance to survive that long.
Is this for males or females? See over. (The comparison of scores is the same for both.) How many males and females did we observe?

lung.complete %>% count(sex)

more males, so ph.ecog graph was for males.

Sex and `ph.ecog` score

plot_predictions(lung.3, condition = c("time", "ph.ecog", "sex"), type = "survival")

Comments

The previous graph was (closer to) males. There were more males in the dataset (sex of 1).
This pair of graphs shows the effect of ph.ecog score (above and below on each facet), and the effect of males (left) vs. females (right).
The difference between males and females is about the same as 1 point on the ph.ecog scale (compare the red curve on the left facet with the green curve on the right facet).

A second go

Can we put the male and female survival curves on the same plot?

Try building plot ourselves:

plot_predictions(lung.3, condition = c("time", 
                                       "ph.ecog", "sex"), 
                 type = "survival", draw = FALSE) %>% 
  select(estimate, time, ph.ecog, sex) -> dd
dd

The plot, setup

ggplot(dd, aes(x = time, y = estimate, 
               colour = ph.ecog, 
               linetype = sex)) + geom_line()

The plot

The summary again

summary(lung.3)

Call:
coxph(formula = Surv(time, status == "dead") ~ sex + ph.ecog, 
    data = lung.complete)

  n= 167, number of events= 120 

          coef exp(coef) se(coef)     z Pr(>|z|)    
sexmale 0.5101    1.6655   0.1969 2.591 0.009579 ** 
ph.ecog 0.4825    1.6201   0.1323 3.647 0.000266 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

        exp(coef) exp(-coef) lower .95 upper .95
sexmale     1.665     0.6004     1.132      2.45
ph.ecog     1.620     0.6172     1.250      2.10

Concordance= 0.641  (se = 0.031 )
Likelihood ratio test= 19.48  on 2 df,   p=6e-05
Wald test            = 19.35  on 2 df,   p=6e-05
Score (logrank) test = 19.62  on 2 df,   p=5e-05

Comments

Males have a higher hazard of death than (baseline) females (positive coef). That is, females are more likely to survive longer than males, all else (ph.ecog score) equal.
A higher ph.ecog score goes with a higher hazard of death (positive coef). So patients with a lower score are more likely to survive longer, all else (sex) equal.
Consistent with graphs we drew.

Assess proportionality of hazards

Explanatory variables are numbers (ph.ecog) and categorical with two levels (sex):

lung.3.zph <- cox.zph(lung.3)
lung.3.zph

        chisq df     p
sex      1.24  1 0.266
ph.ecog  5.27  1 0.022
GLOBAL   6.19  2 0.045

ph.ecog does not show proportional hazards.

Two graphs, one for each explanatory variable:

par(mfrow = c(1, 2)) # make 1x2 grid of base R plots
par(mar = c(4, 4, 0.1, 0.1)) # reduce blank areas around them
plot(lung.3.zph)

Comments on graphs

sex graph (left) goes more or less straight across (good).
ph.ecog graph (right) seems to dip after about 300 days, so hazard of a higher ph.ecog score is higher earlier (not constant over time).
A better analysis would treat ph.ecog score as having a varying effect over time (for example, it only has an effect up to 300 days).

Survival Analysis