Model-free Feature Screening and FDR Control with Knockoff Features

TL;DR

This paper introduces a model-free, data-adaptive feature screening method for high-dimensional data using projection correlation, with an innovative two-step FDR control approach leveraging knockoff features, demonstrating strong empirical results.

Contribution

It develops a novel model-free screening method based on projection correlation and a two-step FDR control procedure using knockoff features, applicable to complex data scenarios.

Findings

Method achieves sure screening and rank consistency.

FDR control is effective when FDR level ≥ 1/s.

Numerical experiments show superior empirical performance.

Abstract

This paper proposes a model-free and data-adaptive feature screening method for ultra-high dimensional datasets. The proposed method is based on the projection correlation which measures the dependence between two random vectors. This projection correlation based method does not require specifying a regression model and applies to the data in the presence of heavy-tailed errors and multivariate response. It enjoys both sure screening and rank consistency properties under weak assumptions. Further, a two-step approach is proposed to control the false discovery rate (FDR) in feature screening with the help of knockoff features. It can be shown that the proposed two-step approach enjoys both sure screening and FDR control if the pre-specified FDR level $\alpha$ is greater or equal to $1/s$, where $s$ is the number of active features. The superior empirical performance of the proposed methods is justified by various numerical experiments and real data applications.

Tables

Table 1. Table 1: The quantiles of minimum model size for linear and Poisson models in Example 1 over 200 replications. The true model size is 5.

		Model 1.a						Model 1.b
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		5.0	5.0	5.0	5.0	8.0		5.0	5.0	6.0	14.0	125.0
DC-SIS		5.0	5.0	5.0	5.0	6.0		5.0	5.0	10.0	81.2	1305.0
bcDC-SIS		5.0	5.0	5.0	5.0	6.0		5.0	5.0	6.0	20.5	231.4
SIS		5.0	5.0	5.0	5.0	5.0		6.0	238.0	1833.0	3878.5	4915.0
$p=10000$
PC-Screen		5.0	5.0	5.0	5.0	7.0		5.0	5.0	8.0	23.0	233.5
DC-SIS		5.0	5.0	5.0	5.0	6.0		5.0	6.8	21.0	204.2	3511.2
bcDC-SIS		5.0	5.0	5.0	5.0	6.0		5.0	5.0	10.0	43.2	718.8
SIS		5.0	5.0	5.0	5.0	5.0		16.0	703.2	3418.5	7432.0	9651.0
		Model 1.c						Model 1.d
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		5.0	5.0	6.0	8.0	50.9		5.0	5.0	6.0	10.2	139.9
DC-SIS		5.0	8.0	42.5	143.0	722.5		5.0	16.0	54.0	189.0	701.6
bcDC-SIS		5.0	5.0	6.0	14.2	80.4		5.0	5.0	8.0	21.8	156.8
SIS		5.0	39.0	81.5	374.0	3244.4		5.0	45.8	130.5	523.5	3241.0
$p=10000$
PC-Screen		5.0	5.0	6.0	8.0	92.4		5.0	6.0	6.0	13.0	176.2
DC-SIS		5.0	18.2	78.0	277.5	1567.5		5.0	30.8	113.0	410.0	2036.5
bcDC-SIS		5.0	5.0	7.0	19.2	179.1		5.0	5.0	10.0	27.0	412.6
SIS		6.0	58.8	180.5	773.2	4189.0		8.9	73.0	244.5	959.8	5777.7
		Model 1.e						Model 1.f
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		5.0	5.0	5.0	5.0	17.2		5.0	5.0	5.0	5.0	7.0
DC-SIS		90.3	396.0	897.5	1762.5	3787.3		76.3	396.0	893.5	1764.8	3471.8
bcDC-SIS		22.7	82.2	259.5	669.2	2358.7		15.0	87.2	266.5	821.2	2471.5
SIS		178.6	604.8	1137.0	2319.8	4303.4		186.0	606.2	1210.5	2253.0	4261.6
$p=10000$
PC-Screen		5.0	5.0	5.0	6.0	23.0		5.0	5.0	5.0	5.0	12.0
DC-SIS		138.0	788.8	1878.5	3301.0	6831.9		154.2	729.2	1725.5	3424.0	6832.8
bcDC-SIS		45.7	163.5	534.5	1520.5	5415.2		30.0	175.8	509.0	1513.2	5580.0
SIS		462.8	1276.8	2460.0	4164.5	8622.5		512.6	1271.2	2484.5	4281.0	8416.3

Table 2. Table 2: The quantiles of minimum model size for nonlinear models in Example 2 over 200 replications. The true model size is 4.

		Model 2.a						Model 2.b
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		4.0	4.0	4.0	5.2	19.1		4.0	5.0	9.5	26.5	261.9
DC-SIS		600.0	1880.2	2994.5	4052.2	4701.3		4.0	7.8	61.0	541.5	2310.9
bcDC-SIS		488.4	1480.8	2863.0	3967.8	4855.9		4.0	6.0	21.5	88.0	893.8
SIS		709.0	2065.8	3062.5	4160.0	4869.4		54.0	658.5	2692.5	4213.0	4829.1
$p=10000$
PC-Screen		4.0	4.0	4.0	5.0	31.0		4.0	5.8	13.0	48.8	393.9
DC-SIS		664.0	3162.5	5655.5	7605.8	9490.1		4.0	13.8	86.0	843.2	5529.2
bcDC-SIS		663.4	2605.8	5578.5	7745.2	9208.4		4.0	8.0	24.5	150.5	1403.9
SIS		986.2	4048.2	6068.0	8298.5	9752.6		64.5	1193.8	4488.0	8139.5	9725.6
		Model 2.c						Model 2.d
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		4.0	4.0	4.0	6.0	21.3		4.0	5.0	9.5	36.0	169.4
DC-SIS		397.6	1536.8	2750.5	3930.0	4721.8		1639.1	3138.5	3851.5	4434.8	4902.4
bcDC-SIS		196.5	1267.5	2774.0	4236.0	4986.2		605.5	1251.8	2073.5	2972.5	4209.1
SIS		421.2	1615.8	2920.0	4057.0	4761.0		2065.0	3487.8	4083.5	4659.0	4950.3
$p=10000$
PC-Screen		4.0	4.0	4.0	7.2	25.0		4.0	5.0	13.0	62.5	297.4
DC-SIS		668.8	3628.2	6243.5	7920.5	9531.6		3409.6	6380.8	7705.0	8756.5	9790.6
bcDC-SIS		288.2	2006.5	4714.5	7370.2	9869.9		776.4	2567.8	4154.5	5517.8	8251.0
SIS		760.9	3609.5	6155.5	8129.8	9577.1		3333.4	6387.8	8007.5	9252.8	9847.2

Table 3. Table 3: The quantiles of minimum model size for multivariate response models in Example 3 over 200 replications. The true model size is 4.

		Model 3.a						Model 3.b
		5%	25%	50%	75%	95%		5%	25%	50%	75%	95%
$p=5000$
PC-Screen		4.0	4.0	4.0	4.0	6.0		4.0	4.0	6.0	18.0	114.4
DC-SIS		53.0	463.0	1211.0	2349.0	3774.2		641.1	2308.8	3307.0	4257.8	4838.0
bcDC-SIS		24.0	215.5	758.5	1965.0	3999.8		225.2	1270.2	2494.5	3596.8	4709.0
$p=10000$
PC-Screen		4.0	4.0	4.0	4.0	9.0		4.0	4.0	8.0	30.2	302.9
DC-SIS		136.8	978.5	2237.5	4506.0	8106.6		1804.7	4510.5	6445.5	8153.5	9707.6
bcDC-SIS		83.0	546.0	1828.5	4264.5	7711.3		444.4	2217.0	4695.0	7122.5	9076.4

Table 4. Table 4: FDR control for linear and Poisson models in Example 4. The true model size is 10.

$\alpha$	$\left|\widehat{𝒜}\right|$	$X_1$	$X_2$	$X_3$	$X_4$	$X_5$	$X_6$	$X_7$	$X_8$	$X_9$	$X_{10}$	All	$\widehat{\text{FDR}}$
	Model 4.a
0.10	11.245	0.995	0.995	1.000	0.995	0.995	0.995	0.995	0.995	0.995	0.990	0.990	0.097
0.15	11.935	0.995	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.990	0.130
0.20	12.860	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.188
0.25	14.095	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.244
0.30	15.270	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.285
	Model 4.b
0.10	10.545	0.945	0.945	0.960	0.950	0.955	0.950	0.955	0.955	0.955	0.945	0.915	0.079
0.15	11.350	0.985	0.990	0.995	0.995	0.995	1.000	0.995	1.000	0.990	0.970	0.920	0.100
0.20	12.460	0.990	0.995	1.000	0.995	0.995	1.000	0.995	1.000	1.000	0.985	0.955	0.164
0.25	13.570	0.995	1.000	1.000	0.995	0.995	1.000	0.995	1.000	1.000	0.985	0.965	0.216
0.30	14.705	1.000	1.000	1.000	0.995	0.995	1.000	1.000	1.000	1.000	0.990	0.980	0.260
	Model 4.c
0.10	11.165	0.995	0.995	0.995	1.000	0.995	0.995	0.995	0.995	0.995	0.995	0.995	0.092
0.15	11.660	1.000	0.995	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.116
0.20	13.075	1.000	0.995	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.193
0.25	14.420	1.000	0.995	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.254
0.30	15.570	1.000	0.995	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.995	0.995	0.299
	Model 4.d
0.10	10.630	0.940	0.945	0.940	0.945	0.955	0.940	0.950	0.940	0.950	0.940	0.925	0.092
0.15	11.660	0.985	0.980	0.995	0.985	0.995	0.990	1.000	1.000	0.990	0.990	0.930	0.122
0.20	12.870	0.990	0.990	1.000	1.000	1.000	0.995	1.000	1.000	0.990	0.990	0.965	0.193
0.25	13.830	0.990	0.995	1.000	1.000	1.000	0.995	1.000	1.000	0.995	0.995	0.975	0.242
0.30	15.020	0.995	0.995	1.000	1.000	1.000	0.995	1.000	1.000	0.995	0.995	0.980	0.288
	Model 4.e
0.10	9.790	0.855	0.870	0.860	0.860	0.865	0.870	0.875	0.870	0.865	0.865	0.815	0.086
0.15	11.590	0.955	0.985	0.975	0.965	0.980	0.995	0.980	0.970	0.965	0.970	0.825	0.123
0.20	12.895	0.970	0.990	0.995	0.995	1.000	1.000	0.995	0.990	0.985	0.980	0.935	0.189
0.25	14.025	0.980	0.995	0.995	0.995	1.000	1.000	0.995	0.990	0.990	0.990	0.935	0.244
0.30	15.040	0.985	1.000	0.995	0.995	1.000	1.000	0.995	0.990	0.990	0.990	0.945	0.283

Table 5. Table 5: Comparison between 𝒜 ^ 1 subscript ^ 𝒜 1 \widehat{\cal A}_{1} and 𝒜 ^ ​ ( T α ) ^ 𝒜 subscript 𝑇 𝛼 \widehat{\cal A}(T_{\alpha}) . The number active features is s = 10 𝑠 10 s=10 and pre-specified FDR level for PC-Knockoff is α = 0.2 . 𝛼 0.2 \alpha=0.2.

	PC-Screen				PC-Knockoff
	$\left|{\widehat{𝒜}}_1\right|$	All	$\widehat{\mathrm{FDR}}$		$\left|\widehat{𝒜}\left(T_{\alpha }\right)\right|$	All	$\widehat{\mathrm{FDR}}$
Model 4.a	100	0.995	0.90		12.9	0.995	0.188
Model 4.b	100	0.995	0.90		12.5	0.955	0.164
Model 4.c	100	0.995	0.90		13.0	0.995	0.193
Model 4.d	100	0.995	0.90		12.9	0.965	0.193
Model 4.e	100	0.995	0.90		12.9	0.935	0.189

Table 6. Table 6: Sample mean and sample standard deviation of R 2 superscript 𝑅 2 R^{2} for training set and test set over 200 bootstrap replications for supermarket data.

	Training set			Test set
	Mean	SD		Mean	SD
PC-Knockoff	0.8681	0.0042		0.8608	0.0281
DC-SIS with ${\chi }^2$-test	0.8534	0.0043		0.8488	0.0276

Equations

$$H_0:\mathbf{x}\text{ and }\mathbf{y}\text{are independent}\text{versus}{H}_1:\text{ otherwise.}$$

$$\mathrm{Pcov}(\mathbf{x},\mathbf{y}{)}^2$$

$$=\iiint {\mathrm{cov}}^2\{I({\mathbf{\alpha }}^{\mathrm{T}}\mathbf{x}\le u),I({\mathbf{\beta }}^{\mathrm{T}}\mathbf{y}\le v)\}d{F}_{U,V}(u,v)d\mathbf{\alpha }d\mathbf{\beta },$$

$$\mathrm{PC}(\mathbf{x},\mathbf{y}{)}^2=\frac{\mathrm{Pcov}(\mathbf{x},\mathbf{y}{)}^2}{\mathrm{Pcov}(\mathbf{x},\mathbf{x})\mathrm{Pcov}(\mathbf{y},\mathbf{y})},$$

$$\mathrm{Pcov}(\mathbf{x},\mathbf{y}{)}^2=$$

$$=$$

$$+$$

$$-2$$

$$\widehat{\mathrm{Pcov}}(\mathbf{X},\mathbf{Y}{)}^2$$

$$\widehat{\mathrm{Pcov}}(\mathbf{X},\mathbf{X}{)}^2$$

$$a_{klr}$$

$${\bar{a}}_{k\cdot r}$$

$$A_{klr}$$

$$b_{klr}$$

$${\bar{b}}_{k\cdot r}$$

$$B_{klr}$$

$$\widehat{\mathrm{PC}}(\mathbf{X},\mathbf{Y}{)}^2=\frac{\widehat{\mathrm{Pcov}}(\mathbf{X},\mathbf{Y}{)}^2}{\widehat{\mathrm{Pcov}}(\mathbf{X},\mathbf{X})\widehat{\mathrm{Pcov}}(\mathbf{Y},\mathbf{Y})}.$$

$$\mathrm{Pr}\left\{|\widehat{\mathrm{Pcov}}(\mathbf{X},\mathbf{Y})^{2}-\mathrm{Pcov}(\mathbf{x},\mathbf{y})^{2}|>\varepsilon\right\}\leq c_{1}\exp\{-c_{2}n\varepsilon^{2}\},\ \text{and}\$$

$$\mathrm{Pr}\left\{|\widehat{\mathrm{PC}}(\mathbf{X},\mathbf{Y}{)}^2-\mathrm{PC}(\mathbf{x},\mathbf{y}{)}^2|>\epsilon \right\}\le 5c_1\exp \{-c_2\sigma n{\epsilon }^2\},$$

$$\mathcal{A}=\{k:F(\mathbf{y}|\mathbf{x})\text{functionally depends on}{X}_k,k=1,\dots ,p\}.$$

$$\hat{\mathcal{A}}(\delta )=\{k:\widehat{\mathrm{PC}}({\mathbf{X}}_k,\mathbf{Y}{)}^2\ge \delta ,1\le k\le p\},$$

$$\mathrm{Pr}\Big{(}{\cal A}\subseteq\widehat{{\cal A}}(\delta)\Big{)}\geq 1-O\Big{(}s\exp\{-c_{4}n^{1-2\kappa}\}\Big{)},$$

$$\mathrm{Pr}\Big{(}{\cal A}\subseteq\widehat{{\cal A}}(c_{3}n^{-\kappa})\Big{)}\geq 1-O(s\exp\{-c_{4}n^{1-2\kappa}\}).$$

$$\mathrm{Pr}\Big{(}{\cal A}\subseteq\widehat{{\cal A}}(c_{3}n^{-\kappa})\Big{)}\geq 1-O(s\exp\{-c_{4}^{\prime}n^{1-2(\kappa+\eta)}\}+n\exp\{-c_{4}^{\prime\prime}n^{\eta}\}),$$

$$\mathrm{Pr}\left(\underset{k\in \mathcal{A}}{\min }{\hat{\omega }}_k-\underset{k\in {\mathcal{A}}^c}{\max }{\hat{\omega }}_k>0\right)>1-O(p\exp \{-c_5{n}^{1-2\kappa }\}),$$

$${\lim \inf }_{n\to \infty }\left(\underset{k\in \mathcal{A}}{\min }{\hat{\omega }}_k-\underset{k\in {\mathcal{A}}^c}{\max }{\hat{\omega }}_k\right)>0,\text{almost surely}.$$

$$\mathrm{cov}(\mathbf{x},\tilde{\mathbf{x}})=\mathbf{G},\text{where}\mathbf{G}=\left[\begin{array}{cc}\mathbf{\Sigma }& \mathbf{\Sigma }-\mathrm{diag}\{\mathbf{h}\}\\ \mathbf{\Sigma }-\mathrm{diag}\{\mathbf{h}\}& \mathbf{\Sigma }\end{array}\right],$$

$$h_j=2{\lambda }_{\min }(\mathbf{\Sigma })\vee 1\text{for}j=1,\dots ,p,$$

$$\sum _j|1-h_j|,$$

$$h_j\ge 0,\text{diag}\{\mathbf{h}\}⪯2\mathbf{\Sigma }.$$

$$W_j=\mathrm{PC}(X_j,\mathbf{y}){}^2-\mathrm{PC}({\tilde{X}}_j,\mathbf{y}){}^2,j=1,\dots ,p,$$

$${\hat{W}}_j=\widehat{\mathrm{PC}}({\mathbf{X}}_j,\mathbf{Y}){}^2-\widehat{\mathrm{PC}}({\tilde{\mathbf{X}}}_j,\mathbf{Y}){}^2,j=1,\dots ,p.$$

$$\mathrm{FDP}(t)=\frac{\#\{j\in {\mathcal{A}}^c:{\hat{W}}_j\ge t\}}{\#\{j:{\hat{W}}_j\ge t\}},$$

$$\#\{j\in {\mathcal{A}}^c:{\hat{W}}_j\ge t\}\approx \#\{j\in {\mathcal{A}}^c:{\hat{W}}_j\le -t\}\le \#\{j:{\hat{W}}_j\le -t\},$$

$$\widehat{\mathrm{FDP}}(t)=\frac{\#\{j:{\hat{W}}_j\le -t\}}{\#\{j:{\hat{W}}_j\ge t\}}.$$

$$T_{\alpha }=\min \left\{t\in \mathcal{W}:\frac{1+\#\{j:{\hat{W}}_j\le -t\}}{\#\{j:{\hat{W}}_j\ge t\}}\le \alpha \right\},$$

Full text

Keywords: Data adaptive, multivariate response model, nonlinear model, projection correlation, rank consistency, sample splitting, sure screening

1 Introduction

The technological development has made extensive data collection and storage feasible in diverse fields. Datasets with ultra-high dimensional features characterize many contemporary research problems in machine learning, computer science, statistics, engineering, social science, finance and so on. When the features contain redundant or noisy information, estimating their functional relationship with the response may become quite challenging in terms of computational expediency, statistical accuracy, and algorithmic stability (Fan et al., 2009; Hall & Miller, 2009; Lv & Liu, 2014). To overcome such challenges caused by ultra-high dimensionality, Fan & Lv (2008) proposed a sure independence screening (SIS) method, which aims to screen out the redundant features by ranking their marginal Pearson correlations. The SIS method is named after the sure independence screening property, which states the selected subset of features contains all the active ones with probability approaching one. The promising numerical performance soon made SIS popular among ultra-high dimensional studies (Liu et al., 2015; Liu & Li, 2020). The sure screening idea has been applied to many important statistical problems including generalized linear model (Fan & Song, 2010), multi-index semi-parametric model (Zhu et al., 2011), nonparametric model (Fan et al., 2011; Liu et al., 2014), quantile regression (He et al., 2013; Wu & Yin, 2015) and compressed sensing (Xue & Zou, 2011) among others.

In addition to the sure screening property, we argue an appealing screening method should satisfy the following two properties. First, the screening method should be model-free in the sense that it can be implemented without specifying a regression model. In the ultra-high dimensional regime, it is challenging, if not impossible, to specify a correct regression model before removing a huge number of redundant features. Hence, the model-free property is desired as it guarantees the effectiveness of the screening method in the presence of model misspecification. The model-free screening method becomes a hot research topic in recent years (Zhu et al., 2011; Li et al., 2012; Mai & Zou, 2015). The second property is data-adaptive which means the screening method should not be sensitive to assumptions like independence, sub-Gaussianity, and univariate response. Such assumptions are usually not satisfied under ultra-high dimensional settings. Even met on the population level, they can be violated in the observed sample due to ultra-high dimensionality. Therefore, the screening methods that are sensitive to such assumptions may perform poorly in real applications. The data-adaptive screening methods also draw a certain amount of attention recently. For instance, He et al. (2013), Wu & Yin (2015) and Ma et al. (2017), among others, considered quantile-based screening which adapts to heavy-tailed data. In addition, Wang (2012) and Fan, Ke & Wang (2020) developed screening methods for strongly correlated features.

Unfortunately, none of the aforementioned screening methods enjoys sure screening, model-free and data-adaptive properties simultaneously. For example, the SIS is tailored to the linear regression and depends on the independence assumption. Li et al. (2012) developed a model-free sure independence screening procedure based on the distance correlation. However, its sure screening property requires sub-exponential assumption for features and response. The Kolmogorov distance based screening method proposed in Mai & Zou (2012) is robust against heavy-tailed data, but it only works for binary classification problems. Pan et al. (2016) proposed a pairwise sure screening procedure for linear discriminant analysis, which is not sensitive to the tail behavior of the features, but requires balanced categories.

In this paper, we propose a model-free and data-adaptive feature screening method named PC-Screen. The PC-Screen is based on ranking the projection correlations between features and response variables. The projection correlation, proposed by Zhu et al. (2017), is a measure of dependence between two random vectors which enjoys several nice probability properties. The PC-Screen does not require specifying any regression model and is insensitive to the correlations and moment conditions of the dataset. As the projection correlation is dimension-free to both random vectors, the PC-Screen can be applied to multi-task learning problems (Caruana, 1997). For instance, we can find a parsimonious set of features that are jointly dependent on the multivariate response. As existing asymptotic results are not suitable for the “large $p$ small $n$ problem” studied in this paper, we also develop non-asymptotic concentration-type inequalities for empirical projection correlation. Further, we demonstrate the proposed PC-Screen enjoys not only the sure screening property but also a stronger result called rank consistency property. The only condition required is a minimum signal strength gap between active and inactive features. The extensive numerical examples show the proposed method wins the horse racing against its competitors in various scenarios.

Most feature screening methods depend on some threshold parameter that controls the cut-off between active and inactive features. The optimal choice of the threshold typically depends on unknown parameters. Under a specific model assumption, the threshold can be selected by cross-validation or information criteria approaches. However, when no model assumption is imposed, such parameter selection approaches are not applicable as the loss function which measures the goodness of fit is not well defined. In addition, existing screening methods tend to sacrifice the false discovery rate for sure screening property by choosing a conservative threshold parameter, leading to an inflated model size. Recently, Zhao & Li (2012) studied the sure independence screening for Cox models with a principled method to select the threshold aiming at controlling the false positive rate. Song et al. (2014) proposed a censored rank independence screening for survival data and chose the threshold by estimating the proportion of active features. The validity of this procedure relies on the independence assumption of multiple test statistics which may easily get violated in the model-free setup. In this paper, we tackle the issue of threshold selection with a two-step procedure named PC-Knockoff. In the first step, we apply the PC-Screen method to obtain an over-fitted subset of moderate size from the ultra-high dimensional features. In the second step, we construct knockoff counterparts for the features which survive in the first step. Conditioning on the sure screening of the first step, we further select a parsimonious model with FDR controlled below a pre-specified level with a statistic that utilizing the knockoff features. Theoretical analysis shows that when the pre-specified level is not too aggressive, the PC-Knockoff procedure enjoys sure screening property, as well as conditional FDR control simultaneously with high probability. We also validate the theoretical findings with various numerical examples.

The rest of the paper is organized as follows. In Section 2, we briefly review the definition and properties of projection correlation and demonstrate its non-asymptotic properties. Then we propose the PC-Screen procedure and show its sure screening and rank consistency properties under very mild conditions. Section 3 studies the PC-Knockoff procedure which selects the threshold with FDR control as well as sure screening property. Section 4 assesses the finite sample performance of the PC-Screen and PC-Knockoff methods with several simulated examples and a real data application. We briefly summarize the paper in Section 5. Due to the limited space, we provide the proofs of some theoretical results in the Appendix, and relegate the proofs of remaining results as well as some additional numerical examples to the a supplementary file.

2 Model-free and data-adaptive screening procedure

2.1 Projection correlation

To pave the way for the proposed screening procedure, we first provide some background on the projection correlation and its properties introduced in Zhu et al. (2017). Let $\mathbf{x}\in\mathbb{R}^{p}$ and $\mathbf{y}\in\mathbb{R}^{q}$ be two random vectors. The projection correlation is elicited by the following independence testing problem,

[TABLE]

The null hypothesis holds if and only if $U=\bm{\alpha}^{\mathrm{\scriptstyle T}}\mathbf{x}$ and $V=\bm{\beta}^{\mathrm{\scriptstyle T}}\mathbf{y}$ are independent for all unit vectors $\bm{\alpha}$ and $\bm{\beta}$. Let $F_{U,V}(u,v)$ be the joint distribution of $(U,\ V)$, and $F_{U}(u)$ and $F_{V}(v)$ be the marginal distributions of $U$ and $V$. The squared projection covariance is defined as

[TABLE]

where $I(\cdot)$ is the indicator function. Furthermore, the projection correlation between $\mathbf{x}$ and $\mathbf{y}$ is defined as the square root of

[TABLE]

and we follow the convention $0/0=0$.

In general $0\leq\mathrm{PC}(\mathbf{x},\mathbf{y})\leq 1$, testing whether $\mathbf{x}$ and $\mathbf{y}$ are independent amounts to testing whether $\mathrm{PC}(\mathbf{x},\mathbf{y})=0$. The projection correlation is a measure of dependence between two random vectors and enjoys some appealing properties. Let $\mathbf{x}$ and $\mathbf{y}$ be two random vectors with continuous marginal and joint probability distributions, $\mathrm{PC}(\mathbf{x},\mathbf{y})=0$ if and only if $\mathbf{x}$ and $\mathbf{y}$ are independent. We remark here this property does not hold in general without the assumption that $(\mathbf{x},\mathbf{y})$ is jointly continuous. When $\mathbf{x}$ and $\mathbf{y}$ are two dependent discrete random variables that are constructed in a similar fashion as in Hoeffding (1948), it is possible that we have $\mathrm{PC}(\mathbf{x},\mathbf{y})=0$.

Zhu et al. (2017) gives an explicit formula for the squared projection covariance in (2.1). Let $(\mathbf{x}_{1},\mathbf{y}_{1}),\ \ldots,\ (\mathbf{x}_{5},\mathbf{y}_{5})$ be 5 independent random copies of $(\mathbf{x},\mathbf{y})$, then

[TABLE]

where $\|\cdot\|$ is the $L_{2}$ norm. Equation (2.3) shows that the projection covariance only depends on the vectors through forms $(\mathbf{x}_{k}-\mathbf{x}_{l})/\|\mathbf{x}_{k}-\mathbf{x}_{l}\|$ and $(\mathbf{y}_{k}-\mathbf{y}_{l})/\|\mathbf{y}_{k}-\mathbf{y}_{l}\|$ whose second moments are unity. This gives us the intuition that the projection covariance is free of the moment conditions on $(\mathbf{x},\mathbf{y})$ which are usually required by some other measurements, such as distance correlation (Li et al., 2012).

Let $\mathbf{X}=(\mathbf{x}_{1},\ \ldots\ ,\mathbf{x}_{n})^{\mathrm{\scriptstyle T}}$ and $\mathbf{Y}=(\mathbf{y}_{1},\ \ldots\ ,\mathbf{y}_{n})^{\mathrm{\scriptstyle T}}$ be an observed sample of $(\mathbf{x},\mathbf{y})$. Equation (2.3) leads to a straightforward estimator of $\mathrm{Pcov}(\mathbf{x},\mathbf{y})^{2}$ based on a $U$-statistic, yet it is difficult to calculate (Székely & Rizzo, 2010). An equivalent form of the $U$-statistic is given in Zhu et al. (2017). In particular, the squared sample projection variance and covariance between $\mathbf{X}$ and $\mathbf{Y}$ can be calculated as

[TABLE]

where $k,l,r=1,\dots,n,$

[TABLE]

Then the sample projection correlation between $\mathbf{X}$ and $\mathbf{Y}$ is defined as the square root of

[TABLE]

Based on (2.4), the sample projection correlation can be computed in $O(n^{3})$.

We first provide exponential-type deviation inequalities for squared sample projection covariance and correlation.

Theorem 1**.**

For any $0<\varepsilon<1$ satisfying $n\geq 10\pi^{2}/\varepsilon$, there exists positive constants $c_{1}$ and $c_{2}$, such that

[TABLE]

[TABLE]

where $\sigma=\min\{\sigma_{x}^{3}\sigma_{y}^{3}/4M^{4},\sigma_{x}^{2}\sigma_{y}^{2}/4M^{4}\}$, $\sigma_{x}=\mathrm{Pcov}(\mathbf{x},\mathbf{x})^{2}$, $\sigma_{y}=\mathrm{Pcov}(\mathbf{y},\mathbf{y})^{2}$ and $M=4\pi^{2}$.

The proof of Theorem 1 is based on an exponential-type deviation inequality for $U$-statistic and can be found in the supplementary material. The above exponential inequalities do not depend on the dimensionality and moment conditions of both random vectors. The exception probability decays exponentially with sample size $n$ which guarantees good finite sample performance of the proposed estimator.

2.2 PC-Screen procedure

In this subsection, we propose a model-free and data-adaptive screening procedure utilizing the nice properties of projection correlation. Let $\mathbf{y}=(Y_{1},\ldots,Y_{q})^{\mathrm{\scriptstyle T}}$ be the response vector of $q$ variables and $\mathbf{x}=(X_{1},\ldots,X_{p})^{\mathrm{\scriptstyle T}}$ be the vector of $p$ features. To avoid the trivial discussion, we restrict ourselves to the non-degenerate case. That is, $\min_{1\leq k\leq p}\mathrm{Pcov}(X_{k},X_{k})^{2}\geq\sigma_{0}^{2}$ and $\min_{1\leq k\leq q}\mathrm{Pcov}(Y_{k},Y_{k})^{2}\geq\sigma_{0}^{2}$ for some $\sigma_{0}>0$.

Denote $F(\mathbf{y}|\mathbf{x})$ the conditional distribution function of $\mathbf{y}$ given $\mathbf{x}$. Without specifying any regression model of $\mathbf{y}$ on $\mathbf{x}$, we define the index set of active features by

[TABLE]

The number of active features is $s=|{\cal A}|$, the cardinality of ${\cal A}$. The features that do not belong to ${\cal A}$ are inactive features. We use ${\cal A}^{c}$, the complement of ${\cal A}$, to denote the index set of inactive features. The above setting abstracts a large number of sparse regression problems including linear model, generalized linear model, additive model, semi-parametric model, non-linear model and so on. Moreover, it also allows the multivariate response and grouped predictors.

Suppose that $\{(\mathbf{x}_{i},\mathbf{y}_{i}),\ i=1,\ldots,n\}$ is a random sample of $(\mathbf{x},\mathbf{y})$. Denote $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\mathrm{\scriptstyle T}}$ and $\mathbf{Y}=(\mathbf{y}_{1},\ldots,\mathbf{y}_{n})^{\mathrm{\scriptstyle T}}$. In the ultra-high dimensional regime, it is natural to assume that the number of features $p$ greatly exceeds the sample size $n$, but the number of active features $s$ is smaller than $n$. For a given feature $X_{k}$, $k=1,\ldots,p$, a sufficient condition for $X_{k}$ to be an inactive feature is the independence between $X_{k}$ and $\mathbf{y}$. This, together with Theorem 1, motivates us to screen out the features whose projection correlations with $\mathbf{y}$ are close to [math]. As a result, we estimate the set of active features by

[TABLE]

where $\delta$ is a pre-specified positive threshold, and $\mathbf{X}_{k}$ is the $k$th column of $\mathbf{X}$. With a proper choice of $\delta$, we show that the proposed feature screening procedure enjoys the sure screening property, which states that with probability approaching 1, all active features are included in $\widehat{\cal A}(\delta)$. We name this feature screening procedure as projection correlation based screening, or PC-Screen.

Denote $\omega_{k}\equiv\mathrm{PC}(X_{k},\mathbf{y})^{2}$ and $\widehat{\omega}_{k}\equiv\widehat{\mathrm{PC}}(\mathbf{X}_{k},\mathbf{Y})^{2}$ to be the squared population and sample projection correlations between the $k$th feature and response, respectively. To analyze the property of PC-Screen, we impose the following minimum signal strength condition.

Condition 1. (Minimum signal strength)

• (a)

For some $c_{3}>0$ and $0\leq\kappa<1/2$, $\min_{k\in{\cal A}}\omega_{k}\geq 2c_{3}n^{-\kappa}.$

• (b)

For some $c_{3}>0$ and $0\leq\kappa<1/2$, $\min_{k\in{\cal A}}\omega_{k}-\max_{k\in{\cal A}^{c}}\omega_{k}\geq 2c_{3}n^{-\kappa}.$

Remark 1**.**

Condition 1 (a) is a minimum signal strength condition that assumes the squared projection correlations between active features and response are uniformly bounded below, and cannot converge to zero too fast as $n$ diverges. Condition 1 (b) imposes an assumption on the gap of signal strength between active and inactive features. Condition 1 (a) is weaker than Condition 1 (b), since $\omega_{k}$ is always non-negative. Such minimum signal strength condition can be viewed as a sparsity assumption that guarantees active features to be distinguished from the inactive ones. Condition 1 is a very mild condition as we allow the minimum signal strength to converge to zero as the sample size diverges. Empirically, one may verify Condition 1 by plotting the sorted empirical projection correlations of features in descending order and visually identifying whether there exists an “elbow” shape, or whether the sequence can be segmented into two groups. A more rigorous verification approach is to consider a multiple-testing procedure. For each feature, one can calculate the projection correlation based test statistic proposed in Zhu et al. (2017). Since it is not the main focus of this paper, we do not pursue more details.

The following two theorems state the sure screening property and the rank consistency property of PC-Screen.

Theorem 2** (Sure screening).**

Under Condition 1 (a), take $\delta\leq\min_{k\in{\cal A}}\omega_{k}/2$, we have

[TABLE]

where $c_{4}$ is a positive constant.

In Theorem 2, if we set $\delta=c_{3}n^{-\kappa}$, which satisfies the condition $\delta\leq\min_{k\in{\cal A}}\omega_{k}/2$, we have

[TABLE]

From (2.7), we know that if $\delta=c_{3}n^{-\kappa}$, all active features are selected with probability approaching 1 as $n\rightarrow\infty$. In fact, any choice of $\delta\leq c_{3}n^{-\kappa}$ leads to the sure screening property. With the same choice of $\delta$, Li et al. (2012) showed that the distance correlation based screening method (DC-SIS) satisfies

[TABLE]

where $c_{4}^{\prime}$, $c_{4}^{\prime\prime}$ and $\eta$ are positive constants. Thus, PC-Screen achieves a faster rate than the DC-SIS since (1) we do not have the term $n\exp\{-c_{4}^{\prime\prime}n^{\eta}\}$ and (2) we do not have an extra $\eta$ in the power of the first term. The faster rate of PC-Screen is due to the fact that projection correlation is not sensitive to dimensionality and free of moment conditions.

Theorem 3** (Rank consistency).**

Under Condition 1 (b), we have

[TABLE]

where $c_{5}$ is some positive constant. If $\log p=o(n^{1-2\kappa})$ with $0\leq\kappa<1/2$, then we have

[TABLE]

The rank consistency in Theorem 3 is a stronger result than the sure screening property. When the signal strength gap between active and inactive features satisfies Condition 1 (b), the active features are always ranked ahead of the inactive ones with high probability. In other words, there exists a choice of $\delta$ on the solution path that can perfectly separate the active and inactive sets with high probability.

3 Screening with FDR control

3.1 Motivation

In the PC-Screen procedure, the threshold $\delta$ controls the cut-off between active and inactive features. Theorem 2 suggests choosing $\delta=cn^{-\kappa}$ for some positive constants $c$ and $\kappa<1/2$. With certain model assumptions, the threshold $\delta$ (or equivalently $c$ and $\kappa$) can be selected by cross-validation or information criterion approaches. However, in the model-free setup, such approaches are not directly applicable as the loss functions that measure the goodness of fit are not well defined. More recently, Zhao & Li (2012) and Song et al. (2014) suggested to select the threshold by controlling the false positive rate and estimating the proportion of active features $s/p$, respectively. These two methods are tailored for survival analysis and are not directly applicable to model-free setup.

In practice, one can arbitrarily select a conservative threshold to ensure that all active features are included with high probability. However, it will include too many inactive features and inflate the FDR. Therefore, selecting the threshold parameter for model-free screening methods will inevitably cause the issue of balancing the trade-off between the sure screening property and FDR control. In this section, we propose a two-step procedure to address this issue by utilizing knockoff features. The construction of knockoff features does not depend on the regression model, and hence is suitable for the model-free setting. Also, the second-order knockoff features can be easily computed by an equicorrelated construction or solving a semidefinite program (Barber & Candès, 2015). The proposed procedure enjoys sure screening property and controls FDR simultaneously with high probability.

3.2 Knockoff features: a brief review

Recently, knockoff has drawn huge attention due to its success in variable selection and many important applications such as genome-wide association studies. The concept of knockoff was first proposed in Barber & Candès (2015) for the fixed design matrix problem and then extended to the random design matrix setting known as Model-X knockoffs (Candès et al., 2018). For a more detailed development of knockoff, see Barber & Candès (2015); Candès et al. (2018); Fan, Demirkaya, Li & Lv (2020) and references therein. In this subsection, we briefly review the notations and definitions of knockoffs for further discussions.

Let $\mathbf{y}=(Y_{1},\ldots,Y_{q})^{\mathrm{\scriptstyle T}}$ be a response vector of $q$ variables and $\mathbf{x}=(X_{1},\ldots,X_{p})^{\mathrm{\scriptstyle T}}$ be a covariate of $p$ features. We say $\widetilde{\mathbf{x}}=(\widetilde{X}_{1},\ldots,\widetilde{X}_{p})$ is a knockoff copy of $\mathbf{x}$ if it satisfies the following conditions.

Condition 2. (Exact knockoff features)

• (a)

Swap $X_{j}$ with its knockoff counterpart $\widetilde{X}_{j}$ does not change the joint distribution of $(\mathbf{x},\widetilde{\mathbf{x}})$ for $j=1,\ldots,p$;

• (b)

Given $\mathbf{x}$, $\widetilde{\mathbf{x}}$ is independent of $\mathbf{y}$, i.e., \widetilde{\mathbf{x}}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\mathbf{y}|\mathbf{x}.

Remark 2**.**

Condition 2 (a) requires the original and knockoff features to be pairwise exchangeable. Condition 2 (b) indicates that knockoff features are conditionally independent of response variables; this is trivially satisfied if $\widetilde{\mathbf{x}}$ is generated without using the information of $\mathbf{y}$.

Constructing knockoff features that exactly follow Condition 2 is challenging, especially when the dimensionality of features $p$ is large. In general, generating exact knockoff features requires the knowledge of the underlying distribution of $\mathbf{x}$, which is usually not available in practice. Barber & Candès (2015) studied the variable selection problem with knockoffs for fixed design matrix $\mathbf{X}\in\mathbb{R}^{n\times p}$. The construction of exact knockoffs is not necessary if assuming the response $\mathbf{y}$ follows a linear regression model with Gaussian error and $n\geq 2p$. A more recent study (Candès et al., 2018) proposed to construct exact knockoff features in a manner of generating sequential conditional independent pairs under the assumption that the distribution of $\mathbf{x}$ is known.

Without the knowledge of the distribution of $\mathbf{x}$, one can construct approximate second-order knockoff features such that $(\mathbf{x},\widetilde{\mathbf{x}})$ is pairwise exchangeable with respect to the first two moments. In other words, the mean vector and covariance matrix of $(\mathbf{x},\widetilde{}\mathbf{x})$ is invariant if we swap $X_{j}$ and $\widetilde{X}_{j}$ for any $j=1,\ldots,p$. The invariant of mean is trivial, and can be achieved by forcing $\text{E}(\mathbf{x})=\text{E}(\widetilde{\mathbf{x}})$. Suppose $\mathrm{cov}(\mathbf{x})=\bm{\Sigma}$, the second-order pairwise exchangeable condition is equivalent to

[TABLE]

where $\mathbf{h}=(h_{1},\ldots,h_{p})^{\mathrm{\scriptstyle T}}$ is a vector that makes $\mathbf{G}$ a positive semidefinite covariance matrix.

Barber & Candès (2015) introduced two approaches to construct the second-order knockoffs. The first approach is known as the equicorrelated construction, which sets

[TABLE]

where $\lambda_{\min}(\bm{\Sigma})$ is the minimum eigenvalue of $\bm{\Sigma}$, and $a\vee b$ denotes the larger one between $a$ and $b$. The second method, named semidefinite programme, finds $h_{j}$ by solving a semidefinite program of the following form

[TABLE]

However, both methods are not directly applicable in the high-dimensional scenarios since both (3.2) and (3.3) require $2p<n$.

Remark 3**.**

When $(\mathbf{x},\widetilde{\mathbf{x}})$ is jointly Gaussian, the equivalence of the first two moments implies the equivalence of the joint distribution and hence (3.1) constructs exact knockoff features. However, when the Gaussian assumption does not hold, the accuracy of the second-order approximation depends on the impact of ignoring higher-order moments of $\mathbf{x}$. Also, the covariance matrix $\bm{\Sigma}$ is usually unknown and needs to be estimated. Hence the estimation accuracy of $\bm{\Sigma}$ also affects the validation of the second-order approximation of exact knockoff features. In addition to sample covariance estimator, more sophisticated estimators of $\bm{\Sigma}$ can be obtained under additional structure or moment conditions, see Fan et al. (2013) and Ke et al. (2019), among others. As the above two issues are not of key interest in this paper, we do not pursue further in these directions.

3.3 FDR control with knockoff features

Suppose $\widetilde{\mathbf{x}}=(\widetilde{X}_{1},\ldots,\widetilde{X}_{p})^{\mathrm{\scriptstyle T}}$ is a knockoff copy of $\mathbf{x}=(X_{1},\ldots,X_{p})^{\mathrm{\scriptstyle T}}$, we propose to measure the population level dependence between $X_{j}$ and response $\mathbf{y}$ by the following quantity

[TABLE]

where ${\mathrm{PC}}(X_{j},\mathbf{y})$ and ${\mathrm{PC}}(\widetilde{X}_{j},\mathbf{y})$ are the projection correlations as defined in (2.2). When $\widetilde{X}_{j}$ is an exact knockoff feature of $X_{j}$, $W_{j}$ will be a non-negative quantity. Further, $W_{j}>0$ implies the distribution of $\mathbf{y}$ depends on $X_{j}$ and $W_{j}=0$ if $\mathbf{y}$ is independent with $X_{j}$ conditioning on active features.

Given a random sample $\{\mathbf{Y},\mathbf{X},\widetilde{\mathbf{X}}\}$ drawn from $\{\mathbf{y},\mathbf{x},\widetilde{\mathbf{x}}\}$, we estimate ${W}_{j}$ by

[TABLE]

where $\widehat{\mathrm{PC}}(\mathbf{X}_{j},\mathbf{Y})$ and $\widehat{\mathrm{PC}}(\widetilde{\mathbf{X}}_{j},\mathbf{Y})$ are sample projection correlations defined in (2.5). Intuitively, a large positive value of $\widehat{W}_{j}$ provides some evidence that the distribution of $\mathbf{y}$ depends on $X_{j}$. On the other hand, if $X_{j}$ is an inactive feature, $|\widehat{W}_{j}|$ is likely be be small and $\widehat{W}_{j}$ is equally likely to be positive or negative. This intuition is justified by the following lemma.

Lemma 1**.**

Let $\widetilde{}\mathbf{x}$ be an exact knockoff copy of $\mathbf{x}$ and ${\cal A}^{c}=\{j_{1},\dots,j_{r}\}$. Then

1. (i)

$W_{j_{k}}=0$ for all $j_{k}\in{\cal A}^{c}$. 2. (ii)

Conditioning on $|\widehat{}\mathbf{w}|=(|\widehat{W}_{1}|,\dots,|\widehat{W}_{p}|)^{\mathrm{\scriptstyle T}}$, $I_{j_{1}},\dots,I_{j_{r}}$ follow i.i.d. $Bernoulli(0.5)$, where $I_{j_{k}}=1$ if $\widehat{W}_{j_{k}}>0$ and 0 otherwise.

For a fixed threshold $t>0$, the false discovery proportion (FDP) is defined as

[TABLE]

where $\#\{\cdot\}$ is the cardinality of a set and we follow the convention that $0/0=0$. The false discovery rate (FDR) is defined as the expectation of FDP, i.e., $\mathrm{FDR}(t)=\text{E}\left[{\mathrm{FDP}}(t)\right]$. According to Lemma 1, $\widehat{W}_{j}$ is equally likely to be positive or negative if $X_{j}$ is an inactive feature. Therefore, we have

[TABLE]

which leads to a conservative estimation of ${\mathrm{FDP}}(t)$,

[TABLE]

To control FDR at a pre-specified level $\alpha$, we follow the knockoff+ procedure (Barber & Candès, 2015) to choose the threshold $T_{\alpha}$ as

[TABLE]

where ${\cal W}=\{|\widehat{W}_{j}|:1\leq j\leq p\}/\{0\}$. The extra term $1$ in the numerator makes the choice of $T_{\alpha}$ slightly more conservative. Then, the active set is selected as

[TABLE]

3.4 PC-Knockoff procedure

The knockoff feature construction methods discussed in Section 3.2 require $2p<n$ and hence are not applicable to high-dimensional scenarios. To address this issue, we propose a two-step procedure, named PC-Knockoff, to screen active features and control the FDR. To avoid mathematical challenges caused by the reuse of sample, we follow the simple sample splitting idea, which has been widely used in statistics and recent examples include hypothesis testing (e.g. Fan et al., 2019), error variance estimation (e.g. Chen et al., 2018), variable selection (e.g. Barber & Candès, 2019), large scale inference (e.g. Fan, Demirkaya, Li & Lv, 2020), and so on. We partition the full sample into two disjoint subsamples with sample sizes $n_{1}$ and $n_{2}=n-n_{1}$. More specifically, let $\mathbf{X}^{(1)}\in\mathbb{R}^{n_{1}\times p}$ and $\mathbf{X}^{(2)}\in\mathbb{R}^{n_{2}\times p}$ be a random partition of $\mathbf{X}$, and let $\mathbf{Y}$ follow the same partition. Without loss of generality, we write

[TABLE]

The two steps of PC-Knockoff procedure are introduced as follows:

(1) Screening step: We rank all $p$ features in descending order based on the sample projection correlation $\widehat{\mathrm{PC}}(\mathbf{X}_{j}^{(1)},\mathbf{Y}^{(1)})$. Then, we select the top $d$ features such that $2d<n_{2}$. Denote the set of selected $d$ features by $\widehat{\cal A}_{1}$. In practice, one can set $d$ to be a relatively large value as long as it satisfies $2d<n_{2}$.

(2) Knockoff step: Let

[TABLE]

We construct knockoff features for $\mathbf{X}^{(2)}_{\widehat{\cal A}_{1}}$ by either the equicorrelated construction as in (3.2) or the semidefinite programme as in (3.3). Denoted $\widetilde{}\mathbf{X}^{(2)}_{\widehat{\cal A}_{1}}$ the constructed knockoff features for $\mathbf{X}^{(2)}_{\widehat{\cal A}_{1}}$. Then, we calculate

[TABLE]

where $\mathbf{X}^{(2)}_{\widehat{\cal A}_{1},j}$ and $\widetilde{\mathbf{X}}^{(2)}_{\widehat{\cal A}_{1},j}$ are the $j$th columns of $\mathbf{X}^{(2)}_{\widehat{\cal A}_{1}}$ and $\widetilde{}\mathbf{X}^{(2)}_{\widehat{\cal A}_{1}}$, respectively. For a pre-specified FDR level $\alpha$, we use (3.6) to choose the threshold $T_{\alpha}$ and the set of selected active features is given by

[TABLE]

The proposed two-step approach has two advantages. First, we apply the PC-Screen procedure to reduce the number of features from $p$ to $d$, which allows us to construct second-order knockoff features. Second, by ruling out $p-d$ inactive features in the screening step, we reduce the total computation cost from $O(n^{3}p)$ to $O(n_{1}^{3}p+n_{2}^{3}d)$. We summarize the PC-Knockoff procedure in Algorithm 1. In fact, the Algorithm 1 provides a general framework for feature screening with FDR control. One can easily modify Algorithm 1 by replacing PC with other measurement statistics such as Pearson correlation, distance correlation, etc.

Remark 4**.**

We want to add a note that Algorithm 1 is not tuning free as one needs to specify the subsample size $n_{1}$ and target dimension $d$ for the screening step. However, Algorithm 1 is not sensitive to the choice of these two hyper-parameters as Theorem 2 guarantees the sure screening property of the screening step under mild conditions. In practice, we suggest small $n_{1}$. As a result, we leave a relatively large subsample for the knockoff step. A relatively large $n_{2}$ allows more features to be selected in the screening step (larger $d$) and more accurate second-order knockoff features constructed in the knockoff step.

Denote ${\cal E}$ the event that the sure screening property is satisfied in the screening step, i.e.,

[TABLE]

The sure screening property in Theorem 2 ensures that, with a relatively large choice of $d$, the event ${\cal E}$ holds with high probability. To be specific, let $\widehat{\omega}^{(1)}_{(1)}\geq\widehat{\omega}^{(1)}_{(2)}\geq\dots\geq\widehat{\omega}^{(1)}_{(p)}$ be the order statistics of sample squared projection correlations based on $(\mathbf{X}^{(1)},\mathbf{Y}^{(1)})$. If Condition 1 (a) holds and $\widehat{\omega}^{(1)}_{(d)}\leq c_{3}n_{1}^{-\kappa}$, then event ${\cal E}$ holds with probability at least $1-O(s\exp\{c_{4}n_{1}^{1-2\kappa}\})$. Conditioning on $\cal E$, the following theorem states that the PC-Knockoff procedure can control the FDR of selected features under the pre-specified level of $\alpha$.

Theorem 4**.**

Let $\widetilde{}\mathbf{X}_{\widehat{\cal A}_{1}}$ be a knockoff copy of $\mathbf{X}_{\widehat{\cal A}_{1}}$ satisfying Condition 2. For any $\alpha\in[0,1]$, the set of selected features $\widehat{\cal A}(T_{\alpha})$ given by Algorithm 1 satisfies

[TABLE]

Theorem 4 states that, with exact knockoff features and conditioning on $\cal E$, the PC-Knockoff procedure can control the FDR under the pre-specified level $\alpha\in[0,1]$. As $\cal E$ occurs with probability close to 1, the $\mathrm{FDR}$ can be controlled even without conditioning on $\cal E$ (Barber & Candès, 2019). We refer to Barber & Candès (2019) and Fan, Demirkaya, Li & Lv (2020) for more discussion regarding $\mathrm{FDR}$ control in two-step procedures.

We hope the PC-Knockoff procedure can maintain sure screening property and control FDR simultaneously. This task is challenging as the procedure needs to balance the trade-off between type I and type II errors. To guarantee the sure screening property, the procedure is likely to select an over-fitted model which leads to higher type I errors as well as FDR. On the other hand, the FDR control forces a parsimonious but possibly under-fitted model that can increase type II errors and ruin the sure screening property. In the following, we study the conditions under which this challenging task can be achieved by the PC-Knockoff procedure. Conditioning on ${\cal E}$, the following theorem states that the simultaneous achievement of sure screening property and FDR control under level $\alpha$, which depends on the relationship between $s$ and $\alpha$.

Theorem 5**.**

Under the conditions of Theorem 4 and further assume $\min_{k\in{\cal A}}W_{j}\geq 4c_{3}n_{2}^{-\kappa}$ for some $c_{3}>0$ and $0\leq\kappa<1/2$.

1. (i)

If $\alpha\geq 1/s$, we have $\text{Pr}({\cal A}\subseteq\widehat{\cal A}(T_{\alpha})|{\cal E})\geq 1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$. 2. (ii)

If $\alpha<1/s$, we have $\text{Pr}({\cal A}\subseteq\widehat{\cal A}(T_{\alpha})\cup\widehat{\cal A}(T_{\alpha})=\varnothing|{\cal E})\geq 1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$. Furthermore, if $s>2$, we have $\text{Pr}({\cal A}\subseteq\widehat{\cal A}(T_{\alpha})|{\cal E})\leq C(s)+O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$, where $0<C(s)<1$ is a constant that only depends on $s$.

The part (i) of Theorem 5 together with Theorem 4 suggest that if $\alpha$ is chosen to be greater or equal to $1/s$, the PC-Knockoff procedure enjoys the sure screening property and controls FDR under $\alpha$ with high probability. When $\alpha$ is chosen to be smaller than $1/s$, we either recover the active set or end up with an empty set with high probability. The probability of recovering the active set is upper bounded by some constant $C(s)$ depending on $s$. Therefore there is no guarantee that PC-Knockoff can select all active features while controlling FDR under $\alpha$. As we know, the value of $\alpha$ controls the amount of type I errors that we can tolerate. One can imagine that the smaller the $\alpha$ is, the more challenging the task is to achieve sure screening and FDR control simultaneously. It is because we allow fewer and fewer inactive features to be selected. There is a phase transition that happens when $\alpha$ goes below $1/s$. In order to satisfy the sure screening property, the procedure will fail to control the FDR at $\alpha$ with non-negligible probability if any inactive feature is selected. In other words, the procedure has to exactly recover the true active set to satisfy FDR control and sure screening simultaneously. We numerically validate this phase transition phenomenon through a simulated example in the supplementary file.

Theorem 5 discourages us to pursue a too aggressive $\alpha$ in practice as the PC-Knockoff procedure may lose the sure screening property. The phase transition between part (i) and part (ii) can also be used as a rule of thumb guideline to estimate $s$ in practice. For a sequence of grid points of $\alpha$ in $(0,1)$, we find the largest grid point $\alpha^{*}$ such that the PC-Knockoff procedure selects an empty set. Then, we can roughly estimate $\widehat{s}$ as the integer part of $1/\alpha^{*}$.

Remark 5**.**

As a two-step procedure, the power of the PC-Knockoff is conditioned on the event that the sure screening property is satisfied in the screening step, i.e., ${\cal E}$. This together with the results in Theorem 5 yield that the probability that the PC-Knockoff procedure does not lose any power is at least $1-O(s\exp\{-c_{4}n_{1}^{1-2\kappa}\})-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$. Empirically, one may follow the “data recycling ” idea proposed in Barber & Candès (2019) to improve the power of the PC-Knockoff procedure.

As a byproduct of the PC-Knockoff procedure, the statistic $\widehat{W}_{j}$ defined in (3.5) can also be used to measure the marginal dependence between the response and the $j$th feature. Given a threshold $\delta$, we can screen a model based on $\widehat{W}_{j}$ by considering the following set

[TABLE]

The next theorem states that the screening procedure based on $\widehat{W}_{j}$ also enjoys the sure screening and rank consistency properties.

Theorem 6**.**

Suppose $\widetilde{}\mathbf{X}$ is an exact knockoff copy for $\mathbf{X}$ and $\min_{k\in{\cal A}}W_{j}\geq 2c_{3}n^{-\kappa}$ for some positive constants $c_{3}>0$ and $0<\kappa<1/2$, then

1. (i)

\mathrm{Pr}\Big{(}{\cal A}\subseteq\widehat{\cal A}_{W}(\delta)\Big{)}\geq 1-O\Big{(}s\exp\{-c_{4}n^{1-2\kappa}\}\Big{)} for $\delta\leq c_{3}n^{-\kappa}$. 2. (ii)

$\mathrm{Pr}\left(\underset{k\in{\cal A}}{\min}\ \widehat{W}_{k}-\underset{k\in{\cal A}^{c}}{\max}\ \widehat{W}_{k}>0\right)\geq 1-O(p\exp\{-c_{4}n^{1-2\kappa}\})$.

Recall that the rank consistency property in Theorem 3 requires a minimum signal gap between active and inactive sets, that is $\min_{k\in{\cal A}}\omega_{k}-\max_{k\in{\cal A}^{c}}{\omega_{k}}>2c_{3}n^{-\kappa}$. However, the rank consistency result in Theorem 6 only requires a minimum signal strength of active features and no condition is imposed on the inactive set. Due to the construction of $W_{j}$, signals of inactive features are canceled out by their knockoff counterparts. As a result, the rank consistency property holds even when some inactive features are spuriously correlated with the response. If we can construct high-quality knockoff features, the screening procedure based on $\widehat{W}_{j}$ can be more powerful than PC-Screen.

4 Numerical examples

4.1 Screening performance

In this subsection, we use simulated examples to assess the finite sample performance of the proposed projection correlation based feature screening procedure (PC-Screen) and compare it with sure independence screening (Fan & Lv, 2008, SIS), distance correlation based screening (Li et al., 2012, DC-SIS) and bias-corrected distance correlation based screening (Székely & Rizzo, 2014, bcDC-SIS). Within each replication, we rank the features in descending order by the above four screening criteria and record the minimum model size that contains all active features. The screening performance is measured by the $5\%,25\%,50\%,75\%$ and $95\%$ quantiles of the minimum model size over 200 replications. Throughout this subsection, we denote $\bm{\Sigma}=(\sigma_{ij})_{p\times p}$ with $\sigma_{ij}=0.5^{|i-j|}$. To mimic ultra-high dimensional scenario, we set $n=100$ and $p=5,000,10,000$ for each example.

Example 1: Linear and Poisson models

Consider the linear model $Y=\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\beta}+\varepsilon$ with $\bm{\beta}=(\mathbf{1}_{5}^{\mathrm{\scriptstyle T}},\mathbf{0}_{p-5}^{\mathrm{\scriptstyle T}})^{\mathrm{\scriptstyle T}}$. We generate covariates $\mathbf{x}$ and $\varepsilon$ independently from the following 4 scenarios.

**  Model 1.a:**

$\mathbf{x}\sim N(\mathbf{0},\mathbf{\Sigma})$ and $\varepsilon\sim N(0,1)$.

**  Model 1.b:**

$\mathbf{x}\sim N(\mathbf{0},\mathbf{\Sigma})$ and $\varepsilon\sim\text{Cauchy}(0,1)$.

**  Model 1.c:**

$\mathbf{u}\sim\text{Cauchy}(\mathbf{0},\mathbf{I}_{p})$, $\mathbf{x}=\mathbf{\Sigma}^{1/2}\mathbf{u}$ and $\varepsilon\sim N(0,1)$.

**  Model 1.d:**

$\mathbf{u}\sim\text{Cauchy}(\mathbf{0},\mathbf{I}_{p})$, $\mathbf{x}=\mathbf{\Sigma}^{1/2}\mathbf{u}$ and $\varepsilon\sim\text{Cauchy}(0,1)$.

In above models, $\text{Cauchy}(\mathbf{0},\mathbf{I}_{p})$ stands for the $p$-dimensional standard Cauchy distribution which is heavy-tailed. Hence, in Models 1.b – 1.d, at least one of $\mathbf{x}$ and $\varepsilon$ is heavy-tailed. We also consider the following two Poisson regression models

**  Model 1.e:**

(Continuous) $Y=\exp\{\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\beta}\}+\varepsilon$, where $\varepsilon\sim N(0,1)$.

**  Model 1.f:**

(Discrete) $Y\sim\text{Poisson}(\exp\{\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\beta}\})$.

Let $\bm{\beta}=(\mathbf{2}_{5}^{\mathrm{\scriptstyle T}},\mathbf{0}_{p-5}^{\mathrm{\scriptstyle T}})^{\mathrm{\scriptstyle T}}$ and we draw $\mathbf{x}$ from $N(\mathbf{0},\mathbf{\Sigma})$. Model 1.e is the Poisson regression model with continuous response while the response in Model 1.f is discrete.

The quantiles of the minimum model size that includes all 5 active features are summarized in Table 1. In the linear benchmark Model 1.a, all four competitors perform well. For Models 1.b – 1.d, the SIS completely fails at the presence of heavy-tailed features and errors. Both distance correlation-based methods struggle to maintain a reasonable model size at $75\%$ and $95\%$ quantiles. In contrast, PC-Screen works reasonably well in all scenarios and outperforms the other three methods by big margins. Similarly, for the Poisson models 1.e and 1.f, PC-Screen can recover the true active set with a model size close to 5 while the other three methods can perform as bad as random guesses at $75\%$ and $95\%$ quantiles.

Example 2: Nonlinear models

Consider the following four non-linear data generating models

**  Model 2.a:**

$Y=5X_{1}+2\sin(\pi X_{2}/2)+2X_{3}\mathbf{1}\{X_{3}>0\}+2\exp\{5X_{4}\}+\varepsilon.$

**  Model 2.b:**

$Y=3X_{1}+3X_{2}^{3}+3X_{3}^{-1}+5\mathbf{1}\{X_{4}>0\}+\varepsilon.$

**  Model 2.c:**

$Y=1-5(X_{2}+X_{3})^{3}\exp\{-5(X_{1}+X_{4}^{2})\}+\varepsilon.$

**  Model 2.d:**

$Y=1-5(X_{2}+X_{3})^{-3}\exp\{1+10\sin(\pi X_{1}/2)+5X_{4}\}+\varepsilon.$

Models 2.a and 2.b admit additive structure while Models 2.c and 2.d have more challenging nonlinear structures. In addition, we generate $\mathbf{x}\sim N(\mathbf{0},\mathbf{\Sigma})$ and $\varepsilon\sim N(0,1)$. Hence, for each model above, the active set contains the first 4 covariates in $\mathbf{x}$.

The quantiles of the minimum model size that includes all 4 active features are presented in Table 2. Again, we observe that PC-Screen significantly outperforms the other three methods in all scenarios. For Model 2.a, the $50\%$ quantile of the minimum model size of PC-Screen is exactly four while the other three methods need much larger model sizes to recover the active set. For Model 2.b, DC-SIS and bcDC-SIS preform comparably to PC-Screen at $5\%$ and $25\%$ quantiles but much worse for other higher quantiles. For Models 2.c and 2.d, PC-Screen performs reasonably well while the other methods fail to effectively screen out the inactive features. The 95% quantiles of SIS, DC-SIS, and bcDC-SIS are almost as large as $p$. This indicates, in the worst-case scenario, SIS, DC-SIS, and bcDC-SIS are hopeless to effectively reduce the dimensionality without missing any active feature.

Example 3: Multivariate response models

In this experiment, we investigate the performance of PC-Screen for multivariate response models. We omit SIS method in this experiment as it is not applicable to multivariate response problems. We generate $\mathbf{y}=(Y_{1},Y_{2})^{\mathrm{\scriptstyle T}}$ from a bivariate normal distribution with conditional mean $\bm{\mu}_{Y|X}=(\mu_{1}(\mathbf{x}),\mu_{2}(\mathbf{x}))^{\mathrm{\scriptstyle T}}$ and covariance matrix $\bm{\Sigma}_{Y|X}=(\sigma_{ij})_{2\times 2}$, where $\sigma_{11}=\sigma_{22}=1$ and $\sigma_{12}=\sigma_{21}=\sigma(\mathbf{x})$. Following the setting in Li et al. (2012), we set $\bm{\theta}=(\mathbf{2}^{\mathrm{\scriptstyle T}}_{4},\mathbf{0}^{\mathrm{\scriptstyle T}}_{p-4})^{\mathrm{\scriptstyle T}}$ and generate $\mu_{1}(\mathbf{x})$, $\mu_{2}(\mathbf{x})$ and $\sigma(\mathbf{x})$ from the two models below.

**  Model 3.a:**

$\mu_{1}(\mathbf{x})=\exp\{2(X_{1}+X_{2})\}$, $\mu_{2}(\mathbf{x})=X_{3}+X_{4}$ and $\sigma(\mathbf{x})=\sin(\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\theta})$.

**  Model 3.b:**

$\mu_{1}(\mathbf{x})=2\sin(\pi X_{1}/2)+X_{3}+\exp\{1+X_{4}\}$, $\mu_{2}(\mathbf{x})=X_{1}^{-2}+X_{2}$ and

$\sigma(\mathbf{x})=(\exp\{\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\theta}\}-1)/(\exp\{\mathbf{x}^{\mathrm{\scriptstyle T}}\bm{\theta}\}+1)$.

The union of the active sets of $\mu_{1}(\mathbf{x})$, $\mu_{2}(\mathbf{x})$ and $\sigma_{1}(\mathbf{x})$ contains the first four covariates in $\mathbf{x}$. The simulation results are summarized in Table 3. Again, the PC-Screen method performs strikingly well compared to the other two methods.

4.2 FDR control performance

In this subsection, we use simulated examples to numerically assess the FDR control as well as sure screening property of the proposed PC-Knockoff procedure. We refer to Algorithm 1 for implementation details.

Example 4: FDR control for linear and Poisson models

Consider the following five regression models with ten active variables.

**  Model 4.a:**

Same as Model 1.a except that $\bm{\beta}=(\mathbf{1}_{s},\mathbf{0}_{p-s})$ with $s=10$.

**  Model 4.b:**

Same as Model 4.a except that $\varepsilon\sim t_{2}$, the $t$ distribution with degrees of freedom 2.

**  Model 4.c:**

Same as Model 4.a except that $\mathbf{x}=0.9\mathbf{x}_{1}+0.1\mathbf{x}_{2}$, where $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ are independently drawn from $\mathbf{x}_{1}\sim N(\mathbf{0},\bm{\Sigma})$ and $\mathbf{x}_{2}\sim t_{2}(\mathbf{0},\bm{\Sigma})$, respectively.

**  Model 4.d:**

Same as Model 1.e except that $\bm{\beta}=2\cdot(\mathbf{1}_{s},\mathbf{0}_{p-s})$ with $s=10$.

**  Model 4.e:**

Same as Model 1.f except that $\bm{\beta}=2\cdot(\mathbf{1}_{s},\mathbf{0}_{p-s})$ with $s=10$.

In this example, we set $n=1000$, $p=5000$ and repeat 200 replications for each scenario. In each replication, we randomly divide the sample into two non-overlapping subsamples. The sample size and target dimension in the screening step are set to be $n_{1}=250$ and $d=100$, and the sample size used to construct knockoff features is $n_{2}=750$. In addition, the covariance matrix is set to be $\bm{\Sigma}=(\sigma_{ij})$ with $\sigma_{ij}=0.5^{|i-j|}$. The performance of FDR control is examined under a sequence of specified levels: $\alpha=0.1,0.15,0.20,0.25$ and $0.30$.

We summarize the results in Table 4 in which $\alpha$ is the pre-specified FDR level, $|\widehat{\cal A}|$ is the average number of selected variables, the column ‘$X_{j}$’ represents the probability that the active variable $X_{j}$ is selected, the column ‘All’ represents the sure screening probability (i.e., the probability that all active variables are selected) and the column ‘$\widehat{\text{FDR}}$’ is the empirical FDR (i.e., the average of empirical FDP). According to Table 4, the proposed PC-Knockoff procedure controls the empirical FDR under the pre-specified level $\alpha$ for almost all scenarios. The only exception is in Model 4.c we have $\widehat{\text{FDR}}=0.254$ at $\alpha=0.25$. In Model 4.c, $\mathbf{x}$ is drawn from a mixture of multivariate normal and multivariate $t_{2}$ distributions. As a result, the second-order knockoffs may not approximate the exact ones very well. For the other four models, the second-order knockoffs perform well as features are normally distributed. Besides FDR control, the PC-Knockoff procedure maintains the sure screening property reasonably well in all scenarios. For Models 4.a – 4.d, the probabilities of selecting all active variables are ranging from $91.5\%$ to $99.5\%$ regardless of the pre-specified level $\alpha$. Even for the more challenging scenario Model 4.e, the probability of selecting all active variables simultaneously is greater than $93.5\%$ when $\alpha\geq 0.2$. Please notice that this is achieved under very small model sizes (i.e., small $|\widehat{\cal A}|$), thanks to the knockoff features.

To illustrate how the two steps in PC-Knockoff work together, we also compare the active set $\widehat{\cal A}_{1}$ selected by the screening step without knockoff features (i.e. selected by PC-Screen) and the active set $\widehat{\cal A}(T_{\alpha})$ selected by PC-Knockoff. Table 5 summarizes the average number of selected features $|\widehat{\cal A}|$, the sure screening probability ‘All’ and empirical false discovery rate $\widehat{\mathrm{FDR}}$ for Models 4.a – 4.e. For the screening step, we conservatively select active features by including the top $d=100$ features. The results show that $\widehat{\cal A}_{1}$ enjoys the sure screening property but has a very high empirical FDR ($\widehat{\mathrm{FDR}}=0.90$). On the contrary, the active set $\widehat{\cal A}(T_{\alpha})$ selected by the PC-Knockoff is able to control the FDR below the pre-specified level with only a little sacrifice of power.

4.3 Supermarket data

In this subsection, we apply the PC-Knockoff procedure to study a supermarket dataset (Wang, 2009; Chen et al., 2018). The dataset consists of $n=464$ observations of daily records from a supermarket. The response $Y$ is the number of customers visited the supermarket on that day. The covariates are sale volumes of $p=6398$ products. Due to data privacy, the detailed product codes are not released in the dataset. Instead, we name the covariates by their indices, i.e., $X_{1},\ldots,X_{6398}$. Both the response and predictors have been standardized to have zero mean and unit variance. The goal is to screen a parsimonious set of products whose sale volumes significantly contribute to the daily number of customers. Meanwhile, we want to control the FDR at level $\alpha=0.1$.

The implementation of PC-Knockoff follows Algorithm 1. We set $n_{1}=200$ and $d=50$ in the screening step. That is, we randomly choose $200$ observations and use PC-Screen to pre-screen 50 features in the screening step. Then we construct second-order knockoffs for the pre-screened 50 features using the remaining 264 observations in the knockoff step. Under the pre-specified FDR level $\alpha=0.1$, the PC-Knockoff procedure selects a model of $12$ variables: $X_{3},X_{6},X_{10},X_{11},X_{30},X_{42},X_{48},X_{71},X_{129},X_{139},X_{176}$, and $X_{400}$. For each selected variable, we draw a scatter plot between this variable and the response. We observe obvious outliers in the scatter plots of $X_{11},X_{71}$ and $X_{400}$. We report this observation in Figure 1, where red triangles denote potential outliers and blue dots denote the rest data points. Besides, red dashed curves and blue solid curves are the fitted local polynomial regression curves with and without potential outliers, receptively. The gray shaded areas are the $95\%$ confidence regions. By comparing blue and red curves, we find the existence of outliers visually alters the fitted curves. This justifies the PC-Knockoff procedure is insensitive to the presence of outliers. The scatter plots between the response and other nine variables ($X_{3}$, $X_{6}$, $X_{10}$, $X_{30}$, $X_{42}$, $X_{48}$, $X_{129}$, $X_{139}$ and $X_{176}$) are presented in Figure 2, which includes the local polynomial regression curves with the $95\%$ confidence regions. Figure 2 indicates that PC-Knockoff can detect various functional relationships.

This supermarket dataset has also been analyzed by Chen et al. (2018). They first apply DC-SIS to screen variables and then fit an additive model with selected variables. They further employ the Wald’s $\chi^{2}$-test with the refitted cross-validation error variance estimate to determine if the selected features are significant at a pre-specified level. As a result, Chen et al. (2018) selected seven significant variables $X_{3},X_{6},X_{11},X_{39},X_{42},X_{62}$ and $X_{139}$. Next, we compare the in-sample fitting and out-of-sample prediction performance between the model selected by PC-Knockoff and the one in Chen et al. (2018) through a bootstrap experiment. In each replication, we randomly split the dataset into a training set of size 400 and a test set of size 64. We fit two additive models with the features selected by the two competitive methods, respectively. Then we calculate and record the training and testing $R^{2}$s for the two models. We repeat it for 200 replications. The sample mean and sample standard deviation of $R^{2}$s are reported in Table 6, which show that the model selected by PC-Knockoff yields higher sample means of $R^{2}$ for both training set and test set.

Due to the limited space, we present an additional real application to microarray data in Section S.2 of the supplementary file.

5 Conclusion

We study the feature screening problems under a general setup. The proposed PC-Screen method developed in Section 2 is novel in the sense that it does not impose any regression model assumption and is robust against possibly heavy-tailed data. The response variable is also allowed to be multidimensional. Regarding the theoretical aspect, the asymptotic results developed in Zhu et al. (2017) do not apply to the “large $p$ small $n$ problem”, and hence we develop non-asymptotic results for the empirical projection correlation. The theoretical analysis shows that the PC-Screen method satisfies sure screening and rank consistency properties. Numerically, we show PC-Screen outperforms popular competitors over various data generative processes. Moreover, this paper tackles the false discovery rate control problem in feature screening. In Section 3, we propose a two-step procedure named PC-Knockoff to determine the threshold of active features. With the sample splitting idea, we first screen the ultra-high dimensional features to a moderate model size and then further select a model with FDR control using a statistic built upon knockoff features. The PC-Knockoff procedure controls FDR under a pre-specified level while maintaining sure screening property with high probability. In practice, the PC-Knockoff procedure works well in various scenarios.

Appendix: Proof of results in Section 3

In the appendix, we provide the proofs of Theorem 4 and Theorem 5 in Section 3. We relegate the proofs of remaining results as well as some additional numerical examples to a supplementary file.

A.1 Proof the Theorem 4

Denote $\widehat{\cal A}_{1}$, the set of features that are selected in the screening step. Throughout this proof, we restrict ourselves to the event ${\cal E}=\{{{\cal A}}\subset\widehat{\cal A}_{1}\}$. According to Theorem 2, the probability of $\cal E$ is at least $1-O(s\exp\{-c_{4}n_{1}^{1-2\kappa}\})$. To simplify the notations, we omit the condition on $\cal E$.

Denote by ${\cal B}={\cal A}\cap\widehat{\cal A}_{1}$ and ${\cal B}^{c}={\cal A}^{c}\cap\widehat{\cal A}_{1}$ the intersections between active and inactive sets with $\widehat{\cal A}_{1}$, respectively. Notice that any feature with $|\widehat{W}_{j}|=0$ will not be included in the final selected set $\widehat{\cal A}(T_{\alpha})$. Hence, without loss of generality, we assume $\widehat{\cal A}_{1}=\{1,2,\dots,d\}$ and $|\widehat{W}_{1}|\geq|\widehat{W}_{2}|\geq\cdots\geq|\widehat{W}_{d}|>0$. For ease of presentation, we further define $\widehat{W}_{d+1}=0$. Then we have

[TABLE]

The first inequality holds since $\{j:j\in{\cal B}^{c}\ \text{and}\ \widehat{W}_{j}\leq-T_{\alpha}\}\subset\{j:\in\widehat{\cal A}_{1}\ \text{and}\ \widehat{W}_{j}\leq-T_{\alpha}\}$ and the second inequality is due to the definition of $T_{\alpha}$ in (3.6). In order to find $T_{\alpha}$, one can simply try different values of $t$ starting from the smallest value $t=|\widehat{W}_{d+1}|=0$, then move to the second smallest value $t=|\widehat{W}_{d}|$, then move to $t=|\widehat{W}_{d-1}|$, and so on. The procedure stops as soon as it finds a value of $t$ satisfying (3.6). In this process, $T_{\alpha}$ can be regarded as a stopping time. More rigorously, for $k=d+1,d,d-1,\dots,1$, we define

[TABLE]

where $V_{+}(k)=\#\{j:j\in{\cal B}^{c},j\leq k,\widehat{W}_{j}>0\}$ and $V_{-}(k)=\#\{j:j\in{\cal B}^{c},j\leq k,\widehat{W}_{j}\leq 0\}$. Let ${\cal F}_{k}$ be the $\sigma$-algebra generated by $\{V_{\pm}(d+1),V_{\pm}(d),\dots,V_{\pm}(k),Z_{d+1},Z_{d}\dots,Z_{k}\}$ where $Z_{d+1}=0$ and

[TABLE]

As a result, given ${\cal F}_{k}$, we know whether $k$ is in the active set ${\cal B}$ or not.

Next we show that the process $M(d+1),M(d),\dots,M(1)$ is a super-martingale running backward with respect to the filtration ${\cal F}_{d+1}\subset{\cal F}_{d}\subset\dots\subset{\cal F}_{1}$. On one hand, if $k\in{\cal B}$, we have $V_{+}(k)=V_{+}(k-1)$, $V_{-}(k)=V_{-}(k-1)$ and thus $M(k)=M(k-1)$. On the other hand, if $k\in{\cal B}^{c}$, we have

[TABLE]

where $I_{k}=I\{W_{k}>0\}$ and $I\{\cdot\}$ is the indicator function. Let $d_{0}=\#\{j:j\in{\cal B}^{c}\}$ with the inactive set ${\cal B}^{c}=\{j_{1},j_{2},\dots,j_{d_{0}}\}$. From Lemma 1 part (ii), we know that $I_{j_{1}},I_{j_{2}}\dots,I_{j_{d_{0}}}$ are i.i.d. Bernoulli(0.5) random variables. Thus conditioning on ${\cal F}_{k}$, (hence $V_{+}(k)$, $V_{-}(k)$ are known), we have

[TABLE]

Thus in the case $k\in{\cal B}^{c}$,

[TABLE]

Therefore, $\text{E}[M(k-1)|{\cal F}_{k}]\leq M(k)$, implying that $M(k),k=d+1,\dots,1$ is a super-martingale with respect to $\{{\cal F}_{k}\}$. By definition, $T_{\alpha}$ is a stopping time with respect to the backward filtration $\{{\cal F}_{k}\}$. According to the optional stopping theorem for super-martingale, we know

[TABLE]

where $X=\#\{j:j\in{\cal B}^{c},W_{j}>0\}$. Since $X\sim\text{Binomial}(d_{0},1/2)$, we have

[TABLE]

Therefore $\text{E}[M(k_{T_{\alpha}})]\leq 1$. From (A.1), we have

[TABLE]

Since $\widehat{\cal A}(T_{\alpha})\subset\widehat{\cal A}_{1}$, we have $\#\{j:j\in{\cal A}^{c}\cap\widehat{\cal A}(T_{\alpha})\}=\#\{j:j\in{\cal B}^{c}\cap\widehat{\cal A}(T_{\alpha})\}$. Then, we conclude

[TABLE]

A.2 Proof of Theorem 5

Now we restrict ourselves to the subset $(\mathbf{X}^{(2)},\mathbf{Y}^{(2)})$. For $j\in\widehat{\cal A}_{1}$, define $\widehat{\omega}_{j}=\widehat{\mathrm{PC}}(\mathbf{X}^{(2)}_{j},\mathbf{Y}^{(2)})^{2}$, $\widetilde{\omega}_{j}=\widehat{\mathrm{PC}}(\widetilde{}\mathbf{X}^{(2)}_{j},\mathbf{Y}^{(2)})^{2}$, $\omega_{j}=\mathrm{PC}(X_{j},\mathbf{y})^{2}$, and $\omega_{j}^{\prime}=\mathrm{PC}(\widetilde{X}_{j},\mathbf{y})^{2}$. Theorem 1 implies that

[TABLE]

Thus we have

[TABLE]

On the other hand, part (i) in Lemma 1 implies that $\omega_{j}=\omega_{j}^{\prime}$ for all $j\in{\cal A}^{c}$. We have

[TABLE]

where the last equality is implied by Theorem 1 and the assumption $n_{2}>d$. Since $\min_{k\in{\cal B}}W_{k}\geq 4c_{3}n_{2}^{-\kappa}$ and $d<n_{2}$, (A.2) implies that $\min_{k\in{\cal B}}\widehat{W}_{k}\geq 2c_{4}n_{2}^{-\kappa}$ with probability at least $1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$. Together with (A.3), we know that $\min_{j\in{\cal B}}\widehat{W}_{j}>\max_{j\in{\cal B}^{c}}|\widehat{W}_{j}|$ with probability $1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$. In other words, with probability at least $1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$, the active features will be ranked ahead of the inactive features. Recall that in order to find the cutoff value $T_{\alpha}$, we start from the smallest value $t=|\widehat{W}_{d+1}|=0$, then move to the second smallest value $t=|\widehat{W}_{d}|$, then move to $t=|\widehat{W}_{d-1}|$, and so on. The procedure stops once it finds a value of $t$ satisfying (3.6). Restrict on the event ${\cal E}^{\prime}=\{\min_{j\in{\cal B}}\widehat{W}_{j}>\max_{j\in{\cal B}^{c}}|\widehat{W}_{j}|\}$, which holds with probability at least $1-O(n_{2}\exp\{-c_{4}n_{2}^{1-2\kappa}\})$ and let $t_{\min}=\min_{j\in{\cal B}}|\widehat{W}_{j}|=\min_{j\in{\cal B}}\widehat{W}_{j}$. If $1/s\leq\alpha$, then

[TABLE]

From the inequality above, we known this process must stop no later than $t$ reaches $t_{\min}$ and hence $T_{\alpha}\leq t_{\min}$. As a result,

[TABLE]

If $1/s>\alpha$, in order to satisfy (3.6), we must have $T_{\alpha}<t_{\min}$ or $T_{\alpha}>\max_{j}|\widehat{W}_{j}|$. $T_{\alpha}<t_{\min}$ indicates all active features are selected and $T_{\alpha}>\max_{j}|\widehat{W}_{j}|$ leads to an empty $\widehat{\cal A}(T_{\alpha})$, as stated in the second part. Notice that $|\widehat{W}_{s+k}|$ satisfies the rule (3.6) is equivalent to that at most $\lfloor(k-1)/(s+1)\rfloor$ of the signs of $\{\widehat{W}_{s+1},\dots,\widehat{W}_{s+k}\}$ are negative. Let $a_{k}$ be the probability that $|\widehat{W}_{s+k}|$ satisfies (3.6) and $b_{k}$ be the probability that the process stops at $|\widehat{W}_{s+k}|$ for $k=1,\dots,d-s$. Let $a_{0}=0$, we have

[TABLE]

It is easy to verify that for $s>2$,

[TABLE]

where $0<C(s)<1$ is some constant only depending on $s$. This completes the proof.

Acknowledgment: We would like to thank the AE and reviewers for their constructive comments, which lead to significant improvement of this work. Liu’s research was supported by National Natural Science Foundation of China (No. 11771361, 11671334, 11871409), Basic Scientific Center Project 71988101 of National Science Foundation of China and JAS14007. Li’s research was supported by NSF grants DMS 1820702 and 1953196. All authors equally contributed to this work, and the authors are listed in seniority. Both Yuan Ke and Jingyuan Liu are designated to be the corresponding authors of this paper.