Sufficient variable screening via directional regression with censored response
Menghao Xu Zhou Yu Jun Shao

#1 1	





Sufficient variable screening  via directional regression with censored response
    Menghao Xu^1, Zhou Yu^1, Jun Shao^1,2

	^1 School of Statistics, East China Normal University, Shanghai, China 

	^2Department of Statistics, University of Wisconsin-Madison, WI, 53706, USA
    December 30, 2023
=================================================================================================================================================================================================





	We in this paper propose a directional regression based approach for ultrahigh dimensional sufficient variable screening with censored responses. The new method is designed in a model-free manner and thus can be adapted to various complex model structures. Under some commonly used assumptions, we show that the proposed method enjoys the sure screening property when the dimension p diverges at an exponential rate of the sample size n. To improve the marginal screening method, the corresponding iterative screening algorithm and stability screening algorithm are further equipped.
	We demonstrate the effectiveness of the proposed method through simulation studies and a real data analysis.
	
	Key Words: Sufficient dimension reduction, Sufficient Variable Selection; Sure independence screening; Ultrahigh dimensional covariates.
	


1.4



§ INTRODUCTION


Data sets collected in many contemporary scientific areas are ultrahigh dimensional and too complex to be analyzed through classical statistical methods.




Consider data observed from a random sample of size n from the 
distribution of (, ), where 
 is a scalar response,  =(_1,…,_p)
is a p-dimensional column vector of covariates, and 
the joint distribution of ( ,  ) is fully nonparametric. 
With an ultrahigh dimension p >> n, it is of great interest to identify 
⊂{1,...,p} such that _ = { x_k: k ∈} 
is truly related to the response. To fulfill the goal of model-free variable selection
based on the training data, 
<cit.> introduced the concept of sufficient variable selection as finding the smallest covariate set _ with ⊂{1,...,p} satisfying 

    |_,

where  stands for independence and | stands for conditioning. For convenience, in what follows we name both  ⊆{1,...,p} and  _= {x_k: k ∈} as covariate set.   
If it is too hard to find the smallest covariate set  satisfying (<ref>) especially when p >> n, a weaker goal is to find a covariate set  
containing  with size as small as possible, which is referred to as sufficient variable screening and is the focus of this paper.  

Research on sufficient variable screening in ultrahigh dimensional setting has gained considerable momentum in recent years. <cit.> and <cit.> proposed to use marginal distance correlation and marginal martingale difference divergence for sufficient variable screening.  Noticing the close relationship between sufficient variable selection and sufficient dimension reduction (; ), <cit.>, <cit.> and <cit.> developed different dimension reduction based screening methods.

In many biomedical studies, the response are often censored rather than fully observed. We consider survival data in which  is the true lifetime,  is the censoring time and 
we only observe ^o= min{,} and the censoring indicator δ=I(≤). 
Sufficient variable selection with censored response is finding  in (<ref>) with Y replaced by (, ), i.e., 

    (,) |_.

While our focus is (<ref>), we can only  observe (^o , δ ), instead of (T,C).

There exists very limited amount of work on model-free variable screening with censored responses.  Assuming |, the quantile  adaptive sure independence screening procedure proposed by 
<cit.> can be naturally extended to survival analysis. <cit.> proposed a survival impact index, which characterizes the impacts of a covariate on the distribution of true lifetime  by evaluating the absolute deviation of the covariate-stratified survival distribution from the unstratified survival distribution. The proposed survival impact index based screening seems to take some advantages over quantile adaptive sure independence screening when dealing with censored responses.

We in this paper give a modification of the directional regression () capable for sufficient dimension reduction with censored response
, and then characterize a suitable modified directional regression index for sufficient variable screening. The sure screening property is established in the ultrahigh dimensional setting, i.e., with probability tending to one, the smallest covariate set  is contained in the set of covariates selected by our proposed procedure. 
We also discuss the limitations of such modified directional regression index and propose a refined iterative procedure of our screening approach. To further enhance the stability of variable screening, we follow  <cit.> and <cit.> to  integrate the resampling scheme into our proposal. 
After screening,  the selected covariate set may contain some unrelevant covariates,   but its size  is much smaller than n so that we may apply variable selection or dimension reduction using an existing method to further reduce the size or dimension of the selected covariate set. 
Our approach are examined through simulation studies and an application to the diffuse large-B-cell lymphoma microarray data <cit.>.




§ MODIFIED DIRECTIONAL REGRESSION INDEX



To derive an index for covariate screening, we first reveal a relationship between sufficient variable selection and sufficient dimension reduction, another perspective in reducing covariate dimension. As a by-product,  we extend one method in sufficient dimension reduction, the directional regression, to survival data with censoring, which leads to an index for sufficient variable screening.



Sufficient dimension reduction aims to identify a linear function of  with dimension lower than p,  without losing information. To be specific, we seek a p× d matrix  with the smallest d such that  

    (,C) |.

The linear space generated by columns of  is called the central subspace and denoted as 𝒮_(,)|. The following result reveals a deep connection between sufficient variable selection (<ref>) and sufficient dimension reduction (<ref>) for censored responses. 	

	Let _1,…,_d be columns of  in (<ref>) and _k be the p× 1 vector whose kth element is 1 all other elements are 0. 		 Then, 
	∑_j=1^d|_k_j|>0 for k ∈ and  ∑_j=1^d|_k_j|=0 for k ∉, where 
	 is given in (<ref>). 




This result tells us that  in sufficient dimension reduction can be also used for sufficient variable selection.  
Inspired by this, in the following we first extend the directional regression <cit.> to find the central space 𝒮_(,)| using survival data with censoring. 

Let 
=^-1/2(-) be the standardized covariate, where 
=E() and =(). Then, 
𝒮_(,)|=^-1/2𝒮_(,)|.   
A key result for the success of  the directional regression in  <cit.> is that, 
if (T,C) is observed,  the  column space  of E[ 2_p-
E{(-)(-)|,,,}]^2 is equal to   
𝒮_(,)|, where _p is the identity matrix of order p and 
(,,) is an independent copy of (,,).
However, in survival analysis  (,) is unobservable; instead, we observe  (^o,δ). The next proposition extends the result in <cit.> to the survival data with censoring. 


	Let =E[2_p-E{(-)(-)|^o ,^o,δ,δ}]^2, where 
	(,^o,δ) is an independent copy of (,^o,δ ).
	

(i) Suppose that 

	

	
  (A1) For any ∈^p and 𝒮_(,)|, E(| P) is a linear function of  for any  projection P onto 𝒮_(,)|;
	

	
  (A2) For any ∈^p and 𝒮_(,)|, (| P) is nonrandom for any  projection P onto 𝒮_(,)|.
	


Then column space of  is contained in 𝒮_(,)|.

(ii) Suppose further that 

	

	
  (A3) For any ∈𝒮_(,)|, ≠ 0, the random variable E{[(-)]^2|^o ,^o,δ,δ} is  not equal  to a constant almost surely.
	


Then column space of  is equal to 𝒮_(,)|. 


Conditions (A1) and (A2) are known as linear conditional mean condition and constant conditional variance condition in the the sufficient dimension reduction literature; see <cit.> and <cit.> for more discussions. Condition (A3) is generally considered to be very mild. See <cit.> for more details.

Proposition <ref> suggests that we can utilize  for estimating 𝒮_(,)|. In applications  we use , a discretized version of , to estimate 𝒮_(,)|. We partition the sample space of the uncensored observations with δ=1 into H_1 non-overlapping intervals  I_11,…,I_1H_1, and the sample space of censoring time C with δ=0 into H_0 non-overlapping intervals I_01,…,I_0H_0. Let p_lj=E[I(δ=l,^o ∈ I_lj)] and _ijlm=E[(-)(-)|δ=i,^o ∈ I_ij,δ=l,^o∈ I_lm)], where (i,j,l,m)∈{i,j,l,m: i, l=0 or 1, j=1,…, H_i, m=1,…, H_l}. Then  is expressed as follows

    =∑_ij∑_lmp_ijp_lm(2_p-_ijlm)^2.

We can recover 𝒮_(,)| through the  eigen-decomposition
_i=λ_i _i, where λ_i's are scalars and _i's are p × 1 vectors, and obtain  =(^-1/2_1,…,^-1/2_d).

As ^-1/2 is involved in , , and , the classical sufficient dimension reduction methods fail to work when p>n unless we have a good estimator of  ^-1/2. 

For sufficient variable selection, we do not need the entire matrix  in (<ref>). Proposition <ref> below shows that the following marginal utility of , 

    g_k^*=∑_ij∑_lmp_i,jp_l,m[_k^-1/2(2_p-_ijlm)^-1/2_k]^2,

is a perfect index for sufficient variable selection.  Note that _k^-1/2(2_p-_ijlm)^-1/2_k is not the kth diagonal element of the matrix 2_p-_ijlm
in (<ref>), but the kth diagonal element of the matrix ^-1/2(2_p-_ijlm)^-1/2.


	If conditions  (A1)-(A3) hold, then g_k^*>0 if k∈ and g_k^*=0 if k∉.


The next result gives an alternative expression of g_k^*, which is useful for our derivation.


	Let _lj=E[ I(δ=l,^o ∈ I_lj)] and _lj=E[ I(δ=l,^o ∈ I_lj)].
	Then	
	
    g_k^*=2∑_lj p_lj[_k^-1/2(p_lj^-1_lj-_p)^-1/2_k]^2+4(∑_lj p_lj^-1_k^-1/2_lj_lj^-1/2_k)^2



However, g_k^* in (<ref>) or (<ref>) still involves ^-1/2 which is hard to estimate when p is bigger than or comparable to n. 
We then follow the idea in  independence variable screening (, , and ), i.e., we replace  ^-1/2 in 
(<ref>) by _p and obtain the following  modified directional regression index, 

    g_k=2∑_lj p_lj(V_ljk/p_lj-1)^2+4(∑_lj U_ljk^2/p_lj)^2

where 
U_ljk=E[z_kI(δ=l,^o ∈ I_lj)] and V_ljk=E[z_k^2 I(δ=l,^o ∈ I_lj)], 
z_k=(x_k-μ_k)/σ_k, μ_k=E(x_k), and σ_k^2= Var(x_k).
Although g_k in (<ref>) is not a prefect index for sufficient variable selection as ^-1/2 may be incorrectly treated as _p, it is good enough for sufficient variable screening, i.e., finding a set containing  in (<ref>)
under some conditions. The following result is an example, in which the conditions 
are similar to those in  <cit.> and <cit.>.


	Assume conditions (A1)-(A3). Suppose also that  Cov(x_i,x_j) has the same sign for i,j∈, and that there exists h∈{1,…,d} such that the (j,h)th element of  in (<ref>) 
	have the same sign for all j∈. Then g_k>0 if k∈.
	
	


















§ SURE INDEPENDENCE SCREENING


In this section we show that variable screening by using the index g_k in (<ref>)
holds some asymptotic properties under some conditions. Procedures with weaker conditions are considered in the next section. 
Let (x_ki, t^o_i, δ_i), i=1,...,n, k=1,...,p, be observations from the random sample from (, ^o, δ ),  
μ̂_k = ∑_i=1^n x_k,i/n, σ̂_k={∑_i=1^n(x_k,i-μ̂_k)^2/n}^1/2, 
z_ki = (x_ki- μ̂_k)/σ̂_k,
p_lj=∑_i=1^n I(δ_i=l,t^o_i  ∈ I_lj)/n, U_ljk=∑_i=1^n z_ki I(δ_i=l,t^o_i  ∈ I_lj)/n, and V_ljk=∑_i=1^n z_ki^2  I(δ_i=l,t^o_i  ∈ I_lj)/n. A sample estimator of g_k in (<ref>) is g_k defined by (<ref>) with p_lj, U_ljk, and V_ljk replaced by 
p_lj, U_ljk, and V_ljk, respectively. 
We select the set of covariates such that g_k is large enough. Define

    ={k: g_k ≥γ, 1≤ k ≤ p},

where γ is a threshold to be specified later. 
To study the theoretical property of  in 
(<ref>), we consider the following conditions:

	

	
  (C1) p>n and log p=O(n^ξ) for some ξ∈ (0,1-2κ), where κ is given in condition (C3);
	

	
  (C2) There exist some
	0< ς <1/4 such that  E{exp(tz_k^2)}≤ K_0 for 1≤ k ≤ p and all |t|≤ς, where K_0 is a fixed constant;
	

	
  (C3) min_k∈ g_k>2c_0 n^-κ for some constants c_0>0 and 0 ≤κ≤ 1/2.
	


Condition (C1) was also used by <cit.> and <cit.>, which allows p to be as large as an exponential of the sample size n. Condition (C2) assumes that  all covariates have an exponential-type tails, which is a common technique condition in ultrahigh dimensional data analysis; see, for example, <cit.>.
Condition (C3) is naturally motivated from Proposition <ref>, and requires that the index g_k for k∈ is not too small, which is also a common condition in the literature of sure independence screening <cit.>.

The next theorem confirms the sure screening property of  .


	(i) Assume conditions (C1) and (C2). Then
	
    {max_1≤ k≤ p |g_k-g_k|≥ C_0(log p/n)^1/2}≤ 72p^-τ-1,

	where τ>0 is a constant and C_0 is defined in (<ref>) in the Appendix.

	(ii)		 Additionally, if condition (C3) also holds and γ≤ c_0 n^-κ, then 
	
    {⊆}≥ 1-72p^-τ-1,

	where  is given by (<ref>).


Since g_k  in (<ref>) is a  modified directional regression index and 
Theorem <ref> indicates that the probability in (<ref>) converges to one as n  diverges to infinity,  
we name the proposed covariate screening procedure as the modified directional regression-sure independence screening () method.
Note that ^-1/2= _p is assumed in the derivation of g_k, but it is not needed in establishing the result in Theorem <ref>, as long as (C1)-(C3) hold true. In the next section we obtain some further results in the case where (C3) may be violated.

The threshold value γ depends on constants c_0 and κ in (C3), which is unknown in real applications.
We follow the convention developed in <cit.> 
and define the  screened covariate set as 

    ^*={k :
    g_k≥g_d_n},

where
g_d_n is the d_nth largest ranked index among all
g_k's. Following <cit.>, d_n can be set as ⌊ n/log n⌋, where ⌊ a ⌋ denotes the integer part of a. Theorem <ref> together with Theorem 1 in <cit.> guarantee ( ⊆^* ) converges to  one as n→∞.

Let _ = {_k : k ∈} be the smallest covariate set related to  the life time  , i.e.,  satisfies  |_. 
Sometimes we are interested in identifying   instead of .
For example, if  we assume |, which is typically needed for 
many survival analysis methods although it is not needed  for the asymptotic property of , then  | and  |_
imply |_ so that survival analysis can be carried out using _. However, identifying  may result in a more efficient analysis if information on C | is useful. 

Since ⊆, the sure screening property ( ⊆^* )→ 1 can still be achieved based on Theorem <ref>. 
Unless |_,   is a strict subset of .
Even If we focus on  only, it is unnecessary to do covariate screening to find 
a ^* with ( ⊆^* )→ 1, because both ^* and ^* are screening methods aimed to reduce the size of  covariate set to a manageable number < n and a further dimension reduction or variable selection can be applied to ^* 

as the size of ^* is much smaller than n, 
i.e., d_n /n → 0.








§ ENHANCED SCREENING WITH ITERATION AND RESAMPLING



 §.§ Iterative variable screening

Condition (C3) plays a key role for the sure independence screening property of . 
However, (C3) may be violated since g_k ignores information contained in . The next result identifies a situation where (C3) does not hold.


	Let _1,…,_d be columns of  in (<ref>). For any , if there exists k ∈ such that  ∑_i=1^d|_k_i|=0, then g_k=0 and, hence, (C3) is violated. 


In the situation described by Proposition <ref>, the sure screening property can not be guaranteed. To circumvent this issue, we should handle the correlations among covariates and consider iterative screening. Suppose that we have already selected a   covariate set  _={_k:
k∈}, where  ⊂{ 1,...,p }. 
Define  _=E(_) and _=(_).  For any  e
∉, let _,e=(_,x_e) and  _e|=_e-_,e_^-1_ be the residual of _e regressed on _. Then,  (_e|,_)=0, which suggests that we can adopt  

the marginal utility of modified directional regression based on (_e|,^o,δ) as an index for iterative screening. Define μ_e|=μ_e-_,e_^-1_, σ^2_e|=σ^2_e-_,e_^-1_,e, and z_e|=(_e|-μ_e|)/σ_e| as the standardized version of _e|. Then we define the following iterative modified directional regression index:

    g_e|=2∑_lj p_lj(V_lje|/p_lj-1)^2+4(∑_lj U_lje|^2/p_lj) ^2,      e∉ ,

where 
U_lje|=E[z_e|I(δ=l,^o ∈ I_lj)] and V_lje|=E[z_e|^2 I(δ=l, ^o ∈ I_lj)]. The next proposition illustrates the advantage of the proposed iterative screening method.


	Let  be a nonempty subset of 
	{1,...,p}. Suppose that 
	
		

		
  (C4) min_k∈,i=1,…,d|β_ik|>c_1n^-θ for some constants c_1>0 and 0< θ≤ 1/8, where  β_ik is the (i,k)th element of 
		 in (<ref>);
		

		
  (C5) σ^2_e|≥ c_2 for some constant c_2>0, where e ∉.
		

	
	Then g_e|>2c_0n^-κ for e∈ with some constants c_0>0 and 0 ≤κ≤ 1/2.

Condition (C4) is a mild condition previously used by <cit.>. Condition (C5) means that the eth relevant covariate missed in the previous steps should not be expressed only by the set of covariates selected by previous steps, which is a general condition under iteration construction.   

The result of this proposition illustrates that utilizing the index g_e| is able to identify the informative predictors missed by . To illustrate, suppose that 
^*_1 =^* as define by (<ref>) is selected by . Suppose that we carry out one iteration to obtain 
a covariate set ^*_2 ={e: e ∉^*_1, 
g_e|^*_1≥g_e|^*_1,q}, where g_e|^*_(1),q is the qth largest ranked index among all g_e|^*_1's.
By Proposition <ref>,  ^*_2 may recover some relevant covariates missed by  ^*_1 selected by , with an appropriate choice of q. The covariate set after iteration is 
^*_1∪^*_2. 
Numerical studies show that q can be much smaller than ⌊ n/log n⌋. 

Although ^*_1∪^*_2 is better than 
^*_1 = ^* in terms of containing relevant covariates,  
its size is always larger than the size of ^*. Hence, to apply the iterative variable screening, we do not have to start with  ^*_1 = ^*, especially when we doubt about whether can select all relevant covariates. 
Instead, we may start with  a ^*_1 smaller than 
^* and set the size of final covariate sets selected after iterations
to be the same as that of ^*. This leads to the following general iterative procedure for covariate screening. 

	

	
  * Step 1. Based on  g_k (k=1,…,p), we select p_1 covariates by . Denote the set of indices of selected covariates by ^*_1.
	

	
  * Step 2.  For e ∉^*_1, we estimate g_e|^*_1 by a sample estimator g_e|^*_1. Based on  g_e|^*_1, we select p_2 covariates by  with the resulting covariate set denoted by ^*_2.
	

	
  * Step 3. Repeat Step 2  until the total selected
	number of covariates is d_n. The final selected covariate set is then _I=^*_1∪⋯∪^*_S, where ∑_v=1^S p_v = d_n. 

We name 
this iterative procedure as the modified directional regression-iterative sure independent screening  () method. 
Under conditions (C1)-(C2) and (C4)-(C5), it can be shown similarly to Theorem 3.1 that 
( ⊂_I) → 1 as n →∞. 
Some simulation results are presented in Section 5 for the selection of p_v's and the results show that S=2 works well under our simulated models.








 §.§ Stability Screening

While  is used to improve  in including all relevant covariates, 
the stability selection approach introduced in 
<cit.> 
is designed to  reduce the number of falsely selected covariates 
through combining resampling with high dimensional variable selection.
<cit.>  adapted this resampling mechanism to iterative sure independence screening for genome-wide association studies. Along with their
developments, we further propose the following procedure to improve . 
The algorithm is based on B independent subsamples of size n_s< n without replacement from the training data set. 
For the bth subsample, we apply  to select a candidate covariate set _I^(b).
The stability screened covariate set based on this  procedure is

    _S={k: π_k≥π_0},     π_k=1/B∑_i=1^B I(k∈_I^(i)),

where I( · ) is the indicator function. Following <cit.>, we prespecify threshold value π_0 to be 0.3 or 0.4 in practical use.
We name this procedure as the modified directional regression-stability sure independence screening  () method. In Section 5, we compare  with  in simulations. 




§ NUMERICAL RESULTS

In this section, we assess the performance of the proposed ,  and  by Monte Carlo simulation. We further examine the proposed screening procedure with an empirical analysis of a real-data example. 
 



 §.§ Simulation study
 
The covariate vector  is generated from the multivariate normal distribution with mean 0 and  covariance matrix  whose (i,j)th element is ρ^|i-j| with  ρ= 0, 0.4, or 0.8 throughout our simulations. Let ϵ∼ N(0,1) be an error term  independent of  and the censoring time C.  We consider the following five models representing various types of covariate functions with different degree of nonlinearity, and multiple failure and censoring distributions. 
 
 	
  * M1. =(2_1)^2+12sin(3_2/7)+0.2ϵ, ∼ N(0, 4) - N(5, 1) + N(15, 1), where _1 and  _2 are p× 1 vectors with their first six components being (1,0,1,0,0,0) and (0,0,0,0,1,1), respectively, and rest components being zeros. 
 	
  * M2. =(2_1)^2+|8_2|+0.2ϵ, ∼ N(0, 4)- N(5, 1) + N(30, 1), where _1 and  _2 are same as those in (M1). 
 	
  * M3. =10sin(_1/4)+4|_2|+0.2ϵ, ∼ N(0, 4)- N(5, 1) + N(15, 1), where _1 and  _2 are same as those in (M1).
 	
  * M4. =exp(_1)+|(_2)^3|+0.2ϵ, ∼ N(0, 4)- N(5, 1) + 4N(30, 1), where _1 and  _2 are p× 1 vectors with their first six components being (-4,4,3,0,0,0) and (0,0,1,0,1,0), respectively, and rest components being zeros.
 	
  * M5. =1.5(_1)^2+exp(_2)+0.2ϵ, =_3+8, where _1, _2 and  _3 are p× 1 vectors with their first six components being (1,0,0,0,0,0), (1,2,2,0,0,0) and (0,0,1,0,0,1), and rest components being zeros.
 
 In all models, T  C |. 
 In M1-M3, 4 relevant covariates for T are _1,_3,_5, and _6. In M4, 4 relevant covariates for T are _1,_2,_3, and _5. In M5, 3 relevant covariates for T are _1,_2, and _3, and 2 relevant covariates for C are _3 and _6. 
 
We first fixed the sample size n to be 200 and the dimension p to be 400, and compare our method   with SII <cit.> and QaSIS <cit.>. To evaluate the performance of the 3  methods, we ran 500 simulations and, for each of the 3 methods, we computed the proportion that an individual relevant predictor was selected and the proportion that all relevant predictors were selected. The simulation results reported in 
Table <ref> show that our method is the best among all methods in most cases. And for all the cases in which the other two methods perform well, our method performs at least better than the other two methods. 
 
Table <ref> reports the average computing time of three method with p=200 and different values of n, or n=200 and various values of p. 
The computations are performed using R on ECNU IBM Platform Application Center 9.1.3. We can see that our method is the most computational efficient among the three methods and is increasingly more efficient as n and p are larger. Also, SII is computational intensive, which may lead problems in applications with large n and/or p. 
  

	



	


Next,  we consider n=300 and p= 2000. As SII is very time consuming for  p=2000 in simulation, we only compare  with QaSIS under this setting. To assess the effect of  α in QaSIS, we obtain results for QaSIS with α=05̇ and 07̇. The results reported in Table <ref> show that   overwhelms QaSIS regardless of the choices of α. Moreover, the performance of QaSIS can be influenced by the choice of α.
	