Survival Analysis: Method

Survival Analysis

Censored times

Survival data record, for each subject, the time from a starting point to an event (death, relapse, failure, churn) or, when the event has not happened by the end of observation, the time for which the subject was followed. The second kind of time is right-censored: it is known only that the event lies beyond it. Discarding censored subjects or treating their times as event times both bias the answer; every method on this page instead uses a censored subject for as long as it was observed, by counting it among those at risk up to its time and not after. The methods assume that censoring is uninformative: a subject censored at t has the same future prospects as the subjects still under observation at t. The status column says which kind each time is, and the Event Code names the value that marks an event.

Kaplan–Meier estimate

Let the distinct event times be t₁ < t₂ < …, with d_j events at t_j among the n_j subjects at risk just before it (those whose time is t_j or more). The probability of surviving past t is the product of the conditional probabilities of surviving each event time up to t: S(t) = Π_{tj ≤ t} (1 − d_j/n_j). The curve is a step function that drops at event times only; a censored time lowers the later numbers at risk. Greenwood's formula gives the variance of ln S as Σ d_j/[n_j(n_j − d_j)], so the standard error of S is S times the square root of that sum; without censoring it reduces to the binomial S(1 − S)/n. Three confidence limits are offered: log, S exp(± z se_{ln S}) truncated at 1 (the default of R's survfit); log-log, exp{−exp[ln(−ln S) ± z se_{ln S}/ln S]}, which stays inside 0 to 1 by construction and behaves best in small samples; and plain, S ± z se, truncated at 0 and 1. When the curve reaches zero the standard error is undefined and the limits are left blank.

The median survival time is the first time at which S falls to 0.5 or below; if the curve is exactly 0.5 on a flat stretch, the midpoint of that stretch is taken. The quartiles are found in the same way at 0.75 and 0.25. Their confidence limits follow Brookmeyer and Crowley: the lower limit is the time at which the lower limit curve of S reaches the level, the upper limit the time at which the upper limit curve does; a limit is blank when its curve never gets there, which is usual for the upper limit of a median in a small or heavily censored group. The mean survival time is the area under S, but it is undefined when the last observation is censored, so the restricted mean (RMST) integrates only to a stated limit τ; its standard error sums the squared remaining areas times the Greenwood increments. The Nelson–Aalen estimate H(t) = Σ d_j/n_j of the cumulative hazard, with variance Σ d_j/n_j², is tabulated beside S; exp(−H) is close to S except where few subjects remain.

Comparing groups: the log-rank family

Under the hypothesis that the groups share one survival curve, the d_j events at t_j should fall on the groups in proportion to their numbers at risk, so group k expects E_kj = d_j n_kj/n_j of them. The statistic of group k is U_k = Σ w_j(d_kj − E_kj), the weighted sum of observed minus expected events, and the covariance of the U follows from the hypergeometric distribution of each event time's table: V_kl = Σ w_j² d_j(n_j − d_j)/(n_j − 1) [δ_kl n_kj/n_j − n_kjn_lj/n_j²]. The U sum to zero, so one group is dropped and U′V⁻¹U is compared with χ² on (number of groups − 1) degrees of freedom. With w = 1 this is the log-rank (Mantel–Haenszel) test, the most powerful when the hazards are proportional. The other weights shift the emphasis: Gehan–Breslow w = n_j and Tarone–Ware w = √n_j weight early times, where many are at risk (Gehan's weights depend on the censoring pattern and can mislead when the groups are censored differently); Peto–Peto uses a survival estimate S̃(t_j) = Π(1 − d/(n + 1)) instead; the Fleming–Harrington family uses S(t_j−)^ρ[1 − S(t_j−)]^γ with the pooled Kaplan–Meier estimate just before each time, so ρ = 1, γ = 0 is the Peto–Prentice form of the Wilcoxon test (R's survdiff with rho = 1) and ρ = 0, γ = 1 looks for late differences. The page reports the chosen test in detail and the whole family beside it; the weight should be chosen before the data are seen.

With a strata column the expected events and the covariance are computed within each stratum and summed, which removes the effect of the stratifying factor from the comparison. For three or more ordered groups the trend statistic (s′U)²/(s′Vs) with scores s = 1, 2, … tests a steady rise or fall of the hazard on one degree of freedom. For two groups the log-rank quantities also give quick hazard-ratio estimates: the ratio of the observed-to-expected ratios, and Peto's one-step value exp(U/V) with limits exp(U/V ± z/√V). Both shrink a large ratio toward 1; the Cox estimate is preferred.

Cox proportional-hazards model

The model writes the hazard of a subject with covariates x as h(t | x) = h₀(t) exp(β′x): an unspecified baseline hazard multiplied by a factor that does not depend on time. exp(β) is the hazard ratio per unit of a numeric covariate, or of a level of a categorical one against its reference (categorical columns enter as dummy variables). The coefficients maximise the partial likelihood, in which each event contributes the ratio of its subject's risk score r = exp(β′x) to the sum of the scores of everyone at risk at that time; the baseline hazard cancels. When d events share a time the partial likelihood must be approximated: Breslow's form uses the full risk-set sum for each of them, Efron's subtracts k/d of the tied events' scores from the sum for the k-th of them (k = 0 … d − 1) and is the more accurate; they coincide without ties. The page solves the score equations by Newton–Raphson with step halving, starting from β = 0, and takes the standard errors from the inverse of the information matrix (minus the second derivative) at the solution.

Three tests of the hypothesis that all coefficients are zero are reported: the likelihood-ratio statistic 2(ℓ₁ − ℓ₀), the Wald statistic β′Iβ and the score statistic U′I⁻¹U evaluated at β = 0, each χ² with as many degrees of freedom as coefficients. The score test of a model holding only the group dummies is the log-rank test when no event times are tied. If a coefficient runs off to infinity (a group with no events, or a covariate that orders the event times perfectly) the likelihood has no finite maximum; the page says so, and the likelihood-ratio and score tests remain usable. With a strata column each stratum has its own baseline hazard and risk sets, and the coefficients are common. Subjects with a missing covariate are left out of the model only.

Harrell's concordance C is the fraction of usable pairs in which the subject with the higher risk score had the earlier event. A pair is usable when the shorter time is an event; pairs tied in the score count one half, pairs with two events at the same time are set aside, and in a stratified model only pairs within a stratum are compared. Its standard error is the infinitesimal-jackknife value of R's concordance function. The baseline cumulative hazard H₀(t) is estimated in the same tie form as the fit (Breslow: d_j/Σr at each event time), and the predicted survival of a subject with covariates x is exp[−H₀(t) e^β′x]. The martingale residual of a subject is its status minus its expected number of events H₀(t_i) e^β′x_i; the residuals sum to zero and are skewed, so the deviance residual sign(M)√{−2[M + δ ln(δ − M)]} is used to look for outliers.

Checking proportional hazards

The Schoenfeld residual of an event is its covariate vector minus the risk-weighted mean of the covariates over the risk set at that time. If the hazard ratio of a covariate is constant the residuals show no trend in time. The page performs the test of Grambsch and Therneau in the form of R's cox.zph: for each term it computes the score test of adding the term multiplied by a function g(t) of time (centred at its mean over the events) to the fitted model, which is the test for a slope of the scaled Schoenfeld residuals against g(t); the GLOBAL line adds all the interactions at once. g may be one minus the Kaplan–Meier estimate (the default, robust to outlying times), time, log time or the rank of time. A small p-value for a term means its effect changes during follow-up: stratify on it, split the follow-up, or read its coefficient as an average effect. The complementary log-log plot gives the same information by eye for the groups: under proportional hazards the curves of ln[−ln S(t)] against ln t are a constant distance apart.

Parametric fits

Each group is also fitted by maximum likelihood with three lifetime distributions, an event contributing its density and a censored time its survival probability: the exponential S = exp(−λt), whose estimate is events divided by total time at risk; the Weibull S = exp[−(t/η)^k]; and the lognormal S = 1 − Φ[(ln t − μ)/σ]. The Weibull and lognormal are fitted in the log-location-scale form ln T = μ + σW of R's survreg by Newton–Raphson (Weibull shape k = 1/σ, scale η = e^μ). AIC = −2ℓ + 2 × (number of parameters) compares the three on one group; on the complementary log-log plot a Weibull is a straight line of slope k and an exponential a line of slope 1. A parametric curve extrapolates beyond the data and gives smooth quantiles, at the price of an assumption the Kaplan–Meier estimate does not make.

Planning a log-rank trial

The power of the log-rank test is set by the number of events. For a hazard ratio HR, two-sided level α, power 1 − β and fractions p_c, p_x of the patients in the control and experimental arms, Schoenfeld's formula requires d = (z_1−α/2 + z_1−β)²/[p_c p_x (ln HR)²] events; Freedman's requires d = (z_1−α/2 + z_1−β)² (1 + φ HR)²/[φ (1 − HR)²] with allocation ratio φ = p_x/p_c, somewhat more for hazard ratios far from 1. Dividing by the average probability that a patient has the event during the study, p_c P_c + p_x [1 − (1 − P_c)^HR], gives the number of patients; increase it for losses to follow-up.

Verification

The built-in data sets are the ones distributed with R: the remission times of the 6-mercaptopurine trial of Freireich and others, the AML maintenance trial, an ovarian cancer trial, the NCCTG lung cancer study and the Veterans Administration lung cancer trial. On them the page reproduces, to the printed digits, the output that the R packages survival and MASS publish as their reference results: Kaplan–Meier tables and medians, the log-rank and G(ρ) tests, Cox coefficients, standard errors, tests, residuals and concordance, the proportional-hazards tests, and the Weibull fits. For the default data the log-rank χ² is 16.79 and the Cox hazard ratio of control to 6-MP is 4.82 with Efron's ties (4.52 with Breslow's).