More Powerful and General Selective Inference for Stepwise Feature Selection Using the Homotopy Method

Traditional statistical inference presumes that a statistical model and a statistical target for which we seek to conduct inference are determined before observing the dataset. Therefore, if we apply traditional methods to hypotheses selected after observing the dataset, the inferential results are no longer valid. This problem has been extensively discussed in the context of feature selection. In fact, even in commonly used feature selection methods such as SFS, correct assessment of the statistical reliability of selected features has long been difficult. Several approaches have been suggested in the literature toward addressing this problem (Benjamini and Yekutieli, 2005; Leeb and Pötscher, 2005; Leeb et al., 2006; Benjamini et al., 2009; Pötscher et al., 2010; Berk et al., 2013; Lockhart et al., 2014; Taylor et al., 2014). A particularly notable approach is conditional SI introduced in the seminal paper by Lee et al. (2016). In their work, the authors showed that, for a set of features selected by Lasso, the selection event can be characterized as a polytope by conditioning on the selected set of features as well as additional events on their signs. Furthermore, Tibshirani et al. (2016) showed that polytope-based SI is also applicable to SFS by additionally conditioning on the history and signs of sequential feature selection. Conditional SI has been actively studied in the past few years and extended to various directions (Fithian et al., 2015; Choi et al., 2017; Tian et al., 2018; Chen and Bien, 2019; Hyun et al., 2018; Loftus and Taylor, 2014, 2015; Panigrahi et al., 2016; Tibshirani et al., 2016; Yang et al., 2016; Suzumura et al., 2017; Tanizaki et al., 2020; Duy et al., 2020a).

In conditional SIs, it is typically preferred to condition on as little as possible so that the inference can be more powerful (Fithian et al., 2014). Namely, the main limitations of current polytope-based SI methods are that excessive conditioning is required to represent the selection event with a single polytope. In the seminal paper by Lee et al. (2016), the authors already discussed the problem of over-conditioning, and explained how extra conditioning on signs can be omitted by an exhaustive enumeration of all possible signs and by taking the union over the resulting polyhedra. However, such an exhaustive enumeration of exponentially increasing number of sign combinations is feasible only when the number of selected features is fairly small. Several other approaches were proposed to circumvent the drawbacks and restrictions of polytope-based SI. Loftus and Taylor (2015) extended polytope-based SI such that selection events characterized by quadratic inequalities can be handled, but this inherits the same over-conditioning problem. To improve the power, Tian et al. (2018) proposed an approach to randomize the algorithm in order to condition on less. Terada and Shimodaira (2019) proposed to use bootstrap re-sampling to characterize the selection event more generally. The disadvantage of these approaches is that additional randomness is introduced into the algorithm and/or the inference.

This study is motivated by Liu et al. (2018) and Duy and Takeuchi (2020). The former studied Lasso SI for full-model parameters, whereas the latter extended the basic idea of the former so that it can be also applied to Lasso SI in more general settings. These two studies go beyond polytope-based SI for more powerful Lasso SI without conditioning on signs. Because Lasso is formulated as a convex optimization problem, the selection event for conditional SI can be characterized based on the optimality conditions of the convex optimization problem. These two methods rely on the optimality conditions to characterize the minimum conditions for the Lasso feature selection event.

In contrast, the SFS algorithm cannot be formulated as an optimization problem. Therefore, existing conditional SI methods for SFS are based on conditioning on the history rather than the output of the SFS algorithm. Namely, existing SI methods actually consider conditions not only on which features are selected but also in what order they are selected. Moreover, when the SFS algorithm is extended to more complex procedures such as FB-SFS, the over-conditioning problem has more pronounced negative effects. Because features can be added and removed in various orders in FB-SFS, if we condition on the history rather than the output of the algorithm, we will suffer from extensive over-conditioning on the entire history of feature addition and removal.

Our contributions are as follows. First, we propose a novel conditional SI method for the SFS algorithm using the homotopy continuation approach and show that the proposed method overcomes the over-conditioning problem in existing methods and enables us to conduct minimally conditioned more powerful conditional SI on the selected features by the SFS algorithm. Second, we develop a conditional SI method for the FB-SFS algorithm with AIC-based stopping criteria by using the homotopy continuation approach. Although the conventional polytope-based SI method is more adversely affected by the over-conditioning problem when SFS is extended to FB-SFS, we show that our homotopy-based method still enables us to conduct minimally-conditioned SI and retains sufficiently high power. We demonstrate the validity and power of the proposed homotopy-based SI methods through numerical experiments on synthetic and real benchmark datasets. The codes will be released when the submission of this draft is accepted.

For a natural number $n$, we denote $[n]$ to be the set $\{1,\ldots,n\}$. The $n\times n$ identity matrix is denoted by $I_{n}$. For a matrix $X$ and a set of column indices ${\mathcal{M}}$, $X_{\mathcal{M}}$ indicates the matrix with the set of columns corresponding to ${\mathcal{M}}$.

We consider the standard forward SFS method as studied in Tibshirani et al. (2016) in §2 and §3. At each step of the SFS method, the feature that most improves the fit is newly added. When each feature has unit length, it is equivalent to the feature which is most correlated with the residual of the least-square regression model fitted with the previously selected features. For a response vector $\bm{y}\in\mathbb{R}^{n}$ and a set of features ${\mathcal{M}}\subseteq[p]$, let $\bm{r}(\bm{y},X_{\mathcal{M}})$ be the residual vector obtained by regressing $\bm{y}$ onto $X_{\mathcal{M}}$ for a set of features ${\mathcal{M}}$, i.e.,

where $P_{X_{\mathcal{M}}}=X_{\mathcal{M}}(X_{\mathcal{M}}^{\top}X_{\mathcal{M}})^{-1}X_{\mathcal{M}}^{\top}$. Let $K$ be the number of the selected features by the SFS method¹¹1 We discuss a situation where the number of features $K$ is selected by cross-validation in Appendix E. In other parts of the paper, we assume that $K$ is determined before looking at the data. . We denote the index of the feature selected at step $k$ as $j_{k}$ and the set of selected features up to step $k$ as ${\mathcal{M}}_{k}=\{j_{1},\ldots,j_{k}\}$ for $k\in[K]$ (Note, however, that, if there is no ambiguity, we simply denote the final set of the selected features as ${\mathcal{M}}~{}(={\mathcal{M}}_{K})$). Feature $j_{k}$ is selected as

where ${\mathcal{M}}_{0}(\bm{y}):=\emptyset$ and $\bm{r}(\bm{y},{\mathcal{M}}_{0}):=\bm{y}$.

In order to quantify the statistical significance of the relation between the selected features and the response, we consider a statistical test for each coefficient of the selected model parameters

Note that the $j^{\rm th}$ coefficient is written as $\hat{\beta}_{{\mathcal{M}},j}=\bm{\eta}^{\top}\bm{y}$ by defining

where $\bm{e}_{j}\in\mathbb{R}^{|{\mathcal{M}}|}$ is a unit vector, $j^{\rm th}$ element of which is 1 and 0 otherwise. Note that $\bm{\eta}$ depends on ${\mathcal{M}}$ and $j$, but we omit the dependence for notational simplicity. We consider the following statistical test for each coefficient $\beta_{{\mathcal{M}},j}=\bm{\eta}^{\top}\bm{\mu}$

where $\beta_{{\mathcal{M}},j}$ is the $j^{\rm th}$ element of the population least squares $\bm{\beta}_{{\mathcal{M}}}=P_{{\mathcal{M}}}\bm{\mu}$, i.e., the projection of $\bm{\mu}$ onto the column space of $X_{{\mathcal{M}}}$. Throughout the paper, we do not assume that the linear model is correctly specified, i.e., it is possible that $\bm{\mu}\neq X_{{\mathcal{M}}}\bm{\beta}$ for any $\bm{\beta}\in\mathbb{R}^{p}$. Even when the linear model is not correctly specified, $\bm{\beta}_{{\mathcal{M}}}$ is still a well-defined best linear approximation.

Because the target of the inference is selected by observing the data $\bm{y}$, if we naively apply a traditional statistical test to the problem in (4) as if the inference target is pre-determined, the result will be invalid (type-I error cannot be controlled at the desired significance level) owing to selection bias. To address the selection bias problem, we consider conditional SI introduced in Lee et al. (2016) and Tibshirani et al. (2016). Let us write the SFS algorithm as the following function:

which maps a response vector $\bm{y}\in\mathbb{R}^{n}$ to the set of the selected feature ${\mathcal{M}}$ by the SFS algorithm.

In conditional SI, the inference is conducted based on the following conditional sampling distribution of the test-statistic:

where

is the nuisance parameter, which is independent of the test statistic. The first condition ${\mathcal{A}}(\bm{Y})={\mathcal{A}}({\bm{y}})$ in (5) indicates the event that the set of the selected features by the $K$-step SFS method with a random response vector $\bm{Y}$ is ${\mathcal{M}}$, i.e., the same as those selected with the observed response vector $\bm{y}$. The second condition ${\bm{q}}({\bm{Y}})={\bm{q}}({\bm{y}})$ indicates that the nuisance parameter for a random response vector $\bm{Y}$ is the same as that for the observed vector $\bm{y}$²²2 The ${\bm{q}}(\bm{Y})$ corresponds to the component $\bm{z}$ in the seminal paper (see Lee et al. (2016), §5, Eq. 5.2, and Theorem 5.2). .

To conduct the conditional inference for (5), the main task is to identify the conditional data space

Once ${\mathcal{Y}}$ is identified, we can easily compute the pivotal quantity

where $F_{m,s^{2}}^{\mathcal{Z}}$ is the c.d.f. of the truncated Normal distribution with mean $m$, variance $s^{2}$, and truncation region ${\mathcal{Z}}$. Later, we will explain how ${\mathcal{Z}}$ is defined in (7). The pivotal quantity is crucial for calculating $p$-value or obtaining confidence interval. Based on the pivotal quantity, we can obtain selective $p$-value (Fithian et al., 2014) in the form of

where $\pi_{j}=1-F^{{\mathcal{Z}}}_{0,{\bm{\eta}}^{\top}\Sigma{\bm{\eta}}}(\bm{\eta}^{\top}{\bm{y}})$. Furthermore, to obtain $1-\alpha$ confidence interval for any $\alpha\in[0,1]$, by inverting the pivotal quantity in (7), we can find the smallest and largest values of $\bm{\eta}^{\top}\bm{\mu}$ such that the value of the pivotal quantity remains in the interval $\left[\frac{\alpha}{2},1-\frac{\alpha}{2}\right]$ (Lee et al., 2016).

Using the second condition in (6), the data in ${\mathcal{Y}}$ are restricted to a line (see Liu et al. (2018) and Fithian et al. (2014)). Therefore, the set ${\mathcal{Y}}$ can be re-written, using a scalar parameter $z\in\mathbb{R}$, as

where ${\bm{a}}={\bm{q}}(\bm{y})$, ${\bm{b}}=\Sigma{\bm{\eta}}({\bm{\eta}}^{\top}\Sigma{\bm{\eta}})^{-1}$, and

Here, $\bm{a}$, $\bm{b}$ and ${\mathcal{Z}}$ depend on ${\mathcal{M}}$ and $j$, but we omit the subscripts for notational simplicity. Now, let us consider a random variable $Z\in\mathbb{R}$ and its observation $z\in\mathbb{R}$ such that they respectively satisfy ${\bm{Y}}={\bm{a}}+{\bm{b}}Z$ and ${\bm{y}}={\bm{a}}+{\bm{b}}z$. The conditional inference in (5) is re-written as the problem of characterizing the sampling distribution of

Because $Z\sim\mathbb{N}(\bm{\eta}^{\top}\bm{\mu},{\bm{\eta}}^{\top}\Sigma{\bm{\eta}})$, $Z\mid Z\in{\mathcal{Z}}$ follows a truncated Normal distribution. Once the truncation region ${\mathcal{Z}}$ is identified, the pivotal quantity in (7) is obtained as $F^{{\mathcal{Z}}}_{\bm{\eta}^{\top}\bm{\mu},{\bm{\eta}}^{\top}\Sigma{\bm{\eta}}}(z)$. Thus, the remaining task is reduced to the characterization of ${\mathcal{Z}}$.

Unfortunately, full identification of the truncation region ${\mathcal{Z}}$ in conditional SI for the SFS algorithm is considered computationally infeasible. Therefore, in existing conditional SI studies such as Tibshirani et al. (2016), the authors circumvent the computational difficulty by over conditioning. Note that over-conditioning is not harmful for selective type-I error control, but it leads to the loss of power (Fithian et al., 2014). In fact, the decrease in the power due to over-conditioning is not unique problem for SFS in Tibshirani et al. (2016) but is a common major problem in many existing conditional SIs (Lee et al., 2016). In the next section, we propose a method to overcome this difficulty.