Quantile Regression with a New Exponentiated Odd Log-Logistic Weibull Distribution

Abstract

We define a new quantile regression model based on a reparameterized exponentiated odd log-logistic Weibull distribution, and obtain some of its structural properties. It includes as sub-models some known regression models that can be utilized in many areas. The maximum likelihood method is adopted to estimate the parameters, and several simulations are performed to study the finite sample properties of the maximum likelihood estimators. The applicability of the proposed regression model is well justified by means of a gastric carcinoma dataset.

1. Introduction

For survival data, the outcome variable is usually the time until the occurrence of an event of interest. A common characteristic of this type of data is the presence of censoring, that is, when the event of interest is not observed for some subjects before the study is finished. Furthermore, this variable depends on one or more explanatory variables (covariables), which have characteristics of the sample under study. Cox’s proportional hazards and accelerated failure time (AFT) models are two common tools in time-to-event modeling. The first class of models has the strong assumption of proportional risks, which is often invalid, so the effects of the covariables on the risk function are examined which can lead to difficult interpretations. The second class assumes that an association exists between the predictors and the survival time, permitting a direct interpretation of the effects of the covariables on lifetimes.

Nevertheless, these methods can fail to capture the heterogeneity of the effects of the covariables. In this respect, the quantile regression (QR) (Koenker and Bassett [1]) can be an alternative to these models, enabling evaluation of the heterogeneous effects of the predictors via analysis of different quantiles. This method involves modeling the quantiles of the survival time and links them to the covariables, providing some advantages, such as:

• Possible identification and inference under the heterogeneous effects of the covariables for different quantiles, thus supplying more complete information about the covariables and more flexibly controlling for the heterogeneity caused by them;
• Flexibility regarding the assumption of proportional risks;
• Provision of a direct interpretation of the results, that is, between the survival time and the covariables of interest;
• Possible analysis of different quantiles, allowing identification of the different effects of the covariables on individuals with different risks; and
• Robustness with respect to outliers in the regression models.

Originally, the QR methods are based on minimizing weighted absolute residuals [1] without any probability distribution, and the estimation of the parameters occurs by means of linear programming algorithms.

Although this approach is very flexible, some challenges such as: (i) the quantile crossing, that is, when two or more estimated quantile curves cross or overlap, causing difficulty in interpretability; and (ii) the drawback of the inability to apply parametric inference tools led to the search for other methods. Regarding the quantile crossing problem, we can verify alternative methods, such as: semiparametric models [2], the support vector (SV) regression approach [3], and a joint quantile estimation approach [4,5]. The weighted absolute residuals estimators coincide with the maximum likelihood estimators (MLEs), when the response follows a skewed Laplace distribution, so the initial association of a continuous distribution to the QR models was based on it (Koenker and Machado [6]).

In the context of censored data, an extensive bibliography can be mentioned, for example: Peng and Huang [7] developed a QR approach for survival data subject to conditionally independent censoring, Wang and Wang [8] proposed a locally weighted censored QR approach following the redistribution-of-mass idea and employed a local reweighing scheme. Zarean et al. [9] used the censored QR for determining overall survival and risk factors in esophageal cancer. Yang [10] presented a new approach for censored QR estimation, and Du et al. [11] developed estimation procedures for partially linear QR models, where some of the responses were censored by another random variable. Further, Xue et al. [12] addressed these limitations by using both simulated examples and data from National Wilms Tumor clinical trials to illustrate proper interpretation of the censored QR model and the differences and advantages of the model compared to the Cox proportional hazard model. Hong et al. [13] provided a practical guide for using QR for right-censored outcome data with covariates of low or high dimensionality, and De Backer et al. [14] studied a novel approach for the estimation of quantiles when facing potential right-censoring of the responses. Recently, De Backer et al. [15] investigated a new procedure for estimating a linear QR with possibly right-censored responses; Qiu et al. [16] considered the QR model for survival data with missing censoring indicators; Yazdani et al. [17] introduced the QR approach for modelling failure time and investigated the covariate effects for different quantiles; Peng [18] provided a comprehensive review of statistical methods for performing QR with different types of survival data; Hsu et al. [19] studied regression models for interval censored data using quantile coefficient functions via a set of parametric basis functions; He et al. [20] provided a unified analysis of the smoothed sequential estimator and its penalized counterpart for increasing dimensions in censored QR; and Wei [21] introduced a discussion about QR for censored data in haematopoietic cell transplant research. Note that all these articles cited in the QR with censored data did not use parametric models or use the skewed Laplace distribution (see [17]), whose estimators coincide.

Subsequently, other distributions were proposed by re-parameterizing them in terms of the quantile function (qf). Recent papers involving models for non-censored data based on other distributions can be mentioned: log-extended exponential-geometric [22]; Birnbaum–Saunders [23,24]; discrete generalized half-normal [25]; transmuted unit-Rayleigh [26]; unit-Burr-XII [27]; unit-Chen [28]; log-symmetric [29]; arcsecant hyperbolic Weibull [30]; and Dagum and Singh–Maddala [31] distributions. However, there is a relative lack in the literature of models for censored data in the parametric context: generalized Gompertz [32] and skew-t [33].

It is well known that the hazard rate function can assume different forms, which has led to the proposal of a large number of new distributions with the purpose of obtaining greater flexibility of data modeling, for example, Ref. [34]. In this sense, we introduce a QR regression model based on a reparameterized, exponentiated, odd log-logistic Weibull (EOLLW) distribution. It has two extra shape parameters, thus enabling the modeling of different forms of hazard rate functions, as well as data with positive or negative symmetric or asymmetric bimodal shapes, making it an alternative to the mixture models commonly used in the presence of bimodality. Another important feature of the new QR model is that it has as special cases: the exponentiated Weibull and odd log-logistic Weibull QR models. A detailed discussion of the theoretical foundations is given in analysis of survival data with concrete applications. The maximum likelihood method is adopted, and several simulations evaluate the behavior of these estimators under some scenarios. Additionally, we show that the model can establish functional relations of the covariables with other parameters, including scale and kurtosis, besides the quantile parameter.

The paper is structured as follows. Section 2 introduces a reparametrization of the EOLLW distribution based on quantiles. Section 3 addresses some mathematical properties. The proposed QR regression model, and some classic inference methods to estimate the parameters are addressed in Section 4. Some simulations are reported in Section 5. Section 6 provides a real application for the new regression model. Section 7 ends with a brief conclusion.

2. The Reparameterized EOLLW Distribution

Let $G(x;\mathbf{\eta })$ be a parent cumulative distribution function (cdf) and $g(x;\mathbf{\eta })=\mathrm{d}G(x;\mathbf{\eta })/\mathrm{d}x$ be its associated probability density function (pdf), both functions of a parameter vector $\mathbf{\eta }$. The cdf of the exponentiated odd log-logistic (EOLL-G) family is given by (Alizadeh et al. [35]) (for $x\in \mathbb{R}$)

$$\begin{array}{c}\hfill F(x;\nu ,\lambda ,\mathbf{\eta })=\frac{G{(x;\mathbf{\eta })}^{\nu \lambda }}{{\left\{G{(x;\mathbf{\eta })}^{\nu }+{[1-G(x;\mathbf{\eta })]}^{\nu }\right\}}^{\lambda }},\end{array}$$

(1)

where $\nu >0$ and $\lambda >0$ are two extra shape parameters.

The pdf corresponding to Equation (1) has the form

$$\begin{array}{c}\hfill f(x;\nu ,\lambda ,\mathbf{\eta })=\frac{\nu \lambda g(x;\mathbf{\eta })G{(x;\mathbf{\eta })}^{\nu \lambda -1}{[1-G(x;\mathbf{\eta })]}^{\nu -1}}{{\left\{G{(x;\mathbf{\eta })}^{\nu }+{[1-G(x;\mathbf{\eta })]}^{\nu }\right\}}^{\lambda +1}}.\end{array}$$

(2)

Henceforth, let $X\sim \mathrm{EOLL}\text{-}G(\nu ,\lambda ,\mathbf{\eta })$ be a random variable with density function (2).

The EOLL-G family reduces to the OLL-G class when $\lambda =1$ (Gleaton and Lynch [36]), and to the exponentiated (Exp-G) family (Mudholkar et al., [37]) when $\nu =1$. Clearly, it becomes the parent $G(x;\mathbf{\eta })$ when $\nu =\lambda =1$.

The EOLLW distribution is defined from (2) by taking the parent Weibull

$$G(x;\gamma ,\sigma )=1-exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\mathrm{and}g(x;\gamma ,\sigma )=\frac{\sigma }{{\gamma }^{\sigma }}{x}^{\sigma -1}exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right],x>0,$$

(3)

respectively, where $\mathbf{\eta }=(\gamma ,\sigma )$, $\gamma >0$ is a scale parameter, and $\sigma >0$ is a shape parameter.

The cdf of the random variable $X\sim \mathrm{EOLLW}(\nu ,\lambda ,\gamma ,\sigma )$ follows from Equations (1) and (3)

$$F(x;\nu ,\lambda ,\gamma ,\sigma )=\frac{{\left\{1-exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu \lambda }}{{\left[{\left\{1-exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu }+{\left\{exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu }\right]}^{\lambda }},x>0.$$

(4)

Based on Equations (2) and (3), the pdf of X becomes

$$f(x;\nu ,\lambda ,\gamma ,\sigma )=\frac{\nu \lambda \sigma x^{\sigma -1}{\left\{exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu }{\left\{1-exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu \lambda -1}}{{\gamma }^{\sigma }{\left[{\left\{1-exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu }+{\left\{exp\left[-{\left(\frac{x}{\gamma }\right)}^{\sigma }\right]\right\}}^{\nu }\right]}^{\lambda +1}},x>0.$$

(5)

The hazard rate function corresponding to (5) is $h(x;\nu ,\lambda ,\gamma ,\sigma )=f(x;\nu ,\lambda ,\gamma ,\sigma )/[1-F(x;\nu ,\lambda ,\gamma ,\sigma \left)\right]$.

By inverting (1), the qf of X reduces to

$$\begin{array}{c}\hfill x=Q\left(q\right)=Q_W\left\{\frac{q^{1/\left(\nu \lambda \right)}}{q^{1/\left(\nu \lambda \right)}+{(1-q^{1/\lambda })}^{1/\nu }}\right\},0<q<1,\end{array}$$

(6)

where $Q_W\left(q\right)=G^{-1}(p;\gamma ,\sigma )$ ($p\in (0,1)$) is the qf of the Weibull distribution, namely

$$\begin{array}{c}\hfill G^{-1}(p;\gamma ,\sigma )=\gamma {[-log(1-p)]}^{1/\sigma }.\end{array}$$

(7)

Thus, we rewrite the $\tau$th quantile (6) as

$$\begin{array}{c}\hfill x=Q\left(q\right)=\gamma {\left\{-log\left[\frac{{(1-q^{1/\lambda })}^{1/\nu }}{q^{1/\left(\nu \lambda \right)}+{(1-q^{1/\lambda })}^{1/\nu }}\right]\right\}}^{1/\sigma }.\end{array}$$

(8)

We can easily obtain the quartiles: first quartile (Q(0.25)), median (Q(0.5)), and third quartile (Q(0.75)).

We define a reparametrization of the pdf (5) as a function of the $\tau$th quantile (6), where the scale $\gamma$ becomes

$$\gamma =\mu {\left\{-log\left[\frac{{(1-{\tau }^{1/\lambda })}^{1/\nu }}{{\tau }^{1/\left(\nu \lambda \right)}+{(1-{\tau }^{\frac{1}{\lambda }})}^{1/\nu }}\right]\right\}}^{-1/\sigma },$$

(9)

$\mu >0$ is the location, and $\tau \in (0,1)$th is the quantile of X (assumed known).

By substituting (9) into Equation (4), the reparameterized cdf of X reduces to

$$F(x;\nu ,\lambda ,\mu ,\sigma )=\frac{{\left[1-exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu \lambda }}{{\left\{{\left[1-exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu }+{\left[exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu }\right\}}^{\lambda }},$$

(10)

where $w(\tau ,\lambda ,\nu )=-log\left\{{(1-{\tau }^{\frac{1}{\lambda }})}^{1/\nu }/\left[{\tau }^{1/\left(\nu \lambda \right)}+{(1-{\tau }^{1/\lambda })}^{1/\nu }\right]\right\}$.

By simple differentiation, the reparameterized pdf of X has the form

$$\begin{array}{c}\hfill f(x;\nu ,\lambda ,\mu ,\sigma )=\frac{\nu \lambda \sigma x^{\sigma -1}wexp{\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}}^{\nu }{\left[1-exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu \lambda -1}}{{\mu }^{\sigma }{\left\{{\left[1-exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu }+{\left[exp\left\{-w{\left(\frac{x}{\mu }\right)}^{\sigma }\right\}\right]}^{\nu }\right\}}^{\lambda +1}}.\end{array}$$

(11)

Henceforth, we redefine $X\sim \mathrm{EOLLW}(\nu ,\lambda ,\mu ,\sigma ,\tau )$ as a random variable with pdf (11), where $\tau \in (0,1)$ is fixed. Figure 1 displays plots of the pdf of X for some $\tau$ values, thus showing its asymmetry and bimodality.

The qf of X is obtained by replacing (9) in Equation (8)

$$x=Q\left(q\right)=\mu {\left[\begin{array}{c}log\left(\frac{{(1-q^{1/\lambda })}^{1/\nu }}{q^{1/\left(\nu \lambda \right)}+{(1-q^{1/\lambda })}^{1/\nu }}\right)/log\left(\frac{{(1-{\tau }^{1/\lambda })}^{1/\nu }}{{\tau }^{1/\left(\nu \lambda \right)}+{(1-{\tau }^{1/\lambda })}^{1/\nu }}\right)\end{array}\right]}^{1/\sigma },0<q<1.$$

(12)

3. Structural Properties

Some known properties of the reparameterized EOLLW distribution are given below:

(A)

Equation (11) gives ${lim}_{x\to \infty }f(x;\nu ,\lambda ,\gamma ,\sigma )=0$. Furthermore (Rodrigues et al. [38]),

$$\begin{array}{c}\hfill \underset{x\to 0^{+}}{lim}f(x;\nu ,\lambda ,\gamma ,\sigma )=\left\{\begin{array}{cc}\infty ,\hfill & \nu \lambda <1/\sigma ,\hfill \\ w/{\mu }^{\sigma },\hfill & \nu \lambda =1/\sigma ,\hfill \\ 0,\hfill & \nu \lambda >1/\sigma ,\hfill \end{array}\right.\end{array}$$

where $w(\tau ,\lambda ,\nu )$ is defined in Section 2.

(B)

The point $(x,f(x\left)\right)$ is called a critical point of f if x is in the domain of f, and either $f^{\prime }\left(x\right)=0$ or $f^{\prime }\left(x\right)$ do not exist. Since f in (11) is differentiable, $f^{\prime }\left(x\right)$ always exists in its domain. Differentiating $f\left(x\right)$ in (11) and setting it equal to zero, its critical points are the roots of (Rodrigues et al. [38])

$$\begin{array}{c}\hfill \frac{T^{″}\left(x\right)}{{\left[T^{\prime }\left(x\right)\right]}^2}=\frac{(\nu +1)T^{\nu }\left(x\right)+(1-\nu \lambda )}{T\left(x\right)[1+T^{\nu }\left(x\right)]},\end{array}$$

(13)

where $T\left(x\right)=G\left(x\right){[1-G\left(x\right)]}^{-1}$, $T^{\prime }\left(x\right)=g\left(x\right){[1-G\left(x\right)]}^{-2}$, $T^{″}\left(x\right)=g^{\prime }\left(x\right){[1-G\left(x\right)]}^{-2}+2g^2\left(x\right){[1-G\left(x\right)]}^{-3}$, and $g\left(x\right)$ and $G\left(x\right)$ are given in (3). Every critical point where f reaches a maximum (resp., minimum) value is called the mode (resp., minimum point). Using the intermediate value theorem, it is simple to verify that, for $\nu \lambda >1/\sigma$, Equation (13) has at least one root in $(0,\infty )$ (see Appendix A).

(C)

If $\sigma =1$ and $\nu >0$ is an integer, the pdf of X is (Rodrigues et al. [38])

• decreasing or decreasing–increasing–decreasing for $\nu \lambda <1$;
• unimodal for $\nu \lambda ⩾1$.

Note that Figure 1a ($\nu =2,\lambda =0.4$) shows the unimodality of the OELLW pdf when $\nu \lambda <1$ and $\sigma =4>1$. We emphasize that the values of the parameters in Figure 1a do not satisfy the hypothesis of the result C(2), so it does not contradict this one.

(D)

For $0<\nu ⩽1$, the pdf of X is (Rodrigues et al. [38])

• decreasing or decreasing–increasing–decreasing for $\nu \lambda <{\sigma }^{-1}$;
• decreasing or uni/bimodal or decreasing–increasing–decreasing for $\nu \lambda ={\sigma }^{-1}$.

(E)

If $X\sim \mathrm{EOLLW}(\nu ,\lambda ,\gamma ,\sigma )$ and $0<\nu ⩽1$, the pdf of X is uni- or bimodal for $\nu \lambda >{\sigma }^{-1}$ (Rodrigues et al. [38]).

Note that Figure 1b $(\nu =0.4,\lambda =0.90)$ and Figure 1c $(\nu =0.4,\lambda =1.5)$ are consistent with this result, because bimodality is obtained and $\nu \lambda >{\sigma }^{-1}$ is satisfied for both cases.

(F)

If D has the Type I Dagum distribution (Dagum [39]), say $D\sim \mathrm{DAGUM}(\nu ,1,\lambda )$, the cdf of $X\sim \mathrm{EOLL}\text{-}G(\nu ,\lambda ,\gamma ,\sigma )$ in (10) can be written as

$$\begin{array}{c}\hfill F(x;\nu ,\lambda ,\mu ,\sigma )=\mathbb{P}[D⩽T\left(x\right)]=\mathbb{P}\left[G^{-1}\left(\frac{D}{1+D};\mathbf{\eta }\right)⩽x\right].\end{array}$$

(14)

Consequently,

$$\begin{array}{c}\hfill X=G^{-1}\left(\frac{D}{1+D}\right)\end{array}$$

is a stochastic representation for X.

(G)

The cdf (10) satisfies the identity

$$\begin{array}{c}\hfill F\left(x\right)=\mathbb{P}\left[G^{-1}\left(\frac{1}{1+B}\right)⩾x\right],\end{array}$$

where $B=1/D$, $D\sim \mathrm{DAGUM}(\nu ,1,\lambda )$ and $G\left(x\right)$ is as given in Item (B).

(H)

We write $\mathcal{A}=\{(\nu ,\lambda )\in {(0,\infty )}^2:(1+\nu )t^{2\nu }+{\nu }^2(1+\lambda )t^{\nu +1}+[{\nu }^2(1+\lambda )+\nu (1-\lambda )+2]t^{\nu }+(1-\nu \lambda )>0,\forall t>0\}$ (see Rodrigues et al. [38]).

• If $(\nu ,\lambda )\in \mathcal{A}$ and $\sigma =1$, the hrf of X is increasing.
• Let ${\nu }^2(1+\lambda )+\nu (1-\lambda )+2⩾0$, $\nu \lambda >1$, and $\nu >0$ be an integer. For example, take the $\nu ⩾1$ integer and $\lambda >{\nu }^{-1}$.

  (a)

  If there exists $0<x^{*}<\nu /\gamma$ such that $h^{\prime }\left(x^{*}\right)=0$, the hrf of X has a bathtub (BT) shape.

  (b)

  If there does not exist $0<x^{*}<\nu /\gamma$ such that $h^{\prime }\left(x^{*}\right)=0$, the hrf of X is increasing.

• Let ${\nu }^2+3\nu -1>0$, $\nu \lambda =1$, and $\nu >0$ be an integer. For example, take $0<\nu <(\sqrt{13}-3)/2$ and $\lambda >n{u}^{-1}$. Under the conditions of Item (a) [Item (b)], the hrf of X has a BT (increasing) shape.

(I)

The EOLLW density transitions from heavy-tailed to light-tailed (Rodrigues et al. [38]).

3.1. Other Properties

3.1.1. Existence of Moments

The tail of the density of X follows from (14),

$$\begin{array}{c}\hfill \mathbb{P}(X⩾x)=\mathbb{P}[D⩾T(x\left)\right],D\sim \mathrm{DAGUM}(\nu ,1,\lambda ),\end{array}$$

where $T\left(x\right)=exp\left[{(x/\gamma )}^{\sigma }\right]-1$. Markov’s inequality gives

$$\begin{array}{c}\hfill \mathbb{P}[D⩾T\left(x\right)]⩽\frac{\mathbb{E}\left(D\right)}{T\left(x\right)}=\frac{\lambda \mathsf{\Gamma }(1-\frac{1}{\nu })\mathsf{\Gamma }(1+\frac{1}{\nu })}{T\left(x\right)},\nu >1.\end{array}$$

Hence,

$$\begin{array}{c}\hfill \mathbb{P}(X⩾x)⩽\frac{\lambda \mathsf{\Gamma }(1-\frac{1}{\nu })\mathsf{\Gamma }(1+\frac{1}{\nu })}{T\left(x\right)},\nu >1.\end{array}$$

(15)

Having an upper bound on the tail of the distribution, we proceed to bound the moments of X. This will prove its existence.

Based on the known formula $\mathbb{E}\left(X^p\right)=p{\int }_0^{\infty }{x}^{p-1}\mathbb{P}(X⩾x)\mathrm{d}x$ (for $X⩾0$ and $p>0$), the inequality holds:

$$\begin{array}{cc}\hfill \mathbb{E}\left(X^p\right)& \stackrel{\left(15\right)}{⩽}p\lambda \mathsf{\Gamma }\left(1-\frac{1}{\nu }\right)\mathsf{\Gamma }\left(1+\frac{1}{\nu }\right){\int }_0^{\infty }\frac{x^{p-1}}{T\left(x\right)}\mathrm{d}x,p>0,\nu >1,\hfill \\ \hfill & =p\lambda \mathsf{\Gamma }\left(1-\frac{1}{\nu }\right)\mathsf{\Gamma }\left(1+\frac{1}{\nu }\right)\frac{{\gamma }^p}{\sigma }{\int }_0^{\infty }\frac{y^{(p/\sigma )-1}}{exp\left(y\right)-1}\mathrm{d}y,\hfill \end{array}$$

where $y=exp\left[{(x/\gamma )}^{\sigma }\right]$.

From the Riemann zeta function $\zeta \left(s\right)={\sum }_{n=1}^{\infty }{n}^{-s}={\left[\mathsf{\Gamma }\left(s\right)\right]}^{-1}{\int }_0^{\infty }{z}^{s-1}{[exp\left(z\right)-1]}^{-1}\mathrm{d}z$, $\mathrm{Re}\left(s\right)>1,$ where $\mathsf{\Gamma }\left(s\right)={\int }_0^{\infty }{x}^{s-1}{\mathrm{e}}^{-\mathrm{x}}\mathrm{d}\mathrm{x}$ is the gamma function, we obtain

$$\begin{array}{c}\hfill \mathbb{E}\left(X^p\right)⩽\frac{p\lambda {\gamma }^p}{\sigma }\mathsf{\Gamma }\left(1-\frac{1}{\nu }\right)\mathsf{\Gamma }\left(1+\frac{1}{\nu }\right)\mathsf{\Gamma }\left(\frac{p}{\sigma }\right)\zeta \left(\frac{p}{\sigma }\right)<\infty ,p>\sigma .\end{array}$$

Thus, for $\nu >1$ and $p>\sigma$, the existence of the pth moment of X (for any $p>0$) is guaranteed.

3.1.2. Gini’s Mean Difference

Given the random variables $X_1,\cdots ,X_n$, the Gini mean difference (GMD) is defined as

$${\mathrm{GMD}}_n={\displaystyle \frac{1}{{\displaystyle \left(\genfrac{}{}{0pt}{}{n}{2}\right)}}}\sum _{1⩽i<j⩽n}\mathbb{E}\left(\right|X_i-X_j\left|\right),$$

(16)

provided the involved expectations exist. The GMD is a very useful measure of variability in the presence of non-normality.

(a)

If $X_1,\cdots ,X_n$ is a sequence of independent and identically distributed (iid) random variables, the classical GMD (La Haye and Zizler [40]) is $\mathrm{GMD}=\mathbb{E}\left(\right|X_1-X_2\left|\right).$

From Proposition 3 of Vila et al. [41], the GMD for a random sample $X_1,\cdots ,X_n$ of the EOLLW model is

$$\begin{array}{c}\hfill \mathrm{GMD}={\int }_0^1(2u-1)F_{X_1}^{-1}\left(u\right)\mathrm{d}u,\end{array}$$

(17)

where

$$F_{X_1}^{-1}\left(u\right)=\mu {\left[\begin{array}{c}\ln \left(\frac{{(1-u^{1/\lambda })}^{1/\nu }}{u^{1/\left(\nu \lambda \right)}+{(1-u^{1/\lambda })}^{1/\nu }}\right)/\ln \left(\frac{{(1-{\tau }^{1/\lambda })}^{1/\nu }}{{\tau }^{1/\left(\nu \lambda \right)}+{(1-{\tau }^{1/\lambda })}^{1/\nu }}\right)\end{array}\right]}^{1/\sigma }.$$

Note that analytically, the GMD (17) for the OELLW model is difficult to obtain. Vila et al. [41] provided the following upper bound $(2/\sqrt{3})\sqrt{\mathrm{Var}\left(X_1\right)}$ for the GMD.

(b)

If $X\sim \mathrm{EOLLW}({\nu }_i,{\lambda }_i,{\gamma }_i,{\sigma }_i)$, and $X_1,\cdots ,X_n$ is a sample not necessarily independent nor identically distributed, the following inequality for the GMD (16) follows from Vila et al. [41]

$$\begin{array}{c}\hfill {\mathrm{GMD}}_n⩽{\displaystyle \frac{1}{{\displaystyle \left(\genfrac{}{}{0pt}{}{n}{2}\right)}}}\sum _{1⩽i<j⩽n}\left[\sqrt{{(\sqrt{\mathrm{Var}\left(X_i\right)}-\sqrt{\mathrm{Var}\left(X_j\right)}{\rho }_{i,j})}^2+\mathrm{Var}\left(X_j\right)(1-{\rho }_{i,j}^2)}+\left|\varsigma \right|\right],\end{array}$$

where $\varsigma =\mathbb{E}\left(X_i\right)-\mathbb{E}\left(X_j\right)$ and ${\rho }_{i,j}=\mathrm{Corr}(X_i,X_j)$, for $i,j=1,\dots ,n$.

Under constraints ${\nu }_i>1$ and ${\sigma }_i<1$, the moments of $X_i\sim \mathrm{EOLLW}({\nu }_i,{\lambda }_i,{\gamma }_i,{\sigma }_i)$ ($i=1,\dots ,n$) always exist (see Section 3.1.1). Then, the mean $\mathbb{E}\left(X_i\right)$, variance $\mathrm{Var}\left(X_i\right)$, and correlations ${\rho }_{i,j}=\mathrm{Corr}(X_i,X_j)$ (for $i,j=1,\dots ,n$) also exist. Hence, for both cases (a) and (b), we can deduce non-trivial upper bounds (then its existence) of the GMD for the EOLLW model.

4. The EOLLW QR Model for Censored Data

A new regression model is defined from the reparametrized EOLLW density (11), and two systematic components for the parameters ${\mu }_i$ and ${\sigma }_i$ (for $i=1,\dots ,n$)

$$\begin{array}{c}\hfill {\mu }_i\left(\tau \right)=exp\left\{{\mathbf{v}}_i^{\top }{\mathbf{\beta }}_1\left(\tau \right)\right\}\mathrm{and}{\sigma }_i\left(\tau \right)=exp\left\{{\mathbf{v}}_i^{\top }{\mathbf{\beta }}_2\left(\tau \right)\right\},\end{array}$$

(18)

where ${\mathbf{\beta }}_1\left(\tau \right)={({\beta }_{10},{\beta }_{11},\cdots ,{\beta }_{1p})}^{\top }$ and ${\mathbf{\beta }}_2\left(\tau \right)={({\beta }_{20},{\beta }_{21},\cdots ,{\beta }_{2p})}^{\top }$ are unknown parameter vectors, and ${\mathbf{v}}_i^{\top }=(v_{i1},\cdots ,v_{ip})$ is the explanatory variable vector. Thus, the heteroscedasticity is modeled via $\sigma$.

The EOLLW QR model is defined by Equations (11) and (18), where $\nu$ and $\lambda$ are unknown constants, and it has as special models:

• Rhe exponentiated Weibull (EW) QR model for $\nu =1$;
• the odd log-logistic Weibull (OLLW) QR model for $\lambda =1$;
• and the Weibull QR model for $\nu =\lambda =1$.

Consider a sample $(x_1,{\delta }_1,{\mathbf{v}}_1),\cdots ,(x_n,{\delta }_n,{\mathbf{v}}_n)$ of independent observations, where each random response is defined by $x_i=min\{X_i,C_i\}$, ${\delta }_i=I_{X_i\le C_i}$ (censoring indicator), where $I(·)$ denotes the indicator function. We consider non-informative censoring and the observed lifetimes and censoring times are independent given ${\mathbf{v}}_i$. Let F and C be the sets of individuals for which $x_i$ is the lifetime or censoring time, respectively. Conventional likelihood estimation techniques can be applied here. The log-likelihood function for the vector $\mathbf{\theta }={({\mathbf{\beta }}_1^{\top }\left(\tau \right),{\mathbf{\beta }}_2^{\top }\left(\tau \right),\nu ,\lambda )}^{\top }$ from model (18) has the form

$$\ell \left(\mathbf{\theta }\right)=\sum _{i\in F}{\ell }_i\left(\mathbf{\theta }\right)+\sum _{i\in C}{\ell }_i^{\left(c\right)}\left(\mathbf{\theta }\right),$$

where ${\ell }_i\left(\mathbf{\theta }\right)=log\left[f\left(x_i\right)\right]$, ${\ell }_i^{\left(c\right)}\left(\mathbf{\theta }\right)=log\left[S\left(x_i\right)\right]$, $f\left(x_i\right)$ is the density (11), $S\left(x_i\right)=1-F\left(x_i\right)$ is the survival function, and $F\left(x_i\right)$ is the cdf (10) of $X_i$. The total log-likelihood function for $\mathbf{\theta }$ can be expressed as

$$\begin{array}{ccc}\hfill \ell \left(\mathbf{\theta }\right)& =& r^{†}log\left(\nu \lambda w\right)+\sum _{i\in F}log\left[{\sigma }_i\left(\tau \right)\right]+\sum _{i\in F}[{\sigma }_i\left(\tau \right)-1]log\left(x_i\right)+\nu \sum _{i\in F}log\left(u_i\right)+\hfill \\ & & (\nu \lambda -1)\sum _{i\in F}log(1-u_i)-\sum _{i\in F}{\sigma }_i\left(\tau \right)log\left[{\mu }_i\left(\tau \right)\right]-\hfill \\ & & (\lambda +1)\sum _{i\in F}log\left[{(1-u_i)}^{\nu }+u_i^{\nu }\right]+\sum _{i\in C}log\left[1-\frac{{(1-u_i)}^{\nu \lambda }}{{[{(1-u_i)}^{\nu }+u_i^{\nu }]}^{\lambda }}\right],\hfill \end{array}$$

(19)

where

$$w(\nu ,\lambda ,\tau )=-log\left(\frac{{(1-{\tau }^{\frac{1}{\lambda }})}^{\frac{1}{\nu }}}{{\tau }^{\frac{1}{\nu \lambda }}+{(1-{\tau }^{\frac{1}{\lambda }})}^{\frac{1}{\nu }}}\right),u_i=exp\left\{-w(\nu ,\lambda ,\tau ){\left[\frac{x_i}{{\mu }_i\left(\tau \right)}\right]}^{{\sigma }_i\left(\tau \right)}\right\},$$

and $r^{†}$ is the number of uncensored observations (failures).

The gamlss package in R [42] is used to find the maximum likelihood estimate $\widehat{\mathbf{\theta }}$ of $\mathbf{\theta }$. This package comes from the general class of generalized additive models for location, scale and shape (GAMLSS) (Rigby and Stasinopoulos [43]). These models allow all parameters of a distribution to be modeled as a function of covariates, such as non-parametric, parametric and/or additive smooth functions. Furthermore, they do not have the restriction that the response distribution belongs to a given family such as the exponential family. The package basically has two algorithms: CG (Cole and Green [44]) and RS (Rigby and Stasinopoulos [43]), whose acronyms come from the names of the authors. These algorithms are stable and do not require precise initial values to guarantee convergence. For this reason, we work with the RS algorithm with initial values for ${\mathbf{\beta }}_1\left(\tau \right)$ and ${\mathbf{\beta }}_2\left(\tau \right)$ obtained from the fitted Weibull QR model ($\nu =\lambda =2$). Compared to the CG algorithm, RS is faster for larger datasets and does not use the expected value of cross derivatives, which can be useful when these values are equal to zero. For more details of the algorithms, see [43].

The codes for the reparametric EOLLW distribution in the GAMLSS framework are available at https://github.com/gabrielamrodrigues/EOLLW_quantiles (accessed on 10 February 2023). Following this approach, different regression models can be constructed by incorporating non-parametric smoothing functions, random effects, or other additive terms to the predictors.

Under conditions that are fulfilled for the parameter vector $\mathbf{\theta }$ in the interior of the parameter space but not on the boundary, the asymptotic distribution of $\sqrt{n}(\widehat{\mathbf{\theta }}-\mathbf{\theta })$ is multivariate normal $N_{2p+2}(0,K{\left(\mathbf{\theta }\right)}^{-1})$, where $K\left(\mathbf{\theta }\right)$ is the information matrix. The asymptotic covariance matrix $K{\left(\mathbf{\theta }\right)}^{-1}$ of $\widehat{\mathbf{\theta }}$ can be approximated by the inverse of the $(2p+2)\times$$(2p+2)$ observed information matrix $-\ddot{\mathbf{L}}\left(\mathbf{\theta }\right)$. The approximate multivariate normal distribution $N_{2p+2}(0,-\ddot{\mathbf{L}}{\left(\mathbf{\theta }\right)}^{-1})$ for $\widehat{\mathbf{\theta }}$ can be used in the classical way to construct approximate confidence regions for some parameters in $\mathbf{\theta }$.

We can use the likelihood ratio (LR) statistic for comparing some sub-models with the EOLLW QR model. We consider the partition $\mathbf{\theta }={({\mathbf{\theta }}_1^T,{\mathbf{\theta }}_2^T)}^T$, where ${\mathbf{\theta }}_1$ is the subset of parameters of interest and ${\mathbf{\theta }}_2$ is the subset of remaining parameters. The LR statistic for testing the null hypothesis $H_0:{\mathbf{\theta }}_1={\mathbf{\theta }}_1^{\left(0\right)}$ versus the alternative hypothesis $H_1:{\mathbf{\theta }}_1\ne {\mathbf{\theta }}_1^{\left(0\right)}$ is given by $w^{*}=2\{\ell \left(\widehat{\mathbf{\theta }}\right)-\ell \left(\tilde{\mathbf{\theta }}\right)\}$, where $\tilde{\mathbf{\theta }}$ and $\widehat{\mathbf{\theta }}$ are the estimates under the null and alternative hypotheses, respectively. The statistic w is asymptotically (as $n\to \infty$) distributed as ${\chi }_k^2$, where k is the dimension of the subset of parameters ${\mathbf{\theta }}_1$ of interest.

The standard maximum likelihood techniques can be adopted for the proposed regression, such as the quantile residuals ($q{r}_i$) (Dunn and Smyth [45]), namely

$$q{r}_i={\mathsf{\Phi }}^{-1}\left\{\frac{{(1-{\widehat{u}}_i)}^{\widehat{\nu }\widehat{\lambda }}}{{[{(1-{\widehat{u}}_i)}^{\widehat{\nu }}+{\widehat{u}}_i^{\widehat{\nu }}]}^{\widehat{\lambda }}}\right\},$$

(20)

where

$${\widehat{u}}_i=exp\left\{-\widehat{w}(\widehat{\nu },\widehat{\lambda },\tau ){\left[\frac{x_i}{{\widehat{\mu }}_i\left(\tau \right)}\right]}^{{\widehat{\sigma }}_i\left(\tau \right)}\right\},\widehat{w}(\widehat{\nu },\widehat{\lambda },\widehat{\tau })=-log\left(\frac{{(1-{\widehat{\tau }}^{\frac{1}{\widehat{\lambda }}})}^{\frac{1}{\widehat{\nu }}}}{{\widehat{\tau }}^{\frac{1}{\widehat{\nu }\widehat{\lambda }}}+{(1-{\widehat{\tau }}^{\frac{1}{\widehat{\lambda }}})}^{\frac{1}{\widehat{\nu }}}}\right),$$

$${\widehat{\mu }}_i\left(\tau \right)=exp\left\{{\mathbf{v}}_i^{\top }{\widehat{\mathbf{\beta }}}_1\left(\tau \right)\right\},{\widehat{\sigma }}_i\left(\tau \right)=exp\left\{{\mathbf{v}}_i^{\top }{\widehat{\mathbf{\beta }}}_2\left(\tau \right)\right\},$$

and $\mathsf{\Phi }{(·)}^{-1}$ is the inverse cumulative standard normal distribution.

5. Simulation Study

A simulation study is carried out to verify the accuracy of the MLEs in the EOLLW QR model for the quartiles $\tau =0.25$, $0.50$ and $0.75$, and approximate censoring percentages 0%, 10% and 50%. Just one covariate $v_1\sim$Binomial $(1,0.5)$ is included in the systematic components:

$${\mu }_i=exp({\beta }_{10}+{\beta }_{11}{v}_{1i}),{\sigma }_i=exp({\beta }_{20}+{\beta }_{21}{v}_{1i}),{\nu }_i=exp\left({\beta }_{30}\right),\mathrm{and}{\lambda }_i=exp\left({\beta }_{40}\right),$$

For each combination, $N=1000$ replicas of sizes $n=100$, 300 and 500 are generated. The true values used are: ${\beta }_{10}=1.5$, ${\beta }_{11}=-1.32$, ${\beta }_{20}=0.5$, ${\beta }_{21}=0.2$, ${\beta }_{30}=1.1$ and ${\beta }_{40}=1.4$.

The inverse transformation method is used to generate the lifetimes $x_1,\cdots ,x_n$ from the EOLLW$({\mu }_i,{\sigma }_i,\nu ,\lambda ,\tau )$ distribution, and the censoring times $c_1,\cdots ,c_n$ are determined from a uniform distribution $(0,k)$, where k controls the censoring percentages. For each scenario, the Average Estimates (AEs), Biases and Mean Square Errors (MSEs) of the MLEs are calculated from:

$$\mathrm{AE}\left(\widehat{\mathbf{\theta }}\right)=\frac{1}{N}\sum _{i=1}^N\widehat{{\theta }_i},\mathrm{Bias}\left(\widehat{\mathbf{\theta }}\right)=\frac{1}{N}\sum _{i=1}^N(\widehat{{\theta }_i}-{\theta }_i),\mathrm{MSE}\left(\widehat{\mathbf{\theta }}\right)=\frac{1}{N}\sum _{i=1}^N{(\widehat{{\theta }_i}-{\theta }_i)}^2,$$

(21)

where ${\widehat{\mathbf{\theta }}}^{\top }=({\widehat{\beta }}_{10},{\widehat{\beta }}_{11},{\widehat{\beta }}_{20},{\widehat{\beta }}_{21},{\widehat{\beta }}_{30},{\widehat{\beta }}_{40})$. The software $\mathtt{R}$ is used and Algorithm 1 presents the simulation steps.

Algorithm 1: Simulation study

Table 1. Simulation results from the fitted EOLLW QR model for $\tau =0.25$.

%	$\mathit{\theta }$	True Value		$\mathit{n}=100$			$\mathit{n}=300$			$\mathit{n}=500$
			AEs	Biases	MSEs	AEs	Biases	MSEs	AEs	Biases	MSEs
0%	${\beta }_{10}$	1.50	1.5017	0.0017	0.0006	1.5010	0.0010	0.0002	1.5015	0.0015	0.0001
	${\beta }_{11}$	−1.32	−1.3208	−0.0008	0.0009	−1.3195	0.0005	0.0003	−1.3204	−0.0004	0.0002
	${\beta }_{20}$	0.50	0.4713	−0.0287	0.3744	0.4689	−0.0311	0.0811	0.4661	−0.0339	0.0647
	${\beta }_{21}$	0.20	0.2026	0.0026	0.0232	0.1983	−0.0017	0.0072	0.1996	−0.0004	0.0043
	${\beta }_{30}$	1.10	0.9322	−0.1678	0.6286	1.1410	0.0410	0.1523	1.1606	0.0606	0.1208
	${\beta }_{40}$	1.40	2.1725	0.7725	2.2507	1.3820	−0.0180	0.2594	1.3228	−0.0772	0.1726
10%	${\beta }_{10}$	1.50	1.5019	0.0019	0.0006	1.5007	0.0007	0.0002	1.5012	0.0012	0.0001
	${\beta }_{11}$	−1.32	−1.3213	−0.0013	0.0010	−1.3199	0.0001	0.0003	−1.3205	-0.0005	0.0002
	${\beta }_{20}$	0.50	0.5107	0.0107	0.3516	0.4687	−0.0313	0.0775	0.4664	−0.0336	0.0663
	${\beta }_{21}$	0.20	0.1998	−0.0002	0.0251	0.1992	−0.0008	0.0080	0.1955	−0.0045	0.0046
	${\beta }_{30}$	1.10	0.8796	−0.2204	0.6288	1.1319	0.0319	0.1517	1.1640	0.0640	0.1210
	${\beta }_{40}$	1.40	2.1030	0.7030	1.7695	1.4245	0.0245	0.3137	1.3198	−0.0802	0.1713
50%	${\beta }_{10}$	1.50	1.5083	0.0083	0.0024	1.5015	0.0015	0.0006	1.5005	0.0005	0.0003
	${\beta }_{11}$	−1.32	−1.3262	−0.0062	0.0029	−1.3208	−0.0008	0.0007	−1.3202	−0.0002	0.0004
	${\beta }_{20}$	0.50	0.7797	0.2797	0.3705	0.5016	0.0016	0.1765	0.4869	−0.0131	0.0922
	${\beta }_{21}$	0.20	0.1276	−0.0724	0.0767	0.1830	−0.0170	0.0197	0.1865	−0.0135	0.0127
	${\beta }_{30}$	1.10	0.5845	−0.5155	0.6681	1.0430	−0.0570	0.3344	1.1165	0.0165	0.1563
	${\beta }_{40}$	1.40	2.1757	0.7757	1.6396	1.6131	0.2131	0.6815	1.4305	0.0305	0.2805

,

Table 2. Simulation results from the fitted EOLLW QR model for $\tau =0.50$.

%	$\mathit{\theta }$	True Value		$\mathit{n}=100$			$\mathit{n}=300$			$\mathit{n}=500$
			AEs	Biases	MSEs	AEs	Biases	MSEs	AEs	Biases	MSEs
0%	${\beta }_{10}$	1.50	1.5018	0.0018	0.0005	1.5007	0.0007	0.0002	1.5010	0.0010	0.0001
	${\beta }_{11}$	−1.32	−1.3213	−0.0013	0.0008	−1.3197	0.0003	0.0003	−1.3205	−0.0005	0.0002
	${\beta }_{20}$	0.50	0.5026	0.0026	0.3890	0.4834	−0.0166	0.0697	0.4664	−0.0336	0.0700
	${\beta }_{21}$	0.20	0.2076	0.0076	0.0202	0.2042	0.0042	0.0066	0.2006	0.0006	0.0039
	${\beta }_{30}$	1.10	0.9111	−0.1889	0.6434	1.1457	0.0457	0.1317	1.1694	0.0694	0.1291
	${\beta }_{40}$	1.40	2.1463	0.7463	2.2571	1.3320	−0.0680	0.2332	1.3050	−0.0950	0.1874
10%	${\beta }_{10}$	1.50	1.5016	0.0016	0.0007	1.5004	0.0004	0.0002	1.5005	0.0005	0.0001
	${\beta }_{11}$	−1.32	−1.3203	−0.0003	0.0010	−1.3197	0.0003	0.0004	−1.3200	−0.0000	0.0002
	${\beta }_{20}$	0.50	0.5788	0.0788	0.2982	0.4814	−0.0186	0.0787	0.4676	−0.0324	0.0666
	${\beta }_{21}$	0.20	0.2037	0.0037	0.0236	0.2032	0.0032	0.0082	0.1994	−0.0006	0.0046
	${\beta }_{30}$	1.10	0.8421	−0.2579	0.5646	1.1418	0.0418	0.1473	1.1686	0.0686	0.1223
	${\beta }_{40}$	1.40	1.9985	0.5985	1.5804	1.3647	−0.0353	0.2598	1.3109	−0.0891	0.1954
50%	${\beta }_{10}$	1.50	1.5040	0.0040	0.0033	1.4995	−0.0005	0.0007	1.4996	−0.0004	0.0004
	${\beta }_{11}$	−1.32	−1.3213	−0.0013	0.0036	−1.3185	0.0015	0.0008	−1.3191	0.0009	0.0005
	${\beta }_{20}$	0.50	0.8284	0.3284	0.4059	0.5002	0.0002	0.2775	0.4885	−0.0115	0.1185
	${\beta }_{21}$	0.20	0.1409	−0.0591	0.0730	0.1852	−0.0148	0.0195	0.1931	−0.0069	0.0127
	${\beta }_{30}$	1.10	0.5795	−0.5205	0.6660	1.0593	−0.0407	0.4554	1.1395	0.0395	0.1775
	${\beta }_{40}$	1.40	2.0447	0.6447	1.4838	1.5840	0.1840	0.6743	1.3576	−0.0424	0.2460

and

Table 3. Simulation results from the fitted EOLLW QR model for $\tau =0.75$.

%	$\mathit{\theta }$	True Value		$\mathit{n}=100$			$\mathit{n}=300$			$\mathit{n}=500$
			AEs	Biases	MSEs	AEs	Biases	MSEs	AEs	Biases	MSEs
0%	${\beta }_{10}$	1.50	1.4978	−0.0022	0.0006	1.4974	−0.0026	0.0002	1.4982	−0.0018	0.0001
	${\beta }_{11}$	−1.32	−1.3202	−0.0002	0.0011	−1.3189	0.0011	0.0004	−1.3200	−0.0000	0.0002
	${\beta }_{20}$	0.50	0.4999	−0.0001	0.4778	0.5031	0.0031	0.0630	0.4692	−0.0308	0.0629
	${\beta }_{21}$	0.20	0.2057	0.0057	0.0187	0.1994	−0.0006	0.0064	0.2003	0.0003	0.0038
	${\beta }_{30}$	1.10	0.9645	−0.1355	0.7290	1.1483	0.0483	0.1149	1.1851	0.0851	0.1136
	${\beta }_{40}$	1.40	2.0572	0.6572	2.0673	1.3201	−0.0799	0.2256	1.2860	−0.1140	0.1613
10%	${\beta }_{10}$	1.50	1.4968	−0.0032	0.0008	1.4973	−0.0027	0.0003	1.4981	−0.0019	0.0002
	${\beta }_{11}$	−1.32	−1.3207	−0.0007	0.0011	−1.3192	0.0008	0.0004	−1.3200	−0.0000	0.0003
	${\beta }_{20}$	0.50	0.5707	0.0707	0.3092	0.5012	0.0012	0.0666	0.4817	−0.0183	0.0553
	${\beta }_{21}$	0.20	0.2045	0.0045	0.0233	0.2046	0.0046	0.0069	0.2004	0.0004	0.0043
	${\beta }_{30}$	1.10	0.9004	−0.1996	0.5399	1.1452	0.0452	0.1179	1.1760	0.0760	0.0991
	${\beta }_{40}$	1.40	1.9765	0.5765	1.7309	1.3391	−0.0609	0.2326	1.2787	−0.1213	0.1641
50%	${\beta }_{10}$	1.50	1.4909	−0.0091	0.0037	1.4949	−0.0051	0.0010	1.4953	−0.0047	0.0007
	${\beta }_{11}$	−1.32	−1.3134	0.0066	0.0043	−1.3177	0.0023	0.0011	−1.3172	0.0028	0.0007
	${\beta }_{20}$	0.50	0.8283	0.3283	0.5621	0.5214	0.0214	0.2210	0.5234	0.0234	0.0790
	${\beta }_{21}$	0.20	0.1300	−0.0700	0.0778	0.1920	−0.0080	0.0184	0.1904	−0.0096	0.0120
	${\beta }_{30}$	1.10	0.6175	−0.4825	0.8154	1.0842	−0.0158	0.3476	1.1419	0.0419	0.1129
	${\beta }_{40}$	1.40	2.0680	0.6680	1.6198	1.5014	0.1014	0.6340	1.3051	−0.0949	0.1893

report the findings. For all scenarios, the AEs converge to the true parameter values, and the biases and MSEs decrease when n increases. These facts indicate that the consistency of the estimators hold. In addition, this behavior is verified even for high censoring percentages. We also found the empirical coverage probabilities (CPs) corresponding to the 95% confidence intervals calculated from the simulations.

Table 4. CPs for the 95% nominal level from the fitted EOLLW QR regression model when $\tau =0.25$, $0.50$ and $0.75$ and approximate censoring percentages 0%, 10% and 50%.

$\mathit{\tau }$	$\mathit{\theta }$	$0\%$$\left(\mathit{n}\right)$			$10\%$$\left(\mathit{n}\right)$			$50\%$$\left(\mathit{n}\right)$
		(100)	(300)	(500)	(100)	(300)	(500)	(100)	(300)	(500)
0.25	${\beta }_{10}$	0.939	0.946	0.957	0.948	0.951	0.947	0.922	0.949	0.953
	${\beta }_{11}$	0.954	0.949	0.959	0.946	0.942	0.966	0.937	0.951	0.946
	${\beta }_{20}$	0.973	0.974	0.965	0.972	0.981	0.969	0.962	0.975	0.981
	${\beta }_{21}$	0.948	0.955	0.956	0.946	0.957	0.954	0.937	0.961	0.954
	${\beta }_{30}$	0.980	1.000	1.000	0.978	0.998	1.000	0.954	0.981	0.993
	${\beta }_{40}$	0.990	0.999	0.991	0.993	0.999	0.993	0.996	0.998	1.000
0.50	${\beta }_{10}$	0.950	0.947	0.952	0.940	0.935	0.951	0.907	0.939	0.968
	${\beta }_{11}$	0.950	0.939	0.956	0.947	0.937	0.950	0.913	0.933	0.996
	${\beta }_{20}$	0.969	0.972	0.959	0.955	0.977	0.962	0.959	0.975	0.973
	${\beta }_{21}$	0.959	0.967	0.964	0.953	0.948	0.957	0.945	0.956	0.998
	${\beta }_{30}$	0.986	0.999	1.000	0.970	0.999	1.000	0.941	0.976	0.995
	${\beta }_{40}$	0.986	0.999	0.987	0.995	0.999	0.995	0.995	0.996	0.989
0.75	${\beta }_{10}$	0.950	0.956	0.963	0.941	0.959	0.958	0.895	0.955	0.938
	${\beta }_{11}$	0.946	0.956	0.959	0.954	0.958	0.963	0.914	0.964	0.949
	${\beta }_{20}$	0.952	0.968	0.966	0.952	0.974	0.969	0.939	0.972	0.991
	${\beta }_{21}$	0.968	0.969	0.965	0.962	0.968	0.968	0.936	0.970	0.959
	${\beta }_{30}$	0.980	0.999	1.000	0.979	1.000	1.000	0.952	0.987	0.996
	${\beta }_{40}$	0.980	0.999	0.989	0.991	0.998	0.996	0.986	0.993	0.999

reports CPs values which approach to the nominal level.

6. Application to Gastric Cancer Data

Gastric cancer is the 5th most common cancer worldwide. There are more than one million new cases of this cancer every year, and it ranked as the 2nd leading cause of mortality from cancer in the world. We consider a survival dataset of patients suffering from gastric adenocarcinoma treated by surgery at Helsinki University Hospital in Finland [46] (available at https://doi.org/10.5061/dryad.hb62394, accessed on 29 November 2022 [47]), which contains 301 individuals with approximate censoring of 60%. Here we consider two covariables. The first corresponds to the classification of Lauren (Figure 2a). Various pathological classifications of the disease exist, but that of Lauren is the most common. Originally developed in the 1960s, the classification system adopted cell structural components to separate the patients in three types: well differentiated (non-cardia/intestinal), poorly differentiated (cardia/diffuse), and mixed disease [48]. Based on histology, the two leading types of gastric cancer are diffuse and intestinal [49]. These two types are reflected in the dataset. The second covariable corresponds to the presence of distant metastasis (M1 disease) (Figure 2b). Many patients diagnosed with gastric cancer present distant metastasis, implying a very poor prognosis, generally indicating prophylactic rather than curative treatment ([50,51]). The objective here is to verify the effects of the covariables in different quantiles, so as to obtain a more complete view of this dataset.

Table 5. Descriptive analysis of gastric cancer data.

	0.25	0.50	0.75	Mean
Lauren 1	4.33	8.57	10.18	7.73
Lauren 2	8.32	11.08	13.30	10.39
Pm 0	7.15	9.66	11.90	9.24
Pm 1	10.06	10.06	10.06	8.48

gives a descriptive summary, which includes the mean times, median times and times for the first and third quartiles. We can observe differences for the Lauren classification covariate: between the quantiles, the average time, and the Lauren 1 and Lauren 2 levels. However, we note subtle differences for the presence of distant metastases covariate. Then, the variables considered are ($i=1,\dots ,301$):

• $x_i$ survival time (in years);
• ${\mathrm{cens}}_i$: censoring indicator (0 = censored, 1 = observed);
• $v_{1i}$: Lauren classification (1 = intestinal, 2 = diffuse), defined by a dummy variable (0 = intestinal, 1 = diffuse);
• $v_{2i}$: Presence of distant metastases (pm) (1 = yes, 0 = no)

Regression Model

We compare the EOLLW QR model with the nested OLLW, Exp-W and Weibull models under three systematic components:

${\mathcal{M}}_0=\left\{\begin{array}{c}\mu \left(\tau \right)=exp\left[{\beta }_{10}\left(\tau \right)\right]\hfill \\ {\sigma }_i\left(\tau \right)=exp\left[{\beta }_{20}\left(\tau \right)\right];\hfill \end{array}\right.$

${\mathcal{M}}_1=\left\{\begin{array}{c}\mu \left(\tau \right)=\mu \left(\tau \right)=exp[{\beta }_{10}\left(\tau \right)+{\beta }_{11}\left(\tau \right)v_{1i}+{\beta }_{12}\left(\tau \right)v_{2i}]\hfill \\ {\sigma }_i\left(\tau \right)=exp\left[{\beta }_{20}\left(\tau \right)\right];\hfill \end{array}\right.$

${\mathcal{M}}_2=\left\{\begin{array}{c}\mu \left(\tau \right)=exp[{\beta }_{10}\left(\tau \right)+{\beta }_{11}\left(\tau \right)v_{1i}+{\beta }_{12}\left(\tau \right)v_{2i}]\hfill \\ {\sigma }_i\left(\tau \right)=exp[{\beta }_{20}\left(\tau \right)+{\beta }_{21}\left(\tau \right)v_{1i}+{\beta }_{22}\left(\tau \right)v_{2i}].\hfill \end{array}\right.$

We consider the following quantiles: $\tau =0.10$, $0.25$, $0.50$, $0.75$ and $0.90$.

Table 6. AIC values for some fitted QR models to gastric cancer data.

Model		$\mathit{\tau }$
		0.10	0.25	0.50	0.75	0.90
${\mathcal{M}}_0$	EOLLW	755.3104	760.6965	755.1712	755.1715	755.1759
	OLLW	773.1060	773.1021	773.1010	773.1017	773.1011
	Exp-W	759.9144	759.0376	758.4132	757.2295	757.9724
	Weibull	813.5743	813.5743	813.5744	813.5744	813.5746
${\mathcal{M}}_1$	EOLLW	755.2259	755.1844	755.1742	755.1753	755.1778
	OLLW	774.2566	774.2433	774.2239	774.2288	774.2294
	Exp-W	762.3150	761.8464	761.2808	760.9295	938.5430
	Weibull	811.0754	811.0752	811.0754	811.0755	811.0756
${\mathcal{M}}_2$	EOLLW	750.3085	750.1898	750.1712	750.1812	750.1930
	OLLW	769.1151	769.1514	769.2313	769.2891	769.3191
	Exp-W	755.8550	755.7825	755.7881	755.8164	768.3938
	Weibull	797.4298	797.4352	797.4666	797.5064	797.5421

reports the Akaike information criterion (AIC) values for the fitted QR regression models. The EOLLW QR model under structure ${\mathcal{M}}_2$ gives the lowest values for these quantiles.

Table 7. LR statistics for the ELLOW QR model under structure ${\mathcal{M}}_2$ and some $\tau$ values for the gastric cancer data.

		$\mathit{\tau }$
Models	Hypotheses	0.10	0.25	0.50	0.75	0.90
EOLLW vs. OLLW	$H_0:\lambda =1$ vs. $H_1:H_0\mathrm{is}\mathrm{false}$	20.80(<0.001)	20.95(<0.001)	21.06(<0.001)	21.10(<0.001)	21.12(<0.001)
EOLLW vs. Exp-W	$H_0:\nu =1$ vs. $H_1:H_0\mathrm{is}\mathrm{false}$	7.54(0.006)	7.58(0.005)	7.61(0.005)	7.63(0.005)	20.20(<0.001)
EOLLW vs. Weibull	$H_0:\lambda =\nu =1$ vs. $H_1:H_0\mathrm{is}\mathrm{false}$	51.12(<0.001)	51.24(<0.001)	51.29(<0.001)	51.28(<0.001)	51.34(<0.001)

gives three likelihood ratio (LR) statistics (p-values in parentheses), thus indicating that the EOLLW QR model under structure ${\mathcal{M}}_2$ is better than the others. Thus, we can consider this model as the predictive model.

Figure 3 displays the MLEs and the corresponding confidence intervals along with the interval $[0.01,0.99]$, and

Table 8. Estimation findings from the ELLOW QR model under structure ${\mathcal{M}}_2$ and $\tau$ = $0.10$, $0.25$, $0.50$, $0.75$ and $0.90$ for the current data.

$\mathit{\tau }$	$\mathit{\theta }$	MLEs	SEs	p-Values
0.10	${\beta }_{10}$	0.372	0.050	<0.01
	${\beta }_{11}$	0.691	0.066	<0.01
	${\beta }_{12}$	−1.100	0.206	<0.01
	${\beta }_{20}$	−0.165	0.022	<0.01
	${\beta }_{21}$	0.269	0.036	<0.01
	${\beta }_{22}$	−0.652	0.113	<0.01
	${\beta }_{30}$	2.092	0.016	<0.01
	${\beta }_{40}$	−2.043	0.017	<0.01
0.25	${\beta }_{10}$	1.322	0.050	<0.01
	${\beta }_{11}$	0.468	0.066	<0.01
	${\beta }_{12}$	−0.350	0.210	0.096
	${\beta }_{20}$	−0.226	0.037	<0.01
	${\beta }_{21}$	0.271	0.057	<0.01
	${\beta }_{22}$	−0.670	0.176	<0.01
	${\beta }_{30}$	2.159	0.026	<0.01
	${\beta }_{40}$	−2.054	0.028	<0.01
0.50	${\beta }_{10}$	1.990	0.050	<0.01
	${\beta }_{11}$	0.308	0.066	<0.01
	${\beta }_{12}$	0.211	0.212	0.320
	${\beta }_{20}$	−0.245	0.061	<0.01
	${\beta }_{21}$	0.269	0.078	<0.01
	${\beta }_{22}$	−0.679	0.131	<0.01
	${\beta }_{30}$	2.178	0.034	<0.01
	${\beta }_{40}$	−2.056	0.037	<0.01
0.75	${\beta }_{10}$	2.362	0.050	<0.01
	${\beta }_{11}$	0.220	0.066	<0.01
	${\beta }_{12}$	0.533	0.213	0.013
	${\beta }_{20}$	−0.239	0.077	<0.01
	${\beta }_{21}$	0.266	0.088	<0.01
	${\beta }_{22}$	−0.682	0.105	<0.01
	${\beta }_{30}$	2.172	0.034	<0.01
	${\beta }_{40}$	−2.054	0.036	<0.01
0.90	${\beta }_{10}$	2.565	0.050	<0.010
	${\beta }_{11}$	0.172	0.066	0.010
	${\beta }_{12}$	0.710	0.213	<0.01
	${\beta }_{20}$	−0.231	0.077	<0.01
	${\beta }_{21}$	0.264	0.087	<0.01
	${\beta }_{22}$	−0.684	0.094	<0.01
	${\beta }_{30}$	2.165	0.033	<0.01
	${\beta }_{40}$	−2.053	0.034	<0.01

gives the MLEs and their standard errors (SEs) for the quantiles $\tau =0.10$, $0.25$, $0.50$, $0.75$ and $0.90$ at the significance level of 5%. The following facts can be mentioned:

• The effect of the Lauren classification 2 in comparison with 1 is decreasing along the quantiles and its confidence interval shows significant estimates for all quantiles. These results corroborate with those point quantiles reported in Table 8.
• The effect of the presence of distant metastasis is rising along the quantiles. Its confidence interval includes zero in the interval $[0.25,0.75)$, thus indicating that the covariable is not significant for these quantiles. These results can be noted by the non-significant p-values for $\tau =0.25$ and $0.50$.
• For the parameters ${\beta }_{21}$ and ${\beta }_{22}$, the estimates are significant for both quantiles, thus indicating that those covariables influence the variability of the survival times.
• The estimates corresponding to the shape parameters ${\beta }_{30}$ and ${\beta }_{40}$ are also significant for all quantiles.

Residual Analysis

Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 provide the normal probability plots of the $q{r}_i$’s in Equation (20) under structure ${\mathcal{M}}_2$ for some quantiles. They reveal that the EOLLW QR model is the best among the fitted models. Further, they approximately follow a standard normal distribution, thus indicating adequate fits. Figure 9 shows the index plot of the $q{r}_i$’s for the EOLLW QR model under structure ${\mathcal{M}}_2$. There are few points outside the interval $[-3,3]$ for both quantiles, and a random pattern around zero which show that these models are very adequate to the current data.

7. Concluding Remarks

We introduced a new quantile regression model for censored data based on the reparametrization of the exponentiated log-logistic odd Weibull (EOLLW) distribution in terms of quantiles with two systematic components. We presented some mathematical properties of the reparametrized EOLLW distribution. The proposed quantile regression model is an important extension of other regression models and can be a valuable addition to the survival analysis area. The new regression model also serves as a good alternative for the analysis of lifetime data and may be more flexible than the exponentiated Weibull, odd log-logistic Weibull and Weibull models. Several simulations were performed for different parameter settings, sample sizes and censoring percentages, to assess the accuracy of the maximum likelihood estimators. The usefulness of the new model was also proved by means of a gastric cancer dataset.