Variable Selection for Length-Biased and Interval-Censored Failure Time Data

Abstract

Length-biased failure time data occur often in various biomedical fields, including clinical trials, epidemiological cohort studies and genome-wide association studies, and their analyses have been attracting a surge of interest. In practical applications, because one may collect a large number of candidate covariates for the failure event of interest, variable selection becomes a useful tool to identify the important risk factors and enhance the estimation accuracy. In this paper, we consider Cox’s proportional hazards model and develop a penalized variable selection technique with various popular penalty functions for length-biased data, in which the failure event of interest suffers from interval censoring. Specifically, a computationally stable and reliable penalized expectation-maximization algorithm via two-stage data augmentation is developed to overcome the challenge in maximizing the intractable penalized likelihood. We establish the oracle property of the proposed method and present some simulation results, suggesting that the proposed method outperforms the traditional variable selection method based on the conditional likelihood. The proposed method is then applied to a set of real data arising from the Prostate, Lung, Colorectal and Ovarian cancer screening trial. The analysis results show that African Americans and having immediate family members with prostate cancer significantly increase the risk of developing prostate cancer, while having diabetes exhibited a significantly lower risk of developing prostate cancer.

1. Introduction

This work proposes a general variable selection method for the length-biased and interval-censored failure time data with the classical proportional hazards (PH) model. Interval-censored data arise when each failure time of interest cannot be measured accurately but is only known to lie in a certain time interval formed by periodical follow-ups [1]. Such data are frequently encountered in many scientific studies including clinical trials and epidemiological surveys, and their regression analysis has been discussed extensively in the literature, see [2,3,4,5,6,7,8] for details. Specifically, Zeng et al. [4], Wang et al. [6] and Zeng et al. [7] investigated inference procedures for the additive hazards, PH and transformation models, respectively.

In addition to interval censoring, left truncation is also frequently encountered in prospective cohort studies, inducing non-randomly selected samples from the target population. A typical example of left truncation occurs in the PLCO study where individuals with any of the PLCO cancers at the onset of the study were not enrolled [6,9]. In particular, when the truncation times follow the uniform distribution (also known as the length-biased or stationarity assumption), the left-truncated data reduce to the length-biased data discussed by many authors, including but not limited to Wang [10], Shen et al. [11] and Ning et al. [12].

The analysis of length-biased data under right censoring has been investigated extensively in the literature [11,13,14,15,16]. To name a few examples, Shen et al. [11] presented unbiased estimating equation approaches for the transformation and accelerated failure time models. Qin and Shen [13] developed an inverse weighted estimating equation approach for the PH model. Qin et al. [14] developed new expectation-maximization (EM) algorithms to estimate the survival function of the failure time. For the length-biased and interval-censored data, Gao and Chan [15] developed an EM algorithm for the PH model via two-stage data augmentation. Further, Shen et al. [16] considered the mixture PH model with a nonsusceptible or cured fraction.

In many practical applications, one may collect a large number of candidate covariates, but in general, only a few covariates are useful to model the failure time of interest. In such a case, penalized variable selection provides a useful tool to eliminate irrelevant variables and further enhance the estimation accuracy. Popular penalty functions include LASSO [17], SCAD [18], adaptive LASSO (ALASSO) [19], SICA [20], SELO [21], MCP [22] and BAR [23,24]. In particular, Fan et al. [25] provided a comprehensive review of variable selection methods and the corresponding algorithms. In recent years, machine learning-based methods have also gained considerable attention due to their great ability in identifying relevant features. To name a few examples, Garavand et al. [26] used clinical examination features and compared different machine learning algorithms in developing a model for the early diagnosis of coronary artery disease. Hosseini et al. [27] used blood microscopic images and a convolutional neural network algorithm for detecting and classifying B-acute lymphoblastic leukemia. Garavand et al. [28] conducted a systematic review of advanced techniques to facilitate the rapid diagnosis of coronary artery disease. Ghaderzadeh and Aria [29] conducted a systematic review of Artificial Intelligence techniques for COVID-19 detection.

Regarding the left-truncated failure time data, a number of variable selection methods have been proposed. In particular, Chen [30] considered the right-censored data and developed a variable selection method for the additive hazards model with covariate measurement errors. He et al. [31] also considered the right-censored data and performed variable selection with penalized estimating equations for the accelerated failure time model. Li et al. [32] developed a conditional likelihood-based variable selection method for left-truncated and interval-censored data under the PH model. However, it is worth noting that the work of Li et al. [32] only involved the ALASSO penalty and their method can be anticipated to lose some efficiency due to ignoring the distribution information of the truncation times.

In this paper, we offer an efficient penalized likelihood method to achieve variable selection in the PH model with length-biased and interval-censored data. Compared with the traditional conditional likelihood method in Li et al. [32], the proposed method yields an efficiency gain via fully taking into account the distribution information of the truncation times. In particular, to optimize the penalized likelihood function with an intractable form, we develop a penalized EM algorithm by introducing pseudo-left-truncated data and Poisson random variables. The proposed method is easy to implement and computationally stable and has desirable advantages over the variable selection method based on the penalized conditional likelihood [32]. An application to a real data set arising from the Prostate, Lung, Colorectal and Ovarian (PLCO) cancer screening trial demonstrates the practical usefulness of the proposed method.

The PLCO cancer screening trial is a large-scale multicenter trial conducted to screen the PLCO cancers and investigate cancer-related mortality. To date, motivated by the rich data structure in the PLCO database, various statistical methods have already been proposed in the literature. To name a few examples, Wang et al. [6] developed an EM algorithm to estimate the spline-based PH model with interval-censored data. Sun et al. [33] considered variable selection in a semiparametric nonmixture cure model with interval-censored data. Li and Peng [34] investigated instrumental variable estimation of complier causal treatment effect with interval-censored data. Withana Gamage et al. [35] considered the estimation of the PH model with left-truncated and arbitrarily interval-censored data.

The remainder of this paper is organized as follows. In Section 2, we first introduce the notation, assumption and corresponding likelihood. Section 3 presents the proposed penalized EM algorithm, and Section 4 establishes the oracle property of the proposed estimators. In Section 5, a simulation study is conducted to assess the variable selection performance and the estimation accuracy of the proposed method, followed by an application in Section 6. Some discussions and conclusions are given in Section 7. Section 8 provides several potential future research directions.

2. Notation, Model and Penalized Likelihood

For the target population, let $\tilde{T}$, $\tilde{A}$ and $\tilde{\mathit{Z}}$ denote the failure time of interest (e.g., the time to the onset of the failure event), the truncation time (e.g., the time to the study enrollment) and the p-dimensional vector of covariates, respectively. Given $\tilde{\mathit{Z}}$, the PH model specifies that the conditional cumulative hazard function of $\tilde{T}$ takes the form

(1)$\Lambda (t\mid \tilde{\mathit{Z}})=\Lambda \left(t\right)exp\left({\mathit{\beta }}^{\top }\tilde{\mathit{Z}}\right),$

where $\mathit{\beta }$ is the p-dimensional vector of the unknown regression coefficients and $\Lambda (·)$ is an unknown increasing cumulative baseline hazard function. Let d denote the number of nonzero components in $\mathit{\beta }$, and ${\mathit{\beta }}_{10}$ and ${\mathit{\beta }}_{20}$ denote all d nonzero and $p-d$ zero-component vectors, respectively.

Under the left-truncation scheme, only individuals with $\tilde{T}\ge \tilde{A}$ are enrolled in the study, and the failure time, truncation time and covariate vector that satisfy $\tilde{T}\ge \tilde{A}$ are denoted by T, A and $\mathit{Z}$, respectively. Then, we know that $(T,A,\mathit{Z})$ has the same joint distribution as $(\tilde{T},\tilde{A},\tilde{\mathit{Z}})$ given $\tilde{T}\ge \tilde{A}$ [36]. As mentioned above, if the truncation time $\tilde{A}$ is further assumed to follow the uniform distribution on $(0,\tau )$, where $(0,\tau )$ is the support of $\tilde{T}$, we have the length-biased sampling mechanism [10,14]. Let $f(t\mid \mathit{z})$ and $S(t\mid \mathit{z})$ denote the density and survival functions of $\tilde{T}$ given $\tilde{\mathit{Z}}=\mathit{z}$, respectively. Under the assumption that $\tilde{A}$ is independent of $\tilde{T}$ given $\tilde{\mathit{Z}}$, the joint density function of $(T,A)$ given $\mathit{Z}=\mathit{z}$ evaluated at $(t,a)$ is

(2)$\frac{h\left(a\right)f(t\mid \mathit{z})}{P(\tilde{T}\ge \tilde{A})}=\frac{h\left(a\right)f(t\mid \mathit{z})}{{\int }_0^{\tau }h\left(a\right)S(a\mid \mathit{z})da},$

where $h\left(a\right)$ denotes the density function of $\tilde{A}$ at time a and equals $1/\tau$ under the length-biased sampling scheme.

Consider a failure time study that recruits n subjects and each failure time suffers from interval censoring due to the periodical examinations for the occurrence of the failure event. For $i=1,\dots ,n$, denote by $T_i$, $A_i$ and ${\mathit{Z}}_i$ the failure time, truncation time and covariate vector of the ith subject in the study, respectively. We assume that there exists a sequence of examination times $A_i\le U_{i,1}<\cdots <U_{i,M_i}<U_{i,M_i+1}$ for subject i, where $M_i$ is a random positive integer, and define $U_{i,M_i+1}=\infty$. Let $L_i$ and $R_i$ denote the endpoints of the smallest interval that brackets $T_i$. That is, $T_i\in (L_i,R_i]$. Clearly, $T_i$ is left-censored if $L_i=A_i$, and $T_i$ is right-censored if $R_i=\infty$. Then, the observed data consist of $\mathcal{O}=\{{\mathcal{O}}_i=(L_i,R_i,A_i,{\mathit{Z}}_i);i=1,\cdots ,n\}$. Under the assumption that examination times are independent of the failure and truncation times given the covariates, the likelihood function based on the observed data can be written as

(3)$L(\mathit{\beta },\Lambda )=\prod _{i=1}^n\frac{S(L_i\mid {\mathit{Z}}_i)-I(R_i<\infty )S(R_i\mid {\mathit{Z}}_i)}{{\int }_0^{\tau }exp\left\{-\Lambda \left(a\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}da/\tau },$

where $S(t\mid {\mathit{Z}}_i)=exp\{-\Lambda \left(t\right)exp\left({\beta }^{\top }{\mathit{Z}}_i\right)\}$.

Essentially, the likelihood (3) is the product of the marginal likelihood and the conditional likelihood, that is, $L(\mathit{\beta },\Lambda )=L_M(\mathit{\beta },\Lambda )\times L_C(\mathit{\beta },\Lambda )$, where $L_M(\mathit{\beta },\Lambda )={\prod }_{i=1}^{n}S(A_i\mid {\mathit{Z}}_i)/{\mu }_{\mathit{\beta },\Lambda }\left({\mathit{Z}}_i\right),$ and $L_C(\mathit{\beta },\Lambda )={\prod }_{i=1}^n\{S(L_i\mid {\mathit{Z}}_i)-I(R_i<\infty )S(R_i\mid {\mathit{Z}}_i)\}/S(A_i\mid {\mathit{Z}}_i).$ In the above, $I(·)$ denotes the indicator function, ${\mu }_{\mathit{\beta },\Lambda }\left({\mathit{Z}}_i\right)={\tau }^{-1}{\int }_0^{\tau }exp\left\{-\Lambda \left(a\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}da$; $L_M(\mathit{\beta },\Lambda )$ is the marginal likelihood of $\{A_i;i=1,\cdots ,n\}$ given $\{{\mathit{Z}}_i;i=1,\cdots ,n\}$; and $L_C(\mathit{\beta },\Lambda )$ is the conditional likelihood given the $A_i$’s. Notably, the commonly used conditional likelihood method only utilizes $L_C(\mathit{\beta },\Lambda )$ for inference, which can be anticipated to lose some estimation efficiency because $L_M(\mathit{\beta },\Lambda )$ also involves the parameters in model (1).

For the nuisance function $\Lambda$, we propose to approximate it with a step function that has non-negative jumps at the unique examination times. Specifically, let $0<t_1<\cdots <t_{K_n}<\infty$ denote the ordered unique values of $\{(L_i,R_{i}I(R_i<\infty )),i=1,\cdots ,n\}$, where $K_n$ is an integer determined by the observed data. For $k=1,\cdots ,K_n$, denote by $f_k$ the non-negative jump size of $\Lambda$ at $t_k$. Then, we have $\Lambda \left(t\right)={\displaystyle \sum _{t_k\le t}}{f}_k$, and the likelihood function (3) can be rewritten as

(4)$L(\mathit{\beta },\mathit{F})=\prod _{i=1}^n\frac{exp\left\{-{\sum }_{t_k\le L_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}-I\left(R_i<\infty \right)exp\left\{-{\sum }_{t_k\le R_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}}{{\int }_0^{\tau }exp\left\{-{\sum }_{t_k\le a}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}da/\tau },$

where $\mathit{F}={(f_1,\dots ,f_{K_n})}^{\top }$.

To accomplish variable selection and estimate the nonzero parameters simultaneously, we propose to maximize the following penalized log-likelihood:

(5)$log{L}_P(\mathit{\beta },\mathit{F})=logL(\mathit{\beta },\mathit{F})-n\sum _{j=1}^p{p}_{\lambda }\left(\left|{\beta }_j\right|\right),$

where $p_{\lambda }\left(\right|{\beta }_j\left|\right)$ denotes a penalty function that depends on the tuning parameter $\lambda$. In what follows, we provide a general maximization procedure for (5) under various commonly adopted penalty functions, such as LASSO, ALASSO, SCAD, SELO, SICA, MCP and BAR [17,18,19,20,21,22,23,24]. Because the penalized log-likelihood (5) has an intractable form, performing direct maximization with some of the existing software is extremely difficult and unstable. This is also the case even without the penalty term as shown in Gao and Chan [15]. In the next section, we will propose a reliable and stable penalized EM algorithm to overcome this computation challenge.

3. Estimation Procedure

The proposed penalized EM algorithm involves two layers of data augmentation, which aims at simplifying the form of (4) and obtaining a tractable objective function. In the first stage of data augmentation, for the ith subject, we introduce a set of independent pseudo-truncated data, ${\mathcal{O}}^{\ast }=\{(T_{im}^{\ast },A_{im}^{\ast },{\mathit{Z}}_i):T_{im}^{\ast }<A_{im}^{\ast },m=1,\cdots ,N_i\}$, which are also referred to as “ghost data” [37], and the random integer $N_i$ follows a negative binomial distribution $N(1,{\pi }_i)$ with $E(N_i\mid {\mathcal{O}}_i)=(1-{\pi }_i)/{\pi }_i$, where

$\begin{array}{cc}\hfill {\pi }_i& =P(T_{im}^{\ast }\ge A_{im}^{\ast }\mid {\mathit{Z}}_i)=1-\sum _{k=0}^{K_n}P(T_{im}^{\ast }=t_k,A_{im}^{\ast }>t_k\mid {\mathit{Z}}_i)\hfill \\ \hfill & =1-\sum _{k=1}^{K_n}(1-t_k/\tau )f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)exp\left\{-\sum _{l=1}^k{f}_{l}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}.\hfill \end{array}$

Given $N_i=n_i$, let $N_{ik}={\sum }_{m=1}^{n_i}I(T_{im}^{\ast }=t_k)$ for each $k=1,\cdots ,K_n$, and then $(N_{i1},\cdots ,N_{i{K}_n})$ follows a multinomial distribution with probabilities $(p_{i1},\cdots ,p_{i{K}_n})$, where $p_{ik}=P(T_{im}^{\ast }=t_k\mid T_{im}^{\ast }<A_{im}^{\ast },{\mathit{Z}}_i)=q_{ik}/(1-{\pi }_i),$ and $q_{ik}=(1-t_k/\tau )f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)exp\{-{\sum }_{l=1}^k{f}_{l}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\}$. In the above, $\tau$ is the finite upper bound of the support of $\tilde{T}$ and can be specified as $t_{K_n}$ in practice [14]. After deleting some constants that are irrelevant to the parameters to be estimated, the augmented likelihood function based on $\{\mathcal{O},{\mathcal{O}}^{\ast }\}$ is

(6)$\begin{array}{cc}\hfill L_1(\mathit{\beta },\mathit{F})=& \prod _{i=1}^n[\left(exp\left\{-\sum _{t_k\le L_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}-I(R_i<\infty )exp\left\{-\sum _{t_k\le R_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}\right)\hfill \\ \hfill & \times \prod _{k=1}^{K_n}{q}_{ik}^{N_{ik}}].\hfill \end{array}$

In the second stage, for the ith subject, we introduce the independent latent variables $W_{ik}$ with $i=1,\cdots ,n$ and $k=1,\cdots ,K_n$, where $W_{ik}$ is the Poisson random variable with mean $f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)$. Then, the likelihood function (6) can be re-expressed with Poisson variables as

$\begin{array}{c}\hfill L_2(\mathit{\beta },\mathit{F})=\prod _{i=1}^n\left\{P\left(\sum _{t_k\le L_i}{W}_{ik}=0\right)P{\left(\sum _{L_i<t_k\le R_i}{W}_{ik}>0\right)}^{I(R_i<\infty )}\times \prod _{k=1}^{K_n}{q}_{ik}^{N_{ik}}\right\}.\end{array}$

Let $P(W_{ik}\mid f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right))$ be the probability mass function of $W_{ik}$ and $R_i^{\ast }=L_{i}I(R_i=\infty )+R_{i}I(R_i<\infty )$. By treating the latent variables $W_{ik}$’s and $N_{ik}$’s as observable, we have the following complete data likelihood

$L_c(\mathit{\beta },\mathit{F})=\prod _{i=1}^n\prod _{k=1}^{K_n}\left\{P{(W_{ik}\mid f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right))}^{I(t_k\le R_i^{\ast })}\times q_{ik}^{N_{ik}}\right\},$

where we require that ${\sum }_{t_k\le L_i}{W}_{ik}=0$ and ${\sum }_{L_i<t_k\le R_i}{W}_{ik}>0$ if $R_i<\infty$, and ${\sum }_{t_k\le L_i}{W}_{ik}=0$ if $R_i=\infty$.

Let $\mathit{\theta }={({\mathit{\beta }}^{\top },f_1,\cdots ,f_{K_n})}^{\top }$, and let ${\mathit{\theta }}^{\left(m\right)}$ be the update of $\mathit{\theta }$ at the mth iteration with $m\ge 0$. Based on $L_c(\mathit{\beta },\mathit{F})$, we can present the expectation step (E-step) and maximization step (M-step) of the proposed algorithm. In the E-step, we calculate the conditional expectations of $W_{ik}$ and $N_{ik}$ given the observed data and ${\mathit{\theta }}^{\left(m\right)}$ in the $log{L}_c(\mathit{\beta },\mathit{F})$. This step yields

(7)$\begin{array}{cc}\hfill Q(\mathit{\theta };{\mathit{\theta }}^{\left(m\right)})& =\sum _{i=1}^n\sum _{k=1}^{K_n}I(t_k\le R_i^{\ast })\left\{E\left(W_{ik}\right)log{f}_k+E\left(W_{ik}\right){\mathit{\beta }}^{\top }{\mathit{Z}}_i-f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}\hfill \\ \hfill & +\sum _{i=1}^n\sum _{k=1}^{K_n}E\left(N_{ik}\right)\left\{log{f}_k+{\mathit{\beta }}^{\top }{\mathit{Z}}_i-\sum _{l=1}^k{f}_{l}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}.\hfill \end{array}$

In particular, at the mth iteration of the algorithm, the expressions of the conditional expectations are given by

$\begin{array}{cc}\hfill E\left(W_{ik}\right)& =I(L_i<t_k\le R_i,R_i<\infty )\frac{f_k^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)}{1-exp\left\{-{\sum }_{L_i<t_l\le R_i}{f}_l^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)\right\}}\hfill \end{array}$

and

$E\left(N_{ik}\right)=\frac{(1-t_k/\tau )f_k^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)exp\left\{-{\sum }_{l=1}^k{f}_l^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)\right\}}{{\pi }_i^{\left(m\right)}},$

where

${\pi }_i^{\left(m\right)}=1-\sum _{k=1}^{K_n}(1-t_k/\tau )f_k^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)exp\left\{-\sum _{l=1}^k{f}_l^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)\right\}.$

For notational simplicity, we omitted the conditional arguments including the observed data and current estimates of the parameters in the above conditional expectations.

In the M-step of the algorithm, by solving $\partial Q(\mathit{\theta };{\mathit{\theta }}^{\left(m\right)})/\partial f_k=0$, we have a closed-form expression for $f_k$, which is given by

(8)$f_k^{(m+1)}=\frac{{\sum }_{i=1}^n\{I(t_k\le R_i^{\ast })E\left(W_{ik}\right)+E\left(N_{ik}\right)\}}{{\sum }_{i=1}^n\{I(t_k\le R_i^{\ast })+{\sum }_{l=k}^{K_n}E\left(N_{il}\right)\}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)},k=1,\cdots ,K_n.$

Next, by plugging (8) into (7), we have the following objective function that only involves the unknown parameter $\mathit{\beta }$:

$\begin{array}{cc}\hfill Q_{\mathrm{new}}(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})=\sum _{i=1}^n\sum _{k=1}^{K_n}& \left\{I(t_k\le R_i^{\ast })E\left(W_{ik}\right)+E\left(N_{ik}\right)\right\}\hfill \\ \hfill & \times \left\{{\mathit{\beta }}^{\top }{\mathit{Z}}_i-log\left[\sum _{i^{{}^{\prime }}=1}^n\left\{I(t_k\le R_{i^{{}^{\prime }}}^{\ast })+\sum _{l=k}^{K_n}E\left(n_{i^{{}^{\prime }}l}\right)\right\}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_{i^{{}^{\prime }}}\right)\right]\right\}.\hfill \end{array}$

To obtain the sparse estimator of $\beta$, we propose to minimize the following objective function

(9)$H(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})+\sum _{j=1}^p{p}_{\lambda }\left(\right|{\beta }_j\left|\right),$

where $H(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})=-\frac{1}{n}{Q}_{\mathrm{new}}(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})$.

For LASSO and ALASSO, the modified shooting algorithm given in Zhang and Lu [38] and others can be adopted to minimize (9). For the BAR, the closed-form solution for $\mathit{\beta }$ is available [39]. For other penalties, after using local linear approximation for $p_{\lambda }\left(\right|{\beta }_j\left|\right)$ [40], one can also adopt the modified shooting algorithm to minimize (9).

In summary, for the given $\lambda$ and initial estimator ${\mathit{\theta }}^{\left(0\right)}$, we repeat the E-step and M-step until the convergence criterion is satisfied, rendering the sparse estimators of the regression parameters. It is worth pointing out that the proposed algorithm is insensitive to the choices of the initial value ${\mathit{\theta }}^{\left(0\right)}$. In practice, one can simply set the initial value of each component in $\mathit{\beta }$ to 0 and the initial value of $f_k$ to $1/K_n$, for $k=1,\cdots ,K_n$. The proposed algorithm is declared to achieve convergence if the summation of the absolute differences in the estimates between two successive iterations is less than a small positive number, such as ${10}^{-3}$.

To select the optimal $\lambda$, we follow Li et al. [39] and others and adopt the BIC criterion, which is defined as

${\mathrm{BIC}}_{\lambda }=-2l(\widehat{\mathit{\beta }},\widehat{\mathit{F}})+d{f}_{\lambda }\times log\left(n\right),$

where $\widehat{\mathit{\beta }}$ is the final estimator of $\mathit{\beta }$, $\widehat{\mathit{F}}$ is the final estimator of $\mathit{F}$, $l(\widehat{\mathit{\beta }},\widehat{\mathit{F}})$ denotes the logarithm of (4) and $d{f}_{\lambda }$ is the total number of the nonzero estimates in $\widehat{\mathit{\beta }}$ and $\widehat{\mathit{F}}$. For a given set including candidate values of $\lambda$, the optimal $\lambda$ can be set as the one that yields the smallest BIC.

4. Asymptotic Properties

Without a loss of generality, we write $\mathit{\beta }={({\mathit{\beta }}_1^{\top },{\mathit{\beta }}_2^{\top })}^{\top }$, where ${\mathit{\beta }}_1$ includes the first d components of $\mathit{\beta }$ that are nonzero and ${\mathit{\beta }}_2$ consists of the remaining zero components. Denote the true value of $\mathit{\beta }$ by ${\mathit{\beta }}_0={({\mathit{\beta }}_{10}^{\top },{\mathit{\beta }}_{20}^{\top })}^{\top }$, where ${\mathit{\beta }}_{10}$ is the true value of ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_{20}$ is the true value of ${\mathit{\beta }}_2$. Let $\widehat{\mathit{\beta }}={({\widehat{\mathit{\beta }}}_1^{\top },{\widehat{\mathit{\beta }}}_2^{\top })}^{\top }$ be the estimator of $\mathit{\beta }$ obtained from the method proposed above, where ${\widehat{\mathit{\beta }}}_1$ denotes the estimate of ${\mathit{\beta }}_1$ and ${\widehat{\mathit{\beta }}}_2$ denotes the estimate of ${\mathit{\beta }}_2$. In what follows, we establish the asymptotic properties of $\widehat{\mathit{\beta }}$.

For any penalty function $p_{\lambda }\left(\right|{\beta }_j\left|\right)$ with tuning parameter $\lambda$, we let $\rho (x;\lambda )={\lambda }^{-1}{p}_{\lambda }\left(x\right)$ and assume that $p_{\lambda }\left(\right|{\beta }_j\left|\right)$ belongs to the function class $\mathcal{P}$ as considered in Lv and Fan [20]:

$\begin{array}{cc}\hfill \mathcal{P}=\{p_{\lambda }(·):\rho (x;\lambda )& \mathrm{is}\mathrm{an}\mathrm{increasing}\mathrm{function}\mathrm{of}x\in [0,\infty )\mathrm{and}\mathrm{has}\mathrm{a}\mathrm{continuous}\mathrm{derivative}\hfill \\ \hfill & {\rho }^{{}^{\prime }}(x;\lambda )\mathrm{for}\mathrm{any}x\in (0,\infty ),\mathrm{where}{\rho }^{{}^{\prime }}(x;\lambda )\in (0,\infty )\mathrm{is}\mathrm{independent}\mathrm{of}\lambda \}.\hfill \end{array}$

The class $\mathcal{P}$ is quite general and includes the penalty functions considered in this work. To establish the asymptotic properties of $\widehat{\mathit{\beta }}$, we need the following regularity conditions.

• (C1). The true regression parameter ${\mathit{\beta }}_0$ lies in a compact set $\mathcal{B}$ in ${\mathcal{R}}^p$, and the true cumulative baseline hazard function ${\Lambda }_0(·)$ is continuously differentiable and positive with ${\Lambda }_0^{{}^{\prime }}\left(t\right)>0$ for all $t\in [{\tau }_1,{\tau }_2]$, where $[{\tau }_1,{\tau }_2]$ is the union of the support of $(U_1,\cdots ,U_M)$ and $0<{\tau }_1<{\tau }_2<\infty$. In addition, we assume that $0<{\Lambda }_0\left({\tau }_2\right)<\infty$.

• (C2). The covariate vector $\mathit{Z}$ is bounded with probability one and the covariance matrix of $\mathit{Z}$ is positive definite.

• (C3). The number of examination times, M, is positive and $E\left(M\right)<\infty$. Additionally, there exists a positive constant $\eta$ such that $P(U_{m+1}-U_m\ge \eta \mid \mathit{Z},M)=1(m=1,\cdots ,M-1)$. Furthermore, there exists a probability measure $\mu$ in $[{\tau }_1,{\tau }_2]$ such that the bivariate distribution function of $(U_m,U_{m+1})$ conditional on $(M,\mathit{Z})$ is dominated by $\mu \times \mu$, and its Radon–Nikodym derivative, denoted by ${\tilde{f}}_m(u,v;M,\mathit{Z})(m=1,\cdots ,m-1)$, can be expanded to a positive and has twice-continuous derivatives with respect to u and v when $v-u>\eta$ and are continuously differentiable with respect to $\mathit{Z}$.

Conditions (C1) and (C2) are standard conditions in a failure time data analysis [7]. Condition (C3) pertains to the joint distribution of the examination times, which ensures that two adjacent examination times should be separated by at least $\eta$. Otherwise, the data may contain exactly observed failure times, requiring different theoretical treatments. Note that the conditions (C1)–(C3) are used in establishing the root-n consistency of the unpenalized maximum likelihood estimator of the regression vector [7], which is required for the penalty term in the penalized likelihood. Also, the conditions (C1)–(C3) ensure that the log profile likelihood ${\ell }_p\left(\mathit{\beta }\right)=\underset{\Lambda }{sup}logL(\mathit{\beta },\Lambda )$ has a quadratic expansion around ${\mathit{\beta }}_0$ [7].

Theorem 1 

(root-n consistency). Under conditions (C1) to (C3), if $\sqrt{n}\lambda =O_p\left(1\right)$, then $\parallel \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\parallel =O_p\left(n^{-1/2}\right)$, where $\parallel ·\parallel$ denotes the Euclidean norm for a given vector.

Theorem 2 

(oracle property). Under conditions (C1) to (C3), if $\sqrt{n}\lambda \to 0$ and $n\lambda \to \infty$, then $\widehat{\mathit{\beta }}$ has the following properties:

• 1. 

  (Sparsity) $\underset{n\to \infty }{lim}P\left({\widehat{\mathit{\beta }}}_2=\mathbf{0}\right)=1$;

• 2. 

  (Asymptotic normality) $\sqrt{n}\left({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10}\right)⇝N\left(\mathbf{0},{\tilde{I}}_{10}^{-1}\right)$, where ${\tilde{I}}_{10}$ is the upper-left $d\times d$ sub-matrix of the efficient Fisher information matrix for β, denoted by ${\tilde{I}}_0$.

Theorem 1 indicates that $\widehat{\mathit{\beta }}$ is consistent, and Theorem 2 (i) implies that $\widehat{\mathit{\beta }}$ is sparse and has the selection consistency property, that is, $\underset{n\to \infty }{lim}P(\{j:{\widehat{\beta }}_j\ne 0\}=\{1,\cdots ,d\})=1$. Theorem 2 (ii) implies that the estimators of the nonzero regression parameters are semiparametrically efficient. The detailed proofs of Theorems 1 and 2 under the $L_1$ penalty ALASSO are given in the Appendix A. For other penalty functions that belong to $\mathcal{P}$, one can prove the above two theorems with analogous techniques, which are omitted in this paper.

5. A Simulation Study

We conducted a simulation study to evaluate the finite-sample performance of the proposed penalized EM algorithm. We first assumed that there exist 10 covariates, following the marginal standard normal distribution with the pairwise correlation $corr(Z_j,Z_k)=0.5^{|j-k|}(j,k=1,\cdots ,10)$. We set the true value of $\mathit{\beta }$ denoted by ${\mathit{\beta }}_0$ to be ${(0.7,0.7,0.7,0,0,0,0,0,0,0)}^{\top }$ (large effects) or ${(0.4,0.4,0.4,0,0,0,0,0,0,0)}^{\top }$ (weak effects). The truncation time $\tilde{A}$ followed the uniform distribution on $(0,\tau )$ with $\tau =15$. The failure time of interest $\tilde{T}$ was generated from model (1) with $\Lambda \left(t\right)=0.3t$. Because we considered the length-biased sampling, only pairs that satisfy $\tilde{T}>\tilde{A}$ were kept in the simulated data, which were denoted as $\{(T_i,A_i,{\mathit{Z}}_i),i=1,\cdots ,n\}$.

To construct interval censoring, for subject i, we generated a series of potential examination times $U_m\sim U_{m-1}+0.1+Uniform(0,2)$ with $m=1,\cdots ,n_c$, $U_0=\tilde{A}$ and $U_{n_c}\le \tau$. Then, $[L_i,R_i)$ was defined as the smallest interval that brackets $T_i$. On average, we had about 30–44% left-censored observations and 19–30% right-censored ones. We considered some classical penalty functions, including LASSO, ALASSO, SCAD, SELO, SICA, MCP and BAR [17,18,20,21,22,24]. In order to find the optimal $\lambda$ for each penalty, we considered 20 equally spaced points over the interval $(a,b)$ and selected the one that minimizes the BIC. In particular, b was chosen to guarantee that all the regression parameter estimates were penalized to zero, while a was chosen to ensure that all the covariates were selected. The following results are based on $n=200$ or 400 and 100 replications.

To assess the variable selection performance, we calculated the FP and TP, which are defined as the average number of selected covariates whose true coefficients are zero and the average number of selected covariates whose true coefficients are not zero, respectively. To measure the estimation accuracy, we reported the median of the mean squared errors (MMSE) and the standard deviation of the mean squared errors (SD), where the mean squared error is defined as ${(\widehat{\mathit{\beta }}-\mathit{\beta })}^{\top }{\Sigma }_{\mathit{Z}}(\widehat{\mathit{\beta }}-\mathit{\beta })$ and ${\Sigma }_{\mathit{Z}}$ denotes the population covariance matrix of the covariate vector $\mathit{Z}$.

Table 1 and Table 2 present the results obtained by the proposed method with ${\mathit{\beta }}_0={(0.7,0.7,0.7,0,0,0,0,0,0,0)}^{\top }$ and ${\mathit{\beta }}_0={(0.4,0.4,0.4,0,0,0,0,0,0,0)}^{\top }$, respectively. In both tables, we also present the results of the oracle estimation and the analysis method without conducting variable selection. The results given in Table 1 and Table 2 show that almost all the penalty functions gave a similar variable selection performance. The only exception is LASSO, which yielded a slightly larger FP. This observation can be anticipated because LASSO often selects more noises than the other penalty functions [39]. In addition, one can see from the tables that except LASSO, the MMSEs yielded by the other penalty functions were close to the oracle estimation and smaller than the analysis method without conducting variable selection. As the sample size increased, the performance of the variable selection and estimation accuracy improved for all the penalty functions.

For comparison, we also include in Table 1 and Table 2 the results obtained by the variable selection method based on the penalized conditional likelihood (PCL). The detailed implementation of the PCL method is given in the Appendix B. Notably, Li et al. [32] also maximized the PCL to conduct variable selection but only considered the ALASSO penalty. It is clear that, compared with the PCL approach, the proposed method yielded smaller MMSEs, implying a more accurate estimation performance. Furthermore, the SDs obtained by the proposed method are smaller than those of the PCL method, because the PCL method ignores the distribution information of the truncation times and thus loses some estimation efficiency.

In this study, we also considered other simulation settings with $p=30$ or 50. Specifically, we set ${\beta }_j=0.7$ for $j=1,\cdots ,3$ and ${\beta }_j=0$ for $j=4,\cdots ,p$ and let other simulation specifications be the same as above. The simulation results are presented in Table 3 and Table 4 and show similar conclusions as above. In particular, the proposed method with ALASSO, SCAD, SELO, SICA, MCP and BAR yielded much smaller MMSEs than the analysis method without conducting variable selection when p increased. This clearly demonstrated the necessity of conducting variable selection in the presence of a large number of covariates.

6. An Application

6.1. Background and Analysis Methods

We applied the proposed method to a set of real data arising from the Prostate, Lung, Colorectal and Ovarian (PLCO) cancer screening trial [9,41]. Sponsored by the National Cancer Institute, the PLCO cancer screening trial was initiated in 1993 and recruited participants who had not previously taken part in any other cancer screening trials at ten screening centers nationwide. The recruited participants were aged from 55 to 74. In particular, the participants who were randomly assigned to the screening group received the Prostate-Specific Antigen (PSA) test periodically over 13 years. If abnormally high PSA levels were found, a prostate biopsy was conducted to determine the occurrence status of prostate cancer. In this study, we focused on the prostate cancer screening data in the screening group and aimed at identifying the important risk factors of the development of prostate cancer. This is because prostate cancer in general causes no signs or symptoms in the early stages, but as the disease progresses, it can cause serious complications, such as urination problems and anemia. Therefore, exploring the risk factors of prostate cancer exhibits a pressing need and is also beneficial to conduct early prevention for males. To this end, the failure time of interest was defined as the age at onset of prostate cancer. Because the participants were only examined intermittently, only interval-censored observations could be obtained for the onset of the prostate cancer. In addition, because the study excluded individuals who had already developed prostate cancer at the study recruitment, the age at the onset of prostate cancer suffered from left truncation with the truncation time being the age the individual enrolled in the study.

We considered seven potential risk factors, including Race (1 for African American and 0 otherwise), Education (1 for at least college and 0 otherwise), Cancer (1 for having an immediate family member with any PLCO cancer and 0 otherwise), ProsCancer (1 for having an immediate family member with prostate cancer and 0 otherwise), Diabetes (1 for having diabetes and 0 otherwise), Stroke (1 for having strokes before and 0 otherwise) and Gallblad (1 for having gall bladder stones and 0 otherwise). The sample size was n = 32,897, and the left and right censoring rates were about $8.67\%$ and $87.69\%$, respectively.

To achieve variable selection, we implemented the proposed method with LASSO, ALASSO, SCAD, SELO, SICA, MCP and BAR as in the simulation study. To select the optimal $\lambda$ for each penalty, we utilized a two-step method. In the first step, we examined a range of points over $(a,b)$ to roughly identify a narrower interval that includes the optimal tuning parameter. Here, as in the simulation study, a was selected to ensure that all the covariates were selected, while b was chosen to ensure that all the regression parameter estimates were penalized to zero. Next, we further considered 20 equally spaced points within the narrower interval and selected the optimal $\lambda$ that minimizes the BIC. In addition, we also employed the BIC to select the best penalty among all the penalties considered for the data, and it turned out that the PH model with SCAD and MCP yielded the smallest BIC value. To calculate the standard error, we used the nonparametric bootstrap method with 100 bootstrap samples. In addition to the proposed method, we also considered the variable method based on the penalized conditional likelihood (PCL) for comparison.

6.2. Results

We summarize in Table 5 the results obtained by the proposed and PCL methods. The results indicated that, except LASSO, the proposed method with other penalties recognized Race, ProsCancer and Diabetes as significant risk factors for prostate cancer. Specifically, being African American and having an immediate family member with prostate cancer increased the risk of developing prostate cancer, while having diabetes exhibited a lower risk of developing prostate cancer. These findings are in accordance with the conclusions obtained by Meister [42], Pierce [43] and others. In contrast, the results given in Table 5 show that the method based on PCL yielded relatively larger standard error estimates compared with the proposed method. This finding again demonstrates the advantage or efficiency gain of the proposed method when taking into account the distribution information of the truncation times in the inference procedure.

7. Discussion and Conclusions

In this article, we considered the length-biased and interval-censored data and developed a penalized analysis procedure to choose important variables among the large number of covariates in the PH model. The main contribution of this work is the development of a novel penalized EM algorithm via introducing two-stage data augmentation, which can greatly simplify the penalized nonparametric maximum likelihood estimation. Specifically, by introducing pseudo-truncated data and Poisson random variables, the possible high-dimensional parameters involved in $\Lambda \left(t\right)$ have explicit solutions, making the proposed algorithm simple and computationally stable. In contrast to the work of Li et al. [32] that only involved the ALASSO penalty, we proposed to jointly utilize the local linear approximation and the modified shooting algorithm, yielding the sparse estimators of the regression parameters under various popular penalty functions. Thus, the proposed method can offer flexible options for the data analyst. The numerical results obtained from a simulation study showed the satisfactory performance and desirable advantage of the proposed method in finite samples. Moreover, by legitimately taking into account the distribution information of the truncation times, the proposed method is more efficient than the traditional penalized conditional likelihood approach (e.g., Li et al.’s method [32]).

Notably, the findings of our prostate cancer data analysis may have certain public health implications. Specifically, African Americans and individuals who have immediate family members with prostate cancer are specific population groups and need to receive early prevention (e.g., cancer screening) in order to reduce the risk of developing prostate cancer.

8. Suggestions for Future Work

Notably, in the proposed method, we only investigated the variable selection technique in the setting of $n>p$. Obviously, in some practical applications, such as the gene expression studies, p is usually much larger than n, and future efforts will be devoted to extend the proposed method to handle the case of $p\gg n$. In addition, generalizations of the proposed method to other regression models (e.g., transformation and additive hazards models [7,44]) and taking into account multivariate interval censoring [45] warrant further research.

Author Contributions

Methodology, F.F.; Writing—review & editing, G.C.; Supervision, J.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data used in the paper are not publicly available but can be requested from https://prevention.cancer.gov/major-programs/prostate-lung-colorectal-and-ovarian-cancer-screening-trial/cancer-data-access-system.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

EM	expectation-maximization
LASSO	the least absolute shrinkage and selection operator penalty
ALASSO	the adaptive LASSO penalty
SCAD	the smoothly clipped absolute deviation penalty
SELO	the seamless-$L_0$ penalty
SICA	the smooth integration of counting and absolute deviation penalty
MCP	the minimax concave penalty
BAR	the broken adaptive ridge penalty

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.

Tables

Table 1

Simulation results for the variable selection and estimation accuracy with large effects.

Method	Penalty	$\mathit{n}=200$			$\mathit{n}=400$
		TP	FP	MMSE (SD)	TP	FP	MMSE (SD)
Proposed method	LASSO	3	1.27	0.163 (0.107)	3	1.15	0.092(0.060)
	ALASSO	3	0.12	0.051 (0.058)	3	0.11	0.025 (0.022)
	SCAD	3	0.09	0.025 (0.046)	3	0.07	0.012 (0.014)
	SELO	3	0.14	0.030 (0.042)	3	0.07	0.013 (0.016)
	SICA	3	0.16	0.030 (0.042)	3	0.07	0.013 (0.015)
	MCP	3	0.16	0.025 (0.042)	3	0.07	0.012 (0.015)
	BAR	3	0.13	0.032 (0.040)	3	0.11	0.014 (0.016)
	Oracle	-	-	0.024 (0.041)	-	-	0.011 (0.014)
	Without VS	-	-	0.057 (0.055)	-	-	0.026 (0.020)
PCL method	LASSO	3	1.42	0.163 (0.162)	3	1.29	0.089 (0.074)
	ALASSO	3	0.20	0.081 (0.108)	3	0.13	0.033 (0.041)
	SCAD	3	0.12	0.056 (0.116)	3	0.08	0.025 (0.041)
	SELO	3	0.18	0.065 (0.087)	3	0.10	0.021 (0.038)
	SICA	3	0.18	0.066 (0.087)	3	0.10	0.021 (0.037)
	MCP	3	0.13	0.054 (0.116)	3	0.10	0.025 (0.043)
	BAR	3	0.19	0.065 (0.089)	3	0.13	0.023 (0.035)

“Proposed method” denotes the proposed penalized EM algorithm; “Without VS” denotes the analysis method without conducting variable selection; and “PCL method” denotes the variable selection method based on the penalized conditional likelihood.

Table 2

Simulation results for the variable selection and estimation accuracy with weak covariate effects.

Method	Penalty	$\mathit{n}=200$			$\mathit{n}=400$
		TP	FP	MMSE (SD)	TP	FP	MMSE (SD)
Proposed method	LASSO	3	0.80	0.047 (0.040)	3	0.73	0.024(0.024)
	ALASSO	3	0.25	0.027 (0.024)	3	0.12	0.009 (0.010)
	SCAD	3	0.36	0.030 (0.028)	3	0.17	0.010 (0.011)
	SELO	3	0.23	0.017 (0.019)	3	0.10	0.008 (0.009)
	SICA	2.99	0.22	0.017 (0.021)	3	0.07	0.008 (0.008)
	MCP	2.98	0.25	0.018 (0.025)	3	0.10	0.007 (0.009)
	BAR	3	0.22	0.015 (0.018)	3	0.12	0.008 (0.009)
	Oracle	-	-	0.013 (0.014)	-	-	0.007 (0.007)
	Without VS	-	-	0.043 (0.031)	-	-	0.020 (0.013)
PCL method	LASSO	2.93	1.18	0.065 (0.075)	3	0.72	0.039 (0.040)
	ALASSO	2.85	0.50	0.060 (0.062)	2.98	0.05	0.015 (0.028)
	SCAD	2.86	0.55	0.076 (0.063)	2.99	0.15	0.017 (0.026)
	SELO	2.87	0.47	0.053 (0.060)	2.98	0.06	0.018 (0.028)
	SICA	2.88	0.48	0.053 (0.059)	2.98	0.07	0.015 (0.027)
	MCP	2.84	0.43	0.072 (0.068)	2.98	0.05	0.015 (0.028)
	BAR	2.87	0.37	0.052 (0.059)	2.97	0.08	0.014 (0.028)

“Proposed method” denotes the proposed penalized EM algorithm; “Without VS” denotes the analysis method without conducting variable selection; and “PCL method” denotes the variable selection method based on the penalized conditional likelihood.

Table 3

Simulation results for the variable selection and estimation accuracy with $p=30$.

Method	Penalty	$\mathit{n}=200$			$\mathit{n}=400$
		TP	FP	MMSE (SD)	TP	FP	MMSE (SD)
Proposed method	LASSO	3	1.59	0.307 (0.143)	3	1.55	0.184 (0.082)
	ALASSO	3	0.45	0.085 (0.075)	3	0.13	0.028 (0.030)
	SCAD	3	0.26	0.041 (0.073)	3	0.15	0.011 (0.017)
	SELO	3	0.35	0.046 (0.051)	3	0.12	0.014 (0.019)
	SICA	3	0.31	0.045 (0.051)	3	0.12	0.015 (0.019)
	MCP	3	0.14	0.030 (0.043)	3	0.13	0.012 (0.017)
	BAR	3	0.46	0.042 (0.044)	3	0.23	0.015 (0.018)
	Oracle	-	-	0.026 (0.035)	-	-	0.010 (0.017)
	Without VS	-	-	0.216 (0.148)	-	-	0.087 (0.043)
PCL method	LASSO	3	1.87	0.362 (0.223)	3	1.33	0.220 (0.129)
	ALASSO	3	0.85	0.114 (0.124)	3	0.24	0.027 (0.046)
	SCAD	2.97	0.41	0.080 (0.134)	3	0.14	0.028 (0.038)
	SELO	2.99	0.80	0.094 (0.106)	3	0.20	0.025 (0.035)
	SICA	3	0.50	0.091 (0.112)	3	0.18	0.024 (0.035)
	MCP	2.98	0.44	0.073 (0.146)	3	0.15	0.021 (0.044)
	BAR	2.99	0.77	0.100 (0.136)	3	0.37	0.026 (0.037)

“Proposed method” denotes the proposed penalized EM algorithm; “Without VS” denotes the analysis method without conducting variable selection; and “PCL method” denotes the variable selection method based on the penalized conditional likelihood.

Table 4

Simulation results for the variable selection and estimation accuracy with $p=50$.

Method	Penalty	$\mathit{n}=200$			$\mathit{n}=400$
		TP	FP	MMSE (SD)	TP	FP	MMSE (SD)
Proposed method	LASSO	3	1.54	0.394 (0.117)	3	1.33	0.269 (0.091)
	ALASSO	3	0.45	0.130 (0.076)	3	0.21	0.037 (0.033)
	SCAD	3	0.29	0.038 (0.061)	3	0.03	0.014 (0.023)
	SELO	3	0.27	0.049 (0.041)	3	0.14	0.018 (0.022)
	SICA	3	0.27	0.048 (0.047)	3	0.20	0.016 (0.020)
	MCP	3	0.14	0.028 (0.046)	3	0.06	0.014 (0.017)
	BAR	3	0.58	0.045 (0.044)	3	0.37	0.017 (0.020)
	Oracle	-	-	0.021 (0.030)	-	-	0.011 (0.016)
	Without VS	-	-	0.637 (0.453)	-	-	0.169 (0.070)
PCL method	LASSO	3	1.74	0.469 (0.263)	3	1.58	0.297 (0.139)
	ALASSO	2.99	0.82	0.155 (0.160)	3	0.38	0.037 (0.045)
	SCAD	2.99	0.61	0.101 (0.152)	3	0.20	0.024 (0.035)
	SELO	2.99	0.70	0.088 (0.105)	3	0.54	0.023 (0.051)
	SICA	2.99	0.65	0.089 (0.103)	3	0.57	0.022 (0.051)
	MCP	2.99	0.56	0.093 (0.157)	3	0.29	0.025 (0.041)
	BAR	3	1.02	0.097 (0.164)	3	0.62	0.036 (0.042)

“Proposed method” denotes the proposed penalized EM algorithm; “Without VS” denotes the analysis method without conducting variable selection; and “PCL method” denotes the variable selection method based on the penalized conditional likelihood.

Table 5

Analysis results of the prostate cancer screening data.

Method	Covariate	LASSO	ALASSO	SCAD	SELO	SICA	MCP	BAR	Without VS
Proposed	Race	0.332 (0.084 )	0.364 * (0.088)	0.444 * (0.072)	0.408 * (0.082)	0.408 * (0.083)	0.444 * (0.072)	0.412 * (0.081)	0.458 * (0.064)
method	Education	0.047 (0.027 )	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0.070 * (0.031)
	Cancer	0.032 (0.029 )	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0.047 (0.032 )
	ProsCancer	0.346 * (0.053 )	0.377 * (0.058)	0.421 * (0.049)	0.397 * (0.053)	0.398 * (0.056)	0.421 * (0.049)	0.405 * (0.052)	0.394 * (0.049)
	Diabetes	−0.310 * (0.060 )	−0.336 * (0.068)	−0.416 * (0.059)	−0.373 * (0.065)	−0.374 * (0.068)	−0.416 * (0.059)	−0.384 * (0.062)	−0.402 * (0.064)
	Stroke	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	−0.244 * (0.109 )
	Gallblad	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	−0.064 (0.058 )
PCL	Race	0.307 * (0.087 )	0.399 * (0.114)	0.420 * (0.102)	0.385 * (0.102)	0.394 * (0.105)	0.420 * (0.102)	0.377 * (0.114)	0.430 * (0.067)
method	Education	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	−0.005 (0.032 )
	Cancer	0.075 * (0.033 )	0.065 (0.049)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0.089 * (0.033 )
	ProsCancer	0.360 * (0.062 )	0.407 * (0.062)	0.453 * (0.062)	0.436 * (0.063)	0.441* (0.064)	0.453 * (0.062)	0.437 * (0.063)	0.407 * (0.051)
	Diabetes	−0.216 * (0.073 )	−0.286 * (0.078)	−0.316 * (0.077)	−0.266 * (0.089)	−0.279 * (0.093)	−0.316 * (0.077)	−0.266 * (0.087)	−0.320 * (0.065 )
	Stroke	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	−0.011 (0.109)
	Gallblad	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0 (-)	0.090 (0.060)

“Proposed method” denotes the proposed penalized EM algorithm; “PCL method” denotes the variable selection method based on the penalized conditional likelihood; “Without VS” denotes the analysis method without conducting variable selection; and “*” indicates that the covariate effect is significant at the level of 0.05.

Appendix A. Proofs of The Asymptotic Results

Appendix A.1. Proof of Theorem 1

Proof. 

The penalized log-likelihood function can be written as $$\begin{array}{c}\hfill log{L}_P(\mathit{\beta },\mathit{F})=logL(\mathit{\beta },\mathit{F})-n\lambda \sum _{j=1}^p\frac{|{\beta }_j|}{|{\tilde{\beta }}_j|},\end{array}$$ According to the definition of the profile likelihood function, for any given $\mathit{\beta }$, the log profile likelihood ${\ell }_p\left(\mathit{\beta }\right)=\underset{\mathit{F}}{sup}logL(\mathit{\beta },\mathit{F})=logL(\mathit{\beta },\mathit{F}\left(\mathit{\beta }\right))$, where $\mathit{F}\left(\mathit{\beta }\right)=arg\underset{\mathit{F}}{sup}logL(\mathit{\beta },\mathit{F})$, which can be obtained by the simplified EM algorithm. Set $logL(\tilde{\mathit{\beta }},\tilde{\mathit{F}})=\underset{\mathit{\beta },\mathit{F}}{max}logL(\mathit{\beta },\mathit{F})$ denotes the maximum value of the likelihood function. Obviously, for any $\mathit{\beta }$, we have $logL(\tilde{\mathit{\beta }},\tilde{\mathit{F}})=\underset{\mathit{\beta },\mathit{F}}{max}logL(\mathit{\beta },\mathit{F})\ge {\ell }_p\left(\mathit{\beta }\right)$. On the other hand, ${\ell }_p\left(\tilde{\mathit{\beta }}\right)=\underset{\mathit{F}}{sup}logL(\tilde{\mathit{\beta }},\mathit{F})\ge logL(\tilde{\mathit{\beta }},\tilde{\mathit{F}})$. Thus, we have ${\ell }_p\left(\tilde{\mathit{\beta }}\right)=logL(\tilde{\mathit{\beta }},\tilde{\mathit{F}})$. Therefore, the estimator $\widehat{\mathit{\beta }}$ obtained by maximizing the penalized log-likelihood is equivalent to the maximum of the following penalized log profile likelihood: (A1)$$Q_n\left(\mathit{\beta }\right)={\ell }_p\left(\mathit{\beta }\right)-n\lambda \sum _{j=1}^p\frac{|{\beta }_j|}{|{\tilde{\beta }}_j|}.$$ Under the regularity conditions (C1)–(C4), it follows from Theorem 1 of Murphy and van der Vaart [46], for any random sequence ${\mathit{\beta }}_n^{\ast }\stackrel{p}{\to }{\mathit{\beta }}_0$, the log profile likelihood ${\ell }_{pn}\left({\mathit{\beta }}_n^{\ast }\right)$ has the following quadratic expansion around ${\mathit{\beta }}_0$: (A2)$$\begin{array}{cc}\hfill {\ell }_p\left({\mathit{\beta }}_n^{\ast }\right)& ={\ell }_p\left({\mathit{\beta }}_0\right)+{\left({\mathit{\beta }}_n^{\ast }-{\mathit{\beta }}_0\right)}^{\top }\sum _{i=1}^n{\tilde{S}}_0\left({\mathcal{O}}_i\right)-\frac{1}{2}n{\left({\mathit{\beta }}_n^{\ast }-{\mathit{\beta }}_0\right)}^{\top }{\tilde{I}}_0\left({\mathit{\beta }}_n^{\ast }-{\mathit{\beta }}_0\right)\hfill \\ \hfill & +o_{P_{{\mathit{\beta }}_0,{\Lambda }_0}}{\left(\sqrt{n}∥{\mathit{\beta }}_n^{\ast }-{\mathit{\beta }}_0∥+1\right)}^2,\hfill \end{array}$$ where ${\tilde{S}}_0$ and ${\tilde{I}}_0$ are the efficient score function for $\mathit{\beta }$ and the efficient Fisher information matrix, respectively.

According to Proposition 3.1 and the discussion in Section 4.4 given by Huang and Wellner [47], ${\ell }_p$ is concave in $(\mathit{\beta },\mathit{F})$. The proof of Theorem 1 in Zeng et al. [7] implies that $\widehat{\Lambda }\left(\tau \right)$ is bounded almost surely when n is large. Hence, without sacrificing generality, we can restrict the parameter $\mathit{F}$ to a compact space. Then, ${\ell }_p\left(\mathit{\beta }\right)=\underset{\mathit{F}}{sup}logL(\mathit{\beta },\mathit{F})$ is concave in $\mathit{\beta }$. On the other hand, combined with the concavity of $-n\lambda {\sum }_{j=1}^p\frac{|{\beta }_j|}{|{\tilde{\beta }}_j|}$, we can obtain that $Q_n\left(\mathit{\beta }\right)$ is concave. In order to establish the root-n consistency of $\widehat{\mathit{\beta }}$, it is enough to prove that, for any $ϵ>0$, there exists a large constant C such that (A3)$$\underset{n\to \infty }{lim\; inf}P\left\{\underset{\parallel \mathbf{u}\parallel \le C}{sup}{Q}_n\left({\mathit{\beta }}_0+n^{-1/2}\mathbf{u}\right)<Q_n\left({\mathit{\beta }}_0\right)\right\}\ge 1-ϵ.$$ This implies that when n is sufficiently large, it can be inferred that there exists a local maximizer of $Q_n\left(\mathit{\beta }\right)$ in the interior of the ball $\{{\mathit{\beta }}_0+n^{-1/2}\mathbf{u}:\parallel \mathbf{u}\parallel \le C\}$ with a probability of at least $1-ϵ$. By the concavity of $Q_n\left(\mathit{\beta }\right)$, the local maximizer must be $\widehat{\mathit{\beta }}$ and thus $∥\widehat{\mathit{\beta }}-{\mathit{\beta }}_0∥=O_p\left(n^{-1/2}\right)$.

To prove (A3), we have $$\begin{array}{cc}\hfill \frac{1}{n}\left\{Q_n({\mathit{\beta }}_0+n^{-1/2}\mathbf{u})-Q_n\left({\mathit{\beta }}_0\right)\right\}& =\frac{1}{n}\left\{{\ell }_p({\mathit{\beta }}_0+n^{-1/2}\mathbf{u})-{\ell }_p\left({\mathit{\beta }}_0\right)\right\}\hfill \\ \hfill & -\lambda \sum _{j=1}^p\left\{\frac{|{\beta }_{j0}+n^{-1/2}{u}_j|-|{\beta }_{j0}|}{|{\tilde{\beta }}_j|}\right\}\hfill \\ \hfill & \le \frac{1}{n}\left\{{\ell }_p({\mathit{\beta }}_0+n^{-1/2}\mathbf{u})-{\ell }_p\left({\mathit{\beta }}_0\right)\right\}\hfill \\ \hfill & -\lambda \sum _{j=1}^d\left\{\frac{|{\beta }_{j0}+n^{-1/2}{u}_j|-|{\beta }_{j0}|}{|{\tilde{\beta }}_j|}\right\}.\hfill \end{array}$$ By the quadratic expansion (A2) and using the fact that $n^{-1/2}{\sum }_{i=1}^n{\tilde{S}}_0\left({\mathcal{O}}_i\right)=O_p\left(1\right)$, we have, when n is large, (A4)$$\begin{array}{cc}\hfill & \frac{1}{n}\left\{{\ell }_p({\beta }_0+n^{-1/2}\mathbf{u})-{\ell }_p\left({\beta }_0\right)\right\}\hfill \\ \hfill & =\frac{1}{n}{\mathbf{u}}^{\top }{n}^{-1/2}\sum _{i=1}^n{\tilde{S}}_0\left({\mathcal{O}}_i\right)-\frac{1}{2n}{\mathbf{u}}^{\top }{\tilde{I}}_0\mathbf{u}+o_{p_{{\beta }_0,{\Lambda }_0}}(n^{-1/2}\parallel \mathbf{u}\parallel +n^{-1/2}{)}^2\hfill \\ \hfill & =\frac{1}{n}{O}_p\left(1\right)\sum _{j=1}^p\left|u_j\right|-\frac{1}{2n}{\mathbf{u}}^{\top }{\tilde{I}}_0\mathbf{u}+o_{p_{{\beta }_0,{\Lambda }_0}}\left\{\frac{{(\parallel \mathbf{u}\parallel +1)}^2}{n}\right\}.\hfill \end{array}$$ On the other hand, we have $$\begin{array}{cc}\hfill & -\lambda \sum _{j=1}^p\left\{\frac{|{\beta }_{j0}+n^{-1/2}{u}_j|-|{\beta }_{j0}|}{|{\tilde{\beta }}_j|}\right\}\hfill \\ \hfill & \le n^{-1/2}\lambda \sum _{j=1}^d\frac{|u_j|}{|\widetilde{{\beta }_j}|}.\hfill \end{array}$$ Note that the unpenalized maximum likelihood estimator $\tilde{\mathit{\beta }}$ satisfies $\parallel \tilde{\mathit{\beta }}-{\mathit{\beta }}_0\parallel =O_p\left(n^{-1/2}\right)$. By the Taylor expansion, we have, for $1\le j\le p$, $$\frac{1}{\left|{\tilde{\beta }}_j\right|}=\frac{1}{\left|{\beta }_{j0}\right|}-\frac{sign\left({\beta }_{j0}\right)}{{\beta }_{j0}^2}\left({\tilde{\beta }}_j-{\beta }_{j0}\right)+o_p\left(\left|{\tilde{\beta }}_j-{\beta }_{j0}\right|\right)=\frac{1}{\left|{\beta }_{j0}\right|}+\frac{O_p\left(1\right)}{\sqrt{n}}.$$ Then, under the condition $\sqrt{n}\lambda =O_p\left(1\right)$, we have $$\begin{array}{cc}\hfill & n^{-1/2}\lambda \sum _{j=1}^d\frac{|u_j|}{|{\tilde{\beta }}_j|}\hfill \\ \hfill & =n^{-1/2}\lambda \sum _{j=1}^d\left(\frac{|u_j|}{|{\beta }_{j0}|}+\frac{|u_j|}{\sqrt{n}}{O}_p\left(1\right)\right)\hfill \\ \hfill & \le C{n}^{-1/2}\lambda O_p\left(1\right)\hfill \\ \hfill & =C{n}^{-1}{O}_p\left(1\right).\hfill \end{array}$$ Therefore, (A5)$$\begin{array}{cc}\hfill & \frac{1}{n}\left\{Q_n({\mathit{\beta }}_0+n^{-1/2}\mathbf{u})-Q_n\left({\mathit{\beta }}_0\right)\right\}\hfill \\ \hfill & \le \frac{1}{n}{O}_p\left(1\right)\sum _{j=1}^p\left|u_j\right|-\frac{1}{2n}{\mathbf{u}}^{\top }{\tilde{I}}_0\mathbf{u}+o_{p_{{\beta }_0,{\Lambda }_0}}\left\{\frac{{(\parallel \mathbf{u}\parallel +1)}^2}{n}\right\}+C{n}^{-1}{O}_p\left(1\right)\hfill \end{array}$$ Thus, in (A5), the first and fourth terms are of the order $C{n}^{-1}$, and the second term is of the order $C^2{n}^{-1}$ because ${\tilde{I}}_0$ is positive definite (see the proof in Zeng et al. [7]). And the third term is of an order smaller than $C^2{n}^{-1}$. Hence, if we choose a sufficiently large C, the second term in (A5) dominates the other term. As a result, (A3) holds, which completes the proof. □

Appendix A.2. Proof of Theorem 2

Proof. 

(i) According to theorem 1, we have $\parallel \widehat{\mathit{\beta }}-{\mathit{\beta }}_0\parallel =O_p\left(n^{-1/2}\right)$. To prove the sparsity property of ${\mathit{\beta }}_2$, it is sufficient to show that for any sequence ${\mathit{\beta }}_1$ such that $\parallel {\mathit{\beta }}_1-{\mathit{\beta }}_{10}\parallel =O_p\left(n^{-1/2}\right)$ and for any positive constant C, (A6)$$\underset{n\to \infty }{lim}P\left\{Q_n\left({\mathit{\beta }}_1,\mathbf{0}\right)=\underset{∥{\mathit{\beta }}_2∥\le C{n}^{-1/2}}{max}{Q}_n\left({\mathit{\beta }}_1,{\mathit{\beta }}_2\right)\right\}=1.$$ To prove (A6), we only need to show that for any such sequence ${\mathit{\beta }}_1$ and any positive constant C, $\partial Q_n\left(\mathit{\beta }\right)/\partial {\beta }_j$ have opposite signs for any ${\beta }_j\in (-C{n}^{-1/2},0)\cup (0,C{n}^{-1/2})(j=d+1,\cdots ,p)$ with the probability tending to 1.

For any ${\mathit{\beta }}_1$ satisfying $\parallel {\mathit{\beta }}_1-{\mathit{\beta }}_{10}\parallel =O_p\left(n^{-1/2}\right)$ and any ${\mathit{\beta }}_2$ such that $\parallel {\mathit{\beta }}_2\parallel \le C{n}^{-1/2}$, when n is large, ${\ell }_p\left\{\mathit{\beta }\right\}$ has the quadratic expansion (A2). It is easy to see that $\partial {\ell }_p\left\{\mathit{\beta }\right\}/\partial {\beta }_j=O_p\left(n^{-1/2}\right)(j=d+1,\cdots ,p)$. Hence, for ${\beta }_j\in (-C{n}^{-1/2},0)\cup (0,C{n}^{-1/2})(j=d+1,\cdots ,p)$, $$\begin{array}{c}\hfill \frac{\partial Q_n\left(\beta \right)}{\partial {\beta }_j}=\frac{\partial {\ell }_p\left(\beta \right)}{\partial {\beta }_j}-n\lambda \frac{sign\left({\beta }_j\right)}{\left|{\tilde{\beta }}_j\right|}=O_p\left(n^{1/2}\right)-\left(n\lambda \right)n^{1/2}\frac{sign\left({\beta }_j\right)}{\left|n^{1/2}{\tilde{\beta }}_j\right|}.\end{array}$$ Because $n^{1/2}({\tilde{\beta }}_j-0)=O_p\left(1\right)$, we have (A7)$$\frac{\partial Q_n\left(\mathit{\beta }\right)}{\partial {\beta }_j}=n^{1/2}\left\{O_p\left(1\right)-n\lambda \frac{sign\left({\beta }_j\right)}{\left|O_p\left(1\right)\right|}\right\}.$$ Because $n{\lambda }_n\to \infty$, it is easy to see that the sign of $\partial Q_n\left(\mathit{\beta }\right)/\partial {\beta }_j$ is determined by the sign of ${\beta }_j$ in (A7) when n is large, and they have opposite signs. Theorem 2(i) is proved.

(ii) We next prove the asymptotic normality of ${\widehat{\mathit{\beta }}}_1$. As shown in (i), it holds that $\underset{n\to \infty }{lim}P\left({\widehat{\mathit{\beta }}}_2=\mathbf{0}\right)=1$. Thus, we only need to derive the asymptotic representation of ${\widehat{\mathit{\beta }}}_1$ in the event $\{{\widehat{\mathit{\beta }}}_2=0\}$.

Let $\widehat{h}=\sqrt{n}({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10})$ and ${\Delta }_{1n}=n^{-1/2}{\sum }_{i=1}^n{\tilde{S}}_{10}\left({\mathcal{O}}_i\right)$, where ${\tilde{S}}_{10}$ is the sub-vector of ${\tilde{S}}_0$ corresponding to ${\mathit{\beta }}_1$. For any random sample in the event $\{{\widehat{\mathit{\beta }}}_2=0\}$, we have $Q_n\left({\widehat{\mathit{\beta }}}_1,\mathbf{0}\right)\ge Q_n\left({\mathit{\beta }}_0+n^{-1/2}\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n},\mathbf{0}\right)\right)$, where ${\tilde{I}}_{10}$ is the upper-left $d\times d$ sub-matrix of the efficient Fisher information matrix ${\tilde{I}}_0$.

Because $\parallel {\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10}\parallel =O_p\left(n^{-1/2}\right)$ according to Theorem 1 and ${\tilde{I}}_{10}^{-1}{\Delta }_{1n}=O_p\left(1\right)$. It follows by the quadratic expansion of the log profile likelihood (A2) that we have $$\begin{array}{cc}\hfill Q_n\left({\widehat{\mathit{\beta }}}_1,\mathbf{0}\right)& ={\ell }_p\left({\mathit{\beta }}_0\right)+{\widehat{h}}^T{\Delta }_{1n}-\frac{1}{2}{\widehat{h}}^T{\tilde{I}}_{10}\widehat{h}+o_{p_{{\mathit{\beta }}_0,{\Lambda }_0}}(\parallel \widehat{h}{\parallel +1)}^2-n\lambda \sum _{j=1}^d\frac{\left|{\widehat{\beta }}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \\ \hfill & ={\ell }_p\left({\mathit{\beta }}_0\right)+{\widehat{h}}^T{\Delta }_{1n}-\frac{1}{2}{\widehat{h}}^T{\tilde{I}}_{10}\widehat{h}+o_p\left(1\right)-n\lambda \sum _{j=1}^d\frac{\left|{\widehat{\beta }}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \end{array}$$ and $$\begin{array}{cc}\hfill & Q_n\left({\mathit{\beta }}_0+n^{-1/2}\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n},\mathbf{0}\right)\right)\hfill \\ \hfill & ={\ell }_p\left({\mathit{\beta }}_0\right)+{\Delta }_{1n}^T{\tilde{I}}_{10}^{-1}{\Delta }_{1n}-\frac{1}{2}{\Delta }_{1n}^T{\tilde{I}}_{10}^{-1}{\Delta }_{1n}+o_{p_{{\beta }_0,{\Lambda }_0}}{\left(∥{\tilde{I}}_{10}^{-1}{\Delta }_{1n}∥+1\right)}^2\hfill \\ \hfill & -n\lambda \sum _{j=1}^d\frac{\left|{\beta }_{j0}+n^{-1/2}{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \\ \hfill & ={\ell }_p\left({\mathit{\beta }}_0\right)+\frac{1}{2}{\Delta }_{1n}^T{\tilde{I}}_{10}^{-1}{\Delta }_{1n}+o_p\left(1\right)-n\lambda \sum _{j=1}^d\frac{\left|{\beta }_{j0}+n^{-1/2}{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}.\hfill \end{array}$$ Thus, (A8)$$\begin{array}{cc}\hfill & Q_n\left({\widehat{\mathit{\beta }}}_1,\mathbf{0}\right)-Q_n\left({\mathit{\beta }}_0+n^{-1/2}\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n},\mathbf{0}\right)\right)\hfill \\ \hfill & ={\widehat{h}}^T{\Delta }_{1n}-\frac{1}{2}{\widehat{h}}^T{\tilde{I}}_{10}\widehat{h}-\frac{1}{2}{\Delta }_{1n}^T{\tilde{I}}_{10}^{-1}{\Delta }_{1n}-n\lambda \sum _{j=1}^d\frac{\left|{\widehat{\beta }}_j\right|-\left|{\beta }_{j0}+n^{-1/2}{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \\ \hfill & =-\frac{1}{2}{\left(\widehat{h}-{\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}^{\top }{\tilde{I}}_{10}\left(\widehat{h}-{\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)-n\lambda \sum _{j=1}^d\frac{\left|{\widehat{\beta }}_j\right|-\left|{\beta }_{j0}+n^{-1/2}{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \\ \hfill & \le -c{∥\widehat{h}-{\tilde{I}}_{10}^{-1}{\Delta }_{1n}∥}^2+\sqrt{n}\lambda \sum _{j=1}^d\frac{\left|\sqrt{n}\left({\widehat{\beta }}_j-{\beta }_{j0}\right)-{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}\hfill \end{array}$$ for some positive constant c because ${\tilde{I}}_{10}$ is positive-definite.

On the other hand, because $Q_n\left({\widehat{\mathit{\beta }}}_1,\mathbf{0}\right)-Q_n\left({\mathit{\beta }}_0+n^{-1/2}\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n},\mathbf{0}\right)\right)\ge 0$, combining (A8) we have that $${∥\widehat{h}-{\tilde{I}}_{10}^{-1}{\Delta }_{1n}∥}^2\le \frac{1}{c}\sqrt{n}\lambda \sum _{j=1}^d\frac{\left|\sqrt{n}\left({\widehat{\beta }}_j-{\beta }_{j0}\right)-{\left({\tilde{I}}_{10}^{-1}{\Delta }_{1n}\right)}_j\right|}{\left|{\tilde{\beta }}_j\right|}+o_p\left(1\right).$$ Note that $\left|{\tilde{\beta }}_j\right|\stackrel{p}{\to }\left|{\beta }_{j0}\right|>0$ for $j=1,\cdots ,d$. Thus, under the condition $\sqrt{n}\lambda \to 0$, $\parallel \widehat{h}-{\tilde{I}}_{10}^{-1}{\Delta }_{1n}\parallel =o_p\left(1\right)$, and then $\sqrt{n}({\widehat{\mathit{\beta }}}_1-{\mathit{\beta }}_{10})={\tilde{I}}_{10}^{-1}{\Delta }_{1n}+o_p\left(1\right)⇝N\left(\mathbf{0},{\tilde{I}}_{10}^{-1}\right).$ □

Appendix B. The Penalized Conditional Likelihood Method

In this section, we provide the EM algorithm for implementing the maximum conditional likelihood approach. The conditional likelihood function can be written as $$L_C(\mathit{\beta },\Lambda )=\prod _{i=1}^n\frac{exp\{-\Lambda \left(L_i\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\}-I(R_i<\infty )exp\{-\Lambda \left(R_i\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\}}{exp\{-\Lambda \left(A_i\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\}}.$$ To maximize the above conditional likelihood function, we treat $\Lambda$ as a step function. Specifically, let $0<t_1<\cdots <t_{K_n}$ denote the ordered distinct time point of all $L_i,R_{i}I(R_i<\infty )$, and $A_i>0$. Let $f_k$ denote the non-negative jump size of the estimator for $\Lambda$ at $t_k$ for $k=1,\cdots ,K_n$. Then, we have $\Lambda \left(t\right)={\displaystyle \sum _{t_k\le t}}{f}_k$ and the above conditional likelihood function can be re-expressed as (A9)$$L_{C1}(\mathit{\beta },\mathit{F})=\prod _{i=1}^{n}exp\left\{-\sum _{A_i\le t_k\le L_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}\times {\left[1-exp\left\{-\sum _{L_i<t_k\le R_i}{f}_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)\right\}\right]}^{I(R_i<\infty )}$$ where $\mathit{F}={(f_1,f_2,\cdots ,f_{K_n})}^{\top }$.

To present the EM algorithm, for the ith subject, we introduce a set of new independent latent variables $\{W_{ik},k=1,\cdots ,K_n,A_i\le t_k\le R_i^{\ast }\}$, where $R_i^{\ast }=L_{i}I(R_i=\infty )+R_{i}I(R_i<\infty )$ and $W_{ik}$ is a Poisson latent variable with the mean $f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right)$. Then, the likelihood function (A9) can be equivalently expressed as $$L_{C2}(\mathit{\beta },\mathit{F})=\prod _{i=1}^{n}P\left(\sum _{A_i\le t_k\le L_i}{W}_{ik}=0\right)P{\left(\sum _{L_i<t_k\le R_i}{W}_{ik}>0\right)}^{I(R_i<\infty )}.$$ Let $P\left(W_{ik}\right|f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right))$ be the probability mass function of $W_{ik}$, and if we assume that the latent variables $W_{ik}$’s were known, we would have the following complete data likelihood function $$L_C^c(\mathit{\beta },\mathit{F})=\prod _{i=1}^n\prod _{k=1}^{K_n}P{\left(W_{ik}\right|f_{k}exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_i\right))}^{I(A_i\le t_k\le R_i^{\ast })},$$ where ${\displaystyle \sum _{t_k\le L_i}}{W}_{ik}=0$, ${\displaystyle \sum _{L_i<t_k\le R_i}}{W}_{ik}>0$ if $R_i<\infty$; ${\displaystyle \sum _{t_k\le L_i}}{W}_{ik}=0$, if $R_i=\infty$.

In the E-step of the EM algorithm, we need to determine the conditional expectation of the log-likelihood function $log{L}_C^c(\mathit{\beta },\mathit{F})$ with respect to all the latent variables and obtain the following objective function: (A10)$$Q_C(\mathit{\theta };{\mathit{\theta }}^{\left(m\right)})=\sum _{i=1}^n\sum _{k=1}^{K_n}I(A_i\le t_k\le R_i^{\ast })\left\{E\left(W_{ik}\right)log{f}_k+E\left(W_{ik}\right){\mathit{\beta }}^{\top }{Z}_i-f_{k}exp\left({\beta }^T{Z}_i\right)\right\}$$ We need to calculate the following conditional expectations at the mth iteration of the algorithm. $$E\left(W_{ik}\right)=I(L_i<t_k\le R_i,R_i<\infty )\times \frac{f_k^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)}{1-exp\left\{-\sum _{L_i<t_l\le R_i}{f}_l^{\left(m\right)}exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)\right\}}.$$ In the M-step, we can update each $f_k$ by solving $\partial Q(\mathit{\theta };{\mathit{\theta }}^{\left(m\right)})/\partial f_k=0$ as follows: (A11)$$f_k=\frac{{\displaystyle \sum _{i=1}^n}I(A_i\le t_k\le R_i^{\ast })E\left(W_{ik}\right)}{{\displaystyle \sum _{i=1}^n}I(A_i\le t_k\le R_i^{\ast })exp\left({\mathit{\beta }}^{\left(m\right)\top }{\mathit{Z}}_i\right)}.$$ Plugging (A11) into (A10), we have the following objective function $$H(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})+\sum _{j=1}^p{P}_{\lambda }\left(\right|{\beta }_j\left|\right)$$ where $$\begin{array}{cc}\hfill H(\mathit{\beta };{\mathit{\theta }}^{\left(m\right)})& =-\frac{1}{n}\sum _{i=1}^n\sum _{k=1}^{K_n}I(A_i\le t_k\le R_i^{\ast })E\left(W_{ik}\right)\{{\mathit{\beta }}^{\top }{\mathit{Z}}_i\hfill \\ \hfill & -log\left[\sum _{i^{{}^{\prime }}=1}^{n}I\left(A_{i^{{}^{\prime }}}\le t_k\le R_{i^{{}^{\prime }}}^{\ast }\right)exp\left({\mathit{\beta }}^{\top }{\mathit{Z}}_{i^{{}^{\prime }}}\right)\right]\}.\hfill \end{array}$$ The sparse estimator of $\mathit{\beta }$ based on $L_C(\mathit{\beta },\mathit{F})$ can be obtained by minimizing the above objective function, which can be accomplished with the same algorithm as in Section 3 of the paper.

References

1. Sun, J. The Statistical Analysis of Interval-Censored Failure Time Data; Springer: New York, NY, USA, 2006.

2. Huang, J. Efficient estimation for the proportional hazards model with interval censoring. Ann. Stat.; 1996; 24, pp. 540-568. [DOI: https://dx.doi.org/10.1214/aos/1032894452]

3. Shen, X. Proportional odds regression and sieve maximum likelihood estimation. Biometrika; 1998; 85, pp. 165-177. [DOI: https://dx.doi.org/10.1093/biomet/85.1.165]

4. Zeng, D.; Cai, J.; Shen, Y. Semiparametric additive risks model for interval-censored data. Stat. Sin.; 2006; 16, pp. 287-302.

5. Zhang, Y.; Hua, L.; Huang, J. A spline-based semiparametric maximum likelihood estimation method for the Cox model with interval-censored data. Scand. J. Stat.; 2010; 37, pp. 338-354. [DOI: https://dx.doi.org/10.1111/j.1467-9469.2009.00680.x]

6. Wang, L.; McMahan, C.S.; Hudgens, M.G.; Qureshi, Z.P. A flexible, computationally efficient method for fitting the proportional hazards model to interval-censored data. Biometrics; 2016; 72, pp. 222-231. [DOI: https://dx.doi.org/10.1111/biom.12389] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26393917]

7. Zeng, D.; Mao, L.; Lin, D.Y. Maximum likelihood estimation for semiparametric transformation models with interval-censored data. Biometrika; 2016; 103, pp. 253-271. [DOI: https://dx.doi.org/10.1093/biomet/asw013]

8. Zhou, Q.; Hu, T.; Sun, J. A sieve semiparametric maximum likelihood approach for regression analysis of bivariate interval-censored failure time data. J. Am. Stat. Assoc.; 2017; 112, pp. 664-672. [DOI: https://dx.doi.org/10.1080/01621459.2016.1158113]

9. Prorok, P.C.; Andriole, G.L.; Bresalier, R.S.; Buys, S.S.; Chia, D.; Crawford, E.D.; Fogel, R.; Gelmann, E.P.; Gilbert, F.; Gohagan, J.K. Design of the prostate, lung, colorectal and ovarian (PLCO) cancer screening trial. Control. Clin. Trials; 2000; 21, pp. 273S-309S. [DOI: https://dx.doi.org/10.1016/S0197-2456(00)00098-2]

10. Wang, M.C. Nonparametric estimation from cross-sectional survival data. J. Am. Stat. Assoc.; 1991; 86, pp. 130-143. [DOI: https://dx.doi.org/10.1080/01621459.1991.10475011]

11. Shen, Y.; Ning, J.; Qin, J. Analyzing length-biased data with semiparametric transformation and accelerated failure time models. J. Am. Stat. Assoc.; 2009; 104, pp. 1192-1202. [DOI: https://dx.doi.org/10.1198/jasa.2009.tm08614]

12. Ning, J.; Qin, J.; Shen, Y. Semiparametric accelerated failure time model for length-biased data with application to dementia study. Stat. Sin.; 2014; 24, pp. 313-333. [DOI: https://dx.doi.org/10.5705/ss.2011.197] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24478570]

13. Qin, J.; Shen, Y. Statistical methods for analyzing right-censored length-biased data under Cox model. Biometrics; 2010; 66, pp. 382-392. [DOI: https://dx.doi.org/10.1111/j.1541-0420.2009.01287.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19522872]

14. Qin, J.; Ning, J.; Liu, H.; Shen, Y. Maximum likelihood estimations and EM algorithms with length-biased data. J. Am. Stat. Assoc.; 2011; 106, pp. 1434-1449. [DOI: https://dx.doi.org/10.1198/jasa.2011.tm10156] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22323840]

15. Gao, F.; Chan, K.C.G. Semiparametric regression analysis of length-biased interval-censored data. Biometrics; 2019; 75, pp. 121-132. [DOI: https://dx.doi.org/10.1111/biom.12970]

16. Shen, P.S.; Peng, Y.; Chen, H.J.; Chen, C.M. Maximum likelihood estimation for length-biased and interval-censored data with a nonsusceptible fraction. Lifetime Data Anal.; 2022; 28, pp. 68-88. [DOI: https://dx.doi.org/10.1007/s10985-021-09536-2]

17. Tibshirani, R. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B (Methodol.); 1996; 58, pp. 267-288. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1996.tb02080.x]

18. Fan, J.; Li, R. 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]

19. 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]

20. Lv, J.; Fan, Y. A unified approach to model selection and sparse recovery using regularized least squares. Ann. Stat.; 2009; 37, pp. 3498-3528. [DOI: https://dx.doi.org/10.1214/09-AOS683]

21. Dicker, L.; Huang, B.; Lin, X. Variable selection and estimation with the seamless-L-0 penalty. Stat. Sin.; 2013; 23, pp. 929-962.

22. Zhang, C.H. Nearly unbiased variable selection under minimax concave penalty. Ann. Stat.; 2010; 38, pp. 894-942. [DOI: https://dx.doi.org/10.1214/09-AOS729]

23. Liu, Z.; Li, G. Efficient regularized regression with penalty for variable selection and network construction. Comput. Math. Methods Med.; 2016; 2016, 3456153. [DOI: https://dx.doi.org/10.1155/2016/3456153]

24. Dai, L.; Chen, K.; Sun, Z.; Liu, Z.; Li, G. Broken adaptive ridge regression and its asymptotic properties. J. Multivar. Anal.; 2018; 168, pp. 334-351. [DOI: https://dx.doi.org/10.1016/j.jmva.2018.08.007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30911202]

25. Fan, J.; Li, R.; Zhang, C.H.; Zou, H. Statistical Foundations of Data Science; Chapman and Hall/CRC: New York, NY, USA, 2020.

26. Garavand, A.; Salehnasab, C.; Behmanesh, A.; Aslani, N.; Zadeh, A.; Ghaderzadeh, M. Efficient Model for Coronary Artery Disease Diagnosis: A Comparative Study of Several Machine Learning Algorithms. J. Healthc. Eng.; 2022; 2022, 5359540. [DOI: https://dx.doi.org/10.1155/2022/5359540] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36304749]

27. Hosseini, A.; Eshraghi, M.A.; Taami, T.; Sadeghsalehi, H.; Hoseinzadeh, Z.; Ghaderzadeh, M.; Rafiee, M. A mobile application based on efficient lightweight CNN model for classification of B-ALL cancer from non-cancerous cells: A design and implementation study. Inform. Med. Unlocked; 2023; 39, 101244. [DOI: https://dx.doi.org/10.1016/j.imu.2023.101244]

28. Garavand, A.; Behmanesh, A.; Aslani, N.; Sadeghsalehi, H.; Ghaderzadeh, M. Towards Diagnostic Aided Systems in Coronary Artery Disease Detection: A Comprehensive Multiview Survey of the State of the Art. Int. J. Intell. Syst.; 2023; 2023, 6442756. [DOI: https://dx.doi.org/10.1155/2023/6442756]

29. Ghaderzadeh, M.; Aria, M. Management of Covid-19 Detection Using Artificial Intelligence in 2020 Pandemic. Proceedings of the ICMHI ’21: 5th International Conference on Medical and Health Informatics; Kyoto, Japan, 14–16 May 2021; pp. 32-38.

30. Chen, L.P. Variable selection and estimation for the additive hazards model subject to left-truncation, right-censoring and measurement error in covariates. J. Stat. Comput. Simul.; 2020; 90, pp. 3261-3300. [DOI: https://dx.doi.org/10.1080/00949655.2020.1800705]

31. He, D.; Zhou, Y.; Zou, H. High-dimensional variable selection with right-censored length-biased data. Stat. Sin.; 2020; 30, pp. 193-215. [DOI: https://dx.doi.org/10.5705/ss.202018.0089]

32. Li, C.; Pak, D.; Todem, D. Adaptive lasso for the Cox regression with interval censored and possibly left truncated data. Stat. Methods Med. Res.; 2020; 29, pp. 1243-1255. [DOI: https://dx.doi.org/10.1177/0962280219856238]

33. Withana Gamage, P.; McMahan, C.; Wang, L. Variable selection in semiparametric nonmixture cure model with interval-censored failure time data. Stat. Med.; 2019; 38, pp. 3026-3039.

34. Li, S.; Peng, L. Instrumental Variable Estimation of Complier Causal Treatment Effect with Interval-Censored Data. Biometrics; 2023; 79, pp. 253-263. [DOI: https://dx.doi.org/10.1111/biom.13565]

35. Withana Gamage, P.; McMahan, C.; Wang, L. A flexible parametric approach for analyzing arbitrarily censored data that are potentially subject to left truncation under the proportional hazards model. Lifetime Data Anal.; 2023; 29, pp. 188-212. [DOI: https://dx.doi.org/10.1007/s10985-022-09579-z]

36. Huang, C.Y.; Qin, J. Semiparametric estimation for the additive hazards model with left-truncated and right-censored data. Biometrika; 2013; 100, pp. 877-888. [DOI: https://dx.doi.org/10.1093/biomet/ast039]

37. Turnbull, B.W. The empirical distribution function with arbitrarily grouped, censored and truncated data. J. R. Stat. Soc. Ser. (Methodol.); 1976; 38, pp. 290-295. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1976.tb01597.x]

38. 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]

39. Li, S.; Wu, Q.; Sun, J. Penalized estimation of semiparametric transformation models with interval-censored data and application to Alzheimer’s disease. Stat. Methods Med. Res.; 2020; 29, pp. 2151-2166. [DOI: https://dx.doi.org/10.1177/0962280219884720] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31718478]

40. Zou, H.; Li, R. One-step sparse estimates in nonconcave penalized likelihood models. Ann. Stat.; 2008; 36, pp. 1509-1533.

41. Andriole, G.L.; Crawford, E.D.; Grubb, R.L.; Buys, S.S.; Chia, D.; Church, T.R.; Fouad, M.N.; Isaacs, C.; Prorok, P. Prostate cancer screening in the randomized prostate, lung, colorectal, and ovarian cancer screening trial: Mortality results after 13 years of follow-up. J. Natl. Cancer Inst.; 2012; 104, pp. 125-132. [DOI: https://dx.doi.org/10.1093/jnci/djr500] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22228146]

42. Meister, K. Risk Factors for Prostate Cancer; American Council on Science and Health: New York, NY, USA, 2002.

43. Pierce, B.L. Why are diabetics at reduced risk for prostate cancer? A review of the epidemiologic evidence. Urol. Oncol. Semin. Orig. Investig.; 2012; 30, pp. 735-743. [DOI: https://dx.doi.org/10.1016/j.urolonc.2012.07.008]

44. Lu, T.; Li, S.; Sun, L. Combined estimating equation approaches for the additive hazards model with left-truncated and interval-censored data. Lifetime Data Anal.; 2023; 29, pp. 672-697. [DOI: https://dx.doi.org/10.1007/s10985-023-09596-6]

45. 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] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34407218]

46. Murphy, S.A.; Van Der Vaart, A.W. On profile likelihood. J. Am. Stat. Assoc.; 2000; 95, pp. 449-465. [DOI: https://dx.doi.org/10.1080/01621459.2000.10474219]

47. Huang, J.; Wellner, J.A. Interval censored survival data: A review of recent progress. Proceedings of the First Seattle Symposium in Biostatistics; Lin, D.Y.; Fleming, T.R. Springer: New York, NY, USA, 1997; pp. 123-169.