A Transparent and Nonlinear Method for Variable Selection ¹¹footnotemark: 1

Let $\mathbf{Y}=(y_{1},\cdots,y_{n})^{\prime}$ be the response, and $\mathbf{X}=(\mathbf{X}_{1},\cdots,\mathbf{X}_{p})$ be the $p$-dimensional input predictors, where $\mathbf{X}_{j}=(x_{1j},\cdots,x_{nj})^{\prime}$, $j=1,\cdots,p$, and $n$ is the sample size. Let $\mathcal{X}=\{1,\ 2,\cdots,p\}$ be the index set of the input predictors. For any subset of indices $\mathcal{S}\subseteq\mathcal{X}$, let $\mathbf{X}_{\mathcal{S}}$ be the data matrix composed of all the $\mathbf{X}_{j}$ that satisfies $j\in\mathcal{S}$. For an $\mathcal{S}$, if there is a function $f(\cdot)$ that makes $\mathbf{Y}=f(\mathbf{X}_{\mathcal{S}})+\mathbf{\varepsilon}$, where $\mathbf{\varepsilon}$ is a stochastic error, then $\mathcal{S}$ stands for the index set of the true model, which is named the relevant subset in this paper. In high-dimensional data, $p$ is close to or even larger than $n$. In such cases, sparsity assumption usually holds that the dimension of $\mathbf{X}_{\mathcal{S}}$ is much less than $p$.

In real problems, $\mathbf{Y}$ and $\mathbf{X}_{\mathcal{S}}$ can be linearly correlated, or they can also be associated in the form of a variety of complex nonlinear functions. To date, handling complex nonlinear correlations remains a great challenge (Fan et al., 2020). The primary goal of the proposed variable selection method is to obtain the correlated subset $\mathcal{S}$ from the input set $\mathcal{X}$. For a given group of predictors $\mathbf{X}_{\mathcal{G}}$, $\mathcal{G}\subseteq\mathcal{X}$, and an unknown predictor $\mathbf{X}_{j}$, $j\in\mathcal{X}\ \backslash\ \mathcal{G}$, the Conditional Dependence Coefficient (CODEC) (Azadkia & Chatterjee, 2021) can measure the nonlinear conditional correlation between $\mathbf{Y}$ and $\mathbf{X}_{j}$ given $\mathbf{X}_{\mathcal{G}}$, as well as the interaction between $\mathbf{X}_{j}$ and $\mathbf{X}_{\mathcal{G}}$ in explaining $\mathbf{Y}$. One of the most important features of the CODEC is that it converges to a limit in $[0,~{}1]$. Given $\mathbf{X}_{\mathcal{G}}$, the limit is $0$ if and only if $\mathbf{Y}$ and $\mathbf{X}_{j}$ are conditionally independent, and is $1$ if and only if $\mathbf{Y}$ is almost surely equal to a measurable function of $\mathbf{X}_{j}$. If $\mathcal{G}\neq\emptyset$, let $\mathcal{G}=\{g_{1},\cdots,g_{q}\}$, where $q\geq 1$. Let $\mathbf{x}_{i}^{\mathcal{G}}=(x_{ig_{1}},\cdots,x_{ig_{q}})$ be the $i$-th observation of $\mathbf{X}_{\mathcal{G}}$, where $i=1,\cdots,n$. CODEC is calculated as

where $R_{h}$ denotes the rank of observation $y_{h}$, i.e., the number of $i$ such that $y_{i}\leq y_{h}$. $M(h)$ denotes the index $i$ of the nearest neighbor $(\mathbf{x}_{i}^{\mathcal{G}},x_{ij})$ of $(\mathbf{x}_{h}^{\mathcal{G}},x_{hj})$ with respect to the Euclidean metric on $\mathbb{R}^{q+1}$, $N(h)$ denotes index $i$ such that $\mathbf{x}_{i}^{\mathcal{G}}$ is the nearest neighbor of $\mathbf{x}_{h}^{\mathcal{G}}$ in $\mathbb{R}^{q}$, and the ties are broken uniformly at random for both $M(h)$ and $N(h)$. $R_{M(h)}$ denotes the rank of $y_{M(h)}$, i.e., the number of $i$ such that $y_{M(i)}\leq y_{M(h)}$, and $R_{N(h)}$ denotes the rank of $y_{N(h)}$.

The CODEC can also measure the unconditional correlation between $\mathbf{Y}$ and a predictor $\mathbf{X}_{j}$, where $j\subseteq\mathcal{X}$, in the absence of any given predictors, i.e., when $\mathcal{G}=\emptyset$. In this case, the CODEC is interpreted as an unconditional dependence coefficient, and is calculated as

where $R_{h}$ denotes the rank of $y_{h}$. $M(h)$ denotes the index $i$ of the nearest neighbor $x_{ij}$ of $x_{hj}$, and the ties are broken uniformly at random for $M(h)$. $R_{M(h)}$ denotes the rank of $y_{M(h)}$. $L_{h}$ denotes the number of $i$ such that $y_{i}\geq y_{h}$.

The calculations above require continuous $\mathbf{Y}$, $\mathbf{X}_{\mathcal{G}}$ and $\mathbf{X}_{j}$, but the CODEC can also be applied to measure the correlations between discrete predictors if the ties are unbonded randomly. If the denominator of $T_{n}(\mathbf{Y},\mathbf{X}_{j}\mid\mathbf{X}_{\mathcal{G}})$ is $0$, the CODEC is undefined. At this point, $Y$ is almost surely equal to a measurable function of $\boldsymbol{X}_{\mathcal{G}}$, and $\boldsymbol{X}_{\mathcal{G}}$ can be regarded as sufficient for predicting $\mathbf{Y}$.

High-dimensional data usually contains uninformative predictors because of some restrictions in data collection. Removing these predictors has limited influence on explaining the response (Li et al., 2017a). For example, in the study of handwriting digits, pixels in the marginal areas may have little explanatory power and can be ignored in the construction of deep neural networks (Chen et al., 2021). Real datasets do not always contain uninformative predictors, since the predictors are usually carefully chosen based on expert experience before being collected to save storage space. However, if we have little prior knowledge, and obtain a dataset with a large number of uninformative predictors, removing them beforehand can greatly improve the efficiency of variable selection.

In this paper, Shannon entropy is used as the measure to distinguish the uninformative predictors from the remaining predictors. Shannon entropy is the probability of all possible values of a predictor, which represents the expectation of the amount of information contained in the predictor (Gray, 2011). Consider a discrete predictor $\mathbf{X}_{j}$ that takes a finite number of $c$ possible values $x_{j,k}\in\{x_{j,1},\cdots,x_{j,c}\}$ with corresponding probabilities $p_{j,k}\in\{p_{j,1},\cdots,p_{j,c}\}$. Its entropy $H(\mathbf{X}_{j})$ is defined as

In general, there is $0\leq H(\mathbf{X}_{j})\leq\ln{c}$. If the distribution of $\mathbf{X}_{j}$ is highly biased toward one of the possible value $x_{i,k}$, $H(\mathbf{X}_{j})$ is the lowest. At this point, if $H(\mathbf{X}_{j})=0$, $\mathbf{X}_{j}$ is defined as an uninformative predictor, i.e., the quantity of information contained in $\mathbf{X}_{j}$ is $0$. Shannon entropy can be applied only to discrete predictors, and data discretization is required beforehand for continuous predictors (Brown et al., 2012).

If a candidate predictor and a subset of relevant predictors are collinear, keeping both of them may affect the robustness of model estimation (Yu & Liu, 2004). Although there is a consensus on the adverse effects of multicollinearity (Fan & Lv, 2010), most variable selection methods manage to avoid discussing the issue. In this paper, such predictors are named redundant predictors, and are measured separately from other predictors to be deleted.

Gram‒Schmidt Orthogonalization (GSO) is adopted in this paper to decompose and identify the information contained in the predictor set, and further to measure multicollinearity among predictors (Wang et al., 2020; Lyu et al., 2017). The Gram‒Schmidt theorem in Euclidean space indicates that for an index subset of predictors $\mathcal{G}=\{g_{1},\cdots,g_{q}\}\subseteq\mathcal{X}$ and the corresponding matrix $\mathbf{X}_{\mathcal{G}}=(\mathbf{X}_{g_{1}},\cdots,\mathbf{X}_{g_{q}})$, if $\mathbf{X}_{\mathcal{G}}$ is linearly independent, one can always construct an orthogonal basis $\mathbf{Z}_{\mathcal{G}}=(\mathbf{Z}_{g_{1}},\cdots,\mathbf{Z}_{g_{q}})$ via GSO, where $\mathbf{Z}_{\mathcal{G}}$ is the linear combination of $\mathbf{X}_{\mathcal{G}}$, and it spans the same space as $\mathbf{X}_{\mathcal{G}}$.

For each $\mathbf{X}_{j}$, $j\in\mathcal{G}$, the orthogonalized variable $\mathbf{Z}_{j}$ is

where $\langle\mathbf{X}_{j},\mathbf{Z}_{k}\rangle$ is the inner product of $\mathbf{X}_{j}$ and $\mathbf{Z}_{k}$, and $\|\mathbf{Z}_{k}\|_{2}$ is the $L^{2}$ norm of $\mathbf{Z}_{k}$.

Let $q$ be the rank of $\mathbf{X}=(\mathbf{X}_{1},\cdots,\mathbf{X}_{p})$. Suppose $q<p$, and $\mathbf{X}_{\mathcal{G}}$ is linearly independent, where the index subset $\mathcal{G}\subseteq\mathcal{X}$; then, an orthogonal basis $\mathbf{Z}=\mathbf{Z}_{\mathcal{G}}$ of $\mathbf{X}$ can be obtained using GSO. For any $j\in\mathcal{X}\ \backslash\ \mathcal{G}$, one can obtain the orthogonal variable $\mathbf{Z}_{j}$ of $\mathbf{X}_{j}$ with GSO based on $\mathbf{Z}$. If the variance of $\mathbf{Z}_{j}$ satisfies $Var(\mathbf{Z}_{j})=0$, all the information of $\mathbf{X}_{j}$ is already contained in $\mathbf{X}_{\mathcal{G}}$, and accordingly, $\mathbf{X}_{j}$ is a redundant predictor, which is collinear with $\mathbf{X}_{\mathcal{G}}$.

In addition to uninformative and redundant predictors, there is another group of predictors that needs to be deleted. One of the most obvious features is that $\mathbf{Y}$ is conditionally independent of these predictors given the relevant subset $\mathcal{S}$. In this paper, these predictors are called conditionally independent predictors. Modeling with the conditionally independent predictors does more harm than good. Not only will it damage interpretability, but it may also reduce prediction accuracy. The CODEC can be used to measure conditionally independent predictors. In particular, given $\mathcal{S}$, if there is $T_{n}(\mathbf{Y},\mathbf{X}_{j}\mid\mathbf{X}_{\mathcal{S}})$ close to $0$ for a predictor $\mathbf{X}_{j}$, where $j\in\mathcal{X}\ \backslash\ \mathcal{S}$, $\mathbf{X}_{j}$ can be considered a conditionally independent predictor.

We generate a regression simulation problem with both monotonic and nonmonotonic functional dependence, interactions, uninformative predictors, multicollinearity, and conditionally independent predictors. In the problem, $90\%$ of the $p$-dimensional input predictors have information, and the remaining $10\%$ are uninformative. The first $90\%$ predictors are consisted of 9 signals $\mathbf{X}_{t_{1}},\cdots,\mathbf{X}_{t_{9}}$, where $t_{g}=(g-1)\cdot p/10+1,\ g=1,\cdots,9$, and the rest are redundant predictors that are collinear with these $9$ signals. The signals are independently sampled from the standard normal distribution $\mathcal{N}(0,~{}1)$. The redundant predictors are highly correlated with one of the signals, and there is $\mathbf{X}_{t_{g}+j}=\mathbf{X}_{t_{g}}+\lambda\mathbf{\epsilon}_{t_{g}+j},\ g=1,\cdots,9,\ j=1,\cdots,(p/10)$, where $\mathbf{\epsilon}_{t_{g}+j}$ is the stochastic error that follows $\mathcal{N}(0,~{}1)$, and $\lambda$ is set as $0.01$. The last $10\%$ predictors are uninformative with only a few nonzero observations. The proportion of nonzero observations is set as $0.001$, and nonzero observations are randomly generated from $\mathcal{N}(0,~{}0.1^{2})$. The response $\mathbf{Y}$ is a nonlinear measurable function of some predictors

where $\mathbf{\varepsilon}$ is the stochastic error that follows $\mathcal{N}(0,~{}0.1^{2})$. $\mathbf{Y}$ is nonlinearly correlated with these four predictors or any of their collinear predictors. $\mathbf{Y}$ is oscillatory functionally dependent on $\mathbf{X}_{t_{3}}$ and $\mathbf{X}_{t_{4}}$, and there are interactions between $\mathbf{X}_{t_{1}}$ and $\mathbf{X}_{t_{2}}$, and $\mathbf{X}_{t_{3}}$ and $\mathbf{X}_{t_{4}}$. The problem has the following three settings by specifying $n$ and $p$. For each setting, $10$ datasets are generated.

• Setting 1: $n=2\ 000,~{}p=1\ 000$;

• Setting 2: $n=2\ 000,~{}p=2\ 000$;

• Setting 3: $n=2\ 000,~{}p=5\ 000$.

In our experiments, TNVS is evaluated against six representative baselines, including two up-to-date model-free and nonlinear methods, namely FOCI (Azadkia & Chatterjee, 2021) and WLS (Zhong et al., 2021), two classic model-free and nonlinear feature screening methods, namely DC-SIS (Li et al., 2012) and SIRS (Zhu et al., 2011), and two stepwise methods based on linear partial correlation, namely PC-simple (Buhlmann et al., 2010) and TPC (Li et al., 2017b).

The parameter settings of the variable selection methods are listed as follows. For TNVS, the uninformative threshold $\alpha_{1}$, relevant threshold $\alpha_{2}$, and redundant threshold $\alpha_{3}$ are empirically set as $0.01,\ -0.01,\ 0.01$, as the default setting in line with previous research (Liu et al., 2018). $d_{\max}$ is set as $\lceil n/\log{n}\rceil$, which represents the smallest positive integer no less than $n/\log{n}$. The termination criterion of FOCI is that the CODECs of all candidates are no larger than 0, or the number of selected predictors reaches $\lceil n/\log{n}\rceil$. The significance levels of PC-simple and TPC are set as $0.05$ (Buhlmann et al., 2010; Li et al., 2017b). The number of selected predictors of WLS, DC-SIS and SIRS are empirically set as $\lceil n/\log{n}\rceil$ in line with the literature. Considering the computing capability, we set the upper bound of the CPU time as $3\ 600$ seconds for all stepwise methods, including TNVS, FOCI, PC-Simple and TPC. In addition, predictors with zero variance are removed beforehand for FOCI. If a selected subset is empty, i.e., no predictor is considered to be correlated with the response, the mean of the response is used as the predictive value for regression problems, and a random category of the response is chosen as the predictive value for classification problems.

10-fold cross-validation is performed on each dataset, where the $n$ samples of the dataset are randomly divided into $10$ equal folds. Each unique fold of samples is used for testing and the remaining $9n/10$ samples are used for training. The performance in each setting is the average result of all folds on all the $10$ corresponding datasets, i.e., each variable selection method is tested for $100$ times in each setting to evaluate their variable selection capability.

All the variable selection methods are implemented in R 4.0.3 for a fair comparison²²2Our code is available at https://github.com/kywang95/TNVS.. The experiments are performed using an Intel Core i7 3.4 GHz processor with 16 GB of RAM.

The following four criteria are adopted to evaluate the effectiveness and efficiency of TNVS and baselines on high-dimensional simulations with complex nonlinear functional dependencies. The results of each simulation in $100$ repetitions are summarized in Table 2.

1. 1.

   $P_{a}$: probability of the discovery of all truly important predictors over all repetitions. The closer $P_{a}$ is to $1$, the more robust the variable selection method is to retain all true predictors. For simulations, true predictors are any of the relevant predictors and the redundant predictors that are collinear with the relvant predictors. For example, for Setting $1$, the indices of true predictors are any combinations of four elements from the subsets $\{1,~{}2,\cdots,100\}$, $\{101,~{}102,\cdots,200\}$, $\{201,~{}202,\cdots,300\}$, and $\{301,~{}302,\cdots,400\}$, and $P_{a}$ is the frequency that contains at least one element from each of the four sets.

2. 2.

   $\mathcal{M}$: the minimum model size required to include all true predictors in all repetitions in which none of the true predictors are omitted. The closer $\mathcal{M}$ is to the true model size, the more concise the selected subset is. The average and the standard deviation of $\mathcal{M}$ are obtained.

3. 3.

   Coverage: The number of true predictors covered by the selected subset. Coverage is less than or equal to the number of true predictors. The average and standard deviation of the coverage in all repetitions are obtained. In the simulated datasets, the closer the coverage is to $4$, the fewer true predictors are omitted.

4. 4.

   Time: the average and standard deviation of running time (in seconds) of variable selection in $100$ repetitions.

Table 2 shows that our TNVS, and the recently proposed FOCI and WLS, outperform the other traditional methods, and TNVS is the most effective and efficient one among these three methods. First, the-close to-1 $P_{a}$ in all settings indicates that these methods can reserve all true predictors in most cases. Among these three methods, TNVS has the smallest minimum model size, which means that the subset selected by TNVS is the most concise. In addition, computational time of TNVS is the least among these three methods, which shows that it is the most efficient. The framework of TNVS is most similar to that of the FOCI. Comparing these two methods, TNVS has a smaller minimum model size and a larger coverage, indicating that TNVS is more accurate and reliable than FOCI.

Table 2 also shows that these nonlinear simulation problems are extremely difficult to solve with the other four prevailing methods, i.e., SIRS, DC-SIS, PC-simple and TPC. Their $P_{a}$ are close to 0, indicating that these methods omit some true predictors and are invalid. PC-simple and TPC are designed for linear regression models, and close-to-0 coverages indicate their weak abilities to identify nonlinearly relevant predictors. SIRS and DC-SIS can only identify the monotonic nonlinear relevance, and having difficulty in removing the multicollinearity, and close-to-0 coverages indicate they cannot retain all important predictors.

Experiments on simulations demonstrate that every module of the proposed TNVS is effective and indispensable. First, compared with PC-simple and TPC, these two stepwise methods with linear partial correlation, TNVS adopts the nonlinear CODEC. This leads to the unsatisfying $P_{a}$ of these two baselines and the close-to-1 $P_{a}$ of TNVS, which indicates that introducing CODEC ensures that TNVS can effectively identify the complex nonlinear associations, such as interacted predictors, and oscillatory functional dependence with interactions. Second, FOCI is a stepwise method with CODEC, but ignoring to sorted disposal the uninformative and redundant predictors. Compared with FOCI, TNVS avoids retaining the large percentages of these unnecessary predictors, and achieves larger $P_{a}$, smaller model size, larger coverage, and less computational time. This implies that prefiltering and batch deletion inherently improve effectiveness and efficiency of TNVS.

We further demonstrate the model interpretability of the proposed variable selection method. Since in simulations, the truly relevant, uninformative, redundant, and conditionally independent predictors are known in advance, we can obtain the proportions of four types of predictors in the selected subset for all methods. Taking a closer look at the three competitive methods mentioned in Table 2, Table 3 shows the proportions of the ground truth relevant ($\textrm{Rel}_{GT}$), uninformative ($\textrm{Uin}_{GT}$), redundant ($\textrm{Red}_{GT}$), and conditionally independent ($\textrm{Cind}_{GT}$) predictors in the selected subset of TNVS, FOCI, and WLS. The average of each quantity in all $100$ repetitions are summarized. The truly relevant predictors in the selected subset are true positive, whereas the rest three types are false positive. Thus, $\textrm{Rel}_{GT}$ is actually the precision of each methods. The larger $\textrm{Rel}_{GT}$ is, and the smaller the proportions of $\textrm{Rel}_{GT}$, $\textrm{Uin}_{GT}$, and $\textrm{Red}_{GT}$ are, the more concise the selected subset is.

Table 3 shows that TNVS has the largest proportion of $\textrm{Red}_{GT}$ in its selected subsets, or equivalently, the largest precision. Nearly all the predictors selected by TNVS are relevant predictors. Uninformative, redundant, and conditionally independent predictors are rarely included in the selected subset of TNVS. The subset selected by FOCI contains a large percentage of uninformative predictors and some redundant predictors. Although WLS can select all relevant predictors, they only take a small part of its selected subset, and many unimportant predictors are falsely picked before the truly important predictors. These findings demonstrate that TNVS enhances the capability of FOCI to handle the lack of information and multicollinearity, and it also achieves a better balance between accuracy and conciseness than WLS.

The most significant difference of TNVS against other methods is that TNVS can divide the predictors into four disjoint subsets, instead of only the selected and the deleted ones. Table 4 shows the proportions of the ground truth types of predictors in the four subsets obtained by TNVS (the relevant subset $\textrm{Rel}_{pred}$, the uninformative subset $\textrm{Uin}_{pred}$, the redundant subset $\textrm{Red}_{pred}$, and the conditionally independent subset $\textrm{Cind}_{pred}$). The average of each quantity in all $100$ repetitions are summarized. The truly relevant predictors grouped into the selected relevant subset are true positive, whereas those grouped into the other three types are false negative. Thus, elements in the diagonal are actually the recall of TNVS. The larger elements in the diagonal are, and the smaller the other elements are, the more accurate the variable selection is.

Table 4 demonstrates the unique strength of TNVS. Besides the selected relevant subset, it further divides the predictors to be removed into the uninformative, redundant, and conditionally independent subsets. The diagonal elements show that the recall of TNVS is high on all four subsets. All uninformative predictors are correctly deleted. A very few relevant predictors are falsely identified as conditionally independent in rare cases, and the misclassification of relevant predictors hinders TNVS from distinguishing the redundancy and conditionally independence. Sufficient evidence shows that when tackling simulations with complex nonlinear functional dependencies, uninformative predictors and collinearity, TNVS provides credible explanations on why certain predictors are removed. To sum up, TNVS outperforms the baselines on model interpretability.

Four prevailing real datasets from various domains are chosen to demonstrate the performance of the proposed TNVS (Li et al., 2017a). The descriptive statistics of the datasets are listed in Table 5. In this paper, the original probes, features, or pixels are directly used as input predictors without any preprocessing such as feature extraction, since the aim of these experiments is to demonstrate the predictive capability and interpretability of the proposed variable selection method in high-dimensional data, not to design a distinct method for specific domain, such as face images recognition.

The reasons for choosing these datasets are as follows. First, these biological or image datasets are high-dimensional data with a large number of predictors. The response and predictors are often nonlinearly correlated in such datasets, and the inputs usually contain a large number of uninformative, redundant, or conditionally independent predictors. Moreover, the face images are considered as the preference since we can easily visualize the variable selection results on the original images, and literature on facial landmark detection (Köstinger et al., 2011) has provided some key reference points. By checking the relevance between our selected predictors and these key reference points, we can easily observe whether the selection results have evident post hoc interpretability. Two ORL datasets with different resolutions are selected to examine whether there are more redundant predictors in images with higher resolution, i.e., whether the images with lower resolution is distinct enough for classification. The datasets are described in detail.

1. 1.

   arcene is a binary classification dataset to discriminate between patients with and without cancer. The data are collected from mass spectrometry analysis of blood serum, and the predictors are continuous. The sample contains $88$ observations with ovarian or prostate cancer and $112$ observations of healthy patients. Each observation has $7\ 000$ real features and $3\ 000$ random probes. These random probes are permuted least-informative feats, and some of them can probably be identified as uninformative predictors.

2. 2.

   isolet is an audio dataset for predicting what letter name is spoken. The dataset contains $1\ 560$ observations, which are audio recordings of $30$ speakers. Each subject speaks names of $26$ letters of the alphabet twice. Each observation contains $617$ features. In this experiment, we generate a high-dimensional dataset with $300$ observations from the raw isolet: the first $6$ categories are considered, and $50$ observations are randomly sampled from each category.

3. 3.

   warpAR10P is a facial recognition dataset with 10 distinct subjects. For each subject, $13$ images are captured under various lighting conditions, with different facial expressions and occlusions such as scarfs and sunglasses. The size of each image is $60\times 40$ pixels, with $256$ gray levels per pixel.

4. 4.

   ORL is another facial recognition dataset containing $400$ face images collected from $40$ individuals, and each image contains $64\times 64$ pixels with $256$ gray levels per pixel. These images exhibit multiple changes, such as various picture-taking times, illumination changes, different expressions, and facial detail changes.

Note that Z-score normalization is adopted for input predictors in all datasets, i.e., $\tilde{\mathbf{X}}_{j}=(\mathbf{X}_{j}–\mu_{j})/\sigma_{j}$, where $\mu_{j}$ and $\sigma_{j}$ are the mean and standard deviation of predictor $\mathbf{X}_{j}$, respectively.

The proposed TNVS is inherently model-free and the selection process relies on no assumptions of the learning model. The effectiveness and generality of the selected subsets are consequently evaluated in terms of their performance on multiple learning models. Specifically, five prevailing predictive models commonly adopted in literature are considered (Salesi et al., 2021; Wan et al., 2022), including Support Vector Regression (SVR) / Support Vector Machine (SVM), Random Forest (RF), Decision Tree (DT), Light Gradient Boosting Machine (lightGBM) and Multi-layer Perceptron (MLP). To examine the predictive capabilities of these variable selection methods, each learning model is trained with the selected predictors, and the results are then compared. These experiments aim to demonstrate the superior generality of variable selection methods rather than the enhanced prediction performance of well-tuned models. Thus, the configurations of the learning models are set to be identical for all variable selection methods to avoid favoring any individual method by combining it with more competitive models. Most parameters are kept as the default settings in scikit-learn (Buitinck et al., 2013). The only adjustment is that the number of iterations is set as $5\ 000$ for MLP. All the learning models are implemented in Python 3.9 using the application programming interface (API) of scikit-learn.

For real datasets, repeated $10$-fold cross-validation is adopted to evaluate each combination of variable selection methods and learning models. The $10$-fold cross-validation is repeated $10$ times on each dataset. All variable selection methods and learning models share the same cross-validation splits so that unbiased estimates are obtained. For all splits, the ratio of classes in the training set is kept the same as that in the full sample. The predictive results are evaluated with out-of-sample accuracy (ACC), recall, $F_{1}$, and the Cohen Kappa Score (Kappa). The ranges of ACC, recall and $F_{1}$ are within $[0,~{}1]$; the larger these indicators are, the better the prediction is. The range of Kappa is within $[-1,~{}1]$; the closer the score is to 1, the better the prediction is.

Since the true models of real datasets are unknown, the predictive results of selected subsets are used to indirectly demonstrate the effectiveness of variable selection. The predictive results obtained using multiple combinations of variable selections and learning models are demonstrated as follows. First, Fig. 3 shows the boxplot of the four predictive indicators obtained with lightGBM on warpAR10P, which compares the predictive capability of the proposed TNVS with the baselines in this case. The predictive results of TNVS are better than those of the baselines, especially PC-Simple and TPC. Other statistical results of the average ACC, recall, $F_{1}$, and Kappa on the real datasets are listed in Table LABEL:tab:appendix_A1 of Appendix B. The results of other datasets and other learning models substantially agree with those in Fig. 3, which demonstrates that TNVS is robustly superior to the baselines.

A statistical test is performed to demonstrate the predictive capability of the proposed TNVS. For every dataset, the Friedman test is performed on statistical results of the four indicators obtained by all learning models to determine the overall performance of a certain variable selection method. The Friedman test provides a mean rank of these methods, and indicates whether there are significant differences among these methods. A higher mean rank means that the method has a better-ranked performance over the competitors. Table 6 shows the result of Friedman test, including the mean rank of these methods, the order of the mean rank, and the p-value of the Friedman test. As seen, the predictors selected by TNVS attain the top two among the seven methods on real datasets under the condition of retaining less predictors than feature screening methods, and there are significant differences among these methods, which further confirms the strength and robustness of TNVS in prediction enhancement.

TNVS transparently groups the predictors into four subsets, i.e., relevant, uninformative, redundant, and conditionally independent, which improves the model interpretability of TNVS. The average cardinalities of the four subsets are counted to check whether these types of predictors exist in the real datasets. Besides, the average compression rate is calculated, which is defined as the proportion of relevant predictors in the inputs. The results of $100$ repetitions on each dataset are shown in Table 7.

As shown in Table 7, the percentage of relevant predictors is lower than $6\%$ in all datasets, and is even lower than $1\%$ in arcene, warpAR10P, and ORL, which indicates that TNVS compresses the input set of predictors to a large degree. Redundant predictors are identified in all datasets, which shows the usefulness of GSO. TNVS identifies the uninformative predictors in arcene as expected.

In addition to model interpretability, we also visualize variable selection results of the ORL dataset as examples to demonstrate post hoc interpretability of TNVS. Under the premise of the transparent search, the results of TNVS are intuitionistic. The literature on facial landmark detection has highlighted some key reference points on the face that could help identity an individual, e.g., the Annotated Facial Landmarks in the Wild (AFLW) marked $21$ facial feature points (Köstinger et al., 2011), including $12$ points on the eyes and eyebrows, $3$ on the mouth, $3$ on the nose, $3$ on each ear, and $1$ on the chin. Checking whether our selected pixels cover most if not all of these key reference points, we can observe qualitatively whether the selection results provide insight to capture the key predictors of these observations.

One of the variable selection results is marked on the original images of ORL. For TNVS, we color the relevant pixels in yellow, the redundant pixels in red, and the unmarked pixels are conditionally independent pixels, as shown in Figure 4. For other methods, we color the selected pixels in yellow, and the unselected pixels are left unmarked. For ORL, TNVS selects $62$ relevant pixels, and removes $66$ redundant pixels, and $3\ 968$ conditionally independent pixels. FOCI selects $26$ pixels, WLS, DC-SIS, and SIRS select $62$ pixels, and PC-simple and TPC select $5$ pixels. Compared with the predictive results of baselines in Table LABEL:tab:appendix_A1 of Appendix B, TNVS achieves comparable or even better performance with more concise subset of predictors.

As is shown in Table LABEL:tab:appendix_A1 of Appendix B, TNVS and WLS are the best two methods on ORL dataset. Their selected pixels both cover the facial points of AFLW, which mainly lie in eye corners, eyebrow corners, tips of nose, corners of mouth, and chins. The accurate detection of key reference points intuitively explains the superior facial discrimination of TNVS and WLS. TNVS further reveals that the redundant pixels of ORL are concentrated mainly in the foreheads, indicating that pixels in these areas are linearly correlated with the relevant pixels. Most regions belong to the conditionally independent subset, indicating that only a very small number of pixels is adequate to explain the response to a certain extent. Some pixels around the eye corners are usually symmetric, and TNVS only select pixels on one side, since given these areas, the pixels in symmetric areas will be conditionally independent. FOCI selects less pixels than TNVS, and these two methods select similar regions of the faces, but FOCI cannot identify the redundant pixels. SIRS and DC-SIS select a large number of redundant pixels in the forehead, and they ignore some feature points in the eye corners, eyebrow corners, and tips of the nose. PC-simple and TPC omit many important facial feature points. Together with the results in Table LABEL:tab:appendix_A1 of Appendix B, this example shows that the predictors selected by TNVS can precisely predict the response, and the reserved or removed predictors can be explicitly interpreted.

We further demonstrate stability of TNVS as the three thresholds vary. These ranges should be tuned within certain ranges. If these thresholds are too large, the predictors in each subset will violate their definitions, i.e., the uninformative subset may include informative predictors, and the redundant subset may include predictors which are not linearly independent. If the method is only robust within small ranges, it will be easily influnced by noise in real data. Here, we take a large enough range for each threshold to show that the outperformance is persistently significant, and the proposed method is robust against hyperparameter variations.

First, the uninformative threshold is tuned within range $\alpha_{1}\in\{0,~{}0.01,\cdots,0.05\}$, and nonlinear correlations are evaluated between the response and the uninformative predictors $T_{n}(\mathbf{Y},\mathbf{X}_{\mathcal{A}_{1}})$, and the nonlinear partial correlations between the response and uninformative predictors given the selected predictors $T_{n}(\mathbf{Y},\mathbf{X}_{\mathcal{A}_{1}}\mid\mathbf{X}_{\mathcal{S}})$ for the four real datasets. If $T_{n}(\mathbf{Y},\mathbf{X}_{\mathcal{A}_{1}}\mid\mathbf{X}_{\mathcal{S}})$ is undefined, it can be regarded that $\mathbf{Y}$ is almost surely a measurable function of $\mathbf{X}_{\mathcal{S}}$, and we set $T_{n}(\mathbf{Y},\mathbf{X}_{\mathcal{A}_{1}}\mid\mathbf{X}_{\mathcal{S}})=0$. For warpAR10P and ORL, the numbers of uninformative predictors are $0$. For arcene and isolet, the numbers of identified uninformative predictors, their correlations and partial correlations are shown in Table 8. The number of uninformative predictors in arcene is the same when $\alpha_{1}\in[0,~{}0.03]$, and that in isolet is the same when $\alpha_{1}\in[0,~{}0.05]$. The correlations and partial correlations between the response and the uninformative predictors are close to $0$ in both datasets, which illustrates that the uninformative predictors have no positive effect on predicting the response. Identifying the uninformative subsets does not heavily depend on the settings of threshold $\alpha_{1}$, which shows the robustness of TNVS at the prefiltering step.

We next tune the remaining two parameters by taking the value of the relevant threshold $\alpha_{2}\in\{-0.05,~{}-0.04,\cdots,0\}$ and the value of the redundant threshold $\alpha_{3}\in\{0,~{}0.01,\cdots,0.05\}$. For each combination of $\alpha_{2}$ and $\alpha_{3}$, a 10-fold cross-validation is performed on the four real datasets to demonstrate the stability of TNVS as these two parameters vary. Fig. 5 shows the average cardinality of the selected subsets obtained by TNVS under different values of $\alpha_{2}$ and $\alpha_{3}$. The number of selected predictors fluctuates slightly in response to the change in $\alpha_{2}$ and $\alpha_{3}$, and there is no obvious and unified trend in different datasets. Fig. C1 of Appendix C further describes the prediction accuracy using the above selected predictors. As expected, the changes in both thresholds have little influence on the predictive results. In summary, the proposed TNVS is insensitive to changes in the thresholds within certain ranges.