Modelling long term survival with non-proportional hazards
Perperoglou, A.
Citation
Perperoglou, A. (2006, October 18). Modelling long term survival with non-proportional hazards. Retrieved from
https://hdl.handle.net/1887/4918
Version:
Corrected Publisher’s Version
License:
Licence agreement concerning inclusion of doctoral thesis in the Institutional Repository of the University of
Leiden
Downloaded from:
https://hdl.handle.net/1887/4918
Chapter
5
Approaches in modelling long
term survival
Abstract
Several modelling techniques have been proposed for non-proportional haz-ards. In this work we consider different models which can be classified into three wide categories: models with time varying effects of the covariates, frailty models and cure rate models. We present those different extensions of the proportional hazards model on an application of 2433 breast cancer pa-tients with a long follow up. We comment on the differences and similarities among the models and evaluate their performance using survival and hazard plots, Brier scores and pseudo-observations.
5.1
Introduction
The Cox proportional hazards [24] model is one of the most popular methods to analyze survival data. In short follow up studies the assumption of a constant hazard ratio is very reasonable. However, in long follow up studies it is more appropriate to assume that time somehow influences the risk ratios, for instance in clinical trials where a treatment yields better survival probabilities, but with time the effect washes away. When the assumption of proportionality is violated then the results from a Cox model are not reliable and other modelling approaches should be considered instead.
varying effects of the coefficients. We apply such a model to the data set and consider Reduced-Rank models as a more general technique.
Another approach would be frailty models. We present the use of a gamma frailty model to the data, and work with the marginal hazard, which is the Burr model. We also present an extension of the latter, the relaxed Burr model. Finally, a cure mixture model which uses proportional hazards to model the survival function and logistic regression to identify cure rates is also used to fit the data.
A common practice for researchers is to report relative risks and the es-timated beta coefficients coming from a Cox model and over-emphasize their importance. The baseline hazard and as a consequence the survival functions get lost along the way and are not reported. However, the survival functions can reveal very important aspects of the data, that are of great clinical signifi-cance. In this work we fit different models, and report the coefficients, but the emphasis is given in plots of hazards and survival. In the extensions of the Cox model considered here, the estimated coefficients are only interpretable at time t = 0, that is, the start of the study. In this work the emphasis is given on the estimated survival probabilities under different modelling techniques rather than how to quantify the covariates effects.
This study reviews some of the extensions of the proportional hazards model and illustrates the fitting procedure in a real data set. Our aim is to investigate the interpretations of different approaches and illustrate the steps of building good survival models. To conclude which model is more appropriate under different cases we demonstrate techniques for evaluating each model. To assess the different modelling strategies, survival and hazard plots are given, we es-timate Brier scores andR2and build pseudo-observations for testing potential lack of fit. As a new method for evaluation we suggest the use of Brier scores on pseudo-observations.
5.2. IASO breast cancer data and proportional hazards analysis 0 5 10 15 20 0.0 0.2 0.4 0.6 0.8 1.0
Years since surgery
Survival
Figure 5.1: Survival function (Kaplan-Meier) of breast cancer patients. Dotted line presents the censoring distribution.
5.2
IASO breast cancer data and proportional hazards
analysis
Table 5.1: Results from Cox regression, age in years, tumor size (in mm), LN+ number of positive lymph nodes, tumor Grade (Bloom-Richardson), Chemotherapy (0 no, 1 yes), Radiotherapy (0 no, 1 yes) and Hormonal treat-ment (0 no, 1 yes)
coef se(coef) z p-value Age 0.017 0.003 4.804 1.6e-06 Tumor 0.011 0.002 4.539 5.6e-06 Ln+ 0.032 0.003 8.703 0.0e-00 Grade 0.373 0.070 5.299 1.2e-07 Chemo 0.006 0.099 0.060 9.5e-01 Radio -0.208 0.095 -2.179 2.9e-02 Horm -0.623 0.087 -7.113 1.1e-12
and axilla, and 1810 cases received hormonal treatment which consist of 20mg of Tamoxifen per day, for at least 5 years following the operation. The mean age of the patients was 56 years (from 23 up to 98) and 602 patient died within the follow up period. A Kaplan-Meier plot is presented in Figure 5.1 along with the survival function of the censoring. More information about patient characteristics can be found in [76].
5.3. Cox models with time varying effects of the covariates
log relative risk. That might not seem as an important increase in a patient risk, but patients with 10 positive lymph nodes have an increased risk of 3 units compared with those with negative lymph nodes. Tumor grading is the second most important factor (z-statistic 5.299), with an increase in risk of 0.373 mov-ing from patients with grade 1 to grade 2. Although gradmov-ing is a categorical variable it was entered as a continuous since its effect is linear. Each additional year of age is associated with 0.017 increase, while the risk difference for the youngest patient (23 years old) to the oldest (98 years old) is 1.275(75 ∗ 0.017), or if expressed in relative risk the oldest patient is 3.58 (exp(1.275)) times more likely to die out of breast cancer than the younger patient. The tumor size is also associated with an increased risk of 0.011 for each additional millimeter. Hormonal treatment is very important to the progress of the disease (z-statistic -7.113) as it decreases the risk of death by 0.623 units, while radiotherapy is also associated with a small decrease in risk. However, from a first look at the table chemotherapy seems to increase the risk of dying from the disease.
The Cox model is a very powerful and useful tool for survival analysis. However, if the assumption of proportionality does not hold, the results might be misleading. In a data set with large follow up period as this one, it is natural to assume that the assumption of proportionality will not hold. Although radiotherapy can be an important factor for survival, it is natural to assume that the effect of the treatment would wash away after some time, and the same might hold for other predictors entered the model. The proportionality assumption says that the risk will remain the same throughout all the follow up period, but in practice it is hard to believe that a treatment given to a patient will remain effective 15 years after. A simple test of proportionality is based on Schoenfeld residuals [39]. For the model fitted, the overallp-value of the test estimated with the function cox.zph was 1.56e-07, a clear indication of violation of the proportionality assumption.
5.3
Cox models with time varying effects of the covariates
An approach for modelling time varying effects of the covariates is to cut the follow up period into smaller time intervals and analyze each one of them using a simple Cox model. This approach is common in very large data sets where computational issues arise. Some special care should be taken for intervals with a small number of events. Time varying effects of the covariates can also be observed using Aalen’s additive regression model [2]. For a short discussion on how to deal with time varying effects see Therneau, pages 142-147 [92].
interaction term of the explanatory variables with a time function. This gives the Cox model with time varying coefficients:
h(t|X) = h0(t) exp(Xβ(t)) (5.1) whereh0(t) is the unspecified baseline hazard, X is a row vector of p covariates and β(t) = β0+ β1f1(t) models the interaction of the covariates with a time function, in this case f1(t). A popular choice for a time function is f (t) = log(t) or f (t) = log(t + 1), and other choices can be time itself, p(t), t2 or other polynomials. Although this extension is very easily implemented in standard software, it has the major drawback that the functional form of the time-varying effect is determined by the choice of the interaction term. A few studies tried to overcome this problem by suggesting fractional polynomials [13], interactions with B-splines [40] or models which use sequential analysis based on a factorization of the likelihood over the time intervals [33]. Another approach is to plot the Schoenfeld residuals for each covariate versus time, and let a spline fit indicate the functional form for that specific covariate and then fit a new model with interactions that mimic the functional form as it was indicated by the smoother. Although this approach is guaranteed to give very significant terms, it is data driven and thus the results can not be safely generalized in other cases.
Perperoglou et al. [77] addressed this problem of choice of time functions and parsimonious modelling by introducing reduced rank regression models in survival analysis. To explain the motivation behind this method start by re-writing equation (5.1) as
h(t|X) = h0(t) exp(XΘF0) (5.2) whereF is a row vector of q time functions including the constant, andΘ is a p × q matrix of estimable coefficients. For instance consider the model: h(t|X) = h0(t) exp(θ11x1+ θ12x1t + θ21x2+ θ22x2t) thenΘ = θ11 θ12 θ21 θ22
5.3. Cox models with time varying effects of the covariates
Table 5.2: Estimated coefficients from reduced-rank hazard regression, rank=3 fixed coef SE f2(t) f3(t) f4(t) f5(t) f6(t) Age 0.066 0.023 -0.073 -0.058 -0.039 -0.062 0.139 Tumor -0.009 0.015 0.036 0.020 0.015 0.015 0.019 Ln+ 0.007 0.028 0.015 0.032 0.021 0.046 0.199 Grade -0.325 0.449 1.671 0.713 0.525 0.303 0.027 Chemo 0.664 0.587 -0.600 -0.813 -0.524 -1.083 0.003 Radio -0.924 0.536 1.227 0.864 0.580 0.841 -1.978 Horm -1.696 0.620 0.358 1.520 0.927 2.472 -1.635
Table 5.3: Partial log-likelihood, number of estimated parameters, and AIC for different rank models.
log-likelihood parameters AIC rank=1 -4118.57 12 -4130.57 rank=2 -4092.22 22 -4114.22 rank=3 -4083.92 30 -4113.92 rank=4 -4081.40 36 -4117.40 rank=5 -4079.69 40 -4119.69 rank=6 -4079.62 42 -4121.62
be preferred. The motivation behind reduced rank models is to choose a matrix F with many flexible time functions, such as B-splines transformations of time, let the rank of the model be determined from a criterion such as AIC, and then let the model describe the time varying effects depending on the choice of the rank.
full rank model. All possible reduced rank models were fitted, from the very rigid rank=1 model to the very flexible full rank model. The rank=3 model was the one with the best AIC, (table 5.3) and it was chosen to analyze the data in the previous study. A rank=3 model has 30 degrees of freedom, much less than the full model that needs 42, however the differences to the full rank model are very subtle. The estimated coefficients under the reduced rank model are presented in table 5.2. The first column of the table corresponds to the fixed effects of the covariates, that is the effects right at the start of the follow up period. Note that these are quite different from the estimates given by the sim-ple Cox model, as they express different things. In a dynamic survival model the effects change with time and it is more interesting to study this covariate pattern through time, instead of looking at the estimated coefficients. In Fig-ure 5.2 the effects of covariates over time are presented. In the upper plot we illustrate the effects of the continuous covariates, while in the lower one, the effects of the categorical covariates. It can be seen that the effect of radiother-apy has a clear time varying behavior, as it decreases the risk of experiencing an event for the first year and then goes to zero after the second year, where its effect diminishes. The benefit of hormonal treatment is also clear for the first couple of years but then it reduces after the fourth year, although it remains an affective treatment for many years following the surgery. On the other hand having a chemotherapy looks like it has an adverse effect in the survival of the patients. However, a clear look at the data will reveal that the patient that had chemotherapy had already been of bad prognosis, and if someone wanted to study the effect of chemotherapy should compare group of similar patients. As it can be seen from the graphs, the rank=3 model is rather flexible for these data, as it approximates the full rank model very closely, and it looks like it overfits the data. If a rank=2 model was chosen instead the result would be more reasonable.
5.3. Cox models with time varying effects of the covariates 0 2 4 6 8 10 0.00 0.02 0.04 0.06
years since surgery
covariate effects age tumor Ln 0 2 4 6 8 10 −1.0 −0.5 0.0 0.5 1.0
years since surgery
covariate effects
grade chemo radio horm
0 2 4 6 8 10
0.00
0.02
0.04
0.06
years since surgery
covariate effects age tumor Ln 0 2 4 6 8 10 −1.0 −0.5 0.0 0.5 1.0
years since surgery
covariate effects
grade chemo radio horm
5.4. Frailty models
number of positive lymph nodes is more important for the disease outcome in the first years, and although the effect drops it always remains an important factor. The size of the tumor affects the survival of the patients for about 5 years, but after that it seems that is not so important for the survival, showing a similar pattern as the tumor grading.
Although the rank=3 model has a slightly better AIC than the rank=2 the difference is very small. The covariate effects, as illustrated in Figure 5.2, have a less reasonable behavior than the ones given under the rank=2 model. In [76] the preferred model was the rank=3 based on the AIC alone and the fact that the model gives a good approximation of the full rank model with less parameters. However, here we are not after an approximation of the full rank model, but we would rather prefer a more parsimonious model that avoids overfitting the data. The rank=2 model has slightly worse AIC (-4114.22 versus -4113.92 of the rank=3 model), however, plots of the covariate effects reveal more (biologically) reasonable patterns. Thus, we would conclude that the rank=2 model is the preferred one, and from now on we use this model for the analysis.
5.4
Frailty models
Frailty models are used in survival analysis to describe individual heterogeneity. The idea is that each person has his own characteristics and thus frailty, and the people who are more frail will tend to experience the event earlier. Since this individual heterogeneity is not observable, a random effect is entered into to the model, which acts multiplicatively on the hazard, and is in fact a random variable that describes excess risk. The frailty model can be written as:
h(t|X, Z) = Zh0(t) exp(Xβ) (5.3) A very popular choice for the random effect is to follow a gamma distribution. LetZ ∼Γ(1/ξ, 1/ξ) with E(Z) = 1 and var(Z) = ξ. Since the frailty is not observable one has to work with the marginal hazard that comes from model (5.3) by taking the expectation of the random effects. Let H0(t) denote the cumulative baseline hazard, then the marginal hazard at the population level is written as:
h(t|X) = h0(t) exp(Xβ) 1 + ξ H0(t) exp(Xβ)
Table 5.4: Estimated coefficients and standard errors from a gamma frailty (Burr) model and Relaxed Burr (RBurr) model.
Burr RBurr
coef se(coef) coef se(coef) Age 0.022 0.004 0.026 0.004 Tumor 0.018 0.003 0.020 0.004 Ln+ 0.045 0.006 0.064 0.012 Grade 0.595 0.094 0.597 0.105 Chemo 0.026 0.132 0.024 0.132 Radio -0.169 0.126 -0.229 0.134 Horm -0.927 0.118 -0.999 0.157
which is the Burr model [19]. Frailty models can be used also to account for the omission of an important covariate from the model, or in some cases to adjust for possible lack of fit and departure from proportionality.
In the breast cancer data, we had evidence that the assumption of propor-tionality is violated, and we would like to see how a more sophisticated model will fit the data. We fitted a gamma frailty model to the data with the same covariates as before, the results are shown in Table 5.4, with a significant frailty term (p-value < 0.001) and estimated variance of 1.53. With the inclusion of the random effect the estimated coefficients are larger than in the simple Cox model. Another main difference is that the hazards under the frailty model are not constant anymore, but they converge. The people that are more frail are supposed to experience an event rather early in time, and then the group of the remaining cases are the stronger individuals, with smaller and constant frailty for the rest of the time period. This setting is useful for describing the mech-anism behind the heterogeneity that is present in the data, and the departure from proportionality. However, in a study with such a long follow up, it will be reasonable to assume that the frailty does not remain constant throughout the follow up period, but changes with time. A reasonable assumption would be for weaker, more frail patients, to experience the event early in time, and then the survivors that do not experience an event for a long time to have a frailty that tends to 1 as time passes by. That assumption leads to the idea of time varying frailties.
5.4. Frailty models
and Henderson [78] as an approximation to a frailty model with time varying effects. Since estimating a model with time varying frailties can get complicated from a mathematical point of view, the authors suggested a generalization of the Burr model (5.4) that would be more flexible and able to approximate a model with autocorrelated frailties. The relaxed Burr model emerges if the link between the estimation of the covariate effects and the baseline hazard in (5.4) is loosened. This can be achieved by replacing the cumulative baseline hazard with a set of time functions F(t), and the variance of the random effect with unknown but estimable coefficients for the time functions δ. Then model (5.4) takes the form:
h(t|X) = h0(t) exp(Xβ)
1 + F(t|δ) exp(Xβ) (5.5) whereF(t|δ) can be any continuous non-negative function starting at F(0|δ)=0. We will always use a linear modelF(t|δ) = f (t)δ where f (t) can be a simple time function multiplied by an unknown but estimable coefficient δ. This gen-eralization of the Burr model allows for more flexible forms of hazard functions, depending on the functional form ofF(t|δ). For estimation of the model pa-rameters β and δ the Newton-Raphson algorithm can be used to maximize the partial likelihood: L = n
∏
i=1 " eXiβ(1 + F(t i|δ)eXiβ)−1 ∑j∈R(ti)e Xjβ(1 + F(ti|δ)eXjβ)−1) #di (5.6) wherediis the event indicator.We fitted the relaxed Burr model on the breast cancer data set, using a simple algorithm written in Fortran 77. We used cubic B-splines with no interior knot as time functions, which is actually just simple cubic polynomials passing through the origin. The choice of time functions was to be flexible enough to allow the model to reveal the appropriate functional form of the time varying pattern. As it can be seen in table (5.4) the estimates of the coefficients and their standard errors are very close to those given by the simple frailty model. However, the merit of the relaxed Burr model lies in its flexibility. In Figure 5.4 the estimated variance of the random effect times the cumulative baseline hazard ( ˆξ ˆH0(t)) is plotted along the estimated coefficients of the time functions of the relaxed Burr model ( ˆδ ˆF(t)). The graph reveals that for the first
0 5 10 15 0.00 0.05 0.10 0.15 0.20 time in years Burr R Burr
Figure 5.4: Plot of ˆξ ˆH0(t) for the Burr model and ˆδ ˆF(t) for the relaxed Burr.
particular ˆδ ˆF(t) reaches a maximum after a decade has elapsed and from then
on it decreases towards zero.
5.5. Cure models 0 5 10 15 −9 −8 −7 −6 −5 time in years log H0(t) PH Burr R Burr
Figure 5.5: Plot of log hazards for two different patient groups. shows hazards that converge after approximately 10 years. The relaxed Burr model gives almost identical hazards to the Burr model for the first 5 years, as is was already seen in graph (5.4) but from that point on the hazards start to diverge and come closer to the ones given by the simple Cox model. That reveals a very interesting pattern for the data. After the elapse of some time, the group of people that remain do not have constant frailties anymore, instead the variance of the frailty drops slowly to zero, and the model comes close to proportional hazards. That strengthen the idea of time varying frailties, and it could point to cure/mixed model, where in long follow up studies a proportion of cases could be consider cured.
5.5
Cure models
thus can be considered cured. In such cases mixture models can be employed to identify and model cured patients. In mixture models the overall survival of the patients consists of two parts, a survival functionSu(t|X) that models the survival of not-cured patients, denoted by the subscriptu, and a probability of a patient been cured π(X∗) which depends on some covariates X∗ and takes the logistic formlog[π(X∗)/{1 − π(X∗)}]. Overall the survival function at timet for patients with covariates X and X∗is given as:
S(t|X, X∗) = [1 − π(X∗)]Su(t|X) + π(X∗) (5.7) NowS(t|X, X∗) is the unconditional survival function of the entire population,
π(X∗) = P(C = 1|X∗) where X∗is a covariate vector that may include exactly
the same covariates asX, and C is an indicator of cured patients, i.e. C = 1 if the patient is cured andC = 0 otherwise. In that setting the logistic function logit(π) = η, where η = X∗β∗, is commonly used to describe the effect of
co-variates on the probability of being cured. For the survival of uncured patients many parametric forms have been suggested by several authors, ([32, 75]). However, staying close to our approach so far we choose the proportional haz-ards model to describe the survival of uncured patients. If the hazard function at timet for a patient with covariates X is given as h(t|X) = h0(t) exp(Xβ) then equation (5.7) is written as:
S(t|X, X∗) = [1 − π(X∗)] exp{−Hu0(t) exp(Xβ)} + π(X∗) (5.8) with unknown parameters for the model β, β∗ and Hu0. Note that the coef-ficients here have an opposite meaning, a positive β means that the risk of experiencing an event is higher, while a positive β∗means that the probability of being cured is higher. For estimation several approaches have been proposed [59, 90, 65]. However, we choose to estimate the model based on the method suggested by Peng [74] in which the M-step of the EM algorithm consists of fitting a proportional hazards model with an offset term and a logistic model with some fixed covariates. The software used here was develop by Peng as an S library [51] and is available for download at www.math.mun.ca/∼ypeng the authors website. We slightly modified this software to work under R 2.0.1.
5.6. Assessment of model fitting and evaluation of uncertainty
Table 5.5: Estimated coefficients, standard errors and Wald statistics (W) from a semiparametric cure model.
β se(β) W β∗ se(β∗) W Intercept 7.441 1.866 3.988 Age 0.017 0.003 4.746 -0.032 0.021 -1.549 Tumor 0.013 0.002 5.234 -0.036 0.013 -2.777 Ln+ 0.051 0.004 12.612 -0.173 0.021 -8.318 Grade 0.292 0.072 4.068 0.564 0.389 1.452 Chemo -0.063 0.100 -0.631 -0.633 0.986 -0.643 Radio -0.289 0.097 -2.975 0.821 0.622 1.320 Horm -0.750 0.087 -8.628 3.786 0.712 5.312
increased risk of being cured. From looking at the estimated coefficients of the logistic model (β∗) it can be concluded that patients with a small number of (or no) infected lymph nodes tend to have a higher probability of being cured (as indicated by the coefficient -0.173 and the large Wald statistic associated with that) with the next most important factor for cure being the hormonal treatment. Following that, is the effect of tumor size (large tumors are more likely to lead to relapse) while the rest of the covariates are associated with small Wald statistics. The estimated coefficients from the Cox model are very similar to the results coming from the simple model presented in table (5.1) with the exception of the chemotherapy effect, which in this model shows that better chances of survival are associated with the treatment. Although the Cox parameter changed sign, it is compensated by the negative coefficient estimate in the cure part of the model, making the net effect of chemotherapy positive.
5.6
Assessment of model fitting and evaluation of
uncer-tainty
Use of pseudo-observations
time pointt the pseudo-observations are defined as: ˆ
P.obsi(t) = nKM(t) − (n − 1)ˆ KMˆ−i(t) (5.9) where KM is the Kaplan-Meier survival estimate at timet and KM−i(t) is the estimate at the same time, leaving individuali out. Each case i is assigned with a pseudo-observation at timet and these values can be as the response in a generalized linear model with link functiong(.)
g(P.obs) = βXˆ
The regression coefficients can now be estimated using generalized estimating equations [62] U(β) = n
∑
i=1 Ui(β) = n∑
i=1 ( ∂ ∂βg −1(βX i)Vi−1(P.obsˆ i− g−1(βXi)) = 0 withVithe covariance matrix for P.obs. A natural link function to use will be the logit function. Once the pseudo-observations are computed, estimation of the coefficients and the covariance matrix can be done using standard statis-tical software that fit generalized linear models. One can regress the pseudo-observations on the covariates and use the estimations that come out a logistic regression model, instead of a Cox model since using this approach there is no worry about censoring.Brier Scores
A measure of uncertainty in survival models was proposed by Graf et al [38] based on Brier scores. Originally Brier scores emerged for judging the inac-curacy of probabilistic weather forecasts based on quadratic scores but have adopted for survival models by Graf in her Ph.D thesis [37] and used as a method of incorporating model uncertainty [11]. A Brier score measures aver-age discrepancies between the true disease outcome and the predictive values from the model. The expected Brier score can be defined as:
5.7. Comparison
up to timet, i.e. ti> t are weighted with 1/G(t). For events before time t, ti< t the weights are ui= 1/G(ti). The average Brier score can be seen as a mean square error of prediction, thus the lower the score the better the model. Based on these scores aR2measure, which is similar to the one used in linear models, can be defined asR2= 1 − Bs(t)m/Bs(t)K, where the indexm denotes the fitted models andK the Kaplan-Meier.
When comparing the Brier scores for a set of values, one problem that has to be solved it the censoring. Using the weightsui defined above is one way to overcome that problem, while an alternative will be to use the pseudo-observations which try to estimateI(Ti> t). Here we combine the use of Brier scores and pseudo-observations, in order to obtain a modified score that would not be heave influenced on the amount of censoring. For that purpose define the Brier pseudo-observation score:
Bp.s(t) =1 n n
∑
i=1 {[P.O(t) − ˆS(t|Xi)]2}where P.O(t) are the pseudo-observations evaluated at time t. When using the pseudo-observations to estimate the Brier scores a positive bias should be expected because of the extra variation. However, since the bias is the same for all models, these scores can also be used for comparison. Nevertheless, it is harder to turn these scores into aR2measure.
5.7
Comparison
Survival plots
0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 time in years Survival PH RR B RB Cure 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 time in years Survival 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 time in years Survival 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 time in years Survival
5.7. Comparison
tumor grading 3. Finally, the last group consists of the worst prognosis patients, being 80 years old with a tumor of 50mm, 15 positive lymph nodes and tumor grading 3.
In Figure 5.6 the survival functions for the four different groups under five different models are presented. The solid line is the survival function out of a Cox proportional hazards model, and it serves as a reference line to the other fitted models. For the first group of patients, the ones with the best prognosis, the relaxed Burr model gives the most optimistic survival function, although the differences among the models are subtle. In the second group of patients, hardly any difference can be seen up to five years of follow up among the rank=2 and the proportional hazards model. From the sixth year and on, the reduced rank model is giving larger probabilities of survival that the rest. As we have commented earlier, the effects of the covariates are weakened under this model after some time has elapsed, thus it is reasonable the effects of the covariates here to weaken after 6 years of follow up and the survival function of group 2 patients to be higher than the rest. The same holds for the cure model, after six years some of the long term survivors are considered cured by the model, thus the probability of survival is higher. The same pattern of the reduced rank model is also present for group 3 and group 4 patients. The differences among the Burr model and relaxed Burr model are hardly visible for up to 6 years in all patient groups, while from the seventh year and on are very small. In the third and fourth group, the Burr model gives slightly larger probabilities of survival at the last couple of years shown in the graph.
Assessment of predictive value
Table 5.6: Brier scores for different models, best model in bold. Time (t) in months.
t KM PH Burr RBurr Rank=2 Cure 24 0.056 0.053 0.053 0.053 0.051 0.053 36 0.098 0.089 0.088 0.089 0.088 0.089 48 0.124 0.113 0.113 0.113 0.113 0.114 60 0.146 0.131 0.130 0.130 0.130 0.132 72 0.166 0.149 0.148 0.147 0.149 0.150 84 0.186 0.166 0.165 0.164 0.167 0.168 96 0.198 0.176 0.175 0.174 0.179 0.179 120 0.218 0.193 0.192 0.190 0.197 0.197
Table 5.7: R2 measure for different models, best model in bold. Time (t) in months.
5.7. Comparison
Table 5.8: Brier scores for different models estimated on pseudo-observations, best model in bold. Time (t) in months.
t KM PH Burr RBurr Rank=2 Cure 24 0.0598 0.0563 0.0558 0.0561 0.0554 0.0563 36 0.1102 0.1013 0.1004 0.1006 0.1000 0.1016 48 0.1472 0.1359 0.1352 0.1352 0.1350 0.1364 60 0.1809 0.1657 0.1648 0.1646 0.1649 0.1675 72 0.2215 0.2041 0.2032 0.2027 0.2041 0.2064 84 0.2690 0.2478 0.2471 0.2463 0.2489 0.2507 96 0.3032 0.2797 0.2788 0.2780 0.2814 0.2831 120 0.3858 0.3602 0.3592 0.3581 0.3633 0.3654
the proportionality assumption does not hold in the data. We will comment on that in the discussion. An important aspect to consider here is how much the Brier score is affected by the amount of censoring going on in the data. As it can be seen from Figure 5.1 on the fifth year more than 50% of the patients are already censored. That means that the Brier score evaluated at timet = 5 years will have 1083 zeros in the vector of u’s. In order to check whether censoring leads to biased Brier scores we can use pseudo-observations instead of the indicator functionI(T > t) in equation (5.10).
We computed the Brier pseudo-observation scores at the same time points as before. The results are presented in table (5.8). The results show that the reduced rank model is the best up to the fifth year, and after that the relaxed Burr model is still the preferred one.
