Cox Proportional Hazards › Formula

The Cox model makes two intertwined claims at once. First, covariates multiply the hazard through an exponential function of a linear predictor. Second, the effect of every covariate is the same at every time. Both claims appear in the same equation, so let us unpack them separately before we introduce any optimization machinery.

We start from the hazard function, which we introduced in the previous chapter. The hazard at time $t$ is the instantaneous rate of an event among survivors:

$$h(t \mid x) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t, x)}{\Delta t}.$$

Here $T$ is the random event time, $x \in \mathbb{R}^p$ is a covariate vector for the individual, and $\Delta t$ is a small time increment. The hazard tells us how many events per unit time occur among those still at risk; the survival function $S(t)$, in contrast, tells us what fraction survive past $t$. They are related by $h(t) = -d \log S(t)/dt$.

Why model the hazard rather than the survival probability directly? A hazard is a nonnegative rate, so a multiplicative model keeps every fitted hazard nonnegative automatically. It also has a natural "events per unit time among survivors" interpretation, which is often exactly the quantity a clinician, engineer, or product manager wants to reason about. Cox's claim is that this rate factors into a piece that depends only on $t$ and a piece that depends only on $x$:

$$h(t \mid x) = h_0(t) \exp(\beta^\top x), \qquad \beta \in \mathbb{R}^p.$$

The function $h_0(t) \ge 0$ is the baseline hazard, defined as the hazard of an individual with $x = 0$, so every new symbol here is accounted for. The exponential does two jobs at once: it maps any real linear predictor to a positive multiplier, and it turns additive changes in $x$ into multiplicative changes in the hazard. A one-unit increase in covariate $x_p$ multiplies the hazard by $e^{\beta_p}$ at every time.

where:

• $t$: time since the start of follow-up, measured in the same units as $T$
• $x \in \mathbb{R}^p$: the covariate vector for the individual
• $\beta \in \mathbb{R}^p$: the coefficient vector, with $\beta_p$ the coefficient on covariate $x_p$
• $h_0(t)$: the baseline hazard, the hazard of an individual whose covariates are all zero
• $\exp(\beta^\top x)$: the individual-specific multiplier that scales the baseline hazard up or down

The model becomes clear when we compare two individuals at the same instant. Taking the ratio of their hazards gives

$$\frac{h(t \mid x_i)}{h(t \mid x_j)} = \frac{h_0(t) \exp(\beta^\top x_i)}{h_0(t) \exp(\beta^\top x_j)} = \exp\bigl(\beta^\top (x_i - x_j)\bigr).$$

The baseline $h_0(t)$ cancels. That cancellation is the entire reason the Cox model is practical: the ratio of hazards depends only on $\beta$ and covariate differences, not on time. This is exactly the proportional hazards assumption, and the quantity $\exp(\beta^\top (x_i - x_j))$ is called the hazard ratio.

where:

• $x_i$: the covariate vector for individual $i$
• $x_j$: the covariate vector for individual $j$
• $\beta$: the coefficient vector shared by all individuals
• $h(t \mid x_i)/h(t \mid x_j)$: the ratio of the two individuals' hazards at time $t$

For fixed covariates and times with positive baseline hazard, for a single binary covariate with $\beta = \log 0.5$, a unit change in the covariate halves the hazard at every point in time. The "proportional" in the name refers to this ratio: the two hazard curves may have any common shape $h_0(t)$, but their vertical ratio is fixed.

Let us be precise about what a hazard ratio is not. For the Cox model the survival function has a very specific form. In the previous chapter we saw that $S(t) = \exp(-H(t))$, where $H(t) = \int_0^t h(s) \, ds$ is the cumulative hazard. Substituting the Cox hazard gives

$$S(t \mid x) = \exp\left(-H_0(t) \exp(\beta^\top x)\right) = S_0(t)^{\exp(\beta^\top x)},$$

where $S_0(t) = \exp(-H_0(t))$ is the baseline survival function and $H_0(t) = \int_0^t h_0(s)\, ds$ is the baseline cumulative hazard.

where:

• $t$: time since the start of follow-up
• $x$: the covariate vector for the individual
• $h_0(s)$: the baseline hazard at time $s$
• $H_0(t) = \int_0^t h_0(s)\, ds$: the baseline cumulative hazard, the total accumulated hazard up to time $t$
• $S_0(t) = \exp(-H_0(t))$: the baseline survival function, the survival probability when all covariates are zero
• $\exp(\beta^\top x)$: the individual-specific multiplier that scales the cumulative hazard
• $S(t \mid x)$: the survival probability of an individual with covariates $x$

The ratio of survival probabilities between two individuals is therefore

$$\frac{S(t \mid x_1)}{S(t \mid x_0)} = \exp\Bigl(-H_0(t) \bigl(e^{\beta^\top x_1} - e^{\beta^\top x_0}\bigr)\Bigr),$$

which depends on $t$ through $H_0(t)$.

where:

• $x_1, x_0$: the covariate vectors of the two individuals being compared
• $H_0(t)$: the baseline cumulative hazard up to time $t$
• $\beta$: the coefficient vector shared by all individuals
• $S(t \mid x_1)/S(t \mid x_0)$: the ratio of the two individuals' survival probabilities at time $t$

A hazard ratio of $0.5$ does not mean "survival is 50% higher" or "the event is twice as likely"; it means "the instantaneous event rate is half as large at every moment." The survival gap can change over time. If the cumulative baseline hazard grows without bound, both survival probabilities tend to zero and their difference eventually vanishes; finite follow-up or a finite limiting cumulative hazard need not show that full pattern. Use the survival curves to explain how a hazard ratio differs from a survival-probability comparison.

The proportional hazards equation describes the model, but it still contains the unknown function $h_0(t)$. If we tried to maximize an ordinary likelihood, we would need to specify and estimate that function at every time. Cox's second insight avoids this entirely by asking a conditional question at each event time. The object that makes the question precise is the risk set. At any time $t$, the risk set $R(t)$ is the set of individuals who have not yet experienced an event and have not yet been censored. Equivalently,

$$R(t) = \{\, j \in \{1, \dots, n\} : \tilde T_j \ge t \,\}.$$

Here $\tilde T_j=\min(T_j,C_j)$ is the observed duration, $T_j$ is the latent event time and $C_j$ is the censoring time. Everyone enters at time zero in this implementation. Individuals censored at an event time are included in the risk set immediately before that time, with events processed before censoring at the same recorded time.

where:

• $t$: a time point at which we evaluate who is still under observation
• $j$: an index running over the $n$ subjects in the study
• $\tilde T_j$: the observed event or censoring duration of subject $j$
• $R(t)$: the risk set, the set of subjects still at risk (not yet failed and not yet censored) at time $t$

A censored subject remains in every risk set up to their censoring time, even though they never contribute a numerator term. This is why survival data cannot be reduced to ordinary classification: dropping censored rows would discard the very information that tells us who was still at risk when later events occurred.

Suppose the observed event times are $t_{(1)} < t_{(2)} < \cdots < t_{(K)}$, with $K$ the number of distinct event times (for now we assume no ties). Let $i_k$ be the individual who failed at $t_{(k)}$. Cox's idea is to ask a conditional question at each event time: given that some individual in the risk set failed at $t_{(k)}$, what is the probability that it was individual $i_k$ in particular? The probability is

$$P(i_k \mid t_{(k)}) = \frac{h(t_{(k)} \mid x_{i_k})}{h(t_{(k)} \mid x_{i_k}) + \sum_{j \in R(t_{(k)}) \setminus \{i_k\}} h(t_{(k)} \mid x_j)}.$$

where:

• $t_{(k)}$: the $k$-th distinct event time
• $i_k$: the individual who fails at $t_{(k)}$
• $x_{i_k}$: the covariate vector of the failing individual
• $R(t_{(k)})$: the risk set at event time $t_{(k)}$
• $h(t_{(k)} \mid x_j)$: the hazard of individual $j$ at time $t_{(k)}$
• $P(i_k \mid t_{(k)})$: the conditional probability that the failing individual was $i_k$, given that exactly one individual in the risk set failed at $t_{(k)}$

Pulling out the factor $h_0(t_{(k)})$, which appears in every term in the numerator and denominator, gives

$$P(i_k \mid t_{(k)}) = \frac{h_0(t_{(k)}) e^{\beta^\top x_{i_k}}}{h_0(t_{(k)}) \sum_{j \in R(t_{(k)})} e^{\beta^\top x_j}} = \frac{e^{\beta^\top x_{i_k}}}{\sum_{j \in R(t_{(k)})} e^{\beta^\top x_j}}.$$

The baseline hazard vanished. This is the partial likelihood contribution at time $t_{(k)}$. It is a "partial" likelihood because it models only the identity of the failure within the risk set, not the absolute time at which the failure occurred. Multiplying contributions across event times gives the partial likelihood

$$L(\beta) = \prod_{k=1}^{K} \frac{\exp(\beta^\top x_{i_k})}{\sum_{j \in R(t_{(k)})} \exp(\beta^\top x_j)}.$$

where:

• $K$: the number of distinct event times
• $t_{(k)}$: the $k$-th distinct event time, with $k = 1, \dots, K$
• $i_k$: the individual who fails at $t_{(k)}$
• $x_{i_k}$: the covariate vector of the failing individual at $t_{(k)}$
• $R(t_{(k)})$: the risk set at event time $t_{(k)}$
• $\beta$: the coefficient vector to be estimated
• $L(\beta)$: the partial likelihood, a function of $\beta$ after the baseline hazard has cancelled

Taking logs turns the product into a sum and the fractions into differences:

$$\ell(\beta) = \log L(\beta) = \sum_{k=1}^{K} \left[ \beta^\top x_{i_k} - \log \sum_{j \in R(t_{(k)})} \exp(\beta^\top x_j) \right].$$

Each event contributes a numerator term (the log weight of the individual that failed) minus a denominator term (the log of the total weight of everyone who was at risk). Maximizing $\ell$ therefore means finding the $\beta$ that makes the observed failures look as likely as possible given who was at risk when they happened. Increasing $\beta_p$, holding the other coefficients fixed, raises the relative weight of higher-$x_p$ individuals. An event’s derivative with respect to this coefficient is positive when the failing individual exceeds the current weighted risk-set mean. Information also depends on covariate spread and weights: a one-subject risk set contributes none, but the number of alternatives alone does not determine the strength of evidence.

The partial likelihood is not the joint density of the observed data: it models the identities of failures conditional on the risk sets, rather than the full distribution of observed event and censoring times. Cox (1975) outlined why partial likelihood supports likelihood-like large-sample inference. Under suitable regularity and sampling conditions, the Cox partial-likelihood estimator is consistent and asymptotically normal (Andersen and Gill, 1982), and inverse observed information provides an approximate model-based covariance. To derive the score equations, treat $\ell$ exactly as you would treat a log-likelihood.

In matrix form, if we let $d_i \in \{0,1\}$ be the event indicator for individual $i$ and $Y_j(t) = \mathbb{1}\{\tilde T_j \ge t\}$ the at-risk indicator, the log partial likelihood can be written

$$\ell(\beta) = \sum_{i=1}^{n} d_i \left[ \beta^\top x_i - \log \sum_{j=1}^{n} Y_j(\tilde T_i) \exp(\beta^\top x_j) \right].$$

where:

• $n$: the number of individuals in the study
• $d_i \in \{0,1\}$: the event indicator for individual $i$ (1 if the event was observed, 0 if censored)
• $Y_j(t) = \mathbb{1}\{\tilde T_j \ge t\}$: the at-risk indicator for individual $j$, equal to 1 if $j$ is still at risk at time $t$
• $\tilde T_j$: the observed event or censoring time of individual $j$
• $x_i, x_j$: the covariate vectors of individuals $i$ and $j$
• $\beta$: the coefficient vector to be estimated

This form makes it clear that only individuals who have $d_i = 1$ contribute numerator terms, while censored individuals still appear in the denominator of every risk set that includes them. The risk-set denominator is the mechanism that lets censored observations speak: they do not need an event to influence the fit, because their presence tells us which covariate profiles were available to fail.

To maximize $\ell(\beta)$ we take derivatives. For the $p$-th component of $\beta$, differentiate the numerator $\beta^\top x_i$ to get $x_{i,p}$, and for the log-sum-exp use the chain rule: the derivative of $\log \sum_j Y_j e^{\beta^\top x_j}$ is the weighted average of $x_{j,p}$ using the exponential weights. The result is the score function

$$\frac{\partial \ell}{\partial \beta_p} = \sum_{k=1}^{K} \left[ x_{i_k,p} - \frac{\sum_{j \in R(t_{(k)})} x_{j,p} \exp(\beta^\top x_j)}{\sum_{j \in R(t_{(k)})} \exp(\beta^\top x_j)} \right].$$

Each event's contribution is exactly "the observed covariate of the failing individual minus the risk-set average of that covariate." Setting the score to zero is what fitting means in the Cox model, and it is a clean "observed minus expected" condition that mirrors ordinary MLE. An event contributes positively when the failing individual has a larger covariate value than the weighted risk-set mean. A positive total score component means that an infinitesimal increase in that coefficient, holding the others fixed, raises the log partial likelihood. A multivariable Newton update couples coordinates through the Hessian, so its individual coefficient changes need not follow the signs of the score components.

The Hessian is obtained by differentiating once more. Writing $\bar x(t_{(k)}) = \sum_{j \in R(t_{(k)})} w_j x_j / \sum_j w_j$ for the risk-set weighted mean with weights $w_j = \exp(\beta^\top x_j)$, we get

$$\frac{\partial^2 \ell}{\partial \beta_p \partial \beta_q} = -\sum_{k=1}^{K} \left[ \overline{(x_p x_q)}(t_{(k)}) - \bar x_p(t_{(k)}) \, \bar x_q(t_{(k)}) \right].$$

where:

• $\ell(\beta)$: the log partial likelihood
• $\beta_p, \beta_q$: the $p$-th and $q$-th components of the coefficient vector
• $k$: an index running over the $K$ distinct event times
• $t_{(k)}$: the $k$-th distinct event time
• $R(t_{(k)})$: the risk set at event time $t_{(k)}$
• $w_j = \exp(\beta^\top x_j)$: the weight of individual $j$ in the risk set
• $\bar x_p(t_{(k)}) = \sum_{j \in R(t_{(k)})} w_j x_{j,p} / \sum_j w_j$: the risk-set weighted mean of covariate $p$
• $\overline{(x_p x_q)}(t_{(k)})$: the risk-set weighted mean of the product $x_{j,p} x_{j,q}$

The bracket is a weighted covariance matrix and is therefore positive semidefinite, which makes the Hessian negative semidefinite. That is, $\ell(\beta)$ is concave. Concavity rules out a stationary point with a lower objective than an attainable global maximum. Strict concavity can give uniqueness when a finite maximum exists, but full design rank alone does not establish finite estimability. Separation or lack of risk-set variation can prevent a finite unique estimate, and unmodified Newton steps have no arbitrary-start convergence guarantee. Unlike ordinary least squares, this fit generally requires iterative optimization. Concavity does not establish an OLS-style closed-form solution, equal computational cost, or convergence of arbitrary Newton steps.

Real data often produce ties, in which several individuals record an event at the same time. With $d_k$ events at time $t_k$ and event set $D_k$, the untied expression no longer specifies an ordering within the tied group. Exact tie treatments depend on whether times are modeled as discrete or as coarsened continuous observations and can be expensive. Two approximations are common. The Breslow approximation treats each tied event as if the others were still at risk:

$$L_k^{\text{Breslow}}(\beta) = \frac{\prod_{i \in D_k} \exp(\beta^\top x_i)} {\left(\sum_{j \in R(t_k)} \exp(\beta^\top x_j)\right)^{d_k}}.$$

where:

• $k$: an index over the distinct event times $t_k$ that carry ties
• $t_k$: a distinct event time at which $d_k$ events are recorded
• $D_k$: the set of individuals who fail at time $t_k$, with $|D_k| = d_k$
• $d_k$: the number of tied events at time $t_k$
• $R(t_k)$: the risk set at time $t_k$
• $x_i, x_j$: the covariate vectors of individuals $i$ and $j$
• $\beta$: the coefficient vector

The Efron approximation accounts for depletion of the tied event set because it gradually removes events from the denominator as we imagine the tie broken:

$$L_k^{\text{Efron}}(\beta) = \frac{\prod_{i \in D_k} \exp(\beta^\top x_i)} {\prod_{\ell=0}^{d_k - 1} \left[ \sum_{j \in R(t_k)} \exp(\beta^\top x_j) - \frac{\ell}{d_k} \sum_{j \in D_k} \exp(\beta^\top x_j) \right]}.$$

where:

• $k$: an index over the distinct event times $t_k$ that carry ties
• $t_k$: a distinct event time at which $d_k$ events are recorded
• $D_k$: the set of individuals who fail at time $t_k$, with $|D_k| = d_k$
• $d_k$: the number of tied events at time $t_k$
• $R(t_k)$: the risk set at time $t_k$
• $\ell$: an index running from $0$ to $d_k - 1$ that progressively removes tied failures from the denominator
• $x_i, x_j$: the covariate vectors of individuals $i$ and $j$
• $\beta$: the coefficient vector

With one event at a time the two formulas coincide. With multiple tied events they can differ, and the difference depends on the event weights as well as the tie count. Efron is often used for coarsened continuous times; genuinely discrete-time data can call for a different model. We specify the tie method explicitly rather than relying on a software default.

The final ingredient is the baseline hazard. Cox's partial likelihood deliberately hides $h_0(t)$, but to turn coefficients into survival curves we need to recover it. The standard estimator is Breslow's:

$$\hat H_0(t) = \sum_{t_k \le t} \frac{d_k}{\sum_{j \in R(t_k)} \exp(\hat\beta^\top x_j)}.$$

Each event time contributes the number of events $d_k$ divided by the total at-risk weight at that time. Summing these increments accumulates an estimate of the cumulative baseline hazard. Given $\hat H_0(t)$, individual survival follows from

$$\hat S(t \mid x) = \exp\Bigl(-\hat H_0(t) \exp(\hat\beta^\top x)\Bigr).$$

where:

• $t$: time since the start of follow-up
• $x$: the covariate vector for the individual
• $\hat H_0(t)$: the Breslow estimate of the baseline cumulative hazard up to time $t$
• $\hat\beta$: the fitted coefficient vector
• $\hat S(t \mid x)$: the estimated survival probability for an individual with covariates $x$

Because $\hat H_0$ is a step function, $\hat S$ is a step function too; between event times, survival is treated as constant.

Three properties guide the fit. First, the log partial likelihood is concave, but a finite identifiable optimum and adequate risk-set covariate variation still matter; full column rank alone does not guarantee a finite estimate or global Newton convergence. Second, the time-invariant hazard multiplier applies to fixed covariates under proportional hazards. Third, inverse observed information gives an approximate model-based covariance under regularity conditions, rather than an exact finite-sample uncertainty guarantee.

There is also an important invariance we will use implicitly later. Because the partial likelihood depends on covariate values only through differences $x_i - x_j$ within each risk set, adding a constant shift to any covariate in $X$ has no effect on $\hat\beta$, only on the baseline, which absorbs the shift. A constant intercept is not separately identified from the Cox baseline hazard. This differs from ordinary logistic regression, where an intercept usually is identifiable.