Chapter 9 Survival analysis
A common form of data in clinical trials is survival (time to event) data. Survival outcomes can often be considered as continuous, but are usually quite skewed. The issue of censoring requires special techniques for statistical analysis.
9.1 Analysis of survival outcomes
9.1.1 Censoring
A survival time \(T\) is called (right-)censored, if it is known that the individual has survived up to time \(T\), but whose survival status past that point is not known. For example, the study might be finished before all patients have died, a patient might decide to withdraw from the study or a patient might die from some other reason, for example in a car crash.
A central assumption in survival analysis is that the censoring time is independent from the survival time (non-informative censoring).
Example 9.1 The Chemotherapy-Study is a randomized study to compare two treatments of head and neck cancer on \(224\) patients: Combination of radiotherapy and chemotherapy (RT+CT) (\(n=112\)) versus radiotherapy alone (RT) (\(n=112\)). The outcome is survival time (in years) from start of treatment. Figure 9.1 shows data for the survival times of the first 10 patients. The survival time is measured from the start of treatment (starting point \(0\)) and it ends if a patient died or is censored.
Figure 9.1: Survival times for a selection of patients from Example 9.1
Table 9.1 shows a selection of variables in the data for the first five patients. The notion overall survival (OS) denotes survival without any further conditions (in contrast to, , progression free survival). In this example, OS is the survival time in years and the \(+\) stands for “censored”. The format of the variable OS (with \(+\) for censored) is the output of the R function Surv(). It combines three variables from the data: an event indicator, time to event and time to last follow-up. Patients with no time to event are censored at their last follow-up. The function Surv() needs the times (combined in one variable) and the events for each patient.
| ID | Sex | Age | Treatment | LK | Stage | PS | OS | |
|---|---|---|---|---|---|---|---|---|
| 1 | 1 | Female | 57.4 | RT+CT | 2c | 4 | 1 | 15.3+ |
| 5 | 2 | Male | 66.9 | RT+CT | 0 | 3 | 1 | 10.5 |
| 9 | 3 | Male | 68.3 | RT+CT | 0 | 3 | 1 | 7.9+ |
| 10 | 4 | Male | 51.6 | RT | 0 | 2 | 0 | 5.4 |
| 2 | 5 | Male | 38.8 | RT | 2c | 4 | 0 | 3.1 |
9.1.2 Survival function and curve
The survival function
\[\begin{equation*} S(t) = \Pr(T > t) \end{equation*}\]
gives the probability of survival at least up to time \(t\). The survival curve is a plot of \(S(t)\) versus \(t\), see Figure 9.2.
Figure 9.2: Examples of survival curves for three different survival functions.
9.1.2.1 Life table method
The life table method is a classical approach to estimate the survival function \(S(t)\). It is based on a partition of the time axis into time intervals \(i=1,\ldots,I\) of typically equal length. The risk of death \(r_i\) in interval \(i\) is estimated as
\[\begin{equation*} \hat r_i = \frac{\mbox{\# patients which died in interval $i$}}{\mbox{\# patients at risk in interval $i$}} = \frac{d_i}{n_i} \end{equation*}\]
By convention, observations censored in interval \(i\) contribute half to the number of patients at risk. Note that the denominator \(n_i\) is monotonically decreasing as a function of time.
The survival function is then estimated as
\[\begin{equation*} \hat S(i) = \prod_{t=1} ^i(1 - \hat r_t) = (1 - \hat r_1) \times \ldots \times (1 - \hat r_i). \end{equation*}\]
Wald confidence intervals can be calculated for \(S(i)\), \(\log \hat S(i)\) or \(\log \{-\log \hat S(i)\}\) with appropriate back-transformation. The corresponding standard errors can be found in textbooks, Collett (2003) (Section 2.2).
Example 9.2 Table 9.3 shows the estimated survival function (values for each interval/year) in the RT+CT group from Example 9.1. Figure 9.3 shows the estimated survival curve with 95% confidence intervals for every year (in red).
| Interval (year) | # patients at beginning of interval | # deaths during interval | # censored | # patients at risk | \(r_i\) | \(1-r_i\) | \(\hat S(i)\) |
| 1 | 112 | 25 | 0 | 112 | 0.22 | 0.78 | 0.78 |
| 2 | 87 | 14 | 1 | 86.5 | 0.16 | 0.84 | 0.65 |
| 3 | 72 | 10 | 1 | 71.5 | 0.14 | 0.86 | 0.56 |
| 4 | 61 | 8 | 1 | 60.5 | 0.13 | 0.87 | 0.49 |
| 5 | 52 | 4 | 4 | 50 | 0.08 | 0.92 | 0.45 |
| 6 | 44 | 3 | 5 | 41.5 | 0.07 | 0.93 | 0.41 |
| 7 | 36 | 4 | 3 | 34.5 | 0.12 | 0.88 | 0.37 |
| 8 | 29 | 4 | 3 | 27.5 | 0.15 | 0.85 | 0.31 |
| 9 | 22 | 1 | 3 | 20.5 | 0.05 | 0.95 | 0.3 |
| 10 | 18 | 1 | 5 | 15.5 | 0.06 | 0.94 | 0.28 |
| 11 | 12 | 2 | 0 | 12 | 0.17 | 0.83 | 0.23 |
| 12 | 10 | 0 | 1 | 9.5 | 0 | 1 | 0.23 |
| 13 | 9 | 2 | 2 | 8 | 0.25 | 0.75 | 0.17 |
| 14 | 5 | 1 | 2 | 4 | 0.25 | 0.75 | 0.13 |
| 15 | 2 | 0 | 1 | 1.5 | 0 | 1 | 0.13 |
| 16 | 1 | 0 | 1 | 0.5 | 0 | 1 | 0.13 |
9.1.2.2 Kaplan-Meier estimate
In many studies, the exact follow-up time for each individual is known. Aggregation in arbitrary time intervals can then be avoided. The Kaplan-Meier estimate of the survival function uses intervals which contain only one event, so \(d_i \in \{0, 1\}\) (if there are no ties). Application of the life table method to these intervals yields the Kaplan-Meier estimate, a step function with jumps of relative size \[ 1/\mbox{number of subjects at risk} \] at the time of each event.
There are some critical assumptions in Kaplan-Meier estimation:
- At any time, patients who are censored have the same survival prospects as those who continue to be followed. If for some patients censoring is related to survival, the estimated survival probabilities would be biased.
- Survival probabilities are the same for subjects recruited early and late in the study. In a long term study, the case mix may change over the period of recruitment, or there may be an innovation in ancillary treatment.
- The event happens at the time specified. Could be a problem, for example, if the event were recurrence of a tumor which would be detected at a regular examination (interval censoring).
Example 9.3 Figure 9.3 compares the Kaplan-Meier estimate of the survival curve to the estimate obtained with the life table method using yearly intervals, in Example 9.1 for the RT+CT group. The Kaplan-Meier estimate can be calculated and plotted with the following R code. The option conf.int = TRUE shows a dotted confidence band and the option mark.time = TRUE would show censoring events (see Figure 9.5).
sakk.surv <- survfit(Surv(os, death) ~ 1, data = sakk.rtct,
conf.type = "log-log")
plot(sakk.surv, conf.int = TRUE, mark.time = FALSE, lwd = 2,
col = col.blue, log = FALSE, xlab = "Survival time (in years)",
ylab = "Survival probability", xlim = c(0, Imax))
Figure 9.3: Comparison of Kaplan-Meier method (blue) with life table method (red, with yearly 95% confidence intervals) to estimate the survival curve in the RT+CT group from Example 9.1.
With survminer::ggsurvplot, a Kaplan-Meier estimate with additional risk table can be obtained as shown in Figure 9.4.
Figure 9.4: Kaplan-Meier estimate of the survival curve in the RT+CT group from Example 9.1, with risk table.
We need to be careful with censored data. Censored observations cannot be simply ignored, i.e. excluded, and they cannot be treated as non-censored. Doing so would lead to biased estimation of the survival function as illustrated by Figure 9.5. Excluding censored observations means losing data on patients without events, so will underestimate the survival curve. Treating them as non-censored, as death events, also leads to underestimation of the survival curve.
Figure 9.5: Kaplan-Meier curves in the RT+CT group from Example 9.1 where censored observations are ignored (top right) and where censored observations are treated as non-censored (bottom right). On the left side (top and bottom) the correct Kaplan-Meier estimate is shown, the dashes indicate censoring events.
9.1.3 Median survival time
The median survival time \(t_{\small{\mbox{Med}}}\) (and other quantiles) can be easily read off the Kaplan-Meier estimate. Also a confidence interval can be derived in this way. The ``square-and-add’’ method, see Appendix A.2.3, can be used to compute a confidence interval for the difference in median survival time between two groups.
Example 9.4 Figure 9.6 is a graphical illustration of the median survival time in the RT+CT group from Example 9.1.
Figure 9.6: Median survival time in the RT+CT group from Example 9.1.
| \(t_{\mbox{\tiny Med}}\) | 95%-CI | |||
|---|---|---|---|---|
| 1 | \(n\) | Number of Deaths | (in years) | (in years) |
| group = RT+CT | 112 | 79 | 3.9 | from 2.2 to 6.0 |
| group = RT | 112 | 88 | 2.4 | from 1.8 to 3.8 |
Figure 9.7: Median survival times with 95% confidence intervals of RT and RT+CT groups from Example 9.1.
The difference in median survival time is 1.5 years (95%-CI: from \(-0.7\) to \(3.7\) years). Both medians with CI are shown in Table 9.4 and graphically in Figure 9.7.
9.2 Comparison of survival curves
We have to be careful about what groups we can compare. Breslow’s first rule of survival analysis is that ’’Selection into the study cohort, or into subgroups to be compared in the analysis, must not depend on events that occur after the start of follow-up``. Survival by tumor response analysis in oncology is an example that leads to biased estimates of the survival distributions, invalid statistical tests, and misleading conclusions. This bias results from the fact that responders must live long enough for response to be observed; there is no such requirement for non-responders (Anderson, Cain, and Gelber 1983, Anderson2008).
Example 9.5 Another example where Breslow’s rule is violated is a study reporting that Oscar winners live longer (Redelmeier and Singh 2001). As other scientists have noticed, this study was subject to survival bias or immortal time bias (Sylvestre, Huszti, and Hanley 2006). See Figure 9.8 for abstracts of both publications.
Figure 9.8: Study about survival in academy award-winning actors and actresses and a reanalysis thereof that raises awareness for immortal time bias.
We now discuss how survival curves can be compared between groups. For example, the survival curves of the RT and RT+CT groups shown in Figure 9.9. Later, we show R code how to generate such a plot with function ggsurvplot.
Figure 9.9: Kaplan-Meier curves for both groups in Example 9.1.
9.2.1 Log-rank test
The log-rank test, also called Mantel-Cox test, quantifies the evidence against the null hypothesis of equal survival curves \(S_A(t) = S_B(t)\) for all times \(t\) with a \(p\)-value. The method compares the observed number of events with the corresponding expected number of events (under the null hypothesis \(H_0\)) in each group.
## Call:
## survdiff(formula = Surv(time, death) ~ treat, data = sakk)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## treat=RT+CT 112 79 87.5 0.822 1.74
## treat=RT 112 88 79.5 0.904 1.74
##
## Chisq= 1.7 on 1 degrees of freedom, p= 0.2The two survival curves can be shown in one plot with the following code, see Figure 9.10:
sakk.surv1 <- survfit(Surv(time, death) ~ treat, data = sakk,
conf.type = "log-log")
ggsurvplot(sakk.surv1, risk.table = TRUE, conf.int = TRUE, pval = TRUE)
Figure 9.10: Kaplan-Meier curves for both groups in Example 9.1 obtained with ggsurvplot.
9.2.2 Hazard function
We can also use the mortality rate to describe the risk of death in a certain interval:
\[\begin{eqnarray*} \frac{\mbox{Death risk in interval}}{\mbox{Length of interval}} \end{eqnarray*}\]
The mortality rate is conditional on having survived up to the interval of interest. The hazard function is obtained by making the length of the intervals very small:
\[\begin{equation*} h(t)= \lim\limits_{h \rightarrow 0} \frac{\operatorname{\mathsf{Pr}}(t < T < t+h \,\vert\,T>t)}{h} \end{equation*}\]
The hazard function and the survival function have a one-to-one relationship:
\[\begin{equation*} S(t) = \exp\left(- \int_0^t h(u) du \right). \end{equation*}\]
Figure 9.11 shows the survival function once for a constant hazard rate and once for an increasing hazard rate.
Figure 9.11: Survival function and hazard ratio (HR) for different hazard rates: constant (top), linearly increasing (middle) and non-linearly increasing (bottom)
9.2.3 Hazard ratio
Let \(h_A(t)\) and \(h_B(t)\) denote the hazard functions in groups A and B. The proportional hazards assumption implies a time-constant hazard ratio%
\[\begin{equation*} \mbox{HR} = {h_A(t)}/{h_B(t)}. \end{equation*}\]
The hazard ratio can be interpreted as
- relative risk at any given time \(t\):
\[\begin{equation*} \mbox{HR} = \frac{h_A(t)}{h_B(t)}\approx \frac{{\operatorname{\mathsf{Pr}}(t < T_A < t+h \,\vert\,T_A>t)}/{h}}{{\operatorname{\mathsf{Pr}}(t < T_B < t+h \,\vert\,T_B>t)}/{h}} \equiv \frac{{\operatorname{\mathsf{Pr}}(t < T_A < t+h \,\vert\,T_A>t)}}{{\operatorname{\mathsf{Pr}}(t < T_B < t+h \,\vert\,T_B>t)}} \end{equation*}\]
- power transformation of the survival function at any time \(t\): % \[\begin{equation*} S_A(t) = S_B(t)^{\scriptsize \mbox{HR}}. \end{equation*}\]
Figure 9.11 shows two different settings of hazard functions (constant and non-constant) that both have the same hazard ratio of 2.
The hazard ratio can be calculated from the output of the log-rank test. Only the observed and expected number of cases are required.
lrTest <- survdiff(Surv(time, death) ~ treat, data = sakk)
observed <- lrTest$obs
expected <- lrTest$exp
ratio <- observed/expected
HR <- ratio[2]/ratio[1]
SElogHR <- sqrt(sum(1/expected))
printWaldCI(log(HR), SElogHR, FUN = exp)## Effect 95% Confidence Interval P-value
## [1,] 1.225 from 0.904 to 1.660 0.199.2.3.1 The Cox model
The (semiparametric) Cox model assumes that the hazard function \(h_i(t)\) of individual \(i\) with explanatory variables \(\boldsymbol{ x}_i\) has the form % \[\begin{equation*} h_i(t; \boldsymbol{ x}_i) = \exp(\boldsymbol{ x}_i^{\top} \boldsymbol{ \beta}) \cdot h_0(t), \end{equation*}\] % where \(h_0(t)\) is the baseline hazard function. Ratios of hazard functions are time-constant: % \[\begin{equation*} \frac{h_i(t; \boldsymbol{ x}_i)}{h_j(t; \boldsymbol{ x}_j)} = \frac{\exp(\boldsymbol{ x}_i^{\top} \boldsymbol{ \beta})} {\exp(\boldsymbol{ x}_j^{\top} \boldsymbol{ \beta})}, \end{equation*}\] % we have proportional hazards. A profile likelihood approach allows to estimate the log hazard ratios \(\boldsymbol{ \beta}\), no matter which form \(h_0(t)\) has. To estimate the hazard ratio, we use Cox-Regression as in the following example.
Example 9.6 Cox regression in Example 9.1:
model.cox0 <- coxph(Surv(os, death) ~ trt2, data = dat.tmp)
knitr::kable(tableRegression(model.cox0, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| trt2RT | 1.23 | from 0.90 to 1.66 | 0.19 |
The comparison RT vs. RT+CT gives \(\mbox{HR}=1.23\) with relatively large confidence interval. There is no evidence for a treatment effect (\(p=0.19\)).
9.2.3.2 Adjusting in survival analysis
There may be differences between treatment groups in important baseline characteristics. Table 9.6 shows the distribution of gender in the two treatment groups in Example 9.1, which is surprisingly unbalanced (percentage female 10% resp. 21%). It is recommended to adjust the treatment effect for potential confounders at baseline, which can be done with Cox regression.
| RT + CT | RT | |
|---|---|---|
| F | 11 | 23 |
| M | 101 | 89 |
Example 9.7 Adjusting for baseline values with Cox regression in Example 9.1:
model.cox1 <- coxph(Surv(os, death) ~ trt
+ sex, data = dat.tmp)
knitr::kable(tableRegression(model.cox1, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| trt (RT) | 1.26 | from 0.92 to 1.71 | 0.14 |
| sex (F) | 0.81 | from 0.52 to 1.26 | 0.34 |
The adjusted hazard of death (instantaneous death risk) is increased by \(26\)% for RT relative to RT + CT (\(p=0.14\)). Women have a \(19\)% reduced hazard of death (\(p=0.34\)) compared to men. However, there is no evidence for both effects.
9.3 Sample size in survival studies
9.3.1 Total number of events
There is a simple formula for the required total number of events Collett (2003) (Ch. 10),
\[\begin{equation*} d = \frac{ 4 (u + v)^2}{ \Delta^2}, \end{equation*}\]
where
- \(\Delta\) is the log hazard ratio to be detected with
- Type I error rate \(\alpha\) (\(\rightarrow\) \(u = z_{1-\alpha/2}\)) and
- power \(1-\beta\) (\(\rightarrow\) \(v = z_{1-\beta}\)).
It can be calculated with the biostatUZH package:
## [1] 2279.3.2 Total number of patients
The required total number of patients \(n\) can be found from
\[\begin{equation*} n = d / \Pr(\mbox{event}), \end{equation*}\]
where \(\Pr(\mbox{event})\) is calculated based on the assumed survival functions in the two groups and assumptions on
- the accrual period \(a\)
- and the follow-up period \(f\).
Example 9.8 Survival of females in lung cancer data is shown in Figure 9.12.
data(cancer, package = "survival")
cancer.female <- subset(cancer, sex == 2)
survObj <- Surv(time = cancer.female$time,
event = (cancer.female$status == 2))
fit <- survfit(survObj ~ 1)
ggsurvplot(fit, data = cancer.female)
Figure 9.12: Kaplan-Meier estimate of the survival curve of females in lung cancer data.
In the function biostatUZH::sampleSizeSurvival, we need to specify:
- HR: Hazard Ratio
- power: power
- sig.level: significance level
- a.length: length of accrual period
- f.length: length of follow-up period
- dist: “nonp” will take Kaplan-Meier survival curve
- survfit.ref: Kaplan-Meier survival curve
res <- sampleSizeSurvival(HR = 0.65, power=.9, sig.level=0.05,
a.length = 400, f.length = 400,
dist = "non-parametric", survfit.ref = fit)
print(t(res))## n HR power sig.level alternative distribution
## [1,] 379 0.65 0.9 0.05 "two.sided" "non-parametric"
## PrEvent Events
## [1,] 0.5996061 227Instead of a non-parametric fit of the survival curves, the function can also be given a parametric input for the survival curves such as a Weibull or exponential distribution with specified parameters.
9.4 Time-varying explanatory variables
Time-varying explanatory variables, e.g. exposure yes/no, can be dealt with by
- splitting the follow-up into parts with constant covariate values and
- treating the parts as distinct subjects.
This approach requires allowance for late entries and is valid if changes in explanatory variables are unrelated to the underlying risk of failure. We will now see two examples of time-dependent analyses.
Example 9.9 The PBC study is about treatment of primary biliary cirrhosis (PBC) of the liver. We will use the first 312 cases in this data set from the package since they participated in the randomized trial and contain largely complete data. The at endpoint can be 0 = censored, 1 = transplant or 2 = dead.
library(survival)
data("pbc")
sel <- c("id", "time", "status", "trt", "age", "sex", "bili", "protime")
## baseline data
head(pbc[,sel])## id time status trt age sex bili protime
## 1 1 400 2 1 58.76523 f 14.5 12.2
## 2 2 4500 0 1 56.44627 f 1.1 10.6
## 3 3 1012 2 1 70.07255 m 1.4 12.0
## 4 4 1925 2 1 54.74059 f 1.8 10.3
## 5 5 1504 1 2 38.10541 f 3.4 10.9
## 6 6 2503 2 2 66.25873 f 0.8 11.0Variables bili and protime are time-dependent measurements. The baseline data shown above, in contrast, only include baseline measurements of these covariates that have multiple measurements over time.
sel <- c("id", "futime", "status", "trt", "age", "sex", "day", "bili",
"protime")
head(pbcseq[,sel], 14)## id futime status trt age sex day bili protime
## 1 1 400 2 1 58.76523 f 0 14.5 12.2
## 2 1 400 2 1 58.76523 f 192 21.3 11.2
## 3 2 5169 0 1 56.44627 f 0 1.1 10.6
## 4 2 5169 0 1 56.44627 f 182 0.8 11.0
## 5 2 5169 0 1 56.44627 f 365 1.0 11.6
## 6 2 5169 0 1 56.44627 f 768 1.9 10.6
## 7 2 5169 0 1 56.44627 f 1790 2.6 11.3
## 8 2 5169 0 1 56.44627 f 2151 3.6 11.5
## 9 2 5169 0 1 56.44627 f 2515 4.2 11.5
## 10 2 5169 0 1 56.44627 f 2882 3.6 11.5
## 11 2 5169 0 1 56.44627 f 3226 4.6 11.5
## 12 3 1012 2 1 70.07255 m 0 1.4 12.0
## 13 3 1012 2 1 70.07255 m 176 1.1 12.0
## 14 3 1012 2 1 70.07255 m 364 1.5 12.0In a baseline analysis (including only baseline measurements), first without and then with adjustment for covariates, we obtain the following results:
## baseline fit for first 312 patients
fit0 <- coxph(Surv(time, status == 2) ~ trt,
data = pbc, subset = (id <= 312))
fit1 <- coxph(Surv(time, status == 2) ~ trt + log(bili) + log(protime),
data = pbc, subset = (id <= 312))
knitr::kable(tableRegression(fit0, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| trt | 0.94 | from 0.66 to 1.34 | 0.75 |
| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| trt | 0.91 | from 0.64 to 1.29 | 0.59 |
| log(bili) | 2.64 | from 2.20 to 3.17 | < 0.0001 |
| log(protime) | 78.76 | from 12.60 to 492.44 | < 0.0001 |
A time-dependent analysis includes all measurements. First, we need to prepare the data accordingly. The data preparation below is adopted from vignette in library .
temp <- subset(pbc, id <= 312, select = c(id:sex, stage)) # baseline
pbc2 <- tmerge(temp, temp, id = id, death = event(time, status)) #set range
pbc2 <- tmerge(pbc2, pbcseq, id = id, ascites = tdc(day, ascites),
bili = tdc(day, bili), albumin = tdc(day, albumin),
protime = tdc(day, protime), alk.phos = tdc(day, alk.phos))
head(pbc2)## id time status trt age sex stage tstart tstop death
## 1 1 400 2 1 58.76523 f 4 0 192 0
## 2 1 400 2 1 58.76523 f 4 192 400 2
## 3 2 4500 0 1 56.44627 f 3 0 182 0
## 4 2 4500 0 1 56.44627 f 3 182 365 0
## 5 2 4500 0 1 56.44627 f 3 365 768 0
## 6 2 4500 0 1 56.44627 f 3 768 1790 0
## ascites bili albumin protime alk.phos
## 1 1 14.5 2.60 12.2 1718
## 2 1 21.3 2.94 11.2 1612
## 3 0 1.1 4.14 10.6 7395
## 4 0 0.8 3.60 11.0 2107
## 5 0 1.0 3.55 11.6 1711
## 6 0 1.9 3.92 10.6 1365## time-dependent fit for first 312 patients
fit2 <- coxph(Surv(tstart, tstop, death == 2) ~ trt + log(bili) +
log(protime), data = pbc2)
knitr::kable(tableRegression(fit2, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| trt | 0.94 | from 0.65 to 1.34 | 0.72 |
| log(bili) | 3.46 | from 2.86 to 4.19 | < 0.0001 |
| log(protime) | 53.39 | from 22.78 to 125.14 | < 0.0001 |
Example 9.10 Survival by tumor response revisited, see Section 9.2. This example is about survival in patients with advanced lung cancer shown in Figure 9.13.
KMfit <- survfit(Surv(time/365.25, status) ~ 1, data = lung)
ggsurvplot(KMfit, xlab = "Years post diagnosis", ylab = "Survival",
conf.int = TRUE)
Figure 9.13: Kaplan-Meier estimate of the survival curve of patients with advanced lung cancer.
We now randomly generate tumor “responders” with a probability of \(0.05\) to become a responder at every month within the first year.
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## [1,] 0 0 0 0 0 0 0 0 0 0 0 NA NA
## [2,] 0 0 0 0 0 1 1 1 1 1 1 1 1
## [3,] 0 0 0 0 0 0 1 1 1 1 1 1 1## responder
## FALSE TRUE
## 145 83For illustration, we first show a wrong analysis comparing responder groups as in Figure 9.14:
badfit <- survfit(Surv(time/365.25, status) ~ responder, data = lung)
ggsurvplot(badfit, xlab = "Years post diagnosis", ylab = "Survival",
conf.int = TRUE, pval = TRUE)
Figure 9.14: Kaplan-Meier estimate of the survival curves of patients with advanced lung cancer, by responder status (wrong analysis).
badcoxfit <- coxph(Surv(time/365.25, status) ~ responder, data = lung)
knitr::kable(tableRegression(badcoxfit, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| responderTRUE | 0.54 | from 0.39 to 0.74 | 0.0002 |
There appears to be an effect of tumor response on survival, but it should not have an effect as responders were generated randomly. The analysis is flawed as it violates ``Breslow’s first rule’’ and is subject to time-dependent bias. More precisely, the bias is that the probability to be a responder increases with time.
At the beginning, we do not know who becomes a responder. We should therefore treat it as a time-dependent covariate. For the correct analysis with time-dependent covariate, we need to split the data by responder status. This augmented data set looks as follows:
## id tstart tstop response status time endpt
## 1 1 0.0 306.0 0 2 306 2
## 2 2 0.0 152.5 0 2 455 0
## 3 2 152.5 455.0 1 2 455 2
## 4 3 0.0 183.0 0 1 1010 0
## 5 3 183.0 1010.0 1 1 1010 1
## 6 4 0.0 91.5 0 2 210 0goodcoxfit <- coxph(Surv(tstart, tstop, endpt==2) ~ response, data=newdata)
knitr::kable(tableRegression(goodcoxfit, latex = FALSE))| Hazard Ratio | 95%-confidence interval | \(p\)-value | |
|---|---|---|---|
| response | 0.87 | from 0.61 to 1.23 | 0.42 |
In the correct analysis, there is no evidence for an effect of tumor response on survival, as it should be for randomly generated responders.
9.5 Additional references
Additional references are M. Bland (2015) (Ch. 16 Time to Event Data), J. N. S. Matthews (2006) (Sec. 7.5 Survival Analysis) and Collett (2003). Studies where the methods from this chapter are applied in practice are for example Blue et al. (2001), Colorectal Cancer Collaborative Group (2000), Zar et al. (2007).