Model-free controlled variable selection via data splitting

Sufficient dimension reduction (SDR) is a powerful technique to extract relevant information from high-dimensional data (Li, 1991; Cook and Weisberg, 1991; Xia et al., 2002; Li and Wang, 2007). We use $Y$ with support $\Omega_{Y}$ to denote the univariate response, and let ${\bf X}=(X_{1},\ldots,X_{p})^{\top}\in\mathbb{R}^{p}$ be the $p$-dimensional vector of all covariates. The basic idea of SDR is to replace the predictor vector with its projection onto a subspace of the predictor space without loss of information on the conditional distribution of $Y$ given ${\bf X}$. In practice, a large number of features in high-dimensional data are typically collected, but only a small portion of them are truly associated with the response variable. However, while grasping important features or patterns in the data, the reduction subspace from SDR usually includes all original variables which makes it difficult to interpret. Therefore, in this paper, we aim at developing a model-free variable selection procedure to screen out truly non-contributing variables with certain error rate control, thus making the subsequent model building feasible or simplified and helping reduce the computational cost caused by high-dimensional data.

Let $F(Y\mid{\bf X})$ denote the conditional distribution function of $Y$ given ${\bf X}$. The index sets of the active and inactive variables are defined respectively as

Many prevalent variable selection procedures have been developed under the paradigm of linear models or generalized linear models, such as LASSO (Tibshirani, 1996), SCAD (Fan and Li, 2001), or adaptive LASSO (Zou, 2006). See the review of Fan and Lv (2010) and the book of Fan et al. (2020) for a fuller list of references. In contrast, model-free variable selection can be achieved by SDR since it does not require complete knowledge of the underlying model, thus researchers can avoid disposing of model misspecification.

SDR methods with variable selection aim to find the active set $\mathcal{A}$ such that

where “$\bot\!\!\!\bot$” stands for independence, ${\bf X}_{\mathcal{A}}=\left\{{\bf X}_{j}:j\in\mathcal{A}\right\}$ denotes the vector containing all active variables and ${\bf X}_{\mathcal{A}}^{c}$ is the complementary set of ${\bf X}_{\mathcal{A}}$. Condition (1) implies that ${\bf X}_{\mathcal{A}}$ contains all the relevant information in terms of predicting $Y$. Li et al. (2005) proposed to combine sufficient dimension reduction and variable selection. Chen et al. (2010) proposed a coordinate-independent sparse estimation that can simultaneously achieve sparse SDR and screen out irrelevant variables efficiently. Wu and Li (2011) focused on the model-free variable selection with a diverging number of predictors. A marginal coordinate hypothesis is proposed by Cook (2004) for model-free variable selection under low-dimensional settings, and then is promoted by Shao et al. (2007) and Yu and Dong (2016). Yu et al. (2016a) constructed marginal coordinate tests for sliced inverse regression (SIR) and Yu et al. (2016b) suggested a trace-pursuit-based utility for ultrahigh-dimensional feature selection. See Li et al. (2020) and Zhu (2020) for a comprehensive review.

However, those existing approaches do not account for uncertainty quantification of the variable selection, i.e., the global error rate control in the selected subset of important covariates in high-dimensional situations. In general high-dimensional nonlinear settings, Candes et al. (2018) developed a Model-X Knockoff framework for controlling false discovery rate (FDR, Benjamini and Hochberg, 1995), which was motivated by the pioneering Knockoff filter (Barber and Candès, 2015). Their statistics constructed via “Knockoff copies” would satisfy (or roughly) joint exchangeability and thus can yield finite-sample FDR control. However, the Model-X Knockoff requires knowing the joint distribution of the covariates, which is typically difficult in high-dimensional settings. Recently, Guo et al. (2024) improved the line of marginal tests (Cook, 2004; Yu and Dong, 2016) by using decorrelated score type statistics to make inferences for a specific predictor which is of interest in advance. They further leveraged the standard Benjamini and Hochberg (1995) on $p$-values to control FDR, but the intensive computation of the decorrelated process may limit its application to high-dimensional situations. In a different direction, Du et al. (2023) proposed a data splitting strategy, named symmetrized data aggregation (SDA), to construct a series of statistics with global symmetry property and then utilize the symmetry to derive a data-driven threshold for error rate control. Specifically, Du et al. (2023) aggregated the dependence structure into a linear model with a pseudo response and a fixed covariate, making the dependence structure become a blessing for power improvement. Similar to the Knockoff method, the SDA is also free of $p$-values and its construction does not rely on contingent assumptions, which motivates us to employ it in sufficient dimension reduction problems.

In this paper, we propose a model-free variable selection procedure that could achieve an effective FDR control. We first recast the problem of conducting variable selection in sufficient dimension reduction into making inferences on regression coefficients in a set of linear regressions with several response transformations. A variable selection procedure is subsequently developed via error rate control for low-dimensional and high-dimensional settings, respectively. Our main contributions include: (1) This novel data-driven selection procedure can control the FDR while being combined with different existing SDR methods for model-free variable selection by choosing different response transformation functions. (2) Our method does not need to estimate any nuisance parameters such as the structural dimension in SDR. (3) Notably, the proposed procedure is computationally efficient and easy to implement since it only involves a one-time split of the data and the calculation of the product of two dimension reduction matrices obtained from two splits. (4) Furthermore, this method can achieve finite-sample and asymptotic FDR control under some mild conditions. (5) Numerical experiments indicate that our procedure exhibits satisfactory FDR control and higher power compared with existing methods.

The rest of this paper is organized as follows. In section 2, we present the problem and model formulation. In section 3, we propose a low-dimensional variable selection procedure with error rate control and then discuss its extension in high-dimensional situations. The finite-sample and asymptotic theories for controlling the FDR are developed in Section 4. Simulation studies and a real-data investigation are conducted in Section 5 to demonstrate the superior performance of the proposed method. Section 6 concludes the paper with several further topics. The main theoretical proofs are given in Appendix. More detailed proofs and additional numerical results are delineated in the Supplementary Material.

Notations. Let $\lambda_{\min}({\bf B})$ and $\lambda_{\max}({\bf B})$ denote the smallest and largest eigenvalues of square matrix ${\bf B}=(b_{ij})$. Write $\|{\bf B}\|_{2}=(\sum\nolimits_{i}\sum\nolimits_{j}b_{ij}^{2})^{1/2}$ and $\|{\bf{B}}\|_{\infty}=\max_{i}\sum\nolimits_{j}|b_{ij}|$. Denote $\|\bm{\mu}\|_{1}=\sum\nolimits_{i}|\mu_{i}|$ and $\|\bm{\mu}\|_{2}=(\sum\nolimits_{i}\mu_{i}^{2})^{1/2}$ be the $L_{1}$ and $L_{2}$ norm of vector $\bm{\mu}$. Denote $\mathbb{E}({\bf X})$ and $\text{cov}({\bf X})$ be the expectation and covariance for random vector ${\bf X}$, respectively. Let $A_{n}\approx B_{n}$ denote that two quantities $A_{n}$ and $B_{n}$ are asymptotically equivalent, in the sense that there is a constant $C>1$ such that $B_{n}/C\leq A_{n}\leq B_{n}C$ with probability tending to 1. The “$\gtrsim$” and “$\lesssim$” are similarly defined.

The variable selection in (1) can be framed as a multiple testing problem

This is known as the marginal coordinate hypothesis described in Cook (2004) and Yu et al. (2016a). Some related works include Li et al. (2005), Shao et al. (2007), and Yu and Dong (2016). This type of selection procedure usually uses some nonnegative marginal utility statistics $W_{j}$’s to measure the importance of $X_{j}$’s to $Y$ in certain sense. However, the global error rate control within those methods is still challenging because the determination of selection thresholds generally involves the approximation to the distribution of $W_{j}$, and the accuracy of asymptotic distributions heavily affects the error rate control.

As a remedy, we consider a reformulation for (2). Let $\bm{\Sigma}=\mathrm{cov}({\bf X})>0$ and assume $\mathbb{E}({\bf X})=0$. Denote $\mathcal{C}$ is dimension reduction subspace. For any function $f(Y)$ satisfying $\mathbb{E}\left\{f(Y)\right\}=0$, it has been demonstrated by Yin and Cook (2002) and Wu and Li (2011) that

under the linearity condition Li (1991), which is usually satisfied when ${\bf X}$ is elliptical distribution. The transformation $f(\cdot)$ is used in a way different from its traditional role of being a mechanism for improving the goodness of model fitting. It serves as an intermediate tool for performing dimension reduction. Consequently, different transformation functions correspond to different SDR methods (Dong, 2021). One can choose a series of transformation functions, $f_{1}(Y),\ldots,f_{H}(Y)$, whose forms do not depend on data. The $H\left(>d\right)$, a pre-specified integer, is usually called a working dimension and $d$ is the true structural dimension of the subspace $\mathcal{C}$. Given the working dimension $H$, at the population level, define

Write $\textbf{B}_{0}=\left(\bm{\beta}_{1}^{0},\ldots,\bm{\beta}_{H}^{0}\right)$, then $\text{Span}(\textbf{B}_{0})\subseteq\mathcal{C}$, and $\text{Span}(\textbf{B}_{0})$ represents the subspace spanned by the column vector of $\textbf{B}_{0}$. By the following usual protocol in the literature of sufficient dimension reduction, we take one step further by assuming the coverage condition $\text{Span}(\textbf{B}_{0})=\mathcal{C}$ whenever $\text{Span}(\textbf{B}_{0})\subseteq\mathcal{C}$. This condition often holds in practice; see Cook and Ni (2006) for further discussion.

For $j=1,\ldots,p$, let $\beta_{hj}^{0}$ be the $j$th element of $\bm{\beta}_{h}^{0}\in\mathbb{R}^{p}$, $h=1,\ldots,H$. If the $j$th variable is unimportant, $\textbf{B}_{j}=\textbf{0}\in\mathbb{R}^{H}$, where $\textbf{B}_{j}=(\beta_{1j}^{0},\beta_{2j}^{0},\cdots,\beta_{Hj}^{0})^{\top}$ denotes the $j$th row of $\textbf{B}_{0}$. Further, $Y\bot\!\!\!\bot{\bf X}\mid{\bf X}_{\mathcal{A}}$ implies that $\sum\nolimits_{h=1}^{H}|\beta_{hj}^{0}|>0$ for $j\in\mathcal{A}$ and $\sum\nolimits_{h=1}^{H}|\beta_{hj}^{0}|=0$ for $j\in\mathcal{A}^{c}$. In other words, if $j$ belongs to the active set $\mathcal{A}$, response $Y$ must depend on $X_{j}$ through at least one of the $H$ linear combinations. If $j$ belongs to the inactive set $\mathcal{A}^{c}$, none of the $H$ linear combinations involve $X_{j}$ (Yu et al., 2016a). Accordingly, the testing problem (2) is equivalent to

Based on the above discussion, selecting active variables in a model-free framework is equivalent to selecting important variables in multiple response linear model. Assume that there are independent and identical distributed data $\mathcal{D}=\left\{{\bf X}_{i},Y_{i}\right\}_{i=1}^{2n}$. Denote $\widehat{\bm{\beta}}_{1},\ldots,\widehat{\bm{\beta}}_{H}$ are the estimators of $\bm{\beta}_{1}^{0},\ldots,\bm{\beta}_{H}^{0}$, and $W_{j}$’s as the marginal statistics based on the sample $\mathcal{D}$ associated with the variants of $\widehat{\bm{\beta}}_{1},\ldots,\widehat{\bm{\beta}}_{H}$. Its explicit form would be given in the next section. A selection procedure with a threshold $L$ is formed as

where $\widehat{\mathcal{A}}(L)$ is the estimate of $\mathcal{A}$ with threshold $L$. Obviously, $L$ plays an important role in variable selection to control the model complexity. We will construct an appropriate threshold by controlling the FDR to achieve model-free variable selection in SDR.

Denote $p_{0}=\left|\mathcal{A}^{c}\right|$, $p_{1}=\left|\mathcal{A}\right|$ and assume that $p_{1}$ is dominated by $p$, i.e., $p_{1}=o(p)$. The false discovery proportion (FDP) associated with the selection procedure (5) is

where $a\vee b=\max\left\{a,b\right\}$ and $\#\{\}$ stands for the cardinality of an event. The FDR is defined as the expectation of the FDP, i.e., $\text{FDR}(L)=\mathbb{E}\left(\text{FDP}(L)\right)$. Our main goal is to find a data-driven threshold $L$ that controls the asymptotic FDR at a target level $\alpha$,

In this section, we first provide a data-driven variable selection procedure in Subsection 3.1 to control the FDR via data splitting technique in a model-free context when $p<n$, and the high-dimensional version is deferred in Subsection 3.2.

In this section, we entirely focus on controlling FDR. We begin by imposing a mild restriction on the response transformation function $f(Y)$.

Assumption 1 distinguishes our approach from most model-based selection methods by transforming a general model into a multivariate response linear problem, thereby achieving model-free variable selection (Wu and Li, 2011). The transformed errors are not independent of the covariates and thus we need more effort for the theoretical analysis. Our first theorem is a finite sample theory for FDR control.

Theorem 4.1 holds regardless of the unknown relationship between variables ${\bf X}$ and response $Y$. This result can be established using the techniques developed in Barber et al. (2020). The quantity $\Delta_{j}$ is interpreted as a measure to investigate the effect of both the asymmetry of $W_{j}$ and the dependence between $W_{j}$ and $\textbf{W}_{-j}$ on FDR. In asymmetric cases, it is still expected that $\Delta_{j}$ will be small, given that both $\widehat{\beta}^{(1)}_{hj}$ and $\widehat{\beta}^{(2)}_{hj}$ converge to normal distributions if $n$ is not too small. Theorem 4.1 implies that tight control of $\Delta_{j}$’s under asymmetric cases also results in effective FDR control.

For asymptotic FDR control of the proposed procedure, we require the following technical assumptions, which are not the weakest one but facilitate the technical proofs in the Supplementary Material. Let $\|{\bf A}-\bm{\Sigma}_{\mathcal{S}}\|_{\infty}=O_{p}(a_{np})$ with $a_{np}\to 0$, where ${\bf A}=n^{-1}\sum\nolimits_{i=1}^{n}{\bf X}_{2\mathcal{S}i}{\bf X}_{2\mathcal{S}i}^{\top}$ and $\bm{\Sigma}_{\mathcal{S}}={\mathbb{E}}({\bf X}_{\mathcal{S}}{\bf X}_{\mathcal{S}}^{\top})$. Define $\upsilon_{n}=\max\left\{\|\bm{\Sigma}_{\mathcal{S}}\|_{\infty},\|\bm{\Sigma}^{-1}_{\mathcal{S}}\|_{\infty}\right\}$ and $\textbf{B}_{0\mathcal{S}}=\left\{\textbf{B}_{j}:j\in\mathcal{S}\right\}$. Denote $d_{n}=|\mathcal{A}|$, $q_{n}=|\mathcal{S}|$ and $q_{0n}=|\mathcal{S}\cap\mathcal{A}^{c}|$. Assume that $q_{n}$ is uniformly bounded above by some non-random sequence $\bar{q}_{n}$.

Theorem 4.2 implies that the variable selection procedure with the data-driven threshold $L$ can control the FDR at the target level asymptotically. Further investigations are needed to better understand the condition $c_{np}a_{np}\upsilon_{n}\sqrt{n\bar{q}_{n}}(\log{\bar{q}_{n}})^{3/2+\gamma}\to 0$. The conventional result of $\|{\bf A}-\bm{\Sigma}_{\mathcal{S}}\|_{\infty}$ indicates that $a_{np}=O_{p}(\upsilon_{n}\sqrt{\log\bar{q}_{n}/n})$. With $c_{np}=d_{n}\sqrt{\log p/n}$ of LASSO selector, the condition degenerates to $d_{n}\upsilon_{n}^{2}\sqrt{\bar{q}_{n}/n}\to 0$ if $p$ is of a polynomial rate of $n$. The above condition basically imposes restrictions on the rate of $d_{n}$, $\upsilon_{n}$, and $\bar{q}_{n}$. Accordingly, the screening stage on split $\mathcal{D}_{1}$ must satisfy $\bar{q}_{n}=o(n)$ if we assume that $d_{n}$ and $\upsilon_{n}$ are bounded. Alternatively, if we only assume that $\upsilon_{n}$ is bounded, then a sufficient requirement for the condition in Theorem 4.2 is $\bar{q}_{n}=o(n^{1/2})$ since $d_{n}\leq\bar{q}_{n}$. This is a reasonable rate in the problem with a diverging number of parameters, such as Fan and Peng (2004) and Wu and Li (2011).

We evaluate the performance of our proposed procedure on several simulated datasets and a real-data example under low-dimensional and high-dimensional settings.

Identifying the truly contributory variables is a critical task in statistical inference. In this paper, we proposed a novel model-free variable selection procedure with a data-driven threshold in sufficient dimension reduction framework for a family of sliced methods via data splitting. The proposed MFSDS is computationally efficient in high-dimensional settings. Theoretical and numerical results show that the proposed MFSDS can asymptotically control the FDR at the target level.

Our work raises several intriguing open questions that deserve further investigation. First, it is interesting to adopt multiple data splitting in MFSDS to improve the stability and robustness but it may require intensive computations. Equal size data splitting is used in our simulation and theory for simplicity. Second, we have checked that a larger sample size for $\mathcal{D}_{1}$ will improve TPR to some extent. In addition, we sacrifice a little detection power to obtain symmetry property via data splitting, while how to construct symmetry property without data splitting or deriving complex asymptotic distribution may be another promising direction. Third, benefiting from the use of the response transformation, the proposed method can work with many classical SDR methods (such as sliced inverse regression (Li, 1991), covariance inverse regression estimator (Cook and Ni, 2006)) and many types of response (such as continuous, discrete, or categorical data). Finally, our method may also be able to dispose of the order determination problem in many situations which merits further investigation, such as Luo and Li (2016) if the importance of $H$ slices has been determined in advance with a diverging number of $H$. It is also worth further exploring the optimal selection of working dimension $H$ in practice.