--- title: "Survival Analysis" editor: markdown: wrap: 72 --- ## 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. ```{r} 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 \small ```{r} my_url <- "http://datafiles.ritsokiguess.site/larynx.csv" larynx <- read_csv(my_url) larynx ``` \normalsize ## Effect of age at diagnosis on survival time? ```{r} #| fig-height: 4.5 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)? ```{r} #| eval: false ggplot(larynx, aes(x = stage, y = time, fill = status)) + geom_boxplot() ``` ## Effect of stage on survival time ```{r} #| echo: false ggplot(larynx, aes(x = stage, y = time, fill = status)) + geom_boxplot() ``` ## 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`: \small ```{r} #| echo: false w <- getOption("width") options(width = 60) ``` ```{r} with(larynx, Surv(time, status == "dead")) ``` ```{r} #| echo: false options(width = w) ``` \normalsize Patients with a `+` by survival time were still alive at end. ## Fit Cox model ```{r} 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`: \small ```{r} drop1(larynx.1, test = "Chisq") ``` \normalsize ## 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: \small ```{r} larynx.2 <- coxph(Surv(time, status == "dead") ~ stage, data = larynx) drop1(larynx.2, test = "Chisq") ``` \normalsize ## Understanding the `stage` effect \footnotesize ```{r} summary(larynx.2) ``` \normalsize ## 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`: ```{r} #| eval: false plot_predictions(larynx.2, condition = c("time", "stage"), type = "survival", vcov = FALSE) ``` ## The plot ```{r} #| echo: false plot_predictions(larynx.2, condition = c("time", "stage"), type = "survival", vcov = FALSE) ``` ## 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 ```{r} #| echo: false #| fig-height: 5 plot_predictions(larynx.1, condition = c("time", "age"), type = "survival", vcov = FALSE) ``` - 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: \small ```{r} new <- datagrid(model = larynx.2, time = c(2, 4, 6), stage = c("C", "D")) new ``` \normalsize ## The predictions ```{r} 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): ```{r} 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 ```{r} 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) ```{r} kidney %>% slice_sample(n = 10) ``` ## Fit Cox model \small ```{r} kidney.1 <- coxph(Surv(time, infected == "yes") ~ type, data = kidney) drop1(kidney.1, test = "Chisq") ``` \normalsize - 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 ```{r} kidney.1.zph <- cox.zph(kidney.1) kidney.1.zph ``` Proportional hazards is rejected (P-value 0.0033). ## A plot ```{r} #| fig-height: 5 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: ```{r} 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: ```{r} larynx.2.zph <- cox.zph(larynx.2, terms = FALSE) larynx.2.zph ``` - No problems with proportional hazards here, though stage C is on the edge. ## The graphs ```{r} #| fig-height: 5 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 ![](lung-cancer-data.png) ## Uh oh, missing values \scriptsize ```{r bSurvival-13} lung %>% select(meal.cal, wt.loss) ``` \normalsize ## A closer look ```{r bSurvival-14, echo=F} options(width = 70) ``` \tiny ```{r bSurvival-15} summary(lung) ``` \normalsize ## Remove obs with *any* missing values ```{r bSurvival-16} lung %>% drop_na() -> lung.complete lung.complete %>% select(meal.cal:wt.loss) %>% slice(1:10) ``` ## Check! \tiny ```{r bSurvival-17} summary(lung.complete) ``` \normalsize No missing values left. ## Recoding All these variables are recorded as numbers, though `status`, `sex` are actually categorical: ```{r} lung.complete %>% mutate(status = recode_values(status, 1 ~ "censored", 2 ~ "dead")) %>% mutate(sex = recode_values(sex, 1 ~ "male", 2 ~ "female")) -> lung.complete ``` ## Check recoding ```{r} #| echo: false w <- getOption("width") options(width = 110) ``` \footnotesize ```{r} lung.complete %>% select(-inst) %>% slice_sample(n = 15) ``` ```{r} #| echo: false options(width = w) ``` \normalsize ## Model 1: use everything except `inst` \footnotesize ```{r bSurvival-18} names(lung.complete) ``` \normalsize - Event was `status` being `dead`: ```{r bSurvival-19 } 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 ```{r bSurvival-20} 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) ```{r} #| results: hide lung.2 <- step(lung.1, test = "Chisq") ``` ## What is left ```{r} tidy(lung.2) ``` ## Fit model 3 Try taking out those marginally significant ones: \small ```{r} lung.3 <- update(lung.2, . ~ . - ph.karno - pat.karno - wt.loss) ``` \normalsize and test whether that was sensible: \small ```{r} anova(lung.3, lung.2) ``` \normalsize Only just. ## Look at model 3 ```{r} summary(lung.3) ``` ## 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`: ```{r} #| fig.height: 4 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: ```{r} #| fig.height: 5 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? ```{r} lung.complete %>% count(sex) ``` - more males, so `ph.ecog` graph was for males. ## Sex and `ph.ecog` score ```{r} #| fig.height = 4.5 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: ```{r} 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 ```{r} #| eval: false ggplot(dd, aes(x = time, y = estimate, colour = ph.ecog, linetype = sex)) + geom_line() ``` ## The plot ```{r} #| echo: false ggplot(dd, aes(x = time, y = estimate, colour = ph.ecog, linetype = sex)) + geom_line() ``` ## The summary again \small ```{r} summary(lung.3) ``` \normalsize ## 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`): ```{r} lung.3.zph <- cox.zph(lung.3) lung.3.zph ``` `ph.ecog` does not show proportional hazards. ## Two graphs, one for each explanatory variable: ```{r} #| fig-height: 5 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).