The XGTDL family of survival distributions

Non-PH parametric survival modelling is developed within the framework of the multiple logistic function. The family considered comprises three basic models: (a) a PH model, (b) an accelerated life model and (c) a model which is non-proportional hazards and non-accelerated life. The last model, the generalised time-dependent logistic model was described first by the author in 1996 and this model gives its name to the entire family. The family is generalised by means of a Gamma frailty extension which is shown to accommodate crossing hazards data. A further generalisation is the inclusion of a dispersion model. These extensions lead naturally to the concept of a multi-parameter regression model described by Burke and MacKenzie in which the scale and shape parameters are modelled simultaneously as functions of covariates. Where possible, we include the MPR extension in the XGTDL family. Following a simulation study, the new models are used to analyse two sets survival data and the methods are discussed.


Introduction
Prior to the development of the proportional hazards (PH) survival model, epidemiologists worked with the multiple logistic function in prospective longitudinal studies where the follow-up duration was fixed, e.g. in relation to disease incidence (Cox 1970). Not all incidence data conform to the PH assumption and MacKenzie (1996) introduced a family of survival models based on the multiple logistic function which specifically recognised this fact. The generalised time-dependent logistic family (GTDL) comprised three models (a) a PH model, (b) an accelerated life (AL) model and (c) a model which is non-PH and non-AL. Development was focussed on model (c), designated the GTDL model, and a Gamma extension demonstrated its ability to deal with crossing hazards survival data (MacKenzie and Ha 2007). This example involved the creation of a covariate-dependent shape parameter. Accordingly, the idea of modelling the shape parameter more generally intrudes and this has motivated the development of multi-parameter regression (MPR) survival models (Burke andMacKenzie 2013, 2017). In this paper, we extend the GTDL family to the XGTDL family ('X' for extended) by including: (a) Gamma frailty, (b) a dispersion model for the frailty variance, and where possible, (c) a MPR model to model non-PH survival data.

The GTDL family
The GTDL family of regression survival models was predicated on the multiple logistic function and the defining hazard functions of the three models in the family are as follows: (a) The logistic PH model (LPH) where π(·) = exp(·)/[1 + exp(·)] and the hazard function is λ(·) > 0 and the following scalers are defined as: α, γ 0 ∈ (−∞, +∞), and ϕ = exp(x β) > 0. There is no intercept term in the linear predictor in ϕ, in keeping with Cox's PH model (Cox 1972).
(c) The GTDL model (TDL when λ 0 = 1 ) where π(·) and α are defined above and λ 0 > 0 is a scalar. Here, γ = x β and x includes an intercept term x 0 . This model is neither PH nor AL.

Some remarks
Remark 1 Whilst we are principally interested in non-PH survival, model (a) is included for the sake of completeness. It has a very simple interpretation in which the baseline hazard is a linear logistic function of time and the hazard function is clearly unbounded above.
Remark 2 Model (b) is rather more complicated and it was developed specifically to avoid having a bounded hazard function (Al-tawarah and MacKenzie, 2003). It was derived by accelerating the baseline survivor function, which is obtained from equation (2), by first setting x = zero vector (whence ϕ = 1) and integrating to obtain the baseline integrated hazard. In particular, the scalar α parameter was introduced to allow for negative trends.

Remark 3
In many ways, model (c) is the most interesting, but its development was not straightforward. The constant λ 0 , introduced to avoid having a bounded hazard function, transpired to be inestimable in maximum likelihood estimation (MLE) analysis as it is aliased with the intercept term β 0 in the linear predictor γ . One approach was to drop the intercept term (MacKenzie 2002;Blagojevic et al. 2003), in which case all of the parameters are estimable in MLE. This model was termed the canonical form and designated CTDL. This idea seemed to work, but in fact it led to a more insidious problem, namely the lack of invariance of the model to the choice of reference subclass when fitting categorical covariates. Moreover, this problem was undetectable when the total sample size was less than 1000 and might well have been overlooked completely. Surprisingly, the solution to this problem (which does not occur with continuous covariates) is to permit λ 0 to vary, for example, with individuals as a Gamma frailty term. Prior to this solution, the GTDL model was used successfully with λ 0 = 1, when it was designated the TDL model. However, today, the Gamma frailty extension is preferred (Blagojevic et al. 2003;MacKenzie and Ha 2007;Ha and MacKenzie 2010).
Later, we extend the entire family to incorporate Gamma frailty in which case all the models will be non-PH.
(c) The GTDL model where γ = x β (including an intercept) is a scale parameter. Again when fitting this model we shall specialise λ 0 = 1, the TDL.

Cure fractions
It may be verified that the behaviour of the survival functions is determined by the value of α. Typically, for α > 0, S(0|x) = 1 and S(∞|x) = lim in other words, the survival function is proper, but for α < 0, S(0|x) = 1 and lim t−→∞ S(∞|x) = 0. Then, the survival function is improper, and in each case, we have a natural tail-deficit or cure rate model. We do not pursue the technical details of cure models in the XGTDL family here, but draw attention to related work by Milani et al. (2015), who studied the tail-deficit GTDL and the tail-deficit GTDL with Gamma frailty, while Rasouli et al. (2016) studied a GTDL model as a classical cure rate mixture model and Calsavara et al. (2020) analysed a melanoma dataset. 1

Frailty extension
Standard arguments involving an unobservable, multiplicative, random effect, u (> 0), on the hazard function yields the general formulae for the marginal survivor and hazard functions, viz. and We use these formulae to incorporate a Gamma frailty (GF) extension in all three members of the GTDL family in turn.

(b) LAL with Gamma frailty
The marginal survivor function is given by with corresponding marginal hazard function where, as before, the regression is held in ϕ = exp(x β).
(c) GTDL with Gamma frailty The marginal survivor function is given by with corresponding marginal hazard function where R(t; x) = [1 + exp(tα + γ )]/[1 + exp(γ )] and γ = x β. Here λ 0 does not appear as it was replaced by the frailty term, u, and integrated out (see Remark 3 in Sect. 2.1 for more explanation).
Thus, all three members of the GTDL family now incorporate frailty and are fully parametric, non-PH, survival models capable of further extension.

Frailty dispersion model
The use of Gamma frailty can be justified on the grounds of mathematical convenience (the probability density function of the marginal distribution can always be found analytically) and on the flexible range of shapes allowed by the Gamma distribution. However, it is sometimes argued that a single Gamma distribution does not afford sufficient flexibility. Accordingly, the flexibility can be extended by allowing the frailty variance to become a person-specific quantity, via another regression model, i.e.
Thus, the dispersion model (DM) is regression-based and hence driven by the data. When a DM is required, it relaxes the assumption that a single Gamma frailty distribution is sufficient. For example, if a binary covariate is supported by the data in the DM we shall have two GF distributions and when the frailty variance does not depend on covariates the model reduces to the standard frailty model viz: The concept of modelling the structure of the dispersion has been addressed in the joint mean-dispersion modelling literature (Smyth 1989;Lee and Nelder 2001;Pan and MacKenzie 2003), but its adoption in the frailty paradigm, in survival analysis, is more recent (Lynch and MacKenzie 2014). Furthermore, the use of frailty dispersion modelling in combination with an underlying GTDL family MPR model (below) is novel in the survival literature.
With the dispersion model in place we have in effect created the core of the XGTDL family. In the next section we extend the family further using the ideas of MPR modelling.

MPR modelling
Following Burke and MacKenzie (2017), the key idea is to model the scale and shape parameters simultaneously in terms of the covariates. Again, this idea is relatively new in the survival literature and is best illustrated by a familiar example. We start with the MPR Weibull (WB) model which is non-PH and used later as a comparator for the MPR models which we develop later in the XGTDL family.

The Weibull MPR model
The Weibull hazard function (Collett 2003) is λ(t) = λγ t γ −1 with (λ > 0 and > 0) whence the MPR form is simply: log e γ = z α (shape) and log e λ = x β (scale) resulting in the MPR hazard function . When x and z contain the same covariates, including intercepts, q = p. When α = (α 0 , α 1 , α 2 , . . . , α q ) = 0 the model reduces to an Exponential model (PH) with parameter ϕ. And when α = (α 0 , 0, 0 . . . , 0) the model is a standard Weibull model (PH) in which α 0 is the shape parameter. Inference in the MPR Weibull model is discussed in detail by Burke and MacKenzie (2017) and a wide range of two-parameter MPR survival models can be fitted using the mpr package developed by Burke (2016) in R (Core 2018).
With these arrangements in place, we can apply the MPR concept to existing models in the XGTDL family. It is worth noting that any model with more than one linear predictor is an MPR model, e.g. a Weibull regression model with a constant shape parameter and a Gamma frailty dispersion model.

A motivating example
We can gain some insight into the flexibility of MPR models by analysing the data provided by the Gastrointestinal Tumor Study Group (Schein and Lombard et al. 1982). This was a randomised controlled trial of the effects of two treatment regimens: chemotherapy vs combined chemotherapy and radiotherapy on the survival times of 90 (45 per regimen) gastric cancer patients. The resulting Kaplan and Meier (1958) survival curves are shown in Fig. 1. The survival curves cross indicating the presence of crossing hazards (violating the Cox (1972) PH assumption) and a significant difference in survival between groups (MacKenzie and Ha 2007). Early survival (0-2 years) is better in the chemotherapy group while longer term survival is better in the combined therapy group. Table 1 shows a number of different parametric models (one PH and the rest non-PH) fitted to the data. Inspecting the treatment effects (theβ 1 s) in Table 1, the TDL + MPR + GF model with Gamma frailty and separate shape parameters (theαs) shows a clear and statistically significant treatment effect. The WB + MPR + GF model shows a marginal result (t = −0.869/0.432 = −2.01) for the treatment effect, but the AIC for the TDL + MPR + GF model is much superior, suggesting that a non-PH base survival function is also required.
It also clear that Gamma frailty is important-the WB + GF model almost detects a difference in treatment effect (t = −1.006/0.508 = −1.98) but the TDL + GF result is unequivocal (t = −3.288/0.822 = −4.0). However, neither of these models can fit the crossing hazards data in Fig. 1. That requires an MPR model with separate tracks for each group.
Looking at the other models in Table 1 would lead one to conclude, erroneously, that there was no statistically significant difference between survival in the two treatment groups. This finding highlights the need for graphical confirmation of tableau results. Chemo. = lower-tailed curve (red) and Chemo.+ radiotherapy = higher tailed curve (black) Figure 2 shows the resulting goodness-of-fit of the TDL + MPR + GF model. The model fits the survival data in both treatment groups well and captures the crossing point of the survival functions. It is apparent that the goodness-of-fit lies, in large part, in modelling the shape parameters in each treatment arm. Accordingly, it is but a short step to consider that modelling the shape parameters may be a useful stratagem in general.

XGTDL-MPR GF models
We have already seen above that the standard two-parameter Weibull survival model can be extended to MPR form (Burke and MacKenzie 2017) and that the TDL + GF model can be extended (MacKenzie and Ha 2007). However, in general, it is still an open question as to whether all two-parameter survival models permit this extension. Accordingly, we turn now to systematically develop the MPR models in the XGTDL family. For convenience we assume that x contains the same covariates as z, whence q = p.

(a) The LPH + MPR + GF model
The LPH model with Gamma frailty is easily extended to MPR form by modelling the α parameter in (9). We choose where, in the shape parameter, α 0 can be negative, representing a declining hazard function (MacKenzie 1996). Note that the scale parameter ϕ is unchanged.
(b) The LAL + MPR + GF model It may be thought that the LAL model with Gamma frailty could be extended to MPR form in a similar fashion, but in this AL model the parameter ϕ = exp(β 1 x 1 + β 2 x 2 +, . . . , +β p x p ) ( shape and scale) influences both shape and scale. From (10), we see that the only other free parameter which could be given a linear predictor is λ, but this, too, influences the scale.
Attempts to implement such a MPR model produced unstable results and we omit this development.
(c) The TDL + MPR + GF model The TDL model with Gamma frailty is extended to MPR form as follows and λ 0 has been replaced by the frailty term and integrated out so we refer to this model as TDL rather than GTDL. Note also that the definition of γ is unchanged.
With these arrangements, the α parameter has become person-specific via the regression so that each subject has a separate covariate-dependent shape parameter. In this way, shape is modelled symmetrically with scale. Clearly, the same scheme applies equally to MPR models without Gamma frailty e.g. the LPH + MPR model.

Estimation and fitting
Model formulation may be quite detailed in the XGTDL family, but, fortunately, Maximum Likelihood (ML) estimation in parametric survival models is classical and straightforward. We assume that the observed survival time, t = min(t 1 , t 2 ), where t 1 is distributed as a XGTDL random variable, T 1 , and T 2 is a chosen censoring distribution such that T 1 is statistically independent of T 2 . For data (t i , d i , x i ), where i = 1, 2, . . . , n subjects and with independent censoring, the likelihood for the parameters, say θ , is where d i is the censoring indicator with d i = 1 for an event and d i = 0 for a censored observation. The hazard functions and survival functions for all of the models in the XGTDL family have been provided in foregoing sections. We fitted the models using a specially written R programme (Core 2018) which called the numerical optimisation procedure nlm to maximise the log-likelihood function, compute the maximum likelihood estimators and their standard errors from the observed information matrix.

Model selection
In the XGTDL family, there are three model classes, LPH, LAL and TDL. For the purpose of selecting amongst models within a class, standard information criteria may be used. We opted for the use of the Akaike (1974) information criterion, AIC = −2 log(θ) + 2 dim(θ), whereθ is the maximum likelihood estimator and dim(θ) is the number of parameters to be estimated. There are two levels of model selection, both of which can be handled by the AIC, namely: (a) covariate selection within a single model type (within a class) and (b) model type selection across classes.
In the three classes, there are 6 + 3 + 6 = 15 model types in all. To see this, let the base = LPH or LAL or TDL basic regression survival models. Then, for the LPH and TDL classes, we have: (1) base, (2) base + GF, (3) base + DM, (4) base + MPR, (5) base + MPR + GF, (6) base + MPR + DM. And for the LAL class, we have the first three model types only. To this, we can add another 6 models with base = WB (Weibull) as comparators making 21 models in all, and not space enough to present the findings for the whole family in the paper. Accordingly, we present only an example simulation.
The use of the AIC may pre-suppose that one is prepared, á priori, to fit several model classes. However, a great deal may be done, initially, by classical plotting methods: for example, by plotting model-based survival functions and/or the cumulative hazard functions against time and comparing with their non-parametric counterparts. In this way, the adequacy of the functional form of the base regression model may be assessed. Next, a check for possible shape effects would seem reasonable. Finally, the addition of gamma frailty might be considered on the grounds that the current regression may not be completely comprehensive. At each stage the AIC and the fit of the survival function can be checked for improvement.

Simulation study
We conducted the simulation with: two covariates x 1 (continuous , N (0, 1)) and x 2 (binary, 50% split); three sample sizes (n = 200, 500 and 1000); three censoring rates (10, 20 and 50%); three values of the frailty variance φ (0.5, 1 and 1.5) and finally for three different parameter vectors. Thus, there were 81 scenarios in principle for each model and the number of replications was m = 1000.
In the simulations, the proportion of right (random) censoring, p c , was controlled using an exponential distribution where the value of the controlling parameter,θ, was obtained from the J (·) function approach. Suppose the independent censoring times follow an exponential distribution with pdf g(t; ϑ) = ϑ exp(−ϑt). Let for fixed θ and a required censoring proportion, p c . Then, we findθ = arg min[J (ϑ)]. These minimizations were carried out in the R statistical package (Core (2018)) using the optimize procedure. The function S(t; θ) is, XGTDL family, model-specific.
We emphasise that the purpose of these simulations is simply to ensure that the parameters are recoverable in the scenarios studied and to ascertain the level of bias (if any). The true values of the parameters are shown at the bottom of the table. Table 2 shows the results for TDL + MPR + GF model. The bias is generally small and decreases with increasing sample size in the scale and shape parameters. The frailty variance is φ = exp(ψ 0 = −0.7) = 0.5. However, for the smallest sample size Table 2 Simulation: ML estimates, estimated standard errors (S.E.) and percentage bias for the TDL, MPR, Gamma frailty model with different sample sizes and censoring rates (n = 200), some estimates can be biased, but the bias is much reduced for increasing sample sizes and small when n = 1000. Of course, we expect that more complicated models tend to require larger sample sizes, but, overall, our simulation suggests that the model parameters are recoverable for reasonable sample sizes.
We also computed the probabilities that the intervals I j = [a j , b j ], for j = 1, . . . , m(= 1000) replications, covered the true value of the corresponding parameter, say θ j , where a j , b j =θ j ∓ 2 × se(θ j ), respectively. Setting c j =1 when I j covers θ j and c j = 0, otherwise, the estimated coverage probability for θ j is then Table 2 the coverage percentage lay between 92.5 and 97.5% (not shown).

Some remarks on model fitting
Remark 1. Throughout the simulation study (and subsequent data analysis), a userwritten R script was used which called the nlm procedure to fit the model. This is a numerical optimization routine with which we are familiar. It can use analytical first and second derivatives, but we allowed the routine to compute numerical derivatives and an approximation to the hessian. This approach has worked well for us in the past. As with any new model, it is important to be careful in the selection of the true values, to be simulated from, which must generate 0 < t < ∞ via the inverse distribution function transform (t = F −1 (u), where u ∼ Uniform(0, 1)). Similarly, with the selection of starting values which must generate a valid interior point in the parameter space. A poor choice of starting values may lead nlm to a false start with no progress from the starting values. This situation was trapped in the user-written R script. At the time of writing, we do not have experience of other methods of fitting these models.

Remark 2.
In the case of the TDL + MPR + GF models, the simulation study suggests that at least 20 observations per parameter are required for successful fitting and estimation of the corresponding variance covariance matrix. In general, studies of survival data with small sample sizes are unlikely to be of great interest. The minimum TDL + MPR + GF model with a single covariate in the scale and shape parameters requires 5 parameters. In the Gastric Tumour Study (1982), there were 90 observations (45 per group) and the 5 parameter TDL + MPR + GF model fitted easily. However, the censoring rates were low: 4.4% in one group and 13.2% in the other.
Remark 3. The distribution of the estimated parameters appeared as a uni-modal, symmetric, histogram-approximately Normal in shape, regardless of the sample size studied. Exceptions occurred when the simulation produced outliers due to nonconvergence of the numerical optimization routine. Such outliers were produced more often in scenarios where the sample size was small and the censoring rate high. When the outliers were removed, the Normal-like distribution re-emerged.
The simulation study included all of the models listed in Table 3 and the results for other models (not shown) follow a similar pattern.

Data analysed
We turn now to analyse the Lung Cancer dataset (Wilkinson 1995). This was a multisource, population-based, study of 855 incident cases diagnosed in Northern Ireland, UK, between October 1st, 1991 and September 30th, 1992. The patients were followed for almost 2 years and survival time was computed as the time from first diagnosis to death or censoring. Some 673 (77%) patients had died by the censoring date (30th May, 1993), leaving 182 (23%) censored. We adopt the conventional 5% level of statistical significance in what follows. Half the patients were dead by 6 months and more than half (51.6%) were allocated to the palliative care group at the time of diagnosis, indicating that no curative treatment was planned. Typically, such poor prognosis results from late diagnosis.

Analysis of treatment
The oncologists were principally interested in the effect of treatment and in which model(s) fitted the treatment survival data, as later, these would be adjusted for the effects of the other measured covariates.
The right panel of Fig. 3 shows the KM estimator by treatment category (Surgery, Chemotherapy, Radiotherapy, Chemo. + Radio., Palliative care). Surgical treatment has the best survival pattern and palliative care the worst. The observed survival tracks do not appear to follow a proportional hazards model-notice the late crossing hazards between chemotherapy and radiotherapy and the convergent behaviour of the combined therapy group.
As a first step, we illustrate fitting the models with only the treatment covariate. Table 3 shows the results of fitting all 21 models in the XGTDL family in terms of the type of parameters fitted, the log-likelihood, the AICs obtained, and the number of parameters fitted. Within a model class (e.g. Weibull) the various model types form a hierarchy in terms of the number of parameters fitted, which is the same across the classes (Weibull, LPH and TDL); LAL being the exception.
Inspection of Table 3 shows that including Gamma frailty in the models results in statistically significant decreases in the observed AICs throughout. For example, in the TDL class a decrease (dAIC) of 22.99 units is obtained between the basic Weibull model (3932.32) and the same model with Gamma frailty included (3888.34), for the addition of just 1 additional parameter. Increasing complexity, by including a dispersion model, also improves model fit throughout the table, although the gains, while still statistically significant, are more modest.
Moving to the MPR modelling scheme (without Gamma frailty) yields models which are close to the basic models with Gamma frailty. For example, in the TDL class the TDL + DM model has an AIC of 3887.20 while the AIC of the TDL + MPR model is 3888.95. In the table, only in the LPH class, is the LPH + MPR model superior to the LPH + DM (dAIC = 2.78). Both types model shape: the MPR model by person-specific covariates influencing the shape parameter and the DM model by Gamma frailty distributions stratified by covariates. That these are different measures of shape is evidenced in the table by the models which support Gamma frailty and MPR parameters. In fact the best model in 3 of the 4 classes is the MPR + GF model. Note that the addition of a dispersion model to a non-MPR model reduces the AICs by an significant quantity. But when the model is MPR, the reduction is not statistically significant. Overall, the conclusion is that the WB + MPR + GF model is a good working assumption.
In these data, the optimum AIC (= 3875.89) is provided by the WB + MPR + GF model which is substantially better than the AIC (= 3896.24) produced by the WB + MPR model which was used to analyse treatment in Burke and MacKenzie (2017). We note that none of the new models proposed fit the data as well as the WB + MPR + GF comparator. The nearest contender is the TDL + MPR + GF model with an AIC of 3883.70 which, however, is also significantly better than the WB + MPR model. We conclude that while the Gastric Cancer data required a non-PH base hazard function the Lung Cancer data do not.

Covariate adjustment
The NI Lung Cancer Study is an observational, population-based, study, and therefore, we adjust the treatment effects for the other measured covariates, using the models which performed well in the previous section, namely, the WB + MPR + GF and TDL + MPR + GF models.
The results are shown in Table 4 for fitting the full and reduced MPR + GF models in the Weibull and TDL classes. There are several interesting findings. First, the two models are re-assuringly broadly similar in relation to the magnitude and interpretation of the effects in the scale and shape regressions. Second, we note the large improvement in AIC when compared to any model involving treatment alone (cf. Table 4), indicating that non-treatment covariates are clearly important. Third, we note that more covariates influence the scale than the shape-twice as many in the Weibull case. This is a common finding in MPR modelling. On eliminating the redundant effects from the full models in Table 5, we reach dAICs of 23.36 and 20.42 for the Weibull and TDL models, respectively, indicating that the reduced models have significantly improved fit. We note, however, that in these reduced models the shape parameters cannot be ignored, dAIC = 9.80 for the Weibull and dAIC = 4.11 for the TDL, respectively, and thus the shape effects are less influential in the TDL + MPR + GF model. Table 3 AICs for XGTDL family and Weibull comparator models fitted using the treatment covariate Table 4 Comparison of Weibull and TDL (MPR + GF) models fitting all nine covariates  Table 5 compares unadjusted and adjusted treatment effects in the MPR + GF models. The adjustment is for the effect of the covariates (case-mix) included in the reduced models (see Table 4). The pattern is the same for both models, namely, adjustment for the effects of case-mix reduces the beneficial effects of surgery, chemotherapy and radiotherapy observed in the unadjusted analysis, but enhances the effect of the combined therapy. However, since only 34 cases received combined therapy the overall effect of adjustment is to reduce the benefit attributable to treatment. That the large benefit in the combined therapy group is an early effect, attenuated later as follow-up progresses, is evidenced by the large positive effect in the shape regression (coefficients not shown).

Discussion
The non-PH GTDL family is of course not the only way to model non-PH data, but it was the use of a member of this family in the crossing hazards example that led, in part, to the broader idea of modelling the shape parameter systematically with the usual scale parameter in survival distributions. A parallel idea was to routinely incorporate frailty because, realistically, survival data do not always conform to the proportional hazards assumption and important covariate information may be missed, even in randomised controlled trials using minimization. Accordingly, the extension to the XGTDL family is a natural statistical modelling development.
Whilst frailty modelling is now relatively well established, the idea of modelling shape symmetrically with scale (the MPR idea) is less well accepted in the survival literature. This paper demonstrates the potential benefits of adopting a MPR approach and corroborates our previous work in this area (Burke and MacKenzie 2017;Peng et al. 2020). Moreover, we have also demonstrated the importance of including both the frailty concept and the MPR approach, a finding which suggests that the MPR model captures a form of time dependence which classical frailty models cannot.
Although this is only one set of data, our experience with other data leads us to suggest that the addition of Gamma frailty and MPR regression to a base model will often be a good general starting place from which to analyse non-PH survival data.
When we embarked on this extension the question, as to whether all survival distributions could be extended easily to MPR form, was open. With the logistic accelerated life (LAL) class, we encountered problems when developing the MPR extension. The difficulty may result simply from the functional form of the hazard in AL models and we continue to investigate this topic in our laboratory. Meanwhile, we hope the models developed thus far in the XGTDL family will be of use to the survival modelling community.