Proportional hazards models are a class of survival models in statistics. Survival models relate the time that passes, before some event occurs, to one or more covariates that may be associated with that quantity of time. In a proportional hazards model, the unique effect of a unit increase in a covariate is multiplicative with respect to the hazard rate. The hazard rate at time is the probability per short time dt that an event will occur between and given that up to time no event has occurred yet. For example, taking a drug may halve one's hazard rate for a stroke occurring, or, changing the material from which a manufactured component is constructed, may double its hazard rate for failure. Other types of survival models such as accelerated failure time models do not exhibit proportional hazards. The accelerated failure time model describes a situation where the biological or mechanical life history of an event is accelerated (or decelerated).
Background
Survival models can be viewed as consisting of two parts: the underlying baseline hazard function, often denoted , describing how the risk of event per time unit changes over time at baseline levels of covariates; and the effect parameters, describing how the hazard varies in response to explanatory covariates. A typical medical example would include covariates such as treatment assignment, as well as patient characteristics such as age at start of study, gender, and the presence of other diseases at start of study, in order to reduce variability and/or control for confounding.
The proportional hazards condition states that covariates are multiplicatively related to the hazard. In the simplest case of stationary coefficients, for example, a treatment with a drug may, say, halve a subject's hazard at any given time , while the baseline hazard may vary. Note however, that this does not double the lifetime of the subject; the precise effect of the covariates on the lifetime depends on the type of . The covariate is not restricted to binary predictors; in the case of a continuous covariate , it is typically assumed that the hazard responds exponentially; each unit increase in results in proportional scaling of the hazard.
The Cox model
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 below, without any consideration of the full hazard function. This approach to survival data is called application of the Cox proportional hazards model, sometimes abbreviated to Cox model or to proportional hazards model. However, Cox also noted that biological interpretation of the proportional hazards assumption can be quite tricky.
Let be the realized values of the p covariates for subject i. The hazard function for the Cox proportional hazards model has the form
This expression gives the hazard function at time t for subject i with covariate vector (explanatory variables) Xi. Note that between subjects, the baseline hazard is identical (has no dependency on i). The only difference between subjects' hazards comes from the baseline scaling factor .
Why it is called "proportional"
To start, suppose we only have a single covariate, , and therefore a single coefficient, . Our model looks like:
Consider the effect of increasing by 1:
We can see that increasing a covariate by 1 scales the original hazard by the constant . Rearranging things slightly, we see that:
The right-hand-side is constant over time (no term has a in it). This relationship, , is called a proportional relationship.
More generally, consider two subjects, i and j, with covariates and respectively. Consider the ratio of their hazards:
The right-hand-side isn't dependent on time, as the only time-dependent factor, , 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, , takes the place of it. Let's see what would happen if we did include an intercept term anyways, denoted :
where we've redefined to be a new baseline hazard, . 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). In other words, adding an intercept term would make the model unidentifiable.
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". It is simply the probability for subjects to have experienced events in the order that they actually have occurred, given the set of times of occurrences and given the subjects' covariates. 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.
To calculate the partial likelihood, the probability for the order of events, let us index the M samples for which events have already occurred by increasing time of occurrence, Y1 < Y2 < ... < YM. Covariates of all other subjects for which no event has occurred get indices M+1,.., N. The partial likelihood can be factorized into one factor for each event that has occurred. The i 'th factor is the probability that out of all subjects (i,i+1,..., N) for which no event has occurred before time Yi, the one that actually occurred at time Yi is the event for subject i:
where ) and the summation is over the set of subjects j where the event has not occurred before time Yi (including subject i itself). Obviously 0 < Li(β) ≤ 1.
Treating the subjects as statistically independent of each other, the partial likelihood for the order of events "Each failure contributes to the likelihood function", Cox (1972), page 191. is
where the subjects for which an event has occurred are indicated by Ci = 1 and all others by Ci = 0. The corresponding log partial likelihood is
where we have written using the indexing introduced above in a more general way, as . Crucially, the effect of the covariates can be estimated without the need to specify the hazard function over time. The partial likelihood can be maximized over β to produce maximum partial likelihood estimates of the model parameters.
The partial score function is
and the Hessian matrix of the partial log likelihood is
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. Let tj denote the unique times, let Hj denote the set of indices i such that Yi = tj and Ci = 1, and let mj = |Hj|. Efron's approach maximizes the following partial likelihood.
The corresponding log partial likelihood is
the score function is
and the Hessian matrix is
where
Note that when Hj is empty (all observations with time tj are censored), the summands in these expressions are treated as zero.
Examples
Below are some worked examples of the Cox model in practice.
A single binary covariate
Suppose the endpoint we are interested in is patient survival during a 5-year observation period after a surgery. Patients can die within the 5-year period, and we record when they died, or patients can live past 5 years, and we only record that they lived past 5 years. The surgery was performed at one of two hospitals, A or B, and we would like to know if the hospital location is associated with 5-year survival. Specifically, we would like to know the relative increase (or decrease) in hazard from a surgery performed at hospital A compared to hospital B. Provided is some (fake) data, where each row represents a patient: T is how long the patient was observed for before death or 5 years (measured in months), and C denotes if the patient died in the 5-year period. We have encoded the hospital as a binary variable denoted X: 1 if from hospital A, 0 from hospital B.
hospital X T C B 0 60 False B 0 32 True B 0 60 False B 0 60 False B 0 60 False A 1 4 True A 1 18 True A 1 60 False A 1 9 True A 1 31 True A 1 53 True A 1 17 True
Our single-covariate Cox proportional model looks like the following, with representing the hospital's effect, and i indexing each patient:
Using statistical software, we can estimate to be 2.12. The hazard ratio is the exponential of this value, . To see why, consider the ratio of hazards, specifically:
Thus, the hazard ratio of hospital A to hospital B is . Putting aside statistical significance for a moment, we can make a statement saying that patients in hospital A are associated with a 8.3x higher risk of death occurring in any short period of time compared to hospital B.
There are important caveats to mention about the interpretation:
a 8.3x higher risk of death does not mean that 8.3x more patients will die in hospital A: survival analysis examines how quickly events occur, not simply whether they occur. More specifically, "risk of death" is a measure of a rate. A rate has units, like meters per second. However, a relative rate does not: a bicycle can go two times faster than another bicycle (the reference bicycle), without specifying any units. Likewise, the risk of death (comparable to the speed of a bike) in hospital A is 8.3 times higher (faster) than the risk of death in hospital B (the reference group). the inverse quantity, is the hazard ratio of hospital B relative to hospital A. We haven't made any inferences about probabilities of survival between the hospitals. This is because we would need an estimate of the baseline hazard rate, , as well as our estimate. However, standard estimation of the Cox proportional hazard model does not directly estimate the baseline hazard rate. Because we have ignored the only time varying component of the model, the baseline hazard rate, our estimate is timescale-invariant. For example, if we had measured time in years instead of months, we would get the same estimate. It is tempting to say that the hospital caused the difference in hazards between the two groups, but since our study is not causal (that is, we do not know how the data was generated), we stick with terminology like "associated".
A single continuous covariate To demonstrate a less traditional use case of survival analysis, the next example will be an economics question: what is the relationship between a company's price-to-earnings ratio (P/E) on their first IPO anniversary and their future survival? More specifically, if we consider a company's "birth event" to be their first IPO anniversary, and any bankruptcy, sale, going private, etc. as a "death" event the company, we'd like to know the influence of the companies' P/E ratio at their "birth" (first IPO anniversary) on their survival.