Concave Group Selection of Nonparameter Additive Accelerated Failure Time Model

Ling Zhu

Jinan University, Guangzhou, China.

DOI: 10.4236/ojs.2021.111008                

Abstract

In this paper, we have studied the nonparameter accelerated failure time (AFT) additive regression model, whose covariates have a nonparametric effect on high-dimensional censored data. We give the asymptotic property of the penalty estimator based on GMCP in the nonparameter AFT model.

Keywords

Accelerated Failure Time Model, Nonparameter Model, Group Minimax Concave Penalty, Weighted Least Squares Estimation

1. Introduction

With the development of the Internet, high-dimensional data has been widely collected in life, especially in the field of medical research and finance, the results or responses of data are censored, so the study of high-dimensional censored data is meaningful. However, due to the impact of “disaster of dimension”, the study of high-dimensional data becomes extremely difficult, and some special methods must be adopted to deal with it. As the number of data dimensions increases, the performance of high-dimensional data structures declines rapidly. In low-dimensional spaces, we often use Euclidean distance to measure the similarity between data; but in high-dimensional spaces, this kind of similarity no longer exists, which makes the data mining of high-dimensional data very severely challenging. On the one hand, the performance of the data mining algorithm based on the index structure is reduced; on the other hand, many mining methods based on the entire spatial distance function will fail. By reducing the number of dimensions, the data can be reduced from high to low dimensions, and then using low-dimensional data processing methods. Therefore, the study of effective dimensionality reduction methods becomes significant in statistics.

In many studies, the main results or responses of survival data are censored. Survival analysis is another important theme of statistics, and it has been widely used in medical research and finance. Therefore, the study of survival data has attracted a lot of attention. The Cox model [1] is the most commonly used regression model for survival data. The alternative method of the PH model is the accelerated failure time model, which directly correlates the logarithm of the failure time with the covariate, and is similar to the traditional linear model, which is easier to explain than the PH model. [2] takes into account both Lasso and threshold gradient oriented regularization for high-dimensional AFT model estimation and variable selection. [3] uses partial least squares (PLS) and Lasso methods to select variables in AFT models with high-dimensional covariates. [4] proposed a robust weighted minimum absolute deviation method to estimate the high-dimensional AFT model. [5] uses COSSO penalty in the nonparameter AFT model for variable selection [6] in the high-dimensional nonparameter AFT model, using the reproduction kernel Hilbert norm penalty for estimation, a new enhanced algorithm is proposed for censoring time data. The algorithm is suitable for fitting parameter accelerated failure time models. [7] studied the elastic net method for variable selection under the Cox proportional hazard model and the AFT model with high-dimensional covariates. [8] developed a robust prediction model for event time results through LASSO regularization. This model is aimed at the Gehan estimation of high-dimensional prediction variables accelerated failure time AFT model. [9] extends rank-based Lasso estimation to the estimation and variable selection in the high-dimensional partial linear acceleration failure time model. [10] uses the bridge penalty for regular estimation and parameter selection of high-dimensional AFT models. Based on the high-dimensional semi-parameter accelerated failure time model, [11] proposed the Buckley-James method of double penalty, which can perform variable selection and parameter estimation at the same time. [12] has developed a method for quickly predicting variable selection and contraction estimation of high-dimensional predictive variable AFT models. The model is related to the correlation vector machine (RVM), which relies on maximum posterior estimation to get sparse estimates quickly. [13] proposes a semiparametric regression model whose covariate effect contains parametric and nonparametric parts. The selection of parametric covariates is achieved by iterative LASSO method, and the nonparametric components are estimated using the sieve method [14], and based on kullback-leibler geometry [15], an empirical model selection tool for nonparameter components was obtained. However, they leave behind some theoretical issues that have not yet been resolved. [16] takes into account the estimation and variable selection of LASSO and MCP in AFT models with high covariates. [17] implements regularization in the high-dimensional AFT model L_1/2 for variable selection. [18] proposed a covariate adjustment screening and variable selection procedure under the accelerated failure time model. It also appropriately adjusted the low-dimensional confounding factors to achieve a more accurate estimation of regression coefficients. [19] proposed an adaptive elastic net and weighted elastic net with censored data and high-dimensional variable selection in the AFT model. [20] proposed to apply a tensor recursive neural network architecture to extract latent representations from the entire patient medical record of the high-dimensional AFT model. [21] considers a novel Sparse L₂ Boosting algorithm, which is based on a semiparameter variable coefficient accelerated failure time model of right-censored survival data with high-dimensional covariates model prediction and variable selection. [22] developed a variable selection method in an AFT model with high-dimensional predictive variables, which consists of a set of algorithms based on two widely used techniques in the field of variable selection in survival analysis synthesis: Buckley-James method and Dantzig selector.

In this article, based on potential predictors, we applied the GMCP (Group Minimax Concave Penalty) penalty method for the first time to the study of a high-dimensional nonparametric accelerated failure time additive regression model (2.1) (MCP, [23] ). The weighted least squares solution of the model based on GMCP penalty is given. We also derived the group coordinate descent algorithm used to calculate the GMCP estimate in this model. Our simulation results show that the weighted least squares estimation based on GMCP penalty works well in the high-dimensional nonparameter accelerated failure time additive regression model, and is superior to the GLasso (Group Least Absolute Shrinkage and Selection Operator) penalty method.

The rest of the paper is organized as follows. In Section 2, we describe the nonparameter accelerated failure time additive regression (NP-AFT-AR) model and our research methods. In Section 3, we give the asymptotic oracle property of GMCP estimation. The simulation results are given in Section 4. Verification of actual data is given in Section 5. The conclusion is given in Section 6.

2. Models and Methods

2.1. Model

In this paper, we study the following nonparametric accelerated failure time additive regression (NP-AFT-AR) model to describe the relationship between the independent predictors or covariates X_j’s and the failure time T:

$T=\exp \left({\eta }_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{f}_j\left(X_j\right)+\epsilon \right)$ (2.1)

where ${\eta }_0$ is the intercept, $X=\left(X_1\mathrm{,}\cdots \mathrm{,}{X}_p\right)$ is a $p\times 1$ vector of covariates, f_j’s are unknown smooth functions with zero means, i.e., $E{f}_j\left(X_j\right)=0$ and $\epsilon$ is the random error term with mean zero and a finite variance ${\sigma }^2$. We consider sample size is small $n<p$, assuming that some additive components $f_j$ are zero, the main purpose of our research is to find the non-zero components and zero components; the second goal is to find the specific functional form of the non-zero components in order to propose a more parsimonious model. In this study, we apply the GMCP penalty in the proposed NP-AFT-AR model for component selection and estimation. We use B-splines to parameterize the nonparameter components, then invoke the inverse probability-of-censoring weighted least squares method to achieve the goals. We treat the spline approximation for each component as a group of variables subject to selection. By the GMCP penalty approach, we show that the proposed method can select significant component functions by choosing the nonzero spline basis functions.

2.2. Weighted Least Squares Estimation

We define $T_i$ as the i^th subject’s survival time, and let $C_i$ denote the censoring time and ${\delta }_i$ denote the event indicator, i.e., ${\delta }_i=I\left(T_i\le C_i\right)$ ; which takes value 1 if the event time is observed, or 0 if the event time is censored. Define $Y_i$ as the minimum of the survival time and the censoring time, i.e., $Y_i=\log \left(\min \left(T_i\mathrm{,}{C}_i\right)\right)$ : Then, the observed data are in the form $\left(Y_i\mathrm{,}{\delta }_i\mathrm{,}{X}_i\right)$, $i=\mathrm{1,}\cdots \mathrm{,}n$. which are assumed to be an independent and identically distributed (i.i.d.) sample from $\left(Y\mathrm{,}\delta \mathrm{,}X\right)$.

Let $Y_{\left(1\right)}\le \cdots \le Y_{\left(n\right)}$ be the order statistics of Y_i’s, ${\delta }_{\left(1\right)}\mathrm{,}\cdots ,{\delta }_{\left(n\right)}$ and $X_{\left(1\right)}\mathrm{,}\cdots ,X_{\left(n\right)}$ are the associated censoring indicators and covariates. Let F be the distribution of T and $\widehat{F_n}$ be its Kaplan-Meier estimator $\widehat{F_n}\left(y\right)={\displaystyle {\sum }_{i=1}^n}{\omega }_{ni}1\left(Y_{\left(i\right)}\le y\right)$, where the ${\omega }_{ni}$ ’s are Kaplan-Meier weights ( [24] ) calculated by

${\omega }_{n1}=\frac{{\delta }_{\left(1\right)}}{n},{\omega }_{ni}=\frac{{\delta }_{\left(i\right)}}{n-i+1}{\displaystyle \underset{j=1}{\overset{i-1}{\prod }}}{\left(\frac{n-j}{n-j+1}\right)}^{{\delta }_{\left(j\right)}},i=2,\cdots ,n$

[4] showed that the weights, ${\omega }_{ni}$ ’s, are the jumps in the Kaplan-Meier estimator. These are equivalent to the inverse probability-of-censoring weights ( [25] [26] ), ${\omega }_{ni}={\delta }_{\left(i\right)}/\widehat{G_n}\left(Y_{\left(i\right)}-\right)$ ; where $\widehat{G_n}$ is the Kaplan-Meier estimator of G, the distribution function of C. The Stute’s weighted least squares loss function for the NP-AFT-AR model (2.1) is defined as

$Q_n=\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}n{\omega }_{ni}{\left\{Y_{\left(i\right)}-{\eta }_0-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{f}_j\left(X_{\left(i\right)j}\right)\right\}}^2$ (2.2)

Here, we use B-spline basis functions to approximated unknown functions f_j’s. For every function component, assuming that $X_j$ is bounded; and $E\left\{f_j\left(X_j\right)\right\}=0,j=1,\cdots ,p$ ; The basis functions are determined by the order $\left(p+1\right)$ and the number of interior knots $\kappa$. The total number of B-spline basis functions for each function component would be $p+\kappa +1$ : For identifiability, satisfy $E{f}_j\left(X_j\right)=0$ ; we take the total number of basis functions to be $M_n=p+\kappa$ only and center all the basis functions at their means. Then the B-splines approximation for each function component, $f_j\left(X_j\right)\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}p$ ; is given by

$f_j\left(X_j\right)\approx {\displaystyle \underset{k\mathrm{<mo>\; =\; </mo>1}}{\overset{M_n}{\sum }}}{\beta }_{jk}{B}_{jk}(\; X\; j\; )$

where $B_{jk}\left(X_j\right)$ are the B-spline basis functions and ${\beta }_j={\left({\beta }_{j1}\mathrm{,}\cdots \mathrm{,}{\beta }_{j{M}_n}\right)}^{\text{T}}$ is the corresponding coefficient parameter vector. Let $B_j$ denote the $n\times M_n$ design matrix of B-spline basis of the j^th predictor and $B_{j\left(i\right)}$ be its i^th row vector corresponding to the sorted data. Denote the $n\times p{M}_n$ design matrix as $B=\left(B_1\mathrm{,}{B}_2\mathrm{,}\cdots \mathrm{,}{B}_p\right)$ ; the i^th row of $B$ as $B_{\left(i\right)}$ ; and the corresponding parameter vector as $\beta ={\left({\beta }_1^{\text{T}}\mathrm{,}\cdots \mathrm{,}{\beta }_p^{\text{T}}\right)}^{\text{T}}$. Then we have

$\begin{array}{c}{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{f}_j\left(X_{\left(i\right)j}\right)=f_1\left(X_{\left(i\right)1}\right)+\cdots +f_p\left(X_{\left(i\right)p}\right)\\ \approx {\displaystyle \underset{k=1}{\overset{M_n}{\sum }}}{\beta }_{1k}{B}_{1k}\left(X_{\left(i\right)1}\right)+\cdots +{\displaystyle \underset{k=1}{\overset{M_n}{\sum }}}{\beta }_{pk}{B}_{pk}\left(X_{\left(i\right)p}\right)\\ ={\displaystyle \underset{j=1}{\overset{p}{\sum }}}{B}_{\left(i\right)j}{\beta }_j\end{array}$ (2.3)

By plugging Equation (2.3) into Equation (2.2), we will get the new loss function as following:

$Q_n\left({\eta }_0\mathrm{,}\beta \right)=\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}n{\omega }_{ni}{\left\{Y_{\left(i\right)}-{\eta }_0-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{B}_{\left(i\right)j}{\beta }_j\right\}}^2$ (2.4)

By centering $B_{\left(i\right)j}$ and $Y_{\left(i\right)}$ with their ${\omega }_{ni}$ -weighted means, the intercept becomes 0. Denote ${\tilde{B}}_{\left(i\right)j}={\left(n{\omega }_{ni}\right)}^{1/2}\left({\hat{B}}_{\left(i\right)j}-{\bar{B}}_{j\omega }\right)$ and $\tilde{Y}={\left(n{\omega }_{ni}\right)}^{1/2}\left(Y_{\left(i\right)}-{\bar{Y}}_{\omega }\right)$ ; where ${\bar{B}}_{j\omega }={\displaystyle {\sum }_{i=1}^n}{\omega }_{ni}{B}_{\left(i\right)j}/{\displaystyle {\sum }_{i=1}^n}{\omega }_{ni}$ and ${\bar{Y}}_{\omega }={\displaystyle {\sum }_{i=1}^n}{\omega }_{ni}{Y}_{\left(i\right)}/{\displaystyle {\sum }_{i=1}^n}{\omega }_{ni}$ Let ${\Vert a\Vert }_2={\left({\displaystyle {\sum }_{j=1}^m}{\left|a_j\right|}^2\right)}^{1/2}$ denote the L₂ norm of any vector $a\in R^m$. For simplicity, we use ${\tilde{B}}_j={\left({\tilde{B}}_{\left(1\right)j}\mathrm{,}\cdots \mathrm{,}{\tilde{B}}_{\left(n\right)j}\right)}^{\text{T}}$ and $\tilde{Y}={\left({\tilde{Y}}_{\left(1\right)}\mathrm{,}\cdots \mathrm{,}{\tilde{Y}}_{\left(n\right)}\right)}^{\text{T}}$. Then we can rewrite the Stute’s weighted least squares loss function Equation (2.4) as

$Q_n\left(\beta \right)=\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left\{{\tilde{Y}}_{\left(i\right)}-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\tilde{B}}_{\left(i\right)j}{\beta }_j\right\}}^2=\frac{1}{2}{\Vert \tilde{Y}-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\tilde{B}}_j{\beta }_j\Vert }_2^2$ (2.5)

2.3. Weighted Least Square Estimation of GMCP Penalty

B-splines approximation is used on the unknown functions, which transforms the nonparameter regression into a parameter regression that makes variable selection and parameter estimation easier to solve. Meanwhile, the grouped variables in ${\tilde{B}}_j$ ; i.e., ${\tilde{B}}_{jk};k=1,\cdots ,M_n$ ; for each $j=\mathrm{1,}\cdots \mathrm{,}p$, are all related to the variable $X_j$ ; so we can consider B-spline basis functions for each nonparameter function $f_j$ to be a group. Instead of selecting the significant nonparameter functions, our task converts to choosing the significant B-spline basis functions from ${\tilde{B}}_j$ or nonzero coefficients from ${\beta }_j$.

In order to carry out variable selection at the group and individual variable levels simultaneously. In our case, the GMCP penalty function is

${\rho }_{\gamma }\left({\Vert {\beta }_j\Vert }_{A_j}\mathrm{,}\lambda \right)=\lambda {\displaystyle {\int }_0^{{\Vert {\beta }_j\Vert }_{A_j}}}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$ (2.6)

where $\gamma$ is a parameter that controls the concavity of $\rho$ and $\lambda$ is the penalty parameter. Here $x_{+}=x{1}_{\left\{x\ge 0\right\}}$. We require $\lambda \ge 0$ and $\gamma >1$. The term MCP comes from the fact that it minimizes the maximum concavity measure defined at (2.2) of [23], subject to conditions on unbiasedness and selection feature. The MCP can be easily understood by considering its derivative

${\stackrel{\dot{}}{\rho }}_{\gamma }\left({\Vert {\beta }_j\Vert }_{A_j}\mathrm{,}\lambda \right)=\lambda {\left(1-\frac{{\Vert {\beta }_j\Vert }_{A_j}}{\gamma \lambda }\right)}_{+}$ (2.7)

where for any $m\times 1$ vector $a$, ${\Vert a\Vert }_1$ is the L₁ norm: ${\Vert a\Vert }_1=\left|a_1+\cdots +a_m\right|$, $\lambda >0$ is the penalty tuning parameter and $A_j=\left\{k\mathrm{<mo>\; :\; </mo>}{\beta }_{jk}\in {\beta }_j\right\}$. In our case, each $A_j$ represents the j^th group of basis functions, i.e., ${\tilde{B}}_{jk}\mathrm{,}k=\mathrm{1,}\cdots \mathrm{,}{M}_n$ ; the values of the basis functions for each nonparameter function $f_j$ may be different from those for another function $f_{j^{\prime }}$ ; and when $j\ne j^{\prime }$ ; we assume there is no overlap between groups. Now combining the objective function in Equation (2.5) and the penalty function in Equation (2.6), we have the penalized weighted least squares objective function for the proposed NP-AFT-AR model as follows:

$Q_{n\lambda }\left(\beta \right)=\frac{1}{2}{\Vert \tilde{Y}-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\tilde{B}}_j{\beta }_j\Vert }_2^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}\lambda {\displaystyle {\int }_0^{{\Vert {\beta }_j\Vert }_{A_j}}}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$ (2.8)

We can conduct group or component selection and estimation by minimizing $Q_{n\lambda }\left(\beta \right)$ : If ${\Vert {\beta }_j\Vert }_{A_j}=0$ ; it implies that the function component $f_j$ is deleted, otherwise, it is selected, further, the individual basis functions within a group can be selected.

2.4. Computation

We derive a group coordinate descent algorithm for computing $\beta$. This algorithm is a natural extension of the standard coordinate descent algorithm ( [27] ). It has also been used in calculating the penalized estimates based on concave penalty functions ( [28] ).

The group coordinate descent algorithm optimizes a target function with respect to a single group at a time, iteratively cycling through all groups until convergence is reached. It is particularly suitable for computing $\beta$, since it has a simple closed form expression for a single-group model, see (2.11) below.

We write $A_j=R_j$ for an $M_n\times M_n$ upper triangular matrix $R_j$ via the Cholesky decomposition. Let ${\theta }_j=R_j{\beta }_j$ and ${\hat{B}}_j={\tilde{B}}_j{R}_j^{-1}$. Simple algebra shows that

$Q\left(\theta \mathrm{,}\lambda \mathrm{,}\gamma \right)=\frac{1}{2}{\Vert \tilde{Y}-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\hat{B}}_j{\theta }_j\Vert }_2^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}\lambda {\displaystyle {\int }_0^{\Vert {\theta }_j\Vert }}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$ (2.9)

Note that $n^{-1}{{\hat{B}}^{\prime }}_j{\hat{B}}_j={R^{\prime }}_j\left(n^{-1}{{\hat{B}}^{\prime }}_j{\hat{B}}_j\right)R_j^{-1}=I_{m_n}$. ${\hat{Y}}_j=\hat{Y}-{\displaystyle {\sum }_{k\ne j}^p}{\hat{B}}_k{\theta }_k$ and

$Q_j\left({\theta }_j\mathrm{,}\lambda \mathrm{,}\gamma \right)=\frac{1}{2}{\Vert {\tilde{Y}}_j-{\hat{B}}_j{\theta }_j\Vert }_2^2+\lambda {\displaystyle {\int }_0^{\Vert {\theta }_j\Vert }}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$ (2.10)

Let ${\eta }_j={\hat{B}}_j{\left({{\hat{B}}^{\prime }}_j{\hat{B}}_j\right)}^{-1}{\hat{Y}}_j$. For $\gamma >1$, it can be verified that the value that minimizes $Q_j\left(\theta \mathrm{,}\lambda \mathrm{,}\gamma \right)$ is

${\tilde{\theta }}_{j\mathrm{,}GM}\left(\lambda \mathrm{,}\gamma \right)=M\left({\eta }_j\mathrm{<mo>\; ;\; </mo>}\lambda \mathrm{,}\gamma \right)\equiv \{\begin{array}{ll}0\hfill & \text{if}\Vert {\eta }_j\Vert \le \lambda \hfill \\ \frac{\gamma }{\gamma -1}\left(1-\frac{\lambda }{\Vert {\eta }_j\Vert }\right){\eta }_j\hfill & \text{if}\lambda <\Vert {\eta }_j\Vert \le \gamma \lambda \hfill \\ {\eta }_j\hfill & \text{if}\Vert {\eta }_j\Vert >\gamma \lambda \hfill \end{array}$ (2.11)

In particular, when $\gamma =\infty$, we have

${\tilde{\theta }}_{j\mathrm{,}GL}={\left(1-\frac{\lambda }{\Vert {\eta }_j\Vert }\right)}_{+}{\eta }_j\mathrm{,}$

which is the GLasso estimate for a single-group model ( [29] ).

The group coordinate descent algorithm can now be implemented as follows. Suppose the current values for the group parameter ${\tilde{\theta }}_k^{\left(s\right)}\mathrm{,}k\ne j$ are given. We want to minimize $Q\left(\theta \mathrm{,}\lambda \mathrm{,}\gamma \right)$ with respect to ${\theta }_j$. Let

$Q_j\left({\theta }_j\mathrm{,}\lambda \mathrm{,}\gamma \right)=\frac{1}{2}{\Vert \tilde{Y}-{\displaystyle \underset{k\ne j}{\sum }}{\hat{B}}_k{\tilde{\theta }}_k^{\left(s\right)}-{\hat{B}}_j{\theta }_j\Vert }_2^2+\lambda {\displaystyle {\int }_0^{\Vert {\theta }_j\Vert }}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$ (2.12)

and write ${\tilde{Y}}_j={\displaystyle {\sum }_{k\ne j}}{\hat{B}}_k{\tilde{\theta }}_k^{\left(s\right)}$ and ${\eta }_j=n^{-1}{{\hat{B}}^{\prime }}_j\left(\tilde{Y}-{\tilde{Y}}_j\right)$. Let ${\tilde{\theta }}_j$ be the minimizer of $Q_j\left({\theta }_j\mathrm{,}\lambda \mathrm{,}\gamma \right)$. When $\gamma >1$, we have ${\tilde{\theta }}_j=M\left({\eta }_j\mathrm{,}\lambda \mathrm{,}\gamma \right)$, where M is defined in (2.11).

For any given $\left(\lambda \mathrm{,}\gamma \right)$, we use (2.11) to cycle through one component at a time. Let ${\tilde{\theta }}^{\left(0\right)}={\left({\tilde{\theta }}_1^{\left(0\right)}{}^{\prime }\mathrm{,}\cdots \mathrm{,}{\tilde{\theta }}_p^{\left(0\right)}{}^{\prime }\right)}^{\prime }$ be the initial value. The proposed coordinate descent algorithm is as follows.

Initial vector of residuals $r=Y-\tilde{Y}$, where $\tilde{Y}={\displaystyle {\sum }_{j=1}^p}{\hat{B}}_j{\theta }_j^{\left(0\right)}$, For $s=\mathrm{0,1,}\cdots$, carry out the following calculation until convergence. For $j=\mathrm{1,}\cdots \mathrm{,}p$, repeat the following steps.

Step 1: Calculate ${\tilde{\eta }}_j=n^{-1}{{\hat{B}}^{\prime }}_{j}r+{\tilde{\theta }}_j^{\left(s\right)}$.

Step 2: Update ${\theta }_j^{\left(s+1\right)}=M\left({\tilde{\eta }}_j\mathrm{<mo>\; ;\; </mo>}\lambda \mathrm{,}\gamma \right)$.

Step 3: Update $r=r-{\hat{B}}_j\left({\theta }_j^{\left(s+1\right)}-{\theta }_j^{\left(s\right)}\right)$ and $j=j+1$.

The last step ensures that r holds the current values of the residuals. Although the objective function is not necessarily convex, it is convex with respect to a single group when the coefficients of all the other groups are fixed.

3. Asymptotic Oracle Properties of GMCP

Let $\left|A\right|$ denote the cardinality of any set $A\in \left\{\mathrm{1,}\cdots \mathrm{,}p\right\}$ and $d_A=\left|A\right|M_n$. Define

$B_A=\left({\tilde{B}}_{jk}\mathrm{,}k=\mathrm{1,}\cdots \mathrm{,}{M}_n\mathrm{<mo>\; ;\; </mo>}j\in A\right)\text{and}{\Sigma }_A=\frac{{B^{\prime }}_A{B}_A}{n}$

Here $B_A$ is $n\times n{d}_A$ dimensional sub-design matrix corresponding to the variables in A, Denote ${\Vert f_j\Vert }_2=\left[E{f}_j^2\left(X_j\right)\right]$ We make the following assumptions.

Similar to [5], we assume:

(C1) T and C are independent.

(C2) $Pr\left(T\le C\mathrm{<mo>\; |\; </mo>}T\mathrm{,}X\right)=Pr\left(T\le C\mathrm{<mo>\; |\; </mo>}T\right)$.

(C3) $E\left(T^2\right)<\infty$ and $E\left(\epsilon \mathrm{<mo>\; |\; </mo>}X\right)=0$.

(C4) Denote ${\tau }_T$ and ${\tau }_C$ as the least upper bounds of T and C, respectively. Then ${\tau }_T<{\tau }_C$ or ${\tau }_T={\tau }_C=\infty$.

(C5) $f^2$ has finite envelope function.

(C6) $E{\left\{f\left(X\right)-f^0\left(X\right)\right\}}^2$ for $f\ne f^0$.

These assumptions correspond to the conditions in [30]. In the random censorship model, (C1) is a basic assumption. (C2) given the failure time T, the censoring indicator is independent of the $X$. (C3) in least-squares estimation, we need the second moment. (C4) assumes the probability of an event being observed is greater than zero, which guarantees the consistency of the estimator. (C5) is a fundamental condition for the consistency and convergence rate in the proofs, and is used in the entropy calculation. (C6) guarantees that $f_0$ is identifiable.

(C7) There exist constant $q^{*}>0,c_1>0$ and $c_2>0$ where $0<c_1\le c_2<\infty$ such that

$c_1\le \frac{{\Vert B_A\nu \Vert }_2^2}{n}\le c_2\mathrm{,}\forall \left|A\right|=q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\Vert \nu \Vert }_2=1\text{and}\nu \in R^{d_A}$

(C8) There is a small constant ${\eta }_1\ge 0$ such that ${\Sigma }_{k\in A_0}{\Vert f_j\Vert }_2\le {\eta }_1$.

(C9) The random errors ${\epsilon }_i,i=1,\cdots ,n$ are independent and identically distributed as $\epsilon$, where $E\left(\epsilon \right)=0$ and $E\left({\epsilon }^2\right)={\sigma }^2<\infty$ ; moreover, the tail probabilities satisfy $P\left(\left|\epsilon \right|>x\right)\le K\exp \left(-C{x}^2\right)$ for $x>0$ and some constantsC and K.

(C10) There exists a positive constant M such that $\left|x_{ik}\right|\le M,i=1,\cdots ,n;k=1,\cdots ,p$.

(C7) is the sparse Riesz condition (SRC) formulated for the nonparameter AFT model (2.1), which controls the range of eigenvalues of the matrix Z. This condition was introduced to study the properties of Lasso for the linear regression model by [31]. (C8) assumes that the unimportant predictors are small in the $L_2$ sense, but do not need to be exactly zero. If ${\eta }_1=0$, (C8) becomes $f_j=0$ for all $k\in A_0$. The problem of variable selection is equivalent to distinguishing nonzero functions from zero functions. (C9) assumes that the distribution of the error terms has sub-Gaussian tails. This condition holds when the error distribution is normal. (C10) assumes that all the predictors are uniformly bounded, which is satisfied in many practical situations.

In this subsection, we simply write ${\hat{f}}_j\left(X_j\right)={\displaystyle {\sum }_{k=1}^{M_n}}{\hat{\beta }}_{jk}{B}_{jk}$ is GMCP estimator. Let ${\beta }_{\ast }^o=\min \left\{{\Vert {\beta }_j^o\Vert }_2\mathrm{,}j\in A_0^c\right\}$ and set ${\beta }_{\ast }^o=\infty$ if $A_0^c$ is empty. Define

${\hat{\beta }}^o=\arg {\min }_b\left\{\frac{1}{2}{\Vert \tilde{Y}-{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\tilde{B}}_j{\beta }_j^o\Vert }_2^2\mathrm{<mo>\; ;\; </mo>}{\Vert {\beta }_j^o\Vert }_2=\mathrm{0,}\forall j\notin A_0^c\right\}$ (3.1)

and

$f_j^o\left(X\right)={\displaystyle \underset{k\mathrm{<mo>\; =\; </mo>1}}{\overset{M_n}{\sum }}}{\hat{\beta }}_{jk}^o{B}_{jk}$

This is the oracle least squares estimator. Of course, it is not a real estimator, since the oracle set is unknown.

We first consider the case where the 2-norm GMCP objective function is convex. This necessarily requires $c_{\min }>0$ where $c_{\min }$ be the smallest eigenvalue of $\Sigma$, and recall $\Sigma =n^{-1}{B}^{\prime }B$. As in [32], define the function

$h\left(t,k\right)=\exp \left(-k{\left(\sqrt{2t-1}-1\right)}^2/4\right),t>1,k=1,2,\cdots$ (3.2)

This function arises from an upper bound for the tail probabilities of the chi-square distributions given in Lemma A.2 in Appendix. This is derived from an exponential inequality for chi-square random variables of [33].

Theorem 3.1. Suppose ${\epsilon }_1\mathrm{,}\cdots \mathrm{,}{\epsilon }_n$ are independent and identically distributed as $N\left(\mathrm{0,}{\sigma }^2\right)$ and (C1)-(C10). Then for any $\left(\lambda \mathrm{,}\gamma \right)$ statisfying $\gamma >1/c_{\min }$ , ${\beta }_{\ast }^o>\lambda \gamma$ and $n{\lambda }^2>{\sigma }^2$ , we have

$P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}(\; \lambda \; )$

and

$P\left(\hat{f}\ne {\hat{f}}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}(\; \lambda \; )$

where ${\eta }_{1n}\left(\lambda \right)=\left(p-q\right)h\left(n{\lambda }^2/{\sigma }^2\mathrm{,}{M}_n\right)$ and ${\eta }_{2n}\left(\lambda \right)=qh\left(c_{1}n{\left({\beta }_{\ast }^o-\gamma \lambda \right)}^2/{\sigma }^2\mathrm{,}{M}_n\right)$.

We give the proof of Theorem 3.1 in Appendix. It provides an upper bound on the probability that $\hat{f}$ is not equal to the oracle estimator in terms of the tail probability functionh in (3.2). The key condition $\gamma >1/c_{\min }$ ensures that the 2-norm GMCP criterion is strictly convex. Nonetheless, this result is a starting for a similar result in $p>n$ case. The following corollary is an immediate consequence of Theorem 3.1.

Corollary 1. suppose that the condition of Therorm 3.1 are satisfied. Also suppose that ${\beta }_{\ast }^o\ge \gamma \lambda +a_n{\tau }_n$ for $a_n\to \infty$ as $n\to \infty$. If $\lambda \ge a_n{\lambda }_n$ , then

$P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\to 0\text{as}n\to \infty$

and

$P\left(\hat{f}\ne {\hat{f}}^o\right)\to 0\text{as}n\to \infty$

where

${\lambda }_n=\sigma \sqrt{2\log \left(\max \left\{p-q\mathrm{,1}\right\}\right)/\left(n{M}_n\right)}\text{and}{\tau }_n=\sigma \sqrt{2\log \left(\max \left(\left\{q\mathrm{,1}\right\}\right)/\left(n{c}_1{M}_n\right)\right)}$

By Corollary 1, the 2-norm GMCP estimator equals the oracle least squares estimator with probability converging to one. This implies it is group selection consistent. We now consider the high-dimensional case where $p>n$. Under condition (C7), let $K_{\ast }=\bar{c}-1/2$, $m_{\ast }=K_{\ast }q$ and $\xi =1/\left(4c^{*}{M}_n\right)$. Define

${\eta }_{3n}\left(\lambda \right)={\left(p-q\right)}^{m_{\ast }}\frac{e^{m_{\ast }}}{m_{\ast }^{m_{\ast }}}h\left(\xi n{\lambda }^2{\sigma }^{-2}/M_n,m_{\ast }{M}_n\right)$

Theorem 3.2. suppose ${\epsilon }_1\mathrm{,}\cdots \mathrm{,}{\epsilon }_n$ are independent and identically distributed as $N\left(\mathrm{0,}{\sigma }^2\right)$ and B satisfies the $SRC\left(q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{c}_1\mathrm{,}{c}_2\right)$ in (C7) with $q^{\mathrm{<mo>\; *\; </mo>}}\ge \left(\bar{c}-1/2\right)$ , $m_{\ast }=K_{\ast }q$ and $\xi =1/\left(4c^{\mathrm{<mo>\; *\; </mo>}}{M}_n\right)$ , we have

$P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}\left(\lambda \right)+{\eta }_{3n}(\; \lambda \; )$

and

$P\left(\hat{f}\ne {\hat{f}}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}(\; \lambda \; )$

where ${\eta }_{1n}\left(\lambda \right)=\left(p-q\right)h\left(n{\lambda }^2/{\sigma }^2\mathrm{,}{M}_n\right)$ and ${\eta }_{2n}\left(\lambda \right)=qh\left(c_{1}n{\left({\beta }_{\ast }^o-\gamma \lambda \right)}^2/{\sigma }^2\mathrm{,}{M}_n\right)$.

Corollary 2. suppose that the condition of Therorm 3.2 are satisfied. Also suppose that ${\beta }_{\ast }^o\ge \gamma \lambda +a_n{\tau }_n$ for $a_n\to \infty$ as $n\to \infty$ . If $\lambda \ge a_n{\lambda }_n^{\mathrm{<mo>\; *\; </mo>}}$ , then

$P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\to 0\text{as}n\to \infty$

and

$P\left(\hat{f}\ne {\hat{f}}^o\right)\to 0\text{as}n\to \infty$

where ${\lambda }_n^{\mathrm{<mo>\; *\; </mo>}}=2\sigma \sqrt{2c_2{M}_n\log \left(p-q\right)/n}$.

Theorem 3.2 and Corollary 2 provide sufficient conditions for the asymptotic oracle property of the global 2-norm GMCP estimator in the $p>n$ situations. Here we allow $p-\left|A_0^c\right|=\exp \left\{O\left(n/\left(c_2{M}_n\right)\right)\right\}$. So p can be greater than n. The condition $n{\lambda }^2\xi >{\sigma }^2{M}_n$ is stronger than the corresponding condition $n{\lambda }^2>{\sigma }^2$ in Theorem 3.5 ( [34] ). The condition $\gamma \ge c_1^{-1}\sqrt{4+\bar{c}}$ ensures that the GMCP criterion is convex in any q-dimensional subspace. It is stronger than the minimal sufficient condition $\gamma >c_1^{-1}$ for convexity in q-dimensional subspaces. This is the price we need to pay in search for a lower-dimensional space that contains the true model.

4. Numerical Simulation

In this section, we conduct simulation studies to evaluate the performance of the GMCP and GLasso penalties in a high-dimensional NP-AFT-AR model with limited samples. We therefore focus on the comparisons of the group selection methods with only the BIC ( [35] ) selected tuning parameter $\left(\lambda \mathrm{,}{M}_n\right)$, is given:

$\text{BIC}\left(\lambda \mathrm{,}{M}_n\right)=\log \left({\text{RSS}}_{\lambda \mathrm{,}{M}_n}\right)+\log \left(n\right)\frac{d{f}_{\lambda \mathrm{,}{M}_n}}{n}$

Where RSS is the sum of squared residuals, df is the number of selected variables given $\left(\lambda \mathrm{,}{M}_n\right)$. We choose $M_n$ from the increasing sequence in Section 5, for any given value of $M_n$, We choose from a sequence of 100 values $\lambda$, from $0.01{\lambda }_{\max }$ to ${\lambda }_{\max }$, Where ${\lambda }_{\max }={\max }_{1\le j\le p}{\Vert {{\tilde{B}}^{\prime }}_j\tilde{Y}\Vert }_2/\sqrt{M_n}{{\tilde{B}}^{\prime }}_j$ is corresponding to the covariate $X_j,j=1,\cdots ,p$ with $n\times M_n$ “design” matrix. ${\lambda }_{\max }$ is the maximum penalty value, which compresses all estimated coefficients to zero.

We compute the empirical prediction mean square error (MSE) to reveal the estimation accuracy. Let ${\hat{f}}_j$ be the estimator of $f_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}p$ ; and we define MSE as

${\text{MSE}}_{f_j}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left|{\hat{f}}_j\left(X_{ij}\right)-f_p\left(X_{ij}\right)\right|}^2$

Three scenarios are considered in the following, where some nonzero components are linear and the response variable is subject to various censoring rates. The sample size $n=400,200$ and a total of 100 simulation runs are used. The logarithm of censoring time $C_i$ is generated from a uniform distribution $U\left(c_1,c_2\right),c_1>0;c_2>0$, where $c_1$ and $c_2$ are determined by a Monte-Carlo method to achieve the censoring rates of 35% and 40% respectively. For example, the censoring rate $cr=Pr\left(T>C\right)$ is approximated by $\widehat{cr}={\displaystyle {\sum }_{i=1}^M}I\left(T_i\mathrm{<mo>\; >\; </mo>}{C}_i\right)/M$ where $T_i$ is drawn from the proposed model (2.1) and $C_i$ is drawn from $U\left(-c_1,c_2\right),c_1>0;c_2>0$, M is the Monte-Carlo simulation runs used to compute cr. When we chose $c_1=0,c_2=4,\widehat{cr}\approx 40\%$, which is considered to be the desired censoring rate. To take account of the computational efficiency and accuracy, we use the cubic B-spline with five evenly distributed interior knots for all the functions $f_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}p$, which gives the number of $3+1+5=9$ basis functions for each nonparametric component. Due to the identifiability constraint, $E\left\{f_j\left(X_j\right)\right\}=0$ ; the actual number of basis functions used is 8. This choice is made because our simulation studies indicated that using a larger number of knots does not improve the finite sample performance (results are not shown).

4.1. Scenario 1 (Covariates Are Independent)

In this scenario, we consider independent covariates and set the intercept ${\eta }_0=0$ : The logarithm of failure times, $T_i\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n$, are generated from

$\begin{array}{l}T=\exp (f_1\left(X_1\right)+f_2\left(X_2\right)+f_3\left(X_3\right)+f_4\left(X_4\right)+f_5\left(X_5\right)\begin{array}{c}\\ \end{array}\\ +f_6\left(X_6\right)+{\displaystyle \underset{j=7}{\overset{p}{\sum }}}{f}_j\left(X_j\right)+\epsilon )\end{array}$

where

$f_1\left(X_1\right)=2{\left(\sin \left(0.25\pi X_1\right)\right)}^3,f_2\left(X_2\right)=2\sin \left(2X_2\right),f_3\left(X_3\right)=X_3^2-\frac{3}{4},$

$f_4\left(X_4\right)=1.2X_4,f_5\left(X_5\right)=\exp \left(-X_5\right)-\frac{25}{12},$

$f_6\left(X_6\right)=\frac{1}{4}{X}_6^3,f_7\left(X_7\right)=\cdots =f_p\left(X_p\right)\equiv 0.$

The predictors are sampled from the $N\left(\mathrm{0,1}\right)$.

We set $p=500$ and consider the cases where $n=400,200$, respectively to see the performance of our proposed methods as the sample size increases. The penalty parameters are selected using CV as described above.

The results for the the GMCP, GSCAD and GLasso methods are given in Table 1 and Table 2 based on 100 replications. The columns in Table 1 include the average number of variables (NV) being selected, model error (ER), percentage of occasions on which correct variables are included in the selected model (%IN) and percentage of occasions on which the exactly correct variables are selected (%CS) with standard error in parentheses. Table 2 summarizes the mean square errors for the six important functions $n^{-1}{\displaystyle {\sum }_{i=1}^n}{\left|{\hat{f}}_j\left(X_{ji}\right)-f_j\left(X_{ji}\right)\right|}^2,j=1,\cdots ,6$ with standard error in parentheses.

Several observations can be obtained from Tables 1-4. The model that was selected by the GMCP and is better than the one selected by the GLasso in terms of model error, the percentage of occasions on which the true variables being selected and the mean square errors for the important coefficient functions. The GMCP includes the correct variables with high probability. When the sample size increases, the performance of both methods becomes better as expected. To examine the estimated nonparametric functions from Concave group Selection methods, we plot GMCP along with the true function components in Figure 1 and Figure 2. The estimated nonparametric coefficient functions are from GMCP method in one run when 100. From the graph, the estimators of the

Table 1. Simulation results. NV, number of selected variables; ER, model error; IN%, percentage of occasions on which the correct variables are included in the selected model; CS%, percentage of occasions on which exactly correct variables are selected, averaged over 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 2. Simulation results. Mean Square errors for the important coefficient functions based on 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 3. Simulation results. NV, number of selected variables; ER, model error; IN%, percentage of occasions on which the correct variables are included in the selected model; CS%, percentage of occasions on which exactly correct variables are selected, averaged over 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 4. Simulation results. Mean Square errors for the important coefficient functions based on 100 replications. Enclosed in parentheses are the corresponding standard errors.

Figure 1. $n=200$, the solid black line is the real function, the dotted red line is the GMCP estimation, CR = 35%.

Figure 2. $n=200$, the solid black line is the real function, the dotted red line is the GMCP estimation, CR = 40%.

nonparameter $f_j\left(X_j\right)\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,6}$, fit the true functions well, which are consistent with the mean square errors for the functions reported in Table 2, Table 4.

4.2. Scenario 2 (Covariates Are Correlated)

In this scenario, we consider correlated covariates and set the intercept ${\eta }_0=0$ : The logarithm of failure times, $T_i,i=1,\cdots ,n$, are generated from?

$\begin{array}{l}T=\exp (f_1\left(X_1\right)+f_2\left(X_2\right)+f_3\left(X_3\right)+f_4\left(X_4\right)+f_5\left(X_5\right)\begin{array}{c}\\ \end{array}\\ +f_6\left(X_6\right)+{\displaystyle \underset{j=7}{\overset{p}{\sum }}}{f}_j\left(X_j\right)+\epsilon )\end{array}$

$f_1\left(X_1\right)=1.2X_1,f_2\left(X_2\right)=2\sin \left(2X_2\right),f_3\left(X_3\right)=\left(X_3^2-\frac{3}{4}\right),$

$f_4\left(X_4\right)=\exp \left(-X_5\right)-\frac{25}{12},f_5\left(X_5\right)=\sin \left(0.5\pi X_5\right),$

$f_6\left(X_6\right)=2{\left(\sin \left(0.25\pi X_6\right)\right)}^3,f_7\left(X_7\right)=\cdots =f_p\left(X_p\right)\equiv 0.$

where the covariates $X=\left(X_1\mathrm{,}{X}_2\mathrm{,}\cdots \mathrm{,}{X}_p\right)$ are generated from $X_p=\left(W_p+0.5U\right)/1.5$ where $W_1\mathrm{,}\cdots ,W_p$ and U are i.i.d. $N\left(\mathrm{0,1}\right)$. This provides a design with a correlation coefficient of 0.5 between all of the covariates.

The simulation study results are reported in Tables 5-8. The conclusions for Scenario 2 are very similar to those for Scenario 1. When the censoring rate increases, the estimation and selection performance decreases for all methods. The results in Table 6, Table 8 show that the GMCP estimator is more accurate than the GLasso estimator for both the individual component functions and the full

Table 5. Simulation results. NV, number of selected variables; ER, model error; IN%, percentage of occasions on which the correct variables are included in the selected model; CS%, percentage of occasions on which exactly correct variables are selected, averaged over 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 6. Simulation results. Mean Square errors for the important coefficient functions based on 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 7. Simulation results. NV, number of selected variables; ER, model error; IN%, percentage of occasions on which the correct variables are included in the selected model; CS%, percentage of occasions on which exactly correct variables are selected, averaged over 100 replications. Enclosed in parentheses are the corresponding standard errors.

Table 8. Simulation results. Mean Square errors for the important coefficient functions based on 100 replications. Enclosed in parentheses are the corresponding standard errors.

model, since the MSE under the GMCP approach is always smaller than that under the GLasso approach. The results in Table 5, Table 7 show that the GMCP method conducts component selection more precisely than the GLasso method, while the GLasso method chooses many zero component functions as nonzero functions. To examine the estimated nonparametric functions from the GMCP, we plot them along with the true function components in Figure 3, Figure 4. The estimated functions are from the GMCP method in one run when $n=200$. The estimation and selection accuracy decrease when covariates are correlated, we can still see that the estimated curves under the GMCP method are close to the true curves compared with the estimated curves under the GLasso method.

5. Application in NA-AFT-Model

In this section, we will use Shedden 2008 (for short) to conduct an empirical analysis of part of the collected lung adenocarcinoma data to illustrate the proposed method. For more information, see [36]. Retrospective data of 442 lung adenocarcinoma patients were collected at multiple locations, including their survival time, some other clinical and demographic data, and the expression level of the 22,283 gene from the following genes: tumor samples. However, most samples have small changes. Therefore, in our application, we randomly select 321 samples, the first 500, 1000 genes. Therefore, $n=321,p=1000$, and the survival rate is 35.8%.

Figure 3. $n=200$, the solid black line is the real function, the dotted red line is the group MCP estimation, the dotted blue line is the group SCAD estimation, and the black line is the group lasso estimation, CR = 35%.

Figure 4. $n=200$, the solid black line is the real function, the dotted red line is the group MCP estimation, the dotted blue line is the group SCAD estimation, and the black line is the group lasso estimation, CR = 40%.

Figure 5. The estimation function graph based on GLasso and GMCP approximates the corresponding 200746_s_at by using the same covariate, where the red dotted line is the GMCP estimate and the gray dotted line is the GLasso estimate.

Here, we are interested in the effect of tumor gene expression levels on the survival time of lung adenocarcinoma patients. Since the linear assumption is always latent in high dimensions, the proposed method may be more suitable for analyzing feature selection problems considering nonlinear effects. In our analysis, we set the spline base $M_n=5$ for each gene. The proposed method selects 1 gene locus under GMCP (ie 200746_s_at). However, when $p=500,1000$, the method under GLasso penalized regression alone selected the 6, 10 gene.

From Figure 5, we find that the larger the dimension, the worse the GLasso method estimation, but it has little effect on the GMCP estimation. Therefore, the verification of the actual data shows that the GMCP penalty is better than the GLasso penalty, and the accuracy is higher, and the calculation cost of the two is the same. Under the same conditions, the GMCP method is more suitable than the GLasso.

6. Concluding Remarks

In this paper, we study the weighted least squares estimation and selection attributes of GMCP in the NP-AFT-AR model with high-dimensional data. For the GMCP method, our simulation results show that GLasso tends to select some unimportant variables. In contrast, GMCP has progressive predictability, which shows that it also has selection consistency.

Appendix Proof

Lemma 1. Let ${\chi }_k^2$ be a random variable with chi-square distribution with k degrees of freedom. For $t>1$, $P\left({\chi }_k^2\ge kt\right)\le h\left(t\mathrm{,}k\right)$, where $h\left(t\mathrm{,}k\right)$ is defined in (3.2).

This lemma is a restatement of the exponential inequality for chi-square distributions of [33].

proofof Theorem 3.1. Since ${\hat{\beta }}^o$ is the oracle least squares estimator, we have ${\hat{\beta }}_j^o\mathrm{,}j\in A_0$ and

$-{{\tilde{B}}^{\prime }}_j\left(\tilde{Y}-\tilde{B}{\hat{\beta }}^o\right)/n=\mathrm{0,}\forall j\in A_0^c$

If ${\Vert {\hat{\beta }}_j^o\Vert }_2/\sqrt{M_n}\ge \lambda \gamma$, then by the definition of the MCP, $\stackrel{\dot{}}{\rho }\left({\Vert {\hat{\beta }}_j^o\Vert }_2\mathrm{<mo>\; ;\; </mo>}\sqrt{M_n}\lambda \mathrm{,}\gamma \right)=0$. Since $c_{\min }>1/\gamma$, the criterion (2.8) is strictly convex. By the Karush-Kuhn-Tucker (KKT) conditions, the equality $\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)={\hat{\beta }}^o$ holds in the intersection of the events

${\Omega }_1\left(\lambda \right)=\left\{\underset{j\in A_0}{\max }{\Vert n^{-1}{{\tilde{B}}^{\prime }}_j\left(\tilde{Y}-\tilde{B}\right){\hat{\beta }}^o\Vert }_2/\sqrt{M_n}\le \lambda \right\}$

and

${\Omega }_2\left(\lambda \right)=\left\{\underset{j\in A_0^c}{\min }{\Vert {\hat{\beta }}_j^o\Vert }_2\ge \gamma \lambda \right\}$

We first bound $1-P\left({\Omega }_1\left(\lambda \right)\right)$. Let ${\hat{\beta }}_{A_0^c}={\left({\hat{\beta }}_j\mathrm{,}j\in A_0^c\right)}^{\prime }$. By (A.1) [34] and using $\tilde{Y}={\tilde{B}}_{A_0^c}{\beta }_{A_0^c}^o+\epsilon$

${\hat{\beta }}_{A_0^c}^o={\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}\tilde{Y}/n={\beta }_{A_0^c}^o+{\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon /n$

It follows that $n^{-1}{\tilde{B}}_j\left(\tilde{Y}-\tilde{B}{\hat{\beta }}^o\right)=n^{-1}{{\tilde{B}}^{\prime }}_k\left(I_n-P_{A_0^c}\right)\epsilon$, where $P_{A_0^c}=n^{-1}{\tilde{B}}_{A_0^c}{\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}$, Because ${{\tilde{B}}^{\prime }}_j{\tilde{B}}_j=I_{M_n}$, ${\Vert {\tilde{B}}_j\left(I_n-P_{A_0^c}\right)\epsilon \Vert }_2^2/{\sigma }^2$ is distributed as a ${\chi }^2$ distribution with $M_n$ degrees of freedom. We have, for $n{\lambda }^2/{\sigma }^2\ge 1$

$\begin{array}{c}1-P\left({\Omega }_1\left(\lambda \right)\right)=P\left(\underset{j\in A_0}{\max }{\Vert n^{1/2}{{\tilde{B}}^{\prime }}_j\left(I_n-P_{A_0^c}\right)\epsilon \Vert }_2^2/\left(M_n{\sigma }^2\right)>n{\lambda }^2/{\sigma }^2\right)\\ \le {\displaystyle \underset{j\in A_0}{\sum }}P\left({\Vert n^{-1/2}{{\tilde{B}}^{\prime }}_j\left(I_n-P_{A_0^c}\right)\epsilon \Vert }_2^2/\left(\sqrt{M_n}{\sigma }^2\right)>M_{n}n{\lambda }^2/{\sigma }^2\right)\\ \le {\displaystyle \underset{j\in A_0}{\sum }}h\left(n{\lambda }^2/{\sigma }^2\mathrm{,}{M}_n\right)\\ \le \left(p-q\right)h\left(n{\lambda }^2/{\sigma }^2\mathrm{,}{M}_n\right)\\ ={\eta }_{1n}\left(\lambda \right)\end{array}$ (6.1)

where we used lemma 1 in the third line. Now consider ${\Omega }_2\left(\lambda \right)$, Recall ${\beta }_{\ast }^o={\min }_{j\in A_0^c}{\Vert {\beta }_j^o\Vert }_2$. If ${\Vert {\hat{\beta }}_j^o-{\beta }_j^o\Vert }_2/\sqrt{M_n}\le {\beta }_{\ast }^o-\gamma \lambda$ for all $j\in A_0^c$, then ${\min }_{j\in A_0^c}{\Vert {\hat{\beta }}_j^o\Vert }_2/\sqrt{M_n}\ge \gamma \lambda$. This implies

$1-P\left({\Omega }_2\left(\lambda \right)\right)\le P\left(\underset{j\in A_0^c}{\max }{\Vert {\hat{\beta }}_j^o-{\beta }_j^o\Vert }_2/\sqrt{M_n}>{\beta }_{\ast }^o-\gamma \lambda \right)$

Let ${\tilde{B}}_j$ be a $M_n\times M_{n}q$ matrix with a $M_n\times M_n$ identity matrix $I_{M_n}$ in the pth block and 0’s elsewhere. Then $n^{1/2}\left({\hat{\beta }}_j^o-{\beta }_j^o\right)=n^{-1/2}{\tilde{B}}_j{\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon$. Note that

${\Vert n^{-1/2}{\tilde{B}}_j{\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2\le {\Vert {\tilde{B}}_j\Vert }_2{\Vert {\Sigma }_{A_0^c}^{-1/2}\Vert }_2{\Vert n^{-1/2}{\Sigma }_{A_0^c}^{-1/2}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2\le c_1^{-1/2}{\Vert n^{-1/2}{\Sigma }_{A_0^c}^{-1/2}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2$

and ${\Vert n^{-1/2}{\Sigma }_{A_0^c}^{-1/2}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2^2/{\sigma }^2$ id distributed as a $\chi$ distribution with q degrees of freedom. Therefore, similar to ${\eta }_{1n}\left(\lambda \right)$, we have, for $c_{1}n\left({\beta }_{\ast }^o-\gamma \lambda \right)/{\sigma }^2>1$,

$\begin{array}{c}1-P\left({\Omega }_2\left(\lambda \right)\right)=P\left(\underset{j\in A_0^c}{\max }{n}^{-1/2}{\Vert {\tilde{B}}_j{\Sigma }_{A_0^c}^{-1}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2/\sqrt{M_n}>\sqrt{n}\left({\beta }_{\ast }^o-\gamma \lambda \right)\right)\\ \le P\left(\underset{j\in A_0^c}{\max }{\Vert n^{-1/2}{\Sigma }_{A_0^c}^{-1/2}{{\tilde{B}}^{\prime }}_{A_0^c}\epsilon \Vert }_2^2/\left(M_n{\sigma }^2\right)>c_{1}n{\left({\beta }_{\ast }^o-\gamma \lambda \right)}^2/{\sigma }^2\right)\\ \le qh\left(c_{1}n{\left({\beta }_{\ast }^o-\gamma \lambda \right)}^2/{\sigma }^2\mathrm{,}{M}_n\right)\\ ={\eta }_{2n}\left(\lambda \right)\end{array}$ (6.2)

Combining ${\eta }_{1n}\left(\lambda \right)$ and ${\eta }_{2n}\left(\lambda \right)$, we have

$P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\le 1-P\left({\Omega }_1\left(\lambda \right)\right)+1-P\left({\Omega }_2\left(\lambda \right)\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}(\; \lambda \; )$

Since $\hat{f}\left(x\right)=B\hat{\beta }$, we can obtain $P\left(\hat{f}\ne {\hat{f}}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}\left(\lambda \right)$. This completes the proof.

For any $Q\subset \left\{\mathrm{1,}\cdots \mathrm{,}p\right\}$ and $m\ge 1$, define

$\zeta \left(\nu \mathrm{<mo>\; ;\; </mo>}m\mathrm{,}B\right)=\max \left\{\frac{{\Vert \left(P_A-P_B\right)\nu \Vert }_2}{{\left(mn\right)}^{1/2}}\mathrm{<mo>\; :\; </mo>}Q\subseteq A\subseteq \left\{\mathrm{1,}\cdots \mathrm{,}p\right\}\mathrm{,}{d}_A=m+d_B\right\}$

Lemma 2. Suppose $\xi n{\lambda }^2>{\sigma }^2{M}_n$. We have

$P\left(2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y};m,A_0^c\right)>\lambda \right)\le {\left(p-q\right)}^m\frac{e^m}{m^m}\exp \left(-m\xi n{\lambda }^2/16\right)$

proof. For any $A\supseteq A_0^c$. We have $\left(P_a-P_{A_0^c}\right){\tilde{B}}_{A_0^c}{\hat{\beta }}_{A_0^c}=0$. Thus

$\left(P_A-P_{A_0^c}\right)\tilde{Y}=\left(P_A-P_{A_0^c}\right)\left({\tilde{B}}_{A_0^c}{\beta }_{A_0^c}+\epsilon \right)=\left(P_A-P_{A_0^c}\right)\epsilon$

Therefore,

$P\left(2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y}\mathrm{<mo>\; ;\; </mo>}m\mathrm{,}{A}_0^c\right)>\lambda \right)=P\left(\underset{A\supseteq A_0^c,\left|A\right|=m+q}{\max }{\Vert \left(P_A-P_{A_0^c}\right)\epsilon \Vert }_2^2/{\sigma }^2>\xi mn{\lambda }^2\right)$

Since $P_A-P_B$ is a projection matrix, ${\Vert \left(P_A-P_{A_0^c}\right)\epsilon \Vert }_2^2/{\sigma }^2\sim {\chi }_{m_A}^2$, where $m_A={\displaystyle {\sum }_{j\in A-A_0^c,A\supseteq A_0^c}}{M}_n\le m{M}_n$. Since there $\left(\begin{array}{c}p-q\\ m\end{array}\right)$ are ways to choose A from $\left\{\mathrm{1,}\cdots \mathrm{,}p\right\}$, we have

$P\left(2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y};m,A_0^c\right)>\lambda \right)\le \left(\begin{array}{c}p-q\\ m\end{array}\right)P\left({\chi }_{m{M}_n}^2>\xi mn{\lambda }^2\right).$

This and Lemma A.2 imply that

$\begin{array}{c}P\left(2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y};m,A_0^c\right)>\lambda \right)\le \left(\begin{array}{c}p-q\\ m\end{array}\right)h\left(\xi n{\lambda }^2/M_n,m{M}_n\right)\\ \le {\left(p-q\right)}^m\frac{e^m}{m^m}h\left(\xi n{\lambda }^2/M_n,m{M}_n\right)\end{array}$ (6.3)

here we used the inequality $\left(\begin{array}{c}p-q\\ m\end{array}\right)\le {\left(p-q\right)}^m\frac{e^m}{m^m}$, this completes the proof.

Define I as any set that satisfies

$\begin{array}{l}{A}_0^c\cup \left\{j\mathrm{<mo>\; :\; </mo>}{\Vert {\hat{\beta }}_j\Vert }_2\ne 0\right\}\subseteq I\\ \subseteq A_0^c\cup \left\{j\mathrm{<mo>\; :\; </mo>}{n}^{-1}{\tilde{B}}^{\prime }\left(\tilde{Y}-\tilde{B}\hat{\beta }\right)=\stackrel{\dot{}}{\rho }\left({\Vert {\hat{\beta }}_j\Vert }_2\mathrm{<mo>\; ;\; </mo>}\sqrt{M_n}\lambda \mathrm{,}\gamma \right)\sqrt{M_n}{\hat{\beta }}_j/{\Vert {\hat{\beta }}_j\Vert }_2\right\}\end{array}$

Lemma 3. Suppose that $\tilde{B}$ satisfies that $SRC\left(q^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{c}_1\mathrm{,}{c}_2\right)$ , $q^{\mathrm{<mo>\; *\; </mo>}}\ge \left(K_{\ast }+1\right)m_{n}q$, and $\gamma \ge c_1^{-1}\sqrt{4+\bar{c}}$ . Let $m_{\ast }=K_{\ast }q$ . Then for any $\tilde{Y}\in {\mathcal{R}}^n$ with $\lambda \ge 2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y}\mathrm{<mo>\; ;\; </mo>}{m}_{\ast }\mathrm{,}{A}_0^c\right)$ , we have $\left|I\right|\le \left(K_{\ast }+1\right)q$ .

proof This lemma can be proved along the line of the proof of Lemma 1 of [23] and is omitted. proofof Theorem 3.2. By Lemma 3, in the event

$2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y}\mathrm{<mo>\; ;\; </mo>}{m}_{\ast }\mathrm{,}{A}_0^c\right)\le \lambda$ (6.4)

we have $\left|I\right|\le \left(K_{\ast }+1\right)q$, Thus in the event (6.4), the original model with p groups reduces a model with at most $\left(K_{\ast }+1\right)q$ groups, in this reduced model, the condition of Theorem 3.2 implies that the conditions of Theorem 3.2. By Lemma 2,

$P\left(2\sqrt{c_2{M}_n}\zeta \left(\tilde{Y}\mathrm{<mo>\; ;\; </mo>}{m}_{\ast }\mathrm{,}{A}_0^c\right)\le \lambda \right)\le {\eta }_{3n}\left(\lambda \right)$ (6.5)

Therefore, combining (6.5) and Theorem 3.1, we have $P\left(\hat{\beta }\left(\lambda \mathrm{,}\gamma \right)\ne {\hat{\beta }}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}\left(\lambda \right)+{\eta }_{3n}\left(\lambda \right)$, since $\hat{f}\left(x\right)=B\hat{\beta }$, we can obtain $P\left(\hat{f}\ne {\hat{f}}^o\right)\le {\eta }_{1n}\left(\lambda \right)+{\eta }_{2n}\left(\lambda \right)+{\eta }_{3n}\left(\lambda \right)$. This proves Theorem 3.2.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Cox, D.R. (1972) Regression Models and Life Tables (with Discussion). Journal of the Royal Statistical Society: Series B, 2, 187-220. https://doi.org/10.1111/j.2517-6161.1972.tb00899.x
[2]	Huang, J., Ma, S. and Xie, H.L. (2006) Regularized Estimation in the Accelerated Failure Time Model with High-Dimensional Covariates. Biometrics, 62, 813-820. https://doi.org/10.1111/j.1541-0420.2006.00562.x
[3]	Datta, S., Jennifer, L.R. and Datta, S. (2007) Predicting Patient Survival from Microarray Data by Accelerated Failure Time Modeling Using Partial Least Squares and LASSO. Biometrics, 63, 259-271. https://doi.org/10.1111/j.1541-0420.2006.00660.x
[4]	Huang, J., Ma, S. and Xie, H.L. (2007) Least Absolute Deviations Estimation for the Accelerated Failure Time Model. Statistica Sinica, 17, 1533-1548.
[5]	Leng, C. and Ma, S. (2007) Accelerated Failure Time Models with Nonlinear Covariates Effects. Australian and New Zealand Journal of Statistics, 49, 155-172. https://doi.org/10.1111/j.1467-842X.2007.00470.x
[6]	Schmid, M. and Hothorn, T. (2008) Flexible Boosting of Accelerated Failure Time Models. BMC Bioinformatics, 9, Article No. 269. https://doi.org/10.1186/1471-2105-9-269
[7]	David, E. and Li, Y. (2009) Survival Analysis with High-Dimensional Covariates: An Application in Microarray Studies. Statistical Applications in Genetics and Molecular Biology, 8, 1-122. https://doi.org/10.2202/1544-6115.1423
[8]	Cai, T., Huang, J. and Tian, L. (2009) Regularized Estimator for the Accelerated Failure Time Model. Biometric, 65, 394-404. https://doi.org/10.1111/j.1541-0420.2008.01074.x
[9]	Johnson, B.A.(2009) Rank-Based Estimation in the l1-Regularized Partly Linear Model for Censored Outcomes with Application to Integrated Analyses of Clinical Predictors and Gene Expression data. Biostatistics, 10, 3659-466.
[10]	Huang, J. and Ma, S. (2010) Variable selection in the accelerated failure time model via the bridge method. Lifetime Data Analysis, 16, 176-195. https://doi.org/10.1007/s10985-009-9144-2
[11]	Wang, Z. and Wang, C.Y. (2010) Buckley-James Boosting for Survival Analysis with High-Dimensional Biomarker Data. Statistical Applications in Genetics and Molecular Biology, 9, Article No. 24. https://doi.org/10.2202/1544-6115.1550
[12]	Liu, F., Dunson, D. and Zou, F. (2011) High-Dimensional Variable Selection in Meta-Analysis for Censored Data. Biostatistics, 67, 504-512. https://doi.org/10.1111/j.1541-0420.2010.01466.x
[13]	Ma, S. and Du, P. (2012) Variable Selection in Partly Linear Regression Model with Diverging Dimensions for Right Censored Data. Statistica Sinica, 22, 1003-1020. https://doi.org/10.5705/ss.2010.267
[14]	Schumaker, L. (1981) Spline Functions: Basic Theory. Cambridge University Press, Cambridge.
[15]	Gu, C. (2004) Model Diagnostics for Smoothing Spline ANOVA Models. Canadian Journal of Statistic, 32, 347-358. https://doi.org/10.2307/3316020
[16]	Hu, J.W. and Chai, H. (2013) Adjusted Regularized Estimation in the Accelerated Failure Time Model with High Dimensional Covariates. Journal of Multivariate Analysis, 122, 96-114. https://doi.org/10.1016/j.jmva.2013.07.011
[17]	Chai, H., Liang, Y. and Liu, X.Y. (2015) The L1/2 Regularization Approach for Survival Analysis in the Accelerated Failure Time Model. Computers in Biology and Medicine, 64, 283-290. https://doi.org/10.1016/j.compbiomed.2014.09.002
[18]	Xia, X.C., Jiang, B.Y., Li, J.L. and Zhang, W.Y. (2016) Low-Dimensional Confounder Adjustment and High-Dimensional Penalized Estimation for Survival Analysis. Lifetime Data Analysis, 22, 547-569. https://doi.org/10.1007/s10985-015-9350-z
[19]	Khan, M.H.R. and Shaw, J.E.H. (2016) Variable Selection for Survival Data with a Class of Adaptive Elastic Net Techniques. Statistics and Computing, 26, 725-741. https://doi.org/10.1007/s11222-015-9555-8
[20]	Yang, Y.C., Fasching, P.A. and Tresp, V. (2017) Modeling Progression Free Survival in Breast Cancer with Tensorized Recurrent Neural Networks and Accelerated Failure Time Models. Proceedings of Machine Learning Research, 68, 164-176.
[21]	Yue, M., Li, J.L. and Ma, S.G. (2018) Sparse Boosting for High-Dimensional Survival Data with Varying Coefficients. Statistics in Medicine, 37, 789-800. https://doi.org/10.1002/sim.7544
[22]	Khan, M.H.R. and Shaw, J.E.H. (2019) Variable Selection for Accelerated Lifetime Models with Synthesized Estimation Techniques. Statistical Methods in Medical Research, 28, 937-952. https://doi.org/10.1177%2F0962280217739522
[23]	Zhang, C.H. (2010) Nearly Unbiased Variable Selection under Minimax Concave Penalty. The Annals of Statistics, 38, 894-942. https://doi.org/10.1214/09-AOS729
[24]	Stute, W. (1993) Almost Sure Representations of the Product-Limit Estimator for Truncated Datay. The Annals of Statistics, 21, 146-156. https://doi.org/10.1214/aos/1176349019
[25]	Zhou, M. (1992) M-Estimation in Censored Linear Models. Biometrika, 79, 837-841.
[26]	Satten, G.A. and Datta, S. (2001) The Kaplan-Meier Estimator as an Inverse-Probability-of-Censoring Weighted Average. The American Statistician 55, 207-210. https://doi.org/10.1198/000313001317098185
[27]	Fu, W.J. (1998) Penalized Regressions: The Bridge versus the LASSO. Journal of Computational and Graphical Statistics, 7, 397-416. https://doi.org/10.1080/10618600.1998.10474784
[28]	Breheny, P. and Huang, J. (2011) Coordinate Descent Algorithms for Nonconvex Penalized Regression, with Applications to Biological Feature Selection. The Annals of Statistics, 5, 232-253. https://doi.org/10.1214/10-AOAS388
[29]	Yuan, M. and Lin, Y. (2006) Model Selection and Estimation in Regression with Grouped Variables. Journal of the Royal Statistical Society: Series B, 68, 49-67. https://doi.org/10.1111/j.1467-9868.2005.00532.x
[30]	Stute, W. (1999) Nonlinear Censored Regression. Statistica Sinica, 9, 1089-1102.
[31]	Zhang, C.H. and Huang, J. (2008) The Sparsity and Bias of the Lasso Selection in High-Dimensional Linear Regression. The Annals of Statistics, 36, 1567-1594. https://doi.org/10.1214/07-AOS520
[32]	Huang, J., Breheny, P. and Ma, S. (2012) A Selective Review of Group Selection in High-Dimensional Models. Statistical Science, 27, 481-499. https://doi.org/10.1214/12-STS392
[33]	Laurent, B. and Laurent, P. (2000) Adaptive Estimation of a Quadratic Functional by Model Selection. The Annals of Statistics, 28, 1302-1338. https://doi.org/10.1214/aos/1015957395
[34]	Yang, G.R., Huang, J. and Zhou, Y. (2014) Concave Group Methods for Variable Selection and Estimation in High-Dimensional Varying Coefficient Models. Science china Mathematics, 57, 2073-2090. https://doi.org/10.1007/s11425-014-4842-y
[35]	Schwarz, G. (1978) Estimating the Dimension of a Model. The Annals of Statistics, 6, 461-464. https://doi.org/10.1214/aos/1176344136
[36]	Director’s Challenge Consortium for the Molecular Classification of Lung Adenocarcinoma (2008) Gene Expression-Based Survival Prediction in Lung Adenocarcinoma: A Multi-Site, Blinded Validation Study. Nature Medicine, 14, 822-827. https://doi.org/10.1038/nm.1790