A Flexible Parametric Modelling Framework for Survival Analysis

We introduce a general, flexible, parametric survival modelling framework which encompasses key shapes of hazard function (constant, increasing, decreasing, up-then-down, down-then-up), various common survival distributions (log-logistic, Burr type XII, Weibull, Gompertz), and includes defective distributions (i.e., cure models). This generality is achieved using four basic distributional parameters: two scale-type parameters and two shape parameters. Generalising to covariate dependence, the scale-type regression components correspond to accelerated failure time (AFT) and proportional hazards (PH) models. Therefore, this general formulation unifies the most popular survival models which allows us to consider the practical value of possible modelling choices for survival data. Furthermore, in line with our proposed flexible baseline distribution, we advocate the use of multi-parameter regression in which more than one distributional parameter depends on covariates - rather than the usual convention of having a single covariate-dependent (scale) parameter. While many choices are available, we suggest introducing covariates through just one or other of the two scale parameters, which covers AFT and PH models, in combination with a `power' shape parameter, which allows for more complex non-AFT/non-PH effects, while the other shape parameter remains covariate-independent, and handles automatic selection of the baseline distribution. We explore inferential issues in simulations, both with and without a covariate, with particular focus on evidence concerning the need, or otherwise, to include both AFT and PH parameters. We illustrate the efficacy of our modelling framework by investigating differences between treatment groups using data from a lung cancer study and a melanoma study. Censoring is accommodated throughout.

Authors

• 14 publications
• 2 publications
• 2 publications
11/16/2021

Multi-Parameter Regression Survival Modelling with Random Effects

We consider a parametric modelling approach for survival data where cova...
01/10/2019

Multi-Parameter Regression Survival Modelling: An Alternative to Proportional Hazards

It is standard practice for covariates to enter a parametric model throu...
05/21/2018

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

The proportional hazards model represents the most commonly assumed haza...
01/10/2019

A Bivariate Power Generalized Weibull Distribution: a Flexible Parametric Model for Survival Analysis

We are concerned with the flexible parametric analysis of bivariate surv...
06/26/2021

Parmsurv: a SAS Macro for Flexible Parametric Survival Analysis with Long-Term Predictions

Health economic evaluations often require predictions of survival rates ...
10/18/2018

Impact of model misspecification in shared frailty survival models

Survival models incorporating random effects to account for unmeasured h...
04/29/2021

Bivariate Analysis of Birth Weight and Gestational Age Depending on Environmental Exposures: Bayesian Distributional Regression with Copulas

In this article, we analyze perinatal data with birth weight (BW) as pri...
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

This article is concerned with both theoretical and practical aspects of parametric survival analysis with a view to providing an attractive and flexible general modelling approach to analysing survival data in areas such as medicine, population health, and disease modelling. In particular, focus will be on the choice of an appropriate flexible form for the distribution of the survival outcome and the efficient use of multi-parameter regression to understand the effects of covariates on survival.

We consider a univariate lifetime random variable,

, the primary survival outcome, whose cumulative hazard function (c.h.f.), , is, atypically perhaps, modelled using a flexible parametric form which we take to be

 H(t)=λH0((ϕt)γ;κ),     t>0. (1)

Here, is an underlying c.h.f. with shape parameter , and , and are further parameters with the following distinct interpretations: controls the horizontal scaling of the hazard function, and is well known as the accelerated failure time parameter (also, is the usual distributional scale parameter); controls the vertical scaling of the hazard function, and is well known as the proportional hazards parameter; and is a second shape parameter which is explicitly defined as a power parameter (unlike

which can enter in potentially more complicated ways, and might even represent a vector of parameters). Were

to be modelled as a location-scale distribution on , then and would be the location and scale of that distribution, respectively, these relationships driving our preference to specify as a power parameter rather than as a more general shape parameter. As will be clear in the sequel, we intend that only one of the scale parameters be present in the model in order to avoid identifiability issues (i.e., we fix or ). However, we write the model in a general way for the purpose of unification of sub-models.

In this article, we also propose a specific choice for , namely

 HA(t;γ,κ)=κ+1κ{(1+tγκ+1)κ−1},     t>0. (2)

corresponding to an adapted form of the ‘power generalised Weibull’ (PGW) distribution introduced by Bagdonaviçius and Nikulin (2002); we will use APGW to stand for ‘adapted PGW’. This choice has some major advantages: with just two shape parameters, the full range of simplest hazard shapes, namely, constant, increasing, decreasing, up-then-down or down-then-up (and no others), are available, the parameters and controlling this through the way they control behaviour of the hazard function near zero and at infinity. Here, we use the simple descriptive terms ‘up-then-down’ and ‘down-then-up’ to avoid the term ‘bathtub-shaped’, which is down-then-up but with a flat valley, the clumsy term ‘upside-down-bathtub-shaped’, and the terms ‘unimodal/uniantimodal’ which also encompass monotone hazards. Our adaptation of the PGW distribution also allows to control distributional choice within the family: for , log-logistic and Burr Type XII distributions are the heaviest tailed members, Weibull distributions () are ‘central’ within the family, and Gompertz-related distributions are the most lightly tailed. See Section 2 for details of this model, which also include its cure model special cases when .

Any one or more of the four distributional parameters in model (1) can be made to depend, typically log-linearly, on covariates; such “multi-parameter regression” is one of the focusses of this work. Indeed, this general formulation covers the most popular survival models, e.g., the accelerated failure time (AFT) model when depends on covariates, the proportional hazards (PH) model when depends on covariates, and semi-parametric versions when is an unspecified function. In particular, an advantage of considering (1) is that one may evaluate the breadth of possible modelling choices. Our primary focus in this respect is to consider which distributional parameters should depend on covariates to assess, for example, whether an AFT model ( regression) is, in general, likely to provide a superior fit when compared with a PH model ( regression), the utility of a simultaneous AFT-PH model (simultaneous and regression components; Chen and Jewell (2001)) when , and the merits of a shape regression component ( or ) in addition to the, more standard, AFT and PH components. One might also consider whether or not non-parametric components should be introduced either for functions of covariates within the regression equations, or for the baseline c.h.f., , or both. However, this is beyond the scope of the current paper.

The reason for our focus on the core model structure rather than the development of non-/semi-parametric approaches is that, within the survival literature, there is a general over-emphasis placed on semi-parametric models – compared with other fields of statistics – to the extent that many useful parametric alternatives do not receive the attention they deserve. In particular, practitioners are often content with the “flexibility” afforded by a non-parametric baseline function without concerning themselves with the possibly inflexible structural assumptions of the model at hand. Indeed, a structurally flexible parametric framework has the potential to outperform a less flexible semi-parametric model; for example, there might be more to be gained by contemplating the extension of a PH model ( regression) to include a shape regression, than by extension to a non-parametric . Of course, this is not to downplay the importance of a sufficiently flexible baseline function, and our proposed choice for , (2), is quite general as it covers a wide variety of popular survival distributions.

After Section 2, in which we justify our choice of baseline distribution and develop its properties, we consider the extension to regression modelling in Section 3

, including model interpretation and estimation. Then, the properties of estimation within this general framework, and further practical aspects, are explored using simulated and real data in Sections

4 and 5, respectively. Finally, we close with some discussion in Section 6.

2 The Specific Model for H0

2.1 Basic Definition and Properties

We recommend for general use the APGW distribution with c.h.f. given by (2) and hazard function is

 hA(t;γ,κ)=γtγ−1(1+tγκ+1)κ−1,     t>0. (3)

This is a tractable distribution with readily available formulae for its (unimodal) density, survivor and quantile functions also. It comes about by a particular vertical and horizontal rescaling of the original PGW distribution which has c.h.f.

(see Bagdonaviçius and Nikulin (2002), Nikulin and Haghighi (2009) and Dimitrakopoulou et al. (2007); the special case of

is the extended exponential distribution of

Nadarajah and Haghighi (2011)). This resulting APGW distribution then retains attractive shape properties of the PGW distribution’s hazard function, includes important survival distributions as special and limiting cases and extends to cure models, as we now show.

First, for fixed ,

 hA(t;γ,κ)∼γtγ−1ast→0   and   hA(t;γ,κ)∼(κ+1)1−κγtκγ−1ast→∞.

The power parameter controls the behaviour of the hazard function at zero: it goes to as As , the hazard function goes to as In fact, the APGW hazard function joins these tails smoothly in such a way that its hazard shapes are readily shown to be as listed in Table 1. Whenever the hazard function is non-monotone, its mode/antimode is at

Defining by (2) allows us to identify an especially large number of special and limiting cases, many important and well known, some less so, as listed in Table 2. (For the ‘Weibull extension’ distribution, see Chen (2000) and Xie et al. (2002).) The shapes of their hazard functions, which are also given in Table 2, reflect the general shape properties of Table 1, of course.

It can also be shown that the APGW distribution retains membership of the log-location-scale-log-concave family of distributions of Jones and Noufaily (2015) and therefore, inter alia, unimodality of densities. We also now note, for future reference, the attractive form of the quantile function associated with , namely

The new adaptation can also be used to widen the family of PGW distributions by taking . For clarity, temporarily define so that . The APGW c.h.f. can then be written as

 HA(t;γ,ψ)=ψ1−ψ(1−1{1+(tγ/ψ)}1−ψ).

This corresponds to a cure model with cure probability

Since the (improper) survival function is in this case of the form , this cure model has an interpretation as the distribution of the minimum of a Poisson number of random variables (e.g. cancer cells, tumours), each following the lifetime distribution with survival function (e.g. Tsodikov et al. (2003)); here, the Poisson parameter is and is the survival function of a scaled Burr Type XII distribution. The hazard functions follow the shape of their limit — the log-logistic — being decreasing for and up-then-down otherwise.

The PGW distribution, and in slightly more complicated form the APGW distribution, exhibit interesting frailty relationships between members of the families with different values of . We defer consideration of these frailty links to Jones et al. (2018)

where we exploit them to obtain a useful bivariate shared frailty model with PGW/APGW marginal distributions. In addition, PGW and APGW distributions are written as linear transformation models in the Appendix.

2.2 Why This Particular Choice for H0?

The PGW/APGW distribution shares the set of hazard behaviours listed in Table 1 with two other established two-shape-parameter lifetime distributions centred on the Weibull distribution, namely, the generalised gamma (GG) and exponentiated Weibull (EW) distributions; see Jones and Noufaily (2015). See Figure 1 for many examples of just how similar the hazard shapes of all three distributions are; in Figure 1, we have chosen the scale parameter such that each distribution has median one, used the PGW vertical scaling and otherwise specified shape parameters only so that all three hazard functions behave as as and as as .

Further effort to choose shape parameters to match hazard functions or other aspects of the distributions even more closely is possible and has been pursued for the EW and GG distributions by Cox and Matheson (2014) and extended to the PGW distribution (what they call the generalised Weibull distribution) by Matheson et al. (2017). Cox and Matheson (2014) state that the “agreement between the two distributions [GG and EW] in our various comparisons, both graphically and in terms of the K–L [Kullback-Leibler] distance, is striking”; after a similar K–L matching exercise, Matheson et al. (2017)

state that “the survival and hazard functions of the [PGW] distribution and its matched GG are visually indistinguishable.” It remains, therefore, to choose between APGW, GG and EW distributions on other grounds. The GG distribution includes the Weibull and gamma distributions as special cases and the lognormal as a limiting one; the EW distribution includes the Weibull and exponentiated exponential distributions. However: we have been unable to match the number of APGW’s advantageous properties as in the previous subsections by similarly adapting either the GG or EW distributions; we prefer the breadth of/difference between the wide range of distributions encompassed by the APGW distribution; and we appreciate the greater tractability of the APGW distribution both mathematically and computationally (for instance, its hazard function has a simpler form compared with the GG — which involves an incomplete gamma function — and the EW).

3 Regression

3.1 Modelling Choices

Within our proposed APGW modelling framework, there are four parameters, , and , which could potentially depend on covariates. Note that most classical modelling approaches are based on having a single covariate-dependent distributional parameter, which we refer to as single parameter regression (SPR), where, understandably, there is a particular emphasis on scale-type parameters, e.g., the accelerated failure time (AFT) model ( regression) and the proportional hazards (PH) model ( regression). However, in line with the flexibility of the APGW distribution itself, we also consider taking a flexible multi-parameter regression (MPR) approach in which more than one parameter may depend on covariates (cf. Burke and MacKenzie (2017), and references therein, for details of multi-parameter regression). The most general linear APGW-MPR is, therefore, given by

 log(ϕ) =xTτ, log(λ) =xTβ, log(γ) =xTα, log(κ+1) =xTν,

where log-link functions are used to respect the positivity of the parameters , and , with a slightly different link function for to accommodate the fact that, within our APGW, it can take values in the range (see Section 2.1), is a vector of covariates, and , , , and are the corresponding vectors of regression coefficients. In practice, we may not necessarily have the same set of covariates appearing in all regression components, and, in our current notation, this can be handled by setting various regression coefficients to zero.

As mentioned in Section 1, we could extend the above regression specification via non-parametric regression functions of , but this is beyond the scope of this paper and, indeed, the MPR approach is, in itself, already flexible without this added complexity. Furthermore, although the general APGW-MPR model offers the opportunity of four regression components simultaneously, this full flexibility is unlikely to be required in practice. In particular, it is well known that and coincide in the Weibull distribution so that only one scale parameter is needed in this case, i.e., and are non-identifiable when . Moreover, our numerical studies (Sections 4 and 5) suggest that, outside of the Weibull, this is effectively much more generally true. Specifically, although and are theoretically identifiable in the non-Weibull cases, the parameters are nonetheless nearly non-identifiable in finite samples, which is an apparently new finding in the literature. Thus, in general, we should fix either or , but we would not simultaneously fix as a scale parameter is a core component for most statistical models.

A good practical choice is composed of the following pieces: (a) only one scale parameter ( or ) depends on covariates, while the other is fixed at one as mentioned above, (b) the shape parameter may depend on covariates, and (c) the shape parameter is constant, i.e., only the intercept, , is non-zero in the vector. This choice provides a useful framework which incorporates, depending on the choice of scale regression, either an AFT () or PH () component, allows for non-AFT/non-PH effects via the coefficients associated with the power parameter (Section 3.2), and automatically selects the underlying baseline distribution via from a range of popular survival distributions (Table 2) including defective distributions, i.e., cure models (Section 2.1).

3.2 Suggested models: M(τ,α) and M(β,α)

Let and denote the two models suggested in the previous paragraph, e.g., the latter is the model with and regression components, along with the shape parameter (but where is a vector of zeros). More generally, beyond these two models, we will use this notation throughout the paper where the arguments of indicate which regression components are present in the model, the absence of either or indicating that this is a vector of zeros. Irrespective of the presence of or regression components in this notation, we assume that and are always present since these are needed to characterise the baseline distribution and the shape of its hazard function (see Tables 1 and 2). Thus, for example, and are, respectively, AFT and PH models with two shape parameters ( and ), is a model which extends the suggested model so that the regression component is also present, and is the most general APGW-MPR model.

We first consider model which extends the basic AFT model, , via the incorporation of the regression component. Now suppose that is a binary covariate and let be the covariate vector with set to zero so that we may write and . As this model extends the AFT model, it is natural to consider its quantile function which is given by

 Q(u|x)=exp(−xTτ)QA1(u;κ)exp(−xTα)

where is the “baseline” quantile function defined in Section 2.1. We can then inspect the quantile ratio

 QRj(u)=Q(u|xj=1)Q(u|xj=0)=exp(−τj)QA1(u;κ)exp(−xT(−j)α){exp(−αj)−1}

where we see that is the key parameter in determining the -dependence. In particular, since is an increasing function of , increases when , decreases when , and is constant (i.e., the usual AFT case) when . Hence, the coefficient characterises the nature of the effect of the binary covariate , and provides a test of the AFT property for that covariate.

Now consider the model which extends the PH model, , and whose hazard function is given by

 h(t|x)=exp(xTβ)hA(t;exp(xTα),κ)

where is defined in (3). The hazard ratio for the binary covariate is

 HRj(t)=h(t|xj=1)h(t|xj=0)=exp(βj+αj)texp(xT(−j)α){exp(αj)−1}g(t;αj,xT(−j)α,κ)

where

 g(t;αj,xT(−j)α,κ)=⎛⎝texp(αj+xT(−j)α)+κ+1texp(xT(−j)α)+κ+1⎞⎠κ−1.

Clearly, characterises departures from proportional hazards as is a constant when . For , we have that and when , while and when . Furthermore, it can be shown that varies monotonically in in the following cases: (i) , or (ii) and . (We do not know about monotonicity or otherwise in the remaining cases.)

From the above we see that, within the suggested and models, the parameters play the following roles: the scale coefficients ( or ) control the overall size of the effect where negative coefficients correspond to longer lifetimes; the shape coefficients describe how covariate effects vary over the lifetime, i.e., permitting non-AFT and non-PH effects; and the shape parameter characterises the baseline distribution. Note that we could, alternatively, achieve non-constant covariate effects via the regression component rather than the component, i.e., using rather than . However, in this case, the interpretation is that such non-constant effects are due to populations which arise from structurally different distributions, rather than different shapes within a given baseline distribution. The latter is arguably more natural as it creates a separation of parameters whereas, in the former, distribution selection and non-constant covariate effects are intertwined. Of course, this is not to say that models with components instead of, or in combination with, components will never be useful in practice. However, we are highlighting practical merits of the and models and, indeed, the general use of these models is motivated by the numerical studies of Sections 4 and 5.

3.3 Estimation

Consider the model formulation given in (1) with all four regression components, i.e., the model. (While we advocate the special cases or , we write the estimation equations in a general form below so as to unify all potential model structures. In particular, estimation of both and is not recommended in practice.) Let , , and be the covariate-dependent distributional parameters for the th individual with covariate vector , and , , , and are the associated vectors of regression coefficients. Allow independent censoring by attaching to each individual an indicator which equals one if the response is observed, and zero if it is right-censored. The log-likelihood function is then given by

 ℓ(θ) =n∑i=1[δi{log(λiγiziti)+m0(zi;κi)}−λiH0(zi;κi)]

where , and, in our proposed APGW case,

 H0(t;κ)=HA(t;1,κ)=κ+1κ{(1+tκ+1)κ−1},
 m0(t;κ)=logh0(t;κ)=(κ−1)log(1+tκ+1).

As usual, the log-likelihood function can be maximised by solving the score equations

 (UTτX, UTβX, UTαX, UTνX)T=04p×1

where is an matrix whose th row is , is a vector of zeros and , , , and are vectors whose th elements are as follows:

 Uτ,i =δi{γi+γizim′0(zi;κi)}−λiγizih0(zi;κi), Uβ,i =δi−λiH0(zi;κi), Uα,i =δi[1+log(zi){1+zim′0(zi;κi)}]−λizilog(zi)h0(zi;κi), Uν,i =[δi{κi+1κi−1m0(zi;κi)−zim′0(zi;κi)} −λiκi+1κ2i(κi−1){κi−1+(1+tκi+1)κia0(zi;κi)}],

where .

Note that the vectors , and are written generically so that they apply to any model of the form given in (1), i.e., they are not specific to the APGW case; the form of , on the other hand, uses the way in which and hence and depend on Thus, although the APGW is certainly a flexible choice (see Section 2), the first three score components extend immediately to other baseline distributions by replacing (and, consequently, and ). Estimation then proceeds once the functional form of has been re-evaluated.

Furthermore, one may, alternatively, prefer to maintain an unspecified baseline distribution, whereby represents an infinite-dimensional (possibly covariate-independent) parameter vector. In this case, estimation equations for the regression coefficients , , and can be based on where is replaced with an appropriate non-parametric estimator (and, similarly, for and ). However, while non-parametric estimation of is reasonably straightforward (say, using a Nelson-Aalen-type estimator), it is well known that terms such as and , which involve and , are more difficult to estimate consistently. We note that semi-parametric versions of the and models have respectively been developed by Chen and Jewell (2001) and Burke et al. (2018). However, we are unaware of a semi-parametric model in the literature. In any case, such semi-parametric models are beyond the scope of the current paper and, indeed, a flexible parametric framework can cover a wide variety of applications as previously discussed in Section 1.

4 Simulation Studies

4.1 Without Covariates

Before considering estimation in the presence of covariates, we first investigate estimation in the context of the APGW model with no covariates. Thus, we simulated data from the APGW distribution parameterised in terms of the following unconstrained parameters: , , and . The values of the first three parameters were fixed at , , , respectively, while was varied such that (rounded to two decimal places); note that , and correspond, respectively, to the log-logistic (), Weibull () and Gompertz () distributions. Furthermore, the sample size was fixed at 1000 and censoring times were generated from an exponential distribution such that, for each value, the censoring rate was fixed at approximately 30%. Within each of the seven simulation scenarios (i.e., varying ), we fitted three different models with the aim of understanding the identifiability of parameters in a finite, but reasonably large, sample: (i) estimate all parameters, (ii) fix at its true value, and (iii) fix at zero. Thus, , and are estimated in all three models. We also considered additional scenarios where but the results are similar and, therefore, are not shown here.

Each scenario was replicated 1000 times, and the results are contained in Table 3

. Clearly, estimation is unstable in model (i), i.e., standard errors are large. This instability arises as a consequence of attempting to estimate the scale parameters,

and , simultaneously. Indeed, in all cases where these parameters are estimated simultaneously, we have found that corr. Of course, it is well known that and are perfectly co-linear in the Weibull case (), but it is interesting to find that this extends (approximately) beyond the Weibull distribution. This appears to be a new finding in survival modelling and implies that these parameters play somewhat similar roles across a range of popular lifetime distributions (it also explains the large standard errors observed in Table 2 of Chen and Jewell (2001)). This instability vanishes once is fixed. In particular, when is set to its true value of (i.e., model (ii)), the estimates display very little bias. Moreover, when is set to an incorrect value, (i.e., model (iii)), converges consistently to a value in the range 1.4–1.5 which compensates for the incorrect specification of and varies smoothly with ; the value of changes somewhat from its value in model (ii), but, interestingly, does not. Furthermore, the fitted survivor curves for both models (not shown) are close to the truth, i.e., there is no reduction in quality of model fit as a consequence of fixing to zero. Similarly (but not shown here), estimation is also stable if is fixed and is estimated, and the fitted survivor curves are again close to the truth. Therefore, the choice of scale — either or (fixing the other to zero) — behaves, approximately, as a model reparameterisation (which it is, exactly, in the Weibull case). We note that, for both models (ii) and (iii), the standard error of can be large when is large. However, this is not a concern as it is a consequence of the fact that the APGW distribution changes very little over a range of large values.

4.2 With a Covariate

We simulated survival times according to the APGW distribution with parameters , , , and where , was varied according to the set , and the remaining parameter values were fixed at , , , , and ; these values were selected to yield realistic survival times. In the notation of Section 3.2, the true model is . As in Section 4.1, the sample size and censored proportion were, respectively, set at 1000 and 30% (with censoring times generated from an exponential distribution). Within each of the seven scenarios (i.e., varying ), we fitted the following three regression models: the more general , the true , and the misspecified , respectively. The results, based on 1000 simulation replicates, are given in Table 4.

Mirroring the case with no covariates (Section 4.1), we find that estimation is unstable when attempting to estimate and coefficients simultaneously in , whereas estimation is stable in both the true and the misspecified models. In the latter, coefficients converge consistently to values varying smoothly with . Note that the results are broadly similar for smaller sample sizes of and (see Appendix for details).

We now consider model fit by inspecting the estimated baseline survivor curves, i.e., the survivor curve for an individual with which we denote by . In particular, we focus on this estimated baseline survivor function evaluated at three true quantiles, namely, , , since is an estimate of the probability . Boxplots of these estimates over simulation replicates arising from the true model, , and the misspecified model, , are shown in Figure 2. We also display the estimates from two simpler (misspecified) models, and , wherein has been dropped from the component (specifically, is set to zero in these simpler models). Clearly both and fit the data very well (apart from a little bias in when ), i.e., the choice of using a or regression component does not alter the model fit much (again mirroring the findings of Section 4.1). On the other hand, when the regression is dropped, the quality of the model fit decreases considerably as this represents a model misspecification in a much stronger sense than switching from to .

5 Data Analysis

5.1 Lung Cancer

We now consider our modelling framework in the context of a lung cancer study which was the subject of a 1995 Queen’s University Belfast PhD thesis by P. Wilkinson (see Burke and MacKenzie (2017)). This study concerns 855 individuals who were diagnosed with lung cancer between 1st October 1991 and 30th September 30 1992, who were then followed up until 30th May 1993 (approximately 20% of survival times were right-censored). The main aim of this study was to investigate differences between the following treatment groups: palliative care, surgery, chemotherapy, radiotherapy, and a combined treatment of chemotherapy and radiotherapy. In our analysis we take palliative care (which is a non-curative treatment providing pain relief) as the reference category. Note that, while various other covariates were captured for each individual, our main aim here is to explore the flexibility of our general modelling scheme in the context of the treatment variable.

As discussed in Section 3, we advocate the use of and since they offer a flexible extension of the popular PH and AFT models (i.e., and , respectively) in which the coefficients indicate non-PH/non-AFT effects (see Section 3.2), and where the baseline distribution is selected via the parameter . Thus, we fitted these two models, and their simpler PH and AFT counterparts, to the lung cancer data. We also fitted and for comparison, albeit we have argued in Section 3.2 that these are perhaps somewhat less natural. These six fitted models are summarised in Table 5.

We immediately see that the largest AIC/BIC values are associated with the simpler single component (i.e., and only) models which suggests that these models are not sufficiently flexible to capture the more complex non-PH/non-AFT effects observed here. Although, in this particular application, the AFT ( only) model has lower AIC/BIC values than that of the PH ( only) model, the fit can be greatly improved by modelling shape (either or ) in addition to scale. Although the best-fitting model here is , the difference is negligible compared with the models we favour, and . Interestingly, these latter two models have very close AIC/BIC values, indicating that the choice of or component is not at all important here (in line with the findings of Section 4.2). The use of more than two regression components did not yield further improvements in fit (models not shown), and, moreover, estimation of such models tends to be unstable — particularly, of course, those with two scale regression components (see also Section 4). Note that we have also avoided shape-only regression models, i.e., , , and , as, typically, models without scale components are not of interest, and, as we would expect, these models fit the current data very poorly indeed (with ).

We now consider the PH-APGW model, , and the two associated shape-regression extensions, and , in more detail. The advantage, in terms of model fit, of shape regression components is clear from Figure 3, while the and models themselves are virtually indistinguishable. Table 6 displays the estimated regression coefficients. We can see that both and suggest a baseline distribution which is between a log-logistic () and a Weibull (), while assumes a separate baseline distribution for each treatment group. Interestingly, in all three models, all shape parameters ( and ) are positive which indicates that the hazards are increasing with time in each treatment group (Table 1). While all three models are in agreement when it comes to the overall effectiveness of each treatment as viewed in terms of the scale coefficients (albeit the chemotherapy effect is only statistically significant in ), the positive shape coefficients in suggest that the effectiveness of each treatment reduces to some extent over time (see Section 3.2) – especially in the case of the combined treatment of chemotherapy and radiotherapy.

The hazard ratios for the models in Table 6 are shown in Figure 4 where those of and are quite similar. They suggest that while the various treatments reduce the hazard in the first few months, their effect is weakened over time and, perhaps, even become inferior to palliative care in the longer term (however, note that very few individuals remain in the sample beyond 15 months). Clearly SPR models, such as , cannot account for covariate effects of this sort.

It is worth highlighting the fact that the basic findings here are qualitatively similar to those of Burke and MacKenzie (2017) who analysed this lung cancer dataset using , i.e., a Weibull MPR model. However, the framework of the current paper permits us to consider a much wider range of model structures and distributions in which appears as a special case. In particular, from Table 6

yields a 95% confidence interval for

, , which does not support the Weibull () baseline distribution. Furthermore, , and we can confirm that the improvement in quality of fit is most evident in the palliative care group (which does not capture so well). Thus, although the basic findings are unaltered in this particular application, the APGW MPR approach yields a better solution in which uncertainty in selecting the baseline distribution is accounted for. Of course, the APGW MPR model will readily adapt to other applications which might differ significantly (both qualitatively and quantitatively) from .

5.2 Melanoma

The Eastern Cooperative Oncology Group (ECOG) trial “EST 1684” was a randomized controlled trial to investigate the adjuvant (i.e., post-surgery) chemotherapy drug “IFN-2b” in treating melanoma (Kirkwood et al., 1996). The outcome variable was relapse-free survival, i.e., time from randomization until the earlier of cancer relapse or death. Patients were recruited to the study between 1984 and 1990, and the study ended in 1993. In total, 284 patients were recruited of which 140 were assigned to the control group, and 144 were assigned to the treatment group. This dataset is available in the R package smcure (Chao et al., 2012), and variations of it have appeared in the cure model literature (Chen et al., 1999; Ibrahim et al., 2001).

As in Section 5.1, we fitted the following models: , , , , , and ; the results are summarised in Table 7. In this case, the two-component models do not provide a large improvement over the one-component models, and has the lowest AIC and the second-lowest BIC ( has the lowest BIC); the AIC and BIC values for are not much larger than for . The models , , and are compared to the Kaplan-Meier (KM) curves in Figure 5. The fitted and curves are similar, and are close to the KM curves. The fitted curves converge later in time, which, visually, look worse compared to the KM curves. However, note that there is very little data in the right tail so that converging curves are plausible when viewed with the level of uncertainty in the tail.

The parameter estimates for , , and are shown in Table 8. Firstly note that the scale coefficient of treatment is negative which indicates that treatment improves survival. Furthermore, for and , note that the parameter is negative, and, similarly, for , the estimated is negative in each treatment group. Thus, all three models point towards a cure proportion. The fact that the power shape parameter, , is greater than one for all three models means that the hazard function has an up-then-down shape. This type of hazard function is commonly observed in the cure literature since the population becomes increasingly composed of cured individuals (i.e., zero hazard) over time.

When , the APGW-MPR cure proportion is given by, i.e., the cure proportion depends on and . Therefore, and suggest that the cure proportion depends on treatment, while suggests that it does not. The estimated cure proportions, along with 95% confidence intervals, are shown in Table 9. Had we fixed to a log-logistic baseline (i.e., ), the resulting and models (not shown) provide an extremely poor fit to the data. This is noteworthy as even the heaviest-tail non-cure APGW model is not supported by the data (and, of course, a Weibull baseline is worse still). The heaviness of tail here can only be supported within the APGW family by a cure model.

6 Discussion

Our proposed APGW-MPR modelling framework is highly flexible and can adapt readily to a wide variety of applications in survival analysis and reliability. In particular, this framework includes the practically important AFT and PH models, and generalises them through shape regression components. Furthermore, the APGW baseline model covers the primary shapes of hazard function (constant, increasing, decreasing, up-then-down, down-then-up) within some of the most popular survival distributions (log-logistic, Burr type XII, Weibull, Gompertz) using only two shape parameters.

In practice, the full four-component APGW-MPR model is likely to be more flexible than is required for most purposes. In fact, we suggest that covariates should appear via just one scale-type component ( or ), along with the shape component which permits survivor functions with differing shapes and indicates departures from more basic AFT or PH effects, while is a covariate-independent parameter which allows us to choose among distributions within one unified framework. We have found that the scale-type parameters ( and ) are highly intertwined in the sense that they cannot be estimated simultaneously within the same model reliably, and are highly correlated. This is true across the full range of distributions (varying ), going well beyond the well-known Weibull case in which the two scale components are equivalent. The implication of this is that, in terms of performance gain, the movement from AFT to PH modelling (or vice versa) might not be very large, whereas we have found that modelling the shape is a more fruitful alteration to the regression specification.

Finally, the perspective of this paper has been to investigate survival modelling generally, to cover some of the most popular models, and to discover some of the better modelling choices that can be made within this framework. Although we have developed these ideas in a fully parametric context, non-parametric equivalents, while possible, are beyond the scope of the present paper (but are investigated in a separate line of work (Burke et al., 2018)). However, it is worth highlighting that perhaps too much emphasis is placed on non-/semi-parametric approaches in survival analysis whereby undue weight is attached to the flexibility of the baseline distribution in comparison to the flexibility of its regression structure. Our general approach to survival modelling provides a framework within which one can consider the most important components of survival modelling (including which might potentially be modelled non-parametrically), and we believe that this insight can lead to better modelling practice in general.

Acknowledgements

The first author would like to thank the Irish Research Council (www.research.ie) for supporting this work (New Foundations award).

References

• Bagdonaviçius and Nikulin (2002) Bagdonaviçius, V. and Nikulin, M. (2002). Accelerated Life Model; Modeling and Statistical Analysis. Chapman & Hall/CRC, Boca Raton, FL.
• Box and Cox (1964) Box, G. and Cox, D. (1964). An analysis of transformations (with discussion). Journal of the Royal Statistical Society Series B, 26:211–252.
• Burke et al. (2018) Burke, K., Eriksson, F., and Pipper, C. B. (2018). Semiparametric multi-parameter regression survival modelling. Submitted.
• Burke and MacKenzie (2017) Burke, K. and MacKenzie, G. (2017). Multi-parameter regression survival modelling: An alternative to proportional hazards. Biometrics, 73:678–686.
• Chao et al. (2012) Chao, C., Yubo, Z., Yingwei, P., and Jiajia, Z. (2012). smcure: fit semiparametric mixture cure models. R package version 2.0.
• Chen et al. (1999) Chen, M., Ibrahim, J., and Sinha, D. (1999). A new Bayesian model for survival data with a surviving fraction. J. Am. Statist. Ass., 94:909–919.
• Chen and Jewell (2001) Chen, Y. and Jewell, N. (2001). On a general class of semiparametric hazards regression models. Biometrika, 88:687–702.
• Chen (2000) Chen, Z. (2000). A new two-parameter lifetime distribution with bathtub shape or increasing failure rate function. Statistical Probability Letters, 49:155–161.
• Cox and Matheson (2014) Cox, C. and Matheson, M. (2014). A comparison of the generalized gamma and exponentiated Weibull distributions. Statistics in Medicine, 33:3772–3780.
• Dimitrakopoulou et al. (2007) Dimitrakopoulou, T., Adamidis, K., and Loukas, S. (2007). A lifetime distribution with an upside-down bathtub-shaped hazard function. IEEE Transactions on Reliability, 56:308–311.
• Ibrahim et al. (2001) Ibrahim, J., Chen, M., and Sinha, D. (2001). Bayesian semiparametric models for survival data with a cure fraction. Biometrics, 57:383–388.
• Jones and Noufaily (2015) Jones, M. and Noufaily, A. (2015). Log-location-scale-log-concave distributions for survival and reliability analysis. Electronic Journal of Statistics, 9:2732–2750.
• Jones et al. (2018) Jones, M., Noufaily, A., and Burke, K. (2018). A bivariate power generalized Weibull distribution: a flexible parametric model for survival analysis. Submitted.
• Kalbfleisch and Prentice (2002) Kalbfleisch, J. and Prentice, R. (2002). The Statistical Analysis of Failure Time Data. Wiley, Hoboken, NJ., 2 edition.
• Kirkwood et al. (1996) Kirkwood, J., Strawderman, M., Ernstoff, M., Smith, T., Borden, E., and Blum, R. (1996). Interferon alfa-2b adjuvant therapy of high-risk resected cutaneous melanoma: the Eastern Cooperative Oncology Group Trial EST 1684. J. Clinical Oncology, 14:7–17.
• Matheson et al. (2017) Matheson, M., Muñoz, A., and Cox, C. (2017). Describing the flexibility of the generalized gamma and related distributions. Journal of Statistical Distributions and Applications, 4:Article 15.
• Nadarajah and Haghighi (2011) Nadarajah, S. and Haghighi, F. (2011). An extension of the exponential distribution. Statistics, 45:543–558.
• Nikulin and Haghighi (2009) Nikulin, M. and Haghighi, F. (2009). On the power generalized Weibull family: model for cancer censored data. Metron, 67:75–86.
• Tsodikov et al. (2003) Tsodikov, A., Ibrahim, J., and Yakovlev, A. (2003). Estimating cure rates from survival data: an alternative to two-component mixture models. Journal of the American Statistical Association, 98:1063–1078.
• Xie et al. (2002) Xie, M., Tang, Y., and Goh, T. (2002). A modified Weibull extension with bathtub-shaped failure rate function. Reliability Engineering and System Safety, 73:279–285.
• Yeo and Johnson (2000) Yeo, I. and Johnson, R. (2000). A new family of power transformations to improve normality or symmetry. Biometrika, 87:954–959.
• Zeng and Lin (2007) Zeng, Z. and Lin, D. (2007). Maximum likelihood estimation in semiparametric regression models with censored data (with discussion). Journal of the Royal Statistical Society Series B, 69:507–564.

Appendix A The (A)PGW Distribution as a Transformation Model

Linear transformation models (c.f. Kalbfleisch and Prentice (2002, sec. 7.5); Zeng and Lin (2007)) are concerned with c.h.f.’s of the form

 HT(t)=w(θH0(t)) (A:1)

where depends log-linearly on covariates and both the transformation function and baseline function are c.h.f.’s. It is easy to see that the (A)PGW c.h.f.’s and are of the form (A:1). Either is a transformation model with the Weibull c.h.f. ; in terms of the overall model, whenever the Weibull is used as baseline, is the horizontal scale parameter so that when covariates are included only in it, the transformation model is an accelerated failure time model. The transformation is a version of the Box-Cox transformation given, in the simpler case of , by (Box and Cox, 1964; Yeo and Johnson, 2000).

If is the lifetime random variable following the transformation model with, for simplicity, c.h.f. , then the model can also be written in the form where follows the distribution with c.h.f. .

Appendix B Simulation Studies with Smaller Samples

This section displays analogous results to those of Table 4 from the main paper but for and respectively.