library(RwR)
library(ggplot2)
library(tibble)
library(dplyr)
library(tidyr)
library(skimr)
library(broom)
library(lubridate)
library(survival)
library(splines)11 Survival analysis
Warning
This chapter is still a draft and minor changes are made without notice.
This chapter deals with regression models for time-to-event data with possible right censoring, which is common for survival data but also in other applications of statistics1 such as reliability theory and credit risk.
1 Time to failure of an electronic component, say, or time to default of a debtor.
Some of the notation in this chapter will deviate a bit from the notation used for other regression models, which is deliberate and intended to follow some of the standard conventions in the survival literature. The outcome variable will, for example, be denoted \(T\) instead of \(Y,\) which is also a natural notation for a time-to-event observation. The regression models are formed in terms of a linear predictor \(\eta = X^T \beta\) based on the vector of predictors \(X \in \mathbb{R}^p\) and a parameter \(\beta \in \mathbb{R}^p,\) and the effect of the predictors on the survival distribution is then always through \(\eta.\)
11.1 Outline and prerequisites
Section 11.2 covers the basic definitions specific to survival analysis and motivates why survival (regression) analysis is special. It also covers the basic nonparametric estimators of survival distributions. Section 11.3 covers parametric survival models and likelihood based methods for estimation and statistical inference. The accelerated failure time models and the proportional hazards models are introduced, which show two different ways that the linear predictor can affect the survival distribution. Section 11.4 primarily treats Cox’s proportional hazards model, which is the most widely used survival regression model. As the name suggests, it is a proportional hazards model, but the methodology is semiparametric. In this case it means that a baseline hazard function is left unconstrained. Section 11.5 introduces survival residuals, which differ a bit from residuals considered for linear and generalized linear models. These differences reflect that survival models do not target conditional expectations, but otherwise the purposes of the residuals for model diagnostics are the same.
The chapter primarily relies on the following R packages.
11.2 Introduction to survival analysis
A typical application of survival analysis is to the study of survival of patients after an initial diagnosis. In such a study, patients are enrolled whenever they are diagnosed with a given (serious, life threatening) disease. We refer to the time of enrollment2 as the baseline, and we note that the baseline for different patients can be different calendar times, and that the survival time is the time from baseline to death. The baseline is thus time zero on the time axis of the survival time.
2 For some applications, all individuals or units can have the same baseline, e.g., in an experiment on the longevity of electronic components, but in many applications the baseline calendar times differ among units. In credit risk models of time to default, say, the baseline could be the date of the loan application.
The patients are followed for some period of time, and if they die during this followup period, their survival time is registered. At a planned calendar time the study ends and the status of each patient is registered, that is, we register if they died, in which case we know their survival time, or if they are still alive. All patients alive at the end of the study have right censored survival times. This form of censoring is sometimes called administrative censoring. Some patients might also have dropped out of the study earlier, because they moved or did not want to participate anymore, say. At the time they drop out, we say that they are lost-to-followup, and this also results in right censoring of their survival time.
In addition to the, possibly right censored, survival times of the patients, we may collect additional data on the patients. We can then use some or all of such additional variables to construct a regression model of the survival distributions. But it is not clear that we want to model expected survival times, and the conditional expectation is often not3 the correct target of a survival regression model.
3 One reason is technical. Expectations rely on the tails of distributions, and with right censoring we do not get that information.
Regression models for survival data thus differ from the linear and generalized linear regression models in at least the following two ways:
- There is almost always a censoring mechanism, and certain aspects of the data are consequently missing. We need to deal with this in the modeling.
- Survival distributions are skewed distributions on the positive half line. It is often the shape of the distribution rather than the location of the distribution that is of interest.
The overall objective of a survival regression model is the same as for any regression model: how are the different predictors associated with the survival time. We might then have different more specific questions we want to answer. For example, to determine if the survival distribution depends on a choice of treatment (testing). Or we might want to quantify how strongly a predictor is associated with the survival outcome.
We might, of course, also be interested in prediction, which we often call survival prognosis. It is rarely of interest to make a point prediction of a survival time for an individual, but we might report a prediction interval (the interquartile range, say), or we might report probabilistic predictions in terms of an entire (conditional) survival distribution.
Example 11.1 We will throughout this chapter use a dataset on prostate cancer survival, which is included in the RwR package as the dataset prostate. The data is from a randomized clinical trial on prostate cancer patients with four different treatments (one placebo and three different dosages of estrogen). There are 502 patients in the trial, and we will in this chapter focus on four of the variables.
rx: Treatment encoded as a factor.status: Either the cause of death or indication that the patient was alive at the end of study.dtime: The censored survival time since baseline in months.age: The baseline age of the patient in years.
There are several other variables included in the dataset, and among these are sdate: the date that the patient entered the study. We will not use sdate for modeling, but it is used in Figure 11.1 to visualize the data. For later purposes, the data is converted to a tibble.
Code
data(prostate, package = "RwR") # version >= 0.3.5
prostate <- prostate |> as_tibble()# A tibble: 502 × 5
rx status dtime age sdate
<fct> <fct> <dbl> <dbl> <date>
1 0.2 mg estrogen alive 72 75 1977-08-10
2 0.2 mg estrogen dead - other ca 1 54 1977-09-21
3 5.0 mg estrogen dead - cerebrovascular 40 69 1978-01-12
4 0.2 mg estrogen dead - cerebrovascular 20 75 1978-03-19
5 placebo alive 65 67 1978-03-22
6 0.2 mg estrogen dead - prostatic ca 24 71 1978-06-14
7 placebo dead - heart or vascular 46 75 1978-06-27
8 placebo alive 62 73 1978-06-28
9 1.0 mg estrogen alive 61 60 1978-08-01
10 1.0 mg estrogen alive 60 78 1978-08-07
# ℹ 492 more rows
Figure 11.1 illustrates the data for the first 100 patients that enter the study. The figure shows their entry date (baseline), their survival time from that date and the status indicator when follow-up ends, indicating if they are dead or if the survival time is censored.
Since the treatments were randomly assigned, differences in the survival distributions between the four treatment groups can be ascribed to the treatment alone. That is, the treatments can be interpreted as the causes of any such differences. We will first of all be interested in testing if there are effects of treatment as well as estimating the magnitudes of such effects.
There are other variables included in the dataset, such as age, and we could use these variables for two different purposes. We could use the variables for prognostic purposes, that is, for fitting regression models that give more accurate survival prognoses for individual patients. Another purpose could be to investigate heterogeneity of treatment effects, that is, to investigate if the treatment works differently on different subgroups of the patients.
11.2.1 Survival distributions and censoring
A survival time \(T^*\) is simply a non-negative random variable. We will denote its distribution function by \(F,\) and if the distribution has density \(f\) w.r.t. Lebesgue measure, we have \[ F(t) = \mathbf{P}(T^* \leq t) = \int_{0}^{t} f(s) \mathrm{d}s. \] The corresponding survival function of \(T^*\) is \[ S(t) \coloneqq 1 - F(t) = \mathbf{P}(T^* > t) = \int_t^{\infty} f(s) \mathrm{d} s. \]
In terms of \(S\) we define the cumulative hazard function as \[ \Lambda(t) \coloneqq - \log(S(t)), \tag{11.1}\] so that \[ S(t) = \exp\left(-\Lambda(t)\right). \tag{11.2}\] If \(S\) is continuously differentiable, \(-S'(t) = f(t)\) is the density of the survival distribution, and the hazard function is defined as \[ \lambda(t) \coloneqq \Lambda'(t) = - (\log S(t))' = \frac{f(t)}{S(t)}. \tag{11.3}\]
There are one-to-one relations between the distribution function, the survival function and the cumulative hazard function.
We can observe that \[ \begin{align*} \lambda(t) & = \lim_{\varepsilon \to 0+} \frac{1}{\varepsilon} \frac{F(t+\varepsilon) - F(t)}{S(t)} \\ & = \lim_{\varepsilon \to 0+} \frac{1}{\varepsilon} \mathbf{P}(T^* \in (t,t+\varepsilon] \, | \, T^* > t), \end{align*} \] and \(\lambda(t)\) thus has the interpretation as the instantaneous rate of death at time \(t\) given survival beyond time \(t.\) We also note that \[ \Lambda(t) = \int_0^t \lambda(s) \mathrm{d} s \] which explains why \(\Lambda\) is called the cumulative hazard function.
We will not observe \(T^*\) but only \[ T = \min\{T^*, C\} \quad \text{ and } \quad \delta = 1(T^* \leq C) \] for a censoring variable \(C.\) We call \(T\) the right censored survival time and \(\delta\) the status indicator. The status indicator is \(1\) if \(T\) is uncensored and \(0\) otherwise. To distinguish the distribution of \(C\) clearly from the distribution of the survival time, we will generally denote its survival function by \(H\) and its density by \(h\) (if it exists). That is, \[ H(t) = \mathbf{P}(C > t) = \int_t^\infty h(s) \mathrm{d} s. \]
11.2.2 Model assumptions
We will generally consider \(n\) observations and let \(T^*_1, \ldots, T^*_n\) denote the survival times and \(C_1,\ldots,C_n\) the censoring times. We then observe the right censored survival times \[ T_i = \min\{T_i^*, C_i\} \] and the status indicators \[ \delta_i = 1(T_i^* \leq C_i) \] for \(i = 1, \ldots, n.\)
The vector of predictor variables is as elsewhere denoted \(X_i \in \mathbb{R}^p,\) and we are generally interested in models of the conditional distributions4 \[ T_i^* \mid X_i = x. \]
4 Note that the target is the conditional distribution of the unobserved and uncensored \(T^*_i,\) not the observed and censored \(T_i.\)
An issue of particular importance for survival modeling is to understand when a predictor variable is observed. In some studies, some potential predictor variables are monitored throughout the study. Such variables as BMI and blood pressure might be recorded on a monthly basis, say. It is possible to develop a general time-to-event theory where predictors also change with time, but it is both technically beyond this chapter, and it is also conceptually challenging to understand and interpret such models correctly.
We will throughout assume that all predictors are observed5 at or before baseline, and we refer to such predictors as baseline predictors.
5 Or, as a minimum, observable at time zero.
The modeling assumptions across the different individuals are now straightforward modifications of A5 or A6 to the survival setup.
SA5: The conditional distribution of \((T^*_i, C_i)\) given \(\mathbf{X}\) depends only on \(X_i\) and \((T^*_1, C_1), \ldots, (T^*_n, C_n)\) are conditionally independent given \(\mathbf{X}.\)
SA6: The triples \((T^*_1, C_1, X_1), \ldots, (T^*_n, C_n, X_n)\) are i.i.d.
Since \((T_i, \delta_i)\) is a function of \((T^*_i, C_i),\) then either assumption implies the same distributional statement about the \((T_i, \delta_i)\)-s as about the \((T^*_i, C_i)\)-s.
An additional delicate survival modeling question is what to assume about the joint distribution of \(T^*_i\) and \(C_i\) (possibly, conditionally on \(X_i\)). Throughout, we will make the following assumption.
Independent censoring: \(T^*_i\) and \(C_i\) are conditionally independent given \(X_i.\)
Even though the assumption above is a conditional independence statement, it is known simply as “independent censoring”. In the literature, “independent censoring” may refer to a slightly more permissive assumption (Definition III.2.1 in Andersen et al. (1993)), stating that information about censoring status at any given timepoint is not allowed to affect the residual survival distribution beyond that timepoint. For this chapter we will make the stronger, but simpler, independent censoring assumption as stated above.
Andersen, Per Kragh, Ørnulf Borgan, Richard D. Gill, and Niels Keiding. 1993. Statistical Models Based on Counting Processes. Springer Series in Statistics. Springer-Verlag.
11.2.3 The Kaplan-Meier and Nelson-Aalen estimators
If \(T_1^*, \ldots, T_n^*\) are i.i.d. and if they were observed, we could define the empirical survival function \[ \hat{S}^*_n(t) = \frac{1}{n} \sum_{i=1}^n 1(T_i^* > t) \] as a nonparametric estimator of \(S.\) Clearly, we cannot compute \(\hat{S}^*_n\) if we only have censored observations. The purpose of this section is to introduce two nonparametric estimators of a survival distribution when the survival times are i.i.d. and independently censored6.
6 In this section there are no predictors, and we effectively assume that \((T_1^*, C_1), \ldots, (T_n^*, C_n)\) are i.i.d. and that \(T_i^*\) and \(C_i\) are independent.
To handle the censoring we introduce the individuals at risk at time \(t\): \[ Y(t) = \sum_{i=1}^n 1(t \leq T_i). \] It can be computed based on the censored survival times, and it counts how many individuals have survived to time \(t\) and are uncensored at that time.
Definition 11.1 below gives the general definition of the Kaplan-Meier estimator of the survival function, but to motivate this definition, we first suppose that there are no ties among the observed survival times, that is, where we have strict inequalities among the ordered survival times \[ T_{(1)} < T_{(2)} < \ldots < T_{(n-1)} < T_{(n)}. \] Here \(T_{(k)}\) denotes the survival time with rank \(k,\) and we will denote the corresponding status indicator by \(\delta_{(k)}.\)
Without censoring, the Kaplan-Meier estimator equals \(\hat{S}^*,\) cf. Exercise 11.2
When there are no ties, and \(T_{(k)} \leq t < T_{(k+1)},\) we can factorize the survival function as \[ \begin{align*} S(t) & = \mathbf{P}(T^* > t) \\ & = \mathbf{P}(T^* > t \mid T^* \geq T_{(k)}) \times \mathbf{P}(T^* \geq T_{(k)} \mid T^* \geq T_{(k-1)}) \\ & \hspace{5mm} \ldots \times \mathbf{P}(T^* \geq T_{(2)} \mid T^* \geq T_{(1)}) \times \mathbf{P}(T^* \geq T_{(1)}) \\ \end{align*} \] where \(\mathbf{P}(T^* \geq T_{(i+1)} \mid T^* \geq T_{(i)})\) is the conditional probability of surviving from \(T_{(i)}\) to \(T_{(i+1)}\) given survival beyond time \(T_{(i)}.\)
Since \(Y(T_{(i)})\) is the number of individuals at risk7 at time \(T_{(i)}\) and \(Y(T_{(i)}) - \delta_{(i)}\) is the number of individuals that survive the interval \([T_{(i)}, T_{(i+1)}),\) we can regard \[ \frac{Y(T_{(i)}) - \delta_{(i)}}{Y(T_{(i)})} = 1-\frac{\delta_{(i)}}{Y(T_{(i)})} \] as an estimator of the conditional probability \(\mathbf{P}(T^* \geq T_{(i+1)} \mid T^* \geq T_{(i)}).\) This gives the Kaplan-Meier estimator for \(t \in [T_{(k)}, T_{(k+1)})\) \[ \hat{S}(t) = \left(1-\frac{\delta_{(k)}}{Y(T_{(k)})}\right) \cdots \left(1-\frac{\delta_{(2)}}{Y(T_{(2)})}\right) \left(1-\frac{\delta_{(1)}}{Y(T_{(1)})}\right). \]
7 Alive and uncensored.
We can write the Kaplan-Meier estimator in a slightly more abstract notation as \[ \hat{S}(t) = \prod_{i: T_i \leq t} \left(1-\frac{\delta_i}{Y(T_i)}\right). \tag{11.4}\]
To define the Kaplan-Meier estimator that allows for multiple deaths at the same time (ties), we introduce the counting process of deaths (non-censored events) \[N(t) = \sum_{i=1}^n 1(T_i \leq t, \delta_i=1),\] and the jumps for the counting process given as \[\Delta N(t) = N(t) - N(t-) = N(t) - \lim_{\varepsilon \to 0+} N(t-\varepsilon).\]
Definition 11.1 With the above notation, the Kaplan-Meier estimator is defined as \[ \hat{S}(t) = \prod_{s \leq t} \left(1-\frac{\Delta N(s)}{Y(s)}\right). \tag{11.5}\]
In the infinite product above, all factors where \(\Delta N(s) = 0\) are equal to \(1.\) There is thus only a finite number of factors different from \(1.\) When \(\Delta N(s) \in \{0, 1\}\) we see that Equation 11.5 coincides with Equation 11.4. The definition allows for \(\Delta N(s) > 1\) to correctly account for multiple deaths at the same time and thus ties among the observed survival times.
The function
survfit() from the R package survival can be used to compute the Kaplan-Meier estimator.Since the cumulative hazard function is \(\Lambda = - \log(S),\) we can estimate it nonparametrically as \[ - \log(\hat{S}(t)) = - \sum_{s \leq t} \log \left(1-\frac{\Delta N(s)}{Y(s)}\right). \tag{11.6}\] But there is also a different, and direct, nonparametric estimator of the cumulative hazard function. The logic behind this estimator is that for a differentiable survival function \(S\) and \(t > s\) (but \(t - s\) small),
\[ \begin{align*} \mathbf{P}(T^* \in [s, t) \mid T^* \geq s) & = \frac{S(s) - S(t)}{S(s)} \\ & \approx - \frac{S'(s)}{S(s)}(t - s) = \lambda(s)(t - s). \end{align*} \]
By a Riemann sum approximation we see that \[ \Lambda(t) \approx \sum_{i: T_{(i)} \leq t} \mathbf{P}(T^* \in [T_{(i)}, T_{(i + 1)}) \mid T^* \geq T_{(i)}), \] and if there are no ties, the conditional probabilities in the sum above can be estimated8 as \(\delta_{(i)} / Y(T_{(i)}).\) This leads to the Nelson-Aalen estimator \[ \hat{\Lambda}(t) = \sum_{i: T_{(i)} \leq t} \frac{\delta_{(i)}}{Y(T_{(i)})}. \tag{11.7}\]
8 Either \(\delta_{(i)} = 1\) and one individual out of \(Y(T_{(i)})\) dies in the interval (at time \(T_{(i)}\)), or \(\delta_{(i)} = 0\) and no individuals die in the interval.
The general definition below allows for multiple deaths at the same time.
Definition 11.2 The Nelson-Aalen estimator of the cumulative hazard function is
\[
\hat{\Lambda}(t) = \sum_{s \leq t} \frac{\Delta N(s)}{Y(s)}.
\tag{11.8}\]
Since the survival function is given by \(S(t) = \exp(-\Lambda(t)),\) we can estimate the survival function by plugging the Nelson-Aalen estimator of the cumulative hazard, which gives \[ \exp\left(-\hat{\Lambda}(t)\right) = \prod_{s \leq t} \exp\left(- \frac{\Delta N(s)}{Y(s)}\right). \tag{11.9}\]
We note that \(\exp(-\hat{\Lambda})\) is not the Kaplan-Meier estimator, nor is \(- \log(\hat{S})\) the Nelson-Aalen estimator.
The function
survfit() with argument stype = 2 estimates the survival function in terms of the Nelson-Aalen estimator via Equation 11.9.The mismatch between the Kaplan-Meier and the Nelson-Aalen estimators is a small nuisance, but we may note that \[ \exp\left(- \frac{\Delta N(s)}{Y(s)}\right) \approx \left(1 - \frac{\Delta N(s)}{Y(s)}\right) \] for large \(Y(s)\) by using the Taylor expansion \(\exp(x) = 1 + x + O(x^2).\) Thus differences between the factors are small unless \(s\) is so large that the number of individuals at risk, \(Y(s),\) is small.
One can argue that the way ties are handled above presumes that survival times can be exactly equal, which is at odds with any continuous survival distribution. If we regard ties as a consequence of lack of precision in our registration or due to rounding of the survival times, it makes sense to replace \(\frac{\Delta N(s)}{Y(s)}\) by \[ \frac{1}{Y(s)} + \frac{1}{Y(s) - 1} + \ldots + \frac{1}{Y(s) - \Delta N(s) + 1}. \] This is known as the Fleming-Harrington correction for ties. It often does not make much of a practical difference though.
The function
survfit() with argument ctype = 2 will use the Fleming-Harrington correction for ties.Example 11.2 We continue Example 11.1 by first considering a summary of the marginal distributions of the four variables.
Code
select(prostate, c(rx, status, dtime, age)) |> skim()Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| rx | 0 | 1 | FALSE | 4 | pla: 127, 1.0: 126, 5.0: 125, 0.2: 124 |
| status | 0 | 1 | FALSE | 10 | ali: 148, dea: 130, dea: 96, dea: 31 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| dtime | 0 | 1 | 36.13 | 23.32 | 0 | 14.25 | 34 | 57.75 | 76 | ▇▆▅▆▆ |
| age | 1 | 1 | 71.46 | 7.08 | 48 | 70.00 | 73 | 76.00 | 89 | ▁▂▅▇▁ |
prostate dataset.
Table 11.1 shows 148 patients (out of 502) were alive at the end of the study. Thus 148 of the survival times are right censored. Table 11.2 shows a more detailed tabulation of the values of the status variable, which specifies one out of nine different causes of death—in addition to specifying if the patient is still alive. We may also note that the survival time (dtime) is registered in months since baseline, which is a very coarse scale. Consequently there are many ties in the data.
| status | Freq |
|---|---|
| alive | 148 |
| dead - cerebrovascular | 31 |
| dead - heart or vascular | 96 |
| dead - other ca | 25 |
| dead - other specific non-ca | 28 |
| dead - prostatic ca | 130 |
| dead - pulmonary embolus | 14 |
| dead - respiratory disease | 16 |
| dead - unknown cause | 7 |
| dead - unspecified non-ca | 7 |
We will fit a survival distribution, using either of the nonparametric estimators above, where we ignore the information about different causes of death. Thus the status variable will only be used to determine the censoring status.
We fit survival functions for all patients using the survfit() function from the survival package. See ?survfit.formula for details about how to specify a model using the formula interface. Note that the left hand side of the formula, that specifies the outcome, uses the special Surv() function syntax to include information about the censoring status. We fit the survival functions below using both the Kaplan-Meier and the Nelson-Aalen estimators, and we also fit the survival function using the Nelson-Aalen estimator with the Fleming-Harrington correction for ties.
Code
prostate_KM <- survfit(
Surv(dtime, status != "alive") ~ 1,
data = prostate,
)
prostate_NA <- survfit(
Surv(dtime, status != "alive") ~ 1,
data = prostate,
stype = 2
)
prostate_NA_FH <- survfit(
Surv(dtime, status != "alive") ~ 1,
data = prostate,
stype = 2,
ctype = 2
)Figure 11.2 shows the fitted survival curves—including a plot that zooms in on a part of the curves. The differences are small and barely visible, and they are most pronounced in the tail of the distribution.
Code
p1 <- summary(prostate_KM, data.frame = TRUE) |>
ggplot(aes(time, surv)) +
geom_step(data = summary(prostate_NA, data.frame = TRUE), color = "blue") +
geom_step(data = summary(prostate_NA_FH, data.frame = TRUE), color = "red") +
geom_step() +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability")
p2 <- p1 + coord_cartesian(xlim = c(50, 75), ylim = c(0.20, 0.40))
gridExtra::grid.arrange(p1, p2, ncol = 2)
We then fit a separate survival function for each treatment group. We can think about this as a simple regression model with a single categorical predictor variable (the treatment variable rx), and the formula specification below is similar to other regression model specifications. The main difference is that using survfit() we fit a survival function for each treatment group instead of a mean value, say.
Code
prostate_KM_rx <- survfit(
Surv(dtime, status != "alive") ~ rx,
data = prostate
)Code
summary(prostate_KM_rx, data.frame = TRUE) |>
ggplot(aes(time, surv, color = strata)) +
geom_step() +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability") +
theme(legend.position = "top")
11.3 Parametric survival models
In this section we focus on parametric models of the survival distribution, and we will first derive a generic expression for the likelihood function. We then turn to regression models and introduce the accelerated failure time models and the proportional hazards models as two different ways for the predictors to affect the outcome distribution.
11.3.1 The likelihood
Recall that in terms of \(T^*\) and \(C\) we have \(T = \min\{T^*, C\}\) and \(\delta = 1(T^* \leq C).\) Our first result gives the density of the joint distribution of \((T, \delta)\) in terms of the marginal distributions of \(T^*\) and \(C\) when they are independent.
Theorem 11.1 If \(T^*\) and \(C\) are independent the joint distribution of \((T, \delta)\) has density \[ g(t, \delta) = f(t)^\delta S(t)^{1-\delta} h(t)^{1-\delta} H(t)^\delta \] w.r.t. the product measure \(m \otimes \tau\) (Lebesgue measure times the counting measure).
Proof. First note that \[ \begin{align*} \mathbf{P}(T \leq t, \delta = 1) & = \mathbf{P}(T^* \leq t, \delta = 1) \\ & = \mathbf{P}(T^* \leq t,C \geq T^*) \\ & = \int_0^t f(s) \int_s^{\infty} h(u) \mathrm{d} u \mathrm{d} s \\ & = \int_0^t f(s) H(s) \mathrm{d} s. \end{align*} \] Likewise, \[ \mathbf{P}(T \leq t, \delta = 0) = \int_0^t h(s) S(s) \mathrm{d} s, \] and we conclude that the density is \[ g(t, \delta) = \left\{ \begin{array}{ll} f(t) H(t) & \text{if } \delta = 1 \\ h(t) S(t) & \text{if } \delta = 0. \end{array} \right. \]
Based on the density result above we can introduce the full likelihood for right censored observations. With \((T_1, \delta_1), \ldots, (T_n, \delta_n)\) independent, \(F_i\) the distribution function for \(T_i^*\) and \(H_i\) the survival function for the censoring variable \(C_i\), the full likelihood is \[ L(F, H) = \prod_{i=1}^n f_i(T_i)^{\delta_i} S_i(T_i)^{1-{\delta_i}} h_i(T_i)^{1-{\delta_i}} H_i(T_i)^{\delta_i}, \] where \(F\) and \(H\) are shorthand for \(F_1, \ldots, F_n\) and \(H_1, \ldots, H_n\).
If the survival distribution can vary independently of the censoring distribution, the full likelihood can be seen as a function of \(F\) only, \[ F \mapsto \prod_{i=1}^n f_i(T_i)^{\delta_i} S_i(T_i)^{1-{\delta_i}} K_i = \left(\prod_{i=1}^n f_i(T_i)^{\delta_i} S_i(T_i)^{1-{\delta_i}}\right) \left(\prod_{i=1}^n K_i\right), \] with \(K_i\) some constants depending on the observations but not on the survival distribution \(F.\) For estimation and inference purposes, we can drop the factor \(\prod_{i=1}^n K_i\) since it is only the first factor that depends on \(F\).
Definition 11.3 The survival likelihood is \[ L(F) \coloneqq \prod_{i=1}^n f_{i}(T_i)^{\delta_i} S_{i}(T_i)^{1-{\delta_i}}. \tag{11.10}\]
For a parametrized model of the distributions of \(T^*_1, \ldots, T^*_n\) we have that \(f_i = f_{i, \beta}\) for some parameter \(\beta.\) This gives the survival likelihood as a function of \(\beta\), \[ L(\beta) = \prod_{i=1}^n f_{i, \beta}(T_i)^{\delta_i} S_{i, \beta}(T_i)^{1-{\delta_i}}. \]
The survival likelihood is the likelihood we will use throughout this chapter. Its derivation made it explicit how this is a likelihood w.r.t. a dominating measure up to a factor that depends on the censoring mechanism only. We can ignore this factor when it does not depend on the unknown parameter \(\beta.\) The argument for that is not independent censoring9 as such but that the \(H_i\)-s do not depend upon \(\beta.\) Technically, the censoring distribution and the parameter \(\beta\) are assumed variation independent.
9 Though independent censoring plays a role here for our particular derivation and form of the full likelihood.
It is possible to take a slightly different point of view and condition on the censoring variables instead. Then the dominating measure for the \(i\)-th censored observation, \(T_i,\) becomes Lebesgue measure plus a unit point mass in \(C_i.\) One arrives at the same likelihood as above though.
We can also express the likelihood in terms of the hazard and cumulative hazard functions.
Proposition 11.1 If \(F_i\) is the distribution function of \(T_i^*\) with hazard function \(\lambda_{i}\) and cumulative hazard function \(\Lambda_{i},\) the survival likelihood equals \[ L(F) = L(\Lambda) \coloneqq \prod_{i=1}^n \lambda_{i}(T_i)^{\delta_i} e^{- \Lambda_{i}(T_i)} \tag{11.11}\]
Proof. Recall that \[ \lambda_{i}(t) = \frac{f_{i}(t)}{S_{i}(t)} \] and that \[ S_{i}(t) = e^{- \Lambda_{i}(t)}. \] Hence \[ \begin{align*} L(F) & = \prod_{i=1}^n f_{i}(T_i)^{\delta_i} S_{i}(T_i)^{1-{\delta_i}} \\ & = \prod_{i=1}^n \lambda_{i}(T_i)^{\delta_i} S_{i}(T_i)^{\delta_i} S_{i}(T_i)^{1-{\delta_i}} \\ & = \prod_{i=1}^n \lambda_{i}(T_i)^{\delta_i} S_{i}(T_i) \\ & = \prod_{i=1}^n \lambda_{i}(T_i)^{\delta_i} e^{- \Lambda_{i}(T_i)} = L(\Lambda). \end{align*} \]
If \(f_{i, \beta}\) is a parametrized density for the survival distribution with hazard function \(\lambda_{i, \beta}\) and cumulative hazard function \(\Lambda_{i, \beta},\) the likelihood as a function of \(\beta\) can then be expressed as \[ L(\beta) = \prod_{i=1}^n \lambda_{i, \beta}(T_i)^{\delta_i} e^{- \Lambda_{i, \beta}(T_i)}. \tag{11.12}\]
Example 11.3 The Weibull distribution has hazard function \[ \lambda(t) = \alpha \gamma t^{\gamma - 1} \] for \(\alpha, \gamma > 0\) two parameters. The special case \(\gamma = 1\) gives us the exponential distribution with constant hazard function. The Weibull distribution has increasing hazard functions over time when \(\gamma > 1\) and decreasing hazard function when \(\gamma < 1.\)
The cumulative hazard function is \[ \Lambda(t) = \alpha \int_0^t \gamma s^{\gamma - 1} \mathrm{d} s = \alpha t^\gamma. \]
For right censored i.i.d. observations from the Weibull model, the likelihood becomes \[ L(\alpha, \gamma) = \prod_{i=1}^n (\alpha \gamma)^{\delta_i} T_i^{\delta_i(\gamma - 1)} e^{- \alpha T_i^\gamma}, \] and the log-likelihood is \[ \begin{align*} \ell(\alpha, \gamma) & = \log L(\alpha, \gamma) \\ & = \sum_{i=1}^n \delta_i (\log(\alpha) + \log(\gamma) + (\gamma - 1)\log(T_i)) - \alpha T_i^\gamma. \end{align*} \]
Example 11.4 (MLE for the exponential distribution) Continuing Example 11.3 and fixing \(\gamma = 1,\) corresponding to an exponential distribution with hazard rate \(\alpha,\) we find that the log-likelihood equals \[ \ell(\alpha) = \sum_{i=1}^n \delta_i \log(\alpha) - \alpha T_i \] up to additive terms that do not depend on \(\alpha.\) Provided that at least one observation is uncensored, this log-likelihood has a unique stationary point \[ \hat{\alpha} = \frac{\sum_{i=1}^n \delta_i}{\sum_{i=1}^n T_i} = \frac{n_u}{\sum_{i=1}^n T_i} \] where \(n_u = \sum_{i=1}^n \delta_i\) is the number of deaths (the number of uncensored events). The second derivative of \(\ell\) is \(-n_u / \alpha^2 < 0,\) so \(\ell\) is strictly concave and \(\hat{\alpha}\) is the MLE.
If we had treated the observations as if there were no censoring, we would get the MLE \[ \frac{n}{\sum_{i=1}^n T_i} > \hat{\alpha}, \] which overestimates the rate.
If we, on the other hand, discard the censored observations, we would get the MLE \[ \frac{n_u}{\sum_{i=1}^n \delta_i T_i} > \hat{\alpha}, \] which likewise overestimates the rate.
11.3.2 Regression models
We introduce two ways that predictor variables can affect the survival outcome distribution: the accelerated failure time models, exemplified by the log-logistic model, and the proportional hazards models, exemplified by the Weibull model. In the latter case, it turns out that the likelihood can be computed and optimized using methods for Poisson generalized linear models.
Accelerated failure time models
Example 11.5 (Log-logistic distribution) Suppose \(Z\) has a logistic distribution with distribution function \[ G(z) = \frac{e^{\lambda z}}{1 + e^{\lambda z}}. \] The distribution of \(T^* = e^Z\) is then the log-logistic distribution, whose distribution function is \[ F(t) = G(\log(t)) = \frac{t^{\lambda}}{1 + t^{\lambda}}, \] with corresponding density \[ f(t) = F'(t) = \frac{\lambda t^{\lambda-1}}{(1+t^{\lambda})^2} \] for \(t > 0.\)
Definition 11.4 Let \(G\) be a distribution function (on \(\mathbb{R}\)). The corresponding accelerated failure time (AFT) model has survival function given as \[ S_{\eta} (t) = 1 - G( (\log(t) - \eta)/\sigma) \] with \(\eta = X^T \beta\) the linear predictor and \(\sigma > 0\) a scale parameter.
The AFT model is a scale-location model of the log-survival time.
The interpretation of the AFT model is that a unit change of the \(j\)-th predictor, with the other predictors fixed, multiplies the failure time by a factor \(e^{\beta_j}.\) The predictor thus changes the time scale, and it accelerates (\(\beta_j < 0\)) or decelerates (\(\beta_j > 0\)) the time to failure.
Example 11.6 (Log-logistic AFT model) If \(G(z) = e^z/(1 + e^z)\) is the standard logistic distribution function and the scale parameter is \(\sigma = \lambda^{-1},\) the survival function for the corresponding AFT model is \[ \begin{aligned} S_{\eta}(t) & = \frac{1}{1 + e^{\lambda(\log(t) - \eta)}} \\ & = \frac{1}{1 + e^{-\lambda \eta} t^{\lambda}}. \end{aligned} \tag{11.13}\] where \(\eta = X^T \beta\) is the linear predictor. This is a log-logistic distribution as introduced in Example 11.5. The effect of \(\eta\) is10 a scale transformation of the log-logistic distribution with scale parameter \(e^{\eta}.\) The density is \[ \begin{aligned} f_{\eta}(t) & = \frac{\lambda e^{-\lambda \eta} t^{\lambda-1}}{(1+e^{-\lambda \eta} t^{\lambda})^2} \\ & = \frac{\lambda e^{\lambda (\log(t) - \eta)} t^{-1}}{(1+e^{\lambda(\log(t) - \eta)})^2}. \end{aligned} \tag{11.14}\]
10 If \(T\) has log-logistic distribution with density \(f_0,\) then \(e^{\eta} T\) is log-logistic distributed with density \(f_{\eta},\) cf. Exercise 11.1.
For maximum likelihood estimation we need to compute and optimize the likelihood in Equation 11.10. Efficient computation requires access to the density as well as the survival function. Both are explicitly available for the log-logistic distribution, which makes this distribution attractive for AFT models.
Example 11.7 We implement the negative log-likelihood (nll()) for the AFT model with a log-logistic survival distribution below. The implementation is based on Equation 11.10 using the expressions (11.13) and (11.14) for the survival function and density, respectively, of the log-logistic distribution.
The nll() function takes four arguments. Besides a vector of parameters (par), the censored survival times (time) and the status indicator (status), it takes a single additional argument x representing a categorical variable with \(k\) values. The parameter vector par is of the form \[
(\sigma, \beta_0, \ldots, \beta_{k-1}) = (\lambda^{-1}, \beta_0, \ldots, \beta_{k-1})
\] with \(\sigma\) the scale parameter, \(\beta_0\) the intercept and \(\beta_1, \ldots, \beta_{k-1}\) are contrasts. The survival function is also implemented below for later usage.
Code
# Argument `x` is an integer representing the k levels of a single
# categorical predictor variable. `par[2:(k+1)]` are the k
# beta-parameters corresponding to the intercept and contrasts.
# Its default value corresponds to no predictor.
nll <- function(par, time, status, x = rep(1, length(time))) {
lambda <- 1 / par[1]
eta <- par[2] + c(0, par[-c(1, 2)])[x] # beta_x
uncens <- status == 1
lt <- log(time)
logf <- (log(lambda) - lambda * eta[uncens] + (lambda - 1) * lt[uncens] -
2 * log(1 + exp(lambda * (lt[uncens] - eta[uncens])))) |> sum()
neg_logS <- log(1 + exp(lambda * (lt[!uncens] - eta[!uncens]))) |> sum()
- logf + neg_logS
}We use optim() with method L-BFGS-B to compute the maximum-likelihood estimate of a single survival distribution. Since some of the patients have survived less than one month, some survival times are registered as \(0,\) which is numerically problematic (due to \(\log(0)\)). For this reason we add \(0.5\) to all times, corresponding to shifting all survival and censoring times by half a month.
Code
prostate_aft <- optim(
c(2, 2),
nll,
method = "L-BFGS-B",
lower = c(0.01, -Inf),
time = prostate$dtime + 0.5, # to avoid log(0)
status = prostate$status != "alive"
)
prostate_aft[1:2]$par
[1] 0.8219375 3.5378103
$value
[1] 1758.922
We see from the optim() result that the parameter estimates are \(\hat{\sigma} =\) 0.82 and \(\hat{\beta}_0 =\) 3.54, and we can also see that \(-\ell(\hat{\sigma}, \hat{\beta}_0) =\) 1758.92. The corresponding fitted survival function is compared to the Kaplan-Meier estimator in Figure 11.4.
We can also use our implementation to fit the regression model where the survival distribution depends on treatment.
Code
prostate_aft_rx <- optim(
c(2, 2, 0, 0, 0),
nll,
method = "L-BFGS-B",
lower = c(0.1, -Inf, -Inf, -Inf, -Inf),
time = prostate$dtime + 0.5, # to avoid log(0)
status = prostate$status != "alive",
x = as.numeric(prostate$rx)
)
prostate_aft_rx[1:2]$par
[1] 0.81992315 3.48305068 -0.04002854 0.33958359 -0.05824063
$value
[1] 1755.851Code
# The log-logistic survival function in the scale-location parametrization.
time0 <- 0:80
fit <- tibble(
time = rep(time0, 4),
surv = c(
S(time0 + 0.5, prostate_aft_rx$par[1], prostate_aft_rx$par[2]),
S(time0 + 0.5, prostate_aft_rx$par[1], prostate_aft_rx$par[2] + prostate_aft_rx$par[3]),
S(time0 + 0.5, prostate_aft_rx$par[1], prostate_aft_rx$par[2] + prostate_aft_rx$par[4]),
S(time0 + 0.5, prostate_aft_rx$par[1], prostate_aft_rx$par[2] + prostate_aft_rx$par[5])
),
strata = rep(paste("rx=", levels(prostate$rx), sep = ""), each = 81)
)
summary(prostate_KM_rx, data.frame = TRUE) |>
ggplot(aes(time + 0.5, surv, color = strata)) +
geom_step(linetype = 2) +
geom_line(data = fit) +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability") +
theme(legend.position = "top")
We can fit AFT models, and a range of other parametric regression models, via the survreg() function from the survival package.
Code
survreg(
Surv(dtime + 0.5, status != "alive") ~ rx,
data = prostate,
dist = "loglogistic"
)Call:
survreg(formula = Surv(dtime + 0.5, status != "alive") ~ rx,
data = prostate, dist = "loglogistic")
Coefficients:
(Intercept) rx0.2 mg estrogen rx1.0 mg estrogen rx5.0 mg estrogen
3.48297715 -0.03983952 0.33947524 -0.05802741
Scale= 0.8199158
Loglik(model)= -1755.9 Loglik(intercept only)= -1758.9
Chisq= 6.14 on 3 degrees of freedom, p= 0.105
n= 502 Note that the parametrization used by survreg() is exactly the same as we implemented, and that the estimated parameter values and computed log-likelihood are the same as well (up to small numerical differences).
The survreg() function even gives a likelihood-ratio test statistic (Chisq) of the hypothesis that there are no differences between the treatments. The corresponding \(p\)-value is computed using an approximating \(\chi^2\)-distribution with \(3\) degrees of freedom. It is above \(0.1\) and we cannot reject that there is no treatment effect.
As for lm and glm objects there are methods for survreg objects implemented in the broom package for the glance() and tidy() functions and in the survival package for anova(). We will illustrate some of their usages in subsequent examples.
Proportional hazards models
A different model class is given by letting the linear predictor have a multiplicative effect on the hazard function. The resulting proportional hazards models are among the most widely used survival models.
Definition 11.5 With \(\lambda_0\) a baseline hazard function, the proportional hazards model is given by the hazard function \[ \lambda(t) = \lambda_0(t) e^{\eta} \] with \(\eta\) the linear predictor.
It follows that for the cumulative hazard function, \[ \Lambda(t) = \Lambda_0(t) e^{\eta}, \] and the proportionality holds for the cumulative hazard as well. The factor \(e^{\beta_j}\) is the hazard ratio between two individuals whose \(j\)-th predictors differ by one unit while all other predictors are equal.
The proportional hazards model is sometimes called the Cox model. Cox (1972) discovered that for the proportional hazards model, the parameters that enter into the linear predictor can be estimated without specifying the baseline hazard function. This will be covered in Section 11.4. Here we will instead show how the parameters can be estimated via maximum-likelihood when the baseline is assumed to be a Weibull distribution.
Cox, D. R. 1972. “Regression Models and Life-Tables.” Journal of the Royal Statistical Society. Series B (Methodological) 34 (2): pp. 187–220.
Example 11.8 The Weibull baseline hazard function and cumulative hazard function are
\[
\lambda_0(t) = \gamma t^{\gamma-1} \quad \text{and} \quad \Lambda_0(t) = t^{\gamma}.
\] The log-likelihood is \[
\begin{align*}
\ell & = \sum_{i=1}^n \delta_i \log(\gamma T_i^{\gamma-1}e^{\eta_i}) - T_i^{\gamma}e^{\eta_i} \\
& = \underbrace{\sum_{i=1}^n \delta_i \log(T_i^{\gamma}e^{\eta_i}) - T_i^{\gamma}e^{\eta_i}}_{\text{Poisson log-likelihood}} +
\sum_{i=1}^n \delta_i \log(\gamma T_i^{-1}).
\end{align*}
\] This is (surprisingly) up to a constant the log-likelihood for a Poisson model of the \(\delta_i\)-s with logarithmic link function and mean value \(T_i^{\gamma}e^{\eta_i}\) for fixed \(\gamma.\) Note that the \(\alpha\) parameter appearing in Example 11.3 is captured by an intercept in the linear predictor.
For a fixed value of the linear predictor we also find that \[ \partial_{\gamma} \ell = \sum_{i=1}^n (\delta_i- T_i^{\gamma} e^{\eta_i}) \log T_i + \frac{\delta_i}{\gamma}. \] Thus \(\gamma\) solves the equation \[ \gamma = \frac{n_u}{\sum_{i=1}^n (T_i^{\gamma} e^{\eta_i} -\delta_i)\log T_i} \tag{11.15}\] with \(n_u\) the number of uncensored observations.
A practical algorithm can alternate between: i) estimate the linear predictor as \(\hat{\eta}^{(k)}\) by optimizing the Poisson log-likelihood (for fixed \(\gamma^{(k-1)}\)); and ii) update \[
\gamma^{(k)} = \frac{n_u}{\sum_{i=1}^n (T_i^{\gamma^{(k-1)}} e^{\hat{\eta}_i^{(k)}} -\delta_i)\log T_i}.
\]
The example above shows how implementations of generalized linear models— and Poisson regression models, in particular— can be used to fit the Weibull proportional hazards model (for fixed \(\gamma\)) with the survival times entering as an offset term. This was of historical importance as it made computations practical early on using some standard algorithms available. Today, the connection between the survival likelihood for Weibull distributions and the Poisson likelihood is mostly of theoretical and historical interest.
To fit the model we will instead show that the proportional hazards model with a Weibull baseline hazard function is also an accelerated failure time model, which makes it possible to use survreg() to fit the model.
Example 11.9 If we take \(G\) in Definition 11.4 to be the extreme value distribution \[ G(z) = 1 - \exp(-e^{z}), \] then the corresponding AFT model has survival function \[ S_\eta(t) = \exp\left(-(e^{-\eta}t)^{\frac{1}{\sigma}}\right) = \exp\left(-e^{-\gamma \eta}t^\gamma\right) \] where \(\gamma = \frac{1}{\sigma}.\) That is, the cumulative hazard function is \(\Lambda_\eta(t) = e^{- \gamma \eta}t^\gamma,\) which is the proportional hazards model with a Weibull baseline, but reparametrized in terms of the scale parameter \(\sigma = \frac{1}{\gamma},\) and the linear predictor is multiplied by \(-\gamma.\)
We can use the survreg() function to fit the model.
Code
prostate_survreg <- survreg(
Surv(dtime + 0.5, status != "alive") ~ rx,
data = prostate,
dist = "weibull"
)
prostate_survreg |> tidy()# A tibble: 5 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.87 0.105 36.9 1.05e-297
2 rx0.2 mg estrogen -0.0314 0.148 -0.212 8.32e- 1
3 rx1.0 mg estrogen 0.392 0.161 2.43 1.49e- 2
4 rx5.0 mg estrogen -0.00534 0.149 -0.0358 9.71e- 1
5 Log(scale) 0.0222 0.0470 0.472 6.37e- 1We can compute the estimated survival functions for the Weibull model much in the same way as we did above for the log-logistic model, where we used the MLE of the parameters and the explicit formula for the survival function. However, this is a bit cumbersome and does not generalize well. We can instead use the predict() method for survreg objects to predict quantiles for the survival functions, which implicitly also gives the estimated survival functions.
Code
fit <- tibble(rx = factor(1:4, labels = levels(prostate$rx)))
p_seq <- seq(0.01, 0.90, 0.01)
pred <- predict(prostate_survreg, newdata = fit, type = "quantile", p = p_seq)
fit <- cbind(fit, as_tibble(pred, .name_repair = "minimal")) |>
pivot_longer(cols = !rx, values_to = "time") |>
mutate(
surv = 1 - rep(p_seq, 4),
strata = paste("rx=", rx, sep = "")
) |>
filter(time <= 80)
summary(prostate_KM_rx, data.frame = TRUE) |>
ggplot(aes(time + 0.5, surv, color = strata)) +
geom_step(linetype = 2) +
geom_line(data = fit, aes(time)) +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability") +
theme(legend.position = "top") +
xlim(0, 82)
Figure 11.6 compares the fitted Weibull distributions to the Kaplan-Meier estimators for the four treatment groups, and Table 11.4 reports the test statistic and \(p\)-value for testing if there are no differences between the four treatment groups. The \(p\)-value is \(0.02\) and a bit smaller than when using the log-logistic distribution.
Code
prostate_survreg |>
anova() |>
select(-4) |>
tidy() |>
filter(term == "rx") |>
knitr::kable(digits = 3)| term | df | deviance | df.residual | p.value |
|---|---|---|---|---|
| rx | 3 | 9.619 | 497 | 0.022 |
To corroborate the technical relation between the log-linear Poisson regression model and the proportional hazards model with a Weibull baseline, we refit the \(\beta\)-parameters using glm() (fixing the scale parameter at its estimated value).
Code
sigma <- prostate_survreg$scale
glm((status != "alive") ~ offset(log(dtime + 0.5) / sigma) + rx,
data = prostate,
family = "poisson"
) |> coefficients() * (- sigma) (Intercept) rx0.2 mg estrogen rx1.0 mg estrogen rx5.0 mg estrogen
3.870053562 -0.031447138 0.392042914 -0.005335812 glm() and identical to those in Table 11.3.
To conclude our analyses of the prostate survival times, we have shown how one treatment group (1.0 mg estrogen) has a slightly different survival distribution than the other three with a better survival prognosis. For the statistical tests of any treatment effect, based on two parametric models, only one is significant. The signal is therefore weak and we cannot draw a clear conclusion on whether there is a treatment effect.
There does not seem to be any detectable difference between the other three treatment groups. It is particularly noteworthy that the treatment group with the highest dose (5.0 mg estrogen) has a survival distribution similar to the placebo group. If there were an effect, we would often expect it to be monotone in dosage, and it would be peculiar if the effect is only there for the medium dose.
11.4 Semiparametric survival models
For the parametric proportional hazards models in Definition 11.5, the baseline hazard function \(\lambda_0\) is either known or parametrized itself. Since any specific parametrized baseline model might be misspecified, we would like to avoid modeling the baseline hazard if possible. This leads to the semiparametric proportional hazards model, where an individual’s hazard is proportional to a common baseline hazard \(\lambda_0\) with the individual relative risk as proportionality factor, but where we do not model11 the baseline \(\lambda_0.\)
11 This is similar to not making distributional assumptions such as A3 or GA3 for the linear and generalized linear models.
12 The \(i\)-th relative risk could be \(e^{\eta_i}\) with \(\eta_i = X_i^T \beta\) the linear predictor, but it does not need to be for the general arguments.
The general proportional hazards model assumes that the dependence on the predictors is via a multiplicative factor on the baseline, and the \(i\)-th individual thus has hazard function \[ \lambda_i(t) = \psi_i \lambda_0(t) \] for \(\psi_i > 0\) the relative risk12 of the \(i\)-th individual. The cumulative hazard function is consequently \(\Lambda_i(t) = \psi_i \Lambda_0(t).\)
In the following sections we discuss parametrization and estimation of the relative risks via maximization of Cox’s partial likelihood. This is an estimation criterion based on the full likelihood but where the baseline hazard has been profiled out.
11.4.1 Cox’s partial likelihood
To ease notation we leave out limits on generic indices \(i\) and \(j\) in products and sums. They will always run through observations from \(1\) to \(n.\)
Definition 11.6 In terms of the observations \(T_1, \ldots, T_n\) and relative risks \(\psi_1, \ldots, \psi_n > 0\) we define the cumulative relative risk process \[ \Psi(t) = \sum_{j: t \leq T_j} \psi_j \] and the cumulative relative risks \(\Psi_i = \Psi(T_i).\)
Definition 11.7 Cox’s partial likelihood for the proportional hazards model13 is \[ L_{\mathrm{par}} (\Psi) = \prod_{i: \delta_i = 1} \frac{\psi_i}{\Psi_i} \] where the cumulative relative risks \(\Psi_i\) are given by Definition 11.6.
13 If there are no ties among the outcome observations in the data, the relative risks \(\psi_i\) are determined uniquely by the jumps of \(\Psi\). If there are ties, we cannot disentangle the individual relative risks completely from \(\Psi\) alone. We write \(L_{\mathrm{par}}\) as a function of \(\Psi\) anyway for notational convenience, remembering that when there are ties, \(\Psi\) must be annotated with information about how the jumps decompose into relative risk contributions for each individual.
The partial likelihood is a function of the relative risks only, but its definition may appear somewhat unmotivated. The following proposition shows how it is related to the survival likelihood.
Proposition 11.2 When \(\Lambda_i(t) = \psi_i \Lambda_0(t)\) the survival likelihood equals \[ L(\Psi, \Lambda_0) \coloneqq L_{\mathrm{par}} (\Psi) \left(\prod_{i} (\Psi_i \lambda_0(T_i))^{\delta_i}\right) e^{-\int_0^{\infty} \Psi(t) \lambda_0(t) \mathrm{d} t}. \]
Proof. First note that with \(\Psi(t) = \sum_{j: t \leq T_j} \psi_j,\) \[
\sum_{i} \Lambda_i(T_i) = \int_0^{\infty} \Psi(t) \lambda_0(t) \mathrm{d} t.
\] The survival likelihood can therefore be factorized as
\[
\begin{align*}
L(\Lambda) & = \prod_{i} \lambda_i(T_i)^{\delta_i} e^{-\Lambda_i(T_i)} \\
& = \left(\prod_{i: \delta_i = 1} (\psi_i \lambda_0(T_i))^{\delta_i}\right) e^{-\int_0^{\infty} \Psi(t) \lambda_0(t)
\mathrm{d} t} \\
& = \prod_{i: \delta_i = 1} \frac{\psi_i}{\Psi_i} \left(\prod_{i}
(\Psi_i \lambda_0(T_i))^{\delta_i}\right) e^{-\int_0^{\infty} \Psi(t) \lambda_0(t)
\mathrm{d} t} = L(\Psi, \Lambda_0).
\end{align*}
\]
Proposition 11.2 decomposes the likelihood into two factors \[ L(\Psi, \Lambda_0) = L_{\mathrm{par}} (\Psi) \tilde{L}(\Psi, \Lambda_0), \] where the second factor is \[ \tilde{L}(\Psi, \Lambda_0) = \left(\prod_{i} (\Psi_i \lambda_0(T_i))^{\delta_i}\right) e^{-\int_0^{\infty} \Psi(t) \Lambda_0(\mathrm{d} t)}. \] It is only this second factor that depends on the baseline hazard. A standard idea is then to compute the profile likelihood by maximizing \(\tilde{L}(\Psi, \Lambda_0)\) over \(\Lambda_0\) for a fixed cumulative relative risk process \(\Psi.\) This is, unfortunately, not directly possible, since \(\tilde{L}(\Psi, \Lambda_0)\) can be driven to \(+ \infty\) if we maximize over \(\Lambda_0\)—even if \(\Lambda_0\) is constrained to be given by a continuous hazard function \(\lambda_0\). It is, however, possible to renormalize the factor such that the profiling operation can be carried out in a meaningful way.
Given \(\varepsilon > 0\), the renormalized likelihood factor is defined as \[ \tilde L_\varepsilon(\Psi, \Lambda_0) = \prod_i \Big(\Psi_i\,\big[\Lambda_0(T_i) - \Lambda_0(T_i-\varepsilon)\big]\Big)^{\delta_i} e^{-\int_0^\infty \Psi(t)\,\Lambda_0(\mathrm{d} t)}. \] For a continuously differentiable \(\Lambda_0\), \[ \Lambda_0(T_i)-\Lambda_0(T_i-\varepsilon) = \varepsilon \lambda_0(T_i) + o(\varepsilon), \] whence if \(\lambda_0(T_i) > 0\) for all uncensored observations, \[ \tilde L_\varepsilon(\Psi, \Lambda_0) = \varepsilon^{n_u}\,\tilde L(\Psi, \Lambda_0)\,(1+o(1)), \] where \(n_u \coloneqq \sum_{i} \delta_i\) is the number of uncensored observations.
Proposition 11.3 If there are no ties in the data and \(\varepsilon > 0\) is smaller than the smallest gap between distinct observation times, then \[ \sup_{\Lambda_0} \tilde L_\varepsilon(\Psi, \Lambda_0) = e^{-n_u} \tag{11.16}\] and the supremum is attained in \[ \hat{\Lambda}_0(t) = \sum_{i: T_i \leq t} \frac{\delta_i}{\Psi_i}. \]
The supremum in Equation 11.16 is taken over all cumulative hazard functions \(\Lambda_0\), regarded as distribution functions for measures, including jump cumulative hazard functions of the form \[
\Lambda_0(t) = \sum_{i : T_i \leq t} \lambda_i
\] for \(\lambda_1, \ldots, \lambda_n \geq 0\). These are measures with point masses \(\lambda_i\) in \(T_i,\) that is, \(\Delta \Lambda_0(T_i) = \lambda_i.\)
Proof. The cumulative relative risk process \(\Psi\) is a left-continuous step function that is constant between observation times. So \(\Psi(s) = \Psi_i\) for \(s \in (T_i-\varepsilon,T_i]\).
With \(\omega_i = \Lambda_0(T_i) - \Lambda_0(T_i - \varepsilon)\), \[ \int_0^\infty \Psi(t) \,\Lambda_0(\mathrm{d}t) \ge \sum_i \Psi_i \omega_i, \] with equality exactly when \(\Lambda_0\) puts no mass outside the windows \((T_i-\varepsilon,T_i]\) on the set where \(\Psi > 0\).
This shows that \[ \log\tilde L_\varepsilon(\Psi, \Lambda_0) \le \sum_i \big(\delta_i\log(\Psi_i \omega_i) - \Psi_i \omega_i\big). \] The right-hand side is maximized term by term when \[ \hat{\omega}_i = \left\{ \begin{array}{ll} \Psi_i^{-1} & \quad \text{if } \delta_i = 1 \\ 0 & \quad \text{if } \delta_i = 0. \end{array}\right. \] This maximum is attained for \(\hat{\Lambda}_0\) with \[ \log\tilde L_\varepsilon(\Psi, \hat{\Lambda}_0) = - \sum_{i} \delta_i = - n_u. \]
We conclude from Proposition 11.2 and Proposition 11.3 that for \(\varepsilon\) small, \[ \begin{align*} L(\Psi, \Lambda_0) & = L_{\mathrm{par}} (\Psi) \tilde{L}(\Psi, \Lambda_0) \\ & \approx L_{\mathrm{par}} (\Psi) \tilde{L}_{\varepsilon}(\Psi, \Lambda_0) \ \varepsilon^{-n_u}, \end{align*} \] with the approximation holding for any fixed \(\Lambda_0(t) = \int_0^t \lambda_0(s) \mathrm{d} s\) as \(\varepsilon \to 0\). Fixing (a small) \(\varepsilon > 0\) and maximizing the renormalized likelihood factor over \(\Lambda_0\) gives \[ L_{\mathrm{par}} (\Psi) e^{-n_u} \varepsilon^{-n_u}, \] which is Cox’s partial likelihood up to the factor \(e^{-n_u} \varepsilon^{-n_u}\) that depends on neither the baseline nor the relative risks.
Note that the supremum in Equation 11.16 is, for any sufficiently small \(\varepsilon > 0\), also attained by any \(\Lambda_0^{\varepsilon}\) with hazard function \(\lambda_0^{\varepsilon}\) fulfilling \(\lambda_0^{\varepsilon}(s) = 0\) for \(s \not \in \cup_i (T_i-\varepsilon,T_i]\) and \[ \int_{T_i - \varepsilon}^{T_i} \lambda_0^{\varepsilon}(s) \mathrm{d} s = \hat{\omega}_i. \] Each of these converges pointwise to \(\hat{\Lambda}_0\) for \(\varepsilon \to 0\), and \(\hat{\Lambda}_0\) is sensibly regarded as the nonparametric maximum-likelihood estimate of \(\Lambda_0\).
11.4.2 Estimation and inference
The arguments following Proposition 11.2 and Proposition 11.3 show that for estimation of parameters that enter into the relative risks14, and not into \(\lambda_0,\) we can ignore the factor \(\tilde{L}(\Psi, \Lambda_0)\) in the likelihood and simply optimize the partial likelihood.
14 And are variation independent of \(\lambda_0.\)
In practice, the relative risks are typically linked to a vector of predictors \(X\) via an exponential function, that is, they are of the form \[
\psi_i = e^{X^T_i \beta} = e^{\eta_i}.
\tag{11.17}\] In this case we find that the partial log-likelihood equals
\[
\ell_{\mathrm{par}}(\beta) = \log L_{\mathrm{par}}(\beta) =
\sum_{i: \delta_i = 1} X_i^T \beta - \log \Bigg(\underbrace{\sum_{j: T_i \leq T_j} \psi_j}_{\Psi_i} \Bigg),
\tag{11.18}\] with corresponding score function \[
\begin{aligned}
\mathcal{U}_{\mathrm{par}}(\beta) & = \nabla \ell_{\mathrm{par}}(\beta) =
\sum_{i: \delta_i = 1} \left( X_i - \sum_{j: T_i \leq T_j} \frac{\psi_j}{\Psi_i} X_j \right) \\
& = \sum_{i: \delta_i = 1} \left( X_i - \bar{X}_i \right),
\end{aligned}
\tag{11.19}\] where \[
\bar{X}_i = \sum_{j: T_i \leq T_j} \frac{\psi_j}{\Psi_i} X_j.
\] The information function becomes
\[ \begin{aligned} \mathcal{J}_{\mathrm{par}}(\beta) & = - D \mathcal{U}_{\mathrm{par}}(\beta) \\ & = \sum_{i: \delta_i = 1} \left( \sum_{j: T_i \leq T_j} \frac{\psi_j}{\Psi_i} X_j X_j^T - \bar{X}_i \bar{X}_i^T \right) \\ & = \sum_{i: \delta_i = 1} \sum_{j: T_i \leq T_j} \frac{\psi_j}{\Psi_i} (X_j - \bar{X}_i)(X_j - \bar{X}_i)^T, \end{aligned} \tag{11.20}\] where we have used that the sum of the weights \(\psi_j / \Psi_i\) for \(j\) with \(T_i \leq T_j\) is one.
The estimate, \(\hat{\beta},\) is computed by maximizing \(\ell_{\mathrm{par}},\) and the observed information, \(\mathcal{J}_{\mathrm{par}}(\hat{\beta}),\) is used for quadratic approximations of \(\ell_{\mathrm{par}},\) for standard error estimates and Wald tests. In practice, the computations proceed just as if \(\ell_{\mathrm{par}}\) were a log-likelihood function.
The particular score test \[ \mathcal{U}_{\mathrm{par}}(0)^T \mathcal{J}_{\mathrm{par}}(0)^{-1} \mathcal{U}_{\mathrm{par}}(0) \] for testing \(\beta = 0\) reduces to a classical test known as the log-rank test in the special case with only a single categorical predictor and no ties.
Example 11.10 Cox’s proportional hazards model can be fitted to data via the function coxph() from the survival R package.
Code
prostate_cox <- coxph(
Surv(dtime + 0.5, status != "alive") ~ rx,
data = prostate,
)
prostate_cox |> tidy()# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 rx0.2 mg estrogen 0.0321 0.145 0.221 0.825
2 rx1.0 mg estrogen -0.386 0.157 -2.46 0.0140
3 rx5.0 mg estrogen 0.00291 0.146 0.0200 0.984 Code
anova(prostate_cox) |>
tidy() |>
filter(term == "rx") |>
knitr::kable(digits = 3)| term | logLik | statistic | df | p.value |
|---|---|---|---|---|
| rx | -2005.899 | 9.709 | 3 | 0.021 |
The parameter estimates are comparable to those obtained by survreg() with the Weibull distribution, by appropriately rescaling, since that model is parametrized as an AFT model. In particular, the Cox model regression parameters have the opposite sign. This is noteworthy, since a negative value of a Cox model regression parameter means that the corresponding predictor has a positive association with the survival time. That is, if \(\beta_j < 0\) then a unit increase of the \(j\)-th predictor multiplies the hazard by the factor \(e^{\beta_j} < 1,\) which results in a survival distribution with longer survival times.
Note that the (partial) likelihood-ratio test reported in Table 11.7 results in a very similar conclusion as using the Weibull baseline.
Once we have estimated \(\beta\) via the partial likelihood, we can estimate the baseline via Proposition 11.3. When \(\psi_i = e^{X_i^T \beta},\) Proposition 11.3 results in the estimator \[ \hat{\Lambda}_0(t) = \sum_{i: T_i \leq t} \frac{\delta_i}{\sum_{j:T_i \leq T_j} e^{X_j^T \hat{\beta}}} \] of the baseline cumulative hazard function in terms of the survival times and the estimate \(\hat{\beta}.\) This estimator is known as the Breslow estimator, which can be seen as an extension of the Nelson-Aalen estimator of the baseline cumulative hazard function.
Example 11.11 We continue Example 11.10 and compute the Breslow estimator of the baseline cumulative hazard using survfit(), see ?survfit.coxph for details. In fact, survfit() returns an estimate of the baseline survival function as \[
\hat{S}_0(t) = e^{-\hat{\Lambda}_0(t)}
\] where \(\hat{\Lambda}_0\) is a Breslow-type estimator.
Code
survfit(prostate_cox) |>
summary(data.frame = TRUE) |>
ggplot(aes(time, surv)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = "gray90") +
geom_step() +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability") 
Figure 11.7 shows the Breslow estimate including a pointwise 95% confidence band, which is also computed by the survfit() function. Figure 11.8 shows the fitted survival functions for the four different treatment groups, using the Breslow estimator of the baseline cumulative hazard function, compared to using a Weibull parametric model as baseline. The survival curves are based on the Cox model and in principle computed as \[
\hat{S}_k = e^{-e^{\hat{\beta}_k} \hat{\Lambda}_0}, \quad k = 0, 1, 2, 3.
\tag{11.21}\] Note that for Cox’s proportional hazards model, \(\beta_0 = 0\) always. That is, there is no intercept parameter. Since the baseline hazard is nonparametric, we cannot identify a separate intercept parameter in this model, and the intercept is effectively absorbed into the nonparametric baseline. As pointed out in the survfit() documentation, this also means that there is no such thing as the baseline hazard. What the baseline hazard is depends on a choice of baseline value of \(X,\) and thus also on how the predictors are encoded. In this example, the baseline computed by survfit() corresponds to the placebo treatment group, but in more complicated setups the baseline corresponds to a pseudo individual—see ?survfit.coxph.
It can be instructive to extract the estimated regression parameters and use (11.21) explicitly, but it is safer and more convenient to compute the estimated survival functions shown in Figure 11.8 using survfit() with a new dataset. The code below shows how, and it gives a survival function for each row of the new data.
Code
pred_frame <- tibble(
data = 1:4,
rx = factor(1:4, labels = levels(prostate$rx)),
strata = factor(1:4, labels = paste("rx=", rx, sep = ""))
)
survfit(prostate_cox, newdata = pred_frame) |>
summary(data.frame = TRUE) |>
left_join(pred_frame) |>
ggplot(aes(time, surv, color = strata)) +
geom_line(data = fit, linetype = 2) +
geom_step() +
ylim(0, 1) +
xlab("time (months)") +
ylab("survival probability") +
theme(legend.position = "top")
11.4.3 How ties are handled
Some of the results and procedures in the previous sections required no ties in the data. In this section we briefly discuss how to handle ties, which occur frequently in practice.
If there are ties among the survival times, the summation over \(j: T_i \leq T_j\) in the definition of the cumulative relative risks is still well defined. If it happens that \(T_i = T_k,\) the corresponding cumulative relative risks \(\Psi_i = \Psi_k\) also become identical. And if both observations are uncensored, they enter into the partial likelihood as the product \[ \frac{\psi_i}{\Psi_i} \frac{\psi_k}{\Psi_k} = \frac{\psi_i \psi_k}{\Psi_i^2}. \] This implicit way of handling ties is known as Breslow’s method. It treats ties as if they were really due to events that happened at the exact same time.
The renormalization argument extends to handle ties in this way. To see this, suppose there are \(d_k\) tied deaths at the \(k\)-th distinct time, where \(\Psi\) takes the value \(\Psi_{(k)}\). The corresponding window then contributes \((\Psi_{(k)}\omega_k)^{d_k} e^{-\Psi_{(k)}\omega_k}\), which is maximized at \(\omega_k = d_k/\Psi_{(k)}\). This gives \[ \sup_{\Lambda_0} \tilde L_\varepsilon = \prod_k d_k^{d_k}\, e^{-n_u}, \] which is still independent of \(\Psi\). The profile is therefore \(L_{\mathrm{par}}\) with Breslow’s tie handling, and the formula for \(\hat\Lambda_0\) is unchanged.
If the tie \(T_i = T_k\) was due to rounding or grouping, the contribution to the partial likelihood from those two observations should have been one of the two factors: \[ \frac{\psi_i \psi_k}{\Psi_i (\Psi_i - \psi_i)} \quad \text{or} \quad \frac{\psi_i \psi_k}{\Psi_k (\Psi_k - \psi_k)}. \tag{11.22}\]
We could try to break the tie, e.g., by letting the R function order() break ties in whatever arbitrary way it chooses depending on the sorting algorithm and the original order of the data. This naive solution will break all ties but in an unsystematic way, choosing one of the two factors above arbitrarily. Breaking ties in this way appears difficult to justify, and we should probably only break ties if we actually know and use the correct ordering prior to the rounding.
Efron’s method15 represents instead the “correct but unknown” factor by an approximation. Efron’s approximation for two ties is \[ \frac{\psi_i \psi_k}{\Psi_i (\Psi_i - (\psi_i + \psi_k)/2)} = \frac{\psi_i \psi_k}{\Psi_k (\Psi_k - (\psi_i + \psi_k)/2)} \] with similar formulas for more than two ties. The approximation is close to the geometric mean of the factors (11.22).
15 The default in coxph() for handling ties is Efron’s method.
Results rarely depend much on whether Efron’s or Breslow’s method (or even naive tie breaking) is used, but Efron’s method is usually preferred if the ties are due to a lack of precision and not to a true discrete nature of the survival times.
11.5 Survival residuals
Model fit can be assessed for survival regression models just as for generalized linear models by means of residuals. There are, of course, several different possible choices.
11.5.1 Cox-Snell residuals
Cox-Snell residuals are not really residuals in a classical sense, but a transformation of the observed survival times to a set of random variables whose distribution should be exponential if the model is correct. They are based on the following result.
Lemma 11.1 If \(T\) has continuous cumulative hazard function \(\Lambda\) then \(\Lambda(T)\) is exponentially distributed with rate 1.
Proof. The survival function \(S(t) = \exp(-\Lambda(t))\) is continuous, whence
\[
S(T) \sim \text{unif}([0,1]).
\] From this it follows that \[
\mathbf{P}(\Lambda(T) \leq t) = \mathbf{P}(S(T) \geq \exp(-t)) = 1 - \exp(-t)
\] for \(t \geq 0,\) or \(\Lambda(T) \sim \text{Exp}(1).\)
Definition 11.8 With \(\hat{\Lambda}_i\) an estimator of the cumulative hazard function for \(T_i\) the Cox-Snell residuals are defined as \[ \hat{\Lambda}_i(T_i) \tag{11.23}\] for \(i = 1, \ldots, n.\)
From the lemma above we conclude that if \((T_1, \delta_1), \ldots, (T_n, \delta_n)\) are independently right censored survival times with corresponding cumulative hazard functions \(\Lambda_1, \ldots, \Lambda_n,\) then \[ (\Lambda_1(T_1), \delta_1), \ldots, (\Lambda_n(T_n), \delta_n) \] are independently right censored i.i.d. exponentially distributed variables.
If there were no censoring, we would expect the Cox-Snell residuals to follow approximately an exponential distribution. With censoring the Cox-Snell residuals are expected to be approximately independently censored exponentially distributed. We use this in practice by estimating the cumulative hazard function for the Cox-Snell residuals using, e.g., the Nelson-Aalen estimator or the negative log of the Kaplan-Meier estimator. The resulting estimate is then compared to the cumulative hazard for the exponential distribution, which is \(\Lambda(t) = t.\) In practice we plot the estimated cumulative hazards against the Cox-Snell residuals, which should fall on a straight line with slope \(1.\)
11.5.2 Martingale residuals
The martingale residuals are defined as \[ \delta_i - \hat{\Lambda}_i(T_i), \] and are thus the status indicators minus the Cox-Snell residuals. Their usefulness is partly justified by the following result.
Lemma 11.2 If \(T = \min\{T^*, C\}\) with \(T^*\) and \(C\) independent and both having continuously differentiable survival functions it holds that \[ \mathbf{E}(\Lambda(T)) = \mathbf{P}(T^* \leq C), \] where \(\Lambda\) is the cumulative hazard function for \(T^*.\)
Proof. Let \(S\) and \(H\) denote the survival functions of \(T^*\) and \(C,\) respectively. By integration by parts \[ \begin{align*} \mathbf{E}(\Lambda(T)) & = - \int_0^{\infty} \Lambda(t) (SH)'(t) \mathrm{d} t = \int_0^{\infty} \Lambda'(t) S(t) H(t) \mathrm{d} t \\ & = \int_0^{\infty} f(t) H(t) \mathrm{d} t = \mathbf{P}(T^* \leq C). \end{align*} \]
The survival function for the minimum of the two independent random variables \(T^*\) and \(C\) is the product \(SH,\) whence the density of their distribution is the negative derivative \(-(SH)'.\)
Lemma 11.2 holds conditionally on \(X\) under independent censoring, and it follows, possibly using the tower property of conditional expectations, that with \(\delta = 1(T^* \leq C),\) \[ \mathbf{E} (\delta - \Lambda(T)) = 0. \tag{11.24}\]
Thus if \(\hat{\Lambda}_i\) were the true cumulative hazard function, the martingale residuals would have mean zero.
Now \(\hat{\Lambda}_i\) is only an estimate, but if the model fits the data, the martingale residuals should have approximately mean \(0.\) In the more advanced theory of survival analysis they also enjoy a martingale property, which explains their name. Here we will not pursue this, but just observe that we can investigate the model fit by investigating if the martingale residuals approximately have mean \(0.\) The residuals can be used much as regression residuals to judge overall fit as well as the fit of how the individual predictor variables enter into the model. It should be noted, though, that they by definition have a left-skewed distribution16 concentrated on \((-\infty, 1).\)
16 The martingale residuals are never approximately normally distributed!
In addition to the martingale residuals, we also define the squared-martingale residuals as \[ (\delta_i - \hat{\Lambda}_i(T_i))^2 - \hat{\Lambda}_i(T_i). \] Exercise 11.8 shows that \(\mathbf{E}((\delta - \Lambda(T))^2) = \mathbf{E}(\Lambda(T))\) under independent censoring, if \(\Lambda\) is the cumulative hazard function for \(T^*\), and the squared-martingale residuals should thus also roughly have mean zero.
Example 11.12 As for other regression models, we can use augment() from the broom package to extract fitted values (the linear predictor) and residuals. Figure 11.9 shows the martingale residuals and squared-martingale residuals plotted against the treatment (which is the only predictor in the model at this point). Since the dataset is rather small, it makes sense here to use geom_beeswarm() from the ggbeeswarm package to visualize the full distribution of the residuals as well as the individual points. Note the upper cap at \(1\) of all martingale residual distributions.
Code
prostate_cox_aug <- augment(prostate_cox, data = prostate) |>
mutate(
Lambda = (status != "alive") - .resid, # Cox-Snell residuals
v = .resid^2 - Lambda # Squared-martingale residuals
)
p1 <- ggplot(prostate_cox_aug, aes(rx, .resid)) +
ggbeeswarm::geom_beeswarm(alpha = 0.4) +
ylab("Martingale residuals") +
theme(
axis.text.x = element_text(angle = 15, hjust = 1, vjust = 1)
)
p2 <- ggplot(prostate_cox_aug, aes(rx, v)) +
ggbeeswarm::geom_beeswarm(alpha = 0.4) +
geom_point(
stat = "summary",
fun = mean,
color = "blue",
size = 2.5
) +
geom_errorbar(
stat = "summary",
fun.data = mean_se,
fun.args = list(mult = 1.96),
color = "blue",
width = 0.15,
linewidth = 0.7
) +
ylab("Squared-martingale residuals") +
theme(
axis.text.x = element_text(angle = 15, hjust = 1, vjust = 1)
)
gridExtra::grid.arrange(p1, p2, ncol = 2)
Since the model only includes the single categorical predictor rx, the martingale residual plot is not informative about misspecifications, because the empirical mean of the martingale residuals is forced to be exactly zero within each group by the estimation procedure. The plot is only shown here to demonstrate what the distribution of martingale residuals looks like. The plot of the squared-martingale residuals in Figure 11.9 is, however, informative as a check of the proportional hazards assumption. The figure includes means and confidence intervals for each group of squared-martingale residuals, and these intervals should include zero for all groups, as they do in the figure.
We can also plot the residuals against another potential predictor, such as age, as in Figure 11.10. We see a clear association between age and the residuals suggesting that age (unsurprisingly) holds predictive information about the survival times.
age, which was not included as a predictor in the model.
We want to include age in the model, and we will first try to include it as an additive term.
Code
prostate_cox_add <- coxph(
Surv(dtime + 0.5, status != "alive") ~ age + rx,
data = prostate,
)Figure 11.11 shows two residual plots. The first compares the Cox-Snell residuals to the exponential distribution, and the second is a standard type of residual plot, this time of the martingale residuals against the fitted values. The fitted values are the linear predictors. For the computation of the Cox-Snell residuals, we used that type.predict = "expected" in the call of augment() gives \(\hat{\Lambda}_i(T_i)\) as the predicted values (.fitted column in the output).
Code
p1 <- augment(prostate_cox_add, data = prostate, type.predict = "expected") |>
survfit(Surv(.fitted, status != "alive") ~ 1, data = _) |>
summary(data.frame = TRUE) |>
ggplot(aes(time, cumhaz)) +
geom_abline(slope = 1, color = "blue") +
geom_point() +
xlab("Cox-Snell residuals") +
ylab("Cumulative hazard")
p2 <- augment(prostate_cox_add) |>
ggplot(aes(.fitted, .resid)) +
geom_point(alpha = 0.4) +
geom_smooth() +
xlab("Fitted values") +
ylab("Martingale residuals")
gridExtra::grid.arrange(p1, p2, ncol = 2)
We can see minor model misspecifications from the residual plots in Figure 11.11. First, there are deviations from a straight line in the Cox-Snell plot, suggesting some form of model misspecification. The smoothed curve in the martingale residual plot also suggests some minor model misspecification. We could suspect that age is not correctly modeled.
11.5.3 Deviance residuals
Alternative residuals include the deviance residuals and the Schoenfeld residuals. We will not treat the latter, but the former are based on the following observation. The parametric log-likelihood with a continuous baseline as well as the logarithm of \(L_{\mathrm{par}}(\Psi) \tilde{L}_{\varepsilon}(\Psi, \Lambda_0)\), evaluated in jump cumulative hazard functions, can be written as \[ \sum_{i=1}^n \underbrace{\delta_i \log (\psi_i \Lambda_{0,i}) - \psi_i \Lambda_{0,i}}_{\text{Poisson log-likelihood term}} + \sum_{i=1}^n \delta_i \log \left(\frac{\lambda_{0,i}}{\Lambda_{0,i}}\right), \] where \(\Lambda_{0,i} = \Lambda_0(T_i)\), \(\lambda_{0,i} = \lambda_0(T_i)\) in the parametric case and \(\lambda_{0,i} = \Delta \Lambda_0(T_i)\) in the jump case. The first term is identical to the Poisson log-likelihood, and the second term does not depend on parameters that enter into the relative risks. Based on this we define the deviance residuals in terms of the Poisson deviances \[ d(\delta_i, \hat{\Lambda}_{i}) = 2\left( \delta_i \log (\delta_i/\hat{\Lambda}_{i}) - \delta_i + \hat{\Lambda}_{i}\right), \] where \(\hat{\Lambda}_{i} = \hat{\psi}_i \hat{\Lambda}_{0,i}.\)
The result is the deviance residuals \[ \text{sign}(\delta_i - \hat{\Lambda}_i) \sqrt{d(\delta_i, \hat{\Lambda}_i)} = \left\{\begin{array}{ll} - \sqrt{2 \hat{\Lambda}_i} & \quad \text{if } \delta_i = 0 \\ \text{sign}( 1 - \hat{\Lambda}_i)\sqrt{ 2( \hat{\Lambda}_i - 1 - \log \hat{\Lambda}_i)} & \quad \text{if } \delta_i = 1 \end{array} \right. \]
Observe that the deviance residuals can be expressed as \[ \text{sign}(\hat{r}_i) \sqrt{-2(\hat{r}_i + \delta_i \log (\delta_i - \hat{r}_i))} \] with \(\hat{r}_i = \delta_i - \hat{\Lambda}_i\) the martingale residuals. Compared to the martingale residuals, the distribution of the deviance residuals is less skewed.
11.5.4 Prostate cancer continued
We conclude this section with a final analysis of the prostate cancer data based on the insights we have obtained from the previous sections. We include both age and treatment rx in the model, now with age nonlinearly expanded using natural cubic splines and rx interacting (linearly) with age.
Code
prostate_cox_int <- coxph(
Surv(dtime + 0.5, status != "alive") ~ ns(age, df = 4) + age*rx,
data = prostate,
)We first reconsider the residual plots in Figure 11.12. In contrast to Figure 11.11 there are no indications in these residual plots of any model misspecification.
Code
p1 <- augment(prostate_cox_int, data = prostate, type.predict = "expected") |>
survfit(Surv(.fitted, status != "alive") ~ 1, data = _) |>
summary(data.frame = TRUE) |>
ggplot(aes(time, cumhaz)) +
geom_abline(slope = 1, color = "blue") +
geom_point() +
xlab("Cox-Snell residuals") +
ylab("Cumulative hazard")
p2 <- augment(prostate_cox_int) |>
ggplot(aes(.fitted, .resid)) +
geom_point() +
geom_smooth() +
xlab("Fitted values") +
ylab("Martingale residuals")
gridExtra::grid.arrange(p1, p2, ncol = 2)
To better understand the interaction model we will compute and visualize the estimated relative risks \[ e^{\hat{\eta}} \] for some combinations of treatment and age as well as the corresponding fitted survival curves.
Code
pred_frame <- expand.grid(
age = c(55, 62.5, 70, 77.5),
rx = factor(1:4, labels = levels(prostate$rx))
)
pred_frame <- mutate(
pred_frame,
risk = predict(
prostate_cox_int,
newdata = pred_frame,
type = "risk"
),
risk = risk / risk[1],
data = 1:16
) | age | rx | risk | data |
|---|---|---|---|
| 55.0 | placebo | 1.000 | 1 |
| 62.5 | placebo | 0.989 | 2 |
| 70.0 | placebo | 0.894 | 3 |
| 77.5 | placebo | 1.313 | 4 |
| 55.0 | 0.2 mg estrogen | 1.438 | 5 |
| 62.5 | 0.2 mg estrogen | 1.231 | 6 |
| 70.0 | 0.2 mg estrogen | 0.963 | 7 |
| 77.5 | 0.2 mg estrogen | 1.224 | 8 |
| 55.0 | 1.0 mg estrogen | 0.391 | 9 |
| 62.5 | 1.0 mg estrogen | 0.499 | 10 |
| 70.0 | 1.0 mg estrogen | 0.582 | 11 |
| 77.5 | 1.0 mg estrogen | 1.102 | 12 |
| 55.0 | 5.0 mg estrogen | 0.494 | 13 |
| 62.5 | 5.0 mg estrogen | 0.679 | 14 |
| 70.0 | 5.0 mg estrogen | 0.852 | 15 |
| 77.5 | 5.0 mg estrogen | 1.738 | 16 |
The computations above give a data frame with estimated relative risks for 16 combinations of age and treatment. Since relative risks are only determined up to a common factor, we have normalized them so that the placebo treatment for age 55 is the reference. Table 11.8 shows the table of the resulting relative risks. Recall that the smaller the relative risk is, the smaller is the risk of dying and the longer is the survival time.
Figure 11.13 shows a plot of the estimated relative risks, from which the nonlinear and non-additive nature of the model is clear. For the high treatment dosage there is a notable change: it reduces the risk relative to placebo for patients younger than about 70 years, while it increases the risk for older patients.
Figure 11.14 shows the fitted survival functions corresponding to the relative risks in Table 11.8. It gives a different perspective on how the estimated survival curves change with age, and how the interaction effect makes the treatment with the high dosage appear beneficial for younger patients while harmful for older patients.
Code
survfit(prostate_cox_int, newdata = pred_frame) |>
summary(data.frame = TRUE) |>
left_join(pred_frame) |>
ggplot(aes(time, surv, color = rx)) +
geom_step(
data = summary(prostate_KM, data.frame = TRUE),
aes(x = time + 0.5, color = NULL),
linetype = 2
) +
geom_step() +
facet_wrap("age", ncol = 4) +
theme(legend.position = "top")
We have not quantified the uncertainty in the above figures. One way to globally test if there is a treatment effect in this interaction model is via a likelihood-ratio test (using the partial likelihood). Table 11.9 shows that the \(p\)-value is about \(0.003\) and thus fairly small. Based on this test we therefore reject the hypothesis that there is no treatment effect.
Code
prostate_cox_age <- coxph(
Surv(dtime + 0.5, status != "alive") ~ ns(age, df = 4),
data = prostate,
)
anova(prostate_cox_age, prostate_cox_int) |>
tidy() |>
knitr::kable()| term | logLik | statistic | df | p.value |
|---|---|---|---|---|
| ~ ns(age, df = 4) | -1999.203 | NA | NA | NA |
| ~ ns(age, df = 4) + age * rx | -1989.345 | 19.7167 | 6 | 0.0031098 |
Today, estrogen treatment of prostate cancer has largely been replaced by better alternatives. Though it is documented to be effective against the cancer, a serious known side effect is increased risk of cardiovascular disease. This knowledge is in line with our findings for this particular study.
11.6 Discrete time survival models
We have in previous sections mostly considered survival distributions that are continuous, that is, distributions on \([0, \infty)\) with a continuous survival function \(S.\) It was typically even differentiable with \(-S'\) the density of the distribution, in which case the hazard function is well defined. In this section we will briefly touch on discrete survival distributions, and almost exclusively on distributions with point masses in the observations \(T_1, \ldots, T_n.\) We assume throughout this section that there are no ties.
Definition 11.9 With observations \(T_1, \ldots, T_n\) and \(P\) a probability measure on \([0, \infty)\) with survival function \(S,\) the empirical likelihood with right censoring is \[ L(P) = \prod_i P(T_i)^{\delta_i} S(T_i)^{1-\delta_i}. \]
The factors in the empirical likelihood are simply the probabilities \[ P(T_i)^{\delta_i} S(T_i)^{1-\delta_i} = \left\{ \begin{array}{ll} \mathbf{P}(T^* = T_i) & \quad \text{if } \delta_i = 1 \\ \mathbf{P}(T^* > T_i) & \quad \text{if } \delta_i = 0 \end{array} \right. \]
Recall also that for the discrete probability measure \(P,\) the survival function is \[ S(t) = \mathbf{P}(T^* > t) = \sum_{i: T_i > t} P(T_i). \] We will show in this section that the Kaplan-Meier estimator can be derived as the maximizer of the empirical likelihood.
11.6.1 Discrete hazards
In the previous section we only introduced hazards for continuous distributions. The discrete hazards are by definition the conditional probabilities \[ \lambda_i \coloneqq \frac{P(T_i)}{S(T_i-)} = \mathbf{P}(T^* = T_i \mid T^* \geq T_i). \] Note that the discrete hazards are probabilities and thus constrained to lie in the interval \([0,1].\)
Lemma 11.3 For a discrete probability measure with point masses in \(T_1, \ldots, T_n\) and corresponding hazards \(\lambda_i\) the survival function is \[ S(t) = \prod_{i: T_i \leq t} (1 - \lambda_i). \]
Proof. Note that \(S(T_{(j)}-) = S(T_{(j-1)})\) for \(j \geq 2\) and that it is sufficient to show the identity for \(t = T_{(j)},\) \(j = 1, \ldots, n.\)
The proof is by induction. For \(T_{(1)}\) we have \(S(T_{(1)}-) = 1\) and \[ S(T_{(1)}) = 1 - P(T_{(1)}) = 1 - \lambda_{(1)}, \] which gives the induction start.
For the induction step, observe that the definition of the hazard implies that \[ \begin{align*} S(T_{(j)}) & = S(T_{(j)}-) - P(T_{(j)}) \\ & = S(T_{(j)}-) - S(T_{(j)}-)\lambda_{(j)} \\ & = S(T_{(j-1)})(1-\lambda_{(j)}), \end{align*} \] which gives the induction step.
Lemma 11.3 shows that the discrete hazards determine the survival function and thus the distribution, but in a slightly different way than for continuous distributions. As a direct corollary we can express the empirical likelihood in terms of the discrete hazards.
Corollary 11.1 With \(Y_i = Y(T_i) = \sum_{j} 1(T_i \leq T_j)\) the number of individuals at risk at time \(T_i,\) the empirical likelihood equals \[ L(P) = \prod_{i} \lambda_i^{\delta_i} (1 - \lambda_i)^{Y_i - \delta_i}. \tag{11.25}\]
Proof. From the definition of the hazards and Lemma 11.3 \[ P(T_i) = \lambda_i S(T_i-) = \lambda_i \prod_{j: T_j < T_i} (1 - \lambda_j). \]
We then plug this expression into the empirical likelihood function and interchange the product order to get \[ \begin{align*} L(P) & = \prod_i \lambda_i ^{\delta_i} \prod_{j: T_j < T_i} (1 - \lambda_j)^{\delta_i} \prod_{j: T_j \leq T_i} (1 - \lambda_j)^{1-\delta_i} \\ & = \prod_{i} \lambda_i^{\delta_i} (1 - \lambda_i)^{Y_i - \delta_i}. \end{align*} \]
11.6.2 Empirical likelihood estimators
In terms of the empirical likelihood we define the nonparametric maximum-likelihood estimator (NPMLE) of the survival distribution as the maximizer of \(L.\) To compute it we maximize \(L\) as a function of the discrete hazards. To this end we need the following lemma.
Lemma 11.4 With \(y \geq 1\) the function \[ \lambda \mapsto \lambda (1 - \lambda)^{y - 1} \] attains a unique maximum over \([0, 1]\) in \(1/y.\)
Proof. With \(f(\lambda) = \lambda (1 - \lambda)^{y - 1}\) we have \[ \begin{align*} f'(\lambda) & = (1 - \lambda)^{y - 1} - \lambda (y-1) (1 - \lambda)^{y - 2} \\ & = (1 - \lambda)^{y - 2}(1 - \lambda - \lambda (y - 1)) \\ & = (1 - \lambda)^{y - 2} (1 - \lambda y). \end{align*} \] For \(y = 1,\) the function is the identity, which attains its maximum in \(1 = 1/y.\) For \(y > 1\) we have \(f(\lambda) \geq 0\) for \(\lambda \in [0,1],\) \(f(0) = f(1) = 0\) and the function has a unique stationary point in \((0, 1)\) at \(1/y\) with \(f(1/y) > 0.\) Thus \(1/y\) is the unique maximizer in \([0,1]\) for \(y \geq 1.\)
Theorem 11.2 The NPMLE of the hazards—the hazards for the probability measure that maximizes the empirical likelihood—is given by \[ \hat{\lambda}_i = \frac{\delta_i}{Y_i}. \]
Proof. The factors of \(L\) in (11.25) equal \[ \lambda_i^{\delta_i} (1 - \lambda_i)^{Y_i - \delta_i} = \left\{ \begin{array}{cc} (1-\lambda_i)^{Y_i} & \quad \text{if } \delta_i = 0 \\ \lambda_i (1-\lambda_i)^{Y_i - 1} & \quad \text{if } \delta_i = 1 \\ \end{array} \right. \] In the first case, when \(\delta_i = 0,\) the factor is maximized in \(0,\) and in the second case, when \(\delta_i = 1,\) the factor is maximized in \(1/Y_i.\)
By plugging in the NPMLE of the hazards in the formula in Lemma 11.3 for the survival function we find that the NPMLE of the survival function is \[ \hat{S}(t) = \prod_{i: T_i \leq t} \left(1 - \frac{\delta_i}{Y_i}\right). \] We recognize this as the Kaplan-Meier estimator.
In terms of the discrete hazards the corresponding jump cumulative hazard function is \[
\Lambda(t) = \sum_{i : T_i \leq t} \lambda_i.
\] If we plug in the NPMLE of the hazards in this formula, we get \[
\hat{\Lambda}(t) = \sum_{i : T_i \leq t} \frac{\delta_i}{Y_i},
\] which we recognize as the Nelson-Aalen estimator. The previously observed
discrepancy between the Kaplan-Meier and Nelson-Aalen estimators was due to the relation in (11.9) assuming a continuous distribution.
The two main points from this section are:
- The Kaplan-Meier estimator has a justification as a nonparametric maximum-likelihood estimator
- When we use the correct relation between hazards and the survival function for discrete distributions, as given by Lemma 11.3, there is no discrepancy between the Kaplan-Meier and Nelson-Aalen estimators.
11.6.3 The discrete proportional hazards model
It is legitimate to consider the proportional hazards model for discrete probability distributions, and use the empirical likelihood function \(L\) for profiling. Contrary to the derivations in Section 11.4, the maximizer is then not explicit. Moreover, the proportionality assumption for the discrete model does not really provide an approximation of the proportionality assumption for the continuous model.
The appropriate generalization of the proportional hazards model to the discrete case is actually the model \[ \lambda_i = 1 - (1 - \lambda_{0,i})^{\psi_i}. \tag{11.26}\]
The empirical likelihood becomes \[ L = \prod_{i} ( 1 - (1 - \lambda_{0,i})^{\psi_i})^{\delta_i} (1 - \lambda_{0,i})^{\Psi_i - \psi_i \delta_i}, \] and similar arguments as in the proof of Theorem 11.2 will show that the NPMLE of \(\lambda_{0,i}\) for given relative risks is \[ \hat{\lambda}_{0,i} = \delta_i\left(1 - \left(1 - \frac{\psi_i}{\Psi_i}\right)^{1/\psi_i}\right). \] The corresponding profile likelihood is \[ \prod_{i: \delta_i = 1} \frac{\psi_i}{\Psi_i} \left(1 - \frac{\psi_i}{\Psi_i}\right)^{\Psi_i/\psi_i - \delta_i}. \tag{11.27}\]
The discrete proportional hazards model (11.26) is arguably the “correct” generalization of the continuous proportional hazards model, see Jacobsen (1984), and using the empirical likelihood yields the profile likelihood (11.27). The resulting profile likelihood is not Cox’s partial likelihood, but it is for most practical purposes very close. In the literature, Cox’s partial likelihood has caught on as the default choice for semiparametric estimation of parameters that enter the relative risks.
Jacobsen, Martin. 1984. “Maximum Likelihood Estimation in the Multiplicative Intensity Model: A Survey.” International Statistical Review 52 (2): 193–207.
Exercises
Exercise 11.1 Show that if \(T\) has a log-logistic distribution with density \(f = f_0\) then \[ e^{\eta} T \] has density \(f_{\eta}.\) Show that if \(T\) has a log-logistic distribution with density \(f_{\eta}\) then \[ \log T - \eta \] has a logistic distribution.
Exercise 11.2 Show that without censoring the Kaplan-Meier estimator equals the empirical survival function \[ \hat{S}(t) = \frac{1}{n} \sum_{i=1}^n 1(T^*_i > t). \]
Exercise 11.3 In this exercise you should consider the prostate dataset from the RwR package and this time include the variable hx in the analysis. The variable hx indicates if the patient has a history of cardiovascular disease. You should ignore the treatment variable rx.
- Fit a proportional hazards model using the two variables
ageandhxas predictors. - Make one or more residual plots for the model you have fitted. Does it fit the data?
- Test the hypothesis that survival is not associated with
hxgivenage.
Exercise 11.4 Let \(T_1^*, \ldots, T^*_n\) be independent survival times for \(n\) patients and let \(X_i \in \{0, 1\}\) indicate if the \(i\)-th patient is treated \((X_i = 1)\) or untreated \((X_i = 0).\) We are interested in estimating the treatment effect, which we summarize in the following parameter of interest \[ \rho := \frac{\text{median}(T_1^* \mid X_1 = 1)}{\text{median}(T_1^* \mid X_1 = 0)}, \] that is, \(\rho\) is the ratio of the median survival time for treated vs. untreated patients.
In this exercise it is assumed that \(T_i^*\) follows a Weibull distribution with shape parameter \(\gamma > 0\) and scale parameter \(\alpha_{X_i},\) that is, the hazard rate is \[ \lambda_i(t) = \alpha_{X_i} \gamma t^{\gamma-1}. \]
Show that \[ \rho = \sqrt[\gamma]{\frac{\alpha_0}{\alpha_1}} \] in terms of the Weibull parameters.
Show that the Weibull model can be formulated as a proportional hazards model in terms of the linear predictor \[ \eta_i = \beta_0 + \beta_1 X_i \] and relate the parameters \(\beta_0\) and \(\beta_1\) to \(\alpha_0\) and \(\alpha_1.\) Express \(\rho\) in terms of the \(\beta_0\) and \(\beta_1\) parameters.
In the remaining part of this exercise we will study two different estimators of \(\rho\) on simulated data. To this end introduce i.i.d. censoring times \(C_1, \ldots, C_n\) and the censored observed survival times \[ T_i = \min\{T_i^*, C_i\}. \]
The data set testWei available from the R package RwR is a simulated data set with \(n = 200,\) \(X_1 = \ldots = X_{100} = 0\) and \(X_{101} = \ldots = X_{200} = 1,\) with \(\gamma = 2,\) \(\beta_0 = 0\) and \(\beta_1 = -1\) (so that \(\rho \approx 1.649\)). The censoring times are uniformly distributed on \([0, 4].\) The time column contains the \(T_i\)-s and the status column is the status indicator (\(1\) is uncensored, \(0\) is censored).
Compute and plot the Kaplan-Meier estimates of the survival function within each of the two treatment groups for the simulated
testWeidata set.Use the Weibull model with fixed \(\gamma = 2\) to estimate \(\beta_0\) and \(\beta_1.\) Compare the resulting estimated survival functions from the Weibull model with the Kaplan-Meier estimates.
Compute the MLE \(\hat{\rho}_{\mathrm{MLE}}\) of \(\rho\) based on the Weibull model (with fixed \(\gamma = 2\)) and compute a 95% confidence interval for \(\rho.\) Compute also, based on the Kaplan-Meier estimates, a direct estimate, \(\hat{\rho}_{\mathrm{KM}},\) of \(\rho.\)
Write an R program that simulates a data set (including the censoring) using the same setup as for the
testWeidata set described above. Simulate \(B = 200\) data sets and compute for each data set the estimates \(\hat{\rho}_{\mathrm{KM}}\) and \(\hat{\rho}_{\mathrm{MLE}}\) as above and compare their distributions.
Hint: Use the R function rweibull() to simulate from the Weibull distribution, but pay attention to the parametrization.
Exercise 11.5 For this exercise we consider the data set lung on lung cancer survival available from the R package RwR. It contains four variables:
| Variable | Description |
|---|---|
| time | survival time in days |
| status | censored = 1, dead = 2 |
| age | age in years |
| sex | male = 1, female = 2 |
- Plot the Kaplan-Meier estimates of the survival functions for males and females separately.
- Use the Cox proportional hazards model to investigate if the survival time depends upon age.
- Check if the Cox proportional hazards model fits the data.
Exercise 11.6 This exercise is about the proportional hazards model with cumulative hazard function \[ \Lambda_i(t) = e^{\eta_i} \Lambda_0(t) \] for the \(i\)-th individual and linear predictor \[ \eta_i = \beta_0 + \beta_{\texttt{age}} X_{i, \texttt{age}} + \beta_{\texttt{sex}} X_{i, \texttt{sex}}. \] Thus the two predictor variables are the age and sex of the individual.
You will need the data set simSurv from the R package RwR that contains the four variables, surv, status, age and sex being the censored survival times, the status indicator, the age and the indicator of gender (1 is male), respectively, of the individuals. There is a total of 200 individuals.
- Estimate, using Cox’s proportional hazards model, the parameters \(\beta_{\texttt{age}}\) and \(\beta_{\texttt{sex}}.\) Compute a 95% confidence interval for \(\beta_{\texttt{sex}}\) and determine whether males have a higher risk of dying than females.
For the remaining questions the baseline cumulative hazard function is assumed to be \[ \Lambda_0(t) = e^{\gamma t} - 1, \quad t \geq 0. \] The cumulative hazard function for the \(i\)-th individual is thus \[ \Lambda_i(t) = e^{\eta_i} ( e^{\gamma t} - 1). \]
We assume that we have independent censored observations \((T_1, \delta_1), \ldots, (T_n, \delta_n)\) with \(T_i \geq 0\) the \(i\)-th censored survival time and \(\delta_i\) the indicator of whether the \(i\)-th observation is uncensored.
Compute the baseline hazard rate \(\lambda_0(t)\) and show that the log-likelihood function is \[ \sum_{i} \delta_i \eta_i + \delta_i \log (e^{\gamma T_i} - 1) - e^{\eta_i} (e^{\gamma T_i} - 1) + \delta_i \log \frac{\gamma e^{\gamma T_i}}{e^{\gamma T_i} - 1}. \]
Plot the profile log-likelihood for the \(\gamma\)-parameter. Compute the MLE of \(\gamma,\) \(\beta_0,\) \(\beta_{\texttt{age}}\) and \(\beta_{\texttt{sex}}.\)
Hint: Identify the first three terms of the log-likelihood as the log-likelihood for a Poisson generalized linear model. Use the glm()-function for the profiling.
- Use the glm framework to compute a 95% confidence interval for \(\beta_{\texttt{sex}}\) by plugging in the MLE, \(\hat{\gamma},\) of \(\gamma\) and pretending that it is fixed and known. Compute next a 95% confidence interval using non-parametric bootstrapping. Compare the results.
Exercise 11.7 This exercise deals with the numerical algorithms associated with the optimization of Cox’s partial likelihood. Consider censored survival data \((T_i, \delta_i, X_i)\) for \(i = 1, \ldots, n\) with \(T_i\) the possibly censored survival time, \(\delta_i\) the indicator of whether \(T_i\) is uncensored, and \(X_i \in \mathbb{R}\) a univariate real valued predictor. The proportional hazards model is given by the individual relative risks \[ \psi_i = e^{\beta X_i} \] for an unknown parameter \(\beta.\) In the following we will make the assumption that the survival times are ordered and with no ties, that is \[ T_1 < T_2 < \ldots < T_n. \]
- Show that Cox’s partial log-likelihood equals \[ \ell(\beta) = \sum_{i: \delta_i = 1} \beta X_i - \log \left(\sum_{j \geq i} e^{\beta X_j} \right), \] and compute the first and second derivative of \(\ell\) w.r.t. \(\beta.\)
The estimator \(\hat{\beta}\) is the maximizer of \(\ell\) provided that \(\ell\) attains a maximum. It then solves the equation \[ \ell'(\beta) = 0. \]
- Show that \(\ell\) is concave. Find a condition on the \(X_i\)-s that ensures that \(\ell\) has a unique maximum.
For the remaining part of this problem you first need to simulate a data set from the model above. You may experiment with different choices of \(\beta,\) but try also \(\beta = 0.\)
- Implement Newton’s algorithm for computing \(\hat{\beta}.\) Apply the algorithm with start value \(\beta_{1} = 0\) to the simulated data set and investigate the values of \(\beta_m,\) \(\ell(\beta_m),\) \(\ell'(\beta_m)\) and \(\ell''(\beta_m)\) in the iterations of the algorithm.
A naive implementation of the computation of e.g. \(\ell\) scales like \(n^2.\) Make sure that your implementation scales like \(n.\)
Compute and plot an estimate of the baseline cumulative hazard function.
Compute \(\hat{\beta}\) using the
coxph()function in R. Compute also a standard 95% confidence interval for \(\beta.\) Compare this interval with a (partial) likelihood interval based on \(\ell.\)
Exercise 11.8 This exercise is on theoretical properties of the survival analysis martingale residuals. Let \(T^*\) and \(C\) be independent random variables on \([0,\infty),\) whose distributions have continuous densities \(f\) and \(h,\) respectively. The corresponding survival functions are denoted \(S\) and \(H.\) Note that \(f = - S'\) and \(h = - H'.\) Then \(T = \min\{T^*, C\}\) denotes the observed (censored) survival time and \(\delta = 1(T^* \leq C)\) is the status indicator.
- Show that for any function \(K: [0,\infty) \times \{0,1\} \to \mathbb{R},\) continuously differentiable in the first variable, it holds that \[ \begin{align*} \mathbf{E}( K(T, \delta)) & = \int_0^\infty \Big(K(0,1) + (K(s, 0) - K(s,1)) S(s) \\ & \quad + \int_0^s \partial_tK(t, 1) S(t) \mathrm{d} t \Big) h(s) \mathrm{d} s \end{align*} \] provided that the expectation is well defined. Use this to show that17 \[ \mathbf{E} (\delta - \Lambda(T)) = 0 \qquad \text{and} \qquad \mathbf{V}(\delta - \Lambda(T)) = \mathbf{E}(\Lambda(T)) \] where \(\Lambda\) denotes the cumulative hazard function for the distribution given by \(f.\)
17 That the expectation is \(0\) is derived in the main text, but a slightly different derivation is provided here. The variance formula is used in Section 11.5 to justify the squared-martingale residuals.
Recall that for the proportional hazards model, the \(i\)-th individual has cumulative hazard function \(\psi_i \Lambda_0\) where \(\psi_i\) is the individual relative risk and \(\Lambda_0\) is the baseline cumulative hazard function.
Assume that \(\hat{\psi}_i\) is the estimate of the \(i\)-th relative risk and that \(\hat{\Lambda}_0\) is the estimate of the baseline. With \(\hat{\Lambda}_i = \hat{\psi}_i \hat{\Lambda}_0(T_i)\) the fitted cumulative hazard function for the \(i\)-th individual evaluated in \(T_i,\) \[ \hat{r}_i = \delta_i - \hat{\Lambda}_i \] is the \(i\)-th martingale residual. Let \(\hat{d}_i\) denote the \(i\)-th deviance residual18.
18 The deviance residuals are sometimes regarded as distribution stabilized transforms of the martingale residuals. This exercise gives some insight into this perspective.
- Give an argument for considering the following standardization of the martingale residuals: \[ \tilde{r}_i = \frac{\hat{r}_i}{\sqrt{\hat{\Lambda}_i}}. \]
- Show that if \(\delta_i = 1\) then \[ \hat{d}_i^2 = \tilde{r}_i^2 + o((1 - \hat{\Lambda}_i)^2). \]
- Explore graphically how good the approximation \[ \hat{d}_i \approx \tilde{r}_i \] is as a function of \(\hat{\Lambda}_i\) when \(\delta_i = 1.\) Does a similar approximation hold if \(\delta_i = 0\)?