Biostat3.net

British Journal of Cancer (2003) 89, 431 – 436& 2003 Cancer Research UK All rights reserved 0007 – 0920/03 $25.00 Survival Analysis Part II: Multivariate data analysis – an introductionto concepts and methods




MJBradburn*,1,TGClark1,SBLove1 andDGAltman1 Cancer Research UK/NHS Centre for Statistics in Medicine, Institute of Health Sciences, Old Road, Oxford OX3 7LF, UK British Journal of Cancer (2003) 89, 431 – 436. doi:10.1038/sj.bjc.6601119 & 2003 Cancer Research UK Keywords: survival analysis; Cox model; AFT model; model selection groups, it offers no estimate of the actual effect size; in otherwords, it offers a statistical, but not a clinical, assessment of the Survival analysis involves the consideration of the time between a factor's impact. The use of a statistical model improves on these fixed starting point (e.g. diagnosis of cancer) and a terminating methods by allowing survival to be assessed with respect to several event (e.g. death). The key feature that distinguishes such data factors simultaneously, and in addition, offers estimates of the from other types is that the event will not necessarily have strength of effect for each constituent factor. Therefore, statistical occurred in all individuals by the time the study ends, and for models are important and frequently used tools which, when these patients, their full survival times are unknown. For instance, constructed appropriately, offer valuable insight into the survival in studies that measure the length of survival after diagnosis of cancer, it is common for a proportion of individuals to remain Several statistical methods have been proposed for modelling alive and disease-free at the end of the follow-up period, and for survival analysis data. We will describe the most important models these patients, we know only a lower limit on their actual time to and illustrate their application using example datasets described in event. Thus, special methods are required for these type of data.
the previous paper (Clark et al, 2003). As before, we will assume The explanation and demonstration of some of the methods throughout that all survival times are independent of each other proposed to analyse such data are the basis of this series.
and that censoring occurs solely as right-censoring and is In the first paper of this series (Clark et al, 2003), we described uninformative. The focus is on covariates that are measured at initial methods for analysing and summarising survival data the time of entry to the study, that may be continuous (e.g. the including the definition of hazard and survival functions, and patient age or tumour size), binary (e.g. gender), unordered testing for a difference between two groups. We continue here by categorical (e.g. histology) or ordered categorical or ordinal (e.g.
considering various statistical models and, in particular, how to performance status or FIGO stage). In the next paper in this series, estimate the effect of one or more factors that may predict survival.
we will discuss the statistical assumptions made when usingstatistical models, and provide advice on choosing the appropriatemodel and covariates therein. We will also consider how to model THE NEED FOR MULTIVARIATE STATISTICAL covariates that change values over time (called ‘time-dependent' or The methods we present here may be divided into two broad The previous paper demonstrated the construction of (Kaplan – categories: proportional hazard approaches (including the semi- Meier) survival curves for different patient groups, and introduced parametric Cox model and fully parametric approaches) and the logrank test to investigate differences between them. Both these accelerated failure time models. These methods have different methods are examples of univariate analysis; they describe the properties and interpretations, but all may be used to summarise survival with respect to the factor under investigation, but survival data.
necessarily ignore the impact of any others. It is more common,at least in clinical investigations, to have a situation where several(known) quantities or covariates, potentially affect patient prog-nosis. For example, suppose two groups of patients are compared: THE COX (‘SEMI-PARAMETRIC') PROPORTIONAL those with and those without a specific genotype. If one of the groups also contains older individuals, any difference in survival The Cox (proportional hazards or PH) model (Cox, 1972) is the may be attributable to genotype or age or indeed both. Hence, most commonly used multivariate approach for analysing survival when investigating survival in relation to any one factor, it is often time data in medical research. It is a survival analysis regression desirable to adjust for the impact of others. Moreover, while the model, which describes the relation between the event incidence, logrank test provides a P-value for the differences between the as expressed by the hazard function and a set of covariates. A fullerexplanation of the hazard function was given in the previous *Correspondence: Mr M Bradburn; E-mail: [email protected] article (Clark et al, 2003). Put briefly, the hazard is the Received 6 December 2002; accepted 30 April 2003 instantaneous event probability at a given time, or the probability Multivariate data analysis MJ Bradburn et al that an individual under observation experiences the event in a stage (an ordinal covariate taking values of 1, 2 3 or 4), histology period centred around that point in time.
(one of seven subtypes), grade (1, 2 or 3), ascites (yes/no) and Mathematically, the Cox model is written as patient age.
Table 1 shows the effect sizes (given as hazard ratios), 95% hðtÞ ¼ h0ðtÞ expfb1x1 þ b2x2 þ    þ bpxpg confidence intervals (CI), regression coefficients and statisticalsignificance for each of these in relation to overall survival. Each where the hazard function h(t) is dependent on (or determined by) factor is assessed through separate univariate Cox regressions a set of p covariates (x1, x2, y, xp), whose impact is measured by (left-hand columns). However, the aim of the database is to the size of the respective coefficients (b1, b2, y, bp). The term h0 is describe how the factors jointly impact on survival, and so all five called the baseline hazard, and is the value of the hazard if all the xi factors were incorporated into the multivariate model (right-hand are equal to zero (the quantity exp(0) equals 1). The ‘t' in h(t) columns). It may be seen that higher FIGO stage, higher grade, reminds us that the hazard may (and probably will) vary over time.
presence of ascites and increased age impaired survival to varying An appealing feature of the Cox model is that the baseline hazard (and statistically significant) degrees. The histology was also of function is estimated nonparametrically, and so unlike most other importance: the figures describe the survival of patients with each statistical models, the survival times are not assumed to follow a histology type in comparison with the serous type. In principle, particular statistical distribution.
any type with a reasonable number of patients could be chosen as The Cox model is essentially a multiple linear regression of the the baseline of comparison. On multivariate analysis Mucinous logarithm of the hazard on the variables xi, with the baseline and serous were the tumour types with the best prognosis, whereas hazard being an ‘intercept' term that varies with time. The undifferentiated and mixed mesodermal were the worst. It is covariates then act multiplicatively on the hazard at any point in possible to present P-values for the comparisons between each time, and this provides us with the key assumption of the PH type and serous, but we have given an overall likelihood ratio test model: the hazard of the event in any group is a constant multiple for the differences between the categories as a whole. The FIGO of the hazard in any other. This assumption implies that the hazard stage could be modelled as a categorical variable in the same curves for the groups should be proportional and cannot cross (see manner as grade and histology, but assuming it is a continuous Figure 1 for examples of each). Proportionality implies that the variable with a linear trend across the four categories performed quantities exp(bi) are called hazard ratios. A value of bi greater sufficiently well.
than zero, or equivalently a hazard ratio greater than one, indicatesthat as the value of the ith covariate increases, the event hazardincreases and thus the length of survival decreases. Put another PARAMETRIC PH MODELS way, a hazard ratio above 1 indicates a covariate that is positivelyassociated with the event probability, and thus negatively Parametric PH models are a class of models similar in concept and associated with the length of survival. This proportionality interpretation to the Cox (PH) model. The key difference between assumption is often appropriate for survival time data but it is the two is that the hazard is assumed to follow a specific statistical important to verify that it holds. We discuss methods for assessing distribution when a fully parametric PH model is fitted to the data, proportionality in the next paper in this series.
whereas the Cox model enforces no such constraint. Other thanthis, the two model types are equivalent. Hazard ratios have thesame interpretation, whether derived from a Cox or a fully The Cox PH model fitted to the ovarian cancer data parametric regression model, and the proportionality of hazards is This large database, as described in the previous paper of this still assumed.
series (Clark et al, 2003), was used to derive a prognostic index for A number of different parametric PH models may be derived by overall survival among ovarian cancer patients in Clark et al choosing different hazard functions. As shown previously, there is (2001). Their analysis included 10 variables, but for simplicity we a direct link between the survival and hazard, and the choice of will consider five, all of which were measured at diagnosis: FIGO hazard distribution determines that of the survival. In fact, themodels commonly applied, such as the Exponential, Weibull orGompertz models, take their names from the distribution that thesurvival times are assumed to follow, but the most distinguishing features between them are in the hazard function. Examples of survival and hazard functions derived from some of theseparametric models were presented in the previous paper (Clark et al, 2003). Figure 1 shows increasing and decreasing Weibullhazard functions, as well as two groups with the latter that are proportional to each other.
Parametric models fitted to the ovarian cancer data The estimated hazard function of the ovarian cancer data asdisplayed in the previous paper (Clark et al, 2003) may be consistent with that derived from a Weibull PH model with decreasing hazard. Fitting this to the ovarian cancer database gives similar results as the Cox model (see Table 2), and may beinterpreted in the same manner. Methods to check for the appropriateness of the Weibull distribution will be discussed in the next paper of this series.
Example of (non-) proportional hazards (groups (c) and (d) only have proportional hazards) using the Weibull distribution. For theWeibull survival model, the hazard function h(t) ¼ ls(lt)s 1 for l, s40: (a) COMPARISON OF THE TWO PH APPROACHES increasing hazard (l ¼ 0.5, s ¼ 1.25); (b) decreasing hazard (l ¼ 0.25,s ¼ 0.75); (c) decreasing hazard (l ¼ 0.5, s ¼ 0.5); (d) decreasing hazard The main drawback of parametric models is the need to specify the (l ¼ 0.25, s ¼ 0.5).
distribution that most appropriately mirrors that of the actual British Journal of Cancer (2003) 89(3), 431 – 436 & 2003 Cancer Research UK Multivariate data analysisMJ Bradburn et al Hazard ratios from the Cox PH model for the ovarian dataset Univariate analysis Multivariate analysis Absence of ascites Age (per 5-year increase) HR ¼ hazard ratio, CI ¼ confidence interval.
Hazard ratios from the Weibull PH model for the ovarian dataset Univariate analysis Multivariate analysis Absence of ascites Age (per 5-year increase) HR ¼ hazard ratio, CI ¼ confidence interval.
survival times. This is an important requirement that needs INTERPRETING THE PH MODEL: BEYOND THE to be verified and an appropriate distribution may be difficult to identify. Where a suitable distribution can be found,however, the parametric model is more informative than the In addition to the ratio of two hazards, it is possible to obtain other Cox model. It is straightforward to derive the hazard function and information from a PH regression model. One simple (and to obtain predicted survival times when using a parametric model, possibly underused) quantity that may be derived from a survival which is not the case in the Cox framework (the use of such model is the predicted survival proportion at any given point in quantities is discussed in the next section). Additionally, the time for a particular risk group. The survival proportion for a appropriate use of these models offers the advantage of being given risk group at any time, S(t), is equal to slightly more efficient; they yield more precise estimates (i.e.
smaller standard errors).
SðtÞ ¼ S0ðtÞexpðgÞ The results from the Cox or parametric PH models may be compared directly, as the model types are merely different where S0(t) is the baseline survival (the survival proportion when approaches to assessing the same quantity. For either method to be valid: (a) the covariate effect needs to be at least approximately b1x1 þ b2x2 þ     þbpxp. Once the value of the baseline survival constant throughout the duration of the study, and (b) the at a given time is derived, then the predicted survival probabilities proportionality assumption must hold. These important issues will for patients with any specified covariate values xi are easily be addressed in the subsequent paper in this series.
obtained. This information could then be displayed via tabular or & 2003 Cancer Research UK British Journal of Cancer (2003) 89(3), 431 – 436 Multivariate data analysis MJ Bradburn et al Survival proportion Survival probability Patient age (years) Illustration of the AFT model: (——), S0(t) the baseline survival function; (       ), S(t1) ¼ S0(jt) for jo1; (- - - - - -), S(t2) ¼ S(jt) for Predicted 5-year survival of ovarian cancer patients by age and graphical displays. Figure 2 illustrates this by giving predicted 5- The AFT model is commonly rewritten as being log-linear with year survival according to patient age and FIGO stage. Further respect to time, giving examples are demonstrated by Christensen (1987) based on theCox model, but can also be used when fitting fully parametric logðTÞ ¼ b0 þ b1x1 þ b2x2 þ    þ bpxp þ e models. In a previous analysis that involved some of the patients inthe present data, Clark et al (2001) produced a nomogram to where e is a measure of (residual) variability in the survival times.
summarise the impact of these and other covariates, and thus Thus, the survival times can be seen to be multiplied by a constant allows the reader to predict the median survival and the 2- and 5- effect under this model specification, and the exponentiated year survival probabilities for patients with given prognostic coefficients, exp(bi), are referred to as time ratios. A time ratio above 1 for the covariate implies that this ‘slows down', or The advantage of fitting a parametric survival model is that prolongs the time to the event, while a time ratio below 1 indicates predictions of the event survival, event hazard, mean and median that an earlier event is more likely.
survival times are readily available. For FIGO stages I – IV, the When the survival times follow a Weibull distribution, it can be median survival times are estimated to be 7.8, 4.0, 2.0 and 1.0 shown that the AFT and PH models are the same. However, the AFT family of models differs crucially from the PH model types interms of their interpretation of effect sizes as time ratios asopposed to hazard ratios.
ACCELERATED FAILURE TIME MODELS The survival times are usually assumed to follow a specific The accelerated failure time (AFT) model is a different type of distributional form in the AFT framework. Distributions such as model that may be used for the analysis of survival time data. For a the Log-Normal, Log-Logistic, Generalised Gamma and Weibull group of patients with covariates (x may be used to represent such survival data. Alternative methods 1, x2, y xp), the model is written mathematically as include the method of Buckley and James (1979), which isdiscussed by Stare et al (2000), and semiparametric AFT models, SðtÞ ¼ S0ðjtÞ in which the baseline survivor function is estimated nonparame-trically (see Wei, 1992, for an overview), but have not yet been where S0(t) is the baseline survivor function and j is an widely implemented in statistical software.
‘acceleration factor' that depends on the covariates according to As with the PH approach, other quantities such as projected survival probabilities may be derived. Also in keeping with PH models is the fact that AFT models make assumptions; the 1x1þb2x2 þ    þ bpxpÞg: appropriate choice of statistical distribution needs to be made, andalso the covariate effects are assumed to be constant and The principle here is that the effect of a covariate is to stretch or multiplicative on the timescale, that is, that the covariate impacts shrink the survival curve along the time axis by a constant relative on survival by a constant factor.
amount j. Figure 3 demonstrates this for the case of a singlecovariate (x1) with two levels, for example, x1 ¼ 0 for a placebo Parametric AFT models fitted to the lung cancer trial data group and x1 ¼ 1 for a new treatment group. The survivalprobabilities, S(t), for the placebo and new treatment groups are We use the non-small cell lung cancer dataset to illustrate the AFT S0(t) and S0(jt), respectively. The proportion of patients who are model, focusing on the relapse-free survival (i.e. , the time from event-free in the placebo group at any time point t1 is the same as diagnosis to the reappearance of cancer, with patients censored at the proportion of those who are event-free in the new treatment time of death if no recurrence had appeared). Again, we present group at a time t2 ¼ jt1. Figure 3 shows the cases where j41 and both the univariate and multivariate effect sizes in Table 3. The jo1, which represent situations where the length of survival is specific comparison of interest was the effect of adjuvant increased and decreased in the new treatment group compared (platinum-based) chemotherapy and radiotherapy compared with with the placebo, respectively.
radiotherapy alone. The unadjusted treatment effect may be British Journal of Cancer (2003) 89(3), 431 – 436 & 2003 Cancer Research UK Multivariate data analysisMJ Bradburn et al Time ratios from the generalised gamma AFT model for the lung cancer trial Univariate analysis Multivariate analysis Treatment (RT+CAP vs RT alone) Cell type (Sq vs non-Sq) Performance status (8 – 10 vs 5 – 7) Nodal involvement Age at diagnosis (/years) Gender (male vs female) Weight loss (X10 vs o10%) Race (white vs non-white) TR ¼ time ratio, CI ¼ confidence interval, RT ¼ radiotherapy, CAP ¼ cytoxan, doxorubicin and platinum-based chemotherapy, Sq ¼ squamous.
summarised by a time ratio of 1.91 (95% CI: 1.21 – 3.01; P ¼ 0.005), The strength of this method is in its simplicity: as the logrank which, having allowed for other covariates increased slightly to test is nonparametric, few distributional assumptions are made, 2.05. Therefore, we can conclude that the time to recurrence was and its interpretation is straightforward. Its main limitation is that significantly prolonged (approximately doubled) among patents it is only applicable when the covariate is categorical (or with given adjuvant chemotherapy in comparison with those who were continuous variables that have been arbitrarily categorised).
Further, this method does not perform well with several covariates, Again, we can derive model-based predictions: overall, patients as the number of individuals in each stratum quickly becomes too allocated to receive adjuvant chemotherapy had a predicted small to allow reasonable comparisons. In addition, it does not median survival time of approximately 16 months, as opposed to quantify the strength of effect of each variable, or even offer a P- 8 months among those treated with radiotherapy alone. Other value for factors other than the one of primary interest. This factors are also significant and would influence these times, but method is not generally regarded as a formal statistical model, but these are of less importance in the context of the comparative trial.
is of use where a very small number of covariates are to be We will return to this example in the next paper of this series.
considered, if only as an exploratory method of analysis.
Aalen's additive model WHICH MODEL SHOULD WE USE: PH OR AFT? Another approach to modelling the relationship between survival From a statistical viewpoint, an obvious way to choose between the and covariates is to assume that the covariates act additively on the two model types is to fit a type that is in keeping with the data. If hazard. Aalen's additive hazard model (Aalen, 1989) is one method the AFT model clearly fits the data better than the PH model, or that has been suggested for this, but its properties are rather unlike vice versa, this model may be adopted as being the more any other model described in this paper. The covariates are appropriate. However, in some cases, either type of model may assumed to impact additively upon a (unknown) baseline hazard, appear to fit the data adequately. In such instances, the choice of but the effects are not constrained to be constant. The impact is model may be influenced by other factors. For instance, if other therefore allowed to vary freely over time according to the studies of a similar nature had all used the Cox regression and underlying equation reported the results as hazard ratios, one may be tempted to followsuit to aid comparability. Against this, the parametric approach hðtÞ ¼ h0ðtÞ þ b1ðtÞx1 þ b2ðtÞx2 þ    þ bpðtÞxp offers more in the way of predictions, and the AFT formulationallows the derivation of a time ratio, which is arguably more where h(t) is the hazard, h0(t) is the baseline hazard and the bi(t) interpretable than a ratio of two hazards. As yet, however, AFT are coefficients, which may change in magnitude and even sign models are relatively unfamiliar and seen rarely in medical with time. Compare this with the Cox regression, where h0(t) is research papers (see Kay and Kinnersley, 2002).
also estimated nonparametrically, but the bi quantify the multi-plicative effect of covariate i on the hazard and are assumedconstant at all times.
As it is not straightforward to estimate h0(t) nonparametrically, the cumulative baseline hazard is used and the regression Stratified survival analysis coefficients that are actually estimated from the data are also the A more straightforward way to incorporate covariates into a cumulative (additional) hazard survival analysis is to use a stratified survival analysis. For example, suppose the covariate of primary interest is treatment, but we wish to control for the clinical stage of the tumour whenmaking the comparison. Here, the survival in each treatment group can be compared within each stage of disease (the ‘strata') The usual method of representing these effects is to graph them by the logrank or some other method, and the differences within against time. The further Bi(t) is from zero at time t, the greater the each stratum are then combined to give an overall comparison of effect the covariate has had on the hazard over the course of the treatments that has been adjusted for the stage.
study up to t. The values of bi(t), the absolute increase in hazard at & 2003 Cancer Research UK British Journal of Cancer (2003) 89(3), 431 – 436 Multivariate data analysis MJ Bradburn et al time t, are not actually observed, but their relative size may be relevant one, because it needs to be kept in mind that all the inferred from the slope of the line. These plots are sometimes models introduced here make certain distributional assumptions called Aalen plots, and they are also used to provide an informal of the survival times that will not always be met.
assessment of the adequacy of the proportional hazards assump- We have focused on the Cox model, the class of parametric PH tion in the Cox model, although Aalen considered its primary role models and AFT models as tools with which to analyse survival as an alternative model in its own right (Aalen, 1993).
time data. Other models exist (see, e.g., Collett (1994) for a more The flexibility of this approach is tempered by the lack of an practical demonstration of some alternatives and Bagdonavicˇius easy interpretation. The Bi(t) coefficients are not easy to under- and Nikulin (2001) for the theoretical background), but many are stand, and as they change repeatedly over time, can offer no single similar to, if not extensions of, the approaches we have discussed.
quantifiable effect size. Formal tests of statistically significant The use of the Cox model offers greater flexibility than parametric covariate effects may be carried out, but Aalen plots are essentially alternatives and, in particular, does not require the direct the only manner with which to interpret the effect sizes. These estimation of the baseline hazard function (i.e. it avoids the need reasons, together with the relative lack of statistical software, are to specify the distribution of the survival times). However, the probably the deciding factors in the relatively minimal use of assumption of proportional hazards is a crucial one that needs to Aalen's model.
be fulfilled for the results to be meaningful, and will not always besatisfied. Further, while the Cox PH model may be valid, other Classification trees and artificial neural networks parametric models will produce more precise estimates where thedistribution is specified correctly.
Two relatively recent developments are classification trees and A further concern is that the choice of covariates to include is artificial neural networks. These methods differ substantially in also far from simple. In the third paper of this series, we will their complexity and interpretation to the methods presented here consider ways to choose between the various model types, to and to each other. Both approaches are described in more detail in identify and assess the importance of covariates, and to verify that a later paper of this series.
the final model is adequate.
The principal strength of statistical models is their ability to assessseveral covariates simultaneously. The strengths of the stratified We wish to thank John Smyth for providing the ovarian cancer logrank test and other such methods are their obvious simplicity data, and Victoria Cornelius and Peter Sasieni for invaluable and the fact that they make fewer parametric assumptions of the comments on an earlier manuscript. Cancer Research UK data. Although these reasons are usually insufficient to suggest that supported all authors. Taane Clark holds a National Health Service the stratified method be used more widely, this second feature is a (UK) Research Training Fellowship.
Aalen OO (1989) A linear regression model for the analysis of life times.
Collett D (1994) Modelling Survival Data in Medical Research. London: Stat Med 8: 907 – 925 Chapman and Hall/CRC Aalen OO (1993) Further results on the non-parametric linear regression Cox DR (1972) Regression models and life tables (with discussion). J R model in survival analysis. Stat Med 12: 1569 – 1588 Statist Soc B 34: 187 – 220 Bagdonavicˇius V, Nikulin M (2001) Accelerated Life Models: Modeling and Kay R, Kinnersley N (2002) On the use of the accelerated failure time Statistical Analysis. London: Chapman & Hall/CRC model as an alternative to the proportional hazards model in the Buckley K, James I (1979) Linear regression with censored data. Biometrics treatment of time to event data: a case study in influenza. Drug Inf J 36: Christensen E (1987) Multivariate survival analysis using Cox's regression Stare J, Heinzl H, Harrell F (2000) On the use of Buckley and James least model. Hepatology 7: 1346 – 1358 squares regression for survival data. New approaches in applied Clark TG, Bradburn MJ, Love SB, Altman DG (2003) Survival analysis. Part statistics: Metodolosˇki zvezki 16 (http://mrvar.fdv.uni-lj.si/pub/mz/ I: basic concepts and first analyses. Br J Cancer 89: 232–238 Clark TG, Stewart ME, Altman DG, Gabra H, Smyth J (2001) A prognostic Wei LJ (1992) The accelerated failure time model: a useful alternative to the model for ovarian cancer. Br J Cancer 85: 944 – 952 Cox regression model in survival analysis. Stat Med 11: 1871 – 1879 British Journal of Cancer (2003) 89(3), 431 – 436 & 2003 Cancer Research UK

Source: http://biostat3.net/download/articles/BJC2003_Survival_analysis_part2.pdf