Notes 3 — Comparing Survival Curves, Continued
Module 1 · Comparing survival curves, continued
Lecture 3.1 onward. Reading: Lecture Notes Part 1, §1.5.
Lecture 3.1 — Comparing survival curves, continued
Multigroup comparison
Now suppose there are $J$ groups rather than two, and we want to assess differences across all of them at once. The hypothesis generalizes the two-group log-rank hypothesis:
\[H_0 : S_1(t) = \cdots = S_J(t) \;\;\text{for all } t > 0\] \[\text{vs}\] \[H_A : S_j(t) \ne S_{j'}(t) \;\;\text{for some } t > 0, \;\;\text{for at least some } j \ne j'\]Read the quantifiers carefully — they stack:
- The null demands all $J$ curves coincide at every time point.
- The alternative is satisfied by any single pair of groups $j \ne j’$ differing at any single time $t$.
So rejecting $H_0$ says only that the curves are not all identical — it does not say which groups differ, or where. This is exactly the ANOVA situation: an omnibus test, with any pairwise follow-up a separate question.
The J×2 table
As in the two-group case, pool the data — now across all $J$ groups — and order the unique event times $t_{(1)} < \cdots < t_{(K)}$. At the $k$-th event time the 2×2 table simply grows a row per group, becoming a $J \times 2$ table:
| No event | Event | ||
|---|---|---|---|
| Group 1 | $n_{1k} - d_{1k}$ | $d_{1k}$ | $n_{1k}$ |
| Group 2 | $n_{2k} - d_{2k}$ | $d_{2k}$ | $n_{2k}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
| Group J | $n_{Jk} - d_{Jk}$ | $d_{Jk}$ | $n_{Jk}$ |
| $n_k - d_k$ | $d_k$ | $n_k$ |
with, for group $j \in {1, \dots, J}$:
\[\begin{aligned} n_{jk} &= \text{\# of individuals at risk in the } j\text{th group at time } t_{(k)} \\[4pt] d_{jk} &= \text{\# of events in the } j\text{th group at time } t_{(k)} \end{aligned}\]The margins pool over the rows exactly as before: $n_k = \sum_j n_{jk}$ at risk in total, $d_k = \sum_j d_{jk}$ events in total at $t_{(k)}$. Again one table per event time, $K$ tables in all.
Observed and expected vectors
With two groups, one number per table ($U_k$) carried all the information. With $J$ groups the per-table summary becomes a vector — the observed event counts in groups $1$ through $J-1$:
\[\mathbf{O}_k = \big( d_{1k},\; d_{2k},\; \dots,\; d_{(J-1)k} \big)\]Why only $J - 1$ entries? The margins of the table are fixed: the $d_{jk}$ must sum to $d_k$, so the last group’s count is determined by the others — $d_{Jk} = d_k - \sum_{j<J} d_{jk}$. Including it would add no information (and would make the variance matrix singular). Same reason the two-group test tracked only Group 1.
Under the null, conditional on the margins, the expected counts are each group’s proportional share of the $d_k$ events — the multigroup version of $E(D_{1k}) = n_{1k} d_k / n_k$:
\[\mathbf{E}_k = \left( \frac{n_{1k}\, d_k}{n_k},\; \frac{n_{2k}\, d_k}{n_k},\; \dots,\; \frac{n_{(J-1)k}\, d_k}{n_k} \right)\]And the per-table variance $V_k$ becomes the $(J-1) \times (J-1)$ variance–covariance matrix of $\mathbf{O}_k - \mathbf{E}_k$ — covariances appear because the group counts within one table compete for the same $d_k$ events: one group getting more forces the others to get fewer.
The multigroup log-rank statistic
Pool each piece across the $K$ event times:
\[\mathbf{O} = \sum_k \mathbf{O}_k \qquad \mathbf{E} = \sum_k \mathbf{E}_k \qquad \mathbf{V} = \sum_k \mathbf{V}_k\]and combine them in a quadratic form:
\[T_{LR} = (\mathbf{O} - \mathbf{E})^T\, \mathbf{V}^{-1}\, (\mathbf{O} - \mathbf{E}) \;\sim\; \chi^2_{J-1}\]Under $H_0$ of no differences across the $J$ groups, $T_{LR}$ has a chi-square distribution with $J - 1$ degrees of freedom — one per free group.
This is the same statistic as before, grown up. With $J = 2$ the vectors collapse to scalars: $\mathbf{O} - \mathbf{E} = \sum_k U_k$, $\mathbf{V} = \sum_k V_k$, and the quadratic form reduces to $(\sum_k U_k)^2 / \sum_k V_k \sim \chi^21$ — exactly the two-group $T{LR}$. The quadratic form $(\cdot)^T \mathbf{V}^{-1} (\cdot)$ is the matrix generalization of “square the deviation and divide by its variance.”
Stratified comparison
Back to comparing two treatment groups — but now adjusting for another categorical variable, with the comparison carried out within its levels.
Why stratify:
- Observational studies — stratification is one approach to control for potential confounding bias.
- Randomized trials — randomization is often stratified to ensure the treatment allocation is balanced across levels of an important prognostic factor (e.g. “randomization was stratified according to site of metastasis, age, and estrogen-receptor status”).
Suppose the stratifying variable has $M$ strata. Running example — the smoking study stratified by gender:
\[\text{gender} \in \{M, F\} \qquad\Longrightarrow\qquad M = 2\]Let $S_{1,m}(t)$ and $S_{0,m}(t)$ denote the survivor functions of treatment groups 1 and 0 within stratum $m$. The null:
\[H_0 : S_{1,m}(t) = S_{0,m}(t) \;\;\text{for all } m, \;\;\text{for all } t > 0\]That is: within every stratum, the two treatment curves coincide everywhere. The treatment comparison happens within strata — men are compared to men, women to women — so a difference between the strata themselves (women faring differently from men overall) does not count as evidence against $H_0$. Only a treatment effect inside a stratum does.
The stratified log-rank statistic
The machinery is the two-group log-rank test, run per stratum. At the $k$-th event time, build the 2×2 table within stratum $m$ — every count picks up an extra subscript:
\[\begin{aligned} n_{1k,m} &= \text{\# at risk in treatment group 1 at } t_{(k)} \text{ from stratum } m \\[4pt] d_{1k,m} &= \text{\# of events in treatment group 1 at } t_{(k)} \text{ from stratum } m \end{aligned}\]and likewise $n_{0k,m}$, $d_{0k,m}$, with within-stratum margins $n_{k,m}$ and $d_{k,m}$.
Each stratum’s table gets its own hypergeometric null — expectation and deviation computed against that stratum’s margins:
\[E(D_{1k,m}) = \frac{n_{1k,m}\, d_{k,m}}{n_{k,m}}\] \[U_{k,m} = d_{1k,m} - \frac{n_{1k,m}\, d_{k,m}}{n_{k,m}} \qquad\qquad \text{Var}[U_{k,m}] = V_{k,m}\]This is where the adjustment lives: Group 1’s “fair share” of the events at $t_{(k)}$ is computed from its share of the risk set within the stratum, so a stratum with worse overall survival cannot masquerade as a treatment effect.
Now pool over both event times and strata — sum everything, then standardize once:
\[T_S = \frac{\left[ \sum_{m=1}^{M} \sum_{k=1}^{K} U_{k,m} \right]^2} {\sum_{m=1}^{M} \sum_{k=1}^{K} V_{k,m}} \;\sim\; \chi^2_1\]Still $\chi^2_1$ — not $\chi^2_M$. The strata are not being compared; they each contribute their own observed-minus-expected evidence about the same single question (does treatment 1 differ from treatment 0?), and that evidence is summed into one deviation before squaring. One comparison, one degree of freedom.
Compare with the unstratified $T_{LR}$: the only change is that expectations and variances are computed within strata and the sum gains an $m$ index. If $M = 1$ the two statistics are identical.
When to use it. The stratified log-rank test has good power when the direction of the treatment effect is the same within each stratum — implicitly, it estimates a common effect across strata. If the interest is in stratum-specific effects (interaction / effect modification — treatment helping women but not men), the stratified log-rank test is the wrong tool: summing $U_{k,m}$ across strata would let opposite-signed effects cancel.
Weighted log-rank tests
The log-rank statistic gives equal weight to every risk set. Sometimes different time points deserve different emphasis. Let $w_k$ be a weight attached to the risk set at $t_{(k)}$; the weighted log-rank statistic is
\[T_W = \frac{\left[ \sum_{k=1}^{K} w_k\, U_k \right]^2} {\sum_{k=1}^{K} w_k^2\, V_k} \;\sim\; \chi^2_1\]Note where the square lands in the denominator: scaling a random variable by $w_k$ scales its variance by $w_k^2$, so the weights enter the numerator linearly but the denominator squared — that keeps the standardization honest. With $w_k = 1$ for all $k$, $T_W$ is exactly $T_{LR}$.
Strategy. Choose weights that emphasize the parts of the time scale where differences are expected, and de-emphasize the parts where there are none.
Pre-specify the weights. Chosen after seeing the data, weights are a fishing license — pick the $w_k$ that inflate your effect. For randomized trials the weights typically must be pre-registered in the analysis plan.
Proposed tests differ only in their choice of $w_k$:
| Test statistic | Weight $w_k$ |
|---|---|
| Log-rank | $1$ |
| Gehan–Breslow | $n_k$ |
| Peto/Prentice | $n \hat{S}(t_{(k)})$ |
| Generalized Wilcoxon | $\hat{S}(t_{(k)})$ |
| Fleming–Harrington | $[\hat{S}(t_{(k)})]^{\rho}\, [1 - \hat{S}(t_{(k)})]^{\gamma}$ |
E.g. Gehan–Breslow weights each risk set by its size $n_k$ — since risk sets shrink over time, that up-weights early differences (and counters the plain log-rank test’s sensitivity to noisy tails). The Fleming–Harrington family tunes the emphasis via $(\rho, \gamma)$: $\rho > 0, \gamma = 0$ leans early, $\rho = 0, \gamma > 0$ leans late, $\rho = \gamma = 0$ recovers the log-rank test.
Modeling survival data: regression
Recap of the observed data — each subject contributes the triple:
\[(X_i,\; \delta_i,\; \mathbf{Z}_i)\] \[\begin{aligned} X_i &= \text{censored event time} \\[4pt] \delta_i &= \text{event indicator} \\[4pt] \mathbf{Z}_i &= \text{set of covariates or characteristics} \end{aligned}\]When $\mathbf{Z}_i$ is a single binary or categorical variable, the log-rank machinery above suffices. But when $\mathbf{Z}_i$ is continuous, or holds multiple variables, group-wise comparison breaks down — there are no groups to compare, or too many. Instead we need to model the distribution of $T \mid \mathbf{Z}$.
Why not linear regression?
Outside survival analysis, the reflex would be ordinary linear regression on the event time:
\[E(T \mid Z) = \beta_0 + \beta_1 Z\]Three reasons this fails for time-to-event outcomes:
- $T > 0$ — event times are non-negative, but nothing stops a linear model from predicting a negative survival time.
- $T$ tends to be right-skewed — linear regression’s machinery (and its inference) leans on symmetric, roughly normal errors; survival times pile up early with a long right tail.
- The mean is not a good measure of centrality — with skewed distributions (and censoring hiding the right tail) the mean is a poor summary; survival work prefers the median and the curve itself.
(And beyond these: linear regression has no way to use a censored observation — a subject known only to have $T_i > x_i$ fits nowhere in least squares.)
Parametric regression
Extend the parametric estimation from Lecture 1.3 to incorporate $\mathbf{Z}$. Assume exponential survival — but now give each subject their own rate:
\[T_i \sim \text{Exponential}(\lambda_i)\]What used to be a single population parameter $\lambda$ is now specific to each study unit: $\lambda_i$. The covariates enter through the rate — but the exponential parameter must be $> 0$, and a linear combination $\beta^T \mathbf{Z}_i$ can be any real number. The fix is to model the log of the rate:
\[\log(\lambda_i) = \beta^T \mathbf{Z}_i \qquad\Longleftrightarrow\qquad \lambda_i = \exp\big( \beta^T \mathbf{Z}_i \big)\]where $\beta$ is a vector of regression coefficients and
\[\beta^T \mathbf{Z}_i = \beta_1 Z_{i1} + \beta_2 Z_{i2} + \cdots + \beta_p Z_{ip}\]is the familiar linear predictor. Whatever value it takes, $\exp(\cdot)$ maps it to a positive rate — the same trick as the log-normal’s $T = \exp(N) > 0$, applied to a parameter instead of a time.
Reading the signs (recall $\mu_i = 1/\lambda_i$ for the exponential): a positive $\beta_j$ means $Z_j$ raises the hazard, hence shortens expected survival.
This is parametric regression under an exponential assumption for $T$ — and it can be rewritten as an accelerated failure time (AFT) model.
Accelerated failure time (AFT) models
The general representation of an AFT model puts the regression on the time scale rather than the rate scale:
\[\log(T_i) = \beta^T \mathbf{Z}_i + \epsilon_i\]coupled with some specification of the distribution of the error $\epsilon_i$. It looks exactly like linear regression — but on $\log T$, which resolves the earlier objections: $T_i = \exp(\beta^T \mathbf{Z}_i + \epsilon_i)$ is automatically positive, and the log tames the right skew.
The choice of the $\epsilon$ distribution is the important one. It:
- defines the distribution of the time-to-event outcome $T_i$ — e.g. extreme-value errors give exponential/Weibull $T$ (recall the lecture-1 note that the log of an exponential is extreme-value), normal errors give log-normal $T$;
- dictates how flexibly the AFT model can represent complex data.
Why “accelerated failure time”? Exponentiating, the covariates multiply time itself: $T_i = e^{\beta^T \mathbf{Z}_i} \cdot e^{\epsilon_i}$ — a subject’s clock runs faster or slower by the factor $e^{\beta^T \mathbf{Z}_i}$. Here a positive coefficient lengthens survival — note the sign flip relative to the rate-scale model $\log \lambda_i = \beta^T \mathbf{Z}_i$, where a positive coefficient raised the hazard and shortened survival. The two parametrizations describe the same exponential model with $\beta$’s of opposite sign — always check which scale software reports.
The Cox proportional hazards model
The alternative to parametric regression — and the most common model used for survival data. Instead of modeling $T$ (or $\log T$) directly, the Cox proportional hazards model specifies the hazard function:
\[\lambda(t \mid \mathbf{Z}) = \lambda_0(t)\, \exp\big( \beta^T \mathbf{Z} \big)\]The pieces:
| Term | Meaning |
|---|---|
| $\lambda(t \mid \mathbf{Z})$ | hazard function at time $t$, given covariates $\mathbf{Z}$ |
| $\lambda_0(t)$ | baseline hazard function — the hazard for the population with $\mathbf{Z} = \mathbf{0}$ |
| $\exp(\beta^T \mathbf{Z})$ | multiplicative effect of the covariates on the hazard |
Setting $\mathbf{Z} = \mathbf{0}$ makes the exponential factor $e^0 = 1$, leaving $\lambda(t \mid \mathbf{0}) = \lambda_0(t)$ — that is why it is the baseline. Every subject’s hazard is their covariates’ constant multiplier applied to the same shared baseline curve.
Two structural points:
- Modeling the hazard is a less direct choice — the hazard is not as straightforward to interpret as a mean or a survival time; the payoff comes in the interpretation of $\beta$ (next).
- $\lambda_0(t)$ is left unspecified — no assumption at all about its shape. This is the model’s defining move: the covariate effect $\exp(\beta^T \mathbf{Z})$ is parametric, the baseline hazard is nonparametric, making the Cox model semiparametric — “half” parametric, “half” nonparametric. It buys the flexibility of the Kaplan–Meier world (no distributional assumption on time) while still delivering regression coefficients for arbitrary covariates.
Interpreting $\beta$: the hazard ratio
To read the $j$-th coefficient, compare two covariate vectors that differ only in coordinate $j$, by exactly one unit:
\[\mathbf{z} = (z_1, \dots, z_j, \dots, z_p) \qquad \mathbf{z}' = (z_1, \dots, z_j + 1, \dots, z_p)\]Take the ratio of the two hazards the model assigns:
\[\begin{aligned} \frac{\lambda(t \mid \mathbf{z}')}{\lambda(t \mid \mathbf{z})} &= \frac{\cancel{\lambda_0(t)}\, \exp\{\beta^T \mathbf{z}'\}} {\cancel{\lambda_0(t)}\, \exp\{\beta^T \mathbf{z}\}} \\[6pt] &= \frac{\exp\{\beta_1 z_1 + \beta_2 z_2 + \cdots + \beta_j (z_j + 1) + \cdots\}} {\exp\{\beta_1 z_1 + \beta_2 z_2 + \cdots + \beta_j z_j + \cdots\}} \\[6pt] &= \frac{e^{\beta_1 z_1}\, e^{\beta_2 z_2} \cdots e^{\beta_j z_j}\, e^{\beta_j} \cdots} {e^{\beta_1 z_1}\, e^{\beta_2 z_2} \cdots e^{\beta_j z_j} \cdots} \\[6pt] &= e^{\beta_j} \end{aligned}\]Two cancellations do all the work:
- The baseline $\lambda_0(t)$ cancels — this is the “proportional hazards” property. The ratio does not depend on $t$: whatever the baseline does over time, the two hazards stay in the same constant proportion. It is also why $\beta$ can be interpreted without ever knowing $\lambda_0$.
- Every matching term cancels — the exponential of a sum factors into a product ($e^{a+b} = e^a e^b$), and all factors are identical except the extra $e^{\beta_j}$ from the $+1$.
So:
\[e^{\beta_j} = \text{the \textbf{hazard ratio} (HR)}\]— the relative change in the hazard at time $t$ associated with a one-unit increase in $Z_j$, holding everything else constant. Note it is the relative change: $e^{\beta_j} = 1.5$ means 50% higher hazard, whatever the baseline level.
Reading the sign (assuming the outcome is “bad” — death, relapse):
| Coefficient | Hazard ratio | Impact of $Z_j$ | |
|---|---|---|---|
| $\beta_j < 0$ | $e^{\beta_j} < 1$ | decreases the magnitude of the hazard | $\Rightarrow$ benefit |
| $\beta_j = 0$ | $e^{\beta_j} = 1$ | none | no effect |
| $\beta_j > 0$ | $e^{\beta_j} > 1$ | increases the magnitude of the hazard | $\Rightarrow$ harm |
The mapping works because $e^x$ is increasing with $e^0 = 1$: the sign of $\beta_j$ lands the hazard ratio on one side of $1$ or the other. Note this is the same direction convention as the rate-scale exponential regression earlier (higher hazard = worse), and the opposite of the AFT time-scale convention.
Why “proportional hazards”
Look again at what the derivation produced:
\[\underbrace{\frac{\lambda(t \mid \mathbf{Z} = \mathbf{z}')}{\lambda(t \mid \mathbf{Z} = \mathbf{z})}}_{\text{a function of } t\text{?}} = \underbrace{e^{\beta_j}}_{\text{no } t \text{ involved}}\]The left side is a ratio of two functions of $t$; the right side is a constant. Regardless of the value of $t$, the two hazards stay in the same proportion — that is the proportional hazards property, and it is an assumption baked into the model’s form. Its implications:
- one covariate’s impact is summarized by a single number — no need for a different hazard ratio at each time;
- time never needs to be invoked when interpreting a study’s results.
Whether proportional hazards actually holds for a given study is a real question (a treatment that helps early but not late violates it) — to be revisited.
Fitting the Cox model: the partial likelihood idea
To estimate $\beta$, Cox (1972) proposed the partial likelihood — different from the usual likelihood specification, but built from two ideas already familiar from Kaplan–Meier and the log-rank test:
- Contributions to the likelihood are made only by people who have the event — nothing is contributed at censoring times, or at times between events.
- Consider the risk set — at each event time, who is still at risk?
The full likelihood would need $\lambda_0(t)$ — which the model deliberately left unspecified. The partial likelihood sidesteps it by asking, at each event time, a question whose answer doesn’t involve the baseline: given that someone in the risk set has the event right now, what is the probability it is this particular person?
The partial likelihood contribution
At a particular event time $t_j$, let $\mathbf{Z}_j$ be the covariates of the individual who had the event there, and $\mathcal{R}_j$ the risk set. The contribution to the partial likelihood is that individual’s hazard as a share of the total hazard among everyone at risk:
\[\begin{aligned} L_j(\beta) &= \frac{\lambda(t \mid \mathbf{Z}_j)} {\sum_{l \in \mathcal{R}_j} \lambda(t \mid \mathbf{Z}_l)} \\[8pt] &= \frac{\cancel{\lambda_0(t_j)}\, \exp(\beta^T \mathbf{Z}_j)} {\sum_{l \in \mathcal{R}_j} \cancel{\lambda_0(t_j)}\, \exp(\beta^T \mathbf{Z}_l)} \\[8pt] &= \frac{\exp(\beta^T \mathbf{Z}_j)} {\sum_{l \in \mathcal{R}_j} \exp(\beta^T \mathbf{Z}_l)} \end{aligned}\]The middle step substitutes the Cox model for every hazard. The baseline $\lambda_0(t_j)$ appears in the numerator and in every term of the denominator’s sum — it is a common factor at the fixed time $t_j$, so it factors out of the sum and cancels completely. What remains depends only on $\beta$ and on who was at risk — the unspecified baseline has vanished, which is exactly what makes estimation possible without ever modeling it.
Read $L_j(\beta)$ as a probability: of all the people at risk at $t_j$, each carries relative hazard weight $\exp(\beta^T \mathbf{Z}_l)$, and the observed event went to individual $j$ with probability proportional to their weight. A $\beta$ that gives high weights to the people who actually had the events will make these fractions large — that is the $\beta$ the data favors.