Statistical methods for censored survival data.

Methods of statistical analysis of censored survival times are briefly reviewed and illustrated by application to clinical trials data. These include estimation of the survival curce, nonparametric tests to compare several survival curves, tests for trend, and regression analysis. Extensions of the methodology are made for application to epidemiologic case-control studies. These are used to estimate relative risks for leukemia associated with radiation exposures. A final section provides some annotated references to the recent literature.


Introduction
Censored survival data arise in a wide variety of statistical investigations. In clinical trials one measures duration of response from start of treatment until relapse or death due to disease. Observations on response time are censored for those subjects still in remission at the study's end, as they are for patients lost to followup during the course of the study. Animal carcinogenesis studies, such as used by the United States Food and Drug Administration to determine the safety offood additives, provide another example. Here the endpoint is the age at diagnosis of a particular kind of cancer, censorship being imposed by death due to other causes, natural or artificial. In tests of the reliability of airplane components, failure times are measured from the start of testing until failure of the component, with censorship imposed by the failure of other components or the necessity of analyzing the data before all items have failed. Figure 1 illustrates the results for the control group in a clinical trial designed to investigate the effects of combined chemotherapy as an adjunct to surgery and radiation in the treatment of childhood rhabdomyosarcoma (1). The endpoint for analysis was the reappearance of tumor, whether at the site of original treatment or through distant metastasis. Children who remained disease-free at the time the *Department of Biostatistics SC-32, University of Washington, Seattle, Washington 98195. data were analyzed had censored observations. In addition to the control arm IA, there were two groups of children who received the drugs actinomycin-D (AMD) and vincristine (VCR): group IB patients were concurrently randomized with the controls, both these groups having apparently had their tumors completely resected; group IIA consisted of patients with microscopic residual disease at the margin of surgical resection.
Interim data from all three arms are presented in Table 1. Note that the censored observations for arm IA, those in the column labeled "disease-free", are smaller in the table than they are in the figure. This is because the figure was drawn from data computed at a later point in time, when additional follow-up was available for patients who had not already died.
Analysis of such data has several goals. For each of the comparison groups one wants an estimate of the survival curve, the probability of surviving t units of time. Statistical tests are required to determine whether the observed differences between the curves are real, or are simply chance effects. If real, a method of quantifying the nature of the differences is desirable. Finally there may be available concomitant observations, including continuous measurements such as age at diagnosis, whose joint effects on survival are important to determine.

Estimation of Survival Curves
When analyzing several groups of survival times the first step is to form a series of 2 x r tables as  9  12  37  25  3  15  16  19  28  9  18  19  20  29  10  24  20  38  10  36  24  42  15  40  24  45  16  45  30  47  30  31  48  34  50  42  52  44  53  59  62 refer to the total number of subjects in the ijh group who remain "at risk", i.e., alive and under observation, just prior to time tk. The tabular entries dik and Sik denote the numbers of those who die at tk, and survive tk, respectively. Table 3 illustrates the calculation of the first three such tables for the data in Table 1. Here r = 3 and ti = 2, t2 = 3, and t3 = 9 months. Note that the tables for increasing tk refer to a constantly diminishing population "at risk" as additional subjects die or are withdrawn (censored) from further observation. Kaplan and Meier (2) derived the maximum likelihood nonparametric estimate of the survival curve based on censored data. This may be calculated recursively, starting from P(to) = 1, and by using the formula (1):  for k = 1, 2, . . ., K. (The group index i has been suppressed for clarity.) In other words, the probability of surviving past tk is estimated as the probability of surviving past tk-l times the conditional probability of surviving past tk given survival to tk. Because of this multiplicative structure, Kaplan and Meier refer to their estimate as the product limit (PL) estimate. In case there is no censorship in the data, it reduces to the familiar empirical distribution function. Table 4 shows the calculation of the relapse-free survival curve from the interim data in Table 1 for treatment group IA. The corresponding curves calculated from final study data for all three treatment curves are shown in Figure 2. Numbers above each curve at annual intervals in this figure refer to numbers of patients still at risk in each group. These are an important means of judging the stability of the estimates, which can in fact be quite unstable in the "tail" ofthe survival distribution where few subjects remain at risk.
The variance of the PL estimate may also be cal- with the understanding that it is applied se;quentially to tied observations. In large samples, P(t) is approximately normally distributed with mean equal to the true survival function P(t) and a variance estimated as shown above (2,3). Note that neither P (t) nor V{P (t)} will change after the last uncensored response time in each group, even though additional subjects continue to be withdrawn from observation. In this region the estimated variance often does not accurately reflect the true variability in the survival curve, which will be substantial unless large numbers of subjects remain on study.

Comparison of Survival Curves
A simple but powerful non-parametric test for the comparison of r survival curves with censored data may also be calculated from the basic data shown in Table 2. This test exploits the fact that, under the null hypothesis of no difference in the underlying survival distributions and conditional upon fixed values for the marginal totals in each 2 x r table, the vector dk = (dlk, . . ., drk)' of observed deaths at tk has an r-dimensional hypergeometric distribution. Consequently the null expectation of the number of deaths in group i at tk is eik = E(dik) = nik (Dk,Nk) (3) i.e., the number at risk in the i-th group times the death rate for all r groups combined (see Table 3 for October 1979 an illustration of this calculation). The covariance matrix Vk ofdk has, under the null hypothesis, an (i,j) component equal to The main idea behind the test is to sum up the statistics calulated from each of the K tables into a vector 0 = Ykdk of observed numbers of deaths in each group, a vector E = :kek of expected numbers of deaths, and a summary covariance matrix V = IkVk. Since the 2 x r tables refer to overlapping sets of subjects they are not, strictly speaking, statistically independent. Nevertheless Cox (4) has shown that V is an appropriate large sample covariance matrix for O-E. Since YOi = ;Ei, i.e. the totals of observed and expected deaths in all r groups agree, V is singular. However, defining 0* and E* to be the first r -1 components of 0 and E, and V* to be the (r -1) x (r -1) upper left hand corner of V, a test statistic for testing equality of the r survival curves is obtained as This is approximately distributed as chi-square on r -1 degrees of freedom under the null hypothesis.
The test Ti was first proposed for survival data by Mantel (5). Cox (6) later derived it from likelihood theory under the proportional hazards (PH) model, in which the instantaneous death rates in the r groups are assumed to be in constant ratio throughout the follow-up period (see below). Peto and Peto (7) argued that it was an asymptotically efficient test under Cox's model and named it the "log rank" test.
A conservative approximation to Ti which requires no matrix inversion is given by the familiar chi-square formula While T2 ' Ti, in fact the two statistics will be quite close, provided that there are few ties among the uncensored survival times (i.e., most of the Dk in Table 2 are unity) and that the patterns of censorship operating in the r groups are not grossly different (8,9). Note that the ½ continuity correction should not be used with survival data. Table 5 illustrates the manner of presentation of the summary and test statistics. Note the calculation of the ratio O/E of observed to expected numbers of deaths in each treatment group. These are very useful as measures of treatment effect, since their ratios, e.g., OiIEl . 021E2, approximate the ratios of death rates in the respective treatment groups (10).

Alternate Weighting Schemes
The summary statistics 0 -E weight the observed differences dkek in each table in a manner which is appropriate to the PH model already mentioned. However this is not the only possible weighting scheme. Multiplying the observed differences by Nk, the total number of subjects in the k-th table, and then summing, gives more weight to the earlier times tk when larger numbers are at risk. This leads to the scores and test statistic T3 = W*'Vw*-lW* (7) (8) (9) where asterisks (*) denote the corresponding r -1 dimensional quantities. A conservative approximation to T3 not requiring matrix inversion is Environmental Health Perspectives ( 1 1) The scores Wi may also be obtained from a pairwise comparison of the observations in the i-th treatment group with those in the remaining r -1 groups. Each such pair is assigned the value + I (or -1) according as the true survival time for the first pair member is known to be smaller than (or larger than) that for the second member. Ties or indeterminate comparisons due to censorship are assigned 0 values. Gehan (11) suggested the use of such scores for the comparison of two samples (r = 2). In this case T4 reduces to the familiar Wilcoxon rank sum test in the absence of ties and censorship. Breslow (12) extended this work to r > 2 samples, proposing the covariance matrix Vw and the statistic T3. These latter statistics, like V and Ti, are valid for situations where the patterns of censorship operative in the r treatment groups are unequal, as in animal carcinogenesis studies where there is differential toxic mortality. The conservative approximation T4, like T2, is strictly valid only where there is equality of censorship.
In practice, the tests Ti and T3 often yield rather similar numerical values (see Table 5). However, this is not always true, and some comment on the proper interpretation when only one statistic is significant is in order. Since T3 weights early values more heavily, it may achieve significance when there is an early separation between the survival curves which later come together or even cross over. Ti gives more weight to the later appearing deaths. A large discrepancy between Ti and T3 generally indicates an interaction between treatment and time on the instantaneous death rates, which is worthy of investigation in its own right.

Testing for Trend
Often the r comparison groups correspond to r levelsxl <X2 < . . . Xr of a quantitative variable such as dose. Global chi-square tests such as Ti through T4 lack statistical power in such situations since they take no account of the natural order of the groups. One needs a single degree of freedom test for trend in survival with increasing dose.
Fortunately such a test is readily calculated from the summary statistics already at hand. In the case of the log rank analysis Vx is a single degree of freedom chi-square for a linear trend of O-E withx. Tarone (13) has suggested using T6 = Ti -T5 as a chi-square on r-2 degrees of freedom for deviations from linearity. An approximation to T5 which only requires calculation of the 0 and E vectors is given by Similarly, when using the W scores, provides a test for linear trend and T9 = T3 -T7 a test for deviations from linearity.
To illustrate these calculations by using the summary data in Table 5, make the fictitious assumption that the three treatment groups IA, IB, and IIA correspond to dose levels xi = 0, X2 = 1, X3 = 2. The statistics for trend are T5 = 9.17 (p = 0.002), T7 = 8.72 (p = 0.003), and T8 = 10.02 (p = 0.002). The deviation chi squares are T6 = 1.93 (NS) and Ts = 1.39 (NS). Hence, this would be a case where the already significant differences are largely explained on the basis of an apparent linear trend in survival with increasing dose.

Adjustment by Stratification
When the r comparison groups differ with respect to factors which influence survival, an analysis which corrects for their possible confounding effects is needed. This may be carried out very simply by dividing the population into strata which are more or less homogeneous internally with respect to the confounding factors. (Of course the number of confounders which may be accommodated simultaneously in this fashion is limited, since if strata become very large in number and small in size a substantial loss of comparative information may result.) Separate survival analyses are performed within each stratum by calculating the summary statistics 0, E, and V defined earlier. These are cumulated by simple addition over strata and used to calculate adjusted test statistics Ti, T2, T5, T6 and T7, in which the cumulated summary statistics replace the stratum specific ones. Likewise adjusted versions of T3, T8, and Ts use the cumulated W and Vw statistics.
Such a stratified analysis was used for a trial of maintenance chemotherapy for children with acute lymphocytic leukemia (14). For this disease it is well known that the diagnostic white blood count (WBC) is an important prognostic factor. The treatment group, consisting of 152 children who received actinomycin (AMD) in addition to standard maintenance drugs, had a median WBC of 10,067; whereas the control group, consisting of 116 children who did not receive AMD, had a median WBC of 14,280. An analysis ignoring this difference in WBC compared the observed number of relapses in the treated and control groups, O1 = 100 and 02 = 81, with expected numbers ofEl = 113.00 and E2 = 68.00. This yielded an (unadjusted) chi-square of T2 = 3.98, which isjust on the borderline of 5% statistical significance.
In order to determine whether the apparent effectiveness of AMD was due, at least in part, to the generally lower WBC's among treated patients, the entire sample of 268 children was divided into three strata as shown in Table 6. A separate calculation of the observed and expected numbers of relapses was made within each stratum, so the totals of O's and E's taken across each row of Table 6 agree. The adjusted expected numbers, El = 110.96 and E2 = 70.04, are now closer to the observed numbers and give an adjusted chi-square of T2 = 2.80, which is no longer statistically significant.
The estimated ratio of relapse rates in the treated vs. control group is 100/113.0 . 81/68.0 = 0.74. When adjusted for WBC in three strata, it is 100/ 110.96 . 81/70.04 = 0.78, again indicating less of a difference between treatment and control group.

Regression Analysis of Survival Data
If the number of confounding variables is large, stratification breaks down since there will be many strata with just one or a few subjects. It may also be of interest to quantify the relationship between survival and several discrete or continuous and con- The usual regression model specifies that the survival times, or some transform such as their logarithm, are equal to a linear combination of the concomitant variables plus a random error term. Unfortunately, to generalize such models for use with censored data is awkward and computationally involved. Thus considerable interest was aroused by Cox (6) when he proposed a model formulated in terms of the effect of the regression variables on instantaneous death rates rather than on times of deathper se. This model turned out to be quite tractable computationally and, as an added benefit, avoided any parametric assumptions about the shape of the underlying survival curve.
Cox's model is defined in terms of the time t specific death rate or hazard function X(tlz) for an individual having ap-vector of covariates z. Specifically he assumes X(tlz) = exp(,f'z)Xo(t), where ,3 is an unknown p-vector of parameters (regression coefficients), while Xo(t) is the unknown hazard or death rate function for an individual with a standard (z = 0) set of convariates. A consequence of this model is that the ratio of hazard functions for two individuals with different sets of covariates, X(tlzi) -= exp {,B' (zl -Z2)} X(tIz2) (15) does not depend on time. Thus it is called the proportional hazards (PH) model.
Several authors (4, 6, 10, [15][16][17] have developed the likelihood analysis of the PH model from rather distinct points of view. Providing that there are no ties in the uncensored data, all derive for the lnlikelihood function of the expression Environmental Health Perspectives where R(tk) is the risk set of subjects still alive and under observation just prior to tk; Zk is the covariate vector for the individual who dies at tk; and the outer summation is over all K true (uncensored) times of death. Different likelihoods arise in the case of ties.
Taking the vector of first partial derivatives of L, setting equal to 0 and solving the resulting nonlinear equations yields a maximum likelihood estimate 3 for the regression coefficients. A covariance matrix for this estimate is obtained in the usual fashion by inversion of the negative of the matrix of second partials of L. The integral Notice that when ,3 = 0 this reduces to the PL estimate of Kaplan and Meier, calculated from the entire set of observations considered as one homogeneous sample. Table 7 illustrates the computer fitting of the PH model to the data on 268 leukemic children. Four regression variables were considered: zi = log (WBC), Z2 = age at diagnosis (years), Z3 = Z2, and Z4 = 1 or 0, according as the patient was treated with AMD or was a control. These four variables were entered into the regression equation sequentially in order to demonstrate their effects on remission duration after adjustment for the preceding variables. A quadratic term in the age variable was required: children in the mid ranges from 2 to 6 years have a better prognosis than at either extreme. The multiplicative effect of treatment on the relapse rate is given by exp (J34) = exp (-0.220) = 0.80, which is quite comparable with the approximate value 0.78 obtained from the simpler stratified analysis. The likelihood ratio test for treatment effectiveness yields a chi-square of 2(-876.95 + 877.99) = 2.09; squaring the standard-ized regression coefficient gives a similar value (_1.45)2 = 2. 11. These are both even smaller than the value of 2.80 obtained after adjustment for WBC in three strata, so that the regression approach has in this case led to an even greater reduction in the statistical significance of the treatment comparison.
As a means of providing a graphical display of the fit of the model, the regression coefficients were used to calculate a prognostic score for each child using the formula S = Z1f31 + Z2f32 + Z3133 + Z4/34 (21) Each of the covariates z was first normalized by subtracting off the mean, so that a score S = 0 represents a "typical" patient. The scores were then used to divide the sample into four groups within each of which PL estimates of the remission duration were calculated. These are plotted in Figure 3 together with predicted remission duration curves, estimated from the model, for specified values ofS. Notice that the predicted curves lie further apart than do the "observed" curves for later days, while the reverse is true for earlier days. This behavior indicates a certain lack of fit of the model, namely that the baseline covariates have more of an effect on early rates of relapse than they do on later ones. The fit could be improved by use of time-dependent covariates z(t) as discussed by Cox (6).

The PH Model in Epidemiology: Applications to RERF Data
The methodology of survival analysis, especially the PH model developed in the last section, is also useful with epidemiologic studies of risk factors for chronic disease. In this context t often represents age, and the endpoint is diagnosis of, or death from, a particular disease in a previously disease-free individual. Xo(t) may then be interpreted as the agespecific incidence or mortality rate for a standard covariate set, while exp (B'z) represents the relative risk or rate ratio (RR) for a subject with covariates z. These are computed from the presumed risk factors under investigation and may themselves depend on age. Application of the previously discussed methodology based on the PH model is straightforward, at least in principle. A slight modification is that the risk sets R(tk) may change with age not only due to the loss of individuals from further observation, as in clinical trials, but also from the enrollment of new subjects in the study at later ages.
Prospective epidemiologic studies involve a cohort of disease-free persons who are enrolled at various ages and kept under continuous surveillance until they either develop the disease, or else are lost to further observation or die from another cause. In order to collect enough cases of rare chronic diseases, tens or even hundreds of thousands of persons may have to be followed over several years. For example, the Radiation Effects Research Foundation (RERF) has kept nearly 100,000 survivors of the Hiroshima and Nagasaki atomic bomb blasts under surveillance for nearly three decades, losing fewer than 0.1% to emigration. With such large cohorts the previously discussed iterative methods for estima-tion and testing of the parameters ,8 and Xo are neither feasible nor necessary.
One method of dealing with large cohorts is to partition the time or age axis into intervals, say ten or twenty, and to postulate a discrete time model for the conditional probabilities of failure within each one. Suppose there areK such intervals (0,tl[, (tl, t2], (tK1, tK], and set Pk(Z) = P[T S tkIT < tk-1, Z] (22) for k = 1, 2, . . ., K, where T denotes the random failure time (age at diagnosis or death), for a subject Environmental Health Perspectives with covariates z. Such an individual contributes a termpk(z) to the likelihood if he develops the disease in the k-th interval, and a term {1 -Pk(z)} if he survives it. If he dies of other causes (is censored) in the interval, Thompson (18) makes the sensible recommendation that he contribute a term 71 -Pk(Z), to reflect his survival disease-free over only a part of the interval.
Cox (6) suggested the use of the linear logistic model logit Pk(Z) = ak + /3'Z (23) where logit (p) = In in/(lp)} (24) This reduces to the PH model in the limit as K increases and the time or age intervals become infinitely small. The term exp (,3' z) in Eq. (23) represents the odds ratio of disease occurrence in each interval, rather than the ratio of instantaneous failure rates.
A similar model was one of several used by Otake (19) to explore the effects of radiation on causespecific mortality using RERF data. He divided the sample into five groups according to age at the time of bomb (ATB) and six classes according to estimated total radiation doses, and considered a variety of causes of death as the endpoint. If nij denotes the number of subjects in age group i and radiation class j, while dij denotes the number of deaths due to specific cause, then his model states logit E(dijln ,) = ai + 1j subject to appropriate constraints on the parameters. This is a linear logistic model for the unconditional probabilities of death over a defined period (in this case 1950-1972) and so does not explicitly account for competing causes ofdeath or differential losses to further observation which may be taking place in the 30 age x dose cells.
Cox's model based on the conditional probabilities would further divide the follow-up period into K intervals, say 1950-1954, 1955-1959, etc. If nik then denotes the number of subjects who were age i ATB, received radiation dosej, and who remained alive at the midpoint of the k-th interval, the model is logit E(diJkln/Jk) = ai +/ + Yk (27) One might include also an interaction term (ay)ik so as to allow age and calendar time to have essentially arbitrary effects on disease incidence. The numbers of deaths dijk may be formally considered to have independent binomial distributions conditionally on the nik (4).
A different approach to this problem, termed by Mantel a "synthetic retrospective study," is to draw a small sample of disease-free "controls" from each of the risk sets R(tk) for comparison with the diseased cases as they arise. The very same theory and methods may be applied also to actual case/control studies in which cases are ascertained as they occur, for example through a population based disease register. The controls, instead of being obtained from a computer file, are samples from the population in which the case arose. They may be patients with different diagnoses in the same hospital, chosen from the same family or neighborhood, or simply sampled at random from the population.
While it is impossible to estimate the incidence rates Xo(t) without studying the full cohort, the PH model assumed for the prospective study implies a probability structure for the sample cases and controls which may be used to estimate the parameter,8 describing the RR associated with the covariates (20). Specifically, suppose m = m(t) cases having covariate vectors zl, ..., Zm are diagnosed (or die) with the particular disease at age t. Suppose also that n disease-free controls with covariates Zm+i, * * * Zm+n are drawn at random from the risk set. The analysis is performed conditionally on fixed values for m, n and the combined set of n + m observed covariate vectors. Using Cox This expression is formally identical to that given earlier for the PH model in survival analysis. However now the risk sets R(tk), instead of containing everyone still alive and disease-free at tk, are replaced by the much smaller sets of mk + nk subjects actually sampled for the retrospective analysis.
In practice, the sets of m + n cases and controls may be matched on other variables x besides age, for example on risk factors already known to be associated with the disease under investigation. The PH model may then be generalized to X(t|z,x) = exp {f,'z} Xo(tlx) (31) which allows the underlying incidence function Xo to depend in an arbitrary way on x as well as t. Interactions between the matching variables and the risk factors continue to be modelled in z. Retrospective sampling carried out within risk sets having similar x and t values leads to the same conditional likelihood given above (20).
A special case of the conditional likelihood occurs when each case is individually matched to exactly R controls, a situation found often in practice. Suppose there are K such sets and denote by ZOk the covariate vector for the case and by Zlk, .. ., ZRk the covariates This approach was taken with an RERF data file consisting of records on 7078 males who were 50+ years old ATB. Thirteen deaths from leukemia occurred (Table 8). For each of these a search was made to determine the risk sets of potential controls who had the same age ATB (exact year), were residing in the same city (Hiroshima or Nagasaki), and who were alive at the end of the same year in which the case died. These ranged in size from 30 to 508 individuals.
In order to illustrate the application of the case/ control methodology, a series of 1, 2, 5, 10, and 20 controls was drawn at random from each of the risk sets. Two covariates were computed from the estimated radiation dose (in rads) for each subject: zl = ln (rads + 1) and Z2 = JiaiTs. The choice of ln (rads + 1) as a covariate in the PH model, which implies there is a linear increase in ln RR with In (rads + 1), was based on two facts: the highly skewed distribution of radiation doses; and Figure 2 of Otake (19), which shows a nonlinear increase in In RR with dose after about 200 rad. The alternate transformation Vrads was used for comparison.
Results of the matched analyses based on the dif- Table 8. Summary data on records used for matched case-control analysis of RERF files.

Risk set Age
Year of Dose, Test for quadratic effect Xi = 7.34 ferent numbers of controls are shown in Table 9. Note the decrease in the estimated standard errors as additional controls are used; however the gain afforded by using twenty rather than ten controls is not great. The regression coefficients of 0.4-0.5 for ln (rads + 1) indicate that, roughly speaking, the risk of leukemia increases by about 5% for every 10% increase in radiation dose. The disparity between the relative risks fitted under the model z = ln (rads + 1) vs. those fitted by z = rads is quite noticeable, especially when one recalls that most doses fall in the 0 -200 range. In order to try to understand better why this might occur, unmatched logistic regression analyses were carried out on the data consisting of all cases and the sets of twenty controls. The estimated slopes and standard errors were similar to those obtained with the matched analysis (Table 10). However the intercepts, corresponding to the estimated log odds of leukemia in the sample at a dose of zero rads, differed markedly between the two models: a = -3.97 for z = ln (rads + 1) vs. a = -3.33 for z = 7rads.
Differences in the estimated relative risk between the two models were thus explained largely by the differences in the absolute risks estimated for the baseline value of the covariate.
A potential hazard of the regression modelling of relative risks is its sensitivity to the choice of scale on which the covariate is measured. Goodness of fit of the model to the data is essential for proper interpretation, and should be explored thoroughly. When both linear and quadratic covariate terms were fitted with the models above, for example, the agreement between the estimated values for a improved substantially. Moreover the fact that the quadratic term was highly significant for the /ads model, and not at all significant for the In (rads + 1) model, showed that the latter gave much better fit to these data (Table 10).

Further Reading
Much of the above material is presented in greater detail in a review article (10) on the PH model and its applications to survival data. Some additional applications of this model to epidemiologic studies are given by Breslow et al. (21,22). Peto et al. (23) present a thorough discussion of its use in the design and analysis of clinical trials.
A computer program for calculating the PL estimate and all the test statistics presented in sections 2-5 above is available (24).
Several authors have pointed out that the W scores defined do not lead to the most efficient generalization of Wilcoxon's test to censored data. They all propose essentially the same statistic as an alternate generalization (7,25,26).
A comparison of the efficiencies of the test statistics using Monte Carlo techniques is made by Lee et al. (27). Efron (17) discusses the efficiency of the likelihood function used with the PH model from a more abstract viewpoint; see also Kalbfleisch (28).
Additional extensions of the PH regression model for use with grouped or heavily tied data are discussed by Cox (6), Kalbfleisch and Prentice (15), Thompson (18), and Prentice and Gloeckler (29).