A novel quantile regression for fractiles based on unit logistic exponential distribution

1.   Introduction

Unit distributions can be obtained using various transformations of well-known continuous distributions, including the negative exponential function transformation. Unit distributions give more flexibility to the original distribution along the unit interval without adding new parameters. Several unit distributions have been used to model data for percentages in many areas such as biology, mortality, recovery rates, economics, health, and reliability analysis, among others. New distributions have been proposed on the unit interval, such as the unit generalized half normal by Korkmaz ^[1], unit inverse Gaussian by Ghitany et al. ^[2], unit Gamma distribution by Consul and Jain ^[3], unit Weibull distribution by Mazucheli et al. ^[4], unit Gompertz distribution by Mazucheli et al. ^[5], unit Omega distribution by Abd El-Monsef et al. ^[6], unit Burr-XII distribution by Korkmaz and Chesneau ^[7], the unit-half normal distribution by Bakouch et al. ^[8], and unit exponential Pareto distribution by Haj Ahmad et al. ^[9]. The unit-exponentiated half-logistic distribution, and the bounded power Lomax were introduced by Hassan et al. (^[10] and ^[11], respectively). Fayomi et al. ^[12] and ^[13] studied the unit-exponentiated Lomax distribution and the unit–power Burr X distribution.

The exponential distribution is widely utilized for modeling real-life data, primarily due to its memoryless feature and analytical simplicity. However, its utility is somewhat constrained because of its constant hazard rate and decreasing density function. To address these limitations and enhance its flexibility, numerous researchers have developed various modifications of the exponential distribution; for example, Hassan et al. ^[14] introduced a new four-parameter extended exponential distribution. Recent extensions include the exponentiated exponential ^[15], beta exponential ^[16], transmuted generalized exponential ^[17], Kumaraswamy transmuted-G family of distributions by ^[18]. A new three-parameter extension of the exponential distribution was introduced by^[19]; also, ^[20] studied a new method for generating distributions with an application to the exponential distribution. Alpha power exponential distribution was discussed by^[21], and extended exponential distribution and its applications were discussed by ^[22]. The Marshall-Olkin alpha power family of distributions was studied by ^[23], the new three-parameter exponential distribution by ^[24], and finally, Topp-leone exponential distribution was studied by ^[25].

Lan and Leemis created a new generalized exponential distribution ^[26]. It was called the logistic-exponential (LE) distribution with explicit-form density and distribution functions. LE is a distribution with two parameters that possess constant, increasing, decreasing, bathtub, and upside-down bathtub failure rate shapes. The only disadvantage of this distribution is that it does not yield explicit-form expressions for the moments, so it must be computed numerically using software like Mathematica, R, or Matlab.

For the LE distribution, the cumulative distribution function $(cdf)$ is given by

with scale parameter $\beta > 0$ and shape parameter $\theta > 0$. The probability density function (pdf) is given by

This paper focuses on a unique special case of the LE distribution, known as the unit-logistic-exponential distribution, confined to the unit interval. The ULE distribution is particularly valuable for its capability to model long-tailed characteristics observed in many real-world datasets. Beyond its utility in capturing long-tailed behavior, the ULE distribution also plays a crucial role in uncovering patterns and trends within the data. By examining the distribution's shape, one can derive insights into the underlying factors influencing the data, which can guide resource allocation decisions or improve the accuracy of future predictions.

Traditional regression methods typically assess differences in outcome variables between populations by focusing on the mean (such as in ordinary least squares regression) or by evaluating a population's average effect (as seen in logistic regression models), after accounting for other explanatory variables. Quantile regression, however, offers the flexibility to examine how the slopes of the regression line vary across different quantiles of the data distribution. For instance, while the median line might remain stable, the 90th quantile prediction line may show a significant upward trend, whereas the 10th quantile prediction line could indicate a notable downward trend.

The primary objective of the research is to explore the new ULE distribution, examine its statistical characteristics, and demonstrate its effectiveness by implementing it in two real-world data examples. The motivations for studying the ULE distribution are as follows: (ⅰ) it is a unique model because it is defined on the unit interval $[0, 1]$ instead of the positive real numbers; (ⅱ) it exhibits significant flexibility in tail behavior, making it useful in risk assessment with relatively better outcomes; and (ⅲ) the strength of this new distribution lies in its ability to model and fit various real-world data with lower error rates compared to other competing models.

This new model demonstrates its effectiveness in fitting the recovery rate for the CD34+ cells, which is a crucial indicator of the sufficiency of peripheral blood stem cells. A more recent purging technique involves the positive selection of CD34+ cells. By focusing on positive CD34+ cell selection, the variability in antigen expression on tumor cells becomes irrelevant, allowing this method to be applied for purging across almost all tumor types, assuming that the tumor cells do not express the CD34+ antigen. The most advanced method currently available for the positive selection of CD34+ cells, enabling their clinical isolation with extremely high purity, is magnetic-activated cell sorting (MACS).

Beyond its applications in modeling the recovery rates, it has a high impact on modeling failure times of components in economic and engineering contexts. Cost-effectiveness analysis (CEA) offers a systematic framework for regulators to compare the quantified benefits of legislative or regulatory decisions against their associated costs. When applied correctly, CEA compels policymakers and regulators to rigorously quantify the health or environmental benefits of potential government actions aimed at reducing risks. It also provides clear metrics for decision-makers, facilitating comparisons among different strategies for addressing the same issue, such as mitigating risks to human health, public safety, or the environment.

The ULE model, while promising and flexible, has some limitations that should be acknowledged as follows: (a) The estimation of parameters in the ULE model, particularly using methods like Bayesian approaches and maximum product spacings, can be computationally intensive. It requires unlimited computational resources and expertise in advanced numerical methods. (b) Although the ULE model was derived from the LE distribution, finding certain statistical properties, such as moments, often requires numerical solutions, such as the quasi-Newton and Markov chain Monte Carlo techniques, that can increase the computational burden and may introduce numerical stability issues in certain scenarios. (c) The ULE model is specifically designed for data within the unit interval $(0, 1)$, limiting its applicability to datasets that can be naturally transformed or confined to this range.

The remainder of this work is arranged as follows: In Section 2, the ULE distribution is presented and its properties are discussed. In Section 3, the maximum likelihood, Bayes, maximum product spacings, least squares estimate of $\theta$ and $\beta$, and the Fisher information matrix are formulated. In Section 4, simulation studies are introduced. In Section 5, a quantile regression model based on ULE distribution is obtained. In Section 6, the two datasets are studied. Finally, the concluding remarks are provided in Section 7.

2.   Unit logistic exponential distribution

Suppose $X$ to be a continuous random variable with LE distribution having two parameters $\theta$ and $\beta\;$, then the cdf of the random variable $Y = e^{-X}$ can be represented by

The pdf for the ULE distribution is defined by

The hazard rate function (HRF) of ULE distribution is as follows

More statistical properties and related functions of the ULE distribution are presented in the following subsections.

2.1.   Shape of density function and hazard rate function

The first derivative of the pdf of ULE distribution is

where $\Psi_1(x) = -\beta \theta +(\beta \theta -1) \left(x^{-\beta }-1\right)^{\theta }+(\beta +1) x^{\beta } \left(\left(x^{-\beta }-1\right)^{\theta }+1\right)-1$.

The above equations require a numerical evaluation to obtain the local maximum, minimum, and inflection points.

Figure 1 illustrates various possible pdf shapes for specific values of the parameters $\theta$ and $\beta$. The pdf can exhibit different forms, such as decreasing (D), increasing (I), unimodal, bathtub-shaped, or increasing-decreasing-increasing, depending on the chosen parameter values. This shape versatility makes the ULE distribution well-suited for analyzing data within the unit interval. Additionally, the left panel of Figure 2 provides a visual representation of the cdf of the ULE distribution, as defined in Eq (2.1).

For the hazard rate function of ULE distribution, we have

and

The first derivative for $h_{ULE}(y)$ is given by

where $\Psi_2(y) = (\beta \theta -1) \left(y^{-\beta }-1\right)^{\theta }+(\beta +1) y^{\beta } \left(\left(y^{-\beta }-1\right)^{\theta }+1\right)-1$. Clearly, the sign of $\frac{\partial\; h_{ULE}(y)}{\partial y}$ depends on the sign of $\Psi(y)$. The HRF (2.3) is increasing for $\beta > 1$ and $\beta \theta > 1$. Additionally, the right panel of Figure 2 shows that for different values of the parameters $\beta$ and $\theta$, the HRF for ULE can accommodate bathtub and IDI shapes.

2.2.   Quantile function and moments

The $\tau_{th}$ quantile of the ULE distribution has a simple closed form and will be utilized later to define a new quantile regression model based on ULE distribution. The quantile regression model for ULE distribution is given by

A random variate can be generated via the inverse transformation method by

To analyze the effect of the parameters $\beta$ and $\theta$ on the skewness and kurtosis of the ULE distribution, we check the trend of the Galton skewness ($sk$) and the Moors kurtosis ($ku$), as described in the following expressions:

and

where $Q(.)$ is the quantile function (2.4). Plots for both measures as functions of $\beta$ and $\theta$ are given in Figure 3. It is observed that the skewness decreases with an increase of $\beta$. For $\beta = 1$, the ULE distribution tends to be symmetric, and for $0 < \beta < 1$, it is positively skewed and decreases toward $0$ as $\theta$ increases. Conversely, for $\beta > 1$, the ULE distribution has negative skewness, which increases toward $0$ as $\theta$ increases. Additionally, the kurtosis increases with $\theta$, while it initially decreases and then increases as $\beta$ increases.

As stated by ^[26], the moment for the logistic exponential distribution exists but cannot be expressed in closed form. Also, the ULE distribution does not have a closed-form expression for its moment and must be computed numerically via

where $\mu_r'$ denotes the $r_{th}$ moment for ULE distribution.

2.3.   Order statistics

When simulating different lifetime systems with some component structures, order statistics, also called ordered random variables, must be considered. Here, along with their typical distributional properties, the order statistics of the ULE distribution are presented. Let $y_{1}, y_{2}, ..., y_{n}$ denote a random sample of size $n$ from the ULE having cdf $F(y)$ and pdf $f(y)$. Then, the order statistics are denoted by, $Y_{(1)}\leq Y_{(2)}\leq...\leq Y_{(n)}$, $Y_{(1)} = min(Y_{1}, Y_{2}, ..., Y_{n}),$ and $Y_{(n)} = max(Y_{1}, Y_{2}, ..., Y_{n})$. The $k^{th}$ order statistic's probability density function is given by

where $y \in (0, 1)$.

In particular, the pdf of $Y_{1}$, can be expressed as

Also, the pdf of $Y_{n}$, can be expressed as

2.4.   Stochastic ordering

Let $Y_1$ and $Y_2$ be two random variables with pdfs $f_1(y)$ and $f_2(y)$, respectively. Let us assume that $Y_1$ is stochastically less than $Y_2$ in terms of the likelihood ratio order ($Y_1\le_{lr}Y_2$) if the ratio $f_2(y)/f_1(y)$ is non-decreasing concerning $y$.

The ULE distribution has a stochastic ordering as illustrated in the next preposition.

Preposition 2.1. Let $Y_1$ and $Y_2$ be two random variables, such that $Y_1 \sim ULE(\theta_1, \beta)$ and $Y_2 \sim ULE(\theta_2, \beta)$. If $\theta_1\ge \theta_2$, then $Y_1\le_{lr}Y_2$.

Proof. The proof is straightforward by taking the first derivative of $f_2(u; \theta_2, \beta)/f_1(y; \theta_1, \beta)$. □

Likelihood ratio order implies a hazard rate order ($Y_1\le_{hr}Y_2$), which in turn implies the usual stochastic order ($Y_1 \le_{st}Y_2$); for further details on stochastic orders, see ^[27].

3.   Methods of estimation

In this research, four estimation methods—namely maximum likelihood estimation, Bayesian estimation, maximum product spacings, and the least squares estimate—are employed to ensure robust and comprehensive parameter estimation. Each method offers unique statistical properties and benefits, contributing to a whole evaluation of the model's performance. The maximum likelihood estimation is a widely used method known for its asymptotic efficiency and consistency, producing parameter estimates that maximize the likelihood of observing the given data. This makes it particularly useful when seeking estimates with desirable large-sample properties. Bayesian estimation, on the other hand, incorporates prior information along with the observed data, providing a flexible framework for parameter inference and enabling the calculation of credible intervals that offer a probabilistic interpretation of parameter uncertainty. The maximum product spacings method is used as an alternative to the maximum likelihood method, particularly effective when dealing with datasets that include extreme or non-uniformly distributed observations. It often provides more stable estimates in cases where traditional likelihood methods may face convergence issues. The least squares estimation method is selected for its simplicity and its ability to minimize the sum of squared differences between observed and predicted values, making it valuable for reducing overall estimation error.

3.1.   Maximum likelihood estimation

The maximum likelihood (ML) estimation method is the most well-known classical inference in statistics. The ML estimates are obtained based on maximizing the log-likelihood function of the ULE distribution. Variance-covariance matrix of the unknown population parameters is utilized to construct the asymptotic confidence intervals for $\theta$ and $\beta$. Numerical methods are used to determine the required estimators, specifically adopting the well-known quasi-Newton method.

Suppose $y_{1}, y_{2}, ..., y_{n}$ represent a random sample of size $n$ with ULE distribution. The likelihood function for ULE distribution is given as

and log-likelihood function of the random sample as

Taking the partial derivatives of Eq (3.2) with respect to $\beta$ and $\theta$, the ML estimates ($\hat{\beta}$ and $\hat{\theta}$) are obtained by solving the following nonlinear systems of equations, as follows:

and

The log-likelihood equations defined in Eqs (3.3) and (3.4) do not provide a closed-form solution. Since obtaining a solution is challenging, it is often more practical to employ nonlinear optimization techniques, such as the quasi-Newton algorithm, to numerically maximize the log-likelihood function. Additionally, explicit confidence intervals for these parameters cannot be directly constructed, necessitating the use of approximate methods to determine the confidence intervals for $\theta$ and $\beta$. To achieve this, we first calculate the second-order partial derivatives, which are essential for deriving the Fisher information matrix.

So, the Fisher information matrix is given by

where

and

The $(1-\zeta)100\%$ confidence interval for the parameters $\beta$ and $\theta$ can be presented as

where $\hat{\theta}$ and $\hat{\beta}$ are the maximum likelihood estimators of $\theta$ and $\beta$, $z_{1-\frac{\zeta}{2}}$ is the percent of the standard normal distribution, and $var(\hat{\theta}), var(\hat{\beta})$ are the asymptotic variances of ML computed utilizing the inverse for the information matrix.

3.2.   Bayes estimation

This section addresses the challenge of deriving Bayes estimators (BE) for the shape and scale parameters of the ULE distribution. Depending on the information known about the parameter, one can choose between using informative or non- informative priors. However, when sample information about the parameters is available, it is preferable to utilize an informative prior. The Gamma prior is preferred in many Bayesian analyses because it is mathematically convenient, conceptually attractive, and computationally friendly. Its flexibility and conjugacy properties make it a powerful tool for representing prior knowledge and updating beliefs in light of new data; see Zellner et al. ^[28]. The following Gamma prior distributions are applied when $\theta$ and $\beta$ are independent parameters.

and

where $a_{i}$ and $b_{i}$, $i = 1, 2$ are supposed to be known hyper-parameters and selected to reflect the prior distribution of the unknown parameters. Thus, we propose to use piecewise independent Gamma priors for both the shape and scale parameters for the ULE distribution because the Gamma distribution is very flexible; see Dey et al.^[29] and Kundu and Howlader ^[30]. The joint prior distribution for $\theta$ and $\beta$ is

The posterior distribution is derived from the likelihood function Eq (3.1) and the prior distribution Eq (3.6). The joint posterior distribution of $\theta$ and $\beta$ is denoted as

The conditional posterior densities of $\theta$ and $\beta$ is presented as

and

The Metropolis-Hastings sampler is necessary for implementing the MCMC technique because the conditional posteriors of $\theta$ and $\beta$ in the previous equations do not conform to any standard distribution. Tierney ^[31] introduced the Metropolis-Hastings (M-H) algorithm within Gibbs sampling to generate posterior samples, and the process is outlined as follows:

$1)$ Suggest initial values to be $(\theta^{(0)}, \beta^{(0)})$.

$2)$ Let $j = 1$.

$3)$ Using the M-H algorithm, create $\theta^{(j)}$ and $\beta^{(j)}$ by using M-H algorithm below using Eqs (3.7) and (3.8) with the normal distribution $N(\theta^{j-1}, var(\theta))$ and $N(\beta^{j-1}, var(\beta))$ and the inverse Fisher information matrix is used to calculate $var(\theta)$ and $var(\beta)$.

(a) We set a $N(\theta^{j-1}, var(\theta))$ and $N(\beta^{j-1}, var(\beta))$ as a proposal distribution for $\theta^{\ast}$ and $\beta^{\ast}$ respectively.

(b) The acceptance rule is defined by the probabilities $\rho_{\theta}$ and $\rho_{\beta}$, which compare the proposed samples against the current samples. We use the acceptance probability:

(c) Set $u_{1}$ and $u_{2}$ from a uniform $(0, 1)$ distribution.

(d) If $u_{1} < \rho_{\theta}$, the request is approved, so make a decision $\theta^{j} = \theta^{\ast}$ or else specify $\theta^{j} = \theta^{j-1}$.

(e) If $u_{2} < \rho_{\beta}$, the request is approved, so make a decision $\beta^{j} = \beta^{\ast}$ or else specify $\beta^{j} = \beta^{j-1}$.

4) Calculating $\theta^{(j)}$ and $\beta^{(j)}$.

5) Set $j = j+1$.

6) Steps from $3$ to $5$ must be repeated $N$ times.

7) To find the credible intervals (CRIs) for $\alpha$ and $\beta$ for $\psi_{k}^{j}$, $j = 1, 2, .., N$, $k = 1, .., 4$ and $(\psi_{1}, \psi_{2}) = (\theta, \beta)$ as $\psi_{k}^{1} < \psi_{k}^{2} < ... < \psi_{k}^{N}$ then the $(1-\gamma)100\%$ CRIs of $\psi_{k}$ is

For practical application, we employ a burn-in period to allow the Markov chain to stabilize. Additionally, thinning reduces autocorrelation in the sampled chains, ensuring that the final posterior estimates represent independent samples. Specifically, after an initial burn-in of $M$ samples, we retain every $k-th$ sample to construct the posterior distribution, as indicated by $\psi_{k}^{(j)}, j = M+1, ..., N$, for sufficiently large $N$.

The variances of the proposal distributions, $var(\theta)$ and $var(\beta)$, are adjusted based on pilot runs to reach an acceptance rate that optimizes the performance of the M-H algorithm. Typically, an acceptance rate between 20% and 30% is targeted to balance convergence speed and sampling efficiency ^[32,33].

In Bayesian analysis, there are four basic elements: the data, the model, the prior, and the loss function. A good Bayesian estimator must minimize the loss function. So, in practical studies, using loss functions in Bayesian inference is important for several reasons, namely: (a) it helps in making decisions that are optimal under uncertainty, and (b) it quantifies the consequences of making incorrect inferences or decisions. By trying to minimize the error itself and the expected cost of errors, which is often more relevant in practice, (c) it can lead to more robust statistical methods. Several authors discussed many loss functions in Bayesian inference, while the squared error loss function was widely used.

The squared error loss (SEL) function is used to get the approximate Bayes estimates of $\hat{\theta}$ and $\hat{\beta}$, being defined as $L\left(\hat{\phi}-\phi \right) = \left(\hat{ \phi}-\phi \right) ^{2}$. The SEL function is symmetric, meaning that it gives equal weights to both over and under-estimation. In real life, we encounter many situations where over-estimation may be more serious than under-estimation, or vice versa. The Bayesian estimation under the SEL is given as:

and

3.3.   Maximum product spacings

Cheng and Amin ^[34] and Ranneby ^[35], introduced the maximum product of spacings (MPS) estimation approach as an alternative ML estimation. They examined the properties of MPS estimators and demonstrated that in certain situations where ML estimators fail to provide consistent and asymptotically efficient estimates, MPS succeeds. These situations include instances where the likelihood function is unbounded from above, heavy-tailed distributions with unspecified scale and location parameters ^[36] and mixture distributions. Consequently, the MPS method overcomes the limitations of the ML method while preserving nearly all of its desirable properties in large sample sizes ^[37]. Assume the parameter vector $\upsilon = (\theta, \beta)$, then the MPS function is presented as

By maximizing $MPS(\upsilon)$ concerning $\upsilon$, one can obtain MPS estimates. The following non-linear equations' simultaneous solutions are also provided for them:

and

where

and

The asymptotic confidence intervals for the parameters $\beta$ and $\theta$, derived from the MPS estimation, rely on the established asymptotic equivalence between the ML and the MPS method, as noted by ^[34,38,39]. Thus, the $100(1-\zeta)\%$ confidence interval for the parameters using MPS is formulated as follows:

where $\hat{\theta}_{MPS}$ and $\hat{\beta}_{MPS}$ are the MPS estimates of $\theta$ and $\beta$, $z_{1-\frac{\zeta}{2}}$ is the percent of the standard normal distribution, and $var(\hat{\theta}_{MPS}), var(\hat{\beta}_{MPS})$ are the asymptotic variances computed utilizing the inverse of the observed information matrix (3.5).

3.4.   Least squares estimation

Least-squares methods produce the estimated parameters with the highest probability (maximum likelihood) of being true if several critical conditions are guaranteed. Least-squares parameter estimation is a basic procedure for evaluating the confidence intervals for the unknown model's parameters. It discusses the practical ways of applying least-squares techniques to experimental data.

The least square estimates (LSE) for $\hat{\theta}_{LSE}$ and $\hat{\beta}_{LSE}$ of $\theta$ and $\beta$, respectively, are observed by minimizing the function:

with respect to $\upsilon$, where $E[F(y_{i}, \theta, \beta)] = \frac{i}{(n+1)}$ for $i = 1, 2, ..., n.$. Then, $\hat{\theta}_{LSE}$ and $\hat{\beta}_{LSE}$ are solutions of the following equations:

and

where $F(y_{i}, \theta, \beta)^{'; \theta}$ and $F(y_{i}, \theta, \beta)^{'; \beta}$ as mentioned in Eqs (3.9) and (3.10).

The bootstrap percentile method (PCI) is commonly used to construct confidence intervals. The following procedure describes the steps to construct PCI:

ⅰ) Calculate (LSE) $\hat{\beta}{LSE}$ and $\hat{\theta}{LSE}$ for the ULE distribution.

ⅱ) Generate a bootstrap sample by using the estimates $(\hat{\beta}{LSE}, \hat{\theta}{LSE})$, and then obtain the bootstrap estimate, denoted as $(\hat{\beta}b, \hat{\theta}b)$, based on the bootstrap sample.

ⅲ) Repeat step $(ii)$ $B$ times to obtain the estimates $(\hat{\beta}{b_1}, \hat{\theta}{b_1})$, $(\hat{\beta}{b_2}, \hat{\theta}{b_2})$, $\dots$ $(\hat{\beta}{b_B}, \hat{\theta}{b_B})$.

ⅳ) Construct the $100(1-\zeta)\%$ bootstrap confidence intervals for the parameters $\beta$ and $\theta$ by using the $\frac{\zeta}{2}$ and $1-\frac{\zeta}{2}$ quantiles of the empirical distribution of the bootstrap estimates for $\beta$ and $\theta$, respectively.

4.   Simulation

The simulation study provides a comprehensive assessment of the estimation methods for estimating the parameters $\theta$ and $\beta$. By generating synthetic datasets with known parameter values, one can evaluate how well these estimation methods converge to the true values as the dataset size increases. Measures such as mean squared errors (MSEs) and average biases (ABs) are used to assess the convergence and accuracy of the point estimators. The performance of $95\%$ confidence intervals is measured based on average widths (AW) and coverage probabilities (CP).

To generate observations, we used the ULE distribution and applied the inverse transformation method with various values of the parameters $(\theta, \beta) = (1, 3),$ $(0.45, 0.5),$ $(3, 0.4),$ $(0.25, 1),$ and $(1.9, 0.8)$. These cases represent various trends of the pdf (2.2) for the ULE distribution. The simulation is redone 1000 times, with varying sample sizes of $n = ,$ 100, $and$ 150 $. For each simulation iteration, we generate 1000 bootstrap samples and accumulate 10000 observations for the Markov chain Monte Carlo (MCMC) algorithm. We used the first 1000 observations as burn-in to minimize initial distribution effects. Additionally, in the thinning process, every fifth observation is selected to reduce dependence and enhance the robustness of the analysis. The simulation results are presented in Tables 1–5. Mathematica scripts were used for all simulation results.

The findings from these tables regarding accuracy and biased-ness are summarized as follows:

● The MSE for all estimation methods decreases as the sample size increases from 50 to 150, indicating improved accuracy with larger sample sizes.

● The ABs generally reduce with larger sample sizes, suggesting that estimators become less biased.

For comparison of estimators from Tables 1–5, the following points are indicated:

● ML estimates show relatively low MSEs and consistent CPs, making it a reliable choice across sample sizes.

● MPSE has competitive MSEs and ABs, but its CPs for $\theta$ can be slightly lower, especially at smaller $n$.

● Bayesian estimation provides good MSEs and high CPs, often outperforming other methods in maintaining high coverage.

● The average width of the confidence intervals decreases as the sample size increases for all estimation techniques indicating more precise estimates with larger datasets. Also, it is observed that the approximate confidence interval (ACI) for MPS has smaller interval widths than the other estimation approaches.

● The coverage probabilities are close to the nominal level, showing that the confidence intervals produced by these methods effectively contain the true parameter values. It is observed that LSE is showing coverage probabilities higher than the nominal level.

Overall, Bayesian estimation offers strong performance with high CPs and reasonable AWs. ML estimation consistently performs well, with solid MSEs, ABs, and CPs, making it a dependable choice. While small-sample scenarios are not the primary focus of this study, we acknowledge their practical importance. Bootstrap or Bayesian approaches will be more appropriate for small-sample analyses rather than the asymptotic confidence intervals.

5.   Quantile regression model for ULE distribution

This section proposes a new quantile regression model based on ULE distribution for component response variables, offering a preference for commonly used beta, Kumaraswamy, and unit-Weibull regression models. When incorporating covariate information into regression analysis of a probability distribution, it is common to connect the mean response to the covariates. In contrast to modeling the mean of the response, quantile regression models were pioneered by ^[40] to model the conditional quantiles of the response variable as a function of the covariates without any distributional assumptions on the error term, and these models are also considered robust regressions.

Since the mean of the ULE distribution does not have a tractable closed-form expression, quantile regression models are a solid substitute. The quantile function of the ULE distribution has a simple treatable form, allowing its cdf and pdf to be expressed using a re-parametrization established on the quantile function in Eq (2.4). Let $\mu = Q(\tau, \theta, \beta)$ and the parameter $\beta = -\frac{\log \left[\left(\frac{1}{\tau}-1\right)^{1/\theta }+1\right]}{\log(\mu)}$. Then, the cdf and pdf of the re-parametrized ULE distribution are given by

and

where $0 < y < 1, 0 < \mu < 1$ and $0 < \tau < 1$.

Given independent random variables $X_i; i = 1, \dots, n$ from the re-parameterized ULE distribution with pdf in Eq (5.2), it is feasible to establish the ULE quantile regression. This regression assumes the functional relation of the median of $Y_i$

where $\boldsymbol{x_i} = (1, x_{1i}, x_{2i}, \dots, x_{ki})$ represents the covariates vector; also, $\boldsymbol{\delta} = (\delta_{0}, \delta_{1}, \dots, \delta_{k})^{T}$ is the coefficients vector in the regression model. Various selections for the link function $g(.)$ can be treated; however, in this context, we only consider the logit link function:

Thus we have

When $\tau = 0.5$, it becomes obvious that the conditional median response is obtained.

5.1.   Maximum likelihood estimation

Referring to Eqs (5.2) and (5.3), the log-likelihood function of the ULE quantile regression model is given by

where ${\bf{\Omega}} = (\theta, \boldsymbol{\delta}^T)^T$ is the vector of the unknown parameters. The ML $\hat{{\bf{\Omega}}} = (\hat{\theta}, \hat{\boldsymbol{\delta}}^T)^T$ can be obtained by directly maximizing the log-likelihood function $\ell({\bf{\Omega}})$ using mathematical packages such as Mathematica or R.

In the classical approach, the confidence intervals for the parameters were constructed based on asymptotic normality, i.e., $\hat{{\bf{\Omega}}}$ can be approximated by a $(k+1)-$variate normal distribution with zero means and covariance matrix $I^{-1}(\hat{{\bf{\Omega}}})$, where $I({\bf{\Omega}})$ is the observed information matrix defined by

The asymptotic $(1-\alpha)100\%$ confidence interval for the vector parameter $\Omega_{j}; j = 1, \dots, k+1$ is $\hat{\Omega}_{j}\pm{{z}_{\frac{\alpha }{2}}}se(\hat{\Omega}_{j})$, where ${{z}_{\frac{\alpha }{2}}}$ represents the upper ${\frac{\alpha }{2}}$ percentile with a standard normal distribution and $se(\hat{\Omega}_{j})$ is the square root of the $j^{th}$ diagonal entry of the Fisher matrix $I^{-1}({\hat{\bf{\Omega}}})$.

5.2.   Model sufficiency

Residual analysis plays a vital role in testing the goodness of fit for a regression model. For this purpose, two types of residuals are employed: randomized quantile residual (RQR) ^[41] and Cox-Snell residual (CSR) ^[42]. The RQRs are computed as follows:

where $G(y_i, \mu_i, \theta)$ denotes the cdf of the ULE quantile regression model in Eq (5.1), and $\Phi^{-1}[.]$ represents the inverse cdf of the standard normal distribution. The adequacy of the fitted model is indicated by the RQRs ($\hat{r}_i$) following a standard normal distribution.

On the other side, the Cox-Snell residuals are written as:

where $\bar{G}(y_i, \mu_i, \theta) = 1 - G(x_i, \mu_i, \theta)$ represents the survival function of the ULE quantile model of regression. The Cox-Snell residuals ($\hat{e}_i$) are expected to follow an exponential distribution with a scale parameter of $1$.

In summary, both the randomized quantile residual and Cox-Snell residual are valuable tools to evaluate the goodness of fit for the ULE quantile regression model and can provide insights into the appropriateness of the model for a given dataset.

5.3.   Simulation

This subsection outlines a simulation study to assess the effectiveness of maximum likelihood estimators (MLEs) and their corresponding asymptotic confidence intervals (CIs) for the parameters in the ULE quantile regression model. The evaluation metrics include mean squared error (MSE) and average bias to assess the point estimates of the parameters, while coverage probability (CPs) along with average width (AWs) to assess the efficiency of the constructed $95\%$ CIs. Two distinct scenarios are considered within the simulation design:

(a) The case of one covariate: the simulation study is performed as

In this case, we consider $(\delta_0, \delta_1)$ $= (-2, 2)$.

(b) The case of two covariates: the simulation study is performed as

In this case, we consider $(\delta_0, \delta_1, \delta_2)$ $= (2, -2, 2)$.

In both cases, the true value for the parameter $\theta$ is $0.5$ and the quantile values are $\tau = 0.1,$ $0.25,$ $0.5,$ $0.75,$ and $0.9$. The covariate $x_{i1}$ is generated from Bernoulli distribution with parameter $0.5$ while the covariate $x_{i2}$ is generated from the standard normal distribution for $n =$$50,$$100$ and $150$. For each scenario, the simulation is repeated $1000$ times.

The simulation results are given in Tables 6 and 7. The results in these tables confirm desirable properties for the maximum likelihood estimates for the parameters of the ULE quantile regression. From Tables 6 and 7, we observe that the MSE, ABs, and AWs decrease with increasing sample size. These metrics are consistently attaining their lowest values for the parameter $\theta$. The MSEs and ABs attain their highest values when $\tau = 0.9$, particularly for small sample sizes. The coverage probabilities remain close to the nominal $95\%$ confidence level.

6.   Applications

This section explores the advantages of the practical use of ULE distribution and its regression model compared to other models commonly used for modeling data within unit intervals. The unit models used in this comparison are the Beta, Kumaraswamy, Topp-Leon, unit-Burr Ⅻ, and unit half-logistic geometric (UHLG) distributions.

To assess and compare the fitted distributions, we employ various goodness-of-fit measures. These include the Kolmogorov-Smirnov distance (KS) with corresponding p-value, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Anderson-Darling statistic ($A^\star$), and Cramer-von Mises statistic ($W^\star$).

Additionally, two graphical methods are used to assess the goodness of fit for the regression models:

$i)$ The quantile-quantile (QQ) plot, which compares randomized quantile residuals against normal quantiles with a simulated enclosure, as proposed by Atkinson ^[43].

$ii)$ The probability-probability (PP) plot compares the empirical probabilities of Cox-Snell residuals with those from the standard exponential distribution.

The model is considered the best fit for the data if its values are closest to the diagonal bar in both plots. Mathematica scripts were used for all numerical calculations and to generate the figures.

6.1.   Recovery rate applications

The recovery rate of CD34+ cells is a crucial factor in assessing the sufficiency of peripheral blood stem cell (PBSC) collection for bone marrow transplantation. This procedure facilitates rapid hematologic recovery following myeloablative therapy for various malignant hematological diseases.

This study analyzed data from 239 patients who underwent autologous PBSC transplantation, focusing on CD34+ cell recovery rates as documented in ^[44]. Additionally, we incorporated relevant covariates that could potentially influence CD34+ cell viability recovery:

● Gender ($y_1$): 0 denotes female and 1 denotes male;

● Chemotherapy ($y_2$): 0 for a single-day chemotherapy regimen and 1 for a three-day regimen;

● Age ($y_3$): Adjusted patient's age, calculated as the current age minus 40.

A univariate modeling approach was employed to assess the suitability of the ULE distribution for modeling CD34+ cell recovery rates. The ML estimates of the ULE distribution's parameters and various goodness-of-fit statistics were calculated and are presented in Table 8. The results from this table strongly suggest that the ULE distribution provides the best fit for the given dataset and it generally outperforms other models across multiple criteria (AIC, BIC, -Loglik, $W^\star$, $A^\star$), indicating that it best fits the CD34+ cell recovery data. The Unit-Omega and Unit-BurrXII models perform well but do not surpass the ULE model according to AIC, BIC, or -Loglik. To further visualize this relationship, Figure 4 was generated. This figure includes histogram plots of the data, overlaid with the estimated densities of the ULE distribution, as well as the Kaplan-Meier curve and the estimated survival curves.

To derive Bayesian estimates and credible intervals, we run the MCMC algorithm 55,000 times, discarding the first 5,000 iterations as burn-in and selecting every 5th value for thinning. The burn-in period is determined based on pilot runs of the Metropolis-Hastings algorithm. Specifically, we monitor the convergence behavior of the chains for the parameters $\theta$ and $\beta$, examining trace plots to identify when the Markov chain stabilizes around the posterior distribution Figure 5 shows a smooth histogram of the marginal posterior density, a trace plot where each chain remains within a similar region without distinct trends, and results from the first 100 lags of the autocorrelation function (ACF). These observations reflect the efficient convergence of the MCMC chain.

Percentile bootstrap confidence intervals were computed from 1,000 bootstrap samples of the parameters. Table 9 displays the maximum likelihood estimates (MLE), bootstrap standard errors (BSE), and the associated $95\%$ confidence intervals. It is evident that the MLE, BSE, Mean Posterior Standard Errors (MPSE), and LSE for the parameter $\theta$ are closely aligned, with confidence intervals that are nearly of the same length.

A multivariate regression model is now employed to examine the effects of gender, chemotherapy, and age on the recovery rate. The proposed regression model is:

Table 10 presents the fitting results for beta, Kumaraswamy, and ULE regression models. Table 10 shows MLE for $\alpha$, and $\delta_i, i = 0, 1, 2, 3$, along with their respective standard errors (SE) and p-values. Additionally, each regression model is accompanied by AIC and BIC statistics.

According to Table 10, there is a statistically significant positive effect of the type of chemotherapy and age on the recovery rate of CD34+ cells, whereas the gender parameter does not significantly affect the recovery rate in any of the regression models.

Additionally, among the beta and Kumaraswamy regression models, the ULE regression model offers the best fit to the data, as indicated by its lower AIC and BIC values. This finding is further explored in Figure 6, which displays the P-P plot of Cox-Snell residuals and Q-Q plots of the randomized quantile residuals with a simulated envelope for all the fitted regression models.

6.2.   Data on the cost-effectiveness of risk management

The dataset in question is presented by Schmit and Roth ^[45], which includes 73 responses from a survey distributed to 374 risk managers of major North American corporations. Schmit and Roth ^[45] sought to assess the cost-effectiveness of different management strategies in reducing a company's risk of property losses and accidents, taking into account specific company characteristics like size and industry. The response variable, y (Firmcost), represents the firm-specific ratio of premiums plus uninsured losses relative to total assets. The covariates associated with this response variable are as follows:

● Suppose ($y_1$): the firm-specific ratio of the total per-occurrence retention levels, as assessed by the corporate risk manager.

● Cap ($y_2$): 1 if the firm utilizes a captive insurance, and zero if it does not.

● Sizelog ($y_3$): the logarithm of the firm's total asset value.

● Indcost ($y_4$): the industry average of premiums and uninsured losses relative to total assets, as reported in the 1985 cost of risk survey (a risk measurement).

● Central ($y_5$): the role of the local manager in determining local retention levels, as evaluated by the corporate risk manager.

● Sizelog ($y_6$): the significance of analytical tools in risk management decision-making, as assessed by the corporate risk manager.

Initially, univariate modeling is applied to the response variable, risk management cost-effectiveness data, to evaluate the performance of the ULE distribution. Table 11 presents the ML estimates for the parameters along with the goodness-of-fit statistics. The findings in Table 11 indicate that the ULE model provides the best overall fit to the risk management cost-effectiveness data, as evidenced by its favorable AIC, BIC, -Loglik, and goodness-of-fit test results (K-S, Cramér-von Mises, and Anderson-Darling). The unit-Omega model performs reasonably well but does not surpass the ULE model across the proposed criteria. This conclusion is further illustrated in Figure 7, which displays histogram plots of the data with the estimated densities, as well as the Kaplan-Meier curve with the estimated survival curves.

To compute the Bootstrap errors (BE) and confidence intervals (CI), we ran the MCMC algorithm 55,000 times, discarding the initial 5,000 iterations as burn-in and selecting every 5th value for thinning. As mentioned in the previous section, the pilot runs of the Metropolis-Hastings algorithm determine the burn-in period through monitoring convergence trace, as shown in Figure 8. This figure illustrates a smooth and well-behaved histogram of the marginal posterior density, a trace plot showing that each chain consistently explores the same region without evident trends, and the results from the first 100 lags of the autocorrelation function (ACF). These findings illustrate fast convergence performance for the MCMC chain. Percentile bootstrap confidence intervals were obtained from 1,000 bootstrap samples of the parameters. Table 12 presents the ML estimates, bootstrap standard errors (BSE), and the corresponding $95\%$ confidence intervals. It is noted that the MLE, BSE, mean posterior standard errors (MPSE), and LSE of parameter $\theta$ are similar, with their confidence intervals nearly the same length. The regression model for $\mu_i$ is presented as

Table 13 presents the fitting results for beta, Kumaraswamy, and ULE regression models. It shows ML estimates for $\alpha$, and $\delta_i, i= 0, 1, 2, 3, 4, 5, 6$, along with their respective standard errors (SE) and p-values. Additionally, each regression model is accompanied by AIC and BIC statistics.

Based on Table 13, we can notice that there is a statistically significant positive impact of the type of chemotherapy and age on the risk management cost-effectiveness data, while the parameter gender has no statistically significant impact on the response variable recovery rate at the usual level for all regression models.

Furthermore, compared to the beta and Kumaraswamy regression models, the ULE regression model provides the most optimal fit to the data as it has the least AIC and BIC. Figure 9 provides a deeper analysis of this result, displaying the P-P plot of Cox-Snell residuals and the Q-Q plots of randomized quantile residuals, complete with a simulated envelope, for all the regression models fitted.

7.   Conclusions

This paper introduced a versatile two-parameter unit distribution called the ULE distribution, designed for modeling datasets confined to the interval $(0, 1)$. We derived and analyzed the survival function and hazard rate function of the ULE distribution, presenting their behavior through graphical representations. We investigated various methods for parameter estimation of the ULE distribution, including maximum likelihood estimation, maximum product spacing, least squares, and Bayesian estimation. An extensive simulation study reveals that Bayesian estimation demonstrates robust performance, characterized by high CPs and acceptable AWs. ML estimate consistently shows strong results with reliable MSEs, ABs, and CPs, making it a trustworthy option. LSE and MPSE may exhibit greater variability with smaller sample sizes, but their performance improves and remains satisfactory as $n$ increases. Additionally, we proposed a novel quantile regression model based on the ULE distribution for unit response variables, providing an alternative to the commonly used beta, Kumaraswamy, and Unit-Weibull regression models. It was found that the randomized quantile residual and Cox-Snell residual are useful tools for assessing the goodness of fit of the ULE quantile regression model, offering valuable insights into the model's suitability for a specific dataset. Our empirical analysis focused on datasets related to the recovery rate of CD34+ cells and the cost-effectiveness of risk management. The numerical analysis with the help of goodness of fit tests emphasizes the suitability of the ULE distribution to model these data. One of the limitations of this model is the time and cost constraints, as we consider a complete data set, units may take a long time to fail, or it can be expensive to perform a complete sample experiment, especially in medical and biological fields. Censoring schemes can be employed to adjust this limitation. Hence, our plan for future work is to implement censored samples instead of complete ones. Another limitation is that even though the ULE distribution is proposed as effective for fitting real-world data, its performance might vary depending on the nature and characteristics of the datasets. It may not always outperform existing models in all scenarios, especially if the data do not exhibit the specific long-tailed behavior for which the ULE distribution is suited. The ULE distribution is defined on the unit interval $(0, 1)$. This limitation implies that it may not be suitable for datasets that extend beyond this interval or where transformations to fit this range introduce bias or distortions.

Author contributions

Hanan Haj Ahmad: Methodology, validation, investigation, resources, writing–review & editing, funding acquisition; Kariema A. Elnagar: Methodology, conceptualization, investigation, resources, data curation, writing–original draft, Writing–review & editing. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [GRANT No. KFU242542].

Conflicts of interest

The authors declare no conflicts of interest.