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

TL;DR

This paper introduces multi-parameter regression (MPR) models for survival analysis, allowing covariates to influence multiple distributional parameters simultaneously, leading to more flexible models that relax proportional hazards assumptions.

Contribution

It proposes a novel MPR modeling approach in survival analysis, including a new variable selection strategy and implementation in an R package.

Findings

MPR models provide greater flexibility than traditional proportional hazards models.

Time-dependent hazard ratios emerge naturally from two-parameter Weibull models.

A new test for proportionality is motivated by MPR models.

Abstract

It is standard practice for covariates to enter a parametric model through a single distributional parameter of interest, for example, the scale parameter in many standard survival models. Indeed, the well-known proportional hazards model is of this kind. In this paper we discuss a more general approach whereby covariates enter the model through more than one distributional parameter simultaneously (e.g., scale and shape parameters). We refer to this practice as "multi-parameter regression" (MPR) modelling and explore its use in a survival analysis context. We find that multi-parameter regression leads to more flexible models which can offer greater insight into the underlying data generating process. To illustrate the concept, we consider the two-parameter Weibull model which leads to time-dependent hazard ratios, thus relaxing the typical proportional hazards assumption and motivating a new test of proportionality. A novel variable selection strategy is introduced for such multi-parameter regression models. It accounts for the correlation arising between the estimated regression coefficients in two or more linear predictors -- a feature which has not been considered by other authors in similar settings. The methods discussed have been implemented in the mpr package in R.

Tables

Table 1. Table 1: Single factor analysis results: standard errors are in brackets.

	MPR				PH
	Scale		Shape		Scale		Shape
Intercept	-1.28	(0.08)	-0.19	(0.04)	-1.48	(0.08)	-0.07	(0.03)
Palliative	0.00	——	0.00	——	0.00	——	——	——
Surgery	-3.91	(0.83)	0.59	(0.20)	-2.22	(0.23)	——	——
Chemo	-0.50	(0.32)	0.07	(0.14)	-0.40	(0.17)	——	——
Radio	-1.26	(0.19)	0.34	(0.07)	-0.57	(0.09)	——	——
C+R	-4.06	(0.88)	0.97	(0.16)	-0.87	(0.20)	——	——
$\ell \left(\hat{\theta }\right)$	-1938.1				-1960.8
AIC	-3896.2				-3933.5
BIC	-3943.8				-3962.0

Table 2. Table 2: Likelihood Ratio Tests and AIC differences

	L.R. Test (p-value)			${\Delta }_{\text{AIC}}$
	Scale	Shape	Joint	Scale	Shape	Joint
Treatment				32.4	14.8	60.4
Age Group	0.227	0.039	0.233	-2.3	2.1	-5.5
WHO Status		0.447		37.2	-4.3	69.2
Sex	0.835	0.627	0.851	-2.0	-1.8	-3.7
Smoker	0.684	0.214	0.037	-4.5	-1.5	1.4
Cell Type	0.019	0.835	0.002	4.0	-5.1	8.6
Metastases		0.056		29.6	1.8	50.8
Sodium	0.009	0.154	0.002	5.4	-0.3	9.4
Albumen		0.251		13.4	-1.2	18.5
Note: ${\Delta }_{\text{AIC}}={\text{AIC}}_{\text{reduced}}-{\text{AIC}}_{\text{full}}$.

Table 3. Table 3: Selection Procedure Results

	AIC Selection		BIC Selection
	PH	MPR	PH	MPR
Treatment	$\beta$	$\beta ,\alpha$	$\beta$	$\beta ,\alpha$
Age Group	—	—	—	—
WHO Status	$\beta$	$\beta$	$\beta$	$\beta$
Sex	—	—	—	—
Smoker	$\beta$	$\alpha$	—	—
Cell Type	$\beta$	$\beta$	$\beta$	$\beta$
Metastases	$\beta$	$\beta ,\alpha$	$\beta$	$\beta$
Sodium	$\beta$	$\beta$	$\beta$	$\beta$
Albumen	$\beta$	$\beta ,\alpha$	$\beta$	$\beta$
AIC	$3723.1$	$3679.7$	—	—
BIC	—	—	$3816.7$	$3802.4$
Note: $\beta$ = “selected in scale” and $\alpha$ = “selected in shape”

Equations

$$\lambda (t)$$

$$\log (\lambda )$$

$$\frac{\lambda (t|c=1)}{\lambda (t|c=0)}=\exp ({\beta }_c+{\alpha }_c)t^{\exp ({\tilde{z}}^{\mathrm{T}}\alpha )\{\exp ({\alpha }_c)-1\}},$$

$$[({\hat{\beta }}_c,{\hat{\alpha }}_c{)}^T-({\beta }_c,{\alpha }_c{)}^T{]}^T{\Sigma }_{{\hat{\beta }}_c,{\hat{\alpha }}_c}^{-1}[({\hat{\beta }}_c,{\hat{\alpha }}_c{)}^T-({\beta }_c,{\alpha }_c{)}^T]\sim {\chi }_2^2.$$

$$\log \lambda$$

$$\log \gamma$$

$${\hat{\pi }}_{x_k\in x}$$

$$\log (\lambda )$$

Full text

Keywords. Crossing hazards; Converging hazards; Diverging hazards; Multi-parameter regression; Non-proportional hazards; Survival analysis; Time-dependent effects; Variable selection.

1 Introduction

The standard approach to (parametric) regression is to relate covariates to one parameter of specific interest, for example, consider generalized linear models (McCullagh and Nelder,, 1989) where covariates enter through the location parameter while the dispersion parameter is simply a constant. We will refer to this standard approach as single parameter regression (SPR). The question we ask is why any one parameter should be deemed the “interest” parameter in terms of covariate analysis? This standard SPR approach neglects the potential influence of covariates on other distributional parameters which may also be important in describing the phenomenon under study. A more flexible approach is to allow covariates to enter the model through multiple distributional parameters, e.g., location and dispersion simultaneously. Hereafter, we refer to this approach as multi-parameter regression (MPR).

Multi-parameter regression modelling has appeared several times in the literature. In the context of normal linear regression, Park, (1966) and Harvey, (1976) modelled the dispersion parameter as a function of covariates to address heteroscedasticity, while Smyth, (1989) modelled dispersion in the more general case of generalized linear models. More recently, structured dispersion models have been considered by Lee and Nelder, (2001, 2006). The generalized additive models for location, scale and shape (GAMLSS) framework (Rigby and Stasinopoulos,, 2005; Stasinopoulos and Rigby,, 2007) goes beyond the exponential family to include a large variety of distributions in which covariates may enter through multiple parameters simultaneously; see also http://www.gamlss.org. The MPR approach has also been considered in longitudinal data analysis (joint mean-covariance models) (Pan and MacKenzie,, 2003, 2006, 2007) and in quantile regression (Noufaily and Jones,, 2013).

In the setting of survival analysis, fully parametric regression models have been much less popular than Cox’s (1972) semi-parametric proportional hazards (PH) model where the hazard function is given by $\lambda(t\,|\,x)=\exp(x^{\mathrm{\scriptscriptstyle T}}\beta)\lambda_{0}(t)$ where $x=(x_{1},\ldots,x_{p})^{\mathrm{\scriptscriptstyle T}}$ and $\beta=(\beta_{1},\ldots,\beta_{p})^{\mathrm{\scriptscriptstyle T}}$ are vectors of covariates and regression coefficients, respectively, and $\lambda_{0}(t)$ is a non-parametric function. Due to its wide acceptance of and availability in standard software, the Cox model is often fitted to data where the proportional hazards assumption does not hold (Schemper,, 1992). Although the Cox model can be generalized in various ways (McGilchrist and Aisbett,, 1991; Gray,, 1992; Klein,, 1992; Nielsen et al.,, 1992; Grambsch and Therneau,, 1994; Therneau and Grambsch,, 2000; Martinussen et al.,, 2002; Tian et al.,, 2005), our focus will be on fully parametric multi-parameter regression models.

We note that the basic Cox model is in fact a single parameter regression model since the scale of $\lambda(t\,|\,x)$ is modelled via the component $\exp(x^{\mathrm{\scriptscriptstyle T}}\beta)$ whereas its shape has no specific structure, being absorbed entirely by $\lambda_{0}(t)$. Although the stratified Cox model is somewhat analogous to a multi-parameter regression model in the sense that covariates enter the shape of $\lambda_{0}(t)$, this is limited to categorical covariates. Moreover, since $\lambda_{0}(t)$ is not modelled directly within the Cox model, we cannot estimate the effects of such stratification variables on the hazard. However, if instead we consider a fully parametric hazard, with parameters controlling its scale and shape, it is straightforward to develop a model where covariates enter through these scale and shape parameters. Indeed, many popular survival distributions have two such parameters, for example, Weibull, log-normal, gamma and log-logistic (Kalbfleisch and Prentice,, 2002; Lawless,, 2003). Within these two-parameter survival distributions, covariates are typically introduced via the scale parameter leading to proportional hazards (PH), accelerated failure time (AFT) or proportional odds (PO) models depending on the distribution. However, further flexibility could clearly be achieved by modelling the scale and shape parameters (i.e., multi-parameter regression), thereby accommodating a wider variety of survival data.

Multi-parameter regression is not wholly novel in a survival setting; early references include Taulbee, (1979) and Gaynor, (1987) who used polynomial and Gompertz hazard models respectively. Anderson, (1991) considered a log-linear model for the survival time where both the location and dispersion parameters depended on covariates (generalizing the AFT model) and, as an example, applied a Weibull MPR model (albeit different from the one in this paper) to the Framingham Heart Study data; interestingly, this model has been adopted in coronary heart disease literature where it is known as a “Framingham equation” (McEwan et al.,, 2004). More recently, the MPR approach has been used in the context of the inverse Gaussian model (Lee and Whitmore,, 2006; Aalen et al.,, 2008; Lee and Whitmore,, 2010). Furthermore, in their development of semi-parametric maximum likelihood estimation, Zeng and Lin, (2007) considered heteroscedastic transformation models for survival data encompassing multi-parameter regression. It is also noteworthy that cure models (Farewell,, 1982; Maller and Zhou,, 1995; Lambert et al.,, 2007; De Castro et al.,, 2010) fall within the MPR framework, i.e., covariates typically enter via the scale parameter and the cure probability.

Notwithstanding these developments, we note that multi-parameter regression in survival has not previously received much attention as a general procedure and is certainly not in mainstream use. We therefore explore and discuss the consequences of the multi-parameter regression approach using the popular Weibull model and show its flexibility in real data analysis. In particular, the MPR model provides a useful alternative to the proportional hazards assumption and produces a new test of deviation from proportionality. It may be fitted using our new mpr package (Burke,, 2016) in R which implements the techniques developed in this paper. See also the gamlss.cens package (Stasinopoulos et al.,, 2016) which extends the GAMLSS framework to survival data.

The paper is organized as follows. In §2 we introduce and develop the Weibull MPR model. Hypothesis testing and variable selection are discussed, including a simulation study to assess our proposed MPR selection procedure, in §3. An application of the methodology to lung cancer data appears in §4 and, finally, we conclude with some remarks in §5.

2 The Weibull multi-parameter regression model

The Weibull distribution is one of the most popular survival distributions and its hazard function is given by

[TABLE]

where $\lambda>0$, the scale parameter, controls the overall magnitude of the hazard and $\gamma>0$, the shape parameter, controls its time-evolution; it can increase ($\gamma>1$), decrease ($\gamma<1$) or remain constant ($\gamma=1$) over time. The Weibull MPR model is generated by introducing covariates into both parameters as follows:

[TABLE]

where the log-link is used to ensure positivity of both parameters, $x=(1,x_{1},\ldots,x_{p})^{\mathrm{\scriptscriptstyle T}}$ and $z=(1,z_{1},\ldots,z_{q})^{\mathrm{\scriptscriptstyle T}}$ are scale and shape covariate vectors which may or may not contain covariates in common and $\beta=(\beta_{0},\beta_{1},\ldots,\beta_{p})^{\mathrm{\scriptscriptstyle T}}$ and $\alpha=(\alpha_{0},\alpha_{1},\ldots,\alpha_{q})^{\mathrm{\scriptscriptstyle T}}$ are the corresponding regression coefficients. Maximum likelihood estimation of the regression coefficients is straightforward and has been implemented in the mpr package (Burke,, 2016).

Under the Weibull MPR model, the hazard ratio for a binary covariate, $c$, common to both scale and shape regression components, is

[TABLE]

where $\beta_{c}$ and $\alpha_{c}$ are the scale and shape coefficients of $c$ respectively and $\tilde{z}^{\mathrm{\scriptscriptstyle T}}\alpha=z^{\mathrm{\scriptscriptstyle T}}\alpha-c\,\alpha_{c}$ represents all of the other terms in the shape linear predictor excluding $c$. The hazard ratio is time-dependent where the time-evolution depends on the other shape-covariates via $\tilde{z}$. This dependence on $\tilde{z}$ can be dealt with by inserting a “typical” shape covariate profile or by integrating (2.3) over the empirical distribution of $\tilde{z}$ in the dataset (cf. Nelder and Lane, (1982); Karrison, (1987); Shen and Fleming, (1997); Martinussen and Pipper, (2013)). Most importantly, it is the sign of $\alpha_{c}$ which solely determines the path of the hazard ratio. When $\alpha_{c}=0$, the ratio of hazards reduces to the usual PH constant, $\exp(\beta_{c})$. Thus, the Weibull MPR model directly extends the proportional hazards model, providing a new test of proportionality adjusted for other scale and shape covariates. Furthermore, for $\alpha_{c}\neq 0$, setting (2.3) equal to one and solving yields an explicit solution for the time-point at which crossing hazards occur: $t_{c}=\exp\left\{(\beta_{c}+\alpha_{c})[\exp(\tilde{z}^{\mathrm{\scriptscriptstyle T}}\alpha)\{1-\exp(\alpha_{c})\}]^{-1}\right\}$.

We note that Anderson, (1991) considered a different Weibull MPR model. Anderson,’s model is based on an alternative hazard parametrization, $\lambda(t)=\lambda\gamma(\lambda t)^{\gamma-1}$, more suited to AFT modelling. One easily finds that it implies a hazard ratio whose functional form is somewhat more complicated than (2.3) and which does not reduce to $\exp(\beta_{c})$ when $\alpha_{c}=0$. Thus, the model presented here provides a more direct extension to the PH model. However, just as our model produces a PH test for a given covariate, Anderson,’s model produces an analogous AFT test – but this was not explored in Anderson, (1991). Furthermore, the primary model emphasized by Anderson, (1991) is a restricted and asymmetric parametrization in which the shape component is related to the scale component via $\log\gamma=\theta_{0}+\theta_{1}\log\lambda$ (although the unrestricted parametrization was also developed and applied further in O’Dell et al., (1994)). Such parametrizations do not permit variable selection in both regression components, a valuable procedure which enables one to determine which covariates have non-PH effects. The identification of such variables is important in practice (see Section 3). Finally, it is worth highlighting that, although both we and Anderson, (1991) consider the Weibull distribution as an example, the MPR approach is not limited to any one distribution and, indeed, the accompanying mpr package (Burke,, 2016) includes a variety of popular survival distributions.

3 Hypothesis testing and variable selection

We now consider the implications of correlation between estimated regression coefficients within MPR models on hypothesis testing and variable selection strategies. These issues do not appear to have been considered in depth by other authors.

3.1 Hypothesis testing

In standard single-parameter regression (SPR) survival models, covariates only appear in the scale and, therefore, the importance of a particular covariate, $c$, is determined by testing if its regression coefficient, $\beta_{c}$, differs significantly from zero. However, in the MPR model there are two hypotheses of interest, namely: (i) $H_{0}:\beta_{c}=0$ and (ii) $H_{0}:\alpha_{c}=0$. Furthermore, the estimated coefficients for the same covariate in different regression components (i.e., $\hat{\beta}_{c}$ and $\hat{\alpha}_{c}$) are typically quite correlated (see Appendix A). The consequence of this higher correlation is that the scale effect, $\hat{\beta}_{c}$, is most heavily adjusted for the shape effect, $\hat{\alpha}_{c}$, and vice versa. Thus, we recommend that covariates should be considered in both regression components simultaneously, since we can imagine scenarios where a covariate only becomes significant when it is present in both components.

It is also of interest to determine the total effect of $c$, i.e., by testing (iii) $H_{0}:\beta_{c}=\alpha_{c}=0$. Using the asymptotic result $(\hat{\beta}_{c},\hat{\alpha}_{c})^{T}\sim N[(\beta_{c},\alpha_{c})^{T},\Sigma_{\hat{\beta}_{c},\hat{\alpha}_{c}}]$, where $\Sigma_{\hat{\beta}_{c},\hat{\alpha}_{c}}$ is the relevant $2\times 2$ covariance sub-matrix of $\Sigma(\hat{\theta})$, we have that

[TABLE]

Thus, a p-value can be obtained by setting $(\beta_{c},\alpha_{c})^{T}=(0,0)^{T}$ in (3.4) and comparing this statistic to the $\chi^{2}_{2}$ distribution. Furthermore, a $(1-\alpha)100\%$ confidence ellipse for $(\beta_{c},\alpha_{c})^{T}$ is given by the set of $(\beta_{c},\alpha_{c})^{T}$ points defining the contour line such that (3.4) is equal to $\chi^{2}_{2,\,1-\alpha}$ (see Friendly et al., (2013) for a detailed account of confidence ellipses).

3.2 Selection procedure

We now define $M(x,z)$ to be a model with scale and shape covariate vectors $x$ and $z$ respectively. Hence, let $M_{0}=M(\tilde{x},\tilde{z})$ be the model where $c$ does not appear (i.e., $c\not\in\tilde{x},\tilde{z}$). Now, there are three models which include $c$: $M_{{\beta_{c}}}=M(\tilde{x}\,\cup\,c,\tilde{z})$, $M_{{\alpha_{c}}}=M(\tilde{x},\tilde{z}\,\cup\,c)$ and $M_{{\beta_{c}\alpha_{c}}}=M(\tilde{x}\,\cup\,c,\tilde{z}\,\cup\,c)$, respectively. Therefore, we can test hypotheses (i), (ii) and (iii) above by comparing $M_{\beta_{c}\alpha_{c}}$ with $M_{\alpha_{c}}$, $M_{\beta_{c}\alpha_{c}}$ with $M_{\beta_{c}}$ and $M_{\beta_{c}\alpha_{c}}$ with $M_{0}$, respectively. This can be done by means of likelihood ratio tests or information criteria (see Burnham and Anderson, (2002)); this approach is more in line with variable selection procedures (Miller,, 2002). Of course inferential problems associated with variable selection methods are well documented (Miller,, 1984; Hurvich and Tsai,, 1990; Zhang,, 1992) and automatic selection procedures are much criticized. Nonetheless, such procedures are useful when the number of covariates is large, so that methods for MPR models are required. Accordingly, an algorithm for MPR forward selection is described in Appendix B which also includes a comparison with the gamlss stepwise procedure.

3.3 Simulation study

We now evaluate the performance of a stagewise MPR variable selection procedure which involves starting from the null model, contains both forward and backward (scale, shape and simultaneous) steps and uses information criteria (AIC and BIC) as the basis of selection.

We simulated data from the Weibull MPR model with

[TABLE]

where $x=z=(1,x_{1},\ldots,x_{10})^{T}$ is a vector of independent binary variables. Thus, the first four covariates, $x_{1}$, $x_{2}$, $x_{3}$ and $x_{4}$, affect both the scale and the shape, $x_{5}$ and $x_{6}$ affect the shape only, $x_{7}$ and $x_{8}$ affect the scale only and, finally, $x_{9}$ and $x_{10}$ have no affect on either component. This setup gives a good variety of scale and shape effects of different strengths where the coefficient values were chosen to lead to realistic survival data. Furthermore, we varied the sample size ($n=$ 100, 500 and 1000) and censored proportion ($p=$ 20%, 50% and 80%) giving 9 scenarios in total, each of which was repeated 500 times.

At the $j$th repetition of a particular scenario, $j=1,\ldots,500$, the selection procedure was applied: let $x^{(j)}$ and $z^{(j)}$ denote the sets of scale and shape covariates selected in this repetition. Over these 500 simulation replicates, the relative (scale and shape) selection frequencies for the covariate $x_{k}$, $k=1,\ldots,10$, are given by

[TABLE]

respectively. The relative frequencies, for selection based on AIC and BIC, respectively, in each of the 9 simulation scenarios, are displayed in Fig. 1.

We first consider the results for selection based on AIC (Fig. 1, top panel). It is clearly difficult to determine which covariates affect survival when the information available is low, i.e., $(n,p)=(100,80\%)$. However, even when the sample size is small, the covariates with stronger effects (i.e., those with larger regression coefficients) are still selected 75% - 85% of the time provided that censoring is low ($p=20\%$). When a moderate level of information is available, $(n,p)=(500,50\%)$, the covariates with stronger effects are selected 100% of the time while those with weaker effects are selected over 80% of the time. Furthermore, as the sample size increases, the relative frequency of selection of any covariate with nonzero effect approaches one. Although important covariates were identified quite successfully, covariates with no effect were selected approximately 20% of the time across all scenarios using AIC. This can easily be combatted by using BIC (Fig. 1, bottom panel) but it comes at the expense of identifying important covariates less often than when using AIC (particularly in low information scenarios) which is not an unexpected finding.

4 Analysis of lung cancer data

Here we analyse a lung cancer dataset which appeared in a 1995 Queen’s University Belfast PhD thesis by P. Wilkinson and was further analysed by MacKenzie, (1996). The data contains individuals diagnosed with lung cancer in Northern Ireland during the period October 1st 1991 to September 30th 1992. Time was measured, in months, from date of diagnosis until the earlier of the occurrence of death or the study end date, which was May 30th 1993. In this observational study there were 855 individuals, of whom 673 died and the remainder were right-censored. A number of (categorical) covariates were measured on each individual, namely: treatment type, age group, WHO status, sex, smoking status, cancer cell type, metastases, and albumen and sodium levels in the blood.

We analyse the data using the Weibull MPR model from the mpr package. Single factor and multi-factor analyses are presented to illustrate both the capabilities of the MPR approach and the application of MPR variable selection. It is worth noting that a number of the categorical variables have missing values; in these cases, an extra category labelled “Missing” was created. Although this ad-hoc approach retains observations in the analysis, it can lead to bias in the estimated parameters (Hortan and Kleinman,, 2007; Jones,, 1996; Vach and Blettner,, 1991). However, addressing missingness is not our aim here.

4.1 Single factor analysis: treatment

We first investigate the unadjusted effect of treatment which has five levels, namely: palliative care (the reference category), surgery, chemotherapy, radiotherapy and chemotherapy + radiotherapy combined. Two models were fitted to the data: the multi-parameter regression Weibull model, and the standard, but more restrictive, shape-constant single parameter regression model, i.e., the proportional hazards Weibull model. These models (labelled MPR and PH) are summarized in Table 1.

In the PH model, the $\hat{\beta}$ coefficients are negative which means that treatment reduces the hazard (compared to palliative care). All treatments are statistically significant (at the 5% level) and, furthermore, the relative merit of each treatment is easily determined from the magnitude of the coefficients. However, the situation is more complex in the MPR model. The $\hat{\beta}$ coefficients are also negative in this case but the $\hat{\alpha}$ coefficients are positive. Thus, although the hazard is reduced through treatment, it is increasing relative to palliative care over time. In other words, the effectiveness of each treatment reduces over time. All treatments are statistically significant in both the scale and shape components apart from chemotherapy (neither component). This latter statement is based on considering the scale and shape effects separately; joint effects can be assessed using (3.4). To this end, joint confidence ellipses for the regression coefficients are shown in Fig. 2 along with the individual confidence intervals for comparison. Although the chemotherapy effect is non-significant, its confidence ellipse only just covers the $(0,0)$. Furthermore, while zero is well within the $\hat{\alpha}_{\text{C}}$ confidence interval, the $\hat{\beta}_{\text{C}}$ confidence interval only just includes zero. This motivates consideration of a reduced MPR model, fixing $\alpha_{\text{C}}=0$. Indeed, we find that this reduced MPR model (not shown) has improved AIC (3894.5) and BIC (3937.2), respectively and, interestingly, we then get $\hat{\beta}_{\text{C}}=-0.37~{}(\text{S.E. }=0.17)$, i.e., very close to the PH analysis, while all other coefficients are virtually unchanged from the original MPR model. This latter finding highlights that the MPR analyis permits the inclusion of simpler PH effects (as expected from Section 2), however, for the sake of brevity, we will not consider the reduced model further.

Although the $\beta$- and $\alpha$-coefficients in the MPR model provide some preliminary insights, the hazard ratios, (2.3), combine information from both of these components which allows us to determine the overall treatment effects. Thus, Fig. 3 shows the hazard ratio for each treatment relative to palliative care along with along with 95% confidence intervals calculated using the delta method. For comparison, the PH hazard ratios and confidence intervals are also shown. Based on the MPR analysis, the only treatment superior to palliative care over the full range of time is surgery. Chemotherapy appears to have a slight early effect in the first 4 months. Radiotherapy and the combined treatment are both superior to palliative care initially, but their effectiveness wears off after about 7 months. Although there appears to be some evidence that the combined treatment may become more hazardous than palliative care later in time, we note that there is much less information in the tails (15 - 20 months) and, indeed, this could also be a consequence of different frailties in the groups. Of course, being an observational study, these are not the true treatment effects since they depend on the net interplay of other factors, for example, many of the individuals selected for surgery were younger, healthier and free from metastases.

Figure 4 shows the Kaplan-Meier survivor curves (Kaplan and Meier,, 1958) for each treatment group along with the model-based (MPR and PH) predicted survivor curves calculated using the well-known relationship $S(t)=\exp(-\int_{0}^{t}\lambda(u)du)$ where $S(t)$ and $\lambda(t)$ are the survivor and hazard functions respectively. We can see that the MPR model fits the data more closely than the PH model. More formally, this improvement in fit is confirmed by carrying out a likelihood ratio test $(p\ll 0.001)$ and the fact that the MPR model has much lower AIC and BIC values. It is noteworthy that the regression coefficients, standard errors and predicted survivor curves arising from the semi-parametric Cox model (not shown) are almost identical to those of the PH Weibull model.

4.2 Multi-factor analysis: variable selection

We now consider a full covariate analysis of the lung cancer data. In order to initially determine the importance of each covariate, first the full model (scale and shape components saturated) was fitted to the data; this model had $\ell(\hat{\theta})=-1798.44$ and AIC $=3700.88$ respectively. Then, 27 reduced models were fitted which arise through excluding each of the nine covariates from scale, the shape and simultaneously from the scale and shape. These models were compared to the full model using likelihood ratio tests and AIC differences.

The results are shown in Table 2. Both the p-values arising from the likelihood ratio tests and the AIC differences agree, i.e., those with smaller p-values have positive $\Delta_{\text{AIC}}$ values and vice versa. It is interesting to note that although age is significant in the shape component, the overall (joint) effect of age is non-significant. Conversely, smoking status is not significant in either component, when judged individually, but the joint effect is significant. These findings are a consequence of the correlation that exists between scale and shape coefficients (discussed in Section 3 and Appendix A) and, furthermore, the need for considering scale, shape and joint effects is clear. The joint $\Delta_{\text{AIC}}$ values give us a sense of the relative importance of each covariate and, in particular, WHO status, treatment and the presence of metastases appear to be the most important factors. Cancer cell type, sodium level and albumen level are also significant, but to a lesser degree. Finally, smoking status has a weak effect whereas age and sex have little impact on survival.

Our stagewise selection procedure (based on both AIC and BIC with scale, shape and simultaneous steps) was applied to the full MPR model and, for comparison, the full PH model (which, of course, only had scale steps). The results are summarized in Table 3. Both the PH and MPR models are in agreement on the significant covariates where, in particular, age and sex are not selected; smoking status is also not selected using BIC. Although the PH and MPR models are in agreement on which covariate effects are significant, the nature of the effects will differ under the two models. Specifically, any covariate with an $\alpha$ effect in the MPR model is judged to have a non-PH effect. Thus, using AIC, treatment, smoking status, metastases and albumen level are judged to be non-PH whereas, using BIC, only treatment is judged to be non-PH. Thus, there is strong evidence that the effect of treatment is non-PH which is also supported by the large $\Delta_{\text{AIC}}$ value for its shape effect in Table 2. The adjusted treatment hazard ratios from this multi-factor MPR model (not shown) are quite similar to the unadjusted hazard ratios in Fig. 3 although the effects of surgery and radiotherapy are reduced (i.e., the hazard ratios are closer to one).

5 Discussion

Multi-parameter regression produces flexible parametric models which can, among other things, relax the proportional hazards assumption by allowing time-dependent hazard ratios. Furthermore, this is achieved without the need for specialized estimation procedures. Naturally, the process of variable selection is more involved in this multi-component setting due to correlation between estimated regression coefficients – an aspect which does not appear to have been considered in detail by other authors. However, our proposed selection procedure generalizes standard approaches (see Appendix B). Although we have focussed only on the two-parameter Weibull model in this paper, the use of alternative models is advisable in practice. In particular, more general parametric models can facilitate more general functional forms for the hazard ratio. Accordingly, the mpr package (Burke,, 2016) includes a variety of distributions. Finally, we believe that it is important for practitioners to diversify beyond standard Cox analyses and we suggest that multi-parameter regression models provide a flexible alternative.

Acknowledgements

We greatly appreciate the insightful comments made by the referees which allowed us to improve on an earlier draft of this paper. The first author would like to thank the Irish Research Council (www.research.ie) for supporting this work.

Appendix A: Correlation between estimated regression coefficients

In this simulation study we investigate the correlation structure for estimated regression coefficients. The data were simulated from the Weibull MPR model with

[TABLE]

where $x_{1}$ and $x_{2}$ are independent binary covariates. We fix regression coefficients and sample size ($n=1000$) here as the correlation structure is similar for other values. However, we varied the proportion censored at three different levels, 20%, 50% and 80%, and each of these simulation scenarios was repeated 500 times. At each repetition the Weibull MPR model was fitted to the simulated data giving a $500\times 6$ matrix of estimated coefficients for each of the three simulation scenarios.

Figure A.1 below shows scatter matrices of the estimated regression coefficients for both models in each of the three scenarios. We have not included the intercepts in these scatter matrices as we are primarily interested in coefficients of covariates. We can see that the “same-covariate-pairs”, $(\hat{\beta}_{1},\hat{\alpha}_{1})$ and $(\hat{\beta}_{2},\hat{\alpha}_{2})$, i.e., pairs corresponding to the same covariate, are highly correlated supporting the discussion in Section 3. All other coefficients are relatively uncorrelated, owing to the fact that $x_{1}$ and $x_{2}$ are independent in this simulation.

Appendix B: Multi-parameter regression selection algorithm

Using the notation of Section 3, we define $M(x,z)$ to be a model with scale and shape covariates $x=(1,x_{1},\ldots,x_{p})^{\mathrm{\scriptscriptstyle T}}$ and $z=(1,z_{1},\ldots,z_{q})^{\mathrm{\scriptscriptstyle T}}$ respectively. Figure B.1 below describes a forward selection procedure starting from the null model $M(x=\emptyset,z=\emptyset)$. At each iteration of the while loop the algorithm searches for a model which is better than the current model; “better” in this context may be based on likelihood ratio tests, Akaike information criterion, bayesian information criterion, etc. Of course, at each such iteration, there are candidate covariates to be considered for inclusion in the scale $\{c\,|\,c\not\in x\}$, the shape $\{c\,|\,c\not\in z\}$ and for simultaneous inclusion $\{c\,|\,c\not\in x,z\}$. Hence, this main loop is composed of three for loops where all of these new models are fitted and compared to the current model. If, at the end of the while loop, a better model is found, then this becomes the current model and another iteration of the loop commences. Otherwise, the loop ends as the best model has been found. A more general stagewise procedure has been implemented in the mpr package (the stepmpr function) which includes backward steps in addition to forward steps within the main while loop, i.e., there are three additional for loops carrying out backward steps in the scale, shape and simultaneously in the scale and shape.

We note that our selection procedure is more general than that which appears in the gamlss package. Within gamlss there are two strategies: “A” and “B”. Strategy A considers only one parameter at a time as follows: first apply forward selection to the scale, then (given the scale model) apply forward selection to the shape and, finally, return to the scale to apply backward selection. Strategy “B” in gamlss considers simultaneous steps only: apply forward and backward selection simultaneously to both the scale and shape. Our procedure generalizes both of these strategies as, at all points in our selection process, forward and backward steps are applied to both the scale and shape parameters individually and simultaneously.