Identification of influential observations in high-dimensional survival data through robust penalized Cox regression based on trimming

1.   Introduction

In several lifetime tests, including, industrial, lifetime and clinical applications, progressive censoring is very useful. Progressive censoring permits the removal of the experimental units surviving until the test finishes. Let an experiment of experiment with $n$ independent units in which it is not desirable to detect all failure times under the cost and time limitations, so only part of failures of the units are observed and the other part are removed from the experiment, such a sample is called a censored sample. Assume that one of the units was broken by accident after the test began, but before all of the units had burned out. If the experiment is still ongoing, this unit must be removed from the life test. The progressive censoring scheme gives a methodology for analyzing this type of data in this case. Some of the most important works on this subject are Balakrishnan and Aggarwala ^[1], Balakrishnan ^[2], and Cramer and Iliopoulos ^[3].

The experimentation time can be very long if the units are very reliable, which is a disadvantage of progressive Type-II censored schemes. Kundu and Joarder ^[4] and Childs et al. ^[5] address this problem by proposing a new type of censoring in which the stopping time of the experiment is minimum value of $\{X_{m:m:n}, T\}$, where the time $T$ is fixed time before the start of the test. This type of censored sampling is called a progressive hybrid censoring sample (PHCS). The total time of the experiment under a PHCS will not exceed $T$. Several authors have studied PHCSs. See, for example, Panahi in ^[6], Alshenawy et al. in ^[7], Hemmati and Khorram in ^[8], and Lin and Huang in ^[9].

However, the weakness of a PHCS is that it cannot be implemented when a few failures can be detected before time $T$. For this reason, Cho et al. ^[10] proposed a general type of censoring, called a generalized Type-I PHCS, in which a smaller number of failures is predetermined. A lifetime test experiment would save the time and costs of failures using this censoring scheme. Moreover, the estimates of the statistical efficiency are improved by the experiment having more failures. In the following section, the generalized Type-I PHCS and its advantages are explained. For recent work on this topic, see, for example, Moihe El-Din et al. ^[11], Mohie El-Din et al. ^[12], and Nagy et al. ^[13].

The Weibull distribution is one of the most important in reliability and life testing, and it is widely utilized in various domains such as reliability theory and clinical trials. For this reason, we used this distribution to express truly real data. The Weibull distribution has the probability density (PD), cumulative distribution (CD), survival (S), and hazard (H) functions given as follows.

For Bayesian inference on the Weibull distribution, see, for example, Mohie El-Din and Nagy ^[14], and Lin et al. in ^[15].

In this paper, we address the development of point and interval estimation and classical and Bayesian inference for the Weibull distribution based on the generalized Type-I PHCS. The Bayesian estimate for any parameter $\beta$, denoted by $\widehat{\beta }_{ BS }$, in terms of the squared error loss function (SELF), is the expected value of the posterior distribution and given by

The LINEX loss function (LLF) can be expressed as follows.

The Bayesian estimator of $\beta$, denoted by $\widehat{\beta }_{BL}$ under the (LLF), the value $\widehat{\beta }_{BL}$ that minimizes $E_{\beta |\underline{X}}\left[ L_{BL}\left(\widehat{\beta }, \beta \right) \right]$ is given by

Calabria and Pulcini ^[16] considered the question of the choice of the value of parameter $v$.

The general entropy loss function (GELF) is another widely used asymmetric loss function. It is given by

The Bayesian estimate $\widehat{\beta }_{BE}$ relative to the GE loss function is given by

The remainder of this article is organized as ollows. Section 2 summarizes the model of the generalized Type-I PHCS. Section 3 extracts the maximum likelihood estimates (ML) and the Bayesian estimates for the unknown parameters and SF and HF under three loss functions. Section 4 derives the Bayesian one-sample prediction for all censoring stage failure times of all withdrawn units. In Section 5, we derive the Bayesian prediction for all withdrawn units in the censoring stage $\left\{ R_{i}, i = 1, ..., m\right\}$, which is called one-sample Bayesian prediction; and in Section 6, we derive the Bayesian prediction of an unobserved future progressive sample from the same distribution, which is called two-sample Bayesian prediction. In Section 7, simulation studies are conducted to compare the efficiency of the proposed inference techniques. In Section 8, a real-life data set is used to demonstrate the theoretical findings. Finally, the paper is concluded in Section 9.

2.   The model explanation

Consider lifetime testing in which $n$ equivalent units are tested. The generalized Type-I PHCS is as follows. Let $T > 0$ and $k, m\in \{1, 2, ..., n\}$ be prefixed integers in which $k < m$ with the predetermined censoring scheme $R = \left(R_{1}, R_{2}, ..., R_{m}\right)$ satisfying $n = m+R_{1}+\ldots +R_{m}$. When the first failure occurs, $R_{1}$ of the remaining units are randomly eliminated. When the second failure occurs $R_{2},$ of the surviving units are eliminated from the experiment. This process repeats until the termination time $T^{\ast } = \max \{X_{k:m:n}, \min \{X_{m:m:n}, T\}\}$ is reached, at which moment the reset surviving units are eliminated from the test. The "generalised Type-I PHCS" modifies the PHCS by allowing the experiment to continue beyond $T$ if only a few failures are observed up to $T$. Ideally, the experimenters would like to observe $m$ failures within this system, but they will observe at least $k$ failures. $D$ is the number of failures observed up to $T$ (see Figure 1).

As mentioned earlier, one of observations from the following types is given under the generalized Type-I PHCS:

1. Suppose the $k^{th}$ failure time occurs after $T$. Then, experiment is terminated at $X_{k:m:n}\;$ and the observations are $\{X_{1:m:n} < ... < X_{k:m:n}\}$.

2. Suppose that $T$ is reached after the $k^{th}$ failure and before the $m^{th}$ failure. In this case, the termination time is $T$ and we observe $\{X_{1:m:n} < ... < X_{k:m:n} < X_{k+1:m:n} < ... < X_{D:m:n}\}.$

3. Suppose that the $m^{th}$ fault was discovered after the $k^{th}$ failure and before $T$. Then, the termination time is $X_{m:m:n}$, and we will find $\{X_{1:m:n} < ... < X_{k:m:n} < X_{k+1:m:n} < ... < X_{m:m:n}\}.$

The joint PDF based on the generalized Type-I PHCS for all cases is now given by:

where $R_{j}^{\ast }$ is the $j^{th}$ value of the vector $R^{\ast }$,

$R_{\tau }^{\ast }$ is the number of units eliminated at time $T$, as determined by

and

The likelihood function of $\lambda, \mu$ under the generalized Type-I PHCS can be derived using (1.1) and (1.2) in (2.1), as

where $W\left(\mu |\underline{\bf{x}}\right) = \sum\limits_{i = 1}^{D^{\ast }}\left(R_{i}^{\ast }+1\right) x_{i}^{\mu }+R_{\tau }^{\ast }T^{\mu }$ and $x_{i} = x_{i:D^{\ast }:n}$ for simplicity of notation.

3.   Maximum Likelihood Estimation

From Equation (2.6), the related log-likelihood function can be found as

equating the first derivatives of (3.1) with respect to $\mu$ and $\lambda$ to zero, we obtain

The ML estimators of $lambda$and $mu$ are then obtained by

By using the numerical technique with the Newton-Raphson iteration method, the ML estimates $\widehat{\lambda }_{ ML }$ and $\widehat{\mu }_{ ML }$ can be obtained by solving (3.2) and (3.3), respectively. Due to the invariance property, the related ML estimations of the SF and HF are therefore given by

3.1.   Approximate confidence intervals for $\lambda$ and $\mu$

The observed Fisher information matrix of parameters $lambda$ and $mu$ for large $D^{\ast },$ is given by

where

and a $100(1-\gamma)\%$ two-sided approximate confidence intervals for the parameters $\lambda$ and $\mu$ are then

and

respectively, where $V\left(\widehat{\lambda }\right)$ and $V\left(\widehat{\mu }\right) \;$are the estimated variances of $\widehat{\lambda } _{ML}$ and $\widehat{\mu }_{ML}$, which are given by the first and the second diagonal element of $I^{-1}(\widehat{\lambda }, \widehat{\mu })$ and $z_{\gamma /2}$ is the upper $\left(\gamma /2\right)$ percentile of the standard normal distribution.

3.2.   Approximate confidence intervals for $S(t)$ and $H(t)$

Greene ^[17] used the delta method to construct the approximate confidence intervals for the SF and HF as a function of the MLEs. This method is used in this subsection to determine the variance of the simpler linear function that can be utilized for inference from large samples, as well as the linear approximation of this function. See Greene ^[17] and Agresti ^[18].

where

and

The approximate estimates of $V\left(\widehat{S}\left(t\right) \right)$ and $V\left(\widehat{H}\left(t\right) \right)$ are then supplied, respectively, by

where $G_{i}^{t}$ is the transpose of $G_{i}$, $i = 1, 2.$ These results provide the approximate confidence intervals for $S(t)$ and $H(t)$ are

and

4.   Bayesian estimations

Assuming that both $\lambda$ and $\mu$ are unknown parameters, a natural choice for the prior distributions of $\lambda$ and $\mu$ is to assume that they are independent gamma distributions $G(a_{1}, b_{1})$ and $G(a_{2}, b_{2})$, respectively. As a result, the following is the joint prior distribution.

$a_{1}, \ b_{1}, \ a_{2}, \ b_{2}$ are positive constants. If hyperparameters $a_{1}, \ b_{1}, \ a_{2}, \ b_{2}$ are set as zero, then the informative priors are reduced to the noninformative priors.

Upon combining (2.6) and (4.1), given the generalized Type-I PHCS, the posterior density function of $\lambda, \mu$ is obtained as

where

Thus, from (1.4), the Bayesian estimates of $\lambda$ and $\mu$ under the SELF are as follows.

From (1.6), we obtain the Bayesian estimator of $\lambda$ and $\mu$ under the LLF,

From (1.8), one obtains the Bayesian estimator of $\lambda$ and $\mu$ under the GELF as follows:

Since the integrals in (4.4), (4.5), (4.6), (4.7), (4.8), and (4.9) cannot be computed analytically, the Markov chain Monte Carlo method (MCMC) is used to evaluate these integrals. Depending on the posterior distribution in (4.2), the conditional posterior distributions $\pi _{1}^{\ast }(\lambda |\mu; \underline{\bf{x}})$ and $\pi _{2}^{\ast }(\mu |\lambda; \underline{ \bf{x}})$ of parameters $\lambda$ and $\mu$ can now be computed and written as follows.

and

It is clear that, the posterior density function $\pi _{1}^{\ast }(\lambda |\mu; \underline{\bf{x}})$ is a gamma density, therefore, samples of $\lambda$ can be easily generated. However, the posterior density function $\pi _{2}^{\ast }(\mu |\lambda; \underline{\bf{x}})$ is not a specific distribution; therefore, it is not possible to generate samples directly by standard methods. From theorem 2 of Kundu ^[19], $\pi _{2}^{\ast }(\mu |\lambda; \underline{\bf{x}})$ is a log-concave function; therefore, to generate random samples from these distributions, we use the Metropolis-Hastings ^[20]. The MCMC algorithm can be described as follows.

Assuming $g\left(\lambda, \mu \right)$ is an arbitrary function in $\lambda$ and $\mu$, the Bayesian estimates of $g$ are obtained using the MCMC values as follows.

Based on $SELF$, $LLF$, and $GELF$, the Bayesian estimates of $g$ are then, respectively, given by

The $100(1-\gamma)\%$ Bayesian confidence interval or credible interval $\left(L, U\right)$ for parameter $\beta$ ($\beta$ is $\lambda$ or $\mu$) if

Since the integration in (4.15) cannot be solved analytically, the $100(1-\gamma)$ MCMC-approximated credibility intervals for $\lambda$ and $\mu$ using the $(N-B)$ using the (N - B) generated values after sorting in ascending order, $\left(\lambda ^{(1)}, \lambda ^{(2)}, ..., \lambda ^{(N-B)}\right)$ and $\left(\mu ^{(1)}, \mu ^{(2)}, ..., \mu ^{(N-B)}\right)$, are given as follows,

The absolute difference between the upper and lower bounds determines the length of the credible intervals.

5.   One-Sample Bayesian prediction

For $\rho = 1, 2, ..., R_{j}^{\ast }$, let $Z_{\rho :R_{j}^{\ast }}$ denote the $\rho ^{th}$ order statistic out of $R_{j}^{\ast }$ removed units at stage $j$. Then, the conditional DF of $Z_{\rho :R_{j}^{\ast }}$, given the observed generalized Type-I PHCS, is given, as in Basak et al.^[21], by

where

with $z_{\tau } = T$.

By using (1.1) and (1.2) in (5.1), given a generalized Type-I PHCS, the conditional DF of $Z_{\rho :R_{j}^{\ast }}$ is then given as follows:

where $C_{q} = \frac{(-1)^{q}\binom{\rho -1}{q}R_{j}^{\ast }!}{(\rho -1)!(R_{j}^{\ast }-\rho)!}$ and $\varpi _{q} = q+R_{j}^{\ast }-\rho +1$ for $q = 0, ..., \rho -1.$

Upon combining (4.2) and (5.2) and using the MCMC technique, the Bayesian predictive DF of $Z_{\rho :R_{j}^{\ast }}$, given a generalized Type-I PHCS, is obtained as

The Bayesian predictive SF of $Z_{\rho :R_{j}^{\ast }}$, given generalized Type-I PHCS, is given as

The Bayesian point predictor of $Z_{\rho :R_{j}^{\ast }}$ under the SELF is the mean of the predictive DF, given by

6.   Two-Sample Bayesian prediction

Let $W_{1:\ell :N}\leq W_{2:\ell :N}\leq \ldots \leq W_{\ell :\ell :N}$ be a future independent progressive Type-II censored sample from the same population with censoring scheme $S = (S_{1}, ..., S_{\ell })$. In this section, we develop a general procedure for deriving the point and interval predictions for $W_{s:\ell :N}$, $1\leq s\leq \ell$, based on the observed generalized Type-I PHCS. The marginal DF of $W_{s:\ell :N}$ is given by Balakrishnan et al. ^[22] as

where $1\leq s\leq \rho,$ $c\left(N, s\right) = N\left(N-S_{1}-1\right)...\left(N-S_{1}...-S_{s-1}+1\right), \; M_{q, s} = N-S_{1}-...-S_{s-q-1}-s+q+1, \;$and $c_{q, s-1} = (-1)^{q}\left\{ \left[ \prod\limits_{u = 1}^{q}\sum\limits_{ \upsilon = s-q}^{s-q+u-1}\left(S_{\upsilon }+1\right) \right] \left[ \prod\limits_{u = 1}^{s-q-1}\sum\limits_{\upsilon = u}^{s-q-1}\left(S_{\upsilon }+1\right) \right] \right\} ^{-1}.$

Upon substituting (1.1) and (1.2) in (6.1), the marginal DF of $W_{s:\ell :N}$ is then obtained as

Upon combining (4.2)and (6.2) and using the MCMC method, given a generalized Type-I PHCS, the Bayesian predictive DF of $W_{s:\ell :N}$ is obtained as

From (6.3), we simply obtain the predictive SF function of $W_{s:\ell :N}$, given a generalized Type-I PHCS, as

The Bayesian point predictor of $W_{s:\ell :N}$, $1\leq s\leq m$, under the SELF is the mean of the predictive DF, given by

where $g_{W_{s:\ell :N}}^{\ast }(w_{s}|\underline{\bf{x}})$ is given as in (6.3).

By solving the following two equations, the Bayesian predictive bounds of the $100(1-\gamma)\%$ equi-tailed (ET)interval for $Z_{\rho :R_{j}^{\ast }}$ and $W_{s:\ell :N}$, $1\leq s\leq m$ can be obtained respectively,

where $G^{\ast }(t|\underline{\bf{x}})$ is given as in (5.4) and (6.4), where $L_{ET}$ and $U_{ET}$ denote the lower and upper bounds, respectively. Furthermore, for the highest posterior density (HPD) method, the following two equations need to be solved:

and

where $g^{\ast }(z|\underline{\bf{x}})$ is as in (5.3) and (6.3), where $L_{HPD}$ and $U_{HPD}$ denote the HPD lower and upper bounds, respectively.

7.   Simulation study

In this section, a Monte Carlo simulation study was conducted to compare the efficiency of ML and Bayesian estimates. Using different values of $n, \; m, \; k$ and $T$, $5000$ generalized Type-I PHCSs were generated from the Weibull distribution (with $\lambda = 1$ and $\mu = 2$). The values of $T$ are chosen such that the three cases of generalized Type-I PHCS occur. Thus, in the first case, a $T$ that lies in the first quarter of the data such that $T^{\ast } = X_{k:m:n}$ is chosen. In the second case, a $T$ that lies in the third quarter such that $T^{\ast } = T$ is chosen. Finally, a $T$ that is sufficiently large such that $T^{\ast } = X_{m:m:n}$ is chosen. We computed the ML estimate and the Bayesian estimates of $\lambda$, $\mu$, $S(t)$, and $H(t)$ (with $t = 0.5$) under the SELF, LLF (with $\upsilon$ = 0.5) and GELF (with $\kappa$ = 0.5) using IP and NIP. We also calculated the mean squared error (MSE) and the expected bias (EB) for each estimate.

The $90\%$ and $95\%$ asymptotic and Bayesian credible confidence intervals with the average length (AL) and the estimated coverage probabilities (CPs) for $\widehat{\lambda }$, $\widehat{\mu }$, $\widehat{S(t)}$, and $\widehat{H(t)}$ are computed.

Different samples of size $(n)$ with different effective sample sizes $(m, k)$ are used to conduct the simulation study. The process of removing the SF units is performed with these censoring schemes.

1. Scheme 1: $R_{i} = \frac{2\left(n-m\right) }{m}$ for odd integers $i$ and $R_{i} = 0$ for even integers of $i$.

2. Scheme 2: $R_{i} = \frac{2\left(n-m\right) }{m}$ for even integers $i$ and $R_{i} = 0$ for add integers of $i$..

3. Scheme 3: $R_{i} = 0$ for $i = 1, 2, ..., m-1$, $R_{i} = n-m$ for $i = m$.

All these cases have been assumed according to the case of generalized Type-I progressive censoring and all Bayesian results are computed based on two different choices for the hyperparameters $(a_{1}, b_{1}, a_{2}, b_{2})$.

1. For the case of $IP:\; a_{1} = 200$, $b_{1} = 200$, $a_{2} = 200$ and $b_{2} = 400$ (by putting the marginal prior distribution of $\lambda$ with mean $\frac{a_{1}}{b_{1}} = 1$ and small variance $\frac{a_{1}}{b_{1}^{2}} = 0.005$ and the marginal prior distribution of $\mu$ with mean $\frac{a_{2}}{b_{2}} = 1$ and variance $\frac{a_{2}}{b_{2}^{2}} = 0.005$).

2. For the case of $NIP:a_{1} = b_{1} = a_{2} = b_{2} = 0$.

The simulated results are displayed in the Appendix of this paper.

8.   Numerical example

To illustrate all conclusions reached for the Weibull distribution, we used a real data consists of 19 values. These data refer to breakthrough times of an offending liquid between electrodes at a voltage of 34 kilovolts, as prepared by Viveros and Balakrishnan in ^[23] from Table 6.1 of Nelson (^[24], p.228). We will use these real data to consider the following progressively censored schemes.

Suppose $m = 10,$ $R = (0, 0, 3, 0, 0, 3, 0, 0, 3, 0)$, Then, we have the following progressive data: 0.19, 0.78, 0.96, 2.78, 3.16, 4.15, 4.85, 7.35, 8.01, and 31.75. If we consider a different $T$, then we have three different generalized Type-I PHCSs.

1. Scheme I: Suppose $T = 4$. Since $T < X_{7:10:19} < X_{10:10:19}$, then the experiment would have terminated at $X_{7:7:19},$ with $R^{\ast} = (0, 0, 3, 0, 0, 0, 9)$ and $R_{\tau }^{\ast } = 0$ and we would have the following data: 0.19, 0.78, 0.96, 2.78, 3.16, 4.15, and 4.67.

2. Scheme II: Suppose $T = 7.5$. Since $X_{7:10:19} < T < X_{10:10:19}$, then the experiment would have terminated at $T = 8,$ with $R^{\ast}$ = (0, 0, 3, 0, 0, 3, 0, 0) and $R_{\tau }^{\ast } = 5$ and we would have the following data: 0.19, 0.78, 0.96, 2.78, 3.16, 4.15, 4.85, and 7.35.

3. Scheme III: Suppose $k = 7$ and $T = 35$. Since $X_{7:10:19} < X_{10:10:19} < T$, then the experiment would have terminated at $X_{10:10:19},$ with $R^{\ast } = R$ and $R_{\tau }^{\ast } = 0$ and we would have the following data: 0.19, 0.78, 0.96, 2.78, 3.16, 4.15, 4.85, 7.35, 8.01, and 31.75.

Based on the generated generalized Type-I PHCS and two different choices of hyperparameters $(a_{1}, b_{1}, a_{2}, b_{2})$ as in the Monte Carlo simulation, Table 1 shows the point predictor and $95\%$ Bayesian prediction bounds of $Z_{\rho :R_{k}^{\ast }}$ for three different censoring schemes, and Table 2 shows the point predictor and $95\%$ Bayesian prediction bounds of $W_{s:\ell :N}$ from the future progressively censored sample of size $\ell = 10$ from a sample of size $N = 20$ with progressive censoring scheme $S = (0, 0, 3, 0, 0, 3, 0, 0, 3, 1)$ for the previous four censoring schemes.

9.   Conclusion

The Bayesian and ML estimates of the unknown parameters and the SF and HF of the Weibull distribution when the observed sample is a generalized Type-I PHCS sample are obtained. In the Bayesian approach, the SELF, LLF and GELF based on IP and NIP distributions are considered. The $90\%$ and $95\%$ asymptotic and credible confidence intervals for the parameters and for the SF and HF are also constructed. The Bayesian point and interval predictions of future order statistics samples from the same population for a progressive Type-II of an unpredictable future sample were also developed. From the numerical results, we derive the following conclusions:

1. From Tables 1–2, the HPD prediction intervals appear to be more accurate than the ET prediction intervals, and the means of the Bayesian point predictor inside the Bayesian prediction intervals.

2. From Tables 3–6 in the appendix, the Bayesian estimates using the IP are better than the MLEs. Furthermore, the results of the ML estimates are similar to the Bayesian estimators with NIP. Thus, when we have no prior knowledge of the unknown parameters, it is often easier to use the ML instead of the Bayesian estimators, since the computation of the Bayesian estimator is more complicated. Moreover, in most cases, the MSE decreases as $n$ and $m$ increase.

3. From Tables 7–10 in the appendix, the AL of confidence intervals decreases as $T$ increases, and the credible intervals perform well compared to the asymptotic confidence intervals. Finally, in all cases AL of the confidence intervals, the $95\%$ intervals are larger than the $90\%$ intervals.

Funding

This work is supported by Researchers Supporting Project number (RSP-2021/323), King Saud University, Riyadh, Saudi Arabia.

Acknowledgments

We are grateful to the referees and the editor for their careful reading and their constructive comments, which leads to this greatly improved paper.

The authors acknowledge financial support from the Researchers Supporting Project number (RSP-2021/323), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare there is no conflict of interest.