On a general structure for hazard-based regression models: an application to population-based cancer research

05/21/2018 ∙ by Francisco J. Rubio, et al. ∙ London School of Hygiene & Tropical Medicine 0

The proportional hazards model represents the most commonly assumed hazard structure when analysing time to event data using regression models. We study a general hazard structure which contains, as particular cases, proportional hazards, accelerated hazards, and accelerated failure time structures, as well as combinations of these. We propose an approach to apply these different hazard structures, based on a flexible parametric distribution (Exponentiated Weibull) for the baseline hazard. This distribution allows us to cover the basic hazard shapes of interest in practice: constant, bathtub, increasing, decreasing, and unimodal. In an extensive simulation study, we evaluate our approach in the context of excess hazard modelling, which is the main quantity of interest in descriptive cancer epidemiology. This study exhibits good inferential properties of the proposed model, as well as good performance when using the Akaike Information Criterion for selecting the hazard structure. An application on lung cancer data illustrates the usefulness of the proposed model.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The analysis of time-to-event data has been dominated by the use of the semi-parametric Cox model during the last decades. While extensions to remove the assumption of proportional hazards (PH) were already discussed in Cox’s original paper, a lot of work has been devoted to improve the flexibility of hazard-based regression models using flexible functions for both the baseline hazard and the inclusion of time-dependent parameters, mainly using splines or fractional polynomials (Durrleman and Simon, 1989; Sleeper and Harrington, 1990; Hess, 1994; Abrahamowicz et al., 1996; Royston and Sauerbrei, 2008). However, the general structure of those models remained the same: the hazard function was expressed with a baseline hazard multiplied by the exponential of a flexible function of time and covariates : . Alternative hazard structures that directly account for time-dependent effects of covariates have been proposed; one of the earliest alternative to the PH model is the accelerated failure time (AFT) model, in which the variables have a direct effect on the time to event, in contrast to the PH model where the covariates affect the hazard function (Kalbfleisch and Prentice, 2011). More recently, in a series of papers (Chen and Wang, 2000b; Chen and Jewell, 2001; Chen et al., 2003, 2014), Y. Chen and co-authors proposed and studied semiparametric accelerated hazards (AH) models, where the covariates have a time-scaling effect on the hazard function, thus allowing for a time-varying effect.

In cancer epidemiology using population-based registry data, the quantity of interest is usually the excess hazard (of death) instead of the overall hazard (Hakulinen and Tenkanen, 1987; Esteve et al., 1990; Remontet et al., 2007; Mariotto et al., 2014). The basic idea behind excess hazard models consists of decomposing the hazard function of death associated with an individual, , as the sum of an excess hazard, , and the general population hazard . The general population hazard is supposed to correctly reflect the other-causes hazard of death in our population of interest, assuming that the contribution in the general population hazard of the disease of interest (e.g. a specific cancer type) is small compared to all others. This excess hazard is interpreted as the hazard that could be due, directly or indirectly, to the cancer under study. Formally, this can be written as:

(1)

where and

are vectors of covariates,

typically corresponding to a subset of covariates of , “age” is the age at diagnosis and “year” is the year of diagnosis (thus “” and “” are respectively the age and the year at time after diagnosis). The population hazard is typically obtained from national life tables based on the available sociodemographic characteristics ( in addition to age and year). A survival function derived from the excess hazard is called the net survival

. Net survival represents a useful way of reporting the probability of survival of cancer patients since this allows for a fairer comparison of survival rates between different populations or countries

(Gamel and Vogel, 2001; Perme et al., 2012; Danieli et al., 2012; Perme et al., 2016)

, and a non-parametric estimator of net survival has been proposed recently

(Perme et al., 2012). Although the interest in cancer epidemiology is on a slightly different quantity (the excess hazard instead of the overall hazard), the methodological developments of regression models have followed the same path than those described above (Bolard et al., 2002; Giorgi et al., 2003; Nelson et al., 2007; Remontet et al., 2007; Mahboubi et al., 2011; Charvat et al., 2016) in terms of the hazard structure adopted.

Our aim is to provide a valuable supplement in the available toolbox for analysing survival data, and which is applicable for both overall and excess hazard regression models. The proposed approach builds on top of two recent developments: (i) the general parametrisation of hazard functions (Chen and Wang, 2000b; Chen and Jewell, 2001), combined with (ii) the use of a flexible parametric distribution for modelling times to event, the exponentiated Weibull distribution (Mudholkar et al., 1996). The paper is organised as follows. In Section 2, we describe the general hazard structure, discuss how the accelerated failure time model, the proportional hazards model, and the accelerated hazards model are nested within this general structure, and discuss the parameter interpretations for these models. In this section, we also introduce the exponentiated Weibull distribution, which will be used to model the “baseline” hazard in the general hazard structure (and the corresponding sub-models) based on its flexibility and ease of implementation. We also discuss maximum likelihood estimation and inference associated with these models. In Section 3, we present an extensive simulation study where we illustrate the inferential properties of the proposed models. In Section 4, we present a data example from lung cancer epidemiology. Here, we illustrate the usefulness of the proposed models and compare them against appropriate competitor approaches. We conclude with a general discussion and point to possible extensions in Section 5. Additional simulations and results are provided in the Supplementary Material.

2 Methods

2.1 The different model structures

We considered the following excess hazard structures (see Chen and Wang (2000b), Chen and Jewell (2001) and Chen et al. (2003) for extensive discussion). We express the different structures below via the hazard function and the cumulative hazard function , respectively, according to time and a vector of covariates . We assume that the vector of covariates does not contain an intercept, in order to avoid identifiability issues. The survival function can be obtained from the well-known relationship . The vector denotes the unknown regression parameters.

  1. Proportional hazards model (PH).

    (2)
  2. Accelerated hazards model (AH).

    (3)
  3. Hybrid hazards model (HH).

    (4)

    where , and .

  4. Accelerated failure time model (AFT).

    (5)
  5. General hazards model (GH).

    (6)

The assumptions behind each of these models are different in nature. The basic idea is to include effects that affect the level of the hazard (time-fixed effects) and the time scale (time-dependent effects) separately, as follows. In the PH model (2), a unit change in a covariate value have a multiplicative effect on the hazard, thus leading to a level change on the y-axis of the hazard. In the AH model (3), the effect of a unit change in a covariate affects the time scale of the baseline hazard (x-axis), thus assuming that there is a time-dependent effect of each covariate. In other words, could be seen as a factor of how much more (or less) time is needed to reach the same hazard level, as compared to baseline, when the covariate is increased by one unit. In addition, the AH model does not assume that the parameters affect the hazard immediately at time , but “gradually”, which is not the case with the other hazard structures Chen et al. (2014). While a “gradual” effect can be useful when estimating a treatment effect, it may not be justified for other variables; it highlights the usefulness of the HH structure (4). Indeed, the HH relaxes the assumption of the AH model by allowing some variables to have a proportional hazards effect rather than time-dependent “gradual” effects. In the AFT model (5

), the effect is on the survival time directly as it can be, in fact, formulated as a log-linear regression model on survival times

(Kalbfleisch and Prentice, 2011). In such models, the estimated regression parameter for a unit change in one covariate accelerates or decelerates the event time, thus leading to a time-dependent effect. Notice that the AFT, PH, and AH models coincide for the case when the baseline hazard is Weibull (Chen and Jewell, 2001). The GH model (6) represents a general hazard structure that contains, as particular cases, the PH, AH, HH, and AFT models. More specifically, if , then ; if , then ; and if , then .

2.2 The Exponentiated Weibull distribution

We propose modelling the baseline hazard

using the Exponentiated Weibull (EW) distribution. The Exponentiated Weibull distribution is simply obtained by exponentiating the Weibull cumulative distribution function to an unspecified positive power

(Mudholkar et al., 1995). This simple transformation adds a second shape parameter that, interestingly, induces considerable flexibility to the hazard function. The hazard function of the Exponentiated Weibull distribution can capture several basic shapes: constant, increasing, decreasing, bathtub, and unimodal (See Figure 1), making it appealing for survival models (Cox and Matheson, 2014). This distribution has been recently used in the context of AFT models Khan (2017), while alternative families of flexible parametric AFT models have also been studied in Rubio and Genton (2016); Rubio and Hong (2016). The Exponentiated Weibull probability density and cumulative distribution functions with scale, shape, and power parameters are given, respectively, by:

(7)

where . This distribution reduces to the Weibull distribution for . The corresponding hazard function is obtained, by definition, as . The EW distribution is identifiable (See Proposition 1 from the Appendix), and general results about the identifiability of the different hazard-based models (2)–(6) are presented in Chen and Jewell (2001). It was shown that the GH model is identifiable except for the case when the baseline hazard is Weibull Chen and Jewell (2001). This is not an issue since, as already mentioned, it corresponds to the case when the AFT, PH, and AH models coincide. Moreover, it has also been shown that the AFT and AH models are not identifiable when the baseline hazard is flat Chen and Jewell (2001) (i.e. exponential), which corresponds to a case when the shape of the baseline hazard does not change with time, a case of little interest in practice.

2.3 Parameter interpretation

Some clarifications on the interpretation of the parameters estimated from the AH and GH models seem appropriate as these models are not as well known as PH and AFT models. Notice that these interpretations directly translate to the HH model as this model is a special case of the GH model.

We start by interpreting the parameters for the AH model (3

), as this will facilitate the interpretation for the GH model. The parameter interpretation depends on the shape of the baseline hazard, which we classify here as monotone (increasing/decreasing) or not (bathtub/unimodal).

For monotonic baseline hazard, a positive value of for one unit change in the covariate means that has a harmful (beneficial) effect if the baseline hazard is increasing (decreasing) (see panels (a) and (b) in Figures S3–S4. A negative value of means that has a beneficial (harmful) effect if the baseline hazard is increasing (decreasing) (see panels (a) and (b) in Figures S1–S2). In other words, a positive value of accelerates the progression of the hazard, which is beneficial for the patients if the hazard decreases and harmful if it increases. On the contrary, a negative value of decelerates the progression of the hazard, which is beneficial for the patients if the hazard increases and harmful if it decreases.

If the shape of the baseline hazard is unimodal (or bathtub), a positive value of accelerates the evolution of the hazard, thus the maximum (minimum) will be reached sooner (see panels (c) and (d) in Figures S3–S4). A negative value of decelerates the evolution of the hazard, thus the maximum (minimum) will be reached later (see panels (c) and (d) in Figures S1–S2).

From the AH model, it is worth noticing that the general shape will not change (a unimodal shape will remain unimodal) and that the peak/minimum (if any) reached by the hazards defined by different subgroups will be at the same level. On the other hand, when using the GH model (6) the parameter is directly multiplying the time , thus re-scaling the timescale (accelerating or decelerating, i.e. they play a role on the x-axis, changing the pace of the hazard’s progression), while the parameter is modifying the level of the hazard (role on the y-axis, changing the magnitude of the hazard, Figure 1). This can be clearly seen in the right panels of Figure 1 for the unimodal shape: in panels (a) and (b) is negative, thus the peak is reached later for patients’ group with (dot-dashed grey lines) compared to the baseline group (solid grey lines), while the parameter changes the magnitude of the hazard (thus the level of the peak). In panels (c) and (d), is positive, thus the peak is reached sooner for patients’ group with compared to the reference group, while the peak level is changed according to . The same interpretation applies for bathtub hazards (black lines in the right panels in Figure 1), and for monotonic hazards (left panels in Figure 1), even though for these later the interplay of both parameters and is less obvious to see on the graphs.

The AH and GH models allow a crossover of the hazards (and also of the survival functions). In some cases, this advantage may lead to difficulties in interpreting clearly the parameters (Chen and Wang, 2000a). Thus, plotting the hazards according to different covariate patterns is recommended to help clarifying the time-to-event process compared to reporting only the survival functions. Indeed, the survival probability, being a cumulative measure, does not help visualising particular features of the instantaneous process (Hess and Levin, 2014).

Figure 1: Hazard shapes obtained with different values of the parameters defining the Exponentiated Weibull distribution , and different values of the regression coefficients and for a binary covariate . Monotonic hazard (increasing in black and decreasing in grey) are displayed on the left column, and bathtub or unimodal (black or grey, respectively) on the right column, where solid lines represent the baseline hazard (e.g. unexposed ) and dot-dashed lines represent the hazard for the exposed ones (). The values of are: , , , for monotonic increasing, monotonic decreasing, bathtub and unimodal shapes, respectively.

2.4 The likelihood function

Let , be a sample of times to event from a population of cancer patients , with covariates , and vital status indicators (1-death, 0-censored). Let also be the corresponding population hazard rates obtained from the national life tables based on the variables , , where represents the age at diagnosis and the year of diagnosis of patient . The likelihood function of the full vector of parameters is then given by

Notice that we have removed the term from the likelihood as it does not depend on the parameters .

In order to obtain confidence intervals for the parameters, we reparameterise the baseline EW hazard in terms of

. Appealing to the consistency and asymptotic normality of the maximum likelihood estimators (MLEs) of the EW distribution, and the availability of large samples in the context of cancer epidemiology, we propose the use of asymptotic confidence intervals of the type , where is the negative of the Hessian matrix of the log-likelihood function under the appropriate parameterisation, and is the confidence level.

Confidence intervals for the net survival curve at specific time-points are obtained using a simulation-based algorithm (Mandel, 2013)

. The idea is to simulate from the asymptotic (multivariate normal) distribution of the parameters in order to obtain a Monte Carlo sample of the net survival at specific time-points, which is used to construct the corresponding confidence intervals.

To measure relative goodness-of-fit of a hazard-based model structure among the ones detailed previously, we employ the Akaike Information Criterion (AIC):

where is the number of parameters (the dimension of ), and is the MLE of .

3 Simulation study

In this section, we present an extensive simulation study where we illustrate the good frequentist properties of the proposed models as well as the ability of the AIC criterion to select models that properly capture the underlying hazard structure. The parameter values are chosen in order to produce scenarios that resemble cancer population studies concerning an aggressive type of cancer (relatively low 5-year survival), such as lung cancer.

3.1 Data generation and simulations designs

We simulated data sets of size , assuming the additive hazard decomposition given in (1

). The variable “age” was simulated using a mixture of uniform distributions with

probability on , probability on , and probability on

years old. The binary variables “sex” and “W” were both simulated from a binomial distribution with probability

(the variable “W” could be viewed as “treatment” or “comorbidity”). In all the scenarios, we simulated the “other-causes” time to event using the UK life tables based on “age” and “sex”, assuming all patients were diagnosed on the same year. The time to event from the excess hazard (cancer event time) was generated using the inverse transform method Ross (2006), and assuming effects of the 3 variables “age”, “sex” and “W”. We assumed either (i) only administrative censoring at years, which induced approximately

censoring in all cases, or (ii) an additional independent random censoring (drop-out) using an exponential distribution with rate parameter

. In this latter case, we choose values for to induce around censoring. We considered 6 generating mechanisms for the excess hazard, all of them with EW baseline hazards: (i) Proportional Hazards (Scenario PH), (ii) Accelerated Failure time (Scenario AFT), (iii) Accelerated Hazard (Scenario AH), (iv) Hybrid Hazard (Scenario HH), (v) General Hazard (Scenario GH), and (vi) a hybrid hazard model where the hazards and the survival curves associated to the variable “W” cross (Scenario CH). The values of the model parameters used for each scenario are summarised in Table S1 from the Appendix. Moreover, tables S2–S7 from the Appendix present the 1-year and 5-year net survival implied by the different combinations of hazard structures and covariate values.

3.2 Analysis of the simulated data

To assess the frequentist properties of the models in all simulation scenarios, we fitted the GH model with EW baseline hazard (GHEW), as compared to the corresponding true generating model. Additionally, to measure the ability of the AIC to select the true generating model, we also fitted the remaining models in each scenario. In all cases, we used the true parameter values of the generating model as initial points in the optimisation step because our aim was more to study the estimator’s properties rather than the properties of the optimisation process. However, as the interest of analysts is also on the overall properties (i.e. including the optimisation process properties in the full process of estimation), we present in the Appendix a thorough study where we present three alternative “automatic” choices for the initial points which we will later use in the real data example (see section 5 in the Appendix). We present results about the convergence using these initial points as well as the resulting estimators, which are virtually the same as those obtained with the chosen initial points, as expected, with the expense of a higher computing time.

We performed the analysis using R software. The optimisation step was conducted using the R commands ‘nlminb()’ and ‘optim()’, while the Hessian matrix used in the construction of the asymptotic confidence intervals are approximated using the command ‘hessian()’ from the R package ‘numDeriv’, which represents an efficient method to calculate the Hessian matrix. The cases where the command ‘hessian()’ produced ‘Inf’ or ‘NaN’ values were excluded as this merely represents a numerical problem associated to numerical differentiation, which mainly affects the most difficult cases (i.e. with and censoring).

3.3 Measures of performance

We report the mean and median of the MLEs for the corresponding models, as well as the empirical standard deviation, the mean (estimated) standard error, the root-mean-square error, and coverage proportions of asymptotic confidence intervals. In addition, we report the proportion of times the different fitted models are selected using AIC. The excess hazard functions associated with the best models selected using AIC are plotted in Figures S5–S14, for different covariate patterns, and compared to the true generating model.

3.4 Simulation results

The results are presented in Tables S8–S29 in the Appendix. For illustrative purposes, we present in Tables 12 and Figure 2 the results for Scenario GH with . In this scenario, we observe very good performance of our approach, with very small bias and a coverage close to the nominal value of 95% (Table 1 and Figure 2). Regarding performances of model selection using the AIC, the true model was selected in around 1/3 of the simulated datasets in the worst case (sample size of 1000 and 50% of censoring), but this proportion increases to 90% or more in situations with a larger sample size (Table 2). More generally, from this simulation study, we can conclude that the MLEs are close to the true values for moderate samples, even in the case of a high censoring rate (Tables S8, S12, S16, S20, S24, S26) and that the coverage is usually quite close to . As expected, the RMSE decreases as the sample size increases (or with identical sample size but lower censoring percentage). For

, we observe low bias and variance of the MLEs, and that the AIC selects the correct model with a high proportion. For

, the AIC selects the true model in more than of cases for all scenarios and whatever the censoring level. Interestingly, in the GH scenario, as long as the sample size was equal to 5000 or larger, the AIC selected the true model in or more samples, regardless of the censoring level. Moreover, even in cases where the incorrect model is selected using AIC, we see that the baseline hazard is close to the true generating model, reflecting that the AIC selects the model closest to the true generating model (see figures S5–S14).

Model Parameter MMLE mMLE ESD Mean Std Error RMSE Coverage
censoring
GHEW (1.75) 1.747 1.726 0.405 0.409 0.405 0.950
(0.6) 0.601 0.600 0.063 0.064 0.063 0.946
(2.5) 2.576 2.523 0.461 0.459 0.467 0.947
(0.1) 0.100 0.100 0.013 0.013 0.013 0.961
(0.1) 0.112 0.109 0.233 0.227 0.233 0.949
(0.1) 0.099 0.104 0.230 0.228 0.230 0.956
(0.05) 0.050 0.050 0.003 0.003 0.003 0.958
(0.2) 0.201 0.202 0.041 0.042 0.041 0.951
(0.25) 0.251 0.251 0.043 0.042 0.043 0.943
censoring
GHEW (1.75) 1.741 1.736 0.457 0.465 0.457 0.957
(0.6) 0.600 0.600 0.074 0.076 0.074 0.948
(2.5) 2.601 2.508 0.556 0.539 0.565 0.963
(0.1) 0.102 0.101 0.016 0.016 0.016 0.959
(0.1) 0.090 0.095 0.267 0.266 0.267 0.957
(0.1) 0.092 0.093 0.269 0.268 0.269 0.948
(0.05) 0.050 0.050 0.004 0.004 0.004 0.945
(0.2) 0.201 0.200 0.045 0.047 0.045 0.956
(0.25) 0.252 0.251 0.049 0.048 0.049 0.949
Table 1: Simulation results for the scenario GH with , , , and . Mean of the MLEs (MMLE), median of the MLEs (mMLE), empirical standard deviation (ESD), mean (estimated) standard error, root-mean-square error (RMSE), and coverage proportions (Coverage).
Model censoring censoring
PHEW 1.6 2.7
AHEW 0 0
AFTEW 56.5 64.5
GHEW 41.9 32.8
PHEW 0 0
AHEW 0 0
AFTEW 3.2 10
GHEW 96.8 90
PHEW 0 0
AHEW 0 0
AFTEW 0 0.5
GHEW 100 99.95
Table 2: Percentage of models selected with AIC in the scenario GH with , , and .
Figure 2: Mean of the best fitted hazards in terms of AIC (dashed lines), compared to the true generating hazard (continuous lines), and 1000 sample-specific fitted hazards (grey lines) in the scenario GH for and censoring. Panels from left to right correspond to covariate values (age, sex, comorbidity)=, respectively. The dashed and continuous lines are virtually the same.

4 Real data example

To illustrate the new proposed models, we analysed a dataset obtained from population-based national cancer registry of lung cancer patients diagnosed in 2012 in the United-Kingdom. We linked these data to administrative data (Hospital Episode Statistics -HES- and Lung Cancer Audit data -LUCADA-) and applied specific algorithms to derive information on stage at diagnosis and presence of comorbidities at the time of diagnosis Benitez-Majano et al. (2016); Maringe et al. (2017). To retrieve information on comorbidity, we used a 6-year period up to 6 months before diagnosis where we checked for the presence of any of 18 comorbidities that are used to define the Charlson Comorbidity Index (CCI), in addition to obesity Maringe et al. (2017). The information was then dichotomised into 2 categories: “no comorbidity” vs. “at least 1 comorbidity (comorbidity indicator = 1)” in our illustrative example. We measured deprivation using the Income Domain from the 2010 England Indices of Multiple Deprivation, defined at the Lower Super Output Area level (mean population ). The Income Domain measures the proportion of the population in an area experiencing deprivation related to low income, and ranges from to (https://www.gov.uk/government/statistics/english-indices-of-deprivation-2010). Follow-up was assessed on the 31st of December 2015, at which time patients alived were censored (so the maximum follow-up was 4 years). We restricted our analysis to women with no missing data, and applied the PH model (2), the AH model (3), the AFT model (5) and the GH model (6), with an Exponentiated Weibull distribution for the baseline hazard. The variables included in the models are ‘agediagc’ (centred age at diagnosis), ‘Istage2’, ‘Istage3’, and ‘Istage4’ (stage at diagnosis), ‘INCOME_SCORE_2015c’ (centred Income Domain), ‘comorbidity’ (presence of comorbidity). The time dependent effects are indicated in Table 3 with the subindex ‘’.

We observed patients with complete cases among which died before the 31st of December 2015. The median follow-up among patients censored was years. The , and quantiles of the patients’ age at diagnosis was , , while the mean was . Among the patients, were Stage I, were Stage II, were Stage III, and were Stage IV. Finally, patients were classified with comorbidity indicator 1. Results of this analysis are presented in Table 3. From this table, we observe that the GHEW is clearly favoured by the AIC (followed by the AFTEW model), thus suggesting the need for including time-dependent as well as proportional effects. The signs of all the estimates are positive in this model. This implies that an increase of one unit in the value of any of the covariates leads to an acceleration of the time to reach the maximum of the hazard as well as an increase of the maximum of the hazard. The magnitude of such acceleration is given by the value of the corresponding estimated parameter, and we noticed that these two effects are different for each variable thus explaining the better fit of the GH model compared to the AFT model. We have compared the GHEW model against alternative models with fewer covariates but they were not favoured by the AIC. For comparison, we have also used the ‘mexhaz’ R package Charvat et al. (2016); Charvat and Belot (2017) to produce models with time-dependent effects (using the command nph) for some or all the variables and we have found that none of these models provided a better fit than the GHEW in terms of AIC (with ‘mexhaz’, the model with the lowest AIC () was the model assuming time-dependent effects for all variables). Figures 3a and 3b show the shapes of the hazard functions associated to the different models for patients aged 70 years at diagnosis, an income score value of 0.15 and with either a stage II or a stage IV cancer, and without comorbidity (panel (a)) or with comorbidity (panel (b)). These figures also show the piecewise excess hazard estimated separately on each of the 4 groups defined by the following values of covariate: , , Stage II or Stage IV, and . Those piecewise hazards represent a way to check the quality of the fit of the different models (Remontet et al., 2007). These figures suggest a good fit of the GH model overall, while the excess hazards obtained from the AH model are always over-estimated after 6 months for the group of stage IV patients. For the group of stage IV patients with comorbidity (panel (b)), the PH model does not capture the high excess hazard just after the diagnosis and the sharp decrease that follows.

Table 4 presents the net survival at years, for the total population (Comorbidity = 0, 1) and for a subgroup (age-group 55-65, Stage I, and Comorbidity = 0, 1, ), using the GHEW model and the Pohar-Perme nonparametric estimator Perme et al. (2012). The confidence intervals for the net survival in the GHEW model were obtained using the simulation-based algorithm described in Section 2.4 using

samples from the asymptotic distribution of the parameters. The Pohar-Perme estimator and its corresponding confidence intervals were calculated using the ‘relsurv’ R package. We observe that the results with the parametric and nonparametric approaches are very close, and that the confidence intervals for the parametric model are slightly shorter (as expected), which shows that the proposed parametric model can accurately capture the underlying hazard structure.

Model PHEW AHEW AFTEW GHEW
scale 0.059 (0.038) 8.482 (0.724) 1.190 (0.175) 1.838 (0.374)
shape 0.188 (0.014) 0.539 (0.046) 0.385 (0.012) 0.442 (0.033)
power 9.175 (1.420) 1.483 (0.129) 4.387 (0.312) 3.593 (0.368)
-0.112 (0.006) 0.032 (0.001) 0.041 (0.004)
-2.977 (0.282) 0.881 (0.065) 0.691 (0.311)
-6.680 (0.337) 1.909 (0.050) 1.707 (0.229)
-10.469 (0.416) 3.003 (0.046) 3.413 (0.226)
-2.668 (0.416) 0.744 (0.115) 0.822 (0.448)
-1.021 (0.106) 0.289 (0.029) 0.539 (0.114)
agediagc 0.022 (0.001) – () 0.034 (0.001)
Istage2 0.721 (0.056) – () 0.845 (0.069)
Istage3 1.473 (0.043) – () 1.849 (0.053)
Istage4 2.211 (0.041) – () 3.073 (0.050)
INCOME_SCORE_2015c 0.527 (0.085) – () 0.750 (0.150)
comorbidity 0.192 (0.021) – () 0.349 (0.039)
AIC 20523.141 20855.753 20189.124 20164.911
Table 3: Maximum likelihood estimates of the parameters (standard errors) for the different excess hazard models, with their corresponding AIC. The time dependent effects are indicated with the subindex ‘’. () By construction, the effects of covariates are constrained to be the same for the time-dependent and time fixed effects in the AFT model (5). See equations (2)-(6) for more details on the different hazard structures.
(a)
(b)
Figure 3: Shapes of the EW excess hazard for PH, AFT, AH, and GH models with baseline covariate values (solid lines) vs. Piecewise excess hazard (grey line): (a) comorbidity = 0 and Stage II (the 2 lowest curves/step functions in each cell), comorbidity = 0 and Stage IV (the 2 highest curves/step functions in each cell); and (b) comorbidity = 1 and Stage II (the 2 lowest curves/step functions in each cell), comorbidity = 1 and Stage IV (the 2 highest curves/step functions in each cell).
Total population by comorbidity
GHEW Pohar-Perme
Comorb. year NS lower upper NS lower upper
0 1 0.407 0.401 0.416 0.408 0.398 0.418
2 0.270 0.265 0.278 0.268 0.259 0.277
3 0.204 0.199 0.211 0.209 0.201 0.218
3.9 0.167 0.162 0.174 0.182 0.172 0.191
1 1 0.380 0.371 0.391 0.370 0.356 0.385
2 0.254 0.246 0.264 0.238 0.225 0.252
3 0.193 0.185 0.203 0.169 0.158 0.182
3.9 0.158 0.150 0.168 0.138 0.126 0.152
Age group 55-65 at Stage I by comorbidity
Comorb. year NS lower upper NS lower upper
0 1 0.924 0.914 0.935 0.928 0.897 0.960
2 0.842 0.827 0.860 0.837 0.794 0.883
3 0.770 0.752 0.791 0.797 0.749 0.847
3.9 0.712 0.692 0.737 0.762 0.705 0.823
1 1 0.881 0.869 0.896 0.901 0.851 0.953
2 0.772 0.754 0.793 0.744 0.674 0.821
3 0.684 0.662 0.710 0.670 0.595 0.755
3.9 0.618 0.595 0.647 0.627 0.541 0.725
Table 4: Net survival (NS) at years and confidence intervals, estimated using the GHEW model and the non parametric Pohar-Perme estimator.

5 Discussion

We have studied a general parametric hazard structure that can capture the basic shapes of the baseline hazard and that contains, as particular cases, the main models of interest in survival analysis: proportional hazards, accelerated failure time, and accelerated hazards models. PH and AFT models already enjoy popularity in survival analysis. However, the limited application of the AH model is mainly due to the lack of efficient and reliable estimation methods Chen et al. (2014). We have shown that, by assuming a flexible parametric distribution function such as the EW, it is possible to conduct classical likelihood inference using already available optimisation algorithms. The combination of such flexible parametric hazard function with the GH structure represents a powerful tool for modelling survival times. Although we have focused on the context of excess hazard models, the proposed flexible parametric GH model is also applicable in the context of overall survival. We have employed the EW distribution for modelling the baseline hazard in the GH model (6), however, there exist other flexible parametric distributions, such as the generalised Gamma (Cox and Matheson, 2014) and the generalised Weibull (Mudholkar et al., 1996) distributions, that can also capture the basic shapes of the hazard function. Here, we only employed the EW distribution as this choice allows for a parsimonious implementation, facilitates the interpretation of the parameters, and the corresponding maximum likelihood estimators (MLEs) are consistent and asymptotically normal in the presence of censored observations (Qian, 2012). The simulation study shows that the proposed model has good frequentist properties and that the selection of an appropriate model structure is feasible for large enough samples. The proposed model can also capture cases where the hazard and the survival curves cross, a case of great interest in practice. Cases when the AH and AFT structures produce crossing hazards have been studied in Zhang and Peng (2009), and these are also illustrated in our simulation study. Despite the flexibility of the EW distribution, there still exist baseline hazard shapes that cannot be captured by this model (such as multimodal hazards). We have conducted additional simulation studies in order to assess the effect of model misspecification on the estimation of net survival. We found that in cases where the EW distribution cannot capture the true shape of the baseline hazard, the net survival functions associated to the fitted models tend to be relatively close to the true model, and that the AIC selects the correct hazard structure with high proportion, provided that the departures from the shapes the EW can capture are moderate. If the departures are severe, this may, unsurprisingly, affect the selection of the correct hazard structure but the selected model tend to resemble the shape of the true generating model (see Hjort, 1992 for a general study of fitting parametric survival models under possible misspecification). These additional studies are available from the authors’ websites, as well as the R code for fitting the models detailed in this paper.

Although we have centered our attention on the GH structure (6), we point out the more general representation discussed in Chen et al. (2003):

where is a known non-negative function that defines the relationship between the baseline hazard and the covariates. This representation contains, for instance, the additive hazards model (Lin and Ying, 1994) , which we do not consider here as it requires additional conditions to guarantee that ; as well as structures including non-linear relationships. We also point out that the use of splines for including non-linear effects can be coupled with any of the aforementioned hazard structures. A potentially useful choice are B-splines, which are implemented in R.

We have employed the AIC for selecting the hazard structure that better fits the data since this criterion can account for possible model misspecification, as this tool asymptotically selects the model that minimises the Kullback-Leibler divergence between the fitted models and the true generating model, under mild regularity conditions. We point out that other criteria such as the Bayesian Information Criterion (BIC) or cross validation can be employed instead. An interesting feature of the GH model is that, by selecting the active variables, we automatically select the hazard structure that better fits the data. The study of efficient variable selection methods in the general structure (

6) would allow the identification of active variables as well as the underlying structure (PH, AH, HH, AFT, or GH) of the hazard function. This points out an interesting research line. The simulation study also illustrates the importance of accounting for sparsity, as the inclusion of spurious variables may bias the estimates of the active variables and the parameters of the baseline hazard. Since, in recent years, more variables such as comorbidities and types of treatments have become available at the population level, we believe this topic will become relevant in cancer epidemiology.

Acknowledgements

Funding: This research was partly supported by Cancer Research UK grant number C7923/A18525. The findings and conclusions in this report are those of the authors and do not necessarily represent the views of Cancer Research UK.

Ethical approvals: We obtained the ethical and statutory approvals required for this research (PIAG 1-05(c)/2007; ECC 1-05(a)/2010); ethical approval updated 6 April 2017 (REC 13/LO/0610). We attest that we have obtained appropriate permissions and paid any required fees for use of copyright protected materials.

References

  • Abrahamowicz et al. [1996] M. Abrahamowicz, T. Mackenzie, and J.M. Esdaile. Time-Dependent Hazard Ratio: Modeling and Hypothesis Testing with Application in Lupus Nephritis. Journal of the American Statistical Association, 91:1432–1439, 1996.
  • Benitez-Majano et al. [2016] S. Benitez-Majano, H. Fowler, C. Maringe, C. Di Girolamo, and B. Rachet. Deriving stage at diagnosis from multiple population-based sources: colorectal and lung cancer in England. British Journal of Cancer, 115(3):391, 2016.
  • Bolard et al. [2002] P. Bolard, C. Quantin, M. Abrahamowicz, J. Estève, R. Giorgi, H. Chadha-Boreham, C. Binquet, and J. Faivre. Assessing time-by-covariate interactions in relative survival models using restrictive cubic spline functions. Journal of Cancer Epidemiology and Prevention, 7:113–122, 2002.
  • Charvat and Belot [2017] H. Charvat and A. Belot. mexhaz: Mixed Effect Excess Hazard Models, 2017. URL https://CRAN.R-project.org/package=mexhaz. R package version 1.4.
  • Charvat et al. [2016] H. Charvat, L. Remontet, N. Bossard, L. Roche, O. Dejardin, B. Rachet, G. Launoy, and A. Belot. A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates. Statistics in Medicine, 35(18):3066–3084, 2016.
  • Chen et al. [2014] Y. Chen, T. Hanson, and J. Zhang. Accelerated hazards model based on parametric families generalized with Bernstein polynomials. Biometrics, 70(1):192–201, 2014.
  • Chen and Wang [2000a] Y. Q. Chen and M. Wang. Estimating a treatment effect with the accelerated hazards models. Controlled Clinical Trials, 21(4):369–380, Aug 2000a.
  • Chen and Jewell [2001] Y.Q. Chen and N.P. Jewell. On a general class of semiparametric hazards regression models. Biometrika, 88(3):687–702, 2001.
  • Chen and Wang [2000b] Y.Q. Chen and M.C. Wang. Analysis of accelerated hazards models. Journal of the American Statistical Association, 95(450):608–618, 2000b.
  • Chen et al. [2003] Y.Q. Chen, N.P. Jewell, and J. Yang. Accelerated hazards model: method, theory and applications. Handbook of Statistics, 23:431–441, 2003.
  • Cox and Matheson [2014] C. Cox and M. Matheson. A comparison of the generalized Gamma and exponentiated Weibull distributions. Statistics in Medicine, 33(21):3772–3780, 2014.
  • Danieli et al. [2012] C. Danieli, L. Remontet, N. Bossard, L. Roche, and A. Belot. Estimating net survival: the importance of allowing for informative censoring. Statistics in Medicine, 31(8):775–786, 2012.
  • Durrleman and Simon [1989] S. Durrleman and R. Simon. Flexible regression models with cubic splines. Statistics in Medicine, 8:551–561, 1989.
  • Esteve et al. [1990] J. Esteve, E. Benhamou, M. Croasdale, and L. Raymond. Relative survival and the estimation of net survival: elements for further discussion. Statistics in Medicine, 9(5):529–538, 1990.
  • Gamel and Vogel [2001] J.W. Gamel and R.L. Vogel. Non-parametric comparison of relative versus cause-specific survival in surveillance, epidemiology and end results (SEER) programme breast cancer patients. Statistical Methods in Medical Research, 10(5):339–352, 2001.
  • Giorgi et al. [2003] R. Giorgi, M. Abrahamowicz, C. Quantin, P. Bolard, J. Estève, J. Gouvernet, and J. Faivre. A relative survival regression model using B-spline functions to model non-proportional hazards. Statistics in Medicine, 22:2767–2784, 2003.
  • Hakulinen and Tenkanen [1987] T. Hakulinen and L. Tenkanen. Regression analysis of relative survival rates. Applied Statistics, pages 309–317, 1987.
  • Hess and Levin [2014] K. R Hess and V.A. Levin. Getting more out of survival data by using the hazard function. Clinical Cancer Research, 20(6):1404–1409, 2014.
  • Hess [1994] K.R. Hess. Assessing time-by-covariate interactions in proportional hazards regression models using cubic spline functions. Statistics in Medicine, 13:1045–1062, 1994.
  • Hjort [1992] N.L. Hjort. On inference in parametric survival data models. International Statistical Review, pages 355–387, 1992.
  • Kalbfleisch and Prentice [2011] J.D. Kalbfleisch and R.L. Prentice. The Statistical Analysis of Failure Time Data, volume 360. John Wiley & Sons, 2011.
  • Khan [2017] S.A. Khan. Exponentiated Weibull regression for time-to-event data. Lifetime Data Analysis, pages 1–27, 2017.
  • Lin and Ying [1994] D.Y. Lin and Z. Ying. Semiparametric analysis of the additive risk model. Biometrika, 81(1):61–71, 1994.
  • Mahboubi et al. [2011] M. Mahboubi, A.and Abrahamowicz, R. Giorgi, C. Binquet, C. Bonithon-Kopp, and C. Quantin. Flexible modeling of the effects of continuous prognostic factors in relative survival. Statistics in Medicine, 30(12):1351–1365, 2011.
  • Mandel [2013] M. Mandel. Simulation-based confidence intervals for functions with complicated derivatives. The American Statistician, 67(2):76–81, 2013.
  • Maringe et al. [2017] C. Maringe, H. Fowler, B. Rachet, and M.A. Luque-Fernandez. Reproducibility, reliability and validity of population-based administrative health data for the assessment of cancer non-related comorbidities. PloS one, 12(3):e0172814, 2017.
  • Mariotto et al. [2014] A.B. Mariotto, A.M. Noone, N. Howlader, H. Cho, J. Keel, G.E .and Garshell, S. Woloshin, and L.M. Schwartz. Cancer survival: an overview of measures, uses, and interpretation. Journal of the National Cancer Institute Monographs, 2014(49):145–186, 2014.
  • Mudholkar et al. [1995] G.S. Mudholkar, D.K. Srivastava, and M. Freimer. The exponentiated Weibull family: a reanalysis of the bus-motor-failure data. Technometrics, 37(4):436–445, 1995.
  • Mudholkar et al. [1996] G.S. Mudholkar, D.K. Srivastava, and G.D. Kollia. A generalization of the Weibull distribution with application to the analysis of survival data. Journal of the American Statistical Association, 91(436):1575–1583, 1996.
  • Nelson et al. [2007] C.P. Nelson, P.C. Lambert, I.B. Squire, and D.R. Jones. Flexible parametric models for relative survival with application in coronary heart disease. Statistics in Medicine, 26(30):5486–5498, 2007.
  • Perme et al. [2012] M.P. Perme, J. Stare, and J. Estève. On estimation in relative survival. Biometrics, 68(1):113–120, 2012.
  • Perme et al. [2016] M.P. Perme, J. Estève, and B. Rachet. Analysing population-based cancer survival–settling the controversies. BMC cancer, 16(1):933, 2016.
  • Qian [2012] L. Qian. The Fisher information matrix for a three-parameter exponentiated Weibull distribution under type II censoring. Statistical Methodology, 9(3):320–329, 2012.
  • Remontet et al. [2007] L. Remontet, N. Bossard, A. Belot, and J. Esteve. An overall strategy based on regression models to estimate relative survival and model the effects of prognostic factors in cancer survival studies. Statistics in Medicine, 26(10):2214–2228, 2007.
  • Ross [2006] S.M. Ross. Simulation. Statistical Modeling and Decision Science. Elsevier Science, 2006. ISBN 9780080517223. URL http://books.google.fr/books?id=TMUt5OXVvY0C.
  • Royston and Sauerbrei [2008] P. Royston and W. Sauerbrei. Multivariable model-building: a pragmatic approach to regression anaylsis based on fractional polynomials for modelling continuous variables, volume 777. John Wiley & Sons, 2008.
  • Rubio and Genton [2016] F.J. Rubio and M.G. Genton.

    Bayesian linear regression with skew-symmetric error distributions with applications to survival analysis.

    Statistics in Medicine, 35(14):2441–2454, 2016.
  • Rubio and Hong [2016] F.J. Rubio and Y. Hong. Survival and lifetime data analysis with a flexible class of distributions. Journal of Applied Statistics, 43(10):1794–1813, 2016.
  • Sleeper and Harrington [1990] L.A. Sleeper and D.P. Harrington. Regression Splines in the Cox Model with Application to Covariate Effects in Liver Disease. Journal of the American Statistical Association, 85:941–949, 1990.
  • Zhang and Peng [2009] J. Zhang and Y. Peng. Crossing hazard functions in common survival models. Statistics & Probability Letters, 79(20):2124–2130, 2009.