A Partitioning-Based Approach to Variable Selection in WLW Model for Multivariate Survival Data

Abstract

In this paper, we propose a new variable selection method using a partitioning-based estimating equation for multivariate survival data to simultaneously perform variable selection and parameter estimation. The main idea of the partitioning-based estimating equation is to partition the score function into small blocks. We construct our method using the SCAD penalty function and achieve the purpose of directly selecting variables through the estimating equation. We further establish asymptotic normality and prove that our method achieves the oracle property. Moreover, we use a simple approximation of the penalty function such that our method can be implemented algorithmically. We conducted simulation studies to validate the performance of our method and analyzed the dataset from the Colon Cancer Study.

1. Introduction

Biomedical research often involves multivariate survival data, such as cancer patients facing local recurrence and repeated hospitalizations in the context of chronic disease management. A significant characteristic of these data is that the survival times of the same individual are correlated. Theoretical development is relatively slow due to the complexity of the dependency structure among survival times, and researchers have mainly focused on modeling the marginal distributions of survival times (see Liang [1], Lin [2], and Spiekerman [3]), where the dependence is not specified.

The most widely used approach is the WLW method [4]. Suppose that there is a random sample of n subjects. Let $T_{ki}$ be the kth survival time of the ith subject, where $k=1,\dots ,K$ and $i=1,\dots ,n$. In the WLW method, the marginal hazard function of $T_{ki}$ is assumed to take the form

${\lambda }_{ki}\left(t\right|Z_{ki})={\lambda }_{k0}\left(t\right)exp\left({\mathit{\beta }}_{k0}^T{Z}_{ki}\left(t\right)\right),t⩾0,$

where $Z_{ki}$ is a p-dimensional possibly covariate vector, ${\lambda }_{k0}\left(t\right)$ is the unspecified baseline hazard function, and ${\mathit{\beta }}_{k0}$ represents the true values of unknown regression coefficients. Let ${\mathit{\beta }}_0={({\mathit{\beta }}_{10}^T,\dots ,{\mathit{\beta }}_{K0}^T)}^T$ and $C_{ki}$ be the censoring time. Let ${\tilde{T}}_{ki}=T_{ki}\wedge C_{ki}$ and ${\delta }_{ki}=I(T_{ki}⩽C_{ki})$. Assume that $T_{ki}$ and $C_{ki}$ are given independent covariates of $Z_{ki}$. Then, the marginal partial likelihood ([5,6]) is

$L_k\left({\beta }_k\right)=\prod _{i=1}^n{\left({\displaystyle \frac{exp\left({\beta }_k^T{Z}_{ki}\left({\tilde{T}}_{ki}\right)\right)}{{\Sigma }_{l\in R_k\left({\tilde{T}}_{ki}\right)}exp\left({\beta }_k^T{Z}_{kl}\left({\tilde{T}}_{ki}\right)\right)}}\right)}^{{\delta }_{ki}},$

and the corresponding score function is

(1)${\displaystyle \frac{\partial log{L}_k\left({\beta }_k\right)}{\partial {\beta }_k}}=\sum _{i=1}^n{\int }_0^{\infty }\left(Z_{ki}\left(t\right)-{\overline{Z}}_{ki}({\beta }_k,t)\right)d{N}_{ki}\left(t\right),$

where $R_k\left(t\right)=\{i:{\tilde{T}}_{ki}⩾t\}$ is the risk set at t, $Y_{ki}\left(t\right)=I({\tilde{T}}_{ki}⩾t)$, $N_{ki}\left(t\right)={\delta }_{ki}I({\tilde{T}}_{ki}⩽t)$, and

${\overline{Z}}_k({\beta }_k,t)={\displaystyle \frac{\sum Z_{ki}\left(t\right)Y_{ki}\left(t\right)exp\left({\beta }_k^T{Z}_{ki}\left(t\right)\right)}{\sum Y_{ki}\left(t\right)exp\left({\beta }_k^T{Z}_{ki}\left(t\right)\right)}}.$

Thus, we have K sets of estimating equations as follows:

(2)$\sum _{i=1}^n{\int }_0^{\infty }(Z_{ki}\left(t\right)-{\overline{Z}}_k({\beta }_k,t))d{N}_{ki}\left(t\right)=0,k=1,2,\dots ,K.$

We can obtain the estimators ${\widehat{\mathit{\beta }}}_1,\dots ,{\widehat{\mathit{\beta }}}_K$. Let ${\widehat{\mathit{\beta }}}_W={({\widehat{\mathit{\beta }}}_1^T,\dots ,{\widehat{\mathit{\beta }}}_K^T)}^T$, which is known as the WLW estimator.

Cui [7] proposed a new method to improve the efficiency of the WLW method, which is called the partition method or partition-estimating equation in this paper. The main idea of their method is to partition the score function into small blocks. He further explored the situation where the number of blocks increases with the sample size, and he established the asymptotic normality of the estimators obtained using their method. His method described and made use of the dependency information among survival times, and the simulation results showed that their method performed better than the WLW method.

In practice, it is always difficult for investigators to identify significant covariates when the number of covariates is large, and variable selection studies increasingly involve the analysis of survival data with high-dimensional covariates to solve this difficulty. Tibshirani [8] proposed the application of the $L_1$ penalty function (LASSO) in the Cox model; Zou [9] proposed the adaptive Lasso (AdpLASSO), which Zhang [10] studied in the Cox model; and Zhang [11] proposed the minimax concave penalty. Several studies have focused on variable selection for multivariate survival data. For example, Cai [12] proposed a variable selection method for a growing number of regression coefficients; Liu [13] proposed a multivariate varying-coefficient hazard model; and Sun [14] developed a variable selection technique for multivariate interval-censored data.

Fan and Li [15] proposed a new type of penalty function called Smoothly Clipped Absolute Deviation (SCAD). The SCAD method combines characteristics of the LASSO and least squares. It compresses the coefficients of the model through a penalty function such that some coefficients are compressed to 0, thereby achieving variable selection, and the larger coefficients can also achieve asymptotically unbiased estimates. Moreover, Fan and Li proposed the oracle property and then introduced the SCAD penalty function into the Cox model [16].

In this paper, we aim to further improve the partition method to propose a new variable selection method for multivariate survival data. Based on the partition method, we make better use of the dependency information among survival times compared to the WLW method. Moreover, we directly achieve the purpose of variable selection using estimating equations. We construct our method with the SCAD penalty function and prove that the obtained estimators possess the oracle property. Numerical studies show that the proposed method performs well.

The rest of this paper is organized as follows. In Section 2, we present the used notation and assumptions. Then, we introduce our method and present its asymptotic and oracle properties in Section 3. We address implementation issues in Section 4, while simulations and an application of our method to real data are given in Section 5 and Section 6. We leave the proofs in the Appendix A.

2. Notation and Assumptions

$T_{ki},Z_{ki}\left(t\right),C_{ki},{\tilde{T}}_{ki},N_{ki}\left(t\right),Y_{ki}\left(t\right),{\lambda }_{ki}\left(t\right|Z_{ki})$, and ${\lambda }_{k0}\left(t\right)$ are the same as in Section 1. To facilitate the notation, let

(3)$M_{ki}\left(t\right)=N_{ki}\left(t\right)-{\int }_0^t{Y}_{ki}\left(u\right){\lambda }_{k0}\left(u\right)exp\left({\beta }_{k0}^T{Z}_{ki}\left(u\right)\right)du,i=1,\dots ,n,k=1,\dots ,K,$

which is a martingale with respect to $\sigma$-filtration ${\mathcal{F}}_{t,ki}=\sigma \{N_{ki}\left(s\right),Y_{ki}\left(s\right),Z_{ki}\left(s\right):0⩽s⩽t\}$.

For $k=1,\dots ,K$, let

$\begin{array}{cc}\hfill & S_k^{\left(d\right)}({\beta }_k,t)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{Y}_{ki}\left(t\right)Z_{ki}{\left(t\right)}^{\otimes d}exp\left({\beta }_k^T{Z}_{ki}\left(t\right)\right),s_k^{\left(d\right)}=E{S}_k^{\left(d\right)}({\beta }_k,t),\hfill \\ \hfill & {\overline{Z}}_k({\beta }_k,t)={\displaystyle \frac{S_k^{\left(1\right)}({\beta }_k,t)}{S_k^{\left(0\right)}({\beta }_k,t)}},{\mu }_k({\beta }_k,t)={\displaystyle \frac{s_k^{\left(1\right)}({\beta }_k,t)}{s_k^{\left(0\right)}({\beta }_k,t)}},\hfill \\ \hfill & V_k({\beta }_k,t)={\displaystyle \frac{S_k^{\left(2\right)}({\beta }_k,t)}{S_k^{\left(0\right)}({\beta }_k,t)}}-{\overline{Z}}_k{({\beta }_k,t)}^{\otimes 2},v_k({\beta }_k,t)={\displaystyle \frac{s_k^{\left(2\right)}({\beta }_k,t)}{s_k^{\left(0\right)}({\beta }_k,t)}}-{\mu }_k{({\beta }_k,t)}^{\otimes 2},\hfill \end{array}$

where E denotes expectation. For a column vector $\alpha$, ${\alpha }^{\otimes 0}=1$, ${\alpha }^{\otimes 1}=\alpha$, and ${\alpha }^{\otimes 2}=\alpha {\alpha }^T$.

The preliminary estimators of the baseline cumulative hazard functions are given by

${\widehat{\Lambda }}_{k0}\left(t\right)={\int }_0^t{\displaystyle \frac{d{N}_{k·}\left(u\right)}{n{S}_k^{\left(0\right)}({\widehat{\beta }}_k,u)}},$

where $N_{k·}\left(t\right)={\sum }_{i=1}^n{N}_{ki}\left(t\right)$. We obtain the following estimated martingales:

${\widehat{M}}_{ki}\left(t\right)=N_{ki}\left(t\right)-{\int }_0^t{Y}_{ki}\left(u\right)exp({\widehat{\beta }}_k^T{Z}_{ki}\left(u\right))d{\widehat{\Lambda }}_{k0}\left(u\right),i=1,\dots ,n,k=1,\dots ,K,$

where ${\widehat{\beta }}_1,\dots ,{\widehat{\beta }}_K$ are the WLW estimators.

We consider only events up to $\tau$ such that $Pr({\tilde{T}}_k⩾\tau )=0$ for all k. Let

$U_k({\beta }_k,t)=\sum _{i=1}^n{\int }_0^t(Z_{ki}\left(u\right)-{\overline{Z}}_k({\beta }_k,u))d{N}_{ki}\left(u\right),t\in [0,\tau ].$

Cui [7] introduced a partition. For the kth event, partition $[0,\tau ]$ into $L_k$ intervals is expressed as follows:

$0=t_0^{\left(k\right)}<t_1^{\left(k\right)}<\cdots <t_{L_k}^{\left(k\right)}=\tau .$

Let $L=L_1+\cdots +L_K$ define partition $\Pi$, as follows:

(4)$\Pi =\{(t_{i_1}^{\left(1\right)},\cdots ,t_{i_K}^{\left(K\right)}):t_{i_k}^{\left(k\right)}⩾0,i=0,1,\dots ,L_k,1⩽k⩽K\}.$

Following Cui [7], we break $U_k({\beta }_k,t)$ into $L_k$ pieces as follows:

$\Delta U_{\Pi }^{\left(k\right)}\left({\beta }_k\right)={(\Delta U_{\Pi }^{\left(k\right)}{({\beta }_k,t_1^{\left(k\right)})}^T,\dots ,\Delta U_{\Pi }^{\left(k\right)}{({\beta }_k,t_{L_k}^{\left(k\right)})}^T)}^T,$

where, for $l=1,\dots ,L_k$,

$\Delta U_{\Pi }^{\left(k\right)}({\beta }_k,t_l^{\left(k\right)})=U_{\Pi }^{\left(k\right)}({\beta }_k,t_l^{\left(k\right)})-U_{\Pi }^{\left(k\right)}({\beta }_k,t_{l-1}^{\left(k\right)})=\sum _{i=1}^n{\int }_{t_{l-1}^{\left(k\right)}}^{t_l^{\left(k\right)}}(Z_{ki}\left(t\right)-{\overline{Z}}_k({\beta }_k,t))d{N}_{ki}\left(t\right).$

Let $\Delta U_{\Pi }={(\Delta U_{\Pi }^{\left(1\right)}{\left({\beta }_1\right)}^T,\dots ,\Delta U_{\Pi }^{\left(K\right)}{\left({\beta }_K\right)}^T)}^T$. Cui [7] introduced the following notations:

$\begin{array}{cc}\hfill \Delta {\widehat{\Psi }}_{\Pi }\left(\widehat{\beta }\right)& =-{\displaystyle \frac{1}{n}}{\displaystyle \frac{\partial \Delta U_{\Pi }\left(\beta \right)}{\partial {\beta }^T}}{|}_{\beta =\widehat{\beta }}\hfill \\ \hfill & ={\left(\begin{array}{cccc}\Delta {\widehat{\Psi }}_{11}\left({\widehat{\beta }}_1\right)& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ \Delta {\widehat{\Psi }}_{1L_1}\left({\widehat{\beta }}_1\right)& 0& \cdots & 0\\ 0& \Delta {\widehat{\Psi }}_{21}\left({\widehat{\beta }}_2\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& \Delta {\widehat{\Psi }}_{2L_2}\left({\widehat{\beta }}_2\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \Delta {\widehat{\Psi }}_{K1}\left({\widehat{\beta }}_K\right)\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \Delta {\widehat{\Psi }}_{K{L}_K}\left({\widehat{\beta }}_K\right)\end{array}\right)}_{pL\times pK},\hfill \end{array}$

$\begin{array}{cc}\hfill & \Delta {\widehat{\Psi }}_{kl}({\widehat{\beta }}_k)=\Delta {\widehat{\Psi }}_{\Pi }({\widehat{\beta }}_k,t_l^{\left(k\right)})-\Delta {\widehat{\Psi }}_{\Pi }({\widehat{\beta }}_k,t_{l-1}^{\left(k\right)})={\int }_{t_{l-1}^{\left(k\right)}}^{t_l^{\left(k\right)}}{V}_k({\widehat{\beta }}_k,u)S_k^{\left(0\right)}({\widehat{\beta }}_k,u)d{\widehat{\Lambda }}_{k0}\left(u\right),\hfill \\ \hfill & \Delta {\widehat{\Sigma }}_{\Pi }(\widehat{\beta })={({\widehat{b}}_{l{l}^{\prime }}^{\left(kj\right)}(\widehat{\beta }))}_{1⩽l⩽L_k,1⩽l^{\prime }⩽L_j,1⩽k,j⩽K},\hfill \\ \hfill & ({\widehat{b}}_{l{l}^{\prime }}^{\left(kj\right)}(\widehat{\beta }))={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left[{\int }_{t_{l-1}^{\left(k\right)}}^{t_l^{\left(k\right)}}(Z_{ki}\left(t\right)-{\overline{Z}}_k({\widehat{\beta }}_k,t))d{\widehat{M}}_{ki}\left(t\right)\right]\left[{\int }_{t_{l^{\prime }-1}^{\left(j\right)}}^{t_{l^{\prime }}^{\left(j\right)}}{(Z_{ji}\left(t\right)-{\overline{Z}}_j({\widehat{\beta }}_j,t))}^{T}d{\widehat{M}}_{j}i\left(t\right)\right],\hfill \\ \hfill & \widehat{\beta }={\widehat{\beta }}_W.\hfill \end{array}$

Cui [7] focused on the situation in which ${\Pi }_n$ varies with sample size n. Let

$\Delta t_{{\Pi }_n}^{\left(k\right)}=\underset{1⩽l⩽L_k}{max}\{t_l^{\left(k\right)}-t_{l-1}^{\left(k\right)}\},$

such that

$\underset{n\to \infty }{lim}\Delta t_{{\Pi }_n}^{\left(k\right)}⟶0,k=1,\dots ,K.$

In addition, he proposed the following estimating equations:

(5)$\Delta {\widehat{\Psi }}_{{\Pi }_n}{({\widehat{\beta }}_W)}^T\Delta {\widehat{\Sigma }}_{{\Pi }_n}{({\widehat{\beta }}_W)}^{-1}\Delta U_{{\Pi }_n}\left(\beta \right)=0.$

We impose the following conditions:

$Pr(Y_{ki}\left(t\right)=1,t\in [0,\tau ])>0,i=1,\dots ,n,k=1,\dots ,K.$

$\left|Z_{kij}\left(0\right)\right|+{\int }_0^t\left|d{Z}_{kij}\left(t\right)\right|<B$ a.s. for $j=1,\dots ,p,i=1,\dots ,n,k=1,\dots ,K$, and some constant $B<\infty .$

${\lambda }_{k0}\left(t\right)$ is a continuous function of $t\in [0,\tau ]$, and there exist constants $C_1>0$ and $\gamma >0$ such that ${\lambda }_{k0}\left(t\right)⩾C{t}^{\gamma },t\in [0,\tau ].$

There exists a neighborhood $\mathcal{B}$ of ${\mathit{\beta }}_{\mathbf{0}}$ such that, for $d=0,1,2$ and $k=1,\dots ,K$,

$\underset{t\in [0,\tau ],\mathit{\beta }\in \mathcal{B}}{sup}\parallel S_k^{\left(d\right)}({\mathit{\beta }}_{\mathit{k}},t)-s_k^{\left(d\right)}({\mathit{\beta }}_{\mathit{k}},t)\parallel =O_p\left({\displaystyle \frac{1}{\sqrt{n}}}\right).$

$s_k^{\left(d\right)}({\mathit{\beta }}_{\mathit{k}},t)$ (k = 1, …, K; d = 0, 1, 2) is a continuous function of $\mathit{\beta }\in \mathcal{B}$ uniformly in $t\in [0,\tau ]$ and is bounded on $\mathcal{B}\times [0,\tau ]$, $s_k^{\left(0\right)}⩾s_0>0(k=1,\dots ,K)$, $s_0$ is a constant, $\mathit{\beta }\in \mathcal{B}$, $t\in [0,\tau ]$, and

$s_k^{\left(1\right)}({\mathit{\beta }}_k,t)={\displaystyle \frac{\partial }{\partial {\mathit{\beta }}_{\mathit{k}}}}{s}_k^{\left(0\right)}({\mathit{\beta }}_k,t),s_k^{\left(2\right)}({\mathit{\beta }}_k,t)={\displaystyle \frac{{\partial }^2}{\partial {\mathit{\beta }}_{\mathit{k}}^{\mathbf{2}}}}{s}_k^{\left(0\right)}({\mathit{\beta }}_k,t),$

for $k=1,\dots ,K,\mathit{\beta }\in \mathcal{B}$, and $t\in [0,\tau ]$.

$v_k({\beta }_k,t)(k=1,\dots ,K)$ is a positive definite matrix on $\mathcal{B}\times [0,\tau ]$,

$\underset{\mathcal{B}\times [0,\tau ]}{inf}{\lambda }_{min}\left(v_k({\beta }_k,t)\right)={\lambda }_{min}^{\left(k\right)}>0,$

and ${\lambda }_{min}(·)$ represents the minimum eigenvalue of the matrix.

For all sufficiently large n, there exists $\Delta {\widehat{\Sigma }}_{{\Pi }_n}{\left({\widehat{\beta }}_W\right)}^{-1}$. We use ${\eta }_n$ to denote the partition index for partition ${\Pi }_n$, ${\eta }_n$ is an increasing positive sequence, and there exists a constant $C_3>C_2>0$ such that

$C_2{\eta }_n⩽min\{L_1,\dots ,L_K\}⩽max\{L_1,\dots ,L_K\}⩽C_3{\eta }_n,$

where we assume that ${\eta }_n\to \infty$, and ${\eta }_n^{4+2\gamma }/\sqrt{n}\to 0$.

Let $\Delta {\Psi }_{\Pi }^T\left(\beta \right)=-{\displaystyle \frac{1}{n}}{\displaystyle \frac{\partial \Delta U_{\Pi }\left(\beta \right)}{\partial {\beta }^T}}$ and $\Delta {\Sigma }_{\Pi }\left(\beta \right)={(b_{l{l}^{\prime }}^{\left(kj\right)}\left(\beta \right))}_{1⩽l⩽L_k,1⩽l^{\prime }⩽L_j,1⩽k,j⩽K}$. Let

$A\left(\mathit{\beta }\right)=\Delta {\Psi }_{{\Pi }_n}^T\left(\mathit{\beta }\right)\Delta {\Sigma }_{{\Pi }_n}^{-1}\left(\mathit{\beta }\right)\Delta {\Psi }_{{\Pi }_n}\left(\mathit{\beta }\right),$

and assume that there exists a constant $C_4$ such that

${\lambda }_{max}\left\{A\left(\mathit{\beta }\right)\right\}<C_4<\infty ,$

Assume that the penalty function $p_{\lambda }\left(\left|{\beta }_j\right|\right)$ satisfies

$\underset{n\to \infty }{lim}\underset{{\beta }_j\to 0+}{lim\; inf}{p}_{\lambda }^{\prime }\left({\beta }_j\right)/\lambda >0$

for all $j=1,\dots ,d_n$. Furthermore, we assume that there exists a constant $C_4$ such that, for nonzero ${\theta }_1$ and ${\theta }_2$, $|p_{\lambda }^{″}\left({\theta }_1\right)-p_{\lambda }^{″}\left({\theta }_2\right)|⩽C_5|{\theta }_1-{\theta }_2|$.

Remark 1.

Conditions 1 and 2 are also used by Andersen et al. [17]. Conditions 4–8, which are adapted from Cai et al. [18] and Cui et al. [7], guarantee the asymptotic normality of the penalized partition estimator. Condition 3 is satisfied for most commonly used distributions in survival analysis [19]. Condition 9 is also used by Cai et al. [12]

3. Main Results

3.1. Construction of Estimators

We introduce the penalized partition-estimating Equation (PPEE) as follows:

(6)$\Delta {\widehat{\Psi }}_{{\Pi }_n}{({\widehat{\beta }}_W)}^T\Delta {\widehat{\Sigma }}_{{\Pi }_n}{({\widehat{\beta }}_W)}^{-1}\Delta U_{{\Pi }_n}\left(\mathit{\beta }\right)-n{\mathit{p}}_{\mathit{\lambda }}^{\mathbf{\prime }}\left(\left|\mathit{\beta }\right|\right)sgn\left(\mathit{\beta }\right)=0$

where ${\mathit{p}}_{\mathit{\lambda }}^{\mathbf{\prime }}\left(\left|\mathit{\beta }\right|\right)={(p_{\lambda }^{\prime }\left(\left|{\beta }_1\right|\right),\dots ,p_{\lambda }^{\prime }\left(\left|{\beta }_d\right|\right))}^T$, and $p_{\lambda }(·)$ is a penalty function. We consider the differential form of the SCAD penalty proposed by Fan and Li [15,20] defined by

$p_{\lambda }^{\prime }\left(\left|\theta \right|\right)=\lambda \left\{I(\left|\theta \right|<\lambda )+{\displaystyle \frac{{(a\lambda -\left|\theta \right|)}_{+}}{(a-1)\lambda }}I(\left|\theta \right|⩾\lambda )\right\},$

for some $a>2$.

Then, we can obtain the estimators $\widehat{\mathit{\beta }}$ by solving the penalized partition-estimating Equation (6). The SCAD penalty function involves two tuning parameters, a and $\lambda$. We will explain how to choose the parameters and how to obtain the estimators $\widehat{\mathit{\beta }}$ in Section 4.

3.2. Asymptotic and Oracle Properties of the Proposed Estimator

Fan and Li proposed the oracle property [21], which means that the estimator has the same limiting distribution as an estimator that knows the true model a priori. In this section, we will provide the asymptotic properties for $\widehat{\mathit{\beta }}$ and show that it achieves the oracle property.

We consider the situation mentioned by Cai [12], where the regression coefficient varies with the sample size n, that is, $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_{d_n})}^T$, where $d_n$ tends to ∞ as $n\to \infty$ and $d_n^4/n\to 0$. Let ${\mathit{\beta }}_{\mathbf{0}}\equiv {({\beta }_{01},\dots ,{\beta }_{0d_n})}^T$ denote the true value of $\mathit{\beta }$. Furthermore, we let ${\mathit{\beta }}_{\mathbf{10}}$ and ${\mathit{\beta }}_{\mathbf{20}}$, respectively, denote the nonzero and zero components of ${\mathit{\beta }}_{\mathbf{0}}$. Without loss of generality, we write ${\mathit{\beta }}_{\mathbf{0}}={({\mathit{\beta }}_{\mathbf{10}}^{\mathit{T}},{\mathit{\beta }}_{\mathbf{20}}^{\mathit{T}})}^T$ and suppose that ${\beta }_{0j}\ne 0$ for $j⩽s_n$ and ${\beta }_{0j}=0$ for $j>s_n$, which means that ${\mathit{\beta }}_{\mathbf{10}}$ consists of the $s_n$ nonzero components of ${\mathit{\beta }}_{\mathbf{0}}$. Let

(7)$a_n=max\{\left|p_{{\lambda }_n}^{\prime }\left(\left|{\beta }_{0j}\right|\right)\right|:{\beta }_{0j}\ne 0\},$

and

(8)$b_n=max\{\left|p_{{\lambda }_n}^{″}\left(\left|{\beta }_{0j}\right|\right)\right|:{\beta }_{0j}\ne 0\}.$

In this section, we will use ${\lambda }_n$ instead of $\lambda$ to emphasize that $\lambda$ depends on n. For simplicity, we let $\mathit{f}\left(\mathit{\beta }\right)=\Delta {\widehat{\Psi }}_{{\Pi }_n}{\left({\widehat{\beta }}_W\right)}^T\Delta {\widehat{\Sigma }}_{{\Pi }_n}{\left({\widehat{\beta }}_W\right)}^{-1}\Delta U_{{\Pi }_n}\left(\mathit{\beta }\right)$ and $\mathit{g}\left(\mathit{\beta }\right)=\mathit{f}\left(\mathit{\beta }\right)-n{\mathit{p}}_{{\mathit{\lambda }}_{\mathit{n}}}^{\mathbf{\prime }}\left(\left|\mathit{\beta }\right|\right)sgn\left(\mathit{\beta }\right)$.

Remark 2.

Here, we discuss the compatibility of the given parameters. We first impose the condition that the SCAD penalty function possesses the following property: for nonzero fixed θ, $lim\sqrt{n/d_n}{p}_{{\lambda }_n}^{\prime }\left(\right|\theta \left|\right)=0$ and $lim{p}_{{\lambda }_n}^{″}\left(\right|\theta \left|\right)=0$. This can be satisfied by appropriately choosing the regularization parameter ${\lambda }_n$. If we choose ${\lambda }_n\to 0$ and $\sqrt{n/d_n}{\lambda }_n\to \infty$, the above property holds because $\sqrt{n/d_n}{p}_{{\lambda }_n}^{\prime }=p_{{\lambda }_n}^{″}\left(\right|\theta \left|\right)=0$ for $\theta \ne 0$; thus, $a_n\to 0$ and $b_n\to 0$ can be satisfied. Furthermore, we can obtain that, for any given constant M, ${inf}_{\left|\theta \right|⩽M{n}^{-1/2}}{p}_{{\lambda }_n}^{\prime }\left(\right|\theta \left|\right)={\lambda }_n$, which means that condition (I) can be satisfied. Therefore, the conditions given in the upcoming theorem will not contradict each other.

Based on Cui’s research [7], we have the following theorem and lemma:

Theorem 1.

Under conditions (1)–(8), let $\left\{{\Pi }_n\right\}$ be the partition sequence that satisfies the following condition:

$\underset{n\to \infty }{lim}\Delta t_{{\Pi }_n}^{\left(k\right)}⟶0,k=1,\dots ,K.$

Then, there exists a matrix W such that

$\underset{n\to \infty }{lim}\Delta {\Psi }_{{\Pi }_n}{\left({\beta }_0\right)}^T\Delta {\Sigma }_{{\Pi }_n}{\left({\beta }_0\right)}^{-1}\Delta {\Psi }_{{\Pi }_n}\left({\beta }_0\right)=W,$

and

${\displaystyle \frac{1}{\sqrt{n}}}\mathit{f}\left({\mathit{\beta }}_{\mathbf{0}}\right)\stackrel{d}{\to }N(0,W).$

Lemma 1.

Under conditions 1-8, we can obtain that

$\parallel {\displaystyle \frac{1}{n}}{\mathit{f}}^{\mathbf{\prime }}\left(\mathit{\beta }\right)+A\left(\mathit{\beta }\right)\parallel =o_p\left(1\right).$

The above theorem and lemma have been proven by Cui [7] and will not be repeated in this paper. With Lemma 1, we can obtain the following theorem, which shows that there exists a penalized partition estimate that converges at the rate $O_p\left\{\sqrt{d_n}(n^{-1/2}+a_n)\right\}$.

Theorem 2.

Under conditions 1–8, if $a_n\to 0$, $b_n\to 0$, $d_n^4/n\to 0$, and ${\eta }_n^{4+2\gamma }/\sqrt{n}\to 0$ as $n\to \infty$, there exists an approximate zero-crossing $\widehat{\mathit{\beta }}$ of $\mathit{g}\left(\mathit{\beta }\right)$ such that $\parallel \widehat{\mathit{\beta }}-{\mathit{\beta }}_{\mathbf{0}}\parallel =O_p\left\{\sqrt{d_n}(n^{-1/2}+a_n)\right\}$.

From Theorem 2, if $a_n=O\left(n^{-1/2}\right)$, which can be achieved by selecting the appropriate ${\lambda }_n$, then there exists a $root-(n/d_n)$-consistent approximate zero-crossing of $\mathit{g}\left(\mathit{\beta }\right)$. Let

$\Sigma =diag\left\{p_{{\lambda }_n}^{″}\left(\left|{\beta }_{01}\right|\right),\cdots ,p_{{\lambda }_n}^{″}\left(\left|{\beta }_{0s_n}\right|\right)\right\},$

and

$\mathit{b}=\left(p_{{\lambda }_n}^{\prime }\left(\left|{\beta }_{01}\right|\right)sgn\left({\beta }_{01}\right),\cdots ,p_{{\lambda }_n}^{\prime }\left(\left|{\beta }_{0s_n}\right|\right)sgn\left({\beta }_{0s_n}\right)\right).$

Then, we can obtain Theorem 3, which shows that the proposed estimator achieves the oracle property.

Theorem 3.

Under conditions 1-9, if $a_n\to 0$, $b_n\to 0$, $d_n^4/n\to 0$, and ${\eta }_n^{4+2\gamma }/\sqrt{n}\to 0$ as $n\to \infty$, and if ${\lambda }_n\to 0$, $\sqrt{n/d_n}{\lambda }_n\to \infty$, and $a_n=O\left(n^{-1/2}\right)$, then under the conditions of Theorem 2, with probability tending to 1, there exists a $root-(n/d_n)$-consistent approximate zero-crossing $\widehat{\mathit{\beta }}={({\widehat{\mathit{\beta }}}_1^T,{\widehat{\mathit{\beta }}}_2^T)}^T$ in Theorem 2 such that

1. 

(Sparsity)${\widehat{\mathit{\beta }}}_2=0$;

2. 

(Asymptotic normality)

(9) $\sqrt{n}(A_{11}+\Sigma )\left\{{\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10}+{(A_{11}+\Sigma )}^{-1}\mathit{b}\right\}\stackrel{d}{\to }N(0,W_{11}),$

where $A_{11}$ and $W_{11}$ are the first $s_n\times s_n$ sub-matrices of $A({\mathit{\beta }}_{10},\mathbf{0})$ and $W({\mathit{\beta }}_{10},\mathbf{0})$, respectively.

Remark 3.

The two theorems above indicate that with the SCAD penalty function, which means $a_n=0$, $\mathit{b}=0$, and $\Sigma =0$ for sufficiently large n, and we have

$\sqrt{n}{A}_{11}({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10})\stackrel{d}{\to }N(0,W_{11}).$

Thus, ${\widehat{\mathit{\beta }}}_1$ possesses the same sampling property as the oracle estimate. The oracle estimator knows $\mathit{\beta }$ in advance, and ${\widehat{\mathit{\beta }}}_2=0$ is the same as that. Hence, the penalized partition estimator that we propose achieves the oracle property.

4. Implementation

4.1. Solution of Penalized Partition-Estimating Equation

In Section 3, we construct the penalized partition-estimating equation. Here, we provide the method for solving the equation. First, we need to establish a reasonable partition ${\Pi }_n$. We make the number of partitions corresponding to each event the same, that is, $L_1=\cdots =L_K=L_{\ast }$. To ensure that the penalized partition-estimating equation is effective, we need to ensure that each interval of the partition contains a certain number of event times or failure times. Hence, a reasonable partition method is as follows: for the $k$–th event, sort the failure times of each subject from small to large, and use the $n/L_{\ast },2n/L_{\ast },\dots ,(L_{\ast }-1)n/L_{\ast }$-th failure time as the cut-point. If $kn/L_{\ast }$ is not an integer ($k=1,\dots ,L_{\ast }-1$), then round it to an integer.

As the derivative function of the SCAD penalty ${\left(p_{\lambda }\left(\left|\beta \right|\right)\right)}^{\prime }=p_{\lambda }^{\prime }\left(\left|\beta \right|\right)sgn\left(\beta \right)$ is discontinuous near 0, we need to obtain an approximate differential penalty function ${\tilde{p}}_{\lambda }^{\prime }\left(\theta \right)$. We rewrite the penalized partition-estimating equation as follows:

(10)$\mathit{f}\left(\mathit{\beta }\right)-n{\tilde{\mathit{p}}}_{\mathit{\lambda }}^{\prime }\left(\left|\mathit{\beta }\right|\right)sgn\left(\mathit{\beta }\right)=0.$

Then, we need to solve the equation above to obtain the penalized partition estimator. In this study, we use the gradient descent method to solve (10).

4.2. Abnormal Condition Handling Within Zero Neighborhood

We need to address the issue mentioned above, that is, the approximation of the differential penalty function when a component ${\beta }_j$ of $\mathit{\beta }$ approaches 0 (abnormal condition handling within the zero neighborhood). In practice, for a very small $ϵ,0<ϵ\ll \lambda$, it can be assumed that $\beta =0$ if $\beta \in (-ϵ,ϵ)$. For the SCAD penalty $p_{\lambda }\left(\left|\beta \right|\right)$, its derivative function is discontinuous near 0. We use a linear function in a small neighborhood $(-ϵ,ϵ)$ near 0 to obtain the approximate differential penalty function ${\tilde{p}}_{\lambda }^{\prime }$:

${\tilde{p}}_{\lambda }^{\prime }\left(\left|\beta \right|\right)=\lambda \left\{{\displaystyle \frac{\left|\beta \right|}{ϵ}}I(\left|\beta \right|<ϵ)+I(ϵ⩽\left|\beta \right|<\lambda )+{\displaystyle \frac{{(a\lambda -\left|\beta \right|)}_{+}}{(a-1)\lambda }}I(\left|\beta \right|⩾\lambda )\right\}.$

This approximation has little impact on the SCAD penalty function. Since there are random errors in the data and model, it will not have a significant impact on the model’s performance as long as the impact of the small neighborhood $(-ϵ,ϵ)$ on the SCAD penalty does not exceed the random errors.

However, there are still issues with this approximation. When ${\beta }_j\in (-ϵ,0)\cup (0,ϵ)$, though ${\beta }_j$ can be approximated as 0 in the analysis of the estimation, the variation in ${\beta }_j$ may cause significant variation in $g\left(\beta \right)=f\left(\beta \right)-n{\tilde{p}}_{\lambda }^{\prime }\left(\left|\beta \right|\right)sgn\left(\beta \right)$, thus resulting in variation in other components of $\beta$, except for ${\beta }_j$, which will make the zero point unstable. If we let all ${\beta }_j\in (-ϵ,ϵ)$ be 0 and prohibit their variation, then although $g\left(\beta \right)$ will remain stable, this method will cause another deficiency, that is, all ${\beta }_j$ close to 0 can no longer escape $(-ϵ,ϵ)$ once they fall into this interval. This can cause a deviation in the zero point and even make it impossible to solve for the zero point. Therefore, we introduce the following Algorithm 1 to solve this problem.

Algorithm 1 Abnormal condition handling within zero neighborhood

Input:   $\beta ,g$ Output:   $\beta ,g$   1:. $I_j=I(\beta \in (-ϵ,ϵ))$,   $j=1,\dots ,d$,   $I=(I_1,\dots ,I_d)$ $//$I denotes a vector representing whether $\beta$ is within a small neighborhood around 0   2:. $\Delta \beta =\alpha g$ $//$$\Delta \beta$ denotes the variation in $\beta$ in the gradient descent method.   3:. $M=$ which $(I==1)$ $//$M denotes the index of the components of $\beta$ that fall into the small neighborhood around 0.   4:. if $M\ne NULL$ then   5:.    if ${max}_{j\in M}|\beta -\Delta {\beta }_j|<ϵ$ then   6:.         ${\beta }_j=0$,   $j\in M$   7:.         $g\left(\beta \right)=f\left(\beta \right)-n{\tilde{p}}_{\lambda }^{\prime }\left(\beta \right)$   8:.         $g_j=0$,   $j\in M$   9:.    end if 10:. end if 11:. return $\beta ,g$

It can be seen that this method can set a small neighborhood with “attraction” near 0. When ${\beta }_j$ falls into this neighborhood, it will temporarily be “attracted” to 0. If $g\left(\beta \right)$ has a large value in this component and can leave this special zero neighborhood in the next step, then this component will not be approximated as 0. Otherwise, when all the components “attracted” to 0 cannot leave this special zero neighborhood in the next step, we will approximate these components to 0, calculate $g\left(\beta \right)$, and perform gradient descent on these components. The next transformation in the algorithm is set to 0.

Condition handling within the zero neighborhood was also mentioned by Fan and Li [21]. They used a quadratic function to approximate the penalty function, resulting in the situation where coefficients are forced to be 0 and can no longer leave the small neighborhood [21]. Our proposed method for abnormal condition handling within the zero neighborhood can solve this problem by setting an appropriate threshold for the “attraction” near 0.

4.3. Tuning Parameter Selection

In order to achieve effective variable selection in the solution of (10), it is necessary to first select appropriate regularization parameters. In their simulations, Fan and Li found that, when $a\approx 3.7$, the SCAD method provided the best variable selection and coefficient estimation performance in penalty least-squares estimation [21]. Subsequently, many scholars continued to use $a\approx 3.7$ in the field of multivariate survival analysis (see Cai [12], liu [13], and Cai [22]). They pointed out that, from a Bayesian statistical point of view, it is suggested that $a\approx 3.7$ be used, as the Bayes risk cannot be reduced much with other choices of a. In this paper, we also take $a\approx 3.7$. Therefore, in the following parameter selection, the value of $\lambda$ needs to be chosen.

Next, we use k-fold cross-validation (k-fold CV) to select $\lambda$. Specifically, we use 10 folds. As it is hard to determine the concrete representation of the objective function of our method, we cannot establish the cross-validation statistic. However, we can apply the 10-fold CV method to the WLW method with the SCAD penalty and select the appropriate regularization parameter from it. The basis for doing so is that the WLW method is a special case of partition estimation when the number of partitions $L_{\ast }$ is set to 1. Therefore, the estimators obtained using the WLW method and the partition-estimating equation have values and errors on the same scale. Hence, it is reasonable to assume that the corresponding optimal regularization parameters have values of similar scale and will not have significant differences.

5. Simulation Study

In this section, we describe our evaluation of the performance of the proposed method based on the results of simulation studies. Raftery [23,24] introduced a bivariate exponential model in 1984. Based on their model, we generate bivariate survival data as follows:

$\begin{array}{cc}\hfill T_1& =\left((1-q_1)Y_1-I_{1}log\left(U\right)\right)exp(-{\beta }_1^{T}Z),\hfill \\ \hfill T_2& =\left((1-q_2)Y_2-I_{2}log(1-U)\right)exp(-{\beta }_2^{T}Z),\hfill \end{array}$

where the covariate Z is a p-dimensional 0–1 vector that indicates the presence or absence of features corresponding to each variable, $T_1$ and $T_2$ are the bivariate failure times, $Y_1$ and $Y_2$ follow a standard exponential distribution, $U\sim U(0,1)$, $I_1$ and $I_2$ are independent random variables $P(I_i=1)=q_i(i=1,2)$, and $q_1$$q_2$ are the pre-set parameters used to adjust the correlation between binary exponential distributions. In addition, $Y_1$, $Y_2$, U, $I_1$, and $I_2$ are independent.

5.1. Different Numbers of Partitions

This section presents simulation experiments to evaluate the performance of the proposed method in parameter estimation for different numbers of partitions $L_{\ast }$. We considered settings with $p=8$, $n=1000$, and $q_1=q_2=q=0.98$. In each simulation, each component of ${\beta }_1$ and ${\beta }_2$ had a probability of $({\displaystyle \frac{1}{6}},{\displaystyle \frac{1}{6}},{\displaystyle \frac{1}{3}},{\displaystyle \frac{1}{6}},{\displaystyle \frac{1}{6}})$ of taking on a value of $(-0.5,-0.25,0,0.25,0.5)$. This experiment did not set a censor time, and the simulation was repeated 1000 times. Figure 1 shows the simulation results for different numbers of partitions. When the number of partitions is 1, it is the WLW method. The results show that the larger the number of partitions, the smaller the MSE of parameter estimation, indicating that the larger the number of partitions, the better the performance of our method. Initially, as the number of partitions increases, the MSE of parameter estimation will significantly decrease, and the advantages of the proposed method will rapidly expand. When the number of partitions is large, as the number of partitions increases, the rate of MSE reduction in parameter estimation slows down, and the improvement in the proposed method’s performance in parameter estimation will become inapparent.

5.2. Different Correlations Between Various Events in Multivariate Survival Data

This section presents simulation experiments to evaluate the performance of the proposed method in the parameter estimation and variable selection for different correlations between various events in multivariate survival data. We considered settings with p = 8, n = 200, 1000, and $L_{\ast }=8$. Let the true values be ${\mathit{\beta }}_{\mathbf{0}}={(0.5,-0.25,0,0,0,0,0.25,0,0,-0.5,0,0.25,0,0,0,0)}^T$. In a simulation study, Cui observed that the partition-estimating function method was superior to the WLW method [7] when there was a significant correlation between various events in multivariate survival data. We let $q_1=q_2=q=0.98,0.8,0.25$. Censoring times were generated from a uniform distribution over $(0,c)$, and we chose $c=1,5$ to change the censoring rate. Each configuration had 1000 replications.

We assess the performance of our method using the model error (ME), similar to Fan and Li [16].

(11)$ME=E{\{exp(-\widehat{\mathit{\beta }}Z)-exp(-{\mathit{\beta }}_{\mathbf{0}}Z)\}}^2.$

In addition, we use the oracle estimator ${\widehat{\mathit{\beta }}}_{OR}$ to define the relative model error (RME) as follows:

(12)$RME={\displaystyle \frac{ME({\widehat{\mathit{\beta }}}_{OR})}{ME(\widehat{\mathit{\beta }})}}.$

We compared our method (PPEE) with the WLW-with-LASSO method and the WLW-with-SCAD method. The results of 1000 simulated datasets are given in Table 1. The column labeled “RME” provides the median of 1000 RMEs, and the column labeled “C” reports the average number of coefficients correctly estimated as 0, while the column “IC” presents the average number of the coefficients erroneously estimated as 0. From Table 1, we can see that our method performs well in terms of variable selection. When there is a high correlation between $T_1$ and $T_2$, our method performs particularly well. While when there is a weak correlation, the performance is slightly better than that of WLW+SCAD but still much better than that of WLW+LASSO. This indicates that our method performs better when the correlation between $T_1$ and $T_2$ is greater. In practical analysis, there is a certain correlation between failure times in multivariate survival data, which is also the advantage of our method over classical methods. In addition, none of the methods set the coefficients of significant variables to 0; thus, there is no under-fitting situation. Moreover, when the sample size n is large, none of the methods set the coefficients of significant variables to 0, and there is no under-fitting situation.

6. The Colon Cancer Study

In this section, we report the results of applying our method to the dataset collected in the Colon Cancer Study [25]. This study was initiated in the 1980s and included 929 patients with stage C disease randomly assigned to observation, levamisole alone, or levamisole combined with fluorouracil. The observation (Obs), levamisole alone (Lev), and levamisole combined with fluorouracil (Lev + 5-FU) groups comprise 315, 310, and 304 patients, respectively. There are multiple failure outcomes, such as the time to cancer recurrence and the survival time. By the end of the study, 155 patients in Obs, 144 in Lev, and 103 in Lev + 5-FU had experienced recurrences, and 114, 109, and 78 had died, respectively. We are interested in the following risk factors: sex, where 1 is male and 0 is female; age; obstruction of the colon by the tumor (obstruct); adherence to nearby organs (adhere); differentiation of the tumor (differ); perforation of the colon (perfor); number of lymph nodes with detectable cancer (nodes); extent of local spread (extent); more than four positive lymph nodes (node4), coded as 1 if true and 0 otherwise; and time from surgery to registration (surg), coded as 1 for long and 0 for short. In addition, similar to Lin [2], we created two dummy variables as follows: Lev, coded as 1 for the levamisole-alone treatment group and 0 for others, and Lev + 5FU, coded as 1 for the levamisole combined with the fluorouracil treatment group and 0 otherwise.

Table 2 shows the estimated coefficients and standard errors for the Colon Cancer Study using different methods, including the unpenalized method (UNM), LASSO, and our method (PPEE). From Table 2, we can see that our method only keeps five significant variables. In the column “LASSO”, we observe that certain variables may have an impact on one failure event but not on another. According to our method, “Lev + 5FU”, “Extent”, and “Node4” all have a significant impact on the two failure events. “Lev + 5FU” has a negative impact on both death and recurrence, which is consistent with Moertel’s study [26]; that is, levamisole combined with fluorouracil is effective in reducing the mortality rate of colon cancer. In addition, Lin [2] found that “extent” and “node4” increased the risk of colon cancer, and our method is consistent with this.

7. Discussion and Conclusions

We proposed a penalized partition-estimating equation for variable selection in multivariate survival data. The partition-estimating equation was originally proposed by Cui [7]. We developed Cui’s method and proposed the penalized partition-estimating equation to simultaneously estimate the parameters and select variables. Compared with the classic Cox regression method, our method can more effectively select variables and estimate coefficients, which is reflected in our simulation experiments. Moreover, our method makes use of the dependency information among survival times and performs better when the correlation between failure times is greater, which is also considered as an advantage of our method over classical methods, as there is a certain correlation between failure times in multivariate survival data in practical analysis. Moreover, we proved the asymptotic and oracle properties of the proposed method.

Future studies can supplement this work. In this study, our method performed well when failure events have a strong correlation. Therefore, one future task is to consider whether our method still performs well when there is no obvious correlation between events and, if not, explore how this problem can be solved. Another interesting task is ultrahigh-dimensional variable selection.

Author Contributions

Methodology, W.C.; Validation, W.C. and W.T.; formal analysis, W.C. and W.T.; writing—original draft preparation, W.T.; writing—review and editing, W.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figure and Tables

Figure 1 Mean squared error (MSE) of parameter estimates for $\beta$ for different number of partitions.

Table 1

Relative model error.

Method	q	RME	Number of Zeros		RME	Number of Zeros
			C	IC		C	IC
		c = 1			c = 5
n = 200
LASSO		0.473	10.458	0.036	0.531	10.513	0.006
SCAD	0.98	0.677	10.703	0.024	0.745	10.783	0.004
PPEE		0.733	10.712	0.024	0.801	10.785	0.004
LASSO		0.469	10.471	0.030	0.518	10.543	0.005
SCAD	0.8	0.691	10.710	0.030	0.753	10.770	0.006
PPEE		0.703	10.705	0.036	0.783	10.776	0.004
LASSO		0.479	10.443	0.036	0.520	10.412	0.004
SCAD	0.25	0.682	10.702	0.026	0.747	10.773	0.002
PPEE		0.691	10.719	0.024	0.756	10.775	0.002
n = 1000
LASSO		0.583	10.498	0.000	0.612	10.532	0.000
SCAD	0.98	0.752	10.811	0.000	0.796	10.848	0.000
PPEE		0.813	10.819	0.000	0.857	10.855	0.000
LASSO		0.590	10.473	0.000	0.593	10.528	0.000
SCAD	0.8	0.743	10.802	0.000	0.790	10.842	0.000
PPEE		0.792	10.811	0.000	0.824	10.847	0.000
LASSO		0.591	10.482	0.000	0.588	10.541	0.000
SCAD	0.25	0.755	10.807	0.000	0.781	10.837	0.000
PPEE		0.761	10.811	0.000	0.792	10.841	0.000

Table 2

Estimated coefficients and standard errors for the Colon Cancer Study using different methods.

Effect	UNM	LASSO	SCAD	PPEE
	Recurrence
Lev	−0.026 (0.111)	0.000	0.000	0.000
Lev + 5FU	−0.499 (0.122)	−0.441 (0.108)	−0.416 (0.108)	−0.428 (0.107)
Sex	−0.138 (0.096)	0.000	0.000	0.000
Age	−0.003 (0.004)	0.000	0.000	0.000
Obstruct	0.194 (0.119)	0.061 (0.095)	0.050 (0.134)	0.048 (0.103)
Perfor	0.211 (0.257)	0.000	0.000	0.000
Adhere	0.161 (0.130)	0.028 (0.137)	0.028 (0.136)	0.000
Nodes	0.038 (0.015)	0.037 (0.017)	0.000	0.000
Differ	0.153 (0.098)	0.118 (0.108)	0.036 (0.106)	0.024 (0.105)
Extent	0.451 (0.119)	0.414 (0.120)	0.393 (0.119)	0.532 (0.116)
Surg	0.240 (0.104)	0.072 (0.110)	0.084 (0.108)	0.000
Node4	0.591 (0.141)	0.641 (0.103)	0.772 (0.146)	0.751 (0.106)
	Death
Lev	−0.041 (0.114)	0.000	0.000	0.000
Lev + 5FU	−0.362 (0.122)	−0.294 (0.109)	−0.209 (0.108)	−0.226 (0.107)
Sex	0.007 (0.097)	0.000	0.000	0.000
Age	0.008 (0.004)	0.006 (0.004)	0.002 (0.004)	0.000
Obstruct	0.269 (0.120)	0.118 (0.135)	0.098 (0.135)	0.094 (0.131)
Perfor	0.017 (0.270)	0.000	0.000	0.000
Adhere	0.170 (0.131)	0.138 (0.145)	0.130 (0.141)	0.135 (0.126)
Nodes	0.044 (0.015)	0.043 (0.014)	0.000	0.000
Differ	0.138 (0.101)	0.106 (0.110)	0.007 (0.106)	0.003 (0.110)
Extent	0.446 (0.118)	0.420 (0.114)	0.377 (0.111)	0.427 (0.112)
Surg	0.240 (0.106)	0.021 (0.113)	0.079 (0.110)	0.000
Node4	0.667 (0.143)	0.641 (0.128)	0.657 (0.143)	0.899 (0.153)

Appendix A. Proofs

For simplicity, we let $F^{\prime }\left(\mathit{\beta }\right)=\mathit{f}\left(\mathit{\beta }\right)$ and $G\left(\mathit{\beta }\right)=F\left(\mathit{\beta }\right)-n{\sum }_{j=1}^{d_n}{p}_{\lambda }\left(\left|{\beta }_j\right|\right)$.

Proof of Theorem 2.

To prove Theorem 2, it is sufficient to show that $\widehat{\mathit{\beta }}$ is a local maximizer of $\mathit{G}\left(\mathit{\beta }\right)$. To prove this, let ${\alpha }_n=\sqrt{d_n}(n^{-1/2}+a_n)$; it is sufficient to show that for any $ϵ>0$, there exists a large constant C such that(A1)$$P\left\{\underset{\parallel \mathit{u}\parallel =C}{sup}\mathit{G}({\mathit{\beta }}_{\mathbf{0}}+{\alpha }_n\mathit{u})<\mathit{G}\left({\mathit{\beta }}_{\mathbf{0}}\right)\right\}⩾1-ϵ.$$

This means that the probability of a local maximum existing in the ball $\{{\mathit{\beta }}_{\mathbf{0}}+{\alpha }_n\mathit{u}:\parallel \mathit{u}\parallel ⩽C\}$ is at least $1-ϵ$. Therefore, a local maximizer exists such that $\parallel \widehat{\mathit{\beta }}-{\mathit{\beta }}_{\mathbf{0}}\parallel =O_p\left({\alpha }_n\right)$, and $\widehat{\mathit{\beta }}$ is a zero-crossing of $\mathit{g}\left(\mathit{\beta }\right)$.

Note that $p_{\lambda }\left(0\right)=0$ and $p_{\lambda }(·)⩾0$; then,$$\begin{array}{cc}\hfill L_n\left(\mathit{u}\right)& =\mathit{G}({\mathit{\beta }}_{\mathbf{0}}+{\alpha }_n\mathit{u})-\mathit{G}\left({\mathit{\beta }}_{\mathbf{0}}\right)\hfill \\ \hfill & ⩽\{F({\mathit{\beta }}_{\mathbf{0}}+{\alpha }_n\mathit{u})-F\left({\mathit{\beta }}_{\mathbf{0}}\right)\}-n\sum _{j=1}^{s_n}\{p_{{\lambda }_n}\left(|{\beta }_{0j}+{\alpha }_n{u}_j|\right)-p_{{\lambda }_n}\left(\left|{\beta }_{0j}\right|\right)\}\hfill \\ \hfill & \triangleq I_1+I_2.\hfill \end{array}$$

Firstly, we consider $I_1$. Let$${\mathit{f}}^{\mathbf{\prime }}\left(\mathit{\beta }\right)={\displaystyle \frac{\partial }{\partial {\mathit{\beta }}^T}}\mathit{f}\left(\mathit{\beta }\right).$$

According to the Taylor expansion,$$\begin{array}{cc}\hfill I_1& ={\alpha }_n{\mathit{u}}^T\mathit{f}\left({\mathit{\beta }}_{\mathbf{0}}\right)+{\displaystyle \frac{1}{2}}{\alpha }_n^2{\mathit{u}}^T{\mathit{f}}^{\mathbf{\prime }}\left({\mathit{\beta }}_{\mathit{n}}^{\mathbf{\ast }}\right)\mathit{u}\hfill \\ \hfill & \triangleq I_{11}+I_{12},\hfill \end{array}$$ where ${\mathit{\beta }}_{\mathit{n}}^{\mathbf{\ast }}$ lies between ${\mathit{\beta }}_{\mathbf{0}}$ and ${\mathit{\beta }}_{\mathbf{0}}+{\alpha }_n\mathit{u}$. According to the Cauchy–Schwarz inequality,$$I_{11}={\alpha }_n{\mathit{u}}^T\mathit{f}\left({\mathit{\beta }}_{\mathbf{0}}\right)⩽{\alpha }_n\parallel \mathit{f}\left({\mathit{\beta }}_{\mathbf{0}}\right)\parallel \parallel \mathit{u}\parallel =O_p\left({\alpha }_n\sqrt{n{d}_n}\right)\parallel \mathit{u}\parallel =O_p\left(n{\alpha }_n^2\right)\parallel \mathit{u}\parallel .$$

Next, we consider $I_{12}$. We can obtain$${\mathit{f}}^{\mathbf{\prime }}\left(\mathit{\beta }\right)=-n\Delta {\widehat{\Psi }}_{{\Pi }_n}{\left(\widehat{{\beta }_W}\right)}^T\Delta {\widehat{\Sigma }}_{{\Pi }_n}{\left(\widehat{{\beta }_W}\right)}^{-1}\Delta {\widehat{\Psi }}_{{\Pi }_n}\left(\mathit{\beta }\right).$$

From Lemma 1, we can obtain(A2)$$\parallel {\displaystyle \frac{1}{n}}{\mathit{f}}^{\prime }\left(\mathit{\beta }\right)+A\left(\mathit{\beta }\right)\parallel =o_p\left(1\right).$$

Consequently,$$I_{12}=-{\displaystyle \frac{1}{2}}n{\alpha }_n^2{\mathit{u}}^{T}A\left({\mathit{\beta }}_{\mathbf{0}}\right)\mathit{u}\{1+op\left(1\right)\}.$$

Thus, $I_{12}$ dominates $I_{11}$ uniformly in $\parallel \mathit{u}\parallel =C$ when choosing a sufficiently large C. Then, using the Taylor expansion, we can obtain$$\begin{array}{cc}\hfill I_2& =-\sum _{j=1}^{s_n}\left\{n{\alpha }_n{p}_{{\lambda }_n}^{\prime }\left(\left|{\beta }_{0j}\right|\right)sgn\left({\beta }_{0j}\right)u_j+n{\alpha }_n^2{p}_{{\lambda }_n}^{″}\left(\left|{\beta }_{0j}\right|\right)u_j^2\{1+o\left(1\right)\}\right\}\hfill \\ \hfill & \triangleq -I_{21}-I_{22}.\hfill \end{array}$$

Using the Cauchy–Schwarz inequality,$$\left|I_{21}\right|⩽\sqrt{s_n}n{\alpha }_n{a}_n\parallel \mathit{u}\parallel ⩽\sqrt{d_n}n{\alpha }_n{a}_n\parallel \mathit{u}\parallel ⩽n{\alpha }_n^2\parallel \mathit{u}\parallel .$$

Furthermore,$$\left|I_{22}\right|=n{\alpha }_n^2\sum _{j=1}^{s_n}\left|p_{{\lambda }_n}^{″}\left(\left|{\beta }_{0j}\right|\right)\right|u_j^2\{1+o\left(1\right)\}⩽2b_{n}n{\alpha }_n^2{\parallel \mathit{u}\parallel }^2.$$

$b_n\to 0$, and therefore, we can choose a sufficiently large C so that $I_{12}$ also dominates both $I_{21}$ and $I_{22}$. This means that (A1) holds, and the proof is complete. □

To prove Theorem 3, we introduce a lemma first.

Lemma A1.

Under the conditions of Theorem 3, with a probability tending to 1, for any given ${\mathit{\beta }}_1$ satisfying $\parallel {\mathit{\beta }}_1-{\mathit{\beta }}_{\mathbf{10}}\parallel =O_p\left(\sqrt{d_n/n}\right)$ and any constant C, the following holds: (A3) $$G\{{({\mathit{\beta }}_1^T,\mathbf{0})}^T\}=\underset{\parallel {\mathit{\beta }}_{\mathbf{2}}\parallel ⩽C\sqrt{d_n/n}}{max}G\{{({\mathit{\beta }}_1^T,{\mathit{\beta }}_2^T)}^T\}.$$

Proof. 

It is sufficient to show that for any ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$ satisfying the conditions above, $\partial G\left(\mathit{\beta }\right)/\partial {\beta }_j$ and ${\beta }_j$ have different signs for all ${\beta }_j\in (-C\sqrt{d_n/n},C\sqrt{d_n/n}),j=s_n+1,\dots ,d_n$. We have$${\displaystyle \frac{\partial G\left(\mathit{\beta }\right)}{\partial {\beta }_j}}={\displaystyle \frac{\partial F\left(\mathit{\beta }\right)}{\partial {\beta }_j}}-n{p}_{{\lambda }_n}^{\prime }\left(\left|{\beta }_j\right|\right)sgn\left({\beta }_j\right).$$

Using the Taylor expansion, we have$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial F\left(\mathit{\beta }\right)}{\partial {\beta }_j}}& ={\displaystyle \frac{\partial F\left({\mathit{\beta }}_{\mathbf{0}}\right)}{\partial {\beta }_j}}-n\sum _{l=1}^{d_n}{A}_{jl}\left({\mathit{\beta }}_{\mathbf{0}}\right)({\beta }_l-{\beta }_{0l})\{1+o_p\left(1\right)\}\hfill \\ \hfill & =I-II,\hfill \end{array}$$ where $A_{jl}\left({\mathit{\beta }}_0\right)$ is the $(j,l)$-element of $A\left({\mathit{\beta }}_{\mathbf{0}}\right)$. According to the standard argument, it follows that$$I=O_p\left(\sqrt{n}\right)=O_p\left(\sqrt{n{d}_n}\right).$$

Next, we consider $II$. According to the Cauchy–Schwarz inequality, it follows that$$\parallel II\parallel ⩽n{\left\{\sum _{l=1}^{d_n}{A}_{jl}^2\left({\mathit{\beta }}_0\right)\right\}}^{1/2}\parallel \mathit{\beta }-{\mathit{\beta }}_{\mathbf{0}}\parallel ,$$ and then we have$$\sum _{j=1}^{d_n}{A}_{jl}^2\left({\mathit{\beta }}_{\mathbf{0}}\right)={\mathit{e}}_j^{T}A{\left({\mathit{\beta }}_{\mathbf{0}}\right)}^{T}A\left(\mathit{\beta }\right){\mathit{e}}_{\mathit{j}}⩽{\lambda }_{max}^2\left\{A\left({\mathit{\beta }}_{\mathbf{0}}\right)\right\}=O\left(1\right),$$ where ${\mathit{e}}_j$ is a $d_n\times 1$ vector whose jth element is 1 and others are 0. Since ${\mathit{\beta }}_{\mathbf{1}}-{\mathit{\beta }}_{\mathbf{10}}=O_p\left(\sqrt{d_n/n}\right)$ and $\parallel {\mathit{\beta }}_{\mathbf{2}}\parallel ⩽C\sqrt{d_n/n}$, we can obtain $\parallel \mathit{\beta }-{\mathit{\beta }}_{\mathbf{0}}\parallel =O_p\left(\sqrt{d_n/n}\right)$. Hence,$$II=O_p\left(\sqrt{n{d}_n}\right),{\displaystyle \frac{\partial F\left(\mathit{\beta }\right)}{\partial {\beta }_j}}=O_p\left(\sqrt{n{d}_n}\right).$$

Thus,$${\displaystyle \frac{\partial G\left(\mathit{\beta }\right)}{\partial {\beta }_j}}=n{\lambda }_n\left[-{\lambda }_n^{-1}{p}_{{\lambda }_n}^{\prime }\left(\left|{\beta }_j\right|\right)sgn\left({\beta }_j\right)+O_p\{1/\left({\lambda }_n\sqrt{n/d_n}\right)\}\right].$$

Given the assumption that ${\lambda }_n\sqrt{n/d_n}\to \infty$, $O_p\{1/\left({\lambda }_n\sqrt{n/d_n}\right)\}=o_p\left(1\right)$, the sign of ${\displaystyle \frac{\partial G\left(\mathit{\beta }\right)}{\partial {\beta }_j}}$ is determined by the sign of ${\beta }_j$. Thus, we complete the proof. □

Proof of Theorem 3.

We can immediately prove part 1 through Lemma A1. To prove part 2, it is sufficient to show that(A4)$${\mathit{f}}_{\mathbf{1}}\left({\mathit{\beta }}_{\mathbf{0}}\right)-n(A_{11}+\Sigma )\left\{\widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}}+{(A_{11}+\Sigma )}^{-1}\mathit{b}\right\}+o_p\left(\sqrt{n}\right)=0,$$ where ${\mathit{f}}_{\mathbf{1}}\left({\mathit{\beta }}_{\mathbf{0}}\right)$ denotes the first $s_n$ components of $\mathit{f}\left({\mathit{\beta }}_{\mathbf{0}}\right)$. Then,$$\sqrt{n}(A_{11}+\Sigma )\left\{{\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10}+{(A_{11}+\Sigma )}^{-1}\mathit{b}\right\}={\displaystyle \frac{1}{\sqrt{n}}}{\mathit{f}}_{\mathbf{1}}\left({\mathit{\beta }}_{\mathbf{0}}\right)+o_p\left(1\right)$$

According to condition (F) and Cui [7], we have$${\displaystyle \frac{1}{\sqrt{n}}}{\mathit{f}}_{\mathbf{1}}\left({\mathit{\beta }}_{\mathbf{0}}\right)\stackrel{d}{\to }N(0,W_{11}).$$

Thus, we need to prove that (A4) is valid. There exists a $\widehat{{\mathit{\beta }}_{\mathbf{1}}}$ in Theorem 2 that is a $root-(n/d_n)$-consistent approximate zero-crossing of $\mathit{g}\left\{{({\mathit{\beta }}_{\mathbf{1}},\mathbf{0})}^T\right\}$, and it satisfies$$\partial G\{{({\widehat{\mathit{\beta }}}_1,\mathbf{0})}^T\}/\partial {\beta }_j=0,j=1,\dots ,s_n.$$

Using the Taylor expansion to $\partial G\{{({\widehat{\mathit{\beta }}}_1,\mathbf{0})}^T\}/\partial {\mathit{\beta }}_{\mathbf{1}}$ at ${\mathit{\beta }}_{\mathbf{10}}$, it follows that$${\mathit{f}}_{\mathbf{1}}\left({\mathit{\beta }}_{\mathbf{0}}\right)+{\mathit{f}}_{\mathbf{1}}^{\mathbf{\prime }}\left({\mathit{\beta }}_{\mathbf{1}}^{\mathbf{\ast }}\right)(\widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}}-n\mathit{b}-n{\Sigma }^{\ast }(\widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}}))=0,$$ where ${\mathit{\beta }}_{\mathbf{1}}^{\mathbf{\ast }}$ lies between $\widehat{{\mathit{\beta }}_{\mathbf{1}}}$ and ${\mathit{\beta }}_{\mathbf{10}}$, and ${\Sigma }^{\ast }=diag\left\{p_{{\lambda }_n}^{″}\left(\left|{\beta }_1^{\ast }\right|\right),\cdots ,p_{{\lambda }_n}^{″}\left(\left|{\beta }_{s_n}^{\ast }\right|\right)\right\}$. From (A2), we have$$\begin{array}{cc}\hfill \parallel \{{\mathit{f}}_1^{\prime }\left({\mathit{\beta }}_{\mathbf{1}}^{\mathbf{\ast }}\right)+n{A}_{11}\}(\widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}})\parallel & ⩽n\parallel n^{-1}{\mathit{f}}_1^{\prime }\left({\mathit{\beta }}_{\mathbf{1}}^{\mathbf{\ast }}\right)+A_{11}\parallel \parallel \widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}}\parallel \hfill \\ \hfill & ⩽n{o}_p\left(1\right)O_p\left(\sqrt{d_n/n}\right)\hfill \\ \hfill & =o_p\left(\sqrt{n/d_n}\right)\hfill \\ \hfill & =o_p\left(\sqrt{n}\right).\hfill \end{array}$$

From condition (I), it follows that$$\parallel n({\Sigma }^{\ast }-\Sigma )(\widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}})\parallel ⩽n{C}_5{\parallel \widehat{{\mathit{\beta }}_{\mathbf{1}}}-{\mathit{\beta }}_{\mathbf{10}}\parallel }^2=C_{4}n{O}_p(d_n/n)=O_p\left(d_n\right)=o_p\left(\sqrt{n}\right).$$

Therefore, (A4) holds. Thus, we complete the proof. □

References

1. Liang, K.Y.; Self, S.G.; Chang, Y.C. Modelling marginal hazards in multivariate failure time data. J. R. Stat. Soc. Ser. B Stat. Methodol.; 1993; 55, pp. 441-453. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1993.tb01914.x]

2. Lin, D.Y. Cox regression analysis of multivariate failure time data: The marginal approach. Stat. Med.; 1994; 13, pp. 2233-2247. [DOI: https://dx.doi.org/10.1002/sim.4780132105] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/7846422]

3. Spiekerman, C.F.; Lin, D.Y. Marginal regression models for multivariate failure time data. J. Am. Stat. Assoc.; 1998; 93, pp. 1164-1175. [DOI: https://dx.doi.org/10.1080/01621459.1998.10473777]

4. Wei, L.J.; Lin, D.Y.; Weissfeld, L. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. J. Am. Stat. Assoc.; 1989; 84, pp. 1065-1073. [DOI: https://dx.doi.org/10.1080/01621459.1989.10478873]

5. Cox, D.R. Regression models and life-tables. J. R. Stat. Soc. Ser. B Methodol.; 1972; 34, pp. 187-202. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1972.tb00899.x]

6. Cox, D.R. Partial likelihood. Biometrika; 1975; 62, pp. 269-276. [DOI: https://dx.doi.org/10.1093/biomet/62.2.269]

7. Cui, W.Q.; Ying, Z.L.; Zhao, L.C. A simple construction of optimal estimation in multivariate marginal Cox regression. Sci. China Math.; 2012; 55, pp. 1827-1857. [DOI: https://dx.doi.org/10.1007/s11425-012-4400-4]

8. Tibshirani, R. The lasso method for variable selection in the Cox model. Stat. Med.; 1997; 16, pp. 385-395. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3]

9. Zou, H. The adaptive lasso and its oracle properties. J. Am. Stat. Assoc.; 2006; 101, pp. 1418-1429. [DOI: https://dx.doi.org/10.1198/016214506000000735]

10. Zhang, H.H.; Lu, W. Adaptive Lasso for Cox’s proportional hazards model. Biometrika; 2007; 94, pp. 691-703. [DOI: https://dx.doi.org/10.1093/biomet/asm037]

11. Zhang, C.H. Nearly unbiased variable selection under minimax concave penalty. Ann. Statist.; 2010; 38, pp. 894-942. [DOI: https://dx.doi.org/10.1214/09-AOS729] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17244211]

12. Cai, J.W.; Fan, J.Q.; Li, R.Z.; Zhou, H.B. Variable selection for multivariate failure time data. Biometrika; 2005; 92, pp. 303-316. [DOI: https://dx.doi.org/10.1093/biomet/92.2.303] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19458784]

13. Liu, J.C.; Zhang, R.Q.; Zhao, W.H.; Lv, Y.Z. Variable selection in semiparametric hazard regression for multivariate survival data. J. Multivar. Anal.; 2015; 142, pp. 26-40. [DOI: https://dx.doi.org/10.1016/j.jmva.2015.07.015]

14. Sun, L.; Li, S.; Wang, L.; Song, X.; Sui, X. Simultaneous variable selection in regression analysis of multivariate interval-censored data. Biometrics; 2022; 78, pp. 1402-1413. [DOI: https://dx.doi.org/10.1111/biom.13548]

15. Fan, J.Q.; Li, R.Z. Variable Selection via Penalized Likelihood; Department of Statistics, UCLA: Los Angeles, CA, USA, 1999.

16. Fan, J.Q.; Li, R.Z. Variable selection for Cox’s proportional hazards model and frailty model. Ann. Stat.; 2002; 30, pp. 74-99. [DOI: https://dx.doi.org/10.1214/aos/1015362185]

17. Andersen, P.K.; Gill, R.D. Cox’s regression model for counting processes: A large sample study. Ann. Stat.; 1982; 10, pp. 1100-1120. [DOI: https://dx.doi.org/10.1214/aos/1176345976]

18. Cai, J.W. Hypothesis testing of hazard ratio parameters in marginal models for multivariate failure time data. Lifetime Data Anal.; 1999; 5, pp. 39-53. [DOI: https://dx.doi.org/10.1023/A:1009679032314]

19. Cui, W.Q. Analysis of Multivariate Survival Data by Marginal Proportional Hazards Regression Models. Ph.D. Thesis; University of Science and Technology of China: Hefei, China, 2004; (In Chinese)

20. Fan, J.Q.; Li, R.Z. New estimation and model selection procedures for semiparametric modeling in longitudinal data analysis. J. Am. Stat. Assoc.; 2004; 99, pp. 710-723. [DOI: https://dx.doi.org/10.1198/016214504000001060]

21. Fan, J.Q.; Li, R.Z. Variable selection via nonconcave penalized likelihood and its oracle properties. J. Am. Stat. Assoc.; 2001; 96, pp. 1348-1360. [DOI: https://dx.doi.org/10.1198/016214501753382273]

22. Cai, K. Bi-level Variable Selection and Dimension-reduction Methods in Complex Lifetime Data Analytics. Ph.D. Thesis; University of Calgary: Calgary, AB, Canada, 2019.

23. Raftery, A.E. A continuous multivariate exponential distribution. Commun. Stat. Theory Methods; 1984; 13, pp. 947-965. [DOI: https://dx.doi.org/10.1080/03610928408828733]

24. Raftery, A.E. Some properties of a new continuous bivariate exponential distribution. Stat. Decis. Suppl. Issue; 1985; 2, pp. 53-58.

25. Moertel, C.G.; Fleming, T.R.; Macdonald, J.S.; Haller, D.G.; Laurie, J.A.; Goodman, P.J.; Ungerleider, J.S.; Emerson, W.A.; Tormey, D.C.; Glick, J.H. . Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma. N. Engl. J. Med.; 1990; 322, pp. 352-358. [DOI: https://dx.doi.org/10.1056/NEJM199002083220602] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/2300087]

26. Moertel, C.G.; Fleming, T.R.; Macdonald, J.S.; Haller, D.G.; Laurie, J.A.; Tangen, C.M.; Ungerleider, J.S.; Emerson, W.A.; Tormey, D.C.; Glick, J.H. . Fluorouracil plus levamisole as effective adjuvant therapy after resection of stage III colon carcinoma: A final report. Ann. Intern. Med.; 1995; 122, pp. 321-326. [DOI: https://dx.doi.org/10.7326/0003-4819-122-5-199503010-00001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/7847642]