Introduction

Sir David Cox observed that if the proportional hazards assumption holds (or, is assumed to hold) then it is possible to estimate the effect parameter(s), denoted ${\displaystyle \beta _{i}}$ below, without any consideration of the full hazard function. This approach to survival data is called application of the Cox proportional hazards model,^[2] sometimes abbreviated to Cox model or to proportional hazards model.^[3] However, Cox also noted that biological interpretation of the proportional hazards assumption can be quite tricky.^[4][5]

Let X_i = (X_i1, … , X_ip) be the realized values of the p covariates for subject i. The hazard function for the Cox proportional hazards model has the form

${\displaystyle {\begin{aligned}\lambda (t|X_{i})&=\lambda _{0}(t)\exp(\beta _{1}X_{i1}+\cdots +\beta _{p}X_{ip})\\&=\lambda _{0}(t)\exp(X_{i}\cdot \beta )\end{aligned}}}$

This expression gives the hazard function at time t for subject i with covariate vector (explanatory variables) X_i. Note that between subjects, the baseline hazard ${\displaystyle \lambda _{0}(t)}$ is identical (has no dependency on i). The only difference between subjects' hazards comes from the baseline scaling factor ${\displaystyle \exp(X_{i}\cdot \beta )}$.

Why it is called "proportional"

To start, suppose we only have a single covariate, ${\displaystyle x}$, and therefore a single coefficient, ${\displaystyle \beta _{1}}$. Consider the effect of increasing ${\displaystyle x}$ by 1:

${\displaystyle {\begin{aligned}\lambda (t|x+1)&=\lambda _{0}(t)\exp(\beta _{1}(x+1))\\&=\lambda _{0}(t)\exp(\beta _{1}x+\beta _{1})\\&={\Bigl (}\lambda _{0}(t)\exp(\beta _{1}x){\Bigr )}\exp(\beta _{1})\\&=\lambda (t|x)\exp(\beta _{1})\end{aligned}}}$

We can see that increasing a covariate by 1 scales the original hazard by the constant ${\displaystyle \exp(\beta _{1})}$. Rearranging things slightly, we see that:

${\displaystyle {\frac {\lambda (t|x+1)}{\lambda (t|x)}}=\exp(\beta _{1})}$

The right-hand-side is constant over time (no term has a ${\displaystyle t}$ in it). This relationship, ${\displaystyle x/y={\text{constant}}}$, is called a proportional relationship.

More generally, consider two subjects, i and j, with covariates ${\displaystyle X_{i}}$ and ${\displaystyle X_{j}}$ respectively. Consider the ratio of their hazards:

${\displaystyle {\begin{aligned}{\frac {\lambda (t|X_{i})}{\lambda (t|X_{j})}}&={\frac {\lambda _{0}(t)\exp(X_{i}\cdot \beta )}{\lambda _{0}(t)\exp(X_{j}\cdot \beta )}}\\&={\frac {{\cancel {\lambda _{0}(t)}}\exp(X_{i}\cdot \beta )}{{\cancel {\lambda _{0}(t)}}\exp(X_{j}\cdot \beta )}}\\&=\exp((X_{i}-X_{j})\cdot \beta )\end{aligned}}}$

The right-hand-side isn't dependent on time, as the only time-dependent factor, ${\displaystyle \lambda _{0}(t)}$, was cancelled out. Thus the ratio of hazards of two subjects is a constant, i.e. the hazards are proportional.

Absence of an intercept term

Often there is an intercept term (also called a constant term or bias term) used in regression models. The Cox model lacks one because the baseline hazard, ${\displaystyle \lambda _{0}(t)}$, takes the place of it. Let's see what would happen if we did include an intercept term anyways, denoted ${\displaystyle \beta _{0}}$:

${\displaystyle {\begin{aligned}\lambda (t|X_{i})&=\lambda _{0}(t)\exp(\beta _{1}X_{i1}+\cdots +\beta _{p}X_{ip}+\beta _{0})\\&=\lambda _{0}(t)\exp(X_{i}\cdot \beta )\exp(\beta _{0})\\&=\left(\exp(\beta _{0})\lambda _{0}(t)\right)\exp(X_{i}\cdot \beta )\\&=\lambda _{0}^{*}(t)\exp(X_{i}\cdot \beta )\end{aligned}}}$

where we've redefined ${\displaystyle \exp(\beta _{0})\lambda _{0}(t)}$ to be a new baseline hazard, ${\displaystyle \lambda _{0}^{*}(t)}$. Thus, the baseline hazard incorporates all parts of the hazard that are not dependent on the subjects' covariates, which includes any intercept term (which is constant for all subjects, by definition).

Likelihood for unique times

The Cox partial likelihood, shown below, is obtained by using Breslow's estimate of the baseline hazard function, plugging it into the full likelihood and then observing that the result is a product of two factors. The first factor is the partial likelihood shown below, in which the baseline hazard has "canceled out". The second factor is free of the regression coefficients and depends on the data only through the censoring pattern. The effect of covariates estimated by any proportional hazards model can thus be reported as hazard ratios.

The likelihood of the event to be observed occurring for subject i at time Y_i can be written as:

${\displaystyle L_{i}(\beta )={\frac {\lambda (Y_{i}\mid X_{i})}{\sum _{j:Y_{j}\geq Y_{i}}\lambda (Y_{i}\mid X_{j})}}={\frac {\lambda _{0}(Y_{i})\theta _{i}}{\sum _{j:Y_{j}\geq Y_{i}}\lambda _{0}(Y_{i})\theta _{j}}}={\frac {\theta _{i}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}},}$

where θ_j = exp(X_j ⋅ β) and the summation is over the set of subjects j where the event has not occurred before time Y_i (including subject i itself). Obviously 0 < L_i(β) ≤ 1. This is a partial likelihood: the effect of the covariates can be estimated without the need to model the change of the hazard over time.

Treating the subjects as if they were statistically independent of each other, the joint probability of all realized events^[6] is the following partial likelihood, where the occurrence of the event is indicated by C_i = 1:

${\displaystyle L(\beta )=\prod _{i:C_{i}=1}L_{i}(\beta ).}$

The corresponding log partial likelihood is

${\displaystyle \ell (\beta )=\sum _{i:C_{i}=1}\left(X_{i}\cdot \beta -\log \sum _{j:Y_{j}\geq Y_{i}}\theta _{j}\right).}$

This function can be maximized over β to produce maximum partial likelihood estimates of the model parameters.

The partial score function is

${\displaystyle \ell ^{\prime }(\beta )=\sum _{i:C_{i}=1}\left(X_{i}-{\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}\right),}$

and the Hessian matrix of the partial log likelihood is

${\displaystyle \ell ^{\prime \prime }(\beta )=-\sum _{i:C_{i}=1}\left({\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}X_{j}^{\prime }}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}-{\frac {\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}\right]\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}^{\prime }\right]}{\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}\right]^{2}}}\right).}$

Using this score function and Hessian matrix, the partial likelihood can be maximized using the Newton-Raphson algorithm. The inverse of the Hessian matrix, evaluated at the estimate of β, can be used as an approximate variance-covariance matrix for the estimate, and used to produce approximate standard errors for the regression coefficients.

Likelihood when there exist tied times

Several approaches have been proposed to handle situations in which there are ties in the time data. Breslow's method describes the approach in which the procedure described above is used unmodified, even when ties are present. An alternative approach that is considered to give better results is Efron's method.^[7] Let t_j denote the unique times, let H_j denote the set of indices i such that Y_i = t_j and C_i = 1, and let m_j = |H_j|. Efron's approach maximizes the following partial likelihood.

${\displaystyle L(\beta )=\prod _{j}{\frac {\prod _{i\in H_{j}}\theta _{i}}{\prod _{\ell =0}^{m_{j}-1}\left[\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}\right]}}.}$

The corresponding log partial likelihood is

${\displaystyle \ell (\beta )=\sum _{j}\left(\sum _{i\in H_{j}}X_{i}\cdot \beta -\sum _{\ell =0}^{m_{j}-1}\log \left(\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}\right)\right),}$

the score function is

${\displaystyle \ell ^{\prime }(\beta )=\sum _{j}\left(\sum _{i\in H_{j}}X_{i}-\sum _{\ell =0}^{m_{j}-1}{\frac {\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}X_{i}}{\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}}}\right),}$

and the Hessian matrix is

${\displaystyle \ell ^{\prime \prime }(\beta )=-\sum _{j}\sum _{\ell =0}^{m_{j}-1}\left({\frac {\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}X_{i}^{\prime }-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}X_{i}X_{i}^{\prime }}{\phi _{j,\ell ,m_{j}}}}-{\frac {Z_{j,\ell ,m_{j}}Z_{j,\ell ,m_{j}}^{\prime }}{\phi _{j,\ell ,m_{j}}^{2}}}\right),}$

where

${\displaystyle \phi _{j,\ell ,m_{j}}=\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}}$

${\displaystyle Z_{j,\ell ,m_{j}}=\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}-{\frac {\ell }{m_{j}}}\sum _{i\in H_{j}}\theta _{i}X_{i}.}$

Note that when H_j is empty (all observations with time t_j are censored), the summands in these expressions are treated as zero.

Examples

Below are some worked examples of the Cox model in practice.