Cox Proportional Hazards

Kaplan-Meier and Nelson-Aalen summarize survival and cumulative hazard for a population or stratum. They can describe a mixture of people with different hazards when their censoring assumptions hold; they do not themselves specify how predictors change an individual's hazard. For equipment with different maintenance histories, a regression model lets us ask how the hazard is associated with those histories while retaining censored observations in the appropriate risk sets.

The Cox proportional hazards model, introduced by David Cox in 1972, is the standard regression model for censored outcomes. Its defining feature is simple. Each individual's hazard is written as a shared, unspecified baseline hazard $h_0(t)$ multiplied by an individual-specific factor $\exp(\beta^\top x)$ that depends only on that individual's covariates. The covariate factor scales the hazard up or down by the same amount at every point in time, which is precisely the proportional-hazards assumption. Because the baseline hazard is left entirely unspecified, Cox's model is semiparametric: it makes a strong parametric claim about how covariates act, and no parametric claim about the shape of time.

The method's second insight is what makes it practical. Cox showed that $\beta$ can be estimated by maximizing a partial likelihood that uses only the order in which events occur and who was still at risk when they did. The unknown baseline $h_0(t)$ cancels at every event time, so we never have to specify it. That cancellation is the core idea of the chapter. We will derive it from first principles, build the fitting machinery from scratch with NumPy, and reconstruct the baseline and survival curves afterwards. Along the way, you will see how each formula translates into code and how to check the result against a mature library.

A hazard ratio is not a survival-probability ratio, and it is not an effect size on the probability scale. Even when hazards are strictly proportional, the ratio of survival probabilities at time $t$ depends on $t$ and on the baseline hazard itself. Holding this distinction straight is essential for correctly communicating results.

The Cox model leaves the baseline hazard shape unspecified. A parametric survival model instead chooses a family for that shape, which can add information if the choice is appropriate but can introduce misspecification. The Cox model trades that distributional choice for the proportional-hazards assumption and a model for the covariate effects. Its efficiency relative to a parametric fit depends on the data-generating process and censoring; there is no universal few-percent information-loss bound.

The hazard-ratio interpretation also suits many applied questions well. Clinicians want to know how much a treatment multiplies a patient's hazard at any moment; engineers want to know how much a maintenance schedule multiplies a failure rate; product managers want to know how much a pricing change multiplies a cancellation hazard. In each case, the Cox coefficient reports the log hazard ratio per unit of covariate, and its exponentiated form $\exp(\beta)$ reports the multiplicative change in hazard.

Statistical inference for the Cox model is well established. The partial likelihood behaves asymptotically like an ordinary likelihood (Andersen and Gill, 1982), so model-based standard errors, Wald intervals, likelihood-ratio tests and score tests are available under their regularity and sampling assumptions. Clustering, misspecification and sparse events can require additional methods. Time-varying covariates, stratification, and penalization all extend naturally. These tools support a familiar regression workflow, while censoring-aware risk sets, proportional-hazards diagnostics, and survival-specific evaluation remain necessary.

The most serious disadvantage is the proportional-hazards assumption itself. The entire model depends on the claim that $\exp(\beta^\top x)$ multiplies the hazard by the same amount at every point in time. If that claim fails, for example, if a treatment is protective early and harmful later, a fitted coefficient can become a population- and follow-up-dependent summary, but their interpretation as a constant hazard ratio becomes misleading, and the survival predictions that follow from the model will be systematically off. We will discuss the standard diagnostics in detail later, because in practice it is the assumption that fails, not the mathematics.

Its relative-hazard interpretation must also be separated from absolute survival and from causal effects. An exponentiated Cox coefficient gives the multiplier of the instantaneous hazard for a one-unit covariate increase, holding the other covariates fixed. The hazard is a rate, not a probability. It says nothing directly about how much survival probability changes, because that depends on the cumulative baseline hazard and on how far into the future you are looking. Reporting "a hazard ratio of 0.5" as if it were "half as likely to die" is a common and consequential misinterpretation.

Finally, the model is not free of practical limitations. The simple model assumes noninformative censoring conditional on the included covariates. If dropout still depends on prognosis after conditioning on those covariates, the fitted coefficients can be biased. It assumes the log hazard is linear in the included covariates; nonlinear effects require their own functional-form checks. Our treatment of tied event times uses approximations. Efron accounts for fractional depletion of the tied-event weight in coarsened continuous times. The recording process and risk-set composition also matter when choosing a tie method. It also does not give you an absolute risk prediction unless you re-estimate the baseline hazard, which is an extra step you may be tempted to skip.

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.

We work through a minimal dataset of four subjects so every number can be checked by hand. Each subject has a single covariate $x \in \{0, 1\}$, an observed time, and an event indicator.

There are three events, at times $2$, $4$, and $6$. Subject 4 is censored at time 8, so they are at risk at every event time but never contribute a numerator term. The three risk sets are:

• at $t=2$: $\{1, 2, 3, 4\}$ (everyone is still at risk)
• at $t=4$: $\{2, 3, 4\}$ (subject 1 has already failed)
• at $t=6$: $\{3, 4\}$ (subject 2 has already failed)

The partial likelihood is

$$L(\beta) = \underbrace{\frac{e^{0}}{e^{0} + e^{\beta} + e^{0} + e^{\beta}}}_{t=2} \cdot \underbrace{\frac{e^{\beta}}{e^{\beta} + e^{0} + e^{\beta}}}_{t=4} \cdot \underbrace{\frac{e^{0}}{e^{0} + e^{\beta}}}_{t=6}.$$

Writing $q = e^{\beta}$ and simplifying,

$$L(\beta) = \frac{1}{2 + 2q} \cdot \frac{q}{1 + 2q} \cdot \frac{1}{1 + q}.$$

Taking logs,

$$\ell(\beta) = -\log(2 + 2q) + \beta - \log(1 + 2q) - \log(1 + q).$$

Let us evaluate this at three values of $\beta$ to see the shape of the surface. At $\beta = 0$ ($q=1$):

$$\ell(0) = -\log 4 + 0 - \log 3 - \log 2 \approx -1.386 - 1.099 - 0.693 = -3.178.$$

At $\beta = -1$ ($q \approx 0.368$):

$$\ell(-1) = -\log(2.736) - 1 - \log(1.736) - \log(1.368) \approx -1.006 - 1 - 0.552 - 0.313 = -2.871.$$

At $\beta = -1.5$ ($q \approx 0.223$):

$$\ell(-1.5) = -\log(2.446) - 1.5 - \log(1.446) - \log(1.223) \approx -2.965.$$

The displayed evaluations bracket the maximum: $\ell(-1.5) \approx -2.965$, $\ell(-1) \approx -2.871$, and $\ell(0) \approx -3.178$, so the largest value is at $\beta=-1$ and the optimum lies between $-1.5$ and $0$. The score can be simplified to $1-2q/(1+q)-2q/(1+2q)$. At $\beta=-1$ it is positive, so increasing $\beta$ locally increases the log partial likelihood. Multiplying the zero-score equation $1-2q/(1+q)-2q/(1+2q)=0$ by its positive denominators yields $1-q-4q^2=0$, whose positive root $q=(\sqrt{17}-1)/8$ determines the finite optimum in this example. We compute it and an approximate Wald interval below; with only three events, the interval is a teaching calculation rather than a precise effect estimate.

The lesson from the hand calculation is not the specific coefficient. It is that the optimization operates on observed covariate values within risk sets. At an event where the failing individual has $x=1$, increasing $\beta$ raises its numerator weight; when the failing individual has $x=0$, that numerator stays at one. In either case, increasing $\beta$ raises the denominator weights of the subjects with $x=1$ who remain at risk. A positive $\beta$ helps only when the failing individual has high covariate values relative to those still at risk. Fitting the Cox model means finding the $\beta$ that balances those two forces across all events.

We now build the partial likelihood from scratch. We use a small equipment-failure dataset with a single binary covariate: maintained equals 1 if the equipment received regular maintenance and 0 otherwise. Our goal is to estimate the effect of maintenance on the failure hazard and then to recover the implied survival curves. The dataset is deliberately small so that the connection between the formulas and the code stays visible.

The hidden setup imports the numerical libraries and shared figure helpers. The table below is a small constructed teaching dataset, with follow-up measured in arbitrary time units; it is not a collected maintenance study.

The two coefficient columns and the two standard-error columns should line up to several decimal places in the displayed table; small residual differences are compatible with optimizer and information-matrix-solve tolerance. Larger discrepancies call for checking the inputs, objective conventions, derivatives, and numerical settings. This is the expected outcome: PHReg and the scratch implementation maximize the same concave log partial likelihood.

Our example has no tied event times, so the untied likelihood is exact. Real data often produce ties; exact methods and approximations depend on how those recorded times are modeled. We can demonstrate the difference by coarsening the observation times into four-unit bins, which forces multiple events to land at the same recording time. We implement both the Breslow and Efron approximations and compare them. The Breslow version treats all tied failures as if they occurred with the full risk set available; the Efron version fractionally reduces the total weight of the tied failures across successive denominators.

The executed tied-data table reports the Breslow and Efron coefficients and their hazard ratios alongside the explicitly selected statsmodels Efron fit. Our Efron coefficient matches the statsmodels Efron coefficient to within $10^{-5}$. These are different approximations applied to this coarsened dataset, not a universal statement about the direction or size of tie bias.

Estimating the baseline hazard lets us convert coefficients into absolute survival predictions rather than relative risk scores. We implement Breslow's estimator and predict survival for both groups at a grid of times. The estimator accumulates one increment at each event time: the number of events divided by the total at-risk weight. The prediction function then applies $\hat S(t \mid x) = \exp(-\hat H_0(t)e^{\hat\beta^\top x})$.