Cox Proportional-Hazards Regression
The Kaplan–Meier estimator describes survival in a single group or compares two groups, but it cannot adjust for the influence of multiple covariates simultaneously. Cox proportional-hazards regression (Cox, 1972) extends survival analysis to the multivariate setting, modelling the relationship between survival time and one or more predictor variables while accounting for censoring. It is the most widely used regression model for time-to-event outcomes in clinical research.
The Hazard Function
The foundation of the Cox model is the hazard function $h(t)$ — the instantaneous rate of experiencing the event at time $t$, conditional on having survived up to that moment:
$$ h(t) = \lim_{\Delta t \to 0} \frac{\Pr(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} $$The hazard is not a probability; it can exceed 1. It describes the risk intensity at each instant. A rising hazard function means the event becomes more likely over time (as in post-operative mortality); a falling one means the risk decreases (as in many infectious diseases after the acute phase).
The Cox Model
The Cox model specifies how the hazard function changes with a set of covariates $X_1, X_2, \ldots, X_k$:
$$ h(t \mid \mathbf{X}) = h_0(t) \cdot \exp(b_1 X_1 + b_2 X_2 + \cdots + b_k X_k) $$The function $h_0(t)$ is the baseline hazard — the hazard when all covariates are zero. It can take any shape and is left completely unspecified. The exponential term multiplies the baseline hazard by a factor that depends on the covariate values but not on time. This structure gives the model two distinctive properties:
- Semi-parametric: the baseline hazard is non-parametric (no distribution is assumed for the event times), while the covariate effects are estimated parametrically.
- Proportional hazards: the ratio of hazards between any two subjects is constant over time, equal to $\exp(\mathbf{b}'(\mathbf{X}_1 - \mathbf{X}_2))$.
Taking logarithms, the model becomes:
$$ \ln\!\left[\frac{h(t \mid \mathbf{X})}{h_0(t)}\right] = b_1 X_1 + b_2 X_2 + \cdots + b_k X_k $$The left-hand side is the log hazard ratio relative to the baseline — a linear function of the covariates, analogous to the linear predictor in logistic regression.
Partial Likelihood Estimation
Cox (1972, 1975) showed that the regression coefficients $b_1, \ldots, b_k$ can be estimated without specifying the baseline hazard $h_0(t)$. At each observed event time $t_j$, only the ordering of subjects matters: the subject who actually experienced the event contributes a term that is the ratio of their exponential score to the sum of scores over all subjects still at risk at $t_j$:
$$ \mathcal{L}_j = \frac{\exp(\mathbf{b}' \mathbf{X}_j)}{\displaystyle\sum_{i \in \mathcal{R}(t_j)} \exp(\mathbf{b}' \mathbf{X}_i)} $$where $\mathcal{R}(t_j)$ is the set of subjects at risk just before $t_j$. Maximising the product of these terms over all event times — the partial likelihood — yields estimates of the regression coefficients. Because $h_0(t)$ cancels entirely, censored observations contribute only through the risk set $\mathcal{R}(t_j)$: they remain in the denominator until they are censored or experience the event, correctly using the information they provide without introducing bias.
Hazard Ratios
The primary output of a Cox model is a set of hazard ratios $\exp(b_j)$ with 95% confidence intervals:
- For a continuous covariate, $\exp(b_j)$ is the hazard ratio associated with a one-unit increase in $X_j$, holding all other covariates constant. If age is entered in years and $\exp(b) = 1.04$, a one-year increase in age is associated with a 4% higher instantaneous event rate.
- For a binary covariate (coded 0/1), $\exp(b_j)$ directly compares the two groups. If $\exp(b) = 0.55$ for treatment versus control, the treated group has a 45% lower hazard at every point in time.
- For a categorical covariate with more than two levels, one level serves as the reference and a separate hazard ratio is estimated for each other level.
The forest plot above illustrates a typical Cox model output. Each row shows the estimated hazard ratio and its 95% CI for one predictor. The reference line at HR = 1 marks the null hypothesis of no effect: estimates to the right indicate increased hazard, estimates to the left indicate reduced hazard. Confidence intervals that do not cross the reference line are statistically significant.
The Proportional-Hazards Assumption
The fundamental assumption of the Cox model is that the hazard ratio between any two subjects remains constant over time. This is not a requirement that the baseline hazard be constant — it can rise, fall, or fluctuate in any way — only that the multiplicative effect of each covariate does not change as time passes. A treatment whose benefit diminishes over time, or a biomarker that predicts early but not late recurrence, violates this assumption.
The assumption can be examined in two complementary ways:
- Log-log plot: plot $\log(-\log \hat{S}(t))$ against $\log(t)$ for each group. Under proportional hazards the curves should be approximately parallel. Systematic convergence or divergence indicates non-proportionality.
- Schoenfeld residuals: plot the scaled Schoenfeld residuals for each covariate against time. A horizontal scatter around zero indicates proportional hazards; a trend over time indicates a time-varying effect (Grambsch & Therneau, 1994).
When the assumption is violated, possible remedies include stratifying on the offending covariate (allowing a separate baseline hazard per stratum) or adding a time-interaction term $b_j(t) \cdot X_j$.
Model Fit and Discrimination
Three test statistics assess the overall significance of the model — whether the predictors jointly improve on the null model containing no covariates:
- The likelihood ratio test compares $-2[\ln \mathcal{L}(\mathbf{0}) - \ln \mathcal{L}(\hat{\mathbf{b}})]$ to a $\chi^2$ distribution with $k$ degrees of freedom. It is generally the most reliable of the three.
- The Wald test evaluates each coefficient as $(b_j / \text{SE}_{b_j})^2$, chi-squared with 1 df. It can be unreliable for large coefficients.
- The score test (log-rank test generalisation) is computed at the null and is numerically equivalent to the log-rank test when there is a single binary covariate.
Harrell's C-index (concordance index) measures how well the model discriminates between subjects who experience the event sooner versus later (Harrell et al., 1996). It is the probability that, for a randomly chosen pair of subjects where one has a shorter survival time, that subject also has a higher predicted hazard. A C-index of 0.5 indicates no discrimination beyond chance; 1.0 is perfect discrimination. Values above 0.70 are generally considered useful in clinical prediction models. MedCalc computes the C-index confidence interval using the method of Pencina & D'Agostino (2004).
Predicted Survival
For a subject with covariate vector $\mathbf{X}$, the predicted survival function is:
$$ \hat{S}(t \mid \mathbf{X}) = \hat{S}_0(t)^{\exp(\hat{\mathbf{b}}' \mathbf{X})} $$where $\hat{S}_0(t) = \exp(-\hat{H}_0(t))$ is the estimated baseline survival function. In practice, the prognostic index $\text{PI} = \hat{b}_1 X_1 + \cdots + \hat{b}_k X_k$ summarises the covariate effect for an individual: a higher PI corresponds to higher predicted hazard and lower predicted survival.
Sample Size
The key determinant of precision in a Cox model is the number of events, not the number of subjects. The rule of thumb recommended by Peduzzi et al. (1995) requires at least 10 events per predictor variable (EPV ≥ 10). With fewer events per variable, coefficient estimates become unstable, confidence intervals widen substantially, and stepwise selection procedures are unreliable. The minimum recommended total is 100 events.
The required sample size can be computed as $N = 10k/p$, where $k$ is the number of predictors and $p$ is the expected proportion of subjects experiencing the event during follow-up.
Example: A study plans to fit a Cox model with 5 predictors. If 40% of subjects are expected to experience the event, the required sample size is $N = 10 \times 5 / 0.40 = 125$ subjects, yielding 50 events — just at the 10 EPV threshold. Increasing to $N = 250$ gives 100 events (20 EPV) and substantially more reliable estimates.
MedCalc's Cox proportional-hazards regression procedure estimates the model coefficients and hazard ratios with 95% CIs, reports the likelihood ratio, Wald, and score tests for overall model significance, computes Harrell's C-index with CI, and produces the predicted survival curve for a subject with specified covariate values.
See the Kaplan–Meier survival analysis chapter for the non-parametric foundation, including the log-rank test and median survival time.