Surviving and non surviving fraction regression models based on the beta modified Weibull distribution

Abstract

We propose two regression models based on the beta modified Weibull distribution. The first one is the long-term mixture lifetime regression applied to survival data when some individuals may never experience the event of interest for possible presence of long-term survivors in the data. This regression attempts to estimate the effects of covariates on the surviving fraction. The second one is regression model based on the log-beta modified Weibull distribution as an alternative to the log-modified Weibull regression. This model aims to estimate the effects of covariates on the survival times. These new models generalize some existing regressions in the literature. For both cases, the model parameters are estimated by the method of maximum likelihood for censored data. We derive the appropriate matrices for assessing the local influence on parameter estimates under different perturbation schemes and present a global sensitivity analysis. A model check based on the quantile residuals is performed to select the appropriate regressions. We reanalyze two data sets available in the literature, one for each regression.

Keywords

Beta modified Weibull distribution censored data local influence long-term survivors regression

1. Introduction

In the last decade, new classes of distributions have been proposed for modeling lifetime data by extending the Weibull distribution. One important extension is the beta modified Weibull (BMW) distribution (Silva et al., 2010). It includes as special cases some important distributions such as the generalized modified Weibull (GMW), beta Weibull (BW), exponentiated Weibull (EW), beta exponential (BE), modified Weibull (MW) and Weibull distributions, among several others. The BMW distribution has attracted several proposals in the area of reliability. The main motivation for using the BMW model is due to its flexibility in accommodating four types of the hazard rate function (hrf) (i.e. increasing, decreasing, unimodal and bathtub) depending on its parameters.

In this paper, we define two new regressions based on this distribution as feasible alternatives for modeling survival times that present or not cure fraction. Cure rate models have been used for modeling time-to-event data for various types of cancers, including breast cancer, non-Hodgkins lymphoma, leukemia, prostate cancer and melanoma.

Perhaps the most popular cure rate regressions is the mixture models (MMs) (Farewell, 1982). Li et al. (2001) investigated MMs in the presence of dependent censoring from the perspective of competing risks and model the dependence between the censoring and survival times using a class of Archimedean copula models. Zeng et al. (2006) studied a class of transformation models for survival data with a cure fraction. The class of transformation models is motivated by biological aspects, and it includes as special cases both the proportional hazards and proportional odds cure models. Further, Rizzato et al. (2009) introduced the generalized log-gamma mixture regression with covariates, Lanjoni et al. (2016) formula a new cure rate survival model, where the time to this event has the Burr XII geometric distribution and recently Ortega et al. (2017) proposed two new models with a cure rate called the odd Birnbaum-Saunders mixture and odd Birnbaum-Saunders geometric models and Ramires et al. (2018) propose a flexible semi-parametric cure rate survival model called the sinh Cauchy cure rate distribution. In this paper, we propose a new regression called the beta modified Weibull mixture (BMWM) regression for survival data with cure rate and long-term survivors. The proposed regression includes as special cases the traditional cure rate models.

In the absence of long term survivors, the location-scale regression model (Lawless, 2003) is distinguished since it is frequently used in clinical trials, when lifetimes are affected by explanatory variables. We define a log-linear regression called the log-beta modified Weibull (LBMW) regression for modeling data whit bathtub shaped, monotonically increasing or decreasing and upside-down bathtub. This new regression model can be an good alternative model for analysis survival data. The main motivation ao inves de motivação poderiamos usar vantagem for using the LBMW regression is that several models listed above, for example, the log-generalized modified Weibull and log-beta Weibull regressions, can be embedded in the LBMW regression. So, this model can be used too for select adequate model.

We consider a frequentist analysis for both regressions. The inferential part is carried out using the asymptotic distribution of the maximum likelihood estimators (MLEs). After modeling, it is important to check the assumptions in the regression and to conduct a robustness study to detect influential or extreme observations that can distort the results of the analysis. We discuss the diagnostic influence based on case-deletion introduced by Cook (1977). But, when case-deletion is used, all information from a single subject is deleted at once and therefore it is hard to tell whether that subject has some influence on a specific aspect of the model. A solution to this problem can be found in the local influence approach, where one investigates how the results of an analysis changed under small perturbations of the model, but where these perturbations can have specific interpretations. Cook (1986) introduced a general framework to detect the influence of observations which indicate how sensitive the analysis is when small perturbations are provoked in the data or in the model. We also used this methodology to detect influential subjects in the BMWM and LBMW regressions. Additionally, quantile residuals (Dunn & Smyth, 1996) are adopted to check classical assumptions in the regressions.

In recent years, several authors proposed various regression models with present or not cure fraction. See, for example, Korkmaz et al. (2019), Lanjoni et al. (2019), Pescim et al. (2019), Scudilio et al. (2019) and Yousof et al. (2019).

The rest of the paper is organized as follows. In Section 2, we review the BMW distribution and define the BMWM regression based on this distribution. In Section 3, we define the BMWM regression for censored data. We estimate the model parameters by the method of maximum likelihood and discuss some goodness-of-fit statistics in Section 4. Several diagnostic measures are addressed in Section 5 by considering case-deletion and normal curvatures of local influence under various perturbation schemes with censored observations. In Section 6, the quantile residuals (qrs) are used to assess departures from the underlying BMWM and LBMW regressions to detect outliers. In Section 7, two real data sets are analyzed, which prove empirically the flexibility, practical relevance and applicability of the proposed regressions. Section 8 ends with some concluding remarks.

2. The BMWM regression

2.1 The BMW distribution

Let $G(t)=1-\exp\{-\alpha t^{\gamma}\exp(\lambda t)\}$ be the cumulative distribution function (cdf) of the MW distribution with parameters $\alpha>0$ , $\gamma>0$ and $\lambda\geqslant 0$ . The survival function of the BMW distribution pioneered by Silva et al. (2010) can be expressed as (for $t>0$ )

$\displaystyle S(t)=1-I_{G(t)}(a,b)=1-\frac{1}{B(a,b)}\int_{0}^{G(t)}w^{a-1}(1-% w)^{b-1}dw,$ (1)

where $a,b,\lambda$ and $\gamma$ are shape parameters, $\alpha$ is a scale parameter, $B(a,b)$ is the complete beta function, $I_{x}(a,b)=B_{x}(a,b)/B(a,b)$ is the incomplete beta function ratio and $B_{x}(a,b)=\int_{0}^{x}w^{a-1}(1-w)^{b-1}dw$ is the incomplete beta function.

The density function and hazard rate function (hrf) corresponding to Eq. (1) are

$\displaystyle f(t)=\frac{\alpha t^{\gamma-1}(\gamma+\lambda t)\exp(\lambda t)}% {B(a,b)}[1-\exp\{-\alpha t^{\gamma}\exp(\lambda t)\}]^{a-1}\exp\{-b\alpha t^{% \gamma}\exp(\lambda t)\},$ (2)

and (by using $I_{z}(b,a)=1-I_{1-z}(a,b)$ )

$\displaystyle h(t)=\frac{\alpha t^{\gamma-1}(\gamma+\lambda t)\exp(\lambda t)}% {B(a,b)I_{\exp\{-\alpha t^{\gamma}\exp(\lambda t)\}}(b,a)}[1-\exp\{-\alpha t^{% \gamma}\exp(\lambda t)\}]^{a-1}\exp\{-b\alpha t^{\gamma}\exp(\lambda t)\}.$

Henceforth, $T\sim\text{BMW}(a,b,\alpha,\gamma,\lambda)$ denotes the random variable $T$ having the density Eq. (2).

According to Silva et al. (2010), a characteristic of the BMW distribution is that its hrf can be bathtub shaped, monotonically increasing or decreasing and upside-down bathtub depending basically on the parameter values. Another important characteristic of this distribution is that it contains as special cases several well-known distributions. For example, it simplifies to the BW distribution when $\lambda=0$ . If $\gamma=1$ , in addition to $\lambda=0$ , it reduces to the BE distribution. The GMW distribution is also a special case when $b=1$ . If $a=1$ in addition to $b=1$ , it gives the MW distribution. For $\lambda=0$ , $\gamma=b=1$ and $\lambda=0$ , the BMW model is equal to the EW and exponentiated exponential (EE) distributions, respectively. For $\gamma=2$ , $\lambda=0$ and $b=1$ , the BMW distribution is the generalized Rayleigh distribution. The Weibull distribution is clearly a basic exemplar when $a=b=1$ and $\lambda=0$ .

The $n$ th ordinary moment $\mu^{\prime}_{n}=E(T^{n})$ of $T$ can be expressed in terms of $a,b,\alpha,\gamma$ and $\lambda$ as (Silva et al., 2010)

$\displaystyle\mu^{\prime}_{n}=\sum_{j=0}^{\infty}\frac{(-1)^{j}\Gamma(a)}{B(a,% b)\Gamma(a-j)(b+j)j!}\sum_{i_{1},\ldots,i_{n}=1}^{\infty}\frac{A_{i_{1},\ldots% ,i_{n}}\Gamma(s_{n}/\gamma+1)}{[\alpha(b+j)]^{s_{n}/\gamma}},$

where $s_{n}=i_{1}+\ldots+i_{n}$ , $A_{i_{1},\ldots,i_{n}}=a_{i_{1}}\ldots a_{i_{n}}$ and

$\displaystyle a_{i}=\frac{(-1)^{i+1}i^{i-2}}{(i-1)!}\left(\frac{\lambda}{% \gamma}\right)^{i-1}.$

2.2 The new mixture regression

In a sample of censored survival times, the proportion of individuals who are immune and not subject to death, failure or relapse may be indicated by a relatively high number of individuals with large censored survival times. In this section, the BMW model is modified to become a new mixture cure rate model.

For formulating the BMWM regression, we consider that the studied population is a mixture of susceptible (uncured) individuals, who may experience the event of interest, and non-susceptible (cured) individuals, who will experience it (Maller & Zhou, 1996). This approach allows to estimate simultaneously whether the event of interest will occur, which is called incidence, and when it will occur, given that it can occur, which is called latency. Let $N_{i}$ , $i=1,\dots,n$ , be the indicator denoting that the $i$ th individual is susceptible ( $N_{i}=1$ ) or non-susceptible $(N_{i}=0)$ . The mixture model has the form

$\displaystyle S_{pop}(t_{i}|{\bf x}_{i})=\pi({\bf x}_{i})+\big{[}1-\pi({\bf x}% _{i})\big{]}S(t_{i}|N_{i}=1,{\bf x}_{i}),$ (3)

where $S_{pop}(t_{i}|{\bf x}_{i})$ is the unconditional survival function of $T_{i}$ for the entire population, $S(t_{i}|N_{i}=1,{\bf x}_{i})$ is the survival function for susceptible individuals given a covariate vector ${\bf x}_{i}^{\top}=(x_{i1},x_{i2},\ldots,x_{ip})$ and $\pi({\bf x}_{i})=P(N_{i}=0|{\bf x}_{i})$ is the probability of cure, which varies from individual to individual. We adopt a logistic link to the covariates, and then the probability that the $i$ th individual is cured can be expressed as

$\displaystyle\pi({\bf x}_{i})=\frac{\exp\big{(}{\bf x}_{i}^{\top}{\mbox{% \boldmath$\beta$}}\big{)}}{1+\exp\big{(}{\bf x}_{i}^{\top}{\mbox{\boldmath$% \beta$}}\big{)}},$ (4)

where ${\mbox{\boldmath$\beta$}}=(\beta_{1},\ldots,\beta_{p})^{\top}$ is the vector of unknown parameters. The BMWM regression follows by assuming that the survival function for the susceptible individuals in Eq. (3) is

$\displaystyle S(t_{i}|N_{i}=1,{\bf x}_{i})=I_{\exp\{-\alpha t_{i}^{\gamma}\exp% (\lambda t_{i})\}}(b,a).$

The BMWM regression when $\pi({\bf x}_{i})=0$ (for all ${\bf x}_{i}$ ) is identical to the BMW regression. For $a=b=1$ , we obtain the MW mixture (MWM) regression with long-term survivors and, for $a=b=1$ and $\lambda=0$ , we obtain the Weibull mixture (WM) regression with cure fraction.

In the BMWM regression, the parameters of interest are ${{\mbox{\boldmath$\theta$}}}=(a,b,\alpha,\gamma,\lambda,{\mbox{\boldmath$\beta% $}}^{\top})^{\top}$ , where $\beta$ indicates the long-term effects. The identifiability between the parameters in the cure fraction and those in the latency distribution for the mixture model has been discussed by Li et al. (2001). The MM is not identifiable when the cure fraction $\pi({\bf x})$ is a constant $\pi$ , but it is identifiable when $\pi({\bf x})$ is modeled by a logistic regression with covariates ${\bf x}$ . So, we have to include some covariates in the cure fraction to ensure identifiability. The estimate of the the proportion of cured individuals follows from Eq. (4) as

$\displaystyle\hat{\pi}({\bf x}_{i})=\frac{\exp\big{(}{\bf x}_{i}^{\top}% \widehat{{\mbox{\boldmath$\beta$}}}\big{)}}{1+\exp\big{(}{\bf x}_{i}^{\top}% \widehat{{\mbox{\boldmath$\beta$}}}\big{)}}.$

Then, the estimated mean cure fraction is $\hat{\pi}=\frac{1}{n}\sum_{i=1}^{n}\hat{\pi}({\bf x}_{i})$ .

2.3 The log-likelihood function

Consider data in the form $(t_{i},\mathbf{x}_{i})$ , $i=1,\dots,n$ , where $t_{i}$ denotes the observed time for the $i$ th subject, i.e., $t_{i}=\min\big{\{}T_{i},C_{i}\big{\}}$ , $T_{i}$ is the lifetime for the $i$ th individual, $C_{i}$ is the censoring time for the $i$ th individual and $\mathbf{x}_{i}$ is the covariate vector. Under these assumptions, the contribution of an individual that failed at $t_{i}$ to the likelihood function is

$\displaystyle\frac{\big{[}1-\pi({\bf x}_{i})\big{]}\alpha t_{i}^{\gamma-1}(% \gamma+\lambda t_{i})\exp(\lambda t_{i})}{B(a,b)}[1-\exp\{-\alpha t_{i}^{% \gamma}\exp(\lambda t_{i})\}]^{a-1}\exp\{-b\alpha t_{i}^{\gamma}\exp(\lambda t% _{i})\},$

and the contribution of an individual that is at risk at $t_{i}$ is

$\displaystyle\pi({\bf x}_{i})+\big{[}1-\pi({\bf x}_{i})\big{]}I_{\exp\{-\alpha t% _{i}^{\gamma}\exp(\lambda t_{i})\}}(b,a).$

Combining the last two equations, the log-likelihood function for the parameter vector ${\mbox{\boldmath$\theta$}}=(a,b,\alpha,\gamma,\lambda,{\mbox{\boldmath$\beta$}% }^{\top})^{\top}$ can be expressed as

$\displaystyle l({\mbox{\boldmath$\theta$}})=\sum_{i\in F}l_{1}\big{[}{\mbox{% \boldmath$\theta$}},t_{i},\pi({\bf x}_{i})\big{]}+\sum_{i\in C}l_{2}\big{[}{% \mbox{\boldmath$\theta$}},t_{i},\pi({\bf x}_{i})\big{]}.$ (5)

Here,

$\displaystyle l_{1}\big{[}{\mbox{\boldmath$\theta$}},t_{i},\pi({\bf x}_{i})% \big{]}=\log\left[\frac{\alpha}{B(a,b)}\right]+\log\big{[}1-\pi({\bf x}_{i})% \big{]}+(\gamma-1)\log(t_{i})+\log(\gamma+\lambda t_{i})+\lambda t_{i}+(a-1)% \log\{1-\exp[-\alpha t_{i}^{\gamma}\exp(\lambda t_{i})]\}-b\alpha t_{i}^{% \gamma}\exp(\lambda t_{i})$

and

$\displaystyle l_{2}\big{[}{\mbox{\boldmath$\theta$}},t_{i},\pi({\bf x}_{i})% \big{]}=\log\left\{\pi({\bf x}_{i})+\big{[}1-\pi({\bf x}_{i})\big{]}I_{\exp\{-% \alpha t_{i}^{\gamma}\exp(\lambda t_{i})\}}(b,a)\right\},$

where $\pi({\bf x}_{i})=\exp({\bf x}_{i}^{\top}{\mbox{\boldmath$\beta$}})/[1+\exp({% \bf x}_{i}^{\top}{\mbox{\boldmath$\beta$}})]$ and $F$ and $C$ denote the sets of individuals corresponding to lifetime and censoring times, respectively.

3. The LBMW regression for censored data

3.1 The log-regression model

We propose a log-regression model where the response variable $Y=\log(T)$ is related to the covariate vector ${\bf x}$ . It is thus important to study the distribution of $Y$ having the LBMW distribution.

Considering the transformations $\rm{\gamma}=\sigma^{-1}$ and $\rm{\alpha}={\rm e}^{-\mu/\sigma}$ , the density function of $Y$ is

$\displaystyle f(y;a,b,\lambda,\mu,\sigma)=\frac{u(\sigma^{-1}+\lambda{\rm e}^{% y})}{B(a,b)}(1-u)^{a-1}{\rm e}^{-bu},y\in\textbf{R},$ (6)

where $u=\exp[(y-\mu)/\sigma+\lambda{\rm e}^{y}]$ , $a>0$ , $b>0$ , $\lambda\geqslant 0$ , $\sigma>0$ , and $\mu\in\textbf{R}$ . We refer to Eq. (6) as the LBMW distribution. Thus,

$\displaystyle\text{if}\quad T\sim\text{BMW}(a,b,\lambda,\gamma,\alpha)\quad% \text{then}\quad Y=\log(T)\sim\text{LBMW}(a,b,\lambda,\mu,\sigma).$

The survival function of $Y$ can be expressed as

$\displaystyle S(y)=1-I_{G(y)}(a,b)=1-\frac{1}{B(a,b)}\int_{0}^{G(y)}w^{a-1}(1-% w)^{b-1}dw,$

where $G(y)=1-\exp\Big{[}-\exp\Big{(}\frac{y-\mu}{\sigma}\Big{)}\exp(\lambda\rm{e}^{y% })\Big{]}$ .

By expanding $Y^{s}=\log(T)^{s}$ in Taylor series around $\mu_{1}^{\prime}$ , the $s$ th moment of $Y$ can be expanded as

$\displaystyle E(Y^{s})=\log(\mu_{1}^{\prime})^{s}+\sum_{i=2}^{\infty}\frac{H^{% (i)}(\mu_{1}^{\prime})\mu_{i}}{i!},$

where $H^{(i)}(\mu_{1}^{\prime})$ is the $i$ th derivative of $H(\mu_{1}^{\prime})=\log(\mu_{1}^{\prime})^{s}$ with respect to $\mu_{1}^{\prime}$ and $\mu_{i}=E(T-\mu_{1}^{\prime})^{i}$ is the $i$ th central moment of $T$ .

The ordinary moments of $T$ can be given as an infinite sum of products of powers of two ordinary moments of $T$ , namely

$\displaystyle E(Y^{s})=\log(\mu_{1}^{\prime})^{r}+\sum_{i=2}^{\infty}\sum^{i}_% {k=0}\frac{(-1)^{k}H^{(i)}(\mu_{1}^{\prime})\mu_{i-k}^{\prime}\mu_{1}^{\prime k% }}{(i-k)!k!},$ (7)

where the moments $\mu_{i-k}^{\prime}$ and $\mu_{1}^{\prime}$ follow directly from the general formula for $\mu^{\prime}_{r}$ given in Section 2.1. Equation (7) is the main result of this section. The derivatives of $H(\mu_{1}^{\prime})=\log(\mu_{1}^{\prime})^{s}$ are easily found in Maple up to any order. Hence, the ordinary moments of the LBMW distribution are functions of the parameters $a,b,\lambda,\sigma$ and $\mu$ .

There are many real situations where the response variable $y_{i}$ depends on a vector of explanatory variables (or covariates) ${\bf x}_{i}^{T}=(x_{i1},x_{i2},\ldots,x_{ip})$ , and can be related to treat-ments, intrinsic traits of sample units, exogenous variables, interactions or time-dependent variables, among others.

Further, we consider that the location parameter $\mu$ of the LBMW regression depends on a matrix of explanatory variables ${\bf X}$ , say $\mu_{i}={\bf x}_{i}^{T}{\mbox{\boldmath$\beta$}}$ . In other words, a regression based on the LBMW distribution Eq. (6) linking the response $Y$ and the covariate vector ${\bf x}$ is

$\displaystyle y_{i}={\bf x}_{i}^{T}{\mbox{\boldmath$\beta$}}+\sigma z_{i},∼{}∼% {}∼{}i=1,\ldots,n,$ (8)

where ${\mbox{\boldmath$\beta$}}=(\beta_{1},\ldots,\beta_{p})^{T}$ and $\sigma>0$ are unknown parameters, ${\bf x}_{i}$ is the vector of the explanatory variables and $Z$ has the density function

$\displaystyle f(z)=\frac{[\sigma^{-1}+\lambda\exp(\mu+\sigma z)]u}{B(a,b)}(1-u% )^{a-1}\exp(-bu),z\in\textbf{R},$ (9)

where $u=u(z)=\exp[z+\lambda\exp(\mu+\sigma z)]$ , $a>0$ , $b>0$ , $\lambda>0$ , $\sigma>0$ , and $\mu\in\textbf{R}$ .

In this case, the survival function of $Y|{\bf x}$ is $S(y|{\bf x})=1-I_{G(y)}(a,b)$ , where $G(y)=1-\exp\{-\exp(\frac{y-{\bf x}^{T}{\mbox{\boldmath$\beta$}}}{\sigma})$ $\exp(\lambda\rm{e}^{y})\}$ .

The LBMW regression Eq. (8) opens new possibilities for fitting several different types of data. It contains as special cases the following well-known regressions and some other new regressions: for $a=b=\lambda=1$ , Eq. (8) is equal to the log-Weibull regression (Lawless, 2003); for $a=b=1$ , it is the same of the LMW regression (Carrasco et al., 2008); for $b=1$ , it reduces to the log-generalized modified Weibull regression (Ortega et al., 2011a); for $\lambda=0$ , it becomes the log-beta Weibull regression (Ortega et al., 2011b); for $\lambda=0$ and $b=1$ , it gives the log-exponentiated Weibull (LEW) regression.

For the interpretation of the estimated coefficients, a possible proposal is based on the ratio of median times (see Hosmer & Lemeshow, 1999). Hence, when the covariable is binary (1 or 0), and considering the ratio of median times with $x=1$ in the numerator, if $\hat{\beta}$ is negative (positive), it implies that the individuals with $x=1$ decrease (increase) the median survival time in $\exp(\hat{\beta})\times$ 100% as compared to those individuals in the group with $x=0$ assumed fixed the other covariables. This interpretation can be extended to continuous or categorical covariables.

3.2 The log-likelihood function

Consider the sample $(y_{1},{\bf x}_{1}),\ldots,(y_{n},{\bf x}_{n})$ of $n$ observations, where ${\bf x}_{i}$ is the explanatory variables vector associated with the $i$ th individual, $y_{i}=\min\big{[}\log(T_{i}),\log(C_{i})\big{]}$ , $T_{i}$ is the lifetime for the $i$ th individual and $C_{i}$ is the censoring for the $i$ th individual (for $i=1,\dots,n$ ). We assume no informative censoring and independence of the observed lifetime and censoring time. The log-likelihood function for the parameter vector ${\mbox{\boldmath$\theta$}}=(a,b,\lambda,\sigma,{\mbox{\boldmath$\beta$}}^{T})^% {T}$ in model Eq. (8) has the form

$\displaystyle l({\mbox{\boldmath$\theta$}})=\sum\limits_{i\in F}l_{1}({\mbox{% \boldmath$\theta$}},z_{i})+\sum\limits_{i\in C}l_{2}({\mbox{\boldmath$\theta$}% },z_{i}),$ (10)

where $F$ denotes the set of uncensored observations, $C$ denotes the set of censored observations,

$\displaystyle l_{1}({\mbox{\boldmath$\theta$}},z_{i})=z_{i}+\log(\sigma^{-1}+u% _{i})+u_{i}-\log[B(a,b)]+(a-1)\log\{1-\exp[-\exp(z_{i}+u_{i})]\}-b\exp(z_{i}+u% _{i}),$ $\displaystyle l_{2}({\mbox{\boldmath$\theta$}},z_{i})=\log\{I_{\exp[-\exp(z_{i% }+u_{i})]}(b,a)\},$

$z_{i}=(y_{i}-{\bf x}_{i}^{T}{\mbox{\boldmath$\beta$}})/\sigma$ and $u_{i}=\lambda\exp({\bf x}_{i}^{T}{\mbox{\boldmath$\beta$}}+\sigma z_{i})$ .

4. Inference

The MLE of the parameter vector $\theta$ can be determined by maximizing the log-likelihood function. We use the NLMixed procedure in the SAS software (see, Littel et al., 2006) or the MaxBFGS routine in the matrix programming language Ox (see, Doornik, 2007) to evaluate the estimate $\widehat{{\mbox{\boldmath$\theta$}}}$ . The regression models codes are available upon request from the interested readers. The inference procedure for $\theta$ can be based on the standard asymptotic likelihood techniques. Approximate confidence intervals and hypothesis tests on the model parameters can be constructed using the multivariate normal approximation given by $\widehat{{\mbox{\boldmath$\theta$}}}\sim N_{m}({\mbox{\boldmath$\theta$}},% \ddot{{\bf L}}({\mbox{\boldmath$\theta$}})^{-1})$ , where $m$ is number of parameters in the regression, $\ddot{{\bf L}}({\mbox{\boldmath$\theta$}})=-\partial^{2}l({\mbox{\boldmath$% \theta$}})/\partial{\mbox{\boldmath$\theta$}}\partial{\mbox{\boldmath$\theta$}% }^{T}$ is the $m\times m$ observed information matrix, whose elements can be computed numerically.

For testing nested models, the likelihood ratio test (LR) can be used to discriminate such models. We can compute the maximum values of the unrestricted and restricted log-likelihoods to construct the LR statistics for testing some sub-models. We consider the partition ${\mbox{\boldmath$\theta$}}=({\mbox{\boldmath$\theta$}}_{1}^{T},{\mbox{% \boldmath$\theta$}}_{2}^{T})^{T}$ , where ${\mbox{\boldmath$\theta$}}_{1}$ is a subset of the parameters of interest and ${\mbox{\boldmath$\theta$}}_{2}$ is a subset of the remaining parameters. The LR statistic for testing the null hypothesis $H_{0}:{\mbox{\boldmath$\theta$}}_{1}={\mbox{\boldmath$\theta$}}_{1}^{(0)}$ versus the alternative hypothesis $H_{1}:{\mbox{\boldmath$\theta$}}_{1}\neq{\mbox{\boldmath$\theta$}}_{1}^{(0)}$ is given by $w=2\{\ell(\widehat{{\mbox{\boldmath$\theta$}}})-\ell(\widetilde{{\mbox{% \boldmath$\theta$}}})\}$ , where $\widetilde{\mbox{\boldmath$\theta$}}$ and $\widehat{\mbox{\boldmath$\theta$}}$ are the MLEs under the null and alternative hypotheses, respectively. The statistic $w$ is asymptotically distributed as $\chi_{k}^{2}$ when $n\to\infty$ , where the degrees of freedom is equal to the difference between the numbers of parameters in the two models. For example, the test of $H_{0}:a=b=1$ versus $H:H_{0}$ is not true is equivalent to compare the BMWM (or LBMW) model with the WMM (or LWM) model.

Further, we can select the best model based on the values of the Akaike information criterion (AIC), corrected Akaike information criterion (AICc) and global deviance (GD) statistics. Given a data set, a researcher has an onerous task to select the best model among some possible alternatives. The model with the smallest AIC, AICc and GD values can be selected as the best one. We can also use the minimum value of the-log-likelihood function to choose the best fitted model.

5. Sensitivity analysis

After fitting the model, it is important to check its assumptions and to conduct a robustness study to detect influential or extreme observations that can cause distortions to the results of the analysis. The first tool to perform sensitivity analysis is the global influence starting from case deletion (see, Cook, 1977). Case deletion is a common approach to study the effect of dropping the $i$ th case from the data set (see, for example, Cook & Weisberg, 1982; Chatterjee & Hadi, 1988; Cook et al., 1988).

Another approach is suggested by Cook (1986), where instead of removing observations, weights are given to them. This is a local influence approach. In survival analysis, several authors have investigated the assessment of local influence as, for instance, Pettit and Bin Daud (1989) and Weissfeld (1990), who investigated local influence in proportional hazard regressions. Escobar and Meeker (1992) adapted local influence methods to regression analysis with censoring, Weissfeld and Schneider (1990) compared local influence and deletion case for the Weibull regression. Besides, there are various papers exploring new distributions such as the followings. Ortega et al. (2003) considered the problem of assessing local influence in generalized log-gamma regression with censored observations, Silva et al. (2008) addressed local influence in the log-Burr regression with censored data, Ortega et al. (2010) investigated local influence for the generalized log-gamma regression with cure fraction, Hashimoto et al. (2012) discussed local influence for the log-Burr XII regression for grouped survival data and Hashimoto et al. (2013) studied local influence in the new Neyman type A beta Weibull regression. Further, Fachini et al. (2014) worked with local influence for the bivariate regression with cure fraction and Ortega et al. (2014) derived the appropriate matrices for assessing local influence on the parameter estimates under different perturbation schemes for the odd Weibull regression.

5.1 Global influence

A first tool to perform sensitivity analysis, as stated before, is by means of global influence starting from case-deletion (see Cook, 1977; Cook et al., 1988). Case-deletion is a common approach to study the effect of dropping the $i$ th case from the data set.

In the following, a quantity with subscript “(i)” means the original quantity with the $i$ th case deleted. Let $\hat{{\mbox{\boldmath$\theta$}}}_{(i)}$ be the MLE of $\theta$ from $l_{(i)}({\mbox{\boldmath$\theta$}})$ . To assess the influence of the $i$ th case on the MLE $\hat{{\mbox{\boldmath$\theta$}}}$ , the basic idea is to compare the difference between $\hat{{\mbox{\boldmath$\theta$}}}_{(i)}$ and $\hat{{\mbox{\boldmath$\theta$}}}$ . If deletion of a case seriously influences the estimates, more attention should be paid to that case. Hence, if $\hat{{\mbox{\boldmath$\theta$}}}_{(i)}$ is far from $\hat{{\mbox{\boldmath$\theta$}}}$ , then the $i$ th case can be regarded as an influential observation. A first measure of the global influence is defined by the standardized norm of $\hat{{\mbox{\boldmath$\theta$}}}_{(i)}-\hat{{\mbox{\boldmath$\theta$}}}$ (generalized Cook distance)

$\displaystyle\textit{GD}_{i}({\mbox{\boldmath$\theta$}})=(\hat{{\mbox{% \boldmath$\theta$}}}_{(i)}-\hat{{\mbox{\boldmath$\theta$}}})^{T}\big{[}\ddot{{% \bf L}}({\mbox{\boldmath$\theta$}})\big{]}(\hat{{\mbox{\boldmath$\theta$}}}_{(% i)}-\hat{{\mbox{\boldmath$\theta$}}}),$

where $\ddot{{\bf L}}({\mbox{\boldmath$\theta$}})$ is the observed information matrix that can be computed numerically.

Another popular measure of the difference between $\hat{{\mbox{\boldmath$\theta$}}}_{(i)}$ and $\hat{{\mbox{\boldmath$\theta$}}}$ is the likelihood distance

$\displaystyle\textit{LD}_{i}({\mbox{\boldmath$\theta$}})=2\Big{\{}l(\hat{{% \mbox{\boldmath$\theta$}}})-l(\hat{{\mbox{\boldmath$\theta$}}}_{(i)})\Big{\}}.$

For the BMWM and LBMW regressiona, the log-likelihood functions are given in Sections 2.3 and 3.2, respectively. The index plots of GD and LD can be used to assess influence of the $i$ th observation.

5.2 Local influence

Another approach is suggested by Cook (1986), where instead of removing observations, weights are given to them. If likelihood displacement $\textit{LD}({\mbox{\boldmath$\omega$}})=2\{l(\hat{{\mbox{\boldmath$\theta$}}})% -l(\hat{{\mbox{\boldmath$\theta$}}}_{{\mbox{\boldmath$\omega$}}})\}$ is used, where $\hat{{\mbox{\boldmath$\theta$}}}_{{\mbox{\boldmath$\omega$}}}$ denotes the MLE under the perturbed model, the normal curvature for $\theta$ at direction ${\bf d},\parallel{\bf d}\parallel=1$ , is given by $C_{{\bf d}}({\mbox{\boldmath$\theta$}})=2|{\bf d}^{T}{\mbox{\boldmath$\Delta$}% }^{T}\big{[}\ddot{{\bf L}}({\mbox{\boldmath$\theta$}})\big{]}^{-1}{\mbox{% \boldmath$\Delta$}}{\bf d}|$ , where $\Delta$ is a $k\times n$ matrix that depends on the perturbation scheme, and whose elements are given by $\Delta_{ji}=\partial^{2}l({\mbox{\boldmath$\theta$}}|{\mbox{\boldmath$\omega$}% })/\partial\theta_{j}\partial\omega_{i}$ , for $i=1,2,\ldots,n$ and $j=1,2,\ldots,k$ , evaluated at $\hat{{\mbox{\boldmath$\theta$}}}$ and ${\mbox{\boldmath$\omega$}}_{0}$ , where ${\mbox{\boldmath$\omega$}}_{0}$ is the no perturbation vector. We can also determine normal curvatures for various index plots, for instance, the index plot of ${\bf d}_{\max}$ , the eigenvector corresponding to $C_{{\bf d}_{\max}}$ , the largest eigenvalue of the matrix ${\bf B}=-{\mbox{\boldmath$\Delta$}}^{T}\big{[}\ddot{{\bf L}}({\mbox{\boldmath$% \theta$}})\big{]}^{-1}{\mbox{\boldmath$\Delta$}}$ and the index plots of $C_{{\bf d}_{i}}$ , called the total local influence (see, for example, Lesaffre & Verbeke, 1998), where ${\bf d}_{i}$ denotes an $n\times 1$ vector of zeros with one at the $i$ th position. Thus, the curvature at direction ${\bf d}_{i}$ takes the form $C_{i}=2|{\mbox{\boldmath$\Delta$}}_{i}^{T}\big{[}\ddot{{\bf L}}({\mbox{% \boldmath$\theta$}})\big{]}^{-1}{\mbox{\boldmath$\Delta$}}_{i}|$ , where ${\mbox{\boldmath$\Delta$}}_{i}^{T}$ denotes the $i$ th row of $\Delta$ .

Next, for three perturbation schemes, we obtain the matrix

$\displaystyle{\mbox{\boldmath$\Delta$}}=({\mbox{\boldmath$\Delta$}}_{ji})_{(p+% 2)\times n}=\Bigg{(}\frac{\partial^{2}l({\mbox{\boldmath$\theta$}}|{\mbox{% \boldmath$\omega$}})}{\partial\theta_{i}{\mbox{\boldmath$\omega$}}_{j}}\Bigg{)% }_{k\times n},j=1,2,\ldots,k\text{ and }i=1,2,\ldots,n,$

by considering the log-likelihood Eqs (5) and (10) and the vector of weights ${\mbox{\boldmath$\omega$}}=(\omega_{1},\omega_{2},\ldots,\omega_{n})^{T}$ .

5.3 Case-weight perturbation

In this case, for the BMWM regression, the log-likelihood function has the form

$\displaystyle l({\mbox{\boldmath$\theta$}}|{\mbox{\boldmath$\omega$}})=\sum% \limits_{i\in F}\omega_{i}l_{1}\big{[}{\mbox{\boldmath$\theta$}},t_{i},\pi({% \bf x}_{i})\big{]}+\sum\limits_{i\in C}\omega_{i}l_{2}\big{[}{\mbox{\boldmath$% \theta$}},t_{i},\pi({\bf x}_{i})\big{]},$

and, for the LBMW regression, it becomes

$\displaystyle l({\mbox{\boldmath$\theta$}}|{\mbox{\boldmath$\omega$}})=\sum% \limits_{i\in F}\omega_{i}l_{1}({\mbox{\boldmath$\theta$}},z_{i})+\sum\limits_% {i\in C}\omega_{i}l_{2}({\mbox{\boldmath$\theta$}},z_{i}),$

where $0\leqslant\omega_{i}\leqslant 1$ and ${\mbox{\boldmath$\omega$}}_{0}=(1,\ldots,1)^{T}$ . The matrix $\Delta$ can be evaluated numerically.

5.4 Response perturbation

For the BMWM regression, each $t_{i}$ is perturbed as $t_{iw}=t_{i}+\omega_{i}S_{t}$ , where $S_{t}$ is a scale factor that can be estimated by the standard deviation of $T$ and $\omega_{i}\in{\bf R}$ . For the LBMW regression, each $y_{i}$ is perturbed as $y_{iw}=y_{i}+\omega_{i}S_{y}$ , where $S_{y}$ is a scale factor that can be estimated by the standard deviation of $Y$ . Here, ${\mbox{\boldmath$\omega$}}_{0}=(0,\ldots,0)^{T}$ . The matrix $\Delta$ is evaluated numerically.

5.5 Explanatory variable perturbation

Consider now an additive perturbation on a particular continuous explanatory va-riable, say $X_{t}$ , by setting $x_{it\omega}=x_{it}+\omega_{i}S_{t}$ , where $S_{t}$ is a scaled factor, $\omega_{i}\in{\bf R}$ . Here, ${\mbox{\boldmath$\omega$}}_{0}=(0,\ldots,0)^{T}$ . The matrix $\Delta$ is evaluated numerically.

6. Residual analysis

The assessment of the fitted regression is an important part of data analysis, particularly with censored data, and especially plots of the residuals play a central role in the checking of statistical regressions. By examining the residuals, we can detect the presence of outlying observations, the absence of components in the systematic part of the regression and departures from the error and variance assumptions.

We use randomized qrs (Dunn & Smyth, 1996) given by

$\displaystyle\widehat{r}_{i}=\Phi^{-1}(\hat{u}_{i}),$

where $\Phi^{-1}(u)$ is the standard normal quantile function and $\hat{u}_{i}$ is defined as if $y_{i}$ is an observation from a continuous response. For the proposed regressions, $\hat{u}_{i}$ is defined as follows:

•
The BMWM regression with cure rate

$\displaystyle\hat{u}_{i}=1-\log\left\{\hat{\pi}({\bf x}_{i})+\big{[}1-\hat{\pi% }({\bf x}_{i})\big{]}I_{\exp\{-\hat{\alpha}t_{i}^{\gamma}\exp(\hat{\lambda}t_{% i})\}}(\hat{a},\hat{b})\right\},$

where

$\displaystyle\hat{\pi}({\bf x}_{i})=\frac{\exp({\bf x}_{i}^{\top}\widehat{{% \mbox{\boldmath$\beta$}}})}{[1+\exp({\bf x}_{i}^{\top}\widehat{{\mbox{% \boldmath$\beta$}}})]}.$
•
The LBMW regression for censored data

$\displaystyle\hat{u}_{i}=I_{\widehat{G}(y_{i})}(\hat{a},\hat{b}),$

where $I_{\cdot}(\cdot)$ is the incomplete beta function ratio defined in Section 2.1 and

$\displaystyle\widehat{G}(y_{i})=1-\exp\Big{[}-\exp\Big{(}\frac{y_{i}-{\bf x}_{% i}^{\top}\widehat{{\mbox{\boldmath$\beta$}}}}{\hat{\sigma}}\Big{)}\exp(\hat{% \lambda}\rm{e}^{y_{i}})\Big{]}.$

For a right censored continuous response, $\hat{u}_{i}$ is defined as a random value from a uniform distribution on the interval $[\hat{u}_{i},1]$ . The true residuals $r_{i}$ ’s have a standard normal distribution if the model is correct (even when the model is not normal). The worm plots, introduced by van Buuren and Fredriks (2001), or QQ-plots of the normalized randomized qrs, can be used to test if the model is appropriate; see, for example, de Castro et al. (2010).
7. Applications

In this section, we prove the flexibility and applicability of the proposed regressions by means of two real data sets.

7.1 Application 1: The BMWM regression

In the first example, we fit the BMWM regression to a data set on cutaneous melanoma described and analyzed by Ibrahim et al. (2001) and Mizoi et al. (2007). The data set was from phase III of the melanoma clinical trial (E1690) conducted by Eastern Cooperative Oncology Group. The objective of the study was to compare a treatment with high-dose of interferon against observations from a control group. Relapse-free survival (years) was the outcome of interest, which was defined as the time from randomization to progression of tumor or death. There are $n=427$ patients on these combined treatment arms. After deleting observations with incomplete data and observed times equal to zero, we obtain a data set of $n=408$ patients with approximately 43% of censored cases. The following variables are associated with each individual (for $i=1,2,\dots,408$ ): $t_{i}$ : observed time (in years); $x_{i1}$ : treatment (0 $=$ control, 1 $=$ interferon dose); $x_{i2}$ : age (in years); $x_{i3}$ : nodule (nodule category: 1 to 4); $x_{i4}$ : sex (0 $=$ male, 1 $=$ female); $x_{i5}$ : p.s. (performance status-patient’s functional capacity scale as regards his/her daily activities: 0 $=$ fully active, 1 $=$ other); $x_{i6}$ : tumor (tumor thickness in mm).

Figure 1 displays the estimated survival curves for interferon group and the control group. An obvious plateau can be noted after about a 5 years’ follow-up, which offers empirical evidence for a cure possibility in E1690.

Figure 1.

Kaplan-meier curves of high-dose interferon and observation groups in E1690.

First, we consider the BMWM regression with long-term survivors given in Eq. (4) with all explanatory variables, i.e. ${\bf x}_{i}^{T}{\mbox{\boldmath$\beta$}}=\beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_% {i2}+\beta_{3}x_{i3}+\beta_{4}x_{i4}+\beta_{5}x_{i5}+\beta_{6}x_{i6}$ .

Table 1

MLEs for the parameters of the BMWM regression with long-term survivors fitted to the cutaneous melanoma data

Parameter	Estimate	Standard error	$p$ -value
$a$	4.5915	0.6744	–
$b$	0.6862	0.3052	–
$\alpha$	1.6001	0.4520	–
$\gamma$	0.7362	0.0789	–
$\lambda$	1E-8	1E-9	–
$\beta_{0}$	2.4161	0.6573	0.0003
$\beta_{1}$	$-$ 0.1599	0.2425	0.5100
$\beta_{2}$	$-$ 0.0159	0.0093	0.0891
$\beta_{3}$	$-$ 0.6026	0.1421	$<$ 0.0010
$\beta_{4}$	0.1991	0.2492	0.4248
$\beta_{5}$	$-$ 0.1358	0.3659	0.7107
$\beta_{6}$	$-$ 0.0767	0.0485	0.1144

We determine the MLEs of the parameters using the NLMixed procedure in SAS. Table 1 lists these estimates for the BMWM regression. The figures in Table 1 reveal that the variable nodule is significant (at 5% level) for the cure fraction. A summary of the AIC, AICc and GD values to compare the BMWM, MWM and WM regressions with long-term survivors is reported in Table 2. The BMWM regression gives the best fit according to these statistics.

Table 2

Some statistics for comparing regression models with long-term survivors

Model	BMWM	MWM	WM
AIC	1043.5	1047.6	1045.6
AICc	1044.2	1048.1	1046.0
GD	1019.5	1027.6	1027.6

An analysis under the BMWM regression with long-term survivors provides a check on the appropriateness of the MWM and WM sub-models and then indicates the extent for which inferences depend upon the regression. For example, the LR statistic for testing the hypotheses $H_{0}:a=b=1$ versus $H_{1}:H_{0}\text{is not true}$ , i.e. to compare the BMWM and MWM regressions, is $w=8.1$ ( $p$ -value $=$ 0.0174), which is favorable toward to the BMWM regression with long-term survivors.

The worm plots (wp function in R) and qqnorm of the qrs in Fig. 2 reveal that the BMWM regression is appropriate to these data.

Figure 2.

(a) worm plots and (b) qqnorm of the qrs for the cutaneous melanoma data.

In summary, we recommend the BMWM regression based on the analysis above. Thus, all subsequent analyzes are based on this regression. In order to detect possible outlying observations, we use the Ox computational program to compute the case-deletion measures $\textit{GD}_{i}({\mbox{\boldmath$\theta$}})$ and $\textit{LD}_{i}({\mbox{\boldmath$\theta$}})$ presented in Section 5.1. The influence measures index plots are displayed in Fig. 3. Based on these plots, the cases $\sharp$ 199 and $\sharp$ 341 are possible influential observations.

Table 3

MLEs for the BMWM regression fitted to the cutaneous melanoma data.

Parameter	Estimate	Standard Error	$p$ -value
$a$	4.1494	0.2819	–
$b$	0.7174	0.2388	–
$\alpha$	1.4825	0.4062	–
$\gamma$	0.7802	0.0403	–
$\lambda$	1E-8	1E-9	–
$\beta_{0}$	1.1141	0.2788	$<$ 0.001
$\beta_{3}$	$-$ 0.4890	0.1104	$<$ 0.001

Figure 3.

(a) index plot of $\textit{GD}_{i}({\mbox{\boldmath$\theta$}})$ for $\theta$ on the cutaneous melanoma data. (b) Index plot of $\textit{LD}_{i}({\mbox{\boldmath$\theta$}})$ for $\theta$ on the cutaneous melanoma data.

Next, we perform an analysis of local influence for the data using the BMWM regression. Local influence figures are omitted.

The cases $\sharp$ 279 and $\sharp$ 341 are very distinguished in relation to the others. All these cases are related to the uncensored observations. The cases $\sharp$ 279 and $\sharp$ 341 have the lowest times. The plots corresponding to the explanatory variables age and tumor are omitted since they indicate no influential observation. Then, the diagnostic analysis detected as potentially influential the following cases: $\sharp$ 199, $\sharp$ 279 and $\sharp$ 341. In order to reveal the impact of these observations on the parameter estimates, we refit the regression under some situations. First, we eliminate individually each one of the three cases. Second, we remove from the original data all the potentially influential observations. We conclude that these three observations do not cause a strong impact on the parameter estimates.

Further, we turn to a simplified model retaining only nodule category as an explanatory variable. The estimates for the BMWM regression with long-term survivors fitted to the cutaneous melanoma data are listed in Table 3, where the only significant variable is $x_{3}$ (nodule). Equation (4) gives the estimate of the proportion of cured individuals

$\displaystyle\hat{\pi}({\bf x}_{i})=\frac{\exp(1.1141-0.4890x_{i3})}{1+\exp(1.% 1141-0.4890x_{i3})}\quad\text{and}\quad\hat{\pi}=\frac{1}{417}\sum_{i=1}^{417}% \hat{\pi}({\bf x}_{i}).$

The estimated mean cure fraction is $\hat{\pi}=0.494$ .

In order to assess if the final fitted regression is appropriate, Fig. 4a displays the empirical survival function and the estimated marginal survival functions given by Eq. (3). Also, the estimated survival functions and the Kaplan-Meier estimate of the cure rate patients stratified by nodule category are displayed in Fig. 4b, from which we verify significant fractions of survivors for all nodule categories. It can be checked that the BMWM regression with long-term survivors provides a good fit to these data. Also, no observation appears as a possible outlier.

Table 4

MLEs of the parameters for the LBMW and LMW regressions and the corresponding SEs for the complete lung cancer data

	Beta modified Weibull			Modified Weibull
Parameters	Estimate	SE	$p$ -value	Estimate	SE	$p$ -value
$a$	6.5235	8.2885	0.4313	–	–	–
$b$	0.9285	1.2420	0.4547	–	–	–
$\lambda$	0.0004	0.0003	0.2029	0.0000	0.0002	–
$\sigma$	2.3655	1.5793	0.1342	0.9281	0.0615	–
$\beta_{0}$	2.7435	1.7814	0.1235	4.9769	0.2117	0.0000
$\beta_{1}$	$-$ 0.1335	0.2070	0.5188	$-$ 0.2285	0.1868	0.2213
$\beta_{2}$	0.0888	0.3479	0.7985	0.3977	0.2548	0.1185
$\beta_{3}$	$-$ 0.7163	0.3191	0.0248	$-$ 0.4285	0.2433	0.0782
$\beta_{4}$	$-$ 0.8545	0.3319	0.0100	$-$ 0.7351	0.2741	0.0073
$\beta_{5}$	0.0390	0.0070	0.0000	0.0301	0.0048	0.0000
$\beta_{6}$	$-$ 0.0020	0.0106	0.8541	$-$ 0.0005	0.0084	0.9553
$\beta_{7}$	0.0092	0.0103	0.3677	0.0061	0.0086	0.4758
$\beta_{8}$	$-$ 0.0062	0.0245	0.7996	$-$ 0.0044	0.0212	0.8362

Figure 4.

(a) kaplan-meier curves (solid lines), the estimated BMWM survival functions and the estimate of the cure fraction for the cutaneous melanoma data and (b) Kaplan-Meier curves (solid lines) and estimated BMWM survival functions stratified by nodule category (1–4, from top to bottom).

7.2 Application 2: The LBMW regression

As a second example, we consider the data from the Veterans Administration lung cancer trial given in Prentice (1973) and reported in Kalbfleisch & Prentice (2002). These data considered males with advanced inoperable lung cancer that received chemotherapy. The survival time is the time from the start of the treatment. The main purpose of the study was to compare the effects of two chemotherapy treatments in prolonging survival time. The explanatory variables in the study are: performance status at diagnosis ( $x_{5}$ ), a measure of general fitness on a scale from 0 to 100, the age in years of the patient ( $x_{6}$ ), the number of months from diagnosis of cancer ( $x_{7}$ ) and priori therapy ( $x_{8}$ ), 0: no priori therapy and 1: priori therapy. In addition, each patient was assigned one of two chemotherapy treatments (standard or test) and the tumors are classified into four types: large, adeno, small and squamous. The data contains $n=137$ observations of which 9 are censored.

Kalbfleisch and Prentice (2002) fitted a generalized F regression to analyze these data and used only the most important covariates: tumor type and performance status (PS). Lawless (2003) and Silva et al. (2009) used only one of the subgroup of patients- no priori therapy or priori therapy- and some covariates. Following Lawless (2003), we center only the explanatory variables $x_{5},x_{6}$ and $x_{7}$ and define the following regression (for $i=1,\ldots,137$ )

$\displaystyle y_{i}=\beta_{0}+\sum_{p=1}^{4}\beta_{p}x_{ip}+\beta_{5}(x_{i5}-% \bar{x}_{5})+\beta_{6}(x_{i6}-\bar{x}_{6})+\beta_{7}(x_{i7}-\bar{x}_{7})+\beta% _{8}x_{i8}+\sigma z_{i},$ (11)

where $y_{i}$ is the logarithm of the survival time $t_{i}$ , the random error $z_{i}$ is given by Eq. (9) and $x_{i1}=0$ if treatment is test, 1 otherwise, $x_{i2}=1$ if tumor type is squamous, 0 otherwise, $x_{i3}=1$ if tumor type is small, 0 otherwise and $x_{i4}=1$ if tumor type is adeno, 0 otherwise.

We fit the LBMW and LMW regressions to the current data. Table 4 gives the estimates of the parameters and their standard errors (SEs) for the two regressions. The required numerical evaluations are implemented in a Quasi-Newton algorithm through the MaxBFGS sub-routine in Ox (see, Doornik, 2007). The MLEs of the $\beta_{r}$ ’s and their SEs are different for the regressions. Further, the SEs of the estimates of the parameters for the LBMW regression are larger than those for the LMW regression. At the 5% level, the variables $x_{3}$ , $x_{4}$ and $x_{5}$ are significant for the LBMW regression, whereas the variables $x_{4}$ and $x_{5}$ are significant for the LMW regression.

Table 5

AIC, AICc and GD statistics for comparing the fitted LBMW and LMW regressions

Model	LBMW	LMW
AIC	410.5	414.3
AICc	413.5	416.4
GD	384.6	392.8

The plots comparing the empirical survival function and the fitted survival functions of the BMW and MW models are displayed in Fig. 5. These plots indicate that both corresponding regressions provide satisfactory fits. However, the BMW regression presents a better fit to these data.

Figure 5.

Empirical and estimated survival functions for the BMW and MW models for lung cancer data.

Further, we select the best fitted regression based on the AIC, AICc and GD values. A summary of these values to compare the LBMW and LMW regressions are given in Table 5. These results indicate that the LBMW regression is more appropriate for fitting these data. In addition, the LR statistic for testing the hypothesis $H_{0}:a=b=1$ (or the LMW regression is adequate) versus $H_{1}:H_{0}$ is not true (or the LBMW regression is adequate) is $w=7.7366$ ( $p$ -value $=$ 0.0209), and thus provides favorable support to the LBMW regression.

The worm plots (wp function in R) and qqnorm of the qrs in Fig. 6 for the fitted LBMW regression reveal that it is acceptable.

In summary, we recommend the LBMW regression based on the analysis above. Thus, all subsequent analyzes are based on this regression. In order to detect possible outlying observations, we use Ox to compute the case-deletion measures $\textit{GD}_{i}({\mbox{\boldmath$\theta$}})$ and $\textit{LD}_{i}({\mbox{\boldmath$\theta$}})$ . The influence measures index plots are displayed in Fig. 7. These plots indicate that the cases $\sharp$ 70 and $\sharp$ 85 are possible influential observations.

Table 6

MLEs for the LBMW regression fitted to the cancer data

Parameter	Estimate	SE	$p$ -value
$a$	6.1985	6.4081	–
$b$	0.90587	1.1425	–
$\lambda$	0.00045	0.00033	–
$\sigma$	2.3412	1.3396	–
$\beta_{0}$	2.7437	1.4537	0.0591
$\beta_{3}$	$-$ 0.7190	0.2446	0.0033
$\beta_{4}$	$-$ 0.9129	0.2861	0.0014
$\beta_{5}$	0.0390	0.0069	$<$ 0.0010

Figure 6.

(a) worm plots and (b) qqnorm of the qrs on the complete lung cancer data.

Figure 7.

(a) index plot of $\textit{GD}_{i}({\mbox{\boldmath$\theta$}})$ for $\theta$ on the cancer data. (b) Index plot of $\textit{LD}_{i}({\mbox{\boldmath$\theta$}})$ for $\theta$ on the cancer data.

Further, we perform an analysis of local influence for the data using the LBMW regression. Local influence figures are omitted. The observations $\sharp$ 44, $\sharp$ 70, $\sharp$ 75 and $\sharp$ 85 are very distinguished in relation to the others. All these cases are related to the uncensored observations. The cases $\sharp$ 70 and $\sharp$ 75 have longer times. On the other hand, the case $\sharp$ 85 has the lowest time and the case $\sharp$ 44 has the longest time when the performance $s t a t u s$ variable corresponds to the value 40. The plots corresponding to the explanatory variables diagnosis and age are omitted since they do not yield influential observations. Therefore, the diagnostic analysis detect as potentially influential the following four cases: $\sharp$ 44, $\sharp$ 70, $\sharp$ 75 and $\sharp$ 85. In order to reveal the impact of these observations on the parameter estimates, we refit the regression under two situations. First, we eliminate individually each one of these four cases. Second, we remove from the original data the totality of potentially influential observations. Hence, these results reveal that all cases do not cause a strong impact on the parameter estimates.

So, the final regression has the form

$\displaystyle y_{i}=\beta_{0}+\beta_{3}x_{i3}+\beta_{4}x_{i4}+\beta_{5}x_{i5}+% \sigma z_{i},∼{}i=1,2,\ldots,1207.$ (12)

The MLEs of the parameters in the final fitted regression are given in Table 6. In summary, we recommend the LBMW regression to fit these data on the basis of the above analysis. We interpret the estimated coefficients of the regression as follows: the expected survival time should increase approximately $4\%[\rm{e}^{0.0390}\times 100\%]$ as the performance status increases one unit, keeping the other variables fixed. In addition, the treatment does not appear to have sizeable effects, but the tumor type is important. Since the values of $\beta_{3}$ and $\beta_{4}$ are negative, the patients whose tumor type is adeno or small present survival smaller probabilities than do the patients with large tumor types.

8. Concluding remarks

In this paper, the log-beta modified Weibull and the beta modified Weibull mixture regressions with the presence of censored data are proposed as good alternatives to model survival lifetime data. We use the Quasi-Newton algorithm to calculate the maximum likelihood estimates. We adopt asymptotic likelihood ratio statistics for testing the model parameters. We discuss the local influence theory for the parameter estimates. So, we propose general attractive regressions for modeling censored and uncensored lifetime data. We prove empirically good adjustment of both regressions through sensitivity analysis in two applications to real data sets. We hope these generalizations may attract wider applications in survival analysis.

Footnotes

Acknowledgments

Special thanks to Fundação de Amparo à Pesquisa da Bahia- FAPESB and CNPq.

References

Buuren

, & Fredriks

(2001). Worm plot: A simple diagnostic device for modelling growth reference curves. Statistics in Medicine, 20, 1259-1277.

Carrasco

J. M. F.

Ortega

E. M. M.

, & Paula

G. A.

(2008). Log-modified weibull regression models with censored data: Sensitivity and residual analysis. Computational Statistics and Data Analysis, 52, 4021-4039.

Chatterjee

, & Hadi

A. S.

(1988). Sensitivity analysis in linear regression. New York: Wiley.

Cook

R. D.

(1977). Detection of influential observations in linear regression. Technometrics, 19, 15-18.

Cook

R. D.

(1986). Assesment of local influence (with discussion). Journal of the Royal Statistical Society, 48, 133-169.

Cook

R. D.

Peña

, & Weisberg

(1988). The likelihood displacement: A unifying principle for influence. Communications in Statistics: Theory and Methods, 17, 623-640.

Cook

R. D.

, & Weisberg

(1982). Residuals and influence in regression. New York: Chapman and Hill.

de Castro

Cancho

V. G.

, & Rodrigues

(2010). A hands-on approach for fitting long-term survival models under the GAMLSS framework. Computer Methods and Programs in Biomedicine, 97, 168-177.

Doornik

(2007). Ox 5: Object-oriented matrix programming language. London: Timberlake Consultants.

10.

Dunn

P. K.

, & Smyth

G. K.

(1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236-244.

11.

Escobar

L. A.

, & Meeker

W. Q.

(1992). Assessing influence in regression analysis with censored data. Biometrics, 48, 507-528.

12.

Fachini

J. B.

Ortega

E. M. M.

, & Cordeiro

G. M.

(2014). A bivariate regression model with cure fraction. Journal of Statistical Computation and Simulation, 84, 1580-1595.

13.

Farewell

V. T.

(1982). The use of mixture models for the analysis of survival data with long-term survivors. Biometrics, 38, 1041-1046.

14.

Ibrahim

J. G.

Chen

M. H.

, & Sinha

(2001). Bayesian survival analysis. Springer-Verlag: New York.

15.

Hashimoto

E. M.

Ortega

E. M. M.

Cordeiro

G. M.

, & Barreto

M. L.

(2012). The log-burr XII regression model for grouped survival data. Journal of Biopharmaceutical Statistics, 22, 141-159.

16.

Hashimoto

E. M.

Cordeiro

G. M.

, & Ortega

E. M. M.

(2013). The new neyman type a beta weibull model with long-term survivors. Computational Statistics, 28, 933-954.

17.

Hosmer

D. W.

, & Lemeshow

(1999). Applied survival analysis. New York: John Wiley and Sons.

18.

Kalbfleisch

J. D.

, & Prentice

R. L.

(2002). The Statistical Analysis of Failure Time Data. 2nd ed. New York: John Wiley.

19.

Korkmaz

M. Ç.

Cordeiro

G. M.

Yousof

H. M.

Pescim

R. R.

Afify

A. Z.

, & Nadarajah

(2019). The Weibull Marshall-Olkin family: Regression model and application to censored data. Communications in Statistics – Theory and Methods, 48, 4171-4194.

20.

Lanjoni

B. R.

Ortega

E. M. M.

, & Cordeiro

G. M.

(2016). Extended burr xii regression models: Theory and applications. Journal of Agricultural Biological and Environmental Statistics, 21, 203-224.

21.

Lawless

J. F.

(2003). Statistical models and methods for lifetime data. 2nd ed. New York: Wiley.

22.

Lesaffre

, & Verbeke

, (1998). Local influence in linear mixed models. Biometrics, 54, 570-582.

23.

C. S.

Taylor

J. M. G.

, & Sy

J. P.

(2001). Identifiability of cure models. Statistics and Probability Letters, 54, 389-395.

24.

Tiwari

R. C.

, & Guha

(2005). Mixture cure survival models with dependent censoring. Harvard University Biostatistics working Paper Series, 41, 211-224.

25.

Littell

R. C.

Milliken

G. A.

Stroup

W. W.

Wolfinger

R. D.

, & Schabernberger

(2006). SAS for Mixed Model (2nd edn). SAS Institute Inc.: Cary, NC, 2006.

26.

Maller

, & Zhou

(1996). Survival analysis with long-term survivors. New York: Wiley.

27.

Mizoi

M. F.

Bolfarine

, & Pedroso-De-Lima

A. C.

(2007). Cure rate model with measurement error. Communications in Statistics Simulation and Computation, 36, 185-196.

28.

Ortega

E. M. M.

Bolfarine

, & Paula

G. A.

(2003). Influence diagnostics in generalized log-gamma regression models. Computational Statistics and Data Analysis, 42, 165-186.

29.

Ortega

E. M. M.

Cancho

V. G.

, & Paula

G. A.

(2010). Generalized log-gamma regression models with cure fraction. Lifetime Data Analysis, 15, 79-106.

30.

Ortega

E. M. M.

Cordeiro

G. M.

, & Carrasco

J. M. F.

(2011a). The log-generalized modified weibull regression model. Brazilian Journal of Probability and Statistics, 25, 64-89.

31.

Ortega

E. M. M.

Cordeiro

G. M.

, & Hashimoto

M. E.

(2011b). A log-linear regression model for the beta-weibull distribution. Communications in Statistics. Simulation and Computation, 40, 1206-1235.

32.

Ortega

E. M. M.

Cordeiro

G. M.

Hashimoto

M. E.

, & Cooray

(2014). A log-linear regression model for the odd weibull distribution with censored data. Journal of Applied Statistics, 41, 1859-1880.

33.

Ortega

E. M. M.

Cordeiro

G. M.

Suzuki

A. K.

, & Ramires

T. G.

(2017). A new extended birnbaum-saunders model with cure fraction: Classical and Bayesian approach. Communications for Statistical Applications and Methods, 24, 397-419.

34.

Pescim

R. R.

Ortega

E. M. M.

Suzuki

A. K.

Cancho

V. G.

, & Cordeiro

G. M.

(2019). A new destructive poisson odd log-logistic generalized half-normal cure rate model. Communications in Statistics- Theory and Methods, 48, 2113-2128.

35.

Pettit

A. N.

, & Bin Daud

(1989). Case-weight measures of influence for proportional hazards regression. Applied Statistics, 38, 51-67.

36.

Prentice

R. L.

(1973). Exponential survival with censoring and explanatory variables. Biometrika, 60, 279-288.

37.

Ramires

T. G.

Hens

Cordeiro

G. M.

, & Ortega

E. M. M.

(2009). Estimating nonlinear effects in the presence of cure fraction using a semi-parametric regression model. Computational Statistics, 33, 709-730.

38.

Rizzato

F. B.

Ortega

E. M. M.

Demétrio

C. G. B.

(2009). The generalized log-gamma mixture model with covariates: local influence and residual analysis. Statistics Methods and Applications, 18, 305-331.

39.

Scudilio

Calsavara

V. F.

Rocha

Louzada

Tomazella

, & Rodrigues

A. S.

(2019). Defective models induced by gamma frailty term for survival data with cured fraction. Journal of Applied Statistics, 46, 484-507.

40.

Silva

G. O.

Ortega

E. M. M.

Garibay

V. C.

, & Barreto

M. L.

(2008). Log-Burr XII regression models with censored data. Computational Statistics and Data Analysis, 52, 3820-3842.

41.

Silva

G. O.

Ortega

E. M. M.

, & Cordeiro

G. M.

(2009). A log-extended Weibull regression model. Computational Statistics and Data Analysis, 53, 4482-4489.

42.

Silva

G. O.

Ortega

E. M. M.

, & Cordeiro

G. M.

(2010). The beta modified weibull distribution. Lifetime Data Analysis, 16, 409-430.

43.

Zeng

Yin

, & Ibrahim

J. G.

(2006). Semiparametric transformation models for survival data with a cure fraction. Journal of the American Statistical Association, 101, 670-684.

44.

Weissfeld

L. A.

(1990). Influence diagnostics for the proportional hazards models. Statistics and Probability Letters, 10, 411-417.

45.

Weissfeld

L. A.

, & Schneider

(1990). Influence diagnostics for the weibull model fit to censored data. Statistics and Probability Letters, 9, 67-73.

46.

Yousof

H. M.

Altun

Rasekhi

Alizadeh

Hamedani

G. G.

, & Ali

M. M.

(2019). A new lifetime model with regression models, characterizations and applications. Communications in Statistics – Simulation and Computation, 48, 264-286.